source: trunk/SRC/Interpolation/compute_fromirr_bilinear_weigaddr.pro @ 238

;+
;
; @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}{type=2d array}
; longitude of the input data
;
; @param olat {in}{required}{type=2d array}
; latitude of the input data
;
; @param omsk {in}{required}{type=2d array or -1}
; land/sea mask of the input data
; put -1 if input data are not masked
;
; @param alonin {in}{required}{type=2d array}
; longitude of the output data
;
; @param alat {in}{required}{type=2d array}
; latitude of the output data
;
; @param amsk {in}{required}{type=2d array or -1}
; land/sea mask of the output data
; put -1 if output data are not masked
;
; @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, 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 same value as their closest neighbor
;
; @history
;  June 2006: Sebastien Masson (smasson\@lodyc.jussieu.fr)
;
; @version
; $Id$
;
;-
;
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]
; mask check
  IF n_elements(omsk) EQ 1 AND omsk[0] EQ -1 THEN omsk = replicate(1b, jpio, jpjo)
  IF n_elements(amsk) EQ 1 AND amsk[0] EQ -1 THEN amsk = replicate(1b, jpia, jpja)
  IF n_elements(omsk) NE jpio*jpjo THEN BEGIN
    ras = report('input grid mask do not have the good size')
    stop
  ENDIF
  IF n_elements(amsk) NE jpia*jpja THEN BEGIN
    ras = report('output grid mask do not have the good size')
    stop
  ENDIF
;
; 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 ouput grid water points
  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
; longitude 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]
; now, we morph the quadrilateral ocean cell into the reference square (0 -> 1)
; in addition we get the coordinates 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 checks...
        tmpmsk = omsk[celladdr[*, ind]]
        CASE 1 OF
          xy[0] LT 0 OR xy[0] GT 1 : stop
          xy[1] LT 0 OR xy[1] GT 1 : stop
          xy[0] EQ 0 AND xy[1] EQ 0 AND tmpmsk[0] EQ 0 : foraddr[n] = -1
          xy[0] EQ 1 AND xy[1] EQ 0 AND tmpmsk[1] EQ 0 : foraddr[n] = -1
          xy[0] EQ 1 AND xy[1] EQ 1 AND tmpmsk[2] EQ 0 : foraddr[n] = -1
          xy[0] EQ 0 AND xy[1] EQ 1 AND tmpmsk[3] EQ 0 : foraddr[n] = -1
          xy[0] EQ 0 AND (tmpmsk[0]+tmpmsk[3]) EQ 0    : foraddr[n] = -1
          xy[0] EQ 1 AND (tmpmsk[1]+tmpmsk[2]) EQ 0    : foraddr[n] = -1
          xy[1] EQ 0 AND (tmpmsk[0]+tmpmsk[1]) EQ 0    : foraddr[n] = -1
          xy[1] EQ 1 AND (tmpmsk[2]+tmpmsk[3]) EQ 0    : foraddr[n] = -1
          ELSE: BEGIN
; we keep its address
        foraddr[n] = ind
; keep the new coordinates
            forweight[n, 0] = xy[0]
            forweight[n, 1] = xy[1]
          END
        ENDCASE



      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)
  IF min(totalweig, max = ma) EQ 0 then stop
  IF ma GT 1 then stop
  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