source: trunk/SRC/Interpolation/compute_fromreg_bilinear_weigaddr.pro @ 327

;+
;
; @file_comments
; compute the weight and address needed to interpolate data from a
; "regular grid" to any grid using the bilinear method
;
; @categories
; Interpolation
;
; @param alonin{in}{required}{type=2d array}
; longitude of the input data
;
; @param alatin {in}{required}{type=2d array}
; latitude of the input data
;
; @param olonin {in}{required}{type=2d array}
; longitude of the output data
;
; @param olat {in}{required}{type=2d array}
; latitude of the output data
;
; @keyword NONORTHERNLINE {type=scalar 0 or 1}{default=0}
; put 1 if you don't want to take into
; account the northern line of the input data when performing the interpolation.
;
; @keyword NOSOUTHERNLINE {type=scalar 0 or 1}{default=0}
; put 1 if you don't want to take into
; account the southern line of the input data when performing the interpolation.
;
; @param weig {out}{type=2d array}
; (see ADDR)
;
; @param addr {out}{type=2d array}
; 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 grid", defined as a grid for which each
;  longitude lines have the same latitude and each latitude columns have the
;  same longitude.
;  - We supposed the data are located on a sphere, with a periodicity along
;  the longitude.
;  - points located out of the southern and northern boundaries are interpolated
;  using a linear interpolation only along the longitudinal direction.
;
; @history
;  November 2005: Sebastien Masson (smasson\@lodyc.jussieu.fr)
;
; @version
; $Id$
;
;-
PRO compute_fromreg_bilinear_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
; for longitude periodic boundary condition we add the fist
; column on the right side of the array and
  alon = [alon, alon[0]+360.]
  jpia = jpia+1L
; 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
;
  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 and 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
;
; we want to do bilinear interpolation => for each ocean point, we must
; find in which atm cell it is located.
; if the ocean point is out of the atm grid, we use closest neighbor
; interpolation
;
; for each T point of oce grid, we find in which atmospheric cell it is
; located.
; As the atmospheric grid is regular, we can use inrecgrid instead
; of inquad.
  pos = inrecgrid(olon, olat, alon[0:jpia-2L], alat[0:jpja-2L] $
                  , checkout = [alon[jpia-1L], alat[jpja-1L]], /output2d)
; checks...
; for longitude, each ocean point must be located in atm cell.
  IF (where(pos[0, *] EQ -1))[0] NE -1 THEN stop
; no ocean point should be located westward of the left boundary of the
; atm cell in which it is supposed to be located
  IF total(olon LT alon[pos[0, *]]) NE 0 THEN stop
; no ocean point should be located eastward of the right boundary of the
; atm cell in which it is supposed to be located
  IF total(olon GT alon[pos[0, *]+1]) NE 0 THEN stop
;
; we use bilinear interpolation
;
; we change the coordinates of each ocean point to fit into a
; rectangle defined by:
;
;  y2 *------------*
;     |            |
;     |            |
;     |            |
;  y1 *------------*
;     x1          x2
;
;    X = (x-x1)/(x2-x1)
;    Y = (y-y1)/(y2-y1)
;
  indx = pos[0, *]
  indy = (temporary(pos))[1, *]
; points located out of the atmospheric grid...(too much northward or southward)
  bad = where(indy EQ -1)
  indy = 0 > indy
;
  IF max(indx) GT jpia-2 THEN stop ; checks...
  IF max(indy) GT jpja-2 THEN stop ; checks...
; x coordinates of the atm cell
  x1 = alon[indx]
  x2 = alon[indx+1]
; new x coordinates of the ocean points in each cell
  divi = temporary(x2)-x1
  glamnew = (olon-x1)/temporary(divi)
  x1 = -1 ; free memory
  olon = -1 ; free memory
; y coordinates of the atm cell
  y1 = alat[indy]
  y2 = alat[indy+1]             ;
; new y coordinates of the ocean points in each cell
  divi = temporary(y2)-y1
  zero = where(divi EQ 0)
  IF zero[0] NE -1 THEN divi[zero] = 1.
  gphinew = (olat-y1)/temporary(divi)
  y1 = -1 ; free memory
; checks...
  IF min(glamnew) LT 0 THEN stop
  IF max(glamnew) GT 1 THEN stop
;
; weight and address array used for bilinear interpolation.
  xaddr = lonarr(4, jpio*jpjo)
  xaddr[0, *] = indx
  xaddr[1, *] = indx + 1L
  xaddr[2, *] = indx + 1L
  xaddr[3, *] = indx
;
  yaddr = lonarr(4, jpio*jpjo)
  yaddr[0, *] = indy
  yaddr[1, *] = indy
  yaddr[2, *] = indy + 1L
  yaddr[3, *] = indy + 1L
; compute the weight for the bilinear interpolation.
  weig = fltarr(4, jpio*jpjo)
  weig[0, *] = (1.-glamnew) * (1.-gphinew)
  weig[1, *] =     glamnew  * (1.-gphinew)
  weig[2, *] =     glamnew  *     gphinew
  weig[3, *] = (1.-glamnew) *     gphinew
; free memory
  gphinew = -1
  IF bad[0] EQ -1 THEN glamnew = -1 ELSE glamnew = (temporary(glamnew))[bad]
; we work now on the "bad" points
; linear interpolation only along the longitudinal direction
  IF bad[0] NE -1 THEN BEGIN
    ybad = olat[bad]
; the ocean points that are not located into an atm cell should be
; located northward of the northern boundary of the atm grid
;      or southward of the southern boundary of the atm grid
    IF total(ybad GE min(alat) AND ybad LE max(alat)) GE 1 THEN stop
;
    weig[0, bad] = (1.-glamnew)
    weig[1, bad] = temporary(glamnew)
    weig[2, bad] = 0.
    weig[3, bad] = 0.
    south = where(ybad LT alat[0])
    IF south[0] NE -1 THEN yaddr[*, bad[temporary(south)]] = 0L
    north = where(ybad GT alat[jpja-1])
    IF north[0] NE -1 THEN yaddr[*, bad[temporary(north)]] = jpja-1
    ybad = -1 & bad = -1 ; free memory
  ENDIF
; check totalweight = 1
  totalweig = abs(1.d - total(weig, 1, /double))
  IF (where(temporary(totalweig) GE 1.e-5))[0] NE -1 THEN stop
;
; come back to the original atm grid without longitudinal overlap.
;
  jpia = jpia-1L
  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