wmscrpmd.F90

Go to the documentation of this file.

 
 
#include "w3macros.h"
!/ ------------------------------------------------------------------- /
MODULE wmscrpmd
  !/
  !/                  +-----------------------------------+
  !/                  | WAVEWATCH III                     |
  !/                  |   E. Rogers, M. Dutour,           |
  !/                  |   A. Roland, F. Ardhuin           |
  !/                  |                        FORTRAN 90 |
  !/                  | Last update :         10-Dec-2014 |
  !/                  +-----------------------------------+
  !/
  !/    06-Sep-2012 : Origination, transfer from WMGRIDMD ( version 4.08 )
  !/    10-Dec-2014 : Add checks for allocate status      ( version 5.04 )
  !/
  !/    Copyright 2012 National Weather Service (NWS),
  !/       National Oceanic and Atmospheric Administration.  All rights
  !/       reserved.  WAVEWATCH III is a trademark of the NWS.
  !/       No unauthorized use without permission.
  !/
  !  1. Purpose :
  !
  !     Routines to determine and process grid dependencies in the
  !     multi-grid wave model.
  !
  !  2. Variables and types :
  !
  !  3. Subroutines and functions :
  !
  !      Name                          Type  Scope      Description
  !     ----------------------------------------------------------------
  !      scrip_wrapper                  Subr. Public   as the name says ...
  !      get_scrip_info_structured      Subr. Public   as the name says ...
  !      get_scrip_info_unstructured    Subr. Public   as the name says ...
  !      get_scrip_info                 Subr. Public   as the name says ...
  !      scrip_info_renormalization     Subr. Public   as the name says ...
  !      TRIANG_INDEXES                 Subr. Public   as the name says ...
  !      get_unstructured_vertex_degree Subr. Public   as the name says ...
  !      GET_BOUNDARY                   Subr. Public   as the name says ...
  !     ----------------------------------------------------------------
  !
  !  4. Subroutines and functions used :
  !
  !      Name                             Type  Module   Description
  !     ----------------------------------------------------------------
  !     get_unstructured_vertex_degree    Subr. W3TRIAMD Manage data structures
  !     ----------------------------------------------------------------
  !
  !  5. Remarks :
  !
  !     - The subroutines TRIANG_INDEXES, get_unstructured_vertex_degree, and
  !       GET_BOUNDARY should probably be renamed and moved to the module w3triamd
  !
  !  6. Switches :
  !
  !
  !     !/S    Enable subroutine tracing.
  !     !/Tn   Enable test output.
  !
  !  7. Source code :
  !
  !/ ------------------------------------------------------------------- /
  !/
  !/ Specify default accessibility
  !/
  PUBLIC
  !/
  !/ Module private variable for checking error returns
  !/
  INTEGER, PRIVATE        :: ISTAT
  !/
CONTAINS
  !/ ------------------------------------------------------------------- /
  SUBROUTINE scrip_wrapper (ID_SRC, ID_DST,                         &
       MAPSTA_SRC,MAPST2_SRC,FLAGLL,GRIDSHIFT,L_MASTER,    &
       L_READ,L_TEST)
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III                     |
    !/                  | E. Rogers, M. Dutour, A. Roland   |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         10-Dec-2014 !
    !/                  +-----------------------------------+
    !/
    !  1. Original author :
    !
    !     Erick Rogers, NRL
    !
    !  2. Last update :
    !
    !     See revisions.
    !
    !  3. Revisions :
    !
    !     29-Apr-2011 : Origination                         ( version 4.01 )
    !     20-Feb-2012 : Mathieu Dutour Sikiric, Aron Roland Z&P
    !                   Modification is to split the code into several
    !                   subroutines
    !                   ---get_scrip_info_structured (the structured case)
    !                   ---get_scrip_info (chooses according to FD/FE)
    !                   ---scrip_info_renormalization (conv_x and all that)
    !     11-Apr-2013   Kevin Lind
    !                   Added code for reading/writing SCRIP remap files
    !/    10-Dec-2014 : Add checks for allocate status      ( version 5.04 )
    !
    !  4. Copyright :
    !
    !  5. Purpose :
    !
    !     Compute grid information required by SCRIP
    !
    !  6. Method :
    !
    !  7. Parameters, Variables and types :
    !
    !  8. Called by :
    !
    !     Subroutine WMGHGH
    !
    !  9. Subroutines and functions used :
    !
    !     Subroutine SCRIP
    !
    ! 10. Error messages:
    !
    ! 11. Remarks :
    !
    ! 12. Structure :
    !
    ! 13. Switches :
    !
    ! 14. Source code :
 
    USE scrip_grids
    USE scrip_remap_vars
    USE scrip_constants
    USE scrip_kindsmod
    USE scrip_interface
    USE w3servmd, ONLY: extcde
    !      USE W3GDATMD, ONLY : GRIDS
 
    IMPLICIT NONE
    INTEGER(SCRIP_I4), INTENT(IN)      :: ID_SRC, ID_DST
    INTEGER(SCRIP_I4), INTENT(IN)      :: MAPSTA_SRC(:,:)
    INTEGER(SCRIP_I4), INTENT(IN)      :: MAPST2_SRC(:,:)
    LOGICAL(SCRIP_LOGICAL), INTENT(IN) :: FLAGLL
    REAL   (SCRIP_R8), INTENT(IN)      :: GRIDSHIFT
    LOGICAL(SCRIP_LOGICAL), INTENT(IN) :: L_MASTER  ! Am I the master processor (do I/O)?
    LOGICAL(SCRIP_LOGICAL), INTENT(IN) :: L_READ    ! Do I read the remap file?
    LOGICAL(SCRIP_LOGICAL), INTENT(IN) :: L_TEST    ! Whether to include test output
    ! in subroutines
    !/ ------------------------------------------------------------------- /
    !/ local variables
    !/
    INTEGER(SCRIP_I4)                 :: IREC,I,J,NI,NJ,IDUM,NK,K, &
         ILINK,IW,ICORNER,NGOODPNTS, &
         NBADPNTS
    INTEGER(SCRIP_I4)                 :: ISRC,JSRC,KSRC,IPNT,KDST, &
         NI_SRC
    REAL   (SCRIP_R8)                 :: LAT_CONVERSION,OFFSET
    REAL   (SCRIP_R8)                 :: CONV_DX,CONV_DY,WEIGHT
    REAL   (SCRIP_R8)                 :: WTSUM
#ifdef W3_T38
    CHARACTER (LEN=10)      :: CDATE_TIME(3)
    INTEGER                 :: DATE_TIME(8)
    INTEGER                 :: ELAPSED_TIME, BEG_TIME, &
         END_TIME
#endif
 
    ! test output for input variables
 
#ifdef W3_T38
    if(l_master)write(*,*)'flagll = ',flagll
    if(l_master)write(*,*)'gridshift = ',gridshift
#endif
 
    !
    ! START: universal settings
    !
 
    !     Set variables for converting to degrees
    !     notes: SCRIP only operates on spherical coordinates, so for the case
    !         where the problem is specified by the user as in a
    !         meters/cartesian coordinate system, it is necessary to make
    !         a temporary conversion to a "fake" spherical coordinate grid,
    !         to keep SCRIP happy. The good news here is that multi-grid
    !         meters-grid simulations will be very rare: we will probably only
    !         encounter them in the context of simple test cases. Strictly
    !         speaking, this conversion does not even need to be physically
    !         correct, e.g. we could say that 1000 km is 1 deg....as long as
    !         we are consistent between grids.
    !     Potential future improvement: make conv_dy and offset such that dst grid
    !         covers a specific longitude range, e.g. 1 deg east to 2 deg east
 
    !
    ! START: set up src grid
    !
 
    !notes: when we work out how to interface with an unstructured grid,
    !       we will need to revisit this issue of how to set grid1_rank, etc.
    !       strategy: declare variables in grid module, but allocate them here.
    !
    grid1_units='degrees' ! the other option is radians...we don't use this
    grid1_name='src' ! this is not used, except for netcdf output
    CALL get_scrip_info(id_src,                                       &
         &   grid1_center_lon, grid1_center_lat,                            &
         &   grid1_corner_lon, grid1_corner_lat, grid1_mask,                &
         &   grid1_dims, grid1_size, grid1_corners, grid1_rank)
    grid2_units='degrees'
    grid2_name='dst'
    CALL get_scrip_info(id_dst,                                       &
         &   grid2_center_lon, grid2_center_lat,                            &
         &   grid2_corner_lon, grid2_corner_lat, grid2_mask,                &
         &   grid2_dims, grid2_size, grid2_corners, grid2_rank)
 
    IF(flagll)THEN
      conv_dx=one
      conv_dy=one
      offset=zero
    ELSE
      lat_conversion=zero ! lat_conversion
      ! notes:    this is the latitude used for conversion everywhere
      !           in the grid (approximation)
      !           (in radians)
      !           conv_dy=92.6*1200.0 ! physical, =92.6/(3/3600)=111000 m = 111 km
      conv_dy=1.0e+6_scrip_r8 ! non-physical, 1e+6=1 deg
      conv_dx=cos(lat_conversion)*conv_dy
      ! notes:    offset (in meters), is necessary so that our grid does
      !           not lie on the branch cut
      offset=75000.0_scrip_r8-min(minval(grid1_center_lon),          &
           &                               minval(grid2_center_lon))
    ENDIF
 
    !.....test output
#ifdef W3_T38
    write(*,*)'l_master = ',l_master
    if(l_master)then
      write(*,*)'conv_dx=', conv_dx
      write(*,*)'conv_dy=', conv_dy
      write(*,*)'offset = ',offset
      write(*,*)'grid1_size=', grid1_size
      write(*,*)'grid2_size=', grid2_size
      write(*,*)'l_read = ',l_read
      write(*,*)'minval(grid1_center_lon) = ',minval(grid1_center_lon)
      write(*,*)'maxval(grid1_center_lon) = ',maxval(grid1_center_lon)
      write(*,*)'minval(grid1_center_lat) = ',minval(grid1_center_lat)
      write(*,*)'maxval(grid1_center_lat) = ',maxval(grid1_center_lat)
      write(*,*)'minval(grid2_center_lon) = ',minval(grid2_center_lon)
      write(*,*)'maxval(grid2_center_lon) = ',maxval(grid2_center_lon)
      write(*,*)'minval(grid2_center_lat) = ',minval(grid2_center_lat)
      write(*,*)'maxval(grid2_center_lat) = ',maxval(grid2_center_lat)
    endif
#endif
 
    CALL scrip_info_renormalization(                                  &
         &   grid1_center_lon, grid1_center_lat,                            &
         &   grid1_corner_lon, grid1_corner_lat, grid1_mask,                &
         &   grid1_dims, grid1_size, grid1_corners, grid1_rank,             &
         &   conv_dx, conv_dy, offset, gridshift)
    CALL scrip_info_renormalization(                                  &
         &   grid2_center_lon, grid2_center_lat,                            &
         &   grid2_corner_lon, grid2_corner_lat, grid2_mask,                &
         &   grid2_dims, grid2_size, grid2_corners, grid2_rank,             &
         &   conv_dx, conv_dy, offset, zero)
 
    !.....Set constants for thresholding weights:
    frac_lowest =1.e-3_scrip_r8
    frac_highest=one+1.e-3_scrip_r8
    wt_lowest   =zero
    wt_highest  =one+1.e-7_scrip_r8
 
#ifdef W3_T38
    call date_and_time (cdate_time(1), cdate_time(2), cdate_time(3), date_time)
    beg_time = ((date_time(5)*60 + date_time(6))*60 +date_time(7))*1000 + date_time(8)
#endif
    CALL scrip(id_src, id_dst, l_master, l_read, l_test)
#ifdef W3_T38
    call date_and_time (cdate_time(1), cdate_time(2), cdate_time(3), date_time)
    end_time = ((date_time(5)*60 + date_time(6))*60 +date_time(7))*1000 + date_time(8)
    elapsed_time = end_time - beg_time
    write(0,*) "SCRIP: ", elapsed_time, " MSEC"
#endif
 
#ifdef W3_T38
    if(l_master)write(*,*)'new minval(grid1_center_lon) = ',minval(grid1_center_lon)
    if(l_master)write(*,*)'new maxval(grid1_center_lon) = ',maxval(grid1_center_lon)
#endif
 
    !.....notes: at this point we have the following useful variables:
    !       num_wts, e.g. num_wts=3....for first order conservative remapping,
    !       only the first one is relevant, the other two are for second-order
    !       remapping.
    !       max_links_map1, e.g. max_links_map1=1369,
    !       grid2_size, e.g. grid2_size=1849,
    !       wts_map1(num_wts,max_links_map1), e.g. wts_map1(3,1369),
    !       grid1_add_map1(max_links_map1),  e.g. grid1_add_map1(1369),
    !       grid2_add_map1(max_links_map1), e.g. grid2_add_map1(1369),
    !       grid2_frac(grid2_size), e.g. grid2_frac(1849),
    !       (see earlier versions for notes re: equivalency in netcdf/matlab)
    !
    !.....test output (optional)
    !
    !.....note re: notation: I use k for the combined i/j array, similar to isea,
    !            but not necessarily the same as isea since some points may
    !            be land etc.
#ifdef W3_T38
    if(l_master)then
      do k=1,grid2_size
        write(403,*)grid2_frac(k)
      end do
      do ilink=1,max_links_map1
        write(405,'(999(1x,f20.7))')(wts_map1(iw,ilink),iw=1,num_wts)
      end do
      do ilink=1,max_links_map1
        write(406,'(i20)')grid1_add_map1(ilink) ! equivalent to
        ! my "src_address"
        write(407,'(i20)')grid2_add_map1(ilink) ! equivalent to
        ! my "dst_address"
      end do
    endif
#endif
 
    !.....organize results and return to wmghgh.
 
    ! what we need, for purpose of feeding back to ww3, for each dst grid node:
    !        a) the set of src grid nodes, in terms of isea, for which
    !           weights are available
    !        b) the corresponding set of weights
 
    ALLOCATE ( wgtdata(grid2_size), stat=istat )
    check_alloc_status( istat )
 
    !.....step 1: count up NR0, NR1, NR2, NRL, NLOC (NR1 and NLOC are denoted
    !             "n" here)
    !     It is especially important to determine how large npnts gets,
    !              so that we can allocate arrays properly
    wgtdata%NR0=0
    wgtdata%NR2=0
    wgtdata%NRL=0
    wgtdata%N=0
 
    ni_src=grid1_dims(1)
    ngoodpnts=0
    nbadpnts=0
 
    DO ilink=1,max_links_map1
 
      !........note that this pair of if-thens *must* be consistent with the
      !........single if-then below
      IF ((grid2_frac(grid2_add_map1(ilink))>frac_lowest) .AND.     &
           (grid2_frac(grid2_add_map1(ilink))<frac_highest)) THEN
        !     then consider this link either to include, or to print a warning
        IF( (wts_map1(1,ilink)>=wt_lowest) .AND.                   &
             (wts_map1(1,ilink)<=wt_highest) )    THEN
          ksrc=grid1_add_map1(ilink)
          jsrc=int((ksrc-1)/ni_src)+1
          isrc=ksrc-(jsrc-1)*ni_src
 
          IF (mapsta_src(jsrc,isrc).EQ.0) THEN ! excluded point
            wgtdata(grid2_add_map1(ilink))%NR0    =              &
                 wgtdata(grid2_add_map1(ilink))%NR0 + 1
            IF (mapst2_src(jsrc,isrc).EQ.0)THEN
              wgtdata(grid2_add_map1(ilink))%NRL =              &
                   wgtdata(grid2_add_map1(ilink))%NRL + 1
            ENDIF
          ELSE IF (abs(mapsta_src(jsrc,isrc)).EQ.1) THEN
            ! sea point
            wgtdata(grid2_add_map1(ilink))%N=                    &
                 wgtdata(grid2_add_map1(ilink))%N+1
          ELSE IF (abs(mapsta_src(jsrc,isrc)).EQ.2) THEN
            ! bnd point
            wgtdata(grid2_add_map1(ilink))%NR2    =              &
                 wgtdata(grid2_add_map1(ilink))%NR2 + 1
          END IF
          ngoodpnts=ngoodpnts+1
        ELSEIF ( (grid1_frac(grid1_add_map1(ilink))>frac_lowest)   &
             .AND. (grid1_frac(grid1_add_map1(ilink))<frac_highest)&
             ) THEN
          ! note that we don't bother giving a warning for the
          ! cases where grid1_frac is strange.
          nbadpnts=nbadpnts+1
          ! we also don't bother giving a warning if we've already
          ! given a lot of warnings (we keep counting, though)
          IF((nbadpnts.LT.5).AND.l_master)THEN ! print warnings
            WRITE(*,'(A)')'%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%'
            WRITE(*,'(A)')'WARNING: SCRIP weight problem '
            WRITE(*,'(4x,A,I7,A,I7)')'ilink = ',ilink,' out of ',&
                 max_links_map1
            WRITE(*,'(4x,A,I7)')'grid1_add_map1(ilink) = ',      &
                 grid1_add_map1(ilink)
            WRITE(*,'(4x,A,I7)')'grid2_add_map1(ilink) = ',      &
                 grid2_add_map1(ilink)
            WRITE(*,'(4x,A,E12.4)')'wts_map1(1,ilink) = ',       &
                 wts_map1(1,ilink)
            WRITE(*,'(4x,A,F10.5)')                              &
                 'grid1_frac(grid1_add_map1(ilink)) = ',              &
                 grid1_frac(grid1_add_map1(ilink))
            WRITE(*,'(4x,A,F10.5)')                              &
                 'grid2_frac(grid2_add_map1(ilink)) = ',              &
                 grid2_frac(grid2_add_map1(ilink))
            WRITE(*,'(4x,A,F10.5)')'grid1_center_lat = ',        &
                 grid1_center_lat(grid1_add_map1(ilink))
            WRITE(*,'(4x,A,F10.5)')'grid1_center_lon = ',        &
                 grid1_center_lon(grid1_add_map1(ilink))
            WRITE(*,'(4x,A,F10.5)')'grid2_center_lat = ',        &
                 grid2_center_lat(grid2_add_map1(ilink))
            WRITE(*,'(4x,A,F10.5)')'grid2_center_lon = ',        &
                 grid2_center_lon(grid2_add_map1(ilink))
            WRITE(*,'(A)')'%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%'
          ENDIF
        ENDIF
      ENDIF
    END DO
    IF((nbadpnts.GT.0).AND.l_master)THEN
      WRITE(*,'(4x,A,I5,A)')'We had problems in ',nbadpnts, &
           ' points.'
      WRITE(*,'(4x,I8,A)')ngoodpnts,' points appear to be OK.'
    ENDIF
    IF( (nbadpnts.GT.(ngoodpnts/30)) .AND.l_master )THEN
      WRITE(*,'(4x,A)')'Error: excessive SCRIP failure. Stopping.'
      stop 'wmscrpmd, case 1'
    ENDIF
 
    !.....step 2: allocate according to the size "n" determined above
    DO kdst=1,grid2_size
      ALLOCATE ( wgtdata(kdst)%W(wgtdata(kdst)%N), stat=istat )
      check_alloc_status( istat )
      ALLOCATE ( wgtdata(kdst)%K(wgtdata(kdst)%N), stat=istat )
      check_alloc_status( istat )
      wgtdata(kdst)%N=0
    END DO
 
    !.....step 3: save weights
    DO ilink=1,max_links_map1
 
      ksrc=grid1_add_map1(ilink)
      jsrc=int((ksrc-1)/ni_src)+1
      isrc=ksrc-(jsrc-1)*ni_src
 
      !........note that this single if-then *must* be consistent with the
      !........pair of if-thens above
      IF ((grid2_frac(grid2_add_map1(ilink))>frac_lowest) .AND.      &
           &       (grid2_frac(grid2_add_map1(ilink))<frac_highest) .AND.     &
           &       (wts_map1(1,ilink)>=wt_lowest) .AND.                       &
           &       (wts_map1(1,ilink)<=wt_highest))THEN
        IF (abs(mapsta_src(jsrc,isrc)).EQ.1) THEN ! sea point
          wgtdata(grid2_add_map1(ilink))%N=                        &
               &             wgtdata(grid2_add_map1(ilink))%N+1
          wgtdata(grid2_add_map1(ilink))%W(wgtdata(                &
               &             grid2_add_map1(ilink))%N)=wts_map1(1,ilink)
          wgtdata(grid2_add_map1(ilink))%K(wgtdata(                &
               &             grid2_add_map1(ilink))%N)=grid1_add_map1(ilink)
        ENDIF
      ENDIF
    END DO
 
    !.....step 4: re-normalize weights. This is necessary because we called
    !.....scrip without the mask. Now that we have the mask in place, we need
    !.....to re-normalize the weights.
    DO kdst=1,grid2_size
      IF (wgtdata(kdst)%N > 0) THEN
        wtsum=zero
        DO ipnt=1,wgtdata(kdst)%N
          wtsum=wtsum+wgtdata(kdst)%W(ipnt)
        END DO
        DO ipnt=1,wgtdata(kdst)%N
          wgtdata(kdst)%W(ipnt)=wgtdata(kdst)%W(ipnt)/wtsum
        END DO
      END IF
    END DO
 
    CALL scrip_clear
  END SUBROUTINE scrip_wrapper
 
  !/ ------------------------------------------------------------------- /
  SUBROUTINE get_scrip_info_unstructured (ID_GRD,                  &
       &   GRID_CENTER_LON, GRID_CENTER_LAT,                             &
       &   GRID_CORNER_LON, GRID_CORNER_LAT, GRID_MASK,                  &
       &   GRID_DIMS, GRID_SIZE, GRID_CORNERS, GRID_RANK)
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III                     |
    !/                  | M. Dutour, A. Roland              |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         10-Dec-2014 !
    !/                  +-----------------------------------+
    !/
    !  1. Original author :
    !
    !     Mathieu Dutour Sikiric, IRB & Aron Roland, Z&P
    !
    !  2. Last update :
    !
    !     See revisions.
    !
    !  3. Revisions :
    !
    !     20-Feb-2012 : Origination
    !/    10-Dec-2014 : Add checks for allocate status      ( version 5.04 )
    !
    !  4. Copyright :
    !
    !  5. Purpose :
    !
    !     Compute grid arrays for scrip for a specific unstructured grid
    !     For interior vertices, we select for every triangle the barycenter
    !     of the triangle. So to every node contained in N triangles we associate
    !     a domain with N corners. Every one of those corners is contained
    !     in 3 different domains.
    !     For nodes that are on the boundary, we have to proceed differently.
    !     For every such node, we have NEIGHBOR_PREV and NEIGHBOR_NEXT which
    !     are the neighbor on each side of the boundary.
    !     We put a corner on the middle of the edge. We also put a corner
    !     on the node itself.
    !     Note that instead of taking barycenters, we could have taken
    !     Voronoi vertices, but this carries danger since Voronoi vertices
    !     can be outside of the triangle. And it leaves open how to treat
    !     the boundary.
    !
    !  6. Method :
    !
    !  7. Parameters, Variables and types :
    !
    !  8. Called by :
    !
    !     Subroutine get_scrip_info
    !
    !  9. Subroutines and functions used :
    !
    ! 10. Error messages:
    !
    ! 11. Remarks :
    !
    ! 12. Structure :
    !
    ! 13. Switches :
    !
    ! 14. Source code :
    USE w3servmd, ONLY: extcde
    USE w3gdatmd, ONLY : grids
    IMPLICIT NONE
    INTEGER, INTENT(IN)      :: ID_GRD
    real*8, INTENT(OUT), ALLOCATABLE :: grid_center_lon(:)
    real*8, INTENT(OUT), ALLOCATABLE :: grid_center_lat(:)
    LOGICAL, INTENT(OUT), ALLOCATABLE :: GRID_MASK(:)
    real*8, INTENT(OUT), ALLOCATABLE :: grid_corner_lon(:,:)
    real*8, INTENT(OUT), ALLOCATABLE :: grid_corner_lat(:,:)
    INTEGER, INTENT(OUT), ALLOCATABLE :: GRID_DIMS(:)
    INTEGER, INTENT(OUT) :: GRID_SIZE, GRID_CORNERS, GRID_RANK
 
    INTEGER DIRAPPROACH, DUALAPPROACH, THEAPPROACH
    INTEGER MNE, MNP, IE, IP, I
    INTEGER NBPLUS, NBMINUS
    INTEGER I1, I2, I3
    real*8 :: elon1, elon2, elon3, elon, elonc
    real*8 :: elat1, elat2, elat3, elat, elatc
    REAL *8 :: DELTALON12, DELTALON13, DELTALAT12, DELTALAT13
    REAL *8 :: THEDET
    real*8 :: pt(3,2)
    INTEGER, POINTER :: IOBP(:), TRIGINCD(:)
    INTEGER, POINTER :: NEIGHBOR_PREV(:), NEIGHBOR_NEXT(:)
    INTEGER, POINTER :: NBASSIGNEDCORNER(:), LISTNBCORNER(:)
    INTEGER, POINTER :: STATUS(:), NEXTVERT(:), PREVVERT(:), FINALVERT(:)
    INTEGER :: MAXCORNER, NBCORNER
    INTEGER :: IDX, IPNEXT, IPPREV, NB, INEXT, IPREV
    real*8, POINTER :: lon_cent_trig(:), lat_cent_trig(:)
    real*8 :: elonip, elonnext, elonprev, elonn, elonp
    real*8 :: elatip, elatnext, elatprev, elatn, elatp
    INTEGER :: ISFINISHED, ZPREV
    INTEGER :: DODEBUG
    grid_rank=2
    dirapproach=1
    dualapproach=2
    theapproach=dualapproach
    mne=grids(id_grd)%NTRI
    mnp=grids(id_grd)%NX
    IF (theapproach .EQ. dirapproach) THEN
      ALLOCATE(grid_center_lon(mne), stat=istat)
      check_alloc_status( istat )
      ALLOCATE(grid_center_lat(mne), stat=istat)
      check_alloc_status( istat )
      ALLOCATE(grid_corner_lon(3,mne), stat=istat)
      check_alloc_status( istat )
      ALLOCATE(grid_corner_lat(3,mne), stat=istat)
      check_alloc_status( istat )
      ALLOCATE(grid_mask(mne), stat=istat)
      check_alloc_status( istat )
      DO ie=1,mne
        i1=grids(id_grd)%TRIGP(1,ie)
        i2=grids(id_grd)%TRIGP(2,ie)
        i3=grids(id_grd)%TRIGP(3,ie)
        elon1=grids(id_grd)%XGRD(1,i1)
        elon2=grids(id_grd)%XGRD(1,i2)
        elon3=grids(id_grd)%XGRD(1,i3)
        elat1=grids(id_grd)%YGRD(1,i1)
        elat2=grids(id_grd)%YGRD(1,i2)
        elat3=grids(id_grd)%YGRD(1,i3)
        elon=(elon1 + elon2 + elon3)/3
        elat=(elat1 + elat2 + elat3)/3
        grid_center_lon(ie)=elon
        grid_center_lat(ie)=elat
        grid_corner_lon(1,ie)=elon1
        grid_corner_lon(2,ie)=elon2
        grid_corner_lon(3,ie)=elon3
        grid_corner_lat(1,ie)=elat1
        grid_corner_lat(2,ie)=elat2
        grid_corner_lat(3,ie)=elat3
        grid_mask(ie)=.true.
      END DO
      grid_corners=3
    END IF
    IF (theapproach .EQ. dualapproach) THEN
      ALLOCATE(trigincd(mnp), stat=istat)
      check_alloc_status( istat )
      ALLOCATE(iobp(mnp), stat=istat)
      check_alloc_status( istat )
      ALLOCATE(neighbor_next(mnp), stat=istat)
      check_alloc_status( istat )
      ALLOCATE(neighbor_prev(mnp), stat=istat)
      check_alloc_status( istat )
      ALLOCATE(nbassignedcorner(mnp), stat=istat)
      check_alloc_status( istat )
      ALLOCATE(listnbcorner(mnp), stat=istat)
      check_alloc_status( istat )
 
      ALLOCATE(status(mnp), stat=istat)
      check_alloc_status( istat )
      ALLOCATE(nextvert(mnp), stat=istat)
      check_alloc_status( istat )
      ALLOCATE(prevvert(mnp), stat=istat)
      check_alloc_status( istat )
      ALLOCATE(finalvert(mnp), stat=istat)
      check_alloc_status( istat )
      ALLOCATE(lon_cent_trig(mne), stat=istat)
      check_alloc_status( istat )
      ALLOCATE(lat_cent_trig(mne), stat=istat)
      check_alloc_status( istat )
 
      CALL get_unstructured_vertex_degree (mnp, mne,             &
           grids(id_grd)%TRIGP, trigincd)
      CALL get_boundary(mnp, mne, grids(id_grd)%TRIGP, iobp,     &
           neighbor_prev, neighbor_next)
 
      ! Find max number of corners
      maxcorner=0
      DO ip=1,mnp
        IF (neighbor_next(ip) .EQ. 0) THEN
          nbcorner=trigincd(ip)
        ELSE
          nbcorner=trigincd(ip) + 3
        END IF
        listnbcorner(ip)=nbcorner
        IF (nbcorner .GT. maxcorner) THEN
          maxcorner=nbcorner
        END IF
      END DO
      grid_corners=maxcorner
 
      ALLOCATE(grid_center_lon(mnp), stat=istat)
      check_alloc_status( istat )
      ALLOCATE(grid_center_lat(mnp), stat=istat)
      check_alloc_status( istat )
      ALLOCATE(grid_corner_lon(maxcorner,mnp), stat=istat)
      check_alloc_status( istat )
      ALLOCATE(grid_corner_lat(maxcorner,mnp), stat=istat)
      check_alloc_status( istat )
      ALLOCATE(grid_mask(mnp), stat=istat)
      check_alloc_status( istat )
 
      ! Add first three corners for boundaries
      nbassignedcorner(:)=0
      DO ip=1,mnp
        grid_mask(ip)=.true.
        IF (neighbor_next(ip) .GT. 0) THEN
          ipnext=neighbor_next(ip)
          ipprev=neighbor_prev(ip)
          elonip=dble(grids(id_grd)%XGRD(1,ip))
          elatip=dble(grids(id_grd)%YGRD(1,ip))
          elonnext=dble(grids(id_grd)%XGRD(1,ipnext))
          elatnext=dble(grids(id_grd)%YGRD(1,ipnext))
          elonprev=dble(grids(id_grd)%XGRD(1,ipprev))
          elatprev=dble(grids(id_grd)%YGRD(1,ipprev))
 
          ! Periodicity fix for corner node
          IF ( abs(elonip - elonnext) .GT. 180.0 ) THEN
            elonnext = elonnext -sign(360.0d0,(elonip - elonnext))
          ENDIF
          IF ( abs(elonip - elonprev) .GT. 180.0 ) THEN
            elonprev = elonprev -sign(360.0d0,(elonip - elonprev))
          ENDIF
 
          elonn=(elonip+elonnext)/2.0
          elatn=(elatip+elatnext)/2.0
          elonp=(elonip+elonprev)/2.0
          elatp=(elatip+elatprev)/2.0
 
 
          grid_corner_lon(1,ip)=elonn
          grid_corner_lat(1,ip)=elatn
          grid_corner_lon(2,ip)=elonip
          grid_corner_lat(2,ip)=elatip
          grid_corner_lon(3,ip)=elonp
          grid_corner_lat(3,ip)=elatp
          nbassignedcorner(ip)=3
        END IF
      END DO
 
      ! Compute centers
      DO ip=1,mnp
        grid_center_lon(ip)=dble(grids(id_grd)%XGRD(1,ip))
        grid_center_lat(ip)=dble(grids(id_grd)%YGRD(1,ip))
      END DO
 
      ! Check triangle node orientation
      ! Compute triangle centers
      nbplus=0
      nbminus=0
      DO ie=1,mne
        i1=grids(id_grd)%TRIGP(1,ie)
        i2=grids(id_grd)%TRIGP(2,ie)
        i3=grids(id_grd)%TRIGP(3,ie)
        pt(1,1)=dble(grids(id_grd)%XGRD(1,i1))
        pt(2,1)=dble(grids(id_grd)%XGRD(1,i2))
        pt(3,1)=dble(grids(id_grd)%XGRD(1,i3))
        pt(1,2)=dble(grids(id_grd)%YGRD(1,i1))
        pt(2,2)=dble(grids(id_grd)%YGRD(1,i2))
        pt(3,2)=dble(grids(id_grd)%YGRD(1,i3))
 
        CALL fix_periodcity(pt)
 
        elon1 = pt(1,1)
        elon2 = pt(2,1)
        elon3 = pt(3,1)
        elat1 = pt(1,2)
        elat2 = pt(2,2)
        elat3 = pt(3,2)
 
        deltalon12=elon2 - elon1
        deltalon13=elon3 - elon1
        deltalat12=elat2 - elat1
        deltalat13=elat3 - elat1
        thedet=deltalon12*deltalat13 - deltalon13*deltalat12
        IF (thedet.GT.0) THEN
          nbplus=nbplus+1
        END IF
        IF (thedet.LT.0) THEN
          nbminus=nbminus+1
        END IF
        elon=(elon1 + elon2 + elon3)/3.0
        elat=(elat1 + elat2 + elat3)/3.0
 
 
        lon_cent_trig(ie)=elon
        lat_cent_trig(ie)=elat
 
      END DO
      dodebug=0
      IF (dodebug.EQ.1) THEN
        print *, 'nbplus=', nbplus, ' nbminus=', nbminus
      END IF
 
      status(:) = 0
      nextvert(:) = 0
      prevvert(:) = 0
      DO ie=1,mne
        DO i=1,3
          CALL triang_indexes(i, inext, iprev)
          ip=grids(id_grd)%TRIGP(i,ie)
          ipnext=grids(id_grd)%TRIGP(inext,ie)
          ipprev=grids(id_grd)%TRIGP(iprev,ie)
          IF (status(ip).EQ.0) THEN
            IF (neighbor_next(ip).EQ.0) THEN
              status(ip)=1
              finalvert(ip)=ipprev
              prevvert(ip)=ipprev
              nextvert(ip)=ipnext
            ELSE
              IF (neighbor_prev(ip).EQ.ipprev) THEN
                status(ip)=1
                prevvert(ip)=ipprev
                nextvert(ip)=ipnext
                finalvert(ip)=neighbor_next(ip)
              END IF
            END IF
          END IF
        END DO
      END DO
      status(:)=0
      DO
        isfinished=1
        DO ie=1,mne
          elon=lon_cent_trig(ie)
          elat=lat_cent_trig(ie)
          DO i=1,3
            CALL triang_indexes(i, inext, iprev)
            ip=grids(id_grd)%TRIGP(i,ie)
            ipnext=grids(id_grd)%TRIGP(inext,ie)
            ipprev=grids(id_grd)%TRIGP(iprev,ie)
            IF (status(ip).EQ.0) THEN
              isfinished=0
              zprev=prevvert(ip)
              IF (zprev.EQ.ipprev) THEN
                idx=nbassignedcorner(ip)
                idx=idx+1
                grid_corner_lon(idx,ip)=elon
                grid_corner_lat(idx,ip)=elat
                nbassignedcorner(ip)=idx
                prevvert(ip)=ipnext
                IF (ipnext.EQ.finalvert(ip)) THEN
                  status(ip)=1
                END IF
              END IF
            END IF
          END DO
        END DO
        IF (isfinished.EQ.1) THEN
          EXIT
        END IF
      END DO
      DO ip=1,mnp
        IF (nbassignedcorner(ip).NE.listnbcorner(ip)) THEN
          WRITE(*,*) 'Incoherent number at IP=', ip
          WRITE(*,*) '  NbAssignedCorner(IP)=', nbassignedcorner(ip)
          WRITE(*,*) '  ListNbCorner(IP)=', listnbcorner(ip)
          WRITE(*,*) '  N_N=', neighbor_next(ip), 'N_P=', neighbor_prev(ip)
          WRITE(*,*) '  TrigIncd=', trigincd(ip)
          stop 'wmscrpmd, case 2'
        END IF
      END DO
 
      ! if the number of corner is below threshold, we have to
      ! add some more.
      DO ip=1,mnp
        nb=nbassignedcorner(ip)
        IF (nb .LT. maxcorner) THEN
          elon=grid_corner_lon(nb,ip)
          elat=grid_corner_lat(nb,ip)
          DO idx=nb+1,maxcorner
            grid_corner_lon(idx,ip)=elon
            grid_corner_lat(idx,ip)=elat
          END DO
        END IF
      END DO
      DEALLOCATE(nbassignedcorner, stat=istat)
      check_dealloc_status( istat )
      DEALLOCATE(listnbcorner, stat=istat)
      check_dealloc_status( istat )
      DEALLOCATE(trigincd, stat=istat)
      check_dealloc_status( istat )
      DEALLOCATE(iobp, stat=istat)
      check_dealloc_status( istat )
      DEALLOCATE(neighbor_prev, stat=istat)
      check_dealloc_status( istat )
      DEALLOCATE(neighbor_next, stat=istat)
      check_dealloc_status( istat )
      DEALLOCATE(status, stat=istat)
      check_dealloc_status( istat )
      DEALLOCATE(nextvert, stat=istat)
      check_dealloc_status( istat )
      DEALLOCATE(prevvert, stat=istat)
      check_dealloc_status( istat )
      DEALLOCATE(finalvert, stat=istat)
      check_dealloc_status( istat )
      DEALLOCATE(lon_cent_trig, stat=istat)
      check_dealloc_status( istat )
      DEALLOCATE(lat_cent_trig, stat=istat)
      check_dealloc_status( istat )
 
      ALLOCATE(grid_dims(2), stat=istat)
      check_alloc_status( istat )
      grid_dims(1)=mnp
      grid_dims(2)=1
      grid_size=mnp
    END IF
  END SUBROUTINE get_scrip_info_unstructured
 
  !/ ------------------------------------------------------------------- /
  SUBROUTINE get_scrip_info_structured (ID_GRD,                     &
       &   GRID_CENTER_LON, GRID_CENTER_LAT,                              &
       &   GRID_CORNER_LON, GRID_CORNER_LAT, GRID_MASK,                   &
       &   GRID_DIMS, GRID_SIZE, GRID_CORNERS, GRID_RANK)
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III                     |
    !/                  | M. Dutour, A. Roland              |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         10-Dec-2014 !
    !/                  +-----------------------------------+
    !/
    !  1. Original author :
    !
    !     Mathieu Dutour Sikiric, IRB & Aron Roland, Z&P
    !
    !  2. Last update :
    !
    !     See revisions.
    !
    !  3. Revisions :
    !
    !     20-Feb-2012 : Origination, this is adapted from Erick Rogers
    !                   code by splitting the code into sections.
    !/    10-Dec-2014 : Add checks for allocate status      ( version 5.04 )
    !
    !  4. Copyright :
    !
    !  5. Purpose :
    !
    !     Compute grid arrays needed for scrip for a specific structured grid.
    !     This is adapted from Erick Rogers original code by splitting
    !     the original scrip_wrapper function
    !
    !  6. Method :
    !
    !  7. Parameters, Variables and types :
    !
    !  8. Called by :
    !
    !     Subroutine get_scrip_info
    !
    !  9. Subroutines and functions used :
    !
    ! 10. Error messages:
    !
    ! 11. Remarks :
    !
    ! 12. Structure :
    !
    ! 13. Switches :
    !
    ! 14. Source code :
    USE w3servmd, ONLY: extcde
    USE w3gdatmd, ONLY : grids
    USE scrip_constants, ONLY : half
    IMPLICIT NONE
    INTEGER, INTENT(IN)      :: ID_GRD
    real*8, INTENT(OUT), ALLOCATABLE :: grid_center_lon(:)
    real*8, INTENT(OUT), ALLOCATABLE :: grid_center_lat(:)
    LOGICAL, INTENT(OUT), ALLOCATABLE :: GRID_MASK(:)
    real*8, INTENT(OUT), ALLOCATABLE :: grid_corner_lon(:,:)
    real*8, INTENT(OUT), ALLOCATABLE :: grid_corner_lat(:,:)
    INTEGER, INTENT(OUT), ALLOCATABLE :: GRID_DIMS(:)
    INTEGER, INTENT(OUT) :: GRID_SIZE, GRID_CORNERS, GRID_RANK
    !
    real*8, ALLOCATABLE    :: xin_grd(:,:), yin_grd(:,:)
    real*8, ALLOCATABLE    :: dxdp_grd(:,:), dxdq_grd(:,:)
    real*8, ALLOCATABLE    :: dydp_grd(:,:), dydq_grd(:,:)
    INTEGER :: N1, N2, NI, NJ
    INTEGER :: IREC, J, I
    grid_rank=2
    n1=SIZE(grids(id_grd)%XGRD,1)
    n2=SIZE(grids(id_grd)%XGRD,2)
    ALLOCATE(xin_grd(n1,n2), stat=istat)
    check_alloc_status( istat )
    ALLOCATE(yin_grd(n1,n2), stat=istat)
    check_alloc_status( istat )
    ALLOCATE(dxdp_grd(n1,n2), stat=istat)
    check_alloc_status( istat )
    ALLOCATE(dxdq_grd(n1,n2), stat=istat)
    check_alloc_status( istat )
    ALLOCATE(dydp_grd(n1,n2), stat=istat)
    check_alloc_status( istat )
    ALLOCATE(dydq_grd(n1,n2), stat=istat)
    check_alloc_status( istat )
 
    xin_grd=dble(grids(id_grd)%XGRD)
    yin_grd=dble(grids(id_grd)%YGRD)
    dxdp_grd=dble(grids(id_grd)%DXDP)
    dxdq_grd=dble(grids(id_grd)%DXDQ)
    dydp_grd=dble(grids(id_grd)%DYDP)
    dydq_grd=dble(grids(id_grd)%DYDQ)
 
    ALLOCATE(grid_dims(grid_rank), stat=istat)
    check_alloc_status( istat )
    grid_dims(1)=n2
    ni=n2
    grid_dims(2)=n1
    nj=n1
 
    grid_size=ni*nj ! hardwired: logically rectangular grid
    grid_corners=4  ! hardwired: each cell has 4 corners
 
    !.....notes: unfortunately, scrip only works for spherical coordinates.
    !         thus, if we want to have a multi-grid case in meters, we have to
    !         fake it. fortunately, it should be pretty rare to have a multi-grid
    !         case in meters.
 
    ALLOCATE(grid_center_lon(ni*nj), stat=istat)
    check_alloc_status( istat )
    ALLOCATE(grid_center_lat(ni*nj), stat=istat)
    check_alloc_status( istat )
    ALLOCATE(grid_corner_lon(4,ni*nj), stat=istat)
    check_alloc_status( istat )
    ALLOCATE(grid_corner_lat(4,ni*nj), stat=istat)
    check_alloc_status( istat )
    ALLOCATE(grid_mask(ni*nj), stat=istat)
    check_alloc_status( istat )
 
    !.....notes: this "gridshift" variable is included because SCRIP sometimes
    !            has trouble when grids cell locations are identical between
    !            the two grids. Thus we apply this to one of the two grids.
 
    !.....notes: The following block of code could be converted to a subroutine.
    !            Since it is called twice, this would save a little space.
    irec=0
    DO j=1,nj
      DO i=1,ni
        irec=irec+1
        grid_center_lon(irec)=xin_grd(j,i)
        grid_center_lat(irec)=yin_grd(j,i)
        grid_mask(irec)=.true.
 
        !..notes: normally, we'd apply the mask like this:
        !           if(abs(mapsta_src(j,i)).eq.1)then
        !              grid1_mask(irec)=.true.
        !           else
        !              grid1_mask(irec)=.false.
        !           endif
        !..but unfortunately, WMGHGH needs information about the overlaying high-res
        !           cells, even those that are masked, for calculating
        !           NRL, NR0, NR1, NR2.
 
        !...........corner 1 : halfway to i-1,j-1
        grid_corner_lon(1,irec)=grid_center_lon(irec)-              &
             &        half*dxdp_grd(j,i)-half*dxdq_grd(j,i)
        grid_corner_lat(1,irec)=grid_center_lat(irec)-              &
             &        half*dydp_grd(j,i)-half*dydq_grd(j,i)
 
        !...........corner 2: halfway to i+1,j-1
        grid_corner_lon(2,irec)=grid_center_lon(irec)+              &
             &        half*dxdp_grd(j,i)-half*dxdq_grd(j,i)
        grid_corner_lat(2,irec)=grid_center_lat(irec)+              &
             &        half*dydp_grd(j,i)-half*dydq_grd(j,i)
 
        !...........corner 3: halfway to i+1,j+1
        grid_corner_lon(3,irec)=grid_center_lon(irec)+              &
             &        half*dxdp_grd(j,i)+half*dxdq_grd(j,i)
        grid_corner_lat(3,irec)=grid_center_lat(irec)+              &
             &        half*dydp_grd(j,i)+half*dydq_grd(j,i)
 
        !...........corner 4: halfway to i-1,j+1
        grid_corner_lon(4,irec)=grid_center_lon(irec)-              &
             &        half*dxdp_grd(j,i)+half*dxdq_grd(j,i)
        grid_corner_lat(4,irec)=grid_center_lat(irec)-              &
             &        half*dydp_grd(j,i)+half*dydq_grd(j,i)
      END DO
    END DO
  END SUBROUTINE get_scrip_info_structured
 
  !/ ------------------------------------------------------------------- /
  SUBROUTINE get_scrip_info(ID_GRD,                                 &
       &   GRID_CENTER_LON, GRID_CENTER_LAT,                              &
       &   GRID_CORNER_LON, GRID_CORNER_LAT, GRID_MASK,                   &
       &   GRID_DIMS, GRID_SIZE, GRID_CORNERS, GRID_RANK)
    !  1. Original author :
    !
    !     Mathieu Dutour Sikiric, IRB & Aron Roland, Z&P
    !
    !  2. Last update :
    !
    !     See revisions.
    !
    !  3. Revisions :
    !
    !     20-Feb-2012 : Origination, this is adapted from Erick Rogers
    !                   code by splitting the code into sections.
    !
    !  4. Copyright :
    !
    !  5. Purpose :
    !
    !     Compute grid for scrip for a specific structured grid.
    !     This is adapted from Erick Rogers code by making it cleaner.
    !
    !  6. Method :
    !
    !  7. Parameters, Variables and types :
    !
    !  8. Called by :
    !
    !     Subroutine scrip_wrapper
    !
    !  9. Subroutines and functions used :
    !
    ! 10. Error messages:
    !
    ! 11. Remarks :
    !
    ! 12. Structure :
    !
    ! 13. Switches :
    !
    ! 14. Source code :
    USE w3servmd, ONLY: extcde
    USE w3gdatmd, ONLY : grids, ungtype
    IMPLICIT NONE
    INTEGER, INTENT(IN)      :: ID_GRD
    real*8, INTENT(OUT), ALLOCATABLE :: grid_center_lon(:)
    real*8, INTENT(OUT), ALLOCATABLE :: grid_center_lat(:)
    LOGICAL, INTENT(OUT), ALLOCATABLE :: GRID_MASK(:)
    real*8, INTENT(OUT), ALLOCATABLE :: grid_corner_lon(:,:)
    real*8, INTENT(OUT), ALLOCATABLE :: grid_corner_lat(:,:)
    INTEGER, INTENT(OUT), ALLOCATABLE :: GRID_DIMS(:)
    INTEGER, INTENT(OUT) :: GRID_SIZE, GRID_CORNERS, GRID_RANK
    real*8 :: dlon1, dlat1, dlon2, dlat2, thedet
    INTEGER :: I, J
    INTEGER :: IC, JC, IP, CHECKSIGNS, NBPLUS, NBMINUS, NBZERO
    INTEGER :: PRINTDATA, PRINTMINMAX
    real*8 :: minlon, minlat, maxlon, maxlat
    real*8 :: minloncorner, maxloncorner, minlatcorner, maxlatcorner
    real*8 :: pt(3,2)
    IF (grids(id_grd)%GTYPE .EQ. ungtype) THEN
      CALL get_scrip_info_unstructured (id_grd,                       &
           &   grid_center_lon, grid_center_lat,                              &
           &   grid_corner_lon, grid_corner_lat, grid_mask,                   &
           &   grid_dims, grid_size, grid_corners, grid_rank)
    ELSE
      CALL get_scrip_info_structured (id_grd,                         &
           &   grid_center_lon, grid_center_lat,                              &
           &   grid_corner_lon, grid_corner_lat, grid_mask,                   &
           &   grid_dims, grid_size, grid_corners, grid_rank)
    END IF
    checksigns=1
    IF (checksigns.EQ.1) THEN
      nbplus=0
      nbminus=0
      nbzero=0
      DO ip=1,grid_size
        DO ic=1,grid_corners
          IF (ic.EQ.grid_corners) THEN
            jc=1
          ELSE
            jc=ic+1
          END IF
 
          pt(1,1) = grid_center_lon(ip)
          pt(1,2) = grid_center_lat(ip)
          pt(2,1) = grid_corner_lon(ic,ip)
          pt(2,2) = grid_corner_lat(ic,ip)
          pt(3,1) = grid_corner_lon(jc,ip)
          pt(3,2) = grid_corner_lat(jc,ip)
 
          CALL fix_periodcity(pt)
 
          dlon1=pt(2,1)-pt(1,1)
          dlon2=pt(3,1)-pt(1,1)
          dlat1=pt(2,2)-pt(1,2)
          dlat2=pt(3,2)-pt(1,2)
 
          thedet=dlon1*dlat2 - dlon2*dlat1
          IF (thedet.GT.1d-8) THEN
            nbplus=nbplus+1
          ELSE IF (thedet.LT.-1d-8) THEN
            nbminus=nbminus+1
          ELSE
            nbzero=nbzero+1
          END IF
        END DO
      END DO
 
      WRITE(*,*) 'SI nbPlus=', nbplus, ' nbMinus=', nbminus, ' nbZero=', nbzero
 
    END IF
  END SUBROUTINE get_scrip_info
 
  !/ ------------------------------------------------------------------- /
  SUBROUTINE scrip_info_renormalization(                            &
       &   GRID_CENTER_LON, GRID_CENTER_LAT,                              &
       &   GRID_CORNER_LON, GRID_CORNER_LAT, GRID_MASK,                   &
       &   GRID_DIMS, GRID_SIZE, GRID_CORNERS, GRID_RANK,                 &
       &   CONV_DX, CONV_DY, OFFSET, GRIDSHIFT)
    !  1. Original author :
    !
    !     Mathieu Dutour Sikiric, IRB & Aron Roland, Z&P
    !     Adapted from Erick Rogers scrip_wrapper
    !
    !  2. Last update :
    !
    !     See revisions.
    !
    !  3. Revisions :
    !
    !     20-Feb-2012 : Origination
    !
    !  4. Copyright :
    !
    !  5. Purpose :
    !
    !     This is adapted from Erick Rogers scrip_wrapper
    !     Purpose is to rescale according to whether the grid is spherical
    !     or not and to adjust by some small shift to keep SCRIP happy
    !     in situations where nodes of different grids overlay
    !
    !  6. Method :
    !
    !    We apply various transformations to the longitude latitude
    !    following here the transformations that were done only in
    !    finite difference meshes.
    !
    !  7. Parameters, Variables and types :
    !
    !  8. Called by :
    !
    !     Subroutine WMGHGH
    !
    !  9. Subroutines and functions used :
    !
    ! 10. Error messages:
    !
    ! 11. Remarks :
    !
    ! 12. Structure :
    !
    ! 13. Switches :
    !
    ! 14. Source code :
    IMPLICIT NONE
    real*8, INTENT(INOUT) :: grid_center_lon(:)
    real*8, INTENT(INOUT) :: grid_center_lat(:)
    LOGICAL, INTENT(IN) :: GRID_MASK(:)
    real*8, INTENT(INOUT) :: grid_corner_lon(:,:)
    real*8, INTENT(INOUT) :: grid_corner_lat(:,:)
    INTEGER, INTENT(IN) :: GRID_DIMS(:)
    INTEGER, INTENT(IN) :: GRID_SIZE, GRID_CORNERS, GRID_RANK
    real*8 :: conv_dx, conv_dy, offset, gridshift
    real*8 deg2rad
    !
    INTEGER :: I, J, IP
    real*8 :: minlon, minlat, maxlon, maxlat, hlon, hlat
    real*8 :: minloncorner, maxloncorner, minlatcorner, maxlatcorner
 
    DO i=1,grid_size
      grid_center_lon(i)=(grid_center_lon(i)+offset)/conv_dx +        &
           &       gridshift
      grid_center_lat(i)=grid_center_lat(i)/conv_dy +                 &
           &       gridshift
      IF(grid_center_lon(i)>360.0) THEN
        grid_center_lon(i)=grid_center_lon(i)-360.0
      END IF
      IF(grid_center_lon(i)<000.0) THEN
        grid_center_lon(i)=grid_center_lon(i)+360.0
      END IF
      DO j=1,grid_corners
        grid_corner_lon(j, i)=(grid_corner_lon(j, i)+offset)/conv_dx+ &
             &       gridshift
        grid_corner_lat(j, i)=grid_corner_lat(j, i)/conv_dy +         &
             &       gridshift
        IF(grid_corner_lon(j,i)>360.0) THEN
          grid_corner_lon(j,i)=grid_corner_lon(j,i)-360.0
        END IF
        IF(grid_corner_lon(j,i)<000.0) THEN
          grid_corner_lon(j,i)=grid_corner_lon(j,i)+360.0
        END IF
      END DO
    END DO
 
  END SUBROUTINE scrip_info_renormalization
 
  !/ ------------------------------------------------------------------- /
  SUBROUTINE triang_indexes(I, INEXT, IPREV)
    !  1. Original author :
    !
    !     Mathieu Dutour Sikiric, IRB & Aron Roland, Z&P
    !
    INTEGER, INTENT(IN)  :: I
    INTEGER, INTENT(OUT) :: INEXT, IPREV
    IF (i.EQ.1) THEN
      inext=3
    ELSE
      inext=i-1
    END IF
    IF (i.EQ.3) THEN
      iprev=1
    ELSE
      iprev=i+1
    END IF
  END SUBROUTINE triang_indexes
 
  !/ ------------------------------------------------------------------- /
  SUBROUTINE get_unstructured_vertex_degree (MNP, MNE, TRIGP,         &
       &          TRIGINCD)
    !  Written:
    !
    !    20-Feb-2012
    !
    !  Author:
    !
    !    Mathieu Dutour Sikiric, IRB & Aron Roland, Z&P
    !
    !  Parameters:
    !   Input:
    !    MNP: number of nodes
    !    INE: list of nodes
    !    MNE: number of triangles
    !   Output:
    !    TrigIncd (number of triangles contained by vertices
    !
    !  Description:
    !    this function returns the list of incidences
    !
    IMPLICIT NONE
    INTEGER, INTENT(IN) :: MNP, MNE
    INTEGER, INTENT(IN) :: TRIGP(:,:)
    INTEGER, INTENT(OUT) :: TRIGINCD(:)
    INTEGER :: IP, IE, I
    trigincd=0
    DO ie=1,mne
      DO i=1,3
        ip=trigp(i,ie)
        trigincd(ip)=trigincd(ip) + 1
      END DO
    END DO
  END SUBROUTINE get_unstructured_vertex_degree
 
  !/ ------------------------------------------------------------------- /
  SUBROUTINE get_boundary(MNP, MNE, TRIGP, IOBP, NEIGHBOR_PREV,       &
       &   NEIGHBOR_NEXT)
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III                     |
    !/                  | M. Dutour, A. Roland              |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         10-Dec-2014 !
    !/                  +-----------------------------------+
    !/
    !  Written:
    !
    !    20-Feb-2012
    !/    10-Dec-2014 : Add checks for allocate status      ( version 5.04 )
    !
    !  Author:
    !
    !    Mathieu Dutour Sikiric, IRB & Aron Roland, Z&P
    !
    !  Parameters:
    !   Input:
    !    MNP: number of nodes
    !    TRIGP: list of nodes
    !    MNE: number of triangles
    !   Output:
    !    NEIGHBOR
    !
    !  Description:
    !    if a node belong to a boundary, the function
    !    returns the neighbor of this point on one side.
    !    if the point is interior then the value 0 is set.
    !
    USE w3servmd, ONLY: extcde
    IMPLICIT NONE
 
    INTEGER, INTENT(IN)             :: MNP, MNE, TRIGP(3,MNE)
    INTEGER, INTENT(INOUT)          :: IOBP(MNP)
    INTEGER, INTENT(INOUT)          :: NEIGHBOR_PREV(MNP)
    INTEGER, INTENT(INOUT)          :: NEIGHBOR_NEXT(MNP)
 
    INTEGER, POINTER :: STATUS(:)
    INTEGER, POINTER :: COLLECTED(:)
    INTEGER, POINTER :: NEXTVERT(:)
    INTEGER, POINTER :: PREVVERT(:)
 
    INTEGER :: IE, I, IP, IP2, IP3
    INTEGER :: ISFINISHED, INEXT, IPREV
    INTEGER :: IPNEXT, IPPREV, ZNEXT, ZPREV
 
    ALLOCATE(status(mnp), stat=istat)
    check_alloc_status( istat )
    ALLOCATE(collected(mnp), stat=istat)
    check_alloc_status( istat )
    ALLOCATE(prevvert(mnp), stat=istat)
    check_alloc_status( istat )
    ALLOCATE(nextvert(mnp), stat=istat)
    check_alloc_status( istat )
    iobp = 0
    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 'wmscrpmd, case 3'
        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. 0) THEN
        iobp(ip)=1
      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 fix_periodcity(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
    real*8, INTENT(INOUT) :: pt(3,2)
    !     ----------------------------------------------------------------
    !
    !     Local variables.
    !     ----------------------------------------------------------------
    INTEGER :: I
    INTEGER :: R1GT180, R2GT180, R3GT180
    !     ----------------------------------------------------------------
    !
    !  4. Subroutines used :
    !
 
    !  5. Called by :
    !
    !      Name         Type  Module   Description
    !     ----------------------------------------------------------------
    !      GET_SCRIP_INFO_UNSTRUCTURED Subr. WMSCRPMD  Element center calculation
    !      GET_SCRIP_INFO              Subr. WMSCRPMD  Check signs
    !     ----------------------------------------------------------------
    !
    !  6. Error messages :
    !
    !       None.
    !
    !  7. Remarks :
    !
    !  8. Structure :
    !
    !  9. Switches :
    !
    ! 10. Source code :
    !/ ------------------------------------------------------------------- /
 
    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 of module WMSCRPMD -------------------------------------------- /
  !/
END MODULE wmscrpmd