Index: palm/trunk/SOURCE/date_and_time_mod.f90
===================================================================
--- palm/trunk/SOURCE/date_and_time_mod.f90	(revision 2718)
+++ palm/trunk/SOURCE/date_and_time_mod.f90	(revision 3298)
@@ -25,4 +25,12 @@
 ! -----------------
 ! $Id$
+! - Minor formatting (kanani)
+! - Added Routines for DEFAULT mode of chemistry emissions (Russo)
+! - Added routine for reading-in real dates: format has to be in DDMMYYYY and 
+!   passed in the namelist parameter date_init (Russo)
+! - Added calculation of several time indices useful in several routines 
+!   of the model (Russo)
+!
+! 2718 2018-01-02 08:49:38Z maronga
 ! Corrected "Former revisions" section
 ! 
@@ -46,38 +54,115 @@
 !> This routine calculates all needed information on date and time used by
 !> other modules
+!> @todo Further testing and revision of routines for updating indices of 
+!>       emissions in the default mode 
+!> @todo Add routine for recognizing leap years 
+!> @todo Add recognition of exact days of week (Monday, Tuesday, etc.)
+!> @todo Reconsider whether to remove day_of_year_init from the namelist: we 
+!>       already implemented changes for calculating it from date_init in 
+!>       calc_date_and_time
 !------------------------------------------------------------------------------!
  MODULE date_and_time_mod
- 
-    USE control_parameters,                                                    &
-        ONLY: time_since_reference_point 
- 
+
+    USE control_parameters,                                                   &
+        ONLY:  coupling_start_time, days_since_reference_point,               &
+               message_string, simulated_time, time_since_reference_point
+
     USE kinds
 
+
     IMPLICIT NONE
 
     PRIVATE
-    
-    PUBLIC   calc_date_and_time, d_hours_day, d_seconds_hour, d_seconds_year,  &
-             day_of_year, day_of_year_init, time_utc, time_utc_init
-
-
-    INTEGER(iwp) ::  day_of_year              !< day of the year (DOY)
-    INTEGER(iwp) ::  day_of_year_init = 172   !< DOY at model start (default: 21 June)
-
-    REAL(wp) ::  time_utc                     !< current model time in UTC
-    REAL(wp) ::  time_utc_init = 43200.0_wp   !< UTC time at model start
-
+
+!-- Variables Declaration
+
+    INTEGER(iwp)        ::  day_of_year      = 0   !< day of the year (DOY)
+    INTEGER(iwp)        ::  day_of_year_init = 0   !< DOY at model start (default: 0)
+
+    ! --- Most of these indices are updated by the routine calc_date_and_time according to the current date and time of the simulation
+    INTEGER(iwp)        ::  hour_of_year = 1                        !< hour of the current year (1:8760(8784))
+    INTEGER(iwp)        ::  hour_of_day=1                           !< hour of the current day (1:24)
+    INTEGER(iwp)        ::  day_of_month=0                          !< day of the current month (1:31)
+    INTEGER(iwp)        ::  month_of_year=0                         !< month of the current year (1:12)
+    INTEGER(iwp)        ::  current_year=0                          !< current year
+    INTEGER(iwp)        ::  hour_call_emis=0                        !< index used to call the emissions just once every hour
+
+    INTEGER(iwp)        ::  ihour            
+    INTEGER(iwp)        ::  index_mm                                !< index months of the default emission mode
+    INTEGER(iwp)        ::  index_dd                                !< index days of the default emission mode
+    INTEGER(iwp)        ::  index_hh                                !< index hours of the default emission mode
+
+    REAL(wp)            ::  time_utc                     !< current model time in UTC
+    REAL(wp)            ::  time_utc_emis                !< current emission module time in UTC
+    REAL(wp)            ::  time_utc_init = 43200.0_wp   !< UTC time at model start
+    REAL(wp)            ::  time_update                  !< used to calculate actual second of the simulation
+ 
     REAL(wp), PARAMETER ::  d_hours_day    = 1.0_wp / 24.0_wp       !< inverse of hours per day (1/24)
     REAL(wp), PARAMETER ::  d_seconds_hour = 1.0_wp / 3600.0_wp     !< inverse of seconds per hour (1/3600)
     REAL(wp), PARAMETER ::  d_seconds_year = 1.0_wp / 31536000.0_wp !< inverse of the seconds per year (1/(365*86400))
     
+    CHARACTER(len=8)    ::  date_init = "31122018"                  !< Starting date of simulation: We selected this because it was a monday
+ 
+    !> --- Parameters
+    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)
+
     SAVE
 
+!-- INTERFACES PART
+    !-- Read initial day and time of simulation
+    INTERFACE init_date_and_time
+       MODULE PROCEDURE init_date_and_time
+    END INTERFACE init_date_and_time
+
+    !-- Get hour index in the DEAFULT case of chemistry emissions :
+    INTERFACE time_default_indices
+       MODULE PROCEDURE time_mdh_indices
+       MODULE PROCEDURE time_hour_indices
+    END INTERFACE time_default_indices
+
+    !-- Calculate current date and time
+    INTERFACE calc_date_and_time
+       MODULE PROCEDURE calc_date_and_time
+    END INTERFACE
+
+
+    !-- Public Interfaces
+    PUBLIC calc_date_and_time, time_default_indices, init_date_and_time
+
+    !-- Public Variables
+    PUBLIC date_init, d_hours_day, d_seconds_hour, d_seconds_year,               &
+           day_of_year, day_of_year_init, time_utc, time_utc_init, day_of_month, &
+           month_of_year, index_mm, index_dd, index_hh, hour_of_day, hour_of_year, &
+           hour_call_emis
+
  CONTAINS
+
+
+ !------------------------------------------------------------------------------!
+ !> Reads starting date from namelist
+ !------------------------------------------------------------------------------!
+  
+    SUBROUTINE init_date_and_time
+
+       IMPLICIT NONE
+
+       IF  (day_of_year_init == 0) THEN
+          ! Day of the month at starting time
+          READ(UNIT=date_init(1:2),fmt=*)day_of_month
+
+          ! Month of the year at starting time
+          READ(UNIT=date_init(3:4),fmt=*)month_of_year
+
+          ! Year at starting time
+          READ(UNIT=date_init(5:8),fmt=*)current_year
+       
+       ENDIF
+
+    END SUBROUTINE init_date_and_time
 
 !------------------------------------------------------------------------------!
 ! Description:
 ! ------------
-!> Calculate current day of the year and time in UTC
+!> Calculate current date and time of the simulation 
 !------------------------------------------------------------------------------!
  
@@ -86,16 +171,255 @@
        IMPLICIT NONE
 
-!
-!--    Calculate current day of the year 
-       day_of_year = day_of_year_init + INT(FLOOR( (time_utc_init + time_since_reference_point)&
-                               / 86400.0_wp ), KIND=iwp)
-                        
+!--    Variables Definition
+       INTEGER		   		:: i_mon       !< Index for going through the different months
+
+       !> Update simulation time in seconds
+       time_update = simulated_time-coupling_start_time
+
+!--    Calculate current day of the simulated time
+       days_since_reference_point=INT(FLOOR( (time_utc_init + time_update) &
+                               / 86400.0_wp ) )
+
+!--    Calculate actual UTC time                        
        time_utc = MOD((time_utc_init + time_since_reference_point), 86400.0_wp)
        
+!sB    PRILIMINARY workaround for time_utc bug concerning time_since_reference_point:
+       time_utc_emis = MOD((time_utc_init + time_update), 86400.0_wp)     
+
+!--    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. 
+
+       IF ( day_of_year_init == 0 ) THEN
+
+          !> 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. 
+          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
+        
+             IF ( month_of_year == 1 ) THEN  !!month of year is read in input
+
+                day_of_year_init = day_of_month
+
+             ELSE
+
+                day_of_year_init= SUM(days( 1:(month_of_year-1) )) + day_of_month  !day_of_month is read in input in this case
+
+             ENDIF
+!kanani: Revise, we cannot force users to provide date_init, maybe set a default value?
+!           ELSE
+! 
+!              message_string = 'date_init not provided in the namelist or'             //          &
+!                               ' given in the wrong format: MUST BE DDMMYYYY'                        
+!              CALL message( 'calc_date_and_time', 'DT0100', 2, 2, 0, 6, 0 ) 
+     
+          ENDIF
+
+       ENDIF
+
+      !-- Calculate actual hour of the day: the first hour of the day is from 00:00:00 to 00:59:59. 
+
+      hour_of_day = INT( FLOOR( time_utc_emis/3600.0_wp ) ) + 1
+
+!--    Calculate current day of the year !TBD: considetr leap years
+       IF ( (day_of_year_init + days_since_reference_point)  .GT. 365 ) THEN
+
+          day_of_year=INT(MOD((day_of_year_init + days_since_reference_point), 365.0_wp))
+
+       ELSE
+          
+          day_of_year = day_of_year_init + days_since_reference_point
+
+       ENDIF
+
+!
+!--    Calculate current hour of the year
+       hour_of_year = ( (day_of_year-1) * 24 ) + hour_of_day  !> actual hour of the year
+        
+
+!
+!--    UPDATE actual day of the month and month of the year
+       !> --------------------------------------------------------------------------------
+       !> The first case is when date_init is not provided: we only know day_of_year_init      
+       IF ( month_of_year == 0 .AND. day_of_month == 0) THEN
+
+          !> The first case is when date_init is not provided: we only know day_of_year_init
+          !DO i_mon=1,12
+             !IF (day_of_year .LE. SUM(days(1:i_mon))) THEN
+          IF ( day_of_year .LE. 31 ) THEN
+
+             month_of_year=1
+             day_of_month=day_of_year
+
+          ELSE
+
+             DO i_mon=2,12   !january is considered in the first case
+                IF ( day_of_year .LE. SUM(days(1:i_mon)) .AND. day_of_year .GT. SUM(days(1:(i_mon-1))) ) THEN
+           
+                   month_of_year=i_mon
+
+                   day_of_month=INT(MOD(day_of_year, SUM(days(1:(i_mon-1))))) 
+
+                   GOTO 38
+
+                ENDIF
+
+             38 ENDDO
+          ENDIF
+       !> --------------------------------------------------------------------------------
+       !> 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
+       !>TBD: something to calculate the current year is missing
+       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
+ 
+          !> calculate month_of_year. TBD: test the condition when day_of_year==31
+ 
+          IF (day_of_year==1) THEN  !> this allows to turn from december to January when passing from a year to another
+ 
+             month_of_year = 1
+        
+          ELSE IF (day_of_year .GT. 1 .AND. day_of_year .GT. SUM(days(1:month_of_year))) THEN
+
+             month_of_year = month_of_year + 1
+
+          ENDIF
+
+          !> calculate day_of_month
+          IF ( month_of_year == 1 ) THEN
+           
+            day_of_month=day_of_year
+
+          ELSE
+
+            day_of_month=INT(MOD(day_of_year, SUM(days(1:(month_of_year-1)))))
+
+          ENDIF
+
+
+
+
+      ELSE
+
+          !> Condition when date_init is provided but it is given in the wrong format
+          message_string = 'date_init not provided in the namelist or'            //          &
+                              ' given in the wrong format: MUST BE DDMMYYYY'                  
+          CALL message( 'calc_date_and_time', 'DT0101', 2, 2, 0, 6, 0 ) 
+
+      ENDIF       
+      
+    END SUBROUTINE calc_date_and_time
+
+!TBD: check the routines used for update of emission indices in the DEFAULT MODE 
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> This routine determines the time factor index in the mdh case of the DEFAULT chemistry emissions mode.
+!------------------------------------------------------------------------------!
+
+ SUBROUTINE time_mdh_indices(daytype_mdh,mo, dd, hh, index_mm, index_dd, index_hh)
+
+    USE indices
+
+    IMPLICIT NONE
+
+    !> IN/OUTPUT
+    INTEGER, INTENT(INOUT) :: mo          !> Month of year
+    INTEGER, INTENT(INOUT) :: dd          !> Day of month
+    INTEGER, INTENT(INOUT) :: hh          !> Hour of day
+    INTEGER, INTENT(INOUT) :: index_mm    !> Index Month
+    INTEGER, INTENT(INOUT) :: index_dd    !> Index Day
+    INTEGER, INTENT(INOUT) :: index_hh    !> Index Hour
+
+    CHARACTER(len=80), INTENT(INOUT) :: daytype_mdh !> type of the day in mdh mode: one of 1-WORKDAY
+                                                    !                                      2-WEEKEND
+                                                    !                                      3-HOLIDAY
+
+    REAL(wp)                         :: frac_day=0  
+
+    !> ------------------------------------------------------------------------
+
+    INTEGER                	     :: weekday
+
+    !> CONSTANTS
+    INTEGER, PARAMETER               :: nmonth = 12
+    INTEGER, PARAMETER               :: nday = 7
+    INTEGER, PARAMETER               :: nhour = 24
+
+    frac_day= (dd-1)/nday    !> indicates the week of the month, supposing the month starts on monday
+
+   ! 1:7      1:31  7                   (0:30)/7
+    weekday = dd-( nday * (INT( CEILING( frac_day ) ) ) ) ! for now we let the year start on Monday. 
+
+    !TBD: set weekday correct based on date
+    index_mm = mo
+    index_dd = nmonth + weekday  !> Index of the days in the mdh mode (13:20)
+
+    SELECT CASE(TRIM(daytype_mdh))
+
+		   CASE ("workday")
+	       
+	           index_hh = nmonth+ nday + hh 
+
+		   CASE ("weekend")
+	       
+	           index_hh = nmonth+ nday + nhour + hh 
+
+		   CASE ("holiday")
+	       
+	           index_hh = nmonth+ nday + 2*nhour + hh 
+    END SELECT
+
+
+ END SUBROUTINE time_mdh_indices
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> This routine determines the time factor index in the HOURLY case of the DEFAULT emissions mode.
+!------------------------------------------------------------------------------!
+
+ SUBROUTINE time_hour_indices(mo,dd,hh,index_hh)
+
+    USE indices
+
+    IMPLICIT NONE
+
+    !> IN/OUTPUT
+    INTEGER, INTENT(INOUT) 		:: mo          !> Month
+    INTEGER, INTENT(INOUT) 		:: hh          !> Hour
+    INTEGER, INTENT(INOUT) 		:: dd          !> Day
+    INTEGER, INTENT(INOUT) 		:: index_hh    !> Index Hour
+ 
+    !> Additional Variables for calculateing indices
+    INTEGER		   		:: index_mm    !> Index Month
+    INTEGER		   		:: index_dd    !> Index Day
+    INTEGER		   		:: i_mon       !> Index for going through the different months
+    INTEGER		   		:: sum_dd      !> Sum days
+
+    !> CONSTANTS
+    INTEGER, PARAMETER 			:: nmonth=12
+    INTEGER, PARAMETER 			:: nday = 7
+    INTEGER, PARAMETER 			:: nhour = 24
+    INTEGER, PARAMETER, DIMENSION(12)   :: days = (/31,28,31,30,31,30,31,31,30,31,30,31/) ! no leap year
+    
+
+    index_mm = mo-1
+    index_dd = dd-1
+    sum_dd=0
+
+    IF (mo == 1) THEN
+
+       index_hh=(index_dd*nhour)+hh
+
+    ELSE 
+
+    	DO i_mon=1,index_mm
        
-
-    END SUBROUTINE calc_date_and_time
-
-
-
- END MODULE date_and_time_mod
+            sum_dd=sum_dd+days(i_mon)
+
+        ENDDO
+      
+    index_hh=(sum_dd*nhour)+(index_dd*nhour)+(hh)
+
+    ENDIF
+ 
+
+ END SUBROUTINE time_hour_indices
+
+
+END MODULE date_and_time_mod
