w3pro1md Module Reference

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

Functions/Subroutines
subroutine	w3map1 (MAPSTA)
	Generate 'map' arrays for the first order upstream scheme. More...
subroutine	w3xyp1 (ISP, DTG, MAPSTA, FIELD, VGX, VGY)
	Propagation in physical space for a given spectral component. More...
subroutine	w3ktp1 (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 first order propagation scheme in single module.

Author

H. L. Tolman

Date

05-Jun-2018

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

w3ktp1()

subroutine w3pro1md::w3ktp1	(	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.

Parameters

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

Author

H. L. Tolman

Date

20-Dec-2004

Definition at line 900 of file w3pro1md.F90.

    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III           NOAA/NCEP |
    !/                  |           H. L. Tolman            |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         20-Dec-2004 |
    !/                  +-----------------------------------+
    !/
    !/    29-Aug-1997 : Final FORTRAN 77                    ( version 1.18 )
    !/    04-Feb-2000 : Upgrade to FORTRAN 90               ( version 2.00 )
    !/    20-Dec-2004 : Multiple grid version.              ( version 3.06 )
    !/
    !  1. Purpose :
    !
    !     Propagation in spectral space.
    !
    !  2. Method :
    !
    !     First order scheme.
    !
    !  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.
    !       CG      R.A.   I   Local group velocities.
    !       WN      R.A.   I   Local wavenumbers.
    !       DEPTH   Real   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.
    !       FRK, FRG, FKC
    !               R.A.  Partial velocity terms.
    !       DWNI    R.A.  Inverse band width.
    !       CTH-WN  R.A.  Propagation velocities of local fluxes.
    !       FLTH-WN R.A.  Propagation fluxes.
    !       AA      R.A.  Extracted spectrum
    !     ----------------------------------------------------------------
    !
    !  4. Subroutines used :
    !
    !     See module documentation.
    !
    !  5. Called by :
    !
    !      Name      Type  Module   Description
    !     ----------------------------------------------------------------
    !      W3WAVE    Subr. W3WAVEMD Wave model routine.
    !     ----------------------------------------------------------------
    !
    !  6. Error messages :
    !
    !       None.
    !
    !  8. Structure :
    !
    !     -----------------------------------------------------------------
    !       1.  Preparations
    !         a Calculate DSDD
    !         b Extract spectrum
    !       2.  Refraction velocities
    !         a Filter level depth reffraction.
    !         b Depth refratcion velocity.
    !         c Current refraction velocity.
    !       3.  Wavenumber shift velocities
    !         a Prepare directional arrays
    !         b Calcuate velocity.
    !       4.  Refraction
    !         a Discrete fluxes.
    !         b Propagation fluxes.
    !         c Refraction.
    !       5.  Wavenumber shifts.
    !         a Discrete fluxes.
    !         b Propagation fluxes.
    !         c Refraction.
    !     -----------------------------------------------------------------
    !
    !  9. Switches :
    !
    !     C/S   Enable subroutine tracing.
    !     C/T   Enable general test output.
    !
    ! 10. Source code :
    !
    !/ ------------------------------------------------------------------- /
    USE constants
    USE w3gdatmd, ONLY: nk, nth, nspec, sig, dsip, ecos, esin, es2, &
         esc, ec2, fachfa, mapwn, flcth, flck, ctmax
    USE w3adatmd, ONLY: is0, is2
    USE w3idatmd, ONLY: flcur
    USE w3odatmd, ONLY: ndst
#ifdef W3_DEBUG
    USE w3odatmd, only : iaproc
#endif
#ifdef W3_S
    USE w3servmd, ONLY: strace
#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, ITH0
    REAL                    :: FDDMAX, FDG, DCYX, DCXXYY, DCXY,     &
         DCXX, DCXYYX, DCYY, FKD, FKD0, CTHB, &
         CWNB
    REAL                    :: VCTH(NSPEC), VCWN(1-NTH:NSPEC+NTH),  &
         VAA(1-NTH:NSPEC+NTH), VFLTH(NSPEC),  &
         VFLWN(1-NTH:NSPEC), DSDD(0:NK+1),    &
         FRK(NK), FRG(NK), FKC(NTH), DWNI(NK)
#ifdef W3_S
    INTEGER, SAVE           :: IENT = 0
    CALL strace (ient, 'W3KTP1')
#endif
    !/
    !/ ------------------------------------------------------------------- /
    !/
#ifdef W3_T
    WRITE (ndst,9000) flcth, flck, facth, fack, ctmax
    WRITE (ndst,9001) isea, depth, cx, cy,                       &
         dddx, dddy, dcxdx, dcxdy, dcydx, dcydy
#endif
    !
    ! 1.  Preparations --------------------------------------------------- *
    ! 1.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,9010)
    DO ik=1, nk+1
      WRITE (ndst,9011) ik, tpi/sig(ik), tpi/wn(ik),             &
           cg(ik), dsdd(ik)
    END DO
#endif
    !
    ! 1.b Extract spectrum
    !
    DO isp=1, nspec
      vaa(isp) = va(isp)
    END DO
    !
    ! 2.  Refraction velocities ------------------------------------------ *
    !
    IF ( flcth ) THEN
      !
      ! 2.a Set slope filter for depth refraction
      !
      fddmax = 0.
      fdg    = facth * cthg0
      !
      DO ith=1, nth
        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
      !
      ! 2.b Depth refraction and great-circle propagation
      !
      DO isp=1, nspec
        vcth(isp) = frg(mapwn(isp)) * ecos(isp)                   &
             + frk(mapwn(isp)) * ( esin(isp)*dddx - ecos(isp)*dddy )
      END DO
      !
#ifdef W3_REFRX
      ! 2.c @C/@x refraction and great-circle propagation
      vcth = 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)
        frk(ik) = frk(ik) / max( 1. , frk(ik)*fddmax/ctmax )
        frg(ik) = fdg * cg(ik)
      END DO
      DO isp=1, nspec
        vcth(isp) = frg(mapwn(isp)) * ecos(isp)            &
             + frk(mapwn(isp)) * ( esin(isp)*dcdx(mapwn(isp)) &
             - ecos(isp)*dcdy(mapwn(isp)) )
      END DO
#endif
      !
      ! 2.d Current refraction
      !
      IF ( flcur ) THEN
        !
        dcyx   = facth *   dcydx
        dcxxyy = facth * ( dcxdx - dcydy )
        dcxy   = facth *   dcxdy
        !
        DO isp=1, nspec
          vcth(isp) = vcth(isp) + es2(isp)*dcyx                 &
               + esc(isp)*dcxxyy - ec2(isp)*dcxy
        END DO
        !
      END IF
      !
    END IF
    !
    ! 3.  Wavenumber shift velocities ------------------------------------ *
    !
    IF ( flck ) THEN
      !
      dcxx   =  - fack *   dcxdx
      dcxyyx =  - fack * ( dcxdy + dcydx )
      dcyy   =  - fack *   dcydy
      fkd    =    fack * ( cx*dddx + cy*dddy )
      !
      DO ith=1, nth
        fkc(ith) = ec2(ith)*dcxx +                                &
             esc(ith)*dcxyyx + es2(ith)*dcyy
      END DO
      !
      isp    = -nth
      DO ik=0, nk+1
        fkd0   = fkd / cg(ik) * dsdd(ik)
        DO ith=1, nth
          isp    = isp + 1
          vcwn(isp) = fkd0 + wn(ik)*fkc(ith)
        END DO
      END DO
      !
      ith0   = nspec - nth
      DO ith=1, nth
        vaa(ith+nspec) = fachfa * vaa(ith+ith0)
        vaa(ith- nth ) = 0.
      END DO
      !
      DO ik=1, nk
        dwni(ik) = cg(ik) / dsip(ik)
      END DO
      !
    END IF
    !
    ! 4.  Refraction ----------------------------------------------------- *
    !
    IF ( flcth ) THEN
      !
      ! 4.a Boundary velocity and fluxes
      !
      DO isp=1, nspec
        cthb       = 0.5 * ( vcth(isp) + vcth(is2(isp)) )
        vflth(isp) = max( cthb , 0. ) * vaa(isp)                 &
             + min( cthb , 0. ) * vaa(is2(isp))
      END DO
      !
      ! 4.b Propagation
      !
      DO isp=1, nspec
        va(isp) = va(isp) + vflth(is0(isp)) - vflth(isp )
      END DO
      !
    END IF
    !
    ! 5.  Wavenumber shifts ---------------------------------------------- *
    !
    IF ( flck ) THEN
      !
      ! 5.a Boundary velocity and fluxes
      !
      DO isp=1-nth, nspec
        cwnb       = 0.5 * ( vcwn(isp) + vcwn(isp+nth) )
        vflwn(isp) = max( cwnb , 0. ) * vaa(  isp  )             &
             + min( cwnb , 0. ) * vaa(isp+nth)
      END DO
      !
      ! 5.c Propagation
      !
      DO isp=1, nspec
        va(isp) = va(isp) + dwni(mapwn(isp)) *                    &
             ( vflwn(isp-nth) - vflwn(isp) )
      END DO
      !
    END IF
    !
    RETURN
    !
    ! Formats
    !
#ifdef W3_T
9000 FORMAT (' TEST W3KTP1 : FLCTH-K, FACTH-K, CTMAX  :',         &
         2l2,2e10.3,f7.3)
9001 FORMAT (' TEST W3KTP1 : LOCAL DATA :',i7,f7.1,2f6.2,1x,      &
         6e10.3)
9010 FORMAT (' TEST W3KTP1 : IK, T, L, CG, DSDD : ')
9011 FORMAT ('              ',i3,f7.2,f7.1,f7.2,e11.3)
#endif
    !/
    !/ End of W3KTP1 ----------------------------------------------------- /
    !/

References w3gdatmd::ctmax, w3gdatmd::dsip, w3gdatmd::ec2, w3gdatmd::ecos, w3gdatmd::es2, w3gdatmd::esc, w3gdatmd::esin, w3gdatmd::fachfa, w3gdatmd::flck, w3gdatmd::flcth, w3idatmd::flcur, w3odatmd::iaproc, w3adatmd::is0, w3adatmd::is2, w3gdatmd::mapwn, w3odatmd::ndst, w3gdatmd::nk, w3gdatmd::nspec, w3gdatmd::nth, w3gdatmd::sig, w3servmd::strace(), and constants::tpi.

Referenced by w3wavemd::w3wave().

w3map1()

subroutine w3pro1md::w3map1

(

integer, dimension(ny*nx), intent(in)

MAPSTA

)

Generate 'map' arrays for the first order upstream scheme.

Parameters

MAPSTA

Status map

Author

H. L. Tolman

Date

06-Dec-2010

Definition at line 108 of file w3pro1md.F90.

    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III           NOAA/NCEP |
    !/                  |           H. L. Tolman            |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         06-Dec-2010 |
    !/                  +-----------------------------------+
    !/
    !/    19-Dec-1996 : Final FORTRAN 77                    ( version 1.18 )
    !/    14-Dec-1999 : Upgrade to FORTRAN 90               ( version 2.00 )
    !/    20-Dec-2004 : Multiple grid version.              ( version 3.06 )
    !/    10-Jan-2007 : Clean-up FACVX/Y compute.           ( version 3.10 )
    !/    06-Dec-2010 : Change from GLOBAL (logical) to ICLOSE (integer) to
    !/                  specify index closure for a grid.   ( version 3.14 )
    !/                  (T. J. Campbell, NRL)
    !/
    !  1. Purpose :
    !
    !     Generate 'map' arrays for the first order upstream scheme.
    !
    !  2. Method :
    !
    !     See section 3.
    !
    !  3. Parameters :
    !
    !     Parameter list
    !     ----------------------------------------------------------------
    !       MAPSTA  I.A.   I   Status map.
    !     ----------------------------------------------------------------
    !
    !  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.   Initialize arrays.
    !      2.   Fill arrays.
    !      3.   Invert arrays.
    !     ------------------------------------------------------
    !
    !  9. Switches :
    !
    !     !/S   Enable subroutine tracing.
    !
    ! 10. Source code :
    !
    !/ ------------------------------------------------------------------- /
    USE w3gdatmd, ONLY: nth, nspec, nx, ny, iclose,                 &
         iclose_none, iclose_smpl, iclose_trpl
    USE w3adatmd, ONLY: is0, is2, facvx, facvy
    USE w3odatmd, ONLY: ndse, iaproc, naperr
    USE w3servmd, ONLY: extcde
#ifdef W3_S
    USE w3servmd, ONLY: strace
#endif
    !/
    IMPLICIT NONE
    !/
    !/ ------------------------------------------------------------------- /
    !/ Parameter list
    !/
    INTEGER, INTENT(IN)     :: MAPSTA(NY*NX)
    !/
    !/ ------------------------------------------------------------------- /
    !/ Local parameters
    !/
    INTEGER                 :: IX, IY, IXY, ISP, IXNEXT
#ifdef W3_S
    INTEGER, SAVE           :: IENT = 0
#endif
    !/
    !/ ------------------------------------------------------------------- /
    !/
#ifdef W3_S
    CALL strace (ient, 'W3MAP1')
#endif
    !
 
    ! 1.  Initialize x-y arrays ------------------------------------------ *
    !
    facvx = 0.
    facvy = 0.
    !
    ! 2.  Fill x-y arrays ------------------------------------------------ *
    !
    !.....FACVY
    DO ix=1, nx
      DO iy=1, ny-1
        ixy    = iy +(ix-1)*ny
        IF ( mapsta( ixy ) .NE. 0 ) facvy(ixy) = facvy(ixy) + 1.
        !.........next point : j+1 : increment IXY by 1
        IF ( mapsta(ixy+1) .NE. 0 ) facvy(ixy) = facvy(ixy) + 1.
      END DO
    END DO
    !
    !.....FACVY for IY=NY
    IF ( iclose.EQ.iclose_trpl ) THEN
      iy=ny
      DO ix=1, nx
        ixy    = iy +(ix-1)*ny
        IF ( mapsta( ixy ) .NE. 0 ) facvy(ixy) = facvy(ixy) + 1.
        !...........next point: j+1: tripole: j==>j+1==>j and i==>ni-i+1
        ixnext=nx-ix+1
        ixy    = iy +(ixnext-1)*ny
        IF ( mapsta( ixy ) .NE. 0 ) facvy(ixy) = facvy(ixy) + 1.
      END DO
      !BGR: Adding the following lines to compute FACVX over all
      !      IX for IY=NY (this allows along-seam propagation).
      !      Located here since already inside "TRPL" if-block.
      !{
      DO ix=1, nx-1
        ixy    = iy +(ix-1)*ny
        IF ( mapsta( ixy ) .NE. 0 ) facvx(ixy) = facvx(ixy) + 1.
        IF ( mapsta(ixy+ny) .NE. 0 ) facvx(ixy) = facvx(ixy) + 1.
      END DO
      !}
    END IF
    !
    !.....FACVX
    DO ix=1, nx-1
      DO iy=2, ny-1
        ixy    = iy +(ix-1)*ny
        IF ( mapsta( ixy  ) .NE. 0 ) facvx(ixy) = facvx(ixy) + 1.
        !.........next point : i+1 : increment IXY by NY
        IF ( mapsta(ixy+ny) .NE. 0 ) facvx(ixy) = facvx(ixy) + 1.
      END DO
    END DO
    !
    !.....FACVX for IX=NX
    IF ( iclose.NE.iclose_none ) THEN
      DO iy=2, ny-1
        ixy    = iy +(nx-1)*ny
        IF ( mapsta(ixy) .NE. 0 ) facvx(ixy) = facvx(ixy) + 1.
        !...........next point : i+1 : increment IXY by NY
        !...........IXY+NY=IY+(IX-1)*NY+NY = IY+IX*NY = IY+NX*NY ==> wrap to IY
        IF ( mapsta(iy ) .NE. 0 ) facvx(ixy) = facvx(ixy) + 1.
      END DO
    END IF
    !
    ! 3.  Invert x-y arrays ---------------------------------------------- *
    !
    DO ixy=1, nx*ny
      IF ( facvx(ixy) .NE. 0. ) facvx(ixy) = 1. / facvx(ixy)
      IF ( facvy(ixy) .NE. 0. ) facvy(ixy) = 1. / facvy(ixy)
    END DO
    !
    ! 4.  Fill theta arrays ---------------------------------------------- *
    !
    DO isp=1, nspec
      is2(isp) = isp + 1
      is0(isp) = isp - 1
    END DO
    !
    DO isp=nth, nspec, nth
      is2(isp) = is2(isp) - nth
    END DO
    !
    DO isp=1, nspec, nth
      is0(isp) = is0(isp) + nth
    END DO
    !
    RETURN
    !/
    !/ End of W3MAP1 ----------------------------------------------------- /
    !/

References w3servmd::extcde(), w3adatmd::facvx, w3adatmd::facvy, w3odatmd::iaproc, w3gdatmd::iclose, w3gdatmd::iclose_none, w3gdatmd::iclose_smpl, w3gdatmd::iclose_trpl, w3adatmd::is0, w3adatmd::is2, w3odatmd::naperr, w3odatmd::ndse, w3gdatmd::nspec, w3gdatmd::nth, w3gdatmd::nx, w3gdatmd::ny, and w3servmd::strace().

Referenced by w3wavemd::w3wave().

w3xyp1()

subroutine w3pro1md::w3xyp1	(	integer, intent(in)	ISP,
		real, intent(in)	DTG,
		integer, dimension(ny*nx), intent(in)	MAPSTA,
		real, dimension(1-ny:ny*(nx+2)), intent(inout)	FIELD,
		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,out]	FIELD	Wave action spectral densities on full grid.
[in]	VGX	Speed of grid.
[in]	VGY	Speed of grid.

Author

H. L. Tolman

Date

29-May-2014

Definition at line 303 of file w3pro1md.F90.

    !/
    !/                  +-----------------------------------+
    !/                  | WAVEWATCH III           NOAA/NCEP |
    !/                  |           H. L. Tolman            |
    !/                  |                        FORTRAN 90 |
    !/                  | Last update :         29-May-2014 |
    !/                  +-----------------------------------+
    !/
    !/    07-Jul-1998 : Final FORTRAN 77                    ( version 1.18 )
    !/    14-Dec-1999 : Upgrade to FORTRAN 90               ( version 2.00 )
    !/    28-Mar-2001 : Partial time step bug fix.          ( version 2.10 )
    !/    02-Apr-2001 : Sub-grid obstructions.              ( version 2.10 )
    !/    26-Dec-2002 : Moving grid version.                ( version 3.02 )
    !/    20-Dec-2004 : Multiple grid version.              ( version 3.06 )
    !/    07-Sep-2005 : Improved XY boundary conditions.    ( 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)
    !/    29-May-2014 : Adding OMPH switch.                 ( version 5.02 )
    !/
    !  1. Purpose :
    !
    !     Propagation in physical space for a given spectral component.
    !
    !  2. Method :
    !
    !     First order scheme with flux formulation.
    !     Curvilinear grid implementation: Fluxes are computed in index space
    !       and then transformed back into 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.
    !       FIELD   R.A.  I/O  Wave action spectral densities on full
    !                          grid.
    !       VGX/Y   Real   I   Speed of grid.
    !     ----------------------------------------------------------------
    !
    !     Local variables.
    !     ----------------------------------------------------------------
    !       NTLOC   Int.  Number of local steps.
    !       DTLOC   Real  Local propagation time step.
    !       VCX     R.A.  Propagation velocities in index space.
    !       VCY     R.A.
    !       CXTOT   R.A.  Propagation velocities in physical space.
    !       CYTOT   R.A.
    !       VFLX    R.A.  Discrete fluxes between grid points in index space.
    !       VFLY    R.A.
    !     ----------------------------------------------------------------
    !
    !  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 :
    !
    !     - The local work arrays are initialized on the first entry to
    !       the routine.
    !     - 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.
    !     - Standard VCB calculation for Y is:
    !               VCB       = FACVY(IXY) * ( VCY2D(IY,IX) + VCY2D(IY+1,IX) )
    !       This is to calculate the flux VCY(IY+0.5). For the tripole grid,
    !       we cannot do it this way, since the sign of VCY flips as we jump
    !       over the seam. If we were to do it this way, VCY(IY) and VCY(IY+1)
    !       are two numbers of similar magnitude and opposite sign, so the
    !       average of the two gives something close to zero, so energy does
    !       not leave via VCY(IY+0.5). One alternative is:
    !               VCB       = VCY2D(IY,IX)
    !       Another alternative is :
    !               VCB       = FACVY(IXY) * ( VCY2D(IY,IX) - VCY2D(IY+1,IX) )
    !       Both appear to give correct results for ww3_tp2.13. We use the
    !       second alternative.
    !
    !  8. Structure :
    !
    !     ---------------------------------------
    !       1.  Preparations
    !         a Set constants
    !         b Initialize arrays
    !       2.  Calculate local discrete fluxes
    !       3.  Calculate propagation fluxes
    !       4.  Propagate
    !       5.  Update boundary conditions
    !     ---------------------------------------
    !
    !  9. Switches :
    !
    !     !/S   Enable subroutine tracing.
    !
    !     !/OMPH  Hybrid OpenMP directives.
    !
    !     !/T   Enable general test output.
    !     !/T1  Test output local fluxes (V)FX-YL.
    !     !/T2  Test output propagation fluxes (V)FLX-Y.
    !     !/T3  Test output propagation.
    !
    ! 10. Source code :
    !
    !/ ------------------------------------------------------------------- /
    USE constants
    !
    USE w3timemd, ONLY: dsec21
    !
    USE w3gdatmd, ONLY: nk, nth, sig, ecos, esin, nx, ny, nsea,     &
         mapsf, dtcfl, iclose, clats, flcx, flcy,    &
         iclose_none, iclose_smpl, iclose_trpl,      &
         flagll, dpdx, dpdy, dqdx, dqdy, gsqrt
    USE w3wdatmd, ONLY: time
    USE w3adatmd, ONLY: cg, cx, cy, atrnx, atrny, facvx, facvy
    USE w3idatmd, ONLY: flcur
    USE w3odatmd, ONLY: ndst, flbpi, nbi, tbpi0, tbpin, isbpi,      &
         bbpi0, bbpin, ndse, iaproc, naperr
    USE w3servmd, ONLY: extcde
#ifdef W3_S
    USE w3servmd, ONLY: strace
#endif
    !/
    IMPLICIT NONE
    !/
    !/ ------------------------------------------------------------------- /
    !/ Parameter list
    !/
    INTEGER, INTENT(IN)     :: ISP, MAPSTA(NY*NX)
    REAL, INTENT(IN)        :: DTG, VGX, VGY
    REAL, INTENT(INOUT)     :: FIELD(1-NY:NY*(NX+2))
    !/
    !/ ------------------------------------------------------------------ /
    !/ Local parameters
    !/
    INTEGER                 :: IK, ITH, NTLOC, ITLOC, ISEA, IXY,    &
         IY0, IX, IY, JXN, JXP, JYN, JYP,     &
         IBI, NYMAX
#ifdef W3_T3
    INTEGER                 ::  IXF, IYF
#endif
#ifdef W3_S
    INTEGER, SAVE           :: IENT = 0
#endif
    REAL                    :: CG0, CGL, CGA, CC, CGN
    REAL                    :: DTLOC,DTRAD, VCB
    REAL                    :: RD1, RD2
    REAL                    :: CP, CQ
#ifdef W3_T3
    REAL                    :: AOLD
#endif
    !/
    !/ Automatic work arrays
    !/
    REAL                    :: CXTOT2D(NY,NX)
    REAL                    :: CYTOT2D(NY,NX)
    REAL                    :: FLD2D(NY+1,NX+1)
    REAL                    :: VCX2D(NY,NX+1)
    REAL                    :: VCY2D(NY+1,NX)
    REAL                    :: VFLX2D(1:NY,0:NX)
    REAL                    :: VFLY2D(NY,NX)
 
    !/
    !/ ------------------------------------------------------------------- /
    !/
#ifdef W3_S
    CALL strace (ient, 'W3XYP1')
#endif
    !
    ! 1.  Preparations --------------------------------------------------- *
 
    ! 1.a Set constants
    !
    ith    = 1 + mod(isp-1,nth)
    ik     = 1 + (isp-1)/nth
    !
    cg0    = 0.575 * grav / sig(1)
    cgl    = 0.575 * grav / sig(ik)
    !
    IF ( flcur ) THEN
      cga    = sqrt(maxval((cgl*ecos(ith)+cx(1:nsea))**2          &
           +(cgl*esin(ith)+cy(1:nsea))**2))
      cc     = sqrt(maxval(cx(1:nsea)**2+cy(1:nsea)**2))
#ifdef W3_MGP
      cga    = sqrt(maxval((cgl*ecos(ith)+cx(1:nsea)-vgx)**2      &
           +(cgl*esin(ith)+cy(1:nsea)-vgy)**2))
      cc     = sqrt(maxval((cx(1:nsea)-vgx)**2+(cy(1:nsea)-vgy)**2))
#endif
    ELSE
      cga    = cgl
#ifdef W3_MGP
      cga    = sqrt((cgl*ecos(ith)-vgx)**2+(cgl*esin(ith)-vgy)**2)
#endif
      cc     = 0.
    END IF
    !
    cgn    = 0.9999 * max( cga, cc, 0.001*cg0 )
    !
    ntloc  = 1 + int(dtg/(dtcfl*cg0/cgn))
    dtloc  = dtg / real(ntloc)
    dtrad  = dtloc
    IF ( flagll ) dtrad=dtrad/(dera*radius)
 
    !
#ifdef W3_T
    WRITE (ndst,9000) ntloc
    WRITE (ndst,9001) isp, ith, ik
#endif
    !
    ! ====================== Loop partial ================================ *
    !
    DO itloc=1, ntloc
      !
      ! 1.b Initialize arrays
      !
#ifdef W3_T1
      WRITE (ndst,9010) itloc
#endif
      !
      vcx2d = 0.
      vcy2d = 0.
      cxtot2d  = 0.
      cytot2d  = 0.
      fld2d  = 0.
      vflx2d  = 0.
      vfly2d  = 0.
      !
      ! 2.  Calculate field and velocities --------------------------------- *
      !
      !     FIELD = A / CG * CLATS
      !     VCX   = COS*CG / CLATS
      !     VCY   = SIN*CG
      !
#ifdef W3_T1
      WRITE (ndst,9020)
#endif
      !
#ifdef W3_OMPH
      !$OMP PARALLEL DO PRIVATE (ISEA, IXY, IX, IY)
#endif
      !
      DO isea=1, nsea
        ix     = mapsf(isea,1)
        iy     = mapsf(isea,2)
        ixy    = mapsf(isea,3)
 
        fld2d(iy,ix) = field(ixy) / cg(ik,isea) * clats(isea)
 
        cxtot2d(iy,ix) = ecos(ith) * cg(ik,isea) / clats(isea)
        cytot2d(iy,ix) = esin(ith) * cg(ik,isea)
#ifdef W3_MGP
        cxtot2d(iy,ix) = cxtot2d(iy,ix) - vgx/clats(isea)
        cytot2d(iy,ix) = cytot2d(iy,ix) - vgy
#endif
 
#ifdef W3_T1
        WRITE (ndst,9021) isea, ixy, fld2d(iy,ix), &
             cxtot2d(iy,ix), cytot2d(iy,ix)
#endif
      END DO
      !
#ifdef W3_OMPH
      !$OMP END PARALLEL DO
#endif
      !
      IF ( flcur ) THEN
        DO isea=1, nsea
          ix     = mapsf(isea,1)
          iy     = mapsf(isea,2)
 
          cxtot2d(iy,ix) = cxtot2d(iy,ix) + cx(isea)/clats(isea)
          cytot2d(iy,ix) = cytot2d(iy,ix) + cy(isea)
 
        END DO
      END IF
 
      IF ( flcx ) THEN
        DO isea=1, nsea
          ix     = mapsf(isea,1)
          iy     = mapsf(isea,2)
          cp=cxtot2d(iy,ix)*dpdx(iy,ix)+cytot2d(iy,ix)*dpdy(iy,ix)
          vcx2d(iy,ix) = cp*dtrad
        END DO
      ELSE
        vcx2d=0.0
      ENDIF
 
      IF ( flcy ) THEN
        DO isea=1, nsea
          ix     = mapsf(isea,1)
          iy     = mapsf(isea,2)
          cq=cxtot2d(iy,ix)*dqdx(iy,ix)+cytot2d(iy,ix)*dqdy(iy,ix)
          vcy2d(iy,ix) = cq*dtrad
        END DO
      ELSE
        vcy2d=0.0
      ENDIF
 
      ! Transform FIELD to index space, i.e. straightened space
      ! Bugfix: This is now done *before* adding the ghost row, so that ghost
      !   row will be in index space (bug applied only to global, irregular
      !   grids, so it did not apply to any test case that existed w/v4.18)
      fld2d(1:ny,1:nx)=fld2d(1:ny,1:nx)*gsqrt(1:ny,1:nx)
 
      !
      ! Deal with longitude closure by duplicating one row *to the right*
      !   in FIELD/FLD2D, VCX
      IF ( iclose.NE.iclose_none ) THEN
#ifdef W3_T1
        WRITE (ndst,9024)
#endif
        DO iy=1, ny
          fld2d(iy,nx+1)=fld2d(iy,1)
          vcx2d(iy,nx+1)=vcx2d(iy,1)
#ifdef W3_T1
          WRITE (ndst,9025) iy, fld2d(iy,nx+1), vcx2d(iy,nx+1)
#endif
        END DO
      END IF
 
      ! Deal with tripole closure by duplicating one row *at the top*
      !   in FIELD/FLD2D, VCY
      IF ( iclose.EQ.iclose_trpl ) THEN
        DO ix=1,nx
          !...........next point: j+1: tripole: j==>j+1==>j and i==>ni-i+1
          fld2d(ny+1,ix)=fld2d(ny,nx-ix+1)
          vcy2d(ny+1,ix)=vcy2d(ny,nx-ix+1)
        END DO
      END IF
 
      !
      ! 3.  Calculate propagation fluxes ----------------------------------- *
      !
      nymax=ny-1
      IF ( iclose.EQ.iclose_trpl ) nymax=ny
      !
#ifdef W3_OMPH
      !$OMP PARALLEL DO PRIVATE (IX, IY, IXY, VCB)
#endif
      !
      DO ix=1, nx
        DO iy=1, nymax
          ixy    = iy +(ix-1)*ny
          vcb       = facvx(ixy) * ( vcx2d(iy,ix) + vcx2d(iy,ix+1) )
          vflx2d(iy,ix) = max( vcb , 0. ) * fld2d(iy,ix)          &
               + min( vcb , 0. ) * fld2d(iy,ix+1)
        END DO
      END DO
      !
#ifdef W3_OMPH
      !$OMP END PARALLEL DO
#endif
      !
      ! Deal with longitude closure by duplicating one row *to the left*
      !    in VFLX. Note that a similar action is not take for tripole grid,
      !    since tripole seam is only: IY=NY communicating with other points
      !    at IY=NY, not a case of IY=NY communicating with IY=1
      IF ( iclose.NE.iclose_none ) THEN
#ifdef W3_T2
        WRITE (ndst,9032)
#endif
        DO iy=1, ny
          vflx2d(iy,0) = vflx2d(iy,nx)
#ifdef W3_T2
          WRITE (ndst,9033) iy, vflx2d(iy,0)
#endif
        END DO
      END IF
      !
#ifdef W3_OMPH
      !$OMP PARALLEL DO PRIVATE (IX, IY, IXY, VCB)
#endif
      !
      DO ix=1, nx
        DO iy=1, ny-1
          ixy    = iy +(ix-1)*ny
          vcb       = facvy(ixy) * ( vcy2d(iy,ix) + vcy2d(iy+1,ix) )
          vfly2d(iy,ix) = max( vcb , 0. ) * fld2d(iy,ix)          &
               + min( vcb , 0. ) * fld2d(iy+1,ix)
        END DO
      END DO
      !
#ifdef W3_OMPH
      !$OMP END PARALLEL DO
#endif
      !
 
      ! For tripole grid, include IY=NY in calculation. VCB is handled
      !    differently. See notes in Section "7. Remarks" above.
      IF ( iclose.EQ.iclose_trpl ) THEN
        iy=ny
        !
#ifdef W3_OMPH
        !$OMP PARALLEL DO PRIVATE (IXY, VCB, IX)
#endif
        !
        DO ix=1, nx
          ixy    = iy +(ix-1)*ny
          vcb       = facvy(ixy) * ( vcy2d(iy,ix) - vcy2d(iy+1,ix) )
          vfly2d(iy,ix) = max( vcb , 0. ) * fld2d(iy,ix)          &
               + min( vcb , 0. ) * fld2d(iy+1,ix)
        END DO
        !
#ifdef W3_OMPH
        !$OMP END PARALLEL DO
#endif
        !
      END IF
 
      ! 4.  Propagate ------------------------------------------------------ *
      !
#ifdef W3_T3
      WRITE (ndst,9040)
#endif
      !
#ifdef W3_OMPH
      !$OMP PARALLEL DO PRIVATE (ISEA, IXY, JXN, JXP, JYN, JYP, IX, IY)
#endif
      !
      DO isea=1, nsea
        !
        ix    = mapsf(isea,1)
        iy    = mapsf(isea,2)
        ixy   = mapsf(isea,3)
 
#ifdef W3_T3
        aold   = fld2d(iy,ix) * cg(ik,isea) / clats(isea)
#endif
        !
        IF (mapsta(ixy).EQ.1) THEN
          !
          IF ( vflx2d(iy,ix-1) .GT. 0. ) THEN
            jxn   = -1
          ELSE
            jxn   =  0
          END IF
          IF ( vflx2d(iy,ix) .LT. 0. ) THEN
            jxp   =  1
          ELSE
            jxp   =  0
          END IF
          IF ( vfly2d(iy-1,ix) .GT. 0. ) THEN
            jyn   = -1
          ELSE
            jyn   =  0
          END IF
          IF ( vfly2d(iy,ix) .LT. 0. ) THEN
            jyp   =  1
          ELSE
            jyp   =  0
          END IF
          !
          fld2d(iy,ix) =   fld2d(iy,ix)                            &
               + atrnx(ixy,jxn) * vflx2d(iy,ix-1)   &
               - atrnx(ixy,jxp) * vflx2d(iy,ix)     &
               + atrny(ixy,jyn) * vfly2d(iy-1,ix)   &
               - atrny(ixy,jyp) * vfly2d(iy,ix)
 
#ifdef W3_T3
          WRITE (ndst,9041) isea, ixy, ixy-ny, ixy-1,         &
               vflx2d(iy,ix-1), vflx2d(iy,ix), vfly2d(iy-1,ix), &
               vfly2d(iy,ix) , cg(ik,isea)/clats(isea),aold,    &
               fld2d(iy,ix)
#endif
          !
          !
#ifdef W3_T3
          WRITE (ndst,9042) isea, mapsta(ixy), aold,fld2d(iy,ix)
#endif
          !
        END IF !  IF (MAPSTA(IXY).EQ.1) THEN
 
        fld2d(iy,ix) = cg(ik,isea) / clats(isea) * fld2d(iy,ix)
 
      END DO !  DO ISEA=1, NSEA
      !
#ifdef W3_OMPH
      !$OMP END PARALLEL DO
#endif
      !
 
      ! Transform FIELD back to physical space, i.e. may be curvilinear
      fld2d(1:ny,1:nx)=fld2d(1:ny,1:nx)/gsqrt(1:ny,1:nx)
      !
      ! 5.  Update boundary conditions ------------------------------------- *
      !
      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
          ix    = mapsf(isbpi(ibi),1)
          iy    = mapsf(isbpi(ibi),2)
          fld2d(iy,ix) = rd1*bbpi0(isp,ibi) + rd2*bbpin(isp,ibi)
        END DO
      END IF
      !
      ! 6.  Put back in 1d shape ------------------------------------------- *
      !
      DO isea=1, nsea
        ix     = mapsf(isea,1)
        iy     = mapsf(isea,2)
        ixy    = mapsf(isea,3)
        field(ixy) = fld2d(iy,ix)
      END DO
      !
      ! ... End of partial time step loop
      !
    END DO !   DO ITLOC=1, NTLOC
    !
    RETURN
    !
    ! Formats
    !
#ifdef W3_T
9000 FORMAT (' TEST W3XYP1 : NTLOC :',i4)
9001 FORMAT (' TEST W3XYP1 : ISP, ITH, IK :',i8,2i4)
#endif
    !
#ifdef W3_T1
9010 FORMAT (' TEST W3XYP1 : INIT. VFX-YL, ITLOC =',i3)
9020 FORMAT (' TEST W3XYP1 : ISEA, IXY, FIELD, VCX, VCY')
9021 FORMAT ('           ',2i8,3e12.4)
9024 FORMAT (' TEST W3XYP1 : GLOBAL CLOSURE: IY, FIELD, VCX ')
9025 FORMAT ('               ',i4,2e12.4)
#endif
    !
#ifdef W3_T2
9032 FORMAT (' TEST W3XYP1 : CLOSE. : IY, VFLX')
9033 FORMAT ('            ',i4,e12.4)
#endif
    !
#ifdef W3_T3
9040 FORMAT (' TEST W3XYP1 : PROPAGATION '/                      &
         '      ISEA, IXY(3), , FLX(2), FLY(2), FAC, A(2)')
9041 FORMAT (2x,4i5,1x,4e10.3,1x,e10.3,1x,2e10.3)
9042 FORMAT (2x,i5,'( MAP = ',i2,' )',56x,2e10.3)
#endif
    !/
    !/ End of W3XYP1 ----------------------------------------------------- /
    !/

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::ecos, w3gdatmd::esin, w3servmd::extcde(), w3adatmd::facvx, w3adatmd::facvy, 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, w3gdatmd::mapsf, w3odatmd::naperr, w3odatmd::nbi, w3odatmd::ndse, w3odatmd::ndst, w3gdatmd::nk, w3gdatmd::nsea, w3gdatmd::nth, w3gdatmd::nx, w3gdatmd::ny, constants::radius, w3gdatmd::sig, w3servmd::strace(), w3odatmd::tbpi0, w3odatmd::tbpin, and w3wdatmd::time.

Referenced by w3wavemd::w3wave().