w3pro3md Module Reference

Bundles routines for third order propagation scheme in single module. More...

Functions/Subroutines
subroutine	w3map3
	Generate 'map' arrays for the ULTIMATE QUICKEST scheme. More...
subroutine	w3mapt
	Generate 'map' arrays for the ULTIMATE QUICKEST scheme to combine GSE alleviation with obstructions. More...
subroutine	w3xyp3 (ISP, DTG, MAPSTA, MAPFS, VQ, VGX, VGY)
	Propagation in phyiscal space for a given spectral component. More...
subroutine	w3ktp3 (ISEA, FACTH, FACK, CTHG0, CG, WN, DW, DDDX, DDDY, CX, CY, DCXDX, DCXDY, DCYDX, DCYDY, DCDX, DCDY, VA, CFLTHMAX, CFLKMAX)
	Propagation in spectral space. More...
subroutine	w3cflxy (ISEA, DTG, MAPSTA, MAPFS, CFLXYMAX, VGX, VGY)
	Computes the maximum CFL number for spatial advection. More...

Detailed Description

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.

Function/Subroutine Documentation

w3cflxy()

subroutine w3pro3md::w3cflxy	(	integer, intent(in)	ISEA,
		real, intent(in)	DTG,
		integer, dimension(ny*nx), intent(in)	MAPSTA,
		integer, dimension(ny*nx), intent(in)	MAPFS,
		real, intent(inout)	CFLXYMAX,
		real, intent(in)	VGX,
		real, intent(in)	VGY
	)

Computes the maximum CFL number for spatial advection.

Used for diagnostic purposes (Could be used to define a local time step ...).

Parameters

[in]	ISEA	Index of grid point.
[in]	DTG	Total time step.
[in]	MAPSTA	Grid point status map.
[in]	MAPFS	Storage map.
[in,out]	CFLXYMAX	Maximum CFL number for XY propagation.
[in]	VGX	Speed of grid.
[in]	VGY	Speed of grid.

Author

F. Ardhuin

Date

31-Oct-2010

Definition at line 1971 of file w3pro3md.F90.

    !/
    !/                  +-----------------------------------+
    !/                  | 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 ----------------------------------------------------- /
    !/

References w3adatmd::atrnx, w3adatmd::atrny, w3odatmd::bbpi0, w3odatmd::bbpin, w3adatmd::cg, w3gdatmd::clats, w3adatmd::cx, w3adatmd::cy, w3gdatmd::dpdx, w3gdatmd::dpdy, w3gdatmd::dqdx, w3gdatmd::dqdy, w3timemd::dsec21(), w3gdatmd::dtcfl, w3gdatmd::dth, w3gdatmd::ecos, w3gdatmd::esin, w3gdatmd::flagll, w3odatmd::flbpi, w3idatmd::flcur, w3gdatmd::flcx, w3gdatmd::flcy, w3gdatmd::gsqrt, w3odatmd::isbpi, w3adatmd::itime, w3adatmd::mapaxy, w3adatmd::mapcxy, w3gdatmd::mapsf, w3adatmd::maptrn, w3adatmd::mapx2, w3adatmd::mapy2, w3adatmd::nact, w3odatmd::nbi, w3adatmd::ncent, w3odatmd::ndse, w3odatmd::ndst, w3gdatmd::nk, w3adatmd::nmx0, w3adatmd::nmx1, w3adatmd::nmx2, w3adatmd::nmy0, w3adatmd::nmy1, w3adatmd::nmy2, w3gdatmd::nsea, w3gdatmd::nth, w3gdatmd::nx, w3gdatmd::ny, w3gdatmd::pfmove, w3gdatmd::sig, w3servmd::strace(), w3odatmd::tbpi0, w3odatmd::tbpin, w3wdatmd::time, w3gdatmd::wdcg, w3gdatmd::wdth, and w3gdatmd::xfr.

Referenced by w3wavemd::w3wave().

w3ktp3()

subroutine w3pro3md::w3ktp3	(	integer, intent(in)	ISEA,
		real, intent(in)	FACTH,
		real, intent(in)	FACK,
		real, intent(in)	CTHG0,
		real, dimension(0:nk+1), intent(in)	CG,
		real, dimension(0:nk+1), intent(in)	WN,
		real, intent(in)	DW,
		real, intent(in)	DDDX,
		real, intent(in)	DDDY,
		real, intent(in)	CX,
		real, intent(in)	CY,
		real, intent(in)	DCXDX,
		real, intent(in)	DCXDY,
		real, intent(in)	DCYDX,
		real, intent(in)	DCYDY,
		real, dimension(0:nk+1), intent(in)	DCDX,
		real, dimension(0:nk+1), intent(in)	DCDY,
		real, dimension(nspec), intent(inout)	VA,
		real, intent(out)	CFLTHMAX,
		real, intent(out)	CFLKMAX
	)

Propagation in spectral space.

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.

Parameters

[in]	ISEA	Number of sea point.
[in]	FACTH	Factor in propagation velocity.
[in]	FACK	Factor in propagation velocity.
[in]	CTHG0	Factor in great circle refracftion term.
[in]	CG	Local group velocities.
[in]	WN	Local wavenumbers.
[in]	DW	Depth.
[in]	DDDX	Depth gradients.
[in]	DDDY	Depth gradients.
[in]	CX	Current components.
[in]	CY	Current components.
[in]	DCXDX	Current gradients.
[in]	DCXDY	Current gradients.
[in]	DCYDX	Current gradients.
[in]	DCYDY	Current gradients.
[in]	DCDX	Phase speed gradients.
[in]	DCDY	Phase speed gradients.
[in,out]	VA	Spectrum.
[out]	CFLTHMAX
[out]	CFLKMAX

Author

H. L. Tolman

Date

01-Jul-2013

Definition at line 1512 of file w3pro3md.F90.

    !/
    !/    *** 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 ----------------------------------------------------- /
    !/

References w3gdatmd::ctmax, w3gdatmd::dmin, w3gdatmd::dsip, w3gdatmd::ec2, w3gdatmd::ecos, w3gdatmd::es2, w3gdatmd::esc, w3gdatmd::esin, w3gdatmd::fachfa, w3gdatmd::flck, w3gdatmd::flcth, w3idatmd::flcur, w3adatmd::itime, w3adatmd::mapth2, w3gdatmd::mapwn, w3adatmd::mapwn2, w3odatmd::ndse, w3odatmd::ndst, w3gdatmd::nk, w3gdatmd::nk2, w3gdatmd::nspec, w3gdatmd::nth, w3gdatmd::sig, w3servmd::strace(), constants::tpi, w3uqckmd::w3qck1(), w3uqckmd::w3qck2(), w3uno2md::w3uno2(), and w3uno2md::w3uno2r().

Referenced by w3wavemd::w3wave().

w3map3()

subroutine w3pro3md::w3map3

Generate 'map' arrays for the ULTIMATE QUICKEST scheme.

Author

H. L. Tolman

Date

01-Apr-2008

Definition at line 140 of file w3pro3md.F90.

    !/
    !/                  +-----------------------------------+
    !/                  | 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 ----------------------------------------------------- /
    !/

References w3gdatmd::gtype, w3adatmd::mapaxy, w3adatmd::mapcxy, w3gdatmd::mapsf, w3gdatmd::mapsta, w3adatmd::mapth2, w3adatmd::mapwn2, w3adatmd::mapx2, w3adatmd::mapy2, w3adatmd::nact, w3adatmd::ncent, w3odatmd::ndst, w3gdatmd::nk, w3adatmd::nmx0, w3adatmd::nmx1, w3adatmd::nmx2, w3adatmd::nmy0, w3adatmd::nmy1, w3adatmd::nmy2, w3gdatmd::nsea, w3gdatmd::nspec, w3gdatmd::nth, w3gdatmd::nx, w3gdatmd::ny, and w3servmd::strace().

Referenced by w3wavemd::w3wave().

w3mapt()

subroutine w3pro3md::w3mapt

Generate 'map' arrays for the ULTIMATE QUICKEST scheme to combine GSE alleviation with obstructions.

Author

H. L. Tolman

Date

17-Dec-2004

Definition at line 525 of file w3pro3md.F90.

    !/
    !/                  +-----------------------------------+
    !/                  | 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 ----------------------------------------------------- /
    !/

References w3adatmd::atrnx, w3adatmd::atrny, w3gdatmd::mapsf, w3adatmd::maptrn, w3gdatmd::nsea, w3gdatmd::nx, w3gdatmd::ny, and w3servmd::strace().

Referenced by w3wavemd::w3wave().

w3xyp3()

subroutine w3pro3md::w3xyp3	(	integer, intent(in)	ISP,
		real, intent(in)	DTG,
		integer, dimension(ny*nx), intent(in)	MAPSTA,
		integer, dimension(ny*nx), intent(in)	MAPFS,
		real, dimension(1-ny:ny*(nx+2)), intent(inout)	VQ,
		real, intent(in)	VGX,
		real, intent(in)	VGY
	)

Propagation in phyiscal space for a given spectral component.

Parameters

[in]	ISP	Number of spectral bin (IK-1)*NTH+ITH.
[in]	DTG	Total time step.
[in]	MAPSTA	Grid point status map.
[in]	MAPFS	Storage map.
[in,out]	VQ	Field to propagate.
[in]	VGX	Speed of grid.
[in]	VGY	Speed of grid.

Author

H. L. Tolman

Date

27-May-2014

Definition at line 639 of file w3pro3md.F90.

    !/
    !/                  +-----------------------------------+
    !/                  | 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 ----------------------------------------------------- /
    !/

References w3adatmd::atrnx, w3adatmd::atrny, w3odatmd::bbpi0, w3odatmd::bbpin, w3adatmd::cg, w3gdatmd::clats, w3adatmd::cx, w3adatmd::cy, constants::dera, w3gdatmd::dpdx, w3gdatmd::dpdy, w3gdatmd::dqdx, w3gdatmd::dqdy, w3timemd::dsec21(), w3gdatmd::dtcfl, w3gdatmd::dth, w3gdatmd::ecos, w3gdatmd::esin, w3servmd::extcde(), w3gdatmd::flagll, w3odatmd::flbpi, w3idatmd::flcur, w3gdatmd::flcx, w3gdatmd::flcy, constants::grav, w3gdatmd::gsqrt, w3odatmd::iaproc, w3gdatmd::iclose, w3gdatmd::iclose_none, w3gdatmd::iclose_smpl, w3gdatmd::iclose_trpl, w3odatmd::isbpi, w3adatmd::itime, w3adatmd::mapaxy, w3adatmd::mapcxy, w3gdatmd::mapsf, w3adatmd::maptrn, w3adatmd::mapx2, w3adatmd::mapy2, w3adatmd::nact, w3odatmd::naperr, w3odatmd::nbi, w3adatmd::ncent, w3odatmd::ndse, w3odatmd::ndst, w3gdatmd::nk, w3adatmd::nmx0, w3adatmd::nmx1, w3adatmd::nmx2, w3adatmd::nmy0, w3adatmd::nmy1, w3adatmd::nmy2, w3gdatmd::nsea, w3gdatmd::nth, w3gdatmd::nx, w3gdatmd::ny, w3gdatmd::pfmove, constants::radius, w3gdatmd::sig, w3servmd::strace(), w3odatmd::tbpi0, w3odatmd::tbpin, w3wdatmd::time, w3uqckmd::w3qck3(), w3uno2md::w3uno2s(), w3gdatmd::wdcg, w3gdatmd::wdth, and w3gdatmd::xfr.

Referenced by w3wavemd::w3wave().