/[lmdze]/trunk/IOIPSL/Calendar/calendar.f
ViewVC logotype

Annotation of /trunk/IOIPSL/Calendar/calendar.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 62 - (hide annotations)
Thu Jul 26 14:37:37 2012 UTC (11 years, 10 months ago) by guez
Original Path: trunk/libf/IOIPSL/calendar.f90
File size: 16004 byte(s)
Changed handling of compiler in compilation system.

Removed the prefix letters "y", "p", "t" or "z" in some names of variables.

Replaced calls to NetCDF by calls to NetCDF95.

Extracted "ioget_calendar" procedures from "calendar.f90" into a
separate file.

Extracted to a separate file, "mathop2.f90", procedures that were not
part of the generic interface "mathop" in "mathop.f90".

Removed computation of "dq" in "bilan_dyn", which was not used.

In "iniadvtrac", removed schemes 20 Slopes and 30 Prather. Was not
compatible with declarations of array sizes.

In "clcdrag", "ustarhb", "vdif_kcay", "yamada4" and "coefkz", changed
the size of some arrays from "klon" to "knon".

Removed possible call to "conema3" in "physiq".

Removed unused argument "cd" in "yamada".

1 guez 30 MODULE calendar
2 guez 36
3     ! From IOIPSL/src/calendar.f90, version 2.0 2004/04/05 14:47:47
4    
5 guez 62 ! This is the calendar used to do all calculations on time. Three
6     ! types of calendars are possible:
7 guez 36
8 guez 62 ! - Gregorian:
9     ! The normal calendar. The time origin for the julian day in this
10     ! case is 24 Nov -4713.
11    
12     ! - No leap:
13     ! A 365 day year without leap years. The origin for the julian days
14     ! is in this case 1 Jan 0.
15    
16     ! - xxxd:
17     ! Year of xxx days with months of equal length. The origin for the
18     ! julian days is then also 1 Jan 0.
19    
20     ! As one can see it is difficult to go from one calendar to the
21     ! other. All operations involving julian days will be wrong. This
22     ! calendar will lock as soon as possible the length of the year and
23     ! forbid any further modification.
24    
25     ! For the non leap-year calendar the method is still brute force.
26     ! We need to find an integer series which takes care of the length
27     ! of the various month. (Jan)
28    
29     USE strlowercase_m, ONLY: strlowercase
30     USE errioipsl, ONLY: histerr
31    
32     IMPLICIT NONE
33    
34 guez 30 PRIVATE
35 guez 62 PUBLIC ymds2ju, ju2ymds, isittime, ioconf_calendar, itau2date, lock_unan, &
36     calendar_used, un_an, un_jour
37 guez 30
38 guez 62 REAL, PARAMETER:: un_jour = 86400. ! one day in seconds
39     LOGICAL:: lock_startdate = .FALSE.
40    
41     CHARACTER(LEN=30):: time_stamp = 'XXXXXXXXXXXXXXXX'
42    
43     ! Description of calendar
44    
45     CHARACTER(LEN=20):: calendar_used = "gregorian"
46     LOGICAL:: lock_unan = .FALSE.
47     REAL:: un_an = 365.2425 ! one year in days
48     INTEGER:: mon_len(12) = (/31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31/)
49    
50     CHARACTER(LEN=3), PARAMETER:: cal(12) = (/'JAN', 'FEB', 'MAR', 'APR', &
51     'MAY', 'JUN', 'JUL', 'AUG', 'SEP', 'OCT', 'NOV', 'DEC'/)
52    
53     REAL, SAVE:: start_day, start_sec
54    
55 guez 30 CONTAINS
56    
57 guez 62 SUBROUTINE ymds2ju (year, month, day, sec, julian)
58 guez 30
59 guez 62 INTEGER, INTENT(IN):: year, month, day
60     REAL, INTENT(IN):: sec
61     REAL, INTENT(OUT):: julian
62 guez 30
63 guez 62 INTEGER:: julian_day
64     REAL:: julian_sec
65 guez 30
66 guez 62 !--------------------------------------------------------------------
67 guez 30
68 guez 62 CALL ymds2ju_internal(year, month, day, sec, julian_day, julian_sec)
69 guez 36 julian = julian_day + julian_sec / un_jour
70 guez 62
71 guez 30 END SUBROUTINE ymds2ju
72    
73     !===
74    
75 guez 62 SUBROUTINE ymds2ju_internal (year, month, day, sec, julian_day, julian_sec)
76 guez 30
77 guez 62 ! Converts year, month, day and seconds into a julian day
78 guez 30
79 guez 62 ! In 1968 in a letter to the editor of Communications of the ACM
80     ! (CACM, volume 11, number 10, October 1968, p.657) Henry F. Fliegel
81     ! and Thomas C. Van Flandern presented such an algorithm.
82 guez 30
83 guez 62 ! See also: http://www.magnet.ch/serendipity/hermetic/cal_stud/jdn.htm
84 guez 30
85 guez 62 ! In the case of the Gregorian calendar we have chosen to use
86     ! the Lilian day numbers. This is the day counter which starts
87     ! on the 15th October 1582.
88     ! This is the day at which Pope Gregory XIII introduced the
89     ! Gregorian calendar.
90     ! Compared to the true Julian calendar, which starts some
91     ! 7980 years ago, the Lilian days are smaler and are dealt with
92     ! easily on 32 bit machines. With the true Julian days you can only
93     ! the fraction of the day in the real part to a precision of
94     ! a 1/4 of a day with 32 bits.
95 guez 30
96 guez 62 INTEGER, INTENT(IN):: year, month, day
97     REAL, INTENT(IN):: sec
98 guez 30
99 guez 62 INTEGER, INTENT(OUT):: julian_day
100     REAL, INTENT(OUT):: julian_sec
101    
102     INTEGER:: jd, m, y, d, ml
103     !--------------------------------------------------------------------
104 guez 30 lock_unan = .TRUE.
105    
106     m = month
107     y = year
108     d = day
109    
110 guez 62 ! We deduce the calendar from the length of the year as it
111     ! is faster than an INDEX on the calendar variable.
112 guez 30
113 guez 62 ! Gregorian
114 guez 30 IF ( (un_an > 365.0).AND.(un_an < 366.0) ) THEN
115     jd = (1461*(y+4800+INT(( m-14 )/12)))/4 &
116     & +(367*(m-2-12*(INT(( m-14 )/12))))/12 &
117     & -(3*((y+4900+INT((m-14)/12))/100))/4 &
118     & +d-32075
119     jd = jd-2299160
120 guez 62 ! No leap or All leap
121 guez 30 ELSE IF (ABS(un_an-365.0) <= EPSILON(un_an) .OR. &
122     & ABS(un_an-366.0) <= EPSILON(un_an)) THEN
123     ml = SUM(mon_len(1:m-1))
124     jd = y*INT(un_an)+ml+(d-1)
125 guez 62 ! Calendar with regular month
126 guez 30 ELSE
127     ml = INT(un_an)/12
128     jd = y*INT(un_an)+(m-1)*ml+(d-1)
129     ENDIF
130    
131     julian_day = jd
132     julian_sec = sec
133 guez 62
134 guez 30 END SUBROUTINE ymds2ju_internal
135 guez 62
136 guez 30 !===
137    
138 guez 62 SUBROUTINE ju2ymds (julian, year, month, day, sec)
139 guez 30
140 guez 62 REAL, INTENT(IN):: julian
141 guez 30
142 guez 62 INTEGER, INTENT(OUT):: year, month, day
143     REAL, INTENT(OUT):: sec
144    
145     INTEGER:: julian_day
146     REAL:: julian_sec
147     !--------------------------------------------------------------------
148 guez 30 julian_day = INT(julian)
149     julian_sec = (julian-julian_day)*un_jour
150    
151 guez 62 CALL ju2ymds_internal(julian_day, julian_sec, year, month, day, sec)
152    
153 guez 30 END SUBROUTINE ju2ymds
154 guez 62
155 guez 30 !===
156    
157 guez 62 SUBROUTINE ju2ymds_internal (julian_day, julian_sec, year, month, day, sec)
158 guez 30
159 guez 62 ! This subroutine computes from the julian day the year,
160     ! month, day and seconds
161 guez 30
162 guez 62 ! In 1968 in a letter to the editor of Communications of the ACM
163     ! (CACM, volume 11, number 10, October 1968, p.657) Henry F. Fliegel
164     ! and Thomas C. Van Flandern presented such an algorithm.
165 guez 30
166 guez 62 ! See also: http://www.magnet.ch/serendipity/hermetic/cal_stud/jdn.htm
167 guez 30
168 guez 62 ! In the case of the Gregorian calendar we have chosen to use
169     ! the Lilian day numbers. This is the day counter which starts
170     ! on the 15th October 1582. This is the day at which Pope
171     ! Gregory XIII introduced the Gregorian calendar.
172     ! Compared to the true Julian calendar, which starts some 7980
173     ! years ago, the Lilian days are smaler and are dealt with easily
174     ! on 32 bit machines. With the true Julian days you can only the
175     ! fraction of the day in the real part to a precision of a 1/4 of
176     ! a day with 32 bits.
177 guez 30
178 guez 62 INTEGER, INTENT(IN):: julian_day
179     REAL, INTENT(IN):: julian_sec
180    
181     INTEGER, INTENT(OUT):: year, month, day
182     REAL, INTENT(OUT):: sec
183    
184     INTEGER:: l, n, i, jd, j, d, m, y, ml
185     INTEGER:: add_day
186     !--------------------------------------------------------------------
187 guez 30 lock_unan = .TRUE.
188    
189     jd = julian_day
190     sec = julian_sec
191     IF (sec > un_jour) THEN
192     add_day = INT(sec/un_jour)
193     sec = sec-add_day*un_jour
194     jd = jd+add_day
195     ENDIF
196    
197 guez 62 ! Gregorian
198 guez 30 IF ( (un_an > 365.0).AND.(un_an < 366.0) ) THEN
199     jd = jd+2299160
200    
201     l = jd+68569
202     n = (4*l)/146097
203     l = l-(146097*n+3)/4
204     i = (4000*(l+1))/1461001
205     l = l-(1461*i)/4+31
206     j = (80*l)/2447
207     d = l-(2447*j)/80
208     l = j/11
209     m = j+2-(12*l)
210     y = 100*(n-49)+i+l
211 guez 62 ! No leap or All leap
212 guez 30 ELSE IF (ABS(un_an-365.0) <= EPSILON(un_an) .OR. &
213     & ABS(un_an-366.0) <= EPSILON(un_an) ) THEN
214     y = jd/INT(un_an)
215     l = jd-y*INT(un_an)
216     m = 1
217     ml = 0
218     DO WHILE (ml+mon_len(m) <= l)
219     ml = ml+mon_len(m)
220     m = m+1
221     ENDDO
222     d = l-ml+1
223 guez 62 ! others
224 guez 30 ELSE
225     ml = INT(un_an)/12
226     y = jd/INT(un_an)
227     l = jd-y*INT(un_an)
228     m = (l/ml)+1
229     d = l-(m-1)*ml+1
230     ENDIF
231    
232     day = d
233     month = m
234     year = y
235 guez 62
236 guez 30 END SUBROUTINE ju2ymds_internal
237 guez 62
238 guez 30 !===
239    
240 guez 62 REAL FUNCTION itau2date (itau, date0, deltat)
241 guez 30
242 guez 62 ! This function transforms itau into a date. The date whith which
243     ! the time axis is going to be labeled
244 guez 30
245 guez 62 ! INPUT
246     ! itau: current time step
247     ! date0: Date at which itau was equal to 0
248     ! deltat: time step between itau s
249    
250     ! OUTPUT
251     ! itau2date: Date for the given itau
252    
253     INTEGER:: itau
254     REAL:: date0, deltat
255     !--------------------------------------------------------------------
256 guez 30 itau2date = REAL(itau)*deltat/un_jour+date0
257 guez 62
258 guez 30 END FUNCTION itau2date
259 guez 62
260 guez 30 !===
261 guez 62
262 guez 30 SUBROUTINE isittime &
263 guez 62 & (itau, date0, dt, freq, last_action, last_check, do_action)
264 guez 30
265 guez 62 ! This subroutine checks the time has come for a given action.
266     ! This is computed from the current time-step(itau).
267     ! Thus we need to have the time delta (dt), the frequency
268     ! of the action (freq) and the last time it was done
269     ! (last_action in units of itau).
270     ! In order to extrapolate when will be the next check we need
271     ! the time step of the last call (last_check).
272 guez 30
273 guez 62 ! The test is done on the following condition:
274     ! the distance from the current time to the time for the next
275     ! action is smaller than the one from the next expected
276     ! check to the next action.
277     ! When the test is done on the time steps simplifactions make
278     ! it more difficult to read in the code.
279     ! For the real time case it is easier to understand !
280 guez 30
281 guez 62 INTEGER, INTENT(IN):: itau
282     REAL, INTENT(IN):: dt, freq
283     INTEGER, INTENT(IN):: last_action, last_check
284     REAL, INTENT(IN):: date0
285 guez 30
286 guez 62 LOGICAL, INTENT(OUT):: do_action
287    
288     REAL:: dt_action, dt_check
289     REAL:: date_last_act, date_next_check, date_next_act, &
290     & date_now, date_mp1, date_mpf
291     INTEGER:: year, month, monthp1, day, next_check_itau, next_act_itau
292     INTEGER:: yearp, dayp
293     REAL:: sec, secp
294     LOGICAL:: check = .FALSE.
295     !--------------------------------------------------------------------
296 guez 30 IF (check) THEN
297 guez 62 WRITE(*, *) &
298     & "isittime 1.0 ", itau, date0, dt, freq, last_action, last_check
299 guez 30 ENDIF
300    
301     IF (last_check >= 0) THEN
302     dt_action = (itau-last_action)*dt
303     dt_check = (itau-last_check)*dt
304     next_check_itau = itau+(itau-last_check)
305    
306 guez 62 !- We are dealing with frequencies in seconds and thus operation
307     !- can be done on the time steps.
308 guez 30
309     IF (freq > 0) THEN
310     IF (ABS(dt_action-freq) <= ABS(dt_action+dt_check-freq)) THEN
311     do_action = .TRUE.
312     ELSE
313     do_action = .FALSE.
314     ENDIF
315    
316 guez 62 !--- Here we deal with frequencies in month and work on julian days.
317 guez 30
318     ELSE
319 guez 62 date_now = itau2date (itau, date0, dt)
320     date_last_act = itau2date (last_action, date0, dt)
321     CALL ju2ymds (date_last_act, year, month, day, sec)
322 guez 35 monthp1 = month - freq
323 guez 30 yearp = year
324    
325 guez 62 !--- Here we compute what logically should be the next month
326 guez 30
327     IF (month >= 13) THEN
328     yearp = year+1
329     monthp1 = monthp1-12
330     ENDIF
331 guez 62 CALL ymds2ju (year, monthp1, day, sec, date_mpf)
332 guez 30
333 guez 62 !--- But it could be that because of a shorter month or a bad
334     !--- starting date that we end up further than we should be.
335     !--- Thus we compute the first day of the next month.
336     !--- We can not be beyond this date and if we are close
337     !--- then we will take it as it is better.
338 guez 30
339     monthp1 = month+ABS(freq)
340     yearp=year
341     IF (monthp1 >= 13) THEN
342     yearp = year+1
343     monthp1 = monthp1 -12
344     ENDIF
345     dayp = 1
346     secp = 0.0
347 guez 62 CALL ymds2ju (yearp, monthp1, dayp, secp, date_mp1)
348 guez 30
349 guez 62 !--- If date_mp1 is smaller than date_mpf or only less than 4 days
350     !--- larger then we take it. This needed to ensure that short month
351     !--- like February do not mess up the thing !
352 guez 30
353     IF (date_mp1-date_mpf < 4.) THEN
354     date_next_act = date_mp1
355     ELSE
356     date_next_act = date_mpf
357     ENDIF
358 guez 62 date_next_check = itau2date (next_check_itau, date0, dt)
359 guez 30
360 guez 62 !--- Transform the dates into time-steps for the needed precisions.
361 guez 30
362     next_act_itau = &
363     & last_action+INT((date_next_act-date_last_act)*(un_jour/dt))
364 guez 62
365 guez 30 IF ( ABS(itau-next_act_itau) &
366     & <= ABS( next_check_itau-next_act_itau)) THEN
367     do_action = .TRUE.
368     IF (check) THEN
369 guez 62 WRITE(*, *) &
370     & 'ACT-TIME: itau, next_act_itau, next_check_itau: ', &
371     & itau, next_act_itau, next_check_itau
372     CALL ju2ymds (date_now, year, month, day, sec)
373     WRITE(*, *) 'ACT-TIME: y, m, d, s: ', year, month, day, sec
374     WRITE(*, *) &
375     & 'ACT-TIME: date_mp1, date_mpf: ', date_mp1, date_mpf
376 guez 30 ENDIF
377     ELSE
378     do_action = .FALSE.
379     ENDIF
380     ENDIF
381    
382     IF (check) THEN
383 guez 62 WRITE(*, *) "isittime 2.0 ", &
384     & date_next_check, date_next_act, ABS(dt_action-freq), &
385     & ABS(dt_action+dt_check-freq), dt_action, dt_check, &
386     & next_check_itau, do_action
387 guez 30 ENDIF
388     ELSE
389     do_action=.FALSE.
390     ENDIF
391 guez 62
392 guez 30 END SUBROUTINE isittime
393 guez 62
394 guez 30 !===
395 guez 62
396 guez 30 SUBROUTINE ioconf_calendar (str)
397    
398 guez 62 ! This routine allows to configure the calendar to be used.
399     ! This operation is only allowed once and the first call to
400     ! ymds2ju or ju2ymsd will lock the current configuration.
401     ! the argument to ioconf_calendar can be any of the following:
402     ! - gregorian: This is the gregorian calendar (default here)
403     ! - noleap: A calendar without leap years = 365 days
404     ! - xxxd: A calendar of xxx days (has to be a modulo of 12)
405     ! with 12 month of equal length
406 guez 30
407 guez 62 CHARACTER(LEN=*), INTENT(IN):: str
408 guez 30
409 guez 62 INTEGER:: leng, ipos
410     CHARACTER(LEN=10):: str10
411     !--------------------------------------------------------------------
412    
413 guez 30 ! 1.0 Clean up the sring !
414    
415     CALL strlowercase (str)
416    
417     IF (.NOT.lock_unan) THEN
418 guez 62
419 guez 30 lock_unan=.TRUE.
420 guez 62
421 guez 30 SELECT CASE(str)
422     CASE('gregorian')
423     calendar_used = 'gregorian'
424     un_an = 365.2425
425 guez 62 mon_len(:)=(/31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31/)
426 guez 30 CASE('standard')
427     calendar_used = 'gregorian'
428     un_an = 365.2425
429 guez 62 mon_len(:)=(/31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31/)
430 guez 30 CASE('proleptic_gregorian')
431     calendar_used = 'gregorian'
432     un_an = 365.2425
433 guez 62 mon_len(:)=(/31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31/)
434 guez 30 CASE('noleap')
435     calendar_used = 'noleap'
436     un_an = 365.0
437 guez 62 mon_len(:)=(/31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31/)
438 guez 30 CASE('365_day')
439     calendar_used = 'noleap'
440     un_an = 365.0
441 guez 62 mon_len(:)=(/31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31/)
442 guez 30 CASE('365d')
443     calendar_used = 'noleap'
444     un_an = 365.0
445 guez 62 mon_len(:)=(/31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31/)
446 guez 30 CASE('all_leap')
447     calendar_used = 'all_leap'
448     un_an = 366.0
449 guez 62 mon_len(:)=(/31, 29, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31/)
450 guez 30 CASE('366_day')
451     calendar_used = 'all_leap'
452     un_an = 366.0
453 guez 62 mon_len(:)=(/31, 29, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31/)
454 guez 30 CASE('366d')
455     calendar_used = 'all_leap'
456     un_an = 366.0
457 guez 62 mon_len(:)=(/31, 29, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31/)
458 guez 30 CASE DEFAULT
459 guez 62 ipos = INDEX(str, 'd')
460 guez 30 IF (ipos == 4) THEN
461 guez 62 READ(str(1:3), '(I3)') leng
462     IF ( (MOD(leng, 12) == 0).AND.(leng > 1) ) THEN
463 guez 30 calendar_used = str
464     un_an = leng
465     mon_len(:) = leng
466     ELSE
467 guez 62 CALL histerr (3, 'ioconf_calendar', &
468 guez 30 & 'The length of the year as to be a modulo of 12', &
469     & 'so that it can be divided into 12 month of equal length', &
470     & str)
471     ENDIF
472     ELSE
473 guez 62 CALL histerr (3, 'ioconf_calendar', &
474     & 'Unrecognized input, please ceck the man pages.', str, ' ')
475 guez 30 ENDIF
476     END SELECT
477     ELSE
478 guez 62 WRITE(str10, '(f10.4)') un_an
479     CALL histerr (2, 'ioconf_calendar', &
480 guez 30 & 'The calendar was already used or configured. You are not', &
481     & 'allowed to change it again. '// &
482 guez 62 & 'The following length of year is used:', str10)
483 guez 30 ENDIF
484 guez 62
485 guez 30 END SUBROUTINE ioconf_calendar
486    
487     END MODULE calendar

  ViewVC Help
Powered by ViewVC 1.1.21