SCHISM coupling with wave model : WWM vs WW3

[dkobashi]

Hi,

This is not an issue, but an input I want to seek. I want to couple SCHISM with a wave model (for nearshore, shelf, and beyond) and I understand that SCHISM can be coupled with WWM. I also heard that SCHISM can be coupled with WaveWatch III (WW3) with an external coupler developed by a group in NOAA (ufs-weather-model: https://github.com/oceanmodeling/ufs-weather-model)

In order to decide which wave model to couple, I want to know the following.

1. Is WWM development/improvement still active and ongoing?
2. Is SCHISM coupling with WW3 getting mature and working well?
3. Are there any tools available to process input data (wind forcing and boundary conditions) for both WWM and WW3? (maybe for WW3, it could be more like a question to WW3 users...)

Since Pyschism and Pylibs do not seem to have a routine for WWM, I searched for documents/pre-processing tools and I found some documentations below.

SCHISM_WWM_Tips_and_Tricks.pdf
SCHISM_WWM3.pdf
Roland_2008_PHD_thesis.pdf

In SCHISM_WWM_Tips_and_Tricks.pdf, there are some python scripts so I can start trying to use them to see if they work to create input files. I checked SCHISM test cases relevant to WWM, but not all input files (e.g. hs, t02, fs, dir netcdf files) are available.

My current preference is to couple SCHISM with WW3. However, after a first glance of ufs-weather-model, the coupling looks very complicated and no instruction of SCHISM-WW3 coupling in detail.

Coupling with WWM appears to be easier than coupling with WW3 and no significant difference in performance between the two. As far as python scripts for preprocessing work, I may want to try WWM, but I am curious to know if anyone has experience in both and can share your thoughts with me.

Thanks in advance.

-DJ

[josephzhang8]

WWM is being maintained, and there is a rudimentary manual in the online manual. We are still trying to write a more comprehensive one with the help from community.

WW3 coupling via UFS is working for 2D radiation stress formulation, and NOAA folks have written some guide for that. @felicio93 could you please point him to the right github place. Thx.

[felicio93]

@dkobashi SCHISM+WW3 via the UFS is under development and the tickets can be found here. If you have specific questions I recommend submitting a ticket there.

Also, if you want to go the SCHISM+WWM route, I recently went over the process of setting up the model and I can help you with it.

[dkobashi]

Thanks @felicio93, @josephzhang8 for your comments,

I start with the SCHISM-WWM coupling.
For now, I am trying to run a python script to create wwmbnd.gr3 using wwmbnd.py I obtained from https://github.com/Krys1202/Guide_SCHISM/blob/main/Python_scripts/wwmbnd.py.

However, the python kernel died when I ran the script. It looks like the script uses lots of exit(), which seems to terminate the program. I thought the right way is to use exit(1) instead of exit()? Any comments?
I am also not 100% sure if I understand the structure of gen_wwmbnd.in used in the script.

4       !ncond, i.e. number of conditions
1 2 2   !cond1 #seg flag ob
2 2 1   !cond2 #seg flag land
3 22 -1 !cond3 #seg flag island
4 99 0  !cond4 #dummy flag other nodes

I have no issue to get netcdf files of hs, t02, spr, fp, and dir from CAWCR and set wwminput.nml. So as soon as I am able to create wwmbnd.gr3, I can test SCHISM-WWM coupling assuming WWM uses the same surface forcing files SHICMS uses with full coupling.

Thanks,

-DJ

[felicio93]

@dkobashi I created this tutorial on how to setup SCHISM+WWM using WW3 Spec bound. I don't think all of this apply to your case, but I am sure there are things you can use. I also have a section on how to setup the bnd file. Please let me know if you need anything else.

SCHISM+WWM Setup with WWIII Spectral Boundaries

Introduction

This is the step-by-step guide for setting up a SCHISM+WWM and force it with WaveWatchIII Global

This guide was prepared by Felicio Cassalho (felicio.cassalho@noaa.gov),

Justification:
SCHISM+WWM documentation is scarce. Adding an active wave boundary is a must in most of the open-ocean domains as it can have a large impact on the Hs, Elevations, etc.
The following tutorial covers:

• How to run WWIII Global and output station data that can be used as spectral boundaries for SCHISM+WWM
• How to setup SCHISM+WWM with spectral open wave boundaries.
  WWM can take two kinds of open wave boundaries: Bulk Parameters (dir, fp, hs, spr, t02) [NOT RECOMMENDED DUE TO LARGE ERRORS], which uses IBOUNDFORMAT = 3, or Wave Spectra [BEST OPTION], which uses IBOUNDFORMAT = 6.
  See the sample wwminput.nml files in the SCHISM repo for reference.

If you need more background on WWIII, here are some resources:

• WWIII Online Tutorial
• WWIII Compilation Tutorial (WW3 Workshop ERDC)
• WWIII/UFS Compilation Tutorial
• If you are running your model on Hercules and do not want recompile WW3 feel free to borrow model setup from: '/work2/noaa/nos-surge/felicioc/BeringSea/WW3/R03_Global_05/'

Step-by-step guide on how to run WaveWatchIII for SCHISM+WWM:

WWIII Inputs:

Grid definition:

• You can find a global WWIII grid (.msh) and the ww3_grid.nml in the sample run on Hercules: '/work2/noaa/nos-surge/felicioc/BeringSea/WW3/R03_Global_05/'
• You might want the WWIII spectral parameterization (&SPECTRUM_NML) to match was you are using in SCHISM+WWM (check wwminput.nml: MDC, MSC; and param.nml: msc2, mdc2)

Wind Forcing:

• Download ERA5 wind data:
  This is very standard, you can use something like:

import cdsapi

dataset = "reanalysis-era5-single-levels"
request = {
    "product_type": ["reanalysis"],
    "variable": [
        "10m_u_component_of_wind",
        "10m_v_component_of_wind",
    ],
    "year": ["2019"],
    "month": ["05","06","07","08","09","10","11"],
    "day": [
        "01", "02", "03",
        "04", "05", "06",
        "07", "08", "09",
        "10", "11", "12",
        "13", "14", "15",
        "16", "17", "18",
        "19", "20", "21",
        "22", "23", "24",
        "25", "26", "27",
        "28", "29", "30",
        "31"
    ],
    "time": [
        "00:00", "01:00", "02:00",
        "03:00", "04:00", "05:00",
        "06:00", "07:00", "08:00",
        "09:00", "10:00", "11:00",
        "12:00", "13:00", "14:00",
        "15:00", "16:00", "17:00",
        "18:00", "19:00", "20:00",
        "21:00", "22:00", "23:00"
    ],
    "data_format": "netcdf",
    "download_format": "unarchived"
}

client = cdsapi.Client()
client.retrieve(dataset, request).download()

• Process the wind data:

  • The global model won't run for long periods (more than +-2 months). So you have to prepare the runs accordingly (for instance by monthly periods). Here is how you can slice the wind data accordingly:
    ncks -d valid_time,0,1488 ORIGINAL_ERA5_FILE.nc wind_201905_to_2019_06.nc.

• Update the ww3_prnc.inp (or .nml):

  • Make sure you ww3_prnc.inp is calling the correct wind file name, dimensions (latitude, longitude, valid_time), and variables (u10, v10). If you are using something other than ERA5 these names might change.

Create the mod_def.ww3 and wind.ww3 executables:

• Submit a job where you call: srun /your/path/WW3/build/bin/ww3_grid ww3_grid.nml > ww3_grid.out and srun /your/path/WW3/build_time/bin/ww3_prnc > ww3_prnc.out
• You can find a sample .job file here (Hercules): /work2/noaa/nos-surge/felicioc/BeringSea/WW3/R03_Global_05/

Prepare station file (ww3_points.list):

This is the file that defines where you will export the spectral data that will be used to force SCHISM+WWM.
You want to place these output stations along your SCHISM+WWM open boundary. The SCHISM+WWM boundary is defined by wwmbnd.gr3. So first we have to create this file.

• wwmbnd.gr3:
  There are 3 kinds of boundaries: not on boundary/inactive (flag=0), Neumann (flag=3), active boundary (flag=2). What we want to do is to set every node to 0, except for where we want to add the wave boundary (set to 2).
  To create wwmbnd.gr3 you can run this script. You will also need a gen_wwmbnd.in, where the first line defines the number of conditions, the second the number of open boundaries (in my case 2) and the kind of boundary they are (in my case 2, which means active). The third line defines the other boundaries (i.e. land), and the kind of boundary they are (0 = inactive). The last line sets values for all the remaning nodes (set to 0).

• ww3_points.list:
  You can set the WWIII to output spectral data at your SCHISM+WWM boundary nodes (flag=2) with:

import numpy as np
import ocsmesh

def convert_longitude(lon, mode=1):
    """
    Converts longitude between formats:
    mode 1: [-180, 180] → [0, 360] (0 in Greenwich)
    mode 2: [0, 360] → [-180, 180] (0 in Greenwich)
    mode 3: [-180, 180] → [0, 360] (180 in Greenwich)
    """
    lon = np.asarray(lon)
    if mode == 1:
        return lon % 360
    elif mode == 2:
        return np.where(lon > 180, lon - 360, lon)
    elif mode == 3:
        return lon + 180
    return lon


rundir = r'/your/path/'
wwmbnd = ocsmesh.Mesh.open(rundir + 'wwmbnd.gr3', crs=4326)
lon, lat, b_flag, node = wwmbnd.coord[:,0], wwmbnd.coord[:,1], wwmbnd.value, wwmbnd.nodes.id()

lon = convert_longitude(lon, mode=2)

b_lat,b_lon,b_id = [],[],[]
for n, i in enumerate(b_flag):
    if abs(int(i[0])) == 2:
        b_lat.append(lat[n])
        b_lon.append(lon[n])
        b_id.append(str(node[n]))

with open(f'{rundir}/ww3_points.list', 'w') as f:
    for n, i in enumerate(b_id):
        f.write(f" {b_lon[n]:.5f} {b_lat[n]:.5f} '{b_id[n]}' \n")

Note that in my case my SCHISM+WWM model was referenced with longitudes from 0 to 360 with 0 in Greenwich, while WWIII was set with longitudes from 0 to 360 with 180 at Greenwich. So I had to apply a longitude conversion.

It might be the case that you don't need a spectral station at every boundary node, specially if your SCHISM+WWM has a much higher resolution than WWIII. I wrote this shell script that rewrites the ww3_points.list selecting every one out of every n (set to 7) stations (note: this only makes sense if your ww3_point.list stations are order correctly. There might be better ways to reduce the size of the ww3_point.lis, this is just what worked for my case).

Run WWIII:

Prepare ww3_shel.nml:

• In &DOMAIN_NML set your start and end dates
• In &INPUT_NML set the wind forcing to true (that will look for wind.ww3)
• In &OUTPUT_TYPE_NML define what you want to output, and the name of the point file (ww3_points.list)
• &OUTPUT_DATE_NML define the start and end dates, and the frequency you want the Field, Point, and Restart outputs. Set the frequency of point outputs to the same frequency you want for your SCHISM+WWM boundary forcing (e.g., 3600sec). If you want to generate a hotstart for a subsequent run, set the frequency of the restart files accordingly.
• see the ww3_shell.nml in /work2/noaa/nos-surge/felicioc/BeringSea/WW3/R03_Global_05/ for reference.

Submit Job:

• Submit a job where you call: srun /your/path/WW3/build/bin/ww3_shel > ww3_shel.out
• You can find a sample .job file here (Hercules): /work2/noaa/nos-surge/felicioc/BeringSea/WW3/R03_Global_05/

PostProcess WWIII Spectral Stations:

Prepare the spectral boundary files for SCHISM+WWM:

This is the tricky part. There is a mismatch in the way WWIII and WWM define direction. You have to "fix" the WWIII directions and concatenate all the files (over time and station) before they can be used in SCHISM+WWM.

• Concatenate by time (and fix direction): If you ran the model for multiple months, move all the *_spec.nc files to an directory and ran the following shell script:

#!/bin/bash

### First move all spec files to a directory (indir)
indir="/work2/noaa/nos-surge/felicioc/BeringSea/WW3/P03/spec_files"
outdir="/work2/noaa/nos-surge/felicioc/BeringSea/WW3/P03"
tmpdir="$outdir/tmp"
mkdir -p "$tmpdir"

# extract station id from file name
station_ids=$(ls $indir/ww3.*_spec.nc | sed -E 's/.*ww3\.([0-9]+)_[0-9]{6}_spec\.nc/\1/' | sort -u)

for station in $station_ids; do
    echo "Processing station: ${station}..."
    # list files for the station sorted by date
    files=$(ls $indir/ww3.${station}_*_spec.nc | sort)
    # concat over time
    ncrcat $files -O "$tmpdir/${station}_time.nc"
    # fix values. lon for Alaska, and direction for WWM
    ncap2 -O -s 'longitude=(longitude<0) ? longitude+360 : longitude' "$tmpdir/${station}_time.nc" "$tmpdir/${station}_time.nc"
    ncap2 -O -s "direction=direction+180." "$tmpdir/${station}_time.nc" "$tmpdir/${station}_time.nc"
    ncap2 -O -s "where(direction>=360) direction=direction-360." "$tmpdir/${station}_time.nc" "$tmpdir/${station}_time.nc"

None: This script reads the files, concatenates over time, and then fixes the longitude (only needed for my alaska model, if it is not needed for you, comment out the line) and the mismatch in direction.

• Concatenate by station: under the tmp direction you will have the spectral boundary for each station (concatenated by time and with fixed directions). Now we have to concatenate all these files. I tried to do that using nco, but I kept getting errors, so I suggest using xarray:

import os
import xarray as xr

path = r"/work2/noaa/nos-surge/felicioc/BeringSea/WW3/P03/"
files = os.listdir(path+r"/tmp/")

concat = xr.concat([xr.open_dataset(path+r"/tmp/"+ff) for ff in files],dim='station')
concat.to_netcdf(path+"ww3_all_station.nc")

This is the final file you need to run SCHISM+WWM with WWIII as spectral boundary

Step-by-step guide on how to run SCHISM+WWM with WaveWatchIII Spectral Boundary Forcing:

SCHISM+WWM Inputs:

param.nml:

• Check msc2 and mdc2 and make sure they match what you defined in WWIII
• Set iwbl to be different than zero
• Review the icou_elfe_wwm and set it according to your case
• Set the iof_wwm you need

wwminput.nml:

• &GRID: check MDC and MSC, make sure it matches SCHISM and WWIII
• &GRID: Following are the parameters I had to change, LBCSE = .true., LBCWA = .false., LBCSP = .true., LINHOM = .true., LBSP1D = .false., LBSP2D = .true., LBINTER = .true., IBOUNDFORMAT = 6, FILEWAVE = ww3_all_station.nc, LINDSPRDEG = .true. If you want to use the same setup, check the input files at /work2/noaa/nos-surge/felicioc/BeringSea/R10*

Run the SCHISM+WWM:

• If your run takes longer than the wallclock time limit, checkout this automated hotstart script I wrote for SCHISM+WWM.

[dkobashi]

@felicio93,

Many thanks for detailed instruction of SCHISM-WWM coupling. It is very helpful.

I totally agree with you in that spectral boundary conditions are better than bulk-wave boundary conditions. However, spectral data are often not easily available so my default option would be to use bulk wave parameters despite its caveats. It seems that you suggest running WW3 at basin or global scale and output spectral points along the boundary to run WWM. This may be ok for short-term run, but may not be for long-term hindcast. So, I will start with bulk wave option, but I will also try to run WW3 as well later.

Actually, I was able to run SCHISM-WWM before I saw your instruction, but the performance varied (good for some stations, bad for others). I think I made a mistake for gen_wwmbnd.in. I used gen_wwmbnd.in in the python script, wwmbnd.py without any modifications.

Your gen_wwmbnd.in is

3 #!ncond, i.e. number of conditions
1 2 2 #!cond1 #seg flag ob
2 59 0 #!cond2 #seg flag land
3 0 0

What's the second column for? The first column is just order number (1,2,3) and the third column is for type of OBC as stated in your instruction. But not sure about the second column. Is that the total number of open/land boundaries or boundary nodes?

Also, when I ran the coupled model, the model was really slow even with a 2D SCHISM run (fully coupled). I also tried an uncoupled run (WWM and SCHISM run separately, but in parallel i.e. icou_elfe_wwm = 0) and the computational speed did not change noticeably. It appears that WWM ran slow, not coupling was slow.
Are there ways to improve computational speed? I understand it would depend on multiple factors e.g. time step (which must be same as SCHISM's), spectral segments, physics, I/O etc.

Bests,

-DJ

[feiye-vims]

@felicio93, do you have a link to your tutorial? We can add the link to the WWM part of the SCHISM manual.

[felicio93]

Hi Fei, this tutorial is on our team's private repo. We are considering moving some of the tutorials we have to a public repo at some point. For now, feel free to copy the sections you think are important and add that to the SCHISM manual.

[felicio93]

@dkobashi the second column is the number of boundaries. In my case I have 2 open boundaries in my mesh.