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

              Line data    Source code
           1             : !------------------------------------------------------------------------------
           2             : !                   Harmonized Emissions Component (HEMCO)                    !
           3             : !------------------------------------------------------------------------------
           4             : !BOP
           5             : !
           6             : ! !MODULE: hco_geotools_mod.F90
           7             : !
           8             : ! !DESCRIPTION: Module HCO\_GeoTools\_Mod contains a collection of
           9             : ! helper routines for extracting geographical information. These
          10             : ! routines are based upon GEOS-5 data and may need to be revised
          11             : ! for other met. fields!
          12             : !\\
          13             : !\\
          14             : ! !INTERFACE:
          15             : !
          16             : MODULE HCO_GeoTools_Mod
          17             : !
          18             : ! !USES:
          19             : !
          20             :   USE HCO_Error_Mod
          21             :
          22             :   IMPLICIT NONE
          23             :   PRIVATE
          24             : !
          25             : ! !PUBLIC MEMBER FUNCTIONS:
          26             : !
          27             :   PUBLIC :: HCO_LandType
          28             :   PUBLIC :: HCO_ValidateLon
          29             :   PUBLIC :: HCO_GetSUNCOS
          30             :   PUBLIC :: HCO_GetHorzIJIndex
          31             :   PUBLIC :: HCO_CalcVertGrid
          32             :   PUBLIC :: HCO_SetPBLm
          33             : !  PUBLIC :: HCO_CalcPBLlev
          34             :
          35             :   INTERFACE HCO_LandType
          36             :      MODULE PROCEDURE HCO_LandType_Dp
          37             :      MODULE PROCEDURE HCO_LandType_Sp
          38             :   END INTERFACE HCO_LandType
          39             :
          40             :   INTERFACE HCO_ValidateLon
          41             :      MODULE PROCEDURE HCO_ValidateLon_Dp
          42             :      MODULE PROCEDURE HCO_ValidateLon_Sp
          43             :   END INTERFACE HCO_ValidateLon
          44             :
          45             : !  INTERFACE HCO_CalcPBLlev
          46             : !     MODULE PROCEDURE HCO_CalcPBLlev2D
          47             : !     MODULE PROCEDURE HCO_CalcPBLlev3D
          48             : !  END INTERFACE HCO_CalcPBLlev
          49             : !
          50             : ! !PRIVATE MEMBER FUNCTIONS:
          51             : !
          52             :   PRIVATE:: HCO_LandType_Dp
          53             :   PRIVATE:: HCO_LandType_Sp
          54             :   PRIVATE:: HCO_ValidateLon_Dp
          55             :   PRIVATE:: HCO_ValidateLon_Sp
          56             : !
          57             : ! !REVISION HISTORY:
          58             : !  18 Dec 2013 - C. Keller   - Initialization
          59             : !  See https://github.com/geoschem/hemco for complete history
          60             : !EOP
          61             : !------------------------------------------------------------------------------
          62             : !BOC
          63             : CONTAINS
          64             : !EOC
          65             : !------------------------------------------------------------------------------
          66             : !                   Harmonized Emissions Component (HEMCO)                    !
          67             : !------------------------------------------------------------------------------
          68             : !BOP
          69             : !
          70             : ! !IROUTINE: HCO_LandType_Sp
          71             : !
          72             : ! !DESCRIPTION: Function HCO\_LANDTYPE returns the land type based upon
          73             : ! GMAO land type fractions. Result is 0=water, 1=land, 2=ice, where ice includes
          74             : ! over both ocean and land, and land includes lakes. Inputs are in single precision.
          75             : !\\
          76             : !\\
          77             : ! !INTERFACE:
          78             : !
          79           0 :   FUNCTION HCO_LandType_Sp( FRLAND, FRLANDIC, FROCEAN, FRSEAICE, FRLAKE ) &
          80             :                           Result ( LandType )
          81             : !
          82             : ! !INPUT PARAMETERS:
          83             : !
          84             :     REAL(sp), INTENT(IN) :: FRLAND    ! Fraction of grid-box covered in land
          85             :     REAL(sp), INTENT(IN) :: FRLANDIC  ! Fraction of grid-box covered in land ice
          86             :     REAL(sp), INTENT(IN) :: FROCEAN   ! Fraction of grid-box covered in ocean
          87             :     REAL(sp), INTENT(IN) :: FRSEAICE  ! Fraction of grid-box covered in sea ice
          88             :     REAL(sp), INTENT(IN) :: FRLAKE    ! Fraction of grid-box covered in lake
          89             : !
          90             : ! !RETURN VALUE
          91             : !
          92             :     INTEGER              :: LandType  ! Land type: 0=ocean, 1=land, 2=ice (ocn,land)
          93             : !
          94             : ! !REMARKS:
          95             : ! Land type index is based on GMAO field LWI, with modification
          96             : ! to classify land ice as ice.
          97             : !
          98             : ! !REVISION HISTORY:
          99             : !  18 Dec 2013 - C. Keller - Initialization!
         100             : !  See https://github.com/geoschem/hemco for complete history
         101             : !EOP
         102             : !------------------------------------------------------------------------------
         103             : !BOC
         104             : !
         105             : ! !DEFINED PARAMETERS:
         106             : !
         107             :     ! Threshold at which a grid-box is considered ice
         108             :     REAL(sp), PARAMETER :: frac_classify_land_ice_as_ice = 0.5_sp
         109             :
         110             :     !--------------------------
         111             :     ! HCO_LANDTYPE begins here
         112             :     !--------------------------
         113             :
         114             :     ! Start with the surface type category definitions based on the GMAO
         115             :     ! definition for land-water-ice index (LWI), which is 0=ocean (non-ice),
         116             :     ! 1=land (includes lakes and ice), 2=ice (over ocean only). Qualify
         117             :     ! land type as location of maximum fraction.
         118             :     LandType = MAXLOC( (/ FRLAND + FRLANDIC + FRLAKE, &
         119             :                           FRSEAICE,                   &
         120           0 :                           FROCEAN - FRSEAICE /), 1 )
         121           0 :     IF ( LandType == 3 ) LandType = 0
         122             :
         123             :     ! Change land type to ice if sufficient ice over land
         124           0 :     IF ( FRLANDIC >= frac_classify_land_ice_as_ice ) THEN
         125           0 :        LandType = 2
         126             :     ENDIF
         127             :
         128           0 :   END FUNCTION HCO_LandType_Sp
         129             : !EOC
         130             : !------------------------------------------------------------------------------
         131             : !                   Harmonized Emissions Component (HEMCO)                    !
         132             : !------------------------------------------------------------------------------
         133             : !BOP
         134             : !
         135             : ! !IROUTINE: HCO_LandType_Dp
         136             : !
         137             : ! !DESCRIPTION: Function HCO\_LANDTYPE returns the land type based upon
         138             : ! GMAO land type fractions. Result is 0=water, 1=land, 2=ice, where ice includes
         139             : ! over both ocean and land, and land includes lakes. Inputs are in double precision.
         140             : !\\
         141             : !\\
         142             : ! !INTERFACE:
         143             : !
         144           0 :   FUNCTION HCO_LandType_Dp( FRLAND, FRLANDIC, FROCEAN, FRSEAICE, FRLAKE ) &
         145             :                           Result ( LandType )
         146             : !
         147             : ! !INPUT PARAMETERS:
         148             : !
         149             :     REAL(dp), INTENT(IN) :: FRLANDIC  ! Fraction of grid-box covered in land ice
         150             :     REAL(dp), INTENT(IN) :: FRLAND    ! Fraction of grid-box covered in land
         151             :     REAL(dp), INTENT(IN) :: FROCEAN   ! Fraction of grid-box covered in ocean
         152             :     REAL(dp), INTENT(IN) :: FRSEAICE  ! Fraction of grid-box covered in sea ice
         153             :     REAL(dp), INTENT(IN) :: FRLAKE    ! Fraction of grid-box covered in lake
         154             : !
         155             : ! !RETURN VALUE:
         156             : !
         157             :     INTEGER              :: LandType  ! Land type: 0=ocean, 1=land, 2=ice (ocn,land)
         158             : !
         159             : ! !REMARKS:
         160             : ! Land type index is based on GMAO field LWI, with modification
         161             : ! to classify land ice as ice.
         162             : !
         163             : ! !REVISION HISTORY:
         164             : !  18 Dec 2013 - C. Keller - Initialization
         165             : !  See https://github.com/geoschem/hemco for complete history
         166             : !EOP
         167             : !------------------------------------------------------------------------------
         168             : !BOC
         169             : !
         170             : ! !DEFINED PARAMETERS::
         171             : !
         172             :     ! Threshold at which a grid-box with land ice is considered ice
         173             :     REAL(dp), PARAMETER :: frac_classify_land_ice_as_ice = 0.5_dp
         174             :
         175             :     !--------------------------
         176             :     ! HCO_LANDTYPE begins here
         177             :     !--------------------------
         178             :
         179             :     ! Start with the surface type category definitions based on the GMAO
         180             :     ! definition for land-water-ice index (LWI), which is 0=ocean (non-ice),
         181             :     ! 1=land (includes lakes and ice), 2=ice (over ocean only). Qualify
         182             :     ! land type as location of maximum fraction.
         183             :     LandType = MAXLOC( (/ FRLAND + FRLANDIC + FRLAKE, &
         184             :                           FRSEAICE,                   &
         185           0 :                           FROCEAN - FRSEAICE /), 1 )
         186           0 :     IF ( LandType == 3 ) LandType = 0
         187             :
         188             :     ! Change land type to ice if sufficient ice over land
         189           0 :     IF ( FRLANDIC >= frac_classify_land_ice_as_ice ) THEN
         190           0 :        LandType = 2
         191             :     ENDIF
         192             :
         193           0 :   END FUNCTION HCO_LandType_Dp
         194             : !EOC
         195             : !------------------------------------------------------------------------------
         196             : !                   Harmonized Emissions Component (HEMCO)                    !
         197             : !------------------------------------------------------------------------------
         198             : !BOP
         199             : !
         200             : ! !IROUTINE: HCO_ValidateLon_Sp
         201             : !
         202             : ! !DESCRIPTION: Subroutine HCO\_ValidateLon\_Sp ensures that the passed
         203             : ! single precision longitude axis LON is steadily increasing.
         204             : !\\
         205             : !\\
         206             : ! !INTERFACE:
         207             : !
         208           0 :   SUBROUTINE HCO_ValidateLon_Sp ( HcoState, NLON, LON, RC )
         209             : !
         210             : ! !USES:
         211             : !
         212             :     USE HCO_STATE_MOD,   ONLY : HCO_STATE
         213             : !
         214             : ! !INPUT/OUTPUT PARAMETERS:
         215             : !
         216             :     TYPE(HCO_State), POINTER       :: HcoState       ! HEMCO state object
         217             :     INTEGER,         INTENT(IN   ) :: NLON        ! # of lons
         218             : !
         219             : ! !INPUT/OUTPUT PARAMETERS:
         220             : !
         221             :     REAL(sp), INTENT(INOUT) :: LON(NLON)   ! longitude axis
         222             :     INTEGER,  INTENT(INOUT) :: RC          ! Return code
         223             : !
         224             : ! !REVISION HISTORY:
         225             : !  16 Jul 2014 - C. Keller - Initialization
         226             : !  See https://github.com/geoschem/hemco for complete history
         227             : !EOP
         228             : !------------------------------------------------------------------------------
         229             : !BOC
         230             : !
         231             : ! !LOCAL VARIABLES:
         232             : !
         233             :     INTEGER             :: I, CNT
         234             :     LOGICAL             :: REDO
         235             :     INTEGER, PARAMETER  :: MAXIT = 10
         236             :
         237             :     !--------------------------
         238             :     ! HCO_ValidateLon_Sp begins here
         239             :     !--------------------------
         240             :
         241           0 :     REDO = .TRUE.
         242           0 :     CNT  = 0
         243             :
         244           0 :     DO WHILE ( REDO )
         245             :
         246             :        ! Exit w/ error after 10 iterations
         247           0 :        CNT = CNT + 1
         248           0 :        IF ( CNT > MAXIT ) THEN
         249             :           CALL HCO_ERROR ( '>10 iterations', RC, &
         250           0 :                            THISLOC='HCO_ValidateLon (HCO_GEOTOOLS_MOD.F90)' )
         251           0 :           RETURN
         252             :        ENDIF
         253             :
         254           0 :        DO I = 1, NLON
         255             :
         256             :           ! If we reach the last grid box, all values are steadily
         257             :           ! increasing (otherwise, the loop would have been exited).
         258           0 :           IF ( I == NLON ) THEN
         259             :              REDO = .FALSE.
         260             :              EXIT
         261             :           ENDIF
         262             :
         263             :           ! Check if next lon value is lower, in which case we subtract
         264             :           ! a value of 360 (degrees) from all longitude values up to
         265             :           ! this point. Then repeat the lookup (from the beginning).
         266           0 :           IF ( LON(I+1) < LON(I) ) THEN
         267           0 :              LON(1:I) = LON(1:I) - 360.0_sp
         268             :              EXIT
         269             :           ENDIF
         270             :
         271             :        ENDDO !I
         272             :     ENDDO ! REDO
         273             :
         274             :     ! Leave w/ success
         275           0 :     RC = HCO_SUCCESS
         276             :
         277             :   END SUBROUTINE HCO_ValidateLon_Sp
         278             : !EOC
         279             : !------------------------------------------------------------------------------
         280             : !                   Harmonized Emissions Component (HEMCO)                    !
         281             : !------------------------------------------------------------------------------
         282             : !BOP
         283             : !
         284             : ! !IROUTINE: HCO_ValidateLon_Dp
         285             : !
         286             : ! !DESCRIPTION: Subroutine HCO\_ValidateLon\_Sp ensures that the passed
         287             : ! double precision longitude axis LON is steadily increasing.
         288             : !\\
         289             : !\\
         290             : ! !INTERFACE:
         291             : !
         292           0 :   SUBROUTINE HCO_ValidateLon_Dp ( HcoState, NLON, LON, RC )
         293             : !
         294             : ! !USES:
         295             : !
         296             :     USE HCO_STATE_MOD,   ONLY : HCO_STATE
         297             : !
         298             : ! !INPUT/OUTPUT PARAMETERS:
         299             : !
         300             :     TYPE(HCO_State), POINTER       :: HcoState       ! HEMCO state object
         301             :     INTEGER,         INTENT(IN   ) :: NLON        ! # of lons
         302             : !
         303             : !
         304             : ! !INPUT/OUTPUT PARAMETERS:
         305             : !
         306             :     REAL(dp), INTENT(INOUT) :: LON(NLON)   ! longitude axis
         307             :     INTEGER,  INTENT(INOUT) :: RC          ! Return code
         308             : !
         309             : ! !REVISION HISTORY:
         310             : !  16 Jul 2014 - C. Keller - Initialization
         311             : !  See https://github.com/geoschem/hemco for complete history
         312             : !EOP
         313             : !------------------------------------------------------------------------------
         314             : !BOC
         315             :     INTEGER             :: I, CNT
         316             :     LOGICAL             :: REDO
         317             :     INTEGER, PARAMETER  :: MAXIT = 10
         318             :
         319             :     !--------------------------
         320             :     ! HCO_ValidateLon_Dp begins here
         321             :     !--------------------------
         322             :
         323           0 :     REDO = .TRUE.
         324           0 :     CNT  = 0
         325             :
         326           0 :     DO WHILE ( REDO )
         327             :
         328             :        ! Exit w/ error after 10 iterations
         329           0 :        CNT = CNT + 1
         330           0 :        IF ( CNT > MAXIT ) THEN
         331             :           CALL HCO_ERROR ( '>10 iterations', RC, &
         332           0 :                            THISLOC='HCO_ValidateLon (HCO_GEOTOOLS_MOD.F90)' )
         333           0 :           RETURN
         334             :        ENDIF
         335             :
         336           0 :        DO I = 1, NLON
         337             :
         338             :           ! If we reach the last grid box, all values are steadily
         339             :           ! increasing (otherwise, the loop would have been exited).
         340           0 :           IF ( I == NLON ) THEN
         341             :              REDO = .FALSE.
         342             :              EXIT
         343             :           ENDIF
         344             :
         345             :           ! Check if next lon value is lower, in which case we subtract
         346             :           ! a value of 360 (degrees) from all longitude values up to
         347             :           ! this point. Then repeat the lookup (from the beginning).
         348           0 :           IF ( LON(I+1) < LON(I) ) THEN
         349           0 :              LON(1:I) = LON(1:I) - 360.0_dp
         350             :              EXIT
         351             :           ENDIF
         352             :
         353             :        ENDDO !I
         354             :     ENDDO ! REDO
         355             :
         356             :     ! Leave w/ success
         357           0 :     RC = HCO_SUCCESS
         358             :
         359             :   END SUBROUTINE HCO_ValidateLon_Dp
         360             : !EOC
         361             : !------------------------------------------------------------------------------
         362             : !                   Harmonized Emissions Component (HEMCO)                    !
         363             : !------------------------------------------------------------------------------
         364             : !BOP
         365             : !
         366             : ! !IROUTINE: HCO_GetSUNCOS
         367             : !
         368             : ! !DESCRIPTION: Subroutine HCO\_GetSUNCOS calculates the solar zenith angle
         369             : ! for the given date.
         370             : !\\
         371             : !\\
         372             : ! !INTERFACE:
         373             : !
         374           0 :   SUBROUTINE HCO_GetSUNCOS( HcoState, SUNCOS, DT, RC )
         375             : !
         376             : ! !USES
         377             : !
         378             :     USE HCO_STATE_MOD,   ONLY : HCO_STATE
         379             :     USE HCO_CLOCK_MOD,   ONLY : HcoClock_Get
         380             :     USE HCO_CLOCK_MOD,   ONLY : HcoClock_GetLocal
         381             : !
         382             : ! !INPUT/OUTPUT PARAMETERS:
         383             : !
         384             :     TYPE(HCO_State), POINTER        :: HcoState       ! HEMCO state object
         385             :     INTEGER,         INTENT(IN   )  :: DT             ! Time shift relative to current date [hrs]
         386             : !
         387             : ! !OUTPUT PARAMETERS:
         388             : !
         389             :     REAL(hp),        INTENT(  OUT)  :: SUNCOS(HcoState%NX,HcoState%NY)
         390             : !
         391             : ! !INPUT/OUTPUT PARAMETERS:
         392             : !
         393             :     INTEGER,         INTENT(INOUT)  :: RC             ! Return code
         394             : !
         395             : ! !REVISION HISTORY:
         396             : !  22 May 2015 - C. Keller   - Initial version, based on GEOS-Chem's dao_mod.F.
         397             : !  See https://github.com/geoschem/hemco for complete history
         398             : !EOP
         399             : !------------------------------------------------------------------------------
         400             : !BOC
         401             : !
         402             : ! !LOCAL VARIABLES:
         403             : !
         404             :     INTEGER              :: I, J,   DOY, HOUR
         405             :     LOGICAL              :: ERR
         406             :     REAL(hp)             :: YMID_R, S_YMID_R,  C_YMID_R
         407             :     REAL(hp)             :: R,      DEC
         408             :     REAL(hp)             :: S_DEC,  C_DEC
         409             :     REAL(hp)             :: SC,     LHR
         410             :     REAL(hp)             :: AHR
         411             :
         412             :     ! Coefficients for solar declination angle
         413             :     REAL(hp),  PARAMETER :: A0 = 0.006918e+0_hp
         414             :     REAL(hp),  PARAMETER :: A1 = 0.399912e+0_hp
         415             :     REAL(hp),  PARAMETER :: A2 = 0.006758e+0_hp
         416             :     REAL(hp),  PARAMETER :: A3 = 0.002697e+0_hp
         417             :     REAL(hp),  PARAMETER :: B1 = 0.070257e+0_hp
         418             :     REAL(hp),  PARAMETER :: B2 = 0.000907e+0_hp
         419             :     REAL(hp),  PARAMETER :: B3 = 0.000148e+0_hp
         420             :
         421             :     CHARACTER(LEN=255) :: LOC
         422             :
         423             :     !-------------------------------
         424             :     ! HCO_GetSUNCOS starts here!
         425             :     !-------------------------------
         426           0 :     LOC = 'HCO_GetSUNCOS (HCO_GEOTOOLS_MOD.F90)'
         427             :
         428             :     ! Get current time information
         429           0 :     CALL HcoClock_Get( HcoState%Clock, cDOY=DOY, cH=HOUR, RC=RC )
         430           0 :     IF ( RC /= HCO_SUCCESS ) THEN
         431           0 :         CALL HCO_ERROR( 'ERROR 0', RC, THISLOC=LOC )
         432           0 :         RETURN
         433             :     ENDIF
         434             :
         435             :     ! Add time adjustment
         436           0 :     HOUR = HOUR + DT
         437             :
         438             :     ! Make sure HOUR is within valid range (0-24)
         439           0 :     IF ( HOUR < 0 ) THEN
         440           0 :        HOUR = HOUR + 24
         441           0 :        DOY  = DOY  - 1
         442           0 :     ELSEIF ( HOUR > 23 ) THEN
         443           0 :        HOUR = HOUR - 24
         444           0 :        DOY  = DOY  + 1
         445             :     ENDIF
         446             :
         447             :     ! Make sure DOY is within valid range of 1 to 365
         448           0 :     DOY = MAX(MIN(DOY,365),1)
         449             :
         450             :     ! Path length of earth's orbit traversed since Jan 1 [radians]
         451           0 :     R        = ( 2e+0_hp * HcoState%Phys%PI / 365e+0_hp ) * DBLE( DOY - 1 )
         452             :
         453             :     ! Solar declination angle (low precision formula) [radians]
         454             :     DEC      = A0 - A1*COS(         R ) + B1*SIN(         R ) &
         455             :                   - A2*COS( 2e+0_hp*R ) + B2*SIN( 2e+0_hp*R ) &
         456           0 :                   - A3*COS( 3e+0_hp*R ) + B3*SIN( 3e+0_hp*R )
         457             :
         458             :     ! Pre-compute sin & cos of DEC outside of DO loops (for efficiency)
         459           0 :     S_DEC    = SIN( DEC )
         460           0 :     C_DEC    = COS( DEC )
         461             :
         462             :     ! Init
         463           0 :     ERR = .FALSE.
         464             :
         465             :     !=================================================================
         466             :     ! Compute cosine of solar zenith angle
         467             :     !=================================================================
         468             : !$OMP PARALLEL DO                                         &
         469             : !$OMP DEFAULT( SHARED )                                   &
         470             : !$OMP PRIVATE( I,      J,   YMID_R, S_YMID_R,  C_YMID_R ) &
         471             : !$OMP PRIVATE( LHR,    AHR, SC,     RC                  )
         472           0 :     DO J = 1, HcoState%NY
         473           0 :     DO I = 1, HcoState%NX
         474             :
         475             :          ! Latitude of grid box [radians]
         476           0 :          YMID_R     = HcoState%Grid%YMID%Val(I,J) * HcoState%Phys%PI_180
         477             :
         478             :          ! Pre-compute sin & cos of DEC outside of I loop (for efficiency)
         479           0 :          S_YMID_R   = SIN( YMID_R )
         480           0 :          C_YMID_R   = COS( YMID_R )
         481             :
         482             :          !==============================================================
         483             :          ! Compute cosine of SZA at the midpoint of the chem timestep
         484             :          ! Required for photolysis, chemistry, emissions, drydep
         485             :          !==============================================================
         486             :
         487             : !-----------------------------------------------------------------------------
         488             : ! Prior to 3/2/17:
         489             : ! Seb Eastham suggested to comment out the call to HcoClock_GetLocal.  If
         490             : ! the Voronoi timezones are used, this will compute the timezones on political
         491             : ! boundaries and not strictly on longitude.  This will cause funny results.
         492             : ! Replace this with a strict longitudinal local time. (bmy, 3/27/17)
         493             : !         ! Local time [hours] at box (I,J) at the midpt of the chem timestep
         494             : !         CALL HcoClock_GetLocal ( HcoState, I, J, cH=LHR, RC=RC )
         495             : !
         496             : !         IF ( RC /= HCO_SUCCESS ) THEN
         497             : !            ERR = .TRUE.
         498             : !            EXIT
         499             : !         ENDIF
         500             : !-----------------------------------------------------------------------------
         501             : ! Prior to 3/2/17:
         502             : ! Seb Eastham says that HOUR (in the new formula below) already contains DT.
         503             : ! so we need to comment this out and just use HOUR + LONGITUDE/15.
         504             : ! (bmy, 3/2/17)
         505             : !         ! Adjust for time shift
         506             : !         LHR = LHR + DT
         507             : !----------------------------------------------------------------------------
         508             :
         509             :          ! Compute local time as UTC + longitude/15 (bmy, 3/2/17)
         510           0 :          LHR = HOUR + ( HcoState%Grid%XMid%Val(I,J) / 15.0_hp )
         511             :
         512           0 :          IF ( LHR <   0.0_hp ) LHR = LHR + 24.0_hp
         513           0 :          IF ( LHR >= 24.0_hp ) LHR = LHR - 24.0_hp
         514             :
         515             :          ! Hour angle at box (I,J) [radians]
         516           0 :          AHR = ABS( LHR - 12.0_hp ) * 15.0_hp * HcoState%Phys%PI_180
         517             :
         518             :          ! Corresponding cosine( SZA ) at box (I,J) [unitless]
         519             :          SC = ( S_YMID_R * S_DEC              ) &
         520           0 :             + ( C_YMID_R * C_DEC * COS( AHR ) )
         521             :
         522             :          ! COS(SZA) at the current time
         523           0 :          SUNCOS(I,J) = SC
         524             :
         525             :       ENDDO
         526             :       ENDDO
         527             : !$OMP END PARALLEL DO
         528             :
         529             :     ! Check error status
         530             :     IF ( ERR ) THEN
         531             :        CALL HCO_ERROR ( &
         532             :          'Cannot calculate SZA', RC, &
         533             :           THISLOC='HCO_GetSUNCOS (hco_geotools_mod.F90)' )
         534             :        RETURN
         535             :     ENDIF
         536             :
         537             :     ! Leave w/ success
         538           0 :     RC = HCO_SUCCESS
         539             :
         540             :   END SUBROUTINE HCO_GetSUNCOS
         541             : !EOC
         542             : #if defined(ESMF_)
         543             : !------------------------------------------------------------------------------
         544             : !                   Harmonized Emissions Component (HEMCO)                    !
         545             : !------------------------------------------------------------------------------
         546             : !BOP
         547             : !
         548             : ! !IROUTINE: HCO_GetHorzIJIndex
         549             : !
         550             : ! !DESCRIPTION: Function HCO\_GetHorzIJIndex returns the grid box index for
         551             : !  the given longitude (deg E, -180...180), and latitude (deg N, -90...90).
         552             : !\\
         553             : !\\
         554             : ! !INTERFACE:
         555             : !
         556             :   SUBROUTINE HCO_GetHorzIJIndex( HcoState, N, Lon, Lat, idx, jdx, RC )
         557             : !
         558             : ! !USES
         559             : !
         560             : #include "MAPL_Generic.h"
         561             :     USE ESMF
         562             :     USE MAPLBase_Mod
         563             :     USE HCO_STATE_MOD,   ONLY : HCO_STATE
         564             : !
         565             : ! !INPUT PARAMETERS:
         566             : !
         567             :     TYPE(HCO_State), POINTER        :: HcoState       ! HEMCO state object
         568             :     INTEGER,         INTENT(IN   )  :: N
         569             :     REAL(hp),        INTENT(IN   )  :: Lon(N)
         570             :     REAL(hp),        INTENT(IN   )  :: Lat(N)
         571             : !
         572             : ! !INPUT/OUTPUT PARAMETERS:
         573             : !
         574             :     INTEGER,         INTENT(INOUT)  :: RC
         575             : !
         576             : ! !OUTPUT PARAMETERS:
         577             : !
         578             :     INTEGER,         INTENT(  OUT)  :: IDX(N), JDX(N)
         579             : !
         580             : ! !REVISION HISTORY:
         581             : !  04 Jun 2015 - C. Keller - Initial version
         582             : !  See https://github.com/geoschem/hemco for complete history
         583             : !EOP
         584             : !------------------------------------------------------------------------------
         585             : !BOC
         586             : !
         587             : ! !LOCAL VARIABLES:
         588             : !
         589             :     REAL             :: LonR(N), LatR(N)
         590             :     REAL, PARAMETER  :: radToDeg = 57.2957795
         591             :     TYPE(ESMF_Grid)  :: Grid
         592             :
         593             :     ! Defined Iam and STATUS
         594             :     __Iam__("HCO_GetHorzIJIndex (hco_geotools_mod.F90)")
         595             :
         596             :     !-------------------------------
         597             :     ! HCO_GetHorzIJIndex begins here
         598             :     !-------------------------------
         599             :
         600             :     ! Get grid
         601             :     ASSERT_(ASSOCIATED(HcoState%GridComp))
         602             :     CALL ESMF_GridCompGet( HcoState%GridComp, Grid=Grid, __RC__ )
         603             :
         604             :     ! Shadow variables
         605             :     LonR(:) = Lon / radToDeg
         606             :     LatR(:) = Lat / radToDeg
         607             :
         608             :     ! Get indeces
         609             :     CALL MAPL_GetHorzIJIndex( npts=N,   II=idx,   JJ=jdx,    &
         610             :                               lon=LonR, lat=LatR, Grid=Grid, &
         611             :                               __RC__)
         612             :
         613             : !!! old version of MAPL:
         614             : !    CALL MAPL_GetHorzIJIndex(N,idx,jdx,LonR,LatR,Grid=Grid,__RC__)
         615             :
         616             :     ! Return w/ success
         617             :     RC =  HCO_SUCCESS
         618             :
         619             :   END SUBROUTINE HCO_GetHorzIJIndex
         620             : !EOC
         621             : #else
         622             : !------------------------------------------------------------------------------
         623             : !                   Harmonized Emissions Component (HEMCO)                    !
         624             : !------------------------------------------------------------------------------
         625             : !BOP
         626             : !
         627             : ! !IROUTINE: HCO_GetHorzIJIndex
         628             : !
         629             : ! !DESCRIPTION: Function HCO\_GetHorzIJIndex returns the grid box index for
         630             : !  the given longitude (deg E, -180...180), and latitude (deg N, -90...90).
         631             : !\\
         632             : !\\
         633             : ! !INTERFACE:
         634             : !
         635           0 :   SUBROUTINE HCO_GetHorzIJIndex( HcoState, N, Lon, Lat, idx, jdx, RC )
         636             : !
         637             : ! !USES:
         638             : !
         639             :     USE HCO_STATE_MOD,   ONLY : HCO_STATE
         640             : !
         641             : ! !INPUT PARAMETERS:
         642             : !
         643             :     TYPE(HCO_State), POINTER        :: HcoState       ! HEMCO state object
         644             :     INTEGER,         INTENT(IN   )  :: N
         645             :     REAL(hp),        INTENT(IN   )  :: Lon(N)
         646             :     REAL(hp),        INTENT(IN   )  :: Lat(N)
         647             : !
         648             : ! !INPUT/OUTPUT PARAMETERS:
         649             : !
         650             :     INTEGER,         INTENT(INOUT)  :: RC
         651             : !
         652             : ! !OUTPUT PARAMETERS:
         653             : !
         654             :     INTEGER,         INTENT(  OUT)  :: IDX(N), JDX(N)
         655             : !
         656             : ! !REVISION HISTORY:
         657             : !  04 Jun 2015 - C. Keller - Initial version
         658             : !  See https://github.com/geoschem/hemco for complete history
         659             : !EOP
         660             : !------------------------------------------------------------------------------
         661             : !BOC
         662             : !
         663             : ! !LOCAL VARIABLES:
         664             : !
         665             :     INTEGER  :: I, J, L, FOUND
         666             :     REAL(hp) :: iLon1, iLon2
         667             :     REAL(hp) :: iLat1, iLat2
         668             :     REAL(hp) :: delta
         669             :
         670             :     !-------------------------------
         671             :     ! HCO_GetHorzIJIndex begins here
         672             :     !-------------------------------
         673             :
         674             :     ! Initialize
         675           0 :     IDX(:) = -1
         676           0 :     JDX(:) = -1
         677           0 :     FOUND  =  0
         678             :
         679             :     ! do for every grid box
         680           0 :     DO J = 1, HcoState%NY
         681           0 :     DO I = 1, HcoState%NX
         682             :
         683             :        ! Get grid edges for this box
         684             :
         685             :        ! Longitude edges
         686           0 :        IF ( ASSOCIATED(HcoState%Grid%XEDGE%Val) ) THEN
         687           0 :           iLon1 = HcoState%Grid%XEDGE%Val(I,  J)
         688           0 :           iLon2 = HcoState%Grid%XEDGE%Val(I+1,J)
         689             :
         690             :        ! Approximate from mid points
         691             :        ELSE
         692           0 :           iLon1 = HcoState%Grid%XMID%Val(I,J)
         693           0 :           IF ( I < HcoState%NX ) THEN
         694           0 :              delta = HcoState%Grid%XMID%Val(I+1,J)-iLon1
         695             :           ELSE
         696           0 :              delta = iLon1 - HcoState%Grid%XMID%Val(I-1,J)
         697             :           ENDIF
         698           0 :           iLon2 = iLon1 + (delta / 2.0_hp)
         699           0 :           iLon1 = iLon1 - (delta / 2.0_hp)
         700             :        ENDIF
         701             :
         702             :        ! Latitude edges
         703           0 :        IF ( ASSOCIATED(HcoState%Grid%YEDGE%Val) ) THEN
         704           0 :           iLat1 = HcoState%Grid%YEDGE%Val(I,J)
         705           0 :           iLat2 = HcoState%Grid%YEDGE%Val(I,J+1)
         706             :
         707             :        ! Approximate from mid points
         708             :        ELSE
         709           0 :           iLat1 = HcoState%Grid%YMID%Val(I,J)
         710           0 :           IF ( J < HcoState%NY ) THEN
         711           0 :              delta = HcoState%Grid%YMID%Val(I,J+1)-iLat1
         712             :           ELSE
         713           0 :              delta = iLat1 - HcoState%Grid%YMID%Val(I,J-1)
         714             :           ENDIF
         715           0 :           iLat2 = iLat1 + (delta / 2.0_hp)
         716           0 :           iLat1 = iLat1 - (delta / 2.0_hp)
         717             :        ENDIF
         718             :
         719             :        ! Check if it's within this box
         720           0 :        DO L = 1, N
         721           0 :           IF ( IDX(L) > 0 ) CYCLE
         722             :
         723           0 :           IF ( Lon(L) >= HcoState%Grid%XEDGE%Val(I,  J  ) .AND. &
         724           0 :                Lon(L) <= HcoState%Grid%XEDGE%Val(I+1,J  ) .AND. &
         725           0 :                Lat(L) >= HcoState%Grid%YEDGE%Val(I  ,J  ) .AND. &
         726           0 :                Lat(L) <= HcoState%Grid%YEDGE%Val(I  ,J+1)        ) THEN
         727           0 :              IDX(L) = I
         728           0 :              JDX(L) = J
         729           0 :              FOUND  = FOUND + 1
         730           0 :              IF ( FOUND == N ) EXIT
         731             :           ENDIF
         732             :        ENDDO
         733             :     ENDDO
         734             :     ENDDO
         735             :
         736             :     ! Return w/ success
         737           0 :     RC = HCO_SUCCESS
         738             :
         739           0 :   END SUBROUTINE HCO_GetHorzIJIndex
         740             : !EOC
         741             : #endif
         742             : !------------------------------------------------------------------------------
         743             : !                   Harmonized Emissions Component (HEMCO)                    !
         744             : !------------------------------------------------------------------------------
         745             : !BOP
         746             : !
         747             : ! !IROUTINE: HCO_CalcVertGrid
         748             : !
         749             : ! !DESCRIPTION: Function HCO\_CalcVertGrid calculates the vertical grid
         750             : !  quantities surface pressure PSFC [Pa], surface geopotential height ZSFC
         751             : !  [m], grid box height BXHEIGHT [m], and pressure edges PEDGE [Pa]. Any of
         752             : !  these fields can be passed explicitly to the routine, in which case these
         753             : !  fields are being used. If not passed through the routine (i.e. if the
         754             : !  corresponding input argument pointer is nullified), the field is searched
         755             : !  in the HEMCO configuration file. If not found in the configuration file,
         756             : !  the field is approximated from other quantities (if possible). For example,
         757             : !  if surface pressures are provided (either passed as argument or in the
         758             : !  HEMCO configuration file as field PSFC), pressure edges are calculated
         759             : !  from PSFC and the vertical grid coordinates (Ap and Bp for a hybrid sigma
         760             : !  coordinate system). The temperature field TK [K] is needed to approximate
         761             : !  box heights and/or geopotential height (via the hydrostatic equation).
         762             : !\\
         763             : !\\
         764             : ! !INTERFACE:
         765             : !
         766           0 :   SUBROUTINE HCO_CalcVertGrid ( HcoState, PSFC, ZSFC, TK, BXHEIGHT, PEDGE, RC )
         767             : !
         768             : ! !USES
         769             : !
         770             :     USE HCO_Arr_Mod,      ONLY : HCO_ArrAssert, HCO_ArrCleanup
         771             :     USE HCO_STATE_MOD,    ONLY : HCO_STATE
         772             :     USE HCO_CALC_MOD,     ONLY : HCO_EvalFld
         773             : !
         774             : ! !INPUT PARAMETERS:
         775             : !
         776             :     TYPE(HCO_State), POINTER        :: HcoState        ! HEMCO state object
         777             :     REAL(hp),        POINTER        :: PSFC(:,:)       ! surface pressure (Pa)
         778             :     REAL(hp),        POINTER        :: ZSFC(:,:)       ! surface geopotential height (m)
         779             :     REAL(hp),        POINTER        :: TK  (:,:,:)     ! air temperature (K)
         780             :     REAL(hp),        POINTER        :: BXHEIGHT(:,:,:) ! grid box height (m)
         781             :     REAL(hp),        POINTER        :: PEDGE(:,:,:)    ! pressure edges (Pa)
         782             : !
         783             : ! !INPUT/OUTPUT PARAMETERS:
         784             : !
         785             :     INTEGER,         INTENT(INOUT)  :: RC
         786             : !
         787             : ! !REVISION HISTORY:
         788             : !  28 Sep 2015 - C. Keller - Initial version
         789             : !  See https://github.com/geoschem/hemco for complete history
         790             : !EOP
         791             : !------------------------------------------------------------------------------
         792             : !BOC
         793             : !
         794             : ! !LOCAL VARIABLES:
         795             : !
         796             :     INTEGER                       :: NX, NY, NZ
         797             :     INTEGER                       :: I,  J,  L
         798             :     LOGICAL                       :: Verb
         799             :     LOGICAL                       :: FoundPSFC
         800             :     LOGICAL                       :: FoundZSFC
         801             :     LOGICAL                       :: FoundTK
         802             :     LOGICAL                       :: FoundPEDGE
         803             :     LOGICAL                       :: FoundBXHEIGHT
         804             :     LOGICAL                       :: ERRBX, ERRZSFC
         805             :     REAL(hp)                      :: P1, P2
         806           0 :     REAL(hp), ALLOCATABLE, TARGET :: TmpTK(:,:,:)
         807           0 :     REAL(hp), POINTER             :: ThisTK(:,:,:)
         808             :     CHARACTER(LEN=255)            :: MSG
         809             :     CHARACTER(LEN=255)            :: LOC = 'HCO_CalcVertGrid (hco_geotools_mod.F90)'
         810             :
         811             :     LOGICAL, SAVE                 :: FIRST         = .TRUE.
         812             :     LOGICAL, SAVE                 :: EVAL_PSFC     = .TRUE.
         813             :     LOGICAL, SAVE                 :: EVAL_ZSFC     = .TRUE.
         814             :     LOGICAL, SAVE                 :: EVAL_TK       = .TRUE.
         815             :     LOGICAL, SAVE                 :: EVAL_PEDGE    = .TRUE.
         816             :     LOGICAL, SAVE                 :: EVAL_BXHEIGHT = .TRUE.
         817             :
         818             :     !-------------------------------
         819             :     ! HCO_CalcVertGrid begins here
         820             :     !-------------------------------
         821             :
         822             :     ! Init
         823           0 :     Verb          = .FALSE.
         824           0 :     FoundPSFC     = .FALSE.
         825           0 :     FoundZSFC     = .FALSE.
         826           0 :     FoundTK       = .FALSE.
         827           0 :     FoundPEDGE    = .FALSE.
         828           0 :     FoundBXHEIGHT = .FALSE.
         829           0 :     ThisTK        => NULL()
         830             :
         831             :     ! Verbose statements
         832           0 :     IF ( HcoState%amIRoot .AND. FIRST .AND. &
         833             :          HCO_IsVerb(HcoState%Config%Err) ) THEN
         834             :        Verb = .TRUE.
         835             :     ENDIF
         836             :     IF ( Verb ) THEN
         837           0 :        MSG = 'Details about vertical grid calculations (only shown on first time step):'
         838           0 :        CALL HCO_MSG(HcoState%Config%Err,MSG,SEP1='-')
         839           0 :        MSG = '1. Input data availability: '
         840           0 :        CALL HCO_MSG(HcoState%Config%Err,MSG,SEP1=' ')
         841             :     ENDIF
         842             :
         843             :     ! ------------------------------------------------------------------
         844             :     ! TK
         845             :     ! ------------------------------------------------------------------
         846             :
         847             :     ! If associated, make sure that array size is correct
         848             :     ! and pass to HEMCO surface pressure field
         849           0 :     IF ( ASSOCIATED(TK) ) THEN
         850           0 :        NX = SIZE(TK,1)
         851           0 :        NY = SIZE(TK,2)
         852           0 :        NZ = SIZE(TK,3)
         853           0 :        IF ( NX /= HcoState%NX .OR. NY /= HcoState%NY .OR. &
         854             :             NZ /= HcoState%NZ ) THEN
         855           0 :           WRITE(MSG,*) 'Wrong TK array size: ', NX, NY, NZ, &
         856           0 :                        '; should be: ', HcoState%NX, HcoState%NY, HcoState%NZ
         857           0 :           CALL HCO_ERROR( MSG, RC, THISLOC=LOC )
         858           0 :           RETURN
         859             :        ENDIF
         860             :
         861             :        ! TK is not a field in grid, so don't pass
         862           0 :        ThisTK  => TK
         863           0 :        FoundTK = .TRUE.
         864             :
         865             :        ! Verbose
         866           0 :        IF ( Verb ) THEN
         867           0 :           WRITE(MSG,*) ' - Temperature field TK obtained from model interface (min,max): ', MINVAL(ThisTK), MAXVAL(ThisTK)
         868           0 :           CALL HCO_MSG(HcoState%Config%Err,MSG)
         869             :        ENDIF
         870             :
         871           0 :     ELSEIF ( EVAL_TK ) THEN
         872           0 :        ALLOCATE(TmpTK(HcoState%NX,HcoState%NY,HcoState%NZ))
         873           0 :        CALL HCO_EvalFld( HcoState, 'TK', TmpTK, RC, FOUND=FoundTK )
         874           0 :        IF ( RC /= HCO_SUCCESS ) THEN
         875           0 :            CALL HCO_ERROR( 'ERROR 1', RC, THISLOC=LOC )
         876           0 :            RETURN
         877             :        ENDIF
         878             :
         879             :        ! TK is sometimes listed as TMPU, so look for that too (bmy, 3/5/21)
         880           0 :        IF ( .not. FoundTK ) THEN
         881           0 :           CALL HCO_EvalFld( HcoState, 'TMPU', TmpTK, RC, FOUND=FoundTK )
         882           0 :           IF ( RC /= HCO_SUCCESS ) THEN
         883           0 :               CALL HCO_ERROR( 'ERROR 2', RC, THISLOC=LOC )
         884           0 :               RETURN
         885             :           ENDIF
         886             :        ENDIF
         887             :
         888           0 :        EVAL_TK = FoundTK
         889           0 :        IF ( FoundTK ) ThisTK => TmpTk
         890             :
         891             :        ! Verbose
         892           0 :        IF ( Verb ) THEN
         893           0 :           IF ( FoundTK ) THEN
         894           0 :              WRITE(MSG,*) ' - Temperature field TK [K] obtained from configuration file (min,max): ', MINVAL(ThisTK), MAXVAL(ThisTK)
         895           0 :              CALL HCO_MSG(HcoState%Config%Err,MSG)
         896             :           ELSE
         897           0 :              WRITE(MSG,*) ' - No temperature field TK found - some vertical grid calculations may not be performed...'
         898           0 :              CALL HCO_MSG(HcoState%Config%Err,MSG)
         899             :           ENDIF
         900             :        ENDIF
         901             :     ENDIF
         902             :
         903             :     ! ------------------------------------------------------------------
         904             :     ! PSFC
         905             :     ! ------------------------------------------------------------------
         906           0 :     CALL HCO_ArrAssert( HcoState%Grid%PSFC, HcoState%NX, HcoState%NY, RC )
         907           0 :     IF ( RC /= HCO_SUCCESS ) THEN
         908           0 :         CALL HCO_ERROR( 'ERROR 3', RC, THISLOC=LOC )
         909           0 :         RETURN
         910             :     ENDIF
         911             :
         912             :     ! If associated, make sure that array size is correct
         913             :     ! and pass to HEMCO surface pressure field
         914           0 :     IF ( ASSOCIATED(PSFC) ) THEN
         915           0 :        NX = SIZE(PSFC,1)
         916           0 :        NY = SIZE(PSFC,2)
         917           0 :        IF ( NX /= HcoState%NX .OR. NY /= HcoState%NY ) THEN
         918           0 :           WRITE(MSG,*) 'Wrong PSFC array size: ', NX, NY, &
         919           0 :                        '; should be: ', HcoState%NX, HcoState%NY
         920           0 :           CALL HCO_ERROR( MSG, RC, THISLOC=LOC )
         921           0 :           RETURN
         922             :        ENDIF
         923             :
         924             :        ! Pass to HcoState array
         925           0 :        HcoState%Grid%PSFC%Val = PSFC
         926           0 :        FoundPSFC              = .TRUE.
         927             :
         928             :        ! Verbose
         929           0 :        IF ( Verb ) THEN
         930           0 :           WRITE(MSG,*) ' - Surface pressure PSFC [Pa] obtained from model interface (min, max): ', MINVAL(HcoState%Grid%PSFC%Val), MAXVAL(HcoState%Grid%PSFC%VAL)
         931           0 :           CALL HCO_MSG(HcoState%Config%Err,MSG)
         932             :        ENDIF
         933             :
         934             :     ! Otherwise, try to read from HEMCO configuration file
         935           0 :     ELSEIF ( EVAL_PSFC ) THEN
         936             :        CALL HCO_EvalFld( HcoState, 'PSFC', HcoState%Grid%PSFC%Val, RC, &
         937           0 :             FOUND=FoundPSFC )
         938           0 :        IF ( RC /= HCO_SUCCESS ) THEN
         939           0 :            CALL HCO_ERROR( 'ERROR 4', RC, THISLOC=LOC )
         940           0 :            RETURN
         941             :        ENDIF
         942             :
         943             :        ! PSFC is sometimes listed as PS, so look for that too (bmy, 3/4/21)
         944           0 :        IF ( .not. FoundPSFC ) THEN
         945             :           CALL HCO_EvalFld( HcoState, 'PS', HcoState%Grid%PSFC%Val, RC, &
         946           0 :                FOUND=FoundPSFC )
         947           0 :           IF ( RC /= HCO_SUCCESS ) THEN
         948           0 :               CALL HCO_ERROR( 'ERROR 5', RC, THISLOC=LOC )
         949           0 :               RETURN
         950             :           ENDIF
         951             :        ENDIF
         952             :
         953           0 :        EVAL_PSFC = FoundPSFC
         954             :
         955             :        ! Verbose
         956           0 :        IF ( Verb ) THEN
         957           0 :           IF ( FoundPSFC ) THEN
         958           0 :              WRITE(MSG,*) ' - Surface pressure PSFC [Pa] obtained from configuration file (min, max): ', MINVAL(HcoState%Grid%PSFC%Val), MAXVAL(HcoState%Grid%PSFC%VAL)
         959           0 :              CALL HCO_MSG(HcoState%Config%Err,MSG)
         960             :           ELSE
         961           0 :              MSG = ' - Surface pressure PSFC not found. Will attempt to calculate it.'
         962           0 :              CALL HCO_MSG(HcoState%Config%Err,MSG)
         963             :           ENDIF
         964             :        ENDIF
         965             :     ENDIF
         966             :
         967             :     ! ------------------------------------------------------------------
         968             :     ! ZSFC
         969             :     ! ------------------------------------------------------------------
         970           0 :     CALL HCO_ArrAssert( HcoState%Grid%ZSFC, HcoState%NX, HcoState%NY, RC )
         971           0 :     IF ( RC /= HCO_SUCCESS ) THEN
         972           0 :         CALL HCO_ERROR( 'ERROR 6', RC, THISLOC=LOC )
         973           0 :         RETURN
         974             :     ENDIF
         975             :
         976             :     ! If associated, make sure that array size is correct
         977             :     ! and pass to HEMCO surface pressure field
         978           0 :     IF ( ASSOCIATED(ZSFC) ) THEN
         979           0 :        NX = SIZE(ZSFC,1)
         980           0 :        NY = SIZE(ZSFC,2)
         981           0 :        IF ( NX /= HcoState%NX .OR. NY /= HcoState%NY ) THEN
         982           0 :           WRITE(MSG,*) 'Wrong ZSFC array size: ', NX, NY, &
         983           0 :                        '; should be: ', HcoState%NX, HcoState%NY
         984           0 :           CALL HCO_ERROR( MSG, RC, THISLOC=LOC )
         985           0 :           RETURN
         986             :        ENDIF
         987             :
         988             :        ! Pass to HcoState array
         989           0 :        HcoState%Grid%ZSFC%Val = ZSFC
         990           0 :        FoundZSFC              = .TRUE.
         991             :
         992             :        ! Verbose
         993           0 :        IF ( Verb ) THEN
         994           0 :           WRITE(MSG,*) ' - Surface geopotential height ZSFC [m] obtained from model interface (min, max): ', MINVAL(HcoState%Grid%ZSFC%Val), MAXVAL(HcoState%Grid%ZSFC%VAL)
         995           0 :           CALL HCO_MSG(HcoState%Config%Err,MSG)
         996             :        ENDIF
         997             :
         998             :     ! Otherwise, try to read from HEMCO configuration file
         999           0 :     ELSEIF ( EVAL_ZSFC ) THEN
        1000           0 :        CALL HCO_EvalFld ( HcoState, 'ZSFC', HcoState%Grid%ZSFC%Val, RC, FOUND=FoundZSFC )
        1001           0 :        IF ( RC /= HCO_SUCCESS ) THEN
        1002           0 :            CALL HCO_ERROR( 'ERROR 7', RC, THISLOC=LOC )
        1003           0 :            RETURN
        1004             :        ENDIF
        1005           0 :        EVAL_ZSFC = FoundZSFC
        1006             :
        1007             :        ! Verbose
        1008           0 :        IF ( Verb ) THEN
        1009           0 :           IF ( FoundZSFC ) THEN
        1010           0 :              WRITE(MSG,*) ' - Surface geopotential height ZSFC [m] obtained from configuration file (min, max): ', MINVAL(HcoState%Grid%ZSFC%Val), MAXVAL(HcoState%Grid%ZSFC%VAL)
        1011           0 :              CALL HCO_MSG(HcoState%Config%Err,MSG)
        1012             :           ELSE
        1013           0 :              MSG = ' - Surface geopotential height ZSFC not found. Will attempt to calculate it.'
        1014           0 :              CALL HCO_MSG(HcoState%Config%Err,MSG)
        1015             :           ENDIF
        1016             :        ENDIF
        1017             :     ENDIF
        1018             :
        1019             :     ! ------------------------------------------------------------------
        1020             :     ! PEDGE
        1021             :     ! ------------------------------------------------------------------
        1022             :     CALL HCO_ArrAssert( HcoState%Grid%PEDGE, HcoState%NX, &
        1023           0 :                         HcoState%NY,         HcoState%NZ+1, RC )
        1024           0 :     IF ( RC /= HCO_SUCCESS ) THEN
        1025           0 :         CALL HCO_ERROR( 'ERROR 8', RC, THISLOC=LOC )
        1026           0 :         RETURN
        1027             :     ENDIF
        1028             :
        1029           0 :     IF ( ASSOCIATED( PEDGE ) ) THEN
        1030             :
        1031           0 :        NX = SIZE(PEDGE,1)
        1032           0 :        NY = SIZE(PEDGE,2)
        1033           0 :        NZ = SIZE(PEDGE,3)
        1034           0 :        IF ( NX /= HcoState%NX .OR. NY /= HcoState%NY .OR. &
        1035             :             NZ /= (HcoState%NZ + 1) ) THEN
        1036           0 :           WRITE(MSG,*) 'Wrong PEDGE array size: ', NX, NY, NZ, &
        1037           0 :                        '; should be: ', HcoState%NX, HcoState%NY, HcoState%NZ+1
        1038           0 :           CALL HCO_ERROR( MSG, RC, THISLOC=LOC )
        1039           0 :           RETURN
        1040             :        ENDIF
        1041             :
        1042           0 :        HcoState%Grid%PEDGE%Val = PEDGE
        1043           0 :        FoundPEDGE = .TRUE.
        1044             :
        1045             :        ! Verbose
        1046           0 :        IF ( Verb ) THEN
        1047           0 :           WRITE(MSG,*) ' - Pressure edges PEDGE obtained from model interface (min, max): ', MINVAL(HcoState%Grid%PEDGE%Val), MAXVAL(HcoState%Grid%PEDGE%VAL)
        1048           0 :           CALL HCO_MSG(HcoState%Config%Err,MSG)
        1049             :        ENDIF
        1050             :
        1051           0 :     ELSEIF ( EVAL_PEDGE ) THEN
        1052             :        CALL HCO_EvalFld ( HcoState, 'PEDGE', &
        1053           0 :                           HcoState%Grid%PEDGE%Val, RC, FOUND=FoundPEDGE )
        1054           0 :        IF ( RC /= HCO_SUCCESS ) THEN
        1055           0 :            CALL HCO_ERROR( 'ERROR 9', RC, THISLOC=LOC )
        1056           0 :            RETURN
        1057             :        ENDIF
        1058           0 :        EVAL_PEDGE = FoundPEDGE
        1059             :
        1060             :        ! Verbose
        1061           0 :        IF ( Verb ) THEN
        1062           0 :           IF ( FoundPEDGE ) THEN
        1063           0 :              WRITE(MSG,*) ' - Pressure edges PEDGE obtained from configuration file (min, max): ', MINVAL(HcoState%Grid%PEDGE%Val), MAXVAL(HcoState%Grid%PEDGE%VAL)
        1064           0 :              CALL HCO_MSG(HcoState%Config%Err,MSG)
        1065             :           ELSE
        1066           0 :              MSG = ' - Pressure edges PEDGE not found. Will attempt to calculate it.'
        1067           0 :              CALL HCO_MSG(HcoState%Config%Err,MSG)
        1068             :           ENDIF
        1069             :        ENDIF
        1070             :     ENDIF
        1071             :
        1072             :     ! ------------------------------------------------------------------
        1073             :     ! BXHEIGHT
        1074             :     ! ------------------------------------------------------------------
        1075             :     CALL HCO_ArrAssert( HcoState%Grid%BXHEIGHT_M, HcoState%NX, &
        1076           0 :                         HcoState%NY,              HcoState%NZ, RC )
        1077           0 :     IF ( RC /= HCO_SUCCESS ) THEN
        1078           0 :         CALL HCO_ERROR( 'ERROR 10', RC, THISLOC=LOC )
        1079           0 :         RETURN
        1080             :     ENDIF
        1081             :
        1082           0 :     IF ( ASSOCIATED( BXHEIGHT ) ) THEN
        1083             :
        1084           0 :        NX = SIZE(BXHEIGHT,1)
        1085           0 :        NY = SIZE(BXHEIGHT,2)
        1086           0 :        NZ = SIZE(BXHEIGHT,3)
        1087           0 :        IF ( NX /= HcoState%NX .OR. NY /= HcoState%NY .OR. &
        1088             :             NZ /= HcoState%NZ ) THEN
        1089           0 :           WRITE(MSG,*) 'Wrong BXHEIGHT array size: ', NX, NY, NZ, &
        1090           0 :                        '; should be: ', HcoState%NX, HcoState%NY, HcoState%NZ
        1091           0 :           CALL HCO_ERROR( MSG, RC, THISLOC=LOC )
        1092           0 :           RETURN
        1093             :        ENDIF
        1094             :
        1095           0 :        HcoState%Grid%BXHEIGHT_M%Val = BXHEIGHT
        1096           0 :        FoundBXHEIGHT = .TRUE.
        1097             :
        1098             :        ! Verbose
        1099           0 :        IF ( Verb ) THEN
        1100           0 :           WRITE(MSG,*) ' - Boxheights BXHEIGHT_M obtained from model interface (min, max): ', MINVAL(HcoState%Grid%BXHEIGHT_M%Val), MAXVAL(HcoState%Grid%BXHEIGHT_M%VAL)
        1101           0 :           CALL HCO_MSG(HcoState%Config%Err,MSG)
        1102             :        ENDIF
        1103             :
        1104             :     ! Otherwise, try to read from HEMCO configuration file
        1105           0 :     ELSEIF ( EVAL_BXHEIGHT ) THEN
        1106             :        CALL HCO_EvalFld ( HcoState, 'BXHEIGHT_M', &
        1107           0 :                           HcoState%Grid%BXHEIGHT_M%Val, RC, FOUND=FoundBXHEIGHT )
        1108           0 :        IF ( RC /= HCO_SUCCESS ) THEN
        1109           0 :            CALL HCO_ERROR( 'ERROR 11', RC, THISLOC=LOC )
        1110           0 :            RETURN
        1111             :        ENDIF
        1112           0 :        EVAL_BXHEIGHT = FoundBXHEIGHT
        1113             :
        1114             :        ! Verbose
        1115           0 :        IF ( Verb ) THEN
        1116           0 :           IF ( FoundBXHEIGHT ) THEN
        1117           0 :              WRITE(MSG,*) ' - Boxheights BXHEIGHT_M obtained from configuration file (min, max): ', MINVAL(HcoState%Grid%BXHEIGHT_M%Val), MAXVAL(HcoState%Grid%BXHEIGHT_M%VAL)
        1118           0 :              CALL HCO_MSG(HcoState%Config%Err,MSG)
        1119             :           ELSE
        1120           0 :              MSG = ' - Boxheights BXHEIGHT_M not found. Will attempt to calculate it.'
        1121           0 :              CALL HCO_MSG(HcoState%Config%Err,MSG)
        1122             :           ENDIF
        1123             :        ENDIF
        1124             :     ENDIF
        1125             :
        1126             :     ! If Box height isn't available, free its memory
        1127             :     ! NOTE: Using NULL instead of deallocate here causes a memory leak!
        1128           0 :     IF ( .NOT. EVAL_BXHEIGHT .OR. .NOT. FoundBXHEIGHT ) THEN
        1129           0 :        CALL HCO_ArrCleanup( HcoState%Grid%BXHEIGHT_M, DeepClean=.TRUE. )
        1130             :     ENDIF
        1131             :
        1132             :     ! ------------------------------------------------------------------
        1133             :     ! Calculate various quantities: the goal is to have the following
        1134             :     ! quantities defined: ZSFC, PSFC, PEDGE, BXHEIGHT_M.
        1135             :     ! - If PEDGE is not yet defined, it is calculated from surface
        1136             :     !   pressure (PSFC).
        1137             :     ! - If BXHEIGHT is not yet defined, it is calculated from PEDGE
        1138             :     !   and TK.
        1139             :     ! - If PSFC is not yet defined, it is set to the first level of
        1140             :     !   PEDGE (if defined) or uniformly set to 101325 Pa.
        1141             :     ! - If ZSFC is not yet defined, it is calculated from PSFC and TK.
        1142             :     ! - If TK is not defined, no attempt is made to initialize it to
        1143             :     !   a useful quantity. Calculations that require TK are omitted.
        1144             :     ! ------------------------------------------------------------------
        1145             :
        1146             :     ! Verbose
        1147           0 :     IF ( Verb ) THEN
        1148           0 :        MSG = '2. Grid calculations: '
        1149           0 :        CALL HCO_MSG(HcoState%Config%Err,MSG,SEP1=' ')
        1150             :     ENDIF
        1151             :
        1152             :     ! Set PSFC
        1153           0 :     IF ( .NOT. FoundPSFC ) THEN
        1154           0 :        IF ( FoundPEDGE ) THEN
        1155           0 :           HcoState%Grid%PSFC%Val(:,:) = HcoState%Grid%PEDGE%Val(:,:,1)
        1156             :
        1157             :           ! Verbose
        1158           0 :           IF ( Verb ) THEN
        1159           0 :              MSG = ' - Surface pressure set to surface pressure edge.'
        1160           0 :              CALL HCO_MSG(HcoState%Config%Err,MSG)
        1161             :           ENDIF
        1162             :        ELSE
        1163           0 :           HcoState%Grid%PSFC%Val(:,:) = 101325.0_hp
        1164           0 :           IF ( FIRST .AND. HcoState%amIRoot ) THEN
        1165             :              MSG = 'Surface pressure PSFC uniformly set to 101325 Pa! ' // &
        1166             :                    'This may affect the accuracy of vertical grid '     // &
        1167             :                    'quantities. It is recommended you provide PSFC via '// &
        1168           0 :                    'the model-HEMCO interface or the HEMCO configuration file!'
        1169           0 :              CALL HCO_WARNING( HcoState%Config%Err,MSG, RC, THISLOC=LOC )
        1170             :           ENDIF
        1171             :
        1172             :           ! Verbose
        1173           0 :           IF ( Verb ) THEN
        1174           0 :              MSG = ' - Surface pressure uniformly set to 101325.0 Pa.'
        1175           0 :              CALL HCO_MSG(HcoState%Config%Err,MSG)
        1176             :           ENDIF
        1177             :        ENDIF
        1178           0 :        FoundPSFC = .TRUE.
        1179             :     ENDIF
        1180             :
        1181             :     ! Set PEDGE
        1182           0 :     IF ( .NOT. FoundPEDGE ) THEN
        1183             :        !$OMP PARALLEL DO        &
        1184             :        !$OMP DEFAULT( SHARED  ) &
        1185             :        !$OMP PRIVATE( I, J, L )
        1186           0 :        DO L = 1, HcoState%NZ+1
        1187           0 :        DO J = 1, HcoState%NY
        1188           0 :        DO I = 1, HcoState%NX
        1189           0 :           HcoState%Grid%PEDGE%Val(I,J,L) &
        1190           0 :            = HcoState%Grid%zGrid%AP(L)   &
        1191           0 :            + ( HcoState%Grid%zGrid%BP(L) &
        1192           0 :              * HcoState%Grid%PSFC%Val(I,J) )
        1193             :        ENDDO
        1194             :        ENDDO
        1195             :        ENDDO
        1196             :        !$OMP END PARALLEL DO
        1197           0 :        FoundPEDGE = .TRUE.
        1198             :
        1199             :        ! Verbose
        1200           0 :        IF ( Verb ) THEN
        1201           0 :           WRITE(MSG,*) ' - PEDGE calculated from PSFC, Ap, and Bp (min, max): ', &
        1202           0 :              MINVAL(HcoState%Grid%PEDGE%Val), MAXVAL(HcoState%Grid%PEDGE%Val)
        1203           0 :           CALL HCO_MSG(HcoState%Config%Err,MSG)
        1204             :        ENDIF
        1205             :     ENDIF
        1206             :
        1207             :     ! Set surface height and/or grid box height
        1208           0 :     IF ( .NOT. FoundZSFC .OR. .NOT. FoundBXHEIGHT ) THEN
        1209           0 :        IF ( FoundTK .AND. FoundPEDGE ) THEN
        1210             :
        1211             :           ! Initialize error flags
        1212           0 :           ERRZSFC = .FALSE.
        1213           0 :           ERRBX   = .FALSE.
        1214             :
        1215             :           ! Make sure box height is initialized if it will be calculated
        1216             :           CALL HCO_ArrAssert( HcoState%Grid%BXHEIGHT_M, HcoState%NX, &
        1217           0 :                               HcoState%NY,              HcoState%NZ, RC )
        1218           0 :           IF ( RC /= HCO_SUCCESS ) THEN
        1219           0 :               CALL HCO_ERROR( 'ERROR 12', RC, THISLOC=LOC )
        1220           0 :               RETURN
        1221             :           ENDIF
        1222             :
        1223             :           !$OMP PARALLEL DO                &
        1224             :           !$OMP DEFAULT( SHARED          ) &
        1225             :           !$OMP PRIVATE( I, J, L, P1, P2 )
        1226           0 :           DO L = 1, HcoState%NZ
        1227           0 :           DO J = 1, HcoState%NY
        1228           0 :           DO I = 1, HcoState%NX
        1229             :
        1230             :              ! BOXHEIGHT (hydrostatic equation)
        1231           0 :              IF ( .NOT. FoundBXHEIGHT ) THEN
        1232           0 :                 P1 = HcoState%Grid%PEDGE%Val(I,J,L)
        1233           0 :                 P2 = HcoState%Grid%PEDGE%Val(I,J,L+1)
        1234           0 :                 IF ( P2 == 0.0_hp ) THEN
        1235             :                    ERRBX = .TRUE.
        1236             :                 ELSE
        1237           0 :                    HcoState%Grid%BXHEIGHT_M%Val(I,J,L) = HcoState%Phys%Rdg0 &
        1238           0 :                                                        * ThisTK(I,J,1)      &
        1239           0 :                                                        * LOG( P1 / P2 )
        1240             :                 ENDIF
        1241             :              ENDIF
        1242             :
        1243             :              ! ZSFC
        1244           0 :              IF ( L == 1 .AND. .NOT. FoundZSFC ) THEN
        1245           0 :                 P1 = 101325.0_hp
        1246           0 :                 P2 = HcoState%Grid%PEDGE%Val(I,J,1)
        1247           0 :                 IF ( P2 == 0.0_hp ) THEN
        1248             :                    ERRZSFC = .TRUE.
        1249             :                 ELSE
        1250           0 :                    HcoState%Grid%ZSFC%Val(I,J) = HcoState%Phys%Rdg0 &
        1251           0 :                                                * ThisTK(I,J,1)      &
        1252           0 :                                                * LOG( P1 / P2 )
        1253             :                 ENDIF
        1254             :              ENDIF
        1255             :           ENDDO
        1256             :           ENDDO
        1257             :           ENDDO
        1258             :           !$OMP END PARALLEL DO
        1259             :
        1260           0 :           IF ( ERRZSFC ) THEN
        1261             :              MSG = 'Cannot calculate surface geopotential heights - at least one ' // &
        1262             :                    'surface pressure value is zero! You can either provide an '    // &
        1263             :                    'updated pressure edge field (PEDGE) or add a field with the '  // &
        1264           0 :                    'surface geopotential height to your configuration file (ZSFC)'
        1265           0 :              CALL HCO_ERROR( MSG, RC, THISLOC=LOC )
        1266           0 :              RETURN
        1267             :           ELSE
        1268           0 :              FoundZSFC = .TRUE.
        1269             :
        1270             :              ! Verbose
        1271           0 :              IF ( Verb ) THEN
        1272           0 :                 WRITE(MSG,*) ' - ZSFC calculated from PSFC and T (min, max): ', &
        1273           0 :                    MINVAL(HcoState%Grid%ZSFC%Val), MAXVAL(HcoState%Grid%ZSFC%Val)
        1274           0 :                 CALL HCO_MSG(HcoState%Config%Err,MSG)
        1275             :              ENDIF
        1276             :           ENDIF
        1277             :
        1278           0 :           IF ( ERRBX ) THEN
        1279             :              MSG = 'Cannot calculate grid box heights - at least one ' // &
        1280             :                    'pressure value is zero! You can either provide an '    // &
        1281             :                    'updated pressure edge field (PEDGE) or add a field with the '  // &
        1282           0 :                    'grid box heights to your configuration file (BOXHEIGHT_M)'
        1283           0 :              CALL HCO_ERROR( MSG, RC, THISLOC=LOC )
        1284           0 :              RETURN
        1285             :           ELSE
        1286           0 :              FoundZSFC = .TRUE.
        1287             :
        1288             :              ! Verbose
        1289           0 :              IF ( Verb ) THEN
        1290           0 :                 WRITE(MSG,*) ' - Boxheights calculated from PEDGE and T (min, max): ', &
        1291           0 :                    MINVAL(HcoState%Grid%BXHEIGHT_M%Val), MAXVAL(HcoState%Grid%BXHEIGHT_M%Val)
        1292           0 :                 CALL HCO_MSG(HcoState%Config%Err,MSG)
        1293             :              ENDIF
        1294             :           ENDIF
        1295             :
        1296             :        ! PEDGE and/or TK not defined
        1297             :        ELSE
        1298           0 :           IF ( .NOT. FoundZSFC .AND. FIRST .AND. HcoState%amIRoot ) THEN
        1299             :              MSG = 'Cannot set surface height ZSFC. This may cause '      // &
        1300             :                    'some extensions to fail. HEMCO tries to calculate '   // &
        1301             :                    'ZSFC from surface pressure and air temperature, but ' // &
        1302           0 :                    'at least one of these variables seem to be missing.'
        1303           0 :              CALL HCO_WARNING( HcoState%Config%Err,MSG, RC, THISLOC=LOC )
        1304             :           ENDIF
        1305           0 :           IF ( .NOT. FoundBXHEIGHT .AND. FIRST .AND. HcoState%amIRoot ) THEN
        1306             :              MSG = 'Cannot set boxheights BXHEIGHT_M. This may cause '      // &
        1307             :                    'some extensions to fail. HEMCO tries to calculate '     // &
        1308             :                    'BXHEIGHT from pressure edges and air temperature, but ' // &
        1309           0 :                    'at least one of these variables seem to be missing.'
        1310           0 :              CALL HCO_WARNING( HcoState%Config%Err,MSG, RC, THISLOC=LOC )
        1311             :           ENDIF
        1312             :        ENDIF
        1313             :     ENDIF
        1314             :
        1315             :     ! ------------------------------------------------------------------
        1316             :     ! Wrap up and leave
        1317             :     ! ------------------------------------------------------------------
        1318             :
        1319             :     ! Verbose
        1320           0 :     IF ( Verb ) THEN
        1321           0 :        WRITE(MSG,*) 'Vertical grid calculations done.'
        1322           0 :        CALL HCO_MSG(HcoState%Config%Err,MSG,SEP2='-')
        1323             :     ENDIF
        1324             :
        1325             :     ! Cleanup
        1326           0 :     ThisTK => NULL()
        1327           0 :     IF ( ALLOCATED(TmpTK) ) DEALLOCATE(TmpTK)
        1328             :
        1329             :     ! Update first flag
        1330           0 :     FIRST = .FALSE.
        1331             :
        1332             :     ! Return w/ success
        1333           0 :     RC = HCO_SUCCESS
        1334             :
        1335           0 :   END SUBROUTINE HCO_CalcVertGrid
        1336             : !EOC
        1337             : !------------------------------------------------------------------------------
        1338             : !                   Harmonized Emissions Component (HEMCO)                    !
        1339             : !------------------------------------------------------------------------------
        1340             : !BOP
        1341             : !
        1342             : ! !IROUTINE: HCO_SetPBLm
        1343             : !
        1344             : ! !DESCRIPTION: Subroutine HCO\_SetPBLm sets the HEMCO PBL mixing height in
        1345             : ! meters. It first tries to read it from field 'FldName' (from the HEMCO data
        1346             : ! list), then to fill it from field 'PBLM', and then assigns the default value
        1347             : ! 'DefVal' to it.
        1348             : !\\
        1349             : !\\
        1350             : ! !INTERFACE:
        1351             : !
        1352           0 :   SUBROUTINE HCO_SetPBLm( HcoState, FldName, PBLM, DefVal, RC )
        1353             : !
        1354             : ! !USES
        1355             : !
        1356             :     USE HCO_Arr_Mod,      ONLY : HCO_ArrAssert
        1357             :     USE HCO_STATE_MOD,    ONLY : HCO_STATE
        1358             :     USE HCO_CALC_MOD,     ONLY : HCO_EvalFld
        1359             : !
        1360             : ! !INPUT PARAMETERS:
        1361             : !
        1362             :     TYPE(HCO_State),  POINTER                  :: HcoState    ! HEMCO state object
        1363             :     CHARACTER(LEN=*), OPTIONAL, INTENT(IN   )  :: FldName     ! field name
        1364             :     REAL(hp),         OPTIONAL, POINTER        :: PBLM(:,:)   ! pbl mixing height
        1365             :     REAL(hp),         OPTIONAL, INTENT(IN   )  :: DefVal      ! default value
        1366             : !
        1367             : ! !INPUT/OUTPUT PARAMETERS:
        1368             : !
        1369             :     INTEGER,          INTENT(INOUT)  :: RC
        1370             : !
        1371             : ! !REVISION HISTORY:
        1372             : !  28 Sep 2015 - C. Keller - Initial version
        1373             : !  See https://github.com/geoschem/hemco for complete history
        1374             : !EOP
        1375             : !------------------------------------------------------------------------------
        1376             : !BOC
        1377             : !
        1378             : ! !LOCAL VARIABLES:
        1379             : !
        1380             :     INTEGER                       :: NX, NY
        1381             :     LOGICAL                       :: FOUND
        1382             :     CHARACTER(LEN=255)            :: MSG
        1383             :     CHARACTER(LEN=255)            :: LOC = 'HCO_SetPBLm (hco_geotools_mod.F90)'
        1384             :
        1385             :     !-------------------------------
        1386             :     ! HCO_SetPBLm begins here
        1387             :     !-------------------------------
        1388             :
        1389             :     ! Init
        1390           0 :     FOUND = .FALSE.
        1391             :
        1392             :     ! Make sure size is ok. Allocate if unallocated.
        1393             :     CALL HCO_ArrAssert( HcoState%Grid%PBLHEIGHT, &
        1394           0 :                         HcoState%NX, HcoState%NY, RC )
        1395           0 :     IF ( RC /= HCO_SUCCESS ) THEN
        1396           0 :         CALL HCO_ERROR( 'ERROR 13', RC, THISLOC=LOC )
        1397           0 :         RETURN
        1398             :     ENDIF
        1399             :
        1400             :     ! Try to read from file first
        1401           0 :     IF ( PRESENT( FldName ) ) THEN
        1402             :        CALL HCO_EvalFld ( HcoState, FldName, &
        1403           0 :           HcoState%Grid%PBLHEIGHT%Val, RC, FOUND=FOUND )
        1404           0 :        IF ( RC /= HCO_SUCCESS ) THEN
        1405           0 :            CALL HCO_ERROR( 'ERROR 14', RC, THISLOC=LOC )
        1406           0 :            RETURN
        1407             :        ENDIF
        1408             :
        1409             :        ! Verbose
        1410           0 :        IF ( HcoState%amIRoot .AND. HCO_IsVerb(HcoState%Config%Err ) ) THEN
        1411           0 :           IF ( FOUND ) THEN
        1412           0 :              WRITE(MSG,*) 'HEMCO PBL heights obtained from field ',TRIM(FldName)
        1413           0 :              CALL HCO_MSG(HcoState%Config%Err,MSG,SEP2='-')
        1414             :           ENDIF
        1415             :        ENDIF
        1416             :     ENDIF
        1417             :
        1418             :     ! Pass 2D field if available
        1419           0 :     IF ( .not. FOUND ) THEN
        1420           0 :        IF ( PRESENT( PBLM ) ) THEN
        1421           0 :           IF ( ASSOCIATED(PBLM) ) THEN
        1422           0 :              NX = SIZE(PBLM,1)
        1423           0 :              NY = SIZE(PBLM,2)
        1424           0 :              IF ( NX /= HcoState%NX .OR. NY /= HcoState%NY ) THEN
        1425           0 :                 WRITE(MSG,*) 'Wrong PBLM array size: ', NX, NY, &
        1426           0 :                              '; should be: ', HcoState%NX, HcoState%NY
        1427           0 :                 CALL HCO_ERROR( MSG, RC, THISLOC=LOC )
        1428           0 :                 RETURN
        1429             :              ENDIF
        1430             :
        1431             :              ! Pass data
        1432           0 :              HcoState%Grid%PBLHEIGHT%Val = PBLM
        1433           0 :              FOUND                       = .TRUE.
        1434             :
        1435             :              ! Verbose
        1436           0 :              IF ( HcoState%amIRoot .AND. HCO_IsVerb(HcoState%Config%Err ) ) THEN
        1437           0 :                 WRITE(MSG,*) 'HEMCO PBL heights obtained from provided 2D field.'
        1438           0 :                 CALL HCO_MSG(HcoState%Config%Err,MSG,SEP2='-')
        1439             :              ENDIF
        1440             :           ENDIF
        1441             :        ENDIF
        1442             :     ENDIF
        1443             :
        1444             :     ! Finally, assign default value if field not yet set
        1445             : !    IF ( .NOT. FOUND .AND. PRESENT(DefVal) ) THEN
        1446           0 :     IF ( .NOT. FOUND ) THEN
        1447           0 :        IF ( PRESENT(DefVal) ) THEN
        1448             :           ! Make sure size is ok
        1449             :           CALL HCO_ArrAssert( HcoState%Grid%PBLHEIGHT, &
        1450           0 :                               HcoState%NX, HcoState%NY, RC )
        1451           0 :           IF ( RC /= HCO_SUCCESS ) THEN
        1452           0 :               CALL HCO_ERROR( 'ERROR 15', RC, THISLOC=LOC )
        1453           0 :               RETURN
        1454             :           ENDIF
        1455             :
        1456             :           ! Pass data
        1457           0 :           HcoState%Grid%PBLHEIGHT%Val = DefVal
        1458           0 :           FOUND                       = .TRUE.
        1459             :
        1460             :           ! Verbose
        1461           0 :           IF ( HcoState%amIRoot .AND. HCO_IsVerb(HcoState%Config%Err ) ) THEN
        1462           0 :              WRITE(MSG,*) 'HEMCO PBL heights uniformly set to ', DefVal
        1463           0 :              CALL HCO_MSG(HcoState%Config%Err,MSG,SEP2='-')
        1464             :           ENDIF
        1465             :        ENDIF
        1466             :     ENDIF
        1467             :
        1468             :     ! Error check
        1469           0 :     IF ( .NOT. FOUND ) THEN
        1470           0 :        WRITE(MSG,*) 'Cannot set PBL height: a valid HEMCO data field, ', &
        1471           0 :           'an explicit 2D field or a default value must be provided!'
        1472           0 :        CALL HCO_ERROR( MSG, RC, THISLOC=LOC )
        1473           0 :        RETURN
        1474             :     ENDIF
        1475             :
        1476             :     ! Return w/ success
        1477           0 :     RC = HCO_SUCCESS
        1478             :
        1479             :   END SUBROUTINE HCO_SetPBLm
        1480             : !EOC
        1481             : !!------------------------------------------------------------------------------
        1482             : !!                   Harmonized Emissions Component (HEMCO)                    !
        1483             : !!------------------------------------------------------------------------------
        1484             : !!BOP
        1485             : !!
        1486             : !! !SUBROUTINE: HCO_CalcPBLlev3D
        1487             : !!
        1488             : !! !DESCRIPTION:
        1489             : !!\\
        1490             : !!\\
        1491             : !! !INTERFACE:
        1492             : !!
        1493             : !  SUBROUTINE HCO_CalcPBLlev3D ( HcoState, PBLFRAC, RC )
        1494             : !!
        1495             : !! !USES
        1496             : !!
        1497             : !    USE HCO_Arr_Mod,      ONLY : HCO_ArrAssert
        1498             : !    USE HCO_STATE_MOD,    ONLY : HCO_STATE
        1499             : !!
        1500             : !! !INPUT PARAMETERS:
        1501             : !!
        1502             : !    TYPE(HCO_State), POINTER        :: HcoState        ! HEMCO state object
        1503             : !    REAL(hp),        POINTER        :: PBLFRAC(:,:,:)  ! planetary PBL fraction
        1504             : !!
        1505             : !! !INPUT/OUTPUT PARAMETERS:
        1506             : !!
        1507             : !    INTEGER,         INTENT(INOUT)  :: RC
        1508             : !!
        1509             : !! !REVISION HISTORY:
        1510             : !!  05 May 2016 - C. Keller - Initial version
        1511             : !!  See https://github.com/geoschem/hemco for complete history
        1512             : !!EOP
        1513             : !!------------------------------------------------------------------------------
        1514             : !!BOC
        1515             : !!
        1516             : !! !LOCAL VARIABLES:
        1517             : !!
        1518             : !    INTEGER             :: I,  J,  L
        1519             : !    CHARACTER(LEN=255)  :: MSG
        1520             : !    CHARACTER(LEN=255)  :: LOC = 'HCO_CalcPBLlev3D (hco_geotools_mod.F90)'
        1521             : !
        1522             : !    !-------------------------------
        1523             : !    ! HCO_CalcPBLlev3D begins here
        1524             : !    !-------------------------------
        1525             : !
        1526             : !    ! Check input array size
        1527             : !    IF ( SIZE(PBLFRAC,1) /= HcoState%NX .OR. &
        1528             : !         SIZE(PBLFRAC,2) /= HcoState%NY .OR. &
        1529             : !         SIZE(PBLFRAC,3) /= HcoState%NZ        ) THEN
        1530             : !       WRITE(MSG,*) 'Input array PBLFRAC has wrong horiz. dimensions: ', &
        1531             : !                     SIZE(PBLFRAC,1),SIZE(PBLFRAC,2),SIZE(PBLFRAC,3)
        1532             : !       CALL HCO_ERROR( MSG, RC, THISLOC=LOC )
        1533             : !       RETURN
        1534             : !    ENDIF
        1535             : !
        1536             : !    ! Make sure array is associated
        1537             : !    CALL HCO_ArrAssert( HcoState%Grid%PBL, HcoState%NX, HcoState%NY, RC )
        1538             : !    IF ( RC /= HCO_SUCCESS ) RETURN
        1539             : !
        1540             : !    ! Initialize values
        1541             : !    HcoState%Grid%PBL%Val = 1.0
        1542             : !
        1543             : !!$OMP PARALLEL DO                                                      &
        1544             : !!$OMP DEFAULT( SHARED )                                                &
        1545             : !!$OMP PRIVATE( I, J, L )                                               &
        1546             : !!$OMP SCHEDULE( DYNAMIC )
        1547             : !    DO J = 1, HcoState%NY
        1548             : !    DO I = 1, HcoState%NX
        1549             : !       ! Search for first level where PBL fraction is zero
        1550             : !       DO L = 1, HcoState%NZ
        1551             : !          IF ( PBLFRAC(I,J,L) > 0.0_hp ) CYCLE
        1552             : !          HcoState%Grid%PBL%Val(I,J) = MAX(L-1,1)
        1553             : !          EXIT
        1554             : !       ENDDO
        1555             : !    ENDDO
        1556             : !    ENDDO
        1557             : !!$OMP END PARALLEL DO
        1558             : !
        1559             : !    ! Return w/ success
        1560             : !    RC = HCO_SUCCESS
        1561             : !
        1562             : !  END SUBROUTINE HCO_CalcPBLlev3D
        1563             : !!EOC
        1564             : !!------------------------------------------------------------------------------
        1565             : !!                   Harmonized Emissions Component (HEMCO)                    !
        1566             : !!------------------------------------------------------------------------------
        1567             : !!BOP
        1568             : !!
        1569             : !! !SUBROUTINE: HCO_CalcPBLlev2D
        1570             : !!
        1571             : !! !DESCRIPTION:
        1572             : !!\\
        1573             : !!\\
        1574             : !! !INTERFACE:
        1575             : !!
        1576             : !  SUBROUTINE HCO_CalcPBLlev2D ( HcoState, PBLlev, RC )
        1577             : !!
        1578             : !! !USES
        1579             : !!
        1580             : !    USE HCO_Arr_Mod,      ONLY : HCO_ArrAssert
        1581             : !    USE HCO_STATE_MOD,    ONLY : HCO_STATE
        1582             : !!
        1583             : !! !INPUT PARAMETERS:
        1584             : !!
        1585             : !    TYPE(HCO_State), POINTER        :: HcoState        ! HEMCO state object
        1586             : !    INTEGER,         POINTER        :: PBLlev(:,:)     ! planetary PBL lev
        1587             : !!
        1588             : !! !INPUT/OUTPUT PARAMETERS:
        1589             : !!
        1590             : !    INTEGER,         INTENT(INOUT)  :: RC
        1591             : !!
        1592             : !! !REVISION HISTORY:
        1593             : !!  05 May 2016 - C. Keller - Initial version
        1594             : !!  See https://github.com/geoschem/hemco for complete history
        1595             : !!EOP
        1596             : !!------------------------------------------------------------------------------
        1597             : !!BOC
        1598             : !!
        1599             : !! !LOCAL VARIABLES:
        1600             : !!
        1601             : !    INTEGER             :: I,  J,  L
        1602             : !    CHARACTER(LEN=255)  :: MSG
        1603             : !    CHARACTER(LEN=255)  :: LOC = 'HCO_CalcPBLlev2D (hco_geotools_mod.F90)'
        1604             : !
        1605             : !    !-------------------------------
        1606             : !    ! HCO_CalcPBLlev2D begins here
        1607             : !    !-------------------------------
        1608             : !
        1609             : !    ! Make sure array is associated
        1610             : !    CALL HCO_ArrAssert( HcoState%Grid%PBL, HcoState%NX, HcoState%NY, RC )
        1611             : !    IF ( RC /= HCO_SUCCESS ) RETURN
        1612             : !
        1613             : !    ! Check input array size
        1614             : !    IF ( SIZE(PBLlev,1) /= HcoState%NX .OR. SIZE(PBLlev,2) /= HcoState%NY ) THEN
        1615             : !       WRITE(MSG,*) 'Input array PBLlev has wrong horiz. dimensions: ', &
        1616             : !              SIZE(PBLlev,1),SIZE(PBLlev,2)
        1617             : !       CALL HCO_ERROR( MSG, RC, THISLOC=LOC )
        1618             : !       RETURN
        1619             : !    ENDIF
        1620             : !
        1621             : !    ! Set values from PBLlev
        1622             : !!$OMP PARALLEL DO                                                      &
        1623             : !!$OMP DEFAULT( SHARED )                                                &
        1624             : !!$OMP PRIVATE( I, J )                                                  &
        1625             : !!$OMP SCHEDULE( DYNAMIC )
        1626             : !    DO J = 1, HcoState%NY
        1627             : !    DO I = 1, HcoState%NX
        1628             : !       HcoState%Grid%PBL%Val(I,J) = MAX(PBLlev(I,J),1)
        1629             : !    ENDDO
        1630             : !    ENDDO
        1631             : !!$OMP END PARALLEL DO
        1632             : !
        1633             : !    ! Return w/ success
        1634             : !    RC = HCO_SUCCESS
        1635             : !
        1636             : !  END SUBROUTINE HCO_CalcPBLlev2D
        1637             : !!EOC
        1638             : END MODULE HCO_GeoTools_Mod
        1639             : !EOM

Generated by: LCOV version 1.14