#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
      REAL*8 , ALLOCATABLE :: GRID_AREA(:)
!/
!/ 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
!/T38      CHARACTER (LEN=10)      :: CDATE_TIME(3)
!/T38      INTEGER                 :: DATE_TIME(8)
!/T38      INTEGER                 :: ELAPSED_TIME, BEG_TIME, &
!/T38                                 END_TIME

! test output for input variables

!/T38      if(l_master)write(*,*)'flagll = ',flagll
!/T38      if(l_master)write(*,*)'gridshift = ',gridshift

! 
! 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
!/T38      write(*,*)'l_master = ',l_master
!/T38      if(l_master)then
!/T38      write(*,*)'conv_dx=', conv_dx
!/T38      write(*,*)'conv_dy=', conv_dy
!/T38      write(*,*)'offset = ',offset
!/T38      write(*,*)'grid1_size=', grid1_size
!/T38      write(*,*)'grid2_size=', grid2_size
!/T38      write(*,*)'l_read = ',l_read
!/T38      write(*,*)'minval(grid1_center_lon) = ',minval(grid1_center_lon)
!/T38      write(*,*)'maxval(grid1_center_lon) = ',maxval(grid1_center_lon)
!/T38      write(*,*)'minval(grid1_center_lat) = ',minval(grid1_center_lat)
!/T38      write(*,*)'maxval(grid1_center_lat) = ',maxval(grid1_center_lat)
!/T38      write(*,*)'minval(grid2_center_lon) = ',minval(grid2_center_lon)
!/T38      write(*,*)'maxval(grid2_center_lon) = ',maxval(grid2_center_lon)
!/T38      write(*,*)'minval(grid2_center_lat) = ',minval(grid2_center_lat)
!/T38      write(*,*)'maxval(grid2_center_lat) = ',maxval(grid2_center_lat)
!/T38      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

!/T38      call date_and_time (CDATE_TIME(1), CDATE_TIME(2), CDATE_TIME(3), DATE_TIME)
!/T38      beg_time = ((date_time(5)*60 + date_time(6))*60 +date_time(7))*1000 + date_time(8)
      CALL SCRIP(ID_SRC, ID_DST, L_MASTER, L_READ, L_TEST)
!/T38      call date_and_time (CDATE_TIME(1), CDATE_TIME(2), CDATE_TIME(3), DATE_TIME)
!/T38      end_time = ((date_time(5)*60 + date_time(6))*60 +date_time(7))*1000 + date_time(8)
!/T38      elapsed_time = end_time - beg_time
!/T38      write(0,*) "SCRIP: ", elapsed_time, " MSEC"

!/T38      if(l_master)write(*,*)'new minval(grid1_center_lon) = ',minval(grid1_center_lon)
!/T38      if(l_master)write(*,*)'new maxval(grid1_center_lon) = ',maxval(grid1_center_lon)

!.....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.
!/T38     if(l_master)then
!/T38     do k=1,grid2_size
!/T38        write(403,*)grid2_frac(k)
!/T38     end do
!/T38     do ilink=1,max_links_map1
!/T38        write(405,'(999(1x,f20.7))')(wts_map1(iw,ilink),iw=1,num_wts)
!/T38     end do
!/T38     do ilink=1,max_links_map1
!/T38        write(406,'(i20)')grid1_add_map1(ilink) ! equivalent to
!/T38                                                ! my "src_address"
!/T38        write(407,'(i20)')grid2_add_map1(ilink) ! equivalent to 
!/T38                                                ! my "dst_address"
!/T38     end do
!/T38     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(IE,1)
          I2=GRIDS(ID_GRD)%TRIGP(IE,2)
          I3=GRIDS(ID_GRD)%TRIGP(IE,3)
          ELON1=GRIDS(ID_GRD)%XYB(I1,1)
          ELON2=GRIDS(ID_GRD)%XYB(I2,1)
          ELON3=GRIDS(ID_GRD)%XYB(I3,1)
          ELAT1=GRIDS(ID_GRD)%XYB(I1,2)
          ELAT2=GRIDS(ID_GRD)%XYB(I2,2)
          ELAT3=GRIDS(ID_GRD)%XYB(I3,2)
          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)%XYB(IP,1))
            ELATIP=DBLE(GRIDS(ID_GRD)%XYB(IP,2))
            ELONNEXT=DBLE(GRIDS(ID_GRD)%XYB(IPNEXT,1))
            ELATNEXT=DBLE(GRIDS(ID_GRD)%XYB(IPNEXT,2))
            ELONPREV=DBLE(GRIDS(ID_GRD)%XYB(IPPREV,1))
            ELATPREV=DBLE(GRIDS(ID_GRD)%XYB(IPPREV,2))
 
            ! 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)%XYB(IP,1))
          GRID_CENTER_LAT(IP)=DBLE(GRIDS(ID_GRD)%XYB(IP,2))
        END DO

        ! Check triangle node orientation
        ! Compute triangle centers
        NBPLUS=0
        NBMINUS=0
        DO IE=1,MNE
          I1=GRIDS(ID_GRD)%TRIGP(IE,1)
          I2=GRIDS(ID_GRD)%TRIGP(IE,2)
          I3=GRIDS(ID_GRD)%TRIGP(IE,3)

          PT(1,1)=DBLE(GRIDS(ID_GRD)%XYB(I1,1))
          PT(2,1)=DBLE(GRIDS(ID_GRD)%XYB(I2,1))
          PT(3,1)=DBLE(GRIDS(ID_GRD)%XYB(I3,1))
          PT(1,2)=DBLE(GRIDS(ID_GRD)%XYB(I1,2))
          PT(2,2)=DBLE(GRIDS(ID_GRD)%XYB(I2,2))
          PT(3,2)=DBLE(GRIDS(ID_GRD)%XYB(I3,2))

          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(IE,I)
            IPNEXT=GRIDS(ID_GRD)%TRIGP(IE,INEXT)
            IPPREV=GRIDS(ID_GRD)%TRIGP(IE,IPREV)
            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(IE,I)
              IPNEXT=GRIDS(ID_GRD)%TRIGP(IE,INEXT)
              IPPREV=GRIDS(ID_GRD)%TRIGP(IE,IPREV)
              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

        ! Calculate area of cell
        ALLOCATE(GRID_AREA(MNP))
        DO I1=1,MNP
          GRID_AREA(I1) = 0.0
          NB=NBASSIGNEDCORNER(I1)
          DO I2=1,NB
            I3 = MOD(I2,NB)+1
            PT(1,1)=GRID_CENTER_LON(I1)
            PT(2,1)=GRID_CORNER_LON(I2,I1)
            PT(3,1)=GRID_CORNER_LON(I3,I1)
            PT(1,2)=GRID_CENTER_LAT(I1)
            PT(2,2)=GRID_CORNER_LAT(I2,I1)
            PT(3,2)=GRID_CORNER_LON(I3,I1)

            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

            GRID_AREA(I1) = GRID_AREA(I1) + ABS(THEDET)
          ENDDO
        ENDDO

        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

!/ ------------------------------------------------------------------- /
      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(IE,I)
          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(MNE,3)
        INTEGER, INTENT(INOUT)          :: IOBP(MNP)
        INTEGER, INTENT(INOUT)          :: NEIGHBOR_PREV(MNP)
        INTEGER, INTENT(INOUT)          :: NEIGHBOR_NEXT(MNP)
        
        INTEGER, POINTER :: STATUS(:)
        INTEGER, POINTER :: COLLECTED(:)
        INTEGER, POINTER :: NEXTVERT(:)
        INTEGER, POINTER :: PREVVERT(:)
        
        INTEGER :: IE, I, IP, IP2, IP3
        INTEGER :: ISFINISHED, INEXT, IPREV
        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(IE,I)
            IPNEXT=TRIGP(IE,INEXT)
            IPPREV=TRIGP(IE,IPREV)
            IF (STATUS(IP).EQ.0) THEN
              STATUS(IP)=1
              PREVVERT(IP)=IPPREV
              NEXTVERT(IP)=IPNEXT
            END IF
          END DO
        END DO
        STATUS(:)=0
        DO
          COLLECTED(:)=0
          DO IE=1,MNE
            DO I=1,3
              CALL TRIANG_INDEXES(I, INEXT, IPREV)
              IP=TRIGP(IE,I)
              IPNEXT=TRIGP(IE,INEXT)
              IPPREV=TRIGP(IE,IPREV)
              IF (STATUS(IP).EQ.0) THEN
                ZNEXT=NEXTVERT(IP)
                IF (ZNEXT.EQ.IPPREV) THEN
                  COLLECTED(IP)=1
                  NEXTVERT(IP)=IPNEXT
                  IF (NEXTVERT(IP).EQ.PREVVERT(IP)) THEN
                    STATUS(IP)=1
                  END IF
                END IF
              END IF
            END DO
          END DO

          ISFINISHED=1
          DO IP=1,MNP
            IF ((COLLECTED(IP).EQ.0).AND.(STATUS(IP).EQ.0)) THEN
              STATUS(IP)=-1
              NEIGHBOR_NEXT(IP)=NEXTVERT(IP)
            END IF
            IF (STATUS(IP).EQ.0) THEN
              ISFINISHED=0
            END IF
          END DO
          IF (ISFINISHED.EQ.1) THEN
            EXIT
          END IF
        END DO

!  Now computing the prev items
        STATUS = 0
        NEXTVERT = 0
        PREVVERT = 0
        DO IE=1,MNE
          DO I=1,3
            CALL TRIANG_INDEXES(I, INEXT, IPREV)
            IP=TRIGP(IE,I)
            IPNEXT=TRIGP(IE,INEXT)
            IPPREV=TRIGP(IE,IPREV)
            IF (STATUS(IP).EQ.0) THEN
              STATUS(IP)=1
              PREVVERT(IP)=IPPREV
              NEXTVERT(IP)=IPNEXT
            END IF
          END DO
        END DO
        STATUS(:)=0
        DO
          COLLECTED(:)=0
          DO IE=1,MNE
            DO I=1,3
              CALL TRIANG_INDEXES(I, INEXT, IPREV)
              IP=TRIGP(IE,I)
              IPNEXT=TRIGP(IE,INEXT)
              IPPREV=TRIGP(IE,IPREV)
              IF (STATUS(IP).EQ.0) THEN
                ZPREV=PREVVERT(IP)
                IF (ZPREV.EQ.IPNEXT) THEN
                  COLLECTED(IP)=1
                  PREVVERT(IP)=IPPREV
                  IF (PREVVERT(IP).EQ.NEXTVERT(IP)) THEN
                    STATUS(IP)=1
                  END IF
                END IF
              END IF
            END DO
          END DO

          ISFINISHED=1
          DO IP=1,MNP
            IF ((COLLECTED(IP).EQ.0).AND.(STATUS(IP).EQ.0)) THEN
              STATUS(IP)=-1
              NEIGHBOR_PREV(IP)=PREVVERT(IP)     ! new code
            END IF
            IF (STATUS(IP).EQ.0) THEN
              ISFINISHED=0
            END IF
          END DO
          IF (ISFINISHED.EQ.1) THEN
            EXIT
          END IF
        END DO
!  Now making checks
        DO IP=1,MNP
          IP2=NEIGHBOR_NEXT(IP)
          IF (IP2.GT.0) THEN
            IP3=NEIGHBOR_PREV(IP2)
            IF (ABS(IP3 - IP).GT.0) THEN
              WRITE(*,*) 'IP=', IP, ' IP2=', IP2, ' IP3=', IP3
              WRITE(*,*) 'We have a dramatic inconsistency'
              STOP '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
!/ ------------------------------------------------------------------- /
      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
