Global mesh for WAVEWATCH3

[Liu-Jincan]

First of all, thank you so much for developing this tool, it is very useful!

Is your feature request related to a problem? Please describe.

I want to create a globalized unstructured grid for WAVEWATCH3 to run, and the code to generate the grid is as follows (using the Example_7 program as a template):

min_el    = 20e3;                % pdf 3.1, minimum resolution in meters.
bbox      = [-180 180;-89 90];   % pdf 2.4, lon min lon max; lat min lat max
max_el    = 120e3;               % pdf 3.1, maximum resolution in meters. 
wl        = 30;                  % pdf 3.4, 30 elements resolve M2 wavelength.
dt        = 0;                   % pdf 4.2, Only reduces res away from coast
grade     = 0.35;                % pdf 4.1, mesh grade in decimal percent. 
R         = 3;                   % pdf 3.3.2, Number of elements to resolve feature.
slp       = 10;                  % pdf 3.5, slope of 10
dem       = 'SRTM15+.nc';        % 
coastline = {'GSHHS_f_L1','GSHHS_f_L6'};
gdat1 = geodata('shp',coastline,'dem',dem,...
                    'bbox',bbox,'h0',min_el);           
fh1 = edgefx('geodata',gdat1,...
            'fs',R,'wl',wl,'max_el',max_el,...
            'slp',slp,'dt',dt,'g',grade); 
mshopts = meshgen('ef',fh1,'bou',gdat1,'nscreen',10,'plot_on',1,'proj','stereo');
mshopts = mshopts.build; 
m = mshopts.grd; % get out the msh object    
m = interp(m,gdat1,'nan','fill','mindepth',1);    

visualization of the m,

When I run write(m,'Global','ww3') , it shows an error in the function writeww3, this is because the variable opedat may be empty, so I modified part of the code of writeww3, (is this a bug?)

if(isempty(opedat))
    fprintf(fid,'%d\n', length(EToV(:,1)));
    m=0;
else
    fprintf(fid,'%d\n', length(EToV(:,1))+opedat.neta);
    m=0;
    for i=1:opedat.nope
        for j=1:opedat.nvdll(i)
            m=m+1;
            fprintf(fid,['%d %s %d %s %d %s %d %s %d %s %d\n'], m,'',15,'',2,'',0,'',0,'',opedat.nbdv(j,i));
        end
    end
end

I use the obtained mesh as input grid for WAVEWATCH3, the final output Netcdf file of WAVEWATCH3 shows the mesh is problematic,

This may be due to the fact that the generated grid surrounds the whole Earth. Steven R. Brus (https://gmd.copernicus.org/articles/14/2917/2021/) also points out the need to modify elements that straddle the −180–+180∘ boundary.

Describe the solution you'd like Is it possible to remove the grid elements around -180 (or 180)?

Describe alternatives you've considered

• When the grid range is determined as [-179.9,179.9;-89.9,89.9], a similar situation as above is encountered.
• When the grid range is determined as [-179,179;-78,78], the final Netcdf shows the mesh is good.
• when the grid range is determined as [-179,179;-89,89], the Mediterranean Sea is not taken into account.
• Replace the projection from stereo to Robinson. Due to the following code in meshgen.m, you need to set h0 to very small for global meshes (stereo does not need to set h0 to very small to get global meshes).

  % kjr Oct 2018 use outer coarsest box for
  % multiscale meshing
  lon_mi = obj.bbox{1}(1,1)-obj.h0(1)/1110;
  lon_ma = obj.bbox{1}(1,2)+obj.h0(1)/1110;
  lat_mi = obj.bbox{1}(2,1)-obj.h0(1)/1110;
  lat_ma = obj.bbox{1}(2,2)+obj.h0(1)/1110;

Additional context This is the first time I submit an issue on github, there may be problems with the format of the question, and I'm looking forward to getting a reply, thanks~~

[krober10nd]

Hello @Liu-Jincan, I'll try to answer all these questions.

1. There are no open boundaries in a truly global model since there's no ocean boundary to force at.
2. The elements are indeed modified to straddle the meridian. Following Steven's paper are you additionally using the modified version of the WW3 code he posted? Additional changes were necessary in the source code of WW3.
3. For global meshes, the stereo option has been extensively tested with both ADCIRC and less so with WW3. I would stick with that option for no and focus on 2).

[krober10nd]

@sbrus89 any ideas what may be wrong? I recall you used global meshes created with OceanMesh at one point to force WW3 models.

[krober10nd]

Oh and in regard to the Mediterranean disappearing please keep in mind you can change the parameter dj_cutoff in the msh.clean routine which is automatically called during generation. See the docstring for setting that value.

[sbrus89]

I also started with the Example 7 script. One that I used successfully can be found here: https://github.com/sbrus89/ww3_utils/blob/master/grid/WW3_global_constant.m

I don't think I knew about the writeww3 method at that time, so I used a separate script to convert from the fort.14 format to the ww3 grid format. It looks like writeww3 should do the exact same thing.

@krober10nd, is correct that there were -180/180 modifications necessary in the WW3 code. However, the mesh file doesn't require any modification. So be sure to use a recent version of WW3. The global unstructured changes were merged in this PR: https://github.com/NOAA-EMC/WW3/pull/335 in March 2021.

[Liu-Jincan]

@krober10nd,

Oh and in regard to the Mediterranean disappearing please keep in mind you can change the parameter dj_cutoff in the msh.clean routine which is automatically called during generation. See the docstring for setting that value.

The default value of dj_cutoff is 0.25, when set to 0.4, the Mediterranean Sea is successfully displayed.

[Liu-Jincan]

@krober10nd,

2. The elements are indeed modified to straddle the meridian. Following Steven's paper are you additionally using the modified version of the WW3 code he posted? Additional changes were necessary in the source code of WW3.

At first, I was using WW3 version 6.07, and when preprocessing the mesh, I would get the STOP WRONG TRIANGLE error. Then I tried to use the modified version of WW3 provided by Steven, but I ended up with this error too (probably because there were some settings I didn't notice). Finally I tried to use the WW3 version downloaded on May 6, 2022, which successfully preprocesses the mesh, as @sbrus89 mentioned above.