Changeset 336 for trunk

Timestamp:

03/20/08 22:36:46 (16 years ago)

Author:

smasson

Message:

improve boxzoom usage in read_ncdf

File:

• trunk/SRC/ToBeReviewed/LECTURE/read_ncdf.pro (modified) (12 diffs)

trunk/SRC/ToBeReviewed/LECTURE/read_ncdf.pro

@@ -77 +77 @@
 ; For ROMS outputs. To define zeta to 0. instead of reading it
 ;
+; @keyword ZINVAR {type=named variable}
+; Set this keyword to a named variable in which 1 is returned if a
+; vertical dimension is found in the variable. Returns 0 otherwise
+;
 ; @keyword _EXTRA
 ; Used to pass keywords to <pro>isafile</pro>, <pro>initncdf</pro>,
@@ -98 +102 @@
 ;-
 FUNCTION read_ncdf, name, beginning, ending, compatibility $
-                  , BOXZOOM=boxzoom, FILENAME=filename $
-                  , PARENTIN=parentin, TIMESTEP=timestep $
-                  , ADDSCL_BEFORE=addscl_before $
-                  , TOUT=tout, NOSTRUCT=nostruct, CONT_NOFILL=CONT_NOFILL, INIT=init $
-                  , GRID=grid, CALLITSELF=callitself, DIREC=direc $
-                  , ZETAFILENAME=zetafilename, ZETAZERO=zetazero $
-                  , _EXTRA=ex
+                    , BOXZOOM = boxzoom, FILENAME = filename $
+                    , PARENTIN = parentin, TIMESTEP = timestep $
+                    , ADDSCL_BEFORE = addscl_before $
+                    , TOUT = tout, NOSTRUCT = nostruct, CONT_NOFILL = CONT_NOFILL, INIT = init $
+                    , GRID = grid, CALLITSELF = callitself, DIREC = direc $
+                    , ZETAFILENAME = zetafilename, ZETAZERO = zetazero $
+                    , ZINVAR = zinvar, _EXTRA = ex
 ;
   compile_opt idl2, strictarrsubs
@@ -129 +133 @@
 ;------------------------------------------------------------
   IF size(filename, /type) NE 7 THEN $
-    return, report('read_ncdf cancelled')
+     return, report('read_ncdf cancelled')
   IF !version.OS_FAMILY EQ 'unix' THEN spawn, '\file '+filename+' > /dev/null'
   cdfid = ncdf_open(filename)
@@ -148 +152 @@
     dimnames[i] = tmp
   ENDFOR
+;------------------------------------------------------------
+; we check if the variable has a vertical dimension and/or 
+; a record dimension. This is useful to define boxzoom 
+; keyword in domdef
+;------------------------------------------------------------
+  dummy = where(varinq.dim EQ inq.recdim, tinvar) 
+; check the presence of a vertical dimension according to the
+; number of dimensions and the presence of a record dimension
+  zinvar = (varinq.ndims EQ 3 AND tinvar NE 1) OR varinq.ndims EQ 4
 ;------------------------------------------------------------
 ; shall we redefine the grid parameters
@@ -175 +188 @@
       date1 = date2jul(beginning[0])
       if n_elements(ending) NE 0 then date2 = date2jul(ending[0]) ELSE date2 = date1
-      firsttps = (where(abs(time - date1) LT 0.9d/86400.d))[0]   ; beware of rounding errors...
+      firsttps = (where(abs(time - date1) LT 0.9d/86400.d))[0] ; beware of rounding errors...
       lasttps = (where(abs(time - date2) LT 0.9d/86400.d))[0]
       jpt = lasttps-firsttps+1
@@ -218 +231 @@
 ; are we in one of the case corresponding to ROMS conventions?
       CASE 1 OF
-        dimnames[2 <(varinq.ndims-1)] EQ 's_w':vargrid = 'W'
+        dimnames[2 < (varinq.ndims-1)] EQ 's_w':vargrid = 'W'
         dimnames[0] EQ 'xi_rho' AND dimnames[1] EQ 'eta_rho':vargrid = 'T'
         dimnames[0] EQ 'xi_u'   AND dimnames[1] EQ 'eta_u'  :vargrid = 'U'
@@ -241 +254 @@
         END
       ENDCASE
-     ENDIF
+    ENDIF
   ENDELSE
 ;---------------------------------------------------------------
@@ -264 +277 @@
     endcase
   ENDIF ELSE BEGIN
-    if keyword_set(boxzoom) then BEGIN
-      Case 1 Of
-        N_Elements(Boxzoom) Eq 1:bte = [lon1, lon2, lat1, lat2, 0., boxzoom[0]]
-        N_Elements(Boxzoom) Eq 2:bte = [lon1, lon2, lat1, lat2, boxzoom[0], boxzoom[1]]
-        N_Elements(Boxzoom) Eq 4:bte = [Boxzoom, vert1, vert2]
-        N_Elements(Boxzoom) Eq 5:bte = [Boxzoom[0:3], 0, Boxzoom[4]]
-        N_Elements(Boxzoom) Eq 6:bte = Boxzoom
-        Else: BEGIN
+    IF keyword_set(boxzoom) then BEGIN
+      CASE N_Elements(Boxzoom) Of
+        1:bte = [lon1, lon2, lat1, lat2, 0., boxzoom[0]]
+        2:bte = [lon1, lon2, lat1, lat2, boxzoom[0], boxzoom[1]]
+        4:bte = [Boxzoom, vert1, vert2]
+        5:bte = [Boxzoom[0:3], 0, Boxzoom[4]]
+        6:bte = Boxzoom
+        ELSE: BEGIN
           ncdf_close, cdfid
           return, report('Wrong Definition of Boxzoom')
-        end
+        END
       ENDCASE
-      savedbox = 1b
-      saveboxparam, 'boxparam4rdncdf.dat'
-      domdef, bte, GRIDTYPE = ['T', vargrid], _extra = ex
+      IF zinvar EQ 1 THEN domdef, bte, GRIDTYPE = ['T', vargrid], _extra = ex $
+      ELSE domdef, bte[0:3], GRIDTYPE = ['T', vargrid], _extra = ex
     ENDIF
     grille, mask, glam, gphi, gdep, nx, ny, nz, firstx, firsty, firstz, lastx, lasty, lastz
@@ -325 +337 @@
   vardate = string(format = '(C(CMoA))', 31*(month-1))+' '+strtrim(day, 1)+', '+strtrim(year, 1)
   varexp = file_basename(filename)
-
 ; We apply the value valmask on land points.
   if NOT keyword_set(cont_nofill) then begin
@@ -344 +355 @@
   IF size(missing_value, /type) NE 7 THEN BEGIN
     IF finite(missing_value) EQ 1 THEN BEGIN
-        CASE 1 OF
-          missing_value GT 1.e6:missing = where(res GT missing_value/10.)
-          missing_value LT -1.e6:missing = where(res LT missing_value/10.)
-          abs(missing_value) LT 1.e-6:missing = where(abs(res) LT 1.e-6)
-          ELSE:missing = where(res EQ missing_value)
-        ENDCASE
+      CASE 1 OF
+        missing_value GT 1.e6:missing = where(res GT missing_value/10.)
+        missing_value LT -1.e6:missing = where(res LT missing_value/10.)
+        abs(missing_value) LT 1.e-6:missing = where(abs(res) LT 1.e-6)
+        ELSE:missing = where(res EQ missing_value)
+      ENDCASE
     ENDIF ELSE missing = where(finite(res) EQ 0)
   ENDIF ELSE missing = -1
@@ -358 +369 @@
   IF missing[0] NE -1 THEN res[temporary(missing)] = !values.f_nan
   IF earth[0] NE -1 THEN BEGIN 
-    IF varinq.ndims eq 3 AND (where(varinq.dim EQ inq.recdim))[0] NE -1 THEN $   ;xyt array    
+    IF varinq.ndims eq 3 AND (where(varinq.dim EQ inq.recdim))[0] NE -1 THEN $ ;xyt array    
        earth = earth#replicate(1, jpt)+replicate(nx*ny, n_elements(earth))#lindgen(jpt)
     IF varinq.ndims eq 4 THEN earth = earth#replicate(1, jpt)+replicate(nx*ny*nz, n_elements(earth))#lindgen(jpt)
@@ -405 +416 @@
 ;---------------------------------------------------------------------
   ncdf_close, cdfid
-
-  IF keyword_set(savedbox) THEN restoreboxparam, 'boxparam4rdncdf.dat'
-
   IF keyword_set(nostruct) THEN return, res
   IF keyword_set(key_forgetold) THEN BEGIN