source: trunk/SRC/Interpolation/compute_fromirr_bilinear_weigaddr.pro

;+
;
; @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
;  -  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 weights are 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 neighbors
;
; @history
;
; - pinsardf 20120813T123457Z curie51.c-curie.tgcc.ccc.cea.fr (Linux)
;
;   * less crypted control print
;
; - June 2006: Sebastien Masson (smasson\@lodyc.jussieu.fr)
;
; @version
; $Id$
;
;-
PRO compute_fromirr_bilinear_weigaddr, olonin, olat, omsk, alonin, alat, amsk $
                                     , weig, addr, lonref1 = lonref1, onplane = onplane, notuse_neighbor = notuse_neighbor
;
  compile_opt idl2, strictarrsubs
;
  s = size(SCOPE_TRACEBACK(/STRUCTURE),/DIMENSION)
  routine = (SCOPE_TRACEBACK(/STRUCTURE))[s-1].Routine
;
  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 BEGIN
    omsk = replicate(1b, jpio, jpjo)
; if this is ORCA2 grid...
    IF (jpio EQ 180 OR jpio EQ 182) AND (jpjo EQ 149  OR jpjo EQ 148  OR jpjo EQ 147 ) THEN BEGIN
; we look for ill defined cells.
      IF jpio EQ 182 THEN lontmp = olonin[1:180, *] ELSE lontmp = olonin
; longitudinal size of the cells...
      a = (lontmp-shift(lontmp, 1, 0))
      d1 = 0.5 * max( [[[a]], [[360+a]]] MOD 360, dimension = 3)
      a = (lontmp-shift(lontmp, 1, 1))
      d2 = 0.5 * max( [[[a]], [[360+a]]] MOD 360, dimension = 3)
      a = (shift(lontmp, 0, 1)-shift(lontmp, 1, 0))
      d3 = 0.5 * max( [[[a]], [[360+a]]] MOD 360, dimension = 3)
      a = (shift(lontmp, 0, 1)-shift(lontmp, 1, 1))
      d4 = 0.5 * max( [[[a]], [[360+a]]] MOD 360, dimension = 3)
      md = [[[d1]], [[d2]], [[d3]], [[d4]]]
      bad = max(md, dimension = 3) GE 1.5
      bad = (bad + shift(bad, -1, -1) + shift(bad, 0, -1) + shift(bad, -1, 0)) < 1
      IF jpio EQ 182 THEN BEGIN
        omsk[1:180, 80:120] = 1b - bad[*, 80:120]
        omsk[0, *] = omsk[180, *]
        omsk[181, *] = omsk[1, *]
      ENDIF ELSE omsk[*, 80:120] = 1b - bad[*, 80:120]
    ENDIF
  ENDIF ELSE BEGIN
    IF n_elements(omsk) NE jpio*jpjo THEN BEGIN
      ras = report('input grid mask do not have the good size')
      stop
    ENDIF
  ENDELSE
  IF n_elements(amsk) EQ 1 AND amsk[0] EQ -1 THEN amsk = replicate(1b, jpia, jpja) ELSE BEGIN
    IF n_elements(amsk) NE jpia*jpja THEN BEGIN
      ras = report('output grid mask do not have the good size')
      stop
    ENDIF
  ENDELSE
;
  IF n_elements(lonref1) EQ 0 THEN lonref1 = 0.
  lonref2 = lonref1 + 360.
; longitude, between 0 and 360
  alon = (alonin + 720. ) MOD 360.
; longitude, between lonref1 and lonref2
  out = where(alon LT lonref1)
  IF out[0] NE -1 THEN alon[out] = alon[out] + 360.
  out = where(alon GT lonref2)
  IF out[0] NE -1 THEN alon[out] = alon[out] - 360.
; longitude, between 0 and 360
  olon = ( olonin + 720. ) MOD 360.
; longitude, between lonref1 and lonref2
  out = where(olon LT lonref1)
  IF out[0] NE -1 THEN olon[out] = olon[out] + 360.
  out = where(olon GT lonref2)
  IF out[0] NE -1 THEN olon[out] = olon[out] - 360.
;
; we work only on the output 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)
  alladdrm1 = shift(alladdr, -1)
  IF ((jpio EQ 92  ) AND (jpjo EQ 76   OR jpjo EQ 75   OR jpjo EQ 74  )) OR $
     ((jpio EQ 182 ) AND (jpjo EQ 149  OR jpjo EQ 148  OR jpjo EQ 147 )) OR $
     ((jpio EQ 722 ) AND (jpjo EQ 522  OR jpjo EQ 521  OR jpjo EQ 520 )) OR $
     ((jpio EQ 1442) AND (jpjo EQ 1021 OR jpjo EQ 1020 OR jpjo EQ 1019)) THEN alladdrm1[jpio-1, *] = alladdr[2, *]

  celladdr = lonarr(4, jpio*(jpjo-1))
  celladdr[0, *] = alladdr
  celladdr[1, *] = alladdrm1
  celladdr[2, *] = temporary(alladdrm1) + jpio
  celladdr[3, *] = temporary(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)])* 3.
  print,  routine, ' total loop length = ' , nawater
  FOR n = 0L, nawater-1 DO BEGIN
; control print
    IF (n MOD 5000) EQ 0 THEN print, routine, ' n = ' , 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 north 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 (lonref1+delta) OR xx GE (lonref2-delta):BEGIN
        lat1 = yy-delta
        lat2 = yy+delta
        good = where((minxcell LE (lonref1+2*delta) OR maxxcell GE (lonref2-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 90   OR jpio EQ 92  ) AND (jpjo EQ 76   OR jpjo EQ 75   OR jpjo EQ 74  ):onsphere = yy GT 40
          (jpio EQ 180  OR jpio EQ 182 ) AND (jpjo EQ 149  OR jpjo EQ 148  OR jpjo EQ 147 ):onsphere = yy GT 40
          (jpio EQ 720  OR jpio EQ 722 ) AND (jpjo EQ 522  OR jpjo EQ 521  OR jpjo EQ 520 ):onsphere = yy GT 40
          (jpio EQ 1440 OR jpio EQ 1442) AND (jpjo EQ 1021 OR jpjo EQ 1020 OR jpjo EQ 1019):onsphere = yy GT 40
          ELSE:onsphere = 1 - keyword_set(onplane) 
        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, delta = delta, /double)
; 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], /double)
        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, /double)
        ENDELSE
; take care of rounding errors...
        zero = where(abs(xy) LT 1e-6)
        IF zero[0] NE -1 THEN xy[zero] = 0.d
        one = where(abs(1.d - xy) LT 1e-6)
        IF one[0] NE -1 THEN xy[one] = 1.d
; 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
  
  IF keyword_set(notuse_neighbor) THEN BEGIN  ; don't do anything for bad points -> they will have no addr/weig like land points
    ok = where(foraddr NE -1, nawater)   ; keep only the points that are above a water oceanic cell.
    foraddr = foraddr[ok]
    forweight = forweight[ok, *]
    awater = awater[temporary(ok)]
  ENDIF
  
; 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))
    onsphere = 1
    print,  routine, ' total loop length = ', n_elements(bad)
    FOR i = 0L, n_elements(bad)-1L DO BEGIN
      IF (i MOD 500) EQ 0 THEN print, routine, ' i = ', i
      xtmp = alon[badaddr[i]]
      ytmp = alat[badaddr[i]]
      CASE 1 OF
; if we are near the north pole
        ytmp GE (90-delta):keep = where(goodlat GE 90-3*delta, cnt)
; if we are near the longitude periodicity area
        xtmp LE (lonref1+delta) OR xtmp GE (lonref2-delta):keep = where((goodlon LE (lonref1+3*delta) OR goodlon GE (lonref2-3*delta)) $
                                                          AND goodlat GE (ytmp-3*delta) AND goodlat LE (ytmp+3*delta), cnt)
; other cases
        ELSE:BEGIN
          keep = where(goodlon GE (xtmp-3*delta) AND goodlon LE (xtmp+3*delta) $
                       AND goodlat GE (ytmp-3*delta) AND goodlat le (ytmp+3*delta), cnt)
; ORCA cases : orca grid is irregular only northward of 40N
          CASE 1 OF
            (jpio EQ 90   OR jpio EQ 92  ) AND (jpjo EQ 76   OR jpjo EQ 75   OR jpjo EQ 74  ):onsphere = yy GT 40
            (jpio EQ 180  OR jpio EQ 182 ) AND (jpjo EQ 149  OR jpjo EQ 148  OR jpjo EQ 147 ):onsphere = yy GT 40
            (jpio EQ 720  OR jpio EQ 722 ) AND (jpjo EQ 522  OR jpjo EQ 521  OR jpjo EQ 520 ):onsphere = yy GT 40
            (jpio EQ 1440 OR jpio EQ 1442) AND (jpjo EQ 1021 OR jpjo EQ 1020 OR jpjo EQ 1019):onsphere = yy GT 40
            ELSE:
          ENDCASE
        END
      ENDCASE
      IF cnt NE 0 THEN BEGIN
        neig[i] = keep[(neighbor(xtmp, ytmp, goodlon[keep], goodlat[keep], sphere = onsphere))[0]]
      ENDIF ELSE neig[i] = (neighbor(alon[badaddr[i]], alat[badaddr[i]], goodlon, goodlat, sphere = onsphere))[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.d - a)*(1.d - b), (1.d - b)*a, a*b, (1.d - a)*b]
  a = -1 &  b = -1              ; free memory
; mask the weight to suppress the corner located on land
  newaweig = newaweig*(omsk[newaaddr])
; for cell with some land corner,
; we have to redistribute the weight on the remaining water corners
; weights normalization
  totalweig = total(newaweig, 1, /double)
;;   IF min(totalweig, max = ma) LE 0.d then stop
;;   IF ma- 1.d GT 1.e-6 then stop
  newaweig = newaweig/(replicate(1.d, 4)#totalweig)
; weights
  weig = dblarr(4, jpia*jpja)
  weig[*, awater] = temporary(newaweig)
; address
  addr = lonarr(4, jpia*jpja)
  addr[*, awater] = temporary(newaaddr)
;
  RETURN
END