Set a land mask on the sea in SeaGrid

Message

Scarlett · #1 Unread post by **Scarlett** » Wed Feb 07, 2018 5:15 pm

Hi ROMS users!
I'm working with SeaGrid, in my map I want to mask a part of the ocean, like it was a land mask. Already I'm doing it clicking on a each ocean cell where that I want do the mask, but I have many cells (because my resolution is 1/24°) and my masking will be very slow.
I wonder if exist an other method to do the mask in a specific area, if I don't want to click cell by cell?
Somebody know how? Thanks in advance.
I add my grid, and a select the part that i want to mask.

Mar.Mo.

kate · #2 Unread post by **kate** » Wed Feb 07, 2018 6:20 pm

Hmm, first, the model requires all angles to be right angles. From what I know of seagrid, your orthogonality errors could be substantial with that grid because you don't ask it to converge.

The matlab editmask has the capability of changing whole rectangles at a time.

I do this stuff outside matlab. Just figure out the patch I want to change and code it up. This is Python, but you can use your language of comfort (you have one, right???):

Code: Select all

import numpy as np
import netCDF4
import sys

ncfile = 'grid_Arctic_2.nc'
nc = netCDF4.Dataset(ncfile, 'a', format='NETCDF3_CLASSIC')

mask_rho[828:980,640:] = 0
mask_rho[837:932,571:] = 0
mask_rho[932:1011,658:] = 0

nc.variables['mask_rho'][:] = mask_rho

nc.close()

Scarlett · #3 Unread post by **Scarlett** » Sat Feb 10, 2018 1:04 am

Hi Kate:
So many thanks for your help, I had not known that the model requires all angles to be right angles. I already corrected this, for the other hand I could mask the sea in matlab, based on your code. Thanks.

But I have a question, Do you know how can I save the netcdf file with the new variable (already edited)?
Thanks in advance.

Mar.Mo.

This is my code:

%% Load Netcdf file
ncfile ='/home/ahumada/roms-3.7/matlab/seagrid/GT4_roms_grd.nc'

nc=netcdf(ncfile)
ncdump(nc)

% extract variable

latr=nc{'lat_rho'}(:,:);
lonr=nc{'lon_rho'}(:,:);
mask_rho=nc{'mask_rho'}(:);
h=nc{'h'}(:,:);

%% mask

% nombre de la variable(renglon i:renglon n, columna i: columna n)= O
% cero es tierra

mask_rho(387:432,203:347) = 0
mask_rho(327:432,386:720) = 0
mask_rho(139:327,641:720) = 0
mask_rho(173:326,534:640) = 0
mask_rho(161:180,543:557) = 0
mask_rho(161:174,558:597) = 0
mask_rho(182:326,506:533) = 0
mask_rho(318:326,500:505) = 0

latr=nc{'lat_rho'}(:,:);
lonr=nc{'lon_rho'}(:,:);
h=nc{'h'}(:,:);
h=h.*mask_rho;

figure(1)
pcolor(lonr, latr,h);
colormap('jet');
hold on
shading interp

c=colorbar('SouthOutside');

shchepet · #4 Unread post by **shchepet** » Sun Feb 11, 2018 11:27 am

Is this the kind of grid you want to generate?

shchepet · #5 Unread post by **shchepet** » Sun Feb 11, 2018 11:53 am

The procedure for generating ROMS land mask from USGS coastline data described here:
viewtopic.php?f=23&t=3878&p=14908

Scarlett · #6 Unread post by **Scarlett** » Wed Feb 14, 2018 5:46 pm

Hi Alexander F. Shchepetkin:
So many thanks for your help. Yes, it is the kind of grid that I want to create.
I followed your advice to do it with the new tools. But I had have some problems with the mpc compile.

I followed the README.compiling file, but when I compile mpc I get:

mpc.f:335:18: Error: Syntax error in IF-expression at (1)
mpc.f:336:19: Error: Syntax error in IF-expression at (1)
mpc.f:337:20: Error: Syntax error in IF-expression at (1)
mpc.f:338:21: Error: Syntax error in IF-expression at (1)
mpc.f:339:22: Error: Syntax error in IF-expression at (1)
mpc.f:340:23: Error: Syntax error in IF-expression at (1)
mpc.f:343:19:

endif
1
Error: Expecting END PROGRAM statement at (1)
mpc.f:344:18:

endif
1
Error: Expecting END PROGRAM statement at (1)
mpc.f:345:17:

endif
1
Error: Expecting END PROGRAM statement at (1)
mpc.f:346:16:

endif
1
Error: Expecting END PROGRAM statement at (1)
mpc.f:347:15:

endif ! After this moment "i" is
1
Error: Expecting END PROGRAM statement at (1)
mpc.f:348:14:

endif ! index of the last character
1
Error: Expecting END PROGRAM statement at (1)
mpc.f:349:13:

endif ! of Fortran type declaration
1
Error: Expecting END PROGRAM statement at (1)
mpc.f:350:12:

endif ! which may be either "real",
1
Error: Expecting END PROGRAM statement at (1)
mpc.f:351:11:

endif ! "integer", or "character"..
1
Error: Expecting END PROGRAM statement at (1)
mpc.f:355:10:

do while (j < length .and. str(j)==' ')
1
Error: Unclassifiable statement at (1)
mpc.f:357:13:

enddo ! first non-blanc str
1
Error: Expecting END IF statement at (1)
mpc.f:374:17: Error: Syntax error in IF-expression at (1)
mpc.f:377:12:

do while( str(is) == ' ' .and. is < length)
1
Error: Unclassifiable statement at (1)
mpc.f:379:15:

enddo
1
Error: Expecting END IF statement at (1)
mpc.f:381:12:

do while( str(ie+1) /= ' ' .and. str(ie+1) /= ','
1
Error: Unclassifiable statement at (1)
mpc.f:384:15:

enddo
1
Error: Expecting END IF statement at (1)
mpc.f:386:19: Error: Syntax error in IF-expression at (1)
mpc.f:391:14:

bffr(k)=str(k+is-1)
1
Error: Unclassifiable statement at (1)
mpc.f:397:16:

str(k+isft)=str(k) ! "isft", then move the rest
1
Error: Unclassifiable statement at (1)
mpc.f:400:14:

str(i+1)='(' ! insertion, adjust the length
1
Error: Unclassifiable statement at (1)
mpc.f:401:14:

str(i+2)='l' ! of the line accordingly,
1
Error: Unclassifiable statement at (1)
mpc.f:402:14:

str(i+3)='e' ! then insert the new-style
1
Error: Unclassifiable statement at (1)
mpc.f:403:14:

str(i+4)='n' ! size specification.
1
Error: Unclassifiable statement at (1)
mpc.f:404:14:

str(i+5)='='
1
Error: Unclassifiable statement at (1)
mpc.f:406:16:

str(i+k+5)=bffr(k)
1
Error: Unclassifiable statement at (1)
mpc.f:408:14:

str(i+ie-is+7)=')'
1
Error: Unclassifiable statement at (1)
mpc.f:416:16:

str(k+isft)=str(k)
1
Error: Unclassifiable statement at (1)
mpc.f:420:14:

str(i+1)='('
1
Error: Unclassifiable statement at (1)
mpc.f:421:14:

str(i+2)='k'
1
Error: Unclassifiable statement at (1)
mpc.f:422:14:

str(i+3)='i'
1
Error: Unclassifiable statement at (1)
mpc.f:423:14:

str(i+4)='n'
1
Error: Unclassifiable statement at (1)
mpc.f:424:14:

str(i+5)='d'
1
Error: Unclassifiable statement at (1)
mpc.f:425:14:

str(i+6)='='
1
Error: Unclassifiable statement at (1)
mpc.f:427:16:

str(i+k+6)=bffr(k)
1
Error: Unclassifiable statement at (1)
mpc.f:429:14:

str(i+ie-is+8)=')'
1
Error: Unclassifiable statement at (1)
mpc.f:459:16:

elseif ( str(j) /= '(' .and. ( str(j) == ' ' .or.
1
Error: Unexpected junk after ELSE statement at (1)
mpc.f:490:16:

str(k+isft)=str(k) ! isft symbols to the right
1
Error: Unclassifiable statement at (1)
mpc.f:493:16:

str(k)=' ' ! and increase length of
1
Error: Unclassifiable statement at (1)
mpc.f:498:14:

str(i+1)='('
1
Error: Unclassifiable statement at (1)
mpc.f:499:14:

str(i+2)='k'
1
Error: Unclassifiable statement at (1)
mpc.f:500:14:

str(i+3)='i'
1
Error: Unclassifiable statement at (1)
mpc.f:501:14:

str(i+4)='n'
1
Error: Unclassifiable statement at (1)
mpc.f:502:14:

str(i+5)='d'
1
Error: Unclassifiable statement at (1)
mpc.f:503:14:

str(i+6)='='
1
Error: Unclassifiable statement at (1)
mpc.f:504:14:

str(i+7)=type
1
Error: Unclassifiable statement at (1)
mpc.f:505:14:

str(i+8)=')'
1
Error: Unclassifiable statement at (1)
mpc.f:507:13:

endif
1
Error: Expecting END PROGRAM statement at (1)
mpc.f:508:11:

endif
1
Error: Expecting END PROGRAM statement at (1)
mpc.f:613:17: Error: Syntax error in IF-expression at (1)
mpc.f:615:16:

elseif (str(i) == '.' .and. lswtch) then
1
Error: Unexpected junk after ELSE statement at (1)
mpc.f:618:13:

endif
1
Error: Expecting END DO statement at (1)
mpc.f:626:19: Error: Syntax error in IF-expression at (1)
mpc.f:628:18:

elseif (str(m) /= ' ') then
1
Error: Unexpected junk after ELSE statement at (1)
mpc.f:630:15:

endif
1
Error: Expecting END DO statement at (1)
mpc.f:633:30:

if ( lswtch .or. str(m) < 'A' .or. ! Cond.(3)
1
Error: Syntax error in IF-expression at (1)
mpc.f:641:21: Error: Syntax error in IF-expression at (1)
mpc.f:643:20:

elseif (str(m) /= ' ') then
1
Error: Unexpected junk after ELSE statement at (1)
mpc.f:645:17:

endif
1
Error: Expecting END DO statement at (1)
mpc.f:650:21: Error: Syntax error in IF-expression at (1)
mpc.f:653:16:

do while (str(i)==' ' .and. i<length)
1
Error: Unclassifiable statement at (1)
mpc.f:655:19:

enddo
1
Error: Expecting END IF statement at (1)
mpc.f:656:24: Error: Syntax error in IF-expression at (1)
mpc.f:660:20:

elseif (str(m) == '_') then
1
Error: Unexpected junk after ELSE statement at (1)
mpc.f:670:23: Error: Syntax error in IF-expression at (1)
mpc.f:671:25: Error: Syntax error in IF-expression at (1)
mpc.f:672:20:

str(m )='D'
1
Error: Unclassifiable statement at (1)
mpc.f:673:20:

str(m+1)='0'
1
Error: Unclassifiable statement at (1)
mpc.f:676:22:

str(i)=str(i+1)
1
Error: Unclassifiable statement at (1)
mpc.f:678:24:

elseif (str(m+1) == 'e') then
1
Error: Unexpected junk after ELSE statement at (1)
mpc.f:679:20:

str(m)='D'
1
Error: Unclassifiable statement at (1)
mpc.f:682:22:

str(i)=str(i+2)
1
Error: Unclassifiable statement at (1)
mpc.f:685:19:

endif
1
Error: Expecting END DO statement at (1)
mpc.f:686:20:

elseif ( str(m) < 'A' .or. ( str(m) > 'Z'. and.
1
Error: Unexpected junk after ELSE statement at (1)
mpc.f:689:16:

do while(str(m-1) == ' ')
1
Error: Unclassifiable statement at (1)
mpc.f:693:18:

str(i+2)=str(i)
1
Error: Unclassifiable statement at (1)
mpc.f:695:16:

str(m)='_'
1
Error: Unclassifiable statement at (1)
mpc.f:696:16:

str(m+1)='8'
1
Error: Unclassifiable statement at (1)
mpc.f:698:17:

endif
1
Error: Expecting END PROGRAM statement at (1)
mpc.f:699:15:

endif
1
Error: Expecting END PROGRAM statement at (1)
mpc.f:700:13:

endif
1
Error: Expecting END PROGRAM statement at (1)
mpc.f:701:11:

enddo
1
Error: Expecting END PROGRAM statement at (1)
mpc.f:706:17: Error: Syntax error in IF-expression at (1)
mpc.f:707:39:

scratch(i:i)=char(ichar(str(i))+case_fold)
1
Error: Syntax error in argument list at (1)
mpc.f:708:72: Error: Unexpected ELSE statement at (1)
mpc.f:709:12:

scratch(i:i)=str(i)
1
Error: Unclassifiable statement at (1)
mpc.f:710:13:

endif
1
Error: Expecting END DO statement at (1)
mpc.f:794:15: Error: Syntax error in IF-expression at (1)
mpc.f:796:10:

do while (str(i) == ' ' .and. i < length)
1
Error: Unclassifiable statement at (1)
mpc.f:798:13:

enddo
1
Error: Expecting END PROGRAM statement at (1)
mpc.f:799:17: Error: Syntax error in IF-expression at (1)
mpc.f:802:12:

do while (str(i) == ' ' .and. i < length)
1
Error: Unclassifiable statement at (1)
mpc.f:804:15:

enddo
1
Error: Expecting END PROGRAM statement at (1)
mpc.f:805:19: Error: Syntax error in IF-expression at (1)
mpc.f:807:14:

do while (str(i) == ' ' .and. i < length)
1
Error: Unclassifiable statement at (1)
mpc.f:809:17:

enddo
1
Error: Expecting END PROGRAM statement at (1)
mpc.f:811:21: Error: Syntax error in IF-expression at (1)
mpc.f:820:18:

str(i+j-1)=scratch(j:j)
1
Error: Unclassifiable statement at (1)
mpc.f:823:17:

endif
1
Error: Expecting END PROGRAM statement at (1)
mpc.f:824:15:

endif
1
Error: Expecting END PROGRAM statement at (1)
mpc.f:825:13:

endif
1
Error: Expecting END PROGRAM statement at (1)
mpc.f:826:11:

endif
1
Error: Expecting END PROGRAM statement at (1)
mpc.f:842:37: Error: Expected a right parenthesis in expression at (1)
mpc.f:844:37: Error: Expected a right parenthesis in expression at (1)
mpc.f:855:19: Error: Syntax error in IF-expression at (1)
mpc.f:857:18:

elseif (str(k) == ',') then
1
Error: Unexpected junk after ELSE statement at (1)
mpc.f:859:18:

elseif (str(k) == '=') then
1
Error: Unexpected junk after ELSE statement at (1)
mpc.f:861:18:

elseif (str(k) == '/') then
1
Error: Unexpected junk after ELSE statement at (1)
mpc.f:863:18:

elseif (str(k) == '+' .or.
1
Error: Unexpected junk after ELSE statement at (1)
mpc.f:867:15:

endif
1
Error: Expecting END DO statement at (1)
mpc.f:909:37: Error: Expected a right parenthesis in expression at (1)
mpc.f:911:38: Error: Expected a right parenthesis in expression at (1)
mpc.f:914:12:

str(i)=' '
1
Error: Unclassifiable statement at (1)
mpc.f:916:10:

str(6)='&'
1
Error: Unclassifiable statement at (1)
mpc.f:919:12:

str(i-m)=str(i)
1
Error: Unclassifiable statement at (1)
mpc.f:922:37: Error: Expected a right parenthesis in expression at (1)
mpc.f:924:37: Error: Expected a right parenthesis in expression at (1)
make: *** [mpc] Error 1

Why do I get this error?
Thanks in advance for your help.

Mar.Mo.

shchepet · #7 Unread post by **shchepet** » Wed Feb 14, 2018 11:25 pm

The reason why you have these compiling problems is a faulty version of mpc.F supplied
with the package at some point in the past. This was corrected, but it appears that the
old version came back because at some point the web page was restored from the backup
after disk failure and I did not realize that this happened. At this moment I do not
have access to the web server to replace tar file, so use file "tools.tar" attached to
this response.

Note that a Matlab grid generation tool which was used to make the grid in the previous
post is also available from the same web suite
http://people.atmos.ucla.edu/alex/ROMS/grid_gui.tar
although it is relatively new and lacks instructions - a README file is being written,
but not included into the package yet. It is basically a tool to create analytical
grid using staged conformal transformations. Setting all weird parameters to zero

cent_lat = EWtpr = NStpr = cushn1 = cushn2 = cushn3 = cushn4 = 0

reverts it to easy_grid type of tool: it creates a patch of Mercator grid transferred
to user-specified location and rotated by user-specified azimuth angle.

Optionally the initial rectangular patch on Mercator plane can be deformed into
a desired shape by two kind of transforms:

(i) "tapering" in either direction - rectangle turns into a sector of polar coordinates
with simultaneous exponential stretching along the radial direction (delta r is
proportional to r), and

(ii) "cushioning" - rectangle becomes convex in one direction, concave in the other
(or vice versa, depending whether parameter cushn is positive or negative).
This is actually transform into elliptic coordinates: imagine a family of confocal ellipses
and a family of hyperbolae with the same foci. It can be proven that ellipses and
hyperbolae intersect each other at right angle. This is the basis.

Parameter cushn controls the location of foci:

cushn=0 means that they are infinitely
far away - rectangle remains rectangle and nothing happens;

cushn > 0 makes them approach from infinity along x-axis - actually places the foci
on the horizontal line passing through point (xctr,yctr) and at equal distance from
the point, so the rectangle become concave in eastern and western side (they become
hyperbolae) and convex on northen and southern (they become cutoffs of ellipses).
Now exactly depends on (xctr,yctr).

cushn < 0 along y-axis -- similarly, but everything is rotated
90 degrees.

It is a bit tricky to control, but after some practice you can fit very weird shapes.

Couple tips:

1. The boundaries of the map are determined automatically based on current grid
configuration. Once the desired geographical region is captured unclick the "redraw map"
box on top -- this would speed up the drawing by a factor of 100 LITERALLY(!) because
map does not need to be redrawn every time when you hit "update" button.

2. Leaving one of the four boxes, nx, ny, size_x, size_y, blank makes the tool to calculate
the missing value from the condition of isotropy of the grid: dx = dy EVERYWHERE (same
as pm == pn EVERYWHERE). This is very useful. Of course, if size_x, size_y are both
specified, but one of nx,ny is not, then isotropy is not exact because of rounding
to the nearest integer; but is nx,ny are both specified, while one of size_x,size_y left
blank, then isotropy is perfect. THIS IS HOW I ALWAYS RECOMMEND to do it at the end:
at first specify both size_x,size_y, and leave nx blank, play with setting until
you like it, then see Matlab command window to figure nx (the tool always prints the
computed value), put nx into the box, but erase size_x or size_y -- the tool will
adjust its value to make the grid be perfectly isotropic.

3. UPDATE button draws the new contour on top of the old (alternating between red and
magenta colors). This is useful to visualize incremental change, so you can accept or
reverse it. Hitting UPDATE button second time without changing any of the parameters
makes single contour.

4. The contours are actually double lines exactly 1dx or 1dy apart from each other.
this is useful because one can use Matlab zoom facility to check how the edge of
the grid is placed relative to the coastline: ideally the coastline should be in
between the double lines - in this case there will be EXACTLY 1 row of mask points.

5. BREXIT means "Break and Exit" this button just quits the tool without doing anything
else, i.e., it does not trigger saving into the file.

shchepet · #8 Unread post by **shchepet** » Wed Feb 14, 2018 11:35 pm

you can fit weird shapes after some practice...

: Black Sea grid.

Scarlett · #9 Unread post by **Scarlett** » Thu Feb 15, 2018 7:21 pm

Dear Alexander:
Thanks a lot for your help, it was and will be so helpful.
Already with the new tools, I could compile mpc. But now when I type make, I get:
--------------------------------------------------------------
[ahumada@master tools]$ make
gfortran -fopenmp -fno-second-underscore -c -I/usr/local/lib -O3 def_clim_file.f
def_clim_file.f:538:40:

ierr=nf_create(trgfile, nf_netcdf4, nctarg)
1
Error: Symbol ‘nf_netcdf4’ at (1) has no IMPLICIT type
make: *** [def_clim_file.o] Error 1
--------------------------------------------------------------
I don't konow if is cause by my netcdf libreries or if is def_clim_file.f file problem.
How can I solve this problem?. Thanks in advance.
Best regards.

Scarlett

shchepet · #10 Unread post by **shchepet** » Fri Feb 16, 2018 9:24 am

Error: Symbol ‘nf_netcdf4’ at (1) has no IMPLICIT type

This one you should be able to overcome yourself. At first, my package comes with several
README.something files: just look through them. Among them is README.compiling which
tells what to do.

Then, it is quite obvious that your netcdf.inc file does not contain definitions for netCDF 4
functions, so you can either

(i) preferred way, install and use netCDF library with full netDF4/hdf5 capabilities, or

(ii) handicapped way (ok as a temporary fix): exclude all functionality from the package
which requires netCDF4 - basically keep trying to compile everything by "make", see which
ones do not compile with errors line this, and exclude them from INST_LIST in Makefile.

Some people find compiling netCDF 4 with hdf5 and zlib too tedious. In reality it takes
about 20 minutes, using mostly computer mouse for copy-pasting and almost without
ever touching keyboard. See README.compiling_sequence inside the attached
netcdf4_help.tar.

Scarlett · #11 Unread post by **Scarlett** » Mon May 28, 2018 4:21 pm

Dear Alexander:
I tried install this mask tool, but was necessary change some of my libraries and dependencies, and this could affect the dependencies that ROMS uses, I chose not install. For a time I decided work with the normal grid (without mask the Atlantic) but is demanding in time way.

So what I did was create a grid file with a mask in the Atlantic from which I created in seagrid. The input ROMS files were creating without problems. But ROMS do not work with this new grid.

Could you tell me why ROMS do not work with this new grid?
Some suggestion? What can I do?

So many thanks for your time and suggestion.
-
Scarlett Mar.Mo.

wilkin · #12 Unread post by **wilkin** » Tue May 29, 2018 8:26 am

Try the tool editmask which is in the myroms/landmask folder of Matlab tools

shchepet · Tue May 29, 2018 12:05 pm

Try the tool editmask which is in the myroms/landmask folder of Matlab tools

John,

Here is the problem: At first many thanks to Andrey Scherbina for creating editmask tool back in 1999.
I checked editmask in myroms/landmask today by getting it from via svn,
svn checkout --username shchepet https://www.myroms.org/svn/src/matlab matlab
and found that despite updated svn $Id: editmask.m 895 2018-02-11 23:15:37Z arango $
and $ Copyright (c) 2002-2018 The ROMS/TOMS Group, the tool itself is pretty old and
needs to be updated in many respects. For example, editmask.m calls read_mask which uses modern
Matlab native netCDF support, but then writes it back via write_mask which still relies on mexnc.

Then this particular version of editmask lacks native support of GSHHS data directly from editmask - hence
user must provide his/her own coastlines - this is a very very old fashion way of doing this (editmask
AGRIF/UCLA use gshhs data for 10+ years already, thanks to m_map package from Rich Pawlowicz it
is done natively - there is no need to convert data into *.mat files). In Rutgers matlab scripts there are two
places where GSHHS data is handled landmask/landsea.m (generates initial version
of mask and puts it into ROMS grid netCDF file) and coastlines/extract_coast.m which
simply reads the desired GSHHS dataset and converts it into *.mat or *.cst file to be used later.

Generating land mask via landsea ---> r_gshhs ---> inpolygon procedure takes considerable
time, and I doubt that it is even practical for more or less reasonable large grid - it may take hours or
even days. Much better alternatives are available for several years already.

Appearance of the tool on the screen is too poor not optimized to use computer screen in any ergonomic
manner, especially given the fact that in nearly 20 years if existence of editmask resolution (number
of pixels) of the monitors changed dramatically, and so do dimensions of ROMS grids.

I maintain my own version of editmask which is designed to work in concert with automatic mask
generation, hence minimizing the need and time consumed by manual editing,

http://people.atmos.ucla.edu/alex/ROMS/grid_gui.tar

http://people.atmos.ucla.edu/alex/ROMS/tools.tar

http://people.atmos.ucla.edu/alex/ROMS/editmask.tar

however my tools are incompatible with Rutgers ROMS for various legacy reasons.

Mainly I do not store meaningless and/or obsolete variables in the grid file.

What is the point of storing domain sizes xl,el if grid is spherical?

Why any one needs x_rho,y_rho, x_psi,y_psi,x_y,y_u, x_v,y_v ?

Why do you need dndx,dmde, if ROMS computes them internally any way?

As for mask arrays, I keep only mask_rho, and do not store three others: they
are computed from mask_rho by model itself in both UCLA and AGRIF/CROCO cases,
besides matlab script landsea.m

Code: Select all

mask_psi(1:Jm-1,1:Im-1)=mask_rho(2:Jm  ,2:Im  ).*                       ...
                        mask_rho(2:Jm  ,1:Im-1).*                       ...
                        mask_rho(1:Jm-1,2:Im  ).*                       ...
                        mask_rho(1:Jm-1,1:Im-1);

does it in not the best way, see viewtopic.php?f=19&t=471 for the discussion.

I no longer need Coriolis frequency "f" in grid netCSF file, since it can be computed from sin(lat_rho), and besides
occasionally I need horizontal Coriolis parameters, f_h ~ cos(lat_rho), so I decided not to add extra variables, but
simply calculate all three of them. (thought "f" is still stored).

Rutgers ROMS users bump into these silly problems, and usually give up without looking into more detail, complaining
about that "model does not run" or something of this sort.

Is it so hard to modify Rutgers ROMS in such a way that it no longer attempts to read mask_u,
mask_v, mask_psi, dmde, and dndx, but rather calculates them internally?

Specifically, to rely only on what is only mathematically necessary?

Here is example:

Code: Select all

netcdf GOM1k1_grd {
dimensions:
        xi_rho = 1652 ;
        eta_rho = 1190 ;
        xi_u = 1651 ;
        eta_v = 1189 ;
variables:
        char spherical ;
                spherical:long_name = "Grid type logical switch" ;
                spherical:option_T = "spherical" ;
        double lon_rho(eta_rho, xi_rho) ;
                lon_rho:long_name = "longitude of RHO-points" ;
                lon_rho:units = "degree East" ;
        double lat_rho(eta_rho, xi_rho) ;
                lat_rho:long_name = "latitude of RHO-points" ;
                lat_rho:units = "degree North" ;
        double lon_psi(eta_v, xi_u) ;
                lon_psi:long_name = "longitude of PSI-points" ;
                lon_psi:units = "degree East" ;
        double lat_psi(eta_v, xi_u) ;
                lat_psi:long_name = "latitude of PSI-points" ;
                lat_psi:units = "degree North" ;
        double pm(eta_rho, xi_rho) ;
                pm:long_name = "curvilinear coordinate metric in XI-direction" ;
                pm:units = "meter-1" ;
        double pn(eta_rho, xi_rho) ;
                pn:long_name = "curvilinear coordinate metric in ETA-direction" ;
                pn:units = "meter-1" ;
        double angle(eta_rho, xi_rho) ;
                angle:long_name = "angle between east and XI-directions" ;
                angle:units = "degrees" ;
        double f(eta_rho, xi_rho) ;
                f:long_name = "Coriolis parameter at RHO-points" ;
                f:units = "second-1" ;
        double hraw(eta_rho, xi_rho) ;
                hraw:long_name = "raw bathymetry at RHO-points" ;
                hraw:units = "meter" ;
                hraw:topo_data_file = "/home/alex/srtm30_plus/" ;
                hraw:smoothing_radius = 4. ;
        double hsmth(eta_rho, xi_rho) ;
                hsmth:long_name = "model bathymetry at RHO-points" ;
                hsmth:units = "meter" ;
                hsmth:hmin = 5. ;
                hsmth:hmax = 4950. ;
                hsmth:r_max = 0.15 ;
                hsmth:method = "mask-dependent rx-conditioned log-smoothing" ;
                hsmth:iters_cond = 1000 ;
        double mask_rho(eta_rho, xi_rho) ;
                mask_rho:long_name = "land/sea mask at RHO-points" ;
                mask_rho:units = "land/water (0/1)" ;
        double h(eta_rho, xi_rho) ;
                h:long_name = "smooth topography merged with parent grid at open boundaries" ;
                h:parent_grid = "GOM5k_grd.nc" ;
                h:merging_flags = "ES" ;
                h:merging_width = 64 ;
        double wgt(eta_rho, xi_rho) ;
                wgt:long_name = "parent-to-child topography merging weights" ;

// global attributes:
                :Title = "ROMS grid produced by Easy Grid" ;
                :Settings = "nx=1650 ny=1188 size_x=1817.0729 size_y=1306 cent_lat=0 tapering=0 Lat=24.665 Lon=-89.23 rotate=7.5 flip_xy=0" ;
                :Date = "06-May-2017" ;
}

where "hraw" and "hsmth" are needed only by tools processing topography, but are not used by ROMS code. The reason for keeping them
is because if topography needs to be "re-smoothed" or "re-merged" with parent grid in can be done stage-by-stage rather than starting
from the very beginning.

Variable "wgt" is the merging function, h = (1-wgt)*hsmth + wgt*h_parent, where wgt=0 everywhere inside the domain, except some areas
(bands) near the open boundaries; h_parent is parent-grid topography interpolated to child grid (usually by splines, see r2r_match_topo.F
from tools.tar). In the case of complicated coastline - open boundaries are partially blocked by narrow strips of land, its setting is may
be nontrivial, see below. Overall "wgt" is/may be used by ROMS model to control sponge layers of enhanced viscosity.

son-goku · #14 Unread post by **son-goku** » Fri Nov 09, 2018 4:25 pm

Hi Alexander F. Shchepetkin:

I try to use the UCLA grid, i downloaded tools.tar from http://people.atmos.ucla.edu/alex/ROMS/tools.tar
I installed all the prerequisites on my machine and followed the steps of the README.compiling file but in step "d" I write make in the terminal
there is an error.

&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
vort.f

70:

C$OMP PARALLEL SHARED(Lm,Mm, pm,pn, f,f_q, fnd_rmask,rmask, dm_u,dn_v,
1
Error: Syntax error in OpenMP variable list at (1)
vort.f

6:

& iA_q)
1
Error: Bad continuation line at (1)
vort.f

6:

& iA_q)
1
Error: Unclassifiable statement at (1)
vort.f

72: Error: Unexpected !$OMP END PARALLEL statement at (1)

&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&

Annex more information in error.zip

I hope you can help me, thank you.

shchepet · #15 Unread post by **shchepet** » Fri Nov 09, 2018 6:23 pm

Obviously this should be

Code: Select all

C$OMP PARALLEL SHARED(Lm,Mm, pm,pn, f,f_q, fnd_rmask,rmask, dm_u,dn_v,
C$OMP&                                                           iA_q)
        call init_vort_thread(Lm,Mm, pm,pn, f,f_q, fnd_rmask,rmask,
     &                                                dm_u,dn_v, iA_q)
C$OMP END PARALLEL

i.e., change "& into "C$OMP&" in continuation line for OpenMP directive
in the original vort.F file.

...interestingly Intel Ifort compiler does not care.

son-goku · #16 Unread post by **son-goku** » Tue Nov 13, 2018 9:17 pm

Thank you so much, I have solved the problem. I think it's because I use Gfortran.

shchepet · Wed Nov 14, 2018 11:23 am

No, you should be able to compile the package using any compiler.

I updated the package (see attachment at the end of this post) specifically paying attention to
the warnings issued by GCC compiler, which was forced into checking conformance to Fortran 2008
standards in pedantic mode. No bugs were discovered, and Gfortran tends to complain excessively
that a variable can be uninitialized where, in fact, I can mathematically prove that it is impossible.

It did, however issued warnings about overloading function standard violation in a couple places:

Code: Select all

integer(kind=2) :: mss(nx,ny), mrg(nx,ny)
....
do j=...
  do i=...
       mrg(i,j)=max(mss(i,j), 1-mss(i,j+1))
  enddo
enddo

where,strictly speaking, in the above expression 1-mss(i,j+1) is of default integer type, that
is (kind=4). Here constant 1 is of(kind=4) by default, so mss(i,j+1) is promoted from kind=2
to kind=4 and subtracted from 1, resulting in a kind=4 second argument of max.
However all arguments of max should be of the same kind, hence warning message.

to get rid of warning message one needs to make 1 be of kind=2 as well

Code: Select all

integer(kind=2) :: mss(nx,ny), mrg(nx,ny)
integer(kind=2), parameter :: ione=1
....
do j=...
  do i=...
       mrg(i,j)=max(mss(i,j), ione-mss(i,j+1))
  enddo
enddo

...it is as silly as this.

son-goku · #18 Unread post by **son-goku** » Wed Nov 21, 2018 5:05 pm

thanks you, I can compile correctly. Can you help me with a grid?

I want to make a grid with sides A, B, C completely right, similar to the attached image

shchepet · Thu Nov 22, 2018 10:21 am

The answer in "no", all three of them, A,B,C, cannot be made straight at the same time.

This is the nature of conformal mapping: once A and C are straight, there is no way to fit a
curve which would go along the coast, and be perpendicular to A and C at the junction points.

B can me made straight easily. but in this case A and C will be curved.

A can be made straight, but then B and C be curved, and the remaining line on the top would
have to enter deeply into the land and get a lot of masking.

Making C straight? I see no point as it is inside the mask any way.

Besides, there is no wisdom even trying to the make side boundaries straight? Even if some
of them are straight, once you depart from the boundary go into the domain, the coordinate
lines are no longer straight any way.

It has to be recognized that this analytical curvilinear grid is optimal by the design in sense
tat it is:

(1) exactly othrogonal;

(2) perfectly isotropic, pm=pn at every grid cell (same as local dx=dy at every cell, despite
that fact that both dx, dy change from cell to cell;

(3) with some diligence has minimal possible non-uniformity of the resolution: the ratio
min(dx)/max(dx) kept as small as possible for a given situation of the particular shape
of the domain and coastline;

The fundamental flaw of seagrid, as well as the more recent package GridBuilder which
replaces it, is that they provide no means to enforce perpenducular junction of the side
boundaries. Both packages allow you to define points at the side and move them in and
out to achieve the desired shape, but the junction points end up be not perpendicular,
and it is basically left up to the user to make them as close to perpendicular as possible --
it is possible in theory, but is very hard to achieve in practice. This makes it very hard to
construct a satisfactory grid: everybody complains about corner problems.

This also leads to the very long-standing confusion in ROMS community (and in POM
community before) that one can fit curvilinear grid to virtually any shape -- as for the
contour yes, you can fit it, but if the corner is not exactly 90 degrees, once you start
filling up the grid points by solving Dirichet elliptic equation, the grid points tend to
either concentrate toward the corner (if >90 degrees) or, conversely, avoid it (< 90).

son-goku · #20 Unread post by **son-goku** » Tue Nov 27, 2018 1:38 am

Thank you

I understand.

I made a useful grid with this tool.

jivica · #21 Unread post by **jivica** » Tue Nov 27, 2018 7:08 am

Dear Sasha,

really like you answers, always getting to the bottom of the problem.

Regarding mesh, I can see that many wants to nest inside global models
(Hycom or NEMO) and do not want/need to download the whole grid domain just to get stripes along boundary. In other words, getting stripe along constant lon, and/or lat of global model saves disk/time using opendap subsettings if your ROMS grid is along that line.

I am curious about your opinion for orthogonality errors, how they plaque solution, what do you think is acceptable limit, is it only causing problem along boundary, how do they impact solution (bad conservation?) etc.

Thanks for reply in advance,

Best
Ivica

shchepet · #22 Unread post by **shchepet** » Wed Nov 28, 2018 9:22 am

What is the acceptable level of orthogonality errors? I would say none. Or, more
precisely, the question is irrelevant, because orthogonality errors can be kept at
the level roundoff errors of double precision -- it is technologically achievable,
and therefore should be done this way. Then there is no point to think about how
orthogonality errors corrupt the solution.

Regarding mesh: depending on the shape of the domain and coastline curvilinear grids
may be very useful, and are presently underutilized. The main reason, in my view, is
embarrassing lack of necessary software tools to generate such grids and to manipulate
data, pre- and post-processing.

If we ever to build a working sea-grid kind tool, it should be built around Laplace
equation solver with Mehrstellen discretization: it is a 9-point scheme, it is
compact fourth-order accurate. These is no point to have anything less accurate.
Laplace equations for both horizontal coordinates are coupled through boundary
conditions. Currently they are not.

Current/previous/existing approaches try to adapt ready-to-use algorithms, whatever
is available in Matlab or elsewhere. But if nothing fits precisely, then basically
settle for surrogates or whatever means to have it half-way, or simply give up.

An example? bilinear interpolation problem. Suppose point (xc,yc) is somewhere inside
the trapezoidal element defined by 4 points (i,j), (i+1,j), (i+1,j+1) (i,j+1) of the
"parent" (generally speaking orthogonal curvilinear) grid; then there must exist two
number p,q bounded by [0,1] each, such that

Code: Select all

xc=(1-p)*(1-q)*x(i,j) + p*(1-q)*x(i+1,j) + (1-p)*q*x(i,j+1) +p*q*x(i+1,j+1)
yc=(1-p)*(1-q)*y(i,j) + p*(1-q)*y(i+1,j) + (1-p)*q*y(i,j+1) +p*q*y(i+1,j+1)

assuming that xc,yc, and all x(i,j),y(i,j) are known, find p, and q.

How to solve it?

ROMS community answer: it cannot be done.

It is nonlinear problem.

Matlab does not provide it.

Therefore we use griddata - we have no other choice.

Griddata splits the trapezoidal element into two triangles. How does it split? By one
of the diagonals. Which one? Either one, 50% chance (b.t.w, the diagonals are exactly
equal to each other, if it is orthogonal curvilinear grid). Then in each triangle it can
uniquely define a linear function of the coordinates x,y and solve for the coefficients.

This is not exactly what we want, but this is what we settle for.

And before one gets to the above problem of finding p,q one should find indices (i,j)
of the parent grid cell which contains point (xc,yc). This requires search. A not so
trivial problem if the grid is curvilinear.

Again, Matlab does not provide it.

so, griddata.

Regarding mesh, I can see that many wants to nest inside global models (Hycom or
NEMO) and do not want/need to download the whole grid domain just to get stripes
along the boundary. In other words, getting stripe along constant lon, and/or lat of global
model saves disk/time using opendap subsettings if your ROMS grid is along that line

So desire to have straight lines along latitude or longitude comes from the fact that it
is easy to cut out data from a global dataset, which which is in lon-lat coordinates.

Getting the whole global data is out of question because of its size.

Frankly I would cut-out a logically-rectangular region which fully contains my ROMS model
domain and use 2D interpolation with rotation of velocity vectors as necessary.

In addition, in terms of read-write performance of netCDF files, getting 1D strip vs. getting
2D patch takes comparable time:

netCDF-3: assuming that longitude corresponds to 1st netCDF dimension, latitude to the 2nd,
reading 1D strip along longitude (1st dimension) takes almost no time: you read from the disk
only what you need; in contrast, getting 1D strip along 2nd dimension is comparable to reading
the entire 2D data bounded by the minimum and maximum latitudes of the desired strip.

netCDF4: assuming the same about 1st and 2nd dimensions: getting 1D strips in lon and lat
direction takes about the same time and requires reading the set of blocks which contain
all the data. Simply put, your reading becomes reading of 2D data regardless along which
direction you need your 1D strip. Depending on the netCDF/HDF5 block size vs. your ROMS
grid geographical extents, it may take exactly the same time as to cut out the entire 2D
region, or less.

CharlesJames · Mon Mar 25, 2019 12:41 am

shchepet wrote:The answer in "no", all three of them, A,B,C, cannot be made straight at the same time.

This is the nature of conformal mapping: once A and C are straight, there is no way to fit a
curve which would go along the coast, and be perpendicular to A and C at the junction points.

B can me made straight easily. but in this case A and C will be curved.

A can be made straight, but then B and C be curved, and the remaining line on the top would
have to enter deeply into the land and get a lot of masking.

Making C straight? I see no point as it is inside the mask any way.

Besides, there is no wisdom even trying to the make side boundaries straight? Even if some
of them are straight, once you depart from the boundary go into the domain, the coordinate
lines are no longer straight any way.

It has to be recognized that this analytical curvilinear grid is optimal by the design in sense
tat it is:

(1) exactly othrogonal;

(2) perfectly isotropic, pm=pn at every grid cell (same as local dx=dy at every cell, despite
that fact that both dx, dy change from cell to cell;

(3) with some diligence has minimal possible non-uniformity of the resolution: the ratio
min(dx)/max(dx) kept as small as possible for a given situation of the particular shape
of the domain and coastline;

The fundamental flaw of seagrid, as well as the more recent package GridBuilder which
replaces it, is that they provide no means to enforce perpenducular junction of the side
boundaries. Both packages allow you to define points at the side and move them in and
out to achieve the desired shape, but the junction points end up be not perpendicular,
and it is basically left up to the user to make them as close to perpendicular as possible --
it is possible in theory, but is very hard to achieve in practice. This makes it very hard to
construct a satisfactory grid: everybody complains about corner problems.

This also leads to the very long-standing confusion in ROMS community (and in POM
community before) that one can fit curvilinear grid to virtually any shape -- as for the
contour yes, you can fit it, but if the corner is not exactly 90 degrees, once you start
filling up the grid points by solving Dirichet elliptic equation, the grid points tend to
either concentrate toward the corner (if >90 degrees) or, conversely, avoid it (< 90).

I did add a new grid format to GridBuilder for purely orthogonal grids last year at the request of John Wilkin, it is similar to the rectangular format but it remains orthogonal (to working precision) under expansion, rotation and translation. In orthogonal mode you can't use control points to generate exotic curvilinear grids (for that you have to use the Free mode which replicates the orignal SeaGrid functionality) but it does still allow telescoping with sliders.

Ocean Modeling Discussion

Set a land mask on the sea in SeaGrid

Set a land mask on the sea in SeaGrid

Re: Set a land mask on the sea in SeaGrid

Re: Set a land mask on the sea in SeaGrid

Re: Set a land mask on the sea in SeaGrid

Re: Set a land mask on the sea in SeaGrid

Re: Set a land mask on the sea in SeaGrid

Re: Set a land mask on the sea in SeaGrid

Re: Set a land mask on the sea in SeaGrid

Re: Set a land mask on the sea in SeaGrid

Re: Set a land mask on the sea in SeaGrid

Re: Set a land mask on the sea in SeaGrid

Re: Set a land mask on the sea in SeaGrid

Re: Set a land mask on the sea in SeaGrid

Re: Set a land mask on the sea in SeaGrid

Re: Set a land mask on the sea in SeaGrid

Re: Set a land mask on the sea in SeaGrid

Re: Set a land mask on the sea in SeaGrid

Re: Set a land mask on the sea in SeaGrid

Re: Set a land mask on the sea in SeaGrid

Re: Set a land mask on the sea in SeaGrid

Re: Set a land mask on the sea in SeaGrid

Re: Set a land mask on the sea in SeaGrid

Re: Set a land mask on the sea in SeaGrid