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 22:39:59 |  | 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