#include "w3macros.h" !/ ------------------------------------------------------------------- / MODULE W3TRIAMD !/ ------------------------------------------------------------------- !/ +-----------------------------------+ !/ | WAVEWATCH III NOAA/NCEP | !/ | F. Ardhuin and A. Roland | !/ | FORTRAN 90 | !/ | Last update : 26-Jan-2014| !/ +-----------------------------------+ !/ !/ 15-Mar-2007 : Origination. ( version 3.13 ) !/ 25-Aug-2011 : Modification of boundary treatment ( version 4.04 ) !/ 30-Aug-2012 : Automatic detection of open BC ( version 4.08 ) !/ 02-Sep-2012 : Clean up of open BC for UG grids ( version 4.08 ) !/ 14-Oct-2013 : Correction of latitude factor ( version 4.12 ) !/ 26-Jan-2014 : Correction interpolation weights ( version 4.18 ) !/ 21-Apr-2016 : New algorithm to detect boundary ( version 5.12 ) !/ ! ! 1. Purpose : ! ! Reads triangle and unstructured grid information ! ! 2. Method : ! ! Look for namelist with name NAME in unit NDS and read if found. ! ! 3. Parameters : ! ! 4. Subroutines used : ! ! Name Type Module Description ! ------------------------------------------------------------------------------------ ! READTRI Subr. Internal Read unstructured grid data from .grd .tri formatted files. ! READMSH Subr. Id. Read unstructured grid data from MSH format ! COUNT Subr. Internal Count connection. ! SPATIAL_GRID Subr. Id. Calculate surfaces. ! NVECTRI Subr. Id. Define cell normals and angles and edge length ! COORDMAX Subr. Id. Calculate useful grid elements ! AREA_SI Subr. Id. Define Connections ! ------------------------------------------------------------------------------------ ! ! 5. Called by : ! ! Program in which it is contained. ! ! 6. Error messages : ! ! 7. Remarks : ! The only point index which is needed is IX and NX stands for the total number of grid point. ! IY and NY are not needed anymore, they are set to 1 in the unstructured case ! Some noticeable arrays are: ! XYB : give the 2D coordinates of all grid points ! TRIGP : give the vertices of each triangle ! 8. Structure : ! ! 9. Switches : ! !/PR3 : Enables unstructured meshes (temporary, will be replace by Unstructured switch) ! 10. Source code : ! !/ ------------------------------------------------------------------- / PUBLIC ! USE CONSTANTS ! USE W3GDATMD, ONLY: W3NMOD, W3SETG ! USE W3ODATMD, ONLY: W3NOUT, W3SETO, W3DMO5 ! USE W3IOGRMD, ONLY: W3IOGR ! USE W3SERVMD, ONLY: ITRACE, NEXTLN, EXTCDE !!/S USE W3SERVMD, ONLY: STRACE ! USE W3ARRYMD, ONLY: INA2R, INA2I !!/T USE W3ARRYMD, ONLY: PRTBLK ! USE W3DISPMD, ONLY: DISTAB ! USE W3GDATMD ! USE W3ODATMD, ONLY: NDSE, NDST, NDSO ! USE W3ODATMD, ONLY: NBI, NBI2, NFBPO, NBO, NBO2, FLBPI, FLBPO, & ! IPBPO, ISBPO, XBPO, YBPO, RDBPO, FNMPRE !--------------------------------------------------------------------- ! !C integer :: node_num integer :: dim_num integer :: triangle_order integer :: triangle_num integer :: bound_edge_num integer :: bound_num !C logical,save, allocatable :: edge_boundary(:) logical,save, allocatable :: node_boundary(:) integer,save, allocatable :: edge_nums(:) integer,save, allocatable :: boundary_node_index(:) !C integer,save, allocatable :: triangle_node(:,:) integer,save, allocatable :: edge(:,:) integer,save, allocatable :: edge_index(:,:) integer,save, allocatable :: iobp_aux(:) INTEGER, SAVE :: N_OUTSIDE_BOUNDARY INTEGER, SAVE, ALLOCATABLE :: OUTSIDE_BOUNDARY(:) real (kind = 8), save, allocatable :: node_xy(:,:) real (kind = 8), save, allocatable :: edge_angle(:,:) CONTAINS !/ -------------------------------------------------------------------/ SUBROUTINE READMSH(NDS,FNAME) !/ ------------------------------------------------------------------- !/ +-----------------------------------+ !/ | WAVEWATCH III NOAA/NCEP | !/ | F. Ardhuin | !/ | FORTRAN 90 | !/ | Last update : 06-Jun-2018| !/ +-----------------------------------+ !/ !/ 15-Feb-2008 : Origination. ( version 3.13 ) !/ 25-Aug-2011 : Change of method for IOBPD ( version 4.04 ) !/ 06-Jun-2018 : Add DEBUGINIT/PDLIB/DEBUGSTP/DEBUGSETIOBP !/ ( version 6.04 ) !/ ! ! 1. Purpose : ! ! Reads triangle and unstructured grid information from GMSH files ! Calls the subroutines needed to compute grid connectivity ! ! 2. Method : ! ! Look for namelist with name NAME in unit NDS and read if found. ! ! 3. Parameters : ! ! Parameter list ! ---------------------------------------------------------------- ! NDS Int. I Data set number used for search. ! NAME C*4 I Name of namelist. ! STATUS C*20 O Status at end of routine, ! '(default values) ' if no namelist found. ! '(user def. values)' if namelist read. ! ---------------------------------------------------------------- ! ! 4. Subroutines used : ! ! Name Type Module Description ! ------------------------------------------------------------------------------------ ! NEXTLN Subr. ! COUNT Subr. Internal Count connection. ! SPATIAL_GRID Subr. Id. Calculate surfaces. ! NVECTRI Subr. Id. Define cell normals and angles and edge length ! COORDMAX Subr. Id. Calculate useful grid elements ! AREA_SI Subr. Id. Define Connections ! ---------------------------------------------------------------- ! ! 5. Called by : ! Name Type Module Description ! ---------------------------------------------------------------- ! W3GRID Prog. Model configuration program ! ---------------------------------------------------------------- ! ! 6. Error messages : ! ! 7. Remarks : ! The only point index which is needed is IX and NX stands for the total number of grid point. ! IY and NY are not needed anymore, they are set to 1 in the unstructured case ! Some noticeable arrays are: ! XYB : give the 2D coordinates of all grid points ! TRIGP : give the vertices of each triangle ! GMSH file gives too much information that is not necessarily required so data processing is needed (data sort and nesting). ! 8. Structure : ! ! 9. Switches : ! ! 10. Source code : ! !/ ------------------------------------------------------------------- / USE W3ODATMD, ONLY: NDSE, NDST, NDSO USE W3GDATMD USE W3SERVMD, ONLY: ITRACE, NEXTLN, EXTCDE USE CONSTANTS, only: LPDLIB USE W3ODATMD, ONLY: IAPROC ! IMPLICIT NONE !/ !/ Parameter list !/ INTEGER, INTENT(IN) :: NDS CHARACTER(60), INTENT(IN) :: FNAME !/ !/ local parameters !/ INTEGER :: i,j,k, NODES, NELTS, ID, KID INTEGER :: ID1, ID2, KID1, ITMP(3) INTEGER :: I1, I2, I3 INTEGER(KIND=4) :: Ind,eltype,ntag, INode CHARACTER :: COMSTR*1, SPACE*1 = ' ', CELS*64 REAL, ALLOCATABLE :: TAGS(:) CHARACTER(LEN=64), ALLOCATABLE :: ELS(:) CHARACTER(LEN=120) :: LINE CHARACTER(LEN=50) :: CHTMP CHARACTER(LEN=10) :: A, B, C INTEGER,ALLOCATABLE :: NELS(:), TRIGPTMP1(:,:), TRIGPTMP2(:,:) INTEGER(KIND=4),ALLOCATABLE :: IFOUND(:), VERTEX(:), BOUNDTMP(:) DOUBLE PRECISION, ALLOCATABLE :: XYBTMP1(:,:),XYBTMP2(:,:) REAL :: z OPEN(NDS,FILE = FNAME,STATUS='old') READ (NDS,'(A)') COMSTR IF (COMSTR.EQ.' ') COMSTR = '$' CALL NEXTLN(COMSTR, NDS, NDSE) READ(NDS,*) i,j,k CALL NEXTLN(COMSTR, NDS, NDSE) LPDLIB = .FALSE. ! ! read number of nodes and nodes from Gmsh files ! READ(NDS,*) NODES ALLOCATE(XYBTMP1(NODES,3)) DO I= 1, NODES READ(NDS,*) j, XYBTMP1(I,1), XYBTMP1(I,2), XYBTMP1(I,3) END DO ! ! read number of elements and elements from Gmsh files ! ALLOCATE(BOUNDTMP(NODES)) N_OUTSIDE_BOUNDARY = 0 CALL NEXTLN(COMSTR, NDS, NDSE) READ(NDS,*) NELTS ALLOCATE(TRIGPTMP1(NELTS, 3)) J = 0 DO I= 1, NELTS READ(NDS,'(A100)') LINE READ(LINE,*) Ind,eltype,ntag ALLOCATE(TAGS(ntag)) SELECT CASE (eltype) ! ! eltype = 15 : boundary points (this is used to make the difference ! between the outside polygon and islands) ! CASE(15) READ(LINE,*) Ind,eltype,ntag,TAGS,INODE N_OUTSIDE_BOUNDARY = N_OUTSIDE_BOUNDARY +1 BOUNDTMP(N_OUTSIDE_BOUNDARY)=INODE ! ! eltype = 2 : triangles ! CASE (2) J = J + 1 READ(LINE,*) Ind,eltype,ntag,tags,ITMP TRIGPTMP1(J,1:3) = ITMP END SELECT DEALLOCATE(TAGS) END DO ! ! organizes the grid data structure ! ALLOCATE(OUTSIDE_BOUNDARY(N_OUTSIDE_BOUNDARY)) OUTSIDE_BOUNDARY(:)=BOUNDTMP(1:N_OUTSIDE_BOUNDARY) NTRI = J ALLOCATE(IFOUND(NODES)) IFOUND = 0 ! ! Verifies that the nodes are used in at least one triangle ! DO K = 1, NTRI I1 = TRIGPTMP1(K,1) I2 = TRIGPTMP1(K,2) I3 = TRIGPTMP1(K,3) IFOUND(I1)= IFOUND(I1) + 1 IFOUND(I2)= IFOUND(I2) + 1 IFOUND(I3)= IFOUND(I3) + 1 END DO J = 0 ALLOCATE(TRIGPTMP2(NTRI,3),VERTEX(NODES),XYBTMP2(NODES,3)) VERTEX(:)=0 XYBTMP2 = 0 DO I = 1, NODES IF( IFOUND(I) .GT. 0) THEN J = J+1 XYBTMP2(J,:) = XYBTMP1(I,:) VERTEX(I) = J END IF END DO ! ! Number of nodes after clean up ! NX = J ! DO I = 1, NTRI I1 = TRIGPTMP1(I,1) I2 = TRIGPTMP1(I,2) I3 = TRIGPTMP1(I,3) TRIGPTMP2(I,1) = VERTEX(I1) TRIGPTMP2(I,2) = VERTEX(I2) TRIGPTMP2(I,3) = VERTEX(I3) END DO ! DEALLOCATE(XYBTMP1, IFOUND,TRIGPTMP1) DEALLOCATE(VERTEX) ! !count points connections to allocate array in W3DIMUG ! CALL COUNT(TRIGPTMP2) CALL W3DIMUG ( 1, NTRI, NX, COUNTOT, NNZ, NDSE, NDST ) ! ! fills arrays ! DO I = 1, NX XYB(I,1) = XYBTMP2(I,1) XYB(I,2) = XYBTMP2(I,2) XYB(I,3) = XYBTMP2(I,3) END DO ! DO I=1, NTRI ITMP = TRIGPTMP2(I,:) TRIGP(I,:) = ITMP END DO ! DEALLOCATE(TRIGPTMP2,XYBTMP2) ! ! call the various routines which define the point spotting strategy ! CALL SPATIAL_GRID CALL NVECTRI CALL COORDMAX ! !READMSH is only called by ww3_grid, thus the grid index is always 1. ! CALL AREA_SI(1) ! CLOSE(NDS) END SUBROUTINE READMSH SUBROUTINE GET_BOUNDARY_STATUS(STATUS) !/ !/ +-----------------------------------+ !/ | WAVEWATCH III NOAA/NCEP | !/ | | !/ | Mathieu Dutour-Sikiric (IRB) | !/ | Aron Roland (BGS IT&E GmbH) | !/ | | !/ | FORTRAN 90 | !/ | Last update : 01-Mai-2018 | !/ +-----------------------------------+ !/ !/ 01-Mai-2018 : Origination. ( version 6.04 ) !/ ! 1. Purpose : boundary status (code duplication) ! 2. Method : ! 3. Parameters : ! ! Parameter list ! ---------------------------------------------------------------- ! ---------------------------------------------------------------- ! ! 4. Subroutines used : ! ! Name Type Module Description ! ---------------------------------------------------------------- ! STRACE Subr. W3SERVMD Subroutine tracing. ! ---------------------------------------------------------------- ! ! 5. Called by : ! ! Name Type Module Description ! ---------------------------------------------------------------- ! ---------------------------------------------------------------- ! ! 6. Error messages : ! 7. Remarks ! 8. Structure : ! 9. Switches : ! ! !/S Enable subroutine tracing. ! ! 10. Source code : ! !/ ------------------------------------------------------------------- / ! USE W3GDATMD, ONLY : TRIGP, NTRI, NX IMPLICIT NONE !/ !/ ------------------------------------------------------------------- / !/ Parameter list !/ !/ ------------------------------------------------------------------- / !/ Local PARAMETERs !/ !/ !/ ------------------------------------------------------------------- / !/ ! integer, intent(out) :: STATUS(NX) INTEGER :: COLLECTED(NX), NEXTVERT(NX), PREVVERT(NX) INTEGER :: ISFINISHED, INEXT, IPREV INTEGER :: IPNEXT, IPPREV, ZNEXT, IP, I, IE STATUS(:) = 0 DO IE=1,NTRI DO I=1,3 IF (I.EQ.1) THEN IPREV=3 ELSE IPREV=I-1 END IF IF (I.EQ.3) THEN INEXT=1 ELSE INEXT=I+1 END IF IP=TRIGP(IE,I) IPNEXT=TRIGP(IE,INEXT) IPPREV=TRIGP(IE,IPREV) IF (STATUS(IP).EQ.0) THEN STATUS(IP)=1 PREVVERT(IP)=IPPREV NEXTVERT(IP)=IPNEXT END IF END DO END DO STATUS(:)=0 DO COLLECTED(:)=0 DO IE=1,NTRI DO I=1,3 IF (I.EQ.1) THEN IPREV=3 ELSE IPREV=I-1 END IF IF (I.EQ.3) THEN INEXT=1 ELSE INEXT=I+1 END IF IP=TRIGP(IE,I) IPNEXT=TRIGP(IE,INEXT) IPPREV=TRIGP(IE,IPREV) IF (STATUS(IP).eq.0) THEN ZNEXT=NEXTVERT(IP) IF (ZNEXT.eq.IPPREV) THEN COLLECTED(IP)=1 NEXTVERT(IP)=IPNEXT IF (NEXTVERT(IP).eq.PREVVERT(IP)) THEN STATUS(IP)=1 END IF END IF END IF END DO END DO ISFINISHED=1 DO IP=1,NX IF ((COLLECTED(IP).eq.0).and.(STATUS(IP).eq.0)) THEN STATUS(IP)=-1 END IF IF (STATUS(IP).eq.0) THEN ISFINISHED=0 END IF END DO IF (ISFINISHED.eq.1) THEN EXIT END IF END DO END SUBROUTINE !/ -------------------------------------------------------------------/ SUBROUTINE READMSHOBC(NDS, FNAME, TMPSTA, UGOBCOK) !/ ------------------------------------------------------------------- !/ +-----------------------------------+ !/ | WAVEWATCH III NOAA/NCEP | !/ | F. Ardhuin | !/ | FORTRAN 90 | !/ | Last update : 14-Mar-2018| !/ +-----------------------------------+ !/ !/ 14-Mar-2018 : Origination. ( version 6.02 ) !/ ! ! 1. Purpose : ! ! Reads open boundary information for UNST grids following GMESH type format ! ! 2. Method : ! ! Reads an ASCII file ! ! 3. Parameters : ! ! Parameter list ! ---------------------------------------------------------------- ! NDS Int. I Data set number used for search. ! FNAME Char*60 I File name ! TMPSTA Char*60 I/O status map to be updated (for OBC, TMPSTA = 2) ! UGOBCOK Logical O flag for proper reading of OBC file ! ---------------------------------------------------------------- ! ! 4. Subroutines used : NONE ! ! 5. Called by : ! Name Type Module Description ! ---------------------------------------------------------------- ! W3GRID Prog. Model configuration program ! ---------------------------------------------------------------- ! ! 6. Error messages : ! ! 7. Remarks : ! ! 8. Structure : ! ! 9. Switches : ! ! 10. Source code : ! !/ ------------------------------------------------------------------- / USE W3GDATMD, ONLY: NX, NY, CCON , COUNTCON USE W3ODATMD, ONLY: NDSE, NDST, NDSO USE W3SERVMD, ONLY: ITRACE, NEXTLN, EXTCDE IMPLICIT NONE !/ ------------------------------------------------------------------- / !/ Parameter list !/ INTEGER, INTENT(IN) :: NDS CHARACTER(60), INTENT(IN) :: FNAME INTEGER, INTENT(INOUT) :: TMPSTA(NY,NX) LOGICAL, INTENT(OUT) :: UGOBCOK !/ !/ local parameters !/ INTEGER :: I, IERR INTEGER(KIND=4) :: Ind,ntag, INode CHARACTER :: COMSTR*1, SPACE*1 = ' ', CELS*64 REAL, ALLOCATABLE :: TAGS(:) CHARACTER(LEN=120) :: LINE UGOBCOK=.FALSE. OPEN(NDS,FILE = FNAME,STATUS='old') READ (NDS,'(A)') COMSTR IF (COMSTR.EQ.' ') COMSTR = '$' CALL NEXTLN(COMSTR, NDS, NDSE) IERR = 0 DO WHILE (IERR.EQ.0) READ (NDS,'(A100)',END=2001,ERR=2002,IOSTAT=IERR) LINE READ(LINE,*,IOSTAT=IERR) Ind,ntag IF (IERR.EQ.0) THEN ALLOCATE(TAGS(ntag)) READ(LINE,*,IOSTAT=IERR) Ind,ntag,TAGS,INODE IF (IERR.EQ.0) THEN TMPSTA(1,INODE)=2 DEALLOCATE(TAGS) ELSE GOTO 2001 END IF END IF END DO CLOSE(NDS) UGOBCOK=.TRUE. RETURN ! 2001 CONTINUE WRITE (NDSE,1001) CALL EXTCDE ( 61 ) ! 2002 CONTINUE WRITE (NDSE,1002) IERR CALL EXTCDE ( 62 ) 1001 FORMAT (/' *** WAVEWATCH III ERROR IN READMSHOBC : '/ & ' PREMATURE END OF FILE IN READING ',A/) 1002 FORMAT (/' *** WAVEWATCH III ERROR IN READMSHOBC : '/ & ' ERROR IN READING ',A,' IOSTAT =',I8/) END SUBROUTINE READMSHOBC !/ ------------------------------------------------------------------- / !/ ------------------------------------------------------------------- / SUBROUTINE UG_GETOPENBOUNDARY(TMPSTA,ZBIN,ZLIM) !/ +-----------------------------------+ !/ | WAVEWATCH III NOAA/NCEP | !/ | F. Ardhuin | !/ | FORTRAN 90 | !/ | Last update : 30-Aug-2012| !/ +-----------------------------------+ !/ !/ 30-Aug-2012 : Adpatation from SHOM-Ifremer program( version 4.07 ) !/ ! ! 1. purpose: defines open boundary points based on depth ! 2. Method : a boundary node has more node around it than triangles ! ! 3. Parameters : ! TMPSTA: status map to be updated (for OBC, TMPSTA = 2) ! ! 4. Subroutines used : ! ! 5. Called by : ! ! Name Type Module Description ! ---------------------------------------------------------------- ! w3GRID Prog. Model configuration program ! ---------------------------------------------------------------- ! ! 6. Error messages : ! ! 7. Remarks : ! ! ! 8. Structure : ! ! 9. Switches : ! ! 10. Source code : USE W3GDATMD, ONLY: NX, NY, CCON, COUNTCON, IOBP IMPLICIT NONE !/ ------------------------------------------------------------------- / !/ Parameter list !/ INTEGER, INTENT(INOUT) :: TMPSTA(NY,NX) REAL , INTENT(IN) :: ZBIN(NY,NX) REAL , INTENT(IN) :: ZLIM !/ !/ ------------------------------------------------------------------- / !/ Local parameters !/ INTEGER :: IBC, IX INTEGER :: MASK(NX), STATUS(NX) ! MASK(:)=1 CALL SET_IOBP (MASK, STATUS) ! DO IBC = 1, N_OUTSIDE_BOUNDARY IX = OUTSIDE_BOUNDARY(IBC) !write(*,*) 'TEST1', IX, TMPSTA(1,IX), CCON(IX), COUNTCON(IX), ZBIN(1,IX), ZLIM IF (IX.NE.0) THEN ! There was a bug in the mesh conversion, OUTSIDE_BOUNDARY(IBC) should not be zero IF ((TMPSTA(1,IX).EQ.1).AND.(STATUS(IX).EQ.0) & .AND.(ZBIN(1,IX).LT.ZLIM)) TMPSTA(1,IX)=2 END IF END DO ! END SUBROUTINE UG_GETOPENBOUNDARY !/ ------------------------------------------------------------------- / SUBROUTINE SPATIAL_GRID !/ ------------------------------------------------------------------- !/ +-----------------------------------+ !/ | WAVEWATCH III NOAA/NCEP | !/ | A. Roland (BGS IT&E GbmH) | !/ | F. Ardhuin (IFREMER) | !/ | FORTRAN 90 | !/ | Last update : 31-Aug-2011| !/ +-----------------------------------+ !/ !/ 15-May-2007 : Origination: adjustment from the WWM code ( version 3.13 ) !/ 31-Aug-2011 : Simplfies the cross products ( version 4.05 ) !/ ! ! 1. Purpose : ! ! Calculates triangle areas and reorders the triangles to have them ! oriented counterclockwise ! ! 2. Method : ! ! The triangle surface calculation is based on cross product. ! ! 3. Parameters : ! ! 4. Subroutines used : ! ! 5. Called by : ! ! Name Type Module Description ! ---------------------------------------------------------------- ! READTRI Subr. Internal Unstructured mesh definition. ! ---------------------------------------------------------------- ! ! 6. Error messages : ! ! 7. Remarks : ! ! This part of code is adapted from the WWM wave model develop at the Darmstadt University ! (Aaron Roland) ! ! 8. Structure : ! ! 9. Switches : ! ! 10. Source code : ! !/ ------------------------------------------------------------------- / USE W3GDATMD IMPLICIT NONE ! !local parameters ! REAL :: TL1, TL2, TL3, TMPTRIGP INTEGER :: I1, I2, I3 INTEGER :: K !/ ------------------------------------------------------------------- / DO K = 1, NTRI I1 = TRIGP(K,1) I2 = TRIGP(K,2) I3 = TRIGP(K,3) ! ! cross product of edge-vector (orientated anticlockwise) ! TRIA(K) = REAL( (XYB(I2,2)-XYB(I1,2)) & ! (Y2-Y1) *(XYB(I1,1)-XYB(I3,1)) & ! *(X1-X3) +(XYB(I3,2)-XYB(I1,2)) & ! (Y3-Y1)*(X2-X1) *(XYB(I2,1)-XYB(I1,1)) )*0.5 ! ! test on negative triangle area, which means that the orientiation is not as assumed to be anticw. ! therefore we swap the nodes !!! ! IF (TRIA(K) .lt. TINY(1.)) THEN TMPTRIGP = TRIGP(K,2) TRIGP(K,2) = TRIGP(K,3) TRIGP(K,3) = TMPTRIGP I2 = TRIGP(K,2) I3 = TRIGP(K,3) TRIA(K) = -1.d0*TRIA(K) STOP 'WRONG TRIANGLE' END IF END DO END SUBROUTINE ! SUBROUTINE NVECTRI !/ ------------------------------------------------------------------- !/ +-----------------------------------+ !/ | WAVEWATCH III NOAA/NCEP | !/ | A. Roland | !/ | FORTRAN 90 | !/ | Last update : 15-May-2008| !/ +-----------------------------------+ !/ !/ 15-May-2007 : Origination: adjustment from the WWM code ( version 3.13 ) !/ ! ! 1. Purpose : ! ! Calculate cell tools: inward normal, angles and length of edges. ! ! 2. Method : ! To get inward pointing normals, triangle are glanced through anti-clockwisely ! ! 3. Parameters : ! ! 4. Subroutines used : ! ! 5. Called by : ! ! Name Type Module Description ! ---------------------------------------------------------------- ! READTRI Subr. Internal Unstructured mesh definition. ! ---------------------------------------------------------------- ! ! 6. Error messages : ! ! 7. Remarks : ! ! 8. Structure : ! ! 9. Switches : ! ! 10. Source code : ! !/ ------------------------------------------------------------------- / USE W3GDATMD USE CONSTANTS IMPLICIT NONE ! !local parameter ! INTEGER :: IP, IE INTEGER :: I1, I2, I3, I11, I22, I33 REAL*8 :: P1(2), P2(2), P3(2) REAL*8 :: R1(2), R2(2), R3(2) REAL*8 :: N1(2), N2(2), N3(2) REAL*8 :: TMP(3) REAL*8 :: TMPINV(3) !/ ------------------------------------------------------------------- / DO IE = 1, NTRI ! ! vertices ! I1 = TRIGP(IE,1) I2 = TRIGP(IE,2) I3 = TRIGP(IE,3) P1(1) = XYB(I1,1) P1(2) = XYB(I1,2) P2(1) = XYB(I2,1) P2(2) = XYB(I2,2) P3(1) = XYB(I3,1) P3(2) = XYB(I3,2) ! ! I1 -> I2, I2 -> I3, I3 -> I1 (anticlockwise orientation is preserved) ! R1 = P3-P2 R2 = P1-P3 R3 = P2-P1 N1(1) = (-R1(2)) N1(2) = ( R1(1)) N2(1) = (-R2(2)) N2(2) = ( R2(1)) N3(1) = (-R3(2)) N3(2) = ( R3(1)) ! ! edges length ! LEN(IE,1) = DSQRT(R1(1)**2+R1(2)**2) LEN(IE,2) = DSQRT(R2(1)**2+R2(2)**2) LEN(IE,3) = DSQRT(R3(1)**2+R3(2)**2) ! ! inward normal used for propagation (not normalized) ! IEN(IE,1) = N1(1) IEN(IE,2) = N1(2) IEN(IE,3) = N2(1) IEN(IE,4) = N2(2) IEN(IE,5) = N3(1) IEN(IE,6) = N3(2) END DO END SUBROUTINE SUBROUTINE COUNT(TRIGPTEMP) !/ ------------------------------------------------------------------- !/ +-----------------------------------+ !/ | WAVEWATCH III NOAA/NCEP | !/ | A. Roland | !/ | F. Ardhuin | !/ | FORTRAN 90 | !/ | Last update : 15-May-2008| !/ +-----------------------------------+ !/ !/ 15-May-2007 : Origination. ( version 3.13 ) !/ ! ! 1. Purpose : ! ! Calculate global and maximum number of connection for array allocations . ! ! 2. Method : ! ! 3. Parameters : ! Parameter list ! ---------------------------------------------------------------- ! NTRI Int. I Total number of triangle. ! TRIGPTEMP Int I Temporary array of triangle vertices ! COUNTRI Int O Maximum number of connected triangle ! for a given points ! COUNTOT Int O Global number of triangle connection ! for the whole grid. ! ---------------------------------------------------------------- ! 4. Subroutines used : ! ! 5. Called by : ! ! Name Type Module Description ! ---------------------------------------------------------------- ! READTRI Subr. Internal Unstructured mesh definition. ! ---------------------------------------------------------------- ! ! 6. Error messages : ! ! 7. Remarks : ! ! 8. Structure : ! ! 9. Switches : ! ! 10. Source code : ! !/ ------------------------------------------------------------------- / USE W3GDATMD IMPLICIT NONE !/ parameter list INTEGER,INTENT(IN) :: TRIGPTEMP(:,:) !/ ------------------------------------------------------------------- / !/ local parameter INTEGER :: CONN(NX) INTEGER :: COUNTER, IP, IE, I, J, N(3) COUNTRI=0 COUNTOT=0 CONN(:)= 0 ! !calculate the number of connected triangles for a given point. ! DO IE = 1,NTRI N(:) = 0. N(1) = TRIGPTEMP(IE,1) N(2) = TRIGPTEMP(IE,2) N(3) = TRIGPTEMP(IE,3) CONN(N(1)) = CONN(N(1)) + 1 CONN(N(2)) = CONN(N(2)) + 1 CONN(N(3)) = CONN(N(3)) + 1 ENDDO COUNTRI = MAXVAL(CONN) ! ! calculate the global number of connections available through the mesh ! J=0 DO IP=1,NX DO I=1,CONN(IP) J=J+1 ENDDO ENDDO COUNTOT=J END SUBROUTINE SUBROUTINE COORDMAX !/ ------------------------------------------------------------------- !/ +-----------------------------------+ !/ | WAVEWATCH III NOAA/NCEP | !/ | F. Ardhuin | !/ | FORTRAN 90 | !/ | Last update : 15-May-2008| !/ +-----------------------------------+ !/ !/ 15-May-2007 : Origination. ( version 3.13 ) !/ ! 1. Purpose : ! ! Calculate first point and last point coordinates, and minimum and maximum edge length. ! ! 2. Method : ! ! 3. Parameters : ! ! 4. Subroutines used : ! ! 5. Called by : ! ! Name Type Module Description ! ---------------------------------------------------------------- ! READTRI Subr. Internal Unstructured mesh definition. ! ---------------------------------------------------------------- ! ! 6. Error messages : ! ! 7. Remarks : ! ! 8. Structure : ! ! 9. Switches : ! ! 10. Source code : ! !/ ------------------------------------------------------------------- / USE W3GDATMD IMPLICIT NONE ! ! maximum of coordinates s ! MAXX = MAXVAL(XYB(:,1)) MAXY = MAXVAL(XYB(:,2)) ! ! minimum of coordinates ! X0 = MINVAL(XYB(:,1)) Y0 = MINVAL(XYB(:,2)) ! !maximum and minimum length of edges ! DXYMAX = MAXVAL(LEN(:,:)) SX = MINVAL(LEN(:,:)) SY = SX ! END SUBROUTINE !------------------------------------------------------------------------- SUBROUTINE AREA_SI(IMOD) !/ ------------------------------------------------------------------- !/ +-----------------------------------+ !/ | WAVEWATCH III NOAA/NCEP | !/ | A. Roland | !/ | FORTRAN 90 | !/ | Last update : 23-Aug-2011| !/ +-----------------------------------+ !/ !/ 15-May-2007 : Origination: adjustment from the WWM code ( version 3.13 ) !/ 23-Aug-2011 : Removes double entries in VNEIGH ( version 4.04 ) !/ ! ! 1. Purpose : ! ! Define optimized connection arrays (points and triangles) for spatial propagation schemes. ! ! 2. Method : ! ! 3. Parameters : ! ! 4. Subroutines used : ! ! 5. Called by : ! ! Name Type Module Description ! ---------------------------------------------------------------- ! READTRI Subr. Internal Unstructured mesh definition. ! ---------------------------------------------------------------- ! ! 6. Error messages : ! ! 7. Remarks : ! ! The storage is optimize especially considering the iterative solver used. ! The schemes used are vertex-centered, a point has to be considered within its ! median dual cell. For a given point, the surface of the dual cell is one third ! of the sum of the surface of connected triangles. ! This routine is from WWM developped in Darmstadt(Aaron Roland) ! ! 8. Structure : ! ! 9. Switches : ! ! 10. Source code : ! !/ ------------------------------------------------------------------- / USE W3GDATMD IMPLICIT NONE !/ input INTEGER, INTENT(IN) :: IMOD !/ local parameters INTEGER :: COUNTER,ifound,alreadyfound INTEGER :: I, J, K, II INTEGER :: IP, IE, POS, POS_I, POS_J, POS_K, IP_I, IP_J, IP_K INTEGER :: I1, I2, I3, IP2, CHILF(NX) INTEGER :: TMP(NX), CELLVERTEX(NX,COUNTRI,2) INTEGER :: COUNT_MAX DOUBLE PRECISION :: TRIA03 INTEGER, ALLOCATABLE :: PTABLE(:,:) !/ ------------------------------------------------------------------- / WRITE(*,'("+TRACE......",A)') 'COMPUTE SI, TRIA und CCON' SI(:) = 0.D0 ! CCON(:) = 0 ! Number of connected Elements DO IE = 1 , NTRI I1 = TRIGP(IE,1) I2 = TRIGP(IE,2) I3 = TRIGP(IE,3) CCON(I1) = CCON(I1) + 1 CCON(I2) = CCON(I2) + 1 CCON(I3) = CCON(I3) + 1 TRIA03 = 1./3. * TRIA(IE) SI(I1) = SI(I1) + TRIA03 SI(I2) = SI(I2) + TRIA03 SI(I3) = SI(I3) + TRIA03 ENDDO CELLVERTEX(:,:,:) = 0 ! Stores for each node the Elementnumbers of the connected Elements ! and the Position of the position of the Node in the Element Index WRITE(*,'("+TRACE......",A)') 'COMPUTE CELLVERTEX' CHILF = 0 DO IE = 1, NTRI DO J=1,3 I = TRIGP(IE,J)!INE(J,IE) CHILF(I) = CHILF(I)+1 CELLVERTEX(I,CHILF(I),1) = IE CELLVERTEX(I,CHILF(I),2) = J END DO ENDDO WRITE(*,'("+TRACE......",A)') 'COMPUTE IE_CELL and POS_CELL' ! ! Second step in storage, the initial 3D array CELLVERTEX, is transformed in a 1D array ! the global index is J . From now, all the computation step based on these arrays must ! abide by the conservation of the 2 loop algorithm (points + connected triangles) ! AR: I will change this now to pointers in order to omit fix loop structure for the LTS stuff ... ! INDEX_CELL(1)=1 J = 0 DO IP = 1, NX DO I = 1, CCON(IP) J = J + 1 IE_CELL(J) = CELLVERTEX(IP,I,1) POS_CELL(J) = CELLVERTEX(IP,I,2) END DO INDEX_CELL(IP+1)=J+1 END DO IF (.NOT. FSNIMP) RETURN J = 0 DO IP = 1, NX DO I = 1, CCON(IP) J = J + 1 END DO END DO COUNT_MAX = J ALLOCATE(PTABLE(COUNT_MAX,7)) J = 0 PTABLE(:,:) = 0. DO IP = 1, NX DO I = 1, CCON(IP) J = J + 1 IE = IE_CELL(J) POS = POS_CELL(J) I1 = TRIGP(IE,1) I2 = TRIGP(IE,2) I3 = TRIGP(IE,3) IF (POS == 1) THEN POS_J = 2 POS_K = 3 ELSE IF (POS == 2) THEN POS_J = 3 POS_K = 1 ELSE POS_J = 1 POS_K = 2 END IF IP_I = IP IP_J = TRIGP(IE,POS_J) IP_K = TRIGP(IE,POS_K) PTABLE(J,1) = IP_I PTABLE(J,2) = IP_J PTABLE(J,3) = IP_K PTABLE(J,4) = POS PTABLE(J,5) = POS_J PTABLE(J,6) = POS_K PTABLE(J,7) = IE END DO END DO ! WRITE(*,'("+TRACE......",A)') 'SET UP SPARSE MATRIX POINTER ... COUNT NONZERO ENTRY' J = 0 K = 0 DO IP = 1, NX TMP(:) = 0 DO I = 1, CCON(IP) J = J + 1 IP_J = PTABLE(J,2) IP_K = PTABLE(J,3) POS = PTABLE(J,4) TMP(IP) = 1 TMP(IP_J) = 1 TMP(IP_K) = 1 END DO K = K + SUM(TMP) END DO NNZ => GRIDS(IMOD)%NNZ NNZ = K ! WRITE(*,'("+TRACE......",A)') 'SET UP SPARSE MATRIX POINTER ... SETUP POINTER' ALLOCATE (GRIDS(IMOD)%JAA(NNZ)) ALLOCATE (GRIDS(IMOD)%IAA(NX+1)) ALLOCATE (GRIDS(IMOD)%POSI(3,COUNT_MAX)) JAA => GRIDS(IMOD)%JAA IAA => GRIDS(IMOD)%IAA POSI => GRIDS(IMOD)%POSI J = 0 K = 0 IAA(1) = 1 JAA = 0 DO IP = 1, NX ! Run through all rows TMP(:)=0 DO I = 1, CCON(IP) ! Check how many entries there are ... J = J + 1 ! this is the same J index as in IE_CELL IP_J = PTABLE(J,2) IP_K = PTABLE(J,3) TMP(IP) = 1 TMP(IP_J) = 1 TMP(IP_K) = 1 END DO DO I = 1, NX ! Run through all columns IF (TMP(I) .GT. 0) THEN ! this is true only for the connected points K = K + 1 JAA(K) = I END IF END DO IAA(IP + 1) = K + 1 END DO POSI = 0 J = 0 DO IP = 1, NX DO I = 1, CCON(IP) J = J + 1 IP_J = PTABLE(J,2) IP_K = PTABLE(J,3) DO K = IAA(IP), IAA(IP+1) - 1 IF (IP == JAA(K)) POSI(1,J) = K IF (IP_J == JAA(K)) POSI(2,J) = K IF (IP_K == JAA(K)) POSI(3,J) = K IF (K == 0) THEN WRITE(*,*) 'ERROR IN AREA_SI K .EQ. 0' STOP END IF END DO END DO END DO DEALLOCATE(PTABLE) END SUBROUTINE SUBROUTINE IS_IN_UNGRID(IMOD, XTIN, YTIN, ITOUT, IS, JS, RW) !/ ------------------------------------------------------------------- !/ +-----------------------------------+ !/ | WAVEWATCH III NOAA/NCEP | !/ | Mathieu Dutour Sikiric, IRB | !/ | Aron Roland, Z&P | !/ | Fabrice Ardhuin | !/ | FORTRAN 90 | !/ | Last update : 26-Jan-2014| !/ +-----------------------------------+ !/ !/ Adapted from other subroutine !/ 15-Oct-2007 : Origination. ( version 3.13 ) !/ 21-Sep-2012 : Uses same interpolation as regular ( version 4.08 ) !/ 26-Jan-2014 : Correcting bug in RW ( version 4.18 ) !/ ! 1. Purpose : ! ! Determine whether a point is inside or outside an unstructured grid, ! and returns index of triangle and interpolation weights ! This is the analogue for triangles of the FUNCTION W3GRMP ! ! 2. Method : ! ! Using barycentric coordinates defined as the ratio of triangle algebric areas ! which are positive or negative. ! Computes the 3 interpolation weights for each triangle until they are all positive ! ! 3. Parameters : ! ! Parameter list ! ---------------------------------------------------------------- ! IMOD Int. I Model number to point to. ! XTIN Real I X-coordinate of target point. ! YTIN Real I Y-coordinate of target point. ! ITOUT Int. I Model number to point to. ! IS,JS I.A. O (I,J) indices of vertices of enclosing grid cell. ! RW R.A. O Array of interpolation weights. ! ---------------------------------------------------------------- ! ! 4. Subroutines used : ! ! None ! ! 5. Called by : ! ! WMGLOW, W3IOPP, WMIOPP, WW3_GINT ! ! 6. Error messages : ! ! - Error checks on previous setting of variable. ! ! 7. Remarks : ! ! 8. Structure : ! ! 9. Switches : ! ! !/S Enable subroutine tracing. ! !/T Enable test output ! ! 10. Source code : ! ! 2. Method : ! ! Using barycentric coordinates. Each coefficient depends on the mass of its related point in the interpolation. ! ! 3. Parameters : ! ! 4. Subroutines used : ! ! 5. Called by : ! ! Name Type Module Description ! ---------------------------------------------------------------- ! W3IOPP Subr. Internal Preprocessing of point output. ! ---------------------------------------------------------------- ! ! 6. Error messages : ! ! 7. Remarks : ! ! This subroutine is adjusted from CREST code (Fabrice Ardhuin) ! For a given output point, the algorithm enable to glance through all the triangles ! to find the one the point belong to, and then make interpolation. ! ! 8. Structure : ! ! 9. Switches : ! ! !/LLG Spherical grid. ! !/XYG Carthesian grid. ! ! 10. Source code : ! !/ ------------------------------------------------------------------- / USE W3GDATMD USE W3SERVMD, ONLY: EXTCDE USE W3ODATMD, ONLY: NDSE IMPLICIT NONE !/ ------------------------------------------------------------------- / ! Parameter list INTEGER, INTENT(IN) :: IMOD REAL , INTENT(IN) :: XTIN, YTIN INTEGER, INTENT(OUT) :: itout INTEGER, INTENT(OUT) :: IS(4), JS(4) REAL, INTENT(OUT) :: RW(4) !/ ------------------------------------------------------------------- / !local parameters DOUBLE PRECISION :: x1, x2, x3 DOUBLE PRECISION :: y1, y2, y3 DOUBLE PRECISION :: s1, s2, s3, sg1, sg2, sg3 INTEGER :: ITRI INTEGER :: I1, I2, I3 INTEGER :: nbFound ! itout = 0 nbFound=0 ITRI = 0 DO WHILE (nbFound.EQ.0.AND.ITRI.LT.GRIDS(IMOD)%NTRI) ITRI = ITRI +1 I1=GRIDS(IMOD)%TRIGP(ITRI,1) I2=GRIDS(IMOD)%TRIGP(ITRI,2) I3=GRIDS(IMOD)%TRIGP(ITRI,3) ! coordinates of the first vertex A x1=GRIDS(IMOD)%XYB(I1,1) y1=GRIDS(IMOD)%XYB(I1,2) ! coordinates of the 2nd vertex B x2=GRIDS(IMOD)%XYB(I2,1) y2=GRIDS(IMOD)%XYB(I2,2) !coordinates of the 3rd vertex C x3=GRIDS(IMOD)%XYB(I3,1) y3=GRIDS(IMOD)%XYB(I3,2) !with M = (XTIN,YTIN) the target point ... !vector product of AB and AC sg3=(y3-y1)*(x2-x1)-(x3-x1)*(y2-y1) !vector product of AB and AM s3=(YTIN-y1)*(x2-x1)-(XTIN-x1)*(y2-y1) !vector product of BC and BA sg1=(y1-y2)*(x3-x2)-(x1-x2)*(y3-y2) !vector product of BC and BM s1=(YTIN-y2)*(x3-x2)-(XTIN-x2)*(y3-y2) !vector product of CA and CB sg2=(y2-y3)*(x1-x3)-(x2-x3)*(y1-y3) !vector product of CA and CM s2=(YTIN-y3)*(x1-x3)-(XTIN-x3)*(y1-y3) IF ((s1*sg1.GE.0).AND.(s2*sg2.GE.0).AND.(s3*sg3.GE.0)) THEN itout=ITRI nbFound=nbFound+1 IS(1)=I1 IS(2)=I2 IS(3)=I3 IS(4)=1 JS(:)=1 RW(1)=s1/sg1 RW(2)=s2/sg2 RW(3)=1.-RW(1)-RW(2) !s3/sg3 RW(4)=0. END IF ENDDO END SUBROUTINE IS_IN_UNGRID !/ ------------------------------------------------------------------- SUBROUTINE IS_IN_UNGRID2(IMOD, XTIN, YTIN, FORCE, ITOUT, IS, JS, RW) !/ ------------------------------------------------------------------- !/ +-----------------------------------+ !/ | WAVEWATCH III NOAA/NCEP | !/ | Mathieu Dutour Sikiric, IRB | !/ | Aron Roland, Z&P | !/ | Fabrice Ardhuin | !/ | FORTRAN 90 | !/ | Last update : 26-Jan-2014| !/ +-----------------------------------+ !/ !/ Adapted from other subroutine !/ 15-Oct-2007 : Origination. ( version 3.13 ) !/ 21-Sep-2012 : Uses same interpolation as regular ( version 4.08 ) !/ 26-Jan-2014 : Correcting bug in RW ( version 4.18 ) !/ ! 1. Purpose : ! ! Determine whether a point is inside or outside an unstructured grid, ! and returns index of triangle and interpolation weights ! This is the analogue for triangles of the FUNCTION W3GRMP ! ! 2. Method : ! ! Using barycentric coordinates defined as the ratio of triangle algebric areas ! which are positive or negative. ! Computes the 3 interpolation weights for each triangle until they are all positive ! ! 3. Parameters : ! ! Parameter list ! ---------------------------------------------------------------- ! IMOD Int. I Model number to point to. ! XTIN Real I X-coordinate of target point. ! YTIN Real I Y-coordinate of target point. ! ITOUT Int. I Model number to point to. ! IS,JS I.A. O (I,J) indices of vertices of enclosing grid cell. ! RW R.A. O Array of interpolation weights. ! ---------------------------------------------------------------- ! ! 4. Subroutines used : ! ! None ! ! 5. Called by : ! ! WMGLOW, W3IOPP, WMIOPP, WW3_GINT ! ! 6. Error messages : ! ! - Error checks on previous setting of variable. ! ! 7. Remarks : ! ! 8. Structure : ! ! 9. Switches : ! ! !/S Enable subroutine tracing. ! !/T Enable test output ! ! 10. Source code : ! ! 2. Method : ! ! Using barycentric coordinates. Each coefficient depends on the mass of its related point in the interpolation. ! ! 3. Parameters : ! ! 4. Subroutines used : ! ! 5. Called by : ! ! Name Type Module Description ! ---------------------------------------------------------------- ! W3IOPP Subr. Internal Preprocessing of point output. ! ---------------------------------------------------------------- ! ! 6. Error messages : ! ! 7. Remarks : ! ! This subroutine is adjusted from CREST code (Fabrice Ardhuin) ! For a given output point, the algorithm enable to glance through all the triangles ! to find the one the point belong to, and then make interpolation. ! ! 8. Structure : ! ! 9. Switches : ! ! !/LLG Spherical grid. ! !/XYG Carthesian grid. ! ! 10. Source code : ! !/ ------------------------------------------------------------------- / USE W3GDATMD USE W3SERVMD, ONLY: EXTCDE USE W3ODATMD, ONLY: NDSE IMPLICIT NONE !/ ------------------------------------------------------------------- / ! Parameter list INTEGER, INTENT(IN) :: IMOD, FORCE REAL , INTENT(IN) :: XTIN, YTIN INTEGER, INTENT(OUT) :: itout INTEGER, INTENT(OUT) :: IS(4), JS(4) REAL, INTENT(OUT) :: RW(4) !/ ------------------------------------------------------------------- / !local parameters DOUBLE PRECISION :: x1, x2, x3, D1, D2, D3, DISTMIN, DDMIN DOUBLE PRECISION :: s1, s2, s3, sg1, sg2, sg3, smin, ssum DOUBLE PRECISION :: y1, y2, y3 INTEGER :: ITRI, ITRIS INTEGER :: I1, I2, I3 INTEGER :: nbFound LOGICAL :: MAPSTAOK ! itout = 0 nbFound=0 ITRI = 0 ITRIS = 1 ssum = 0 smin = 0 DO WHILE (nbFound.EQ.0.AND.ITRI.LT.GRIDS(IMOD)%NTRI) ITRI = ITRI +1 I1=GRIDS(IMOD)%TRIGP(ITRI,1) I2=GRIDS(IMOD)%TRIGP(ITRI,2) I3=GRIDS(IMOD)%TRIGP(ITRI,3) ! coordinates of the first vertex A x1=GRIDS(IMOD)%XYB(I1,1) y1=GRIDS(IMOD)%XYB(I1,2) ! coordinates of the 2nd vertex B x2=GRIDS(IMOD)%XYB(I2,1) y2=GRIDS(IMOD)%XYB(I2,2) !coordinates of the 3rd vertex C x3=GRIDS(IMOD)%XYB(I3,1) y3=GRIDS(IMOD)%XYB(I3,2) !with M = (XTIN,YTIN) the target point ... !vector product of AB and AC sg3=(y3-y1)*(x2-x1)-(x3-x1)*(y2-y1) !vector product of AB and AM s3=(YTIN-y1)*(x2-x1)-(XTIN-x1)*(y2-y1) !vector product of BC and BA sg1=(y1-y2)*(x3-x2)-(x1-x2)*(y3-y2) !vector product of BC and BM s1=(YTIN-y2)*(x3-x2)-(XTIN-x2)*(y3-y2) !vector product of CA and CB sg2=(y2-y3)*(x1-x3)-(x2-x3)*(y1-y3) !vector product of CA and CM s2=(YTIN-y3)*(x1-x3)-(XTIN-x3)*(y1-y3) ! ssum = ABS(s1*sg1)+ABS(s2*sg2)+ABS(s3*sg3) MAPSTAOK = ((GRIDS(IMOD)%MAPSTA(1,I1).GE.1).AND. & (GRIDS(IMOD)%MAPSTA(1,I2).GE.1).AND.(GRIDS(IMOD)%MAPSTA(1,I3).GE.1)) IF (FORCE.LT.2) MAPSTAOK =.TRUE. ssum = (XTIN-(x1+x2+x3)/3.)**2+(YTIN-(y1+y2+y2)/3.)**2 IF (smin.EQ.0.AND. MAPSTAOK ) smin=ssum !WRITE(6,*) 'ssum',ITRI,MAPSTAOK,ssum,smin IF (ssum.LT.smin .AND. MAPSTAOK ) THEN smin=ssum ITRIS=ITRI ENDIF IF ((s1*sg1.GE.0).AND.(s2*sg2.GE.0).AND.(s3*sg3.GE.0)) THEN itout=ITRI nbFound=nbFound+1 IS(1)=I1 IS(2)=I2 IS(3)=I3 IS(4)=1 JS(:)=1 RW(1)=s1/sg1 RW(2)=s2/sg2 RW(3)=1.-RW(1)-RW(2) !s3/sg3 RW(4)=0. END IF ENDDO IF (itout.EQ.0.AND.FORCE.GT.0) THEN ITRI=ITRIS I1=GRIDS(IMOD)%TRIGP(ITRI,1) I2=GRIDS(IMOD)%TRIGP(ITRI,2) I3=GRIDS(IMOD)%TRIGP(ITRI,3) ! coordinates of the first vertex A x1=GRIDS(IMOD)%XYB(I1,1) y1=GRIDS(IMOD)%XYB(I1,2) ! coordinates of the 2nd vertex B x2=GRIDS(IMOD)%XYB(I2,1) y2=GRIDS(IMOD)%XYB(I2,2) !coordinates of the 3rd vertex C x3=GRIDS(IMOD)%XYB(I3,1) y3=GRIDS(IMOD)%XYB(I3,2) D1=(XTIN-X1)**2+(YTIN-Y1)**2 D2=(XTIN-X2)**2+(YTIN-Y2)**2 D3=(XTIN-X3)**2+(YTIN-Y3)**2 IF (D1.LE.D2.AND.D1.LE.D3) IS(1)=I1 IF (D2.LE.D1.AND.D2.LE.D3) IS(1)=I2 IF (D3.LE.D2.AND.D3.LE.D1) IS(1)=I3 IS(2:4)=1 JS(:)=1 RW(1)=1 RW(2:4)=0. ITOUT=ITRI ENDIF END SUBROUTINE IS_IN_UNGRID2 !/ ------------------------------------------------------------------- / SUBROUTINE UG_GRADIENTS (PARAM, DIFFX, DIFFY) !/ +-----------------------------------+ !/ | WAVEWATCH III NOAA/NCEP | !/ | F. Ardhuin | !/ | A. Roland | !/ | FORTRAN 90 | !/ | Last update : 14-Oct-2013| !/ +-----------------------------------+ !/ !/ 15-Nov-2007 : Origination. ( version 3.13 ) !/ 31-Oct-2010 : Merging of 4.03 with 3.14-Ifremer ( version 4.04 ) !/ 08-Nov-2011 : Correction for zero grad. on contour( version 4.04 ) !/ 14-Oct-2013 : Correction of latitude factor ( version 4.12 ) !/ 01-Mai-2018 : Using linear shape function for gradients [ version 6.04) !/ ! ! 1. purpose: calculate gradients at a point via its connection. ! 2. Method : using linear shape function this is a basis on which ! all advection schemes in Roland (2008) are checked. ! ! 3. Parameters : ! PARAM : depth or current field (indices 0 to NSEA) ! DIFFX : x gradient (indices 1 to NX) ! DIFFY : y gradient (indices 1 to NX) ! ! 4. Subroutines used : ! ! 5. Called by : ! ! Name Type Module Description ! ---------------------------------------------------------------- ! W3WAVE Subr. Actual wind wave routine ! ---------------------------------------------------------------- ! ! 6. Error messages : ! ! 7. Remarks : ! ! This subroutine is adjusted from WWM code (Aaron Roland) ! ! 8. Structure : ! ! 9. Switches : ! ! 10. Source code : USE CONSTANTS USE W3GDATMD, ONLY : TRIGP, NTRI, NX, NSEA, MAPFS, CLATIS, & MAPSTA, ANGLE, FLAGLL, IOBP, IEN, TRIA, NSEAL, NTRI USE W3ADATMD, ONLY : NSEALM IMPLICIT NONE REAL, INTENT(IN) :: PARAM(0:NSEA) REAL, INTENT(OUT) :: DIFFX(:,:), DIFFY(:,:) ! local parameters INTEGER :: VERTICES(3), NI(3), NI_GL(3) REAL :: TMP1(3), TMP2(3) INTEGER :: I, IX, IE, IE_GL REAL :: VAR(3), FACT, LATMEAN REAL :: DIFFXTMP, DIFFYTMP REAL :: DEDX(3), DEDY(3) REAL :: DVDXIE, DVDYIE REAL :: WEI(NX), WEI_LOCAL(NSEAL) DIFFX = 0. DIFFY = 0. ! IF (FLAGLL) THEN FACT=1./(DERA*RADIUS) ELSE FACT=1. END IF WEI = 0. DO IE = 1, NTRI NI = TRIGP(IE,:) LATMEAN = 1./3. * SUM(CLATIS(MAPFS(1,NI))) WEI(NI) = WEI(NI) + 2.*TRIA(IE) IF (TRIA(IE) .LT. TINY(1.)) THEN WRITE(*,*) 'AREA SMALLER ZERO', IE, NTRI, TRIA(IE) STOP ENDIF DEDX(1) = IEN(IE,1) DEDX(2) = IEN(IE,3) DEDX(3) = IEN(IE,5) DEDY(1) = IEN(IE,2) DEDY(2) = IEN(IE,4) DEDY(3) = IEN(IE,6) VAR = PARAM(MAPFS(1,NI)) * FACT DVDXIE = DOT_PRODUCT( VAR,DEDX) DVDYIE = DOT_PRODUCT( VAR,DEDY) DIFFX(1,NI) = DIFFX(1,NI) + DVDXIE * LATMEAN DIFFY(1,NI) = DIFFY(1,NI) + DVDYIE END DO DIFFX(1,:) = DIFFX(1,:)/WEI DIFFY(1,:) = DIFFY(1,:)/WEI ! END SUBROUTINE UG_GRADIENTS !/ ------------------------------------------------------------------- / SUBROUTINE W3NESTUG(DISTMIN,FLOK) !/ !/ +-----------------------------------+ !/ | WAVEWATCH III NOAA/NCEP | !/ | | !/ | Aron Roland (BGS IT&E GmbH) | !/ | Mathieu Dutour-Sikiric (IRB) | !/ | | !/ | FORTRAN 90 | !/ | Last update : 01-June-2018 | !/ +-----------------------------------+ !/ !/ 01-June-2018 : Origination. ( version 6.04 ) !/ ! 1. Purpose : UGTYPE nesting initialization ! 2. Method : ! 3. Parameters : ! ! Parameter list ! ---------------------------------------------------------------- ! ---------------------------------------------------------------- ! ! 4. Subroutines used : ! ! Name Type Module Description ! ---------------------------------------------------------------- ! STRACE Subr. W3SERVMD Subroutine tracing. ! ---------------------------------------------------------------- ! ! 5. Called by : ! ! Name Type Module Description ! ---------------------------------------------------------------- ! ---------------------------------------------------------------- ! ! 6. Error messages : ! 7. Remarks ! 8. Structure : ! 9. Switches : ! ! !/S Enable subroutine tracing. ! ! 10. Source code : ! !/ ------------------------------------------------------------------- / ! USE W3ODATMD, ONLY: NBI, NDSE, ISBPI, XBPI, YBPI USE W3GDATMD, ONLY: NX, XYB, XGRD, YGRD, MAPSTA, MAPFS, MAPSF REAL, INTENT(IN) :: DISTMIN LOGICAL, INTENT(INOUT) :: FLOK INTEGER :: I, J, JMEMO, IS, IX, N, IX1(NBI) REAL :: DIST, DIST0 ! N = 0 ! !1. look for input boundary point index ! warning: if land points are included as boundary points to abide by the nest ! file, their status should be -2. ! IX1 = 0 ISBPI = 1 DO IX = 1, NX IF (ABS(MAPSTA (1,IX)) .EQ. 2) THEN N = N + 1 IF (N.GT.NBI) THEN WRITE(NDSE,*) 'Error: boundary node index > NBI ... nest.ww3 file is not consistent with mod_def.ww3' STOP ENDIF IX1(N) = IX END IF END DO ! !2. Matches the model grid points (where MAPSTA = 2) with the points in nest.ww3 ! For this, we use the nearest point in the nest file. ! DO I = 1, NBI !FA: This will not work with FLAGLL=.F. (XY grid) DIST0 = 360**2 IS=1 DO J = 1, N DIST=(XBPI(I)-XYB(IX1(J),1))**2+(YBPI(I)-XYB(IX1(J),2))**2 IF (DIST.LT.DIST0) THEN IS = MAPFS(1,IX1(J)) DIST0=DIST JMEMO=J END IF END DO DIST0=SQRT(DIST0) IF (DIST0.LE.DISTMIN) THEN ISBPI(I)=IS ELSE FLOK=.TRUE. END IF END DO IF ( N .NE. NBI) THEN WRITE(NDSE ,900) N, NBI DO J=1,N WRITE(6,*) 'THIS POINT HAS MAPSTA=2:',ISBPI(J) END DO ISBPI(N+1:NBI)=ISBPI(1) END IF 900 FORMAT (/' *** WAVEWATCH III ERROR IN W3IOBC : '/ & ' NUMBER OF MAPSTA=2 DIFFERS FROM NUMBER IN nest.ww3 '/ & ' CHECK nest.ww3 AND ww3_grid.inp ',2I8/) END SUBROUTINE !/ ------------------------------------------------------------------- / SUBROUTINE SET_IOBP (MASK, STATUS) !/ !/ +-----------------------------------+ !/ | WAVEWATCH III NOAA/NCEP | !/ | | !/ | Aron Roland (BGS IT&E GmbH) | !/ | Mathieu Dutour-Sikiric (IRB) | !/ | | !/ | FORTRAN 90 | !/ | Last update : 01-June-2018 | !/ +-----------------------------------+ !/ !/ 01-June-2018 : Origination. ( version 6.04 ) !/ ! 1. Purpose : setup boundary pointer ! 2. Method : ! 3. Parameters : ! ! Parameter list ! ---------------------------------------------------------------- ! ---------------------------------------------------------------- ! ! 4. Subroutines used : ! ! Name Type Module Description ! ---------------------------------------------------------------- ! STRACE Subr. W3SERVMD Subroutine tracing. ! ---------------------------------------------------------------- ! ! 5. Called by : ! ! Name Type Module Description ! ---------------------------------------------------------------- ! ---------------------------------------------------------------- ! ! 6. Error messages : ! 7. Remarks ! 8. Structure : ! 9. Switches : ! ! !/S Enable subroutine tracing. ! ! 10. Source code : ! !/ ------------------------------------------------------------------- / ! !/ ! USE CONSTANTS ! USE W3GDATMD, ONLY: NX, NTRI, TRIGP USE W3ODATMD, ONLY: IAPROC IMPLICIT NONE !/ !/ ------------------------------------------------------------------- / !/ Parameter list !/ INTEGER, INTENT(IN) :: MASK(NX) INTEGER, INTENT(OUT) :: STATUS(NX) ! INTEGER :: COLLECTED(NX), NEXTVERT(NX), PREVVERT(NX) INTEGER :: ISFINISHED !, INEXT, IPREV INTEGER :: INEXT(3), IPREV(3) INTEGER :: ZNEXT, IP, I, IE, IPNEXT, IPPREV, COUNT integer nb0, nb1, nbM1 STATUS = -1 INEXT=(/ 2, 3, 1 /) !IPREV=1+MOD(I+1,3) IPREV=(/ 3, 1, 2 /) !INEXT=1+MOD(I,3) DO IE=1,NTRI ! If one of the points of the triangle is masked out (land) then do as if triangle does not exist... ! IF ((MASK(TRIGP(IE,1)).GT.0).AND.(MASK(TRIGP(IE,2)).GT.0).AND.(MASK(TRIGP(IE,3)).GT.0)) THEN DO I=1,3 IP=TRIGP(IE,I) CALL TRIANG_INDEXES(I, IPNEXT, IPPREV) !IPNEXT=TRIGP(IE,INEXT(I)) !IPPREV=TRIGP(IE,IPREV(I)) IF (STATUS(IP).EQ.-1) THEN STATUS(IP)=1 PREVVERT(IP)=IPPREV NEXTVERT(IP)=IPNEXT END IF END DO ! ENDIF END DO STATUS(:)=-1 ! COUNT = 0 DO COUNT = COUNT + 1 COLLECTED(:)=0 DO IE=1,NTRI ! IF ((MASK(TRIGP(IE,1)).GT.0).AND.(MASK(TRIGP(IE,2)).GT.0).AND.(MASK(TRIGP(IE,3)).GT.0)) THEN DO I=1,3 IP=TRIGP(IE,I) CALL TRIANG_INDEXES(I, IPNEXT, IPPREV) !IPNEXT=TRIGP(IE,INEXT(I)) !IPPREV=TRIGP(IE,IPREV(I)) IF (STATUS(IP).EQ.-1) THEN ZNEXT=NEXTVERT(IP) IF (ZNEXT.EQ.IPPREV) THEN COLLECTED(IP)=1 NEXTVERT(IP)=IPNEXT IF (NEXTVERT(IP).EQ.PREVVERT(IP)) THEN STATUS(IP)=1 END IF END IF END IF END DO ! END IF ! end of test on MASK END DO ! ! Checks that all nodes have been treated ... ! ISFINISHED=1 DO IP=1,NX IF (MASK(IP).LE.0) THEN STATUS(IP)=0 ELSE IF ((COLLECTED(IP).EQ.0).AND.(STATUS(IP).EQ.-1)) THEN STATUS(IP)=0 END IF IF (STATUS(IP).eq.-1) THEN ISFINISHED=0 END IF ENDIF END DO IF (ISFINISHED.EQ.1) THEN EXIT END IF END DO STATUS = 1 CALL GET_BOUNDARY(NX, NTRI, TRIGP, STATUS, PREVVERT, NEXTVERT) ! DO IP= 1, NX ! WRITE(12000,*) IP, STATUS(IP) ! ENDDO !#ifdef MPI_PARALL_GRID ! CALL exchange_p2di(STATUS) !#endif END SUBROUTINE SET_IOBP !/ ------------------------------------------------------------------- / SUBROUTINE GET_BOUNDARY(MNP, MNE, TRIGP, IOBP, NEIGHBOR_PREV, & & NEIGHBOR_NEXT) !/ !/ +-----------------------------------+ !/ | WAVEWATCH III NOAA/NCEP | !/ | | !/ | Aron Roland (BGS IT&E GmbH) | !/ | Mathieu Dutour-Sikiric (IRB) | !/ | | !/ | FORTRAN 90 | !/ | Last update : 01-June-2018 | !/ +-----------------------------------+ !/ !/ 01-June-2018 : Origination. ( version 6.04 ) !/ ! 1. Purpose : find boundary points ! 2. Method : ! 3. Parameters : ! ! Parameter list ! ---------------------------------------------------------------- ! ---------------------------------------------------------------- ! ! 4. Subroutines used : ! ! Name Type Module Description ! ---------------------------------------------------------------- ! STRACE Subr. W3SERVMD Subroutine tracing. ! ---------------------------------------------------------------- ! ! 5. Called by : ! ! Name Type Module Description ! ---------------------------------------------------------------- ! ---------------------------------------------------------------- ! ! 6. Error messages : ! 7. Remarks ! 8. Structure : ! 9. Switches : ! ! !/S Enable subroutine tracing. ! ! 10. Source code : ! !/ ------------------------------------------------------------------- / ! USE W3SERVMD, ONLY: EXTCDE IMPLICIT NONE !/ !/ ------------------------------------------------------------------- / !/ Parameter list !/ !/ ------------------------------------------------------------------- / !/ Local PARAMETERs !/ INTEGER, INTENT(IN) :: MNP, MNE, TRIGP(MNE,3) INTEGER, INTENT(INOUT) :: IOBP(MNP) INTEGER, INTENT(INOUT) :: NEIGHBOR_PREV(MNP) INTEGER, INTENT(INOUT) :: NEIGHBOR_NEXT(MNP) INTEGER, POINTER :: STATUS(:) INTEGER, POINTER :: COLLECTED(:) INTEGER, POINTER :: NEXTVERT(:) INTEGER, POINTER :: PREVVERT(:) INTEGER :: IE, I, IP, IP2, IP3 INTEGER :: ISFINISHED, INEXT, IPREV, ISTAT INTEGER :: IPNEXT, IPPREV, ZNEXT, ZPREV ALLOCATE(STATUS(MNP), STAT=ISTAT) CHECK_ALLOC_STATUS ( ISTAT ) ALLOCATE(COLLECTED(MNP), STAT=ISTAT) CHECK_ALLOC_STATUS ( ISTAT ) ALLOCATE(PREVVERT(MNP), STAT=ISTAT) CHECK_ALLOC_STATUS ( ISTAT ) ALLOCATE(NEXTVERT(MNP), STAT=ISTAT) CHECK_ALLOC_STATUS ( ISTAT ) NEIGHBOR_NEXT = 0 NEIGHBOR_PREV = 0 ! Now computing the next items STATUS = 0 NEXTVERT = 0 PREVVERT = 0 DO IE=1,MNE DO I=1,3 CALL TRIANG_INDEXES(I, INEXT, IPREV) IP=TRIGP(IE,I) IPNEXT=TRIGP(IE,INEXT) IPPREV=TRIGP(IE,IPREV) IF (STATUS(IP).EQ.0) THEN STATUS(IP)=1 PREVVERT(IP)=IPPREV NEXTVERT(IP)=IPNEXT END IF END DO END DO STATUS(:)=0 DO COLLECTED(:)=0 DO IE=1,MNE DO I=1,3 CALL TRIANG_INDEXES(I, INEXT, IPREV) IP=TRIGP(IE,I) IPNEXT=TRIGP(IE,INEXT) IPPREV=TRIGP(IE,IPREV) IF (STATUS(IP).EQ.0) THEN ZNEXT=NEXTVERT(IP) IF (ZNEXT.EQ.IPPREV) THEN COLLECTED(IP)=1 NEXTVERT(IP)=IPNEXT IF (NEXTVERT(IP).EQ.PREVVERT(IP)) THEN STATUS(IP)=1 END IF END IF END IF END DO END DO ISFINISHED=1 DO IP=1,MNP IF ((COLLECTED(IP).EQ.0).AND.(STATUS(IP).EQ.0)) THEN STATUS(IP)=-1 NEIGHBOR_NEXT(IP)=NEXTVERT(IP) END IF IF (STATUS(IP).EQ.0) THEN ISFINISHED=0 END IF END DO IF (ISFINISHED.EQ.1) THEN EXIT END IF END DO ! Now computing the prev items STATUS = 0 NEXTVERT = 0 PREVVERT = 0 DO IE=1,MNE DO I=1,3 CALL TRIANG_INDEXES(I, INEXT, IPREV) IP=TRIGP(IE,I) IPNEXT=TRIGP(IE,INEXT) IPPREV=TRIGP(IE,IPREV) IF (STATUS(IP).EQ.0) THEN STATUS(IP)=1 PREVVERT(IP)=IPPREV NEXTVERT(IP)=IPNEXT END IF END DO END DO STATUS(:)=0 DO COLLECTED(:)=0 DO IE=1,MNE DO I=1,3 CALL TRIANG_INDEXES(I, INEXT, IPREV) IP=TRIGP(IE,I) IPNEXT=TRIGP(IE,INEXT) IPPREV=TRIGP(IE,IPREV) IF (STATUS(IP).EQ.0) THEN ZPREV=PREVVERT(IP) IF (ZPREV.EQ.IPNEXT) THEN COLLECTED(IP)=1 PREVVERT(IP)=IPPREV IF (PREVVERT(IP).EQ.NEXTVERT(IP)) THEN STATUS(IP)=1 END IF END IF END IF END DO END DO ISFINISHED=1 DO IP=1,MNP IF ((COLLECTED(IP).EQ.0).AND.(STATUS(IP).EQ.0)) THEN STATUS(IP)=-1 NEIGHBOR_PREV(IP)=PREVVERT(IP) ! new code END IF IF (STATUS(IP).EQ.0) THEN ISFINISHED=0 END IF END DO IF (ISFINISHED.EQ.1) THEN EXIT END IF END DO ! Now making checks DO IP=1,MNP IP2=NEIGHBOR_NEXT(IP) IF (IP2.GT.0) THEN IP3=NEIGHBOR_PREV(IP2) IF (ABS(IP3 - IP).GT.0) THEN WRITE(*,*) 'IP=', IP, ' IP2=', IP2, ' IP3=', IP3 WRITE(*,*) 'We have a dramatic inconsistency' STOP END IF END IF END DO ! Now assigning the boundary IOBP array DO IP=1,MNP IF (STATUS(IP).EQ.-1 .AND. IOBP(IP) .EQ. 1) THEN IOBP(IP)=0 END IF END DO DEALLOCATE(STATUS, STAT=ISTAT) CHECK_DEALLOC_STATUS ( ISTAT ) DEALLOCATE(COLLECTED, STAT=ISTAT) CHECK_DEALLOC_STATUS ( ISTAT ) DEALLOCATE(NEXTVERT, STAT=ISTAT) CHECK_DEALLOC_STATUS ( ISTAT ) DEALLOCATE(PREVVERT, STAT=ISTAT) CHECK_DEALLOC_STATUS ( ISTAT ) END SUBROUTINE !/ ------------------------------------------------------------------- / SUBROUTINE TRIANG_INDEXES(I, INEXT, IPREV) !/ !/ +-----------------------------------+ !/ | WAVEWATCH III NOAA/NCEP | !/ | | !/ | Aron Roland (BGS IT&E GmbH) | !/ | Mathieu Dutour-Sikiric (IRB) | !/ | | !/ | FORTRAN 90 | !/ | Last update : 01-June-2018 | !/ +-----------------------------------+ !/ !/ 01-June-2018 : Origination. ( version 6.04 ) !/ ! 1. Purpose : set indices of the triangle ! 2. Method : ! 3. Parameters : ! ! Parameter list ! ---------------------------------------------------------------- ! ---------------------------------------------------------------- ! ! 4. Subroutines used : ! ! Name Type Module Description ! ---------------------------------------------------------------- ! STRACE Subr. W3SERVMD Subroutine tracing. ! ---------------------------------------------------------------- ! ! 5. Called by : ! ! Name Type Module Description ! ---------------------------------------------------------------- ! ---------------------------------------------------------------- ! ! 6. Error messages : ! 7. Remarks ! 8. Structure : ! 9. Switches : ! ! !/S Enable subroutine tracing. ! ! 10. Source code : ! !/ ------------------------------------------------------------------- / IMPLICIT NONE ! !/ ------------------------------------------------------------------- / !/ Parameter list !/ !/ ------------------------------------------------------------------- / !/ Local PARAMETERs !/ !/ !/ ------------------------------------------------------------------- / !/ INTEGER, INTENT(IN) :: I INTEGER, INTENT(OUT) :: INEXT, IPREV IF (I.EQ.1) THEN INEXT=3 ELSE INEXT=I-1 END IF IF (I.EQ.3) THEN IPREV=1 ELSE IPREV=I+1 END IF END SUBROUTINE !/ ------------------------------------------------------------------- / SUBROUTINE GET_INTERFACE() !/ !/ +-----------------------------------+ !/ | WAVEWATCH III NOAA/NCEP | !/ | | !/ | Aron Roland (BGS IT&E GmbH) | !/ | Mathieu Dutour-Sikiric (IRB) | !/ | | !/ | FORTRAN 90 | !/ | Last update : 01-June-2018 | !/ +-----------------------------------+ !/ !/ 01-June-2018 : Origination. ( version 6.04 ) !/ ! 1. Purpose : get interface nodes of the wetting and drying part ! 2. Method : ! 3. Parameters : ! ! Parameter list ! ---------------------------------------------------------------- ! ---------------------------------------------------------------- ! ! 4. Subroutines used : ! ! Name Type Module Description ! ---------------------------------------------------------------- ! STRACE Subr. W3SERVMD Subroutine tracing. ! ---------------------------------------------------------------- ! ! 5. Called by : ! ! Name Type Module Description ! ---------------------------------------------------------------- ! ---------------------------------------------------------------- ! ! 6. Error messages : ! 7. Remarks ! 8. Structure : ! 9. Switches : ! ! !/S Enable subroutine tracing. ! ! 10. Source code : ! !/ ------------------------------------------------------------------- / ! USE CONSTANTS, ONLY : LPDLIB USE W3GDATMD, ONLY : NX, IOBP, CCON, NSEAL, IOBDP, IE_CELL, IOBDP, TRIGP IMPLICIT NONE !/ !/ ------------------------------------------------------------------- / !/ Parameter list !/ !/ ------------------------------------------------------------------- / !/ Local PARAMETERs !/ !/ !/ ------------------------------------------------------------------- / !/ INTEGER :: I, J, IP, IE J = 0 DO IP = 1, NSEAL IF (IOBP(IP) .NE. 1 .OR. IOBDP(IP) .EQ. 0) CYCLE DO I = 1, CCON(IP) J = J + 1 IE = IE_CELL(J) IF (ANY(IOBDP(TRIGP(IE,:)) .EQ. 0)) THEN IOBDP(IP) = -1 ! Set this node as a wet node adjacent to a dry one ... now what's next? Here on this points we want to c EXIT ENDIF ENDDO ENDDO END SUBROUTINE !/ ------------------------------------------------------------------- / SUBROUTINE SETUGIOBP ( ) !/ !/ +-----------------------------------+ !/ | WAVEWATCH III NOAA/NCEP | !/ | Fabrice Ardhuin | !/ | Aron Roland | !/ | FORTRAN 90 | !/ | Last update : 17-Apr-2016 | !/ +-----------------------------------+ !/ !/ 23-Aug-2011 : Origination. ( version 4.04 ) !/ 17-Apr-2016 : Uses optimized boundary detection ( version 5.10 ) !/ ! 1. Purpose : ! ! Redefines the values of the boundary points and angle pointers ! based on the MAPSTA array ! ! 2. Method : ! ! Adapted boundary detection from A. Roland and M. Dutour (WWM code) ! ! 3. Parameters : ! ! Parameter list ! ---------------------------------------------------------------- ! ---------------------------------------------------------------- ! ! Local variables. ! ---------------------------------------------------------------- ! ---------------------------------------------------------------- ! ! 4. Subroutines used : ! ! 5. Called by : ! ! Name Type Module Description ! ---------------------------------------------------------------- ! WW3_GRID Prog. WW3_GRID Grid preprocessor ! W3ULEV Subr. W3UPDTMD Water level update ! ---------------------------------------------------------------- ! ! 6. Error messages : ! ! None. ! ! 7. Remarks : ! ! 8. Structure : ! ! 9. Switches : ! ! !/S Enable subroutine tracing. ! ! 10. Source code : !/ ------------------------------------------------------------------- / !/ ! USE CONSTANTS ! USE W3GDATMD, ONLY: NX, NY, NSEA, MAPFS, & NK, NTH, DTH, XFR, MAPSTA, COUNTRI, & ECOS, ESIN, IEN, NTRI, TRIGP, & IOBP,IOBPD, IOBPA, & REFPARS, REFLC, REFLD, & ANGLE0, ANGLE USE W3ODATMD, ONLY: TBPI0, TBPIN, FLBPI USE W3ADATMD, ONLY: CG, CX, CY, ATRNX, ATRNY, ITIME, CFLXYMAX USE W3IDATMD, ONLY: FLCUR USE W3ODATMD, only : IAPROC IMPLICIT NONE !/ ------------------------------------------------------------------- / !/ Parameter list !/ !/ !/ ------------------------------------------------------------------- / !/ Local parameters !/ INTEGER :: ITH, IX, I, J, IP, IE, NDIRSUM REAL (KIND = 8) :: COSSUM, SINSUM REAL (KIND = 8) :: DIRMIN, DIRMAX, SHIFT, TEMPO, DIRCOAST REAL (KIND = 8) :: X1, X2, Y1, Y2, DXP1, DXP2, DXP3 REAL (KIND = 8) :: DYP1, DYP2, DYP3, eDet1, eDet2, EVX, EVY REAL(KIND=8), PARAMETER :: THR = TINY(1.) INTEGER :: I1, I2, I3 INTEGER :: ITMP(NX), NEXTVERT(NX), PREVVERT(NX) !/ ------------------------------------------------------------------- / ! ! 1. Preparations --------------------------------------------------- * ! 1.a Set constants ! ! !--- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! 2. Searches for boundary points ! ITMP = MAPSTA(1,:) CALL SET_IOBP(ITMP, IOBP) ! !--- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! 3. Defines directions pointing into land or sea ! IOBPD(:,:) = 0 IOBPA(:) = 0 ! DO IP=1,NX IF ((MAPSTA(1,IP).EQ.2).AND.(IOBP(IP).EQ.0)) IOBPA(IP)=1 !WRITE(600,*) IP,ITMP(IP),IOBP(IP),IOBPA(IP) END DO DO IE=1,NTRI I1 = TRIGP(IE,1) I2 = TRIGP(IE,2) I3 = TRIGP(IE,3) DXP1 = IEN(IE,6) DYP1 = - IEN(IE,5) DXP2 = IEN(IE,2) DYP2 = - IEN(IE,1) DXP3 = IEN(IE,4) DYP3 = - IEN(IE,3) DO ITH=1,NTH EVX=ECOS(ITH) EVY=ESIN(ITH) DO I=1,3 IF (I.eq.1) THEN x1= DXP1 y1= DYP1 x2= - DXP3 y2= - DYP3 IP= I1 END IF IF (I.eq.2) THEN x1 = DXP2 y1 = DYP2 x2 = - DXP1 y2 = - DYP1 IP = I2 END IF IF (I.eq.3) THEN x1 = DXP3 y1 = DYP3 x2 = - DXP2 y2 = - DYP2 IP = I3 END IF IF (IOBP(IP) .eq. 0) THEN ! physical boundary eDet1 = THR-x1*EVY+y1*EVX eDet2 = THR+x2*EVY-y2*EVX IF ((eDet1.gt.0.).and.(eDet2.gt.0.)) THEN ! this is the case of waves going towards the boundary ... IOBPD(ITH,IP)=1 ENDIF ELSE ! water ... IOBPD(ITH,IP)=1 END IF END DO END DO END DO DO IP = 1, NX IF ( IOBPA(IP) .eq. 1 .OR. IOBP(IP) .eq. 3 .OR. IOBP(IP) .eq. 4) IOBPD(:,IP) = 1 END DO !2do: recode for mpi ! IF (LBCWA .OR. LBCSP) THEN ! IF (.NOT. ANY(IOBP .EQ. 2)) THEN ! CALL WWM_ABORT('YOU IMPOSED BOUNDARY CONDITIONS BUT IN THE BOUNDARY FILE ARE NO NODES WITH FLAG = 2') ! ENDIF ! ENDIF !#ifdef MPI_PARALL_GRID ! CALL exchange_p2di(IOBWB) ! DO ID = 1, MDC ! iwild = IOBPD(ID,:) ! CALL exchange_p2di(iwild) ! IOBPD(ID,:) = iwild ! ENDDO !#endif !--- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! 3. Updates the reflection direction and sharp / flat shoreline angle ! ! Finds the shoreline direction from IOBPD ! REFLC(1,:)= 0. REFLD(:,:)= 1 DO IP=1,NX IF (IOBP(IP).EQ.0.AND.MAPSTA(1,IP).EQ.1) THEN COSSUM=0. SINSUM=0. NDIRSUM=0. DO ITH=1,NTH COSSUM=COSSUM+IOBPD(ITH,IP)*ECOS(ITH) SINSUM=SINSUM+IOBPD(ITH,IP)*ESIN(ITH) NDIRSUM=NDIRSUM+IOBPD(ITH,IP) END DO DIRCOAST=ATAN2(SINSUM, COSSUM) REFLD(1,MAPFS(1,IP)) = 1+MOD(NTH+NINT(DIRCOAST/DTH),NTH) REFLD(2,MAPFS(1,IP)) = 4-MAX(2,NINT(4.*REAL(NDIRSUM)/REAL(NTH))) REFLC(1,MAPFS(1,IP))= REFPARS(1) END IF END DO !DO IX=1,NX !DO ITH=1,NTH ! WRITE(500+IAPROC,*) IX,ITH,IOBP(IX),IOBPA(IX),IOBPD(ITH,IX) !,REFLD(1:2,MAPFS(1,IX)) !ENDDO !ENDDO ! ! Recomputes the angles used in the gradients estimation ! RETURN END SUBROUTINE SETUGIOBP !/ ------------------------------------------------------------------- / END MODULE W3TRIAMD