source: trunk/Interpolation/compute_fromreg_bilinear_weigaddr.pro @ 59

;+
; NAME: compute_fromreg_bilinear_weigaddr
;
; PURPOSE: compute the weight and address neede to interpolate data from a
;          "regular grid" to any grid using the bilinear method
;   
; CATEGORY:interpolation
;
; CALLING SEQUENCE: 
;     compute_fromreg_bilinear_weigaddr, alon, alat, olon, olat, weig, addr
;
; INPUTS:
;     lonin and latin: longitude/latitude of the input data 
;     lonout and latout: longitude/latitude of the output data 
;
; KEYWORD PARAMETERS: 
;
;     /NONORTHERNLINE and /NOSOUTHERNLINE: activate if you don't whant to take into
;          account the northen/southern line of the input data when perfoming the
;          interpolation.
;
; OUTPUTS: 
;     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, jpio, jpjo, /over)
;
; COMMON BLOCKS: none
;
; SIDE EFFECTS: ?
;
; RESTRICTIONS:
;  -  the input grid must be a "regular 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 out of the southern and northern boundaries are interpolated
;     using a linear interpolation only along the longitudinal direction.
; 
; EXAMPLE: 
;
; MODIFICATION HISTORY:
;  November 2005: Sebastien Masson (smasson@lodyc.jussieu.fr) 
; 
;-
;
;----------------------------------------------------------
;----------------------------------------------------------
;
PRO compute_fromreg_bilinear_weigaddr, alonin, alatin, olonin, olat, weig, addr $
  , NONORTHERNLINE = nonorthernline, NOSOUTHERNLINE = nosouthernline
;
  compile_opt strictarr, 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 bondary 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 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
;
; we want to do biliear 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 armospheric 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 points 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 bondary 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 bondary 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 points 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 boudary of the atm grid 
;      or southward of the southern boudary 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)]] = 0L
    ybad = -1 & bad = -1 ; free memory
  ENDIF
; check totalweight = 1
  totalweig = abs(1-total(weig, 1))
  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