source: trunk/SRC/Interpolation/compute_fromreg_imoms3_weigaddr.pro @ 157

;+
;
; @file_comments
; compute the weight and address neede to interpolate data from a
; "regular grid" to any grid using the imoms3 method
;
; @categories 
; Interpolation
;
; @param alonin {in}{required}
; longitude of the input data
;
; @param alatin {in}{required}
; latitude of the input data
;
; @param olonin {in}{required}
; longitude of the output data
; @param olat {in}{required}
; latitude of the output data
;
; @keyword NONORTHERNLINE
; @keyword NOSOUTHERNLINE
; activate if you don't want to take into account the northen/southern line
; of the input data when perfoming the interpolation.
;
; @param weig {out}
; @param addr {out}
; 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)
;
; @restrictions
;  - 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.
;  -  points located between the first/last 2 lines are interpolated
;  using a imoms3 interpolation along the longitudinal direction and linear
;  interpolation along the latitudinal direction
;  - points located out of the southern and northern boundaries are interpolated
;  using a imoms3 interpolation only along the longitudinal direction.
;
; @history
;  November 2005: Sebastien Masson (smasson\@lodyc.jussieu.fr)
;  March 2006: works for rectangular grids
;
; @version $Id$
;
;-
;
;----------------------------------------------------------
;----------------------------------------------------------
;
PRO compute_fromreg_imoms3_weigaddr, alonin, alatin, olonin, olat, weig, addr $
                                     , NONORTHERNLINE = nonorthernline, NOSOUTHERNLINE = nosouthernline
;
  compile_opt idl2, strictarrsubs
;
  alon = alonin
  alat = alatin
  olon = olonin
;
  jpia = n_elements(alon)
  jpja = n_elements(alat)
;
  jpio = (size(olon, /dimensions))[0]
  jpjo = (size(olon, /dimensions))[1]
;
; alon
  minalon = min(alon,  max = maxalon)
  IF maxalon-minalon GE 360. THEN stop
; alon must be monotonically increasing
  IF array_equal(sort(alon), lindgen(jpia)) NE 1 THEN BEGIN
    shiftx = -(where(alon EQ min(alon)))[0]
    alon = shift(alon, shiftx)
    IF array_equal(sort(alon), lindgen(jpia)) NE 1 THEN stop
  ENDIF ELSE shiftx = 0
; 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 noregx = 1
; we extend the longitude range of alon (-> easy interpolation even
; near minalon et maxalon)
  toadd = 10*jpia/360+1
  alon = [alon[jpia-toadd:jpia-1]-360., alon[*], alon[0:toadd-1]+360.]
  jpia = jpia+2*toadd
; alat
  revy = alat[0] GT alat[1]
  IF revy THEN alat = reverse(alat)
; alat must be monotonically increasing
  IF array_equal(sort(alat), lindgen(jpja)) NE 1 THEN stop
; 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 THEN noregy = 1
;
  if keyword_set(nonorthernline) then BEGIN
    jpja = jpja - 1L
    alat = alat[0: jpja-1L]
  ENDIF
  if keyword_set(nosouthernline) then BEGIN
    alat = alat[1: jpja-1L]
    jpja = jpja - 1L
  ENDIF
; olon between minalon et minalon+360
  out = where(olon LT minalon)
  WHILE out[0] NE -1 DO BEGIN
    olon[out] = olon[out]+360.
    out = where(olon LT minalon)
  ENDWHILE
  out = where(olon GE minalon+360.)
  WHILE out[0] NE -1 DO BEGIN
    olon[out] = olon[out]- 360.
    out = where(olon GE minalon+360.)
  ENDWHILE
; make sure that all values of olon are located within values of alon
  IF min(olon, max = ma) LT minalon THEN stop
  IF ma GE minalon+360. THEN stop
;
  xaddr = lonarr(16, jpio*jpjo)
  yaddr = lonarr(16, jpio*jpjo)
  weig = fltarr(16, jpio*jpjo)
;
  indexlon = value_locate(alon, olon)
  IF total(alon[indexlon] GT olon) NE 0 THEN stop
  IF total(alon[indexlon + 1L] LE olon) NE 0 THEN stop
  IF (where(indexlon LE 1L     ))[0] NE -1 THEN stop
  IF (where(indexlon GE jpia-3L))[0] NE -1 THEN stop
  indexlat = value_locate(alat, olat)
;
; for the ocean points located below the atm line
; jpja-2 and above the line 1
; for those points we can always find 16 neighbors
; imoms interpolation along longitude and latitude
;
  short = where(indexlat LT jpja-2L AND indexlat GE 1L)
  ilon = indexlon[short]
  ilat = indexlat[short]
;
  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)
  IF abs(mi-1) GE 1.e-6 THEN stop
  IF abs(ma-1) GE 1.e-6 THEN stop
;
  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)
  IF abs(mi-1) GE 1.e-6 THEN stop
  IF abs(ma-1) GE 1.e-6 THEN stop
;
; line 0
  xaddr[0, short] = ilon - 1L
  xaddr[1, short] = ilon
  xaddr[2, short] = ilon + 1L
  xaddr[3, short] = ilon + 2L
  yaddr[0, short] = ilat - 1L
  yaddr[1, short] = yaddr[0, short]
  yaddr[2, short] = yaddr[0, short]
  yaddr[3, short] = yaddr[0, short]
  weig[0, short] = wx0 * wy0
  weig[1, short] = wx1 * wy0
  weig[2, short] = wx2 * wy0
  weig[3, short] = wx3 * wy0
; line 1
  xaddr[4, short] = ilon - 1L
  xaddr[5, short] = ilon
  xaddr[6, short] = ilon + 1L
  xaddr[7, short] = ilon + 2L
  yaddr[4, short] = ilat
  yaddr[5, short] = ilat
  yaddr[6, short] = ilat
  yaddr[7, short] = ilat
  weig[4, short] = wx0 * wy1
  weig[5, short] = wx1 * wy1
  weig[6, short] = wx2 * wy1
  weig[7, short] = wx3 * wy1
; line 2
  xaddr[8, short] = ilon - 1L
  xaddr[9, short] = ilon
  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
  weig[10, short] = wx2 * wy2
  weig[11, short] = wx3 * wy2
; line 3
  xaddr[12, short] = ilon - 1L
  xaddr[13, short] = ilon
  xaddr[14, short] = ilon + 1L
  xaddr[15, short] = ilon + 2L
  yaddr[12, short] = ilat + 2L
  yaddr[13, short] = yaddr[12, short]
  yaddr[14, short] = yaddr[12, short]
  yaddr[15, short] = yaddr[12, short]
  weig[12, short] = wx0 * wy3
  weig[13, short] = wx1 * wy3
  weig[14, short] = wx2 * wy3
  weig[15, short] = wx3 * wy3
;
  mi = min(total(weig[*, short], 1), max = ma)
  IF abs(mi-1) GE 1.e-6 THEN stop
  IF abs(ma-1) GE 1.e-6 THEN stop
;
; for the ocean points located between the atm lines
; jpja-2 and jpja-1 or between the atm lines 0 and 1
; linear interpolation between line 1 and line 2
;
  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+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+1L]-olat[short])/delta
    IF min(d2, max = ma) LE 0 THEN stop
    IF ma GT 1 THEN stop
    wy2 = 1.- temporary(d2)
;
    mi = min(wy1+wy2, max = ma)
    IF abs(mi-1) GE 1.e-6 THEN stop
    IF abs(ma-1) GE 1.e-6 THEN stop
; but imoms3 along the longitude
    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)
    IF abs(mi-1) GE 1.e-6 THEN stop
    IF abs(ma-1) GE 1.e-6 THEN stop
; line 1
    xaddr[0, short] = ilon - 1L
    xaddr[1, short] = ilon
    xaddr[2, short] = ilon + 1L
    xaddr[3, short] = ilon + 2L
    yaddr[0, short] = ilat
    yaddr[1, short] = ilat
    yaddr[2, short] = ilat
    yaddr[3, short] = ilat
    weig[0, short] = wx0 * wy1
    weig[1, short] = wx1 * wy1
    weig[2, short] = wx2 * wy1
    weig[3, short] = wx3 * wy1
; line 2
    xaddr[4, short] = ilon - 1L
    xaddr[5, short] = ilon
    xaddr[6, short] = ilon + 1L
    xaddr[7, short] = ilon + 2L
    yaddr[4, short] = ilat + 1L
    yaddr[5, short] = yaddr[4, short]
    yaddr[6, short] = yaddr[4, short]
    yaddr[7, short] = yaddr[4, short]
    weig[4, short] = wx0 * wy2
    weig[5, short] = wx1 * wy2
    weig[6, short] = wx2 * wy2
    weig[7, short] = wx3 * wy2
;
    mi = min(total(weig[*, short], 1), max = ma)
    IF abs(mi-1) GE 1.e-6 THEN stop
    IF abs(ma-1) GE 1.e-6 THEN stop
;
  ENDIF
;
; for the ocean points located below the line 0
; Interpolation only along the longitude
;
  short = where(indexlat EQ -1)
  IF short[0] NE -1 THEN BEGIN
    ilon = indexlon[short]
;
    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)
    IF abs(mi-1) GE 1.e-6 THEN stop
    IF abs(ma-1) GE 1.e-6 THEN stop
; line 1
    xaddr[0, short] = ilon - 1L
    xaddr[1, short] = ilon
    xaddr[2, short] = ilon + 1L
    xaddr[3, short] = ilon + 2L
    yaddr[0:3, short] = 0
    weig[0, short] = wx0
    weig[1, short] = wx1
    weig[2, short] = wx2
    weig[3, short] = wx3
;
    mi = min(total(weig[*, short], 1), max = ma)
    IF abs(mi-1) GE 1.e-6 THEN stop
    IF abs(ma-1) GE 1.e-6 THEN stop
;
  ENDIF
;
; for the ocean points located above jpia-1
; Interpolation only along the longitude
;
  short = where(indexlat EQ jpja-1L)
  IF short[0] NE -1 THEN BEGIN
    ilon = indexlon[short]
;
    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)
    IF abs(mi-1) GE 1.e-6 THEN stop
    IF abs(ma-1) GE 1.e-6 THEN stop
; line 1
    xaddr[0, short] = ilon-1L
    xaddr[1, short] = ilon
    xaddr[2, short] = ilon+1L
    xaddr[3, short] = ilon+2L
    yaddr[0:3, short] = jpja-1L
    weig[0, short] = wx0
    weig[1, short] = wx1
    weig[2, short] = wx2
    weig[3, short] = wx3
;
    mi = min(total(weig[*, short], 1), max = ma)
    IF abs(mi-1) GE 1.e-6 THEN stop
    IF abs(ma-1) GE 1.e-6 THEN stop
;
  ENDIF
;
; Come back to the original index of atm grid without longitudinal overlap.
;
;
  xaddr = temporary(xaddr) - toadd
  jpia = jpia - 2*toadd
; make sure all values are ge 0
  xaddr = temporary(xaddr) + jpia
; range the values between 0 and jpia-1
  xaddr = temporary(xaddr) mod jpia
;
; take into account shiftx if needed
  IF shiftx NE 0 THEN xaddr = (temporary(xaddr) - shiftx) MOD jpia
; take into account nosouthernline and nonorthernline
  if keyword_set(nosouthernline) then BEGIN
    yaddr = temporary(yaddr) + 1L
    jpja = jpja + 1L
  ENDIF
  if keyword_set(nonorthernline) then jpja = jpja + 1L
; take into account revy if needed
  IF revy EQ 1 THEN yaddr = jpja - 1L - temporary(yaddr)
;                        ;
  addr = temporary(yaddr)*jpia+temporary(xaddr)
;
  RETURN
END