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

          Line data    Source code
       1             : !BOC
       2             : #if defined ( MODEL_GCCLASSIC ) || defined( MODEL_WRF ) || defined( MODEL_CESM ) || defined( HEMCO_STANDALONE )
       3             : ! The 'standard' HEMCO I/O module is used for:
       4             : ! - HEMCO Standalone (HEMCO_STANDALONE)
       5             : ! - GEOS-Chem 'Classic' (MODEL_GCCLASSIC)
       6             : ! - WRF-GC (MODEL_WRF)
       7             : ! - CESM-GC and CAM-Chem / HEMCO-CESM (MODEL_CESM)
       8             : !EOC
       9             : !------------------------------------------------------------------------------
      10             : !                   Harmonized Emissions Component (HEMCO)                    !
      11             : !------------------------------------------------------------------------------
      12             : !BOP
      13             : !
      14             : ! !MODULE: hcoio_read_std_mod.F90
      15             : !
      16             : ! !DESCRIPTION: Module HCOIO\_read\_mod controls data processing
      17             : ! (file reading, unit conversion, regridding) for HEMCO in the
      18             : ! 'standard' environment (i.e. non-ESMF).
      19             : !
      20             : ! This module implements the 'standard' environment (i.e. non-ESMF).
      21             : !\\
      22             : !\\
      23             : ! !INTERFACE:
      24             : !
      25             : MODULE HCOIO_Read_Mod
      26             : !
      27             : ! !USES:
      28             : !
      29             :   USE HCO_Types_Mod
      30             :   USE HCO_Error_Mod
      31             :   USE HCO_CharTools_Mod
      32             :   USE HCO_State_Mod,       ONLY : Hco_State
      33             :   USE HCOIO_Util_Mod
      34             : 
      35             :   IMPLICIT NONE
      36             :   PRIVATE
      37             : !
      38             : ! !PUBLIC MEMBER FUNCTIONS:
      39             : !
      40             :   PUBLIC  :: HCOIO_Read
      41             :   PUBLIC  :: HCOIO_CloseAll
      42             : !
      43             : ! !REMARKS:
      44             : !  Beginning with HEMCO 3.0.0, all I/O modules use the same module names,
      45             : !  and their compilation depends on pre-processor flags defined at the top
      46             : !  of the file.
      47             : !
      48             : !  This is to streamline the implementation of one unified Data Input Layer,
      49             : !  that can be switched in and out at compile time, and reduce branching of
      50             : !  code paths elsewhere.
      51             : !
      52             : ! !REVISION HISTORY:
      53             : !  22 Aug 2013 - C. Keller   - Initial version
      54             : !  See https://github.com/geoschem/hemco for complete history
      55             : !EOP
      56             : !------------------------------------------------------------------------------
      57             : !BOC
      58             : !
      59             : ! !DEFINED PARAMETERS
      60             : !
      61             :   ! Parameter used for difference testing of floating points
      62             :   REAL(dp), PRIVATE, PARAMETER :: EPSILON = 1.0e-5_dp
      63             : 
      64             : #if defined( MODEL_CESM ) || defined( MODEL_WRF )
      65             :   REAL(hp), PRIVATE            :: GC_72_EDGE_SIGMA(73) = (/ &
      66             :     1.000000E+00, 9.849998E-01, 9.699136E-01, 9.548285E-01, 9.397434E-01, 9.246593E-01, &
      67             :     9.095741E-01, 8.944900E-01, 8.794069E-01, 8.643237E-01, 8.492406E-01, 8.341584E-01, &
      68             :     8.190762E-01, 7.989697E-01, 7.738347E-01, 7.487007E-01, 7.235727E-01, 6.984446E-01, &
      69             :     6.733175E-01, 6.356319E-01, 5.979571E-01, 5.602823E-01, 5.226252E-01, 4.849751E-01, &
      70             :     4.473417E-01, 4.097261E-01, 3.721392E-01, 3.345719E-01, 2.851488E-01, 2.420390E-01, &
      71             :     2.055208E-01, 1.746163E-01, 1.484264E-01, 1.261653E-01, 1.072420E-01, 9.115815E-02, &
      72             :     7.748532E-02, 6.573205E-02, 5.565063E-02, 4.702097E-02, 3.964964E-02, 3.336788E-02, &
      73             :     2.799704E-02, 2.341969E-02, 1.953319E-02, 1.624180E-02, 1.346459E-02, 1.112953E-02, &
      74             :     9.171478E-03, 7.520355E-03, 6.135702E-03, 4.981002E-03, 4.023686E-03, 3.233161E-03, &
      75             :     2.585739E-03, 2.057735E-03, 1.629410E-03, 1.283987E-03, 1.005675E-03, 7.846040E-04, &
      76             :     6.089317E-04, 4.697755E-04, 3.602270E-04, 2.753516E-04, 2.082408E-04, 1.569208E-04, &
      77             :     1.184308E-04, 8.783617E-05, 6.513694E-05, 4.737232E-05, 3.256847E-05, 1.973847E-05, &
      78             :     9.869233E-06/)
      79             : #endif
      80             : 
      81             : CONTAINS
      82             : !EOC
      83             : !------------------------------------------------------------------------------
      84             : !                   Harmonized Emissions Component (HEMCO)                    !
      85             : !------------------------------------------------------------------------------
      86             : !BOP
      87             : !
      88             : ! !IROUTINE: HCOIO_Read
      89             : !
      90             : ! !DESCRIPTION: Reads a netCDF file and returns the regridded array in proper
      91             : ! units. This routine uses the HEMCO generic data reading and regridding
      92             : ! routines.
      93             : !\\
      94             : !\\
      95             : ! Two different regridding algorithm are used: NCREGRID for 3D data with
      96             : ! vertical regridding, and map\_a2a for all other data. map\_a2a also
      97             : ! supports index-based remapping, while this feature is currently not
      98             : ! possible in combination with NCREGRID.
      99             : !\\
     100             : !\\
     101             : ! 3D data is vertically regridded onto the simulation grid on the sigma
     102             : ! interface levels. In order to calculate these levels correctly, the netCDF
     103             : ! vertical coordinate description must adhere to the CF - conventions. See
     104             : ! routine NC\_Get\_Sigma\_Levels in Ncdf\_Mod for more details.
     105             : !\\
     106             : !\\
     107             : ! A simpler vertical interpolation scheme is used if (a) the number of
     108             : ! vertical levels of the input data corresponds to the number of levels
     109             : ! on the simulation grid (direct mapping, no remapping), (b) the vertical
     110             : ! level variable name (long\_name) contains the word "GEOS-Chem level". In
     111             : ! the latter case, the vertical levels of the input data is interpreted as
     112             : ! GEOS vertical levels and mapped onto the simulation grid using routine
     113             : ! ModelLev\_Interpolate.
     114             : !\\
     115             : !\\
     116             : ! !INTERFACE:
     117             : !
     118           0 :   SUBROUTINE HCOIO_Read( HcoState, Lct, RC )
     119             : !
     120             : ! !USES:
     121             : !
     122             :     USE HCO_Ncdf_Mod,       ONLY : NC_Open
     123             :     USE HCO_Ncdf_Mod,       ONLY : NC_Close
     124             :     USE HCO_Ncdf_Mod,       ONLY : NC_Read_Var
     125             :     USE HCO_Ncdf_Mod,       ONLY : NC_Read_Arr
     126             :     USE HCO_Ncdf_Mod,       ONLY : NC_Get_Grid_Edges
     127             :     USE HCO_Ncdf_Mod,       ONLY : NC_Get_Sigma_Levels
     128             :     USE HCO_CHARPAK_MOD,    ONLY : TRANLC
     129             :     USE HCO_Unit_Mod,       ONLY : HCO_Unit_Change
     130             :     USE HCO_Unit_Mod,       ONLY : HCO_Unit_ScalCheck
     131             :     USE HCO_Unit_Mod,       ONLY : HCO_IsUnitless
     132             :     USE HCO_Unit_Mod,       ONLY : HCO_IsIndexData
     133             :     USE HCO_Unit_Mod,       ONLY : HCO_UnitTolerance
     134             :     USE HCO_GeoTools_Mod,   ONLY : HCO_ValidateLon
     135             :     USE HCO_FileData_Mod,   ONLY : FileData_ArrCheck
     136             :     USE HCO_FileData_Mod,   ONLY : FileData_ArrInit
     137             :     USE HCO_FileData_Mod,   ONLY : FileData_Cleanup
     138             :     USE HCOIO_MESSY_MOD,    ONLY : HCO_MESSY_REGRID
     139             :     USE HCO_INTERP_MOD,     ONLY : REGRID_MAPA2A
     140             :     USE HCO_INTERP_MOD,     ONLY : ModelLev_Check
     141             :     USE HCO_CLOCK_MOD,      ONLY : HcoClock_Get
     142             :     USE HCO_DIAGN_MOD,      ONLY : Diagn_Update
     143             :     USE HCO_EXTLIST_MOD,    ONLY : HCO_GetOpt
     144             :     USE HCO_TIDX_MOD,       ONLY : tIDx_IsInRange
     145             : 
     146             :     include "netcdf.inc"
     147             : !
     148             : ! !INPUT PARAMETERS:
     149             : !
     150             :     TYPE(HCO_State),  POINTER        :: HcoState   ! HEMCO state object
     151             :     TYPE(ListCont),   POINTER        :: Lct        ! HEMCO list container
     152             : !
     153             : ! !INPUT/OUTPUT PARAMETERS:
     154             : !
     155             :     INTEGER,          INTENT(INOUT)  :: RC         ! Success or failure?
     156             : !
     157             : ! !REVISION HISTORY:
     158             : !  13 Mar 2013 - C. Keller   - Initial version
     159             : !  See https://github.com/geoschem/hemco for complete history
     160             : !EOP
     161             : !------------------------------------------------------------------------------
     162             : !BOC
     163             : !
     164             : ! !LOCAL VARIABLES:
     165             : !
     166             :     CHARACTER(LEN=255)            :: thisUnit, LevUnit, LevName
     167             :     CHARACTER(LEN=1023)           :: MSG, LOC
     168             :     CHARACTER(LEN=1023)           :: srcFile, srcFile2
     169             :     INTEGER                       :: NX, NY
     170             :     INTEGER                       :: NCRC, Flag, AS
     171             :     INTEGER                       :: ncLun, ncLun2
     172             :     INTEGER                       :: ierr,   v_id
     173             :     INTEGER                       :: nlon,   nlat,  nlev, nTime
     174             :     INTEGER                       :: lev1,   lev2,  dir
     175             :     INTEGER                       :: tidx1,  tidx2,  ncYr,  ncMt
     176             :     INTEGER                       :: tidx1b, tidx2b, ncYr2, ncMt2
     177             :     INTEGER                       :: HcoID
     178             :     INTEGER                       :: ArbIdx
     179             :     INTEGER                       :: nlatEdge, nlonEdge
     180             :     INTEGER                       :: Direction
     181             :     REAL(hp)                      :: MW_g
     182             :     REAL(sp)                      :: wgt1,   wgt2
     183           0 :     REAL(sp), POINTER             :: ncArr(:,:,:,:)
     184           0 :     REAL(sp), POINTER             :: ncArr2(:,:,:,:)
     185           0 :     REAL(hp), POINTER             :: SigEdge(:,:,:)
     186           0 :     REAL(hp), POINTER             :: SigLev (:,:,:)
     187           0 :     REAL(hp), POINTER             :: LonMid   (:)
     188           0 :     REAL(hp), POINTER             :: LatMid   (:)
     189           0 :     REAL(hp), POINTER             :: LevMid   (:)
     190           0 :     REAL(hp), POINTER             :: LonEdge  (:)
     191           0 :     REAL(hp), POINTER             :: LatEdge  (:)
     192             :     REAL(hp)                      :: UnitFactor
     193             :     LOGICAL                       :: KeepSpec
     194             :     LOGICAL                       :: FOUND
     195             :     LOGICAL                       :: IsModelLevel
     196             :     LOGICAL                       :: DoReturn
     197             :     INTEGER                       :: UnitTolerance
     198             :     INTEGER                       :: AreaFlag, TimeFlag
     199             :     REAL(dp)                      :: YMDhma,  YMDhmb, YMDhm1
     200             :     REAL(dp)                      :: oYMDhm1, oYMDhm2
     201             :     INTEGER                       :: cYr, cMt, cDy, cHr, Yr1, Yr2
     202             :     INTEGER                       :: nYears, iYear
     203             :     INTEGER                       :: I, J
     204             : 
     205             :     ! Use MESSy regridding routines?
     206             :     LOGICAL                       :: UseMESSy
     207             : 
     208             :     ! SAVEd scalars
     209             :     LOGICAL, SAVE                 :: doPrintWarning = .TRUE.
     210             : 
     211             :     !=================================================================
     212             :     ! HCOIO_READ begins here
     213             :     !=================================================================
     214           0 :     LOC = 'HCOIO_READ (HCOIO_READ_STD_MOD.F90)'
     215             : 
     216             :     ! Do not try to read a mask file where the mask bounding box limits
     217             :     ! are given in the srcFile location, as there is no file to read.
     218             :     ! This should fix https://github.com/geoschem/HEMCO/issues/153.
     219             :     !   -- Bob Yantosca (12 Jul 2022)
     220           0 :     IF ( Lct%Dct%DctType == HCO_DCTTYPE_MASK ) THEN
     221           0 :        IF ( .not. Lct%Dct%Dta%ncRead ) RETURN
     222             :     ENDIF
     223             : 
     224             :     ! Enter
     225           0 :     CALL HCO_ENTER( HcoState%Config%Err, LOC, RC )
     226           0 :     IF ( RC /= HCO_SUCCESS ) THEN
     227           0 :         CALL HCO_ERROR( 'ERROR 0', RC, THISLOC=LOC )
     228           0 :         RETURN
     229             :     ENDIF
     230             : 
     231             :     ! Initialize pointers
     232           0 :     ncArr   => NULL()
     233           0 :     ncArr2  => NULL()
     234           0 :     SigEdge => NULL()
     235           0 :     SigLev  => NULL()
     236           0 :     LonMid  => NULL()
     237           0 :     LatMid  => NULL()
     238           0 :     LevMid  => NULL()
     239           0 :     LonEdge => NULL()
     240           0 :     LatEdge => NULL()
     241             : 
     242             :     ! Zero local variables for safety's sake
     243           0 :     dir     =  0
     244           0 :     lev1    =  0
     245           0 :     lev2    =  0
     246           0 :     ncYr    =  0
     247           0 :     ncMt    =  0
     248           0 :     ncYr2   =  0
     249           0 :     ncMt2   =  0
     250           0 :     nLon    =  0
     251           0 :     nLat    =  0
     252           0 :     nLev    =  0
     253           0 :     nTime   =  0
     254           0 :     tIdx1   =  0
     255           0 :     tIdx2   =  0
     256           0 :     tidx1b  =  0
     257           0 :     tidx2b  =  0
     258           0 :     wgt1    =  0.0_sp
     259           0 :     wgt2    =  0.0_sp
     260             : 
     261             :     ! Get unit tolerance set in configuration file
     262           0 :     UnitTolerance = HCO_UnitTolerance( HcoState%Config )
     263             : 
     264             :     ! For convenience, copy horizontal grid dimensions from HEMCO
     265             :     ! state object
     266           0 :     NX = HcoState%NX
     267           0 :     NY = HcoState%NY
     268             : 
     269             :     ! Verbose
     270           0 :     IF ( HCO_IsVerb( HcoState%Config%Err ) ) THEN
     271           0 :        WRITE(MSG,*) 'Processing container: ', TRIM(Lct%Dct%cName)
     272           0 :        CALL HCO_MSG( HcoState%Config%Err, MSG, SEP1='-' )
     273             :     ENDIF
     274             : 
     275             :     ! If the file has cycle flag "E" (e.g. it's a restart file), then we will
     276             :     ! read it only once and then never again.  If the file has already been
     277             :     ! read on a previous call, then don't call HCOIO_READ. (bmy, 10/4/18)
     278             :     !
     279             :     ! Moved this handling from hcoio_dataread_mod, as it is non-MAPL specific
     280             :     ! (hplin, 4/5/21)
     281             :     IF ( Lct%Dct%Dta%CycleFlag == HCO_CFLAG_EXACT .and.                      &
     282           0 :          Lct%Dct%Dta%UpdtFlag  == HCO_UFLAG_ONCE  .and.                      &
     283             :          Lct%Dct%Dta%isTouched                          ) THEN
     284             : 
     285             :        ! Print a warning message only once
     286           0 :        IF ( doPrintWarning ) THEN
     287           0 :           doPrintWarning = .FALSE.
     288             :           MSG = 'No further attempts will be made to read file: ' //         &
     289           0 :                 TRIM( Lct%Dct%Dta%NcFile )
     290           0 :           CALL HCO_WARNING ( HcoState%Config%Err, MSG, RC )
     291             :        ENDIF
     292             : 
     293             :        ! Return without reading
     294           0 :        CALL HCO_LEAVE( HcoState%Config%Err, RC )
     295           0 :        RETURN
     296             :     ENDIF
     297             : 
     298             :     ! ----------------------------------------------------------------
     299             :     ! Parse source file name. This will replace all tokens ($ROOT,
     300             :     ! ($YYYY), etc., with valid values.
     301             :     ! ----------------------------------------------------------------
     302           0 :     CALL SrcFile_Parse( HcoState, Lct, srcFile, FOUND, RC )
     303           0 :     IF ( RC /= HCO_SUCCESS ) THEN
     304             :        MSG = 'Error encountered in routine "SrcFile_Parse", located '     // &
     305           0 :              'module src/Core/hcoio_read_std_mod.F90!'
     306           0 :         CALL HCO_ERROR( MSG, RC )
     307           0 :         RETURN
     308             :     ENDIF
     309             : 
     310             :     ! Handle found or not in the standard way if HEMCO is in regular run mode.
     311           0 :     IF ( .NOT. HcoState%Options%isDryRun ) THEN
     312             : 
     313             :        !====================================================================
     314             :        ! HEMCO is in regular simulation mode (not dry-run)!
     315             :        !====================================================================
     316             : 
     317             :        ! If file not found, return w/ error. No error if cycling attribute is
     318             :        ! select to range. In that case, just make sure that array is empty.
     319           0 :        IF ( .NOT. FOUND ) THEN
     320           0 :           IF ( ( Lct%Dct%Dta%CycleFlag == HCO_CFLAG_RANGE ) .OR.      &
     321             :                ( Lct%Dct%Dta%CycleFlag == HCO_CFLAG_EXACT )     ) THEN
     322             : 
     323             :              ! If MustFind flag is enabled, return with error if field is not
     324             :              ! found
     325           0 :              IF ( Lct%Dct%Dta%MustFind ) THEN
     326             :                 MSG = 'Cannot find file for current simulation time: ' // &
     327             :                      TRIM(srcFile) // ' - Cannot get field ' // &
     328             :                      TRIM(Lct%Dct%cName) // '. Please check file name ' // &
     329           0 :                      'and time (incl. time range flag) in the config. file'
     330           0 :                 CALL HCO_ERROR( MSG, RC )
     331           0 :                 RETURN
     332             : 
     333             :              ! If MustFind flag is not enabled, ignore this field and return
     334             :              ! with a warning.
     335             :              ELSE
     336           0 :                 CALL FileData_Cleanup( Lct%Dct%Dta, DeepClean=.FALSE. )
     337             :                 MSG = 'No valid file found for current simulation time - data '// &
     338           0 :                      'will be ignored for time being - ' // TRIM(Lct%Dct%cName)
     339           0 :                 CALL HCO_WARNING ( HcoState%Config%Err, MSG, RC )
     340           0 :                 CALL HCO_LEAVE ( HcoState%Config%Err,  RC )
     341           0 :                 RETURN
     342             :              ENDIF
     343             : 
     344             :           ELSE
     345             :              MSG = 'Cannot find file for current simulation time: ' // &
     346             :                   TRIM(srcFile) // ' - Cannot get field ' // &
     347             :                   TRIM(Lct%Dct%cName) // '. Please check file name ' // &
     348           0 :                   'and time (incl. time range flag) in the config. file'
     349           0 :              CALL HCO_ERROR( MSG, RC )
     350           0 :              RETURN
     351             :           ENDIF
     352             :        ENDIF
     353             : 
     354             :     ELSE
     355             : 
     356             :        !====================================================================
     357             :        ! HEMCO is in a "dry-run" mode!
     358             :        !====================================================================
     359             : 
     360             :        ! Simulate file read buffer
     361           0 :        IF ( TRIM(HcoState%ReadLists%FileInArchive) == TRIM(srcFile) ) THEN
     362           0 :           CALL HCO_LEAVE ( HcoState%Config%Err,  RC )
     363           0 :           RETURN
     364             :        ENDIF
     365             : 
     366             :        ! If file exists, print the result. If NOT, then handle accordingly
     367             :        ! But NEVER error out (HCO_ERROR), as we want to get a list of all
     368             :        ! files. (hplin, 11/2/19)
     369           0 :        IF ( .NOT. FOUND ) THEN
     370           0 :           IF ( ( Lct%Dct%Dta%CycleFlag == HCO_CFLAG_RANGE )       .OR.       &
     371             :                ( Lct%Dct%Dta%CycleFlag == HCO_CFLAG_EXACT )     ) THEN
     372             : 
     373             :              ! If MustFind flag is enabled, return with error if field is not
     374             :              ! found
     375           0 :              IF ( Lct%Dct%Dta%MustFind ) THEN
     376             :                 MSG = 'Cannot find file for current simulation time: '    // &
     377             :                      TRIM(srcFile) // ' - Cannot get field '              // &
     378             :                      TRIM(Lct%Dct%cName) // '. Please check file name '   // &
     379           0 :                      'and time (incl. time range flag) in the config. file'
     380           0 :                 CALL HCO_Warning( HcoState%Config%Err, MSG, RC )
     381             : 
     382             :                 ! Write a msg to stdout (NOT FOUND)
     383           0 :                 WRITE( 6, 300 ) TRIM( srcFile )
     384             :  300            FORMAT( 'HEMCO: REQUIRED FILE NOT FOUND ', a )
     385             : 
     386             :              ! If MustFind flag is not enabled, ignore this field and return
     387             :              ! with a warning.
     388             :              ELSE
     389           0 :                 CALL FileData_Cleanup( Lct%Dct%Dta, DeepClean=.FALSE. )
     390             :                 MSG = 'No valid file found for current simulation time - '// &
     391             :                      'data will be ignored for time being - '             // &
     392           0 :                      TRIM(Lct%Dct%cName)
     393           0 :                 CALL HCO_WARNING ( HcoState%Config%Err, MSG, RC )
     394             : 
     395             :                 ! Write a msg to stdout (OPTIONAL)
     396           0 :                 WRITE( 6, 310 ) TRIM( srcFile )
     397             :  310            FORMAT( 'HEMCO: OPTIONAL FILE NOT FOUND ', a )
     398             : 
     399             :              ENDIF
     400             : 
     401             :           ! Not range or exact
     402             :           ELSE
     403             :              MSG = 'Cannot find file for current simulation time: '       // &
     404             :                   TRIM(srcFile) // ' - Cannot get field '                 // &
     405             :                   TRIM(Lct%Dct%cName) // '. Please check file name '      // &
     406           0 :                   'and time (incl. time range flag) in the config. file'
     407           0 :              CALL HCO_WARNING ( HcoState%Config%Err, MSG, RC )
     408             : 
     409             :              ! Write a msg to stdout (NOT FOUND)
     410           0 :              WRITE( 6, 300 ) TRIM(srcFile)
     411             : 
     412             :           ENDIF
     413             :        ELSE
     414             : 
     415             :           ! Write a mesage to stdout (HEMCO: Opening...)
     416           0 :           WRITE( 6, 100 ) TRIM( srcFile )
     417             : 
     418             :        ENDIF
     419             : 
     420             :        ! It is safe to leave now, we do not need to handle opening the file.
     421             :        ! This may be changed in the future if the "dry-run" mode requires
     422             :        ! a check of the file contents... this may be beyond the scope for now.
     423             : 
     424             :        ! Simulate the "reading" in netCDF to prevent duplicate entries
     425             :        ! in the log
     426           0 :        HcoState%ReadLists%FileInArchive = TRIM(srcFile)
     427             : 
     428             :        ! Skip further processing
     429           0 :        CALL HCO_LEAVE ( HcoState%Config%Err,  RC )
     430           0 :        RETURN
     431             :     ENDIF ! End of dry-run mode else clause
     432             : 
     433             :     ! ----------------------------------------------------------------
     434             :     ! Open netCDF
     435             :     ! ----------------------------------------------------------------
     436             : 
     437             :     ! Check if file is already in buffer. In that case use existing
     438             :     ! open stream. Otherwise open new file. At any given time there
     439             :     ! can only be one file in buffer.
     440           0 :     ncLun = -1
     441           0 :     IF ( HcoState%ReadLists%FileLun > 0 ) THEN
     442           0 :        IF ( TRIM(HcoState%ReadLists%FileInArchive) == TRIM(srcFile) ) THEN
     443           0 :           ncLun = HcoState%ReadLists%FileLun
     444             :        ELSE
     445           0 :           CALL NC_CLOSE ( HcoState%ReadLists%FileLun )
     446           0 :           HcoState%ReadLists%FileLun = -1
     447             :        ENDIF
     448             :     ENDIF
     449             : 
     450             :     ! To read from existing stream:
     451           0 :     IF ( ncLun > 0 ) THEN
     452             : 
     453             :        ! Verbose mode
     454           0 :        IF ( HCO_IsVerb( HcoState%Config%Err ) ) THEN
     455           0 :           WRITE(MSG,*) 'Reading from existing stream: ', TRIM(srcFile)
     456           0 :           CALL HCO_MSG( HcoState%Config%Err, MSG )
     457             :        ENDIF
     458             : 
     459             :     ! To open a new file:
     460             :     ELSE
     461           0 :        CALL NC_OPEN ( TRIM(srcFile), ncLun )
     462             : 
     463             :        ! Verbose mode
     464           0 :        IF ( HCO_IsVerb( HcoState%Config%Err ) ) THEN
     465           0 :           WRITE(MSG,*) 'Opening file: ', TRIM(srcFile)
     466           0 :           CALL HCO_MSG( HcoState%Config%Err, MSG )
     467             :        ENDIF
     468             : 
     469             :        ! Also write to standard output
     470           0 :        WRITE( 6, 100 ) TRIM( srcFile )
     471             :  100   FORMAT( 'HEMCO: Opening ', a )
     472             : 
     473             :        ! This is now the file in archive
     474           0 :        HcoState%ReadLists%FileInArchive = TRIM(srcFile)
     475           0 :        HcoState%ReadLists%FileLun       = ncLun
     476             :     ENDIF
     477             : 
     478             :     ! ----------------------------------------------------------------
     479             :     ! Extract time slice information
     480             :     ! This determines the lower and upper time slice index (tidx1
     481             :     ! and tidx2) to be read based upon the time slice information
     482             :     ! extracted from the file and the time stamp settings set in the
     483             :     ! HEMCO configuration file. Multiple time slices are only selected
     484             :     ! for weekdaily data or for 'autodetected' hourly data (using the
     485             :     ! wildcard character in the configuration file time attribute) or
     486             :     ! if data shall be interpolated between two (consecutive) time
     487             :     ! slices. The weights to be assigned to those two time slices is
     488             :     ! also calculated in GET_TIMEIDX and returned as variables wgt1
     489             :     ! and wgt2, respectively.
     490             :     ! ----------------------------------------------------------------
     491             :     CALL GET_TIMEIDX ( HcoState,  Lct,                &
     492             :                        ncLun,     tidx1,    tidx2,    &
     493             :                        wgt1,      wgt2,     oYMDhm1,  &
     494           0 :                        YMDhma,    YMDhm1,   RC        )
     495           0 :     IF ( RC /= HCO_SUCCESS ) THEN
     496           0 :         CALL HCO_ERROR( 'ERROR 1', RC, THISLOC=LOC )
     497           0 :         RETURN
     498             :     ENDIF
     499             : 
     500             :     !-----------------------------------------------------------------
     501             :     ! Check for negative tidx1. tidx1 can still be negative if:
     502             :     ! (a) CycleFlag is set to range and the current simulation
     503             :     ! time is outside of the data time range. In this case, we
     504             :     ! prompt a warning and make sure that there is no data
     505             :     ! associated with this FileData container.
     506             :     ! (b) CycleFlag is set to exact and none of the data time
     507             :     ! stamps matches the current simulation time exactly. Return
     508             :     ! with error!
     509             :     !-----------------------------------------------------------------
     510           0 :     IF ( tidx1 < 0 ) THEN
     511           0 :        DoReturn = .FALSE.
     512           0 :        IF ( Lct%Dct%Dta%CycleFlag == HCO_CFLAG_CYCLE ) THEN
     513           0 :           MSG = 'Invalid time index in ' // TRIM(srcFile)
     514           0 :           CALL HCO_ERROR( MSG, RC )
     515             :           DoReturn = .TRUE.
     516           0 :        ELSEIF ( ( Lct%Dct%Dta%CycleFlag == HCO_CFLAG_RANGE ) .OR.      &
     517             :                 ( Lct%Dct%Dta%CycleFlag == HCO_CFLAG_EXACT )     ) THEN
     518           0 :           IF ( Lct%Dct%Dta%MustFind ) THEN
     519             :              MSG = 'Cannot find field with valid time stamp in ' // &
     520             :                    TRIM(srcFile) // ' - Cannot get field ' // &
     521             :                    TRIM(Lct%Dct%cName) // '. Please check file name ' // &
     522           0 :                    'and time (incl. time range flag) in the config. file'
     523           0 :              CALL HCO_ERROR( MSG, RC )
     524             :              DoReturn = .TRUE.
     525             :           ELSE
     526           0 :              CALL FileData_Cleanup( Lct%Dct%Dta, DeepClean=.FALSE.)
     527             :              MSG = 'Simulation time is outside of time range provided for '//&
     528           0 :                   TRIM(Lct%Dct%cName) // ' - field is ignored for the time being!'
     529           0 :              CALL HCO_WARNING ( HcoState%Config%Err, MSG, RC )
     530           0 :              DoReturn = .TRUE.
     531           0 :              CALL HCO_LEAVE ( HcoState%Config%Err,  RC )
     532             :           ENDIF
     533             :        ENDIF
     534             : 
     535             :        ! Eventually return here
     536             :        IF ( DoReturn ) THEN
     537           0 :           RETURN
     538             :        ENDIF
     539             :     ENDIF
     540             : 
     541             :     ! ----------------------------------------------------------------
     542             :     ! Check if variable is in file
     543             :     ! ----------------------------------------------------------------
     544           0 :     ierr = Nf_Inq_Varid( ncLun, Lct%Dct%Dta%ncPara, v_id )
     545           0 :     IF ( ierr /= NF_NOERR ) THEN
     546             : 
     547             :        ! If MustFind flag is enabled, return with error if field is not
     548             :        ! found
     549           0 :        IF ( Lct%Dct%Dta%MustFind ) THEN
     550             :           MSG = 'Cannot find field ' // TRIM(Lct%Dct%cName) // &
     551           0 :                 '. Please check variable name in the config. file'
     552           0 :           CALL HCO_ERROR( MSG, RC )
     553           0 :           RETURN
     554             : 
     555             :        ! If MustFind flag is not enabled, ignore this field and return
     556             :        ! with a warning.
     557             :        ELSE
     558           0 :           CALL FileData_Cleanup( Lct%Dct%Dta, DeepClean=.FALSE. )
     559             :           MSG = 'Cannot find field ' // TRIM(Lct%Dct%cName) // &
     560           0 :                 '. Will be ignored for time being.'
     561           0 :           CALL HCO_WARNING ( HcoState%Config%Err, MSG, RC )
     562           0 :           CALL HCO_LEAVE ( HcoState%Config%Err,  RC )
     563           0 :           RETURN
     564             :        ENDIF
     565             :     ENDIF
     566             : 
     567             :     ! ----------------------------------------------------------------
     568             :     ! Read grid
     569             :     ! ----------------------------------------------------------------
     570             : 
     571             :     ! Extract longitude midpoints
     572           0 :     CALL NC_READ_VAR ( ncLun, 'lon', nlon, thisUnit, LonMid, NCRC )
     573           0 :     IF ( NCRC /= 0 ) THEN
     574           0 :        CALL HCO_ERROR( 'NC_READ_VAR: lon', RC )
     575           0 :        RETURN
     576             :     ENDIF
     577             : 
     578           0 :     IF ( nlon == 0 ) THEN
     579           0 :        CALL NC_READ_VAR ( ncLun, 'longitude', nlon, thisUnit, LonMid, NCRC )
     580             :     ENDIF
     581           0 :     IF ( NCRC /= 0 ) THEN
     582           0 :        CALL HCO_ERROR( 'NC_READ_VAR: longitude', RC )
     583           0 :        RETURN
     584             :     ENDIF
     585             : 
     586           0 :     IF ( nlon == 0 ) THEN
     587           0 :        CALL NC_READ_VAR ( ncLun, 'Longitude', nlon, thisUnit, LonMid, NCRC )
     588             :     ENDIF
     589           0 :     IF ( NCRC /= 0 ) THEN
     590           0 :        CALL HCO_ERROR( 'NC_READ_LON: Longitude', RC )
     591           0 :        RETURN
     592             :     ENDIF
     593             : 
     594           0 :     IF ( nlon == 0 ) THEN
     595             :        MSG = 'Cannot find longitude variable in ' // TRIM(srcFile) // &
     596           0 :              ' - Must be one of `lon`, `longitude`, `Longitude`'
     597           0 :        CALL HCO_ERROR( MSG, RC )
     598           0 :        RETURN
     599             :     ENDIF
     600             : 
     601             :     ! Unit must be degrees_east
     602           0 :     CALL TRANLC( thisUnit)
     603           0 :     IF ( INDEX( thisUnit, 'degrees_east' ) == 0 ) THEN
     604             :        MSG = 'illegal longitude unit in ' // TRIM(srcFile) // &
     605           0 :              ' - Must be `degrees_east`.'
     606           0 :        CALL HCO_ERROR( MSG, RC )
     607           0 :        RETURN
     608             :     ENDIF
     609             : 
     610             :     ! Make sure longitude is steadily increasing.
     611           0 :     CALL HCO_ValidateLon( HcoState, nlon, LonMid, RC )
     612           0 :     IF ( RC /= HCO_SUCCESS ) THEN
     613           0 :         CALL HCO_ERROR( 'ERROR 2', RC, THISLOC=LOC )
     614           0 :         RETURN
     615             :     ENDIF
     616             : 
     617             :     ! Extract latitude midpoints
     618           0 :     CALL NC_READ_VAR ( ncLun, 'lat', nlat, thisUnit, LatMid, NCRC )
     619           0 :     IF ( NCRC /= 0 ) THEN
     620           0 :        CALL HCO_ERROR( 'NC_READ_LON: lat', RC )
     621           0 :        RETURN
     622             :     ENDIF
     623             : 
     624           0 :     IF ( nlat == 0 ) THEN
     625           0 :        CALL NC_READ_VAR ( ncLun, 'latitude', nlat, thisUnit, LatMid, NCRC )
     626             :     ENDIF
     627           0 :     IF ( NCRC /= 0 ) THEN
     628           0 :        CALL HCO_ERROR( 'NC_READ_LON: latitude', RC )
     629           0 :        RETURN
     630             :     ENDIF
     631             : 
     632           0 :     IF ( nlat == 0 ) THEN
     633           0 :        CALL NC_READ_VAR ( ncLun, 'Latitude', nlat, thisUnit, LatMid, NCRC )
     634             :     ENDIF
     635           0 :     IF ( NCRC /= 0 ) THEN
     636           0 :        CALL HCO_ERROR( 'NC_READ_LON: Latitude', RC )
     637           0 :        RETURN
     638             :     ENDIF
     639             : 
     640           0 :     IF ( nlat == 0 ) THEN
     641             :        MSG = 'Cannot find latitude variable in ' // TRIM(srcFile) // &
     642           0 :              ' - Must be one of `lat`, `latitude`, `Latitude`'
     643           0 :        CALL HCO_ERROR( MSG, RC )
     644           0 :        RETURN
     645             :     ENDIF
     646             : 
     647             :     ! Unit must be degrees_north
     648           0 :     CALL TRANLC( thisUnit)
     649           0 :     IF ( INDEX( thisUnit, 'degrees_north' ) == 0 ) THEN
     650             :        MSG = 'illegal latitude unit in ' // TRIM(srcFile) // &
     651           0 :              ' - Must be `degrees_north`.'
     652           0 :        CALL HCO_ERROR( MSG, RC )
     653           0 :        RETURN
     654             :     ENDIF
     655             : 
     656             :     ! Get level index if we are dealing with 3D data
     657           0 :     IF ( Lct%Dct%Dta%SpaceDim == 3 ) THEN
     658             : 
     659             :        ! Try to extract level midpoints
     660           0 :        LevName = 'lev'
     661           0 :        CALL NC_READ_VAR ( ncLun, LevName, nlev, LevUnit, LevMid, NCRC )
     662           0 :        IF ( NCRC /= 0 ) THEN
     663           0 :           CALL HCO_ERROR( 'NC_READ_VAR: lev', RC )
     664           0 :           RETURN
     665             :        ENDIF
     666           0 :        IF ( nlev == 0 ) THEN
     667           0 :           LevName = 'height'
     668           0 :           CALL NC_READ_VAR ( ncLun, LevName, nlev, LevUnit, LevMid, NCRC )
     669           0 :           IF ( NCRC /= 0 ) THEN
     670           0 :              CALL HCO_ERROR( 'NC_READ_VAR: height', RC )
     671           0 :              RETURN
     672             :           ENDIF
     673             :        ENDIF
     674           0 :        IF ( nlev == 0 ) THEN
     675           0 :           LevName = 'level'
     676           0 :           CALL NC_READ_VAR ( ncLun, LevName, nlev, LevUnit, LevMid, NCRC )
     677           0 :           IF ( NCRC /= 0 ) THEN
     678           0 :              CALL HCO_ERROR( 'NC_READ_VAR: level', RC )
     679           0 :              RETURN
     680             :           ENDIF
     681             :        ENDIF
     682             : 
     683             :        ! Error check
     684           0 :        IF ( nlev == 0 ) THEN
     685             :           MSG = 'Cannot find vertical coordinate variable in ' // &
     686           0 :                  TRIM(SrcFile) // ' - Must be one of `lev`, `level`, `height`.'
     687           0 :           CALL HCO_ERROR( MSG, RC )
     688           0 :           RETURN
     689             :        ENDIF
     690             : 
     691             :        ! Are these model levels? This will only return true if 
     692             :        ! (1) the variable is on 72/73 levels and you are going to 47
     693             :        ! levels, (2) if you are on 102/103 levels and you are going
     694             :        ! to 74 levels, (3) if you are on 47/48 levels and you are
     695             :        ! going to 72 levels. Otherwise, use MESSy (nbalasus, 8/24/2023).
     696           0 :        IF ( Lct%Dct%Dta%Levels == 0 ) THEN
     697             : 
     698             : #if defined( MODEL_CESM ) || defined( MODEL_WRF )
     699             : 
     700             :           ! In WRF/CESM, IsModelLevel has a different meaning of "GEOS-Chem levels"
     701             :           ! because the models in WRF and CESM are user-defined and thus fixed input
     702             :           ! files would never be on the model level. In this case, a check is added
     703             :           ! in order to match the file with known GEOS-Chem levels, and if so, the
     704             :           ! data will be handled later in this file accordingly to be vertically
     705             :           ! regridded to the runtime model levels using MESSy.
     706             :           ! This fixes a regression from the vertical regridding fixes in 3.7.1.
     707             :           !
     708             :           ! The meaning of "is model levels" in WRF and CESM are different.
     709             :           ! Model levels can be changed and thus data is never on the model level.
     710             :           ! In this case, IsModelLevel means that the data is on standard
     711             :           ! GEOS-Chem levels, and if so, the data will be handled accordingly
     712             :           ! using a hard-coded set of GEOS-Chem levels to be interpolated using MESSy.
     713             :           ! (hplin, 10/15/23)
     714           0 :           IF ( TRIM(LevUnit) == "level" .or. TRIM(LevUnit) == "GEOS-Chem level" ) THEN
     715             :           ! the below check will be obsolete and is unmaintainable, but would be consistent with ModelLev_Check.
     716             :           ! it is more robust to check for the explicit intention of LevUnit
     717             :           !     nlev == 47 .or. nlev == 48 .or. nlev == 36 .or. nlev == 72 .or. nlev == 73 ) THEN
     718           0 :               IsModelLevel = .true.
     719             :           ENDIF
     720             : 
     721             : #else
     722             : 
     723             :           CALL ModelLev_Check( HcoState, nlev, IsModelLevel, RC )
     724             :           IF ( RC /= HCO_SUCCESS ) THEN
     725             :               CALL HCO_ERROR( 'ERROR 3', RC, THISLOC=LOC )
     726             :               RETURN
     727             :           ENDIF
     728             : 
     729             : #endif
     730             : 
     731             :           ! Set level indeces to be read
     732           0 :           lev1 = 1
     733           0 :           lev2 = nlev
     734             : 
     735             :        ! If levels are explicitly given:
     736             :        ELSE
     737             : 
     738             :           ! Number of levels to be read must be smaller or equal to total
     739             :           ! number of available levels
     740           0 :           IF ( ABS(Lct%Dct%Dta%Levels) > nlev ) THEN
     741           0 :              WRITE(MSG,*) Lct%Dct%Dta%Levels, ' levels requested but file ', &
     742           0 :                 'has only ', nlev, ' levels: ', TRIM(Lct%Dct%cName)
     743           0 :              CALL HCO_ERROR( MSG, RC )
     744           0 :              RETURN
     745             :           ENDIF
     746             : 
     747             :           ! Set levels to be read
     748           0 :           IF ( Lct%Dct%Dta%Levels > 0 ) THEN
     749           0 :              lev1 = 1
     750           0 :              lev2 = Lct%Dct%Dta%Levels
     751             : 
     752             :           ! Reverse axis!
     753             :           ELSE
     754           0 :              lev1 = nlev
     755           0 :              lev2 = nlev + Lct%Dct%Dta%Levels + 1
     756             :           ENDIF
     757             : 
     758             :           ! Use MESSy regridding
     759           0 :           IsModelLevel = .FALSE.
     760             : 
     761             :        ENDIF
     762             : 
     763             :        ! Verbose
     764           0 :        IF ( HCO_IsVerb( HcoState%Config%Err ) ) THEN
     765           0 :           WRITE(MSG,*) 'Will read vertical levels ', lev1, ' to ', lev2
     766           0 :           CALL HCO_MSG(HcoState%Config%Err,MSG)
     767             :        ENDIF
     768             : 
     769           0 :        IF ( HCO_IsVerb( HcoState%Config%Err ) .AND. IsModelLevel ) THEN
     770           0 :           WRITE(MSG,*) 'Data is assumed to already be on the model level grid'
     771           0 :           CALL HCO_MSG(HcoState%Config%Err,MSG)
     772             :        ENDIF
     773             : 
     774             :     ! For 2D data, set lev1 and lev2 to zero. This will ignore
     775             :     ! the level dimension in the netCDF reading call that follows.
     776             :     ELSE
     777           0 :        nlev        = 0
     778           0 :        lev1        = 0
     779           0 :        lev2        = 0
     780           0 :        IsModelLevel = .FALSE.
     781             :     ENDIF
     782             : 
     783             :     ! ----------------------------------------------------------------
     784             :     ! Check for arbitrary additional dimension. Will return -1 if not
     785             :     ! set.
     786             :     ! ----------------------------------------------------------------
     787           0 :     CALL GetArbDimIndex( HcoState, ncLun, Lct, ArbIdx, RC )
     788           0 :     IF ( RC /= HCO_SUCCESS ) THEN
     789           0 :         CALL HCO_ERROR( 'ERROR 4', RC, THISLOC=LOC )
     790           0 :         RETURN
     791             :     ENDIF
     792             : 
     793             :     ! ----------------------------------------------------------------
     794             :     ! Read data
     795             :     ! ----------------------------------------------------------------
     796             : 
     797             :     ! Verbose mode
     798           0 :     IF ( HCO_IsVerb( HcoState%Config%Err ) ) THEN
     799           0 :        WRITE(MSG,*) 'Reading variable ', TRIM(Lct%Dct%Dta%ncPara)
     800           0 :        CALL HCO_MSG(HcoState%Config%Err,MSG)
     801             :     ENDIF
     802             : 
     803             :     CALL NC_READ_ARR( fID     = ncLun,              &
     804             :                       ncVar   = Lct%Dct%Dta%ncPara, &
     805             :                       lon1    = 1,                  &
     806             :                       lon2    = nlon,               &
     807             :                       lat1    = 1,                  &
     808             :                       lat2    = nlat,               &
     809             :                       lev1    = lev1,               &
     810             :                       lev2    = lev2,               &
     811             :                       time1   = tidx1,              &
     812             :                       time2   = tidx2,              &
     813             :                       ncArr   = ncArr,              &
     814             :                       varUnit = thisUnit,           &
     815             :                       wgt1    = wgt1,               &
     816             :                       wgt2    = wgt2,               &
     817             :                       MissVal = HCO_MISSVAL,        &
     818             :                       ArbIdx  = ArbIdx,             &
     819           0 :                       RC      = NCRC                 )
     820             : 
     821           0 :     IF ( NCRC /= 0 ) THEN
     822           0 :        CALL HCO_ERROR( 'NC_READ_ARRAY', RC )
     823           0 :        RETURN
     824             :     ENDIF
     825             : 
     826             :     ! Check for missing values: set base emissions and masks to 0, and
     827             :     ! scale factors to 1. This will make sure that these entries will
     828             :     ! be ignored.
     829             : !!! CALL CheckMissVal ( Lct, ncArr )
     830             : 
     831             :     !-----------------------------------------------------------------
     832             :     ! Eventually do interpolation between files. This is a pretty
     833             :     ! crude implementation for data interpolation between different
     834             :     ! files. It is only applied to data that is marked as interpolated
     835             :     ! data and if no appropriate interpolation date could be found in
     836             :     ! the first file. This will only be the case if the preferred date-
     837             :     ! time is outside the file range.
     838             :     !-----------------------------------------------------------------
     839           0 :     IF ( Lct%Dct%Dta%CycleFlag == HCO_CFLAG_INTER .AND. wgt1 < 0.0_sp ) THEN
     840             : 
     841             :        ! No need to read another file if the previous file had exactly the
     842             :        ! time stamp we were looking for.
     843           0 :        IF ( oYMDhm1 == YMDhma ) THEN
     844           0 :           FOUND = .FALSE.
     845             :        ELSE
     846             :           ! Check if there exists another file for a future/previous date
     847             :           ! oYMDhm1 is the originally preferred date, YMDhma is the date read
     848             :           ! from the file. If YMDhma is in the future, we need to look for a
     849             :           ! file in the past. Else, need to look for a file in the future.
     850           0 :           IF ( oYMDhm1 < YMDhma ) THEN
     851           0 :              Direction = -1
     852             :           ELSE
     853           0 :              Direction = +1
     854             :           ENDIF
     855             :           CALL SrcFile_Parse ( HcoState,  Lct, srcFile2, &
     856           0 :                                FOUND, RC, Direction = Direction )
     857           0 :           IF ( RC /= HCO_SUCCESS ) THEN
     858           0 :               CALL HCO_ERROR( 'ERROR 5', RC, THISLOC=LOC )
     859           0 :               RETURN
     860             :           ENDIF
     861             :        ENDIF
     862             : 
     863             :        ! If found, read data. Assume that all meta-data is the same.
     864           0 :        IF ( FOUND ) THEN
     865             : 
     866             :           ! Open file
     867           0 :           CALL NC_OPEN ( TRIM(srcFile2), ncLun2 )
     868             : 
     869             :           ! Define time stamp to be read. Use this call only
     870             :           ! to get the datetime of the first time slice (YMDhm1).
     871             :           ! All other values will be ignored and reset below.
     872             :           CALL GET_TIMEIDX ( HcoState,  Lct,               &
     873             :                              ncLun2,    tidx1,    tidx2,   &
     874             :                              wgt1,      wgt2,     oYMDhm2, &
     875           0 :                              YMDhmb,    YMDhm1,   RC       )
     876           0 :           IF ( RC /= HCO_SUCCESS ) THEN
     877           0 :               CALL HCO_ERROR( 'ERROR 6', RC, THISLOC=LOC )
     878           0 :               RETURN
     879             :           ENDIF
     880             : 
     881             :           ! Always read first time slice
     882           0 :           tidx1 = 1
     883           0 :           tidx2 = 1
     884           0 :           wgt1  = -1.0_sp
     885           0 :           wgt2  = -1.0_sp
     886             : 
     887             :           ! Read data and write into array ncArr2
     888             :           CALL NC_READ_ARR( fID     = ncLun2,             &
     889             :                             ncVar   = Lct%Dct%Dta%ncPara, &
     890             :                             lon1    = 1,                  &
     891             :                             lon2    = nlon,               &
     892             :                             lat1    = 1,                  &
     893             :                             lat2    = nlat,               &
     894             :                             lev1    = lev1,               &
     895             :                             lev2    = lev2,               &
     896             :                             time1   = tidx1,              &
     897             :                             time2   = tidx2,              &
     898             :                             ncArr   = ncArr2,             &
     899             :                             varUnit = thisUnit,           &
     900             :                             wgt1    = wgt1,               &
     901             :                             wgt2    = wgt2,               &
     902             :                             MissVal = HCO_MISSVAL,        &
     903             :                             ArbIdx  = ArbIdx,             &
     904           0 :                             RC      = NCRC                 )
     905           0 :           IF ( NCRC /= 0 ) THEN
     906           0 :              CALL HCO_ERROR( 'NC_READ_ARRAY (2)', RC )
     907           0 :              RETURN
     908             :           ENDIF
     909             : 
     910             :           ! Eventually fissing values
     911             : !!!          CALL CheckMissVal ( Lct, ncArr2 )
     912             : 
     913             :           ! Calculate weights to be applied to ncArr2 and ncArr1. These
     914             :           ! weights are calculated based on the originally preferred
     915             :           ! datetime oYMDh1 and the selected datetime of file 1 (YMDhma)
     916             :           ! and file 2 (YMDhm1)
     917             :           ! If date on file 1 < date on file 2:
     918           0 :           IF ( YMDhma < YMDhm1 ) THEN
     919           0 :              CALL GetWeights ( YMDhma, YMDhm1, oYMDhm1, wgt1, wgt2 )
     920             :           ! If date on file 1 > date on file 2:
     921           0 :           ELSEIF ( YMDhma > YMDhm1 ) THEN
     922           0 :              CALL GetWeights ( YMDhm1, YMDhma, oYMDhm1, wgt2, wgt1 )
     923             :           ! If both datetimes are for some reason the same (this should
     924             :           ! not happen!)
     925             :           ELSE
     926           0 :              wgt1 = 0.5_sp
     927           0 :              wgt2 = 0.5_sp
     928             :           ENDIF
     929             : 
     930             :           ! Apply weights
     931           0 :           ncArr = (wgt1 * ncArr) + (wgt2 * ncArr2)
     932             : 
     933             :           ! Verbose
     934           0 :           IF ( HCO_IsVerb( HcoState%Config%Err ) ) THEN
     935           0 :              MSG = 'Interpolated data between two files:'
     936           0 :              CALL HCO_MSG(HcoState%Config%Err,MSG)
     937           0 :              MSG = '- File 1: ' // TRIM(srcFile)
     938           0 :              CALL HCO_MSG(HcoState%Config%Err,MSG)
     939           0 :              WRITE(MSG,*) '   Time stamp used: ', YMDhma
     940           0 :              CALL HCO_MSG(HcoState%Config%Err,MSG)
     941           0 :              WRITE(MSG,*) '   Applied weight: ', wgt1
     942           0 :              CALL HCO_MSG(HcoState%Config%Err,MSG)
     943           0 :              MSG = '- File 2: ' // TRIM(srcFile2)
     944           0 :              CALL HCO_MSG(HcoState%Config%Err,MSG)
     945           0 :              WRITE(MSG,*) '   Time stamp used: ', YMDhm1
     946           0 :              CALL HCO_MSG(HcoState%Config%Err,MSG)
     947           0 :              WRITE(MSG,*) '   Applied weight: ', wgt2
     948           0 :              CALL HCO_MSG(HcoState%Config%Err,MSG)
     949             :           ENDIF
     950             : 
     951             :           ! Cleanup
     952           0 :           IF ( ASSOCIATED(ncArr2) ) DEALLOCATE(ncArr2)
     953             : 
     954             :           ! Close file
     955           0 :           CALL NC_CLOSE ( ncLun2 )
     956             :        ENDIF !FOUND
     957             : 
     958             :     !-----------------------------------------------------------------
     959             :     ! Eventually calculate averages. Currently, averages are only
     960             :     ! calculated on the year dimension, e.g. over years.
     961             :     !-----------------------------------------------------------------
     962           0 :     ELSEIF ( Lct%Dct%Dta%CycleFlag == HCO_CFLAG_AVERG    .OR. &
     963             :              Lct%Dct%Dta%CycleFlag == HCO_CFLAG_RANGEAVG       ) THEN
     964             : 
     965             :        ! cYr is the current simulation year
     966             :        CALL HcoClock_Get( HcoState%Clock, cYYYY=cYr, cMM=cMt, cDD=cDy, &
     967           0 :                           cH=cHr, RC=RC )
     968           0 :        IF ( RC /= HCO_SUCCESS ) THEN
     969           0 :            CALL HCO_ERROR( 'ERROR 7', RC, THISLOC=LOC )
     970           0 :            RETURN
     971             :        ENDIF
     972             : 
     973             :        ! Determine year range to be read:
     974             :        ! By default, we would like to average between the year range given
     975             :        ! in the time attribute
     976           0 :        Yr1 = Lct%Dct%Dta%ncYrs(1)
     977           0 :        Yr2 = Lct%Dct%Dta%ncYrs(2)
     978             : 
     979             :        ! If averaging shall only be performed if outside the given
     980             :        ! range, check if current simulation date is within the range
     981             :        ! provied in the configuration file. If so, set year range to
     982             :        ! be read to current year only.
     983           0 :        IF ( ( Lct%Dct%Dta%CycleFlag == HCO_CFLAG_RANGEAVG ) ) THEN
     984           0 :           IF ( tIDx_IsInRange(Lct,cYr,cMt,cDy,cHr) ) THEN
     985           0 :              Yr1 = cYr
     986           0 :              Yr2 = cYr
     987             :           ENDIF
     988             :        ENDIF
     989             : 
     990             :        ! Total number of years to be read
     991           0 :        nYears = Yr2 - Yr1 + 1
     992             : 
     993             :        ! Read and add annual data if there is more than one year to be
     994             :        ! used.
     995           0 :        IF ( nYears > 1 ) THEN
     996             : 
     997             :           ! Cleanup ncArr. This is refilled again
     998           0 :           ncArr = 0.0_sp
     999             : 
    1000           0 :           DO iYear = Yr1, Yr2
    1001             : 
    1002             :              ! Get file name for this year
    1003             :              CALL SrcFile_Parse ( HcoState, Lct, srcFile2, &
    1004           0 :                                   FOUND, RC, Year=iYear )
    1005           0 :              IF ( RC /= HCO_SUCCESS ) THEN
    1006           0 :                  CALL HCO_ERROR( 'ERROR 8', RC, THISLOC=LOC )
    1007           0 :                  RETURN
    1008             :              ENDIF
    1009             : 
    1010             :              ! If found, read data. Assume that all meta-data is the same.
    1011           0 :              IF ( .NOT. FOUND ) THEN
    1012           0 :                 WRITE(MSG,*) 'Cannot find file for year ', iYear, ' - needed ', &
    1013           0 :                    'to perform time-averaging on file ', TRIM(Lct%Dct%Dta%ncFile)
    1014           0 :                 CALL HCO_ERROR( MSG, RC )
    1015           0 :                 RETURN
    1016             :              ENDIF
    1017             : 
    1018             :              ! Open file
    1019           0 :              CALL NC_OPEN ( TRIM(srcFile2), ncLun2 )
    1020             : 
    1021             :              ! Define time stamp to be read.
    1022             :              CALL GET_TIMEIDX ( HcoState,  Lct,               &
    1023             :                                 ncLun2,    tidx1,    tidx2,   &
    1024             :                                 wgt1,      wgt2,     oYMDhm2, &
    1025             :                                 YMDhmb,    YMDhm1,   RC,      &
    1026           0 :                                 Year=iYear                    )
    1027           0 :              IF ( RC /= HCO_SUCCESS ) THEN
    1028           0 :                  CALL HCO_ERROR( 'ERROR 9', RC, THISLOC=LOC )
    1029           0 :                  RETURN
    1030             :              ENDIF
    1031             : 
    1032             :              ! Do not perform weights
    1033           0 :              wgt1  = -1.0_sp
    1034           0 :              wgt2  = -1.0_sp
    1035             : 
    1036             :              ! Read data and write into array ncArr2
    1037             :              CALL NC_READ_ARR( fID     = ncLun2,             &
    1038             :                                ncVar   = Lct%Dct%Dta%ncPara, &
    1039             :                                lon1    = 1,                  &
    1040             :                                lon2    = nlon,               &
    1041             :                                lat1    = 1,                  &
    1042             :                                lat2    = nlat,               &
    1043             :                                lev1    = lev1,               &
    1044             :                                lev2    = lev2,               &
    1045             :                                time1   = tidx1,              &
    1046             :                                time2   = tidx2,              &
    1047             :                                ncArr   = ncArr2,             &
    1048             :                                varUnit = thisUnit,           &
    1049             :                                wgt1    = wgt1,               &
    1050             :                                wgt2    = wgt2,               &
    1051             :                                MissVal = HCO_MISSVAL,        &
    1052             :                                ArbIdx  = ArbIdx,             &
    1053           0 :                                RC      = NCRC                 )
    1054           0 :              IF ( NCRC /= 0 ) THEN
    1055           0 :                 CALL HCO_ERROR( 'NC_READ_ARRAY (3)', RC )
    1056           0 :                 RETURN
    1057             :              ENDIF
    1058             : 
    1059             :              ! Eventually fissing values
    1060             : !!!             CALL CheckMissVal ( Lct, ncArr2 )
    1061             : 
    1062             :              ! Add all values to ncArr
    1063           0 :              ncArr = ncArr + ncArr2
    1064             : 
    1065             :              ! Cleanup
    1066           0 :              IF ( ASSOCIATED(ncArr2) ) DEALLOCATE(ncArr2)
    1067             : 
    1068             :              ! Close file
    1069           0 :              CALL NC_CLOSE ( ncLun2 )
    1070             : 
    1071             :           ENDDO !iYear
    1072             : 
    1073             :           ! Now calculate average
    1074           0 :           ncArr = ncArr / REAL(nYears,sp)
    1075             : 
    1076             :           ! Verbose
    1077           0 :           IF ( HcoState%amIRoot .AND. HCO_IsVerb( HcoState%Config%Err ) ) THEN
    1078           0 :              WRITE(MSG,110) TRIM(Lct%Dct%cName), Yr1, Yr2
    1079           0 :              CALL HCO_MSG(HcoState%Config%Err,MSG)
    1080             :           ENDIF
    1081             :  110      FORMAT( 'Field ', a, ': Average data over years ', I4.4, ' to ', I4.4 )
    1082             : 
    1083             :        ENDIF !nYears>1
    1084             :     ENDIF !Averaging
    1085             : 
    1086             :     !-----------------------------------------------------------------
    1087             :     ! Convert to HEMCO units
    1088             :     ! HEMCO data are all in kg/m2/s for fluxes and kg/m3 for
    1089             :     ! concentrations. Unit conversion is performed based on the
    1090             :     ! unit on the input file and the srcUnit attribute given in the
    1091             :     ! configuration file. By default, HEMCO will attempt to convert
    1092             :     ! the units found in the input file to the standard quantities
    1093             :     ! for mass (kg), area (m2 or m3), and time (s). For instance,
    1094             :     ! g/cm2/hr will be converted to kg/m2/s. The exceptions to this
    1095             :     ! rule are:
    1096             :     ! 1. If srcUnit is set to '1', the input data are expected to
    1097             :     !    be unitless. If the units string on the input file is none
    1098             :     !    of the units recognized by HEMCO as unitless, an error is
    1099             :     !    returned if the unit tolerance setting is set to zero, or
    1100             :     !    a warning is prompted if unit tolerance is greater than zero.
    1101             :     ! 2. If srcUnit is set to 'count', no unit conversion is performed
    1102             :     !    and data will be treated as 'index' data, e.g. regridding will
    1103             :     !    preserve the absolute values.
    1104             :     !-----------------------------------------------------------------
    1105             : 
    1106             :     ! If OrigUnit is set to wildcard character: use unit from source file
    1107           0 :     IF ( TRIM(Lct%Dct%Dta%OrigUnit) == &
    1108             :          TRIM(HCO_GetOpt(HcoState%Config%ExtList,'Wildcard')) ) THEN
    1109           0 :        Lct%Dct%Dta%OrigUnit = TRIM(thisUnit)
    1110             :     ENDIF
    1111             : 
    1112             :     ! If OrigUnit is set to '1' or to 'count', perform no unit
    1113             :     ! conversion.
    1114           0 :     IF ( HCO_IsUnitLess(Lct%Dct%Dta%OrigUnit)  .OR. &
    1115             :          HCO_IsIndexData(Lct%Dct%Dta%OrigUnit)       ) THEN
    1116             : 
    1117             :        ! Check if file unit is also unitless. This will return 0 for
    1118             :        ! unitless, 1 for HEMCO emission unit, 2 for HEMCO conc. unit,
    1119             :        ! -1 otherwise.
    1120           0 :        Flag = HCO_UNIT_SCALCHECK( thisUnit )
    1121             : 
    1122             :        ! Return with error if: (1) thisUnit is recognized as HEMCO unit and
    1123             :        ! unit tolerance is set to zero; (2) thisUnit is neither unitless nor
    1124             :        ! a HEMCO unit and unit tolerance is set to zero or one.
    1125             :        ! The unit tolerance is defined in the configuration file.
    1126           0 :        IF ( Flag /= 0 .AND. UnitTolerance == 0 ) THEN
    1127             :           MSG = 'Illegal unit: ' // TRIM(thisUnit) // '. File: ' // &
    1128           0 :                 TRIM(srcFile)
    1129           0 :           CALL HCO_ERROR( MSG, RC )
    1130           0 :           RETURN
    1131             :        ENDIF
    1132             : 
    1133             :        ! Prompt a warning if thisUnit is not recognized as unitless.
    1134           0 :        IF ( Flag /= 0 ) THEN
    1135             :           MSG = 'Data is treated as unitless, but file attribute suggests ' // &
    1136           0 :                 'it is not: ' // TRIM(thisUnit) // '. File: ' // TRIM(srcFile)
    1137           0 :           CALL HCO_WARNING( HcoState%Config%Err, MSG, RC )
    1138             :        ENDIF
    1139             : 
    1140             :        ! Verbose mode
    1141           0 :        IF ( HCO_IsVerb( HcoState%Config%Err ) ) THEN
    1142           0 :           WRITE(MSG,*) 'Based on srcUnit attribute (', TRIM(Lct%Dct%Dta%OrigUnit), &
    1143           0 :                        '), no unit conversion is performed.'
    1144           0 :           CALL HCO_MSG(HcoState%Config%Err,MSG)
    1145             :        ENDIF
    1146             : 
    1147             :     ! Convert to HEMCO units in all other cases.
    1148             :     ELSE
    1149             : 
    1150             :        ! For zero unit tolerance, make sure that thisUnit matches
    1151             :        ! with unit set in configuration file. For higher unit
    1152             :        ! tolerances, prompt a level 3 warning.
    1153           0 :        IF ( TRIM(Lct%Dct%Dta%OrigUnit) /= TRIM(thisUnit) ) THEN
    1154             :           MSG = 'File units do not match: ' // TRIM(thisUnit) // &
    1155             :                 ' vs. ' // TRIM(Lct%Dct%Dta%OrigUnit)    // &
    1156           0 :                 '. File: ' // TRIM(srcFile)
    1157             : 
    1158           0 :           IF ( UnitTolerance == 0 ) THEN
    1159           0 :              CALL HCO_ERROR( MSG, RC )
    1160           0 :              RETURN
    1161             :           ELSE
    1162           0 :              CALL HCO_WARNING( HcoState%Config%Err, MSG, RC )
    1163             :           ENDIF
    1164             :        ENDIF
    1165             : 
    1166             :        ! Shadow species properties needed for unit conversion
    1167           0 :        HcoID = Lct%Dct%HcoID
    1168           0 :        IF ( HcoID > 0 ) THEN
    1169           0 :           MW_g = HcoState%Spc(HcoID)%MW_g
    1170             :        ELSE
    1171           0 :           MW_g = -999.0_hp
    1172             :        ENDIF
    1173             : 
    1174             :        ! Now convert to HEMCO units. This attempts to convert mass,
    1175             :        ! area/volume and time to HEMCO standards (kg, m2/m3, s).
    1176           0 :        ncYr  = FLOOR( MOD( oYMDhm1, 1.0e12_dp ) / 1.0e8_dp )
    1177           0 :        ncMt  = FLOOR( MOD( oYMDhm1, 1.0e8_dp  ) / 1.0e6_dp )
    1178             : 
    1179           0 :        IF ( ncYr == 0 ) THEN
    1180           0 :           CALL HcoClock_Get( HcoState%Clock, cYYYY = ncYr, RC=RC )
    1181           0 :           IF ( RC /= HCO_SUCCESS ) THEN
    1182           0 :               CALL HCO_ERROR( 'ERROR 10', RC, THISLOC=LOC )
    1183           0 :               RETURN
    1184             :           ENDIF
    1185             :        ENDIF
    1186           0 :        IF ( ncMt == 0 ) THEN
    1187           0 :           CALL HcoClock_Get( HcoState%Clock, cMM   = ncMt, RC=RC )
    1188           0 :           IF ( RC /= HCO_SUCCESS ) THEN
    1189           0 :               CALL HCO_ERROR( 'ERROR 11', RC, THISLOC=LOC )
    1190           0 :               RETURN
    1191             :           ENDIF
    1192             :        ENDIF
    1193             : 
    1194             :        ! Verbose mode
    1195           0 :        IF ( HcoState%amIRoot .and. HCO_IsVerb( HcoState%Config%Err ) ) THEN
    1196           0 :           WRITE(MSG,*) 'Unit conversion settings: '
    1197           0 :           CALL HCO_MSG(HcoState%Config%Err,MSG)
    1198           0 :           WRITE(MSG,*) '- Year, month        : ', ncYr, ncMt
    1199           0 :           CALL HCO_MSG(HcoState%Config%Err,MSG)
    1200             :        ENDIF
    1201             : 
    1202             :        CALL HCO_UNIT_CHANGE(                 &
    1203             :             HcoConfig     = HcoState%Config, &
    1204             :             Array         = ncArr,           &
    1205             :             Units         = thisUnit,        &
    1206             :             MW            = MW_g,            &
    1207             :             YYYY          = ncYr,            &
    1208             :             MM            = ncMt,            &
    1209             :             AreaFlag      = AreaFlag,        &
    1210             :             TimeFlag      = TimeFlag,        &
    1211             :             FACT          = UnitFactor,      &
    1212           0 :             RC            = RC                )
    1213           0 :        IF ( RC /= HCO_SUCCESS ) THEN
    1214           0 :           MSG = 'Cannot convert units for ' // TRIM(Lct%Dct%cName)
    1215           0 :           CALL HCO_ERROR( MSG , RC )
    1216           0 :           RETURN
    1217             :        ENDIF
    1218             : 
    1219             :        ! Verbose mode
    1220           0 :        IF ( UnitFactor /= 1.0_hp ) THEN
    1221           0 :           IF ( HCO_IsVerb( HcoState%Config%Err ) ) THEN
    1222           0 :              WRITE(MSG,*) 'Data was in units of ', TRIM(thisUnit), &
    1223           0 :                           ' - converted to HEMCO units by applying ', &
    1224           0 :                           'scale factor ', UnitFactor
    1225           0 :              CALL HCO_MSG(HcoState%Config%Err,MSG)
    1226             :           ENDIF
    1227             :        ELSE
    1228           0 :           IF ( HCO_IsVerb( HcoState%Config%Err ) ) THEN
    1229           0 :              WRITE(MSG,*) 'Data was in units of ', TRIM(thisUnit), &
    1230           0 :                           ' - unit conversion factor is ', UnitFactor
    1231           0 :              CALL HCO_MSG(HcoState%Config%Err,MSG)
    1232             :           ENDIF
    1233             :        ENDIF
    1234             : 
    1235             :        ! Check for valid unit combinations, i.e. emissions must be kg/m2/s,
    1236             :        ! concentrations kg/m3. Eventually multiply by emission time step
    1237             :        ! or divide by area to obtain those values.
    1238             : 
    1239             :        ! Concentration data
    1240           0 :        IF ( AreaFlag == 3 .AND. TimeFlag == 0 ) THEN
    1241           0 :           Lct%Dct%Dta%IsConc = .TRUE.
    1242             : 
    1243             :        ! If concentration data is per second (kg/m3/s), multiply by emission
    1244             :        ! time step to get concentration (kg/m3).
    1245           0 :        ELSEIF ( AreaFlag == 3 .AND. TimeFlag == 1 ) THEN
    1246           0 :           Lct%Dct%Dta%IsConc = .TRUE.
    1247             : 
    1248           0 :           ncArr = ncArr * HcoState%TS_EMIS
    1249             :           MSG = 'Data converted from kg/m3/s to kg/m3: ' // &
    1250           0 :                 TRIM(Lct%Dct%cName) // ': ' // TRIM(thisUnit)
    1251           0 :           CALL HCO_WARNING( HcoState%Config%Err, MSG, RC )
    1252             : 
    1253             :        ! Unitless data
    1254           0 :        ELSEIF ( AreaFlag == -1 .AND. TimeFlag == -1 ) THEN
    1255             :           ! nothing do to
    1256             : 
    1257             :        ! Emission data
    1258           0 :        ELSEIF ( AreaFlag == 2 .AND. TimeFlag == 1 ) THEN
    1259             :           ! nothing do to
    1260             : 
    1261             :        ! Emission data that is not per time (kg/m2): convert to kg/m2/s
    1262           0 :        ELSEIF ( AreaFlag == 2 .AND. TimeFlag == 0 ) THEN
    1263           0 :           ncArr = ncArr / HcoState%TS_EMIS
    1264             :           MSG = 'Data converted from kg/m2 to kg/m2/s: ' // &
    1265           0 :                 TRIM(Lct%Dct%cName) // ': ' // TRIM(thisUnit)
    1266           0 :           CALL HCO_WARNING( HcoState%Config%Err, MSG, RC )
    1267             : 
    1268             :        ! Emission data that is not per area (i.e. kg/s) needs to be converted
    1269             :        ! to per area manually.
    1270           0 :        ELSEIF ( AreaFlag == 0 .AND. TimeFlag == 1 ) THEN
    1271             : 
    1272             :           ! Get lat edges: those are read from file if possible, otherwise
    1273             :           ! calculated from the lat midpoints.
    1274             :           ! ==> Sine of lat is needed. Do conversion right here.
    1275             :           CALL NC_GET_GRID_EDGES ( ncLun, 2, LatMid,   nlat, &
    1276           0 :                                    LatEdge,  nlatEdge, NCRC   )
    1277           0 :           IF ( NCRC /= 0 ) THEN
    1278           0 :              MSG = 'Cannot read lat edge of ' // TRIM(srcFile)
    1279           0 :              CALL HCO_ERROR( MSG, RC )
    1280           0 :              RETURN
    1281             :           ENDIF
    1282             : 
    1283             :           ! Now normalize data by area calculated from lat edges.
    1284             :           CALL NORMALIZE_AREA( HcoState, ncArr,   nlon, &
    1285           0 :                                LatEdge,  srcFile, RC     )
    1286           0 :           IF ( RC /= HCO_SUCCESS ) THEN
    1287           0 :               CALL HCO_ERROR( 'ERROR 12', RC, THISLOC=LOC )
    1288           0 :               RETURN
    1289             :           ENDIF
    1290             : 
    1291             :        ! All other combinations are invalid
    1292             :        ELSE
    1293             :           MSG = 'Unit must be unitless, emission or concentration: ' // &
    1294           0 :                 TRIM(Lct%Dct%cName) // ': ' // TRIM(thisUnit)
    1295           0 :           CALL HCO_ERROR( MSG, RC )
    1296           0 :           RETURN
    1297             :        ENDIF
    1298             :     ENDIF ! Unit conversion
    1299             : 
    1300             :     !-----------------------------------------------------------------
    1301             :     ! Get horizontal grid edges
    1302             :     !-----------------------------------------------------------------
    1303             : 
    1304             :     ! Get longitude edges and make sure they are steadily increasing.
    1305             :     CALL NC_GET_GRID_EDGES ( ncLun, 1, LonMid,   nlon, &
    1306           0 :                              LonEdge,  nlonEdge, NCRC   )
    1307           0 :     IF ( NCRC /= 0 ) THEN
    1308           0 :        MSG = 'Cannot read lon edge of ' // TRIM(srcFile)
    1309           0 :        CALL HCO_ERROR( MSG, RC )
    1310           0 :        RETURN
    1311             :     ENDIF
    1312           0 :     CALL HCO_ValidateLon( HcoState, nlonEdge, LonEdge, RC )
    1313           0 :     IF ( RC /= HCO_SUCCESS ) THEN
    1314           0 :         CALL HCO_ERROR( 'ERROR 13', RC, THISLOC=LOC )
    1315           0 :         RETURN
    1316             :     ENDIF
    1317             : 
    1318             :     ! Get latitude edges (only if they have not been read yet
    1319             :     ! for unit conversion)
    1320           0 :     IF ( .NOT. ASSOCIATED( LatEdge ) ) THEN
    1321             :        CALL NC_GET_GRID_EDGES ( ncLun, 2, LatMid,   nlat, &
    1322           0 :                                 LatEdge,  nlatEdge, NCRC   )
    1323           0 :        IF ( NCRC /= 0 ) THEN
    1324           0 :           MSG = 'Cannot read lat edge of ' // TRIM(srcFile)
    1325           0 :           CALL HCO_ERROR( MSG, RC )
    1326           0 :           RETURN
    1327             :        ENDIF
    1328             :     ENDIF
    1329             : 
    1330             :     !-----------------------------------------------------------------
    1331             :     ! Determine regridding algorithm to be applied: use NCREGRID from
    1332             :     ! MESSy only if we need to regrid vertical levels. For all other
    1333             :     ! fields, use the much faster map_a2a.
    1334             :     ! Perform no vertical regridding if the vertical levels are model
    1335             :     ! levels. Model levels are assumed to start at the surface, i.e.
    1336             :     ! the first input level must correspond to the surface level. The
    1337             :     ! total number of vertical levels must not match the number of
    1338             :     ! vertical levels on the simulation grid. Data is not extrapolated
    1339             :     ! beyond the existing levels.
    1340             :     ! Vertical regridding based on NCREGRID will always map the input
    1341             :     ! data onto the entire simulation grid (no extrapolation beyond
    1342             :     ! the vertical input coordinates).
    1343             :     ! Index-based remapping can currently not be done with the MESSy
    1344             :     ! routines, i.e. it is not possible to vertically regrid index-
    1345             :     ! based data.
    1346             :     !-----------------------------------------------------------------
    1347             : 
    1348           0 :     UseMESSy = .FALSE.
    1349           0 :     IF ( nlev > 1 .AND. .NOT. IsModelLevel ) THEN
    1350           0 :        UseMESSy = .TRUE.
    1351             :     ENDIF
    1352             : 
    1353             : #if defined( MODEL_CESM ) || defined( MODEL_WRF )
    1354             :     ! If in WRF or the CESM environment, the vertical grid is arbitrary.
    1355             :     ! MESSy regridding ALWAYS has to be used.
    1356           0 :     IF ( nlev > 1 ) THEN
    1357           0 :       UseMESSy = .TRUE.
    1358             : 
    1359           0 :       IF ( HCO_IsVerb( HcoState%Config%Err ) ) THEN
    1360           0 :         WRITE(MSG,*) '  ==> WRF/CESM: Always forcing MESSy regridding for number of verticals', nlev, IsModelLevel
    1361           0 :         CALL HCO_MSG(HcoState%Config%Err,MSG)
    1362             :       ENDIF
    1363             :     ENDIF
    1364             : #endif
    1365             : 
    1366           0 :     IF ( HCO_IsIndexData(Lct%Dct%Dta%OrigUnit) .AND. UseMESSy ) THEN
    1367             :        MSG = 'Cannot do MESSy regridding for index data: ' // &
    1368           0 :              TRIM(srcFile)
    1369           0 :        CALL HCO_ERROR( MSG, RC )
    1370           0 :        RETURN
    1371             :     ENDIF
    1372             : 
    1373             :     !-----------------------------------------------------------------
    1374             :     ! Use MESSy regridding
    1375             :     !-----------------------------------------------------------------
    1376           0 :     IF ( UseMESSy ) THEN
    1377           0 :        IF ( HCO_IsVerb( HcoState%Config%Err ) ) THEN
    1378           0 :           WRITE(MSG,*) '  ==> Use MESSy regridding (NCREGRID)'
    1379           0 :           CALL HCO_MSG(HcoState%Config%Err,MSG)
    1380             :        ENDIF
    1381             : 
    1382             : #if !defined( MODEL_CESM ) && !defined( MODEL_WRF )
    1383             :        ! If we do MESSy regridding, we can only do one time step
    1384             :        ! at a time at the moment!
    1385             :        IF ( tidx1 /= tidx2 ) THEN
    1386             :           MSG = 'Cannot do MESSy regridding for more than one time step; ' &
    1387             :                 // TRIM(srcFile)
    1388             :           CALL HCO_ERROR( MSG, RC )
    1389             :           RETURN
    1390             :        ENDIF
    1391             : 
    1392             :        ! Note: This seems to be a soft restriction - removing this
    1393             :        ! does not conflict with MESSy regridding. Need to check (hplin, 5/30/20)
    1394             :        ! This has to be used for WRF-GC and CESM so ifdefd out
    1395             : #endif
    1396             : 
    1397             : #if defined( MODEL_WRF ) || defined( MODEL_CESM )
    1398             :        !--------------------------------------------------------------
    1399             :        ! Eventually get sigma levels
    1400             :        ! For files that have hardcoded GEOS-Chem "index"-based levels,
    1401             :        ! translate these levels back into a sigma representation
    1402             :        ! of the GEOS-Chem levels (sigma = p/ps on INTERFACE)
    1403             :        !
    1404             :        ! There are caveats with this. This is essentially a copy of the
    1405             :        ! hardcoded hPa lists from
    1406             :        ! http://wiki.seas.harvard.edu/geos-chem/index.php/GEOS-Chem_vertical_grids
    1407             :        ! hard-coded by hand, and we only assume that the data is either
    1408             :        ! 47-levels or 72-levels.
    1409             :        !
    1410             :        ! Parse the 72 list using regex like so: ^ ?\d{1,2} then remove the lines
    1411             :        ! Then you have the 73 edges.
    1412             :        !
    1413             :        ! psfc = PEDGE(0) = 1013.250 hPa
    1414             :        !
    1415             :        ! Ported from the original WRF-GC implementation (hplin, 5/27/20)
    1416             :        !--------------------------------------------------------------
    1417           0 :        IF ( nlev > 1 .AND. IsModelLevel ) THEN
    1418           0 :           IF ( HCO_IsVerb( HcoState%Config%Err ) ) THEN
    1419           0 :             WRITE(MSG,*) '  ==> WRF/CESM: Writing in fixed sigma coordinates for GEOS-Chem levels', nlon, nlat
    1420           0 :             CALL HCO_MSG(HcoState%Config%Err,MSG)
    1421             :           ENDIF
    1422             : 
    1423           0 :           ALLOCATE(SigEdge(nlon, nlat, nlev))
    1424             : 
    1425           0 :           DO I = 1, nlon
    1426           0 :              DO J = 1, nlat
    1427             :                ! Fill with pre-defined, hard coded sigma levels computed.
    1428           0 :                SigEdge(I, J, :) = GC_72_EDGE_SIGMA(1:nlev)
    1429             :              ENDDO
    1430             :            ENDDO
    1431             :        ENDIF
    1432             : #endif
    1433             : 
    1434             :        !--------------------------------------------------------------
    1435             :        ! Eventually get sigma levels
    1436             :        ! Vertical regridding is performed on sigma interface levels:
    1437             :        ! sigma(i,j,l) = p(i,j,l) / ps(i,j)
    1438             :        ! NC_Get_Sigma_Levels attempts to create the sigma levels from
    1439             :        ! the content of the netCDF file.
    1440             :        ! For now, it is assumed that all input data is on vertical
    1441             :        ! mid-point levels, and the interface values are calculated
    1442             :        ! by linear interpolation of the mid-point values in a second
    1443             :        ! step.
    1444             :        ! For model levels, the sigma levels don't need to be known
    1445             :        ! as vertical interpolation will be done based on subroutine
    1446             :        ! ModelLev_Interpolate (within HCO_MESSY_REGRID).
    1447             :        !--------------------------------------------------------------
    1448           0 :        IF ( nlev > 1 .AND. .NOT. IsModelLevel ) THEN
    1449             : 
    1450             :           ! Get sigma levels
    1451             :           CALL NC_Get_Sigma_Levels ( fID     = ncLun,   &
    1452             :                                      ncFile  = srcFile, &
    1453             :                                      levName = LevName, &
    1454             :                                      lon1    = 1,       &
    1455             :                                      lon2    = nlon,    &
    1456             :                                      lat1    = 1,       &
    1457             :                                      lat2    = nlat,    &
    1458             :                                      lev1    = 1,       &
    1459             :                                      lev2    = nlev,    &
    1460             :                                      time    = tidx1,   &
    1461             :                                      SigLev  = SigLev,  &
    1462             :                                      Dir     = dir,     &
    1463           0 :                                      RC      = NCRC      )
    1464           0 :           IF ( NCRC /= 0 ) THEN
    1465           0 :              CALL HCO_ERROR( 'Cannot read sigma levels of '//TRIM(srcFile), RC )
    1466           0 :              RETURN
    1467             :           ENDIF
    1468             : 
    1469             :           ! Interpolate onto edges
    1470           0 :           CALL SigmaMidToEdges ( HcoState, SigLev, SigEdge, RC )
    1471           0 :           IF ( RC /= HCO_SUCCESS ) THEN
    1472           0 :               CALL HCO_ERROR( 'ERROR 14', RC, THISLOC=LOC )
    1473           0 :               RETURN
    1474             :           ENDIF
    1475             : 
    1476             :           ! Sigma levels are not needed anymore
    1477           0 :           IF ( ASSOCIATED(SigLev) ) DEALLOCATE(SigLev)
    1478             : 
    1479             :           !-----------------------------------------------------------
    1480             :           ! Flip vertical axis if positive axis is 'down', i.e. level
    1481             :           ! index 1 is the top of the atmosphere
    1482             :           !-----------------------------------------------------------
    1483           0 :           IF ( dir == -1 ) THEN
    1484           0 :              SigEdge(:,:,:  ) = SigEdge(:,:,nlev+1:1:-1  )
    1485           0 :              NcArr  (:,:,:,:) = NcArr  (:,:,nlev  :1:-1,:)
    1486             :           ENDIF
    1487             : 
    1488             :        ENDIF ! nlev>1
    1489             : 
    1490             : #if defined( MODEL_WRF ) || defined( MODEL_CESM )
    1491             :        ! Input data is "never" on model levels because model levels can change! (hplin, 5/29/20)
    1492           0 :        IsModelLevel = .false.
    1493             : #endif
    1494             : 
    1495             :        ! Now do the regridding
    1496             :        CALL HCO_MESSY_REGRID ( HcoState,  NcArr,                 &
    1497             :                                LonEdge,   LatEdge,      SigEdge, &
    1498           0 :                                Lct,       IsModelLevel, RC        )
    1499           0 :        IF ( RC /= HCO_SUCCESS ) THEN
    1500           0 :            CALL HCO_ERROR( 'ERROR 15', RC, THISLOC=LOC )
    1501           0 :            RETURN
    1502             :        ENDIF
    1503             : 
    1504             :        ! Cleanup
    1505           0 :        IF ( ASSOCIATED(SigEdge) ) DEALLOCATE(SigEdge)
    1506             : 
    1507             :     !-----------------------------------------------------------------
    1508             :     ! Use map_a2a regridding
    1509             :     !-----------------------------------------------------------------
    1510             :     ELSE
    1511           0 :        IF ( HCO_IsVerb( HcoState%Config%Err ) ) THEN
    1512           0 :           WRITE(MSG,*) '  ==> Use map_a2a regridding'
    1513           0 :           CALL HCO_MSG(HcoState%Config%Err,MSG)
    1514             :        ENDIF
    1515             : 
    1516           0 :        CALL REGRID_MAPA2A ( HcoState, NcArr, LonEdge, LatEdge, Lct, RC )
    1517           0 :        IF ( RC /= HCO_SUCCESS ) THEN
    1518           0 :            CALL HCO_ERROR( 'ERROR 16', RC, THISLOC=LOC )
    1519           0 :            RETURN
    1520             :        ENDIF
    1521             : 
    1522             :     ENDIF
    1523             : 
    1524             :     !-----------------------------------------------------------------
    1525             :     ! Add to diagnostics (if it exists)
    1526             :     !-----------------------------------------------------------------
    1527           0 :     IF ( HcoState%Options%Field2Diagn ) THEN
    1528           0 :        IF ( Lct%Dct%Dta%SpaceDim == 3 .AND. ASSOCIATED(Lct%Dct%Dta%V3) ) THEN
    1529           0 :           IF ( ASSOCIATED(Lct%Dct%Dta%V3(1)%Val) ) THEN
    1530             :              CALL Diagn_Update ( HcoState, cName=TRIM(Lct%Dct%cName), &
    1531           0 :                                  Array3D=Lct%Dct%Dta%V3(1)%Val, COL=-1, RC=RC )
    1532           0 :              IF ( RC /= HCO_SUCCESS ) THEN
    1533           0 :                  CALL HCO_ERROR( 'ERROR 17', RC, THISLOC=LOC )
    1534           0 :                  RETURN
    1535             :              ENDIF
    1536             :           ENDIF
    1537           0 :        ELSEIF ( Lct%Dct%Dta%SpaceDim == 2 .AND. ASSOCIATED(Lct%Dct%Dta%V2) ) THEN
    1538           0 :           IF ( ASSOCIATED(Lct%Dct%Dta%V2(1)%Val) ) THEN
    1539             :              CALL Diagn_Update ( HcoState, cName=TRIM(Lct%Dct%cName), &
    1540           0 :                                  Array2D=Lct%Dct%Dta%V2(1)%Val, COL=-1, RC=RC )
    1541           0 :              IF ( RC /= HCO_SUCCESS ) THEN
    1542           0 :                  CALL HCO_ERROR( 'ERROR 18', RC, THISLOC=LOC )
    1543           0 :                  RETURN
    1544             :              ENDIF
    1545             :           ENDIF
    1546             :        ENDIF
    1547             :     ENDIF
    1548             : 
    1549             :     !-----------------------------------------------------------------
    1550             :     ! Cleanup and leave
    1551             :     !-----------------------------------------------------------------
    1552           0 :     IF ( ASSOCIATED ( ncArr   ) ) DEALLOCATE ( ncArr   )
    1553           0 :     IF ( ASSOCIATED ( LonMid  ) ) DEALLOCATE ( LonMid  )
    1554           0 :     IF ( ASSOCIATED ( LatMid  ) ) DEALLOCATE ( LatMid  )
    1555           0 :     IF ( ASSOCIATED ( LevMid  ) ) DEALLOCATE ( LevMid  )
    1556           0 :     IF ( ASSOCIATED ( LonEdge ) ) DEALLOCATE ( LonEdge )
    1557           0 :     IF ( ASSOCIATED ( LatEdge ) ) DEALLOCATE ( LatEdge )
    1558             : 
    1559             :     ! Return w/ success
    1560           0 :     CALL HCO_LEAVE ( HcoState%Config%Err,  RC )
    1561             : 
    1562           0 :   END SUBROUTINE HCOIO_Read
    1563             : !EOC
    1564             : !------------------------------------------------------------------------------
    1565             : !                   Harmonized Emissions Component (HEMCO)                    !
    1566             : !------------------------------------------------------------------------------
    1567             : !BOP
    1568             : !
    1569             : ! !IROUTINE: HCOIO_CloseAll
    1570             : !
    1571             : ! !DESCRIPTION: Subroutine HCOIO\_CloseAll makes sure that there is no open
    1572             : ! netCDF file left in the stream.
    1573             : !\\
    1574             : !\\
    1575             : ! !INTERFACE:
    1576             : !
    1577           0 :   SUBROUTINE HCOIO_CloseAll( HcoState, RC )
    1578             : !
    1579             : ! !USES:
    1580             : !
    1581             :     USE HCO_Ncdf_Mod,   ONLY : NC_CLOSE
    1582             : !
    1583             : ! !INPUT PARAMTERS:
    1584             : !
    1585             :     TYPE(HCO_State), POINTER          :: HcoState    ! HEMCO state
    1586             : !
    1587             : ! !INPUT/OUTPUT PARAMETERS:
    1588             : !
    1589             :     INTEGER,          INTENT(INOUT)   :: RC
    1590             : !
    1591             : ! !REVISION HISTORY:
    1592             : !  24 Mar 2016 - C. Keller: Initial version
    1593             : !  See https://github.com/geoschem/hemco for complete history
    1594             : !EOP
    1595             : !------------------------------------------------------------------------------
    1596             : !BOC
    1597             : !
    1598             : ! !LOCAL VARIABLES:
    1599             : !
    1600             :     !======================================================================
    1601             :     ! HCOIO_CloseAll begins here
    1602             :     !======================================================================
    1603           0 :     IF ( HcoState%ReadLists%FileLun > 0 ) THEN
    1604           0 :        CALL NC_CLOSE( HcoState%ReadLists%FileLun )
    1605           0 :        HcoState%ReadLists%FileLun = -1
    1606             :     ENDIF
    1607             : 
    1608             :     ! Return w/ success
    1609           0 :     RC = HCO_SUCCESS
    1610             : 
    1611           0 :   END SUBROUTINE HCOIO_CloseAll
    1612             : !EOC
    1613             : END MODULE HCOIO_Read_Mod
    1614             : #endif

Generated by: LCOV version 1.14