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/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');
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');
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 ready-to-go version of the same mask post-processed by single_connect and commit_rmask.
- EastPac_curv_mask1.png (4.05 KiB) Viewed 53219 times
-
- This is "raw" mask generated by gshhs_to_roms_mask from USGS coastline data as seen in grid-index coordinates (basically by ncview)
- EastPac_curv_mask.png (5 KiB) Viewed 53219 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 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.
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.
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 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.
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.
- Attachments
-
- tools.tar
- source code
- (1.48 MiB) Downloaded 805 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 -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
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
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 copy-pasting and almost without
ever touching keyboard. See README.compiling_sequence inside the attached
netcdf4_help.tar.
- Attachments
-
- netcdf4_help.tar
- help for zlib-hdf5-netcdf compiling
- (20 KiB) Downloaded 761 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 697 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 08901-8521, USA. ph: 609-630-0559 jwilkin@rutgers.edu
71 Dudley Rd, New Brunswick, NJ 08901-8521, USA. ph: 609-630-0559 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 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);
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" ;
}
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.
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 639 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 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
...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), 1-mss(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), ione-mss(i,j+1))
enddo
enddo
- Attachments
-
- tools.tar
- updated November2018
- (1.54 MiB) Downloaded 609 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 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).
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).
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: 172
- Joined: Mon May 05, 2003 2:41 pm
- Location: The University of Western Australia, Perth, Australia
- Contact:
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 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
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 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.
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)
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 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
- Posts: 43
- 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 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).