[1] | 1 | %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
| 2 | % |
---|
| 3 | % Build a ROMS grid file |
---|
| 4 | % |
---|
| 5 | % |
---|
| 6 | % Pierrick Penven, IRD, 2002. |
---|
| 7 | % |
---|
| 8 | % Version of 10-Oct-2002 |
---|
| 9 | %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
| 10 | clear all |
---|
| 11 | close all |
---|
| 12 | %%%%%%%%%%%%%%%%%%%%% USERS DEFINED VARIABLES %%%%%%%%%%%%%%%%%%%%%%%% |
---|
| 13 | % |
---|
| 14 | % Title |
---|
| 15 | % |
---|
| 16 | title='TEST'; |
---|
| 17 | config='test'; |
---|
| 18 | % |
---|
| 19 | % Grid file name |
---|
| 20 | % |
---|
| 21 | grdname='test_grd.nc'; |
---|
| 22 | % |
---|
| 23 | % Slope parameter (r=grad(h)/h) maximum value for topography smoothing |
---|
| 24 | % |
---|
| 25 | rtarget=0.15; |
---|
| 26 | % |
---|
| 27 | % Grid dimensions: |
---|
| 28 | % lonmin : Minimum longitude [degree east] |
---|
| 29 | % lonmax : Maximum longitude [degree east] |
---|
| 30 | % latmin : Minimum latitude [degree north] |
---|
| 31 | % latmax : Maximum latitude [degree north] |
---|
| 32 | if config=='test' |
---|
| 33 | % |
---|
| 34 | % Test |
---|
| 35 | lon1=43.9; |
---|
| 36 | lat1=8.7; |
---|
| 37 | lon2=34.1; |
---|
| 38 | lat2=31.2; |
---|
| 39 | lon3=46.9; |
---|
| 40 | lat3=10.7; |
---|
| 41 | |
---|
| 42 | end |
---|
| 43 | % |
---|
| 44 | % Grid resolution [degree] |
---|
| 45 | % |
---|
| 46 | dl=1/4; |
---|
| 47 | % |
---|
| 48 | % Minimum depth [m] |
---|
| 49 | % |
---|
| 50 | hmin=10; |
---|
| 51 | % |
---|
| 52 | % Topography netcdf file name (ETOPO 2) |
---|
| 53 | % |
---|
| 54 | topofile='../Topo/etopo2.nc'; |
---|
| 55 | % |
---|
| 56 | % GSHSS user defined coastline (see m_map) |
---|
| 57 | % |
---|
| 58 | coastfileplot=''; |
---|
| 59 | % |
---|
| 60 | % |
---|
| 61 | %%%%%%%%%%%%%%%%%%% END USERS DEFINED VARIABLES %%%%%%%%%%%%%%%%%%%%%%% |
---|
| 62 | % |
---|
| 63 | % Title |
---|
| 64 | % |
---|
| 65 | disp(' ') |
---|
| 66 | disp([' Making the grid: ',grdname]) |
---|
| 67 | disp(' ') |
---|
| 68 | disp([' Title: ',title]) |
---|
| 69 | disp(' ') |
---|
| 70 | disp([' Resolution: 1/',num2str(1/dl),' deg']) |
---|
| 71 | % |
---|
| 72 | % Get the Longitude and latitude |
---|
| 73 | % |
---|
| 74 | theta=atan((lat3-lat1)/(lon3-lon1)) |
---|
| 75 | l=(lon2-lon1)*cos(theta)+(lat2-lat1)*sin(theta) |
---|
| 76 | L=(lat2-lat1)*cos(theta)-(lon2-lon1)*sin(theta) |
---|
| 77 | x=(0:dl:l); |
---|
| 78 | y=(0:dl:L); |
---|
| 79 | [xr,yr]=meshgrid(x,y); |
---|
| 80 | Lonr=lon1+xr*cos(theta)-yr*sin(theta); |
---|
| 81 | Latr=lat1+xr*sin(theta)+yr*cos(theta); |
---|
| 82 | [Lonu,Lonv,Lonp]=rho2uvp(Lonr); |
---|
| 83 | [Latu,Latv,Latp]=rho2uvp(Latr); |
---|
| 84 | % |
---|
| 85 | % Create the grid file |
---|
| 86 | % |
---|
| 87 | disp(' ') |
---|
| 88 | disp(' Create the grid file...') |
---|
| 89 | [M,L]=size(Latp); |
---|
| 90 | disp([' LLm = ',num2str(L-1)]) |
---|
| 91 | disp([' MMm = ',num2str(M-1)]) |
---|
| 92 | create_grid(L,M,grdname,title) |
---|
| 93 | % |
---|
| 94 | % Fill the grid file |
---|
| 95 | % |
---|
| 96 | disp(' ') |
---|
| 97 | disp(' Fill the grid file...') |
---|
| 98 | nc=netcdf(grdname,'write'); |
---|
| 99 | nc{'lat_u'}(:)=Latu; |
---|
| 100 | nc{'lon_u'}(:)=Lonu; |
---|
| 101 | nc{'lat_v'}(:)=Latv; |
---|
| 102 | nc{'lon_v'}(:)=Lonv; |
---|
| 103 | nc{'lat_rho'}(:)=Latr; |
---|
| 104 | nc{'lon_rho'}(:)=Lonr; |
---|
| 105 | nc{'lat_psi'}(:)=Latp; |
---|
| 106 | nc{'lon_psi'}(:)=Lonp; |
---|
| 107 | result=close(nc); |
---|
| 108 | % |
---|
| 109 | % Compute the metrics |
---|
| 110 | % |
---|
| 111 | disp(' ') |
---|
| 112 | disp(' Compute the metrics...') |
---|
| 113 | [pm,pn,dndx,dmde]=get_metrics(grdname); |
---|
| 114 | xr=0.*pm; |
---|
| 115 | yr=xr; |
---|
| 116 | for i=1:L |
---|
| 117 | xr(:,i+1)=xr(:,i)+2./(pm(:,i+1)+pm(:,i)); |
---|
| 118 | end |
---|
| 119 | for j=1:M |
---|
| 120 | yr(j+1,:)=yr(j,:)+2./(pn(j+1,:)+pn(j,:)); |
---|
| 121 | end |
---|
| 122 | [xu,xv,xp]=rho2uvp(xr); |
---|
| 123 | [yu,yv,yp]=rho2uvp(yr); |
---|
| 124 | dx=1./pm; |
---|
| 125 | dy=1./pn; |
---|
| 126 | dxmax=max(max(dx/1000)); |
---|
| 127 | dxmin=min(min(dx/1000)); |
---|
| 128 | dymax=max(max(dy/1000)); |
---|
| 129 | dymin=min(min(dy/1000)); |
---|
| 130 | disp(' ') |
---|
| 131 | disp([' Min dx=',num2str(dxmin),' km - Max dx=',num2str(dxmax),' km']) |
---|
| 132 | disp([' Min dy=',num2str(dymin),' km - Max dy=',num2str(dymax),' km']) |
---|
| 133 | % |
---|
| 134 | % Add topography from topofile |
---|
| 135 | % |
---|
| 136 | disp(' ') |
---|
| 137 | disp(' Add topography...') |
---|
| 138 | h=add_topo(grdname,topofile); |
---|
| 139 | % |
---|
| 140 | % Compute the mask |
---|
| 141 | % |
---|
| 142 | maskr=h>0; |
---|
| 143 | maskr=process_mask(maskr); |
---|
| 144 | if config=='kago' |
---|
| 145 | maskr(Latr>45.5)=0; |
---|
| 146 | end |
---|
| 147 | if config=='test' |
---|
| 148 | maskr(Lonr<41.5 & Latr<14.4)=0.; |
---|
| 149 | maskr(Lonr<40.5 & Latr<14.6)=0.; |
---|
| 150 | end |
---|
| 151 | [masku,maskv,maskp]=uvp_mask(maskr); |
---|
| 152 | % |
---|
| 153 | % Smooth the topography |
---|
| 154 | % |
---|
| 155 | h = smoothgrid(h,hmin,rtarget); |
---|
| 156 | % |
---|
| 157 | % Angle between XI-axis and the direction |
---|
| 158 | % to the EAST at RHO-points [radians]. |
---|
| 159 | % |
---|
| 160 | angle=get_angle(Latu,Lonu); |
---|
| 161 | % |
---|
| 162 | % Coriolis parameter |
---|
| 163 | % |
---|
| 164 | f=4*pi*sin(pi*Latr/180)/(24*3600); |
---|
| 165 | % |
---|
| 166 | % Write it down |
---|
| 167 | % |
---|
| 168 | disp(' ') |
---|
| 169 | disp(' Write it down...') |
---|
| 170 | nc=netcdf(grdname,'write'); |
---|
| 171 | nc{'h'}(:)=h; |
---|
| 172 | nc{'pm'}(:)=pm; |
---|
| 173 | nc{'pn'}(:)=pn; |
---|
| 174 | nc{'dndx'}(:)=dndx; |
---|
| 175 | nc{'dmde'}(:)=dmde; |
---|
| 176 | nc{'mask_u'}(:)=masku; |
---|
| 177 | nc{'mask_v'}(:)=maskv; |
---|
| 178 | nc{'mask_psi'}(:)=maskp; |
---|
| 179 | nc{'mask_rho'}(:)=maskr; |
---|
| 180 | nc{'x_u'}(:)=xu; |
---|
| 181 | nc{'y_u'}(:)=yu; |
---|
| 182 | nc{'x_v'}(:)=xv; |
---|
| 183 | nc{'y_v'}(:)=yv; |
---|
| 184 | nc{'x_rho'}(:)=xr; |
---|
| 185 | nc{'y_rho'}(:)=yr; |
---|
| 186 | nc{'x_psi'}(:)=xp; |
---|
| 187 | nc{'y_psi'}(:)=yp; |
---|
| 188 | nc{'angle'}(:)=angle; |
---|
| 189 | nc{'f'}(:)=f; |
---|
| 190 | nc{'spherical'}(:)='T'; |
---|
| 191 | result=close(nc); |
---|
| 192 | disp(' ') |
---|
| 193 | disp([' Size of the grid: L = ',... |
---|
| 194 | num2str(L),' - M = ',num2str(M)]) |
---|
| 195 | % |
---|
| 196 | % make a plot |
---|
| 197 | % |
---|
| 198 | disp(' ') |
---|
| 199 | disp(' Do a plot...') |
---|
| 200 | themask=ones(size(maskr)); |
---|
| 201 | themask(maskr==0)=NaN; |
---|
| 202 | domaxis=[min(min(Lonr)) max(max(Lonr)) min(min(Latr)) max(max(Latr))]; |
---|
| 203 | colaxis=[min(min(h)) max(max(h))]; |
---|
| 204 | fixcolorbar([0.25 0.05 0.5 0.03],colaxis,... |
---|
| 205 | 'Topography',10) |
---|
| 206 | width=1; |
---|
| 207 | height=0.8; |
---|
| 208 | subplot('position',[0. 0.14 width height]) |
---|
| 209 | m_proj('mercator',... |
---|
| 210 | 'lon',[domaxis(1) domaxis(2)],... |
---|
| 211 | 'lat',[domaxis(3) domaxis(4)]); |
---|
| 212 | m_pcolor(Lonr,Latr,h.*themask); |
---|
| 213 | %shading interp |
---|
| 214 | caxis(colaxis) |
---|
| 215 | hold on |
---|
| 216 | m_contour(Lonr,Latr,h,[hmin 100 200 500 1000 2000 4000],'k--'); |
---|
| 217 | if ~isempty(coastfileplot) |
---|
| 218 | m_usercoast(coastfileplot,'patch',[.9 .9 .9]); |
---|
| 219 | else |
---|
| 220 | % m_gshhs_l('patch',[.9 .9 .9]) |
---|
| 221 | end |
---|
| 222 | m_grid('box','fancy',... |
---|
| 223 | 'xtick',5,'ytick',5,'tickdir','out',... |
---|
| 224 | 'fontsize',7); |
---|
| 225 | hold off |
---|
| 226 | % |
---|
| 227 | % End |
---|
| 228 | % |
---|
| 229 | %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
---|