Line data    Source code

       1             : !> \file    mo_edk_get_nc_time.f90
       2             : !> \brief   \copybrief mo_edk_get_nc_time
       3             : !> \details \copydetails mo_edk_get_nc_time
       4             : 
       5             : !> \brief   Module to get time vector from NetCDF file.
       6             : !> \author  Matthias Zink
       7             : !> \date    Oct 2012
       8             : !> \authors Matthias Cuntz, Juliane Mai
       9             : !> \date    Nov 2014
      10             : !!          - time int or double
      11             : !> \author  Stephan Thober
      12             : !> \date    Sep 2015
      13             : !!          - added read for hourly data
      14             : !> \author  Robert Schweppe
      15             : !> \date    Nov 2017
      16             : !!          - restructured routine, reads vector now
      17             : !> \author  Maren Kaluza
      18             : !> \date    May 2018
      19             : !!          - fixed bug in time reading
      20             : !> \copyright Copyright 2005-\today, the CHS Developers, Sabine Attinger: All rights reserved.
      21             : !! EDK is released under the LGPLv3+ license \license_note
      22             : module mo_edk_get_nc_time
      23             : 
      24             :   implicit none
      25             : 
      26             :   private
      27             : 
      28             :   public :: get_time_vector_and_select
      29             : 
      30             : contains
      31             : 
      32             :   !> \brief   Extract time vector
      33             :   !> \details Extract time vector in unit julian hours and get supposed time step in hours
      34           0 :   subroutine get_time_vector_and_select(var, fname, inctimestep, time_start, time_cnt, target_period)
      35             : 
      36             :     use mainVar, only : period, dayhours, daysecs, yeardays
      37             :     use mo_julian, only : caldat, dec2date, julday
      38             :     use mo_kind, only : dp, i4, i8
      39             :     use mo_message, only : message
      40             :     use mo_netcdf, only : NcVariable
      41             :     use mo_string_utils, only : DIVIDE_STRING
      42             : 
      43             :     implicit none
      44             : 
      45             :     !> variable of interest
      46             :     type(NcVariable), intent(in) :: var
      47             : 
      48             :     !> fname of ncfile for error message
      49             :     character(256), intent(in) :: fname
      50             : 
      51             :     !> flag for requested time step
      52             :     integer(i4), intent(in) :: inctimestep
      53             : 
      54             :     !> time_start index of time selection
      55             :     integer(i4), intent(out) :: time_start
      56             : 
      57             :     !> time_count of indexes of time selection
      58             :     integer(i4), intent(out) :: time_cnt
      59             : 
      60             :     !> reference period
      61             :     type(period), intent(in), optional :: target_period
      62             : 
      63             :     ! reference time of NetCDF
      64             :     integer(i4) :: yRef, dRef, mRef, hRef, jRef
      65             : 
      66             :     ! netcdf attribute values
      67             :     character(256) :: AttValues
      68             : 
      69             :     ! dummies for netcdf attribute handling
      70           0 :     character(256), dimension(:), allocatable :: strArr, date, time
      71             : 
      72             :     ! native time step converter in ncfile
      73             :     integer(i8) :: time_step_seconds
      74             : 
      75             :     ! time vector
      76           0 :     integer(i8), allocatable, dimension(:) :: time_data
      77             : 
      78             :     ! period of ncfile, for clipping
      79             :     type(period) :: nc_period, clip_period
      80             : 
      81             :     integer(i4) :: ncJulSta1, dd, n_time
      82             : 
      83             :     integer(i4) :: mmcalstart, mmcalend, yycalstart, yycalend
      84             : 
      85             :     integer(i4) :: mmncstart, yyncstart
      86             : 
      87             :     ! helper variable for error output
      88             :     integer(i4) :: hstart_int, hend_int
      89             : 
      90             :     ! helper variable for error output
      91             :     character(256) :: error_msg
      92             : 
      93             : 
      94           0 :     call var%getAttribute('units', AttValues)
      95             :     ! AttValues looks like "<unit> since YYYY-MM-DD[ HH:MM:SS]"
      96             :     ! split at space
      97           0 :     call DIVIDE_STRING(trim(AttValues), ' ', strArr)
      98             : 
      99             :     ! determine reference time at '-' and convert to integer
     100           0 :     call DIVIDE_STRING(trim(strArr(3)), '-', date)
     101           0 :     read(date(1), *) yRef
     102           0 :     read(date(2), *) mRef
     103           0 :     read(date(3), *) dRef
     104             : 
     105           0 :     jRef = julday(dd = dRef, mm = mRef, yy = yRef)
     106             : 
     107             :     ! if existing also read in the time (only hour so far)
     108           0 :     hRef = 0
     109           0 :     if(size(strArr) .gt. 3) then
     110           0 :       call DIVIDE_STRING(trim(strArr(4)), ':', time)
     111           0 :       read(time(1), *) hRef
     112             :     end if
     113             : 
     114             :     ! determine the step_size
     115           0 :     if (strArr(1) .EQ. 'days') then
     116             :       time_step_seconds = int(DaySecs)
     117           0 :     else if (strArr(1) .eq. 'hours') then
     118             :       time_step_seconds = int(DaySecs / DayHours)
     119           0 :     else if (strArr(1) .eq. 'minutes') then
     120             :       time_step_seconds = int(DaySecs / DayHours / 60._dp)
     121           0 :     else if (strArr(1) .eq. 'seconds') then
     122             :       time_step_seconds = 1_i8
     123             :     else
     124             :       call message('***ERROR: Please provide the input data in (days, hours, minutes, seconds) ', &
     125           0 :               'since YYYY-MM-DD[ HH:MM:SS] in the netcdf file. Found: ', trim(AttValues))
     126           0 :       stop
     127             :     end if
     128             : 
     129             :     ! get the time vector
     130           0 :     call var%getData(time_data)
     131             :     ! convert array from units since to seconds
     132           0 :     time_data = time_data * time_step_seconds
     133             : 
     134             :     ! check for length of time vector, needs to be at least of length 2, otherwise step width check fails
     135           0 :     if (size(time_data) .le. 1) then
     136           0 :       call message('***ERROR: length of time dimension needs to be at least 2 in file: ' // trim(fname))
     137           0 :       stop
     138             :     end if
     139             : 
     140             :     ! check for equal timesteps and timestep must not be multiple of native timestep
     141             :     error_msg = '***ERROR: time_steps are not equal over all times in file and/or do not conform to' // &
     142           0 :             ' requested timestep in file (' // trim(fname) // ') : '
     143             : 
     144             :     ! compare the read period from ncfile to the period required
     145             :     ! convert julian second information back to date via conversion to float
     146             :     ! the 0.5_dp is for the different reference of fractional julian days, hours are truncated
     147           0 :     n_time = size(time_data)
     148           0 :     call dec2date(time_data(1) / DaySecs - 0.5_dp + jRef + hRef / 24._dp, nc_period%dStart, nc_period%mStart, &
     149           0 :             nc_period%yStart, hstart_int)
     150           0 :     nc_period%julStart = int(time_data(1) / DaySecs + jRef + hRef / 24._dp)
     151           0 :     call dec2date(time_data(n_time) / DaySecs - 0.5_dp + jRef + hRef / 24._dp, nc_period%dEnd, nc_period%mEnd, &
     152           0 :             nc_period%yEnd, hend_int)
     153           0 :     nc_period%julEnd = int(time_data(n_time) / DaySecs + jRef + hRef / 24._dp)
     154             : 
     155             :     ! if no target period is present, use the whole time period
     156           0 :     if (present(target_period)) then
     157           0 :       clip_period = target_period
     158             :     else
     159           0 :       clip_period = nc_period
     160             :     end if
     161             : 
     162             :     ! prepare the selection and check for required time_step
     163           0 :     select case(inctimestep)
     164             :     case(-1) ! daily
     165             :       ! difference must be 1 day
     166           0 :       if (.not. all(abs((time_data(2 : n_time) - time_data(1 : n_time - 1)) / DaySecs - 1._dp) .lt. 1.e-6)) then
     167           0 :         call message(error_msg // trim('daily'))
     168           0 :         stop
     169             :       end if
     170           0 :       ncJulSta1 = nc_period%julStart
     171           0 :       time_start = clip_period%julStart - ncJulSta1 + 1_i4
     172           0 :       time_cnt = clip_period%julEnd - clip_period%julStart + 1_i4
     173             :     case(-2) ! monthly
     174             :       ! difference must be between 28 and 31 days
     175           0 :       if (any(abs((time_data(2 : n_time) - time_data(1 : n_time - 1)) / DaySecs) .gt. 31._dp) .or. &
     176             :               any(abs((time_data(2 : n_time) - time_data(1 : n_time - 1)) / DaySecs) .lt. 28._dp)) then
     177           0 :         call message(error_msg // trim('monthly'))
     178           0 :         stop
     179             :       end if
     180             : 
     181           0 :       call caldat(clip_period%julStart, dd, mmcalstart, yycalstart)
     182           0 :       call caldat(nc_period%julStart, dd, mmncstart, yyncstart)
     183             :       ! monthly timesteps are usually set by month end, so for beginning, we need 1st of month
     184           0 :       ncJulSta1 = julday(1, mmncstart, yyncstart)
     185           0 :       call caldat(clip_period%julEnd, dd, mmcalend, yycalend)
     186           0 :       time_start = (yycalstart * 12 + mmcalstart) - (yyncstart * 12 + mmncstart) + 1_i4
     187           0 :       time_cnt = (yycalend * 12 + mmcalend) - (yycalstart * 12 + mmcalstart) + 1_i4
     188             :     case(-3) ! yearly
     189             :       ! difference must be between 365 and 366 days
     190           0 :       if (any(abs((time_data(2 : n_time) - time_data(1 : n_time - 1)) / DaySecs) .gt. (YearDays + 1._dp)) .or. &
     191             :               any(abs((time_data(2 : n_time) - time_data(1 : n_time - 1)) / DaySecs) .lt. YearDays)) then
     192           0 :         call message(error_msg // 'yearly')
     193           0 :         stop
     194             :       end if
     195           0 :       call caldat(clip_period%julStart, dd, mmcalstart, yycalstart)
     196           0 :       call caldat(nc_period%julStart, dd, mmncstart, yyncstart)
     197             :       ! yearly timesteps are usually set by year end, so for beginning, we need 1st of year
     198           0 :       ncJulSta1 = julday(1, 1, yyncstart)
     199           0 :       call caldat(clip_period%julEnd, dd, mmcalend, yycalend)
     200           0 :       time_start = yycalstart - yyncstart + 1_i4
     201           0 :       time_cnt = yycalend - yycalstart + 1_i4
     202             :     case(-4) ! hourly
     203             :       ! difference must be 1 hour
     204           0 :       if (.not. all(abs((time_data(2 : n_time) - time_data(1 : n_time - 1)) / 3600._dp - 1._dp) .lt. 1.e-6)) then
     205           0 :         call message(error_msg // 'hourly')
     206           0 :         stop
     207             :       end if
     208           0 :       ncJulSta1 = nc_period%julStart
     209           0 :       time_start = (clip_period%julStart - ncJulSta1) * 24_i4 + 1_i4 ! convert to hours; always starts at one
     210           0 :       time_cnt = (clip_period%julEnd - clip_period%julStart + 1_i4) * 24_i4 ! convert to hours
     211             :     case default ! no output at all
     212           0 :       call message('***ERROR: read_forcing_nc: unknown nctimestep switch.')
     213           0 :       stop
     214             :     end select
     215             : 
     216             :     ! Check if time steps in file cover simulation period
     217           0 :     if (.not. ((ncJulSta1 .LE. clip_period%julStart) .AND. (nc_period%julEnd .GE. clip_period%julEnd))) then
     218             :       call message('***ERROR: read_forcing_nc: time period of input data: ', trim(fname), &
     219           0 :               '          is not matching modelling period.')
     220           0 :       stop
     221             :     end if
     222             : 
     223           0 :   end subroutine get_time_vector_and_select
     224             : 
     225             : end module mo_edk_get_nc_time