w3psmcmd.F90

Go to the documentation of this file.

 
 
#include "w3macros.h"
!/ ------------------------------------------------------------------- /
MODULE w3psmcmd
  !/
  !/                  +------------------------------------+
  !/                  | Spherical Multiple-Cell (SMC) grid |
  !/                  | Adv, GCT, Rfr, Dif subroutines.    |
  !/                  |           Jian-Guo Li              |
  !/                  | First created:     8 Nov 2010      |
  !/                  | Last modified:    23 Mar 2020      |
  !/                  +------------------------------------+
  !/
  !/    08-Nov-2010 : Coding started by adapting w3pro2md.ftn.
  !/    18-Nov-2010 : Refraction and GCT by rotation and k-shift.
  !/    12-Apr-2011 : Restore x-advective flux for intermediate update.
  !/     3-Jun-2011 : New refraction formulation using Cg only.
  !/     8-Jun-2011 : Optimise classic refraction formulation.
  !/    16-Jun-2011 : Add refraction limter to gradient direction.
  !/     1-Jul-2011 : New refraction using Cg and gradient limiter.
  !/    28-Jul-2011 : Finalise with old refraction scheme and gradient limiter.
  !/     4-Nov-2011 : Separate x and y obstruction coefficients.
  !/     5-Jan-2012 : Update to multi-resolution SMC grid with sub-time-steps.
  !/     2-Feb-2012 : Separate single- and multi-resolution advection.
  !/     6-Mar-2012 : Tidy up code and minor adjustments, CLATF.
  !/    12-Mar-2012 : Remove net flux bug and optimise upstream code.
  !/    16-Jan-2013 : Adapted for Version 4.08, removing FACX/Y.
  !/    16-Sep-2013 : Add Arctic part for SMC grid in WW3 V4.11
  !/     3-Jan-2014 : Remove bug in SMCDHXY for AU/V as cell size.
  !/     7-Jan-2014 : Remove bug in SMCGtCrfr for K definition.
  !/    28-Jan-2014 : Move Arctic boundary condition update out.
  !/    18-Aug-2015 : New gradient, average and 3rd order advection subs.
  !/     3-Sep-2015 : UNO3 advection scheme by logical option FUNO3.
  !/    14-Sep-2015 : Modify DHDX/Y for Arctic part refraction term.
  !/     8-Aug-2017 : Update SMCGradn for 0 or 1 boundary conditions.
  !/     9-Jan-2018 : Parallelization by adding OpenMP directives.
  !/    19-Feb-2020 : Additions for OMP bit-reproducability (C.Bunney)
  !/    23-Mar-2020 : Add extra parenthese for single ATOMIC line update.
  !/    22-Oct-2020 : Two new subs for lat-lon points mapping to cells.
  !/
  !  1. Purpose :
  !
  !     Bundles routines for SMC advection (UNO2) and diffusion schemes in
  !     single module, including GCT and refraction rotation schemes.
  !
  !  2. Variables and types :
  !
  !      Name      Type  Scope    Description
  !     ----------------------------------------------------------------
  !      TRNMIN    R.P.  Private   Minimum transparancy for local
  !     ----------------------------------------------------------------
  !
  !  3. Subroutines and functions :
  !
  !      Name      Type  Scope    Description
  !     ----------------------------------------------------------------
  !      W3PSMC    Subr. Public   Spatial propagation on SMC grid.
  !      W3KRTN    Subr. Public   Spectral modification by GCT and refraction.
  !      SMCxUNO2  Subr. Public   Irregular grid mid-flux on U-faces by UNO2.
  !      SMCyUNO2  Subr. Public   Irregular grid mid-flux on V-faces by UNO2.
  !      SMCxUNO2r Subr. Public   Regular grid mid-flux on U-faces by UNO2.
  !      SMCyUNO2r Subr. Public   Regular grid mid-flux on V-faces by UNO2.
  !      SMCkUNO2  Subr. Public   Shift in k-space due to refraction by UNO2.
  !      SMCxUNO3  Subr. Public   Irregular grid mid-flux on U-faces by UNO3.
  !      SMCyUNO3  Subr. Public   Irregular grid mid-flux on V-faces by UNO3.
  !      SMCxUNO3r Subr. Public   Regular grid 3rd order U-mid-flux by UNO3.
  !      SMCyUNO3r Subr. Public   Regular grid 3rd order V-mid-flux by UNO3.
  !      SMCGtCrfr Subr. Public   Refraction and GCT rotation in theta.
  !      SMCDHXY   Subr. Public   Evaluate depth gradient and refraction limiter.
  !      SMCGradn  Subr. Public   Evaluate local gradient for sea points.
  !      SMCAverg  Subr. Public   Numerical 1-2-1 average of sea point field.
  !      W3GATHSMC W3SCATSMC      Gather and scatter spectral components.
  !      W3SMCELL  Subr. Public   Calculate cell centre lat-lon for given ids.
  !      W3SMCGMP  Subr. Public   Map lat-lon points to SMC grid cells.
  !     ----------------------------------------------------------------
  !
  !  4. Subroutines and functions used :
  !
  !      Name      Type  Module   Description
  !     ----------------------------------------------------------------
  !      STRACE    Subr. W3SERVMD Subroutine tracing.
  !      W3ACTURN  Subr. W3SERVMD Subroutine rotating action spectrum.
  !     ----------------------------------------------------------------
  !
  !  5. Remarks :
  !
  !  6. Switches :
  !
  !       !/MGP   Correct for motion of grid.
  !       !/S     Enable subroutine tracing.
  !       !/Tn    Enable test output.
  !
  !  7. Source code :
  !
  !/ ------------------------------------------------------------------- /
  !/
  !/ Use omp_lib for OpenMP functions if switched on.   JGLi10Jan2018
  !$       USE omp_lib
  !/
  !/ Public variables
  !/
  PUBLIC
  !/
  !/ Private data !/
  REAL, PRIVATE, PARAMETER:: TRNMIN = 0.95
  !/
CONTAINS
 
  !/ ------------------------------------------------------------------- /
  SUBROUTINE w3psmc ( ISP, DTG, VQ )
    !/
    !/                  +------------------------------------+
    !/                  | Spherical Multiple-Cell (SMC) grid |
    !/                  | Advection and diffusion sub.       |
    !/                  | First created:   JG Li  8 Nov 2010 |
    !/                  | Last modified:   JG Li 18 Apr 2018 |
    !/                  +------------------------------------+
    !/
    !/    08-Nov-2010 : Origination.                JGLi    ( version 1.00 )
    !/    16-Dec-2010 : Check U/V CFL values.       JGLi    ( version 1.10 )
    !/    18-Mar-2011 : Check MPI communication.    JGLi    ( version 1.20 )
    !/    16-May-2011 : Tidy up diagnosis lines.    JGLi    ( version 1.30 )
    !/     4 Nov-2011 : Separate x and y obstruc.   JGLi    ( version 1.40 )
    !/     5 Jan-2012 : Multi-resolution SMC grid.  JGLi    ( version 1.50 )
    !/     2 Feb-2012 : Separate single multi adv.  JGLi    ( version 1.60 )
    !/     6 Mar-2012 : Minor adjustments of CLATF. JGLi    ( version 1.70 )
    !/    12 Feb-2012 : Remove net flux bug.        JGLi    ( version 1.80 )
    !/    16 Sep-2013 : Add Arctic part.            JGLi    ( version 2.00 )
    !/     3 Sep-2015 : Gradient, UNO3 and Average. JGLi    ( version 2.10 )
    !/    26 Feb-2016 : Update boundary spectra.    JGLi    ( version 2.20 )
    !/    23 Mar-2016 : Add current option.         JGLi    ( version 2.30 )
    !/    18 Apr-2018 : Refined sub-grid blocking.  JGLi    ( version 2.40 )
    !/
    !  1. Purpose :
    !
    !     Propagation in phyiscal space for a given spectral component.
    !
    !  2. Method :
    !
    !     Unstructured SMC grid, point-oriented face and cell loops.
    !     UNO2 advection scheme and isotropic FTCS diffusion scheme.
    !
    !  3. Parameters :
    !
    !     Parameter list
    !     ----------------------------------------------------------------
    !       ISP     Int.   I   Number of spectral bin (IK-1)*NTH+ITH
    !       FACX/Y  Real   I   Factor in propagation velocity (1 or 0 *DT/DX)
    !       DTG     Real   I   Total time step.
    !       MAPSTA  I.A.   I   Grid point status map.
    !       MAPFS   I.A.   I   Storage map.
    !       VQ      R.A.  I/O  Field to propagate.
    !     ----------------------------------------------------------------
    !
    !     Local variables.
    !     ----------------------------------------------------------------
    !       NTLOC   Int   Number of local time steps.
    !       DTLOC   Real  Local propagation time step.
    !       CGD     Real  Deep water group velocity.
    !       DSSD, DNND    Deep water diffusion coefficients.
    !       ULCFLX  R.A.  Local courant numbers in 'x' (norm. velocities)
    !       VLCFLY  R.A.  Id. in 'y'.
    !       CXTOT   R.A.  Propagation velocities in physical space.
    !       CYTOT   R.A.
    !       DFRR    Real  Relative frequency increment.
    !       DX0I    Real  Inverted grid incremenent in meters (longitude, eq.).
    !       DY0I    Real  Inverted grid incremenent in meters (latitude).
    !     ----------------------------------------------------------------
    !
    !  4. Subroutines used :
    !
    !     See module documentation.
    !
    !  5. Called by :
    !
    !      Name      Type  Module   Description
    !     ----------------------------------------------------------------
    !      W3WAVE    Subr. W3WAVEMD Wave model routine.
    !     ---------------------------------------------------------------
    !
    !  6. Error messages :
    !
    !  7. Remarks :
    !
    !  8. Structure :
    !
    !     ---------------------------------------------
    !       1.  Preparations
    !         a Set constants
    !         b Initialize arrays
    !       2.  Prepare arrays
    !         a Velocities and 'Q'
    !         b diffusion coefficients
    !       3.  Loop over sub-steps
    !       ----------------------------------------
    !         a Propagate
    !         b Update boundary conditions
    !         c Diffusion correction
    !       ----------------------------------------
    !       4.  Store Q field in spectra
    !     ---------------------------------------------
    !
    !  9. Switches :
    !
    !       !/MGP   Correct for motion of grid.
    !
    !       !/TDYN  Dynamic increase of DTME
    !       !/DSS0  Disable diffusion in propagation direction
    !       !/XW0   Propagation diffusion only.
    !       !/XW1   Growth diffusion only.
    !
    !       !/S     Enable subroutine tracing.
    !
    !       !/T     Enable general test output.
    !       !/T1    Dump of input field and fluxes.
    !       !/T2    Dump of output field.
    !
    ! 10. Source code :
    !
    !/ ------------------------------------------------------------------- /
    USE constants
    !
    USE w3timemd, ONLY: dsec21
    !
    USE w3gdatmd, ONLY: nk, nth, dth, xfr, esin, ecos, sig, nx, ny,  &
         nsea, sx, sy, mapsf, funo3, fverg,           &
         ijkcel, ijkufc, ijkvfc, ncel, nufc, nvfc,    &
         ijkcel3, ijkcel4,                            &
         ijkvfc5, ijkvfc6,ijkufc5,ijkufc6,            &
         nlvcel, nlvufc, nlvvfc, nrlv, mrfct,         &
         dtcfl, clats, dtms, ctrnx, ctrny
    USE w3gdatmd, ONLY: nglo, angarc, arctc
    USE w3wdatmd, ONLY: time
    USE w3adatmd, ONLY: cg, wn, u10, cx, cy, atrnx, atrny, itime
    !
    USE w3idatmd, ONLY: flcur
    USE w3odatmd, ONLY: ndse, ndst, flbpi, nbi, tbpi0, tbpin,       &
         isbpi, bbpi0, bbpin
    !
#ifdef W3_S
    USE w3servmd, ONLY: strace
#endif
    !/
    IMPLICIT NONE
    !/
    !/ ------------------------------------------------------------------- /
    !/ Parameter list
    !/
    INTEGER, INTENT(IN)     :: ISP
    REAL,    INTENT(IN)     :: DTG
    REAL,    INTENT(INOUT)  :: VQ(NSEA)
    !/
    !/ ------------------------------------------------------------------- /
    !/ Local parameters
    !/
    INTEGER                 :: ITH, IK, NTLOC, ITLOC, ISEA, IXY,    &
         IY, IY0, IP, IBI, LvR
    INTEGER                 :: i, j, k, L, M, N, LL, MM, NN, LMN,   &
         iuf, juf, ivf, jvf, icl, jcl
#ifdef W3_S
    INTEGER, SAVE           :: IENT = 0
#endif
    REAL                    :: CG0, CGA, CGN, CGX, CGY, FMR, RD1,   &
         RD2, CXMIN, CXMAX, CYMIN, CYMAX,     &
         CXC, CYC, DTLDX, DTLDY
    REAL                    :: DTLOC, CGCOS, CGSIN, FUTRN, FVTRN,   &
         DFRR, DX0I, DY0I, CGD, DSSD,         &
         DNND, DCELL, XWIND, TFAC, DSS, DNN
    REAL                    :: PCArea, ARCTH
    LOGICAL                 :: YFIRST
    !/
    !/ Automatic work arrays
    !
    REAL, Dimension(-9:NCel) ::  FCNt, AFCN, BCNt, UCFL, VCFL, CQ,  &
         CQA, CXTOT, CYTOT
    REAL, Dimension(   NUFc) ::  FUMD, FUDIFX, ULCFLX
    REAL, Dimension(   NVFc) ::  FVMD, FVDIFY, VLCFLY
    !/
    !/ ------------------------------------------------------------------- /
    !/
#ifdef W3_S
    CALL strace (ient, 'W3PSMC')
#endif
    !
    ! 1.  Preparations --------------------------------------------------- *
 
    !!Li  Spectral bin direction and frequency indices
    ith    = 1 + mod(isp-1,nth)
    ik     = 1 + (isp-1)/nth
 
    !!Li  Maximum group speed for 1st and the transported frequency bin
    !!Li  A factor of 1.2 is added to account for the shallow water peak.
    cg0    = 0.6 * grav / sig(1)
    cga    = 0.6 * grav / sig(ik)
 
    !!Li  Maximum group speed for given spectral bin. First bin speed is
    !!Li  used to avoid zero speed component.
    !     CGX    = ABS( CGA * ECOS(ITH) )
    !     CGY    = ABS( CGA * ESIN(ITH) )
    cgx    =  cga * ecos(ith)
    cgy    =  cga * esin(ith)
 
    !!Li  Add maximum current components to maximum group components.
    IF ( flcur ) THEN
      cxmin  = minval( cx(1:nsea) )
      cxmax  = maxval( cx(1:nsea) )
      cymin  = minval( cy(1:nsea) )
      cymax  = maxval( cy(1:nsea) )
      IF ( abs(cgx+cxmin) .GT. abs(cgx+cxmax) ) THEN
        cgx    = cgx + cxmin
      ELSE
        cgx    = cgx + cxmax
      END IF
      IF ( abs(cgy+cymin) .GT. abs(cgy+cymax) ) THEN
        cgy    = cgy + cymin
      ELSE
        cgy    = cgy + cymax
      END IF
      cxc    = max( abs(cxmin) , abs(cxmax) )
      cyc    = max( abs(cymin) , abs(cymax) )
    ELSE
      cxc    = 0.
      cyc    = 0.
    END IF
 
    !!Li  Base-cell grid lenth at Equator (size-4 on SMC625 grid).
    dx0i  = 1.0/(sx * dera * radius)
    dy0i  = 1.0/(sy * dera * radius)
 
    !!Li  Miminum time step determined by Courant number < 0.8
    !!Li  Note, minimum x grid length is half the Equator value.
    !!Li  Minimum time step should not be less than sub w3init requirement,
    !!Li  where IAPPRO array is initialised for propagation parallization.
    cgn   = 0.9999 * max( abs(cgx), abs(cgy), cxc, cyc, 0.001*cg0 )
    dtloc = dtcfl*cg0/cgn
    ntloc  = 1 + int(dtg/dtloc - 0.001)
    dtloc  = dtg / real(ntloc)
 
    !!Li  Group speed component common factors, FACX=DTG*DX0I
    !!Li  FACX and FACY are evaluated here directly.  JGLi16Jan2013
    !     CGCOS   = FACX * ECOS(ITH) / REAL(NTLOC)
    !     CGSIN   = FACY * ESIN(ITH) / REAL(NTLOC)
    cgcos   = ecos(ith)
    cgsin   = esin(ith)
    dtldx   = dtloc * dx0i
    dtldy   = dtloc * dy0i
    !
    yfirst = mod(itime,2) .EQ. 0
    !
    !Li   Homogenous diffusion Fourier number DNND and DSSD will be used.
    !Li   They have to be divided by base-cell size for size-1 stability.
    !Li   So they are equivalent to the Fourier number in size-1 cell at
    !Li   the sub-time step DTLOC/MRFct.
    IF ( dtms .GT. 0. ) THEN
      dfrr   = xfr - 1.
      cgd    = 0.5 * grav / sig(ik)
      dnn    = ((dth*cgd)**2)*dtms / 12.
      dnnd   = dnn*dtloc*(dx0i*dx0i)
      dssd   = dnn*dtloc*(dy0i*dy0i)
    ELSE
      dssd = 0.0
      dnnd = 0.0
    END IF
    !
    ! 1.b Initialize arrays
    !
#ifdef W3_T
    WRITE (ndst,9010)
#endif
    !
    ulcflx = 0.
    vlcfly = 0.
 
    !Li    Pass spectral element VQ to CQ and define size-1 cell CFL
#ifdef W3_OMPG
    !$OMP Parallel DO Private(ISEA)
#endif
    DO isea=1, nsea
      !Li  Transported variable is divided by CG as in WW3.
      cq(isea) = vq(isea)/cg(ik,isea)
      !Li  Resetting NaNQ VQ to zero if any.   JGLi18Mar2013
      IF( .NOT. (cq(isea) .EQ. cq(isea)) )  cq(isea) = 0.0
    END DO
#ifdef W3_OMPG
    !$OMP END Parallel DO
#endif
 
    !Li  Add current components if any to wave velocity.
    IF ( flcur ) THEN
#ifdef W3_OMPG
      !$OMP Parallel DO Private(ISEA)
#endif
      DO isea=1, nsea
        cxtot(isea) = (cgcos * cg(ik,isea) + cx(isea))
        cytot(isea) = (cgsin * cg(ik,isea) + cy(isea))
      ENDDO
#ifdef W3_OMPG
      !$OMP END Parallel DO
#endif
    ELSE
      !Li   No current case use group speed only.
#ifdef W3_OMPG
      !$OMP Parallel DO Private(ISEA)
#endif
      DO isea=1, nsea
        cxtot(isea) =  cgcos * cg(ik,isea)
        cytot(isea) =  cgsin * cg(ik,isea)
      END DO
#ifdef W3_OMPG
      !$OMP END Parallel DO
#endif
      !Li   End of IF( FLCUR ) block.
    ENDIF
 
    !Li   Arctic cell velocity components need to be rotated
    !Li   back to local east referenence system for propagation.
    IF( arctc ) THEN
      DO isea=nglo+1, nsea
        arcth = angarc(isea-nglo)*dera
        cxc = cxtot(isea)*cos(arcth) + cytot(isea)*sin(arcth)
        cyc = cytot(isea)*cos(arcth) - cxtot(isea)*sin(arcth)
        cxtot(isea) = cxc
        cytot(isea) = cyc
      END DO
      !Li   Polar cell area factor for V-flux update
      pcarea = dy0i/(mrfct*pi*dx0i*float(ijkcel(4,nsea)))
      !Li   V-component is reset to zero for Polar cell as direction
      !Li   is undefined there.
      cytot(nsea) = 0.0
    ENDIF
 
 
    !Li     Convert velocity components into CFL factors.
#ifdef W3_OMPG
    !$OMP Parallel DO Private(ISEA)
#endif
    DO isea=1, nsea
      ucfl(isea) = dtldx*cxtot(isea)/clats(isea)
      vcfl(isea) = dtldy*cytot(isea)
    ENDDO
#ifdef W3_OMPG
    !$OMP END Parallel DO
#endif
 
    !Li  Initialise boundary cell CQ and Velocity values.
    cq(-9:0)=0.0
    ucfl(-9:0)=0.0
    vcfl(-9:0)=0.0
    !
    ! 3.  Loop over frequency-dependent sub-steps -------------------------*
    !
    DO itloc=1, ntloc
      !
      !     Initialise net flux arrays.
      fcnt(-9:ncel) = 0.0
      afcn(-9:ncel) = 0.0
      bcnt(-9:ncel) = 0.0
      !
      !     Single-resolution SMC grid uses regular grid advection with
      !     partial blocking enabled when NRLv = 1
      IF ( nrlv .EQ. 1 ) THEN
        IF( funo3 ) THEN
          !  Use 3rd order UNO3 scheme.  JGLi20Aug2015
          CALL smcxuno3r(1, nufc, cq, ucfl, ulcflx, dnnd, fumd, fudifx)
        ELSE
          !  Call SMCxUNO2 to calculate MFx value
          CALL smcxuno2r(1, nufc, cq, ucfl, ulcflx, dnnd, fumd, fudifx)
        ENDIF
 
        !  Store conservative flux in FCNt advective one in AFCN
#ifdef W3_OMPG
        !$OMP Parallel DO Private(i, M, N, FUTRN)
#endif
        DO i=1, nufc
          m=ijkufc5(i)
          n=ijkufc6(i)
          futrn = fumd(i)*ulcflx(i) - fudifx(i)
 
          !! Add sub-grid transparency for input flux update.  JGLi16May2011
          !! Transparency is also applied on diffusion flux.   JGLi12Mar2012
          !! Replace CRITICAL with ATOMIC.  JGLi15Jan2019
          !! !$OMP CRITICAL
          !! Remove boundary cell flux update or M N > 0.  JGLi28Mar2019
          IF( m > 0 ) THEN
            IF( (ctrnx(m)+ctrnx(n)) .GE. 1.96 )  THEN
#ifdef W3_OMPG
              !$OMP ATOMIC
#endif
              fcnt(m) = fcnt(m) - futrn
            ELSE IF( ulcflx(i) .GE. 0.0 )  THEN
#ifdef W3_OMPG
              !$OMP ATOMIC
#endif
              fcnt(m) = fcnt(m) - futrn*ctrnx(m)
            ELSE
#ifdef W3_OMPG
              !$OMP ATOMIC
#endif
              fcnt(m) = fcnt(m) - futrn*ctrnx(n)*ctrnx(m)
            ENDIF
            !  Also divided by another cell length as UCFL is in basic unit.
#ifdef W3_OMPG
            !$OMP ATOMIC
#endif
            ! ChrisB: Re-arranged the RHS term below to make it
            ! valid for OMP ATMOIC directive.
            afcn(m) = afcn(m) - (fumd(i)*ucfl(m) - fudifx(i))
          ENDIF
 
          IF( n > 0 ) THEN
            IF( (ctrnx(m)+ctrnx(n)) .GE. 1.96 )  THEN
#ifdef W3_OMPG
              !$OMP ATOMIC
#endif
              fcnt(n) = fcnt(n) + futrn
            ELSE IF( ulcflx(i) .GE. 0.0 )  THEN
#ifdef W3_OMPG
              !$OMP ATOMIC
#endif
              fcnt(n) = fcnt(n) + futrn*ctrnx(m)*ctrnx(n)
            ELSE
#ifdef W3_OMPG
              !$OMP ATOMIC
#endif
              fcnt(n) = fcnt(n) + futrn*ctrnx(n)
            ENDIF
            !  Also divided by another cell length as UCFL is in basic unit.
#ifdef W3_OMPG
            !$OMP ATOMIC
#endif
            afcn(n) = afcn(n) + (fumd(i)*ucfl(n) - fudifx(i))
          ENDIF
          !! !$OMP END CRITICAL
 
        ENDDO
#ifdef W3_OMPG
        !$OMP END Parallel DO
#endif
 
        !  Store conservative update in CQA and advective update in CQ
        !  The side length in MF value has to be cancelled with cell length
        !  Note ULCFLX has been divided by the cell size inside SMCxUNO2.
#ifdef W3_OMPG
        !$OMP Parallel DO Private(n)
#endif
        DO n=1, nsea
          cqa(n)=cq(n) + fcnt(n)/float(ijkcel3(n))
          cq(n)=cq(n) + afcn(n)/float(ijkcel3(n))
        ENDDO
#ifdef W3_OMPG
        !$OMP END Parallel DO
#endif
 
        !  Call advection subs.
        IF( funo3 ) THEN
          !  Use 3rd order UNO3 scheme.  JGLi20Aug2015
          CALL smcyuno3r(1, nvfc, cq, vcfl, vlcfly, dssd, fvmd, fvdify)
        ELSE
          !  Call SMCyUNO2 to calculate MFy value
          CALL smcyuno2r(1, nvfc, cq, vcfl, vlcfly, dssd, fvmd, fvdify)
        ENDIF
 
#ifdef W3_OMPG
        !$OMP Parallel DO Private(j, M, N, FVTRN)
#endif
        DO j=1, nvfc
          m=ijkvfc5(j)
          n=ijkvfc6(j)
          fvtrn = fvmd(j)*vlcfly(j) - fvdify(j)
 
          !! Add sub-grid transparency for input flux update.  JGLi16May2011
          !! Transparency is also applied on diffusion flux.   JGLi12Mar2012
          !! Replace CRITICAL with ATOMIC.  JGLi15Jan2019
          !! !$OMP CRITICAL
          !! Remove boundary cell flux update or M N > 0.  JGLi28Mar2019
          IF( m > 0 ) THEN
            IF( (ctrny(m)+ctrny(n)) .GE. 1.96 )  THEN
#ifdef W3_OMPG
              !$OMP ATOMIC
#endif
              bcnt(m) = bcnt(m) - fvtrn
            ELSE IF( vlcfly(j) .GE. 0.0 )  THEN
#ifdef W3_OMPG
              !$OMP ATOMIC
#endif
              bcnt(m) = bcnt(m) - fvtrn*ctrny(m)
            ELSE
#ifdef W3_OMPG
              !$OMP ATOMIC
#endif
              bcnt(m) = bcnt(m) - fvtrn*ctrny(n)*ctrny(m)
            ENDIF
          ENDIF
          IF( n > 0 ) THEN
            IF( (ctrny(m)+ctrny(n)) .GE. 1.96 )  THEN
#ifdef W3_OMPG
              !$OMP ATOMIC
#endif
              bcnt(n) = bcnt(n) + fvtrn
            ELSE IF( vlcfly(j) .GE. 0.0 )  THEN
#ifdef W3_OMPG
              !$OMP ATOMIC
#endif
              bcnt(n) = bcnt(n) + fvtrn*ctrny(m)*ctrny(n)
            ELSE
#ifdef W3_OMPG
              !$OMP ATOMIC
#endif
              bcnt(n) = bcnt(n) + fvtrn*ctrny(n)
            ENDIF
          ENDIF
          !! !$OMP END CRITICAL
        ENDDO
#ifdef W3_OMPG
        !$OMP END Parallel DO
#endif
 
        !  Store conservative update of CQA in CQ
        !  The v side length in MF value has to be cancelled with cell length
        !! One cosine factor is also needed to be divided for SMC grid
#ifdef W3_OMPG
        !$OMP Parallel DO Private(n)
#endif
        DO n=1, nsea
          cq(n)=cqa(n) + bcnt(n)/( clats(n)*float(ijkcel3(n)) )
        ENDDO
#ifdef W3_OMPG
        !$OMP END Parallel DO
#endif
        !   Polar cell needs a special area factor, one-level case.
        IF( arctc ) cq(nsea) = cqa(nsea) + bcnt(nsea)*pcarea
 
        !   End of single-resolution advection and diffusion.
      ELSE
 
        !     Multi-resolution SMC grid uses irregular grid advection scheme
        !     without partial blocking when NRLv > 1
        !
        ! 3.a    Multiresolution sub-steps depend on refined levels MRFct
        DO lmn = 1, mrfct
          !
          ! 3.b    Loop over all levels, starting from the finest level.
          DO ll= 1, nrlv
 
            !        Cell size of this level
            lvr=2**(ll - 1)
            fmr=float( lvr )
            !
            ! 3.c    Calculate this level only if size is factor of LMN
            IF( mod(lmn, lvr) .EQ. 0 ) THEN
              !
              ! 3.d    Select cell and face ranges
              icl=nlvcel(ll-1)+1
              iuf=nlvufc(ll-1)+1
              ivf=nlvvfc(ll-1)+1
              jcl=nlvcel(ll)
              juf=nlvufc(ll)
              jvf=nlvvfc(ll)
              !
              !  Use 3rd order UNO3 scheme.  JGLi03Sep2015
              IF( funo3 ) THEN
                CALL smcxuno3(iuf, juf, cq, ucfl, ulcflx, dnnd, fumd, fudifx, fmr)
              ELSE
                !  Call SMCxUNO2 to calculate finest level (size-1) MFx value
                CALL smcxuno2(iuf, juf, cq, ucfl, ulcflx, dnnd, fumd, fudifx, fmr)
              ENDIF
 
              !  Store fineset level conservative flux in FCNt advective one in AFCN
#ifdef W3_OMPG
              !$OMP Parallel DO Private(i, L, M, FUTRN)
#endif
              DO i=iuf, juf
                l=ijkufc5(i)
                m=ijkufc6(i)
                futrn = fumd(i)*ulcflx(i) - fudifx(i)
                !! Replace CRITICAL with ATOMIC.  JGLi15Jan2019
                !! !$OMP CRITICAL
                !! Remove boundary cell flux update or L M > 0.  JGLi28Mar2019
                IF( l > 0 ) THEN
                  !! Add sub-grid blocking for refined cells.   JGLi18Apr2018
                  IF( (ctrnx(m)+ctrnx(l)) .GE. 1.96 )  THEN
#ifdef W3_OMPG
                    !$OMP ATOMIC
#endif
                    fcnt(l) = fcnt(l) - futrn
                  ELSE IF( ulcflx(i) .GE. 0.0 ) THEN
#ifdef W3_OMPG
                    !$OMP ATOMIC
#endif
                    fcnt(l) = fcnt(l) - futrn*ctrnx(l)
                  ELSE
#ifdef W3_OMPG
                    !$OMP ATOMIC
#endif
                    fcnt(l) = fcnt(l) - futrn*ctrnx(l)*ctrnx(m)
                  ENDIF
#ifdef W3_OMPG
                  !$OMP ATOMIC
#endif
                  ! ChrisB: Re-arranged the RHS term below to make it
                  ! valid for OMP ATMOIC directive.
                  afcn(l) = afcn(l) - (fumd(i)*ucfl(l)*fmr - fudifx(i))
                ENDIF
                IF( m > 0 ) THEN
                  !! Add sub-grid blocking for refined cells.   JGLi18Apr2018
                  IF( (ctrnx(m)+ctrnx(l)) .GE. 1.96 )  THEN
#ifdef W3_OMPG
                    !$OMP ATOMIC
#endif
                    fcnt(m) = fcnt(m) + futrn
                  ELSE IF( ulcflx(i) .GE. 0.0 ) THEN
#ifdef W3_OMPG
                    !$OMP ATOMIC
#endif
                    fcnt(m) = fcnt(m) + futrn*ctrnx(m)*ctrnx(l)
                  ELSE
#ifdef W3_OMPG
                    !$OMP ATOMIC
#endif
                    fcnt(m) = fcnt(m) + futrn*ctrnx(m)
                  ENDIF
#ifdef W3_OMPG
                  !$OMP ATOMIC
#endif
                  afcn(m) = afcn(m) + (fumd(i)*ucfl(m)*fmr - fudifx(i))
                ENDIF
                !! !$OMP END CRITICAL
              ENDDO
#ifdef W3_OMPG
              !$OMP END Parallel DO
#endif
 
              !  Store conservative update in CQA and advective update in CQ
              !  The side length in MF value has to be cancelled with cell y-length.
              !  Also divided by another cell x-size as UCFL is in size-1 unit.
#ifdef W3_OMPG
              !$OMP Parallel DO Private(n)
#endif
              DO n=icl, jcl
                cqa(n)=cq(n) + fcnt(n)/float( ijkcel3(n)*ijkcel4(n) )
                cq(n)=cq(n) + afcn(n)/float( ijkcel3(n)*ijkcel4(n) )
                fcnt(n)=0.0
                afcn(n)=0.0
              ENDDO
#ifdef W3_OMPG
              !$OMP END Parallel DO
#endif
              !
              !  Use 3rd order UNO3 scheme.  JGLi03Sep2015
              IF( funo3 ) THEN
                CALL smcyuno3(ivf, jvf, cq, vcfl, vlcfly, dssd, fvmd, fvdify, fmr)
              ELSE
                !  Call SMCyUNO2 to calculate MFy value
                CALL smcyuno2(ivf, jvf, cq, vcfl, vlcfly, dssd, fvmd, fvdify, fmr)
              ENDIF
              !
              !  Store conservative flux in BCNt
#ifdef W3_OMPG
              !$OMP Parallel DO Private(j, L, M, FVTRN)
#endif
              DO j=ivf, jvf
                l=ijkvfc5(j)
                m=ijkvfc6(j)
                fvtrn = fvmd(j)*vlcfly(j) - fvdify(j)
                !! Replace CRITICAL with ATOMIC.  JGLi15Jan2019
                !! !$OMP CRITICAL
                !! Remove boundary cell flux update or L M > 0.  JGLi28Mar2019
                IF( l > 0 ) THEN
                  !! Add sub-grid blocking for refined cells.   JGLi18Apr2018
                  IF( (ctrny(m)+ctrny(l)) .GE. 1.96 )  THEN
#ifdef W3_OMPG
                    !$OMP ATOMIC
#endif
                    bcnt(l) = bcnt(l) - fvtrn
                  ELSE IF( vlcfly(j) .GE. 0.0 )  THEN
#ifdef W3_OMPG
                    !$OMP ATOMIC
#endif
                    bcnt(l) = bcnt(l) - fvtrn*ctrny(l)
                  ELSE
#ifdef W3_OMPG
                    !$OMP ATOMIC
#endif
                    bcnt(l) = bcnt(l) - fvtrn*ctrny(l)*ctrny(m)
                  ENDIF
                ENDIF
                IF( m > 0 ) THEN
                  !! Add sub-grid blocking for refined cells.   JGLi18Apr2018
                  IF( (ctrny(m)+ctrny(l)) .GE. 1.96 )  THEN
#ifdef W3_OMPG
                    !$OMP ATOMIC
#endif
                    bcnt(m) = bcnt(m) + fvtrn
                  ELSE IF( vlcfly(j) .GE. 0.0 )  THEN
#ifdef W3_OMPG
                    !$OMP ATOMIC
#endif
                    bcnt(m) = bcnt(m) + fvtrn*ctrny(m)*ctrny(l)
                  ELSE
#ifdef W3_OMPG
                    !$OMP ATOMIC
#endif
                    bcnt(m) = bcnt(m) + fvtrn*ctrny(m)
                  ENDIF
                ENDIF
                !! !$OMP END CRITICAL
              ENDDO
#ifdef W3_OMPG
              !$OMP END Parallel DO
#endif
 
              !  Store conservative update of CQA in CQ
              !  The v side length in MF value has to be cancelled with x-size.
              !  Also divided by cell y-size as VCFL is in size-1 unit.
              !! One cosine factor is also needed to be divided for SMC grid.
#ifdef W3_OMPG
              !$OMP Parallel DO Private(n)
#endif
              DO n=icl, jcl
                cq(n)=cqa(n) + bcnt(n)/( clats(n)*            &
                     &             float( ijkcel3(n)*ijkcel4(n) ) )
                bcnt(n)=0.0
              ENDDO
#ifdef W3_OMPG
              !$OMP END Parallel DO
#endif
              !Li  Polar cell needs a special area factor, multi-level case.
              IF( arctc .AND. jcl .EQ. nsea ) THEN
                cq(nsea) = cqa(nsea) + bcnt(nsea)*pcarea
              ENDIF
              !
              !  End of refine level if block  MOD(LMN, LvR) .EQ. 0
            ENDIF
 
            !  End of refine level loop LL=1, NRLv
          ENDDO
          !!
          !!    END of multi-resolution sub-step loop LMN = 1, MRFct
        ENDDO
 
        !  End of multi-resolution advection ELSE block of NRLv > 1
      ENDIF
 
      !!   Update boundary spectra if any.  JGLi26Feb2016
      !
      IF ( flbpi ) THEN
        rd1 = dsec21(tbpi0, time)-dtg*real(ntloc-itloc)/real(ntloc)
        rd2 = dsec21(tbpi0, tbpin)
        IF ( rd2 .GT. 0.001 ) THEN
          rd2    = min(1.,max(0.,rd1/rd2))
          rd1    = 1. - rd2
        ELSE
          rd1    = 0.
          rd2    = 1.
        END IF
        DO ibi=1, nbi
          isea     = isbpi(ibi)
          cq(isea) = (rd1*bbpi0(isp,ibi) + rd2*bbpin(isp,ibi))   &
               /cg(ik,isea)
        END DO
      ENDIF
      !
      !!    End of ITLOC DO
    ENDDO
 
    !  Average with 1-2-1 scheme.  JGLi20Aug2015
    IF(fverg) CALL smcaverg(cq)
 
    !
    ! 4.  Store results in VQ in proper format --------------------------- *
    !
#ifdef W3_OMPG
    !$OMP Parallel DO Private(ISEA)
#endif
    DO isea=1, nsea
      vq(isea) =  max( 0. , cq(isea)*cg(ik,isea) )
    END DO
#ifdef W3_OMPG
    !$OMP END Parallel DO
#endif
    !
    RETURN
    !
    ! Formats
    !
#ifdef W3_T
9001 FORMAT (' TEST W3PSMC : ISP, ITH, IK, COS-SIN :',i8,2i4,2f7.3)
9003 FORMAT (' TEST W3PSMC : NO DISPERSION CORRECTION ')
9010 FORMAT (' TEST W3PSMC : INITIALIZE ARRAYS')
9020 FORMAT (' TEST W3PSMC : CALCULATING LCFLX/Y AND DSS/NN (NSEA=', &
         i6,')')
#endif
#ifdef W3_T1
9021 FORMAT (1x,i6,2i5,e12.4,2f7.3)
#endif
#ifdef W3_T
9022 FORMAT (' TEST W3PSMC : CORRECTING FOR CURRENT')
9040 FORMAT (' TEST W3PSMC : FIELD AFTER PROP. (NSEA=',i6,')')
#endif
#ifdef W3_T2
9041 FORMAT (1x,i6,2i5,e12.4)
#endif
    !/
    !/ End of W3PSMC ----------------------------------------------------- /
    !/
  END SUBROUTINE w3psmc
  !/
  !/ ------------------------------------------------------------------- /
  SUBROUTINE w3krtn ( ISEA, FACTH, FACK, CTHG0, CG, WN, DEPTH,    &
       DDDX, DDDY, ALFLMT, CX, CY, DCXDX, DCXDY,   &
       DCYDX, DCYDY, DCDX, DCDY, VA )
    !/
    !/                  +------------------------------------+
    !/                  | Spherical Multiple-Cell (SMC) grid |
    !/                  | Refraction and great-cirle turning |
    !/                  |           Jian-Guo Li              |
    !/                  | First created:     8 Nov 2010      |
    !/                  | Last modified:    06-Jun-2018      |
    !/                  +------------------------------------+
    !/
    !/    08-Nov-2010 : Origination.                        ( version 1.00 )
    !/    10-Jun-2011 : New refraction formulation.         ( version 1.10 )
    !/    16-Jun-2011 : Add refraction limiter to gradient. ( version 1.20 )
    !/    21-Jul-2011 : Old refraction formula + limiter.   ( version 1.30 )
    !/    26-Jul-2011 : Tidy up refraction schemes.         ( version 1.40 )
    !/    28-Jul-2011 : Finalise with old refraction.       ( version 1.50 )
    !/    23-Mar-2016 : Add current option in refraction.   ( version 2.30 )
    !/    06-Jun-2018 : Add DEBUGDCXDX                      ( version 6.04 )
    !/
    !/
    !  1. Purpose :
    !
    !     Refraction and great-circle turning by spectral rotation.
    !
    !  2. Method :
    !
    !     Linear interpolation equivalent to 1st order upstream scheme
    !     but without restriction on rotation angle.  However, refraction
    !     is limited towards the depth gradient direction (< 90 degree).
    !     Refraction induced spectral shift in the k-space will remain
    !     to be advected using the UNO2 scheme.
    !
    !  3. Parameters :
    !
    !     Parameter list
    !     ----------------------------------------------------------------
    !       ISEA    Int.   I   Number of sea point.
    !       FACTH/K Real   I   Factor in propagation velocity.
    !       CTHG0   Real   I   Factor in great circle refracftion term.
    !       MAPxx2  I.A.   I   Propagation and storage maps.
    !       CG      R.A.   I   Local group velocities.
    !       WN      R.A.   I   Local wavenumbers.
    !       DEPTH   R.A.   I   Depth.
    !       DDDx    Real   I   Depth gradients.
    !       CX/Y    Real   I   Current components.
    !       DCxDx   Real   I   Current gradients.
    !       DCDx    Real   I   Phase speed gradients.
    !       VA      R.A.  I/O  Spectrum.
    !     ----------------------------------------------------------------
    !
    !     Local variables.
    !     ----------------------------------------------------------------
    !       DPH2K   R.A.  2*Depth*Wave_number_K
    !       SNH2K   R.A.  SINH(2*Depth*Wave_number_K)
    !       FDD, FDU, FGC, FCD, FCU
    !               R.A.  Directionally varying part of depth, current and
    !                     great-circle refraction terms and of consit.
    !                     of Ck term.
    !       CFLT-K  R.A.  Propagation velocities of local fluxes.
    !       DB      R.A.  Wavenumber band widths at cell centers.
    !       DM      R.A.  Wavenumber band widths between cell centers and
    !                     next cell center.
    !       Q       R.A.  Extracted spectrum
    !     ----------------------------------------------------------------
    !
    !  4. Subroutines used :
    !
    !       SMCGtCrfr Refraction and GCT rotation in theta.
    !       SMCkUNO2  Refraction shift in k-space by UNO2.
    !       STRACE    Service routine.
    !
    !  5. Called by :
    !
    !       W3WAVE   Wave model routine.
    !
    !  6. Error messages :
    !
    !       None.
    !
    !  8. Structure :
    !
    !     -----------------------------------------------------------------
    !       1.  Preparations
    !         a Initialize arrays
    !         b Set constants and counters
    !       2.  Point  preparations
    !         a Calculate SNH2K
    !         b Extract spectrum
    !       3.  Refraction velocities
    !         a Filter level depth reffraction.
    !         b Depth refratcion velocity.
    !         c Current refraction velocity.
    !       4.  Wavenumber shift velocities
    !         a Prepare directional arrays
    !         b Calcuate velocity.
    !       5.  Propagate.
    !       6.  Store results.
    !     -----------------------------------------------------------------
    !
    !  9. Switches :
    !
    !       !/S     Enable subroutine tracing.
    !       !/T     Enable general test output.
    !
    ! 10. Source code :
    !
    !/ ------------------------------------------------------------------- /
    USE constants
    USE w3gdatmd, ONLY: nk, nth, nspec, sig, dsip, ecos, esin, &
         ec2, esc, es2, flcth, flck, ctmax, dth
    USE w3adatmd, ONLY: itime
    USE w3idatmd, ONLY: flcur
    USE w3odatmd, ONLY: ndse, ndst
#ifdef W3_S
    USE w3servmd, ONLY: strace
#endif
    !/
    IMPLICIT NONE
    !/
    !/ ------------------------------------------------------------------- /
    !/ Parameter list
    !/
    INTEGER, INTENT(IN) :: ISEA
#ifdef W3_S
    INTEGER, SAVE       :: IENT = 0
#endif
    REAL, INTENT(IN)    :: FACTH, FACK, CTHG0, CG(0:NK+1),      &
         WN(0:NK+1), DEPTH, DDDX, DDDY,       &
         ALFLMT(NTH), CX, CY, DCXDX, DCXDY,   &
         DCYDX, DCYDY, DCDX(0:NK+1), DCDY(0:NK+1)
    REAL, INTENT(INOUT) :: VA(NSPEC)
    !/
    !/ ------------------------------------------------------------------- /
    !/ Local parameters
    !/
    INTEGER    :: ITH, IK, ISP
    REAL       :: FGC, FKD, FKS, FRK(NK), FRG(NK), DDNorm(NTH),     &
         FKC(NTH), VQ(NSPEC), VCFLT(NSPEC), DEPTH30,       &
         DB(0:NK+1), DM(-1:NK+1), CFLK(NTH,0:NK),          &
         !Li   For new refraction scheme using Cg.  JGLi26Jul2011
         !                   DPH2K(0:NK+1), SNH2K(0:NK+1)
         !Li   For old refraction scheme using phase speed.  JGLi26Jul2011
         sigsnh(0:nk+1)
    !/
    !/ ------------------------------------------------------------------- /
    !/
#ifdef W3_S
    CALL strace (ient, 'W3KRTN')
#endif
    !
    ! 1.  Preparation for point ------------------------------------------ *
    !     Array with partial derivative of sigma versus depth
    !Li   Use of minimum depth 30 m for refraction factor.  JGLi12Feb2014
    depth30=max(30.0, depth)
 
    DO ik=0, nk+1
      !Li   For old refraction scheme using phase speed.  JGLi8Jun2011
      !        DPH2K(IK) = 2.0*WN(IK)*DEPTH
      !        SNH2K(IK) = SINH( DPH2K(IK) )
      !Li   For new refraction scheme using Cg.  JGLi3Jun2011
      !!        SIGSNH(IK) = SIG(IK)/SINH(2.0*WN(IK)*DEPTH)
      !!AC  Replacing SIGSINH with a delimiter to prevent the SINH value from
      !!AC  becoming significantly large. Right now set to a max around 1E21
      !         SIGSNH(IK) = SIG(IK)/SINH(MIN(2.0*WN(IK)*DEPTH,50.0))
      !Li   Refraction factor uses minimum depth of 30 m.  JGLi12Feb2014
      sigsnh(ik) = sig(ik)/sinh(min(2.0*wn(ik)*depth30,50.0))
    END DO
    !
    ! 2.  Extract spectrum without mapping
    !
    vq = va
 
    ! 3.  Refraction velocities ------------------------------------------ *
    !
    !
    IF ( flcth ) THEN
      !
      ! 3.a Set slope filter for depth refraction
      !
      !Li   Lift theta-refraction limit and use new formulation.   25 Nov 2010
      !Li   FACTH  = DTG / DTH / REAL(NTLOC), CTHG0 = - TAN(DERA*Y) / RADIUS
      fgc    = facth * cthg0
      !
      DO ik=1, nk
        !Li   New refraction formulation using Cg only.  JGLi3Jun2011
        !         FRK(IK) = FACTH*2.0*SIG(IK)*(1.-DEPTH*SIG(IK)*SIG(IK)/GRAV)  &
        !     &                                /(DPH2K(IK)+SNH2K(IK))
        !Li   Old refraction formulation using phase speed.  JGLi8Jun2011
        frk(ik) = facth * sigsnh(ik)
        !Li
        frg(ik) = fgc * cg(ik)
      END DO
      !
      !Li   Current induced refraction stored in FKC.    JGLi30Mar2016
      IF ( flcur ) THEN
        DO ith=1, nth
          !Li   Put a CTMAX limit on current theta rotation.  JGLi02Mar2017
          !               FKC(ITH) = FACTH*( DCYDX*ES2(ITH) - DCXDY*EC2(ITH) +  &
          fgc      = facth*( dcydx*es2(ith) - dcxdy*ec2(ith) +  &
               (dcxdx - dcydy)*esc(ith) )
          fkc(ith) = sign( min(abs(fgc), ctmax), fgc )
        END DO
      ELSE
        fkc(:)=0.0
      END IF
      !
      ! 3.b Depth refraction and great-circle turning.
      !
      DO ith=1, nth
        ddnorm(ith)=esin(ith)*dddx-ecos(ith)*dddy
        DO ik=1, nk
          isp = (ik-1)*nth + ith
          !Li   Apply depth gradient limited refraction, current and GCT term
          vcflt(isp)=frg(ik)*ecos(ith) + fkc(ith) +          &
               sign( min(abs(frk(ik)*ddnorm(ith)), alflmt(ith)),  &
               !Li   For new refraction scheme using Cg.  JGLi26Jul2011
               !                             FRK(IK)*DDNorm(ITH) )
               !Li   For old refraction scheme using phase speed.  JGLi26Jul2011
               ddnorm(ith) )
        END DO
      END DO
 
    END IF
    !
    ! 4.  Wavenumber shift velocities due to current refraction ---------- *
    !
    IF ( flck ) THEN
      !
      ! 4.a Directionally dependent part
      !
      DO ith=1, nth
        !Li   Depth induced refraction is suspended as it is absorbed in
        !Li   the fixed frequency bin used for wave spectrum.  JGLi30Mar2016
        !           FKC(ITH) = ( ECOS(ITH)*DDDX + ESIN(ITH)*DDDY )
        fkc(ith) = -dcxdx*ec2(ith) -dcydy*es2(ith)               &
             -(dcxdy + dcydx)*esc(ith)
      END DO
      fkd = cx*dddx + cy*dddy
      !
      ! 4.b Band widths
      !
      !Li  Cell and side indices for k-dimension are arranged as
      !    Cell:    | -1 | 0 | 1 | 2 | ... | NK | NK+1 | NK+2 |
      !    Side:        -1   0   1   2 ...     NK     NK+1
      !Li  DSIP = SIG(K+1) - SIG(K), radian frequency increment
      !
      DO ik=0, nk
        db(ik) = dsip(ik) / cg(ik)
        dm(ik) = wn(ik+1) - wn(ik)
      END DO
      db(nk+1) = dsip(nk+1) / cg(nk+1)
      dm(nk+1) = dm(nk)
      dm(  -1) = dm( 0)
 
      ! 4.c Courant number of k-velocity without dividing by dk
      !!Li  FACK = DTG / REAL(NTLOC)
      !
      DO ik=0, nk
        !Li   For new refraction scheme using Cg.  JGLi3Jun2011
        !           FKS   = - FACK*WN(IK)*SIG(IK)/SNH2K(IK)
        !Li   Old refraction formulation using phase speed.  JGLi8Jun2011
        !           FKS   = - FACK*WN(IK)*SIGSNH(IK)
        !Li   Current induced k-shift.   JGLi30Mar2016
        fks = max( 0.0, cg(ik)*wn(ik)-0.5*sig(ik) )*fkd /    &
             ( depth30*cg(ik) )
        DO ith=1, nth
          cflk(ith,ik) = fack*( fks + fkc(ith)*wn(ik) )
        END DO
      END DO
      !Li   No CFL limiter is required here as it is applied in SMCkUNO2.
      !
    END IF
    !
    ! 5.  Propagate ------------------------------------------------------ *
    !
    IF ( mod(itime,2) .EQ. 0 ) THEN
      IF ( flck ) THEN
        !!Li  Refraction caused k-space shift.
        CALL smckuno2(cflk, vq, db, dm)
      END IF
      IF ( flcth ) THEN
        !!Li  GCT and refraction by rotation in theta direction.
        CALL smcgtcrfr(vcflt, vq)
      END IF
    ELSE
      IF ( flcth ) THEN
        !!Li  GCT and refraction by rotation in theta direction.
        CALL smcgtcrfr(vcflt, vq)
      END IF
      IF ( flck )  THEN
        !!Li  Refraction caused k-space shift.
        CALL smckuno2(cflk, vq, db, dm)
      END IF
    END IF
    !
    ! 6.  Store reults --------------------------------------------------- *
    !
    va = vq
    !
    RETURN
    !
    !/ End of W3KRTN ----------------------------------------------------- /
    !/
  END SUBROUTINE w3krtn
 
 
  ! Subroutine that calculate mid-flux values for x dimension
  SUBROUTINE smcxuno2(NUA, NUB, CF, UC, UFLX, AKDif, FU, FX, FTS)
 
    USE constants
    USE w3gdatmd, ONLY: ncel, mrfct, nufc, ijkcel, ijkufc, clats, &
         ijkcel3, ijkcel4
    USE w3odatmd, ONLY: ndse, ndst
 
    IMPLICIT NONE
    INTEGER, INTENT( IN):: NUA, NUB
    REAL,    INTENT( IN):: CF(-9:NCel), UC(-9:NCel), AKDif, FTS
    REAL,    INTENT(Out):: UFLX(NUFc), FU(NUFc), FX(NUFc)
    !
    INTEGER ::  i, j, k, L, M, N, ij
    REAL:: CNST, CNST0, CNST1, CNST2, CNST3, CNST4, CNST5, CNST6, CNST8, CNST9
 
    !     Two layer of boundary cells are added to each boundary cell face
    !     with all boundary cell values CF(-9:0)=0.0.
 
    !     Diffusion Fourier no. at sub-time-step, proportional to face size,
    !     which is also equal to the sub-time-step factor FTS.
    cnst0=akdif*fts*fts
    !     Uniform diffusion coefficient for all sizes.  JGLi24Feb2012
    !        CNST0=AKDif*MRFct*FTS
 
#ifdef W3_OMPG
    !$OMP Parallel Default(Shared), Private(i, ij, K, L, M, N),  &
    !$OMP& Private(CNST,CNST1,CNST2,CNST3,CNST4,CNST5,CNST6,CNST8,CNST9)
#endif
 
    !    Notice an extra side length L is multiplied to mid-flux to give correct
    !    proportion of flux into the cells.  This length will be removed by the
    !    cell length when the tracer concentration is updated.
 
#ifdef W3_OMPG
    !$OMP DO
#endif
 
    DO i=nua, nub
 
      !    Select Upstream, Central and Downstream cells
      k=ijkufc(4,i)
      l=ijkufc(5,i)
      m=ijkufc(6,i)
      n=ijkufc(7,i)
 
      !    Face bounding cell lengths and central gradient
      cnst2=float( ijkcel3(l) )
      cnst3=float( ijkcel3(m) )
      cnst5=(cf(m)-cf(l))/( cnst2 + cnst3 )
 
      !    Courant number in local size-1 cell, arithmetic average.
      cnst6=0.5*( uc(l)+uc(m) )*fts
      uflx(i) = cnst6
 
      !    Multi-resolution SMC grid requires flux multiplied by face factor.
      cnst8 = float( ijkufc(3,i) )
 
      !    Diffusion factor in local size-1 cell, plus the cosine factors.
      !    2.0 factor to cancel that in gradient CNST5.  JGLi08Mar2012
      !    The maximum cell number is used to avoid the boundary cell number
      !    in selection of the cosine factor.
      ij= max(l, m)
      cnst9 = 2.0/( clats( ij )*clats( ij ) )
 
      !    For positive velocity case
      IF(cnst6 >= 0.0)  THEN
 
        !    Use central cell velocity for boundary flux.  JGLi06Apr2011
        IF( m .LE. 0) uflx(i) = uc(l)*fts
 
        !    Upstream cell length and gradient, depending on UFLX sign.
        cnst1=float( ijkcel3(k) )
        cnst4=(cf(l)-cf(k))/( cnst2 + cnst1 )
 
        !    Use minimum gradient all region.
        cnst=sign(1.0, cnst5)*min( abs(cnst4), abs(cnst5) )
 
        !    Mid-flux value inside central cell
        fu(i)=(cf(l) + cnst*(cnst2 - uflx(i)))*cnst8
 
        !    For negative velocity case
      ELSE
 
        !    Use central cell velocity for boundary flux.  JGLi06Apr2011
        IF( l .LE. 0) uflx(i) = uc(m)*fts
 
        !    Upstream cell length and gradient, depending on UFLX sign.
        cnst1=float( ijkcel3(n) )
        cnst4=(cf(n)-cf(m))/( cnst1 + cnst3 )
 
        !    Use minimum gradient outside monotonic region.
        cnst=sign(1.0, cnst5)*min( abs(cnst4), abs(cnst5) )
 
        !    Mid-flux value inside central cell M
        fu(i)=(cf(m) - cnst*(cnst3+uflx(i)))*cnst8
 
      ENDIF
 
      !    Diffusion flux by face gradient x DT
      fx(i)=cnst0*cnst5*cnst8*cnst9
 
    END DO
 
#ifdef W3_OMPG
    !$OMP END DO
    !$OMP END Parallel
#endif
 
    RETURN
  END SUBROUTINE smcxuno2
 
 
  SUBROUTINE smcyuno2(NVA, NVB, CF, VC, VFLY, AKDif, FV, FY, FTS)
 
    USE constants
    USE w3gdatmd, ONLY: ncel, mrfct, nvfc, ijkcel, ijkvfc, clatf, ijkcel4
    USE w3odatmd, ONLY: ndse, ndst
 
    IMPLICIT NONE
    INTEGER, INTENT( IN):: NVA, NVB
    REAL,    INTENT( IN):: CF(-9:NCel), VC(-9:NCel), AKDif, FTS
    REAL,    INTENT(Out):: VFLY(NVFc), FV(NVFc), FY(NVFc)
    INTEGER ::  i, j, k, L, M, N, ij
    REAL:: CNST, CNST0, CNST1, CNST2, CNST3, CNST4, CNST5, CNST6, CNST8
 
    !     Notice an extra side length L is multiplied to mid-flux to give correct
    !     proportion of flux into the cells.  This length will be removed by the
    !     cell length when the tracer concentration is updated.
 
    !     Diffusion Fourier no. at sub-time-step, proportional to face size,
    !     which is also equal to the sub-time-step factor FTS.
    !        CNST0=AKDif*FTS*FTS
    !     2.0 factor to cancel that in gradient CNST5.  JGLi08Mar2012
    cnst0=akdif*fts*fts*2.0
    !     Uniform diffusion coefficient for all sizes.  JGLi24Feb2012
    !        CNST0=AKDif*MRFct*FTS
 
#ifdef W3_OMPG
    !$OMP Parallel Default(Shared), Private(j, K, L, M, N),  &
    !$OMP& Private(CNST,CNST1,CNST2,CNST3,CNST4,CNST5,CNST6,CNST8)
    !$OMP DO
#endif
 
    DO j=nva, nvb
 
      !    Select Upstream, Central and Downstream cells
      k=ijkvfc(4,j)
      l=ijkvfc(5,j)
      m=ijkvfc(6,j)
      n=ijkvfc(7,j)
 
      !    Face bounding cell lengths and gradient
      cnst2=float( ijkcel4(l) )
      cnst3=float( ijkcel4(m) )
      cnst5=(cf(m)-cf(l))/( cnst2 + cnst3 )
 
      !    Courant number in local size-1 cell unit
      !    Multiply by multi-resolution time step factor  FTS
      cnst6=0.5*( vc(l)+vc(m) )*fts
      vfly(j) = cnst6
 
      !    Face size integer and cosine factor.
      !    CLATF is defined on V-face for SMC grid.  JGLi28Feb2012
      cnst8=clatf(j)*float( ijkvfc(3,j) )
 
      !    For positive velocity case
      IF(cnst6 >= 0.0)  THEN
 
        !    Boundary cell y-size is set equal to central cell y-size
        !    as y-boundary cell sizes are not proportional to refined
        !    inner cells but constant of the base cell y-size, and
        !    Use central cell speed for face speed.  JGLi06Apr2011
        IF( m .LE. 0 ) THEN
          vfly(j) = vc(l)*fts
          cnst3   = cnst2
        ENDIF
 
        !    Upstream cell size and irregular grid gradient, depending on VFLY.
        cnst1=float( ijkcel4(k) )
        cnst4=(cf(l)-cf(k))/( cnst2 + cnst1 )
 
        !    Use minimum gradient outside monotonic region
        cnst=sign(1.0, cnst5)*min( abs(cnst4), abs(cnst5) )
 
        !    Mid-flux value multiplied by face width and cosine factor
        fv(j)=( cf(l) + cnst*(cnst2 - vfly(j)) )*cnst8
 
        !    For negative velocity case
      ELSE
 
        !    Set boundary cell y-size equal to central cell y-size and
        !    Use central cell speed for flux face speed.  JGLi06Apr2011
        IF( l .LE. 0 ) THEN
          vfly(j) = vc(m)*fts
          cnst2   = cnst3
        ENDIF
 
        !    Upstream cell size and gradient, depending on VFLY sign.
        !    Side gradients for central cell includs 0.5 factor.
        cnst1=float( ijkcel4(n) )
        cnst4=(cf(n)-cf(m))/( cnst1 + cnst3 )
 
        !    Use minimum gradient outside monotonic region
        cnst=sign(1.0, cnst5)*min( abs(cnst4), abs(cnst5) )
 
        !    Mid-flux value multiplied by face width and cosine factor
        fv(j)=( cf(m) - cnst*(cnst3 + vfly(j)) )*cnst8
 
      ENDIF
 
      !    Diffusion flux by face gradient x DT x face_width x cos(lat)
      !    Multiply by multi-resolution time step factor FTS
      fy(j)=cnst0*cnst5*cnst8
 
    END DO
 
#ifdef W3_OMPG
    !$OMP END DO
    !$OMP END Parallel
#endif
 
    RETURN
  END SUBROUTINE smcyuno2
 
 
  SUBROUTINE smcxuno2r(NUA, NUB, CF, UC, UFLX, AKDif, FU, FX)
 
    USE constants
    USE w3gdatmd, ONLY: nsea, ny, ncel, nufc, ijkcel, ijkufc, clats
    USE w3gdatmd, ONLY: ijkcel3
    USE w3odatmd, ONLY: ndse, ndst
 
    IMPLICIT NONE
    INTEGER, INTENT( IN):: NUA, NUB
    REAL,    INTENT( IN):: CF(-9:NCel), UC(-9:NCel), AKDif
    REAL,    INTENT(Out):: UFLX(NUFc), FU(NUFc), FX(NUFc)
    !
    INTEGER ::  i, j, k, L, M, N, ij
    REAL:: CNST, CNST0, CNST1, CNST2, CNST3, CNST4, CNST5, CNST6
 
    !      Two layer of boundary cells are added to each boundary cell face
    !      with all boundary cell values CF(-9:0)=0.0.
 
    !      Notice an extra side length L is multiplied to mid-flux to give correct
    !      proportion of flux into the cells.  This length will be removed by the
    !      cell length when the tracer concentration is updated.
 
#ifdef W3_OMPG
    !$OMP Parallel Default(Shared), Private(i, ij, K, L, M, N),  &
    !$OMP& Private(CNST,CNST0,CNST1,CNST2,CNST3,CNST4,CNST5,CNST6)
    !$OMP DO
#endif
 
    DO i=nua, nub
 
      !    Select Upstream, Central and Downstream cells
      k=ijkufc(4,i)
      l=ijkufc(5,i)
      m=ijkufc(6,i)
      n=ijkufc(7,i)
 
      !    Face bounding cell lengths and gradient
      cnst2=float( ijkcel3(l) )
      cnst3=float( ijkcel3(m) )
      cnst5=(cf(m)-cf(l))
 
      !    Averaged Courant number for base-level cell face
      cnst6= 0.5*( uc(l)+uc(m) )
      uflx(i) = cnst6
 
      !    Diffusion Fourier number in local cell size
      !    To avoid boundary cell number, use maximum of L and M.
      ij= max(l, m)
      cnst0 = 2.0/( clats(ij)*clats(ij) )
 
      !    For positive velocity case
      IF(cnst6 >= 0.0)  THEN
 
        !    Use central cell velocity for boundary flux.  JGLi06Apr2011
        IF( m .LE. 0) uflx(i) = uc(l)
 
        !    Side gradient for upstream cell as regular grid.
        cnst4=(cf(l)-cf(k))
 
        !    Use minimum gradient all region with 0.5 factor
        cnst=sign(0.5, cnst5)*min( abs(cnst4), abs(cnst5) )
 
        !    Mid-flux value inside central cell
        fu(i)=(cf(l) + cnst*(1.0-uflx(i)/cnst2))
 
        !    For negative velocity case
      ELSE
 
        !    Use central cell velocity for boundary flux.  JGLi06Apr2011
        IF( l .LE. 0) uflx(i) = uc(m)
 
        !    Side gradient for upstream cell, depneding on UFLX sign.
        cnst4=(cf(n)-cf(m))
 
        !    Use minimum gradient outside monotonic region, include 0.5 factor
        cnst=sign(0.5, cnst5)*min( abs(cnst4), abs(cnst5) )
 
        !    Mid-flux value inside central cell M
        fu(i)=(cf(m) - cnst*(1.0+uflx(i)/cnst3))
 
      ENDIF
 
      !    Diffusion flux by face gradient x DT
      fx(i)=akdif*cnst0*cnst5/(cnst2 + cnst3)
 
    END DO
 
#ifdef W3_OMPG
    !$OMP END DO
    !$OMP END Parallel
#endif
 
    RETURN
  END SUBROUTINE smcxuno2r
 
 
  SUBROUTINE smcyuno2r(NVA, NVB, CF, VC, VFLY, AKDif, FV, FY)
 
    USE constants
    USE w3gdatmd, ONLY: nsea, ny, ncel, nvfc, ijkcel, ijkvfc, clatf
    USE w3odatmd, ONLY: ndse, ndst
 
    IMPLICIT NONE
    INTEGER, INTENT( IN):: NVA, NVB
    REAL,    INTENT( IN):: CF(-9:NCel), VC(-9:NCel), AKDif
    REAL,    INTENT(Out):: VFLY(NVFc), FV(NVFc), FY(NVFc)
    INTEGER ::  i, j, k, L, M, N, ij
    REAL:: CNST, CNST0, CNST1, CNST2, CNST3, CNST4, CNST5, CNST6, CNST8
 
    !     Notice an extra side length L is multiplied to mid-flux to give correct
    !     proportion of flux into the cells.  This length will be removed by the
    !     cell length when the tracer concentration is updated.
 
#ifdef W3_OMPG
    !$OMP Parallel Default(Shared), Private(j, K, L, M, N),  &
    !$OMP& Private(CNST,CNST4,CNST5,CNST6,CNST8)
    !$OMP DO
#endif
 
    DO j=nva, nvb
 
      !    Select Upstream, Central and Downstream cells
      k=ijkvfc(4,j)
      l=ijkvfc(5,j)
      m=ijkvfc(6,j)
      n=ijkvfc(7,j)
 
      !    Central face gradient.
      cnst5=(cf(m)-cf(l))
 
      !    Courant number in basic cell unit as dy is constant
      cnst6=0.5*( vc(l)+vc(m) )
      vfly(j) = cnst6
 
      !    Face size integer and cosine factor
      !    CLATF is defined on V-face for SMC grid.  JGLi28Feb2012
      cnst8=clatf(j)*float( ijkvfc(3,j) )
 
      !    For positive velocity case
      IF(cnst6 >= 0.0)  THEN
 
        !    Use central cell speed for flux face speed.  JGLi06Apr2011
        IF( m .LE. 0 ) vfly(j) = vc(l)
 
        !    Upstream face gradient, depending on VFLY sign.
        cnst4=(cf(l)-cf(k))
 
        !    Use minimum gradient, including 0.5 factor and central sign.
        cnst=sign(0.5, cnst5)*min( abs(cnst4), abs(cnst5) )
 
        !    Mid-flux value multiplied by face width and cosine factor
        fv(j)=( cf(l) + cnst*(1.0-vfly(j)) )*cnst8
 
        !    For negative velocity case
      ELSE
 
        !    Use central cell speed for flux face speed.  JGLi06Apr2011
        IF( l .LE. 0 ) vfly(j) = vc(m)
 
        !    Side gradients for upstream face, depending on VFLY sign.
        cnst4=(cf(n)-cf(m))
 
        !    Use minimum gradient, including 0.5 factor and central sign.
        cnst=sign(0.5, cnst5)*min( abs(cnst4), abs(cnst5) )
 
        !    Mid-flux value multiplied by face width and cosine factor
        fv(j)=( cf(m) - cnst*(1.0+vfly(j)) )*cnst8
 
      ENDIF
 
      !    Diffusion flux by face gradient x DT x face_width x cos(lat)
      fy(j)=akdif*cnst5*cnst8
 
    END DO
 
#ifdef W3_OMPG
    !$OMP END DO
    !$OMP END Parallel
#endif
 
    RETURN
  END SUBROUTINE smcyuno2r
 
 
  SUBROUTINE smcxuno3(NUA, NUB, CF, UC, UFLX, AKDif, FU, FX, FTS)
 
    USE constants
    USE w3gdatmd, ONLY: ncel, mrfct, nufc, ijkcel, ijkufc, clats
    USE w3gdatmd, ONLY: ijkcel3
    USE w3odatmd, ONLY: ndse, ndst
 
    IMPLICIT NONE
    INTEGER, INTENT( IN):: NUA, NUB
    REAL,    INTENT( IN):: CF(-9:NCel), UC(-9:NCel), AKDif, FTS
    REAL,    INTENT(Out):: UFLX(NUFc), FU(NUFc), FX(NUFc)
    !
    INTEGER ::  i, j, k, L, M, N, ij
    REAL    :: CNST, CNST0, CNST1, CNST2, CNST3, CNST4, CNST5, CNST6,  &
         CNST7, CNST8, CNST9
 
    !     Two layer of boundary cells are added to each boundary cell face
    !     with all boundary cell values CF(-9:0)=0.0.
 
    !     Diffusion Fourier no. at sub-time-step, proportional to face size,
    !     which is also equal to the sub-time-step factor FTS.
    !         CNST0=AKDif*FTS*FTS
    !     2.0 factor to cancel that in gradient CNST5.  JGLi03Sep2015
    cnst0=akdif*fts*fts*2.0
 
    !     Notice an extra side length L is multiplied to mid-flux to give correct
    !     proportion of flux into the cells.  This length will be removed by the
    !     cell length when the tracer concentration is updated.
 
#ifdef W3_OMPG
    !$OMP Parallel Default(Shared), Private(i, ij, K, L, M, N),  &
    !$OMP& Private(CNST,CNST1,CNST2,CNST3,CNST4,CNST5,CNST6,CNST7,CNST8,CNST9)
    !$OMP DO
#endif
 
    DO i=nua, nub
 
      !    Select Upstream, Central and Downstream cells
      k=ijkufc(4,i)
      l=ijkufc(5,i)
      m=ijkufc(6,i)
      n=ijkufc(7,i)
 
      !    Face bounding cell lengths and central gradient
      cnst2=float( ijkcel3(l) )
      cnst3=float( ijkcel3(m) )
      cnst5=(cf(m)-cf(l))/( cnst2 + cnst3 )
 
      !    Courant number in local size-1 cell, arithmetic average.
      cnst6=0.5*( uc(l)+uc(m) )*fts
      uflx(i) = cnst6
 
      !    Multi-resolution SMC grid requires flux multiplied by face factor.
      cnst8 = float( ijkufc(3,i) )
 
      !    Diffusion factor in local size-1 cell, plus the cosine factors.
      !    2.0 factor to cancel that in gradient CNST5.  JGLi08Mar2012
      !    The maximum cell number is used to avoid the boundary cell number
      !    in selection of the cosine factor.
      ij= max(l, m)
 
      !    For positive velocity case
      IF(cnst6 >= 0.0)  THEN
 
        !    Use central cell velocity for boundary flux.  JGLi06Apr2011
        IF( m .LE. 0) uflx(i) = uc(l)*fts
 
        !    Upstream cell length and gradient, depending on UFLX sign.
        cnst1=float( ijkcel3(k) )
        cnst4=(cf(l)-cf(k))/( cnst2 + cnst1 )
 
        !    Second order gradient
        cnst7 = cnst5 - cnst4
        cnst9 = 2.0/( cnst3+cnst2+cnst2+cnst1 )
 
        !    Use 3rd order scheme
        IF( abs(cnst7) .LT. 0.6*cnst9*abs(cf(m)-cf(k)) ) THEN
          cnst= cnst5 - ( cnst3+uflx(i) )*cnst7*cnst9/1.5
 
          !    Use doubled UNO2 scheme
        ELSE IF( dble(cnst4)*dble(cnst5) .GT. 0.d0 ) THEN
          cnst=sign(2.0, cnst5)*min( abs(cnst4), abs(cnst5) )
 
        ELSE
          !    Use minimum gradient UNO2 scheme
          cnst=sign(1.0, cnst5)*min( abs(cnst4), abs(cnst5) )
 
        ENDIF
 
        !    Mid-flux value inside central cell
        fu(i)=(cf(l) + cnst*(cnst2 - uflx(i)))*cnst8
 
        !    For negative velocity case
      ELSE
 
        !    Use central cell velocity for boundary flux.  JGLi06Apr2011
        IF( l .LE. 0) uflx(i) = uc(m)*fts
 
        !    Upstream cell length and gradient, depending on UFLX sign.
        cnst1=float( ijkcel3(n) )
        cnst4=(cf(n)-cf(m))/( cnst1 + cnst3 )
 
        !    Second order gradient
        cnst7 = cnst4 - cnst5
        cnst9 = 2.0/( cnst2+cnst3+cnst3+cnst1 )
 
        !    Use 3rd order scheme
        IF( abs(cnst7) .LT. 0.6*cnst9*abs(cf(n)-cf(l)) ) THEN
          cnst= cnst5 + ( cnst2-uflx(i) )*cnst7*cnst9/1.5
 
          !    Use doubled UNO2 scheme
        ELSE IF( dble(cnst4)*dble(cnst5) .GT. 0.d0 ) THEN
          cnst=sign(2.0, cnst5)*min( abs(cnst4), abs(cnst5) )
 
        ELSE
          !    Use minimum gradient UNO2 scheme.
          cnst=sign(1.0, cnst5)*min( abs(cnst4), abs(cnst5) )
 
        ENDIF
 
        !    Mid-flux value inside central cell M
        fu(i)=(cf(m) - cnst*(cnst3+uflx(i)))*cnst8
 
      ENDIF
 
      !    Diffusion flux by face gradient x DT
      fx(i)=cnst0*cnst5*cnst8/( clats( ij )*clats( ij ) )
 
    END DO
 
#ifdef W3_OMPG
    !$OMP END DO
    !$OMP END Parallel
#endif
 
    RETURN
  END SUBROUTINE smcxuno3
 
 
  SUBROUTINE smcyuno3(NVA, NVB, CF, VC, VFLY, AKDif, FV, FY, FTS)
 
    USE constants
    USE w3gdatmd, ONLY: ncel, mrfct, nvfc, ijkcel, ijkvfc, clatf
    USE w3gdatmd, ONLY: ijkcel4
    USE w3odatmd, ONLY: ndse, ndst
 
    IMPLICIT NONE
    INTEGER, INTENT( IN):: NVA, NVB
    REAL,    INTENT( IN):: CF(-9:NCel), VC(-9:NCel), AKDif, FTS
    REAL,    INTENT(Out):: VFLY(NVFc), FV(NVFc), FY(NVFc)
    INTEGER ::  i, j, k, L, M, N, ij
    REAL:: CNST, CNST0, CNST1, CNST2, CNST3, CNST4, CNST5, CNST6,  &
         CNST7, CNST8, CNST9
 
    !     Notice an extra side length L is multiplied to mid-flux to give correct
    !     proportion of flux into the cells.  This length will be removed by the
    !     cell length when the tracer concentration is updated.
 
    !     Diffusion Fourier no. at sub-time-step, proportional to face size,
    !     which is also equal to the sub-time-step factor FTS.
    !        CNST0=AKDif*FTS*FTS
    !     2.0 factor to cancel that in gradient CNST5.  JGLi08Mar2012
    cnst0=akdif*fts*fts*2.0
    !     Uniform diffusion coefficient for all sizes.  JGLi24Feb2012
    !        CNST0=AKDif*MRFct*FTS
 
#ifdef W3_OMPG
    !$OMP Parallel Default(Shared), Private(j, K, L, M, N),  &
    !$OMP& Private(CNST,CNST1,CNST2,CNST3,CNST4,CNST5,CNST6,CNST7,CNST8,CNST9)
    !$OMP DO
#endif
 
    DO j=nva, nvb
 
      !    Select Upstream, Central and Downstream cells
      k=ijkvfc(4,j)
      l=ijkvfc(5,j)
      m=ijkvfc(6,j)
      n=ijkvfc(7,j)
 
      !    Face bounding cell lengths and gradient
      cnst2=float( ijkcel4(l) )
      cnst3=float( ijkcel4(m) )
      cnst5=(cf(m)-cf(l))/( cnst2 + cnst3 )
 
      !    Courant number in local size-1 cell unit
      !    Multiply by multi-resolution time step factor  FTS
      cnst6=0.5*( vc(l)+vc(m) )*fts
      vfly(j) = cnst6
 
      !    Face size integer and cosine factor.
      !    CLATF is defined on V-face for SMC grid.  JGLi28Feb2012
      cnst8=clatf(j)*float( ijkvfc(3,j) )
 
      !    For positive velocity case
      IF(cnst6 >= 0.0)  THEN
 
        !    Boundary cell y-size is set equal to central cell y-size
        !    as y-boundary cell sizes are not proportional to refined
        !    inner cells but constant of the base cell y-size, and
        !    Use central cell speed for face speed.  JGLi06Apr2011
        IF( m .LE. 0 ) THEN
          vfly(j) = vc(l)*fts
          cnst3   = cnst2
        ENDIF
 
        !    Upstream cell size and irregular grid gradient, depending on VFLY.
        cnst1=float( ijkcel4(k) )
        cnst4=(cf(l)-cf(k))/( cnst2 + cnst1 )
 
        !    Second order gradient
        cnst7 = cnst5 - cnst4
        cnst9 = 2.0/( cnst3+cnst2+cnst2+cnst1 )
 
        !    Use 3rd order scheme
        IF( abs(cnst7) .LT. 0.6*cnst9*abs(cf(m)-cf(k)) ) THEN
          cnst= cnst5 - ( cnst3+vfly(j) )*cnst7*cnst9/1.5
 
          !    Use doubled UNO2 scheme
        ELSE IF( dble(cnst4)*dble(cnst5) .GT. 0.d0 ) THEN
          cnst=sign(2.0, cnst5)*min( abs(cnst4), abs(cnst5) )
 
        ELSE
 
          !    Use minimum gradient outside monotonic region
          cnst=sign(1.0, cnst5)*min( abs(cnst4), abs(cnst5) )
 
        ENDIF
 
        !    Mid-flux value multiplied by face width and cosine factor
        fv(j)=( cf(l) + cnst*(cnst2 - vfly(j)) )*cnst8
 
        !    For negative velocity case
      ELSE
 
        !    Set boundary cell y-size equal to central cell y-size and
        !    Use central cell speed for flux face speed.  JGLi06Apr2011
        IF( l .LE. 0 ) THEN
          vfly(j) = vc(m)*fts
          cnst2   = cnst3
        ENDIF
 
        !    Upstream cell size and gradient, depending on VFLY sign.
        !    Side gradients for central cell includs 0.5 factor.
        cnst1=float( ijkcel4(n) )
        cnst4=(cf(n)-cf(m))/( cnst1 + cnst3 )
 
        !    Second order gradient
        cnst7 = cnst4 - cnst5
        cnst9 = 2.0/( cnst2+cnst3+cnst3+cnst1 )
 
        !    Use 3rd order scheme
        IF( abs(cnst7) .LT. 0.6*cnst9*abs(cf(n)-cf(l)) ) THEN
          cnst= cnst5 + ( cnst2-vfly(j) )*cnst7*cnst9/1.5
 
          !    Use doubled UNO2 scheme
        ELSE IF( dble(cnst4)*dble(cnst5) .GT. 0.d0 ) THEN
          cnst=sign(2.0, cnst5)*min( abs(cnst4), abs(cnst5) )
 
        ELSE
 
          !    Use minimum gradient outside monotonic region
          cnst=sign(1.0, cnst5)*min( abs(cnst4), abs(cnst5) )
 
        ENDIF
 
        !    Mid-flux value multiplied by face width and cosine factor
        fv(j)=( cf(m) - cnst*(cnst3 + vfly(j)) )*cnst8
 
      ENDIF
 
      !    Diffusion flux by face gradient x DT x face_width x cos(lat)
      !    Multiply by multi-resolution time step factor FTS
      fy(j)=cnst0*cnst5*cnst8
 
    END DO
 
#ifdef W3_OMPG
    !$OMP END DO
    !$OMP END Parallel
#endif
 
    RETURN
  END SUBROUTINE smcyuno3
 
 
  SUBROUTINE smcxuno3r(NUA, NUB, CF, UC, UFLX, AKDif, FU, FX)
 
    USE constants
    USE w3gdatmd, ONLY: nsea, ny, ncel, nufc, ijkcel, ijkufc, clats
    USE w3gdatmd, ONLY: ijkcel3
    USE w3odatmd, ONLY: ndse, ndst
 
    IMPLICIT NONE
    INTEGER, INTENT( IN):: NUA, NUB
    REAL,    INTENT( IN):: CF(-9:NCel), UC(-9:NCel), AKDif
    REAL,    INTENT(Out):: UFLX(NUFc), FU(NUFc), FX(NUFc)
    !
    INTEGER ::  i, j, k, L, M, N, ij
    REAL:: CNST, CNST0, CNST1, CNST2, CNST3, CNST4, CNST5, CNST6, &
         CNST7, CNST8, CNST9
    !    Two layer of boundary cells are added to each boundary cell face
    !    with all boundary cell values CF(-9:0)=0.0.
 
    !    Notice an extra side length L is multiplied to mid-flux to give correct
    !    proportion of flux into the cells.  This length will be removed by the
    !    cell length when the tracer concentration is updated.
 
#ifdef W3_OMPG
    !$OMP Parallel Default(Shared), Private(i, ij, K, L, M, N),  &
    !$OMP& Private(CNST,CNST0,CNST1,CNST2,CNST3,CNST4,CNST5,CNST6,CNST7,CNST8,CNST9)
    !$OMP DO
#endif
 
    DO i=nua, nub
 
      !    Select Upstream, Central and Downstream cells
      k=ijkufc(4,i)
      l=ijkufc(5,i)
      m=ijkufc(6,i)
      n=ijkufc(7,i)
 
      !    Face bounding cell lengths and gradient
      cnst2=float( ijkcel3(l) )
      cnst3=float( ijkcel3(m) )
      cnst5=(cf(m)-cf(l))
 
      !    Averaged Courant number for base-level cell face
      cnst6= 0.5*( uc(l)+uc(m) )
      uflx(i) = cnst6
 
      !    Diffusion Fourier number in local cell size
      !    To avoid boundary cell number, use maximum of L and M.
      ij= max(l, m)
      cnst0 = 2.0/( clats(ij)*clats(ij) )
 
      !    For positive velocity case
      IF(cnst6 >= 0.0)  THEN
 
        !    Use central cell velocity for boundary flux.  JGLi06Apr2011
        IF( m .LE. 0) uflx(i) = uc(l)
 
        !    Side gradient for upstream cell as regular grid.
        cnst4=(cf(l)-cf(k))
        cnst8=(cf(m)-cf(k))
        cnst9=(cnst5-cnst4)
 
        IF( abs(cnst9) <= 0.6*abs(cnst8) )  THEN
          !    Use 3rd order scheme in limited zone, note division by 2 grid sizes
          cnst=0.5*cnst5 - (1.0+uflx(i)/cnst2)*cnst9/6.0
        ELSE IF( dble(cnst4)*dble(cnst5) .GT. 0.d0 ) THEN
          !    Use doubled minimum gradient in rest of monotonic region
          cnst=sign(1.0, cnst5)*min( abs(cnst4), abs(cnst5) )
        ELSE
          !    Use minimum gradient all region with 0.5 factor
          cnst=sign(0.5, cnst5)*min( abs(cnst4), abs(cnst5) )
        ENDIF
        !    Mid-flux value inside central cell
        fu(i)=(cf(l) + cnst*(1.0-uflx(i)/cnst2))
 
        !    For negative velocity case
      ELSE
 
        !    Use central cell velocity for boundary flux.  JGLi06Apr2011
        IF( l .LE. 0) uflx(i) = uc(m)
 
        !    Side gradient for upstream cell, depneding on UFLX sign.
        cnst4=(cf(n)-cf(m))
        cnst8=(cf(n)-cf(l))
        cnst9=(cnst4-cnst5)
 
        IF( abs(cnst9) <= 0.6*abs(cnst8) )  THEN
          !    Use 3rd order scheme in limited zone, note division by 2 grid sizes
          cnst=0.5*cnst5 + (1.0-uflx(i)/cnst3)*cnst9/6.0
        ELSE IF( dble(cnst4)*dble(cnst5) .GT. 0.d0 ) THEN
          !    Use doubled minimum gradient in rest of monotonic region
          cnst=sign(1.0, cnst5)*min( abs(cnst4), abs(cnst5) )
        ELSE
          !    Use minimum gradient outside monotonic region, include 0.5 factor
          cnst=sign(0.5, cnst5)*min( abs(cnst4), abs(cnst5) )
        ENDIF
 
        !    Mid-flux value inside central cell M
        fu(i)=(cf(m) - cnst*(1.0+uflx(i)/cnst3))
 
      ENDIF
 
      !    Diffusion flux by face gradient x DT
      fx(i)=akdif*cnst0*cnst5/(cnst2 + cnst3)
 
    END DO
 
#ifdef W3_OMPG
    !$OMP END DO
    !$OMP END Parallel
#endif
 
    ! 999  PRINT*, ' Sub SMCxUNO3r ended.'
 
    RETURN
  END SUBROUTINE smcxuno3r
 
 
  SUBROUTINE smcyuno3r(NVA, NVB, CF, VC, VFLY, AKDif, FV, FY)
 
    USE constants
    USE w3gdatmd, ONLY: nsea, ny, ncel, nvfc, ijkcel, ijkvfc, clatf
    USE w3odatmd, ONLY: ndse, ndst
 
    IMPLICIT NONE
    INTEGER, INTENT( IN):: NVA, NVB
    REAL,    INTENT( IN):: CF(-9:NCel), VC(-9:NCel), AKDif
    REAL,    INTENT(Out):: VFLY(NVFc), FV(NVFc), FY(NVFc)
    INTEGER ::  i, j, k, L, M, N, ij
    REAL    :: CNST, CNST0, CNST1, CNST2, CNST3, CNST4, CNST5, CNST6, &
         CNST7, CNST8, CNST9
 
    !     Notice an extra side length L is multiplied to mid-flux to give correct
    !     proportion of flux into the cells.  This length will be removed by the
    !     cell length when the tracer concentration is updated.
 
#ifdef W3_OMPG
    !$OMP Parallel Default(Shared), Private(j, K, L, M, N),  &
    !$OMP& Private(CNST,CNST4,CNST5,CNST6,CNST7,CNST8,CNST9)
    !$OMP DO
#endif
 
    DO j=nva, nvb
 
      !    Select Upstream, Central and Downstream cells
      k=ijkvfc(4,j)
      l=ijkvfc(5,j)
      m=ijkvfc(6,j)
      n=ijkvfc(7,j)
 
      !    Central face gradient.
      cnst5=(cf(m)-cf(l))
 
      !    Courant number in basic cell unit as dy is constant
      cnst6=0.5*( vc(l)+vc(m) )
      vfly(j) = cnst6
 
      !    Face size integer and cosine factor
      !    CLATF is defined on V-face for SMC grid.  JGLi28Feb2012
      cnst7=clatf(j)*float( ijkvfc(3,j) )
 
      !    For positive velocity case
      IF(cnst6 >= 0.0)  THEN
 
        !    Use central cell speed for flux face speed.  JGLi06Apr2011
        IF( m .LE. 0 ) vfly(j) = vc(l)
 
        !    Upstream face gradient, depending on VFLY sign.
        cnst4=(cf(l)-cf(k))
 
        !    Second gradient for 3rd scheme
        cnst8=(cf(m)-cf(k))
        cnst9=(cnst5-cnst4)
 
        IF( abs(cnst9) <= 0.6*abs(cnst8) )  THEN
          !    Use 3rd order scheme in limited zone, note division by 2 grid sizes
          cnst=0.5*cnst5-(1.0+vfly(j))*cnst9/6.0
        ELSE IF( dble(cnst4)*dble(cnst5) .GT. 0.d0 ) THEN
          !    Use doubled minimum gradient in rest of monotonic region
          cnst=sign(1.0, cnst5)*min( abs(cnst4), abs(cnst5) )
        ELSE
          !    Use minimum gradient, including 0.5 factor and central sign.
          cnst=sign(0.5, cnst5)*min( abs(cnst4), abs(cnst5) )
        ENDIF
 
        !    Mid-flux value multiplied by face width and cosine factor
        fv(j)=( cf(l) + cnst*(1.0-vfly(j)) )*cnst7
 
        !    For negative velocity case
      ELSE
 
        !    Use central cell speed for flux face speed.  JGLi06Apr2011
        IF( l .LE. 0 ) vfly(j) = vc(m)
 
        !    Side gradients for upstream face, depending on VFLY sign.
        cnst4=(cf(n)-cf(m))
 
        !    Second gradient for 3rd scheme
        cnst8=(cf(n)-cf(l))
        cnst9=(cnst4-cnst5)
 
        IF( abs(cnst9) <= 0.6*abs(cnst8) )  THEN
          !    Use 3rd order scheme in limited zone, note division by 2 grid sizes
          cnst=0.5*cnst5+(1.0-vfly(j))*cnst9/6.0
        ELSE IF( dble(cnst4)*dble(cnst5) .GT. 0.d0 ) THEN
          !    Use doubled minimum gradient in rest of monotonic region
          cnst=sign(1.0, cnst5)*min( abs(cnst4), abs(cnst5) )
        ELSE
          !    Use minimum gradient, including 0.5 factor and central sign.
          cnst=sign(0.5, cnst5)*min( abs(cnst4), abs(cnst5) )
        ENDIF
        !    Mid-flux value multiplied by face width and cosine factor
        fv(j)=( cf(m) - cnst*(1.0+vfly(j)) )*cnst7
 
      ENDIF
 
      !    Diffusion flux by face gradient x DT x face_width x cos(lat)
      fy(j)=akdif*cnst5*cnst7
 
    END DO
 
#ifdef W3_OMPG
    !$OMP END DO
    !$OMP END Parallel
#endif
 
    RETURN
  END SUBROUTINE smcyuno3r
 
  !
  ! Add optional zero-gradient bounday conditions.    JGLi08Aug2017
  !
  SUBROUTINE smcgradn(CVQ, GrdX, GrdY, L0r1)
 
    USE constants
    USE w3gdatmd, ONLY: nsea,   nufc,   nvfc,   mrfct,       &
         ijkcel, ijkufc, ijkvfc, clats, sx, sy
    USE w3gdatmd, ONLY: arctc
    USE w3odatmd, ONLY: ndse, ndst
 
    IMPLICIT NONE
    !!    New boundary conditions depending on user defined L0r1.
    !!    L0r1 = 0 will set zero at land points while L0r1 > 0 invokes
    !!    the zero-gradient boundary condition.    JGLi08Aug2017
    REAL,    INTENT( IN)::  CVQ(NSEA)
    REAL,    INTENT(Out):: GrdX(NSEA), GrdY(NSEA)
    INTEGER, INTENT( IN):: L0r1
    !
    INTEGER :: I, J, K, L, M, N
    REAL:: CNST, CNST0, CNST1, CNST2, CNST3, CNST4, CNST5, CNST6
    REAL :: DX0I, DY0I
 
    !     Use a few working arrays
    REAL,  Dimension(-9:NSEA):: CVF, AUN, AVN
 
    !     Two layer of boundary cells are added to each boundary cell face
    !     with all boundary cell default values CVF(-9:0)= 0.0.
    cvf(-9:0)  = 0.0
    cvf(1:nsea)=cvq(1:nsea)
 
    !!    Initialize arrays
    aun = 0.
    avn = 0.
    grdx = 0.
    grdy = 0.
 
    !!    Multi-resolution base-cell size defined by refined levels.
    !!    So the MRFct converts the base cell SX, SY into size-1 cell lenth.
    !!    Constant size-1 dy=DY0 and dx on Equator DX0, inverted.
    dx0i   = mrfct/ ( sx * dera * radius )
    dy0i   = mrfct/ ( sy * dera * radius )
 
#ifdef W3_OMPG
    !$OMP Parallel Default(Shared), Private(i, j, K, L, M, N),  &
    !$OMP& Private(CNST,CNST0,CNST1,CNST2,CNST3,CNST4,CNST5,CNST6)
    !$OMP DO
#endif
 
    !!    Calculate x-gradient by averaging U-face gradients.
    DO i=1, nufc
 
      !     Select Upstream, Central and Downstream cells
      l=ijkufc(5,i)
      m=ijkufc(6,i)
 
      !!      For zero-gradient boundary conditions, simply skip boundary faces.
      IF( l0r1 .EQ. 0 .OR. (l > 0 .AND. m > 0) ) THEN
 
        !         Multi-resolution SMC grid requires flux multiplied by face factor.
        cnst1=float( ijkufc(3,i) )
 
        !         Face bounding cell lengths and central gradient
        cnst2=float( ijkcel(3,l) )
        cnst3=float( ijkcel(3,m) )
 
        !         Side gradients over 2 cell lengths for central cell.
        !         Face size factor is also included for average.
        cnst5=cnst1*(cvf(m)-cvf(l))/(cnst2+cnst3)
#if defined W3_B4B && defined W3_OMPG
        cnst5=int(cnst5 * 1.0e6)   ! CB: B4B
#endif
 
        !! Replace CRITICAL with ATOMIC.  JGLi15Jan2019
        !! !$OMP CRITICAL
        !    Store side gradient in two neighbouring cells
        !! Remove boundary cell flux update or L M > 0.  JGLi28Mar2019
        IF( l > 0 ) THEN
#ifdef W3_OMPG
          !$OMP ATOMIC
#endif
          aun(l) = aun(l) + cnst5
        ENDIF
        IF( m > 0 ) THEN
#ifdef W3_OMPG
          !$OMP ATOMIC
#endif
          aun(m) = aun(m) + cnst5
        ENDIF
        !! !$OMP END CRITICAL
 
      ENDIF
    END DO
 
#ifdef W3_OMPG
    !$OMP END DO
#endif
 
#if defined W3_B4B && defined W3_OMPG
    !$OMP SINGLE
    aun = aun / 1.0e6  ! CB B4B
    !$OMP END SINGLE
#endif
 
    !  Assign averaged side-gradient to GrdX, plus latitude factor
    !  Note averaging over 2 times of cell y-width factor but AUN
    !  has already been divied by two cell lengths.
 
#ifdef W3_OMPG
    !$OMP DO
#endif
 
    DO n=1, nsea
      !       Cell y-size IJKCel(4,i) is used to cancel the face size-factor in AUN.
      !       Plus the actual physical length scale for size-1 cell.
      !       Note polar cell (if any) AUN = 0.0 as it has no U-face.
      grdx(n)=dx0i*aun(n)/( clats(n)*ijkcel(4,n) )
 
    ENDDO
 
#ifdef W3_OMPG
    !$OMP END DO
    !$OMP DO
#endif
 
    !!    Calculate y-gradient by averaging V-face gradients.
    DO j=1, nvfc
 
      !       Select Central and Downstream cells
      l=ijkvfc(5,j)
      m=ijkvfc(6,j)
 
      !!      For zero-gradient boundary conditions, simply skip boundary faces.
      IF( l0r1 .EQ. 0 .OR. (l > 0 .AND. m > 0) ) THEN
 
        !       Face size is required for multi-resolution grid.
        cnst1=real( ijkvfc(3,j) )
 
        !         Cell y-length of UCD cells
        cnst2=real( ijkcel(4,l) )
        cnst3=real( ijkcel(4,m) )
 
        !         Side gradients over 2 cell lengths for central cell.
        !         Face size factor is also included for average.
        cnst6=cnst1*(cvf(m)-cvf(l))/(cnst2+cnst3)
#if defined W3_B4B && defined W3_OMPG
        cnst6 = int(cnst6 * 1.0e6) ! CB B4B
#endif
 
        !! Replace CRITICAL with ATOMIC.  JGLi15Jan2019
        !! !$OMP CRITICAL
        !! Remove boundary cell flux update or L M > 0.  JGLi28Mar2019
        IF( l > 0 ) THEN
          !    Store side gradient in two neighbouring cells
#ifdef W3_OMPG
          !$OMP ATOMIC
#endif
          avn(l) = avn(l) + cnst6
        ENDIF
        IF( m > 0 ) THEN
#ifdef W3_OMPG
          !$OMP ATOMIC
#endif
          avn(m) = avn(m) + cnst6
        ENDIF
        !! !$OMP END CRITICAL
 
      ENDIF
    END DO
 
#ifdef W3_OMPG
    !$OMP END DO
#endif
 
#if defined W3_B4B && defined W3_OMPG
    !$OMP SINGLE
    avn = avn / 1.0e6  !CB B4B
    !$OMP END SINGLE
#endif
 
#ifdef W3_OMPG
    !$OMP DO
#endif
 
    !  Assign averaged side-gradient to GrdY.
    DO n=1, nsea
      !  AV is divided by the cell x-size IJKCel(3,i) to cancel face
      !  size-factor, and physical y-distance of size-1 cell.
      grdy(n)=dy0i*avn(n)/real( ijkcel(3,n) )
 
    END DO
 
#ifdef W3_OMPG
    !$OMP END DO
    !$OMP END Parallel
#endif
 
    !!Li  Y-gradient for polar cell in Arctic part is set to zero.
    IF( arctc ) grdy(nsea) = 0.0
 
    RETURN
  END SUBROUTINE smcgradn
 
 
  SUBROUTINE smcaverg(CVQ)
 
    USE constants
    USE w3gdatmd, ONLY: nsea,   nufc,   nvfc,    &
         ijkcel, ijkufc, ijkvfc,  &
         ijkufc5, ijkufc6
    USE w3gdatmd, ONLY: arctc
    USE w3odatmd, ONLY: ndse, ndst
 
    IMPLICIT NONE
    REAL, INTENT(INOUT) :: CVQ(-9:NSEA)
    !
    INTEGER :: I, J, K, L, M, N
    REAL :: CNST, CNST0, CNST1, CNST2, CNST3, CNST4, CNST5, CNST6
 
    !     Use a few working arrays
    REAL, Dimension(-9:NSEA) :: CVF, AUN, AVN
 
    !     Two layer of boundary cells are added to each boundary cell face
    !     with all boundary cell values stored in CVF(-9:0).
    cvf=cvq
 
    !!    Initialize arrays
    aun = 0.
    avn = 0.
 
    !!Li  Save polar cell value if any.
    cnst0 = cvq(nsea)
 
#ifdef W3_OMPG
    !$OMP Parallel Default(Shared), Private(i, j, L, M, n),  &
    !$OMP& Private(CNST3,CNST4,CNST5,CNST6)
    !$OMP DO
#endif
 
    !!    Calculate x-gradient by averaging U-face gradients.
    DO i=1, nufc
 
      !    Select Upstream, Central and Downstream cells
      l=ijkufc5(i)
      m=ijkufc6(i)
 
      !       Multi-resolution SMC grid requires flux multiplied by face factor.
      cnst5=real( ijkufc(3,i) )*(cvf(m)+cvf(l))
#ifdef W3_B4B
      !OMPG           CNST5=int(CNST5 * 1.0e6)
#endif
 
      !! Replace CRITICAL with ATOMIC.  JGLi15Jan2019
      !! !$OMP CRITICAL
      !    Store side gradient in two neighbouring cells
      !! Remove boundary cell flux update or L M > 0.  JGLi28Mar2019
      IF( l > 0 ) THEN
#ifdef W3_OMPG
        !$OMP ATOMIC
#endif
        aun(l) = aun(l) + cnst5
      ENDIF
      IF( m > 0 ) THEN
#ifdef W3_OMPG
        !$OMP ATOMIC
#endif
        aun(m) = aun(m) + cnst5
      ENDIF
      !! !$OMP END CRITICAL
 
    END DO
 
#ifdef W3_OMPG
    !$OMP END DO
#endif
 
#if defined W3_B4B && defined W3_OMPG
    !$OMP SINGLE
    aun = aun / 1.0e6  !CB B4B
    !$OMP END SINGLE
#endif
 
#ifdef W3_OMPG
    !$OMP DO
#endif
 
    !!    Calculate y-gradient by averaging V-face gradients.
    DO j=1, nvfc
 
      !       Select Central and Downstream cells
      l=ijkvfc(5,j)
      m=ijkvfc(6,j)
 
      !       Face size is required for multi-resolution grid.
      cnst6=real( ijkvfc(3,j) )*(cvf(m)+cvf(l))
#if defined W3_B4B && defined W3_OMPG
      cnst6=int(cnst6 * 1e6)
#endif
 
      !! Replace CRITICAL with ATOMIC.  JGLi15Jan2019
      !! !$OMP CRITICAL
      !    Store side gradient in two neighbouring cells
      !! Remove boundary cell flux update or L M > 0.  JGLi28Mar2019
      IF( l > 0 ) THEN
#ifdef W3_OMPG
        !$OMP ATOMIC
#endif
        avn(l) = avn(l) + cnst6
      ENDIF
      IF( m > 0 ) THEN
#ifdef W3_OMPG
        !$OMP ATOMIC
#endif
        avn(m) = avn(m) + cnst6
      ENDIF
      !! !$OMP END CRITICAL
 
    END DO
 
#ifdef W3_OMPG
    !$OMP END DO
#endif
 
#if defined W3_B4B && defined W3_OMPG
    !$OMP SINGLE
    avn = avn / 1.0e6  !CB B4B
    !$OMP END SINGLE
#endif
 
#ifdef W3_OMPG
    !$OMP DO
#endif
 
    !  Assign averaged value back to CVQ.
    DO n=1, nsea
 
      cnst3=0.125/real( ijkcel(3,n) )
      cnst4=0.125/real( ijkcel(4,n) )
      !       AUN is divided by the cell y-size IJKCel(4,n) and AVN by
      !       the cell x-size IJKCel(3,n) to cancel face size factors.
      cvq(n)= aun(n)*cnst4 + avn(n)*cnst3
 
    END DO
 
#ifdef W3_OMPG
    !$OMP END DO
    !$OMP END Parallel
#endif
 
    !!Li  Polar cell (if any) keep original value.
    IF( arctc ) cvq(nsea) = cnst0
 
    ! 999  PRINT*, ' Sub SMCAverg ended.'
 
    RETURN
  END SUBROUTINE smcaverg
 
 
  SUBROUTINE smcgtcrfr(CoRfr, SpeTHK)
    USE constants
    USE w3gdatmd, ONLY: nk, nth, dth, ctmax
 
    IMPLICIT NONE
    REAL, INTENT(IN)   ::  CoRfr(NTH, NK)
    REAL, INTENT(INOUT):: SpeTHK(NTH, NK)
    INTEGER ::  I, J, K, L, M, N
    REAL, Dimension(NTH):: SpeGCT, Spectr
    REAL:: CNST, CNST0, CNST1, CNST2, CNST3, CNST4, CNST5, CNST6
 
    !     Loop through NK spectral bins.
    DO n=1, nk
 
      !!      Asign cell spectrum to temporary variable Spcetr
      spectr=spethk(1:nth,n)
      spegct=0.0
 
      !!      Loop through NTH directional bins for each cell spectrum
      DO j=1, nth
 
        !         GCT + refraction Courant number for this dirctional bin
        cnst6=corfr(j,n)
 
        !         Work out integer number of bins to be skipped.
        !         If K is great than NTH, full circle turning is removed.
        cnst5=abs( cnst6 )
        k= mod( int(cnst5), nth )
 
        !         Partitioning faraction of the spectral component
        cnst1=cnst5 - float( int(cnst5) )
        cnst2=1.0 - cnst1
 
        !         For positive turning case
        IF(cnst6 > 0.0)  THEN
 
          !           Select the upstream and downstream bins to rotate in, wrap at end
          l=j+k
          m=j+k+1
          IF( l .GT. nth ) l = l - nth
          IF( m .GT. nth ) m = m - nth
 
          !!          Divide the j bin energy by fraction of CNST6 and store in SpeGCT
          spegct(l)=spegct(l)+spectr(j)*cnst2
          spegct(m)=spegct(m)+spectr(j)*cnst1
 
          !         For negative or no turning case
        ELSE
 
          !           Select the upstream and downstream bins to rotate in, wrap at end
          l=j-k
          m=j-k-1
          IF( l .LT. 1 ) l = l + nth
          IF( m .LT. 1 ) m = m + nth
 
          !!          Divide the bin energy by fraction of CNST6 and store in SpeGCT
          spegct(l)=spegct(l)+spectr(j)*cnst2
          spegct(m)=spegct(m)+spectr(j)*cnst1
 
        ENDIF
 
        !!      End of directional loop j
      END DO
 
      !!      Store GCT spectrum
      spethk(1:nth,n) = spegct
 
      !!    End of cell loop n
    END DO
 
    ! 999  PRINT*, ' Sub SMCGtCrfr ended.'
 
    RETURN
  END SUBROUTINE smcgtcrfr
 
 
  !
  ! Fix bug on CFL limiter and add positive filter.  JGLi28Jun2017
  !
  SUBROUTINE smckuno2(CoRfr, SpeTHK, DKC, DKS)
 
    USE constants
    USE w3gdatmd, ONLY: nk, nk2, nth, dth, xfr, ctmax
 
    IMPLICIT NONE
    REAL, INTENT(IN)   ::  CoRfr(NTH, 0:NK), DKC(0:NK+1), DKS(-1:NK+1)
    REAL, INTENT(INOUT):: SpeTHK(NTH, NK)
    INTEGER ::  I, J, K, L, M, N
    REAL, Dimension(-1:NK+2):: SpeRfr, Spectr, SpeFlx
    REAL:: CNST, CNST0, CNST1, CNST2, CNST3, CNST4, CNST5, CNST6
 
    cnst=xfr**(-7)
 
    DO n=1, nth
 
      !!      Asign cell spectrum to temporary variable Spcetr
      spectr(-1)  =0.0
      spectr( 0)  =0.0
      spectr(1:nk)=spethk(n,1:nk)
      spectr(nk+1)=spectr(nk  )*cnst
      spectr(nk+2)=spectr(nk+1)*cnst
 
      !!      Calculate k-space gradient for NK+2 faces by UNO2 scheme
      sperfr(-1)= 0.0
      DO j=-1, nk+1
        sperfr(j)=(spectr(j+1)-spectr(j))/dks(j)
      ENDDO
 
      !!      Calculate k-space fluxes for NK+1 faces by UNO2 scheme
      DO j=0, nk
 
        !         Final refraction Courant number for this k-space face
        cnst6=corfr(n,j)
        !!        Note CoRfr is CFL for k but without dividing dk.
 
        !         For positive case
        IF(cnst6 > 0.0)  THEN
 
          cnst0 = min( ctmax*dkc(j), cnst6 )
          speflx(j) = cnst0*( spectr(j) + sign(0.5, sperfr(j))*(dkc(j)-cnst0)  &
               *min( abs(sperfr(j-1)), abs(sperfr(j)) ) )
 
          !         For negative or no turning case
        ELSE
 
          cnst0 = min( ctmax*dkc(j+1), -cnst6 )
          speflx(j) = -cnst0*( spectr(j+1) - sign(0.5, sperfr(j))*(dkc(j+1)-cnst0) &
               *min( abs(sperfr(j+1)), abs(sperfr(j)) ) )
 
        ENDIF
 
        !!      End of flux loop j
      END DO
 
      !!      Update spectrum for the given direction
      DO j=1, nk
        !         Final refraction Courant number for this k-space face
        !         SpeTHK(n, j) = Spectr(j) + (SpeFlx(j-1) - SpeFlx(j))/DKC(j)
        !         Add positive filter in case negative values.  JGLi27Jun2017
        spethk(n, j) = max( 0.0, spectr(j)+(speflx(j-1)-speflx(j))/dkc(j) )
      END DO
 
      !!    End of directional loop n
    END DO
 
    ! 999   PRINT*, ' Sub SMCkUNO2 ended.'
 
    RETURN
  END SUBROUTINE smckuno2
 
 
  SUBROUTINE smcdhxy
    USE constants
    USE w3gdatmd, ONLY: nx, ny, nsea, mapsta, mapfs, mrfct, ijkcel,  &
         clats, nth, dth, esin, ecos, refran, dmin
    USE w3gdatmd, ONLY: nglo, angarc, arctc
    USE w3adatmd, ONLY: dw, dddx, dddy, dhdx, dhdy, dhlmt
    USE w3odatmd, ONLY: ndse, ndst
 
    IMPLICIT NONE
 
    INTEGER :: I, J, K, L, M, N
    REAL :: CNST, CNST0, CNST1, CNST2, CNST3, CNST4, CNST5, CNST6
    REAL, Dimension(NSEA) :: HCel, GrHx, GrHy
    !     REAL, Dimension(-9:NSEA) :: HCel
 
    !!    Assign water depth to HCel from DW integer values.
    !!    set half the minimum depth DMIN for negative cells.
    !     HCel(-9:0) = 0.5*DMIN
    hcel(1:nsea)= dw(1:nsea)
 
    !!    Reset shallow water depth with minimum depth
#ifdef W3_OMPG
    !$OMP Parallel DO Private(k)
#endif
    DO k=1, nsea
      IF(dw(k) .LT. dmin)  hcel(k)=dmin
    ENDDO
#ifdef W3_OMPG
    !$OMP END Parallel DO
#endif
 
    !!    Initialize full grid gradient arrays
    dddx = 0.
    dddy = 0.
 
    !!    Use zero-gradient boundary condition or L0r1 > 0
    l = 1
 
    !!    Calculate sea point water depth gradient
    CALL smcgradn(hcel, grhx, grhy, l)
 
    !!    Pass gradient values to DHDX, DHDY
    dhdx(1:nsea) = grhx
    dhdy(1:nsea) = grhy
 
    !!   Apply limiter to depth-gradient and copy to full grid.
#ifdef W3_OMPG
    !$OMP Parallel DO Private(i,j,k,m,n, CNST0, CNST1, CNST2)
#endif
    DO n=1,nsea
 
      !  A limiter of gradient <= 0.1 is applied.
      IF( abs( dhdx(n) ) .GT.  0.1) dhdx(n)=sign( 0.1, dhdx(n) )
      IF( abs( dhdy(n) ) .GT.  0.1) dhdy(n)=sign( 0.1, dhdy(n) )
 
      !! Asign DHDX value to full grid variable DDDX
      i= ijkcel(1,n)/mrfct + 1
      j= ijkcel(2,n)/mrfct + 1
      k= max(1, ijkcel(3,n)/mrfct)
      m= max(1, ijkcel(4,n)/mrfct)
      dddx(j:j+m-1,i:i+k-1)  = dhdx(n)
      dddy(j:j+m-1,i:i+k-1)  = dhdy(n)
 
      !Li  Depth gradient in the Arctic part has to be rotated into
      !Li  the map-east system for calculation of refraction.
      IF( arctc .AND. n .GT. nglo ) THEN
        cnst0 = angarc(n - nglo)*dera
        cnst1 = dhdx(n)*cos(cnst0) - dhdy(n)*sin(cnst0)
        cnst2 = dhdx(n)*sin(cnst0) + dhdy(n)*cos(cnst0)
        dhdx(n) = cnst1
        dhdy(n) = cnst2
      ENDIF
 
    END DO
#ifdef W3_OMPG
    !$OMP END Parallel DO
#endif
 
    !! Calculate the depth gradient limiter for refraction.
#ifdef W3_T
    l = 0 !CB - added T switch
#endif
 
#ifdef W3_OMPG
    !$OMP Parallel DO Private(i, n, CNST4, CNST6)
#endif
    DO n=1,nsea
 
      !Li   Work out magnitude of depth gradient
      cnst4 = 1.0001*sqrt(dhdx(n)*dhdx(n) + dhdy(n)*dhdy(n))
      !
      !Li   Directional depedent depth gradient limiter.  JGLi16Jun2011
      IF ( cnst4 .GT. 1.0e-5 ) THEN
 
#if defined W3_T && defined W3_OMPG
        !$OMP ATOMIC Update   !CB - added T switch
#endif
        l = l + 1       !CB - added T switch
#if defined W3_T && defined W3_OMPG
        !$OMP END ATOMIC      !CB - added T switch
#endif
 
        DO i=1, nth
          !Li       Refraction is done only when depth gradient is non-zero.
          !Li       Note ACOS returns value between [0, Pi), always positive.
          cnst6 = acos(-(dhdx(n)*ecos(i)+dhdy(n)*esin(i))/cnst4 )
          !Li   User-defined refraction limiter added.   JGLi09Jan2012
          dhlmt(i,n)=min(refran, 0.75*min(cnst6,abs(pi-cnst6)))/dth
        END DO
        !Li   Output some values for inspection.  JGLi22Jul2011
#ifdef W3_T
        IF( mod(n, 1000) .EQ. 0 )   &
             WRITE(ndst,'(i8,18F5.1)' ) n, (dhlmt(i,n), i=1,18)
#endif
 
      ELSE
        dhlmt(:,n) = 0.0
      ENDIF
 
    ENDDO
#ifdef W3_OMPG
    !$OMP END Parallel DO
#endif
 
#ifdef W3_T
    WRITE(ndst,*) ' No. Refraction points =', l
#endif
 
#ifdef W3_T
999 print*, ' Sub SMCDHXY ended.'
#endif
 
    RETURN
  END SUBROUTINE smcdhxy
 
 
  SUBROUTINE smcdcxy
    USE constants
    USE w3gdatmd, ONLY: nx, ny, nsea, mapsta, mapfs, mrfct, ijkcel
    USE w3gdatmd, ONLY: nglo, angarc, arctc
    USE w3adatmd, ONLY: cx, cy, dcxdx, dcxdy, dcydx, dcydy
    USE w3odatmd, ONLY: ndse, ndst
 
    IMPLICIT NONE
 
    INTEGER :: I, J, K, L, M, N
    REAL :: CNST, CNST0, CNST1, CNST2, CNST3, CNST4, CNST5, CNST6
    REAL, Dimension(NSEA) :: CXCY, GrHx, GrHy
    !     REAL, Dimension(-9:NSEA) :: CXCY
 
    !!    Assign current CX speed to CXCY and set negative cells.
    !     CXCY(-9:0) = 0.0
    !!    Use zero-gradient boundary condition or L0r1 > 0
    l = 1
    cxcy(1:nsea)= cx(1:nsea)
 
    !!   Initialize full grid gradient arrays
    dcxdx = 0.0
    dcxdy = 0.0
 
    !!    Calculate sea point water depth gradient
    CALL smcgradn(cxcy, grhx, grhy, l)
 
    !!   Apply limiter to CX-gradient and copy to full grid.
#ifdef W3_OMPG
    !$OMP Parallel DO Private(i, j, k, m, n, CNST0, CNST1, CNST2)
#endif
    DO n=1,nsea
 
      !       A limiter of gradient <= 0.01 is applied.
      IF( abs( grhx(n) ) .GT.  0.01) grhx(n)=sign( 0.01, grhx(n) )
      IF( abs( grhy(n) ) .GT.  0.01) grhy(n)=sign( 0.01, grhy(n) )
 
      !Li     Current gradient in the Arctic part has to be rotated into
      !Li     the map-east system for calculation of refraction.
      IF( arctc .AND. n .GT. nglo ) THEN
        cnst0 = angarc(n - nglo)*dera
        cnst1 = grhx(n)*cos(cnst0) - grhy(n)*sin(cnst0)
        cnst2 = grhx(n)*sin(cnst0) + grhy(n)*cos(cnst0)
        grhx(n) = cnst1
        grhy(n) = cnst2
      ENDIF
 
      !! Asign CX gradients to full grid variable DCXDX/Y
      i= ijkcel(1,n)/mrfct + 1
      j= ijkcel(2,n)/mrfct + 1
      k= max(1, ijkcel(3,n)/mrfct)
      m= max(1, ijkcel(4,n)/mrfct)
      dcxdx(j:j+m-1,i:i+k-1)  = grhx(n)
      dcxdy(j:j+m-1,i:i+k-1)  = grhy(n)
 
    END DO
#ifdef W3_OMPG
    !$OMP END Parallel DO
#endif
 
    !!    Assign current CY speed to CXCY and set negative cells.
    !     CXCY(-9:0) = 0.0
    !!    Use zero-gradient boundary condition or L0r1 > 0
    l = 1
    cxcy(1:nsea)= cy(1:nsea)
 
    !!    Initialize full grid gradient arrays
    dcydx = 0.0
    dcydy = 0.0
 
    !!    Calculate sea point water depth gradient
    CALL smcgradn(cxcy, grhx, grhy, l)
 
    !!    Apply limiter to CX-gradient and copy to full grid.
#ifdef W3_OMPG
    !$OMP Parallel DO Private(i, j, k, m, n, CNST0, CNST1, CNST2)
#endif
    DO n=1,nsea
 
      !!      A limiter of gradient <= 0.1 is applied.
      IF( abs( grhx(n) ) .GT.  0.01) grhx(n)=sign( 0.01, grhx(n) )
      IF( abs( grhy(n) ) .GT.  0.01) grhy(n)=sign( 0.01, grhy(n) )
 
      !!      Current gradient in the Arctic part has to be rotated into
      !!      the map-east system for calculation of refraction.
      IF( arctc .AND. n .GT. nglo ) THEN
        cnst0 = angarc(n - nglo)*dera
        cnst1 = grhx(n)*cos(cnst0) - grhy(n)*sin(cnst0)
        cnst2 = grhx(n)*sin(cnst0) + grhy(n)*cos(cnst0)
        grhx(n) = cnst1
        grhy(n) = cnst2
      ENDIF
 
      !!      Asign CX gradients to full grid variable DCXDX/Y
      i= ijkcel(1,n)/mrfct + 1
      j= ijkcel(2,n)/mrfct + 1
      k= max(1, ijkcel(3,n)/mrfct)
      m= max(1, ijkcel(4,n)/mrfct)
      dcydx(j:j+m-1,i:i+k-1)  = grhx(n)
      dcydy(j:j+m-1,i:i+k-1)  = grhy(n)
 
    END DO
#ifdef W3_OMPG
    !$OMP END Parallel DO
#endif
 
#ifdef W3_T
999 print*, ' Sub SMCDCXY ended.'
#endif
 
    RETURN
  END SUBROUTINE smcdcxy
 
  !/
  !/ ------------------------------------------------------------------- /
  !/ ------------------------------------------------------------------- /
  !/
  SUBROUTINE w3gathsmc ( ISPEC, FIELD )
    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH-III           NOAA/NCEP |
    !/                  |           H. L. Tolman            |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         13-Jun-2006 |
    !/                  +-----------------------------------+
    !/
    !/    04-Jan-1999 : Distributed FORTRAN 77 version.     ( version 1.18 )
    !/    13-Jan-2000 : Upgrade to FORTRAN 90               ( version 2.00 )
    !/                  Major changes to logistics.
    !/    29-Dec-2004 : Multiple grid version.              ( version 3.06 )
    !/    13-Jun-2006 : Split STORE in G/SSTORE             ( version 3.09 )
    !/     9-Dec-2010 : Adapted for SMC grid propagtion. JGLi
    !/
    !  1. Purpose :
    !
    !     Gather spectral bin information into a propagation field array.
    !
    !  2. Method :
    !
    !     Direct copy or communication calls (MPP version).
    !
    !  3. Parameters :
    !
    !     Parameter list
    !     ----------------------------------------------------------------
    !       ISPEC   Int.   I   Spectral bin considered.
    !       FIELD   R.A.   O   Full field to be propagated.
    !     ----------------------------------------------------------------
    !
    !  4. Subroutines used :
    !
    !      Name      Type  Module   Description
    !     ----------------------------------------------------------------
    !      STRACE    Subr. W3SERVMD Subroutine tracing.
    !
    !      MPI_STARTALL, MPI_WAITALL
    !                Subr. mpif.h   MPI persistent comm. routines (!/MPI).
    !     ----------------------------------------------------------------
    !
    !  5. Called by :
    !
    !      Name      Type  Module   Description
    !     ----------------------------------------------------------------
    !      W3WAVE    Subr. W3WAVEMD Actual wave model routine.
    !     ----------------------------------------------------------------
    !
    !  6. Error messages :
    !
    !       None.
    !
    !  7. Remarks :
    !
    !     - The field is extracted but not converted.
    !     - Array FIELD is not initialized.
    !     - MPI version requires posing of send and receive calls in
    !       W3WAVE to match local calls.
    !     - MPI version does not require an MPI_TESTALL call for the
    !       posted gather operation as MPI_WAITALL is mandatory to
    !       reset persistent communication for next time step.
    !     - MPI version allows only two new pre-fetch postings per
    !       call to minimize chances to be slowed down by gathers that
    !       are not yet needed, while maximizing the pre-loading
    !       during the early (low-frequency) calls to the routine
    !       where the amount of calculation needed for proagation is
    !       the largest.
    !
    !  8. Structure :
    !
    !     See source code.
    !
    !  9. Switches :
    !
    !     !/SHRD  Switch for message passing method.
    !     !/MPI   Id.
    !
    !     !/S     Enable subroutine tracing.
    !     !/MPIT  MPI test output.
    !
    ! 10. Source code :
    !
    !/ ------------------------------------------------------------------- /
#ifdef W3_S
    USE w3servmd, ONLY: strace
#endif
    !/
    USE w3gdatmd, ONLY: nspec, nx, ny, nsea, nseal, ncel, mapsf
    USE w3wdatmd, ONLY: a => va
#ifdef W3_MPI
    USE w3adatmd, ONLY: mpibuf, bstat, ibfloc, isploc, bispl, &
         nsploc, nrqsg2, irqsg2, gstore
    USE w3odatmd, ONLY: ndst, iaproc, naproc
#endif
    !/
    IMPLICIT NONE
    !
#ifdef W3_MPI
    include "mpif.h"
#endif
    !/
    !/ ------------------------------------------------------------------- /
    !/ Parameter list
    !/
    INTEGER, INTENT(IN)     :: ISPEC
    REAL,    INTENT(OUT)    :: FIELD(NCel)
    !/
    !/ ------------------------------------------------------------------- /
    !/ Local parameters
    !/
#ifdef W3_SHRD
    INTEGER                 :: ISEA, IXY
#endif
#ifdef W3_MPI
    INTEGER                 :: STATUS(MPI_STATUS_SIZE,NSPEC),  &
         IOFF, IERR_MPI, JSEA, ISEA,     &
         IXY, IS0, IB0, NPST, J
#endif
#ifdef W3_S
    INTEGER, SAVE           :: IENT
#endif
    !/
    !/ ------------------------------------------------------------------- /
    !/
#ifdef W3_S
    CALL strace (ient, 'W3GATH')
#endif
    !
    !      FIELD  = 0.
    !
    ! 1.  Shared memory version ------------------------------------------ /
    !
#ifdef W3_SHRD
    DO isea=1, nsea
      field(isea) = a(ispec,isea)
    END DO
    !
    RETURN
#endif
    !
    ! 2.  Distributed memory version ( MPI ) ----------------------------- /
    ! 2.a Update counters
    !
#ifdef W3_MPI
    isploc = isploc + 1
    ibfloc = ibfloc + 1
    IF ( ibfloc .GT. mpibuf ) ibfloc = 1
    !
    ! 2.b Check status of present buffer
    ! 2.b.1 Scatter (send) still in progress, wait to end
    !
    IF ( bstat(ibfloc) .EQ. 2 ) THEN
      ioff =  1 + (bispl(ibfloc)-1)*nrqsg2
      IF ( nrqsg2 .GT. 0 ) CALL                              &
           mpi_waitall ( nrqsg2, irqsg2(ioff,2),             &
           status, ierr_mpi )
      bstat(ibfloc) = 0
    END IF
    !
    ! 2.b.2 Gather (recv) not yet posted, post now
    !
    IF ( bstat(ibfloc) .EQ. 0 ) THEN
      bstat(ibfloc) = 1
      bispl(ibfloc) = isploc
      ioff =  1 + (isploc-1)*nrqsg2
      IF ( nrqsg2 .GT. 0 ) CALL                              &
           mpi_startall ( nrqsg2, irqsg2(ioff,1), ierr_mpi )
    END IF
    !
    ! 2.c Put local spectral densities in store
    !
    DO jsea=1, nseal
      gstore(iaproc+(jsea-1)*naproc,ibfloc) = a(ispec,jsea)
    END DO
    !
    ! 2.d Wait for remote spectral densities
    !
    ioff =  1 + (bispl(ibfloc)-1)*nrqsg2
    IF ( nrqsg2 .GT. 0 ) CALL                                  &
         mpi_waitall ( nrqsg2, irqsg2(ioff,1), status, ierr_mpi )
    !
    ! 2.e Convert storage array to field.
    !
    DO isea=1, nsea
      field(isea) = gstore(isea,ibfloc)
    END DO
    !
    ! 2.f Pre-fetch data in available buffers
    !
    is0    = isploc
    ib0    = ibfloc
    npst   = 0
    DO j=1, mpibuf-1
      is0    = is0 + 1
      IF ( is0 .GT. nsploc ) EXIT
      ib0    = 1 + mod(ib0,mpibuf)
      IF ( bstat(ib0) .EQ. 0 ) THEN
        bstat(ib0) = 1
        bispl(ib0) = is0
        ioff       = 1 + (is0-1)*nrqsg2
        IF ( nrqsg2 .GT. 0 ) CALL                            &
             mpi_startall ( nrqsg2, irqsg2(ioff,1), ierr_mpi )
        npst       = npst + 1
      END IF
      IF ( npst .GE. 2 ) EXIT
    END DO
    RETURN
#endif
    !
    !/ End of W3GATHSMC ----------------------------------------------------- /
    !/
  END SUBROUTINE w3gathsmc
  !
  !/ ------------------------------------------------------------------- /
  SUBROUTINE w3scatsmc ( ISPEC, MAPSTA, FIELD )
    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH-III           NOAA/NCEP |
    !/                  |           H. L. Tolman            |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         13-Jun-2006 |
    !/                  +-----------------------------------+
    !/
    !/    04-Jan-1999 : Distributed FORTRAN 77 version.     ( version 1.18 )
    !/    13-Jan-2000 : Upgrade to FORTRAN 90               ( version 2.00 )
    !/                  Major changes to logistics.
    !/    28-Dec-2004 : Multiple grid version.              ( version 3.06 )
    !/    07-Sep-2005 : Updated boundary conditions.        ( version 3.08 )
    !/    13-Jun-2006 : Split STORE in G/SSTORE             ( version 3.09 )
    !/     9-Dec-2010 : Adapted for SMC grid propagtion.     JGLi09Dec2010
    !/    16-Jan-2012 : Remove MAPSTA checking for SMC grid. JGLi16Jan2012
    !/
    !/
    !  1. Purpose :
    !
    !     'Scatter' data back to spectral storage after propagation.
    !
    !  2. Method :
    !
    !     Direct copy or communication calls (MPP version).
    !     See also W3GATH.
    !
    !  3. Parameters :
    !
    !     Parameter list
    !     ----------------------------------------------------------------
    !       ISPEC   Int.   I   Spectral bin considered.
    !       MAPSTA  I.A.   I   Status map for spatial grid.
    !       FIELD   R.A.   I   SMC grid field to be propagated.
    !     ----------------------------------------------------------------
    !
    !  4. Subroutines used :
    !
    !      Name      Type  Module   Description
    !     ----------------------------------------------------------------
    !      STRACE    Subr. W3SERVMD Subroutine tracing.
    !
    !      MPI_STARTALL, MPI_WAITALL, MPI_TESTALL
    !                Subr. mpif.h   MPI persistent comm. routines (!/MPI).
    !     ----------------------------------------------------------------
    !
    !  5. Called by :
    !
    !      Name      Type  Module   Description
    !     ----------------------------------------------------------------
    !      W3WAVE    Subr. W3WAVEMD Actual wave model routine.
    !     ----------------------------------------------------------------
    !
    !  6. Error messages :
    !
    !     None.
    !
    !  7. Remarks :
    !
    !     - The field is put back but not converted !
    !     - MPI persistent communication calls initialize in W3MPII.
    !     - See W3GATH and W3MPII for additional comments on data
    !       buffering.
    !
    !  8. Structure :
    !
    !     See source code.
    !
    !  9. Switches :
    !
    !     !/SHRD  Switch for message passing method.
    !     !/MPI   Id.
    !
    !     !/S     Enable subroutine tracing.
    !
    ! 10. Source code :
    !
    !/ ------------------------------------------------------------------- /
#ifdef W3_S
    USE w3servmd, ONLY: strace
#endif
    !/
    USE w3gdatmd, ONLY: nspec, nx, ny, nsea, ncel, nseal, mapsf
    USE w3wdatmd, ONLY: a => va
#ifdef W3_MPI
    USE w3adatmd, ONLY: mpibuf, bstat, ibfloc, isploc, bispl, &
         nsploc, nrqsg2, irqsg2, sstore
    USE w3odatmd, ONLY: iaproc, naproc
#endif
    USE w3odatmd, ONLY: ndst
    !/
    IMPLICIT NONE
    !
#ifdef W3_MPI
    include "mpif.h"
#endif
    !/
    !/ ------------------------------------------------------------------- /
    !/ Parameter list
    !/
    INTEGER, INTENT(IN)     :: ISPEC, MAPSTA(NY*NX)
    REAL,    INTENT(IN)     :: FIELD(NCel)
    !/
    !/ ------------------------------------------------------------------- /
    !/ Local parameters
    !/
#ifdef W3_SHRD
    INTEGER                 :: ISEA, IXY
#endif
#ifdef W3_MPI
    INTEGER                 :: ISEA, IXY, IOFF, IERR_MPI, J,   &
         STATUS(MPI_STATUS_SIZE,NSPEC),  &
         JSEA, IB0
#endif
#ifdef W3_S
    INTEGER, SAVE           :: IENT
#endif
#ifdef W3_MPI
    LOGICAL                 :: DONE
#endif
    !/
    !/ ------------------------------------------------------------------- /
    !/
#ifdef W3_S
    CALL strace (ient, 'W3SCAT')
#endif
    !
    ! 1.  Shared memory version ------------------------------------------ *
    !
#ifdef W3_SHRD
    DO isea=1, nsea
      ixy           = mapsf(isea,3)
      IF ( mapsta(ixy) .GE. 1 ) a(ispec,isea) = field(isea)
    END DO
    !
    RETURN
#endif
    !
    ! 2.  Distributed memory version ( MPI ) ----------------------------- *
    ! 2.a Initializations
    !
    ! 2.b Convert full grid to sea grid, active points only
    !
#ifdef W3_MPI
    DO isea=1, nsea
      ixy    = mapsf(isea,3)
      IF ( mapsta(ixy) .GE. 1 ) sstore(isea,ibfloc) = field(isea)
    END DO
    !
    ! 2.c Send spectral densities to appropriate remote
    !
    ioff   = 1 + (isploc-1)*nrqsg2
    IF ( nrqsg2 .GT. 0 ) CALL                                  &
         mpi_startall ( nrqsg2, irqsg2(ioff,2), ierr_mpi )
    bstat(ibfloc) = 2
    !
    ! 2.d Save locally stored results
    !
    DO jsea=1, nseal
      !!Li   ISEA   = IAPROC+(JSEA-1)*NAPROC
      isea   = min( iaproc+(jsea-1)*naproc, nsea )
      a(ispec,jsea) = sstore(isea,ibfloc)
    END DO
    !
    ! 2.e Check if any sends have finished
    !
    ib0    = ibfloc
    !
    DO j=1, mpibuf
      ib0    = 1 + mod(ib0,mpibuf)
      IF ( bstat(ib0) .EQ. 2 ) THEN
        ioff   = 1 + (bispl(ib0)-1)*nrqsg2
        IF ( nrqsg2 .GT. 0 ) THEN
          CALL mpi_testall ( nrqsg2, irqsg2(ioff,2), done,  &
               status, ierr_mpi )
        ELSE
          done   = .true.
        END IF
        IF ( done .AND. nrqsg2.GT.0 ) CALL                   &
             mpi_waitall ( nrqsg2, irqsg2(ioff,2),       &
             status, ierr_mpi )
        IF ( done ) THEN
          bstat(ib0) = 0
        END IF
      END IF
    END DO
#endif
    !
    ! 2.f Last component, finish message passing, reset buffer control
    !
#ifdef W3_MPI
    IF ( isploc .EQ. nsploc ) THEN
      DO ib0=1, mpibuf
        IF ( bstat(ib0) .EQ. 2 ) THEN
          ioff   = 1 + (bispl(ib0)-1)*nrqsg2
          IF ( nrqsg2 .GT. 0 ) CALL                        &
               mpi_waitall ( nrqsg2, irqsg2(ioff,2),       &
               status, ierr_mpi )
          bstat(ib0) = 0
        END IF
      END DO
      isploc = 0
      ibfloc = 0
    END IF
    RETURN
#endif
    !
    ! Formats
    !
    !/
    !/ End of W3SCATSMC ----------------------------------------------------- /
    !/
  END SUBROUTINE w3scatsmc
  !/
  !/ End of two new subs for SMC grid.    JGLi 15Mar2011
  !/
 
  SUBROUTINE w3smcell( IMOD, NC, IDCl, XLon, YLat )
    !! -------------------------------------------------------------------
    !!
    !!    Generated for WW3 Multi-grid boundary matching.  JGLi19Oct2020
    !!
    !! 1. Purpose:
    !
    !     Calculate the cell centre longitude and latitude in degree for a given
    !     list of cell identity or sequential numbers in the IMOD sub-grid.
    !
    !! 2. Method:
    !
    !     Regular grid SX, SY, X0, Y0 and SMC grid MRFct and IJKCel arrays
    !     in W3GDATMD are used to work out SMC grid origin and increments.
    !     Then given cell centre coordinates are calculated.  Longitude is
    !     wrapped into [0, 360) range, latitude in in (-90, 90) range.
    !     The polar cell centre is off the N-Pole to avoid singularity but
    !     its centre values are not used for propagation schemes.
    !
    !! 3. Parameters:
    !     ----------------------------------------------------------------
    !     IMOD    Int.   I   Model number to point to.
    !     NC      Int.   I   Numcer of cells to be calculated.
    !     IDCl    Int.   I   List of cell id or sequential numbers.
    !     XLon    Real   O   X-Longitude in degree of listed cells.
    !     YLat    Real   O   Y-Latitude in degree of listed cells.
    !     ----------------------------------------------------------------
    !
    !! 4. Subroutines used:
    !
    !     None
    !
    !! 5. Called by:
    !
    !     WMGLOW, W3IOPP, WMIOPP, WW3_GINT
    !
    !! 6. Error messages:
    !
    !     - Error checks on previous setting of variable.
    !
    !! 7. Remarks:
    !
    !! 8. Structure:
    !
    !! 9. Switches:
    !
    !     !/S     Enable subroutine tracing.
    !     !/T     Enable test output
    !
    ! 10. Source code:
    !
    !/ ------------------------------------------------------------------- /
    USE constants
    USE w3gdatmd
    USE w3servmd, ONLY: extcde
    USE w3odatmd, ONLY: ndse, ndst
#ifdef W3_S
    USE w3servmd, ONLY: strace
#endif
 
    IMPLICIT NONE
 
    !/ ------------------------------------------------------------------- /
    ! Input/Output variables
 
    INTEGER, INTENT(IN)               :: IMOD, NC
    INTEGER, INTENT(in), Dimension(NC):: IDCl        ! Automatic array
    REAL   , INTENT(out),Dimension(NC):: XLon, YLat
    !/ ------------------------------------------------------------------- /
    ! Local variables.
 
    REAL       :: XI0, YJ0, DXG, DYG, DX1, DY1
    INTEGER    :: I1, I3, J2, J4, MRF, ij, ijp, NSEM
#ifdef W3_S
    INTEGER    ::  IENT = 0
    CALL strace (ient, 'W3SMCELL')
#endif
 
    !! 1. Convert regular grid parameters into SMC grid origin and increments.
    dxg = grids(imod)%SX
    dyg = grids(imod)%SY
    xi0 = grids(imod)%X0 - 0.5*dxg
    yj0 = grids(imod)%Y0 - 0.5*dyg
    mrf = grids(imod)%MRFct
    dx1 = dxg/real(mrf)
    dy1 = dyg/real(mrf)
    nsem = grids(imod)%NSEA
 
    !! 2. Loop over listed cells and work out their centre coordinates.
 
#ifdef W3_OMPG
    !$OMP Parallel DO Private(ij, ijp, I1, J2, I3, J4 )
#endif
    DO ij = 1, nc
      ijp = idcl(ij)
      !!Li  Return South Pole lon-lat values for any ids < 1 or > NSEA
      !!    so these out of range points will not be matched to any cell.
      IF( ijp < 1 .OR. ijp > nsem ) THEN
        xlon(ij) = 0.0
        ylat(ij) = -90.0
      ELSE
        !!    Fetch cell array indexes from given sub-grid.
        i1=grids(imod)%IJKCel(1, ijp)
        j2=grids(imod)%IJKCel(2, ijp)
        i3=grids(imod)%IJKCel(3, ijp)
        j4=grids(imod)%IJKCel(4, ijp)
 
        !!    Calculate its cell centre lon-lat values.
        xlon(ij) = xi0 + ( float(i1) + 0.5*float(i3) )*dx1
        ylat(ij) = yj0 + ( float(j2) + 0.5*float(j4) )*dy1
      ENDIF
    END DO
#ifdef W3_OMPG
    !$OMP END Parallel DO
#endif
 
    !! 3. Wrap negative logitudes into [0, 360) range.
    WHERE( xlon < 0.0 ) xlon = xlon + 360.0
    !
    RETURN
  END SUBROUTINE w3smcell
  !!
 
  !!
  SUBROUTINE w3smcgmp( IMOD, NC, XLon, YLat, IDCl )
    !! -------------------------------------------------------------------
    !!
    !!    Generated for WW3 Multi-grid boundary matching.  JGLi22Oct2020
    !!
    !! 1. Purpose:
    !
    !     Determine whether a list of points are inside the IMOD SMC sub-grid
    !     and return the IMOD sub-grid cell indexes, if any.
    !
    !! 2. Method:
    !
    !     Convert point XLon and YLat values into cell indices i, j.
    !     Match with cell ranges (i,i+di) and (j,j+dj) to see i,j in
    !     which cell.  Return the matched cell number.  Otherwide,
    !     return an index of 0, or no matching cell found.
    !
    !! 3. Parameters:
    !     ----------------------------------------------------------------
    !     IMOD    Int.   I   Model number to point to.
    !     XLon    Real   I   X-Longitude in degree of search points.
    !     YLat    Real   I   Y-Latitude  in degree of search points.
    !     NC      Int.   I   Number of points to be searched.
    !     IDCl    Int.   O   Model number to point to.
    !     ----------------------------------------------------------------
    !
    !! 4. Subroutines used:
    !
    !     None
    !
    !! 5. Called by:
    !
    !     WMGLOW, W3IOPP, WMIOPP, WW3_GINT
    !
    !! 6. Error messages:
    !
    !     - Error checks on previous setting of variable.
    !
    !! 7. Remarks:
    !
    !! 8. Structure:
    !
    !! 9. Switches:
    !
    !     !/S     Enable subroutine tracing.
    !     !/T     Enable test output
    !
    ! 10. Source code:
    !
    !/ ------------------------------------------------------------------- /
    USE constants
    USE w3gdatmd
    USE w3servmd, ONLY: extcde
    USE w3odatmd, ONLY: ndse, ndst
#ifdef W3_S
    USE w3servmd, ONLY: strace
#endif
 
    IMPLICIT NONE
 
    !/ ------------------------------------------------------------------- /
    ! Iuput/Output variables
    INTEGER, INTENT(IN)                :: IMOD, NC
    REAL   , INTENT(in),  Dimension(NC):: XLon, YLat
    INTEGER, INTENT(out), Dimension(NC):: IDCl
    !/ ------------------------------------------------------------------- /
    ! Local variables
    INTEGER, Dimension(NC) :: IX1, JY1
    REAL       :: XI0, YJ0, DXG, DYG, DX1, DY1, XLow(NC)
    INTEGER    :: I1, I3, J2, J4, ij, ijp, MRF, NSEM, NFund
#ifdef W3_S
    INTEGER    ::  IENT = 0
    CALL strace (ient, 'W3SMCGMP')
#endif
 
    !! 1. Convert XLon YLat into SMC grid indexes in present SMC grid.
    dxg = grids(imod)%SX
    dyg = grids(imod)%SY
    xi0 = grids(imod)%X0 - 0.5*dxg
    yj0 = grids(imod)%Y0 - 0.5*dyg
    mrf = grids(imod)%MRFct
    dx1 = dxg/real(mrf)
    dy1 = dyg/real(mrf)
    nsem = grids(imod)%NSEA
 
    !!  Wrap longitude so they are great than XI0.
    xlow = xlon
    WHERE( xlow < xi0 ) xlow = xlow + 360.0
 
    !!  Convert XLon and YLat into SMC indexes.
    ix1 = floor( (xlow - xi0)/dx1 )
    jy1 = floor( (ylat - yj0)/dy1 )
 
    !!  Initialise IDCl to be all 0
    idcl = 0
 
    !! 2. Loop over all cells until all input points are found.
    nfund = 0
    ij = 0
    DO WHILE( ij < nsem .AND. nfund < nc )
      ij = ij + 1
      i1=grids(imod)%IJKCel(1, ij)
      j2=grids(imod)%IJKCel(2, ij)
      i3=grids(imod)%IJKCel(3, ij)
      j4=grids(imod)%IJKCel(4, ij)
      lpnbis: DO ijp = 1, nc
        IF( idcl(ijp) .EQ. 0 ) THEN
          !!  Check if IX1 and JY1 fall inside the cell i,j range.
          IF((ix1(ijp) .GE. i1) .AND. (ix1(ijp) .LT. i1+i3) .AND.  &
               (jy1(ijp) .GE. j2) .AND. (jy1(ijp) .LT. j2+j4)) THEN
            nfund = nfund + 1
            idcl(ijp) = ij
            EXIT  lpnbis
          ENDIF
        ENDIF
      END DO  lpnbis
    END DO
 
    !!  If any IDCl element remians to be 0, it means no cell is found
    !!  covering this point.  So check IDCl(ij) > 0 to ensure in grid.
 
    RETURN
 
  END SUBROUTINE w3smcgmp
  !!
 
  !/ End of module W3PSMCMD -------------------------------------------- /
  !/
END MODULE w3psmcmd