w3pro2md Module Reference

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

Functions/Subroutines
subroutine	w3map2
	Generate 'map' arrays for the ULTIMATE QUICKEST scheme. More...
subroutine	w3xyp2 (ISP, DTG, MAPSTA, MAPFS, VQ, VGX, VGY)
	Propagation in physical space for a given spectral component. More...
subroutine	w3ktp2 (ISEA, FACTH, FACK, CTHG0, CG, WN, DEPTH, DDDX, DDDY, CX, CY, DCXDX, DCXDY, DCYDX, DCYDY, DCDX, DCDY, VA)
	Propagation in spectral space. More...

Detailed Description

Bundles routines for third order porpagation scheme in single module.

Author

H. L. Tolman

Date

29-May-2014

Copyright

Copyright 2009-2022 National Weather Service (NWS), National Oceanic and Atmospheric Administration. All rights reserved. WAVEWATCH III is a trademark of the NWS. No unauthorized use without permission.

Function/Subroutine Documentation

w3ktp2()

subroutine w3pro2md::w3ktp2	(	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)	DEPTH,
		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
	)

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]	DEPTH	Depth.
[in]	DDDX	Depth gradient.
[in]	DDDY	Depth gradient.
[in]	CX	Current component.
[in]	CY	Current component.
[in]	DCXDX	Current gradients.
[in]	DCXDY	Current gradients.
[in]	DCYDX	Current gradients.
[in]	DCYDY	Current gradients.
[in]	DCDX	Phase speed gradient.
[in]	DCDY	Phase speed gradient.
[in,out]	VA	Spectrum.

Author

H. L. Tolman

Date

01-Jul-2013

Definition at line 1305 of file w3pro2md.F90.

    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III           NOAA/NCEP |
    !/                  |           H. L. Tolman            |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         01-Jul-2013 |
    !/                  +-----------------------------------+
    !/
    !/    14-Feb-2000 : Origination.                        ( version 2.08 )
    !/    17-Dec-2004 : Multiple grid version.              ( version 3.06 )
    !/    01-Jul-2013 : Adding UQ and UNO switches to chose between third
    !/                  and second order schemes.           ( version 4.12 )
    !/
    !  1. Purpose :
    !
    !     Propagation in spectral space.
    !
    !  2. Method :
    !
    !     Third order QUICKEST scheme with ULTIMATE limiter.
    !
    !     As with the spatial propagation, the two spaces are considered
    !     independently, but the propagation is performed in a 2-D space.
    !     Compared to the propagation in physical space, the directions
    !     rerpesent a closed space and are therefore comparable to the
    !     longitudinal or 'X' propagation. The wavenumber space has to be
    !     extended to allow for boundary treatment. Using a simple first
    !     order boundary treatment at both sided, two points need to
    !     be added. This implies that the spectrum needs to be extended,
    !     shifted and rotated, as is performed using MAPTH2 as set
    !     in W3MAP3.
    !
    !  3. Parameters :
    !
    !     Parameter list
    !     ----------------------------------------------------------------
    !       ISEA    Int.   I   Number of sea point.
    !       FACTH/K Real   I   Factor in propagation velocity.
    !       CTHG0   Real   I   Factor in great circle refracftion term.
    !       MAPxx2  I.A.   I   Propagation and storage maps.
    !       CG      R.A.   I   Local group velocities.
    !       WN      R.A.   I   Local wavenumbers.
    !       DEPTH   R.A.   I   Depth.
    !       DDDx    Real   I   Depth gradients.
    !       CX/Y    Real   I   Current components.
    !       DCxDx   Real   I   Current gradients.
    !       DCDX-Y  Real   I   Phase speed gradients.
    !       VA      R.A.  I/O  Spectrum.
    !     ----------------------------------------------------------------
    !
    !     Local variables.
    !     ----------------------------------------------------------------
    !       DSDD    R.A.  Partial derivative of sigma for depth.
    !       FDD, FDU, FDG, FCD, FCU
    !               R.A.  Directionally varying part of depth, current and
    !                     great-circle refraction terms and of consit.
    !                     of Ck term.
    !       CFLT-K  R.A.  Propagation velocities of local fluxes.
    !       DB      R.A.  Wavenumber band widths at cell centers.
    !       DM      R.A.  Wavenumber band widths between cell centers and
    !                     next cell center.
    !       Q       R.A.  Extracted spectrum
    !     ----------------------------------------------------------------
    !
    !  4. Subroutines used :
    !
    !       W3QCK1   Actual propagation routine.
    !       W3QCK2   Actual propagation routine.
    !       STRACE   Service routine.
    !
    !  5. Called by :
    !
    !       W3WAVE   Wave model routine.
    !
    !  6. Error messages :
    !
    !       None.
    !
    !  8. Structure :
    !
    !     -----------------------------------------------------------------
    !       1.  Preparations
    !         a Initialize arrays
    !         b Set constants and counters
    !       2.  Point  preparations
    !         a Calculate DSDD
    !         b Extract spectrum
    !       3.  Refraction velocities
    !         a Filter level depth reffraction.
    !         b Depth refratcion velocity.
    !         c Current refraction velocity.
    !       4.  Wavenumber shift velocities
    !         a Prepare directional arrays
    !         b Calcuate velocity.
    !       5.  Propagate.
    !       6.  Store results.
    !     -----------------------------------------------------------------
    !
    !  9. Switches :
    !
    !       !/S     Enable subroutine tracing.
    !       !/T     Enable general test output.
    !
    ! 10. Source code :
    !
    !/ ------------------------------------------------------------------- /
    USE constants
    USE w3gdatmd, ONLY: nk, nk2, nth, nspec, sig, dsip, ecos, esin, &
         ec2, esc, es2, fachfa, mapwn, flcth, flck,  &
         ctmax
    USE w3adatmd, ONLY: mapth2, mapwn2, itime
    USE w3idatmd, ONLY: flcur
    USE w3odatmd, ONLY: ndse, ndst
#ifdef W3_S
    USE w3servmd, ONLY: strace
#endif
#ifdef W3_UQ
    USE w3uqckmd
#endif
#ifdef W3_UNO
    USE w3uno2md
#endif
    !/
    IMPLICIT NONE
    !/
    !/ ------------------------------------------------------------------- /
    !/ Parameter list
    !/
    INTEGER, INTENT(IN)     :: ISEA
    REAL, INTENT(IN)        :: FACTH, FACK, CTHG0, CG(0:NK+1),      &
         WN(0:NK+1), DEPTH, DDDX, DDDY,       &
         CX, CY, DCXDX, DCXDY, DCYDX, DCYDY
    REAL, INTENT(IN)        :: DCDX(0:NK+1), DCDY(0:NK+1)
    REAL, INTENT(INOUT)     :: VA(NSPEC)
    !/
    !/ ------------------------------------------------------------------- /
    !/ Local parameters
    !/
    INTEGER                 :: ITH, IK, ISP
#ifdef W3_S
    INTEGER, SAVE           :: IENT = 0
#endif
    REAL                    :: FDDMAX, FDG, FKD, FKD0, DCYX,        &
         DCXXYY, DCXY, DCXX, DCXYYX, DCYY
    REAL                    :: DSDD(0:NK+1), FRK(NK), FRG(NK),      &
         FKC(NTH), VQ(-NK-1:NK2*(NTH+2)),     &
         DB(NK2,NTH+1), DM(NK2,0:NTH+1),      &
         VCFLT(NK2*(NTH+1)), CFLK(NK2,NTH)
    !/
    !/ ------------------------------------------------------------------- /
    !/
#ifdef W3_S
    CALL strace (ient, 'W3KTP2')
#endif
    !
    ! 1.  Preparations --------------------------------------------------- *
    ! 1.a Initialize arrays
    !
    IF ( flck ) vq = 0.
    !
#ifdef W3_T
    WRITE (ndst,9000) flcth, flck, facth, fack, ctmax
    WRITE (ndst,9010) isea, depth, cx, cy, dddx, dddy,           &
         dcxdx, dcxdy, dcydx, dcydy
#endif
    !
    ! 2.  Preparation for point ------------------------------------------ *
    ! 2.a Array with partial derivative of sigma versus depth
    !
    DO ik=0, nk+1
      IF ( depth*wn(ik) .LT. 5. ) THEN
        dsdd(ik) = max( 0. ,                                     &
             cg(ik)*wn(ik)-0.5*sig(ik) ) / depth
      ELSE
        dsdd(ik) = 0.
      END IF
    END DO
    !
#ifdef W3_T
    WRITE (ndst,9020)
    DO ik=1, nk+1
      WRITE (ndst,9021) ik, tpi/sig(ik), tpi/wn(ik),             &
           cg(ik), dsdd(ik)
    END DO
#endif
    !
    ! 2.b Extract spectrum
    !
    DO isp=1, nspec
      vq(mapth2(isp)) = va(isp)
    END DO
    !
    ! 3.  Refraction velocities ------------------------------------------ *
    !
    IF ( flcth ) THEN
      !
      ! 3.a Set slope filter for depth refraction
      !
      fddmax = 0.
      fdg    = facth * cthg0
      !
      DO ith=1, nth/2
        fddmax = max(fddmax,abs(esin(ith)*dddx-ecos(ith)*dddy))
      END DO
      !
      DO ik=1, nk
        frk(ik) = facth * dsdd(ik) / wn(ik)
        frk(ik) = frk(ik) / max( 1. , frk(ik)*fddmax/ctmax )
        frg(ik) = fdg * cg(ik)
      END DO
      !
      ! 3.b Depth refraction and great-circle propagation
      !
      DO isp=1, nspec
        vcflt(mapth2(isp)) = frg(mapwn(isp)) * ecos(isp)          &
             + frk(mapwn(isp)) * ( esin(isp)*dddx - ecos(isp)*dddy )
      END DO
      !
#ifdef W3_REFRX
      ! 3.c @C/@x refraction and great-circle propagation
      vcflt  = 0.
      frk    = 0.
      fddmax = 0.
      !
      DO isp=1, nspec
        fddmax = max( fddmax , abs(                      &
             esin(isp)*dcdx(mapwn(isp)) - ecos(isp)*dcdy(mapwn(isp)) ) )
      END DO
      !
      DO ik=1, nk
        frk(ik) = facth * cg(ik) * wn(ik) / sig(ik)
      END DO
      DO isp=1, nspec
        vcflt(mapth2(isp)) = frg(mapwn(isp)) * ecos(isp)   &
             + frk(mapwn(isp)) * ( esin(isp)*dcdx(mapwn(isp)) &
             - ecos(isp)*dcdy(mapwn(isp)) )
      END DO
#endif
      !
      ! 3.d Current refraction
      !
      IF ( flcur ) THEN
        !
        dcyx   = facth *   dcydx
        dcxxyy = facth * ( dcxdx - dcydy )
        dcxy   = facth *   dcxdy
        !
        DO isp=1, nspec
          vcflt(mapth2(isp)) = vcflt(mapth2(isp)) +             &
               es2(isp)*dcyx  + esc(isp)*dcxxyy - ec2(isp)*dcxy
        END DO
        !
      END IF
      !
    END IF
    !
    ! 4.  Wavenumber shift velocities ------------------------------------ *
    !     FACK is just the time step, which is accounted for in W3QCK2
    !
    IF ( flck ) THEN
      !
      ! 4.a Directionally dependent part
      !
      dcxx   =  -   dcxdx
      dcxyyx =  - ( dcxdy + dcydx )
      dcyy   =  -   dcydy
      fkd    =    ( cx*dddx + cy*dddy )
      !
      DO ith=1, nth
        fkc(ith) = ec2(ith)*dcxx +                                &
             esc(ith)*dcxyyx + es2(ith)*dcyy
      END DO
      !
      ! 4.b Velocities
      !
      DO ik=0, nk+1
        fkd0   = fkd / cg(ik) * dsdd(ik)
        DO ith=1, nth
          cflk(ik+1,ith) = fkd0 + wn(ik)*fkc(ith)
        END DO
      END DO
      !
      ! 4.c Band widths
      !
      DO ik=0, nk
        db(ik+1,1) = dsip(ik) / cg(ik)
        dm(ik+1,1) = wn(ik+1) - wn(ik)
      END DO
      db(nk+2,1) = dsip(nk+1) / cg(nk+1)
      dm(nk+2,1) = 0.
      !
      DO ith=2, nth
        DO ik=1, nk+2
          db(ik,ith) = db(ik,1)
          dm(ik,ith) = dm(ik,1)
        END DO
      END DO
      !
    END IF
    !
    ! 5.  Propagate ------------------------------------------------------ *
    !
    IF ( mod(itime,2) .EQ. 0 ) THEN
      IF ( flck ) THEN
        DO ith=1, nth
          vq(nk+2+(ith-1)*nk2) = fachfa * vq(nk+1+(ith-1)*nk2)
        END DO
        !
#ifdef W3_UQ
        CALL w3qck2 ( nth, nk2, nth, nk2, cflk, fack, db, dm,   &
             vq, .false., 1, mapth2, nspec,            &
             mapwn2, nspec-nth, nspec, nspec+nth,      &
             ndse, ndst )
#endif
        !
#ifdef W3_UNO
        CALL w3uno2 ( nth, nk2, nth, nk2, cflk, fack, db, dm,  &
             vq, .false., 1, mapth2, nspec,           &
             mapwn2, nspec-nth, nspec, nspec+nth,     &
             ndse, ndst )
#endif
      END IF
      IF ( flcth ) THEN
        !
#ifdef W3_UQ
        CALL w3qck1 ( nth, nk2, nth, nk2, vcflt, vq, .true.,    &
             nk2, mapth2, nspec, mapth2, nspec, nspec, &
             nspec, ndse, ndst )
#endif
        !
#ifdef W3_UNO
        CALL w3uno2r( nth, nk2, nth, nk2, vcflt, vq, .true.,   &
             nk2, mapth2, nspec, mapth2, nspec, nspec,&
             nspec, ndse, ndst )
#endif
        !
      END IF
    ELSE
      IF ( flcth ) THEN
        !
#ifdef W3_UQ
        CALL w3qck1 ( nth, nk2, nth, nk2, vcflt, vq, .true.,    &
             nk2, mapth2, nspec, mapth2, nspec, nspec, &
             nspec, ndse, ndst )
#endif
        !
#ifdef W3_UNO
        CALL w3uno2r( nth, nk2, nth, nk2, vcflt, vq, .true.,   &
             nk2, mapth2, nspec, mapth2, nspec, nspec,&
             nspec, ndse, ndst )
#endif
        !
      END IF
      IF ( flck )  THEN
        DO ith=1, nth
          vq(nk+2+(ith-1)*nk2) = fachfa * vq(nk+1+(ith-1)*nk2)
        END DO
        !
#ifdef W3_UQ
        CALL w3qck2 ( nth, nk2, nth, nk2, cflk, fack, db, dm,   &
             vq, .false., 1, mapth2, nspec,            &
             mapwn2, nspec-nth, nspec, nspec+nth,      &
             ndse, ndst )
#endif
        !
#ifdef W3_UNO
        CALL w3uno2 ( nth, nk2, nth, nk2, cflk, fack, db, dm,  &
             vq, .false., 1, mapth2, nspec,           &
             mapwn2, nspec-nth, nspec, nspec+nth,     &
             ndse, ndst )
#endif
        !
      END IF
    END IF
    !
    ! 6.  Store reults --------------------------------------------------- *
    !
    DO isp=1, nspec
      va(isp) = vq(mapth2(isp))
    END DO
    !
    RETURN
    !
    ! Formats
    !
#ifdef W3_T
9000 FORMAT ( ' TEST W3KTP2 : FLCTH-K, FACTH-K, CTMAX  :',        &
         2l2,2e10.3,f7.3)
9010 FORMAT ( ' TEST W3KTP2 : LOCAL DATA :',i7,f7.1,2f6.2,1x,6e10.2)
9020 FORMAT ( ' TEST W3KTP2 : IK, T, L, CG, DSDD : ')
9021 FORMAT ( '              ',i3,f7.2,f7.1,f7.2,e11.3)
#endif
    !
#ifdef W3_T0
9040 FORMAT (/' TEST W3KTP2 : NORMALIZED ',a/)
9041 FORMAT (1x,60(1x,i2))
9042 FORMAT (1x,60i3)
#endif
    !/
    !/ End of W3KTP2 ----------------------------------------------------- /
    !/

References w3gdatmd::ctmax, 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().

w3map2()

subroutine w3pro2md::w3map2

Generate 'map' arrays for the ULTIMATE QUICKEST scheme.

Author

H. L. Tolman

Date

09-Nov-2005

Definition at line 136 of file w3pro2md.F90.

    !/
    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III           NOAA/NCEP |
    !/                  |           H. L. Tolman            |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         09-Nov-2005 |
    !/                  +-----------------------------------+
    !/
    !/    19-Dec-1996 : Final FORTRAN 77                    ( version 1.18 )
    !/    15-Dec-1999 : Upgrade to FORTRAN 90               ( version 2.00 )
    !/    09-Feb-2001 : Clean up of parameter lists         ( version 2.08 )
    !/    20-Dec-2004 : Multiple grid version.              ( version 3.06 )
    !/    09-Nov-2005 : Removing soft boundary option.      ( version 3.08 )
    !/
    !  1. Purpose :
    !
    !     Generate 'map' arrays for the ULTIMATE QUICKEST scheme.
    !
    !  2. Method :
    !
    !     MAPX2, MAPY2, MAPTH2 and MAPWN2 contain consecutive 1-D spatial
    !     grid counters (e.g., IXY = (IX-1)*MY + IY). The arrays are
    !     devided in three parts. For MAPX2, these ranges are :
    !
    !         1    - NMX0  Counters IXY for which grid point (IX,IY) and
    !                      (IX+1,IY) both are active grid points.
    !       NMX0+1 - NMX1  Id. only (IX,IY) active.
    !       NMX1+1 - NMX2  Id. only (IX+1,IY) active.
    !
    !     The array MAPY2 has a similar layout varying IY instead of IX.
    !     MAPXY contains similar information for the cross term in the
    !     diffusion correction (counter NMXY only).
    !
    !  3. Parameters :
    !
    !     Parameter list
    !     ----------------------------------------------------------------
    !     ----------------------------------------------------------------
    !
    !  4. Subroutines used :
    !
    !     See module documentation.
    !
    !  5. Called by :
    !
    !      Name      Type  Module   Description
    !     ----------------------------------------------------------------
    !      W3WAVE    Subr. W3WAVEMD Wave model routine.
    !     ---------------------------------------------------------------
    !
    !  6. Error messages :
    !
    !  7. Remarks :
    !
    !  8. Structure :
    !
    !     ------------------------------------------------------
    !      1.   Map MAPX2
    !        a  Range 1 to NMX0
    !        b  Range NMX0+1 to NMX1
    !        c  Range NMX1+1 to NMX2
    !      2.   Map MAPY2
    !        a  Range 1 to NMY0
    !        b  Range NMY0+1 to NMY1
    !        c  Range NMY1+1 to NMY2
    !      3.   Map MAPAXY and MAPXY
    !      4.   Maps for intra-spectral propagation
    !        a  MAPTH2, MAPATK
    !        b  MAPWN2
    !     ------------------------------------------------------
    !
    !  9. Switches :
    !
    !     !/S   Enable subroutine tracing.
    !     !/T   Enable test output.
    !
    ! 10. Source code :
    !/ ------------------------------------------------------------------- /
    USE w3gdatmd, ONLY: nk, nth, nspec, nx, ny, nsea, mapsta
    USE w3adatmd, ONLY: nmx0, nmx1, nmx2, nmy0, nmy1, nmy2, nact,   &
         nmxy, mapx2, mapy2, mapaxy, mapxy,          &
         mapth2, mapwn2
    USE w3odatmd, ONLY: ndst
#ifdef W3_S
    USE w3servmd, ONLY: strace
#endif
    !/
    IMPLICIT NONE
    !/
    !/ ------------------------------------------------------------------- /
    !/ Parameter list
    !/
    !/ ------------------------------------------------------------------- /
    !/ Local parameters
    !/
    INTEGER                 :: IX, IY, IXY0, IX2, IY2, IX0, IY0,    &
         IK, ITH, ISP, ISP0, ISP2
#ifdef W3_S
    INTEGER, SAVE           :: IENT = 0
#endif
#ifdef W3_T
    INTEGER                 :: MAPTXY(NY,NX), I, IXY
    INTEGER                 :: MAPTST(NK+2,NTH)
#endif
    !/
    !/ ------------------------------------------------------------------- /
    !/
#ifdef W3_S
    CALL strace (ient, 'W3MAP2')
#endif
    !
    ! 1.  Map MAPX2 ------------------------------------------------------ *
    ! 1.a Range 1 to NMX0
    !
#ifdef W3_T
    maptxy = 0.
#endif
    !
    nmx0   = 0
    DO ix=1, nx
      ixy0   = (ix-1)*ny
      ix2    = 1 + mod(ix,nx)
      DO iy=2, ny-1
        IF ( mapsta(iy,ix).EQ.1 .AND. mapsta(iy,ix2).EQ.1 ) THEN
          nmx0   = nmx0 + 1
          mapx2(nmx0) = ixy0 + iy
#ifdef W3_T
          maptxy(iy,ix) = maptxy(iy,ix) + 1
#endif
        END IF
      END DO
    END DO
    !
    ! 1.b Range NMX0+1 to NMX1
    !
    nmx1   = nmx0
    DO ix=1, nx
      ixy0   = (ix-1)*ny
      ix2    = 1 + mod(ix,nx)
      DO iy=2, ny-1
        IF ( mapsta(iy,ix).EQ.1 .AND. mapsta(iy,ix2).NE.1 ) THEN
          nmx1   = nmx1 + 1
          mapx2(nmx1) = ixy0 + iy
#ifdef W3_T
          maptxy(iy,ix) = maptxy(iy,ix) + 2
#endif
        END IF
      END DO
    END DO
    !
    ! 1.c Range NMX1+1 to NMX2
    !
    nmx2   = nmx1
    DO ix=1, nx
      ixy0   = (ix-1)*ny
      ix2    = 1 + mod(ix,nx)
      DO iy=2, ny-1
        IF ( mapsta(iy,ix).NE.1 .AND. mapsta(iy,ix2).EQ.1 ) THEN
          nmx2   = nmx2 + 1
          mapx2(nmx2) = ixy0 + iy
#ifdef W3_T
          maptxy(iy,ix) = maptxy(iy,ix) + 4
#endif
        END IF
      END DO
    END DO
    !
 
#ifdef W3_T
    WRITE (ndst,9000) 'MAPX2', nmx0, nmx1-nmx0,                  &
         nmx2-nmx1, nmx2
    DO iy=ny, 1, -1
      WRITE (ndst,9001) (maptxy(iy,ix),ix=1, nx)
    END DO
#endif
    !
    ! 2.  Map MAPY2 ------------------------------------------------------ *
    ! 2.a Range 1 to NMY0
    !
#ifdef W3_T
    maptxy = 0.
#endif
    !
    nmy0   = 0
    DO ix=1, nx
      ixy0   = (ix-1)*ny
      DO iy=1, ny-1
        iy2    = iy + 1
        IF ( mapsta(iy,ix).EQ.1 .AND. mapsta(iy2,ix).EQ.1 ) THEN
          nmy0   = nmy0 + 1
          mapy2(nmy0) = ixy0 + iy
#ifdef W3_T
          maptxy(iy,ix) = maptxy(iy,ix) + 1
#endif
        END IF
      END DO
    END DO
    !
    ! 2.b Range NMY0+1 to NMY1
    !
    nmy1   = nmy0
    DO ix=1, nx
      ixy0   = (ix-1)*ny
      DO iy=1, ny-1
        iy2    = iy + 1
        IF ( mapsta(iy,ix).EQ.1 .AND. mapsta(iy2,ix).NE.1 ) THEN
          nmy1   = nmy1 + 1
          mapy2(nmy1) = ixy0 + iy
#ifdef W3_T
          maptxy(iy,ix) = maptxy(iy,ix) + 2
#endif
        END IF
      END DO
    END DO
    !
    ! 2.c Range NMY1+1 to NMY2
    !
    nmy2   = nmy1
    DO ix=1, nx
      ixy0   = (ix-1)*ny
      DO iy=1, ny-1
        iy2    = iy + 1
        IF ( mapsta(iy,ix).NE.1 .AND. mapsta(iy2,ix).EQ.1 ) THEN
          nmy2   = nmy2 + 1
          mapy2(nmy2) = ixy0 + iy
#ifdef W3_T
          maptxy(iy,ix) = maptxy(iy,ix) + 4
#endif
        END IF
      END DO
    END DO
    !
#ifdef W3_T
    WRITE (ndst,9000) 'MAPY2', nmy0, nmy1-nmy0,                  &
         nmy2-nmy1, nmy2
    DO iy=ny, 1, -1
      WRITE (ndst,9001) (maptxy(iy,ix),ix=1, nx)
    END DO
#endif
    !
    ! 3.  Map MAPAXY and MAPXY ------------------------------------------- *
    !
    nact   = 0
    DO ix=1, nx
      iy0    = (ix-1)*ny
      DO iy=2, ny-1
        IF ( mapsta(iy,ix).EQ.1 ) THEN
          nact         = nact + 1
          mapaxy(nact) = iy0 + iy
        END IF
      END DO
    END DO
    !
    nmxy   = 0
    DO ix=1, nx
      ixy0   = (ix-1)*ny
      ix2    = ix+1
      IF ( ix .EQ. nx ) ix2 =  1
      DO iy=2, ny-1
        IF ( mapsta( iy ,ix ).GE.1 .AND.                            &
             mapsta( iy ,ix2).GE.1 .AND.                            &
             mapsta(iy+1,ix ).GE.1 .AND.                            &
             mapsta(iy+1,ix2).GE.1 ) THEN
          nmxy   = nmxy + 1
          mapxy(nmxy) = ixy0 + iy
        END IF
      END DO
    END DO
    !
    ! 4.  Maps for intra-spectral propagation ---------------------------- *
    !
    IF ( mapth2(1) .NE. 0 ) RETURN
    !
#ifdef W3_T
    maptst = 0
#endif
    !
    ! 4.a MAPTH2 and MAPBTK
    !
    DO ik=1, nk
      DO ith=1, nth
        isp    = ith + (ik-1)*nth
        isp2   = (ik+1) + (ith-1)*(nk+2)
        mapth2(isp) = isp2
#ifdef W3_T
        maptst(ik+1,ith) = maptst(ik+1,ith) + 1
#endif
      END DO
    END DO
    !
#ifdef W3_T
    WRITE (ndst,9000) 'MAPTH2', isp, 0, 0, isp
    DO ik=nk+2, 1, -1
      WRITE (ndst,9001) (maptst(ik,ith),ith=1, nth)
    END DO
    maptst = 0
#endif
    !
    ! 4.b MAPWN2
    !
    isp0   = 0
    DO ik=1, nk-1
      DO ith=1, nth
        isp0   = isp0 + 1
        isp2   = (ik+1) + (ith-1)*(nk+2)
        mapwn2(isp0) = isp2
#ifdef W3_T
        maptst(ik+1,ith) = maptst(ik+1,ith) + 1
#endif
      END DO
    END DO
    !
    DO ith=1, nth
      isp0   = isp0 + 1
      isp2   = nk+1 + (ith-1)*(nk+2)
      mapwn2(isp0) = isp2
#ifdef W3_T
      maptst(nk+1,ith) = maptst(nk+1,ith) + 2
#endif
    END DO
    !
    DO ith=1, nth
      isp0   = isp0 + 1
      isp2   = 1 + (ith-1)*(nk+2)
      mapwn2(isp0) = isp2
#ifdef W3_T
      maptst(1,ith) = maptst(1,ith) + 4
#endif
    END DO
    !
#ifdef W3_T
    WRITE (ndst,9000) 'MAPWN2', nspec-nth, nth, nth, nspec+nth
    DO ik=nk+2, 1, -1
      WRITE (ndst,9001) (maptst(ik,ith),ith=1, nth)
    END DO
#endif
    !
    RETURN
    !
    ! Formats
    !
#ifdef W3_T
9000 FORMAT (/' TEST W3MAP2 : TEST MAP FOR PROPAGATION'/          &
         '      MAP     : ',a/                               &
         '      CENTRAL : ',i6/                              &
         '      ABOVE   : ',i6/                              &
         '      BELOW   : ',i6/                              &
         '      TOTAL   : ',i6/)
9001 FORMAT (1x,130i1)
9010 FORMAT (' TEST W3MAP2 : COMPOSITE MAPS TH2, WN2 AND BTK')
9011 FORMAT (2x,60i2)
#endif
    !/
    !/ End of W3MAP2 ----------------------------------------------------- /
    !/

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

Referenced by w3wavemd::w3wave().

w3xyp2()

subroutine w3pro2md::w3xyp2	(	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 physical 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
[in]	VGY

Author

H. L. Tolman

Date

29-May-2014

Definition at line 509 of file w3pro2md.F90.

    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III           NOAA/NCEP |
    !/                  |           H. L. Tolman            |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         29-May-2014 |
    !/                  +-----------------------------------+
    !/
    !/    07-Jul-1998 : Final FORTRAN 77                    ( version 1.18 )
    !/    15-Dec-1999 : Upgrade to FORTRAN 90               ( version 2.00 )
    !/    24-Jan-2001 : Flat grid version                   ( version 2.06 )
    !/    09-Feb-2001 : Clean up of parameter lists         ( version 2.08 )
    !/    14-Feb-2001 : Unit numbers in UQ routines         ( version 2.08 )
    !/    13-Nov-2001 : Sub-grid obstructions.              ( version 2.14 )
    !/    26-Dec-2002 : Moving grid option,                 ( version 3.02 )
    !/    20-Dec-2004 : Multiple grid version.              ( version 3.06 )
    !/    07-Sep-2005 : Improved boundary conditions.       ( version 3.08 )
    !/    09-Nov-2005 : Removing soft boundary option.      ( version 3.08 )
    !/    05-Mar-2008 : Added NEC sxf90 compiler directives.
    !/                  (Chris Bunney, UK Met Office)       ( version 3.13 )
    !/    30-Oct-2009 : Implement curvilinear grid type.    ( version 3.14 )
    !/                  (W. E. Rogers & T. J. Campbell, NRL)
    !/    06-Dec-2010 : Change from GLOBAL (logical) to ICLOSE (integer) to
    !/                  specify index closure for a grid.   ( version 3.14 )
    !/                  (T. J. Campbell, NRL)
    !/    23-Dec-2010 : Fix HPFAC and HQFAC by including the COS(YGRD)
    !/                  factor with DXDP and DXDQ terms.    ( version 3.14 )
    !/                  (T. J. Campbell, NRL)
    !/    01-Jul-2013 : Adding UQ and UNO switches to chose between third
    !/                  and second order schemes.           ( version 4.12 )
    !/    29-May-2014 : Adding OMPH switch.                 ( version 5.02 )
    !/
    !  1. Purpose :
    !
    !     Propagation in physical space for a given spectral component.
    !
    !  2. Method :
    !
    !     Third-order ULTIMATE QUICKEST scheme and diffusion correction
    !     for linear dispersion (see manual).
    !     Curvilinear grid implementation: Fluxes are computed in index space
    !       and then transformed back into physical space.  The diffusion term
    !       is handled in physical space.
    !
    !  3. Parameters :
    !
    !     Parameter list
    !     ----------------------------------------------------------------
    !       ISP     Int.   I   Number of spectral bin (IK-1)*NTH+ITH
    !       DTG     Real   I   Total time step.
    !       MAPSTA  I.A.   I   Grid point status map.
    !       MAPFS   I.A.   I   Storage map.
    !       VQ      R.A.  I/O  Field to propagate.
    !     ----------------------------------------------------------------
    !
    !     Local variables.
    !     ----------------------------------------------------------------
    !       NTLOC   Int   Number of local time steps.
    !       DTLOC   Real  Local propagation time step.
    !       CGD     Real  Deep water group velocity.
    !       DSSD, DNND    Deep water diffusion coefficients.
    !       VLCFLX  R.A.  Local courant numbers in index space (norm. velocities)
    !       VLCFLY  R.A.
    !       CXTOT   R.A.  Propagation velocities in physical space.
    !       CYTOT   R.A.
    !       CELLP   Real  Cell Reynolds/Peclet number used to calculate
    !                     diffusion coefficient for growing spectral
    !                     components.
    !       DFRR    Real  Relative frequency increment.
    !     ----------------------------------------------------------------
    !
    !  4. Subroutines used :
    !
    !     See module documentation.
    !
    !  5. Called by :
    !
    !      Name      Type  Module   Description
    !     ----------------------------------------------------------------
    !      W3WAVE    Subr. W3WAVEMD Wave model routine.
    !     ---------------------------------------------------------------
    !
    !  6. Error messages :
    !
    !       None.
    !
    !  7. Remarks :
    !
    !     - Note that the ULTIMATE limiter does not guarantee non-zero
    !       energies.
    !     - The present scheme shows a strong distortion when propaga-
    !       ting a field under an angle with the grid in a truly 2-D
    !       fashion. Propagation is therefore split along the two
    !       axes.
    !     - Two boundary treatments are available. The first uses real
    !       boundaries in each space. In this case waves will not
    !       penetrate in narrow straights under an angle with the grid.
    !       This behavior is improved by using a 'soft' option, in
    !       which the 'X' or 'Y' sweep allows for energy to go onto
    !       the land. This improves the above behavior, but implies
    !       that X-Y connenctions are required in barriers for them
    !       to become inpenetrable.
    !     - If TDYN is set to zero, ALL diffusion is skipped. Set TDYN
    !       to a small positive number to have growth diffusion only.
    !     - Curvilinear grid implementation. Variables FACX, FACY, CCOS, CSIN,
    !       CCURX, CCURY are not needed and have been removed.  FACX is accounted
    !       for as approriate in this subroutine.  FACX is also accounted for in
    !       the case of .NOT.FLCX.  Since FACX is removed, there is now a check for
    !       .NOT.FLCX in this subroutine.  In CFL calcs dx and dy are omitted,
    !       since dx=dy=1 in index space.  Curvilinear grid derivatives
    !       (DPDY, DQDX, etc.) and metric (GSQRT) are brought in via W3GDATMD.
    !     - Factors VFDIFX_FAC VFDIFY_FAC VFDIFC_FAC are introduced so that results
    !       match for test case tp2.3.  Use of these factors is optional and removal
    !       can significantly reduce size/cost of code. These variants are marked as
    !       CURV1 or CURV2. NCEP will make final decision re: which version to adopt.
    !       CURV1 is the shorter version and results do not match the original code
    !       for all test cases. CURV2 is the longer version and results do match the
    !       original. DETAILED EXPLANATION: Discrepancies occur at the boundaries.
    !       This is because, at the boundaries, the pre-curvilinear version zeros out
    !       some terms in the diffusion calculation. Since they are zeroed out,
    !       they aren't there to *cancel* certain other terms: these "other terms"
    !       affect the result, so they have to be retained in the long vesion (CURV2)
    !       to get an exact match. In the short version, the "canceling out" is
    !       performed prior to coding the scheme, so both the canceled and canceling
    !       terms are always omitted.
    !
    !  8. Structure :
    !
    !     ---------------------------------------------
    !       1.  Preparations
    !         a Set constants
    !         b Initialize arrays
    !       2.  Prepare arrays
    !         a Velocities and 'Q'
    !         b diffusion coefficients
    !       3.  Loop over sub-steps
    !       ----------------------------------------
    !         a Propagate
    !         b Update boundary conditions
    !         c Diffusion correction
    !       ----------------------------------------
    !       4.  Store Q field in spectra
    !     ---------------------------------------------
    !
    !  9. Switches :
    !
    !       !/MGP   Correct for motion of grid.
    !
    !       !/TDYN  Dynamic increase of DTME
    !       !/DSS0  Disable diffusion in propagation direction
    !       !/XW0   Propagation diffusion only.
    !       !/XW1   Growth diffusion only.
    !
    !       !/OMPH  Hybrid OpenMP directives.
    !
    !       !/S     Enable subroutine tracing.
    !
    !       !/T     Enable general test output.
    !       !/T1    Dump of input field and fluxes.
    !       !/T2    Dump of output field.
    !
    ! 10. Source code :
    !
    !/ ------------------------------------------------------------------- /
    USE constants
    !
    USE w3timemd, ONLY: dsec21
    !
    USE w3gdatmd, ONLY: nk, nth, dth, xfr, esin, ecos, sig, nx, ny, &
         nsea, sx, sy, mapsf, iclose, flcx, flcy,    &
         iclose_none, iclose_smpl, iclose_trpl,      &
         dtcfl, clats, dtme, clatmn, flagll,         &
         hpfac, hqfac, dpdx, dpdy, dqdx, dqdy, gsqrt
    USE w3wdatmd, ONLY: time
    USE w3adatmd, ONLY: cg, wn, u10, cx, cy, atrnx, atrny, itime,   &
         nmx0, nmx1, nmx2, nmy0, nmy1, nmy2, nact,   &
         nmxy, mapx2, mapy2, mapaxy, mapxy
    USE w3idatmd, ONLY: flcur
    USE w3odatmd, ONLY: ndse, ndst, flbpi, nbi, tbpi0, tbpin,       &
         isbpi, bbpi0, bbpin, iaproc, naperr
    USE w3servmd, ONLY: extcde
#ifdef W3_S
    USE w3servmd, ONLY: strace
#endif
#ifdef W3_UQ
    USE w3uqckmd
#endif
#ifdef W3_UNO
    USE w3uno2md
#endif
    !/
    IMPLICIT NONE
    !/
    !/ ------------------------------------------------------------------- /
    !/ Parameter list
    !/
    INTEGER, INTENT(IN)     :: ISP, MAPSTA(NY*NX), MAPFS(NY*NX)
    REAL, INTENT(IN)        :: DTG, VGX, VGY
    REAL, INTENT(INOUT)     :: VQ(1-NY:NY*(NX+2))
    !/
    !/ ------------------------------------------------------------------- /
    !/ Local parameters
    !/
    INTEGER                 :: ITH, IK, NTLOC, ITLOC, ISEA, IXY,    &
         IX,IY, IY0, IP, IBI
    INTEGER                 :: TTEST(2),DTTST
#ifdef W3_S
    INTEGER, SAVE           :: IENT = 0
#endif
    REAL                    :: CG0, CGA, CGN, CGX, CGY, CXC, CYC,  &
         CXMIN, CXMAX, CYMIN, CYMAX
    REAL                    :: DTLOC, DTRAD,               &
         DFRR, CELLP,  CGD, DSSD,    &
         DNND, DCELL, XWIND, TFAC, DSS, DNN
    REAL                    :: RD1, RD2
    REAL                    :: RFAC, DFAC, DVQ, QXX, QXY, QYY
    REAL                    :: CP, CQ
    LOGICAL                 :: YFIRST
    LOGICAL                 :: GLOBAL
    !/
    !/ Automatic work arrays
    !/
    REAL                    :: VLCFLX((NX+1)*NY), VLCFLY((NX+1)*NY),&
         VFDIFX(1-NY:NX*NY), VFDIFY(NX*NY),   &
         VFDIFC(1-NY:NX*NY), VDXX((NX+1)*NY), &
         VDYY(NX*NY), VDXY((NX+1)*NY)
 
    REAL                    :: CXTOT((NX+1)*NY), CYTOT(NX*NY)
    REAL                    :: VQ_OLD(1-NY:NY*(NX+2))
    !CURV2: BEGIN -----------------------------------------------------------------
    REAL                    :: VFDIFX_FAC(1-NY:NX*NY),              &
         VFDIFY_FAC(1-NY:NX*NY),              &
         VFDIFC_FAC(1-NY:NX*NY)
    !CURV2: END -------------------------------------------------------------------
    !/
    !/ ------------------------------------------------------------------- /
    !/
#ifdef W3_S
    CALL strace (ient, 'W3XYP2')
#endif
    !
    !     IF ( MAXVAL(VQ) .EQ. 0. ) THEN
    !         IF ( NBI .EQ. 0 ) THEN
    !             RETURN
    !           ELSE
    !             IF ( MAXVAL(BBPI0(ISP,:)) .EQ. 0. .AND.            &
    !                  MAXVAL(BBPIN(ISP,:)) .EQ. 0. ) RETURN
    !           END IF
    !       END IF
    !
    ! 1.  Preparations --------------------------------------------------- *
 
    IF ( iclose .EQ. iclose_trpl ) THEN
      IF (iaproc .EQ. naperr) &
           WRITE(ndse,*)'SUBROUTINE W3XYP2 IS NOT YET ADAPTED FOR '//    &
           'TRIPOLE GRIDS. STOPPING NOW.'
      CALL extcde ( 1 )
    END IF
 
    ! 1.a Set constants
    !
    global = iclose.NE.iclose_none
    ith    = 1 + mod(isp-1,nth)
    ik     = 1 + (isp-1)/nth
    !
    cg0    = 0.575 * grav / sig(1)
    cga    = 0.575 * grav / sig(ik)
    cgx    = cga * ecos(ith)
    cgy    = cga * esin(ith)
#ifdef W3_MGP
    cgx    = cgx - vgx
    cgy    = cgy - vgy
#endif
    !
    IF ( flcur ) THEN
      cxmin  = minval( cx(1:nsea) )
      cxmax  = maxval( cx(1:nsea) )
      cymin  = minval( cy(1:nsea) )
      cymax  = maxval( cy(1:nsea) )
      IF ( abs(cgx+cxmin) .GT. abs(cgx+cxmax) ) THEN
        cgx    = cgx + cxmin
      ELSE
        cgx    = cgx + cxmax
      END IF
      IF ( abs(cgy+cymin) .GT. abs(cgy+cymax) ) THEN
        cgy    = cgy + cymin
      ELSE
        cgy    = cgy + cymax
      END IF
      cxc    = max( abs(cxmin) , abs(cxmax) )
      cyc    = max( abs(cymin) , abs(cymax) )
#ifdef W3_MGP
      cxc    = max( abs(cxmin-vgx) , abs(cxmax-vgx) )
      cyc    = max( abs(cymin-vgy) , abs(cymax-vgy) )
#endif
    ELSE
      cxc    = 0.
      cyc    = 0.
    END IF
    !
    cgn    = 0.9999 * max( abs(cgx) , abs(cgy) , cxc, cyc, 0.001*cg0 )
    !
    ntloc  = 1 + int(dtg/(dtcfl*cg0/cgn))
    dtloc  = dtg / real(ntloc)
    dtrad  = dtloc
    IF ( flagll ) dtrad=dtrad/(dera*radius)
    !
    IF ( flagll ) THEN
      rfac = dera * radius
      dfac = 1. / rfac**2
    ELSE
      rfac = 1.
      dfac = 1.
    END IF
    !
    ttest(1) = time(1)
    ttest(2) = 0
    dttst = dsec21(ttest,time)
    yfirst = mod(nint(dttst/dtg),2) .EQ. 0
    !
#ifdef W3_T
    WRITE (ndst,9000) yfirst
    WRITE (ndst,9001) isp, ith, ik, ecos(ith), esin(ith)
#endif
    !
#ifdef W3_TDYN
    IF ( isp .EQ. 1 ) dtme = dtme + dtg
#endif
    !
    IF ( dtme .NE. 0. ) THEN
      dfrr   = xfr - 1.
      cellp  = 10.
      cgd    = 0.5 * grav / sig(ik)
      dssd   = ( dfrr * cgd )**2 * dtme / 12.
      dnnd   = ( cgd * dth )**2 * dtme / 12.
#ifdef W3_T
      WRITE (ndst,9002) dfrr, cellp, dtme
    ELSE
      WRITE (ndst,9003)
#endif
    END IF
 
    !
    ! 1.b Initialize arrays
    !
#ifdef W3_T
    WRITE (ndst,9010)
#endif
    !
    vlcflx = 0.
    vlcfly = 0.
    vfdifx = 0.
    vfdify = 0.
    vfdifc = 0.
    vdxx   = 0.
    vdyy   = 0.
    vdxy   = 0.
    cxtot  = 0.
    cytot  = 0.
    !
    ! 2.  Calculate velocities and diffusion coefficients ---------------- *
    ! 2.a Velocities
    !
    !     Q     = ( A / CG * CLATS )
    !     LCFLX = ( COS*CG / CLATS ) * DT / DX
    !     LCFLY = (     SIN*CG )     * DT / DX
    !
#ifdef W3_T
    WRITE (ndst,9020) nsea
#endif
    !
#ifdef W3_OMPH
    !$OMP PARALLEL DO PRIVATE (ISEA, IXY)
#endif
    !
    DO isea=1, nsea
      ixy         = mapsf(isea,3)
      vq(ixy) = vq(ixy) / cg(ik,isea) * clats(isea)
      cxtot(ixy) = ecos(ith) * cg(ik,isea) / clats(isea)
      cytot(ixy) = esin(ith) * cg(ik,isea)
#ifdef W3_MGP
      cxtot(ixy) = cxtot(ixy) - vgx/clats(isea)
      cytot(ixy) = cytot(ixy) - vgy
#endif
#ifdef W3_T1
      IF ( .NOT. flcur )                                        &
           WRITE (ndst,9021) isea, mapsf(isea,1), mapsf(isea,2), &
           vq(ixy), cxtot(ixy), cytot(ixy)
#endif
    END DO
    !
#ifdef W3_OMPH
    !$OMP END PARALLEL DO
#endif
    !
    IF ( flcur ) THEN
#ifdef W3_T
      WRITE (ndst,9022)
#endif
      DO isea=1, nsea
        ixy         = mapsf(isea,3)
        cxtot(ixy) = cxtot(ixy) + cx(isea)/clats(isea)
        cytot(ixy) = cytot(ixy) + cy(isea)
#ifdef W3_T1
        WRITE (ndst,9021) isea, mapsf(isea,1), mapsf(isea,2), &
             vq(ixy), cxtot(ixy), cytot(ixy)
#endif
      END DO
    END IF
 
    !
#ifdef W3_OMPH
    !$OMP PARALLEL DO PRIVATE (ISEA, IX, IY, IXY, CP, CQ)
#endif
    !
    DO isea=1, nsea
      ix     = mapsf(isea,1)
      iy     = mapsf(isea,2)
      ixy    = mapsf(isea,3)
      cp=cxtot(ixy)*dpdx(iy,ix)+cytot(ixy)*dpdy(iy,ix)
      cq=cxtot(ixy)*dqdx(iy,ix)+cytot(ixy)*dqdy(iy,ix)
      vlcflx(ixy) = cp*dtrad
      vlcfly(ixy) = cq*dtrad
    END DO
    !
#ifdef W3_OMPH
    !$OMP END PARALLEL DO
#endif
    !
    ! 2.b Diffusion coefficients
    !
    IF ( dtme .NE. 0. ) THEN
      !
#ifdef W3_OMPH
      !$OMP PARALLEL DO PRIVATE (ISEA, IX, IY, IXY, &
      !$OMP&                     DCELL, XWIND, TFAC, DSS, DNN)
#endif
      !
      DO isea=1, nsea
        ix        = mapsf(isea,1)
        iy        = mapsf(isea,2)
        ixy       = mapsf(isea,3)
        IF ( min( atrnx(ixy,1) , atrnx(ixy,-1) ,                 &
             atrny(ixy,1) , atrny(ixy,-1) ) .GT. trnmin ) THEN
          dcell     = cgd * min( hpfac(iy,ix)*rfac, &
               hqfac(iy,ix)*rfac ) / cellp
          xwind     = 3.3 * u10(isea)*wn(ik,isea)/sig(ik) - 2.3
          xwind     = max( 0. , min( 1. , xwind ) )
#ifdef W3_XW0
          xwind     = 0.
#endif
#ifdef W3_XW1
          xwind     = 1.
#endif
          tfac      = min( 1. , (clats(isea)/clatmn)**2 )
          dss       = xwind * dcell + (1.-xwind) * dssd * tfac
#ifdef W3_DSS0
          dss       = 0.
#endif
          dnn       = xwind * dcell + (1.-xwind) * dnnd * tfac
 
          vdxx(ixy) = dtloc * (dss*ecos(ith)**2+dnn*esin(ith)**2)
          vdyy(ixy) = dtloc * (dss*esin(ith)**2+dnn*ecos(ith)**2) &
               / clats(isea)**2
          vdxy(ixy) = dtloc * (dss-dnn) * esin(ith)*ecos(ith)     &
               / clats(isea)
 
        END IF
      END DO
      !
#ifdef W3_OMPH
      !$OMP END PARALLEL DO
#endif
      !
    END IF
    !
    ! 3.  Loop over sub-steps -------------------------------------------- *
    !
    DO itloc=1, ntloc
      !
      ! 3.a Propagate fields
      !
#ifdef W3_OMPH
      !$OMP PARALLEL DO PRIVATE (ISEA, IX, IY, IXY )
#endif
      !
      DO isea=1, nsea
        ix     = mapsf(isea,1)
        iy     = mapsf(isea,2)
        ixy    = mapsf(isea,3)
        vq(ixy)= vq(ixy) * gsqrt(iy,ix)
      END DO
      !
#ifdef W3_OMPH
      !$OMP END PARALLEL DO
#endif
      !
      IF ( yfirst ) THEN
        !
#ifdef W3_UQ
        IF ( flcy ) CALL w3qck3                               &
             (nx, ny, nx, ny, vlcfly, atrny, vq,       &
             .false., 1, mapaxy, nact, mapy2, nmy0,   &
             nmy1, nmy2, ndse, ndst )
        IF ( flcx ) CALL w3qck3                               &
             (nx, ny, nx, ny, vlcflx, atrnx, vq,       &
             global, ny, mapaxy, nact, mapx2, nmx0,   &
             nmx1, nmx2, ndse, ndst )
#endif
        !
#ifdef W3_UNO
        IF ( flcy ) CALL w3uno2s                             &
             (nx, ny, nx, ny, vlcfly, atrny, vq,      &
             .false., 1, mapaxy, nact, mapy2, nmy0,  &
             nmy1, nmy2, ndse, ndst )
        IF ( flcx ) CALL w3uno2s                             &
             (nx, ny, nx, ny, vlcflx, atrnx, vq,      &
             global, ny, mapaxy, nact, mapx2, nmx0,  &
             nmx1, nmx2, ndse, ndst )
#endif
        !
      ELSE
        !
#ifdef W3_UQ
        IF ( flcx ) CALL w3qck3                               &
             (nx, ny, nx, ny, vlcflx, atrnx, vq,       &
             global, ny, mapaxy, nact, mapx2, nmx0,   &
             nmx1, nmx2, ndse, ndst )
        IF ( flcy ) CALL w3qck3                               &
             (nx, ny, nx, ny, vlcfly, atrny, vq,       &
             .false., 1, mapaxy, nact, mapy2, nmy0,   &
             nmy1, nmy2, ndse, ndst )
#endif
        !
#ifdef W3_UNO
        IF ( flcx ) CALL w3uno2s                             &
             (nx, ny, nx, ny, vlcflx, atrnx, vq,      &
             global, ny, mapaxy, nact, mapx2, nmx0,  &
             nmx1, nmx2, ndse, ndst )
        IF ( flcy ) CALL w3uno2s                             &
             (nx, ny, nx, ny, vlcfly, atrny, vq,      &
             .false., 1, mapaxy, nact, mapy2, nmy0,  &
             nmy1, nmy2, ndse, ndst )
#endif
        !
      END IF
      !
#ifdef W3_OMPH
      !$OMP PARALLEL DO PRIVATE (ISEA, IX, IY, IXY )
#endif
      !
      DO isea=1, nsea
        ix     = mapsf(isea,1)
        iy     = mapsf(isea,2)
        ixy    = mapsf(isea,3)
        vq(ixy)= vq(ixy) / gsqrt(iy,ix)
      END DO
      !
#ifdef W3_OMPH
      !$OMP END PARALLEL DO
#endif
      !
      ! 3.b Update boundaries
      !
      IF ( flbpi ) THEN
        rd1    = dsec21( tbpi0, time ) - dtg *                   &
             REAL(NTLOC-ITLOC)/REAL(NTLOC)
        rd2    = dsec21( tbpi0, tbpin )
        IF ( rd2 .GT. 0.001 ) THEN
          rd2    = min(1.,max(0.,rd1/rd2))
          rd1    = 1. - rd2
        ELSE
          rd1    = 0.
          rd2    = 1.
        END IF
        DO ibi=1, nbi
          isea    = isbpi(ibi)
          ixy     = mapsf(isbpi(ibi),3)
          vq(ixy) = ( rd1*bbpi0(isp,ibi) + rd2*bbpin(isp,ibi) )   &
               / cg(ik,isea) * clats(isea)
        END DO
      END IF
      !
      ! 3.c Diffusion correction
      !
      IF ( dtme .NE. 0. ) THEN
 
        IF ( global ) THEN
          DO iy=1, ny
            vq(iy+nx*ny) = vq(iy)
          END DO
        END IF
 
        !CURV2: BEGIN -----------------------------------------------------------------
        vfdifx_fac=0.0
        DO ip=1, nmx0
          ixy         = mapx2(ip)
          vfdifx_fac(ixy) = 1.0
        END DO
        vfdify_fac=0.0
        DO ip=1, nmy0
          ixy         = mapy2(ip)
          vfdify_fac(ixy) = 1.0
        END DO
        vfdifc_fac=0.0
        DO ip=1, nmxy
          ixy         = mapxy(ip)
          vfdifc_fac(ixy) = 1.0
        END DO
        IF ( global ) THEN
          iy0     = (nx-1)*ny
          DO iy=1, ny
            vfdifx_fac(iy-ny) = vfdifx_fac(iy+iy0)
            vfdifc_fac(iy-ny) = vfdifc_fac(iy+iy0)
          END DO
        END IF
        !CURV2: END -------------------------------------------------------------------
 
        vq_old = vq
        !
#ifdef W3_OMPH
        !$OMP PARALLEL DO PRIVATE (ISEA, IX, IY, IXY, &
        !$OMP&                     QXX, QYY, QXY, DVQ )
#endif
        !
        DO ip=1, nact
          ixy    = mapaxy(ip)
          isea   = mapfs(ixy)
          ix     = mapsf(isea,1)
          iy     = mapsf(isea,2)
 
          !CURV1: DOES NOT GIVE EXACT MATCH TO EARLIER WW3 VERSION FOR TEST CASE TP2.3
          !CURV1: NEAR THE BOUNDARY, NOTE THAT THIS VERSION USES NON-ACTIVE GRID POINTS
          !CURV1: IN ITS CALCS. THIS IS NO PROBLEM, AS LONG AS THE NON-ACTIVE GRID POINTS
          !CURV1: EXIST IN THE ARRAY AND ARE VQ=0. ALSO NOTE THAT WITH THE SHORT VERSION
          !CURV1: OF THE CODE, VFDIF?_FAC VARIABLES CAN BE REMOVED AND 3-4 DO LOOPS ABOVE
          !CURV1: CAN BE REMOVED.
          !CURV1: BEGIN -----------------------------------------------------------------
          !                QXX = VQ_OLD(IXY+NY) - 2.0*VQ_OLD(IXY) + VQ_OLD(IXY-NY)
          !                QYY = VQ_OLD(IXY+1)  - 2.0*VQ_OLD(IXY) + VQ_OLD(IXY-1)
          !                QXY = VQ_OLD(IXY+NY+1) - VQ_OLD(IXY-NY+1) &
          !                    - VQ_OLD(IXY+NY-1) + VQ_OLD(IXY-NY-1)
          !CURV1: END -------------------------------------------------------------------
 
          !CURV2: DOES GIVE EXACT MATCH TO EARLIER WW3 VERSION FOR TEST CASE TP2.3. NOTE
          !CURV2: THAT IF VFDIFC_FAC VARIABLES ARE ALL UNITY, MANY TERMS CANCEL OUT.
          !CURV2: HOWEVER, VFDIFC_FAC IS ZERO WHEN A RELATED VQ POINT IS NOT AN ACTIVE
          !CURV2: GRID POINT
          !CURV2: BEGIN -----------------------------------------------------------------
          qxx = vfdifx_fac(ixy)   *vq_old(ixy+ny)     &
               - vfdifx_fac(ixy)   *vq_old(ixy)        &
               - vfdifx_fac(ixy-ny)*vq_old(ixy)        &
               + vfdifx_fac(ixy-ny)*vq_old(ixy-ny)
          qyy = vfdify_fac(ixy)   *vq_old(ixy+1)      &
               - vfdify_fac(ixy)   *vq_old(ixy)        &
               - vfdify_fac(ixy-1) *vq_old(ixy)        &
               + vfdify_fac(ixy-1) *vq_old(ixy-1)
          qxy = vfdifc_fac(ixy)     *vq_old(ixy)      &
               + vfdifc_fac(ixy-ny-1)*vq_old(ixy)      &
               - vfdifc_fac(ixy-1)   *vq_old(ixy)      &
               - vfdifc_fac(ixy-ny)  *vq_old(ixy)      &
               + vfdifc_fac(ixy-ny)  *vq_old(ixy+1)    &
               - vfdifc_fac(ixy)     *vq_old(ixy+1)    &
               + vfdifc_fac(ixy-1)   *vq_old(ixy-1)    &
               - vfdifc_fac(ixy-ny-1)*vq_old(ixy-1)    &
               + vfdifc_fac(ixy-1)   *vq_old(ixy+ny)   &
               - vfdifc_fac(ixy)     *vq_old(ixy+ny)   &
               + vfdifc_fac(ixy-ny)  *vq_old(ixy-ny)   &
               - vfdifc_fac(ixy-ny-1)*vq_old(ixy-ny)   &
               + vfdifc_fac(ixy)     *vq_old(ixy+ny+1) &
               - vfdifc_fac(ixy-1)   *vq_old(ixy+ny-1) &
               + vfdifc_fac(ixy-ny-1)*vq_old(ixy-ny-1) &
               - vfdifc_fac(ixy-ny)  *vq_old(ixy-ny+1)
          !CURV2: END -------------------------------------------------------------------
          !
          qxy = 0.25*qxy
          !
          dvq = vdxx(ixy)*(  dpdx(iy,ix)*dpdx(iy,ix)*qxx          &
               + 2.0*dpdx(iy,ix)*dqdx(iy,ix)*qxy      &
               + dqdx(iy,ix)*dqdx(iy,ix)*qyy )        &
               + vdyy(ixy)*(  dpdy(iy,ix)*dpdy(iy,ix)*qxx          &
               + 2.0*dpdy(iy,ix)*dqdy(iy,ix)*qxy      &
               + dqdy(iy,ix)*dqdy(iy,ix)*qyy)         &
               + 2.0*vdxy(ixy)*(  dpdx(iy,ix)*dpdy(iy,ix)*qxx      &
               + dqdx(iy,ix)*dqdy(iy,ix)*qyy      &
               +      (  dpdx(iy,ix)*dqdy(iy,ix)  &
               + dqdx(iy,ix)*dpdy(iy,ix) )*qxy )
          !
          vq(ixy) = vq_old(ixy) + dvq * dfac
          !
        END DO
        !
#ifdef W3_OMPH
        !$OMP END PARALLEL DO
#endif
        !
      END IF
      !
      yfirst = .NOT. yfirst
    END DO
    !
    ! 4.  Store results in VQ in proper format --------------------------- *
    !
#ifdef W3_T
    WRITE (ndst,9040) nsea
#endif
    !
#ifdef W3_OMPH
    !$OMP PARALLEL DO PRIVATE (ISEA, IXY )
#endif
    !
    DO isea=1, nsea
      ixy    = mapsf(isea,3)
      IF ( mapsta(ixy) .GT. 0 ) THEN
#ifdef W3_T2
        WRITE (ndst,9041) isea, mapsf(isea,1), mapsf(isea,2), vq(ixy)
#endif
        vq(ixy) =  max( 0. , cg(ik,isea) / clats(isea) * vq(ixy) )
        !         ELSE
        !           VQ(IXY) =  0.
      END IF
    END DO
    !
#ifdef W3_OMPH
    !$OMP END PARALLEL DO
#endif
    !
    RETURN
    !
    ! Formats
    !
#ifdef W3_T
9000 FORMAT (' TEST W3XYP2 : YFIRST  :',l2)
9001 FORMAT (' TEST W3XYP2 : ISP, ITH, IK, COS-SIN :',i8,2i4,2f7.3)
9002 FORMAT (' TEST W3XYP2 : DFRR, CELLP, DTME     :',3e10.3)
9003 FORMAT (' TEST W3XYP2 : NO DISPERSION CORRECTION ')
9010 FORMAT (' TEST W3XYP2 : INITIALIZE ARRAYS')
9020 FORMAT (' TEST W3XYP2 : CALCULATING LCFLX/Y AND DSS/NN (NSEA=', &
         i6,')')
9022 FORMAT (' TEST W3XYP2 : CORRECTING FOR CURRENT')
9040 FORMAT (' TEST W3XYP2 : FIELD AFTER PROP. (NSEA=',i6,')')
#endif
#ifdef W3_T1
9021 FORMAT (1x,i6,2i5,e12.4,2f7.3)
#endif
#ifdef W3_T2
9041 FORMAT (1x,i6,2i5,e12.4)
#endif
    !/
    !/ End of W3XYP2 ----------------------------------------------------- /
    !/

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

Referenced by w3wavemd::w3wave().