!> @file
!> @brief Bundles routines for third order porpagation scheme in single
!>  module.
!>
!> @author H. L. Tolman
!> @date   29-May-2014
!>

#include "w3macros.h"
!/ ------------------------------------------------------------------- /

!>
!> @brief Bundles routines for third order porpagation scheme in single
!>  module.
!>
!> @author H. L. Tolman
!> @date   29-May-2014
!>
!> @copyright Copyright 2009-2022 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.
!>
MODULE W3PRO2MD
  !/
  !/                  +-----------------------------------+
  !/                  | WAVEWATCH III           NOAA/NCEP |
  !/                  |           H. L. Tolman            |
  !/                  |                        FORTRAN 90 |
  !/                  | Last update :         29-May-2014 |
  !/                  +-----------------------------------+
  !/
  !/    04-Feb-2000 : Origination.                        ( version 2.00 )
  !/    24-Jan-2001 : Flat grid version                   ( version 2.06 )
  !/    08-Feb-2001 : UQ routines moved to own module     ( version 2.08 )
  !/    09-Feb-2001 : Clean up of parameter lists         ( version 2.08 )
  !/    14-Feb-2001 : Unit numbers in UQ routines         ( version 2.08 )
  !/    13-Nov-2001 : Sub-grid obstacles added.           ( version 2.14 )
  !/    26-Dec-2002 : Moving grid option.                 ( version 3.02 )
  !/    20-Dec-2004 : Multiple grid version.              ( version 3.06 )
  !/    07-Sep-2005 : Improved 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 )
  !/    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)
  !/    06-Dec-2010 : Change from GLOBAL (logical) to ICLOSE (integer) to
  !/                  specify index closure for a grid.   ( version 3.14 )
  !/                  (T. J. Campbell, NRL)
  !/    23-Dec-2010 : Fix HPFAC and HQFAC by including the COS(YGRD)
  !/                  factor with DXDP and DXDQ terms.    ( version 3.14 )
  !/                  (T. J. Campbell, NRL)
  !/    01-Jul-2013 : Adding UQ and UNO switches to chose between third
  !/                  and second order schemes.           ( version 4.12 )
  !/    29-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 porpagation scheme in single
  !     module.
  !
  !  2. Variables and types :
  !
  !      Name      Type  Scope    Description
  !     ----------------------------------------------------------------
  !      TRNMIN    R.P.  Private   Minimum transparancy for local
  !     ----------------------------------------------------------------
  !
  !  3. Subroutines and functions :
  !
  !      Name      Type  Scope    Description
  !     ----------------------------------------------------------------
  !      W3MAP2    Subr. Public   Set up auxiliary maps.
  !      W3XYP2    Subr. Public   Third order spatial propagation.
  !      W3KTP2    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 :
  !
  !  6. Switches :
  !
  !       !/UQ    3rd order UQ propagation scheme.
  !       !/UNO   2nd order UNO propagation scheme.
  !
  !       !/MGP   Correct for motion of grid.
  !
  !       !/OMPH  Hybrid OpenMP directives.
  !
  !       !/TDYN, !/DSS0, !/XW0, !/XW1
  !               Diffusion options in W3XYP2
  !
  !       !/S     Enable subroutine tracing.
  !       !/Tn    Enable test output.
  !
  !  7. Source code :
  !
  !/ ------------------------------------------------------------------- /
  !/
  !/ Public variables
  !/
  PUBLIC
  !/
  !/ Private data
  !/
  REAL, PRIVATE, PARAMETER:: TRNMIN = 0.95
  !/
CONTAINS
  !/ ------------------------------------------------------------------- /
  !>
  !> @brief Generate 'map' arrays for the ULTIMATE QUICKEST scheme.
  !>
  !> @author H. L. Tolman
  !> @date   09-Nov-2005
  !>
  SUBROUTINE W3MAP2
    !/
    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III           NOAA/NCEP |
    !/                  |           H. L. Tolman            |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         09-Nov-2005 |
    !/                  +-----------------------------------+
    !/
    !/    19-Dec-1996 : Final FORTRAN 77                    ( version 1.18 )
    !/    15-Dec-1999 : Upgrade to FORTRAN 90               ( version 2.00 )
    !/    09-Feb-2001 : Clean up of parameter lists         ( version 2.08 )
    !/    20-Dec-2004 : Multiple grid version.              ( version 3.06 )
    !/    09-Nov-2005 : Removing soft boundary option.      ( version 3.08 )
    !/
    !  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.
    !     MAPXY contains similar information for the cross term in the
    !     diffusion correction (counter NMXY only).
    !
    !  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 and MAPXY
    !      4.   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
    USE W3ADATMD, ONLY: NMX0, NMX1, NMX2, NMY0, NMY1, NMY2, NACT,   &
         NMXY, MAPX2, MAPY2, MAPAXY, MAPXY,          &
         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,    &
         IK, ITH, ISP, ISP0, ISP2
#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, 'W3MAP2')
#endif
    !
    ! 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 and MAPXY ------------------------------------------- *
    !
    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
    !
    NMXY   = 0
    DO IX=1, NX
      IXY0   = (IX-1)*NY
      IX2    = IX+1
      IF ( IX .EQ. NX ) IX2 =  1
      DO IY=2, NY-1
        IF ( MAPSTA( IY ,IX ).GE.1 .AND.                            &
             MAPSTA( IY ,IX2).GE.1 .AND.                            &
             MAPSTA(IY+1,IX ).GE.1 .AND.                            &
             MAPSTA(IY+1,IX2).GE.1 ) THEN
          NMXY   = NMXY + 1
          MAPXY(NMXY) = IXY0 + IY
        END IF
      END DO
    END DO
    !
    ! 4.  Maps for intra-spectral propagation ---------------------------- *
    !
    IF ( MAPTH2(1) .NE. 0 ) RETURN
    !
#ifdef W3_T
    MAPTST = 0
#endif
    !
    ! 4.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
    !
    ! 4.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 W3MAP2 : TEST MAP FOR PROPAGATION'/          &
         '      MAP     : ',A/                               &
         '      CENTRAL : ',I6/                              &
         '      ABOVE   : ',I6/                              &
         '      BELOW   : ',I6/                              &
         '      TOTAL   : ',I6/)
9001 FORMAT (1X,130I1)
9010 FORMAT (' TEST W3MAP2 : COMPOSITE MAPS TH2, WN2 AND BTK')
9011 FORMAT (2X,60I2)
#endif
    !/
    !/ End of W3MAP2 ----------------------------------------------------- /
    !/
  END SUBROUTINE W3MAP2
  !/ ------------------------------------------------------------------- /
  !>
  !> @brief Propagation in physical space for a given spectral component.
  !>
  !> @param[in]    ISP     Number of spectral bin (IK-1)*NTH+ITH.
  !> @param[in]    DTG     Total time step.
  !> @param[in]    MAPSTA  Grid point status map.
  !> @param[in]    MAPFS   Storage map.
  !> @param[inout] VQ      Field to propagate.
  !> @param[in]    VGX
  !> @param[in]    VGY
  !>
  !> @author H. L. Tolman
  !> @date   29-May-2014
  !>
  SUBROUTINE W3XYP2 ( ISP, DTG, MAPSTA, MAPFS, VQ, VGX, VGY )
    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III           NOAA/NCEP |
    !/                  |           H. L. Tolman            |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         29-May-2014 |
    !/                  +-----------------------------------+
    !/
    !/    07-Jul-1998 : Final FORTRAN 77                    ( version 1.18 )
    !/    15-Dec-1999 : Upgrade to FORTRAN 90               ( version 2.00 )
    !/    24-Jan-2001 : Flat grid version                   ( version 2.06 )
    !/    09-Feb-2001 : Clean up of parameter lists         ( version 2.08 )
    !/    14-Feb-2001 : Unit numbers in UQ routines         ( version 2.08 )
    !/    13-Nov-2001 : Sub-grid obstructions.              ( version 2.14 )
    !/    26-Dec-2002 : Moving grid option,                 ( version 3.02 )
    !/    20-Dec-2004 : Multiple grid version.              ( version 3.06 )
    !/    07-Sep-2005 : Improved 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)
    !/    06-Dec-2010 : Change from GLOBAL (logical) to ICLOSE (integer) to
    !/                  specify index closure for a grid.   ( version 3.14 )
    !/                  (T. J. Campbell, NRL)
    !/    23-Dec-2010 : Fix HPFAC and HQFAC by including the COS(YGRD)
    !/                  factor with DXDP and DXDQ terms.    ( version 3.14 )
    !/                  (T. J. Campbell, NRL)
    !/    01-Jul-2013 : Adding UQ and UNO switches to chose between third
    !/                  and second order schemes.           ( version 4.12 )
    !/    29-May-2014 : Adding OMPH switch.                 ( version 5.02 )
    !/
    !  1. Purpose :
    !
    !     Propagation in physical space for a given spectral component.
    !
    !  2. Method :
    !
    !     Third-order ULTIMATE QUICKEST scheme and diffusion correction
    !     for linear dispersion (see manual).
    !     Curvilinear grid implementation: Fluxes are computed in index space
    !       and then transformed back into physical space.  The diffusion term
    !       is handled in physical space.
    !
    !  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.
    !     ----------------------------------------------------------------
    !
    !     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.
    !       VLCFLX  R.A.  Local courant numbers in index space (norm. velocities)
    !       VLCFLY  R.A.
    !       CXTOT   R.A.  Propagation velocities in physical space.
    !       CYTOT   R.A.
    !       CELLP   Real  Cell Reynolds/Peclet number used to calculate
    !                     diffusion coefficient for growing spectral
    !                     components.
    !       DFRR    Real  Relative frequency increment.
    !     ----------------------------------------------------------------
    !
    !  4. Subroutines used :
    !
    !     See module documentation.
    !
    !  5. Called by :
    !
    !      Name      Type  Module   Description
    !     ----------------------------------------------------------------
    !      W3WAVE    Subr. W3WAVEMD 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 distortion 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.
    !     - If TDYN is set to zero, ALL diffusion is skipped. Set TDYN
    !       to a small positive number to have growth diffusion only.
    !     - 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.
    !     - Factors VFDIFX_FAC VFDIFY_FAC VFDIFC_FAC are introduced so that results
    !       match for test case tp2.3.  Use of these factors is optional and removal
    !       can significantly reduce size/cost of code. These variants are marked as
    !       CURV1 or CURV2. NCEP will make final decision re: which version to adopt.
    !       CURV1 is the shorter version and results do not match the original code
    !       for all test cases. CURV2 is the longer version and results do match the
    !       original. DETAILED EXPLANATION: Discrepancies occur at the boundaries.
    !       This is because, at the boundaries, the pre-curvilinear version zeros out
    !       some terms in the diffusion calculation. Since they are zeroed out,
    !       they aren't there to *cancel* certain other terms: these "other terms"
    !       affect the result, so they have to be retained in the long vesion (CURV2)
    !       to get an exact match. In the short version, the "canceling out" is
    !       performed prior to coding the scheme, so both the canceled and canceling
    !       terms are always omitted.
    !
    !  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.
    !
    !       !/OMPH  Hybrid OpenMP directives.
    !
    !       !/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, ICLOSE, FLCX, FLCY,    &
         ICLOSE_NONE, ICLOSE_SMPL, ICLOSE_TRPL,      &
         DTCFL, CLATS, DTME, CLATMN, FLAGLL,         &
         HPFAC, HQFAC, DPDX, DPDY, DQDX, DQDY, GSQRT
    USE W3WDATMD, ONLY: TIME
    USE W3ADATMD, ONLY: CG, WN, U10, CX, CY, ATRNX, ATRNY, ITIME,   &
         NMX0, NMX1, NMX2, NMY0, NMY1, NMY2, NACT,   &
         NMXY, MAPX2, MAPY2, MAPAXY, MAPXY
    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,    &
         IX,IY, IY0, IP, IBI
    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                    :: DTLOC, DTRAD,               &
         DFRR, CELLP,  CGD, DSSD,    &
         DNND, DCELL, XWIND, TFAC, DSS, DNN
    REAL                    :: RD1, RD2
    REAL                    :: RFAC, DFAC, DVQ, QXX, QXY, QYY
    REAL                    :: CP, CQ
    LOGICAL                 :: YFIRST
    LOGICAL                 :: GLOBAL
    !/
    !/ Automatic work arrays
    !/
    REAL                    :: VLCFLX((NX+1)*NY), VLCFLY((NX+1)*NY),&
         VFDIFX(1-NY:NX*NY), VFDIFY(NX*NY),   &
         VFDIFC(1-NY:NX*NY), VDXX((NX+1)*NY), &
         VDYY(NX*NY), VDXY((NX+1)*NY)

    REAL                    :: CXTOT((NX+1)*NY), CYTOT(NX*NY)
    REAL                    :: VQ_OLD(1-NY:NY*(NX+2))
    !CURV2: BEGIN -----------------------------------------------------------------
    REAL                    :: VFDIFX_FAC(1-NY:NX*NY),              &
         VFDIFY_FAC(1-NY:NX*NY),              &
         VFDIFC_FAC(1-NY:NX*NY)
    !CURV2: END -------------------------------------------------------------------
    !/
    !/ ------------------------------------------------------------------- /
    !/
#ifdef W3_S
    CALL STRACE (IENT, 'W3XYP2')
#endif
    !
    !     IF ( MAXVAL(VQ) .EQ. 0. ) THEN
    !         IF ( NBI .EQ. 0 ) THEN
    !             RETURN
    !           ELSE
    !             IF ( MAXVAL(BBPI0(ISP,:)) .EQ. 0. .AND.            &
    !                  MAXVAL(BBPIN(ISP,:)) .EQ. 0. ) RETURN
    !           END IF
    !       END IF
    !
    ! 1.  Preparations --------------------------------------------------- *

    IF ( ICLOSE .EQ. ICLOSE_TRPL ) THEN
      IF (IAPROC .EQ. NAPERR) &
           WRITE(NDSE,*)'SUBROUTINE W3XYP2 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
    !
    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    = 0.9999 * 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)
    !
    IF ( FLAGLL ) THEN
      RFAC = DERA * RADIUS
      DFAC = 1. / RFAC**2
    ELSE
      RFAC = 1.
      DFAC = 1.
    END IF
    !
    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
    !
#ifdef W3_TDYN
    IF ( ISP .EQ. 1 ) DTME = DTME + DTG
#endif
    !
    IF ( DTME .NE. 0. ) THEN
      DFRR   = XFR - 1.
      CELLP  = 10.
      CGD    = 0.5 * GRAV / SIG(IK)
      DSSD   = ( DFRR * CGD )**2 * DTME / 12.
      DNND   = ( CGD * DTH )**2 * DTME / 12.
#ifdef W3_T
      WRITE (NDST,9002) DFRR, CELLP, DTME
    ELSE
      WRITE (NDST,9003)
#endif
    END IF

    !
    ! 1.b Initialize arrays
    !
#ifdef W3_T
    WRITE (NDST,9010)
#endif
    !
    VLCFLX = 0.
    VLCFLY = 0.
    VFDIFX = 0.
    VFDIFY = 0.
    VFDIFC = 0.
    VDXX   = 0.
    VDYY   = 0.
    VDXY   = 0.
    CXTOT  = 0.
    CYTOT  = 0.
    !
    ! 2.  Calculate velocities and diffusion coefficients ---------------- *
    ! 2.a Velocities
    !
    !     Q     = ( A / CG * CLATS )
    !     LCFLX = ( COS*CG / CLATS ) * DT / DX
    !     LCFLY = (     SIN*CG )     * DT / DX
    !
#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
      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
    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
    !
    ! 2.b Diffusion coefficients
    !
    IF ( DTME .NE. 0. ) THEN
      !
#ifdef W3_OMPH
      !$OMP PARALLEL DO PRIVATE (ISEA, IX, IY, IXY, &
      !$OMP&                     DCELL, XWIND, TFAC, DSS, DNN)
#endif
      !
      DO ISEA=1, NSEA
        IX        = MAPSF(ISEA,1)
        IY        = MAPSF(ISEA,2)
        IXY       = MAPSF(ISEA,3)
        IF ( MIN ( ATRNX(IXY,1) , ATRNX(IXY,-1) ,                 &
             ATRNY(IXY,1) , ATRNY(IXY,-1) ) .GT. TRNMIN ) THEN
          DCELL     = CGD * MIN ( HPFAC(IY,IX)*RFAC, &
               HQFAC(IY,IX)*RFAC ) / CELLP
          XWIND     = 3.3 * U10(ISEA)*WN(IK,ISEA)/SIG(IK) - 2.3
          XWIND     = MAX ( 0. , MIN ( 1. , XWIND ) )
#ifdef W3_XW0
          XWIND     = 0.
#endif
#ifdef W3_XW1
          XWIND     = 1.
#endif
          TFAC      = MIN ( 1. , (CLATS(ISEA)/CLATMN)**2 )
          DSS       = XWIND * DCELL + (1.-XWIND) * DSSD * TFAC
#ifdef W3_DSS0
          DSS       = 0.
#endif
          DNN       = XWIND * DCELL + (1.-XWIND) * DNND * TFAC

          VDXX(IXY) = DTLOC * (DSS*ECOS(ITH)**2+DNN*ESIN(ITH)**2)
          VDYY(IXY) = DTLOC * (DSS*ESIN(ITH)**2+DNN*ECOS(ITH)**2) &
               / CLATS(ISEA)**2
          VDXY(IXY) = DTLOC * (DSS-DNN) * ESIN(ITH)*ECOS(ITH)     &
               / CLATS(ISEA)

        END IF
      END DO
      !
#ifdef W3_OMPH
      !$OMP END PARALLEL DO
#endif
      !
    END IF
    !
    ! 3.  Loop over sub-steps -------------------------------------------- *
    !
    DO ITLOC=1, NTLOC
      !
      ! 3.a Propagate fields
      !
#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
      !
#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.b Update boundaries
      !
      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)
          IXY     = MAPSF(ISBPI(IBI),3)
          VQ(IXY) = ( RD1*BBPI0(ISP,IBI) + RD2*BBPIN(ISP,IBI) )   &
               / CG(IK,ISEA) * CLATS(ISEA)
        END DO
      END IF
      !
      ! 3.c Diffusion correction
      !
      IF ( DTME .NE. 0. ) THEN

        IF ( GLOBAL ) THEN
          DO IY=1, NY
            VQ(IY+NX*NY) = VQ(IY)
          END DO
        END IF

        !CURV2: BEGIN -----------------------------------------------------------------
        VFDIFX_FAC=0.0
        DO IP=1, NMX0
          IXY         = MAPX2(IP)
          VFDIFX_FAC(IXY) = 1.0
        END DO
        VFDIFY_FAC=0.0
        DO IP=1, NMY0
          IXY         = MAPY2(IP)
          VFDIFY_FAC(IXY) = 1.0
        END DO
        VFDIFC_FAC=0.0
        DO IP=1, NMXY
          IXY         = MAPXY(IP)
          VFDIFC_FAC(IXY) = 1.0
        END DO
        IF ( GLOBAL ) THEN
          IY0     = (NX-1)*NY
          DO IY=1, NY
            VFDIFX_FAC(IY-NY) = VFDIFX_FAC(IY+IY0)
            VFDIFC_FAC(IY-NY) = VFDIFC_FAC(IY+IY0)
          END DO
        END IF
        !CURV2: END -------------------------------------------------------------------

        VQ_OLD = VQ
        !
#ifdef W3_OMPH
        !$OMP PARALLEL DO PRIVATE (ISEA, IX, IY, IXY, &
        !$OMP&                     QXX, QYY, QXY, DVQ )
#endif
        !
        DO IP=1, NACT
          IXY    = MAPAXY(IP)
          ISEA   = MAPFS(IXY)
          IX     = MAPSF(ISEA,1)
          IY     = MAPSF(ISEA,2)

          !CURV1: DOES NOT GIVE EXACT MATCH TO EARLIER WW3 VERSION FOR TEST CASE TP2.3
          !CURV1: NEAR THE BOUNDARY, NOTE THAT THIS VERSION USES NON-ACTIVE GRID POINTS
          !CURV1: IN ITS CALCS. THIS IS NO PROBLEM, AS LONG AS THE NON-ACTIVE GRID POINTS
          !CURV1: EXIST IN THE ARRAY AND ARE VQ=0. ALSO NOTE THAT WITH THE SHORT VERSION
          !CURV1: OF THE CODE, VFDIF?_FAC VARIABLES CAN BE REMOVED AND 3-4 DO LOOPS ABOVE
          !CURV1: CAN BE REMOVED.
          !CURV1: BEGIN -----------------------------------------------------------------
          !                QXX = VQ_OLD(IXY+NY) - 2.0*VQ_OLD(IXY) + VQ_OLD(IXY-NY)
          !                QYY = VQ_OLD(IXY+1)  - 2.0*VQ_OLD(IXY) + VQ_OLD(IXY-1)
          !                QXY = VQ_OLD(IXY+NY+1) - VQ_OLD(IXY-NY+1) &
          !                    - VQ_OLD(IXY+NY-1) + VQ_OLD(IXY-NY-1)
          !CURV1: END -------------------------------------------------------------------

          !CURV2: DOES GIVE EXACT MATCH TO EARLIER WW3 VERSION FOR TEST CASE TP2.3. NOTE
          !CURV2: THAT IF VFDIFC_FAC VARIABLES ARE ALL UNITY, MANY TERMS CANCEL OUT.
          !CURV2: HOWEVER, VFDIFC_FAC IS ZERO WHEN A RELATED VQ POINT IS NOT AN ACTIVE
          !CURV2: GRID POINT
          !CURV2: BEGIN -----------------------------------------------------------------
          QXX = VFDIFX_FAC(IXY)   *VQ_OLD(IXY+NY)     &
               - VFDIFX_FAC(IXY)   *VQ_OLD(IXY)        &
               - VFDIFX_FAC(IXY-NY)*VQ_OLD(IXY)        &
               + VFDIFX_FAC(IXY-NY)*VQ_OLD(IXY-NY)
          QYY = VFDIFY_FAC(IXY)   *VQ_OLD(IXY+1)      &
               - VFDIFY_FAC(IXY)   *VQ_OLD(IXY)        &
               - VFDIFY_FAC(IXY-1) *VQ_OLD(IXY)        &
               + VFDIFY_FAC(IXY-1) *VQ_OLD(IXY-1)
          QXY = VFDIFC_FAC(IXY)     *VQ_OLD(IXY)      &
               + VFDIFC_FAC(IXY-NY-1)*VQ_OLD(IXY)      &
               - VFDIFC_FAC(IXY-1)   *VQ_OLD(IXY)      &
               - VFDIFC_FAC(IXY-NY)  *VQ_OLD(IXY)      &
               + VFDIFC_FAC(IXY-NY)  *VQ_OLD(IXY+1)    &
               - VFDIFC_FAC(IXY)     *VQ_OLD(IXY+1)    &
               + VFDIFC_FAC(IXY-1)   *VQ_OLD(IXY-1)    &
               - VFDIFC_FAC(IXY-NY-1)*VQ_OLD(IXY-1)    &
               + VFDIFC_FAC(IXY-1)   *VQ_OLD(IXY+NY)   &
               - VFDIFC_FAC(IXY)     *VQ_OLD(IXY+NY)   &
               + VFDIFC_FAC(IXY-NY)  *VQ_OLD(IXY-NY)   &
               - VFDIFC_FAC(IXY-NY-1)*VQ_OLD(IXY-NY)   &
               + VFDIFC_FAC(IXY)     *VQ_OLD(IXY+NY+1) &
               - VFDIFC_FAC(IXY-1)   *VQ_OLD(IXY+NY-1) &
               + VFDIFC_FAC(IXY-NY-1)*VQ_OLD(IXY-NY-1) &
               - VFDIFC_FAC(IXY-NY)  *VQ_OLD(IXY-NY+1)
          !CURV2: END -------------------------------------------------------------------
          !
          QXY = 0.25*QXY
          !
          DVQ = VDXX(IXY)*(  DPDX(IY,IX)*DPDX(IY,IX)*QXX          &
               + 2.0*DPDX(IY,IX)*DQDX(IY,IX)*QXY      &
               + DQDX(IY,IX)*DQDX(IY,IX)*QYY )        &
               + VDYY(IXY)*(  DPDY(IY,IX)*DPDY(IY,IX)*QXX          &
               + 2.0*DPDY(IY,IX)*DQDY(IY,IX)*QXY      &
               + DQDY(IY,IX)*DQDY(IY,IX)*QYY)         &
               + 2.0*VDXY(IXY)*(  DPDX(IY,IX)*DPDY(IY,IX)*QXX      &
               + DQDX(IY,IX)*DQDY(IY,IX)*QYY      &
               +      (  DPDX(IY,IX)*DQDY(IY,IX)  &
               + DQDX(IY,IX)*DPDY(IY,IX) )*QXY )
          !
          VQ(IXY) = VQ_OLD(IXY) + DVQ * DFAC
          !
        END DO
        !
#ifdef W3_OMPH
        !$OMP END PARALLEL DO
#endif
        !
      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_T2
        WRITE (NDST,9041) ISEA, MAPSF(ISEA,1), MAPSF(ISEA,2), VQ(IXY)
#endif
        VQ(IXY) =  MAX ( 0. , CG(IK,ISEA) / CLATS(ISEA) * VQ(IXY) )
        !         ELSE
        !           VQ(IXY) =  0.
      END IF
    END DO
    !
#ifdef W3_OMPH
    !$OMP END PARALLEL DO
#endif
    !
    RETURN
    !
    ! Formats
    !
#ifdef W3_T
9000 FORMAT (' TEST W3XYP2 : YFIRST  :',L2)
9001 FORMAT (' TEST W3XYP2 : ISP, ITH, IK, COS-SIN :',I8,2I4,2F7.3)
9002 FORMAT (' TEST W3XYP2 : DFRR, CELLP, DTME     :',3E10.3)
9003 FORMAT (' TEST W3XYP2 : NO DISPERSION CORRECTION ')
9010 FORMAT (' TEST W3XYP2 : INITIALIZE ARRAYS')
9020 FORMAT (' TEST W3XYP2 : CALCULATING LCFLX/Y AND DSS/NN (NSEA=', &
         I6,')')
9022 FORMAT (' TEST W3XYP2 : CORRECTING FOR CURRENT')
9040 FORMAT (' TEST W3XYP2 : FIELD AFTER PROP. (NSEA=',I6,')')
#endif
#ifdef W3_T1
9021 FORMAT (1X,I6,2I5,E12.4,2f7.3)
#endif
#ifdef W3_T2
9041 FORMAT (1X,I6,2I5,E12.4)
#endif
    !/
    !/ End of W3XYP2 ----------------------------------------------------- /
    !/
  END SUBROUTINE W3XYP2
  !/
  !/ ------------------------------------------------------------------- /
  !>
  !> @brief Propagation in spectral space.
  !>
  !> @details 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
  !>  represent 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.
  !>
  !> @param[in]    ISEA   Number of sea point.
  !> @param[in]    FACTH  Factor in propagation velocity.
  !> @param[in]    FACK   Factor in propagation velocity.
  !> @param[in]    CTHG0  Factor in great circle refracftion term.
  !> @param[in]    CG     Local group velocities.
  !> @param[in]    WN     Local wavenumbers.
  !> @param[in]    DEPTH  Depth.
  !> @param[in]    DDDX   Depth gradient.
  !> @param[in]    DDDY   Depth gradient.
  !> @param[in]    CX     Current component.
  !> @param[in]    CY     Current component.
  !> @param[in]    DCXDX  Current gradients.
  !> @param[in]    DCXDY  Current gradients.
  !> @param[in]    DCYDX  Current gradients.
  !> @param[in]    DCYDY  Current gradients.
  !> @param[in]    DCDX   Phase speed gradient.
  !> @param[in]    DCDY   Phase speed gradient.
  !> @param[inout] VA     Spectrum.
  !>
  !> @author H. L. Tolman
  !> @date   01-Jul-2013
  !>
  SUBROUTINE W3KTP2 ( ISEA, FACTH, FACK, CTHG0, CG, WN, DEPTH,    &
       DDDX, DDDY, CX, CY, DCXDX, DCXDY,           &
       DCYDX, DCYDY, DCDX, DCDY, VA )
    !/
    !/                  +-----------------------------------+
    !/                  | 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 )
    !/    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.
    !       DEPTH   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
    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), DEPTH, DDDX, DDDY,       &
         CX, CY, DCXDX, DCXDY, DCYDX, DCYDY
    REAL, INTENT(IN)        :: DCDX(0:NK+1), DCDY(0:NK+1)
    REAL, INTENT(INOUT)     :: VA(NSPEC)
    !/
    !/ ------------------------------------------------------------------- /
    !/ 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
    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, 'W3KTP2')
#endif
    !
    ! 1.  Preparations --------------------------------------------------- *
    ! 1.a Initialize arrays
    !
    IF ( FLCK ) VQ = 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
      !
      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)
        FRK(IK) = FRK(IK) / MAX ( 1. , FRK(IK)*FDDMAX/CTMAX )
        FRG(IK) = FDG * CG(IK)
      END DO
      !
      ! 3.b Depth refraction and great-circle propagation
      !
      DO ISP=1, NSPEC
        VCFLT(MAPTH2(ISP)) = FRG(MAPWN(ISP)) * ECOS(ISP)          &
             + FRK(MAPWN(ISP)) * ( ESIN(ISP)*DDDX - ECOS(ISP)*DDDY )
      END DO
      !
#ifdef W3_REFRX
      ! 3.c @C/@x refraction and great-circle propagation
      VCFLT  = 0.
      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
      DO ISP=1, NSPEC
        VCFLT(MAPTH2(ISP)) = FRG(MAPWN(ISP)) * ECOS(ISP)   &
             + FRK(MAPWN(ISP)) * ( ESIN(ISP)*DCDX(MAPWN(ISP)) &
             - ECOS(ISP)*DCDY(MAPWN(ISP)) )
      END DO
#endif
      !
      ! 3.d Current refraction
      !
      IF ( FLCUR ) THEN
        !
        DCYX   = FACTH *   DCYDX
        DCXXYY = FACTH * ( DCXDX - DCYDY )
        DCXY   = FACTH *   DCXDY
        !
        DO ISP=1, NSPEC
          VCFLT(MAPTH2(ISP)) = VCFLT(MAPTH2(ISP)) +             &
               ES2(ISP)*DCYX  + ESC(ISP)*DCXXYY - EC2(ISP)*DCXY
        END DO
        !
      END IF
      !
    END IF
    !
    ! 4.  Wavenumber shift velocities ------------------------------------ *
    !     FACK is just the time step, which is accounted for in W3QCK2
    !
    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 Velocities
      !
      DO IK=0, NK+1
        FKD0   = FKD / CG(IK) * DSDD(IK)
        DO ITH=1, NTH
          CFLK(IK+1,ITH) = FKD0 + WN(IK)*FKC(ITH)
        END DO
      END DO
      !
      ! 4.c 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
      !
    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 W3KTP2 : FLCTH-K, FACTH-K, CTMAX  :',        &
         2L2,2E10.3,F7.3)
9010 FORMAT ( ' TEST W3KTP2 : LOCAL DATA :',I7,F7.1,2F6.2,1X,6E10.2)
9020 FORMAT ( ' TEST W3KTP2 : IK, T, L, CG, DSDD : ')
9021 FORMAT ( '              ',I3,F7.2,F7.1,F7.2,E11.3)
#endif
    !
#ifdef W3_T0
9040 FORMAT (/' TEST W3KTP2 : NORMALIZED ',A/)
9041 FORMAT (1X,60(1X,I2))
9042 FORMAT (1X,60I3)
#endif
    !/
    !/ End of W3KTP2 ----------------------------------------------------- /
    !/
  END SUBROUTINE W3KTP2
  !/
  !/ End of module W3PRO2MD -------------------------------------------- /
  !/
END MODULE W3PRO2MD
