LCOV - code coverage report
Current view: | top level - hemco/HEMCO/src/Core - hco_calc_mod.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 0 | 703 | 0.0 %
Date: | 2025-01-13 21:54:50 |  | Functions: | 0 | 15 | 0.0 %

              Line data    Source code
           1             : !------------------------------------------------------------------------------
           2             : !                   Harmonized Emissions Component (HEMCO)                    !
           3             : !------------------------------------------------------------------------------
           4             : !BOP
           5             : !
           6             : ! !MODULE: hco_calc_mod.F90
           7             : !
           8             : ! !DESCRIPTION: Module HCO\_Calc\_Mod contains routines to calculate
           9             : ! HEMCO core emissions based on the content of the HEMCO EmisList
          10             : ! object. All emissions are in [kg/m2/s].
          11             : !\\
          12             : !\\
          13             : ! Emissions for the current datetime are calculated by multiplying base
          14             : ! emissions fields with the associated scale factors. Different
          15             : ! inventories are merged/overlayed based upon the category and hierarchy
          16             : ! attributes assigned to the individual base fields. Within the same
          17             : ! category, fields of higher hierarchy overwrite lower-hierarchy fields.
          18             : ! Emissions of different categories are always added.
          19             : !\\
          20             : !\\
          21             : ! The assembled emission array is written into the corresponding emission
          22             : ! rates array of the HEMCO state object: HcoState%Spc(HcoID)%Emis, where
          23             : ! HcoID denotes the corresponding species ID. Emis covers dimension lon,
          24             : ! lat, lev on the HEMCO grid, i.e. unlike the emission arrays in EmisList
          25             : ! that only cover the levels defined in the source files, Emis extends
          26             : ! over all vertical model levels.
          27             : !\\
          28             : !\\
          29             : ! Negative emissions are not supported and are ignored. Surface
          30             : ! deposition velocities are stored in HcoState%Spc(HcoID)%Depv and can
          31             : ! be added therein.
          32             : !\\
          33             : !\\
          34             : ! In addition to emissions and surface deposition rates, HEMCO also
          35             : ! supports concentrations (kg/m3). Data is automatically written into
          36             : ! the concentration array HcoState%Spc(HcoID)%Conc if the source data
          37             : ! is marked as being concentration data (i.e. if Dta%IsConc is .TRUE.).
          38             : ! The latter is automatically determined by HEMCO based upon the data
          39             : ! units.
          40             : !\\
          41             : !\\
          42             : ! All emission calculation settings are passed through the HcoState
          43             : ! options object (HcoState%Options). These include:
          44             : !
          45             : ! \begin{itemize}
          46             : !  \item ExtNr: extension number to be considered.
          47             : !  \item SpcMin: lower species ID (HEMCO ID) to be considered.
          48             : !  \item SpcMax: upper species ID (HEMCO ID) to be considered. If set
          49             : !        to -1, all species above or equal to SpcMin are considered.
          50             : !  \item CatMin: lower emission category to be considered.
          51             : !  \item CatMax: upper emission category to be considered. If set to
          52             : !        -1, all categories above or equal to CatMin are considered.
          53             : !  \item FillBuffer: if set to TRUE, the emissions will be written into
          54             : !        buffer array HcoState%Buffer3D instead of HcoState%Spc(ID)%Emis.
          55             : !        If this option is enabled, only one species can be calculated at
          56             : !        a time (by setting SpcMin/SpcMax, CatMin/CatMax and/or ExtNr
          57             : !        accordingly). This option is useful for extensions, e.g. if
          58             : !        additional scalings need to be done on some emission fields
          59             : !        assembled by HEMCO (e.g. PARANOX extension).
          60             : ! \end{itemize}
          61             : !
          62             : ! !INTERFACE:
          63             : !
          64             : MODULE HCO_Calc_Mod
          65             : !
          66             : ! !USES:
          67             : !
          68             :   USE HCO_Diagn_Mod
          69             :   USE HCO_Error_Mod
          70             :   USE HCO_Types_Mod
          71             :   USE HCO_DataCont_Mod, ONLY : Pnt2DataCont
          72             :
          73             :   IMPLICIT NONE
          74             :   PRIVATE
          75             : !
          76             : ! !PUBLIC MEMBER FUNCTIONS:
          77             : !
          78             :   PUBLIC  :: HCO_CalcEmis
          79             :   PUBLIC  :: HCO_CheckDepv
          80             :   PUBLIC  :: HCO_EvalFld
          81             :   PUBLIC  :: HCO_MaskFld
          82             : #ifdef ADJOINT
          83             :   PUBLIC  :: GET_CURRENT_EMISSIONS_ADJ
          84             : #endif
          85             : !
          86             : ! !PRIVATE MEMBER FUNCTIONS:
          87             : !
          88             :   PRIVATE :: GET_CURRENT_EMISSIONS
          89             :   PRIVATE :: GetMaskVal
          90             :   PRIVATE :: GetDilFact
          91             :   PRIVATE :: GetVertIndx
          92             :   PRIVATE :: GetIdx
          93             : !
          94             : ! !PARAMETER
          95             : !
          96             :   ! Mask threshold. All mask values below this value will be evaluated
          97             :   ! as zero (= outside of mask), and all values including and above this
          98             :   ! value as inside the mask. This is only of relevance if the MaskFractions
          99             :   ! option is false. If MaskFractions is true, the fractional mask values are
         100             :   ! considered, e.g. a grid box can contribute 40% to a mask region, etc.
         101             :   ! The MaskFractions toggle can be set in the settings section of the HEMCO
         102             :   ! configuration file (Use mask fractions: true/false). It defaults to false.
         103             :   REAL(sp), PARAMETER  :: MASK_THRESHOLD = 0.5_sp
         104             : !
         105             : ! ============================================================================
         106             : !
         107             : ! !REVISION HISTORY:
         108             : !  25 Aug 2012 - C. Keller   - Initial version.
         109             : !  See https://github.com/geoschem/hemco for complete history
         110             : !EOP
         111             : !------------------------------------------------------------------------------
         112             : !BOC
         113             :   INTERFACE HCO_EvalFld
         114             :      MODULE PROCEDURE HCO_EvalFld_2D
         115             :      MODULE PROCEDURE HCO_EvalFld_3D
         116             :   END INTERFACE HCO_EvalFld
         117             :
         118             : CONTAINS
         119             : !EOC
         120             : !------------------------------------------------------------------------------
         121             : !                   Harmonized Emissions Component (HEMCO)                    !
         122             : !------------------------------------------------------------------------------
         123             : !BOP
         124             : !
         125             : ! !IROUTINE: HCO_CalcEmis
         126             : !
         127             : ! !DESCRIPTION: Subroutine HCO\_CalcEmis calculates the 3D emission
         128             : ! fields at current datetime for the specified species, categories, and
         129             : ! extension number.
         130             : !\\
         131             : !\\
         132             : ! !INTERFACE:
         133             : !
         134           0 :   SUBROUTINE HCO_CalcEmis( HcoState, UseConc, RC )
         135             : !
         136             : ! !USES:
         137             : !
         138             :     USE HCO_STATE_MOD,    ONLY : HCO_State
         139             :     USE HCO_ARR_MOD,      ONLY : HCO_ArrAssert
         140             :     USE HCO_DATACONT_MOD, ONLY : ListCont_NextCont
         141             :     USE HCO_FILEDATA_MOD, ONLY : FileData_ArrIsDefined
         142             :     USE HCO_Scale_Mod,    ONLY : HCO_ScaleArr
         143             : !
         144             : ! !INPUT PARAMETERS:
         145             : !
         146             :     LOGICAL,         INTENT(IN   )  :: UseConc    ! Use concentration fields?
         147             : !
         148             : ! !INPUT/OUTPUT PARAMETERS:
         149             : !
         150             :     TYPE(HCO_State), POINTER        :: HcoState   ! HEMCO state object
         151             :     INTEGER,         INTENT(INOUT)  :: RC         ! Return code
         152             : !
         153             : ! !REVISION HISTORY:
         154             : !  25 Aug 2012 - C. Keller   - Initial Version
         155             : !  See https://github.com/geoschem/hemco for complete history
         156             : !EOP
         157             : !------------------------------------------------------------------------------
         158             : !BOC
         159             : !
         160             : ! !LOCAL VARIABLES:
         161             : !
         162             :     ! Working pointers: list and data container
         163             :     TYPE(ListCont), POINTER :: Lct
         164             :     TYPE(DataCont), POINTER :: Dct
         165             :
         166             :     ! Temporary emission arrays
         167             :     REAL(hp), POINTER       :: OutArr(:,:,:) => NULL()
         168             :     REAL(hp), TARGET        :: SpcFlx( HcoState%NX, &
         169             :                                        HcoState%NY, &
         170           0 :                                        HcoState%NZ   )
         171             :     REAL(hp), TARGET        :: CatFlx( HcoState%NX, &
         172             :                                        HcoState%NY, &
         173           0 :                                        HcoState%NZ   )
         174             :     REAL(hp), TARGET        :: TmpFlx( HcoState%NX, &
         175             :                                        HcoState%NY, &
         176           0 :                                        HcoState%NZ   )
         177             :     REAL(hp)                :: Mask  ( HcoState%NX, &
         178             :                                        HcoState%NY, &
         179           0 :                                        HcoState%NZ   )
         180             :     REAL(hp)                :: HirFlx( HcoState%NX, &
         181             :                                        HcoState%NY, &
         182           0 :                                        HcoState%NZ   )
         183             :     REAL(hp)                :: HirMsk( HcoState%NX, &
         184             :                                        HcoState%NY, &
         185           0 :                                        HcoState%NZ   )
         186             :
         187             :     ! Integers
         188             :     INTEGER             :: ThisSpc, PrevSpc ! current and previous species ID
         189             :     INTEGER             :: ThisCat, PrevCat ! current and previous category
         190             :     INTEGER             :: ThisHir, PrevHir ! current and previous hierarchy
         191             :     INTEGER             :: SpcMin,  SpcMax  ! range of species to be considered
         192             :     INTEGER             :: CatMin,  CatMax  ! range of categories to be considered
         193             :     INTEGER             :: ExtNr            ! Extension Nr to be used
         194             :     INTEGER             :: nI, nJ, nL
         195             :     INTEGER             :: nnSpec, FLAG
         196             :
         197             :     LOGICAL             :: Found, DoDiagn, EOL, UpdateCat
         198             :
         199             :     ! For error handling & verbose mode
         200             :     CHARACTER(LEN=255)  :: MSG, LOC
         201             :
         202             :     ! testing / debugging
         203             :     integer :: ix,iy
         204             :
         205             :     !=================================================================
         206             :     ! HCO_CalcEmis begins here!
         207             :     !=================================================================
         208             :
         209             :     ! testing only
         210           0 :     ix = 30
         211           0 :     iy = 34
         212             :
         213             :     ! Initialize
         214           0 :     LOC = 'HCO_CalcEmis (HCO_CALC_MOD.F90)'
         215           0 :     Lct => NULL()
         216           0 :     Dct => NULL()
         217             :
         218             :     ! Enter routine
         219           0 :     CALL HCO_ENTER (HcoState%Config%Err, LOC, RC )
         220           0 :     IF(RC /= HCO_SUCCESS) RETURN
         221             :
         222             :     !-----------------------------------------------------------------
         223             :     ! Initialize variables
         224             :     !-----------------------------------------------------------------
         225             :
         226             :     ! Initialize
         227           0 :     SpcFlx(:,:,:)    = 0.0_hp
         228           0 :     CatFlx(:,:,:)    = 0.0_hp
         229           0 :     HirFlx(:,:,:)    = 0.0_hp
         230           0 :     HirMsk(:,:,:)    = 0.0_hp
         231           0 :     PrevSpc          = -1
         232           0 :     PrevHir          = -1
         233           0 :     PrevCat          = -1
         234           0 :     nnSpec           = 0
         235             :
         236             :     ! Pass emission grid dimensions
         237           0 :     nI = HcoState%NX
         238           0 :     nJ = HcoState%NY
         239           0 :     nL = HcoState%NZ
         240             :
         241             :     ! Pass calculation options
         242           0 :     SpcMin  = HcoState%Options%SpcMin        !Lower species ID
         243           0 :     SpcMax  = HcoState%Options%SpcMax        !Upper species ID
         244           0 :     CatMin  = HcoState%Options%CatMin        !Lower emission category
         245           0 :     CatMax  = HcoState%Options%CatMax        !Upper emission category
         246           0 :     ExtNr   = HcoState%Options%ExtNr         !Extension number
         247           0 :     DoDiagn = HcoState%Options%AutoFillDiagn !Write AutoFill diagnostics?
         248             :
         249             :     ! Verbose mode
         250           0 :     IF ( HCO_IsVerb(HcoState%Config%Err ) ) THEN
         251           0 :        WRITE (MSG, *) 'Run HEMCO calculation w/ following options:'
         252           0 :        CALL HCO_MSG ( HcoState%Config%Err, MSG )
         253           0 :        WRITE (MSG, "(A20,I5)")    'Extension number:', ExtNr
         254           0 :        CALL HCO_MSG ( HcoState%Config%Err, MSG )
         255           0 :        WRITE (MSG, "(A20,I5,I5)") 'Tracer range:', SpcMin, SpcMax
         256           0 :        CALL HCO_MSG ( HcoState%Config%Err, MSG )
         257           0 :        WRITE (MSG, "(A20,I5,I5)") 'Category range:', CatMin, CatMax
         258           0 :        CALL HCO_MSG ( HcoState%Config%Err, MSG )
         259           0 :        WRITE (MSG, *) 'Auto diagnostics: ', DoDiagn
         260           0 :        CALL HCO_MSG ( HcoState%Config%Err, MSG )
         261             :     ENDIF
         262             :
         263             :     !=================================================================
         264             :     ! Walk through all containers of EmisList and determine the
         265             :     ! emissions for all containers that qualify for calculation.
         266             :     ! The containers in EmisList are sorted by species, category and
         267             :     ! hierarchy. This enables a straightforward, piece-by-piece
         268             :     ! assembly of the final emission array (start with lowest
         269             :     ! hierarchy emissions, then overwrite piece-by-piece with higher
         270             :     ! hierarchy values).
         271             :     !=================================================================
         272             :
         273             :     ! Point to the head of the emissions linked list
         274           0 :     EOL = .FALSE. ! End of list
         275           0 :     Lct => NULL()
         276           0 :     CALL ListCont_NextCont ( HcoState%EmisList, Lct, FLAG )
         277             :
         278             :     ! Do until end of EmisList (==> loop over all emission containers)
         279             :     DO
         280             :        ! Have we reached the end of the list?
         281           0 :        EOL = .FALSE.
         282           0 :        IF ( FLAG /= HCO_SUCCESS ) EOL = .TRUE.
         283             :
         284             :        ! ------------------------------------------------------------
         285             :        ! Select container and update all working variables & arrays.
         286             :        ! ------------------------------------------------------------
         287             :        IF ( .NOT. EOL ) THEN
         288             :
         289             :           ! Dct is the current data container
         290           0 :           Dct => Lct%Dct
         291             :
         292             :           ! Check if this is a base field
         293           0 :           IF ( Dct%DctType /= HCO_DCTTYPE_BASE ) THEN
         294           0 :              CALL ListCont_NextCont ( HcoState%EmisList, Lct, FLAG )
         295           0 :              CYCLE
         296             :           ENDIF
         297             :
         298             :           ! Sanity check: Make sure this container holds data.
         299             :           ! 'Empty' containers are possible if the simulation time
         300             :           ! is outside of the specified data time range and time
         301             :           ! slice cycling is deactivated (CycleFlag > 1).
         302           0 :           IF( .NOT. FileData_ArrIsDefined(Lct%Dct%Dta) ) THEN
         303           0 :              CALL ListCont_NextCont ( HcoState%EmisList, Lct, FLAG )
         304           0 :              CYCLE
         305             :           ENDIF
         306             :
         307             :           ! Check if this is the specified extension number
         308           0 :           IF ( Dct%ExtNr /= ExtNr ) THEN
         309           0 :              CALL ListCont_NextCont ( HcoState%EmisList, Lct, FLAG )
         310           0 :              CYCLE
         311             :           ENDIF
         312             :
         313             :           ! Advance to next container if the species ID is outside
         314             :           ! the specified species range (SpcMin - SpcMax). Consider
         315             :           ! all species above SpcMin if SpcMax is negative!
         316           0 :           IF( (  Dct%HcoID < SpcMin                     ) .OR. &
         317             :               ( (Dct%HcoID > SpcMax) .AND. (SpcMax > 0) ) ) THEN
         318           0 :              CALL ListCont_NextCont ( HcoState%EmisList, Lct, FLAG )
         319           0 :              CYCLE
         320             :           ENDIF
         321             :
         322             :           ! Advance to next emission field if the emission category of
         323             :           ! the current container is outside of the specified species
         324             :           ! range (CatMin - CatMax). Consider all categories above CatMin
         325             :           ! if CatMax is negative!
         326           0 :           IF( (  Dct%Cat < CatMin                     ) .OR. &
         327             :               ( (Dct%Cat > CatMax) .AND. (CatMax > 0) ) ) THEN
         328           0 :              CALL ListCont_NextCont ( HcoState%EmisList, Lct, FLAG )
         329           0 :              CYCLE
         330             :           ENDIF
         331             :
         332             :           ! Check if this container holds data in the desired unit format,
         333             :           ! i.e. concentration data if UseConc is enabled, emission data
         334             :           ! otherwise.
         335           0 :           IF ( UseConc .NEQV. Dct%Dta%IsConc ) THEN
         336           0 :              CALL ListCont_NextCont ( HcoState%EmisList, Lct, FLAG )
         337           0 :              CYCLE
         338             :           ENDIF
         339             :
         340             :           ! Update working variables
         341           0 :           ThisSpc = Dct%HcoID
         342           0 :           ThisCat = Dct%Cat
         343           0 :           ThisHir = Dct%Hier
         344             :
         345             :        ! If end of list, use dummy values for ThisSpc, ThisCat and ThisHir
         346             :        ! to make sure that emissions are added to HEMCO in the section
         347             :        ! below!
         348             :        ELSE
         349           0 :           ThisSpc = -1
         350           0 :           ThisCat = -1
         351           0 :           ThisHir = -1
         352             :        ENDIF
         353             :
         354             :        !--------------------------------------------------------------------
         355             :        ! Before computing emissions of current data container make sure that
         356             :        ! emissions of previous container are properly archived.
         357             :        !--------------------------------------------------------------------
         358             :
         359             :        ! Add emissions on hierarchy level to the category flux array. Do
         360             :        ! this only if this is a new species, a new category or a new
         361             :        ! hierarchy level.
         362             :        ! Note: no need to add to diagnostics because hierarchy level
         363             :        ! diagnostics are filled right after computing the emissions of
         364             :        ! a given data container (towards the end of the DO loop).
         365             :        IF ( (ThisHir /= PrevHir) .OR. &
         366           0 :             (ThisSpc /= PrevSpc) .OR. &
         367             :             (ThisCat /= PrevCat)        ) THEN
         368             :
         369             :           ! Add hierarchy level emissions to category array over the
         370             :           ! covered regions.
         371           0 :           CatFlx = ( (1.0_hp - HirMsk) * CatFlx ) + HirFlx
         372             :
         373             :           ! Reset
         374           0 :           HirFlx = 0.0_hp
         375           0 :           HirMsk = 0.0_hp
         376             :        ENDIF
         377             :
         378             :        !--------------------------------------------------------------------
         379             :        ! If this is a new species or category, pass the previously collected
         380             :        ! emissions to the species array. Update diagnostics at category level.
         381             :        ! Skip this step for first species, i.e. if PrevSpc is still -1.
         382             :        !--------------------------------------------------------------------
         383           0 :        UpdateCat = .FALSE.
         384           0 :        IF ( ThisCat /= PrevCat ) UpdateCat = .TRUE.
         385           0 :        IF ( ThisSpc /= PrevSpc ) UpdateCat = .TRUE.
         386           0 :        IF ( PrevCat <= 0 .OR. PrevSpc <= 0 ) UpdateCat = .FALSE.
         387           0 :        IF ( UpdateCat ) THEN
         388             :
         389             :           ! CatFlx holds the emissions for this category. Pass this to
         390             :           ! the species array SpcFlx.
         391           0 :           SpcFlx(:,:,:) = SpcFlx(:,:,:) + CatFlx(:,:,:)
         392             :
         393             :           ! verbose
         394           0 :           IF ( HCO_IsVerb(HcoState%Config%Err ) ) THEN
         395           0 :              WRITE(MSG,*) 'Added category emissions to species array: '
         396           0 :              CALL HCO_MSG(HcoState%Config%Err,MSG)
         397           0 :              WRITE(MSG,*) 'Species       : ', PrevSpc
         398           0 :              CALL HCO_MSG(HcoState%Config%Err,MSG)
         399           0 :              WRITE(MSG,*) 'Category      : ', PrevCat
         400           0 :              CALL HCO_MSG(HcoState%Config%Err,MSG)
         401           0 :              WRITE(MSG,*) 'Cat. emissions: ', SUM(CatFlx)
         402           0 :              CALL HCO_MSG(HcoState%Config%Err,MSG)
         403           0 :              WRITE(MSG,*) 'Spc. emissions: ', SUM(SpcFlx)
         404           0 :              CALL HCO_MSG(HcoState%Config%Err,MSG)
         405             :           ENDIF
         406             :
         407             :           ! Add category emissions to diagnostics at category level
         408             :           ! (only if defined in the diagnostics list).
         409           0 :           IF ( Diagn_AutoFillLevelDefined(HcoState%Diagn,3) .AND. DoDiagn ) THEN
         410             :              ! Bug fix: Make sure to pass COL=-1 to ensure all HEMCO diagnostics
         411             :              ! are updated, including those manually defined in other models
         412             :              ! (mps, 11/30/21)
         413             :              CALL Diagn_Update( HcoState,    ExtNr=ExtNr,   &
         414             :                                 Cat=PrevCat, Hier=-1,  HcoID=PrevSpc, &
         415           0 :                                 AutoFill=1,  Array3D=CatFlx, COL=-1, RC=RC )
         416           0 :              IF ( RC /= HCO_SUCCESS ) THEN
         417           0 :                  CALL HCO_ERROR( 'ERROR 0', RC, THISLOC=LOC )
         418           0 :                  RETURN
         419             :              ENDIF
         420             : #ifdef ADJOINT
         421             :              IF (HcoState%IsAdjoint) THEN
         422             :                 CALL Diagn_Update( HcoState,    ExtNr=ExtNr,             &
         423             :                                    Cat=PrevCat, Hier=-1,  HcoID=PrevSpc, &
         424             :                                    AutoFill=1,  Array3D=CatFlx,          &
         425             :                                    COL=HcoState%Diagn%HcoDiagnIDAdjoint, RC=RC )
         426             :                 IF ( RC /= HCO_SUCCESS ) THEN
         427             :                     CALL HCO_ERROR( 'ERROR 1', RC, THISLOC=LOC )
         428             :                     RETURN
         429             :                 ENDIF
         430             :              ENDIF
         431             : #endif
         432             :           ENDIF
         433             :
         434             :           ! Reset CatFlx array and the previously used hierarchy
         435             :           ! ==> Emission hierarchies are only important within the
         436             :           ! same category, hence always start over at lowest hierarchy
         437             :           ! when entering a new category.
         438           0 :           CatFlx(:,:,:)  = 0.0_hp
         439             :           PrevHir        = -1
         440             :        ENDIF
         441             :
         442             :        !--------------------------------------------------------------------
         443             :        ! If this is a new species, pass previously calculated emissions
         444             :        ! to the final emissions array in HcoState.
         445             :        ! Update diagnostics at extension number level.
         446             :        ! Don't do before first emission calculation, i.e. if PrevSpc
         447             :        ! is still the initialized value of -1!
         448             :        !--------------------------------------------------------------------
         449           0 :        IF ( ThisSpc /= PrevSpc .AND. PrevSpc > 0 ) THEN
         450             :
         451             :           ! Add to OutArr
         452           0 :           OutArr(:,:,:) = OutArr(:,:,:) + SpcFlx(:,:,:)
         453             :
         454             :           ! testing only
         455           0 :           IF ( HCO_IsVerb(HcoState%Config%Err ) ) THEN
         456           0 :              WRITE(MSG,*) 'Added total emissions to output array: '
         457           0 :              CALL HCO_MSG(HcoState%Config%Err,MSG)
         458           0 :              WRITE(MSG,*) 'Species: ', PrevSpc
         459           0 :              CALL HCO_MSG(HcoState%Config%Err,MSG)
         460           0 :              WRITE(MSG,*) 'SpcFlx : ', SUM(SpcFlx)
         461           0 :              CALL HCO_MSG(HcoState%Config%Err,MSG)
         462           0 :              WRITE(MSG,*) 'OutArr : ', SUM(OutArr)
         463           0 :              CALL HCO_MSG(HcoState%Config%Err,MSG)
         464             :           ENDIF
         465             :
         466             :           ! Add to diagnostics at extension number level.
         467             :           ! The same diagnostics may be updated multiple times during
         468             :           ! the same time step, continuously adding emissions to it.
         469           0 :           IF ( Diagn_AutoFillLevelDefined(HcoState%Diagn,2) .AND. DoDiagn ) THEN
         470             :              ! Bug fix: Make sure to pass COL=-1 to ensure all HEMCO diagnostics
         471             :              ! are updated, including those manually defined in other models
         472             :              ! (mps, 11/30/21)
         473             :              CALL Diagn_Update( HcoState,  ExtNr=ExtNr,  &
         474             :                                 Cat=-1,    Hier=-1,  HcoID=PrevSpc, &
         475           0 :                                AutoFill=1,Array3D=SpcFlx, COL=-1, RC=RC )
         476           0 :              IF ( RC /= HCO_SUCCESS ) THEN
         477           0 :                  CALL HCO_ERROR( 'ERROR 2', RC, THISLOC=LOC )
         478           0 :                  RETURN
         479             :              ENDIF
         480             : #ifdef ADJOINT
         481             :              IF (HcoState%IsAdjoint) THEN
         482             :                 CALL Diagn_Update( HcoState,  ExtNr=ExtNr,             &
         483             :                                    Cat=-1,    Hier=-1,  HcoID=PrevSpc, &
         484             :                                    AutoFill=1,Array3D=SpcFlx,          &
         485             :                                    COL=HcoState%Diagn%HcoDiagnIDAdjoint, RC=RC )
         486             :                 IF ( RC /= HCO_SUCCESS ) THEN
         487             :                     CALL HCO_ERROR( 'ERROR 3', RC, THISLOC=LOC )
         488             :                     RETURN
         489             :                 ENDIF
         490             :              ENDIF
         491             : #endif
         492             :           ENDIF
         493             :
         494             :           ! Reset arrays and previous hierarchy.
         495           0 :           SpcFlx(:,:,:)  =  0.0_hp
         496           0 :           PrevCat        =  -1
         497           0 :           PrevHir        =  -1
         498           0 :           OutArr         => NULL()
         499             :        ENDIF
         500             :
         501             :        !--------------------------------------------------------------------
         502             :        ! Exit DO loop here if end of list
         503             :        !--------------------------------------------------------------------
         504           0 :        IF ( EOL ) EXIT
         505             :
         506             :        !--------------------------------------------------------------------
         507             :        ! Update/archive information on species level if needed
         508             :        !--------------------------------------------------------------------
         509           0 :        IF ( ThisSpc /= PrevSpc .AND. ThisSpc > 0 ) THEN
         510             :
         511             :           ! Update number of species for which emissions have been
         512             :           ! calculated.
         513           0 :           nnSpec = nnSpec + 1
         514             :
         515             :           ! To write emissions into temporary array, make OutArr point
         516             :           ! to the buffer array HcoState%Buffer3D.
         517           0 :           IF ( HcoState%Options%FillBuffer ) THEN
         518             :
         519             :              ! Cannot use temporary array for more than one species!
         520           0 :              IF ( nnSpec > 1 ) THEN
         521           0 :                 MSG = 'Cannot fill buffer for more than one species!'
         522           0 :                 CALL HCO_ERROR( MSG, RC )
         523           0 :                 RETURN
         524             :              ENDIF
         525             :
         526             :              ! Point to array and check allocation status as well as
         527             :              ! array size.
         528           0 :              OutArr => HcoState%Buffer3D%Val
         529           0 :              IF ( .NOT. ASSOCIATED( OutArr ) ) THEN
         530           0 :                 MSG = 'Buffer array is not associated'
         531           0 :                 CALL HCO_ERROR( MSG, RC )
         532           0 :                 RETURN
         533             :              ENDIF
         534             :              IF ( (SIZE(OutArr,1) /= nI) .OR. &
         535           0 :                   (SIZE(OutArr,2) /= nJ) .OR. &
         536             :                   (SIZE(OutArr,3) /= nL)       ) THEN
         537           0 :                 MSG = 'Buffer array has wrong dimension!'
         538           0 :                 CALL HCO_ERROR( MSG, RC )
         539           0 :                 RETURN
         540             :              ENDIF
         541             :
         542             :           ! To write emissions directly into HcoState, make OutArr
         543             :           ! point to current species' array in HcoState. Use emission
         544             :           ! array for emissions, and concentration array for concentrations.
         545             :           ELSE
         546             :
         547             :              ! For concentrations:
         548           0 :              IF ( UseConc ) THEN
         549           0 :                 CALL HCO_ArrAssert( HcoState%Spc(ThisSpc)%Conc, &
         550           0 :                                     nI, nJ, nL, RC             )
         551           0 :                 IF ( RC /= HCO_SUCCESS ) THEN
         552           0 :                     CALL HCO_ERROR( 'ERROR 4', RC, THISLOC=LOC )
         553           0 :                     RETURN
         554             :                 ENDIF
         555           0 :                 OutArr => HcoState%Spc(ThisSpc)%Conc%Val
         556             :
         557             :              ! For emissions:
         558             :              ELSE
         559           0 :                 CALL HCO_ArrAssert( HcoState%Spc(ThisSpc)%Emis, &
         560           0 :                                     nI, nJ, nL, RC             )
         561           0 :                 IF ( RC /= HCO_SUCCESS ) THEN
         562           0 :                     CALL HCO_ERROR( 'ERROR 5', RC, THISLOC=LOC )
         563           0 :                     RETURN
         564             :                 ENDIF
         565           0 :                 OutArr => HcoState%Spc(ThisSpc)%Emis%Val
         566             :              ENDIF
         567             :
         568             :           ENDIF
         569             :
         570             :           ! verbose mode
         571           0 :           IF ( HCO_IsVerb(HcoState%Config%Err ) ) THEN
         572           0 :              WRITE(MSG,*) 'Calculating emissions for species ', &
         573           0 :                            TRIM(HcoState%Spc(ThisSpc)%SpcName)
         574           0 :              CALL HCO_MSG( HcoState%Config%Err, MSG, SEP1='-', SEP2='-' )
         575             :           ENDIF
         576             :        ENDIF
         577             :
         578             :        !--------------------------------------------------------------------
         579             :        ! Get current emissions and write into TmpFlx array. The array Mask
         580             :        ! denotes all valid grid boxes for this inventory.
         581             :        !--------------------------------------------------------------------
         582           0 :        TmpFlx(:,:,:) = 0.0_hp
         583           0 :        CALL GET_CURRENT_EMISSIONS( HcoState, Dct, nI, nJ, nL, TmpFlx, Mask, RC )
         584           0 :        IF ( RC /= HCO_SUCCESS ) THEN
         585           0 :            CALL HCO_ERROR( 'ERROR 6', RC, THISLOC=LOC )
         586           0 :            RETURN
         587             :        ENDIF
         588             :
         589             :        ! Eventually add universal scale factor
         590           0 :        CALL HCO_ScaleArr( HcoState, ThisSpc, TmpFlx, RC )
         591           0 :        IF ( RC /= HCO_SUCCESS ) THEN
         592           0 :            CALL HCO_ERROR( 'ERROR 7', RC, THISLOC=LOC )
         593           0 :            RETURN
         594             :        ENDIF
         595             :
         596             :        ! Check for negative values according to the corresponding setting
         597             :        ! in the configuration file: 2 means allow negative values, 1 means
         598             :        ! set to zero and prompt a warning, else return with error.
         599           0 :        IF ( HcoState%Options%NegFlag /= 2 ) THEN
         600             :
         601           0 :           IF ( ANY(TmpFlx < 0.0_hp) ) THEN
         602             :
         603             :              ! Set to zero and prompt warning
         604           0 :              IF ( HcoState%Options%NegFlag == 1 ) THEN
         605           0 :                 WHERE ( TmpFlx < 0.0_hp ) TmpFlx = 0.0_hp
         606           0 :                 MSG = 'Negative emissions set to zero: '// TRIM(Dct%cName)
         607           0 :                 CALL HCO_WARNING( HcoState%Config%Err, MSG, RC )
         608             :
         609             :              ! Return with error
         610             :              ELSE
         611             :                 MSG = 'Negative emissions in: '// TRIM(Dct%cName) // '. ' // &
         612           0 :                 'To allow negatives, edit settings in the configuration file.'
         613           0 :                 CALL HCO_ERROR( MSG, RC )
         614           0 :                 RETURN
         615             :              ENDIF
         616             :           ENDIF
         617             :        ENDIF
         618             :
         619             :        ! ------------------------------------------------------------
         620             :        ! Collect all emissions of the same category (and species) on
         621             :        ! the hierarchy level into array HirFlx. HirMsk contains the
         622             :        ! combined covered region. That is, if there are two regional
         623             :        ! inventories with the same hierarchy HirMsk will cover both
         624             :        ! of these regions.
         625             :        ! The specified field hierarchies determine whether the
         626             :        ! temporary emissions are added (if hierarchy is the same
         627             :        ! as the previously used hierarchy), or if they overwrite the
         628             :        ! previous values in HirFlx (if hierarchy is higher than the
         629             :        ! previous hierarchy).
         630             :        ! ------------------------------------------------------------
         631             :
         632             :        ! Add emissions to the hierarchy array HirFlx if this hierarchy
         633             :        ! is the same as previous hierarchy
         634           0 :        IF ( ThisHir == PrevHir ) THEN
         635           0 :           HirFlx = HirFlx + TmpFlx
         636           0 :           HirMsk = HirMsk + Mask
         637             :
         638             :           ! Make sure mask values do not exceed 1.0
         639           0 :           WHERE(HirMsk > 1.0 ) HirMsk = 1.0
         640             :
         641             :        ! If hierarchy is larger than those of the previously used
         642             :        ! fields, overwrite HirFlx with new values.
         643             :        ELSE
         644             :
         645           0 :           HirFlx = TmpFlx
         646           0 :           HirMsk = Mask
         647             :
         648             :        ENDIF
         649             :
         650             :        ! Update diagnostics at hierarchy level. Make sure that only
         651             :        ! positive values are used.
         652             :        ! The same diagnostics may be updated multiple times
         653             :        ! during the same time step, continuously adding
         654             :        ! emissions to it.
         655             :        ! Now remove PosOnly flag. TmpFlx is initialized to zero, so it's
         656             :        ! ok to keep negative values (ckeller, 7/12/15).
         657           0 :        IF ( Diagn_AutoFillLevelDefined(HcoState%Diagn,4) .AND. DoDiagn ) THEN
         658             :              ! Bug fix: Make sure to pass COL=-1 to ensure all HEMCO diagnostics
         659             :              ! are updated, including those manually defined in other models
         660             :              ! (mps, 11/30/21)
         661             :           CALL Diagn_Update( HcoState,       ExtNr=ExtNr,   &
         662             :                              Cat=ThisCat,Hier=ThisHir,   HcoID=ThisSpc, &
         663             :                              AutoFill=1, Array3D=TmpFlx, &
         664           0 :                              COL=-1, RC=RC )
         665           0 :           IF ( RC /= HCO_SUCCESS ) THEN
         666           0 :               CALL HCO_ERROR( 'ERROR 8', RC, THISLOC=LOC )
         667           0 :               RETURN
         668             :           ENDIF
         669             : #ifdef ADJOINT
         670             :           IF (HcoState%IsAdjoint) THEN
         671             :              ! I don't know why I chose collection=-1 instead of
         672             :              ! collection=HcoState%Diagn%HcoDiagnIDAdjoint like in the other
         673             :              ! parts of the adjoint code here, but it's what worked in the
         674             :              ! old repo so I'm keeping it for now. May need to change
         675             :              CALL Diagn_Update( HcoState,       ExtNr=ExtNr,               &
         676             :                                 Cat=ThisCat,Hier=ThisHir,   HcoID=ThisSpc, &
         677             :                                 AutoFill=1, Array3D=TmpFlx,                &
         678             :                                 COL=-1, RC=RC )
         679             :              IF ( RC /= HCO_SUCCESS ) THEN
         680             :                  CALL HCO_ERROR( 'ERROR 9', RC, THISLOC=LOC )
         681             :                  RETURN
         682             :              ENDIF
         683             :           ENDIF
         684             :
         685             : #endif
         686             :        ENDIF
         687             :
         688             :        ! Update previously used species, category and hierarchy
         689           0 :        PrevSpc = ThisSpc
         690           0 :        PrevCat = ThisCat
         691           0 :        PrevHir = ThisHir
         692             :
         693             :        ! Advance to next emission container
         694           0 :        CALL ListCont_NextCont( HcoState%EmisList, Lct, FLAG )
         695             :
         696             :     ENDDO ! Loop over EmisList
         697             :
         698             :     ! Make sure internal pointers are nullified
         699           0 :     Lct    => NULL()
         700           0 :     Dct    => NULL()
         701           0 :     OutArr => NULL()
         702             :
         703             :     ! verbose
         704           0 :     IF ( HCO_IsVerb(HcoState%Config%Err ) ) THEN
         705           0 :        WRITE (MSG, *) 'HEMCO emissions successfully calculated!'
         706           0 :        CALL HCO_MSG ( HcoState%Config%Err, MSG )
         707             :     ENDIF
         708             :
         709             :     ! Leave w/ success
         710           0 :     CALL HCO_LEAVE ( HcoState%Config%Err, RC )
         711             :
         712             :   END SUBROUTINE HCO_CalcEmis
         713             : !EOC
         714             : !------------------------------------------------------------------------------
         715             : !                   Harmonized Emissions Component (HEMCO)                    !
         716             : !------------------------------------------------------------------------------
         717             : !BOP
         718             : !
         719             : ! !IROUTINE: HCO_CheckDepv
         720             : !
         721             : ! !DESCRIPTION: Subroutine HCO\_CheckDepv is a simple routine to check the
         722             : ! dry deposition frequency value. This is to avoid unrealistically high
         723             : ! deposition frequencies that may occur if grid box concentrations are very
         724             : ! low. The deposition frequency is limited to a value that will make sure
         725             : ! that the drydep exponent ( exp( -depfreq * dt ) ) is still small enough to
         726             : ! remove all species mass. The maximum limit of depfreq * dt can be defined
         727             : ! as a HEMCO option (MaxDepExp). Its default value is 20.0.
         728             : !\\
         729             : !\\
         730             : ! !INTERFACE:
         731             : !
         732           0 :   SUBROUTINE HCO_CheckDepv( HcoState, Depv, RC )
         733             : !
         734             : ! !USES:
         735             : !
         736             :     USE HCO_STATE_MOD,    ONLY : HCO_State
         737             : !
         738             : ! !INPUT/OUTPUT PARAMETERS:
         739             : !
         740             :     TYPE(HCO_State), POINTER        :: HcoState   ! HEMCO state object
         741             :     REAL(hp),        INTENT(INOUT)  :: Depv       ! Deposition velocity
         742             :     INTEGER,         INTENT(INOUT)  :: RC         ! Return code
         743             : !
         744             : ! !REVISION HISTORY:
         745             : !  11 Mar 2015 - C. Keller   - Initial Version
         746             : !  See https://github.com/geoschem/hemco for complete history
         747             : !EOP
         748             : !------------------------------------------------------------------------------
         749             : !BOC
         750             : !
         751             : ! !LOCAL VARIABLES:
         752             : !
         753             :     REAL(hp)  :: ExpVal
         754             :
         755             :     !=================================================================
         756             :     ! HCO_CheckDepv begins here!
         757             :     !=================================================================
         758             :
         759           0 :     ExpVal = Depv * HcoState%TS_EMIS
         760           0 :     IF ( ExpVal > HcoState%Options%MaxDepExp ) THEN
         761           0 :        Depv = HcoState%Options%MaxDepExp / HcoState%TS_EMIS
         762             :     ENDIF
         763             :
         764           0 :   END SUBROUTINE HCO_CheckDepv
         765             : !EOC
         766             : !------------------------------------------------------------------------------
         767             : !                   Harmonized Emissions Component (HEMCO)                    !
         768             : !------------------------------------------------------------------------------
         769             : !BOP
         770             : !
         771             : ! !IROUTINE: Get_Current_Emissions
         772             : !
         773             : ! !DESCRIPTION: Subroutine Get\_Current\_Emissions calculates the current
         774             : !  emissions for the specified emission container.
         775             : !  This subroutine is only called by HCO\_CalcEmis and for base emission
         776             : !  containers, i.e. containers of type 1.
         777             : !\\
         778             : !\\
         779             : ! !INTERFACE:
         780             : !
         781           0 :   SUBROUTINE Get_Current_Emissions( HcoState, BaseDct,   nI,   nJ,           &
         782           0 :                                     nL,       OUTARR_3D, MASK, RC, UseLL    )
         783             : !
         784             : ! !USES:
         785             : !
         786             :     USE HCO_State_Mod,    ONLY : HCO_State
         787             :     USE HCO_tIdx_MOD,     ONLY : tIDx_GetIndx
         788             :     USE HCO_FileData_Mod, ONLY : FileData_ArrIsDefined
         789             : !
         790             : ! !INPUT PARAMETERS:
         791             : !
         792             :     INTEGER,           INTENT(IN)  :: nI                  ! # of lons
         793             :     INTEGER,           INTENT(IN)  :: nJ                  ! # of lats
         794             :     INTEGER,           INTENT(IN)  :: nL                  ! # of levs
         795             : !
         796             : ! !INPUT/OUTPUT PARAMETERS:
         797             : !
         798             :
         799             :     TYPE(HCO_State), POINTER       :: HcoState            ! HEMCO state object
         800             :     TYPE(DataCont),  POINTER       :: BaseDct             ! base emission
         801             :                                                           !  container
         802             :     REAL(hp),        INTENT(INOUT) :: outArr_3D(nI,nJ,nL) ! output array
         803             :     REAL(hp),        INTENT(INOUT) :: mask     (nI,nJ,nL) ! mask array
         804             :     INTEGER,         INTENT(INOUT) :: RC
         805             : !
         806             : ! !OUTPUT PARAMETERS:
         807             : !
         808             :     INTEGER,         INTENT(  OUT), OPTIONAL :: UseLL
         809             : !
         810             : ! !REMARKS:
         811             : !  This routine uses multiple loops over all grid boxes (base emissions
         812             : !  and scale factors use separate loops). In an OMP environment, this approach
         813             : !  seems to be faster than using only one single loop (but repeated calls to
         814             : !  point to containers, etc.). The alternative approach is used in routine
         815             : !  Get\_Current\_Emissions\_B at the end of this module and may be employed
         816             : !  on request.
         817             : !
         818             : ! !REVISION HISTORY:
         819             : !  25 Aug 2012 - C. Keller   - Initial Version
         820             : !  See https://github.com/geoschem/hemco for complete history
         821             : !EOP
         822             : !------------------------------------------------------------------------------
         823             : !BOC
         824             : !
         825             : ! !LOCAL VARIABLES:
         826             : !
         827             :     ! Pointers
         828             :     TYPE(DataCont), POINTER :: ScalDct
         829             :     TYPE(DataCont), POINTER :: MaskDct
         830             :     TYPE(DataCont), POINTER :: LevDct1
         831             :     TYPE(DataCont), POINTER :: LevDct2
         832             :
         833             :     ! Scalars
         834             :     REAL(sp)                :: tmpVal, MaskScale
         835             :     REAL(hp)                :: outData
         836             :     REAL(hp)                :: ScalFact
         837             :     INTEGER                 :: tIDx,  IDX
         838             :     INTEGER                 :: totLL, nnLL
         839             :     INTEGER                 :: I,     J,     L,      N
         840             :     INTEGER                 :: LowLL, UppLL, ScalLL, TmpLL
         841             :     INTEGER                 :: EC,    ERROR
         842             :     CHARACTER(LEN=255)      :: MSG,   thisLoc
         843             :     LOGICAL                 :: NegScalExist
         844             :     LOGICAL                 :: MaskFractions
         845             :     LOGICAL                 :: isLevDct1
         846             :     LOGICAL                 :: isLevDct2
         847             :     LOGICAL                 :: isMaskDct
         848             :     LOGICAL                 :: isPblHt
         849             :     LOGICAL                 :: isBoxHt
         850             :     INTEGER                 :: LevDct1_Unit
         851             :     INTEGER                 :: LevDct2_Unit
         852             :
         853             :     ! testing only
         854             :     INTEGER, PARAMETER      :: IX=25, IY=25
         855             :
         856             :     !=================================================================
         857             :     ! GET_CURRENT_EMISSIONS begins here
         858             :     !=================================================================
         859             :
         860             :     ! Initialize
         861           0 :     ScalDct => NULL()
         862           0 :     MaskDct => NULL()
         863           0 :     msg     = ''
         864           0 :     thisLoc = 'GET_CURRENT_EMISSIONS (hco_calc_mod.F90)'
         865             :
         866             :     ! Enter
         867           0 :     CALL HCO_Enter( HcoState%Config%Err, thisLoc, RC )
         868           0 :     IF ( RC /= HCO_SUCCESS ) RETURN
         869             :
         870             :     ! Check if container contains data
         871           0 :     IF ( .NOT. FileData_ArrIsDefined(BaseDct%Dta) ) THEN
         872           0 :        msg = 'Array not defined: ' // TRIM(BaseDct%cName)
         873           0 :        CALL HCO_Error( msg, RC, thisLoc )
         874           0 :        RETURN
         875             :     ENDIF
         876             :
         877             :     ! Initialize mask. By default, assume that we use all grid boxes.
         878           0 :     MASK          = 1.0_hp
         879           0 :     MaskFractions = HcoState%Options%MaskFractions
         880             :
         881             :     ! Verbose output
         882           0 :     IF ( HCO_IsVerb(HcoState%Config%Err ) ) THEN
         883           0 :        WRITE(msg,*) 'Evaluate field ', TRIM(BaseDct%cName)
         884           0 :        CALL HCO_Msg( HcoState%Config%Err, msg, sep1=' ' )
         885             :     ENDIF
         886             :
         887             : #if !defined ( ESMF_ )
         888             :     ! Put check for PBLHEIGHT here, if not running in ESMF/MAPL
         889           0 :     IF ( .NOT. ASSOCIATED( HcoState%Grid%PBLHEIGHT%Val ) ) THEN
         890           0 :        msg = 'PBLHEIGHT (in meters) is missing in HEMCO state'
         891           0 :        CALL HCO_Error( msg, RC, thisLoc )
         892           0 :        RETURN
         893             :     ENDIF
         894             : #endif
         895             :
         896             :     !========================================================================
         897             :     ! Set base emissions
         898             :     !========================================================================
         899             :
         900             :     ! Check for level index containers
         901             :     ! Move error checks here, outside of the parallel DO loop
         902             :     ! Remove ELSE blocks for efficiency
         903           0 :     LevDct1 => NULL()
         904           0 :     IF ( BaseDct%levScalID1 > 0 ) THEN
         905           0 :        CALL Pnt2DataCont( HcoState, BaseDct%levScalID1, LevDct1, RC )
         906           0 :        IF ( RC /= HCO_SUCCESS ) THEN
         907           0 :           msg = 'Could not get a pointer to LevDct1 from HcoState!'
         908           0 :           CALL HCO_Error( msg, RC, thisLoc )
         909           0 :           RETURN
         910             :        ENDIF
         911             :     ENDIF
         912             :
         913           0 :     LevDct2 => NULL()
         914           0 :     IF ( BaseDct%levScalID2 > 0 ) THEN
         915           0 :        CALL Pnt2DataCont( HcoState, BaseDct%levScalID2, LevDct2, RC )
         916           0 :        IF ( RC /= HCO_SUCCESS ) THEN
         917           0 :           msg = 'Could not get a pointer to LevDct2 in HcoState!'
         918           0 :           CALL HCO_Error( msg, RC, thisLoc )
         919           0 :           RETURN
         920             :        ENDIF
         921             :     ENDIF
         922             :
         923             :     ! Test whether LevDct1 and LevDct2 are associated
         924           0 :     isLevDct1 = ASSOCIATED( LevDct1 )
         925           0 :     isLevDct2 = ASSOCIATED( LevDct2 )
         926             :
         927             :     ! Get the units of LevDct1 (if it exists)
         928           0 :     LevDct1_Unit = -1
         929           0 :     IF ( isLevDct1 ) THEN
         930           0 :        LevDct1_Unit = GetEmisLUnit( HcoState, LevDct1 )
         931           0 :        IF ( LevDct1_Unit < 0 ) THEN
         932           0 :           MSG = 'LevDct1 units are not defined!'
         933           0 :           CALL HCO_Error( msg, RC, thisLoc )
         934           0 :           RETURN
         935             :        ENDIF
         936             :     ENDIF
         937             :
         938             :     ! Get the units of LevDct2 (if it exists)
         939           0 :     LevDct2_Unit = -1
         940           0 :     IF ( isLevDct2 ) THEN
         941           0 :        LevDct2_Unit = GetEmisLUnit( HcoState, LevDct2 )
         942           0 :        IF ( LevDct2_Unit < 0 ) THEN
         943           0 :           MSG = 'LevDct2_Units are not defined!'
         944           0 :           CALL HCO_Error( msg, RC, thisLoc )
         945           0 :           RETURN
         946             :        ENDIF
         947             :     ENDIF
         948             :
         949             :     ! Throw an error if boxheight is missing and the units are in meters
         950           0 :     IF ( LevDct1_Unit == HCO_EMISL_M  .or.                                  &
         951             :          LevDct2_Unit == HCO_EMISL_M ) THEN
         952           0 :        IF ( .not. ASSOCIATED( HcoState%Grid%BXHEIGHT_M%Val ) ) THEN
         953           0 :           MSG = 'Boxheight (in meters) is missing in HEMCO state'
         954           0 :           CALL HCO_Error( msg, RC, thisLoc )
         955           0 :           RETURN
         956             :        ENDIF
         957             :     ENDIF
         958             :
         959             :     ! Throw an error if boxheight is missing and the units are in PBL frac
         960           0 :     IF ( LevDct1_Unit == HCO_EMISL_PBL  .or.                                &
         961             :          LevDct2_Unit == HCO_EMISL_PBL ) THEN
         962           0 :        IF ( .not. ASSOCIATED( HcoState%Grid%PBLHEIGHT%Val ) ) THEN
         963           0 :           MSG = 'Boundary layer height is missing in HEMCO state'
         964           0 :           CALL HCO_Error( msg, RC, thisLoc )
         965           0 :           RETURN
         966             :        ENDIF
         967             :     ENDIF
         968             :
         969             :     ! Initialize non-private loop variables here
         970           0 :     error = 0
         971           0 :     totLL = 0.0
         972           0 :     nnLL  = 0.0
         973             :
         974             :     ! Loop over all latitudes and longitudes
         975             :     !$OMP PARALLEL DO                                                        &
         976             :     !$OMP DEFAULT( SHARED                                                   )&
         977             :     !$OMP PRIVATE( I, J, L, tIdx, tmpVal, LowLL, UppLL, EC                  )&
         978             :     !$OMP REDUCTION( +:totLL                                                )&
         979             :     !$OMP REDUCTION( +:nnLL                                                 )&
         980             :     !$OMP COLLAPSE( 2                                                       )&
         981             :     !$OMP SCHEDULE( DYNAMIC, 8                                              )
         982           0 :     DO J = 1, nJ
         983           0 :     DO I = 1, nI
         984             :
         985             :        ! Continue to end of loop if an error has occurred
         986             :        ! (we cannot exit from a parallel loop)
         987           0 :        IF ( error > 0 ) CYCLE
         988             :
         989             :        ! Zero private variables for safety's sake
         990           0 :        tmpVal  = 0.0_hp
         991           0 :        lowLL   = 0
         992           0 :        uppLL   = 0
         993           0 :        EC      = HCO_SUCCESS
         994             :
         995             :        !---------------------------------------------------------------------
         996             :        ! Get current time index for this container and at this location
         997             :        !---------------------------------------------------------------------
         998           0 :        tIDx = tIDx_GetIndx( HcoState, BaseDct%Dta, I, J )
         999           0 :        IF ( tIDx < 1 ) THEN
        1000           0 :           WRITE( msg, * ) 'Cannot get time slice index at location ',I,J,    &
        1001           0 :                           ': ', TRIM(BaseDct%cName), tIDx
        1002           0 :           error = 1
        1003           0 :           CYCLE
        1004             :        ENDIF
        1005             :
        1006             :        !---------------------------------------------------------------------
        1007             :        ! Get lower and upper vertical index
        1008             :        !---------------------------------------------------------------------
        1009             :        CALL GetVertIndx( HcoState,     BaseDct,   isLevDct1, LevDct1,        &
        1010             :                          LevDct1_Unit, isLevDct2, LevDct2,   LevDct2_Unit,   &
        1011             :                          I,            J,         LowLL,     UppLL,          &
        1012           0 :                          RC=EC                                              )
        1013             :
        1014             :        ! Trap error
        1015           0 :        IF ( EC /= HCO_SUCCESS ) THEN
        1016           0 :           WRITE(MSG,*) 'Error getting vertical index at location ',I,J,&
        1017           0 :                        ': ', TRIM(BaseDct%cName)
        1018           0 :           error = 1
        1019           0 :           CYCLE
        1020             :        ENDIF
        1021             :
        1022             :        !---------------------------------------------------------------------
        1023             :        ! Apply vertical dilution factor (if necessary)
        1024             :        !---------------------------------------------------------------------
        1025             :
        1026             :        ! Update variables for computing the average level
        1027           0 :        totLL = totLL + UppLL
        1028           0 :        nnLL  = nnLL  + 1
        1029             :
        1030             :        ! Loop over all levels
        1031           0 :        DO L = LowLL, UppLL
        1032             :
        1033             :           ! Get base value. Use uniform value if scalar field.
        1034           0 :           tmpVal = Get_Value_From_DataCont( I, J, L, tIdx, BaseDct )
        1035             :
        1036             :           ! If it's a missing value, mask box as unused
        1037             :           ! and set value to zero.
        1038           0 :           IF ( tmpVal == HCO_MISSVAL ) THEN
        1039           0 :              mask(I,J,:)      = 0.0_hp
        1040           0 :              outArr_3D(I,J,L) = 0.0_hp
        1041           0 :              CYCLE
        1042             :           ENDIF
        1043             :
        1044             :           ! Otherwise, apply the vertical dilution factor (if necessary)
        1045             :           CALL Apply_Dilution_Factor( I        = I,                          &
        1046             :                                       J        = J,                          &
        1047             :                                       L        = L,                          &
        1048             :                                       LowLL    = LowLL,                      &
        1049             :                                       UppLL    = UppLL,                      &
        1050             :                                       HcoState = HcoState,                   &
        1051             :                                       BaseDct  = BaseDct,                    &
        1052             :                                       inData   = tmpVal,                     &
        1053           0 :                                       outData  = outArr_3D(I,J,L),           &
        1054           0 :                                       RC       = EC                         )
        1055             :
        1056             :           ! Trap errors
        1057           0 :           IF ( EC /= HCO_SUCCESS ) THEN
        1058           0 :              error = 1
        1059           0 :              msg   = 'Error encountered in routine "Apply_Dilution_Factor"!'
        1060           0 :              CYCLE
        1061             :           ENDIF
        1062             :
        1063             :        ENDDO !L
        1064             :
        1065             :     ENDDO !I
        1066             :     ENDDO !J
        1067             :     !$OMP END PARALLEL DO
        1068             :
        1069             :     !------------------------------------------------------------------------
        1070             :     ! Return if an error occurred in the parallel loop above
        1071             :     !------------------------------------------------------------------------
        1072           0 :     IF ( error == 1 ) THEN
        1073           0 :        CALL HCO_Error( msg, RC, thisLoc )
        1074           0 :        RETURN
        1075             :     ENDIF
        1076             :
        1077             :     !========================================================================
        1078             :     ! Apply scale factors to base emissions
        1079             :     !
        1080             :     ! The container IDs of all scale factors associated with this base
        1081             :     ! container are stored in vector Scal_cID.
        1082             :     !========================================================================
        1083           0 :     IF ( BaseDct%nScalID > 0 ) THEN
        1084             :
        1085             :        ! Loop over all scale factors for this base emissions
        1086           0 :        DO N = 1, BaseDct%nScalID
        1087             :
        1088             :           !------------------------------------------------------------------
        1089             :           ! Get a pointer to the data container holding this scale factor
        1090             :           !------------------------------------------------------------------
        1091             :
        1092             :           ! Get the Nth scale factor container ID
        1093           0 :           IDX = BaseDct%Scal_cID(N)
        1094             :
        1095             :           ! Point to data container with the given container ID
        1096           0 :           CALL Pnt2DataCont( HcoState, IDX, ScalDct, RC )
        1097           0 :           IF ( RC /= HCO_SUCCESS ) THEN
        1098             :              msg = 'Could not get scale factor for base emissions field ' // &
        1099           0 :                    TRIM( BaseDct%cName )
        1100           0 :              CALL HCO_Error( msg, RC, thisLoc )
        1101           0 :              RETURN
        1102             :           ENDIF
        1103             :
        1104             :           ! Sanity check: scale field cannot be a base field
        1105           0 :           IF ( ScalDct%DctType == HCO_DCTTYPE_BASE ) THEN
        1106           0 :              msg = 'Wrong scale field type: ' // TRIM(ScalDct%cName)
        1107           0 :              CALL HCO_Error( msg, RC, thisLoc )
        1108           0 :              RETURN
        1109             :           ENDIF
        1110             :
        1111             :           !------------------------------------------------------------------
        1112             :           ! Skip this scale factor if no data defined. This is possible
        1113             :           ! if scale factors are only defined for a given time range and
        1114             :           ! the simulation datetime is outside of this range.
        1115             :           !------------------------------------------------------------------
        1116           0 :           IF ( .not. FileData_ArrIsDefined( ScalDct%Dta ) ) THEN
        1117           0 :              IF ( HCO_IsVerb( HcoState%Config%Err ) ) THEN
        1118             :                 msg = 'Skip scale factor ' // TRIM( ScalDct%cName )       // &
        1119           0 :                      ' because it is not defined for this datetime.'
        1120           0 :                 CALL HCO_MSG( HcoState%Config%Err, msg )
        1121             :              ENDIF
        1122             :              CYCLE
        1123             :           ENDIF
        1124             :
        1125             :           ! Verbose printout
        1126           0 :           IF ( HCO_IsVerb(HcoState%Config%Err ) ) THEN
        1127           0 :              MSG = 'Applying scale factor ' // TRIM(ScalDct%cName)
        1128           0 :              CALL HCO_MSG(HcoState%Config%Err,MSG)
        1129             :           ENDIF
        1130             :
        1131             :           !------------------------------------------------------------------
        1132             :           ! Get vertical extension of this scale factor array.
        1133             :           !------------------------------------------------------------------
        1134           0 :           ScalLL = 1
        1135           0 :           IF ( ScalDct%Dta%SpaceDim == 3 ) THEN
        1136           0 :              ScalLL = SIZE( ScalDct%Dta%V3(1)%Val, 3 )
        1137             :           ENDIF
        1138             :
        1139             :           !------------------------------------------------------------------
        1140             :           ! Check if there is a mask field associated with this scale
        1141             :           ! factor. In this case, get a pointer to the corresponding
        1142             :           ! mask field and evaluate scale factors only inside the mask
        1143             :           ! region.
        1144             :           !------------------------------------------------------------------
        1145           0 :           IF ( ASSOCIATED(ScalDct%Scal_cID) ) THEN
        1146           0 :              CALL Pnt2DataCont( HcoState, ScalDct%Scal_cID(1), MaskDct, RC )
        1147           0 :              IF ( RC /= HCO_SUCCESS ) THEN
        1148             :                 msg = 'Could not get mask field for scale factor: ' //       &
        1149           0 :                       TRIM( ScalDct%cName )
        1150           0 :                 CALL HCO_Error( msg, RC, thisLoc )
        1151           0 :                 RETURN
        1152             :              ENDIF
        1153             :
        1154             :              ! Sanity check: The data container MaskDct must be a mask!
        1155           0 :              IF ( MaskDct%DctType /= HCO_DCTTYPE_MASK ) THEN
        1156             :                 MSG = 'Invalid mask for scale factor: '                   // &
        1157             :                       TRIM( ScalDct%cName ) // '; mask: '                 // &
        1158           0 :                       TRIM( MaskDct%cName )
        1159           0 :                 CALL HCO_Error( msg, RC, thisLoc )
        1160           0 :                 RETURN
        1161             :              ENDIF
        1162             :           ENDIF
        1163             :
        1164             :           ! Set a flag to denote whether MaskDct is associated
        1165             :           ! This can be done outside of the parallel loops below
        1166           0 :           isMaskDct = ASSOCIATED( MaskDct )
        1167             :
        1168             :           !------------------------------------------------------------------
        1169             :           ! Apply the mask to the scale factor
        1170             :           !------------------------------------------------------------------
        1171             :
        1172             :           ! Reinitialize error flag. Will be set to > 0 if error occurs,
        1173             :           ! and to -1 if negative scale factor is ignored.
        1174           0 :           error = 0
        1175             :
        1176             :           ! Loop over all latitudes and longitudes
        1177             :           !$OMP PARALLEL DO                                                  &
        1178             :           !$OMP DEFAULT( SHARED                                             )&
        1179             :           !$OMP PRIVATE( I,     J,     tIdx,  tmpVal,    L                  )&
        1180             :           !$OMP PRIVATE( LowLL, UppLL, tmpLL, MaskScale, EC                 )&
        1181             :           !$OMP COLLAPSE( 2                                                 )&
        1182             :           !$OMP SCHEDULE( DYNAMIC, 8                                        )
        1183           0 :           DO J = 1, nJ
        1184           0 :           DO I = 1, nI
        1185             :
        1186             :              ! Continue to end of loop if an error has occurred
        1187             :              ! (we cannot exit from a parallel loop)
        1188           0 :              IF ( error > 0 ) CYCLE
        1189             :
        1190             :              !---------------------------------------------------------------
        1191             :              ! If there is a mask associated with this scale factors, check
        1192             :              ! if this grid box is within or outside of the mask region.
        1193             :              ! Values that partially fall into the mask region are either
        1194             :              ! treated as binary (100% inside or outside), or partially
        1195             :              ! (using the real grid area fractions), depending on the
        1196             :              ! HEMCO options.
        1197             :              !---------------------------------------------------------------
        1198             :
        1199             :              ! Default mask scaling is 1.0 (no mask applied)
        1200           0 :              maskScale = 1.0_sp
        1201             :
        1202             :              ! If there is a mask applied to this scale factor ...
        1203           0 :              IF ( isMaskDct ) THEN
        1204             :                 CALL GetMaskVal ( MaskDct,   I,             J,               &
        1205           0 :                                   MaskScale, MaskFractions, EC              )
        1206           0 :                 IF ( EC /= HCO_SUCCESS ) THEN
        1207             :                    error = 4
        1208             :                    CYCLE
        1209             :                 ENDIF
        1210             :              ENDIF
        1211             :
        1212             :              ! We continue an skip this grid box if mask is completely zero
        1213           0 :              IF ( maskScale <= 0.0_sp ) CYCLE
        1214             :
        1215             :              ! Get current time index for this container and at this location
        1216           0 :              tIDx = tIDx_GetIndx( HcoState, ScalDct%Dta, I, J )
        1217           0 :              IF ( tIDx < 1 ) THEN
        1218           0 :                 WRITE(*,*) 'Cannot get time slice index at location ',I,J,   &
        1219           0 :                            ': ', TRIM(ScalDct%cName), tIDx
        1220           0 :                 error = 3
        1221           0 :                 CYCLE
        1222             :              ENDIF
        1223             :
        1224             :              !---------------------------------------------------------------
        1225             :              ! Check if this is a mask. If so, add mask values to the MASK
        1226             :              ! array. For now, we assume masks to be binary, i.e. 0 or 1.
        1227             :              ! We may want to change that in future to also support values
        1228             :              ! in between. This is especially important when regridding
        1229             :              ! high resolution masks onto coarser grids!
        1230             :              !---------------------------------------------------------------
        1231           0 :              IF ( ScalDct%DctType == HCO_DCTTYPE_MASK ) THEN
        1232             :
        1233             :                 ! Get mask value
        1234           0 :                 CALL GetMaskVal( ScalDct, I, J, TMPVAL, MaskFractions, EC )
        1235           0 :                 IF ( EC /= HCO_SUCCESS ) THEN
        1236             :                    error = 4
        1237             :                    CYCLE
        1238             :                 ENDIF
        1239             :
        1240             :                 ! Pass to output mask
        1241           0 :                 mask(I,J,:) = mask(I,J,:) * TMPVAL
        1242             :
        1243             :                 ! Verbose printout
        1244             :                 IF ( HCO_IsVerb(HcoState%Config%Err) .and.                   &
        1245           0 :                      I==1 .AND. J==1                       ) THEN
        1246             :                    msg = 'Mask field ' // TRIM( ScalDct%cName )           // &
        1247           0 :                          ' found and added to temporary mask.'
        1248           0 :                    CALL HCO_Msg( HcoState%Config%Err, msg )
        1249             :                 ENDIF
        1250             :
        1251             :                 ! Advance to next grid box
        1252             :                 CYCLE
        1253             :              ENDIF
        1254             :
        1255             :              !---------------------------------------------------------------
        1256             :              ! For non-mask fields, apply scale factors to all levels
        1257             :              ! of the base field individually. If the scale factor
        1258             :              ! field has more than one vertical level, use the
        1259             :              ! vertical level closest to the corresponding vertical
        1260             :              ! level of the base emission field
        1261             :              !---------------------------------------------------------------
        1262             :
        1263             :              ! Get lower and upper vertical index
        1264             :              CALL GetVertIndx( HcoState, BaseDct,       isLevDct1,           &
        1265             :                                LevDct1,  LevDct1_Unit,  isLevDct2,           &
        1266             :                                LevDct2,  LevDct2_Unit,  I,                   &
        1267           0 :                                J,        LowLL,         UppLL,     EC       )
        1268             :
        1269             :              ! Trap errors
        1270           0 :              IF ( EC /= HCO_SUCCESS ) THEN
        1271             :                 error = 1
        1272             :                 CYCLE
        1273             :              ENDIF
        1274             :
        1275             :              ! Loop over all vertical levels of the base field
        1276           0 :              DO L = LowLL,UppLL
        1277             :
        1278             :                 ! If the vertical level exceeds the number of available
        1279             :                 ! scale factor levels, use the highest available level.
        1280             :                 ! Otherwise use the same vertical level index.
        1281           0 :                 TmpLL = L
        1282           0 :                 IF ( L > ScalLL ) TmpLL = ScalLL
        1283             :
        1284             :                 !------------------------------------------------------------
        1285             :                 ! Get scale factor for this grid box. Use same uniform
        1286             :                 ! value if it's a scalar field.
        1287             :                 !------------------------------------------------------------
        1288           0 :                 tmpVal = Get_Value_From_DataCont( I, J, L, tIdX, ScalDct )
        1289             :
        1290             :                 ! Set missing value to one
        1291           0 :                 IF ( tmpVal == HCO_MISSVAL ) tmpVal = 1.0_sp
        1292             :
        1293             :                 ! Eventually apply mask scaling
        1294           0 :                 IF ( maskScale /= 1.0_sp ) THEN
        1295           0 :                    tmpVal = tmpVal * MaskScale
        1296             :                 ENDIF
        1297             :
        1298             :                 !------------------------------------------------------------
        1299             :                 ! Negative scale factors: proceed according to the "negative
        1300             :                 ! value" setting specified in the HEMCO configuration file
        1301             :                 ! NegFlag = 2: Use this value (i.e. pass thru this IF stmt)
        1302             :                 ! NegFlag = 1: Ignore ands how warning
        1303             :                 ! NegFlag = 0: Return with error
        1304             :                 !------------------------------------------------------------
        1305           0 :                 IF ( tmpVal < 0.0_sp .and. HcoState%Options%NegFlag /= 2 ) THEN
        1306             :
        1307             :                    ! NegFlag = 1: ignore and show warning
        1308             :                    ! Otherwise return with error
        1309           0 :                    IF ( HcoState%Options%NegFlag == 1 ) THEN
        1310             :                       error = -1
        1311             :                       CYCLE
        1312             :                    ELSE
        1313           0 :                       WRITE( 6, * ) 'Negative scale factor at ',             &
        1314           0 :                                      I, J, TmpLL, tidx, ': ',                &
        1315           0 :                                      TRIM(ScalDct%cName), TMPVAL
        1316           0 :                       error = 1
        1317           0 :                       CYCLE
        1318             :                    ENDIF
        1319             :                 ENDIF
        1320             :
        1321             :                 !------------------------------------------------------------
        1322             :                 ! Apply scale factor in accordance to field operator
        1323             :                 !------------------------------------------------------------
        1324             :                 CALL Apply_Scale_Factor( I       = I,                        &
        1325             :                                          J       = J,                        &
        1326             :                                          L       = L,                        &
        1327             :                                          ScalDct = ScalDct,                  &
        1328             :                                          scalFac = tmpVal,                   &
        1329           0 :                                          dataVal = outArr_3D(I,J,L),         &
        1330           0 :                                          RC      = EC                       )
        1331             :
        1332             :                 ! Trap errors
        1333           0 :                 IF ( EC /= HCO_SUCCESS ) THEN
        1334           0 :                    error = 2
        1335           0 :                    CYCLE
        1336             :                 ENDIF
        1337             :
        1338             :              ENDDO
        1339             :
        1340             :              !---------------------------------------------------------------
        1341             :              ! Verbose printout
        1342             :              !----------------------------------------------------------------
        1343             :              IF ( HCO_IsVerb(HcoState%Config%Err)  .and.                     &
        1344           0 :                   I == ix .and. J == iy           ) THEN
        1345           0 :                 write(MSG,*) 'Scale field ', TRIM(ScalDct%cName)
        1346           0 :                 CALL HCO_MSG(HcoState%Config%Err,MSG)
        1347           0 :                 write(MSG,*) 'Time slice: ', tIdx
        1348           0 :                 CALL HCO_MSG(HcoState%Config%Err,MSG)
        1349           0 :                 write(MSG,*) 'IX, IY: ', IX, IY
        1350           0 :                 CALL HCO_MSG(HcoState%Config%Err,MSG)
        1351           0 :                 write(MSG,*) 'Scale factor (IX,IY,L1): ', TMPVAL
        1352           0 :                 CALL HCO_MSG(HcoState%Config%Err,MSG)
        1353           0 :                 write(MSG,*) 'Mathematical operation : ', ScalDct%Oper
        1354           0 :                 CALL HCO_MSG(HcoState%Config%Err,MSG)
        1355             :              ENDIF
        1356             :
        1357             :           ENDDO !I
        1358             :           ENDDO !J
        1359             :           !$OMP END PARALLEL DO
        1360             :
        1361             :           !------------------------------------------------------------------
        1362             :           ! Return with error if an error was encountered in the loop above
        1363             :           !------------------------------------------------------------------
        1364           0 :           IF ( error > 0 ) THEN
        1365             :
        1366             :              ! Construct the appropriate error message
        1367           0 :              SELECT CASE( error )
        1368             :                 CASE( 1 )
        1369           0 :                    msg = 'Negative scale factor found (aborted): '
        1370             :                 CASE( 2 )
        1371           0 :                    msg = 'Illegal mathematical operator for scale factor: '
        1372             :                 CASE( 3 )
        1373           0 :                    msg = 'Encountered negative time index for scale factor: '
        1374             :                 CASE( 4 )
        1375           0 :                    msg = 'Error applying mask to scale factor: '
        1376             :                 CASE DEFAULT
        1377           0 :                    msg = 'Error when applying scale factor: '
        1378             :              END SELECT
        1379             :
        1380             :              ! Append name of scale factor to error message
        1381           0 :              msg = TRIM( msg ) // TRIM( ScalDct%cName )
        1382             :
        1383             :              ! Exit with error message
        1384           0 :              ScalDct => NULL()
        1385           0 :              CALL HCO_Error( msg, RC, thisLoc )
        1386           0 :              RETURN
        1387             :           ENDIF
        1388             :
        1389             :           ! eventually prompt warning for negative values
        1390           0 :           IF ( ERROR == -1 ) THEN
        1391             :              msg = 'Negative scale factor found (ignored): '              // &
        1392           0 :                     TRIM( ScalDct%cName )
        1393           0 :              CALL HCO_WARNING( HcoState%Config%Err, msg, RC )
        1394             :           ENDIF
        1395             :
        1396             :           ! Free pointer
        1397           0 :           MaskDct => NULL()
        1398             :
        1399             :        ENDDO ! N
        1400             :     ENDIF    ! N > 0
        1401             :
        1402             :     !========================================================================
        1403             :     ! Update optional variables
        1404             :     !========================================================================
        1405           0 :     IF ( PRESENT( UseLL) ) THEN
        1406           0 :        UseLL = 1
        1407           0 :        IF ( nnLL > 0 ) UseLL = NINT(REAL(TotLL,kind=sp)/REAL(nnLL,kind=sp))
        1408             :     ENDIF
        1409             :
        1410             :     ! Weight output emissions by mask
        1411           0 :     outArr_3D = outArr_3D * mask
        1412             :
        1413             :     ! Cleanup and leave w/ success
        1414           0 :     ScalDct => NULL()
        1415           0 :     CALL HCO_Leave( HcoState%Config%Err, RC )
        1416             :
        1417             :   END SUBROUTINE Get_Current_Emissions
        1418             : !EOC
        1419             : !------------------------------------------------------------------------------
        1420             : !                  GEOS-Chem Global Chemical Transport Model                  !
        1421             : !------------------------------------------------------------------------------
        1422             : !BOP
        1423             : !
        1424             : ! !IROUTINE: Get_Value_From_DataCont
        1425             : !
        1426             : ! !DESCRIPTION: Returns a data value stored in a data container at a given
        1427             : !  grid box and time.
        1428             : !\\
        1429             : !\\
        1430             : ! !INTERFACE:
        1431             : !
        1432           0 :   FUNCTION Get_Value_From_DataCont( I, J, L, tIdx, Dct ) RESULT( dataVal )
        1433             : !
        1434             : ! !INPUT PARAMETERS:
        1435             : !
        1436             :     INTEGER,        INTENT(IN) :: I        ! Longitude (or X-dim) index
        1437             :     INTEGER,        INTENT(IN) :: J        ! Latitude  (or Y-dim) index
        1438             :     INTEGER,        INTENT(IN) :: L        ! Vertical level index
        1439             :     INTEGER,        INTENT(IN) :: tIdx     ! Time index
        1440             :     TYPE(DataCont), POINTER    :: Dct      ! Data Container object
        1441             : !
        1442             : ! !RETURN VALUE:
        1443             : !
        1444             :     REAL(sp)                   :: dataVal  ! Data at this grid box and time
        1445             : !
        1446             : ! !REMARKS:
        1447             : !  This code was abstracted out of Get_Current_Emisssions for clarity.
        1448             : !  We have refactored the code to remove ELSE blocks for better efficiency.
        1449             : !EOP
        1450             : !------------------------------------------------------------------------------
        1451             : !BOC
        1452             :
        1453             :     ! Data is a 1-D scaler: Return a uniform value
        1454           0 :     IF ( Dct%Dta%SpaceDim == 1 ) THEN
        1455           0 :        dataVal = Dct%Dta%V2(tIDx)%Val(1,1)
        1456           0 :        RETURN
        1457             :     ENDIF
        1458             :
        1459             :     ! Data is a 2-D array: Return value at (I,J)
        1460           0 :     IF ( Dct%Dta%SpaceDim == 2 ) THEN
        1461           0 :        dataVal = Dct%Dta%V2(tIDx)%Val(I,J)
        1462           0 :        RETURN
        1463             :     ENDIF
        1464             :
        1465             :     ! Data is a 3-D array: Return value at (I,J,L)
        1466           0 :     dataVal = Dct%Dta%V3(tIDx)%Val(I,J,L)
        1467             :
        1468           0 :   END FUNCTION Get_Value_From_DataCont
        1469             : !EOC
        1470             : !------------------------------------------------------------------------------
        1471             : !                  GEOS-Chem Global Chemical Transport Model                  !
        1472             : !------------------------------------------------------------------------------
        1473             : !BOP
        1474             : !
        1475             : ! !IROUTINE: Apply_Dilution_Factor
        1476             : !
        1477             : ! !DESCRIPTION: Applies the dilution factor to input data.  This algorithm
        1478             : !  has been abstracted out of Get_Current_Emissions for computational
        1479             : !  efficiency.
        1480             : !\\
        1481             : !\\
        1482             : ! !INTERFACE:
        1483             : !
        1484           0 :   SUBROUTINE Apply_Dilution_Factor( I,       J,         L,      LowLL,       &
        1485             :                                     UppLL,   HcoState, BaseDct, inData,      &
        1486             :                                     outData, RC                             )
        1487             : !
        1488             : ! !USES:
        1489             : !
        1490             :     USE HCO_Error_Mod
        1491             :     USE HCO_State_Mod, ONLY : HCO_State
        1492             : !
        1493             : ! !INPUT PARAMETERS:
        1494             : !
        1495             :     INTEGER,         INTENT(IN)    :: I         ! Longitude (or X-dim) index
        1496             :     INTEGER,         INTENT(IN)    :: J         ! Latitude  (or Y-dim) index
        1497             :     INTEGER,         INTENT(IN)    :: L         ! Vertical level index
        1498             :     INTEGER,         INTENT(IN)    :: LowLL     ! Lower level for emissions
        1499             :     INTEGER,         INTENT(IN)    :: UppLL     ! Upper level for emissions
        1500             :     TYPE(DataCont),  POINTER       :: BaseDct   ! Base data container
        1501             :     REAL(sp),        INTENT(IN)    :: inData    ! Input data
        1502             : !
        1503             : ! !INPUT/OUTPUT PARAMETERS:
        1504             : !
        1505             :     TYPE(HCO_State), POINTER       :: HcoState  ! HEMCO state object
        1506             :     REAL(hp),        INTENT(INOUT) :: outData   ! Data w/ dil. factor applied
        1507             : !
        1508             : ! !OUTPUT PARAMETERS:
        1509             : !
        1510             :     INTEGER,         INTENT(OUT)   :: RC        ! Success or failure
        1511             : !
        1512             : !
        1513             : ! !REMARKS:
        1514             : !  This code was abstracted out of Get_Current_Emisssions for clarity.  We
        1515             : !  have also refactored the code to remove ELSE blocks for better efficiency.
        1516             : !EOP
        1517             : !------------------------------------------------------------------------------
        1518             : !BOC
        1519             : !
        1520             : ! !LOCAL VARIABLES:
        1521             : !
        1522             :     ! Scalars
        1523             :     REAL(hp)           :: dilFact
        1524             :
        1525             :     ! Strings
        1526             :     CHARACTER(LEN=255) :: errMsg, thisLoc
        1527             :
        1528             :     !========================================================================
        1529             :     ! Apply_Dilution_Factor begins here!
        1530             :     !========================================================================
        1531             :
        1532             :     ! Initialize
        1533           0 :     RC      = HCO_SUCCESS
        1534           0 :     errMsg  = ''
        1535           0 :     thisLoc = 'Apply_Dilution_Factor (in src/Core/hco_calc_mod.F90)'
        1536           0 :     dilFact = 0.0_hp
        1537             :
        1538             :     !========================================================================
        1539             :     ! Get dilution factor. Never dilute 3D emissions.
        1540             :     !========================================================================
        1541           0 :     IF ( BaseDct%Dta%SpaceDim == 3 ) THEN
        1542           0 :        outData = inData
        1543           0 :        RETURN
        1544             :     ENDIF
        1545             :
        1546             :     !========================================================================
        1547             :     ! If emission level mode is 2, copy emissions to all level
        1548             :     ! A separate scale factor should be used to distribute vertically
        1549             :     !========================================================================
        1550           0 :     IF ( BaseDct%Dta%EmisLmode == 2 ) THEN
        1551           0 :        outData = inData
        1552           0 :        RETURN
        1553             :     ENDIF
        1554             :
        1555             :     !========================================================================
        1556             :     ! Otherwise, compute the vertical dilution factor
        1557             :     ! and apply it to the input data.
        1558             :     !========================================================================
        1559             :     CALL GetDilFact( I          = I,                                         &
        1560             :                      J          = J,                                         &
        1561             :                      L          = L,                                         &
        1562             :                      LowLL      = LowLL,                                     &
        1563             :                      UppLL      = UppLL,                                     &
        1564             :                      HcoState   = HcoState,                                  &
        1565             :                      EmisL1     = BaseDct%Dta%EmisL1,                        &
        1566             :                      EmisL1Unit = BaseDct%Dta%EmisL1Unit,                    &
        1567             :                      EmisL2     = BaseDct%Dta%EmisL2,                        &
        1568             :                      EmisL2Unit = BaseDct%Dta%EmisL2Unit,                    &
        1569             :                      DilFact    = dilFact,                                   &
        1570           0 :                      RC         = RC                                        )
        1571             :
        1572             :     ! Trap error
        1573           0 :     IF ( RC /= HCO_SUCCESS ) THEN
        1574           0 :        WRITE(*,*) 'Error getting dilution factor at ',I,J,&
        1575           0 :             ': ', TRIM(BaseDct%cName)
        1576           0 :        RC = 1
        1577           0 :        RETURN
        1578             :     ENDIF
        1579             :
        1580             :     ! Scale base emission by dilution factor
        1581           0 :     outData = dilFact * inData
        1582             :
        1583             :   END SUBROUTINE Apply_Dilution_Factor
        1584             : !EOC
        1585             : !------------------------------------------------------------------------------
        1586             : !                  GEOS-Chem Global Chemical Transport Model                  !
        1587             : !------------------------------------------------------------------------------
        1588             : !BOP
        1589             : !
        1590             : ! !IROUTINE: Apply_Scale_Factor
        1591             : !
        1592             : ! !DESCRIPTION: Applies scale factors to the input
        1593             : !\\
        1594             : !\\
        1595             : ! !INTERFACE:
        1596             : !
        1597           0 :   SUBROUTINE Apply_Scale_Factor( I, J, L, ScalDct, scalFac, dataVal, RC )
        1598             : !
        1599             : ! !USES:
        1600             : !
        1601             :     USE HCO_Error_Mod
        1602             : !
        1603             : ! !INPUT PARAMETERS:
        1604             : !
        1605             :     INTEGER,        INTENT(IN)    :: I        ! Longitude (or X-dim) index
        1606             :     INTEGER,        INTENT(IN)    :: J        ! Latitude  (or Y-dim) index
        1607             :     INTEGER,        INTENT(IN)    :: L        ! Vertical level index
        1608             :     TYPE(DataCont), POINTER       :: ScalDct  ! Scale Factor data container
        1609             :     REAL(sp),       INTENT(IN)    :: scalFac  ! Scale factor
        1610             : !
        1611             : ! !INPUT/OUTPUT PARAMETERS:
        1612             : !
        1613             :     REAL(hp),       INTENT(INOUT) :: dataVal  ! Data to be scaled
        1614             : !
        1615             : ! !OUTPUT PARAMETERS:
        1616             : !
        1617             :     INTEGER,        INTENT(OUT)   :: RC       ! Success or failure
        1618             : !
        1619             : ! !REMARKS:
        1620             : !  This code was abstracted out of Get_Current_Emisssions for clarity.  We
        1621             : !  have also refactored the code to remove ELSE blocks for better efficiency.
        1622             : !EOP
        1623             : !------------------------------------------------------------------------------
        1624             : !BOC
        1625             : !
        1626             : ! !LOCAL VARIABLES:
        1627             : !
        1628             :     ! Strings
        1629             :     CHARACTER(LEN=255) :: errMsg, thisLoc
        1630             :
        1631             :     !========================================================================
        1632             :     ! Apply_Scale_Factor begins here!
        1633             :     !========================================================================
        1634             :
        1635             :     ! Initialize
        1636           0 :     RC      = HCO_SUCCESS
        1637           0 :     errMsg  = ''
        1638           0 :     thisLoc = 'Apply_Scale_Factor (src/Core/hco_calc_mod.F90)'
        1639             :
        1640             :     !------------------------------------------------------------------------
        1641             :     ! Operation = 1: multiply
        1642             :     !------------------------------------------------------------------------
        1643           0 :     IF ( ScalDct%Oper == 1 ) THEN
        1644           0 :        dataVal = dataVal * scalFac
        1645           0 :        RETURN
        1646             :     ENDIF
        1647             :
        1648             :     !------------------------------------------------------------------------
        1649             :     ! Operation = -1: divide
        1650             :     !------------------------------------------------------------------------
        1651           0 :     IF ( ScalDct%Oper == -1 ) THEN
        1652           0 :        IF ( scalFac /= 0.0_sp ) THEN
        1653           0 :            dataVal = dataVal / scalFac
        1654             :        ENDIF
        1655           0 :        RETURN
        1656             :     ENDIF
        1657             :
        1658             :     !------------------------------------------------------------------------
        1659             :     ! Operation = 2: square
        1660             :     !------------------------------------------------------------------------
        1661           0 :     IF ( ScalDct%Oper == 2 ) THEN
        1662           0 :        dataVal = dataVal * scalFac * scalFac
        1663           0 :        RETURN
        1664             :     ENDIF
        1665             :
        1666             :     !------------------------------------------------------------------------
        1667             :     ! Return w/ error otherwise (Oper 3 is only allowed for masks!)
        1668             :     !------------------------------------------------------------------------
        1669           0 :     WRITE( errMsg, * ) 'Illegal operator for: ', TRIM( ScalDct%cName ),      &
        1670           0 :                        ' operation: ', ScalDct%Oper
        1671           0 :     CALL HCO_Error( ErrMsg, RC, thisLoc )
        1672           0 :     RC = 2
        1673             :
        1674             :   END SUBROUTINE Apply_Scale_Factor
        1675             : !EOC
        1676             : !------------------------------------------------------------------------------
        1677             : !                   Harmonized Emissions Component (HEMCO)                    !
        1678             : !------------------------------------------------------------------------------
        1679             : !BOP
        1680             : !
        1681             : ! !IROUTINE: HCO_EvalFld_3D
        1682             : !
        1683             : ! !DESCRIPTION: Subroutine HCO\_EvalFld\_3D returns the 3D data field belonging
        1684             : !  to the emissions list data container with field name 'cName'. The returned
        1685             : !  data field is the completely evaluated field, e.g. the base field multiplied
        1686             : !  by all scale factors and with all masking being applied (as specified in the
        1687             : !  HEMCO configuration file). This distinguished this routine from HCO\_GetPtr
        1688             : !  in hco\_emislist\_mod.F90, which returns a reference to the unevaluated data
        1689             : !  field.
        1690             : !\\
        1691             : !\\
        1692             : ! !INTERFACE:
        1693             : !
        1694           0 :   SUBROUTINE HCO_EvalFld_3D( HcoState, cName, Arr3D, RC, FOUND )
        1695             : !
        1696             : ! !USES:
        1697             : !
        1698             :     USE HCO_STATE_MOD,    ONLY : HCO_State
        1699             :     USE HCO_DATACONT_MOD, ONLY : ListCont_Find
        1700             : !
        1701             : ! !INPUT PARAMETERS:
        1702             : !
        1703             :     CHARACTER(LEN=*), INTENT(IN   )  :: cName
        1704             : !
        1705             : ! !INPUT/OUTPUT PARAMETERS:
        1706             : !
        1707             :     TYPE(HCO_State),  POINTER        :: HcoState     ! HEMCO state object
        1708             :     REAL(hp),         INTENT(INOUT)  :: Arr3D(:,:,:) ! 3D array
        1709             :     INTEGER,          INTENT(INOUT)  :: RC           ! Return code
        1710             : !
        1711             : ! !OUTPUT PARAMETERS:
        1712             : !
        1713             :     LOGICAL,          INTENT(  OUT), OPTIONAL  :: FOUND
        1714             : !
        1715             : ! !REVISION HISTORY:
        1716             : !  11 May 2015 - C. Keller   - Initial Version
        1717             : !  See https://github.com/geoschem/hemco for complete history
        1718             : !EOP
        1719             : !------------------------------------------------------------------------------
        1720             : !BOC
        1721             : !
        1722             : ! !LOCAL VARIABLES:
        1723             : !
        1724             :     ! Scalars
        1725             :     LOGICAL                 :: FND
        1726             :     INTEGER                 :: AS, nI, nJ, nL, FLAG
        1727             :
        1728             :     ! Arrays
        1729           0 :     REAL(hp), ALLOCATABLE   :: Mask(:,:,:)
        1730             :
        1731             :     ! Working pointers: list and data container
        1732             :     TYPE(ListCont), POINTER :: Lct
        1733             :
        1734             :     ! For error handling & verbose mode
        1735             :     CHARACTER(LEN=255)      :: MSG
        1736             :     CHARACTER(LEN=255)      :: LOC = "HCO_EvalFld_3d (HCO_calc_mod.F90)"
        1737             :
        1738             :     !=================================================================
        1739             :     ! HCO_EvalFld_3D begins here!
        1740             :     !=================================================================
        1741             :
        1742             :     ! Init
        1743           0 :     RC    = HCO_SUCCESS
        1744           0 :     Lct   => NULL()
        1745           0 :     IF ( PRESENT(FOUND) ) FOUND = .FALSE.
        1746             :
        1747             :     ! Search for base container
        1748           0 :     CALL ListCont_Find ( HcoState%EmisList, TRIM(cName), FND, Lct )
        1749           0 :     IF ( PRESENT(FOUND) ) FOUND = FND
        1750             :
        1751             :     ! If not found, return here
        1752           0 :     IF ( .NOT. FND ) THEN
        1753           0 :        IF ( PRESENT(FOUND) ) THEN
        1754           0 :           RETURN
        1755             :        ELSE
        1756           0 :           MSG = 'Cannot find in EmisList: ' // TRIM(cName)
        1757           0 :           CALL HCO_ERROR( MSG, RC, THISLOC=LOC )
        1758           0 :           RETURN
        1759             :        ENDIF
        1760             :     ENDIF
        1761             :
        1762             :     ! Init
        1763           0 :     Arr3D = 0.0_hp
        1764             :
        1765             :     ! Define output dimensions
        1766           0 :     nI = SIZE(Arr3D,1)
        1767           0 :     nJ = SIZE(Arr3D,2)
        1768           0 :     nL = SIZE(Arr3D,3)
        1769             :
        1770             :     ! Sanity check: horizontal grid dimensions are expected to be on HEMCO grid
        1771           0 :     IF ( nI /= HcoState%NX .OR. nJ /= HcoState%nY ) THEN
        1772           0 :        WRITE(MSG,*) "Horizontal dimension error: ", TRIM(cName), nI, nJ
        1773           0 :        CALL HCO_ERROR( MSG, RC, THISLOC=LOC )
        1774           0 :        RETURN
        1775             :     ENDIF
        1776             :
        1777             :     ! Make sure mask array is defined
        1778           0 :     ALLOCATE(MASK(nI,nJ,nL),STAT=AS)
        1779           0 :     IF ( AS /= 0 ) THEN
        1780           0 :        CALL HCO_ERROR( 'Cannot allocate MASK', RC, THISLOC=LOC )
        1781           0 :        RETURN
        1782             :     ENDIF
        1783           0 :     mask = 0.0_hp
        1784             :
        1785             :     ! Calculate emissions for base container
        1786             :     CALL GET_CURRENT_EMISSIONS( HcoState, Lct%Dct, nI,   nJ,                 &
        1787           0 :                                 nL,       Arr3D,   Mask, RC                 )
        1788           0 :     IF ( RC /= HCO_SUCCESS ) THEN
        1789           0 :         CALL HCO_ERROR( 'ERROR 14', RC, THISLOC=LOC )
        1790           0 :         RETURN
        1791             :     ENDIF
        1792             :
        1793             :     ! All done
        1794           0 :     IF (ALLOCATED(MASK) ) DEALLOCATE(MASK)
        1795           0 :     Lct => NULL()
        1796             :
        1797           0 :   END SUBROUTINE HCO_EvalFld_3D
        1798             : !EOC
        1799             : !------------------------------------------------------------------------------
        1800             : !                   Harmonized Emissions Component (HEMCO)                    !
        1801             : !------------------------------------------------------------------------------
        1802             : !BOP
        1803             : !
        1804             : ! !IROUTINE: HCO_EvalFld_2D
        1805             : !
        1806             : ! !DESCRIPTION: Subroutine HCO\_EvalFld\_2D returns the 2D data field belonging
        1807             : !  to the emissions list data container with field name 'cName'. The returned
        1808             : !  data field is the completely evaluated field, e.g. the base field multiplied
        1809             : !  by all scale factors and with all masking being applied (as specified in the
        1810             : !  HEMCO configuration file). This distinguished this routine from HCO\_GetPtr
        1811             : !  in hco\_emislist\_mod.F90, which returns a reference to the unevaluated data
        1812             : !  field.
        1813             : !\\
        1814             : !\\
        1815             : !\\
        1816             : ! !INTERFACE:
        1817             : !
        1818           0 :   SUBROUTINE HCO_EvalFld_2D( HcoState, cName, Arr2D, RC, FOUND )
        1819             : !
        1820             : ! !USES:
        1821             : !
        1822             :     USE HCO_STATE_MOD,    ONLY : HCO_State
        1823             :     USE HCO_DATACONT_MOD, ONLY : ListCont_Find
        1824             : !
        1825             : ! !INPUT PARAMETERS:
        1826             : !
        1827             :     CHARACTER(LEN=*), INTENT(IN   )  :: cName
        1828             : !
        1829             : ! !INPUT/OUTPUT PARAMETERS:
        1830             : !
        1831             :     TYPE(HCO_State),  POINTER        :: HcoState     ! HEMCO state object
        1832             :     REAL(hp),         INTENT(INOUT)  :: Arr2D(:,:)   ! 2D array
        1833             :     INTEGER,          INTENT(INOUT)  :: RC           ! Return code
        1834             : !
        1835             : ! !OUTPUT PARAMETERS:
        1836             : !
        1837             :     LOGICAL,          INTENT(  OUT), OPTIONAL  :: FOUND
        1838             : !
        1839             : ! !REVISION HISTORY:
        1840             : !  11 May 2015 - C. Keller   - Initial Version
        1841             : !  See https://github.com/geoschem/hemco for complete history
        1842             : !EOP
        1843             : !------------------------------------------------------------------------------
        1844             : !BOC
        1845             : !
        1846             : ! !LOCAL VARIABLES:
        1847             : !
        1848             :     ! Scalars
        1849             :     LOGICAL                 :: FND
        1850             :     INTEGER                 :: AS, nI, nJ, nL, UseLL, FLAG
        1851             :
        1852             :     ! Arrays
        1853           0 :     REAL(hp), ALLOCATABLE   :: Mask (:,:,:)
        1854           0 :     REAL(hp), ALLOCATABLE   :: Arr3D(:,:,:)
        1855             :
        1856             :     ! Working pointers: list and data container
        1857             :     TYPE(ListCont), POINTER :: Lct
        1858             :
        1859             :     ! For error handling & verbose mode
        1860             :     CHARACTER(LEN=255)      :: MSG
        1861             :     CHARACTER(LEN=255)      :: LOC = "HCO_EvalFld_2d (HCO_calc_mod.F90)"
        1862             :
        1863             :     !=================================================================
        1864             :     ! HCO_EvalFld_2D begins here!
        1865             :     !=================================================================
        1866             :
        1867             :     ! Init
        1868           0 :     RC    = HCO_SUCCESS
        1869           0 :     Lct   => NULL()
        1870           0 :     IF ( PRESENT(FOUND) ) FOUND = .FALSE.
        1871             :
        1872             :     ! Search for base container
        1873           0 :     CALL ListCont_Find ( HcoState%EmisList, TRIM(cName), FND, Lct )
        1874           0 :     IF ( PRESENT(FOUND) ) FOUND = FND
        1875             :
        1876             :     ! If not found, return here
        1877           0 :     IF ( .NOT. FND ) THEN
        1878           0 :        IF ( PRESENT(FOUND) ) THEN
        1879           0 :           RETURN
        1880             :        ELSE
        1881           0 :           MSG = 'Cannot find in EmisList: ' // TRIM(cName)
        1882           0 :           CALL HCO_ERROR( MSG, RC, THISLOC=LOC )
        1883           0 :           RETURN
        1884             :        ENDIF
        1885             :     ENDIF
        1886             :
        1887             :     ! Init Arr2D
        1888           0 :     Arr2D = 0.0_hp
        1889             :
        1890             :     ! Define output dimensions
        1891           0 :     nI = SIZE(Arr2D,1)
        1892           0 :     nJ = SIZE(Arr2D,2)
        1893           0 :     nL = 1
        1894             :
        1895             :     ! Sanity check: horizontal grid dimensions are expected to be on HEMCO grid
        1896           0 :     IF ( nI /= HcoState%NX .OR. nJ /= HcoState%nY ) THEN
        1897           0 :        WRITE(MSG,*) "Horizontal dimension error: ", TRIM(cName), nI, nJ
        1898           0 :        CALL HCO_ERROR( MSG, RC, THISLOC=LOC )
        1899           0 :        RETURN
        1900             :     ENDIF
        1901             :
        1902             :     ! Make sure mask array is defined
        1903           0 :     ALLOCATE(MASK(nI,nJ,nL),Arr3D(nI,nJ,nL),STAT=AS)
        1904           0 :     IF ( AS /= 0 ) THEN
        1905           0 :        CALL HCO_ERROR( 'Cannot allocate MASK', RC, THISLOC=LOC )
        1906           0 :        RETURN
        1907             :     ENDIF
        1908           0 :     Arr3D = 0.0_hp
        1909           0 :     Mask  = 0.0_hp
        1910             :
        1911             :     ! Calculate emissions for base container
        1912             :     CALL GET_CURRENT_EMISSIONS( HcoState, Lct%Dct, nI,   nJ,                 &
        1913           0 :                                 nL,       Arr3D,   Mask, RC, UseLL=UseLL    )
        1914           0 :     IF ( RC /= HCO_SUCCESS ) THEN
        1915           0 :         CALL HCO_ERROR( 'ERROR 15', RC, THISLOC=LOC )
        1916           0 :         RETURN
        1917             :     ENDIF
        1918             :
        1919             :     ! Place 3D array into 2D array. UseLL returns the vertical level into which
        1920             :     ! emissions have been added within GET_CURRENT_EMISSIONS. This should be
        1921             :     ! level 1 for most cases but it can be another level if specified so.
        1922             :     ! Return a warning if level is not 1 (ckeller, 11/1/16).
        1923           0 :     UseLL = MIN( MAX(useLL,1), SIZE(Arr3D,3) )
        1924           0 :     IF ( UseLL /= 1 ) THEN
        1925           0 :        WRITE(MSG,*) "2D data was emitted above surface - this information might be lost: " , TRIM(cName), UseLL
        1926           0 :        CALL HCO_WARNING( HcoState%Config%Err, MSG, RC, THISLOC=LOC )
        1927             :     ENDIF
        1928             :
        1929             :     ! Pass 3D data to 2D array
        1930           0 :     Arr2D = Arr3D(:,:,UseLL)
        1931             :
        1932             :     ! All done
        1933           0 :     IF (ALLOCATED(MASK ) ) DEALLOCATE(MASK )
        1934           0 :     IF (ALLOCATED(Arr3D) ) DEALLOCATE(Arr3D)
        1935           0 :     Lct => NULL()
        1936             :
        1937           0 :   END SUBROUTINE HCO_EvalFld_2D
        1938             : !EOC
        1939             : !------------------------------------------------------------------------------
        1940             : !                   Harmonized Emissions Component (HEMCO)                    !
        1941             : !------------------------------------------------------------------------------
        1942             : !BOP
        1943             : !
        1944             : ! !IROUTINE: GetMaskVal
        1945             : !
        1946             : ! !DESCRIPTION: Subroutine GetMaskVal is a helper routine to get the mask
        1947             : !  value at a given location.
        1948             : !\\
        1949             : !\\
        1950             : ! !INTERFACE:
        1951             : !
        1952           0 :   SUBROUTINE GetMaskVal ( Dct, I, J, MaskVal, Fractions, RC )
        1953             : !
        1954             : ! !USES:
        1955             : !
        1956             : !
        1957             : ! !INPUT PARAMETERS:
        1958             : !
        1959             :     INTEGER,         INTENT(IN   ) :: I                   ! # of lons
        1960             :     INTEGER,         INTENT(IN   ) :: J                   ! # of lats
        1961             :     LOGICAL,         INTENT(IN   ) :: Fractions           ! Use fractions?
        1962             : !
        1963             : ! !INPUT/OUTPUT PARAMETERS:
        1964             : !
        1965             :     TYPE(DataCont),  POINTER       :: Dct                 ! Mask container
        1966             :     REAL(sp),        INTENT(INOUT) :: MaskVal
        1967             :     INTEGER,         INTENT(INOUT) :: RC
        1968             : !
        1969             : ! !REVISION HISTORY:
        1970             : !  09 Apr 2015 - C. Keller   - Initial Version
        1971             : !  See https://github.com/geoschem/hemco for complete history
        1972             : !EOP
        1973             : !------------------------------------------------------------------------------
        1974             : !BOC
        1975             : !
        1976             : ! !LOCAL VARIABLES:
        1977             : !
        1978             :
        1979             :     !=================================================================
        1980             :     ! GetMaskVal begins here
        1981             :     !=================================================================
        1982             :
        1983             :     ! Mask value over this grid box
        1984           0 :     MaskVal = Dct%Dta%V2(1)%Val(I,J)
        1985             :
        1986             :     ! Negative mask values are treated as zero (exclude).
        1987           0 :     IF ( (MaskVal <= 0.0_sp) .OR. (MaskVal == HCO_MISSVAL) ) THEN
        1988           0 :        MaskVal = 0.0_sp
        1989           0 :     ELSEIF ( MaskVal > 1.0_sp ) THEN
        1990           0 :        MaskVal = 1.0_sp
        1991             :     ENDIF
        1992             :
        1993             :     ! For operator set to 3, mirror value
        1994             :     ! MaskVal=1 becomes 0 and MaskVal=0/missing becomes 1
        1995           0 :     IF ( Dct%Oper == 3 ) THEN
        1996           0 :        IF ( (MaskVal == 0.0_sp) .OR. (MaskVal == HCO_MISSVAL) ) THEN
        1997           0 :           MaskVal = 1.0_sp
        1998           0 :        ELSEIF ( MaskVal == 1.0_sp ) THEN
        1999           0 :           MaskVal = 1.0_sp - MaskVal
        2000             :        ENDIF
        2001             :     ENDIF
        2002             :
        2003             :     ! Treat as binary?
        2004           0 :     IF ( .NOT. Fractions ) THEN
        2005           0 :        IF ( MaskVal < MASK_THRESHOLD ) THEN
        2006           0 :           MaskVal = 0.0_sp
        2007             :        ELSE
        2008           0 :           MaskVal = 1.0_sp
        2009             :        ENDIF
        2010             :     ENDIF
        2011             :
        2012             :     ! Return w/ success
        2013           0 :     RC = HCO_SUCCESS
        2014             :
        2015           0 :   END SUBROUTINE GetMaskVal
        2016             : !EOC
        2017             : !------------------------------------------------------------------------------
        2018             : !                   Harmonized Emissions Component (HEMCO)                    !
        2019             : !------------------------------------------------------------------------------
        2020             : !BOP
        2021             : !
        2022             : ! !IROUTINE: HCO_MaskFld
        2023             : !
        2024             : ! !DESCRIPTION: Subroutine HCO\_MaskFld is a helper routine to get the mask
        2025             : ! field with the given name. The returned mask field is fully evaluated,
        2026             : ! e.g. the data operation flag associated with this mask field is already
        2027             : ! taken into account. For instance, if the data operator of a mask field is
        2028             : ! set to 3, the returned array contains already the mirrored mask values.
        2029             : !\\
        2030             : !\\
        2031             : ! !INTERFACE:
        2032             : !
        2033           0 :   SUBROUTINE HCO_MaskFld ( HcoState, MaskName, Mask, RC, FOUND )
        2034             : !
        2035             : ! !USES:
        2036             : !
        2037             :     USE HCO_STATE_MOD,    ONLY : HCO_State
        2038             :     USE HCO_DATACONT_MOD, ONLY : ListCont_Find
        2039             : !
        2040             : ! !INPUT PARAMETERS:
        2041             : !
        2042             :     TYPE(HCO_STATE), POINTER                 :: HcoState
        2043             :     CHARACTER(LEN=*),INTENT(IN   )           :: MaskName
        2044             : !
        2045             : ! !INPUT/OUTPUT PARAMETERS:
        2046             : !
        2047             :     REAL(sp),        INTENT(INOUT)           :: Mask(:,:)
        2048             :     INTEGER,         INTENT(INOUT)           :: RC
        2049             : !
        2050             : ! !OUTPUT PARAMETERS:
        2051             : !
        2052             :     LOGICAL,         INTENT(  OUT), OPTIONAL :: FOUND
        2053             : !
        2054             : ! !REVISION HISTORY:
        2055             : !  11 Jun 2015 - C. Keller   - Initial Version
        2056             : !  See https://github.com/geoschem/hemco for complete history
        2057             : !EOP
        2058             : !------------------------------------------------------------------------------
        2059             : !BOC
        2060             : !
        2061             : ! !LOCAL VARIABLES:
        2062             : !
        2063             :     INTEGER                 :: I, J, FLAG
        2064             :
        2065             :     LOGICAL                 :: FND, ERR
        2066             :     LOGICAL                 :: Fractions
        2067             :
        2068             :     TYPE(ListCont), POINTER :: MaskLct
        2069             :
        2070             :     CHARACTER(LEN=255)      :: MSG
        2071             :     CHARACTER(LEN=255)      :: LOC = 'HCO_MaskFld (hco_calc_mod.F90)'
        2072             :
        2073             :     !=================================================================
        2074             :     ! HCO_MaskFld begins here
        2075             :     !=================================================================
        2076             :
        2077             :     ! Nullify
        2078           0 :     MaskLct  => NULL()
        2079             :
        2080             :     ! Init: default is mask value of 1
        2081           0 :     MASK = 1.0_sp
        2082           0 :     ERR  = .FALSE.
        2083           0 :     FND  = .FALSE.
        2084             :
        2085             :     ! Search for mask field within EmisList
        2086           0 :     CALL ListCont_Find ( HcoState%EmisList, TRIM(MaskName), FND, MaskLct )
        2087             :
        2088           0 :     IF ( .NOT. FND .AND. .NOT. PRESENT(FOUND) ) THEN
        2089           0 :        MSG = 'Cannot find mask field ' // TRIM(MaskName)
        2090           0 :        CALL HCO_MSG(HcoState%Config%Err,MSG,SEP1='!')
        2091             :        MSG = 'Make sure this field is listed in the mask section '  // &
        2092             :            'of the HEMCO configuration file. You may also need to ' // &
        2093           0 :            'set the optional attribute `ReadAlways` to `yes`, e.g.'
        2094           0 :        CALL HCO_MSG(HcoState%Config%Err,MSG)
        2095           0 :        MSG = '5000 TESTMASK     -140/10/-40/90 - - - xy 1 1 -140/10/-40/90 yes'
        2096           0 :        CALL HCO_MSG(HcoState%Config%Err,MSG)
        2097             :        CALL HCO_ERROR ( &
        2098           0 :                         'Error reading mask '//TRIM(MaskName), RC, THISLOC=LOC )
        2099           0 :        RETURN
        2100             :     ENDIF
        2101           0 :     IF ( PRESENT(FOUND) ) FOUND = FND
        2102             :
        2103             :     ! Do only if found
        2104           0 :     IF ( FND ) THEN
        2105             :
        2106             :        ! Use mask fractions?
        2107           0 :        Fractions = HcoState%Options%MaskFractions
        2108             :
        2109             :        ! Make sure mask array has correct dimensions
        2110           0 :        IF ( SIZE(MASK,1) /= HcoState%NX .OR. SIZE(MASK,2) /= HcoState%NY ) THEN
        2111           0 :           WRITE(MSG,*) 'Input mask array has wrong dimensions. Must be ', &
        2112           0 :              HcoState%NX, HcoState%NY, ' but found ', SIZE(MASK,1), SIZE(MASK,2)
        2113           0 :           CALL HCO_ERROR ( MSG, RC, THISLOC=LOC )
        2114           0 :           RETURN
        2115             :        ENDIF
        2116             :
        2117             :        ! Do for every grid box
        2118             :        !$OMP PARALLEL DO            &
        2119             :        !$OMP DEFAULT( SHARED      ) &
        2120             :        !$OMP PRIVATE( I, J        )
        2121           0 :        DO J = 1, HcoState%NY
        2122           0 :        DO I = 1, HcoState%NX
        2123           0 :           CALL GetMaskVal( MaskLct%Dct, I, J, Mask(I,J), Fractions, RC )
        2124           0 :           IF ( RC /= HCO_SUCCESS ) THEN
        2125             :              ERR = .TRUE.
        2126             :              EXIT
        2127             :           ENDIF
        2128             :        ENDDO
        2129             :        ENDDO
        2130             :        !$OMP END PARALLEL DO
        2131             :
        2132             :        ! Error check
        2133           0 :        IF ( ERR ) THEN
        2134           0 :           MSG = 'Error in GetMaskVal'
        2135           0 :           CALL HCO_ERROR ( MSG, RC, THISLOC=LOC )
        2136           0 :           RETURN
        2137             :        ENDIF
        2138             :
        2139             :     ENDIF
        2140             :
        2141             :     ! Free pointer
        2142           0 :     MaskLct  => NULL()
        2143             :
        2144             :     ! Return w/ success
        2145           0 :     RC = HCO_SUCCESS
        2146             :
        2147             :   END SUBROUTINE HCO_MaskFld
        2148             : !EOC
        2149             : !------------------------------------------------------------------------------
        2150             : !                   Harmonized Emissions Component (HEMCO)                    !
        2151             : !------------------------------------------------------------------------------
        2152             : !BOP
        2153             : !
        2154             : ! !IROUTINE: GetVertIndx
        2155             : !
        2156             : ! !DESCRIPTION: Subroutine GetVertIndx is a helper routine to get the vertical
        2157             : !  index range of the given data field.
        2158             : !\\
        2159             : !\\
        2160             : ! !INTERFACE:
        2161             : !
        2162           0 :   SUBROUTINE GetVertIndx( HcoState, Dct,           isLevDct1,                &
        2163             :                           LevDct1,  LevDct1_Unit,  isLevDct2,                &
        2164             :                           LevDct2,  LevDct2_Unit,  I,                        &
        2165             :                           J,        LowLL,         UppLL,     RC            )
        2166             : !
        2167             : ! !USES:
        2168             : !
        2169             :     USE HCO_State_Mod, ONLY : HCO_State
        2170             : !
        2171             : ! !INPUT PARAMETERS:
        2172             : !
        2173             :     TYPE(HCO_State), POINTER       :: HcoState      ! HEMCO state object
        2174             :     LOGICAL,         INTENT(IN)    :: isLevDct1     ! Is LevDct1 not null?
        2175             :     TYPE(DataCont),  POINTER       :: LevDct1       ! Level index 1 container
        2176             :     INTEGER,         INTENT(IN)    :: LevDct1_Unit  ! LevDct1 unit code
        2177             :     LOGICAL,         INTENT(IN)    :: isLevDct2     ! Is LevDct2 not null?
        2178             :     TYPE(DataCont),  POINTER       :: LevDct2       ! Level index 2 container
        2179             :     INTEGER,         INTENT(IN)    :: LevDct2_Unit  ! LevDct2 unit code
        2180             :     INTEGER,         INTENT(IN)    :: I             ! lon index
        2181             :     INTEGER,         INTENT(IN)    :: J             ! lat index
        2182             : !
        2183             : ! !INPUT/OUTPUT PARAMETERS:
        2184             : !
        2185             :     TYPE(DataCont),  POINTER       :: Dct           ! Mask container
        2186             :     INTEGER,         INTENT(INOUT) :: LowLL         ! lower level index
        2187             :     INTEGER,         INTENT(INOUT) :: UppLL         ! upper level index
        2188             :     INTEGER,         INTENT(INOUT) :: RC
        2189             : !
        2190             : ! !REVISION HISTORY:
        2191             : !  06 May 2016 - C. Keller   - Initial Version
        2192             : !  See https://github.com/geoschem/hemco for complete history
        2193             : !EOP
        2194             : !------------------------------------------------------------------------------
        2195             : !BOC
        2196             : !
        2197             : ! !LOCAL VARIABLES:
        2198             : !
        2199             :     INTEGER             :: EmisLUnit
        2200             :     REAL(hp)            :: EmisL
        2201             :     CHARACTER(LEN=255)  :: errMsg, thisLoc
        2202             :
        2203             :     !=======================================================================
        2204             :     ! GetVertIndx begins here
        2205             :     !=======================================================================
        2206             :
        2207             :     ! Initialize
        2208           0 :     RC      = HCO_SUCCESS
        2209           0 :     errMsg  = ''
        2210           0 :     thisLoc = 'GetVertIndx (src/Core/hco_calc_mod.F90)'
        2211             :
        2212             :     !-----------------------------------------------------------------------
        2213             :     ! Get vertical extension of base emission array.
        2214             :     !
        2215             :     ! Unlike the output array OUTARR_3D, the data containers do not
        2216             :     ! necessarily extent over the entire troposphere but only cover
        2217             :     ! the effectively filled vertical levels. For most inventories,
        2218             :     ! this is only the first model level.
        2219             :     !-----------------------------------------------------------------------
        2220           0 :     IF ( Dct%Dta%SpaceDim == 3 ) THEN
        2221           0 :        LowLL = 1
        2222           0 :        UppLL = SIZE( Dct%Dta%V3(1)%Val, 3 )
        2223           0 :        RC    = HCO_SUCCESS
        2224           0 :        RETURN
        2225             :     ENDIF
        2226             :
        2227             :     !-----------------------------------------------------------------------
        2228             :     ! For 2D field, check if it shall be spread out over multiple
        2229             :     ! levels. Possible to go from PBL to max. specified level.
        2230             :     !-----------------------------------------------------------------------
        2231             :
        2232             :     ! Lower level
        2233             :     ! --> Check if scale factor is used to determine lower and/or
        2234             :     !     upper level
        2235             :     !
        2236             :     ! NOTE: Get rid of ELSE block for efficiency (bmy, 09 Mar 2023)
        2237           0 :     EmisL     = Dct%Dta%EmisL1
        2238           0 :     EmisLUnit = Dct%Dta%EmisL1Unit
        2239             :
        2240           0 :     IF ( isLevDct1 ) THEN
        2241           0 :        EmisL = GetEmisL( HcoState, LevDct1, I, J )
        2242           0 :        IF ( EmisL < 0.0_hp ) THEN
        2243           0 :           RC = HCO_FAIL
        2244           0 :           RETURN
        2245             :        ENDIF
        2246           0 :        EmisLUnit = LevDct1_Unit
        2247             :     ENDIF
        2248             :
        2249           0 :     CALL GetIdx( HcoState, I, J, EmisL, EmisLUnit, LowLL, RC )
        2250           0 :     IF ( RC /= HCO_SUCCESS ) THEN
        2251           0 :        errMsg = 'Error encountered in rouitine "GetIdx" (for LevDct1)!'
        2252           0 :        CALL HCO_ERROR( errMsg, RC, thisLoc )
        2253           0 :        RETURN
        2254             :     ENDIF
        2255             :
        2256             :     ! Upper level
        2257             :     ! NOTE: Get rid of ELSE block for efficiency (bmy, 09 Mar 2023)
        2258           0 :     EmisL     = Dct%Dta%EmisL2
        2259           0 :     EmisLUnit = Dct%Dta%EmisL2Unit
        2260             :
        2261           0 :     IF ( isLevDct2 ) THEN
        2262           0 :        EmisL = GetEmisL( HcoState, LevDct2, I, J )
        2263           0 :        IF ( EmisL < 0.0_hp ) THEN
        2264           0 :           RC = HCO_FAIL
        2265           0 :           RETURN
        2266             :        ENDIF
        2267           0 :        EmisLUnit = LevDct2_Unit
        2268             :     ENDIF
        2269             :
        2270           0 :     CALL GetIdx( HcoState, I, J, EmisL, EmisLUnit, UppLL, RC )
        2271           0 :     IF ( RC /= HCO_SUCCESS ) THEN
        2272           0 :        errMsg = 'Error encountered in rouitine "GetIdx" (for LevDct2)!'
        2273           0 :        CALL HCO_ERROR( errMsg, RC, thisLoc )
        2274           0 :        RETURN
        2275             :     ENDIF
        2276             :
        2277             :     ! Upper level must not be lower than lower level
        2278           0 :     UppLL = MAX(LowLL, UppLL)
        2279             :
        2280             :     ! Return w/ success
        2281           0 :     RC = HCO_SUCCESS
        2282             :
        2283             :   END SUBROUTINE GetVertIndx
        2284             : !EOC
        2285             : !------------------------------------------------------------------------------
        2286             : !                   Harmonized Emissions Component (HEMCO)                    !
        2287             : !------------------------------------------------------------------------------
        2288             : !BOP
        2289             : !
        2290             : ! !FUNCTION: GetEmisL
        2291             : !
        2292             : ! !DESCRIPTION: Returns the emission level read from a scale factor.
        2293             : !\\
        2294             : !\\
        2295             : ! !INTERFACE:
        2296             : !
        2297           0 :   FUNCTION GetEmisL( HcoState, LevDct, I, J ) RESULT ( EmisL )
        2298             : !
        2299             : ! !USES:
        2300             : !
        2301             :     USE HCO_TYPES_MOD
        2302             :     USE HCO_STATE_MOD,    ONLY : HCO_State
        2303             :     USE HCO_tIdx_MOD,     ONLY : tIDx_GetIndx
        2304             : !
        2305             : ! !INPUT PARAMETERS:
        2306             : !
        2307             :     TYPE(HCO_State), POINTER        :: HcoState       ! HEMCO state object
        2308             :     TYPE(DataCont),  POINTER        :: LevDct         ! Level index 1 container
        2309             :     INTEGER,         INTENT(IN   )  :: I, J           ! horizontal index
        2310             : !
        2311             : ! !RETURN VALUE:
        2312             : !
        2313             :     REAL(hp)                        :: EmisL
        2314             : !
        2315             : ! !REVISION HISTORY:
        2316             : !  26 Jan 2018 - C. Keller - Initial version
        2317             : !  See https://github.com/geoschem/hemco for complete history
        2318             : !EOP
        2319             : !------------------------------------------------------------------------------
        2320             : !BOC
        2321             : !
        2322             : ! !LOCAL VARIABLES:
        2323             : !
        2324             :     INTEGER  :: levtidx
        2325             :
        2326             :     !=================================================================
        2327             :     ! GetEmisL begins here
        2328             :     !=================================================================
        2329           0 :     levtidx = tIDx_GetIndx( HcoState, LevDct%Dta, I, J )
        2330           0 :     IF ( levtidx <= 0 ) THEN
        2331             :        WRITE(*,*)' Cannot get time slice for field '//&
        2332           0 :        TRIM(LevDct%cName)//': GetEmisL (hco_calc_mod.F90)'
        2333           0 :        EmisL = -1.0
        2334           0 :        RETURN
        2335             :     ENDIF
        2336             :
        2337           0 :     IF ( LevDct%Dta%SpaceDim == 1 ) THEN
        2338           0 :        EmisL = LevDct%Dta%V2(levtidx)%Val(1,1)
        2339           0 :     ELSEIF ( LevDct%Dta%SpaceDim == 2 ) THEN
        2340           0 :        EmisL = LevDct%Dta%V2(levtidx)%Val(I,J)
        2341           0 :     ELSEIF ( LevDct%Dta%SpaceDim == 3 ) THEN
        2342           0 :        EmisL = LevDct%Dta%V3(levtidx)%Val(I,J,1)
        2343             :     ENDIF
        2344             :
        2345           0 :     IF ( EmisL == HCO_MISSVAL ) EmisL = 0.0_hp
        2346             :
        2347             : END FUNCTION GetEmisL
        2348             : !EOC
        2349             : !------------------------------------------------------------------------------
        2350             : !                   Harmonized Emissions Component (HEMCO)                    !
        2351             : !------------------------------------------------------------------------------
        2352             : !BOP
        2353             : !
        2354             : ! !FUNCTION: GetEmisLUnit
        2355             : !
        2356             : ! !DESCRIPTION: Returns the emission level unit read from a scale factor.
        2357             : !\\
        2358             : !\\
        2359             : ! !INTERFACE:
        2360             : !
        2361           0 :   FUNCTION GetEmisLUnit( HcoState, LevDct ) RESULT( EmisLUnit )
        2362             : !
        2363             : ! !USES:
        2364             : !
        2365             :     USE HCO_TYPES_MOD
        2366             :     USE HCO_STATE_MOD, ONLY : HCO_State
        2367             : !
        2368             : ! !INPUT PARAMETERS:
        2369             : !
        2370             :     TYPE(HCO_State), POINTER        :: HcoState       ! HEMCO state object
        2371             :     TYPE(DataCont),  POINTER        :: LevDct         ! Level index 1 container
        2372             : !
        2373             : ! !RETURN VALUE:
        2374             : !
        2375             :     INTEGER                         :: EmisLUnit
        2376             : !
        2377             : ! !REVISION HISTORY:
        2378             : !  26 Jan 2018 - C. Keller - Initial version
        2379             : !  See https://github.com/geoschem/hemco for complete history
        2380             : !EOP
        2381             : !------------------------------------------------------------------------------
        2382             : !BOC
        2383             : !
        2384             : ! !LOCAL VARIABLES:
        2385             : !
        2386             :     !=================================================================
        2387             :     ! GetEmisLUnit begins here
        2388             :     !=================================================================
        2389             :
        2390             :     ! For now, only meters are supported
        2391           0 :     EmisLUnit = HCO_EMISL_M
        2392             :
        2393             :     ! Dummy check that units on field are actually in meters
        2394           0 :     IF ( TRIM(LevDct%Dta%OrigUnit) /= 'm' .AND. &
        2395             :          TRIM(LevDct%Dta%OrigUnit) /= '1'        ) THEN
        2396             :        WRITE(*,*) TRIM(LevDct%cName)// &
        2397             :        ' must have units of `m`, instead found '//&
        2398           0 :        TRIM(LevDct%Dta%OrigUnit)//': GetEmisLUnit (hco_calc_mod.F90)'
        2399           0 :        EmisLUnit = -1
        2400             :     ENDIF
        2401             :
        2402           0 : END FUNCTION GetEmisLUnit
        2403             : !EOC
        2404             : !------------------------------------------------------------------------------
        2405             : !                   Harmonized Emissions Component (HEMCO)                    !
        2406             : !------------------------------------------------------------------------------
        2407             : !BOP
        2408             : !
        2409             : ! !IROUTINE: GetIdx
        2410             : !
        2411             : ! !DESCRIPTION: Subroutine GetIdx is a helper routine to return the vertical
        2412             : !  level index for a given altitude. The altitude can be provided in level
        2413             : !  coordinates, in units of meters or as the 'PBL mixing height'.
        2414             : !\\
        2415             : !\\
        2416             : ! !INTERFACE:
        2417             : !
        2418           0 :   SUBROUTINE GetIdx( HcoState, I, J, alt, altu, lidx, RC )
        2419             : !
        2420             : ! !USES:
        2421             : !
        2422             :     USE HCO_TYPES_MOD
        2423             :     USE HCO_STATE_MOD,   ONLY : HCO_STATE
        2424             : !
        2425             : ! !INPUT PARAMETERS:
        2426             : !
        2427             :     TYPE(HCO_State), POINTER        :: HcoState       ! HEMCO state object
        2428             :     INTEGER,         INTENT(IN   )  :: I, J           ! horizontal index
        2429             :     INTEGER,         INTENT(IN   )  :: altu           ! altitude unit
        2430             : !
        2431             : ! !OUTPUT PARAMETERS:
        2432             : !
        2433             :     INTEGER,         INTENT(  OUT)  :: lidx           ! level index
        2434             : !
        2435             : ! !INPUT/OUTPUT PARAMETERS:
        2436             : !
        2437             :     REAL(hp),        INTENT(INOUT)  :: alt            ! altitude
        2438             :     INTEGER,         INTENT(INOUT)  :: RC
        2439             : !
        2440             : ! !REMARKS:
        2441             : !  The code was refactored to avoid ELSE blocks (which are computationally
        2442             : !  expensive) following the "never-nesting" technique. (Bob Y., 10 Mar 2023)
        2443             : !
        2444             : ! !REVISION HISTORY:
        2445             : !  09 May 2016 - C. Keller - Initial version
        2446             : !  See https://github.com/geoschem/hemco for complete history
        2447             : !EOP
        2448             : !------------------------------------------------------------------------------
        2449             : !BOC
        2450             : !
        2451             : ! !LOCAL VARIABLES:
        2452             : !
        2453             :     INTEGER                 :: L
        2454             :     REAL(hp)                :: altb, altt
        2455             :     CHARACTER(LEN=255)      :: MSG
        2456             :     CHARACTER(LEN=255)      :: LOC = 'GetIdx (hco_calc_mod.F90)'
        2457             :
        2458             :     !=================================================================
        2459             :     ! HCO_GetIdx begins here
        2460             :     !=================================================================
        2461             :
        2462             :     ! Init
        2463           0 :     RC = HCO_SUCCESS
        2464             :
        2465             :     ! Input data is in level coordinates;
        2466             :     ! Return the corresponding level index
        2467           0 :     IF ( altu == HCO_EMISL_LEV ) THEN
        2468           0 :        lidx = INT( alt )
        2469           0 :        RETURN
        2470             :     ENDIF
        2471             :
        2472             :     ! Input specifies the top-of-atmosphere;
        2473             :     ! Return the top-of-atmosphere level index
        2474           0 :     IF ( altu == HCO_EMISL_TOP ) THEN
        2475           0 :        lidx = HCOState%NZ
        2476           0 :        RETURN
        2477             :     ENDIF
        2478             :
        2479             :     ! Input data is in meters or specifies the PBL top;
        2480             :     ! Find the corresponding level index
        2481           0 :     IF ( altu == HCO_EMISL_M .OR. altu == HCO_EMISL_PBL ) THEN
        2482             :
        2483             :        ! Eventually get altitude from PBL height
        2484           0 :        IF ( altu == HCO_EMISL_PBL ) THEN
        2485           0 :           alt = HcoState%Grid%PBLHEIGHT%Val(I,J)
        2486             :        ENDIF
        2487             :
        2488             :        ! Special case of negative height
        2489           0 :        IF ( alt <= 0.0_hp ) THEN
        2490           0 :           lidx = 1
        2491           0 :           RETURN
        2492             :        ENDIF
        2493             :
        2494             :        ! Loop over data until we are within desired level
        2495             :        ! NOTE: This can be rewritten more efficiently (bmy, 3/5/21)
        2496           0 :        altt = 0.0_hp
        2497           0 :        altb = 0.0_hp
        2498           0 :        lidx = -1
        2499           0 :        DO L = 1, HcoState%NZ
        2500           0 :           altt = altb + HcoState%Grid%BXHEIGHT_M%Val(I,J,L)
        2501           0 :           IF ( alt >= altb .AND. alt < altt ) THEN
        2502           0 :              lidx = L
        2503           0 :              RETURN
        2504             :           ENDIF
        2505           0 :           altb = altt
        2506             :        ENDDO
        2507             :
        2508             :        ! If altitude is above maximum level
        2509           0 :        IF ( lidx == -1 .AND. alt >= altt ) THEN
        2510           0 :           lidx = HcoState%NZ
        2511             :           WRITE(MSG,*)                                                       &
        2512           0 :                'Level is above max. grid box level - use top level ', alt
        2513           0 :           CALL HCO_WARNING ( HcoState%Config%Err, MSG, RC, THISLOC=LOC )
        2514           0 :           RETURN
        2515             :        ENDIF
        2516             :
        2517             :        RETURN
        2518             :     ENDIF
        2519             :
        2520             :     ! Error if we drop down to here
        2521           0 :     MSG = 'Illegal altitude unit'
        2522           0 :     CALL HCO_ERROR ( MSG, RC, THISLOC=LOC )
        2523             :
        2524             :   END SUBROUTINE GetIdx
        2525             : !EOC
        2526             : !------------------------------------------------------------------------------
        2527             : !                   Harmonized Emissions Component (HEMCO)                    !
        2528             : !------------------------------------------------------------------------------
        2529             : !BOP
        2530             : !
        2531             : ! !IROUTINE: GetDilFact
        2532             : !
        2533             : ! !DESCRIPTION: Subroutine GetDilFact returns the vertical dilution factor,
        2534             : ! that is the factor that is to be applied to distribute emissions into
        2535             : ! multiple vertical levels. If grid box height information are available,
        2536             : ! these are used to compute the distribution factor. Otherwise, equal weight
        2537             : ! is given to all vertical levels.
        2538             : !\\
        2539             : !\\
        2540             : ! !TODO: Dilution factors are currently only weighted by grid box heights
        2541             : ! (if these information are available) but any pressure information are
        2542             : ! ignored.
        2543             : !\\
        2544             : !\\
        2545             : ! !INTERFACE:
        2546             : !
        2547           0 :   SUBROUTINE GetDilFact( HcoState,   EmisL1, EmisL1Unit, EmisL2,             &
        2548             :                          EmisL2Unit, I,      J,          L,                  &
        2549             :                          LowLL,      UppLL,  DilFact,    RC                 )
        2550             : !
        2551             : ! !USES:
        2552             : !
        2553             :     USE HCO_STATE_MOD,    ONLY : HCO_State
        2554             : !
        2555             : ! !INPUT PARAMETERS:
        2556             : !
        2557             :     TYPE(HCO_State), POINTER       :: HcoState    ! HEMCO state object
        2558             :     INTEGER,         INTENT(IN)    :: I           ! lon index
        2559             :     INTEGER,         INTENT(IN)    :: J           ! lat index
        2560             :     INTEGER,         INTENT(IN)    :: L           ! lev index
        2561             :     INTEGER,         INTENT(IN)    :: LowLL       ! lower level index
        2562             :     INTEGER,         INTENT(IN)    :: UppLL       ! upper level index
        2563             : !
        2564             : ! !OUTPUT PARAMETERS:
        2565             : !
        2566             :     REAL(hp),        INTENT(OUT)   :: DilFact     ! Dilution factor
        2567             : !
        2568             : ! !INPUT/OUTPUT PARAMETERS:
        2569             : !
        2570             :     REAL(hp),        INTENT(INOUT) :: EmisL1
        2571             :     INTEGER,         INTENT(INOUT) :: EmisL1Unit
        2572             :     REAL(hp),        INTENT(INOUT) :: EmisL2
        2573             :     INTEGER,         INTENT(INOUT) :: EmisL2Unit
        2574             :     INTEGER,         INTENT(INOUT) :: RC
        2575             : !
        2576             : ! !REVISION HISTORY:
        2577             : !  06 May 2016 - C. Keller   - Initial Version
        2578             : !  See https://github.com/geoschem/hemco for complete history
        2579             : !EOP
        2580             : !------------------------------------------------------------------------------
        2581             : !BOC
        2582             : !
        2583             : ! !LOCAL VARIABLES:
        2584             : !
        2585             :     INTEGER            :: L1
        2586             :     CHARACTER(LEN=255) :: MSG
        2587             :     CHARACTER(LEN=255) :: LOC = 'GetDilFact (hco_calc_mod.F90)'
        2588             :     REAL(hp)           :: h1, h2, dh, dh1, dh2
        2589             :     REAL(hp)           :: UppLLR, LowLLR
        2590             :
        2591             :     !=================================================================
        2592             :     ! GetDilFact begins here
        2593             :     !=================================================================
        2594             :
        2595             :     ! Init
        2596           0 :     DilFact = 1.0_hp
        2597           0 :     RC = HCO_SUCCESS
        2598             :
        2599             :     ! Nothing to do if it's only one level
        2600           0 :     IF ( LowLL == UppLL ) RETURN
        2601             :
        2602             :     ! Compute dilution factor based on boxheights if this information
        2603             :     ! is available
        2604           0 :     IF ( ASSOCIATED( HcoState%Grid%BXHEIGHT_M%Val ) ) THEN
        2605             :
        2606             :        ! Get height of bottom level LowLL (in m)
        2607           0 :        IF ( EmisL1Unit == HCO_EMISL_M ) THEN
        2608           0 :           h1 = EmisL1
        2609           0 :        ELSEIF ( EmisL1Unit == HCO_EMISL_PBL ) THEN
        2610           0 :           h1 = HcoState%Grid%PBLHEIGHT%Val(I,J)
        2611             :        ELSE
        2612           0 :           IF ( LowLL > 1 ) THEN
        2613           0 :              h1 = SUM(HcoState%Grid%BXHEIGHT_M%Val(I,J,1:(LowLL-1)))
        2614             :           ELSE
        2615             :              h1 = 0.0_hp
        2616             :           ENDIF
        2617             :        ENDIF
        2618             :
        2619             :        ! Get height of top level UppLL (in m)
        2620           0 :        IF ( EmisL2Unit == HCO_EMISL_M ) THEN
        2621           0 :           h2 = EmisL2
        2622           0 :        ELSEIF ( EmisL2Unit == HCO_EMISL_PBL ) THEN
        2623           0 :           h2 = HcoState%Grid%PBLHEIGHT%Val(I,J)
        2624             :        ELSE
        2625           0 :           h2 = SUM(HcoState%Grid%BXHEIGHT_M%Val(I,J,1:UppLL))
        2626             :        ENDIF
        2627             :
        2628             :        ! If vertical weight option is enabled, calculate vertical
        2629             :        ! distribution factor relative to the grid cell heights. This
        2630             :        ! is the default (and recommended) option as this makes sure
        2631             :        ! that the same amount of mass is emitted into each layer.
        2632           0 :        IF ( HcoState%Options%VertWeight ) THEN
        2633             :
        2634             :           ! Height of grid box of interest (in m)
        2635           0 :           dh = HcoState%Grid%BXHEIGHT_M%Val(I,J,L)
        2636             :
        2637             :           ! Adjust dh if we are in lowest level
        2638           0 :           IF ( L == LowLL ) THEN
        2639           0 :              dh = SUM(HcoState%Grid%BXHEIGHT_M%Val(I,J,1:LowLL)) - h1
        2640             :           ENDIF
        2641             :
        2642             :           ! Adjust dh if we are in top level
        2643           0 :           IF ( L == UppLL ) THEN
        2644           0 :              dh = h2 - SUM(HcoState%Grid%BXHEIGHT_M%Val(I,J,1:(UppLL-1)))
        2645             :           ENDIF
        2646             :
        2647             :           ! compute dilution factor: the new flux should emit the same mass per
        2648             :           ! volume, i.e. flux_total/column_total = flux_level/column_level
        2649             :           ! --> flux_level = fluxtotal * column_level / column_total.
        2650           0 :           IF ( h2 > h1 ) THEN
        2651           0 :              DilFact = dh / ( h2 - h1 )
        2652             :           ELSE
        2653           0 :              MSG = 'GetDilFact h2 not greater than h1'
        2654           0 :              CALL HCO_ERROR ( HcoState%Config%Err, MSG, RC, THISLOC=LOC )
        2655           0 :              RETURN
        2656             :           ENDIF
        2657             :
        2658             :        ! If VertWeight option is turned off, emit the same flux in each layer.
        2659             :        ! Since model layers have different depths, this will result in differnt
        2660             :        ! total emissions per layer.
        2661             :        ELSE
        2662             :
        2663             :           ! Get fractional layer indeces for lower and upper level. This makes
        2664             :           ! sure that only fractions of the lower and upper level are being
        2665             :           ! considered, so that double-counting is avoided if a model layer
        2666             :           ! serves both as the top layer and the bottom layer (e.g., wildfire
        2667             :           ! emissions emitted from bottom to the top of PBL, and from the top
        2668             :           ! of PBL to 5000m).
        2669           0 :           LowLLR = REAL(LowLL,hp) - 1.0_hp
        2670           0 :           UppLLR = REAL(UppLL,hp)
        2671           0 :           dh1 = 0.0_hp
        2672           0 :           DO L1 = 1, HcoState%NZ
        2673           0 :              dh2 = SUM(HcoState%Grid%BXHEIGHT_M%Val(I,J,1:L1))
        2674           0 :              IF ( h1 >= dh1 .AND. h1 < dh2 ) THEN
        2675           0 :                 LowLLR = REAL(L1,hp) - ( (dh2-h1)/(dh2-dh1) )
        2676             :              ENDIF
        2677           0 :              IF ( h2 > dh1 .AND. h2 <= dh2 ) THEN
        2678           0 :                 UppLLR = REAL(L1,hp) - ( (dh2-h2)/(dh2-dh1) )
        2679             :              ENDIF
        2680             :              ! top layer is bottom layer in next loop
        2681           0 :              dh1 = dh2
        2682             :           ENDDO
        2683             :
        2684             :           ! Dilution factor using fractional levels
        2685           0 :           IF ( UppLLR <= LowLLR ) THEN
        2686           0 :              DilFact = 1.0_hp / REAL(UppLL-LowLL+1,hp)
        2687             :           ELSE
        2688           0 :              DilFact = 1.0_hp / (UppLLR-LowLLR)
        2689             :           ENDIF
        2690             :
        2691             :        ENDIF
        2692             :
        2693             :     ! Approximate dilution factor otherwise
        2694             :     ELSE
        2695             :
        2696           0 :        DilFact = 1.0_hp / REAL(UppLL-LowLL+1,hp)
        2697             :     ENDIF
        2698             :
        2699             :     ! Return w/ success
        2700           0 :     RC = HCO_SUCCESS
        2701             :
        2702             :   END SUBROUTINE GetDilFact
        2703             : #ifdef ADJOINT
        2704             : !BOP
        2705             : !
        2706             : ! !IROUTINE: Get_Current_Emissions
        2707             : !
        2708             : ! !DESCRIPTION: Subroutine Get\_Current\_Emissions calculates the current
        2709             : !  emissions for the specified emission container.
        2710             : !  This subroutine is only called by HCO\_CalcEmis and for base emission
        2711             : !  containers, i.e. containers of type 1.
        2712             : !\\
        2713             : !\\
        2714             : ! !INTERFACE:
        2715             : !
        2716             :   SUBROUTINE Get_Current_Emissions_Adj( HcoState, BaseDct,   &
        2717             :                                     nI, nJ, nL, OUTARR_3D, MASK, RC, UseLL )
        2718             : !
        2719             : ! !USES:
        2720             : !
        2721             :     USE HCO_State_Mod,    ONLY : HCO_State
        2722             :     USE HCO_tIdx_MOD,     ONLY : tIDx_GetIndx
        2723             :     USE HCO_FileData_Mod, ONLY : FileData_ArrIsDefined
        2724             : !
        2725             : ! !INPUT PARAMETERS:
        2726             : !
        2727             :     INTEGER,           INTENT(IN)  :: nI                  ! # of lons
        2728             :     INTEGER,           INTENT(IN)  :: nJ                  ! # of lats
        2729             :     INTEGER,           INTENT(IN)  :: nL                  ! # of levs
        2730             : !
        2731             : ! !INPUT/OUTPUT PARAMETERS:
        2732             : !
        2733             :
        2734             :     TYPE(HCO_State), POINTER       :: HcoState            ! HEMCO state object
        2735             :     TYPE(DataCont),  POINTER       :: BaseDct             ! base emission
        2736             :                                                           !  container
        2737             :     REAL(hp),        INTENT(INOUT) :: OUTARR_3D(nI,nJ,nL) ! output array
        2738             :     REAL(hp),        INTENT(INOUT) :: MASK     (nI,nJ,nL) ! mask array
        2739             :     INTEGER,         INTENT(INOUT) :: RC
        2740             : !
        2741             : ! !OUTPUT PARAMETERS:
        2742             : !
        2743             :     INTEGER,         INTENT(  OUT), OPTIONAL :: UseLL
        2744             : !
        2745             : ! !REMARKS:
        2746             : !  This routine uses multiple loops over all grid boxes (base emissions
        2747             : !  and scale factors use separate loops). In an OMP environment, this approach
        2748             : !  seems to be faster than using only one single loop (but repeated calls to
        2749             : !  point to containers, etc.). The alternative approach is used in routine
        2750             : !  Get\_Current\_Emissions\_B at the end of this module and may be employed
        2751             : !  on request.
        2752             : !
        2753             : ! !REVISION HISTORY:
        2754             : !  25 Aug 2012 - C. Keller   - Initial Version
        2755             : !  09 Nov 2012 - C. Keller   - MASK update. Masks are now treated
        2756             : !                              separately so that multiple masks can be
        2757             : !                              added.
        2758             : !  06 Jun 2014 - R. Yantosca - Cosmetic changes in ProTeX headers
        2759             : !  07 Sep 2014 - C. Keller   - Mask update. Now set mask to zero as soon as
        2760             : !                              on of the applied masks is zero.
        2761             : !  03 Dec 2014 - C. Keller   - Now calculate time slice index on-the-fly.
        2762             : !  29 Dec 2014 - C. Keller   - Added scale factor masks.
        2763             : !  02 Mar 2015 - C. Keller   - Now check for missing values. Missing values are
        2764             : !                              excluded from emission calculation.
        2765             : !  26 Oct 2016 - R. Yantosca - Don't nullify local ptrs in declaration stmts
        2766             : !  11 May 2017 - C. Keller   - Added universal scaling
        2767             : !EOP
        2768             : !------------------------------------------------------------------------------
        2769             : !BOC
        2770             : !
        2771             : ! !LOCAL VARIABLES:
        2772             : !
        2773             :     ! Pointers
        2774             :     TYPE(DataCont), POINTER :: ScalDct
        2775             :     TYPE(DataCont), POINTER :: MaskDct
        2776             :     TYPE(DataCont), POINTER :: LevDct1
        2777             :     TYPE(DataCont), POINTER :: LevDct2
        2778             :
        2779             :     ! Scalars
        2780             :     REAL(sp)                :: TMPVAL, MaskScale
        2781             :     REAL(hp)                :: DilFact
        2782             :     REAL(hp)                :: ScalFact
        2783             :     INTEGER                 :: tIDx, IDX
        2784             :     INTEGER                 :: I, J, L, N
        2785             :     INTEGER                 :: LowLL, UppLL, ScalLL, TmpLL
        2786             :     INTEGER                 :: ERROR
        2787             :     INTEGER                 :: TotLL, nnLL
        2788             :     CHARACTER(LEN=255)      :: MSG, LOC
        2789             :     LOGICAL                 :: NegScalExist
        2790             :     LOGICAL                 :: MaskFractions
        2791             :     LOGICAL                 :: isLevDct1
        2792             :     LOGICAL                 :: isLevDct2
        2793             :     LOGICAL                 :: isMaskDct
        2794             :     LOGICAL                 :: isPblHt
        2795             :     LOGICAL                 :: isBoxHt
        2796             :     INTEGER                 :: LevDct1_Unit
        2797             :     INTEGER                 :: LevDct2_Unit
        2798             :
        2799             :     ! testing only
        2800             :     INTEGER                 :: IX, IY
        2801             :
        2802             :     !=================================================================
        2803             :     ! GET_CURRENT_EMISSIONS begins here
        2804             :     !=================================================================
        2805             :
        2806             :     ! Initialize
        2807             :     ScalDct => NULL()
        2808             :     MaskDct => NULL()
        2809             :     LOC     = 'GET_CURRENT_EMISSIONS_ADJ (hco_calc_mod.F90)'
        2810             :
        2811             :     ! Enter
        2812             :     CALL HCO_ENTER(HcoState%Config%Err, LOC, RC )
        2813             :     IF(RC /= HCO_SUCCESS) RETURN
        2814             :
        2815             :     ! testing only:
        2816             :     IX = 3 !-1
        2817             :     IY = 8 !-1
        2818             :
        2819             :     ! Check if container contains data
        2820             :     IF ( .NOT. FileData_ArrIsDefined(BaseDct%Dta) ) THEN
        2821             :        MSG = 'Array not defined: ' // TRIM(BaseDct%cName)
        2822             :        CALL HCO_ERROR( MSG, RC )
        2823             :        RETURN
        2824             :     ENDIF
        2825             :
        2826             :     ! Initialize mask. By default, assume that we use all grid boxes.
        2827             :     MASK(:,:,:)  = 1.0_hp
        2828             :     MaskFractions = HcoState%Options%MaskFractions
        2829             :
        2830             :     ! Verbose
        2831             :     IF ( HCO_IsVerb(HcoState%Config%Err) ) THEN
        2832             :        WRITE(MSG,*) 'Evaluate field ', TRIM(BaseDct%cName)
        2833             :        CALL HCO_MSG(HcoState%Config%Err,MSG,SEP1=' ')
        2834             :     ENDIF
        2835             :
        2836             :     ! ----------------------------------------------------------------
        2837             :     ! Set base emissions
        2838             :     ! ----------------------------------------------------------------
        2839             :
        2840             :     ! Initialize ERROR. Will be set to 1 if error occurs below
        2841             :     ERROR = 0
        2842             :
        2843             :     ! Initialize variables to compute average vertical level index
        2844             :     totLL = 0
        2845             :     nnLL  = 0
        2846             :
        2847             :     ! Check for level index containers
        2848             :     IF ( BaseDct%levScalID1 > 0 ) THEN
        2849             :        CALL Pnt2DataCont( HcoState, BaseDct%levScalID1, LevDct1, RC )
        2850             :        IF ( RC /= HCO_SUCCESS ) THEN
        2851             :            CALL HCO_ERROR( 'ERROR 18', RC, THISLOC=LOC )
        2852             :            RETURN
        2853             :        ENDIF
        2854             :     ELSE
        2855             :        LevDct1 => NULL()
        2856             :     ENDIF
        2857             :     IF ( BaseDct%levScalID2 > 0 ) THEN
        2858             :        CALL Pnt2DataCont( HcoState, BaseDct%levScalID2, LevDct2, RC )
        2859             :        IF ( RC /= HCO_SUCCESS ) THEN
        2860             :            CALL HCO_ERROR( 'ERROR 19', RC, THISLOC=LOC )
        2861             :            RETURN
        2862             :        ENDIF
        2863             :     ELSE
        2864             :        LevDct2 => NULL()
        2865             :     ENDIF
        2866             :
        2867             :     ! Test whether LevDct1 and LevDct2 are associated
        2868             :     isLevDct1 = ASSOCIATED( LevDct1 )
        2869             :     isLevDct2 = ASSOCIATED( LevDct2 )
        2870             :
        2871             :     ! Get the units of LevDct1 (if it exists)
        2872             :     IF ( isLevDct1 ) THEN
        2873             :        LevDct1_Unit = GetEmisLUnit( HcoState, LevDct1 )
        2874             :        IF ( LevDct1_Unit < 0 ) THEN
        2875             :           MSG = 'LevDct1 units are not defined!'
        2876             :           CALL HCO_ERROR ( HcoState%Config%Err, MSG, RC, THISLOC=LOC )
        2877             :           RC = HCO_FAIL
        2878             :           RETURN
        2879             :        ENDIF
        2880             :     ELSE
        2881             :        LevDct1_Unit = -1
        2882             :     ENDIF
        2883             :
        2884             :     ! Get the units of LevDct2 (if it exists)
        2885             :     IF ( isLevDct2 ) THEN
        2886             :        LevDct2_Unit = GetEmisLUnit( HcoState, LevDct2 )
        2887             :        IF ( LevDct2_Unit < 0 ) THEN
        2888             :           MSG = 'LevDct2_Units are not defined!'
        2889             :           CALL HCO_ERROR ( HcoState%Config%Err, MSG, RC, THISLOC=LOC )
        2890             :           RETURN
        2891             :        ENDIF
        2892             :     ELSE
        2893             :        LevDct2_Unit = -1
        2894             :     ENDIF
        2895             :
        2896             :     ! Throw an error if boxheight is missing and the units are in meters
        2897             :     IF ( LevDct1_Unit == HCO_EMISL_M  .or.                                  &
        2898             :          LevDct2_Unit == HCO_EMISL_M ) THEN
        2899             :        IF ( .NOT. ASSOCIATED(HcoState%Grid%BXHEIGHT_M%Val) ) THEN
        2900             :           MSG = 'Boxheight (in meters) is missing in HEMCO state'
        2901             :           CALL HCO_ERROR ( HcoState%Config%Err, MSG, RC, THISLOC=LOC )
        2902             :           RETURN
        2903             :        ENDIF
        2904             :     ENDIF
        2905             :
        2906             :     ! Throw an error if boxheight is missing and the units are in PBL frac
        2907             :     IF ( LevDct1_Unit == HCO_EMISL_PBL  .or.                                &
        2908             :          LevDct2_Unit == HCO_EMISL_PBL ) THEN
        2909             :        IF ( .NOT. ASSOCIATED(HcoState%Grid%PBLHEIGHT%Val) ) THEN
        2910             :           MSG = 'Boundary layer height is missing in HEMCO state'
        2911             :           CALL HCO_ERROR ( HcoState%Config%Err, MSG, RC, THISLOC=LOC )
        2912             :           RETURN
        2913             :        ENDIF
        2914             :     ENDIF
        2915             :
        2916             :     ! Loop over all latitudes and longitudes
        2917             : !$OMP PARALLEL DO                                                            &
        2918             : !$OMP DEFAULT( SHARED                                                       )&
        2919             : !$OMP PRIVATE( I, J, L, tIdx, TMPVAL, DilFact, LowLL, UppLL                 )&
        2920             : !$OMP COLLAPSE( 2                                                           )&
        2921             : !$OMP SCHEDULE( DYNAMIC, 4                                                  )&
        2922             : !$OMP REDUCTION( +:totLL                                                    )&
        2923             : !$OMP REDUCTION( +:nnLL                                                     )
        2924             :     DO J = 1, nJ
        2925             :     DO I = 1, nI
        2926             :
        2927             :        ! Zero for safety's sake
        2928             :        totLL = 0
        2929             :        nnLL  = 0
        2930             :
        2931             :        ! Get current time index for this container and at this location
        2932             :        tIDx = tIDx_GetIndx( HcoState, BaseDct%Dta, I, J )
        2933             :        IF ( tIDx < 1 ) THEN
        2934             :           WRITE(MSG,*) 'Cannot get time slice index at location ',I,J,&
        2935             :                        ': ', TRIM(BaseDct%cName), tIDx
        2936             :           ERROR = 1
        2937             :           EXIT
        2938             :        ENDIF
        2939             :
        2940             :        ! Get lower and upper vertical index
        2941             :        CALL GetVertIndx ( HcoState,     BaseDct,   isLevDct1, LevDct1,       &
        2942             :                           LevDct1_Unit, isLevDct2, LevDct2,   LevDct2_Unit,  &
        2943             :                           I,            J,         LowLL,     UppLL,         &
        2944             :                           RC                                                )
        2945             :        IF ( RC /= HCO_SUCCESS ) THEN
        2946             :           WRITE(MSG,*) 'Error getting vertical index at location ',I,J,&
        2947             :                        ': ', TRIM(BaseDct%cName)
        2948             :           ERROR = 1 ! Will cause error
        2949             :           EXIT
        2950             :        ENDIF
        2951             :
        2952             :        ! average upper level
        2953             :        totLL = totLL + UppLL
        2954             :        nnLL  = nnLL + 1
        2955             :
        2956             :        ! Loop over all levels
        2957             :        DO L = LowLL, UppLL
        2958             :
        2959             :           ! Get base value. Use uniform value if scalar field.
        2960             :           IF ( BaseDct%Dta%SpaceDim == 1 ) THEN
        2961             :              TMPVAL = BaseDct%Dta%V2(tIDx)%Val(1,1)
        2962             :           ELSEIF ( BaseDct%Dta%SpaceDim == 2 ) THEN
        2963             :              TMPVAL = BaseDct%Dta%V2(tIDx)%Val(I,J)
        2964             :           ELSE
        2965             :              TMPVAL = BaseDct%Dta%V3(tIDx)%Val(I,J,L)
        2966             :           ENDIF
        2967             :
        2968             :           ! If it's a missing value, mask box as unused and set value to zero
        2969             :           IF ( TMPVAL == HCO_MISSVAL ) THEN
        2970             :              MASK(I,J,:)      = 0.0_hp
        2971             :              OUTARR_3D(I,J,L) = 0.0_hp
        2972             :
        2973             :           ! Pass base value to output array
        2974             :           ELSE
        2975             :
        2976             :              ! Get dilution factor. Never dilute 3D emissions.
        2977             :              IF ( BaseDct%Dta%SpaceDim == 3 ) THEN
        2978             :                 DilFact = 1.0_hp !1.0
        2979             :
        2980             :              ! 2D dilution factor
        2981             :              ELSE
        2982             :                 CALL GetDilFact ( HcoState,    BaseDct%Dta%EmisL1, &
        2983             :                                   BaseDct%Dta%EmisL1Unit, BaseDct%Dta%EmisL2,  &
        2984             :                                   BaseDct%Dta%EmisL2Unit, I, J, L, LowLL,  &
        2985             :                                   UppLL, DilFact, RC )
        2986             :                 IF ( RC /= HCO_SUCCESS ) THEN
        2987             :                    WRITE(MSG,*) 'Error getting dilution factor at ',I,J,&
        2988             :                                 ': ', TRIM(BaseDct%cName)
        2989             :                    ERROR = 1
        2990             :                    EXIT
        2991             :                 ENDIF
        2992             :              ENDIF
        2993             :
        2994             :              ! Scale base emission by dilution factor
        2995             :              OUTARR_3D(I,J,L) = DilFact * TMPVAL
        2996             :           ENDIF
        2997             :        ENDDO !L
        2998             :
        2999             :     ENDDO !I
        3000             :     ENDDO !J
        3001             : !$OMP END PARALLEL DO
        3002             :
        3003             :     ! Check for error
        3004             :     IF ( ERROR == 1 ) THEN
        3005             :        CALL HCO_ERROR( MSG, RC )
        3006             :        RETURN
        3007             :     ENDIF
        3008             :
        3009             :     ! ----------------------------------------------------------------
        3010             :     ! Apply scale factors
        3011             :     ! The container IDs of all scale factors associated with this base
        3012             :     ! container are stored in vector Scal_cID.
        3013             :     ! ----------------------------------------------------------------
        3014             :
        3015             :     ! Loop over scale factors
        3016             :     IF ( BaseDct%nScalID > 0 ) THEN
        3017             :
        3018             :     DO N = 1, BaseDct%nScalID
        3019             :
        3020             :        ! Get the scale factor container ID for the current slot
        3021             :        IDX = BaseDct%Scal_cID(N)
        3022             :
        3023             :        ! Point to data container with the given container ID
        3024             :        CALL Pnt2DataCont( HcoState, IDX, ScalDct, RC )
        3025             :        IF ( RC /= HCO_SUCCESS ) THEN
        3026             :            CALL HCO_ERROR( 'ERROR 20', RC, THISLOC=LOC )
        3027             :            RETURN
        3028             :        ENDIF
        3029             :
        3030             :        ! Sanity check: scale field cannot be a base field
        3031             :        IF ( (ScalDct%DctType == HCO_DCTTYPE_BASE) ) THEN
        3032             :           MSG = 'Wrong scale field type: ' // TRIM(ScalDct%cName)
        3033             :           CALL HCO_ERROR( MSG, RC )
        3034             :           RETURN
        3035             :        ENDIF
        3036             :
        3037             :        ! Skip this scale factor if no data defined. This is possible
        3038             :        ! if scale factors are only defined for a given time range and
        3039             :        ! the simulation datetime is outside of this range.
        3040             :        IF ( .NOT. FileData_ArrIsDefined(ScalDct%Dta) ) THEN
        3041             :           IF ( HCO_IsVerb( HcoState%Config%Err ) ) THEN
        3042             :              MSG = 'Skip scale factor '//TRIM(ScalDct%cName)// &
        3043             :                    ' because it is not defined for this datetime.'
        3044             :              CALL HCO_MSG(HcoState%Config%Err,MSG)
        3045             :           ENDIF
        3046             :           CYCLE
        3047             :        ENDIF
        3048             :
        3049             :        ! Verbose mode
        3050             :        IF ( HCO_IsVerb(HcoState%Config%Err ) ) THEN
        3051             :           MSG = 'Applying scale factor ' // TRIM(ScalDct%cName)
        3052             :           CALL HCO_MSG(HcoState%Config%Err,MSG)
        3053             :        ENDIF
        3054             :
        3055             :        ! Get vertical extension of this scale factor array.
        3056             :        IF( (ScalDct%Dta%SpaceDim<=2) ) THEN
        3057             :           ScalLL = 1
        3058             :        ELSE
        3059             :           ScalLL = SIZE(ScalDct%Dta%V3(1)%Val,3)
        3060             :        ENDIF
        3061             :
        3062             :        ! Check if there is a mask field associated with this scale
        3063             :        ! factor. In this case, get a pointer to the corresponding
        3064             :        ! mask field and evaluate scale factors only inside the mask
        3065             :        ! region.
        3066             :        IF ( ASSOCIATED(ScalDct%Scal_cID) ) THEN
        3067             :           CALL Pnt2DataCont( HcoState, ScalDct%Scal_cID(1), MaskDct, RC )
        3068             :           IF ( RC /= HCO_SUCCESS ) THEN
        3069             :               CALL HCO_ERROR( 'ERROR 21', RC, THISLOC=LOC )
        3070             :               RETURN
        3071             :           ENDIF
        3072             :
        3073             :           ! Must be mask field
        3074             :           IF ( MaskDct%DctType /= HCO_DCTTYPE_MASK ) THEN
        3075             :              MSG = 'Invalid mask for scale factor: '//TRIM(ScalDct%cName)
        3076             :              MSG = TRIM(MSG) // '; mask: '//TRIM(MaskDct%cName)
        3077             :              CALL HCO_ERROR( MSG, RC )
        3078             :              RETURN
        3079             :           ENDIF
        3080             :        ENDIF
        3081             :
        3082             :        ! Reinitialize error flag. Will be set to 1 or 2 if error occurs,
        3083             :        ! and to -1 if negative scale factor is ignored.
        3084             :        ERROR = 0
        3085             :
        3086             :        ! Loop over all latitudes and longitudes
        3087             : !$OMP PARALLEL DO                                                            &
        3088             : !$OMP DEFAULT( SHARED )                                                      &
        3089             : !$OMP PRIVATE( I, J, tIdx, TMPVAL, L, LowLL, UppLL, tmpLL, MaskScale        )&
        3090             : !$OMP COLLAPSE( 2                                                           )&
        3091             : !$OMP SCHEDULE( DYNAMIC, 4                                                  )
        3092             :        DO J = 1, nJ
        3093             :        DO I = 1, nI
        3094             :
        3095             :           ! ------------------------------------------------------------
        3096             :           ! If there is a mask associated with this scale factors, check
        3097             :           ! if this grid box is within or outside of the mask region.
        3098             :           ! Values that partially fall into the mask region are either
        3099             :           ! treated as binary (100% inside or outside), or partially
        3100             :           ! (using the real grid area fractions), depending on the
        3101             :           ! HEMCO options.
        3102             :           ! ------------------------------------------------------------
        3103             :
        3104             :           ! Default mask scaling is 1.0 (no mask applied)
        3105             :           MaskScale = 1.0_sp
        3106             :
        3107             :           ! If there is a mask applied to this scale factor ...
        3108             :           IF ( ASSOCIATED(MaskDct) ) THEN
        3109             :              CALL GetMaskVal ( MaskDct, I, J, &
        3110             :                                MaskScale, MaskFractions, RC )
        3111             :              IF ( RC /= HCO_SUCCESS ) THEN
        3112             :                 ERROR = 4
        3113             :                 EXIT
        3114             :              ENDIF
        3115             :           ENDIF
        3116             :
        3117             :           ! We can skip this grid box if mask is completely zero
        3118             :           IF ( MaskScale <= 0.0_sp ) CYCLE
        3119             :
        3120             :           ! Get current time index for this container and at this location
        3121             :           tIDx = tIDx_GetIndx( HcoState, ScalDct%Dta, I, J )
        3122             :           IF ( tIDx < 1 ) THEN
        3123             :              WRITE(*,*) 'Cannot get time slice index at location ',I,J,&
        3124             :                           ': ', TRIM(ScalDct%cName), tIDx
        3125             :              ERROR = 3
        3126             :              EXIT
        3127             :           ENDIF
        3128             :
        3129             :           ! Check if this is a mask. If so, add mask values to the MASK
        3130             :           ! array. For now, we assume masks to be binary, i.e. 0 or 1.
        3131             :           ! We may want to change that in future to also support values
        3132             :           ! in between. This is especially important when regridding
        3133             :           ! high resolution masks onto coarser grids!
        3134             :           ! ------------------------------------------------------------
        3135             :           IF ( ScalDct%DctType == HCO_DCTTYPE_MASK ) THEN
        3136             :
        3137             :              ! Get mask value
        3138             :              CALL GetMaskVal ( ScalDct, I, J, &
        3139             :                                TMPVAL,    MaskFractions, RC )
        3140             :              IF ( RC /= HCO_SUCCESS ) THEN
        3141             :                 ERROR = 4
        3142             :                 EXIT
        3143             :              ENDIF
        3144             :
        3145             :              ! Pass to output mask
        3146             :              MASK(I,J,:) = MASK(I,J,:) * TMPVAL
        3147             :
        3148             :              ! testing only
        3149             :              IF ( HCO_IsVerb(HcoState%Config%Err) .AND. I==1 .AND. J==1 ) THEN
        3150             :                 write(MSG,*) 'Mask field ', TRIM(ScalDct%cName),   &
        3151             :                      ' found and added to temporary mask.'
        3152             :                 CALL HCO_MSG(HcoState%Config%Err,MSG)
        3153             :              ENDIF
        3154             :
        3155             :              ! Advance to next grid box
        3156             :              CYCLE
        3157             :           ENDIF! DctType=MASK
        3158             :
        3159             :           ! ------------------------------------------------------------
        3160             :           ! For non-mask fields, apply scale factors to all levels
        3161             :           ! of the base field individually. If the scale factor
        3162             :           ! field has more than one vertical level, use the
        3163             :           ! vertical level closest to the corresponding vertical
        3164             :           ! level of the base emission field
        3165             :           ! ------------------------------------------------------------
        3166             :
        3167             :           ! Get lower and upper vertical index
        3168             :           CALL GetVertIndx( HcoState, BaseDct,       isLevDct1,              &
        3169             :                             LevDct1,  LevDct1_Unit,  isLevDct2,              &
        3170             :                             LevDct2,  LevDct2_Unit,  I,                      &
        3171             :                             J,        LowLL,         UppLL,      RC         )
        3172             :           IF ( RC /= HCO_SUCCESS ) THEN
        3173             :              ERROR = 1 ! Will cause error
        3174             :              EXIT
        3175             :           ENDIF
        3176             :
        3177             :           ! Loop over all vertical levels of the base field
        3178             :           DO L = LowLL,UppLL
        3179             :              ! If the vertical level exceeds the number of available
        3180             :              ! scale factor levels, use the highest available level.
        3181             :              IF ( L > ScalLL ) THEN
        3182             :                 TmpLL = ScalLL
        3183             :              ! Otherwise use the same vertical level index.
        3184             :              ELSE
        3185             :                 TmpLL = L
        3186             :              ENDIF
        3187             :
        3188             :              ! Get scale factor for this grid box. Use same uniform
        3189             :              ! value if it's a scalar field
        3190             :              IF ( ScalDct%Dta%SpaceDim == 1 ) THEN
        3191             :                 TMPVAL = ScalDct%Dta%V2(tidx)%Val(1,1)
        3192             :              ELSEIF ( ScalDct%Dta%SpaceDim == 2 ) THEN
        3193             :                 TMPVAL = ScalDct%Dta%V2(tidx)%Val(I,J)
        3194             :              ELSE
        3195             :                 TMPVAL = ScalDct%Dta%V3(tidx)%Val(I,J,TmpLL)
        3196             :              ENDIF
        3197             :
        3198             :              ! Set missing value to one
        3199             :              IF ( TMPVAL == HCO_MISSVAL ) TMPVAL = 1.0_sp
        3200             :
        3201             :              ! Eventually apply mask scaling
        3202             :              IF ( MaskScale /= 1.0_sp ) THEN
        3203             :                 TMPVAL = TMPVAL * MaskScale
        3204             :              ENDIF
        3205             :
        3206             :              ! For negative scale factor, proceed according to the
        3207             :              ! negative value setting specified in the configuration
        3208             :              ! file (NegFlag = 2: use this value):
        3209             :              IF ( TMPVAL < 0.0_sp .AND. HcoState%Options%NegFlag /= 2 ) THEN
        3210             :
        3211             :                 ! NegFlag = 1: ignore and show warning
        3212             :                 IF ( HcoState%Options%NegFlag == 1 ) THEN
        3213             :                    ERROR = -1 ! Will prompt warning
        3214             :                    CYCLE
        3215             :
        3216             :                 ! Return w/ error otherwise
        3217             :                 ELSE
        3218             :                    WRITE(*,*) 'Negative scale factor at ',I,J,TmpLL,tidx,&
        3219             :                               ': ', TRIM(ScalDct%cName), TMPVAL
        3220             :                    ERROR = 1 ! Will cause error
        3221             :                    EXIT
        3222             :                 ENDIF
        3223             :              ENDIF
        3224             :
        3225             :              ! -------------------------------------------------------
        3226             :              ! Apply scale factor in accordance to field operator
        3227             :              ! -------------------------------------------------------
        3228             :
        3229             :              ! Oper 1: multiply
        3230             :              IF ( ScalDct%Oper == 1 ) THEN
        3231             :                 OUTARR_3D(I,J,L) = OUTARR_3D(I,J,L) * TMPVAL
        3232             :
        3233             :              ! Oper -1: divide
        3234             :              ELSEIF ( ScalDct%Oper == -1 ) THEN
        3235             :                 ! Ignore zeros to avoid NaN
        3236             :                 IF ( TMPVAL /= 0.0_sp ) THEN
        3237             :                    OUTARR_3D(I,J,L) = OUTARR_3D(I,J,L) / TMPVAL
        3238             :                 ENDIF
        3239             :
        3240             :              ! Oper 2: square
        3241             :              ELSEIF ( ScalDct%Oper == 2 ) THEN
        3242             :                 OUTARR_3D(I,J,L) = OUTARR_3D(I,J,L) * TMPVAL * TMPVAL
        3243             :
        3244             :              ! Return w/ error otherwise (Oper 3 is only allowed for masks!)
        3245             :              ELSE
        3246             :                 WRITE(*,*) 'Illegal operator for ', TRIM(ScalDct%cName), ScalDct%Oper
        3247             :                 ERROR = 2 ! Will cause error
        3248             :                 EXIT
        3249             :              ENDIF
        3250             :
        3251             :           ENDDO !LL
        3252             :
        3253             :           ! Verbose mode
        3254             :           if ( HCO_IsVerb(HcoState%Config%Err) .and. i == ix .and. j == iy ) then
        3255             :              write(MSG,*) 'Scale field ', TRIM(ScalDct%cName)
        3256             :              CALL HCO_MSG(HcoState%Config%Err,MSG)
        3257             :              write(MSG,*) 'Time slice: ', tIdx
        3258             :              CALL HCO_MSG(HcoState%Config%Err,MSG)
        3259             :              write(MSG,*) 'IX, IY: ', IX, IY
        3260             :              CALL HCO_MSG(HcoState%Config%Err,MSG)
        3261             :              write(MSG,*) 'Scale factor (IX,IY,L1): ', TMPVAL
        3262             :              CALL HCO_MSG(HcoState%Config%Err,MSG)
        3263             :              write(MSG,*) 'Mathematical operation : ', ScalDct%Oper
        3264             :              CALL HCO_MSG(HcoState%Config%Err,MSG)
        3265             : !             write(lun,*) 'Updt (IX,IY,L1): ', OUTARR_3D(IX,IY,1)
        3266             :           endif
        3267             :
        3268             :        ENDDO !I
        3269             :        ENDDO !J
        3270             : !$OMP END PARALLEL DO
        3271             :
        3272             :        ! error check
        3273             :        IF ( ERROR > 0 ) THEN
        3274             :           IF ( ERROR == 1 ) THEN
        3275             :              MSG = 'Negative scale factor found (aborted): ' // TRIM(ScalDct%cName)
        3276             :           ELSEIF ( ERROR == 2 ) THEN
        3277             :              MSG = 'Illegal mathematical operator for scale factor: ' // TRIM(ScalDct%cName)
        3278             :           ELSEIF ( ERROR == 3 ) THEN
        3279             :              MSG = 'Encountered negative time index for scale factor: ' // TRIM(ScalDct%cName)
        3280             :           ELSEIF ( ERROR == 4 ) THEN
        3281             :              MSG = 'Mask error in ' // TRIM(ScalDct%cName)
        3282             :           ELSE
        3283             :              MSG = 'Error when applying scale factor: ' // TRIM(ScalDct%cName)
        3284             :           ENDIF
        3285             :           ScalDct => NULL()
        3286             :           CALL HCO_ERROR( MSG, RC )
        3287             :           RETURN
        3288             :        ENDIF
        3289             :
        3290             :        ! eventually prompt warning for negative values
        3291             :        IF ( ERROR == -1 ) THEN
        3292             :           MSG = 'Negative scale factor found (ignored): ' // TRIM(ScalDct%cName)
        3293             :           CALL HCO_WARNING( HcoState%Config%Err, MSG, RC )
        3294             :        ENDIF
        3295             :
        3296             :        ! Free pointer
        3297             :        MaskDct => NULL()
        3298             :
        3299             :     ENDDO ! N
        3300             :     ENDIF ! N > 0
        3301             :
        3302             :     ! Update optional variables
        3303             :     IF ( PRESENT(UseLL) ) THEN
        3304             :        UseLL = 1
        3305             :        IF ( nnLL > 0 ) UseLL = NINT(REAL(TotLL,kind=sp)/REAL(nnLL,kind=sp))
        3306             :     ENDIF
        3307             :
        3308             :     ! Weight output emissions by mask
        3309             :     OUTARR_3D = OUTARR_3D * MASK
        3310             :
        3311             :     ! Cleanup and leave w/ success
        3312             :     ScalDct => NULL()
        3313             :     CALL HCO_LEAVE ( HcoState%Config%Err, RC )
        3314             :
        3315             :   END SUBROUTINE Get_Current_Emissions_Adj
        3316             : !EOC
        3317             : #endif
        3318             : !EOC
        3319             : END MODULE HCO_Calc_Mod

Generated by: LCOV version 1.14