w3pro3md.F90

Go to the documentation of this file.

 
 
#include "w3macros.h"
!/ ------------------------------------------------------------------- /
MODULE w3pro3md
  !/
  !/                  +-----------------------------------+
  !/                  | WAVEWATCH III           NOAA/NCEP |
  !/                  |           H. L. Tolman            |
  !/                  |                        FORTRAN 90 |
  !/                  | Last update :         27-May-2014 |
  !/                  +-----------------------------------+
  !/
  !/    27-Feb-2000 : Origination.                        ( version 2.08 )
  !/    17-Sep-2000 : Clean-up.                           ( version 2.13 )
  !/    10-Dec-2001 : Sub-grid obstructions.              ( version 2.14 )
  !/    16-Oct-2002 : Change INTENT for ATRN in W3XYP3.   ( version 3.00 )
  !/    26-Dec-2002 : Moving grid version.                ( version 3.02 )
  !/    01-Aug-2003 : Moving grid GSE correction.         ( version 3.03 )
  !/    17-Dec-2004 : Multiple grid version.              ( version 3.06 )
  !/    07-Sep-2005 : Upgrade XY boundary conditions.     ( version 3.08 )
  !/    09-Nov-2005 : Removing soft boundary option.      ( version 3.08 )
  !/    05-Mar-2008 : Added NEC sxf90 compiler directives.
  !/                  (Chris Bunney, UK Met Office)       ( version 3.13 )
  !/    01-Apr-2008 : Bug fix W3MAP3 MAPSTA range check.  ( version 3.13 )
  !/    29-May-2009 : Preparing distribution version.     ( version 3.14 )
  !/    30-Oct-2009 : Implement curvilinear grid type.    ( version 3.14 )
  !/                  (W. E. Rogers & T. J. Campbell, NRL)
  !/    17-Aug-2010 : Add test output W3XYP3.           ( version 3.14.5 )
  !/    06-Dec-2010 : Change from GLOBAL (logical) to ICLOSE (integer) to
  !/                  specify index closure for a grid.   ( version 3.14 )
  !/                  (T. J. Campbell, NRL)
  !/    26-Dec-2012 : More initializations.               ( version 4.11 )
  !/    01-Jul-2013 : Adding UQ and UNO switches to chose between third
  !/                  and second order schemes.           ( version 4.12 )
  !/    12-Sep-2013 : Add documentation for global clos.  ( version 4.12 )
  !/    27-May-2014 : Adding OMPH switch.                 ( version 5.02 )
  !/
  !/    Copyright 2009-2014 National Weather Service (NWS),
  !/       National Oceanic and Atmospheric Administration.  All rights
  !/       reserved.  WAVEWATCH III is a trademark of the NWS.
  !/       No unauthorized use without permission.
  !/
  !  1. Purpose :
  !
  !     Bundles routines for third order propagation scheme in single
  !     module.
  !
  !  2. Variables and types :
  !
  !      Name      Type  Scope    Description
  !     ----------------------------------------------------------------
  !      TRNMIN    R.P.  Private   Minimum transparancy for local
  !                                switching off of averaging.
  !     ----------------------------------------------------------------
  !
  !     Also work arrays for W3KTP3 (private).
  !
  !  3. Subroutines and functions :
  !
  !      Name      Type  Scope    Description
  !     ----------------------------------------------------------------
  !      W3MAP3    Subr. Public   Set up auxiliary maps.
  !      W3MAPT    Subr. Public   Set up transparency map for GSE.
  !      W3XYP3    Subr. Public   Third order spatial propagation.
  !      W3KTP3    Subr. Public   Third order spectral propagation.
  !     ----------------------------------------------------------------
  !
  !  4. Subroutines and functions used :
  !
  !      Name      Type  Module   Description
  !     ----------------------------------------------------------------
  !      STRACE    Subr. W3SERVMD Subroutine tracing.
  !      W3QCK1    Subr. W3UQCKMD Regular grid UQ scheme.
  !      W3QCK2    Subr.   Id.    Irregular grid UQ scheme.
  !      W3QCK3    Subr.   Id.    Regular grid UQ scheme + obstructions.
  !      W3UNO2    Subr. W3UNO2MD UNO2 scheme for irregular grid.
  !      W3UNO2r   Subr.   Id.    UNO2 scheme reduced to regular grid.
  !      W3UNO2s   Subr.   Id.    UNO2 regular grid with subgrid
  !                               obstruction.
  !     ----------------------------------------------------------------
  !
  !  5. Remarks :
  !
  !     - The averaging is not performed around semi-transparent grid
  !       points to avoid that leaking through barriers occurs.
  !
  !  6. Switches :
  !
  !       !/UQ    3rd order UQ propagation scheme.
  !       !/UNO   2nd order UNO propagation scheme.
  !
  !       !/MGP   Moving grid corrections.
  !       !/MGG   Moving grid corrections.
  !
  !       !/OMPH  Hybrid OpenMP directives.
  !
  !       !/S     Enable subroutine tracing.
  !       !/Tn    Enable test output.
  !
  !  7. Source code :
  !
  !/ ------------------------------------------------------------------- /
  !/
  !/ Public variables
  !/
  PUBLIC
  !/
  !/ Private data
  !/
  REAL, PRIVATE, PARAMETER:: TRNMIN = 0.95
  !/
CONTAINS
  !/ ------------------------------------------------------------------- /
  SUBROUTINE w3map3
    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III           NOAA/NCEP |
    !/                  |           H. L. Tolman            |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         01-Apr-2008 |
    !/                  +-----------------------------------+
    !/
    !/    27-Feb-2000 : Origination.                        ( version 2.08 )
    !/    10-Dec-2001 : Sub-grid obstructions.              ( version 2.14 )
    !/                  (array allocation only.)
    !/    17-Dec-2004 : Multiple grid version.              ( version 3.06 )
    !/    09-Nov-2005 : Removing soft boundary option.      ( version 3.08 )
    !/    01-Apr-2008 : Bug fix sec. 4 MAPSTA range check.  ( version 3.13 )
    !/
    !  1. Purpose :
    !
    !     Generate 'map' arrays for the ULTIMATE QUICKEST scheme.
    !
    !  2. Method :
    !
    !     MAPX2, MAPY2, MAPTH2 and MAPWN2 contain consecutive 1-D spatial
    !     grid counters (e.g., IXY = (IX-1)*MY + IY). The arrays are
    !     devided in three parts. For MAPX2, these ranges are :
    !
    !         1    - NMX0  Counters IXY for which grid point (IX,IY) and
    !                      (IX+1,IY) both are active grid points.
    !       NMX0+1 - NMX1  Id. only (IX,IY) active.
    !       NMX1+1 - NMX2  Id. only (IX+1,IY) active.
    !
    !     The array MAPY2 has a similar layout varying IY instead of IX.
    !
    !  3. Parameters :
    !
    !     Parameter list
    !     ----------------------------------------------------------------
    !     ----------------------------------------------------------------
    !
    !  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.   Map MAPX2
    !        a  Range 1 to NMX0
    !        b  Range NMX0+1 to NMX1
    !        c  Range NMX1+1 to NMX2
    !      2.   Map MAPY2
    !        a  Range 1 to NMY0
    !        b  Range NMY0+1 to NMY1
    !        c  Range NMY1+1 to NMY2
    !      3.   Map MAPAXY
    !      4.   Map MAPCXY
    !      5.   Maps for intra-spectral propagation
    !        a  MAPTH2, MAPATK
    !        b  MAPWN2
    !     ------------------------------------------------------
    !
    !  9. Switches :
    !
    !     !/S   Enable subroutine tracing.
    !     !/T   Enable test output.
    !
    ! 10. Source code :
    !/ ------------------------------------------------------------------- /
    USE w3gdatmd, ONLY: nk, nth, nspec, nx, ny, nsea, mapsta, mapsf,&
         gtype
    USE w3adatmd, ONLY: nmx0, nmx1, nmx2, nmy0, nmy1, nmy2, nact,   &
         ncent, mapx2, mapy2, mapaxy, mapcxy,        &
         mapth2, mapwn2
    USE w3odatmd, ONLY: ndst
#ifdef W3_S
    USE w3servmd, ONLY: strace
#endif
    !/
    IMPLICIT NONE
    !/
    !/ ------------------------------------------------------------------- /
    !/ Parameter list
    !/
    !/ ------------------------------------------------------------------- /
    !/ Local parameters
    !/
    INTEGER                 :: IX, IY, IXY0, IX2, IY2, IX0, IY0,    &
         ISEA, IK, ITH, ISP, ISP0, ISP2, NCENTC
#ifdef W3_S
    INTEGER, SAVE           :: IENT = 0
#endif
#ifdef W3_T
    INTEGER                 :: MAPTXY(NY,NX), I, IXY
    INTEGER                 :: MAPTST(NK+2,NTH)
#endif
    !/
    !/ ------------------------------------------------------------------- /
    !/
#ifdef W3_S
    CALL strace (ient, 'W3MAP3')
#endif
    !
    IF (gtype .LT. 3) THEN
      ! 1.  Map MAPX2 ------------------------------------------------------ *
      ! 1.a Range 1 to NMX0
      !
#ifdef W3_T
      maptxy = 0.
#endif
      !
      nmx0   = 0
      DO ix=1, nx
        ixy0   = (ix-1)*ny
        ix2    = 1 + mod(ix,nx)
        DO iy=2, ny-1
          IF ( mapsta(iy,ix).EQ.1 .AND. mapsta(iy,ix2).EQ.1 ) THEN
            nmx0   = nmx0 + 1
            mapx2(nmx0) = ixy0 + iy
#ifdef W3_T
            maptxy(iy,ix) = maptxy(iy,ix) + 1
#endif
          END IF
        END DO
      END DO
      !
      ! 1.b Range NMX0+1 to NMX1
      !
      nmx1   = nmx0
      DO ix=1, nx
        ixy0   = (ix-1)*ny
        ix2    = 1 + mod(ix,nx)
        DO iy=2, ny-1
          IF ( mapsta(iy,ix).EQ.1 .AND. mapsta(iy,ix2).NE.1 ) THEN
            nmx1   = nmx1 + 1
            mapx2(nmx1) = ixy0 + iy
#ifdef W3_T
            maptxy(iy,ix) = maptxy(iy,ix) + 2
#endif
          END IF
        END DO
      END DO
      !
      ! 1.c Range NMX1+1 to NMX2
      !
      nmx2   = nmx1
      DO ix=1, nx
        ixy0   = (ix-1)*ny
        ix2    = 1 + mod(ix,nx)
        DO iy=2, ny-1
          IF ( mapsta(iy,ix).NE.1 .AND. mapsta(iy,ix2).EQ.1 ) THEN
            nmx2   = nmx2 + 1
            mapx2(nmx2) = ixy0 + iy
#ifdef W3_T
            maptxy(iy,ix) = maptxy(iy,ix) + 4
#endif
          END IF
        END DO
      END DO
      !
#ifdef W3_T
      WRITE (ndst,9000) 'MAPX2', nmx0, nmx1-nmx0,                     &
           nmx2-nmx1, nmx2
      DO iy=ny, 1, -1
        WRITE (ndst,9001) (maptxy(iy,ix),ix=1, nx)
      END DO
#endif
      !
      ! 2.  Map MAPY2 ------------------------------------------------------ *
      ! 2.a Range 1 to NMY0
      !
#ifdef W3_T
      maptxy = 0.
#endif
      !
      nmy0   = 0
      DO ix=1, nx
        ixy0   = (ix-1)*ny
        DO iy=1, ny-1
          iy2    = iy + 1
          IF ( mapsta(iy,ix).EQ.1 .AND. mapsta(iy2,ix).EQ.1 ) THEN
            nmy0   = nmy0 + 1
            mapy2(nmy0) = ixy0 + iy
#ifdef W3_T
            maptxy(iy,ix) = maptxy(iy,ix) + 1
#endif
          END IF
        END DO
      END DO
      !
      ! 2.b Range NMY0+1 to NMY1
      !
      nmy1   = nmy0
      DO ix=1, nx
        ixy0   = (ix-1)*ny
        DO iy=1, ny-1
          iy2    = iy + 1
          IF ( mapsta(iy,ix).EQ.1 .AND. mapsta(iy2,ix).NE.1 ) THEN
            nmy1   = nmy1 + 1
            mapy2(nmy1) = ixy0 + iy
#ifdef W3_T
            maptxy(iy,ix) = maptxy(iy,ix) + 2
#endif
          END IF
        END DO
      END DO
      !
      ! 2.c Range NMY1+1 to NMY2
      !
      nmy2   = nmy1
      DO ix=1, nx
        ixy0   = (ix-1)*ny
        DO iy=1, ny-1
          iy2    = iy + 1
          IF ( mapsta(iy,ix).NE.1 .AND. mapsta(iy2,ix).EQ.1 ) THEN
            nmy2   = nmy2 + 1
            mapy2(nmy2) = ixy0 + iy
#ifdef W3_T
            maptxy(iy,ix) = maptxy(iy,ix) + 4
#endif
          END IF
        END DO
      END DO
      !
#ifdef W3_T
      WRITE (ndst,9000) 'MAPY2', nmy0, nmy1-nmy0,                     &
           nmy2-nmy1, nmy2
      DO iy=ny, 1, -1
        WRITE (ndst,9001) (maptxy(iy,ix),ix=1, nx)
      END DO
#endif
      !
      ! 3.  Map MAPAXY ----------------------------------------------------- *
      !
      nact   = 0
      DO ix=1, nx
        iy0    = (ix-1)*ny
        DO iy=2, ny-1
          IF ( mapsta(iy,ix).EQ.1 ) THEN
            nact         = nact + 1
            mapaxy(nact) = iy0 + iy
          END IF
        END DO
      END DO
      !
      ! 4.  Map MAPCXY ----------------------------------------------------- *
      !
      ncent  = 0
      ncentc = nsea
      mapcxy = 0
      !
      DO isea=1,  nsea
        ix      = mapsf(isea,1)
        ix0    = ix-1
        ix2    = ix+1
        iy      = mapsf(isea,2)
        iy0    = iy-1
        iy2    = iy+1
        IF ( ix .EQ. nx ) ix2 = 1
        IF ( ix .EQ. 1 ) ix0 = nx
        IF ( mapsta(iy,ix).EQ.2 .OR. mapsta(iy,ix).LT.0 ) THEN
          mapcxy(ncentc) = isea
          ncentc = ncentc - 1
        ELSE IF ( mapsta(iy0,ix0).GE.1 .AND.                     &
             mapsta(iy0,ix ).GE.1 .AND.                     &
             mapsta(iy0,ix2).GE.1 .AND.                     &
             mapsta(iy ,ix0).GE.1 .AND.                     &
             mapsta(iy ,ix2).GE.1 .AND.                     &
             mapsta(iy2,ix0).GE.1 .AND.                     &
             mapsta(iy2,ix ).GE.1 .AND.                     &
             mapsta(iy2,ix2).GE.1 ) THEN
          ncent  = ncent + 1
          mapcxy(ncent) = isea
        ELSE
          mapcxy(ncentc) = isea
          ncentc = ncentc - 1
        END IF
      END DO
    END IF
    !
    ! 5.  Maps for intra-spectral propagation ---------------------------- *
    !
    IF ( mapth2(1) .NE. 0 ) RETURN
    !
#ifdef W3_T
    maptst = 0
#endif
    !
    ! 5.a MAPTH2 and MAPBTK
    !
    DO ik=1, nk
      DO ith=1, nth
        isp    = ith + (ik-1)*nth
        isp2   = (ik+1) + (ith-1)*(nk+2)
        mapth2(isp) = isp2
#ifdef W3_T
        maptst(ik+1,ith) = maptst(ik+1,ith) + 1
#endif
      END DO
    END DO
    !
#ifdef W3_T
    WRITE (ndst,9000) 'MAPTH2', isp, 0, 0, isp
    DO ik=nk+2, 1, -1
      WRITE (ndst,9001) (maptst(ik,ith),ith=1, nth)
    END DO
    maptst = 0
#endif
    !
    ! 5.b MAPWN2
    !
    isp0   = 0
    DO ik=1, nk-1
      DO ith=1, nth
        isp0   = isp0 + 1
        isp2   = (ik+1) + (ith-1)*(nk+2)
        mapwn2(isp0) = isp2
#ifdef W3_T
        maptst(ik+1,ith) = maptst(ik+1,ith) + 1
#endif
      END DO
    END DO
    !
    DO ith=1, nth
      isp0   = isp0 + 1
      isp2   = nk+1 + (ith-1)*(nk+2)
      mapwn2(isp0) = isp2
#ifdef W3_T
      maptst(nk+1,ith) = maptst(nk+1,ith) + 2
#endif
    END DO
    !
    DO ith=1, nth
      isp0   = isp0 + 1
      isp2   = 1 + (ith-1)*(nk+2)
      mapwn2(isp0) = isp2
#ifdef W3_T
      maptst(1,ith) = maptst(1,ith) + 4
#endif
    END DO
    !
#ifdef W3_T
    WRITE (ndst,9000) 'MAPWN2', nspec-nth, nth, nth, nspec+nth
    DO ik=nk+2, 1, -1
      WRITE (ndst,9001) (maptst(ik,ith),ith=1, nth)
    END DO
#endif
    !
    RETURN
    !
    ! Formats
    !
#ifdef W3_T
9000 FORMAT (/' TEST W3MAP3 : TEST MAP FOR PROPAGATION'/          &
         '      MAP     : ',a/                               &
         '      CENTRAL : ',i6/                              &
         '      ABOVE   : ',i6/                              &
         '      BELOW   : ',i6/                              &
         '      TOTAL   : ',i6/)
9001 FORMAT (1x,130i1)
9010 FORMAT (' TEST W3MAP3 : COMPOSITE MAPS TH2, WN2 AND BTK')
9011 FORMAT (2x,60i2)
#endif
    !/
    !/ End of W3MAP3 ----------------------------------------------------- /
    !/
  END SUBROUTINE w3map3
  !/ ------------------------------------------------------------------- /
  SUBROUTINE w3mapt
    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III           NOAA/NCEP |
    !/                  |           H. L. Tolman            |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         17-Dec-2004 |
    !/                  +-----------------------------------+
    !/
    !/    10-Dec-2001 : Origination.                        ( version 2.14 )
    !/    17-Dec-2004 : Multiple grid version.              ( version 3.06 )
    !/
    !  1. Purpose :
    !
    !     Generate 'map' arrays for the ULTIMATE QUICKEST scheme to combine
    !     GSE alleviation with obstructions.
    !
    !  2. Method :
    !
    !  3. Parameters :
    !
    !     Parameter list
    !     ----------------------------------------------------------------
    !     ----------------------------------------------------------------
    !
    !  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 :
    !
    !     See source code.
    !
    !  9. Switches :
    !
    !     !/S   Enable subroutine tracing.
    !
    ! 10. Source code :
    !/ ------------------------------------------------------------------- /
    USE w3gdatmd, ONLY: nx, ny, nsea, mapsf
    USE w3adatmd, ONLY: atrnx, atrny, maptrn
#ifdef W3_S
    USE w3servmd, ONLY: strace
#endif
    !/
    IMPLICIT NONE
    !/
    !/ ------------------------------------------------------------------- /
    !/ Parameter list
    !/
    !/ ------------------------------------------------------------------- /
    !/ Local parameters
    !/
    INTEGER                 :: ISEA, IXY
#ifdef W3_S
    INTEGER, SAVE           :: IENT = 0
#endif
    !/
    !/ ------------------------------------------------------------------- /
    !/
#ifdef W3_S
    CALL strace (ient, 'W3MAPT')
#endif
    !
    ! 1.  Map MAPTRN ----------------------------------------------------- *
    !
    DO isea=1, nsea
      ixy         = mapsf(isea,3)
 
      !notes: Oct 22 2012: I changed this because it *looks* like a bug.
      !       I have not confirmed that it is a bug.
      !       Old code is given (2 lines). Only the first line is
      !       changed.
 
      !old    MAPTRN(IXY) = MIN( ATRNX(IXY,1) ,ATRNY(IXY,-1) ,              &
      !old                       ATRNY(IXY,1), ATRNY(IXY,-1) ) .LT. TRNMIN
 
      maptrn(ixy) = min( atrnx(ixy,1) ,atrnx(ixy,-1) ,              &
           atrny(ixy,1), atrny(ixy,-1) ) .LT. trnmin
    END DO
    !
    RETURN
    !
    ! Formats
    !/
    !/ End of W3MAPT ----------------------------------------------------- /
    !/
  END SUBROUTINE w3mapt
  !/ ------------------------------------------------------------------- /
  SUBROUTINE w3xyp3 ( ISP, DTG, MAPSTA, MAPFS, VQ, VGX, VGY )
    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III           NOAA/NCEP |
    !/                  |           H. L. Tolman            |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         27-May-2014 |
    !/                  +-----------------------------------+
    !/
    !/    27-Feb-2000 : Origination.                        ( version 2.08 )
    !/    17-Sep-2000 : Clean-up.                           ( version 2.13 )
    !/    10-Dec-2001 : Sub-grid obstructions.              ( version 2.14 )
    !/    16-Oct-2002 : Change INTENT for ATRNRX/Y.         ( version 3.00 )
    !/    26-Dec-2002 : Moving grid version.                ( version 3.02 )
    !/    01-Aug-2003 : Moving grid GSE correction.         ( version 3.03 )
    !/    17-Dec-2004 : Multiple grid version.              ( version 3.06 )
    !/    07-Sep-2005 : Upgrade XY boundary conditions.     ( version 3.08 )
    !/    09-Nov-2005 : Removing soft boundary option.      ( version 3.08 )
    !/    05-Mar-2008 : Added NEC sxf90 compiler directives.
    !/                  (Chris Bunney, UK Met Office)       ( version 3.13 )
    !/    30-Oct-2009 : Implement curvilinear grid type.    ( version 3.14 )
    !/                  (W. E. Rogers & T. J. Campbell, NRL)
    !/    17-Aug-2010 : Add test output.                  ( version 3.14.5 )
    !/    30-Oct-2010 : Implement unstructured grid         ( version 3.14 )
    !/    06-Dec-2010 : Change from GLOBAL (logical) to ICLOSE (integer) to
    !/                  specify index closure for a grid.   ( version 3.14 )
    !/                  (T. J. Campbell, NRL)
    !/    01-Jul-2013 : Adding UQ and UNO switches to chose between third
    !/                  and second order schemes.           ( version 4.12 )
    !/    12-Sep-2013 : Add documentation for global clos.  ( version 4.12 )
    !/    27-May-2014 : Adding OMPH switch.                 ( version 5.02 )
    !/
    !  1. Purpose :
    !
    !     Propagation in phyiscal space for a given spectral component.
    !
    !  2. Method :
    !
    !     Third-order ULTIMATE QUICKEST scheme with averaging.
    !     Curvilinear grid implementation: Fluxes are computed in index space
    !       and then transformed back into physical space.  The diffusion term
    !       is handled in physical space.  The averaging scheme is adapted for
    !       curvilinear grids by applying the appropriate local rotations and
    !       adjustments to interpolation weights which control the strength of
    !       the averaging in axial directions.
    !
    !  3. Parameters :
    !
    !     Parameter list
    !     ----------------------------------------------------------------
    !       ISP     Int.   I   Number of spectral bin (IK-1)*NTH+ITH
    !       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.
    !       VGX/Y   Real   I   Speed of grid.
    !     ----------------------------------------------------------------
    !
    !     Local variables.
    !     ----------------------------------------------------------------
    !       NTLOC   Int   Number of local time steps.
    !       DTLOC   Real  Local propagation time step.
    !       VCFL0X  R.A.  Local courant numbers for absolute group vel.
    !                     using local X-grid step.
    !       VCFL0Y  R.A.  Id. in Y.
    !     ----------------------------------------------------------------
    !
    !  4. Subroutines used :
    !
    !       W3QCK3   Actual propagation algorithm
    !
    !       STRACE   Service routine.
    !
    !  5. Called by :
    !
    !       W3WAVE   Wave model routine.
    !
    !  6. Error messages :
    !
    !       None.
    !
    !  7. Remarks :
    !
    !     - Note that the ULTIMATE limiter does not guarantee non-zero
    !       energies.
    !     - The present scheme shows a strong distorsion when propaga-
    !       ting a field under an angle with the grid in a truly 2-D
    !       fashion. Propagation is therefore split along the two
    !       axes.
    !     - Two boundary treatments are available. The first uses real
    !       boundaries in each space. In this case waves will not
    !       penetrate in narrow straights under an angle with the grid.
    !       This behavior is improved by using a 'soft' option, in
    !       which the 'X' or 'Y' sweep allows for energy to go onto
    !       the land. This improves the above behavior, but implies
    !       that X-Y connenctions are required in barriers for them
    !       to become inpenetrable.
    !     - Note that unlike in W3XYP2, isotropic diffusion is never
    !       used for growth.
    !     - Curvilinear grid implementation. Variables FACX, FACY, CCOS, CSIN,
    !       CCURX, CCURY are not needed and have been removed.  FACX is accounted
    !       for as approriate in this subroutine.  FACX is also accounted for in
    !       the case of .NOT.FLCX.  Since FACX is removed, there is now a check for
    !       .NOT.FLCX in this subroutine.  In CFL calcs dx and dy are omitted,
    !       since dx=dy=1 in index space.  Curvilinear grid derivatives
    !       (DPDY, DQDX, etc.) and metric (GSQRT) are brought in via W3GDATMD.
    !     - The strength of the averaging scheme is dependent on grid resolution.
    !       Since grid resolution is non-uniform for curvilinear grids, this means
    !       that the strength of the averaging is also non-uniform. This may not be
    !       a desirable effect. A potential future upgrade would be to add an
    !       additional term/factor that balances the effect of the spatial
    !       variation of grid resolution.
    !
    !  8. Structure :
    !
    !     ---------------------------------------------
    !       1.  Preparations
    !         a Set constants
    !         b Initialize arrays
    !       2.  Prepare arrays
    !         a Velocities and 'Q'
    !       3.  Loop over sub-steps
    !       ----------------------------------------
    !         a Average
    !         b Propagate
    !         c Update boundaries.
    !       ----------------------------------------
    !       4.  Store Q field in spectra
    !     ---------------------------------------------
    !
    !  9. Switches :
    !
    !       !/S     Enable subroutine tracing.
    !
    !       !/OMPH  Hybrid OpenMP directives.
    !
    !       !/MGP   Moving grid corrections.
    !       !/MGG   Moving grid corrections.
    !
    !       !/T     Enable general test output.
    !       !/T0    Dump of precalcaulted interpolation data.
    !       !/T1    Dump of input field and fluxes.
    !       !/T2    Dump of output field (before boundary update).
    !       !/T3    Dump of output field (final).
    !
    ! 10. Source code :
    !
    !/ ------------------------------------------------------------------- /
    USE constants
    !
    USE w3timemd, ONLY: dsec21
    !
    USE w3gdatmd, ONLY: nx, ny, nsea, mapsf, dtcfl, clats,      &
         iclose, flcx, flcy, nk, nth, dth, xfr,  &
         iclose_none, iclose_smpl, iclose_trpl,  &
         ecos, esin, sig, wdcg, wdth, pfmove,    &
         flagll, dpdx, dpdy, dqdx, dqdy, gsqrt
    USE w3wdatmd, ONLY: time
    USE w3adatmd, ONLY: nmx0, nmx1, nmx2, nmy0, nmy1, nmy2, nact,   &
         ncent, mapx2, mapy2, mapaxy, mapcxy,        &
         maptrn, cg, cx, cy, atrnx, atrny, itime
    USE w3idatmd, ONLY: flcur
    USE w3odatmd, ONLY: ndse, ndst, flbpi, nbi, tbpi0, tbpin,       &
         isbpi, bbpi0, bbpin, iaproc, naperr
    USE w3servmd, ONLY: extcde
#ifdef W3_S
    USE w3servmd, ONLY: strace
#endif
#ifdef W3_UQ
    USE w3uqckmd
#endif
#ifdef W3_UNO
    USE w3uno2md
#endif
    !/
    IMPLICIT NONE
    !/
    !/ ------------------------------------------------------------------- /
    !/ Parameter list
    !/
    INTEGER, INTENT(IN)     :: ISP, MAPSTA(NY*NX), MAPFS(NY*NX)
    REAL, INTENT(IN)        :: DTG, VGX, VGY
    REAL, INTENT(INOUT)     :: VQ(1-NY:NY*(NX+2))
    !/
    !/ ------------------------------------------------------------------- /
    !/ Local parameters
    !/
    INTEGER                 :: ITH, IK, NTLOC, ITLOC, ISEA, IXY, IP
    INTEGER                 :: IX, IY, IXC, IYC, IBI
    INTEGER                 :: IIXY1(NSEA), IIXY2(NSEA),            &
         IIXY3(NSEA), IIXY4(NSEA)
    INTEGER                 :: TTEST(2),DTTST
#ifdef W3_S
    INTEGER, SAVE           :: IENT = 0
#endif
    REAL                    :: CG0, CGA, CGN, CGX, CGY, CXC, CYC,   &
         CXMIN, CXMAX, CYMIN, CYMAX
    REAL                    :: CGC, FGSE = 1.
    REAL                    :: FTH, FTHX, FTHY, FCG, FCGX, FCGY
    REAL                    :: DTLOC, DTRAD,                        &
         DXCGN, DYCGN, DXCGS, DYCGS, DXCGC,   &
         DYCGC
    REAL                    :: RDI1(NSEA), RDI2(NSEA),              &
         RDI3(NSEA), RDI4(NSEA)
    REAL                    :: TMPX, TMPY, RD1, RD2, RD3, RD4
    LOGICAL                 :: YFIRST
    LOGICAL                 :: GLOBAL
    REAL                    :: CP, CQ
    !/
    !/ Automatic work arrays
    !/
    INTEGER                 :: MAPSTX(1-2*NY:NY*(NX+2))
    REAL                    :: VLCFLX((NX+1)*NY), VLCFLY((NX+1)*NY),&
         AQ(1-NY:NY*(NX+2))
    REAL                    :: CXTOT((NX+1)*NY), CYTOT(NX*NY)
    !/
    !/ ------------------------------------------------------------------- /
    !/
#ifdef W3_S
    CALL strace (ient, 'W3XYP3')
#endif
    !
    ! 1.  Preparations --------------------------------------------------- *
 
    IF ( iclose .EQ. iclose_trpl ) THEN
      !/ ------------------------------------------------------------------- /
      IF (iaproc .EQ. naperr) &
           WRITE(ndse,*)'SUBROUTINE W3XYP3 IS NOT YET ADAPTED FOR '//    &
           'TRIPOLE GRIDS. STOPPING NOW.'
      CALL extcde ( 1 )
    END IF
 
    ! 1.a Set constants
    !
    global = iclose.NE.iclose_none
    ith    = 1 + mod(isp-1,nth)
    ik     = 1 + (isp-1)/nth
    !
    cg0    = 0.575 * grav / sig(1)
    cga    = 0.575 * grav / sig(ik)
    cgx    = cga * ecos(ith)
    cgy    = cga * esin(ith)
#ifdef W3_MGP
    cgx    = cgx - vgx
    cgy    = cgy - vgy
#endif
    cgc    = sqrt( cgx**2 + cgy**2 )
#ifdef W3_MGG
    fgse   = ( cga / max(0.001*cga,cgc) )**pfmove
#endif
    !
    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) )
#ifdef W3_MGP
      cxc    = max( abs(cxmin-vgx) , abs(cxmax-vgx) )
      cyc    = max( abs(cymin-vgy) , abs(cymax-vgy) )
#endif
    ELSE
      cxc    = 0.
      cyc    = 0.
    END IF
    !
    cgn    = max( abs(cgx) , abs(cgy) , cxc, cyc, 0.001*cg0 )
    !
    ntloc  = 1 + int(dtg/(dtcfl*cg0/cgn))
    dtloc  = dtg / real(ntloc)
    dtrad  = dtloc
    IF ( flagll ) dtrad=dtrad/(dera*radius)
    !
    ttest(1) = time(1)
    ttest(2) = 0
    dttst = dsec21(ttest,time)
    yfirst = mod(nint(dttst/dtg),2) .EQ. 0
    !
#ifdef W3_T
    WRITE (ndst,9000) yfirst
    WRITE (ndst,9001) isp, ith, ik, ecos(ith), esin(ith)
#endif
    !
    ! 1.b Initialize arrays
    !
#ifdef W3_T
    WRITE (ndst,9010)
#endif
    !
    vlcflx = 0.
    vlcfly = 0.
    cxtot  = 0.
    cytot  = 0.
    !
    mapstx(1:nx*ny) = mapsta(1:nx*ny)
    !
    IF ( global ) THEN
      mapstx(1-2*ny:0)            = mapsta((nx-2)*ny+1:nx*ny)
      mapstx(nx*ny+1:nx*ny+2*ny)  = mapsta(1:2*ny)
    ELSE
      mapstx(1-2*ny:0)            = 0
      mapstx(nx*ny+1:nx*ny+2*ny)  = 0
    END IF
    !
    ! 1.c Pre-calculate interpolation info
    !
    fth = fgse * wdth * dth * dtloc
    fcg = fgse * wdcg * 0.5 * (xfr-1./xfr) * dtloc
    IF ( flagll ) THEN
      fth = fth / radius / dera
      fcg = fcg / radius / dera
    END IF
    fcg = fcg / real(ntloc) !TJC: required to match original (is this correct?)
    fthx = - fth * esin(ith)
    fthy =   fth * ecos(ith)
    fcgx =   fcg * ecos(ith)
    fcgy =   fcg * esin(ith)
    !
#ifdef W3_T0
    WRITE (ndst,9011)
#endif
    !
#ifdef W3_OMPH
    !$OMP PARALLEL DO PRIVATE (ISEA, IX, IY, TMPX, TMPY, &
    !$OMP&                     DXCGN, DYCGN, DXCGS, DYCGS, DXCGC, DYCGC, &
    !$OMP&                     IXC, IYC)
#endif
    !
    DO isea=1, nsea
      !
      ix  = mapsf(isea,1)
      iy  = mapsf(isea,2)
      !
      ! 1.c.1 Normal and parallel width ...
      !
      tmpx   = fthx / clats(isea)
      tmpy   = fthy
      dxcgn  = dpdx(iy,ix)*tmpx + dpdy(iy,ix)*tmpy
      dycgn  = dqdx(iy,ix)*tmpx + dqdy(iy,ix)*tmpy
      tmpx   = fcgx / clats(isea)
      tmpy   = fcgy
      dxcgs  = dpdx(iy,ix)*tmpx + dpdy(iy,ix)*tmpy
      dycgs  = dqdx(iy,ix)*tmpx + dqdy(iy,ix)*tmpy
      !
      ! 1.c.2 "Sum" corner (and mirror image) ...
      !
      dxcgc  = dxcgn + dxcgs
      dycgc  = dycgn + dycgs
      !
      ixc    = ny
      IF ( dxcgc .LT. 0. ) ixc = - ixc
      iyc    = 1
      IF ( dycgc .LT. 0. ) iyc = - iyc
      !
      iixy1(isea) = ixc + iyc
      IF ( abs(dxcgc) .GT. abs(dycgc) ) THEN
        iixy2(isea) = ixc
        rdi1(isea) = abs(dycgc/dxcgc)
        rdi2(isea) = abs(dxcgc)
      ELSE
        iixy2(isea) = iyc
        IF ( abs(dycgc) .GT. 1.e-5 ) THEN
          rdi1(isea) = abs(dxcgc/dycgc)
        ELSE
          rdi1(isea) = 1.
        END IF
        rdi2(isea) = abs(dycgc)
      END IF
      !
#ifdef W3_T0
      WRITE (ndst,9012) isea, ith, iixy1(isea), iixy2(isea), &
           rdi1(isea), rdi2(isea)*cg(ik,1)
#endif
      !
      ! 1.c.2 "Difference" corner (and mirror image) ...
      !
      dxcgc  = dxcgn - dxcgs
      dycgc  = dycgn - dycgs
      !
      ixc    = ny
      IF ( dxcgc .LT. 0. ) ixc = - ixc
      iyc    = 1
      IF ( dycgc .LT. 0. ) iyc = - iyc
      !
      iixy3(isea) = ixc + iyc
      IF ( abs(dxcgc) .GT. abs(dycgc) ) THEN
        iixy4(isea) = ixc
        rdi3(isea) = abs(dycgc/dxcgc)
        rdi4(isea) = abs(dxcgc)
      ELSE
        iixy4(isea) = iyc
        IF ( abs(dycgc) .GT. 1.e-5 ) THEN
          rdi3(isea) = abs(dxcgc/dycgc)
        ELSE
          rdi3(isea) = 1.
        END IF
        rdi4(isea) = abs(dycgc)
      END IF
      !
#ifdef W3_T0
      WRITE (ndst,9013) iixy3(isea), iixy4(isea), rdi3(isea), &
           rdi4(isea)*cg(ik,1)
#endif
      !
    END DO
    !
#ifdef W3_OMPH
    !$OMP END PARALLEL DO
#endif
    !
    ! 2.  Calculate velocities and diffusion coefficients ---------------- *
    ! 2.a Velocities
    !
    !     Q     = ( A / CG * CLATS )
    !     LCFLX = ( COS*CG / CLATS ) * DT / DX
    !     LCFLY = (     SIN*CG )     * DT / DY
    !
#ifdef W3_T
    WRITE (ndst,9020) nsea
#endif
    !
#ifdef W3_OMPH
    !$OMP PARALLEL DO PRIVATE (ISEA, IXY)
#endif
    !
    DO isea=1, nsea
      ixy         = mapsf(isea,3)
      vq(ixy) = vq(ixy) / cg(ik,isea) * clats(isea)
      cxtot(ixy) = ecos(ith) * cg(ik,isea) / clats(isea)
      cytot(ixy) = esin(ith) * cg(ik,isea)
#ifdef W3_MGP
      cxtot(ixy) = cxtot(ixy) - vgx/clats(isea)
      cytot(ixy) = cytot(ixy) - vgy
#endif
#ifdef W3_T1
      IF ( .NOT. flcur )                                        &
           WRITE (ndst,9021) isea, mapsf(isea,1), mapsf(isea,2), &
           vq(ixy), cxtot(ixy), cytot(ixy)
#endif
    END DO
    !
#ifdef W3_OMPH
    !$OMP END PARALLEL DO
#endif
    !
    IF ( flcur ) THEN
#ifdef W3_T
      WRITE (ndst,9022)
#endif
      !
#ifdef W3_OMPH
      !$OMP PARALLEL DO PRIVATE (ISEA,IXY)
#endif
      !
      DO isea=1, nsea
        ixy         = mapsf(isea,3)
        cxtot(ixy) = cxtot(ixy) + cx(isea)/clats(isea)
        cytot(ixy) = cytot(ixy) + cy(isea)
#ifdef W3_T1
        WRITE (ndst,9021) isea, mapsf(isea,1), mapsf(isea,2), &
             vq(ixy), cxtot(ixy), cytot(ixy)
#endif
      END DO
      !
#ifdef W3_OMPH
      !$OMP END PARALLEL DO
#endif
      !
    END IF
    !
#ifdef W3_OMPH
    !$OMP PARALLEL DO PRIVATE (ISEA,IX, IY, IXY, CP, CQ)
#endif
    !
    DO isea=1, nsea
      ix     = mapsf(isea,1)
      iy     = mapsf(isea,2)
      ixy    = mapsf(isea,3)
      cp = cxtot(ixy)*dpdx(iy,ix) + cytot(ixy)*dpdy(iy,ix)
      cq = cxtot(ixy)*dqdx(iy,ix) + cytot(ixy)*dqdy(iy,ix)
      vlcflx(ixy) = cp*dtrad
      vlcfly(ixy) = cq*dtrad
    END DO
    !
#ifdef W3_OMPH
    !$OMP END PARALLEL DO
#endif
    !
    ! 3.  Loop over sub-steps -------------------------------------------- *
    !
    DO itloc=1, ntloc
      !
      ! 3.a Average
      !
      aq     = vq
      vq     = 0.
      !
      ! 3.a.1 Central points
      !
      DO ip=1, ncent
        isea    = mapcxy(ip)
        ixy     = mapsf(isea,3)
        IF ( maptrn(ixy) ) THEN
          vq(ixy) = aq(ixy)
        ELSE
          rd1     = rdi1(isea)
          rd2     = min( 1. , rdi2(isea) * cg(ik,isea) )
          rd3     = rdi3(isea)
          rd4     = min( 1. , rdi4(isea) * cg(ik,isea) )
          vq(ixy          ) = vq(ixy          )                   &
               + aq(ixy) * (3.-rd2-rd4)/3.
          vq(ixy+iixy1(isea)) = vq(ixy+iixy1(isea))               &
               + aq(ixy) * rd2*rd1/6.
          vq(ixy+iixy2(isea)) = vq(ixy+iixy2(isea))               &
               + aq(ixy) * (1.-rd1)*rd2/6.
          vq(ixy+iixy3(isea)) = vq(ixy+iixy3(isea))               &
               + aq(ixy) * rd4*rd3/6.
          vq(ixy+iixy4(isea)) = vq(ixy+iixy4(isea))               &
               + aq(ixy) * (1.-rd3)*rd4/6.
          vq(ixy-iixy1(isea)) = vq(ixy-iixy1(isea))               &
               + aq(ixy) * rd2*rd1/6.
          vq(ixy-iixy2(isea)) = vq(ixy-iixy2(isea))               &
               + aq(ixy) * (1.-rd1)*rd2/6.
          vq(ixy-iixy3(isea)) = vq(ixy-iixy3(isea))               &
               + aq(ixy) * rd4*rd3/6.
          vq(ixy-iixy4(isea)) = vq(ixy-iixy4(isea))               &
               + aq(ixy) * (1.-rd3)*rd4/6.
        END IF
      END DO
      !
      ! 3.a.2 Near-coast points
      !
      DO ip=ncent+1, nsea
        isea    = mapcxy(ip)
        ix      = mapsf(isea,1)
        ixy     = mapsf(isea,3)
        IF ( mapsta(ixy) .LE. 0 ) cycle
        IF ( maptrn(ixy) ) THEN
          vq(ixy) = aq(ixy)
        ELSE
          rd1     = rdi1(isea)
          rd3     = rdi3(isea)
          rd2     = min( 1. , rdi2(isea) * cg(ik,isea) )
          rd4     = min( 1. , rdi4(isea) * cg(ik,isea) )
          vq(ixy          ) = vq(ixy          )                   &
               + aq(ixy) * (3.-rd2-rd4)/3.
          !
          ixc    = sign(ny,iixy1(isea))
          iyc    = iixy1(isea) - ixc
          IF ( mapstx(ixy+iixy1(isea)) .GE. 1 .AND.               &
               .NOT. ( mapstx(ixy+ixc).LE.0 .AND.                 &
               mapstx(ixy+iyc).LE.0 ) ) THEN
            vq(ixy+iixy1(isea)) = vq(ixy+iixy1(isea))      &
                 + aq(ixy) * rd2*rd1/6.
          ELSE
            vq(ixy          ) = vq(ixy          )          &
                 + aq(ixy) * rd2*rd1/6.
          END IF
          IF ( mapstx(ixy-iixy1(isea)) .GE. 1 .AND.               &
               .NOT. ( mapstx(ixy-ixc).LE.0 .AND.                 &
               mapstx(ixy-iyc).LE.0 ) ) THEN
            vq(ixy-iixy1(isea)) = vq(ixy-iixy1(isea))      &
                 + aq(ixy) * rd2*rd1/6.
          ELSE
            vq(ixy          ) = vq(ixy          )          &
                 + aq(ixy) * rd2*rd1/6.
          END IF
 
          IF ( mapstx(ixy+iixy2(isea)) .GE. 1 ) THEN
            vq(ixy+iixy2(isea)) = vq(ixy+iixy2(isea))      &
                 + aq(ixy) * (1.-rd1)*rd2/6.
          ELSE
            vq(ixy          ) = vq(ixy          )          &
                 + aq(ixy) * (1.-rd1)*rd2/6.
          END IF
          IF ( mapstx(ixy-iixy2(isea)) .GE. 1 ) THEN
            vq(ixy-iixy2(isea)) = vq(ixy-iixy2(isea))      &
                 + aq(ixy) * (1.-rd1)*rd2/6.
          ELSE
            vq(ixy          ) = vq(ixy          )          &
                 + aq(ixy) * (1.-rd1)*rd2/6.
          END IF
          !
          ixc    = sign(ny,iixy3(isea))
          iyc    = iixy3(isea) - ixc
          IF ( mapstx(ixy+iixy3(isea)) .GE. 1 .AND.               &
               .NOT. ( mapstx(ixy+ixc).LE.0 .AND.                 &
               mapstx(ixy+iyc).LE.0 ) ) THEN
            vq(ixy+iixy3(isea)) = vq(ixy+iixy3(isea))      &
                 + aq(ixy) * rd4*rd3/6.
          ELSE
            vq(ixy          ) = vq(ixy          )          &
                 + aq(ixy) * rd4*rd3/6.
          END IF
          IF ( mapstx(ixy-iixy3(isea)) .GE. 1 .AND.               &
               .NOT. ( mapstx(ixy-ixc).LE.0 .AND.                 &
               mapstx(ixy-iyc).LE.0 ) ) THEN
            vq(ixy-iixy3(isea)) = vq(ixy-iixy3(isea))      &
                 + aq(ixy) * rd4*rd3/6.
          ELSE
            vq(ixy          ) = vq(ixy          )          &
                 + aq(ixy) * rd4*rd3/6.
          END IF
          !
          IF ( mapstx(ixy+iixy4(isea)) .GE. 1 ) THEN
            vq(ixy+iixy4(isea)) = vq(ixy+iixy4(isea))      &
                 + aq(ixy) * (1.-rd3)*rd4/6.
          ELSE
            vq(ixy          ) = vq(ixy          )          &
                 + aq(ixy) * (1.-rd3)*rd4/6.
          END IF
          IF ( mapstx(ixy-iixy4(isea)) .GE. 1 ) THEN
            vq(ixy-iixy4(isea)) = vq(ixy-iixy4(isea))      &
                 + aq(ixy) * (1.-rd3)*rd4/6.
          ELSE
            vq(ixy          ) = vq(ixy          )          &
                 + aq(ixy) * (1.-rd3)*rd4/6.
          END IF
          !
        END IF
        !
      END DO
      !
      ! 3.a.3 Restore boundary data
      !
      DO ixy=1, nx*ny
        IF ( mapsta(ixy).EQ.2 ) vq(ixy) = aq(ixy)
      END DO
      !
      ! 3.a.4 Global closure (averaging only, propagation is closed in W3QCK3).
      !
      IF ( global ) THEN
        DO iy=1, ny
          vq(iy          ) = vq(iy          ) + vq(nx*ny+iy)
          vq((nx-1)*ny+iy) = vq((nx-1)*ny+iy) + vq(iy-ny)
        END DO
      END IF
      !
      ! 3.b Propagate fields
      !
      !     Transform VQ to straightened space
      !
#ifdef W3_OMPH
      !$OMP PARALLEL DO PRIVATE (ISEA, IX, IY, IXY)
#endif
      !
      DO isea=1, nsea
        ix     = mapsf(isea,1)
        iy     = mapsf(isea,2)
        ixy    = mapsf(isea,3)
        vq(ixy)= vq(ixy) * gsqrt(iy,ix)
      END DO
      !
#ifdef W3_OMPH
      !$OMP END PARALLEL DO
#endif
      !
      IF ( yfirst ) THEN
        !
#ifdef W3_UQ
        IF ( flcy ) CALL w3qck3                               &
             (nx, ny, nx, ny, vlcfly, atrny, vq,       &
             .false., 1, mapaxy, nact, mapy2, nmy0,   &
             nmy1, nmy2, ndse, ndst )
        IF ( flcx ) CALL w3qck3                               &
             (nx, ny, nx, ny, vlcflx, atrnx, vq,       &
             global, ny, mapaxy, nact, mapx2, nmx0,   &
             nmx1, nmx2, ndse, ndst )
#endif
        !
#ifdef W3_UNO
        IF ( flcy ) CALL w3uno2s                             &
             (nx, ny, nx, ny, vlcfly, atrny, vq,      &
             .false., 1, mapaxy, nact, mapy2, nmy0,  &
             nmy1, nmy2, ndse, ndst )
        IF ( flcx ) CALL w3uno2s                             &
             (nx, ny, nx, ny, vlcflx, atrnx, vq,      &
             global, ny, mapaxy, nact, mapx2, nmx0,  &
             nmx1, nmx2, ndse, ndst )
#endif
        !
      ELSE
        !
#ifdef W3_UQ
        IF ( flcx ) CALL w3qck3                               &
             (nx, ny, nx, ny, vlcflx, atrnx, vq,       &
             global, ny, mapaxy, nact, mapx2, nmx0,   &
             nmx1, nmx2, ndse, ndst )
        IF ( flcy ) CALL w3qck3                               &
             (nx, ny, nx, ny, vlcfly, atrny, vq,       &
             .false., 1, mapaxy, nact, mapy2, nmy0,   &
             nmy1, nmy2, ndse, ndst )
#endif
        !
#ifdef W3_UNO
        IF ( flcx ) CALL w3uno2s                             &
             (nx, ny, nx, ny, vlcflx, atrnx, vq,      &
             global, ny, mapaxy, nact, mapx2, nmx0,  &
             nmx1, nmx2, ndse, ndst )
        IF ( flcy ) CALL w3uno2s                             &
             (nx, ny, nx, ny, vlcfly, atrny, vq,      &
             .false., 1, mapaxy, nact, mapy2, nmy0,  &
             nmy1, nmy2, ndse, ndst )
#endif
        !
      END IF
      !
      !     Transform VQ back to normal space
      !
#ifdef W3_OMPH
      !$OMP PARALLEL DO PRIVATE (ISEA, IX, IY, IXY)
#endif
      !
      DO isea=1, nsea
        ix     = mapsf(isea,1)
        iy     = mapsf(isea,2)
        ixy    = mapsf(isea,3)
        vq(ixy)= vq(ixy) / gsqrt(iy,ix)
      END DO
      !
#ifdef W3_OMPH
      !$OMP END PARALLEL DO
#endif
      !
      ! 3.c Update boundaries
      !
#ifdef W3_T
      WRITE (ndst,9030) nsea
#endif
      !
#ifdef W3_T2
      DO isea=1, nsea
        ixy    = mapsf(isea,3)
        IF ( mapsta(ixy) .GT. 0 ) THEN
          WRITE(ndst,9031)isea,mapsf(isea,1),mapsf(isea,2),vq(ixy)
          vq(ixy) =  max( 0. , cg(ik,isea)/clats(isea)*vq(ixy) )
        END IF
      END DO
#endif
      !
      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
          ixy     = mapsf(isbpi(ibi),3)
          vq(ixy) = ( rd1*bbpi0(isp,ibi) + rd2*bbpin(isp,ibi) )   &
               / cg(ik,isbpi(ibi)) * clats(isbpi(ibi))
        END DO
      END IF
      !
      yfirst = .NOT. yfirst
    END DO
    !
    ! 4.  Store results in VQ in proper format --------------------------- *
    !
#ifdef W3_T
    WRITE (ndst,9040) nsea
#endif
    !
#ifdef W3_OMPH
    !$OMP PARALLEL DO PRIVATE (ISEA, IXY)
#endif
    !
    DO isea=1, nsea
      ixy    = mapsf(isea,3)
      IF ( mapsta(ixy) .GT. 0 ) THEN
#ifdef W3_T3
        WRITE (ndst,9041) isea, mapsf(isea,1), mapsf(isea,2), vq(ixy)
#endif
        vq(ixy) =  max( 0. , cg(ik,isea)/clats(isea)*vq(ixy) )
      END IF
    END DO
    !
#ifdef W3_OMPH
    !$OMP END PARALLEL DO
#endif
    !
    RETURN
    !
    ! Formats
    !
#ifdef W3_T
9000 FORMAT (' TEST W3XYP3 : YFIRST  :',l2)
9001 FORMAT (' TEST W3XYP3 : ISP, ITH, IK, COS-SIN :',i8,2i4,2f7.3)
9010 FORMAT (' TEST W3XYP3 : INITIALIZE ARRAYS')
9020 FORMAT (' TEST W3XYP3 : CALCULATING VCFL0X/Y (NSEA=',i6,')')
9022 FORMAT (' TEST W3XYP3 : CALCULATING VCFLUX/Y')
9030 FORMAT (' TEST W3XYP3 : FIELD BEFORE BPI. (NSEA=',i6,')')
9040 FORMAT (' TEST W3XYP3 : FIELD AFTER PROP. (NSEA=',i6,')')
#endif
#ifdef W3_T0
9011 FORMAT (' TEST W3XYP3 : PREPARE AVERAGING')
9012 FORMAT ('         ',4i4,2f7.3)
9013 FORMAT ('         ',8x,2i4,2f7.3)
#endif
    !
#ifdef W3_T1
9021 FORMAT (1x,i6,2i5,e12.4,2f7.3)
#endif
    !
#ifdef W3_T2
9031 FORMAT (1x,i6,2i5,e12.4)
#endif
    !
#ifdef W3_T3
9041 FORMAT (1x,i6,2i5,e12.4)
#endif
    !/
    !/ End of W3XYP3 ----------------------------------------------------- /
    !/
  END SUBROUTINE w3xyp3
  !/ ------------------------------------------------------------------- /
  SUBROUTINE w3ktp3 ( ISEA, FACTH, FACK, CTHG0, CG, WN, DW,       &
       DDDX, DDDY, CX, CY, DCXDX, DCXDY,           &
       DCYDX, DCYDY, DCDX, DCDY, VA, CFLTHMAX, CFLKMAX )
    !/
    !/    *** THIS ROUTINE SHOULD BE IDENTICAL TO W3KTP2 ***
    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III           NOAA/NCEP |
    !/                  |           H. L. Tolman            |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         01-Jul-2013 |
    !/                  +-----------------------------------+
    !/
    !/    14-Feb-2000 : Origination.                        ( version 2.08 )
    !/    17-Dec-2004 : Multiple grid version.              ( version 3.06 )
    !/    06-Mar-2011 : Output of maximum CFL  (F. Ardhuin) ( version 3.14 )
    !/    24-Aug-2011 : Limiter on k advection (F. Ardhuin) ( version 4.04 )
    !/    25-Aug-2011 : DEPTH  = MAX ( DMIN, DW(ISEA) )     ( version 4.04 )
    !/    26-Dec-2012 : More initializations.               ( version 4.11 )
    !/    01-Jul-2013 : Adding UQ and UNO switches to chose between third
    !/                  and second order schemes.           ( version 4.12 )
    !/
    !  1. Purpose :
    !
    !     Propagation in spectral space.
    !
    !  2. Method :
    !
    !     Third order QUICKEST scheme with ULTIMATE limiter.
    !
    !     As with the spatial propagation, the two spaces are considered
    !     independently, but the propagation is performed in a 2-D space.
    !     Compared to the propagation in physical space, the directions
    !     rerpesent a closed space and are therefore comparable to the
    !     longitudinal or 'X' propagation. The wavenumber space has to be
    !     extended to allow for boundary treatment. Using a simple first
    !     order boundary treatment at both sided, two points need to
    !     be added. This implies that the spectrum needs to be extended,
    !     shifted and rotated, as is performed using MAPTH2 as set
    !     in W3MAP3.
    !
    !  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.
    !       DW      R.A.   I   Depth.
    !       DDDx    Real   I   Depth gradients.
    !       CX/Y    Real   I   Current components.
    !       DCxDx   Real   I   Current gradients.
    !       DCDX-Y  Real   I   Phase speed gradients.
    !       VA      R.A.  I/O  Spectrum.
    !     ----------------------------------------------------------------
    !
    !     Local variables.
    !     ----------------------------------------------------------------
    !       DSDD    R.A.  Partial derivative of sigma for depth.
    !       FDD, FDU, FDG, 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 :
    !
    !       W3QCK1   Actual propagation routine.
    !       W3QCK2   Actual propagation routine.
    !       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 DSDD
    !         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, nk2, nth, nspec, sig, dsip, ecos, esin, &
         ec2, esc, es2, fachfa, mapwn, flcth, flck,  &
         ctmax, dmin
    USE w3adatmd, ONLY: mapth2, mapwn2, itime
    USE w3idatmd, ONLY: flcur
    USE w3odatmd, ONLY: ndse, ndst
#ifdef W3_S
    USE w3servmd, ONLY: strace
#endif
#ifdef W3_UQ
    USE w3uqckmd
#endif
#ifdef W3_UNO
    USE w3uno2md
#endif
    !/
    IMPLICIT NONE
    !/
    !/ ------------------------------------------------------------------- /
    !/ Parameter list
    !/
    INTEGER, INTENT(IN)     :: ISEA
    REAL, INTENT(IN)        :: FACTH, FACK, CTHG0, CG(0:NK+1),      &
         WN(0:NK+1), DW, DDDX, DDDY,       &
         CX, CY, DCXDX, DCXDY, DCYDX, DCYDY
    REAL, INTENT(IN)        :: DCDX(0:NK+1), DCDY(0:NK+1)
    REAL, INTENT(INOUT)     :: VA(NSPEC)
    REAL, INTENT(OUT)       :: CFLTHMAX, CFLKMAX
    !/
    !/ ------------------------------------------------------------------- /
    !/ Local parameters
    !/
    INTEGER                 :: ITH, IK, ISP
#ifdef W3_S
    INTEGER, SAVE           :: IENT = 0
#endif
    REAL                    :: FDDMAX, FDG, FKD, FKD0, DCYX,        &
         DCXXYY, DCXY, DCXX, DCXYYX, DCYY,    &
         VELNOFILT, VELFAC, DEPTH
    REAL                    :: DSDD(0:NK+1), FRK(NK), FRG(NK),      &
         FKC(NTH), VQ(-NK-1:NK2*(NTH+2)),     &
         DB(NK2,NTH+1), DM(NK2,0:NTH+1),      &
         VCFLT(NK2*(NTH+1)), CFLK(NK2,NTH)
    !/
    !/ ------------------------------------------------------------------- /
    !/
#ifdef W3_S
    CALL strace (ient, 'W3KTP3')
#endif
    !
    ! 1.  Preparations --------------------------------------------------- *
    ! 1.a Initialize arrays
    !
    depth  = max( dmin, dw )
    vq       = 0.
    IF ( flcth ) vcflt    = 0.
    IF ( flck  ) cflk     = 0.
    cflthmax = 0.
    cflkmax  = 0.
    !
#ifdef W3_T
    WRITE (ndst,9000) flcth, flck, facth, fack, ctmax
    WRITE (ndst,9010) isea, depth, cx, cy, dddx, dddy,           &
         dcxdx, dcxdy, dcydx, dcydy
#endif
    !
    ! 2.  Preparation for point ------------------------------------------ *
    ! 2.a Array with partial derivative of sigma versus depth
    !
    DO ik=0, nk+1
      IF ( depth*wn(ik) .LT. 5. ) THEN
        dsdd(ik) = max( 0. ,                                     &
             cg(ik)*wn(ik)-0.5*sig(ik) ) / depth
      ELSE
        dsdd(ik) = 0.
      END IF
    END DO
    !
#ifdef W3_T
    WRITE (ndst,9020)
    DO ik=1, nk+1
      WRITE (ndst,9021) ik, tpi/sig(ik), tpi/wn(ik),             &
           cg(ik), dsdd(ik)
    END DO
#endif
    !
    ! 2.b Extract spectrum
    !
    DO isp=1, nspec
      vq(mapth2(isp)) = va(isp)
    END DO
    !
    ! 3.  Refraction velocities ------------------------------------------ *
    !
    IF ( flcth ) THEN
      !
      ! 3.a Set slope filter for depth refraction
      !
      ! N.B.:  FACTH = DTG / DTH / REAL(NTLOC)  (value set in w3wavemd)
      !        namely, FACTH*VC=1 corresponds to CFL=1
      !
      fddmax = 0.
      fdg    = facth * cthg0
      !
      DO ith=1, nth/2
        fddmax = max(fddmax,abs(esin(ith)*dddx-ecos(ith)*dddy))
      END DO
      !
      DO ik=1, nk
        frk(ik) = facth * dsdd(ik) / wn(ik)
        !
        ! Removes the filtering that was done at that stage (F. Ardhuin 2011/03/06)
        !
        !            FRK(IK) = FRK(IK) / MAX ( 1. , FRK(IK)*FDDMAX/CTMAX )
        frg(ik) = fdg * cg(ik)
      END DO
      !
      ! 3.b Current refraction
      !
      IF ( flcur ) THEN
        !
        dcyx   = facth *   dcydx
        dcxxyy = facth * ( dcxdx - dcydy )
        dcxy   = facth *   dcxdy
        !
        DO isp=1, nspec
          vcflt(mapth2(isp)) = es2(isp)*dcyx  +     &
               esc(isp)*dcxxyy - ec2(isp)*dcxy
        END DO
        !
      ELSE
        vcflt(:)=0.
      END IF
      !
      ! 3.c Depth refraction and great-circle propagation
      !
#ifdef W3_REFRX
      ! 3.d @C/@x refraction and great-circle propagation
      frk    = 0.
      fddmax = 0.
      !
      DO isp=1, nspec
        fddmax = max( fddmax , abs(                      &
             esin(isp)*dcdx(mapwn(isp)) - ecos(isp)*dcdy(mapwn(isp)) ) )
      END DO
      DO ik=1, nk
        frk(ik) = facth * cg(ik) * wn(ik) / sig(ik)
      END DO
#endif
      !
      DO isp=1, nspec
        velnofilt = vcflt(mapth2(isp))       &
             + frg(mapwn(isp)) * ecos(isp)              &
             + frk(mapwn(isp)) * ( esin(isp)*dddx - ecos(isp)*dddy )
        !
#ifdef W3_REFRX
        ! 3.d @C/@x refraction and great-circle propagation
        velnofilt = vcflt(mapth2(isp))       &
             + frg(mapwn(isp)) * ecos(isp)            &
             + frk(mapwn(isp)) * ( esin(isp)*dcdx(mapwn(isp)) &
             - ecos(isp)*dcdy(mapwn(isp)) )
#endif
        cflthmax = max(cflthmax, abs(velnofilt))
        !
        ! Puts filtering on total velocity (including currents and great circle effects)
        ! the filtering limits VCFLT to be less than CTMAX
        ! this modification was proposed by F. Ardhuin 2011/03/06
        !
        vcflt(mapth2(isp))=sign(min(abs(velnofilt),ctmax),velnofilt)
      END DO
    END IF
    !
    ! 4.  Wavenumber shift velocities ------------------------------------ *
    ! N.B.:  FACK = DTG / REAL(NTLOC)  (value set in w3wavemd)
    !        namely, FACK*VC/DK=1 corresponds to CFL=1
    !
    IF ( flck ) THEN
      !
      ! 4.a Directionally dependent part
      !
      dcxx   =  -   dcxdx
      dcxyyx =  - ( dcxdy + dcydx )
      dcyy   =  -   dcydy
      fkd    =    ( cx*dddx + cy*dddy )
      !
      DO ith=1, nth
        fkc(ith) = ec2(ith)*dcxx +                                &
             esc(ith)*dcxyyx + es2(ith)*dcyy
      END DO
      !
      ! 4.b Band widths
      !
      DO ik=0, nk
        db(ik+1,1) = dsip(ik) / cg(ik)
        dm(ik+1,1) = wn(ik+1) - wn(ik)
      END DO
      db(nk+2,1) = dsip(nk+1) / cg(nk+1)
      dm(nk+2,1) = 0.
      !
      DO ith=2, nth
        DO ik=1, nk+2
          db(ik,ith) = db(ik,1)
          dm(ik,ith) = dm(ik,1)
        END DO
      END DO
      !
      ! 4.c Velocities
      !
      DO ik=0, nk+1
        fkd0   = fkd / cg(ik) * dsdd(ik)
        velfac =  fack/db(ik+1,1)
        DO ith=1, nth
          !
          ! Puts filtering on velocity (needs the band widths)
          !
          velnofilt = ( fkd0 + wn(ik)*fkc(ith) ) * velfac   ! this is velocity * DT / DK
          cflkmax = max(cflkmax, abs(velnofilt))
          cflk(ik+1,ith) = sign(min(abs(velnofilt),ctmax),velnofilt)/velfac
          !CFLK(IK+1,ITH) = FKD0 + WN(IK)*FKC(ITH)          ! this was without the limiter ...
        END DO
      END DO
      !
    END IF
    !
    ! 5.  Propagate ------------------------------------------------------ *
    !
    IF ( mod(itime,2) .EQ. 0 ) THEN
      IF ( flck ) THEN
        DO ith=1, nth
          vq(nk+2+(ith-1)*nk2) = fachfa * vq(nk+1+(ith-1)*nk2)
        END DO
        !
#ifdef W3_UQ
        CALL w3qck2 ( nth, nk2, nth, nk2, cflk, fack, db, dm,   &
             vq, .false., 1, mapth2, nspec,            &
             mapwn2, nspec-nth, nspec, nspec+nth,      &
             ndse, ndst )
#endif
        !
#ifdef W3_UNO
        CALL w3uno2 ( nth, nk2, nth, nk2, cflk, fack, db, dm,  &
             vq, .false., 1, mapth2, nspec,           &
             mapwn2, nspec-nth, nspec, nspec+nth,     &
             ndse, ndst )
#endif
        !
      END IF
      IF ( flcth ) THEN
        !
#ifdef W3_UQ
        CALL w3qck1 ( nth, nk2, nth, nk2, vcflt, vq, .true.,    &
             nk2, mapth2, nspec, mapth2, nspec, nspec, &
             nspec, ndse, ndst )
#endif
        !
#ifdef W3_UNO
        CALL w3uno2r( nth, nk2, nth, nk2, vcflt, vq, .true.,   &
             nk2, mapth2, nspec, mapth2, nspec, nspec,&
             nspec, ndse, ndst )
#endif
        !
      END IF
    ELSE
      IF ( flcth ) THEN
        !
#ifdef W3_UQ
        CALL w3qck1 ( nth, nk2, nth, nk2, vcflt, vq, .true.,    &
             nk2, mapth2, nspec, mapth2, nspec, nspec, &
             nspec, ndse, ndst )
#endif
        !
#ifdef W3_UNO
        CALL w3uno2r( nth, nk2, nth, nk2, vcflt, vq, .true.,   &
             nk2, mapth2, nspec, mapth2, nspec, nspec,&
             nspec, ndse, ndst )
#endif
        !
      END IF
      IF ( flck ) THEN
        DO ith=1, nth
          vq(nk+2+(ith-1)*nk2) = fachfa * vq(nk+1+(ith-1)*nk2)
        END DO
        !
#ifdef W3_UQ
        CALL w3qck2 ( nth, nk2, nth, nk2, cflk, fack, db, dm,   &
             vq, .false., 1, mapth2, nspec,            &
             mapwn2, nspec-nth, nspec, nspec+nth,      &
             ndse, ndst )
#endif
        !
#ifdef W3_UNO
        CALL w3uno2 ( nth, nk2, nth, nk2, cflk, fack, db, dm,  &
             vq, .false., 1, mapth2, nspec,           &
             mapwn2, nspec-nth, nspec, nspec+nth,     &
             ndse, ndst )
#endif
        !
      END IF
    END IF
    !
    ! 6.  Store reults --------------------------------------------------- *
    !
    DO isp=1, nspec
      va(isp) = vq(mapth2(isp))
    END DO
    !
    RETURN
    !
    ! Formats
    !
#ifdef W3_T
9000 FORMAT ( ' TEST W3KTP3 : FLCTH-K, FACTH-K, CTMAX  :',        &
         2l2,2e10.3,f7.3)
9010 FORMAT ( ' TEST W3KTP3 : LOCAL DATA :',i7,f7.1,2f6.2,1x,6e10.2)
9020 FORMAT ( ' TEST W3KTP3 : IK, T, L, CG, DSDD : ')
9021 FORMAT ( '              ',i3,f7.2,f7.1,f7.2,e11.3)
#endif
    !
#ifdef W3_T0
9040 FORMAT (/' TEST W3KTP3 : NORMALIZED ',a/)
9041 FORMAT (1x,60(1x,i2))
9042 FORMAT (1x,60i3)
#endif
    !/
    !/ End of W3KTP3 ----------------------------------------------------- /
    !/
  END SUBROUTINE w3ktp3
  !/ ------------------------------------------------------------------- /
  SUBROUTINE w3cflxy ( ISEA, DTG, MAPSTA, MAPFS, CFLXYMAX, VGX, VGY )
    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III           NOAA/NCEP |
    !/                  |           F. Ardhuin            |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         31-Oct-2010 |
    !/                  +-----------------------------------+
    !/
    !/    07-Mar-2011 : Origination.                        ( version 3.14 )
    !/
    !  1. Purpose :
    !
    !     Computes the maximum CFL number for spatial advection. Used for diagnostic
    !     purposes. (Could be used to define a local time step ...)
    !
    !  2. Method :
    !
    !
    !  3. Parameters :
    !
    !     Parameter list
    !     ----------------------------------------------------------------
    !       ISEA    Int.   I   Index of grid point.
    !       DTG     Real   I   Total time step.
    !       MAPSTA  I.A.   I   Grid point status map.
    !       MAPFS   I.A.   I   Storage map.
    !       CFLXYMAX Real  O   Maximum CFL number for XY propagation.
    !       VGX/Y   Real   I   Speed of grid.
    !     ----------------------------------------------------------------
    !
    !     Local variables.
    !     ----------------------------------------------------------------
    !       NTLOC   Int   Number of local time steps.
    !       DTLOC   Real  Local propagation time step.
    !       VCFL0X  R.A.  Local courant numbers for absolute group vel.
    !                     using local X-grid step.
    !       VCFL0Y  R.A.  Id. in Y.
    !     ----------------------------------------------------------------
    !
    !  4. Subroutines used :
    !
    !       STRACE   Service routine.
    !
    !  5. Called by :
    !
    !       W3WAVE   Wave model routine.
    !
    !  6. Error messages :
    !
    !       None.
    !
    !  7. Remarks :
    !
    !     - Curvilinear grid implementation. Variables FACX, FACY, CCOS, CSIN,
    !       CCURX, CCURY are not needed and have been removed.  FACX is accounted
    !       for as approriate in this subroutine.  FACX is also accounted for in
    !       the case of .NOT.FLCX.  Since FACX is removed, there is now a check for
    !       .NOT.FLCX in this subroutine.  In CFL calcs dx and dy are omitted,
    !       since dx=dy=1 in index space.  Curvilinear grid derivatives
    !       (DPDY, DQDX, etc.) and metric (GSQRT) are brought in via W3GDATMD.
    !
    !  8. Structure :
    !
    !     ---------------------------------------------
    !     ---------------------------------------------
    !
    !  9. Switches :
    !
    !       !/S     Enable subroutine tracing.
    !
    !       !/MGP   Moving grid corrections.
    !       !/MGG   Moving grid corrections.
    !
    !       !/T     Enable general test output.
    !
    ! 10. Source code :
    !
    !/ ------------------------------------------------------------------- /
    USE constants
    !
    USE w3timemd, ONLY: dsec21
    !
    USE w3gdatmd, ONLY: nx, ny, nsea, mapsf, dtcfl, clats,      &
         flcx, flcy, nk, nth, dth, xfr,          &
         ecos, esin, sig, wdcg, wdth, pfmove,    &
         flagll, dpdx, dpdy, dqdx, dqdy, gsqrt
    USE w3wdatmd, ONLY: time
    USE w3adatmd, ONLY: nmx0, nmx1, nmx2, nmy0, nmy1, nmy2, nact,   &
         ncent, mapx2, mapy2, mapaxy, mapcxy,        &
         maptrn, cg, 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)     :: ISEA, MAPSTA(NY*NX), MAPFS(NY*NX)
    REAL, INTENT(IN)        :: DTG, VGX, VGY
    REAL, INTENT(INOUT)     :: CFLXYMAX
    !/
    !/ ------------------------------------------------------------------- /
    !/ Local parameters
    !/
    INTEGER                 :: ITH, IK, IXY, IP
    INTEGER                 :: IX, IY, IXC, IYC, IBI
#ifdef W3_S
    INTEGER, SAVE           :: IENT = 0
#endif
    REAL                    :: CG0, CGA, CGN, CGX, CGY, CXC, CYC,   &
         CXMIN, CXMAX, CYMIN, CYMAX
    REAL                    :: CGC, FGSE = 1.
    REAL                    :: FTH, FTHX, FTHY, FCG, FCGX, FCGY
    REAL                    :: CP, CQ
    !/
    !/ Automatic work arrays
    !/
    REAL                    :: VLCFLX, VLCFLY
    REAL                    :: CXTOT, CYTOT
    !/
    !/ ------------------------------------------------------------------- /
    !/
#ifdef W3_S
    CALL strace (ient, 'W3XYCFL')
#endif
    !
    ! 1.  Preparations --------------------------------------------------- *
    ! 1.a Set constants
    !
    !
    cflxymax=0.
    ix  = mapsf(isea,1)
    iy  = mapsf(isea,2)
    ixy = mapsf(isea,3)
    DO ik=1,nk
      DO ith=1,nth
        cxtot = ecos(ith) * cg(ik,isea) / clats(isea)
        cytot = esin(ith) * cg(ik,isea)
#ifdef W3_MGP
        cxtot = cxtot - vgx/clats(isea)
        cytot = cytot - vgy
#endif
 
 
        IF ( flcur ) THEN
          cxtot = cxtot + cx(isea)/clats(isea)
          cytot = cytot + cy(isea)
        END IF
 
        cp = cxtot*dpdx(iy,ix) + cytot*dpdy(iy,ix)
        cq = cxtot*dqdx(iy,ix) + cytot*dqdy(iy,ix)
        vlcflx = cp*dtg
        vlcfly = cq*dtg
        cflxymax = max(vlcflx,vlcfly,cflxymax)
      END DO
    END DO
 
    RETURN
    !/
    !/ End of W3XYCFL ----------------------------------------------------- /
    !/
  END SUBROUTINE w3cflxy
 
  !/
  !/ End of module W3PRO3MD -------------------------------------------- /
  !/
END MODULE w3pro3md