Index: trunk/Interpolation/.idlwave_catalog
===================================================================
--- trunk/Interpolation/.idlwave_catalog	(revision 68)
+++ trunk/Interpolation/.idlwave_catalog	(revision 69)
@@ -2,5 +2,5 @@
 ;; IDLWAVE catalog for library saxo
 ;; Automatically Generated -- do not edit.
-;; Created by idlwave_catalog on Tue Jan 24 17:06:19 2006
+;; Created by idlwave_catalog on Thu May  4 08:42:55 2006
 ;;
 (setq idlwave-library-catalog-libname "saxo")
@@ -10,5 +10,5 @@
    ("clickincell" fun nil (lib "clickincell.pro" nil "saxo") "Result = %s" (nil ("_EXTRA") ("CELLTYPE") ("COLOR") ("DRAWCELL") ("IJ") ("ORIGINAL")))
    ("compute_fromreg_bilinear_weigaddr" pro nil (lib "compute_fromreg_bilinear_weigaddr.pro" nil "saxo") "%s, alonin, alatin, olonin, olat, weig, addr" (nil ("NONORTHERNLINE") ("NOSOUTHERNLINE")))
-   ("compute_fromreg_imoms3_weigaddr" pro nil (lib "compute_fromreg_imoms3_weigaddr.pro" nil "saxo") "%s, alonin, alatin, olonin, olat, weig, addr" (nil ("FORCE") ("NONORTHERNLINE") ("NOSOUTHERNLINE")))
+   ("compute_fromreg_imoms3_weigaddr" pro nil (lib "compute_fromreg_imoms3_weigaddr.pro" nil "saxo") "%s, alonin, alatin, olonin, olat, weig, addr" (nil ("NONORTHERNLINE") ("NOSOUTHERNLINE")))
    ("cutpar" fun nil (lib "cutpar.pro" nil "saxo") "Result = %s(x0, y0, x1, y1, x2, y2, x3, y3, n)" (nil ("endpoints") ("onsphere")))
    ("cutsegment" fun nil (lib "cutsegment.pro" nil "saxo") "Result = %s(x0, y0, x1, y1, n)" (nil ("endpoints") ("onsphere")))
@@ -23,4 +23,9 @@
    ("neighbor" fun nil (lib "neighbor.pro" nil "saxo") "Result = %s(p0lon, p0lat, neighlon, neighlat)" (nil ("distance") ("radians") ("sphere")))
    ("quadrilateral2square" fun nil (lib "quadrilateral2square.pro" nil "saxo") "Result = %s(x0in, y0in, x1in, y1in, x2in, y2in, x3in, y3in, xxin, yyin)" (nil ("PERF")))
+   ("spl_fstdrv" fun nil (lib "spl_fstdrv.pro" nil "saxo") "Result = %s(x, y, yscd, x2)" (nil))
+   ("pure_concave" fun nil (lib "spl_incr.pro" nil "saxo") "Result = %s(x1, x2, y1, y2, der2, x)" (nil))
+   ("pure_convex" fun nil (lib "spl_incr.pro" nil "saxo") "Result = %s(x1, x2, y1, y2, der2, x)" (nil))
+   ("spl_incr" fun nil (lib "spl_incr.pro" nil "saxo") "Result = %s(x, y, x2)" (nil ("YP0") ("YPN_1")))
+   ("spl_keep_mean" fun nil (lib "spl_keep_mean.pro" nil "saxo") "Result = %s(x, yin, x2)" (nil ("GE0") ("YP0") ("YPN_1")))
    ("square2quadrilateral" fun nil (lib "square2quadrilateral.pro" nil "saxo") "Result = %s(x0in, y0in, x1in, y1in, x2in, y2in, x3in, y3in, xxin, yyin)" (nil))
    ("testinterp" pro nil (lib "testinterp.pro" nil "saxo") "%s" (nil))))
Index: trunk/Interpolation/clickincell.pro
===================================================================
--- trunk/Interpolation/clickincell.pro	(revision 68)
+++ trunk/Interpolation/clickincell.pro	(revision 69)
@@ -19,5 +19,5 @@
 ;     CELLTYPE = 'T', 'W', 'U', 'V' or 'F': This this the type of point
 ;     that is located in the center of the cell which the click is
-;     located. defalut is T type of cell (with corner defined by F
+;     located. default is T type of cell (with corner defined by F
 ;     points).
 ;
Index: trunk/Interpolation/compute_fromreg_imoms3_weigaddr.pro
===================================================================
--- trunk/Interpolation/compute_fromreg_imoms3_weigaddr.pro	(revision 68)
+++ trunk/Interpolation/compute_fromreg_imoms3_weigaddr.pro	(revision 69)
@@ -31,6 +31,7 @@
 ;
 ; RESTRICTIONS:
-;  -  the input grid must be a "regular grid", defined as a grid for which
-;     lontitudes and latitudes are regularly spaced.
+;  -  the input grid must be a "regular/rectangular grid", defined as a grid for
+;     which each lontitudes lines have the same latitude and each latitudes columns
+;     have the same longitude.
 ;  -  We supposed the data are located on a sphere, with a periodicity along
 ;     the longitude.
@@ -45,5 +46,5 @@
 ; MODIFICATION HISTORY:
 ;  November 2005: Sebastien Masson (smasson@lodyc.jussieu.fr) 
-; 
+;  March 2006: works for rectangular grids
 ;-
 ;
@@ -52,5 +53,5 @@
 ;
 PRO compute_fromreg_imoms3_weigaddr, alonin, alatin, olonin, olat, weig, addr $
-  , NONORTHERNLINE = nonorthernline, NOSOUTHERNLINE = nosouthernline, FORCE = force
+                                     , NONORTHERNLINE = nonorthernline, NOSOUTHERNLINE = nosouthernline
 ;
   compile_opt strictarr, strictarrsubs 
@@ -75,11 +76,8 @@
     IF array_equal(sort(alon), lindgen(jpia)) NE 1 THEN stop
   ENDIF ELSE shiftx = 0
-; alon must be regularly spaced
+; alon is it regularly spaced?
   step = alon-shift(alon, 1)
   step[0] = step[0] + 360.
-  IF total((step-step[0]) GE 1.e-6) NE 0 THEN BEGIN 
-    print,  'input longitude must be regularly spaced, we stop...'
-    stop
-  ENDIF
+  IF total((step-step[0]) GE 1.e-6) NE 0 THEN noregx = 1
 ; we extend the longitude range of alon (-> easy interpolation even
 ; near minalon et maxalon)
@@ -92,11 +90,8 @@
 ; alat must be monotonically increasing
   IF array_equal(sort(alat), lindgen(jpja)) NE 1 THEN stop
-; alat must be regularly spaced
+; alat is it regularly spaced?
   step = alat-shift(alat, 1)
   step = step[1:jpja - 1L]
-  IF total((step-step[0]) GE 1.e-6) NE 0 AND NOT keyword_set(force) THEN BEGIN 
-    print,  'input latitude must be regularly spaced, we stop...'
-    stop
-  ENDIF
+  IF total((step-step[0]) GE 1.e-6) NE 0 THEN noregy = 1
 ;
   if keyword_set(nonorthernline) then BEGIN 
@@ -143,24 +138,42 @@
   ilat = indexlat[short]
 ; 
-  delta = alat[ilat+1L]-alat[ilat]
-  IF max(abs(delta-delta[0])) GE 1.e-6 THEN stop
-  delta = delta[0]
-;
-  d0 = (alat[ilat-1]-olat[short])/delta
-  IF min(d0, max = ma) LE -2 THEN stop
-  IF ma GT -1 THEN stop
-  wy0 = imoms3(temporary(d0))
-  d1 = (alat[ilat  ]-olat[short])/delta
-  IF min(d1, max = ma) LE -1 THEN stop
-  IF ma GT 0 THEN stop
-  wy1 = imoms3(temporary(d1))
-  d2 = (alat[ilat+1]-olat[short])/delta
-  IF min(d2, max = ma) LE 0 THEN stop
-  IF ma GT 1 THEN stop  
-  wy2 = imoms3(temporary(d2))
-  d3 = (alat[ilat+2]-olat[short])/delta
-  IF min(d3, max = ma) LE 1 THEN stop
-  IF ma GT 2 THEN stop  
-  wy3 = imoms3(temporary(d3))
+  IF NOT keyword_set(noregy) THEN BEGIN 
+    delta = alat[ilat+1L]-alat[ilat]
+    IF max(abs(delta-delta[0])) GE 1.e-6 THEN stop
+    delta = delta[0]
+;
+    d0 = (alat[ilat-1L]-olat[short])/delta
+    IF min(d0, max = ma) LE -2 THEN stop
+    IF ma GT -1 THEN stop
+    wy0 = imoms3(temporary(d0))
+    d1 = (alat[ilat   ]-olat[short])/delta
+    IF min(d1, max = ma) LE -1 THEN stop
+    IF ma GT 0 THEN stop
+    wy1 = imoms3(temporary(d1))
+    d2 = (alat[ilat+1L]-olat[short])/delta
+    IF min(d2, max = ma) LE 0 THEN stop
+    IF ma GT 1 THEN stop  
+    wy2 = imoms3(temporary(d2))
+    d3 = (alat[ilat+2L]-olat[short])/delta
+    IF min(d3, max = ma) LE 1 THEN stop
+    IF ma GT 2 THEN stop  
+    wy3 = imoms3(temporary(d3))
+  ENDIF ELSE BEGIN 
+    nele = n_elements(short)
+    wy0 = fltarr(nele)
+    wy1 = fltarr(nele)
+    wy2 = fltarr(nele)
+    wy3 = fltarr(nele)
+    FOR i = 0L, nele-1 DO BEGIN
+      IF i MOD 10000 EQ 0 THEN print, i
+      newlat = spl_incr(alat[ilat[i]-1L:ilat[i]+2L], [-1., 0., 1., 2.], olat[short[i]])
+      IF newlat LE 0 THEN stop
+      IF newlat GT 1 THEN stop
+      wy0[i] = imoms3(newlat+1)
+      wy1[i] = imoms3(newlat)
+      wy2[i] = imoms3(1-newlat)
+      wy3[i] = imoms3(2-newlat)
+    ENDFOR
+  ENDELSE 
 ;
   mi = min(wy0+wy1+wy2+wy3, max = ma)
@@ -168,24 +181,42 @@
   IF abs(ma-1) GE 1.e-6 THEN stop
 ;
-  delta = alon[ilon]-alon[ilon-1]
-  IF max(abs(delta-delta[0])) GE 1.e-6 THEN stop
-  delta = delta[0]
-;
-  d0 = (alon[ilon-1]-olon[short])/delta
-  IF min(d0, max = ma) LE -2 THEN stop
-  IF ma GT -1 THEN stop
-  wx0 = imoms3(temporary(d0))
-  d1 = (alon[ilon  ]-olon[short])/delta
-  IF min(d1, max = ma) LE -1 THEN stop
-  IF ma GT 0 THEN stop
-  wx1 = imoms3(temporary(d1))
-  d2 = (alon[ilon+1]-olon[short])/delta
-  IF min(d2, max = ma) LE 0 THEN stop
-  IF ma GT 1 THEN stop
-  wx2 = imoms3(temporary(d2))
-  d3 = (alon[ilon+2]-olon[short])/delta
-  IF min(d3, max = ma) LE 1 THEN stop
-  IF ma GT 2 THEN stop  
-  wx3 = imoms3(temporary(d3))
+  IF NOT keyword_set(noregx) THEN BEGIN 
+    delta = alon[ilon]-alon[ilon-1L]
+    IF max(abs(delta-delta[0])) GE 1.e-6 THEN stop
+    delta = delta[0]
+;
+    d0 = (alon[ilon-1L]-olon[short])/delta
+    IF min(d0, max = ma) LE -2 THEN stop
+    IF ma GT -1 THEN stop
+    wx0 = imoms3(temporary(d0))
+    d1 = (alon[ilon   ]-olon[short])/delta
+    IF min(d1, max = ma) LE -1 THEN stop
+    IF ma GT 0 THEN stop
+    wx1 = imoms3(temporary(d1))
+    d2 = (alon[ilon+1L]-olon[short])/delta
+    IF min(d2, max = ma) LE 0 THEN stop
+    IF ma GT 1 THEN stop
+    wx2 = imoms3(temporary(d2))
+    d3 = (alon[ilon+2L]-olon[short])/delta
+    IF min(d3, max = ma) LE 1 THEN stop
+    IF ma GT 2 THEN stop  
+    wx3 = imoms3(temporary(d3))
+  ENDIF ELSE BEGIN 
+    nele = n_elements(short)
+    wx0 = fltarr(nele)
+    wx1 = fltarr(nele)
+    wx2 = fltarr(nele)
+    wx3 = fltarr(nele)
+    FOR i = 0L, nele-1 DO BEGIN
+      IF i MOD 10000 EQ 0 THEN print, i
+      newlon = spl_incr(alon[ilon[i]-1L:ilon[i]+2L], [-1., 0., 1., 2.], olon[short[i]])
+      IF newlon LE 0 THEN stop
+      IF newlon GT 1 THEN stop
+      wx0[i] = imoms3(newlon+1)
+      wx1[i] = imoms3(newlon)
+      wx2[i] = imoms3(1-newlon)
+      wx3[i] = imoms3(2-newlon)
+    ENDFOR
+  ENDELSE 
 ;
   mi = min(wx0+wx1+wx2+wx3, max = ma)
@@ -224,10 +255,10 @@
   xaddr[10, short] = ilon + 1L
   xaddr[11, short] = ilon + 2L
-  yaddr[8 , short] = ilat + 1L
-  yaddr[9 , short] = yaddr[8 , short]
-  yaddr[10, short] = yaddr[8 , short]
-  yaddr[11, short] = yaddr[8 , short]
-  weig[8 , short] = wx0 * wy2
-  weig[9 , short] = wx1 * wy2
+  yaddr[8, short] = ilat + 1L
+  yaddr[9, short] = yaddr[8, short]
+  yaddr[10, short] = yaddr[8, short]
+  yaddr[11, short] = yaddr[8, short]
+  weig[8, short] = wx0 * wy2
+  weig[9, short] = wx1 * wy2
   weig[10, short] = wx2 * wy2
   weig[11, short] = wx3 * wy2
@@ -254,18 +285,20 @@
 ; linear interpolation between line 1 and line 2
 ;
-  short = where(indexlat EQ jpja-2 OR indexlat EQ 0)
+  short = where(indexlat EQ jpja-2L OR indexlat EQ 0)
   IF short[0] NE -1 THEN BEGIN
     ilon = indexlon[short]
     ilat = indexlat[short]
 ;
-    delta = alat[ilat+1]-alat[ilat]
-    IF max(abs(delta-delta[0])) GE 1.e-6 THEN stop
-    delta = delta[0]
-;
-    d1 = (alat[ilat  ]-olat[short])/delta
+    delta = alat[ilat+1L]-alat[ilat]
+    IF NOT keyword_set(noregy) THEN BEGIN 
+      IF max(abs(delta-delta[0])) GE 1.e-6 THEN stop
+      delta = delta[0] 
+    ENDIF 
+;
+    d1 = (alat[ilat   ]-olat[short])/delta
     IF min(d1, max = ma) LE -1 THEN stop
     IF ma GT 0 THEN stop
     wy1 = 1.+ temporary(d1)
-    d2 = (alat[ilat+1]-olat[short])/delta
+    d2 = (alat[ilat+1L]-olat[short])/delta
     IF min(d2, max = ma) LE 0 THEN stop
     IF ma GT 1 THEN stop  
@@ -276,24 +309,42 @@
     IF abs(ma-1) GE 1.e-6 THEN stop
 ; but imoms3 along the longitude
-    delta = alon[ilon]-alon[ilon-1]
-    IF max(abs(delta-delta[0])) GE 1.e-6 THEN stop
-    delta = delta[0]
-;
-    d0 = (alon[ilon-1]-olon[short])/delta
-    IF min(d0, max = ma) LE -2 THEN stop
-    IF ma GT -1 THEN stop
-    wx0 = imoms3(temporary(d0))
-    d1 = (alon[ilon  ]-olon[short])/delta
-    IF min(d1, max = ma) LE -1 THEN stop
-    IF ma GT 0 THEN stop
-    wx1 = imoms3(temporary(d1))
-    d2 = (alon[ilon+1]-olon[short])/delta
-    IF min(d2, max = ma) LE 0 THEN stop
-    IF ma GT 1 THEN stop
-    wx2 = imoms3(temporary(d2))
-    d3 = (alon[ilon+2]-olon[short])/delta
-    IF min(d3, max = ma) LE 1 THEN stop
-    IF ma GT 2 THEN stop  
-    wx3 = imoms3(temporary(d3))
+    IF NOT keyword_set(noregx) THEN BEGIN 
+      delta = alon[ilon]-alon[ilon-1L]
+      IF max(abs(delta-delta[0])) GE 1.e-6 THEN stop
+      delta = delta[0]
+;
+      d0 = (alon[ilon-1L]-olon[short])/delta
+      IF min(d0, max = ma) LE -2 THEN stop
+      IF ma GT -1 THEN stop
+      wx0 = imoms3(temporary(d0))
+      d1 = (alon[ilon   ]-olon[short])/delta
+      IF min(d1, max = ma) LE -1 THEN stop
+      IF ma GT 0 THEN stop
+      wx1 = imoms3(temporary(d1))
+      d2 = (alon[ilon+1L]-olon[short])/delta
+      IF min(d2, max = ma) LE 0 THEN stop
+      IF ma GT 1 THEN stop
+      wx2 = imoms3(temporary(d2))
+      d3 = (alon[ilon+2L]-olon[short])/delta
+      IF min(d3, max = ma) LE 1 THEN stop
+      IF ma GT 2 THEN stop  
+      wx3 = imoms3(temporary(d3))
+    ENDIF ELSE BEGIN 
+      nele = n_elements(short)
+      wx0 = fltarr(nele)
+      wx1 = fltarr(nele)
+      wx2 = fltarr(nele)
+      wx3 = fltarr(nele)
+      FOR i = 0L, nele-1 DO BEGIN
+        IF i MOD 10000 EQ 0 THEN print, i
+        newlon = spl_incr(alon[ilon[i]-1L:ilon[i]+2L], [-1., 0., 1., 2.], olon[short[i]])
+        IF newlon LE 0 THEN stop
+        IF newlon GT 1 THEN stop
+        wx0[i] = imoms3(newlon+1)
+        wx1[i] = imoms3(newlon)
+        wx2[i] = imoms3(1-newlon)
+        wx3[i] = imoms3(2-newlon)
+      ENDFOR
+    ENDELSE 
 ;
     mi = min(wx0+wx1+wx2+wx3, max = ma)
@@ -340,24 +391,42 @@
     ilon = indexlon[short]
 ;
-    delta = alon[ilon]-alon[ilon-1]
-    IF max(abs(delta-delta[0])) GE 1.e-6 THEN stop
-    delta = delta[0]
-;
-    d0 = (alon[ilon-1]-olon[short])/delta
-    IF min(d0, max = ma) LE -2 THEN stop
-    IF ma GT -1 THEN stop
-    wx0 = imoms3(temporary(d0))
-    d1 = (alon[ilon  ]-olon[short])/delta
-    IF min(d1, max = ma) LE -1 THEN stop
-    IF ma GT 0 THEN stop
-    wx1 = imoms3(temporary(d1))
-    d2 = (alon[ilon+1]-olon[short])/delta
-    IF min(d2, max = ma) LE 0 THEN stop
-    IF ma GT 1 THEN stop
-    wx2 = imoms3(temporary(d2))
-    d3 = (alon[ilon+2]-olon[short])/delta
-    IF min(d3, max = ma) LE 1 THEN stop
-    IF ma GT 2 THEN stop  
-    wx3 = imoms3(temporary(d3))
+    IF NOT keyword_set(noregx) THEN BEGIN 
+      delta = alon[ilon]-alon[ilon-1L]
+      IF max(abs(delta-delta[0])) GE 1.e-6 THEN stop
+      delta = delta[0]
+;
+      d0 = (alon[ilon-1L]-olon[short])/delta
+      IF min(d0, max = ma) LE -2 THEN stop
+      IF ma GT -1 THEN stop
+      wx0 = imoms3(temporary(d0))
+      d1 = (alon[ilon   ]-olon[short])/delta
+      IF min(d1, max = ma) LE -1 THEN stop
+      IF ma GT 0 THEN stop
+      wx1 = imoms3(temporary(d1))
+      d2 = (alon[ilon+1L]-olon[short])/delta
+      IF min(d2, max = ma) LE 0 THEN stop
+      IF ma GT 1 THEN stop
+      wx2 = imoms3(temporary(d2))
+      d3 = (alon[ilon+2L]-olon[short])/delta
+      IF min(d3, max = ma) LE 1 THEN stop
+      IF ma GT 2 THEN stop  
+      wx3 = imoms3(temporary(d3))
+    ENDIF ELSE BEGIN 
+      nele = n_elements(short)
+      wx0 = fltarr(nele)
+      wx1 = fltarr(nele)
+      wx2 = fltarr(nele)
+      wx3 = fltarr(nele)
+      FOR i = 0L, nele-1 DO BEGIN
+        IF i MOD 10000 EQ 0 THEN print, i
+        newlon = spl_incr(alon[ilon[i]-1L:ilon[i]+2L], [-1., 0., 1., 2.], olon[short[i]])
+        IF newlon LE 0 THEN stop
+        IF newlon GT 1 THEN stop
+        wx0[i] = imoms3(newlon+1)
+        wx1[i] = imoms3(newlon)
+        wx2[i] = imoms3(1-newlon)
+        wx3[i] = imoms3(2-newlon)
+      ENDFOR
+    ENDELSE 
 ;
     mi = min(wx0+wx1+wx2+wx3, max = ma)
@@ -384,28 +453,46 @@
 ; Interpolation only along the longitude 
 ;
-  short = where(indexlat EQ jpja-1)
+  short = where(indexlat EQ jpja-1L)
   IF short[0] NE -1 THEN BEGIN
     ilon = indexlon[short]
 ;
-    delta = alon[ilon]-alon[ilon-1]
-    IF max(abs(delta-delta[0])) GE 1.e-6 THEN stop
-    delta = delta[0]
-;
-    d0 = (alon[ilon-1]-olon[short])/delta
-    IF min(d0, max = ma) LE -2 THEN stop
-    IF ma GT -1 THEN stop
-    wx0 = imoms3(temporary(d0))
-    d1 = (alon[ilon  ]-olon[short])/delta
-    IF min(d1, max = ma) LE -1 THEN stop
-    IF ma GT 0 THEN stop
-    wx1 = imoms3(temporary(d1))
-    d2 = (alon[ilon+1]-olon[short])/delta
-    IF min(d2, max = ma) LE 0 THEN stop
-    IF ma GT 1 THEN stop
-    wx2 = imoms3(temporary(d2))
-    d3 = (alon[ilon+2]-olon[short])/delta
-    IF min(d3, max = ma) LE 1 THEN stop
-    IF ma GT 2 THEN stop  
-    wx3 = imoms3(temporary(d3))
+    IF NOT keyword_set(noregx) THEN BEGIN 
+      delta = alon[ilon]-alon[ilon-1L]
+      IF max(abs(delta-delta[0])) GE 1.e-6 THEN stop
+      delta = delta[0]
+;
+      d0 = (alon[ilon-1L]-olon[short])/delta
+      IF min(d0, max = ma) LE -2 THEN stop
+      IF ma GT -1 THEN stop
+      wx0 = imoms3(temporary(d0))
+      d1 = (alon[ilon   ]-olon[short])/delta
+      IF min(d1, max = ma) LE -1 THEN stop
+      IF ma GT 0 THEN stop
+      wx1 = imoms3(temporary(d1))
+      d2 = (alon[ilon+1L]-olon[short])/delta
+      IF min(d2, max = ma) LE 0 THEN stop
+      IF ma GT 1 THEN stop
+      wx2 = imoms3(temporary(d2))
+      d3 = (alon[ilon+2L]-olon[short])/delta
+      IF min(d3, max = ma) LE 1 THEN stop
+      IF ma GT 2 THEN stop  
+      wx3 = imoms3(temporary(d3))
+    ENDIF ELSE BEGIN 
+      nele = n_elements(short)
+      wx0 = fltarr(nele)
+      wx1 = fltarr(nele)
+      wx2 = fltarr(nele)
+      wx3 = fltarr(nele)
+      FOR i = 0L, nele-1 DO BEGIN
+        IF i MOD 10000 EQ 0 THEN print, i
+        newlon = spl_incr(alon[ilon[i]-1L:ilon[i]+2L], [-1., 0., 1., 2.], olon[short[i]])
+        IF newlon LE 0 THEN stop
+        IF newlon GT 1 THEN stop
+        wx0[i] = imoms3(newlon+1)
+        wx1[i] = imoms3(newlon)
+        wx2[i] = imoms3(1-newlon)
+        wx3[i] = imoms3(2-newlon)
+      ENDFOR
+    ENDELSE 
 ;
     mi = min(wx0+wx1+wx2+wx3, max = ma)
Index: trunk/Interpolation/extrapolate.pro
===================================================================
--- trunk/Interpolation/extrapolate.pro	(revision 68)
+++ trunk/Interpolation/extrapolate.pro	(revision 69)
@@ -42,5 +42,5 @@
   sqrtinv = 1./sqrt(2)
 ;
-  cnt = 0
+  cnt = 1
 ; When we look for the coast line, we don't whant to select the
 ; borderlines of the array. -> we force the value of the mask for
Index: trunk/Interpolation/fromreg.pro
===================================================================
--- trunk/Interpolation/fromreg.pro	(revision 68)
+++ trunk/Interpolation/fromreg.pro	(revision 69)
@@ -2,15 +2,8 @@
 ; NAME: fromreg
 ;
-; PURPOSE: interpolate data from a "regular grid" to any grid.
+; PURPOSE: interpolate data from a "regular/rectangular grid" to any grid.
 ;   2 metods availables: bilinear and imoms3 
-; 
-;   In the case on 'bilinear interpolation':
-;   A "regular grid" is defined as a grid for which each lontitudes lines have 
+;   A "regular/rectangular grid" is defined as a grid for which each lontitudes lines have 
 ;   the same latitude and each latitudes columns have the same longitude.
-;   
-; 
-;   In the case on 'imoms3 interpolation':
-;   A "regular grid" is defined as a grid for which lontitudes and latitudes
-;   are regularly spaced.
 ;   
 ; CATEGORY:interpolation
Index: trunk/Interpolation/get_gridparams.pro
===================================================================
--- trunk/Interpolation/get_gridparams.pro	(revision 68)
+++ trunk/Interpolation/get_gridparams.pro	(revision 69)
@@ -6,5 +6,5 @@
 ;      and make sure it is 1D or 2D arrays 
 ;   
-;   or 2) given longitude and latidude arrays get their dimensions and make sure 
+;   or 2) given longitude and latitude arrays get their dimensions and make sure 
 ;      they are 1D or 2D arrays 
 ;
@@ -13,9 +13,9 @@
 ; CALLING SEQUENCE:
 ; 
-; 1) ncdf_get_gridparams, file, lonname, latname, lon, lat, jpi, jpj, n_dimensions
+; 1) get_gridparams, file, lonname, latname, lon, lat, jpi, jpj, n_dimensions
 ;
 ; or 
 ;
-; 2) ncdf_get_gridparams, lon, lat, jpi, jpj, n_dimensions
+; 2) get_gridparams, lon, lat, jpi, jpj, n_dimensions
 ;
 ; INPUTS:
@@ -29,5 +29,5 @@
 ;
 ; 2) lon and lat: 1d or 2D arrays defining longitudes and latitudes. 
-;    Note that these arrays are also outpust and can therefore be modified. 
+;    Note that these arrays are also outputs and can therefore be modified. 
 
 ; KEYWORD PARAMETERS: none
@@ -35,5 +35,5 @@
 ; OUTPUTS:
 ;  lon the variable that will contain the longitudes
-;  lat the variable that will contain the latidutes
+;  lat the variable that will contain the latitudes
 ;  jpi the number of points in the longitudinal direction
 ;  jpj the number of points in the latitudinal direction
Index: trunk/Interpolation/imoms3.pro
===================================================================
--- trunk/Interpolation/imoms3.pro	(revision 68)
+++ trunk/Interpolation/imoms3.pro	(revision 69)
@@ -7,10 +7,10 @@
   IF test1[0] NE -1 THEN BEGIN
     xtmp = x[test1]
-    y[test1] = .5*xtmp^3 -xtmp^2 -.5*xtmp+1
+    y[test1] = 0.5*xtmp*xtmp*xtmp -xtmp*xtmp -0.5*xtmp + 1
   ENDIF
   test1 = where(x LT 2 AND x GE 1)
   IF test1[0] NE -1 THEN BEGIN
     xtmp = x[test1]
-    y[test1] = (-1./6.)*xtmp^3 + xtmp^2 - (11./6.)*xtmp+1
+    y[test1] = (-1./6.)*xtmp*xtmp*xtmp + xtmp*xtmp - (11./6.)*xtmp + 1
   ENDIF
 
Index: trunk/Interpolation/spl_fstdrv.pro
===================================================================
--- trunk/Interpolation/spl_fstdrv.pro	(revision 69)
+++ trunk/Interpolation/spl_fstdrv.pro	(revision 69)
@@ -0,0 +1,92 @@
+;------------------------------------------------------------
+;------------------------------------------------------------
+;------------------------------------------------------------
+;+
+; NAME:spl_fstdrv
+;
+; PURPOSE: SPL_FSTDRV returns the values of the first derivative of
+; the interpolating function at the points X2i. it is a double
+; precision array.
+;
+; Given the arrays X and Y, which tabulate a function (with the X[i]
+; AND Y[i] in ascending order), and given an input value X2, the
+; SPL_INCR function returns an interpolated value for the given values
+; of X2. The interpolation method is based on cubic spline, corrected
+; in a way that interpolated value are also in ascending order
+;
+; CATEGORY:
+;
+; CALLING SEQUENCE: y2 =  spl_fstdrv(x, y, yscd, x2)
+;
+; INPUTS:
+;
+;    x: An n-element (at least 2) input vector that specifies the
+;    tabulate points in ascending order.
+;
+;    y: f(x) = y. An n-element input vector that specifies the values
+;    of the tabulated function F(Xi) corresponding to Xi.
+;
+;    yscd: The output from SPL_INIT for the specified X and Y.
+;
+;    x2: The input values for which the first derivative values are
+;    desired. X can be scalar or an array of values.
+
+; KEYWORD PARAMETERS: none
+;
+; OUTPUTS: 
+;
+;    y2: f'(x2) = y2. 
+;
+; COMMON BLOCKS: none
+;
+; SIDE EFFECTS: ?
+;
+; RESTRICTIONS: ?
+;
+; EXAMPLE:
+;
+; MODIFICATION HISTORY:
+;  Sebastien Masson (smasson@lodyc.jussieu.fr): May 2005
+;-
+;------------------------------------------------------------
+;------------------------------------------------------------
+;------------------------------------------------------------
+FUNCTION spl_fstdrv, x, y, yscd, x2
+;
+; compute the first derivative of the spline function;
+;
+  nx = n_elements(x)
+  ny = n_elements(y)
+; x must have at least 2 elements
+  IF nx LT 2 THEN stop   
+; y must have the same number of elements than x
+  IF nx NE ny THEN stop
+; define loc in a way that 
+;  if loc[i] eq -1   :                 x2[i] <  x[0]
+;  if loc[i] eq nx2-1:                 x2[i] >= x[nx-1]
+;  else              :    x[loc[i]] <= x2[i] <  x[loc[i]+1]
+  loc = value_locate(x, x2)
+; change loc in order to 
+; use x[0]    and x[1]    even if x2[i] <  x[0]    -> extrapolation
+; use x[nx-2] and x[nx-1] even if x2[i] >= x[nx-1] -> extrapolation
+  loc = 0 > temporary(loc) < (nx-2) 
+; distance between to consecutive x
+  deltax = x[loc+1]-x[loc]
+; distance between to consecutive y
+  deltay = y[loc+1]-y[loc]
+; relative distance between x2[i] and x[loc[i]+1]
+  a = (x[loc+1]-x2)/deltax
+; relative distance between x2[i] and x[loc[i]]
+  b = 1.0d - a
+; compute the first derivative on x (see numerical recipes Chap 3.3)
+  yfrst = temporary(deltay)/deltax $
+    - 1.0d/6.0d * (3.0d*a*a - 1.0d) * deltax * yscd[loc] $
+    + 1.0d/6.0d * (3.0d*b*b - 1.0d) * deltax * yscd[loc+1]
+; beware of the computation precision... 
+; force near zero values to be exactly 0.0
+  zero = where(abs(yfrst) LT 1.e-10)
+  IF zero[0] NE -1 THEN yfrst[zero] = 0.0d
+
+  RETURN, yfrst
+END
+
Index: trunk/Interpolation/spl_incr.pro
===================================================================
--- trunk/Interpolation/spl_incr.pro	(revision 69)
+++ trunk/Interpolation/spl_incr.pro	(revision 69)
@@ -0,0 +1,496 @@
+;------------------------------------------------------------
+;------------------------------------------------------------
+;------------------------------------------------------------
+;+
+; NAME:spl_incr
+;
+; PURPOSE:
+;
+; Given the arrays X and Y, which tabulate a function (with the X[i]
+; AND Y[i] in ascending order), and given an input value X2, the
+; SPL_INCR function returns an interpolated value for the given values
+; of X2. The interpolation method is based on cubic spline, corrected
+; in a way that interpolated values are also monotonically increasing.
+;
+; CATEGORY:
+;
+; CALLING SEQUENCE: y2 =  spl_incr(x, y, x2)
+;
+; INPUTS:
+;
+;    x: An n-element (at least 2) input vector that specifies the
+;    tabulate points in a strict ascending order.
+;
+;    y: f(x) = y. An n-element input vector that specifies the values
+;    of the tabulated function F(Xi) corresponding to Xi. As f is
+;    supposed to be monotonically increasing, y values must be
+;    monotonically increasing. y can have equal consecutive values.
+;
+;    x2: The input values for which the interpolated values are
+;    desired. Its values must be strictly monotonically increasing. 
+;
+; KEYWORD PARAMETERS:
+;
+;    YP0: The first derivative of the interpolating function at the
+;    point X0. If YP0 is omitted, the second derivative at the
+;    boundary is set to zero, resulting in a "natural spline."
+;
+;    YPN_1: The first derivative of the interpolating function at the
+;    point Xn-1. If YPN_1 is omitted, the second derivative at the
+;    boundary is set to zero, resulting in a "natural spline." 
+;
+; OUTPUTS: 
+;
+;    y2: f(x2) = y2. Double precision array
+;
+; COMMON BLOCKS: none
+;
+; SIDE EFFECTS: ?
+;
+; RESTRICTIONS:
+;   It might be possible that y2[i+1]-y2[i] has very small negative
+;   values (amplitude smaller than 1.e-6)...
+;
+; EXAMPLE:
+;
+;     n = 100L
+;     x = (dindgen(n))^2
+;     y = abs(randomn(0, n))
+;     y[n/2:n/2+1] = 0.
+;     y[n-n/3] = 0.
+;     y[n-n/6:n-n/6+5] = 0.
+;     y = total(y, /cumulative, /double)
+;     x2 = dindgen((n-1)^2)
+;     n2 = n_elements(x2)
+;     print, min(y[1:n-1]-y[0:n-2]) LT 0
+;     y2 = spl_incr( x, y, x2)
+;     splot, x, y, xstyle = 1, ystyle = 1, ysurx=.25, petit = [1, 2, 1], /land
+;     oplot, x2, y2, color = 100
+;     c = y2[1:n2-1] - y2[0:n2-2]
+;     print, min(c) LT 0 
+;     print, min(c, max = ma), ma
+;     splot,c,xstyle=1,ystyle=1, yrange=[-.01,.05], ysurx=.25, petit = [1, 2, 2], /noerase
+;     oplot,[0, n_elements(c)], [0, 0], linestyle = 1
+;
+; MODIFICATION HISTORY:
+;  Sebastien Masson (smasson@lodyc.jussieu.fr): May-Dec 2005
+;-
+;------------------------------------------------------------
+;------------------------------------------------------------
+;------------------------------------------------------------
+FUNCTION pure_concave, x1, x2, y1, y2, der2, x
+; X^n type
+  xx = (double(x)-double(x1))/(double(x2)-double(x1))
+  f = (double(x2)-double(x1))/(double(y2)-double(y1))
+  n = der2*temporary(f)
+  res = xx^(n)
+;   IF check_math() GT 0 THEN BEGIN
+;       zero = where(abs(res) LT 1.e-10)
+;       IF zero[0] NE -1 THEN res[zero] = 0.0d
+;   END
+  res = temporary(res)*(double(y2)-double(y1))+y1
+;
+;  IF array_equal(sort(res), lindgen(n_elements(res)) ) NE 1 THEN stop
+  RETURN, res
+END
+
+FUNCTION pure_convex, x1, x2, y1, y2, der2, x
+; 1-(1-X)^n type
+  xx = 1.0d - (double(x)-double(x1))/(double(x2)-double(x1))
+  f = (double(x2)-double(x1))/(double(y2)-double(y1))
+  n = der2*temporary(f)
+  res = xx^(n)
+;   IF check_math() GT 0 THEN BEGIN
+;       zero = where(abs(res) LT 1.e-10)
+;       IF zero[0] NE -1 THEN res[zero] = 0.0d
+;   END
+  res = 1.0d - temporary(res)
+  res = temporary(res)*(y2-y1)+y1
+;
+;  IF array_equal(sort(res), lindgen(n_elements(res)) ) NE 1 THEN stop
+  RETURN, res
+END
+
+FUNCTION spl_incr, x, y, x2, YP0 = yp0, YPN_1 = ypn_1
+;
+;---------------------------------
+; check and initialisation ...
+;---------------------------------
+  nx = n_elements(x)
+  ny = n_elements(y)
+  nx2 = n_elements(x2)
+; x must have at least 2 elements
+  IF nx LT 2 THEN stop 
+; y must have the same number of elements than x
+  IF nx NE ny THEN stop
+; x be monotonically increasing
+  IF min(x[1:nx-1]-x[0:nx-2]) LE 0 THEN stop 
+; x2 be monotonically increasing
+  IF N_ELEMENTS(X2) GE 2 THEN $
+  IF min(x2[1:nx2-1]-x2[0:nx2-2])  LE 0 THEN stop 
+; y be monotonically increasing
+  IF min(y[1:ny-1]-y[0:ny-2]) LT 0 THEN stop 
+;---------------------------------
+; first check: check if two consecutive values are equal
+;---------------------------------
+  bad = where(y[1:ny-1]-y[0:ny-2] EQ 0, cntbad)
+  IF cntbad NE 0 THEN BEGIN
+; define the results: y2
+      y2 = dblarr(nx2)
+; define xinx2: see help of value_locate
+;  if xinx2[i] eq -1   :                 x[bad[i]] <  x2[0]
+;  if xinx2[i] eq nx2-1:                 x[bad[i]] >= x2[nx2-1]
+;  else                : x2[xinx2[i]] <= x[bad[i]] <  x2[xinx2[i]+1]
+    xinx2 = value_locate(x2, x[bad])
+    xinx2_1 = value_locate(x2, x[bad+1])
+;
+; left side ... if there is x2 values smaller that x[bad[0]]. 
+; we force ypn_1 = 0.0d
+    IF xinx2[0] NE -1 THEN BEGIN
+      IF bad[0] EQ 0 THEN BEGIN
+        IF xinx2[0] NE 0 THEN stop
+        y2[0] = y[0]
+      ENDIF ELSE BEGIN 
+        y2[0:xinx2[0]] $
+          = spl_incr(x[0:bad[0]], y[0:bad[0]], x2[0:xinx2[0]] $
+                     , yp0 = yp0, ypn_1 = 0.0d)
+      ENDELSE 
+    ENDIF 
+; flat section
+    IF xinx2_1[0] NE -1 THEN $
+      y2[(xinx2[0]+1) < xinx2_1[0] : xinx2_1[0]] = y[bad[0]]
+; middle pieces ... if cntbad gt 1 then we have to cut spl_incr in
+; more than 2 pieces...
+      IF cntbad GT 1 THEN BEGIN
+; we take care of the piece located wetween bad[ib-1]+1 and bad[ib]
+        FOR ib = 1, cntbad-1 DO BEGIN
+; if there is x2 values smaller that x[bad[ib]], then the x2 values
+; located between bad[ib-1]+1 and bad[ib] are (xinx2_1[ib-1]+1:xinx2[ib]
+; and if we don't have two consecutive flat sections
+          IF xinx2[ib] NE -1 AND (bad[ib-1] NE bad[ib]-1) THEN begin
+            y2[(xinx2_1[ib-1]+1) < xinx2[ib]:xinx2[ib]] $
+              = spl_incr(x[bad[ib-1]+1:bad[ib]], y[bad[ib-1]+1:bad[ib]] $
+                         , x2[(xinx2_1[ib-1]+1) < xinx2[ib]:xinx2[ib]] $
+                         , yp0 = 0.0d, ypn_1 = 0.0d)
+          ENDIF
+; flat section
+          IF xinx2_1[ib] NE -1 THEN $
+            y2[(xinx2[ib]+1) < xinx2_1[ib] : xinx2_1[ib]] = y[bad[ib]]
+        ENDFOR
+      ENDIF
+; right side ... if there is x2 values larger that x[bad[cntbad-1]+1]. 
+; we force yp0 = 0.0d
+      IF xinx2_1[cntbad-1] NE nx2-1 THEN $
+        y2[xinx2_1[cntbad-1]+1:nx2-1] $
+        = spl_incr(x[bad[cntbad-1]+1:nx-1], y[bad[cntbad-1]+1:nx-1] $
+                        , x2[xinx2_1[cntbad-1]+1:nx2-1] $
+                        , yp0 = 0.0d, ypn_1 = ypn_1new)
+
+    RETURN, y2
+
+  ENDIF
+;-----------
+; compute the second derivative of the cubic spline on each x.
+;-----------
+  yscd = spl_init(x, y, yp0 = yp0, ypn_1 = ypn_1, /double)
+;---------------------------------
+; second check: none of the first derivative on x values must be negative.
+;---------------------------------
+;
+; compute the first derivative on x
+;
+  yifrst = spl_fstdrv(x, y, yscd, x)
+;
+; we force the negative first derivative to 0 by calling again
+; spl_incr with the keywords yp0 and ypn_1 to specify the
+; first derivative equal to 0
+;
+  bad = where(yifrst LT 0.0d, cntbad)
+  IF cntbad NE 0 THEN BEGIN
+;
+; we define the new values of the keyword ypn_1:
+; if the first derivative of the last value of x is negative 
+; we define the new values of the keyword ypn_1 to 0.0d0
+    IF bad[cntbad-1] EQ nx-1 THEN BEGIN 
+      ypn_1new = 0.0d
+; we remove this case from the list
+      IF cntbad GE 2 THEN bad = bad[0:cntbad-2]
+      cntbad = cntbad-1
+; else we take the value of ypn_1 if it was already defined
+    ENDIF ELSE IF n_elements(ypn_1) NE 0 THEN ypn_1new = ypn_1
+;
+; we define the new values of the keyword yp0:
+; if the first derivative of the first value of x is negative 
+; we define the new values of the keyword yp0 to 0.0
+    IF bad[0] EQ 0 THEN BEGIN
+      yp0new = 0.0d
+; we remove this case from the list
+      IF cntbad GE 2 THEN bad = bad[1:cntbad-1]
+      cntbad = cntbad-1
+; else we take the value of yp0 if it was already defined
+    ENDIF ELSE IF n_elements(yp0) NE 0 THEN yp0new = yp0
+;
+; if all the negative derivative corresponded to one of the cases above,
+; then we can directly call spl_incr with the new yp0new and ypn_1new
+    IF cntbad LE 0 THEN BEGIN
+      y2 = spl_incr(x, y, x2, yp0 = yp0new, ypn_1 = ypn_1new)
+;
+; else: there is still cases with negative derivative ...
+; we will cut spl_incr in n spl_incr and specify yp0, ypn_1
+; for each of this n spl_incr  
+    ENDIF ELSE BEGIN
+; define xinx2: see help of value_locate
+;  if xinx2[i] eq -1   :                 x[bad[i]] <  x2[0]
+;  if xinx2[i] eq nx2-1:                 x[bad[i]] >= x2[nx2-1]
+;  else                : x2[xinx2[i]] <= x[bad[i]] <  x2[xinx2[i]+1]
+      xinx2 = value_locate(x2, x[bad])
+      y2 = dblarr(nx2)
+; left side ... if there is x2 values smaller that x[bad[0]]. 
+; we force ypn_1 = 0.0d
+      IF xinx2[0] NE -1 THEN $
+        y2[0:xinx2[0]] $
+        = spl_incr(x[0:bad[0]], y[0:bad[0]], x2[0:xinx2[0]] $
+                        , yp0 = yp0new, ypn_1 = 0.0d)
+; middle pieces ... if cntbad gt 1 then we have to cut spl_incr in
+; more than 2 pieces -> we have middle pieces for which 
+; we force yp0 = 0.0d and ypn_1 = 0.0d
+      IF cntbad GT 1 THEN BEGIN
+; we take care of the piece located wetween bad[ib-1] and bad[ib]
+        FOR ib = 1, cntbad-1 DO BEGIN
+; if there is x2 values smaller that x[bad[ib]], then the x2 values
+; located between bad[ib-1] and bad[ib] are (xinx2[ib-1]+1:xinx2[ib]
+          IF xinx2[ib] NE -1 THEN begin
+            y2[(xinx2[ib-1]+1) < xinx2[ib]:xinx2[ib]] $
+              = spl_incr(x[bad[ib-1]:bad[ib]], y[bad[ib-1]:bad[ib]] $
+                              , x2[(xinx2[ib-1]+1) < xinx2[ib]:xinx2[ib]] $
+                              , yp0 = 0.0d, ypn_1 = 0.0d)
+          endif
+        ENDFOR
+      ENDIF
+; right side ... if there is x2 values larger that x[bad[cntbad-1]]. 
+; we force yp0 = 0.0d
+      IF xinx2[cntbad-1] NE nx2-1 THEN $
+        y2[xinx2[cntbad-1]+1:nx2-1] $
+        = spl_incr(x[bad[cntbad-1]:nx-1], y[bad[cntbad-1]:nx-1] $
+                        , x2[xinx2[cntbad-1]+1:nx2-1] $
+                        , yp0 = 0.0d, ypn_1 = ypn_1new)
+    ENDELSE 
+; we return the checked and corrected value of yfrst
+;       FOR i = 0, nx-1 DO BEGIN
+;         same = where(abs(x2- x[i]) LT 1.e-10, cnt)
+; ;        IF cnt NE 0 THEN y2[same] = y[i] 
+;       ENDFOR
+    RETURN, y2
+  ENDIF
+;
+; we can be in this part of the code only if:
+;  (1) spl_incr is called by itself
+;  (2) none are the first derivative in x are negative (because they have been 
+;      checked and corrected by the previous call to spl_incr, see above)  
+;---------------------------------
+; third check: we have to make sure that the first derivative cannot
+;               have negative values between on x[0] and x[nx-1]
+;---------------------------------
+;
+; first we compute the first derivative, next we correct the values
+; where we know that the first derivative can be negative. 
+;
+  y2 = spl_interp(x, y, yscd, x2, /double)
+;
+; between x[i] and x[i+1], the cubic spline is a cubic function:
+; y  =  a*X^3 +  b*X^2 + c*X + d
+; y' = 3a*X^2 + 2b*X   + c
+; y''= 6a*X   + 2b
+; if we take X = x[i+1]-x[i] then
+;    d = y[i]; c = y'[i]; b = 0.5 * y''[i], 
+;    a = 1/6 * (y''[i+1]-y''[i])/(x[i+1]-x[i])
+; 
+; y'[i] and y'[i+1] are positive so y' can be negative
+; between x[i] and x[i+1] only if 
+;   1) a > 0
+;            ==> y''[i+1] > y''[i]
+;   2) y' reach its minimum value between  x[i] and x[i+1] 
+;      -> 0 < - b/(3a) < x[i+1]-x[i] 
+;            ==> y''[i+1] > 0 > y''[i]
+;
+; we do a first selection by looking for those points...
+;
+  loc = lindgen(nx-1)
+  maybebad = where(yscd[loc] LE 0.0d AND yscd[loc+1] GE 0.0d, cntbad)
+;
+  IF cntbad NE 0 THEN BEGIN
+
+    mbbloc = loc[maybebad]
+
+    aaa = (yscd[mbbloc+1]-yscd[mbbloc])/(6.0d*(x[mbbloc+1]-x[mbbloc]))
+    bbb = 0.5d * yscd[mbbloc]
+    ccc = yifrst[mbbloc]
+    ddd = y[mbbloc]
+;
+; definitive selection:
+; y' can become negative if and only if (2b)^2 - 4(3a)c > 0
+; y' can become negative if and only if    b^2  - (3a)c > 0
+;
+    delta = bbb*bbb - 3.0d*aaa*ccc
+;
+    bad = where(delta GT 0, cntbad)
+;
+    IF cntbad NE 0 THEN BEGIN
+      delta = delta[bad]
+      aaa = aaa[bad]
+      bbb = bbb[bad]
+      ccc = ccc[bad]
+      ddd = ddd[bad]
+      bad = maybebad[bad]
+; define xinx2_1: see help of value_locate
+;  if xinx2_1[i] eq -1   :                   x[bad[i]] <  x2[0]
+;  if xinx2_1[i] eq nx2-1:                   x[bad[i]] >= x2[nx2-1]
+;  else                  : x2[xinx2_1[i]] <= x[bad[i]] <  x2[xinx2_1[i]+1]
+      xinx2_1 = value_locate(x2, x[bad])
+; define xinx2_2: see help of value_locate
+;  if xinx2_2[i] eq -1   :                   x[bad[i]+1] <  x2[0]
+;  if xinx2_2[i] eq nx2-1:                   x[bad[i]+1] >= x2[nx2-1]
+;  else                  : x2[xinx2_2[i]] <= x[bad[i]+1] <  x2[xinx2_2[i]+1]
+      xinx2_2 = value_locate(x2, x[bad+1])
+; to avoid the particular case when x2 = x[bad[i]]
+; and there is no other x2 point until x[bad[i]+1]
+      xinx2_1 = xinx2_1 < (xinx2_2-1)
+;
+      FOR ib = 0, cntbad-1 DO BEGIN
+;
+; at least one of the x2 points must be located between
+; x[bad[ib]] and x[bad[ib]+1]
+        IF x2[0] LE x[bad[ib]+1] AND x2[nx2-1] GE x[bad[ib]] THEN BEGIN
+;
+          CASE 1 OF
+            yifrst[bad[ib]+1] EQ 0.0d:BEGIN
+; case pur convex: we use the first derivative of 1-(1-x)^n
+; and ajust n to get the good value: yifrst[bad[ib]] in x[bad[ib]]
+              y2[xinx2_1[ib]+1:xinx2_2[ib]] $
+                = pure_convex(x[bad[ib]], x[bad[ib]+1] $
+                              , y[bad[ib]], y[bad[ib]+1]  $
+                              , yifrst[bad[ib]] $
+                              , x2[xinx2_1[ib]+1:xinx2_2[ib]])
+            END
+            yifrst[bad[ib]] EQ 0.0d:BEGIN
+; case pur concave: we use the first derivative of x^n
+; and ajust n to get the good value: yifrst[bad[ib]+1] in x[bad[ib]+1]
+              y2[xinx2_1[ib]+1:xinx2_2[ib]] $
+                = pure_concave(x[bad[ib]], x[bad[ib]+1] $
+                               , y[bad[ib]], y[bad[ib]+1] $
+                               , yifrst[bad[ib]+1] $
+                               , x2[xinx2_1[ib]+1:xinx2_2[ib]])
+            END
+            ELSE:BEGIN
+; in those cases, the first derivative has 2 zero between
+; x[bad[ib]] and x[bad[ib]+1]. We look for the minimum value of the
+; first derivative that correspond to the inflection point of y              
+              xinfl = -bbb[ib]/(3.0d*aaa[ib])
+; we compute the y value for xinfl
+              yinfl = aaa[ib]*xinfl*xinfl*xinfl + bbb[ib]*xinfl*xinfl $
+                + ccc[ib]*xinfl + ddd[ib]
+;                
+              CASE 1 OF
+; if y[xinfl] smaller than y[bad[ib]] then we conserve y2 until
+; the first zero of y2 and from this point we use x^n and ajust n to
+; get the good value: yifrst[bad[ib]+1] in x[bad[ib]+1]
+                yinfl LT y[bad[ib]]:BEGIN
+; value of the first zero (y'[xzero]=0)
+                  xzero = (-bbb[ib]-sqrt(delta[ib]))/(3.0d*aaa[ib])
+; value of y[xzero]...
+                  yzero = aaa[ib]*xzero*xzero*xzero + bbb[ib]*xzero*xzero $
+                    + ccc[ib]*xzero + ddd[ib]
+; if yzero > y[bad[ib]+1] then we cannot applay the method we want to
+; apply => we use then convex-concave case by changing by hand the
+; value of yinfl and xinfl
+                  IF yzero GT y[bad[ib]+1] THEN BEGIN
+                    yinfl = 0.5d*(y[bad[ib]+1]+y[bad[ib]])
+                    xinfl = 0.5d*(x[bad[ib]+1]-x[bad[ib]])
+                    GOTO, convexconcave
+                  ENDIF
+; define xinx2_3: see help of value_locate
+;  if xinx2_3[ib] eq -1   :                x[bad[ib]]+xzero <  x2[0]
+;  if xinx2_3[ib] eq nx2-1:                x[bad[ib]]+xzero >= x2[nx2-1]
+;  else                   : x2[xinx2_3] <= x[bad[ib]]+xzero <  x2[xinx3_2+1]
+                  xinx2_3 = value_locate(x2, x[bad[ib]]+xzero)
+; to avoid the particular case when x2 = x[bad[ib]]+xzero
+; and there is no other x2 point until x[bad[ib]+1]
+                  xinx2_3 = xinx2_3 < (xinx2_2[ib]-1)
+                  IF xinx2_2[ib] GE xinx2_3+1 THEN BEGIN
+                    y2[xinx2_3+1:xinx2_2[ib]] $
+                      = pure_concave(x[bad[ib]]+xzero, x[bad[ib]+1] $
+                                     , yzero, y[bad[ib]+1] $
+                                     , yifrst[bad[ib]+1] $
+                                     , x2[xinx2_3+1:xinx2_2[ib]])
+                  ENDIF                
+                END
+; if y[xinfl] bigger than y[bad[ib]+1] then we conserve y2 from
+; the second zero of y2 and before this point we use 1-(1-x)^n and
+; ajust n to get the good value: yifrst[bad[ib]] in x[bad[ib]]
+                yinfl GT y[bad[ib]+1]:BEGIN
+; value of the second zero (y'[xzero]=0)
+                  xzero = (-bbb[ib]+sqrt(delta[ib]))/(3.0d*aaa[ib])
+; value of y[xzero]...
+                  yzero = aaa[ib]*xzero*xzero*xzero + bbb[ib]*xzero*xzero $
+                    + ccc[ib]*xzero + ddd[ib]
+; if yzero < y[bad[ib]] then we cannot applay the method we want to
+; apply => we use then convex-concave case by changing by hand the
+; value of yinfl and xinfl
+                  IF yzero lt y[bad[ib]] THEN BEGIN
+                    yinfl = 0.5d*(y[bad[ib]+1]+y[bad[ib]])
+                    xinfl = 0.5d*(x[bad[ib]+1]-x[bad[ib]])
+                    GOTO, convexconcave
+                  ENDIF
+; define xinx2_3: see help of value_locate
+;  if xinx2_3[ib] eq -1   :                x[bad[ib]]+xzero <  x2[0]
+;  if xinx2_3[ib] eq nx2-1:                x[bad[ib]]+xzero >= x2[nx2-1]
+;  else                   : x2[xinx2_3] <= x[bad[ib]]+xzero <  x2[xinx3_2+1]
+                  xinx2_3 = value_locate(x2, x[bad[ib]]+xzero)
+                  IF xinx2_3 ge xinx2_1[ib]+1 THEN BEGIN
+                    y2[xinx2_1[ib]+1:xinx2_3] $
+                      = pure_convex(x[bad[ib]], x[bad[ib]]+xzero  $
+                                    , y[bad[ib]], yzero   $
+                                    , yifrst[bad[ib]] $
+                                    , x2[xinx2_1[ib]+1:xinx2_3])
+                  ENDIF                
+                END
+                ELSE:BEGIN
+convexconcave:
+; define xinx2_3: see help of value_locate
+;  if xinx2_3[ib] eq -1   :                x[bad[ib]]+xzero <  x2[0]
+;  if xinx2_3[ib] eq nx2-1:                x[bad[ib]]+xzero >= x2[nx2-1]
+;  else                   : x2[xinx2_3] <= x[bad[ib]]+xzero <  x2[xinx3_2+1]
+                  xinx2_3 = value_locate(x2, x[bad[ib]]+xinfl)
+
+                  IF xinx2_3 ge xinx2_1[ib]+1 THEN BEGIN
+                    y2[xinx2_1[ib]+1:xinx2_3] $
+                      = pure_convex(x[bad[ib]], x[bad[ib]]+xinfl  $
+                                    , y[bad[ib]], yinfl  $
+                                    , yifrst[bad[ib]] $
+                                    , x2[xinx2_1[ib]+1:xinx2_3])
+                    
+                  ENDIF                
+                  IF xinx2_2[ib] GE xinx2_3+1 THEN BEGIN
+                    y2[xinx2_3+1:xinx2_2[ib]] $
+                      = pure_concave(x[bad[ib]]+xinfl, x[bad[ib]+1] $
+                                     , yinfl, y[bad[ib]+1] $
+                                     , yifrst[bad[ib]+1] $
+                                     , x2[xinx2_3+1:xinx2_2[ib]])
+                  ENDIF                
+                END
+              ENDCASE
+
+            END 
+          ENDCASE
+        ENDIF
+      ENDFOR
+
+    ENDIF
+  ENDIF
+;
+  RETURN, y2
+;
+;------------------------------------------------------------------
+;------------------------------------------------------------------
+;
+END
Index: trunk/Interpolation/spl_keep_mean.pro
===================================================================
--- trunk/Interpolation/spl_keep_mean.pro	(revision 69)
+++ trunk/Interpolation/spl_keep_mean.pro	(revision 69)
@@ -0,0 +1,144 @@
+;------------------------------------------------------------
+;------------------------------------------------------------
+;------------------------------------------------------------
+;+
+; NAME:spl_keep_mean
+;
+; PURPOSE:
+;
+; Given the arrays X and Y, which tabulate a function (with the X[i]
+; AND Y[i] in ascending order), and given an input value X2, the
+; SPL_INCR function returns an interpolated value for the given values
+; of X2. The interpolation method is based on cubic spline, corrected
+; in a way that integral of the interpolated values is the same as the
+; integral of the input values. (-> for exemple to build daily data
+; from monthly mean and keep the monthly mean of the computed daily
+; data equa to the original values)
+;
+; CATEGORY:
+;
+; CALLING SEQUENCE: y2 =  spl_keep_mean(x, y, x2)
+;
+; INPUTS:
+;
+;    x: An n-element (at least 2) input vector that specifies the
+;    tabulate points in a strict ascending order.
+;
+;    y: an array with one element less than x. y[i] represents the
+;    mean value between x[i] and x[i+1]. if /GE0 is activated, y must
+;    have positive values.
+;
+;    x2: The input values for which the interpolated values are
+;    desired. Its values must be strictly monotonically increasing. 
+;
+; KEYWORD PARAMETERS:
+;
+;    /GE0: to force that y2 is always GE than 0. In that case, y must
+;    also be GE than 0.
+;
+;    YP0: The first derivative of the interpolating function at the
+;    point X0. If YP0 is omitted, the second derivative at the
+;    boundary is set to zero, resulting in a "natural spline."
+;
+;    YPN_1: The first derivative of the interpolating function at the
+;    point Xn-1. If YPN_1 is omitted, the second derivative at the
+;    boundary is set to zero, resulting in a "natural spline." 
+;
+; OUTPUTS: 
+;
+;    y2: the meean value between two consecutive values of x2. This
+;    array has one element less than y2. y2 has double precision.
+;
+; COMMON BLOCKS: none
+;
+; SIDE EFFECTS: ?
+;
+; RESTRICTIONS:
+;   It might be possible that y2 has very small negative values
+;   (amplitude smaller than 1.e-6)... 
+;
+;
+; EXAMPLE:
+;
+;    12 monthly values of precipitations into daily values:
+;
+;    yr1 = 1990
+;    yr2 = 1992
+;    nyr = yr2-yr1+1
+;    n1 = 12*nyr+1
+;    x = julday(1+findgen(n1), replicate(1, n1) $
+;           , replicate(yr1, n1), fltarr(n1))
+;    n2 = 365*nyr + total(leapyr(yr1+indgen(nyr))) + 1
+;    x2 = julday(replicate(1, n2), 1+findgen(n2) $
+;               , replicate(yr1, n2), fltarr(n2))
+;    y = abs(randomn(0, n1-1))
+;    y2 = spl_keep_mean(x, y, x2, /ge0)
+
+;    print, min(x, max = ma), ma
+;    print, min(x2, max = ma), ma
+;    print, vairdate([min(x, max = ma), ma])
+;    print, total(y*(x[1:n1-1]-x[0:n1-2]))
+;    print, total(y2*(x2[1:n2-1]-x2[0:n2-2]))
+;
+; MODIFICATION HISTORY:
+;  Sebastien Masson (smasson@lodyc.jussieu.fr): May 2005
+;-
+;------------------------------------------------------------
+;------------------------------------------------------------
+;------------------------------------------------------------
+FUNCTION spl_keep_mean, x, yin, x2, YP0 = yp0, YPN_1 = ypn_1, GE0 = ge0
+;
+;---------------------------------
+; check and initialisation ...
+;---------------------------------
+;
+  nx = n_elements(x)
+  ny = n_elements(yin)
+  nx2 = n_elements(x2)
+; x must have at least 2 elements
+  IF nx LT 2 THEN stop 
+; x2 must have at least 2 elements
+  IF nx2 LT 2 THEN stop 
+; x be monotonically increasing
+  IF min(x[1:nx-1]-x[0:nx-2]) LE 0 THEN stop 
+; x2 be monotonically increasing
+  IF min(x2[1:nx2-1]-x2[0:nx2-2])  LE 0 THEN stop 
+;
+;---------------------------------
+; compute the integral of y
+;---------------------------------
+; if spl_keep_mean is called by the user (and not by itself) we must compute
+; the integral of y. yin must have one element less than x
+  IF nx NE ny+1 THEN stop
+  y = double(yin)*double(x[1:nx-1]-x[0:nx-2])
+  y = [0.0d, temporary(y)]
+  y = total(temporary(y), /cumulative, /double)
+;
+;---------------------------------
+; compute the "spline" interpolation
+;---------------------------------
+;
+  IF keyword_set(ge0) THEN BEGIN
+; if the want that the interpolated values are always >= 0, we must
+; have yin >= 0.0d
+    IF min(yin) LT 0 THEN stop
+; call spl_incr
+    y2 = spl_incr(x, temporary(y), x2, yp0 = yp0, ypn_1 = ypn_1)
+  ENDIF ELSE BEGIN
+    yscd = spl_init(x, y, yp0 = yp0, ypn_1 = ypn_1, /double)
+    y2 = spl_interp(x, y, temporary(yscd), x2, /double)
+  ENDELSE
+;                ;
+;---------------------------------
+; Compute the derivative of y
+;---------------------------------
+;
+  yfrst = (y2[1:nx2-1]-y2[0:nx2-2])/(x2[1:nx2-1]-x2[0:nx2-2])
+; it can happen that we have very small negative values (-1.e-6 for ex)
+;  yfrst = 0.0d > temporary(yfrst)
+  RETURN, yfrst
+  
+;------------------------------------------------------------------
+;------------------------------------------------------------------
+;
+ END
