% Main program to generate the grid file for ROMS. % Please mannually step through this script and edit when % it is requested. You will find a line marked EDIT. % E. Di Lorenzo - May 2003 disp('Please mannually step through this script'); return % This script will help you generate rectangular curvilinear grids only. %========================================================== % Preparing the grid param %========================================================== % EDIT: set a name to recognize the grid nameit='ias40'; % 8 km grid % EDIT: set path where you want to store the grid. outdir=['./']; if ~exist(outdir), unix(['mkdir ',outdir]);end grdfile=[outdir,nameit,'_grid.nc']; % EDIT: insert detailed coastline file if wanted. You can extract it % form the web at: rimmer.ngdc.noaa.gov/coast % After you donwload the file use ConvertCstdat2mat.m % otherwise use default. The format is lon,lat in the mat file. seagridCoastline=which('rgrd_CoastlineWorld.mat'); % default centered on Atlantic/Indian seagridCoastline=which('rgrd_CoastlineWorldPacific.mat'); % default centered on Pacific rgrd_DrawCstLine(seagridCoastline); % zoom in the picture to approximately select the location of % your final grid. [Lons]=get(gca,'xlim'); LonMin=Lons(1); LonMax=Lons(2); [Lats]=get(gca,'ylim'); LatMin=Lats(1); LatMax=Lats(2); close % EDIT: set ax= [LonMin, LonMax, LatMin, LatMax] % so that you include the location of the grid your % are building ax= [LonMin, LonMax, LatMin, LatMax] % EDIT: put the name of the file in which the scirpt will % save the 4 corners of your grid in LON,LAT . CornerFile= [outdir,'GridCorners-',nameit,'.dat']; % yu have two options to generate the corners of the grid: % OPT: 1 % If your corner are the one in ax then just write the % corner file. WriteCornerFileFromAX.m rgrd_WriteCornerFileFromAX(ax,CornerFile); % OPT: 2 % IF you already have a corner file than skip the FindGridCorners % use this routine to design grid if needed. FindGridCorners.m rgrd_FindGridCorners (ax,seagridCoastline,CornerFile); close all % Plot grid rgrd_PlotGridCorners(ax,CornerFile,seagridCoastline); % Find the distance and needed reolution for the grid. % TellMeCornerDist.m [x_spacing,y_spacing]=rgrd_TellMeCornerDist(CornerFile); % remember these spacing, round them to be even. disp('You need to set these in Seagrid'); Cells_Edge_1 = round(y_spacing/2)*2 Cells_Edge_2 = round(x_spacing/2)*2 %========================================================== % Build grid - Running Seagrid %========================================================== % execute seagrid and prepare grid % in seagrid just % 1) load the coastline file and the boundary file % 2) set the spacing in setup using Cells_Edge_1, Cells_Edge_2 % 3) On the menu compute select "Set all water" % 4) then save the seagrid.mat file and exit. seagrid(seagridCoastline) % if you like you can also load the topography and use the % mask computation options. I prefer taking care of this later % in this script. It is up to you. % convert seagrid file to ROMS grid seagrid2roms('seagrid40.mat', grdfile); % EDIT add new grid in the rnt_gridinfo.m configfile=which('rnt_gridinfo'); unix(['vi ',configfile]); grd=rnt_gridload(nameit); %========================================================== % topography %========================================================== load H_julios h grd1=rnt_gridload('ias20'); h=rnt_griddata(grd1.lonr,grd1.latr,h,grd.lonr,grd.latr,'cubic'); %EDIT: set minimum depth for topography. hmin = 20; hraw=-h;hraw(hraw< hmin) = hmin; hmax=5500; hraw(hraw> hmax) = hmax; grd.h=hraw; % now match topography with N. Atlantic grid grd=rnt_gridload(nameit); load hnatl.mat hn=rnt_griddata(grd1.lonr,grd1.latr,hn,grd.lonr,grd.latr,'cubic'); hn(isnan(hn)) = grd.h(isnan(hn)); nbuff = 10; alfa=zeros(size(grd.h)); for i=0:nbuff %north alfa( end-i,1:end-i) = cos( i*pi/nbuff/2 ); %east alfa( 1:end-i, end-i) = cos( i*pi/nbuff/2 ); end beta = 1 - alfa; hnew = alfa.*hn + beta.*grd.h; [I,J]=size(hnew); i=1:I; j=J-60:J; pcolor(grd.lonr(i,j), grd.latr(i,j), grd.maskr(i,j).*(hnew(i,j) - hn(i,j))); shading interp; colorbar % rfac=0.3; % maximum slope of topography allowed in smoothing % tmp = smoothgrid(hraw,hmin,rfac); % rnt_plcm(rvalue(tmp),grd); title 'rvalue h'; % Will try shapiro to preserve depth of straights tmp = shapiro2(hnew,2,2); % once filtered for i=1:9 tmp = shapiro2(tmp,2,2); % once filtered end rnt_plcm(rvalue(tmp),grd); title 'rvalue h'; %save topo nc=netcdf(grd.grdfile,'w'); nc{'h'}(:,:)=tmp'; nc{'hraw'}(1,:,:)=hraw'; close(nc) grd=rnt_gridload(nameit); nat=rnt_gridload('natl'); [I,J] = rgrd_FindInsideIJ(grd.lonr,grd.latr,nat.lonr,nat.latr); Irange = [I(1) I(end)] Jrange = [J(1) J(end)] % end smoothing ---------------------------------------- figure; grd=rnt_gridload(nameit); subplot(2,2,1); rnt_plc3(1./grd.pm/1000,grd,0,5,0,0); title 'DX km' subplot(2,2,2); rnt_plc3(grd.h,grd,0,5,0,0); title (['h (min depth ',num2str(grd.hc),' )']) subplot(2,2,3); rnt_plc3(rvalue(grd.h),grd,0,5,0,0); title 'rvalue h' subplot(2,2,4); rnt_plc3(grd.h-sq(grd.hraw(1,:,:)),grd,0,5,0,0); title 'h - hraw' %========================================================== % Masking %========================================================== grd=rnt_gridload(nameit); mask=grd.maskr; mask(grd.h <= hmin) =0; %/sdb/edl/ROMS-pak/matlib/rgrd/rgrd_GetMask.m % you can also interpolate the mask from an existing grid % before you do the final re-touches grd1=rnt_gridload('ias20'); mask = rgrd_GetMask(grd1,grd); mask(isnan(mask)) =0; %save raw mask nc=netcdf(grdfile,'w'); nc{'mask_rho'}(:)=mask'; close(nc); % utility from Rutgers modified to use rnt_griddata % will produce a file called ijcoast.mat seagridCoastline=grd.cstfile; [C]=rgrd_ijcoast(grdfile,seagridCoastline); editmask(grdfile,'ijcoast.mat'); %by Andrey Scherbina % if you do not want to change anything still % run editmask, press the SAVE button once and % then exit. Thanks. % Your grid is ready. grd=rnt_gridload(nameit);