#include "w3macros.h"
!/ ------------------------------------------------------------------- /
MODULE W3PRO3MD
!/
!/ +-----------------------------------+
!/ | WAVEWATCH III NOAA/NCEP |
!/ | H. L. Tolman |
!/ | FORTRAN 90 |
!/ | Last update : 27-May-2014 |
!/ +-----------------------------------+
!/
!/ 27-Feb-2000 : Origination. ( version 2.08 )
!/ 17-Sep-2000 : Clean-up. ( version 2.13 )
!/ 10-Dec-2001 : Sub-grid obstructions. ( version 2.14 )
!/ 16-Oct-2002 : Change INTENT for ATRN in W3XYP3. ( version 3.00 )
!/ 26-Dec-2002 : Moving grid version. ( version 3.02 )
!/ 01-Aug-2003 : Moving grid GSE correction. ( version 3.03 )
!/ 17-Dec-2004 : Multiple grid version. ( version 3.06 )
!/ 07-Sep-2005 : Upgrade XY boundary conditions. ( version 3.08 )
!/ 09-Nov-2005 : Removing soft boundary option. ( version 3.08 )
!/ 05-Mar-2008 : Added NEC sxf90 compiler directives.
!/ (Chris Bunney, UK Met Office) ( version 3.13 )
!/ 01-Apr-2008 : Bug fix W3MAP3 MAPSTA range check. ( version 3.13 )
!/ 29-May-2009 : Preparing distribution version. ( version 3.14 )
!/ 30-Oct-2009 : Implement curvilinear grid type. ( version 3.14 )
!/ (W. E. Rogers & T. J. Campbell, NRL)
!/ 17-Aug-2010 : Add test output W3XYP3. ( version 3.14.5 )
!/ 06-Dec-2010 : Change from GLOBAL (logical) to ICLOSE (integer) to
!/ specify index closure for a grid. ( version 3.14 )
!/ (T. J. Campbell, NRL)
!/ 26-Dec-2012 : More initializations. ( version 4.11 )
!/ 01-Jul-2013 : Adding UQ and UNO switches to chose between third
!/ and second order schemes. ( version 4.12 )
!/ 12-Sep-2013 : Add documentation for global clos. ( version 4.12 )
!/ 27-May-2014 : Adding OMPH switch. ( version 5.02 )
!/
!/ Copyright 2009-2014 National Weather Service (NWS),
!/ National Oceanic and Atmospheric Administration. All rights
!/ reserved. WAVEWATCH III is a trademark of the NWS.
!/ No unauthorized use without permission.
!/
! 1. Purpose :
!
! Bundles routines for third order propagation scheme in single
! module.
!
! 2. Variables and types :
!
! Name Type Scope Description
! ----------------------------------------------------------------
! TRNMIN R.P. Private Minimum transparancy for local
! switching off of averaging.
! ----------------------------------------------------------------
!
! Also work arrays for W3KTP3 (private).
!
! 3. Subroutines and functions :
!
! Name Type Scope Description
! ----------------------------------------------------------------
! W3MAP3 Subr. Public Set up auxiliary maps.
! W3MAPT Subr. Public Set up transparency map for GSE.
! W3XYP3 Subr. Public Third order spatial propagation.
! W3KTP3 Subr. Public Third order spectral propagation.
! ----------------------------------------------------------------
!
! 4. Subroutines and functions used :
!
! Name Type Module Description
! ----------------------------------------------------------------
! STRACE Subr. W3SERVMD Subroutine tracing.
! W3QCK1 Subr. W3UQCKMD Regular grid UQ scheme.
! W3QCK2 Subr. Id. Irregular grid UQ scheme.
! W3QCK3 Subr. Id. Regular grid UQ scheme + obstructions.
! W3UNO2 Subr. W3UNO2MD UNO2 scheme for irregular grid.
! W3UNO2r Subr. Id. UNO2 scheme reduced to regular grid.
! W3UNO2s Subr. Id. UNO2 regular grid with subgrid
! obstruction.
! ----------------------------------------------------------------
!
! 5. Remarks :
!
! - The averaging is not performed around semi-transparent grid
! points to avoid that leaking through barriers occurs.
!
! 6. Switches :
!
! !/UQ 3rd order UQ propagation scheme.
! !/UNO 2nd order UNO propagation scheme.
!
! !/MGP Moving grid corrections.
! !/MGG Moving grid corrections.
!
! !/OMPH Hybrid OpenMP directives.
!
! !/S Enable subroutine tracing.
! !/Tn Enable test output.
!
! 7. Source code :
!
!/ ------------------------------------------------------------------- /
!/
!/ Public variables
!/
PUBLIC
!/
!/ Private data
!/
REAL, PRIVATE, PARAMETER:: TRNMIN = 0.95
!/
CONTAINS
!/ ------------------------------------------------------------------- /
SUBROUTINE W3MAP3
!/
!/ +-----------------------------------+
!/ | WAVEWATCH III NOAA/NCEP |
!/ | H. L. Tolman |
!/ | FORTRAN 90 |
!/ | Last update : 01-Apr-2008 |
!/ +-----------------------------------+
!/
!/ 27-Feb-2000 : Origination. ( version 2.08 )
!/ 10-Dec-2001 : Sub-grid obstructions. ( version 2.14 )
!/ (array allocation only.)
!/ 17-Dec-2004 : Multiple grid version. ( version 3.06 )
!/ 09-Nov-2005 : Removing soft boundary option. ( version 3.08 )
!/ 01-Apr-2008 : Bug fix sec. 4 MAPSTA range check. ( version 3.13 )
!/
! 1. Purpose :
!
! Generate 'map' arrays for the ULTIMATE QUICKEST scheme.
!
! 2. Method :
!
! MAPX2, MAPY2, MAPTH2 and MAPWN2 contain consecutive 1-D spatial
! grid counters (e.g., IXY = (IX-1)*MY + IY). The arrays are
! devided in three parts. For MAPX2, these ranges are :
!
! 1 - NMX0 Counters IXY for which grid point (IX,IY) and
! (IX+1,IY) both are active grid points.
! NMX0+1 - NMX1 Id. only (IX,IY) active.
! NMX1+1 - NMX2 Id. only (IX+1,IY) active.
!
! The array MAPY2 has a similar layout varying IY instead of IX.
!
! 3. Parameters :
!
! Parameter list
! ----------------------------------------------------------------
! ----------------------------------------------------------------
!
! 4. Subroutines used :
!
! See module documentation.
!
! 5. Called by :
!
! Name Type Module Description
! ----------------------------------------------------------------
! W3WAVE Subr. W3WAVEMD Wave model routine.
! ----------------------------------------------------------------
!
! 6. Error messages :
!
! 7. Remarks :
!
! 8. Structure :
!
! ------------------------------------------------------
! 1. Map MAPX2
! a Range 1 to NMX0
! b Range NMX0+1 to NMX1
! c Range NMX1+1 to NMX2
! 2. Map MAPY2
! a Range 1 to NMY0
! b Range NMY0+1 to NMY1
! c Range NMY1+1 to NMY2
! 3. Map MAPAXY
! 4. Map MAPCXY
! 5. Maps for intra-spectral propagation
! a MAPTH2, MAPATK
! b MAPWN2
! ------------------------------------------------------
!
! 9. Switches :
!
! !/S Enable subroutine tracing.
! !/T Enable test output.
!
! 10. Source code :
!/ ------------------------------------------------------------------- /
USE W3GDATMD, ONLY: NK, NTH, NSPEC, NX, NY, NSEA, MAPSTA, MAPSF,&
GTYPE
USE W3ADATMD, ONLY: NMX0, NMX1, NMX2, NMY0, NMY1, NMY2, NACT, &
NCENT, MAPX2, MAPY2, MAPAXY, MAPCXY, &
MAPTH2, MAPWN2
USE W3ODATMD, ONLY: NDST
#ifdef W3_S
USE W3SERVMD, ONLY: STRACE
#endif
!/
IMPLICIT NONE
!/
!/ ------------------------------------------------------------------- /
!/ Parameter list
!/
!/ ------------------------------------------------------------------- /
!/ Local parameters
!/
INTEGER :: IX, IY, IXY0, IX2, IY2, IX0, IY0, &
ISEA, IK, ITH, ISP, ISP0, ISP2, NCENTC
#ifdef W3_S
INTEGER, SAVE :: IENT = 0
#endif
#ifdef W3_T
INTEGER :: MAPTXY(NY,NX), I, IXY
INTEGER :: MAPTST(NK+2,NTH)
#endif
!/
!/ ------------------------------------------------------------------- /
!/
#ifdef W3_S
CALL STRACE (IENT, 'W3MAP3')
#endif
!
IF (GTYPE .LT. 3) THEN
! 1. Map MAPX2 ------------------------------------------------------ *
! 1.a Range 1 to NMX0
!
#ifdef W3_T
MAPTXY = 0.
#endif
!
NMX0 = 0
DO IX=1, NX
IXY0 = (IX-1)*NY
IX2 = 1 + MOD(IX,NX)
DO IY=2, NY-1
IF ( MAPSTA(IY,IX).EQ.1 .AND. MAPSTA(IY,IX2).EQ.1 ) THEN
NMX0 = NMX0 + 1
MAPX2(NMX0) = IXY0 + IY
#ifdef W3_T
MAPTXY(IY,IX) = MAPTXY(IY,IX) + 1
#endif
END IF
END DO
END DO
!
! 1.b Range NMX0+1 to NMX1
!
NMX1 = NMX0
DO IX=1, NX
IXY0 = (IX-1)*NY
IX2 = 1 + MOD(IX,NX)
DO IY=2, NY-1
IF ( MAPSTA(IY,IX).EQ.1 .AND. MAPSTA(IY,IX2).NE.1 ) THEN
NMX1 = NMX1 + 1
MAPX2(NMX1) = IXY0 + IY
#ifdef W3_T
MAPTXY(IY,IX) = MAPTXY(IY,IX) + 2
#endif
END IF
END DO
END DO
!
! 1.c Range NMX1+1 to NMX2
!
NMX2 = NMX1
DO IX=1, NX
IXY0 = (IX-1)*NY
IX2 = 1 + MOD(IX,NX)
DO IY=2, NY-1
IF ( MAPSTA(IY,IX).NE.1 .AND. MAPSTA(IY,IX2).EQ.1 ) THEN
NMX2 = NMX2 + 1
MAPX2(NMX2) = IXY0 + IY
#ifdef W3_T
MAPTXY(IY,IX) = MAPTXY(IY,IX) + 4
#endif
END IF
END DO
END DO
!
#ifdef W3_T
WRITE (NDST,9000) 'MAPX2', NMX0, NMX1-NMX0, &
NMX2-NMX1, NMX2
DO IY=NY, 1, -1
WRITE (NDST,9001) (MAPTXY(IY,IX),IX=1, NX)
END DO
#endif
!
! 2. Map MAPY2 ------------------------------------------------------ *
! 2.a Range 1 to NMY0
!
#ifdef W3_T
MAPTXY = 0.
#endif
!
NMY0 = 0
DO IX=1, NX
IXY0 = (IX-1)*NY
DO IY=1, NY-1
IY2 = IY + 1
IF ( MAPSTA(IY,IX).EQ.1 .AND. MAPSTA(IY2,IX).EQ.1 ) THEN
NMY0 = NMY0 + 1
MAPY2(NMY0) = IXY0 + IY
#ifdef W3_T
MAPTXY(IY,IX) = MAPTXY(IY,IX) + 1
#endif
END IF
END DO
END DO
!
! 2.b Range NMY0+1 to NMY1
!
NMY1 = NMY0
DO IX=1, NX
IXY0 = (IX-1)*NY
DO IY=1, NY-1
IY2 = IY + 1
IF ( MAPSTA(IY,IX).EQ.1 .AND. MAPSTA(IY2,IX).NE.1 ) THEN
NMY1 = NMY1 + 1
MAPY2(NMY1) = IXY0 + IY
#ifdef W3_T
MAPTXY(IY,IX) = MAPTXY(IY,IX) + 2
#endif
END IF
END DO
END DO
!
! 2.c Range NMY1+1 to NMY2
!
NMY2 = NMY1
DO IX=1, NX
IXY0 = (IX-1)*NY
DO IY=1, NY-1
IY2 = IY + 1
IF ( MAPSTA(IY,IX).NE.1 .AND. MAPSTA(IY2,IX).EQ.1 ) THEN
NMY2 = NMY2 + 1
MAPY2(NMY2) = IXY0 + IY
#ifdef W3_T
MAPTXY(IY,IX) = MAPTXY(IY,IX) + 4
#endif
END IF
END DO
END DO
!
#ifdef W3_T
WRITE (NDST,9000) 'MAPY2', NMY0, NMY1-NMY0, &
NMY2-NMY1, NMY2
DO IY=NY, 1, -1
WRITE (NDST,9001) (MAPTXY(IY,IX),IX=1, NX)
END DO
#endif
!
! 3. Map MAPAXY ----------------------------------------------------- *
!
NACT = 0
DO IX=1, NX
IY0 = (IX-1)*NY
DO IY=2, NY-1
IF ( MAPSTA(IY,IX).EQ.1 ) THEN
NACT = NACT + 1
MAPAXY(NACT) = IY0 + IY
END IF
END DO
END DO
!
! 4. Map MAPCXY ----------------------------------------------------- *
!
NCENT = 0
NCENTC = NSEA
MAPCXY = 0
!
DO ISEA=1, NSEA
IX = MAPSF(ISEA,1)
IX0 = IX-1
IX2 = IX+1
IY = MAPSF(ISEA,2)
IY0 = IY-1
IY2 = IY+1
IF ( IX .EQ. NX ) IX2 = 1
IF ( IX .EQ. 1 ) IX0 = NX
IF ( MAPSTA(IY,IX).EQ.2 .OR. MAPSTA(IY,IX).LT.0 ) THEN
MAPCXY(NCENTC) = ISEA
NCENTC = NCENTC - 1
ELSE IF ( MAPSTA(IY0,IX0).GE.1 .AND. &
MAPSTA(IY0,IX ).GE.1 .AND. &
MAPSTA(IY0,IX2).GE.1 .AND. &
MAPSTA(IY ,IX0).GE.1 .AND. &
MAPSTA(IY ,IX2).GE.1 .AND. &
MAPSTA(IY2,IX0).GE.1 .AND. &
MAPSTA(IY2,IX ).GE.1 .AND. &
MAPSTA(IY2,IX2).GE.1 ) THEN
NCENT = NCENT + 1
MAPCXY(NCENT) = ISEA
ELSE
MAPCXY(NCENTC) = ISEA
NCENTC = NCENTC - 1
END IF
END DO
END IF
!
! 5. Maps for intra-spectral propagation ---------------------------- *
!
IF ( MAPTH2(1) .NE. 0 ) RETURN
!
#ifdef W3_T
MAPTST = 0
#endif
!
! 5.a MAPTH2 and MAPBTK
!
DO IK=1, NK
DO ITH=1, NTH
ISP = ITH + (IK-1)*NTH
ISP2 = (IK+1) + (ITH-1)*(NK+2)
MAPTH2(ISP) = ISP2
#ifdef W3_T
MAPTST(IK+1,ITH) = MAPTST(IK+1,ITH) + 1
#endif
END DO
END DO
!
#ifdef W3_T
WRITE (NDST,9000) 'MAPTH2', ISP, 0, 0, ISP
DO IK=NK+2, 1, -1
WRITE (NDST,9001) (MAPTST(IK,ITH),ITH=1, NTH)
END DO
MAPTST = 0
#endif
!
! 5.b MAPWN2
!
ISP0 = 0
DO IK=1, NK-1
DO ITH=1, NTH
ISP0 = ISP0 + 1
ISP2 = (IK+1) + (ITH-1)*(NK+2)
MAPWN2(ISP0) = ISP2
#ifdef W3_T
MAPTST(IK+1,ITH) = MAPTST(IK+1,ITH) + 1
#endif
END DO
END DO
!
DO ITH=1, NTH
ISP0 = ISP0 + 1
ISP2 = NK+1 + (ITH-1)*(NK+2)
MAPWN2(ISP0) = ISP2
#ifdef W3_T
MAPTST(NK+1,ITH) = MAPTST(NK+1,ITH) + 2
#endif
END DO
!
DO ITH=1, NTH
ISP0 = ISP0 + 1
ISP2 = 1 + (ITH-1)*(NK+2)
MAPWN2(ISP0) = ISP2
#ifdef W3_T
MAPTST(1,ITH) = MAPTST(1,ITH) + 4
#endif
END DO
!
#ifdef W3_T
WRITE (NDST,9000) 'MAPWN2', NSPEC-NTH, NTH, NTH, NSPEC+NTH
DO IK=NK+2, 1, -1
WRITE (NDST,9001) (MAPTST(IK,ITH),ITH=1, NTH)
END DO
#endif
!
RETURN
!
! Formats
!
#ifdef W3_T
9000 FORMAT (/' TEST W3MAP3 : TEST MAP FOR PROPAGATION'/ &
' MAP : ',A/ &
' CENTRAL : ',I6/ &
' ABOVE : ',I6/ &
' BELOW : ',I6/ &
' TOTAL : ',I6/)
9001 FORMAT (1X,130I1)
9010 FORMAT (' TEST W3MAP3 : COMPOSITE MAPS TH2, WN2 AND BTK')
9011 FORMAT (2X,60I2)
#endif
!/
!/ End of W3MAP3 ----------------------------------------------------- /
!/
END SUBROUTINE W3MAP3
!/ ------------------------------------------------------------------- /
SUBROUTINE W3MAPT
!/
!/ +-----------------------------------+
!/ | WAVEWATCH III NOAA/NCEP |
!/ | H. L. Tolman |
!/ | FORTRAN 90 |
!/ | Last update : 17-Dec-2004 |
!/ +-----------------------------------+
!/
!/ 10-Dec-2001 : Origination. ( version 2.14 )
!/ 17-Dec-2004 : Multiple grid version. ( version 3.06 )
!/
! 1. Purpose :
!
! Generate 'map' arrays for the ULTIMATE QUICKEST scheme to combine
! GSE alleviation with obstructions.
!
! 2. Method :
!
! 3. Parameters :
!
! Parameter list
! ----------------------------------------------------------------
! ----------------------------------------------------------------
!
! 4. Subroutines used :
!
! See module documentation.
!
! 5. Called by :
!
! Name Type Module Description
! ----------------------------------------------------------------
! W3WAVE Subr. W3WAVEMD Wave model routine.
! ----------------------------------------------------------------
!
! 6. Error messages :
!
! 7. Remarks :
!
! 8. Structure :
!
! See source code.
!
! 9. Switches :
!
! !/S Enable subroutine tracing.
!
! 10. Source code :
!/ ------------------------------------------------------------------- /
USE W3GDATMD, ONLY: NX, NY, NSEA, MAPSF
USE W3ADATMD, ONLY: ATRNX, ATRNY, MAPTRN
#ifdef W3_S
USE W3SERVMD, ONLY: STRACE
#endif
!/
IMPLICIT NONE
!/
!/ ------------------------------------------------------------------- /
!/ Parameter list
!/
!/ ------------------------------------------------------------------- /
!/ Local parameters
!/
INTEGER :: ISEA, IXY
#ifdef W3_S
INTEGER, SAVE :: IENT = 0
#endif
!/
!/ ------------------------------------------------------------------- /
!/
#ifdef W3_S
CALL STRACE (IENT, 'W3MAPT')
#endif
!
! 1. Map MAPTRN ----------------------------------------------------- *
!
DO ISEA=1, NSEA
IXY = MAPSF(ISEA,3)
!notes: Oct 22 2012: I changed this because it *looks* like a bug.
! I have not confirmed that it is a bug.
! Old code is given (2 lines). Only the first line is
! changed.
!old MAPTRN(IXY) = MIN( ATRNX(IXY,1) ,ATRNY(IXY,-1) , &
!old ATRNY(IXY,1), ATRNY(IXY,-1) ) .LT. TRNMIN
MAPTRN(IXY) = MIN( ATRNX(IXY,1) ,ATRNX(IXY,-1) , &
ATRNY(IXY,1), ATRNY(IXY,-1) ) .LT. TRNMIN
END DO
!
RETURN
!
! Formats
!/
!/ End of W3MAPT ----------------------------------------------------- /
!/
END SUBROUTINE W3MAPT
!/ ------------------------------------------------------------------- /
SUBROUTINE W3XYP3 ( ISP, DTG, MAPSTA, MAPFS, VQ, VGX, VGY )
!/
!/ +-----------------------------------+
!/ | WAVEWATCH III NOAA/NCEP |
!/ | H. L. Tolman |
!/ | FORTRAN 90 |
!/ | Last update : 27-May-2014 |
!/ +-----------------------------------+
!/
!/ 27-Feb-2000 : Origination. ( version 2.08 )
!/ 17-Sep-2000 : Clean-up. ( version 2.13 )
!/ 10-Dec-2001 : Sub-grid obstructions. ( version 2.14 )
!/ 16-Oct-2002 : Change INTENT for ATRNRX/Y. ( version 3.00 )
!/ 26-Dec-2002 : Moving grid version. ( version 3.02 )
!/ 01-Aug-2003 : Moving grid GSE correction. ( version 3.03 )
!/ 17-Dec-2004 : Multiple grid version. ( version 3.06 )
!/ 07-Sep-2005 : Upgrade XY boundary conditions. ( version 3.08 )
!/ 09-Nov-2005 : Removing soft boundary option. ( version 3.08 )
!/ 05-Mar-2008 : Added NEC sxf90 compiler directives.
!/ (Chris Bunney, UK Met Office) ( version 3.13 )
!/ 30-Oct-2009 : Implement curvilinear grid type. ( version 3.14 )
!/ (W. E. Rogers & T. J. Campbell, NRL)
!/ 17-Aug-2010 : Add test output. ( version 3.14.5 )
!/ 30-Oct-2010 : Implement unstructured grid ( version 3.14 )
!/ 06-Dec-2010 : Change from GLOBAL (logical) to ICLOSE (integer) to
!/ specify index closure for a grid. ( version 3.14 )
!/ (T. J. Campbell, NRL)
!/ 01-Jul-2013 : Adding UQ and UNO switches to chose between third
!/ and second order schemes. ( version 4.12 )
!/ 12-Sep-2013 : Add documentation for global clos. ( version 4.12 )
!/ 27-May-2014 : Adding OMPH switch. ( version 5.02 )
!/
! 1. Purpose :
!
! Propagation in phyiscal space for a given spectral component.
!
! 2. Method :
!
! Third-order ULTIMATE QUICKEST scheme with averaging.
! Curvilinear grid implementation: Fluxes are computed in index space
! and then transformed back into physical space. The diffusion term
! is handled in physical space. The averaging scheme is adapted for
! curvilinear grids by applying the appropriate local rotations and
! adjustments to interpolation weights which control the strength of
! the averaging in axial directions.
!
! 3. Parameters :
!
! Parameter list
! ----------------------------------------------------------------
! ISP Int. I Number of spectral bin (IK-1)*NTH+ITH
! DTG Real I Total time step.
! MAPSTA I.A. I Grid point status map.
! MAPFS I.A. I Storage map.
! VQ R.A. I/O Field to propagate.
! VGX/Y Real I Speed of grid.
! ----------------------------------------------------------------
!
! Local variables.
! ----------------------------------------------------------------
! NTLOC Int Number of local time steps.
! DTLOC Real Local propagation time step.
! VCFL0X R.A. Local courant numbers for absolute group vel.
! using local X-grid step.
! VCFL0Y R.A. Id. in Y.
! ----------------------------------------------------------------
!
! 4. Subroutines used :
!
! W3QCK3 Actual propagation algorithm
!
! STRACE Service routine.
!
! 5. Called by :
!
! W3WAVE Wave model routine.
!
! 6. Error messages :
!
! None.
!
! 7. Remarks :
!
! - Note that the ULTIMATE limiter does not guarantee non-zero
! energies.
! - The present scheme shows a strong distorsion when propaga-
! ting a field under an angle with the grid in a truly 2-D
! fashion. Propagation is therefore split along the two
! axes.
! - Two boundary treatments are available. The first uses real
! boundaries in each space. In this case waves will not
! penetrate in narrow straights under an angle with the grid.
! This behavior is improved by using a 'soft' option, in
! which the 'X' or 'Y' sweep allows for energy to go onto
! the land. This improves the above behavior, but implies
! that X-Y connenctions are required in barriers for them
! to become inpenetrable.
! - Note that unlike in W3XYP2, isotropic diffusion is never
! used for growth.
! - Curvilinear grid implementation. Variables FACX, FACY, CCOS, CSIN,
! CCURX, CCURY are not needed and have been removed. FACX is accounted
! for as approriate in this subroutine. FACX is also accounted for in
! the case of .NOT.FLCX. Since FACX is removed, there is now a check for
! .NOT.FLCX in this subroutine. In CFL calcs dx and dy are omitted,
! since dx=dy=1 in index space. Curvilinear grid derivatives
! (DPDY, DQDX, etc.) and metric (GSQRT) are brought in via W3GDATMD.
! - The strength of the averaging scheme is dependent on grid resolution.
! Since grid resolution is non-uniform for curvilinear grids, this means
! that the strength of the averaging is also non-uniform. This may not be
! a desirable effect. A potential future upgrade would be to add an
! additional term/factor that balances the effect of the spatial
! variation of grid resolution.
!
! 8. Structure :
!
! ---------------------------------------------
! 1. Preparations
! a Set constants
! b Initialize arrays
! 2. Prepare arrays
! a Velocities and 'Q'
! 3. Loop over sub-steps
! ----------------------------------------
! a Average
! b Propagate
! c Update boundaries.
! ----------------------------------------
! 4. Store Q field in spectra
! ---------------------------------------------
!
! 9. Switches :
!
! !/S Enable subroutine tracing.
!
! !/OMPH Hybrid OpenMP directives.
!
! !/MGP Moving grid corrections.
! !/MGG Moving grid corrections.
!
! !/T Enable general test output.
! !/T0 Dump of precalcaulted interpolation data.
! !/T1 Dump of input field and fluxes.
! !/T2 Dump of output field (before boundary update).
! !/T3 Dump of output field (final).
!
! 10. Source code :
!
!/ ------------------------------------------------------------------- /
USE CONSTANTS
!
USE W3TIMEMD, ONLY: DSEC21
!
USE W3GDATMD, ONLY: NX, NY, NSEA, MAPSF, DTCFL, CLATS, &
ICLOSE, FLCX, FLCY, NK, NTH, DTH, XFR, &
ICLOSE_NONE, ICLOSE_SMPL, ICLOSE_TRPL, &
ECOS, ESIN, SIG, WDCG, WDTH, PFMOVE, &
FLAGLL, DPDX, DPDY, DQDX, DQDY, GSQRT
USE W3WDATMD, ONLY: TIME
USE W3ADATMD, ONLY: NMX0, NMX1, NMX2, NMY0, NMY1, NMY2, NACT, &
NCENT, MAPX2, MAPY2, MAPAXY, MAPCXY, &
MAPTRN, CG, CX, CY, ATRNX, ATRNY, ITIME
USE W3IDATMD, ONLY: FLCUR
USE W3ODATMD, ONLY: NDSE, NDST, FLBPI, NBI, TBPI0, TBPIN, &
ISBPI, BBPI0, BBPIN, IAPROC, NAPERR
USE W3SERVMD, ONLY: EXTCDE
#ifdef W3_S
USE W3SERVMD, ONLY: STRACE
#endif
#ifdef W3_UQ
USE W3UQCKMD
#endif
#ifdef W3_UNO
USE W3UNO2MD
#endif
!/
IMPLICIT NONE
!/
!/ ------------------------------------------------------------------- /
!/ Parameter list
!/
INTEGER, INTENT(IN) :: ISP, MAPSTA(NY*NX), MAPFS(NY*NX)
REAL, INTENT(IN) :: DTG, VGX, VGY
REAL, INTENT(INOUT) :: VQ(1-NY:NY*(NX+2))
!/
!/ ------------------------------------------------------------------- /
!/ Local parameters
!/
INTEGER :: ITH, IK, NTLOC, ITLOC, ISEA, IXY, IP
INTEGER :: IX, IY, IXC, IYC, IBI
INTEGER :: IIXY1(NSEA), IIXY2(NSEA), &
IIXY3(NSEA), IIXY4(NSEA)
INTEGER :: TTEST(2),DTTST
#ifdef W3_S
INTEGER, SAVE :: IENT = 0
#endif
REAL :: CG0, CGA, CGN, CGX, CGY, CXC, CYC, &
CXMIN, CXMAX, CYMIN, CYMAX
REAL :: CGC, FGSE = 1.
REAL :: FTH, FTHX, FTHY, FCG, FCGX, FCGY
REAL :: DTLOC, DTRAD, &
DXCGN, DYCGN, DXCGS, DYCGS, DXCGC, &
DYCGC
REAL :: RDI1(NSEA), RDI2(NSEA), &
RDI3(NSEA), RDI4(NSEA)
REAL :: TMPX, TMPY, RD1, RD2, RD3, RD4
LOGICAL :: YFIRST
LOGICAL :: GLOBAL
REAL :: CP, CQ
!/
!/ Automatic work arrays
!/
INTEGER :: MAPSTX(1-2*NY:NY*(NX+2))
REAL :: VLCFLX((NX+1)*NY), VLCFLY((NX+1)*NY),&
AQ(1-NY:NY*(NX+2))
REAL :: CXTOT((NX+1)*NY), CYTOT(NX*NY)
!/
!/ ------------------------------------------------------------------- /
!/
#ifdef W3_S
CALL STRACE (IENT, 'W3XYP3')
#endif
!
! 1. Preparations --------------------------------------------------- *
IF ( ICLOSE .EQ. ICLOSE_TRPL ) THEN
!/ ------------------------------------------------------------------- /
IF (IAPROC .EQ. NAPERR) &
WRITE(NDSE,*)'SUBROUTINE W3XYP3 IS NOT YET ADAPTED FOR '// &
'TRIPOLE GRIDS. STOPPING NOW.'
CALL EXTCDE ( 1 )
END IF
! 1.a Set constants
!
GLOBAL = ICLOSE.NE.ICLOSE_NONE
ITH = 1 + MOD(ISP-1,NTH)
IK = 1 + (ISP-1)/NTH
!
CG0 = 0.575 * GRAV / SIG(1)
CGA = 0.575 * GRAV / SIG(IK)
CGX = CGA * ECOS(ITH)
CGY = CGA * ESIN(ITH)
#ifdef W3_MGP
CGX = CGX - VGX
CGY = CGY - VGY
#endif
CGC = SQRT ( CGX**2 + CGY**2 )
#ifdef W3_MGG
FGSE = ( CGA / MAX(0.001*CGA,CGC) )**PFMOVE
#endif
!
IF ( FLCUR ) THEN
CXMIN = MINVAL ( CX(1:NSEA) )
CXMAX = MAXVAL ( CX(1:NSEA) )
CYMIN = MINVAL ( CY(1:NSEA) )
CYMAX = MAXVAL ( CY(1:NSEA) )
IF ( ABS(CGX+CXMIN) .GT. ABS(CGX+CXMAX) ) THEN
CGX = CGX + CXMIN
ELSE
CGX = CGX + CXMAX
END IF
IF ( ABS(CGY+CYMIN) .GT. ABS(CGY+CYMAX) ) THEN
CGY = CGY + CYMIN
ELSE
CGY = CGY + CYMAX
END IF
CXC = MAX ( ABS(CXMIN) , ABS(CXMAX) )
CYC = MAX ( ABS(CYMIN) , ABS(CYMAX) )
#ifdef W3_MGP
CXC = MAX ( ABS(CXMIN-VGX) , ABS(CXMAX-VGX) )
CYC = MAX ( ABS(CYMIN-VGY) , ABS(CYMAX-VGY) )
#endif
ELSE
CXC = 0.
CYC = 0.
END IF
!
CGN = MAX ( ABS(CGX) , ABS(CGY) , CXC, CYC, 0.001*CG0 )
!
NTLOC = 1 + INT(DTG/(DTCFL*CG0/CGN))
DTLOC = DTG / REAL(NTLOC)
DTRAD = DTLOC
IF ( FLAGLL ) DTRAD=DTRAD/(DERA*RADIUS)
!
TTEST(1) = TIME(1)
TTEST(2) = 0
DTTST = DSEC21(TTEST,TIME)
YFIRST = MOD(NINT(DTTST/DTG),2) .EQ. 0
!
#ifdef W3_T
WRITE (NDST,9000) YFIRST
WRITE (NDST,9001) ISP, ITH, IK, ECOS(ITH), ESIN(ITH)
#endif
!
! 1.b Initialize arrays
!
#ifdef W3_T
WRITE (NDST,9010)
#endif
!
VLCFLX = 0.
VLCFLY = 0.
CXTOT = 0.
CYTOT = 0.
!
MAPSTX(1:NX*NY) = MAPSTA(1:NX*NY)
!
IF ( GLOBAL ) THEN
MAPSTX(1-2*NY:0) = MAPSTA((NX-2)*NY+1:NX*NY)
MAPSTX(NX*NY+1:NX*NY+2*NY) = MAPSTA(1:2*NY)
ELSE
MAPSTX(1-2*NY:0) = 0
MAPSTX(NX*NY+1:NX*NY+2*NY) = 0
END IF
!
! 1.c Pre-calculate interpolation info
!
FTH = FGSE * WDTH * DTH * DTLOC
FCG = FGSE * WDCG * 0.5 * (XFR-1./XFR) * DTLOC
IF ( FLAGLL ) THEN
FTH = FTH / RADIUS / DERA
FCG = FCG / RADIUS / DERA
END IF
FCG = FCG / REAL(NTLOC) !TJC: required to match original (is this correct?)
FTHX = - FTH * ESIN(ITH)
FTHY = FTH * ECOS(ITH)
FCGX = FCG * ECOS(ITH)
FCGY = FCG * ESIN(ITH)
!
#ifdef W3_T0
WRITE (NDST,9011)
#endif
!
#ifdef W3_OMPH
!$OMP PARALLEL DO PRIVATE (ISEA, IX, IY, TMPX, TMPY, &
!$OMP& DXCGN, DYCGN, DXCGS, DYCGS, DXCGC, DYCGC, &
!$OMP& IXC, IYC)
#endif
!
DO ISEA=1, NSEA
!
IX = MAPSF(ISEA,1)
IY = MAPSF(ISEA,2)
!
! 1.c.1 Normal and parallel width ...
!
TMPX = FTHX / CLATS(ISEA)
TMPY = FTHY
DXCGN = DPDX(IY,IX)*TMPX + DPDY(IY,IX)*TMPY
DYCGN = DQDX(IY,IX)*TMPX + DQDY(IY,IX)*TMPY
TMPX = FCGX / CLATS(ISEA)
TMPY = FCGY
DXCGS = DPDX(IY,IX)*TMPX + DPDY(IY,IX)*TMPY
DYCGS = DQDX(IY,IX)*TMPX + DQDY(IY,IX)*TMPY
!
! 1.c.2 "Sum" corner (and mirror image) ...
!
DXCGC = DXCGN + DXCGS
DYCGC = DYCGN + DYCGS
!
IXC = NY
IF ( DXCGC .LT. 0. ) IXC = - IXC
IYC = 1
IF ( DYCGC .LT. 0. ) IYC = - IYC
!
IIXY1(ISEA) = IXC + IYC
IF ( ABS(DXCGC) .GT. ABS(DYCGC) ) THEN
IIXY2(ISEA) = IXC
RDI1 (ISEA) = ABS(DYCGC/DXCGC)
RDI2 (ISEA) = ABS(DXCGC)
ELSE
IIXY2(ISEA) = IYC
IF ( ABS(DYCGC) .GT. 1.E-5 ) THEN
RDI1(ISEA) = ABS(DXCGC/DYCGC)
ELSE
RDI1(ISEA) = 1.
END IF
RDI2(ISEA) = ABS(DYCGC)
END IF
!
#ifdef W3_T0
WRITE (NDST,9012) ISEA, ITH, IIXY1(ISEA), IIXY2(ISEA), &
RDI1(ISEA), RDI2(ISEA)*CG(IK,1)
#endif
!
! 1.c.2 "Difference" corner (and mirror image) ...
!
DXCGC = DXCGN - DXCGS
DYCGC = DYCGN - DYCGS
!
IXC = NY
IF ( DXCGC .LT. 0. ) IXC = - IXC
IYC = 1
IF ( DYCGC .LT. 0. ) IYC = - IYC
!
IIXY3(ISEA) = IXC + IYC
IF ( ABS(DXCGC) .GT. ABS(DYCGC) ) THEN
IIXY4(ISEA) = IXC
RDI3 (ISEA) = ABS(DYCGC/DXCGC)
RDI4 (ISEA) = ABS(DXCGC)
ELSE
IIXY4(ISEA) = IYC
IF ( ABS(DYCGC) .GT. 1.E-5 ) THEN
RDI3(ISEA) = ABS(DXCGC/DYCGC)
ELSE
RDI3(ISEA) = 1.
END IF
RDI4(ISEA) = ABS(DYCGC)
END IF
!
#ifdef W3_T0
WRITE (NDST,9013) IIXY3(ISEA), IIXY4(ISEA), RDI3(ISEA), &
RDI4(ISEA)*CG(IK,1)
#endif
!
END DO
!
#ifdef W3_OMPH
!$OMP END PARALLEL DO
#endif
!
! 2. Calculate velocities and diffusion coefficients ---------------- *
! 2.a Velocities
!
! Q = ( A / CG * CLATS )
! LCFLX = ( COS*CG / CLATS ) * DT / DX
! LCFLY = ( SIN*CG ) * DT / DY
!
#ifdef W3_T
WRITE (NDST,9020) NSEA
#endif
!
#ifdef W3_OMPH
!$OMP PARALLEL DO PRIVATE (ISEA, IXY)
#endif
!
DO ISEA=1, NSEA
IXY = MAPSF(ISEA,3)
VQ (IXY) = VQ(IXY) / CG(IK,ISEA) * CLATS(ISEA)
CXTOT(IXY) = ECOS(ITH) * CG(IK,ISEA) / CLATS(ISEA)
CYTOT(IXY) = ESIN(ITH) * CG(IK,ISEA)
#ifdef W3_MGP
CXTOT(IXY) = CXTOT(IXY) - VGX/CLATS(ISEA)
CYTOT(IXY) = CYTOT(IXY) - VGY
#endif
#ifdef W3_T1
IF ( .NOT. FLCUR ) &
WRITE (NDST,9021) ISEA, MAPSF(ISEA,1), MAPSF(ISEA,2), &
VQ(IXY), CXTOT(IXY), CYTOT(IXY)
#endif
END DO
!
#ifdef W3_OMPH
!$OMP END PARALLEL DO
#endif
!
IF ( FLCUR ) THEN
#ifdef W3_T
WRITE (NDST,9022)
#endif
!
#ifdef W3_OMPH
!$OMP PARALLEL DO PRIVATE (ISEA,IXY)
#endif
!
DO ISEA=1, NSEA
IXY = MAPSF(ISEA,3)
CXTOT(IXY) = CXTOT(IXY) + CX(ISEA)/CLATS(ISEA)
CYTOT(IXY) = CYTOT(IXY) + CY(ISEA)
#ifdef W3_T1
WRITE (NDST,9021) ISEA, MAPSF(ISEA,1), MAPSF(ISEA,2), &
VQ(IXY), CXTOT(IXY), CYTOT(IXY)
#endif
END DO
!
#ifdef W3_OMPH
!$OMP END PARALLEL DO
#endif
!
END IF
!
#ifdef W3_OMPH
!$OMP PARALLEL DO PRIVATE (ISEA,IX, IY, IXY, CP, CQ)
#endif
!
DO ISEA=1, NSEA
IX = MAPSF(ISEA,1)
IY = MAPSF(ISEA,2)
IXY = MAPSF(ISEA,3)
CP = CXTOT(IXY)*DPDX(IY,IX) + CYTOT(IXY)*DPDY(IY,IX)
CQ = CXTOT(IXY)*DQDX(IY,IX) + CYTOT(IXY)*DQDY(IY,IX)
VLCFLX(IXY) = CP*DTRAD
VLCFLY(IXY) = CQ*DTRAD
END DO
!
#ifdef W3_OMPH
!$OMP END PARALLEL DO
#endif
!
! 3. Loop over sub-steps -------------------------------------------- *
!
DO ITLOC=1, NTLOC
!
! 3.a Average
!
AQ = VQ
VQ = 0.
!
! 3.a.1 Central points
!
DO IP=1, NCENT
ISEA = MAPCXY(IP)
IXY = MAPSF(ISEA,3)
IF ( MAPTRN(IXY) ) THEN
VQ(IXY) = AQ(IXY)
ELSE
RD1 = RDI1(ISEA)
RD2 = MIN ( 1. , RDI2(ISEA) * CG(IK,ISEA) )
RD3 = RDI3(ISEA)
RD4 = MIN ( 1. , RDI4(ISEA) * CG(IK,ISEA) )
VQ(IXY ) = VQ(IXY ) &
+ AQ(IXY) * (3.-RD2-RD4)/3.
VQ(IXY+IIXY1(ISEA)) = VQ(IXY+IIXY1(ISEA)) &
+ AQ(IXY) * RD2*RD1/6.
VQ(IXY+IIXY2(ISEA)) = VQ(IXY+IIXY2(ISEA)) &
+ AQ(IXY) * (1.-RD1)*RD2/6.
VQ(IXY+IIXY3(ISEA)) = VQ(IXY+IIXY3(ISEA)) &
+ AQ(IXY) * RD4*RD3/6.
VQ(IXY+IIXY4(ISEA)) = VQ(IXY+IIXY4(ISEA)) &
+ AQ(IXY) * (1.-RD3)*RD4/6.
VQ(IXY-IIXY1(ISEA)) = VQ(IXY-IIXY1(ISEA)) &
+ AQ(IXY) * RD2*RD1/6.
VQ(IXY-IIXY2(ISEA)) = VQ(IXY-IIXY2(ISEA)) &
+ AQ(IXY) * (1.-RD1)*RD2/6.
VQ(IXY-IIXY3(ISEA)) = VQ(IXY-IIXY3(ISEA)) &
+ AQ(IXY) * RD4*RD3/6.
VQ(IXY-IIXY4(ISEA)) = VQ(IXY-IIXY4(ISEA)) &
+ AQ(IXY) * (1.-RD3)*RD4/6.
END IF
END DO
!
! 3.a.2 Near-coast points
!
DO IP=NCENT+1, NSEA
ISEA = MAPCXY(IP)
IX = MAPSF(ISEA,1)
IXY = MAPSF(ISEA,3)
IF ( MAPSTA(IXY) .LE. 0 ) CYCLE
IF ( MAPTRN(IXY) ) THEN
VQ(IXY) = AQ(IXY)
ELSE
RD1 = RDI1(ISEA)
RD3 = RDI3(ISEA)
RD2 = MIN ( 1. , RDI2(ISEA) * CG(IK,ISEA) )
RD4 = MIN ( 1. , RDI4(ISEA) * CG(IK,ISEA) )
VQ(IXY ) = VQ(IXY ) &
+ AQ(IXY) * (3.-RD2-RD4)/3.
!
IXC = SIGN(NY,IIXY1(ISEA))
IYC = IIXY1(ISEA) - IXC
IF ( MAPSTX(IXY+IIXY1(ISEA)) .GE. 1 .AND. &
.NOT. ( MAPSTX(IXY+IXC).LE.0 .AND. &
MAPSTX(IXY+IYC).LE.0 ) ) THEN
VQ(IXY+IIXY1(ISEA)) = VQ(IXY+IIXY1(ISEA)) &
+ AQ(IXY) * RD2*RD1/6.
ELSE
VQ(IXY ) = VQ(IXY ) &
+ AQ(IXY) * RD2*RD1/6.
END IF
IF ( MAPSTX(IXY-IIXY1(ISEA)) .GE. 1 .AND. &
.NOT. ( MAPSTX(IXY-IXC).LE.0 .AND. &
MAPSTX(IXY-IYC).LE.0 ) ) THEN
VQ(IXY-IIXY1(ISEA)) = VQ(IXY-IIXY1(ISEA)) &
+ AQ(IXY) * RD2*RD1/6.
ELSE
VQ(IXY ) = VQ(IXY ) &
+ AQ(IXY) * RD2*RD1/6.
END IF
IF ( MAPSTX(IXY+IIXY2(ISEA)) .GE. 1 ) THEN
VQ(IXY+IIXY2(ISEA)) = VQ(IXY+IIXY2(ISEA)) &
+ AQ(IXY) * (1.-RD1)*RD2/6.
ELSE
VQ(IXY ) = VQ(IXY ) &
+ AQ(IXY) * (1.-RD1)*RD2/6.
END IF
IF ( MAPSTX(IXY-IIXY2(ISEA)) .GE. 1 ) THEN
VQ(IXY-IIXY2(ISEA)) = VQ(IXY-IIXY2(ISEA)) &
+ AQ(IXY) * (1.-RD1)*RD2/6.
ELSE
VQ(IXY ) = VQ(IXY ) &
+ AQ(IXY) * (1.-RD1)*RD2/6.
END IF
!
IXC = SIGN(NY,IIXY3(ISEA))
IYC = IIXY3(ISEA) - IXC
IF ( MAPSTX(IXY+IIXY3(ISEA)) .GE. 1 .AND. &
.NOT. ( MAPSTX(IXY+IXC).LE.0 .AND. &
MAPSTX(IXY+IYC).LE.0 ) ) THEN
VQ(IXY+IIXY3(ISEA)) = VQ(IXY+IIXY3(ISEA)) &
+ AQ(IXY) * RD4*RD3/6.
ELSE
VQ(IXY ) = VQ(IXY ) &
+ AQ(IXY) * RD4*RD3/6.
END IF
IF ( MAPSTX(IXY-IIXY3(ISEA)) .GE. 1 .AND. &
.NOT. ( MAPSTX(IXY-IXC).LE.0 .AND. &
MAPSTX(IXY-IYC).LE.0 ) ) THEN
VQ(IXY-IIXY3(ISEA)) = VQ(IXY-IIXY3(ISEA)) &
+ AQ(IXY) * RD4*RD3/6.
ELSE
VQ(IXY ) = VQ(IXY ) &
+ AQ(IXY) * RD4*RD3/6.
END IF
!
IF ( MAPSTX(IXY+IIXY4(ISEA)) .GE. 1 ) THEN
VQ(IXY+IIXY4(ISEA)) = VQ(IXY+IIXY4(ISEA)) &
+ AQ(IXY) * (1.-RD3)*RD4/6.
ELSE
VQ(IXY ) = VQ(IXY ) &
+ AQ(IXY) * (1.-RD3)*RD4/6.
END IF
IF ( MAPSTX(IXY-IIXY4(ISEA)) .GE. 1 ) THEN
VQ(IXY-IIXY4(ISEA)) = VQ(IXY-IIXY4(ISEA)) &
+ AQ(IXY) * (1.-RD3)*RD4/6.
ELSE
VQ(IXY ) = VQ(IXY ) &
+ AQ(IXY) * (1.-RD3)*RD4/6.
END IF
!
END IF
!
END DO
!
! 3.a.3 Restore boundary data
!
DO IXY=1, NX*NY
IF ( MAPSTA(IXY).EQ.2 ) VQ(IXY) = AQ(IXY)
END DO
!
! 3.a.4 Global closure (averaging only, propagation is closed in W3QCK3).
!
IF ( GLOBAL ) THEN
DO IY=1, NY
VQ(IY ) = VQ(IY ) + VQ(NX*NY+IY)
VQ((NX-1)*NY+IY) = VQ((NX-1)*NY+IY) + VQ(IY-NY)
END DO
END IF
!
! 3.b Propagate fields
!
! Transform VQ to straightened space
!
#ifdef W3_OMPH
!$OMP PARALLEL DO PRIVATE (ISEA, IX, IY, IXY)
#endif
!
DO ISEA=1, NSEA
IX = MAPSF(ISEA,1)
IY = MAPSF(ISEA,2)
IXY = MAPSF(ISEA,3)
VQ(IXY)= VQ(IXY) * GSQRT(IY,IX)
END DO
!
#ifdef W3_OMPH
!$OMP END PARALLEL DO
#endif
!
IF ( YFIRST ) THEN
!
#ifdef W3_UQ
IF ( FLCY ) CALL W3QCK3 &
(NX, NY, NX, NY, VLCFLY, ATRNY, VQ, &
.FALSE., 1, MAPAXY, NACT, MAPY2, NMY0, &
NMY1, NMY2, NDSE, NDST )
IF ( FLCX ) CALL W3QCK3 &
(NX, NY, NX, NY, VLCFLX, ATRNX, VQ, &
GLOBAL, NY, MAPAXY, NACT, MAPX2, NMX0, &
NMX1, NMX2, NDSE, NDST )
#endif
!
#ifdef W3_UNO
IF ( FLCY ) CALL W3UNO2s &
(NX, NY, NX, NY, VLCFLY, ATRNY, VQ, &
.FALSE., 1, MAPAXY, NACT, MAPY2, NMY0, &
NMY1, NMY2, NDSE, NDST )
IF ( FLCX ) CALL W3UNO2s &
(NX, NY, NX, NY, VLCFLX, ATRNX, VQ, &
GLOBAL, NY, MAPAXY, NACT, MAPX2, NMX0, &
NMX1, NMX2, NDSE, NDST )
#endif
!
ELSE
!
#ifdef W3_UQ
IF ( FLCX ) CALL W3QCK3 &
(NX, NY, NX, NY, VLCFLX, ATRNX, VQ, &
GLOBAL, NY, MAPAXY, NACT, MAPX2, NMX0, &
NMX1, NMX2, NDSE, NDST )
IF ( FLCY ) CALL W3QCK3 &
(NX, NY, NX, NY, VLCFLY, ATRNY, VQ, &
.FALSE., 1, MAPAXY, NACT, MAPY2, NMY0, &
NMY1, NMY2, NDSE, NDST )
#endif
!
#ifdef W3_UNO
IF ( FLCX ) CALL W3UNO2s &
(NX, NY, NX, NY, VLCFLX, ATRNX, VQ, &
GLOBAL, NY, MAPAXY, NACT, MAPX2, NMX0, &
NMX1, NMX2, NDSE, NDST )
IF ( FLCY ) CALL W3UNO2s &
(NX, NY, NX, NY, VLCFLY, ATRNY, VQ, &
.FALSE., 1, MAPAXY, NACT, MAPY2, NMY0, &
NMY1, NMY2, NDSE, NDST )
#endif
!
END IF
!
! Transform VQ back to normal space
!
#ifdef W3_OMPH
!$OMP PARALLEL DO PRIVATE (ISEA, IX, IY, IXY)
#endif
!
DO ISEA=1, NSEA
IX = MAPSF(ISEA,1)
IY = MAPSF(ISEA,2)
IXY = MAPSF(ISEA,3)
VQ(IXY)= VQ(IXY) / GSQRT(IY,IX)
END DO
!
#ifdef W3_OMPH
!$OMP END PARALLEL DO
#endif
!
! 3.c Update boundaries
!
#ifdef W3_T
WRITE (NDST,9030) NSEA
#endif
!
#ifdef W3_T2
DO ISEA=1, NSEA
IXY = MAPSF(ISEA,3)
IF ( MAPSTA(IXY) .GT. 0 ) THEN
WRITE(NDST,9031)ISEA,MAPSF(ISEA,1),MAPSF(ISEA,2),VQ(IXY)
VQ(IXY) = MAX ( 0. , CG(IK,ISEA)/CLATS(ISEA)*VQ(IXY) )
END IF
END DO
#endif
!
IF ( FLBPI ) THEN
RD1 = DSEC21 ( TBPI0, TIME ) - DTG * &
REAL(NTLOC-ITLOC)/REAL(NTLOC)
RD2 = DSEC21 ( TBPI0, TBPIN )
IF ( RD2 .GT. 0.001 ) THEN
RD2 = MIN(1.,MAX(0.,RD1/RD2))
RD1 = 1. - RD2
ELSE
RD1 = 0.
RD2 = 1.
END IF
DO IBI=1, NBI
IXY = MAPSF(ISBPI(IBI),3)
VQ(IXY) = ( RD1*BBPI0(ISP,IBI) + RD2*BBPIN(ISP,IBI) ) &
/ CG(IK,ISBPI(IBI)) * CLATS(ISBPI(IBI))
END DO
END IF
!
YFIRST = .NOT. YFIRST
END DO
!
! 4. Store results in VQ in proper format --------------------------- *
!
#ifdef W3_T
WRITE (NDST,9040) NSEA
#endif
!
#ifdef W3_OMPH
!$OMP PARALLEL DO PRIVATE (ISEA, IXY)
#endif
!
DO ISEA=1, NSEA
IXY = MAPSF(ISEA,3)
IF ( MAPSTA(IXY) .GT. 0 ) THEN
#ifdef W3_T3
WRITE (NDST,9041) ISEA, MAPSF(ISEA,1), MAPSF(ISEA,2), VQ(IXY)
#endif
VQ(IXY) = MAX ( 0. , CG(IK,ISEA)/CLATS(ISEA)*VQ(IXY) )
END IF
END DO
!
#ifdef W3_OMPH
!$OMP END PARALLEL DO
#endif
!
RETURN
!
! Formats
!
#ifdef W3_T
9000 FORMAT (' TEST W3XYP3 : YFIRST :',L2)
9001 FORMAT (' TEST W3XYP3 : ISP, ITH, IK, COS-SIN :',I8,2I4,2F7.3)
9010 FORMAT (' TEST W3XYP3 : INITIALIZE ARRAYS')
9020 FORMAT (' TEST W3XYP3 : CALCULATING VCFL0X/Y (NSEA=',I6,')')
9022 FORMAT (' TEST W3XYP3 : CALCULATING VCFLUX/Y')
9030 FORMAT (' TEST W3XYP3 : FIELD BEFORE BPI. (NSEA=',I6,')')
9040 FORMAT (' TEST W3XYP3 : FIELD AFTER PROP. (NSEA=',I6,')')
#endif
#ifdef W3_T0
9011 FORMAT (' TEST W3XYP3 : PREPARE AVERAGING')
9012 FORMAT (' ',4I4,2F7.3)
9013 FORMAT (' ',8X,2I4,2F7.3)
#endif
!
#ifdef W3_T1
9021 FORMAT (1X,I6,2I5,E12.4,2f7.3)
#endif
!
#ifdef W3_T2
9031 FORMAT (1X,I6,2I5,E12.4)
#endif
!
#ifdef W3_T3
9041 FORMAT (1X,I6,2I5,E12.4)
#endif
!/
!/ End of W3XYP3 ----------------------------------------------------- /
!/
END SUBROUTINE W3XYP3
!/ ------------------------------------------------------------------- /
SUBROUTINE W3KTP3 ( ISEA, FACTH, FACK, CTHG0, CG, WN, DW, &
DDDX, DDDY, CX, CY, DCXDX, DCXDY, &
DCYDX, DCYDY, DCDX, DCDY, VA, CFLTHMAX, CFLKMAX )
!/
!/ *** THIS ROUTINE SHOULD BE IDENTICAL TO W3KTP2 ***
!/
!/ +-----------------------------------+
!/ | WAVEWATCH III NOAA/NCEP |
!/ | H. L. Tolman |
!/ | FORTRAN 90 |
!/ | Last update : 01-Jul-2013 |
!/ +-----------------------------------+
!/
!/ 14-Feb-2000 : Origination. ( version 2.08 )
!/ 17-Dec-2004 : Multiple grid version. ( version 3.06 )
!/ 06-Mar-2011 : Output of maximum CFL (F. Ardhuin) ( version 3.14 )
!/ 24-Aug-2011 : Limiter on k advection (F. Ardhuin) ( version 4.04 )
!/ 25-Aug-2011 : DEPTH = MAX ( DMIN, DW(ISEA) ) ( version 4.04 )
!/ 26-Dec-2012 : More initializations. ( version 4.11 )
!/ 01-Jul-2013 : Adding UQ and UNO switches to chose between third
!/ and second order schemes. ( version 4.12 )
!/
! 1. Purpose :
!
! Propagation in spectral space.
!
! 2. Method :
!
! Third order QUICKEST scheme with ULTIMATE limiter.
!
! As with the spatial propagation, the two spaces are considered
! independently, but the propagation is performed in a 2-D space.
! Compared to the propagation in physical space, the directions
! rerpesent a closed space and are therefore comparable to the
! longitudinal or 'X' propagation. The wavenumber space has to be
! extended to allow for boundary treatment. Using a simple first
! order boundary treatment at both sided, two points need to
! be added. This implies that the spectrum needs to be extended,
! shifted and rotated, as is performed using MAPTH2 as set
! in W3MAP3.
!
! 3. Parameters :
!
! Parameter list
! ----------------------------------------------------------------
! ISEA Int. I Number of sea point.
! FACTH/K Real I Factor in propagation velocity.
! CTHG0 Real I Factor in great circle refracftion term.
! MAPxx2 I.A. I Propagation and storage maps.
! CG R.A. I Local group velocities.
! WN R.A. I Local wavenumbers.
! DW R.A. I Depth.
! DDDx Real I Depth gradients.
! CX/Y Real I Current components.
! DCxDx Real I Current gradients.
! DCDX-Y Real I Phase speed gradients.
! VA R.A. I/O Spectrum.
! ----------------------------------------------------------------
!
! Local variables.
! ----------------------------------------------------------------
! DSDD R.A. Partial derivative of sigma for depth.
! FDD, FDU, FDG, FCD, FCU
! R.A. Directionally varying part of depth, current and
! great-circle refraction terms and of consit.
! of Ck term.
! CFLT-K R.A. Propagation velocities of local fluxes.
! DB R.A. Wavenumber band widths at cell centers.
! DM R.A. Wavenumber band widths between cell centers and
! next cell center.
! Q R.A. Extracted spectrum
! ----------------------------------------------------------------
!
! 4. Subroutines used :
!
! W3QCK1 Actual propagation routine.
! W3QCK2 Actual propagation routine.
! STRACE Service routine.
!
! 5. Called by :
!
! W3WAVE Wave model routine.
!
! 6. Error messages :
!
! None.
!
! 8. Structure :
!
! -----------------------------------------------------------------
! 1. Preparations
! a Initialize arrays
! b Set constants and counters
! 2. Point preparations
! a Calculate DSDD
! b Extract spectrum
! 3. Refraction velocities
! a Filter level depth reffraction.
! b Depth refratcion velocity.
! c Current refraction velocity.
! 4. Wavenumber shift velocities
! a Prepare directional arrays
! b Calcuate velocity.
! 5. Propagate.
! 6. Store results.
! -----------------------------------------------------------------
!
! 9. Switches :
!
! !/S Enable subroutine tracing.
! !/T Enable general test output.
!
! 10. Source code :
!
!/ ------------------------------------------------------------------- /
USE CONSTANTS
USE W3GDATMD, ONLY: NK, NK2, NTH, NSPEC, SIG, DSIP, ECOS, ESIN, &
EC2, ESC, ES2, FACHFA, MAPWN, FLCTH, FLCK, &
CTMAX, DMIN
USE W3ADATMD, ONLY: MAPTH2, MAPWN2, ITIME
USE W3IDATMD, ONLY: FLCUR
USE W3ODATMD, ONLY: NDSE, NDST
#ifdef W3_S
USE W3SERVMD, ONLY: STRACE
#endif
#ifdef W3_UQ
USE W3UQCKMD
#endif
#ifdef W3_UNO
USE W3UNO2MD
#endif
!/
IMPLICIT NONE
!/
!/ ------------------------------------------------------------------- /
!/ Parameter list
!/
INTEGER, INTENT(IN) :: ISEA
REAL, INTENT(IN) :: FACTH, FACK, CTHG0, CG(0:NK+1), &
WN(0:NK+1), DW, DDDX, DDDY, &
CX, CY, DCXDX, DCXDY, DCYDX, DCYDY
REAL, INTENT(IN) :: DCDX(0:NK+1), DCDY(0:NK+1)
REAL, INTENT(INOUT) :: VA(NSPEC)
REAL, INTENT(OUT) :: CFLTHMAX, CFLKMAX
!/
!/ ------------------------------------------------------------------- /
!/ Local parameters
!/
INTEGER :: ITH, IK, ISP
#ifdef W3_S
INTEGER, SAVE :: IENT = 0
#endif
REAL :: FDDMAX, FDG, FKD, FKD0, DCYX, &
DCXXYY, DCXY, DCXX, DCXYYX, DCYY, &
VELNOFILT, VELFAC, DEPTH
REAL :: DSDD(0:NK+1), FRK(NK), FRG(NK), &
FKC(NTH), VQ(-NK-1:NK2*(NTH+2)), &
DB(NK2,NTH+1), DM(NK2,0:NTH+1), &
VCFLT(NK2*(NTH+1)), CFLK(NK2,NTH)
!/
!/ ------------------------------------------------------------------- /
!/
#ifdef W3_S
CALL STRACE (IENT, 'W3KTP3')
#endif
!
! 1. Preparations --------------------------------------------------- *
! 1.a Initialize arrays
!
DEPTH = MAX ( DMIN, DW )
VQ = 0.
IF ( FLCTH ) VCFLT = 0.
IF ( FLCK ) CFLK = 0.
CFLTHMAX = 0.
CFLKMAX = 0.
!
#ifdef W3_T
WRITE (NDST,9000) FLCTH, FLCK, FACTH, FACK, CTMAX
WRITE (NDST,9010) ISEA, DEPTH, CX, CY, DDDX, DDDY, &
DCXDX, DCXDY, DCYDX, DCYDY
#endif
!
! 2. Preparation for point ------------------------------------------ *
! 2.a Array with partial derivative of sigma versus depth
!
DO IK=0, NK+1
IF ( DEPTH*WN(IK) .LT. 5. ) THEN
DSDD(IK) = MAX ( 0. , &
CG(IK)*WN(IK)-0.5*SIG(IK) ) / DEPTH
ELSE
DSDD(IK) = 0.
END IF
END DO
!
#ifdef W3_T
WRITE (NDST,9020)
DO IK=1, NK+1
WRITE (NDST,9021) IK, TPI/SIG(IK), TPI/WN(IK), &
CG(IK), DSDD(IK)
END DO
#endif
!
! 2.b Extract spectrum
!
DO ISP=1, NSPEC
VQ(MAPTH2(ISP)) = VA(ISP)
END DO
!
! 3. Refraction velocities ------------------------------------------ *
!
IF ( FLCTH ) THEN
!
! 3.a Set slope filter for depth refraction
!
! N.B.: FACTH = DTG / DTH / REAL(NTLOC) (value set in w3wavemd)
! namely, FACTH*VC=1 corresponds to CFL=1
!
FDDMAX = 0.
FDG = FACTH * CTHG0
!
DO ITH=1, NTH/2
FDDMAX = MAX(FDDMAX,ABS(ESIN(ITH)*DDDX-ECOS(ITH)*DDDY))
END DO
!
DO IK=1, NK
FRK(IK) = FACTH * DSDD(IK) / WN(IK)
!
! Removes the filtering that was done at that stage (F. Ardhuin 2011/03/06)
!
! FRK(IK) = FRK(IK) / MAX ( 1. , FRK(IK)*FDDMAX/CTMAX )
FRG(IK) = FDG * CG(IK)
END DO
!
! 3.b Current refraction
!
IF ( FLCUR ) THEN
!
DCYX = FACTH * DCYDX
DCXXYY = FACTH * ( DCXDX - DCYDY )
DCXY = FACTH * DCXDY
!
DO ISP=1, NSPEC
VCFLT(MAPTH2(ISP)) = ES2(ISP)*DCYX + &
ESC(ISP)*DCXXYY - EC2(ISP)*DCXY
END DO
!
ELSE
VCFLT(:)=0.
END IF
!
! 3.c Depth refraction and great-circle propagation
!
#ifdef W3_REFRX
! 3.d @C/@x refraction and great-circle propagation
FRK = 0.
FDDMAX = 0.
!
DO ISP=1, NSPEC
FDDMAX = MAX ( FDDMAX , ABS ( &
ESIN(ISP)*DCDX(MAPWN(ISP)) - ECOS(ISP)*DCDY(MAPWN(ISP)) ) )
END DO
DO IK=1, NK
FRK(IK) = FACTH * CG(IK) * WN(IK) / SIG(IK)
END DO
#endif
!
DO ISP=1, NSPEC
VELNOFILT = VCFLT(MAPTH2(ISP)) &
+ FRG(MAPWN(ISP)) * ECOS(ISP) &
+ FRK(MAPWN(ISP)) * ( ESIN(ISP)*DDDX - ECOS(ISP)*DDDY )
!
#ifdef W3_REFRX
! 3.d @C/@x refraction and great-circle propagation
VELNOFILT = VCFLT(MAPTH2(ISP)) &
+ FRG(MAPWN(ISP)) * ECOS(ISP) &
+ FRK(MAPWN(ISP)) * ( ESIN(ISP)*DCDX(MAPWN(ISP)) &
- ECOS(ISP)*DCDY(MAPWN(ISP)) )
#endif
CFLTHMAX = MAX(CFLTHMAX, ABS(VELNOFILT))
!
! Puts filtering on total velocity (including currents and great circle effects)
! the filtering limits VCFLT to be less than CTMAX
! this modification was proposed by F. Ardhuin 2011/03/06
!
VCFLT(MAPTH2(ISP))=SIGN(MIN(ABS(VELNOFILT),CTMAX),VELNOFILT)
END DO
END IF
!
! 4. Wavenumber shift velocities ------------------------------------ *
! N.B.: FACK = DTG / REAL(NTLOC) (value set in w3wavemd)
! namely, FACK*VC/DK=1 corresponds to CFL=1
!
IF ( FLCK ) THEN
!
! 4.a Directionally dependent part
!
DCXX = - DCXDX
DCXYYX = - ( DCXDY + DCYDX )
DCYY = - DCYDY
FKD = ( CX*DDDX + CY*DDDY )
!
DO ITH=1, NTH
FKC(ITH) = EC2(ITH)*DCXX + &
ESC(ITH)*DCXYYX + ES2(ITH)*DCYY
END DO
!
! 4.b Band widths
!
DO IK=0, NK
DB(IK+1,1) = DSIP(IK) / CG(IK)
DM(IK+1,1) = WN(IK+1) - WN(IK)
END DO
DB(NK+2,1) = DSIP(NK+1) / CG(NK+1)
DM(NK+2,1) = 0.
!
DO ITH=2, NTH
DO IK=1, NK+2
DB(IK,ITH) = DB(IK,1)
DM(IK,ITH) = DM(IK,1)
END DO
END DO
!
! 4.c Velocities
!
DO IK=0, NK+1
FKD0 = FKD / CG(IK) * DSDD(IK)
VELFAC = FACK/DB(IK+1,1)
DO ITH=1, NTH
!
! Puts filtering on velocity (needs the band widths)
!
VELNOFILT = ( FKD0 + WN(IK)*FKC(ITH) ) * VELFAC ! this is velocity * DT / DK
CFLKMAX = MAX(CFLKMAX, ABS(VELNOFILT))
CFLK(IK+1,ITH) = SIGN(MIN(ABS(VELNOFILT),CTMAX),VELNOFILT)/VELFAC
!CFLK(IK+1,ITH) = FKD0 + WN(IK)*FKC(ITH) ! this was without the limiter ...
END DO
END DO
!
END IF
!
! 5. Propagate ------------------------------------------------------ *
!
IF ( MOD(ITIME,2) .EQ. 0 ) THEN
IF ( FLCK ) THEN
DO ITH=1, NTH
VQ(NK+2+(ITH-1)*NK2) = FACHFA * VQ(NK+1+(ITH-1)*NK2)
END DO
!
#ifdef W3_UQ
CALL W3QCK2 ( NTH, NK2, NTH, NK2, CFLK, FACK, DB, DM, &
VQ, .FALSE., 1, MAPTH2, NSPEC, &
MAPWN2, NSPEC-NTH, NSPEC, NSPEC+NTH, &
NDSE, NDST )
#endif
!
#ifdef W3_UNO
CALL W3UNO2 ( NTH, NK2, NTH, NK2, CFLK, FACK, DB, DM, &
VQ, .FALSE., 1, MAPTH2, NSPEC, &
MAPWN2, NSPEC-NTH, NSPEC, NSPEC+NTH, &
NDSE, NDST )
#endif
!
END IF
IF ( FLCTH ) THEN
!
#ifdef W3_UQ
CALL W3QCK1 ( NTH, NK2, NTH, NK2, VCFLT, VQ, .TRUE., &
NK2, MAPTH2, NSPEC, MAPTH2, NSPEC, NSPEC, &
NSPEC, NDSE, NDST )
#endif
!
#ifdef W3_UNO
CALL W3UNO2r( NTH, NK2, NTH, NK2, VCFLT, VQ, .TRUE., &
NK2, MAPTH2, NSPEC, MAPTH2, NSPEC, NSPEC,&
NSPEC, NDSE, NDST )
#endif
!
END IF
ELSE
IF ( FLCTH ) THEN
!
#ifdef W3_UQ
CALL W3QCK1 ( NTH, NK2, NTH, NK2, VCFLT, VQ, .TRUE., &
NK2, MAPTH2, NSPEC, MAPTH2, NSPEC, NSPEC, &
NSPEC, NDSE, NDST )
#endif
!
#ifdef W3_UNO
CALL W3UNO2r( NTH, NK2, NTH, NK2, VCFLT, VQ, .TRUE., &
NK2, MAPTH2, NSPEC, MAPTH2, NSPEC, NSPEC,&
NSPEC, NDSE, NDST )
#endif
!
END IF
IF ( FLCK ) THEN
DO ITH=1, NTH
VQ(NK+2+(ITH-1)*NK2) = FACHFA * VQ(NK+1+(ITH-1)*NK2)
END DO
!
#ifdef W3_UQ
CALL W3QCK2 ( NTH, NK2, NTH, NK2, CFLK, FACK, DB, DM, &
VQ, .FALSE., 1, MAPTH2, NSPEC, &
MAPWN2, NSPEC-NTH, NSPEC, NSPEC+NTH, &
NDSE, NDST )
#endif
!
#ifdef W3_UNO
CALL W3UNO2 ( NTH, NK2, NTH, NK2, CFLK, FACK, DB, DM, &
VQ, .FALSE., 1, MAPTH2, NSPEC, &
MAPWN2, NSPEC-NTH, NSPEC, NSPEC+NTH, &
NDSE, NDST )
#endif
!
END IF
END IF
!
! 6. Store reults --------------------------------------------------- *
!
DO ISP=1, NSPEC
VA(ISP) = VQ(MAPTH2(ISP))
END DO
!
RETURN
!
! Formats
!
#ifdef W3_T
9000 FORMAT ( ' TEST W3KTP3 : FLCTH-K, FACTH-K, CTMAX :', &
2L2,2E10.3,F7.3)
9010 FORMAT ( ' TEST W3KTP3 : LOCAL DATA :',I7,F7.1,2F6.2,1X,6E10.2)
9020 FORMAT ( ' TEST W3KTP3 : IK, T, L, CG, DSDD : ')
9021 FORMAT ( ' ',I3,F7.2,F7.1,F7.2,E11.3)
#endif
!
#ifdef W3_T0
9040 FORMAT (/' TEST W3KTP3 : NORMALIZED ',A/)
9041 FORMAT (1X,60(1X,I2))
9042 FORMAT (1X,60I3)
#endif
!/
!/ End of W3KTP3 ----------------------------------------------------- /
!/
END SUBROUTINE W3KTP3
!/ ------------------------------------------------------------------- /
SUBROUTINE W3CFLXY ( ISEA, DTG, MAPSTA, MAPFS, CFLXYMAX, VGX, VGY )
!/
!/ +-----------------------------------+
!/ | WAVEWATCH III NOAA/NCEP |
!/ | F. Ardhuin |
!/ | FORTRAN 90 |
!/ | Last update : 31-Oct-2010 |
!/ +-----------------------------------+
!/
!/ 07-Mar-2011 : Origination. ( version 3.14 )
!/
! 1. Purpose :
!
! Computes the maximum CFL number for spatial advection. Used for diagnostic
! purposes. (Could be used to define a local time step ...)
!
! 2. Method :
!
!
! 3. Parameters :
!
! Parameter list
! ----------------------------------------------------------------
! ISEA Int. I Index of grid point.
! DTG Real I Total time step.
! MAPSTA I.A. I Grid point status map.
! MAPFS I.A. I Storage map.
! CFLXYMAX Real O Maximum CFL number for XY propagation.
! VGX/Y Real I Speed of grid.
! ----------------------------------------------------------------
!
! Local variables.
! ----------------------------------------------------------------
! NTLOC Int Number of local time steps.
! DTLOC Real Local propagation time step.
! VCFL0X R.A. Local courant numbers for absolute group vel.
! using local X-grid step.
! VCFL0Y R.A. Id. in Y.
! ----------------------------------------------------------------
!
! 4. Subroutines used :
!
! STRACE Service routine.
!
! 5. Called by :
!
! W3WAVE Wave model routine.
!
! 6. Error messages :
!
! None.
!
! 7. Remarks :
!
! - Curvilinear grid implementation. Variables FACX, FACY, CCOS, CSIN,
! CCURX, CCURY are not needed and have been removed. FACX is accounted
! for as approriate in this subroutine. FACX is also accounted for in
! the case of .NOT.FLCX. Since FACX is removed, there is now a check for
! .NOT.FLCX in this subroutine. In CFL calcs dx and dy are omitted,
! since dx=dy=1 in index space. Curvilinear grid derivatives
! (DPDY, DQDX, etc.) and metric (GSQRT) are brought in via W3GDATMD.
!
! 8. Structure :
!
! ---------------------------------------------
! ---------------------------------------------
!
! 9. Switches :
!
! !/S Enable subroutine tracing.
!
! !/MGP Moving grid corrections.
! !/MGG Moving grid corrections.
!
! !/T Enable general test output.
!
! 10. Source code :
!
!/ ------------------------------------------------------------------- /
USE CONSTANTS
!
USE W3TIMEMD, ONLY: DSEC21
!
USE W3GDATMD, ONLY: NX, NY, NSEA, MAPSF, DTCFL, CLATS, &
FLCX, FLCY, NK, NTH, DTH, XFR, &
ECOS, ESIN, SIG, WDCG, WDTH, PFMOVE, &
FLAGLL, DPDX, DPDY, DQDX, DQDY, GSQRT
USE W3WDATMD, ONLY: TIME
USE W3ADATMD, ONLY: NMX0, NMX1, NMX2, NMY0, NMY1, NMY2, NACT, &
NCENT, MAPX2, MAPY2, MAPAXY, MAPCXY, &
MAPTRN, CG, CX, CY, ATRNX, ATRNY, ITIME
USE W3IDATMD, ONLY: FLCUR
USE W3ODATMD, ONLY: NDSE, NDST, FLBPI, NBI, TBPI0, TBPIN, &
ISBPI, BBPI0, BBPIN
#ifdef W3_S
USE W3SERVMD, ONLY: STRACE
#endif
!/
IMPLICIT NONE
!/
!/ ------------------------------------------------------------------- /
!/ Parameter list
!/
INTEGER, INTENT(IN) :: ISEA, MAPSTA(NY*NX), MAPFS(NY*NX)
REAL, INTENT(IN) :: DTG, VGX, VGY
REAL, INTENT(INOUT) :: CFLXYMAX
!/
!/ ------------------------------------------------------------------- /
!/ Local parameters
!/
INTEGER :: ITH, IK, IXY, IP
INTEGER :: IX, IY, IXC, IYC, IBI
#ifdef W3_S
INTEGER, SAVE :: IENT = 0
#endif
REAL :: CG0, CGA, CGN, CGX, CGY, CXC, CYC, &
CXMIN, CXMAX, CYMIN, CYMAX
REAL :: CGC, FGSE = 1.
REAL :: FTH, FTHX, FTHY, FCG, FCGX, FCGY
REAL :: CP, CQ
!/
!/ Automatic work arrays
!/
REAL :: VLCFLX, VLCFLY
REAL :: CXTOT, CYTOT
!/
!/ ------------------------------------------------------------------- /
!/
#ifdef W3_S
CALL STRACE (IENT, 'W3XYCFL')
#endif
!
! 1. Preparations --------------------------------------------------- *
! 1.a Set constants
!
!
CFLXYMAX=0.
IX = MAPSF(ISEA,1)
IY = MAPSF(ISEA,2)
IXY = MAPSF(ISEA,3)
DO IK=1,NK
DO ITH=1,NTH
CXTOT = ECOS(ITH) * CG(IK,ISEA) / CLATS(ISEA)
CYTOT = ESIN(ITH) * CG(IK,ISEA)
#ifdef W3_MGP
CXTOT = CXTOT - VGX/CLATS(ISEA)
CYTOT = CYTOT - VGY
#endif
IF ( FLCUR ) THEN
CXTOT = CXTOT + CX(ISEA)/CLATS(ISEA)
CYTOT = CYTOT + CY(ISEA)
END IF
CP = CXTOT*DPDX(IY,IX) + CYTOT*DPDY(IY,IX)
CQ = CXTOT*DQDX(IY,IX) + CYTOT*DQDY(IY,IX)
VLCFLX = CP*DTG
VLCFLY = CQ*DTG
CFLXYMAX = MAX(VLCFLX,VLCFLY,CFLXYMAX)
END DO
END DO
RETURN
!/
!/ End of W3XYCFL ----------------------------------------------------- /
!/
END SUBROUTINE W3CFLXY
!/
!/ End of module W3PRO3MD -------------------------------------------- /
!/
END MODULE W3PRO3MD