Changeset 172 for trunk/SRC/ToBeReviewed/LECTURE/read_ncdf.pro

Timestamp:

09/11/06 09:11:26 (18 years ago)

Author:

smasson

Message:

bugfix + manage roms outputs

File:

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

trunk/SRC/ToBeReviewed/LECTURE/read_ncdf.pro

@@ -91 +91 @@
                     , PARENTIN = parentin, TIMESTEP = timestep, TIMEVAR = timevar $
                     , TOUT = tout, NOSTRUCT = nostruct, CONT_NOFILL = CONT_NOFILL, INIT = init $
-                    , GRID = grid, FBASE2TBASE = fbase2tbase, _EXTRA = ex
+                    , GRID = grid, FBASE2TBASE = fbase2tbase, _EXTRA = ex, CALLITSELF = callitself
 ;---------------------------------------------------------
 ;
@@ -129 +129 @@
   ENDIF
   varcontient = ncdf_varinq(cdfid, name)
+  IF varcontient.ndims LT 2 THEN return, report('read_ncdf cannot read scalar or 1D data')
+; look for the dimension names
+  dimnames = strarr(varcontient.ndims)
+  FOR i = 0, varcontient.ndims-1 DO BEGIN
+    ncdf_diminq, cdfid, varcontient.dim[i], tmp, dimsize
+    dimnames[i] = tmp
+  ENDFOR 
 ;------------------------------------------------------------
 ; shall we redefine the grid parameters
@@ -144 +151 @@
     if n_elements(ending) NE 0 then lasttps = ending[0] ELSE lasttps = firsttps
     jpt = lasttps-firsttps+1
-    time = julday(1, 1, 1) + lindgen(jpt)
+    IF NOT keyword_set(callitself) then time = julday(1, 1, 1) + lindgen(jpt)
   ENDIF ELSE BEGIN
     if keyword_set(parent) then BEGIN
@@ -236 +243 @@
       mots = str_sep(value, ' ')
       unite = mots[0]
+      IF unite NE 'seconds' AND unite NE 'hours' AND unite NE 'days' $
+         AND unite NE 'months' AND unite NE 'years' THEN BEGIN
+        ncdf_close, cdfid
+        return, report('time units does not start with seconds/hours/days/months/years')
+      ENDIF
+      IF stregex(value, '[^ ]* since ([0-9]){4}-([0-9]){2}-([0-9]){2}.*', /boolean) EQ 0 THEN BEGIN
+        ncdf_close, cdfid
+        return, report('attribut units of time has not the good format: [^ ]* since ([0-9]){4}-([0-9]){2}-([0-9]){2}.*')
+      ENDIF
       depart = str_sep(mots[2], '-')
       ncdf_varget, cdfid, timeid, time
@@ -260 +276 @@
         ELSE:BEGIN
           ncdf_close, cdfid
-          return, report('The "units" attribu of the time axis must be something like: !C "seconds since 0001-01-01 ..." !C "days since 1979-01-01 ..." !C "months since 1979-01-01 ..." !C "years since 1979-01-01 ..."')
+          return, report('The "units" attribute of the time axis must be something like: !C "seconds since 0001-01-01 ..." !C "days since 1979-01-01 ..." !C "months since 1979-01-01 ..." !C "years since 1979-01-01 ..."')
         end
       ENDCASE
@@ -291 +307 @@
     vargrid = 'T'               ; default definition
     IF finite(glamu[0]) EQ 1 THEN BEGIN
-      pattern = ['GRID.', 'GRID_', 'GRID', 'UPID_', '30ID_']
-      gdtype = ['T', 'U', 'V', 'W', 'F']
-      fnametest = strupcase(filename)
-      FOR i = 0, n_elements(pattern)-1 DO BEGIN
-        FOR j = 0, n_elements(gdtype)-1 DO BEGIN
-          substr = pattern[i]+gdtype[j]
-          pos = strpos(fnametest, substr)
-          IF pos NE -1 THEN $
-             vargrid = strmid(fnametest, pos+strlen(substr)-1, 1)
-        ENDFOR
-      ENDFOR
-    ENDIF
+; are we in one of the case corresponding to ROMS conventions?
+      CASE 1 OF
+        dimnames[2 <(varcontient.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'
+        dimnames[0] EQ 'xi_v'   AND dimnames[1] EQ 'eta_v'  :vargrid = 'V'
+        dimnames[0] EQ 'xi_psi' AND dimnames[1] EQ 'eta_psi':vargrid = 'F'
+        dimnames[0] EQ 'xi_rho' AND dimnames[1] EQ 'eta_v'  :vargrid = 'V'
+        dimnames[0] EQ 'xi_u'   AND dimnames[1] EQ 'eta_rho':vargrid = 'U'
+        dimnames[0] EQ 'xi_u'   AND dimnames[1] EQ 'eta_v'  :vargrid = 'F'
+        ELSE:BEGIN 
+; could we define the grid type from the file name??
+          pattern = ['GRID.', 'GRID_', 'GRID', 'UPID_', '30ID_']
+          gdtype = ['T', 'U', 'V', 'W', 'F']
+          fnametest = strupcase(filename)
+          FOR i = 0, n_elements(pattern)-1 DO BEGIN
+            FOR j = 0, n_elements(gdtype)-1 DO BEGIN
+              substr = pattern[i]+gdtype[j]
+              pos = strpos(fnametest, substr)
+              IF pos NE -1 THEN $
+                 vargrid = strmid(fnametest, pos+strlen(substr)-1, 1)
+            ENDFOR
+          ENDFOR
+        END
+      ENDCASE
+     ENDIF
   ENDELSE
 ;---------------------------------------------------------------
@@ -380 +410 @@
   ENDIF
 ;
+  IF n_elements(key_zreverse) EQ 0 THEN key_zreverse = 0
+  IF keyword_set(key_zreverse) THEN BEGIN
+    tmp = jpk-1-firstz
+    firstz = jpk-1-lastz
+    lastz = tmp
+  ENDIF
+;
   IF keyword_set(fbase2tbase) THEN BEGIN
     case strupcase(vargrid) of
@@ -419 +456 @@
 ;---------------------------------------------------------------------
 ; varname
-  varname = name
+  IF NOT keyword_set(callitself) THEN varname = name
 ; varunit
   if varcontient.natts NE 0 then begin
@@ -428 +465 @@
     found = (where(lowattnames EQ 'units'))[0]
     IF found NE -1 then ncdf_attget, cdfid, name, attnames[found], value ELSE value = ''
-    varunit = strtrim(string(value), 2)
+    IF NOT keyword_set(callitself) THEN varunit = strtrim(string(value), 2)
 ;
     found = (where(lowattnames EQ 'add_offset'))[0]
@@ -443 +480 @@
 ;
   ENDIF ELSE BEGIN 
-    varunit = ''
+    IF NOT keyword_set(callitself) THEN varunit = ''
     add_offset = 0.
     scale_factor = 1.
@@ -457 +494 @@
 
 ; we apply reverse
-  if keyword_set(key_yreverse) then res = reverse(temporary(res),  2)
-  if keyword_set(key_zreverse) AND (size(res))[0] EQ 3 AND jpt EQ 1 then res = reverse(temporary(res), 3)
-  if keyword_set(key_zreverse) AND (size(res))[0] EQ 4 THEN res = reverse(temporary(res), 3)
-
+  if keyword_set(key_yreverse) AND ny NE 1 THEN $
+     res = reverse(reform(res, nx, ny, nz, jpt, /overwrite),  2)
+  if keyword_set(key_zreverse) AND nz NE 1 $
+     AND varcontient.ndims - ((where(varcontient.dim EQ contient.recdim))[0] NE -1) EQ 3 THEN $
+        res = reverse(reform(res, nx, ny, nz, jpt, /overwrite),  3)
 ; We apply the value valmask on land points.
   if NOT keyword_set(cont_nofill) then begin
@@ -501 +539 @@
   if add_offset NE 0 then res = temporary(res)+add_offset
   if missing[0] NE -1 then res[temporary(missing)] = !values.f_nan
-  if earth[0] NE -1 then res[temporary(earth)] = 1e20
+  if earth[0] NE -1 then res[temporary(earth)] = 1.e20
+;---------------------------------------------------------------------
+; if it is roms outputs, we need to get additionals infos...
+  IF NOT keyword_set(callitself) THEN BEGIN
+    IF strmid(dimnames[0], 0, 3) EQ 'xi_' AND strmid(dimnames[1], 0, 4) EQ 'eta_' THEN BEGIN
+      ncdf_attget, cdfid, 'theta_s', theta_s, /global
+      ncdf_attget, cdfid, 'theta_b', theta_b, /global
+      ncdf_attget, cdfid, 'hc', hc, /global
+; +++ binder l'exsitance de h et zeta...
+; +++ mettre zeta a 0 par defaut
+      hroms = read_ncdf('h', 0, 0, FILENAME = filename $
+                        , /TIMESTEP, /NOSTRUCT, CONT_NOFILL = CONT_NOFILL $
+                        , GRID = vargrid, /CALLITSELF, _EXTRA = ex)
+      zeta = read_ncdf('zeta', firsttps, lasttps, FILENAME = filename $
+                       , /TIMESTEP, /NOSTRUCT, CONT_NOFILL = CONT_NOFILL $
+                       , GRID = vargrid, /CALLITSELF, _EXTRA = ex)
+      romszinfos = {h:temporary(hroms), zeta:temporary(zeta), theta_s:theta_s, theta_b:theta_b, hc:hc}
+    ENDIF ELSE romszinfos = {h:-1, zeta:-1, theta_s:-1, theta_b:-1, hc:-1}
+  ENDIF
 ;---------------------------------------------------------------------
   ncdf_close, cdfid
-;---------------------------------------------------------------------
-  if keyword_set(savedbox) THEN restoreboxparam, 'boxparam4rdncdf.dat'
-  if keyword_set(nostruct) then return, res $
-  ELSE BEGIN 
-    IF keyword_set(key_forgetold) THEN BEGIN
-      return, {arr:res, grid:vargrid, unit:varunit, experiment:varexp, name:varname} 
-    ENDIF ELSE BEGIN 
-      return, {tab:res, grille:vargrid, unite:varunit, experience:varexp, nom:varname}
-    ENDELSE 
-  ENDELSE 
+
+  IF keyword_set(savedbox) THEN restoreboxparam, 'boxparam4rdncdf.dat'
+
+  IF keyword_set(nostruct) THEN return, res
+  IF keyword_set(key_forgetold) THEN BEGIN
+    return, {arr:temporary(res), grid:vargrid, unit:varunit, experiment:varexp, name:varname} 
+  ENDIF ELSE BEGIN 
+    return, {tab:temporary(res), grille:vargrid, unite:varunit, experience:varexp, nom:varname}
+  ENDELSE
+
 END