Changeset 396

Timestamp:

04/08/09 18:14:31 (15 years ago)

Author:

smasson

Message:

bugfix on leapyr calendar with months le 0 or months gt 12 + add true julian day as unit for calendar

Location:

trunk/SRC

Files:

• Calendar/julday.pro (modified) (1 diff)
• ReadWrite/ncdf_getatt.pro (modified) (1 diff)
• ReadWrite/ncdf_gettime.pro (modified) (5 diffs)

trunk/SRC/Calendar/julday.pro

@@ -257 +257 @@
       L_YEAR = LONG(YEAR)
 ;
-      tochange = where(L_MONTH LT 0)
-      IF tochange[0] NE -1 THEN BEGIN
-        L_YEAR[tochange] = L_YEAR[tochange]+L_MONTH[tochange]/12-1
-        L_MONTH[tochange] =  12 + L_MONTH[tochange] MOD 12
-      ENDIF
-;
-      tochange = where(L_MONTH GT 12)
-      IF tochange[0] NE -1 THEN BEGIN
-        L_YEAR[tochange] = L_YEAR[tochange]+L_MONTH[tochange]/12
-        L_MONTH[tochange] =  L_MONTH[tochange] MOD 12
+      tochange = where(L_MONTH LE 0)
+      IF tochange[0] NE -1 THEN BEGIN
+        tmp = L_MONTH[tochange] - 12L
+        L_YEAR[tochange] = L_YEAR[tochange] + tmp/12L
+        L_MONTH[tochange] =  (temporary(tmp) MOD 12L) + 12L
+      ENDIF
+;
+      tochange = where(L_MONTH GT 12L)
+      IF tochange[0] NE -1 THEN BEGIN
+        tmp = L_MONTH[tochange] - 1L
+        L_YEAR[tochange] = L_YEAR[tochange] + tmp/12L
+        L_MONTH[tochange] =  (temporary(tmp) MOD 12L) + 1L
       ENDIF
 ;

trunk/SRC/ReadWrite/ncdf_getatt.pro

@@ -102 +102 @@
         'units':BEGIN 
           ncdf_attget, cdfid, varid, attname, tmp
-          units = strtrim(tmp, 2)
+          units = strcompress(strtrim(tmp, 2))
         END
         'calendar':BEGIN 

trunk/SRC/ReadWrite/ncdf_gettime.pro

@@ -61 +61 @@
 ;
   IF n_elements(cdfid) EQ 0 THEN BEGIN
-    cdfid = ncdf_open(isafile(filename, title = 'which file must be open by ncdf_gettime?', IODIRECTORY=iodir, _extra = ex))
+    cdfid = ncdf_open(isafile(filename, title = 'which file must be open by ncdf_gettime?', IODIRECTORY = iodir, _extra = ex))
     tobeclosed = 1
   ENDIF
@@ -149 +149 @@
 ;----------------------------------------------------
 ;
+  value = strlowcase(value)
+  IF value NE 'true julian day' THEN BEGIN 
 ; time_counter:units = "seconds since 0001-01-01 00:00:00" ;
 ; time_counter:units = "hours since 0001-01-01 00:00:00" ;
@@ -156 +158 @@
 ;
 ; we decript the "units" attribute to find the time origin
-  value = strtrim(strcompress(string(value)), 2)
-  words = str_sep(value, ' ')
-  units = words[0]
-  units = strlowcase(units)
-  IF strpos(units, 's', strlen(units)-1) NE -1 THEN units = strmid(units, 0, strlen(units)-1)
-  IF strpos(units, 'julian_') NE -1 THEN units = strmid(units, 7)
-  IF units NE 'second' AND units NE 'hour' AND units NE 'day' $
-     AND units NE 'month' AND units NE 'year' THEN BEGIN
-    CASE caller OF
-      'read_ncdf':err = 'time units does not start with seconds/hours/days/months/years !C Use the TIMESTEP keyword'
-      'scanfile':err = 'time units does not start with seconds/hours/days/months/years !C we create a fake calendar ...'
-    ENDCASE
-    IF keyword_set(tobeclosed) THEN ncdf_close, cdfid
-    return, -jpt
-  ENDIF
-  IF stregex(value, '[^ ]* since ([0-9]){1,4}-([0-9]){1,2}-([0-9]){1,2}.*', /boolean) EQ 0 THEN BEGIN
-    CASE caller OF
-      'read_ncdf':err = 'attribute units of time has not the good format: [^ ]* since ([0-9]){1,4}-([0-9]){1,2}-([0-9]){1,2}.*!C Use the TIMESTEP keyword'
-      'scanfile':err = 'attribute units of time has not the good format: [^ ]* since ([0-9]){1,4}-([0-9]){1,2}-([0-9]){1,2}.*!C we create a fake calendar ...'
-    ENDCASE
-    IF keyword_set(tobeclosed) THEN ncdf_close, cdfid
-    return, -jpt
-  ENDIF
-  start = str_sep(words[2], '-')
+    words = str_sep(value, ' ')
+    units = words[0]
+    IF strpos(units, 's', strlen(units)-1) NE -1 THEN units = strmid(units, 0, strlen(units)-1)
+    IF strpos(units, 'julian_') NE -1 THEN units = strmid(units, 7)
+    IF units NE 'second' AND units NE 'hour' AND units NE 'day' $
+       AND units NE 'month' AND units NE 'year' THEN BEGIN
+      CASE caller OF
+        'read_ncdf':err = 'time units does not start with seconds/hours/days/months/years !C Use the TIMESTEP keyword'
+        'scanfile':err = 'time units does not start with seconds/hours/days/months/years !C we create a fake calendar ...'
+      ENDCASE
+      IF keyword_set(tobeclosed) THEN ncdf_close, cdfid
+      return, -jpt
+    ENDIF
+    IF stregex(value, '[^ ]* since ([0-9]){1,4}-([0-9]){1,2}-([0-9]){1,2}.*', /boolean) EQ 0 THEN BEGIN
+      CASE caller OF
+        'read_ncdf':err = 'attribute units of time has not the good format: [^ ]* since ([0-9]){1,4}-([0-9]){1,2}-([0-9]){1,2}.*!C Use the TIMESTEP keyword'
+        'scanfile':err = 'attribute units of time has not the good format: [^ ]* since ([0-9]){1,4}-([0-9]){1,2}-([0-9]){1,2}.*!C we create a fake calendar ...'
+      ENDCASE
+      IF keyword_set(tobeclosed) THEN ncdf_close, cdfid
+      return, -jpt
+    ENDIF
+    start = str_sep(words[2], '-')
+  ENDIF ELSE units = value
 ;----------------------------------------------------
 ; compute time axis
@@ -185 +186 @@
   ncdf_varget, cdfid, timeid, time
   time = double(time)
-  case units of
+  case units OF
+    'true julian day':
     'second':time = julday(start[1], start[2], start[0], 0, 0, 0)+time/86400.d
     'hour':time = julday(start[1], start[2], start[0], 0, 0, 0)+time/24.d
@@ -191 +193 @@
     'month':BEGIN
       if total(fix(time) NE time) NE 0 then $ ; we switch to days with 30d/m
-         time = julday(start[1], start[2], start[0])+round(time*30) $
-      ELSE for t = 0, n_elements(time)-1 DO $
-         time[t] = julday(start[1]+time[t], start[2], start[0])
+         time = julday(start[1], start[2], start[0], 0, 0, 0) + round(time*30) $
+      ELSE time = julday(start[1]+fix(time), replicate(start[2], jpt), replicate(start[0], jpt), 0, 0, 0)
     END
     'year':BEGIN
       if total(fix(time) NE time) NE 0 then $ ; we switch to days with 365d/y
-         time = julday(start[1], start[2], start[0])+round(time*365) $
-      ELSE for t = 0, n_elements(time)-1 do $
-         time[t] = julday(start[1], start[2], start[0]+time[t])
+         time = julday(start[1], start[2], start[0], 0, 0, 0) + round(time*365) $
+      ELSE time = julday(replicate(start[1], jpt), replicate(start[2], jpt), start[0]+fix(time), 0, 0, 0)
     END
   ENDCASE