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

!> @file date_and_time_mod.f90
!------------------------------------------------------------------------------!
! This file is part of the PALM model system.
!
! PALM is free software: you can redistribute it and/or modify it under the
! terms of the GNU General Public License as published by the Free Software
! Foundation, either version 3 of the License, or (at your option) any later
! version.
!
! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
!
! You should have received a copy of the GNU General Public License along with
! PALM. If not, see <http://www.gnu.org/licenses/>.
!
! Copyright 1997-2019 Leibniz Universitaet Hannover
!------------------------------------------------------------------------------!
!
! Current revisions:
! ------------------
! 
! 
! Former revisions:
! -----------------
! $Id: date_and_time_mod.f90 3839 2019-03-28 21:12:25Z suehring $
! further tabs removed, unused variables removed
! 
! 3655 2019-01-07 16:51:22Z knoop
! further tabs removed
! 
! 3467 2018-10-30 19:05:21Z suehring
! Tabs removed
! 
! 3458 2018-10-30 14:51:23Z kanani
! from chemistry branch r3443, banzhafs:
! Added initial hour_of_day, hour_of_year, day_of_year and month_of_year to
! init_date_and_time
!
! 3298 2018-10-02 12:21:11Z kanani
! - 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
! 
! 2701 2017-12-15 15:40:50Z suehring
! Changes from last commit documented
!
! 2698 2017-12-14 18:46:24Z suehring
! Bugfix in definition of d_seconds_year.
!
! 2696 2017-12-14 17:12:51Z kanani
! Change in file header (GPL part)
! 
! 2544 2017-10-13 18:09:32Z maronga
! Initial revision
! 
! 

!
! Description:
! ------------
!> 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
!> @todo time_utc during spin-up  
!------------------------------------------------------------------------------!
 MODULE date_and_time_mod

    USE control_parameters,                                                   &
        ONLY:  coupling_start_time, days_since_reference_point,               &
               message_string, simulated_time, time_since_reference_point

    USE kinds


    IMPLICIT NONE

    PRIVATE

!-- 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)        ::  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 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 = "21062017"                  !< 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

    !-- Get hour index in the PRE-PROCESSED case of chemistry emissions :
    INTERFACE time_preprocessed_indices
       MODULE PROCEDURE time_preprocessed_indices
    END INTERFACE time_preprocessed_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, time_preprocessed_indices

    !-- 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

       !--    Variables Definition
       INTEGER ::  i_mon       !< Index for going through the different months

       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


       !-- Calculate initial 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_init/3600.0_wp ) ) + 1

       !-- Calculate initial day day_of_year_init in case date_init is given or day_of_year_init is given
       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

       
       !-- Initial day of the year
       day_of_year = day_of_year_init

       !-- Initial hour of the year
       hour_of_year = ( (day_of_year-1) * 24 ) + hour_of_day

       !--Initial 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

         
          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( 'init_date_and_time', 'DT0102', 2, 2, 0, 6, 0 ) 

       ENDIF


    END SUBROUTINE init_date_and_time

!------------------------------------------------------------------------------!
! Description:
! ------------
!> Calculate current date and time of the simulation 
!------------------------------------------------------------------------------!
 
    SUBROUTINE calc_date_and_time

       IMPLICIT NONE

!--    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 changes due to changes in time_since_reference_point in
!sB    radiation_model_mod during runtime:
       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

          ! fix the date if the day is 1st and earlier day is needed due to spinup
          IF ( day_of_month < 1 ) THEN
 
             ! if the day is the first day in the year
             IF ( month_of_year  ==  1 ) THEN
                month_of_year = 12
                day_of_month = 31
 
             ! other cases
             ELSE
                month_of_year = month_of_year - 1
                day_of_month = days(month_of_year)
             ENDIF
 
          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


!------------------------------------------------------------------------------!
! Description:
! ------------
!> This routine determines the time factor index in the PRE-PROCESSED emissions mode.
!------------------------------------------------------------------------------!

 SUBROUTINE time_preprocessed_indices(index_hh)

    USE indices

    IMPLICIT NONE

!
!-- In/output
    INTEGER, INTENT(INOUT) ::  index_hh    !> Index Hour
!
!-- Additional Variables for calculateing indices
!-- Constants
    INTEGER, PARAMETER ::  nhour = 24

    IF (days_since_reference_point == 0) THEN

       index_hh=hour_of_day

    ELSE

       index_hh=(days_since_reference_point*nhour)+(hour_of_day)

    ENDIF


 END SUBROUTINE time_preprocessed_indices


!------------------------------------------------------------------------------!
! 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                  :: 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

         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