text: !/ ------------------------------------------------------------------- /
        MODULE W3PRO3MD
  !/
  !/                  +-----------------------------------+
  !/                  | WAVEWATCH III           NOAA/NCEP |
  !/                  |           H. L. Tolman            |
  !/                  |                        FORTRAN 90 |
  !/                  | Last update :         29-May-2009 |
  !/                  +-----------------------------------+
  !/
  !/    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 )
  !/
  !/    Copyright 2009 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
  !                                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. W3UQCKMD Irregular grid UQ scheme.
  !      W3QCK3    Subr. W3UQCKMD Regular grid UQ scheme + obstructions.
  !     ----------------------------------------------------------------
  !
  !  5. Remarks :
  !
  !     - The averaging is not performed around semi-transparent grid
  !       points to avoid that leaking through barriers occurs.
  !
  !  6. Switches :
  !
  !       !/C90   Cray FORTRAN 90 compiler directives.
  !       !/NEC   NEC SXF90 compiler directives.
  !
  !       !/MGP   Moving grid corrections.
  !       !/MGG   Moving grid corrections.
  !
  !       !/S     Enable subroutine tracing.
  !       !/Tn    Enable test output.
  !
  !  7. Source code :
  !
  !/ ------------------------------------------------------------------- /
  !/
  !/ Public variables
  !/
        PUBLIC
  !/
  !/ Private data
  !/
        REAL, PRIVATE, PARAMETER:: TRNMIN = 0.95
  !/
        CONTAINS
  !/ ------------------------------------------------------------------- /
        SUBROUTINE W3MAP3
  !/
  !/                  +-----------------------------------+
  !/                  | WAVEWATCH III           NOAA/NCEP |
  !/                  |           H. L. Tolman            |
  !/                  |                        FORTRAN 90 |
  !/                  | Last update :         01-Apr-2008 |
  !/                  +-----------------------------------+
  !/
  !/    27-Feb-2000 : Origination.                        ( version 2.08 )
  !/    10-Dec-2001 : Sub-grid obstructions.              ( version 2.14 )
  !/                  (array allocation only.)
  !/    17-Dec-2004 : Multiple grid version.              ( version 3.06 )
  !/    09-Nov-2005 : Removing soft boundary option.      ( version 3.08 )
  !/    01-Apr-2008 : Bug fix sec. 4 MAPSTA range check.  ( version 3.13 )
  !/
  !  1. Purpose :
  !
  !     Generate 'map' arrays for the ULTIMATE QUICKEST scheme.
  !
  !  2. Method :
  !
  !     MAPX2, MAPY2, MAPTH2 and MAPWN2 contain consecutive 1-D spatial
  !     grid counters (e.g., IXY = (IX-1)*MY + IY). The arrays are
  !     devided in three parts. For MAPX2, these ranges are :
  !
  !         1    - NMX0  Counters IXY for which grid point (IX,IY) and
  !                      (IX+1,IY) both are active grid points.
  !       NMX0+1 - NMX1  Id. only (IX,IY) active.
  !       NMX1+1 - NMX2  Id. only (IX+1,IY) active.
  !
  !     The array MAPY2 has a similar layout varying IY instead of IX.
  !
  !  3. Parameters :
  !
  !     Parameter list
  !     ----------------------------------------------------------------
  !     ----------------------------------------------------------------
  !
  !  4. Subroutines used :
  !
  !     See module documentation.
  !
  !  5. Called by :
  !
  !      Name      Type  Module   Description
  !     ----------------------------------------------------------------
  !      W3WAVE    Subr. W3WAVEMD Wave model routine.
  !     ----------------------------------------------------------------
  !
  !  6. Error messages :
  !
  !  7. Remarks :
  !
  !  8. Structure :
  !
  !     ------------------------------------------------------
  !      1.   Map MAPX2
  !        a  Range 1 to NMX0
  !        b  Range NMX0+1 to NMX1
  !        c  Range NMX1+1 to NMX2
  !      2.   Map MAPY2
  !        a  Range 1 to NMY0
  !        b  Range NMY0+1 to NMY1
  !        c  Range NMY1+1 to NMY2
  !      3.   Map MAPAXY
  !      4.   Map MAPCXY
  !      5.   Maps for intra-spectral propagation
  !        a  MAPTH2, MAPATK
  !        b  MAPWN2
  !     ------------------------------------------------------
  !
  !  9. Switches :
  !
  !     !/S   Enable subroutine tracing.
  !     !/T   Enable test output.
  !
  ! 10. Source code :
  !/ ------------------------------------------------------------------- /
        USE W3GDATMD, ONLY: NK, NTH, NSPEC, NX, NY, NSEA, MAPSTA, MAPSF
        USE W3ADATMD, ONLY: NMX0, NMX1, NMX2, NMY0, NMY1, NMY2, NACT,   &
                            NCENT, MAPX2, MAPY2, MAPAXY, MAPCXY,        &
                            MAPTH2, MAPWN2
        USE W3ODATMD, ONLY: NDST
  !/
        IMPLICIT NONE
  !/
  !/ ------------------------------------------------------------------- /
  !/ Parameter list
  !/
  !/ ------------------------------------------------------------------- /
  !/ Local parameters
  !/
        INTEGER                 :: IX, IY, IXY0, IX2, IY2, IX0, IY0,    &
                                   ISEA, IK, ITH, ISP, ISP0, ISP2, NCENTC
  !/
  !/ ------------------------------------------------------------------- /
  !/
  !
  ! 1.  Map MAPX2 ------------------------------------------------------ *
  ! 1.a Range 1 to NMX0
  !
        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
              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
              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
              END IF
            END DO
          END DO
  !
  ! 2.  Map MAPY2 ------------------------------------------------------ *
  ! 2.a Range 1 to NMY0
  !
        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
              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
              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
              END IF
             END DO
           END DO
  !
  ! 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
  !
  ! 5.  Maps for intra-spectral propagation ---------------------------- *
  !
        IF ( MAPTH2(1) .NE. 0 ) RETURN
  !
  ! 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
            END DO
          END DO
  !
  ! 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
            END DO
          END DO
  !
        DO ITH=1, NTH
          ISP0   = ISP0 + 1
          ISP2   = NK+1 + (ITH-1)*(NK+2)
          MAPWN2(ISP0) = ISP2
          END DO
  !
        DO ITH=1, NTH
          ISP0   = ISP0 + 1
          ISP2   = 1 + (ITH-1)*(NK+2)
          MAPWN2(ISP0) = ISP2
          END DO
  !
        RETURN
  !
  ! Formats
  !
  !/
  !/ End of W3MAP3 ----------------------------------------------------- /
  !/
        END SUBROUTINE W3MAP3
  !/ ------------------------------------------------------------------- /
        SUBROUTINE W3MAPT
  !/
  !/                  +-----------------------------------+
  !/                  | WAVEWATCH III           NOAA/NCEP |
  !/                  |           H. L. Tolman            |
  !/                  |                        FORTRAN 90 |
  !/                  | Last update :         17-Dec-2004 |
  !/                  +-----------------------------------+
  !/
  !/    10-Dec-2001 : Origination.                        ( version 2.14 )
  !/    17-Dec-2004 : Multiple grid version.              ( version 3.06 )
  !/
  !  1. Purpose :
  !
  !     Generate 'map' arrays for the ULTIMATE QUICKEST scheme to combine
  !     GSE alleviation with obstructions.
  !
  !  2. Method :
  !
  !  3. Parameters :
  !
  !     Parameter list
  !     ----------------------------------------------------------------
  !     ----------------------------------------------------------------
  !
  !  4. Subroutines used :
  !
  !     See module documentation.
  !
  !  5. Called by :
  !
  !      Name      Type  Module   Description
  !     ----------------------------------------------------------------
  !      W3WAVE    Subr. W3WAVEMD Wave model routine.
  !     ----------------------------------------------------------------
  !
  !  6. Error messages :
  !
  !  7. Remarks :
  !
  !  8. Structure :
  !
  !     See source code.
  !
  !  9. Switches :
  !
  !     !/S   Enable subroutine tracing.
  !
  ! 10. Source code :
  !/ ------------------------------------------------------------------- /
        USE W3GDATMD, ONLY: NX, NY, NSEA, MAPSF
        USE W3ADATMD, ONLY: ATRNX, ATRNY, MAPTRN
  !/
        IMPLICIT NONE
  !/
  !/ ------------------------------------------------------------------- /
  !/ Parameter list
  !/
  !/ ------------------------------------------------------------------- /
  !/ Local parameters
  !/
        INTEGER                 :: ISEA, IXY
  !/
  !/ ------------------------------------------------------------------- /
  !/
  !
  ! 1.  Map MAPTRN ----------------------------------------------------- *
  !
        DO ISEA=1, NSEA
          IXY         = MAPSF(ISEA,3)
          MAPTRN(IXY) = MIN( ATRNX(IXY,1) ,ATRNY(IXY,-1) ,              &
                             ATRNY(IXY,1), ATRNY(IXY,-1) ) .LT. TRNMIN
          END DO
  !
        RETURN
  !
  ! Formats
  !/
  !/ End of W3MAPT ----------------------------------------------------- /
  !/
        END SUBROUTINE W3MAPT
  !/ ------------------------------------------------------------------- /
        SUBROUTINE W3XYP3 ( ISP, FACX, FACY, DTG, MAPSTA, MAPFS, VQ,    &
                            VGX, VGY )
  !/
  !/                  +-----------------------------------+
  !/                  | WAVEWATCH III           NOAA/NCEP |
  !/                  |           H. L. Tolman            |
  !/                  |                        FORTRAN 90 |
  !/                  | Last update :         05-Mar-2008 |
  !/                  +-----------------------------------+
  !/
  !/    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 )
  !/
  !  1. Purpose :
  !
  !     Propagation in phyiscal space for a given spectral component.
  !
  !  2. Method :
  !
  !     Third-order ULTIMATE QUICKEST scheme with averaging.
  !
  !  3. Parameters :
  !
  !     Parameter list
  !     ----------------------------------------------------------------
  !       ISP     Int.   I   Number of spectral bin (IK-1)*NTH+ITH
  !       FACX/Y  Real   I   Factor in propagation velocity.
  !                          ( 1 or 0 * DT / DX )
  !       DTG     Real   I   Total time step.
  !       MAPSTA  I.A.   I   Grid point status map.
  !       MAPFS   I.A.   I   Storage map.
  !       VQ      R.A.  I/O  Field to propagate.
  !       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.
  !
  !  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 :
  !
  !       !/NEC   Enable NEC SXF90 compiler directives.
  !
  !       !/S     Enable subroutine tracing.
  !
  !       !/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.
  !
  ! 10. Source code :
  !
  !/ ------------------------------------------------------------------- /
        USE W3CONSTANTS
  !
        USE W3TIMEMD, ONLY: DSEC21
  !
        USE W3GDATMD, ONLY: NX, NY, NSEA, MAPSF, DTCFL, CLAT, CLATS,    &
                            GLOBAL, FLCX, FLCY, NK, NTH, DTH, XFR,      &
                            ECOS, ESIN, SIG, WDCG, WDTH, PFMOVE
        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
        USE W3UQCKMD
  !/
        IMPLICIT NONE
  !/
  !/ ------------------------------------------------------------------- /
  !/ Parameter list
  !/
        INTEGER, INTENT(IN)     :: ISP, MAPSTA(NY*NX), MAPFS(NY*NX)
        REAL, INTENT(IN)        :: FACX, FACY, 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, IIXY1(NY),         &
                                   IIXY2(NY), IIXY3(NY), IIXY4(NY), IBI
        REAL                    :: CG0, CGA, CGN, CGX, CGY, CXC, CYC,   &
                                   CXMIN, CXMAX, CYMIN, CYMAX
        REAL                    :: CGC, FGSE = 1., FTH, FCG
        REAL                    :: DTLOC, CCOS, CSIN, CCURX, CCURY
        REAL                    :: DXCGN, DYCGN, DXCGS, DYCGS, DXCGC,   &
                                   DYCGC, RDI1(NY), RDI2(NY), RDI3(NY), &
                                   RDI4(NY), RD1, RD2, RD3, RD4
        LOGICAL                 :: YFIRST
  !/
  !/ 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))
  !/
  !/ ------------------------------------------------------------------- /
  !/
  !
  ! 1.  Preparations --------------------------------------------------- *
  ! 1.a Set constants
  !
        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)
        CGC    = SQRT ( CGX**2 + CGY**2 )
  !
        IF ( FLCUR ) THEN
            CXMIN  = MINVAL ( CX(1:NSEA) )
            CXMAX  = MAXVAL ( CX(1:NSEA) )
            CYMIN  = MINVAL ( CY(1:NSEA) )
            CYMAX  = MAXVAL ( CY(1:NSEA) )
            IF ( ABS(CGX+CXMIN) .GT. ABS(CGX+CXMAX) ) THEN
                CGX    = CGX + CXMIN
              ELSE
                CGX    = CGX + CXMAX
              END IF
            IF ( ABS(CGY+CYMIN) .GT. ABS(CGY+CYMAX) ) THEN
                CGY    = CGY + CYMIN
              ELSE
                CGY    = CGY + CYMAX
              END IF
            CXC    = MAX ( ABS(CXMIN) , ABS(CXMAX) )
            CYC    = MAX ( ABS(CYMIN) , ABS(CYMAX) )
          ELSE
            CXC    = 0.
            CYC    = 0.
          END IF
  !
        CGN    = MAX ( ABS(CGX) , ABS(CGY) , CXC, CYC, 0.001*CG0 )
  !
        NTLOC  = 1 + INT(DTG/(DTCFL*CG0/CGN))
        DTLOC  = DTG / REAL(NTLOC)
  !
        CCOS   = FACX * ECOS(ITH) / REAL(NTLOC)
        CSIN   = FACY * ESIN(ITH) / REAL(NTLOC)
  !
        CCURX  = FACX / REAL(NTLOC)
        CCURY  = FACY / REAL(NTLOC)
  !
        YFIRST = MOD(ITIME,2) .EQ. 0
  !
  ! 1.b Initialize arrays
  !
        VLCFLX = 0.
        VLCFLY = 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 / REAL(NTLOC)
        FCG    =   FGSE * WDCG * 0.5 * (XFR-1./XFR) / REAL(NTLOC)
  !
        DO IY=1, NY
  !
  ! 1.c.1 Normal and parallel width ...
  !
          DXCGN  = - FTH * FACX * ESIN(ITH) / CLAT(IY)
          DYCGN  =   FTH * FACY * ECOS(ITH)
          DXCGS  =   FCG * CCOS / CLAT(IY)
          DYCGS  =   FCG * CSIN
  !
  ! 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(IY) = IXC + IYC
          IF ( ABS(DXCGC) .GT. ABS(DYCGC) ) THEN
              IIXY2(IY) = IXC
              RDI1 (IY) = ABS(DYCGC/DXCGC)
              RDI2 (IY) = ABS(DXCGC)
            ELSE
              IIXY2(IY) = IYC
              IF ( ABS(DYCGC) .GT. 1.E-5 ) THEN
                  RDI1(IY) = ABS(DXCGC/DYCGC)
                ELSE
                  RDI1(IY) = 1.
                END IF
              RDI2(IY) = ABS(DYCGC)
            END IF
  !
  ! 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(IY) = IXC + IYC
          IF ( ABS(DXCGC) .GT. ABS(DYCGC) ) THEN
              IIXY4(IY) = IXC
              RDI3 (IY) = ABS(DYCGC/DXCGC)
              RDI4 (IY) = ABS(DXCGC)
            ELSE
              IIXY4(IY) = IYC
              IF ( ABS(DYCGC) .GT. 1.E-5 ) THEN
                  RDI3(IY) = ABS(DXCGC/DYCGC)
                ELSE
                  RDI3(IY) = 1.
                END IF
              RDI4(IY) = ABS(DYCGC)
            END IF
  !
          END DO
  !
  ! 2.  Calculate velocities and diffusion coefficients ---------------- *
  ! 2.a Velocities
  !
  !     Q     = ( A / CG * CLATS )
  !     LCFLX = ( COS*CG / CLATS ) * DT / DX
  !     LCFLY = (     SIN*CG )     * DT / DY
  !
        DO ISEA=1, NSEA
          IXY         = MAPSF(ISEA,3)
          VQ    (IXY) = VQ(IXY) / CG(IK,ISEA) * CLATS(ISEA)
          VLCFLX(IXY) = CCOS * CG(IK,ISEA) / CLATS(ISEA)
          VLCFLY(IXY) = CSIN * CG(IK,ISEA)
          END DO
  !
        IF ( FLCUR ) THEN
            DO ISEA=1, NSEA
              IXY         = MAPSF(ISEA,3)
              VLCFLX(IXY) = VLCFLX(IXY) + CCURX*CX(ISEA)/CLATS(ISEA)
              VLCFLY(IXY) = VLCFLY(IXY) + CCURY*CY(ISEA)
             END DO
          END IF
  !
  ! 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)
            IY      = MAPSF(ISEA,2)
            IF ( MAPTRN(IXY) ) THEN
                VQ(IXY) = AQ(IXY)
              ELSE
                RD1     = RDI1(IY)
                RD2     = MIN ( 1. , RDI2(IY) * CG(IK,ISEA) )
                RD3     = RDI3(IY)
                RD4     = MIN ( 1. , RDI4(IY) * CG(IK,ISEA) )
                VQ(IXY          ) = VQ(IXY          )                   &
                                     + AQ(IXY) * (3.-RD2-RD4)/3.
                VQ(IXY+IIXY1(IY)) = VQ(IXY+IIXY1(IY))                   &
                                     + AQ(IXY) * RD2*RD1/6.
                VQ(IXY+IIXY2(IY)) = VQ(IXY+IIXY2(IY))                   &
                                     + AQ(IXY) * (1.-RD1)*RD2/6.
                VQ(IXY+IIXY3(IY)) = VQ(IXY+IIXY3(IY))                   &
                                     + AQ(IXY) * RD4*RD3/6.
                VQ(IXY+IIXY4(IY)) = VQ(IXY+IIXY4(IY))                   &
                                     + AQ(IXY) * (1.-RD3)*RD4/6.
                VQ(IXY-IIXY1(IY)) = VQ(IXY-IIXY1(IY))                   &
                                     + AQ(IXY) * RD2*RD1/6.
                VQ(IXY-IIXY2(IY)) = VQ(IXY-IIXY2(IY))                   &
                                     + AQ(IXY) * (1.-RD1)*RD2/6.
                VQ(IXY-IIXY3(IY)) = VQ(IXY-IIXY3(IY))                   &
                                     + AQ(IXY) * RD4*RD3/6.
                VQ(IXY-IIXY4(IY)) = VQ(IXY-IIXY4(IY))                   &
                                     + 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)
            IY      = MAPSF(ISEA,2)
            IXY     = MAPSF(ISEA,3)
            IF ( MAPSTA(IXY) .LE. 0 ) CYCLE
            IF ( MAPTRN(IXY) ) THEN
                VQ(IXY) = AQ(IXY)
              ELSE
                RD1     = RDI1(IY)
                RD3     = RDI3(IY)
                RD2     = MIN ( 1. , RDI2(IY) * CG(IK,ISEA) )
                RD4     = MIN ( 1. , RDI4(IY) * CG(IK,ISEA) )
                VQ(IXY          ) = VQ(IXY          )                   &
                                     + AQ(IXY) * (3.-RD2-RD4)/3.
  !
                IXC    = SIGN(NY,IIXY1(IY))
                IYC    = IIXY1(IY) - IXC
                IF ( MAPSTX(IXY+IIXY1(IY)) .GE. 1 .AND.                 &
                     .NOT. ( MAPSTX(IXY+IXC).LE.0 .AND.                 &
                             MAPSTX(IXY+IYC).LE.0 ) ) THEN
                         VQ(IXY+IIXY1(IY)) = VQ(IXY+IIXY1(IY))          &
                                              + AQ(IXY) * RD2*RD1/6.
                  ELSE
                         VQ(IXY          ) = VQ(IXY          )          &
                                              + AQ(IXY) * RD2*RD1/6.
                  END IF
                IF ( MAPSTX(IXY-IIXY1(IY)) .GE. 1 .AND.                 &
                     .NOT. ( MAPSTX(IXY-IXC).LE.0 .AND.                 &
                             MAPSTX(IXY-IYC).LE.0 ) ) THEN
                         VQ(IXY-IIXY1(IY)) = VQ(IXY-IIXY1(IY))          &
                                              + AQ(IXY) * RD2*RD1/6.
                  ELSE
                         VQ(IXY          ) = VQ(IXY          )          &
                                              + AQ(IXY) * RD2*RD1/6.
                  END IF
   
                IF ( MAPSTX(IXY+IIXY2(IY)) .GE. 1 ) THEN
                         VQ(IXY+IIXY2(IY)) = VQ(IXY+IIXY2(IY))          &
                                              + AQ(IXY) * (1.-RD1)*RD2/6.
                  ELSE
                         VQ(IXY          ) = VQ(IXY          )          &
                                              + AQ(IXY) * (1.-RD1)*RD2/6.
                  END IF
                IF ( MAPSTX(IXY-IIXY2(IY)) .GE. 1 ) THEN
                         VQ(IXY-IIXY2(IY)) = VQ(IXY-IIXY2(IY))          &
                                              + AQ(IXY) * (1.-RD1)*RD2/6.
                  ELSE
                         VQ(IXY          ) = VQ(IXY          )          &
                                              + AQ(IXY) * (1.-RD1)*RD2/6.
                  END IF
  !
                IXC    = SIGN(NY,IIXY3(IY))
                IYC    = IIXY3(IY) - IXC
                IF ( MAPSTX(IXY+IIXY3(IY)) .GE. 1 .AND.                 &
                     .NOT. ( MAPSTX(IXY+IXC).LE.0 .AND.                 &
                             MAPSTX(IXY+IYC).LE.0 ) ) THEN
                         VQ(IXY+IIXY3(IY)) = VQ(IXY+IIXY3(IY))          &
                                              + AQ(IXY) * RD4*RD3/6.
                  ELSE
                         VQ(IXY          ) = VQ(IXY          )          &
                                              + AQ(IXY) * RD4*RD3/6.
                  END IF
                IF ( MAPSTX(IXY-IIXY3(IY)) .GE. 1 .AND.                 &
                     .NOT. ( MAPSTX(IXY-IXC).LE.0 .AND.                 &
                             MAPSTX(IXY-IYC).LE.0 ) ) THEN
                         VQ(IXY-IIXY3(IY)) = VQ(IXY-IIXY3(IY))          &
                                              + AQ(IXY) * RD4*RD3/6.
                  ELSE
                         VQ(IXY          ) = VQ(IXY          )          &
                                              + AQ(IXY) * RD4*RD3/6.
                  END IF
  !
                IF ( MAPSTX(IXY+IIXY4(IY)) .GE. 1 ) THEN
                         VQ(IXY+IIXY4(IY)) = VQ(IXY+IIXY4(IY))          &
                                              + AQ(IXY) * (1.-RD3)*RD4/6.
                  ELSE
                         VQ(IXY          ) = VQ(IXY          )          &
                                              + AQ(IXY) * (1.-RD3)*RD4/6.
                  END IF
                IF ( MAPSTX(IXY-IIXY4(IY)) .GE. 1 ) THEN
                         VQ(IXY-IIXY4(IY)) = VQ(IXY-IIXY4(IY))          &
                                              + 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
  !
          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
  !
          IF ( YFIRST ) THEN
              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 )
            ELSE
              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 )
            END IF
  !
  ! 3.c 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
                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 --------------------------- *
  !
        DO ISEA=1, NSEA
          IXY    = MAPSF(ISEA,3)
          IF ( MAPSTA(IXY) .GT. 0 ) THEN
              VQ(IXY) =  MAX ( 0. , CG(IK,ISEA)/CLATS(ISEA)*VQ(IXY) )
            END IF
          END DO
  !
        RETURN
  !
  ! Formats
  !
  !/
  !/ End of W3XYP3 ----------------------------------------------------- /
  !/
        END SUBROUTINE W3XYP3
  !/ ------------------------------------------------------------------- /
        SUBROUTINE W3KTP3 ( ISEA, FACTH, FACK, CTHG0, CG, WN, DEPTH,    &
                            DDDX, DDDY, CX, CY, DCXDX, DCXDY,           &
                            DCYDX, DCYDY, VA )
  !/
  !/    *** THIS ROUTINE SHOULD BE IDENTICAL TO W3KTP2 ***
  !/
  !/                  +-----------------------------------+
  !/                  | WAVEWATCH III           NOAA/NCEP |
  !/                  |           H. L. Tolman            |
  !/                  |                        FORTRAN 90 |
  !/                  | Last update :         17-Dec-2004 |
  !/                  +-----------------------------------+
  !/
  !/    14-Feb-2000 : Origination.                        ( version 2.08 )
  !/    17-Dec-2004 : Multiple grid version.              ( version 3.06 )
  !/
  !  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.
  !       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 W3CONSTANTS
        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
        USE W3UQCKMD
  !/
        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(INOUT)     :: VA(NSPEC)
  !/
  !/ ------------------------------------------------------------------- /
  !/ Local parameters
  !/
        INTEGER                 :: ITH, IK, ISP
        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)
  !/
  !/ ------------------------------------------------------------------- /
  !/
  !
  ! 1.  Preparations --------------------------------------------------- *
  ! 1.a Initialize arrays
  !
        IF ( FLCK ) VQ = 0.
  !
  ! 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
  !
  ! 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
  !
  ! 3.c 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
                CALL W3QCK2 ( NTH, NK2, NTH, NK2, CFLK, FACK, DB, DM,   &
                              VQ, .FALSE., 1, MAPTH2, NSPEC,            &
                              MAPWN2, NSPEC-NTH, NSPEC, NSPEC+NTH,      &
                              NDSE, NDST )
              END IF
            IF ( FLCTH ) THEN
                CALL W3QCK1 ( NTH, NK2, NTH, NK2, VCFLT, VQ, .TRUE.,    &
                              NK2, MAPTH2, NSPEC, MAPTH2, NSPEC, NSPEC, &
                              NSPEC, NDSE, NDST )
              END IF
          ELSE
            IF ( FLCTH ) THEN
                CALL W3QCK1 ( NTH, NK2, NTH, NK2, VCFLT, VQ, .TRUE.,    &
                              NK2, MAPTH2, NSPEC, MAPTH2, NSPEC, NSPEC, &
                              NSPEC, NDSE, NDST )
              END IF
            IF ( FLCK )  THEN
                DO ITH=1, NTH
                  VQ(NK+2+(ITH-1)*NK2) = FACHFA * VQ(NK+1+(ITH-1)*NK2)
                  END DO
                CALL W3QCK2 ( NTH, NK2, NTH, NK2, CFLK, FACK, DB, DM,   &
                              VQ, .FALSE., 1, MAPTH2, NSPEC,            &
                              MAPWN2, NSPEC-NTH, NSPEC, NSPEC+NTH,      &
                              NDSE, NDST )
            END IF
          END IF
  !
  ! 6.  Store reults --------------------------------------------------- *
  !
        DO ISP=1, NSPEC
          VA(ISP) = VQ(MAPTH2(ISP))
          END DO
  !
        RETURN
  !
  ! Formats
  !
  !/
  !/ End of W3KTP3 ----------------------------------------------------- /
  !/
        END SUBROUTINE W3KTP3
  !/
  !/ End of module W3PRO3MD -------------------------------------------- /
  !/
        END MODULE W3PRO3MD
  