#include "w3macros.h" !/ ------------------------------------------------------------------- / MODULE W3PRO2MD !/ !/ +-----------------------------------+ !/ | WAVEWATCH III NOAA/NCEP | !/ | H. L. Tolman | !/ | FORTRAN 90 | !/ | Last update : 29-May-2014 | !/ +-----------------------------------+ !/ !/ 04-Feb-2000 : Origination. ( version 2.00 ) !/ 24-Jan-2001 : Flat grid version ( version 2.06 ) !/ 08-Feb-2001 : UQ routines moved to own module ( version 2.08 ) !/ 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 obstacles added. ( 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 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 ) !/ 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) !/ 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 ) !/ !/ 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 porpagation scheme in single ! module. ! ! 2. Variables and types : ! ! Name Type Scope Description ! ---------------------------------------------------------------- ! TRNMIN R.P. Private Minimum transparancy for local ! ---------------------------------------------------------------- ! ! 3. Subroutines and functions : ! ! Name Type Scope Description ! ---------------------------------------------------------------- ! W3MAP2 Subr. Public Set up auxiliary maps. ! W3XYP2 Subr. Public Third order spatial propagation. ! W3KTP2 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 : ! ! 6. Switches : ! ! !/UQ 3rd order UQ propagation scheme. ! !/UNO 2nd order UNO propagation scheme. ! ! !/MGP Correct for motion of grid. ! ! !/OMPH Hybrid OpenMP directives. ! ! !/C90 Cray FORTRAN 90 compiler directives. ! !/NEC NEC SXF90 compiler directives. ! ! !/TDYN, !/DSS0, !/XW0, !/XW1 ! Diffusion options in W3XYP2 ! ! !/S Enable subroutine tracing. ! !/Tn Enable test output. ! ! 7. Source code : ! !/ ------------------------------------------------------------------- / !/ !/ Public variables !/ PUBLIC !/ !/ Private data !/ REAL, PRIVATE, PARAMETER:: TRNMIN = 0.95 !/ CONTAINS !/ ------------------------------------------------------------------- / SUBROUTINE W3MAP2 !/ !/ !/ +-----------------------------------+ !/ | 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 !/S USE W3SERVMD, ONLY: STRACE !/ IMPLICIT NONE !/ !/ ------------------------------------------------------------------- / !/ Parameter list !/ !/ ------------------------------------------------------------------- / !/ Local parameters !/ INTEGER :: IX, IY, IXY0, IX2, IY2, IX0, IY0, & IK, ITH, ISP, ISP0, ISP2 !/S INTEGER, SAVE :: IENT = 0 !/T INTEGER :: MAPTXY(NY,NX), I, IXY !/T INTEGER :: MAPTST(NK+2,NTH) !/ !/ ------------------------------------------------------------------- / !/ !/S CALL STRACE (IENT, 'W3MAP2') ! ! 1. Map MAPX2 ------------------------------------------------------ * ! 1.a Range 1 to NMX0 ! !/T MAPTXY = 0. ! 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 !/T MAPTXY(IY,IX) = MAPTXY(IY,IX) + 1 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 !/T MAPTXY(IY,IX) = MAPTXY(IY,IX) + 2 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 !/T MAPTXY(IY,IX) = MAPTXY(IY,IX) + 4 END IF END DO END DO ! !/T WRITE (NDST,9000) 'MAPX2', NMX0, NMX1-NMX0, & !/T NMX2-NMX1, NMX2 !/T DO IY=NY, 1, -1 !/T WRITE (NDST,9001) (MAPTXY(IY,IX),IX=1, NX) !/T END DO ! ! 2. Map MAPY2 ------------------------------------------------------ * ! 2.a Range 1 to NMY0 ! !/T MAPTXY = 0. ! 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 !/T MAPTXY(IY,IX) = MAPTXY(IY,IX) + 1 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 !/T MAPTXY(IY,IX) = MAPTXY(IY,IX) + 2 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 !/T MAPTXY(IY,IX) = MAPTXY(IY,IX) + 4 END IF END DO END DO ! !/T WRITE (NDST,9000) 'MAPY2', NMY0, NMY1-NMY0, & !/T NMY2-NMY1, NMY2 !/T DO IY=NY, 1, -1 !/T WRITE (NDST,9001) (MAPTXY(IY,IX),IX=1, NX) !/T END DO ! ! 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 ! !/T MAPTST = 0 ! ! 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 !/T MAPTST(IK+1,ITH) = MAPTST(IK+1,ITH) + 1 END DO END DO ! !/T WRITE (NDST,9000) 'MAPTH2', ISP, 0, 0, ISP !/T DO IK=NK+2, 1, -1 !/T WRITE (NDST,9001) (MAPTST(IK,ITH),ITH=1, NTH) !/T END DO ! !/T MAPTST = 0 ! ! 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 !/T MAPTST(IK+1,ITH) = MAPTST(IK+1,ITH) + 1 END DO END DO ! DO ITH=1, NTH ISP0 = ISP0 + 1 ISP2 = NK+1 + (ITH-1)*(NK+2) MAPWN2(ISP0) = ISP2 !/T MAPTST(NK+1,ITH) = MAPTST(NK+1,ITH) + 2 END DO ! DO ITH=1, NTH ISP0 = ISP0 + 1 ISP2 = 1 + (ITH-1)*(NK+2) MAPWN2(ISP0) = ISP2 !/T MAPTST(1,ITH) = MAPTST(1,ITH) + 4 END DO ! !/T WRITE (NDST,9000) 'MAPWN2', NSPEC-NTH, NTH, NTH, NSPEC+NTH !/T DO IK=NK+2, 1, -1 !/T WRITE (NDST,9001) (MAPTST(IK,ITH),ITH=1, NTH) !/T END DO ! RETURN ! ! Formats ! !/T 9000 FORMAT (/' TEST W3MAP2 : TEST MAP FOR PROPAGATION'/ & !/T ' MAP : ',A/ & !/T ' CENTRAL : ',I6/ & !/T ' ABOVE : ',I6/ & !/T ' BELOW : ',I6/ & !/T ' TOTAL : ',I6/) !/T 9001 FORMAT (1X,130I1) ! !/T 9010 FORMAT (' TEST W3MAP2 : COMPOSITE MAPS TH2, WN2 AND BTK') !/T 9011 FORMAT (2X,60I2) !/ !/ End of W3MAP2 ----------------------------------------------------- / !/ END SUBROUTINE W3MAP2 !/ ------------------------------------------------------------------- / SUBROUTINE W3XYP2 ( ISP, DTG, MAPSTA, MAPFS, VQ, VGX, VGY ) !/ !/ +-----------------------------------+ !/ | 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. ! ! !/NEC Enable NEC SXF90 compiler 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 !/S USE W3SERVMD, ONLY: STRACE !/UQ USE W3UQCKMD !/UNO USE W3UNO2MD !/ 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 !/S INTEGER, SAVE :: IENT = 0 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 ------------------------------------------------------------------- !/ !/ ------------------------------------------------------------------- / !/ !/S CALL STRACE (IENT, 'W3XYP2') ! ! 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) !/MGP CGX = CGX - VGX !/MGP CGY = CGY - VGY ! 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) ) !/MGP CXC = MAX ( ABS(CXMIN-VGX) , ABS(CXMAX-VGX) ) !/MGP CYC = MAX ( ABS(CYMIN-VGY) , ABS(CYMAX-VGY) ) 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 ! YFIRST = MOD(ITIME,2) .EQ. 0 ! !/T WRITE (NDST,9000) YFIRST !/T WRITE (NDST,9001) ISP, ITH, IK, ECOS(ITH), ESIN(ITH) ! !/TDYN IF ( ISP .EQ. 1 ) DTME = DTME + DTG ! 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. !/T WRITE (NDST,9002) DFRR, CELLP, DTME !/T ELSE !/T WRITE (NDST,9003) END IF ! ! 1.b Initialize arrays ! !/T WRITE (NDST,9010) ! 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 ! !/T WRITE (NDST,9020) NSEA ! !/OMPH/!$OMP PARALLEL DO PRIVATE (ISEA, IXY) ! 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) !/MGP CXTOT(IXY) = CXTOT(IXY) - VGX/CLATS(ISEA) !/MGP CYTOT(IXY) = CYTOT(IXY) - VGY !/T1 IF ( .NOT. FLCUR ) & !/T1 WRITE (NDST,9021) ISEA, MAPSF(ISEA,1), MAPSF(ISEA,2), & !/T1 VQ(IXY), CXTOT(IXY), CYTOT(IXY) END DO ! !/OMPH/!$OMP END PARALLEL DO ! IF ( FLCUR ) THEN !/T WRITE (NDST,9022) DO ISEA=1, NSEA IXY = MAPSF(ISEA,3) CXTOT(IXY) = CXTOT(IXY) + CX(ISEA)/CLATS(ISEA) CYTOT(IXY) = CYTOT(IXY) + CY(ISEA) !/T1 WRITE (NDST,9021) ISEA, MAPSF(ISEA,1), MAPSF(ISEA,2), & !/T1 VQ(IXY), CXTOT(IXY), CYTOT(IXY) END DO END IF ! !/OMPH/!$OMP PARALLEL DO PRIVATE (ISEA, IX, IY, IXY, CP, CQ) ! 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 ! !/OMPH/!$OMP END PARALLEL DO ! ! 2.b Diffusion coefficients ! IF ( DTME .NE. 0. ) THEN ! !/OMPH/!$OMP PARALLEL DO PRIVATE (ISEA, IX, IY, IXY, & !/OMPH/!$OMP& DCELL, XWIND, TFAC, DSS, DNN) ! 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 ) ) !/XW0 XWIND = 0. !/XW1 XWIND = 1. TFAC = MIN ( 1. , (CLATS(ISEA)/CLATMN)**2 ) DSS = XWIND * DCELL + (1.-XWIND) * DSSD * TFAC !/DSS0 DSS = 0. 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 ! !/OMPH/!$OMP END PARALLEL DO ! END IF ! ! 3. Loop over sub-steps -------------------------------------------- * ! DO ITLOC=1, NTLOC ! ! 3.a Propagate fields ! !/OMPH/!$OMP PARALLEL DO PRIVATE (ISEA, IX, IY, IXY ) ! 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 ! !/OMPH/!$OMP END PARALLEL DO ! IF ( YFIRST ) THEN ! !/UQ IF ( FLCY ) CALL W3QCK3 & !/UQ (NX, NY, NX, NY, VLCFLY, ATRNY, VQ, & !/UQ .FALSE., 1, MAPAXY, NACT, MAPY2, NMY0, & !/UQ NMY1, NMY2, NDSE, NDST ) !/UQ IF ( FLCX ) CALL W3QCK3 & !/UQ (NX, NY, NX, NY, VLCFLX, ATRNX, VQ, & !/UQ GLOBAL, NY, MAPAXY, NACT, MAPX2, NMX0, & !/UQ NMX1, NMX2, NDSE, NDST ) ! !/UNO IF ( FLCY ) CALL W3UNO2s & !/UNO (NX, NY, NX, NY, VLCFLY, ATRNY, VQ, & !/UNO .FALSE., 1, MAPAXY, NACT, MAPY2, NMY0, & !/UNO NMY1, NMY2, NDSE, NDST ) !/UNO IF ( FLCX ) CALL W3UNO2s & !/UNO (NX, NY, NX, NY, VLCFLX, ATRNX, VQ, & !/UNO GLOBAL, NY, MAPAXY, NACT, MAPX2, NMX0, & !/UNO NMX1, NMX2, NDSE, NDST ) ! ELSE ! !/UQ IF ( FLCX ) CALL W3QCK3 & !/UQ (NX, NY, NX, NY, VLCFLX, ATRNX, VQ, & !/UQ GLOBAL, NY, MAPAXY, NACT, MAPX2, NMX0, & !/UQ NMX1, NMX2, NDSE, NDST ) !/UQ IF ( FLCY ) CALL W3QCK3 & !/UQ (NX, NY, NX, NY, VLCFLY, ATRNY, VQ, & !/UQ .FALSE., 1, MAPAXY, NACT, MAPY2, NMY0, & !/UQ NMY1, NMY2, NDSE, NDST ) ! !/UNO IF ( FLCX ) CALL W3UNO2s & !/UNO (NX, NY, NX, NY, VLCFLX, ATRNX, VQ, & !/UNO GLOBAL, NY, MAPAXY, NACT, MAPX2, NMX0, & !/UNO NMX1, NMX2, NDSE, NDST ) !/UNO IF ( FLCY ) CALL W3UNO2s & !/UNO (NX, NY, NX, NY, VLCFLY, ATRNY, VQ, & !/UNO .FALSE., 1, MAPAXY, NACT, MAPY2, NMY0, & !/UNO NMY1, NMY2, NDSE, NDST ) ! END IF ! !/OMPH/!$OMP PARALLEL DO PRIVATE (ISEA, IX, IY, IXY ) ! 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 ! !/OMPH/!$OMP END PARALLEL DO ! ! 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 ! !/OMPH/!$OMP PARALLEL DO PRIVATE (ISEA, IX, IY, IXY, & !/OMPH/!$OMP& QXX, QYY, QXY, DVQ ) ! 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 ! !/OMPH/!$OMP END PARALLEL DO ! END IF ! YFIRST = .NOT. YFIRST END DO ! ! 4. Store results in VQ in proper format --------------------------- * ! !/T WRITE (NDST,9040) NSEA ! !/C90/!DIR$ IVDEP !/NEC/!CDIR NODEP ! !/OMPH/!$OMP PARALLEL DO PRIVATE (ISEA, IXY ) ! DO ISEA=1, NSEA IXY = MAPSF(ISEA,3) IF ( MAPSTA(IXY) .GT. 0 ) THEN !/T2 WRITE (NDST,9041) ISEA, MAPSF(ISEA,1), MAPSF(ISEA,2), VQ(IXY) VQ(IXY) = MAX ( 0. , CG(IK,ISEA) / CLATS(ISEA) * VQ(IXY) ) ! ELSE ! VQ(IXY) = 0. END IF END DO ! !/OMPH/!$OMP END PARALLEL DO ! RETURN ! ! Formats ! !/T 9000 FORMAT (' TEST W3XYP2 : YFIRST :',L2) !/T 9001 FORMAT (' TEST W3XYP2 : ISP, ITH, IK, COS-SIN :',I8,2I4,2F7.3) !/T 9002 FORMAT (' TEST W3XYP2 : DFRR, CELLP, DTME :',3E10.3) !/T 9003 FORMAT (' TEST W3XYP2 : NO DISPERSION CORRECTION ') ! !/T 9010 FORMAT (' TEST W3XYP2 : INITIALIZE ARRAYS') ! !/T 9020 FORMAT (' TEST W3XYP2 : CALCULATING LCFLX/Y AND DSS/NN (NSEA=', & !/T I6,')') !/T1 9021 FORMAT (1X,I6,2I5,E12.4,2f7.3) !/T 9022 FORMAT (' TEST W3XYP2 : CORRECTING FOR CURRENT') ! !/T 9040 FORMAT (' TEST W3XYP2 : FIELD AFTER PROP. (NSEA=',I6,')') !/T2 9041 FORMAT (1X,I6,2I5,E12.4) !/ !/ End of W3XYP2 ----------------------------------------------------- / !/ END SUBROUTINE W3XYP2 !/ !/ ------------------------------------------------------------------- / SUBROUTINE W3KTP2 ( ISEA, FACTH, FACK, CTHG0, CG, WN, DEPTH, & DDDX, DDDY, CX, CY, DCXDX, DCXDY, & DCYDX, DCYDY, DCDX, DCDY, VA ) !/ !/ +-----------------------------------+ !/ | 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 !/S USE W3SERVMD, ONLY: STRACE !/UQ USE W3UQCKMD !/UNO USE W3UNO2MD !/ 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 !/S INTEGER, SAVE :: IENT = 0 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) !/ !/ ------------------------------------------------------------------- / !/ !/S CALL STRACE (IENT, 'W3KTP2') ! ! 1. Preparations --------------------------------------------------- * ! 1.a Initialize arrays ! IF ( FLCK ) VQ = 0. ! !/T WRITE (NDST,9000) FLCTH, FLCK, FACTH, FACK, CTMAX !/T WRITE (NDST,9010) ISEA, DEPTH, CX, CY, DDDX, DDDY, & !/T DCXDX, DCXDY, DCYDX, DCYDY ! ! 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 ! !/T WRITE (NDST,9020) !/T DO IK=1, NK+1 !/T WRITE (NDST,9021) IK, TPI/SIG(IK), TPI/WN(IK), & !/T CG(IK), DSDD(IK) !/T END DO ! ! 2.b Extract spectrum ! DO ISP=1, NSPEC VQ(MAPTH2(ISP)) = VA(ISP) END DO ! ! 3. Refraction velocities ------------------------------------------ * ! IF ( FLCTH ) THEN ! ! 3.a Set slope filter for depth refraction ! FDDMAX = 0. FDG = FACTH * CTHG0 ! DO ITH=1, NTH/2 FDDMAX = MAX(FDDMAX,ABS(ESIN(ITH)*DDDX-ECOS(ITH)*DDDY)) END DO ! DO IK=1, NK FRK(IK) = FACTH * DSDD(IK) / WN(IK) FRK(IK) = FRK(IK) / MAX ( 1. , FRK(IK)*FDDMAX/CTMAX ) FRG(IK) = FDG * CG(IK) END DO ! ! 3.b Depth refraction and great-circle propagation ! DO ISP=1, NSPEC VCFLT(MAPTH2(ISP)) = FRG(MAPWN(ISP)) * ECOS(ISP) & + FRK(MAPWN(ISP)) * ( ESIN(ISP)*DDDX - ECOS(ISP)*DDDY ) END DO ! !/REFRX! 3.c @C/@x refraction and great-circle propagation !/REFRX VCFLT = 0. !/REFRX FRK = 0. !/REFRX FDDMAX = 0. ! !/REFRX DO ISP=1, NSPEC !/REFRX FDDMAX = MAX ( FDDMAX , ABS ( & !/REFRX ESIN(ISP)*DCDX(MAPWN(ISP)) - ECOS(ISP)*DCDY(MAPWN(ISP)) ) ) !/REFRX END DO ! !/REFRX DO IK=1, NK !/REFRX FRK(IK) = FACTH * CG(IK) * WN(IK) / SIG(IK) !/REFRX END DO !/REFRX DO ISP=1, NSPEC !/REFRX VCFLT(MAPTH2(ISP)) = FRG(MAPWN(ISP)) * ECOS(ISP) & !/REFRX + FRK(MAPWN(ISP)) * ( ESIN(ISP)*DCDX(MAPWN(ISP)) & !/REFRX - ECOS(ISP)*DCDY(MAPWN(ISP)) ) !/REFRX END DO ! ! 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 ! !/UQ CALL W3QCK2 ( NTH, NK2, NTH, NK2, CFLK, FACK, DB, DM, & !/UQ VQ, .FALSE., 1, MAPTH2, NSPEC, & !/UQ MAPWN2, NSPEC-NTH, NSPEC, NSPEC+NTH, & !/UQ NDSE, NDST ) ! !/UNO CALL W3UNO2 ( NTH, NK2, NTH, NK2, CFLK, FACK, DB, DM, & !/UNO VQ, .FALSE., 1, MAPTH2, NSPEC, & !/UNO MAPWN2, NSPEC-NTH, NSPEC, NSPEC+NTH, & !/UNO NDSE, NDST ) END IF IF ( FLCTH ) THEN ! !/UQ CALL W3QCK1 ( NTH, NK2, NTH, NK2, VCFLT, VQ, .TRUE., & !/UQ NK2, MAPTH2, NSPEC, MAPTH2, NSPEC, NSPEC, & !/UQ NSPEC, NDSE, NDST ) ! !/UNO CALL W3UNO2r( NTH, NK2, NTH, NK2, VCFLT, VQ, .TRUE., & !/UNO NK2, MAPTH2, NSPEC, MAPTH2, NSPEC, NSPEC,& !/UNO NSPEC, NDSE, NDST ) ! END IF ELSE IF ( FLCTH ) THEN ! !/UQ CALL W3QCK1 ( NTH, NK2, NTH, NK2, VCFLT, VQ, .TRUE., & !/UQ NK2, MAPTH2, NSPEC, MAPTH2, NSPEC, NSPEC, & !/UQ NSPEC, NDSE, NDST ) ! !/UNO CALL W3UNO2r( NTH, NK2, NTH, NK2, VCFLT, VQ, .TRUE., & !/UNO NK2, MAPTH2, NSPEC, MAPTH2, NSPEC, NSPEC,& !/UNO NSPEC, NDSE, NDST ) ! END IF IF ( FLCK ) THEN DO ITH=1, NTH VQ(NK+2+(ITH-1)*NK2) = FACHFA * VQ(NK+1+(ITH-1)*NK2) END DO ! !/UQ CALL W3QCK2 ( NTH, NK2, NTH, NK2, CFLK, FACK, DB, DM, & !/UQ VQ, .FALSE., 1, MAPTH2, NSPEC, & !/UQ MAPWN2, NSPEC-NTH, NSPEC, NSPEC+NTH, & !/UQ NDSE, NDST ) ! !/UNO CALL W3UNO2 ( NTH, NK2, NTH, NK2, CFLK, FACK, DB, DM, & !/UNO VQ, .FALSE., 1, MAPTH2, NSPEC, & !/UNO MAPWN2, NSPEC-NTH, NSPEC, NSPEC+NTH, & !/UNO NDSE, NDST ) ! END IF END IF ! ! 6. Store reults --------------------------------------------------- * ! DO ISP=1, NSPEC VA(ISP) = VQ(MAPTH2(ISP)) END DO ! RETURN ! ! Formats ! !/T 9000 FORMAT ( ' TEST W3KTP2 : FLCTH-K, FACTH-K, CTMAX :', & !/T 2L2,2E10.3,F7.3) !/T 9010 FORMAT ( ' TEST W3KTP2 : LOCAL DATA :',I7,F7.1,2F6.2,1X,6E10.2) !/T 9020 FORMAT ( ' TEST W3KTP2 : IK, T, L, CG, DSDD : ') !/T 9021 FORMAT ( ' ',I3,F7.2,F7.1,F7.2,E11.3) ! !/T0 9040 FORMAT (/' TEST W3KTP2 : NORMALIZED ',A/) !/T0 9041 FORMAT (1X,60(1X,I2)) !/T0 9042 FORMAT (1X,60I3) !/ !/ End of W3KTP2 ----------------------------------------------------- / !/ END SUBROUTINE W3KTP2 !/ !/ End of module W3PRO2MD -------------------------------------------- / !/ END MODULE W3PRO2MD