WW3/w3triamd_8F90_source.html

#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:

  !                     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

  !      USE W3ARRYMD, ONLY: INA2R, INA2I

  !      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(:)


CONTAINS

  !/ -------------------------------------------------------------------/


  SUBROUTINE readmsh(NDS,FNAME)

    !/ -------------------------------------------------------------------

    !/                  +-----------------------------------+

    !/                  | WAVEWATCH III           NOAA/NCEP |

    !/                  |           F. Ardhuin              |

    !/                  |           A. Roland               |

    !/                  |                        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:

    !                     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, ONLY: zb, xgrd, ygrd, ntri, nx, countot, trigp, nnz, w3dimug

    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.

#ifdef W3_PDLIB

    lpdlib = .true.

#endif

    !

    ! read number of nodes and nodes from Gmsh files

    !

    READ(nds,*) nodes

    ALLOCATE(xybtmp1(3,nodes))

    DO i= 1, nodes

      READ(nds,*) j, xybtmp1(1,i), xybtmp1(2,i), xybtmp1(3,i)

    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(3,nelts))

    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(1:3,j) = 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(1,k)

      i2 = trigptmp1(2,k)

      i3 = trigptmp1(3,k)


      ifound(i1)= ifound(i1) + 1

      ifound(i2)= ifound(i2) + 1

      ifound(i3)= ifound(i3) + 1

    END DO


    j = 0


    ALLOCATE(trigptmp2(3,ntri),vertex(nodes),xybtmp2(3,nodes))

    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(1,i)

      i2 = trigptmp1(2,i)

      i3 = trigptmp1(3,i)

      trigptmp2(1,i) = vertex(i1)

      trigptmp2(2,i) = vertex(i2)

      trigptmp2(3,i) = 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

      xgrd(1,i) = xybtmp2(1,i)

      ygrd(1,i) = xybtmp2(2,i)

      zb(i)     = xybtmp2(3,i)

    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


#ifdef W3_PDLIB

    IF(.false.) THEN

#endif

      CALL area_si(1)

#ifdef W3_PDLIB

    ENDIF

#endif

    !

    CLOSE(nds)

  END SUBROUTINE readmsh

  !/--------------------------------------------------------------------/


  SUBROUTINE readmsh_iobp(NDS,FNAME)

    !/ -------------------------------------------------------------------

    !/                  +-----------------------------------+

    !/                  | WAVEWATCH III           NOAA/NCEP |

    !/                  |           A. Roland               |

    !/                  |           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:

    !                     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

    LOGICAL                            :: lfile_exists

    CHARACTER                          :: COMSTR*1, SPACE*1 = ' ', cels*64

    DOUBLE PRECISION, ALLOCATABLE      :: XYBTMP1(:,:)


    INQUIRE(file=fname, exist=lfile_exists)

    IF (.NOT. lfile_exists) RETURN

    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)

    !

    ! read number of nodes and nodes from Gmsh files

    !

    READ(nds,*) nodes

    ALLOCATE(xybtmp1(3,nodes))

    DO i= 1, nodes

      READ(nds,*) j, xybtmp1(1,i), xybtmp1(2,i), xybtmp1(3,i)

      IF (int(xybtmp1(3,i)) .EQ. 3) iobp(i) = 3

    END DO

    !

    CLOSE(nds)

  END SUBROUTINE readmsh_iobp

  !/--------------------------------------------------------------------/


  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 :

    !

    !/ ------------------------------------------------------------------- /

#ifdef W3_S

    USE w3servmd, ONLY: strace

#endif

    !


#ifdef W3_PDLIB

    use yowelementpool, only: ne_global

    use yownodepool, only: np_global

#endif

    USE w3gdatmd, ONLY : trigp, ntri, nx

    IMPLICIT NONE

    !/

    !/ ------------------------------------------------------------------- /

    !/ Parameter list

    !/

    !/ ------------------------------------------------------------------- /

    !/ Local PARAMETERs

    !/

#ifdef W3_S

    INTEGER, SAVE           :: IENT = 0

#endif

    !/

    !/ ------------------------------------------------------------------- /

    !/

    !

    integer*2, intent(out) :: STATUS(NX)

    INTEGER :: COLLECTED(NX), NEXTVERT(NX), PREVVERT(NX)

    INTEGER :: ISFINISHED, INEXT, IPREV

    INTEGER :: IPNEXT, IPPREV, ZNEXT, IP, I, IE

#ifdef W3_S

    CALL strace (ient, 'VA_SETUP_IOBPD')

#endif

    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(i,ie)

        ipnext=trigp(inext,ie)

        ipprev=trigp(iprev,ie)

        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(i,ie)

          ipnext=trigp(inext,ie)

          ipprev=trigp(iprev,ie)

          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 get_boundary_status


  !/ -------------------------------------------------------------------/


  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

#ifdef W3_S

    USE w3servmd, ONLY: strace

#endif


    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


#ifdef W3_S

    USE w3servmd, ONLY: strace

#endif


    IMPLICIT NONE

    !/ ------------------------------------------------------------------- /

    !/ Parameter list

    !/

    INTEGER, INTENT(INOUT)  :: TMPSTA(NY,NX)

#ifdef W3_S

    INTEGER, SAVE           :: IENT = 0

#endif

    REAL   , INTENT(IN)     :: ZBIN(NY,NX)

    REAL   , INTENT(IN)     :: ZLIM

    !/

    !/ ------------------------------------------------------------------- /

    !/ Local parameters

    !/

    INTEGER                 ::  IBC, IX

    INTEGER                 ::  MASK(NX)

    INTEGER*2               ::  STATUS(NX)

    !

    mask(:)=1

    CALL set_iobp (mask, status)

    !

#ifdef W3_S

    CALL strace (ient, 'UG_GETOPENBOUNDARY')

#endif

    DO ibc = 1, n_outside_boundary

      ix = outside_boundary(ibc)

      !write(*,*) 'TEST1', IX, TMPSTA(1,IX), CCON(IX), COUNTCON(IX), ZBIN(1,IX), ZLIM

      ! OUTSIDE_BOUNDARY(IBC) is defined over the full nodes NODES indexes

      ! whereas TMPSTA and ZBIN are defined over the clean up list of nodes NX

      IF ((ix.NE.0).AND.(ix.LE.nx)) THEN

        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

#ifdef W3_S

    USE w3servmd, ONLY: strace

#endif

    USE w3odatmd, ONLY: ndse


    IMPLICIT NONE

    !

    !local parameters

    !

    REAL              :: TL1, TL2, TL3, TMPTRIGP

    INTEGER           :: I1, I2, I3

    INTEGER           :: K

    real*8            :: pt(3,2)


#ifdef W3_S

    INTEGER                      ::  IENT = 0

#endif

    !/ ------------------------------------------------------------------- /

#ifdef W3_S

    CALL strace (ient, 'SPATIAL_GRID')

#endif


    DO k = 1, ntri


      i1 = trigp(1,k)

      i2 = trigp(2,k)

      i3 = trigp(3,k)


!AR: todo call this only for global grid

      CALL fix_periodcity(i1,i2,i3,xgrd,ygrd,pt)

      !

      ! cross product of edge-vector  (orientated anticlockwise)

      !

      tria(k) = real( (pt(2,2)-pt(1,2)) &

           *(pt(1,1)-pt(3,1))      &

           +(pt(3,2)-pt(1,2))      &

           *(pt(2,1)-pt(1,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(2,k)

        trigp(2,k) = trigp(3,k)

        trigp(3,k) = tmptrigp

        i2 = trigp(2,k)

        i3 = trigp(3,k)

        tria(k) = -1.d0*tria(k)

      END IF

    END DO

  END SUBROUTINE spatial_grid

  !/--------------------------------------------------------------------/

  !

  !/--------------------------------------------------------------------/


  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

#ifdef W3_S

    USE w3servmd, ONLY: strace

#endif

    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)

    real*8    :: pt(3,2)

#ifdef W3_S

    INTEGER                      ::  IENT = 0

#endif

    !/ ------------------------------------------------------------------- /

#ifdef W3_S

    CALL strace (ient, 'NVECTRI')

#endif

    !

    DO ie = 1, ntri

      !

      ! vertices

      !

      i1 = trigp(1,ie)

      i2 = trigp(2,ie)

      i3 = trigp(3,ie)


      CALL fix_periodcity(i1,i2,i3,xgrd,ygrd,pt)


      p1(1) = pt(1,1)

      p1(2) = pt(1,2)

      p2(1) = pt(2,1)

      p2(2) = pt(2,2)

      p3(1) = pt(3,1)

      p3(2) = pt(3,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 nvectri

  !/---------------------------------------------------------------------------


  !/------------------------------------------------------------------------


  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

#ifdef W3_S

    USE w3servmd, ONLY: strace

#endif

    IMPLICIT NONE


    !/ parameter list


    INTEGER,INTENT(IN) :: TRIGPTEMP(:,:)

    !/ ------------------------------------------------------------------- /

    !/ local parameter


    INTEGER               :: CONN(NX)

    INTEGER               :: COUNTER, IP, IE, I, J, N(3)

#ifdef W3_S

    INTEGER                      ::  IENT = 0

#endif

    !/------------------------------------------------------------------------


#ifdef W3_S

    CALL strace (ient, 'COUNT')

#endif


    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(1,ie)

      n(2) = trigptemp(2,ie)

      n(3) = trigptemp(3,ie)

      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 count


  !/----------------------------------------------------------------------------


  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

#ifdef W3_S

    USE w3servmd, ONLY: strace

#endif

    IMPLICIT NONE

#ifdef W3_S

    INTEGER                      ::  IENT = 0

#endif


#ifdef W3_S

    CALL strace (ient, 'COORDMAX')

#endif

    !

    ! maximum of coordinates s

    !

    maxx = maxval(xgrd(1,:))

    maxy = maxval(ygrd(1,:))

    !

    ! minimum of coordinates

    !

    x0 = minval(xgrd(1,:))

    y0 = minval(ygrd(1,:))

    !

    !maximum and minimum length of edges

    !

    dxymax = maxval(len(:,:))

    sx = minval(len(:,:))

    sy = sx

    !

  END SUBROUTINE coordmax

  !-------------------------------------------------------------------------


  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

#ifdef W3_S

    USE w3servmd, ONLY: strace

#endif

    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(:,:)


#ifdef W3_S

    INTEGER                      ::  IENT = 0

#endif

    !/ ------------------------------------------------------------------- /


#ifdef W3_S

    CALL strace (ient, 'AREA_SI')

#endif


    si(:) = 0.d0

    !

    ccon(:) = 0     ! Number of connected Elements

    DO ie = 1, ntri

      i1 = trigp(1,ie)

      i2 = trigp(2,ie)

      i3 = trigp(3,ie)

      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 Node in the Element Index


    chilf = 0


    DO ie = 1, ntri

      DO j=1,3

        i = trigp(j,ie)!INE(J,IE)

        chilf(i) = chilf(i)+1

        cellvertex(i,chilf(i),1) = ie

        cellvertex(i,chilf(i),2) = j

      END DO

    ENDDO

    !

    ! 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)

    !

    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(1,ie)

        i2 = trigp(2,ie)

        i3 = trigp(3,ie)

        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(pos_j,ie)

        ip_k = trigp(pos_k,ie)

        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(*,*) .EQ.'ERROR IN AREA_SI K  0'

            stop

          END IF

        END DO

      END DO

    END DO


    DEALLOCATE(ptable)


  END SUBROUTINE area_si


  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

#ifdef W3_S

    USE w3servmd, ONLY: strace

#endif

    USE w3odatmd, ONLY: ndse

    IMPLICIT NONE


    !/ ------------------------------------------------------------------- /

    ! Parameter list


    INTEGER, INTENT(IN)            :: IMOD

    DOUBLE PRECISION, 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

    real*8                       :: pt(3,2)

    INTEGER                      :: ITRI

    INTEGER                      :: I1, I2, I3

    INTEGER                      :: nbFound

#ifdef W3_S

    INTEGER                      ::  IENT = 0

    CALL strace (ient, 'IS_IN_UNGRID')

#endif


    !

    itout = 0

    nbfound=0

    itri = 0

    DO WHILE (nbfound.EQ.0.AND.itri.LT.grids(imod)%NTRI)

      itri = itri +1

      i1=grids(imod)%TRIGP(1,itri)

      i2=grids(imod)%TRIGP(2,itri)

      i3=grids(imod)%TRIGP(3,itri)


      CALL fix_periodcity(i1,i2,i3,grids(imod)%XGRD,grids(imod)%YGRD,pt)

      ! coordinates of the first vertex A

      x1 = pt(1,1)

      y1 = pt(1,2)

      ! coordinates of the 2nd vertex B

      x2 = pt(2,1)

      y2 = pt(2,2)

      !coordinates of the 3rd vertex C

      x3 = pt(3,1)

      y3 = pt(3,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

#ifdef W3_S

    USE w3servmd, ONLY: strace

#endif

    USE w3odatmd, ONLY: ndse

    IMPLICIT NONE


    !/ ------------------------------------------------------------------- /

    ! Parameter list


    INTEGER, INTENT(IN)            :: IMOD, FORCE

    DOUBLE PRECISION, 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

#ifdef W3_S

    INTEGER                      ::  IENT = 0

    CALL strace (ient, 'IS_IN_UNGRID2')

#endif


    !

    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(1,itri)

      i2=grids(imod)%TRIGP(2,itri)

      i3=grids(imod)%TRIGP(3,itri)

      ! coordinates of the first vertex A

      x1=grids(imod)%XGRD(1,i1)

      y1=grids(imod)%YGRD(1,i1)

      ! coordinates of the 2nd vertex B

      x2=grids(imod)%XGRD(1,i2)

      y2=grids(imod)%XGRD(1,i2)

      !coordinates of the 3rd vertex C

      x3=grids(imod)%XGRD(1,i3)

      y3=grids(imod)%YGRD(1,i3)

      !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(1,itri)

      i2=grids(imod)%TRIGP(2,itri)

      i3=grids(imod)%TRIGP(3,itri)

      ! coordinates of the first vertex A

      x1=grids(imod)%XGRD(1,i1)

      y1=grids(imod)%YGRD(1,i1)

      ! coordinates of the 2nd vertex B

      x2=grids(imod)%XGRD(1,i2)

      y2=grids(imod)%YGRD(1,i2)

      !coordinates of the 3rd vertex C

      x3=grids(imod)%XGRD(1,i3)

      y3=grids(imod)%YGRD(1,i3)

      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

#ifdef W3_PDLIB

    USE yowelementpool

    use yownodepool,    only: pdlib_ien, pdlib_tria, npa

    USE yowexchangemodule, only : pdlib_exchange1dreal

#endif


    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)

    real*8               :: rtmp(nseal)


    diffx = 0.

    diffy = 0.

    !

    IF (flagll) THEN

      fact=1./(dera*radius)

    ELSE

      fact=1.

    END IF


#ifdef W3_PDLIB

    IF (.NOT. lpdlib) THEN

#endif

      wei = 0.

      DO ie = 1, ntri

        ni = trigp(:,ie)

        latmean = 1./3. * sum(clatis(mapfs(1,ni)))

        wei(ni) = wei(ni) + 2.*tria(ie)

        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

#ifdef W3_PDLIB

    ELSE

      wei_local = 0.

      DO ie = 1, ne

        ni      = ine(:,ie)

        ie_gl   = ielg(ie)

        ni_gl   = trigp(:,ie_gl)

        latmean = 1./3. * sum(clatis(mapfs(1,ni_gl)))

        wei_local(ni) = wei_local(ni) + 2.*pdlib_tria(ie)

        dedx(1)     = pdlib_ien(1,ie)

        dedx(2)     = pdlib_ien(3,ie)

        dedx(3)     = pdlib_ien(5,ie)

        dedy(1)     = pdlib_ien(2,ie)

        dedy(2)     = pdlib_ien(4,ie)

        dedy(3)     = pdlib_ien(6,ie)

        var         = param(mapfs(1,ni_gl)) * 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_local

      diffy(1,:) = diffy(1,:)/wei_local

    ENDIF

    CALL pdlib_exchange1dreal(diffx(1,:))

    CALL pdlib_exchange1dreal(diffy(1,:))

#endif

    !

  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 :

    !

    !/ ------------------------------------------------------------------- /

#ifdef W3_S

    USE w3servmd, ONLY: strace

#endif

    !

    USE w3odatmd, ONLY: nbi, ndse, isbpi, xbpi, ybpi

    USE w3gdatmd, ONLY: nx, 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

#ifdef W3_T

        WRITE(ndse ,*)'ADDING BOUNDARY POINT:',n,ix

#endif

      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

      dist0 = huge(1.)

      is = 1

      DO j = 1, n

        dist = (xbpi(i) - xgrd(1,ix1(j)))**2 + (ybpi(i) - ygrd(1,ix1(j)))**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

#ifdef W3_T

        WRITE(ndse ,'(A,I6,A,I7,A,I6)') 'MATCHED BOUNDARY POINT:',i,'GRID POINT:', &

             mapsf(is,1),'INDEX IN nest.ww3:', jmemo

#endif

      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 w3nestug


  !/ ------------------------------------------------------------------- /


  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 :

    !

    !/ ------------------------------------------------------------------- /

#ifdef W3_S

    USE w3servmd, ONLY: strace

#endif

    !

    !/

    !

    USE constants

    !

    !

    USE w3gdatmd, ONLY: nx, ntri, trigp

    USE w3odatmd, ONLY: iaproc


    IMPLICIT NONE


    !/

    !/ ------------------------------------------------------------------- /

    !/ Parameter list

    !/

    INTEGER, INTENT(IN)   :: MASK(NX)

    INTEGER*2, 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(1,IE)).GT.0).AND.(MASK(TRIGP(2,IE)).GT.0).AND.(MASK(TRIGP(3,IE)).GT.0)) THEN

      DO i=1,3

        ip=trigp(i,ie)

        CALL triang_indexes(i, ipnext, ipprev)

        !IPNEXT=TRIGP(INEXT(I),IE)

        !IPPREV=TRIGP(IPREV(I),IE)

        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(1,IE)).GT.0).AND.(MASK(TRIGP(2,IE)).GT.0).AND.(MASK(TRIGP(3,IE)).GT.0)) THEN

        DO i=1,3

          ip=trigp(i,ie)

          CALL triang_indexes(i, ipnext, ipprev)

          !IPNEXT=TRIGP(INEXT(I),IE)

          !IPPREV=TRIGP(IPREV(I),IE)

          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)


    !#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 :

    !

    !/ ------------------------------------------------------------------- /

#ifdef W3_S

    USE w3servmd, ONLY: strace

#endif

    !

    USE w3servmd, ONLY: extcde

    IMPLICIT NONE

    !/

    !/ ------------------------------------------------------------------- /

    !/ Parameter list

    !/

    !/ ------------------------------------------------------------------- /

    !/ Local PARAMETERs

    !/

#ifdef W3_S

    INTEGER, SAVE           :: IENT = 0

#endif


    INTEGER, INTENT(IN)             :: MNP, MNE, TRIGP(3,MNE)

    INTEGER*2, 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

#ifdef W3_S

    CALL strace (ient, 'GET_BOUNDARY')

#endif

    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(i,ie)

        ipnext=trigp(inext,ie)

        ipprev=trigp(iprev,ie)

        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(i,ie)

          ipnext=trigp(inext,ie)

          ipprev=trigp(iprev,ie)

          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(i,ie)

        ipnext=trigp(inext,ie)

        ipprev=trigp(iprev,ie)

        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(i,ie)

          ipnext=trigp(inext,ie)

          ipprev=trigp(iprev,ie)

          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 get_boundary


  !/ ------------------------------------------------------------------- /


  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 :

    !

    !/ ------------------------------------------------------------------- /

#ifdef W3_S

    USE w3servmd, ONLY: strace

#endif

    IMPLICIT NONE

    !

    !/ ------------------------------------------------------------------- /

    !/ Parameter list

    !/

    !/ ------------------------------------------------------------------- /

    !/ Local PARAMETERs

    !/

#ifdef W3_S

    INTEGER, SAVE           :: IENT = 0

#endif

    !/

    !/ ------------------------------------------------------------------- /

    !/

    INTEGER, INTENT(IN)  :: I

    INTEGER, INTENT(OUT) :: INEXT, IPREV

#ifdef W3_S

    CALL strace (ient, 'TRIANG_INDEXES')

#endif

    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 triang_indexes


  !/ ------------------------------------------------------------------- /


  SUBROUTINE set_ug_iobp()

    !/

    !/                  +-----------------------------------+

    !/                  | 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,                          &

#ifdef w3_ref1

         refpars, reflc, refld,                      &

#endif

         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

#ifdef W3_S

    USE w3servmd, ONLY: strace

#endif


    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)

    CHARACTER(60) :: FNAME

#ifdef W3_S

    INTEGER, SAVE           :: IENT = 0

#endif

    !/ ------------------------------------------------------------------- /

    !

    ! 1.  Preparations --------------------------------------------------- *

    ! 1.a Set constants

    !

#ifdef W3_S

    CALL strace (ient, 'SETUGIOBP')

#endif

    !

    !--- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    ! 2.  Searches for boundary points

    !

    itmp = mapsta(1,:)

    CALL set_iobp(itmp, iobp)

    fname = 'meshbnd.msh'

    CALL readmsh_iobp(23456,fname)

    !

    !--- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    ! 3. Defines directions pointing into land or sea

    !

    iobpd(:,:) = 0

    iobpa(:)   = 0

    !

    DO ip=1,nx

      IF (mapsta(1,ip).EQ.2) THEN

        iobpa(ip) = 1

        iobp(ip)  = 2

      ENDIF

    END DO


    DO ie = 1,ntri

      i1   =   trigp(1,ie)

      i2   =   trigp(2,ie)

      i3   =   trigp(3,ie)

      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


#ifdef W3_REF1

    !

    ! 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

#endif

    !

    ! Recomputes the angles used in the gradients estimation

    !

    !

    RETURN

  END SUBROUTINE set_ug_iobp

  !/ ------------------------------------------------------------------- /


  SUBROUTINE fix_periodcity(I1,I2,I3,XGRD,YGRD,PT)

    !/

    !/                  +-----------------------------------+

    !/                  | WAVEWATCH III           NOAA/NCEP |

    !/                  |        Steven Brus                |

    !/                  |        Ali Abdolali               |

    !/                  |                        FORTRAN 90 |

    !/                  | Last update :         21-May-2020 |

    !/                  +-----------------------------------+

    !/

    !/    21-May-2020 : Origination.                        ( version 6.07 )

    !/

    !/

    !  1. Purpose :

    !

    !     Adjust element longitude coordinates for elements straddling the

    !     dateline with distance of ~360 degrees

    !

    !  2. Method :

    !

    !     Detect if element has nodes on both sides of dateline and adjust

    !     coordinates so that all nodes have the same sign

    !

    !  3. Parameters :

    !

    !     Parameter list

    !     ----------------------------------------------------------------

    IMPLICIT NONE

    INTEGER, INTENT(IN) :: I1, I2, I3

    DOUBLE PRECISION, INTENT(IN) :: XGRD(:,:), YGRD(:,:)

    real*8, INTENT(OUT) :: pt(3,2)

    !     ----------------------------------------------------------------

    !

    !     Local variables.

    !     ----------------------------------------------------------------

    INTEGER :: I

    INTEGER :: R1GT180, R2GT180, R3GT180

    !     ----------------------------------------------------------------

    !

    !  4. Subroutines used :

    !


    !  5. Called by :

    !

    !      Name         Type  Module   Description

    !     ----------------------------------------------------------------

    !      SPATIAL_GRID Subr. W3TRIAM  Triangle area calculation

    !      NVECTRI      Subr. W3TRIAM  Edge length, angle, normal calcuation

    !      IS_IN_UNGRID Subr. W3TRIAM  Point in element calculation

    !     ----------------------------------------------------------------

    !

    !  6. Error messages :

    !

    !       None.

    !

    !  7. Remarks :

    !

    !  8. Structure :

    !

    !  9. Switches :

    !

    ! 10. Source code :

    !/ ------------------------------------------------------------------- /


    pt(1,1) = xgrd(1,i1)

    pt(1,2) = ygrd(1,i1)

    pt(2,1) = xgrd(1,i2)

    pt(2,2) = ygrd(1,i2)

    pt(3,1) = xgrd(1,i3)

    pt(3,2) = ygrd(1,i3)


    r1gt180 = merge(1, 0, abs(pt(3,1)-pt(2,1)).GT.180)

    r2gt180 = merge(1, 0, abs(pt(1,1)-pt(3,1)).GT.180)

    r3gt180 = merge(1, 0, abs(pt(2,1)-pt(1,1)).GT.180)

    ! if R1GT180+R2GT180+R3GT180 .eq. 0 the element does not cross the dateline

    ! if R1GT180+R2GT180+R3GT180 .eq. 1 the element contains the pole

    ! if R1GT180+R2GT180+R3GT180 .eq. 2 the element crosses the dateline


    IF ( r1gt180 + r2gt180 == 2 ) THEN

      pt(3,1)=pt(3,1)-sign(360.0d0,(pt(3,1)-pt(2,1)))

    ELSE IF ( r2gt180 + r3gt180 == 2 ) THEN

      pt(1,1)=pt(1,1)-sign(360.0d0,(pt(1,1)-pt(2,1)))

    ELSE IF ( r1gt180 + r3gt180 == 2 ) THEN

      pt(2,1)=pt(2,1)-sign(360.0d0,(pt(2,1)-pt(3,1)))

    ENDIF


    RETURN

  END SUBROUTINE fix_periodcity

END MODULE w3triamd