LCOV - code coverage report
Current view: | top level - hemco/HEMCO/src/Core - hcoio_util_mod.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 0 | 1024 | 0.0 %
Date: | 2024-12-17 17:57:11 |  | Functions: | 0 | 21 | 0.0 %

              Line data    Source code
           1             : !------------------------------------------------------------------------------
           2             : !                   Harmonized Emissions Component (HEMCO)                    !
           3             : !------------------------------------------------------------------------------
           4             : !BOP
           5             : !
           6             : ! !MODULE: hcoio_util_mod.F90
           7             : !
           8             : ! !DESCRIPTION: Module HCOIO\_Util\_Mod contains utility functions
           9             : ! for use in data processing including file reading, unit conversions,
          10             : ! and regridding.
          11             : !\\
          12             : !\\
          13             : ! !INTERFACE:
          14             : !
          15             : MODULE HCOIO_Util_Mod
          16             : !
          17             : ! !USES:
          18             : !
          19             :   USE HCO_Types_Mod
          20             :   USE HCO_Error_Mod
          21             :   USE HCO_CharTools_Mod
          22             :   USE HCO_State_Mod,       ONLY : Hco_State
          23             :
          24             :   IMPLICIT NONE
          25             :   PRIVATE
          26             : !
          27             : ! !PUBLIC MEMBER FUNCTIONS:
          28             : !
          29             : #if !defined(ESMF_)
          30             :   PUBLIC :: GET_TIMEIDX
          31             :   PUBLIC :: Check_AvailYMDhm
          32             :   PUBLIC :: prefYMDhm_Adjust
          33             :   PUBLIC :: Set_tIdx2
          34             :   PUBLIC :: IsClosest
          35             :   PUBLIC :: GetIndex2Interp
          36             :   PUBLIC :: GetWeights
          37             :   PUBLIC :: YMDhm2hrs
          38             :   PUBLIC :: Normalize_Area
          39             :   PUBLIC :: SrcFile_Parse
          40             :   PUBLIC :: SigmaMidToEdges
          41             :   PUBLIC :: CheckMissVal
          42             :   PUBLIC :: GetArbDimIndex
          43             : #endif
          44             :   PUBLIC :: HCOIO_ReadOther
          45             :   PUBLIC :: HCOIO_ReadCountryValues
          46             :   PUBLIC :: HCOIO_ReadFromConfig
          47             :   PUBLIC :: GetDataVals
          48             :   PUBLIC :: GetSliceIdx
          49             :   PUBLIC :: FillMaskBox
          50             :   PUBLIC :: ReadMath
          51             : !
          52             : ! !REVISION HISTORY:
          53             : !  12 Jun 2020 - E. Lundgren - Initial version, created from subset of
          54             : !                              hcoio_util_mod.F90
          55             : !  See https://github.com/geoschem/hemco for complete history
          56             : !EOP
          57             : !------------------------------------------------------------------------------
          58             : !BOC
          59             : !
          60             : ! !DEFINED PARAMETERS
          61             : !
          62             :   ! Parameter used for difference testing of floating points
          63             :   REAL(dp), PRIVATE, PARAMETER :: EPSILON = 1.0e-5_dp
          64             :
          65             : CONTAINS
          66             : !EOC
          67             : #if !defined( ESMF_ )
          68             : !------------------------------------------------------------------------------
          69             : !                   Harmonized Emissions Component (HEMCO)                    !
          70             : !------------------------------------------------------------------------------
          71             : !BOP
          72             : !
          73             : ! !IROUTINE: Get_TimeIdx
          74             : !
          75             : ! !DESCRIPTION: Returns the lower and upper time slice index (tidx1
          76             : ! and tidx2, respectively) to be read. These values are determined
          77             : ! based upon the time slice information extracted from the netCDF file,
          78             : ! the time stamp settings set in the config. file, and the current
          79             : ! simulation date.
          80             : !\\
          81             : !\\
          82             : ! Return arguments wgt1 and wgt2 denote the weights to be given to
          83             : ! the two time slices. This is only of relevance for data that shall
          84             : ! be interpolated between two (not necessarily consecutive) time slices.
          85             : ! In all other cases, the returned weights are negative and will be
          86             : ! ignored.
          87             : !\\
          88             : !\\
          89             : ! Also returns the time slice year and month, as these values may be
          90             : ! used for unit conversion.
          91             : !\\
          92             : !\\
          93             : ! !INTERFACE:
          94             : !
          95           0 :   SUBROUTINE GET_TIMEIDX( HcoState,  Lct,               &
          96             :                           ncLun,     tidx1,    tidx2,   &
          97             :                           wgt1,      wgt2,     oYMDhm,  &
          98             :                           YMDhm,     YMDhm1,   RC,      &
          99             :                           Year )
         100             : !
         101             : ! !USES:
         102             : !
         103             :     USE HCO_Ncdf_Mod,  ONLY : NC_Read_Time_YYYYMMDDhhmm
         104             :     USE HCO_tIdx_Mod,  ONLY : HCO_GetPrefTimeAttr
         105             : !
         106             : ! !INPUT PARAMETERS:
         107             : !
         108             :     TYPE(HCO_State),  POINTER                  :: HcoState  ! HcoState object
         109             :     TYPE(ListCont),   POINTER                  :: Lct       ! List container
         110             :     INTEGER,          INTENT(IN   )            :: ncLun     ! open ncLun
         111             :     INTEGER,          INTENT(IN   ), OPTIONAL  :: Year      ! year to be used
         112             : !
         113             : ! !OUTPUT PARAMETERS:
         114             : !
         115             :     INTEGER,          INTENT(  OUT)            :: tidx1  ! lower time idx
         116             :     INTEGER,          INTENT(  OUT)            :: tidx2  ! upper time idx
         117             :     REAL(sp),         INTENT(  OUT)            :: wgt1   ! weight to tidx1
         118             :     REAL(sp),         INTENT(  OUT)            :: wgt2   ! weight to tidx2
         119             :     REAL(dp),         INTENT(  OUT)            :: oYMDhm ! preferred time slice
         120             :     REAL(dp),         INTENT(  OUT)            :: YMDhm  ! selected time slice
         121             :     REAL(dp),         INTENT(  OUT)            :: YMDhm1 ! 1st time slice in file
         122             : !
         123             : ! !INPUT/OUTPUT PARAMETERS:
         124             : !
         125             :     INTEGER,          INTENT(INOUT)            :: RC
         126             : !
         127             : ! !REVISION HISTORY:
         128             : !  13 Mar 2013 - C. Keller - Initial version
         129             : !  See https://github.com/geoschem/hemco for complete history
         130             : !EOP
         131             : !------------------------------------------------------------------------------
         132             : !BOC
         133             : !
         134             : ! !LOcAL VARIABLES:
         135             : !
         136             :     CHARACTER(LEN=255)    :: MSG, LOC
         137             :     CHARACTER(LEN=1023)   :: MSG_LONG
         138             :     INTEGER               :: tidx1a
         139             :     INTEGER               :: nTime,  T, CNT, NCRC
         140             :     INTEGER               :: prefYr, prefMt, prefDy, prefHr, prefMn
         141             :     INTEGER               :: refYear
         142             :     REAL(dp)              :: origYMDhm, prefYMDhm
         143           0 :     REAL(dp),   POINTER   :: availYMDhm(:)
         144             :     LOGICAL               :: ExitSearch
         145             :     LOGICAL               :: verb
         146             :
         147             :     !=================================================================
         148             :     ! GET_TIMEIDX begins here
         149             :     !=================================================================
         150             :
         151             :     ! Initialize
         152           0 :     LOC         = 'GET_TIMEIDX (HCOIO_UTIL_MOD.F90)'
         153             :
         154             :     ! Officially enter Get_TimeIdx
         155           0 :     CALL HCO_ENTER( HcoState%Config%Err, LOC, RC )
         156           0 :     IF ( RC /= HCO_SUCCESS ) THEN
         157           0 :         CALL HCO_ERROR( 'ERROR 0', RC, THISLOC=LOC )
         158           0 :         RETURN
         159             :     ENDIF
         160           0 :     verb = HCO_IsVerb( HcoState%Config%Err )
         161             :
         162             :     ! Initialize local variables for safety's sake
         163           0 :     nTime      =  0
         164           0 :     cnt        =  0
         165           0 :     prefYr     =  0
         166           0 :     prefMt     =  0
         167           0 :     prefDy     =  0
         168           0 :     prefHr     =  0
         169           0 :     prefMn     =  0
         170           0 :     refYear    =  0
         171           0 :     origYMDhm  =  0
         172           0 :     prefYMDhm  =  0
         173           0 :     tidx1      =  0
         174           0 :     tidx2      =  0
         175           0 :     tidx1a     =  0
         176           0 :     wgt1       = -1.0_sp
         177           0 :     wgt2       = -1.0_sp
         178           0 :     oYMDhm     =  0.0_dp
         179           0 :     YMDhm      =  0.0_dp
         180           0 :     YMDhm1     =  0.0_dp
         181           0 :     ExitSearch = .FALSE.
         182           0 :     availYMDhm => NULL()
         183             :
         184             :     ! ----------------------------------------------------------------
         185             :     ! Extract netCDF time slices (YYYYMMDDhhmm)
         186             :     ! ----------------------------------------------------------------
         187             :     CALL NC_READ_TIME_YYYYMMDDhhmm( ncLun, nTime,    availYMDhm,  &
         188           0 :                                     refYear=refYear, RC=NCRC     )
         189           0 :     IF ( NCRC /= 0 ) THEN
         190           0 :        CALL HCO_ERROR( 'NC_READ_TIME_YYYYMMDDhhmm', RC )
         191           0 :        RETURN
         192             :     ENDIF
         193             :
         194             :     ! Return warning if netCDF reference year prior to 1901: it seems
         195             :     ! like there are some problems with that and the time slices can be
         196             :     ! off by one day!
         197           0 :     IF ( (refYear <= 1900) .AND. (nTime > 0) ) THEN
         198             :        MSG = 'ncdf reference year is prior to 1901 - ' // &
         199           0 :             'time stamps may be wrong!'
         200           0 :        CALL HCO_WARNING ( HcoState%Config%Err, MSG, RC )
         201             :     ENDIF
         202             :
         203             :     ! verbose mode
         204           0 :     IF ( verb ) THEN
         205           0 :        write(MSG,*) 'Number of time slices found: ', nTime
         206           0 :        CALL HCO_MSG(HcoState%Config%Err,MSG)
         207           0 :        IF ( nTime > 0 ) THEN
         208           0 :           write(MSG,*) 'Time slice range : ', &
         209           0 :                        availYMDhm(1), availYMDhm(nTime)
         210           0 :           CALL HCO_MSG(HcoState%Config%Err,MSG)
         211             :        ENDIF
         212             :     ENDIF
         213             :
         214             :     ! ----------------------------------------------------------------
         215             :     ! Select time slices to read
         216             :     ! ----------------------------------------------------------------
         217             :
         218             :     ! ----------------------------------------------------------------
         219             :     ! Get preferred time stamp to read based upon the specs set in the
         220             :     ! config. file.
         221             :     ! This can return value -1 for prefHr, indicating that all
         222             :     ! corresponding time slices shall be read.
         223             :     ! This call will return -1 for all date attributes if the
         224             :     ! simulation date is outside of the data range given in the
         225             :     ! configuration file.
         226             :     ! ----------------------------------------------------------------
         227             :     CALL HCO_GetPrefTimeAttr ( HcoState, Lct, &
         228           0 :                                prefYr, prefMt, prefDy, prefHr, prefMn, RC )
         229           0 :     IF ( RC /= HCO_SUCCESS ) THEN
         230             :        MSG = &
         231           0 :          'Error encountered in HCO_GetPrefTimeAttr for ' // TRIM(Lct%Dct%cName)
         232           0 :        CALL HCO_ERROR( MSG, RC )
         233           0 :        IF ( ASSOCIATED(availYMDhm) ) THEN
         234           0 :           DEALLOCATE(availYMDhm)
         235             :           availYMDhm => NULL()
         236             :        ENDIF
         237           0 :        RETURN
         238             :     ENDIF
         239             :
         240             :     ! Eventually force preferred year to passed value
         241           0 :     IF ( PRESENT(Year) ) prefYr = Year
         242             :
         243             :     ! Check if we are outside of provided range
         244           0 :     IF ( prefYr < 0 .OR. prefMt < 0 .OR. prefDy < 0 ) THEN
         245             :
         246             :        ! This should only happen for 'range' data
         247           0 :        IF ( Lct%Dct%Dta%CycleFlag /= HCO_CFLAG_RANGE ) THEN
         248           0 :           MSG = 'Cannot get preferred datetime for ' // TRIM(Lct%Dct%cName)
         249           0 :           CALL HCO_ERROR( MSG, RC )
         250           0 :           IF ( ASSOCIATED(availYMDhm) ) THEN
         251           0 :              DEALLOCATE(availYMDhm)
         252             :              availYMDhm => NULL()
         253             :           ENDIF
         254           0 :           RETURN
         255             :        ENDIF
         256             :
         257             :        ! If this part of the code gets executed, the data associated
         258             :        ! with this container shall not be used at the current date.
         259             :        ! To do so, set the time indeces to -1 and leave right here.
         260           0 :        tidx1 = -1
         261           0 :        tidx2 = -1
         262             :
         263             :        ! Leave w/ success
         264           0 :        CALL HCO_LEAVE( HcoState%Config%Err,  RC )
         265           0 :        RETURN
         266             :     ENDIF
         267             :
         268             :     ! origYMDhm is the preferred datetime. Store into shadow variable
         269             :     ! prefYMDhm. prefYMDhm may be adjusted if origYMDhm is outside of the
         270             :     ! netCDF datetime range.
         271             :     ! Now put origYMDhm, prefYMDhm in YYYYMMDDhhmm format (bmy, 4/10/17)
         272             :     origYMDhm = ( DBLE(      prefYr      ) * 1.0e8_dp ) + &
         273             :                 ( DBLE(      prefMt      ) * 1.0e6_dp ) + &
         274             :                 ( DBLE(      prefDy      ) * 1.0e4_dp ) + &
         275             :                 ( DBLE( MAX( prefHr, 0 ) ) * 1.0e2_dp ) + &
         276           0 :                 ( DBLE( MAX( prefMn, 0 ) )          )
         277           0 :     prefYMDhm = origYMDhm
         278             :
         279             :     ! verbose mode
         280           0 :     IF ( verb ) THEN
         281           0 :        write(MSG,'(A30,f14.0)') 'preferred datetime: ', prefYMDhm
         282           0 :        CALL HCO_MSG(HcoState%Config%Err,MSG)
         283             :     ENDIF
         284             :
         285             :     ! ================================================================
         286             :     ! Case 1: Only one time slice available.
         287             :     ! ================================================================
         288           0 :     IF ( nTime == 1 ) THEN
         289           0 :        tidx1 = 1
         290           0 :        tidx2 = 1
         291             :
         292             :     ! ================================================================
         293             :     ! Case 2: More than one time slice available. Determine lower
         294             :     ! and upper time slice index from file & HEMCO settings.
         295             :     ! ================================================================
         296           0 :     ELSEIF ( nTime > 1 ) THEN
         297             :
         298             :        ! Init
         299           0 :        tidx1   = -1
         300           0 :        tidx2   = -1
         301             :
         302             :        ! -------------------------------------------------------------
         303             :        ! Check if preferred datetime prefYMDhm is within the range
         304             :        ! available time slices, e.g. it falls within the interval
         305             :        ! of availYMDhm. In this case, set tidx1 to the index of the
         306             :        ! closest time slice that is not in the future.
         307             :        ! -------------------------------------------------------------
         308           0 :        CALL Check_AvailYMDhm ( Lct, nTime, availYMDhm, prefYMDhm, tidx1a )
         309             :
         310             :        ! -------------------------------------------------------------
         311             :        ! Check if we need to continue search. Even if the call above
         312             :        ! returned a time slice, it may be possible to continue looking
         313             :        ! for a better suited time stamp. This is only the case if
         314             :        ! there are discontinuities in the time stamps, e.g. if a file
         315             :        ! contains monthly data for 2005 and 2020. In that case, the
         316             :        ! call above would return the index for Dec 2005 for any
         317             :        ! simulation date between 2005 and 2010 (e.g. July 2010),
         318             :        ! whereas it makes more sense to use July 2005 (and eventually
         319             :        ! interpolate between the July 2005 and July 2020 data).
         320             :        ! The IsClosest command checks if there are any netCDF time
         321             :        ! stamps (prior to the selected one) that are closer to each
         322             :        ! other than the difference between the preferred time stamp
         323             :        ! prefYMDhm and the currently selected time stamp
         324             :        ! availYMDhm(tidx1a). In that case, it continues the search by
         325             :        ! updating prefYMDhm so that it falls within the range of the
         326             :        ! 'high-frequency' interval.
         327             :        ! -------------------------------------------------------------
         328           0 :        ExitSearch = .FALSE.
         329           0 :        IF ( Lct%Dct%Dta%CycleFlag == HCO_CFLAG_EXACT ) THEN
         330             :           ExitSearch = .TRUE.
         331           0 :        ELSE IF ( tidx1a > 0 ) THEN
         332           0 :           ExitSearch = IsClosest( prefYMDhm, availYMDhm, nTime, tidx1a )
         333             :        ENDIF
         334             :
         335             :        ! When using the interpolation flag, use the first or last timestep
         336             :        ! when outside of the available date range
         337           0 :        IF ( Lct%Dct%Dta%CycleFlag == HCO_CFLAG_INTER .and. tidx1a < 0 ) THEN
         338           0 :           IF ( prefYMDhm < availYMDhm(1) ) THEN
         339           0 :              tidx1a = 1
         340           0 :           ELSE IF ( prefYMDhm > availYMDhm(nTime) ) THEN
         341           0 :              tidx1a = nTime
         342             :           ENDIF
         343             :        ENDIF
         344             :
         345             :        ! Do not continue search if data is to be interpolated and is
         346             :        ! not discontinuous (mps, 10/23/19)
         347           0 :        IF ( Lct%Dct%Dta%CycleFlag == HCO_CFLAG_INTER .and. &
         348             :             .not. Lct%Dct%Dta%Discontinuous ) THEN
         349             :           ExitSearch = .TRUE.
         350             :        ENDIF
         351             :
         352             :        ! Write to tidx1 if this is the best match.
         353           0 :        IF ( ExitSearch ) THEN
         354           0 :           tidx1 = tidx1a
         355             :
         356             :        ! -------------------------------------------------------------
         357             :        ! If search shall be continued, adjust preferred year, then
         358             :        ! month, then day to the closest available year (month, day)
         359             :        ! in the time slices, and check if this is a better match.
         360             :        ! -------------------------------------------------------------
         361             :        ELSE
         362             :
         363             :           ! Adjust year, month, and day (in this order).
         364           0 :           CNT  = 0
         365             :           DO
         366           0 :              CNT = CNT + 1
         367           0 :              IF ( ExitSearch .OR. CNT > 3 ) EXIT
         368             :
         369             :              ! Adjust prefYMDhm at the given level (1=Y, 2=M, 3=D)
         370           0 :              CALL prefYMDhm_Adjust ( nTime, availYMDhm, prefYMDhm, CNT, tidx1a )
         371             :
         372             :              ! verbose mode
         373           0 :              IF ( verb ) THEN
         374           0 :                 write(MSG,'(A30,f14.0)') 'adjusted preferred datetime: ', &
         375           0 :                      prefYMDhm
         376           0 :                 CALL HCO_MSG(HcoState%Config%Err,MSG)
         377             :              ENDIF
         378             :
         379             :              ! check for time stamp with updated date/time
         380           0 :              CALL Check_AvailYMDhm ( Lct, nTime, availYMDhm, prefYMDhm, tidx1a )
         381             :
         382             :              ! Can we leave now?
         383           0 :              ExitSearch = IsClosest( prefYMDhm, availYMDhm, nTime, tidx1a )
         384           0 :              IF ( ExitSearch ) tidx1 = tidx1a
         385             :
         386             :           ENDDO
         387             :        ENDIF
         388             :
         389             :        ! -------------------------------------------------------------
         390             :        ! If tidx1 still isn't defined, i.e. prefYMDhm is still
         391             :        ! outside the range of availYMDhm, set tidx1 to the closest
         392             :        ! available date. This must be 1 or nTime!
         393             :        ! -------------------------------------------------------------
         394           0 :        IF ( .NOT. ExitSearch ) THEN
         395           0 :           IF ( prefYMDhm < availYMDhm(1) ) THEN
         396           0 :              tidx1 = 1
         397             :           ELSE
         398           0 :              tidx1 = nTime
         399             :           ENDIF
         400             :        ENDIF
         401             :
         402             :        ! -------------------------------------------------------------
         403             :        ! If we are dealing with 3-hourly or hourly data, select all timesteps
         404             :        ! -------------------------------------------------------------
         405             :
         406             :        ! Hour flag is -1: wildcard
         407           0 :        IF ( Lct%Dct%Dta%ncHrs(1) == -1 .AND. nTime == 8 ) THEN
         408           0 :           tidx1 = 1
         409           0 :           tidx2 = nTime
         410             :
         411             :           ! verbose mode
         412           0 :           IF ( verb ) THEN
         413           0 :              WRITE(MSG,*) 'Data is 3-hourly. Entire day will be read.'
         414           0 :              CALL HCO_MSG(HcoState%Config%Err,MSG)
         415             :           ENDIF
         416             :        ENDIF
         417           0 :        IF ( Lct%Dct%Dta%ncHrs(1) == -1 .AND. nTime == 24 ) THEN
         418           0 :           tidx1 = 1
         419           0 :           tidx2 = nTime
         420             :
         421             :           ! verbose mode
         422           0 :           IF ( verb ) THEN
         423           0 :              WRITE(MSG,*) 'Data is hourly. Entire day will be read.'
         424           0 :              CALL HCO_MSG(HcoState%Config%Err,MSG)
         425             :           ENDIF
         426             :        ENDIF
         427             :
         428             :        ! -------------------------------------------------------------
         429             :        ! If we are dealing with weekday data, pick the slice to be
         430             :        ! used based on the current day of week.
         431             :        ! The ncDys flag has been set in subroutine HCO_ExtractTime
         432             :        ! (hco_tidx_mod.F90) based upon the time attributes set in the
         433             :        ! configuration file. It can have the following values:
         434             :        ! >0  : specific days are given.
         435             :        ! -1  : wildcard (autodetect)
         436             :        ! -10 : WD (weekday).
         437             :        ! -999: determine from current simulation day.
         438             :        ! For specific days or if determined from the current datetime
         439             :        ! (flags >0 or -999), the weekday is not taken into account.
         440             :        ! If auto-detection is enabled, days are treated as weekday if
         441             :        ! (and only if) there are exactly 7 time slices. Otherwise, they
         442             :        ! are interpreted as 'regular' day data.
         443             :        ! If flag is set to -10, e.g. time attribute is 'WD', the current
         444             :        ! time index is assumed to hold Sunday data, with the following
         445             :        ! six slices being Mon, Tue, ..., Sat. For weekdaily data, all
         446             :        ! seven time slices will be read into memory so that at any given
         447             :        ! time, the local weekday can be taken (weekdaily data is always
         448             :        ! assumed to be in local time).
         449             :        ! -------------------------------------------------------------
         450             :
         451             :        ! Day flag is -1: wildcard
         452           0 :        IF ( Lct%Dct%Dta%ncDys(1) == -1 .AND. nTime == 7 ) THEN
         453           0 :           tidx1 = 1
         454           0 :           tidx2 = nTime
         455             :
         456             :           ! Make sure data is treated in local time
         457           0 :           Lct%Dct%Dta%IsLocTime = .TRUE.
         458             :
         459             :        ! Day flag is -10: WD
         460           0 :        ELSEIF ( Lct%Dct%Dta%ncDys(1) == -10 ) THEN
         461             :
         462             :           ! There must be at least 7 time slices
         463           0 :           IF ( nTime < 7 ) THEN
         464             :              MSG = 'Data must have exactly 7 time slices '// &
         465           0 :                    'if you set day attribute to WD: '//TRIM(Lct%Dct%cName)
         466           0 :              CALL HCO_ERROR( MSG, RC )
         467           0 :              IF ( ASSOCIATED(availYMDhm) ) THEN
         468           0 :                 DEALLOCATE(availYMDhm)
         469             :                 availYMDhm => NULL()
         470             :              ENDIF
         471           0 :              RETURN
         472             :           ENDIF
         473             :
         474             :           ! If there are exactly seven time slices, interpret them as
         475             :           ! the seven weekdays.
         476           0 :           IF ( nTime == 7 ) THEN
         477           0 :              tidx1 = 1
         478           0 :              tidx2 = 7
         479             :
         480             :           ! If there are more than 7 time slices, interpret the current
         481             :           ! selected index as sunday of the current time frame (e.g. sunday
         482             :           ! data of current month), and select the time slice index
         483             :           ! accordingly. This requires that there are at least 6 more time
         484             :           ! slices following the current one.
         485             :           ELSE
         486           0 :              IF ( tidx1 < 0 ) THEN
         487           0 :                 WRITE(MSG,*) 'Cannot get weekday slices for: ', &
         488           0 :                    TRIM(Lct%Dct%cName), '. Cannot find first time slice.'
         489           0 :                 CALL HCO_ERROR( MSG, RC )
         490           0 :                 IF ( ASSOCIATED(availYMDhm) ) THEN
         491           0 :                    DEALLOCATE(availYMDhm)
         492             :                    availYMDhm => NULL()
         493             :                 ENDIF
         494           0 :                 RETURN
         495             :              ENDIF
         496             :
         497           0 :              IF ( (tidx1+6) > nTime ) THEN
         498           0 :                 WRITE(MSG,*) 'Cannot get weekday for: ',TRIM(Lct%Dct%cName), &
         499           0 :                    '. There are less than 6 additional time slices after ',  &
         500           0 :                    'selected start date ', availYMDhm(tidx1)
         501           0 :                 CALL HCO_ERROR( MSG, RC )
         502           0 :                 IF ( ASSOCIATED(availYMDhm) ) THEN
         503           0 :                    DEALLOCATE(availYMDhm)
         504             :                    availYMDhm => NULL()
         505             :                 ENDIF
         506           0 :                 RETURN
         507             :              ENDIF
         508           0 :              tidx2 = tidx1 + 6
         509             :           ENDIF
         510             :
         511             :           ! Make sure data is treated in local time
         512           0 :           Lct%Dct%Dta%IsLocTime = .TRUE.
         513             :
         514             :        ENDIF
         515             :
         516             :        ! -------------------------------------------------------------
         517             :        ! Now need to set upper time slice index tidx2. This index
         518             :        ! is only different from tidx1 if:
         519             :        ! (1) We interpolate between two time slices, i.e. TimeCycle
         520             :        !     attribute is set to 'I'. In this case, we simply pick
         521             :        !     the next higher time slice index and calculate the
         522             :        !     weights for time1 and time2 based on the current time.
         523             :        ! (2) Multiple hourly slices are read (--> prefHr = -1 or -10,
         524             :        !     e.g. hour attribute in config. file was set to wildcard
         525             :        !     character or data is in local hours). In this case,
         526             :        !     check if there are multiple time slices for the selected
         527             :        !     date (y/m/d).
         528             :        ! tidx2 has already been set to proper value above if it's
         529             :        ! weekday data.
         530             :        ! -------------------------------------------------------------
         531           0 :        IF ( tidx2 < 0 ) THEN
         532             :
         533             :           ! Interpolate between dates
         534           0 :           IF ( Lct%Dct%Dta%CycleFlag == HCO_CFLAG_INTER ) THEN
         535             :
         536             :              CALL GetIndex2Interp( HcoState,   Lct,       nTime,      &
         537             :                                    availYMDhm, prefYMDhm, origYMDhm,  &
         538             :                                    tidx1,      tidx2,     wgt1,       &
         539           0 :                                    wgt2,       RC                   )
         540           0 :              IF ( RC /= HCO_SUCCESS ) THEN
         541             :                 MSG = 'Error encountered in GetIndex2Interp for: '        // &
         542           0 :                      TRIM(Lct%Dct%Cname)
         543           0 :                 CALL HCO_ERROR( MSG, RC )
         544           0 :                 IF ( ASSOCIATED(availYMDhm) ) THEN
         545           0 :                    DEALLOCATE(availYMDhm)
         546             :                    availYMDhm => NULL()
         547             :                 ENDIF
         548           0 :                 RETURN
         549             :              ENDIF
         550             :
         551             :           ! Check for multiple hourly data
         552           0 :           ELSEIF ( tidx1 > 0 .AND. prefHr < 0 ) THEN
         553           0 :              CALL SET_TIDX2 ( nTime, availYMDhm, tidx1, tidx2 )
         554             :
         555             :              ! Denote as local time if necessary
         556           0 :              IF ( Lct%Dct%Dta%ncHrs(1) == -10 ) THEN
         557           0 :                 Lct%Dct%Dta%IsLocTime = .TRUE.
         558             :              ENDIF
         559             :           ELSE
         560           0 :              tidx2 = tidx1
         561             :           ENDIF
         562             :        ENDIF
         563             :
         564             :     ! ================================================================
         565             :     ! Case 3: No time slice available. Set both indeces to zero. Data
         566             :     ! with no time stamp must have CycleFlag 'Cycling'.
         567             :     ! ================================================================
         568             :     ELSE
         569           0 :        IF ( Lct%Dct%Dta%CycleFlag /= HCO_CFLAG_CYCLE ) THEN
         570             :           MSG = 'Field has no time/date variable - cycle flag must' // &
         571             :                 'be set to `C` in the HEMCO configuration file:'    // &
         572           0 :                 TRIM(Lct%Dct%cName)
         573           0 :           CALL HCO_ERROR( MSG, RC )
         574           0 :           IF ( ASSOCIATED(availYMDhm) ) THEN
         575           0 :              DEALLOCATE(availYMDhm)
         576             :              availYMDhm => NULL()
         577             :           ENDIF
         578           0 :           RETURN
         579             :        ENDIF
         580             :
         581           0 :        tidx1 = 0
         582           0 :        tidx2 = 0
         583             :     ENDIF
         584             :
         585             :     !-----------------------------------------------------------------
         586             :     ! Sanity check: if CycleFlag is set to 'Exact', the file time stamp
         587             :     ! must exactly match the current time.
         588             :     !-----------------------------------------------------------------
         589           0 :     IF ( (Lct%Dct%Dta%CycleFlag == HCO_CFLAG_EXACT) .AND. (tidx1 > 0) ) THEN
         590           0 :        IF ( availYMDhm(tidx1) /= prefYMDhm ) THEN
         591           0 :           tidx1 = -1
         592           0 :           tidx2 = -1
         593             :        ENDIF
         594             :     ENDIF
         595             :
         596             :     !-----------------------------------------------------------------
         597             :     ! If multiple time slices are read, extract time interval between
         598             :     ! time slices in memory (in hours). This is to make sure that the
         599             :     ! cycling between the slices will be done at the correct rate
         600             :     ! (e.g. every hour, every 3 hours, ...).
         601             :     !-----------------------------------------------------------------
         602           0 :     IF ( (tidx2>tidx1) .AND. (Lct%Dct%Dta%CycleFlag/=HCO_CFLAG_INTER) ) THEN
         603           0 :        Lct%Dct%Dta%DeltaT = YMDhm2hrs( availYMDhm(tidx1+1) - availYMDhm(tidx1) )
         604             :     ELSE
         605           0 :        Lct%Dct%Dta%DeltaT = 0
         606             :     ENDIF
         607             :
         608             :     ! verbose mode
         609           0 :     IF ( verb ) THEN
         610           0 :        WRITE(MSG,'(A30,I14)') 'selected tidx1: ', tidx1
         611           0 :        CALL HCO_MSG(HcoState%Config%Err,MSG)
         612           0 :        IF ( tidx1 > 0 ) THEN
         613           0 :           WRITE(MSG,'(A30,f14.0)') 'corresponding datetime 1: ', &
         614           0 :                availYMDhm(tidx1)
         615           0 :           CALL HCO_MSG(HcoState%Config%Err,MSG)
         616           0 :           IF ( wgt1 >= 0.0_sp ) THEN
         617           0 :              WRITE(MSG,*) 'weight1: ', wgt1
         618           0 :              CALL HCO_MSG(HcoState%Config%Err,MSG)
         619             :           ENDIF
         620             :        ENDIF
         621             :
         622           0 :        IF ( (tidx2 /= tidx1) ) THEN
         623           0 :           WRITE(MSG,'(A30,I14)') 'selected tidx2: ', tidx2
         624           0 :           CALL HCO_MSG(HcoState%Config%Err,MSG)
         625           0 :           WRITE(MSG,'(A30,f14.0)') 'corresponding datetime 2: ', &
         626           0 :                availYMDhm(tidx2)
         627           0 :           CALL HCO_MSG(HcoState%Config%Err,MSG)
         628           0 :           IF ( wgt1 >= 0.0_sp ) THEN
         629           0 :              WRITE(MSG,*) 'weight2: ', wgt2
         630           0 :              CALL HCO_MSG(HcoState%Config%Err,MSG)
         631             :           ENDIF
         632             :        ENDIF
         633             :
         634           0 :        WRITE(MSG,'(A30,I14)') 'assigned delta t [h]: ', Lct%Dct%Dta%DeltaT
         635           0 :        CALL HCO_MSG(HcoState%Config%Err,MSG)
         636           0 :        WRITE(MSG,*) 'local time? ', Lct%Dct%Dta%IsLocTime
         637           0 :        CALL HCO_MSG(HcoState%Config%Err,MSG)
         638             :     ENDIF
         639             :
         640             :     ! ----------------------------------------------------------------
         641             :     ! TODO: set time brackets
         642             :     ! --> In future, we may want to set time brackets denoting the
         643             :     ! previous and next time slice available in the netCDF file. This
         644             :     ! may become useful for temporal interpolations and more efficient
         645             :     ! data update calls (only update if new time slice is available).
         646             :     ! ----------------------------------------------------------------
         647             :
         648             :     !-----------------------------------------------------------------
         649             :     ! Prepare output, cleanup and leave
         650             :     !-----------------------------------------------------------------
         651             :
         652             :     ! ncYr and ncMt are the year and month fo the time slice to be
         653             :     ! used. These values may be required to convert units to 'per
         654             :     ! seconds'.
         655           0 :     IF ( tidx1 > 0 ) THEN
         656           0 :        YMDhm  = availYMDhm(tidx1)
         657           0 :        YMDhm1 = availYMDhm(1)
         658           0 :        oYMDhm = origYMDhm
         659             :     ENDIF
         660             :
         661             :     ! Deallocate and nullify the pointer
         662           0 :     IF ( ASSOCIATED(availYMDhm) ) THEN
         663           0 :        DEALLOCATE(availYMDhm)
         664             :        availYMDhm => NULL()
         665             :     ENDIF
         666             :
         667             :     ! Return w/ success
         668           0 :     CALL HCO_LEAVE ( HcoState%Config%Err,  RC )
         669             :
         670           0 :   END SUBROUTINE GET_TIMEIDX
         671             : !EOC
         672             : !------------------------------------------------------------------------------
         673             : !                   Harmonized Emissions Component (HEMCO)                    !
         674             : !------------------------------------------------------------------------------
         675             : !BOP
         676             : !
         677             : ! !IROUTINE: Check_AvailYMDhm
         678             : !
         679             : ! !DESCRIPTION: Checks if prefYMDhm is within the range of availYMDhm
         680             : ! and returns the location of the closest vector element that is in
         681             : ! the past (--> tidx1). tidx1 is set to -1 otherwise.
         682             : !\\
         683             : !\\
         684             : ! !INTERFACE:
         685             : !
         686           0 :   SUBROUTINE Check_AvailYMDhm( Lct, N, availYMDhm, prefYMDhm, tidx1 )
         687             : !
         688             : ! !INPUT PARAMETERS:
         689             : !
         690             :     TYPE(ListCont),   POINTER      :: Lct
         691             :     INTEGER,          INTENT(IN)   :: N
         692             :     REAL(dp),         INTENT(IN)   :: availYMDhm(N)
         693             :     REAL(dp),         INTENT(IN)   :: prefYMDhm
         694             : !
         695             : ! !OUTPUT PARAMETERS:
         696             : !
         697             :     INTEGER,          INTENT(OUT)  :: tidx1
         698             : !
         699             : ! !REVISION HISTORY:
         700             : !  13 Mar 2013 - C. Keller   - Initial version
         701             : !  See https://github.com/geoschem/hemco for complete history
         702             : !EOP
         703             : !------------------------------------------------------------------------------
         704             : !BOC
         705             : !
         706             : ! !LOCAL VARIABLES:
         707             : !
         708             :     INTEGER :: I, nTime
         709             :
         710             :     !=================================================================
         711             :     ! Check_availYMDhm begins here
         712             :     !=================================================================
         713             :
         714             :     ! Init
         715           0 :     tidx1 = -1
         716             :
         717             :     ! Return if preferred datetime not within the vector range
         718           0 :     IF ( prefYMDhm < availYMDhm(1) .OR. prefYMDhm > availYMDhm(N) ) RETURN
         719             :
         720             :     ! To avoid out-of-bounds error in the loop below:
         721             :     ! (1) For interpolated data, the upper loop limit should be N;
         722             :     ! (2) Otherwise, the upper loop limit should be N-1.
         723             :     ! (bmy, 4/28/21)
         724           0 :     nTime = N - 1
         725           0 :     IF ( Lct%Dct%Dta%CycleFlag == HCO_CFLAG_INTER ) nTime = N
         726             :
         727             :     ! Get closest index that is not in the future
         728           0 :     DO I = 1, nTime
         729             :
         730             :        ! NOTE: Epsilon test is more robust than an equality test
         731             :        ! for double-precision variables (bmy, 4/11/17)
         732           0 :        IF ( ABS( availYMDhm(I) - prefYMDhm ) < EPSILON ) THEN
         733           0 :           tidx1 = I
         734           0 :           EXIT
         735             :        ENDIF
         736             :
         737             :        ! Check if next time slice is in the future, in which case the
         738             :        ! current slice is selected. Don't do this for a CycleFlag of
         739             :        ! 3 (==> exact match).
         740           0 :        IF ( (availYMDhm(I+1)       >  prefYMDhm        ) .AND. &
         741           0 :             (Lct%Dct%Dta%CycleFlag /= HCO_CFLAG_EXACT) ) THEN
         742           0 :           tidx1 = I
         743           0 :           EXIT
         744             :        ENDIF
         745             :     ENDDO
         746             :
         747             :   END SUBROUTINE Check_AvailYMDhm
         748             : !EOC
         749             : !------------------------------------------------------------------------------
         750             : !                   Harmonized Emissions Component (HEMCO)                    !
         751             : !------------------------------------------------------------------------------
         752             : !BOP
         753             : !
         754             : ! !IROUTINE: prefYMDhm_Adjust
         755             : !
         756             : ! !DESCRIPTION: Adjusts prefYMDhm to the closest available time attribute. Can
         757             : ! be adjusted for year (level=1), month (level=2), or day (level=3).
         758             : !\\
         759             : !\\
         760             : ! !INTERFACE:
         761             : !
         762           0 :   SUBROUTINE prefYMDhm_Adjust( N, availYMDhm, prefYMDhm, level, tidx1 )
         763             : !
         764             : ! !INPUT PARAMETERS:
         765             : !
         766             :     INTEGER   , INTENT(IN)     :: N
         767             :     REAL(dp)  , INTENT(IN)     :: availYMDhm(N)
         768             :     INTEGER   , INTENT(IN)     :: level
         769             :     INTEGER   , INTENT(IN)     :: tidx1
         770             : !
         771             : ! !INPUT/OUTPUT PARAMETERS:
         772             : !
         773             :     REAL(dp)  , INTENT(INOUT)  :: prefYMDhm
         774             : !
         775             : ! !REVISION HISTORY:
         776             : !  13 Mar 2013 - C. Keller   - Initial version
         777             : !  See https://github.com/geoschem/hemco for complete history
         778             : !EOP
         779             : !------------------------------------------------------------------------------
         780             : !BOC
         781             : !
         782             : ! !LOCAL VARIABLES:
         783             : !
         784             :     ! Scalars
         785             :     INTEGER          :: I, IMIN, IMAX
         786             :     REAL(dp)         :: origYr,  origMt,  origDy, origHr, origMi
         787             :     REAL(dp)         :: refAttr, tmpAttr, newAttr
         788             :     REAL(dp)         :: iDiff,   minDiff
         789             :     REAL(dp)         :: modVal
         790             :     REAL(dp)         :: div
         791             :
         792             :     !=================================================================
         793             :     ! prefYMDhm_Adjust begins here!
         794             :     !=================================================================
         795             :
         796             :     ! Get original Yr, Mt, Day, Hr, Mi
         797             :     ! Time values are now in YYYYMMDDhhmm format (bmy, 4/11/17)
         798           0 :     origYr = FLOOR( MOD( prefYMDhm, 1.0e12_dp ) / 1.0e8_dp )
         799           0 :     origMt = FLOOR( MOD( prefYMDhm, 1.0e8_dp  ) / 1.0e6_dp )
         800           0 :     origDy = FLOOR( MOD( prefYMDhm, 1.0e6_dp  ) / 1.0e4_dp )
         801           0 :     origHr = FLOOR( MOD( prefYMDhm, 1.0e4_dp  ) / 1.0e2_dp )
         802           0 :     origMi = FLOOR( MOD( prefYMDhm, 1.0e2_dp  )            )
         803             :
         804             :     ! Extract new attribute from availYMDhm and insert into prefYMDhm. Pick
         805             :     ! closest available value.
         806           0 :     SELECT CASE ( level )
         807             :        ! --- Year
         808             :        CASE ( 1 )
         809             :           modVal  = 1.0e12_dp
         810             :           div     = 1.0e8_dp
         811           0 :           refAttr = origYr
         812             :
         813             :        ! --- Month
         814             :        CASE ( 2 )
         815           0 :           modVal  = 1.0e8_dp
         816           0 :           div     = 1.0e6_dp
         817           0 :           refAttr = origMt
         818             :
         819             :        ! --- Day
         820             :        CASE ( 3 )
         821           0 :           modVal  = 1.0e6_dp
         822           0 :           div     = 1.0e4_dp
         823           0 :           refAttr = origMt
         824             :
         825             :        ! --- Hour
         826             :        CASE ( 4 )
         827           0 :           modval  = 1.0e4_dp
         828           0 :           div     = 1.0e2_dp
         829           0 :           refAttr = origHr
         830             :
         831             :        ! --- Minute
         832             :        CASE ( 5 )
         833           0 :           modVal  = 1.0e2_dp
         834           0 :           div     = 1.0_dp
         835           0 :           refAttr = origMi
         836             :
         837             :        CASE DEFAULT
         838           0 :           RETURN
         839             :     END SELECT
         840             :
         841             :     ! Maximum loop number:
         842             :     ! If tidx1 is already set, only search values in the past.
         843           0 :     IF ( tidx1 > 0 ) THEN
         844             :        IMIN = 1
         845             :        IMAX = tidx1
         846             :
         847             :     ! If tidx1 is not yet set, prefYMDhm must be outside the range of
         848             :     ! availYMDhm. Pick only the closest available time stamp.
         849             :     ELSE
         850           0 :        IF ( prefYMDhm > availYMDhm(1) ) THEN
         851           0 :           IMIN = N
         852           0 :           IMAX = N
         853             :        ELSE
         854             :           IMIN = 1
         855             :           IMAX = 1
         856             :        ENDIF
         857             :     ENDIF
         858             :
         859             :     ! Select current minimum value
         860           0 :     minDiff = 10000000000000000.0_dp
         861           0 :     newAttr = -1d0
         862           0 :     DO I = IMIN, IMAX
         863           0 :        tmpAttr = FLOOR( MOD(availYMDhm(I),modVal) / div )
         864           0 :        iDiff   = ABS( tmpAttr - refAttr )
         865           0 :        IF ( iDiff < minDiff ) THEN
         866           0 :           newAttr = tmpAttr
         867           0 :           minDiff = iDiff
         868             :        ENDIF
         869             :     ENDDO
         870             :
         871             :     ! Just reuse current value if no better value could be found
         872           0 :     IF ( newAttr < 0 ) THEN
         873           0 :        newAttr = refAttr
         874             :     ENDIF
         875             :
         876             :     ! Update variable
         877             :     ! --- Year
         878           0 :     IF ( level == 1 ) THEN
         879             :        prefYMDhm = ( newAttr * 1.0e8_dp ) + &
         880             :                    ( origMt  * 1.0e6_dp ) + &
         881             :                    ( origDy  * 1.0e4_dp ) + &
         882             :                    ( origHr  * 1.0e2_dp ) + &
         883           0 :                    ( origMi          )
         884             :
         885             :     ! --- Month
         886           0 :     ELSEIF ( level == 2 ) THEN
         887             :        prefYMDhm = ( origYr  * 1.0e8_dp ) + &
         888             :                    ( newAttr * 1.0e6_dp ) + &
         889             :                    ( origDy  * 1.0e4_dp ) + &
         890             :                    ( origHr  * 1.0e2_dp ) + &
         891           0 :                    ( origMi             )
         892             :
         893             :     ! --- Day
         894           0 :     ELSEIF ( level == 3 ) THEN
         895             :        prefYMDhm = ( origYr  * 1.0e8_dp  ) + &
         896             :                    ( origMt  * 1.0e6_dp  ) + &
         897             :                    ( newAttr * 1.0e4_dp  ) + &
         898             :                    ( origHr  * 1.0e2_dp  ) + &
         899           0 :                    ( origMi              )
         900             :
         901             :     ! --- Hour
         902           0 :     ELSEIF ( level == 4 ) THEN
         903             :        prefYMDhm = ( origYr  * 1.0e8_dp  ) + &
         904             :                    ( origMt  * 1.0e6_dp  ) + &
         905             :                    ( origDy  * 1.0e4_dp  ) + &
         906             :                    ( newAttr * 1.0e2_dp  ) + &
         907           0 :                    ( origMi              )
         908             :     ! --- Minute
         909           0 :     ELSEIF ( level == 5 ) THEN
         910             :        prefYMDhm = ( origYr  * 1.0e8_dp  ) + &
         911             :                    ( origMt  * 1.0e6_dp  ) + &
         912             :                    ( origDy  * 1.0e4_dp  ) + &
         913             :                    ( origHr  * 1.0e2_dp  ) + &
         914           0 :                    ( newAttr             )
         915             :
         916             :     ENDIF
         917             :
         918             :   END SUBROUTINE prefYMDhm_Adjust
         919             : !EOC
         920             : !------------------------------------------------------------------------------
         921             : !                   Harmonized Emissions Component (HEMCO)                    !
         922             : !------------------------------------------------------------------------------
         923             : !BOP
         924             : !
         925             : ! !IROUTINE: Set_tIdx2
         926             : !
         927             : ! !DESCRIPTION: sets the upper time slice index by selecting the range
         928             : ! of all elements in availYMDhm with the same date (year,month,day) as
         929             : ! availYMDh(tidx1).
         930             : !\\
         931             : !\\
         932             : ! !INTERFACE:
         933             : !
         934           0 :   SUBROUTINE Set_tIdx2( N, availYMDhm, tidx1, tidx2 )
         935             : !
         936             : ! !INPUT PARAMETERS:
         937             : !
         938             :     INTEGER,  INTENT(IN)  :: N               ! Number of times
         939             :     REAL(dp), INTENT(IN)  :: availYMDhm(N)   ! Time stamp vector
         940             :     INTEGER,  INTENT(IN)  :: tidx1           ! Lower time slice index
         941             : !
         942             : ! !INPUT/OUTPUT PARAMETERS:
         943             : !
         944             :     INTEGER,  INTENT(OUT) :: tidx2           ! Upper time slice index
         945             : !
         946             : ! !REVISION HISTORY:
         947             : !  13 Mar 2013 - C. Keller   - Initial version
         948             : !  See https://github.com/geoschem/hemco for complete history
         949             : !EOP
         950             : !------------------------------------------------------------------------------
         951             : !BOC
         952             : !
         953             : ! !LOCAL VARIABLES:
         954             : !
         955             :     INTEGER :: YMD, I, IYMD
         956             :
         957             :     !=================================================================
         958             :     ! SET_TIDX2 begins here!
         959             :     !=================================================================
         960             :
         961             :     ! Init
         962           0 :     tidx2 = tidx1
         963             :
         964             :     ! Sanity check
         965           0 :     IF ( tidx1 == N ) RETURN
         966             :
         967             :     ! Get wanted YMD
         968           0 :     YMD = floor(availYMDhm(tidx1) / 1.0e4_dp)
         969             :
         970             :     ! See how many more tile slices with the same YMD exist from index
         971             :     ! tidx1 onwards.
         972           0 :     DO I = tidx1, N
         973           0 :        iYMD = floor(availYMDhm(I) / 1.0e4_dp)
         974           0 :        IF ( iYMD == YMD ) THEN
         975           0 :           tidx2 = I
         976           0 :        ELSEIF ( iYMD > YMD ) THEN
         977             :           EXIT
         978             :        ENDIF
         979             :     ENDDO
         980             :
         981             :   END SUBROUTINE Set_tIdx2
         982             : !EOC
         983             : !------------------------------------------------------------------------------
         984             : !                   Harmonized Emissions Component (HEMCO)                    !
         985             : !------------------------------------------------------------------------------
         986             : !BOP
         987             : !
         988             : ! !IROUTINE: IsClosest
         989             : !
         990             : ! !DESCRIPTION: function IsClosest returns true if the selected time index
         991             : ! is the 'closest' one. It is defined as being closest if:
         992             : ! (a) the currently selected index exactly matches the preferred one.
         993             : ! (b) the time gap between the preferred time stamp and the currently selected
         994             : ! index is at least as small as any other gap of consecutive prior time stamps.
         995             : !\\
         996             : !\\
         997             : ! !INTERFACE:
         998             : !
         999           0 :   FUNCTION IsClosest ( prefYMDhm, availYMDhm, nTime, ctidx1 ) RESULT ( Closest )
        1000             : !
        1001             : ! !INPUT PARAMETERS:
        1002             : !
        1003             :     INTEGER,    INTENT(IN)  :: nTime
        1004             :     REAL(dp),   INTENT(IN)  :: prefYMDhm
        1005             :     REAL(dp),   INTENT(IN)  :: availYMDhm(nTime)
        1006             :     INTEGER,    INTENT(IN)  :: ctidx1
        1007             : !
        1008             : ! !OUTPUT PARAMETERS:
        1009             : !
        1010             :     LOGICAL              :: Closest
        1011             : !
        1012             : ! !REVISION HISTORY:
        1013             : !  03 Mar 2015 - C. Keller   - Initial version
        1014             : !  See https://github.com/geoschem/hemco for complete history
        1015             : !EOP
        1016             : !------------------------------------------------------------------------------
        1017             : !BOC
        1018             : !
        1019             : ! !LOCAL VARIABLES:
        1020             : !
        1021             :     INTEGER :: N
        1022             :     INTEGER :: diff, idiff
        1023             :
        1024             :     !=================================================================
        1025             :     ! IsClosest begins here!
        1026             :     !=================================================================
        1027             :
        1028             :     ! Init
        1029           0 :     Closest = .TRUE.
        1030             :
        1031             :     ! It's not closest if index is not defined
        1032           0 :     IF ( ctidx1 <= 0 ) THEN
        1033           0 :        Closest = .FALSE.
        1034             :        RETURN
        1035             :     ENDIF
        1036             :
        1037             :     ! It's closest if it is the first index
        1038           0 :     IF ( ctidx1 == 1 ) RETURN
        1039             :
        1040             :     ! It's closest if it matches date exactly
        1041             :     ! NOTE: Epsilon test is more robust than an equality test
        1042             :     ! for double-precision variables (bmy, 4/11/17)
        1043           0 :     IF ( ABS( availYMDhm(ctidx1) - prefYMDhm ) < EPSILON ) RETURN
        1044             :
        1045             :     ! It's closest if current select one is in the future
        1046           0 :     IF ( availYMDhm(ctidx1) > prefYMDhm ) RETURN
        1047             :
        1048             :     ! Check if any of the time stamps in the past have closer intervals
        1049             :     ! than the current select time stamp to it's previous one
        1050           0 :     diff = prefYMDhm - availYMDhm(ctidx1)
        1051           0 :     DO N = 2, ctidx1
        1052           0 :        idiff = availYMDhm(N) - availYMDhm(N-1)
        1053           0 :        IF ( idiff < diff ) THEN
        1054           0 :           Closest = .FALSE.
        1055             :           RETURN
        1056             :        ENDIF
        1057             :     ENDDO
        1058             :
        1059             :   END FUNCTION IsClosest
        1060             : !EOC
        1061             : !------------------------------------------------------------------------------
        1062             : !                   Harmonized Emissions Component (HEMCO)                    !
        1063             : !------------------------------------------------------------------------------
        1064             : !BOP
        1065             : !
        1066             : ! !IROUTINE: GetIndex2Interp
        1067             : !
        1068             : ! !DESCRIPTION: GetIndex2Interp
        1069             : !\\
        1070             : !\\
        1071             : ! !INTERFACE:
        1072             : !
        1073           0 :   SUBROUTINE GetIndex2Interp ( HcoState,  Lct,                   &
        1074           0 :                                nTime,     availYMDhm,            &
        1075             :                                prefYMDhm, origYMDhm,  tidx1,     &
        1076             :                                tidx2,     wgt1,       wgt2,  RC   )
        1077             : !
        1078             : ! !INPUT PARAMETERS:
        1079             : !
        1080             :     TYPE(HCO_State),  POINTER       :: HcoState
        1081             :     TYPE(ListCont),   POINTER       :: Lct
        1082             :     INTEGER,          INTENT(IN)    :: nTime
        1083             :     REAL(dp),         INTENT(IN)    :: availYMDhm(nTime)
        1084             :     REAL(dp),         INTENT(IN)    :: prefYMDhm
        1085             :     REAL(dp),         INTENT(IN)    :: origYMDhm
        1086             :     INTEGER,          INTENT(IN)    :: tidx1
        1087             : !
        1088             : ! !OUTPUT PARAMETERS:
        1089             : !
        1090             :     INTEGER,          INTENT(OUT)   :: tidx2
        1091             : !
        1092             : ! !INPUT/OUTPUT PARAMETERS:
        1093             : !
        1094             :     REAL(sp),         INTENT(INOUT) :: wgt1
        1095             :     REAL(sp),         INTENT(INOUT) :: wgt2
        1096             :     INTEGER,          INTENT(INOUT) :: RC
        1097             : !
        1098             : ! !REVISION HISTORY:
        1099             : !  02 Mar 2015 - C. Keller   - Initial version
        1100             : !  See https://github.com/geoschem/hemco for complete history
        1101             : !EOP
        1102             : !------------------------------------------------------------------------------
        1103             : !BOC
        1104             : !
        1105             : ! !LOCAL VARIABLES:
        1106             : !
        1107             :     ! Scalars
        1108             :     INTEGER             :: I
        1109             :     REAL(dp)            :: tmpYMDhm
        1110             :     LOGICAL             :: verb
        1111             :
        1112             :     ! Strings
        1113             :     CHARACTER(LEN=255)  :: MSG
        1114             :     CHARACTER(LEN=255)  :: LOC = 'GetIndex2Interp (hcoio_util_mod.F90)'
        1115             :
        1116             :     !=================================================================
        1117             :     ! GetIndex2Interp begins here
        1118             :     !=================================================================
        1119             :
        1120             :     ! Verbose mode?
        1121             :     verb = HCO_IsVerb( HcoState%Config%Err )
        1122             :
        1123             :     ! If the originally wanted datetime was below the available data
        1124             :     ! range, set all weights to the first index.
        1125           0 :     IF ( origYMDhm <= availYMDhm(1) ) THEN
        1126           0 :        tidx2 = tidx1
        1127           0 :        wgt1  = 1.0_sp
        1128           0 :        wgt2  = 0.0_sp
        1129             :
        1130             :     ! If the originally wanted datetime is beyond the available data
        1131             :     ! range, set tidx2 to tidx1 but leave weights in their original
        1132             :     ! values (-1.0). The reason is that we will attempt to interpolate
        1133             :     ! between a second file, which is only done if the weights are
        1134             :     ! negative.
        1135           0 :     ELSEIF ( origYMDhm >= availYMDhm(nTime) ) THEN
        1136           0 :        tidx2 = tidx1
        1137             :
        1138             :     ! No interpolation needed if there is a time slices that exactly
        1139             :     ! matches the (originally) preferred datetime.
        1140             :     ! NOTE: An Epsilon test is more robust than an equality test
        1141             :     ! for double-precision variables (bmy, 4/11/17)
        1142           0 :     ELSEIF ( ABS( origYMDhm - availYMDhm(tidx1) ) < EPSILON ) THEN
        1143           0 :        tidx2 = tidx1
        1144           0 :        wgt1  = 1.0_sp
        1145           0 :        wgt2  = 0.0_sp
        1146             :
        1147             :     ! If we are inside the data range but none of the time slices
        1148             :     ! matches the preferred datetime, get the second time slices that
        1149             :     ! shall be used for data interpolation. This not necessarily needs
        1150             :     ! to be the consecutive time slice. For instance, imagine a data
        1151             :     ! set that contains montlhly data for years 2005 and 2010. For
        1152             :     ! Feb 2007, we would want to interpolate between Feb 2005 and Feb
        1153             :     ! 2010 data. The index tidx1 already points to Feb 2005, but the
        1154             :     ! upper index tidx2 needs to be set accordingly.
        1155             :     ELSE
        1156             :
        1157             :        ! Init
        1158           0 :        tidx2 = -1
        1159             :
        1160             :        ! Search for a time slice in the future that has the same
        1161             :        ! month/day/hour as currently selected time slice.
        1162             :        tmpYMDhm = availYMDhm(tidx1)
        1163             :        DO
        1164             :           ! Increase by one year
        1165           0 :           tmpYMDhm = tmpYMDhm + 1.0e8_dp
        1166             :
        1167             :           ! Exit if we are beyond available dates
        1168           0 :           IF ( tmpYMDhm > availYMDhm(nTime) ) EXIT
        1169             :
        1170             :           ! Check if there is a time slice with that date
        1171           0 :           DO I = tidx1,nTime
        1172           0 :              IF ( tmpYMDhm == availYMDhm(I) ) THEN
        1173           0 :                 tidx2 = I
        1174           0 :                 EXIT
        1175             :              ENDIF
        1176             :           ENDDO
        1177           0 :           IF ( tidx2 > 0 ) EXIT
        1178             :        ENDDO
        1179             :
        1180             :        ! Repeat above but now only modify day
        1181           0 :        IF ( tidx2 < 0 ) THEN
        1182             :           tmpYMDhm = availYMDhm(tidx1)
        1183             :           DO
        1184             :              ! Increase by one day
        1185           0 :              tmpYMDhm = tmpYMDhm + 1.0e4_dp
        1186             :
        1187             :              ! Exit if we are beyond available dates
        1188           0 :              IF ( tmpYMDhm > availYMDhm(nTime) ) EXIT
        1189             :
        1190             :              ! Check if there is a time slice with that date
        1191           0 :              DO I = tidx1,nTime
        1192           0 :                 IF ( tmpYMDhm == availYMDhm(I) ) THEN
        1193           0 :                    tidx2 = I
        1194           0 :                    EXIT
        1195             :                 ENDIF
        1196             :              ENDDO
        1197           0 :              IF ( tidx2 > 0 ) EXIT
        1198             :           ENDDO
        1199             :        ENDIF
        1200             :
        1201             :        ! If all of those tests failed, simply get the next time
        1202             :        ! slice.
        1203           0 :        IF ( tidx2 < 0 ) THEN
        1204           0 :           tidx2 = tidx1 + 1
        1205             :
        1206             :           ! Make sure that tidx2 does not exceed nTime, which is
        1207             :           ! the number of time slices in the file. This can cause
        1208             :           ! an out-of-bounds error. (bmy, 3/7/19)
        1209           0 :           IF ( tidx2 > nTime ) tidx2 = nTime
        1210             :
        1211             :           ! Prompt warning
        1212           0 :           WRITE(MSG,*) 'Having problems in finding the next time slice ', &
        1213           0 :                 'to interpolate from, just take the next available ',     &
        1214           0 :                 'slice. Interpolation will be performed from ',           &
        1215           0 :                 availYMDhm(tidx1), ' to ', availYMDhm(tidx2), '. Data ',    &
        1216           0 :                 'container: ', TRIM(Lct%Dct%cName)
        1217           0 :           CALL HCO_WARNING(HcoState%Config%Err, MSG, RC, THISLOC=LOC)
        1218             :        ENDIF
        1219             :
        1220             :        ! Calculate weights wgt1 and wgt2 to be given to slice 1 and
        1221             :        ! slice2, respectively.
        1222           0 :        CALL GetWeights ( availYMDhm(tidx1), availYMDhm(tidx2), origYMDhm, &
        1223           0 :                          wgt1, wgt2 )
        1224             :
        1225             :     ENDIF
        1226             :
        1227             :     ! Return w/ success
        1228           0 :     RC = HCO_SUCCESS
        1229             :
        1230           0 :   END SUBROUTINE GetIndex2Interp
        1231             : !EOC
        1232             : !------------------------------------------------------------------------------
        1233             : !                   Harmonized Emissions Component (HEMCO)                    !
        1234             : !------------------------------------------------------------------------------
        1235             : !BOP
        1236             : !
        1237             : ! !IROUTINE: GetWeights
        1238             : !
        1239             : ! !DESCRIPTION: Helper function to get the interpolation weights between
        1240             : ! two datetime intervals (int1, int2) and for a given time cur.
        1241             : !\\
        1242             : !\\
        1243             : ! !INTERFACE:
        1244             : !
        1245           0 :   SUBROUTINE GetWeights ( int1, int2, cur, wgt1, wgt2 )
        1246             : !
        1247             : ! !INPUT PARAMETERS:
        1248             : !
        1249             :     REAL(dp),         INTENT(IN   )   :: int1, int2, cur
        1250             : !
        1251             : ! !INPUT/OUTPUT PARAMETERS:
        1252             : !
        1253             :     REAL(sp),         INTENT(  OUT)   :: wgt1, wgt2
        1254             : !
        1255             : ! !REVISION HISTORY:
        1256             : !  04 Mar 2015 - C. Keller - Initial version
        1257             : !  See https://github.com/geoschem/hemco for complete history
        1258             : !EOP
        1259             : !------------------------------------------------------------------------------
        1260             : !BOC
        1261             : !
        1262             : ! !LOCAL VARIABLES:
        1263             : !
        1264             :     REAL(dp)              :: diff1, diff2
        1265             :     REAL(dp)              :: jdc, jd1, jd2
        1266             :
        1267             :     !=================================================================
        1268             :     ! GetWeights begins here!
        1269             :     !=================================================================
        1270             :
        1271             :     ! Convert dates to Julian dates
        1272           0 :     jdc = YMDhm2jd ( cur  )
        1273           0 :     jd1 = YMDhm2jd ( int1 )
        1274           0 :     jd2 = YMDhm2jd ( int2 )
        1275             :
        1276             :     ! Check if outside of range
        1277           0 :     IF ( jdc <= jd1 ) THEN
        1278           0 :        wgt1 = 1.0_sp
        1279           0 :     ELSEIF ( jdc >= jd2 ) THEN
        1280           0 :        wgt1 = 0.0_sp
        1281             :     ELSE
        1282           0 :        diff1 = jd2 - jdc
        1283           0 :        diff2 = jd2 - jd1
        1284           0 :        wgt1  = diff1 / diff2
        1285             :     ENDIF
        1286             :
        1287             :     ! second weight is just complement of wgt1
        1288           0 :     wgt2  = 1.0_sp - wgt1
        1289             :
        1290           0 :   END SUBROUTINE GetWeights
        1291             : !EOC
        1292             : !------------------------------------------------------------------------------
        1293             : !                   Harmonized Emissions Component (HEMCO)                    !
        1294             : !------------------------------------------------------------------------------
        1295             : !BOP
        1296             : !
        1297             : ! !IROUTINE: YMDhm2jd
        1298             : !
        1299             : ! !DESCRIPTION: returns the julian date of element YMDhm.
        1300             : !\\
        1301             : !\\
        1302             : ! !INTERFACE:
        1303             : !
        1304           0 :   FUNCTION YMDhm2jd ( YMDhm ) RESULT ( jd )
        1305             : !
        1306             : ! !USES:
        1307             : !
        1308             :     USE HCO_Julday_Mod
        1309             : !
        1310             : ! !INPUT PARAMETERS:
        1311             : !
        1312             :     REAL(dp), INTENT(IN)  :: YMDhm
        1313             : !
        1314             : ! !INPUT/OUTPUT PARAMETERS:
        1315             : !
        1316             :     REAL(hp) :: jd
        1317             : !
        1318             : ! !REVISION HISTORY:
        1319             : !  24 Feb 2019 - C. Keller - Initial version
        1320             : !  See https://github.com/geoschem/hemco for complete history
        1321             : !EOP
        1322             : !------------------------------------------------------------------------------
        1323             : !BOC
        1324             : !
        1325             : ! !LOCAL VARIABLES:
        1326             : !
        1327             :     INTEGER               :: yr, mt, dy, hr, mn
        1328             :     REAL(dp)              :: utc, day
        1329             :
        1330             :     !=================================================================
        1331             :     ! YMDh2jd begins here!
        1332             :     !=================================================================
        1333           0 :     yr  = FLOOR( MOD( YMDhm, 1.0e12_dp ) / 1.0e8_dp )
        1334           0 :     mt  = FLOOR( MOD( YMDhm, 1.0e8_dp  ) / 1.0e6_dp )
        1335           0 :     dy  = FLOOR( MOD( YMDhm, 1.0e6_dp  ) / 1.0e4_dp )
        1336           0 :     hr  = FLOOR( MOD( YMDhm, 1.0e4_dp  ) / 1.0e2_dp )
        1337           0 :     mn  = FLOOR( MOD( YMDhm, 1.0e2_dp  ) )
        1338             :     utc = ( REAL(hr,dp) / 24.0_dp    ) + &
        1339             :           ( REAL(mn,dp) / 1440.0_dp  ) + &
        1340           0 :           ( REAL(0 ,dp) / 86400.0_dp )
        1341           0 :     day = REAL(dy,dp) + utc
        1342           0 :     jd  = JULDAY( yr, mt, day )
        1343             :
        1344           0 :   END FUNCTION YMDhm2jd
        1345             : !EOC
        1346             : !------------------------------------------------------------------------------
        1347             : !                   Harmonized Emissions Component (HEMCO)                    !
        1348             : !------------------------------------------------------------------------------
        1349             : !BOP
        1350             : !
        1351             : ! !IROUTINE: YMDhm2hrs
        1352             : !
        1353             : ! !DESCRIPTION: returns the hours of element YMDhm. For simplicity, 30 days are
        1354             : ! assigned to every month. At the moment, this routine is only called to
        1355             : ! determine the time interval between two emission time slices (DeltaT) and
        1356             : ! this approximation is good enough.
        1357             : !\\
        1358             : !\\
        1359             : ! !INTERFACE:
        1360             : !
        1361           0 :   FUNCTION YMDhm2hrs ( YMDhm ) RESULT ( hrs )
        1362             : !
        1363             : ! !INPUT PARAMETERS:
        1364             : !
        1365             :     REAL(dp), INTENT(IN)  :: YMDhm
        1366             : !
        1367             : ! !INPUT/OUTPUT PARAMETERS:
        1368             : !
        1369             :     INTEGER              :: hrs
        1370             : !
        1371             : ! !REVISION HISTORY:
        1372             : !  26 Jan 2015 - C. Keller - Initial version
        1373             : !  See https://github.com/geoschem/hemco for complete history
        1374             : !EOP
        1375             : !------------------------------------------------------------------------------
        1376             : !BOC
        1377             :
        1378             :     !=================================================================
        1379             :     ! YMDh2hrs begins here!
        1380             :     !=================================================================
        1381             :     hrs = FLOOR( MOD( YMDhm, 1.0e12_dp ) / 1.0e8_dp ) * 8760 + &
        1382             :           FLOOR( MOD( YMDhm, 1.0e8_dp  ) / 1.0e6_dp ) * 720  + &
        1383             :           FLOOR( MOD( YMDhm, 1.0e6_dp  ) / 1.0e4_dp ) * 24   + &
        1384           0 :           FLOOR( MOD( YMDhm, 1.0e4_dp  ) / 1.0e2_dp )
        1385             :
        1386           0 :   END FUNCTION YMDhm2hrs
        1387             : !EOC
        1388             : !------------------------------------------------------------------------------
        1389             : !                   Harmonized Emissions Component (HEMCO)                    !
        1390             : !------------------------------------------------------------------------------
        1391             : !BOP
        1392             : !
        1393             : ! !IROUTINE: Normalize_Area
        1394             : !
        1395             : ! !DESCRIPTION: Subroutine Normalize\_Area normalizes the given array
        1396             : ! by the surface area calculated from the given netCDF file.
        1397             : !\\
        1398             : !\\
        1399             : ! !INTERFACE:
        1400             : !
        1401           0 :   SUBROUTINE Normalize_Area( HcoState, Array, nlon, LatEdge, FN, RC )
        1402             : !
        1403             : ! !INPUT PARAMETERS:
        1404             : !
        1405             :     TYPE(HCO_State),  POINTER         :: HcoState    ! HEMCO state object
        1406             :     INTEGER,          INTENT(IN   )   :: nlon        ! # of lon midpoints
        1407             :     REAL(hp),         POINTER         :: LatEdge(:)  ! lat edges
        1408             :     CHARACTER(LEN=*), INTENT(IN   )   :: FN          ! filename
        1409             : !
        1410             : ! !INPUT/OUTPUT PARAMETERS:
        1411             : !
        1412             :     REAL(sp),         POINTER         :: Array(:,:,:,:) ! Data
        1413             :     INTEGER,          INTENT(INOUT)   :: RC             ! Return code
        1414             : !
        1415             : ! !REVISION HISTORY:
        1416             : !  13 Mar 2013 - C. Keller - Initial version
        1417             : !  See https://github.com/geoschem/hemco for complete history
        1418             : !EOP
        1419             : !------------------------------------------------------------------------------
        1420             : !BOC
        1421             : !
        1422             : ! !LOCAL VARIABLES:
        1423             : !
        1424             :     REAL(hp)              :: DLAT, AREA
        1425             :     INTEGER               :: NLAT, J
        1426             :     CHARACTER(LEN=255)    :: MSG, LOC
        1427             :
        1428             :     !=================================================================
        1429             :     ! NORNALIZE_AREA begins here!
        1430             :     !=================================================================
        1431             :
        1432             :     ! Initialize
        1433           0 :     LOC    = 'NORMALIZE_AREA (hcoio_util_mod.F90 )'
        1434             :
        1435             :     ! Check array size
        1436           0 :     NLAT = SIZE(LatEdge,1) - 1
        1437             :
        1438           0 :     IF ( SIZE(Array,1) /= nlon ) THEN
        1439           0 :        MSG = 'Array size does not agree with nlon: ' // TRIM(FN)
        1440           0 :        CALL HCO_ERROR( MSG, RC, THISLOC=LOC )
        1441           0 :        RETURN
        1442             :     ENDIF
        1443           0 :     IF ( SIZE(Array,2) /= NLAT ) THEN
        1444           0 :        MSG = 'Array size does not agree with nlat: ' // TRIM(FN)
        1445           0 :        CALL HCO_ERROR( MSG, RC, THISLOC=LOC )
        1446           0 :        RETURN
        1447             :     ENDIF
        1448             :
        1449             :     ! Loop over all latitudes
        1450           0 :     DO J = 1, NLAT
        1451             :        ! get grid box area in m2 for grid box with lower and upper latitude
        1452             :        ! llat/ulat:  Area = 2 * PI * Re^2 * DLAT / nlon,
        1453             :        ! where DLAT = abs( sin(ulat) - sin(llat) )
        1454           0 :        DLAT = ABS( SIN(LatEdge(J+1)*HcoState%Phys%PI_180)  &
        1455           0 :                    - SIN(LatEdge(J)*HcoState%Phys%PI_180) )
        1456             :        AREA = ( 2_hp * HcoState%Phys%PI * DLAT * HcoState%Phys%Re**2 ) &
        1457           0 :               / REAL(nlon,hp)
        1458             :
        1459             :        ! convert array data to m-2
        1460           0 :        ARRAY(:,J,:,:) = ARRAY(:,J,:,:) / AREA
        1461             :     ENDDO
        1462             :
        1463             :     ! Prompt a warning
        1464           0 :     WRITE(MSG,*) 'No area unit found in ' // TRIM(FN) // ' - convert to m-2!'
        1465           0 :     CALL HCO_WARNING ( HcoState%Config%Err, MSG, RC, THISLOC=LOC )
        1466             :
        1467             :     ! Leave w/ success
        1468           0 :     RC = HCO_SUCCESS
        1469             :
        1470             :   END SUBROUTINE Normalize_Area
        1471             : !EOC
        1472             : !------------------------------------------------------------------------------
        1473             : !                   Harmonized Emissions Component (HEMCO)                    !
        1474             : !------------------------------------------------------------------------------
        1475             : !BOP
        1476             : !
        1477             : ! !IROUTINE: SrcFile_Parse
        1478             : !
        1479             : ! !DESCRIPTION: Routine SrcFile\_Parse parses the source file name ('ncFile')
        1480             : ! of the provided list container Lct. In particular, it searches for tokens
        1481             : ! such as $ROOT, $YYYY, etc., within the file name and replaces those values
        1482             : ! with the intendend characters. The parsed file name is returned in string
        1483             : ! srcFile, while the original file name is retained in Lct.
        1484             : !\\
        1485             : !\\
        1486             : ! It now also checks if the file exists. If the file does not exist and the
        1487             : ! file name contains date tokens, it tries to adjust the file name to the
        1488             : ! closest available date in the past. The optional flag FUTURE can be used
        1489             : ! to denote that the next available file in the future shall be selected,
        1490             : ! even if there is a file that exactly matches the preferred date time. This
        1491             : ! is useful for interpolation between fields.
        1492             : !\\
        1493             : !\\
        1494             : ! !INTERFACE:
        1495             : !
        1496           0 :   SUBROUTINE SrcFile_Parse ( HcoState, Lct, srcFile, FOUND, RC, &
        1497             :                              Direction, Year )
        1498             : !
        1499             : ! !USES:
        1500             : !
        1501             :     USE HCO_TIDX_MOD,         ONLY : HCO_GetPrefTimeAttr
        1502             :     USE HCO_TIDX_MOD,         ONLY : tIDx_IsInRange
        1503             :     USE HCO_CLOCK_MOD,        ONLY : HcoClock_Get
        1504             :     USE HCO_CLOCK_MOD,        ONLY : Get_LastDayOfMonth
        1505             : !
        1506             : ! !INPUT PARAMETERS:
        1507             : !
        1508             :     TYPE(HCO_State),  POINTER                 :: HcoState   ! HEMCO state object
        1509             :     TYPE(ListCont),   POINTER                 :: Lct        ! HEMCO list
        1510             :     INTEGER,          INTENT(IN   ), OPTIONAL :: Direction  ! Look for file in
        1511             :                                                             ! future (+1) or
        1512             :                                                             ! past (-1)
        1513             :     INTEGER,          INTENT(IN   ), OPTIONAL :: Year       ! To use fixed year
        1514             : !
        1515             : ! !OUTPUT PARAMETERS:
        1516             : !
        1517             :     CHARACTER(LEN=*), INTENT(  OUT)           :: srcFile    ! output string
        1518             :     LOGICAL,          INTENT(  OUT)           :: FOUND      ! Does file exist?
        1519             : !
        1520             : ! !INPUT/OUTPUT PARAMETERS:
        1521             : !
        1522             :     INTEGER,          INTENT(INOUT)           :: RC         ! return code
        1523             : !
        1524             : ! !REVISION HISTORY:
        1525             : !  01 Oct 2014 - C. Keller - Initial version
        1526             : !  See https://github.com/geoschem/hemco for complete history
        1527             : !EOP
        1528             : !------------------------------------------------------------------------------
        1529             : !BOC
        1530             : !
        1531             : ! !LOCAL VARIABLES:
        1532             : !
        1533             :     INTEGER :: INC,     CNT,    TYPCNT, TYP,   NEWTYP
        1534             :     INTEGER :: prefYr,  prefMt, prefDy, prefHr, prefMn
        1535             :     INTEGER :: origYr,  origMt, origDy, origHr
        1536             :     LOGICAL :: hasFile, hasYr,  hasMt,  hasDy, hasHr
        1537             :     LOGICAL :: nextTyp
        1538             :     CHARACTER(LEN=1023) :: MSG, LOC
        1539             :     CHARACTER(LEN=1023) :: srcFileOrig
        1540             :
        1541             :     ! maximum # of iterations for file search
        1542             :     INTEGER, PARAMETER :: MAXIT = 10000
        1543             :
        1544             :     !=================================================================
        1545             :     ! SrcFile_Parse
        1546             :     !=================================================================
        1547             :
        1548             :     ! Initialize
        1549           0 :     LOC     = 'SrcFile_Parse (HCOIO_UTIL_MOD.F90)'
        1550           0 :     RC      = HCO_SUCCESS
        1551           0 :     found   = .FALSE.
        1552           0 :     srcFile = Lct%Dct%Dta%ncFile
        1553             :
        1554             :     ! verbose mode
        1555           0 :     IF ( HCO_IsVerb( HcoState%Config%Err ) ) THEN
        1556           0 :        WRITE(MSG,*) 'Parsing source file and replacing tokens'
        1557           0 :        CALL HCO_MSG(HcoState%Config%Err,MSG)
        1558             :     ENDIF
        1559             :
        1560             :     ! Get preferred dates (to be passed to parser)
        1561             :     CALL HCO_GetPrefTimeAttr ( HcoState, Lct, &
        1562           0 :                                prefYr, prefMt, prefDy, prefHr, prefMn, RC )
        1563           0 :     IF ( RC /= HCO_SUCCESS ) THEN
        1564           0 :         CALL HCO_ERROR( 'ERROR 1', RC, THISLOC=LOC )
        1565           0 :         RETURN
        1566             :     ENDIF
        1567             :
        1568             :     ! Make sure dates are not negative
        1569           0 :     IF ( prefYr <= 0 ) THEN
        1570           0 :        CALL HcoClock_Get( HcoState%Clock, cYYYY = prefYr, RC = RC )
        1571           0 :        IF ( RC /= HCO_SUCCESS ) THEN
        1572           0 :            CALL HCO_ERROR( 'ERROR 2', RC, THISLOC=LOC )
        1573           0 :            RETURN
        1574             :        ENDIF
        1575             :     ENDIF
        1576           0 :     IF ( prefMt <= 0 ) THEN
        1577           0 :        CALL HcoClock_Get( HcoState%Clock, cMM   = prefMt, RC = RC )
        1578           0 :        IF ( RC /= HCO_SUCCESS ) THEN
        1579           0 :            CALL HCO_ERROR( 'ERROR 3', RC, THISLOC=LOC )
        1580           0 :            RETURN
        1581             :        ENDIF
        1582             :     ENDIF
        1583           0 :     IF ( prefDy <= 0 ) THEN
        1584           0 :        CALL HcoClock_Get( HcoState%Clock, cDD   = prefDy, RC = RC )
        1585           0 :        IF ( RC /= HCO_SUCCESS ) THEN
        1586           0 :            CALL HCO_ERROR( 'ERROR 4', RC, THISLOC=LOC )
        1587           0 :            RETURN
        1588             :        ENDIF
        1589             :     ENDIF
        1590           0 :     IF ( prefHr <  0 ) THEN
        1591           0 :        CALL HcoClock_Get( HcoState%Clock, cH    = prefHr, RC = RC )
        1592           0 :        IF ( RC /= HCO_SUCCESS ) THEN
        1593           0 :            CALL HCO_ERROR( 'ERROR 5', RC, THISLOC=LOC )
        1594           0 :            RETURN
        1595             :        ENDIF
        1596             :     ENDIF
        1597             :
        1598             :     ! Eventually replace default preferred year with specified one
        1599           0 :     IF ( PRESENT(Year) ) prefYr = Year
        1600             :
        1601             :     ! Call the parser
        1602             :     CALL HCO_CharParse ( HcoState%Config, srcFile, prefYr, prefMt, &
        1603           0 :                          prefDy, prefHr, prefMn, RC )
        1604           0 :     IF ( RC /= HCO_SUCCESS ) THEN
        1605           0 :         CALL HCO_ERROR( 'ERROR 6', RC, THISLOC=LOC )
        1606           0 :         RETURN
        1607             :     ENDIF
        1608           0 :     srcFileOrig = TRIM(srcFile)
        1609             :
        1610             :     ! Check if file exists
        1611           0 :     INQUIRE( FILE=TRIM(srcFile), EXIST=HasFile )
        1612             :
        1613             :     ! If the direction flag is on, force HasFile to be false.
        1614           0 :     IF ( PRESENT(Direction) ) THEN
        1615           0 :        IF ( Direction /= 0 ) HasFile = .FALSE.
        1616             :     ENDIF
        1617             :
        1618             :     !-----------------------------------------------------------------------
        1619             :     ! If this is a HEMCO dry-run simulation, then do not enter the loop
        1620             :     ! where we will attempt to go back in time until a file is found.
        1621             :     ! For the dry-run we need to report all files, even missing.
        1622             :     ! This fixes Github issue geoschem/geos-chem #312. (bmy, 6/9/20)
        1623             :     !-----------------------------------------------------------------------
        1624           0 :     IF ( HcoState%Options%isDryRun ) THEN
        1625             :
        1626             :        ! Make sure that the year is not 1, this indicates that the
        1627             :        ! preferred year is outside of the years specified in the
        1628             :        ! time range settings in the configuration file, and will
        1629             :        ! lead to files with a year of "0001" in the path.
        1630             :        ! (bmy, 6/9/20)
        1631           0 :        IF ( prefyr == 1 ) THEN
        1632             :           MSG = 'Cannot find file for current simulation time: ' // &
        1633             :                TRIM(srcFile) // ' - Cannot get field ' // &
        1634             :                TRIM(Lct%Dct%cName) // '. Please check file name ' // &
        1635           0 :                'and time (incl. time range flag) in the config. file'
        1636           0 :           CALL HCO_ERROR( MSG, RC )
        1637           0 :           RETURN
        1638             :        ENDIF
        1639             :
        1640             :        ! Otherwise return with success
        1641           0 :        RC    = HCO_SUCCESS
        1642           0 :        Found = HasFile
        1643           0 :        RETURN
        1644             :     ENDIF
        1645             :
        1646             :     ! If file does not exist, check if we can adjust prefYr, prefMt, etc.
        1647           0 :     IF ( .NOT. HasFile .AND. Lct%Dct%DctType /= HCO_CFLAG_EXACT ) THEN
        1648             :
        1649             :        ! Check if any token exist
        1650           0 :        HasYr = ( INDEX(TRIM(Lct%Dct%Dta%ncFile),'YYYY') > 0 )
        1651           0 :        HasMt = ( INDEX(TRIM(Lct%Dct%Dta%ncFile),'MM'  ) > 0 )
        1652           0 :        HasDy = ( INDEX(TRIM(Lct%Dct%Dta%ncFile),'DD'  ) > 0 )
        1653           0 :        HasHr = ( INDEX(TRIM(Lct%Dct%Dta%ncFile),'HH'  ) > 0 )
        1654             :
        1655             :        ! Search for file
        1656           0 :        IF ( HasYr .OR. HasMt .OR. HasDy .OR. HasHr ) THEN
        1657             :
        1658             :           ! Date increments
        1659           0 :           INC = -1
        1660           0 :           IF ( PRESENT(Direction) ) THEN
        1661           0 :              INC = Direction
        1662             :           ENDIF
        1663             :
        1664             :           ! Initialize counters
        1665           0 :           CNT = 0
        1666             :
        1667             :           ! Type is the update type (see below)
        1668           0 :           TYP = 0
        1669             :
        1670             :           ! Mirror preferred variables
        1671           0 :           origYr = prefYr
        1672           0 :           origMt = prefMt
        1673           0 :           origDy = prefDy
        1674           0 :           origHr = prefHr
        1675             :
        1676             :           ! Do until file is found or counter exceeds threshold
        1677           0 :           DO WHILE ( .NOT. HasFile )
        1678             :
        1679             :              ! Inrease counter
        1680           0 :              CNT = CNT + 1
        1681           0 :              IF ( CNT > MAXIT ) EXIT
        1682             :
        1683             :              ! Increase update type if needed:
        1684           0 :              nextTyp = .FALSE.
        1685             :
        1686             :              ! Type 0: Initialization
        1687           0 :              IF ( TYP == 0 ) THEN
        1688             :                 nextTyp = .TRUE.
        1689             :              ! Type 1: update hour only
        1690           0 :              ELSEIF ( TYP == 1 .AND. TYPCNT > 24 ) THEN
        1691             :                 nextTyp = .TRUE.
        1692             :              ! Type 2: update day only
        1693           0 :              ELSEIF ( TYP == 2 .AND. TYPCNT > 31 ) THEN
        1694             :                 nextTyp = .TRUE.
        1695             :              ! Type 3: update month only
        1696           0 :              ELSEIF ( TYP == 3 .AND. TYPCNT > 12 ) THEN
        1697             :                 nextTyp = .TRUE.
        1698             :              ! Type 4: update year only
        1699           0 :              ELSEIF ( TYP == 4 .AND. TYPCNT > 300 ) THEN
        1700             :                 nextTyp = .TRUE.
        1701             :              ! Type 5: update hour and day
        1702           0 :              ELSEIF ( TYP == 5 .AND. TYPCNT > 744 ) THEN
        1703             :                 nextTyp = .TRUE.
        1704             :              ! Type 6: update day and month
        1705           0 :              ELSEIF ( TYP == 6 .AND. TYPCNT > 372 ) THEN
        1706             :                 nextTyp = .TRUE.
        1707             :              ! Type 7: update month and year
        1708           0 :              ELSEIF ( TYP == 7 .AND. TYPCNT > 3600 ) THEN
        1709             :                 EXIT
        1710             :              ENDIF
        1711             :
        1712             :              ! Get next type
        1713             :              IF ( nextTyp ) THEN
        1714           0 :                 NEWTYP = -1
        1715           0 :                 IF     ( hasHr .AND. TYP < 1 ) THEN
        1716             :                    NEWTYP = 1
        1717           0 :                 ELSEIF ( hasDy .AND. TYP < 2 ) THEN
        1718             :                    NEWTYP = 2
        1719           0 :                 ELSEIF ( hasMt .AND. TYP < 3 ) THEN
        1720             :                    NEWTYP = 3
        1721           0 :                 ELSEIF ( hasYr .AND. TYP < 4 ) THEN
        1722             :                    NEWTYP = 4
        1723             :                 ELSEIF ( hasDy .AND. TYP < 2 ) THEN
        1724             :                    NEWTYP = 5
        1725             :                 ELSEIF ( hasDy .AND. TYP < 2 ) THEN
        1726             :                    NEWTYP = 6
        1727             :                 ELSEIF ( hasDy .AND. TYP < 2 ) THEN
        1728             :                    NEWTYP = 7
        1729             :                 ENDIF
        1730             :
        1731             :                 ! Exit if no other type found
        1732             :                 IF ( NEWTYP < 0 ) EXIT
        1733             :
        1734             :                 ! This is the new type, reset type counter
        1735           0 :                 TYP    = NEWTYP
        1736           0 :                 TYPCNT = 0
        1737             :
        1738             :                 ! Make sure we reset all values
        1739           0 :                 prefYr = origYr
        1740           0 :                 prefMt = origMt
        1741           0 :                 prefDy = origDy
        1742           0 :                 prefHr = origHr
        1743             :
        1744             :              ENDIF
        1745             :
        1746             :              ! Update preferred datetimes
        1747           0 :              SELECT CASE ( TYP )
        1748             :                 ! Adjust hour only
        1749             :                 CASE ( 1 )
        1750           0 :                    prefHr = prefHr + INC
        1751             :                 ! Adjust day only
        1752             :                 CASE ( 2 )
        1753           0 :                    prefDy = prefDy + INC
        1754             :                 ! Adjust month only
        1755             :                 CASE ( 3 )
        1756           0 :                    prefMt = prefMt + INC
        1757             :                 ! Adjust year only
        1758             :                 CASE ( 4 )
        1759           0 :                    prefYr = prefYr + INC
        1760             :                 ! Adjust hour and day
        1761             :                 CASE ( 5 )
        1762           0 :                    prefHr = prefHr + INC
        1763           0 :                    IF ( MOD(TYPCNT,24) == 0 ) prefDy = prefDy + INC
        1764             :                 ! Adjust day and month
        1765             :                 CASE ( 6 )
        1766           0 :                    prefDy = prefDy + INC
        1767           0 :                    IF ( MOD(TYPCNT,31) == 0 ) prefMt = prefMt + INC
        1768             :                 ! Adjust month and year
        1769             :                 CASE ( 7 )
        1770           0 :                    prefMt = prefMt + INC
        1771           0 :                    IF ( MOD(TYPCNT,12) == 0 ) prefYr = prefYr + INC
        1772             :                 CASE DEFAULT
        1773           0 :                    EXIT
        1774             :              END SELECT
        1775             :
        1776             :              ! Check if we need to adjust a year/month/day/hour
        1777           0 :              IF ( prefHr < 0 ) THEN
        1778           0 :                 prefHr = 23
        1779           0 :                 prefDy = prefDy - 1
        1780             :              ENDIF
        1781           0 :              IF ( prefHr > 23 ) THEN
        1782           0 :                 prefHr = 0
        1783           0 :                 prefDy = prefDy + 1
        1784             :              ENDIF
        1785           0 :              IF ( prefDy < 1  ) THEN
        1786           0 :                 prefDy = 31
        1787           0 :                 prefMt = prefMt - 1
        1788             :              ENDIF
        1789           0 :              IF ( prefDy > 31 ) THEN
        1790           0 :                 prefDy = 1
        1791           0 :                 prefMt = prefMt + 1
        1792             :              ENDIF
        1793           0 :              IF ( prefMt < 1  ) THEN
        1794           0 :                 prefMt = 12
        1795           0 :                 prefYr = prefYr - 1
        1796             :              ENDIF
        1797           0 :              IF ( prefMt > 12 ) THEN
        1798           0 :                 prefMt = 1
        1799           0 :                 prefYr = prefYr + 1
        1800             :              ENDIF
        1801             :
        1802             :              ! Make sure day does not exceed max. number of days in this month
        1803           0 :              prefDy = MIN( prefDy, Get_LastDayOfMonth( prefMt, prefYr ) )
        1804             :
        1805             :              ! Mirror original file
        1806           0 :              srcFile = Lct%Dct%Dta%ncFile
        1807             :
        1808             :              ! Call the parser with adjusted values
        1809             :              CALL HCO_CharParse ( HcoState%Config, srcFile, prefYr, &
        1810           0 :                                   prefMt, prefDy, prefHr, prefMn, RC )
        1811           0 :              IF ( RC /= HCO_SUCCESS ) THEN
        1812           0 :                  CALL HCO_ERROR( 'ERROR 7', RC, THISLOC=LOC )
        1813           0 :                  RETURN
        1814             :              ENDIF
        1815             :
        1816             :              ! Check if this file exists
        1817           0 :              INQUIRE( FILE=TRIM(srcFile), EXIST=HasFile )
        1818             :
        1819             :              ! Update counter
        1820           0 :              TYPCNT = TYPCNT + 1
        1821             :           ENDDO
        1822             :        ENDIF
        1823             :     ENDIF
        1824             :
        1825             :     ! Additional check for data with a given range: make sure that the selected
        1826             :     ! field is not outside of the given range
        1827           0 :     IF ( HasFile .AND. ( Lct%Dct%Dta%CycleFlag == HCO_CFLAG_RANGE ) ) THEN
        1828           0 :        HasFile = TIDX_IsInRange ( Lct, prefYr, prefMt, prefDy, prefHr )
        1829             :     ENDIF
        1830             :
        1831             :     ! Restore original source file name and date to avoid confusion in log file
        1832           0 :     IF ( .not. HasFile ) THEN
        1833           0 :        srcFile = Trim(srcFileOrig)
        1834             :     ENDIF
        1835             :
        1836             :     ! Return variable
        1837           0 :     FOUND = HasFile
        1838             :
        1839             :     ! Return w/ success
        1840           0 :     RC = HCO_SUCCESS
        1841             :
        1842           0 :   END SUBROUTINE SrcFile_Parse
        1843             : !EOC
        1844             : !------------------------------------------------------------------------------
        1845             : !                   Harmonized Emissions Component (HEMCO)                    !
        1846             : !------------------------------------------------------------------------------
        1847             : !BOP
        1848             : !
        1849             : ! !IROUTINE: SigmaMidToEdges
        1850             : !
        1851             : ! !DESCRIPTION: Helper routine to interpolate sigma mid point values to edges.
        1852             : ! A simple linear interpolation is performed.
        1853             : !\\
        1854             : !\\
        1855             : ! !INTERFACE:
        1856             : !
        1857           0 :   SUBROUTINE SigmaMidToEdges ( HcoState, SigMid, SigEdge, RC )
        1858             : !
        1859             : ! !INPUT PARAMETERS:
        1860             : !
        1861             :     TYPE(HCO_State),  POINTER                 :: HcoState        ! HEMCO state
        1862             :     REAL(hp),         POINTER                 :: SigMid(:,:,:)   ! sigma levels
        1863             : !
        1864             : ! !OUTPUT PARAMETERS:
        1865             : !
        1866             :     REAL(hp),         POINTER                 :: SigEdge(:,:,:)  ! sigma edges
        1867             :     INTEGER,          INTENT(  OUT)           :: RC              ! return code
        1868             : !
        1869             : ! !REVISION HISTORY:
        1870             : !  03 Oct 2013 - C. Keller - Initial version
        1871             : !  See https://github.com/geoschem/hemco for complete history
        1872             : !EOP
        1873             : !------------------------------------------------------------------------------
        1874             : !BOC
        1875             : !
        1876             : ! !LOCAL VARIABLES:
        1877             : !
        1878             :     INTEGER            :: L, AS
        1879             :     INTEGER            :: nx, ny, nz
        1880             :     CHARACTER(LEN=255) :: MSG
        1881             :     CHARACTER(LEN=255) :: LOC = 'SigmaMidToEdges (hcoio_util_mod.F90)'
        1882             :
        1883             :     !=================================================================
        1884             :     ! SigmaMidToEdges begins here!
        1885             :     !=================================================================
        1886             :
        1887             :     ! Allocate space as required
        1888           0 :     nx = SIZE(SigMid,1)
        1889           0 :     ny = SIZE(SigMid,2)
        1890           0 :     nz = SIZE(SigMid,3)
        1891           0 :     IF ( ASSOCIATED(SigEdge) ) DEALLOCATE(SigEdge)
        1892           0 :     ALLOCATE(SigEdge(nx,ny,nz+1),STAT=AS)
        1893           0 :     IF ( AS/=0 ) THEN
        1894             :        CALL HCO_ERROR( 'Allocate SigEdge', RC, &
        1895           0 :                        THISLOC=LOC )
        1896           0 :        RETURN
        1897             :     ENDIF
        1898           0 :     SigEdge = 0.0_hp
        1899             :
        1900             :     ! Calculate sigma edges by linear interpolation (symmetric mid-points)
        1901           0 :     DO L = 1, nz-1
        1902           0 :        SigEdge(:,:,L+1) = ( SigMid(:,:,L) + SigMid(:,:,L+1) ) / 2.0_hp
        1903             :     ENDDO
        1904             :
        1905             :     ! Get outermost values:
        1906           0 :     SigEdge(:,:,1   ) = SigMid(:,:,1 ) - ( SigEdge(:,:,2) - SigMid(:,:,1)   )
        1907           0 :     SigEdge(:,:,nz+1) = SigMid(:,:,nz) + ( SigMid(:,:,nz) - SigEdge(:,:,nz) )
        1908             :
        1909             :     ! Return w/ success
        1910           0 :     RC = HCO_SUCCESS
        1911             :
        1912             :   END SUBROUTINE SigmaMidToEdges
        1913             : !EOC
        1914             : !------------------------------------------------------------------------------
        1915             : !                   Harmonized Emissions Component (HEMCO)                    !
        1916             : !------------------------------------------------------------------------------
        1917             : !BOP
        1918             : !
        1919             : ! !IROUTINE: CheckMissVal
        1920             : !
        1921             : ! !DESCRIPTION: Checks for missing values in the passed array. Missing values
        1922             : ! of base emissions and masks are set to 0, missing values of scale factors
        1923             : ! are set to 1.
        1924             : !\\
        1925             : ! !INTERFACE:
        1926             : !
        1927           0 :   SUBROUTINE CheckMissVal ( Lct, Arr )
        1928             : !
        1929             : ! !INPUT PARAMETERS:
        1930             : !
        1931             :     TYPE(ListCont),   POINTER                 :: Lct
        1932             :     REAL(sp),         POINTER                 :: Arr(:,:,:,:)
        1933             : !
        1934             : ! !REVISION HISTORY:
        1935             : !  04 Mar 2015 - C. Keller - Initial version
        1936             : !  See https://github.com/geoschem/hemco for complete history
        1937             : !EOP
        1938             : !------------------------------------------------------------------------------
        1939             : !BOC
        1940             : !
        1941             : ! !LOCAL VARIABLES:
        1942             : !
        1943             :     !=================================================================
        1944             :     ! CheckMissVal begins here!
        1945             :     !=================================================================
        1946             :
        1947             :     ! Error trap
        1948           0 :     IF ( .NOT. ASSOCIATED(Arr) ) RETURN
        1949             :
        1950           0 :     IF ( ANY(Arr == HCO_MISSVAL) ) THEN
        1951             :        ! Base emissions
        1952           0 :        IF ( Lct%Dct%DctType == HCO_DCTTYPE_BASE ) THEN
        1953           0 :           WHERE(Arr == HCO_MISSVAL) Arr = 0.0_sp
        1954             :        ! Scale factor
        1955           0 :        ELSEIF ( Lct%Dct%DctType == HCO_DCTTYPE_SCAL ) THEN
        1956           0 :           WHERE(Arr == HCO_MISSVAL) Arr = 1.0_sp
        1957             :        ! Mask
        1958           0 :        ELSEIF ( Lct%Dct%DctType == HCO_DCTTYPE_MASK ) THEN
        1959           0 :           WHERE(Arr == HCO_MISSVAL) Arr = 0.0_sp
        1960             :        ENDIF
        1961             :     ENDIF
        1962             :
        1963             :   END SUBROUTINE CheckMissVal
        1964             : !EOC
        1965             : !------------------------------------------------------------------------------
        1966             : !                   Harmonized Emissions Component (HEMCO)                    !
        1967             : !------------------------------------------------------------------------------
        1968             : !BOP
        1969             : !
        1970             : ! !IROUTINE: GetArbDimIndex
        1971             : !
        1972             : ! !DESCRIPTION: Subroutine GetArbDimIndex returns the index of the arbitrary
        1973             : ! file dimension. -1 if no such dimension is defined.
        1974             : !\\
        1975             : ! !INTERFACE:
        1976             : !
        1977           0 :   SUBROUTINE GetArbDimIndex( HcoState, Lun, Lct, ArbIdx, RC )
        1978             : !
        1979             : ! !USES:
        1980             : !
        1981             :     USE HCO_m_netcdf_io_checks
        1982             :     USE HCO_m_netcdf_io_get_dimlen
        1983             :     USE HCO_ExtList_Mod,    ONLY : GetExtOpt
        1984             : !
        1985             : ! !INPUT PARAMETERS:
        1986             : !
        1987             :     TYPE(HCO_State),  POINTER                 :: HcoState
        1988             :     INTEGER,          INTENT(IN   )           :: Lun
        1989             :     TYPE(ListCont),   POINTER                 :: Lct
        1990             : !
        1991             : ! !OUTPUT PARAMETERS:
        1992             : !
        1993             :     INTEGER,          INTENT(  OUT)           :: ArbIdx
        1994             :     INTEGER,          INTENT(  OUT)           :: RC
        1995             : !
        1996             : ! !REVISION HISTORY:
        1997             : !  22 Sep 2015 - C. Keller - Initial version
        1998             : !  See https://github.com/geoschem/hemco for complete history
        1999             : !EOP
        2000             : !------------------------------------------------------------------------------
        2001             : !BOC
        2002             : !
        2003             : ! !LOCAL VARIABLES:
        2004             : !
        2005             :     INTEGER             :: TargetVal, nVal
        2006             :     LOGICAL             :: Found
        2007             :     CHARACTER(LEN=255)  :: ArbDimVal
        2008             :     CHARACTER(LEN=511)  :: MSG
        2009             :     CHARACTER(LEN=255)  :: LOC = 'GetArbDimIndex (hcoio_util_mod.F90)'
        2010             :
        2011             :     !=================================================================
        2012             :     ! GetArbDimIndex
        2013             :     !=================================================================
        2014             :
        2015             :     ! Assume success until otherwise
        2016           0 :     RC = HCO_SUCCESS
        2017             :
        2018             :     ! Init
        2019           0 :     ArbIdx = -1
        2020           0 :     IF ( TRIM(Lct%Dct%Dta%ArbDimName) == 'none' ) RETURN
        2021             :
        2022             :     ! Check if variable exists
        2023           0 :     Found = Ncdoes_Dim_Exist ( Lun, TRIM(Lct%Dct%Dta%ArbDimName) )
        2024           0 :     IF ( .NOT. Found ) THEN
        2025             :        MSG = 'Cannot read dimension ' // TRIM(Lct%Dct%Dta%ArbDimName) &
        2026             :              // ' from file ' // &
        2027           0 :              TRIM(Lct%Dct%Dta%ncFile)
        2028           0 :        CALL HCO_ERROR( MSG, RC, THISLOC=LOC )
        2029           0 :        RETURN
        2030             :     ENDIF
        2031             :
        2032             :     ! Get dimension length
        2033           0 :     CALL Ncget_Dimlen ( Lun, TRIM(Lct%Dct%Dta%ArbDimName), nVal )
        2034             :
        2035             :     ! Get value to look for. This is archived in variable ArbDimVal.
        2036             :     ! Eventually need to extract value from HEMCO settings
        2037           0 :     ArbDimVal = TRIM(Lct%Dct%Dta%ArbDimVal)
        2038             :
        2039             :     ! If string starts with a number, evaluate value directly
        2040             :     IF ( ArbDimVal(1:1) == '0' .OR. &
        2041             :          ArbDimVal(1:1) == '1' .OR. &
        2042             :          ArbDimVal(1:1) == '2' .OR. &
        2043             :          ArbDimVal(1:1) == '3' .OR. &
        2044             :          ArbDimVal(1:1) == '4' .OR. &
        2045             :          ArbDimVal(1:1) == '5' .OR. &
        2046             :          ArbDimVal(1:1) == '6' .OR. &
        2047             :          ArbDimVal(1:1) == '7' .OR. &
        2048           0 :          ArbDimVal(1:1) == '8' .OR. &
        2049             :          ArbDimVal(1:1) == '9'       ) THEN
        2050           0 :        READ(ArbDimVal,*) TargetVal
        2051             :
        2052             :     ! Otherwise, assume this is a HEMCO option (including a token)
        2053             :     ELSE
        2054           0 :        IF ( ArbDimVal(1:1) == '$' ) ArbDimVal = ArbDimVal(2:LEN(ArbDimVal))
        2055             :        CALL GetExtOpt ( HcoState%Config, ExtNr=-999, &
        2056             :                         OptName=TRIM(ArbDimVal), &
        2057           0 :                         OptValInt=TargetVal, FOUND=Found, RC=RC )
        2058           0 :        IF ( RC /= HCO_SUCCESS ) THEN
        2059           0 :            CALL HCO_ERROR( 'ERROR 8', RC, THISLOC=LOC )
        2060           0 :            RETURN
        2061             :        ENDIF
        2062           0 :        IF ( .NOT. Found ) THEN
        2063           0 :           WRITE(MSG,*) 'Cannot evaluate additional dimension value ', &
        2064           0 :              TRIM(ArbDimVal), '. This does not seem to be a number nor ', &
        2065           0 :              'a HEMCO token/setting. This error happened when evaluating ', &
        2066           0 :              'dimension ', TRIM(Lct%Dct%Dta%ArbDimName), ' belonging to ', &
        2067           0 :              'file ', TRIM(Lct%Dct%Dta%ncFile)
        2068           0 :           CALL HCO_ERROR( MSG, RC, THISLOC=LOC )
        2069           0 :           RETURN
        2070             :        ENDIF
        2071             :     ENDIF
        2072             :
        2073           0 :     IF ( TargetVal > nVal ) THEN
        2074           0 :        WRITE(MSG,*) 'Desired dimension value ', TargetVal, &
        2075           0 :           ' exceeds corresponding dimension length on that file: ', nVal, &
        2076           0 :           'This error happened when evaluating ', &
        2077           0 :           'dimension ', TRIM(Lct%Dct%Dta%ArbDimName), ' belonging to ', &
        2078           0 :           'file ', TRIM(Lct%Dct%Dta%ncFile)
        2079           0 :        CALL HCO_ERROR( MSG, RC, THISLOC=LOC )
        2080           0 :        RETURN
        2081             :
        2082             :     ELSE
        2083           0 :        ArbIdx = TargetVal
        2084             :     ENDIF
        2085             :
        2086             :     ! Verbose
        2087           0 :     IF ( HcoState%amIRoot .AND. HCO_IsVerb( HcoState%Config%Err ) ) THEN
        2088           0 :        WRITE(MSG,*) 'Additional dimension ', TRIM(Lct%Dct%Dta%ArbDimName), &
        2089           0 :                     ' in ', TRIM(Lct%Dct%Dta%ncFile), ': use index ',      &
        2090           0 :                     ArbIdx, ' (set: ', Lct%Dct%Dta%ArbDimVal, ')'
        2091           0 :        CALL HCO_MSG(HcoState%Config%Err,MSG)
        2092             :     ENDIF
        2093             :
        2094             :     ! Return w/ success
        2095           0 :     RC = HCO_SUCCESS
        2096             :
        2097             :   END SUBROUTINE GetArbDimIndex
        2098             : !EOC
        2099             : #endif
        2100             : !------------------------------------------------------------------------------
        2101             : !                   Harmonized Emissions Component (HEMCO)                    !
        2102             : !------------------------------------------------------------------------------
        2103             : !BOP
        2104             : !
        2105             : ! !IROUTINE: HCOIO_ReadOther
        2106             : !
        2107             : ! !DESCRIPTION: Subroutine HCOIO\_ReadOther is a wrapper routine to
        2108             : ! read data from sources other than netCDF.
        2109             : !\\
        2110             : !\\
        2111             : ! If a file name is given (ending with '.txt'), the data are assumed
        2112             : ! to hold country-specific values (e.g. diurnal scale factors). In all
        2113             : ! other cases, the data is directly read from the configuration file
        2114             : ! (scalars).
        2115             : !\\
        2116             : !\\
        2117             : ! !INTERFACE:
        2118             : !
        2119           0 :   SUBROUTINE HCOIO_ReadOther( HcoState, Lct, RC )
        2120             : !
        2121             : ! !USES:
        2122             : !
        2123             : !
        2124             : ! !INPUT PARAMTERS:
        2125             : !
        2126             :     TYPE(HCO_State), POINTER          :: HcoState    ! HEMCO state
        2127             : !
        2128             : ! !INPUT/OUTPUT PARAMETERS:
        2129             : !
        2130             :     TYPE(ListCont),   POINTER         :: Lct
        2131             :     INTEGER,          INTENT(INOUT)   :: RC
        2132             : !
        2133             : ! !REVISION HISTORY:
        2134             : !  22 Dec 2014 - C. Keller: Initial version
        2135             : !  See https://github.com/geoschem/hemco for complete history
        2136             : !EOP
        2137             : !------------------------------------------------------------------------------
        2138             : !BOC
        2139             : !
        2140             : ! !LOCAL VARIABLES:
        2141             : !
        2142             :     CHARACTER(LEN=255) :: MSG, LOC
        2143             :
        2144             :     !======================================================================
        2145             :     ! HCOIO_ReadOther begins here
        2146             :     !======================================================================
        2147           0 :     LOC = 'HCOIO_ReadOther (HCOIO_UTIL_MOD.F90)'
        2148             :
        2149             :     ! Error check: data must be in local time
        2150           0 :     IF ( .NOT. Lct%Dct%Dta%IsLocTime ) THEN
        2151             :        MSG = 'Cannot read data from file that is not in local time: ' // &
        2152           0 :              TRIM(Lct%Dct%cName)
        2153           0 :        CALL HCO_ERROR( MSG, RC, THISLOC='HCOIO_ReadOther (hcoio_dataread_mod.F90)' )
        2154           0 :        RETURN
        2155             :     ENDIF
        2156             :
        2157             :     ! Read an ASCII file as country values
        2158           0 :     IF ( INDEX( TRIM(Lct%Dct%Dta%ncFile), '.txt' ) > 0 ) THEN
        2159           0 :        CALL HCOIO_ReadCountryValues( HcoState, Lct, RC )
        2160           0 :        IF ( RC /= HCO_SUCCESS ) THEN
        2161           0 :            CALL HCO_ERROR( 'ERROR 9', RC, THISLOC=LOC )
        2162           0 :            RETURN
        2163             :        ENDIF
        2164             :
        2165             :     ! Directly read from configuration file otherwise
        2166             :     ELSE
        2167           0 :        CALL HCOIO_ReadFromConfig( HcoState, Lct, RC )
        2168           0 :        IF ( RC /= HCO_SUCCESS ) THEN
        2169           0 :            CALL HCO_ERROR( 'ERROR 10', RC, THISLOC=LOC )
        2170           0 :            RETURN
        2171             :        ENDIF
        2172             :     ENDIF
        2173             :
        2174             :     ! Return w/ success
        2175           0 :     RC = HCO_SUCCESS
        2176             :
        2177             :   END SUBROUTINE HCOIO_ReadOther
        2178             : !EOC
        2179             : !------------------------------------------------------------------------------
        2180             : !                   Harmonized Emissions Component (HEMCO)                    !
        2181             : !------------------------------------------------------------------------------
        2182             : !BOP
        2183             : !
        2184             : ! !IROUTINE: HCOIO_ReadCountryValues
        2185             : !
        2186             : ! !DESCRIPTION: Subroutine HCOIO\_ReadCountryValues
        2187             : !\\
        2188             : !\\
        2189             : ! !INTERFACE:
        2190             : !
        2191           0 :   SUBROUTINE HCOIO_ReadCountryValues ( HcoState, Lct, RC )
        2192             : !
        2193             : ! !USES:
        2194             : !
        2195             :     USE HCO_inquireMod,     ONLY : findFreeLUN
        2196             :     USE HCO_CHARTOOLS_MOD,  ONLY : HCO_CMT, HCO_SPC, NextCharPos
        2197             :     USE HCO_EmisList_Mod,   ONLY : HCO_GetPtr
        2198             :     USE HCO_FileData_Mod,   ONLY : FileData_ArrCheck
        2199             : !
        2200             : ! !INPUT PARAMTERS:
        2201             : !
        2202             :     TYPE(HCO_State), POINTER          :: HcoState    ! HEMCO state
        2203             : !
        2204             : ! !INPUT/OUTPUT PARAMETERS:
        2205             : !
        2206             :     TYPE(ListCont),   POINTER         :: Lct
        2207             :     INTEGER,          INTENT(INOUT)   :: RC
        2208             : !
        2209             : ! !REVISION HISTORY:
        2210             : !  22 Dec 2014 - C. Keller: Initial version
        2211             : !  See https://github.com/geoschem/hemco for complete history
        2212             : !EOP
        2213             : !------------------------------------------------------------------------------
        2214             : !BOC
        2215             : !
        2216             : ! !LOCAL VARIABLES:
        2217             : !
        2218             :     INTEGER               :: IUFILE, IOS
        2219             :     INTEGER               :: ID1, ID2, I, NT, CID, NLINE
        2220           0 :     REAL(sp), POINTER     :: CNTR(:,:)
        2221           0 :     INTEGER,  ALLOCATABLE :: CIDS(:,:)
        2222           0 :     REAL(hp), POINTER     :: Vals(:)
        2223             :     LOGICAL               :: Verb
        2224             :     CHARACTER(LEN=2047)   :: LINE
        2225             :     CHARACTER(LEN=255)    :: MSG, DUM, CNT
        2226             :     CHARACTER(LEN=255)    :: LOC = 'HCOIO_ReadCountryValues (hcoio_util_mod.F90)'
        2227             :
        2228             :     !======================================================================
        2229             :     ! HCOIO_ReadCountryValues begins here
        2230             :     !======================================================================
        2231             :
        2232             :     ! Init
        2233           0 :     CNTR => NULL()
        2234           0 :     Vals => NULL()
        2235             :
        2236             :     ! verbose mode?
        2237           0 :     Verb = HCO_IsVerb( HcoState%Config%Err )
        2238             :
        2239             :     ! Verbose
        2240           0 :     IF ( Verb ) THEN
        2241           0 :        MSG = 'Use country-specific values for ' // TRIM(Lct%Dct%cName)
        2242           0 :        CALL HCO_MSG(HcoState%Config%Err,MSG)
        2243           0 :        MSG = '- Source file: ' // TRIM(Lct%Dct%Dta%ncFile)
        2244           0 :        CALL HCO_MSG(HcoState%Config%Err,MSG)
        2245             :     ENDIF
        2246             :
        2247             :     ! Open file
        2248           0 :     IUFILE = FindFreeLun()
        2249           0 :     OPEN ( IUFILE, FILE=TRIM( Lct%Dct%Dta%ncFile ), STATUS='OLD', IOSTAT=IOS )
        2250           0 :     IF ( IOS /= 0 ) THEN
        2251           0 :        MSG = 'Cannot open ' // TRIM(Lct%Dct%Dta%ncFile)
        2252           0 :        CALL HCO_ERROR( MSG, RC, THISLOC=LOC )
        2253           0 :        RETURN
        2254             :     ENDIF
        2255             :
        2256             :     ! Repeat for every line
        2257             :     NLINE = 0
        2258             :     DO
        2259             :
        2260             :        ! Read line
        2261           0 :        READ( IUFILE, '(a)', IOSTAT=IOS ) LINE
        2262             :
        2263             :        ! End of file?
        2264           0 :        IF ( IOS < 0 ) EXIT
        2265             :
        2266             :        ! Error?
        2267           0 :        IF ( IOS > 0 ) THEN
        2268           0 :           MSG = 'Error reading ' // TRIM(Lct%Dct%Dta%ncFile)
        2269           0 :           MSG = TRIM(MSG) // ' - last valid line: ' // TRIM(LINE)
        2270           0 :           CALL HCO_ERROR( MSG, RC, THISLOC=LOC )
        2271           0 :           RETURN
        2272             :        ENDIF
        2273             :
        2274             :        ! Skip commented lines and/or empty lines
        2275           0 :        IF ( TRIM(LINE) == ''      ) CYCLE
        2276           0 :        IF ( LINE(1:1)  == HCO_CMT ) CYCLE
        2277             :
        2278             :        ! First (valid) line holds the name of the mask container
        2279           0 :        IF ( NLINE == 0 ) THEN
        2280             :
        2281             :           ! Get pointer to mask. Convert to integer
        2282           0 :           CALL HCO_GetPtr( HcoState, TRIM(LINE), CNTR, RC )
        2283           0 :           IF ( RC /= HCO_SUCCESS ) THEN
        2284           0 :               CALL HCO_ERROR( 'ERROR 11', RC, THISLOC=LOC )
        2285           0 :               RETURN
        2286             :           ENDIF
        2287           0 :           ALLOCATE( CIDS(HcoState%NX, HcoState%NY), STAT=IOS )
        2288           0 :           IF ( IOS /= 0 ) THEN
        2289           0 :              CALL HCO_ERROR( 'Cannot allocate CIDS', RC, THISLOC=LOC )
        2290           0 :              RETURN
        2291             :           ENDIF
        2292           0 :           CIDS = NINT(CNTR)
        2293             :
        2294             :           ! Verbose
        2295           0 :           IF ( HCO_IsVerb( HcoState%Config%Err ) ) THEN
        2296           0 :              MSG = '- Use ID mask ' // TRIM(LINE)
        2297           0 :              CALL HCO_MSG(HcoState%Config%Err,MSG)
        2298             :           ENDIF
        2299             :
        2300             :           ! Go to next line
        2301           0 :           NLINE = NLINE + 1
        2302           0 :           CYCLE
        2303             :        ENDIF
        2304             :
        2305             :        ! Get first space character to skip country name.
        2306             :        ! We assume here that a country name is given right at the
        2307             :        ! beginning of the line, e.g. 'USA 744 1.05/1.02/...'
        2308           0 :        ID1 = NextCharPos( LINE, HCO_SPC )
        2309           0 :        CNT = LINE(1:ID1)
        2310             :
        2311             :        ! Get country ID
        2312           0 :        DO I = ID1, LEN(LINE)
        2313           0 :           IF ( LINE(I:I) /= HCO_SPC ) EXIT
        2314             :        ENDDO
        2315           0 :        ID1 = I
        2316           0 :        ID2 = NextCharPos( LINE, HCO_SPC, START=ID1 )
        2317             :
        2318           0 :        IF ( ID2 >= LEN(LINE) .OR. ID2 < 0 ) THEN
        2319           0 :           MSG = 'Cannot extract country ID from: ' // TRIM(LINE)
        2320           0 :           CALL HCO_ERROR( MSG, RC, THISLOC=LOC )
        2321           0 :           RETURN
        2322             :        ENDIF
        2323           0 :        DUM = LINE(ID1:ID2)
        2324           0 :        READ( DUM, * ) CID
        2325             :
        2326             :        ! Extract data values
        2327           0 :        ID1  = ID2+1
        2328           0 :        ID2  = LEN(LINE)
        2329           0 :        LINE = LINE(ID1:ID2)
        2330           0 :        CALL GetDataVals( HcoState, Lct, LINE, Vals, RC )
        2331           0 :        IF ( RC /= HCO_SUCCESS ) THEN
        2332           0 :            CALL HCO_ERROR( 'ERROR 12', RC, THISLOC=LOC )
        2333           0 :            RETURN
        2334             :        ENDIF
        2335             :
        2336             :        ! Check data / array dimensions
        2337           0 :        NT = SIZE(Vals,1)
        2338             :        CALL FileData_ArrCheck( HcoState%Config, Lct%Dct%Dta, &
        2339           0 :                                HcoState%NX, HcoState%NY, NT, RC )
        2340           0 :        IF ( RC /= HCO_SUCCESS ) THEN
        2341           0 :            CALL HCO_ERROR( 'ERROR 13', RC, THISLOC=LOC )
        2342           0 :            RETURN
        2343             :        ENDIF
        2344             :
        2345             :        ! Pass to data array. If the country ID is larger than zero, fill
        2346             :        ! only those grid boxes. Otherwise, fill all grid boxes that have
        2347             :        ! not yet been filled.
        2348           0 :        DO I = 1, NT
        2349           0 :           IF ( CID == 0 ) THEN
        2350           0 :              WHERE ( Lct%Dct%Dta%V2(I)%Val <= 0.0_sp )
        2351           0 :                 Lct%Dct%Dta%V2(I)%Val = Vals(I)
        2352             :              ENDWHERE
        2353             :           ELSE
        2354           0 :              WHERE ( CIDS == CID )
        2355           0 :                 Lct%Dct%Dta%V2(I)%Val = Vals(I)
        2356             :              ENDWHERE
        2357             :           ENDIF
        2358             :        ENDDO
        2359             :
        2360             :        ! Verbose
        2361           0 :        IF ( HCO_IsVerb( HcoState%Config%Err ) ) THEN
        2362           0 :           WRITE(MSG,*) '- Obtained values for ',TRIM(CNT),' ==> ID:', CID
        2363           0 :           CALL HCO_MSG(HcoState%Config%Err,MSG)
        2364             :        ENDIF
        2365             :
        2366             :        ! Cleanup
        2367           0 :        IF ( ASSOCIATED(Vals) ) DEALLOCATE( Vals )
        2368           0 :        Vals => NULL()
        2369             :
        2370             :        ! Update # of read lines
        2371           0 :        NLINE = NLINE + 1
        2372             :     ENDDO
        2373             :
        2374             :     ! Close file
        2375           0 :     CLOSE ( IUFILE )
        2376             :
        2377             :     ! Data is 2D
        2378           0 :     Lct%Dct%Dta%SpaceDim  = 2
        2379             :
        2380             :     ! Make sure data is in local time
        2381           0 :     IF ( .NOT. Lct%Dct%Dta%IsLocTime ) THEN
        2382           0 :        Lct%Dct%Dta%IsLocTime = .TRUE.
        2383             :        MSG = 'Data assigned to mask regions will be treated in local time: '//&
        2384           0 :               TRIM(Lct%Dct%cName)
        2385           0 :        CALL HCO_WARNING( HcoState%Config%Err, MSG, RC, THISLOC=LOC )
        2386             :     ENDIF
        2387             :
        2388             :     ! Cleanup
        2389           0 :     Cntr => NULL()
        2390           0 :     IF ( ALLOCATED(CIDS) ) DEALLOCATE ( CIDS )
        2391             :
        2392             :     ! Return w/ success
        2393           0 :     RC = HCO_SUCCESS
        2394             :
        2395           0 :   END SUBROUTINE HCOIO_ReadCountryValues
        2396             : !EOC
        2397             : !------------------------------------------------------------------------------
        2398             : !                   Harmonized Emissions Component (HEMCO)                    !
        2399             : !------------------------------------------------------------------------------
        2400             : !BOP
        2401             : !
        2402             : ! !IROUTINE: HCOIO_ReadFromConfig
        2403             : !
        2404             : ! !DESCRIPTION: Subroutine HCOIO\_ReadFromConfig reads data directly from
        2405             : ! the configuration file (instead of reading it from a netCDF file).
        2406             : ! These data is always assumed to be spatially uniform, but it is possible
        2407             : ! to specify multiple time slices by separating the individual time slice
        2408             : ! values by the HEMCO separator sign ('/' by default). The time dimension
        2409             : ! of these data is either determined from the srcTime attribute or estimated
        2410             : ! from the number of time slices provided. For example, if no srcTime is
        2411             : ! specified and 24 time slices are provided, data is assumed to represent
        2412             : ! hourly data. Similarly, data is assumed to represent weekdaily or monthly
        2413             : ! data for 7 or 12 time slices, respectively.
        2414             : !\\
        2415             : !\\
        2416             : ! If the srcTime attribute is defined, the time slices are determined from
        2417             : ! this attribute. Only one time dimension (year, month, day, or hour) can
        2418             : ! be defined for scalar fields!
        2419             : !\\
        2420             : !\\
        2421             : ! !INTERFACE:
        2422             : !
        2423           0 :   SUBROUTINE HCOIO_ReadFromConfig( HcoState, Lct, RC )
        2424             : !
        2425             : ! !USES:
        2426             : !
        2427             :     USE HCO_FILEDATA_MOD,   ONLY : FileData_ArrCheck
        2428             : !
        2429             : ! !INPUT PARAMTERS:
        2430             : !
        2431             :     TYPE(HCO_State), POINTER          :: HcoState    ! HEMCO state
        2432             : !
        2433             : ! !INPUT/OUTPUT PARAMETERS:
        2434             : !
        2435             :     TYPE(ListCont),   POINTER         :: Lct
        2436             :     INTEGER,          INTENT(INOUT)   :: RC
        2437             : !
        2438             : ! !REVISION HISTORY:
        2439             : !  24 Jul 2014 - C. Keller: Initial version
        2440             : !  See https://github.com/geoschem/hemco for complete history
        2441             : !EOP
        2442             : !------------------------------------------------------------------------------
        2443             : !BOC
        2444             : !
        2445             : ! !LOCAL VARIABLES:
        2446             : !
        2447             :     INTEGER            :: I, NT
        2448           0 :     REAL(hp), POINTER  :: Vals(:)
        2449             :     CHARACTER(LEN=255) :: MSG
        2450             :     CHARACTER(LEN=255) :: LOC = 'HCOIO_ReadFromConfig (hcoio_util_mod.F90)'
        2451             :
        2452             :     !======================================================================
        2453             :     ! HCOIO_ReadFromConfig begins here
        2454             :     !======================================================================
        2455             :
        2456             :     ! Init
        2457           0 :     Vals => NULL()
        2458             :
        2459             :     ! Verbose
        2460           0 :     IF ( HCO_IsVerb( HcoState%Config%Err ) ) THEN
        2461           0 :        WRITE(MSG, *) 'Read from config file: ', TRIM(Lct%Dct%cName)
        2462           0 :        CALL HCO_MSG(HcoState%Config%Err,MSG)
        2463             :     ENDIF
        2464             :
        2465             :     !-------------------------------------------------------------------
        2466             :     ! Get data values for this time step.
        2467             :     !-------------------------------------------------------------------
        2468           0 :     CALL GetDataVals( HcoState, Lct, Lct%Dct%Dta%ncFile, Vals, RC )
        2469           0 :     IF ( RC /= HCO_SUCCESS ) THEN
        2470           0 :         CALL HCO_ERROR( 'ERROR 14', RC, THISLOC=LOC )
        2471           0 :         RETURN
        2472             :     ENDIF
        2473             :
        2474             :     !-------------------------------------------------------------------
        2475             :     ! Copy data into array.
        2476             :     !-------------------------------------------------------------------
        2477             :
        2478             :     ! Number of values
        2479           0 :     NT = SIZE(Vals,1)
        2480             :
        2481             :     ! For masks, interpret data as mask corners (lon1/lat1/lon2/lat2)
        2482             :     ! with no time dimension
        2483           0 :     IF ( Lct%Dct%DctType == HCO_DCTTYPE_MASK ) THEN
        2484             :
        2485             :        ! Make sure data is allocated
        2486             :        CALL FileData_ArrCheck( HcoState%Config, Lct%Dct%Dta, &
        2487           0 :                                HcoState%NX, HcoState%NY, 1, RC )
        2488           0 :        IF ( RC /= HCO_SUCCESS ) THEN
        2489           0 :            CALL HCO_ERROR( 'ERROR 15', RC, THISLOC=LOC )
        2490           0 :            RETURN
        2491             :        ENDIF
        2492             :
        2493             :        ! Fill array: 1.0 within grid box, 0.0 outside.
        2494           0 :        CALL FillMaskBox( HcoState, Lct, Vals, RC )
        2495           0 :        IF ( RC /= HCO_SUCCESS ) THEN
        2496           0 :            CALL HCO_ERROR( 'ERROR 16', RC, THISLOC=LOC )
        2497           0 :            RETURN
        2498             :        ENDIF
        2499             :
        2500             :        ! Data is 2D
        2501           0 :        Lct%Dct%Dta%SpaceDim = 2
        2502             :
        2503             :     ! For base emissions and scale factors, interpret data as scalar
        2504             :     ! values with a time dimension.
        2505             :     ELSE
        2506             :
        2507           0 :        CALL FileData_ArrCheck( HcoState%Config, Lct%Dct%Dta, 1, 1, NT, RC )
        2508           0 :        IF ( RC /= HCO_SUCCESS ) THEN
        2509           0 :            CALL HCO_ERROR( 'ERROR 17', RC, THISLOC=LOC )
        2510           0 :            RETURN
        2511             :        ENDIF
        2512           0 :        DO I = 1, NT
        2513           0 :           Lct%Dct%Dta%V2(I)%Val(1,1) = Vals(I)
        2514             : !==============================================================================
        2515             : ! KLUDGE BY BOB YANTOSCA (05 Jan 2016)
        2516             : !
        2517             : ! This WRITE statement avoids a seg fault in some Intel Fortran Compiler
        2518             : ! versions, such as ifort 12 and ifort 13.  The ADVANCE="no" prevents
        2519             : ! carriage returns from being added to the log file, and the '' character
        2520             : ! will prevent text from creeping across the screen.
        2521             : !
        2522             : ! NOTE: This section only gets executed during the initialization phase,
        2523             : ! when we save data not read from netCDF files into the HEMCO data structure.
        2524             : ! This type of data includes scale factors and mask data specified as vectors
        2525             : ! in the HEMCO configuration file.  Therefore, this section will only get
        2526             : ! executed at startup, so the WRITE statment should not add significant
        2527             : ! overhead to the simulation.
        2528             : !
        2529             : ! The root issue seems to be an optimization bug in the compiler.
        2530             : !==============================================================================
        2531             : #if defined( LINUX_IFORT )
        2532             :           WRITE( 6, '(a)', ADVANCE='no' ) ''
        2533             : #endif
        2534             :
        2535             :        ENDDO
        2536             :
        2537             :        ! Data is 1D
        2538           0 :        Lct%Dct%Dta%SpaceDim  = 1
        2539             :
        2540             :        ! Make sure data is in local time
        2541           0 :        IF ( .NOT. Lct%Dct%Dta%IsLocTime ) THEN
        2542           0 :           Lct%Dct%Dta%IsLocTime = .TRUE.
        2543             :           MSG = 'Scale factors read from file are treated as local time: '// &
        2544           0 :                  TRIM(Lct%Dct%cName)
        2545           0 :           CALL HCO_WARNING( HcoState%Config%Err, MSG, RC, THISLOC=LOC )
        2546             :        ENDIF
        2547             :
        2548             :     ENDIF
        2549             :
        2550             :     ! Cleanup
        2551           0 :     IF ( ASSOCIATED(Vals) ) DEALLOCATE(Vals)
        2552             :
        2553             :     ! Return w/ success
        2554           0 :     RC = HCO_SUCCESS
        2555             :
        2556           0 :   END SUBROUTINE HCOIO_ReadFromConfig
        2557             : !EOC
        2558             : !------------------------------------------------------------------------------
        2559             : !                   Harmonized Emissions Component (HEMCO)                    !
        2560             : !------------------------------------------------------------------------------
        2561             : !BOP
        2562             : !
        2563             : ! !IROUTINE: GetSliceIdx
        2564             : !
        2565             : ! !DESCRIPTION: gets the time slice index to be used for data directly
        2566             : ! read from the HEMCO configuration file. prefDt denotes the preferred
        2567             : ! time attribute (year, month, or day). DtType is used to identify the
        2568             : ! time attribute type (1=year, 2=month, 3=day). The time slice index will
        2569             : ! be selected based upon those two variables. IDX is the selected time
        2570             : ! slice index. It will be set to -1 if the current simulation date
        2571             : ! is outside of the specified time range and the time cycle attribute is
        2572             : ! not enabled for this field.
        2573             : !\\
        2574             : !\\
        2575             : ! !INTERFACE:
        2576             : !
        2577           0 :   SUBROUTINE GetSliceIdx ( HcoState, Lct, DtType, prefDt, IDX, RC )
        2578             : !
        2579             : ! !INPUT PARAMETERS:
        2580             : !
        2581             :     TYPE(HCO_State),  POINTER                 :: HcoState
        2582             :     TYPE(ListCont),   POINTER                 :: Lct
        2583             :     INTEGER,          INTENT(IN   )           :: DtType
        2584             :     INTEGER,          INTENT(IN   )           :: prefDt
        2585             : !
        2586             : ! !INPUT/OUTPUT PARAMETERS:
        2587             : !
        2588             :     INTEGER,          INTENT(INOUT)           :: IDX
        2589             :     INTEGER,          INTENT(INOUT)           :: RC
        2590             : !
        2591             : ! !REVISION HISTORY:
        2592             : !  13 Mar 2013 - C. Keller - Initial version
        2593             : !  See https://github.com/geoschem/hemco for complete history
        2594             : !EOP
        2595             : !------------------------------------------------------------------------------
        2596             : !BOC
        2597             : !
        2598             : ! !LOCAL VARIABLES:
        2599             : !
        2600             :     INTEGER            :: lowDt, uppDt
        2601             :     CHARACTER(LEN=255) :: MSG
        2602             :     CHARACTER(LEN=255) :: LOC = 'GetSliceIdx (hcoio_util_mod.F90)'
        2603             :
        2604             :     !=================================================================
        2605             :     ! GetSliceIdx begins here!
        2606             :     !=================================================================
        2607             :
        2608             :     ! Init
        2609           0 :     RC = HCO_SUCCESS
        2610             :
        2611             :     ! Get upper and lower time range
        2612           0 :     IF ( DtType == 1 ) THEN
        2613           0 :        lowDt = Lct%Dct%Dta%ncYrs(1)
        2614           0 :        uppDt = Lct%Dct%Dta%ncYrs(2)
        2615           0 :     ELSEIF ( DtType == 2 ) THEN
        2616           0 :        lowDt = Lct%Dct%Dta%ncMts(1)
        2617           0 :        uppDt = Lct%Dct%Dta%ncMts(2)
        2618           0 :     ELSEIF ( DtType == 3 ) THEN
        2619           0 :        lowDt = Lct%Dct%Dta%ncDys(1)
        2620           0 :        uppDt = Lct%Dct%Dta%ncDys(2)
        2621             :     ELSE
        2622           0 :        WRITE(MSG,*) "DtType must be one of 1, 2, 3: ", DtType
        2623           0 :        CALL HCO_ERROR( MSG, RC, THISLOC=LOC )
        2624           0 :        RETURN
        2625             :     ENDIF
        2626             :
        2627             :     ! Check for cycle flags:
        2628             :
        2629             :     ! Data cycle set to range or exact date: in these cases, the
        2630             :     ! the preferred date will be equal to the current date, so
        2631             :     ! check if the preferred date is indeed within the available
        2632             :     ! range (lowDt, uppDt).
        2633             :     ! For data only to be used within the specified range, set
        2634             :     ! index to -1. This will force the scale factors to be set to
        2635             :     ! zero!
        2636           0 :     IF ( prefDt < lowDt .OR. prefDt > uppDt ) THEN
        2637           0 :        IF ( ( Lct%Dct%Dta%CycleFlag == HCO_CFLAG_EXACT ) .OR.      &
        2638             :             ( Lct%Dct%Dta%CycleFlag == HCO_CFLAG_RANGE )     ) THEN
        2639           0 :           IDX = -1
        2640           0 :           RETURN
        2641             :        ELSE
        2642             :           ! this here should never happen, since for a cycle flag of 1,
        2643             :           ! the preferred date should always be restricted to the range
        2644             :           ! of available time stamps.
        2645           0 :           MSG = 'preferred date is outside of range: ' // TRIM(Lct%Dct%cName)
        2646           0 :           CALL HCO_ERROR( MSG, RC, THISLOC=LOC )
        2647           0 :           RETURN
        2648             :        ENDIF
        2649             :     ENDIF
        2650             :
        2651             :     ! If the code makes it to here, prefDt is within the available data range
        2652             :     ! and we simply get the wanted index from the current index and the lowest
        2653             :     ! available index.
        2654           0 :     IDX = prefDt - lowDt + 1
        2655             :
        2656             :   END SUBROUTINE GetSliceIdx
        2657             : !EOC
        2658             : !------------------------------------------------------------------------------
        2659             : !                   Harmonized Emissions Component (HEMCO)                    !
        2660             : !------------------------------------------------------------------------------
        2661             : !BOP
        2662             : !
        2663             : ! !IROUTINE: GetDataVals
        2664             : !
        2665             : ! !DESCRIPTION: Subroutine GetDataVals extracts the data values from ValStr
        2666             : ! and writes them into vector Vals. ValStr is typically a character string
        2667             : ! read from an external ASCII file or directly from the HEMCO configuration
        2668             : ! file. Depending on the time specifications provided in the configuration
        2669             : ! file, Vals will be filled with only a subset of the values of ValStr.
        2670             : !\\
        2671             : !\\
        2672             : ! !INTERFACE:
        2673             : !
        2674           0 :   SUBROUTINE GetDataVals ( HcoState, Lct, ValStr, Vals, RC )
        2675             : !
        2676             : ! !USES:
        2677             : !
        2678             :     USE HCO_CHARTOOLS_MOD,  ONLY : HCO_CharSplit
        2679             :     USE HCO_EXTLIST_MOD,    ONLY : HCO_GetOpt
        2680             :     USE HCO_UNIT_MOD,       ONLY : HCO_Unit_Change
        2681             :     USE HCO_tIdx_Mod,       ONLY : HCO_GetPrefTimeAttr
        2682             :     USE HCO_CLOCK_MOD,      ONLY : HcoClock_Get
        2683             : !
        2684             : ! !INPUT PARAMTERS:
        2685             : !
        2686             :     TYPE(HCO_State),  POINTER         :: HcoState    ! HEMCO state
        2687             :     CHARACTER(LEN=*), INTENT(IN   )   :: ValStr
        2688             : !
        2689             : ! !INPUT/OUTPUT PARAMETERS:
        2690             : !
        2691             :     TYPE(ListCont),   POINTER         :: Lct
        2692             :     INTEGER,          INTENT(INOUT)   :: RC
        2693             : !
        2694             : ! !OUTPUT PARAMETERS:
        2695             : !
        2696             :     REAL(hp),         POINTER         :: Vals(:)
        2697             : !
        2698             : ! !REVISION HISTORY:
        2699             : !  22 Dec 2014 - C. Keller: Initial version
        2700             : !  See https://github.com/geoschem/hemco for complete history
        2701             : !EOP
        2702             : !------------------------------------------------------------------------------
        2703             : !BOC
        2704             : !
        2705             : ! !LOCAL VARIABLES:
        2706             : !
        2707             :     INTEGER            :: HcoID
        2708             :     INTEGER            :: I, N, NUSE, AS
        2709             :     INTEGER            :: IDX1, IDX2
        2710             :     INTEGER            :: AreaFlag, TimeFlag, Check
        2711             :     INTEGER            :: prefYr, prefMt, prefDy, prefHr, prefMn
        2712             :     INTEGER            :: cYr,    cMt,    cDy,    cHr
        2713             :     REAL(hp)           :: MW_g
        2714             :     REAL(hp)           :: UnitFactor
        2715             :     REAL(hp)           :: FileVals(100)
        2716           0 :     REAL(hp), POINTER  :: FileArr(:,:,:,:)
        2717             :     LOGICAL            :: IsPerArea
        2718             :     LOGICAL            :: IsMath
        2719             :     CHARACTER(LEN=255) :: MSG
        2720             :     CHARACTER(LEN=255) :: LOC = 'GetDataVals (hcoio_util_mod.F90)'
        2721             :
        2722             :     !======================================================================
        2723             :     ! GetDataVals begins here
        2724             :     !======================================================================
        2725             :
        2726             :     ! Initialize
        2727           0 :     FileArr => NULL()
        2728             :
        2729             :     ! Shadow species properties needed for unit conversion
        2730           0 :     HcoID = Lct%Dct%HcoID
        2731           0 :     IF ( HcoID > 0 ) THEN
        2732           0 :        MW_g = HcoState%Spc(HcoID)%MW_g
        2733             :     ELSE
        2734           0 :        MW_g = -999.0_hp
        2735             :     ENDIF
        2736             :
        2737             :     ! Is this a math expression?
        2738           0 :     IsMath = .FALSE.
        2739           0 :     IF ( LEN(ValStr) > 5 ) THEN
        2740           0 :        IF ( ValStr(1:5)=='MATH:' ) IsMath = .TRUE.
        2741             :     ENDIF
        2742             :
        2743             :     ! Evaluate math expression if string starts with 'MATH:'
        2744             :     IF ( IsMath ) THEN
        2745           0 :        CALL ReadMath ( HcoState, Lct, ValStr, FileVals, N, RC )
        2746           0 :        IF ( RC /= HCO_SUCCESS ) THEN
        2747           0 :            CALL HCO_ERROR( 'ERROR 18', RC, THISLOC=LOC )
        2748           0 :            RETURN
        2749             :        ENDIF
        2750             :
        2751             :     ! Use regular string parser otherwise
        2752             :     ELSE
        2753             :        CALL HCO_CharSplit ( ValStr, &
        2754             :                             HCO_GetOpt(HcoState%Config%ExtList,'Separator'), &
        2755             :                             HCO_GetOpt(HcoState%Config%ExtList,'Wildcard'), &
        2756           0 :                             FileVals, N, RC )
        2757           0 :        IF ( RC /= HCO_SUCCESS ) THEN
        2758           0 :            CALL HCO_ERROR( 'ERROR 19', RC, THISLOC=LOC )
        2759           0 :            RETURN
        2760             :        ENDIF
        2761             :     ENDIF
        2762             :
        2763             :     ! Return w/ error if no scale factor defined
        2764           0 :     IF ( N == 0 ) THEN
        2765             :        MSG = 'Cannot read data: ' // TRIM(Lct%Dct%cName) // &
        2766           0 :              ': ' // TRIM(ValStr)
        2767           0 :        CALL HCO_ERROR( MSG, RC, THISLOC=LOC)
        2768           0 :        RETURN
        2769             :     ENDIF
        2770             :
        2771             :     ! Get the preferred times, i.e. the preferred year, month, day,
        2772             :     ! or hour (as specified in the configuration file).
        2773             :     CALL HCO_GetPrefTimeAttr( HcoState, Lct, &
        2774           0 :                               prefYr, prefMt, prefDy, prefHr, prefMn, RC )
        2775           0 :     IF ( RC /= HCO_SUCCESS ) THEN
        2776           0 :         CALL HCO_ERROR( 'ERROR 20', RC, THISLOC=LOC )
        2777           0 :         RETURN
        2778             :     ENDIF
        2779             :
        2780             :     ! ----------------------------------------------------------------
        2781             :     ! For masks, assume that values represent the corners of the mask
        2782             :     ! box, e.g. there must be four values. Masks are time-independent
        2783             :     ! and unitless
        2784             :     ! ----------------------------------------------------------------
        2785           0 :     IF ( Lct%Dct%DctType == HCO_DCTTYPE_MASK ) THEN
        2786             :
        2787             :        ! There must be exactly four values
        2788           0 :        IF ( N /= 4 ) THEN
        2789             :           MSG = 'Mask values are not lon1/lat1/lon2/lat2: ' // &
        2790           0 :                 TRIM(ValStr) // ' --> ' // TRIM(Lct%Dct%cName)
        2791           0 :           CALL HCO_ERROR( MSG, RC, THISLOC=LOC )
        2792           0 :           RETURN
        2793             :        ENDIF
        2794             :
        2795             :        ! Pass to FileArr array (will be used below)
        2796           0 :        NUSE = 4
        2797           0 :        ALLOCATE( FileArr(1,1,1,NUSE), STAT=AS )
        2798             :        IF ( AS /= 0 ) THEN
        2799           0 :           MSG = 'Cannot allocate FileArr'
        2800           0 :           CALL HCO_ERROR( MSG, RC, THISLOC=LOC )
        2801           0 :           RETURN
        2802             :        ENDIF
        2803           0 :        FileArr(1,1,1,:) = FileVals(1:NUSE)
        2804             :
        2805             :     ! ----------------------------------------------------------------
        2806             :     ! For non-masks, the data is interpreted as uniform values with
        2807             :     ! a time dimension. Need to select the time slices to be used at
        2808             :     ! this time (depending on the provided time attributes), as well
        2809             :     ! as to ensure that values are in the correct units.
        2810             :     ! Use all time slices unless a time interval is provided in
        2811             :     ! attribute srcTime of the configuration file.
        2812             :     ! ----------------------------------------------------------------
        2813             :     ELSE
        2814             :
        2815             :        ! If there is only one value use this one and ignore any time
        2816             :        ! preferences.
        2817           0 :        IF ( N == 1 ) THEN
        2818           0 :           NUSE = 1
        2819           0 :           IDX1 = 1
        2820           0 :           IDX2 = 1
        2821             :
        2822             :        ! If it's a math expression use all passed values
        2823           0 :        ELSEIF ( IsMath ) THEN
        2824           0 :           NUSE = N
        2825           0 :           IDX1 = 1
        2826           0 :           IDX2 = N
        2827             :
        2828             :        ELSE
        2829             :           ! Currently, data read directly from the configuration file can only
        2830             :           ! represent one time dimension, i.e. it can only be yearly, monthly,
        2831             :           ! daily (or hourly data, but this is read all at the same time).
        2832             :
        2833             :           ! Annual data
        2834           0 :           IF ( Lct%Dct%Dta%ncYrs(1) /= Lct%Dct%Dta%ncYrs(2) ) THEN
        2835             :              ! Error check
        2836             :              IF ( Lct%Dct%Dta%ncMts(1) /= Lct%Dct%Dta%ncMts(2) .OR. &
        2837           0 :                   Lct%Dct%Dta%ncDys(1) /= Lct%Dct%Dta%ncDys(2) .OR. &
        2838             :                   Lct%Dct%Dta%ncHrs(1) /= Lct%Dct%Dta%ncHrs(2)       ) THEN
        2839             :                 MSG = 'Data must not have more than one time dimension: ' // &
        2840           0 :                        TRIM(Lct%Dct%cName)
        2841           0 :                 CALL HCO_ERROR( MSG, RC, THISLOC=LOC )
        2842           0 :                 RETURN
        2843             :              ENDIF
        2844             :
        2845           0 :              CALL GetSliceIdx ( HcoState, Lct, 1, prefYr, IDX1, RC )
        2846           0 :              IF ( RC /= HCO_SUCCESS ) THEN
        2847           0 :                  CALL HCO_ERROR( 'ERROR 21', RC, THISLOC=LOC )
        2848           0 :                  RETURN
        2849             :              ENDIF
        2850           0 :              IDX2 = IDX1
        2851           0 :              NUSE = 1
        2852             :
        2853             :           ! Monthly data
        2854           0 :           ELSEIF ( Lct%Dct%Dta%ncMts(1) /= Lct%Dct%Dta%ncMts(2) ) THEN
        2855             :              ! Error check
        2856           0 :              IF ( Lct%Dct%Dta%ncDys(1) /= Lct%Dct%Dta%ncDys(2) .OR. &
        2857             :                   Lct%Dct%Dta%ncHrs(1) /= Lct%Dct%Dta%ncHrs(2)       ) THEN
        2858             :                 MSG = 'Data must only have one time dimension: ' // &
        2859           0 :                       TRIM(Lct%Dct%cName)
        2860           0 :                 CALL HCO_ERROR( MSG, RC, THISLOC=LOC )
        2861           0 :                 RETURN
        2862             :              ENDIF
        2863             :
        2864           0 :              CALL GetSliceIdx ( HcoState, Lct, 2, prefMt, IDX1, RC )
        2865           0 :              IF ( RC /= HCO_SUCCESS ) THEN
        2866           0 :                  CALL HCO_ERROR( 'ERROR 22', RC, THISLOC=LOC )
        2867           0 :                  RETURN
        2868             :              ENDIF
        2869           0 :              IDX2 = IDX1
        2870           0 :              NUSE = 1
        2871             :
        2872             :           ! Daily data
        2873           0 :           ELSEIF ( Lct%Dct%Dta%ncDys(1) /= Lct%Dct%Dta%ncDys(2) ) THEN
        2874             :              ! Error check
        2875           0 :              IF ( Lct%Dct%Dta%ncHrs(1) /= Lct%Dct%Dta%ncHrs(2) ) THEN
        2876             :                 MSG = 'Data must only have one time dimension: ' // &
        2877           0 :                       TRIM(Lct%Dct%cName)
        2878           0 :                 CALL HCO_ERROR( MSG, RC, THISLOC=LOC )
        2879           0 :                 RETURN
        2880             :              ENDIF
        2881             :
        2882           0 :              CALL GetSliceIdx ( HcoState, Lct, 3, prefDy, IDX1, RC )
        2883           0 :              IF ( RC /= HCO_SUCCESS ) THEN
        2884           0 :                  CALL HCO_ERROR( 'ERROR 23', RC, THISLOC=LOC )
        2885           0 :                  RETURN
        2886             :              ENDIF
        2887           0 :              IDX2 = IDX1
        2888           0 :              NUSE = 1
        2889             :
        2890             :           ! All other cases (incl. hourly data): read all time slices).
        2891             :           ELSE
        2892           0 :              IDX1 = 1
        2893           0 :              IDX2 = N
        2894           0 :              NUSE = N
        2895             :           ENDIF
        2896             :        ENDIF
        2897             :
        2898             :        ! ----------------------------------------------------------------
        2899             :        ! Read selected time slice(s) into data array
        2900             :        ! ----------------------------------------------------------------
        2901           0 :        IF ( IDX2 > N ) THEN
        2902           0 :           WRITE(MSG,*) 'Index ', IDX2, ' is larger than number of ', &
        2903           0 :                        'values found: ', TRIM(Lct%Dct%cName)
        2904           0 :           CALL HCO_ERROR( MSG, RC, THISLOC=LOC )
        2905           0 :           RETURN
        2906             :        ENDIF
        2907             :
        2908           0 :        ALLOCATE( FileArr(1,1,1,NUSE), STAT=AS )
        2909             :        IF ( AS /= 0 ) THEN
        2910           0 :           MSG = 'Cannot allocate FileArr'
        2911           0 :           CALL HCO_ERROR( MSG, RC, THISLOC=LOC )
        2912           0 :           RETURN
        2913             :        ENDIF
        2914             :
        2915             :        ! Check for range/exact flag
        2916             :        ! If range is given, the preferred Yr/Mt/Dy/Hr will be negative
        2917             :        ! if we are outside the desired range.
        2918           0 :        IF ( Lct%Dct%Dta%CycleFlag == HCO_CFLAG_RANGE ) THEN
        2919           0 :           IF ( prefYr == -1 .OR. prefMt == -1 .OR. prefDy == -1 ) IDX1 = -1
        2920           0 :           IF ( Lct%Dct%Dta%ncHrs(1) >= 0 .AND. prefHr == -1 )     IDX1 = -1
        2921             :
        2922             :        ! If flag is exact, the preferred date must be equal to the current
        2923             :        ! simulation date.
        2924           0 :        ELSEIF ( Lct%Dct%Dta%CycleFlag == HCO_CFLAG_EXACT ) THEN
        2925           0 :           IF ( Lct%Dct%Dta%ncYrs(1) > 0 ) THEN
        2926           0 :              IF ( prefYr < Lct%Dct%Dta%ncYrs(1) .OR. &
        2927           0 :                   prefYr > Lct%Dct%Dta%ncYrs(2) ) IDX1 = -1
        2928             :           ENDIF
        2929           0 :           IF ( Lct%Dct%Dta%ncMts(1) > 0 ) THEN
        2930           0 :              IF ( prefMt < Lct%Dct%Dta%ncMts(1) .OR. &
        2931           0 :                   prefMt > Lct%Dct%Dta%ncMts(2) ) IDX1 = -1
        2932             :           ENDIF
        2933           0 :           IF ( Lct%Dct%Dta%ncDys(1) > 0 ) THEN
        2934           0 :              IF ( prefDy < Lct%Dct%Dta%ncDys(1) .OR. &
        2935           0 :                   prefDy > Lct%Dct%Dta%ncDys(2) ) IDX1 = -1
        2936             :           ENDIF
        2937           0 :           IF ( Lct%Dct%Dta%ncHrs(1) >= 0 ) THEN
        2938           0 :              IF ( prefHr < Lct%Dct%Dta%ncHrs(1) .OR. &
        2939           0 :                   prefHr > Lct%Dct%Dta%ncHrs(2) ) IDX1 = -1
        2940             :           ENDIF
        2941             :        ENDIF
        2942             :
        2943             :        ! IDX1 becomes -1 for data that is outside of the valid range
        2944             :        ! (and no time cycling enabled). In this case, make sure that
        2945             :        ! scale factor is set to zero.
        2946           0 :        IF ( IDX1 < 0 ) THEN
        2947           0 :           IF ( Lct%Dct%DctType == HCO_DCTTYPE_BASE ) THEN
        2948           0 :              FileArr(1,1,1,:) = 0.0_hp
        2949             :              MSG = 'Base field outside of range - set to zero: ' // &
        2950           0 :                    TRIM(Lct%Dct%cName)
        2951           0 :              CALL HCO_WARNING ( HcoState%Config%Err, MSG, RC, THISLOC=LOC )
        2952             : #if defined( MODEL_GEOS )
        2953             :           ELSEIF ( Lct%Dct%DctType == HCO_DCTTYPE_MASK ) THEN
        2954             :              FileArr(1,1,1,:) = 0.0_hp
        2955             :              MSG = 'Mask outside of range - set to zero: ' // &
        2956             :                    TRIM(Lct%Dct%cName)
        2957             :              CALL HCO_WARNING ( HcoState%Config%Err, MSG, RC, THISLOC=LOC )
        2958             : #endif
        2959             :           ELSE
        2960           0 :              FileArr(1,1,1,:) = 1.0_hp
        2961             :              MSG = 'Scale factor outside of range - set to one: ' // &
        2962           0 :                    TRIM(Lct%Dct%cName)
        2963           0 :              CALL HCO_WARNING ( HcoState%Config%Err, MSG, RC, THISLOC=LOC )
        2964             :           ENDIF
        2965             :        ELSE
        2966           0 :           FileArr(1,1,1,:) = FileVals(IDX1:IDX2)
        2967             :        ENDIF
        2968             :
        2969             :        ! ----------------------------------------------------------------
        2970             :        ! Convert data to HEMCO units
        2971             :        ! ----------------------------------------------------------------
        2972             :        CALL HCO_UNIT_CHANGE( HcoConfig     = HcoState%Config,            &
        2973             :                              Array         = FileArr,                    &
        2974             :                              Units         = TRIM(Lct%Dct%Dta%OrigUnit), &
        2975             :                              MW            = MW_g,                       &
        2976             :                              YYYY          = -999,                       &
        2977             :                              MM            = -999,                       &
        2978             :                              AreaFlag      = AreaFlag,                   &
        2979             :                              TimeFlag      = TimeFlag,                   &
        2980             :                              FACT          = UnitFactor,                 &
        2981           0 :                              RC            = RC                           )
        2982           0 :        IF ( RC /= HCO_SUCCESS ) THEN
        2983           0 :            CALL HCO_ERROR( 'ERROR 24', RC, THISLOC=LOC )
        2984           0 :            RETURN
        2985             :        ENDIF
        2986             :
        2987             :        ! Verbose mode
        2988           0 :        IF ( UnitFactor /= 1.0_hp ) THEN
        2989           0 :           IF ( HCO_IsVerb( HcoState%Config%Err ) ) THEN
        2990           0 :              WRITE(MSG,*) 'Data was in units of ', TRIM(Lct%Dct%Dta%OrigUnit), &
        2991           0 :                           ' - converted to HEMCO units by applying ', &
        2992           0 :                           'scale factor ', UnitFactor
        2993           0 :              CALL HCO_MSG(HcoState%Config%Err,MSG)
        2994             :           ENDIF
        2995             :        ELSE
        2996           0 :           IF ( HCO_IsVerb( HcoState%Config%Err ) ) THEN
        2997           0 :              WRITE(MSG,*) 'Data was in units of ', TRIM(Lct%Dct%Dta%OrigUnit), &
        2998           0 :                           ' - unit conversion factor is ', UnitFactor
        2999           0 :              CALL HCO_MSG(HcoState%Config%Err,MSG)
        3000             :           ENDIF
        3001             :        ENDIF
        3002             :
        3003             :        ! Data must be ...
        3004             :        ! ... concentration ...
        3005           0 :        IF ( AreaFlag == 3 .AND. TimeFlag == 0 ) THEN
        3006           0 :           Lct%Dct%Dta%IsConc = .TRUE.
        3007             :
        3008           0 :        ELSEIF ( AreaFlag == 3 .AND. TimeFlag == 1 ) THEN
        3009           0 :           Lct%Dct%Dta%IsConc = .TRUE.
        3010           0 :           FileArr = FileArr * HcoState%TS_EMIS
        3011             :           MSG = 'Data converted from kg/m3/s to kg/m3: ' // &
        3012           0 :                 TRIM(Lct%Dct%cName) // ': ' // TRIM(Lct%Dct%Dta%OrigUnit)
        3013           0 :           CALL HCO_WARNING ( HcoState%Config%Err, MSG, RC, THISLOC=LOC )
        3014             :
        3015             :        ! ... emissions or unitless ...
        3016           0 :        ELSEIF ( (AreaFlag == -1 .AND. TimeFlag == -1) .OR. &
        3017             :                 (AreaFlag ==  2 .AND. TimeFlag ==  1)       ) THEN
        3018           0 :           Lct%Dct%Dta%IsConc = .FALSE.
        3019             :
        3020             :        ! ... invalid otherwise:
        3021             :        ELSE
        3022             :           MSG = 'Unit must be unitless, emission or concentration: ' // &
        3023           0 :                 TRIM(Lct%Dct%cName) // ': ' // TRIM(Lct%Dct%Dta%OrigUnit)
        3024           0 :           CALL HCO_ERROR( MSG, RC, THISLOC=LOC )
        3025           0 :           RETURN
        3026             :        ENDIF
        3027             :
        3028             :        ! Auto-detect delta t [in hours] between time slices.
        3029             :        ! Scale factors can be:
        3030             :        ! length 1 : constant
        3031             :        ! length 7 : weekday factors: Sun, Mon, ..., Sat
        3032             :        ! length 12: monthly factors: Jan, Feb, ..., Dec
        3033             :        ! length 24: hourly  factors: 12am, 1am, ... 11pm
        3034           0 :        IF ( NUSE == 1 ) THEN
        3035           0 :           Lct%Dct%Dta%DeltaT = 0
        3036           0 :        ELSEIF ( NUSE == 7 ) THEN
        3037           0 :           Lct%Dct%Dta%DeltaT = 24
        3038           0 :        ELSEIF ( NUSE == 12 ) THEN
        3039           0 :           Lct%Dct%Dta%DeltaT = 720
        3040           0 :        ELSEIF ( NUSE == 24 ) THEN
        3041           0 :           Lct%Dct%Dta%DeltaT = 1
        3042             :        ELSE
        3043             :           MSG = 'Factor must be of length 1, 7, 12, or 24!' // &
        3044           0 :                  TRIM(Lct%Dct%cName)
        3045           0 :           CALL HCO_ERROR( MSG, RC, THISLOC=LOC)
        3046           0 :           RETURN
        3047             :        ENDIF
        3048             :
        3049             :     ENDIF ! Masks vs. non-masks
        3050             :
        3051             :     ! Copy data into output array.
        3052           0 :     IF ( ASSOCIATED(Vals) ) DEALLOCATE( Vals )
        3053           0 :     ALLOCATE( Vals(NUSE), STAT=AS )
        3054             :     IF ( AS /= 0 ) THEN
        3055           0 :        MSG = 'Cannot allocate Vals'
        3056           0 :        CALL HCO_ERROR( MSG, RC, THISLOC=LOC )
        3057           0 :        RETURN
        3058             :     ENDIF
        3059           0 :     Vals(:) = FileArr(1,1,1,:)
        3060             :
        3061             :     ! Cleanup
        3062           0 :     IF ( ASSOCIATED(FileArr) ) DEALLOCATE(FileArr)
        3063           0 :     FileArr => NULL()
        3064             :
        3065             :     ! Return w/ success
        3066           0 :     RC = HCO_SUCCESS
        3067             :
        3068           0 :   END SUBROUTINE GetDataVals
        3069             : !EOC
        3070             : !------------------------------------------------------------------------------
        3071             : !                   Harmonized Emissions Component (HEMCO)                    !
        3072             : !------------------------------------------------------------------------------
        3073             : !BOP
        3074             : !
        3075             : ! !IROUTINE: FillMaskBox
        3076             : !
        3077             : ! !DESCRIPTION: Subroutine FillMaskBox fills the data array of the passed list
        3078             : ! container Lct according to the mask region provided in Vals. Vals contains
        3079             : ! the mask region of interest, denoted by the lower left and upper right grid
        3080             : ! box corners: lon1, lat1, lon2, lat2. The data array of Lct is filled such
        3081             : ! that all grid boxes are set to 1 whose mid-point is inside of the given box
        3082             : ! range.
        3083             : !\\
        3084             : !\\
        3085             : ! !INTERFACE:
        3086             : !
        3087           0 :   SUBROUTINE FillMaskBox ( HcoState, Lct, Vals, RC )
        3088             : !
        3089             : ! !USES:
        3090             : !
        3091             : !
        3092             : ! !INPUT PARAMTERS:
        3093             : !
        3094             :     TYPE(HCO_State),  POINTER         :: HcoState    ! HEMCO state
        3095             :     REAL(hp)        , POINTER         :: Vals(:)
        3096             : !
        3097             : ! !INPUT/OUTPUT PARAMETERS:
        3098             : !
        3099             :     TYPE(ListCont),   POINTER         :: Lct
        3100             :     INTEGER,          INTENT(INOUT)   :: RC
        3101             : !
        3102             : ! !REVISION HISTORY:
        3103             : !  29 Dec 2014 - C. Keller - Initial version
        3104             : !  See https://github.com/geoschem/hemco for complete history
        3105             : !EOP
        3106             : !------------------------------------------------------------------------------
        3107             : !BOC
        3108             : !
        3109             : ! !LOCAL VARIABLES:
        3110             : !
        3111             :     LOGICAL            :: GridPoint
        3112             :     INTEGER            :: I, J
        3113             :     REAL(hp)           :: LON1, LON2, LAT1, LAT2
        3114             :     REAL(hp)           :: XDG1, XDG2, YDG1, YDG2
        3115             :     REAL(hp)           :: ILON, ILAT
        3116             :     CHARACTER(LEN=255) :: MSG
        3117             :     CHARACTER(LEN=255) :: LOC = 'FillMaskBox (hcoio_util_mod.F90)'
        3118             :
        3119             :     !=================================================================
        3120             :     ! FillMaskBox begins here!
        3121             :     !=================================================================
        3122             :
        3123             :     ! Extract lon1, lon2, lat1, lat2
        3124           0 :     LON1 = VALS(1)
        3125           0 :     LAT1 = VALS(2)
        3126           0 :     LON2 = VALS(3)
        3127           0 :     LAT2 = VALS(4)
        3128             :
        3129             :     ! Check if this is mask is a point. In this case, we need the grid
        3130             :     ! box edges being defined.
        3131           0 :     GridPoint = .FALSE.
        3132           0 :     IF ( ( LON1 == LON2 ) .AND. ( LAT1 == LAT2 ) ) THEN
        3133           0 :        IF ( .NOT. ASSOCIATED(HcoState%Grid%XEDGE%Val) .OR. &
        3134             :             .NOT. ASSOCIATED(HcoState%Grid%YEDGE%Val)       ) THEN
        3135             :           MSG = 'Cannot evaluate grid point mask - need grid box '   // &
        3136             :                 'edges for this. This error occurs if a mask covers '// &
        3137             :                 'a fixed grid point (e.g. lon1=lon2 and lat1=lat2) ' // &
        3138           0 :                 'but HEMCO grid edges are not defined.'
        3139           0 :           CALL HCO_ERROR( MSG, RC, THISLOC=LOC )
        3140           0 :           RETURN
        3141             :        ENDIF
        3142             :        GridPoint = .TRUE.
        3143             :     ENDIF
        3144             :
        3145             :     ! Check for every grid box if mid point is within mask region.
        3146             :     ! Set to 1.0 if this is the case.
        3147             : !$OMP PARALLEL DO                        &
        3148             : !$OMP DEFAULT( SHARED                 )  &
        3149             : !$OMP PRIVATE( I, J, ILON, ILAT       )  &
        3150             : !$OMP PRIVATE( XDG1, XDG2, YDG1, YDG2 )  &
        3151             : !$OMP SCHEDULE( DYNAMIC               )
        3152           0 :     DO J = 1, HcoState%NY
        3153           0 :     DO I = 1, HcoState%NX
        3154             :
        3155             :        ! If it's a grid point, check if it's within this
        3156             :        ! grid box
        3157           0 :        IF ( GridPoint ) THEN
        3158           0 :           XDG1 = HcoState%Grid%XEDGE%Val(I  ,J  )
        3159           0 :           XDG2 = HcoState%Grid%XEDGE%Val(I+1,J  )
        3160           0 :           YDG1 = HcoState%Grid%YEDGE%Val(I  ,J  )
        3161           0 :           YDG2 = HcoState%Grid%YEDGE%Val(I  ,J+1)
        3162           0 :           IF ( XDG1 >= 180.0_hp ) XDG1 = XDG1 - 360.0_hp
        3163           0 :           IF ( XDG2 >= 180.0_hp ) XDG2 = XDG2 - 360.0_hp
        3164             :
        3165             :           IF ( LON1 >= XDG1 .AND. LON1 <= XDG2 .AND. &
        3166           0 :                LAT1 >= YDG1 .AND. LAT1 <= YDG2        ) THEN
        3167           0 :              Lct%Dct%Dta%V2(1)%Val(I,J) = 1.0_sp
        3168             :           ENDIF
        3169             :
        3170             :        ! Check if mid point is within mask region
        3171             :        ELSE
        3172             :           ! Get longitude and latitude at this grid box
        3173           0 :           ILON = HcoState%Grid%XMID%Val(I,J)
        3174           0 :           ILAT = HcoState%Grid%YMID%Val(I,J)
        3175           0 :           IF ( ILON >= 180.0_hp ) ILON = ILON - 360.0_hp
        3176             :
        3177             :           IF ( ILON >= LON1 .AND. ILON <= LON2 .AND. &
        3178           0 :                ILAT >= LAT1 .AND. ILAT <= LAT2        ) THEN
        3179           0 :              Lct%Dct%Dta%V2(1)%Val(I,J) = 1.0_sp
        3180             :           ENDIF
        3181             :        ENDIF
        3182             :
        3183             :     ENDDO
        3184             :     ENDDO
        3185             : !$OMP END PARALLEL DO
        3186             :
        3187             :     ! Return w/ success
        3188           0 :     RC = HCO_SUCCESS
        3189             :
        3190             :   END SUBROUTINE FillMaskBox
        3191             : !EOC
        3192             : !------------------------------------------------------------------------------
        3193             : !                   Harmonized Emissions Component (HEMCO)                    !
        3194             : !------------------------------------------------------------------------------
        3195             : !BOP
        3196             : !
        3197             : ! !IROUTINE: ReadMath
        3198             : !
        3199             : ! !DESCRIPTION: Subroutine ReadMath reads and evaluates a mathematical
        3200             : ! expression. Mathematical expressions can combine time-stamps with
        3201             : ! mathematical functions, e.g. to yield the sine of current simulation hour.
        3202             : ! Mathematical expressions must start with the identifier 'MATH:', followed
        3203             : ! by the actual expression. Each expression must include at least one
        3204             : ! variable (evaluated at runtime). The following variables are currently
        3205             : ! supported: YYYY (year), MM (month), DD (day), HH (hour), LH (local hour),
        3206             : ! NN (minute), SS (second), WD (weekday), LWD (local weekday),
        3207             : ! DOY (day of year), ELH (elapsed hours), ELS (elapsed seconds).
        3208             : ! In addition, the following variables can be used: PI (3.141...), DOM
        3209             : ! (\# of days of current month).
        3210             : ! For example, the following expression would yield a continuous sine
        3211             : ! curve as function of hour of day: 'MATH:sin(HH/24*PI*2)'.
        3212             : !\\
        3213             : !\\
        3214             : ! For a full list of valid mathematical expressions, see module interpreter.F90.
        3215             : !\\
        3216             : !\\
        3217             : ! !INTERFACE:
        3218             : !
        3219           0 :   SUBROUTINE ReadMath( HcoState, Lct, ValStr, Vals, N, RC )
        3220             : !
        3221             : ! !USES:
        3222             : !
        3223             :     USE HCO_CLOCK_MOD,      ONLY : HcoClock_Get
        3224             :     USE HCO_tIdx_Mod,       ONLY : HCO_GetPrefTimeAttr
        3225             :     USE INTERPRETER
        3226             : !
        3227             : ! !INPUT PARAMTERS:
        3228             : !
        3229             :     TYPE(HCO_State),  POINTER         :: HcoState    ! HEMCO state
        3230             :     TYPE(ListCont),   POINTER         :: Lct
        3231             :     CHARACTER(LEN=*), INTENT(IN   )   :: ValStr
        3232             : !
        3233             : ! !INPUT/OUTPUT PARAMETERS:
        3234             : !
        3235             :     REAL(hp),         INTENT(INOUT)   :: Vals(:)
        3236             :     INTEGER,          INTENT(INOUT)   :: RC
        3237             : !
        3238             : ! !OUTPUT PARAMETERS:
        3239             : !
        3240             :     INTEGER,          INTENT(  OUT)   :: N
        3241             : !
        3242             : ! !REVISION HISTORY:
        3243             : !  11 May 2017 - C. Keller - Initial version
        3244             : !  See https://github.com/geoschem/hemco for complete history
        3245             : !EOP
        3246             : !------------------------------------------------------------------------------
        3247             : !BOC
        3248             : !
        3249             : ! !LOCAL VARIABLES:
        3250             : !
        3251             :     LOGICAL            :: EOS
        3252             :     INTEGER            :: STRL
        3253             :     INTEGER            :: I, NVAL, LHIDX, LWDIDX
        3254             :     INTEGER            :: prefYr, prefMt, prefDy, prefHr, prefMn
        3255             :     INTEGER            :: prefWD, prefDOY, prefS, LMD, cHr
        3256             :     INTEGER            :: nSteps
        3257             :     REAL(hp)           :: ELH, ELS
        3258             :     REAL(hp)           :: Val
        3259             :     CHARACTER(LEN=255) :: MSG
        3260             :     CHARACTER(LEN=255) :: LOC = 'ReadMath (hcoio_util_mod.F90)'
        3261             :
        3262             :     ! Variables used by the evaluator to build and to determine the value
        3263             :     ! of the expressions
        3264             :     character(len = 10) :: all_variables(12)
        3265             :     real(hp)            :: all_variablesvalues(12)
        3266             :
        3267             :     !String variable that will store the function that the evaluator will build
        3268             :     character (len = 275)  :: func
        3269             :
        3270             :     !String variable that will return the building of the expression result
        3271             :     !If everything was ok then statusflag = 'ok', otherwise statusflag = 'error'
        3272             :     character (len = 5)  :: statusflag
        3273             :
        3274             :     !======================================================================
        3275             :     ! ReadMath begins here
        3276             :     !======================================================================
        3277             :
        3278             :     ! Substring (without flag 'MATH:')
        3279           0 :     STRL = LEN(ValStr)
        3280           0 :     IF ( STRL < 6 ) THEN
        3281             :        MSG = 'Math expression is too short - expected `MATH:<expr>`: ' &
        3282           0 :              //TRIM(ValStr)
        3283           0 :        CALL HCO_ERROR( MSG, RC, THISLOC=LOC )
        3284           0 :        RETURN
        3285             :     ENDIF
        3286           0 :     func = ValStr(6:STRL)
        3287             :
        3288             :     ! Get preferred time stamps
        3289             :     CALL HCO_GetPrefTimeAttr( HcoState, Lct, &
        3290           0 :                               prefYr, prefMt, prefDy, prefHr, prefMn, RC )
        3291           0 :     IF ( RC /= HCO_SUCCESS ) THEN
        3292           0 :         CALL HCO_ERROR( 'ERROR 25', RC, THISLOC=LOC )
        3293           0 :         RETURN
        3294             :     ENDIF
        3295             :
        3296             :     ! Get some other current time stamps
        3297             :     CALL HcoClock_Get( HcoState%Clock,  cS=prefS,     cH=cHr, &
        3298             :                        cWEEKDAY=prefWD, cDOY=prefDOY, LMD=LMD,      &
        3299           0 :                        nSteps=nSteps,   RC=RC )
        3300           0 :     IF ( RC /= HCO_SUCCESS ) THEN
        3301           0 :         CALL HCO_ERROR( 'ERROR 26', RC, THISLOC=LOC )
        3302           0 :         RETURN
        3303             :     ENDIF
        3304             :
        3305             :     ! GetPrefTimeAttr can return -999 for hour. In this case set to current
        3306             :     ! simulation hour
        3307           0 :     IF ( prefHr < 0 ) prefHr = cHr
        3308             :
        3309             :     ! Parse function. This will replace any tokens in the function with the
        3310             :     ! actual token values. (ckeller, 7/7/17)
        3311             :     CALL HCO_CharParse ( HcoState%Config, func, &
        3312           0 :                          prefYr, prefMt, prefDy, prefHr, prefMn, RC )
        3313           0 :     IF ( RC /= HCO_SUCCESS ) THEN
        3314           0 :         CALL HCO_ERROR( 'ERROR 27', RC, THISLOC=LOC )
        3315           0 :         RETURN
        3316             :     ENDIF
        3317             :
        3318             :     ! Elapsed hours and seconds since start time
        3319           0 :     ELS = HcoState%TS_DYN * nSteps
        3320           0 :     ELH = ELS / 3600.0_hp
        3321             :
        3322             :     ! Check which variables are in string.
        3323             :     ! Possible variables are YYYY, MM, DD, WD, HH, NN, SS, DOY, ELH, ELS
        3324           0 :     NVAL   = 0
        3325           0 :     LHIDX  = -1
        3326           0 :     LWDIDX = -1
        3327             :
        3328           0 :     IF ( INDEX(func,'YYYY') > 0 ) THEN
        3329           0 :        NVAL                      = NVAL + 1
        3330           0 :        all_variables(NVAL)       = 'yyyy'
        3331           0 :        all_variablesvalues(NVAL) = prefYr
        3332             :     ENDIF
        3333           0 :     IF ( INDEX(func,'MM') > 0 ) THEN
        3334           0 :        NVAL                      = NVAL + 1
        3335           0 :        all_variables(NVAL)       = 'mm'
        3336           0 :        all_variablesvalues(NVAL) = prefMt
        3337             :     ENDIF
        3338           0 :     IF ( INDEX(func,'DD') > 0 ) THEN
        3339           0 :        NVAL                      = NVAL + 1
        3340           0 :        all_variables(NVAL)       = 'dd'
        3341           0 :        all_variablesvalues(NVAL) = prefDy
        3342             :     ENDIF
        3343           0 :     IF ( INDEX(func,'WD') > 0 ) THEN
        3344           0 :        NVAL                      = NVAL + 1
        3345           0 :        all_variables(NVAL)       = 'wd'
        3346           0 :        all_variablesvalues(NVAL) = prefWD
        3347             :     ENDIF
        3348           0 :     IF ( INDEX(func,'LWD') > 0 ) THEN
        3349           0 :        NVAL                      = NVAL + 1
        3350           0 :        all_variables(NVAL)       = 'lwd'
        3351           0 :        all_variablesvalues(NVAL) = prefWD
        3352           0 :        LWDIDX                    = NVAL
        3353             :     ENDIF
        3354           0 :     IF ( INDEX(func,'HH') > 0 ) THEN
        3355           0 :        NVAL                      = NVAL + 1
        3356           0 :        all_variables(NVAL)       = 'hh'
        3357           0 :        all_variablesvalues(NVAL) = prefHr
        3358             :     ENDIF
        3359           0 :     IF ( INDEX(func,'LH') > 0 ) THEN
        3360           0 :        NVAL                      = NVAL + 1
        3361           0 :        all_variables(NVAL)       = 'lh'
        3362           0 :        all_variablesvalues(NVAL) = prefHr
        3363           0 :        LHIDX                     = NVAL
        3364             :     ENDIF
        3365           0 :     IF ( INDEX(func,'NN') > 0 ) THEN
        3366           0 :        NVAL                      = NVAL + 1
        3367           0 :        all_variables(NVAL)       = 'nn'
        3368           0 :        all_variablesvalues(NVAL) = prefMn
        3369             :     ENDIF
        3370           0 :     IF ( INDEX(func,'SS') > 0 ) THEN
        3371           0 :        NVAL                      = NVAL + 1
        3372           0 :        all_variables(NVAL)       = 'ss'
        3373           0 :        all_variablesvalues(NVAL) = prefS
        3374             :     ENDIF
        3375           0 :     IF ( INDEX(func,'DOY') > 0 ) THEN
        3376           0 :        NVAL                      = NVAL + 1
        3377           0 :        all_variables(NVAL)       = 'doy'
        3378           0 :        all_variablesvalues(NVAL) = prefDOY
        3379             :     ENDIF
        3380           0 :     IF ( INDEX(func,'PI') > 0 ) THEN
        3381           0 :        NVAL                      = NVAL + 1
        3382           0 :        all_variables(NVAL)       = 'pi'
        3383           0 :        all_variablesvalues(NVAL) = HcoState%Phys%PI
        3384             :     ENDIF
        3385           0 :     IF ( INDEX(func,'DOM') > 0 ) THEN
        3386           0 :        NVAL                      = NVAL + 1
        3387           0 :        all_variables(NVAL)       = 'dom'
        3388           0 :        all_variablesvalues(NVAL) = LMD
        3389             :     ENDIF
        3390           0 :     IF ( INDEX(func,'ELH') > 0 ) THEN
        3391           0 :        NVAL                      = NVAL + 1
        3392           0 :        all_variables(NVAL)       = 'elh'
        3393           0 :        all_variablesvalues(NVAL) = ELH
        3394             :     ENDIF
        3395           0 :     IF ( INDEX(func,'ELS') > 0 ) THEN
        3396           0 :        NVAL                      = NVAL + 1
        3397           0 :        all_variables(NVAL)       = 'els'
        3398           0 :        all_variablesvalues(NVAL) = ELS
        3399             :     ENDIF
        3400             :
        3401             :     ! Error trap: cannot have local hour and local weekday in
        3402             :     ! same expression
        3403           0 :     IF ( LHIDX > 0 .AND. LWDIDX > 0 ) THEN
        3404             :        MSG = 'Cannot have local hour and local weekday in '//&
        3405           0 :              'same expression: '//TRIM(func)
        3406           0 :        CALL HCO_ERROR( MSG, RC, THISLOC=LOC )
        3407           0 :        RETURN
        3408             :     ENDIF
        3409             :
        3410             :     ! N is the number of expressions.
        3411           0 :     Vals(:) = -999.0_hp
        3412           0 :     IF ( LHIDX > 0 ) THEN
        3413           0 :        N = 24
        3414           0 :     ELSEIF ( LWDIDX > 0 ) THEN
        3415           0 :        N = 7
        3416             :     ELSE
        3417           0 :        N = 1
        3418             :     ENDIF
        3419             :
        3420             :     ! Evaluate expression
        3421             :     !Initialize function
        3422           0 :     call init (func, all_variables(1:NVAL), statusflag)
        3423           0 :     IF(statusflag == 'ok') THEN
        3424           0 :        DO I=1,N
        3425           0 :           IF ( LHIDX  > 0 ) all_variablesvalues(LHIDX)  = I-1
        3426           0 :           IF ( LWDIDX > 0 ) all_variablesvalues(LWDIDX) = I-1
        3427           0 :           Val = evaluate( all_variablesvalues(1:NVAL) )
        3428           0 :           Vals(I) = Val
        3429             :
        3430             :           ! Verbose
        3431           0 :           IF ( HCO_IsVerb( HcoState%Config%Err ) ) THEN
        3432           0 :              WRITE(MSG,*) 'Evaluated function: ',TRIM(func),' --> ', Val
        3433           0 :              CALL HCO_MSG(HcoState%Config%Err,MSG)
        3434             :           ENDIF
        3435             :        ENDDO
        3436             :     ELSE
        3437           0 :        MSG = 'Error evaluation function: '//TRIM(func)
        3438           0 :        CALL HCO_ERROR( MSG, RC, THISLOC=LOC )
        3439           0 :        RETURN
        3440             :     ENDIF
        3441           0 :     call destroyfunc()
        3442             :
        3443             :     ! Return w/ success
        3444           0 :     RC = HCO_SUCCESS
        3445             :
        3446             :   END SUBROUTINE ReadMath
        3447             : !EOC
        3448             : END MODULE HCOIO_Util_Mod

Generated by: LCOV version 1.14