Index: trunk/INTERP2/make_weight.pro
===================================================================
--- trunk/INTERP2/make_weight.pro	(revision 2)
+++ trunk/INTERP2/make_weight.pro	(revision 48)
@@ -6,5 +6,5 @@
 ;
 ; PURPOSE: Creates for a given ORCA grid, and its ORCA_GEO equivalent
-;          grid  the required 
+;          grid the required
 ;
 ; CATEGORY : Subroutine
@@ -36,123 +36,123 @@
 
 
-@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)
+    @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))
+    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
+    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)
-   
+
+        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 
-    
+; Convert them into cartesian
+
     rec_coord=cv_coord(/DEGREES,FROM_SPHERE=sph_coord,/TO_RECT)
 
@@ -162,10 +162,10 @@
 
 
-    
+
 
 ;****************************************************************
-; Finds in which grid cell point stand   
+; Finds in which grid cell point stand
 ;****************************************************************
-    
+
     xx_bis=reform(xx_bis,jpi,jpj)
     yy_bis=reform(yy_bis,jpi,jpj)
@@ -176,5 +176,5 @@
     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),$
@@ -183,6 +183,6 @@
                  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
 
@@ -191,21 +191,21 @@
 
     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 
+    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
@@ -223,5 +223,5 @@
     point_y4=answer_j+1
 
-   
+
 
 ; Now we know which 4 points are located nearby our output T point.
@@ -231,24 +231,24 @@
     glamt_ref=replicate(0.,4)
     gphit_ref=replicate(0.,4)
-  
-  
-   
-    glamt_ref(0)=glamt_full(point_x1+1,point_y1)              
+
+
+
+    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)              
+    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)              
+
+    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)              
+
+    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)
 
@@ -262,5 +262,5 @@
 
    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)
@@ -270,5 +270,5 @@
    y1=rec_coord_ref(1,1)
    z1=rec_coord_ref(2,1)
- 
+
    x2=rec_coord_ref(0,2)
    y2=rec_coord_ref(1,2)
@@ -286,5 +286,5 @@
    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
@@ -292,5 +292,5 @@
 
    XT=rec_coord_ref(*,1:4)
- 
+
    XM=MATRIX_MULTIPLY(XT,XT,/BTRANSPOSE)
    EVAL=HQR(ELMHES(XM),/DOUBLE)
@@ -299,6 +299,6 @@
    EVAL=FLOAT(EVAL)
    EVEC=FLOAT(EVEC)
-   ORDER=SORT(EVAL)  
-   
+   ORDER=SORT(EVAL)
+
 
    xx=total(rec_coord_ref(*,0)*EVEC(*,ORDER(1)))
@@ -316,5 +316,5 @@
    x4=total(rec_coord_ref(*,4)*EVEC(*,ORDER(1)))
    y4=total(rec_coord_ref(*,4)*EVEC(*,ORDER(2)))
-   
+
 
 ;**************************************************
@@ -358,5 +358,5 @@
 
 b_coeff= -(x2-x1)*yj0-(xx-x1)*(yj0-yj1)+(yy-y1)*(xj0-xj1)+(y2-y1)*xj0
-   
+
 c_coeff= (xx-x1)*yj0-(yy-y1)*xj0
 
@@ -432,5 +432,5 @@
 
 b_coeff= -(x4-x1)*yi0-(xx-x1)*(yi0-yi1)+(yy-y1)*(xi0-xi1)+(y4-y1)*xi0
-   
+
 c_coeff= (xx-x1)*yi0-(yy-y1)*xi0
 
@@ -441,5 +441,5 @@
 jP=999.
 
-if abs(a_coeff) GT 1.E-5  then begin
+if abs(a_coeff) GT 1.E-5 then begin
 
 sol1=(-b_coeff+sqrt(delta))/(2.*a_coeff)
@@ -461,6 +461,6 @@
 
 endelse
-   
-   
+
+
    weight_T(ji,jj,0) = (1.-iP)*(1.-jP)
 
@@ -471,5 +471,5 @@
    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
 
@@ -477,5 +477,5 @@
       printf, lun_err,'IP=',iP,' JP=',jP
 
-     
+
       printf, lun_err,'XP=',xx,' YP=',yy
       printf, lun_err,'XA=',x1,' YA=',y1
@@ -483,45 +483,45 @@
       printf, lun_err,'XC=',x3,' YC=',y3
       printf, lun_err,'XD=',x4,' YD=',y4
-      
-   endif
-  
-
-   endif
-
-  
-  end
+
+               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
+    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
+    ; Dimension
    xidout = NCDF_DIMDEF(idout, 'x', jpiglo)
    yidout = NCDF_DIMDEF(idout, 'y', jpjglo)
@@ -531,10 +531,10 @@
 
 
-; Attributes   
+; 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)
@@ -573,5 +573,5 @@
    NCDF_ATTPUT, idout, id0, 'units','no unit'
    NCDF_ATTPUT, idout, id0, 'long_name','W1'
- 
+
 ; Variable 1
    NCDF_ATTPUT, idout, id1, 'units','no unit'
@@ -588,5 +588,5 @@
    NCDF_ATTPUT, idout, id4, 'units','no unit'
    NCDF_ATTPUT, idout, id4, 'long_name','i1'
- 
+
 ; Variable 1
    NCDF_ATTPUT, idout, id5, 'units','no unit'
@@ -604,5 +604,5 @@
    NCDF_ATTPUT, idout, id8, 'units','no unit'
    NCDF_ATTPUT, idout, id8, 'long_name','j1'
- 
+
 ; Variable 1
    NCDF_ATTPUT, idout, id9, 'units','no unit'
@@ -636,5 +636,5 @@
    NCDF_VARPUT, idout, id_3, temps(0)
 
- 
+
    NCDF_VARPUT, idout, id0,weight_full(*,*,0)
    NCDF_VARPUT, idout, id1,weight_full(*,*,1)
@@ -649,5 +649,5 @@
    NCDF_VARPUT, idout, id10,pointer_j_full(*,*,2)
    NCDF_VARPUT, idout, id11,pointer_j_full(*,*,3)
-  
+
 
 
