Add SCHISM-WWM setup for ensembles

docker/pyschism/docker/Dockerfile

@@ -47,10 +47,6 @@ RUN git clone https://github.com/schism-dev/pyschism.git && \
     conda run -p $ENV_PREFIX --no-capture-output \
         pip install "ensembleperturbation>=1.1.2"
 
-RUN mkdir -p $PROJECT_DIR/scripts
-COPY docker/*.py ${PROJECT_DIR}/scripts/
-ENV PYTHONPATH ${PROJECT_DIR}/scripts/
-
 ENV CONDA_DIR /opt/conda
 
 RUN conda clean --all
@@ -61,6 +57,11 @@ RUN chown -R $UID:$GID $HOME
 
 USER $USER
 
+RUN mkdir -p $PROJECT_DIR/scripts
+COPY docker/*.py ${PROJECT_DIR}/scripts/
+COPY docker/refs ${PROJECT_DIR}/refs/
+ENV PYTHONPATH ${PROJECT_DIR}/scripts/
+
 RUN mkdir -p $HOME/.local/share/pyschism
 
 # Ref: https://pythonspeed.com/articles/activate-conda-dockerfile/

docker/pyschism/docker/refs/param.nml

@@ -0,0 +1,69 @@
+&CORE
+  ipre=0
+  ibc=1
+  ibtp=0
+  nspool=24
+  ihfskip=11088
+  dt=150.0
+  rnday=19.25
+
+  msc2 = 24     !same as msc in .nml ... for consitency check between SCHISM and WWM
+  mdc2 = 30     !same as mdc in .nml
+/
+
+&OPT
+  start_year=2018
+  start_month=8
+  start_day=30
+  start_hour=6.0
+  utc_start=-0.0
+  ics=2
+  ihot=1
+  nchi=-1
+  hmin_man=1.0
+  ic_elev=1
+  nws=-1
+  wtiminc=150.0
+
+  icou_elfe_wwm = 1
+  nstep_wwm = 4  !call WWM every this many time steps
+  iwbl = 0 !wave boundary layer formulation (used only if USE_WMM and
+           !icou_elfe_wwm/=0 and nchi=1. If icou_elfe_wwm=0, set iwbl=0):
+           !1-modified Grant-Madsen formulation; 2-Soulsby (1997)
+  hmin_radstress = 1. !min. total water depth used only in radiation stress calculation [m]
+!  nrampwafo = 0       !ramp-up option for the wave forces (1: on; 0: off)
+  drampwafo = 0.      !ramp-up period in days for the wave forces (no ramp-up if <=0)
+  turbinj = 0.15      !% of depth-induced wave breaking energy injected in turbulence
+                      !(default: 0.15 (15%), as proposed by Feddersen, 2012)
+  turbinjds = 1.0     !% of wave energy dissipated through whitecapping injected in turbulence
+                      !(default: 1 (100%), as proposed by Paskyabi et al. 2012)
+  alphaw = 0.5        !for itur=4 : scaling parameter for the surface roughness z0s = alphaw*Hm0.
+                      !If negative z0s = abs(alphaw) e.g. z0s=0.2 m (Feddersen and Trowbridge, 2005)
+                         ! Vortex Force terms (off/on:0/1)
+/
+
+&SCHOUT
+  nhot=1
+  nhot_write=11088
+
+  iof_hydro(14) = 1
+
+  iof_wwm(1)  = 1 !sig. height (m) {sigWaveHeight}   2D
+  iof_wwm(2)  = 0 !Mean average period (sec) - TM01 {meanWavePeriod}  2D
+  iof_wwm(3)  = 0 !Zero down crossing period for comparison with buoy (s) - TM02 {zeroDowncrossPeriod}  2D
+  iof_wwm(4)  = 0 !Average period of wave runup/overtopping - TM10 {TM10}  2D
+  iof_wwm(5)  = 0 !Mean wave number (1/m) {meanWaveNumber}  2D
+  iof_wwm(6)  = 0 !Mean wave length (m) {meanWaveLength}  2D
+  iof_wwm(7)  = 0 !Mean average energy transport direction (degr) - MWD in NDBC? {meanWaveDirection}  2D
+  iof_wwm(8) = 0 !Mean directional spreading (degr) {meanDirSpreading}  2D
+  iof_wwm(9) = 1 !Discrete peak period (sec) - Tp {peakPeriod}  2D
+  iof_wwm(10) = 0 !Continuous peak period based on higher order moments (sec) {continuousPeakPeriod}  2D
+  iof_wwm(11) = 0 !Peak phase vel. (m/s) {peakPhaseVel}  2D
+  iof_wwm(12) = 0 !Peak n-factor {peakNFactor}   2D
+  iof_wwm(13) = 0 !Peak group vel. (m/s) {peakGroupVel}   2D
+  iof_wwm(14) = 0 !Peak wave number {peakWaveNumber}  2D
+  iof_wwm(15) = 0 !Peak wave length {peakWaveLength}  2D
+  iof_wwm(16) = 1 !Peak (dominant) direction (degr) {dominantDirection}  2D
+  iof_wwm(17) = 0 !Peak directional spreading {peakSpreading}  2D
+
+/

docker/pyschism/docker/setup_ensemble.py

@@ -26,6 +26,10 @@
 from pyschism.forcing import NWM
 from ensembleperturbation.perturbation.atcf import perturb_tracks
 
+
+import wwm
+
+
 logger = logging.getLogger(__name__)
 logger.setLevel(logging.INFO)
 
@@ -40,6 +44,7 @@ def main(args):
     nwm_file = EFS_MOUNT_POINT / args.nwm_file
     mesh_dir = EFS_MOUNT_POINT / args.mesh_directory
     hr_prelandfall = args.hours_before_landfall
+    use_wwm = args.use_wwm
 
     workdir = out_dir
     mesh_file = mesh_dir / 'mesh_w_bdry.grd'
@@ -191,6 +196,9 @@ def main(args):
         'parallel': True
     })
 
+    if use_wwm:
+        wwm.setup_wwm(mesh_file, workdir)
+
 
 def parse_arguments():
     argument_parser = ArgumentParser()
@@ -249,6 +257,9 @@ def parse_arguments():
     argument_parser.add_argument(
         "-b", "--hours-before-landfall", type=int
     )
+    argument_parser.add_argument(
+        "--use-wwm", action="store_true"
+    )
 
     argument_parser.add_argument(
         "name", help="name of the storm", type=str)