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

#include "w3macros.h"
!/ ------------------------------------------------------------------- /
!>
!> @brief Bundles routines for third order propagation scheme in single
!>  module.
!>
!> @author H. L. Tolman
!> @date   27-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 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
  !/ ------------------------------------------------------------------- /
  !>
  !> @brief Generate 'map' arrays for the ULTIMATE QUICKEST scheme.
  !>
  !> @author H. L. Tolman
  !> @date   01-Apr-2008
  !>
  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
  !/ ------------------------------------------------------------------- /
  !>
  !> @brief Generate 'map' arrays for the ULTIMATE QUICKEST scheme to combine
  !>  GSE alleviation with obstructions.
  !>
  !> @author H. L. Tolman
  !> @date   17-Dec-2004
  !>
  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
  !/ ------------------------------------------------------------------- /
  !>
  !> @brief Propagation in phyiscal 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     Speed of grid.
  !> @param[in]    VGY     Speed of grid.
  !>
  !> @author H. L. Tolman
  !> @date   27-May-2014
  !>
  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
  !/ ------------------------------------------------------------------- /
!>
!> @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]    DW        Depth.
!> @param[in]    DDDX      Depth gradients.
!> @param[in]    DDDY      Depth gradients.
!> @param[in]    CX        Current components.
!> @param[in]    CY        Current components.
!> @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 gradients.
!> @param[in]    DCDY      Phase speed gradients.
!> @param[inout] VA        Spectrum.
!> @param[out]   CFLTHMAX
!> @param[out]   CFLKMAX
!>
!> @author H. L. Tolman
!> @date   01-Jul-2013
!>
  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
  !/ ------------------------------------------------------------------- /
  !>
  !> @brief Computes the maximum CFL number for spatial advection.
  !>
  !> @details Used for diagnostic purposes (Could be used to define a
  !>  local time step ...).
  !>
  !> @param[in]    ISEA      Index of grid point.
  !> @param[in]    DTG       Total time step.
  !> @param[in]    MAPSTA    Grid point status map.
  !> @param[in]    MAPFS     Storage map.
  !> @param[inout] CFLXYMAX  Maximum CFL number for XY propagation.
  !> @param[in]    VGX       Speed of grid.
  !> @param[in]    VGY       Speed of grid.
  !>
  !> @author F. Ardhuin
  !> @date   31-Oct-2010
  !>
  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
