WW3/wmgridmd_8F90_source.html

#include "w3macros.h"

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

MODULE wmgridmd

  !/

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

  !/                  | WAVEWATCH III           NOAA/NCEP |

  !/                  |           H. L. Tolman            |

  !/                  |           W. E. Rogers            |

  !/                  |                        FORTRAN 90 |

  !/                  | Last update :         10-Dec-2014 |

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

  !/

  !/    28-Dec-2005 : Origination WMGLOW, WMGHGH, WMRSPC. ( version 3.08 )

  !/    09-Mar-2006 : Carry land mask in WMGHGH.          ( version 3.09 )

  !/    24-Apr-2006 : Origination WMGEQL.                 ( version 3.09 )

  !/    25-Jul-2006 : Point output grid in WMRSPC.        ( version 3.10 )

  !/    23-Dec-2006 : Adding group test in WMGEQL.        ( version 3.10 )

  !/    28-Dec-2006 : Simplify NIT for partial comm.      ( version 3.10 )

  !/    22-Jan-2007 : Add saving of NAVMAX in WMGEQL.     ( version 3.10 )

  !/    02-Feb-2007 : Setting FLAGST in WMGEQL.           ( version 3.10 )

  !/    07-Feb-2007 : Setting FLAGST in WMGHGH.           ( version 3.10 )

  !/    15-Feb-2007 : Tweaking MAPODI algorithm in WMGEQL.( version 3.10 )

  !/    11-Apr-2008 : Bug fix active edges WMGEQL.        ( version 3.13 )

  !/    14-Apr-2008 : Bug fix for global grids WMGEQL.    ( version 3.13 )

  !/    26-Mar-2009 : Adding test output !/T9 to WMGLOW.  ( version 3.14 )

  !/    20-May-2009 : Linking FLAGST and FLGHG1.          ( version 3.14 )

  !/    26-May-2009 : Fix erroneous cyclic upd in WMGHGH. ( version 3.14 )

  !/    29-May-2009 : Preparing distribution version.     ( version 3.14 )

  !/    30-Oct-2009 : Implement run-time grid selection.  ( version 3.14 )

  !/                  (W. E. Rogers & T. J. Campbell, NRL)

  !/    06-Dec-2010 : Change from GLOBAL (logical) to ICLOSE (integer) to

  !/                  specify index closure for a grid.   ( version 3.14 )

  !/                  (T. J. Campbell, NRL)

  !/    23-Dec-2010 : Fix HPFAC and HQFAC by including the COS(YGRD)

  !/                  factor with DXDP and DXDQ terms.    ( version 3.14 )

  !/                  (T. J. Campbell, NRL)

  !/    12-Mar-2012 : Use MPI_COMM_NULL in checks.        ( version 3.14 )

  !/    06-Jun-2012 : Porting bugfixes from 3.14 to 4.07  ( version 4.07 )

  !/    05-Sep-2012 : Implementation of UNGTYPE with SCRIP

  !/                       (Mathieu Dutour Sikiric, IRB; Aron Roland, Z&P)

  !/    21-Sep-2012 : Modify WMGHGH to support SCRIP remap( version 4.11 )

  !/                  write/read capabilities (K. Lind, NRL)

  !/    05-Aug-2013 : Change PR2/3 to UQ/UNO in distances.( version 4.12 )

  !/    10-Dec-2014 : Add checks for allocate status      ( version 5.04 )

  !/    28-Oct-2020 : Add SMCTYPE for SMC sub-grid.  JGLi ( version 7.13 )

  !/    26-Jan-2021 : Add WMSMCEQL for SMC sub-grid. JGLi ( version 7.13 )

  !/

  !/    Copyright 2009-2013 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

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

  !      WMGLOW    Subr. Public   Dependencies to lower ranked grids.

  !      WMGHGH    Subr. Public   Dependencies to higher ranked grids.

  !      WMGEQL    Subr. Public   Dependencies to same ranked grids.

  !      WMRSPC    Subr. Public   Make map of flags for spectral

  !                               conversion between grids.

  !      WMSMCEQL  Subr. Public   Dependencies on same ranked SMC grids.

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

  !

  !  4. Subroutines and functions used :

  !

  !      Name      Type  Module   Description

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

  !      W3SETO, W3SETG, W3DMO5, WMSETM

  !                Subr. W3xDATMD Manage data structures.

  !

  !      STRACE    Sur.  W3SERVMD Subroutine tracing.

  !      EXTCDE    Subr.   Id.    Program abort.

  !

  !      MPI_BCAST, MPI_BARRIER

  !                Subr. mpif.h   Comunication routines.

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

  !

  !  5. Remarks :

  !

  !      - WMGLOW and WMGHGH need to be run in this order to

  !        assure proper resolving of cross-dependencies.

  !      - WMGLOW and WMGEQL, idem.

  !

  !  6. Switches :

  !

  !     !/PRn  propagation scheme.

  !     !/UQ   propagation scheme.

  !     !/UNO  propagation scheme.

  !     !/SMC  Enable SMC sub-grids.

  !

  !     !/SHRD Distributed memory approach

  !     !/DIST

  !     !/MPI

  !

  !     !/O12  Removed boundary points output WMGEQL (central).

  !     !/O13  Removed boundary points output WMGEQL (edge).

  !

  !     !/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 wmglow ( FLRBPI )

    !/

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

    !/                  | WAVEWATCH III           NOAA/NCEP |

    !/                  |           H. L. Tolman            |

    !/                  |           W. E. Rogers            |

    !/                  |                        FORTRAN 90 |

    !/                  | Last update :         06-Jun-2018 !

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

    !/

    !/    06-Oct-2005 : Origination.                        ( version 3.08 )

    !/    10-Feb-2006 : Add test on grid resolution.        ( version 3.09 )

    !/    26-Mar-2009 : Adding test output !/T9.            ( version 3.14 )

    !/    30-Oct-2009 : Implement run-time grid selection.  ( version 3.14 )

    !/                  (W. E. Rogers & T. J. Campbell, NRL)

    !/    22-Dec-2010 : Adapt for use with irregular grids  ( version 3.14 )

    !/                  (W. E. Rogers, NRL)

    !/    12-Mar-2012 : Use MPI_COMM_NULL in checks.        ( version 4.07 )

    !/    06-Jun-2012 : Porting bugfixes from 3.14 to 4.07  ( version 4.07 )

    !/    10-Dec-2014 : Add checks for allocate status      ( version 5.04 )

    !/    06-Jun-2018 : Use W3PARALL                        ( version 6.04 )

    !/

    !  1. Purpose :

    !

    !     Determine relations to lower ranked grids for each grid.

    !     On the fly, the opposite relations are also saved.

    !

    !  2. Method :

    !

    !     Map active boundary points to lower ranked grids.

    !

    !  3. Parameters :

    !

    !     Parameter list

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

    !       FLRBPI  L.A.   O   Array with flags for external file use.

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

    !

    !  4. Subroutines used :

    !

    !      Name      Type  Module   Description

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

    !      W3SETO, W3SETG, W3DMO5

    !                Subr. W3xDATMD Manage data structures.

    !

    !      STRACE    Subr. W3SERVMD Subroutine tracing.

    !      EXTCDE    Subr.   Id.    Program abort.

    !

    !      MPI_BCAST, MPI_BARRIER

    !                Subr. mpif.h   Comunication routines.

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

    !

    !  5. Called by :

    !

    !      Name      Type  Module   Description

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

    !      WMINIT    Subr  WMINITMD Multi-grid model initialization.

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

    !

    !  6. Error messages :

    !

    !  7. Remarks :

    !

    !     - For MPI version it is assumed that NX, NY, NSEA, and NSEAL are

    !       properly initialized even if the grid is not run on the local

    !       process.

    !

    !  8. Structure :

    !

    !     See source code.

    !

    !  9. Switches :

    !

    !     !/MPI  Distribbuted memory management.

    !

    !     !/S    Enable subroutine tracing.

    !     !/T    Enable test output.

    !     !/T1   Test output for individual boundary points

    !     !/T2   Test output cross-reference table

    !     !/T9   Test output of map of boundary origine.

    !

    ! 10. Source code :

    !

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

    !

    USE w3servmd, ONLY: extcde

#ifdef W3_S

    USE w3servmd, ONLY: strace

#endif

    !

    USE w3gdatmd

    USE w3odatmd

    USE w3triamd

    USE wmmdatmd

    USE w3parall, ONLY : init_get_jsea_isproc

    !

    IMPLICIT NONE

    !

#ifdef W3_MPI

    include "mpif.h"

#endif

    !/

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

    !/ Parameter list

    !/

    LOGICAL, INTENT(OUT), OPTIONAL :: FLRBPI(NRGRD)

    !/

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

    !/ Local parameters

    !/

    INTEGER                 :: I, IBI, IX, IY, JS, J,       &

         JTOT, I1, J1, I2, J2

#ifdef W3_MPI

    INTEGER                 :: NXYG, IERR_MPI

#endif

#ifdef W3_S

    INTEGER, SAVE           :: IENT = 0

#endif

    INTEGER, ALLOCATABLE    :: TSTORE(:,:)

#ifdef W3_MPI

    LOGICAL                 :: FLBARR

#endif

    REAL                    :: XA, YA

    REAL                    :: FACTOR

    LOGICAL                 :: GRIDD(NRGRD,NRGRD) ! indicates grid-to-grid

    ! dependency

    LOGICAL                 :: RFILE(NRGRD),  FLAGOK

    LOGICAL                 :: INGRID ! indicates whether boundary point

    ! is in lower rank grid

    INTEGER                 :: IVER(4),JVER(4) ! (I,J) indices of vertices

    ! of cell (in lower rank grid J) enclosing

    ! boundary point (in higher rank grid I)

    REAL                    :: RW(4) ! Array of interpolation weights.

    INTEGER                 :: KVER ! counter for 4 vertices


    REAL                    :: DX_MIN_GRIDI,DY_MIN_GRIDI,DX_MAX_GRIDI, &

         DY_MAX_GRIDI

    REAL                    :: DX_MIN_GRIDJ,DY_MIN_GRIDJ,DX_MAX_GRIDJ, &

         DY_MAX_GRIDJ

    INTEGER  :: ITRI, IM1, IM2, IT, JT, ISFIRST, ITOUT, NBRELEVANT

    REAL :: DIST_MIN, DIST_MAX, EDIST

    LOGICAL RESOL_CHECK

    !

#ifdef W3_T9

    CHARACTER(LEN=1), ALLOCATABLE :: TMAP(:,:)

#endif

    !/

#ifdef W3_S

    CALL strace (ient, 'WMGLOW')

#endif

    !

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

    ! 1.  Test grid, Initialize and synchronize grids as needed ( !/MPI )

    !

#ifdef W3_MPI

    flbarr = .false.

#endif

    !

    DO i=1, nrgrd

      !

      IF ( .NOT. grids(i)%GINIT ) THEN

        IF ( improc .EQ. nmperr ) WRITE (mdse,1000) i

        CALL extcde ( 1000 )

      END IF


      CALL w3seto ( i, mdse, mdst )

      CALL w3setg ( i, mdse, mdst )

      !

#ifdef W3_MPI

      flbarr = flbarr .OR. mdatas(i)%FBCAST

      IF ( mdatas(i)%FBCAST .AND.                              &

           mdatas(i)%MPI_COMM_BCT.NE.mpi_comm_null ) THEN

        nxyg   = grids(i)%NX * grids(i)%NY

        CALL mpi_bcast ( grids(i)%MAPSTA(1,1), nxyg,        &

             mpi_integer, 0,                    &

             mdatas(i)%MPI_COMM_BCT, ierr_mpi )

        CALL mpi_bcast ( grids(i)%MAPST2(1,1), nxyg,        &

             mpi_integer, 0,                    &

             mdatas(i)%MPI_COMM_BCT, ierr_mpi )

        CALL mpi_bcast ( grids(i)%MAPFS (1,1), nxyg,        &

             mpi_integer, 0,                    &

             mdatas(i)%MPI_COMM_BCT, ierr_mpi )

        nxyg   = 3*grids(i)%NSEA

        CALL mpi_bcast ( grids(i)%MAPSF (1,1), nxyg,        &

             mpi_integer, 0,                    &

             mdatas(i)%MPI_COMM_BCT, ierr_mpi )

        CALL mpi_bcast ( grids(i)%CLATIS(1), nsea, mpi_real, 0,&

             mdatas(i)%MPI_COMM_BCT, ierr_mpi )

        CALL mpi_bcast ( sgrds(i)%SIG(0), nk+2, mpi_real, 0,&

             mdatas(i)%MPI_COMM_BCT, ierr_mpi )

      END IF

#endif

      !

    END DO

    !

#ifdef W3_MPI

    IF (flbarr) CALL mpi_barrier (mpi_comm_mwave,ierr_mpi)

#endif

    !

#ifdef W3_T

    WRITE (mdst,9010)

#endif

    !

#ifdef W3_SMC

    !!  Check GTYPE for all grids.

    IF( improc.EQ.nmperr )  WRITE(mdse,*) " GTYPES in WMGLOW:", &

         ( grids(i)%GTYPE, i=1, nrgrd )

#endif

    !

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

    ! 2.  Process grids

    !

    IF ( flagll ) THEN

      factor = 1.

    ELSE

      factor = 1.e-3

    END IF

    !

    gridd = .false.

    rfile = .false.

    !

    IF ( .NOT. ALLOCATED(nbi2g) ) THEN

      ALLOCATE ( nbi2g(nrgrd,nrgrd), stat=istat )

      check_alloc_status( istat )

    END IF

    nbi2g  = 0

    !

#ifdef W3_T

    WRITE (mdst,9020)

#endif

    !

    DO i=1, nrgrd

      !

#ifdef W3_T

      WRITE (mdst,9021) i, grank(i), outpts(i)%OUT5%NBI

#endif

      !

      ! 2.a Test for input boundary points

      !

      IF ( outpts(i)%OUT5%NBI .EQ. 0 ) THEN

#ifdef W3_T

        WRITE (mdst,9022) 'NO INPUT BOUNDARY POINTS, SKIPPING'

#endif

        cycle

      END IF

      !

      ! 2.b Test for lowest rank

      !

      IF ( grank(i) .EQ. 1 ) THEN

        rfile(i) = .true.

#ifdef W3_T

        WRITE (mdst,9022) 'RANK = 1, DATA FROM FILE'

#endif

        cycle

      END IF

      !

#ifdef W3_SMC

      !!  SMC grid only appears in same ranked group.  JGLi23Mar2021

      IF( grids(i)%GTYPE .EQ. smctype ) THEN

        IF( improc.EQ.nmperr ) WRITE(mdse,*) ' WMGLOW skip SMC grid', i

        cycle

      END IF

#endif

      !

      ! 2.c Search for input boundary points

      !


#ifdef W3_T

      WRITE (mdst,9022) 'SEARCHING FOR ACTIVE BOUNDARY POINTS'

#endif

      ibi    = 0

      !

      ! ... Set up data structure for grid

      !

      CALL w3seto ( i, mdse, mdst )

      CALL w3setg ( i, mdse, mdst )

      CALL w3dmo5 ( i, mdse, mdst, 1 )

      ALLOCATE ( tstore(nbi,0:4), stat=istat )

      check_alloc_status( istat )

      !

      ! ... Set up loop structure for grid

      !

      DO iy=1, ny

        DO ix=1, nx


          !notes : MAPSTA refers to GRIDS(I)%MAPSTA ...this is set in W3SETG

          IF ( abs(mapsta(iy,ix)) .EQ. 2 ) THEN

            xa     = real(xgrd(iy,ix)) !old code: X0 + REAL(IX-1)*SX

            ya     = real(ygrd(iy,ix)) !old code: Y0 + REAL(IY-1)*SY

            !

            ! ... Loop over previous (lower ranked) grids, going in order from highest

            !          of lower ranked grids (I-1) to lowest of lower ranked grids (1)

            !

            js     = 0

            !

            DO j=i-1, 1, -1

              !

              IF ( grank(j) .GE. grank(i) ) cycle

              !

#ifdef W3_SMC

              !!  SMC grid only suppots same ranked group so far.  JGLi12Apr2021

              IF( grids(j)%GTYPE .EQ. smctype ) THEN

                IF( improc.EQ.nmperr ) WRITE(mdse,*)   &

                     ' WMGLOW skip SMC grid', j

                cycle

              END IF

#endif

              !

              ! ... Check if in grid

              !

              ! notes:

              ! old version (v4.00):

              !        if in grid, return location in grid: a) JX, JY

              !                                   (lower left indices of cell),

              !                                             b) RX, RY

              !                                   (normalized location in cell)

              !        in not in grid, cycle (search next grid)

              ! new version (v4.01):

              !        Check if point within grid and compute interpolation weights using GSU

              !        if not in grid, cycle (search next grid)

              !

              IF (grids(j)%GTYPE .EQ. ungtype) THEN

                !AR: Here we need to take special care in the case that any problem occurs due to the XA, YA beeing 4 byte

                CALL is_in_ungrid(j, dble(xa), dble(ya), itout, iver, jver, rw)

                IF (itout.EQ.0) THEN

                  ingrid=.false.

                ELSE

                  ingrid=.true.

                  flagok =( abs(grids(j)%MAPSTA(jver(1),iver(1))).GE.1 .OR.  &

                       rw(1).LT.0.05 ) .AND.      &

                       ( abs(grids(j)%MAPSTA(jver(2),iver(2))).GE.1 .OR.  &

                       rw(2).LT.0.05 ) .AND.      &

                       ( abs(grids(j)%MAPSTA(jver(3),iver(3))).GE.1 .OR.  &

                       rw(3).LT.0.05 )

                END IF

                nbrelevant=3

              ELSE

                ingrid = w3grmp( grids(j)%GSU,  xa,  ya, iver , jver, rw )

                !          Print *, 'J=', J, 'IX=', IX, 'IY=', IY

                !          Print *, 'IN=', INGRID, 'XA=', XA, 'YA=', YA

                !          Print *, '    1: IVER=', IVER(1), 'JVER=', JVER(1), 'RW=', RW(1)

                !          Print *, '    2: IVER=', IVER(2), 'JVER=', JVER(2), 'RW=', RW(2)

                !          Print *, '    3: IVER=', IVER(3), 'JVER=', JVER(3), 'RW=', RW(3)

                !          Print *, '    4: IVER=', IVER(4), 'JVER=', JVER(4), 'RW=', RW(4)

                IF (ingrid) THEN

                  flagok =( abs(grids(j)%MAPSTA(jver(1),iver(1))).GE.1 .OR.  &

                       rw(1).LT.0.05 ) .AND.      &

                       ( abs(grids(j)%MAPSTA(jver(2),iver(2))).GE.1 .OR.  &

                       rw(2).LT.0.05 ) .AND.      &

                       ( abs(grids(j)%MAPSTA(jver(4),iver(4))).GE.1 .OR.  &

                       rw(4)   .LT.0.05 ) .AND.      &

                       ( abs(grids(j)%MAPSTA(jver(3),iver(3))).GE.1 .OR.  &

                       rw(3)   .LT.0.05 )

                END IF

                nbrelevant=4

              END IF

              !  internal name=                 GSU XTIN YTIN    IS     JS  RW   (notes)

              !  role=out                       in    in   in   out    out out

              !  size=                          ---    1    1     4      4   4

              !

              ! notes:

              !      - organization of IVER(4),JVER(4),RW(4) as returned by W3GRMP are

              !        as follows:

              !              Point 1 : lower i , lower j (JY1,JX1)

              !              Point 2 : upper i , lower j (JY1,JX2)

              !              Point 3 : upper i , upper j (JY2,JX2)

              !              Point 4 : lower i , upper j (JY2,JX1)

              !              (counter-clockwise starting from lower i, lower j)

              !

              !  ... if not in grid, warning message and cycle (search next grid)

              IF ( .NOT.ingrid ) THEN

#ifdef W3_T

                IF ( iaproc .EQ. naperr ) THEN

                  IF ( flagll ) THEN

                    WRITE (ndse,2000) xa, ya

                  ELSE

                    WRITE (ndse,2001) xa, ya

                  END IF

                END IF

#endif

                cycle

              END IF


              !

              ! ... Check against MAPSTA

              !


              ! Notes:

              ! Old code        | becomes | New code

              !-----------------| --------| -------

              ! (1.-RX)*(1.-RY) | becomes | RW(1)

              ! RX*(1.-RY)      | becomes | RW(2)

              ! (1.-RX)*RY      | becomes | RW(4)

              ! RX*RY           | becomes | RW(3)

              ! JX1             | becomes | IVER(1)

              ! JY1             | becomes | JVER(1)

              ! JX2             | becomes | IVER(3)

              ! JY2             | becomes | JVER(3)


              ! Notes:

              ! IVER(1)=IVER(4), IVER(2)=IVER(3)

              ! JVER(1)=JVER(2), JVER(3)=JVER(4)


              ! point 1:

              flagok = ( abs(grids(j)%MAPSTA(jver(1),iver(1))).GE.1 .OR.  &

                   rw(1).LT.0.05 ) .AND.      &

                   ! point 2:

                   ( abs(grids(j)%MAPSTA(jver(2),iver(2))).GE.1 .OR.  &

                   rw(2).LT.0.05 ) .AND.      &

                   ! point 4:

                   ( abs(grids(j)%MAPSTA(jver(4),iver(4))).GE.1 .OR.  &

                   rw(4)   .LT.0.05 ) .AND.      &

                   ! point 3:

                   ( abs(grids(j)%MAPSTA(jver(3),iver(3))).GE.1 .OR.  &

                   rw(3)   .LT.0.05 )

              !

              IF ( .NOT.flagok ) cycle

              !

              ! ... We found interpolation data !

              !

              js     = j

              ibi    = ibi + 1

              gridd(i,js) = .true.

              !

              xbpi(ibi)  = xa

              ybpi(ibi)  = ya

              isbpi(ibi) = mapfs(iy,ix)

              !

              tstore(ibi, 0) = js

              !

              ! notes:

              !   To maintain perfect consistency with old code, we would make code such that:

              !   -  point 1 in GSU goes to point 1 in RDBPI, TSTORE

              !   -  point 2 in GSU goes to point 2 in RDBPI, TSTORE

              !   -  point 4 in GSU goes to point 3 in RDBPI, TSTORE

              !   -  point 3 in GSU goes to point 4 in RDBPI, TSTORE

              !   Instead, here, we map point 4 in GSU goes to point 4 in RDBPI, TSTORE, etc.

              !   Thus the ordering of RDBPI, TSTORE has changed.

              !   I have no reason to believe that the ordering in RDBPI, TSTORE is important.

              !   I have gone through test case mww3_test_02 for gridsets a,b,c,d and found

              !      no change in result vs v4.00.


              DO kver=1,4

                IF (kver .LE. nbrelevant) THEN

                  IF ( abs(grids(j)%MAPSTA(jver(kver),iver(kver))).GE.1 &

                       .AND. rw(kver)   .GT.0.05 ) THEN

                    rdbpi(ibi,kver) = rw(kver)

                    tstore(ibi,kver) = grids(j)%MAPFS(jver(kver),iver(kver))

                  ELSE

                    rdbpi(ibi,kver) = 0.

                    tstore(ibi,kver) = 0

                  END IF

                ELSE

                  rdbpi(ibi,kver) = 0.

                  tstore(ibi,kver) = 0

                END IF


              END DO


              !

              ! .....normalize weights to give sum(R)=1

              rdbpi(ibi,:) = rdbpi(ibi,:) / sum(rdbpi(ibi,:))

              !

              ! Search was successful, so no need to search through other grids, so exit loop

              EXIT

            END DO ! "DO J=..."

            !

            IF ( js.EQ.0 .AND. improc.EQ.nmperr )                 &

                 WRITE (mdse,1020) i, ix, iy, xa, ya

            !

          END IF ! If a boundary point...


        END DO ! "DO IX=..."

      END DO ! "DO IY=..."


      !

      ! 2.d Error checks

      !

      IF ( ibi .EQ. 0 ) THEN

        rfile(i) = .true.

        IF ( improc .EQ. nmperr ) WRITE (mdse,1021)

        DEALLOCATE ( outpts(i)%OUT5%IPBPI, outpts(i)%OUT5%ISBPI,  &

             outpts(i)%OUT5%XBPI,  outpts(i)%OUT5%YBPI,   &

             outpts(i)%OUT5%RDBPI, stat=istat )

        check_dealloc_status( istat )

        cycle

      ELSE IF ( ibi .NE. outpts(i)%OUT5%NBI ) THEN

        CALL extcde ( 1020 )

      ENDIF

      !

      ! 2.e Sort spectra by grid, fill IPBPI, and get NBI2 and ....

      !


      ipbpi  = 0

      nbi2   = 0

      !

      DO j=1, nrgrd

        DO i1=1, nbi

          IF ( tstore(i1,0) .NE. j ) cycle

          DO j1=1, 4

            IF ( tstore(i1,j1).NE.0 .AND. ipbpi(i1,j1).EQ.0 ) THEN

              nbi2         = nbi2 + 1

              ipbpi(i1,j1) = nbi2

              DO i2=i1, nbi

                IF ( tstore(i2,0) .NE. j ) cycle

                DO j2=1, 4

                  IF ( tstore(i2,j2) .EQ. tstore(i1,j1) )         &

                       ipbpi(i2,j2) = nbi2

                END DO

              END DO

            END IF

          END DO

        END DO

      END DO

      !

      ! 2.f Set up spectral storage and cross-grid mapping

      !

      CALL w3dmo5 ( i, mdse, mdst, 3 )

      !

      ALLOCATE ( mdatas(i)%NBI2S(nbi2,2), stat=istat )

      check_alloc_status( istat )

      nbi2s  =>  mdatas(i)%NBI2S

      !

      DO i1=1, nbi

        DO j1=1, 4

          IF ( ipbpi(i1,j1) .NE. 0 ) THEN

            nbi2s(ipbpi(i1,j1),1) = tstore(i1,0)

            nbi2s(ipbpi(i1,j1),2) = tstore(i1,j1)

          END IF

        END DO

      END DO

      !

      DO i1=1, nbi2

        nbi2g(i,nbi2s(i1,1)) = nbi2g(i,nbi2s(i1,1)) + 1

      END DO

      !

      ! 2.g Test output

      !

#ifdef W3_T1

      WRITE (mdst,9023)

      DO j=1, nbi

        WRITE (mdst,9024) j, isbpi(j), factor*xbpi(j),          &

             factor*ybpi(j), ipbpi(j,:), rdbpi(j,:), tstore(j,:)

      END DO

#endif

      !

#ifdef W3_T2

      WRITE (mdst,9025)

      DO j=1, nbi2

        WRITE (mdst,9026) j, nbi2s(j,:)

      END DO

#endif

      !

#ifdef W3_T9

      ALLOCATE ( tmap(nx,ny), stat=istat )

      check_alloc_status( istat )

#endif

      !

#ifdef W3_T9

      DO ix=1, nx

        DO iy=1, ny

          IF ( abs(mapsta(iy,ix)) .EQ. 0 ) then

            tmap(ix,iy) = '/'

          ELSE IF ( abs(mapsta(iy,ix)) .EQ. 1 ) then

            tmap(ix,iy) = '-'

          ELSE IF ( abs(mapsta(iy,ix)) .EQ. 2 ) then

            tmap(ix,iy) = 'X'

          END IF

        END DO

      END DO

#endif

      !

#ifdef W3_T9

      DO j=1, nbi

        ix = mapsf(isbpi(j),1)

        iy = mapsf(isbpi(j),2)

        WRITE (tmap(ix,iy),'(I1)') tstore(j,0)

      END DO

#endif

      !

#ifdef W3_T9

      DO j=1, 1+(nx-1)/130

        WRITE (mdst,9029) i, j

        DO iy=ny, 1, -1

          i1 = j*130-129

          i2 = min( nx , j*130 )

          WRITE (mdst,'(1X,130A1)') tmap(i1:i2,iy)

        END DO

      END DO

#endif

      !

#ifdef W3_T9

      DEALLOCATE ( tmap, stat=istat )

      check_dealloc_status( istat )

#endif

      !

      DEALLOCATE ( tstore, stat=istat )

      check_dealloc_status( istat )

      !

    END DO

    !

#ifdef W3_T

    WRITE (mdst,9027)

    DO i=1, nrgrd

      WRITE (mdst,9028) outpts(i)%OUT5%NBI, outpts(i)%OUT5%NBI2, &

           rfile(i), nbi2g(i,:)

    END DO

#endif

    !

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

    ! 3.  Finalyze grid dependencies in GRDLOW

    ! 3.a Get size of array and dimension

    !


    ! notes:

    !    GRIDD(I,J) indicates whether grid I is dependent on lower ranked grid J

    !    JS counts the number of grids J that grid I is dependent on

    !    GRDLOW is sized to accomodate the grid with the largest JS


    jtot   = 0

    DO i=1, nrgrd

      js     = 0

      DO j=1, nrgrd

        IF ( gridd(i,j) ) js = js + 1

      END DO

      jtot   = max( jtot , js )

    END DO

    !

    IF ( ALLOCATED(grdlow) ) THEN

      DEALLOCATE ( grdlow, stat=istat )

      check_dealloc_status( istat )

    END IF

    ALLOCATE ( grdlow(nrgrd,0:jtot), stat=istat )

    check_alloc_status( istat )

    grdlow = 0

    !

#ifdef W3_T

    WRITE (mdst,9030) jtot

#endif

    !

    ! 3.b Fill array

    !

    flagok = .true.

    !

    DO i=1, nrgrd

      jtot   = 0

      DO j=1, nrgrd

        IF ( gridd(i,j) ) THEN

          jtot   = jtot + 1

          grdlow(i,jtot) = j

          ! ... error checking: catch situation where ranks are inconsistent with

          !                     resolution


          ! notes:

          ! old code:   SXJ=GRIDS(J)%SX

          !             SXI=GRIDS(I)%SX

          !             SYJ=GRIDS(J)%SY

          !             SYI=GRIDS(I)%SY

          !             also, old code did not need to check both min and max,

          !             since they were the same

          ! new code:

          !       SXI(:,:) ==> GRIDS(I)%HPFAC ! resolution in higher rank grid I

          !                    (approximate in case of irregular grids)

          !       SYI(:,:) ==> GRIDS(I)%HQFAC ! viz.

          !       SXJ(:,:) ==> GRIDS(J)%HPFAC ! resolution in lower rank grid J

          !                    (approximate in case of irregular grids)

          !       SYJ(:,:) ==> GRIDS(J)%HQFAC ! viz.


          ! notes:

          ! for irregular grids, we require

          !       1) smallest cell in low rank grid is larger than smallest cell

          !          in high rank grid

          !       2) largest cell in low rank grid is larger than largest cell

          !          in high rank grid

          ! Each dimension (along i/p and j/q axes) is checked separately,

          !     making 4 checks total.

          ! This is strict, and may generate "false positives" in error checking

          !     here. In this case, the user may wish to disable this error checking.

          ! For case of regular grids, we cannot use HPFAC, since it goes to zero

          !     at pole. We instead use good ol' SX and SY


          IF ( grids(i)%GTYPE .EQ. clgtype ) THEN

            dx_min_gridi=minval(grids(i)%HPFAC)

            dy_min_gridi=minval(grids(i)%HQFAC)

            dx_max_gridi=maxval(grids(i)%HPFAC)

            dy_max_gridi=maxval(grids(i)%HQFAC)

          ELSEIF ( grids(i)%GTYPE .EQ. rlgtype .OR.      &

               grids(i)%GTYPE .EQ. smctype ) THEN

            !!Li  SMC grid shares mesh with regular grid.  22Mar2021

            dx_min_gridi=grids(i)%SX

            dy_min_gridi=grids(i)%SY

            dx_max_gridi=grids(i)%SX

            dy_max_gridi=grids(i)%SY

          ELSEIF ( grids(i)%GTYPE .EQ. ungtype ) THEN

            isfirst=1

            dist_max=0

            dist_min=0

            DO itri=1,grids(i)%NTRI

              DO it=1,3

                IF (it.EQ.3) THEN

                  jt=1

                ELSE

                  jt=it+1

                END IF

                im1=grids(i)%TRIGP(it,itri)

                im2=grids(i)%TRIGP(jt,itri)

                edist=w3dist(flagll, real(grids(i)%XGRD(1,im1)), &

                     REAL(GRIDS(I)%YGRD(1,IM1)), REAL(GRIDS(I)%XGRD(1,IM2)), &

                     REAL(GRIDS(I)%YGRD(1,IM2)))

                IF (isfirst.EQ.1) THEN

                  dist_max=edist

                  dist_min=edist

                  isfirst=0

                ELSE

                  IF (edist.GT.dist_max) THEN

                    dist_max=edist

                  END IF

                  IF (edist.LT.dist_min) THEN

                    dist_min=edist

                  END IF

                END IF

              END DO

            END DO

            dx_min_gridi=dist_min

            dy_min_gridi=dist_min

            dx_max_gridi=dist_max

            dy_max_gridi=dist_max

          ELSE

            CALL extcde ( 601 )

          END IF


          IF ( grids(j)%GTYPE .EQ. clgtype ) THEN

            dx_min_gridj=minval(grids(j)%HPFAC)

            dy_min_gridj=minval(grids(j)%HQFAC)

            dx_max_gridj=maxval(grids(j)%HPFAC)

            dy_max_gridj=maxval(grids(j)%HQFAC)

          ELSEIF ( grids(j)%GTYPE .EQ. rlgtype .OR.      &

               grids(j)%GTYPE .EQ. smctype ) THEN

            !!Li  SMC grid shares mesh with regular grid.  22Mar2021

            dx_min_gridj=grids(j)%SX

            dy_min_gridj=grids(j)%SY

            dx_max_gridj=grids(j)%SX

            dy_max_gridj=grids(j)%SY

          ELSEIF ( grids(j)%GTYPE .EQ. ungtype ) THEN

            isfirst=1

            dist_max=0

            dist_min=0

            DO itri=1,grids(j)%NTRI

              DO it=1,3

                IF (it.EQ.3) THEN

                  jt=1

                ELSE

                  jt=it+1

                END IF

                im1=grids(j)%TRIGP(it,itri)

                im2=grids(j)%TRIGP(jt,itri)

                edist=w3dist(flagll, real(grids(j)%XGRD(1,im1)), &

                     REAL(GRIDS(J)%YGRD(1,IM1)), REAL(GRIDS(J)%XGRD(1,IM2)), &

                     REAL(GRIDS(J)%YGRD(1,IM2)))

                IF (isfirst.EQ.1) THEN

                  dist_max=edist

                  dist_min=edist

                  isfirst=0

                ELSE

                  IF (edist.GT.dist_max) THEN

                    dist_max=edist

                  END IF

                  IF (edist.LT.dist_min) THEN

                    dist_min=edist

                  END IF

                END IF

              END DO

            END DO

            dx_min_gridj=dist_min

            dy_min_gridj=dist_min

            dx_max_gridj=dist_max

            dy_max_gridj=dist_max

          ELSE

            CALL extcde ( 602 )

          END IF


          resol_check=.false.

#ifdef W3_T38

          resol_check=.true.

#endif

          IF (resol_check) THEN

            IF ( dx_min_gridj .LT. 0.99*dx_min_gridi .OR.    &

                 dy_min_gridj .LT. 0.99*dy_min_gridi .OR.    &

                 dx_max_gridj .LT. 0.99*dx_max_gridi .OR.    &

                 dy_max_gridj .LT. 0.99*dy_max_gridi  ) THEN

              print *, 'DX_MIN_GRID I=', dx_min_gridi, ' J=', dx_min_gridj

              print *, 'DX_MAX_GRID I=', dx_max_gridi, ' J=', dx_max_gridj

              IF ( improc.EQ.nmperr ) WRITE (mdse,1030)  &

                   j, grank(j), dx_min_gridj, dy_min_gridj, &

                   dx_max_gridj, dy_max_gridj, &

                   i, grank(i), dx_min_gridi, dy_min_gridi, &

                   dx_max_gridi, dy_max_gridi

              flagok = .false.

            END IF

          END IF


        END IF ! IF ( GRIDD(I,J) ) THEN


      END DO ! DO J=...

      grdlow(i,0) = jtot

    END DO ! DO I=...

    !

#ifdef W3_T

    WRITE (mdst,9031)

    DO i=1, nrgrd

      WRITE (mdst,9032) i, grdlow(i,0:grdlow(i,0))

    END DO

#endif

    !

    IF ( .NOT. flagok ) CALL extcde ( 1030 )

    !

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

    ! 4.  Finalyze grid dependencies in GRDHGH

    ! 4.a Get size of array and dimension

    !

    jtot   = 0

    DO i=1, nrgrd

      js     = 0

      DO j=1, nrgrd

        IF ( gridd(j,i) ) js = js + 1

      END DO

      jtot   = max( jtot , js )

    END DO

    !

    IF ( ALLOCATED(grdhgh) ) THEN

      DEALLOCATE ( grdhgh, stat=istat )

      check_dealloc_status( istat )

    END IF

    ALLOCATE ( grdhgh(nrgrd,0:jtot), stat=istat )

    check_alloc_status( istat )

    grdhgh = 0

    !

#ifdef W3_T

    WRITE (mdst,9040) jtot

#endif

    !

    ! 4.b Fill array

    !

    DO i=1, nrgrd ! low rank grid

      jtot   = 0

      DO j=1, nrgrd

        IF ( gridd(j,i) ) THEN ! grid j is of higher rank than grid i

          ! *and* there is dependency

          jtot   = jtot + 1  ! count the number of grids of higher

          ! rank than grid i

          grdhgh(i,jtot) = j ! save the grid number of the higher rank grid

        END IF

      END DO

      grdhgh(i,0) = jtot ! save the count of higher ranked grids

    END DO

    !

#ifdef W3_T

    WRITE (mdst,9041)

    DO i=1, nrgrd

      WRITE (mdst,9042) i, grdhgh(i,0:grdhgh(i,0))

    END DO

#endif

    !

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

    ! 5.  Export file flags

    !

    IF ( PRESENT(flrbpi) ) flrbpi = rfile

    !

    RETURN

    !

    ! Formats

    !

1000 FORMAT (/' *** WAVEWATCH III ERROR IN WMGLOW : *** '/        &

         '     GRID NOT INITIALIZED, GRID NR',i4 /)

    !

1020 FORMAT (/' *** WAVEWATCH III ERROR IN WMGLOW : *** '/       &

         '     CANNOT FIND SOURCE FOR BOUNDARY DATA '/          &

         '     GRID, IX, IY, X, Y:',3i6,2e12.4/)

    !

1021 FORMAT (/' *** WAVEWATCH III ERROR IN WMGLOW : *** '/       &

         '     NONE OF BOUNDARY POINTS CAN BE MAPPED'/          &

         '     READING FROM FILE INSTEAD'/)

    !

1030 FORMAT (/' *** WAVEWATCH III ERROR IN WMGLOW : *** '/       &

         '     RANKS AND RESOLUTIONS INCONSISTENT'/             &

         '        GRID',i4,' RANK',i4,' RESOLUTION :',4e10.3/   &

         '        GRID',i4,' RANK',i4,' RESOLUTION :',4e10.3/)

    !

2000 FORMAT (/' *** WAVEWATCH-III WARNING : BOUNDARY POINT'/  &

         '     NOT FOUND IN LOWER RANK GRID : ',2f10.3/  &

         '     POINT SKIPPED '/)

    !

2001 FORMAT (/' *** WAVEWATCH-III WARNING : BOUNDARY POINT'/  &

         '     NOT FOUND IN LOWER RANK GRID : ',2e10.3/  &

         '     POINT SKIPPED '/)

    !

#ifdef W3_T

9010 FORMAT ( ' TEST WMGLOW : ALL GRIDS INITIALIZED')

#endif

    !

#ifdef W3_T

9020 FORMAT ( ' TEST WMGLOW : STARTING LOOP OVER GRIDS')

9021 FORMAT ( ' TEST WMGLOW : I, RANK, NBI :',2i4,i6)

9022 FORMAT ( '               ',a)

#endif

#ifdef W3_T1

9023 FORMAT (' TEST WMGLOW : POINT DATA ')

9024 FORMAT (i5,i8,2f6.1,4i5,4f5.2,i3,4i8)

#endif

#ifdef W3_T2

9025 FORMAT (' TEST WMGLOW : NBI2S ')

9026 FORMAT ('           ',2i4,2x,i8)

#endif

#ifdef W3_T

9027 FORMAT (' TEST WMGLOW : NBI, NBI2, RFILE, NBI2G ')

9028 FORMAT ('               ',2i5,l2,' : ',20i5)

#endif

#ifdef W3_T9

9029 FORMAT (' TEST WMGLOW : SOURCE MAP GRID',i3,'   PART',i3)

#endif

    !

#ifdef W3_T

9030 FORMAT ( ' TEST WMGLOW : GRDLOW DIMENSIONED AT ',i2)

9031 FORMAT ( ' TEST WMGLOW : GRDLOW :')

9032 FORMAT ( '                 ',2i4,' : ',20i3)

#endif

    !

#ifdef W3_T

9040 FORMAT ( ' TEST WMGLOW : GRDHGH DIMENSIONED AT ',i2)

9041 FORMAT ( ' TEST WMGLOW : GRDHGH :')

9042 FORMAT ( '                 ',2i4,' : ',20i3)

#endif

    !/

    !/ End of WMGLOW ----------------------------------------------------- /

    !/

  END SUBROUTINE wmglow


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

  SUBROUTINE wmghgh

    !/

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

    !/                  | WAVEWATCH III           NOAA/NCEP |

    !/                  |           H. L. Tolman            |

    !/                  |           W. E. Rogers            |

    !/                  |                        FORTRAN 90 |

    !/                  | Last update :         10-Dec-2014 !

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

    !/

    !/    28-Dec-2005 : Origination.                        ( version 3.08 )

    !/    09-Mar-2006 : Carry over land mask.               ( version 3.09 )

    !/    28-Dec-2006 : Simplify NIT for partial comm.      ( version 3.10 )

    !/    07-Feb-2007 : Setting FLAGST.                     ( version 3.10 )

    !/    20-May-2009 : Linking FLAGST and FLGHG1.          ( version 3.14 )

    !/    26-May-2009 : Fix erroneous cyclic updating.      ( version 3.14 )

    !/    30-Oct-2009 : Implement run-time grid selection.  ( version 3.14 )

    !/                  (W. E. Rogers & T. J. Campbell, NRL)

    !/    06-Dec-2010 : Change from GLOBAL (logical) to ICLOSE (integer) to

    !/                  specify index closure for a grid.   ( version 3.14 )

    !/                  (T. J. Campbell, NRL)

    !/    23-Dec-2010 : Fix HPFAC and HQFAC by including the COS(YGRD)

    !/                  factor with DXDP and DXDQ terms.    ( version 3.14 )

    !/                  (T. J. Campbell, NRL)

    !/    07-Jul-2011 : Bug fix for IX bounds with wrapping ( version 3.14+)

    !/                  grids (see use of "IDSTLA" below)

    !/                  (W. E. Rogers, NRL)

    !/    02-Aug-2011 : Adapted for use with irregular      ( version 3.14+)

    !/                  grids (W. E. Rogers, NRL)

    !/    21-Sep-2012 : Modified to implement SCRIP remap   ( version 4.11 )

    !/                  file read and write option

    !/                  (K. R. Lind, NRL)

    !/    05-Aug-2013 : Change PR2/3 to UQ/UNO in distances.( version 4.12 )

    !/    10-Dec-2014 : Add checks for allocate status      ( version 5.04 )

    !/    20-Jan-2017 : Fix SCRIP ALLWGTS allocation error and improve

    !/                  SCRIPNC SCRIP_STOP report and exit. ( version 6.02 )

    !/

    !  1. Purpose :

    !

    !     Determine relation to higher ranked grids for each grid.

    !     Base map set in WMGLOW, supplemental data computed here.

    !

    !  2. Method :

    !

    !     Map averaging information for higher ranked grid to lower

    !     ranked grid.

    !

    !  3. Parameters :

    !

    !  4. Subroutines used :

    !

    !      Name      Type  Module   Description

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

    !      W3SETO, W3SETG, W3DMO5, WMSETM

    !                Subr. W3xDATMD Manage data structures.

    !      STRACE    Sur.  W3SERVMD Subroutine tracing.

    !      EXTCDE    Sur.    Id.    Program abort.

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

    !

    !  5. Called by :

    !

    !  6. Error messages :

    !

    !  7. Remarks :

    !

    !     Regarding the map of distances to the boundary :

    !      - v4.00 : the map of distances to the boundary was intentionally

    !                not an accurate characteristic distance. It was felt that

    !                it was more important that it be 'safe' and quick to compute.

    !                An iterative method was used to compute distance by starting

    !                at boundary and working inwards one grid row layer at a time,

    !                incrementing distance by dx etc. until the distance map was

    !                filled in. This was characterized as "local increment solution

    !                only."

    !      - v4.01 : conversion to work with irregular grids. Author could not

    !                think of any way to retain "local increment solution" method

    !                for situation of irregular grids. Therefore method has been

    !                changed to compute accurate distances. New method is also

    !                more transparent and simpler with much less code, thus

    !                easier to modify or debug. It is expected that this method

    !                could be more expensive to compute. Isolated timings were

    !                not performed. Since the iteration step has been removed,

    !                it is hoped that the expense is at least offset somewhat.

    !

    !     Regarding method of calculating weights :

    !              o If SCRIP software is not compiled into WW3 by user

    !                (i.e. if SCRIP switch is not set, then original method

    !                (denoted "_OM") will be used.

    !              o If SCRIP is activated by user, and all grids are

    !                regular and specified in terms of meters (cartesian),

    !                then WMGHGH will calculate weights using both methods,

    !                and then compare the two, producing an error message

    !                if they do not match (built-in regression testing)

    !                For more info, see Section 0a below.

    !

    !     re: Inconsistent RANK vs NBI (warning message) in Section 1.d :

    !         This was an error, but has been changed to a warning to allow

    !         more flexibility, e.g. having two outer grids with different

    !         rank (latter to avoid handling via WMGEQL).

    !         Change made July 2016.

    !         Old system:

    !           * grid rank > 1 and NBI>0 : do computations

    !           * grid rank > 1 and NBI=0 : error message

    !           * grid rank = 1 and NBI>0 : do nothing

    !           * grid rank = 1 and NBI=0 : do nothing

    !         New system:

    !           * grid rank > 1 and NBI>0 : do computations

    !           * grid rank > 1 and NBI=0 : do nothing w/ warning message

    !           * grid rank = 1 and NBI>0 : do nothing

    !           * grid rank = 1 and NBI=0 : do nothing

    !

    !  8. Structure :

    !

    !     See source code.

    !

    !  9. Switches :

    !

    !     !/SHRD Distributed memory approach

    !     !/DIST

    !

    !     !/PRn  propagation scheme.

    !

    !     !/S    Enable subroutine tracing.

    !     !/T    Enable test output.

    !     !/T3   Test output for received spectra.

    !     !/T4   Test output for sent spectra.

    !

    ! 10. Source code :

    !

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

    !

    USE constants

    USE w3servmd, ONLY: extcde

    USE w3gsrumd, ONLY: w3dist

#ifdef W3_S

    USE w3servmd, ONLY: strace

#endif

    !

    USE w3gdatmd

    USE w3odatmd

    USE wmmdatmd

    USE w3parall, ONLY : init_get_jsea_isproc

    !      USE W3PARALL, ONLY : INIT_GET_JSEA_ISPROC_GLOB

#ifdef W3_SCRIP

    USE wmscrpmd

    USE scrip_interface

#endif

    !/

    IMPLICIT NONE

    !

#ifdef W3_MPI

    include "mpif.h"

#endif

    !

    !/

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

    !/ Parameter list

    !/

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

    !/ Local parameters

    !/


    ! notes re: variable names: During the extension for irregular grids,

    !    some variable were renamed to make the code more readable:

    !         JX==> ISRC

    !         JY==> JSRC

    !         IX==> IDST

    !         IY==> JDST

    !         grid I ==> grid GDST

    !         grid J ==> grid GSRC


    INTEGER                 :: GDST, IJ, IDST, JDST, GSRC, JJ, IB, ISEA,    &

         JSEA, IDSTLA, IDSTHA, JDSTLA, JDSTHA,        &

         ISRC, JSRC, ISRCL, ISRCH, JSRCL, JSRCH, NIT, &

         NRTOT, NROK, JF, JR, NLMAX, ISPROC, ISPRO2,  &

         IREC, ISND, ITMP,ILOC

#ifdef W3_SCRIP

    INTEGER                 ::  NLMAX_SCRIP

#endif


#ifdef W3_DIST

    INTEGER                 :: LTAG0

#endif

#ifdef W3_MPI

    INTEGER                 :: IERR_MPI

#endif

#ifdef W3_S

    INTEGER, SAVE           :: IENT = 0

#endif


    INTEGER, ALLOCATABLE    :: IDSTL(:), IDSTH(:), JDSTL(:), JDSTH(:),      &

         MAPTST(:,:),                                 &

         I1(:,:), I2(:,:), I3(:), I4(:),              &

         INFLND(:,:)

    INTEGER, ALLOCATABLE    :: NX_BEG(:), NX_END(:)

#ifdef W3_MPIBDI

    INTEGER, ALLOCATABLE    :: NX_SIZE(:), IRQ(:), MSTAT(:,:)

#endif

#ifdef W3_MPI

    INTEGER                 :: IM, NX_REM, TAG, NRQ

#endif


    INTEGER, ALLOCATABLE    :: TMPINT_OM(:,:),TMPINT(:,:)

    REAL, ALLOCATABLE       :: TMPRL_OM(:,:) ,TMPRL(:,:)

    REAL, ALLOCATABLE       :: BDIST_OM(:)   ,BDIST(:)

    INTEGER                 :: NR0   , NR1   , NR2   , NRL   , NLOC

    INTEGER                 :: NR0_OM, NR1_OM, NR2_OM, NRL_OM, NLOC_OM


#ifdef W3_DIST

    INTEGER, ALLOCATABLE    :: LTAG(:)

#endif


    REAL                    :: FACTOR, STX, STY, STXY, NEWVAL,      &

         XL, XH, YL, YH, XA, YA, DXC, JD,     &

         WX, WY, WTOT


    LOGICAL                 :: FLGREC


    LOGICAL, ALLOCATABLE    :: GRIDOK(:),                           &

         STMASK(:,:), MASKI(:,:), TMPLOG(:)


    INTEGER                 :: JBND,IBND ! counter for boundary points

    REAL                    :: DD  ! distance to boundary point

    ! (temporary variable)

    REAL                    :: XDST,YDST

    REAL                    :: XSRC,YSRC

    REAL                    :: WXWY

    INTEGER                 :: NJDST,NIDST,KDST

    INTEGER                 :: NJSRC,NISRC,KSRC

    INTEGER                 :: IPNT,ICOUNT,IPNT2

    INTEGER                 :: DST_GRID_SIZE,ISTOP,JTMP


    REAL                    :: DX_MAX_GDST,DY_MAX_GDST

    REAL                    :: DX_MIN_GSRC,DY_MIN_GSRC


#ifdef W3_SCRIP

    TYPE allwgt

      TYPE(weight_data), POINTER :: WGTDATA(:)

    END TYPE allwgt

    TYPE(allwgt), ALLOCATABLE :: ALLWGTS(:)

    LOGICAL              :: L_MASTER = .true.

    LOGICAL              :: L_READ   = .false.

    LOGICAL              :: L_WRITE = .false.

#endif

#ifdef W3_SCRIPNC

    INTEGER              :: IMPROC_ASSIGN

    CHARACTER(LEN=80)    :: interp_file1, interp_file_test

    CHARACTER(LEN=3)     :: cdst, csrc

    LOGICAL, ALLOCATABLE :: LGRDREAD(:,:)

    LOGICAL, ALLOCATABLE :: LGRDWRITE(:,:)

    INTEGER              :: NGRDRANK(2)

#endif


    LOGICAL                 :: LSCRIP=.false.      ! true if SCRIP switch is set,

    ! indicates that SCRIP code has

    ! been compiled into WW3

    LOGICAL                 :: LSCRIPNC=.false.    ! true if SCRIPNC switch is set,

    ! indicates that SCRIP code has

    ! been compiled with netCDF

    ! into WW3

    LOGICAL                 :: L_STOP = .false.    ! true if SCRIPNC switch is set

    ! and STOP_SCRIP file exists

    LOGICAL                 :: T38=.false.         ! true if T38 switch is set.

    ! This logical is necessary

    ! since it isn't possible to

    ! have two switches disabling

    ! the same line of code.

    LOGICAL                 :: ALL_REGULAR=.true.  ! true if all grids are

    ! regular grids

    LOGICAL                 :: DO_CHECKING=.false. ! true if we will be

    ! checking "old method" of

    ! computing weights vs.

    ! SCRIP method of computing

    ! weights.

    LOGICAL                 :: OLD_METHOD=.false.  ! true if we will compute

    ! using "old method" (does

    ! not necessarily mean

    ! that this solution is

    ! utilized)

    LOGICAL                 :: LMPIBDI=.false.     ! true if MPIBDI switch is set

    LOGICAL                 :: CALLED_SCRIP=.false.! true if SCRIP has been

    ! called for this processor


    INTEGER  :: ITRI, IM1, IM2, IT, JT, IsFirst

    REAL :: DIST_MIN, DIST_MAX, eDist


#ifdef W3_T

    CHARACTER(LEN=1), ALLOCATABLE :: MAPST(:,:)

#endif

    !/

#ifdef W3_T38

    CHARACTER (LEN=10) :: CDATE_TIME(3)

    INTEGER            :: DATE_TIME(8)

    INTEGER            :: ELAPSED_TIME, BEG_TIME(10), END_TIME

    INTEGER            :: NMYOUT=42

    CHARACTER (LEN=14) :: CMYOUT="myout00000.lis"

    CHARACTER (LEN=5)  :: CRANK

#endif


#ifdef W3_T38

    WRITE(crank, "(I5.5)") improc-1

    cmyout(6:10) = crank(1:5)

    OPEN (nmyout, file=cmyout, status="REPLACE")

#endif

#ifdef W3_S

    CALL strace (ient, 'WMGHGH')

#endif

    !

#ifdef W3_MPI

    CALL mpi_barrier(mpi_comm_mwave, ierr_mpi)

#endif

#ifdef W3_T38

    CALL date_and_time ( cdate_time(1), cdate_time(2), cdate_time(3), date_time)

    beg_time(1) = ((date_time(5)*60 + date_time(6))*60 + date_time(7))*1000 + date_time(8)

    WRITE(nmyout,*) "WMGHGH: START: 0 MSEC"

#endif


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

    ! 0.  Initializations / tests

    !

    IF ( .NOT. ALLOCATED(grdhgh) ) THEN

      IF ( improc.EQ.nmperr ) WRITE(mdse,1000)

      CALL extcde (1000)

    END IF


#ifdef W3_MPIBDI

    lmpibdi=.true.

#endif

#ifdef W3_SCRIP

    IF (improc .EQ. 1) THEN

      l_master = .true.

      l_write  = .true.

    ELSE

      l_master = .false.

      l_write  = .false.

    ENDIF

#endif

#ifdef W3_SCRIPNC

    INQUIRE(file="SCRIP_STOP", exist=l_stop)

    improc_assign = 1

#endif

    !

    !KRL  Allocate helper arrays to enable bottleneck loop parallelization

    ALLOCATE ( nx_beg(nmproc), nx_end(nmproc), stat=istat )

    check_alloc_status( istat )

#ifdef W3_MPIBDI

    ALLOCATE ( nx_size(nmproc), irq(2*nmproc), &

         mstat(mpi_status_size,2*nmproc), stat=istat )

    check_alloc_status( istat )

#endif

    !

    !!HT:

    !!HT: Set up and initialize storage data structures ....

    !!HT:

#ifdef W3_T38

    CALL date_and_time ( cdate_time(1), cdate_time(2), cdate_time(3), date_time)

    beg_time(2) = ((date_time(5)*60 + date_time(6))*60 + date_time(7))*1000 + date_time(8)

#endif

    DO gdst=1, nrgrd

      DO gsrc=1, nrgrd

        IF ( hgstge(gdst,gsrc)%INIT ) THEN

          IF ( hgstge(gdst,gsrc)%NREC .NE. 0 ) THEN

            DEALLOCATE (                                         &

                 hgstge(gdst,gsrc)%LJSEA , hgstge(gdst,gsrc)%NRAVG, &

                 hgstge(gdst,gsrc)%IMPSRC, hgstge(gdst,gsrc)%ITAG , &

                 hgstge(gdst,gsrc)%WGTH  , hgstge(gdst,gsrc)%SHGH , &

                 stat=istat )

            check_dealloc_status( istat )

          END IF

          IF ( hgstge(gdst,gsrc)%NSND .NE. 0 ) THEN

            DEALLOCATE (                                         &

                 hgstge(gdst,gsrc)%ISEND ,                          &

                 stat=istat )

            check_dealloc_status( istat )

          END IF

          hgstge(gdst,gsrc)%NTOT = 0

          hgstge(gdst,gsrc)%NREC = 0

          hgstge(gdst,gsrc)%NRC1 = 0

          hgstge(gdst,gsrc)%NSND = 0

          hgstge(gdst,gsrc)%NSN1 = 0

          hgstge(gdst,gsrc)%NSMX = 0

          hgstge(gdst,gsrc)%INIT = .false.

        END IF

      END DO

    END DO

    gdst=-999 ! unset grid

    gsrc=-999 ! unset grid


#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(2)

    WRITE(nmyout,*) "WMGHGH, LOOP 1 TOOK ", elapsed_time, " MSEC"

#endif


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

    ! 0.a Plan future behavior by setting logical variables.


#ifdef W3_SCRIP

    lscrip=.true.

#endif

#ifdef W3_SCRIPNC

    lscripnc=.true.

#endif

#ifdef W3_T38

    t38=.true.

#endif


    DO gdst=1, nrgrd

      IF ( grids(gdst)%GTYPE .NE. rlgtype .AND.        &

           grids(gdst)%GTYPE .NE. smctype ) THEN

        !!Li  Add SMCTYPE option into ALL_REGULAR case.  JGLi20Nov2020

        all_regular=.false.

      END IF

    END DO


    ! Notes re: FLAGLL case: Old method calculates overlap area based on deg lat

    !       and deg lon. New method (SCRIP) calculates overlap area based on real

    !       distances. Therefore weights will not match for FLAGLL case, so we

    !       do not perform checking for FLAGLL case.


    IF ( (.NOT.flagll) .AND. all_regular .AND. lscrip ) THEN

      IF ( improc.EQ.nmperr ) &

           WRITE (mdse,'(/2A)')'We will check SCRIP calculations ', &

           'against old method of calculating weights.'

      do_checking=.true.

    END IF


    IF (do_checking .OR. (.NOT.lscrip)) old_method=.true.


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

    ! 0.b Check solution method


    IF ( (.NOT.lscrip) .AND. (.NOT.all_regular) .AND. &

         (nrgrd.GT.1) ) THEN

      IF ( improc.EQ.nmperr ) &

           WRITE (mdse,'(/3A)') ' *** ERROR WMGHGH: ', &

           'IRREGULAR or UNSTRUCTURED grid detected: this requires ', &

           'SCRIP switch.'

      CALL extcde ( 999 )

    END IF


    !

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

    ! 1.  Set boundary distance maps

    ! 1.a Check if needed

    !

    !!HT: FLGBDI is a flag set in WMMDATMD to .FALSE. and is used to identify

    !!HT:        if the boundary distance maps have been initialized

    !!HT:

    !!HT: For each individual grid a map is generated identifying the distance

    !!HT: to open boundaries (MAPBDI, saved in structure MDATA in WMMDATMD).

    !!HT: This map is used later to choose if more that 1 high-res grids

    !!HT: could provide data to a low-res grid. The high-res grid with data

    !!HT: furthest away from its own open boundary will be used.


#ifdef W3_SCRIPNC

    IF (.NOT. l_stop) THEN   ! Do not need MAPBDI if going to stop after generating mappings

#endif

      IF ( .NOT. flgbdi ) THEN

        !

        IF ( flagll ) THEN

          factor = radius * dera

          !notes: was FACTOR = RADIUS / 360. (bug fix)

        ELSE

          factor = 1.

        END IF

        !

#ifdef W3_T

        WRITE (mdst,9010)

#endif

        !

        ! 1.b Loop over grids

        !

#ifdef W3_T38

        CALL date_and_time ( cdate_time(1), cdate_time(2), cdate_time(3), date_time)

        beg_time(3) = ((date_time(5)*60 + date_time(6))*60 + date_time(7))*1000 + date_time(8)

        elapsed_time = beg_time(3) - beg_time(1)

        WRITE(nmyout,*) "WMGHGH, BEGINNING BOTTLENECK LOOP AT ", elapsed_time, " MSEC"

#endif

        DO gdst=1, nrgrd


#ifdef W3_T38

          IF(improc.EQ.nmperr)WRITE(mdse,*)'GDST = ',gdst,' OUT OF ',nrgrd

#endif


          CALL w3seto ( gdst, mdse, mdst )

          CALL w3setg ( gdst, mdse, mdst )

          CALL wmsetm ( gdst, mdse, mdst )


          !            IF ( GTYPE .EQ. UNGTYPE ) THEN

          !               IF ( IMPROC.EQ.NMPERR ) &

          !                    WRITE (MDSE,'(/2A)') ' *** ERROR WMGHGH: ', &

          !                    'UNSTRUCTURED GRID SUPPORT NOT YET IMPLEMENTED ***'

          !               CALL EXTCDE ( 999 )

          !            END IF


          !

#ifdef W3_T

          WRITE (mdst,9011) gdst, grank(gdst), nbi

#endif

          !

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

          ! Inconsistent RANK vs NBI (warning message)

          !   This was an error, now changed to a warning (see notes section)

          IF ( (grank(gdst).NE.1) .AND. (nbi.EQ.0) ) THEN

            IF ( improc.EQ.nmperr ) &

                 WRITE (mdse,'(/2A)') ' WARNING in WMGHGH:   ', &

                 'NBI=0 AND RANK > 1 '

          END IF


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

          ! 1.c NBI=0, so computations not needed (test output only)

          !

          IF ( (nbi.EQ.0) .OR. (grank(gdst).EQ.1) ) THEN

            ! (then do nothing except test output)

#ifdef W3_T

            WRITE (mdst,9012)

#endif


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

            ! 1.d NBI>0, Generate map with distances to boundary.


            !!HT: Initialize MAPBDI

            !!HT:   0. for active boundary points

            !!HT:  -1. for points that are not considered at all (rescaled for test

            !!HT:      output only, only negative value is essentially later).

            !!HT:  -2. for points that still need to be processed.


          ELSE


            IF(improc.EQ.nmperr)WRITE(mdse,'(A)') &

                 '  Generating map with distances to boundary.'

            !  for purposes of screen output, would be useful to wait for other processors to catch up here...(if mpibdi switch used)

            ALLOCATE ( mdatas(gdst)%MAPBDI(ny,nx), stat=istat )

            check_alloc_status( istat )

            mapbdi => mdatas(gdst)%MAPBDI

            !

            !KRL           Set up ranges for X. If not MPIBDI, just 1 to NX

            nx_beg(improc) = 1

            nx_end(improc) = nx

#ifdef W3_MPIBDI

            nx_beg(1) = 1

            IF ( nmproc .EQ. 1 ) THEN

              nx_end(1) = nx

              nx_size(1) = nx

            ELSE

              nx_rem = mod( nx, nmproc )

              nx_size(1) = nx / nmproc

              IF (nx_rem .GT. 0) nx_size(1) = nx_size(1) + 1

              nx_end(1) = nx_beg(1) + nx_size(1) - 1

              DO im = 2, nmproc

                nx_beg(im) = nx_end(im-1) + 1

                nx_size(im) = nx / nmproc

                IF (im .LE. nx_rem) nx_size(im) = nx_size(im) + 1

                nx_end(im) = nx_beg(im) + nx_size(im) - 1

                nx_size(im-1) = nx_size(im-1) * ny

              END DO

              nx_size(nmproc) = nx_size(nmproc) * ny

            END IF

#endif

            !KRL           Setup complete

            !

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

            ! Loop to determine MAPBDI

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


            IF(improc.EQ.nmperr)WRITE(mdse,'(A)') &

                 'Starting MAPBDI 1st loop.'


            DO idst=nx_beg(improc), nx_end(improc)

              IF(mod(idst,250).EQ.0)THEN

                IF(lmpibdi)THEN

                  WRITE(mdse,'(4x,3(A,I5))')&

                       'processing column ',idst,' out of ',nx, &

                       ' on processor ',improc

                ELSEIF(improc.EQ.nmperr)THEN

                  WRITE(mdse,'(4x,2(A,I5))')&

                       'processing column ',idst,' out of ',nx

                ENDIF

              ENDIF

              DO jdst=1, ny

                IF ( mapsta(jdst,idst) .EQ. 0 ) THEN ! (excluded point)

                  mapbdi(jdst,idst) = -1. / sig(1) * dtmax ! new (bug fix)

                ELSE IF ( abs(mapsta(jdst,idst)) .EQ. 2 ) THEN

                  ! (boundary point)

                  mapbdi(jdst,idst) =  0.

                ELSE ! ABS(MAPSTA)=1 (sea point)

                  mapbdi(jdst,idst) = 1.0e+10

                ENDIF ! IF MAPSTA

              END DO ! DO JDST...

            END DO ! DO IDST...


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


            IF(improc.EQ.nmperr)WRITE(mdse,'(A)') &

                 'Starting MAPBDI 2nd loop.'


            DO ibnd=1,nx

              IF ( (mod(ibnd,25).EQ.0) .AND. &

                   (improc.EQ.nmperr) ) THEN

                WRITE(mdse,'(4x,2(A,I5))') &

                     'bnd. point ',ibnd,' out of ',nx

              ENDIF

              DO jbnd=1,ny

                IF ( abs(mapsta(jbnd,ibnd)) .EQ. 2 ) THEN

                  ! (boundary point)

#ifdef W3_OMPH

                  !$OMP PARALLEL DO PRIVATE(IDST,JDST,DD),SCHEDULE(DYNAMIC)

#endif

                  DO idst=nx_beg(improc), nx_end(improc)

                    DO jdst=1, ny

                      IF (abs(mapsta(jdst,idst)) .EQ. 1) THEN

                        !....find distance to this boundary point.

                        dd=factor*w3dist(flagll,real(xgrd(jdst,idst)), &

                             REAL(YGRD(JDST,IDST)),REAL(XGRD(JBND,IBND)), &

                             REAL(YGRD(JBND,IBND)))


                        ! Notes: The origin of "0.58 * GRAV" is to translate from distance (in meters)

                        ! to time (in seconds) required for a wave to travel from the boundary to point

                        ! JDST,IDST based on a specific group velocity 0.58*grav would be the group

                        ! velocity of a 7.3 s wave in deep water. Significance of T=7.3 s is explained

                        ! in notes by HT below.


                        dd=dd/ ( 0.58 * grav )

                        mapbdi(jdst,idst)=min(mapbdi(jdst,idst),dd)

                      ENDIF

                    END DO ! DO JDST

                  END DO ! DO IDST

#ifdef W3_OMPH

                  !$OMP END PARALLEL DO

#endif

                ENDIF ! (if BND point)


              END DO ! DO JBND

            END DO ! DO IBND


            IF(improc.EQ.nmperr)WRITE(mdse,'(A)') &

                 'Finished MAPBDI 2nd loop.'


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


#ifdef W3_MPIBDI

            !KRL          Exchange (Note: for efficiency, post receives first)

            !KRL          MPI_ALLGATHERV would do this, but freezes for PGI and open_mpi

            !KRL          This suggests they use blocking SEND/RECV, so this is faster anyway and less implementation-dependent

            nrq = 0

            DO im = 1, nmproc

              IF ( im .NE. improc ) THEN

                nrq = nrq + 1

                tag = nmproc * im + improc

                CALL mpi_irecv ( mapbdi(1,nx_beg(im)), nx_size(im), mpi_real, im - 1, tag, mpi_comm_mwave, &

                     irq(nrq), ierr_mpi )

              END IF

            END DO

            DO im = 1, nmproc

              IF ( im .NE. improc ) THEN

                nrq = nrq + 1

                tag = nmproc * improc + im

                CALL mpi_isend( mapbdi(1,nx_beg(improc)), nx_size(improc), mpi_real, im - 1, tag, mpi_comm_mwave, &

                     irq(nrq), ierr_mpi )

              END IF

            END DO

            CALL mpi_waitall( nrq, irq, mpi_status_ignore, ierr_mpi )

            CALL mpi_barrier ( mpi_comm_mwave, ierr_mpi )

#endif


            IF(improc.EQ.nmperr)WRITE(mdse,'(A/)') &

                 '  Finished generating map with distances to boundary.'


            !...notes regarding old method of doing what we just did

            !!HT:

            !!HT: (1)

            !!HT:

            !!HT: CHANGE array is used to identify grid points that still need to

            !!HT: be processed, and that are adjacent to points that have been

            !!HT: processed. Only those points can be updated in this step of the

            !!HT: loop started above here. The two loops below set the CHANGE array.

            !!HT:

            !!HT: (2)

            !!HT:

            !!HT: CHANGD identify if more points have been updated

            !!HT:

            !!HT: STX and STY are partial normalized distances, defined as the

            !!HT: physical distance Delta Y ( FACTOR * SY ) and Delta X

            !!HT: ( FACTOR * SX * XLAT(JDST) ), devided by the sistance traveled,

            !!HT: which is CgMAX * DTMAX. CgMAX is approximately 1.15 * CgDEEP,

            !!HT: or 1.15 * 0.5 * C_DEEP = 0.58 * GRAV / SIG(1). Since SIG(1) and

            !!HT: DTMAX may vary, these two factors are not included in MAPBDI.

            !!HT:

            !!HT: This defines MAPBDI similar to an inverse CFL number.

            !!HT:

            !!HT: (3)

            !!HT:

            !!HT: ERROR : Should be CLAT(JDST), not CLATI(JDST) : "STX    = FACTOR * SX * CLATI(JDST) / ( 0.58 * GRAV )"


            ! 1.e Test output

            !

            !!HT: Note that SIG(1) and DTMAX are included here so that the map defines

            !!HT: how many time steps DTMAX it takes to reach this place.


#ifdef W3_T

            WRITE (mdst,9013)

            DO jdst=ny,1 , -1

              WRITE (mdst,9014) nint(mapbdi(jdst,:)*sig(1)/dtmax)

            END DO

#endif

            !

          END IF

        END DO

        flgbdi = .true.

      END IF

#ifdef W3_SCRIPNC

    END IF

#endif

    !

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

    ! 2.  Data sources for reconcilliation

    ! 2.a Loop over grids, processing check

    !


    !!HT: GRDHGH(GDST,0) was set in WMGLOW to identify how many grids may

    !!HT: contribute from higher ranks to the present grid (GDST).


    ALLOCATE ( i1(nrgrd,nmproc), i2(nrgrd,nmproc),                  &

         i3(nrgrd), i4(nrgrd), stat=istat )

    check_alloc_status( istat )


#ifdef W3_DIST

    ltag0  = 0

#endif


#ifdef W3_SCRIPNC

    !    If reading/writing SCRIP files, need to determine in advance which it is to avoid race condition:

    !       Processor writing file before other processor can check for it

    ngrdrank = shape(grdhgh)

    ALLOCATE( lgrdread(ngrdrank(1)-1, ngrdrank(2)), stat=istat )

    check_alloc_status( istat )

    ALLOCATE(lgrdwrite(ngrdrank(1)-1, ngrdrank(2)), stat=istat )

    check_alloc_status( istat )

    DO gdst=1, nrgrd

      DO jj = 1, grdhgh(gdst,0)

        IF ( grdhgh(gdst,0) .EQ. 0 ) THEN

          !             If no remap, then no file

          lgrdread(gdst,jj)  = .false.

          lgrdwrite(gdst,jj) = .false.

        ELSE

          gsrc   = grdhgh(gdst,jj)

          interp_file1 = "rmp_src_to_dst_conserv_xxx_xxx.nc"

          WRITE(cdst, "(I3.3)") gdst

          WRITE(csrc, "(I3.3)") gsrc

          interp_file1(24:26) = csrc

          interp_file1(28:30) = cdst

          INQUIRE(file=interp_file1, exist=l_read)

          !             At this point, file either exists already (L_READ = .TRUE.) or needs to be written

          lgrdread(gdst,jj)  = l_read

          lgrdwrite(gdst,jj) = .NOT. l_read

        END IF

      END DO

    END DO

#endif

#ifdef W3_MPI

    IF (lscripnc) CALL mpi_barrier(mpi_comm_mwave, ierr_mpi)

#endif


    lowrank_grid : DO gdst=1, nrgrd


#ifdef W3_T38

      CALL date_and_time ( cdate_time(1), cdate_time(2), cdate_time(3), date_time)

      beg_time(2) = ((date_time(5)*60 + date_time(6))*60 + date_time(7))*1000 + date_time(8)

      elapsed_time = beg_time(2) - beg_time(1)

      WRITE(nmyout,*) "WMGHGH, LOOP LOWRANK_GRID, GDST= ", gdst, " START: ", elapsed_time, " MSEC"

#endif


      ! Test output

#ifdef W3_T

      WRITE (mdst,9020) gdst, grdhgh(gdst,0)

#endif

#ifdef W3_SCRIP

      IF ( improc.EQ.nmperr.AND.t38 )WRITE(mdst,*)'GDST = ',gdst,' OUT OF ',nrgrd

#endif


      !

      IF ( grdhgh(gdst,0) .EQ. 0 ) THEN ! no grids of higher rank than this

        ! one.

#ifdef W3_T

        WRITE (mdst,9021)

#endif

      ELSE ! processing required


        !

        ! 2.b Process grid

        ! 2.b.1 Preparations

        !

        !!HT: Grid I has higher rank grids covering it, we now set up MAPTST

        !!HT:   MAPTST shows from which gid the data is averages.

        !!HT:   INFLND inferred land points based on land in high-res grids.

        !!HT:

        CALL w3seto ( gdst, mdse, mdst )

        CALL w3setg ( gdst, mdse, mdst )

        CALL wmsetm ( gdst, mdse, mdst )


        !  W3SETG set ICLOSE for us, and we have determined that there is

        !   interaction between high and low rank. So this is a good point

        !   to check the closure type.


        IF ( iclose .EQ. iclose_trpl ) THEN

          IF ( improc.EQ.nmperr ) &

               WRITE(mdse,*)'SUBROUTINE WMGHGH IS'// &

               ' NOT YET ADAPTED FOR TRIPOLE GRIDS. STOPPING NOW.'

          CALL extcde ( 1 )

        END IF


        ALLOCATE ( maptst(ny,nx), inflnd(ny,nx), stat=istat )

        check_alloc_status( istat )

        maptst = 0

        inflnd = 0


        !################################################################

        ! Start new block of code: Calculate weights by calling SCRIP interface

        !################################################################


        ! Notes on grid variables:

        ! GRIDS(GSRC)%{grid variable} (src grid, high rank, high resolution grid)

        ! GRIDS(GDST)%{grid variable} (dst grid, low rank, low resolution grid)


        !   At this point, we are working on a particular low rank (dst) grid.

        !   We will save our weight information in the structure "ALLWGTS".

        !   For this dst grid, it is possible to have many src grids. That is

        !      why we store it this way.

        !   First, we ALLOCATE ALLWGTS from 1 up to the largest value of all

        !      the possible source grids. This will be referenced as "GSRC"

        !      Not every value of GSRC will be filled (e.g. "1" usually isn't filled)

        !      but since we are doing this as a derived data type, we are still

        !      efficient in terms of memory usage.

        !   Inside SCRIP interface, we have:

        !      type weight_data

        !         integer (kind=int_kind)              :: n    ! number of weights for

        !                                               dst cell, formerly npnts(:)

        !         real    (kind=dbl_kind), allocatable :: w(:) ! weights, sized by n,

        !                                               formerly wxwy(:,:)

        !         integer (kind=int_kind), allocatable :: k(:) ! source grid cells,

        !                                               sized by n, formerly KSRC(:,:)

        !      end type weight_data

        !   ....

        !     type(weight_data), allocatable :: WGTDATA(:)

        !   ....

        !     ALLOCATE ( WGTDATA(grid2_size), STAT=ISTAT ) ! grid2=destination grid

        !     CHECK_ALLOC_STATUS ( ISTAT )


#ifdef W3_SCRIP

        njdst=ny

        nidst=nx

        ALLOCATE ( allwgts(maxval(grdhgh)), stat=istat )

        check_alloc_status( istat )

#endif


        ! Next, we loop through the src grids for the dst grid that we are working on.

#ifdef W3_SCRIP

        DO jj=1, grdhgh(gdst,0)

#endif

#ifdef W3_T38

          CALL date_and_time ( cdate_time(1), cdate_time(2), cdate_time(3), date_time)

          beg_time(3) = ((date_time(5)*60 + date_time(6))*60 + date_time(7))*1000 + date_time(8)

          elapsed_time = beg_time(3) - beg_time(1)

          WRITE(nmyout,*) "WMGHGH, LOOP JJ= ", jj, " START: ", elapsed_time, " MSEC"

#endif


#ifdef W3_SCRIP

          gsrc   = grdhgh(gdst,jj)

          nisrc=grids(gsrc)%NX

          njsrc=grids(gsrc)%NY ! only needed for diagnostics

#endif


          ! Next, we call SCRIP for this src grid.

          ! Conditions for calling SCRIP are:

          ! 1) Not using L_STOP: in this case, all processes need all the weight

          !    information, so all processes need to call SCRIP for all grid pairs

          ! OR

          ! 2) Using L_STOP, writing .nc files and not reading .nc files. With

          !    L_STOP, different processors are doing different things, and so

          !    have different settings for L_WRITE. L_READ is the same for all

          !    processors, since it is simply based on whether the file already

          !    exists.


#ifdef W3_SCRIPNC

          interp_file1 = "rmp_src_to_dst_conserv_xxx_xxx.nc"

          WRITE(cdst, "(I3.3)") gdst

          WRITE(csrc, "(I3.3)") gsrc

          interp_file1(24:26) = csrc

          interp_file1(28:30) = cdst

          l_read = lgrdread(gdst, jj)

#endif

#ifdef W3_T38

          CALL date_and_time ( cdate_time(1), cdate_time(2), cdate_time(3), date_time)

          beg_time(4) = ((date_time(5)*60 + date_time(6))*60 + date_time(7))*1000 + date_time(8)

          elapsed_time = beg_time(3) - beg_time(1)

          WRITE(nmyout,*) "WMGHGH, SCRIP WRAPPER START: ", elapsed_time, " MSEC"

#endif

#ifdef W3_SCRIPNC

          IF (l_stop) l_write = (improc .EQ. improc_assign)

#endif


#ifdef W3_SCRIPNC

          IF(l_stop.AND.l_read)THEN

            IF ( improc.EQ.nmperr ) &

                 WRITE(mdse,'(A)')'ERROR: You should either have SCRIP_STOP '// &

                 'or remapping (.nc) files. Not both. We will exit now.'

            CALL extcde (505)

          ENDIF

#endif


#ifdef W3_SCRIP

          called_scrip=.false. ! initialize

#endif


#ifdef W3_SCRIPNC

          IF ((.NOT. l_stop) .OR. ((.NOT. l_read) .AND. l_write)) THEN

#endif

#ifdef W3_SCRIP

            IF (l_stop) THEN ! we are sending different grids to different processors

              WRITE(mdse,'(A,2(I5),A,I5)')'Calling SCRIP for GSRC,GDST = ', &

                   gsrc,gdst,' on processor ',improc

            ELSEIF(improc.EQ.nmperr)THEN

              WRITE(mdse,'(A,2(I5))')'Calling SCRIP interface for GSRC,GDST = ', &

                   gsrc,gdst

            ENDIF

            CALL scrip_wrapper (gsrc, gdst,         &

                 grids(gsrc)%MAPSTA,grids(gsrc)%MAPST2,flagll,      &

                 grids(gsrc)%GRIDSHIFT,l_write,l_read,t38)

            called_scrip=.true.

#endif

#ifdef W3_SCRIPNC

          END IF

#endif

#ifdef W3_SCRIP

          CALL flush(mdse)

#endif

#ifdef W3_SCRIPNC

          IF (l_stop) THEN

            IF (.NOT. l_read) THEN

              improc_assign = improc_assign + 1

              IF (improc_assign .GT. nmproc) improc_assign = 1

              IF(called_scrip)THEN ! we called scrip_wrapper, so we need

                ! to deallocate before leaving

                dst_grid_size=nidst*njdst

                DO kdst=1,dst_grid_size

                  DEALLOCATE(wgtdata(kdst)%W, stat=istat )

                  check_dealloc_status( istat )

                  DEALLOCATE(wgtdata(kdst)%K, stat=istat )

                  check_dealloc_status( istat )

                END DO

                DEALLOCATE(wgtdata, stat=istat )

                check_dealloc_status( istat )

              END IF

              cycle ! cycle out of this loop :  DO JJ=1, GRDHGH(GDST,0)

            END IF

          END IF

#endif

#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(4)

          WRITE(nmyout,*) "WMGHGH, SCRIP WRAPPER, GSRC= ", gsrc, " TOOK ", elapsed_time, " MSEC"

#endif


#ifdef W3_SCRIP

          IF(.NOT.called_scrip)THEN ! we should not be here, since we need

            ! WGTDATA(KDST)%N which is created by SCRIP

            IF ( improc.EQ.nmperr ) WRITE(mdse,'(A)')'ERROR: we '// &

                 'should have cycled out by now. We will exit now.'

            CALL extcde (506)

          ENDIF

#endif


          ! SCRIP has now created the data strucure "WGTDATA" and stored the weights

          ! in it. However, this is only for the present src grid. We want to store the

          ! data for all the src grids. Thus, we use a new data structure of type

          ! "ALLWGT" to store this data. First though, we need to ALLOCATE it:

          ! (note: "k" is equivalent to isea, but includes *all* points)


#ifdef W3_SCRIP

          dst_grid_size=nidst*njdst

          ALLOCATE(allwgts(gsrc)%WGTDATA(dst_grid_size),stat=istat)

          check_alloc_status( istat )

          DO kdst=1,dst_grid_size

            ALLOCATE(allwgts(gsrc)%WGTDATA(kdst) &

                 %W(wgtdata(kdst)%N),stat=istat)

            check_alloc_status( istat )

            ALLOCATE(allwgts(gsrc)%WGTDATA(kdst) &

                 %K(wgtdata(kdst)%N),stat=istat)

            check_alloc_status( istat )

          END DO

#endif


          ! Now that we have it allocated, we can just copy WGTDATA into ALLWGTS


          ! Notes re: short and long way to do this:

          ! pgf90 on IBM Opteron, gfortran, g95, xlf, all tested ok with "short method"

          ! pgf90 on our linux workstations (Intel) requires the "long method"

          ! (possible compiler bug)

          ! ALLWGTS(GSRC)%WGTDATA = WGTDATA                               !short method


          ! BEGIN long method for filling derived data type "ALLWGTS"


#ifdef W3_SCRIP

          DO kdst=1,dst_grid_size

            allwgts(gsrc)%WGTDATA(kdst)%N=wgtdata(kdst)%N

            allwgts(gsrc)%WGTDATA(kdst)%NR0=wgtdata(kdst)%NR0

            allwgts(gsrc)%WGTDATA(kdst)%NR2=wgtdata(kdst)%NR2

            allwgts(gsrc)%WGTDATA(kdst)%NRL=wgtdata(kdst)%NRL

            DO ipnt=1,wgtdata(kdst)%N

              allwgts(gsrc)%WGTDATA(kdst)%W(ipnt) &

                   =wgtdata(kdst)%W(ipnt)

              allwgts(gsrc)%WGTDATA(kdst)%K(ipnt) &

                   =wgtdata(kdst)%K(ipnt)

            END DO

          END DO

#endif


          ! END long method for filling derived data type "ALLWGTS"


          ! We're done with WGTDATA, so we can DEALLOCATE it. This is important,

          ! since it will be allocated again the next time SCRIP is called.


#ifdef W3_SCRIP

          DO kdst=1,dst_grid_size

            DEALLOCATE(wgtdata(kdst)%W, stat=istat )

            check_dealloc_status( istat )

            DEALLOCATE(wgtdata(kdst)%K, stat=istat )

            check_dealloc_status( istat )

          END DO

          DEALLOCATE(wgtdata, stat=istat )

          check_dealloc_status( istat )

#endif


          ! Here's a "test output" block of code to demonstrate how the weights can

          ! be called up from ALLWGTS...and to verify that the data is stored properly.

          ! (again note that "k" is equivalent to isea, but includes *all* points)


#ifdef W3_SCRIP

          IF(t38)THEN

            WRITE(mdst,'(/2A)')'      XDST       YDST     ', &

                 '     XSRC          YSRC          WXWY'

            DO jdst=1,njdst

              DO idst=1,nidst

                kdst=(jdst-1)*nidst+idst

                xdst=real(grids(gdst)%XGRD(jdst,idst))

                ydst=real(grids(gdst)%YGRD(jdst,idst))

                DO ipnt=1,allwgts(gsrc)%WGTDATA(kdst)%N

                  ksrc=allwgts(gsrc)%WGTDATA(kdst)%K(ipnt)

                  jsrc=int((ksrc-1)/nisrc)+1

                  isrc=ksrc-(jsrc-1)*nisrc

                  xsrc=real(grids(gsrc)%XGRD(jsrc,isrc))

                  ysrc=real(grids(gsrc)%YGRD(jsrc,isrc))

                  wxwy=allwgts(gsrc)%WGTDATA(kdst)%W(ipnt)

                  WRITE(mdst,'(5(1X,F12.5))')xdst,ydst,xsrc, &

                       ysrc,wxwy

                END DO

              END DO

            END DO !  DO JDST=1,NJDST

          ENDIF ! IF(T38)THEN

#endif


#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(3)

          WRITE(nmyout,*) "WMGHGH, LOOP JJ, GSRC= ", gsrc, " TOOK ", elapsed_time, " MSEC"

#endif


#ifdef W3_SCRIP

        END DO ! DO JJ=1, GRDHGH(GDST,0)

        gsrc = -999 ! unset grid

#endif


        ! If SCRIPNC and L_STOP, then cycle LOWRANK_GRID loop and deallocate

        ! storage associated with dst grid.


#ifdef W3_SCRIPNC

        IF (l_stop) THEN

          IF ( ALLOCATED(maptst) ) THEN

            DEALLOCATE ( maptst, stat=istat )

            check_dealloc_status( istat )

          END IF

          IF ( ALLOCATED(inflnd) ) THEN

            DEALLOCATE ( inflnd, stat=istat )

            check_dealloc_status( istat )

          END IF

          IF ( ALLOCATED(allwgts) ) THEN

            DO jj=1, grdhgh(gdst,0)

              gsrc = grdhgh(gdst,jj)

              IF ( ASSOCIATED(allwgts(gsrc)%WGTDATA) ) THEN

                DO kdst=1, dst_grid_size


                  !#########################################################################################

                  !menta: for some reason gfortran complains that these lines are too long. Unindenting them

                  IF ( ALLOCATED(allwgts(gsrc)%WGTDATA(kdst)%W) ) THEN

                    DEALLOCATE ( allwgts(gsrc)%WGTDATA(kdst)%W, stat=istat )

                    check_dealloc_status( istat )

                  END IF

                  IF ( ALLOCATED(allwgts(gsrc)%WGTDATA(kdst)%K) ) THEN

                    DEALLOCATE ( allwgts(gsrc)%WGTDATA(kdst)%K, stat=istat )

                    check_dealloc_status( istat )

                  END IF

                  !#########################################################################################

                END DO

                DEALLOCATE ( allwgts(gsrc)%WGTDATA, stat=istat )

                check_dealloc_status( istat )

                NULLIFY ( allwgts(gsrc)%WGTDATA )

              END IF

            END DO

            DEALLOCATE ( allwgts, stat=istat )

            check_dealloc_status( istat )

          END IF

          cycle lowrank_grid

        END IF

#endif


        !################################################################

        ! End new block of code: Calculate weights by calling SCRIP interface

        !################################################################


        ! 2.b.2 Find points used for boundary data in higher ranked grids

        !

        !!HT: These points are marked in MAPTST as negative values to assure

        !!HT: that the grid poits used for boundary data are not getting

        !!HT: averaged values from high-reswolution grids as that will result

        !!HT: in cyclic, possibly non-conservative updating.

        !!HT:

        !!HT: NBI2S has all necessary data set in WMGLOW as called before WMGHGH.

        !!HT:

        !!HT: JJ loop goes over grids that reviously have been identified as

        !!HT: getting data from the grid presently cousidered.

        !

        ! notes: The purpose of this is to identify points

        !        that should not be used in the averaging procedure. It is

        !        related to statement in Tolman (OM, 2008):  "Second, Eq (7) is not

        !        applied to grid points in the low resolution grid that contribute

        !        to boundary data for the high resolution grid. This avoids cyclic

        !        updating of data between grids.


        ! notes:   GRDHGH(GDST,0) is number of grids of higher rank than the present

        !                         grid (GDST)

        !          GRDHGH(GDST,1...etc.) is the grid number


        ! notes: Setting MAPTST=negative here is probably overkill, since it means

        !        we will not have weights for this point. However, to change this,

        !        we would need a new variable to use in its place, since we need

        !        to mark the point for use in STMASK determination.


        DO jj=1, grdhgh(gdst,0)

          gsrc   = grdhgh(gdst,jj)

          DO ib=1, SIZE(mdatas(gsrc)%NBI2S(:,1))

            IF ( mdatas(gsrc)%NBI2S(ib,1) .EQ. gdst ) THEN

              idst     = mapsf(mdatas(gsrc)%NBI2S(ib,2),1)

              jdst     = mapsf(mdatas(gsrc)%NBI2S(ib,2),2)

              maptst(jdst,idst) = - gsrc

            END IF

          END DO

        END DO

        gsrc   = -999 ! unset grid

        !

        ! 2.b.3 Range of coverage per grid


        !!HT:

        !!HT: In this JJ loop, we go over all higher resolution grids to find

        !!HT: ranges that can be averaged to replace data in the 'I' (GDST) grid.!

        !!HT:


        ALLOCATE ( idstl(grdhgh(gdst,0)), idsth(grdhgh(gdst,0)),  &

             jdstl(grdhgh(gdst,0)), jdsth(grdhgh(gdst,0)),  &

             gridok(grdhgh(gdst,0)),bdist(grdhgh(gdst,0)),  &

             stat=istat                                     )

        check_alloc_status( istat )


        IF (old_method) THEN

          ALLOCATE (bdist_om(grdhgh(gdst,0)), stat=istat )

          check_alloc_status( istat )

        END IF


        !

        !       Notes: For case of lower ranked grid GDST being irregular, grid indices

        !              i and j do not correspond to x and y, so optimization

        !              by limiting search in manner of pre-curvilinear versions of

        !              WW3 is not appropriate.

        !

        IF ( (gtype .EQ. clgtype).or.(gtype .EQ. ungtype) ) THEN


          idstla   = 1

          idstha   = nx

          jdstla   = 1

          jdstha   = ny


        ELSE


          ! loop through higher ranked grids GSRC


          DO jj=1, grdhgh(gdst,0)

            gsrc   = grdhgh(gdst,jj)

            !!HT:

            !!HT: XL,XH, and YL,YH and the low and high (X,Y) values of the grid box

            !!HT: in the grid 'I' fro which the high-res data needs to be averaged.

            !!HT: To be efficient, we compute a range of high-res grid point that

            !!HT: could be considered, rather than looking through the whole grid.

            !!HT: This will work only for the old grids, not for the newer curvilinear

            !!HT: and unstructured grids.

            !!HT:

            !!HT: This sets the range in the low-res grid to consider.

            !

            !       Notes (HLT): outer edges already taken off here ...

            !              will not work in a simple way for spherical grids,

            !              so we don't even try ....

            !

            !       Notes: SX and SY are only used in cases where GTYPE .NE. CLGTYPE,

            !              i.e. regular grids. In case of regular grids, SX and SY

            !              can be replaced with HPFAC HQFAC, if desired.

            !

            ! find upper and lower bounds of higher ranks grids


            IF ( (grids(gsrc)%GTYPE .EQ. clgtype) .OR. &

                 (grids(gsrc)%GTYPE .EQ. ungtype) ) THEN


              !       Notes: in case of irregular grids, there is no obvious way to

              !              offset by dx/2 dy/2, so we omit that sliver (thus we increase

              !              search area slightly).


              xl=real(minval(grids(gsrc)%XGRD))

              yl=real(minval(grids(gsrc)%YGRD))

              xh=real(maxval(grids(gsrc)%XGRD))

              yh=real(maxval(grids(gsrc)%YGRD))


            ELSE


              xl     = grids(gsrc)%X0 + 0.5 * grids(gsrc)%SX

              xh     = grids(gsrc)%X0 + ( real(grids(gsrc)%NX) - 1.5 ) &

                   * grids(gsrc)%SX

              yl     = grids(gsrc)%Y0 + 0.5 * grids(gsrc)%SY

              yh     = grids(gsrc)%Y0 + ( real(grids(gsrc)%NY) - 1.5 ) &

                   * grids(gsrc)%SY


            ENDIF ! IF ( GRIDS(GSRC)%GTYPE .EQ. CLGTYPE )


            !

            ! find where this falls in the current (low) ranked grid


            IF ( flagll ) THEN

              idstl(jj) = 1

              idsth(jj) = nx

            ELSE


              ! Notes: from check "IF ( GTYPE .EQ. CLGTYPE ) THEN" above, we know that

              !         GTYPE .NE CLGTYPE....so it is safe to use SX SY etc here


              idstl(jj) = 2 + int( (xl-x0)/sx + 0.49 )

              idsth(jj) = 1 + int( (xh-x0)/sx - 0.49 )

            END IF


            jdstl(jj) = 2 + int( (yl-y0)/sy + 0.49 )

            jdsth(jj) = 1 + int( (yh-y0)/sy - 0.49 )


            idstl(jj) = max(  1 , idstl(jj) )

            idsth(jj) = min( nx , idsth(jj) )

            jdstl(jj) = max(  1 , jdstl(jj) )

            jdsth(jj) = min( ny , jdsth(jj) )

            !

#ifdef W3_T

            WRITE (mdst,9022) gsrc, idstl(jj),idsth(jj), &

                 jdstl(jj),jdsth(jj)

#endif

            !

          END DO  ! end loop through higher ranked grids

          gsrc   = -999 ! unset grid

          !


          ! save the extremities of that set of high-ranked grids

          idstla   = minval(idstl)

          idstha   = maxval(idsth)

          jdstla   = minval(jdstl)

          jdstha   = maxval(jdsth)


        ENDIF ! IF ( (GTYPE .EQ. CLGTYPE ) .or. (GTYPE .EQ. UNGTYPE))


        ! loop through higher ranked grids


        !

        ! 2.b.4 Point by point check

        !

        ! Notes: We loop through all grids of higher rank

        !        GSRC=the grid number of the higher rank grid.

        !        NLMAX is used for dimensioning purposes.

        !        It is apparently using the ratio between the resolution

        !        of the low rank grid (GDST) and high rank grid (GSRC)

        !        Obviously, we cannot use this calculation for irregular grids.


        nlmax  = 0

#ifdef W3_SCRIP

        nlmax_scrip=0

#endif

        DO jj=1, grdhgh(gdst,0)

          gsrc   = grdhgh(gdst,jj)


          ! Notes: NLMAX is used to dimension TMPINT,TMPRL, and to set ITAG and LTAG

          !        (MPI case).

          !        As we remove more of the older code, it may turn out that

          !        NLMAX is no longer needed, in which case we can remove this

          !        block of code. For example, the weights data structure is introduced

          !        to WW3 already dimensioned.


          IF ( grids(gdst)%GTYPE .EQ. clgtype ) THEN

            dx_max_gdst=maxval(grids(gdst)%HPFAC)

            dy_max_gdst=maxval(grids(gdst)%HQFAC)

          ELSEIF ( grids(gdst)%GTYPE .EQ. rlgtype ) THEN

            dx_max_gdst=grids(gdst)%SX

            dy_max_gdst=grids(gdst)%SY

          ELSE

            isfirst=1

            dist_max=0

            dist_min=0

            DO itri=1,grids(gdst)%NTRI

              DO it=1,3

                IF (it.eq.3) THEN

                  jt=1

                ELSE

                  jt=it+1

                END IF

                im1=grids(gdst)%TRIGP(it,itri)

                im2=grids(gdst)%TRIGP(jt,itri)

                edist=w3dist(flagll, real(grids(gdst)%XGRD(1,im1)), &

                     REAL(GRIDS(GDST)%YGRD(1,IM1)), &

                     REAL(GRIDS(GDST)%XGRD(1,IM2)), REAL(GRIDS(GDST)%YGRD(1,IM2)))

                IF (isfirst.eq.1) THEN

                  dist_max=edist

                  dist_min=edist

                  isfirst=0

                ELSE

                  IF (edist.gt.dist_max) THEN

                    dist_max=edist

                  END IF

                  IF (edist.lt.dist_min) THEN

                    dist_min=edist

                  END IF

                END IF

              END DO

            END DO

            dx_max_gdst=dist_max

            dy_max_gdst=dist_max

          END IF


          IF ( grids(gsrc)%GTYPE .EQ. clgtype ) THEN

            dx_min_gsrc=minval(grids(gsrc)%HPFAC)

            dy_min_gsrc=minval(grids(gsrc)%HQFAC)

          ELSEIF ( grids(gsrc)%GTYPE .EQ. rlgtype ) THEN

            dx_min_gsrc=grids(gsrc)%SX

            dy_min_gsrc=grids(gsrc)%SY

          ELSE

            isfirst=1

            dist_max=0

            dist_min=0

            DO itri=1,grids(gsrc)%NTRI

              DO it=1,3

                IF (it.eq.3) THEN

                  jt=1

                ELSE

                  jt=it+1

                END IF

                im1=grids(gsrc)%TRIGP(it,itri)

                im2=grids(gsrc)%TRIGP(jt,itri)

                edist=w3dist(flagll, real(grids(gsrc)%XGRD(1,im1)), &

                     REAL(GRIDS(GSRC)%YGRD(1,IM1)), &

                     REAL(GRIDS(GSRC)%XGRD(1,IM2)), REAL(GRIDS(GSRC)%YGRD(1,IM2)))

                IF (isfirst.eq.1) THEN

                  dist_max=edist

                  dist_min=edist

                  isfirst=0

                ELSE

                  IF (edist.gt.dist_max) THEN

                    dist_max=edist

                  END IF

                  IF (edist.lt.dist_min) THEN

                    dist_min=edist

                  END IF

                END IF

              END DO

            END DO

            dx_min_gsrc=dist_min

            dy_min_gsrc=dist_min

          END IF


          ! notes: original code was much simpler:

          !              NLMAX  = MAX ( NLMAX , (2+INT(SX/GRIDS(J)%SX+0.001)) *  &

          !                                     (2+INT(SY/GRIDS(J)%SY+0.001)) )


          nlmax = max( nlmax ,                             &

               (2+int(dx_max_gdst/dx_min_gsrc+0.001)) *  &

               (2+int(dy_max_gdst/dy_min_gsrc+0.001)) )


#ifdef W3_T38

          WRITE(mdst,*)'ratio 1 = ',(dx_max_gdst/dx_min_gsrc), &

               dx_max_gdst,dx_min_gsrc

          WRITE(mdst,*)'ratio 2 = ',(dy_max_gdst/dy_min_gsrc), &

               dy_max_gdst,dy_min_gsrc

          WRITE(mdse,*)'GSRC, NLMAX = ',gsrc,nlmax

#endif


#ifdef W3_SCRIP

          DO jdst=1, ny

            DO idst=1, nx

              kdst=(jdst-1)*nidst+idst

              nloc=allwgts(gsrc)%WGTDATA(kdst)%N

              nlmax_scrip=max(nlmax_scrip,nloc)

            END DO

          END DO

#endif


        END DO !  DO JJ=1, GRDHGH(GDST,0)

        gsrc=-999 ! unset grid


        ! Notes regarding 3 possible scenarios:

        !        If only using SCRIP, then

        !            * set NLMAX=NLMAX_SCRIP here.

        !            * TMPRL_OM will not be created

        !            * TMPRL will be calculated using SCRIP

        !        If only using old method

        !            * NLMAX is already set, and SCRIP switch does not exist, so

        !               nothing is done here

        !            * both TMPRL and TMPRL_OM will be dimensioned

        !            * TMPRL_OM will be calculated

        !            * TMPRL_OM will be copied to TMPRL for use

        !        If using both methods ("DO_CHECKING")

        !            * set NLMAX=MAX(NLMAX, NLMAX_SCRIP) here.

        !            * both TMPRL_OM and TMPRL will be created

        !            * both will be calculated using the 2 methods, and

        !               checked against each other

        !            * the SCRIP version of weights (TMPRL) will be the ones used.


#ifdef W3_SCRIP

        IF ( improc.EQ.nmperr.AND.t38 ) &

             WRITE(mdse,*) 'NLMAX,NLMAX_SCRIP=',nlmax,nlmax_scrip

        IF(do_checking)THEN

          nlmax = max(nlmax, nlmax_scrip)

        ELSE

          nlmax = nlmax_scrip

        ENDIF

        IF ( improc.EQ.nmperr.AND.t38 ) &

             WRITE(mdse,*) 'NEW NLMAX:',nlmax

#endif


        IF(nlmax.GT.100)THEN

          WRITE(mdse,'(/A,I8)') &

               'WARNING: unusually large value for NLMAX : ',nlmax

        END IF


        nrtot  = 0

        IF(old_method)THEN

          ALLOCATE ( tmpint_om(nx*ny,-4:nlmax), stat=istat )

          check_alloc_status( istat )

          ALLOCATE ( tmprl_om(nx*ny,0:nlmax),   stat=istat )

          check_alloc_status( istat )

        ENDIF

        ALLOCATE ( tmpint(nx*ny,-4:nlmax),       stat=istat )

        check_alloc_status( istat )

        ALLOCATE ( tmprl(nx*ny,0:nlmax),         stat=istat )

        check_alloc_status( istat )

        ALLOCATE ( tmplog(nx*ny),                stat=istat )

        check_alloc_status( istat )

        !

#ifdef W3_DIST

        ALLOCATE ( ltag(nlmax), stat=istat )

        check_alloc_status( istat )

        DO jj=1, nlmax

          ltag(jj) = jj + ltag0

        END DO

#endif

        !

        !!HT:

        !!HT: After the search range is set, we are actually searching in the

        !!HT: high-res grid. IDST, JDST are counters in the grid to which the

        !!HT: averaged data is to go. XA and YA are center locatons of target

        !!HT: grid. Necxt two loops over all relevant point in target grid.

        !!HT:


        ! Notes: This is the start of the large loop through the individual

        !        grid points of the low-rank grid.

        !        The checks below for JDST.LT.JDSTLA , IDST.LT.IDSTLA etc are to save

        !        time but will only be useful for the case of regular grids.


        lowrank_j : DO jdst=1, ny

          IF ( jdst.LT.jdstla .OR. jdst.GT.jdstha ) cycle


          lowrank_i : DO idst=1, nx

            IF ( idst.LT.idstla .OR. idst.GT.idstha ) cycle

            ! check that this is a sea point

            IF ( abs(mapsta(jdst,idst)) .NE. 1 ) cycle

            ! MAPTST: see Section 2.b.2 above

            IF ( maptst(jdst,idst) .LT. 0 ) cycle

            xa     = real(xgrd(jdst,idst)) ! old code: X0 + REAL(IDST-1)*SX

            ya     = real(ygrd(jdst,idst)) ! old code: Y0 + REAL(JDST-1)*SY


            !!HT: For each point in the target grid, loop over all relevant high-res

            !!HT: grid (JJ loop).


            nrok   = 0


            ! notes:   GRDHGH(GDST,0) is number of grids of higher rank than the present

            !          grid (GDST)

            !          GRDHGH(GDST,1...etc.) is the grid number


            DO jj=1, grdhgh(gdst,0)

              gsrc   = grdhgh(gdst,jj)


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

              ! Start counting using old method

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!


              !       Note for LLG: Assumption is made that the higher ranked grid

              !                     cannot be global.

              !

              !!HT: Set search range in [candidate] high-res grid.


              ! Notes: The quantities XL YL XH YH apear to be a bounding box in

              !        index space for later searching within the high rank grid (or

              !        otherwise making computations from the high rank grid). They

              !        are the distance from the coarse grid point to the origin of

              !        the high rank grid, measured in grid cells of the high rank

              !        grid. So, they are the i and j values in the high rank

              !        bounding the low rank grid cell.


              IF (old_method)THEN

                ! ...then  we do the counting using the old method


                ! Notes: Resulting "old method" variables are saved with "_OM" suffix.


                IF ( flagll ) THEN

                  dxc  = mod( 1080.+xa-grids(gsrc)%X0 , 360. )

                  xl   = 1. + (dxc-0.5*sx)/grids(gsrc)%SX

                  xh   = 1. + (dxc+0.5*sx)/grids(gsrc)%SX

                ELSE

                  xl   = 1. + (xa-grids(gsrc)%X0-0.5*sx)/grids(gsrc)%SX

                  xh   = 1. + (xa-grids(gsrc)%X0+0.5*sx)/grids(gsrc)%SX

                END IF

                yl     = 1. + (ya-grids(gsrc)%Y0-0.5*sy)/grids(gsrc)%SY

                yh     = 1. + (ya-grids(gsrc)%Y0+0.5*sy)/grids(gsrc)%SY


                isrcl  = nint(xl+0.01)

                isrch  = nint(xh-0.01)

                jsrcl  = nint(yl+0.01)

                jsrch  = nint(yh-0.01)


                IF ( isrcl.LT.1 .OR. isrch.GT.grids(gsrc)%NX .OR.          &

                     jsrcl.LT.1 .OR. jsrch.GT.grids(gsrc)%NY ) THEN

                  !                            dst point not in src grid, so go to next src grid

                  gridok(jj) = .false.  ! does this get used anywhere?

                  cycle ! leave GSRC loop

                END IF


                !!HT: Loop over search range in high-res grid, ISRC and JSRC loops.

                !!HT:  NR0_OM counts high-res grid points with MAPSTA=0, etc.

                !!HT:  NRL_OM separately identifies explitcit land points.

                !!HT:  BDIST_OM saves the boundary data from the source grid.

                !!HT:


                ! Notes: We appear to be searching for the smallest boundary distance and

                !        doing some counting

                !        Purpose of counting is unknown (for dimensioning?)


                ! Initialize

                nr0_om    = 0 ! counter of MAPSTA=0 (indicates

                ! excluded point)

                nrl_om    = 0 ! counter of MAPSTA=0 (indicates

                ! excluded point) and MAPST2=0

                ! (indicates land)

                nr1_om    = 0 ! counter of |MAPSTA|=1

                ! (indicates sea point)

                nr2_om    = 0 ! counter of |MAPSTA|=2

                ! (indicates boundary point)

                bdist_om(jj) = 9.99e33


                DO isrc=isrcl, isrch

                  DO jsrc=jsrcl, jsrch

                    IF (grids(gsrc)%MAPSTA(jsrc,isrc).EQ.0) THEN

                      ! excluded point

                      nr0_om    = nr0_om + 1


                      ! Notes: Q: What does MAPST2=0 indicate?

                      !        A: MAPST2 is the "second grid status map"

                      !           For disabled points (MAPSTA=0) , MAPST2 indicates land (0) or

                      !           excluded (1). For sea and active boundary points, MAPST2 indicates

                      !           a) ice coverage b) drying out of points c) land in moving grid or

                      !           inferred land in nesting and d) masked in two-way nesting


                      IF (grids(gsrc)%MAPST2(jsrc,isrc).EQ.0) &

                           nrl_om = nrl_om + 1

                    ELSE IF (abs(grids(gsrc)%MAPSTA(jsrc,isrc)) &

                         .EQ.1) THEN ! sea point

                      nr1_om    = nr1_om + 1


                      ! Notes: check against stored "distance to boundary point"

                      !        This BDIST_OM array will be used later, when we select

                      !        the high rank grid to average from.


                      bdist_om(jj) = min( bdist_om(jj) ,     &

                           mdatas(gsrc)%MAPBDI(jsrc,isrc) )

                    ELSE IF (abs(grids(gsrc)%MAPSTA(jsrc,isrc)) &

                         .EQ.2) THEN ! bnd point

                      nr2_om    = nr2_om + 1

                    END IF

                  END DO ! DO JSRC=JSRCL, JSRCH

                END DO ! DO ISRC=ISRCL, ISRCH


              END IF ! (if OLD_METHOD)


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

              ! Done with counting using old method.

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

              ! Start counting using new method

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!


              ! Initialize


#ifdef W3_SCRIP

              bdist(jj) = 9.99e33

#endif


              ! Notes on variables used here:

              ! IDST, JDST given by loop, NIDST set above, the rest we need to set here


#ifdef W3_SCRIP

              nisrc=grids(gsrc)%NX

              kdst=(jdst-1)*nidst+idst

#endif


#ifdef W3_SCRIP

              DO ipnt=1,allwgts(gsrc)%WGTDATA(kdst)%N

                ksrc=allwgts(gsrc)%WGTDATA(kdst)%K(ipnt)

                jsrc=int((ksrc-1)/nisrc)+1

                isrc=ksrc-(jsrc-1)*nisrc

                IF (abs(grids(gsrc)%MAPSTA(jsrc,isrc)).EQ.1) THEN

#endif

                  ! sea point

#ifdef W3_SCRIP

                  bdist(jj) = min( bdist(jj) ,               &

                       mdatas(gsrc)%MAPBDI(jsrc,isrc) )

                ELSE

                  IF ( improc.EQ.nmperr ) &

                       WRITE(mdse,*) &

                       'we masked non-sea points. (coding error)'

                  CALL extcde ( 999 )

                END IF

              END DO

#endif


              ! Pull NR0, etc. from storage...

#ifdef W3_SCRIP

              nr0 = allwgts(gsrc)%WGTDATA(kdst)%NR0

#endif

              !                             counter of MAPSTA=0 (indicates excluded point)

#ifdef W3_SCRIP

              nrl = allwgts(gsrc)%WGTDATA(kdst)%NRL

#endif

              !                             counter of MAPSTA=0 (indicates excluded point)

              !                             and MAPST2=0 (indicates land)

#ifdef W3_SCRIP

              nr1 = allwgts(gsrc)%WGTDATA(kdst)%N

#endif

              !                             counter of |MAPSTA|=1 (indicates sea point)

#ifdef W3_SCRIP

              nr2 = allwgts(gsrc)%WGTDATA(kdst)%NR2

#endif

              !                             counter of |MAPSTA|=2 (indicates boundary point)


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

              ! Finished counting using new method.

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!


              ! Compare results

              IF(do_checking)THEN

                ! then it is OK to compare with the values that we got using the old method

#ifdef W3_T38

                WRITE(mdst,*)'STARTING TEST 1'

#endif

                IF(nr0_om.NE.nr0)THEN

                  IF ( improc.EQ.nmperr )WRITE (mdse,'(/1A,2(I8))') &

                       ' *** ERROR WMGHGH: NR0_OM,NR0 = ',nr0_om,nr0

                  CALL extcde ( 999 )

                ENDIF

                IF(nr1_om.NE.nr1)THEN

                  IF ( improc.EQ.nmperr )WRITE (mdse,'(/1A,2(I8))') &

                       ' *** ERROR WMGHGH: NR1_OM,NR1 = ',nr1_om,nr1

                  CALL extcde ( 999 )

                ENDIF

                IF(nr2_om.NE.nr2)THEN

                  IF ( improc.EQ.nmperr )WRITE (mdse,'(/1A,2(I8))') &

                       ' *** ERROR WMGHGH: NR2_OM,NR2 = ',nr2_om,nr2

                  CALL extcde ( 999 )

                ENDIF

                IF(nrl_om.NE.nrl)THEN

                  IF ( improc.EQ.nmperr )WRITE (mdse,'(/1A,2(I8))') &

                       ' *** ERROR WMGHGH: NRL_OM,NRL = ',nrl_om,nrl

                  CALL extcde ( 999 )

                ENDIF

                IF(bdist_om(jj).NE.bdist(jj))THEN

                  IF ( improc.EQ.nmperr ) &

                       WRITE (mdse,'(/2A,2(F12.5))') &

                       ' *** ERROR WMGHGH: ', &

                       ' BDIST_OM(JJ),BDIST(JJ) = ', &

                       bdist_om(jj),bdist(jj)

                  CALL extcde ( 999 )

                ENDIF

#ifdef W3_T38

                WRITE(mdst,*)'PASSED TEST 1'

#endif

              END IF ! (if DO_CHECKING)


              ! Notes: We are done with the counting. If we didn't use SCRIP to get NR0,

              !  etc., then we need to set them using the _OM variables.


              IF(.NOT.lscrip)THEN

                nr0=nr0_om

                nr1=nr1_om

                nr2=nr2_om

                nrl=nrl_om

                bdist=bdist_om

              END IF


              ! Notes: Potential future improvement: for irregular grids, it would make

              ! more sense to use the overlapped area, rather than simply counting cells

              ! to decide on "inferred land". However, since grid cell size it typically

              ! fairly uniform locally, the current approach will suffice for now.


              ! Notes: This is the only place that the "NRL" "NR0" "NR1" and "NR2" variables

              !        are used directly. They affect MAPST2 indirectly below.

              !        The calculation itself is essentially "50% or more of the grid

              !        cells are land".


              IF ( nrl .GT. (nr0+nr1+nr2)/2 ) THEN


                ! Notes: This is not considered an OK grid (NROK is not incremented) and

                !        it is considered "inferred land"


                inflnd(jdst,idst) = 1

              ELSE

                gridok(jj) = nr1.GT.0 .AND. nr2.EQ.0


                ! Notes: for a grid cell to be considered "OK", we require that there is

                ! at least one sea point being used, and no boundary points being used


                IF ( gridok(jj) ) nrok = nrok + 1

              END IF


            END DO !  GSRC loop

            gsrc=-999 ! unset grid


            IF ( nrok .EQ. 0 ) THEN


              ! Notes:  exit IDST loop since there are no "OK" source grid cells for this

              ! dst point. At this point, INFLND could be 1, but isn't necessarily 1


              cycle


            ELSE


              ! Notes: If any grids are OK for this dst point, then we override any prior

              ! setting of INFLD=1. Apparently this is for the situation of having some src

              ! grids giving INFLD=1 and another giving INFLD=0 for the same dst point.

              ! I wouldn't expect this to happen very often.


              inflnd(jdst,idst) = 0


            END IF

            !

            ! 2.b.5 Select source grid

            !

            ! Notes: It appears that we are selecting the high rank grid from

            !        which we will perform the averaging.

            !        The code is written such that the first higher rank

            !        grid that we find has the rank that we want, but isn't necessarily the

            !        grid that we want.

            !        Are grids necessarily in order of rank? If so, then we want the grid

            !        that is higher rank but of nearest rank to the present grid.

            !        Anyway, once we have decided on the grid rank that we want, we select

            !        the specific grid according to criterion: larger distance to

            !        boundary = better

            !        Keep in mind that this grid is selected for *this* (IDST,JDST) and not

            !        necesssarily for the next...


            jf     = 0


            !!HT: Another loop over all high-res grid to decide which grid will

            !!HT: be used to average data. If more than 1 grids are available,

            !!HT: the boundary distance in the high-res grid, stored in BDIST is

            !!HT: used to make the choice.


            !!ER: Old logic was to select a grid that is of the "next rank up,

            !!ER: for which data is available". This was done by searching

            !!ER: from 1 to GRDHGH (remember that the available source grids

            !!ER: are in order of rank), and exiting when the rank increased.

            !!ER: The problem with selecting the "lowest rank grid with rank

            !!ER: greater than that of GDST" is that at this point, we have

            !!ER: no knowledge of what will be masked in that SRC grid, since

            !!ER: we haven't updated MAPSTA for that GSRC yet, based on what

            !!ER: points are covered by higher rank grids (in case of masking

            !!ER: computations). We can avoid this problem by reversing the

            !!ER: order of GSRC search (from highest rank to lowest rank of

            !!ER: higher rank). For example, grid 1 wants data from grid 2,

            !!ER: but just gets zeros because grid 2 is masked there, because

            !!ER: grid 2 is masked by grid 3 in Section 2.3.2 below. Going

            !!ER: directly to GSRC=3 for GDST=1 (skipping GSRC=2) avoids

            !!ER: this: the highest rank at that location will never be

            !!ER: masked by a higher rank grid.


            DO jj=grdhgh(gdst,0),1,-1

              !old              DO JJ=1, GRDHGH(GDST,0) ! used before Aug 2014


              gsrc   = grdhgh(gdst,jj)


              IF ( gridok(jj) ) THEN

                IF ( jf .EQ. 0 ) THEN ! we haven't already found a grid

                  jf     = gsrc ! now we have found a grid.

                  jr     = grank(gsrc)

                  ! this is the rank that we want....the rank of the first grid that we find

                  jd     = bdist(jj) ! larger distance = better

                ELSE

                  ! we already found a grid, but maybe this one is better

                  IF ( grank(gsrc) .NE. jr ) EXIT

                  ! this is not the rank that we want

                  IF ( bdist(jj) .GT. jd ) THEN

                    ! we like this grid better

                    jf     = gsrc

                    jd     = bdist(jj)

                  END IF

                END IF

              END IF

            END DO

            gsrc=jf

#ifdef W3_T38

            WRITE(mdst,'(A,2(I8),A,I8)')'For grid point IDST,JDST = ',idst,jdst,', we selected GSRC = ',gsrc

#endif


            !!HT: Data for grid point IDST,JDST in the low-res grid will be taken from

            !!HT: high-res grid GSRC.


            maptst(jdst,idst) = gsrc

            !

            ! 2.b.6 Store data (temp)

            !

            ! Notes: This section is for calculations of weights for the

            !        area-weighted averaging.


            nrtot  = nrtot + 1

            tmpint(nrtot,-4) = idst

            tmpint(nrtot,-3) = jdst

            tmpint(nrtot,-2) = mapfs(jdst,idst)

            tmpint(nrtot,-1) = gsrc

            tmprl(nrtot, 0) = jd * sig(1) / dtmax


            ! Notes: Calculation for XL YL XH YH is same as it was in section 2.b.4, so

            !        see notes in that section.


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

            !...Begin block of code for computing weights using old method

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!


            IF (old_method)THEN

              ! it is OK to do the counting using the old method

              ! (These variables are saved with "_OM" suffix)


              DO itmp=-4,-1

                tmpint_om(nrtot,itmp)=tmpint(nrtot,itmp)

              END DO

              tmprl_om(nrtot,0)=tmprl(nrtot,0)


              IF ( flagll ) THEN

                dxc    = mod( 1080.+xa-grids(gsrc)%X0 , 360. )

                xl     = 1. + (dxc-0.5*sx)/grids(gsrc)%SX

                xh     = 1. + (dxc+0.5*sx)/grids(gsrc)%SX

              ELSE

                xl     = 1. + (xa-grids(gsrc)%X0-0.5*sx)/grids(gsrc)%SX

                xh     = 1. + (xa-grids(gsrc)%X0+0.5*sx)/grids(gsrc)%SX

              END IF

              yl     = 1. + (ya-grids(gsrc)%Y0-0.5*sy)/grids(gsrc)%SY

              yh     = 1. + (ya-grids(gsrc)%Y0+0.5*sy)/grids(gsrc)%SY


              ! Notes: Here, we save the search bounds. These bounds are given in terms of

              !        index space of the high rank grid. "L" and "H" here do *not* refer

              !        to grid rank! They refer to lower and upper bounds.


              isrcl    = nint(xl+0.01)

              isrch    = nint(xh-0.01)

              jsrcl    = nint(yl+0.01)

              jsrch    = nint(yh-0.01)


              wtot   = 0.

              nloc_om   = 0

              DO isrc=isrcl, isrch

                wx     = min(xh,real(isrc)+0.5) - max(xl,real(isrc)-0.5)

                DO jsrc=jsrcl, jsrch

                  IF (abs(grids(gsrc)%MAPSTA(jsrc,isrc)).EQ.1) THEN

                    ! sea point

                    wy     = min(yh,real(jsrc)+0.5) -               &

                         max(yl,real(jsrc)-0.5)

                    wtot   = wtot + wx*wy

                    nloc_om   = nloc_om + 1

                    ! Notes: check here that we are sufficiently dimensioned.

                    IF ( nloc_om .GT. nlmax ) THEN

                      IF ( improc.EQ.nmperr ) WRITE (mdse,1020)

                      CALL extcde(1020)

                    END IF

                    tmpint_om(nrtot,nloc_om) =  &

                         grids(gsrc)%MAPFS(jsrc,isrc)

                    tmprl_om(nrtot,nloc_om) = wx*wy

                  END IF

                END DO

              END DO

              tmpint_om(nrtot,0) = nloc_om

              tmprl_om(nrtot,1:nloc_om) = tmprl_om(nrtot,1:nloc_om) &

                   / wtot


            END IF ! (if OLD_METHOD)


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

            !...End block of code for computing weights using old method

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

            !...Begin block of code for "computing" weights using new method

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!


            ! Notes: Weights have already been computed by SCRIP.

            !        We just need to transfer them to TMPINT and TMPRL


#ifdef W3_SCRIP

            kdst=(jdst-1)*nidst+idst

            nloc=allwgts(gsrc)%WGTDATA(kdst)%N

            tmpint(nrtot,0) = nloc

            nisrc=grids(gsrc)%NX

#endif


            ! Test output

#ifdef W3_SCRIP

            IF ( improc.EQ.nmperr.AND.t38 ) THEN

              WRITE(mdst,*)'GSRC,KDST,NLOC = ',gsrc,kdst,nloc

            ENDIF

#endif


            ! Notes: check here that we are sufficiently dimensioned.


#ifdef W3_SCRIP

            IF ( nloc .GT. nlmax ) THEN

              IF ( improc.EQ.nmperr ) THEN

                WRITE (mdse,'(/2A,4(1x,I8))') &

                     ' *** ERROR WMGHGH: ', &

                     ' IDST,JDST,NLOC,NLMAX = ', &

                     idst,jdst,nloc,nlmax

                WRITE(mdse,1021)

              ENDIF

              CALL extcde(1021)

            END IF

            DO ipnt=1,nloc

              ksrc=allwgts(gsrc)%WGTDATA(kdst)%K(ipnt)

              jsrc=int((ksrc-1)/nisrc)+1

              isrc=ksrc-(jsrc-1)*nisrc

              tmpint(nrtot,ipnt) = grids(gsrc)%MAPFS(jsrc,isrc)

              tmprl(nrtot,ipnt)= &

                   allwgts(gsrc)%WGTDATA(kdst)%W(ipnt) ! WX*WY / WTOT

            END DO !  DO IPNT=1,NLOC

#endif


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

            !...End block of code for "computing" weights using new method

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

            !...Begin block of code that is just for testing

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!


            IF (do_checking)THEN

              ! compare with the values that we got using the old method

#ifdef W3_T38

              WRITE(mdst,*)'STARTING TEST 2'

#endif

              if (nloc.NE.nloc_om) THEN

                IF ( improc.EQ.nmperr )WRITE (mdse,'(/1A,2(I8))') &

                     ' *** ERROR WMGHGH: NLOC,NLOC_OM  = ',nloc,nloc_om

                CALL extcde ( 999 )

              END IF

              istop=0

              icount=0

              DO ipnt=1,nloc

                DO ipnt2=1,nloc

                  IF (tmpint_om(nrtot,ipnt).EQ.tmpint(nrtot,ipnt2))THEN

                    ! we found our point

                    icount=icount+1

                    IF(abs(tmprl_om(nrtot,ipnt)-tmprl(nrtot,ipnt2)) &

                         .GT.4.0e-5)then

                      IF ( improc.EQ.nmperr )WRITE &

                           (mdse,'(/2A,2(F12.5))') &

                           ' *** ERROR WMGHGH: ', &

                           ' *** TMPRL_OM(NRTOT,IPNT),TMPRL(NRTOT,IPNT2) = ', &

                           tmprl_om(nrtot,ipnt),tmprl(nrtot,ipnt2)

                      istop=1

                    END IF

                  END IF

                END DO

              END DO

              IF(icount.NE.nloc)THEN

                IF ( improc.EQ.nmperr )WRITE (mdse,'(/1A,2(I8))') &

                     ' *** ERROR WMGHGH: ICOUNT,NLOC = ',icount,nloc

                istop=1

              END IF

              IF(istop.EQ.1)THEN

                CALL extcde ( 999 )

              END IF

#ifdef W3_T38

              WRITE(mdst,*)'PASSED TEST 2'

#endif


            END IF ! (if both grids are regular grids)


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

            !...End block of code that is just for testing

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!


            nrok   = 0


          END DO lowrank_i !  DO IDST=1, NX

        END DO lowrank_j ! DO JDST=1, NY


#ifdef W3_T38

        WRITE(mdst,*)'WMGHGH Section 2.b.6 completed.'

#endif


        ! Notes: We are done with the counting. If we didn't use SCRIP to get

        !     TMPINT, TMPRL, then we need to set them using the _OM variables.


        IF(.NOT.lscrip)THEN

          tmpint=tmpint_om

          tmprl=tmprl_om

        END IF


#ifdef W3_T

        WRITE (mdst,9023) gdst, nrtot

#endif

        !

        ! 2.c Set up masks based on stencil width of scheme and inferred land

        ! 2.c.1 Inferred land

        !

        !!HT: Inferred land from INFLND is added to MAPSTA / MAPST2

        !!HT:

        mapst2 = mapst2 - 4*mod(mapst2/4,2)

        mapst2 = mapst2 + 4*inflnd

        DO idst=1, nx

          DO jdst=1, ny

            IF ( mapst2(jdst,idst).GT.0 ) mapsta(jdst,idst) =   &

                 - abs(mapsta(jdst,idst))

          END DO

        END DO

        !

        ! 2.c.2 Masking

        !

        !!HT: This is masking in the low-res grid to identify where the grid

        !!HT: is covered by high-res grids, and far enough away from the

        !!HT: high-res grid edges so that no dynamic computations are needed

        !!HT: in the low-res grid.

        !!HT:

        ALLOCATE ( stmask(ny,0:nx+1), maski(ny,nx), stat=istat )

        check_alloc_status( istat )

        IF ( mdatas(gdst)%MSKINI ) THEN

          DEALLOCATE ( mdatas(gdst)%MAPMSK, stat=istat )

          check_dealloc_status( istat )

        END IF

        ALLOCATE ( mdatas(gdst)%MAPMSK(ny,nx), stat=istat )

        check_alloc_status( istat )

        mapmsk => mdatas(gdst)%MAPMSK

        mdatas(gdst)%MSKINI = .true.


        mapmsk = mod(mapst2/8,2)

        mapst2 = mapst2 - 8*mapmsk


        !!HT: STMASK (logical) is used to start this up. We first use the point

        !!HT:        MAPTST that have been marked as used for boundary points in

        !!HT:        the corrsponding high-res grids.

        !!HT: NIT sets the stencil width of the propagation scheme, used to see

        !!HT:        how far we need to move in from the boundary points of

        !!HT:        the high-res grid to reach the area in the low-res grid

        !!HT:        where we do not need to compute.


        stmask(:,1:nx) = maptst .LT. 0

        stmask(:,0) = stmask(:,nx)

        stmask(:,nx+1) = stmask(:,1)


#ifdef W3_PR0

        nit    = 0

#endif

#ifdef W3_PR1

        nit    = ( 1 + int(dtmax/dtcfl-0.001) ) * 1

#endif

#ifdef W3_UQ

        nit    = ( 1 + int(dtmax/dtcfl-0.001) ) * 3

#endif

#ifdef W3_UNO

        nit    = ( 1 + int(dtmax/dtcfl-0.001) ) * 3

#endif


        idstla=2

        idstha=nx-1


        ! notes....bug fix: in official release 3.14, the if-then below

        !  was missing. This would produce incorrect results for a global grid that

        !  had a higher rank grid on the branch cut (180 to -180 or 360 to 0). See

        !  treatment of STMASK after the if-then statement. There, it is clear that

        !  it was intended that MASKI be available for i=1 and i=nx, ... but it wasn't

        !  available. Symptoms of bug: when using "T T" for masking options, a strip

        !  of land would be placed along the i-column just east of the branch cut.

        !  This would be seen in the global (low rank) grid.


        IF ( iclose.NE.iclose_none ) THEN

          idstla=1

          idstha=nx

        END IF


        DO jtmp=1, nit

          maski  = .false.

          DO idst=idstla,idstha

            DO jdst=2, ny-1

              IF ( .NOT. stmask(jdst,idst) .AND. (                     &

                   stmask(jdst+1,idst+1) .OR.  stmask(jdst+1,idst  ) .OR. &

                   stmask(jdst+1,idst-1) .OR.  stmask(jdst  ,idst-1) .OR. &

                   stmask(jdst-1,idst-1) .OR.  stmask(jdst-1,idst  ) .OR. &

                   stmask(jdst-1,idst+1) .OR.  stmask(jdst  ,idst+1) ) )  &

                   maski(jdst,idst) = .true.

            END DO

          END DO

          stmask(:,1:nx) = stmask(:,1:nx) .OR. maski

          stmask(:,0) = stmask(:,nx)

          stmask(:,nx+1) = stmask(:,1)

        END DO


        !!HT: Loop over all point from which low-res grid gets data from

        !!HT: high-res grid(s). Comparing to STMASK shows which points can be

        !!HT: masked out for computation.

        !!HT:

        !!HT: MAPMSK is stored in WMMDATMD for use in wave model.


        DO iloc=1, nrtot

          idst     = tmpint(iloc,-4)

          jdst     = tmpint(iloc,-3)

          tmplog(iloc) = stmask(jdst,idst)

          IF ( .NOT. stmask(jdst,idst) ) THEN

            mapmsk(jdst,idst) = 1

            IF ( flghg1 ) mapsta(jdst,idst) = -abs(mapsta(jdst,idst))

            maptst(jdst,idst) = 99

          END IF

        END DO


        IF ( flghg1 ) mapst2 = mapst2 + 8*mapmsk


        DEALLOCATE ( stmask, maski, stat=istat )

        check_dealloc_status( istat )


        !!HT: Now that all temporary data is stored, and all mosks are set, all

        !!HT: can be put from temporaty storage in permanent storage.

        !!HT:

        !!HT: Should require no modifications for newer grids ...

        !!HT:    ... unless more data is needed than for old grids .....


        !

        ! 2.d Set up mapping for staging data

        ! 2.d.1 Set counters / required array sizes

        !

#ifdef W3_SHRD

        isproc = 1

        ispro2 = 1

#endif

        i1     = 0

        i2     = 0

        i3     = 0

        i4     = 0


        DO iloc=1, nrtot


          jj     = tmpint(iloc,-1)

          hgstge(gdst,jj)%NTOT = hgstge(gdst,jj)%NTOT + 1

          isea   = tmpint(iloc,-2)

          CALL init_get_jsea_isproc(isea, jsea, isproc)

#ifdef W3_DIST

          isproc = isproc + croot - 1

#endif

          !

          i1(jj,isproc) = i1(jj,isproc) + 1

          IF ( tmplog(iloc) ) i2(jj,isproc) = i2(jj,isproc) + 1

          IF ( improc .EQ. isproc ) THEN

            hgstge(gdst,jj)%NSMX = max(hgstge(gdst,jj)%NSMX,tmpint(iloc,0))

          END IF


          DO jr=1, tmpint(iloc,0)

            isea   = tmpint(iloc,jr)

            CALL init_get_jsea_isproc_glob(isea, jj, jsea, ispro2)

            IF ( ispro2 .EQ. improc ) THEN

              hgstge(gdst,jj)%NSND = hgstge(gdst,jj)%NSND + 1

              IF ( tmplog(iloc) ) hgstge(gdst,jj)%NSN1 =              &

                   hgstge(gdst,jj)%NSN1 + 1

            END IF

          END DO

          !

        END DO


        hgstge(gdst,:)%NREC = i1(:,improc)

        hgstge(gdst,:)%NRC1 = i2(:,improc)

        !

        ! 2.d.2 ALLOCATE (DEALLOCATE in section 0 as needed)

        !

        DO gsrc=1, nrgrd

          IF ( hgstge(gdst,gsrc)%NREC .GT. 0 ) THEN

            ALLOCATE (                                          &

                 hgstge(gdst,gsrc)%LJSEA (hgstge(gdst,gsrc)%NREC), &

                 hgstge(gdst,gsrc)%NRAVG (hgstge(gdst,gsrc)%NREC), &

                 hgstge(gdst,gsrc)%IMPSRC(hgstge(gdst,gsrc)%NREC,  &

                 hgstge(gdst,gsrc)%NSMX),                          &

                 hgstge(gdst,gsrc)%ITAG  (hgstge(gdst,gsrc)%NREC,  &

                 hgstge(gdst,gsrc)%NSMX),                          &

                 hgstge(gdst,gsrc)%WGTH  (hgstge(gdst,gsrc)%NREC,  &

                 hgstge(gdst,gsrc)%NSMX),                          &

                 hgstge(gdst,gsrc)%SHGH  (sgrds(gsrc)%NSPEC,       &

                 hgstge(gdst,gsrc)%NSMX,                           &

                 hgstge(gdst,gsrc)%NREC), stat=istat               )

            check_alloc_status( istat )

#ifdef W3_T3

            hgstge(gdst,gsrc)%LJSEA  = -1

            hgstge(gdst,gsrc)%NRAVG  = -1

            hgstge(gdst,gsrc)%IMPSRC = -1

            hgstge(gdst,gsrc)%ITAG   = -1

            hgstge(gdst,gsrc)%WGTH   = -1.

#endif

          END IF

          IF ( hgstge(gdst,gsrc)%NSND .GT. 0 ) THEN

            ALLOCATE ( hgstge(gdst,gsrc)%ISEND (hgstge(gdst,gsrc)%NSND,5), &

                 stat=istat )

            check_alloc_status( istat )

#ifdef W3_T4

            hgstge(gdst,gsrc)%ISEND = -1

#endif

          END IF

          hgstge(gdst,gsrc)%INIT = .true.

        END DO

        !

        ! 2.d.3 Fill allocated arrays

        !

        flgrec = .true.

        i2     = i1 + 1

        i1     = 0

        i4     = hgstge(gdst,:)%NSND + 1

        i3     = 0


        DO iloc=1, nrtot


          isea   = tmpint(iloc,-2)

          jj     = tmpint(iloc,-1)

          nr0    = tmpint(iloc, 0)

          CALL init_get_jsea_isproc(isea, jsea, isproc)

#ifdef W3_DIST

          isproc = isproc + croot - 1

          flgrec = isproc .EQ. improc

#endif

          !

          IF ( tmplog(iloc) ) THEN

            i1(jj,isproc) = i1(jj,isproc) + 1

            irec          = i1(jj,isproc)

          ELSE

            i2(jj,isproc) = i2(jj,isproc) - 1

            irec          = i2(jj,isproc)

          END IF


          IF ( flgrec ) THEN

            hgstge(gdst,jj)%LJSEA(irec) = jsea

            hgstge(gdst,jj)%NRAVG(irec) = nr0

            hgstge(gdst,jj)%WGTH(irec,:nr0) = tmprl(iloc,1:nr0)

#ifdef W3_DIST

            hgstge(gdst,jj)%ITAG(irec,:nr0) = ltag(:nr0)

#endif

          END IF


          DO ij=1, nr0


            isea   = tmpint(iloc,ij)

            CALL init_get_jsea_isproc_glob(isea, jj, jsea, ispro2)

            IF ( flgrec ) hgstge(gdst,jj)%IMPSRC(irec,ij) = ispro2


            IF ( ispro2 .EQ. improc ) THEN

              IF ( tmplog(iloc) ) THEN

                i3(jj) = i3(jj) + 1

                isnd   = i3(jj)

              ELSE

                i4(jj) = i4(jj) - 1

                isnd   = i4(jj)

              END IF

              hgstge(gdst,jj)%ISEND(isnd,1) = jsea

#ifdef W3_DIST

              hgstge(gdst,jj)%ISEND(isnd,2) = isproc

#endif

              hgstge(gdst,jj)%ISEND(isnd,3) = irec

              hgstge(gdst,jj)%ISEND(isnd,4) = ij

#ifdef W3_DIST

              hgstge(gdst,jj)%ISEND(isnd,5) = ltag(ij)

#endif

            END IF


          END DO

          !

#ifdef W3_DIST

          ltag   = ltag + nr0

          ltag0  = ltag0 + nr0

#endif

          !

        END DO

        !

        ! 2.e Adjust FLAGST using MAPTST

        !

#ifdef W3_T

        ALLOCATE ( mapst(ny,nx), stat=istat )

        check_alloc_status( istat )

        mapst  = '-'

#endif

        !

        DO isea=1, nsea

          idst     = mapsf(isea,1)

          jdst     = mapsf(isea,2)

          IF ( maptst(jdst,idst) .GT. 0 ) flagst(isea) = .NOT. flghg1

#ifdef W3_T

          IF ( flagst(isea) ) THEN

            mapst(jdst,idst)  = 'O'

          ELSE

            mapst(jdst,idst)  = 'X'

          END IF

#endif

        END DO

        !

        ! 2.f Test output map

        !

#ifdef W3_T

        WRITE (mdst,9025) 'MAPTST'

        DO jdst=ny,1 , -1

          WRITE (mdst,9026) maptst(jdst,:) + 88*inflnd(jdst,:)

        END DO

#endif

        !

#ifdef W3_T

        WRITE (mdst,9025) 'MAPSTA'

        DO jdst=ny,1 , -1

          WRITE (mdst,9026) mapsta(jdst,:)

        END DO

#endif

        !

#ifdef W3_T

        WRITE (mdst,9025) 'MAPST2'

        DO jdst=ny,1 , -1

          WRITE (mdst,9026) mapst2(jdst,:)

        END DO

#endif

        !

#ifdef W3_T

        WRITE (mdst,9025) 'FLAGST'

        DO jdst=ny,1 , -1

          WRITE (mdst,9027) mapst(jdst,:)

        END DO

#endif

        !

        DEALLOCATE ( maptst, inflnd, stat=istat )

        check_dealloc_status( istat )

#ifdef W3_T

        DEALLOCATE ( mapst, stat=istat )

        check_dealloc_status( istat )

#endif

        !

        ! 2.g Test output receiving

        !

#ifdef W3_T3

        DO gsrc=1, nrgrd

          nr0    = hgstge(gdst,gsrc)%NREC

          IF ( nr0 .EQ. 0 ) THEN

            WRITE (mdst,9030) gsrc

          ELSE

            WRITE (mdst,9031) gsrc, nr0

            DO irec=1, nr0

              jsea   = hgstge(gdst,gsrc)%LJSEA(irec)

              nrtot  = hgstge(gdst,gsrc)%NRAVG(irec)

              IF ( nrtot .LE. 15 ) THEN

                WRITE (mdst,9032) jsea, nrtot,           &

                     hgstge(gdst,gsrc)%WGTH(irec,:nrtot)

              ELSE

                WRITE (mdst,9032) jsea, nrtot,           &

                     hgstge(gdst,gsrc)%WGTH(irec,1:15)

                WRITE (mdst,9033)                        &

                     hgstge(gdst,gsrc)%WGTH(irec,16:nrtot)

              END IF

              WRITE (mdst,9034)                            &

                   hgstge(gdst,gsrc)%IMPSRC(irec,1:nrtot)

              WRITE (mdst,9034)                            &

                   hgstge(gdst,gsrc)%ITAG(irec,1:nrtot)

            END DO

          END IF

        END DO

#endif

        !

        ! 2.h Test output sending

        !

#ifdef W3_T4

        DO gsrc=1, nrgrd

          nr0    = hgstge(gdst,gsrc)%NSND

          IF ( nr0 .EQ. 0 ) THEN

            WRITE (mdst,9040) gsrc

          ELSE

            WRITE (mdst,9041) gsrc, nr0

            DO isnd=1, nr0

              WRITE (mdst,9042) hgstge(gdst,gsrc)%ISEND(isnd,:)

            END DO

          END IF

        END DO

#endif

        !

        ! 2.i Final clean up

        !

        DEALLOCATE ( idstl, idsth, jdstl, jdsth, gridok,  bdist, &

             tmpint, tmprl, tmplog, stat=istat )

        check_dealloc_status( istat )


        IF (old_method) THEN

          DEALLOCATE ( bdist_om, tmpint_om, tmprl_om, stat=istat )

          check_dealloc_status( istat )

        END IF


#ifdef W3_DIST

        DEALLOCATE ( ltag, stat=istat )

        check_dealloc_status( istat )

#endif

        !


        ! Notes: We are done with this dst (low rank) grid, so we deallocate ALLWGTS .

        !        This is important because ALLWGTS will be allocated again for the next

        !        dst grid.


#ifdef W3_SCRIP

        DO jj=1, grdhgh(gdst,0)

          gsrc   = grdhgh(gdst,jj)

          DO kdst=1,dst_grid_size

            DEALLOCATE ( allwgts(gsrc)%WGTDATA(kdst)%W, stat=istat )

            check_dealloc_status( istat )

            DEALLOCATE ( allwgts(gsrc)%WGTDATA(kdst)%K, stat=istat )

            check_dealloc_status( istat )

          END DO

          DEALLOCATE ( allwgts(gsrc)%WGTDATA, stat=istat )

          check_dealloc_status( istat )

        END DO

        DEALLOCATE ( allwgts, stat=istat )

        check_dealloc_status( istat )

#endif


      END IF !  IF ( GRDHGH(GDST,0) ...

#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(2)

      WRITE(nmyout,*) "WMGHGH, LOOP LOWRANK_GRID, GDST= ", gdst, " TOOK ", elapsed_time, " MSEC"

#endif

    END DO lowrank_grid


    ! If SCRIPNC and L_STOP, then we are done

    IF ( lscripnc .AND. l_stop ) THEN

      ! WW3 processes wait here till all have finished

#ifdef W3_MPI

      CALL mpi_barrier( mpi_comm_mwave, ierr_mpi )

#endif

      ! This is not a true error, so exit code is zero

      WRITE( mdse, '(A,I4.4,A)' ) 'IMPROC=',improc, &

           ': STOP_SCRIP option invoked: '// &

           'non-error exit after writing remap netcdf files'

      CALL extcde( 0 )

    END IF


    DEALLOCATE ( i1, i2, i3, i4, stat=istat )

    check_dealloc_status( istat )

#ifdef W3_MPIBDI

    DEALLOCATE ( nx_size, irq, mstat, stat=istat )

    check_dealloc_status( istat )

#endif

    DEALLOCATE ( nx_beg, nx_end, stat=istat )

    check_dealloc_status( istat )


    !

    ! 2.j Test output counters

    !

#ifdef W3_T

    WRITE (mdst,9028) 'NTOT'

    DO jj=1, nrgrd

      WRITE (mdst,9029) hgstge(jj,:)%NTOT

    END DO

#endif

    !

#ifdef W3_T

    WRITE (mdst,9028) 'NREC'

    DO jj=1, nrgrd

      WRITE (mdst,9029) hgstge(jj,:)%NREC

    END DO

#endif

    !

#ifdef W3_T

    WRITE (mdst,9028) 'NRC1'

    DO jj=1, nrgrd

      WRITE (mdst,9029) hgstge(jj,:)%NRC1

    END DO

#endif

    !

#ifdef W3_T

    WRITE (mdst,9028) 'NSND'

    DO jj=1, nrgrd

      WRITE (mdst,9029) hgstge(jj,:)%NSND

    END DO

#endif

    !

#ifdef W3_T

    WRITE (mdst,9028) 'NSN1'

    DO jj=1, nrgrd

      WRITE (mdst,9029) hgstge(jj,:)%NSN1

    END DO

#endif

    !

#ifdef W3_T

    WRITE (mdst,9028) 'NSMX'

    DO jj=1, nrgrd

      WRITE (mdst,9029) hgstge(jj,:)%NSMX

    END DO

#endif

    !

#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(1)

    WRITE(nmyout,*) "WMGHGH, ALL TOOK ", elapsed_time, " MSEC"

#endif


    RETURN

    !

    ! Formats

    !

1000 FORMAT (/' *** WAVEWATCH III ERROR IN WMGHGH : *** '/           &

         '     GRDHGH NOT YET ALLOCATED, CALL WMGLOW FIRST'/)

1020 FORMAT (/' *** WAVEWATCH III ERROR IN WMGHGH : *** '/           &

         '     TMPINT AND TMPRL TOO SMALL (w/out SCRIP)'/)

#ifdef W3_SCRIP

1021 FORMAT (/' *** WAVEWATCH III ERROR IN WMGHGH : *** '/           &

         '     TMPINT AND TMPRL TOO SMALL (w/SCRIP) '/)

#endif

    !

#ifdef W3_T

9010 FORMAT ( ' TEST WMGHGH : INITIALIZE BOUNDARY DISTANCE MAPS')

9011 FORMAT ( '               GRID = ',i3,'  RANK = ',i3,         &

         '  NBI = ',i6)

9012 FORMAT ( '                  *** MAP NOT NEEDED ***')

9013 FORMAT ( ' TEST WMGHGH : FINAL MAP ')

9014 FORMAT (2x,65i2)

#endif

    !/

#ifdef W3_T

9020 FORMAT ( ' TEST WMGHGH : GRID',i3,' HAS',i3,' DATA SOURCES')

9021 FORMAT ( '               NO PROCESSING REQUIRED')

9022 FORMAT ( ' TEST WMGHGH : GRID',i3,' COVERS ',4i8)

9023 FORMAT ( ' TEST WMGHGH : GRID',i3,                           &

         ', NR OF POINTS TO PROCESS:',i5)

9025 FORMAT ( ' TEST WMGHGH : FINAL ',a)

9026 FORMAT (2x,65i2)

9027 FORMAT (2x,65a2)

#endif

    !

#ifdef W3_T

9028 FORMAT ( ' TEST WMGHGH : COUNTERS ',a)

9029 FORMAT (2x,20i6)

#endif

    !

#ifdef W3_T3

9030 FORMAT ( ' TEST WMGHG : FROM GRID',i3,', NO DATA TO RECEIVE')

9031 FORMAT ( ' TEST WMGHG : FROM GRID',i3,', RECEIVING ',i6)

9032 FORMAT ( 2x,i10,i6,15f6.2)

9033 FORMAT ( 18x,15f6.2)

9034 FORMAT ( 18x,15i6)

#endif

    !

#ifdef W3_T4

9040 FORMAT ( ' TEST WMGHG : FROM GRID',i3,', NO DATA TO SEND')

9041 FORMAT ( ' TEST WMGHG : FROM GRID',i3,', SENDING ',i6)

9042 FORMAT ( 12x,i10,4i6)

#endif

    !/

    !/ End of WMGHGH ----------------------------------------------------- /

    !/

  END SUBROUTINE wmghgh

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

  SUBROUTINE wmgeql

    !/

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

    !/                  | WAVEWATCH III           NOAA/NCEP |

    !/                  |           H. L. Tolman            |

    !/                  |                        FORTRAN 90 |

    !/                  | Last update :         10-Dec-2014 !

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

    !/

    !/    24-Apr-2006 : Origination.                        ( version 3.09 )

    !/    23-Dec-2006 : Adding group test.                  ( version 3.10 )

    !/    28-Dec-2006 : Simplify NIT for partial comm.      ( version 3.10 )

    !/    22-Jan-2007 : Add saving og NAVMAX.               ( version 3.10 )

    !/    02-Feb-2007 : Setting FLAGST for replaced points. ( version 3.10 )

    !/    15-Feb-2007 : Tweaking MAPODI algorithm in WMGEQL.( version 3.10 )

    !/    11-Apr-2008 : Big fix active edges (MAPSTA=2)     ( version 3.13 )

    !/    14-Apr-2008 : Big fix for global grids.           ( version 3.13 )

    !/    20-May-2009 : Linking FLAGST and FLGHG1.          ( version 3.14 )

    !/    30-Oct-2009 : Implement run-time grid selection.  ( version 3.14 )

    !/                  (W. E. Rogers & T. J. Campbell, NRL)

    !/    06-Dec-2010 : Change from GLOBAL (logical) to ICLOSE (integer) to

    !/                  specify index closure for a grid.   ( version 3.14 )

    !/                  (T. J. Campbell, NRL)

    !/    23-Dec-2010 : Fix HPFAC and HQFAC by including the COS(YGRD)

    !/                  factor with DXDP and DXDQ terms.    ( version 3.14 )

    !/                  (T. J. Campbell, NRL)

    !/    05-Aug-2013 : Change PR2/3 to UQ/UNO in distances.( version 4.12 )

    !/    10-Dec-2014 : Add checks for allocate status      ( version 5.04 )

    !/    28-Oct-2020 : Add SMCTYPE for SMC sub-grid.  JGLi ( version 6.xx )

    !/

    !  1. Purpose :

    !

    !     Determine relations to same ranked grids for each grid.

    !

    !  2. Method :

    !

    !     Cross mapping of grid points, determine boundary distance data

    !     and interpolation weights.

    !

    !  3. Parameters :

    !

    !  4. Subroutines used :

    !

    !      Name      Type  Module   Description

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

    !      W3SETG, W3SETO, WMSETM

    !                Subr. W3GDATMD Manage data structures.

    !      STRACE    Subr. W3SERVMD Subroutine tracing.

    !      EXTCDE    Subr.   Id.    Program abort.

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

    !

    !  5. Called by :

    !

    !  6. Error messages :

    !

    !  7. Remarks :

    !

    !     - In looking for compatable boundary points in overlapping grids

    !       two assumptions hav been made.

    !        a) No active boundary points exist in global grids.

    !        b) For a lower resolution grid an expanded sewarch area is

    !           required for corresponding active grid points. By limiting

    !           the resolution ratio to 2, only one extra grid point needs

    !           to be considered (JXL2 versus JXL etc.).

    !

    !  8. Structure :

    !

    !  9. Switches :

    !

    !     !/PRn  Propagation scheme.

    !

    !     !/O12  Removed boundary points output (central).

    !     !/O13  Removed boundary points output (edge).


    !     !/S    Enable subroutine tracing.

    !     !/T    Enable test output.

    !     !/T5   Detailed test output 'receiving'.

    !     !/T6   Detailed test output 'sending'.

    !     !/T7   Detailed test output all.

    !

    !     !/MPI  Distribbuted memory management.

    !

    ! 10. Source code :

    !

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

    !

    USE constants

    USE w3gdatmd

    USE w3odatmd

    USE w3adatmd

    USE wmmdatmd

    !

    USE w3servmd, ONLY: extcde

    !      USE W3PARALL, ONLY: INIT_GET_JSEA_ISPROC_GLOB

#ifdef W3_S

    USE w3servmd, ONLY: strace

#endif

    !

    IMPLICIT NONE

    !/

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

    !/ Parameter list

    !/

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

    !/ Local parameters

    !/

    INTEGER                  :: I, J, IX, IXL, IXH, IY, IYL, IYH,   &

         JX, JXL, JXH, JXL2, JXH2,           &

         JY, JYL, JYH, JYL2, JYH2,           &

         NR, NT, NA, NTL, JJ, NIT, NG, NOUT, &

         ISEA, JSEA, ISPROC, ITAG, TGRP,     &

         EXTRA, IP, NP

#ifdef W3_T7

    INTEGER                  :: IA

#endif

#ifdef W3_S

    INTEGER, SAVE            :: IENT = 0

#endif

    INTEGER, ALLOCATABLE     :: MAP3D(:,:,:), NREC(:), NSND(:),     &

         NTPP(:), MAPOUT(:,:)

    REAL                     :: FACTOR, XSL, XSH, YSL, YSH, XA, YA, &

         XR, YR, RX(2), RY(2), STX, STY,     &

         STXY, NEWVAL, WGTH

    REAL, PARAMETER          :: TODO   = -9.99e25

    REAL, PARAMETER          :: ODIMAX = 25.

    REAL, PARAMETER          :: FACMAX = 2.001

    REAL, ALLOCATABLE        :: WGT3D(:,:,:)

    LOGICAL                  :: CHANGE, XEXPND, YEXPND

    LOGICAL, ALLOCATABLE     :: SHRANK(:,:), DOGRID(:),             &

         MASKA(:,:), MASKI(:,:)

#ifdef W3_T5

    CHARACTER(LEN=18), ALLOCATABLE :: TSTR(:)

    CHARACTER(LEN=18)              :: DSTR

#endif

    !

    TYPE store

      INTEGER                :: NTOT, NFIN

      INTEGER, POINTER       :: IX(:), IY(:), NAV(:), ISS(:,:),     &

           JSS(:,:), IPS(:,:), ITG(:,:)

      REAL, POINTER          :: AWG(:,:)

      LOGICAL, POINTER       :: FLA(:)

      LOGICAL                :: INIT

    END TYPE store

    !

    TYPE(store), ALLOCATABLE :: STORES(:,:)

    !/

#ifdef W3_S

    CALL strace (ient, 'WMGEQL')

#endif

    !

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

    ! 0.  Initializations

    !


    ALLOCATE ( shrank(nrgrd,nrgrd), stores(nrgrd,nrgrd),            &

         dogrid(nrgrd), stat=istat )

    check_alloc_status( istat )

    !

    shrank = .false.

    !

    DO i=1, nrgrd


      DO j=1, nrgrd

        stores(i,j)%INIT = .false.

        stores(i,j)%NTOT = 0

        stores(i,j)%NFIN = 0

      END DO

    END DO

    !

    IF ( flagll ) THEN

      factor = radius * dera

      !notes: was FACTOR = RADIUS / 360. (I don't know where this came from....

      !       ...maybe it was supposed to be CIRCUMFERENCE/360)

    ELSE

      factor = 1.

    END IF

    itag   = 0

    !

#ifdef W3_SMC

    !!  Check GTYPE for all grids.

    IF( improc.EQ.nmperr )  WRITE (mdse,*) " WMGEQL GTYPE:",   &

         ( grids(i)%GTYPE, i=1, nrgrd )

#endif

    !

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

    ! 1.  Grid point relations and temp data storage

    ! 1.a Outer loop over all grids

    !

#ifdef W3_T

    WRITE (mdst,9010)

#endif

    !

    DO i=1, nrgrd


#ifdef W3_T

      WRITE (mdst,9011) i, grank(i)

#endif

      !

      ! 1.b Find grids with same rank

      !

      nr     = 0

      !

#ifdef W3_SMC

      !!  SMC grids use WMSMCEQL for equal ranked grids.  JGLi23Mar2021

      IF( grids(i)%GTYPE .EQ. smctype ) cycle

#endif

      !

      DO j=1, nrgrd


        IF ( grank(i).NE.grank(j) .OR. i.EQ.j ) cycle

        !

#ifdef W3_SMC

        IF( grids(j)%GTYPE .EQ. smctype ) cycle

#endif

        !

#ifdef W3_T

        WRITE (mdst,9012) j

#endif

        shrank(i,j) = .true.

        nr          = nr + 1

      END DO

      !

      CALL w3setg ( i, mdse, mdst )

      !

      dogrid(i) = nr .GT. 0


      !..notes: we will reach this point even if there are no equal rank grids


#ifdef W3_T

      IF ( nr .EQ. 0 ) WRITE (mdst,9013) 'NO GRIDS WITH SAME RANK'

#endif

      IF ( nr .EQ. 0 ) cycle


      !..notes: we will not reach this point if are no equal rank grids. that makes it a good place to check against grid type


      IF ( iclose .EQ. iclose_trpl ) THEN

        IF ( improc.EQ.nmperr ) WRITE(mdse,*)'SUBROUTINE WMGEQL IS'// &

             ' NOT YET ADAPTED FOR TRIPOLE GRIDS. STOPPING NOW.'

        CALL extcde ( 1 )

      END IF


      ! Unresolved bug: this triggers even for 2 irregular grids that are not overlapping!

      !   We should only be checking for cases of 2 irregular grids of equal rank that

      !   are overlapping. Unfortunately, at this point in the routine, we don't know

      !   whether they are overlapping...requires more code to do this, since all code

      !   in this routine is for regular grids. Fortunately, there is really no

      !   disadvantage to making the two irregular grids to be different rank using

      !   ww3_multi.inp


      IF ( gtype .EQ. ungtype ) THEN

        IF ( improc.EQ.nmperr )WRITE (mdse,'(/3A)') ' *** ERROR ',  &

             'WMGEQL: UNSTRUCTURED GRID SUPPORT NOT YET ',          &

             'IMPLEMENTED ***'

        CALL extcde ( 999 )

      END IF

      IF ( gtype .EQ. clgtype ) THEN

        IF ( improc.EQ.nmperr )WRITE (mdse,'(/3A)') ' *** ERROR ',  &

             'WMGEQL: CURVILINEAR GRID SUPPORT NOT IMPLEMENTED ',   &

             'FOR NRGRD > 1 ***'

        CALL extcde ( 999 )

      END IF


      !

      ! 1.c Fill TMPMAP with raw relational data

      !

      ! 1.c.1 Loop over grids, select same rank

      !

      DO j=1, nrgrd


        IF ( .NOT. shrank(i,j) ) cycle

        !

        ! 1.c.2 Determine shared area

        !       Don't even try for X in LLG

        !


        !       Note: Check is against FLAGLL. Would it be more appropriate

        !                   to check against ICLOSE?

        IF ( flagll ) THEN

          ixl    =  1

          ixh    = nx

        ELSE

          xsl    = ( grids(j)%X0 - x0 ) / sx - 0.01

          xsh    = ( grids(j)%X0 + grids(j)%SX*(grids(j)%NX-1)   &

               - x0 ) / sx + 0.01

          ixl    = max( 1+nint(xsl) ,  1 )

          ixh    = min( 1+nint(xsh) , nx )

        END IF

        !

        ysl    = ( grids(j)%Y0 - y0 ) / sy - 0.01

        ysh    = ( grids(j)%Y0 + grids(j)%SY*(grids(j)%NY-1)   &

             - y0 ) / sy + 0.01

        iyl    = max( 1+nint(ysl) ,  1 )

        iyh    = min( 1+nint(ysh) , ny )

        !

        nt     = (1+ixh-ixl) * (1+iyh-iyl)

        IF ( nt .EQ. 0 ) cycle

        !

        stores(i,j)%INIT = .true.

        ALLOCATE ( stores(i,j)%IX(nt)   , stores(i,j)%IY(nt)   ,    &

             stores(i,j)%NAV(nt)  , stores(i,j)%FLA(nt)  ,    &

             stores(i,j)%ISS(nt,4), stores(i,j)%JSS(nt,4),    &

             stores(i,j)%IPS(nt,4), stores(i,j)%ITG(nt,4),    &

             stores(i,j)%AWG(nt,4), stat=istat )

        check_alloc_status( istat )

        stores(i,j)%NAV  = 0

        stores(i,j)%FLA  = .false.

        stores(i,j)%ISS  = 0

        stores(i,j)%JSS  = 0

        stores(i,j)%IPS  = 0

        stores(i,j)%ITG  = 0

        stores(i,j)%AWG  = 0.

        !

        ! 1.c.3 Loops over shared area

        !

        nt     = 0

        !

        xexpnd = sx .GT. grids(j)%SX

        yexpnd = sy .GT. grids(j)%SY

        !

        DO ix=ixl, ixh

          xa     = x0 + real(ix-1)*sx

          IF ( flagll ) THEN

            xr     = 1. + mod(1080. + xa - grids(j)%X0 , 360. ) &

                 / grids(j)%SX

          ELSE

            xr     = 1. + (xa-grids(j)%X0) / grids(j)%SX

          END IF

          jxl    = int(xr)

          jxh    = jxl + 1

          rx(1)     = 1. - mod(xr,1.)

          IF ( rx(1).GT.0.99 .OR. jxh.EQ.grids(j)%NX+1 ) THEN

            jxh    = jxl

            rx(1)  = 1.

          END IF

          IF ( rx(1).LT.0.01 .OR. jxl.EQ.0 ) THEN

            jxl    = jxh

            rx(1)  = 1.

          END IF

          rx(2)  = 1. - rx(1)

          !

          IF ( jxl.LT.1 .OR. jxh.GT.grids(j)%NX ) cycle

          !

          IF ( xexpnd ) THEN

            jxl2   = max( 1 , jxl-1 )

            jxh2   = min( grids(j)%NX , jxh+1 )

          ELSE

            jxl2   = jxl

            jxh2   = jxh

          END IF

          !

          DO iy=iyl, iyh

            ya     = y0 + real(iy-1)*sy

            yr     = 1. + (ya-grids(j)%Y0) / grids(j)%SY

            jyl    = int(yr)

            jyh    = jyl + 1

            ry(1)  = 1. - mod(yr,1.)

            IF ( ry(1).GT.0.99 .OR. jyh.EQ.grids(j)%NY+1 ) THEN

              jyh    = jyl

              ry(1)  = 1.

            END IF

            IF ( ry(1).LT.0.01 .OR. jyl.EQ.0 ) THEN

              jyl    = jyh

              ry(1)  = 1.

            END IF

            IF ( ry(1) .GT. 0.99 ) jyh = jyl

            ry(2)  = 1. - ry(1)

            !

            IF ( jyl.LT.1 .OR. jyh.GT.grids(j)%NY ) cycle

            !

            IF ( yexpnd ) THEN

              jyl2   = max( 1 , jyl-1 )

              jyh2   = min( grids(j)%NY , jyh+1 )

            ELSE

              jyl2   = jyl

              jyh2   = jyh

            END IF

            !

            ! 1.c.4 Temp storage of raw data

            !

            nt     = nt + 1

            na     = 0

#ifdef W3_SHRD

            isproc = 1

#endif

            stores(i,j)%IX(nt)  = ix

            stores(i,j)%IY(nt)  = iy

            !

            DO jx = jxl, jxh

              DO jy = jyl, jyh

                IF ( grids(j)%MAPSTA(jy,jx) .NE. 0 ) THEN

                  na     = na + 1

                  itag   = itag + 1

                  wgth   = rx(1+jx-jxl) * ry(1+jy-jyl)

                  isea   = grids(j)%MAPFS(jy,jx)

                  IF ( isea .EQ. 0 ) THEN

                    jsea   = 0

#ifdef W3_MPI

                    isproc = 1

#endif

                  ELSE

                    CALL init_get_jsea_isproc_glob(isea, j, jsea, isproc)

                  END IF

                  stores(i,j)%AWG(nt,na) = wgth

                  stores(i,j)%ISS(nt,na) = isea

                  stores(i,j)%JSS(nt,na) = jsea

                  stores(i,j)%IPS(nt,na) = isproc

                  stores(i,j)%ITG(nt,na) = itag

                END IF

              END DO

            END DO

            !

            DO jx = jxl2, jxh2

              DO jy = jyl2, jyh2

                IF ( abs(grids(j)%MAPSTA(jy,jx)) .EQ. 2 )           &

                     stores(i,j)%FLA(nt) = .true.

              END DO

            END DO

            !

            wgth   = sum( stores(i,j)%AWG(nt,1:na) )

            IF ( wgth .LT. 0.499 ) THEN

              na    = 0

            ELSE

              stores(i,j)%AWG(nt,:) = stores(i,j)%AWG(nt,:) / wgth

            END IF

            !

            stores(i,j)%NAV(nt) = na

            !

            ! ... End of loops in 1.c

            !

          END DO

        END DO

        !

        stores(i,j)%NTOT = nt

        !

      END DO

      !

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

      ! 2.  Generate open edge distance maps

      ! 2.a Base map based on MAPSTA only, time step not included.

      !

#ifdef W3_T

      WRITE (mdst,9020) i

#endif

      !

      ALLOCATE ( mdatas(i)%MAPODI(ny,nx), stat=istat )

      check_alloc_status( istat )

      mapodi => mdatas(i)%MAPODI

      mapodi = 0.

      !

      DO ix=1, nx

        DO iy=1, ny

          IF ( abs(mapsta(iy,ix)) .EQ. 1 ) THEN

            mapodi(iy,ix) = todo

          ELSE IF ( abs(mapsta(iy,ix)) .EQ. 2 ) THEN

            mapodi(iy,ix) = -2. / sig(1) * dtmax

          ELSE

            mapodi(iy,ix) = -1. / sig(1) * dtmax

          END IF

        END DO

      END DO

      !

      ! 2.b Add in cross-grid information from STORES

      !

      ALLOCATE ( maska(ny,nx), stat=istat )

      check_alloc_status( istat )

      maska  = .false.

      !

      DO j=1, nrgrd

        IF ( .NOT. shrank(i,j) ) cycle

        DO jj=1, stores(i,j)%NTOT

          ix     = stores(i,j)%IX(jj)

          iy     = stores(i,j)%IY(jj)

          IF ( ix.EQ.1 .OR. ix.EQ.nx .OR. iy.EQ.1 .OR. iy.EQ.ny ) THEN

            maska(iy,ix) = stores(i,j)%FLA(jj) .OR.               &

                 stores(i,j)%NAV(jj).EQ.0

            IF ( abs(mapsta(iy,ix)).EQ.2 .AND.                    &

                 .NOT.stores(i,j)%FLA(jj) .AND.                   &

                 stores(i,j)%NAV(jj).GT.0 ) THEN

              mapodi(iy,ix) = 0.

#ifdef W3_O13

              IF ( improc.EQ.nmperr )                      &

                   WRITE (mdse,1020) i, ix, 1

#endif

            END IF

          ELSE

            maska(iy,ix) = stores(i,j)%FLA(jj)

          END IF

          IF ( mapsta(iy,ix).EQ.0 .AND. mapst2(iy,ix) .EQ.1 .AND.   &

               stores(i,j)%NAV(jj).GT.0 ) mapodi(iy,ix) = 0.

        END DO

      END DO

      !

      ! 2.c Remove incompatable boundary points

      !

      ALLOCATE ( maski(ny,nx), stat=istat )

      check_alloc_status( istat )

      maski  = .false.

      !

      DO ix=2, nx-1

        DO iy=2, ny-1

          IF ( abs(mapsta(iy,ix)) .EQ. 2 .AND.                      &

               .NOT. maska(iy,ix) .AND. (                           &

               mapodi(iy-1,ix  ) .GE. 0. .OR.              &

               mapodi(iy+1,ix  ) .GE. 0. .OR.              &

               mapodi(iy  ,ix-1) .GE. 0. .OR.              &

               mapodi(iy  ,ix+1) .GE. 0. ) ) THEN

            maski(iy,ix) = .true.

#ifdef W3_O12

            IF ( improc.EQ.nmperr ) WRITE (mdse,1020) i, ix, iy

#endif

          END IF

        END DO

      END DO

      !

      DEALLOCATE ( maska, stat=istat )

      check_dealloc_status( istat )

      !

      DO ix=1, nx

        DO iy=1, ny

          IF ( maski(iy,ix) ) mapodi(iy,ix) = 0.

        END DO

      END DO

      !

      ! 2.d Mask out influenced edge

      !

#ifdef W3_PR0

      nit    = 0

#endif

#ifdef W3_PR1

      nit    = ( 1 + int(dtmax/dtcfl-0.001) ) * 1

#endif

#ifdef W3_UQ

      nit    = ( 1 + int(dtmax/dtcfl-0.001) ) * 3

#endif

#ifdef W3_UNO

      nit    = ( 1 + int(dtmax/dtcfl-0.001) ) * 3

#endif

      !

      IF ( iclose.NE.iclose_none ) THEN

        ixl    =  1

        ixh    = nx

      ELSE

        ixl    =  2

        ixh    = nx - 1

      END IF

      !

      DO j=1, nit

        !

        maski  = .false.

        !

        DO ix=ixl, ixh

          IF ( ix .EQ. 1 ) THEN

            jxl    = nx

            jxh    =  2

          ELSE IF ( ix .EQ. nx ) THEN

            jxl    = nx - 1

            jxh    = 1

          ELSE

            jxl    = ix - 1

            jxh    = ix + 1

          END IF

          !

          DO iy=2, ny-1

            IF ( mapodi(iy,ix) .EQ. todo .AND. (                    &

                 mapodi(iy+1,ix ) .GE. 0. .OR.                      &

                 mapodi(iy  ,jxl) .GE. 0. .OR.                      &

                 mapodi(iy-1,ix ) .GE. 0. .OR.                      &

                 mapodi(iy  ,jxh) .GE. 0. .OR.                      &

                 ( mapodi(iy+1,jxh) .GE. 0. .AND. .NOT.               &

                 ( mapsta(iy+1,ix ) .NE. 1 .AND.                      &

                 mapsta(iy  ,jxh) .NE. 1 ) ) .OR.                   &

                 ( mapodi(iy+1,jxl) .GE. 0. .AND. .NOT.               &

                 ( mapsta(iy+1,ix ) .NE. 1 .AND.                      &

                 mapsta(iy  ,jxl) .NE. 1 ) ) .OR.                   &

                 ( mapodi(iy-1,jxl) .GE. 0. .AND. .NOT.               &

                 ( mapsta(iy-1,ix ) .NE. 1 .AND.                      &

                 mapsta(iy  ,jxl) .NE. 1 ) ) .OR.                   &

                 ( mapodi(iy-1,jxh) .GE. 0. .AND. .NOT.               &

                 ( mapsta(iy-1,ix ) .NE. 1 .AND.                      &

                 mapsta(iy  ,jxh) .NE. 1 ) ) ) )                    &

                 maski(iy,ix) = .true.

          END DO

          !

        END DO

        !

        DO ix=ixl, ixh

          DO iy=2, ny-1

            IF ( maski(iy,ix) ) mapodi(iy,ix) = 0.

          END DO

        END DO

        !

      END DO

      !

      ! 2.e Compute distances

      !

      DO

        maski  = .false.

        !

        DO ix=ixl, ixh

          IF ( ix .EQ. 1 ) THEN

            jxl    = nx

            jxh    =  2

          ELSE IF ( ix .EQ. nx ) THEN

            jxl    = nx - 1

            jxh    = 1

          ELSE

            jxl    = ix - 1

            jxh    = ix + 1

          END IF

          DO iy=2, ny-1

            IF ( mapodi(iy,ix) .EQ. todo .AND. (                    &

                 mapodi(iy+1,ix ) .GE. 0. .OR.                      &

                 mapodi(iy-1,ix ) .GE. 0. .OR.                      &

                 mapodi(iy  ,jxh) .GE. 0. .OR.                      &

                 mapodi(iy  ,jxl) .GE. 0. .OR.                      &

                 ( mapodi(iy+1,jxh) .GE. 0. .AND. .NOT.               &

                 ( mapsta(iy+1,ix ) .NE. 1 .AND.                      &

                 mapsta(iy  ,jxh) .NE. 1 ) ) .OR.                   &

                 ( mapodi(iy+1,jxl) .GE. 0. .AND. .NOT.               &

                 ( mapsta(iy+1,ix ) .NE. 1 .AND.                      &

                 mapsta(iy  ,jxl) .NE. 1 ) ) .OR.                   &

                 ( mapodi(iy-1,jxl) .GE. 0. .AND. .NOT.               &

                 ( mapsta(iy-1,ix ) .NE. 1 .AND.                      &

                 mapsta(iy  ,jxl) .NE. 1 ) ) .OR.                   &

                 ( mapodi(iy-1,jxh) .GE. 0. .AND. .NOT.               &

                 ( mapsta(iy-1,ix ) .NE. 1 .AND.                      &

                 mapsta(iy  ,jxh) .NE. 1 ) ) ) )                    &

                 maski(iy,ix) = .true.

          END DO

        END DO

        !

        change = .false.

        DO iy=2, ny-1

          DO ix=ixl, ixh

            IF ( ix .EQ. 1 ) THEN

              jxl    = nx

              jxh    =  2

            ELSE IF ( ix .EQ. nx ) THEN

              jxl    = nx - 1

              jxh    = 1

            ELSE

              jxl    = ix - 1

              jxh    = ix + 1

            END IF

            isea   = mapfs(iy,ix)

            sty    = factor * hqfac(iy,ix) / ( 0.58 * grav )

            stx    = factor * hpfac(iy,ix) &

                 / ( 0.58 * grav )

            stxy   = sqrt( stx**2 + sty**2 )

            IF ( maski(iy,ix) ) THEN

              newval = odimax / sig(1) * dtmax

              IF ( mapodi(iy+1,ix ).GE.0. .AND. .NOT.             &

                   maski(iy+1,ix ) ) newval = min(              &

                   newval , mapodi(iy+1,ix )+sty )

              IF ( mapodi(iy-1,ix ).GE.0. .AND. .NOT.             &

                   maski(iy-1,ix ) ) newval = min(              &

                   newval , mapodi(iy-1,ix )+sty )

              IF ( mapodi(iy  ,jxh).GE.0. .AND. .NOT.             &

                   maski(iy  ,jxh) ) newval = min(              &

                   newval , mapodi(iy  ,jxh)+stx)

              IF ( mapodi(iy  ,jxl).GE.0. .AND. .NOT.             &

                   maski(iy  ,jxl) ) newval = min(              &

                   newval , mapodi(iy  ,jxl)+stx)

              IF ( mapodi(iy+1,jxh).GE.0. .AND. .NOT.             &

                   ( mapsta(iy+1,ix ) .NE. 1 .AND.                  &

                   mapsta(iy  ,jxh) .NE. 1 ) ) newval =           &

                   min( newval , mapodi(iy+1,jxh)+stxy)

              IF ( mapodi(iy+1,jxl).GE.0. .AND. .NOT.             &

                   ( mapsta(iy+1,ix ) .NE. 1 .AND.                  &

                   mapsta(iy  ,jxl) .NE. 1 ) ) newval =           &

                   min( newval , mapodi(iy+1,jxl)+stxy)

              IF ( mapodi(iy-1,jxl).GE.0. .AND. .NOT.             &

                   ( mapsta(iy-1,ix ) .NE. 1 .AND.                  &

                   mapsta(iy  ,jxl) .NE. 1 ) ) newval =           &

                   min( newval , mapodi(iy-1,jxl)+stxy)

              IF ( mapodi(iy-1,jxh).GE.0. .AND. .NOT.             &

                   ( mapsta(iy-1,ix ) .NE. 1 .AND.                  &

                   mapsta(iy  ,jxh) .NE. 1 ) ) newval =           &

                   min( newval , mapodi(iy-1,jxh)+stxy)

              mapodi(iy,ix) = newval

              change = .true.

            END IF

          END DO

        END DO

        !

        IF ( .NOT. change ) EXIT

      END DO

      !

      DO ix=2, nx-1

        DO iy=2, ny-1

          IF ( mapodi(iy,ix) .EQ. todo )                            &

               mapodi(iy,ix) = 2. * odimax / sig(1) * dtmax

        END DO

      END DO

      !

      DEALLOCATE ( maski, stat=istat )

      check_dealloc_status( istat )

      !

      ! 2.f Update FLAGST

      !

      DO isea=1, nsea

        ix     = mapsf(isea,1)

        iy     = mapsf(isea,2)

        IF ( mapodi(iy,ix) .EQ. 0. ) flagst(isea) = .NOT. flghg1

      END DO

      !

      ! 2.g Test output

      !

#ifdef W3_T

      np     = 1 + (nx-1)/65

      DO ip=1, np

        ixl    = 1 + (ip-1)*65

        ixh    = min( nx, ip*65 )

        WRITE (mdst,9024) ixl, ixh

        DO iy=ny,1 , -1

          WRITE (mdst,9025) nint(mapodi(iy,ixl:ixh)*sig(1)/dtmax)

        END DO

      END DO

#endif

      !

      ! ... End of loop in 1.a

      !

    END DO

    !

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

    ! 3.  Final data base (full data base, scratched at end of routine)

    ! 3.a Loop over grids

    !

    ALLOCATE ( nrec(nrgrd), nsnd(nrgrd), ntpp(nmproc), stat=istat )

    check_alloc_status( istat )

    !

    DO i=1, nrgrd

      IF ( .NOT. dogrid(i) ) cycle

#ifdef W3_T

      WRITE (mdst,9030) i

#endif

      !

      CALL w3setg ( i, mdse, mdst )

      CALL w3seto ( i, mdse, mdst )

      CALL wmsetm ( i, mdse, mdst )

      !

      ALLOCATE ( map3d(ny,nx,-4:nrgrd), wgt3d(ny,nx,0:nrgrd), stat=istat )

      check_alloc_status( istat )

      map3d  = 0

      wgt3d  = 0.

      nrec   = 0

      nsnd   = 0

      !

      ! 3.b Filling MAP3D and WGT3D, as well as NREC and NSND

      !

      DO j=1, nrgrd

        IF ( .NOT. shrank(i,j) ) cycle

#ifdef W3_T

        WRITE (mdst,9031) j

#endif

        mapodi => mdatas(j)%MAPODI

        !

        DO jj=1, stores(i,j)%NTOT

          ix     = stores(i,j)%IX(jj)

          iy     = stores(i,j)%IY(jj)

          wgt3d(iy,ix,0) = mdatas(i)%MAPODI(iy,ix)

          map3d(iy,ix,-2) = mapfs(iy,ix)

          IF ( map3d(iy,ix,-2) .NE. 0 ) THEN

            map3d(iy,ix,-3) = 1 + (map3d(iy,ix,-2)-1)/naproc

#ifdef W3_SHRD

            map3d(iy,ix,-4) = 1

#endif

#ifdef W3_MPI

            map3d(iy,ix,-4) = map3d(iy,ix,-2) -              &

                 (map3d(iy,ix,-3)-1)*naproc + croot - 1

#endif

          END IF

          IF ( wgt3d(iy,ix,0).GE.0. .AND. mapsta(iy,ix).NE.0. .AND. &

               stores(i,j)%NAV(jj).GT.0 ) THEN

            wgt3d(iy,ix,j) = odimax / sig(1) * dtmax

            DO na=1, stores(i,j)%NAV(jj)

              jx     = grids(j)%MAPSF(stores(i,j)%ISS(jj,na),1)

              jy     = grids(j)%MAPSF(stores(i,j)%ISS(jj,na),2)

              IF ( mapodi(jy,jx) .GE. 0. ) wgt3d(iy,ix,j) =       &

                   min( wgt3d(iy,ix,j) , mapodi(jy,jx) )

            END DO

            IF ( wgt3d(iy,ix,j) .GT. 0. ) map3d(iy,ix,j) = 1

          END IF

        END DO

        !

        stores(i,j)%NFIN = sum(map3d(:,:,j))

#ifdef W3_T

        WRITE (mdst,9032) stores(i,j)%NFIN, stores(i,j)%NTOT

#endif

        !

      END DO

      !

      mapodi => mdatas(i)%MAPODI

      DO ix=1, nx

        DO iy=1, ny

          map3d(iy,ix, 0) = maxval(map3d(iy,ix,1:))

          map3d(iy,ix,-1) = sum(map3d(iy,ix,1:))

          IF ( map3d(iy,ix,-1) .GT. 0 ) THEN

            IF ( mapodi(iy,ix)*sig(1)/dtmax .GT. 1.5*odimax ) THEN

              wgt3d(iy,ix, 0:) = 0.

              map3d(iy,ix,-1:) = 0

            ELSE

              wgth   = sum(wgt3d(iy,ix,:))

              IF ( wgth .GT. 1.e-25 ) THEN

                wgt3d(iy,ix,:) = wgt3d(iy,ix,:) / wgth

              ELSE

                wgt3d(iy,ix,:) = 0.

              END IF

              IF ( map3d(iy,ix,-4) .EQ. improc ) THEN

                nrec(i) = nrec(i) + 1

                DO jj=1, nrgrd

                  IF ( map3d(iy,ix,jj) .GT. 0 )               &

                       nrec(jj) = nrec(jj) + 1

                END DO

              END IF

            END IF

          END IF

        END DO

      END DO

      !

      DO j=1, nrgrd

        IF ( .NOT. shrank(i,j) ) cycle

        DO jj=1, stores(i,j)%NTOT

          ix     = stores(i,j)%IX(jj)

          iy     = stores(i,j)%IY(jj)

          IF ( map3d(iy,ix,j) .NE. 0 ) THEN

            DO na=1, stores(i,j)%NAV(jj)

              IF ( stores(i,j)%IPS(jj,na) .EQ. improc )           &

                   nsnd(j) = nsnd(j) + 1

            END DO

          END IF

        END DO

      END DO

      !

      ng     = maxval(map3d(:,:,-1))

      ntl    = sum(map3d(:,:,0))

      !

      ! 3.c Check for points with all ODI = 0

      !

      mapodi => mdatas(i)%MAPODI

      nout   = 0

      !

      jxl    = nx

      jxh    =  1

      jyl    = ny

      jyh    =  1

      !

      ALLOCATE ( mapout(ny,nx), stat=istat )

      check_alloc_status( istat )

      mapout = mapsta

      !

      DO ix=1, nx

        DO iy=1, ny

          IF ( abs(mapsta(iy,ix)).EQ. 1 .AND.                       &

               mapodi(iy,ix) .EQ. 0.    .AND.                       &

               map3d(iy,ix,-1) .EQ. 0 ) THEN

            nout = nout + 1

            IF ( improc.EQ.nmperr .AND. nout.EQ. 1 )              &

                 WRITE(mdse,*) ' '

            IF ( improc.EQ.nmperr .AND. nout.LE.25 )              &

                 WRITE(mdse,1001) i, ix, iy

            IF ( improc.EQ.nmperr .AND. nout.EQ.25 )              &

                 WRITE(mdse,1006)

            jxl    = min( ix, jxl )

            jxh    = max( ix, jxh )

            jyl    = min( iy, jyl )

            jyh    = max( iy, jyh )

            mapout(iy,ix) = 999

          END IF

        END DO

      END DO

      !

      ! 3.d Test and error output

      !

#ifdef W3_T

      WRITE (mdst,9033) ntl, ng, nout

      WRITE (mdst,9034) nrec

      WRITE (mdst,9035) nsnd

      WRITE (mdst,9036)

      DO iy=ny,1 , -1

        WRITE (mdst,9037) map3d(iy,:,-1)

      END DO

#endif

      !

      IF ( nout .GT. 0 ) THEN

        IF ( improc.EQ.nmperr ) THEN

          WRITE(mdse,1000) i, nout

          extra  = 2

          jxl = max(  1, jxl - extra )

          jxh = min( nx, jxh + extra )

          jyl = max(  1, jyl - extra )

          jyh = min( ny, jyh + extra )

          WRITE (mdse,1002) jxl, jxh, jyl, jyh

          np     = 1 + (jxh-jxl)/65

          DO ip=1, np

            ixl    = jxl + (ip-1)*65

            ixh    = min( nx, ixl+64 )

            WRITE (mdse,1005) ixl, ixh

            WRITE (mdse,1003) 'STATUS MAP MAPSTA'

            DO iy=jyh, jyl, -1

              WRITE (mdse,1004) mapsta(iy,ixl:ixh)

            END DO

            WRITE (mdse,1003) 'MISSING POINTS IN MAPSTA (**)'

            DO iy=jyh, jyl, -1

              WRITE (mdse,1004) mapout(iy,ixl:ixh)

            END DO

            WRITE (mdse,1003) 'OPEN BOUND. DISTANCE MAP MAPODI'

            DO iy=jyh, jyl, -1

              WRITE (mdse,1004)                                 &

                   nint(mapodi(iy,ixl:ixh)*sig(1)/dtmax)

            END DO

            WRITE (mdse,1003) 'GRID COVERAGE MAP MAP3D'

            DO iy=jyh, jyl, -1

              WRITE (mdse,1004) map3d(iy,ixl:ixh,-1)

            END DO

            WRITE (mdse,*)

          END DO

        END IF

        CALL extcde (1000)

      END IF

      !

      DEALLOCATE ( mapout, stat=istat )

      check_dealloc_status( istat )

      !

#ifdef W3_T7

      WRITE (mdst,9330) i

      DO j=1, nrgrd

        IF ( .NOT. shrank(i,j) ) THEN

          IF ( i .NE. j ) WRITE (mdst,9331) j

          cycle

        END IF

        WRITE (mdst,9332) j, stores(i,j)%NFIN, i, j

        IF ( stores(i,j)%NFIN .EQ. 0 ) cycle

        ntl    = 0

        DO jj=1, stores(i,j)%NTOT

          ix     = stores(i,j)%IX(jj)

          iy     = stores(i,j)%IY(jj)

          IF ( map3d(iy,ix,j) .EQ. 0 ) cycle

          ntl    = ntl + 1

          na     = stores(i,j)%NAV(jj)

          WRITE (mdst,9333) ntl, ix, iy, map3d(iy,ix,-2),       &

               map3d(iy,ix,-3), map3d(iy,ix,-4),   &

               wgt3d(iy,ix,0), wgt3d(iy,ix,j), na, &

               stores(i,j)%ISS(jj,1),              &

               stores(i,j)%JSS(jj,1),              &

               stores(i,j)%IPS(jj,1),              &

               stores(i,j)%AWG(jj,1),              &

               stores(i,j)%ITG(jj,1)

          DO ia=2, na

            WRITE (mdst,9334) stores(i,j)%ISS(jj,ia),           &

                 stores(i,j)%JSS(jj,ia),           &

                 stores(i,j)%IPS(jj,ia),           &

                 stores(i,j)%AWG(jj,ia),           &

                 stores(i,j)%ITG(jj,ia)

          END DO

        END DO

      END DO

#endif

      !

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

      ! 4.  Save data base as needed in EQSTGE

      !

      ! 4.a   ALLOCATE storage

      ! 4.a.1 Local counters, weights and sea counters (grid 'I')

      !

      IF ( eqstge(i,i)%NREC .NE. 0 ) THEN

        DEALLOCATE (eqstge(i,i)%ISEA , eqstge(i,i)%JSEA ,         &

             eqstge(i,i)%WGHT, stat=istat )

        check_dealloc_status( istat )

        eqstge(i,i)%NREC = 0

#ifdef W3_T

        WRITE (mdst,9040) i, i

#endif

      END IF

      !

      IF ( nrec(i) .GT. 0 ) THEN

        ALLOCATE ( eqstge(i,i)%ISEA(nrec(i))  ,                   &

             eqstge(i,i)%JSEA(nrec(i))  ,                   &

             eqstge(i,i)%WGHT(nrec(i)), stat=istat )

        check_alloc_status( istat )

        eqstge(i,i)%NREC = nrec(i)

#ifdef W3_T

        WRITE (mdst,9041) i, i, nrec(i)

#endif

      END IF

      !

      ! 4.a.1 Local counters, arrays weights etc.  (grid 'J' receive)

      !

      DO j=1, nrgrd

        IF ( i .EQ. j ) cycle

        eqstge(i,i)%NTOT = stores(i,j)%NFIN

        !

        IF ( eqstge(i,j)%NREC .NE. 0 ) THEN

          DEALLOCATE ( eqstge(i,j)%ISEA , eqstge(i,j)%JSEA ,      &

               eqstge(i,j)%WGHT , eqstge(i,j)%SEQL ,      &

               eqstge(i,j)%NAVG , eqstge(i,j)%WAVG ,      &

               eqstge(i,j)%RIP  , eqstge(i,j)%RTG, stat=istat )

          check_dealloc_status( istat )

          eqstge(i,j)%NREC   = 0

          eqstge(i,j)%NAVMAX = 1

#ifdef W3_T

          WRITE (mdst,9042) i, j

#endif

        END IF

        !

        IF ( nrec(j) .GT. 0 ) THEN

          na     =  maxval( stores(i,j)%NAV(1:stores(i,j)%NTOT) )

          eqstge(i,j)%NAVMAX = na

          ALLOCATE ( eqstge(i,j)%ISEA(nrec(j))  ,                 &

               eqstge(i,j)%JSEA(nrec(j))  ,                 &

               eqstge(i,j)%WGHT(nrec(j))  ,                 &

               eqstge(i,j)%SEQL(sgrds(j)%NSPEC,nrec(j),na), &

               eqstge(i,j)%NAVG(nrec(j))  ,                 &

               eqstge(i,j)%WAVG(nrec(j),na),                &

               eqstge(i,j)%RIP(nrec(j),na),                 &

               eqstge(i,j)%RTG(nrec(j),na), stat=istat )

          check_alloc_status( istat )

          eqstge(i,j)%NREC = nrec(j)

#ifdef W3_T

          WRITE (mdst,9043) i, j, nrec(j), na

#endif

        END IF

        !

        IF ( eqstge(i,j)%NSND .NE. 0 ) THEN

          DEALLOCATE ( eqstge(i,j)%SIS , eqstge(i,j)%SJS ,        &

               eqstge(i,j)%SI1 , eqstge(i,j)%SI2 ,        &

               eqstge(i,j)%SIP , eqstge(i,j)%STG, stat=istat )

          check_dealloc_status( istat )

          eqstge(i,j)%NSND = 0

#ifdef W3_T

          WRITE (mdst,9044) i, j

#endif

        END IF

        !

        IF ( nsnd(j) .GT. 0 ) THEN

          ALLOCATE ( eqstge(i,j)%SIS(nsnd(j)) ,                   &

               eqstge(i,j)%SJS(nsnd(j)) ,                   &

               eqstge(i,j)%SI1(nsnd(j)) ,                   &

               eqstge(i,j)%SI2(nsnd(j)) ,                   &

               eqstge(i,j)%SIP(nsnd(j)) ,                   &

               eqstge(i,j)%STG(nsnd(j)), stat=istat )

          check_alloc_status( istat )

          eqstge(i,j)%NSND = nsnd(j)

#ifdef W3_T

          WRITE (mdst,9045) i, j, nsnd(j)

#endif

        END IF

        !

      END DO

      !

      ! 4.b   Store data in EQSTGE

      ! 4.b.1 Grid I (JSEA and weight only)

      !

      IF ( eqstge(i,i)%NREC .GT. 0 ) THEN

        ntl    = 0

        DO ix=1, nx

          DO iy=1, ny

            IF ( map3d(iy,ix,0) .EQ. 0 ) cycle

            IF ( map3d(iy,ix,-4) .NE. improc ) cycle

            ntl    = ntl + 1

            eqstge(i,i)%ISEA(ntl) = map3d(iy,ix,-2)

            eqstge(i,i)%JSEA(ntl) = map3d(iy,ix,-3)

            eqstge(i,i)%WGHT(ntl) = wgt3d(iy,ix,0)

          END DO

        END DO

      END IF

      !

      ! 4.b.2 All other grids, info for receiving

      !

      DO j=1, nrgrd

        IF ( .NOT. shrank(i,j) ) cycle

        IF ( eqstge(i,j)%NREC .EQ. 0 ) cycle

        ntl    = 0

        !

        DO jj=1, stores(i,j)%NTOT

          ix     = stores(i,j)%IX(jj)

          iy     = stores(i,j)%IY(jj)

          IF ( map3d(iy,ix,j) .EQ. 0 ) cycle

          IF ( map3d(iy,ix,-4) .NE. improc ) cycle

          ntl    = ntl + 1

          eqstge(i,j)%ISEA(ntl) = map3d(iy,ix,-2)

          eqstge(i,j)%JSEA(ntl) = map3d(iy,ix,-3)

          eqstge(i,j)%WGHT(ntl) = wgt3d(iy,ix,j)

          na     = stores(i,j)%NAV(jj)

          eqstge(i,j)%NAVG(ntl) = na

          eqstge(i,j)%WAVG(ntl,1:na) = stores(i,j)%AWG(jj,1:na)

          eqstge(i,j)%RIP (ntl,1:na) = stores(i,j)%IPS(jj,1:na)

          eqstge(i,j)%RTG (ntl,1:na) = stores(i,j)%ITG(jj,1:na)

        END DO

        !

      END DO

      !

      ! 4.b.3 All other grids, info for sending

      !

      DO j=1, nrgrd

        IF ( .NOT. shrank(i,j) ) cycle

        IF ( eqstge(i,j)%NSND .EQ. 0 ) cycle

        ntpp   = 0

        ntl    = 0

        !

        DO jj=1, stores(i,j)%NTOT

          ix     = stores(i,j)%IX(jj)

          iy     = stores(i,j)%IY(jj)

          IF ( map3d(iy,ix,j) .NE. 0 ) THEN

            ntpp(map3d(iy,ix,-4)) = ntpp(map3d(iy,ix,-4)) + 1

            DO na=1, stores(i,j)%NAV(jj)

              IF ( stores(i,j)%IPS(jj,na) .EQ. improc ) THEN

                ntl    = ntl + 1

                eqstge(i,j)%SIS(ntl) = stores(i,j)%ISS(jj,na)

                eqstge(i,j)%SJS(ntl) = stores(i,j)%JSS(jj,na)

                eqstge(i,j)%SI1(ntl) = ntpp(map3d(iy,ix,-4))

                eqstge(i,j)%SI2(ntl) = na

                eqstge(i,j)%SIP(ntl) = map3d(iy,ix,-4)

                eqstge(i,j)%STG(ntl) = stores(i,j)%ITG(jj,na)

              END IF

            END DO

          END IF

        END DO

        !

      END DO

      !

      ! 4.c Detailed test output

      !

#ifdef W3_T5

      dstr   = '                  '

#endif

      !

#ifdef W3_T5

      IF ( eqstge(i,i)%NREC .EQ. 0 ) THEN

        WRITE (mdst,9140) i

      ELSE

        WRITE (mdst,9141) i

        na     = 0

        DO j=1, nrgrd

          IF ( i.EQ.j .OR. eqstge(i,j)%NREC.EQ.0 ) cycle

          na     = na + 1

          nsnd(na) = j

        END DO

        WRITE (mdst,9142) nsnd(1:na)

        WRITE (mdst,9143)

        ALLOCATE ( tstr(na), stat=istat )

        check_alloc_status( istat )

        DO jj=1, eqstge(i,i)%NREC

          DO ng=1, na

            j      = nsnd(ng)

            tstr(ng) = dstr

            DO ntl=1, eqstge(i,j)%NREC

              IF ( eqstge(i,i)%ISEA(jj) .EQ.                &

                   eqstge(i,j)%ISEA(ntl) ) THEN

                WRITE (tstr(ng),9144) ntl,                &

                     eqstge(i,j)%WGHT(ntl),        &

                     eqstge(i,j)%NAVG(ntl)

                EXIT

              END IF

            END DO

          END DO

          WRITE (mdst,9145) jj, eqstge(i,i)%ISEA(jj),         &

               eqstge(i,i)%JSEA(jj),         &

               eqstge(i,i)%WGHT(jj),         &

               tstr

        END DO

        DEALLOCATE ( tstr, stat=istat )

        check_dealloc_status( istat )

      END IF

#endif

      !

#ifdef W3_T5

      DO j=1, nrgrd

        IF ( i.EQ.j .OR. eqstge(i,j)%NREC.EQ.0 ) cycle

        WRITE (mdst,9146) j

        DO jj=1, eqstge(i,j)%NREC

          WRITE (mdst,9147) jj, eqstge(i,j)%NAVG(jj),           &

               ( eqstge(i,j)%WAVG(jj,na),            &

               eqstge(i,j)%RIP (jj,na),            &

               eqstge(i,j)%RTG (jj,na),            &

               na=1, eqstge(i,j)%NAVG(jj) )

        END DO

      END DO

#endif

      !

#ifdef W3_T6

      DO j=1, nrgrd

        IF ( i .EQ. j ) cycle

        IF ( eqstge(i,j)%NSND .EQ. 0 ) THEN

          WRITE (mdst,9240) j

        ELSE

          WRITE (mdst,9241) j

          DO jj=1, eqstge(i,j)%NSND

            WRITE (mdst,9242) jj, eqstge(i,j)%SIS(jj),        &

                 eqstge(i,j)%SJS(jj),        &

                 eqstge(i,j)%SI1(jj),        &

                 eqstge(i,j)%SI2(jj),        &

                 eqstge(i,j)%SIP(jj),        &

                 eqstge(i,j)%STG(jj)

          END DO

        END IF

      END DO

#endif

      !

      ! ... End of loop started in 3.a

      !

      DEALLOCATE ( map3d, wgt3d, stat=istat )

      check_dealloc_status( istat )

    END DO

    !

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

    ! 5.  Generate GRDEQL

    ! 5.a Size of array

    !

    nrec   = 0

    !

    DO i=1, nrgrd

      DO j=1, nrgrd

        IF ( i.EQ.j .OR. stores(i,j)%NFIN.EQ.0 ) cycle

        nrec(i) = nrec(i) + 1

      END DO

    END DO

    !

    na     = maxval(nrec)

    ALLOCATE ( grdeql(nrgrd,0:na), stat=istat )

    check_alloc_status( istat )

    grdeql = 0

    !

#ifdef W3_T

    WRITE (mdst,9050) na

#endif

    !

    ! 5.b Fill array

    !

    DO i=1, nrgrd

      grdeql(i,0) = nrec(i)

      jj          = 0

      DO j=1, nrgrd

        IF ( i.EQ.j .OR. stores(i,j)%NFIN.EQ.0 ) cycle

        jj           = jj + 1

        grdeql(i,jj) = j

      END DO

    END DO

    !

#ifdef W3_T

    WRITE (mdst,9051)

    DO i=1, nrgrd

      WRITE (mdst,9052) i, grdeql(i,0:grdeql(i,0))

    END DO

#endif

    !

    ! 5.c Resolution test

    !


    IF ( flagll ) THEN

      factor = 1.

    ELSE

      factor = 1.e-3

    END IF

    !

    ! notes: This resolution test, with FACMAX=2, is pretty strict, so

    !        it is not going to be appropriate for irregular grids.

    !        We'll just have to trust the judgement of the user in the

    !        case of irregular grids. But if we change our minds and do

    !        some kind of check for irregular grids, we could make

    !        a check against median(HPFAC) and median(HQFAC).


    DO i=1, nrgrd

      CALL w3setg ( i, mdse, mdst )

      IF ( gtype.EQ.rlgtype ) THEN

        DO jj=1, grdeql(i,0)

          j      = grdeql(i,jj)

          IF ( grids(j)%GTYPE.EQ.rlgtype ) THEN

            IF ( sx/grids(j)%SX .GT. facmax     .OR.                   &

                 sx/grids(j)%SX .LT. 1./facmax  .OR.                   &

                 sy/grids(j)%SY .GT. facmax     .OR.                   &

                 sy/grids(j)%SY .LT. 1./facmax ) THEN

              IF ( improc.EQ.nmperr ) WRITE(mdse,1050) i, factor*sx,  &

                   factor*sy, j, factor*grids(j)%SX, factor*grids(j)%SY

              CALL extcde ( 1050 )

            END IF ! IF ( SX/GR ...

          END IF ! IF ( GRIDS(J)%GTYPE...

        END DO ! DO JJ=...

      END IF ! IF ( GTYPE....

    END DO ! DO I=...

    !

    ! 5.d Group number test

    !

    DO i=1, nrgrd

      IF ( grdeql(i,0) .GE. 2 ) THEN

        tgrp   = grgrp(grdeql(i,1))

        DO j=2, grdeql(i,0)

          IF ( grgrp(grdeql(i,j)) .NE. tgrp ) THEN

            IF ( improc .EQ. nmperr ) WRITE(mdse,1051)          &

                 grdeql(i,1), grgrp(grdeql(i,1)),   &

                 grdeql(i,j), grgrp(grdeql(i,j))

            CALL extcde ( 1051 )

          END IF

        END DO

      END IF

    END DO

    !

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

    ! 6.  Final clean up

    !

    DO i=1, nrgrd

      DO j=1, nrgrd

        IF ( stores(i,j)%INIT ) THEN

          DEALLOCATE ( stores(i,j)%IX  , stores(i,j)%IY  ,        &

               stores(i,j)%NAV , stores(i,j)%FLA ,        &

               stores(i,j)%ISS , stores(i,j)%JSS ,        &

               stores(i,j)%IPS , stores(i,j)%ITG ,        &

               stores(i,j)%AWG , stat=istat )

          check_dealloc_status( istat )

        END IF

      END DO

    END DO

    !

    DEALLOCATE ( shrank, stores, nrec, nsnd, ntpp, stat=istat )

    check_dealloc_status( istat )

    !

    RETURN

    !

    ! Formats

    !

1000 FORMAT (/' *** WAVEWATCH III ERROR IN WMGEQL : *** '/           &

         '     UNCOVERED EDGE POINTS FOR GRID',i4,'  (',i6,')'/)

1001 FORMAT ( '     GRID',i4,'  POINT',2i5,' NOT COVERED (WMGEQL)')

1002 FORMAT ( '     DIAGNOSTICS IX AND IY RANGE:',4i6)

1003 FORMAT (/'     SHOWING ',a/)

1004 FORMAT (2x,65i2)

1005 FORMAT (/'     SHOWING IX RANGE ',2i6)

1006 FORMAT ( '        (WILL NOT PRINT ANY MORE UNCOVERED POINTS)')

    !

1020 FORMAT (/' *** WAVEWATCH III WARNING WMGEQL : *** '/            &

         '     REMOVED BOUNDARY POINT FROM OPEN EDGE DISTANCE MAP'/     &

         '     GRID, IX, IY :',3i6)

    !

1050 FORMAT (/' *** WAVEWATCH III ERROR IN WMGEQL : *** '/           &

         '     GRID INCREMENTS TOO DIFFERENT '/                 &

         '     GRID',i4,'   INCREMENTS ',2f8.2/                 &

         '     GRID',i4,'   INCREMENTS ',2f8.2/)

1051 FORMAT (/' *** WAVEWATCH III ERROR IN WMGEQL : *** '/           &

         '     OVERLAPPING GRIDS NEED TO BE IN SAME GROUP '/    &

         '         GRID',i4,' IN GROUP',i4/                     &

         '         GRID',i4,' IN GROUP',i4/)

    !

#ifdef W3_T

9010 FORMAT ( ' TEST WMGEQL : STARTING LOOP OVER GRIDS')

9011 FORMAT ( ' TEST WMGEQL : I, RANK :',2i4)

9012 FORMAT ( '               GRID ',i3,' HAS SAME RANK')

9013 FORMAT ( '               ',a)

#endif

    !

#ifdef W3_T

9020 FORMAT ( ' TEST WMGEQL : GENERATING DISTANCE MAP GRID ',i3)

9024 FORMAT ( ' TEST WMGEQL : FINAL MAP FOR X RANGE ',2i6)

9025 FORMAT (2x,65i2)

#endif

    !

#ifdef W3_T

9030 FORMAT ( ' TEST WMGEQL : DEPENDENCE INFORMATION GRID ',i3)

9031 FORMAT ( '               CHECKING GRID ',i3)

9032 FORMAT ( '               POINTS USED/AVAIL :',2i6)

9033 FORMAT ( '               TOTAL/GRIDS/OUT   :',3i6)

9034 FORMAT ( '               LOCAL PER GRID    :',15i6)

9035 FORMAT ( '               SENDING PER GRID  :',15i6)

9036 FORMAT ( ' TEST WMGEQL : NUMBER OF CONTRIBUTING GRIDS MAP')

9037 FORMAT (2x,65i2)

#endif

    !

#ifdef W3_T

9040 FORMAT ( ' TEST WMGEQL : GRID ',i2,'-',i2,' CLEAR STORAGE')

9041 FORMAT ( ' TEST WMGEQL : GRID ',i2,'-',i2,' STORAGE SIZE',i6)

9042 FORMAT ( '               RECV ',i2,'-',i2,' CLEAR STORAGE')

9043 FORMAT ( '               RECV ',i2,'-',i2,' STORAGE SIZE',2i6)

9044 FORMAT ( '               SEND ',i2,'-',i2,' CLEAR STORAGE')

9045 FORMAT ( '               SEND ',i2,'-',i2,' STORAGE SIZE',i6)

#endif

    !

#ifdef W3_T

9050 FORMAT ( ' TEST WMGEQL : GRDEQL DIMENSIONED AT ',i2)

9051 FORMAT ( ' TEST WMGEQL : GRDEQL :')

9052 FORMAT ( '                 ',2i4,' : ',20i3)

#endif

    !

#ifdef W3_T5

9140 FORMAT ( ' TEST WMGEQL : NO RECEIVING DATA FOR GRID ',i3,   &

         ' <=====================================')

9141 FORMAT ( ' TEST WMGEQL : RECEIVING DATA GRID ',i3,          &

         ' <=====================================')

9142 FORMAT ( '               RECEIVING FROM GRID(S) ',10i3)

9143 FORMAT (16x,'COUNT, ISEA, JSEA, WEIGHT / ',                 &

         'COUNT WEIGHT NR PER GRID')

9144 FORMAT (i6,f6.2,i6)

9145 FORMAT (12x,3i6,f6.2,10(' - ',a))

9146 FORMAT ( ' TEST WMGEQL : RECEIVING DATA AVG. GRID ',i3)

9147 FORMAT (12x,i6,i2,4(f8.2,i4,i6))

#endif

    !

#ifdef W3_T6

9240 FORMAT ( ' TEST WMGEQL : NO SENDING DATA FOR GRID ',i3,     &

         ' <=====================================')

9241 FORMAT ( ' TEST WMGEQL : SENDING DATA GRID ',i3,            &

         ' <====================================='/         &

         11x,'COUNT, ISEA, JSEA, ARRAY IND., PROC, TAG')

9242 FORMAT ( '              ',4i8,i4,2i8)

#endif

    !

#ifdef W3_T7

9330 FORMAT ( ' TEST WMGEQL : FULL SOURCE INFO GRID ',i3,        &

         ' <=====================================')

9331 FORMAT ( '               GRID ',i3,' IS NOT OF SAME RANK')

9332 FORMAT ( '               GRID ',i3,' CONTRIBUTES TO',i6,    &

         ' GRID POINTS'/                                    &

         18x,'<---------- GRID',i6,' ---------->',          &

         4x,'<----------- GRID',i6,' ----------->'/        &

         18x,'NR   IX   IY  ISEA  JSEA  IP  WGTH',          &

         2x,'  WGTH  NA  ISEA  JSEA  IP  WGTH   TAG' )

9333 FORMAT (15x,3i5,2i6,i4,f6.2,2x,f6.2,i4,2i6,i4,f6.2,i6)

9334 FORMAT (64x,2i6,i4,f6.2,i6)

#endif

    !/

    !/ End of WMGEQL ----------------------------------------------------- /

    !/

  END SUBROUTINE wmgeql

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

  SUBROUTINE wmrspc

    !/

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

    !/                  | WAVEWATCH III           NOAA/NCEP |

    !/                  |           H. L. Tolman            |

    !/                  |                        FORTRAN 90 |

    !/                  | Last update :         10-Dec-2014 !

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

    !/

    !/    22-Sep-2005 : Origination.                        ( version 3.08 )

    !/    25-Jul-2006 : Point output grid added.            ( version 3.10 )

    !/    30-Oct-2009 : Implement run-time grid selection.  ( version 3.14 )

    !/                  (W. E. Rogers & T. J. Campbell, NRL)

    !/    10-Dec-2014 : Add checks for allocate status      ( version 5.04 )

    !/

    !  1. Purpose :

    !

    !     Generate map with flogs for need of spectral grid conversion

    !     between models.

    !

    !  2. Method :

    !

    !     Test of parameters as introduced before in W3IOBC.

    !

    !  3. Parameters :

    !

    !  4. Subroutines used :

    !

    !      Name      Type  Module   Description

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

    !      STRACE    Sur.  W3SERVMD Subroutine tracing.

    !      EXTCDE    Subr.   Id.    Program abort.

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

    !

    !  5. Called by :

    !

    !      Name      Type  Module   Description

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

    !      WMINIT    Subr  WMINITMD Multi-grid model initialization.

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

    !

    !  6. Error messages :

    !

    !  7. Remarks :

    !

    !  8. Structure :

    !

    !     See source code.

    !

    !  9. Switches :

    !

    !     !/S    Enable subroutine tracing.

    !     !/T    Enable test output

    !

    ! 10. Source code :

    !

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

    !

    USE w3servmd, ONLY: extcde

#ifdef W3_S

    USE w3servmd, ONLY: strace

#endif

    !

    USE w3gdatmd

    USE w3odatmd, ONLY: unipts

    USE wmmdatmd

    !

    IMPLICIT NONE

    !/

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

    !/ Parameter list

    !/

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

    !/ Local parameters

    !/

    INTEGER                 :: I, J, LOW

#ifdef W3_S

    INTEGER, SAVE           :: IENT = 0

#endif

    !/

#ifdef W3_S

    CALL strace (ient, 'WMRSPC')

#endif

    !

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

    ! 0.  Initializations

    !

    IF ( unipts ) THEN

      low    = 0

    ELSE

      low    = 1

    END IF

    IF ( .NOT. ALLOCATED(respec) ) THEN

      ALLOCATE ( respec(low:nrgrd,low:nrgrd), stat=istat )

      check_alloc_status( istat )

    END IF

    respec = .false.

    !

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

    ! 1.  Fill map with flags

    !

    DO i=low, nrgrd

      DO j=i+1, nrgrd

        respec(i,j) = sgrds(i)%NK    .NE. sgrds(j)%NK     .OR.      &

             sgrds(i)%NTH   .NE. sgrds(j)%NTH    .OR.      &

             sgrds(i)%XFR   .NE. sgrds(j)%XFR    .OR.      &

             sgrds(i)%FR1   .NE. sgrds(j)%FR1    .OR.      &

             sgrds(i)%TH(1) .NE. sgrds(j)%TH(1)

        respec(j,i) = respec(i,j)

      END DO

    END DO

    !

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

    ! 2.  Test output

    !

#ifdef W3_T

    WRITE (mdst,9000)

    DO i=low, nrgrd

      WRITE (mdst,9001) i, respec(i,:)

    END DO

#endif

    !

    RETURN

    !

    ! Formats

    !

#ifdef W3_T

9000 FORMAT ( 'TEST WMRSPC : MAP RESPEC FILLED ')

9001 FORMAT ( '              ',i4,' : ',20l2)

#endif

    !/

    !/ End of WMRSPC ----------------------------------------------------- /

    !/

  END SUBROUTINE wmrspc

  !/

  !!

  SUBROUTINE wmsmceql

    !!

    !!    Adapted from multi-grid sub WMGEQL for set up equal ranked SMC

    !!    grid boundary points.       JGLi10Aug2020

    !!    Move boundary point matching to sub-grid root PEs and broadcast to

    !!    all other PEs.              JGLi02Dec2020

    !!    Clear bugs for 3 sub-grids case and finalise output messages.

    !!                                JGLi26Jan2021

    !!    Add regular grid to SMC grid same ranked group.

    !!                                JGLi12Apr2021

    !!

    !  1. Purpose :

    !

    !     Determine relations to same ranked SMC grids for each grid.

    !     Set boundary points update for regular grid in same ranked group.

    !

    !  2. Method :

    !

    !     Cross mapping of grid points, use nearest sea points and no

    !     interpolation is required so far.

    !

    !  3. Parameters :

    !

    !  4. Subroutines used :

    !

    !      Name      Type  Module   Description

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

    !      W3SETG, W3SETO, WMSETM

    !                Subr. W3GDATMD Manage data structures.

    !      W3SMCGMP, Subr. W3PSMCMD Mapping Lon-Lat points to SMC grid cells.

    !      W3SMCELL, Subr. W3PSMCMD Find Lon-Lat for SMC grid cell centre.

    !      STRACE    Subr. W3SERVMD Subroutine tracing.

    !      EXTCDE    Subr.   Id.    Program abort.

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

    !

    !  5. Called by :

    !

    !  6. Error messages :

    !

    !  7. Remarks :

    !

    !       SMC J sub-grid has own boundary cell number W3GDATMD's NBSMC and

    !       their ID list are stored in W3GDATMD's GRIDS(J)%ISMCBP(NBSMC),

    !       which are the global ISEA values of the NBSMC boundary cells.

    !       So there is no need to look for boundary points, but just

    !       fetching the boundary cell list from each SMC sub-grid.

    !       No interpolation is required as one to one correspondance is

    !       assumed among SMC sub-grid boundary points.  JGLi06Nov2020

    !       Sub WMIOEG and WMIOES are modified to use the new EQSTGE array

    !       for same ranked SMC sub-grids only.          JGLi26Jan2021

    !

    !  8. Structure :

    !

    !  9. Switches :

    !

    !     !/PRn  Propagation scheme.

    !     !/SMC  For SMC grid.

    !

    !     !/S    Enable subroutine tracing.

    !     !/T    Enable test output.

    !

    !     !/MPI  Distribbuted memory management.

    !     !/SHRD Shared memory case.

    !

    ! 10. Source code :

    !

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

    !

    USE constants

    USE w3gdatmd

    USE w3odatmd

    USE w3adatmd

    USE wmmdatmd

    !

    USE w3servmd, ONLY: extcde

#ifdef W3_SMC

    USE w3psmcmd, ONLY: w3smcgmp, w3smcell

#endif


#ifdef W3_S

    USE w3servmd, ONLY: strace

#endif

    !

    IMPLICIT NONE

    !

#ifdef W3_MPI

    include "mpif.h"

#endif

    !/

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

    !/ Parameter list

    !/

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

    !/ Local parameters

    !/

    INTEGER                :: I, J, IX, IY, IXY, JX, JY, NPJ,      &

         NR, NT, NA, NTL, JJ, NIT, NG, NOUT,  &

         ISEA, JSEA, IPRC, ITAG, TGRP, NPMX,  &

         IP, NP, ICROOT, JCROOT, IEER


#ifdef W3_MPI

    INTEGER, Dimension(MPI_STATUS_SIZE):: MPIState

#endif


#ifdef W3_S

    INTEGER, SAVE            :: IENT = 0

#endif

    INTEGER, ALLOCATABLE   :: NREC(:),  NSND(:),  NTPP(:),  &

         IBPTS(:), JBPTS(:), IPBPT(:)

    REAL, PARAMETER        :: ODIMAX = 25.

    REAL, ALLOCATABLE      :: XLon(:), YLat(:)

    LOGICAL                :: CHANGE

    LOGICAL, ALLOCATABLE   :: SHRANK(:,:), DOGRID(:)

#ifdef W3_T5

    CHARACTER(LEN=18), ALLOCATABLE :: TSTR(:)

    CHARACTER(LEN=18)              :: DSTR

#endif

    !

    TYPE store

      INTEGER              :: NTOT, NFIN

      INTEGER, POINTER     :: ICVBP(:), MSDBP(:), ISS(:), JSS(:),   &

           JCVBP(:), IPCVB(:), IPS(:), ITG(:)

      LOGICAL, POINTER     :: FLA(:)

      LOGICAL              :: INIT

    END TYPE store

    !

    TYPE(store), ALLOCATABLE :: STORES(:,:)

    !/

#ifdef W3_S

    CALL strace (ient, 'WMSMCEQL ')

#endif

    !

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

    ! 0.  Initializations

    !

    ALLOCATE ( shrank(nrgrd,nrgrd), stores(nrgrd,nrgrd),            &

         dogrid(nrgrd), stat=istat )

    check_alloc_status( istat )

    !

    shrank = .false.

    !

    DO i=1, nrgrd


      DO j=1, nrgrd

        stores(i,j)%INIT = .false.

        stores(i,j)%NTOT = 0

        stores(i,j)%NFIN = 0

      END DO

    END DO

    !

    itag   = 0

    !

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

    ! 1.  Grid point relations and temp data storage

    ! 1.a Outer loop over all grids

    !

#ifdef W3_T

    WRITE (mdst,9010)

#endif

    !

    DO i=1, nrgrd


#ifdef W3_T

      WRITE (mdst,9011) i, grank(i)

#endif

      !

      ! 1.b Find sub grids with same rank

      !

      nr = 0

      !

      DO j=1, nrgrd

        IF( (grank(i).NE.grank(j)) .OR. (i.EQ.j) ) cycle

        shrank(i,j) = .true.

        nr = nr + 1

      END DO

      !

      dogrid(i) = nr .GT. 0

      !

      IF( nr .EQ. 0 ) cycle

      !

      CALL w3setg( i, mdse, mdst )

      !

      ! Find local root PE and NAPROC for I grid.

#ifdef W3_SHRD

      icroot = 1

#endif

#ifdef W3_MPI

      icroot = mdatas(i)%CROOT

#endif

      np     = outpts(i)%NAPROC

      !

      ! 1.c Fetch Grid I boundary points.

      !

      nt = 0

      IF( grids(i)%GTYPE .EQ. rlgtype ) THEN

        ! 1.c.1 Regular grid I boundary points are stored in NBI.

        nt = outpts(i)%OUT5%NBI

#ifdef W3_MPI

        IF( improc .EQ. icroot )  THEN

#endif

          WRITE(mdse,*) "ICROOT, NT are", icroot, nt

#ifdef W3_MPI

        ENDIF

#endif

        !

        ! 1.c.2 SMC grid I boundary cell ids are saved in NBSMC.

#ifdef W3_SMC

      ELSEIF( grids(i)%GTYPE .EQ. smctype ) THEN

#endif

#ifdef W3_MPI

#ifdef W3_SMC

        IF( improc .EQ. icroot )  THEN

#endif

#endif

#ifdef W3_SMC

          nt = grids(i)%NBSMC

          WRITE(mdse,*) "ICROOT, NT are", icroot, nt

#endif

#ifdef W3_MPI

#ifdef W3_SMC

        ENDIF

#endif

#endif


#ifdef W3_MPI

#ifdef W3_SMC

        CALL mpi_bcast( nt, 1, mpi_integer, &

             icroot-1, mpi_comm_mwave, ieer)

#endif

#endif


        !  Need to wait for all PEs get these values.

#ifdef W3_MPI

#ifdef W3_SMC

        CALL mpi_barrier (mpi_comm_mwave,ieer)

#endif

#endif

        !

      ENDIF  !! GTYPE .EQ. RLGTYPE

      !

      IF( nt .EQ. 0 ) cycle


      IF( nt > 0 ) THEN

        ALLOCATE( ibpts(nt), jbpts(nt), ipbpt(nt),   &

             xlon(nt), ylat(nt), stat=istat )

        check_alloc_status( istat )


        !  Use saved I-grid boundary cell list.

#ifdef W3_MPI

        IF( improc .EQ. icroot )  THEN

#endif

          IF( grids(i)%GTYPE .EQ. rlgtype ) THEN

            !! Loop over regular grid mesh to find boundary points.

            ixy = 0

            DO isea=1, nsea

              ix = mapsf(isea, 1)

              iy = mapsf(isea, 2)

              IF( abs(mapsta(iy,ix)) .EQ. 2 ) THEN

                ixy = ixy + 1

                xlon(ixy) = real(xgrd(iy,ix))

                ylat(ixy) = real(ygrd(iy,ix))

                ibpts(ixy) = isea

                jbpts(ixy) = 1 + (isea - 1)/np

                ipbpt(ixy) = icroot-1 + isea-(jbpts(ixy)-1)*np

              ENDIF

            ENDDO

            !

#ifdef W3_SMC

          ELSEIF( grids(i)%GTYPE .EQ. smctype ) THEN


            ibpts = grids(i)%ISMCBP(1:nt)

            CALL w3smcell( i, nt, ibpts, xlon, ylat )

            DO ix = 1, nt

#endif

              !  Global processor IPBPT and local JSEA, for ISEA spectrum in I grid.

#ifdef W3_SMC

              isea = ibpts(ix)

              jsea = 1 + (isea - 1)/np

              ipbpt(ix) = icroot - 1 + isea - (jsea - 1)*np

              jbpts(ix) = jsea

            ENDDO

#endif

            !

          ENDIF     !!  RLGTYPE

#ifdef W3_MPI

        ENDIF  !!  ICROOT

#endif

        !

        !  All have to wait for ICROOT finishes conversion of cell ids to XLon-YLat

#ifdef W3_MPI

        CALL mpi_barrier (mpi_comm_mwave,ieer)

#endif

        !

        ! Then broadcast IBPTS, IPBPT, XLon, and YLat to all PEs

#ifdef W3_MPI

        CALL mpi_bcast( ibpts(1), nt, mpi_integer,  &

             icroot-1, mpi_comm_mwave, ieer)

        CALL mpi_bcast( jbpts(1), nt, mpi_integer,  &

             icroot-1, mpi_comm_mwave, ieer)

        CALL mpi_bcast( ipbpt(1), nt, mpi_integer,  &

             icroot-1, mpi_comm_mwave, ieer)

        CALL mpi_bcast( xlon(1), nt, mpi_real,  &

             icroot-1, mpi_comm_mwave, ieer)

        CALL mpi_bcast( ylat(1), nt, mpi_real,  &

             icroot-1, mpi_comm_mwave, ieer)

#endif


        ! 1.d  Loop over J grids, select same rank

        !

        DO j=1, nrgrd


          IF( .NOT. shrank(i,j) ) cycle

          !! Only SMC J-grid provides boundary spectra for I-Grid.

          IF( grids(j)%GTYPE .NE. smctype ) cycle

          !

          ! Find local root PE and NAPROC for J grid.

#ifdef W3_SHRD

          jcroot = 1

#endif

#ifdef W3_MPI

          jcroot = mdatas(j)%CROOT

#endif

          npj    = outpts(j)%NAPROC

          !

          ! Find out whether any I-grid boundary points matched in J-Grid.

          !

          stores(i,j)%INIT = .true.

          ALLOCATE( stores(i,j)%ICVBP(nt), stores(i,j)%MSDBP(nt), &

               stores(i,j)%JCVBP(nt), stores(i,j)%IPCVB(nt), &

               stores(i,j)%ISS(nt),   stores(i,j)%JSS(nt),   &

               stores(i,j)%IPS(nt),   stores(i,j)%ITG(nt),   &

               stores(i,j)%FLA(nt),   stat=istat )

          check_alloc_status( istat )

          stores(i,j)%ICVBP = 0

          stores(i,j)%JCVBP = 0

          stores(i,j)%IPCVB = 0

          stores(i,j)%MSDBP = 0

          stores(i,j)%ISS   = 0

          stores(i,j)%JSS   = 0

          stores(i,j)%IPS   = 0

          stores(i,j)%ITG   = 0

          stores(i,j)%FLA   = .false.

          !

          !  Work out which I-grid bounary points are matched in J-grid on JCROOT.

#ifdef W3_MPI

#ifdef W3_SMC

          IF( improc .EQ. jcroot ) THEN

#endif

#endif

#ifdef W3_SMC

            CALL w3smcgmp( j, nt, xlon, ylat, stores(i,j)%MSDBP )

#endif

#ifdef W3_MPI

#ifdef W3_SMC

          ENDIF

#endif

#endif

          !

          ! Then broadcast the results to all PEs

#ifdef W3_MPI

#ifdef W3_SMC

          CALL mpi_bcast( stores(i,j)%MSDBP(1), nt, mpi_integer,  &

               jcroot-1, mpi_comm_mwave, ieer)

#endif

#endif

          !

          !  Need to wait for all PEs get these values.

#ifdef W3_MPI

#ifdef W3_SMC

          CALL mpi_barrier( mpi_comm_mwave, ieer)

#endif

#endif

          !

#ifdef W3_SMC

          stores(i,j)%ICVBP = ibpts

          stores(i,j)%JCVBP = jbpts

          stores(i,j)%IPCVB = ipbpt

#endif

          !

          !  Check which I-grid boundary points matched inside J-Grid

          ntl= 0

          DO jx=1, nt

            IF( stores(i,j)%MSDBP(jx) .EQ. 0 )  cycle


            !  Process J-grid send point if it matches I-grid boundary point.

            ntl    = ntl + 1

            itag   = itag + 1

            isea   = stores(i,j)%MSDBP(jx)

            !  Find global processor IPRC and local JSEA on J-grid, holding ISEA spectrum.

            jsea = 1 + (isea - 1)/npj

            iprc = jcroot - 1 + isea - (jsea - 1)*npj

            !  Store these spectral location info in STORES.

            stores(i,j)%ISS(jx) = isea

            stores(i,j)%JSS(jx) = jsea

            stores(i,j)%IPS(jx) = iprc

            stores(i,j)%ITG(jx) = itag

            stores(i,j)%FLA(jx) = .true.

          END DO

          !

          !  SMC grid boundary points are supposed to be 1 to 1 correspondant

          !  so there is no need for interpolation.   JGLi03Nov2020

          !

          stores(i,j)%NTOT = nt

          stores(i,j)%NFIN = ntl

          !

#ifdef W3_MPI

#ifdef W3_SMC

          IF( improc .EQ. nmperr )              &

#endif

#endif

#ifdef W3_SMC

               WRITE(mdse,1060) i, nt, j, ntl

#endif

          !

          ! ... End of loops J in 1.c

        END DO

        !

        !!  Free temporary space for I-grid.

        DEALLOCATE( ibpts, jbpts, ipbpt, xlon, ylat, stat=istat )

        check_dealloc_status( istat )


      END IF   ! NT > 0

      !

      ! ... End of 1.a loop I grid.

    END DO

    !

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

    ! 3.  Final data base (full data base, scratched at end of routine)

    ! 3.a Loop over grids

    !

    ALLOCATE( nrec(nrgrd), nsnd(nrgrd), ntpp(nmproc), stat=istat )

    check_alloc_status( istat )

    !

    DO i=1, nrgrd

      IF ( .NOT. dogrid(i) ) cycle

#ifdef W3_T

      WRITE (mdst,9030) i

#endif

      !

      CALL w3setg ( i, mdse, mdst )

      CALL w3seto ( i, mdse, mdst )

      CALL wmsetm ( i, mdse, mdst )

      !

      nrec   = 0

      nsnd   = 0

      !

      ! Find local root PE and maximum PE for I grid.

#ifdef W3_SHRD

      icroot = 1

#endif

#ifdef W3_MPI

      icroot = mdatas(i)%CROOT

#endif

      npmx   = outpts(i)%NAPROC + icroot - 1

      !

      ! 3.b Filling NREC and NSND for grid I

      !

      !!  Work out how many I-grid boundary points to be updated by other grids.

      !!  Use matched J-grid points to selected I-grid points. JGLi26Jan2021

      DO j=1, nrgrd

        IF( .NOT. shrank(i,j) ) cycle

        IF( stores(i,j)%NFIN > 0 ) THEN

          DO ix = 1, stores(i,j)%NTOT

            IF( stores(i,j)%MSDBP(ix) > 0  .AND.  &

                 stores(i,j)%IPCVB(ix) .EQ. improc ) THEN

              nrec(i) = nrec(i) + 1

              nrec(j) = nrec(j) + 1

            END IF

          END DO

        END IF

      END DO


      !  Accumulate all related I-Grid points to be send to other sub-grids.

      !  Add IPRC range check to ensure sending from I-grid.  JGLi22Jan2021

      DO j=1, nrgrd

        IF( .NOT. shrank(j,i) ) cycle

        IF( stores(j,i)%NFIN > 0 ) THEN

          DO iy=1, stores(j,i)%NTOT

            IF( stores(j,i)%MSDBP(iy) > 0  .AND.  &

                 stores(j,i)%IPS(  iy) .EQ. improc ) THEN

              nsnd(j) = nsnd(j) + 1

            ENDIF

          END DO

        END IF

      END DO

      !

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

      ! 4.  Save data base as needed in EQSTGE

      !

      ! 4.a   ALLOCATE storage

      ! 4.a.1 Local counters, weights and sea counters (grid 'I')

      !

      IF( eqstge(i,i)%NREC .NE. 0 ) THEN

        DEALLOCATE (eqstge(i,i)%ISEA, eqstge(i,i)%JSEA ,         &

             eqstge(i,i)%WGHT, stat=istat )

        check_dealloc_status( istat )

        eqstge(i,i)%NREC = 0

#ifdef W3_T

        WRITE (mdst,9040) i, i

#endif

      END IF

      !

      IF( nrec(i) .GT. 0 ) THEN

        ALLOCATE( eqstge(i,i)%ISEA(nrec(i)),                   &

             eqstge(i,i)%JSEA(nrec(i)),                   &

             eqstge(i,i)%WGHT(nrec(i)), stat=istat )

        check_alloc_status( istat )

        eqstge(i,i)%NREC = nrec(i)

#ifdef W3_T

        WRITE (mdst,9041) i, i, nrec(i)

#endif

      END IF

      !

      !! Initial NTOT for grid I before summing over other grids.  JGLi18Jan2021

      eqstge(i,i)%NTOT = 0

      !

      ! 4.a.1 Local counters, arrays weights etc.  (grid 'J' receive)

      !

      DO j=1, nrgrd

        IF( i .EQ. j ) cycle

        !

        !!   Looks strange to store in EQSTGE(I,I) as other J-grid may

        !!   overwrite the value. Should be suspended?  JGLi30Dec2020

        !         EQSTGE(I,I)%NTOT = STORES(I,J)%NFIN

        !!   Changed to summation ove all other J-grids NFIN.  Not sure where

        !!   NTOT is used but keep it anyway.   JGLi18Jan2021

        eqstge(i,i)%NTOT = eqstge(i,i)%NTOT + stores(i,j)%NFIN

        !

        IF( eqstge(i,j)%NREC .NE. 0 ) THEN

          DEALLOCATE( eqstge(i,j)%ISEA , eqstge(i,j)%JSEA ,      &

               eqstge(i,j)%WGHT , eqstge(i,j)%SEQL ,      &

               eqstge(i,j)%NAVG , eqstge(i,j)%WAVG ,      &

               eqstge(i,j)%RIP  , eqstge(i,j)%RTG,        &

               stat=istat )

          check_dealloc_status( istat )

          eqstge(i,j)%NREC   = 0

          eqstge(i,j)%NAVMAX = 1

        END IF

        !

        IF( nrec(j) .GT. 0 ) THEN

          na = 1

          eqstge(i,j)%NAVMAX = na

          ALLOCATE( eqstge(i,j)%ISEA(nrec(j)),                   &

               eqstge(i,j)%JSEA(nrec(j)),                   &

               eqstge(i,j)%WGHT(nrec(j)),                   &

               eqstge(i,j)%SEQL(sgrds(j)%NSPEC,nrec(j),na), &

               eqstge(i,j)%NAVG(nrec(j)),                   &

               eqstge(i,j)%WAVG(nrec(j),na),                &

               eqstge(i,j)%RIP( nrec(j),na),                &

               eqstge(i,j)%RTG( nrec(j),na), stat=istat )

          check_alloc_status( istat )

          eqstge(i,j)%NREC = nrec(j)

        END IF

        !

        IF( eqstge(j,i)%NSND .NE. 0 ) THEN

          DEALLOCATE( eqstge(j,i)%SIS, eqstge(j,i)%SJS ,        &

               eqstge(j,i)%SI1, eqstge(j,i)%SI2 ,        &

               eqstge(j,i)%SIP, eqstge(j,i)%STG, stat=istat)

          check_dealloc_status( istat )

          eqstge(j,i)%NSND = 0

        END IF

        !

        IF( nsnd(j) .GT. 0 ) THEN

          ALLOCATE( eqstge(j,i)%SIS(nsnd(j)),                   &

               eqstge(j,i)%SJS(nsnd(j)),                   &

               eqstge(j,i)%SI1(nsnd(j)),                   &

               eqstge(j,i)%SI2(nsnd(j)),                   &

               eqstge(j,i)%SIP(nsnd(j)),                   &

               eqstge(j,i)%STG(nsnd(j)), stat=istat )

          check_alloc_status( istat )

          eqstge(j,i)%NSND = nsnd(j)

        END IF

        !

      END DO

      !

      ! 4.b   Store data in EQSTGE

      ! 4.b.1 Grid I (JSEA and weight only) also filled in J-Grid loop

      !       but it accumulates all points received by I-grid.

      nt = 0

      !

      ! 4.b.2 Info for I-grid receiving from all other grids

      DO j=1, nrgrd

        IF( .NOT. shrank(i,j) ) cycle

        IF( eqstge(i,j)%NREC .EQ. 0 ) cycle

        ntl = 0

        DO ix=1, stores(i,j)%NTOT

          IF(   stores(i,j)%MSDBP(ix) > 0  .AND.  &

               stores(i,j)%IPCVB(ix) .EQ. improc ) THEN

            ! All points received by I-grid accumulated from each J-grid.

            nt = nt + 1

            eqstge(i,i)%ISEA(nt) = stores(i,j)%ICVBP(ix)

            eqstge(i,i)%JSEA(nt) = stores(i,j)%JCVBP(ix)

            ! No need to alter local spectra for SMC grid.  JGLi08Dec2020

            eqstge(i,i)%WGHT(nt) = 1.0


            ! Boundary points received by I-grid from J-grid.

            ntl = ntl + 1

            eqstge(i,j)%ISEA(ntl) = stores(i,j)%ICVBP(ix)

            eqstge(i,j)%JSEA(ntl) = stores(i,j)%JCVBP(ix)

            !! Boundary spectra will be substituted fully.  JGLi08Dec2020

            eqstge(i,j)%WGHT(ntl) = 1.0

            eqstge(i,j)%NAVG(ntl) = 1

            eqstge(i,j)%WAVG(ntl,1) = 1.0

            eqstge(i,j)%RIP (ntl,1) = stores(i,j)%IPS(ix)

            eqstge(i,j)%RTG (ntl,1) = stores(i,j)%ITG(ix)

          END IF

        END DO


      END DO

      !

      ! 4.b.3 All other grids, info for sending

      !

      DO j=1, nrgrd

        IF ( .NOT. shrank(j,i) ) cycle

        IF ( eqstge(j,i)%NSND .EQ. 0 ) cycle

        ntpp   = 0

        ntl    = 0

        DO iy =1, stores(j,i)%NTOT

          IF(    stores(j,i)%MSDBP(iy) > 0 ) THEN

            iprc=stores(j,i)%IPS(  iy)

            ntpp(iprc) = ntpp(iprc) + 1

            IF( iprc .EQ. improc ) THEN

              ntl    = ntl + 1

              eqstge(j,i)%SIS(ntl) = stores(j,i)%ISS(iy)

              eqstge(j,i)%SJS(ntl) = stores(j,i)%JSS(iy)

              eqstge(j,i)%SI1(ntl) = ntpp(iprc)

              eqstge(j,i)%SI2(ntl) = 1

              eqstge(j,i)%SIP(ntl) = stores(j,i)%IPCVB(iy)

              eqstge(j,i)%STG(ntl) = stores(j,i)%ITG(iy)

            END IF

          END IF

        END DO

        !

      END DO

      !

      !  End of 3.a loop for I grid.

    END DO

    !

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

    ! 5.  Generate GRDEQL

    ! 5.a Size of array

    !

    nrec   = 0

    !

    DO i=1, nrgrd

      DO j=1, nrgrd

        IF ( i.EQ.j .OR. stores(i,j)%NFIN.EQ.0 ) cycle

        nrec(i) = nrec(i) + 1

      END DO

    END DO

    !

    na = maxval(nrec)

    ALLOCATE( grdeql(nrgrd,0:na), stat=istat )

    check_alloc_status( istat )

    grdeql = 0

    !

#ifdef W3_T

    WRITE (mdst,9050) na

#endif

    !

    ! 5.b Fill array

    !

    DO i=1, nrgrd

      grdeql(i,0) = nrec(i)

      jj          = 0

      DO j=1, nrgrd

        IF ( i.EQ.j .OR. stores(i,j)%NFIN.EQ.0 ) cycle

        jj           = jj + 1

        grdeql(i,jj) = j

      END DO

    END DO

    !

#ifdef W3_T

    WRITE (mdst,9051)

    DO i=1, nrgrd

      WRITE (mdst,9052) i, grdeql(i,0:grdeql(i,0))

    END DO

#endif

    !

    ! 5.d Group number test

    !

    DO i=1, nrgrd

      IF( grdeql(i,0) .GE. 2 ) THEN

        tgrp = grgrp(grdeql(i,1))

        DO j=2, grdeql(i,0)

          IF( grgrp(grdeql(i,j)) .NE. tgrp ) THEN

            IF( improc .EQ. nmperr ) WRITE(mdse,1051)          &

                 grdeql(i,1), grgrp(grdeql(i,1)),   &

                 grdeql(i,j), grgrp(grdeql(i,j))

            CALL extcde ( 1051 )

          END IF

        END DO

      END IF

    END DO

    !

    !  Wait all PEs finishing EQSTGE setup before clean up.  JGLi20Jan2021

#ifdef W3_MPI

#ifdef W3_SMC

    CALL mpi_barrier (mpi_comm_mwave,ieer)

#endif

#endif

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

    ! 6.  Final clean up

    !

    DO i=1, nrgrd

      DO j=1, nrgrd

        IF( stores(i,j)%INIT ) THEN

          DEALLOCATE( stores(i,j)%ICVBP, stores(i,j)%MSDBP,   &

               stores(i,j)%JCVBP, stores(i,j)%IPCVB,   &

               stores(i,j)%ISS  , stores(i,j)%JSS  ,   &

               stores(i,j)%IPS  , stores(i,j)%ITG  ,   &

               stores(i,j)%FLA  , stat=istat )

          check_dealloc_status( istat )

        END IF

      END DO

    END DO

    !

    DEALLOCATE( shrank, stores, nrec, nsnd, ntpp, stat=istat )

    check_dealloc_status( istat )

    !

#ifdef W3_MPI

#ifdef W3_SMC

    IF( improc .EQ. nmperr )              &

#endif

#endif

#ifdef W3_SMC

         WRITE(mdse,*) " *** WMSMCEQL completed from PE ", improc

#endif


    RETURN

    !

    ! Formats

    !

1000 FORMAT (/' *** WAVEWATCH III ERROR IN WMSMCEQL : *** '/         &

         '     UNCOVERED EDGE POINTS FOR GRID',i4,'  (',i6,')'/)

1001 FORMAT ( '     GRID',i4,'  POINT',2i5,' NOT COVERED (WMGEQL)')

1002 FORMAT ( '     DIAGNOSTICS IX AND IY RANGE:',4i6)

1003 FORMAT (/'     SHOWING ',a/)

1004 FORMAT (2x,65i2)

1005 FORMAT (/'     SHOWING IX RANGE ',2i6)

1006 FORMAT ( '        (WILL NOT PRINT ANY MORE UNCOVERED POINTS)')

    !

1020 FORMAT (/' *** WAVEWATCH III WARNING WMSMCEQL : *** '/          &

         '     REMOVED BOUNDARY POINT FROM OPEN EDGE DISTANCE MAP'/     &

         '     GRID, IX, IY :',3i6)

    !

1050 FORMAT (/' *** WAVEWATCH III ERROR IN WMSMCEQL : *** '/         &

         '     GRID INCREMENTS TOO DIFFERENT '/                 &

         '     GRID',i4,'   INCREMENTS ',2f8.2/                 &

         '     GRID',i4,'   INCREMENTS ',2f8.2/)

1051 FORMAT (/' *** WAVEWATCH III ERROR IN WMSMCEQL : *** '/         &

         '     OVERLAPPING GRIDS NEED TO BE IN SAME GROUP '/    &

         '         GRID',i4,' IN GROUP',i4/                     &

         '         GRID',i4,' IN GROUP',i4/)

1060 FORMAT (' Grid NBPI from',2i6,' found in',2i6)


    !

#ifdef W3_T

9010 FORMAT ( ' TEST WMSMCEQL : STARTING LOOP OVER GRIDS')

9011 FORMAT ( ' TEST WMSMCEQL : I, RANK :',2i4)

9012 FORMAT ( '               GRID ',i3,' HAS SAME RANK')

9013 FORMAT ( '               ',a)

#endif

    !

#ifdef W3_T

9020 FORMAT ( ' TEST WMSMCEQL : GENERATING DISTANCE MAP GRID ',i3)

9024 FORMAT ( ' TEST WMSMCEQL : FINAL MAP FOR X RANGE ',2i6)

9025 FORMAT (2x,65i2)

#endif

    !

#ifdef W3_T

9030 FORMAT ( ' TEST WMSMCEQL : DEPENDENCE INFORMATION GRID ',i3)

9031 FORMAT ( '               CHECKING GRID ',i3)

9032 FORMAT ( '               POINTS USED/AVAIL :',2i6)

9033 FORMAT ( '               TOTAL/GRIDS/OUT   :',3i6)

9034 FORMAT ( '               LOCAL PER GRID    :',15i6)

9035 FORMAT ( '               SENDING PER GRID  :',15i6)

9036 FORMAT ( ' TEST WMSMCEQL : NUMBER OF CONTRIBUTING GRIDS MAP')

9037 FORMAT (2x,65i2)

#endif

    !

#ifdef W3_T

9040 FORMAT ( ' TEST WMSMCEQL : GRID ',i2,'-',i2,' CLEAR STORAGE')

9041 FORMAT ( ' TEST WMSMCEQL : GRID ',i2,'-',i2,' STORAGE SIZE',i6)

9042 FORMAT ( '               RECV ',i2,'-',i2,' CLEAR STORAGE')

9043 FORMAT ( '               RECV ',i2,'-',i2,' STORAGE SIZE',2i6)

9044 FORMAT ( '               SEND ',i2,'-',i2,' CLEAR STORAGE')

9045 FORMAT ( '               SEND ',i2,'-',i2,' STORAGE SIZE',i6)

#endif

    !

#ifdef W3_T

9050 FORMAT ( ' TEST WMSMCEQL : GRDEQL DIMENSIONED AT ',i2)

9051 FORMAT ( ' TEST WMSMCEQL : GRDEQL :')

9052 FORMAT ( '                 ',2i4,' : ',20i3)

#endif

    !

    !/

    !/ End of WMSMCEQL  -------------------------------------------------- /

    !/

  END SUBROUTINE wmsmceql

  !!


  !/ End of module WMGRIDMD -------------------------------------------- /

  !/

END MODULE wmgridmd