Changeset 303 for trunk/SRC/Interpolation

Timestamp:

10/16/07 14:59:26 (17 years ago)

Author:

smasson

Message:

bugfix + speedup

Location:

trunk/SRC/Interpolation

Files:

• compute_fromirr_bilinear_weigaddr.pro (modified) (6 diffs)
• file_interp.pro (modified) (2 diffs)
• inquad.pro (modified) (4 diffs)

trunk/SRC/Interpolation/compute_fromirr_bilinear_weigaddr.pro

@@ -130 +130 @@
 ;
   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, *] = shift(alladdr, -1)
-  celladdr[2, *] = shift(alladdr + jpio, -1)
-  celladdr[3, *] = alladdr + jpio
+  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
@@ -153 +159 @@
 ; We look for which ocean water cell contains the atmosphere water point
 ;
-  delta = max([(360./jpio), (180./jpjo)])* 4.
+  delta = max([(360./jpio), (180./jpjo)])* 3.
   FOR n = 0L, nawater-1 DO BEGIN
 ; control print
@@ -185 +191 @@
 ; 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 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
@@ -202 +208 @@
                    , xcell[2, good], ycell[2, good] $
                    , xcell[1, good], ycell[1, good] $
-                   , onsphere = onsphere, newcoord = newcoord, /noprint)
+                   , 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]
@@ -266 +272 @@
 ; 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]
+    onsphere = 1
+    FOR i = 0L, n_elements(bad)-1L DO BEGIN
+      IF (i MOD 500) EQ 0 THEN print, 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 delta OR xtmp GE (360-delta):keep = where((goodlon LE 3*delta OR goodlon GE (360-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
@@ -274 +306 @@
     foraddr[bad] = foraddr[neig]
     forweight[bad, *] = forweight[neig, *]
-  endif
+  ENDIF
 ; transform the address of the ocean cell into the address of its 4 corners
   newaaddr = celladdr[*, temporary(foraddr)]

trunk/SRC/Interpolation/file_interp.pro

@@ -402 +402 @@
                 off = [0, 0, k]
                 ncdf_varget, inid, i, data, offset = off, count = incnt
+                IF n_elements(inmasksz) GE 3 THEN BEGIN
+                  IF inmasksz[2] EQ indimsz[varinq.dim[2]] AND varinq.dim[2] NE ininq.recdim THEN tmpmsk = inmask[*, *, k] $
+                  ELSE tmpmsk = inmask[*, *, 0]
+                ENDIF ELSE tmpmsk = inmask[*, *, 0]
                 IF inmasksz[2] EQ indimsz[varinq.dim[2]] AND varinq.dim[2] NE ininq.recdim THEN tmpmsk = inmask[*, *, k] $
                 ELSE tmpmsk = inmask[*, *, 0]
@@ -414 +418 @@
                 FOR k = 0, indimsz[varinq.dim[2]]-1 DO BEGIN
                   incnt = [indimsz[varinq.dim[0: 1]], 1, 1]
-                  outcnt = [outdimsz[varinq.dim[0: 1]], 1]
+                  outcnt = [outdimsz[varinq.dim[0: 1]], 1, 1]
                   off = [0, 0, k, t]
                   ncdf_varget, inid, i, data, offset = off, count = incnt
-                  IF inmasksz[2] EQ indimsz[varinq.dim[2]] THEN tmpmsk = inmask[*, *, k] ELSE tmpmsk = inmask
+                  IF n_elements(inmasksz) GE 3 THEN BEGIN
+                    IF inmasksz[2] EQ indimsz[varinq.dim[2]] THEN tmpmsk = inmask[*, *, k] ELSE tmpmsk = inmask
+                  ENDIF ELSE tmpmsk = inmask[*, *, 0]
                   IF interp THEN data = call_interp2d(temporary(data), inlon, inlat, temporary(tmpmsk), outlon, outlat $
                                                       , INIRR = inirr, METHOD = method, SMOOTH = smooth $

trunk/SRC/Interpolation/inquad.pro

@@ -33 +33 @@
 ; automatically.
 ;
-; @keyword ZOOMRADIUS {default=4}
-; the zoom (circle centered on the (x,y) with a radius of
-; zoomradius degree where we look for the quadrilateral which
-; contains the (x,y) point) used for the satellite projection
-; when /ONSPHERE is activated.
-; 4 seems to be the minimum which can be used.
-; Can be increase if the cell size is larger than 5 degrees.
+; @keyword DELTA {default=4}
+; to speed up the program, we reduce the aera where we look for potential
+; quadrilaterals containing (x,y). Delta defines the limit of the box
+; centred on (x,y) with a zonal and meridional extent of delta degrees.
 ;
 ; @keyword NOPRINT
@@ -85 +82 @@
 ;-
 ;
-FUNCTION inquad, x, y, x1, y1, x2, y2, x3, y3, x4, y4, ONSPHERE = onsphere,  DOUBLE = double, ZOOMRADIUS = zoomradius, NOPRINT = noprint, NEWCOORD = newcoord
+FUNCTION inquad, x, y, x1, y1, x2, y2, x3, y3, x4, y4, ONSPHERE = onsphere,  DOUBLE = double, DELTA = delta, NOPRINT = noprint, NEWCOORD = newcoord
 ;
   compile_opt idl2, strictarrsubs
@@ -113 +110 @@
     save = {map:!map, x:!x, y:!y, z:!z, p:!p}
 ; do a satellite projection
-    IF NOT keyword_set(zoomradius) THEN zoomradius = 4
-    map_set, y[0], x[0], 0, /satellite, sat_p = [1+zoomradius*20/6371.229, 0, 0], /noerase, /iso, /noborder
+    IF NOT keyword_set(delta) THEN delta = 4
+    map_set, y[0], x[0], 0, /stereo, limit = [y[0]-delta/2, x[0]-delta/2, y[0]+delta/2, x[0]+delta/2], /noerase, /iso, /noborder
 ; use normal coordinates to reject cells which are out of the projection.
     tmp  = convert_coord(x, y, /DATA, /TO_NORMAL, DOUBLE = double)
@@ -218 +215 @@
 ; quadrilateral.
 ; row j of test corresponds to all the points localized in cell j
+
+  IF keyword_set(double) THEN one = 1.d ELSE one = 1. 
+  nquad_1 = replicate(one, nquad)
+  ntofind_1 = replicate(one, ntofind)
+  x_nquad = x[*]#replicate(one, nquad)
+  y_nquad = y[*]#replicate(one, nquad)
+
   test = $
-; (x-x1)
-  ((x[*]#replicate(1, nquad)-replicate(1, ntofind)#x1[*]) $
-; *(y2-y1)
-  *(replicate(1, ntofind)#(y2-y1)[*]) $
-; GE (x2-x1)
-  GE ((replicate(1, ntofind)#(x2-x1)[*]) $
-; *(y-y1)
-  *(y[*]#replicate(1, nquad)-replicate(1, ntofind)#y1[*])))
-;-------
-  test = temporary(test) $
-; *(x-x2)
-  *((x[*]#replicate(1, nquad)-replicate(1, ntofind)#x2[*]) $
-; *(y3-y2)
-  *(replicate(1, ntofind)#(y3-y2)[*]) $
-; GE (x3-x2)
-  GE ((replicate(1, ntofind)#(x3-x2)[*]) $
-; *(y-y2)
-  *(y[*]#replicate(1, nquad)-replicate(1, ntofind)#y2[*])))
-;-------
-  test = temporary(test) $
-; *(x-x3)
-  *((x[*]#replicate(1, nquad)-replicate(1, ntofind)#x3[*]) $
-; *(y4-y3)
-  *(replicate(1, ntofind)#(y4-y3)[*]) $
-; GE (x4-x3)
-  GE ((replicate(1, ntofind)#(x4-x3)[*]) $
-; *(y-y3)
-  *(y[*]#replicate(1, nquad)-replicate(1, ntofind)#y3[*])))
-;-------
-  test = temporary(test) $
-; *(x-x4)
-  *((x[*]#replicate(1, nquad)-replicate(1, ntofind)#x4[*]) $
-; *(y1-y4)
-  *(replicate(1, ntofind)#(y1-y4)[*]) $
-; GE (x1-x4)
-  GE ((replicate(1, ntofind)#(x1-x4)[*]) $
-; *(y-y4)
-  *(y[*]#replicate(1, nquad)-replicate(1, ntofind)#y4[*])))
+;  (x-x1)*(y2-y1) GE (x2-x1)*(y-y1)
+  ( (x_nquad - ntofind_1#x1[*]) * (ntofind_1#(y2-y1)[*]) ) GE ( (ntofind_1#(x2-x1)[*]) * (y_nquad - ntofind_1#y1[*]) ) AND $
+; *(x-x2)*(y3-y2) GE (x3-x2)*(y-y2)
+  ( (x_nquad - ntofind_1#x2[*]) * (ntofind_1#(y3-y2)[*]) ) GE ( (ntofind_1#(x3-x2)[*]) * (y_nquad - ntofind_1#y2[*]) ) AND $
+; *(x-x3)*(y4-y3) GE (x4-x3)*(y-y3)
+  ( (x_nquad - ntofind_1#x3[*]) * (ntofind_1#(y4-y3)[*]) ) GE ( (ntofind_1#(x4-x3)[*]) * (y_nquad - ntofind_1#y3[*]) ) AND $
+; *(x-x4)*(y1-y4) GE (x1-x4)*(y-y4)
+  ( (x_nquad - ntofind_1#x4[*]) * (ntofind_1#(y1-y4)[*]) ) GE ( (ntofind_1#(x1-x4)[*]) * (y_nquad - ntofind_1#y4[*]) )
+;
+  nquad_1 = 1 & ntofind_1 = 1 & x_nquad = 1 & y_nquad = 1   ; free memory
 ;
 ; check test if ntofind gt 1