;------------------------------------------------------------ ;------------------------------------------------------------ ;------------------------------------------------------------ ; ; NAME: make_weight.pro ; ; PURPOSE: Creates for a given ORCA grid, and its ORCA_GEO equivalent ; grid the required ; ; CATEGORY : Subroutine ; ; CALLING SEQUENCE : make_weight ; ; INPUTS : ; All the variables of a given ORCA meshmask file and its ORCA_GEO equivalent ; OUTPUTS : ; Weights and pointers required to make interpolation between ; a field on the ORCA T grid and a field on the ORCA_GEO T grid ; ; COMMON BLOCKS: ; None ; ; SIDE EFFECTS: ; ; RESTRICTIONS: ; ; EXAMPLE: ; ; MODIFICATION HISTORY: 08/2002 Robinson Hordoir ; ;------------------------------------------------------------ ;------------------------------------------------------------ ;------------------------------------------------------------ pro make_weight @initorca @naminterp2 @init_path2 ;***************************************** ; Error file opening ;***************************************** openu, lun_err,errorfile,/get_lun ;************************ ; Compute Geogrid ;************************ ; Let's get our reference latitude northernmost=where(gphit EQ max(gphit)) jnorth=northernmost/jpi inorth=min(northernmost-jnorth*jpi) ref_latitude=gphit(inorth,*) ; Let's get the reference longitude equador=where(abs(gphit) EQ min(abs(gphit)) ) ref_longitude=glamt(equador) ; Now we have the new grid glamt_reg=ref_longitude#replicate(1,n_elements(ref_latitude)) gphit_reg=replicate(1,n_elements(ref_longitude))#ref_latitude ; Get the output mask mask_reg=read_ncdf('tmask',FILENAME=outputgrid,/NOSTRUCT,/cont_nofill) ; Get the latitude and longitude with extra boundary lines get_fullmesh, glamt_full,gphit_full sph_coord_full=replicate(9999.,3,jpiglo*jpjglo) sph_coord_full(0,*)=reform(glamt_full,jpiglo*jpjglo) sph_coord_full(1,*)=reform(gphit_full,jpiglo*jpjglo) sph_coord_full(2,*)=replicate(1.,jpiglo*jpjglo) rec_grid_full=cv_coord(/DEGREES,FROM_SPHERE=sph_coord_full,/TO_RECT) xx_grid_full=rec_grid_full(0,*) yy_grid_full=rec_grid_full(1,*) zz_grid_full=rec_grid_full(2,*) sph_coord=replicate(9999.,3,jpi*jpj) radius=replicate(1.,jpi,jpj) radius(0:jpi/2,jpj-1)=1000. sph_coord(0,*)=reform(glamt,jpi*jpj) sph_coord(1,*)=reform(gphit,jpi*jpj) sph_coord(2,*)=reform(radius,jpi*jpj) rec_grid=cv_coord(/DEGREES,FROM_SPHERE=sph_coord,/TO_RECT) xx_grid=rec_grid(0,*) yy_grid=rec_grid(1,*) zz_grid=rec_grid(2,*) ; Init array to stock pointers pointer_Tx=replicate(999,jpi,jpj,4) pointer_Ty=replicate(999,jpi,jpj,4) ; Init weight matric weight_T=double(replicate(0.,jpi,jpj,4)) ; Get the latitude to start from not_found=1 for jj=0, jpj-1 do begin lat_max_ORCA=max( gphit(*,jj) ) lat_max_reg =max( gphit_reg(*,jj) ) if ((lat_max_ORCA NE lat_max_reg) AND (not_found EQ 1)) then begin jpj_start=jj-2 not_found=0 endif end for jj=0,jpj-1 do begin for ji=0,jpi-1 do begin if mask_reg(ji,jj) NE 0 then begin ; Coordinates of our point xx=glamt_reg(ji,jj) yy=gphit_reg(ji,jj) sph_coord=replicate(1.,3) sph_coord(0)=xx sph_coord(1)=yy ; Convert them into cartesian rec_coord=cv_coord(/DEGREES,FROM_SPHERE=sph_coord,/TO_RECT) xx_bis=replicate(rec_coord(0),jpi*jpj) yy_bis=replicate(rec_coord(1),jpi*jpj) zz_bis=replicate(rec_coord(2),jpi*jpj) ;**************************************************************** ; Finds in which grid cell point stand ;**************************************************************** xx_bis=reform(xx_bis,jpi,jpj) yy_bis=reform(yy_bis,jpi,jpj) zz_bis=reform(zz_bis,jpi,jpj) xx_grid_full=reform(xx_grid_full,jpiglo,jpjglo) yy_grid_full=reform(yy_grid_full,jpiglo,jpjglo) zz_grid_full=reform(zz_grid_full,jpiglo,jpjglo) check_inside,xx_bis ,yy_bis ,zz_bis,$ 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),$ xx_grid_full(2:jpiglo-1,0:jpjglo-2),yy_grid_full(2:jpiglo-1,0:jpjglo-2),zz_grid_full(2:jpiglo-1,0:jpjglo-2),$ xx_grid_full(2:jpiglo-1,1:jpjglo-1),yy_grid_full(2:jpiglo-1,1:jpjglo-1),zz_grid_full(2:jpiglo-1,1:jpjglo-1),$ 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 if answer(0) EQ -1 then begin print, 'ERROR IN MAKE WEIGHT : POSITION NOT LOCATED ON GRID' stop endif answer_j=answer/jpi answer_i=answer-answer_j*jpi answer_j=answer_j(where(answer_i EQ min(answer_i))) answer_i=answer_i(where(answer_i EQ min(answer_i))) answer_j=answer_j(0) answer_i=answer_i(0) pointer_Tx(ji,jj,0)=answer_i pointer_Ty(ji,jj,0)=answer_j point_x1=answer_i point_y1=answer_j pointer_Tx(ji,jj,1)=answer_i+1 pointer_Ty(ji,jj,1)=answer_j point_x2=answer_i+1 point_y2=answer_j pointer_Tx(ji,jj,2)=answer_i+1 pointer_Ty(ji,jj,2)=answer_j+1 point_x3=answer_i+1 point_y3=answer_j+1 pointer_Tx(ji,jj,3)=answer_i pointer_Ty(ji,jj,3)=answer_j+1 point_x4=answer_i point_y4=answer_j+1 ; Now we know which 4 points are located nearby our output T point. ; Let's get the coordinates of this points in long and lat glamt_ref=replicate(0.,4) gphit_ref=replicate(0.,4) glamt_ref(0)=glamt_full(point_x1+1,point_y1) gphit_ref(0)=gphit_full(point_x1+1,point_y1) glamt_ref(1)=glamt_full(point_x2+1,point_y2) gphit_ref(1)=gphit_full(point_x2+1,point_y2) glamt_ref(2)=glamt_full(point_x3+1,point_y3) gphit_ref(2)=gphit_full(point_x3+1,point_y3) glamt_ref(3)=glamt_full(point_x4+1,point_y4) gphit_ref(3)=gphit_full(point_x4+1,point_y4) ; Let's make things simple for coordinates sph_coord_ref=replicate(0.,3,5) sph_coord_ref(0,0)=glamt_reg(ji,jj) sph_coord_ref(1,0)=gphit_reg(ji,jj) sph_coord_ref(2,0)=1. sph_coord_ref(0,1:4)=glamt_ref sph_coord_ref(1,1:4)=gphit_ref sph_coord_ref(2,1:4)=1. rec_coord_ref=cv_coord(/DEGREES,FROM_SPHERE=sph_coord_ref,/TO_RECT) xx=rec_coord_ref(0,0) yy=rec_coord_ref(1,0) zz=rec_coord_ref(2,0) x1=rec_coord_ref(0,1) y1=rec_coord_ref(1,1) z1=rec_coord_ref(2,1) x2=rec_coord_ref(0,2) y2=rec_coord_ref(1,2) z2=rec_coord_ref(2,2) x3=rec_coord_ref(0,3) y3=rec_coord_ref(1,3) z3=rec_coord_ref(2,3) x4=rec_coord_ref(0,4) y4=rec_coord_ref(1,4) z4=rec_coord_ref(2,4) xG=total(rec_coord_ref(0,1:4))/4. yG=total(rec_coord_ref(1,1:4))/4. zG=total(rec_coord_ref(2,1:4))/4. rec_coord_ref(0,*)=rec_coord_ref(0,*)-xG rec_coord_ref(1,*)=rec_coord_ref(1,*)-yG rec_coord_ref(2,*)=rec_coord_ref(2,*)-zG XT=rec_coord_ref(*,1:4) XM=MATRIX_MULTIPLY(XT,XT,/BTRANSPOSE) EVAL=HQR(ELMHES(XM),/DOUBLE) EVEC=EIGENVEC(XM,EVAL) EVAL=FLOAT(EVAL) EVEC=FLOAT(EVEC) ORDER=SORT(EVAL) xx=total(rec_coord_ref(*,0)*EVEC(*,ORDER(1))) yy=total(rec_coord_ref(*,0)*EVEC(*,ORDER(2))) x1=total(rec_coord_ref(*,1)*EVEC(*,ORDER(1))) y1=total(rec_coord_ref(*,1)*EVEC(*,ORDER(2))) x2=total(rec_coord_ref(*,2)*EVEC(*,ORDER(1))) y2=total(rec_coord_ref(*,2)*EVEC(*,ORDER(2))) x3=total(rec_coord_ref(*,3)*EVEC(*,ORDER(1))) y3=total(rec_coord_ref(*,3)*EVEC(*,ORDER(2))) x4=total(rec_coord_ref(*,4)*EVEC(*,ORDER(1))) y4=total(rec_coord_ref(*,4)*EVEC(*,ORDER(2))) ;************************************************** ; Computes i coordinate of point in the grid cell * ;************************************************** ; Calculate coordinates of j vector ; Vector j0 is vector AD/(AD) xj0=x4-x1 yj0=y4-y1 ; Normalised vector norm_j0=sqrt(xj0*xj0+yj0*yj0) xj0=xj0/norm_j0 yj0=yj0/norm_j0 ; Vector j1 is vector BC/(BC) xj1=x3-x2 yj1=y3-y2 ; Normalised vector norm_j1=sqrt(xj1*xj1+yj1*yj1) xj1=xj1/norm_j1 yj1=yj1/norm_j1 ; Caculate coefficient for solving 2nd degree equation in i a_coeff= (x2-x1)*(yj0-yj1)-(y2-y1)*(xj0-xj1) b_coeff= -(x2-x1)*yj0-(xx-x1)*(yj0-yj1)+(yy-y1)*(xj0-xj1)+(y2-y1)*xj0 c_coeff= (xx-x1)*yj0-(yy-y1)*xj0 delta=max([b_coeff*b_coeff-4*a_coeff*c_coeff,0]) iP=999. if abs(a_coeff) GT 1.E-5 then begin sol1=(-b_coeff+sqrt(delta))/(2.*a_coeff) sol2=(-b_coeff-sqrt(delta))/(2.*a_coeff) dist_sol1=sol1*sol1+(sol1-1)*(sol1-1) dist_sol2=sol2*sol2+(sol2-1)*(sol2-1) iP=sol1 if dist_sol1 LT dist_sol2 then iP=sol1 if dist_sol2 LT dist_sol1 then iP=sol2 endif else begin iP=-c_coeff/b_coeff endelse ;************************************************** ; Computes j coordinate of point in the grid cell * ;************************************************** ; Calculate coordinates of i vector ; Vector i0 is vector AB/(AB) xi0=x2-x1 yi0=y2-y1 ; Normalised vector norm_i0=sqrt(xi0*xi0+yi0*yi0) if norm_i0 LT 0. then stio xi0=xi0/norm_i0 yi0=yi0/norm_i0 ; Vector i1 is vector DC/(DC) xi1=x3-x4 yi1=y3-y4 ; Normalised vector norm_i1=sqrt(xi1*xi1+yi1*yi1) if norm_i1 LT 0. then stop xi1=xi1/norm_i1 yi1=yi1/norm_i1 ; Caculate coefficient for solving 2nd degree equation in j a_coeff= (x4-x1)*(yi0-yi1)-(y4-y1)*(xi0-xi1) b_coeff= -(x4-x1)*yi0-(xx-x1)*(yi0-yi1)+(yy-y1)*(xi0-xi1)+(y4-y1)*xi0 c_coeff= (xx-x1)*yi0-(yy-y1)*xi0 delta=max([b_coeff*b_coeff-4*a_coeff*c_coeff,0]) jP=999. if abs(a_coeff) GT 1.E-5 then begin sol1=(-b_coeff+sqrt(delta))/(2.*a_coeff) sol2=(-b_coeff-sqrt(delta))/(2.*a_coeff) dist_sol1=sol1*sol1+(sol1-1)*(sol1-1) dist_sol2=sol2*sol2+(sol2-1)*(sol2-1) jP=sol1 if dist_sol1 LT dist_sol2 then jP=sol1 if dist_sol2 LT dist_sol1 then jP=sol2 endif else begin jP=-c_coeff/b_coeff endelse weight_T(ji,jj,0) = (1.-iP)*(1.-jP) weight_T(ji,jj,1) = iP*(1.-jP) weight_T(ji,jj,2) = iP*jP weight_T(ji,jj,3) = (1.-iP)*jP if min(weight_T(ji,jj,*)) LT 0 OR max(weight_T(ji,jj,*)) GT 1 then begin printf, lun_err,'JI=',ji,' JJ=',jj printf, lun_err,'IP=',iP,' JP=',jP printf, lun_err,'XP=',xx,' YP=',yy printf, lun_err,'XA=',x1,' YA=',y1 printf, lun_err,'XB=',x2,' YB=',y2 printf, lun_err,'XC=',x3,' YC=',y3 printf, lun_err,'XD=',x4,' YD=',y4 endif endif end end close, lun_err free_lun,lun_err ; Convert to full format weight_full=replicate(0.,jpiglo,jpjglo,4) pointer_i_full=replicate(0,jpiglo,jpjglo,4) pointer_j_full=replicate(0,jpiglo,jpjglo,4) weight_full(1:jpiglo-2,0:jpjglo-2,*)=weight_T pointer_i_full(1:jpiglo-2,0:jpjglo-2,*)=pointer_Tx pointer_j_full(1:jpiglo-2,0:jpjglo-2,*)=pointer_Ty ; Let's save Weight functions and the pointer arrays ; Creates Netcdf File to save this vargrid = 'T' ; Name idout = NCDF_CREATE(weightfile,/clobber) NCDF_CONTROL, idout, /nofill ; Dimension xidout = NCDF_DIMDEF(idout, 'x', jpiglo) yidout = NCDF_DIMDEF(idout, 'y', jpjglo) didout = NCDF_DIMDEF(idout, 'deptht', 1) tidout = NCDF_DIMDEF(idout, 'time_counter', /unlimited) ; Attributes id_0 = NCDF_VARDEF(idout, 'nav_lon' , [xidout, yidout ], /FLOAT) id_1 = NCDF_VARDEF(idout, 'nav_lat' , [xidout, yidout ], /FLOAT) id_2 = NCDF_VARDEF(idout, 'deptht' , [ didout ], /FLOAT) id_3 = NCDF_VARDEF(idout, 'time_counter', [ tidout], /FLOAT) id0 = NCDF_VARDEF(idout, 'weight_1' , [xidout, yidout,didout,tidout ], /FLOAT) id1 = NCDF_VARDEF(idout, 'weight_2' , [xidout, yidout,didout,tidout ], /FLOAT) id2 = NCDF_VARDEF(idout, 'weight_3' , [xidout, yidout,didout,tidout ], /FLOAT) id3 = NCDF_VARDEF(idout, 'weight_4' , [xidout, yidout,didout,tidout ], /FLOAT) id4 = NCDF_VARDEF(idout, 'pointer_1_i' , [xidout, yidout ,didout,tidout ], /LONG) id5 = NCDF_VARDEF(idout, 'pointer_2_i' , [xidout, yidout ,didout,tidout ], /LONG) id6 = NCDF_VARDEF(idout, 'pointer_3_i' , [xidout, yidout,didout,tidout ], /LONG) id7 = NCDF_VARDEF(idout, 'pointer_4_i' , [xidout, yidout,didout,tidout ], /LONG) id8 = NCDF_VARDEF(idout, 'pointer_1_j' , [xidout, yidout,didout,tidout ], /LONG) id9 = NCDF_VARDEF(idout, 'pointer_2_j' , [xidout, yidout,didout,tidout ], /LONG) id10 = NCDF_VARDEF(idout, 'pointer_3_j' , [xidout, yidout,didout,tidout ], /LONG) id11 = NCDF_VARDEF(idout, 'pointer_4_j' , [xidout, yidout,didout,tidout ], /LONG) ; Variable 0 NCDF_ATTPUT, idout, id_0, 'units', 'degrees_east' NCDF_ATTPUT, idout, id_0, 'long_name', 'Longitude' NCDF_ATTPUT, idout, id_0, 'nav_model', 'Default grid' ; Variable 1 NCDF_ATTPUT, idout, id_1, 'units', 'degrees_north' NCDF_ATTPUT, idout, id_1, 'long_name', 'Latitude' NCDF_ATTPUT, idout, id_1, 'nav_model', 'Default grid' ; Variable 2 NCDF_ATTPUT, idout, id_2, 'units','meters' NCDF_ATTPUT, idout, id_2, 'long_name','Depth' NCDF_ATTPUT, idout, id_2, 'nav_model','Default grid' ; Variable3 NCDF_ATTPUT, idout, id_3, 'units', 'seconds since 0001-01-15 12:00:00 ' NCDF_ATTPUT, idout, id_3, 'calendar','noleap' NCDF_ATTPUT, idout, id_3, 'title', 'Time' NCDF_ATTPUT, idout, id_3, 'long_name', 'Time axis' NCDF_ATTPUT, idout, id_3, 'time_origin','0000-DEC-15 00:00:00' ; Variable 0 NCDF_ATTPUT, idout, id0, 'units','no unit' NCDF_ATTPUT, idout, id0, 'long_name','W1' ; Variable 1 NCDF_ATTPUT, idout, id1, 'units','no unit' NCDF_ATTPUT, idout, id1, 'long_name','W2' ; Variable 0 NCDF_ATTPUT, idout, id2, 'units','no unit' NCDF_ATTPUT, idout, id2, 'long_name','W3' ; Variable 0 NCDF_ATTPUT, idout, id3, 'units','no unit' NCDF_ATTPUT, idout, id3, 'long_name','W4' ; Variable 0 NCDF_ATTPUT, idout, id4, 'units','no unit' NCDF_ATTPUT, idout, id4, 'long_name','i1' ; Variable 1 NCDF_ATTPUT, idout, id5, 'units','no unit' NCDF_ATTPUT, idout, id5, 'long_name','i2' ; Variable 0 NCDF_ATTPUT, idout, id6, 'units','no unit' NCDF_ATTPUT, idout, id6, 'long_name','i3' ; Variable 0 NCDF_ATTPUT, idout, id7, 'units','no unit' NCDF_ATTPUT, idout, id7, 'long_name','i4' ; Variable 0 NCDF_ATTPUT, idout, id8, 'units','no unit' NCDF_ATTPUT, idout, id8, 'long_name','j1' ; Variable 1 NCDF_ATTPUT, idout, id9, 'units','no unit' NCDF_ATTPUT, idout, id9, 'long_name','j2' ; Variable 0 NCDF_ATTPUT, idout, id10, 'units','no unit' NCDF_ATTPUT, idout, id10, 'long_name','j3' ; Variable 0 NCDF_ATTPUT, idout, id11, 'units','no unit' NCDF_ATTPUT, idout, id11, 'long_name','j4' NCDF_CONTROL, idout, /ENDEF ; Writting NCDF_VARPUT, idout, id_0, glamt NCDF_VARPUT, idout, id_1, gphit prof=5. NCDF_VARPUT, idout, id_2, prof temps=lonarr(12) for i=0,0 do temps[i] = 86400l*(julday(1+i,15,1)-julday(1,15,1)) NCDF_VARPUT, idout, id_3, temps(0) NCDF_VARPUT, idout, id0,weight_full(*,*,0) NCDF_VARPUT, idout, id1,weight_full(*,*,1) NCDF_VARPUT, idout, id2,weight_full(*,*,2) NCDF_VARPUT, idout, id3,weight_full(*,*,3) NCDF_VARPUT, idout, id4,pointer_i_full(*,*,0) NCDF_VARPUT, idout, id5,pointer_i_full(*,*,1) NCDF_VARPUT, idout, id6,pointer_i_full(*,*,2) NCDF_VARPUT, idout, id7,pointer_i_full(*,*,3) NCDF_VARPUT, idout, id8,pointer_j_full(*,*,0) NCDF_VARPUT, idout, id9,pointer_j_full(*,*,1) NCDF_VARPUT, idout, id10,pointer_j_full(*,*,2) NCDF_VARPUT, idout, id11,pointer_j_full(*,*,3) NCDF_CLOSE, idout end