Set a land mask on the sea in SeaGrid

 Posts: 48
 Joined: Tue Aug 04, 2015 4:42 pm
 Location: Universidad del Mar (UMAR), Mexico
 Contact:
Set a land mask on the sea in SeaGrid
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.
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.
Re: Set a land mask on the sea in SeaGrid
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???):
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()

 Posts: 48
 Joined: Tue Aug 04, 2015 4:42 pm
 Location: Universidad del Mar (UMAR), Mexico
 Contact:
Re: Set a land mask on the sea in SeaGrid
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/roms3.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');
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/roms3.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');
Re: Set a land mask on the sea in SeaGrid
Is this the kind of grid you want to generate?
Re: Set a land mask on the sea in SeaGrid
The procedure for generating ROMS land mask from USGS coastline data described here:
viewtopic.php?f=23&t=3878&p=14908
viewtopic.php?f=23&t=3878&p=14908
 Attachments

 This is readytogo version of the same mask postprocessed by single_connect and commit_rmask.
 EastPac_curv_mask1.png (4.05 KiB) Viewed 4731 times

 This is "raw" mask generated by gshhs_to_roms_mask from USGS coastline data as seen in gridindex coordinates (basically by ncview)
 EastPac_curv_mask.png (5 KiB) Viewed 4731 times

 Posts: 48
 Joined: Tue Aug 04, 2015 4:42 pm
 Location: Universidad del Mar (UMAR), Mexico
 Contact:
Re: Set a land mask on the sea in SeaGrid
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 IFexpression at (1)
mpc.f:336:19: Error: Syntax error in IFexpression at (1)
mpc.f:337:20: Error: Syntax error in IFexpression at (1)
mpc.f:338:21: Error: Syntax error in IFexpression at (1)
mpc.f:339:22: Error: Syntax error in IFexpression at (1)
mpc.f:340:23: Error: Syntax error in IFexpression 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 nonblanc str
1
Error: Expecting END IF statement at (1)
mpc.f:374:17: Error: Syntax error in IFexpression 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 IFexpression at (1)
mpc.f:391:14:
bffr(k)=str(k+is1)
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 newstyle
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+ieis+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+ieis+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 IFexpression 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 IFexpression 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 IFexpression at (1)
mpc.f:641:21: Error: Syntax error in IFexpression 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 IFexpression 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 IFexpression 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 IFexpression at (1)
mpc.f:671:25: Error: Syntax error in IFexpression 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(m1) == ' ')
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 IFexpression 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 IFexpression 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 IFexpression 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 IFexpression 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 IFexpression at (1)
mpc.f:820:18:
str(i+j1)=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 IFexpression 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(im)=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.
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 IFexpression at (1)
mpc.f:336:19: Error: Syntax error in IFexpression at (1)
mpc.f:337:20: Error: Syntax error in IFexpression at (1)
mpc.f:338:21: Error: Syntax error in IFexpression at (1)
mpc.f:339:22: Error: Syntax error in IFexpression at (1)
mpc.f:340:23: Error: Syntax error in IFexpression 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 nonblanc str
1
Error: Expecting END IF statement at (1)
mpc.f:374:17: Error: Syntax error in IFexpression 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 IFexpression at (1)
mpc.f:391:14:
bffr(k)=str(k+is1)
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 newstyle
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+ieis+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+ieis+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 IFexpression 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 IFexpression 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 IFexpression at (1)
mpc.f:641:21: Error: Syntax error in IFexpression 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 IFexpression 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 IFexpression 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 IFexpression at (1)
mpc.f:671:25: Error: Syntax error in IFexpression 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(m1) == ' ')
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 IFexpression 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 IFexpression 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 IFexpression 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 IFexpression 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 IFexpression at (1)
mpc.f:820:18:
str(i+j1)=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 IFexpression 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(im)=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.
Re: Set a land mask on the sea in SeaGrid
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 userspecified location and rotated by userspecified 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 xaxis  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 yaxis  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.
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 userspecified location and rotated by userspecified 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 xaxis  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 yaxis  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.
 Attachments

 tools.tar
 source code
 (1.48 MiB) Downloaded 157 times
Re: Set a land mask on the sea in SeaGrid
you can fit weird shapes after some practice...
 Attachments
Last edited by shchepet on Thu Nov 22, 2018 10:27 am, edited 1 time in total.

 Posts: 48
 Joined: Tue Aug 04, 2015 4:42 pm
 Location: Universidad del Mar (UMAR), Mexico
 Contact:
Re: Set a land mask on the sea in SeaGrid
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 fnosecondunderscore 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
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 fnosecondunderscore 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
Re: Set a land mask on the sea in SeaGrid
This one you should be able to overcome yourself. At first, my package comes with severalError: Symbol ‘nf_netcdf4’ at (1) has no IMPLICIT type
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 copypasting and almost without
ever touching keyboard. See README.compiling_sequence inside the attached
netcdf4_help.tar.
 Attachments

 netcdf4_help.tar
 help for zlibhdf5netcdf compiling
 (20 KiB) Downloaded 96 times

 Posts: 48
 Joined: Tue Aug 04, 2015 4:42 pm
 Location: Universidad del Mar (UMAR), Mexico
 Contact:
Re: Set a land mask on the sea in SeaGrid
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.
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.
 Attachments

 GT5s_1_12_roms_grd_sponge2_mask.nc
 (17.76 MiB) Downloaded 69 times
Re: Set a land mask on the sea in SeaGrid
Try the tool editmask which is in the myroms/landmask folder of Matlab tools
John Wilkin: DMCS Rutgers University
71 Dudley Rd, New Brunswick, NJ 089018521, USA. ph: 6096300559 jwilkin@rutgers.edu
71 Dudley Rd, New Brunswick, NJ 089018521, USA. ph: 6096300559 jwilkin@rutgers.edu
Re: Set a land mask on the sea in SeaGrid
John,Try the tool editmask which is in the myroms/landmask folder of Matlab tools
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 20180211 23:15:37Z arango $
and $ Copyright (c) 20022018 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:Jm1,1:Im1)=mask_rho(2:Jm ,2:Im ).* ...
mask_rho(2:Jm ,1:Im1).* ...
mask_rho(1:Jm1,2:Im ).* ...
mask_rho(1:Jm1,1:Im1);
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 RHOpoints" ;
lon_rho:units = "degree East" ;
double lat_rho(eta_rho, xi_rho) ;
lat_rho:long_name = "latitude of RHOpoints" ;
lat_rho:units = "degree North" ;
double lon_psi(eta_v, xi_u) ;
lon_psi:long_name = "longitude of PSIpoints" ;
lon_psi:units = "degree East" ;
double lat_psi(eta_v, xi_u) ;
lat_psi:long_name = "latitude of PSIpoints" ;
lat_psi:units = "degree North" ;
double pm(eta_rho, xi_rho) ;
pm:long_name = "curvilinear coordinate metric in XIdirection" ;
pm:units = "meter1" ;
double pn(eta_rho, xi_rho) ;
pn:long_name = "curvilinear coordinate metric in ETAdirection" ;
pn:units = "meter1" ;
double angle(eta_rho, xi_rho) ;
angle:long_name = "angle between east and XIdirections" ;
angle:units = "degrees" ;
double f(eta_rho, xi_rho) ;
f:long_name = "Coriolis parameter at RHOpoints" ;
f:units = "second1" ;
double hraw(eta_rho, xi_rho) ;
hraw:long_name = "raw bathymetry at RHOpoints" ;
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 RHOpoints" ;
hsmth:units = "meter" ;
hsmth:hmin = 5. ;
hsmth:hmax = 4950. ;
hsmth:r_max = 0.15 ;
hsmth:method = "maskdependent rxconditioned logsmoothing" ;
hsmth:iters_cond = 1000 ;
double mask_rho(eta_rho, xi_rho) ;
mask_rho:long_name = "land/sea mask at RHOpoints" ;
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 = "parenttochild 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 = "06May2017" ;
}
is because if topography needs to be "resmoothed" or "remerged" with parent grid in can be done stagebystage rather than starting
from the very beginning.
Variable "wgt" is the merging function, h = (1wgt)*hsmth + wgt*h_parent, where wgt=0 everywhere inside the domain, except some areas
(bands) near the open boundaries; h_parent is parentgrid 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.
Re: Set a land mask on the sea in SeaGrid
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.f70:
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.f6:
& iA_q)
1
Error: Bad continuation line at (1)
vort.f6:
& iA_q)
1
Error: Unclassifiable statement at (1)
vort.f72: Error: Unexpected !$OMP END PARALLEL statement at (1)
&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
Annex more information in error.zip
I hope you can help me, thank you.
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.f70:
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.f6:
& iA_q)
1
Error: Bad continuation line at (1)
vort.f6:
& iA_q)
1
Error: Unclassifiable statement at (1)
vort.f72: Error: Unexpected !$OMP END PARALLEL statement at (1)
&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
Annex more information in error.zip
I hope you can help me, thank you.
 Attachments

 Error.zip
 More images about error and a file.txt whit all output of terminal
 (325.17 KiB) Downloaded 54 times
Re: Set a land mask on the sea in SeaGrid
Obviously this should be
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.
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
in the original vort.F file.
...interestingly Intel Ifort compiler does not care.
Re: Set a land mask on the sea in SeaGrid
Thank you so much, I have solved the problem. I think it's because I use Gfortran.
Re: Set a land mask on the sea in SeaGrid
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:
where,strictly speaking, in the above expression 1mss(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
...it is as silly as this.
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), 1mss(i,j+1))
enddo
enddo
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), ionemss(i,j+1))
enddo
enddo
 Attachments

 tools.tar
 updated November2018
 (1.54 MiB) Downloaded 54 times
Re: Set a land mask on the sea in SeaGrid
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
I want to make a grid with sides A, B, C completely right, similar to the attached image
Re: Set a land mask on the sea in SeaGrid
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 nonuniformity 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 longstanding 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).
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 nonuniformity 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 longstanding 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).
Re: Set a land mask on the sea in SeaGrid
Thank you
I understand.
I made a useful grid with this tool.
I understand.
I made a useful grid with this tool.
 jivica
 Posts: 137
 Joined: Mon May 05, 2003 2:41 pm
 Location: The University of Western Australia, Perth, Australia
Re: Set a land mask on the sea in SeaGrid
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
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
Re: Set a land mask on the sea in SeaGrid
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 postprocessing.
If we ever to build a working seagrid kind tool, it should be built around Laplace
equation solver with Mehrstellen discretization: it is a 9point scheme, it is
compact fourthorder 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 readytouse algorithms, whatever
is available in Matlab or elsewhere. But if nothing fits precisely, then basically
settle for surrogates or whatever means to have it halfway, 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
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.
is easy to cut out data from a global dataset, which which is in lonlat coordinates.
Getting the whole global data is out of question because of its size.
Frankly I would cutout a logicallyrectangular region which fully contains my ROMS model
domain and use 2D interpolation with rotation of velocity vectors as necessary.
In addition, in terms of readwrite performance of netCDF files, getting 1D strip vs. getting
2D patch takes comparable time:
netCDF3: 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.
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 postprocessing.
If we ever to build a working seagrid kind tool, it should be built around Laplace
equation solver with Mehrstellen discretization: it is a 9point scheme, it is
compact fourthorder 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 readytouse algorithms, whatever
is available in Matlab or elsewhere. But if nothing fits precisely, then basically
settle for surrogates or whatever means to have it halfway, 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=(1p)*(1q)*x(i,j) + p*(1q)*x(i+1,j) + (1p)*q*x(i,j+1) +p*q*x(i+1,j+1)
yc=(1p)*(1q)*y(i,j) + p*(1q)*y(i+1,j) + (1p)*q*y(i,j+1) +p*q*y(i+1,j+1)
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.
So desire to have straight lines along latitude or longitude comes from the fact that itRegarding 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
is easy to cut out data from a global dataset, which which is in lonlat coordinates.
Getting the whole global data is out of question because of its size.
Frankly I would cutout a logicallyrectangular region which fully contains my ROMS model
domain and use 2D interpolation with rotation of velocity vectors as necessary.
In addition, in terms of readwrite performance of netCDF files, getting 1D strip vs. getting
2D patch takes comparable time:
netCDF3: 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
 Posts: 33
 Joined: Thu May 24, 2007 12:12 pm
 Location: South Australian Research and Development Institute
Re: Set a land mask on the sea in SeaGrid
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.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 nonuniformity 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 longstanding 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).