source: palm/trunk/SOURCE/date_and_time_mod.f90 @ 4016

Last change on this file since 4016 was 3839, checked in by moh.hefny, 6 years ago

fixing negative day situation when 1st day of month is used with spinup

  • Property svn:keywords set to Id
File size: 21.3 KB
Line 
1!> @file date_and_time_mod.f90
2!------------------------------------------------------------------------------!
3! This file is part of the PALM model system.
4!
5! PALM is free software: you can redistribute it and/or modify it under the
6! terms of the GNU General Public License as published by the Free Software
7! Foundation, either version 3 of the License, or (at your option) any later
8! version.
9!
10! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
13!
14! You should have received a copy of the GNU General Public License along with
15! PALM. If not, see <http://www.gnu.org/licenses/>.
16!
17! Copyright 1997-2019 Leibniz Universitaet Hannover
18!------------------------------------------------------------------------------!
19!
20! Current revisions:
21! ------------------
22!
23!
24! Former revisions:
25! -----------------
26! $Id: date_and_time_mod.f90 3839 2019-03-28 21:12:25Z forkel $
27! further tabs removed, unused variables removed
28!
29! 3655 2019-01-07 16:51:22Z knoop
30! further tabs removed
31!
32! 3467 2018-10-30 19:05:21Z suehring
33! Tabs removed
34!
35! 3458 2018-10-30 14:51:23Z kanani
36! from chemistry branch r3443, banzhafs:
37! Added initial hour_of_day, hour_of_year, day_of_year and month_of_year to
38! init_date_and_time
39!
40! 3298 2018-10-02 12:21:11Z kanani
41! - Minor formatting (kanani)
42! - Added Routines for DEFAULT mode of chemistry emissions (Russo)
43! - Added routine for reading-in real dates: format has to be in DDMMYYYY and
44!   passed in the namelist parameter date_init (Russo)
45! - Added calculation of several time indices useful in several routines
46!   of the model (Russo)
47!
48! 2718 2018-01-02 08:49:38Z maronga
49! Corrected "Former revisions" section
50!
51! 2701 2017-12-15 15:40:50Z suehring
52! Changes from last commit documented
53!
54! 2698 2017-12-14 18:46:24Z suehring
55! Bugfix in definition of d_seconds_year.
56!
57! 2696 2017-12-14 17:12:51Z kanani
58! Change in file header (GPL part)
59!
60! 2544 2017-10-13 18:09:32Z maronga
61! Initial revision
62!
63!
64
65!
66! Description:
67! ------------
68!> This routine calculates all needed information on date and time used by
69!> other modules
70!> @todo Further testing and revision of routines for updating indices of
71!>       emissions in the default mode
72!> @todo Add routine for recognizing leap years
73!> @todo Add recognition of exact days of week (Monday, Tuesday, etc.)
74!> @todo Reconsider whether to remove day_of_year_init from the namelist: we
75!>       already implemented changes for calculating it from date_init in
76!>       calc_date_and_time
77!> @todo time_utc during spin-up 
78!------------------------------------------------------------------------------!
79 MODULE date_and_time_mod
80
81    USE control_parameters,                                                   &
82        ONLY:  coupling_start_time, days_since_reference_point,               &
83               message_string, simulated_time, time_since_reference_point
84
85    USE kinds
86
87
88    IMPLICIT NONE
89
90    PRIVATE
91
92!-- Variables Declaration
93
94    INTEGER(iwp)        ::  day_of_year      = 0   !< day of the year (DOY)
95    INTEGER(iwp)        ::  day_of_year_init = 0   !< DOY at model start (default: 0)
96
97    ! --- Most of these indices are updated by the routine calc_date_and_time according to the current date and time of the simulation
98    INTEGER(iwp)        ::  hour_of_year = 1                        !< hour of the current year (1:8760(8784))
99    INTEGER(iwp)        ::  hour_of_day=1                           !< hour of the current day (1:24)
100    INTEGER(iwp)        ::  day_of_month=0                          !< day of the current month (1:31)
101    INTEGER(iwp)        ::  month_of_year=0                         !< month of the current year (1:12)
102    INTEGER(iwp)        ::  current_year=0                          !< current year
103    INTEGER(iwp)        ::  hour_call_emis=0                        !< index used to call the emissions just once every hour
104
105    INTEGER(iwp)        ::  index_mm                                !< index months of the default emission mode
106    INTEGER(iwp)        ::  index_dd                                !< index days of the default emission mode
107    INTEGER(iwp)        ::  index_hh                                !< index hours of the emission mode
108
109    REAL(wp)            ::  time_utc                     !< current model time in UTC
110    REAL(wp)            ::  time_utc_emis                !< current emission module time in UTC
111    REAL(wp)            ::  time_utc_init = 43200.0_wp   !< UTC time at model start
112    REAL(wp)            ::  time_update                  !< used to calculate actual second of the simulation
113 
114    REAL(wp), PARAMETER ::  d_hours_day    = 1.0_wp / 24.0_wp       !< inverse of hours per day (1/24)
115    REAL(wp), PARAMETER ::  d_seconds_hour = 1.0_wp / 3600.0_wp     !< inverse of seconds per hour (1/3600)
116    REAL(wp), PARAMETER ::  d_seconds_year = 1.0_wp / 31536000.0_wp !< inverse of the seconds per year (1/(365*86400))
117   
118    CHARACTER(len=8)    ::  date_init = "21062017"                  !< Starting date of simulation: We selected this because it was a monday
119 
120    !> --- Parameters
121    INTEGER, PARAMETER, DIMENSION(12) :: days = (/31,28,31,30,31,30,31,31,30,31,30,31/) ! total number of days for each month (no leap year)
122
123    SAVE
124
125!-- INTERFACES PART
126    !-- Read initial day and time of simulation
127    INTERFACE init_date_and_time
128       MODULE PROCEDURE init_date_and_time
129    END INTERFACE init_date_and_time
130
131    !-- Get hour index in the DEAFULT case of chemistry emissions :
132    INTERFACE time_default_indices
133       MODULE PROCEDURE time_mdh_indices
134       MODULE PROCEDURE time_hour_indices
135    END INTERFACE time_default_indices
136
137    !-- Get hour index in the PRE-PROCESSED case of chemistry emissions :
138    INTERFACE time_preprocessed_indices
139       MODULE PROCEDURE time_preprocessed_indices
140    END INTERFACE time_preprocessed_indices
141
142
143    !-- Calculate current date and time
144    INTERFACE calc_date_and_time
145       MODULE PROCEDURE calc_date_and_time
146    END INTERFACE
147
148
149    !-- Public Interfaces
150    PUBLIC calc_date_and_time, time_default_indices, init_date_and_time, time_preprocessed_indices
151
152    !-- Public Variables
153    PUBLIC date_init, d_hours_day, d_seconds_hour, d_seconds_year,               &
154           day_of_year, day_of_year_init, time_utc, time_utc_init, day_of_month, &
155           month_of_year, index_mm, index_dd, index_hh, hour_of_day, hour_of_year, &
156           hour_call_emis
157
158 CONTAINS
159
160
161 !------------------------------------------------------------------------------!
162 !> Reads starting date from namelist
163 !------------------------------------------------------------------------------!
164 
165    SUBROUTINE init_date_and_time
166
167       IMPLICIT NONE
168
169       !--    Variables Definition
170       INTEGER ::  i_mon       !< Index for going through the different months
171
172       IF  (day_of_year_init == 0) THEN
173          ! Day of the month at starting time
174          READ(UNIT=date_init(1:2),fmt=*)day_of_month
175
176          ! Month of the year at starting time
177          READ(UNIT=date_init(3:4),fmt=*)month_of_year
178
179          ! Year at starting time
180          READ(UNIT=date_init(5:8),fmt=*)current_year
181       
182       ENDIF
183
184
185       !-- Calculate initial hour of the day: the first hour of the day is from 00:00:00 to 00:59:59.
186
187       hour_of_day = INT( FLOOR( time_utc_init/3600.0_wp ) ) + 1
188
189       !-- Calculate initial day day_of_year_init in case date_init is given or day_of_year_init is given
190       IF ( day_of_year_init == 0 ) THEN
191
192          !> Condition for printing an error when date_init is not provided when day_of_year_init is not given in the namelist or when the format of the date is not the one required by PALM.
193          IF ( day_of_month .GT. 0 .AND. day_of_month .LE. 31 .AND. month_of_year .GT. 0 .AND. month_of_year .LE. 12) THEN
194       
195             IF ( month_of_year == 1 ) THEN  !!month of year is read in input
196
197                day_of_year_init = day_of_month
198
199             ELSE
200
201                day_of_year_init= SUM(days( 1:(month_of_year-1) )) + day_of_month  !day_of_month is read in input in this case
202
203             ENDIF
204!kanani: Revise, we cannot force users to provide date_init, maybe set a default value?
205!           ELSE
206!
207!              message_string = 'date_init not provided in the namelist or'             //          &
208!                               ' given in the wrong format: MUST BE DDMMYYYY'                       
209!              CALL message( 'calc_date_and_time', 'DT0100', 2, 2, 0, 6, 0 )
210     
211          ENDIF
212
213       ENDIF
214
215       
216       !-- Initial day of the year
217       day_of_year = day_of_year_init
218
219       !-- Initial hour of the year
220       hour_of_year = ( (day_of_year-1) * 24 ) + hour_of_day
221
222       !--Initial day of the month and month of the year
223       !> --------------------------------------------------------------------------------
224       !> The first case is when date_init is not provided: we only know day_of_year_init     
225       IF ( month_of_year == 0 .AND. day_of_month == 0) THEN
226
227         
228          IF ( day_of_year .LE. 31 ) THEN
229
230             month_of_year=1
231             day_of_month=day_of_year
232
233          ELSE
234
235             DO i_mon=2,12   !january is considered in the first case
236                IF ( day_of_year .LE. SUM(days(1:i_mon)) .AND. day_of_year .GT. SUM(days(1:(i_mon-1))) ) THEN
237           
238                   month_of_year=i_mon
239
240                   day_of_month=INT(MOD(day_of_year, SUM(days(1:(i_mon-1))))) 
241
242                   GOTO 38
243
244                ENDIF
245
246             38 ENDDO
247          ENDIF
248       !> --------------------------------------------------------------------------------
249       !> in the second condition both day of month and month_of_year are either given in input (passed to date_init) or we are in some day successive to the initial one, so that day_of_month has already be computed in previous step
250       !>TBD: something to calculate the current year is missing
251       ELSEIF ( day_of_month .GT. 0 .AND. day_of_month .LE. 31 .AND. month_of_year .GT. 0 .AND. month_of_year .LE. 12) THEN
252 
253          !> calculate month_of_year. TBD: test the condition when day_of_year==31
254 
255          IF (day_of_year==1) THEN  !> this allows to turn from december to January when passing from a year to another
256 
257             month_of_year = 1
258       
259          ELSE IF (day_of_year .GT. 1 .AND. day_of_year .GT. SUM(days(1:month_of_year))) THEN
260
261             month_of_year = month_of_year + 1
262
263          ENDIF
264
265          !> calculate day_of_month
266          IF ( month_of_year == 1 ) THEN
267           
268            day_of_month=day_of_year
269
270          ELSE
271
272            day_of_month=INT(MOD(day_of_year, SUM(days(1:(month_of_year-1)))))
273
274          ENDIF
275
276
277       ELSE
278
279          !> Condition when date_init is provided but it is given in the wrong format
280          message_string = 'date_init not provided in the namelist or'            //          &
281                              ' given in the wrong format: MUST BE DDMMYYYY'                 
282          CALL message( 'init_date_and_time', 'DT0102', 2, 2, 0, 6, 0 ) 
283
284       ENDIF
285
286
287    END SUBROUTINE init_date_and_time
288
289!------------------------------------------------------------------------------!
290! Description:
291! ------------
292!> Calculate current date and time of the simulation
293!------------------------------------------------------------------------------!
294 
295    SUBROUTINE calc_date_and_time
296
297       IMPLICIT NONE
298
299!--    Variables Definition
300       INTEGER :: i_mon       !< Index for going through the different months
301
302       !> Update simulation time in seconds
303       time_update = simulated_time-coupling_start_time
304
305!--    Calculate current day of the simulated time
306       days_since_reference_point=INT(FLOOR( (time_utc_init + time_update) &
307                               / 86400.0_wp ) )
308
309!--    Calculate actual UTC time                       
310       time_utc = MOD((time_utc_init + time_since_reference_point), 86400.0_wp)
311       
312!sB    PRILIMINARY workaround for time_utc changes due to changes in time_since_reference_point in
313!sB    radiation_model_mod during runtime:
314       time_utc_emis = MOD((time_utc_init + time_update), 86400.0_wp)     
315
316!--    Calculate initial day of the year: it is calculated only once. In fact, day_of_year_init is initialized to 0 and then a positive value is passed. This condition is also called only when day_of_year_init is not given in the namelist.
317
318       IF ( day_of_year_init == 0 ) THEN
319
320          !> Condition for printing an error when date_init is not provided when day_of_year_init is not given in the namelist or when the format of the date is not the one required by PALM.
321          IF ( day_of_month .GT. 0 .AND. day_of_month .LE. 31 .AND. month_of_year .GT. 0 .AND. month_of_year .LE. 12) THEN
322       
323             IF ( month_of_year == 1 ) THEN  !!month of year is read in input
324
325                day_of_year_init = day_of_month
326
327             ELSE
328
329                day_of_year_init= SUM(days( 1:(month_of_year-1) )) + day_of_month  !day_of_month is read in input in this case
330
331             ENDIF
332!kanani: Revise, we cannot force users to provide date_init, maybe set a default value?
333!           ELSE
334!
335!              message_string = 'date_init not provided in the namelist or'             //          &
336!                               ' given in the wrong format: MUST BE DDMMYYYY'                       
337!              CALL message( 'calc_date_and_time', 'DT0100', 2, 2, 0, 6, 0 )
338     
339          ENDIF
340
341       ENDIF
342
343      !-- Calculate actual hour of the day: the first hour of the day is from 00:00:00 to 00:59:59.
344
345      hour_of_day = INT( FLOOR( time_utc_emis/3600.0_wp ) ) + 1
346
347!--    Calculate current day of the year !TBD: considetr leap years
348       IF ( (day_of_year_init + days_since_reference_point)  .GT. 365 ) THEN
349
350          day_of_year=INT(MOD((day_of_year_init + days_since_reference_point), 365.0_wp))
351
352       ELSE
353         
354          day_of_year = day_of_year_init + days_since_reference_point
355
356       ENDIF
357
358!
359!--    Calculate current hour of the year
360       hour_of_year = ( (day_of_year-1) * 24 ) + hour_of_day  !> actual hour of the year
361       
362
363!
364!--    UPDATE actual day of the month and month of the year
365       !> --------------------------------------------------------------------------------
366       !> The first case is when date_init is not provided: we only know day_of_year_init     
367       IF ( month_of_year == 0 .AND. day_of_month == 0) THEN
368
369          !> The first case is when date_init is not provided: we only know day_of_year_init
370          !DO i_mon=1,12
371             !IF (day_of_year .LE. SUM(days(1:i_mon))) THEN
372          IF ( day_of_year .LE. 31 ) THEN
373
374             month_of_year=1
375             day_of_month=day_of_year
376
377          ELSE
378
379             DO i_mon=2,12   !january is considered in the first case
380                IF ( day_of_year .LE. SUM(days(1:i_mon)) .AND. day_of_year .GT. SUM(days(1:(i_mon-1))) ) THEN
381           
382                   month_of_year=i_mon
383
384                   day_of_month=INT(MOD(day_of_year, SUM(days(1:(i_mon-1))))) 
385
386                   GOTO 38
387
388                ENDIF
389
390             38 ENDDO
391          ENDIF
392       !> --------------------------------------------------------------------------------
393       !> in the second condition both day of month and month_of_year are either given in input (passed to date_init) or we are in some day successive to the initial one, so that day_of_month has already be computed in previous step
394       !>TBD: something to calculate the current year is missing
395       ELSEIF ( day_of_month .GT. 0 .AND. day_of_month .LE. 31 .AND. month_of_year .GT. 0 .AND. month_of_year .LE. 12) THEN
396 
397          !> calculate month_of_year. TBD: test the condition when day_of_year==31
398 
399          IF (day_of_year==1) THEN  !> this allows to turn from december to January when passing from a year to another
400 
401             month_of_year = 1
402       
403          ELSE IF (day_of_year .GT. 1 .AND. day_of_year .GT. SUM(days(1:month_of_year))) THEN
404
405             month_of_year = month_of_year + 1
406
407          ENDIF
408
409          !> calculate day_of_month
410          IF ( month_of_year == 1 ) THEN
411           
412            day_of_month=day_of_year
413
414          ELSE
415
416            day_of_month=INT(MOD(day_of_year, SUM(days(1:(month_of_year-1)))))
417
418          ENDIF
419
420          ! fix the date if the day is 1st and earlier day is needed due to spinup
421          IF ( day_of_month < 1 ) THEN
422 
423             ! if the day is the first day in the year
424             IF ( month_of_year  ==  1 ) THEN
425                month_of_year = 12
426                day_of_month = 31
427 
428             ! other cases
429             ELSE
430                month_of_year = month_of_year - 1
431                day_of_month = days(month_of_year)
432             ENDIF
433 
434          ENDIF
435
436      ELSE
437
438          !> Condition when date_init is provided but it is given in the wrong format
439          message_string = 'date_init not provided in the namelist or'            //          &
440                              ' given in the wrong format: MUST BE DDMMYYYY'                 
441          CALL message( 'calc_date_and_time', 'DT0101', 2, 2, 0, 6, 0 ) 
442
443      ENDIF       
444
445    END SUBROUTINE calc_date_and_time
446
447
448!------------------------------------------------------------------------------!
449! Description:
450! ------------
451!> This routine determines the time factor index in the PRE-PROCESSED emissions mode.
452!------------------------------------------------------------------------------!
453
454 SUBROUTINE time_preprocessed_indices(index_hh)
455
456    USE indices
457
458    IMPLICIT NONE
459
460!
461!-- In/output
462    INTEGER, INTENT(INOUT) ::  index_hh    !> Index Hour
463!
464!-- Additional Variables for calculateing indices
465!-- Constants
466    INTEGER, PARAMETER ::  nhour = 24
467
468    IF (days_since_reference_point == 0) THEN
469
470       index_hh=hour_of_day
471
472    ELSE
473
474       index_hh=(days_since_reference_point*nhour)+(hour_of_day)
475
476    ENDIF
477
478
479 END SUBROUTINE time_preprocessed_indices
480
481
482!------------------------------------------------------------------------------!
483! Description:
484! ------------
485!> This routine determines the time factor index in the mdh case of the DEFAULT chemistry emissions mode.
486!------------------------------------------------------------------------------!
487
488 SUBROUTINE time_mdh_indices(daytype_mdh,mo, dd, hh, index_mm, index_dd, index_hh)
489
490    USE indices
491
492    IMPLICIT NONE
493
494    !> IN/OUTPUT
495    INTEGER, INTENT(INOUT) :: mo          !> Month of year
496    INTEGER, INTENT(INOUT) :: dd          !> Day of month
497    INTEGER, INTENT(INOUT) :: hh          !> Hour of day
498    INTEGER, INTENT(INOUT) :: index_mm    !> Index Month
499    INTEGER, INTENT(INOUT) :: index_dd    !> Index Day
500    INTEGER, INTENT(INOUT) :: index_hh    !> Index Hour
501
502    CHARACTER(len=80), INTENT(INOUT) :: daytype_mdh !> type of the day in mdh mode: one of 1-WORKDAY
503                                                    !                                      2-WEEKEND
504                                                    !                                      3-HOLIDAY
505
506    REAL(wp)                         :: frac_day=0 
507
508    !> ------------------------------------------------------------------------
509
510    INTEGER                          :: weekday
511
512    !> CONSTANTS
513    INTEGER, PARAMETER               :: nmonth = 12
514    INTEGER, PARAMETER               :: nday = 7
515    INTEGER, PARAMETER               :: nhour = 24
516
517    frac_day= (dd-1)/nday    !> indicates the week of the month, supposing the month starts on monday
518
519   ! 1:7      1:31  7                   (0:30)/7
520    weekday = dd-( nday * (INT( CEILING( frac_day ) ) ) ) ! for now we let the year start on Monday.
521
522    !TBD: set weekday correct based on date
523    index_mm = mo
524    index_dd = nmonth + weekday  !> Index of the days in the mdh mode (13:20)
525
526    SELECT CASE(TRIM(daytype_mdh))
527
528       CASE ("workday")
529       
530          index_hh = nmonth+ nday + hh 
531
532       CASE ("weekend")
533       
534          index_hh = nmonth+ nday + nhour + hh 
535
536       CASE ("holiday")
537       
538          index_hh = nmonth+ nday + 2*nhour + hh 
539
540    END SELECT
541
542
543 END SUBROUTINE time_mdh_indices
544
545!------------------------------------------------------------------------------!
546! Description:
547! ------------
548!> This routine determines the time factor index in the HOURLY case of the DEFAULT emissions mode.
549!------------------------------------------------------------------------------!
550
551 SUBROUTINE time_hour_indices(mo,dd,hh,index_hh)
552
553    USE indices
554
555    IMPLICIT NONE
556
557    !> IN/OUTPUT
558    INTEGER, INTENT(INOUT)              :: mo          !> Month
559    INTEGER, INTENT(INOUT)              :: hh          !> Hour
560    INTEGER, INTENT(INOUT)              :: dd          !> Day
561    INTEGER, INTENT(INOUT)              :: index_hh    !> Index Hour
562 
563    !> Additional Variables for calculateing indices
564    INTEGER                             :: index_mm    !> Index Month
565    INTEGER                             :: index_dd    !> Index Day
566    INTEGER                             :: i_mon       !> Index for going through the different months
567    INTEGER                             :: sum_dd      !> Sum days
568
569    !> CONSTANTS
570    INTEGER, PARAMETER                  :: nhour = 24
571    INTEGER, PARAMETER, DIMENSION(12)   :: days = (/31,28,31,30,31,30,31,31,30,31,30,31/) ! no leap year
572   
573
574    index_mm = mo-1
575    index_dd = dd-1
576    sum_dd=0
577
578    IF (mo == 1) THEN
579
580       index_hh=(index_dd*nhour)+hh
581
582    ELSE
583
584       DO i_mon=1,index_mm
585
586         sum_dd=sum_dd+days(i_mon)
587
588       ENDDO
589     
590       index_hh=(sum_dd*nhour)+(index_dd*nhour)+(hh)
591
592    ENDIF
593 
594
595 END SUBROUTINE time_hour_indices
596
597
598END MODULE date_and_time_mod
Note: See TracBrowser for help on using the repository browser.