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

              Line data    Source code
           1             : !------------------------------------------------------------------------------
           2             : !                   Harmonized Emissions Component (HEMCO)                    !
           3             : !------------------------------------------------------------------------------
           4             : !BOP
           5             : !
           6             : ! !MODULE: hco_interp_mod.F90
           7             : !
           8             : ! !DESCRIPTION: Module HCO\_INTERP\_MOD contains routines to interpolate
           9             : ! input data onto the HEMCO grid. This module contains routine for
          10             : ! horizontal regridding between regular grids (MAP\_A2A), as well as
          11             : ! vertical interpolation amongst GEOS model levels (full <--> reduced).
          12             : !\\
          13             : !\\
          14             : ! Regridding is supported for concentration quantities (default) and
          15             : ! index-based values. For the latter, the values in the regridded grid
          16             : ! boxes correspond to the value of the original grid that contrbutes most
          17             : ! to the given box.
          18             : !\\
          19             : !\\
          20             : ! !INTERFACE:
          21             : !
          22             : MODULE HCO_Interp_Mod
          23             : !
          24             : ! !USES:
          25             : !
          26             :   USE HCO_Types_Mod
          27             :   USE HCO_Error_Mod
          28             :   USE HCO_State_Mod, ONLY : Hco_State
          29             :
          30             :   IMPLICIT NONE
          31             :   PRIVATE
          32             : !
          33             : ! !PUBLIC MEMBER FUNCTIONS:
          34             : !
          35             :   PUBLIC  :: ModelLev_Check
          36             :   PUBLIC  :: ModelLev_Interpolate
          37             :   PUBLIC  :: REGRID_MAPA2A
          38             : !
          39             : ! !PRIVATE MEMBER FUNCTIONS:
          40             : !
          41             :   PRIVATE :: COLLAPSE
          42             : !
          43             : ! !REVISION HISTORY:
          44             : !  30 Dec 2014 - C. Keller - Initialization
          45             : !  See https://github.com/geoschem/hemco for complete history
          46             : !EOP
          47             : !------------------------------------------------------------------------------
          48             : !BOC
          49             : !
          50             : ! !PRIVATE VARIABLES:
          51             : !
          52             :   ! AP parameter of native GEOS-5 grid. Needed to remap GEOS-5 data from native
          53             :   ! onto the reduced vertical grid.
          54             :   REAL(hp), TARGET :: G5_EDGE_NATIVE(73) = (/                          &
          55             :               0.000000e+00_hp, 4.804826e-02_hp, 6.593752e+00_hp, 1.313480e+01_hp, &
          56             :               1.961311e+01_hp, 2.609201e+01_hp, 3.257081e+01_hp, 3.898201e+01_hp, &
          57             :               4.533901e+01_hp, 5.169611e+01_hp, 5.805321e+01_hp, 6.436264e+01_hp, &
          58             :               7.062198e+01_hp, 7.883422e+01_hp, 8.909992e+01_hp, 9.936521e+01_hp, &
          59             :               1.091817e+02_hp, 1.189586e+02_hp, 1.286959e+02_hp, 1.429100e+02_hp, &
          60             :               1.562600e+02_hp, 1.696090e+02_hp, 1.816190e+02_hp, 1.930970e+02_hp, &
          61             :               2.032590e+02_hp, 2.121500e+02_hp, 2.187760e+02_hp, 2.238980e+02_hp, &
          62             :               2.243630e+02_hp, 2.168650e+02_hp, 2.011920e+02_hp, 1.769300e+02_hp, &
          63             :               1.503930e+02_hp, 1.278370e+02_hp, 1.086630e+02_hp, 9.236572e+01_hp, &
          64             :               7.851231e+01_hp, 6.660341e+01_hp, 5.638791e+01_hp, 4.764391e+01_hp, &
          65             :               4.017541e+01_hp, 3.381001e+01_hp, 2.836781e+01_hp, 2.373041e+01_hp, &
          66             :               1.979160e+01_hp, 1.645710e+01_hp, 1.364340e+01_hp, 1.127690e+01_hp, &
          67             :               9.292942e+00_hp, 7.619842e+00_hp, 6.216801e+00_hp, 5.046801e+00_hp, &
          68             :               4.076571e+00_hp, 3.276431e+00_hp, 2.620211e+00_hp, 2.084970e+00_hp, &
          69             :               1.650790e+00_hp, 1.300510e+00_hp, 1.019440e+00_hp, 7.951341e-01_hp, &
          70             :               6.167791e-01_hp, 4.758061e-01_hp, 3.650411e-01_hp, 2.785261e-01_hp, &
          71             :               2.113490e-01_hp, 1.594950e-01_hp, 1.197030e-01_hp, 8.934502e-02_hp, &
          72             :               6.600001e-02_hp, 4.758501e-02_hp, 3.270000e-02_hp, 2.000000e-02_hp, &
          73             :               1.000000e-02_hp /)
          74             :
          75             :   ! AP parameter of native 102-layer GISS grid
          76             :   REAL(hp), TARGET :: E102_EDGE_NATIVE(103) = (/                                            &
          77             :                 0.0000000,   2.7871507,   5.5743014,   8.3614521,  11.1486028, 13.9357536,  &
          78             :                16.7229043,  19.5100550,  22.2972057,  25.0843564,  27.8715071, 30.6586578,  &
          79             :                33.4458085,  36.2329593,  39.0201100,  41.8087123,  44.6089278, 47.4534183,  &
          80             :                50.4082336,  53.5662786,  57.0095710,  60.7533531,  64.7323011, 68.8549615,  &
          81             :                73.0567364,  77.2969797,  81.5364973,  85.7346430,  89.8565776, 93.8754457,  &
          82             :                97.7709243, 101.5277712, 105.1350991, 108.5878272, 111.8859556, 115.0302100, &
          83             :               118.0249453, 120.8854039, 123.6326345, 126.2811535, 128.8360417, 131.2987506, &
          84             :               133.6736353, 135.9708571, 138.2013035, 140.3700552, 142.4814670, 144.5457005, &
          85             :               146.5692881, 148.5464231, 150.4712991, 152.3497225, 154.1875000, 144.5468750, &
          86             :               135.1875000, 126.0781250, 117.1914062, 108.5859375, 100.3671875, 92.5898438,  &
          87             :                85.2265625,  78.2226562,  71.5546875,  65.2226562,  59.2226562, 53.5546875,  &
          88             :                48.2226562,  43.2226562,  38.5546875,  34.2226562,  30.2226562, 26.5507812,  &
          89             :                23.1875000,  20.0781250,  17.1896562,  14.5684375,  12.2865742, 10.3573086,  &
          90             :                 8.7353750,   7.3664922,   6.2100156,   5.2343633,   4.4119297, 3.7186797,   &
          91             :                 3.1341479,   2.6404328,   2.2207877,   1.8587369,   1.5477125, 1.2782115,   &
          92             :                 1.0427319,   0.8367716,   0.6514691,   0.4772511,   0.3168814, 0.1785988,   &
          93             :                 0.1000000,   0.0560000,   0.0320000,   0.0180000,   0.0100000, 0.0050000,   &
          94             :                 0.0020000 /)
          95             :
          96             : CONTAINS
          97             : !EOC
          98             : !------------------------------------------------------------------------------
          99             : !                   Harmonized Emissions Component (HEMCO)                    !
         100             : !------------------------------------------------------------------------------
         101             : !BOP
         102             : !
         103             : ! !IROUTINE: Regrid_MAPA2A
         104             : !
         105             : ! !DESCRIPTION: Subroutine Regrid\_MAPA2A regrids input array NcArr onto
         106             : ! the simulation grid and stores the data in list container Lct. Horizontal
         107             : ! regridding is performed using MAP\_A2A algorithm. Vertical interpolation
         108             : ! between GEOS levels (full vs. reduced, GEOS-5 vs. GEOS-4), is also
         109             : ! supported.
         110             : !\\
         111             : !\\
         112             : ! This routine can remap concentrations and index-based quantities.
         113             : !\\
         114             : !\\
         115             : ! !INTERFACE:
         116             : !
         117           0 :   SUBROUTINE REGRID_MAPA2A( HcoState, NcArr, LonE, LatE, Lct, RC )
         118             : !
         119             : ! !USES:
         120             : !
         121             :     USE HCO_REGRID_A2A_Mod, ONLY : MAP_A2A
         122             :     USE HCO_FileData_Mod,   ONLY : FileData_ArrCheck
         123             :     USE HCO_UNIT_MOD,       ONLY : HCO_IsIndexData
         124             : !
         125             : ! !INPUT PARAMETERS:
         126             : !
         127             :     TYPE(HCO_State),  POINTER        :: HcoState          ! HEMCO state object
         128             :     REAL(sp),         POINTER        :: NcArr(:,:,:,:)    ! 4D input data
         129             :     REAL(hp),         POINTER        :: LonE(:)           ! Input grid longitude edges
         130             :     REAL(hp),         POINTER        :: LatE(:)           ! Input grid latitude edges
         131             : !
         132             : ! !INPUT/OUTPUT PARAMETERS:
         133             : !
         134             :     TYPE(ListCont),   POINTER        :: Lct               ! HEMCO list container
         135             :     INTEGER,          INTENT(INOUT)  :: RC                ! Success or failure?
         136             : !
         137             : ! !REVISION HISTORY:
         138             : !  03 Feb 2015 - C. Keller   - Initial version
         139             : !  See https://github.com/geoschem/hemco for complete history
         140             : !EOP
         141             : !------------------------------------------------------------------------------
         142             : !BOC
         143             : !
         144             : ! !LOCAL VARIABLES:
         145             : !
         146             :     INTEGER                 :: nLonEdge, nLatEdge
         147             :     INTEGER                 :: NX, NY, NZ, NLEV, NTIME, NCELLS
         148             :     INTEGER                 :: I, J, L, T, AS, I2
         149             :     INTEGER                 :: nIndex
         150           0 :     REAL(sp), ALLOCATABLE   :: LonEdgeI(:)
         151           0 :     REAL(sp), ALLOCATABLE   :: LatEdgeI(:)
         152           0 :     REAL(sp)                :: LonEdgeO(HcoState%NX+1)
         153           0 :     REAL(sp)                :: LatEdgeO(HcoState%NY+1)
         154             :
         155           0 :     REAL(sp), POINTER       :: ORIG_2D(:,:)
         156           0 :     REAL(sp), POINTER       :: REGR_2D(:,:)
         157           0 :     REAL(sp), POINTER       :: REGR_4D(:,:,:,:)
         158             :
         159           0 :     REAL(sp), ALLOCATABLE, TARGET :: FRACS(:,:,:,:)
         160           0 :     REAL(hp), ALLOCATABLE         :: REGFRACS(:,:,:,:)
         161           0 :     REAL(hp), ALLOCATABLE         :: MAXFRACS(:,:,:,:)
         162           0 :     REAL(hp), ALLOCATABLE         :: INDECES(:,:,:,:)
         163           0 :     REAL(hp), ALLOCATABLE         :: UNIQVALS(:)
         164             :     REAL(hp)                      :: IVAL
         165             :     LOGICAL                       :: IsIndex
         166             :
         167             :     LOGICAL                 :: VERB
         168             :     CHARACTER(LEN=255)      :: MSG
         169             :     CHARACTER(LEN=255)      :: LOC = 'ModelLev_Interpolate (hco_interp_mod.F90)'
         170             :
         171             :     !=================================================================
         172             :     ! REGRID_MAPA2A begins here
         173             :     !=================================================================
         174             :
         175             :     ! Init
         176           0 :     ORIG_2D => NULL()
         177           0 :     REGR_2D => NULL()
         178           0 :     REGR_4D => NULL()
         179             :
         180             :     ! Check for verbose mode
         181           0 :     verb = HCO_IsVerb( HcoState%Config%Err )
         182             :
         183             :     ! get longitude / latitude sizes
         184           0 :     nLonEdge = SIZE(LonE,1)
         185           0 :     nLatEdge = SIZE(LatE,1)
         186             :
         187             :     ! Write input grid edges to shadow variables so that map_a2a accepts them
         188             :     ! as argument.
         189             :     ! Also, for map_a2a, latitudes have to be sines...
         190           0 :     ALLOCATE(LonEdgeI(nlonEdge), LatEdgeI(nlatEdge), STAT=AS )
         191           0 :     IF ( AS /= 0 ) THEN
         192           0 :        CALL HCO_ERROR( 'alloc error LonEdgeI/LatEdgeI', RC, THISLOC=LOC )
         193           0 :        RETURN
         194             :     ENDIF
         195           0 :     LonEdgeI(:) = LonE
         196           0 :     LatEdgeI(:) = SIN( LatE * HcoState%Phys%PI_180 )
         197             :
         198             :     ! Get output grid edges from HEMCO state
         199           0 :     LonEdgeO(:) = HcoState%Grid%XEDGE%Val(:,1)
         200           0 :     LatEdgeO(:) = HcoState%Grid%YSIN%Val(1,:)
         201             :
         202             :     ! Get input array sizes
         203           0 :     NX     = size(ncArr,1)
         204           0 :     NY     = size(ncArr,2)
         205           0 :     NLEV   = size(ncArr,3)
         206           0 :     NTIME  = size(ncArr,4)
         207           0 :     NCELLS = NX * NY * NLEV * NTIME
         208             :
         209             :     ! Are these index-based data? If so, need to remap the fraction (1 or 0)
         210             :     ! of every value independently. For every grid box, the value with the
         211             :     ! highest overlap (closest to 1) is taken.
         212           0 :     IsIndex = HCO_IsIndexData(Lct%Dct%Dta%OrigUnit)
         213             :
         214           0 :     IF ( IsIndex ) THEN
         215             :
         216             :        ! Allocate working arrays:
         217             :        ! - FRACS contains the fractions on the original grid. These are
         218             :        !   binary (1 or 0).
         219             :        ! - MAXFRACS stores the highest used fraction for each output grid
         220             :        !   box. Will be updated continously.
         221             :        ! - INDECES is the output array holding the index-based remapped
         222             :        !   values. Will be updated continuously.
         223             :        ! - UNIQVALS is a vector holding all unique values of the input
         224             :        !   array (NINDEX is the number of unique values).
         225             :        !
         226             :        ! ckeller, 9/24/15: Extend vertical axis of MAXFRACS, REGFRACS, and
         227             :        ! INDECES to HcoState%NZ+1 for fields that are on edges instead of
         228             :        ! mid-points.
         229           0 :        ALLOCATE( FRACS(NX,NY,NLEV,NTIME), STAT=AS )
         230           0 :        IF ( AS /= 0 ) THEN
         231           0 :           CALL HCO_ERROR( 'alloc error FRACS', RC, THISLOC=LOC )
         232           0 :           RETURN
         233             :        ENDIF
         234           0 :        ALLOCATE( MAXFRACS(HcoState%NX,HcoState%NY,HcoState%NZ+1,NTIME), STAT=AS )
         235           0 :        IF ( AS /= 0 ) THEN
         236           0 :           CALL HCO_ERROR( 'alloc error MAXFRACS', RC, THISLOC=LOC )
         237           0 :           RETURN
         238             :        ENDIF
         239           0 :        ALLOCATE( REGFRACS(HcoState%NX,HcoState%NY,HcoState%NZ+1,NTIME), STAT=AS )
         240           0 :        IF ( AS /= 0 ) THEN
         241           0 :           CALL HCO_ERROR( 'alloc error INDECES', RC, THISLOC=LOC )
         242           0 :           RETURN
         243             :        ENDIF
         244           0 :        ALLOCATE( INDECES(HcoState%NX,HcoState%NY,HcoState%NZ+1,NTIME), STAT=AS )
         245           0 :        IF ( AS /= 0 ) THEN
         246           0 :           CALL HCO_ERROR( 'alloc error INDECES', RC, THISLOC=LOC )
         247           0 :           RETURN
         248             :        ENDIF
         249           0 :        ALLOCATE( UNIQVALS(NCELLS), STAT=AS )
         250           0 :        IF ( AS /= 0 ) THEN
         251           0 :           CALL HCO_ERROR( 'alloc error INDECES', RC, THISLOC=LOC )
         252           0 :           RETURN
         253             :        ENDIF
         254           0 :        FRACS    = 0.0_sp
         255           0 :        REGFRACS = 0.0_hp
         256           0 :        MAXFRACS = 0.0_hp
         257           0 :        INDECES  = 0.0_hp
         258           0 :        UNIQVALS = 0.0_hp
         259             :
         260             :        ! Get unique values. Loop over all input data values and add
         261             :        ! them to UNIQVALS vector if UNIQVALS doesn't hold that same value
         262             :        ! yet.
         263           0 :        NINDEX = 0
         264           0 :        DO T = 1, NTIME
         265           0 :        DO L = 1, NLEV
         266           0 :        DO J = 1, NY
         267           0 :        DO I = 1, NX
         268             :
         269             :           ! Current value
         270           0 :           IVAL = NcArr(I,J,L,T)
         271             :
         272             :           ! Check if value already exists in UNIQVALS
         273           0 :           IF ( NINDEX > 0 ) THEN
         274           0 :              IF ( ANY(UNIQVALS(1:NINDEX) == IVAL) ) CYCLE
         275             :           ENDIF
         276             :
         277             :           ! Add to UNIQVALS
         278           0 :           NINDEX = NINDEX + 1
         279           0 :           UNIQVALS(NINDEX) = IVAL
         280             :        ENDDO
         281             :        ENDDO
         282             :        ENDDO
         283             :        ENDDO
         284             :
         285             :        ! Verbose mode
         286           0 :        IF ( verb ) THEN
         287           0 :           MSG = 'Do index based regridding for field ' // TRIM(Lct%Dct%cName)
         288           0 :           CALL HCO_MSG(HcoState%Config%Err,MSG)
         289           0 :           WRITE(MSG,*) '   - Number of indeces: ', NINDEX
         290           0 :           CALL HCO_MSG(HcoState%Config%Err,MSG)
         291             :        ENDIF
         292             :
         293             :     ELSE
         294           0 :        NINDEX = 1
         295             :     ENDIF
         296             :
         297             :     ! Define array to put horizontally regridded data onto. If this
         298             :     ! is 3D data, we first regrid all vertical levels horizontally
         299             :     ! and then pass these data to the list container. In this second
         300             :     ! step, levels may be deflated/collapsed.
         301             :
         302             :     ! 2D data is directly passed to the data container
         303           0 :     IF ( Lct%Dct%Dta%SpaceDim <= 2 ) THEN
         304             :        CALL FileData_ArrCheck( HcoState%Config, Lct%Dct%Dta, &
         305           0 :                                HcoState%NX, HcoState%NY, NTIME, RC )
         306           0 :        IF ( RC /= 0 ) RETURN
         307             :     ENDIF
         308             :
         309             :     ! 3D data and index data is first written into a temporary array,
         310             :     ! REGR_4D.
         311           0 :     IF ( Lct%Dct%Dta%SpaceDim == 3 .OR. IsIndex ) THEN
         312           0 :        ALLOCATE( REGR_4D(HcoState%NX,HcoState%NY,NLEV,NTIME), STAT=AS )
         313             :        IF ( AS /= 0 ) THEN
         314           0 :           CALL HCO_ERROR( 'alloc error REGR_4D', RC, THISLOC=LOC )
         315           0 :           RETURN
         316             :        ENDIF
         317           0 :        REGR_4D = 0.0_hp
         318             :     ENDIF
         319             :
         320             :     ! Do regridding for every index value. If it's not index data, this loop
         321             :     ! is executed only once (NINDEX=1).
         322           0 :     DO I = 1, NINDEX
         323             :
         324             :        ! For index based data, create fractions array for the given index.
         325           0 :        IF ( IsIndex ) THEN
         326           0 :           IVAL = UNIQVALS(I)
         327           0 :           WHERE( ncArr == IVAL )
         328             :              FRACS = 1.0_sp
         329             :           ELSEWHERE
         330             :              FRACS = 0.0_sp
         331             :           END WHERE
         332             :        ENDIF
         333             :
         334             :        ! Regrid horizontally
         335           0 :        DO T = 1, NTIME
         336           0 :        DO L = 1, NLEV
         337             :
         338             :           ! Point to 2D slices to be regridded:
         339             :           ! - Original 2D array
         340           0 :           IF ( IsIndex ) THEN
         341           0 :             ORIG_2D => FRACS(:,:,L,T)
         342             :           ELSE
         343           0 :             ORIG_2D => ncArr(:,:,L,T)
         344             :           ENDIF
         345             :
         346             :           ! - Regridded 2D array
         347           0 :           IF ( Lct%Dct%Dta%SpaceDim <= 2 .AND. .NOT. IsIndex ) THEN
         348           0 :              REGR_2D => Lct%Dct%Dta%V2(T)%Val(:,:)
         349             :           ELSE
         350           0 :              REGR_2D => REGR_4D(:,:,L,T)
         351             :           ENDIF
         352             :
         353             :           ! Do the regridding
         354             :           CALL MAP_A2A( NX,      NY, LonEdgeI,    LatEdgeI, ORIG_2D,  &
         355             :                         HcoState%NX, HcoState%NY, LonEdgeO, LatEdgeO, &
         356           0 :                         REGR_2D, 0, 0, HCO_MISSVAL )
         357           0 :           ORIG_2D => NULL()
         358           0 :           REGR_2D => NULL()
         359             :
         360             :        ENDDO !L
         361             :        ENDDO !T
         362             :
         363             :        ! Eventually collapse levels onto simulation levels.
         364           0 :        IF ( Lct%Dct%Dta%SpaceDim == 3 ) THEN
         365           0 :           CALL ModelLev_Interpolate( HcoState, REGR_4D, Lct, RC )
         366           0 :           IF ( RC /= HCO_SUCCESS ) THEN
         367           0 :               CALL HCO_ERROR( 'ERROR 0', RC, THISLOC=LOC )
         368           0 :               RETURN
         369             :           ENDIF
         370             :        ENDIF
         371             :
         372             :        ! For index based data, map fractions back to corresponding value.
         373             :        ! Array INDECES holds the index-based remapped values. Set INDECES
         374             :        ! to current index value in every grid box where the regridded
         375             :        ! fraction of this index is higher than any previous fraction
         376             :        ! (array MAXFRACS stores the highest used fraction in each grid box).
         377           0 :        IF ( IsIndex ) THEN
         378             :
         379             :           ! Reset
         380           0 :           REGFRACS = 0.0_hp
         381             :
         382             :           ! 3D data written to Lct needs to be mapped back onto REGR_4D.
         383           0 :           IF ( Lct%Dct%Dta%SpaceDim == 3 ) THEN
         384           0 :              DO T = 1, NTIME
         385           0 :                 NZ = SIZE(Lct%Dct%Dta%V3(T)%Val,3)
         386           0 :                 REGFRACS(:,:,1:NZ,T) = Lct%Dct%Dta%V3(T)%Val(:,:,:)
         387             :              ENDDO
         388             :           ELSE
         389           0 :              REGFRACS(:,:,1:NLEV,:) = REGR_4D(:,:,:,:)
         390             :           ENDIF
         391             :
         392             :           ! REGR_4D are the remapped fractions.
         393           0 :           DO T  = 1, NTIME
         394           0 :           DO L  = 1, HcoState%NZ
         395           0 :           DO J  = 1, HcoState%NY
         396           0 :           DO I2 = 1, HcoState%NX
         397           0 :              IF ( REGFRACS(I2,J,L,T) > MAXFRACS(I2,J,L,T) ) THEN
         398           0 :                 MAXFRACS(I2,J,L,T) = REGR_4D(I2,J,L,T)
         399           0 :                 INDECES (I2,J,L,T) = IVAL
         400             :              ENDIF
         401             :           ENDDO
         402             :           ENDDO
         403             :           ENDDO
         404             :           ENDDO
         405             :
         406             : !------------------------------------------------------------------------------
         407             : ! Prior to 9/29/16:
         408             : !          ! This code is preblematic in Gfortran.  Replace it with the
         409             : !          ! explicit DO loops above.  Leave this here for reference.
         410             : !          ! (sde, bmy, 9/21/16)
         411             : !          WHERE ( REGFRACS > MAXFRACS )
         412             : !             MAXFRACS = REGR_4D
         413             : !             INDECES  = IVAL
         414             : !          END WHERE
         415             : !------------------------------------------------------------------------------
         416             :        ENDIF
         417             :
         418             :     ENDDO !I
         419             :
         420             :     ! For index values, pass index data to data container.
         421           0 :     IF ( IsIndex ) THEN
         422           0 :        IF ( Lct%Dct%Dta%SpaceDim == 3 ) THEN
         423           0 :           DO T = 1, NTIME
         424           0 :              NZ = SIZE(Lct%Dct%Dta%V3(T)%Val,3)
         425           0 :              Lct%Dct%Dta%V3(T)%Val(:,:,:) = INDECES(:,:,1:NZ,T)
         426             :           ENDDO
         427             :        ELSE
         428           0 :           DO T = 1, NTIME
         429           0 :              Lct%Dct%Dta%V2(T)%Val(:,:)   = INDECES(:,:,1,T)
         430             :           ENDDO
         431             :        ENDIF
         432             :     ENDIF
         433             :
         434             :     ! Cleanup
         435           0 :     DEALLOCATE(LonEdgeI, LatEdgeI)
         436           0 :     IF ( ASSOCIATED( REGR_4D  ) ) DEALLOCATE( REGR_4D  )
         437           0 :     IF ( ALLOCATED ( FRACS    ) ) DEALLOCATE( FRACS    )
         438           0 :     IF ( ALLOCATED ( REGFRACS ) ) DEALLOCATE( REGFRACS )
         439           0 :     IF ( ALLOCATED ( MAXFRACS ) ) DEALLOCATE( MAXFRACS )
         440           0 :     IF ( ALLOCATED ( INDECES  ) ) DEALLOCATE( INDECES  )
         441           0 :     IF ( ALLOCATED ( UNIQVALS ) ) DEALLOCATE( UNIQVALS )
         442             :
         443             :     ! Return w/ success
         444           0 :     RC = HCO_SUCCESS
         445             :
         446           0 :   END SUBROUTINE REGRID_MAPA2A
         447             : !EOC
         448             : !------------------------------------------------------------------------------
         449             : !                   Harmonized Emissions Component (HEMCO)                    !
         450             : !------------------------------------------------------------------------------
         451             : !BOP
         452             : !
         453             : ! !IROUTINE: ModelLev_Check
         454             : !
         455             : ! !DESCRIPTION: Subroutine ModelLev\_Check checks if the passed number of
         456             : ! vertical levels indicates that these are model levels or not.
         457             : !\\
         458             : !\\
         459             : ! !INTERFACE:
         460             : !
         461           0 :   SUBROUTINE ModelLev_Check( HcoState, nLev, IsModelLev, RC )
         462             : !
         463             : ! !USES:
         464             : !
         465             :     USE HCO_FileData_Mod,   ONLY : FileData_ArrCheck
         466             : !
         467             : ! !INPUT PARAMETERS:
         468             : !
         469             :     TYPE(HCO_State),  POINTER        :: HcoState          ! HEMCO state object
         470             :     INTEGER,          INTENT(IN   )  :: nlev              ! number of levels
         471             : !
         472             : ! !INPUT/OUTPUT PARAMETERS:
         473             : !
         474             :     LOGICAL,          INTENT(INOUT)  :: IsModelLev        ! Are these model levels?
         475             :     INTEGER,          INTENT(INOUT)  :: RC                ! Success or failure?
         476             : !
         477             : ! !REVISION HISTORY:
         478             : !  29 Sep 2015 - C. Keller   - Initial version
         479             : !  See https://github.com/geoschem/hemco for complete history
         480             : !EOP
         481             : !------------------------------------------------------------------------------
         482             : !BOC
         483             : !
         484             : ! !LOCAL VARIABLES:
         485             : !
         486             :     INTEGER                 :: nz
         487             :
         488             :     !=================================================================
         489             :     ! ModelLev_Check begins here
         490             :     !=================================================================
         491             :
         492             :     ! Assume success until otherwise
         493           0 :     RC = HCO_SUCCESS
         494             :
         495             :     ! Shadow number of vertical levels on grid
         496           0 :     nz = HcoState%NZ
         497             :
         498             :     ! Assume model levels if input data levels correspond to # of grid
         499             :     ! levels or levels + 1 (edges)
         500           0 :     IF ( nlev == nz .OR. nlev == nz + 1 ) THEN
         501           0 :        IsModelLev = .TRUE.
         502             :
         503             :    ! If input is 72 layer (or 3/11/36 layer) and output is 47 layer
         504           0 :     ELSEIF ( nz == 47 ) THEN
         505           0 :        IsModelLev = ( nlev == 72 .OR. nlev == 73 .OR. nlev == 36 .OR. nlev == 3 .OR. nlev == 11)
         506             :
         507             :     ! If input is 102 layer and output is 74 layer
         508           0 :     ELSEIF ( nz == 74 ) THEN
         509           0 :        IsModelLev = ( nlev == 102 .OR. nlev == 103 )
         510             :
         511             :     ! If input is 47 layer (or 3/11/36 layer) and output is 72 layer
         512           0 :     ELSEIF ( nz == 72 ) THEN
         513           0 :        IsModelLev = ( nlev == 47 .OR. nlev == 48 .OR. nlev == 36 .OR. nlev == 3 .OR. nlev == 11)
         514             :
         515             :     ELSE
         516           0 :       IsModelLev = .FALSE.
         517             :     ENDIF
         518             :
         519           0 :   END SUBROUTINE ModelLev_Check
         520             : !EOC
         521             : !------------------------------------------------------------------------------
         522             : !                   Harmonized Emissions Component (HEMCO)                    !
         523             : !------------------------------------------------------------------------------
         524             : !BOP
         525             : !
         526             : ! !IROUTINE: ModelLev_Interpolate
         527             : !
         528             : ! !DESCRIPTION: Subroutine ModelLev\_Interpolate puts 3D data from an
         529             : ! arbitrary number of model levels onto the vertical levels of the simulation
         530             : ! grid. Since the input data is already on model levels, this is only to
         531             : ! collapse fields between native/reduced vertical levels, e.g. from
         532             : ! 72 native GEOS-5 levels onto the reduced 47 levels.
         533             : !
         534             : !
         535             : ! The input data (REGR\_4D) is expected to be already regridded horizontally.
         536             : ! The 4th dimension of REGR\_4D denotes time.
         537             : !
         538             : !
         539             : ! The 3rd dimension of REGR\_3D holds the vertical levels. It is assumed that
         540             : ! these are model levels, starting at the surface (level 1). If the input
         541             : ! data holds 72/73 input levels, this is interpreted as native data and will
         542             : ! be collapsed onto the reduced GEOS-5 grid. If the input holds 102/103 input
         543             : ! levels, this is interpreted as native data and will be collapsed onto the
         544             : ! reduced GISS grid. If the input holds 47/48 input levels, this is interpreted
         545             : ! as reduced GEOS-5 data and it will be inflated to the native GEOS-5 grid
         546             : ! (with a warning, as this is not recommended). If the input holds N input levels,
         547             : ! (where N = 3, 11, or 36 to account for NEI and AEIC emissions), this is assumed
         548             : ! to be the first N levels of the GEOS-5 grid, meaning they will be written as
         549             : ! levels 1-N of a 47 or 72 level output grid (with the remaining values left to
         550             : ! be zero) (nbalasus, 8/29/2023).
         551             : !
         552             : !
         553             : ! Currently, this routine can remap the following combinations:
         554             : !
         555             : ! * Native GEOS-5 onto reduced GEOS-5 (72 --> 47 levels, 73 --> 48 edges)
         556             : ! * Native GISS onto reduced GISS (102 --> 74 levels, 103 --> 75 edges)
         557             : ! * Reduced GEOS-5 onto native GEOS-5 (47 --> 72 levels, 48 --> 73 edges)
         558             : ! * N (N = 3/11/36) levels onto native/reduced GEOS-5 (N --> levels 1-N levels of 47/72 level grid, rest are 0)
         559             : !
         560             : ! !INTERFACE:
         561             : !
         562           0 :   SUBROUTINE ModelLev_Interpolate( HcoState, REGR_4D, Lct, RC )
         563             : !
         564             : ! !USES:
         565             : !
         566             :     USE HCO_FileData_Mod,   ONLY : FileData_ArrCheck
         567             : !
         568             : ! !INPUT PARAMETERS:
         569             : !
         570             :     TYPE(HCO_State),  POINTER        :: HcoState          ! HEMCO state object
         571             :     REAL(sp),         POINTER        :: REGR_4D(:,:,:,:)  ! 4D input data
         572             : !
         573             : ! !INPUT/OUTPUT PARAMETERS:
         574             : !
         575             :     TYPE(ListCont),   POINTER        :: Lct               ! HEMCO list container
         576             :     INTEGER,          INTENT(INOUT)  :: RC                ! Success or failure?
         577             : !
         578             : ! !REVISION HISTORY:
         579             : !  30 Dec 2014 - C. Keller   - Initial version
         580             : !  See https://github.com/geoschem/hemco for complete history
         581             : !EOP
         582             : !------------------------------------------------------------------------------
         583             : !BOC
         584             : !
         585             : ! !LOCAL VARIABLES:
         586             : !
         587             :     INTEGER                 :: nx, ny, nz, nt
         588             :     INTEGER                 :: fineIDX, coarseIDX
         589             :     INTEGER                 :: minlev, nlev, nout
         590             :     INTEGER                 :: L, T, NL, I
         591             :     INTEGER                 :: OS
         592             :     LOGICAL                 :: verb, infl, clps
         593             :     LOGICAL                 :: DONE
         594             :     CHARACTER(LEN=255)      :: MSG, LOC
         595             :
         596             :     !=================================================================
         597             :     ! ModelLev_Interpolate begins here
         598             :     !=================================================================
         599           0 :     LOC = 'ModelLev_Interpolate (HCO_INTERP_MOD.F90)'
         600             :
         601             :     ! Enter
         602           0 :     CALL HCO_ENTER (HcoState%Config%Err, LOC, RC )
         603           0 :     IF ( RC /= HCO_SUCCESS ) THEN
         604           0 :         CALL HCO_ERROR( 'ERROR 1', RC, THISLOC=LOC )
         605           0 :         RETURN
         606             :     ENDIF
         607             :
         608             :     ! Check for verbose mode
         609           0 :     verb = HCO_IsVerb( HcoState%Config%Err )
         610           0 :     IF ( verb ) THEN
         611           0 :        MSG = 'Vertically interpolate model levels: '//TRIM(Lct%Dct%cName)
         612           0 :        CALL HCO_MSG(HcoState%Config%Err,MSG)
         613             :     ENDIF
         614             :
         615             :     ! Get HEMCO grid dimensions
         616           0 :     nx = HcoState%NX
         617           0 :     ny = HcoState%NY
         618           0 :     nz = HcoState%NZ
         619             :
         620             :     ! Input data must be on horizontal HEMCO grid
         621           0 :     IF ( SIZE(REGR_4D,1) /= nx ) THEN
         622           0 :        WRITE(MSG,*) 'x dimension mismatch ', TRIM(Lct%Dct%cName), &
         623           0 :           ': ', nx, SIZE(REGR_4D,1)
         624           0 :        CALL HCO_ERROR( MSG, RC )
         625           0 :        RETURN
         626             :     ENDIF
         627           0 :     IF ( SIZE(REGR_4D,2) /= ny ) THEN
         628           0 :        WRITE(MSG,*) 'y dimension mismatch ', TRIM(Lct%Dct%cName), &
         629           0 :           ': ', ny, SIZE(REGR_4D,2)
         630           0 :        CALL HCO_ERROR( MSG, RC )
         631           0 :        RETURN
         632             :     ENDIF
         633             :
         634             :     ! Get vertical and time dimension of input data
         635           0 :     nlev = SIZE(REGR_4D,3)
         636           0 :     nt   = SIZE(REGR_4D,4)
         637             :
         638             :     ! Check to make sure ModelLev_Interpolate should have been called
         639             :     IF ( ( ( nlev == nz ) .OR. ( nlev == nz+1 ) ) .OR. &                          ! write data without doing anything
         640             :          ( ( nz == 47 ) .AND. ( ( nlev == 72 ) .OR. ( nlev == 73 ) ) ) .OR. &     ! collapse native to reduced GEOS-5
         641             :          ( ( nz == 74 ) .AND. ( ( nlev == 102 ) .OR. ( nlev == 103 ) ) ) .OR. &   ! collapse native to reduced GISS
         642           0 :          ( ( nz == 72 ) .AND. ( ( nlev == 47 ) .OR. ( nlev == 48 ) ) ) .OR. &     ! inflate reduced to native GEOS-5
         643             :          ( ( ( nz == 72 ) .OR. ( nz == 47 ) ) .AND. ( ( nlev == 3 ) .OR. ( nlev == 11 ) .OR. ( nlev == 36 ) ) ) ) THEN  ! write 3/11/36 levels to reduced/native GEOS-5
         644             :          ! do nothing
         645             :     ELSE
         646           0 :       WRITE(MSG,*) 'ModelLev_Interpolate was called but MESSy should have been used: ',TRIM(Lct%Dct%cName)
         647           0 :       CALL HCO_ERROR( MSG, RC )
         648           0 :       RETURN
         649             :     ENDIF
         650             :
         651             :     ! Vertical interpolation done?
         652           0 :     DONE = .FALSE.
         653             :
         654             :     !===================================================================
         655             :     ! If no vertical interpolation is needed, then (1) save the 4D
         656             :     ! input data array to to the HEMCO list container object and
         657             :     ! (2) exit this subroutine.
         658             :     !===================================================================
         659           0 :     IF ( ( nlev == nz ) .OR. ( nlev == nz+1 ) ) THEN
         660             :
         661           0 :        CALL FileData_ArrCheck( HcoState%Config, Lct%Dct%Dta, nx, ny, nlev, nt, RC )
         662           0 :        IF ( RC /= HCO_SUCCESS ) THEN
         663           0 :            CALL HCO_ERROR( 'ERROR 2', RC, THISLOC=LOC )
         664           0 :            RETURN
         665             :        ENDIF
         666             :
         667           0 :        DO T = 1, nt
         668           0 :           Lct%Dct%Dta%V3(T)%Val(:,:,:) = REGR_4D(:,:,:,T)
         669             :        ENDDO
         670             :
         671             :        ! Verbose
         672           0 :        IF ( HCO_IsVerb( HcoState%Config%Err ) ) THEN
         673           0 :           MSG = '# of input levels = # of output levels - passed as is.'
         674           0 :           CALL HCO_MSG(HcoState%Config%Err,MSG)
         675             :        ENDIF
         676             :
         677             :        ! Done!
         678             :        DONE = .TRUE.
         679             :     ENDIF
         680             :
         681             :     !===================================================================
         682             :     ! Do vertical regridding:
         683             :     !===================================================================
         684             :     IF ( .NOT. DONE ) THEN
         685             :
         686             :        !----------------------------------------------------------------
         687             :        ! Native to reduced GEOS-5 levels
         688             :        !----------------------------------------------------------------
         689           0 :        IF ( nz == 47 ) THEN
         690             :
         691             :           ! Determine if the variable is on model levels or edges
         692             :           IF ( nlev == 72 ) THEN
         693           0 :              NL   = 36
         694           0 :              nout = 47
         695             :           ELSEIF ( nlev == 73 ) THEN
         696           0 :              NL   = 37
         697           0 :              nout = 48
         698             :           ELSEIF ( nlev == 3 ) THEN
         699           0 :              NL   = 3
         700           0 :              nout = 47
         701             :           ELSEIF ( nlev == 11 ) THEN
         702           0 :              NL   = 11
         703           0 :              nout = 47
         704             :           ELSEIF ( nlev == 36 ) THEN
         705           0 :              NL   = 36
         706           0 :              nout = 47
         707             :           ELSE
         708             :              MSG = 'Can only remap from native onto reduced GEOS-5 if '// &
         709           0 :                    'input data has exactly 72, 73, 3, 11, or 36 levels: '//TRIM(Lct%Dct%cName)
         710           0 :              CALL HCO_ERROR( MSG, RC )
         711           0 :              RETURN
         712             :           ENDIF ! nlev == (72,73,3,11,36ELSE)
         713             :
         714             :           ! Make sure output array is allocated
         715           0 :           CALL FileData_ArrCheck( HcoState%Config, Lct%Dct%Dta, nx, ny, nout, nt, RC )
         716             :
         717             :           ! Do for every time slice
         718           0 :           DO T = 1, nt
         719             :
         720             :              ! Levels that are passed level-by-level.
         721           0 :              DO L = 1, NL
         722           0 :                 Lct%Dct%Dta%V3(T)%Val(:,:,L) = REGR_4D(:,:,L,T)
         723             :              ENDDO !L
         724             :
         725             :             ! If remapping model grid layers, collapse layers
         726           0 :              IF ( nlev == 72 ) THEN
         727             :                ! Collapse two levels (e.g. levels 37-38 into level 37):
         728           0 :                CALL COLLAPSE( Lct, REGR_4D, 37, 37, 2, T, 5, RC )
         729           0 :                CALL COLLAPSE( Lct, REGR_4D, 38, 39, 2, T, 5, RC )
         730           0 :                CALL COLLAPSE( Lct, REGR_4D, 39, 41, 2, T, 5, RC )
         731           0 :                CALL COLLAPSE( Lct, REGR_4D, 40, 43, 2, T, 5, RC )
         732             :                ! Collapse four levels:
         733           0 :                CALL COLLAPSE( Lct, REGR_4D, 41, 45, 4, T, 5, RC )
         734           0 :                CALL COLLAPSE( Lct, REGR_4D, 42, 49, 4, T, 5, RC )
         735           0 :                CALL COLLAPSE( Lct, REGR_4D, 43, 53, 4, T, 5, RC )
         736           0 :                CALL COLLAPSE( Lct, REGR_4D, 44, 57, 4, T, 5, RC )
         737           0 :                CALL COLLAPSE( Lct, REGR_4D, 45, 61, 4, T, 5, RC )
         738           0 :                CALL COLLAPSE( Lct, REGR_4D, 46, 65, 4, T, 5, RC )
         739           0 :                CALL COLLAPSE( Lct, REGR_4D, 47, 69, 4, T, 5, RC )
         740             :              ! If remapping model grid edges, sample at edges
         741           0 :              ELSEIF ( nlev == 73 ) THEN
         742           0 :                Lct%Dct%Dta%V3(T)%Val(:,:,38) = REGR_4D(:,:,39,T)
         743           0 :                Lct%Dct%Dta%V3(T)%Val(:,:,39) = REGR_4D(:,:,41,T)
         744           0 :                Lct%Dct%Dta%V3(T)%Val(:,:,40) = REGR_4D(:,:,43,T)
         745           0 :                Lct%Dct%Dta%V3(T)%Val(:,:,41) = REGR_4D(:,:,45,T)
         746           0 :                Lct%Dct%Dta%V3(T)%Val(:,:,42) = REGR_4D(:,:,49,T)
         747           0 :                Lct%Dct%Dta%V3(T)%Val(:,:,43) = REGR_4D(:,:,53,T)
         748           0 :                Lct%Dct%Dta%V3(T)%Val(:,:,44) = REGR_4D(:,:,57,T)
         749           0 :                Lct%Dct%Dta%V3(T)%Val(:,:,45) = REGR_4D(:,:,61,T)
         750           0 :                Lct%Dct%Dta%V3(T)%Val(:,:,46) = REGR_4D(:,:,65,T)
         751           0 :                Lct%Dct%Dta%V3(T)%Val(:,:,47) = REGR_4D(:,:,69,T)
         752           0 :                Lct%Dct%Dta%V3(T)%Val(:,:,48) = REGR_4D(:,:,73,T)
         753             :             ! If the input is N (N = 3/11/36) levels, levels N+1 to 47 are set to 0
         754           0 :              ELSEIF ( nlev == 3 ) THEN
         755           0 :                DO L = 4,47
         756           0 :                   Lct%Dct%Dta%V3(T)%Val(:,:,L) = 0.0_hp
         757             :                ENDDO !L
         758           0 :              ELSEIF ( nlev == 11 ) THEN
         759           0 :                DO L = 12,47
         760           0 :                   Lct%Dct%Dta%V3(T)%Val(:,:,L) = 0.0_hp
         761             :                ENDDO !L
         762           0 :              ELSEIF ( nlev == 36 ) THEN
         763           0 :                DO L = 37,47
         764           0 :                   Lct%Dct%Dta%V3(T)%Val(:,:,L) = 0.0_hp
         765             :                ENDDO !L
         766             :              ENDIF ! nlev == (72,73,36)
         767             :
         768             :           ENDDO ! T
         769             :
         770             :           ! Verbose
         771           0 :           IF ( HCO_IsVerb( HcoState%Config%Err ) ) THEN
         772           0 :              WRITE(MSG,*) 'Mapped ', nlev, ' levels onto reduced GEOS-5 levels.'
         773           0 :              CALL HCO_MSG(HcoState%Config%Err,MSG)
         774           0 :              IF ( nlev == 3 .OR. nlev == 11 .OR. nlev == 36 ) THEN
         775           0 :                WRITE(MSG,*) 'The input variable has ', nlev, 'L, which were written to be L 1-', nlev, ' on the output 47 L grid (remaining values set to 0).'
         776             :              ELSE
         777           0 :                WRITE(MSG,*) 'Pressure-weighted vertical regridding was done - consider if this is appropriate for the variable units.'
         778           0 :              CALL HCO_MSG(HcoState%Config%Err,MSG)
         779             :              ENDIF
         780             :           ENDIF
         781             :
         782             :           ! Done!
         783             :           DONE = .TRUE.
         784             :
         785             :        !----------------------------------------------------------------
         786             :        ! Native to reduced GISS levels
         787             :        !----------------------------------------------------------------
         788           0 :        ELSEIF ( nz == 74 ) THEN
         789             :
         790             :           ! Determine if the variable is on model levels or edges
         791           0 :           IF ( nlev == 102 ) THEN
         792           0 :              NL   = 60
         793           0 :              nout = 74
         794           0 :           ELSEIF ( nlev == 103 ) THEN
         795           0 :              NL   = 61
         796           0 :              nout = 75
         797             :           ELSE
         798             :              MSG = 'Can only remap from native onto reduced GISS if '// &
         799           0 :                    'input data has exactly 102 or 103 levels: '//TRIM(Lct%Dct%cName)
         800           0 :              CALL HCO_ERROR( MSG, RC )
         801           0 :              RETURN
         802             :           ENDIF ! nlev == (102,103,ELSE)
         803             :
         804             :           ! Make sure output array is allocated
         805           0 :           CALL FileData_ArrCheck( HcoState%Config, Lct%Dct%Dta, nx, ny, nout, nt, RC )
         806             :
         807             :           ! Do for every time slice
         808           0 :           DO T = 1, nt
         809             :
         810             :              ! Levels that are passed level-by-level.
         811           0 :              DO L = 1, NL
         812           0 :                 Lct%Dct%Dta%V3(T)%Val(:,:,L) = REGR_4D(:,:,L,T)
         813             :              ENDDO !L
         814             :
         815             :              ! If remapping model grid layers, collapse layers
         816           0 :              IF ( nlev == 102 ) THEN
         817             :                 ! Collapse two levels (e.g. levels 61-62 into level 61):
         818           0 :                 CALL COLLAPSE( Lct, REGR_4D, 61, 61, 2, T, 22, RC )
         819           0 :                 CALL COLLAPSE( Lct, REGR_4D, 62, 63, 2, T, 22, RC )
         820           0 :                 CALL COLLAPSE( Lct, REGR_4D, 63, 65, 2, T, 22, RC )
         821           0 :                 CALL COLLAPSE( Lct, REGR_4D, 64, 67, 2, T, 22, RC )
         822           0 :                 CALL COLLAPSE( Lct, REGR_4D, 65, 69, 2, T, 22, RC )
         823           0 :                 CALL COLLAPSE( Lct, REGR_4D, 66, 71, 2, T, 22, RC )
         824           0 :                 CALL COLLAPSE( Lct, REGR_4D, 67, 73, 2, T, 22, RC )
         825             :                 ! Collapse four levels:
         826           0 :                 CALL COLLAPSE( Lct, REGR_4D, 68, 75, 4, T, 22, RC )
         827           0 :                 CALL COLLAPSE( Lct, REGR_4D, 69, 79, 4, T, 22, RC )
         828           0 :                 CALL COLLAPSE( Lct, REGR_4D, 70, 83, 4, T, 22, RC )
         829           0 :                 CALL COLLAPSE( Lct, REGR_4D, 71, 87, 4, T, 22, RC )
         830           0 :                 CALL COLLAPSE( Lct, REGR_4D, 72, 91, 4, T, 22, RC )
         831           0 :                 CALL COLLAPSE( Lct, REGR_4D, 73, 95, 4, T, 22, RC )
         832           0 :                 CALL COLLAPSE( Lct, REGR_4D, 74, 99, 4, T, 22, RC )
         833             :              ! If remapping model grid edges, sample at the edges
         834             :              ELSE
         835           0 :                 Lct%Dct%Dta%V3(T)%Val(:,:,62) = REGR_4D(:,:,63,T)
         836           0 :                 Lct%Dct%Dta%V3(T)%Val(:,:,63) = REGR_4D(:,:,65,T)
         837           0 :                 Lct%Dct%Dta%V3(T)%Val(:,:,64) = REGR_4D(:,:,67,T)
         838           0 :                 Lct%Dct%Dta%V3(T)%Val(:,:,65) = REGR_4D(:,:,69,T)
         839           0 :                 Lct%Dct%Dta%V3(T)%Val(:,:,66) = REGR_4D(:,:,71,T)
         840           0 :                 Lct%Dct%Dta%V3(T)%Val(:,:,67) = REGR_4D(:,:,73,T)
         841           0 :                 Lct%Dct%Dta%V3(T)%Val(:,:,68) = REGR_4D(:,:,75,T)
         842           0 :                 Lct%Dct%Dta%V3(T)%Val(:,:,69) = REGR_4D(:,:,79,T)
         843           0 :                 Lct%Dct%Dta%V3(T)%Val(:,:,70) = REGR_4D(:,:,83,T)
         844           0 :                 Lct%Dct%Dta%V3(T)%Val(:,:,71) = REGR_4D(:,:,87,T)
         845           0 :                 Lct%Dct%Dta%V3(T)%Val(:,:,72) = REGR_4D(:,:,91,T)
         846           0 :                 Lct%Dct%Dta%V3(T)%Val(:,:,73) = REGR_4D(:,:,95,T)
         847           0 :                 Lct%Dct%Dta%V3(T)%Val(:,:,74) = REGR_4D(:,:,99,T)
         848           0 :                 Lct%Dct%Dta%V3(T)%Val(:,:,75) = REGR_4D(:,:,103,T)
         849             :              ENDIF ! nlev == (102,ELSE)
         850             :           ENDDO ! T
         851             :
         852             :           ! Verbose
         853           0 :           IF ( HCO_IsVerb(HcoState%Config%Err ) ) THEN
         854           0 :              WRITE(MSG,*) 'Mapped ', nlev, ' levels onto reduced GISS levels.'
         855           0 :              CALL HCO_MSG(HcoState%Config%Err,MSG)
         856           0 :              WRITE(MSG,*) 'Pressure-weighted vertical regridding was done - consider if this is appropriate for the variable units.'
         857           0 :              CALL HCO_MSG(HcoState%Config%Err,MSG)
         858             :           ENDIF
         859             :
         860             :           ! Done!
         861             :           DONE = .TRUE.
         862             :
         863             :        !----------------------------------------------------------------
         864             :        ! Reduced to native GEOS-5 levels
         865             :        !----------------------------------------------------------------
         866           0 :        ELSEIF ( nz == 72 ) THEN
         867             :
         868             :           ! Determine if the variable is on model levels or edges
         869             :           IF ( nlev == 47 ) THEN
         870           0 :              NL   = 36
         871           0 :              nout = 72
         872             :           ELSEIF ( nlev == 48 ) THEN
         873           0 :              NL   = 37
         874           0 :              nout = 73
         875             :           ELSEIF ( nlev == 3 ) THEN
         876           0 :              NL   = 3
         877           0 :              nout = 72
         878             :           ELSEIF ( nlev == 11 ) THEN
         879           0 :              NL   = 11
         880           0 :              nout = 72
         881             :           ELSEIF ( nlev == 36 ) THEN
         882           0 :              NL   = 36
         883           0 :              nout = 72
         884             :           ELSE
         885             :              MSG = 'Can only remap from reduced onto native GEOS-5 if '// &
         886           0 :                    'input data has exactly 47, 48, 3, 11, or 36 levels: '//TRIM(Lct%Dct%cName)
         887           0 :              CALL HCO_ERROR( MSG, RC )
         888           0 :              RETURN
         889             :           ENDIF ! nlev == (48,48,3,11,36,ELSE)
         890             :
         891             :           ! Make sure output array is allocated
         892           0 :           CALL FileData_ArrCheck( HcoState%Config, Lct%Dct%Dta, nx, ny, nout, nt, RC )
         893             :
         894             :           ! Do for every time slice
         895           0 :           DO T = 1, nt
         896             :
         897             :              ! Levels that are passed level-by-level.
         898           0 :              DO L = 1, NL
         899           0 :                 Lct%Dct%Dta%V3(T)%Val(:,:,L) = REGR_4D(:,:,L,T)
         900             :              ENDDO !L
         901             :
         902             :              ! If remapping model grid layers, inflate layers
         903           0 :              IF ( nlev == 47 ) THEN
         904             :
         905             :                 ! Inflate two levels (e.g. levels 37-38 on the fine grid are copies of level 37 on the coarse grid):
         906             :                 coarseIDX = 36
         907           0 :                 DO I = 1,8
         908           0 :                    fineIDX = 36 + I
         909           0 :                    IF ( MOD(I,2) /= 0 ) THEN
         910           0 :                      coarseIDX = coarseIDX + 1
         911             :                    ENDIF
         912           0 :                    Lct%Dct%Dta%V3(T)%Val(:,:,fineIDX) = REGR_4D(:,:,coarseIDX,T)
         913             :                 ENDDO ! I
         914             :
         915             :                 ! Inflate four levels (e.g. levels 45-48 on the fine grid are copies of level 41 on the coarse grid)
         916             :                 coarseIDX = 40
         917           0 :                 DO I = 1,28
         918           0 :                    fineIDX = 44 + i
         919           0 :                    IF ( MOD(I-1,4) == 0 ) THEN
         920           0 :                       coarseIDX = coarseIDX + 1
         921             :                    ENDIF
         922           0 :                    Lct%Dct%Dta%V3(T)%Val(:,:,fineIDX) = REGR_4D(:,:,coarseIDX,T)
         923             :                 ENDDO ! I
         924             :
         925             :              ! If remapping model grid edges, inflate edges
         926           0 :              ELSEIF ( nlev == 48 ) THEN
         927             :
         928             :                 ! Sample every two edges (e.g. edges 38-39 on the fine grid are copies of edge 38 on the coarse grid)
         929             :                 coarseIDX = 37
         930           0 :                 DO I = 1,8
         931           0 :                    fineIDX = 37 + I
         932           0 :                    IF ( MOD(I,2) /= 0 ) THEN
         933           0 :                      coarseIDX = coarseIDX + 1
         934             :                    ENDIF
         935           0 :                    Lct%Dct%Dta%V3(T)%Val(:,:,fineIDX) = REGR_4D(:,:,coarseIDX,T)
         936             :                 ENDDO ! I
         937             :
         938             :                 ! Sample every four edges (e.g. edges 46-49 on the fine grid are copies of edge 42 on the coarse grid)
         939             :                 coarseIDX = 40
         940           0 :                 DO I = 1,28
         941           0 :                    fineIDX = 44 + i
         942           0 :                    IF ( MOD(I-1,4) == 0 ) THEN
         943           0 :                       coarseIDX = coarseIDX + 1
         944             :                    ENDIF
         945           0 :                    Lct%Dct%Dta%V3(T)%Val(:,:,fineIDX) = REGR_4D(:,:,coarseIDX,T)
         946             :                 ENDDO ! I
         947             :
         948           0 :              ELSEIF ( nlev == 3 ) THEN
         949           0 :                DO L = 4,72
         950           0 :                   Lct%Dct%Dta%V3(T)%Val(:,:,L) = 0.0_hp
         951             :                ENDDO !L
         952             :
         953           0 :              ELSEIF ( nlev == 11 ) THEN
         954           0 :                DO L = 12,72
         955           0 :                   Lct%Dct%Dta%V3(T)%Val(:,:,L) = 0.0_hp
         956             :                ENDDO !L
         957             :
         958           0 :              ELSEIF ( nlev == 36 ) THEN
         959           0 :                DO L = 37,72
         960           0 :                   Lct%Dct%Dta%V3(T)%Val(:,:,L) = 0.0_hp
         961             :                ENDDO !L
         962             :
         963             :              ENDIF ! nlev == (47,48,3,11,36)
         964             :           ENDDO ! T
         965             :
         966             :           ! Verbose
         967           0 :           IF ( HCO_IsVerb(HcoState%Config%Err ) ) THEN
         968           0 :              WRITE(MSG,*) 'Mapped ', nlev, ' levels onto native GEOS-5 levels.'
         969           0 :              CALL HCO_MSG(HcoState%Config%Err,MSG)
         970           0 :              IF ( nlev == 3 .OR. nlev == 11 .OR. nlev == 36 ) THEN
         971           0 :                WRITE(MSG,*) 'The input variable has ', nlev, 'L, which were written to be L 1-', nlev, ' on the output 72 L grid (remaining values set to 0).'
         972             :              ELSE
         973           0 :                WRITE(MSG,*) 'Inflating from 47/48 to 72/73 levels is not recommended and is likely not mass-conserving.'
         974             :              ENDIF
         975           0 :              CALL HCO_MSG(HcoState%Config%Err,MSG)
         976             :           ENDIF
         977             :
         978             :           ! Done!
         979             :           DONE = .TRUE.
         980             :
         981             :        ENDIF ! nz == (47,74,72)
         982             :
         983             :     ENDIF ! Vertical regridding required
         984             :
         985             :     !===================================================================
         986             :     ! Error check / verbose mode
         987             :     !===================================================================
         988             :     IF ( DONE ) THEN
         989           0 :       IF ( HCO_IsVerb(HcoState%Config%Err ) ) THEN
         990           0 :           WRITE(MSG,*) 'Did vertical regridding for ',TRIM(Lct%Dct%cName),':'
         991           0 :           CALL HCO_MSG(HcoState%Config%Err,MSG)
         992           0 :           WRITE(MSG,*) 'Number of original levels: ', nlev
         993           0 :           CALL HCO_MSG(HcoState%Config%Err,MSG)
         994           0 :           WRITE(MSG,*) 'Number of output levels: ', SIZE(Lct%Dct%Dta%V3(1)%Val,3)
         995           0 :           CALL HCO_MSG(HcoState%Config%Err,MSG)
         996             :        ENDIF
         997             :     ELSE
         998           0 :        WRITE(MSG,*) 'Vertical regridding failed: ',TRIM(Lct%Dct%cName)
         999           0 :        CALL HCO_ERROR( MSG, RC )
        1000           0 :        RETURN
        1001             :     ENDIF
        1002             :
        1003             :     ! Return w/ success
        1004           0 :     CALL HCO_LEAVE ( HcoState%Config%Err, RC )
        1005             :
        1006             :   END SUBROUTINE ModelLev_Interpolate
        1007             : !EOC
        1008             : !------------------------------------------------------------------------------
        1009             : !                   Harmonized Emissions Component (HEMCO)                    !
        1010             : !------------------------------------------------------------------------------
        1011             : !BOP
        1012             : !
        1013             : ! !IROUTINE: COLLAPSE
        1014             : !
        1015             : ! !DESCRIPTION: Helper routine to collapse input levels onto the output grid.
        1016             : ! The input data is weighted by the grid box thicknesses defined on top of
        1017             : ! this module. The input parameter T determines the time slice to be considered,
        1018             : ! and MET denotes the met field type of the input data (22 = GISS, 5= GEOS-5).
        1019             : !\\
        1020             : !\\
        1021             : ! !INTERFACE:
        1022             : !
        1023           0 :   SUBROUTINE COLLAPSE ( Lct, REGR_4D, OutLev, InLev1, NLEV, T, MET, RC )
        1024             : !
        1025             : ! !INPUT PARAMETERS:
        1026             : !
        1027             :     REAL(sp),         POINTER        :: REGR_4D(:,:,:,:)  ! 4D input data
        1028             :     INTEGER,          INTENT(IN)     :: OutLev
        1029             :     INTEGER,          INTENT(IN)     :: InLev1
        1030             :     INTEGER,          INTENT(IN)     :: NLEV
        1031             :     INTEGER,          INTENT(IN)     :: T
        1032             :     INTEGER,          INTENT(IN)     :: MET               ! 22=GISS E2.2, else GEOS-5
        1033             : !
        1034             : ! !INPUT/OUTPUT PARAMETERS:
        1035             : !
        1036             :     TYPE(ListCont),   POINTER        :: Lct               ! HEMCO list container
        1037             :     INTEGER,          INTENT(INOUT)  :: RC                ! Success or failure?
        1038             : !
        1039             : ! !REVISION HISTORY:
        1040             : !  30 Dec 2014 - C. Keller   - Initial version
        1041             : !  See https://github.com/geoschem/hemco for complete history
        1042             : !EOP
        1043             : !------------------------------------------------------------------------------
        1044             : !BOC
        1045             :     INTEGER               :: I, NZ, ILEV, TOPLEV
        1046             :     REAL(hp)              :: THICK
        1047           0 :     REAL(hp), POINTER     :: EDG(:)
        1048           0 :     REAL(hp), ALLOCATABLE :: WGT(:)
        1049             :     CHARACTER(LEN=255)      :: MSG
        1050             :     CHARACTER(LEN=255)      :: LOC = 'ModelLev_Interpolate (hco_interp_mod.F90)'
        1051             :
        1052             :     !=================================================================
        1053             :     ! COLLAPSE begins here
        1054             :     !=================================================================
        1055             :
        1056             :     ! Init
        1057           0 :     EDG => NULL()
        1058             :
        1059             :     ! Reset
        1060           0 :     Lct%Dct%Dta%V3(T)%Val(:,:,OutLev) = 0.0_hp
        1061             :
        1062             :     ! Don't do anything if there are not enough levels in REGR_4D
        1063           0 :     NZ = SIZE(REGR_4D,3)
        1064           0 :     IF ( NZ < InLev1 ) RETURN
        1065             :
        1066             :     ! Get maximum level to be used for pressure thickness calculations.
        1067           0 :     TOPLEV = InLev1 + NLEV
        1068             :
        1069             :     ! Get pointer to grid edges on the native input grid
        1070           0 :     IF ( Met == 22 ) THEN
        1071           0 :        EDG => E102_EDGE_NATIVE(InLev1:TOPLEV)
        1072           0 :     ELSEIF ( Met == 5 ) THEN
        1073           0 :        EDG => G5_EDGE_NATIVE(InLev1:TOPLEV)
        1074             :     ELSE
        1075           0 :        WRITE(MSG,*) 'The Met value given was not valid: ', Met
        1076           0 :        CALL HCO_ERROR( MSG, RC, THISLOC=LOC )
        1077           0 :        RETURN
        1078             :     ENDIF
        1079             :
        1080             :     ! Thickness of output level
        1081           0 :     THICK = EDG(1) - EDG(NLEV+1)
        1082             :
        1083             :     ! Get level weights
        1084           0 :     ALLOCATE(WGT(NLEV))
        1085           0 :     WGT = 0.0
        1086           0 :     DO I = 1, NLEV
        1087           0 :        WGT(I) = ( EDG(I) - EDG(I+1) ) / THICK
        1088             :     ENDDO
        1089             :
        1090             :     ! Pass levels to output data, one after each other
        1091           0 :     Lct%Dct%Dta%V3(T)%Val(:,:,OutLev) = REGR_4D(:,:,InLev1,T) * WGT(1)
        1092           0 :     DO I = 1, NLEV-1
        1093           0 :        ILEV = InLev1 + I
        1094           0 :        IF ( NZ < ILEV ) EXIT
        1095           0 :        Lct%Dct%Dta%V3(T)%Val(:,:,OutLev) = Lct%Dct%Dta%V3(T)%Val(:,:,OutLev) &
        1096           0 :                                          + ( REGR_4D(:,:,ILEV,T) * WGT(I+1) )
        1097             :     ENDDO
        1098             :
        1099             :     ! Cleanup
        1100           0 :     DEALLOCATE(WGT)
        1101           0 :     EDG => NULL()
        1102             :
        1103           0 :   END SUBROUTINE COLLAPSE
        1104             : !EOC
        1105             : END MODULE HCO_Interp_Mod

Generated by: LCOV version 1.14