WAVEWATCH III  beta 0.0.1
w3srcemd.F90
Go to the documentation of this file.
1 
8 
9 #include "w3macros.h"
10 !/ ------------------------------------------------------------------- /
11 
24 MODULE w3srcemd
25  !/
26  !/ +-----------------------------------+
27  !/ | WAVEWATCH III NOAA/NCEP |
28  !/ | H. L. Tolman |
29  !/ | F. Ardhuin |
30  !/ | FORTRAN 90 |
31  !/ | Last update : 22-Mar-2021 |
32  !/ +-----------------------------------+
33  !/
34  !/ For updates see subroutine.
35  !/
36  ! 1. Purpose :
37  !
38  ! Source term integration routine.
39  !
40  ! 2. Variables and types :
41  !
42  ! Name Type Scope Description
43  ! ----------------------------------------------------------------
44  ! OFFSET R.P. Private Offset in time integration scheme.
45  ! 0.5 in original WAM, now 1.0
46  ! ----------------------------------------------------------------
47  !
48  ! 3. Subroutines and functions :
49  !
50  ! Name Type Scope Description
51  ! ----------------------------------------------------------------
52  ! W3SRCE Subr. Public Calculate and integrate source terms.
53  ! ----------------------------------------------------------------
54  !
55  ! 4. Subroutines and functions used :
56  !
57  ! See corresponding documentation of W3SRCE.
58  !
59  ! 5. Remarks :
60  !
61  ! 6. Switches :
62  !
63  ! See section 9 of W3SRCE.
64  !
65  ! 7. Source code :
66  !
67  !/ ------------------------------------------------------------------- /
68  !/
69  REAL, PARAMETER, PRIVATE:: OFFSET = 1.
70  !/
71 CONTAINS
72  !/ ------------------------------------------------------------------- /
73 
190  SUBROUTINE w3srce ( srce_call, IT, ISEA, JSEA, IX, IY, IMOD, &
191  SPECOLD, SPEC, VSIO, VDIO, SHAVEIO, &
192  ALPHA, WN1, CG1, CLATSL, &
193  D_INP, U10ABS, U10DIR, &
194 #ifdef W3_FLX5
195  TAUA, TAUADIR, &
196 #endif
197  AS, USTAR, USTDIR, &
198  CX, CY, ICE, ICEH, ICEF, ICEDMAX, &
199  REFLEC, REFLED, DELX, DELY, DELA, TRNX, &
200  TRNY, BERG, FPI, DTDYN, FCUT, DTG, TAUWX, &
201  TAUWY, TAUOX, TAUOY, TAUWIX, TAUWIY, TAUWNX,&
202  TAUWNY, PHIAW, CHARN, TWS, PHIOC, WHITECAP, &
203  D50, PSIC, BEDFORM , PHIBBL, TAUBBL, TAUICE,&
204  PHICE, TAUOCX, TAUOCY, WNMEAN, DAIR, COEF)
205  !/
206  !/ +-----------------------------------+
207  !/ | WAVEWATCH III NOAA/NCEP |
208  !/ | H. L. Tolman |
209  !/ | F. Ardhuin |
210  !/ | A. Roland |
211  !/ | M. Dutour Sikiric |
212  !/ | FORTRAN 90 |
213  !/ | Last update : 22-Mar-2021 |
214  !/ +-----------------------------------+
215  !/
216  !/ 06-Dec-1996 : Final FORTRAN 77 ( version 1.18 )
217  !/ 04-Feb-2000 : Upgrade to FORTRAN 90 ( version 2.00 )
218  !/ 14-Feb-2000 : Exact-NL added ( version 2.01 )
219  !/ 04-May-2000 : Non-central integration ( version 2.03 )
220  !/ 02-Feb-2001 : Xnl version 3.0 ( version 2.07 )
221  !/ 09-May-2002 : Switch clean up. ( version 2.21 )
222  !/ 13-Nov-2002 : Add stress vector. ( version 3.00 )
223  !/ 27-Nov-2002 : First version of VDIA and MDIA. ( version 3.01 )
224  !/ 07-Oct-2003 : Output options for NN training. ( version 3.05 )
225  !/ 24-Dec-2004 : Multiple model version. ( version 3.06 )
226  !/ 23-Jun-2006 : Linear input added. ( version 3.09 )
227  !/ 27-Jun-2006 : Adding file name preamble. ( version 3.09 )
228  !/ 04-Jul-2006 : Separation of stress computation. ( version 3.09 )
229  !/ 16-Apr-2007 : Miche style limiter added. ( version 3.11 )
230  !/ (J. H. Alves)
231  !/ 25-Apr-2007 : Battjes-Janssen Sdb added. ( version 3.11 )
232  !/ (J. H. Alves)
233  !/ 09-Oct-2007 : Adding WAM 4+ and SB1 options. ( version 3.13 )
234  !/ (F. Ardhuin)
235  !/ 29-May-2009 : Preparing distribution version. ( version 3.14 )
236  !/ 19-Aug-2010 : Making treatment of 0 water depth ( version 3.14.6 )
237  !/ consistent with the rest of the model.
238  !/ 31-Mar-2010 : Adding ice conc. and reflections ( version 3.14.4 )
239  !/ 15-May-2010 : Adding transparencies ( version 3.14.4 )
240  !/ 01-Jun-2011 : Movable bed bottom friction in BT4 ( version 4.01 )
241  !/ 01-Jul-2011 : Energy and momentum flux, friction ( version 4.01 )
242  !/ 24-Aug-2011 : Uses true depth for depth-induced ( version 4.04 )
243  !/ 16-Sep-2011 : Initialization of TAUWAX, TAUWAY ( version 4.04 )
244  !/ 1-Dec-2011 : Adding BYDRZ source term package ( version 4.04 )
245  !/ ST6 and optional Hwang (2011)
246  !/ stresses FLX4.
247  !/ 14-Mar-2012 : Update of BT4, passing PSIC ( version 4.04 )
248  !/ 13-Jul-2012 : Move GMD (SNL3) and nonlinear filter (SNLS)
249  !/ from 3.15 (HLT). ( version 4.08 )
250  !/ 28-Aug-2013 : Corrected MLIM application ( version 4.11 )
251  !/ 10-Sep-2013 : Special treatment for IG band ( version 4.15 )
252  !/ 14-Nov-2013 : Make orphaned pars in par lst local ( version 4.13 )
253  !/ 17-Nov-2013 : Coupling fraction of ice-free ( version 4.13 )
254  !/ surface to SIN and SDS. (S. Zieger)
255  !/ 01-Avr-2014 : Adding ice thickness and floe size ( version 4.18 )
256  !/ 23-May-2014 : Adding ice fluxes to W3SRCE ( version 5.01 )
257  !/ 27-Aug-2015 : Adding inputs to function W3SIS2 ( version 5.10 )
258  !/ 13-Dec-2015 : Implicit integration of Sice (F.A.) ( version 5.10 )
259  !/ 30-Jul-2017 : Adds TWS in interface ( version 6.04 )
260  !/ 07-Jan-2018 : Allows variable ice scaling (F.A.) ( version 6.04 )
261  !/ 01-Jan-2018 : Add implicit source term integration ( version 6.04)
262  !/ 01-Jan-2018 : within PDLIB (A. Roland, M. Dutour
263  !/ 18-Aug-2018 : S_{ice} IC5 (Q. Liu) ( version 6.06)
264  !/ 26-Aug-2018 : UOST (Mentaschi et al. 2015, 2018) ( version 6.06 )
265  !/ 22-Mar-2021 : Add extra fields used in coupling ( version 7.13 )
266  !/ 07-Jun-2021 : S_{nl5} GKE NL5 (Q. Liu) ( version 7.13 )
267  !/ 19-Jul-2021 : Momentum and air density support ( version 7.14 )
268  !/
269  !/ Copyright 2009-2013 National Weather Service (NWS),
270  !/ National Oceanic and Atmospheric Administration. All rights
271  !/ reserved. WAVEWATCH III is a trademark of the NWS.
272  !/ No unauthorized use without permission.
273  !/
274  ! 1. Purpose :
275  !
276  ! Calculate and integrate source terms for a single grid point.
277  !
278  ! 2. Method :
279  !
280  ! Physics : see manual and corresponding subroutines.
281  !
282  ! Numerics :
283  !
284  ! Dynamic-implicit integration of the source terms based on
285  ! WW-II (Tolman 1992). The dynamic time step is calculated
286  ! given a maximum allowed change of spectral densities for
287  ! frequencies / wavenumbers below the usual cut-off.
288  ! The maximum change is given by the minimum of a parametric
289  ! and a relative change. The parametric change relates to a
290  ! PM type equilibrium range
291  !
292  ! -1 (2pi)**4 1
293  ! dN(k) = Xp alpha pi ---------- ------------
294  ! max g**2 k**3 sigma
295  !
296  ! 1 .
297  ! = FACP ------------ (1)
298  ! k**3 sigma .
299  !
300  ! where
301  ! alpha = 0.62e-4 (set in W3GRID)
302  ! Xp fraction of PM shape (read in W3GRID)
303  ! FACP combined factor (set in W3GRID)
304  !
305  ! The maximum relative change is given as
306  !
307  ! / +- -+ \ .
308  ! dN(k) = Xr max | N(k) , max | Nx , Xfilt N(k) | | (2)
309  ! max \ +- max-+ / .
310  !
311  ! where
312  ! Xr fraction of relative change (read in W3GRID)
313  ! Xfilt filter level (read in W3GRID)
314  ! Nx Maximum parametric change (1)
315  ! for largest wavenumber.
316  !
317  ! 3. Parameters :
318  !
319  ! Parameter list
320  ! ----------------------------------------------------------------
321  ! IX,IY Int. I Discrete grid point counters.
322  ! IMOD Int. I Model number.
323  ! SPEC R.A. I/O Spectrum (action) in 1-D form.
324  ! ALPHA R.A. I/O Nondimenional 1-D spectrum corresponding
325  ! to above full spectra (Phillip's const.).
326  ! Calculated separately for numerical
327  ! economy on vector machine (W3SPR2).
328  ! WN1 R.A. I Discrete wavenumbers.
329  ! CG1 R.A. I Id. group velocities.
330  ! D_INP Real. I Depth. Compared to DMIN to get DEPTH.
331  ! U10ABS Real. I Wind speed at reference height.
332  ! U10DIR Real. I Id. wind direction.
333  ! TAUA Real. I Magnitude of total atmospheric stress ( !/FLX5 )
334  ! TAUADIR Real. I Direction of atmospheric stress ( !/FLX5 )
335  ! AS Real. I Air-sea temp. difference. ( !/ST3 )
336  ! USTAR Real. !/O Friction velocity.
337  ! USTDIR Real !/O Idem, direction.
338  ! CX-Y Real. I Current velocity components. ( !/BS1 )
339  ! ICE Real I Sea ice concentration
340  ! ICEH Real I Sea ice thickness
341  ! ICEF Real I/O Sea ice maximum floe diameter (updated)
342  ! ICEDMAX Real I/O Sea ice maximum floe diameter
343  ! BERG Real I Iceberg damping coefficient ( !/BS1 )
344  ! REFLEC R.A. I reflection coefficients ( !/BS1 )
345  ! REFLED I.A. I reflection direction ( !/BS1 )
346  ! TRNX-Y Real I Grid transparency in X and Y ( !/BS1 )
347  ! DELX Real. I grid cell size in X direction ( !/BS1 )
348  ! DELY Real. I grid cell size in Y direction ( !/BS1 )
349  ! DELA Real. I grid cell area ( !/BS1 )
350  ! FPI Real I/O Peak-input frequency. ( !/ST2 )
351  ! WHITECAP R.A. O Whitecap statisics ( !/ST4 )
352  ! DTDYN Real O Average dynamic time step.
353  ! FCUT Real O Cut-off frequency for tail.
354  ! DTG Real I Global time step.
355  ! D50 Real I Sand grain size ( !/BT4 )
356  ! BEDFORM R.A. I/O Bedform parameters ( !/BT4 )
357  ! PSIC Real I Critical Shields ( !/BT4 )
358  ! PHIBBL Real O Energy flux to BBL ( !/BTx )
359  ! TAUBBL R.A. O Momentum flux to BBL ( !/BTx )
360  ! TAUICE R.A. O Momentum flux to sea ice ( !/ICx )
361  ! PHICE Real O Energy flux to sea ice ( !/ICx )
362  ! TAUOCX-YReal O Total ocean momentum components
363  ! WNMEAN Real O Mean wave number
364  ! DAIR Real I Air density
365  ! ----------------------------------------------------------------
366  ! Note: several pars are set to I/O to avoid compiler warnings.
367  !
368  ! 4. Subroutines used :
369  !
370  ! Name Type Module Description
371  ! ----------------------------------------------------------------
372  ! W3SPRn Subr. W3SRCnMD Mean wave parameters for use in
373  ! source terms.
374  ! W3FLXn Subr. W3FLXnMD Flux/stress computation.
375  ! W3SLNn Subr. W3SLNnMD Linear input.
376  ! W3SINn Subr. W3SRCnMD Input source term.
377  ! W3SNLn Subr. W3SNLnMD Nonlinear interactions.
378  ! W3SNLS Subr. W3SNLSMD Nonlinear smoother.
379  ! W3SDSn Subr. W3SRCnMD Whitecapping source term
380  ! W3SBTn Subr. W3SBTnMD Bottom friction source term.
381  ! W3SDBn Subr. W3SBTnMD Depth induced breaking source term.
382  ! W3STRn Subr. W3STRnMD Triad interaction source term.
383  ! W3SBSn Subr. W3SBSnMD Bottom scattering source term.
384  ! W3REFn Subr. W3REFnMD Reflexions (shore, icebergs ...).
385  ! STRACE Subr. W3SERVMD Subroutine tracing (!/S)
386  ! ----------------------------------------------------------------
387  !
388  ! 5. Called by :
389  !
390  ! Name Type Module Description
391  ! ----------------------------------------------------------------
392  ! W3WAVE Subr. W3WAVEMD Actual wave model routine.
393  ! ----------------------------------------------------------------
394  !
395  ! 6. Error messages :
396  !
397  ! None.
398  !
399  ! 7. Remarks :
400  !
401  ! - No testing is performed on the status of the grid point.
402  !
403  ! 8. Structure :
404  !
405  ! -----------------------------------------------------------------
406  ! 1. Preparations
407  ! a Set maximum change and wavenumber arrays.
408  ! b Prepare dynamic time stepping.
409  ! c Compute mean parameters. ( W3SPRn )
410  ! d Compute stresses (if posible).
411  ! e Prepare cut-off
412  ! f Test output for !/NNT option.
413  ! --start-dynamic-integration-loop---------------------------------
414  ! 2. Calculate source terms
415  ! a Input. ( W3SLNx, W3SINn )
416  ! b Nonlinear interactions. ( W3SNLn )
417  ! c Dissipation ( W3SDSn )
418  ! 1 as included in source terms ( W3SDSn )
419  ! 2 optional dissipation due to different physics ( W3SWLn )
420  ! d Bottom friction. ( W3SBTn )
421  ! 3. Calculate cut-off frequencie(s)
422  ! 4. Summation of source terms and diagonal term and time step.
423  ! 5. Increment spectrum.
424  ! 6. Add tail
425  ! a Mean wave parameters and cut-off ( W3SPRn )
426  ! b 'Seeding' of spectrum. ( !/SEED )
427  ! c Add tail
428  ! 7. Check if integration complete.
429  ! --end-dynamic-integration-loop-----------------------------------
430  ! 8. Save integration data.
431  ! -----------------------------------------------------------------
432  !
433  ! 9. Switches :
434  !
435  ! !/FLX1 Wu (1980) stress computation. ( Choose one )
436  ! !/FLX2 T&C (1996) stress computation.
437  ! !/FLX3 T&C (1996) stress computation with cap.
438  ! !/FLX4 Hwang (2011) stress computation (2nd order).
439  ! !/FLX5 Direct use of stress from atmoshperic model.
440  !
441  ! !/LN0 No linear input. ( Choose one )
442  !
443  ! !/ST0 No input and dissipation. ( Choose one )
444  ! !/ST1 WAM-3 input and dissipation.
445  ! !/ST2 Tolman and Chalikov (1996) input and dissipation.
446  ! !/ST3 WAM 4+ input and dissipation.
447  ! !/ST4 Ardhuin et al. (2009, 2010)
448  ! !/ST6 BYDB source terms after Babanin, Young, Donelan and Banner.
449  !
450  ! !/NL0 No nonlinear interactions. ( Choose one )
451  ! !/NL1 Discrete interaction approximation.
452  ! !/NL2 Exact nonlinear interactions.
453  ! !/NL3 Generalized Multiple DIA.
454  ! !/NL4 Two Scale Approximation
455  ! !/NL5 Generalized Kinetic Equation.
456  ! !/NLS Nonlinear HF smoother.
457  !
458  ! !/BT0 No bottom friction. ( Choose one )
459  ! !/BT1 JONSWAP bottom friction.
460  ! !/BT4 Bottom friction using movable bed roughness
461  ! (Tolman 1994, Ardhuin & al. 2003)
462  ! !/BT8 Muddy bed (Dalrymple & Liu).
463  ! !/BT9 Muddy bed (Ng).
464  !
465  ! !/IC1 Dissipation via interaction with ice according to simple
466  ! methods: 1) uniform in frequency or
467  ! !/IC2 2) Liu et al. model
468  ! !/IC3 Dissipation via interaction with ice according to a
469  ! viscoelastic sea ice model (Wang and Shen 2010).
470  ! !/IC4 Dissipation via interaction with ice as a function of freq.
471  ! (empirical/parametric methods)
472  ! !/IC5 Dissipation via interaction with ice according to a
473  ! viscoelastic sea ice model (Mosig et al. 2015).
474  ! !/DB0 No depth-limited breaking. ( Choose one )
475  ! !/DB1 Battjes-Janssen depth-limited breaking.
476  !
477  ! !/TR0 No triad interactions. ( Choose one )
478  ! !/TR1 Lumped Triad Approximation (LTA).
479  !
480  ! !/BS0 No bottom scattering. ( Choose one )
481  ! !/BS1 Scattering term by Ardhuin and Magne (2007).
482  !
483  ! !/MLIM Miche style limiter for shallow water and steepness.
484  !
485  ! !/SEED 'Seeding' of lowest frequency for suffuciently strong
486  ! winds.
487  !
488  ! !/NNT Write output to file test_data_NNN.ww3 for NN training.
489  !
490  ! !/S Enable subroutine tracing.
491  ! !/T Enable general test output.
492  !
493  ! 10. Source code :
494  !
495  !/ ------------------------------------------------------------------- /
496  USE constants, ONLY: dwat, srce_imp_post, srce_imp_pre, &
497  srce_direct, grav, tpi, tpiinv
498  USE w3gdatmd, ONLY: nk, nth, nspec, sig, th, dmin, dtmax, &
499  dtmin, facti1, facti2, facsd, fachfa, facp, &
500  xfc, xflt, xrel, xft, fxfm, fxpm, dden, &
501  fte, ftf, fhmax, ecos, esin, iicedisp, &
502  icescales, iicesmooth
503  USE w3wdatmd, ONLY: time
504  USE w3odatmd, ONLY: ndse, ndst, iaproc
505  USE w3idatmd, ONLY: inflags2
506  USE w3dispmd
507 #ifdef W3_T
508  USE constants, ONLY: rade
509 #endif
510 #ifdef W3_REF1
511  USE w3gdatmd, ONLY: iobp, iobpd, gtype, ungtype, refpars
512 #endif
513 #ifdef W3_NNT
514  USE w3odatmd, ONLY: iaproc, screen, fnmpre
515 #endif
516 #ifdef W3_FLD1
517  USE w3fld1md, ONLY: w3fld1
518  USE w3gdatmd, ONLY: aalpha
519 #endif
520 #ifdef W3_FLD2
521  USE w3fld2md, ONLY: w3fld2
522  USE w3gdatmd, ONLY: aalpha
523 #endif
524 #ifdef W3_FLX1
525  USE w3flx1md
526 #endif
527 #ifdef W3_FLX2
528  USE w3flx2md
529 #endif
530 #ifdef W3_FLX3
531  USE w3flx3md
532 #endif
533 #ifdef W3_FLX4
534  USE w3flx4md
535 #endif
536 #ifdef W3_FLX5
537  USE w3flx5md
538 #endif
539 #ifdef W3_LN1
540  USE w3sln1md
541 #endif
542 #ifdef W3_ST0
543  USE w3src0md
544 #endif
545 #ifdef W3_ST1
546  USE w3src1md
547 #endif
548 #ifdef W3_ST2
549  USE w3src2md
550  USE w3gdatmd, ONLY : zwind
551 #endif
552 #ifdef W3_ST3
553  USE w3src3md
554  USE w3gdatmd, ONLY : zzwnd, ffxfm, ffxpm
555 #endif
556 #ifdef W3_ST4
557  USE w3src4md, ONLY : w3spr4, w3sin4, w3sds4
558  USE w3gdatmd, ONLY : zzwnd, ffxfm, ffxpm, ffxfa, sintailpar
559 #endif
560 #ifdef W3_ST6
561  USE w3src6md
562  USE w3swldmd, ONLY : w3swl6
563  USE w3gdatmd, ONLY : swl6s6
564 #endif
565 #ifdef W3_NL1
566  USE w3snl1md
567  USE w3gdatmd, ONLY: iqtpe
568 #endif
569 #ifdef W3_NL2
570  USE w3snl2md
571 #endif
572 #ifdef W3_NL3
573  USE w3snl3md
574 #endif
575 #ifdef W3_NL4
576  USE w3snl4md
577 #endif
578 #ifdef W3_NL5
579  USE w3snl5md
580  USE w3timemd, ONLY: tick21
581 #endif
582 #ifdef W3_NLS
583  USE w3snlsmd
584 #endif
585 #ifdef W3_BT1
586  USE w3sbt1md
587 #endif
588 #ifdef W3_BT4
589  USE w3sbt4md
590 #endif
591 #ifdef W3_BT8
592  USE w3sbt8md
593 #endif
594 #ifdef W3_BT9
595  USE w3sbt9md
596 #endif
597 #ifdef W3_IC1
598  USE w3sic1md
599 #endif
600 #ifdef W3_IC2
601  USE w3sic2md
602 #endif
603 #ifdef W3_IC3
604  USE w3sic3md
605 #endif
606 #ifdef W3_IC4
607  USE w3sic4md
608 #endif
609 #ifdef W3_IC5
610  USE w3sic5md
611 #endif
612 #ifdef W3_IS1
613  USE w3sis1md
614 #endif
615 #ifdef W3_IS2
616  USE w3sis2md
617  USE w3gdatmd, ONLY : is2pars
618 #endif
619 #ifdef W3_DB1
620  USE w3sdb1md
621 #endif
622 #ifdef W3_TR1
623  USE w3str1md
624 #endif
625 #ifdef W3_BS1
626  USE w3sbs1md
627 #endif
628 #ifdef W3_REF1
629  USE w3ref1md
630 #endif
631 #ifdef W3_IG1
632  USE w3gdatmd, ONLY : igpars
633 #endif
634 #ifdef W3_S
635  USE w3servmd, ONLY: strace
636 #endif
637 #ifdef W3_NNT
638  USE w3servmd, ONLY: extcde
639 #endif
640 #ifdef W3_UOST
641  USE w3uostmd, ONLY: uost_srctrmcompute
642 #endif
643 #ifdef W3_PDLIB
644  USE pdlib_w3profsmd, ONLY : b_jac, aspar_jac, aspar_diag_all
645  USE yownodepool, ONLY: pdlib_i_diag, pdlib_si
646  USE w3gdatmd, ONLY: b_jgs_limiter, fssource, optioncall
647  USE w3gdatmd, ONLY: iobp_loc, iobpd_loc, b_jgs_limiter_func
648  USE w3wdatmd, ONLY: va
649  USE w3parall, ONLY: imem, lsloc
650 #endif
651  !/
652  IMPLICIT NONE
653  !/
654  !/ ------------------------------------------------------------------- /
655  !/ Parameter list
656  !/
657  INTEGER, INTENT(IN) :: srce_call, IT, ISEA, JSEA, IX, IY, IMOD
658  REAL, intent(in) :: SPECOLD(NSPEC), CLATSL
659  REAL, INTENT(OUT) :: VSIO(NSPEC), VDIO(NSPEC)
660  LOGICAL, INTENT(OUT) :: SHAVEIO
661  REAL, INTENT(IN) :: D_INP, U10ABS, &
662  U10DIR, AS, CX, CY, DTG, D50,PSIC, &
663  ICE, ICEH
664 #ifdef W3_FLX5
665  REAL, INTENT(IN) :: TAUA, TAUADIR
666 #endif
667  INTEGER, INTENT(IN) :: REFLED(6)
668  REAL, INTENT(IN) :: REFLEC(4), DELX, DELY, DELA, &
669  TRNX, TRNY, BERG, ICEDMAX, DAIR
670  REAL, INTENT(INOUT) :: WN1(NK), CG1(NK), &
671  SPEC(NSPEC), ALPHA(NK), USTAR, &
672  USTDIR, FPI, TAUOX, TAUOY, &
673  TAUWX, TAUWY, PHIAW, PHIOC, PHICE, &
674  CHARN, TWS, BEDFORM(3), PHIBBL, &
675  TAUBBL(2), TAUICE(2), WHITECAP(4), &
676  TAUWIX, TAUWIY, TAUWNX, TAUWNY, &
677  ICEF, TAUOCX, TAUOCY, WNMEAN
678  REAL, INTENT(OUT) :: DTDYN, FCUT
679  REAL, INTENT(IN) :: COEF
680  !/
681  !/ ------------------------------------------------------------------- /
682  !/ Local parameters
683  !/
684  INTEGER :: IK, ITH, IS, IS0, NSTEPS, NKH, NKH1, &
685  IKS1, IS1, NSPECH, IDT, IERR, ISP
686  REAL :: DTTOT, FHIGH, DT, AFILT, DAMAX, AFAC, &
687  HDT, ZWND, FP, DEPTH, TAUSCX, TAUSCY, FHIGI
688  ! Scaling factor for SIN, SDS, SNL
689  REAL :: ICESCALELN, ICESCALEIN, ICESCALENL, ICESCALEDS
690  REAL :: EMEAN, FMEAN, AMAX, CD, Z0, SCAT, &
691  SMOOTH_ICEDISP
692  REAL :: WN_R(NK), CG_ICE(NK), ALPHA_LIU(NK), ICECOEF2, R(NK)
693  DOUBLE PRECISION :: ATT, ISO
694  REAL :: EBAND, DIFF, EFINISH, HSTOT, PHINL, &
695  FMEAN1, FMEANWS, &
696  FACTOR, FACTOR2, DRAT, TAUWAX, TAUWAY, &
697  MWXFINISH, MWYFINISH, A1BAND, B1BAND, &
698  COSI(2)
699  REAL :: SPECINIT(NSPEC), SPEC2(NSPEC), FRLOCAL, JAC2
700  REAL :: DAM (NSPEC), DAM2(NSPEC), WN2(NSPEC), &
701  VSLN(NSPEC), &
702  VSIN(NSPEC), VDIN(NSPEC), &
703  VSNL(NSPEC), VDNL(NSPEC), &
704  VSDS(NSPEC), VDDS(NSPEC), &
705  VSBT(NSPEC), VDBT(NSPEC)
706  REAL :: VS(NSPEC), VD(NSPEC), EB(NK)
707 
708  LOGICAL :: SHAVE
709  LOGICAL :: LBREAK
710  LOGICAL, SAVE :: FIRST = .true.
711  LOGICAL :: PrintDeltaSmDA
712  REAL :: eInc1, eInc2, eVS, eVD, JAC
713  REAL :: DeltaSRC(NSPEC)
714 
715  REAL :: FOUT(NK,NTH), SOUT(NK,NTH), DOUT(NK,NTH)
716  REAL, SAVE :: TAUNUX, TAUNUY
717  LOGICAL, SAVE :: FLTEST = .false., flagnn = .true.
718 
719 #ifdef W3_OMPG
720  !$omp threadprivate( TAUNUX, TAUNUY)
721  !$omp threadprivate( FLTEST, FLAGNN )
722  !$omp threadprivate( FIRST )
723 #endif
724 
725  !/
726  !/ ------------------------------------------------------------------- /
727  !/ Local parameters dependent on compile switch
728  !/
729 #ifdef W3_S
730  INTEGER, SAVE :: IENT = 0
731 #endif
732 
733 #ifdef W3_NNT
734  INTEGER, SAVE :: NDSD = 89, ndsd2 = 88, j
735  REAL :: QCERR = 0. !/XNL2 and !/NNT
736 #endif
737 
738 #ifdef W3_NL5
739  INTEGER :: QI5TSTART(2)
740  REAL :: QR5KURT
741  INTEGER, PARAMETER :: NL5_SELECT = 1
742  REAL, PARAMETER :: NL5_OFFSET = 0. ! explicit dyn.
743 #endif
744 
745 #ifdef W3_SEED
746  REAL :: UC, SLEV
747 #endif
748 
749 #ifdef W3_MLIM
750  REAL :: HM, EM
751 #endif
752 
753 #ifdef W3_NNT
754  REAL :: FACNN
755 #endif
756 
757 #ifdef W3_T
758  REAL :: DTRAW
759 #endif
760 
761 #if defined(W3_IC1) || defined(W3_IC2) || defined(W3_IC3) || defined(W3_IC4) || defined(W3_IC5)
762  REAL :: VSIC(NSPEC), VDIC(NSPEC)
763 #endif
764 
765 #ifdef W3_DB1
766  REAL :: VSDB(NSPEC), VDDB(NSPEC)
767 #endif
768 
769 #ifdef W3_TR1
770  REAL :: VSTR(NSPEC), VDTR(NSPEC)
771 #endif
772 
773 #ifdef W3_BS1
774  REAL :: VSBS(NSPEC), VDBS(NSPEC)
775 #endif
776 
777 #ifdef W3_REF1
778  REAL :: VREF(NSPEC)
779 #endif
780 
781 #if defined(W3_IS1) || defined(W3_IS2)
782  REAL :: VSIR(NSPEC), VDIR(NSPEC)
783 #endif
784 
785 #ifdef W3_IS2
786  REAL :: VDIR2(NSPEC)
787  DOUBLE PRECISION :: SCATSPEC(NTH)
788 #endif
789 
790 #ifdef W3_UOST
791  REAL :: VSUO(NSPEC), VDUO(NSPEC)
792 #endif
793 
794 #ifdef W3_ST1
795  REAL :: FH1, FH2
796 #endif
797 
798 #ifdef W3_ST2
799  REAL :: FHTRAN, DFH, FACDIA, FACPAR
800 #endif
801 
802 #ifdef W3_ST3
803  REAL :: FMEANS, FH1, FH2
804 #endif
805 
806 #ifdef W3_ST4
807  REAL :: FMEANS, FH1, FH2, FAGE, DLWMEAN
808  REAL :: BRLAMBDA(NSPEC)
809 #endif
810 
811 #if defined(W3_ST3) || defined(W3_ST4)
812  LOGICAL :: LLWS(NSPEC)
813 #endif
814 
815 #ifdef W3_ST6
816  REAL :: VSWL(NSPEC), VDWL(NSPEC)
817 #endif
818 
819 #ifdef W3_PDLIB
820  REAL :: PreVS, DVS, SIDT, FAKS, MAXDAC
821 #endif
822 
823 #ifdef W3_NNT
824  CHARACTER(LEN=17), SAVE :: FNAME = 'test_data_nnn.ww3'
825 #endif
826  !
827  !/ -- End of variable delclarations
828  !
829 #ifdef W3_S
830  CALL strace (ient, 'W3SRCE')
831 #endif
832 
833 #ifdef W3_T
834  fltest = .true.
835 #endif
836  !
837  iks1 = 1
838 #ifdef W3_IG1
839  ! Does not integrate source terms for IG band if IGPARS(12) = 0.
840  IF (nint(igpars(12)).EQ.0) iks1 = nint(igpars(5))
841 #endif
842  is1=(iks1-1)*nth+1
843 
844  !! Initialise source term arrays:
845  vd = 0.
846  vs = 0.
847  vdio = 0.
848  vsio = 0.
849  vsbt = 0.
850  vdbt = 0.
851 
852 #if defined(W3_LN0) || defined(W3_LN1) || defined(W3_SEED)
853  vsln = 0.
854 #endif
855 
856 #if defined(W3_ST0) || defined(W3_ST3) || defined(W3_ST4)
857  vsin = 0.
858  vdin = 0.
859 #endif
860 
861 #if defined(W3_NL0) || defined(W3_NL1)
862  vsnl = 0.
863  vdnl = 0.
864 #endif
865 
866 #ifdef W3_TR1
867  vstr = 0.
868  vdtr = 0.
869 #endif
870 
871 #if defined(W3_ST0) || defined(W3_ST4)
872  vsds = 0.
873  vdds = 0.
874 #endif
875 
876 #ifdef W3_DB1
877  vsdb = 0.
878  vddb = 0.
879 #endif
880 
881 #if defined(W3_IC1) || defined(W3_IC2) || defined(W3_IC3) || defined(W3_IC4) || defined(W3_IC5)
882  vsic = 0.
883  vdic = 0.
884 #endif
885 
886 #ifdef W3_UOST
887  vsuo = 0.
888  vduo = 0.
889 #endif
890 
891 #if defined(W3_IS1) || defined(W3_IS2)
892  vsir = 0.
893  vdir = 0.
894 #endif
895 
896 #ifdef W3_IS2
897  vdir2 = 0.
898 #endif
899 
900 #ifdef W3_ST6
901  vswl = 0.
902  vdwl = 0.
903 #endif
904 
905 #if defined(W3_ST0) || defined(W3_ST1) || defined(W3_ST6)
906  zwnd = 10.
907 #endif
908 
909 #if defined(W3_ST2)
910  zwnd = zwind
911 #endif
912 
913 #if defined(W3_ST4)
914  zwnd = zzwnd
915 #endif
916  !
917  ! 1. Preparations --------------------------------------------------- *
918  !
919  depth = max( dmin , d_inp )
920  drat = dair / dwat
921  icescaleln = max(0.,min(1.,1.-ice*icescales(1)))
922  icescalein = max(0.,min(1.,1.-ice*icescales(2)))
923  icescalenl = max(0.,min(1.,1.-ice*icescales(3)))
924  icescaleds = max(0.,min(1.,1.-ice*icescales(4)))
925 
926 #ifdef W3_T
927  WRITE (ndst,9000)
928  WRITE (ndst,9001) depth, u10abs, u10dir*rade
929 #endif
930 
931  ! 1.a Set maximum change and wavenumber arrays.
932  !
933  !XP = 0.15
934  !FACP = XP / PI * 0.62E-3 * TPI**4 / GRAV**2
935  !
936  DO ik=1, nk
937  dam(1+(ik-1)*nth) = facp / ( sig(ik) * wn1(ik)**3 )
938  wn2(1+(ik-1)*nth) = wn1(ik)
939  END DO
940  !
941  DO ik=1, nk
942  is0 = (ik-1)*nth
943  DO ith=2, nth
944  dam(ith+is0) = dam(1+is0)
945  wn2(ith+is0) = wn2(1+is0)
946  END DO
947  END DO
948  !
949  ! 1.b Prepare dynamic time stepping
950  !
951  dtdyn = 0.
952  dttot = 0.
953  nsteps = 0
954  phiaw = 0.
955  charn = 0.
956  tws = 0.
957  phinl = 0.
958  phibbl = 0.
959  tauwix = 0.
960  tauwiy = 0.
961  tauwnx = 0.
962  tauwny = 0.
963  tauwax = 0.
964  tauway = 0.
965  tauscx = 0.
966  tauscy = 0.
967  taubbl = 0.
968  tauice = 0.
969  phice = 0.
970  tauocx = 0.
971  tauocy = 0.
972  wnmean = 0.
973 
974  !
975  ! TIME is updated in W3WAVEMD prior to the call of W3SCRE, we should
976  ! move 'TIME' one time step backward (QL)
977 #ifdef W3_NL5
978  qi5tstart = time
979  CALL tick21 (qi5tstart, -1.0 * dtg)
980 #endif
981  !
982 #ifdef W3_DEBUGSRC
983  IF (ix .eq. debug_node) THEN
984  WRITE(740+iaproc,*) 'W3SRCE start sum(SPEC)=', sum(spec)
985  WRITE(740+iaproc,*) 'W3SRCE start sum(SPECOLD)=', sum(specold)
986  WRITE(740+iaproc,*) 'W3SRCE start sum(SPECINIT)=', sum(specinit)
987  WRITE(740+iaproc,*) 'W3SRCE start sum(VSIO)=', sum(vsio)
988  WRITE(740+iaproc,*) 'W3SRCE start sum(VDIO)=', sum(vdio)
989  WRITE(740+iaproc,*) 'W3SRCE start USTAR=', ustar
990  END IF
991 #endif
992 
993 #ifdef W3_ST4
994  dlwmean= 0.
995  brlambda(:)=0.
996  whitecap(:)=0.
997 #endif
998  !
999  ! 1.c Set mean parameters
1000  !
1001 #ifdef W3_ST0
1002  CALL w3spr0 (spec, cg1, wn1, emean, fmean, wnmean, amax)
1003  fp = 0.85 * fmean
1004 #endif
1005 #ifdef W3_ST1
1006  CALL w3spr1 (spec, cg1, wn1, emean, fmean, wnmean, amax)
1007  fp = 0.85 * fmean
1008 #endif
1009 #ifdef W3_ST2
1010  CALL w3spr2 (spec, cg1, wn1, depth, fpi, u10abs, ustar, &
1011  emean, fmean, wnmean, amax, alpha, fp )
1012 #endif
1013 #ifdef W3_ST3
1014  tauwx=0.
1015  tauwy=0.
1016  IF ( it .eq. 0 ) THEN
1017  llws(:) = .true.
1018  ustar=0.
1019  ustdir=0.
1020  CALL w3spr3 (spec, cg1, wn1, emean, fmean, fmeans, wnmean, &
1021  amax, u10abs, u10dir, ustar, ustdir, &
1022  tauwx, tauwy, cd, z0, charn, llws, fmeanws)
1023  ELSE
1024  CALL w3spr3 (spec, cg1, wn1, emean, fmean, fmeans, wnmean, &
1025  amax, u10abs, u10dir, ustar, ustdir, &
1026  tauwx, tauwy, cd, z0, charn, llws, fmeanws)
1027  CALL w3sin3 ( spec, cg1, wn2, u10abs, ustar, drat, as, &
1028  u10dir, z0, cd, tauwx, tauwy, tauwax, tauway, &
1029  ice, vsin, vdin, llws, ix, iy )
1030  END IF
1031  CALL w3spr3 (spec, cg1, wn1, emean, fmean, fmeans, wnmean, &
1032  amax, u10abs, u10dir, ustar, ustdir, &
1033  tauwx, tauwy, cd, z0, charn, llws, fmeanws)
1034  tws = 1./fmeanws
1035 #endif
1036 #ifdef W3_ST4
1037  IF (sintailpar(4).GT.0.5) THEN ! this is designed to keep the bug as an option
1038  tauwx=0.
1039  tauwy=0.
1040  END IF
1041  IF ( it .EQ. 0 ) THEN
1042  llws(:) = .true.
1043  tauwx=0.
1044  tauwy=0.
1045  ustar=0.
1046  ustdir=0.
1047  ELSE
1048  CALL w3spr4 (spec, cg1, wn1, emean, fmean, fmean1, wnmean, &
1049  amax, u10abs, u10dir, &
1050 #ifdef W3_FLX5
1051  taua, tauadir, dair, &
1052 #endif
1053  ustar, ustdir, &
1054  tauwx, tauwy, cd, z0, charn, llws, fmeanws, dlwmean)
1055 #endif
1056 
1057 #if defined(W3_DEBUGSRC) && defined(W3_ST4)
1058  IF (ix == debug_node) THEN
1059  WRITE(740+iaproc,*) '1: out value USTAR=', ustar, ' USTDIR=', ustdir
1060  WRITE(740+iaproc,*) '1: out value EMEAN=', emean, ' FMEAN=', fmean
1061  WRITE(740+iaproc,*) '1: out value FMEAN1=', fmean1, ' WNMEAN=', wnmean
1062  WRITE(740+iaproc,*) '1: out value CD=', cd, ' Z0=', z0
1063  WRITE(740+iaproc,*) '1: out value ALPHA=', charn, ' FMEANWS=', fmeanws
1064  END IF
1065 #endif
1066 
1067 #ifdef W3_ST4
1068  IF (sintailpar(4).GT.0.5) CALL w3sin4 ( spec, cg1, wn2, u10abs, ustar, drat, as, &
1069  u10dir, z0, cd, tauwx, tauwy, tauwax, tauway, &
1070  vsin, vdin, llws, ix, iy, brlambda )
1071  END IF
1072 #endif
1073 #if defined(W3_DEBUGSRC) && defined(W3_ST4)
1074  IF (ix == debug_node) THEN
1075  WRITE(740+iaproc,*) '1: U10DIR=', u10dir, ' Z0=', z0, ' CHARN=', charn
1076  WRITE(740+iaproc,*) '1: USTAR=', ustar, ' U10ABS=', u10abs, ' AS=', as
1077  WRITE(740+iaproc,*) '1: DRAT=', drat
1078  WRITE(740+iaproc,*) '1: TAUWX=', tauwx, ' TAUWY=', tauwy
1079  WRITE(740+iaproc,*) '1: TAUWAX=', tauwax, ' TAUWAY=', tauway
1080  WRITE(740+iaproc,*) '1: min(CG1)=', minval(cg1), ' max(CG1)=', maxval(cg1)
1081  WRITE(740+iaproc,*) '1: W3SIN4(min/max/sum)VSIN=', minval(vsin), maxval(vsin), sum(vsin)
1082  WRITE(740+iaproc,*) '1: W3SIN4(min/max/sum)VDIN=', minval(vdin), maxval(vdin), sum(vdin)
1083  END IF
1084 #endif
1085 
1086 #ifdef W3_ST4
1087  CALL w3spr4 (spec, cg1, wn1, emean, fmean, fmean1, wnmean, &
1088  amax, u10abs, u10dir, &
1089 #ifdef W3_FLX5
1090  taua, tauadir, dair, &
1091 #endif
1092  ustar, ustdir, &
1093  tauwx, tauwy, cd, z0, charn, llws, fmeanws, dlwmean)
1094  tws = 1./fmeanws
1095 #endif
1096 #ifdef W3_ST6
1097  CALL w3spr6 (spec, cg1, wn1, emean, fmean, wnmean, amax, fp)
1098 #endif
1099  !
1100  ! 1.c2 Stores the initial data
1101  !
1102  specinit = spec
1103  !
1104  ! 1.d Stresses
1105  !
1106 #ifdef W3_FLX1
1107  CALL w3flx1 ( zwnd, u10abs, u10dir, ustar, ustdir, z0, cd )
1108 #endif
1109 #ifdef W3_FLX2
1110  CALL w3flx2 ( zwnd, depth, fp, u10abs, u10dir, &
1111  ustar, ustdir, z0, cd )
1112 #endif
1113 #ifdef W3_FLX3
1114  CALL w3flx3 ( zwnd, depth, fp, u10abs, u10dir, &
1115  ustar, ustdir, z0, cd )
1116 #endif
1117 #ifdef W3_FLX4
1118  CALL w3flx4 ( zwnd, u10abs, u10dir, ustar, ustdir, z0, cd )
1119 #endif
1120 #ifdef W3_FLX5
1121  CALL w3flx5 ( zwnd, u10abs, u10dir, taua, tauadir, dair, &
1122  ustar, ustdir, z0, cd, charn )
1123 #endif
1124  !
1125  ! 1.e Prepare cut-off beyond which the tail is imposed with a power law
1126  !
1127 #ifdef W3_ST0
1128  fhigh = sig(nk)
1129 #endif
1130 #ifdef W3_ST1
1131  fh1 = fxfm * fmean
1132  fh2 = fxpm / ustar
1133  fhigh = max( fh1 , fh2 )
1134  IF (fltest) WRITE (ndst,9004) fh1*tpiinv, fh2*tpiinv, fhigh*tpiinv
1135 #endif
1136 #ifdef W3_ST2
1137  fhigh = xfc * fpi
1138 #endif
1139 #ifdef W3_ST3
1140  fhigh = max(ffxfm * max(fmean,fmeanws),ffxpm / ustar)
1141 #endif
1142 #ifdef W3_ST4
1143  ! Introduces a Long & Resio (JGR2007) type dependance on wave age
1144  ! !/ST4 FAGE = FFXFA*TANH(0.3*U10ABS*FMEANWS*TPI/GRAV)
1145  fage = 0.
1146  fhigh = max( (ffxfm + fage ) * max(fmean1,fmeanws), ffxpm / ustar)
1147  fhigi = ffxfa * fmean1
1148 #endif
1149 #ifdef W3_ST6
1150  IF (fxfm .LE. 0) THEN
1151  fhigh = sig(nk)
1152  ELSE
1153  fhigh = max(fxfm * fmean, fxpm / ustar)
1154  ENDIF
1155 #endif
1156  !
1157  ! 1.f Prepare output file for !/NNT option
1158  !
1159 #ifdef W3_NNT
1160  IF ( it .EQ. 0 ) THEN
1161  j = len_trim(fnmpre)
1162  WRITE (fname(11:13),'(I3.3)') iaproc
1163  OPEN (ndsd,file=fnmpre(:j)//fname,form='UNFORMATTED', convert=file_endian, &
1164  err=800,iostat=ierr)
1165  WRITE (ndsd,err=801,iostat=ierr) nk, nth
1166  WRITE (ndsd,err=801,iostat=ierr) sig(1:nk) * tpiinv
1167  OPEN (ndsd2,file=fnmpre(:j)//'time.ww3', &
1168  form='FORMATTED',err=800,iostat=ierr)
1169  END IF
1170 #endif
1171  !
1172  ! ... Branch point dynamic integration - - - - - - - - - - - - - - - -
1173  !
1174  DO
1175  !
1176  nsteps = nsteps + 1
1177  !
1178 #ifdef W3_T
1179  WRITE (ndst,9020) nsteps, dttot
1180 #endif
1181  !
1182  ! 2. Calculate source terms ----------------------------------------- *
1183  !
1184  ! 2.a Input.
1185  !
1186 #ifdef W3_LN1
1187  CALL w3sln1 ( wn1, fhigh, ustar, u10dir , vsln )
1188 #endif
1189  !
1190 #ifdef W3_ST1
1191  CALL w3sin1 ( spec, wn2, ustar, u10dir , vsin, vdin )
1192 #endif
1193 #ifdef W3_ST2
1194  CALL w3sin2 ( spec, cg1, wn2, u10abs, u10dir, cd, z0, &
1195  fpi, vsin, vdin )
1196 #endif
1197 #ifdef W3_ST3
1198  CALL w3sin3 ( spec, cg1, wn2, u10abs, ustar, drat, as, &
1199  u10dir, z0, cd, tauwx, tauwy, tauwax, tauway, &
1200  ice, vsin, vdin, llws, ix, iy )
1201 #endif
1202 #ifdef W3_ST4
1203  CALL w3sin4 ( spec, cg1, wn2, u10abs, ustar, drat, as, &
1204  u10dir, z0, cd, tauwx, tauwy, tauwax, tauway, &
1205  vsin, vdin, llws, ix, iy, brlambda )
1206 #endif
1207 
1208 #if defined(W3_DEBUGSRC) && defined(W3_ST4)
1209  IF (ix == debug_node) THEN
1210  WRITE(740+iaproc,*) '2 : W3SIN4(min/max/sum)VSIN=', minval(vsin), maxval(vsin), sum(vsin)
1211  WRITE(740+iaproc,*) '2 : W3SIN4(min/max/sum)VDIN=', minval(vdin), maxval(vdin), sum(vdin)
1212  END IF
1213 #endif
1214 
1215 #ifdef W3_ST6
1216  CALL w3sin6 ( spec, cg1, wn2, u10abs, ustar, ustdir, cd, dair, &
1217  tauwx, tauwy, tauwax, tauway, vsin, vdin )
1218 #endif
1219  !
1220  ! 2.b Nonlinear interactions.
1221  !
1222 #ifdef W3_NL1
1223  IF (iqtpe.GT.0) THEN
1224  CALL w3snl1 ( spec, cg1, wnmean*depth, vsnl, vdnl )
1225  ELSE
1226  CALL w3snlgqm ( spec, cg1, wn1, depth, vsnl, vdnl )
1227  END IF
1228 #endif
1229 #ifdef W3_NL2
1230  CALL w3snl2 ( spec, cg1, depth, vsnl, vdnl )
1231 #endif
1232 #ifdef W3_NL3
1233  CALL w3snl3 ( spec, cg1, wn1, depth, vsnl, vdnl )
1234 #endif
1235 #ifdef W3_NL4
1236  CALL w3snl4 ( spec, cg1, wn1, depth, vsnl, vdnl )
1237 #endif
1238 #ifdef W3_NL5
1239  CALL w3snl5 ( spec, cg1, wn1, fmean, qi5tstart, &
1240  u10abs, u10dir, jsea, vsnl, vdnl, qr5kurt)
1241 #endif
1242  !
1243 #ifdef W3_PDLIB
1244  IF (.NOT. fssource .or. lsloc) THEN
1245 #endif
1246 #ifdef W3_TR1
1247  CALL w3str1 ( spec, specold, cg1, wn1, depth, ix, vstr, vdtr )
1248 #endif
1249 #ifdef W3_PDLIB
1250  ENDIF
1251 #endif
1252  !
1253  ! 2.c Dissipation... except for ST4
1254  ! 2.c1 as in source term package
1255  !
1256 #ifdef W3_ST1
1257  CALL w3sds1 ( spec, wn2, emean, fmean, wnmean, vsds, vdds )
1258 #endif
1259 #ifdef W3_ST2
1260  CALL w3sds2 ( spec, cg1, wn1, fpi, ustar, alpha,vsds, vdds )
1261 #endif
1262 #ifdef W3_ST3
1263  CALL w3sds3 ( spec, wn1, cg1, emean, fmeans, wnmean, &
1264  ustar, ustdir, depth, vsds, vdds, ix, iy )
1265 #endif
1266 #ifdef W3_ST4
1267  CALL w3sds4 ( spec, wn1, cg1, ustar, ustdir, depth, dair, vsds, &
1268  vdds, ix, iy, brlambda, whitecap, dlwmean )
1269 #endif
1270 #if defined(W3_DEBUGSRC) && defined(W3_ST4)
1271  IF (ix == debug_node) THEN
1272  WRITE(740+iaproc,*) '2 : W3SDS4(min/max/sum)VSDS=', minval(vsds), maxval(vsds), sum(vsds)
1273  WRITE(740+iaproc,*) '2 : W3SDS4(min/max/sum)VDDS=', minval(vdds), maxval(vdds), sum(vdds)
1274  END IF
1275 #endif
1276 
1277 #ifdef W3_ST6
1278  CALL w3sds6 ( spec, cg1, wn1, vsds, vdds )
1279 #endif
1280  !
1281 #ifdef W3_PDLIB
1282  IF (.NOT. fssource .or. lsloc) THEN
1283 #endif
1284 #ifdef W3_DB1
1285  CALL w3sdb1 ( ix, spec, depth, emean, fmean, wnmean, cg1, &
1286  lbreak, vsdb, vddb )
1287 #endif
1288 #ifdef W3_PDLIB
1289  ENDIF
1290 #endif
1291  !
1292  ! 2.c2 optional dissipation parameterisations
1293  !
1294 #ifdef W3_ST6
1295  IF (swl6s6) THEN
1296  CALL w3swl6 ( spec, cg1, wn1, vswl, vdwl )
1297  END IF
1298 #endif
1299  !
1300  ! 2.d Bottom interactions.
1301  !
1302 #ifdef W3_BT1
1303  CALL w3sbt1 ( spec, cg1, wn1, depth, vsbt, vdbt )
1304 #endif
1305 #ifdef W3_BT4
1306  CALL w3sbt4 ( spec, cg1, wn1, depth, d50, psic, taubbl, &
1307  bedform, vsbt, vdbt, ix, iy )
1308 #endif
1309 #ifdef W3_BT8
1310  CALL w3sbt8 ( spec, depth, vsbt, vdbt, ix, iy )
1311 #endif
1312 #ifdef W3_BT9
1313  CALL w3sbt9 ( spec, depth, vsbt, vdbt, ix, iy )
1314 #endif
1315  !
1316 #ifdef W3_BS1
1317  CALL w3sbs1 ( spec, cg1, wn1, depth, cx, cy, &
1318  tauscx, tauscy, vsbs, vdbs )
1319 #endif
1320  !
1321  ! 2.e Unresolved Obstacles Source Term
1322  !
1323 #ifdef W3_UOST
1324  ! UNRESOLVED OBSTACLES
1325  CALL uost_srctrmcompute(ix, iy, spec, cg1, dt, &
1326  u10abs, u10dir, vsuo, vduo)
1327 #endif
1328  !
1329  ! 2.g Dump training data if necessary
1330  !
1331 #ifdef W3_NNT
1332  WRITE (screen,8888) time, dttot, flagnn, qcerr
1333  WRITE (ndsd2,8888) time, dttot, flagnn, qcerr
1334 8888 FORMAT (1x,i8.8,1x,i6.6,f8.1,l2,f8.2)
1335  WRITE (ndsd,err=801,iostat=ierr) ix, iy, time, nsteps, &
1336  dttot, flagnn, depth, u10abs, u10dir
1337  !
1338  IF ( flagnn ) THEN
1339  DO ik=1, nk
1340  facnn = tpi * sig(ik) / cg1(ik)
1341  DO ith=1, nth
1342  is = ith + (ik-1)*nth
1343  fout(ik,ith) = spec(is) * facnn
1344  sout(ik,ith) = vsnl(is) * facnn
1345  dout(ik,ith) = vdnl(is)
1346  END DO
1347  END DO
1348  WRITE (ndsd,err=801,iostat=ierr) fout
1349  WRITE (ndsd,err=801,iostat=ierr) sout
1350  WRITE (ndsd,err=801,iostat=ierr) dout
1351  END IF
1352 #endif
1353  !
1354  ! 3. Set frequency cut-off ------------------------------------------ *
1355  !
1356 #ifdef W3_ST2
1357  fhigh = xfc * fpi
1358  IF ( fltest ) WRITE (ndst,9005) fhigh*tpiinv
1359 #endif
1360  nkh = min( nk , int(facti2+facti1*log(max(1.e-7,fhigh))) )
1361  nkh1 = min( nk , nkh+1 )
1362  nspech = nkh1*nth
1363 #ifdef W3_T
1364  WRITE (ndst,9021) nkh, nkh1, nspech
1365 #endif
1366  !
1367  ! 4. Summation of source terms and diagonal term and time step ------ *
1368  !
1369  dt = min( dtg-dttot , dtmax )
1370  afilt = max( dam(nspec) , xflt*amax )
1371  !
1372  ! For input and dissipation calculate the fraction of the ice-free
1373  ! surface. In the presence of ice, the effective water surface
1374  ! is reduce to a fraction of the cell size free from ice, and so is
1375  ! input :
1376  ! SIN = (1-ICE)**ISCALEIN*SIN and SDS=(1-ICE)**ISCALEDS*SDS ------------------ *
1377  ! INFLAGS2(4) is true if ice concentration was ever read during
1378  ! this simulation
1379  IF ( inflags2(4) ) THEN
1380  vsnl(1:nspech) = icescalenl * vsnl(1:nspech)
1381  vdnl(1:nspech) = icescalenl * vdnl(1:nspech)
1382  vsln(1:nspech) = icescaleln * vsln(1:nspech)
1383  vsin(1:nspech) = icescalein * vsin(1:nspech)
1384  vdin(1:nspech) = icescalein * vdin(1:nspech)
1385  vsds(1:nspech) = icescaleds * vsds(1:nspech)
1386  vdds(1:nspech) = icescaleds * vdds(1:nspech)
1387  END IF
1388 
1389 #ifdef W3_PDLIB
1390  IF (b_jgs_limiter_func == 2) THEN
1391  DO ik=1, nk
1392  jac = cg1(ik)/clatsl
1393  jac2 = 1./tpi/sig(ik)
1394  frlocal = sig(ik)*tpiinv
1395 #ifdef W3_ST6
1396  dam2(1+(ik-1)*nth) = 5e-7 * grav/frlocal**4 * ustar * fmean * dtmin * jac * jac2
1397 #else
1398  dam2(1+(ik-1)*nth) = 5e-7 * grav/frlocal**4 * ustar * max(fmeanws,fmean) * dtmin * jac * jac2
1399 #endif
1400  !FROM WWM: 5E-7 * GRAV/FR(IS)**4 * USTAR * MAX(FMEANWS(IP),FMEAN(IP)) * DT4S/PI2/SPSIG(IS)
1401  END DO
1402  DO ik=1, nk
1403  is0 = (ik-1)*nth
1404  DO ith=2, nth
1405  dam2(ith+is0) = dam2(1+is0)
1406  END DO
1407  END DO
1408  ENDIF
1409 #endif
1410  !
1411  DO is=is1, nspech
1412  vs(is) = vsln(is) + vsin(is) + vsnl(is) &
1413  + vsds(is) + vsbt(is)
1414 #ifdef W3_ST6
1415  vs(is) = vs(is) + vswl(is)
1416 #endif
1417 #if defined(W3_TR1) && !defined(W3_PDLIB)
1418  vs(is) = vs(is) + vstr(is)
1419 #endif
1420 #ifdef W3_BS1
1421  vs(is) = vs(is) + vsbs(is)
1422 #endif
1423 #ifdef W3_UOST
1424  vs(is) = vs(is) + vsuo(is)
1425 #endif
1426  vd(is) = vdin(is) + vdnl(is) &
1427  + vdds(is) + vdbt(is)
1428 #ifdef W3_ST6
1429  vd(is) = vd(is) + vdwl(is)
1430 #endif
1431 #if defined(W3_TR1) && !defined(W3_PDLIB)
1432  vd(is) = vd(is) + vdtr(is)
1433 #endif
1434 #ifdef W3_BS1
1435  vd(is) = vd(is) + vdbs(is)
1436 #endif
1437 #ifdef W3_UOST
1438  vd(is) = vd(is) + vduo(is)
1439 #endif
1440  damax = min( dam(is) , max( xrel*specinit(is) , afilt ) )
1441  afac = 1. / max( 1.e-10 , abs(vs(is)/damax) )
1442 #ifdef W3_NL5
1443  IF (nl5_select .EQ. 1) THEN
1444  dt = min( dt , afac / ( max( 1.e-10, &
1445  1. + nl5_offset*afac*min(0.,vd(is)) ) ) )
1446  ELSE
1447 #endif
1448  dt = min( dt , afac / ( max( 1.e-10, &
1449  1. + offset*afac*min(0.,vd(is)) ) ) )
1450 #ifdef W3_NL5
1451  ENDIF
1452 #endif
1453  END DO ! end of loop on IS
1454 
1455  !
1456  dt = max( 0.5, dt ) ! The hardcoded min. dt is a problem for certain cases e.g. laborotary scale problems.
1457  !
1458  dtdyn = dtdyn + dt
1459 #ifdef W3_T
1460  dtraw = dt
1461 #endif
1462  idt = 1 + int( 0.99*(dtg-dttot)/dt ) ! number of iterations
1463  dt = (dtg-dttot)/real(idt) ! actualy time step
1464  shave = dt.LT.dtmin .AND. dt.LT.dtg-dttot ! limiter check ...
1465  shaveio = shave
1466  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!
1467  !
1468 #ifdef W3_NL5
1469  dt = int(dt) * 1.0
1470 #endif
1471  IF (srce_call .eq. srce_imp_post) dt = dtg ! for implicit part
1472 #ifdef W3_NL5
1473  IF (nl5_select .EQ. 1) THEN
1474  hdt = nl5_offset * dt
1475  ELSE
1476 #endif
1477  hdt = offset * dt
1478 #ifdef W3_NL5
1479  ENDIF
1480 #endif
1481  dttot = dttot + dt
1482 
1483 #ifdef W3_DEBUGSRC
1484  IF (ix == debug_node) WRITE(*,'(A20,2I10,5F20.10,L20)') 'TIMINGS 2', idt, nsteps, dt, dtmin, dtdyn, hdt, dttot, shave
1485  IF (ix == debug_node) THEN
1486  WRITE(740+iaproc,*) '1: min/max/sum(VS)=', minval(vs), maxval(vs), sum(vs)
1487  WRITE(740+iaproc,*) '1: min/max/sum(VD)=', minval(vd), maxval(vd), sum(vd)
1488  WRITE(740+iaproc,*) 'min/max/sum(VSIN)=', minval(vsin), maxval(vsin), sum(vsin)
1489  WRITE(740+iaproc,*) 'min/max/sum(VDIN)=', minval(vdin), maxval(vdin), sum(vdin)
1490  WRITE(740+iaproc,*) 'min/max/sum(VSLN)=', minval(vsln), maxval(vsln), sum(vsln)
1491  WRITE(740+iaproc,*) 'min/max/sum(VSNL)=', minval(vsnl), maxval(vsnl), sum(vsnl)
1492  WRITE(740+iaproc,*) 'min/max/sum(VDNL)=', minval(vdnl), maxval(vdnl), sum(vdnl)
1493  WRITE(740+iaproc,*) 'min/max/sum(VSDS)=', minval(vsds), maxval(vsds), sum(vsds)
1494  WRITE(740+iaproc,*) 'min/max/sum(VDDS)=', minval(vdds), maxval(vdds), sum(vdds)
1495 #ifdef W3_ST6
1496  WRITE(740+iaproc,*) 'min/max/sum(VSWL)=', minval(vswl), maxval(vswl), sum(vswl)
1497  WRITE(740+iaproc,*) 'min/max/sum(VDWL)=', minval(vdwl), maxval(vdwl), sum(vdwl)
1498 #endif
1499 #ifdef W3_DB1
1500  WRITE(740+iaproc,*) 'min/max/sum(VSDB)=', minval(vsdb), maxval(vsdb), sum(vsdb)
1501  WRITE(740+iaproc,*) 'min/max/sum(VDDB)=', minval(vddb), maxval(vddb), sum(vddb)
1502 #endif
1503 #ifdef W3_TR1
1504  WRITE(740+iaproc,*) 'min/max/sum(VSTR)=', minval(vstr), maxval(vstr), sum(vstr)
1505  WRITE(740+iaproc,*) 'min/max/sum(VDTR)=', minval(vdtr), maxval(vdtr), sum(vdtr)
1506 #endif
1507 #ifdef W3_BS1
1508  WRITE(740+iaproc,*) 'min/max/sum(VSBS)=', minval(vsbs), maxval(vsbs), sum(vsbs)
1509  WRITE(740+iaproc,*) 'min/max/sum(VDBS)=', minval(vdbs), maxval(vdbs), sum(vdbs)
1510 #endif
1511  WRITE(740+iaproc,*) 'min/max/sum(VSBT)=', minval(vsbt), maxval(vsbt), sum(vsbt)
1512  WRITE(740+iaproc,*) 'min/max/sum(VDBT)=', minval(vdbt), maxval(vdbt), sum(vdbt)
1513  END IF
1514 #endif
1515 
1516 #ifdef W3_PDLIB
1517  IF (srce_call .eq. srce_imp_pre) THEN
1518  IF (lsloc) THEN
1519  IF (imem == 1) THEN
1520  sidt = pdlib_si(jsea) * dtg
1521  DO ik = 1, nk
1522  jac = clatsl/cg1(ik)
1523  DO ith = 1, nth
1524  isp = ith + (ik-1)*nth
1525  vd(isp) = min(0., vd(isp))
1526  IF (b_jgs_limiter_func == 2) THEN
1527  maxdac = max(dam(isp),dam2(isp))
1528  ELSE
1529  maxdac = dam(isp)
1530  ENDIF
1531  faks = dtg / max( 1. , (1.-dtg*vd(isp)))
1532  dvs = vs(isp) * faks
1533  IF (.NOT. b_jgs_limiter) THEN
1534  dvs = sign(min(maxdac,abs(dvs)),dvs)
1535  ENDIF
1536  prevs = dvs / faks
1537  evs = prevs / cg1(ik) * clatsl
1538  evd = min(0.,vd(isp))
1539  b_jac(isp,jsea) = b_jac(isp,jsea) + sidt * (evs - evd*spec(isp)*jac)
1540  aspar_jac(isp,pdlib_i_diag(jsea)) = aspar_jac(isp,pdlib_i_diag(jsea)) - sidt * evd
1541 #ifdef W3_DB1
1542  evs = vsdb(isp) * jac
1543  evd = min(0.,vddb(isp))
1544  IF (evs .gt. 0.) THEN
1545  evs = 2*evs
1546  evd = -evd
1547  ELSE
1548  evs = -evs
1549  evd = 2*evd
1550  ENDIF
1551 #endif
1552  b_jac(isp,jsea) = b_jac(isp,jsea) + sidt * evs
1553  aspar_jac(isp,pdlib_i_diag(jsea)) = aspar_jac(isp,pdlib_i_diag(jsea)) - sidt * evd
1554 
1555 #ifdef W3_TR1
1556  evs = vstr(isp) * jac
1557  evd = vdtr(isp)
1558  IF (evs .gt. 0.) THEN
1559  evs = 2*evs
1560  evd = -evd
1561  ELSE
1562  evs = -evs
1563  evd = 2*evd
1564  ENDIF
1565 #endif
1566  b_jac(isp,jsea) = b_jac(isp,jsea) + sidt * evs
1567  aspar_jac(isp,pdlib_i_diag(jsea)) = aspar_jac(isp,pdlib_i_diag(jsea)) - sidt * evd
1568  END DO
1569  END DO
1570 
1571  ELSEIF (imem == 2) THEN
1572 
1573  sidt = pdlib_si(jsea) * dtg
1574  DO ik=1,nk
1575  jac = clatsl/cg1(ik)
1576  DO ith=1,nth
1577  isp=ith + (ik-1)*nth
1578  vd(isp) = min(0., vd(isp))
1579  IF (b_jgs_limiter_func == 2) THEN
1580  maxdac = max(dam(isp),dam2(isp))
1581  ELSE
1582  maxdac = dam(isp)
1583  ENDIF
1584  faks = dtg / max( 1. , (1.-dtg*vd(isp)))
1585  dvs = vs(isp) * faks
1586  IF (.NOT. b_jgs_limiter) THEN
1587  dvs = sign(min(maxdac,abs(dvs)),dvs)
1588  ENDIF
1589  prevs = dvs / faks
1590  evs = prevs / cg1(ik) * clatsl
1591  evd = vd(isp)
1592 #ifdef W3_DB1
1593  evs = evs + dble(vsdb(isp)) * jac
1594  evd = evd + min(0.,dble(vddb(isp)))
1595  b_jac(isp,jsea) = b_jac(isp,jsea) + sidt * (evs - evd*va(isp,jsea))
1596  aspar_diag_all(isp,jsea) = aspar_diag_all(isp,jsea) - sidt * evd
1597 #endif
1598  END DO
1599  END DO
1600  ENDIF
1601  ENDIF
1602 
1603  printdeltasmda=.false.
1604  IF (printdeltasmda .eqv. .true.) THEN
1605  DO is=1,nspec
1606  deltasrc(is) = vsin(is) - spec(is)*vdin(is)
1607  END DO
1608  WRITE(740+iaproc,*) 'min/max/sum(VSIN)=', minval(vsin), maxval(vsin), sum(vsin)
1609  WRITE(740+iaproc,*) 'min/max/sum(DeltaIN)=', minval(deltasrc), maxval(deltasrc), sum(deltasrc)
1610  !
1611  DO is=1,nspec
1612  deltasrc(is) = vsnl(is) - spec(is)*vdnl(is)
1613  END DO
1614  WRITE(740+iaproc,*) 'min/max/sum(VSNL)=', minval(vsnl), maxval(vsnl), sum(vsnl)
1615  WRITE(740+iaproc,*) 'min/max/sum(DeltaNL)=', minval(deltasrc), maxval(deltasrc), sum(deltasrc)
1616  !
1617  DO is=1,nspec
1618  deltasrc(is) = vsds(is) - spec(is)*vdds(is)
1619  END DO
1620  WRITE(740+iaproc,*) 'min/max/sum(VSDS)=', minval(vsds), maxval(vsds), sum(vsds)
1621  WRITE(740+iaproc,*) 'min/max/sum(DeltaDS)=', minval(deltasrc), maxval(deltasrc), sum(deltasrc)
1622  !
1623  ! DO IS=1,NSPEC
1624  ! DeltaSRC(IS) = VSIC(IS) - SPEC(IS)*VDIC(IS)
1625  ! END DO
1626  WRITE(740+iaproc,*) 'min/max/sum(DeltaDS)=', minval(deltasrc), maxval(deltasrc), sum(deltasrc)
1627  END IF
1628 
1629  IF (.not. lsloc) THEN
1630  IF (optioncall .eq. 1) THEN
1631  CALL sign_vsd_patankar_ww3(spec,vs,vd)
1632  ELSE IF (optioncall .eq. 2) THEN
1633  CALL sign_vsd_semi_implicit_ww3(spec,vs,vd)
1634  ELSE IF (optioncall .eq. 3) THEN
1635  CALL sign_vsd_semi_implicit_ww3(spec,vs,vd)
1636  ENDIF
1637  vsio = vs
1638  vdio = vd
1639  ENDIF
1640 
1641 #ifdef W3_DEBUGSRC
1642  IF (ix == debug_node) THEN
1643  WRITE(740+iaproc,*) ' srce_imp_pre : SHAVE = ', shave
1644  WRITE(740+iaproc,*) ' srce_imp_pre : DT=', dt, ' HDT=', hdt, 'DTG=', dtg
1645  WRITE(740+iaproc,*) ' srce_imp_pre : sum(SPEC)=', sum(spec)
1646  WRITE(740+iaproc,*) ' srce_imp_pre : sum(VSTOT)=', sum(vs)
1647  WRITE(740+iaproc,*) ' srce_imp_pre : sum(VDTOT)=', sum(min(0. , vd))
1648  END IF
1649 
1650  IF (ix == debug_node) WRITE(44,'(1EN15.4)') sum(vsin)
1651  IF (ix == debug_node) WRITE(44,'(1EN15.4)') sum(vdin)
1652  IF (ix == debug_node) WRITE(44,'(1EN15.4)') sum(vsds)
1653  IF (ix == debug_node) WRITE(44,'(1EN15.4)') sum(vdds)
1654  IF (ix == debug_node) WRITE(44,'(1EN15.4)') sum(vsnl)
1655  IF (ix == debug_node) WRITE(44,'(1EN15.4)') sum(vdnl)
1656  IF (ix == debug_node) WRITE(44,'(1EN15.4)') sum(vsln)
1657  IF (ix == debug_node) WRITE(44,'(1EN15.4)') sum(vsbt)
1658  IF (ix == debug_node) WRITE(44,'(1EN15.4)') sum(vs)
1659  IF (ix == debug_node) WRITE(44,'(1EN15.4)') sum(vd)
1660 #endif
1661  RETURN ! return everything is done for the implicit ...
1662 
1663  END IF ! srce_imp_pre
1664 ! --end W3_PDLIB
1665 #endif
1666  !
1667 #ifdef W3_T
1668  WRITE (ndst,9040) dtraw, dt, shave
1669 #endif
1670  !
1671  ! 5. Increment spectrum --------------------------------------------- *
1672  !
1673  IF (srce_call .eq. srce_direct) THEN
1674  IF ( shave ) THEN
1675  DO is=is1, nspech
1676  einc1 = vs(is) * dt / max( 1. , (1.-hdt*vd(is)))
1677  einc2 = sign( min(dam(is),abs(einc1)) , einc1 )
1678  spec(is) = max( 0. , spec(is)+einc2 )
1679  END DO
1680  ELSE
1681  !
1682  DO is=is1, nspech
1683  einc1 = vs(is) * dt / max( 1. , (1.-hdt*vd(is)))
1684  spec(is) = max( 0. , spec(is)+einc1 )
1685  END DO
1686  END IF
1687  !
1688 #ifdef W3_DB1
1689  DO is=is1, nspech
1690  einc1 = vsdb(is) * dt / max( 1. , (1.-hdt*vddb(is)))
1691  spec(is) = max( 0. , spec(is)+einc1 )
1692  END DO
1693 #endif
1694 #ifdef W3_TR1
1695  DO is=is1, nspech
1696  einc1 = vdtr(is) * dt / max( 1. , (1.-hdt*vdtr(is)))
1697  spec(is) = max( 0. , spec(is)+einc1 )
1698  END DO
1699 #endif
1700 
1701 #ifdef W3_DEBUGSRC
1702  IF (ix == debug_node) WRITE(44,'(1EN15.4)') sum(vsin)
1703  IF (ix == debug_node) WRITE(44,'(1EN15.4)') sum(vdin)
1704  IF (ix == debug_node) WRITE(44,'(1EN15.4)') sum(vsds)
1705  IF (ix == debug_node) WRITE(44,'(1EN15.4)') sum(vdds)
1706  IF (ix == debug_node) WRITE(44,'(1EN15.4)') sum(vsnl)
1707  IF (ix == debug_node) WRITE(44,'(1EN15.4)') sum(vdnl)
1708  IF (ix == debug_node) WRITE(44,'(1EN15.4)') sum(vsln)
1709  IF (ix == debug_node) WRITE(44,'(1EN15.4)') sum(vsbt)
1710  IF (ix == debug_node) WRITE(44,'(1EN15.4)') sum(vs)
1711  IF (ix == debug_node) WRITE(44,'(1EN15.4)') sum(vd)
1712  IF (ix == debug_node) THEN
1713  WRITE(740+iaproc,*) ' srce_direct : SHAVE = ', shave
1714  WRITE(740+iaproc,*) ' srce_direct : DT=', dt, ' HDT=', hdt, 'DTG=', dtg
1715  WRITE(740+iaproc,*) ' srce_direct : sum(SPEC)=', sum(spec)
1716  WRITE(740+iaproc,*) ' srce_direct : sum(VSTOT)=', sum(vs)
1717  WRITE(740+iaproc,*) ' srce_direct : sum(VDTOT)=', sum(min(0. , vd))
1718  END IF
1719 #endif
1720  END IF ! srce_call .eq. srce_direct
1721  !
1722  ! 5.b Computes
1723  ! atmos->wave flux PHIAW-------------------------------- *
1724  ! wave ->BBL flux PHIBBL------------------------------- *
1725  ! wave ->ice flux PHICE ------------------------------- *
1726  !
1727  whitecap(3)=0.
1728  hstot=0.
1729  DO ik=iks1, nk
1730  factor = dden(ik)/cg1(ik) !Jacobian to get energy in band
1731  factor2= factor*grav*wn1(ik)/sig(ik) ! coefficient to get momentum
1732 
1733  ! Wave direction is "direction to"
1734  ! therefore there is a PLUS sign for the stress
1735  DO ith=1, nth
1736  is = (ik-1)*nth + ith
1737  cosi(1)=ecos(is)
1738  cosi(2)=esin(is)
1739  phiaw = phiaw + (vsin(is))* dt * factor &
1740  / max( 1. , (1.-hdt*vdin(is))) ! semi-implict integration scheme
1741 
1742  phibbl= phibbl- (vsbt(is))* dt * factor &
1743  / max( 1. , (1.-hdt*vdbt(is))) ! semi-implict integration scheme
1744  phinl = phinl + vsnl(is)* dt * factor &
1745  / max( 1. , (1.-hdt*vdnl(is))) ! semi-implict integration scheme
1746  IF (vsin(is).GT.0.) whitecap(3) = whitecap(3) + spec(is) * factor
1747  hstot = hstot + spec(is) * factor
1748  END DO
1749  END DO
1750  whitecap(3) = 4. * sqrt(whitecap(3))
1751  hstot =4.*sqrt(hstot)
1752  tauwix = tauwix + tauwx * drat * dt
1753  tauwiy = tauwiy + tauwy * drat * dt
1754  tauwnx = tauwnx + tauwax * drat * dt
1755  tauwny = tauwny + tauway * drat * dt
1756  ! MISSING: TAIL TO BE ADDED ?
1757  !
1758 #ifdef W3_NLS
1759  CALL w3snls ( spec, cg1, wn1, depth, u10abs, dt, aa=spec )
1760 #endif
1761  !
1762  ! 6. Add tail ------------------------------------------------------- *
1763  ! a Mean parameters
1764  !
1765  !
1766 #ifdef W3_ST0
1767  CALL w3spr0 (spec, cg1, wn1, emean, fmean, wnmean, amax)
1768 #endif
1769 #ifdef W3_ST1
1770  CALL w3spr1 (spec, cg1, wn1, emean, fmean, wnmean, amax)
1771 #endif
1772 #ifdef W3_ST2
1773  CALL w3spr2 (spec, cg1, wn1, depth, fpi, u10abs, ustar, &
1774  emean, fmean, wnmean, amax, alpha, fp )
1775 #endif
1776 #ifdef W3_ST3
1777  CALL w3spr3 (spec, cg1, wn1, emean, fmean, fmeans, &
1778  wnmean, amax, u10abs, u10dir, ustar, ustdir, &
1779  tauwx, tauwy, cd, z0, charn, llws, fmeanws)
1780 #endif
1781 #ifdef W3_ST4
1782  CALL w3spr4 (spec, cg1, wn1, emean, fmean, fmean1, wnmean,&
1783  amax, u10abs, u10dir, &
1784 #ifdef W3_FLX5
1785  taua, tauadir, dair, &
1786 #endif
1787  ustar, ustdir, &
1788  tauwx, tauwy, cd, z0, charn, llws, fmeanws, dlwmean)
1789 #endif
1790 #ifdef W3_ST6
1791  CALL w3spr6 (spec, cg1, wn1, emean, fmean, wnmean, amax, fp)
1792 #endif
1793  !
1794 #ifdef W3_FLX2
1795  CALL w3flx2 ( zwnd, depth, fp, u10abs, u10dir, &
1796  ustar, ustdir, z0, cd )
1797 #endif
1798 #ifdef W3_FLX3
1799  CALL w3flx3 ( zwnd, depth, fp, u10abs, u10dir, &
1800  ustar, ustdir, z0, cd )
1801 #endif
1802  !
1803 #ifdef W3_ST1
1804  fh1 = fxfm * fmean
1805  fhigh = min( sig(nk) , max( fh1 , fh2 ) )
1806  nkh = max( 2 , min( nkh1 , &
1807  int( facti2 + facti1*log(max(1.e-7,fhigh)) ) ) )
1808  !
1809  IF ( fltest ) WRITE (ndst,9060) &
1810  fh1*tpiinv, fh2*tpiinv, fhigh*tpiinv, nkh
1811 #endif
1812  !
1813 #ifdef W3_ST2
1814  fhtran = xft*fpi
1815  fhigh = xfc*fpi
1816  dfh = fhigh - fhtran
1817  nkh = max( 1 , &
1818  int( facti2 + facti1*log(max(1.e-7,fhtran)) ) )
1819  !
1820  IF ( fltest ) WRITE (ndst,9061) fhtran, fhigh, nkh
1821 #endif
1822  !
1823 #ifdef W3_ST3
1824  fh1 = ffxfm * fmean
1825  fh2 = ffxpm / ustar
1826  fhigh = min( sig(nk) , max( fh1 , fh2 ) )
1827  nkh = max( 2 , min( nkh1 , &
1828  int( facti2 + facti1*log(max(1.e-7,fhigh)) ) ) )
1829  !
1830  IF ( fltest ) WRITE (ndst,9062) &
1831  fh1*tpiinv, fh2*tpiinv, fhigh*tpiinv, nkh
1832 #endif
1833  !
1834 #ifdef W3_ST4
1835  ! Introduces a Long & Resio (JGR2007) type dependance on wave age
1836  fage = ffxfa*tanh(0.3*u10abs*fmeanws*tpi/grav)
1837  fh1 = (ffxfm+fage) * fmean1
1838  fh2 = ffxpm / ustar
1839  fhigh = min( sig(nk) , max( fh1 , fh2 ) )
1840  nkh = max( 2 , min( nkh1 , &
1841  int( facti2 + facti1*log(max(1.e-7,fhigh)) ) ) )
1842 #endif
1843  !
1844 #ifdef W3_ST6
1845  IF (fxfm .LE. 0) THEN
1846  fhigh = sig(nk)
1847  ELSE
1848  fhigh = min( sig(nk), max(fxfm * fmean, fxpm / ustar) )
1849  ENDIF
1850  nkh = max( 2 , min( nkh1 , &
1851  int( facti2 + facti1*log(max(1.e-7,fhigh)) ) ) )
1852  !
1853  IF ( fltest ) WRITE (ndst,9063) fhigh*tpiinv, nkh
1854 #endif
1855  !
1856  ! 6.b Limiter for shallow water or Miche style criterion
1857  ! Last time step ONLY !
1858  ! uses true depth (D_INP) instead of limited depth
1859  !
1860 #ifdef W3_MLIM
1861  IF ( dttot .GE. 0.9999*dtg ) THEN
1862  hm = fhmax *tanh(wnmean*max(0.,d_inp)) / max(1.e-4,wnmean )
1863  em = hm * hm / 16.
1864  IF ( emean.GT.em .AND. emean.GT.1.e-30 ) THEN
1865  spec = spec / emean * em
1866  emean = em
1867  END IF
1868  END IF
1869 #endif
1870  !
1871  ! 6.c Seeding of spectrum
1872  ! alpha = 0.005 , 0.5 in eq., 0.25 for directional distribution
1873  !
1874 #ifdef W3_SEED
1875  DO ik=min(nk,nkh), nk
1876  uc = facsd * grav / sig(ik)
1877  slev = min( 1. , max( 0. , u10abs/uc-1. ) ) * &
1878  6.25e-4 / wn1(ik)**3 / sig(ik)
1879  IF (inflags2(4)) slev=slev*(1-ice)
1880  DO ith=1, nth
1881  spec(ith+(ik-1)*nth) = max( spec(ith+(ik-1)*nth) , &
1882  slev * max( 0. , cos(u10dir-th(ith)) )**2 )
1883  END DO
1884  END DO
1885 #endif
1886  !
1887  ! 6.d Add tail
1888  !
1889  DO ik=nkh+1, nk
1890 #ifdef W3_ST2
1891  facdia = max( 0. , min( 1., (sig(ik)-fhtran)/dfh) )
1892  facpar = max( 0. , 1.-facdia )
1893 #endif
1894  DO ith=1, nth
1895  spec(ith+(ik-1)*nth) = spec(ith+(ik-2)*nth) * fachfa &
1896 #ifdef W3_ST2
1897  * facdia + facpar * spec(ith+(ik-1)*nth) &
1898 #endif
1899  + 0.
1900  END DO
1901  END DO
1902  !
1903  ! 6.e Update wave-supported stress----------------------------------- *
1904  !
1905 #ifdef W3_ST3
1906  CALL w3sin3 ( spec, cg1, wn2, u10abs, ustar, drat, as, &
1907  u10dir, z0, cd, tauwx, tauwy, tauwax, tauway, &
1908  ice, vsin, vdin, llws, ix, iy )
1909 #endif
1910 #ifdef W3_ST4
1911  CALL w3sin4 ( spec, cg1, wn2, u10abs, ustar, drat, as, &
1912  u10dir, z0, cd, tauwx, tauwy, tauwax, tauway, &
1913  vsin, vdin, llws, ix, iy, brlambda )
1914  IF (sintailpar(4).LT.0.5) CALL w3spr4 (spec, cg1, wn1, emean, fmean, fmean1, wnmean,&
1915  amax, u10abs, u10dir, &
1916 #ifdef W3_FLX5
1917  taua, tauadir, dair, &
1918 #endif
1919  ustar, ustdir, &
1920  tauwx, tauwy, cd, z0, charn, llws, fmeanws, dlwmean)
1921 #endif
1922 
1923  !
1924  ! 7. Check if integration complete ---------------------------------- *
1925  !
1926  ! Update QI5TSTART (Q. Liu)
1927 #ifdef W3_NL5
1928  CALL tick21(qi5tstart, dt)
1929 #endif
1930 
1931  IF (srce_call .eq. srce_imp_post) THEN
1932  EXIT
1933  ENDIF
1934 
1935  IF ( dttot .GE. 0.9999*dtg ) THEN
1936  ! IF (IX == DEBUG_NODE) WRITE(*,*) 'DTTOT, DTG', DTTOT, DTG
1937  EXIT
1938  ENDIF
1939 
1940  END DO ! INTEGRATION LOOP
1941 
1942 #ifdef W3_DEBUGSRC
1943  IF (ix .eq. debug_node) THEN
1944  WRITE(740+iaproc,*) 'NSTEPS=', nsteps
1945  WRITE(740+iaproc,*) '1 : sum(SPEC)=', sum(spec)
1946  END IF
1947  WRITE(740+iaproc,*) 'DT=', dt, 'DTG=', dtg
1948 #endif
1949  !
1950  ! ... End point dynamic integration - - - - - - - - - - - - - - - - - -
1951  !
1952  ! 8. Save integration data ------------------------------------------ *
1953  !
1954  dtdyn = dtdyn / real(max(1,nsteps))
1955  fcut = fhigh * tpiinv
1956  !
1957  GOTO 888
1958  !
1959  ! Error escape locations
1960  !
1961 #ifdef W3_NNT
1962 800 CONTINUE
1963  WRITE (ndse,8000) fname, ierr
1964  CALL extcde (1)
1965  !
1966 801 CONTINUE
1967  WRITE (ndse,8001) ierr
1968  CALL extcde (2)
1969 #endif
1970  !
1971 888 CONTINUE
1972  !
1973  ! 9.a Computes PHIOC------------------------------------------ *
1974  ! The wave to ocean flux is the difference between initial energy
1975  ! and final energy, plus wind input plus the SNL flux to high freq.,
1976  ! minus the energy lost to the bottom boundary layer (BBL)
1977  !
1978 #ifdef W3_DEBUGSRC
1979  IF (ix .eq. debug_node) THEN
1980  WRITE(740+iaproc,*) '2 : sum(SPEC)=', sum(spec)
1981  END IF
1982 #endif
1983  efinish = 0.
1984  mwxfinish = 0.
1985  mwyfinish = 0.
1986  DO ik=1, nk
1987  eband = 0.
1988  a1band = 0.
1989  b1band = 0.
1990  DO ith=1, nth
1991  diff = specinit(ith+(ik-1)*nth)-spec(ith+(ik-1)*nth)
1992  eband = eband + diff
1993  a1band = a1band + diff*ecos(ith)
1994  b1band = b1band + diff*esin(ith)
1995  END DO
1996  efinish = efinish + eband * dden(ik) / cg1(ik)
1997  mwxfinish = mwxfinish + a1band * dden(ik) / cg1(ik) &
1998  * wn1(ik)/sig(ik)
1999  mwyfinish = mwyfinish + b1band * dden(ik) / cg1(ik) &
2000  * wn1(ik)/sig(ik)
2001  END DO
2002  !
2003  ! Transformation in momentum flux in m^2 / s^2
2004  !
2005  tauox=(grav*mwxfinish+tauwix-taubbl(1))/dtg
2006  tauoy=(grav*mwyfinish+tauwiy-taubbl(2))/dtg
2007  tauwix=tauwix/dtg
2008  tauwiy=tauwiy/dtg
2009  tauwnx=tauwnx/dtg
2010  tauwny=tauwny/dtg
2011  taubbl(:)=taubbl(:)/dtg
2012  tauocx=dair*coef*coef*ustar*ustar*cos(ustdir) + dwat*(tauox-tauwix)
2013  tauocy=dair*coef*coef*ustar*ustar*sin(ustdir) + dwat*(tauoy-tauwiy)
2014  !
2015  ! Transformation in wave energy flux in W/m^2=kg / s^3
2016  !
2017  phioc =dwat*grav*(efinish+phiaw-phibbl)/dtg
2018  phiaw =dwat*grav*phiaw /dtg
2019  phinl =dwat*grav*phinl /dtg
2020  phibbl=dwat*grav*phibbl/dtg
2021  !
2022  ! 10.1 Adds ice scattering and dissipation: implicit integration---------------- *
2023  ! INFLAGS2(4) is true if ice concentration was ever read during
2024  ! this simulation
2025  !
2026 #ifdef W3_DEBUGSRC
2027  IF (ix .eq. debug_node) THEN
2028  WRITE(740+iaproc,*) '3 : sum(SPEC)=', sum(spec)
2029  END IF
2030 #endif
2031 
2032  IF ( inflags2(4).AND.ice.GT.0 ) THEN
2033 
2034  IF (iicedisp) THEN
2035  icecoef2 = 1e-6
2036  CALL liu_forward_dispersion (iceh,icecoef2,depth, &
2037  sig,wn_r,cg_ice,alpha_liu)
2038  !
2039  IF (iicesmooth) THEN
2040 #ifdef W3_IS2
2041  DO ik=1,nk
2042  smooth_icedisp=0.
2043  IF (is2pars(14)*(tpi/wn_r(ik)).LT.icef) THEN ! IF ICE IS NOT TOO MUCH BROKEN
2044  smooth_icedisp=tanh((icef-is2pars(14)*(tpi/wn_r(ik)))/(icef*is2pars(13)))
2045  END IF
2046  wn_r(ik)=wn1(ik)*(1-smooth_icedisp)+wn_r(ik)*(smooth_icedisp)
2047  END DO
2048 #endif
2049  END IF
2050  ELSE
2051  wn_r=wn1
2052  cg_ice=cg1
2053  END IF
2054  !
2055  r(:)=1 ! In case IC2 is defined but not IS2
2056  !
2057 #ifdef W3_IC1
2058  CALL w3sic1 ( spec,depth, cg1, ix, iy, vsic, vdic )
2059 #endif
2060 #ifdef W3_IS2
2061  CALL w3sis2 ( spec, depth, ice, iceh, icef, icedmax, ix, iy, &
2062  vsir, vdir, vdir2, wn1, cg1, wn_r, cg_ice, r )
2063 #endif
2064 #ifdef W3_IC2
2065  CALL w3sic2 ( spec, depth, iceh, icef, cg1, wn1,&
2066  ix, iy, vsic, vdic, wn_r, cg_ice, alpha_liu, r)
2067 #endif
2068 #ifdef W3_IC3
2069  CALL w3sic3 ( spec,depth, cg1, wn1, ix, iy, vsic, vdic )
2070 #endif
2071 #ifdef W3_IC4
2072  CALL w3sic4 ( spec,depth, cg1, ix, iy, vsic, vdic )
2073 #endif
2074 #ifdef W3_IC5
2075  CALL w3sic5 ( spec,depth, cg1, wn1, ix, iy, vsic, vdic )
2076 #endif
2077  !
2078 #ifdef W3_IS1
2079  CALL w3sis1 ( spec, ice, vsir )
2080 #endif
2081  spec2 = spec
2082  !
2083  tauice(:) = 0.
2084  phice = 0.
2085  DO ik=1,nk
2086  is = 1+(ik-1)*nth
2087  !
2088  ! First part of ice term integration: dissipation part
2089  !
2090  att=1.
2091 #ifdef W3_IC1
2092  att=exp(ice*vdic(is)*dtg)
2093 #endif
2094 #ifdef W3_IC2
2095  att=exp(ice*vdic(is)*dtg)
2096 #endif
2097 #ifdef W3_IC3
2098  att=exp(ice*vdic(is)*dtg)
2099 #endif
2100 #ifdef W3_IC4
2101  att=exp(ice*vdic(is)*dtg)
2102 #endif
2103 #ifdef W3_IC5
2104  att=exp(ice*vdic(is)*dtg)
2105 #endif
2106 #ifdef W3_IS1
2107  att=att*exp(ice*vdir(is)*dtg)
2108 #endif
2109 #ifdef W3_IS2
2110  att=att*exp(ice*vdir2(is)*dtg)
2111  IF (is2pars(2).EQ.0) THEN ! Reminder : IS2PARS(2) = IS2BACKSCAT
2112  !
2113  ! If there is not re-distribution in directions the scattering is just an attenuation
2114  !
2115  att=att*exp((ice*vdir(is))*dtg)
2116  END IF
2117 #endif
2118  spec(1+(ik-1)*nth:nth+(ik-1)*nth) = att*spec2(1+(ik-1)*nth:nth+(ik-1)*nth)
2119  !
2120  ! Second part of ice term integration: scattering including re-distribution in directions
2121  !
2122 #ifdef W3_IS2
2123  IF (is2pars(2).GE.0) THEN
2124  IF (is2pars(20).GT.0.5) THEN
2125  !
2126  ! Case of isotropic back-scatter: the directional spectrum is decomposed into
2127  ! - an isotropic part (ISO): eigenvalue of scattering is 0
2128  ! - the rest (SPEC-ISO): eigenvalue of scattering is VDIR(IS)
2129  !
2130  scat = exp(vdir(is)*is2pars(2)*dtg)
2131  iso = sum(spec(1+(ik-1)*nth:nth+(ik-1)*nth))/nth
2132  spec(1+(ik-1)*nth:nth+(ik-1)*nth) = iso &
2133  +(spec(1+(ik-1)*nth:nth+(ik-1)*nth)-iso)*scat
2134  ELSE
2135  !
2136  ! General solution with matrix exponentials: same as bottom scattering, see Ardhuin & Herbers (JFM 2002)
2137  !
2138  scatspec(1:nth)=dble(spec(1+(ik-1)*nth:nth+(ik-1)*nth))
2139  spec(1+(ik-1)*nth:nth+(ik-1)*nth) = &
2140  REAL(MATMUL(IS2EIGVEC(:,:), EXP(IS2EIGVAL(:)*VDIR(IS)*DTG*IS2PARS(2)) &
2141  *MATMUL(TRANSPOSE(IS2EIGVEC(:,:)),SCATSPEC)))
2142  END IF
2143  END IF
2144 #endif
2145  !
2146  ! 10.2 Fluxes of energy and momentum due to ice effects
2147  !
2148  factor = dden(ik)/cg1(ik) !Jacobian to get energy in band
2149  factor2= factor*grav*wn1(ik)/sig(ik) ! coefficient to get momentum
2150  DO ith = 1,nth
2151  is = ith+(ik-1)*nth
2152  phice = phice + (spec(is)-spec2(is)) * factor
2153  cosi(1)=ecos(is)
2154  cosi(2)=esin(is)
2155  tauice(:) = tauice(:) - (spec(is)-spec2(is))*factor2*cosi(:)
2156  END DO
2157  END DO
2158  phice =-1.*dwat*grav*phice /dtg
2159  tauice(:)=tauice(:)/dtg
2160  ELSE
2161 #ifdef W3_IS2
2162  IF (is2pars(10).LT.0.5) THEN
2163  icef = 0.
2164  ENDIF
2165 #endif
2166  END IF
2167  !
2168  !
2169  ! - - - - - - - - - - - - - - - - - - - - - -
2170  ! 11. Sea state dependent stress routine calls
2171  ! - - - - - - - - - - - - - - - - - - - - - -
2172  !Note the Sea-state dependent stress calculations are primarily for high-wind
2173  !conditions (>10 m/s). It is not recommended to use these at lower wind
2174  !in their current state.
2175  !
2176 #ifdef W3_DEBUGSRC
2177  IF (ix .eq. debug_node) THEN
2178  WRITE(740+iaproc,*) '4 : sum(SPEC)=', sum(spec)
2179  END IF
2180 #endif
2181 
2182  ! FLD1/2 requires the calculation of FPI:
2183 #ifdef W3_FLD1
2184  CALL calc_fpi(spec, cg1, fpi, vsin )
2185 #endif
2186 #ifdef W3_FLD2
2187  CALL calc_fpi(spec, cg1, fpi, vsin )
2188 #endif
2189  !
2190 #ifdef W3_FLD1
2191  IF (u10abs.GT.10. .and. hstot.gt.0.5) then
2192  CALL w3fld1 ( spec,min(fpi/tpi,2.0),coef*u10abs*cos(u10dir), &
2193  coef*u10abs*sin(u10dir), zwnd, depth, 0.0, &
2194  dair, ustar, ustdir, z0,taunux,taunuy,charn)
2195  ELSE
2196  charn = aalpha
2197  ENDIF
2198 #endif
2199 #ifdef W3_FLD2
2200  IF (u10abs.GT.10. .and. hstot.gt.0.5) then
2201  CALL w3fld2 ( spec,min(fpi/tpi,2.0),coef*u10abs*cos(u10dir), &
2202  coef*u10abs*sin(u10dir), zwnd, depth, 0.0, &
2203  dair, ustar, ustdir, z0,taunux,taunuy,charn)
2204  ELSE
2205  charn = aalpha
2206  ENDIF
2207 #endif
2208  !
2209  ! 12. includes shoreline reflection --------------------------------------------- *
2210  !
2211 #ifdef W3_DEBUGSRC
2212  IF (ix .eq. debug_node) THEN
2213  WRITE(740+iaproc,*) '5 : sum(SPEC)=', sum(spec)
2214  END IF
2215 #endif
2216 
2217 #ifdef W3_REF1
2218  IF (reflec(1).GT.0.OR.reflec(2).GT.0.OR.(reflec(4).GT.0.AND.berg.GT.0)) THEN
2219  CALL w3sref ( spec, cg1, wn1, emean, fmean, depth, cx, cy, &
2220  reflec, refled, trnx, trny, &
2221  berg, dtg, ix, iy, jsea, vref )
2222  IF (gtype.EQ.ungtype.AND.refpars(3).LT.0.5) THEN
2223 #ifdef W3_PDLIB
2224  IF (iobp_loc(jsea).EQ.0) THEN
2225 #else
2226  IF (iobp(ix).EQ.0) THEN
2227 #endif
2228  DO ik=1, nk
2229  DO ith=1, nth
2230  isp = ith+(ik-1)*nth
2231 #ifdef W3_PDLIB
2232  IF (iobpd_loc(ith,jsea).EQ.0) spec(isp) = dtg*vref(isp)
2233 #else
2234  IF (iobpd(ith,ix).EQ.0) spec(isp) = dtg*vref(isp)
2235 #endif
2236  END DO
2237  END DO
2238  ELSE
2239  spec(:) = spec(:) + dtg * vref(:)
2240  ENDIF
2241  ELSE
2242  spec(:) = spec(:) + dtg * vref(:)
2243  END IF
2244  END IF
2245 #endif
2246 
2247  !
2248 #ifdef W3_DEBUGSRC
2249  IF (ix .eq. debug_node) THEN
2250  WRITE(740+iaproc,*) '6 : sum(SPEC)=', sum(spec)
2251  END IF
2252 #endif
2253 
2254  first = .false.
2255 
2256  IF (it.EQ.0) spec = specinit
2257 
2258  spec = max(0., spec)
2259  !
2260  RETURN
2261  !
2262  ! Formats
2263  !
2264 #ifdef W3_NNT
2265 8000 FORMAT (/' *** ERROR W3SRCE : ERROR IN OPENING FILE ',a,' ***'/ &
2266  ' IOSTAT = ',i10/)
2267 8001 FORMAT (/' *** ERROR W3SRCE : ERROR IN WRITING TO FILE ***'/ &
2268  ' IOSTAT = ',i10/)
2269 #endif
2270  !
2271 #ifdef W3_T
2272 9000 FORMAT (' TEST W3SRCE : COUNTERS : NO LONGER AVAILABLE')
2273 9001 FORMAT (' TEST W3SRCE : DEPTH :',f8.1/ &
2274  ' WIND SPEED :',f8.1/ &
2275  ' WIND DIR :',f8.1)
2276 #endif
2277 #ifdef W3_ST1
2278 9004 FORMAT (' TEST W3SRCE : FHIGH (3X) : ',3f8.4/ &
2279  ' ------------- NEW DYNAMIC INTEGRATION LOOP', &
2280  ' ------------- ')
2281 #endif
2282 #ifdef W3_ST2
2283 9005 FORMAT (' TEST W3SRCE : FHIGH : ',f8.4/ &
2284  ' ------------- NEW DYNAMIC INTEGRATION LOOP', &
2285  ' ------------- ')
2286 #endif
2287 #ifdef W3_ST3
2288 9006 FORMAT (' TEST W3SRCE : FHIGH (3X) : ',3f8.4/ &
2289  ' ------------- NEW DYNAMIC INTEGRATION LOOP', &
2290  ' ------------- ')
2291 #endif
2292 #ifdef W3_ST4
2293 9006 FORMAT (' TEST W3SRCE : FHIGH (3X) : ',3f8.4/ &
2294  ' ------------- NEW DYNAMIC INTEGRATION LOOP', &
2295  ' ------------- ')
2296 #endif
2297  !
2298 #ifdef W3_T
2299 9020 FORMAT (' TEST W3SRCE : NSTEP : ',i4,' DTTOT :',f6.1)
2300 9021 FORMAT (' TEST W3SRCE : NKH (3X) : ',2i3,i6)
2301 9040 FORMAT (' TEST W3SRCE : DTRAW, DT, SHAVE :',2f6.1,2x,l1)
2302 #endif
2303  !
2304 #ifdef W3_ST1
2305 9060 FORMAT (' TEST W3SRCE : FHIGH (3X) : ',3f8.4/ &
2306  ' NKH : ',i3)
2307 #endif
2308 #ifdef W3_ST2
2309 9061 FORMAT (' TEST W3SRCE : FHIGH (2X) : ',2f8.4/ &
2310  ' NKH : ',i3)
2311 #endif
2312 #ifdef W3_ST3
2313 9062 FORMAT (' TEST W3SRCE : FHIGH (3X) : ',3f8.4/ &
2314  ' NKH : ',i3)
2315 #endif
2316 #ifdef W3_ST4
2317 9062 FORMAT (' TEST W3SRCE : FHIGH (3X) : ',3f8.4/ &
2318  ' NKH : ',i3)
2319 #endif
2320 #ifdef W3_ST6
2321 9063 FORMAT (' TEST W3SRCE : FHIGH : ',f8.4/ &
2322  ' NKH : ',i3)
2323 #endif
2324  !/
2325  !/ End of W3SRCE ----------------------------------------------------- /
2326  !/
2327  END SUBROUTINE w3srce
2328  !/ ------------------------------------------------------------------- /
2329 
2343  SUBROUTINE calc_fpi( A, CG, FPI, S )
2344  !/
2345  !/ +-----------------------------------+
2346  !/ | WAVEWATCH III NOAA/NCEP |
2347  !/ | Jessica Meixner |
2348  !/ | |
2349  !/ | FORTRAN 90 |
2350  !/ | Last update : 06-Jun-2018 |
2351  !/ +-----------------------------------+
2352  !/
2353  !/ 06-Jul-2016 : Origination ( version 5.12 )
2354  !/ 06-Jul-2016 : Add SUBROUTINE SIGN_VSD_SEMI_IMPLICIT_WW3
2355  !/ Add optional DEBUGSRC/PDLIB ( version 6.04 )
2356  !/
2357  ! 1. Purpose :
2358  !
2359  ! Calculate equivalent peak frequency
2360  !
2361  ! 2. Method :
2362  !
2363  ! Tolman and Chalikov (1996), equivalent peak frequency from source
2364 
2365  ! 3. Parameters :
2366  !
2367  ! Parameter list
2368  ! ----------------------------------------------------------------
2369  ! A R.A. I Action density spectrum (1-D).
2370  ! CG R.A. I Group velocities for k-axis of spectrum.
2371  ! FPI R.A. O Input 'peak' frequency.
2372  ! S R.A. I Source term (1-D version).
2373  ! ----------------------------------------------------------------
2374  !
2375  ! 4. Subroutines used :
2376  !
2377  ! Name Type Module Description
2378  ! ----------------------------------------------------------------
2379  ! STRACE Subr. W3SERVMD Subroutine tracing.
2380  ! ----------------------------------------------------------------
2381  !
2382  ! 5. Called by :
2383  !
2384  ! Name Type Module Description
2385  ! ----------------------------------------------------------------
2386  ! W3SRCE Subr.
2387  ! ----------------------------------------------------------------
2388  !
2389  ! 6. Error messages :
2390  !
2391  ! 7. Remarks :
2392  !
2393  ! 8. Structure :
2394  !
2395  ! See source code.
2396  !
2397  ! 9. Switches :
2398  !
2399  ! !/S Enable subroutine tracing.
2400  !
2401  ! 10. Source code :
2402  !
2403  !/ ------------------------------------------------------------------- /
2404  USE constants
2405  USE w3gdatmd, ONLY: nk, nth, nspec, xfr, dden, sig,fte, fttr
2406 #ifdef W3_S
2407  USE w3servmd, ONLY: strace
2408 #endif
2409  !
2410  IMPLICIT NONE
2411  !/
2412  !/ ------------------------------------------------------------------- /
2413  !/ Parameter list
2414  !/
2415  REAL, INTENT(IN) :: A(NSPEC), CG(NK), S(NSPEC)
2416  REAL, INTENT(OUT) :: FPI
2417  !/
2418  !/ ------------------------------------------------------------------- /
2419  !/ Local parameters
2420  !/
2421  INTEGER :: IS, IK
2422 #ifdef W3_S
2423  INTEGER, SAVE :: IENT = 0
2424 #endif
2425  REAL :: M0, M1, SIN1A(NK)
2426  !/
2427  !/ ------------------------------------------------------------------- /
2428  !/
2429 #ifdef W3_S
2430  CALL strace (ient, 'CALC_FPI')
2431 #endif
2432  !
2433  ! Calculate FPI: equivalent peak frequncy from wind source term
2434  ! input
2435  !
2436  DO ik=1, nk
2437  sin1a(ik) = 0.
2438  DO is=(ik-1)*nth+1, ik*nth
2439  sin1a(ik) = sin1a(ik) + max( 0. , s(is) )
2440  END DO
2441  END DO
2442  !
2443  m0 = 0.
2444  m1 = 0.
2445  DO ik=1, nk
2446  sin1a(ik) = sin1a(ik) * dden(ik) / ( cg(ik) * sig(ik)**3 )
2447  m0 = m0 + sin1a(ik)
2448  m1 = m1 + sin1a(ik)/sig(ik)
2449  END DO
2450  !
2451  sin1a(nk) = sin1a(nk) / dden(nk)
2452  m0 = m0 + sin1a(nk) * fte
2453  m1 = m1 + sin1a(nk) * fttr
2454  IF ( m1 .LT. 1e-20 ) THEN
2455  fpi = xfr * sig(nk)
2456  ELSE
2457  fpi = m0 / m1
2458  END IF
2459 
2460  END SUBROUTINE calc_fpi
2461  !/ ------------------------------------------------------------------- /!
2462 
2474  SUBROUTINE sign_vsd_semi_implicit_ww3(SPEC, VS, VD)
2475  !/
2476  !/ +-----------------------------------+
2477  !/ | WAVEWATCH III NOAA/NCEP |
2478  !/ | |
2479  !/ | Aron Roland (BGS IT&E GmbH) |
2480  !/ | Mathieu Dutour-Sikiric (IRB) |
2481  !/ | |
2482  !/ | FORTRAN 90 |
2483  !/ | Last update : 01-June-2018 |
2484  !/ +-----------------------------------+
2485  !/
2486  !/ 01-June-2018 : Origination. ( version 6.04 )
2487  !/
2488  ! 1. Purpose : Put source term in matrix same as done always
2489  ! 2. Method :
2490  ! 3. Parameters :
2491  !
2492  ! Parameter list
2493  ! ----------------------------------------------------------------
2494  ! ----------------------------------------------------------------
2495  !
2496  ! 4. Subroutines used :
2497  !
2498  ! Name Type Module Description
2499  ! ----------------------------------------------------------------
2500  ! STRACE Subr. W3SERVMD Subroutine tracing.
2501  ! ----------------------------------------------------------------
2502  !
2503  ! 5. Called by :
2504  !
2505  ! Name Type Module Description
2506  ! ----------------------------------------------------------------
2507  ! ----------------------------------------------------------------
2508  !
2509  ! 6. Error messages :
2510  ! 7. Remarks
2511  ! 8. Structure :
2512  ! 9. Switches :
2513  !
2514  ! !/S Enable subroutine tracing.
2515  !
2516  ! 10. Source code :
2517  !
2518  !/ ------------------------------------------------------------------- /
2519 #ifdef W3_S
2520  USE w3servmd, ONLY: strace
2521 #endif
2522  !
2523  USE w3gdatmd, only : nth, nk, nspec
2524  IMPLICIT NONE
2525  !/
2526  !/ ------------------------------------------------------------------- /
2527  !/ Parameter list
2528  !/
2529  !/ ------------------------------------------------------------------- /
2530  !/ Local PARAMETERs
2531  !/
2532 #ifdef W3_S
2533  INTEGER, SAVE :: IENT = 0
2534 #endif
2535  !/
2536  !/ ------------------------------------------------------------------- /
2537  !/
2538 
2539  INTEGER :: ISP, ITH, IK, IS
2540  REAL, INTENT(IN) :: SPEC(NSPEC)
2541  REAL, INTENT(INOUT) :: VS(NSPEC), VD(NSPEC)
2542 #ifdef W3_S
2543  CALL strace (ient, 'SIGN_VSD_SEMI_IMPLICIT_WW3')
2544 #endif
2545  DO is=1,nspec
2546  vd(is) = min(0., vd(is))
2547  END DO
2548  END SUBROUTINE sign_vsd_semi_implicit_ww3
2549  !/ ------------------------------------------------------------------- /
2550 
2562  SUBROUTINE sign_vsd_patankar_ww3(SPEC, VS, VD)
2563  !/
2564  !/ +-----------------------------------+
2565  !/ | WAVEWATCH III NOAA/NCEP |
2566  !/ | |
2567  !/ | Aron Roland (BGS IT&E GmbH) |
2568  !/ | Mathieu Dutour-Sikiric (IRB) |
2569  !/ | |
2570  !/ | FORTRAN 90 |
2571  !/ | Last update : 01-June-2018 |
2572  !/ +-----------------------------------+
2573  !/
2574  !/ 01-June-2018 : Origination. ( version 6.04 )
2575  !/
2576  ! 1. Purpose : Put source term in matrix Patankar style (experimental)
2577  ! 2. Method :
2578  ! 3. Parameters :
2579  !
2580  ! Parameter list
2581  ! ----------------------------------------------------------------
2582  ! ----------------------------------------------------------------
2583  !
2584  ! 4. Subroutines used :
2585  !
2586  ! Name Type Module Description
2587  ! ----------------------------------------------------------------
2588  ! STRACE Subr. W3SERVMD Subroutine tracing.
2589  ! ----------------------------------------------------------------
2590  !
2591  ! 5. Called by :
2592  !
2593  ! Name Type Module Description
2594  ! ----------------------------------------------------------------
2595  ! ----------------------------------------------------------------
2596  !
2597  ! 6. Error messages :
2598  ! 7. Remarks
2599  ! 8. Structure :
2600  ! 9. Switches :
2601  !
2602  ! !/S Enable subroutine tracing.
2603  !
2604  ! 10. Source code :
2605  !
2606  !/ ------------------------------------------------------------------- /
2607 #ifdef W3_S
2608  USE w3servmd, ONLY: strace
2609 #endif
2610  !
2611 
2612  USE w3gdatmd, only : nth, nk, nspec
2613  IMPLICIT NONE
2614  !/
2615  !/ ------------------------------------------------------------------- /
2616  !/ Parameter list
2617  !/
2618  !/ ------------------------------------------------------------------- /
2619  !/ Local PARAMETERs
2620  !/
2621 #ifdef W3_S
2622  INTEGER, SAVE :: IENT = 0
2623 #endif
2624  !/
2625  !/ ------------------------------------------------------------------- /
2626  !/
2627  INTEGER :: ISP, ITH, IK, IS
2628  REAL, INTENT(IN) :: SPEC(NSPEC)
2629  REAL, INTENT(INOUT) :: VS(NSPEC), VD(NSPEC)
2630 #ifdef W3_S
2631  CALL strace (ient, 'SIGN_VSD_PATANKAR_WW3')
2632 #endif
2633  DO is=1,nspec
2634  vd(is) = min(0., vd(is))
2635  vs(is) = max(0., vs(is))
2636  END DO
2637  END SUBROUTINE sign_vsd_patankar_ww3
2638  !/
2639  !/ End of module W3SRCEMD -------------------------------------------- /
2640  !/
2641 END MODULE w3srcemd
pdlib_w3profsmd::aspar_diag_all
real, dimension(:,:), allocatable aspar_diag_all
Definition: w3profsmd_pdlib.F90:114
w3gdatmd::nk
integer, pointer nk
Definition: w3gdatmd.F90:1230
w3gdatmd::xflt
real, pointer xflt
Definition: w3gdatmd.F90:1245
constants::srce_imp_pre
integer, parameter srce_imp_pre
srce_imp_pre
Definition: constants.F90:98
w3sln1md
Definition: w3sln1md.F90:3
w3sln1md::w3sln1
subroutine w3sln1(K, FHIGH, USTAR, USDIR, S)
Definition: w3sln1md.F90:57
w3str1md::w3str1
subroutine w3str1(A, AOLD, CG, WN, DEPTH, IX, S, D)
Triad interaction source term computed using the Lumped Triad Appproximation (LTA) of Eldeberky (1996...
Definition: w3str1md.F90:184
w3src3md::w3sin3
subroutine w3sin3(A, CG, K, U, USTAR, DRAT, AS, USDIR, Z0, CD, TAUWX, TAUWY, TAUWNX, TAUWNY, ICE, S, D, LLWS, IX, IY)
Calculate diagonal and input source term for WAM4+ approach.
Definition: w3src3md.F90:386
w3gdatmd::swl6s6
logical, pointer swl6s6
Definition: w3gdatmd.F90:1338
w3uostmd
Parmeterization of the unresolved obstacles.
Definition: w3uostmd.F90:22
w3sbt9md::w3sbt9
subroutine w3sbt9(AC, H_WDEPTH, S, D, IX, IY)
Compute dissipation by viscous fluid mud using Ng (2000) (adapted from Erick Rogers code by Mark Orze...
Definition: w3sbt9md.F90:119
w3gdatmd::nspec
integer, pointer nspec
Definition: w3gdatmd.F90:1230
w3swldmd
Source term module for swell dissipation.
Definition: w3swldmd.F90:18
w3str1md
Module for inclusion of triad nonlinear interaction according to Eldeberky's (1996) Lumped Triad Inte...
Definition: w3str1md.F90:17
w3flx5md::w3flx5
subroutine w3flx5(ZWND, U10, U10D, TAUA, TAUADIR, RHOAIR, UST, USTD, Z0, CD, CHARN)
Unified process to obtain friction velocity and drag when stresses are an input (from atmospheric mod...
Definition: w3flx5md.F90:119
w3snl1md::w3snl1
subroutine w3snl1(A, CG, KDMEAN, S, D)
Calculate nonlinear interactions and the diagonal term of its derivative.
Definition: w3snl1md.F90:115
w3gdatmd::ungtype
integer, parameter ungtype
Definition: w3gdatmd.F90:626
w3gdatmd::dmin
real, pointer dmin
Definition: w3gdatmd.F90:1183
constants::srce_direct
integer, parameter srce_direct
srce_direct
Definition: constants.F90:96
w3wdatmd
Define data structures to set up wave model dynamic data for several models simultaneously.
Definition: w3wdatmd.F90:18
w3flx4md::w3flx4
subroutine w3flx4(ZWND, U10, U10D, UST, USTD, Z0, CD)
Flux/stress computations according to Hwang (JTECH, 2011).
Definition: w3flx4md.F90:106
w3gdatmd::fxfm
real, pointer fxfm
Definition: w3gdatmd.F90:1245
w3idatmd::inflags2
logical, dimension(:), pointer inflags2
Definition: w3idatmd.F90:260
pdlib_w3profsmd::aspar_jac
real, dimension(:,:), allocatable aspar_jac
Definition: w3profsmd_pdlib.F90:114
yownodepool::pdlib_i_diag
integer, dimension(:), allocatable, target, public pdlib_i_diag
Definition: yownodepool.F90:85
w3gdatmd::iobpd_loc
integer *1, dimension(:,:), pointer iobpd_loc
Definition: w3gdatmd.F90:1116
w3parall::lsloc
logical, parameter lsloc
Definition: w3parall.F90:89
w3uostmd::uost_srctrmcompute
subroutine, public uost_srctrmcompute(IX, IY, SPEC, CG, DT, U10ABS, U10DIR, S, D)
Estimates the UOST source term for a give spectrum.
Definition: w3uostmd.F90:299
w3sis1md::w3sis1
subroutine, public w3sis1(A, ICE, S)
Spectral reflection due to ice.
Definition: w3sis1md.F90:75
w3fld1md::w3fld1
subroutine w3fld1(ASPC, FPI, WNDX, WNDY, ZWND, DEPTH, RIB, DAIR, UST, USTD, Z0, TAUNUX, TAUNUY, CHARN)
Definition: w3fld1md.F90:96
w3gdatmd::sig
real, dimension(:), pointer sig
Definition: w3gdatmd.F90:1234
w3sbt1md::w3sbt1
subroutine w3sbt1(A, CG, WN, DEPTH, S, D)
Bottom friction source term according to the empirical JONSWAP formulation.
Definition: w3sbt1md.F90:89
w3gdatmd::ffxfa
real, pointer ffxfa
Definition: w3gdatmd.F90:1324
w3src2md::w3sds2
subroutine w3sds2(A, CG, K, FPI, USTAR, ALFA, S, D)
Calculate whitecapping source term and diagonal term of derivative.
Definition: w3src2md.F90:583
w3sdb1md
Dummy slot for bottom friction source term.
Definition: w3sdb1md.F90:24
w3odatmd::iaproc
integer, pointer iaproc
Definition: w3odatmd.F90:457
w3wdatmd::time
integer, dimension(:), pointer time
Definition: w3wdatmd.F90:172
w3flx2md::w3flx2
subroutine w3flx2(ZWIND, DEPTH, FP, U, UDIR, UST, USTD, Z0, CD)
FLux/stress computations according Tolman and Chalikov (1996).
Definition: w3flx2md.F90:91
w3ref1md::w3sref
subroutine w3sref(A, CG, WN, EMEAN, FMEAN, DEPTH, CX1, CY1, REFLC, REFLD, TRNX, TRNY, BERG, DT, IX, IY, JSEA, S)
Computes coastal and iceberg/island reflections and adds free IG energy.
Definition: w3ref1md.F90:113
w3gdatmd::ecos
real, dimension(:), pointer ecos
Definition: w3gdatmd.F90:1234
constants::rade
real, parameter rade
RADE Conversion factor from radians to degrees.
Definition: constants.F90:76
constants::srce_imp_post
integer, parameter srce_imp_post
srce_imp_post
Definition: constants.F90:97
w3sic2md::w3sic2
subroutine, public w3sic2(A, DEPTH, ICEH, ICEF, CG, WN, IX, IY, S, D, WN_R, CG_ICE, ALPHA, R)
S_{ice} source term using 5 parameters read from input files.
Definition: w3sic2md.F90:122
w3gdatmd::zzwnd
real, pointer zzwnd
Definition: w3gdatmd.F90:1311
w3gdatmd::fhmax
real, pointer fhmax
Definition: w3gdatmd.F90:1245
w3srcemd::w3srce
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
Calculate and integrate source terms for a single grid point.
Definition: w3srcemd.F90:195
w3odatmd::fnmpre
character(len=80) fnmpre
Definition: w3odatmd.F90:330
w3gdatmd::fssource
logical, pointer fssource
Definition: w3gdatmd.F90:1406
w3fld1md
Definition: w3fld1md.F90:2
w3sbt4md::w3sbt4
subroutine w3sbt4(A, CG, WN, DEPTH, D50, PSIC, TAUBBL, BEDFORM, S, D, IX, IY)
Computes the SHOWEX bottom friction with movable bed effects.
Definition: w3sbt4md.F90:341
w3fld2md
Definition: w3fld2md.F90:2
w3gdatmd::iobp_loc
integer *2, dimension(:), pointer iobp_loc
Definition: w3gdatmd.F90:1117
w3src4md::w3sds4
subroutine w3sds4(A, K, CG, USTAR, USDIR, DEPTH, DAIR, SRHS, DDIAG, IX, IY, BRLAMBDA, WHITECAP, DLWMEAN)
Calculate whitecapping source term and diagonal term of derivative.
Definition: w3src4md.F90:2034
w3wdatmd::va
real, dimension(:,:), pointer va
Definition: w3wdatmd.F90:183
w3src6md::w3sin6
subroutine, public w3sin6(A, CG, WN2, UABS, USTAR, USDIR, CD, DAIR, TAUWX, TAUWY, TAUNWX, TAUNWY, S, D)
Observation-based source term for wind input.
Definition: w3src6md.F90:292
w3gdatmd::th
real, dimension(:), pointer th
Definition: w3gdatmd.F90:1234
w3gdatmd::fxpm
real, pointer fxpm
Definition: w3gdatmd.F90:1245
w3snl4md::w3snl4
subroutine w3snl4(A, CG, WN, DEPTH, S, D)
Interface module for TSA type nonlinear interactions.
Definition: w3snl4md.F90:611
w3snlsmd::w3snls
subroutine w3snls(A, CG, WN, DEPTH, UABS, DT, SNL, AA)
High-frequeny filter based on the nonlinear interactions for an uresolved quadruplet.
Definition: w3snlsmd.F90:132
w3src2md::w3sin2
subroutine w3sin2(A, CG, K, U, UDIR, CD, Z0, FPI, S, D)
Calculate input source term.
Definition: w3src2md.F90:309
w3gdatmd::ffxfm
real, pointer ffxfm
Definition: w3gdatmd.F90:1311
yownodepool::pdlib_si
real(rkind), dimension(:), allocatable, target, public pdlib_si
Definition: yownodepool.F90:80
w3gdatmd::xrel
real, pointer xrel
Definition: w3gdatmd.F90:1245
w3gdatmd::xfc
real, pointer xfc
Definition: w3gdatmd.F90:1245
w3odatmd::ndse
integer, pointer ndse
Definition: w3odatmd.F90:456
w3flx3md
FLux/stress computations according Tolman and Chalikov (1996).
Definition: w3flx3md.F90:23
w3flx4md
Flux/stress computations according to Hwang ( 2011).
Definition: w3flx4md.F90:27
w3flx5md
Unified process to obtain friction velocity and drag when stresses are an input (from atmospheric mod...
Definition: w3flx5md.F90:28
w3gdatmd::fachfa
real, pointer fachfa
Definition: w3gdatmd.F90:1232
w3flx1md
Flux/stress computations according to Wu (1980).
Definition: w3flx1md.F90:21
w3flx2md
FLux/stress computations according Tolman and Chalikov (1996).
Definition: w3flx2md.F90:21
w3src4md::w3sin4
subroutine w3sin4(A, CG, K, U, USTAR, DRAT, AS, USDIR, Z0, CD, TAUWX, TAUWY, TAUWNX, TAUWNY, S, D, LLWS, IX, IY, BRLAMBDA)
Calculate diagonal and input source term for WAM4+ approach.
Definition: w3src4md.F90:426
w3sdb1md::w3sdb1
subroutine w3sdb1(IX, A, DEPTH, EMEAN, FMEAN, WNMEAN, CG, LBREAK, S, D)
Compute depth-induced breaking using Battjes and Janssen bore model approach.
Definition: w3sdb1md.F90:97
w3gdatmd::esin
real, dimension(:), pointer esin
Definition: w3gdatmd.F90:1234
w3gdatmd::iicesmooth
logical, pointer iicesmooth
Definition: w3gdatmd.F90:1217
yownodepool
Has data that belong to nodes.
Definition: yownodepool.F90:39
w3srcemd::calc_fpi
subroutine calc_fpi(A, CG, FPI, S)
Calculate equivalent peak frequency.
Definition: w3srcemd.F90:2344
w3servmd
Definition: w3servmd.F90:3
w3srcemd::sign_vsd_patankar_ww3
subroutine sign_vsd_patankar_ww3(SPEC, VS, VD)
Put source term in matrix Patankar style (experimental).
Definition: w3srcemd.F90:2563
w3gdatmd::facsd
real, pointer facsd
Definition: w3gdatmd.F90:1245
w3timemd::tick21
subroutine tick21(TIME, DTIME)
Definition: w3timemd.F90:84
w3sic4md::w3sic4
subroutine, public w3sic4(A, DEPTH, CG, IX, IY, S, D)
S_{ice} source term using 5 parameters read from input files.
Definition: w3sic4md.F90:121
w3snl5md::w3snl5
subroutine, public w3snl5(A, CG, WN, FMEAN, T1ABS, U10, UDIR, JSEA, S, D, KURT)
Interface to CalcQRSNL subroutine of the GKE module.
Definition: w3snl5md.F90:100
constants::tpiinv
real, parameter tpiinv
TPIINV Inverse of 2*Pi.
Definition: constants.F90:74
w3gdatmd::dtmin
real, pointer dtmin
Definition: w3gdatmd.F90:1183
w3snlsmd
Nonlinear interaction based ‘smoother’ for high frequencies.
Definition: w3snlsmd.F90:21
w3gdatmd::ffxpm
real, pointer ffxpm
Definition: w3gdatmd.F90:1311
w3gdatmd::nth
integer, pointer nth
Definition: w3gdatmd.F90:1230
w3dispmd::liu_forward_dispersion
subroutine liu_forward_dispersion(H_ICE, VISC, H_WDEPTH, SIGMA, K_SOLUTION, CG, ALPHA)
Definition: w3dispmd.F90:688
w3gdatmd::fttr
real, pointer fttr
Definition: w3gdatmd.F90:1232
w3odatmd
Definition: w3odatmd.F90:3
w3srcemd
Source term integration routine.
Definition: w3srcemd.F90:24
w3gdatmd::iicedisp
logical, pointer iicedisp
Definition: w3gdatmd.F90:1217
w3odatmd::screen
integer, pointer screen
Definition: w3odatmd.F90:456
w3parall::imem
integer, parameter imem
Definition: w3parall.F90:90
constants::dwat
real, parameter dwat
DWAT Density of water (kg/m3).
Definition: constants.F90:62
w3sbs1md::w3sbs1
subroutine w3sbs1(A, CG, WN, DEPTH, CX1, CY1, TAUSCX, TAUSCY, S, D)
Bottom scattering source term.
Definition: w3sbs1md.F90:114
w3gdatmd::optioncall
integer optioncall
Definition: w3gdatmd.F90:1110
w3sic5md::w3sic5
subroutine, public w3sic5(A, DEPTH, CG, WN, IX, IY, S, D)
Calculate ice source term S_{ice} according to 3 different sea ice models.
Definition: w3sic5md.F90:169
file
file(STRINGS ${CMAKE_BINARY_DIR}/switch switch_strings) separate_arguments(switches UNIX_COMMAND $
Definition: CMakeLists.txt:3
w3gdatmd::facp
real, pointer facp
Definition: w3gdatmd.F90:1245
w3gdatmd::b_jgs_limiter
logical, pointer b_jgs_limiter
Definition: w3gdatmd.F90:1413
w3gdatmd::iobpd
integer *1, dimension(:,:), pointer iobpd
Definition: w3gdatmd.F90:1130
w3gdatmd::b_jgs_limiter_func
integer, pointer b_jgs_limiter_func
Definition: w3gdatmd.F90:1417
w3gdatmd::xft
real, pointer xft
Definition: w3gdatmd.F90:1245
constants::tpi
real, parameter tpi
TPI 2*Pi.
Definition: constants.F90:72
w3src1md::w3sin1
subroutine w3sin1(A, K, USTAR, USDIR, S, D)
Calculate diagonal of input source (actual source term put together in W3SRCE).
Definition: w3src1md.F90:256
w3servmd::strace
subroutine strace(IENT, SNAME)
Definition: w3servmd.F90:148
w3src6md::w3spr6
subroutine, public w3spr6(A, CG, WN, EMEAN, FMEAN, WNMEAN, AMAX, FP)
Calculate mean wave parameters.
Definition: w3src6md.F90:122
w3gdatmd::facti2
real, pointer facti2
Definition: w3gdatmd.F90:1232
w3snl2md::w3snl2
subroutine w3snl2(A, CG, DEPTH, S, D)
Interface to exact interactions.
Definition: w3snl2md.F90:96
w3gdatmd::gtype
integer, pointer gtype
Definition: w3gdatmd.F90:1094
w3idatmd
Define data structures to set up wave model input data for several models simultaneously.
Definition: w3idatmd.F90:16
w3gdatmd::sintailpar
real, dimension(:), pointer sintailpar
Definition: w3gdatmd.F90:1324
w3gdatmd::xfr
real, pointer xfr
Definition: w3gdatmd.F90:1232
w3sis1md
Diffusion source term.
Definition: w3sis1md.F90:22
w3sis2md
Floe-size dependant scattering of waves in the marginal ice zone.
Definition: w3sis2md.F90:33
w3gdatmd::fte
real, pointer fte
Definition: w3gdatmd.F90:1232
w3sbt8md::w3sbt8
subroutine w3sbt8(AC, H_WDEPTH, S, D, IX, IY)
Compute dissipation by viscous fluid mud using Dalrymple and Liu (1978).
Definition: w3sbt8md.F90:112
w3gdatmd::iobp
integer *2, dimension(:), pointer iobp
Definition: w3gdatmd.F90:1129
w3sbt8md
Contains routines for computing dissipation by viscous fluid mud using Dalrymple and Liu (1978) "Thin...
Definition: w3sbt8md.F90:25
w3src3md::w3spr3
subroutine w3spr3(A, CG, WN, EMEAN, FMEAN, FMEANS, WNMEAN, AMAX, U, UDIR, USTAR, USDIR, TAUWX, TAUWY, CD, Z0, CHARN, LLWS, FMEANWS)
Calculate mean wave parameters for the use in the source term routines.
Definition: w3src3md.F90:137
w3odatmd::ndst
integer, pointer ndst
Definition: w3odatmd.F90:456
w3sbt9md
Contains routines for computing dissipation by viscous fluid mud using Ng (2000).
Definition: w3sbt9md.F90:25
w3src2md::w3spr2
subroutine w3spr2(A, CG, WN, DEPTH, FPI, U, USTAR, EMEAN, FMEAN, WNMEAN, AMAX, ALFA, FP)
Calculate mean wave parameters for the use in the source term routines (Tolman and Chalikov).
Definition: w3src2md.F90:127
w3gdatmd::icescales
real, dimension(:), pointer icescales
Definition: w3gdatmd.F90:1183
w3sbt4md
SHOWEX bottom friction source term (Ardhuin et al.
Definition: w3sbt4md.F90:25
w3sbt1md
JONSWAP bottom friction routine.
Definition: w3sbt1md.F90:21
w3src3md::w3sds3
subroutine w3sds3(A, K, CG, EMEAN, FMEAN, WNMEAN, USTAR, USDIR, DEPTH, S, D, IX, IY)
Calculate whitecapping source term and diagonal term of derivative.
Definition: w3src3md.F90:1255
constants
Define some much-used constants for global use (all defined as PARAMETER).
Definition: constants.F90:20
w3srcemd::sign_vsd_semi_implicit_ww3
subroutine sign_vsd_semi_implicit_ww3(SPEC, VS, VD)
Put source term in matrix same as done always.
Definition: w3srcemd.F90:2475
w3gdatmd::dden
real, dimension(:), pointer dden
Definition: w3gdatmd.F90:1234
w3gdatmd
Definition: w3gdatmd.F90:16
w3swldmd::w3swl6
subroutine, public w3swl6(A, CG, WN, S, D)
Turbulent dissipation of narrow-banded swell.
Definition: w3swldmd.F90:252
w3sbs1md
This module computes a scattering term based on the theory by Ardhuin and Magne (JFM 2007).
Definition: w3sbs1md.F90:22
w3sis2md::w3sis2
subroutine, public w3sis2(A, DEPTH, CICE, ICEH, ICEF, ICEDMAX, IX, IY, S, D, DISSIP, WN, CG, WN_R, CG_ICE, R)
Wave scattering in the MIZ, adapted from Dumont et al.
Definition: w3sis2md.F90:654
w3gdatmd::igpars
real, dimension(:), pointer igpars
Definition: w3gdatmd.F90:1142
w3gdatmd::refpars
real, dimension(:), pointer refpars
Definition: w3gdatmd.F90:1139
w3fld2md::w3fld2
subroutine w3fld2(ASPC, FPI, WNDX, WNDY, ZWND, DEPTH, RIB, DAIR, UST, USTD, Z0, TAUNUX, TAUNUY, CHARN)
Definition: w3fld2md.F90:67
w3servmd::extcde
subroutine extcde(IEXIT, UNIT, MSG, FILE, LINE, COMM)
Definition: w3servmd.F90:736
w3src6md::w3sds6
subroutine, public w3sds6(A, CG, WN, S, D)
Observation-based source term for dissipation.
Definition: w3src6md.F90:547
w3snl3md
Generalized and optimized multiple DIA implementation.
Definition: w3snl3md.F90:24
w3gdatmd::zwind
real, pointer zwind
Definition: w3gdatmd.F90:1304
w3snl4md
Generic shallow-water Boltzmann integral (FBI or TSA).
Definition: w3snl4md.F90:25
w3snl5md
Interface module for GKE (resonant & quasi-resonant four-wave interactions).
Definition: w3snl5md.F90:25
w3src3md
The 'WAM4+' source terms based on P.A.E.M.
Definition: w3src3md.F90:30
w3src4md
The 'SHOM/Ifremer' source terms based on P.A.E.M.
Definition: w3src4md.F90:28
w3gdatmd::aalpha
real, pointer aalpha
Definition: w3gdatmd.F90:1311
w3snl1md
Bundles routines to calculate nonlinear wave-wave interactions according to the Discrete Interaction ...
Definition: w3snl1md.F90:25
w3src6md
Observation-based wind input and dissipation after Donelan et al (2006), and Babanin et al.
Definition: w3src6md.F90:27
w3snl2md
Interface module to exact nonlinear interactions.
Definition: w3snl2md.F90:23
w3src0md
Mean wave parameter computation for case without input and dissipation.
Definition: w3src0md.F90:15
w3src1md
Bundle WAM cycle 3 input and dissipation source terms with their defining parameters.
Definition: w3src1md.F90:15
w3src2md
Tolman and Chalikov (1996) input and dissipation source terms.
Definition: w3src2md.F90:16
w3gdatmd::ftf
real, pointer ftf
Definition: w3gdatmd.F90:1232
w3gdatmd::is2pars
real, dimension(:), pointer is2pars
Definition: w3gdatmd.F90:1161
w3snl3md::w3snl3
subroutine w3snl3(A, CG, WN, DEPTH, S, D)
Multiple Discrete Interaction Parameterization for arbitrary depths with generalized quadruplet layou...
Definition: w3snl3md.F90:181
w3gdatmd::facti1
real, pointer facti1
Definition: w3gdatmd.F90:1232
w3snl1md::w3snlgqm
subroutine w3snlgqm(A, CG, WN, DEPTH, TSTOTn, TSDERn)
Definition: w3snl1md.F90:789
w3flx3md::w3flx3
subroutine w3flx3(ZWIND, DEPTH, FP, U, UDIR, UST, USTD, Z0, CD)
FLux/stress computations according Tolman and Chalikov (1996).
Definition: w3flx3md.F90:96
w3gdatmd::iqtpe
integer, pointer iqtpe
Definition: w3gdatmd.F90:1345
w3ref1md
This module computes shoreline reflection, and unresolved islands and iceberg reflections.
Definition: w3ref1md.F90:22
w3timemd
Definition: w3timemd.F90:3
w3src4md::w3spr4
subroutine w3spr4(A, CG, WN, EMEAN, FMEAN, FMEAN1, WNMEAN, AMAX, U, UDIR, ifdef W3_FLX5
Calculate mean wave parameters for the use in the source term routines.
Definition: w3src4md.F90:145
w3flx1md::w3flx1
subroutine w3flx1(ZWND, U10, U10D, UST, USTD, Z0, CD)
FLux/stress computations according to Wu (1980).
Definition: w3flx1md.F90:89
w3parall
Parallel routines for implicit solver.
Definition: w3parall.F90:22
w3dispmd
Definition: w3dispmd.F90:3
pdlib_w3profsmd::b_jac
real, dimension(:,:), allocatable b_jac
Definition: w3profsmd_pdlib.F90:114
w3sic3md::w3sic3
subroutine, public w3sic3(A, DEPTH, CG, WN, IX, IY, S, D)
Definition: w3sic3md.F90:97
w3src1md::w3sds1
subroutine w3sds1(A, K, EMEAN, FMEAN, WNMEAN, S, D)
Calculate whitecapping source term and diagonal term of derivative.
Definition: w3src1md.F90:433
w3src0md::w3spr0
subroutine w3spr0(A, CG, WN, EMEAN, FMEAN, WNMEAN, AMAX)
Calculate mean wave parameters.
Definition: w3src0md.F90:82
constants::grav
real, parameter grav
GRAV Acc.
Definition: constants.F90:61
pdlib_w3profsmd
Definition: w3profsmd_pdlib.F90:4
w3src1md::w3spr1
subroutine w3spr1(A, CG, WN, EMEAN, FMEAN, WNMEAN, AMAX)
Definition: w3src1md.F90:88
w3sic3md
Definition: w3sic3md.F90:3
w3sic4md
Calculate ice source term S_{ice} according to simple methods.
Definition: w3sic4md.F90:27
w3sic5md
Calculate ice source term S_{ice} according to different ice models:
Definition: w3sic5md.F90:25
w3gdatmd::dtmax
real, pointer dtmax
Definition: w3gdatmd.F90:1183
w3sic1md
Calculate ice source term S_{ice} according to simple methods.
Definition: w3sic1md.F90:23
w3sic2md
Calculate ice dissipation source term S_{ice}.
Definition: w3sic2md.F90:27
w3sic1md::w3sic1
subroutine, public w3sic1(A, DEPTH, CG, IX, IY, S, D)
S_{ice} source term using 5 parameters read from input files.
Definition: w3sic1md.F90:94