![(please configure the [header_logo] section in trac.ini)](site/roms_agrif_sst.png)
Changeset 7
- Timestamp:
- 05/04/11 17:34:07 (13 years ago)
- Location:
- Roms_tools/Aforc_NCEP
- Files:
-
- 2 edited
-
get_NCEP_grid.m (modified) (4 diffs)
-
make_NCEP.m (modified) (25 diffs)
Legend:
- Unmodified
- Added
- Removed
-
Roms_tools/Aforc_NCEP/get_NCEP_grid.m
r1 r7 79 79 lat=lat(j); 80 80 % 81 % If we are in OpenDap (Get_My_Data=0) we need 82 % a shift of decal=1, because the indexes begin at 0 in OpenDap. 83 % If we use local data, Get_My_Data=1, there is no shift needed, 84 % then decal=0; 85 86 if Get_My_Data==1 87 decal=0; 88 else 89 decal=1; 90 end 91 % 81 92 if ~isempty(i1) 82 i1min=min(i1)- 1;83 i1max=max(i1)- 1;93 i1min=min(i1)-decal; 94 i1max=max(i1)-decal; 84 95 else 85 96 i1min=[]; … … 87 98 end 88 99 if ~isempty(i2) 89 i2min=min(i2)- 1;90 i2max=max(i2)- 1;100 i2min=min(i2)-decal; 101 i2max=max(i2)-decal; 91 102 else 92 103 i2min=[]; … … 94 105 end 95 106 if ~isempty(i3) 96 i3min=min(i3)- 1;97 i3max=max(i3)- 1;107 i3min=min(i3)-decal; 108 i3max=max(i3)-decal; 98 109 else 99 110 i3min=[]; … … 101 112 end 102 113 % 103 jmin=min(j)- 1;104 jmax=max(j)- 1;114 jmin=min(j)-decal; 115 jmax=max(j)-decal; 105 116 jrange=['[',num2str(jmin),':',num2str(jmax),']']; 106 117 -
Roms_tools/Aforc_NCEP/make_NCEP.m
r1 r7 2 2 % 3 3 % make_NCEP.m 4 % 4 % 5 5 % Create and fill frc and bulk files with NCEP data. 6 6 % (NCEP Reanalysis) … … 8 8 % The on-line reference to NCEP is at 9 9 % http://www.cdc.noaa.gov/cdc/reanalysis/ 10 % 11 % Further Information: 10 % 11 % Further Information: 12 12 % http://www.brest.ird.fr/Roms_tools/ 13 % 13 % 14 14 % This file is part of ROMSTOOLS 15 15 % … … 29 29 % MA 02111-1307 USA 30 30 % 31 % Copyright (c) 2005-2006 by Pierrick Penven 32 % e-mail:Pierrick.Penven@ird.fr 33 % 31 % Copyright (c) 2005-2006 by Pierrick Penven 32 % e-mail:Pierrick.Penven@ird.fr 33 % 34 34 % Adapted from a previous verions from 35 35 % Alvaro Peliz (U. Aveiro) & Patrick Marchesiello (IRD) - 2005 … … 48 48 % 49 49 if NCEP_version==1 50 frc_prefix=[frc_prefix,'_NCEP1_']; 50 frc_prefix=[frc_prefix,'_NCEP1_']; 51 51 blk_prefix=[blk_prefix,'_NCEP1_']; 52 52 elseif NCEP_version==2 53 frc_prefix=[frc_prefix,'_NCEP2_']; 53 frc_prefix=[frc_prefix,'_NCEP2_']; 54 54 blk_prefix=[blk_prefix,'_NCEP2_']; 55 55 end … … 103 103 disp(['====================']) 104 104 download_NCEP(Ymin,Ymax,Mmin,Mmax,lonmin,lonmax,latmin,latmax,... 105 105 NCEP_dir,NCEP_version,Yorig,Get_My_Data,My_NCEP_dir) 106 106 disp(['=====================']) 107 107 disp(['DOWNLOAD STEP FINISH']) … … 125 125 lat1=nc{'lat'}(:); 126 126 [lon1,lat1]=meshgrid(lon1,lat1); 127 127 128 128 if Get_My_Data~=1 129 129 mask=1-squeeze(nc{'landsfc'}(:)); … … 131 131 mask=1-squeeze(nc{'land'}(:)); 132 132 end 133 mask(mask==0)=NaN; 133 mask(mask==0)=NaN; 134 134 close(nc); 135 135 136 136 % 137 137 %Loop on the years and the months … … 144 144 % 145 145 for Y=Ymin:Ymax 146 if Y==Ymin 146 if Y==Ymin 147 147 mo_min=Mmin; 148 148 else … … 157 157 disp(' ') 158 158 disp(['Processing year ',num2str(Y),... 159 159 ' - month ',num2str(M)]) 160 160 disp(' ') 161 161 % … … 164 164 if Get_My_Data~=1 165 165 nc=netcdf([NCEP_dir,'tmp2m_Y',num2str(Y),'M',num2str(M),'.nc']); 166 elseif Get_My_Data==1 166 elseif Get_My_Data==1 167 167 nc=netcdf([NCEP_dir,'prate_Y',num2str(Y),'M',num2str(M),'.nc']); 168 168 end … … 180 180 disp(['tlen=',num2str(tlen)]) 181 181 disp(['Overlap is ',num2str(itolap_ncep),' it of 6 hours']) 182 disp(['Overlap is ',num2str(itolap_ncep/4),' days before and after']) 182 disp(['Overlap is ',num2str(itolap_ncep/4),' days before and after']) 183 183 time=0*(1:tlen); 184 time(itolap_ncep+1:tlen0+itolap_ncep)=NCEP_time; 184 time(itolap_ncep+1:tlen0+itolap_ncep)=NCEP_time; 185 185 disp(['=====================']) 186 186 disp('Compute time for roms file') 187 187 disp(['=====================']) 188 188 for aa=1:itolap_ncep 189 190 end 191 189 time(aa)=time(itolap_ncep+1)-(itolap_ncep+1-aa)*dt; 190 end 191 192 192 for aa=1:itolap_ncep 193 193 time(tlen0+itolap_ncep+aa)=time(tlen0+itolap_ncep)+aa*dt; … … 199 199 % Create the ROMS forcing files 200 200 blkname=[blk_prefix,'Y',num2str(Y),... 201 201 'M',num2str(M),nc_suffix]; 202 202 frcname=[frc_prefix,'Y',num2str(Y),... 203 203 'M',num2str(M),nc_suffix]; 204 204 if makeblk==1 205 206 205 disp(['Create a new bulk file: ' blkname]) 206 create_bulk(blkname,grdname,ROMS_title,time,0); 207 207 disp([' ']) 208 208 end 209 209 if makefrc==1 210 210 disp(['Create a new forcing file: ' frcname]) 211 211 disp([' ']) 212 213 214 215 212 create_forcing(frcname,grdname,ROMS_title,... 213 time,0,0,... 214 0,0,0,... 215 0,0,0,0,0,0) 216 216 end 217 217 % … … 219 219 % 220 220 if add_tides==1 221 221 add_tidal_data(tidename,grdname,frcname,Y,M) 222 222 end 223 223 % 224 224 % Open the ROMS forcing files 225 225 if makefrc==1 226 226 nc_frc=netcdf(frcname,'write'); 227 227 else 228 228 nc_frc=[]; 229 229 end 230 230 if makeblk==1 231 231 nc_blk=netcdf(blkname,'write'); 232 232 else 233 233 nc_blk=[]; 234 234 end 235 235 % … … 238 238 Ym=Y; 239 239 if Mm==0 240 241 242 end 243 240 Mm=12; 241 Ym=Y-1; 242 end 243 244 244 if Get_My_Data~=1 245 245 fname = [NCEP_dir,'tmp2m_Y',num2str(Ym),'M',num2str(Mm),'.nc']; 246 246 nc=netcdf([NCEP_dir,'tmp2m_Y',num2str(Ym),'M',num2str(Mm),'.nc']); 247 247 248 248 elseif Get_My_Data==1 249 249 fname = [NCEP_dir,'prate_Y',num2str(Ym),'M',num2str(Mm),'.nc']; … … 253 253 disp(' ') 254 254 disp('======================================================') 255 disp('Perform interpolations for the previous month') 255 disp('Perform interpolations for the previous month') 256 256 disp('======================================================') 257 257 disp(' ') 258 258 if exist(fname)==0 259 260 261 262 259 disp(['No data for the previous month: using current month']) 260 tndx=1; 261 Mm=M; 262 Ym=Y; 263 263 else 264 265 264 nc=netcdf(fname); 265 tndx=length(nc('time')); 266 266 if makefrc==1 267 267 for aa=1:itolap_ncep 268 268 nc_frc{'sms_time'}(aa)=nc{'time'}(tndx-(itolap_ncep-aa)); 269 269 end 270 270 end 271 271 % 272 273 272 if makeblk==1 273 for aa=1:itolap_ncep 274 274 nc_blk{'bulk_time'}(aa)=nc{'time'}(tndx-(itolap_ncep-aa)); 275 275 end 276 276 end 277 277 close(nc) 278 278 end … … 281 281 % 282 282 for aa=1:itolap_ncep 283 aa0=itolap_ncep-aa; 283 aa0=itolap_ncep-aa; 284 284 interp_NCEP(NCEP_dir,Ym,Mm,Roa,interp_method,lon1,lat1,... 285 286 end 287 %###################################################################### 288 % 285 mask,tndx-aa0,nc_frc,nc_blk,lon,lat,angle,aa,Get_My_Data) 286 end 287 %###################################################################### 288 % 289 289 disp(' ') 290 290 disp('======================================================') … … 297 297 298 298 for tndx=1:tlen0 299 300 301 299 if mod(tndx,20)==0 300 disp(['Step: ',num2str(tndx),' of ',num2str(tlen0)]) 301 end 302 302 interp_NCEP(NCEP_dir,Y,M,Roa,interp_method,lon1,lat1,... 303 mask,tndx,nc_frc,nc_blk,lon,lat,angle,tndx+itolap_ncep,Get_My_Data) 304 end 305 306 disp(' ') 307 disp('======================================================') 308 disp('Perform interpolations for next month') 303 mask,tndx,nc_frc,nc_blk,lon,lat,angle,tndx+itolap_ncep,Get_My_Data) 304 end 305 306 disp(' ') 307 disp('======================================================') 308 disp('Perform interpolations for next month') 309 309 disp('======================================================') 310 310 disp(' ') … … 314 314 Mp=M+1; 315 315 Yp=Y; 316 % 317 % Perform the interpolations for the next month 318 % 319 disp('Last steps') 316 320 if Mp==13 317 318 319 end 320 321 Mp=1; 322 Yp=Y+1; 323 end 324 321 325 if Get_My_Data~=1 322 326 fname=[NCEP_dir,'tmp2m_Y',num2str(Yp),'M',num2str(Mp),'.nc']; … … 324 328 fname=[NCEP_dir,'prate_Y',num2str(Yp),'M',num2str(Mp),'.nc']; 325 329 end 326 330 327 331 if exist(fname)==0 328 329 330 331 332 disp(['No data for the next month: using current month']) 333 tndx=tlen0; 334 Mp=M; 335 Yp=Y; 332 336 else 333 337 nc=netcdf(fname); 334 338 if makefrc==1 335 339 disp('sms_time') 336 for tndx=tlen0+itolap_ncep :tlen;337 nc_frc{'sms_time'}(tndx)=nc{'time'}(tndx-tlen0- (itolap_ncep-1));340 for tndx=tlen0+itolap_ncep+1:tlen; 341 nc_frc{'sms_time'}(tndx)=nc{'time'}(tndx-tlen0-itolap_ncep); 338 342 end; 339 343 end 340 344 % 341 345 if makeblk==1 342 346 disp('bulk_time') 343 for tndx=tlen0+itolap_ncep:tlen; 344 nc_blk{'bulk_time'}(tndx)=nc{'time'}(tndx-tlen0-(itolap_ncep-1)); 345 end; 346 end 347 close(nc) 348 end 349 % 350 % Perform the interpolations for the next month 351 % 352 disp('Last steps') 347 for tndx=tlen0+itolap_ncep+1:tlen; 348 nc_blk{'bulk_time'}(tndx)=nc{'time'}(tndx-tlen0-itolap_ncep); 349 end; 350 end 351 close(nc) 352 end 353 % 353 354 for tndx=tlen0+itolap_ncep+1:tlen; 354 355 355 disp(['tndx= ',num2str(tndx)]) 356 tout=tndx; 356 357 disp(['tout=tndx ',num2str(tndx)]) 357 358 358 if Mp==M 359 tin=tlen0; % persistency if current month is used 359 360 disp(['tin=',num2str(tin)]) 360 361 361 else 362 tin=tndx-tlen0-itolap_ncep; 362 363 disp(['tin=',num2str(tin)]) 363 364 365 mask,tin,nc_frc,nc_blk,lon,lat,angle,tout,Get_My_Data) 364 end 365 interp_NCEP(NCEP_dir,Yp,Mp,Roa,interp_method,lon1,lat1,... 366 mask,tin,nc_frc,nc_blk,lon,lat,angle,tout,Get_My_Data) 366 367 end; 367 368 % … … 369 370 % 370 371 if ~isempty(nc_frc) 371 372 close(nc_frc); 372 373 end 373 374 if ~isempty(nc_blk) 374 375 close(nc_blk); 375 376 end 376 377 end … … 382 383 % 383 384 disp('======================================================') 384 disp('Add spin up phase') 385 disp('Add spin up phase') 385 386 if SPIN_Long>0 386 387 M=Mmin-1; … … 389 390 M=M+1; 390 391 if M==13 391 M=1; 392 M=1; 392 393 Y=Y+1; 393 394 end … … 401 402 frcname=[frc_prefix,'Y',num2str(Ymin),'M',num2str(M),nc_suffix]; 402 403 frcname2=[frc_prefix,'Y',num2str(Y),'M',num2str(M),nc_suffix]; 403 disp(['Create ',frcname2]) 404 eval(['!cp ',frcname,' ',frcname2]) 404 disp(['Create ',frcname2]) 405 eval(['!cp ',frcname,' ',frcname2]) 405 406 % 406 407 % Change the time … … 424 425 blkname=[blk_prefix,'Y',num2str(Ymin),'M',num2str(M),nc_suffix]; 425 426 blkname2=[blk_prefix,'Y',num2str(Y),'M',num2str(M),nc_suffix]; 426 disp(['Create ',blkname2]) 427 eval(['!cp ',blkname,' ',blkname2]) 427 disp(['Create ',blkname2]) 428 eval(['!cp ',blkname,' ',blkname2]) 428 429 % 429 430 % Change the time … … 470 471 test_forcing(frcname,grdname,'sustr',slides,3,coastfileplot) 471 472 figure 472 test_forcing(frcname,grdname,'svstr',slides,3,coastfileplot) 473 end 474 end 475 473 test_forcing(frcname,grdname,'svstr',slides,3,coastfileplot) 474 end 475 end 476
Note: See TracChangeset
for help on using the changeset viewer.
