Changeset 538 in lmdz_wrf for trunk/tools

Timestamp:

Jun 29, 2015, 7:43:26 PM (10 years ago)

Author:

lfita

Message:

Adding the right '360dqy' calendar calculation

File:

• trunk/tools/nc_var_tools.py (modified) (6 diffs)

trunk/tools/nc_var_tools.py

@@ -2958 +2958 @@
     return timeinf
 
-def date_juliandate(rdate):
-    """ Function to transform from a standard date to a matrix julian day based one
-      rdate= `datetime' object to transform
+def date_juliandate(refyr, TOTsecs):
+    """ Function to transform from a total quantity of seconds since the 
+      beginning of a year to 360 days / year calendar
+      TOTsecs= total number of seconds
     """
-    import datetime as dt
-
     fname = 'date_juliandate'
 
-    jday = int(rdate.timetuple().tm_yday)
-    rldate[0] = int(rdate.year)
-    rldate[1] = jday/30 + 1
-    rldate[2] = jday - (rldate[1]-1)*30
-    rldate[3] = int(rdate.hour)
-    rldate[4] = int(rdate.second)
-    rldate[5] = int(rdate.minute)
-
-    if rldate[2] < 0: 
-        print errormsg
-        print '  ' + fname + ': wrong date', rdate, 'it gives a negative',           \
-          rldate[2],'day!!'
+    secsYear = 3600*24*30*12.
+    secsMonth = 3600*24*30.
+    secsDay = 3600*24.
+    secsHour = 3600.
+    secsMinute = 60.
+
+    rldate = np.zeros((6), dtype=int)
+    rldate[0] = refyr 
+    if TOTsecs > secsYear:
+        rldate[0] = refyr + int(TOTsecs/secsYear) + 1
+        TOTsecs = TOTsecs - int(TOTsecs/secsYear)*secsYear
+
+    if np.mod(TOTsecs,secsMonth) == 0:
+        rldate[1] = int(TOTsecs/secsMonth)
+    else:
+        rldate[1] = int(TOTsecs/secsMonth) + 1
+    TOTsecs = TOTsecs - int(TOTsecs/secsMonth)*secsMonth
+
+    if np.mod(TOTsecs,secsDay) == 0:
+        rldate[2] = int(TOTsecs/secsDay)
+    else:
+        rldate[2] = int(TOTsecs/secsDay) + 1
+    TOTsecs = TOTsecs - int(TOTsecs/secsDay)*secsDay
+
+    if np.mod(TOTsecs,secsHour) == 0:
+        rldate[3] = int(TOTsecs/secsHour)
+    else:
+        rldate[3] = int(TOTsecs/secsHour) + 1
+    TOTsecs = TOTsecs - int(TOTsecs/secsHour)*secsHour
+
+    if np.mod(TOTsecs,secsMinute) == 0:
+        rldate[4] = int(TOTsecs/secsMinute)
+    else:
+        rldate[4] = int(TOTsecs/secsMinute) + 1
+    TOTsecs = TOTsecs - int(TOTsecs/secsMinute)*secsMinute
+
+    rldate[5] = TOTsecs
+
+#    print refyr,TOTsecs,':',rldate
+#    quit()
 
     return rldate
@@ -2994 +3021 @@
     """
     import datetime as dt
+    fname = 'CFtimes_datetime'
 
     times = ncfile.variables[tname]
@@ -3001 +3029 @@
     if not searchInlist(attvar, 'units'):
         print errormsg
-        print '    CFtimes_datetime: "time" does not have attribute: "units"'
+        print '  ' + fname + ": '" + tname + "' does not have attribute: 'units'"
         quit(-1)
     else:
@@ -3022 +3050 @@
     realdates = np.zeros((dimt, 6), dtype=int)
 
+    secsDay=3600*24.
+
 # Checking calendar! 
 ## 
     y360 = False
-    daycal360 = ['earth_360d', '360d', '360days']
+    daycal360 = ['earth_360d', '360d', '360days', '360_day']
     if searchInlist(attvar, 'calendar'):
         calendar = times.getncattr('calendar')
@@ -3045 +3075 @@
         if tunits == 'weeks':
             for it in range(dimt):
-                realdate = refdate + dt.timedelta(weeks=float(times[it]))
-                realdates[it,:] = date_juliandate(realdate)
+                deltat = dt.timedelta(weeks=float(times[it]))
+                Tsecs = deltat.days*secsDay + deltat.seconds + deltat.microseconds/1000.
+                realdates[it,:] = date_juliandate(refdate.year,Tsecs)
         elif tunits == 'days':
             for it in range(dimt):
-                realdate = refdate + dt.timedelta(days=float(times[it]))
-                realdates[it,:] = date_juliandate(realdate)
+                deltat = dt.timedelta(days=float(times[it]))
+                Tsecs = deltat.days*secsDay + deltat.seconds + deltat.microseconds/1000.
+                realdates[it,:] = date_juliandate(refdate.year,Tsecs)
         elif tunits == 'hours':
            for it in range(dimt):
-                realdate = refdate + dt.timedelta(hours=float(times[it]))
-                realdates[it,:] = date_juliandate(realdate)
+                realdate = dt.timedelta(hours=float(times[it]))
+                Tsecs = deltat.days*secsDay + deltat.seconds + deltat.microseconds/1000.
+                realdates[it,:] = date_juliandate(refdate.year,Tsecs)
         elif tunits == 'minutes':
            for it in range(dimt):
-                realdate = refdate + dt.timedelta(minutes=float(times[it]))
-                realdates[it,:] = date_juliandate(realdate)
+                realdate = dt.timedelta(minutes=float(times[it]))
+                Tsecs = deltat.days*secsDay + deltat.seconds + deltat.microseconds/1000.
+                realdates[it,:] = date_juliandate(refdate.year,Tsecs)
         elif tunits == 'seconds':
            for it in range(dimt):
-                realdate = refdate + dt.timedelta(seconds=float(times[it]))
-                realdates[it,:] = date_juliandate(realdate)
+                realdate = dt.timedelta(seconds=float(times[it]))
+                Tsecs = deltat.days*secsDay + deltat.seconds + deltat.microseconds/1000.
+                realdates[it,:] = date_juliandate(refdate.year,Tsecs)
         elif tunits == 'miliseconds':
            for it in range(dimt):
-                realdate = refdate + dt.timedelta(miliseconds=float(times[it]))
-                realdates[it,:] = date_juliandate(realdate)
+                realdate = dt.timedelta(miliseconds=float(times[it]))
+                Tsecs = deltat.days*secsDay + deltat.seconds + deltat.microseconds/1000.
+                realdates[it,:] = date_juliandate(refdate.year,Tsecs)
         else:
               print errormsg
@@ -16202 +16238 @@
                         slicevar.append(ddt)
 
-#                print ddt,d1,d2,mattimes[ddt,:],'var:',slicevar,'new:',slicenewvar
+                print ddt,d1,d2,mattimes[ddt,:],'var:',slicevar,'new:',slicenewvar
 
                 newvarv[tuple(slicenewvar)] = oVn[tuple(slicevar)]