tge.f90

Go to the documentation of this file.

!/===========================================================================/
! Copyright (c) 2007, The University of Massachusetts Dartmouth 
! Produced at the School of Marine Science & Technology 
! Marine Ecosystem Dynamics Modeling group
! All rights reserved.
!
! FVCOM has been developed by the joint UMASSD-WHOI research team. For 
! details of authorship and attribution of credit please see the FVCOM
! technical manual or contact the MEDM group.
!
! 
! This file is part of FVCOM. For details, see http://fvcom.smast.umassd.edu 
! The full copyright notice is contained in the file COPYRIGHT located in the 
! root directory of the FVCOM code. This original header must be maintained
! in all distributed versions.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" 
! AND ANY EXPRESS OR  IMPLIED WARRANTIES, INCLUDING,  BUT NOT  LIMITED TO,
! THE IMPLIED WARRANTIES OF MERCHANTABILITY AND  FITNESS FOR A PARTICULAR
! PURPOSE ARE DISCLAIMED.  
!
!/---------------------------------------------------------------------------/
! CVS VERSION INFORMATION
! $Id$
! $Name$
! $Revision$
!/===========================================================================/

   SUBROUTINE triangle_grid_edge

! NEEDS UPDATE TO FIX ERROR MESSAGES

!==============================================================================!
!  This program is used to define the non-overlapped, unstructured             !
!  triangular meshes used for flux computations. The mesh could be             !
!  created using the commerical software called "sms8.0" or other              !
!  mesh generation programs. The mesh file generated by sms8.0 can             !
!  be directly used for this subroutine, while the mesh file                   !
!  generated using other programs must be coverted the data format             !
!  to meet the required format used here.                                      !
!==============================================================================!
!     variable list:                                   !
!  vx(m)    :: vx(i) = x-coordinate of node i (input from mesh)                !
!  vy(m)    :: vy(i) = y-coordinate of node i (input from mesh)                !
!  nv(n,3)  :: nv(i:1-3) = 3 node numbers of element i                         !
!  xc(n)    :: xc(i) = x-coordinate of element i (calculated from vx)          !
!  yc(n)    :: yc(i) = y-coordinate of element i (calculated from vy)          !
!  cor(n)   :: cor(i) = f plane coriolis at element i                          !
!                                                                              !
!  nbe(n,3) :: nbe(i,1->3) = element index of 1->3 neighbors of element i      !
!  isbce(n) :: flag if element is on the boundary, see below for values        !
!  isonb(m) :: flag is node is on the boundary, see below for values           !
!                                                                              !
!  ntve(m)  :: the number of neighboring elements of node m                    !
!  nbve(m,ntve(m)) :: nbve(i,1->ntve(i)) = ntve elements containing node i     !
!  nbvt(m,ntve(m)) :: nbvt(i,j) = the node number of node i in element         ! 
!                     nbve(i,j) (has a value of 1,2,or 3)                      ! 
!                                                                              !
!   ne       ::   number of unique element edges                               !
!  iec(ne,2) :: nbe(i,1->2) cell number of cell(s) connected to edge i         !
!   isbc(ne) :: flag marking edge property                                     !
!              isbc(i) = 0:  element edge i is in the interior                 !
!              isbc(i) = 1:  element edge i is on the boundary                 !
!ienode(ne,2):: ienode(i,1->2) node numbers at each end of element edge i      ! 
!  xijc(ne)  :: xijc(i) = x-coordinate of mid point of element edge i          ! 
!  yijc(ne)  :: yijc(i) = y-coordinate of mid point of element edge i          ! 
!  dltxyc(ne):: dltxyc(i) = length of element edge i                           !
!  dltxc(ne) :: dltxc(i) = deltax (x-projection) of element edge i             !
!  dltyc(ne) :: dltyc(i) = deltay (y-projection) of element edge i             !
!  sitac(ne) :: sitac(i) =  arctg(dltyc,dltxc) (angle of inclination of edge)  !
!                                                                              !
!==============================================================================!
!     classification of the triangles nodes, and edges                         !
!                                                                              !
!     isonb(i)=0:  node in the interior computational domain                   !
!     isonb(i)=1:  node on the solid boundary                                  !
!     isonb(i)=2:  node on the open boundary                                   !
!                                                                              !
!     isbce(i)=0:  element in the interior computational domain                !
!     isbce(i)=1:  element on the solid boundary                               !
!     isbce(i)=2:  element on the open boundary                                !
!     isbce(i)=3:  element with 2 solid boundary edges                         !
!                                                                              !
!      isbc(i)=0:  element edge in interior                                    !
!      isbc(i)=1:  element edge on boundary                                    !
!==============================================================================!


!==============================================================================|
!   FIND NEIGHBORING ELEMENTS, MARK BOUNDARY NODES AND ELEMENTS                |
!                                          |
!   NBE(N,3) :  NBE(I,J) = IDENTITY OF NBOR ELMNT TO TRI I ON EDGE J           |
!   IBCE(N)  :  DESCRIBED IN SUBROUTINE HEADING                        |    
!   ISONB(M):  DESCRIBED IN SUBROUTINE HEADING                         |    
!==============================================================================|
   USE all_vars
   USE mod_spherical


   USE mod_par

   USE mod_obcs
   IMPLICIT NONE
   INTEGER, ALLOCATABLE, DIMENSION(:,:) :: TEMP,TEMP2,NB_TMP,ISET
   INTEGER I,J,DIF1,DIF2,DIF3,II,JJ,NTMP,NCNT,INEY,NFLAG
   INTEGER ITMP1,ITMP2,ITMP3,JN,JJB,IBCETMP,NCTMP,NCETMP,NPT
   INTEGER, ALLOCATABLE :: CELLS(:,:),CELLCNT(:),NBET(:,:)
   REAL(SP)  DTMP
   INTEGER N1,N2,N3,J1,J2,J3,SBUF, IERR
   
   ! STORAGE FOR ISONB, ISBCE TO RUN EXCHANGE AND SET HALO'S
   REAL(SP), ALLOCATABLE :: FTEMP(:)

   Character(len=8)  :: tmpstr
   Character(len=200):: errstr


   REAL(SP) DELTX,DELTY,ALPHA1,ALPHA2

!----------------------------REPORT--------------------------------------------!
   IF (dbg_set(dbg_log))THEN
      WRITE(ipt,*  )'!'
      WRITE(ipt,*)'!           SETTING UP TRIS/ELEMENTS/EDGES/CVS          '
      WRITE(ipt,*  )'!'
   END IF
!----------------------------INITIALIZE----------------------------------------!
   
   isbce = 0
   isonb = 0
   nbe   = 0

!
!----DETERMINE NBE(i=1:n,j=1:3): INDEX OF 1 to 3 NEIGHBORING ELEMENTS----------!
!
   ALLOCATE(nbet(nt,3)) ; nbet = 0
   ALLOCATE(cells(mt,50)) ; cells = 0
   ALLOCATE(cellcnt(mt))  ; cellcnt = 0
   DO i=1,nt
     n1 = nv(i,1) ; cellcnt(n1) = cellcnt(n1)+1
     n2 = nv(i,2) ; cellcnt(n2) = cellcnt(n2)+1
     n3 = nv(i,3) ; cellcnt(n3) = cellcnt(n3)+1
     cells(nv(i,1),cellcnt(n1)) = i
     cells(nv(i,2),cellcnt(n2)) = i
     cells(nv(i,3),cellcnt(n3)) = i
   END DO
   if(maxval(cellcnt) > 50)write(ipt,*)'bad',maxval(cellcnt)
   DO i=1,nt
     n1 = nv(i,1)
     n2 = nv(i,2)
     n3 = nv(i,3)
     DO j1 = 1,cellcnt(n1) 
     DO j2 = 1,cellcnt(n2) 
       IF((cells(n1,j1) == cells(n2,j2)).AND. cells(n1,j1) /= i)nbe(i,3) = cells(n1,j1)
     END DO
     END DO
     DO j2 = 1,cellcnt(n2) 
     DO j3 = 1,cellcnt(n3) 
       IF((cells(n2,j2) == cells(n3,j3)).AND. cells(n2,j2) /= i)nbe(i,1) = cells(n2,j2)
     END DO
     END DO
     DO j1 = 1,cellcnt(n1) 
     DO j3 = 1,cellcnt(n3) 
       IF((cells(n1,j1) == cells(n3,j3)).AND. cells(n1,j1) /= i)nbe(i,2) = cells(n3,j3)
     END DO
     END DO
   END DO
   DEALLOCATE(cells,cellcnt)

   ! OLD - USED TO MAKE GLOBAL INDEX ARRAYS TO SAVE PARALLEL OUTPUT
   ! REPLACED IN MOD_NCDIO WITH MORE GENERIC TYPE GRID DATA!
!   IF (PAR) THEN
!      ALLOCATE(NBEGL(0:NT,3)); NBEGL      = 0
!      DO I = 1,NT
!         NBEGL(I,:) = EGID_X(NBE(I,:))
!      END DO
!   ELSE
!      NBEGL => NBE
!   END IF

   IF (dbg_set(dbg_log))WRITE(ipt,*)  '!  NEIGHBOR FINDING      :    COMPLETE'
!
!--ENSURE ALL ELEMENTS HAVE AT LEAST ONE NEIGHBOR------------------------------!
!
   nflag = 0
   DO i=1,nt
     IF(sum(nbe(i,1:3))==0)THEN 
       nflag = 1
       WRITE(ipt,*)'ELEMENT ',i,' AT ',xc(i),yc(i),' HAS NO NEIGHBORS'
       CALL pstop
     END IF
   END DO
   IF(nflag == 1) CALL pstop
     
!
!----IF ELEMENT ON BOUNDARY SET ISBCE(I)=1 AND ISONB(J)=1 FOR BOUNDARY NODES J-!
!

   DO i=1,nt 
     IF(min(nbe(i,1),nbe(i,2),nbe(i,3))==0)THEN    !!ELEMENT ON BOUNDARY
       isbce(i) = 1
       IF(nbe(i,1) == 0)THEN 
         isonb(nv(i,2)) = 1 ; isonb(nv(i,3)) = 1
       END IF
       IF(nbe(i,2) ==0) THEN
         isonb(nv(i,1)) = 1 ; isonb(nv(i,3)) = 1
       END IF
       IF(nbe(i,3) ==0) THEN
         isonb(nv(i,1)) = 1 ; isonb(nv(i,2)) = 1
       END IF
     END IF
   END DO
   IF (dbg_set(dbg_log)) WRITE(ipt,*)  '!  ISONB SETTING         :    COMPLETE'
        
!==============================================================================|
!             DEFINE NTVE, NBVE, NBVT                                          !
!                                                                              !
! ntve(1:m):           total number of the surrounding triangles               !
!                      connected to the given node                             !
! nbve(1:m, 1:ntve+1): the identification number of surrounding                !
!                      triangles with a common node (counted clockwise)        !
! nbvt(1:m,ntve(1:m)): the idenfication number of a given node over            !
!                      each individual surrounding triangle(counted            !
!                      clockwise)                                              !
! ntsn(1:m):           total number of surrounding nodes                       !
! nbsn(1:m, ntsn):     the identification number of surrounding nodes          !
!                      (counted clockwise)                                     !
! nbse(1:m,2*ntsn):    the identification number of control volume s           !
!                      edges between two neighbor nodes                        !
!==============================================================================|

!
!----DETERMINE MAX NUMBER OF SURROUNDING ELEMENTS------------------------------!
!
   mx_nbr_elem = 0
   DO i=1,m
     ncnt = 0
     DO j=1,nt
!       IF( (NV(J,1)-I)*(NV(J,2)-I)*(NV(J,3)-I)==0) NCNT = NCNT + 1
!       IF( (NV(J,1)-I) ==0 .OR. (NV(J,2)-I) ==0 .OR. (NV(J,3)-I)==0) NCNT = NCNT + 1
       IF( float(nv(j,1)-i)*float(nv(j,2)-i)*float(nv(j,3)-i) == 0.0_sp) &
         ncnt = ncnt + 1
     END DO
     mx_nbr_elem = max(mx_nbr_elem,ncnt)
   END DO
!   WRITE(IPT,*) 'MAXIMUM NUMBER OF NEIGHBOR ELEMENTS',MX_NBR_ELEM


!
!----ALLOCATE ARRAYS BASED ON MX_NBR_ELEM--------------------------------------!
! 
   ALLOCATE(nbve(m,mx_nbr_elem+1)); nbve = 0
   ALLOCATE(nbvt(m,mx_nbr_elem+1)); nbvt = 0
!   ALLOCATE(NBSN(M,MX_NBR_ELEM+2))
   ALLOCATE(nbsn(m,mx_nbr_elem+3)); nbsn = 0
!
!--DETERMINE NUMBER OF SURROUNDING ELEMENTS FOR NODE I = NTVE(I)---------------!
!--DETERMINE NBVE - INDICES OF NEIGHBORING ELEMENTS OF NODE I------------------!
!--DETERMINE NBVT - INDEX (1,2, or 3) OF NODE I IN NEIGHBORING ELEMENT---------!
!
       
   DO i=1,m
     ncnt=0
     DO j=1,nt
!       IF( (NV(J,1)-I) == 0 .OR.  (NV(J,2)-I) == 0 .OR. (NV(J,3)-I) == 0)THEN 
        IF (float(nv(j,1)-i)*float(nv(j,2)-i)*float(nv(j,3)-i) == 0.0_sp)THEN
         ncnt = ncnt+1
         nbve(i,ncnt)=j
         IF((nv(j,1)-i) == 0) nbvt(i,ncnt)=1
         IF((nv(j,2)-i) == 0) nbvt(i,ncnt)=2
         IF((nv(j,3)-i) == 0) nbvt(i,ncnt)=3
       END IF
     ENDDO
     ntve(i)=ncnt
   ENDDO
!
!--Reorder Order Elements Surrounding a Node to Go in a Cyclical Procession----!
!--Determine NTSN  = Number of Nodes Surrounding a Node (+1)-------------------!
!--Determine NBSN  = Node Numbers of Nodes Surrounding a Node------------------!
!
!Lettmann&JQI   ALLOCATE(NB_TMP(M,MX_NBR_ELEM+1))
   ALLOCATE(nb_tmp(mx_nbr_elem+1,2))
   DO i=1,m
     IF(isonb(i) == 0) THEN
       nb_tmp(1,1)=nbve(i,1)
       nb_tmp(1,2)=nbvt(i,1)
       DO j=2,ntve(i)+1
         ii=nb_tmp(j-1,1)
         jj=nb_tmp(j-1,2)
         nb_tmp(j,1)=nbe(ii,jj+1-int((jj+1)/4)*3)
         jj=nb_tmp(j,1)
         IF((nv(jj,1)-i) == 0) nb_tmp(j,2)=1
         IF((nv(jj,2)-i) == 0) nb_tmp(j,2)=2
         IF((nv(jj,3)-i) == 0) nb_tmp(j,2)=3
       ENDDO

       DO j=2,ntve(i)+1
         nbve(i,j)=nb_tmp(j,1)
       ENDDO

       DO j=2,ntve(i)+1
         nbvt(i,j)=nb_tmp(j,2)
       ENDDO

       ntmp=ntve(i)+1
       IF(nbve(i,1) /= nbve(i,ntmp)) THEN
          print*, myid,ngid(i),ntmp,'NBVE(I,nTMP) NOT CORRECT!!'
          print*, "NBVE(I,:)=",egid(nbve(i,:))
          CALL pstop
       ENDIF

       IF(nbvt(i,1) /= nbvt(i,ntmp)) THEN
          print*, ngid(i),'NBVT(I) NOT CORRECT!!'
          print*, "NBVT(I,:)=",nbvt(i,:)
          CALL pstop
       END IF

       ntsn(i)=ntve(i)

       DO j=1,ntsn(i)
         ii=nbve(i,j)
         jj=nbvt(i,j)
         nbsn(i,j)=nv(ii,jj+1-int((jj+1)/4)*3)
       ENDDO

       ntsn(i)=ntsn(i)+1
       nbsn(i,ntsn(i))=nbsn(i,1)

     ELSE 
           jjb=0

       DO j=1,ntve(i)
         jj=nbvt(i,j)
         IF(nbe(nbve(i,j),jj+2-int((jj+2)/4)*3) == 0) THEN
           jjb=jjb+1
           nb_tmp(jjb,1)=nbve(i,j)
           nb_tmp(jjb,2)=nbvt(i,j)
         END IF
       ENDDO

       IF(jjb /= 1) THEN
         WRITE(ipt,*) 'ERROR IN ISONB !,I,J', i,j
         CALL fatal_error("ERROR IN TGE DETERMINING ISONB")
       END IF

       DO j=2,ntve(i)
         ii=nb_tmp(j-1,1)
         jj=nb_tmp(j-1,2)
         nb_tmp(j,1)=nbe(ii,jj+1-int((jj+1)/4)*3)
         jj=nb_tmp(j,1)
         IF((nv(jj,1)-i) == 0) nb_tmp(j,2)=1
         IF((nv(jj,2)-i) == 0) nb_tmp(j,2)=2
         IF((nv(jj,3)-i) == 0) nb_tmp(j,2)=3
       ENDDO

       DO j=1,ntve(i)
         nbve(i,j)=nb_tmp(j,1)
         nbvt(i,j)=nb_tmp(j,2)
       ENDDO

       nbve(i,ntve(i)+1)=0
       ntsn(i)=ntve(i)+1
       nbsn(i,1)=i

       DO j=1,ntsn(i)-1
         ii=nbve(i,j)
         jj=nbvt(i,j)
         nbsn(i,j+1)=nv(ii,jj+1-int((jj+1)/4)*3)
       ENDDO

       j=ntsn(i)
       ii=nbve(i,j-1)
       jj=nbvt(i,j-1)
       nbsn(i,j+1)=nv(ii,jj+2-int((jj+2)/4)*3)
       ntsn(i)=ntsn(i)+2
       nbsn(i,ntsn(i))=i
     END IF
   END DO
   DEALLOCATE(nb_tmp)
   IF(mx_nbr_elem+3 -maxval(ntsn) < 0)THEN
      WRITE(ipt,*)'CHECK NTSN/NBSN',maxval(ntsn),mx_nbr_elem+3
      CALL pstop
   END IF
  

   ! OLD - USED TO MAKE GLOBAL INDEX ARRAYS TO SAVE PARALLEL OUTPUT
   ! REPLACED IN MOD_NCDIO WITH MORE GENERIC TYPE GRID DATA!
!   IF(PAR) THEN
!      ALLOCATE(NBSNGL(M,MX_NBR_ELEM+3))
!      ALLOCATE(NBVEGL(M,MX_NBR_ELEM+1))
!      DO I = 1,M
!         NBVEGL(I,:)=EGID_X(NBVE(I,:))
!
!         NBSNGL(I,:) = NGID_X(NBSN(I,:))
!      END DO
!   ELSE
!      NBVEGL => NBVE
!
!      NBSNGL =>NBSN
!   END IF

   IF (dbg_set(dbg_log))WRITE(ipt,*)  '!  NBVE/NBVT             :    COMPLETE'
     

!==============================================================================!
!  Define the parameters of each triangular edge                               !
!                                                                              !
!  ne           :    number of unique element edges                            !
!  iec(1:ne,1:2):    counting number identifying two connected cells           !
!  isbc(1:ne):       0: triangle s edge in the interior                        !
!                    1: triangle s edge on the boundary                        !
!  ienode(1:ne,1:2): the identification number of two end points of a          !
!                    edge                                                      !
!  xijc(1:ne):       the x-coordinate location of the middle points            !
!                    of a edge                                                 !
!  yijc(1:ne):       the y-coordinate location of the middle points            !
!                    of a edge                                                 !
!  dltxyc(1:ne):     length of the edge                                        !
!  dltxc(1:ne):      vx(ienode(i,2))-vx(idnode(i,1))                           !
!  dltyc(1:ne):      vy(ienode(i,2))-vy(idnode(i,1))                           !
!  sitac(1:ne):      arctg(dltyc,dltxc)                                        !
!==============================================================================!

   ALLOCATE(iset(nt,3),temp((nt)*3,2),temp2((nt)*3,2))
   iset = 0
   ne = 0
   temp = 0 
   temp2 = 0
   DO i=1,nt
     DO j=1,3
       IF(iset(i,j) == 0)THEN
         ne   = ne + 1
         iney = nbe(i,j)
         iset(i,j) = 1
         DO jn=1,3
           IF(i == nbe(iney,jn)) iset(iney,jn) = 1
         END DO
         temp(ne,1) = i ; temp(ne,2) = iney
         temp2(ne,1) = nv(i,j+1-int((j+1)/4)*3)
         temp2(ne,2) = nv(i,j+2-int((j+2)/4)*3)
       END IF
     END DO
   END DO
   DEALLOCATE(iset)
!
!--ALLOCATE ARRAYS REQUIRING NUMBER OF EDGES-----------------------------------!
!
   ALLOCATE(iec(ne,2))
   ALLOCATE(ienode(ne,2))
   ALLOCATE(xijc(ne))
   ALLOCATE(yijc(ne))
   ALLOCATE(dltxyc(ne))
   ALLOCATE(dltxc(ne))
   ALLOCATE(dltyc(ne))
   ALLOCATE(sitac(ne))
   ALLOCATE(isbc(ne))


   iec(:,1) = temp(1:ne,1)
   iec(:,2) = temp(1:ne,2)
   ienode(:,1) = temp2(1:ne,1)
   ienode(:,2) = temp2(1:ne,2)


   DEALLOCATE(temp,temp2)
         
         
!
!------MARK ELEMENT EDGES THAT ARE ON THE BOUNDARY-----------------------------!
!
   isbc = 0
   DO i=1,ne
     IF((iec(i,1) == 0) .OR. (iec(i,2) == 0)) isbc(i) = 1 
   END DO

!
!------CALCULATE ELEMENT EDGE METRICS------------------------------------------!
!
   DO i=1,ne
     dltxc(i) =  vx(ienode(i,2))-vx(ienode(i,1))
     dltyc(i) =  vy(ienode(i,2))-vy(ienode(i,1))
     xijc(i)  = (vx(ienode(i,1))+vx(ienode(i,2)))/2.0_sp
     yijc(i)  = (vy(ienode(i,1))+vy(ienode(i,2)))/2.0_sp
     dltxyc(i)= sqrt(dltxc(i)**2+dltyc(i)**2)
     sitac(i) = atan2(dltyc(i),dltxc(i))
   END DO
    IF (dbg_set(dbg_log))WRITE(ipt,*)  '!  EDGE SETUP            :    COMPLETE'


!==============================================================================!
!  read triangular mesh parameters on open boundary :                          !
!  iobce:   number of open boundary cells.                                     !
!  isbcn:   number of open boundary nodes.                                     !
!  i_obc_e: counter number of open boundary cells                              !
!  i_obc_n: counter number of open boundary nodes                              !
!==============================================================================!

!
!----TRAVERSE  BOUNDARY NODE NUMBERS AND SET ISONB(NODE)=2---------------------!
!
   DO i=1,iobcn
     isonb(i_obc_n(i))=2
   ENDDO

!
!---- SET HALO VALUES IF PAR
! 


!
!----DETERMINE IF ELEMENT IS ON OPEN BOUNDARY (CONTAINS EDGE ON OPEN BOUNDARY)-!
!
   ibcetmp=0
   DO i=1,n
     itmp1=isonb(nv(i,1))
     itmp2=isonb(nv(i,2))
     itmp3=isonb(nv(i,3))

     IF(sum(isonb(nv(i,1:3))) == 4) THEN
       isbce(i)=2
       ibcetmp =ibcetmp+1
     ELSE IF(sum(isonb(nv(i,1:3))) > 4) THEN
       print*,'SORRY, THE BOUNDARY CELL',i,'IS NOT GOOD FOR MODEL.'
       print*,'IT HAS EITHER TWO SIDES OF OPEN BOUNDARY OR ONE OPEN BOUNDARY'
       print*,'AND ONE SOLID BOUNDARY. PLEASE CHECK AND MODIFIED IT.'
       print*,'THIS MESSAGE IS IN SUBROUTINE TRIANGLE_GRID_EDGE (TGE.F)'
       print*,'STOP RUNNING...'
       CALL pstop

     END IF
   END DO

   DO i=1,nt
     IF((nbe(i,1)+nbe(i,2)+nbe(i,3) == 0).AND.(isbce(i) /= 2)) isbce(i)=3
     IF((nbe(i,1)+nbe(i,2) == 0).AND.(isbce(i) /= 2)) isbce(i)=3
     IF((nbe(i,2)+nbe(i,3) == 0).AND.(isbce(i) /= 2)) isbce(i)=3
     IF((nbe(i,1)+nbe(i,3) == 0).AND.(isbce(i) /= 2)) isbce(i)=3
   ENDDO

! SET HALO VALUES IF PAR


!==============================================================================!
!  xije(1:nc,1:2):  the x coordinate locations of starting and ending          !
!                   points of the control volumes edges                        !
!  yije(1:nc,1:2):  the y coordinate locations of starting and ending          !
!                   points of the control volumes edges                        !
!  niec(1:nc,1:2):  the counting number of left and right nodes                !
!                   conected to this control volumes edge from                 !
!                   starting point to ending point                             !
!  dltxe(1:nc):     the x distance of individual edges                         !
!  dltye(1:nc)      the y distance of individual edges                         !
!  dltxye(1:nc):    the length of individual edges                             !
!  ntrg(1:nc)  :    element associated with this control volume edge           !
!==============================================================================!
   nctmp  = 0
   ncetmp = 0 

   DO i=1,ne
     IF(isbc(i) == 0) THEN
       IF(iec(i,1) <= n)THEN
         nctmp=nctmp+1
         npt  =nctmp
       ELSE
         ncetmp = ncetmp + 1
         npt    = ncetmp+(3*n)
       END IF
       xije(npt,1) = xc(iec(i,1))
       yije(npt,1) = yc(iec(i,1))
       xije(npt,2) = xijc(i)
       yije(npt,2) = yijc(i)
       niec(npt,1) = ienode(i,1)
       niec(npt,2) = ienode(i,2)
       ntrg(npt)   = iec(i,1)
       dltxe(npt)  = xije(npt,2)-xije(npt,1)
       dltye(npt)  = yije(npt,2)-yije(npt,1)
       dtmp        = dltxe(npt)*dltxe(npt)+dltye(npt)*dltye(npt)
       dltxye(npt) = sqrt(dtmp)
       sitae(npt)  = atan2(dltye(npt),dltxe(npt))

       IF(iec(i,2) <= n)THEN
         nctmp=nctmp+1
         npt  =nctmp
       ELSE
         ncetmp = ncetmp + 1
         npt    = ncetmp+(3*n)
       END IF
       xije(npt,1)=xc(iec(i,2))
       yije(npt,1)=yc(iec(i,2))
       xije(npt,2)=xijc(i)
       yije(npt,2)=yijc(i)
       niec(npt,1)=ienode(i,2)
       niec(npt,2)=ienode(i,1)
       ntrg(npt)=iec(i,2)
       dltxe(npt)=xije(npt,2)-xije(npt,1)
       dltye(npt)=yije(npt,2)-yije(npt,1)
       dtmp=dltxe(npt)*dltxe(npt)+dltye(npt)*dltye(npt)
       dltxye(npt)=sqrt(dtmp)
       sitae(npt)=atan2(dltye(npt),dltxe(npt))
     ELSE IF(isbc(i) == 1) THEN
       IF(iec(i,1) <= n)THEN
         nctmp=nctmp+1
         npt  =nctmp
       ELSE
         ncetmp = ncetmp + 1
         npt    = ncetmp+(3*n)
       END IF
       IF(iec(i,1) == 0) THEN
         print*, i,'IEC(I,1)===0'
         CALL pstop
       END IF
       xije(npt,1)=xc(iec(i,1))
       yije(npt,1)=yc(iec(i,1))
       xije(npt,2)=xijc(i)
       yije(npt,2)=yijc(i)
       niec(npt,1)=ienode(i,1)
       niec(npt,2)=ienode(i,2)
       ntrg(npt)=iec(i,1)
       dltxe(npt)=xije(npt,2)-xije(npt,1)
       dltye(npt)=yije(npt,2)-yije(npt,1)
       dtmp=dltxe(npt)*dltxe(npt)+dltye(npt)*dltye(npt)
       dltxye(npt)=sqrt(dtmp)
       sitae(npt)=atan2(dltye(npt),dltxe(npt))
     ELSE
       WRITE(ipt,*) 'ISBC(I) NOT CORRECT, I==',i
       CALL pstop
     END IF
   ENDDO

   ncv_i  = nctmp
   ncv    = ncetmp+nctmp
      
   IF(ncv /= 3*(nt)) THEN
     print*,'NCV IS NOT CORRECT, PLEASE CHECK THE SETUP'
     CALL pstop
   END IF
   IF(ncv_i /= 3*n) THEN
     print*,'NCV_I IS NOT CORRECT, PLEASE CHECK THE SETUP'
      CALL pstop
   END IF

   DO i=1,ncv_i
     IF(niec(i,1) > m .OR. niec(i,2) > m)THEN
       write(ipt,*)'problemas',niec(i,1),niec(i,2),m
       CALL pstop
     END IF
   END DO

!==============================================================================!
!  nisbce_1/nisbce_2/nisbce_3:  number of elements with isbce of 1,2,3         !
!  lisbce_1/lisbce_2/lisbce_3:  list of elements with isbce of 1,2,3           !
!  epor                      :  element porosity (=0 if isbce = 2)             !
!==============================================================================!

!
!--COUNT NUMBER OF ELEMENTS OF EACH TYPE (ISBCE=1,2,3)-------------------------!
!
   nisbce_1 = 0 ; nisbce_2 = 0 ; nisbce_3 = 0
   DO i=1,n 
     IF(isbce(i) == 1) nisbce_1 = nisbce_1 + 1
     IF(isbce(i) == 2) nisbce_2 = nisbce_2 + 1
     IF(isbce(i) == 3) nisbce_3 = nisbce_3 + 1
   END DO

!
!--ALLOCATE ELEMENT TYPE ARRAYS LISBCE_1,LISBCE_2,LISBCE_3---------------------!
!
   IF(nisbce_1 > 0)THEN
     ALLOCATE( lisbce_1(nisbce_1) )
   ELSE
!     WRITE(IPT,*)  '!  WARNING               :    NO ELEMENTS WITH ISBCE=1'
   END IF
     
   IF(nisbce_2 > 0)THEN
     ALLOCATE( lisbce_2(nisbce_2) )
   ELSE
!     WRITE(IPT,*)  '!  WARNING               :    NO ELEMENTS WITH ISBCE=2'
   END IF
     
   IF(nisbce_3 > 0)THEN
     ALLOCATE( lisbce_3(nisbce_3) )
   ELSE
!     WRITE(IPT,*)  '!  WARNING               :    NO ELEMENTS WITH ISBCE=3'
   END IF

!
!--LOAD ELEMENT TYPE ARRAYS LISBCE_1,LISBCE_2,LISBCE_3--------------------------!
!
   nisbce_1 = 0 ; nisbce_2 = 0 ; nisbce_3 = 0
   DO i=1,n 
     IF(isbce(i) == 1) THEN
       nisbce_1 = nisbce_1 + 1
       lisbce_1(nisbce_1) = i
     END IF
     IF(isbce(i) == 2) THEN
       nisbce_2 = nisbce_2 + 1
       lisbce_2(nisbce_2) = i
     END IF
     IF(isbce(i) == 3) THEN
       nisbce_3 = nisbce_3 + 1
       lisbce_3(nisbce_3) = i
     END IF
   END DO

!
!--SET ELEMENT POROSITY---------------------------------------------------------!
!
   ALLOCATE(epor(0:nt)) ; epor = 1.0_sp
   DO i=1,n
     IF(isbce(i) == 2) epor(i) = 0.0_sp
   END DO
     
   IF (dbg_set(dbg_log))WRITE(ipt,*)  '!  NISBCE/LISBCE/EPOR    :    COMPLETE'
   IF (dbg_set(dbg_log))WRITE(ipt,*)  '!  TRIS/EDGES/CVOLS      :    COMPLETE'

   RETURN
   END SUBROUTINE triangle_grid_edge
!==============================================================================!