Changeset 48
- Timestamp:
- 03/16/14 20:38:39 (10 years ago)
- Location:
- trunk
- Files:
-
- 45 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/FIG_bar_transp_write.pro
r2 r48 1 PRO FIG_bar_transp_write,iyear,ian 1 PRO FIG_bar_transp_write,iyear,ian 2 2 @init2 3 3 @initorca2_bab … … 15 15 ;e_year0='10' 16 16 key_portrait = 0 17 ; stockage des fichiers brut 17 ; stockage des fichiers brut 18 18 ioDATA='/usr/work/sur/fvi/OPA/ORCA2/'+e_exp 19 file_U=e_exp+'_5d_'+ian+'0101_'+ian+'1231_grid_U.nc' 19 file_U=e_exp+'_5d_'+ian+'0101_'+ian+'1231_grid_U.nc' 20 20 print, file_U 21 file_V=e_exp+'_5d_'+ian+'0101_'+ian+'1231_grid_V.nc' 21 file_V=e_exp+'_5d_'+ian+'0101_'+ian+'1231_grid_V.nc' 22 22 print, file_V 23 23 … … 27 27 t_bt = 'bar_transp' 28 28 ioORLN2 = '/usr/work/sur/fvi/OPA/ORCA2' 29 ;facteur d'echelle vertical 29 ;facteur d'echelle vertical for partial steps 30 30 e3v3d=read_ncdf('e3v_ps',0,/timestep,iodir=ioORLN2,/nostruct,/tout,filename='meshmask_bab.nc') 31 31 e3u3d=read_ncdf('e3u_ps',0,/timestep,iodir=ioORLN2,/nostruct,/tout,filename='meshmask_bab.nc') … … 34 34 35 35 ; vertical integration: 36 e3t3d=make_array(jpi,jpj,jpk) 36 e3t3d=make_array(jpi,jpj,jpk) 37 37 for k=0, jpk-1 do begin &$ 38 for j=0,jpj-1 do begin &$ 38 for j=0,jpj-1 do begin &$ 39 39 for i=0,jpi-1 do begin &$ 40 40 e3t3d(i,j,k) = e3t(k) &$ … … 44 44 jpt = 73 45 45 46 ;vud = make_array(jpi,jpj,jpt) 46 ;vud = make_array(jpi,jpj,jpt) 47 47 ;vvd = make_array(jpi, jpj, jpt) 48 48 49 ; ouverture des fichiers dans lesquels on va écrire49 ; ouverture des fichiers dans lesquels on va écrire 50 50 id3=NCDF_OPEN('/usr/work/sur/fvi/geomag/U_5d_'+iyear+'_grid_T.nc',/write) 51 51 id4=NCDF_OPEN('/usr/work/sur/fvi/geomag/V_5d_'+iyear+'_grid_T.nc',/write) 52 52 53 FOR jt = 0, jpt-1 DO BEGIN &$ 53 FOR jt = 0, jpt-1 DO BEGIN &$ 54 54 55 55 ; ouverture des fichiers et stockage en memoire partial steps 56 vu=read_ncdf('vozocrtx',jt,jt, /timestep, iodir=ioDATA,/nostruct,/TOUT,filename=file_U) &$ 57 vv=read_ncdf('vomecrty',jt,jt, /timestep,iodir=ioDATA,/nostruct,/TOUT,filename=file_V) &$ 56 vu=read_ncdf('vozocrtx',jt,jt, /timestep, iodir=ioDATA,/nostruct,/TOUT,filename=file_U) &$ 57 vv=read_ncdf('vomecrty',jt,jt, /timestep,iodir=ioDATA,/nostruct,/TOUT,filename=file_V) &$ 58 58 59 59 60 60 61 61 ; Somme sur la verticale partial steps 62 vum=total( vu*e3u3d*umask(),3 ) &$ 63 vvm=total( vv*e3v3d*vmask(),3 ) &$ 62 vum=total( vu*e3u3d*umask(),3 ) &$ 63 vvm=total( vv*e3v3d*vmask(),3 ) &$ 64 64 65 65 ; Shit sur la grille T partial steps 66 vut= (vum+shift(vum,1,0) )*0.5 &$ 66 vut= (vum+shift(vum,1,0) )*0.5 &$ 67 67 vvt= (vvm+shift(vvm,0,1) )*0.5 &$ 68 ; Bande de recouvrement 68 ; Bande de recouvrement 69 69 vut(0, *) = vut(jpi-2, *) 70 70 vvt(*, 0) = 0. 71 ; stockage dans le fichier de sortie 71 ; stockage dans le fichier de sortie 72 72 NCDF_VARPUT, id3,'sossheig',vut, offset = [0, 0, jt] 73 73 NCDF_VARPUT, id4,'sossheig',vvt, offset = [0, 0, jt] … … 75 75 print, jt 76 76 77 ENDFOR 77 ENDFOR 78 78 ; on ferme le NetCDF 79 79 NCDF_CLOSE,id3 … … 107 107 ; tidout = NCDF_DIMDEF(idout, 'time_counter', /unlimited) 108 108 109 109 110 110 ; Attributs globaux 111 111 ; id0 = NCDF_VARDEF(idout, 'nav_lon' , [xidout, yidout ], /FLOAT) … … 141 141 ; NCDF_ATTPUT, idout, id5, 'title', 'vomecrty' 142 142 ; NCDF_ATTPUT, idout, id5, 'long_name', 'meridian current summmed on vertical' 143 143 144 144 ; NCDF_CONTROL, idout, /ENDEF 145 145 … … 148 148 ; Creation de la latitude 149 149 ; NCDF_VARPUT, idout, id1, gphit 150 ; Creation de la profondeur 150 ; Creation de la profondeur 151 151 ; prof=1. 152 ; NCDF_VARPUT, idout, id2, prof 152 ; NCDF_VARPUT, idout, id2, prof 153 153 ; Creation du calendrier 154 154 ; temps = findgen(jpt)*5.*86400. +julday(1, 1, 1) … … 158 158 ; NCDF_VARPUT, idout, id3, temps 159 159 160 160 161 161 ; Ecriture des donnees 162 162 -
trunk/INTERP2/check_inside.pro
r2 r48 93 93 94 94 result= where( ( (AB_Vect_BP LE 0.) AND (BC_Vect_CP LE 0.) AND (CD_Vect_DP LE 0.) AND (DA_Vect_AP LE 0.) ) OR $ 95 ( (AB_Vect_BP GE 0.) AND (BC_Vect_CP GE 0.) AND (CD_Vect_DP GE 0.) AND (DA_Vect_AP GE 0.) ) ) 95 ( (AB_Vect_BP GE 0.) AND (BC_Vect_CP GE 0.) AND (CD_Vect_DP GE 0.) AND (DA_Vect_AP GE 0.) ) ) 96 96 97 97 -
trunk/INTERP2/init_path2.pro
r2 r48 38 38 39 39 ; 40 ; 2. Input information 40 ; 2. Input information 41 41 ; ==================== 42 42 ; … … 44 44 ; ------------------- 45 45 ; 46 46 47 47 48 48 datafile = 'ForcageMOED_5d_1993_grid_T.nc' … … 61 61 62 62 ; 63 ; 3. Output information 63 ; 3. Output information 64 64 ; ===================== 65 65 ; … … 82 82 ; 83 83 ;;;; outputfile = 'DivBustar_5d_1998_reg.nc' 84 84 85 85 86 86 ; … … 90 90 91 91 errorfile='error_wT' 92 92 93 93 ;*************** END OF USER'S PART ******************* 94 94 ;********DO NOT MODIFY ANYTHING BELLOW THIS LINE******* -
trunk/INTERP2/initorca.pro
r2 r48 1 ;2 1 ; 3 2 @common -
trunk/INTERP2/initorca_grand.pro
r2 r48 1 ;2 ;3 1 @common 4 2 @init_path2 -
trunk/INTERP2/initorca_reg.pro
r2 r48 1 ;2 ;3 1 @common 4 2 @init_path2 -
trunk/INTERP2/interp2.pro
r2 r48 13 13 ; 14 14 ; INPUTS : 15 ; A full ORCA grid file + a 1D or 2D field on this grid 15 ; A full ORCA grid file + a 1D or 2D field on this grid 16 16 ; OUTPUTS : 17 17 ; - A grid file for the geographical equivalent to the given 18 ; ORCA grid 18 ; ORCA grid 19 19 ; - The interpolated field on this new grid + the weights and 20 20 ; pointers required to creat it. … … 75 75 76 76 end 77 78 -
trunk/INTERP2/make_geogrid.pro
r2 r48 53 53 ; Output values are set to input values 54 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 55 glamt_reg=glamt 56 glamu_reg=glamu 57 glamv_reg=glamv 58 glamf_reg=glamf 59 60 gphit_reg=gphit 61 gphiu_reg=gphiu 62 gphiv_reg=gphiv 63 gphif_reg=gphif 64 65 e1t_reg=e1t 66 e1u_reg=e1u 67 e1v_reg=e1v 68 e1f_reg=e1f 69 70 e2t_reg=e2t 71 e2u_reg=e2u 72 e2v_reg=e2v 73 e2f_reg=e2f 74 75 tmask_reg=tmask 76 umask_reg=tmask 77 vmask_reg=tmask 78 fmask_reg=tmask 79 79 80 80 … … 96 96 ; Let's get the reference longitude 97 97 98 equador=where(abs(gphit) EQ 98 equador=where(abs(gphit) EQ min(abs(gphit)) ) 99 99 100 100 ref_longitude=glamt(equador) … … 171 171 172 172 ;************************************************ 173 ; 2nd STEP : Define the mask of the outp ût grid173 ; 2nd STEP : Define the mask of the outpût grid 174 174 ;************************************************ 175 175 … … 188 188 ; to avoid any problem 189 189 190 sph_coord(2,jpj*jpi-1-jpi/2:jpj*jpi-1)=10. 190 sph_coord(2,jpj*jpi-1-jpi/2:jpj*jpi-1)=10. 191 191 192 192 sph_coord(0,*)=reform(glamt,jpi*jpj) … … 196 196 197 197 xx_grid=rec_grid(0,*) 198 yy_grid=rec_grid(1,*) 198 yy_grid=rec_grid(1,*) 199 199 zz_grid=rec_grid(2,*) 200 200 201 201 202 202 not_found=1 203 203 for jj=0, jpj-1 do begin 204 204 205 206 207 208 209 205 lat_max_ORCA=max( gphit(*,jj) ) 206 lat_max_reg =max( gphit_reg(*,jj) ) 207 208 209 if ((lat_max_ORCA NE lat_max_reg) AND (not_found EQ 1)) then begin 210 210 jpj_start=jj-2 211 211 not_found=0 212 212 endif 213 213 214 214 end … … 220 220 for jj=jpj/2,jpj-1 do begin 221 221 for ji=0,jpi-1 do begin 222 222 223 223 224 224 ; Coordinates of our point … … 226 226 xx=glamt_reg(ji,jj) 227 227 yy=gphit_reg(ji,jj) 228 228 229 229 sph_coord=replicate(1.,3) 230 230 sph_coord(0)=xx 231 231 sph_coord(1)=yy 232 232 233 ; Convert them into cartesian 234 233 ; Convert them into cartesian 234 235 235 rec_coord=cv_coord(/DEGREES,FROM_SPHERE=sph_coord,/TO_RECT) 236 236 … … 242 242 +(yy_grid-yy_bis)*(yy_grid-yy_bis) $ 243 243 +(zz_grid-zz_bis)*(zz_grid-zz_bis)) 244 244 245 245 246 246 sort_mod=where(modules EQ min(modules)) 247 247 248 answj=sort_mod(0)/jpi 248 answj=sort_mod(0)/jpi 249 249 250 250 answi=sort_mod(0)-answj*jpi … … 252 252 253 253 mask_reg(ji,jj,*)=tmask(answi,answj,*) 254 255 254 255 256 256 endfor 257 257 endfor … … 265 265 ; Create Global Mask file for T point 266 266 267 268 269 267 tmask_reg_glo=replicate(0.,jpiglo,jpjglo,jpk) 268 269 tmask_reg_glo(1:jpiglo-2,0:jpjglo-2, *)=mask_reg 270 270 271 271 … … 274 274 ; North 275 275 276 276 tmask_reg_glo(*,jpjglo-1,*)=tmask_reg_glo(*,jpjglo-2,*) 277 277 ;E/W 278 278 279 280 279 tmask_reg_glo(0,*,*)=tmask_reg_glo(jpiglo-2,*,*) 280 tmask_reg_glo(jpiglo-1,*,*)=tmask_reg_glo(1,*,*) 281 281 282 282 ; Global T mask in created … … 285 285 286 286 287 288 289 290 291 292 287 umask_reg=tmask_reg_glo(1:jpiglo-2,0:jpjglo-2,*)*tmask_reg_glo(2:jpiglo-1,0:jpjglo-2,*) 288 289 vmask_reg=tmask_reg_glo(1:jpiglo-2,0:jpjglo-2,*)*tmask_reg_glo(1:jpiglo-2,1:jpjglo-1,*) 290 291 fmask_reg=tmask_reg_glo(1:jpiglo-2,0:jpjglo-2,*)*tmask_reg_glo(2:jpiglo-1,0:jpjglo-2,*) $ 292 *tmask_reg_glo(1:jpiglo-2,1:jpjglo-1,*)*tmask_reg_glo(2:jpiglo-1,1:jpjglo-1,*) 293 293 294 294 ; Create Global Mask file for U,V & F points 295 295 296 297 298 299 300 301 302 303 296 umask_reg_glo=replicate(0.,jpiglo,jpjglo,jpk) 297 umask_reg_glo(1:jpiglo-2,0:jpjglo-2, *)=umask_reg 298 299 vmask_reg_glo=replicate(0.,jpiglo,jpjglo,jpk) 300 vmask_reg_glo(1:jpiglo-2,0:jpjglo-2, *)=vmask_reg 301 302 fmask_reg_glo=replicate(0.,jpiglo,jpjglo,jpk) 303 fmask_reg_glo(1:jpiglo-2,0:jpjglo-2, *)=fmask_reg 304 304 305 305 ; Data Duplication … … 307 307 ; North 308 308 309 310 311 309 umask_reg_glo(*,jpjglo-1,*)=umask_reg_glo(*,jpjglo-2,*) 310 vmask_reg_glo(*,jpjglo-1,*)=vmask_reg_glo(*,jpjglo-2,*) 311 fmask_reg_glo(*,jpjglo-1,*)=fmask_reg_glo(*,jpjglo-2,*) 312 312 313 313 ;E/W 314 314 315 316 317 318 319 320 321 322 315 umask_reg_glo(0,*,*)=umask_reg_glo(jpiglo-2,*,*) 316 umask_reg_glo(jpiglo-1,*,*)=umask_reg_glo(1,*,*) 317 318 vmask_reg_glo(0,*,*)=vmask_reg_glo(jpiglo-2,*,*) 319 vmask_reg_glo(jpiglo-1,*,*)=vmask_reg_glo(1,*,*) 320 321 fmask_reg_glo(0,*,*)=fmask_reg_glo(jpiglo-2,*,*) 322 fmask_reg_glo(jpiglo-1,*,*)=fmask_reg_glo(1,*,*) 323 323 324 324 ; ***** NAV_LON & NAV_LAT ********** … … 327 327 ; Create Global longitude and latitude 328 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 329 glamt_reg_glo=replicate(0.,jpiglo,jpjglo) 330 331 glamt_reg_glo(1:jpiglo-2,0:jpjglo-2)=glamt_reg 332 333 gphit_reg_glo=replicate(0.,jpiglo,jpjglo) 334 335 gphit_reg_glo(1:jpiglo-2,0:jpjglo-2)=gphit_reg 336 337 nav_lon_reg_glo=replicate(0.,jpiglo,jpjglo) 338 339 nav_lat_reg_glo=replicate(0.,jpiglo,jpjglo) 340 341 nav_lon_reg_glo=glamt_reg_glo 342 343 nav_lat_reg_glo=gphit_reg_glo 344 344 345 345 ; Data Duplication … … 347 347 ; E/W Transition 348 348 349 350 351 352 353 354 355 ; North 356 357 358 349 nav_lon_reg_glo(0,*)=nav_lon_reg_glo(jpiglo-2,*) 350 nav_lon_reg_glo(jpiglo-1,*)=nav_lon_reg_glo(1,*) 351 352 nav_lat_reg_glo(0,*)=nav_lat_reg_glo(jpiglo-2,*) 353 nav_lat_reg_glo(jpiglo-1,*)=nav_lat_reg_glo(1,*) 354 355 ; North 356 357 nav_lon_reg_glo(*,jpjglo-1)= nav_lon_reg_glo(*,jpjglo-2) 358 nav_lat_reg_glo(*,jpjglo-1)= nav_lat_reg_glo(*,jpjglo-2) 359 359 360 360 ; ***** GLAMT ***** 361 ; We don't do it to save memory : same as nav_lon 361 ; We don't do it to save memory : same as nav_lon 362 362 363 363 ; ***** GLAMU ***** 364 364 365 366 367 368 ; North 369 370 371 372 ; E/W 373 374 375 376 365 glamu_reg_glo=replicate(0.,jpiglo,jpjglo) 366 glamu_reg_glo(1:jpiglo-2,0:jpjglo-2)=glamu_reg 367 368 ; North 369 370 glamu_reg_glo(*,jpjglo-1)=glamu_reg_glo(*,jpjglo-2) 371 372 ; E/W 373 374 glamu_reg_glo(0,*)=glamu_reg_glo(jpiglo-2,*) 375 glamu_reg_glo(jpiglo-1,*)=glamu_reg_glo(1,*) 376 377 377 ; ***** GLAMV ***** 378 378 379 380 381 382 ; North 383 384 385 386 ; E/W 387 388 389 379 glamv_reg_glo=replicate(0.,jpiglo,jpjglo) 380 glamv_reg_glo(1:jpiglo-2,0:jpjglo-2)=glamv_reg 381 382 ; North 383 384 glamv_reg_glo(*,jpjglo-1)=glamv_reg_glo(*,jpjglo-2) 385 386 ; E/W 387 388 glamv_reg_glo(0,*)=glamv_reg_glo(jpiglo-2,*) 389 glamv_reg_glo(jpiglo-1,*)=glamv_reg_glo(1,*) 390 390 391 391 392 392 ; ***** GLAMF ***** 393 393 394 395 396 397 ; North 398 399 400 401 ; E/W 402 403 404 394 glamf_reg_glo=replicate(0.,jpiglo,jpjglo) 395 glamf_reg_glo(1:jpiglo-2,0:jpjglo-2)=glamf_reg 396 397 ; North 398 399 glamf_reg_glo(*,jpjglo-1)=glamf_reg_glo(*,jpjglo-2) 400 401 ; E/W 402 403 glamf_reg_glo(0,*)=glamf_reg_glo(jpiglo-2,*) 404 glamf_reg_glo(jpiglo-1,*)=glamf_reg_glo(1,*) 405 405 406 406 407 407 ; ***** GPHIT ***** 408 ; We don't do it to save memory : same as nav_lat 408 ; We don't do it to save memory : same as nav_lat 409 409 410 410 ; ***** GPHIU ***** 411 411 412 413 414 415 ; North 416 417 418 419 ; E/W 420 421 422 423 412 gphiu_reg_glo=replicate(0.,jpiglo,jpjglo) 413 gphiu_reg_glo(1:jpiglo-2,0:jpjglo-2)=gphiu_reg 414 415 ; North 416 417 gphiu_reg_glo(*,jpjglo-1)=gphiu_reg_glo(*,jpjglo-2) 418 419 ; E/W 420 421 gphiu_reg_glo(0,*)=gphiu_reg_glo(jpiglo-2,*) 422 gphiu_reg_glo(jpiglo-1,*)=gphiu_reg_glo(1,*) 423 424 424 ; ***** GPHIV ***** 425 425 426 427 428 429 ; North 430 431 432 433 ; E/W 434 435 436 426 gphiv_reg_glo=replicate(0.,jpiglo,jpjglo) 427 gphiv_reg_glo(1:jpiglo-2,0:jpjglo-2)=gphiv_reg 428 429 ; North 430 431 gphiv_reg_glo(*,jpjglo-1)=gphiv_reg_glo(*,jpjglo-2) 432 433 ; E/W 434 435 gphiv_reg_glo(0,*)=gphiv_reg_glo(jpiglo-2,*) 436 gphiv_reg_glo(jpiglo-1,*)=gphiv_reg_glo(1,*) 437 437 438 438 439 439 ; ***** GPHIF ***** 440 440 441 442 443 444 ; North 445 446 447 448 ; E/W 449 450 451 441 gphif_reg_glo=replicate(0.,jpiglo,jpjglo) 442 gphif_reg_glo(1:jpiglo-2,0:jpjglo-2)=gphif_reg 443 444 ; North 445 446 gphif_reg_glo(*,jpjglo-1)=gphif_reg_glo(*,jpjglo-2) 447 448 ; E/W 449 450 gphif_reg_glo(0,*)=gphif_reg_glo(jpiglo-2,*) 451 gphif_reg_glo(jpiglo-1,*)=gphif_reg_glo(1,*) 452 452 453 453 ; ***** FF ******* 454 454 455 456 457 ff_reg_glo = 2.*omega*sin(gphif_reg_glo*!PI/180.) 458 455 omega = 2.0 * !PI / 86164.0 456 457 ff_reg_glo = 2.*omega*sin(gphif_reg_glo*!PI/180.) 458 459 459 460 460 ; *********** E1T*************** 461 462 463 464 465 ; North 466 467 468 469 ; E/W 470 471 472 461 462 e1t_reg_glo=replicate(0.,jpiglo,jpjglo) 463 e1t_reg_glo(1:jpiglo-2,0:jpjglo-2)=e1t_reg 464 465 ; North 466 467 e1t_reg_glo(*,jpjglo-1)=e1t_reg_glo(*,jpjglo-2) 468 469 ; E/W 470 471 e1t_reg_glo(0,*)=e1t_reg_glo(jpiglo-2,*) 472 e1t_reg_glo(jpiglo-1,*)=e1t_reg_glo(1,*) 473 473 474 474 475 475 476 476 ; *********** E2T*************** 477 478 479 480 481 ; North 482 483 484 485 ; E/W 486 487 488 477 478 e2t_reg_glo=replicate(0.,jpiglo,jpjglo) 479 e2t_reg_glo(1:jpiglo-2,0:jpjglo-2)=e2t_reg 480 481 ; North 482 483 e2t_reg_glo(*,jpjglo-1)=e2t_reg_glo(*,jpjglo-2) 484 485 ; E/W 486 487 e2t_reg_glo(0,*)=e2t_reg_glo(jpiglo-2,*) 488 e2t_reg_glo(jpiglo-1,*)=e2t_reg_glo(1,*) 489 489 490 490 ; *********** E1U*************** 491 492 493 494 495 ; North 496 497 498 499 ; E/W 500 501 502 491 492 e1u_reg_glo=replicate(0.,jpiglo,jpjglo) 493 e1u_reg_glo(1:jpiglo-2,0:jpjglo-2)=e1u_reg 494 495 ; North 496 497 e1u_reg_glo(*,jpjglo-1)=e1u_reg_glo(*,jpjglo-2) 498 499 ; E/W 500 501 e1u_reg_glo(0,*)=e1u_reg_glo(jpiglo-2,*) 502 e1u_reg_glo(jpiglo-1,*)=e1u_reg_glo(1,*) 503 503 504 504 505 505 ; *********** E2U*************** 506 507 508 509 510 ; North 511 512 513 514 ; E/W 515 516 517 506 507 e2u_reg_glo=replicate(0.,jpiglo,jpjglo) 508 e2u_reg_glo(1:jpiglo-2,0:jpjglo-2)=e2u_reg 509 510 ; North 511 512 e2u_reg_glo(*,jpjglo-1)=e2u_reg_glo(*,jpjglo-2) 513 514 ; E/W 515 516 e2u_reg_glo(0,*)=e2u_reg_glo(jpiglo-2,*) 517 e2u_reg_glo(jpiglo-1,*)=e2u_reg_glo(1,*) 518 518 519 519 ; *********** E1V*************** 520 521 522 523 524 ; North 525 526 527 528 ; E/W 529 530 531 520 521 e1v_reg_glo=replicate(0.,jpiglo,jpjglo) 522 e1v_reg_glo(1:jpiglo-2,0:jpjglo-2)=e1v_reg 523 524 ; North 525 526 e1v_reg_glo(*,jpjglo-1)=e1v_reg_glo(*,jpjglo-2) 527 528 ; E/W 529 530 e1v_reg_glo(0,*)=e1v_reg_glo(jpiglo-2,*) 531 e1v_reg_glo(jpiglo-1,*)=e1v_reg_glo(1,*) 532 532 533 533 ; *********** E2V*************** 534 535 536 537 538 ; North 539 540 541 542 ; E/W 543 544 545 534 535 e2v_reg_glo=replicate(0.,jpiglo,jpjglo) 536 e2v_reg_glo(1:jpiglo-2,0:jpjglo-2)=e2v_reg 537 538 ; North 539 540 e2v_reg_glo(*,jpjglo-1)=e2v_reg_glo(*,jpjglo-2) 541 542 ; E/W 543 544 e2v_reg_glo(0,*)=e2v_reg_glo(jpiglo-2,*) 545 e2v_reg_glo(jpiglo-1,*)=e2v_reg_glo(1,*) 546 546 547 547 ; *********** E1F ************** 548 548 549 550 551 552 ; North 553 554 555 556 ; E/W 557 558 559 549 e1f_reg_glo=replicate(0.,jpiglo,jpjglo) 550 e1f_reg_glo(1:jpiglo-2,0:jpjglo-2)=e1f_reg 551 552 ; North 553 554 e1f_reg_glo(*,jpjglo-1)=e1f_reg_glo(*,jpjglo-2) 555 556 ; E/W 557 558 e1f_reg_glo(0,*)=e1f_reg_glo(jpiglo-2,*) 559 e1f_reg_glo(jpiglo-1,*)=e1f_reg_glo(1,*) 560 560 561 561 ; *********** E2F*************** 562 563 564 565 566 ; North 567 568 569 570 ; E/W 571 572 573 562 563 e2f_reg_glo=replicate(0.,jpiglo,jpjglo) 564 e2f_reg_glo(1:jpiglo-2,0:jpjglo-2)=e2f_reg 565 566 ; North 567 568 e2f_reg_glo(*,jpjglo-1)=e2f_reg_glo(*,jpjglo-2) 569 570 ; E/W 571 572 e2f_reg_glo(0,*)=e2f_reg_glo(jpiglo-2,*) 573 e2f_reg_glo(jpiglo-1,*)=e2f_reg_glo(1,*) 574 574 575 575 … … 577 577 ;****************END OF ATTRIBUTES**************** 578 578 579 579 580 580 ;************************************************* 581 581 ; 4th STEP : CREATES NCDF MESHMASK FILE … … 585 585 586 586 587 587 vargrid = 'T' 588 588 589 589 ; Name 590 591 592 590 idout = NCDF_CREATE(outputgrid,/clobber) 591 print, 'Creation du fichier Netcdf' 592 NCDF_CONTROL, idout, /nofill 593 593 ; Dimension 594 594 595 596 597 598 599 600 601 595 xidout = NCDF_DIMDEF(idout, 'x', jpiglo) 596 yidout = NCDF_DIMDEF(idout, 'y', jpjglo) 597 didout = NCDF_DIMDEF(idout, 'z', jpk) 598 x_a = NCDF_DIMDEF(idout, 'x_a',1) 599 y_a = NCDF_DIMDEF(idout, 'y_a',1) 600 z_a = NCDF_DIMDEF(idout, 'z_a',1) 601 tidout = NCDF_DIMDEF(idout, 'time_counter', /unlimited) 602 602 603 603 ; Attributes 604 605 606 607 608 609 610 611 612 613 614 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629 630 631 632 633 634 604 605 id0 = NCDF_VARDEF(idout, 'nav_lon' , [xidout, yidout ],/FLOAT ) 606 id1 = NCDF_VARDEF(idout, 'nav_lat' , [xidout, yidout ],/FLOAT ) 607 id2 = NCDF_VARDEF(idout, 'nav_lev' , [ didout ],/FLOAT ) 608 id3 = NCDF_VARDEF(idout, 'time' , [ tidout],/FLOAT ) 609 id4 = NCDF_VARDEF(idout, 'time_steps' , [ tidout],/LONG ) 610 id5 = NCDF_VARDEF(idout, 'glamt' , [xidout, yidout, z_a , tidout],/DOUBLE) 611 id6 = NCDF_VARDEF(idout, 'glamu' , [xidout, yidout, z_a , tidout],/DOUBLE) 612 id7 = NCDF_VARDEF(idout, 'glamv' , [xidout, yidout, z_a , tidout],/DOUBLE) 613 id8 = NCDF_VARDEF(idout, 'glamf' , [xidout, yidout, z_a , tidout],/DOUBLE) 614 id9 = NCDF_VARDEF(idout, 'gphit' , [xidout, yidout, z_a , tidout],/DOUBLE) 615 id10 = NCDF_VARDEF(idout, 'gphiu' , [xidout, yidout, z_a , tidout],/DOUBLE) 616 id11 = NCDF_VARDEF(idout, 'gphiv' , [xidout, yidout, z_a , tidout],/DOUBLE) 617 id12 = NCDF_VARDEF(idout, 'gphif' , [xidout, yidout, z_a , tidout],/DOUBLE) 618 id13 = NCDF_VARDEF(idout, 'e1t' , [xidout, yidout, z_a , tidout],/DOUBLE) 619 id14 = NCDF_VARDEF(idout, 'e1u' , [xidout, yidout, z_a , tidout],/DOUBLE) 620 id15 = NCDF_VARDEF(idout, 'e1v' , [xidout, yidout, z_a , tidout],/DOUBLE) 621 id16 = NCDF_VARDEF(idout, 'e1f' , [xidout, yidout, z_a , tidout],/DOUBLE) 622 id17 = NCDF_VARDEF(idout, 'e2t' , [xidout, yidout, z_a , tidout],/DOUBLE) 623 id18 = NCDF_VARDEF(idout, 'e2u' , [xidout, yidout, z_a , tidout],/DOUBLE) 624 id19 = NCDF_VARDEF(idout, 'e2v' , [xidout, yidout, z_a , tidout],/DOUBLE) 625 id20 = NCDF_VARDEF(idout, 'e2f' , [xidout, yidout, z_a , tidout],/DOUBLE) 626 id21 = NCDF_VARDEF(idout, 'tmask' , [xidout, yidout, didout, tidout],/DOUBLE) 627 id22 = NCDF_VARDEF(idout, 'umask' , [xidout, yidout, didout, tidout],/DOUBLE) 628 id23 = NCDF_VARDEF(idout, 'vmask' , [xidout, yidout, didout, tidout],/DOUBLE) 629 id24 = NCDF_VARDEF(idout, 'fmask' , [xidout, yidout, didout, tidout],/DOUBLE) 630 id25 = NCDF_VARDEF(idout, 'ff' , [xidout, yidout, z_a , tidout],/DOUBLE) 631 id26 = NCDF_VARDEF(idout, 'gdept' , [x_a , y_a , didout, tidout],/DOUBLE) 632 id27 = NCDF_VARDEF(idout, 'gdepw' , [x_a , y_a , didout, tidout],/DOUBLE) 633 id28 = NCDF_VARDEF(idout, 'e3t' , [x_a , y_a , didout, tidout],/DOUBLE) 634 id29 = NCDF_VARDEF(idout, 'e3w' , [x_a , y_a , didout, tidout],/DOUBLE) 635 635 636 636 ; Variable 0 637 638 639 637 NCDF_ATTPUT, idout, id0, 'units', 'degrees_east' 638 NCDF_ATTPUT, idout, id0, 'long_name', 'Longitude' 639 NCDF_ATTPUT, idout, id0, 'nav_model', 'Default grid' 640 640 ; Variable 1 641 642 643 641 NCDF_ATTPUT, idout, id1, 'units', 'degrees_north' 642 NCDF_ATTPUT, idout, id1, 'long_name', 'Latitude' 643 NCDF_ATTPUT, idout, id1, 'nav_model', 'Default grid' 644 644 ; Variable 2 645 646 647 645 NCDF_ATTPUT, idout, id2, 'units','meters' 646 NCDF_ATTPUT, idout, id2, 'long_name','Depth' 647 NCDF_ATTPUT, idout, id2, 'nav_model','Default grid' 648 648 ; Variable 3 649 650 651 652 653 654 655 656 649 NCDF_ATTPUT, idout, id3, 'units', 'seconds since 0001-01-15 12:00:00 ' 650 NCDF_ATTPUT, idout, id3, 'calendar','noleap' 651 NCDF_ATTPUT, idout, id3, 'title', 'Time' 652 NCDF_ATTPUT, idout, id3, 'long_name', 'Time axis' 653 NCDF_ATTPUT, idout, id3, 'time_origin','0000-DEC-15 00:00:00' 654 655 656 NCDF_CONTROL, idout, /ENDEF 657 657 658 658 … … 661 661 662 662 663 664 665 666 667 668 669 670 671 672 673 674 675 676 677 678 679 680 681 682 683 684 685 686 687 688 689 690 691 692 693 694 695 696 697 698 699 700 701 702 703 704 705 706 707 708 709 710 711 712 713 663 NCDF_VARPUT, idout, id0 , nav_lon_reg_glo 664 NCDF_VARPUT, idout, id1 , nav_lat_reg_glo 665 NCDF_VARPUT, idout, id2 , gdept 666 NCDF_VARPUT, idout, id3 , 1. 667 NCDF_VARPUT, idout, id4 , 1. 668 NCDF_VARPUT, idout, id5 , nav_lon_reg_glo 669 670 NCDF_VARPUT, idout, id6 , glamu_reg_glo 671 NCDF_VARPUT, idout, id7 , glamv_reg_glo 672 NCDF_VARPUT, idout, id8 , glamf_reg_glo 673 NCDF_VARPUT, idout, id9 , nav_lat_reg_glo 674 NCDF_VARPUT, idout, id10, gphiu_reg_glo 675 NCDF_VARPUT, idout, id11, gphiv_reg_glo 676 NCDF_VARPUT, idout, id12, gphif_reg_glo 677 NCDF_VARPUT, idout, id13, e1t_reg_glo 678 NCDF_VARPUT, idout, id14, e1u_reg_glo 679 NCDF_VARPUT, idout, id15, e1v_reg_glo 680 NCDF_VARPUT, idout, id16, e1f_reg_glo 681 NCDF_VARPUT, idout, id17, e2t_reg_glo 682 NCDF_VARPUT, idout, id18, e2u_reg_glo 683 NCDF_VARPUT, idout, id19, e2v_reg_glo 684 NCDF_VARPUT, idout, id20, e2f_reg_glo 685 NCDF_VARPUT, idout, id21, tmask_reg_glo 686 NCDF_VARPUT, idout, id22, umask_reg_glo 687 NCDF_VARPUT, idout, id23, vmask_reg_glo 688 NCDF_VARPUT, idout, id24, fmask_reg_glo 689 NCDF_VARPUT, idout, id25, ff_reg_glo 690 691 dummy=replicate(0.,1,1,jpk,1) 692 dummy(0,0,*,0)=gdept 693 694 NCDF_VARPUT, idout, id26, dummy 695 696 dummy=replicate(0.,1,1,jpk,1) 697 dummy(0,0,*,0)=gdepw 698 699 NCDF_VARPUT, idout, id27, dummy 700 701 dummy=replicate(0.,1,1,jpk,1) 702 dummy(0,0,*,0)=e3t 703 704 NCDF_VARPUT, idout, id28, dummy 705 706 dummy=replicate(0.,1,1,jpk,1) 707 dummy(0,0,*,0)=e3w 708 709 NCDF_VARPUT, idout, id29, dummy 710 711 NCDF_CLOSE, idout 712 713 714 714 715 715 end -
trunk/INTERP2/make_weight.pro
r2 r48 6 6 ; 7 7 ; PURPOSE: Creates for a given ORCA grid, and its ORCA_GEO equivalent 8 ; grid the required8 ; grid the required 9 9 ; 10 10 ; CATEGORY : Subroutine … … 36 36 37 37 38 @initorca39 @naminterp240 @init_path241 42 43 ;*****************************************44 ; Error file opening45 ;*****************************************46 47 openu, lun_err,errorfile,/get_lun48 49 ;************************50 ; Compute Geogrid51 ;************************52 ; Let's get our reference latitude53 54 northernmost=where(gphit EQ max(gphit))55 56 jnorth=northernmost/jpi57 58 inorth=min(northernmost-jnorth*jpi)59 60 ref_latitude=gphit(inorth,*)61 62 ; Let's get the reference longitude63 64 equador=where(abs(gphit) EQ min(abs(gphit)) )65 66 ref_longitude=glamt(equador)67 68 ; Now we have the new grid69 70 71 glamt_reg=ref_longitude#replicate(1,n_elements(ref_latitude))72 73 74 gphit_reg=replicate(1,n_elements(ref_longitude))#ref_latitude75 76 ; Get the output mask77 78 mask_reg=read_ncdf('tmask',FILENAME=outputgrid,/NOSTRUCT,/cont_nofill)79 80 ; Get the latitude and longitude with extra boundary lines81 82 get_fullmesh, glamt_full,gphit_full83 84 sph_coord_full=replicate(9999.,3,jpiglo*jpjglo)85 86 sph_coord_full(0,*)=reform(glamt_full,jpiglo*jpjglo)87 sph_coord_full(1,*)=reform(gphit_full,jpiglo*jpjglo)88 sph_coord_full(2,*)=replicate(1.,jpiglo*jpjglo)89 90 rec_grid_full=cv_coord(/DEGREES,FROM_SPHERE=sph_coord_full,/TO_RECT)91 92 xx_grid_full=rec_grid_full(0,*)93 yy_grid_full=rec_grid_full(1,*) 94 zz_grid_full=rec_grid_full(2,*)95 96 97 sph_coord=replicate(9999.,3,jpi*jpj)98 99 radius=replicate(1.,jpi,jpj)100 radius(0:jpi/2,jpj-1)=1000.101 102 sph_coord(0,*)=reform(glamt,jpi*jpj)103 sph_coord(1,*)=reform(gphit,jpi*jpj)104 sph_coord(2,*)=reform(radius,jpi*jpj)105 106 rec_grid=cv_coord(/DEGREES,FROM_SPHERE=sph_coord,/TO_RECT)107 108 xx_grid=rec_grid(0,*)109 yy_grid=rec_grid(1,*) 110 zz_grid=rec_grid(2,*)111 112 113 114 ; Init array to stock pointers115 116 pointer_Tx=replicate(999,jpi,jpj,4)117 pointer_Ty=replicate(999,jpi,jpj,4)38 @initorca 39 @naminterp2 40 @init_path2 41 42 43 ;***************************************** 44 ; Error file opening 45 ;***************************************** 46 47 openu, lun_err,errorfile,/get_lun 48 49 ;************************ 50 ; Compute Geogrid 51 ;************************ 52 ; Let's get our reference latitude 53 54 northernmost=where(gphit EQ max(gphit)) 55 56 jnorth=northernmost/jpi 57 58 inorth=min(northernmost-jnorth*jpi) 59 60 ref_latitude=gphit(inorth,*) 61 62 ; Let's get the reference longitude 63 64 equador=where(abs(gphit) EQ min(abs(gphit)) ) 65 66 ref_longitude=glamt(equador) 67 68 ; Now we have the new grid 69 70 71 glamt_reg=ref_longitude#replicate(1,n_elements(ref_latitude)) 72 73 74 gphit_reg=replicate(1,n_elements(ref_longitude))#ref_latitude 75 76 ; Get the output mask 77 78 mask_reg=read_ncdf('tmask',FILENAME=outputgrid,/NOSTRUCT,/cont_nofill) 79 80 ; Get the latitude and longitude with extra boundary lines 81 82 get_fullmesh, glamt_full,gphit_full 83 84 sph_coord_full=replicate(9999.,3,jpiglo*jpjglo) 85 86 sph_coord_full(0,*)=reform(glamt_full,jpiglo*jpjglo) 87 sph_coord_full(1,*)=reform(gphit_full,jpiglo*jpjglo) 88 sph_coord_full(2,*)=replicate(1.,jpiglo*jpjglo) 89 90 rec_grid_full=cv_coord(/DEGREES,FROM_SPHERE=sph_coord_full,/TO_RECT) 91 92 xx_grid_full=rec_grid_full(0,*) 93 yy_grid_full=rec_grid_full(1,*) 94 zz_grid_full=rec_grid_full(2,*) 95 96 97 sph_coord=replicate(9999.,3,jpi*jpj) 98 99 radius=replicate(1.,jpi,jpj) 100 radius(0:jpi/2,jpj-1)=1000. 101 102 sph_coord(0,*)=reform(glamt,jpi*jpj) 103 sph_coord(1,*)=reform(gphit,jpi*jpj) 104 sph_coord(2,*)=reform(radius,jpi*jpj) 105 106 rec_grid=cv_coord(/DEGREES,FROM_SPHERE=sph_coord,/TO_RECT) 107 108 xx_grid=rec_grid(0,*) 109 yy_grid=rec_grid(1,*) 110 zz_grid=rec_grid(2,*) 111 112 113 114 ; Init array to stock pointers 115 116 pointer_Tx=replicate(999,jpi,jpj,4) 117 pointer_Ty=replicate(999,jpi,jpj,4) 118 118 ; Init weight matric 119 119 120 weight_T=double(replicate(0.,jpi,jpj,4))120 weight_T=double(replicate(0.,jpi,jpj,4)) 121 121 122 122 ; Get the latitude to start from 123 123 124 not_found=1125 for jj=0, jpj-1 do begin126 127 lat_max_ORCA=max( gphit(*,jj) )128 lat_max_reg =max( gphit_reg(*,jj) )129 130 131 if ((lat_max_ORCA NE lat_max_reg) AND (not_found EQ 1)) then begin132 jpj_start=jj-2133 not_found=0134 endif135 136 end124 not_found=1 125 for jj=0, jpj-1 do begin 126 127 lat_max_ORCA=max( gphit(*,jj) ) 128 lat_max_reg =max( gphit_reg(*,jj) ) 129 130 131 if ((lat_max_ORCA NE lat_max_reg) AND (not_found EQ 1)) then begin 132 jpj_start=jj-2 133 not_found=0 134 endif 135 136 end 137 137 138 138 139 139 for jj=0,jpj-1 do begin 140 141 for ji=0,jpi-1 do begin142 143 144 if mask_reg(ji,jj) NE 0 then begin145 146 ; Coordinates of our point147 148 xx=glamt_reg(ji,jj)149 yy=gphit_reg(ji,jj)150 140 141 for ji=0,jpi-1 do begin 142 143 144 if mask_reg(ji,jj) NE 0 then begin 145 146 ; Coordinates of our point 147 148 xx=glamt_reg(ji,jj) 149 yy=gphit_reg(ji,jj) 150 151 151 sph_coord=replicate(1.,3) 152 152 sph_coord(0)=xx 153 153 sph_coord(1)=yy 154 154 155 ; Convert them into cartesian 156 155 ; Convert them into cartesian 156 157 157 rec_coord=cv_coord(/DEGREES,FROM_SPHERE=sph_coord,/TO_RECT) 158 158 … … 162 162 163 163 164 164 165 165 166 166 ;**************************************************************** 167 ; Finds in which grid cell point stand 167 ; Finds in which grid cell point stand 168 168 ;**************************************************************** 169 169 170 170 xx_bis=reform(xx_bis,jpi,jpj) 171 171 yy_bis=reform(yy_bis,jpi,jpj) … … 176 176 zz_grid_full=reform(zz_grid_full,jpiglo,jpjglo) 177 177 178 178 179 179 check_inside,xx_bis ,yy_bis ,zz_bis,$ 180 180 xx_grid_full(1:jpiglo-2,0:jpjglo-2),yy_grid_full(1:jpiglo-2,0:jpjglo-2),zz_grid_full(1:jpiglo-2,0:jpjglo-2),$ … … 183 183 xx_grid_full(1:jpiglo-2,1:jpjglo-1),yy_grid_full(1:jpiglo-2,1:jpjglo-1),zz_grid_full(1:jpiglo-2,1:jpjglo-1), answer 184 184 185 186 185 186 187 187 if answer(0) EQ -1 then begin 188 188 … … 191 191 192 192 endif 193 194 193 194 195 195 answer_j=answer/jpi 196 196 answer_i=answer-answer_j*jpi 197 198 197 198 199 199 answer_j=answer_j(where(answer_i EQ min(answer_i))) 200 200 answer_i=answer_i(where(answer_i EQ min(answer_i))) 201 201 202 202 answer_j=answer_j(0) 203 answer_i=answer_i(0) 204 205 pointer_Tx(ji,jj,0)=answer_i 203 answer_i=answer_i(0) 204 205 pointer_Tx(ji,jj,0)=answer_i 206 206 pointer_Ty(ji,jj,0)=answer_j 207 207 point_x1=answer_i 208 208 point_y1=answer_j 209 209 210 210 pointer_Tx(ji,jj,1)=answer_i+1 211 211 pointer_Ty(ji,jj,1)=answer_j … … 223 223 point_y4=answer_j+1 224 224 225 225 226 226 227 227 ; Now we know which 4 points are located nearby our output T point. … … 231 231 glamt_ref=replicate(0.,4) 232 232 gphit_ref=replicate(0.,4) 233 234 235 236 glamt_ref(0)=glamt_full(point_x1+1,point_y1) 233 234 235 236 glamt_ref(0)=glamt_full(point_x1+1,point_y1) 237 237 gphit_ref(0)=gphit_full(point_x1+1,point_y1) 238 238 239 glamt_ref(1)=glamt_full(point_x2+1,point_y2) 239 glamt_ref(1)=glamt_full(point_x2+1,point_y2) 240 240 gphit_ref(1)=gphit_full(point_x2+1,point_y2) 241 242 glamt_ref(2)=glamt_full(point_x3+1,point_y3) 241 242 glamt_ref(2)=glamt_full(point_x3+1,point_y3) 243 243 gphit_ref(2)=gphit_full(point_x3+1,point_y3) 244 244 245 246 glamt_ref(3)=glamt_full(point_x4+1,point_y4) 245 246 glamt_ref(3)=glamt_full(point_x4+1,point_y4) 247 247 gphit_ref(3)=gphit_full(point_x4+1,point_y4) 248 249 248 249 250 250 251 251 ; Let's make things simple for coordinates 252 252 253 253 sph_coord_ref=replicate(0.,3,5) 254 254 … … 262 262 263 263 rec_coord_ref=cv_coord(/DEGREES,FROM_SPHERE=sph_coord_ref,/TO_RECT) 264 264 265 265 xx=rec_coord_ref(0,0) 266 266 yy=rec_coord_ref(1,0) … … 270 270 y1=rec_coord_ref(1,1) 271 271 z1=rec_coord_ref(2,1) 272 272 273 273 x2=rec_coord_ref(0,2) 274 274 y2=rec_coord_ref(1,2) … … 286 286 yG=total(rec_coord_ref(1,1:4))/4. 287 287 zG=total(rec_coord_ref(2,1:4))/4. 288 288 289 289 rec_coord_ref(0,*)=rec_coord_ref(0,*)-xG 290 290 rec_coord_ref(1,*)=rec_coord_ref(1,*)-yG … … 292 292 293 293 XT=rec_coord_ref(*,1:4) 294 294 295 295 XM=MATRIX_MULTIPLY(XT,XT,/BTRANSPOSE) 296 296 EVAL=HQR(ELMHES(XM),/DOUBLE) … … 299 299 EVAL=FLOAT(EVAL) 300 300 EVEC=FLOAT(EVEC) 301 ORDER=SORT(EVAL) 302 301 ORDER=SORT(EVAL) 302 303 303 304 304 xx=total(rec_coord_ref(*,0)*EVEC(*,ORDER(1))) … … 316 316 x4=total(rec_coord_ref(*,4)*EVEC(*,ORDER(1))) 317 317 y4=total(rec_coord_ref(*,4)*EVEC(*,ORDER(2))) 318 318 319 319 320 320 ;************************************************** … … 358 358 359 359 b_coeff= -(x2-x1)*yj0-(xx-x1)*(yj0-yj1)+(yy-y1)*(xj0-xj1)+(y2-y1)*xj0 360 360 361 361 c_coeff= (xx-x1)*yj0-(yy-y1)*xj0 362 362 … … 432 432 433 433 b_coeff= -(x4-x1)*yi0-(xx-x1)*(yi0-yi1)+(yy-y1)*(xi0-xi1)+(y4-y1)*xi0 434 434 435 435 c_coeff= (xx-x1)*yi0-(yy-y1)*xi0 436 436 … … 441 441 jP=999. 442 442 443 if abs(a_coeff) GT 1.E-5 443 if abs(a_coeff) GT 1.E-5 then begin 444 444 445 445 sol1=(-b_coeff+sqrt(delta))/(2.*a_coeff) … … 461 461 462 462 endelse 463 464 463 464 465 465 weight_T(ji,jj,0) = (1.-iP)*(1.-jP) 466 466 … … 471 471 weight_T(ji,jj,3) = (1.-iP)*jP 472 472 473 473 474 474 if min(weight_T(ji,jj,*)) LT 0 OR max(weight_T(ji,jj,*)) GT 1 then begin 475 475 … … 477 477 printf, lun_err,'IP=',iP,' JP=',jP 478 478 479 479 480 480 printf, lun_err,'XP=',xx,' YP=',yy 481 481 printf, lun_err,'XA=',x1,' YA=',y1 … … 483 483 printf, lun_err,'XC=',x3,' YC=',y3 484 484 printf, lun_err,'XD=',x4,' YD=',y4 485 486 endif487 488 489 endif490 491 492 end485 486 endif 487 488 489 endif 490 491 492 end 493 493 494 494 end 495 495 496 close, lun_err497 free_lun,lun_err498 499 ; Convert to full format500 501 502 weight_full=replicate(0.,jpiglo,jpjglo,4)503 pointer_i_full=replicate(0,jpiglo,jpjglo,4)504 pointer_j_full=replicate(0,jpiglo,jpjglo,4)505 506 weight_full(1:jpiglo-2,0:jpjglo-2,*)=weight_T507 pointer_i_full(1:jpiglo-2,0:jpjglo-2,*)=pointer_Tx508 pointer_j_full(1:jpiglo-2,0:jpjglo-2,*)=pointer_Ty509 510 511 512 513 ; Let's save Weight functions and the pointer arrays514 515 516 ; Creates Netcdf File to save this517 518 519 vargrid = 'T'520 521 ; Name496 close, lun_err 497 free_lun,lun_err 498 499 ; Convert to full format 500 501 502 weight_full=replicate(0.,jpiglo,jpjglo,4) 503 pointer_i_full=replicate(0,jpiglo,jpjglo,4) 504 pointer_j_full=replicate(0,jpiglo,jpjglo,4) 505 506 weight_full(1:jpiglo-2,0:jpjglo-2,*)=weight_T 507 pointer_i_full(1:jpiglo-2,0:jpjglo-2,*)=pointer_Tx 508 pointer_j_full(1:jpiglo-2,0:jpjglo-2,*)=pointer_Ty 509 510 511 512 513 ; Let's save Weight functions and the pointer arrays 514 515 516 ; Creates Netcdf File to save this 517 518 519 vargrid = 'T' 520 521 ; Name 522 522 idout = NCDF_CREATE(weightfile,/clobber) 523 523 524 524 NCDF_CONTROL, idout, /nofill 525 ; Dimension525 ; Dimension 526 526 xidout = NCDF_DIMDEF(idout, 'x', jpiglo) 527 527 yidout = NCDF_DIMDEF(idout, 'y', jpjglo) … … 531 531 532 532 533 ; Attributes 533 ; Attributes 534 534 id_0 = NCDF_VARDEF(idout, 'nav_lon' , [xidout, yidout ], /FLOAT) 535 535 id_1 = NCDF_VARDEF(idout, 'nav_lat' , [xidout, yidout ], /FLOAT) 536 536 id_2 = NCDF_VARDEF(idout, 'deptht' , [ didout ], /FLOAT) 537 537 id_3 = NCDF_VARDEF(idout, 'time_counter', [ tidout], /FLOAT) 538 538 539 539 id0 = NCDF_VARDEF(idout, 'weight_1' , [xidout, yidout,didout,tidout ], /FLOAT) 540 540 id1 = NCDF_VARDEF(idout, 'weight_2' , [xidout, yidout,didout,tidout ], /FLOAT) … … 573 573 NCDF_ATTPUT, idout, id0, 'units','no unit' 574 574 NCDF_ATTPUT, idout, id0, 'long_name','W1' 575 575 576 576 ; Variable 1 577 577 NCDF_ATTPUT, idout, id1, 'units','no unit' … … 588 588 NCDF_ATTPUT, idout, id4, 'units','no unit' 589 589 NCDF_ATTPUT, idout, id4, 'long_name','i1' 590 590 591 591 ; Variable 1 592 592 NCDF_ATTPUT, idout, id5, 'units','no unit' … … 604 604 NCDF_ATTPUT, idout, id8, 'units','no unit' 605 605 NCDF_ATTPUT, idout, id8, 'long_name','j1' 606 606 607 607 ; Variable 1 608 608 NCDF_ATTPUT, idout, id9, 'units','no unit' … … 636 636 NCDF_VARPUT, idout, id_3, temps(0) 637 637 638 638 639 639 NCDF_VARPUT, idout, id0,weight_full(*,*,0) 640 640 NCDF_VARPUT, idout, id1,weight_full(*,*,1) … … 649 649 NCDF_VARPUT, idout, id10,pointer_j_full(*,*,2) 650 650 NCDF_VARPUT, idout, id11,pointer_j_full(*,*,3) 651 651 652 652 653 653 -
trunk/INTERP2/naminterp2.pro
r2 r48 21 21 ; 22 22 ; MODIFICATION HISTORY: 08/2002 Robinson Hordoir 23 ; 03/2003 Robinson Hordoir & Christian Ethe 23 ; 03/2003 Robinson Hordoir & Christian Ethe 24 24 ; 25 25 ;------------------------------------------------------------ … … 27 27 ;------------------------------------------------------------ 28 28 ;--------------------------------------------------- 29 ; INTERP2 : NAMELIST 29 ; INTERP2 : NAMELIST 30 30 ;--------------------------------------------------- 31 31 … … 36 36 37 37 ; From which ORCA grid 38 ; Please give geographical accuracy of your input grid 38 ; Please give geographical accuracy of your input grid 39 39 ; (ex. for ORCA2, grid_def=2.) 40 40 … … 71 71 ;**************************************** 72 72 73 ; Number of levels 73 ; Number of levels 74 74 75 75 nlevel=1 … … 95 95 do_rempl = 0 96 96 97 ; If extrapolation is required, how many i nterations shall be done ?97 ; If extrapolation is required, how many iterations shall be done ? 98 98 99 99 niterations = 10 -
trunk/INTERP2/newinterp.pro
r2 r48 12 12 ; 13 13 ; INPUTS : 14 ; ORCA_GEO grid file + ORCA_GEO weight file + Input field 14 ; ORCA_GEO grid file + ORCA_GEO weight file + Input field 15 15 ; on ORCA grid 16 16 ; OUTPUTS : … … 108 108 109 109 for jj=0,jpj-1 do begin 110 111 110 111 112 112 if mask_reg(ji,jj,0) NE 0 then begin 113 113 … … 117 117 + weight(ji,jj,2)*input_field(pointer_i(ji,jj,2)+1,pointer_j(ji,jj,2),*,*) $ 118 118 + weight(ji,jj,3)*input_field(pointer_i(ji,jj,3)+1,pointer_j(ji,jj,3),*,*) 119 119 120 120 endif 121 121 … … 131 131 ; Data of the input grid is duplicated for the northern boundary condition 132 132 133 zdata_inv(0:jpi/2 -1 ,*,*,*) = zdata_north(jpi/2:jpi-1,*,*,*) 134 zdata_inv(jpi/2:jpi-1 ,*,*,*) = zdata_north(0:jpi/2-1,*,*,*) 133 zdata_inv(0:jpi/2 -1 ,*,*,*) = zdata_north(jpi/2:jpi-1,*,*,*) 134 zdata_inv(jpi/2:jpi-1 ,*,*,*) = zdata_north(0:jpi/2-1,*,*,*) 135 135 zdata_inv = reverse(zdata_inv,2) 136 136 zdata_north(*,jpj:jpj+jplus-1,*,*) = zdata_inv(*,jpj:jpj+jplus-1,*,*) … … 178 178 179 179 ; Attributes 180 180 181 181 id0 = NCDF_VARDEF(idout, 'nav_lon' , [xidout, yidout ],/FLOAT ) 182 182 id1 = NCDF_VARDEF(idout, 'nav_lat' , [xidout, yidout ],/FLOAT ) … … 207 207 208 208 209 209 210 210 NCDF_CONTROL, idout, /ENDEF 211 211 … … 222 222 NCDF_CLOSE, idout 223 223 224 225 224 225 226 226 227 227 end 228 229 230 231 232 233 234 235 -
trunk/angle.pro
r2 r48 16 16 ;; common 17 17 ;; /comcoh/ glamf, : longitudes and latitudes at F-points 18 ;; gphif 18 ;; gphif 19 19 ;; 20 20 ;; Output : 21 21 ;; ------- 22 22 ;; common 23 ;; /comcoh/ gsinu,gcosu : sinus and cosinus of the angle 24 ;; gsinv,gcosv between north-south direction 25 ;; 23 ;; /comcoh/ gsinu,gcosu : sinus and cosinus of the angle 24 ;; gsinv,gcosv between north-south direction 25 ;; and the j-direction of the mesh 26 26 ;; 27 27 ;; Modifications: … … 37 37 ;;;--------------------------------------------------------------------- 38 38 ; 39 ;40 39 ; I. Compute the cosinus and sinus 41 40 ; ================================ 42 ; 41 ; 43 42 ; ... north pole direction & modulous (at u-point) 44 43 zxnpu = 0. - fsxspp( glamu, gphiu ) 45 44 zynpu = 0. - fsyspp( glamu, gphiu ) 46 45 znnpu = zxnpu*zxnpu + zynpu*zynpu 47 ; 46 ; 48 47 ; ... north pole direction & modulous (at v-point) 49 48 zxnpv = 0. - fsxspp( glamv, gphiv ) 50 49 zynpv = 0. - fsyspp( glamv, gphiv ) 51 50 znnpv = zxnpv*zxnpv + zynpv*zynpv 52 ; 51 ; 53 52 ; ... j-direction: f-point segment direction (u-point) 54 53 zxffu= fsxspp( glamf, gphif ) - fsxspp( shift(glamf, 0, 1), shift(gphif, 0, 1) ) … … 56 55 zmnpfu= sqrt ( znnpu * ( zxffu*zxffu + zyffu*zyffu ) ) 57 56 58 59 ; 57 58 ; 60 59 ; ... i-direction: f-point segment direction (v-point) 61 60 zxffv= fsxspp( glamf, gphif ) - fsxspp( shift(glamf, 1, 0), shift(gphif, 1, 0) ) 62 61 zyffv= fsyspp( glamf, gphif ) - fsyspp( shift(glamf, 1, 0), shift(gphif, 1, 0) ) 63 62 zmnpfv= sqrt ( znnpv * ( zxffv*zxffv + zyffv*zyffv ) ) 64 65 ; 63 64 ; 66 65 ; ... cosinus and sinus using scalar and vectorial products 67 66 gsinu = ( zxnpu*zyffu - zynpu*zxffu ) / zmnpfu … … 69 68 70 69 71 ; 70 ; 72 71 ; ... cosinus and sinus using scalar and vectorial products 73 72 ; (caution, rotation of 90 degres) … … 75 74 gcosv =-( zxnpv*zyffv - zynpv*zxffv ) / zmnpfv 76 75 77 78 ; 79 ; 76 77 ; 78 ; 80 79 ; II. Geographic mesh 81 80 ; =================== 82 ; 83 ; 81 ; 82 ; 84 83 ind = where(abs(glamf-shift(glamf, 0, 1)) LT 1.e-8) 85 84 gsinu(ind) = 0.d 86 85 gcosu(ind) = 1.d 87 ; 86 ; 88 87 ind = where(abs(gphif-shift(gphif, 1, 0)) LT 1.e-8) 89 88 gsinv(ind) = 0.d … … 91 90 ; 92 91 return 93 END 92 END -
trunk/common_interp.pro
r41 r48 1 1 ;+ 2 2 ; 3 ; @history 3 ; @history 4 4 ; 11/99 A. Jouzeau 5 5 ; 10/2001 R. Hordoir, Added variables for … … 8 8 ; 03/2003 R. Hordoir, Added variables for 9 9 ; irregular input grid 10 ; 10 ; 11 11 ;- 12 12 ;----------------------------------------------------------------------- -
trunk/condmag_from_orca.pro
r41 r48 27 27 ; to specify on which grid we do the interpolation T, U, V 28 28 ; must belong to T, U or V 29 ; ++ je ne sais pas comment ce param ètre doit intervenir29 ; ++ je ne sais pas comment ce paramÚtre doit intervenir 30 30 ; 31 31 ; @keyword DRAKKAR_EXP … … 33 33 ; only used when orcares = ORCA025 34 34 ; must be G42 ++ G70 35 ; ++ je ne sais pas comment ce param ètre doit intervenir a priori36 ; la partie mesh est la m ême pour toutes les expériences mais attention37 ; pas forc ément le mask.35 ; ++ je ne sais pas comment ce paramÚtre doit intervenir a priori 36 ; la partie mesh est la même pour toutes les expériences mais attention 37 ; pas forcément le mask. 38 38 ; 39 39 ; @keyword PERF … … 88 88 PRO condmag_from_orca, orcares, method, gridtype, DRAKKAR_EXP = drakkar_exp, PERF = perf 89 89 ; 90 compile_opt idl2, strictarrsubs 91 ; 92 IF keyword_set(perf) EQ 1 THEN BEGIN 93 msg = 'iii : start profiler' 90 compile_opt idl2, strictarrsubs 91 ; 92 IF keyword_set(perf) EQ 1 THEN BEGIN 93 msg = 'iii : start profiler' 94 ras = report(msg) 95 ;PROFILER, /SYSTEM & PROFILER 96 PROFILER, /SYSTEM & PROFILER 97 ; ++ ne tracera que les modules compilés au moment de l'appel à cette commande 98 ; donc ++ regarder si on a bien tout et si on doit ajouter des PROFILER,toto dans 99 ; tous les modules 100 ENDIF 101 ; 102 ;---- 103 ; check input parameters 104 ;---- 105 ; 106 ; check orcares definition 107 ; 108 CASE orcares OF 109 'ORCA2': BEGIN 110 msg = 'iii : valid orcares parameter = ' + orcares 111 ras = report(msg) 112 filename_oce = 'meshmask_bab.nc' 113 IF keyword_set(DRAKKAR_EXP) THEN BEGIN 114 msg = 'www : unused DRAKKAR_EXP keyword = ' + drakkar_exp 115 ras = report(msg) 116 END 117 END 118 'ORCA025': BEGIN 119 msg = 'iii : valid orcares parameter = ' + orcares 120 ras = report(msg) 121 IF keyword_set(DRAKKAR_EXP) THEN BEGIN 122 msg = 'iii : DRAKKAR_EXP keyword set' 123 ras = report(msg) 124 msg = 'iii : DRAKKAR_EXP = ' + drakkar_exp 125 ras = report(msg) 126 CASE drakkar_exp OF 127 'G42' : BEGIN 128 msg = 'iii : valid DRAKKAR_EXP keyword = ' + drakkar_exp 129 ras = report(msg) 130 END 131 'G70' : BEGIN 132 msg = 'iii : valid DRAKKAR_EXP keyword = ' + drakkar_exp 133 ras = report(msg) 134 END 135 ELSE : BEGIN 136 msg = 'eee : invalid DRAKKAR_EXP keyword = ' + drakkar_exp 137 ras = report(msg) 138 RETURN 139 END 140 ENDCASE 141 filename_oce = orcares + '-' + drakkar_exp + '_mesh_hgr.nc' 142 ENDIF ELSE BEGIN 143 msg = 'eee : unset DRAKKAR_EXP keyword' 144 ras = report(msg) 145 msg = 'eee : orcares must be G42 or G70' 146 ras = report(msg) 147 RETURN 148 ENDELSE 149 END 150 ELSE : BEGIN 151 msg = 'eee : invalid orcares parameter = ' + orcares 152 ras = report(msg) 153 msg = 'eee : orcares must be ORCA2 or ORCA025' 154 ras = report(msg) 155 RETURN 156 END 157 ENDCASE 158 ; 159 ; check method definition 160 CASE method OF 161 'bilinear': BEGIN 162 msg = 'iii : valid method parameter = ' + method 163 ras = report(msg) 164 END 165 ELSE : BEGIN 166 msg = 'eee : invalid method parameter = ' + method 167 ras = report(msg) 168 msg = 'eee : method must be bilinear' 169 ras = report(msg) 170 RETURN 171 END 172 ENDCASE 173 ; 174 ; check gridtype definition 175 gridtype = strupcase(gridtype) 176 CASE gridtype OF 177 'T' : BEGIN 178 msg = 'iii : valid gridtype parameter = ' + gridtype 179 ras = report(msg) 180 END 181 'U' : BEGIN 182 msg = 'iii : valid gridtype parameter = ' + gridtype 183 ras = report(msg) 184 END 185 'V' : BEGIN 186 msg = 'iii : valid gridtype parameter = ' + gridtype 187 ras = report(msg) 188 END 189 ELSE : BEGIN 190 msg = 'eee : invalid gridtype parameter = ' + gridtype 191 ras = report(msg) 192 msg = 'eee : gridtype must be T, U or V' 193 ras = report(msg) 194 RETURN 195 END 196 ENDCASE 197 ; 198 ; check for input files 199 ; 200 ; test if ${GEOMAG_ID} defined 201 geomag_id_env=GETENV('GEOMAG_ID') 202 CASE geomag_id_env OF 203 '' : BEGIN 204 msg = 'eee : ${GEOMAG_ID} is not defined' 205 ras = report(msg) 206 RETURN 207 END 208 ELSE: BEGIN 209 msg = 'iii : ${GEOMAG_ID} is ' + geomag_id_env 210 ras = report(msg) 211 END 212 ENDCASE 213 ; 214 iodirin = isadirectory(geomag_id_env) 215 ; 216 ; existence and protection of ${GEOMAG_ID} 217 IF (FILE_TEST(iodirin, /DIRECTORY,/EXECUTABLE , /READ) EQ 0) THEN BEGIN 218 msg = 'eee : the directory' + iodirin + ' is not accessible.' 219 ras = report(msg) 220 RETURN 221 ENDIF 222 ; 223 ; build input filenames 224 filename_cond_sed = 'cond_sed' + '_' + orcares +'.nc' 225 fullfilename_cond_sed = iodirin + filename_cond_sed 226 ; 227 ; check if this file exists 228 fullfilename_cond_sed = isafile(iodirin + filename_cond_sed, NEW=0,/MUST_EXIST) 229 IF fullfilename_cond_sed[0] EQ '' THEN BEGIN 230 msg = 'eee : the file ' + fullfilename_cond_sed + ' was not found.' 231 ras = report(msg) 232 RETURN 233 ENDIF 234 ; 235 ; protection 236 IF (FILE_TEST(fullfilename_cond_sed[0], /READ) EQ 0) THEN BEGIN 237 msg = 'eee : the file ' + fullfilename_cond_sed[0] + ' is not readable.' 238 ras = report(msg) 239 RETURN 240 ENDIF 241 ; 242 filename_Br = 'Br' + '_' + orcares +'.nc' 243 fullfilename_Br = iodirin + filename_Br 244 ; 245 ; check if this file exists 246 fullfilename_Br = isafile(iodirin + filename_Br, NEW=0,/MUST_EXIST) 247 IF fullfilename_Br[0] EQ '' THEN BEGIN 248 msg = 'eee : the file ' + fullfilename_Br + ' was not found.' 249 ras = report(msg) 250 RETURN 251 ENDIF 252 ; 253 ; protection 254 IF (FILE_TEST(fullfilename_Br[0], /READ) EQ 0) THEN BEGIN 255 msg = 'eee : the file ' + fullfilename_Br[0] + ' is not readable.' 256 ras = report(msg) 257 RETURN 258 ENDIF 259 ; 260 ; mesh mask 261 ; check if this file exists 262 fullfilename_oce = isafile(iodirin + filename_oce, NEW=0, /MUST_EXIST, $ 263 RECURSIVE=0) 264 IF fullfilename_oce[0] EQ '' THEN BEGIN 265 msg = 'eee : the file ' + fullfilename_oce + ' was not found.' 266 ras = report(msg) 267 RETURN 268 ENDIF 269 ; 270 ; protection 271 IF (FILE_TEST(fullfilename_oce[0], /READ) EQ 0) THEN BEGIN 272 msg = 'eee : the file ' + fullfilename_oce[0] + ' is not readable.' 273 ras = report(msg) 274 RETURN 275 ENDIF 276 ; 277 filename_condmag = 'condmag.nc' 278 fullfilename_condmag = iodirin + filename_condmag 279 ; 280 ; check if this file exists 281 fullfilename_condmag = isafile(iodirin + filename_condmag, NEW=0,/MUST_EXIST) 282 IF fullfilename_condmag[0] EQ '' THEN BEGIN 283 msg = 'eee : the file ' + fullfilename_condmag + ' was not found.' 284 ras = report(msg) 285 RETURN 286 ENDIF 287 ; 288 ; protection 289 IF (FILE_TEST(fullfilename_condmag[0], /READ) EQ 0) THEN BEGIN 290 msg = 'eee : the file ' + fullfilename_condmag[0] + ' is not readable.' 291 ras = report(msg) 292 RETURN 293 ENDIF 294 ; 295 ; test if ${GEOMAG_OD} defined 296 geomag_od_env=GETENV('GEOMAG_OD') 297 CASE geomag_od_env OF 298 '' : BEGIN 299 msg = 'eee : ${GEOMAG_OD} is not defined' 300 ras = report(msg) 301 RETURN 302 END 303 ELSE: BEGIN 304 msg = 'iii : ${GEOMAG_OD} is ' + geomag_od_env 305 ras = report(msg) 306 END 307 ENDCASE 308 ; 309 ; check if output data will be possible 310 iodirout = isadirectory(geomag_od_env) 311 ; 312 ; existence and protection 313 IF (FILE_TEST(iodirout, /DIRECTORY,/WRITE) EQ 0) THEN BEGIN 314 msg = 'eee : the directory' + iodirout + ' was not found.' 315 ras = report(msg) 316 RETURN 317 ENDIF 318 ; 319 ; build output filename 320 filename = 'condmag_from_' + orcares +'.nc' 321 fullfilename = iodirout + filename 322 ; 323 ; in order to avoid unexpected overwritten 324 IF (FILE_TEST(fullfilename) EQ 1) THEN BEGIN 325 msg = 'eee : the file ' + fullfilename + ' already exists.' 326 ras = report(msg) 327 RETURN 328 ENDIF 329 ; 330 ;---- 331 ; Oceanic grid parameters ie input grid + mask 332 ;---- 333 ; 334 initocemesh, orcares, DRAKKAR_EXP = drakkar_exp 335 ; 336 ;---- 337 ; read data to interpolate on regular grid 338 ;---- 339 ; 340 a=read_ncdf('cond_sed', 0, /timestep, $ 341 file = fullfilename_cond_sed[0], /nostruct) 342 msg = 'iii : ' + fullfilename_cond_sed[0] + ' opened for read' 94 343 ras = report(msg) 95 ;PROFILER, /SYSTEM & PROFILER 96 PROFILER, /SYSTEM & PROFILER 97 ; ++ ne tracera que les modules compilés au moment de l'appel à cette commande 98 ; donc ++ regarder si on a bien tout et si on doit ajouter des PROFILER,toto dans 99 ; tous les modules 100 ENDIF 101 ; 102 ;---- 103 ; check input parameters 104 ;---- 105 ; 106 ; check orcares definition 107 ; 108 CASE orcares OF 109 'ORCA2': BEGIN 110 msg = 'iii : valid orcares parameter = ' + orcares 111 ras = report(msg) 112 filename_oce = 'meshmask_bab.nc' 113 IF keyword_set(DRAKKAR_EXP) THEN BEGIN 114 msg = 'www : unused DRAKKAR_EXP keyword = ' + drakkar_exp 115 ras = report(msg) 116 END 117 END 118 'ORCA025': BEGIN 119 msg = 'iii : valid orcares parameter = ' + orcares 120 ras = report(msg) 121 IF keyword_set(DRAKKAR_EXP) THEN BEGIN 122 msg = 'iii : DRAKKAR_EXP keyword set' 123 ras = report(msg) 124 msg = 'iii : DRAKKAR_EXP = ' + drakkar_exp 125 ras = report(msg) 126 CASE drakkar_exp OF 127 'G42' : BEGIN 128 msg = 'iii : valid DRAKKAR_EXP keyword = ' + drakkar_exp 129 ras = report(msg) 130 END 131 'G70' : BEGIN 132 msg = 'iii : valid DRAKKAR_EXP keyword = ' + drakkar_exp 133 ras = report(msg) 134 END 135 ELSE : BEGIN 136 msg = 'eee : invalid DRAKKAR_EXP keyword = ' + drakkar_exp 137 ras = report(msg) 138 RETURN 139 END 140 ENDCASE 141 filename_oce = orcares + '-' + drakkar_exp + '_mesh_hgr.nc' 142 ENDIF ELSE BEGIN 143 msg = 'eee : unset DRAKKAR_EXP keyword' 144 ras = report(msg) 145 msg = 'eee : orcares must be G42 or G70' 146 ras = report(msg) 147 RETURN 148 ENDELSE 149 END 150 ELSE : BEGIN 151 msg = 'eee : invalid orcares parameter = ' + orcares 152 ras = report(msg) 153 msg = 'eee : orcares must be ORCA2 or ORCA025' 154 ras = report(msg) 155 RETURN 156 END 157 ENDCASE 158 ; 159 ; check method definition 160 CASE method OF 161 'bilinear': BEGIN 162 msg = 'iii : valid method parameter = ' + method 163 ras = report(msg) 164 END 165 ELSE : BEGIN 166 msg = 'eee : invalid method parameter = ' + method 167 ras = report(msg) 168 msg = 'eee : method must be bilinear' 169 ras = report(msg) 170 RETURN 171 END 172 ENDCASE 173 ; 174 ; check gridtype definition 175 gridtype = strupcase(gridtype) 176 CASE gridtype OF 177 'T' : BEGIN 178 msg = 'iii : valid gridtype parameter = ' + gridtype 179 ras = report(msg) 180 END 181 'U' : BEGIN 182 msg = 'iii : valid gridtype parameter = ' + gridtype 183 ras = report(msg) 184 END 185 'V' : BEGIN 186 msg = 'iii : valid gridtype parameter = ' + gridtype 187 ras = report(msg) 188 END 189 ELSE: BEGIN 190 msg = 'eee : invalid gridtype parameter = ' + gridtype 191 ras = report(msg) 192 msg = 'eee : gridtype must be T, U or V' 193 ras = report(msg) 194 RETURN 195 END 196 ENDCASE 197 ; 198 ; check for input files 199 ; 200 ; test if ${GEOMAG_ID} defined 201 geomag_id_env=GETENV('GEOMAG_ID') 202 CASE geomag_id_env OF 203 '' : BEGIN 204 msg = 'eee : ${GEOMAG_ID} is not defined' 205 ras = report(msg) 206 RETURN 207 END 208 ELSE: BEGIN 209 msg = 'iii : ${GEOMAG_ID} is ' + geomag_id_env 210 ras = report(msg) 211 END 212 ENDCASE 213 ; 214 iodirin = isadirectory(geomag_id_env) 215 ; 216 ; existence and protection of ${GEOMAG_ID} 217 IF (FILE_TEST(iodirin, /DIRECTORY,/EXECUTABLE , /READ) EQ 0) THEN BEGIN 218 msg = 'eee : the directory' + iodirin + ' is not accessible.' 219 ras = report(msg) 220 RETURN 221 ENDIF 222 ; 223 ; build input filenames 224 filename_cond_sed = 'cond_sed' + '_' + orcares +'.nc' 225 fullfilename_cond_sed = iodirin + filename_cond_sed 226 ; 227 ; check if this file exists 228 fullfilename_cond_sed = isafile(iodirin + filename_cond_sed, NEW=0,/MUST_EXIST) 229 IF fullfilename_cond_sed[0] EQ '' THEN BEGIN 230 msg = 'eee : the file ' + fullfilename_cond_sed + ' was not found.' 231 ras = report(msg) 232 RETURN 233 ENDIF 234 ; 235 ; protection 236 IF (FILE_TEST(fullfilename_cond_sed[0], /READ) EQ 0) THEN BEGIN 237 msg = 'eee : the file ' + fullfilename_cond_sed[0] + ' is not readable.' 238 ras = report(msg) 239 RETURN 240 ENDIF 241 ; 242 filename_Br = 'Br' + '_' + orcares +'.nc' 243 fullfilename_Br = iodirin + filename_Br 244 ; 245 ; check if this file exists 246 fullfilename_Br = isafile(iodirin + filename_Br, NEW=0,/MUST_EXIST) 247 IF fullfilename_Br[0] EQ '' THEN BEGIN 248 msg = 'eee : the file ' + fullfilename_Br + ' was not found.' 249 ras = report(msg) 250 RETURN 251 ENDIF 252 ; 253 ; protection 254 IF (FILE_TEST(fullfilename_Br[0], /READ) EQ 0) THEN BEGIN 255 msg = 'eee : the file ' + fullfilename_Br[0] + ' is not readable.' 256 ras = report(msg) 257 RETURN 258 ENDIF 259 ; 260 ; mesh mask 261 ; check if this file exists 262 fullfilename_oce = isafile(iodirin + filename_oce, NEW=0, /MUST_EXIST, $ 263 RECURSIVE=0) 264 IF fullfilename_oce[0] EQ '' THEN BEGIN 265 msg = 'eee : the file ' + fullfilename_oce + ' was not found.' 266 ras = report(msg) 267 RETURN 268 ENDIF 269 ; 270 ; protection 271 IF (FILE_TEST(fullfilename_oce[0], /READ) EQ 0) THEN BEGIN 272 msg = 'eee : the file ' + fullfilename_oce[0] + ' is not readable.' 273 ras = report(msg) 274 RETURN 275 ENDIF 276 ; 277 filename_condmag = 'condmag.nc' 278 fullfilename_condmag = iodirin + filename_condmag 279 ; 280 ; check if this file exists 281 fullfilename_condmag = isafile(iodirin + filename_condmag, NEW=0,/MUST_EXIST) 282 IF fullfilename_condmag[0] EQ '' THEN BEGIN 283 msg = 'eee : the file ' + fullfilename_condmag + ' was not found.' 284 ras = report(msg) 285 RETURN 286 ENDIF 287 ; 288 ; protection 289 IF (FILE_TEST(fullfilename_condmag[0], /READ) EQ 0) THEN BEGIN 290 msg = 'eee : the file ' + fullfilename_condmag[0] + ' is not readable.' 291 ras = report(msg) 292 RETURN 293 ENDIF 294 ; 295 ; test if ${GEOMAG_OD} defined 296 geomag_od_env=GETENV('GEOMAG_OD') 297 CASE geomag_od_env OF 298 '' : BEGIN 299 msg = 'eee : ${GEOMAG_OD} is not defined' 300 ras = report(msg) 301 RETURN 302 END 303 ELSE: BEGIN 304 msg = 'iii : ${GEOMAG_OD} is ' + geomag_od_env 305 ras = report(msg) 306 END 307 ENDCASE 308 ; 309 ; check if output data will be possible 310 iodirout = isadirectory(geomag_od_env) 311 ; 312 ; existence and protection 313 IF (FILE_TEST(iodirout, /DIRECTORY,/WRITE) EQ 0) THEN BEGIN 314 msg = 'eee : the directory' + iodirout + ' was not found.' 315 ras = report(msg) 316 RETURN 317 ENDIF 318 ; 319 ; build output filename 320 filename = 'condmag_from_' + orcares +'.nc' 321 fullfilename = iodirout + filename 322 ; 323 ; in order to avoid unexpected overwritten 324 IF (FILE_TEST(fullfilename) EQ 1) THEN BEGIN 325 msg = 'eee : the file ' + fullfilename + ' already exists.' 326 ras = report(msg) 327 RETURN 328 ENDIF 329 ; 330 ;---- 331 ; Oceanic grid parameters ie input grid + mask 332 ;---- 333 ; 334 initocemesh, orcares, DRAKKAR_EXP = drakkar_exp 335 ; 336 ;---- 337 ; read data to interpolate on regular grid 338 ;---- 339 ; 340 a=read_ncdf('cond_sed', 0, /timestep, $ 341 file = fullfilename_cond_sed[0], /nostruct) 342 msg = 'iii : ' + fullfilename_cond_sed[0] + ' opened for read' 343 ras = report(msg) 344 ; mask and limite cond_sed on values below 5000 ++ 345 cond_sedin = extrapsmooth(a, tmask[*, *, 0]*(a le 5000), /x_periodic) 346 ; 347 b = read_ncdf('Br', 0, /timestep, $ 348 file = fullfilename_Br[0], /nostruct) 349 msg = 'iii : ' + fullfilename_Br[0] + ' opened for read' 350 ras = report(msg) 351 ; mask and limit Br on values below/above ++ 352 Brin = extrapsmooth(b, tmask[*, *, 0]*(b), /x_periodic) 353 ; 354 glamin = glamt 355 gphiin = gphit 356 maskin = tmask[*, *, 0] 357 ; 358 ;--- 359 ; conductivity and magnetic field grid parameters ie output grid 360 ;---- 344 ; mask and limite cond_sed on values below 5000 ++ 345 cond_sedin = extrapsmooth(a, tmask[*, *, 0]*(a le 5000), /x_periodic) 346 ; 347 b = read_ncdf('Br', 0, /timestep, $ 348 file = fullfilename_Br[0], /nostruct) 349 msg = 'iii : ' + fullfilename_Br[0] + ' opened for read' 350 ras = report(msg) 351 ; mask and limit Br on values below/above ++ 352 Brin = extrapsmooth(b, tmask[*, *, 0]*(b), /x_periodic) 353 ; 354 glamin = glamt 355 gphiin = gphit 356 maskin = tmask[*, *, 0] 357 ; 358 ;--- 359 ; conductivity and magnetic field grid parameters ie output grid 360 ;---- 361 361 362 initncdf, fullfilename_condmag[0], xaxisname = 'lo', yaxisname = 'la'363 msg = 'iii : ' + fullfilename_condmag[0] + ' opened for read'364 ras = report(msg)365 ;366 glamout = glamt367 gphiout = gphit368 ;362 initncdf, fullfilename_condmag[0], xaxisname = 'lo', yaxisname = 'la' 363 msg = 'iii : ' + fullfilename_condmag[0] + ' opened for read' 364 ras = report(msg) 365 ; 366 glamout = glamt 367 gphiout = gphit 368 ; 369 369 370 dimidxout = jpi371 dimidyout = jpj372 ;373 ;---------------374 ; Interpolate375 ;---------------376 ;++ cond_sedout = fromirr('bilinear', cond_sedin, glamin, gphiin, maskin, glamout, gphiout, -1, WEIG = weights, ADDR = addresses)370 dimidxout = jpi 371 dimidyout = jpj 372 ; 373 ;--------------- 374 ; Interpolate 375 ;--------------- 376 ;++ cond_sedout = fromirr('bilinear', cond_sedin, glamin, gphiin, maskin, glamout, gphiout, -1, WEIG = weights, ADDR = addresses) 377 377 cond_sedout=dindgen(dimidxout,dimidyout) 378 ;Brout=fromirr('bilinear', Brin, weights, addresses) ; ++ pb doc % Variable is undefined: LONOUT379 ;++Brout = fromirr('bilinear', Brin, glamin, gphiin, maskin, glamout, gphiout, -1, WEIG = weights, ADDR = addresses)380 ;378 ;Brout=fromirr('bilinear', Brin, weights, addresses) ; ++ pb doc % Variable is undefined: LONOUT 379 ;++Brout = fromirr('bilinear', Brin, glamin, gphiin, maskin, glamout, gphiout, -1, WEIG = weights, ADDR = addresses) 380 ; 381 381 Brout=dindgen(dimidxout,dimidyout) 382 382 383 ;384 ;---------------385 ; Produce outputs386 ;---------------387 ;388 netcdf_id = NCDF_CREATE(fullfilename, /clobber)389 NCDF_CONTROL, netcdf_id, /NOFILL390 ;391 ; dimension392 dimidxout = NCDF_DIMDEF(netcdf_id, 'la' , jpi)393 dimidyout = NCDF_DIMDEF(netcdf_id, 'lo' , jpj)394 ;395 ; global attributes396 NCDF_ATTPUT, netcdf_id, 'Conventions', 'GDT 1.2', /GLOBAL ; ++ conformite397 ;++ NCDF_ATTPUT, netcdf_id, 'file_name' , filename, /GLOBAL398 ;++ NCDF_ATTPUT, netcdf_id, 'Title' , title, /GLOBAL399 ;400 ; declaration of variables401 varid = lonarr(4)402 ;403 varid[0] = NCDF_VARDEF(netcdf_id, 'la' , [dimidxout], /FLOAT)404 varid[1] = NCDF_VARDEF(netcdf_id, 'lo' , [dimidyout], /FLOAT)405 varid[2] = NCDF_VARDEF(netcdf_id, 'cond_sed' , [dimidxout, dimidyout], /FLOAT)406 NCDF_ATTPUT, netcdf_id, varid[2], 'units' , 'siemens'407 varid[3] = NCDF_VARDEF(netcdf_id, 'Br' , [dimidxout, dimidyout], /FLOAT)408 NCDF_ATTPUT, netcdf_id, varid[3], 'units' , 'siemens'409 ;410 NCDF_CONTROL, netcdf_id, /ENDEF411 ;412 NCDF_VARPUT, netcdf_id, varid[0], dindgen(dimidxout,dimidyout)413 NCDF_VARPUT, netcdf_id, varid[1], dindgen(dimidyout,dimidyout)414 ;415 ; write the data416 ;417 NCDF_VARPUT, netcdf_id, varid[2], cond_sedout418 ;++ NCDF_VARPUT, netcdf_id, 'Br', Brout419 ;---------------------------420 ; close the netcdf files421 NCDF_CLOSE, netcdf_id422 ;423 msg = 'iii : ' + fullfilename + ' created'424 ras = report(msg)425 ;426 IF keyword_set(perf) EQ 1 THEN BEGIN427 msg = 'iii : report profiler results'428 ras = report(msg)429 ; ++ tri par ordre alpha , par ordre de pourcentage croissant ou décroissant430 ; +++ d'utilisation431 profiler,/REPORT432 ; shut down all profiling (according to433 ; http://www.dfanning.com/code_tips/whyslow.html)434 profiler, /CLEAR, /SYSTEM435 profiler, /CLEAR, /RESET383 ; 384 ;--------------- 385 ; Produce outputs 386 ;--------------- 387 ; 388 netcdf_id = NCDF_CREATE(fullfilename, /clobber) 389 NCDF_CONTROL, netcdf_id, /NOFILL 390 ; 391 ; dimension 392 dimidxout = NCDF_DIMDEF(netcdf_id, 'la' , jpi) 393 dimidyout = NCDF_DIMDEF(netcdf_id, 'lo' , jpj) 394 ; 395 ; global attributes 396 NCDF_ATTPUT, netcdf_id, 'Conventions', 'GDT 1.2', /GLOBAL ; ++ conformite 397 ;++ NCDF_ATTPUT, netcdf_id, 'file_name' , filename, /GLOBAL 398 ;++ NCDF_ATTPUT, netcdf_id, 'Title' , title, /GLOBAL 399 ; 400 ; declaration of variables 401 varid = lonarr(4) 402 ; 403 varid[0] = NCDF_VARDEF(netcdf_id, 'la' , [dimidxout], /FLOAT) 404 varid[1] = NCDF_VARDEF(netcdf_id, 'lo' , [dimidyout], /FLOAT) 405 varid[2] = NCDF_VARDEF(netcdf_id, 'cond_sed' , [dimidxout, dimidyout], /FLOAT) 406 NCDF_ATTPUT, netcdf_id, varid[2], 'units' , 'siemens' 407 varid[3] = NCDF_VARDEF(netcdf_id, 'Br' , [dimidxout, dimidyout], /FLOAT) 408 NCDF_ATTPUT, netcdf_id, varid[3], 'units' , 'siemens' 409 ; 410 NCDF_CONTROL, netcdf_id, /ENDEF 411 ; 412 NCDF_VARPUT, netcdf_id, varid[0], dindgen(dimidxout,dimidyout) 413 NCDF_VARPUT, netcdf_id, varid[1], dindgen(dimidyout,dimidyout) 414 ; 415 ; write the data 416 ; 417 NCDF_VARPUT, netcdf_id, varid[2], cond_sedout 418 ;++ NCDF_VARPUT, netcdf_id, 'Br', Brout 419 ;--------------------------- 420 ; close the netcdf files 421 NCDF_CLOSE, netcdf_id 422 ; 423 msg = 'iii : ' + fullfilename + ' created' 424 ras = report(msg) 425 ; 426 IF keyword_set(perf) EQ 1 THEN BEGIN 427 msg = 'iii : report profiler results' 428 ras = report(msg) 429 ; ++ tri par ordre alpha , par ordre de pourcentage croissant ou décroissant 430 ; +++ d'utilisation 431 profiler,/REPORT 432 ; shut down all profiling (according to 433 ; http://www.dfanning.com/code_tips/whyslow.html) 434 profiler, /CLEAR, /SYSTEM 435 profiler, /CLEAR, /RESET 436 436 ENDIF 437 ;437 ; 438 438 END -
trunk/condmag_on_orca.pro
r41 r48 92 92 ;- 93 93 PRO condmag_on_orca, orcares, method, gridtype, DRAKKAR_EXP = drakkar_exp 94 ;95 compile_opt idl2, strictarrsubs96 ;97 ;----98 ; check input parameters99 ;----100 ;101 ; check orcares definition102 ;103 CASE orcares OF104 'ORCA2': BEGIN94 ; 95 compile_opt idl2, strictarrsubs 96 ; 97 ;---- 98 ; check input parameters 99 ;---- 100 ; 101 ; check orcares definition 102 ; 103 CASE orcares OF 104 'ORCA2': BEGIN 105 105 msg = 'iii : valid orcares parameter = ' + orcares 106 106 ras = report(msg) … … 111 111 END 112 112 END 113 'ORCA025': BEGIN114 msg = 'iii : valid orcares parameter = ' + orcares113 'ORCA025': BEGIN 114 msg = 'iii : valid orcares parameter = ' + orcares 115 115 ras = report(msg) 116 116 IF keyword_set(DRAKKAR_EXP) THEN BEGIN … … 136 136 filename_oce = orcares + '-' + drakkar_exp + '_mesh_hgr.nc' 137 137 ENDIF ELSE BEGIN 138 msg = 'eee : unset DRAKKAR_EXP keyword' 138 msg = 'eee : unset DRAKKAR_EXP keyword' 139 139 ras = report(msg) 140 140 msg = 'eee : orcares must be G42 or G70' … … 143 143 ENDELSE 144 144 END 145 ELSE: BEGIN145 ELSE : BEGIN 146 146 msg = 'eee : invalid orcares parameter = ' + orcares 147 147 ras = report(msg) … … 149 149 ras = report(msg) 150 150 RETURN 151 152 ENDCASE153 ;154 ; check method definition155 CASE method OF156 'bilinear': BEGIN151 END 152 ENDCASE 153 ; 154 ; check method definition 155 CASE method OF 156 'bilinear': BEGIN 157 157 msg = 'iii : valid method parameter = ' + method 158 158 ras = report(msg) … … 162 162 ras = report(msg) 163 163 END 164 ELSE : BEGIN164 ELSE : BEGIN 165 165 msg = 'eee : invalid method parameter = ' + method 166 166 ras = report(msg) … … 169 169 RETURN 170 170 END 171 ENDCASE172 ;173 ; check gridtype definition174 gridtype = strupcase(gridtype)175 CASE gridtype OF176 'T' : BEGIN171 ENDCASE 172 ; 173 ; check gridtype definition 174 gridtype = strupcase(gridtype) 175 CASE gridtype OF 176 'T' : BEGIN 177 177 msg = 'iii : valid gridtype parameter = ' + gridtype 178 178 ras = report(msg) … … 186 186 ras = report(msg) 187 187 END 188 ELSE : BEGIN188 ELSE : BEGIN 189 189 msg = 'eee : invalid gridtype parameter = ' + gridtype 190 190 ras = report(msg) … … 193 193 RETURN 194 194 END 195 ENDCASE196 ;197 ; check for input files198 ;199 ; test if ${GEOMAG_ID} defined200 geomag_id_env=GETENV('GEOMAG_ID')201 CASE geomag_id_env OF202 '' : BEGIN195 ENDCASE 196 ; 197 ; check for input files 198 ; 199 ; test if ${GEOMAG_ID} defined 200 geomag_id_env=GETENV('GEOMAG_ID') 201 CASE geomag_id_env OF 202 '' : BEGIN 203 203 msg = 'eee : ${GEOMAG_ID} is not defined' 204 204 ras = report(msg) … … 209 209 ras = report(msg) 210 210 END 211 ENDCASE212 ;213 iodirin = isadirectory(geomag_id_env)214 ;215 ; existence and protection of ${GEOMAG_ID}216 IF (FILE_TEST(iodirin, /DIRECTORY, /EXECUTABLE, /READ) EQ 0) THEN BEGIN217 msg = 'eee : the directory' + iodirin + ' is not accessible.'218 ras = report(msg)219 RETURN220 ENDIF221 ;222 ; conductivity and magnetic field223 filename_condmag = 'condmag.nc'224 ;225 ; check if this file exists226 fullfilename_condmag = isafile(iodirin + filename_condmag, NEW=0, /MUST_EXIST)227 IF fullfilename_condmag[0] EQ '' THEN BEGIN228 msg = 'eee : the file ' + fullfilename_condmag + ' was not found.'229 ras = report(msg)230 RETURN231 ENDIF232 ;233 ; protection234 IF (FILE_TEST(fullfilename_condmag[0], /READ) EQ 0) THEN BEGIN235 msg = 'eee : the file ' + fullfilename_condmag[0] + ' is not readable.'236 ras = report(msg)237 RETURN238 ENDIF239 ;240 ; mesh mask241 ; check if this file exists242 fullfilename_oce = isafile(iodirin + filename_oce, NEW=0, /MUST_EXIST, $243 244 IF fullfilename_oce[0] EQ '' THEN BEGIN245 msg = 'eee : the file ' + fullfilename_oce + ' was not found.'246 ras = report(msg)247 RETURN248 ENDIF249 ;250 ; protection251 IF (FILE_TEST(fullfilename_oce[0], /READ) EQ 0) THEN BEGIN252 msg = 'eee : the file ' + fullfilename_oce[0] + ' is not readable.'253 ras = report(msg)254 RETURN255 ENDIF256 ;257 ; test if ${GEOMAG_OD} defined258 geomag_od_env=GETENV('GEOMAG_OD')259 CASE geomag_od_env OF260 '' : BEGIN261 msg = 'eee : ${GEOMAG_OD} is not defined'262 ras = report(msg)211 ENDCASE 212 ; 213 iodirin = isadirectory(geomag_id_env) 214 ; 215 ; existence and protection of ${GEOMAG_ID} 216 IF (FILE_TEST(iodirin, /DIRECTORY, /EXECUTABLE, /READ) EQ 0) THEN BEGIN 217 msg = 'eee : the directory' + iodirin + ' is not accessible.' 218 ras = report(msg) 219 RETURN 220 ENDIF 221 ; 222 ; conductivity and magnetic field 223 filename_condmag = 'condmag.nc' 224 ; 225 ; check if this file exists 226 fullfilename_condmag = isafile(iodirin + filename_condmag, NEW=0, /MUST_EXIST) 227 IF fullfilename_condmag[0] EQ '' THEN BEGIN 228 msg = 'eee : the file ' + fullfilename_condmag + ' was not found.' 229 ras = report(msg) 230 RETURN 231 ENDIF 232 ; 233 ; protection 234 IF (FILE_TEST(fullfilename_condmag[0], /READ) EQ 0) THEN BEGIN 235 msg = 'eee : the file ' + fullfilename_condmag[0] + ' is not readable.' 236 ras = report(msg) 237 RETURN 238 ENDIF 239 ; 240 ; mesh mask 241 ; check if this file exists 242 fullfilename_oce = isafile(iodirin + filename_oce, NEW=0, /MUST_EXIST, $ 243 RECURSIVE=0) 244 IF fullfilename_oce[0] EQ '' THEN BEGIN 245 msg = 'eee : the file ' + fullfilename_oce + ' was not found.' 246 ras = report(msg) 247 RETURN 248 ENDIF 249 ; 250 ; protection 251 IF (FILE_TEST(fullfilename_oce[0], /READ) EQ 0) THEN BEGIN 252 msg = 'eee : the file ' + fullfilename_oce[0] + ' is not readable.' 253 ras = report(msg) 254 RETURN 255 ENDIF 256 ; 257 ; test if ${GEOMAG_OD} defined 258 geomag_od_env=GETENV('GEOMAG_OD') 259 CASE geomag_od_env OF 260 '' : BEGIN 261 msg = 'eee : ${GEOMAG_OD} is not defined' 262 ras = report(msg) 263 263 RETURN 264 264 END 265 ELSE: BEGIN266 msg = 'iii : ${GEOMAG_OD} is ' + geomag_od_env265 ELSE: BEGIN 266 msg = 'iii : ${GEOMAG_OD} is ' + geomag_od_env 267 267 ras = report(msg) 268 268 END 269 ENDCASE 270 ; 271 ; check if output data will be possible 272 iodirout = isadirectory(geomag_od_env) 273 ; 274 ; existence and protection 275 IF (FILE_TEST(iodirout, /DIRECTORY, /WRITE) EQ 0) THEN BEGIN 276 msg = 'eee : the directory' + iodirout + ' was not found.' 277 ras = report(msg) 278 RETURN 279 ENDIF 280 ; 281 ; d'après ncdump -h /usr/work/sur/fvi/OPA/geomag/condmag.nc 282 condmaglonname = 'lo' 283 condmaglatname = 'la' 284 ; 285 ;---- 286 ; conductivity and magnetic field grid parameters 287 ;---- 288 ; 289 get_gridparams, fullfilename_condmag[0], $ 290 condmaglonname, condmaglatname, $ 291 condmaglon, condmaglat, jpia, jpja, 1, /DOUBLE 292 ; 293 ;---- 294 ; Oceanic grid parameters 295 ;---- 296 ; 297 olonname = 'glam' + STRLOWCASE(gridtype) 298 olatname = 'gphi' + STRLOWCASE(gridtype) 299 get_gridparams, fullfilename_oce[0], $ 300 olonname, olatname, $ 301 olon, olat, jpio, jpjo, 2, /DOUBLE 302 msg = 'iii : ' + fullfilename_oce[0] + ' opened for read' 303 ras = report(msg) 304 ; 305 ;--------------- 306 ; Compute weight and address 307 ;--------------- 308 ; 309 CASE method OF 310 'bilinear': compute_fromreg_bilinear_weigaddr, $ 269 ENDCASE 270 ; 271 ; check if output data will be possible 272 iodirout = isadirectory(geomag_od_env) 273 ; 274 ; existence and protection 275 IF (FILE_TEST(iodirout, /DIRECTORY, /WRITE) EQ 0) THEN BEGIN 276 msg = 'eee : the directory' + iodirout + ' was not found.' 277 ras = report(msg) 278 RETURN 279 ENDIF 280 ; 281 ; d'aprÚs ncdump -h /usr/work/sur/fvi/OPA/geomag/condmag.nc 282 condmaglonname = 'lo' 283 condmaglatname = 'la' 284 ; 285 ;---- 286 ; conductivity and magnetic field grid parameters 287 ;---- 288 ; 289 get_gridparams, fullfilename_condmag[0], $ 290 condmaglonname, condmaglatname, $ 291 condmaglon, condmaglat, jpia, jpja, 1, /DOUBLE 292 ; 293 ;---- 294 ; Oceanic grid parameters 295 ;---- 296 ; 297 olonname = 'glam' + STRLOWCASE(gridtype) 298 olatname = 'gphi' + STRLOWCASE(gridtype) 299 get_gridparams, fullfilename_oce[0], $ 300 olonname, olatname, $ 301 olon, olat, jpio, jpjo, 2, /DOUBLE 302 msg = 'iii : ' + fullfilename_oce[0] + ' opened for read' 303 ras = report(msg) 304 ; 305 ;--------------- 306 ; Compute weight and address 307 ;--------------- 308 ; 309 CASE method OF 310 'bilinear': compute_fromreg_bilinear_weigaddr, $ 311 condmaglon, condmaglat, olon, olat, weig, addr 312 'imoms3' : compute_fromreg_imoms3_weigaddr, $ 311 313 condmaglon, condmaglat, olon, olat, weig, addr 312 'imoms3' : compute_fromreg_imoms3_weigaddr, $ 313 condmaglon, condmaglat, olon, olat, weig, addr 314 ENDCASE 315 ; 316 ; reading condmag file 317 netcdf_id_condmag = NCDF_OPEN(fullfilename_condmag[0], /NOWRITE) 318 msg = 'iii : ' + fullfilename_condmag[0] + ' opened for read' 319 ras = report(msg) 320 ; 321 varname_cond_sed = 'cond_sed' 322 varname_br = 'Br' 323 ; 324 varinq_cond_sed = NCDF_VARINQ(netcdf_id_condmag, varname_cond_sed) 325 NCDF_VARGET, netcdf_id_condmag, varname_cond_sed, varin_cond_sed 326 ; 327 varinq_br = NCDF_VARINQ(netcdf_id_condmag, varname_br) 328 NCDF_VARGET, netcdf_id_condmag, varname_br, varin_br 329 ; 330 NCDF_CLOSE, netcdf_id_condmag 331 ; 332 ;--------------------------- 333 ; do the interpolation 334 varin_cond_sed = TOTAL(weig*varin_cond_sed[addr], 1) 335 varin_cond_sed = REFORM(varin_cond_sed, jpio, jpjo, /OVER) 336 varin_br = TOTAL(weig*varin_br[addr], 1) 337 varin_br = REFORM(varin_br, jpio, jpjo, /OVER) 338 ; 339 varout_cond_sed = TEMPORARY(varin_cond_sed) 340 varout_br = TEMPORARY(varin_br) 341 ; 342 ; put back the masked value 343 ;++ IF bad[0] NE -1 THEN varout_cond_sed[TEMPORARY(bad)] = 32767 344 ;++ IF bad[0] NE -1 THEN varout_br[TEMPORARY(bad)] = 32767 345 ; 346 ; 347 ; produce outputs 348 condmag_output, orcares, $ 349 varinq_cond_sed, 'Conductance', 'Conductance', 'siemens', $ 350 jpio, jpjo, olon, olat, $ 351 varout_cond_sed 352 condmag_output, orcares, $ 353 varinq_br, 'Magnetic field', 'magnetic field', 'tesla', $ 354 jpio, jpjo, olon, olat, $ 355 varout_br 356 ; 314 ENDCASE 315 ; 316 ; reading condmag file 317 netcdf_id_condmag = NCDF_OPEN(fullfilename_condmag[0], /NOWRITE) 318 msg = 'iii : ' + fullfilename_condmag[0] + ' opened for read' 319 ras = report(msg) 320 ; 321 varname_cond_sed = 'cond_sed' 322 varname_br = 'Br' 323 ; 324 varinq_cond_sed = NCDF_VARINQ(netcdf_id_condmag, varname_cond_sed) 325 NCDF_VARGET, netcdf_id_condmag, varname_cond_sed, varin_cond_sed 326 ; 327 varinq_br = NCDF_VARINQ(netcdf_id_condmag, varname_br) 328 NCDF_VARGET, netcdf_id_condmag, varname_br, varin_br 329 ; 330 NCDF_CLOSE, netcdf_id_condmag 331 ; 332 ;--------------------------- 333 ; do the interpolation 334 varin_cond_sed = TOTAL(weig*varin_cond_sed[addr], 1) 335 varin_cond_sed = REFORM(varin_cond_sed, jpio, jpjo, /OVER) 336 varin_br = TOTAL(weig*varin_br[addr], 1) 337 varin_br = REFORM(varin_br, jpio, jpjo, /OVER) 338 ; 339 varout_cond_sed = TEMPORARY(varin_cond_sed) 340 varout_br = TEMPORARY(varin_br) 341 ; 342 ; put back the masked value 343 ;++ IF bad[0] NE -1 THEN varout_cond_sed[TEMPORARY(bad)] = 32767 344 ;++ IF bad[0] NE -1 THEN varout_br[TEMPORARY(bad)] = 32767 345 ; 346 ; 347 ; produce outputs 348 condmag_output, orcares, $ 349 varinq_cond_sed, 'Conductance', 'Conductance', 'siemens', $ 350 jpio, jpjo, olon, olat, $ 351 varout_cond_sed 352 condmag_output, orcares, $ 353 varinq_br, 'Magnetic field', 'magnetic field', 'tesla', $ 354 jpio, jpjo, olon, olat, $ 355 varout_br 356 ; 357 357 END -
trunk/condmag_output.pro
r36 r48 79 79 ;- 80 80 PRO condmag_output, orcares, variable, title, long_name, units, jpio, jpjo, olon, olat, values 81 ;82 compile_opt idl2, strictarrsubs83 ;84 ;----85 ; check input parameters86 ;----87 ;88 ; check orcares definition89 ;90 CASE orcares OF91 'ORCA2': BEGIN81 ; 82 compile_opt idl2, strictarrsubs 83 ; 84 ;---- 85 ; check input parameters 86 ;---- 87 ; 88 ; check orcares definition 89 ; 90 CASE orcares OF 91 'ORCA2': BEGIN 92 92 msg = 'iii : valid orcares parameter = ' + orcares 93 93 ras = report(msg) 94 95 'ORCA025': BEGIN94 END 95 'ORCA025': BEGIN 96 96 msg = 'iii : valid orcares parameter = ' + orcares 97 97 ras = report(msg) 98 99 ELSE: BEGIN98 END 99 ELSE : BEGIN 100 100 msg = 'eee : invalid orcares parameter = ' + orcares 101 101 ras = report(msg) … … 103 103 ras = report(msg) 104 104 RETURN 105 106 ENDCASE107 ;108 ; check variable definition109 ; ++ dimension etc.110 CASE variable.NAME OF111 'cond_sed': BEGIN112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 105 END 106 ENDCASE 107 ; 108 ; check variable definition 109 ; ++ dimension etc. 110 CASE variable.NAME OF 111 'cond_sed': BEGIN 112 msg = 'iii : valid variable.NAME parameter = ' + variable.NAME 113 ras = report(msg) 114 CASE title OF 115 'Conductance': BEGIN 116 msg = 'iii : valid title parameter = ' + title 117 ras = report(msg) 118 END 119 ELSE : BEGIN 120 msg = 'eee : invalid title parameter = ' + title 121 ras = report(msg) 122 msg = 'eee : title must be Conductance' 123 ras = report(msg) 124 RETURN 125 END 126 ENDCASE 127 CASE long_name OF 128 'Conductance': BEGIN 129 msg = 'iii : valid long_name parameter = ' + long_name 130 ras = report(msg) 131 END 132 ELSE : BEGIN 133 msg = 'eee : invalid long_name parameter = ' + long_name 134 ras = report(msg) 135 msg = 'eee : long_name must be Conductance' 136 ras = report(msg) 137 RETURN 138 END 139 ENDCASE 140 CASE units OF 141 'siemens': BEGIN 142 142 msg = 'iii : valid units parameter = ' + units 143 143 ras = report(msg) 144 145 144 END 145 ELSE : BEGIN 146 146 msg = 'eee : invalid units parameter = ' + units 147 147 ras = report(msg) … … 149 149 ras = report(msg) 150 150 RETURN 151 END 152 ENDCASE 153 END 154 'Br' : BEGIN 155 msg = 'iii : valid variable.NAME parameter = ' + variable.NAME 151 END 152 ENDCASE 153 END 154 'Br' : BEGIN 155 msg = 'iii : valid variable.NAME parameter = ' + variable.NAME 156 ras = report(msg) 157 CASE title OF 158 'Magnetic field': BEGIN 159 msg = 'iii : valid title parameter = ' + title 160 ras = report(msg) 161 END 162 ELSE : BEGIN 163 msg = 'eee : invalid title parameter = ' + title 156 164 ras = report(msg) 157 CASE title OF 158 'Magnetic field': BEGIN 159 msg = 'iii : valid title parameter = ' + title 160 ras = report(msg) 161 END 162 ELSE : BEGIN 163 msg = 'eee : invalid title parameter = ' + title 164 ras = report(msg) 165 msg = 'eee : title must be Magnetic field' 166 ras = report(msg) 167 RETURN 168 END 169 ENDCASE 170 CASE long_name OF 171 'magnetic field': BEGIN 172 msg = 'iii : valid long_name parameter = ' + long_name 173 ras = report(msg) 174 END 175 ELSE : BEGIN 176 msg = 'eee : invalid long_name parameter = ' + long_name 177 ras = report(msg) 178 msg = 'eee : long_name must be magnetic field' 179 ras = report(msg) 180 RETURN 181 END 182 ENDCASE 183 CASE units OF 184 'tesla': BEGIN 185 msg = 'iii : valid units parameter = ' + units 186 ras = report(msg) 187 END 188 ELSE : BEGIN 189 msg = 'eee : invalid units parameter = ' + units 190 ras = report(msg) 191 msg = 'eee : units must be tesla' 192 ras = report(msg) 193 RETURN 194 END 195 ENDCASE 196 END 197 ELSE : BEGIN 198 msg = 'eee : invalid variable.NAME parameter = ' + variable.NAME 199 ras = report(msg) 200 msg = 'eee : variable.NAME must be cond_sed or Br' 165 msg = 'eee : title must be Magnetic field' 201 166 ras = report(msg) 202 167 RETURN 203 168 END 204 ENDCASE 205 ; 206 ; check jpio 207 IF (jpio LT 1) THEN BEGIN 208 msg = ' eee : invalid jpio parameter = ' + STRING(jpio) 209 ras = report(msg) 210 msg = ' eee : jpio must be greater or equal to 1' 211 ras = report(msg) 212 RETURN 213 ENDIF ELSE BEGIN 214 msg = 'iii : valid jpio parameter = ' + STRING(jpio) 215 ras = report(msg) 216 ENDELSE 217 ; 218 ; check jpjo 219 IF (jpjo LT 1) THEN BEGIN 220 msg = ' eee : invalid jpjo parameter = ' + STRING(jpjo) 221 ras = report(msg) 222 msg = ' eee : jpjo must be greater or equal to 1' 223 ras = report(msg) 224 RETURN 225 ENDIF ELSE BEGIN 226 msg = 'iii : valid jpjo parameter = ' + STRING(jpjo) 227 ras = report(msg) 228 ENDELSE 229 ; 230 ; check olon ++ what to check 231 ; when decided reproduce with olat ++ 232 IF (N_ELEMENTS(olon) NE jpio*jpjo) THEN BEGIN 233 msg = ' eee : invalid nb of elements of olon parameter = ' + STRING(N_ELEMENTS(olon)) 234 ras = report(msg) 235 msg = ' eee : nb of elements of olon must be equal to jpio*jpjo = ' + STRING(jpio*jpjo) 236 ras = report(msg) 237 RETURN 238 ENDIF ELSE BEGIN 239 msg = 'iii : valid nb of elements olon parameter' 240 ras = report(msg) 241 ENDELSE 242 ; 243 ; check values ++ 244 ; test if ${GEOMAG_OD} defined 245 geomag_od_env=GETENV('GEOMAG_OD') 246 CASE geomag_od_env OF 247 '' : BEGIN 248 msg = 'eee : ${GEOMAG_OD} is not defined' 249 ras = report(msg) 250 RETURN 251 END 252 ELSE: BEGIN 169 ENDCASE 170 CASE long_name OF 171 'magnetic field': BEGIN 172 msg = 'iii : valid long_name parameter = ' + long_name 173 ras = report(msg) 174 END 175 ELSE : BEGIN 176 msg = 'eee : invalid long_name parameter = ' + long_name 177 ras = report(msg) 178 msg = 'eee : long_name must be magnetic field' 179 ras = report(msg) 180 RETURN 181 END 182 ENDCASE 183 CASE units OF 184 'tesla': BEGIN 185 msg = 'iii : valid units parameter = ' + units 186 ras = report(msg) 187 END 188 ELSE : BEGIN 189 msg = 'eee : invalid units parameter = ' + units 190 ras = report(msg) 191 msg = 'eee : units must be tesla' 192 ras = report(msg) 193 RETURN 194 END 195 ENDCASE 196 END 197 ELSE : BEGIN 198 msg = 'eee : invalid variable.NAME parameter = ' + variable.NAME 199 ras = report(msg) 200 msg = 'eee : variable.NAME must be cond_sed or Br' 201 ras = report(msg) 202 RETURN 203 END 204 ENDCASE 205 ; 206 ; check jpio 207 IF (jpio LT 1) THEN BEGIN 208 msg = ' eee : invalid jpio parameter = ' + STRING(jpio) 209 ras = report(msg) 210 msg = ' eee : jpio must be greater or equal to 1' 211 ras = report(msg) 212 RETURN 213 ENDIF ELSE BEGIN 214 msg = 'iii : valid jpio parameter = ' + STRING(jpio) 215 ras = report(msg) 216 ENDELSE 217 ; 218 ; check jpjo 219 IF (jpjo LT 1) THEN BEGIN 220 msg = ' eee : invalid jpjo parameter = ' + STRING(jpjo) 221 ras = report(msg) 222 msg = ' eee : jpjo must be greater or equal to 1' 223 ras = report(msg) 224 RETURN 225 ENDIF ELSE BEGIN 226 msg = 'iii : valid jpjo parameter = ' + STRING(jpjo) 227 ras = report(msg) 228 ENDELSE 229 ; 230 ; check olon ++ what to check 231 ; when decided reproduce with olat ++ 232 IF (N_ELEMENTS(olon) NE jpio*jpjo) THEN BEGIN 233 msg = ' eee : invalid nb of elements of olon parameter = ' + STRING(N_ELEMENTS(olon)) 234 ras = report(msg) 235 msg = ' eee : nb of elements of olon must be equal to jpio*jpjo = ' + STRING(jpio*jpjo) 236 ras = report(msg) 237 RETURN 238 ENDIF ELSE BEGIN 239 msg = 'iii : valid nb of elements olon parameter' 240 ras = report(msg) 241 ENDELSE 242 ; 243 ; check values ++ 244 ; test if ${GEOMAG_OD} defined 245 geomag_od_env=GETENV('GEOMAG_OD') 246 CASE geomag_od_env OF 247 '' : BEGIN 248 msg = 'eee : ${GEOMAG_OD} is not defined' 249 ras = report(msg) 250 RETURN 251 END 252 ELSE : BEGIN 253 253 msg = 'iii : ${GEOMAG_OD} is ' + geomag_od_env 254 254 ras = report(msg) 255 END 256 ENDCASE 257 ; 258 ; check if output data will be possible 259 iodirout = isadirectory(geomag_od_env) 260 ; 261 ; existence and protection 262 IF (FILE_TEST(iodirout, /DIRECTORY,/WRITE) EQ 0) THEN BEGIN 263 msg = 'eee : the directory' + iodirout + ' was not found.' 264 ras = report(msg) 265 RETURN 266 ENDIF 267 ; 268 ; build output filenames 269 filename = variable.NAME + '_' + orcares +'.nc' 270 fullfilename = iodirout + filename 271 ; 272 ; in order to avoid unexpected overwritten 273 IF (FILE_TEST(fullfilename) EQ 1) THEN BEGIN 274 msg = 'eee : the file ' + fullfilename + ' already exists.' 275 ras = report(msg) 276 RETURN 277 ENDIF 278 ; 279 ; d'après ncdump -h /usr/work/sur/fvi/OPA/geomag/condmag.nc 280 ; ces lignes n'ont rien a faire la ++ ces information devrait etre en parameter 281 condmaglonname = 'lo' 282 condmaglatname = 'la' 283 varname = variable.NAME 284 ; 285 ;--------------------------- 286 ; Creation of the NetCdf output file 287 ;--------------------------- 288 ; 289 netcdf_id = NCDF_CREATE(fullfilename, /clobber) 290 NCDF_CONTROL, netcdf_id, /NOFILL 291 ; 292 ; dimension 293 dimidx = NCDF_DIMDEF(netcdf_id, 'x' , jpio) 294 dimidy = NCDF_DIMDEF(netcdf_id, 'y' , jpjo) 295 dimidt = NCDF_DIMDEF(netcdf_id, 'lo', /UNLIMITED) 296 ; 297 ; global attributes 298 NCDF_ATTPUT, netcdf_id, 'Conventions', 'GDT 1.2', /GLOBAL ; ++ conformite 299 NCDF_ATTPUT, netcdf_id, 'file_name' , filename, /GLOBAL 300 NCDF_ATTPUT, netcdf_id, 'Title' , title, /GLOBAL 301 ; 302 ; declaration of variables 303 ; 4 common variables for the two files to produce 304 varid = lonarr(4) 305 ; 306 varid[0] = NCDF_VARDEF(netcdf_id, 'nav_lon' , [dimidx, dimidy], /FLOAT) 307 NCDF_ATTPUT, netcdf_id, varid[0], 'units' , 'degrees_east' 308 NCDF_ATTPUT, netcdf_id, varid[0], 'valid_min', min(olon, max = omax), /FLOAT 309 NCDF_ATTPUT, netcdf_id, varid[0], 'valid_max', omax, /FLOAT 310 NCDF_ATTPUT, netcdf_id, varid[0], 'long_name', 'Longitude at t-point' 311 ; 312 varid[1] = NCDF_VARDEF(netcdf_id, 'nav_lat' , [dimidx, dimidy], /FLOAT) 313 NCDF_ATTPUT, netcdf_id, varid[1], 'units' , 'degrees_north' 314 NCDF_ATTPUT, netcdf_id, varid[1], 'valid_min', min(olat, max = omax), /FLOAT 315 NCDF_ATTPUT, netcdf_id, varid[1], 'valid_max', omax, /FLOAT 316 NCDF_ATTPUT, netcdf_id, varid[1], 'long_name', 'Latitude at t-point' 317 ; 318 varid[2] = NCDF_VARDEF(netcdf_id, 'time' , [dimidt], /FLOAT) 319 ; 320 varid[3] = NCDF_VARDEF(netcdf_id, varname, [dimidx, dimidy, dimidt], /FLOAT) 321 NCDF_ATTPUT, netcdf_id, varid[3], 'long_name', long_name 322 NCDF_ATTPUT, netcdf_id, varid[3], 'units', units 323 ; pour min et max, il faut avoir lu la variable ... cf. plus bas ++ 324 ; donc pour l'instant on les met à valeur manquante 325 NCDF_ATTPUT, netcdf_id, varid[3], 'valid_min', !VALUES.F_NAN, /FLOAT 326 NCDF_ATTPUT, netcdf_id, varid[3], 'valid_max', !VALUES.F_NAN, /FLOAT 327 ; 328 ;--------------------------- 329 ; end of header definition, writing of the NetCdf files 330 ;--------------------------- 331 ; 332 NCDF_CONTROL, netcdf_id, /ENDEF 333 ; 334 NCDF_VARPUT, netcdf_id, 'nav_lon', olon 335 NCDF_VARPUT, netcdf_id, 'nav_lat', olat 336 NCDF_VARPUT, netcdf_id, varid[2], FLOAT(0.5) ; ++ c'est quoi cette valeur 337 ; 338 ; compute min max 339 minarr = min(values, max = maxarr) 340 ; 341 ; put back the masked value 342 ;++ IF bad[0] NE -1 THEN values[TEMPORARY(bad)] = 32767 343 ; 344 ; write the data 345 NCDF_VARPUT, netcdf_id, varname, values, COUNT = [jpio, jpjo, 1], OFFSET = [0, 0, 0] 346 ; 347 ; update min max attributes 348 NCDF_CONTROL, netcdf_id, /REDEF 349 NCDF_ATTPUT, netcdf_id, varid[3], 'valid_min', minarr, /FLOAT 350 NCDF_ATTPUT, netcdf_id, varid[3], 'valid_max', maxarr, /FLOAT 351 NCDF_CONTROL, netcdf_id, /ENDEF 352 ; 353 ; close the netcdf files 354 NCDF_CLOSE, netcdf_id 355 ; 356 msg = 'iii : ' + fullfilename + ' created' 357 ras = report(msg) 358 ; 255 END 256 ENDCASE 257 ; 258 ; check if output data will be possible 259 iodirout = isadirectory(geomag_od_env) 260 ; 261 ; existence and protection 262 IF (FILE_TEST(iodirout, /DIRECTORY,/WRITE) EQ 0) THEN BEGIN 263 msg = 'eee : the directory' + iodirout + ' was not found.' 264 ras = report(msg) 265 RETURN 266 ENDIF 267 ; 268 ; build output filenames 269 filename = variable.NAME + '_' + orcares +'.nc' 270 fullfilename = iodirout + filename 271 ; 272 ; in order to avoid unexpected overwritten 273 IF (FILE_TEST(fullfilename) EQ 1) THEN BEGIN 274 msg = 'eee : the file ' + fullfilename + ' already exists.' 275 ras = report(msg) 276 RETURN 277 ENDIF 278 ; 279 ; d'aprÚs ncdump -h /usr/work/sur/fvi/OPA/geomag/condmag.nc 280 ; ces lignes n'ont rien a faire la ++ ces information devrait etre en parameter 281 condmaglonname = 'lo' 282 condmaglatname = 'la' 283 varname = variable.NAME 284 ; 285 ;--------------------------- 286 ; Creation of the NetCdf output file 287 ;--------------------------- 288 ; 289 netcdf_id = NCDF_CREATE(fullfilename, /clobber) 290 NCDF_CONTROL, netcdf_id, /NOFILL 291 ; 292 ; dimension 293 dimidx = NCDF_DIMDEF(netcdf_id, 'x' , jpio) 294 dimidy = NCDF_DIMDEF(netcdf_id, 'y' , jpjo) 295 dimidt = NCDF_DIMDEF(netcdf_id, 'lo', /UNLIMITED) 296 ; 297 ; global attributes 298 NCDF_ATTPUT, netcdf_id, 'Conventions', 'GDT 1.2', /GLOBAL ; ++ conformite 299 NCDF_ATTPUT, netcdf_id, 'file_name' , filename, /GLOBAL 300 NCDF_ATTPUT, netcdf_id, 'Title' , title, /GLOBAL 301 ; 302 ; declaration of variables 303 ; 4 common variables for the two files to produce 304 varid = lonarr(4) 305 ; 306 varid[0] = NCDF_VARDEF(netcdf_id, 'nav_lon' , [dimidx, dimidy], /FLOAT) 307 NCDF_ATTPUT, netcdf_id, varid[0], 'units' , 'degrees_east' 308 NCDF_ATTPUT, netcdf_id, varid[0], 'valid_min', min(olon, max = omax), /FLOAT 309 NCDF_ATTPUT, netcdf_id, varid[0], 'valid_max', omax, /FLOAT 310 NCDF_ATTPUT, netcdf_id, varid[0], 'long_name', 'Longitude at t-point' 311 ; 312 varid[1] = NCDF_VARDEF(netcdf_id, 'nav_lat' , [dimidx, dimidy], /FLOAT) 313 NCDF_ATTPUT, netcdf_id, varid[1], 'units' , 'degrees_north' 314 NCDF_ATTPUT, netcdf_id, varid[1], 'valid_min', min(olat, max = omax), /FLOAT 315 NCDF_ATTPUT, netcdf_id, varid[1], 'valid_max', omax, /FLOAT 316 NCDF_ATTPUT, netcdf_id, varid[1], 'long_name', 'Latitude at t-point' 317 ; 318 varid[2] = NCDF_VARDEF(netcdf_id, 'time' , [dimidt], /FLOAT) 319 ; 320 varid[3] = NCDF_VARDEF(netcdf_id, varname, [dimidx, dimidy, dimidt], /FLOAT) 321 NCDF_ATTPUT, netcdf_id, varid[3], 'long_name', long_name 322 NCDF_ATTPUT, netcdf_id, varid[3], 'units', units 323 ; pour min et max, il faut avoir lu la variable ... cf. plus bas ++ 324 ; donc pour l'instant on les met à valeur manquante 325 NCDF_ATTPUT, netcdf_id, varid[3], 'valid_min', !VALUES.F_NAN, /FLOAT 326 NCDF_ATTPUT, netcdf_id, varid[3], 'valid_max', !VALUES.F_NAN, /FLOAT 327 ; 328 ;--------------------------- 329 ; end of header definition, writing of the NetCdf files 330 ;--------------------------- 331 ; 332 NCDF_CONTROL, netcdf_id, /ENDEF 333 ; 334 NCDF_VARPUT, netcdf_id, 'nav_lon', olon 335 NCDF_VARPUT, netcdf_id, 'nav_lat', olat 336 NCDF_VARPUT, netcdf_id, varid[2], FLOAT(0.5) ; ++ c'est quoi cette valeur 337 ; 338 ; compute min max 339 minarr = min(values, max = maxarr) 340 ; 341 ; put back the masked value 342 ;++ IF bad[0] NE -1 THEN values[TEMPORARY(bad)] = 32767 343 ; 344 ; write the data 345 NCDF_VARPUT, netcdf_id, varname, values, COUNT = [jpio, jpjo, 1], OFFSET = [0, 0, 0] 346 ; 347 ; update min max attributes 348 NCDF_CONTROL, netcdf_id, /REDEF 349 NCDF_ATTPUT, netcdf_id, varid[3], 'valid_min', minarr, /FLOAT 350 NCDF_ATTPUT, netcdf_id, varid[3], 'valid_max', maxarr, /FLOAT 351 NCDF_CONTROL, netcdf_id, /ENDEF 352 ; 353 ; close the netcdf files 354 NCDF_CLOSE, netcdf_id 355 ; 356 msg = 'iii : ' + fullfilename + ' created' 357 ras = report(msg) 359 358 END -
trunk/correc_angle.pro
r41 r48 12 12 ; zresul_vv, zresul_vu, zresul_u, zresul_v 13 13 ; 14 ; INPUTS: 15 ; zresul_uu : x-component of interpolated field at 16 ; u-point 17 ; zresul_uv : x-component of interpolated field at 18 ; v-point 19 ; zresul_vv : y-component of interpolated field at 20 ; v-point 21 ; zresul_vu : y-component of interpolated field at 22 ; u-point 14 ; INPUTS: 15 ; zresul_uu : x-component of interpolated field at 16 ; u-point 17 ; zresul_uv : x-component of interpolated field at 18 ; v-point 19 ; zresul_vv : y-component of interpolated field at 20 ; v-point 21 ; zresul_vu : y-component of interpolated field at 22 ; u-point 23 23 ; 24 24 ; KEYWORD PARAMETERS: None 25 25 ; 26 ; OUTPUTS: 27 ; zresul_u : corrected x-component of interpolated 26 ; OUTPUTS: 27 ; zresul_u : corrected x-component of interpolated 28 28 ; field at u-point 29 ; zresul_v : corrected y-component of interpolated 29 ; zresul_v : corrected y-component of interpolated 30 30 ; field at v-point 31 31 ; 32 32 ; COMMON BLOCKS: 33 ; 33 ; common_interp.pro 34 34 ; 35 ; SIDE EFFECTS: 35 ; SIDE EFFECTS: 36 36 ; 37 37 ; RESTRICTIONS: … … 55 55 angle, lon_u, lat_u, lon_v, lat_v, lon_f, lat_f, gsinu, gsinv, gcosu, gcosv 56 56 ; 57 ; 2. Correct the two components of the field 57 ; 2. Correct the two components of the field 58 58 ; ========================================== 59 59 ; … … 68 68 ; 69 69 return 70 END 70 END -
trunk/detectbarotropicmode.sh
r47 r48 16 16 # ======== 17 17 # 18 # .. code-bl cok:: bash18 # .. code-block:: bash 19 19 # 20 20 # detectbarotropicmode.sh -l list -r resolution -exp experience -g grid \ … … 31 31 # in the box [-50,30] [295,335]: 32 32 # 33 # .. code-block 33 # .. code-block:: bash 34 34 # 35 35 # detectbarotropicmode.sh -l ${GEOMAG_OD}/list_G42 -r ORCA025 -exp G42 -g gridT \ … … 60 60 command=detectbarotropicmode.sh 61 61 # 62 usage=" 63 # 64 set +u65 while [ ! -z "${1}"]62 usage="Usage : ${command} -l list -r orcares -exp drakkar_exp -g grid -latmin latmin -latmax latmax -longmin longmin -longmax longmax" 63 # 64 set -u 65 while [ ${#} -gt 0 ] 66 66 do 67 case ${1} in68 -l)69 list=${2}70 shift71 ;;72 -r)73 orcares=${2}74 shift75 ;;76 -exp)77 drakkar_exp=${2}78 shift79 ;;80 -g)81 grid=${2}82 shift83 ;;84 -latmin)85 latmin=${2}86 shift87 ;;88 -latmax)89 latmax=${2}90 shift91 ;;92 -longmin)93 longmin=${2}94 shift95 ;;96 -longmax)97 longmax=${2}98 shift99 ;;100 *)101 # other choice102 echo "eee : unknown option ${1}"103 echo "${usage}"104 ;;105 esac106 # next flag107 shift67 case ${1} in 68 -l) 69 list=${2} 70 shift 71 ;; 72 -r) 73 orcares=${2} 74 shift 75 ;; 76 -exp) 77 drakkar_exp=${2} 78 shift 79 ;; 80 -g) 81 grid=${2} 82 shift 83 ;; 84 -latmin) 85 latmin=${2} 86 shift 87 ;; 88 -latmax) 89 latmax=${2} 90 shift 91 ;; 92 -longmin) 93 longmin=${2} 94 shift 95 ;; 96 -longmax) 97 longmax=${2} 98 shift 99 ;; 100 *) 101 # other choice 102 echo "eee : unknown option ${1}" 103 echo "${usage}" 104 ;; 105 esac 106 # next flag 107 shift 108 108 done 109 #110 set -u111 109 # 112 110 # check GEOMAG environement 113 111 if [ ! -d ${GEOMAG_LOG} ] 114 112 then 115 echo " eee : \${GEOMAG_LOG} not found"116 exit 1113 echo " eee : \${GEOMAG_LOG} not found" 114 exit 1 117 115 fi 118 116 # ++ blindage en tout genre 119 117 # 120 118 case ${grid} in 121 gridT)122 echo "iii : grid=${orcares}"123 ;;124 *)125 echo "eee : pb \${grid} = ${grid}"126 exit 1127 ;;119 gridT) 120 echo "iii : grid=${orcares}" 121 ;; 122 *) 123 echo "eee : pb \${grid} = ${grid}" 124 exit 1 125 ;; 128 126 esac 129 127 # 130 128 case ${orcares} in 131 ORCA025)132 echo "iii : orcares=${orcares}"133 ;;134 *)135 echo "eee : pb \${orcares} = ${orcares}"136 exit 1137 ;;129 ORCA025) 130 echo "iii : orcares=${orcares}" 131 ;; 132 *) 133 echo "eee : pb \${orcares} = ${orcares}" 134 exit 1 135 ;; 138 136 esac 139 137 # 140 138 case ${drakkar_exp} in 141 G42)142 echo "iii : drakkar_exp=${drakkar_exp}"143 ;;144 *)145 echo "eee : pb \${drakkar_exp} = ${drakkar_exp}"146 exit 1147 ;;139 G42) 140 echo "iii : drakkar_exp=${drakkar_exp}" 141 ;; 142 *) 143 echo "eee : pb \${drakkar_exp} = ${drakkar_exp}" 144 exit 1 145 ;; 148 146 esac 149 147 # 150 148 # ++ conversion of latitude and longitude in x and y indexes 151 149 case ${orcares} in 152 ORCA025)153 xmin=529154 xmax=689155 ymin=266156 ymax=372157 echo "iii : forcage xmin xmax ymin ymax : ${xmin} ${xmax} ${ymin} ${ymax}"158 ;;159 *)160 echo "eee : pb conversion lat,long to x,y"161 exit 1162 ;;150 ORCA025) 151 xmin=529 152 xmax=689 153 ymin=266 154 ymax=372 155 echo "iii : forcage xmin xmax ymin ymax : ${xmin} ${xmax} ${ymin} ${ymax}" 156 ;; 157 *) 158 echo "eee : pb conversion lat,long to x,y" 159 exit 1 160 ;; 163 161 esac 164 162 # … … 167 165 for file in `cat ${list}` 168 166 do 169 echo "iii : extraction from ${file}" # ++ voir seulement si debug demandé170 filessh=${GEOMAG_OD}/`basename ${file} .nc`_ssh.nc171 rm -f ${filessh} 2> /dev/null172 ncks --dimension x,${xmin},${xmax} --dimension y,${ymin},${ymax} -v sossheig,nav_lon,nav_lat \173 174 status_ncks=${?}175 if [ ${status_ncks} -ne 0 ]176 then177 echo "eee : pb with ncks ${file}"178 exit 1179 else180 list_ssh="${list_ssh} ${filessh}"181 fi167 echo "iii : extraction from ${file}" # ++ voir seulement si debug demandé 168 filessh=${GEOMAG_OD}/`basename ${file} .nc`_ssh.nc 169 rm -f ${filessh} 2> /dev/null 170 ncks --dimension x,${xmin},${xmax} --dimension y,${ymin},${ymax} -v sossheig,nav_lon,nav_lat \ 171 ${file} ${filessh} 172 status_ncks=${?} 173 if [ ${status_ncks} -ne 0 ] 174 then 175 echo "eee : pb with ncks ${file}" 176 exit 1 177 else 178 list_ssh="${list_ssh} ${filessh}" 179 fi 182 180 done 183 181 line=`head -n 1 ${list}` … … 195 193 if [ ${status_ncrcat} -ne 0 ] 196 194 then 197 echo "eee : pb with ncrcat"198 exit 1195 echo "eee : pb with ncrcat" 196 exit 1 199 197 else 200 echo "iii : result in ${filetot}"198 echo "iii : result in ${filetot}" 201 199 fi 202 200 # -
trunk/divfred.pro
r41 r48 25 25 ; - Requires SAXO tools 26 26 ; - les matrices u et v peuvent de 2 a 4 dimensions. 27 ; attention pour distinguer les diff érentes configurations de u et v27 ; attention pour distinguer les différentes configurations de u et v 28 28 ; (xy, xyz, xyt, xyzt), on regarde la variable du common 29 29 ; -time qui contient le calendrier en jour julien d''IDL auquel 30 30 ; se rapportent u et v ansi que la variable 31 ; -jpt qui est le nombre de pas de temps à considérer dans time.32 ; les tableaux u et v sont d écoupés sur le même domaine33 ; g éographique. A cause du décalage des grilles T, U, V et F il est34 ; possible que ces 2 tableaux n''aient pas la m ême taille et se35 ; reportent à des indices différents. Si tel est le cas les tableaux36 ; sont red écoupés sur les indices qu'ils ont en commun et le dommaine37 ; est red éfini pour qu'il colle àces indices communs.38 ; pour éviter ces redécoupes utiliser le mot clé/memeindice dans31 ; -jpt qui est le nombre de pas de temps à considérer dans time. 32 ; les tableaux u et v sont découpés sur le même domaine 33 ; géographique. A cause du décalage des grilles T, U, V et F il est 34 ; possible que ces 2 tableaux n''aient pas la même taille et se 35 ; reportent à des indices différents. Si tel est le cas les tableaux 36 ; sont redécoupés sur les indices qu'ils ont en commun et le dommaine 37 ; est redéfini pour qu'il colle à ces indices communs. 38 ; pour éviter ces redécoupes utiliser le mot clé /memeindice dans 39 39 ; <pro>domdef</pro> 40 40 ; 41 ; les points sur le bord du dessin sont mis à!values.f_nan41 ; les points sur le bord du dessin sont mis à !values.f_nan 42 42 ; 43 43 ; @todo 44 ; ++ pas fini de comprendre, tester (compatibilit ésaxo), adapter, commenter44 ; ++ pas fini de comprendre, tester (compatibilité saxo), adapter, commenter 45 45 ; nettoyer 46 46 ; ++ a comparer et merger avec SAXO_DIR/ToBeReviewed/CALCULS/div.pro … … 111 111 v = v[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, *] 112 112 END 113 ELSE :BEGIN113 ELSE :BEGIN 114 114 zdiv = -1 115 115 GOTO, sortie … … 142 142 ; mise a !values.f_nan de la bordure 143 143 ;------------------------------------------------------------ 144 if 144 if NOT keyword_set(key_periodic) OR nx NE jpi then begin 145 145 zdiv(0, *, *) = !values.f_nan 146 146 zdiv(nx-1, *, *) = !values.f_nan … … 161 161 varunits = '1e6*s-1' 162 162 domdef, (glamt[indice2d])[0, 0], (glamu[indice2d])[nx-1, 0],(gphit[indice2d])[0, 0], (gphiv[indice2d])[0, ny-1], vert1, vert2, grille = ['t'] 163 if keyword_set(direc) then 163 if keyword_set(direc) then zdiv = moyenne(zdiv,direc,/nan, boite = boite) 164 164 END 165 165 ;---------------------------------------------------------------------------- … … 189 189 v = v[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, *] 190 190 END 191 ELSE :return, -1191 ELSE : return, -1 192 192 endcase 193 193 ;------------------------------------------------------------ … … 216 216 ; mise a !values.f_nan de la bordure 217 217 ;------------------------------------------------------------ 218 if 218 if NOT keyword_set(key_periodic) OR nx NE jpi then begin 219 219 zdiv(0, *, *) = !values.f_nan 220 220 zdiv(nx-1, *, *) = !values.f_nan … … 232 232 varunits = '1e6*s-1' 233 233 domdef, (glamt[indice2d])[0, 0], (glamu[indice2d])[nx-1, 0],(gphit[indice2d])[0, 0], (gphiv[indice2d])[0, ny-1], vert1, vert2, grille = ['t'] 234 if keyword_set(direc) then 234 if keyword_set(direc) then zdiv = grossemoyenne(zdiv,direc,/nan, boite = boite) 235 235 END 236 236 ;---------------------------------------------------------------------------- … … 247 247 ;---------------------------------------------------------------------------- 248 248 ;---------------------------------------------------------------------------- 249 ELSE :BEGIN ;xy249 ELSE : BEGIN ;xy 250 250 indice3d = lindgen(jpi, jpj, jpk) 251 251 indice3d = indice3d[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, firstzt] … … 273 273 v = v[indice2d] 274 274 END 275 ELSE :return, -1275 ELSE : return, -1 276 276 endcase 277 277 ;------------------------------------------------------------ … … 289 289 ; mise a !values.f_nan de la bordure 290 290 ;------------------------------------------------------------ 291 if 291 if NOT keyword_set(key_periodic) OR nx NE jpi then begin 292 292 zdiv(0, *) = !values.f_nan 293 293 zdiv(nx-1, *) = !values.f_nan … … 308 308 varunits = '1e6*s-1' 309 309 domdef, (glamt[indice2d])[0, 0], (glamu[indice2d])[nx-1, 0],(gphit[indice2d])[0, 0], (gphiv[indice2d])[0, ny-1], vert1, vert2, grille = ['t'] 310 if keyword_set(direc) then 310 if keyword_set(direc) then zdiv = moyenne(zdiv,direc,/nan, boite = boite) 311 311 312 312 END -
trunk/draw.pro
r2 r48 25 25 IF n_elements(ind) GT 1 THEN field(ind) = 0. 26 26 ; 27 IF size(glam, /n_dimensions) EQ 1 THEN triangles = triangule() ELSE triangulate, glamt, gphit, triangles 27 IF size(glam, /n_dimensions) EQ 1 THEN triangles = triangule() ELSE triangulate, glamt, gphit, triangles 28 28 ; 29 IF keyword_set(land) THEN BEGIN 29 IF keyword_set(land) THEN BEGIN 30 30 if dessin EQ 1 then plt, field, zmin, zmax, nite = 0, /carte, petit = [nbcol, nbdessin, dessin], /remplit, div = 4, /cont_nofill, retain = 2, /nocontour, /land $ 31 31 ELSE plt, field, zmin, zmax, nite = 0, /carte, petit = [nbcol, nbdessin, dessin], $ 32 32 /noerase, /remplit, div = 4, /cont_nofill, /nocontour, /land 33 ENDIF ELSE BEGIN 33 ENDIF ELSE BEGIN 34 34 if dessin EQ 1 then plt, field, zmin, zmax, nite = 0, /carte, petit = [nbcol, nbdessin, dessin], /remplit, div = 4, /cont_nofill, retain = 2, /nocontour, /port $ 35 35 ELSE plt, field, zmin, zmax, nite = 0, /carte, petit = [nbcol, nbdessin, dessin], $ 36 36 /noerase, /remplit, div = 4, /cont_nofill, /nocontour, /port 37 ENDELSE 37 ENDELSE 38 38 return 39 39 end -
trunk/extrap.pro
r41 r48 2 2 ; NAME: extrap.pro (based on remplit.pro) 3 3 ; 4 ; PURPOSE: Extrapolates ocean data values on land points 4 ; PURPOSE: Extrapolates ocean data values on land points 5 5 ; 6 6 ; CATEGORY: Subroutine … … 8 8 ; CALLING SEQUENCE: extrap, zdata, zmask, it, val 9 9 ; 10 ; INPUTS: 10 ; INPUTS: 11 11 ; zdata : field to extrapolate 12 12 ; zmask : field's mask 13 13 ; it : iteration 14 ; val : field value on masked points 14 ; val : field value on masked points 15 15 ; 16 16 ; KEYWORD PARAMETERS: None 17 17 ; 18 ; OUTPUTS: 19 ; zdata : extrapolated field 18 ; OUTPUTS: 19 ; zdata : extrapolated field 20 20 ; zmask : mask after extrapolation 21 21 ; 22 22 ; COMMON BLOCKS: 23 ; common_interp.pro23 ; common_interp.pro 24 24 ; 25 25 ; … … 38 38 ; IF keyword_set(key_performance) THEN print, systime(1)-temp 39 39 ; 40 ; 2. Extrapolating when needed 40 ; 2. Extrapolating when needed 41 41 ; ============================ 42 42 ; … … 71 71 ; IF keyword_set(key_performance) THEN print, systime(1)-temp 72 72 ; IF keyword_set(key_performance) THEN temp = systime(1) 73 ; IF n_elements(ocean2) NE jpiatm*(jpjatm+4) THEN BEGIN 73 ; IF n_elements(ocean2) NE jpiatm*(jpjatm+4) THEN BEGIN 74 74 ; land = where(zmask EQ 0.) 75 75 ; zremplit(land) = val 76 ; ENDIF 76 ; ENDIF 77 77 ; IF keyword_set(key_performance) THEN print, systime(1)-temp 78 78 ; … … 95 95 terre = where(mmmask EQ 0) 96 96 if terre[0] EQ -1 then GOTO, fini 97 ; les points du bord du cadre doivent maintenan etre dans la terre97 ; les points du bord du cadre doivent maintenant etre dans la terre 98 98 mmmask[0, *] = 0 99 99 mmmask[jpiatm+1, *] = 0 … … 119 119 +1./sqrt(2)*(z[(jpiatm+2)+1+cote]+z[(jpiatm+2)-1+cote]+z[-(jpiatm+2)+1+cote]+z[-(jpiatm+2)-1+cote]) 120 120 poids = voisins[cote] 121 121 122 122 z[cote] = zcote/poids 123 123 ;--------------------------------------------------------------- -
trunk/forcage.pro
r2 r48 1 PRO forcagequimarche,iyear,ian 1 PRO forcagequimarche,iyear,ian 2 2 @init2 3 3 @initorca2_bab 4 4 5 6 7 8 5 @common 9 10 11 6 12 7 ;ian='01' … … 16 11 rep_fred='/usr/work/sur/fvi/OPA/geomag/' 17 12 key_portrait = 0 18 ; stockage des fichiers brut 13 ; stockage des fichiers brut 19 14 ioDATA='/usr/work/sur/fvi/OPA/ORCA2/' 20 file_U=e_exp+'_5d_'+ian+'0101_'+ian+'1231_grid_U.nc' 15 file_U=e_exp+'_5d_'+ian+'0101_'+ian+'1231_grid_U.nc' 21 16 print, file_U 22 file_V=e_exp+'_5d_'+ian+'0101_'+ian+'1231_grid_V.nc' 17 file_V=e_exp+'_5d_'+ian+'0101_'+ian+'1231_grid_V.nc' 23 18 print, file_V 24 file_T=e_exp+'_5d_'+ian+'0101_'+ian+'1231_grid_T.nc' 19 file_T=e_exp+'_5d_'+ian+'0101_'+ian+'1231_grid_T.nc' 25 20 print, file_T 26 21 file_Sed= rep_fred+'cond_sed_ORCA2.nc' … … 33 28 t_bt = 'bar_transp' 34 29 ioORLN2 = '/usr/work/sur/fvi/OPA/ORCA2' 35 ;facteur d'echelle vertical 30 ;facteur d'echelle vertical for partial steps 36 31 e3v3d=read_ncdf('e3v_ps',0,/timestep,iodir=ioORLN2,/nostruct,/tout,filename='meshmask_bab.nc') 37 32 e3u3d=read_ncdf('e3u_ps',0,/timestep,iodir=ioORLN2,/nostruct,/tout,filename='meshmask_bab.nc') … … 44 39 45 40 ; vertical integration: 46 e3t3d=make_array(jpi,jpj,jpk) 41 e3t3d=make_array(jpi,jpj,jpk) 47 42 for k=0, jpk-1 do begin &$ 48 for j=0,jpj-1 do begin &$ 43 for j=0,jpj-1 do begin &$ 49 44 for i=0,jpi-1 do begin &$ 50 45 e3t3d(i,j,k) = e3t(k) &$ … … 54 49 jpt = 73 55 50 56 ;vud = make_array(jpi,jpj,jpt) 51 ;vud = make_array(jpi,jpj,jpt) 57 52 ;vvd = make_array(jpi, jpj, jpt) 58 53 divBustar = make_array(jpi, jpj, jpt) 59 54 diver2=fltarr(jpi,jpj,1,73) 60 55 61 ; ouverture des fichiers dans lesquels on va écrire56 ; ouverture des fichiers dans lesquels on va écrire 62 57 ;id3=NCDF_OPEN('/usr/work/sur/fvi/OPA/geomag/U_5d_'+iyear+'_grid_T.nc',/write) 63 58 ;id4=NCDF_OPEN('/usr/work/sur/fvi/OPA/geomag/V_5d_'+iyear+'_grid_T.nc',/write) … … 66 61 ;id4=NCDF_OPEN('/usr/work/sur/fvi/OPA/ORCA2/DivBustar_5d_'+iyear+'_grid_T.nc',/write) 67 62 68 FOR jt = 0, jpt-1 DO BEGIN &$ 63 FOR jt = 0, jpt-1 DO BEGIN &$ 69 64 70 65 ; ouverture des fichiers et stockage en memoire partial steps 71 vu=read_ncdf('vozocrtx',jt,jt, /timestep, iodir=ioDATA,/nostruct,/TOUT,filename=file_U) &$ 72 ;stop 73 vv=read_ncdf('vomecrty',jt,jt, /timestep,iodir=ioDATA,/nostruct,/TOUT,filename=file_V) &$ 66 vu=read_ncdf('vozocrtx',jt,jt, /timestep, iodir=ioDATA,/nostruct,/TOUT,filename=file_U) &$ 67 ;stop 68 vv=read_ncdf('vomecrty',jt,jt, /timestep,iodir=ioDATA,/nostruct,/TOUT,filename=file_V) &$ 74 69 ;stop 75 70 ; lecture salinite & temperature … … 107 102 108 103 conduct_v=(conduct+shift(conduct,0,1,0))/2. 109 104 110 105 u_cond_u=total( vu*conduct_u*e3u3d*umask(),3) 111 106 … … 126 121 ;for jj=0,147 do begin &$ 127 122 ; for ji=0,179 do begin &$ 128 ; if( Diver(ji,jj) GT 1e10 ) then begin &$ 129 ; 123 ; if( Diver(ji,jj) GT 1e10 ) then begin &$ 124 ; Diver(ji,jj) = 0. &$ 130 125 ; endif &$ 131 126 ; endfor &$ 132 ;endfor 127 ;endfor 133 128 134 129 ;stop … … 137 132 ;stop 138 133 ; Somme sur la verticale partial steps 139 ; vum=total( vu*e3u3d*umask(),3 ) &$ 140 ; vvm=total( vv*e3v3d*vmask(),3 ) &$ 134 ; vum=total( vu*e3u3d*umask(),3 ) &$ 135 ; vvm=total( vv*e3v3d*vmask(),3 ) &$ 141 136 142 137 ; Shift sur la grille T partial steps 143 ; vut= (vum+shift(vum,1,0) )*0.5 &$ 138 ; vut= (vum+shift(vum,1,0) )*0.5 &$ 144 139 ; vvt= (vvm+shift(vvm,0,1) )*0.5 &$ 145 ; Bande de recouvrement 140 ; Bande de recouvrement 146 141 ; vut(0, *) = vut(jpi-2, *) 147 142 ; vvt(*, 0) = 0. 148 ; stockage dans le fichier de sortie 143 ; stockage dans le fichier de sortie 149 144 ;NCDF_VARPUT, id3,'sossheig',vut, offset = [0, 0, jt] 150 145 ;NCDF_VARPUT, id4,'sossheig',Diver, offset = [0, 0, jt] … … 153 148 154 149 155 ENDFOR 150 ENDFOR 156 151 ; on ferme le NetCDF 157 152 ;NCDF_CLOSE,id3 … … 218 213 ; Creation de la latitude 219 214 NCDF_VARPUT, idout, id1, gphit 220 ; Creation de la profondeur 221 NCDF_VARPUT, idout, id2, gdept 215 ; Creation de la profondeur 216 NCDF_VARPUT, idout, id2, gdept 222 217 ; Creation du calendrier 223 218 224 219 NCDF_VARPUT, idout, id3, temps 225 220 226 221 227 222 ; Ecriture des donnees 228 223 229 ; ecriture des glam_8 224 ; ecriture des glam_8 230 225 NCDF_VARPUT, idout, id4 , diver2 231 226 -
trunk/forcageforMOED.pro
r2 r48 1 1 PRO forcageforMOED 2 @init2 3 @initorca2_bab 4 5 @common 6 7 ian='01' 8 iyear='1993' 9 10 e_exp='ESS' 11 rep_fred='/usr/work/sur/fvi/OPA/geomag/' 12 key_portrait = 0 13 ; stockage des fichiers brut 14 ioDATA='/usr/work/sur/fvi/OPA/ORCA2/' 15 file_U=e_exp+'_5d_'+ian+'0101_'+ian+'1231_grid_U.nc' 16 print, file_U 17 file_V=e_exp+'_5d_'+ian+'0101_'+ian+'1231_grid_V.nc' 18 print, file_V 19 file_T=e_exp+'_5d_'+ian+'0101_'+ian+'1231_grid_T.nc' 20 print, file_T 21 file_Sed= rep_fred+'cond_sed_ORCA2.nc' 22 file_Br= rep_fred+'Br_ORCA2.nc' 23 24 25 ; title 26 t_exp= e_exp 27 28 t_bt = 'bar_transp' 29 ioORLN2 = '/usr/work/sur/fvi/OPA/ORCA2' 30 ;facteur d'echelle vertical for partial steps 31 e3v3d=read_ncdf('e3v_ps',0,/timestep,iodir=ioORLN2,/nostruct,/tout,filename='meshmask_bab.nc') 32 e3u3d=read_ncdf('e3u_ps',0,/timestep,iodir=ioORLN2,/nostruct,/tout,filename='meshmask_bab.nc') 33 SIGMAsed=read_ncdf('cond_sed',0,/timestep,/nostruct,/tout,filename=file_Sed,/cont_nofill) 34 BR=read_ncdf('Br',0,/timestep,/nostruct,/tout,filename=file_Br,/cont_nofill) 35 36 37 38 39 40 ; vertical integration: 41 e3t3d=make_array(jpi,jpj,jpk) 42 for k=0, jpk-1 do begin &$ 43 for j=0,jpj-1 do begin &$ 2 @init2 3 @initorca2_bab 4 5 @common 6 7 ian='01' 8 iyear='1993' 9 10 e_exp='ESS' 11 rep_fred='/usr/work/sur/fvi/OPA/geomag/' 12 key_portrait = 0 13 ; stockage des fichiers brut 14 ioDATA='/usr/work/sur/fvi/OPA/ORCA2/' 15 file_U=e_exp+'_5d_'+ian+'0101_'+ian+'1231_grid_U.nc' 16 print, file_U 17 file_V=e_exp+'_5d_'+ian+'0101_'+ian+'1231_grid_V.nc' 18 print, file_V 19 file_T=e_exp+'_5d_'+ian+'0101_'+ian+'1231_grid_T.nc' 20 print, file_T 21 file_Sed= rep_fred+'cond_sed_ORCA2.nc' 22 file_Br= rep_fred+'Br_ORCA2.nc' 23 24 25 ; title 26 t_exp= e_exp 27 28 t_bt = 'bar_transp' 29 ioORLN2 = '/usr/work/sur/fvi/OPA/ORCA2' 30 ;facteur d'echelle vertical for partial steps 31 e3v3d=read_ncdf('e3v_ps',0,/timestep,iodir=ioORLN2,/nostruct,/tout,filename='meshmask_bab.nc') 32 e3u3d=read_ncdf('e3u_ps',0,/timestep,iodir=ioORLN2,/nostruct,/tout,filename='meshmask_bab.nc') 33 SIGMAsed=read_ncdf('cond_sed',0,/timestep,/nostruct,/tout,filename=file_Sed,/cont_nofill) 34 BR=read_ncdf('Br',0,/timestep,/nostruct,/tout,filename=file_Br,/cont_nofill) 35 36 ; vertical integration: 37 e3t3d=make_array(jpi,jpj,jpk) 38 for k=0, jpk-1 do begin &$ 39 for j=0,jpj-1 do begin &$ 44 40 for i=0,jpi-1 do begin &$ 45 e3t3d(i,j,k) = e3t(k) &$41 e3t3d(i,j,k) = e3t(k) &$ 46 42 endfor &$ 47 endfor &$ 48 endfor 49 jpt = 73 50 51 ;vud = make_array(jpi,jpj,jpt) 52 ;vvd = make_array(jpi, jpj, jpt) 53 divBustar = make_array(jpi, jpj, jpt) 54 diver3=replicate(0.,182,149,1,73) 55 sigma3=replicate(0.,182,149,1,73) 56 Jpu=replicate(0.,182,149,1,73) 57 Jpv=replicate(0.,182,149,1,73) 58 59 60 FOR jt = 0, jpt-1 DO BEGIN &$ 61 62 ; ouverture des fichiers et stockage en memoire partial steps 63 vu=read_ncdf('vozocrtx',jt,jt, /timestep, iodir=ioDATA,/nostruct,/TOUT,filename=file_U) &$ 64 ;stop 65 vv=read_ncdf('vomecrty',jt,jt, /timestep,iodir=ioDATA,/nostruct,/TOUT,filename=file_V) &$ 66 ;stop 67 ; lecture salinite & temperature 68 temp= read_ncdf('votemper',jt,jt,/timestep,iodir=ioDATA,/nostruct,/tout,filename=file_T) 69 ;stop 70 salin=read_ncdf('vosaline',jt,jt,/timestep,iodir=ioDATA,/nostruct,/tout,filename=file_T) 71 ;stop 72 conduct=0.02047780622061 + 0.00273147624197*temp + 0.00035133182334*temp*temp + 0.09139808809909*salin + 0.00241425798890*salin*temp -0.00023998958774*salin*salin 73 mask_t=where(conduct GT 1.e+19) 74 conduct(mask_t)=0. 75 ; Somme conduct au point T 76 77 78 ; Calcul SIGMA 79 80 SIGMAoc=total(conduct*e3t3d*tmask,3) 81 SIGMA=SIGMAsed+SIGMAoc 82 83 84 85 SIGMA_u=(SIGMA+shift(SIGMA,-1,0))/2. 86 87 SIGMA_v=(SIGMA+shift(SIGMA,0,1))/2. 88 89 ; Calcul B en points u et v 90 91 BR_u=(BR+shift(BR,-1,0))/2. 92 93 BR_v=(BR+shift(BR,0,1))/2. 94 95 96 ; Calcul integrale conduct 97 98 conduct_u=(conduct+shift(conduct,-1,0,0))/2. 99 100 conduct_v=(conduct+shift(conduct,0,1,0))/2. 101 102 u_cond_u=total( vu*conduct_u*e3u3d*umask(),3) 103 104 v_cond_v=total( vv*conduct_v*e3v3d*vmask(),3) 105 106 107 ;calcul de Jp= u_star x Fz k 108 Jpv_star= -1*BR_u*u_cond_u 109 110 Jpu_star= BR_v*v_cond_v 111 112 ; Divergence du champ 113 114 Diver=divfred(Jpu_star,Jpv_star) 115 116 Diver=Diver*1e-6 117 ;stop 118 lecontinent=where(Diver GT 1.E08) 119 Diver(lecontinent)=0. 120 121 122 123 ;stop 124 ;bande de recouvrement:: 125 126 diver3(1:180,0:147,*,jt)=Diver(*,*) 127 diver3(0,*,*,*)=diver3(180,*,*,*) 128 diver3(181,*,*,*)=diver3(1,*,*,*) 129 130 sigma3(1:180,0:147,*,jt)=SIGMA(*,*) 131 sigma3(0,*,*,*)=sigma3(180,*,*,*) 132 sigma3(181,*,*,*)=sigma3(1,*,*,*) 133 134 Jpu(1:180,0:147,*,jt)=Jpu_star(*,*) 135 Jpu(0,*,*,*)=Jpu(180,*,*,*) 136 Jpu(181,*,*,*)=Jpu(1,*,*,*) 137 138 Jpv(1:180,0:147,*,jt)=Jpv_star(*,*) 139 Jpv(0,*,*,*)=Jpv(180,*,*,*) 140 Jpv(181,*,*,*)=Jpv(1,*,*,*) 141 142 143 144 145 print, jt 146 147 148 ENDFOR 149 150 151 152 153 temps=fltarr(73) 154 temps(0)=0. 155 for jt=0,71 do begin &$ 156 temps(jt+1)=temps(jt) +5*86400. &$ 157 endfor 158 print,temps 159 160 vargrid = 'T' 161 iodir = '/usr/work/sur/fvi/OPA/ORCA2/' 162 ; Nom 163 idout = NCDF_CREATE(iodir+'ForcageMOED_5d_'+iyear+'_grid_T.nc',/clobber) 164 print, 'Creation du fichier Netcdf' 165 NCDF_CONTROL, idout, /nofill 166 ; Dimension 167 xidout = NCDF_DIMDEF(idout, 'x',jpiglo) 168 yidout = NCDF_DIMDEF(idout, 'y',jpjglo) 169 didout = NCDF_DIMDEF(idout, 'deptht',1) 170 tidout = NCDF_DIMDEF(idout, 'time_counter', /unlimited) 171 172 didout1 = NCDF_DIMDEF(idout, 'deptht1',jpk) 43 endfor &$ 44 endfor 45 jpt = 73 46 47 ;vud = make_array(jpi,jpj,jpt) 48 ;vvd = make_array(jpi, jpj, jpt) 49 divBustar = make_array(jpi, jpj, jpt) 50 diver3=replicate(0.,182,149,1,73) 51 sigma3=replicate(0.,182,149,1,73) 52 Jpu=replicate(0.,182,149,1,73) 53 Jpv=replicate(0.,182,149,1,73) 54 55 56 FOR jt = 0, jpt-1 DO BEGIN &$ 57 58 ; ouverture des fichiers et stockage en memoire partial steps 59 vu=read_ncdf('vozocrtx',jt,jt, /timestep, iodir=ioDATA,/nostruct,/TOUT,filename=file_U) &$ 60 ;stop 61 vv=read_ncdf('vomecrty',jt,jt, /timestep,iodir=ioDATA,/nostruct,/TOUT,filename=file_V) &$ 62 ;stop 63 ; lecture salinite & temperature 64 temp= read_ncdf('votemper',jt,jt,/timestep,iodir=ioDATA,/nostruct,/tout,filename=file_T) 65 ;stop 66 salin=read_ncdf('vosaline',jt,jt,/timestep,iodir=ioDATA,/nostruct,/tout,filename=file_T) 67 ;stop 68 conduct=0.02047780622061 + 0.00273147624197*temp + 0.00035133182334*temp*temp + 0.09139808809909*salin + 0.00241425798890*salin*temp -0.00023998958774*salin*salin 69 mask_t=where(conduct GT 1.e+19) 70 conduct(mask_t)=0. 71 ; Somme conduct au point T 72 73 74 ; Calcul SIGMA 75 76 SIGMAoc=total(conduct*e3t3d*tmask,3) 77 SIGMA=SIGMAsed+SIGMAoc 78 79 SIGMA_u=(SIGMA+shift(SIGMA,-1,0))/2. 80 81 SIGMA_v=(SIGMA+shift(SIGMA,0,1))/2. 82 83 ; Calcul B en points u et v 84 85 BR_u=(BR+shift(BR,-1,0))/2. 86 87 BR_v=(BR+shift(BR,0,1))/2. 88 89 90 ; Calcul integrale conduct 91 92 conduct_u=(conduct+shift(conduct,-1,0,0))/2. 93 94 conduct_v=(conduct+shift(conduct,0,1,0))/2. 95 96 u_cond_u=total( vu*conduct_u*e3u3d*umask(),3) 97 98 v_cond_v=total( vv*conduct_v*e3v3d*vmask(),3) 99 100 101 ;calcul de Jp= u_star x Fz k 102 Jpv_star= -1*BR_u*u_cond_u 103 104 Jpu_star= BR_v*v_cond_v 105 106 ; Divergence du champ 107 108 Diver=divfred(Jpu_star,Jpv_star) 109 110 Diver=Diver*1e-6 111 ;stop 112 lecontinent=where(Diver GT 1.E08) 113 Diver(lecontinent)=0. 114 115 116 117 ;stop 118 ;bande de recouvrement: 119 120 diver3(1:180,0:147,*,jt)=Diver(*,*) 121 diver3(0,*,*,*)=diver3(180,*,*,*) 122 diver3(181,*,*,*)=diver3(1,*,*,*) 123 124 sigma3(1:180,0:147,*,jt)=SIGMA(*,*) 125 sigma3(0,*,*,*)=sigma3(180,*,*,*) 126 sigma3(181,*,*,*)=sigma3(1,*,*,*) 127 128 Jpu(1:180,0:147,*,jt)=Jpu_star(*,*) 129 Jpu(0,*,*,*)=Jpu(180,*,*,*) 130 Jpu(181,*,*,*)=Jpu(1,*,*,*) 131 132 Jpv(1:180,0:147,*,jt)=Jpv_star(*,*) 133 Jpv(0,*,*,*)=Jpv(180,*,*,*) 134 Jpv(181,*,*,*)=Jpv(1,*,*,*) 135 136 print, jt 137 138 ENDFOR 139 140 temps=fltarr(73) 141 temps(0)=0. 142 for jt=0,71 do begin &$ 143 temps(jt+1)=temps(jt) +5*86400. &$ 144 endfor 145 print,temps 146 147 vargrid = 'T' 148 iodir = '/usr/work/sur/fvi/OPA/ORCA2/' 149 ; Nom 150 idout = NCDF_CREATE(iodir+'ForcageMOED_5d_'+iyear+'_grid_T.nc',/clobber) 151 print, 'Creation du fichier Netcdf' 152 NCDF_CONTROL, idout, /nofill 153 ; Dimension 154 xidout = NCDF_DIMDEF(idout, 'x',jpiglo) 155 yidout = NCDF_DIMDEF(idout, 'y',jpjglo) 156 didout = NCDF_DIMDEF(idout, 'deptht',1) 157 tidout = NCDF_DIMDEF(idout, 'time_counter', /unlimited) 158 159 didout1 = NCDF_DIMDEF(idout, 'deptht1',jpk) 173 160 174 161 175 162 ; Attributs globaux 176 id0 = NCDF_VARDEF(idout, 'nav_lon' , [xidout, yidout ], /FLOAT)177 id1 = NCDF_VARDEF(idout, 'nav_lat' , [xidout, yidout ], /FLOAT)178 id2 = NCDF_VARDEF(idout, 'deptht' , [ didout1 ], /FLOAT)179 id3 = NCDF_VARDEF(idout, 'time_counter', [ tidout], /FLOAT)180 id4 = NCDF_VARDEF(idout, 'DivJp' , [xidout, yidout, didout, tidout], /DOUBLE)181 id5 = NCDF_VARDEF(idout, 'Sigma' , [xidout, yidout, didout, tidout], /DOUBLE)182 id6 = NCDF_VARDEF(idout, 'Jpu' , [xidout, yidout, didout, tidout], /DOUBLE)183 id7 = NCDF_VARDEF(idout, 'Jpv' , [xidout, yidout, didout, tidout], /DOUBLE)163 id0 = NCDF_VARDEF(idout, 'nav_lon' , [xidout, yidout ], /FLOAT) 164 id1 = NCDF_VARDEF(idout, 'nav_lat' , [xidout, yidout ], /FLOAT) 165 id2 = NCDF_VARDEF(idout, 'deptht' , [ didout1 ], /FLOAT) 166 id3 = NCDF_VARDEF(idout, 'time_counter', [ tidout], /FLOAT) 167 id4 = NCDF_VARDEF(idout, 'DivJp' , [xidout, yidout, didout, tidout], /DOUBLE) 168 id5 = NCDF_VARDEF(idout, 'Sigma' , [xidout, yidout, didout, tidout], /DOUBLE) 169 id6 = NCDF_VARDEF(idout, 'Jpu' , [xidout, yidout, didout, tidout], /DOUBLE) 170 id7 = NCDF_VARDEF(idout, 'Jpv' , [xidout, yidout, didout, tidout], /DOUBLE) 184 171 185 172 ; Variable 0 186 NCDF_ATTPUT, idout, id0, 'units', 'degrees_east'187 NCDF_ATTPUT, idout, id0, 'long_name', 'Longitude'188 NCDF_ATTPUT, idout, id0, 'nav_model', 'Default grid'173 NCDF_ATTPUT, idout, id0, 'units', 'degrees_east' 174 NCDF_ATTPUT, idout, id0, 'long_name', 'Longitude' 175 NCDF_ATTPUT, idout, id0, 'nav_model', 'Default grid' 189 176 ; Variable 1 190 NCDF_ATTPUT, idout, id1, 'units', 'degrees_north'191 NCDF_ATTPUT, idout, id1, 'long_name', 'Latitude'192 NCDF_ATTPUT, idout, id1, 'nav_model', 'Default grid'177 NCDF_ATTPUT, idout, id1, 'units', 'degrees_north' 178 NCDF_ATTPUT, idout, id1, 'long_name', 'Latitude' 179 NCDF_ATTPUT, idout, id1, 'nav_model', 'Default grid' 193 180 ; Variable 2 194 NCDF_ATTPUT, idout, id2, 'units','meters'195 NCDF_ATTPUT, idout, id2, 'long_name','Depth'196 NCDF_ATTPUT, idout, id2, 'nav_model','Default grid'181 NCDF_ATTPUT, idout, id2, 'units','meters' 182 NCDF_ATTPUT, idout, id2, 'long_name','Depth' 183 NCDF_ATTPUT, idout, id2, 'nav_model','Default grid' 197 184 ; Variable3 198 NCDF_ATTPUT, idout, id3, 'units', 'seconds since 0001-01-01 00:00:00 '199 NCDF_ATTPUT, idout, id3, 'calendar','noleap'200 NCDF_ATTPUT, idout, id3, 'title', 'Time'201 NCDF_ATTPUT, idout, id3, 'long_name', 'Time axis'202 NCDF_ATTPUT, idout, id3, 'time_origin','0001-JAN-01 00:00:00'185 NCDF_ATTPUT, idout, id3, 'units', 'seconds since 0001-01-01 00:00:00 ' 186 NCDF_ATTPUT, idout, id3, 'calendar','noleap' 187 NCDF_ATTPUT, idout, id3, 'title', 'Time' 188 NCDF_ATTPUT, idout, id3, 'long_name', 'Time axis' 189 NCDF_ATTPUT, idout, id3, 'time_origin','0001-JAN-01 00:00:00' 203 190 ; Variables 204 NCDF_ATTPUT, idout, id4, 'long_name', 'Divergence Jp'191 NCDF_ATTPUT, idout, id4, 'long_name', 'Divergence Jp' 205 192 ; Variables 206 NCDF_ATTPUT, idout, id5, 'long_name', 'Total Conductance'207 NCDF_ATTPUT, idout, id6, 'long_name', 'Jpu=u_star x B (x comp)'208 NCDF_ATTPUT, idout, id7, 'long_name', 'Jpv=u_star x B (y comp)'209 210 211 NCDF_CONTROL, idout, /ENDEF193 NCDF_ATTPUT, idout, id5, 'long_name', 'Total Conductance' 194 NCDF_ATTPUT, idout, id6, 'long_name', 'Jpu=u_star x B (x comp)' 195 NCDF_ATTPUT, idout, id7, 'long_name', 'Jpv=u_star x B (y comp)' 196 197 198 NCDF_CONTROL, idout, /ENDEF 212 199 213 200 ; Creation de la longitude 214 NCDF_VARPUT, idout, id0, glamt201 NCDF_VARPUT, idout, id0, glamt 215 202 ; Creation de la latitude 216 NCDF_VARPUT, idout, id1, gphit217 ; Creation de la profondeur 218 NCDF_VARPUT, idout, id2, gdept203 NCDF_VARPUT, idout, id1, gphit 204 ; Creation de la profondeur 205 NCDF_VARPUT, idout, id2, gdept 219 206 ; Creation du calendrier 220 207 221 NCDF_VARPUT, idout, id3, temps 222 223 224 ; Ecriture des donnees 225 226 ; ecriture des glam_8 227 NCDF_VARPUT, idout, id4 , diver3 228 NCDF_VARPUT, idout, id5 , sigma3 229 NCDF_VARPUT, idout, id6 , Jpu 230 NCDF_VARPUT, idout, id7 , Jpv 231 232 233 234 NCDF_CLOSE, idout 235 208 NCDF_VARPUT, idout, id3, temps 209 210 211 ; Ecriture des donnees 212 213 ; ecriture des glam_8 214 NCDF_VARPUT, idout, id4 , diver3 215 NCDF_VARPUT, idout, id5 , sigma3 216 NCDF_VARPUT, idout, id6 , Jpu 217 NCDF_VARPUT, idout, id7 , Jpv 218 219 NCDF_CLOSE, idout 236 220 237 221 END -
trunk/forcagequimarche.pro
r36 r48 2 2 ; 3 3 ; @file_comments 4 ; calcul d'une matrice 3D de conductivit ésigma=f(T,S) en fonction de T et S4 ; calcul d'une matrice 3D de conductivité sigma=f(T,S) en fonction de T et S 5 5 ; de ORCA, sur la grille T 6 6 ; … … 67 67 ; 68 68 PRO forcagequimarche, orcares, iyear, ian, DRAKKAR_EXP = drakkar_exp 69 ; 70 ;---- 71 ; check input parameters 72 ;---- 73 ; 74 ; check orcares definition 75 ; 76 CASE orcares OF 77 'ORCA2': BEGIN 78 msg = 'iii : valid orcares parameter = ' + orcares 79 ras = report(msg) 80 filename_oce = 'meshmask_bab.nc' 81 IF keyword_set(DRAKKAR_EXP) THEN BEGIN 82 msg = 'www : unused DRAKKAR_EXP keyword = ' + drakkar_exp 83 ras = report(msg) 84 END 69 ; 70 ;---- 71 ; check input parameters 72 ;---- 73 ; 74 ; check orcares definition 75 ; 76 CASE orcares OF 77 'ORCA2': BEGIN 78 msg = 'iii : valid orcares parameter = ' + orcares 79 ras = report(msg) 80 filename_oce = 'meshmask_bab.nc' 81 IF keyword_set(DRAKKAR_EXP) THEN BEGIN 82 msg = 'www : unused DRAKKAR_EXP keyword = ' + drakkar_exp 83 ras = report(msg) 84 END 85 END 86 ELSE : BEGIN 87 msg = 'eee : invalid orcares parameter = ' + orcares 88 ras = report(msg) 89 msg = 'eee : orcares must be ORCA2 or ORCA025++' 90 ras = report(msg) 91 RETURN 92 END 93 ENDCASE 94 95 ; ++ check iyear 96 ; ++ check ian 97 98 ;++@init2 99 ; init grid 100 initocemeshmask, orcares, DRAKKAR_EXP = drakkar_exp 101 102 @common 103 ; 104 ; check for input files 105 ; 106 ; test if ${GEOMAG_ID} defined 107 geomag_id_env=GETENV('GEOMAG_ID') 108 CASE geomag_id_env OF 109 '' : BEGIN 110 msg = 'eee : ${GEOMAG_ID} is not defined' 111 ras = report(msg) 112 RETURN 85 113 END 86 ELSE : BEGIN 87 msg = 'eee : invalid orcares parameter = ' + orcares 88 ras = report(msg) 89 msg = 'eee : orcares must be ORCA2 or ORCA025++' 90 ras = report(msg) 91 RETURN 92 END 114 ELSE : BEGIN 115 msg = 'iii : ${GEOMAG_ID} is ' + geomag_id_env 116 ras = report(msg) 117 END 118 ENDCASE 119 ; 120 iodirin = isadirectory(geomag_id_env) 121 ; 122 ; existence and protection of ${GEOMAG_ID} 123 IF (FILE_TEST(iodirin, /DIRECTORY, /EXECUTABLE, /READ) EQ 0) THEN BEGIN 124 msg = 'eee : the directory' + iodirin + ' is not accessible.' 125 ras = report(msg) 126 RETURN 127 ENDIF 128 ; 129 ; test if ${GEOMAG_OD} defined 130 geomag_od_env=GETENV('GEOMAG_OD') 131 CASE geomag_od_env OF 132 '' : BEGIN 133 msg = 'eee : ${GEOMAG_OD} is not defined' 134 ras = report(msg) 135 RETURN 136 END 137 ELSE : BEGIN 138 msg = 'iii : ${GEOMAG_OD} is ' + geomag_od_env 139 ras = report(msg) 140 END 93 141 ENDCASE 94 95 ; ++ check iyear 96 ; ++ check ian 97 98 ;++@init2 99 ; init grid 100 initocemeshmask, orcares, DRAKKAR_EXP = drakkar_exp 101 102 @common 103 ; 104 ; check for input files 105 ; 106 ; test if ${GEOMAG_ID} defined 107 geomag_id_env=GETENV('GEOMAG_ID') 108 CASE geomag_id_env OF 109 '' : BEGIN 110 msg = 'eee : ${GEOMAG_ID} is not defined' 111 ras = report(msg) 112 RETURN 113 END 114 ELSE: BEGIN 115 msg = 'iii : ${GEOMAG_ID} is ' + geomag_id_env 116 ras = report(msg) 117 END 118 ENDCASE 119 ; 120 iodirin = isadirectory(geomag_id_env) 121 ; 122 ; existence and protection of ${GEOMAG_ID} 123 IF (FILE_TEST(iodirin, /DIRECTORY, /EXECUTABLE, /READ) EQ 0) THEN BEGIN 124 msg = 'eee : the directory' + iodirin + ' is not accessible.' 125 ras = report(msg) 126 RETURN 127 ENDIF 128 ; 129 ; test if ${GEOMAG_OD} defined 130 geomag_od_env=GETENV('GEOMAG_OD') 131 CASE geomag_od_env OF 132 '' : BEGIN 133 msg = 'eee : ${GEOMAG_OD} is not defined' 134 ras = report(msg) 135 RETURN 136 END 137 ELSE: BEGIN 138 msg = 'iii : ${GEOMAG_OD} is ' + geomag_od_env 139 ras = report(msg) 140 END 141 ENDCASE 142 ; 143 ; check if output data will be possible 144 iodirout = isadirectory(geomag_od_env) 145 ; 146 ; existence and protection 147 IF (FILE_TEST(iodirout, /DIRECTORY, /WRITE) EQ 0) THEN BEGIN 148 msg = 'eee : the directory' + iodirout + ' was not found.' 149 ras = report(msg) 150 RETURN 151 ENDIF 152 ; 153 e_exp='ESS' 154 key_portrait = 0 155 ; stockage des fichiers brut 156 ioDATA=geomag_id_env 157 158 CASE orcares OF 159 'ORCA2': BEGIN 160 file_U=e_exp+'_5d_'+ian+'0101_'+ian+'1231_grid_U.nc' 161 file_V=e_exp+'_5d_'+ian+'0101_'+ian+'1231_grid_V.nc' 162 file_T=e_exp+'_5d_'+ian+'0101_'+ian+'1231_grid_T.nc' 163 jpt = 73 ;time counter des fichiers ci-dessus 164 END 165 ENDCASE 166 file_Sed= geomag_id_env+'cond_sed_'+orcares+'.nc' 167 file_Br= geomag_id_env+'Br_'+orcares + '.nc' 168 ; 169 ; title 170 t_exp= e_exp 171 t_bt = 'bar_transp' 172 ioORLN2 = geomag_id_env 173 ; 174 ;facteur d'echelle vertical for partial steps 175 e3v3d=read_ncdf('e3v_ps',0,/timestep,iodir=ioORLN2,/nostruct,/tout,filename='meshmask_bab.nc') 176 e3u3d=read_ncdf('e3u_ps',0,/timestep,iodir=ioORLN2,/nostruct,/tout,filename='meshmask_bab.nc') 177 SIGMAsed=read_ncdf('cond_sed',0,/timestep,/nostruct,/tout,filename=file_Sed,/cont_nofill) 178 ;BR=read_ncdf('Br',0,/timestep,/nostruct,/tout,filename=file_Br) 179 BR=read_ncdf('Br',0,/timestep,/nostruct,/tout,filename=file_Br,/cont_nofill) 180 ; 181 ; vertical integration: 182 e3t3d=make_array(jpi,jpj,jpk) 183 for k=0, jpk-1 do begin &$ 184 for j=0,jpj-1 do begin &$ 185 for i=0,jpi-1 do begin &$ 186 e3t3d(i,j,k) = e3t(k) &$ 187 endfor &$ 188 endfor &$ 189 endfor 190 ; 191 ;vud = make_array(jpi,jpj,jpt) 192 ;vvd = make_array(jpi, jpj, jpt) 193 divBustar = make_array(jpi, jpj, jpt) 194 diver3=replicate(0.,182,149,1,jpt) 195 sigma3=replicate(0.,182,149,1,jpt) 196 ; 197 FOR jt = 0, jpt-1 DO BEGIN &$ 198 ; 199 ; ouverture des fichiers et stockage en memoire partial steps 200 vu=read_ncdf('vozocrtx',jt,jt, /timestep, iodir=ioDATA,/nostruct,/TOUT,filename=file_U) &$ 142 ; 143 ; check if output data will be possible 144 iodirout = isadirectory(geomag_od_env) 145 ; 146 ; existence and protection 147 IF (FILE_TEST(iodirout, /DIRECTORY, /WRITE) EQ 0) THEN BEGIN 148 msg = 'eee : the directory' + iodirout + ' was not found.' 149 ras = report(msg) 150 RETURN 151 ENDIF 152 ; 153 e_exp='ESS' 154 key_portrait = 0 155 ; stockage des fichiers brut 156 ioDATA=geomag_id_env 157 158 CASE orcares OF 159 'ORCA2': BEGIN 160 file_U=e_exp+'_5d_'+ian+'0101_'+ian+'1231_grid_U.nc' 161 file_V=e_exp+'_5d_'+ian+'0101_'+ian+'1231_grid_V.nc' 162 file_T=e_exp+'_5d_'+ian+'0101_'+ian+'1231_grid_T.nc' 163 jpt = 73 ;time counter des fichiers ci-dessus 164 END 165 ENDCASE 166 file_Sed= geomag_id_env+'cond_sed_'+orcares+'.nc' 167 file_Br= geomag_id_env+'Br_'+orcares + '.nc' 168 ; 169 ; title 170 t_exp= e_exp 171 t_bt = 'bar_transp' 172 ioORLN2 = geomag_id_env 173 ; 174 ;facteur d'echelle vertical for partial steps 175 e3v3d=read_ncdf('e3v_ps',0,/timestep,iodir=ioORLN2,/nostruct,/tout,filename='meshmask_bab.nc') 176 e3u3d=read_ncdf('e3u_ps',0,/timestep,iodir=ioORLN2,/nostruct,/tout,filename='meshmask_bab.nc') 177 SIGMAsed=read_ncdf('cond_sed',0,/timestep,/nostruct,/tout,filename=file_Sed,/cont_nofill) 178 ;BR=read_ncdf('Br',0,/timestep,/nostruct,/tout,filename=file_Br) 179 BR=read_ncdf('Br',0,/timestep,/nostruct,/tout,filename=file_Br,/cont_nofill) 180 ; 181 ; vertical integration: 182 e3t3d=make_array(jpi,jpj,jpk) 183 for k=0, jpk-1 do begin &$ 184 for j=0,jpj-1 do begin &$ 185 for i=0,jpi-1 do begin &$ 186 e3t3d(i,j,k) = e3t(k) &$ 187 endfor &$ 188 endfor &$ 189 endfor 190 ; 191 ;vud = make_array(jpi,jpj,jpt) 192 ;vvd = make_array(jpi, jpj, jpt) 193 divBustar = make_array(jpi, jpj, jpt) 194 diver3=replicate(0.,182,149,1,jpt) 195 sigma3=replicate(0.,182,149,1,jpt) 196 ; 197 FOR jt = 0, jpt-1 DO BEGIN &$ 198 ; 199 ; ouverture des fichiers et stockage en memoire partial steps 200 vu=read_ncdf('vozocrtx',jt,jt, /timestep, iodir=ioDATA,/nostruct,/TOUT,filename=file_U) &$ 201 ;stop 202 vv=read_ncdf('vomecrty',jt,jt, /timestep,iodir=ioDATA,/nostruct,/TOUT,filename=file_V) &$ 203 ;stop 204 ; lecture salinite & temperature 205 temp= read_ncdf('votemper',jt,jt,/timestep,iodir=ioDATA,/nostruct,/tout,filename=file_T) 206 ;stop 207 salin=read_ncdf('vosaline',jt,jt,/timestep,iodir=ioDATA,/nostruct,/tout,filename=file_T) 208 ;stop 209 conduct=0.02047780622061 + 0.00273147624197*temp + 0.00035133182334*temp*temp + 0.09139808809909*salin + 0.00241425798890*salin*temp -0.00023998958774*salin*salin 210 mask_t=where(conduct GT 1.e+19) 211 conduct(mask_t)=0. 212 ; Somme conduct au point T 213 ; 214 ; Calcul SIGMA 215 ; 216 SIGMAoc=total(conduct*e3t3d*tmask,3) 217 SIGMA=SIGMAsed+SIGMAoc 218 ; 219 SIGMA_u=(SIGMA+shift(SIGMA,-1,0))/2. 220 SIGMA_v=(SIGMA+shift(SIGMA,0,1))/2. 221 ; 222 ; Calcul B en points u et v 223 ; 224 BR_u=(BR+shift(BR,-1,0))/2. 225 BR_v=(BR+shift(BR,0,1))/2. 226 ; 227 ; Calcul integrale conduct 228 ; 229 conduct_u=(conduct+shift(conduct,-1,0,0))/2. 230 conduct_v=(conduct+shift(conduct,0,1,0))/2. 231 ; 232 u_cond_u=total( vu*conduct_u*e3u3d*umask(),3) 233 v_cond_v=total( vv*conduct_v*e3v3d*vmask(),3) 234 ; 235 Bu_star= BR_u*u_cond_u/SIGMA_u 236 Bv_star= BR_v*v_cond_v/SIGMA_v 237 238 ; Transport horizontal 239 T_u=total( vu*e3u3d*umask(),3) 240 T_v=total( vv*e3v3d*vmask(),3) 241 242 ; 243 ; Divergence du champ 244 Diver=divfred(Bu_star,Bv_star) 245 246 Diver=Diver*1e-6 247 ;stop 248 lecontinent=where(Diver GT 1.E08) 249 Diver(lecontinent)=0. 250 201 251 ;stop 202 vv=read_ncdf('vomecrty',jt,jt, /timestep,iodir=ioDATA,/nostruct,/TOUT,filename=file_V) &$ 203 ;stop 204 ; lecture salinite & temperature 205 temp= read_ncdf('votemper',jt,jt,/timestep,iodir=ioDATA,/nostruct,/tout,filename=file_T) 206 ;stop 207 salin=read_ncdf('vosaline',jt,jt,/timestep,iodir=ioDATA,/nostruct,/tout,filename=file_T) 208 ;stop 209 conduct=0.02047780622061 + 0.00273147624197*temp + 0.00035133182334*temp*temp + 0.09139808809909*salin + 0.00241425798890*salin*temp -0.00023998958774*salin*salin 210 mask_t=where(conduct GT 1.e+19) 211 conduct(mask_t)=0. 212 ; Somme conduct au point T 213 ; 214 ; Calcul SIGMA 215 ; 216 SIGMAoc=total(conduct*e3t3d*tmask,3) 217 SIGMA=SIGMAsed+SIGMAoc 218 ; 219 SIGMA_u=(SIGMA+shift(SIGMA,-1,0))/2. 220 SIGMA_v=(SIGMA+shift(SIGMA,0,1))/2. 221 ; 222 ; Calcul B en points u et v 223 ; 224 BR_u=(BR+shift(BR,-1,0))/2. 225 BR_v=(BR+shift(BR,0,1))/2. 226 ; 227 ; Calcul integrale conduct 228 ; 229 conduct_u=(conduct+shift(conduct,-1,0,0))/2. 230 conduct_v=(conduct+shift(conduct,0,1,0))/2. 231 ; 232 u_cond_u=total( vu*conduct_u*e3u3d*umask(),3) 233 v_cond_v=total( vv*conduct_v*e3v3d*vmask(),3) 234 ; 235 Bu_star= BR_u*u_cond_u/SIGMA_u 236 Bv_star= BR_v*v_cond_v/SIGMA_v 237 238 ; Transport horizontal 239 T_u=total( vu*e3u3d*umask(),3) 240 T_v=total( vv*e3v3d*vmask(),3) 241 242 ; 243 ; Divergence du champ 244 Diver=divfred(Bu_star,Bv_star) 245 246 Diver=Diver*1e-6 247 ;stop 248 lecontinent=where(Diver GT 1.E08) 249 Diver(lecontinent)=0. 250 251 ;stop 252 ;bande de recouvrement:: 252 ;bande de recouvrement: 253 253 254 254 diver3(1:180,0:147,*,jt)=Diver(*,*) … … 342 342 343 343 NCDF_CLOSE, idout 344 ; ++ l' équivalent avec forcage_output344 ; ++ l'équivalent avec forcage_output 345 345 forcage_output, orcares, variable, 'Divergence', long_name, jpiglo, jpjglo, gphit, glamt,diver3 346 346 -
trunk/fsxspp.pro
r2 r48 10 10 ; 11 11 return, z 12 END 12 END -
trunk/fsyspp.pro
r2 r48 10 10 ; 11 11 return, z 12 END 12 END -
trunk/geomag_profile.sh
r47 r48 30 30 # =========== 31 31 # 32 # .. option:: -d 33 # .. option:: -i 34 # .. option:: -o 35 # .. option:: -t 32 # .. option:: -d <directory> 33 # .. option:: -i <indir> 34 # .. option:: -o <outdir> 35 # .. option:: -t <tempdir> 36 36 # 37 37 # define GEOMAG environment … … 73 73 # ++ besoin de posix 74 74 # 75 # ++ pas de MANPATH defini par d éfaut sur zeus75 # ++ pas de MANPATH defini par défaut sur zeus 76 76 # 77 77 # EVOLUTIONS … … 116 116 usage=" Usage : ${command} -d directory -i indir -o outdir -t tempdir" 117 117 # 118 set +u119 while [ ! -z "${1}" ]120 do121 case ${1} in122 -d)123 # directory for application choosen by user (see svn checkout command used)124 directory=${2}125 shift126 ;;127 -i)128 # directory for inputs choosen by user129 indir=${2}130 shift131 ;;132 -o)133 # directory for outputs choosen by user134 outdir=${2}135 shift136 ;;137 -t)138 # directory for temporary outputs choosen by user139 tempdir=${2}140 shift141 ;;142 *)143 # other choice144 echo "eee : unknown option ${1}"145 echo "${usage}"146 # nb : no exit because this file should be launched by login process147 ;;148 esac149 # next flag150 shift151 done152 unset usage153 #154 118 set -u 155 # 156 # check for ${directory} 157 if [ ! -d ${directory} ] 119 pb=0 120 # 121 minargcount=8 122 maxargcount=8 123 narg=${#} 124 # 125 if [ ${narg} -lt ${minargcount} ] 158 126 then 159 echo " eee : ${directory} not found" 160 # nb : no exit because this file should be launched by login process 127 echo "eee : not enought arguments (${narg} vs ${minargcount})" 128 echo "${usage}" 129 # nb : no exit because this file should be launched by login process 130 pb=1 161 131 fi 162 # 163 # check for permission on directory164 if [ ! -x ${directory} ]132 unset minargcount 133 # 134 if [ ${narg} -gt ${maxargcount} ] 165 135 then 166 echo " eee : ${directory} not reachable" 167 # nb : no exit because this file should be launched by login process 136 echo "eee : too many arguments (${narg} vs ${maxargcount})" 137 echo "${usage}" 138 # nb : no exit because this file should be launched by login process 139 pb=1 168 140 fi 169 # 170 set -u 171 system=$(uname) 172 case "${system}" in 173 IRIX64) 174 echo " www : no specific posix checking" 175 ;; 176 *) 177 set -o posix 178 ;; 179 esac 180 unset system 181 # 182 GEOMAG=${directory} 183 export GEOMAG 184 unset directory 185 # 186 # add GEOMAG tools to PATH 187 # if not already done 188 suppath=$(echo ${GEOMAG} | tr -s "/") 189 echo ${PATH} | grep -q "${suppath}:" 190 test_path=${?} 191 if [ ${test_path} -ne 0 ] 141 unset maxargcount 142 unset narg 143 # 144 if [ ${pb} -eq 0 ] 192 145 then 193 PATH=${suppath}:${PATH} 194 export PATH 195 else 196 # option bavarde oui/non pas encore implantée ++ 197 echo "${command} : iii : ${suppath}/ already in \${PATH}" 146 # 147 while [ ${#} -gt 0 ] 148 do 149 case ${1} in 150 -d) 151 # directory for application choosen by user 152 # (see svn checkout command used) 153 directory=${2} 154 shift 155 ;; 156 -i) 157 # directory for inputs choosen by user 158 indir=${2} 159 shift 160 ;; 161 -o) 162 # directory for outputs choosen by user 163 outdir=${2} 164 shift 165 ;; 166 -t) 167 # directory for temporary outputs choosen by user 168 tempdir=${2} 169 shift 170 ;; 171 *) 172 # other choice 173 echo "eee : unknown option ${1}" 174 echo "${usage}" 175 # nb : no exit because this file should be launched by login process 176 ;; 177 esac 178 # next flag 179 shift 180 done 181 unset usage 198 182 fi 199 unset test_path 200 unset suppath 201 # 202 # add GEOMAG manuals to MANPATH 203 # if not already done 204 suppath=$(echo ${GEOMAG}/doc/manuals/man/ | tr -s "/") 205 echo ${MANPATH} | grep -q "${suppath}:" 206 test_manpath=${?} 207 if [ ${test_manpath} -ne 0 ] 183 if [ ${pb} -eq 0 ] 208 184 then 209 MANPATH=${suppath}:${MANPATH} 210 export MANPATH 211 else 212 # option bavarde oui/non pas encore implantée ++ 213 echo "${command} : iii : ${suppath} already in \${MANPATH}" 185 # 186 # check for ${directory} 187 if [ ! -d ${directory} ] 188 then 189 echo " eee : ${directory} not found" 190 # nb : no exit because this file should be launched by login process 191 fi 192 # 193 # check for permission on directory 194 if [ ! -x ${directory} ] 195 then 196 echo " eee : ${directory} not reachable" 197 # nb : no exit because this file should be launched by login process 198 fi 199 # 200 system=$(uname) 201 case "${system}" in 202 IRIX64) 203 echo " www : no specific posix checking" 204 ;; 205 *) 206 set -o posix 207 ;; 208 esac 209 unset system 210 # 211 GEOMAG=${directory} 212 export GEOMAG 213 unset directory 214 # 215 # add GEOMAG tools to PATH 216 # if not already done 217 suppath=$(echo ${GEOMAG} | tr -s "/") 218 echo ${PATH} | grep -q "${suppath}:" 219 test_path=${?} 220 if [ ${test_path} -ne 0 ] 221 then 222 PATH=${suppath}:${PATH} 223 export PATH 224 else 225 # option bavarde oui/non pas encore implantée ++ 226 echo "${command} : iii : ${suppath}/ already in \${PATH}" 227 fi 228 unset test_path 229 unset suppath 230 # 231 # add GEOMAG manuals to MANPATH 232 # if not already done 233 suppath=$(echo ${GEOMAG}/doc/manuals/man/ | tr -s "/") 234 echo ${MANPATH} | grep -q "${suppath}:" 235 test_manpath=${?} 236 if [ ${test_manpath} -ne 0 ] 237 then 238 MANPATH=${suppath}:${MANPATH} 239 export MANPATH 240 else 241 # option bavarde oui/non pas encore implantée ++ 242 echo "${command} : iii : ${suppath} already in \${MANPATH}" 243 fi 244 unset test_manpath 245 unset suppath 246 # 247 GEOMAG_LOG=${tempdir} 248 export GEOMAG_LOG 249 unset tempdir 250 if [ ! -d ${GEOMAG_LOG} ] 251 then 252 mkdir -p ${GEOMAG_LOG} 253 status=${?} 254 if [ ${status} -ne 0 ] 255 then 256 echo "${command} : eee : can not create \${GEOMAG_LOG}" 257 # nb : no exit because this file should be launched by login process 258 else 259 echo "${command} : iii : creation of \${GEOMAG_LOG}" 260 fi 261 unset status 262 fi 263 # check for permission on GEOMAG_LOG 264 if [ ! -x ${GEOMAG_LOG} ] 265 then 266 echo " eee : ${GEOMAG_LOG} not reachable" 267 # nb : no exit because this file should be launched by login process 268 fi 269 # 270 # check for permission on GEOMAG_LOG 271 if [ ! -w ${GEOMAG_LOG} ] 272 then 273 echo " eee : ${GEOMAG_LOG} not writable" 274 # nb : no exit because this file shouldreachable be launched 275 # by login process 276 fi 277 # 278 EDITOR=vi 279 export EDITOR 280 # 281 # io directories 282 GEOMAG_ID=${indir} 283 export GEOMAG_ID 284 unset indir 285 if [ ! -d ${GEOMAG_ID} ] 286 then 287 mkdir -p ${GEOMAG_ID} 288 echo "${command} : iii : creation of \${GEOMAG_ID}" 289 fi 290 # check for permission on GEOMAG_ID 291 if [ ! -x ${GEOMAG_ID} ] 292 then 293 echo " eee : ${GEOMAG_ID} not reachable" 294 # nb : no exit because this file should be launched by login process 295 fi 296 # 297 GEOMAG_OD=${outdir} 298 export GEOMAG_OD 299 unset outdir 300 if [ ! -d ${GEOMAG_OD} ] 301 then 302 mkdir -p ${GEOMAG_OD} 303 echo "${command} : iii : creation of \${GEOMAG_OD}" 304 fi 305 # check for permission on GEOMAG_OD 306 if [ ! -x ${GEOMAG_OD} ] 307 then 308 echo " eee : ${GEOMAG_OD} not reachable" 309 # nb : no exit because this file should be launched by login process 310 fi 311 if [ ! -w ${GEOMAG_OD} ] 312 then 313 echo " eee : ${GEOMAG_OD} not writable" 314 # nb : no exit because this file should be launched by login process 315 fi 316 # 214 317 fi 215 unset test_manpath 216 unset suppath 217 # 218 GEOMAG_LOG=${tempdir} 219 export GEOMAG_LOG 220 unset tempdir 221 if [ ! -d ${GEOMAG_LOG} ] 222 then 223 mkdir -p ${GEOMAG_LOG} 224 status=${?} 225 if [ ${status} -ne 0 ] 226 then 227 echo "${command} : eee : can not create \${GEOMAG_LOG}" 228 # nb : no exit because this file should be launched by login process 229 else 230 echo "${command} : iii : creation of \${GEOMAG_LOG}" 231 fi 232 unset status 233 fi 234 # check for permission on GEOMAG_LOG 235 if [ ! -x ${GEOMAG_LOG} ] 236 then 237 echo " eee : ${GEOMAG_LOG} not reachable" 238 # nb : no exit because this file should be launched by login process 239 fi 240 # 241 # check for permission on GEOMAG_LOG 242 if [ ! -w ${GEOMAG_LOG} ] 243 then 244 echo " eee : ${GEOMAG_LOG} not writable" 245 # nb : no exit because this file shouldreachable be launched by login process 246 fi 247 # 248 EDITOR=vi 249 export EDITOR 250 # 251 # io directories 252 GEOMAG_ID=${indir} 253 export GEOMAG_ID 254 unset indir 255 if [ ! -d ${GEOMAG_ID} ] 256 then 257 mkdir -p ${GEOMAG_ID} 258 echo "${command} : iii : creation of \${GEOMAG_ID}" 259 fi 260 # check for permission on GEOMAG_ID 261 if [ ! -x ${GEOMAG_ID} ] 262 then 263 echo " eee : ${GEOMAG_ID} not reachable" 264 # nb : no exit because this file should be launched by login process 265 fi 266 # 267 GEOMAG_OD=${outdir} 268 export GEOMAG_OD 269 unset outdir 270 if [ ! -d ${GEOMAG_OD} ] 271 then 272 mkdir -p ${GEOMAG_OD} 273 echo "${command} : iii : creation of \${GEOMAG_OD}" 274 fi 275 # check for permission on GEOMAG_OD 276 if [ ! -x ${GEOMAG_OD} ] 277 then 278 echo " eee : ${GEOMAG_OD} not reachable" 279 # nb : no exit because this file should be launched by login process 280 fi 281 if [ ! -w ${GEOMAG_OD} ] 282 then 283 echo " eee : ${GEOMAG_OD} not writable" 284 # nb : no exit because this file should be launched by login process 285 fi 286 # 318 unset pb 287 319 # end 288 320 unset command -
trunk/getmodelout.sh
r47 r48 27 27 # ======== 28 28 # 29 # To get DRAKKAR G42 experiment output on grid T points bet eween year 9 and 10:29 # To get DRAKKAR G42 experiment output on grid T points between year 9 and 10: 30 30 # 31 31 # .. code-block:: bash … … 51 51 usage=" Usage : ${command} -r orcares -exp drakkar_exp -g grid -o list" 52 52 # 53 set +u54 while [ ! -z "${1}"]53 set -u 54 while [ ${#} -gt 0 ] 55 55 do 56 56 case ${1} in … … 89 89 done 90 90 unset usage 91 #92 set -u93 91 # 94 92 # check GEOMAG environement -
trunk/grad.pro
r41 r48 11 11 res = litchamp(field) 12 12 taille=size(res) 13 grille, mask, glam, gphi, gdep, nx, ny,nz,premierx,premiery,premierz,dernierx, derniery, dernierz 13 grille, mask, glam, gphi, gdep, nx, ny,nz,premierx,premiery,premierz,dernierx, derniery, dernierz 14 14 if n_elements(valmask) EQ 0 then valmask = 1e20 15 15 case strupcase(vargrid) of 16 16 'T':BEGIN 17 17 case direc of 18 'x':BEGIN 18 'x':BEGIN 19 19 divi = e1u[premierx:dernierx, premiery:derniery] 20 20 newmask = (umask())[premierx:dernierx, premiery:derniery, premierz:dernierz] … … 35 35 vargrid = 'W' 36 36 END 37 ELSE :return, report('Bad definition of direction argument')37 ELSE : return, report('Bad definition of direction argument') 38 38 ENDCASE 39 39 END … … 47 47 vargrid = 'T' 48 48 END 49 ELSE :return, report('Bad definition of direction argument')49 ELSE : return, report('Bad definition of direction argument') 50 50 endcase 51 51 END … … 71 71 vargrid = 'W' 72 72 END 73 ELSE :return, report('Bad definition of direction argument')73 ELSE : return, report('Bad definition of direction argument') 74 74 endcase 75 75 END … … 95 95 vargrid = 'W' 96 96 END 97 ELSE :return, report('Bad definition of direction argument')97 ELSE : return, report('Bad definition of direction argument') 98 98 endcase 99 99 END … … 103 103 ; 'y':divi = (shift(e2u, 0, -1))[premierx:dernierx, premiery:derniery] 104 104 ; 'z':divi = e3w[premierz:dernierz] 105 ; ELSE :return, report('Bad definition of direction argument')105 ; ELSE : return, report('Bad definition of direction argument') 106 106 ; endcase 107 107 ; END 108 ELSE :return, report('Bad definition of vargrid')108 ELSE : return, report('Bad definition of vargrid') 109 109 ENDCASE 110 110 res = fitintobox(res) … … 117 117 if earth[0] NE -1 then res[earth] = !values.f_nan 118 118 case direc of 119 'x':BEGIN 119 'x':BEGIN 120 120 res = (shift(res, -1, 0)-res)/divi 121 121 if key_periodique EQ 0 OR nx NE jpi THEN res[nx-1, *] = !values.f_nan 122 122 if vargrid EQ 'T' OR vargrid EQ 'V' then res = shift(res, 1, 0) 123 123 END 124 'y':BEGIN 124 'y':BEGIN 125 125 res = (shift(res, 0, -1)-res)/divi 126 126 res[*, ny-1] = !values.f_nan 127 if vargrid EQ 'T' OR vargrid EQ 'U' then res = shift(res, 0, 1) 128 END 129 ELSE :return, report('Bad definition of direction argument for the type of array')127 if vargrid EQ 'T' OR vargrid EQ 'U' then res = shift(res, 0, 1) 128 END 129 ELSE : return, report('Bad definition of direction argument for the type of array') 130 130 ENDCASE 131 131 earth = where(newmask[*, *, premierz] EQ 0) … … 143 143 divi = divi[*]#replicate(1, jpt) 144 144 case direc of 145 'x':BEGIN 145 'x':BEGIN 146 146 res = (shift(res, -1, 0, 0)-res)/divi 147 147 if key_periodique EQ 0 OR nx NE jpi THEN res[nx-1, *, *] = !values.f_nan 148 148 if vargrid EQ 'T' OR vargrid EQ 'V' then res = shift(res, 1, 0, 0) 149 149 END 150 'y':BEGIN 150 'y':BEGIN 151 151 res = (shift(res, 0, -1, 0)-res)/divi 152 152 res[*, ny-1, *] = !values.f_nan 153 if vargrid EQ 'T' OR vargrid EQ 'U' then res = shift(res, 0, 1, 0) 154 END 155 ELSE :return, report('Bad definition of direction argument for the type of array')153 if vargrid EQ 'T' OR vargrid EQ 'U' then res = shift(res, 0, 1, 0) 154 END 155 ELSE : return, report('Bad definition of direction argument for the type of array') 156 156 ENDCASE 157 157 earth = where(newmask[*, *, premierz] EQ 0) 158 if earth[0] NE -1 then BEGIN 158 if earth[0] NE -1 then BEGIN 159 159 earth = earth#replicate(1, jpt)+replicate(1, n_elements(earth))#(nx*ny*lindgen(jpt)) 160 160 res[earth] = valmask … … 168 168 if earth[0] NE -1 then res[earth] = !values.f_nan 169 169 case direc OF 170 'x':BEGIN 170 'x':BEGIN 171 171 divi = divi[*]#replicate(1, nz) 172 172 res = (shift(res, -1, 0, 0)-res)/divi … … 174 174 if vargrid EQ 'T' OR vargrid EQ 'V' then res = shift(res, 1, 0, 0) 175 175 END 176 'y':BEGIN 176 'y':BEGIN 177 177 divi = divi[*]#replicate(1, nz) 178 178 res = (shift(res, 0, -1, 0)-res)/divi 179 179 res[*, ny-1, *] = !values.f_nan 180 if vargrid EQ 'T' OR vargrid EQ 'U' then res = shift(res, 0, 1, 0) 180 if vargrid EQ 'T' OR vargrid EQ 'U' then res = shift(res, 0, 1, 0) 181 181 END 182 182 'z':BEGIN 183 183 divi = reform(replicate(1, nx*ny)#divi, nx, ny, nz) 184 184 if nx EQ 1 OR ny EQ 1 then res = reform(res, nx, ny, nz) 185 if vargrid EQ 'W' THEN BEGIN 185 if vargrid EQ 'W' THEN BEGIN 186 186 res = (shift(res, 0, 0, 1)-res)/divi 187 187 res[*, *, 0] = !values.f_nan … … 206 206 divi = reform(divi[*]#replicate(1, jpt), nx, ny, nz, jpt, /over) 207 207 if nx EQ 1 OR ny EQ 1 then res = reform(res, nx, ny, nz, jpt) 208 if vargrid EQ 'W' THEN BEGIN 208 if vargrid EQ 'W' THEN BEGIN 209 209 res = (shift(res, 0, 0, 1, 0)-res)/divi 210 210 res[*, *, 0, *] = !values.f_nan … … 227 227 return, res 228 228 end 229 230 231 232 233 234 -
trunk/init2.pro
r2 r48 37 37 homedir = '/usr/work/sur/fvi/OPA/geomag/' 38 38 iodir = '/usr/work/sur/fvi/OPA/ORCA2/' 39 animdir =iodir 39 animdir =iodir 40 40 psdir = '/usr/work/sur/fvi/IDL/IDL_PS/' 41 41 imagedir = iodir -
trunk/init_interp.pro
r2 r48 1 1 ;------------------------------------------------------------ 2 2 ; MODIFICATION HISTORY: Sebastien Masson (smasson@lodyc.jussieu.fr) 3 ; 3 ; 26/5/98 4 4 ;------------------------------------------------------------ 5 5 -
trunk/init_path.pro
r2 r48 9 9 output_dir = '/usr/work/sur/fvi/OPA/geomag/' 10 10 ; 11 ; 2. Input information 11 ; 2. Input information 12 12 ; ==================== 13 13 ; … … 15 15 ; ------------------- 16 16 ; 17 IF keyword_set(nscal) THEN BEGIN 18 data_u_file = 'OMIP_east_west_stress.nc'19 data_v_file = 'OMIP_north_south_stress.nc'20 ENDIF ELSE BEGIN 21 datafile = 'condmag.nc'22 ENDELSE 17 IF keyword_set(nscal) THEN BEGIN 18 data_u_file = 'OMIP_east_west_stress.nc' 19 data_v_file = 'OMIP_north_south_stress.nc' 20 ENDIF ELSE BEGIN 21 datafile = 'condmag.nc' 22 ENDELSE 23 23 ; 24 24 ; 2.2. Data mask file name … … 32 32 gridfile = 'meshmask_bab.nc' 33 33 ; 34 ; 3. Output information 34 ; 3. Output information 35 35 ; ===================== 36 36 ; … … 38 38 ; --------------------- 39 39 ; 40 IF keyword_set(nscal) THEN BEGIN 41 outputfile_x = 'OMIP_east_west_stress_ORCA_R2_grid_U.nc'42 outputfile_y = 'OMIP_north_south_stress_ORCA_R2_grid_V.nc'43 ENDIF ELSE BEGIN 44 ;outputfile = 'cond_sed_ORCA2.nc'45 outputfile = 'Br_ORCA2.nc'46 ENDELSE 40 IF keyword_set(nscal) THEN BEGIN 41 outputfile_x = 'OMIP_east_west_stress_ORCA_R2_grid_U.nc' 42 outputfile_y = 'OMIP_north_south_stress_ORCA_R2_grid_V.nc' 43 ENDIF ELSE BEGIN 44 ; outputfile = 'cond_sed_ORCA2.nc' 45 outputfile = 'Br_ORCA2.nc' 46 ENDELSE 47 47 ; 48 48 printf, 40, '' 49 IF keyword_set(nscal) THEN BEGIN 50 printf, 40, 'Data input file (x-component) : ', string(input_dir, '/', data_u_file)51 printf, 40, 'Data input file (y-component) : ', string(input_dir, '/', data_v_file)52 ENDIF ELSE BEGIN 53 printf, 40, 'Data input file : ', string(input_dir, '/', datafile)54 ENDELSE 49 IF keyword_set(nscal) THEN BEGIN 50 printf, 40, 'Data input file (x-component) : ', string(input_dir, '/', data_u_file) 51 printf, 40, 'Data input file (y-component) : ', string(input_dir, '/', data_v_file) 52 ENDIF ELSE BEGIN 53 printf, 40, 'Data input file : ', string(input_dir, '/', datafile) 54 ENDELSE 55 55 IF keyword_set(key_mask) THEN printf, 40, 'Mask input file : ', string(input_dir, '/', maskfile) 56 IF keyword_set(nscal) THEN BEGIN 57 printf, 40, 'Output file (x-component) : ', string(output_dir, '/', outputfile_x)58 printf, 40, 'Output file (y-component) : ', string(output_dir, '/', outputfile_y)59 ENDIF ELSE BEGIN 60 printf, 40, 'Output file : ', string(output_dir, '/', outputfile)61 ENDELSE 56 IF keyword_set(nscal) THEN BEGIN 57 printf, 40, 'Output file (x-component) : ', string(output_dir, '/', outputfile_x) 58 printf, 40, 'Output file (y-component) : ', string(output_dir, '/', outputfile_y) 59 ENDIF ELSE BEGIN 60 printf, 40, 'Output file : ', string(output_dir, '/', outputfile) 61 ENDELSE 62 62 printf, 40, 'Target grid file : ', string(input_dir, '/', gridfile) 63 -
trunk/initocemeshmask.pro
r41 r48 82 82 'ORCA2': BEGIN 83 83 msg = 'iii : valid orcares parameter = ' + orcares 84 84 ras = report(msg) 85 85 filename_oce = 'meshmask_bab.nc' 86 86 IF keyword_set(DRAKKAR_EXP) THEN BEGIN … … 121 121 ENDELSE 122 122 END 123 ELSE 123 ELSE : BEGIN 124 124 msg = 'eee : invalid orcares parameter = ' + orcares 125 125 ras = report(msg) 126 126 msg = 'eee : orcares must be ORCA2 or ORCA025' 127 127 ras = report(msg) 128 128 RETURN 129 129 END 130 130 ENDCASE 131 131 ; … … 140 140 RETURN 141 141 END 142 ELSE : BEGIN142 ELSE : BEGIN 143 143 msg = 'iii : ${GEOMAG_ID} is ' + geomag_id_env 144 144 ras = report(msg) 145 145 END 146 146 ENDCASE 147 147 iodirin = isadirectory(geomag_id_env) … … 197 197 ; 198 198 END 199 ELSE 199 ELSE : BEGIN 200 200 msg = 'eee : invalid orcares parameter = ' + orcares 201 201 ras = report(msg) … … 203 203 ras = report(msg) 204 204 RETURN 205 205 END 206 206 ENDCASE 207 207 jpt = 1 -
trunk/initorca.pro
r2 r48 1 1 ;------------------------------------------------------------ 2 2 ; MODIFICATION HISTORY: Sebastien Masson (smasson@lodyc.jussieu.fr) 3 ; 3 ; 26/5/98 4 4 ; pris le 25 janvier 99 par Maurice Imbard 5 5 ;------------------------------------------------------------ -
trunk/interp.pro
r41 r48 2 2 ; NAME: interp.pro 3 3 ; 4 ; PURPOSE: Performs a bilinear or bicubic interpolation on 4 ; PURPOSE: Performs a bilinear or bicubic interpolation on 5 5 ; data from a regular origin grid to any grid 6 6 … … 9 9 ; CALLING SEQUENCE: interp, znumdta, zdata_name, mask, t, glam, gphi, zresul 10 10 ; 11 ; INPUTS: 11 ; INPUTS: 12 12 ; 13 13 ; znumdta : ID of input data file 14 14 ; zdata_name : input data name 15 15 ; mask : data mask 16 ; t : timestep 16 ; t : timestep 17 17 ; glam : longitude of target grid 18 ; gphi : latitude " " " 18 ; gphi : latitude " " " 19 19 ; 20 20 ; KEYWORD PARAMETERS: None 21 21 ; 22 ; OUTPUTS: 22 ; OUTPUTS: 23 23 ; zresul : interpolated field 24 24 ; 25 25 ; COMMON BLOCKS: 26 ; 27 ; 28 ; SIDE EFFECTS: 26 ; common_interp.pro 27 ; 28 ; SIDE EFFECTS: 29 29 ; 30 30 ; RESTRICTIONS: … … 42 42 ; grid on meridian axis 43 43 ; 44 ;45 44 ;- 46 45 ;------------------------------------------------------------ … … 55 54 ;------------------------------------------------- 56 55 ; 57 ;58 56 ; 0. Extracting data at timestep t 59 57 ; ================================ … … 62 60 zzresul2 = fltarr(jpioce, jpjoce, jpkoce) 63 61 64 FOR jk = 0, jpkatm-1 DO BEGIN 65 66 62 FOR jk = 0, jpkatm-1 DO BEGIN 63 64 67 65 IF keyword_set(ndim) THEN $ 68 66 ncdf_varget, znumdta, zdata_name, zdata, count = [jpiatm, jpjatm, 1, 1], offset = [0, 0, jk, t] ELSE BEGIN … … 70 68 ncdf_varget, znumdta, zdata_name, zdata, count = [jpiatm, jpjatm, 1], offset = [0, 0, t] ELSE $ 71 69 ncdf_varget, znumdta, zdata_name, zdata, count = [jpiatm, jpjatm,1, 1], offset = [0, 0, 0, t] 72 ENDELSE 73 70 ENDELSE 71 74 72 zmask = fltarr(jpiatm, jpjatm) 75 73 zmask = mask(*, *, jk) 76 74 printf, 40, '' 77 75 printf, 40, 'Data reading OK' 78 IF keyword_set(key_vari) THEN BEGIN 76 IF keyword_set(key_vari) THEN BEGIN 79 77 varid = ncdf_varid(znumdta, zdata_name) 80 78 ncdf_attget, znumdta, varid, "scale_factor", scale_factor 81 79 ncdf_attget, znumdta, varid, "add_offset", add_offset 82 80 83 81 printf, 40, '' 84 82 printf, 40, 'Scale factor and add_offset reading OK' 85 ENDIF 86 ; Scale Factor + Add_offset 83 ENDIF 84 ; Scale Factor + Add_offset 87 85 ; Modif Fred 88 ; 86 ; zdata=ncdf_lec('/usr/work/sur/fvi/OPA/geomag/condmag.nc',var='cond_sed') 89 87 zdata=ncdf_lec('/usr/work/sur/fvi/OPA/geomag/condmag.nc',var='Br') 90 88 … … 99 97 100 98 ; zdata = double(zdata) 101 IF keyword_set(nreverse) THEN BEGIN 99 IF keyword_set(nreverse) THEN BEGIN 102 100 zdata = reverse(zdata, 2) 103 101 printf, 40, 'Re-arranging data in latitudes : North-South --> South-North' … … 109 107 ; ============================================= 110 108 ; 111 112 113 114 IF 115 116 IF 117 109 110 111 112 IF keyword_set(nscal) AND keyword_set(north_pole) THEN BEGIN 113 114 IF zdata_name EQ data_u_name then BEGIN 115 118 116 119 117 zdata_u=zdata … … 126 124 ENDIF 127 125 128 IF 129 126 IF zdata_name EQ data_v_name then BEGIN 127 130 128 131 129 zdata_v=zdata … … 139 137 140 138 ; Check is made here 141 139 142 140 global_vect=replicate(0.,jpiatm,jpjatm,2) 143 141 144 142 global_vect(*,*,0)=zdata_u 145 global_vect(*,*,1)=zdata_v 146 143 global_vect(*,*,1)=zdata_v 144 147 145 northwind, global_vect,zdata_name,t 148 146 149 147 zdata_u=global_vect(*,*,0) 150 148 zdata_v=global_vect(*,*,1) … … 155 153 156 154 157 IF 158 155 IF zdata_name EQ data_u_name then BEGIN 156 159 157 zdata=zdata_u 160 158 161 159 ENDIF 162 160 163 IF 164 161 IF zdata_name EQ data_v_name then BEGIN 162 165 163 zdata=zdata_v 166 164 … … 168 166 169 167 ENDIF 170 171 ; 172 ; 173 ; 1. Applying mask on data and plotting 168 169 ; 170 ; 1. Applying mask on data and plotting 174 171 ; ===================================== 175 172 ; 176 173 val = 1.e20 177 IF keyword_set(nmiss) THEN BEGIN 178 174 IF keyword_set(nmiss) THEN BEGIN 175 179 176 ind = where(abs(zdata) GE abs(missing_val)) 180 177 IF ind[0] NE -1 then zmask(ind) = 0. 181 182 ENDIF 178 179 ENDIF 183 180 184 181 ind = where(zmask EQ 0.) 185 186 IF ind[0] NE -1 then 182 183 IF ind[0] NE -1 then zdata(ind) = val 187 184 188 185 ind = where(zmask NE 0.) 189 IF ind[0] NE -1 THEN min = min(zdata(ind)) ELSE BEGIN 186 IF ind[0] NE -1 THEN min = min(zdata(ind)) ELSE BEGIN 190 187 min = missing_val 191 ENDELSE 192 IF ind[0] NE -1 THEN max = max(zdata(ind)) ELSE BEGIN 188 ENDELSE 189 IF ind[0] NE -1 THEN max = max(zdata(ind)) ELSE BEGIN 193 190 max = missing_val 194 191 zdata(*, *) = missing_val 195 ENDELSE 192 ENDELSE 196 193 IF t MOD ndraw EQ 0 AND keyword_set(keydraw) AND jk MOD nlevel EQ 0 THEN BEGIN 197 194 IF keyword_set(key_ps) THEN openps, string(zdata_name, '_controle_iteration_', strtrim(t+1, 2), '.ps') 198 195 draw, zdata, lonatm, latatm, jpiatm, jpjatm, 1, 1, 3, $ 199 196 TITLE = zdata_name+', iteration = '+strtrim(t+1,1), min, max 200 ENDIF 197 ENDIF 201 198 ; 202 199 ; 2. Extrapolating coastal values on land points … … 207 204 tempsun2 = systime(1) ; pour key_performance 208 205 ; 209 ; ... adding a border to overcome shifting problems 206 ; ... adding a border to overcome shifting problems 210 207 ; 211 208 zero = fltarr(jpiatm+2, jpjatm+2) … … 218 215 219 216 ; 220 ; ... filling 221 ; 222 223 if max(zmask) GT 0 224 225 WHILE n LE nfil AND eps LE 1.e-5 DO BEGIN 217 ; ... filling 218 ; 219 220 if max(zmask) GT 0 then begin 221 222 WHILE n LE nfil AND eps LE 1.e-5 DO BEGIN 226 223 227 224 extrap, zdata, zmask, n, val … … 237 234 zmask = zmask[1:jpiatm, 1:jpjatm] 238 235 ; 239 if keyword_set(key_performance) THEN print, 'temps extrap', systime(1)-tempsun2 236 if keyword_set(key_performance) THEN print, 'temps extrap', systime(1)-tempsun2 240 237 ; .. drawing 241 IF t MOD ndraw EQ 0 AND keyword_set(keydraw) AND jk MOD nlevel EQ 0 THEN BEGIN 238 IF t MOD ndraw EQ 0 AND keyword_set(keydraw) AND jk MOD nlevel EQ 0 THEN BEGIN 242 239 draw, zdata, lonatm, latatm, jpiatm, jpjatm, 1, 2, 3, $ 243 240 TITLE = 'extrapolated '+zdata_name+', iteration = '+strtrim(t+1,1), min, max 244 ENDIF 245 241 ENDIF 242 246 243 ind = where(zmask NE 0. AND zdata GE val) 247 IF n_elements(ind) GT 1 THEN BEGIN 244 IF n_elements(ind) GT 1 THEN BEGIN 248 245 print, 'We stop because of insufficient extrapolation on missing values' 249 246 print, n_elements(ind), ' points non extrapolated' 250 247 print, 'Increase nfill in naminterp.pro' 251 248 stop 252 ENDIF 249 ENDIF 253 250 254 251 255 252 ; 3. Treatment of Boundary Conditions 256 253 ; ======================================================== 257 ; Treatment of boundary conditions require extra stripes 254 ; Treatment of boundary conditions require extra stripes 258 255 ; Size of extra stripe is calculated after a given extension 259 256 ; along x and y axis on the OUTPUT grid. iplus and jplus will be … … 266 263 267 264 ; Initialisation of array for the northern boundary condition 268 265 269 266 zdata_north=replicate(0.,jpiatm,jpjatm+jplus) 270 267 zdata_north(*,0:jpjatm-1)=zdata 271 268 272 269 ; Data of the input grid is duplicated for the northern boundary condition 273 270 274 zdata_inv=shift(zdata_north(*,jpjatm-1-jplus+1:jpjatm-1),jpiatm/2,0) 275 271 zdata_inv=shift(zdata_north(*,jpjatm-1-jplus+1:jpjatm-1),jpiatm/2,0) 272 276 273 zdata_inv=reverse(zdata_inv,2) 277 274 … … 281 278 282 279 zdata=zdata_north 283 280 284 281 ; East/West Boundary Condition, zdata_extend will be the array that 285 282 ; contains both northern boundary condition and E/W boundary … … 288 285 289 286 zdata_extend=replicate(0.,jpiatm+2*iplus,jpjatm+jplus) 290 287 291 288 zdata_extend(iplus:jpiatm-1+iplus,*)=zdata 292 289 293 290 zdata_extend(0:iplus-1,*)=zdata(jpiatm-1-iplus+1:jpiatm-1,*) 294 291 295 292 zdata_extend(jpiatm+iplus:jpiatm+2*iplus-1,*)=zdata(0:iplus-1,*) 296 293 297 294 ; 298 295 ; 3.1 Interpolating data on output grid and plotting result … … 304 301 ; grid 305 302 306 print, 'Time Step=',t 303 print, 'Time Step=',t 307 304 print, 'Input Level=',jk 308 305 … … 311 308 gphi2 = replicate(0.,jpioce,jpjoce) 312 309 313 310 314 311 315 312 for ji = 0, jpioce - 1 do begin 316 for jj = 0, jpjoce - 1 do begin 317 318 319 320 diff = - gphi(ji,jj) $ 321 + phi_field 322 323 diff = where( diff GE 0. ) 324 325 326 j_above=min( diff ) 327 313 for jj = 0, jpjoce - 1 do begin 314 315 diff = - gphi(ji,jj) $ 316 + phi_field 317 318 diff = where( diff GE 0. ) 319 320 321 j_above=min( diff ) 322 328 323 if diff(0) eq -1 then begin 329 324 ; print, 'Output grid goes above input grid latitude top value !!!' … … 334 329 gphi2(ji,jj) = j_fine 335 330 endif else begin 336 337 j_bellow = j_above -1 331 332 j_bellow = j_above -1 338 333 339 334 if j_bellow GE 0 then begin … … 345 340 endif 346 341 347 gphi2(ji,jj) = j_fine 342 gphi2(ji,jj) = j_fine 348 343 if j_fine GT j_above then stop 349 344 350 345 endif else begin 351 346 … … 353 348 354 349 gphi2(ji,jj) = j_fine 355 350 356 351 endelse 357 358 359 endelse 360 352 353 354 endelse 355 361 356 endfor 362 357 endfor 363 358 364 365 IF keyword_set(ninterp) THEN BEGIN 359 360 IF keyword_set(ninterp) THEN BEGIN 366 361 367 362 zresul = interpolate(zdata_extend, glam2, gphi2, cubic = -.5) 368 363 369 ENDIF ELSE BEGIN 370 364 ENDIF ELSE BEGIN 365 371 366 zresul = bilinear(zdata_extend , glam2, gphi2) 372 367 373 368 ENDELSE 374 375 IF t MOD ndraw EQ 0 AND keyword_set(keydraw) AND jk MOD nlevel EQ 0 THEN BEGIN 369 370 IF t MOD ndraw EQ 0 AND keyword_set(keydraw) AND jk MOD nlevel EQ 0 THEN BEGIN 376 371 draw, zresul, glam, gphi, jpioce, jpjoce, 1, 3, 3, $ 377 372 TITLE = 'interpolated '+zdata_name+', iteration = '+strtrim(t+1,1), min, max 378 IF keyword_set(key_ps) THEN BEGIN 373 IF keyword_set(key_ps) THEN BEGIN 379 374 closeps 380 ENDIF ELSE BEGIN 375 ENDIF ELSE BEGIN 381 376 print, 'Press return to continue' 382 377 zxc = strarr(1) 383 378 read, zxc 384 ENDELSE 385 ENDIF 379 ENDELSE 380 ENDIF 386 381 zzresul(*, *, jk) = temporary(zresul) 387 ENDFOR 382 ENDFOR 388 383 zresul = zzresul 389 384 390 385 ; 391 ; 4. Interpolating data on the vertical axis if the outdept not equal datadept 386 ; 4. Interpolating data on the vertical axis if the outdept not equal datadept 392 387 ; =========================================================================== 393 388 394 389 395 IF keyword_set(ndim) THEN BEGIN 390 IF keyword_set(ndim) THEN BEGIN 396 391 397 392 nlev_flag = WHERE ( datadept(0:jpkatm-1) ne outdept(0:jpkoce-1) ) 398 393 IF nlev_flag(0) NE -1 THEN BEGIN 399 FOR jk = 0, jpkoce-1 DO BEGIN 400 394 FOR jk = 0, jpkoce-1 DO BEGIN 395 401 396 jn = where(abs(datadept) GE abs(outdept(jk)) $ 402 397 AND abs(shift(datadept, 1)) LT abs(outdept(jk))) 403 404 405 IF min(jn) GE 0 THEN BEGIN 398 399 400 IF min(jn) GE 0 THEN BEGIN 406 401 407 402 ikt=jn-1 … … 409 404 410 405 IF jn NE 0 then begin 411 406 412 407 zt1=abs(outdept(jk))-abs(datadept(ikt)) 413 408 zt2=abs(datadept(ikb))-abs(outdept(jk)) … … 418 413 419 414 ENDIF ELSE BEGIN 420 415 421 416 zzresul2(*, *, jk) = zzresul(*, *, ikb) 422 423 ENDELSE 424 425 ENDIF ELSE BEGIN 417 418 ENDELSE 419 420 ENDIF ELSE BEGIN 426 421 427 422 zzresul2(*, *, jk) = temporary(zzresul2(*, *, jk-1)) 428 423 429 ENDELSE 430 ENDFOR 424 ENDELSE 425 ENDFOR 431 426 zresul = zzresul2 432 427 ENDIF ELSE BEGIN … … 435 430 print, ' ' 436 431 ENDELSE 437 ENDIF 438 439 440 if keyword_set(key_performance) THEN print, 'temps interp', systime(1)-tempsun 432 ENDIF 433 434 435 if keyword_set(key_performance) THEN print, 'temps interp', systime(1)-tempsun 441 436 442 437 return 443 END 438 END -
trunk/interpolation.pro
r41 r48 9 9 ; CALLING SEQUENCE: interpolation 10 10 ; 11 ; INPUTS: Complete the INIT_PATH and NAMINTERP files 11 ; INPUTS: Complete the INIT_PATH and NAMINTERP files 12 12 ; 13 13 ; KEYWORD PARAMETERS: None … … 16 16 ; 17 17 ; COMMON BLOCKS: 18 ; 19 ; 20 ; SIDE EFFECTS: 18 ; common_interp.pro 19 ; 20 ; SIDE EFFECTS: 21 21 ; 22 22 ; RESTRICTIONS: … … 70 70 71 71 ; 72 ;--------------------------------------------------- 72 ;--------------------------------------------------- 73 73 ; READING INPUT MASK & OUTPUT COORDINATES * 74 74 ;--------------------------------------------------- … … 87 87 lonatm(i) = west_lon + i*dx 88 88 ENDFOR 89 FOR j = 0, jpjatm-1 DO BEGIN 89 FOR j = 0, jpjatm-1 DO BEGIN 90 90 latatm(j) = phi_input(j) 91 ENDFOR 91 ENDFOR 92 92 ; 93 93 ;--------------------------------------------------- … … 99 99 printf, 40, '===========================' 100 100 ; 101 IF keyword_set(negmask) THEN BEGIN 101 IF keyword_set(negmask) THEN BEGIN 102 102 mask = abs(temporary(mask)) 103 103 printf, 40, 'Taking the abs of mask (-1 --> +1)' 104 ENDIF 105 IF keyword_set(nreverse) THEN BEGIN 104 ENDIF 105 IF keyword_set(nreverse) THEN BEGIN 106 106 mask = reverse(mask, 2) 107 107 printf, 40, 'Re-arranging mask in latitudes : North-South --> South-North' 108 ENDIF 108 ENDIF 109 109 ind = where(mask NE 0 AND mask NE 1) 110 IF n_elements(ind) NE 1 THEN BEGIN 110 IF n_elements(ind) NE 1 THEN BEGIN 111 111 printf, 40, '' 112 112 printf, 40, '************************************************************' … … 115 115 printf, 40, '' 116 116 printf, 40, 'program stops' 117 stop 118 ENDIF 119 IF keyword_set(invmask) THEN BEGIN 117 stop 118 ENDIF 119 IF keyword_set(invmask) THEN BEGIN 120 120 mask = inv_mask(mask) 121 121 printf, 40, 'Giving value = 1 to mask on ocean points ' … … 141 141 ; ============================================== 142 142 ; 143 IF keyword_set(nscal) THEN BEGIN 143 IF keyword_set(nscal) THEN BEGIN 144 144 printf, 40, '' 145 145 printf, 40, ' V-1- Writing attributes of interpolated field to files' … … 149 149 printf, 40, '' 150 150 printf, 40, 'Writing attributes of x-component of interpolated field' 151 ncdf_varget, nummsh, mask_u_name, umask, count = [jpioce, jpjoce, jpkoce, 1] 151 ncdf_varget, nummsh, mask_u_name, umask, count = [jpioce, jpjoce, jpkoce, 1] 152 152 netcdf_output, mask_u_name, umask, data_u_name, lon_u, lat_u, outputfile_x, $ 153 153 title_x, long_x_name, POINT = "u", $ … … 156 156 printf, 40, '' 157 157 printf, 40, 'Writing attributes of y-component of interpolated field' 158 ncdf_varget, nummsh, mask_v_name, vmask, count = [jpioce, jpjoce, jpkoce, 1] 158 ncdf_varget, nummsh, mask_v_name, vmask, count = [jpioce, jpjoce, jpkoce, 1] 159 159 netcdf_output, mask_v_name, vmask, data_v_name, lon_v, lat_v, outputfile_y, $ 160 160 title_y, long_y_name, POINT = "v", $ 161 numout_v, varid_v, numdta_v 162 ENDIF ELSE BEGIN 161 numout_v, varid_v, numdta_v 162 ENDIF ELSE BEGIN 163 163 printf, 40, '' 164 164 printf, 40, ' V-1- Writing attributes of interpolated field to file' 165 165 printf, 40, ' -----------------------------------------------------' 166 166 ; 167 ncdf_varget, nummsh, mask_t_name, tmask, count = [jpioce, jpjoce, jpkoce, 1] 167 ncdf_varget, nummsh, mask_t_name, tmask, count = [jpioce, jpjoce, jpkoce, 1] 168 168 netcdf_output, mask_t_name, tmask, data_name, lon_t, lat_t, outputfile, $ 169 169 title, long_name, POINT = "t", $ 170 170 numout, varid, numdta 171 ENDELSE 172 173 174 ; 175 ; 1. Interpolation 171 ENDELSE 172 173 174 ; 175 ; 1. Interpolation 176 176 ; ================ 177 177 ; … … 184 184 ; BEGINNING OF TIME LOOP 185 185 ; 186 IF keyword_set(nscal) THEN BEGIN 186 IF keyword_set(nscal) THEN BEGIN 187 187 mi_u = 1.e20 188 188 ma_u = -1.e20 189 189 mi_v = 1.e20 190 190 ma_v = -1.e20 191 ENDIF ELSE BEGIN 191 ENDIF ELSE BEGIN 192 192 mi = 1.e20 193 193 ma = -1.e20 … … 195 195 196 196 197 FOR t = nit000-1, nitend-1 DO BEGIN 197 FOR t = nit000-1, nitend-1 DO BEGIN 198 198 tempsun4 = systime(1) ; pour key_performance 199 199 printf, 40, '' … … 204 204 IF keyword_set(nscal) THEN BEGIN 205 205 printf, 40, 'A vectorial field has been chosen for interpolation' 206 IF keyword_set(north_pole) THEN BEGIN 206 IF keyword_set(north_pole) THEN BEGIN 207 207 printf, 40, 'Vectorial Fied will be Tested for Coherence nearby North Pole' 208 208 ENDIF … … 218 218 printf, 40, 'interpolating ', data_name, ' at t-point' 219 219 interp, numdta, data_name, 1,mask, t, lon_t, lat_t, phi_input, zresul 220 ENDELSE 220 ENDELSE 221 221 ENDELSE 222 222 … … 224 224 225 225 ; 226 ; 2. Angles correction for 2D fields 226 ; 2. Angles correction for 2D fields 227 227 ; ================================== 228 ; 229 IF keyword_set(nscal) THEN BEGIN 228 ; 229 IF keyword_set(nscal) THEN BEGIN 230 230 printf, 40, '' 231 231 printf, 40, ' angles correction' … … 233 233 tempsun2 = systime(1) ; pour key_performance 234 234 correc_angle, zresul_uu, zresul_uv, zresul_vv, zresul_vu, zresul_u, zresul_v 235 if keyword_set(key_performance) THEN print, 'temps angles', systime(1)-tempsun2 235 if keyword_set(key_performance) THEN print, 'temps angles', systime(1)-tempsun2 236 236 IF t MOD ndraw EQ 0 AND keyword_set(keydraw) THEN BEGIN 237 237 ncdf_varget, numdta_u, data_u_name, zdata_u, $ … … 245 245 zu = zdata_u 246 246 zv = zdata_v 247 IF keyword_set(nmiss) THEN BEGIN 247 IF keyword_set(nmiss) THEN BEGIN 248 248 ind = where(abs(zu) GE missing_val) 249 249 zu(ind) = 0. 250 250 ind = where(abs(zv) GE missing_val) 251 251 zv(ind) = 0. 252 ENDIF 252 ENDIF 253 253 ind = where(mask NE 0.) 254 254 minu = min(zu(ind)) … … 266 266 draw, zresul_v, lon_f, lat_f, jpioce, jpjoce, 2, 4, 2, $ 267 267 TITLE = data_v_name+' interpolated, iteration = '+strtrim(t+1,1), minv, maxv, /land 268 IF keyword_set(key_ps) THEN BEGIN 268 IF keyword_set(key_ps) THEN BEGIN 269 269 closeps 270 ENDIF ELSE BEGIN 270 ENDIF ELSE BEGIN 271 271 print, 'Press return to continue' 272 272 zxc = strarr(1) 273 273 read, zxc 274 ENDELSE 274 ENDELSE 275 275 ENDIF 276 ENDIF 277 278 ; 279 ; 3. NetCDF output to file of interpolated field 276 ENDIF 277 278 ; 279 ; 3. NetCDF output to file of interpolated field 280 280 ; ============================================== 281 281 ; 282 282 tempsun3 = systime(1) ; pour key_performance 283 IF keyword_set(nscal) THEN BEGIN 283 IF keyword_set(nscal) THEN BEGIN 284 284 printf, 40, '' 285 285 printf, 40, 'Storing components of interpolated field' 286 286 287 287 288 288 ncdf_varput, numout_u, varid_u, zresul_u, offset = [0, 0, 0, t-nit000+1] 289 289 ind = where(abs(zresul_u) LT 1.e10) 290 IF ind[0] NE -1 THEN BEGIN 290 IF ind[0] NE -1 THEN BEGIN 291 291 IF min(zresul_u) LT mi_u THEN mi_u = min(zresul_u[ind]) 292 292 IF max(zresul_u) GT ma_u THEN ma_u = max(zresul_u[ind]) 293 ENDIF ELSE BEGIN 293 ENDIF ELSE BEGIN 294 294 mi_u = missing_val 295 295 ma_u = missing_val 296 296 ENDELSE 297 297 298 298 ncdf_varput, numout_v, varid_v, zresul_v, offset = [0, 0, 0, t-nit000+1] 299 299 ind = where(abs(zresul_v) LT 1.e10) 300 IF ind[0] NE -1 THEN BEGIN 300 IF ind[0] NE -1 THEN BEGIN 301 301 IF min(zresul_v) LT mi_v THEN mi_v = min(zresul_v[ind]) 302 302 IF max(zresul_v) GT ma_v THEN ma_v = max(zresul_v[ind]) 303 ENDIF ELSE BEGIN 303 ENDIF ELSE BEGIN 304 304 mi_v = missing_val 305 305 ma_v = missing_val 306 ENDELSE 307 ENDIF ELSE BEGIN 306 ENDELSE 307 ENDIF ELSE BEGIN 308 308 printf, 40, '' 309 309 printf, 40, 'Storing field' 310 310 311 311 ncdf_varput, numout, varid, zresul, offset = [0, 0, 0, t-nit000+1] 312 312 ind = where(abs(zresul) LT 1.e10) 313 313 IF min(zresul) LT mi THEN mi = min(zresul[ind]) 314 314 IF max(zresul) GT ma THEN ma = max(zresul[ind]) 315 ENDELSE 316 if keyword_set(key_performance) THEN print, 'temps ecriture variable', systime(1)-tempsun3 317 if keyword_set(key_performance) THEN printf, 40, 'temps ecriture variable', systime(1)-tempsun3 318 if keyword_set(key_performance) THEN print, 'temps boucle', systime(1)-tempsun4 319 if keyword_set(key_performance) THEN printf, 40, 'temps boucle', systime(1)-tempsun4 315 ENDELSE 316 if keyword_set(key_performance) THEN print, 'temps ecriture variable', systime(1)-tempsun3 317 if keyword_set(key_performance) THEN printf, 40, 'temps ecriture variable', systime(1)-tempsun3 318 if keyword_set(key_performance) THEN print, 'temps boucle', systime(1)-tempsun4 319 if keyword_set(key_performance) THEN printf, 40, 'temps boucle', systime(1)-tempsun4 320 320 ENDFOR 321 321 ; … … 330 330 ; storing min and max in attributes 331 331 ; 332 IF keyword_set(nscal) THEN BEGIN 332 IF keyword_set(nscal) THEN BEGIN 333 333 ncdf_control, numout_u, /redef 334 334 ncdf_control, numout_v, /redef … … 337 337 ncdf_attput, numout_v, varid_v, "valid_min", mi_v 338 338 ncdf_attput, numout_v, varid_v, "valid_max", ma_v 339 ENDIF ELSE BEGIN 339 ENDIF ELSE BEGIN 340 340 ncdf_control, numout, /redef 341 341 ncdf_attput, numout, varid, "valid_min", mi 342 342 ncdf_attput, numout, varid, "valid_max", ma 343 ENDELSE 343 ENDELSE 344 344 IF keyword_set(key_mask) THEN ncdf_close, nummsk 345 IF keyword_set(nscal) THEN BEGIN 345 IF keyword_set(nscal) THEN BEGIN 346 346 ncdf_close, numdta_u 347 347 ncdf_close, numdta_v 348 348 ncdf_close, numout_u 349 349 ncdf_close, numout_v 350 ENDIF ELSE BEGIN 350 ENDIF ELSE BEGIN 351 351 ncdf_close, numout 352 352 ncdf_close, numdta 353 ENDELSE 353 ENDELSE 354 354 ncdf_close, nummsh 355 355 ; … … 360 360 ; 361 361 close, 40 362 if keyword_set(key_performance) THEN print, 'temps total', systime(1)-tempsun 363 if keyword_set(key_performance) THEN printf,40, 'temps total', systime(1)-tempsun 362 if keyword_set(key_performance) THEN print, 'temps total', systime(1)-tempsun 363 if keyword_set(key_performance) THEN printf,40, 'temps total', systime(1)-tempsun 364 364 ; 365 365 return 366 END 366 END -
trunk/inv_mask.pro
r41 r48 8 8 ; CALLING SEQUENCE: inv_mask, mask 9 9 ; 10 ; INPUTS: 11 ; mask : mask field composed of 0 and 1 values 10 ; INPUTS: 11 ; mask : mask field composed of 0 and 1 values 12 12 ; 13 13 ; KEYWORD PARAMETERS: None 14 14 ; 15 ; OUTPUTS: 16 ; mask : mask field composed of 0 and 1 values 17 ; inverted 15 ; OUTPUTS: 16 ; mask : mask field composed of 0 and 1 values 17 ; inverted 18 18 ; 19 19 ; COMMON BLOCKS: None 20 20 ; 21 ; SIDE EFFECTS: 21 ; SIDE EFFECTS: 22 22 ; 23 23 ; RESTRICTIONS: … … 43 43 printf, 40, 'Inversion done' 44 44 return, z 45 END 45 END -
trunk/naminterp.pro
r2 r48 13 13 ndim = 0 14 14 ; Case ndim = 0 15 ; -if the field and mask are really 2D (x,y(,t)) then 15 ; -if the field and mask are really 2D (x,y(,t)) then 16 16 ; key_2d = 1 17 17 ; -else if they are (x,y,z_a(,t)) where z_a=1 then … … 35 35 ; 1.4. Mask on points 36 36 ; ------------------- 37 ; key_mask = 1 if some points are masked 37 ; key_mask = 1 if some points are masked 38 38 ; = 0 if no masked points 39 39 ; … … 48 48 ; ------------------------- 49 49 ; 50 keydraw = 0 ; =0 if no drawing =1 else 50 keydraw = 0 ; =0 if no drawing =1 else 51 51 key_ps = 0 ; =0 if no Postscript file wanted =1 if PS wanted (implies keydraw=1) 52 52 ndraw = 1 53 53 nlevel = 33 ; every nlevel levels 54 54 ; 55 ; 56 ; 2. Information about the data file 55 ; 2. Information about the data file 57 56 ; ================================== 58 57 ; … … 60 59 ; ----------------- 61 60 ; ... data 62 IF keyword_set(nscal) THEN BEGIN 61 IF keyword_set(nscal) THEN BEGIN 63 62 data_u_name = 'stress_x' 64 63 data_v_name = 'stress_y' … … 73 72 IF keyword_set(key_mask) THEN mask_name = 'tmask' 74 73 ; Input Latitude Axis 75 IF keyword_set(nscal) THEN BEGIN 74 IF keyword_set(nscal) THEN BEGIN 76 75 y_axis_u = 'nav_lat' 77 76 y_axis_v = 'nav_lat' 78 ENDIF ELSE BEGIN 77 ENDIF ELSE BEGIN 79 78 y_axis = 'la' 80 79 ENDELSE 81 80 82 81 ; ... input grid depth (3D field only) 83 IF keyword_set(ndim) THEN BEGIN 82 IF keyword_set(ndim) THEN BEGIN 84 83 datadept_name = 'deptht' 85 84 ; ... location (ndep =1 if depth in data file, =0 if in mask file) 86 85 ndep = 1 87 86 ENDIF 88 ; ... time 87 ; ... time 89 88 time_name = 'lo' 90 89 ; ... output grid mask(s) 91 IF keyword_set(nscal) THEN BEGIN 90 IF keyword_set(nscal) THEN BEGIN 92 91 mask_u_name = 'umask' 93 92 mask_v_name = 'vmask' 94 ENDIF ELSE BEGIN 93 ENDIF ELSE BEGIN 95 94 mask_t_name = 'tmask' 96 95 ENDELSE 97 ; ... output grid depth 96 ; ... output grid depth 98 97 IF keyword_set(ndim) THEN outdept_name = 'gdept' 99 ; ... output grid coordinates 100 IF keyword_set(nscal) THEN BEGIN 98 ; ... output grid coordinates 99 IF keyword_set(nscal) THEN BEGIN 101 100 lon_u_name = 'glamu' 102 101 lat_u_name = 'gphiu' … … 105 104 lon_f_name = 'glamf' 106 105 lat_f_name = 'gphif' 107 ENDIF ELSE BEGIN 106 ENDIF ELSE BEGIN 108 107 lon_name = 'nav_lon' 109 108 lat_name = 'nav_lat' 110 ENDELSE 109 ENDELSE 111 110 ; ... missing value 112 ; 111 ; 113 112 ; nmiss = 0 --> no missing value 114 113 ; nmiss = 1 --> missing value of value missing_val … … 141 140 nit000 = 1 ; first time step 142 141 nitend = 1 ; last time step 143 nstep = nitend-nit000+1 ; number of steps 142 nstep = nitend-nit000+1 ; number of steps 144 143 ; 145 144 ; 2.4. Pretreatment of input variables 146 145 ; ------------------------------------ 147 146 ; data = data*scale_factor + add_offset 148 ; when handling multiple files with scale_factor and 147 ; when handling multiple files with scale_factor and 149 148 ; add_offset defined, put key_vari=1 in case they change 150 149 ; from one file to another … … 165 164 ; 166 165 ; mask = 1 on OCEAN points --> invmask = 0 167 ; mask = 1 on LAND points --> invmask = 1 166 ; mask = 1 on LAND points --> invmask = 1 168 167 ; 169 168 invmask = 0 … … 174 173 ; 175 174 extend = 0 176 ;177 175 ; 178 176 ; 3. Information about the output grid … … 191 189 pres_y=2. ; Precision of the output grid along y axis 192 190 ext_x=40 ; Number of extension point on the output grid to treat for E/W boundary conditions 193 ext_y=40 ; Number of extension point on the output grid to treat for North fold boundary conditions; 191 ext_y=40 ; Number of extension point on the output grid to treat for North fold boundary conditions; 194 192 ; modlam --> 0 if all longitudes are lower than 360 degrees 195 193 ; --> 1 if some longitudes are greater than 360 degrees 196 ; 194 ; 197 195 modlam = 0 198 196 ; 199 ; 200 ; 4. Title(s) of output file(s) 197 ; 4. Title(s) of output file(s) 201 198 ; ============================= 202 199 ; 203 IF keyword_set(nscal) THEN BEGIN 200 IF keyword_set(nscal) THEN BEGIN 204 201 long_x_name = 'Zonal wind stress' 205 202 title_x = "Daily interpolated OMIP stress x wind on ORCA2 grid" 206 203 long_y_name = 'Meridional wind stress' 207 204 title_y = "Daily interpolated OMIP stress y wind on ORCA2 grid" 208 ENDIF ELSE BEGIN 205 ENDIF ELSE BEGIN 209 206 long_name = 'Salinity' 210 207 title = "Monthly Levitus Sea Salinity corrected" 211 ENDELSE 208 ENDELSE 212 209 ; 213 210 key_performance = 0 … … 217 214 IF keyword_set(ndim) THEN printf, 40, ' --> 3D' ELSE printf, 40, ' --> 2D' 218 215 ; 219 IF keyword_set(nscal) THEN BEGIN220 printf, 40, 'Interpolation of the vector array : ' 221 printf, 40, ' --> ', '[', data_u_name, ',', data_v_name, ']' 222 ENDIF ELSE BEGIN 223 printf, 40, 'Interpolation of the scalar array : ' 216 IF keyword_set(nscal) THEN BEGIN 217 printf, 40, 'Interpolation of the vector array : ' 218 printf, 40, ' --> ', '[', data_u_name, ',', data_v_name, ']' 219 ENDIF ELSE BEGIN 220 printf, 40, 'Interpolation of the scalar array : ' 224 221 printf, 40, ' --> ', data_name 225 ENDELSE 222 ENDELSE 226 223 ; 227 224 printf, 40, 'Type of the interpolation :' -
trunk/netcdf_input.pro
r41 r48 2 2 ; NAME: netcdf_input.pro 3 3 ; 4 ; PURPOSE: Reads coordinates, data and mask 4 ; PURPOSE: Reads coordinates, data and mask 5 5 ; 6 6 ; CATEGORY: Subroutine … … 8 8 ; CALLING SEQUENCE: netcdf_input, mask 9 9 ; 10 ; INPUTS: 10 ; INPUTS: 11 11 ; 12 12 ; maskfile : name of mask file 13 13 ; gridfile : name of grid file 14 ; mask_name : mask name 14 ; mask_name : mask name 15 15 ; SCALAR FIELD : 16 16 ; datafile : name of input data file … … 29 29 ; KEYWORD PARAMETERS: None 30 30 ; 31 ; OUTPUTS: 31 ; OUTPUTS: 32 32 ; mask : mask of the origin grid 33 33 ; nummsk : ID of input mask file 34 34 ; nummsh : ID of output coordinates and mask file 35 ; SCALAR FIELD : 35 ; SCALAR FIELD : 36 36 ; lon_t : longitudes of target grid at t-point 37 ; lat_t : latitudes of target grid at t-point 37 ; lat_t : latitudes of target grid at t-point 38 38 ; numdta : ID of input data file 39 39 ; VECTOR FIELD : 40 40 ; lon_u : longitudes of target grid at u-point 41 ; lat_u : latitudes of target grid at u-point 41 ; lat_u : latitudes of target grid at u-point 42 42 ; lon_v : longitudes of target grid at v-point 43 ; lat_v : latitudes of target grid at v-point 43 ; lat_v : latitudes of target grid at v-point 44 44 ; lon_f : longitudes of target grid at f-point 45 ; lat_f : latitudes of target grid at f-point 45 ; lat_f : latitudes of target grid at f-point 46 46 ; numdta_u : ID of input data x-component file 47 47 ; numdta_v : ID of input data y-component file 48 48 ; 49 49 ; COMMON BLOCKS: 50 ; 50 ; common_interp.pro 51 51 ; 52 ; SIDE EFFECTS: 52 ; SIDE EFFECTS: 53 53 ; 54 54 ; RESTRICTIONS: … … 74 74 IF keyword_set(key_mask) THEN nummsk = ncdf_open(input_dir+ '/'+ maskfile) 75 75 IF keyword_set(key_mask) THEN printf, 40, maskfile, ' opening OK' 76 IF keyword_set(nscal) THEN BEGIN 76 IF keyword_set(nscal) THEN BEGIN 77 77 numdta_u = ncdf_open(string(input_dir, '/', data_u_file)) 78 78 printf, 40, data_u_file, ' opening OK' 79 79 numdta_v = ncdf_open(string(input_dir, '/', data_v_file)) 80 80 IF data_u_file NE data_v_file THEN printf, 40, data_v_file, ' opening OK' 81 ENDIF ELSE BEGIN 81 ENDIF ELSE BEGIN 82 82 numdta = ncdf_open(input_dir+ '/'+ datafile) 83 83 printf, 40, datafile, ' opening OK' 84 ENDELSE 84 ENDELSE 85 85 nummsh = ncdf_open(input_dir+ '/'+ gridfile) 86 86 printf, 40, gridfile, ' opening OK' … … 89 89 ; ===================================== 90 90 ; 91 IF keyword_set(nscal) THEN BEGIN 91 IF keyword_set(nscal) THEN BEGIN 92 92 ncdf_varget, nummsh, lon_u_name, lon_u 93 93 ncdf_varget, nummsh, lat_u_name, lat_u … … 96 96 ncdf_varget, nummsh, lon_v_name, lon_v 97 97 ncdf_varget, nummsh, lat_v_name, lat_v 98 IF keyword_set(modlam) THEN BEGIN 98 IF keyword_set(modlam) THEN BEGIN 99 99 ii=where(lon_f gt 360) 100 100 lon_f(ii)=lon_f(ii)-360 … … 103 103 ii=where(lon_v gt 360) 104 104 lon_v(ii)=lon_v(ii)-360 105 ENDIF 106 ENDIF ELSE BEGIN 105 ENDIF 106 ENDIF ELSE BEGIN 107 107 ncdf_varget, nummsh, lon_name, lon_t 108 108 ncdf_varget, nummsh, lat_name, lat_t 109 IF keyword_set(modlam) THEN BEGIN 109 IF keyword_set(modlam) THEN BEGIN 110 110 ii=where(lon_t gt 360) 111 111 lon_t(ii)=lon_t(ii)-360 112 ENDIF 113 ENDELSE 112 ENDIF 113 ENDELSE 114 114 printf, 40, '' 115 115 printf, 40, 'Coordinates of output grid reading OK' … … 122 122 123 123 IF keyword_set(nscal) THEN BEGIN 124 124 125 125 ncdf_varget, numdta_u,y_axis_u,phi_input_u, count = [jpiatm, jpjatm] 126 126 ncdf_varget, numdta_v,y_axis_v,phi_input_v, count = [jpiatm, jpjatm] 127 127 128 128 phi_input_u=phi_input_u(0,*) 129 129 phi_input_v=phi_input_v(0,*) 130 130 131 131 ENDIF ELSE BEGIN 132 132 133 134 ncdf_varget, numdta ,y_axis ,phi_input , count = [ jpjatm] 135 ; modif FRED133 134 ncdf_varget, numdta ,y_axis ,phi_input , count = [ jpjatm] 135 ; modif FRED 136 136 phi_input=reverse(phi_input) 137 ; fin modif FRED 137 ; fin modif FRED 138 138 ;stop 139 139 ; phi_input=phi_input(0,*) … … 143 143 ; Data Mask 144 144 145 IF keyword_set(key_mask) THEN BEGIN 145 IF keyword_set(key_mask) THEN BEGIN 146 146 IF keyword_set(ndim) THEN $ 147 ncdf_varget, nummsk, mask_name, mask, count = [jpiatm, jpjatm, jpkatm, 1] ELSE BEGIN 147 ncdf_varget, nummsk, mask_name, mask, count = [jpiatm, jpjatm, jpkatm, 1] ELSE BEGIN 148 148 149 149 IF keyword_set(key_2d) THEN $ 150 150 ncdf_varget, nummsk, mask_name, mask, count = [jpiatm, jpjatm] ELSE $ 151 151 ncdf_varget, nummsk, mask_name, mask, count = [jpiatm, jpjatm, 1, 1] 152 ENDELSE 153 152 ENDELSE 153 154 154 printf, 40, '' 155 155 printf, 40, 'Data mask reading OK' … … 157 157 IF keyword_set(ndim) THEN mask = replicate(1., jpiatm, jpjatm, jpkatm) ELSE $ 158 158 mask = replicate(1., jpiatm, jpjatm) 159 ENDELSE 159 ENDELSE 160 160 ; 161 161 ; 3. Reading depths of input and output grids (3D field) 162 162 ; ====================================================== 163 163 ; 164 IF keyword_set(ndim) THEN BEGIN 164 IF keyword_set(ndim) THEN BEGIN 165 165 IF keyword_set(ndep) THEN $ 166 166 ncdf_varget, numdta, datadept_name, datadept ELSE $ … … 169 169 datadept = float(datadept) 170 170 outdept = float(outdept) 171 ENDIF 171 ENDIF 172 172 ; 173 173 return 174 END 174 END -
trunk/netcdf_output.pro
r41 r48 9 9 ; zoutputfile, ztitle, POINT = point, znumout, zvarid 10 10 ; 11 ; INPUTS: 11 ; INPUTS: 12 12 ; zout_mask_name : name of interpolated field mask 13 ; 13 ; zout_mask : mask of interpolated field 14 14 ; zdata_name : name of data 15 15 ; zlon,zlat : coordinates of output grid … … 20 20 ; KEYWORD PARAMETERS: None 21 21 ; 22 ; OUTPUTS: 22 ; OUTPUTS: 23 23 ; NetCDF file containing the attributes of interpolated field 24 24 ; znumout : ID of output file … … 26 26 ; 27 27 ; COMMON BLOCKS: 28 ; 28 ; common_interp.pro 29 29 ; 30 ; SIDE EFFECTS: 30 ; SIDE EFFECTS: 31 31 ; 32 32 ; RESTRICTIONS: … … 66 66 ; -------------- 67 67 ; 68 IF size(zlon, /n_dimensions) EQ 1 THEN BEGIN 68 IF size(zlon, /n_dimensions) EQ 1 THEN BEGIN 69 69 varlon = ncdf_vardef(znumout, 'nav_lon', [xid], /float) 70 70 varlat = ncdf_vardef(znumout, 'nav_lat', [yid], /float) 71 ENDIF ELSE BEGIN 71 ENDIF ELSE BEGIN 72 72 varlon = ncdf_vardef(znumout, 'nav_lon', [xid, yid], /float) 73 73 varlat = ncdf_vardef(znumout, 'nav_lat', [xid, yid], /float) 74 ENDELSE 75 IF keyword_set(ndim) THEN BEGIN 74 ENDELSE 75 IF keyword_set(ndim) THEN BEGIN 76 76 varlev1 = ncdf_varid(nummsh, 'nav_lev') 77 77 vart1 = ncdf_varid(znumdta, time_name) 78 varlev = ncdf_vardef(znumout, 'nav_lev', [zid], /float) 78 varlev = ncdf_vardef(znumout, 'nav_lev', [zid], /float) 79 79 vart = ncdf_vardef(znumout, 'time', [tid], /float) 80 80 ; varmsk = ncdf_vardef(znumout, zout_mask_name, [xid, yid, zid], /float) 81 81 zvarid = ncdf_vardef(znumout, zdata_name, [xid, yid, zid, tid], /float) 82 ENDIF ELSE BEGIN 82 ENDIF ELSE BEGIN 83 83 vart1 = ncdf_varid(znumdta, time_name) 84 84 vart = ncdf_vardef(znumout, 'time', [tid], /float) 85 85 ; varmsk = ncdf_vardef(znumout, zout_mask_name, [xid, yid], /float) 86 86 zvarid = ncdf_vardef(znumout, zdata_name, [xid, yid, tid], /float) 87 ENDELSE 87 ENDELSE 88 88 ; 89 89 ; 2.3. Attributes … … 98 98 ncdf_attput, znumout, varlat, "valid_max", max(zlat) 99 99 ncdf_attput, znumout, varlat, "long_name", string("Latitude at ", point, "-point") 100 IF keyword_set(ndim) THEN BEGIN 100 IF keyword_set(ndim) THEN BEGIN 101 101 rien = ncdf_attcopy(nummsh, varlev1, "units", znumout, varlev) 102 102 rien = ncdf_attcopy(nummsh, varlev1, "valid_min", znumout, varlev) … … 106 106 zatt = ncdf_varinq(znumdta, vart1) 107 107 znbatt = zatt.natts 108 FOR i = 0, znbatt-1 DO BEGIN 108 FOR i = 0, znbatt-1 DO BEGIN 109 109 namatt = ncdf_attname(znumdta, vart1, i) 110 110 rien = ncdf_attcopy(znumdta, vart1, namatt,znumout, vart) 111 ENDFOR 112 ncdf_attput, znumout, zvarid, "units", units 113 ncdf_attput, znumout, zvarid, "long_name", zlong_name 111 ENDFOR 112 ncdf_attput, znumout, zvarid, "units", units 113 ncdf_attput, znumout, zvarid, "long_name", zlong_name 114 114 ; 115 115 ; ... Global attributes … … 119 119 ncdf_attput, znumout,"Title", ztitle, /global 120 120 ; 121 ; 3. Writing data 121 ; 3. Writing data 122 122 ; =============== 123 123 ; … … 131 131 ncdf_varput, znumout, varlon, zlon 132 132 ncdf_varput, znumout, varlat, zlat 133 IF keyword_set(ndim) THEN BEGIN 133 IF keyword_set(ndim) THEN BEGIN 134 134 ncdf_varget, nummsh, 'nav_lev', z, count = [jpkoce] 135 135 ncdf_varput, znumout, varlev, z 136 ENDIF 136 ENDIF 137 137 ; 138 138 ; 3.2. Time and mask … … 144 144 ncdf_varput, znumout, vart, zt 145 145 ;ncdf_varput, znumout, varmsk, zout_mask 146 if keyword_set(key_performance) THEN print, 'temps netcdf_output', systime(1)-tempsun 146 if keyword_set(key_performance) THEN print, 'temps netcdf_output', systime(1)-tempsun 147 147 ; 148 148 return 149 END 149 END -
trunk/northwind.pro
r41 r48 3 3 ; 4 4 ; PURPOSE: Check on in the intput field the level of coherence of winds nearby north pole, 5 ; and corrects it. Thi stask is performed before interpolation 5 ; and corrects it. Thi stask is performed before interpolation 6 6 ; 7 7 ; CATEGORY: Subroutine … … 9 9 ; CALLING SEQUENCE: northwind 10 10 ; 11 ; INPUTS: 12 ; 11 ; INPUTS: 13 12 ; 14 13 ; KEYWORD PARAMETERS: None 15 14 ; 16 ; OUTPUTS: 15 ; OUTPUTS: 17 16 ; zresul : the same wind field with its upward stripe 18 17 ; corrected in order to get coherence of winds on north polte 19 18 ; 20 19 ; COMMON BLOCKS: 21 ; common_interp.pro 22 ; 20 ; common_interp.pro 23 21 ; 24 22 PRO northwind, datglo,zdata_name,t … … 26 24 @common_interp 27 25 28 ; Treatment of north pole stripe 26 ; Treatment of north pole stripe 29 27 30 28 31 29 ; Get from global array, the last stripe and the last but one stripe 32 30 33 31 ; Initialising data arrays 34 32 35 33 datglo_rect = replicate(0.,2,jpiatm) 36 datglo_rectb = replicate(0.,2,jpiatm) 34 datglo_rectb = replicate(0.,2,jpiatm) 37 35 38 36 ; Put data of the last two stripes into arrays 39 37 40 38 41 39 datglo_rect(0,*) = datglo(*,jpjatm-1,0) 42 40 datglo_rect(1,*) = datglo(*,jpjatm-1,1) … … 45 43 datglo_rectb(1,*)= datglo(*,jpjatm-2,1) 46 44 47 45 48 46 49 47 ; Converts data of these two lines into polar coordinates … … 59 57 datglo_polarb=CV_COORD(FROM_RECT=datglo_rect,/TO_POLAR) 60 58 datglo_polarb(0,*)=datglo_polar(0,*)+2*!PI*indgen(jpiatm)/jpiatm 61 59 62 60 63 61 ; Coherence of the last line needs to be corrected … … 69 67 mean_argb=replicate((total(datglo_polarb(0,*))/jpiatm+2*!PI*west_lon/360.+!PI*(jpiatm-1.)/float(jpiatm)) ,jpiatm) 70 68 71 69 72 70 73 71 ; Checks whether module is about the same on the North Pole line : 74 72 75 73 test_mod=abs(mean_mod-datglo_polar(1,*)) 76 74 77 75 IF ( max(test_mod)/mean_mod(0) ) GT 0.01 THEN BEGIN 78 79 print, 'MODULE IS NOT COHERENT AND WILL BE CORRECTED TO A MEAN VALUE OF :', mean_mod(0) 80 76 77 print, 'MODULE IS NOT COHERENT AND WILL BE CORRECTED TO A MEAN VALUE OF :', mean_mod(0) 78 81 79 datglo_polar(1,*)=mean_mod 82 80 83 81 ENDIF 84 82 85 83 86 84 ; Checks whether argument is coherent on the North Pole line : … … 88 86 test_arg=abs(mean_arg-datglo_polar(0,*)) 89 87 90 IF ( max(test_arg)/mean_arg(0)) GT 0.01 THEN BEGIN 88 IF ( max(test_arg)/mean_arg(0)) GT 0.01 THEN BEGIN 91 89 92 90 print, 'ARGUMENT IS NOT COHERENT AND WILL BE CORRECTED TO A MEAN VALUE :', mean_argb(0) 93 91 94 92 datglo_polar(0,*)=mean_argb-!PI/2-(2*!PI*indgen(jpiatm)/float(jpiatm)+replicate(2.*!PI*west_lon/360.,jpiatm)) 95 96 93 94 97 95 98 96 ; Cross Corellation Method for wind coherence on North Pole 99 97 100 98 all_phases=indgen(jpiatm) 101 99 102 100 Ux=mean_mod*cos(all_phases*2.*!PI/jpiatm) 103 Uy=mean_mod*sin(all_phases*2.*!PI/jpiatm) 101 Uy=mean_mod*sin(all_phases*2.*!PI/jpiatm) 104 102 105 103 phyx=C_CORRELATE(Ux,datglo(*,jpjatm-2,0),all_phases) … … 110 108 deph_x=where(phyx EQ max(phyx)) 111 109 deph_y=where(phyy EQ max(phyy)) 112 110 113 111 deph_x=deph_x(0) 114 112 deph_y=deph_y(0) 115 113 116 114 117 115 118 116 … … 122 120 datglo(*,jpjatm-1,0)=Ux 123 121 datglo(*,jpjatm-1,1)=Uy 124 122 125 123 126 124 ENDIF 127 125 128 129 126 127 130 128 131 129 end -
trunk/preproc_mask.pro
r41 r48 3 3 ; NAME: preproc_mask.pro 4 4 ; 5 ; PURPOSE: Extends land-sea mask by one point over ocean 5 ; PURPOSE: Extends land-sea mask by one point over ocean 6 6 ; 7 7 ; CATEGORY: Function … … 9 9 ; CALLING SEQUENCE: preproc_mask, mask 10 10 ; 11 ; INPUTS: 12 ; mask : mask field composed of 0 and 1 values 13 ; with 0 over land points 11 ; INPUTS: 12 ; mask : mask field composed of 0 and 1 values 13 ; with 0 over land points 14 14 ; 15 15 ; KEYWORD PARAMETERS: None 16 16 ; 17 ; OUTPUTS: 18 ; mask : mask field extended 19 ; 17 ; OUTPUTS: 18 ; mask : mask field extended 20 19 ; 21 ; COMMON BLOCKS: 20 ; 21 ; COMMON BLOCKS: 22 22 ; common_interp 23 23 ; 24 ; SIDE EFFECTS: 24 ; SIDE EFFECTS: 25 25 ; 26 26 ; RESTRICTIONS: … … 50 50 ; 51 51 return, z 52 END 52 END -
trunk/rseries_ncdf.pro
r41 r48 65 65 ; "years since 1979-01-01 00:00:00" 66 66 ; 4) The calendar must be the Gregorian calendar 67 ; 5) The name of the file must begin ing by67 ; 5) The name of the file must beginning by 68 68 ; exp_freq_datefirst_datelast_* 69 69 ; 6) The maximum gap between 2 consecutive files must be one day. … … 118 118 if n_elements(freqin) EQ 0 then freq = '5d' ELSE freq = freqin 119 119 ;------------------------------------------------ 120 ; determina ition of the filename beginning with exp_freq_*120 ; determination of the filename beginning with exp_freq_* 121 121 ;------------------------------------------------ 122 122 IF keyword_set(gridtype) EQ 0 THEN gridtype = '' … … 141 141 ; date1 142 142 CASE strmid(freqin, 0, 1, /rever) OF 143 'm' :div = 100144 'y' :div = 10000145 ELSE :div = 1143 'm' : div = 100 144 'y' : div = 10000 145 ELSE : div = 1 146 146 ENDCASE 147 147 goodfile = where(datefirst le date1/div and datelast GE date1/div) … … 243 243 'm':datelast = 100*datelast + daysinmonth(datelast MOD 100, datelast / 100) 244 244 'y':div = 10000*datelast+360+(5+leapyr(datelast+century))*(key_caltype EQ 'greg') 245 ELSE :div = 1245 ELSE : div = 1 246 246 ENDCASE 247 247 if date2 GT datelast THEN BEGIN 248 248 ; if we need to read more than one file, 249 249 ; first we read the first file 250 ;++ print,var,date1+century*1000000L, datelast[0]+century*1000000L 251 ;++READ, B, PROMPT='Enter Name: ' 250 ;++ print,var,date1+century*1000000L, datelast[0]+century*1000000L 251 ;++READ, B, PROMPT='Enter Name: ' 252 252 res1 = read_ncdf(var, date1+century*1000000L, datelast[0]+century*1000000L $ 253 253 , filename = filename, /nostruct, _extra = ex) … … 269 269 res = temporary(res1) 270 270 END 271 ELSE :BEGIN271 ELSE : BEGIN 272 272 ; we glue the result of the first read and the result of the new call 273 273 ; to rseries_ncdf -
trunk/step1_diff.pro
r41 r48 114 114 ENDELSE 115 115 ; 116 ;117 116 ; write output 118 117 ; in order to avoid unexpected overwritten … … 122 121 RETURN 123 122 ENDIF 124 ;125 123 ; 126 124 ;--------------------------- -
trunk/step2_diff.sh
r47 r48 16 16 # ======== 17 17 # 18 # .. code-bl cok:: bash18 # .. code-block:: bash 19 19 # 20 20 # step2_diff.sh -y year … … 25 25 # For test on year 1994: 26 26 # 27 # .. code-bl cok:: bash27 # .. code-block:: bash 28 28 # 29 29 # ./step2_diff.sh -y 1994 … … 59 59 usage=" Usage : ${command} -y year" 60 60 # 61 set +u 62 while [ ! -z "${1}" ] 61 while [ ${#} -gt 0 ] 63 62 do 64 63 case ${1} in
Note: See TracChangeset
for help on using the changeset viewer.