Index: trunk/SRC/Interpolation/compute_fromirr_bilinear_weigaddr.pro
===================================================================
--- trunk/SRC/Interpolation/compute_fromirr_bilinear_weigaddr.pro	(revision 110)
+++ trunk/SRC/Interpolation/compute_fromirr_bilinear_weigaddr.pro	(revision 110)
@@ -0,0 +1,236 @@
+;+
+; @file_comments compute the weight and address needed to interpolate data from
+;          an "irregular 2D grid" (defined as a grid made of quadrilateral cells)
+;          to any grid using the bilinear method
+;   
+; @categories interpolation
+;
+;     @param olonin {in}{required} longitudeof the input data 
+;     @param olat   {in}{required} latitude of the input data 
+;     @param omsk   {in}{required} land/se mask of the input data 
+;     @param alonin {in}{required} longitude of the output data 
+;     @param alat   {in}{required} latitude of the output data 
+;     @param amsk   {in}{required} land/se mask of the output data 
+;
+; @returns 
+;     weig, addr: 2D arrays, weig and addr are the weight and addresses used to
+;     perform the interpolation:
+;          dataout = total(weig*datain[addr], 1)
+;          dataout = reform(dataout, jpia, jpja, /over)
+;
+; @restrictions
+;  -  the input grid must be an "irregular 2D grid", defined as a grid made 
+;     of quadrilateral cells which corners positions are defined with olonin and olat
+;  -  We supposed the data are located on a sphere, with a periodicity along
+;     the longitude
+;  -  to perform the bilinear interpolation within quadrilateral cells, we
+;     first morph the cell into a square cell and then compute the bilinear
+;     interpolation.
+;  -  if some corners of the cell are land points, their weight is set to 0 
+;     and the weight is redistributed on the remaining "water" corners
+;  -  points located out of the southern and northern boundaries or in cells
+;    containing only land points are set the the same value as their closest neighbor l; 
+; @history
+;  June 2006: Sebastien Masson (smasson\@lodyc.jussieu.fr) 
+; 
+;-
+;
+;----------------------------------------------------------
+;----------------------------------------------------------
+;
+PRO compute_fromirr_bilinear_weigaddr, olonin, olat, omsk, alonin, alat, amsk, weig, addr
+;
+  compile_opt idl2, strictarrsubs 
+;
+  jpia = (size(alonin, /dimensions))[0]
+  jpja = (size(alonin, /dimensions))[1]
+;
+  jpio = (size(olonin, /dimensions))[0]
+  jpjo = (size(olonin, /dimensions))[1]
+;
+; longitude, between 0 and 360
+  alon = alonin MOD 360
+  out = where(alon LT 0) 
+  IF out[0] NE -1 THEN alon[out] = alon[out]+360
+  olon = olonin MOD 360
+  out = where(olon LT 0) 
+  IF out[0] NE -1 THEN olon[out] = olon[out]+360
+;
+; we work only on the water points
+  owater = where(omsk EQ 1)
+  nowater = n_elements(owater) 
+  awater = where(amsk EQ 1)
+  nawater = n_elements(awater) 
+;
+; define all cells that have corners located at olon, olat
+; we define the cell with the address of all corners
+;
+;            3        2
+;             +------+
+;             |      |
+;             |      |
+;             |      |
+;             +------+
+;            0        1
+;
+  alladdr = lindgen(jpio, jpjo-1)
+  celladdr = lonarr(4, jpio*(jpjo-1))
+  celladdr[0, *] = alladdr
+  celladdr[1, *] = shift(alladdr, -1)
+  celladdr[2, *] = shift(alladdr + jpio, -1)
+  celladdr[3, *] = alladdr + jpio
+;
+; list the cells that have at least 1 ocean point as corner
+  good = where(total(omsk[celladdr], 1) GT 0)
+; keep only those cells 
+  celladdr = celladdr[*, temporary(good)]
+;
+  xcell = olon[celladdr]
+  minxcell = min(xcell, dimension = 1, max = maxxcell)
+  ycell = olat[celladdr]
+  minycell = min(ycell, dimension = 1, max = maxycell)
+; foraddr: address of the ocean water cell associated to each atmosphere water point
+  foraddr = lonarr(nawater)
+; forweight: x/y position of the atmosphere water point in the ocean water cell
+  forweight = dblarr(nawater, 2)
+;
+; Loop on all the water point of the atmosphere
+; We look for which ocean water cell contains the atmosphere water point
+;
+  delta = max([(360./jpio), (180./jpjo)])* 4.
+  FOR n = 0L, nawater-1 DO BEGIN
+; control print
+    IF (n MOD 5000) EQ 0 THEN print, n
+; longitute and latitude of the atmosphere water point
+    xx = alon[awater[n]]
+    yy = alat[awater[n]]
+; 1) we reduce the number of ocean cells for which we will try to know if
+; xx,yy is inside. 
+    CASE 1 OF
+; if we are near the norh pole
+      yy GE (90-delta):BEGIN
+        lat1 = 90-2*delta
+        good = where(maxycell GE lat1)
+        onsphere = 1
+      END
+; if we are near the longitude periodicity area
+      xx LE delta OR xx GE (360-delta):BEGIN 
+        lat1 = yy-delta
+        lat2 = yy+delta
+        good = where((minxcell LE 2*delta OR maxxcell GE (360-2*delta)) AND maxycell GE lat1 AND minycell LE lat2)
+        onsphere = 1
+      END
+; other cases
+      ELSE:BEGIN
+        lon1 = xx-delta
+        lon2 = xx+delta
+        lat1 = yy-delta
+        lat2 = yy+delta
+        good = where(maxxcell GE lon1 AND minxcell LE lon2 AND maxycell GE lat1 AND minycell le lat2)
+; ORCA cases : orca grid is irregular only northward of 40N
+        CASE 1 OF
+          jpio EQ 92   AND (jpjo EQ 76   OR jpjo EQ 75   OR jpjo EQ 74  ):onsphere = yy GT 40
+          jpio EQ 180  AND (jpjo EQ 149  OR jpjo EQ 148  OR jpjo EQ 147 ):onsphere = yy GT 40
+          jpio EQ 720  AND (jpjo EQ 522  OR jpjo EQ 521  OR jpjo EQ 520 ):onsphere = yy GT 40
+          jpio EQ 1440 AND (jpjo EQ 1021 OR jpjo EQ 1020 OR jpjo EQ 1019):onsphere = yy GT 40
+          ELSE:onsphere = 1
+        ENDCASE
+      END
+    ENDCASE
+; we found a short list of possible ocean water cells containing the atmosphere water point
+    IF good[0] NE -1 THEN BEGIN
+; in which cell is located the atmosphere water point?
+; Warning! inquad use clockwise quadrilateral definition
+      ind = inquad(xx, yy $
+                   , xcell[0, good], ycell[0, good] $
+                   , xcell[3, good], ycell[3, good] $
+                   , xcell[2, good], ycell[2, good] $
+                   , xcell[1, good], ycell[1, good] $
+                   , onsphere = onsphere, newcoord = newcoord, /noprint)
+; keep only the first cell (if the atmospheric point was located in several ocean cells)
+      ind = ind[0]
+; we found one ocean water cell containing the atmosphere water point
+      IF ind NE -1 THEN BEGIN
+        ind = good[ind]
+; we keep its address 
+        foraddr[n] = ind
+; now, we morph the quadrilateral ocean cell into the reference square (0 -> 1)
+; in addition we get the corrdinates of the atmospheric point in this new morphed square
+        IF onsphere THEN BEGIN
+; Warning! quadrilateral2square use anticlockwise quadrilateral definition
+          xy = quadrilateral2square(newcoord[0, 0], newcoord[1, 0] $
+                                    , newcoord[0, 3], newcoord[1, 3] $
+                                    , newcoord[0, 2], newcoord[1, 2] $
+                                    , newcoord[0, 1], newcoord[1, 1] $
+                                    , newcoord[0, 4], newcoord[1, 4])
+        ENDIF ELSE BEGIN
+          xy = quadrilateral2square(xcell[0, ind], ycell[0, ind] $
+                                    , xcell[1, ind], ycell[1, ind] $
+                                    , xcell[2, ind], ycell[2, ind] $
+                                    , xcell[3, ind], ycell[3, ind], xx, yy)
+        ENDELSE 
+; take care of rounding errors...
+        zero = where(abs(xy) LT 1e-4)
+        IF zero[0] NE -1 THEN xy[zero] = 0
+        one = where(abs(1-xy) LT 1e-4)
+        IF one[0] NE -1 THEN xy[one] = 1
+; some (useless) checks...
+        IF  xy[0] LT 0 OR xy[0] GT 1 THEN stop
+        IF  xy[0] LT 0 OR xy[0] GT 1 THEN stop
+; keep the new coordinates
+        forweight[n, 0] = xy[0]
+        forweight[n, 1] = xy[1]
+      ENDIF ELSE foraddr[n] = -1
+    ENDIF ELSE foraddr[n] = -1
+  ENDFOR
+; do we have some water atmospheric points that are not located in an water oceanic cell?
+  bad = where(foraddr EQ -1)
+  IF bad[0] NE -1 THEN BEGIN
+; yes?
+; we look for neighbor water atmospheric point located in water oceanic cell
+    badaddr = awater[bad]
+    good = where(foraddr NE -1)
+; list the atmospheric points located in water oceanic cell
+    goodaddr = awater[good]
+; there longitude and latitude
+    goodlon = alon[goodaddr]
+    goodlat = alat[goodaddr]
+; for all the bad points, look for a neighbor
+    neig = lonarr(n_elements(bad))
+    FOR i = 0, n_elements(bad)-1 DO BEGIN
+      neig[i] = (neighbor(alon[badaddr[i]], alat[badaddr[i]], goodlon, goodlat, /sphere))[0]
+    ENDFOR
+; get the address regarding foraddr
+    neig = good[neig]
+; associate each bad point with its neighbor (get its address and weight)
+    foraddr[bad] = foraddr[neig]
+    forweight[bad, *] = forweight[neig, *]
+  endif
+; transform the address of the ocean cell into the address of its 4 corners
+  newaaddr = celladdr[*, temporary(foraddr)]
+; now we compute the weight to give at each corner
+  newaweig = dblarr(4, nawater)
+  a = reform(forweight[*, 0], 1, nawater)
+  b = reform(forweight[*, 1], 1, nawater)
+  forweight =  -1 ; free memory
+  newaweig = [(1-a)*(1-b), (1-b)*a, a*b, (1-a)*b]
+  a = -1 &  b = -1 ; free memory
+; mask the weight to suppress the corner located on land
+  newaweig = newaweig*((omsk)[newaaddr])
+  totalweig = total(newaweig, 1)
+; for cell with some land corner,
+; we have to redistribute the weight on the reaining water corners
+; weights normalization
+  totalweig = total(newaweig, 1)
+  newaweig = newaweig/(replicate(1., 4)#totalweig)
+  totalweig = total(newaweig, 1)
+
+; weights
+  weig = dblarr(4, jpia*jpja)
+  weig[*, awater] = temporary(newaweig)
+; address
+  addr = dblarr(4, jpia*jpja)
+  addr[*, awater] = temporary(newaaddr)
+;
+  RETURN
+END
Index: trunk/SRC/Interpolation/compute_fromreg_bilinear_weigaddr.pro
===================================================================
--- trunk/SRC/Interpolation/compute_fromreg_bilinear_weigaddr.pro	(revision 109)
+++ trunk/SRC/Interpolation/compute_fromreg_bilinear_weigaddr.pro	(revision 110)
@@ -1,4 +1,4 @@
 ;+
-; @file_comments compute the weight and address neede to interpolate data from a
+; @file_comments compute the weight and address needed to interpolate data from a
 ;          "regular grid" to any grid using the bilinear method
 ;   
Index: trunk/SRC/Interpolation/fromirr.pro
===================================================================
--- trunk/SRC/Interpolation/fromirr.pro	(revision 110)
+++ trunk/SRC/Interpolation/fromirr.pro	(revision 110)
@@ -0,0 +1,95 @@
+;+
+;
+; @file_comments interpolate data from an irregular 2D grid to any 2D grid.
+;   Only 1 metod available: bilinear
+;   
+; @categories interpolation
+;
+; @examples  
+; dataout = fromirr(method, datain [, lonin, latin, mskin, lonout, latout, mskout])
+;
+;    @param method: {in}{required} a string defining the interpolation method. must be 'bilinear'
+;    @param datain: {in}{required} a 2D array the input data to interpolate
+;    @param lonin: {in}{optional} a 2D array defining the longitude of the input data
+;           optionals if WEIG and ADDR keywords used.
+;    @param latin: {in}{optional} a 2D array defining the latitude of the input data.
+;                     optionals if WEIG and ADDR keywords used.
+;    @param mskin: {in}{optional} a 2D array, the land-sea mask of the input data (1 on ocean, 0 on land)
+;    @param lonout: {in}{optional} 1D or 2D array defining the longitude of the output data.
+;                       optionals if WEIG and ADDR keywords used.
+;    @param latout: {in}{optional} 1D or 2D array defining the latitude of the output data.
+;                       optionals if WEIG and ADDR keywords used.
+;    @param mskout: {in}{optional} a 2D array, the land-sea mask of the ouput data (1 on ocean, 0 on land)
+;
+; @keyword WEIG (see ADDR)
+; @keyword ADDR: 2D arrays, weig and addr are the weight and addresses used to
+;     perform the interpolation:
+;          dataout = total(weig*datain[addr], 1)
+;          dataout = reform(dataout, jpio, jpjo, /over)
+;     Those keywords can be set to named variables into which the values will be
+;     copied when the current routine exits. Next, they can be used to perform
+;     the interpolation whithout computing again those 2 parameters. This greatly
+;     speed-up the interpolation! In that case, lonin, latin, lonout and latout are not necessary.
+;
+; @returns 2D array: the interpolated data
+;
+; @restrictions We supposed the data are located on a sphere, with a periodicity along
+;              the longitude.
+;              Note that the input data can contain the same cells several times
+;             (like ORCA grid near the north pole boundary)
+;
+; @examples
+;  
+;  tncep = fromirr('bilinear', topa, glamt, gphit, tmask[*,*,0], lonout, latout, mskout)
+;
+;  or 
+;
+;  t1ncep = fromirr('bilinear', topa, glamt, gphit, tmask[*,*,0], lonout, latout, mskout $
+;                            , WEIG = a, ADDR = b)
+;  help, a, b
+;  t2ncep = fromirr('bilinear', topa, WEIG = a, ADDR = b)
+;
+; @history
+;  June 2006: Sebastien Masson (smasson\@lodyc.jussieu.fr) 
+; 
+;-
+;----------------------------------------------------------
+;----------------------------------------------------------
+;
+FUNCTION fromirr, method, datain, lonin, latin, mskin, lonout, latout, mskout $
+                  , WEIG = weig, ADDR = addr
+;
+  compile_opt strictarr, strictarrsubs 
+;
+  IF NOT (keyword_set(weig) AND keyword_set(addr)) THEN BEGIN
+;---------------
+; atmospheric grid parameters
+;---------------
+    alon = lonin
+    alat = latin
+    get_gridparams, alon, alat, jpia, jpja, 2, /double
+;---------------
+; Oceanic grid parameters
+;---------------
+    olon = lonout
+    olat = latout
+    get_gridparams, olon, olat, jpio, jpjo, 2, /double
+;---------------
+; Compute weight and address
+;---------------
+    CASE method OF
+      'bilinear':compute_fromirr_bilinear_weigaddr, alon, alat, mskin, olon, olat, mskout, weig, addr
+      ELSE:BEGIN 
+        print, ' unknown interpolation method... we stop'
+        stop
+      ENDELSE
+    ENDCASE
+  ENDIF
+;---------------
+; to the interpolation
+;---------------
+  dataout = total(weig*datain[addr], 1)
+  dataout = reform(dataout, jpio, jpjo, /over)
+;
+  RETURN, dataout
+END
Index: trunk/SRC/Interpolation/fromreg.pro
===================================================================
--- trunk/SRC/Interpolation/fromreg.pro	(revision 109)
+++ trunk/SRC/Interpolation/fromreg.pro	(revision 110)
@@ -14,10 +14,15 @@
 ;            must be 'bilinear' or 'imoms3'
 ;    @param datain {in}{required}  a 2D array the input data to interpolate
-;    @param lonin latin {in}{required}  longitude/latitude of the input data. optionals if
+;    @param lonin {in}{optional}  longitude of the input data. optionals if
 ;            WEIG and ADDR keywords used.
-;    @param lonout latout {in}{required}  longitude/latitude of the output data. optionals if
+;    @param latin {in}{optional}  latitude of the input data. optionals if
+;            WEIG and ADDR keywords used.
+;    @param lonout {in}{optional}  longitude of the output data. optionals if
+;            WEIG and ADDR keywords used.
+;    @param lonout {in}{optional}  latitude of the output data. optionals if
 ;            WEIG and ADDR keywords used.
 ;
-; @keyword     WEIG, ADDR 2D arrays, weig and addr are the weight and addresses used to
+; @keyword     WEIG (see ADDR)
+; @keyword     ADDR 2D arrays, weig and addr are the weight and addresses used to
 ;     perform the interpolation:
 ;          dataout = total(weig*datain[addr], 1)
@@ -51,5 +56,4 @@
 ; 
 ;-
-;
 ;----------------------------------------------------------
 ;----------------------------------------------------------
Index: trunk/SRC/Interpolation/inquad.pro
===================================================================
--- trunk/SRC/Interpolation/inquad.pro	(revision 109)
+++ trunk/SRC/Interpolation/inquad.pro	(revision 110)
@@ -285,4 +285,5 @@
   IF keyword_set(onsphere) THEN BEGIN
     IF arg_present(newcoord) THEN BEGIN
+      found = found[0]
       newcoord = [[x1[found], y1[found]] $
                   , [x2[found], y2[found]] $
