Changeset 48 for trunk/INTERP2/make_weight.pro
- Timestamp:
- 03/16/14 20:38:39 (10 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
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
Note: See TracChangeset
for help on using the changeset viewer.