!> @file !> @brief Source term integration routine. !> !> @author H. L. Tolman !> @author F. Ardhuin !> @date 22-Mar-2021 !> #include "w3macros.h" !/ ------------------------------------------------------------------- / !> !> @brief Source term integration routine. !> !> @author H. L. Tolman !> @author F. Ardhuin !> @date 22-Mar-2021 !> !> @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. !> MODULE W3SRCEMD !/ !/ +-----------------------------------+ !/ | WAVEWATCH III NOAA/NCEP | !/ | H. L. Tolman | !/ | F. Ardhuin | !/ | FORTRAN 90 | !/ | Last update : 22-Mar-2021 | !/ +-----------------------------------+ !/ !/ For updates see subroutine. !/ ! 1. Purpose : ! ! Source term integration routine. ! ! 2. Variables and types : ! ! Name Type Scope Description ! ---------------------------------------------------------------- ! OFFSET R.P. Private Offset in time integration scheme. ! 0.5 in original WAM, now 1.0 ! ---------------------------------------------------------------- ! ! 3. Subroutines and functions : ! ! Name Type Scope Description ! ---------------------------------------------------------------- ! W3SRCE Subr. Public Calculate and integrate source terms. ! ---------------------------------------------------------------- ! ! 4. Subroutines and functions used : ! ! See corresponding documentation of W3SRCE. ! ! 5. Remarks : ! ! 6. Switches : ! ! See section 9 of W3SRCE. ! ! 7. Source code : ! !/ ------------------------------------------------------------------- / !/ REAL, PARAMETER, PRIVATE:: OFFSET = 1. !/ CONTAINS !/ ------------------------------------------------------------------- / !> !> @brief Calculate and integrate source terms for a single grid point. !> !> @verbatim !> Physics : see manual and corresponding subroutines. !> !> Numerics : !> !> Dynamic-implicit integration of the source terms based on !> WW-II (Tolman 1992). The dynamic time step is calculated !> given a maximum allowed change of spectral densities for !> frequencies / wavenumbers below the usual cut-off. !> The maximum change is given by the minimum of a parametric !> and a relative change. The parametric change relates to a !> PM type equilibrium range !> !> -1 (2pi)**4 1 !> dN(k) = Xp alpha pi ---------- ------------ !> max g**2 k**3 sigma !> !> 1 . !> = FACP ------------ (1) !> k**3 sigma . !> !> where !> alpha = 0.62e-4 (set in W3GRID) !> Xp fraction of PM shape (read in W3GRID) !> FACP combined factor (set in W3GRID) !> !> The maximum relative change is given as !> !> / +- -+ \ . !> dN(k) = Xr max | N(k) , max | Nx , Xfilt N(k) | | (2) !> max \ +- max-+ / . !> !> where !> Xr fraction of relative change (read in W3GRID) !> Xfilt filter level (read in W3GRID) !> Nx Maximum parametric change (1) !> for largest wavenumber. !> @endverbatim !> !> @param[in] srce_call !> @param[in] IT !> @param[in] ISEA !> @param[in] JSEA !> @param[in] IX Discrete grid point counters. !> @param[in] IY Discrete grid point counters. !> @param[in] IMOD Model number. !> @param[in] SPECOLD !> @param[inout] SPEC Spectrum (action) in 1-D form. !> @param[out] VSIO !> @param[out] VDIO !> @param[out] SHAVEIO !> @param[inout] ALPHA Nondimensional 1-D spectrum corresponding !> to above full spectra (Phillip's const.). !> @param[inout] WN1 Discrete wavenumbers. !> @param[inout] CG1 Id. group velocities. !> @param[in] CLATSL !> @param[in] D_INP Depth, compared to DMIN to get DEPTH. !> @param[in] U10ABS Wind speed at reference height. !> @param[in] U10DIR Id. wind direction. !> @param[inout] TAUA Magnitude of total atmospheric stress. !> @param[inout] TAUADIR Direction of atmospheric stress. !> @param[in] AS Air-sea temperature difference. !> @param[inout] USTAR Friction velocity. !> @param[inout] USTDIR Idem, direction. !> @param[in] CX Current velocity component. !> @param[in] CY Current velocity component. !> @param[in] ICE Sea ice concentration. !> @param[in] ICEH Sea ice thickness. !> @param[inout] ICEF Sea ice floe diameter. !> @param[in] ICEDMAX Sea ice maximum floe diameter !> @param[in] REFLEC Reflection coefficients. !> @param[in] REFLED Reflection direction. !> @param[in] DELX Grid cell size in X direction. !> @param[in] DELY Grid cell size in Y direction. !> @param[in] DELA Grid cell area. !> @param[in] TRNX Grid transparency in X. !> @param[in] TRNY Grid transparency in Y. !> @param[in] BERG Iceberg damping coefficient. !> @param[inout] FPI Peak-input frequency. !> @param[out] DTDYN Average dynamic time step. !> @param[out] FCUT Cut-off frequency for tail. !> @param[in] DTG Global time step. !> @param[inout] TAUWX !> @param[inout] TAUWY !> @param[inout] TAUOX !> @param[inout] TAUWIX !> @param[inout] TAUWIY !> @param[inout] TAUWNX !> @param[inout] TAUWNY !> @param[inout] PHIAW !> @param[inout] CHARN !> @param[inout] TWS !> @param[inout] PHIOC !> @param[inout] WHITECAP Whitecap statistics. !> @param[in] D50 Sand grain size. !> @param[in] PSIC Critical shields. !> @param[inout] BEDFORM Bedform parameters. !> @param[inout] PHIBBL Energy flux to BBL. !> @param[inout] TAUBBL Momentum flux to BBL. !> @param[inout] TAUICE Momentum flux to sea ice. !> @param[inout] PHICE Energy flux to sea ice. !> @param[inout] TAUOCX Total ocean momentum component. !> @param[inout] TAUOCY Total ocean momentum component. !> @param[inout] WNMEAN Mean wave number. !> @param[in] DAIR Air density. !> @param[in] COEF !> !> @author H. L. Tolman !> @author F. Ardhuin !> @author A. Roland !> @author M. Dutour Sikiric !> @date 22-Mar-2021 !> SUBROUTINE W3SRCE ( srce_call, IT, ISEA, JSEA, IX, IY, IMOD, & SPECOLD, SPEC, VSIO, VDIO, SHAVEIO, & ALPHA, WN1, CG1, CLATSL, & D_INP, U10ABS, U10DIR, & #ifdef W3_FLX5 TAUA, TAUADIR, & #endif AS, USTAR, USTDIR, & CX, CY, ICE, ICEH, ICEF, ICEDMAX, & REFLEC, REFLED, DELX, DELY, DELA, TRNX, & TRNY, BERG, FPI, DTDYN, FCUT, DTG, TAUWX, & TAUWY, TAUOX, TAUOY, TAUWIX, TAUWIY, TAUWNX,& TAUWNY, PHIAW, CHARN, TWS, PHIOC, WHITECAP, & D50, PSIC, BEDFORM , PHIBBL, TAUBBL, TAUICE,& PHICE, TAUOCX, TAUOCY, WNMEAN, DAIR, COEF) !/ !/ +-----------------------------------+ !/ | WAVEWATCH III NOAA/NCEP | !/ | H. L. Tolman | !/ | F. Ardhuin | !/ | A. Roland | !/ | M. Dutour Sikiric | !/ | FORTRAN 90 | !/ | Last update : 22-Mar-2021 | !/ +-----------------------------------+ !/ !/ 06-Dec-1996 : Final FORTRAN 77 ( version 1.18 ) !/ 04-Feb-2000 : Upgrade to FORTRAN 90 ( version 2.00 ) !/ 14-Feb-2000 : Exact-NL added ( version 2.01 ) !/ 04-May-2000 : Non-central integration ( version 2.03 ) !/ 02-Feb-2001 : Xnl version 3.0 ( version 2.07 ) !/ 09-May-2002 : Switch clean up. ( version 2.21 ) !/ 13-Nov-2002 : Add stress vector. ( version 3.00 ) !/ 27-Nov-2002 : First version of VDIA and MDIA. ( version 3.01 ) !/ 07-Oct-2003 : Output options for NN training. ( version 3.05 ) !/ 24-Dec-2004 : Multiple model version. ( version 3.06 ) !/ 23-Jun-2006 : Linear input added. ( version 3.09 ) !/ 27-Jun-2006 : Adding file name preamble. ( version 3.09 ) !/ 04-Jul-2006 : Separation of stress computation. ( version 3.09 ) !/ 16-Apr-2007 : Miche style limiter added. ( version 3.11 ) !/ (J. H. Alves) !/ 25-Apr-2007 : Battjes-Janssen Sdb added. ( version 3.11 ) !/ (J. H. Alves) !/ 09-Oct-2007 : Adding WAM 4+ and SB1 options. ( version 3.13 ) !/ (F. Ardhuin) !/ 29-May-2009 : Preparing distribution version. ( version 3.14 ) !/ 19-Aug-2010 : Making treatment of 0 water depth ( version 3.14.6 ) !/ consistent with the rest of the model. !/ 31-Mar-2010 : Adding ice conc. and reflections ( version 3.14.4 ) !/ 15-May-2010 : Adding transparencies ( version 3.14.4 ) !/ 01-Jun-2011 : Movable bed bottom friction in BT4 ( version 4.01 ) !/ 01-Jul-2011 : Energy and momentum flux, friction ( version 4.01 ) !/ 24-Aug-2011 : Uses true depth for depth-induced ( version 4.04 ) !/ 16-Sep-2011 : Initialization of TAUWAX, TAUWAY ( version 4.04 ) !/ 1-Dec-2011 : Adding BYDRZ source term package ( version 4.04 ) !/ ST6 and optional Hwang (2011) !/ stresses FLX4. !/ 14-Mar-2012 : Update of BT4, passing PSIC ( version 4.04 ) !/ 13-Jul-2012 : Move GMD (SNL3) and nonlinear filter (SNLS) !/ from 3.15 (HLT). ( version 4.08 ) !/ 28-Aug-2013 : Corrected MLIM application ( version 4.11 ) !/ 10-Sep-2013 : Special treatment for IG band ( version 4.15 ) !/ 14-Nov-2013 : Make orphaned pars in par lst local ( version 4.13 ) !/ 17-Nov-2013 : Coupling fraction of ice-free ( version 4.13 ) !/ surface to SIN and SDS. (S. Zieger) !/ 01-Avr-2014 : Adding ice thickness and floe size ( version 4.18 ) !/ 23-May-2014 : Adding ice fluxes to W3SRCE ( version 5.01 ) !/ 27-Aug-2015 : Adding inputs to function W3SIS2 ( version 5.10 ) !/ 13-Dec-2015 : Implicit integration of Sice (F.A.) ( version 5.10 ) !/ 30-Jul-2017 : Adds TWS in interface ( version 6.04 ) !/ 07-Jan-2018 : Allows variable ice scaling (F.A.) ( version 6.04 ) !/ 01-Jan-2018 : Add implicit source term integration ( version 6.04) !/ 01-Jan-2018 : within PDLIB (A. Roland, M. Dutour !/ 18-Aug-2018 : S_{ice} IC5 (Q. Liu) ( version 6.06) !/ 26-Aug-2018 : UOST (Mentaschi et al. 2015, 2018) ( version 6.06 ) !/ 22-Mar-2021 : Add extra fields used in coupling ( version 7.13 ) !/ 07-Jun-2021 : S_{nl5} GKE NL5 (Q. Liu) ( version 7.13 ) !/ 19-Jul-2021 : Momentum and air density support ( version 7.14 ) !/ !/ Copyright 2009-2013 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 : ! ! Calculate and integrate source terms for a single grid point. ! ! 2. Method : ! ! Physics : see manual and corresponding subroutines. ! ! Numerics : ! ! Dynamic-implicit integration of the source terms based on ! WW-II (Tolman 1992). The dynamic time step is calculated ! given a maximum allowed change of spectral densities for ! frequencies / wavenumbers below the usual cut-off. ! The maximum change is given by the minimum of a parametric ! and a relative change. The parametric change relates to a ! PM type equilibrium range ! ! -1 (2pi)**4 1 ! dN(k) = Xp alpha pi ---------- ------------ ! max g**2 k**3 sigma ! ! 1 . ! = FACP ------------ (1) ! k**3 sigma . ! ! where ! alpha = 0.62e-4 (set in W3GRID) ! Xp fraction of PM shape (read in W3GRID) ! FACP combined factor (set in W3GRID) ! ! The maximum relative change is given as ! ! / +- -+ \ . ! dN(k) = Xr max | N(k) , max | Nx , Xfilt N(k) | | (2) ! max \ +- max-+ / . ! ! where ! Xr fraction of relative change (read in W3GRID) ! Xfilt filter level (read in W3GRID) ! Nx Maximum parametric change (1) ! for largest wavenumber. ! ! 3. Parameters : ! ! Parameter list ! ---------------------------------------------------------------- ! IX,IY Int. I Discrete grid point counters. ! IMOD Int. I Model number. ! SPEC R.A. I/O Spectrum (action) in 1-D form. ! ALPHA R.A. I/O Nondimenional 1-D spectrum corresponding ! to above full spectra (Phillip's const.). ! Calculated separately for numerical ! economy on vector machine (W3SPR2). ! WN1 R.A. I Discrete wavenumbers. ! CG1 R.A. I Id. group velocities. ! D_INP Real. I Depth. Compared to DMIN to get DEPTH. ! U10ABS Real. I Wind speed at reference height. ! U10DIR Real. I Id. wind direction. ! TAUA Real. I Magnitude of total atmospheric stress ( !/FLX5 ) ! TAUADIR Real. I Direction of atmospheric stress ( !/FLX5 ) ! AS Real. I Air-sea temp. difference. ( !/ST3 ) ! USTAR Real. !/O Friction velocity. ! USTDIR Real !/O Idem, direction. ! CX-Y Real. I Current velocity components. ( !/BS1 ) ! ICE Real I Sea ice concentration ! ICEH Real I Sea ice thickness ! ICEF Real I/O Sea ice maximum floe diameter (updated) ! ICEDMAX Real I/O Sea ice maximum floe diameter ! BERG Real I Iceberg damping coefficient ( !/BS1 ) ! REFLEC R.A. I reflection coefficients ( !/BS1 ) ! REFLED I.A. I reflection direction ( !/BS1 ) ! TRNX-Y Real I Grid transparency in X and Y ( !/BS1 ) ! DELX Real. I grid cell size in X direction ( !/BS1 ) ! DELY Real. I grid cell size in Y direction ( !/BS1 ) ! DELA Real. I grid cell area ( !/BS1 ) ! FPI Real I/O Peak-input frequency. ( !/ST2 ) ! WHITECAP R.A. O Whitecap statisics ( !/ST4 ) ! DTDYN Real O Average dynamic time step. ! FCUT Real O Cut-off frequency for tail. ! DTG Real I Global time step. ! D50 Real I Sand grain size ( !/BT4 ) ! BEDFORM R.A. I/O Bedform parameters ( !/BT4 ) ! PSIC Real I Critical Shields ( !/BT4 ) ! PHIBBL Real O Energy flux to BBL ( !/BTx ) ! TAUBBL R.A. O Momentum flux to BBL ( !/BTx ) ! TAUICE R.A. O Momentum flux to sea ice ( !/ICx ) ! PHICE Real O Energy flux to sea ice ( !/ICx ) ! TAUOCX-YReal O Total ocean momentum components ! WNMEAN Real O Mean wave number ! DAIR Real I Air density ! ---------------------------------------------------------------- ! Note: several pars are set to I/O to avoid compiler warnings. ! ! 4. Subroutines used : ! ! Name Type Module Description ! ---------------------------------------------------------------- ! W3SPRn Subr. W3SRCnMD Mean wave parameters for use in ! source terms. ! W3FLXn Subr. W3FLXnMD Flux/stress computation. ! W3SLNn Subr. W3SLNnMD Linear input. ! W3SINn Subr. W3SRCnMD Input source term. ! W3SNLn Subr. W3SNLnMD Nonlinear interactions. ! W3SNLS Subr. W3SNLSMD Nonlinear smoother. ! W3SDSn Subr. W3SRCnMD Whitecapping source term ! W3SBTn Subr. W3SBTnMD Bottom friction source term. ! W3SDBn Subr. W3SBTnMD Depth induced breaking source term. ! W3STRn Subr. W3STRnMD Triad interaction source term. ! W3SBSn Subr. W3SBSnMD Bottom scattering source term. ! W3REFn Subr. W3REFnMD Reflexions (shore, icebergs ...). ! STRACE Subr. W3SERVMD Subroutine tracing (!/S) ! ---------------------------------------------------------------- ! ! 5. Called by : ! ! Name Type Module Description ! ---------------------------------------------------------------- ! W3WAVE Subr. W3WAVEMD Actual wave model routine. ! ---------------------------------------------------------------- ! ! 6. Error messages : ! ! None. ! ! 7. Remarks : ! ! - No testing is performed on the status of the grid point. ! ! 8. Structure : ! ! ----------------------------------------------------------------- ! 1. Preparations ! a Set maximum change and wavenumber arrays. ! b Prepare dynamic time stepping. ! c Compute mean parameters. ( W3SPRn ) ! d Compute stresses (if posible). ! e Prepare cut-off ! f Test output for !/NNT option. ! --start-dynamic-integration-loop--------------------------------- ! 2. Calculate source terms ! a Input. ( W3SLNx, W3SINn ) ! b Nonlinear interactions. ( W3SNLn ) ! c Dissipation ( W3SDSn ) ! 1 as included in source terms ( W3SDSn ) ! 2 optional dissipation due to different physics ( W3SWLn ) ! d Bottom friction. ( W3SBTn ) ! 3. Calculate cut-off frequencie(s) ! 4. Summation of source terms and diagonal term and time step. ! 5. Increment spectrum. ! 6. Add tail ! a Mean wave parameters and cut-off ( W3SPRn ) ! b 'Seeding' of spectrum. ( !/SEED ) ! c Add tail ! 7. Check if integration complete. ! --end-dynamic-integration-loop----------------------------------- ! 8. Save integration data. ! ----------------------------------------------------------------- ! ! 9. Switches : ! ! !/FLX1 Wu (1980) stress computation. ( Choose one ) ! !/FLX2 T&C (1996) stress computation. ! !/FLX3 T&C (1996) stress computation with cap. ! !/FLX4 Hwang (2011) stress computation (2nd order). ! !/FLX5 Direct use of stress from atmoshperic model. ! ! !/LN0 No linear input. ( Choose one ) ! ! !/ST0 No input and dissipation. ( Choose one ) ! !/ST1 WAM-3 input and dissipation. ! !/ST2 Tolman and Chalikov (1996) input and dissipation. ! !/ST3 WAM 4+ input and dissipation. ! !/ST4 Ardhuin et al. (2009, 2010) ! !/ST6 BYDB source terms after Babanin, Young, Donelan and Banner. ! ! !/NL0 No nonlinear interactions. ( Choose one ) ! !/NL1 Discrete interaction approximation. ! !/NL2 Exact nonlinear interactions. ! !/NL3 Generalized Multiple DIA. ! !/NL4 Two Scale Approximation ! !/NL5 Generalized Kinetic Equation. ! !/NLS Nonlinear HF smoother. ! ! !/BT0 No bottom friction. ( Choose one ) ! !/BT1 JONSWAP bottom friction. ! !/BT4 Bottom friction using movable bed roughness ! (Tolman 1994, Ardhuin & al. 2003) ! !/BT8 Muddy bed (Dalrymple & Liu). ! !/BT9 Muddy bed (Ng). ! ! !/IC1 Dissipation via interaction with ice according to simple ! methods: 1) uniform in frequency or ! !/IC2 2) Liu et al. model ! !/IC3 Dissipation via interaction with ice according to a ! viscoelastic sea ice model (Wang and Shen 2010). ! !/IC4 Dissipation via interaction with ice as a function of freq. ! (empirical/parametric methods) ! !/IC5 Dissipation via interaction with ice according to a ! viscoelastic sea ice model (Mosig et al. 2015). ! !/DB0 No depth-limited breaking. ( Choose one ) ! !/DB1 Battjes-Janssen depth-limited breaking. ! ! !/TR0 No triad interactions. ( Choose one ) ! !/TR1 Lumped Triad Approximation (LTA). ! ! !/BS0 No bottom scattering. ( Choose one ) ! !/BS1 Scattering term by Ardhuin and Magne (2007). ! ! !/MLIM Miche style limiter for shallow water and steepness. ! ! !/SEED 'Seeding' of lowest frequency for suffuciently strong ! winds. ! ! !/NNT Write output to file test_data_NNN.ww3 for NN training. ! ! !/S Enable subroutine tracing. ! !/T Enable general test output. ! ! 10. Source code : ! !/ ------------------------------------------------------------------- / USE CONSTANTS, ONLY: DWAT, srce_imp_post, srce_imp_pre, & srce_direct, GRAV, TPI, TPIINV USE W3GDATMD, ONLY: NK, NTH, NSPEC, SIG, TH, DMIN, DTMAX, & DTMIN, FACTI1, FACTI2, FACSD, FACHFA, FACP, & XFC, XFLT, XREL, XFT, FXFM, FXPM, DDEN, & FTE, FTF, FHMAX, ECOS, ESIN, IICEDISP, & ICESCALES, IICESMOOTH USE W3WDATMD, ONLY: TIME USE W3ODATMD, ONLY: NDSE, NDST, IAPROC USE W3IDATMD, ONLY: INFLAGS2 USE W3DISPMD #ifdef W3_T USE CONSTANTS, ONLY: RADE #endif #ifdef W3_REF1 USE W3GDATMD, ONLY: IOBP, IOBPD, GTYPE, UNGTYPE, REFPARS #endif #ifdef W3_NNT USE W3ODATMD, ONLY: IAPROC, SCREEN, FNMPRE #endif #ifdef W3_FLD1 USE W3FLD1MD, ONLY: W3FLD1 USE W3GDATMD, ONLY: AALPHA #endif #ifdef W3_FLD2 USE W3FLD2MD, ONLY: W3FLD2 USE W3GDATMD, ONLY: AALPHA #endif #ifdef W3_FLX1 USE W3FLX1MD #endif #ifdef W3_FLX2 USE W3FLX2MD #endif #ifdef W3_FLX3 USE W3FLX3MD #endif #ifdef W3_FLX4 USE W3FLX4MD #endif #ifdef W3_FLX5 USE W3FLX5MD #endif #ifdef W3_LN1 USE W3SLN1MD #endif #ifdef W3_ST0 USE W3SRC0MD #endif #ifdef W3_ST1 USE W3SRC1MD #endif #ifdef W3_ST2 USE W3SRC2MD USE W3GDATMD, ONLY : ZWIND #endif #ifdef W3_ST3 USE W3SRC3MD USE W3GDATMD, ONLY : ZZWND, FFXFM, FFXPM #endif #ifdef W3_ST4 USE W3SRC4MD, ONLY : W3SPR4, W3SIN4, W3SDS4 USE W3GDATMD, ONLY : ZZWND, FFXFM, FFXPM, FFXFA #endif #ifdef W3_ST6 USE W3SRC6MD USE W3SWLDMD, ONLY : W3SWL6 USE W3GDATMD, ONLY : SWL6S6 #endif #ifdef W3_NL1 USE W3SNL1MD #endif #ifdef W3_NL2 USE W3SNL2MD #endif #ifdef W3_NL3 USE W3SNL3MD #endif #ifdef W3_NL4 USE W3SNL4MD #endif #ifdef W3_NL5 USE W3SNL5MD USE W3TIMEMD, ONLY: TICK21 #endif #ifdef W3_NLS USE W3SNLSMD #endif #ifdef W3_BT1 USE W3SBT1MD #endif #ifdef W3_BT4 USE W3SBT4MD #endif #ifdef W3_BT8 USE W3SBT8MD #endif #ifdef W3_BT9 USE W3SBT9MD #endif #ifdef W3_IC1 USE W3SIC1MD #endif #ifdef W3_IC2 USE W3SIC2MD #endif #ifdef W3_IC3 USE W3SIC3MD #endif #ifdef W3_IC4 USE W3SIC4MD #endif #ifdef W3_IC5 USE W3SIC5MD #endif #ifdef W3_IS1 USE W3SIS1MD #endif #ifdef W3_IS2 USE W3SIS2MD USE W3GDATMD, ONLY : IS2PARS #endif #ifdef W3_DB1 USE W3SDB1MD #endif #ifdef W3_TR1 USE W3STR1MD #endif #ifdef W3_BS1 USE W3SBS1MD #endif #ifdef W3_REF1 USE W3REF1MD #endif #ifdef W3_IG1 USE W3GDATMD, ONLY : IGPARS #endif #ifdef W3_S USE W3SERVMD, ONLY: STRACE #endif #ifdef W3_NNT USE W3SERVMD, ONLY: EXTCDE #endif #ifdef W3_UOST USE W3UOSTMD, ONLY: UOST_SRCTRMCOMPUTE #endif #ifdef W3_PDLIB USE PDLIB_W3PROFSMD, ONLY : B_JAC, ASPAR_JAC, ASPAR_DIAG_ALL USE yowNodepool, ONLY: PDLIB_I_DIAG, PDLIB_SI USE W3GDATMD, ONLY: B_JGS_LIMITER, FSSOURCE, optionCall USE W3GDATMD, ONLY: IOBP_LOC, IOBPD_LOC, B_JGS_LIMITER_FUNC USE W3WDATMD, ONLY: VA USE W3PARALL, ONLY: IMEM, LSLOC #endif !/ IMPLICIT NONE !/ !/ ------------------------------------------------------------------- / !/ Parameter list !/ INTEGER, INTENT(IN) :: srce_call, IT, ISEA, JSEA, IX, IY, IMOD REAL, intent(in) :: SPECOLD(NSPEC), CLATSL REAL, INTENT(OUT) :: VSIO(NSPEC), VDIO(NSPEC) LOGICAL, INTENT(OUT) :: SHAVEIO REAL, INTENT(IN) :: D_INP, U10ABS, & U10DIR, AS, CX, CY, DTG, D50,PSIC, & ICE, ICEH #ifdef W3_FLX5 REAL, INTENT(IN) :: TAUA, TAUADIR #endif INTEGER, INTENT(IN) :: REFLED(6) REAL, INTENT(IN) :: REFLEC(4), DELX, DELY, DELA, & TRNX, TRNY, BERG, ICEDMAX, DAIR REAL, INTENT(INOUT) :: WN1(NK), CG1(NK), & SPEC(NSPEC), ALPHA(NK), USTAR, & USTDIR, FPI, TAUOX, TAUOY, & TAUWX, TAUWY, PHIAW, PHIOC, PHICE, & CHARN, TWS, BEDFORM(3), PHIBBL, & TAUBBL(2), TAUICE(2), WHITECAP(4), & TAUWIX, TAUWIY, TAUWNX, TAUWNY, & ICEF, TAUOCX, TAUOCY, WNMEAN REAL, INTENT(OUT) :: DTDYN, FCUT REAL, INTENT(IN) :: COEF !/ !/ ------------------------------------------------------------------- / !/ Local parameters !/ INTEGER :: IK, ITH, IS, IS0, NSTEPS, NKH, NKH1, & IKS1, IS1, NSPECH, IDT, IERR, ISP REAL :: DTTOT, FHIGH, DT, AFILT, DAMAX, AFAC, & HDT, ZWND, FP, DEPTH, TAUSCX, TAUSCY, FHIGI ! Scaling factor for SIN, SDS, SNL REAL :: ICESCALELN, ICESCALEIN, ICESCALENL, ICESCALEDS REAL :: EMEAN, FMEAN, AMAX, CD, Z0, SCAT, & SMOOTH_ICEDISP REAL :: WN_R(NK), CG_ICE(NK), ALPHA_LIU(NK), ICECOEF2, R(NK) DOUBLE PRECISION :: ATT, ISO REAL :: EBAND, DIFF, EFINISH, HSTOT, PHINL, & FMEAN1, FMEANWS, & FACTOR, FACTOR2, DRAT, TAUWAX, TAUWAY, & MWXFINISH, MWYFINISH, A1BAND, B1BAND, & COSI(2) REAL :: SPECINIT(NSPEC), SPEC2(NSPEC), FRLOCAL, JAC2 REAL :: DAM (NSPEC), DAM2(NSPEC), WN2(NSPEC), & VSLN(NSPEC), & VSIN(NSPEC), VDIN(NSPEC), & VSNL(NSPEC), VDNL(NSPEC), & VSDS(NSPEC), VDDS(NSPEC), & VSBT(NSPEC), VDBT(NSPEC) REAL :: VS(NSPEC), VD(NSPEC), EB(NK) LOGICAL :: SHAVE LOGICAL :: LBREAK LOGICAL, SAVE :: FIRST = .TRUE. LOGICAL :: PrintDeltaSmDA REAL :: eInc1, eInc2, eVS, eVD, JAC REAL :: DeltaSRC(NSPEC) REAL :: FOUT(NK,NTH), SOUT(NK,NTH), DOUT(NK,NTH) REAL, SAVE :: TAUNUX, TAUNUY LOGICAL, SAVE :: FLTEST = .FALSE., FLAGNN = .TRUE. #ifdef W3_OMPG !$omp threadprivate( TAUNUX, TAUNUY) !$omp threadprivate( FLTEST, FLAGNN ) !$omp threadprivate( FIRST ) #endif !/ !/ ------------------------------------------------------------------- / !/ Local parameters dependent on compile switch !/ #ifdef W3_S INTEGER, SAVE :: IENT = 0 #endif #ifdef W3_NNT INTEGER, SAVE :: NDSD = 89, NDSD2 = 88, J REAL :: QCERR = 0. !/XNL2 and !/NNT #endif #ifdef W3_NL5 INTEGER :: QI5TSTART(2) REAL :: QR5KURT INTEGER, PARAMETER :: NL5_SELECT = 1 REAL, PARAMETER :: NL5_OFFSET = 0. ! explicit dyn. #endif #ifdef W3_SEED REAL :: UC, SLEV #endif #ifdef W3_MLIM REAL :: HM, EM #endif #ifdef W3_NNT REAL :: FACNN #endif #ifdef W3_T REAL :: DTRAW #endif #if defined(W3_IC1) || defined(W3_IC2) || defined(W3_IC3) || defined(W3_IC4) || defined(W3_IC5) REAL :: VSIC(NSPEC), VDIC(NSPEC) #endif #ifdef W3_DB1 REAL :: VSDB(NSPEC), VDDB(NSPEC) #endif #ifdef W3_TR1 REAL :: VSTR(NSPEC), VDTR(NSPEC) #endif #ifdef W3_BS1 REAL :: VSBS(NSPEC), VDBS(NSPEC) #endif #ifdef W3_REF1 REAL :: VREF(NSPEC) #endif #if defined(W3_IS1) || defined(W3_IS2) REAL :: VSIR(NSPEC), VDIR(NSPEC) #endif #ifdef W3_IS2 REAL :: VDIR2(NSPEC) DOUBLE PRECISION :: SCATSPEC(NTH) #endif #ifdef W3_UOST REAL :: VSUO(NSPEC), VDUO(NSPEC) #endif #ifdef W3_ST1 REAL :: FH1, FH2 #endif #ifdef W3_ST2 REAL :: FHTRAN, DFH, FACDIA, FACPAR #endif #ifdef W3_ST3 REAL :: FMEANS, FH1, FH2 #endif #ifdef W3_ST4 REAL :: FMEANS, FH1, FH2, FAGE, DLWMEAN REAL :: BRLAMBDA(NSPEC) #endif #if defined(W3_ST3) || defined(W3_ST4) LOGICAL :: LLWS(NSPEC) #endif #ifdef W3_ST6 REAL :: VSWL(NSPEC), VDWL(NSPEC) #endif #ifdef W3_PDLIB REAL :: PreVS, DVS, SIDT, FAKS, MAXDAC #endif #ifdef W3_NNT CHARACTER(LEN=17), SAVE :: FNAME = 'test_data_nnn.ww3' #endif ! !/ -- End of variable delclarations ! #ifdef W3_S CALL STRACE (IENT, 'W3SRCE') #endif #ifdef W3_T FLTEST = .TRUE. #endif ! IKS1 = 1 #ifdef W3_IG1 ! Does not integrate source terms for IG band if IGPARS(12) = 0. IF (NINT(IGPARS(12)).EQ.0) IKS1 = NINT(IGPARS(5)) #endif IS1=(IKS1-1)*NTH+1 !! Initialise source term arrays: VD = 0. VS = 0. VDIO = 0. VSIO = 0. VSBT = 0. VDBT = 0. #if defined(W3_LN0) || defined(W3_LN1) || defined(W3_SEED) VSLN = 0. #endif #if defined(W3_ST0) || defined(W3_ST3) || defined(W3_ST4) VSIN = 0. VDIN = 0. #endif #if defined(W3_NL0) || defined(W3_NL1) VSNL = 0. VDNL = 0. #endif #ifdef W3_TR1 VSTR = 0. VDTR = 0. #endif #if defined(W3_ST0) || defined(W3_ST4) VSDS = 0. VDDS = 0. #endif #ifdef W3_DB1 VSDB = 0. VDDB = 0. #endif #if defined(W3_IC1) || defined(W3_IC2) || defined(W3_IC3) || defined(W3_IC4) || defined(W3_IC5) VSIC = 0. VDIC = 0. #endif #ifdef W3_UOST VSUO = 0. VDUO = 0. #endif #if defined(W3_IS1) || defined(W3_IS2) VSIR = 0. VDIR = 0. #endif #ifdef W3_IS2 VDIR2 = 0. #endif #ifdef W3_ST6 VSWL = 0. VDWL = 0. #endif #if defined(W3_ST0) || defined(W3_ST1) || defined(W3_ST6) ZWND = 10. #endif #if defined(W3_ST2) ZWND = ZWIND #endif #if defined(W3_ST4) ZWND = ZZWND #endif ! ! 1. Preparations --------------------------------------------------- * ! DEPTH = MAX ( DMIN , D_INP ) DRAT = DAIR / DWAT ICESCALELN = MAX(0.,MIN(1.,1.-ICE*ICESCALES(1))) ICESCALEIN = MAX(0.,MIN(1.,1.-ICE*ICESCALES(2))) ICESCALENL = MAX(0.,MIN(1.,1.-ICE*ICESCALES(3))) ICESCALEDS = MAX(0.,MIN(1.,1.-ICE*ICESCALES(4))) #ifdef W3_T WRITE (NDST,9000) WRITE (NDST,9001) DEPTH, U10ABS, U10DIR*RADE #endif ! 1.a Set maximum change and wavenumber arrays. ! !XP = 0.15 !FACP = XP / PI * 0.62E-3 * TPI**4 / GRAV**2 ! DO IK=1, NK DAM(1+(IK-1)*NTH) = FACP / ( SIG(IK) * WN1(IK)**3 ) WN2(1+(IK-1)*NTH) = WN1(IK) END DO ! DO IK=1, NK IS0 = (IK-1)*NTH DO ITH=2, NTH DAM(ITH+IS0) = DAM(1+IS0) WN2(ITH+IS0) = WN2(1+IS0) END DO END DO ! ! 1.b Prepare dynamic time stepping ! DTDYN = 0. DTTOT = 0. NSTEPS = 0 PHIAW = 0. CHARN = 0. TWS = 0. PHINL = 0. PHIBBL = 0. TAUWIX = 0. TAUWIY = 0. TAUWNX = 0. TAUWNY = 0. TAUWAX = 0. TAUWAY = 0. TAUSCX = 0. TAUSCY = 0. TAUBBL = 0. TAUICE = 0. PHICE = 0. TAUOCX = 0. TAUOCY = 0. WNMEAN = 0. ! ! TIME is updated in W3WAVEMD prior to the call of W3SCRE, we should ! move 'TIME' one time step backward (QL) #ifdef W3_NL5 QI5TSTART = TIME CALL TICK21 (QI5TSTART, -1.0 * DTG) #endif ! #ifdef W3_DEBUGSRC IF (IX .eq. DEBUG_NODE) THEN WRITE(740+IAPROC,*) 'W3SRCE start sum(SPEC)=', sum(SPEC) WRITE(740+IAPROC,*) 'W3SRCE start sum(SPECOLD)=', sum(SPECOLD) WRITE(740+IAPROC,*) 'W3SRCE start sum(SPECINIT)=', sum(SPECINIT) WRITE(740+IAPROC,*) 'W3SRCE start sum(VSIO)=', sum(VSIO) WRITE(740+IAPROC,*) 'W3SRCE start sum(VDIO)=', sum(VDIO) WRITE(740+IAPROC,*) 'W3SRCE start USTAR=', USTAR END IF #endif #ifdef W3_ST4 DLWMEAN= 0. BRLAMBDA(:)=0. WHITECAP(:)=0. #endif ! ! 1.c Set mean parameters ! #ifdef W3_ST0 CALL W3SPR0 (SPEC, CG1, WN1, EMEAN, FMEAN, WNMEAN, AMAX) FP = 0.85 * FMEAN #endif #ifdef W3_ST1 CALL W3SPR1 (SPEC, CG1, WN1, EMEAN, FMEAN, WNMEAN, AMAX) FP = 0.85 * FMEAN #endif #ifdef W3_ST2 CALL W3SPR2 (SPEC, CG1, WN1, DEPTH, FPI, U10ABS, USTAR, & EMEAN, FMEAN, WNMEAN, AMAX, ALPHA, FP ) #endif #ifdef W3_ST3 TAUWX=0. TAUWY=0. IF ( IT .eq. 0 ) THEN LLWS(:) = .TRUE. USTAR=0. USTDIR=0. CALL W3SPR3 (SPEC, CG1, WN1, EMEAN, FMEAN, FMEANS, WNMEAN, & AMAX, U10ABS, U10DIR, USTAR, USTDIR, & TAUWX, TAUWY, CD, Z0, CHARN, LLWS, FMEANWS) ELSE CALL W3SPR3 (SPEC, CG1, WN1, EMEAN, FMEAN, FMEANS, WNMEAN, & AMAX, U10ABS, U10DIR, USTAR, USTDIR, & TAUWX, TAUWY, CD, Z0, CHARN, LLWS, FMEANWS) CALL W3SIN3 ( SPEC, CG1, WN2, U10ABS, USTAR, DRAT, AS, & U10DIR, Z0, CD, TAUWX, TAUWY, TAUWAX, TAUWAY, & ICE, VSIN, VDIN, LLWS, IX, IY ) END IF CALL W3SPR3 (SPEC, CG1, WN1, EMEAN, FMEAN, FMEANS, WNMEAN, & AMAX, U10ABS, U10DIR, USTAR, USTDIR, & TAUWX, TAUWY, CD, Z0, CHARN, LLWS, FMEANWS) TWS = 1./FMEANWS #endif #ifdef W3_ST4 TAUWX=0. TAUWY=0. IF ( IT .eq. 0 ) THEN LLWS(:) = .TRUE. USTAR=0. USTDIR=0. ELSE CALL W3SPR4 (SPEC, CG1, WN1, EMEAN, FMEAN, FMEAN1, WNMEAN, & AMAX, U10ABS, U10DIR, & #ifdef W3_FLX5 TAUA, TAUADIR, DAIR, & #endif USTAR, USTDIR, & TAUWX, TAUWY, CD, Z0, CHARN, LLWS, FMEANWS, DLWMEAN) #endif #if defined(W3_DEBUGSRC) && defined(W3_ST4) IF (IX == DEBUG_NODE) THEN WRITE(740+IAPROC,*) '1: out value USTAR=', USTAR, ' USTDIR=', USTDIR WRITE(740+IAPROC,*) '1: out value EMEAN=', EMEAN, ' FMEAN=', FMEAN WRITE(740+IAPROC,*) '1: out value FMEAN1=', FMEAN1, ' WNMEAN=', WNMEAN WRITE(740+IAPROC,*) '1: out value CD=', CD, ' Z0=', Z0 WRITE(740+IAPROC,*) '1: out value ALPHA=', CHARN, ' FMEANWS=', FMEANWS END IF #endif #ifdef W3_ST4 CALL W3SIN4 ( SPEC, CG1, WN2, U10ABS, USTAR, DRAT, AS, & U10DIR, Z0, CD, TAUWX, TAUWY, TAUWAX, TAUWAY, & VSIN, VDIN, LLWS, IX, IY, BRLAMBDA ) END IF #endif #if defined(W3_DEBUGSRC) && defined(W3_ST4) IF (IX == DEBUG_NODE) THEN WRITE(740+IAPROC,*) '1: U10DIR=', U10DIR, ' Z0=', Z0, ' CHARN=', CHARN WRITE(740+IAPROC,*) '1: USTAR=', USTAR, ' U10ABS=', U10ABS, ' AS=', AS WRITE(740+IAPROC,*) '1: DRAT=', DRAT WRITE(740+IAPROC,*) '1: TAUWX=', TAUWX, ' TAUWY=', TAUWY WRITE(740+IAPROC,*) '1: TAUWAX=', TAUWAX, ' TAUWAY=', TAUWAY WRITE(740+IAPROC,*) '1: min(CG1)=', minval(CG1), ' max(CG1)=', maxval(CG1) WRITE(740+IAPROC,*) '1: W3SIN4(min/max/sum)VSIN=', minval(VSIN), maxval(VSIN), sum(VSIN) WRITE(740+IAPROC,*) '1: W3SIN4(min/max/sum)VDIN=', minval(VDIN), maxval(VDIN), sum(VDIN) END IF #endif #ifdef W3_ST4 CALL W3SPR4 (SPEC, CG1, WN1, EMEAN, FMEAN, FMEAN1, WNMEAN, & AMAX, U10ABS, U10DIR, & #ifdef W3_FLX5 TAUA, TAUADIR, DAIR, & #endif USTAR, USTDIR, & TAUWX, TAUWY, CD, Z0, CHARN, LLWS, FMEANWS, DLWMEAN) TWS = 1./FMEANWS #endif #ifdef W3_ST6 CALL W3SPR6 (SPEC, CG1, WN1, EMEAN, FMEAN, WNMEAN, AMAX, FP) #endif ! ! 1.c2 Stores the initial data ! SPECINIT = SPEC ! ! 1.d Stresses ! #ifdef W3_FLX1 CALL W3FLX1 ( ZWND, U10ABS, U10DIR, USTAR, USTDIR, Z0, CD ) #endif #ifdef W3_FLX2 CALL W3FLX2 ( ZWND, DEPTH, FP, U10ABS, U10DIR, & USTAR, USTDIR, Z0, CD ) #endif #ifdef W3_FLX3 CALL W3FLX3 ( ZWND, DEPTH, FP, U10ABS, U10DIR, & USTAR, USTDIR, Z0, CD ) #endif #ifdef W3_FLX4 CALL W3FLX4 ( ZWND, U10ABS, U10DIR, USTAR, USTDIR, Z0, CD ) #endif #ifdef W3_FLX5 CALL W3FLX5 ( ZWND, U10ABS, U10DIR, TAUA, TAUADIR, DAIR, & USTAR, USTDIR, Z0, CD, CHARN ) #endif ! ! 1.e Prepare cut-off beyond which the tail is imposed with a power law ! #ifdef W3_ST0 FHIGH = SIG(NK) #endif #ifdef W3_ST1 FH1 = FXFM * FMEAN FH2 = FXPM / USTAR FHIGH = MAX ( FH1 , FH2 ) IF (FLTEST) WRITE (NDST,9004) FH1*TPIINV, FH2*TPIINV, FHIGH*TPIINV #endif #ifdef W3_ST2 FHIGH = XFC * FPI #endif #ifdef W3_ST3 FHIGH = MAX(FFXFM * MAX(FMEAN,FMEANWS),FFXPM / USTAR) #endif #ifdef W3_ST4 ! Introduces a Long & Resio (JGR2007) type dependance on wave age ! !/ST4 FAGE = FFXFA*TANH(0.3*U10ABS*FMEANWS*TPI/GRAV) FAGE = 0. FHIGH = MAX( (FFXFM + FAGE ) * MAX(FMEAN1,FMEANWS), FFXPM / USTAR) FHIGI = FFXFA * FMEAN1 #endif #ifdef W3_ST6 IF (FXFM .LE. 0) THEN FHIGH = SIG(NK) ELSE FHIGH = MAX (FXFM * FMEAN, FXPM / USTAR) ENDIF #endif ! ! 1.f Prepare output file for !/NNT option ! #ifdef W3_NNT IF ( IT .EQ. 0 ) THEN J = LEN_TRIM(FNMPRE) WRITE (FNAME(11:13),'(I3.3)') IAPROC OPEN (NDSD,FILE=FNMPRE(:J)//FNAME,form='UNFORMATTED', convert=file_endian, & ERR=800,IOSTAT=IERR) WRITE (NDSD,ERR=801,IOSTAT=IERR) NK, NTH WRITE (NDSD,ERR=801,IOSTAT=IERR) SIG(1:NK) * TPIINV OPEN (NDSD2,FILE=FNMPRE(:J)//'time.ww3', & FORM='FORMATTED',ERR=800,IOSTAT=IERR) END IF #endif ! ! ... Branch point dynamic integration - - - - - - - - - - - - - - - - ! DO ! NSTEPS = NSTEPS + 1 ! #ifdef W3_T WRITE (NDST,9020) NSTEPS, DTTOT #endif ! ! 2. Calculate source terms ----------------------------------------- * ! ! 2.a Input. ! #ifdef W3_LN1 CALL W3SLN1 ( WN1, FHIGH, USTAR, U10DIR , VSLN ) #endif ! #ifdef W3_ST1 CALL W3SIN1 ( SPEC, WN2, USTAR, U10DIR , VSIN, VDIN ) #endif #ifdef W3_ST2 CALL W3SIN2 ( SPEC, CG1, WN2, U10ABS, U10DIR, CD, Z0, & FPI, VSIN, VDIN ) #endif #ifdef W3_ST3 CALL W3SIN3 ( SPEC, CG1, WN2, U10ABS, USTAR, DRAT, AS, & U10DIR, Z0, CD, TAUWX, TAUWY, TAUWAX, TAUWAY, & ICE, VSIN, VDIN, LLWS, IX, IY ) #endif #ifdef W3_ST4 CALL W3SIN4 ( SPEC, CG1, WN2, U10ABS, USTAR, DRAT, AS, & U10DIR, Z0, CD, TAUWX, TAUWY, TAUWAX, TAUWAY, & VSIN, VDIN, LLWS, IX, IY, BRLAMBDA ) #endif #if defined(W3_DEBUGSRC) && defined(W3_ST4) IF (IX == DEBUG_NODE) THEN WRITE(740+IAPROC,*) '2 : W3SIN4(min/max/sum)VSIN=', minval(VSIN), maxval(VSIN), sum(VSIN) WRITE(740+IAPROC,*) '2 : W3SIN4(min/max/sum)VDIN=', minval(VDIN), maxval(VDIN), sum(VDIN) END IF #endif #ifdef W3_ST6 CALL W3SIN6 ( SPEC, CG1, WN2, U10ABS, USTAR, USTDIR, CD, DAIR, & TAUWX, TAUWY, TAUWAX, TAUWAY, VSIN, VDIN ) #endif ! ! 2.b Nonlinear interactions. ! #ifdef W3_NL1 CALL W3SNL1 ( SPEC, CG1, WNMEAN*DEPTH, VSNL, VDNL ) #endif #ifdef W3_NL2 CALL W3SNL2 ( SPEC, CG1, DEPTH, VSNL, VDNL ) #endif #ifdef W3_NL3 CALL W3SNL3 ( SPEC, CG1, WN1, DEPTH, VSNL, VDNL ) #endif #ifdef W3_NL4 CALL W3SNL4 ( SPEC, CG1, WN1, DEPTH, VSNL, VDNL ) #endif #ifdef W3_NL5 CALL W3SNL5 ( SPEC, CG1, WN1, FMEAN, QI5TSTART, & U10ABS, U10DIR, JSEA, VSNL, VDNL, QR5KURT) #endif ! #ifdef W3_PDLIB IF (.NOT. FSSOURCE .or. LSLOC) THEN #endif #ifdef W3_TR1 CALL W3STR1 ( SPEC, SPECOLD, CG1, WN1, DEPTH, IX, VSTR, VDTR ) #endif #ifdef W3_PDLIB ENDIF #endif ! ! 2.c Dissipation... except for ST4 ! 2.c1 as in source term package ! #ifdef W3_ST1 CALL W3SDS1 ( SPEC, WN2, EMEAN, FMEAN, WNMEAN, VSDS, VDDS ) #endif #ifdef W3_ST2 CALL W3SDS2 ( SPEC, CG1, WN1, FPI, USTAR, ALPHA,VSDS, VDDS ) #endif #ifdef W3_ST3 CALL W3SDS3 ( SPEC, WN1, CG1, EMEAN, FMEANS, WNMEAN, & USTAR, USTDIR, DEPTH, VSDS, VDDS, IX, IY ) #endif #ifdef W3_ST4 CALL W3SDS4 ( SPEC, WN1, CG1, USTAR, USTDIR, DEPTH, DAIR, VSDS, & VDDS, IX, IY, BRLAMBDA, WHITECAP, DLWMEAN ) #endif #if defined(W3_DEBUGSRC) && defined(W3_ST4) IF (IX == DEBUG_NODE) THEN WRITE(740+IAPROC,*) '2 : W3SDS4(min/max/sum)VSDS=', minval(VSDS), maxval(VSDS), sum(VSDS) WRITE(740+IAPROC,*) '2 : W3SDS4(min/max/sum)VDDS=', minval(VDDS), maxval(VDDS), sum(VDDS) END IF #endif #ifdef W3_ST6 CALL W3SDS6 ( SPEC, CG1, WN1, VSDS, VDDS ) #endif ! #ifdef W3_PDLIB IF (.NOT. FSSOURCE .or. LSLOC) THEN #endif #ifdef W3_DB1 CALL W3SDB1 ( IX, SPEC, DEPTH, EMEAN, FMEAN, WNMEAN, CG1, & LBREAK, VSDB, VDDB ) #endif #ifdef W3_PDLIB ENDIF #endif ! ! 2.c2 optional dissipation parameterisations ! #ifdef W3_ST6 IF (SWL6S6) THEN CALL W3SWL6 ( SPEC, CG1, WN1, VSWL, VDWL ) END IF #endif ! ! 2.d Bottom interactions. ! #ifdef W3_BT1 CALL W3SBT1 ( SPEC, CG1, WN1, DEPTH, VSBT, VDBT ) #endif #ifdef W3_BT4 CALL W3SBT4 ( SPEC, CG1, WN1, DEPTH, D50, PSIC, TAUBBL, & BEDFORM, VSBT, VDBT, IX, IY ) #endif #ifdef W3_BT8 CALL W3SBT8 ( SPEC, DEPTH, VSBT, VDBT, IX, IY ) #endif #ifdef W3_BT9 CALL W3SBT9 ( SPEC, DEPTH, VSBT, VDBT, IX, IY ) #endif ! #ifdef W3_BS1 CALL W3SBS1 ( SPEC, CG1, WN1, DEPTH, CX, CY, & TAUSCX, TAUSCY, VSBS, VDBS ) #endif ! ! 2.e Unresolved Obstacles Source Term ! #ifdef W3_UOST ! UNRESOLVED OBSTACLES CALL UOST_SRCTRMCOMPUTE(IX, IY, SPEC, CG1, DT, & U10ABS, U10DIR, VSUO, VDUO) #endif ! ! 2.g Dump training data if necessary ! #ifdef W3_NNT WRITE (SCREEN,8888) TIME, DTTOT, FLAGNN, QCERR WRITE (NDSD2,8888) TIME, DTTOT, FLAGNN, QCERR 8888 FORMAT (1X,I8.8,1X,I6.6,F8.1,L2,F8.2) WRITE (NDSD,ERR=801,IOSTAT=IERR) IX, IY, TIME, NSTEPS, & DTTOT, FLAGNN, DEPTH, U10ABS, U10DIR ! IF ( FLAGNN ) THEN DO IK=1, NK FACNN = TPI * SIG(IK) / CG1(IK) DO ITH=1, NTH IS = ITH + (IK-1)*NTH FOUT(IK,ITH) = SPEC(IS) * FACNN SOUT(IK,ITH) = VSNL(IS) * FACNN DOUT(IK,ITH) = VDNL(IS) END DO END DO WRITE (NDSD,ERR=801,IOSTAT=IERR) FOUT WRITE (NDSD,ERR=801,IOSTAT=IERR) SOUT WRITE (NDSD,ERR=801,IOSTAT=IERR) DOUT END IF #endif ! ! 3. Set frequency cut-off ------------------------------------------ * ! #ifdef W3_ST2 FHIGH = XFC * FPI IF ( FLTEST ) WRITE (NDST,9005) FHIGH*TPIINV #endif NKH = MIN ( NK , INT(FACTI2+FACTI1*LOG(MAX(1.E-7,FHIGH))) ) NKH1 = MIN ( NK , NKH+1 ) NSPECH = NKH1*NTH #ifdef W3_T WRITE (NDST,9021) NKH, NKH1, NSPECH #endif ! ! 4. Summation of source terms and diagonal term and time step ------ * ! DT = MIN ( DTG-DTTOT , DTMAX ) AFILT = MAX ( DAM(NSPEC) , XFLT*AMAX ) ! ! For input and dissipation calculate the fraction of the ice-free ! surface. In the presence of ice, the effective water surface ! is reduce to a fraction of the cell size free from ice, and so is ! input : ! SIN = (1-ICE)**ISCALEIN*SIN and SDS=(1-ICE)**ISCALEDS*SDS ------------------ * ! INFLAGS2(4) is true if ice concentration was ever read during ! this simulation IF ( INFLAGS2(4) ) THEN VSNL(1:NSPECH) = ICESCALENL * VSNL(1:NSPECH) VDNL(1:NSPECH) = ICESCALENL * VDNL(1:NSPECH) VSLN(1:NSPECH) = ICESCALELN * VSLN(1:NSPECH) VSIN(1:NSPECH) = ICESCALEIN * VSIN(1:NSPECH) VDIN(1:NSPECH) = ICESCALEIN * VDIN(1:NSPECH) VSDS(1:NSPECH) = ICESCALEDS * VSDS(1:NSPECH) VDDS(1:NSPECH) = ICESCALEDS * VDDS(1:NSPECH) END IF #ifdef W3_PDLIB IF (B_JGS_LIMITER_FUNC == 2) THEN DO IK=1, NK JAC = CG1(IK)/CLATSL JAC2 = 1./TPI/SIG(IK) FRLOCAL = SIG(IK)*TPIINV #ifdef W3_ST6 DAM2(1+(IK-1)*NTH) = 5E-7 * GRAV/FRLOCAL**4 * USTAR * FMEAN * DTMIN * JAC * JAC2 #else DAM2(1+(IK-1)*NTH) = 5E-7 * GRAV/FRLOCAL**4 * USTAR * MAX(FMEANWS,FMEAN) * DTMIN * JAC * JAC2 #endif !FROM WWM: 5E-7 * GRAV/FR(IS)**4 * USTAR * MAX(FMEANWS(IP),FMEAN(IP)) * DT4S/PI2/SPSIG(IS) END DO DO IK=1, NK IS0 = (IK-1)*NTH DO ITH=2, NTH DAM2(ITH+IS0) = DAM2(1+IS0) END DO END DO ENDIF #endif ! DO IS=IS1, NSPECH VS(IS) = VSLN(IS) + VSIN(IS) + VSNL(IS) & + VSDS(IS) + VSBT(IS) #ifdef W3_ST6 VS(IS) = VS(IS) + VSWL(IS) #endif #if defined(W3_TR1) && !defined(W3_PDLIB) VS(IS) = VS(IS) + VSTR(IS) #endif #ifdef W3_BS1 VS(IS) = VS(IS) + VSBS(IS) #endif #ifdef W3_UOST VS(IS) = VS(IS) + VSUO(IS) #endif VD(IS) = VDIN(IS) + VDNL(IS) & + VDDS(IS) + VDBT(IS) #ifdef W3_ST6 VD(IS) = VD(IS) + VDWL(IS) #endif #if defined(W3_TR1) && !defined(W3_PDLIB) VD(IS) = VD(IS) + VDTR(IS) #endif #ifdef W3_BS1 VD(IS) = VD(IS) + VDBS(IS) #endif #ifdef W3_UOST VD(IS) = VD(IS) + VDUO(IS) #endif DAMAX = MIN ( DAM(IS) , MAX ( XREL*SPECINIT(IS) , AFILT ) ) AFAC = 1. / MAX( 1.E-10 , ABS(VS(IS)/DAMAX) ) #ifdef W3_NL5 IF (NL5_SELECT .EQ. 1) THEN DT = MIN ( DT , AFAC / ( MAX ( 1.E-10, & 1. + NL5_OFFSET*AFAC*MIN(0.,VD(IS)) ) ) ) ELSE #endif DT = MIN ( DT , AFAC / ( MAX ( 1.E-10, & 1. + OFFSET*AFAC*MIN(0.,VD(IS)) ) ) ) #ifdef W3_NL5 ENDIF #endif END DO ! end of loop on IS ! DT = MAX ( 0.5, DT ) ! The hardcoded min. dt is a problem for certain cases e.g. laborotary scale problems. ! DTDYN = DTDYN + DT #ifdef W3_T DTRAW = DT #endif IDT = 1 + INT ( 0.99*(DTG-DTTOT)/DT ) ! number of iterations DT = (DTG-DTTOT)/REAL(IDT) ! actualy time step SHAVE = DT.LT.DTMIN .AND. DT.LT.DTG-DTTOT ! limiter check ... SHAVEIO = SHAVE DT = MAX ( DT , MIN (DTMIN,DTG-DTTOT) ) ! override dt with input time step or last time step if it is bigger ... anyway the limiter is on! ! #ifdef W3_NL5 DT = INT(DT) * 1.0 #endif IF (srce_call .eq. srce_imp_post) DT = DTG ! for implicit part #ifdef W3_NL5 IF (NL5_SELECT .EQ. 1) THEN HDT = NL5_OFFSET * DT ELSE #endif HDT = OFFSET * DT #ifdef W3_NL5 ENDIF #endif DTTOT = DTTOT + DT #ifdef W3_DEBUGSRC IF (IX == DEBUG_NODE) WRITE(*,'(A20,2I10,5F20.10,L20)') 'TIMINGS 2', IDT, NSTEPS, DT, DTMIN, DTDYN, HDT, DTTOT, SHAVE IF (IX == DEBUG_NODE) THEN WRITE(740+IAPROC,*) '1: min/max/sum(VS)=', minval(VS), maxval(VS), sum(VS) WRITE(740+IAPROC,*) '1: min/max/sum(VD)=', minval(VD), maxval(VD), sum(VD) WRITE(740+IAPROC,*) 'min/max/sum(VSIN)=', minval(VSIN), maxval(VSIN), sum(VSIN) WRITE(740+IAPROC,*) 'min/max/sum(VDIN)=', minval(VDIN), maxval(VDIN), sum(VDIN) WRITE(740+IAPROC,*) 'min/max/sum(VSLN)=', minval(VSLN), maxval(VSLN), sum(VSLN) WRITE(740+IAPROC,*) 'min/max/sum(VSNL)=', minval(VSNL), maxval(VSNL), sum(VSNL) WRITE(740+IAPROC,*) 'min/max/sum(VDNL)=', minval(VDNL), maxval(VDNL), sum(VDNL) WRITE(740+IAPROC,*) 'min/max/sum(VSDS)=', minval(VSDS), maxval(VSDS), sum(VSDS) WRITE(740+IAPROC,*) 'min/max/sum(VDDS)=', minval(VDDS), maxval(VDDS), sum(VDDS) #ifdef W3_ST6 WRITE(740+IAPROC,*) 'min/max/sum(VSWL)=', minval(VSWL), maxval(VSWL), sum(VSWL) WRITE(740+IAPROC,*) 'min/max/sum(VDWL)=', minval(VDWL), maxval(VDWL), sum(VDWL) #endif #ifdef W3_DB1 WRITE(740+IAPROC,*) 'min/max/sum(VSDB)=', minval(VSDB), maxval(VSDB), sum(VSDB) WRITE(740+IAPROC,*) 'min/max/sum(VDDB)=', minval(VDDB), maxval(VDDB), sum(VDDB) #endif #ifdef W3_TR1 WRITE(740+IAPROC,*) 'min/max/sum(VSTR)=', minval(VSTR), maxval(VSTR), sum(VSTR) WRITE(740+IAPROC,*) 'min/max/sum(VDTR)=', minval(VDTR), maxval(VDTR), sum(VDTR) #endif #ifdef W3_BS1 WRITE(740+IAPROC,*) 'min/max/sum(VSBS)=', minval(VSBS), maxval(VSBS), sum(VSBS) WRITE(740+IAPROC,*) 'min/max/sum(VDBS)=', minval(VDBS), maxval(VDBS), sum(VDBS) #endif WRITE(740+IAPROC,*) 'min/max/sum(VSBT)=', minval(VSBT), maxval(VSBT), sum(VSBT) WRITE(740+IAPROC,*) 'min/max/sum(VDBT)=', minval(VDBT), maxval(VDBT), sum(VDBT) END IF #endif #ifdef W3_PDLIB IF (srce_call .eq. srce_imp_pre) THEN IF (LSLOC) THEN IF (IMEM == 1) THEN SIDT = PDLIB_SI(JSEA) * DTG DO IK = 1, NK JAC = CLATSL/CG1(IK) DO ITH = 1, NTH ISP = ITH + (IK-1)*NTH VD(ISP) = MIN(0., VD(ISP)) IF (B_JGS_LIMITER_FUNC == 2) THEN MAXDAC = MAX(DAM(ISP),DAM2(ISP)) ELSE MAXDAC = DAM(ISP) ENDIF FAKS = DTG / MAX ( 1. , (1.-DTG*VD(ISP))) DVS = VS(ISP) * FAKS IF (.NOT. B_JGS_LIMITER) THEN DVS = SIGN(MIN(MAXDAC,ABS(DVS)),DVS) ENDIF PreVS = DVS / FAKS eVS = PreVS / CG1(IK) * CLATSL eVD = MIN(0.,VD(ISP)) B_JAC(ISP,JSEA) = B_JAC(ISP,JSEA) + SIDT * (eVS - eVD*SPEC(ISP)*JAC) ASPAR_JAC(ISP,PDLIB_I_DIAG(JSEA)) = ASPAR_JAC(ISP,PDLIB_I_DIAG(JSEA)) - SIDT * eVD #ifdef W3_DB1 eVS = VSDB(ISP) * JAC eVD = MIN(0.,VDDB(ISP)) IF (eVS .gt. 0.) THEN evS = 2*evS evD = -evD ELSE evS = -evS evD = 2*evD ENDIF #endif B_JAC(ISP,JSEA) = B_JAC(ISP,JSEA) + SIDT * eVS ASPAR_JAC(ISP,PDLIB_I_DIAG(JSEA)) = ASPAR_JAC(ISP,PDLIB_I_DIAG(JSEA)) - SIDT * eVD #ifdef W3_TR1 eVS = VSTR(ISP) * JAC eVD = VDTR(ISP) IF (eVS .gt. 0.) THEN evS = 2*evS evD = -evD ELSE evS = -evS evD = 2*evD ENDIF #endif B_JAC(ISP,JSEA) = B_JAC(ISP,JSEA) + SIDT * eVS ASPAR_JAC(ISP,PDLIB_I_DIAG(JSEA)) = ASPAR_JAC(ISP,PDLIB_I_DIAG(JSEA)) - SIDT * eVD END DO END DO ELSEIF (IMEM == 2) THEN SIDT = PDLIB_SI(JSEA) * DTG DO IK=1,NK JAC = CLATSL/CG1(IK) DO ITH=1,NTH ISP=ITH + (IK-1)*NTH VD(ISP) = MIN(0., VD(ISP)) IF (B_JGS_LIMITER_FUNC == 2) THEN MAXDAC = MAX(DAM(ISP),DAM2(ISP)) ELSE MAXDAC = DAM(ISP) ENDIF FAKS = DTG / MAX ( 1. , (1.-DTG*VD(ISP))) DVS = VS(ISP) * FAKS IF (.NOT. B_JGS_LIMITER) THEN DVS = SIGN(MIN(MAXDAC,ABS(DVS)),DVS) ENDIF PreVS = DVS / FAKS eVS = PreVS / CG1(IK) * CLATSL eVD = VD(ISP) #ifdef W3_DB1 eVS = eVS + DBLE(VSDB(ISP)) * JAC eVD = evD + MIN(0.,DBLE(VDDB(ISP))) B_JAC(ISP,JSEA) = B_JAC(ISP,JSEA) + SIDT * (eVS - eVD*VA(ISP,JSEA)) ASPAR_DIAG_ALL(ISP,JSEA) = ASPAR_DIAG_ALL(ISP,JSEA) - SIDT * eVD #endif END DO END DO ENDIF ENDIF PrintDeltaSmDA=.FALSE. IF (PrintDeltaSmDA .eqv. .TRUE.) THEN DO IS=1,NSPEC DeltaSRC(IS) = VSIN(IS) - SPEC(IS)*VDIN(IS) END DO WRITE(740+IAPROC,*) 'min/max/sum(VSIN)=', minval(VSIN), maxval(VSIN), sum(VSIN) WRITE(740+IAPROC,*) 'min/max/sum(DeltaIN)=', minval(DeltaSRC), maxval(DeltaSRC), sum(DeltaSRC) ! DO IS=1,NSPEC DeltaSRC(IS) = VSNL(IS) - SPEC(IS)*VDNL(IS) END DO WRITE(740+IAPROC,*) 'min/max/sum(VSNL)=', minval(VSNL), maxval(VSNL), sum(VSNL) WRITE(740+IAPROC,*) 'min/max/sum(DeltaNL)=', minval(DeltaSRC), maxval(DeltaSRC), sum(DeltaSRC) ! DO IS=1,NSPEC DeltaSRC(IS) = VSDS(IS) - SPEC(IS)*VDDS(IS) END DO WRITE(740+IAPROC,*) 'min/max/sum(VSDS)=', minval(VSDS), maxval(VSDS), sum(VSDS) WRITE(740+IAPROC,*) 'min/max/sum(DeltaDS)=', minval(DeltaSRC), maxval(DeltaSRC), sum(DeltaSRC) ! ! DO IS=1,NSPEC ! DeltaSRC(IS) = VSIC(IS) - SPEC(IS)*VDIC(IS) ! END DO WRITE(740+IAPROC,*) 'min/max/sum(DeltaDS)=', minval(DeltaSRC), maxval(DeltaSRC), sum(DeltaSRC) END IF IF (.not. LSLOC) THEN IF (optionCall .eq. 1) THEN CALL SIGN_VSD_PATANKAR_WW3(SPEC,VS,VD) ELSE IF (optionCall .eq. 2) THEN CALL SIGN_VSD_SEMI_IMPLICIT_WW3(SPEC,VS,VD) ELSE IF (optionCall .eq. 3) THEN CALL SIGN_VSD_SEMI_IMPLICIT_WW3(SPEC,VS,VD) ENDIF VSIO = VS VDIO = VD ENDIF #ifdef W3_DEBUGSRC IF (IX == DEBUG_NODE) THEN WRITE(740+IAPROC,*) ' srce_imp_pre : SHAVE = ', SHAVE WRITE(740+IAPROC,*) ' srce_imp_pre : DT=', DT, ' HDT=', HDT, 'DTG=', DTG WRITE(740+IAPROC,*) ' srce_imp_pre : sum(SPEC)=', sum(SPEC) WRITE(740+IAPROC,*) ' srce_imp_pre : sum(VSTOT)=', sum(VS) WRITE(740+IAPROC,*) ' srce_imp_pre : sum(VDTOT)=', sum(MIN(0. , VD)) END IF IF (IX == DEBUG_NODE) WRITE(44,'(1EN15.4)') SUM(VSIN) IF (IX == DEBUG_NODE) WRITE(44,'(1EN15.4)') SUM(VDIN) IF (IX == DEBUG_NODE) WRITE(44,'(1EN15.4)') SUM(VSDS) IF (IX == DEBUG_NODE) WRITE(44,'(1EN15.4)') SUM(VDDS) IF (IX == DEBUG_NODE) WRITE(44,'(1EN15.4)') SUM(VSNL) IF (IX == DEBUG_NODE) WRITE(44,'(1EN15.4)') SUM(VDNL) IF (IX == DEBUG_NODE) WRITE(44,'(1EN15.4)') SUM(VSLN) IF (IX == DEBUG_NODE) WRITE(44,'(1EN15.4)') SUM(VSBT) IF (IX == DEBUG_NODE) WRITE(44,'(1EN15.4)') SUM(VS) IF (IX == DEBUG_NODE) WRITE(44,'(1EN15.4)') SUM(VD) #endif RETURN ! return everything is done for the implicit ... END IF ! srce_imp_pre ! --end W3_PDLIB #endif ! #ifdef W3_T WRITE (NDST,9040) DTRAW, DT, SHAVE #endif ! ! 5. Increment spectrum --------------------------------------------- * ! IF (srce_call .eq. srce_direct) THEN IF ( SHAVE ) THEN DO IS=IS1, NSPECH eInc1 = VS(IS) * DT / MAX ( 1. , (1.-HDT*VD(IS))) eInc2 = SIGN ( MIN (DAM(IS),ABS(eInc1)) , eInc1 ) SPEC(IS) = MAX ( 0. , SPEC(IS)+eInc2 ) END DO ELSE ! DO IS=IS1, NSPECH eInc1 = VS(IS) * DT / MAX ( 1. , (1.-HDT*VD(IS))) SPEC(IS) = MAX ( 0. , SPEC(IS)+eInc1 ) END DO END IF ! #ifdef W3_DB1 DO IS=IS1, NSPECH eInc1 = VSDB(IS) * DT / MAX ( 1. , (1.-HDT*VDDB(IS))) SPEC(IS) = MAX ( 0. , SPEC(IS)+eInc1 ) END DO #endif #ifdef W3_TR1 DO IS=IS1, NSPECH eInc1 = VDTR(IS) * DT / MAX ( 1. , (1.-HDT*VDTR(IS))) SPEC(IS) = MAX ( 0. , SPEC(IS)+eInc1 ) END DO #endif #ifdef W3_DEBUGSRC IF (IX == DEBUG_NODE) WRITE(44,'(1EN15.4)') SUM(VSIN) IF (IX == DEBUG_NODE) WRITE(44,'(1EN15.4)') SUM(VDIN) IF (IX == DEBUG_NODE) WRITE(44,'(1EN15.4)') SUM(VSDS) IF (IX == DEBUG_NODE) WRITE(44,'(1EN15.4)') SUM(VDDS) IF (IX == DEBUG_NODE) WRITE(44,'(1EN15.4)') SUM(VSNL) IF (IX == DEBUG_NODE) WRITE(44,'(1EN15.4)') SUM(VDNL) IF (IX == DEBUG_NODE) WRITE(44,'(1EN15.4)') SUM(VSLN) IF (IX == DEBUG_NODE) WRITE(44,'(1EN15.4)') SUM(VSBT) IF (IX == DEBUG_NODE) WRITE(44,'(1EN15.4)') SUM(VS) IF (IX == DEBUG_NODE) WRITE(44,'(1EN15.4)') SUM(VD) IF (IX == DEBUG_NODE) THEN WRITE(740+IAPROC,*) ' srce_direct : SHAVE = ', SHAVE WRITE(740+IAPROC,*) ' srce_direct : DT=', DT, ' HDT=', HDT, 'DTG=', DTG WRITE(740+IAPROC,*) ' srce_direct : sum(SPEC)=', sum(SPEC) WRITE(740+IAPROC,*) ' srce_direct : sum(VSTOT)=', sum(VS) WRITE(740+IAPROC,*) ' srce_direct : sum(VDTOT)=', sum(MIN(0. , VD)) END IF #endif END IF ! srce_call .eq. srce_direct ! ! 5.b Computes ! atmos->wave flux PHIAW-------------------------------- * ! wave ->BBL flux PHIBBL------------------------------- * ! wave ->ice flux PHICE ------------------------------- * ! WHITECAP(3)=0. HSTOT=0. DO IK=IKS1, NK FACTOR = DDEN(IK)/CG1(IK) !Jacobian to get energy in band FACTOR2= FACTOR*GRAV*WN1(IK)/SIG(IK) ! coefficient to get momentum ! Wave direction is "direction to" ! therefore there is a PLUS sign for the stress DO ITH=1, NTH IS = (IK-1)*NTH + ITH COSI(1)=ECOS(IS) COSI(2)=ESIN(IS) PHIAW = PHIAW + (VSIN(IS))* DT * FACTOR & / MAX ( 1. , (1.-HDT*VDIN(IS))) ! semi-implict integration scheme PHIBBL= PHIBBL- (VSBT(IS))* DT * FACTOR & / MAX ( 1. , (1.-HDT*VDBT(IS))) ! semi-implict integration scheme PHINL = PHINL + VSNL(IS)* DT * FACTOR & / MAX ( 1. , (1.-HDT*VDNL(IS))) ! semi-implict integration scheme IF (VSIN(IS).GT.0.) WHITECAP(3) = WHITECAP(3) + SPEC(IS) * FACTOR HSTOT = HSTOT + SPEC(IS) * FACTOR END DO END DO WHITECAP(3) = 4. * SQRT(WHITECAP(3)) HSTOT =4.*SQRT(HSTOT) TAUWIX = TAUWIX + TAUWX * DRAT * DT TAUWIY = TAUWIY + TAUWY * DRAT * DT TAUWNX = TAUWNX + TAUWAX * DRAT * DT TAUWNY = TAUWNY + TAUWAY * DRAT * DT ! MISSING: TAIL TO BE ADDED ? ! #ifdef W3_NLS CALL W3SNLS ( SPEC, CG1, WN1, DEPTH, U10ABS, DT, AA=SPEC ) #endif ! ! 6. Add tail ------------------------------------------------------- * ! a Mean parameters ! ! #ifdef W3_ST0 CALL W3SPR0 (SPEC, CG1, WN1, EMEAN, FMEAN, WNMEAN, AMAX) #endif #ifdef W3_ST1 CALL W3SPR1 (SPEC, CG1, WN1, EMEAN, FMEAN, WNMEAN, AMAX) #endif #ifdef W3_ST2 CALL W3SPR2 (SPEC, CG1, WN1, DEPTH, FPI, U10ABS, USTAR, & EMEAN, FMEAN, WNMEAN, AMAX, ALPHA, FP ) #endif #ifdef W3_ST3 CALL W3SPR3 (SPEC, CG1, WN1, EMEAN, FMEAN, FMEANS, & WNMEAN, AMAX, U10ABS, U10DIR, USTAR, USTDIR, & TAUWX, TAUWY, CD, Z0, CHARN, LLWS, FMEANWS) #endif #ifdef W3_ST4 CALL W3SPR4 (SPEC, CG1, WN1, EMEAN, FMEAN, FMEAN1, WNMEAN,& AMAX, U10ABS, U10DIR, & #ifdef W3_FLX5 TAUA, TAUADIR, DAIR, & #endif USTAR, USTDIR, & TAUWX, TAUWY, CD, Z0, CHARN, LLWS, FMEANWS, DLWMEAN) #endif #ifdef W3_ST6 CALL W3SPR6 (SPEC, CG1, WN1, EMEAN, FMEAN, WNMEAN, AMAX, FP) #endif ! #ifdef W3_FLX2 CALL W3FLX2 ( ZWND, DEPTH, FP, U10ABS, U10DIR, & USTAR, USTDIR, Z0, CD ) #endif #ifdef W3_FLX3 CALL W3FLX3 ( ZWND, DEPTH, FP, U10ABS, U10DIR, & USTAR, USTDIR, Z0, CD ) #endif ! #ifdef W3_ST1 FH1 = FXFM * FMEAN FHIGH = MIN ( SIG(NK) , MAX ( FH1 , FH2 ) ) NKH = MAX ( 2 , MIN ( NKH1 , & INT ( FACTI2 + FACTI1*LOG(MAX(1.E-7,FHIGH)) ) ) ) ! IF ( FLTEST ) WRITE (NDST,9060) & FH1*TPIINV, FH2*TPIINV, FHIGH*TPIINV, NKH #endif ! #ifdef W3_ST2 FHTRAN = XFT*FPI FHIGH = XFC*FPI DFH = FHIGH - FHTRAN NKH = MAX ( 1 , & INT ( FACTI2 + FACTI1*LOG(MAX(1.E-7,FHTRAN)) ) ) ! IF ( FLTEST ) WRITE (NDST,9061) FHTRAN, FHIGH, NKH #endif ! #ifdef W3_ST3 FH1 = FFXFM * FMEAN FH2 = FFXPM / USTAR FHIGH = MIN ( SIG(NK) , MAX ( FH1 , FH2 ) ) NKH = MAX ( 2 , MIN ( NKH1 , & INT ( FACTI2 + FACTI1*LOG(MAX(1.E-7,FHIGH)) ) ) ) ! IF ( FLTEST ) WRITE (NDST,9062) & FH1*TPIINV, FH2*TPIINV, FHIGH*TPIINV, NKH #endif ! #ifdef W3_ST4 ! Introduces a Long & Resio (JGR2007) type dependance on wave age FAGE = FFXFA*TANH(0.3*U10ABS*FMEANWS*TPI/GRAV) FH1 = (FFXFM+FAGE) * FMEAN1 FH2 = FFXPM / USTAR FHIGH = MIN ( SIG(NK) , MAX ( FH1 , FH2 ) ) NKH = MAX ( 2 , MIN ( NKH1 , & INT ( FACTI2 + FACTI1*LOG(MAX(1.E-7,FHIGH)) ) ) ) #endif ! #ifdef W3_ST6 IF (FXFM .LE. 0) THEN FHIGH = SIG(NK) ELSE FHIGH = MIN ( SIG(NK), MAX(FXFM * FMEAN, FXPM / USTAR) ) ENDIF NKH = MAX ( 2 , MIN ( NKH1 , & INT ( FACTI2 + FACTI1*LOG(MAX(1.E-7,FHIGH)) ) ) ) ! IF ( FLTEST ) WRITE (NDST,9063) FHIGH*TPIINV, NKH #endif ! ! 6.b Limiter for shallow water or Miche style criterion ! Last time step ONLY ! ! uses true depth (D_INP) instead of limited depth ! #ifdef W3_MLIM IF ( DTTOT .GE. 0.9999*DTG ) THEN HM = FHMAX *TANH(WNMEAN*MAX(0.,D_INP)) / MAX(1.E-4,WNMEAN ) EM = HM * HM / 16. IF ( EMEAN.GT.EM .AND. EMEAN.GT.1.E-30 ) THEN SPEC = SPEC / EMEAN * EM EMEAN = EM END IF END IF #endif ! ! 6.c Seeding of spectrum ! alpha = 0.005 , 0.5 in eq., 0.25 for directional distribution ! #ifdef W3_SEED DO IK=MIN(NK,NKH), NK UC = FACSD * GRAV / SIG(IK) SLEV = MIN ( 1. , MAX ( 0. , U10ABS/UC-1. ) ) * & 6.25E-4 / WN1(IK)**3 / SIG(IK) IF (INFLAGS2(4)) SLEV=SLEV*(1-ICE) DO ITH=1, NTH SPEC(ITH+(IK-1)*NTH) = MAX ( SPEC(ITH+(IK-1)*NTH) , & SLEV * MAX ( 0. , COS(U10DIR-TH(ITH)) )**2 ) END DO END DO #endif ! ! 6.d Add tail ! DO IK=NKH+1, NK #ifdef W3_ST2 FACDIA = MAX ( 0. , MIN ( 1., (SIG(IK)-FHTRAN)/DFH) ) FACPAR = MAX ( 0. , 1.-FACDIA ) #endif DO ITH=1, NTH SPEC(ITH+(IK-1)*NTH) = SPEC(ITH+(IK-2)*NTH) * FACHFA & #ifdef W3_ST2 * FACDIA + FACPAR * SPEC(ITH+(IK-1)*NTH) & #endif + 0. END DO END DO ! ! 6.e Update wave-supported stress----------------------------------- * ! #ifdef W3_ST3 CALL W3SIN3 ( SPEC, CG1, WN2, U10ABS, USTAR, DRAT, AS, & U10DIR, Z0, CD, TAUWX, TAUWY, TAUWAX, TAUWAY, & ICE, VSIN, VDIN, LLWS, IX, IY ) #endif #ifdef W3_ST4 CALL W3SIN4 ( SPEC, CG1, WN2, U10ABS, USTAR, DRAT, AS, & U10DIR, Z0, CD, TAUWX, TAUWY, TAUWAX, TAUWAY, & VSIN, VDIN, LLWS, IX, IY, BRLAMBDA ) #endif ! ! 7. Check if integration complete ---------------------------------- * ! ! Update QI5TSTART (Q. Liu) #ifdef W3_NL5 CALL TICK21(QI5TSTART, DT) #endif IF (srce_call .eq. srce_imp_post) THEN EXIT ENDIF IF ( DTTOT .GE. 0.9999*DTG ) THEN ! IF (IX == DEBUG_NODE) WRITE(*,*) 'DTTOT, DTG', DTTOT, DTG EXIT ENDIF END DO ! INTEGRATION LOOP #ifdef W3_DEBUGSRC IF (IX .eq. DEBUG_NODE) THEN WRITE(740+IAPROC,*) 'NSTEPS=', NSTEPS WRITE(740+IAPROC,*) '1 : sum(SPEC)=', sum(SPEC) END IF WRITE(740+IAPROC,*) 'DT=', DT, 'DTG=', DTG #endif ! ! ... End point dynamic integration - - - - - - - - - - - - - - - - - - ! ! 8. Save integration data ------------------------------------------ * ! DTDYN = DTDYN / REAL(MAX(1,NSTEPS)) FCUT = FHIGH * TPIINV ! GOTO 888 ! ! Error escape locations ! #ifdef W3_NNT 800 CONTINUE WRITE (NDSE,8000) FNAME, IERR CALL EXTCDE (1) ! 801 CONTINUE WRITE (NDSE,8001) IERR CALL EXTCDE (2) #endif ! 888 CONTINUE ! ! 9.a Computes PHIOC------------------------------------------ * ! The wave to ocean flux is the difference between initial energy ! and final energy, plus wind input plus the SNL flux to high freq., ! minus the energy lost to the bottom boundary layer (BBL) ! #ifdef W3_DEBUGSRC IF (IX .eq. DEBUG_NODE) THEN WRITE(740+IAPROC,*) '2 : sum(SPEC)=', sum(SPEC) END IF #endif EFINISH = 0. MWXFINISH = 0. MWYFINISH = 0. DO IK=1, NK EBAND = 0. A1BAND = 0. B1BAND = 0. DO ITH=1, NTH DIFF = SPECINIT(ITH+(IK-1)*NTH)-SPEC(ITH+(IK-1)*NTH) EBAND = EBAND + DIFF A1BAND = A1BAND + DIFF*ECOS(ITH) B1BAND = B1BAND + DIFF*ESIN(ITH) END DO EFINISH = EFINISH + EBAND * DDEN(IK) / CG1(IK) MWXFINISH = MWXFINISH + A1BAND * DDEN(IK) / CG1(IK) & * WN1(IK)/SIG(IK) MWYFINISH = MWYFINISH + B1BAND * DDEN(IK) / CG1(IK) & * WN1(IK)/SIG(IK) END DO ! ! Transformation in momentum flux in m^2 / s^2 ! TAUOX=(GRAV*MWXFINISH+TAUWIX-TAUBBL(1))/DTG TAUOY=(GRAV*MWYFINISH+TAUWIY-TAUBBL(2))/DTG TAUWIX=TAUWIX/DTG TAUWIY=TAUWIY/DTG TAUWNX=TAUWNX/DTG TAUWNY=TAUWNY/DTG TAUBBL(:)=TAUBBL(:)/DTG TAUOCX=DAIR*COEF*COEF*USTAR*USTAR*COS(USTDIR) + DWAT*(TAUOX-TAUWIX) TAUOCY=DAIR*COEF*COEF*USTAR*USTAR*SIN(USTDIR) + DWAT*(TAUOY-TAUWIY) ! ! Transformation in wave energy flux in W/m^2=kg / s^3 ! PHIOC =DWAT*GRAV*(EFINISH+PHIAW-PHIBBL)/DTG PHIAW =DWAT*GRAV*PHIAW /DTG PHINL =DWAT*GRAV*PHINL /DTG PHIBBL=DWAT*GRAV*PHIBBL/DTG ! ! 10.1 Adds ice scattering and dissipation: implicit integration---------------- * ! INFLAGS2(4) is true if ice concentration was ever read during ! this simulation ! #ifdef W3_DEBUGSRC IF (IX .eq. DEBUG_NODE) THEN WRITE(740+IAPROC,*) '3 : sum(SPEC)=', sum(SPEC) END IF #endif IF ( INFLAGS2(4).AND.ICE.GT.0 ) THEN IF (IICEDISP) THEN ICECOEF2 = 1E-6 CALL LIU_FORWARD_DISPERSION (ICEH,ICECOEF2,DEPTH, & SIG,WN_R,CG_ICE,ALPHA_LIU) ! IF (IICESMOOTH) THEN #ifdef W3_IS2 DO IK=1,NK SMOOTH_ICEDISP=0. IF (IS2PARS(14)*(TPI/WN_R(IK)).LT.ICEF) THEN ! IF ICE IS NOT TOO MUCH BROKEN SMOOTH_ICEDISP=TANH((ICEF-IS2PARS(14)*(TPI/WN_R(IK)))/(ICEF*IS2PARS(13))) END IF WN_R(IK)=WN1(IK)*(1-SMOOTH_ICEDISP)+WN_R(IK)*(SMOOTH_ICEDISP) END DO #endif END IF ELSE WN_R=WN1 CG_ICE=CG1 END IF ! R(:)=1 ! In case IC2 is defined but not IS2 ! #ifdef W3_IC1 CALL W3SIC1 ( SPEC,DEPTH, CG1, IX, IY, VSIC, VDIC ) #endif #ifdef W3_IS2 CALL W3SIS2 ( SPEC, DEPTH, ICE, ICEH, ICEF, ICEDMAX, IX, IY, & VSIR, VDIR, VDIR2, WN1, CG1, WN_R, CG_ICE, R ) #endif #ifdef W3_IC2 CALL W3SIC2 ( SPEC, DEPTH, ICEH, ICEF, CG1, WN1,& IX, IY, VSIC, VDIC, WN_R, CG_ICE, ALPHA_LIU, R) #endif #ifdef W3_IC3 CALL W3SIC3 ( SPEC,DEPTH, CG1, WN1, IX, IY, VSIC, VDIC ) #endif #ifdef W3_IC4 CALL W3SIC4 ( SPEC,DEPTH, CG1, IX, IY, VSIC, VDIC ) #endif #ifdef W3_IC5 CALL W3SIC5 ( SPEC,DEPTH, CG1, WN1, IX, IY, VSIC, VDIC ) #endif ! #ifdef W3_IS1 CALL W3SIS1 ( SPEC, ICE, VSIR ) #endif SPEC2 = SPEC ! TAUICE(:) = 0. PHICE = 0. DO IK=1,NK IS = 1+(IK-1)*NTH ! ! First part of ice term integration: dissipation part ! ATT=1. #ifdef W3_IC1 ATT=EXP(ICE*VDIC(IS)*DTG) #endif #ifdef W3_IC2 ATT=EXP(ICE*VDIC(IS)*DTG) #endif #ifdef W3_IC3 ATT=EXP(ICE*VDIC(IS)*DTG) #endif #ifdef W3_IC4 ATT=EXP(ICE*VDIC(IS)*DTG) #endif #ifdef W3_IC5 ATT=EXP(ICE*VDIC(IS)*DTG) #endif #ifdef W3_IS1 ATT=ATT*EXP(ICE*VDIR(IS)*DTG) #endif #ifdef W3_IS2 ATT=ATT*EXP(ICE*VDIR2(IS)*DTG) IF (IS2PARS(2).EQ.0) THEN ! Reminder : IS2PARS(2) = IS2BACKSCAT ! ! If there is not re-distribution in directions the scattering is just an attenuation ! ATT=ATT*EXP((ICE*VDIR(IS))*DTG) END IF #endif SPEC(1+(IK-1)*NTH:NTH+(IK-1)*NTH) = ATT*SPEC2(1+(IK-1)*NTH:NTH+(IK-1)*NTH) ! ! Second part of ice term integration: scattering including re-distribution in directions ! #ifdef W3_IS2 IF (IS2PARS(2).GE.0) THEN IF (IS2PARS(20).GT.0.5) THEN ! ! Case of isotropic back-scatter: the directional spectrum is decomposed into ! - an isotropic part (ISO): eigenvalue of scattering is 0 ! - the rest (SPEC-ISO): eigenvalue of scattering is VDIR(IS) ! SCAT = EXP(VDIR(IS)*IS2PARS(2)*DTG) ISO = SUM(SPEC(1+(IK-1)*NTH:NTH+(IK-1)*NTH))/NTH SPEC(1+(IK-1)*NTH:NTH+(IK-1)*NTH) = ISO & +(SPEC(1+(IK-1)*NTH:NTH+(IK-1)*NTH)-ISO)*SCAT ELSE ! ! General solution with matrix exponentials: same as bottom scattering, see Ardhuin & Herbers (JFM 2002) ! SCATSPEC(1:NTH)=DBLE(SPEC(1+(IK-1)*NTH:NTH+(IK-1)*NTH)) SPEC(1+(IK-1)*NTH:NTH+(IK-1)*NTH) = & REAL(MATMUL(IS2EIGVEC(:,:), EXP(IS2EIGVAL(:)*VDIR(IS)*DTG*IS2PARS(2)) & *MATMUL(TRANSPOSE(IS2EIGVEC(:,:)),SCATSPEC))) END IF END IF #endif ! ! 10.2 Fluxes of energy and momentum due to ice effects ! FACTOR = DDEN(IK)/CG1(IK) !Jacobian to get energy in band FACTOR2= FACTOR*GRAV*WN1(IK)/SIG(IK) ! coefficient to get momentum DO ITH = 1,NTH IS = ITH+(IK-1)*NTH PHICE = PHICE + (SPEC(IS)-SPEC2(IS)) * FACTOR COSI(1)=ECOS(IS) COSI(2)=ESIN(IS) TAUICE(:) = TAUICE(:) - (SPEC(IS)-SPEC2(IS))*FACTOR2*COSI(:) END DO END DO PHICE =-1.*DWAT*GRAV*PHICE /DTG TAUICE(:)=TAUICE(:)/DTG ELSE #ifdef W3_IS2 IF (IS2PARS(10).LT.0.5) THEN ICEF = 0. ENDIF #endif END IF ! ! ! - - - - - - - - - - - - - - - - - - - - - - ! 11. Sea state dependent stress routine calls ! - - - - - - - - - - - - - - - - - - - - - - !Note the Sea-state dependent stress calculations are primarily for high-wind !conditions (>10 m/s). It is not recommended to use these at lower wind !in their current state. ! #ifdef W3_DEBUGSRC IF (IX .eq. DEBUG_NODE) THEN WRITE(740+IAPROC,*) '4 : sum(SPEC)=', sum(SPEC) END IF #endif ! FLD1/2 requires the calculation of FPI: #ifdef W3_FLD1 CALL CALC_FPI(SPEC, CG1, FPI, VSIN ) #endif #ifdef W3_FLD2 CALL CALC_FPI(SPEC, CG1, FPI, VSIN ) #endif ! #ifdef W3_FLD1 IF (U10ABS.GT.10. .and. HSTOT.gt.0.5) then CALL W3FLD1 ( SPEC,min(FPI/TPI,2.0),COEF*U10ABS*COS(U10DIR), & COEF*U10ABS*Sin(U10DIR), ZWND, DEPTH, 0.0, & DAIR, USTAR, USTDIR, Z0,TAUNUX,TAUNUY,CHARN) ELSE CHARN = AALPHA ENDIF #endif #ifdef W3_FLD2 IF (U10ABS.GT.10. .and. HSTOT.gt.0.5) then CALL W3FLD2 ( SPEC,min(FPI/TPI,2.0),COEF*U10ABS*COS(U10DIR), & COEF*U10ABS*Sin(U10DIR), ZWND, DEPTH, 0.0, & DAIR, USTAR, USTDIR, Z0,TAUNUX,TAUNUY,CHARN) ELSE CHARN = AALPHA ENDIF #endif ! ! 12. includes shoreline reflection --------------------------------------------- * ! #ifdef W3_DEBUGSRC IF (IX .eq. DEBUG_NODE) THEN WRITE(740+IAPROC,*) '5 : sum(SPEC)=', sum(SPEC) END IF #endif #ifdef W3_REF1 IF (REFLEC(1).GT.0.OR.REFLEC(2).GT.0.OR.(REFLEC(4).GT.0.AND.BERG.GT.0)) THEN CALL W3SREF ( SPEC, CG1, WN1, EMEAN, FMEAN, DEPTH, CX, CY, & REFLEC, REFLED, TRNX, TRNY, & BERG, DTG, IX, IY, JSEA, VREF ) IF (GTYPE.EQ.UNGTYPE.AND.REFPARS(3).LT.0.5) THEN #ifdef W3_PDLIB IF (IOBP_LOC(JSEA).EQ.0) THEN #else IF (IOBP(IX).EQ.0) THEN #endif DO IK=1, NK DO ITH=1, NTH ISP = ITH+(IK-1)*NTH #ifdef W3_PDLIB IF (IOBPD_LOC(ITH,JSEA).EQ.0) SPEC(ISP) = DTG*VREF(ISP) #else IF (IOBPD(ITH,IX).EQ.0) SPEC(ISP) = DTG*VREF(ISP) #endif END DO END DO ELSE SPEC(:) = SPEC(:) + DTG * VREF(:) ENDIF ELSE SPEC(:) = SPEC(:) + DTG * VREF(:) END IF END IF #endif ! #ifdef W3_DEBUGSRC IF (IX .eq. DEBUG_NODE) THEN WRITE(740+IAPROC,*) '6 : sum(SPEC)=', sum(SPEC) END IF #endif FIRST = .FALSE. IF (IT.EQ.0) SPEC = SPECINIT SPEC = MAX(0., SPEC) ! RETURN ! ! Formats ! #ifdef W3_NNT 8000 FORMAT (/' *** ERROR W3SRCE : ERROR IN OPENING FILE ',A,' ***'/ & ' IOSTAT = ',I10/) 8001 FORMAT (/' *** ERROR W3SRCE : ERROR IN WRITING TO FILE ***'/ & ' IOSTAT = ',I10/) #endif ! #ifdef W3_T 9000 FORMAT (' TEST W3SRCE : COUNTERS : NO LONGER AVAILABLE') 9001 FORMAT (' TEST W3SRCE : DEPTH :',F8.1/ & ' WIND SPEED :',F8.1/ & ' WIND DIR :',F8.1) #endif #ifdef W3_ST1 9004 FORMAT (' TEST W3SRCE : FHIGH (3X) : ',3F8.4/ & ' ------------- NEW DYNAMIC INTEGRATION LOOP', & ' ------------- ') #endif #ifdef W3_ST2 9005 FORMAT (' TEST W3SRCE : FHIGH : ',F8.4/ & ' ------------- NEW DYNAMIC INTEGRATION LOOP', & ' ------------- ') #endif #ifdef W3_ST3 9006 FORMAT (' TEST W3SRCE : FHIGH (3X) : ',3F8.4/ & ' ------------- NEW DYNAMIC INTEGRATION LOOP', & ' ------------- ') #endif #ifdef W3_ST4 9006 FORMAT (' TEST W3SRCE : FHIGH (3X) : ',3F8.4/ & ' ------------- NEW DYNAMIC INTEGRATION LOOP', & ' ------------- ') #endif ! #ifdef W3_T 9020 FORMAT (' TEST W3SRCE : NSTEP : ',I4,' DTTOT :',F6.1) 9021 FORMAT (' TEST W3SRCE : NKH (3X) : ',2I3,I6) 9040 FORMAT (' TEST W3SRCE : DTRAW, DT, SHAVE :',2F6.1,2X,L1) #endif ! #ifdef W3_ST1 9060 FORMAT (' TEST W3SRCE : FHIGH (3X) : ',3F8.4/ & ' NKH : ',I3) #endif #ifdef W3_ST2 9061 FORMAT (' TEST W3SRCE : FHIGH (2X) : ',2F8.4/ & ' NKH : ',I3) #endif #ifdef W3_ST3 9062 FORMAT (' TEST W3SRCE : FHIGH (3X) : ',3F8.4/ & ' NKH : ',I3) #endif #ifdef W3_ST4 9062 FORMAT (' TEST W3SRCE : FHIGH (3X) : ',3F8.4/ & ' NKH : ',I3) #endif #ifdef W3_ST6 9063 FORMAT (' TEST W3SRCE : FHIGH : ',F8.4/ & ' NKH : ',I3) #endif !/ !/ End of W3SRCE ----------------------------------------------------- / !/ END SUBROUTINE W3SRCE !/ ------------------------------------------------------------------- / !> !> @brief Calculate equivalent peak frequency. !> !> @details Tolman and Chalikov (1996), equivalent peak frequency from source. !> !> @param[in] A Action density spectrum (1-D). !> @param[in] CG Group velocities for k-axis of spectrum. !> @param[out] FPI Input 'peak' frequency. !> @param[in] S Source term (1-D version). !> !> @author Jessica Meixner !> @date 06-Jun-2018 !> SUBROUTINE CALC_FPI( A, CG, FPI, S ) !/ !/ +-----------------------------------+ !/ | WAVEWATCH III NOAA/NCEP | !/ | Jessica Meixner | !/ | | !/ | FORTRAN 90 | !/ | Last update : 06-Jun-2018 | !/ +-----------------------------------+ !/ !/ 06-Jul-2016 : Origination ( version 5.12 ) !/ 06-Jul-2016 : Add SUBROUTINE SIGN_VSD_SEMI_IMPLICIT_WW3 !/ Add optional DEBUGSRC/PDLIB ( version 6.04 ) !/ ! 1. Purpose : ! ! Calculate equivalent peak frequency ! ! 2. Method : ! ! Tolman and Chalikov (1996), equivalent peak frequency from source ! 3. Parameters : ! ! Parameter list ! ---------------------------------------------------------------- ! A R.A. I Action density spectrum (1-D). ! CG R.A. I Group velocities for k-axis of spectrum. ! FPI R.A. O Input 'peak' frequency. ! S R.A. I Source term (1-D version). ! ---------------------------------------------------------------- ! ! 4. Subroutines used : ! ! Name Type Module Description ! ---------------------------------------------------------------- ! STRACE Subr. W3SERVMD Subroutine tracing. ! ---------------------------------------------------------------- ! ! 5. Called by : ! ! Name Type Module Description ! ---------------------------------------------------------------- ! W3SRCE Subr. ! ---------------------------------------------------------------- ! ! 6. Error messages : ! ! 7. Remarks : ! ! 8. Structure : ! ! See source code. ! ! 9. Switches : ! ! !/S Enable subroutine tracing. ! ! 10. Source code : ! !/ ------------------------------------------------------------------- / USE CONSTANTS USE W3GDATMD, ONLY: NK, NTH, NSPEC, XFR, DDEN, SIG,FTE, FTTR #ifdef W3_S USE W3SERVMD, ONLY: STRACE #endif ! IMPLICIT NONE !/ !/ ------------------------------------------------------------------- / !/ Parameter list !/ REAL, INTENT(IN) :: A(NSPEC), CG(NK), S(NSPEC) REAL, INTENT(OUT) :: FPI !/ !/ ------------------------------------------------------------------- / !/ Local parameters !/ INTEGER :: IS, IK #ifdef W3_S INTEGER, SAVE :: IENT = 0 #endif REAL :: M0, M1, SIN1A(NK) !/ !/ ------------------------------------------------------------------- / !/ #ifdef W3_S CALL STRACE (IENT, 'CALC_FPI') #endif ! ! Calculate FPI: equivalent peak frequncy from wind source term ! input ! DO IK=1, NK SIN1A(IK) = 0. DO IS=(IK-1)*NTH+1, IK*NTH SIN1A(IK) = SIN1A(IK) + MAX ( 0. , S(IS) ) END DO END DO ! M0 = 0. M1 = 0. DO IK=1, NK SIN1A(IK) = SIN1A(IK) * DDEN(IK) / ( CG(IK) * SIG(IK)**3 ) M0 = M0 + SIN1A(IK) M1 = M1 + SIN1A(IK)/SIG(IK) END DO ! SIN1A(NK) = SIN1A(NK) / DDEN(NK) M0 = M0 + SIN1A(NK) * FTE M1 = M1 + SIN1A(NK) * FTTR IF ( M1 .LT. 1E-20 ) THEN FPI = XFR * SIG(NK) ELSE FPI = M0 / M1 END IF END SUBROUTINE CALC_FPI !/ ------------------------------------------------------------------- /! !> !> @brief Put source term in matrix same as done always. !> !> @param[in] SPEC !> @param[inout] VS !> @param[inout] VD !> !> @author Aron Roland !> @author Mathieu Dutour-Sikiric !> @date 01-Jun-2018 !> SUBROUTINE SIGN_VSD_SEMI_IMPLICIT_WW3(SPEC, VS, VD) !/ !/ +-----------------------------------+ !/ | WAVEWATCH III NOAA/NCEP | !/ | | !/ | Aron Roland (BGS IT&E GmbH) | !/ | Mathieu Dutour-Sikiric (IRB) | !/ | | !/ | FORTRAN 90 | !/ | Last update : 01-June-2018 | !/ +-----------------------------------+ !/ !/ 01-June-2018 : Origination. ( version 6.04 ) !/ ! 1. Purpose : Put source term in matrix same as done always ! 2. Method : ! 3. Parameters : ! ! Parameter list ! ---------------------------------------------------------------- ! ---------------------------------------------------------------- ! ! 4. Subroutines used : ! ! Name Type Module Description ! ---------------------------------------------------------------- ! STRACE Subr. W3SERVMD Subroutine tracing. ! ---------------------------------------------------------------- ! ! 5. Called by : ! ! Name Type Module Description ! ---------------------------------------------------------------- ! ---------------------------------------------------------------- ! ! 6. Error messages : ! 7. Remarks ! 8. Structure : ! 9. Switches : ! ! !/S Enable subroutine tracing. ! ! 10. Source code : ! !/ ------------------------------------------------------------------- / #ifdef W3_S USE W3SERVMD, ONLY: STRACE #endif ! USE W3GDATMD, only : NTH, NK, NSPEC IMPLICIT NONE !/ !/ ------------------------------------------------------------------- / !/ Parameter list !/ !/ ------------------------------------------------------------------- / !/ Local PARAMETERs !/ #ifdef W3_S INTEGER, SAVE :: IENT = 0 #endif !/ !/ ------------------------------------------------------------------- / !/ INTEGER :: ISP, ITH, IK, IS REAL, INTENT(IN) :: SPEC(NSPEC) REAL, INTENT(INOUT) :: VS(NSPEC), VD(NSPEC) #ifdef W3_S CALL STRACE (IENT, 'SIGN_VSD_SEMI_IMPLICIT_WW3') #endif DO IS=1,NSPEC VD(IS) = MIN(0., VD(IS)) END DO END SUBROUTINE SIGN_VSD_SEMI_IMPLICIT_WW3 !/ ------------------------------------------------------------------- / !> !> @brief Put source term in matrix Patankar style (experimental). !> !> @param[in] SPEC !> @param[inout] VS !> @param[inout] VD !> !> @author Aron Roland !> @author Mathieu Dutour-Sikiric !> @date 01-Jun-2018 !> SUBROUTINE SIGN_VSD_PATANKAR_WW3(SPEC, VS, VD) !/ !/ +-----------------------------------+ !/ | WAVEWATCH III NOAA/NCEP | !/ | | !/ | Aron Roland (BGS IT&E GmbH) | !/ | Mathieu Dutour-Sikiric (IRB) | !/ | | !/ | FORTRAN 90 | !/ | Last update : 01-June-2018 | !/ +-----------------------------------+ !/ !/ 01-June-2018 : Origination. ( version 6.04 ) !/ ! 1. Purpose : Put source term in matrix Patankar style (experimental) ! 2. Method : ! 3. Parameters : ! ! Parameter list ! ---------------------------------------------------------------- ! ---------------------------------------------------------------- ! ! 4. Subroutines used : ! ! Name Type Module Description ! ---------------------------------------------------------------- ! STRACE Subr. W3SERVMD Subroutine tracing. ! ---------------------------------------------------------------- ! ! 5. Called by : ! ! Name Type Module Description ! ---------------------------------------------------------------- ! ---------------------------------------------------------------- ! ! 6. Error messages : ! 7. Remarks ! 8. Structure : ! 9. Switches : ! ! !/S Enable subroutine tracing. ! ! 10. Source code : ! !/ ------------------------------------------------------------------- / #ifdef W3_S USE W3SERVMD, ONLY: STRACE #endif ! USE W3GDATMD, only : NTH, NK, NSPEC IMPLICIT NONE !/ !/ ------------------------------------------------------------------- / !/ Parameter list !/ !/ ------------------------------------------------------------------- / !/ Local PARAMETERs !/ #ifdef W3_S INTEGER, SAVE :: IENT = 0 #endif !/ !/ ------------------------------------------------------------------- / !/ INTEGER :: ISP, ITH, IK, IS REAL, INTENT(IN) :: SPEC(NSPEC) REAL, INTENT(INOUT) :: VS(NSPEC), VD(NSPEC) #ifdef W3_S CALL STRACE (IENT, 'SIGN_VSD_PATANKAR_WW3') #endif DO IS=1,NSPEC VD(IS) = MIN(0., VD(IS)) VS(IS) = MAX(0., VS(IS)) END DO END SUBROUTINE SIGN_VSD_PATANKAR_WW3 !/ !/ End of module W3SRCEMD -------------------------------------------- / !/ END MODULE W3SRCEMD