wmgridmd Module Reference

Routines to determine and process grid dependencies in the multi-grid wave model. More...

Functions/Subroutines
subroutine	wmglow (FLRBPI)
	Determine relations to lower ranked grids for each grid. More...
subroutine	wmghgh
	Determine relation to higher ranked grids for each grid. More...
subroutine	wmgeql
	Determine relations to same ranked grids for each grid. More...
subroutine	wmrspc
	Generate map with flags for need of spectral grid conversion between models. More...
subroutine	wmsmceql
	Determine relations to same ranked SMC grids for each grid. More...

Detailed Description

Routines to determine and process grid dependencies in the multi-grid wave model.

Author

H. L. Tolman

W. E. Rogers

Date

10-Dec-2014

Function/Subroutine Documentation

wmgeql()

subroutine wmgridmd::wmgeql

Determine relations to same ranked grids for each grid.

Cross mapping of grid points, determine boundary distance data and interpolation weights.

Author

H. L. Tolman

Date

10-Dec-2014

Definition at line 3709 of file wmgridmd.F90.

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

References w3gdatmd::clgtype, wmmdatmd::croot, constants::dera, w3gdatmd::dtcfl, w3gdatmd::dtmax, wmmdatmd::eqstge, w3servmd::extcde(), w3gdatmd::flagll, w3gdatmd::flagst, wmmdatmd::flghg1, wmmdatmd::grank, constants::grav, wmmdatmd::grdeql, wmmdatmd::grgrp, w3gdatmd::grids, w3gdatmd::gtype, w3gdatmd::hpfac, w3gdatmd::hqfac, w3gdatmd::iclose, w3gdatmd::iclose_none, w3gdatmd::iclose_trpl, wmmdatmd::improc, wmmdatmd::init_get_jsea_isproc_glob(), w3gdatmd::mapfs, wmmdatmd::mapodi, w3gdatmd::mapsf, w3gdatmd::mapst2, w3gdatmd::mapsta, wmmdatmd::mdatas, wmmdatmd::mdse, wmmdatmd::mdst, w3odatmd::naproc, wmmdatmd::nmperr, wmmdatmd::nmproc, wmmdatmd::nrgrd, w3gdatmd::nsea, w3gdatmd::nx, w3gdatmd::ny, constants::radius, w3gdatmd::rlgtype, w3gdatmd::sgrds, w3gdatmd::sig, w3gdatmd::smctype, w3servmd::strace(), w3gdatmd::sx, w3gdatmd::sy, w3gdatmd::ungtype, w3gdatmd::w3setg(), w3odatmd::w3seto(), wmmdatmd::wmsetm(), w3gdatmd::x0, and w3gdatmd::y0.

Referenced by wminitmd::wminit(), and wminitmd::wminitnml().

wmghgh()

subroutine wmgridmd::wmghgh

Determine relation to higher ranked grids for each grid.

Base map set in WMGLOW, supplemental data computed here. Map averaging information for higher ranked grid to lower ranked grid.

Author

H. L. Tolman

W. E. Rogers

Date

10-Dec-2014

Definition at line 1100 of file wmgridmd.F90.

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

References w3gdatmd::clgtype, wmmdatmd::croot, constants::dera, w3gdatmd::dtcfl, w3gdatmd::dtmax, w3servmd::extcde(), file(), w3gdatmd::flagll, w3gdatmd::flagst, wmmdatmd::flgbdi, wmmdatmd::flghg1, wmmdatmd::grank, constants::grav, wmmdatmd::grdhgh, w3gdatmd::grids, w3gdatmd::gtype, wmmdatmd::hgstge, w3gdatmd::iclose, w3gdatmd::iclose_none, w3gdatmd::iclose_trpl, wmmdatmd::improc, include(), w3parall::init_get_jsea_isproc(), wmmdatmd::init_get_jsea_isproc_glob(), wmmdatmd::mapbdi, w3gdatmd::mapfs, wmmdatmd::mapmsk, w3gdatmd::mapsf, w3gdatmd::mapst2, w3gdatmd::mapsta, wmmdatmd::mdatas, wmmdatmd::mdse, wmmdatmd::mdst, wmmdatmd::mpi_comm_mwave, w3odatmd::nbi, wmmdatmd::nmperr, wmmdatmd::nmproc, wmmdatmd::nrgrd, w3gdatmd::nsea, w3gdatmd::nx, w3gdatmd::ny, constants::radius, w3gdatmd::rlgtype, wmscrpmd::scrip_wrapper(), w3gdatmd::sgrds, w3gdatmd::sig, w3gdatmd::smctype, w3servmd::strace(), w3gdatmd::sx, w3gdatmd::sy, w3gdatmd::ungtype, w3gdatmd::w3setg(), w3odatmd::w3seto(), scrip_interface::wgtdata, wmmdatmd::wmsetm(), w3gdatmd::x0, w3gdatmd::xgrd, and w3gdatmd::y0.

Referenced by wminitmd::wminit(), and wminitmd::wminitnml().

wmglow()

subroutine wmgridmd::wmglow

(

logical, dimension(nrgrd), intent(out), optional

FLRBPI

)

Determine relations to lower ranked grids for each grid.

On the fly, the opposite relations are also saved. Map active boundary points to lower ranked grids.

Parameters

[out]

FLRBPI

Array with flags for external file use.

Author

H. L. Tolman

W. E. Rogers

Date

06-Jun-2018

Definition at line 152 of file wmgridmd.F90.

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

References w3gdatmd::clgtype, w3servmd::extcde(), w3gdatmd::flagll, wmmdatmd::grank, wmmdatmd::grdhgh, wmmdatmd::grdlow, w3gdatmd::grids, w3odatmd::iaproc, wmmdatmd::improc, include(), w3parall::init_get_jsea_isproc(), w3odatmd::ipbpi, w3triamd::is_in_ungrid(), w3odatmd::isbpi, w3gdatmd::mapfs, w3gdatmd::mapsf, w3gdatmd::mapsta, wmmdatmd::mdatas, wmmdatmd::mdse, wmmdatmd::mdst, wmmdatmd::mpi_comm_mwave, w3odatmd::naperr, w3odatmd::nbi, w3odatmd::nbi2, wmmdatmd::nbi2g, wmmdatmd::nbi2s, w3odatmd::ndse, w3gdatmd::nk, wmmdatmd::nmperr, w3gdatmd::nsea, w3gdatmd::nx, w3gdatmd::ny, w3odatmd::outpts, w3odatmd::rdbpi, w3gdatmd::rlgtype, w3gdatmd::sgrds, w3gdatmd::smctype, w3servmd::strace(), w3gdatmd::ungtype, w3odatmd::w3dmo5(), w3gdatmd::w3setg(), w3odatmd::w3seto(), w3odatmd::xbpi, w3gdatmd::xgrd, w3odatmd::ybpi, and w3gdatmd::ygrd.

Referenced by wminitmd::wminit(), and wminitmd::wminitnml().

wmrspc()

subroutine wmgridmd::wmrspc

Generate map with flags for need of spectral grid conversion between models.

Test of parameters as introduced before in W3IOBC.

Author

H. L. Tolman

Date

10-Dec-2014

Definition at line 5139 of file wmgridmd.F90.

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

References w3servmd::extcde(), wmmdatmd::mdst, wmmdatmd::nrgrd, wmmdatmd::respec, w3gdatmd::sgrds, w3servmd::strace(), and w3odatmd::unipts.

Referenced by wminitmd::wminit(), and wminitmd::wminitnml().

wmsmceql()

subroutine wmgridmd::wmsmceql

Determine relations to same ranked SMC grids for each grid.

Set boundary points update for regular grid in same ranked group.

Cross mapping of grid points, use nearest sea points and no interpolation is required so far.

Author

J G Li

Date

12-Apr-2021

Definition at line 5285 of file wmgridmd.F90.

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

References wmmdatmd::eqstge, w3servmd::extcde(), wmmdatmd::grank, wmmdatmd::grdeql, wmmdatmd::grgrp, w3gdatmd::grids, wmmdatmd::improc, include(), w3gdatmd::mapsf, w3gdatmd::mapsta, wmmdatmd::mdatas, wmmdatmd::mdse, wmmdatmd::mdst, wmmdatmd::mpi_comm_mwave, wmmdatmd::nmperr, wmmdatmd::nmproc, wmmdatmd::nrgrd, w3gdatmd::nsea, w3odatmd::outpts, w3gdatmd::rlgtype, w3gdatmd::sgrds, w3gdatmd::smctype, w3servmd::strace(), w3gdatmd::w3setg(), w3odatmd::w3seto(), w3psmcmd::w3smcell(), w3psmcmd::w3smcgmp(), wmmdatmd::wmsetm(), w3gdatmd::xgrd, and w3gdatmd::ygrd.

Referenced by wminitmd::wminit(), and wminitmd::wminitnml().