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

              Line data    Source code
           1             : !------------------------------------------------------------------------------
           2             : !                   Harmonized Emissions Component (HEMCO)                    !
           3             : !------------------------------------------------------------------------------
           4             : !BOP
           5             : !
           6             : ! !MODULE: hcoio_messy_mod.F90
           7             : !
           8             : ! !DESCRIPTION: Module HCOIO\_MESSY\_MOD interfaces HEMCO with the regridding
           9             : ! tool NCREGRID of the Modular Earth Submodel System (MESSy). This regridding
          10             : ! scheme is used for vertical regridding and/or for index data, i.e. data with
          11             : ! discrete values (e.g. land type integers). This code currently only works for
          12             : ! rectilinear (regular lon-lat) grids but can be extended to support curvilinear
          13             : ! grids.
          14             : !\\
          15             : ! REFERENCES:
          16             : ! \begin{itemize}
          17             : ! \item Joeckel, P. Technical note: Recursive rediscretisation of geo-
          18             : ! scientific data in the Modular Earth Submodel System (MESSy), ACP, 6,
          19             : ! 3557--3562, 2006.
          20             : ! \end{itemize}
          21             : ! !INTERFACE:
          22             : !
          23             : MODULE HCOIO_MESSY_MOD
          24             : !
          25             : ! !USES:
          26             : !
          27             :   USE HCO_ERROR_MOD
          28             :   USE HCO_TYPES_MOD,        ONLY : ListCont
          29             :   USE HCO_STATE_MOD,        ONLY : Hco_State
          30             :   USE MESSY_NCREGRID_BASE,  ONLY : NARRAY, AXIS
          31             :
          32             :   IMPLICIT NONE
          33             :   PRIVATE
          34             : !
          35             : ! !PUBLIC MEMBER FUNCTIONS:
          36             : !
          37             :   PUBLIC  :: HCO_MESSY_REGRID
          38             : !
          39             : ! !PRIVATE MEMBER FUNCTIONS:
          40             : !
          41             :   ! Set this value to TRUE if you want to reduce the output array
          42             :   ! to the minimum required number of vertical levels.
          43             :   LOGICAL, PARAMETER            :: ReduceVert = .FALSE.
          44             : !
          45             : ! !MODULE INTERFACES:
          46             : !
          47             : ! !REVISION HISTORY:
          48             : !  24 Jun 2014 - C. Keller   - Initial version
          49             : !  See https://github.com/geoschem/hemco for complete history
          50             : !EOP
          51             : !------------------------------------------------------------------------------
          52             : !BOC
          53             : !
          54             : ! !MODULE VARIABLES:
          55             :
          56             :   CONTAINS
          57             : !EOC
          58             : !------------------------------------------------------------------------------
          59             : !                   Harmonized Emissions Component (HEMCO)                    !
          60             : !------------------------------------------------------------------------------
          61             : !BOP
          62             : !
          63             : ! !IROUTINE: Hco_Messy_Regrid
          64             : !
          65             : ! !DESCRIPTION: This is the wrapper routine to regrid a 4D input array
          66             : ! NcArr (x,y,z,t) onto the HEMCO emissions grid (defined in HcoState)
          67             : ! using the regridding tool NCREGRID. LonEdge, LatEdge and LevEdge are
          68             : ! the grid point edges of the input grid. The data is written into list
          69             : ! container Lct.
          70             : !\\
          71             : !\\
          72             : ! If the input grid is 2D (horizontal only), LevEdge must not be specified
          73             : ! (null pointer) and the data is regridded in the horizontal only. If the
          74             : ! input grid has only one vertical level, it is assumed that this is the
          75             : ! surface level and the output data is 3D but with only one vertical level.
          76             : !\\
          77             : !\\
          78             : ! For input data with more than one vertical level, the data is mapped
          79             : ! onto the entire 3D grid. The module parameter ReduceVert can be used to
          80             : ! cap the output data at the lowest possible level. For example, if the
          81             : ! input grid only covers three surface levels with a minimum sigma value
          82             : ! of 0.75, vertical regridding is performed within this sigma range (1-0.75)
          83             : ! and the output grid is reduced accordingly. This option is not used in the
          84             : ! standard HEMCO setup because problems can arise if the data array of a
          85             : ! given container suddently changes its size (i.e. when updated data covers
          86             : ! more/less vertical levels than the data beforehand).
          87             : !\\
          88             : !\\
          89             : ! The input argument IsModelLev denotes whether or not the vertical
          90             : ! coordinates of the input data are on model levels. If set to yes and LevEdge
          91             : ! is not provided (i.e. a nullified pointer), the MESSy regridding routines are
          92             : ! only used for the horizontal remapping and subroutine ModelLev\_Interpolate
          93             : ! (module hco\_interp\_mod.F90) is used for the vertical remapping.
          94             : !\\
          95             : !\\
          96             : ! !INTERFACE:
          97             : !
          98           0 :   SUBROUTINE HCO_MESSY_REGRID ( HcoState,  NcArr,               &
          99             :                                 LonEdge,   LatEdge,    LevEdge, &
         100             :                                 Lct,       IsModelLev, RC        )
         101             : !
         102             : ! !USES:
         103             : !
         104             :   USE HCO_FILEDATA_MOD,     ONLY : FileData_ArrCheck
         105             :   USE HCO_UNIT_MOD,         ONLY : HCO_IsIndexData
         106             :   USE HCO_INTERP_MOD,       ONLY : ModelLev_Interpolate
         107             :   USE MESSY_NCREGRID_BASE,  ONLY : RG_INT, RG_IDX
         108             :   USE MESSY_NCREGRID_BASE,  ONLY : NREGRID
         109             :   USE MESSY_NCREGRID_BASE,  ONLY : INIT_NARRAY
         110             : !
         111             : ! !INPUT/OUTPUT PARAMETERS:
         112             : !
         113             :     TYPE(HCO_State),  POINTER        :: HcoState        ! HEMCO obj.
         114             :     REAL(sp),         POINTER        :: ncArr(:,:,:,:)  ! Input array(x,y,z,t)
         115             :     REAL(hp),         POINTER        :: LonEdge(:)      ! lon edges
         116             :     REAL(hp),         POINTER        :: LatEdge(:)      ! lat edges
         117             :     REAL(hp),         POINTER        :: LevEdge(:,:,:)  ! sigma level edges
         118             :     TYPE(ListCont),   POINTER        :: Lct             ! Target list container
         119             :     LOGICAL,          INTENT(IN   )  :: IsModelLev      ! Are these model levels?
         120             :     INTEGER,          INTENT(INOUT)  :: RC              ! Return code
         121             : !
         122             : ! !REVISION HISTORY:
         123             : !  27 Jun 2014 - C. Keller - Initial version
         124             : !  See https://github.com/geoschem/hemco for complete history
         125             : !EOP
         126             : !------------------------------------------------------------------------------
         127             : !BOC
         128             : !
         129             : ! !ROUTINE ARGUMENTS:
         130             : !
         131           0 :     TYPE(narray), POINTER         :: narr_src(:)
         132           0 :     TYPE(narray), POINTER         :: narr_dst(:)
         133           0 :     TYPE(axis),   POINTER         :: axis_src(:)
         134           0 :     TYPE(axis),   POINTER         :: axis_dst(:)
         135           0 :     INTEGER,      POINTER         :: rg_type(:)
         136           0 :     INTEGER,      POINTER         :: rcnt   (:)
         137           0 :     REAL(dp),     POINTER         :: sovl   (:)
         138           0 :     REAL(dp),     POINTER         :: dovl   (:)
         139           0 :     REAL(hp),     POINTER         :: lon    (:)
         140           0 :     REAL(hp),     POINTER         :: lat    (:)
         141           0 :     REAL(hp),     POINTER         :: sigma  (:,:,:)
         142           0 :     REAL(sp),     POINTER         :: ArrIn  (:,:,:,:)
         143           0 :     REAL(sp),     POINTER         :: ArrOut (:,:,:,:)
         144           0 :     REAL(hp), ALLOCATABLE, TARGET :: sigout (:,:,:)
         145             :     REAL(hp)                      :: sigMin
         146             :     INTEGER                       :: NZIN, NZOUT, NTIME
         147             :     INTEGER                       :: NXIN, NYIN
         148             :     INTEGER                       :: I, L, AS
         149             :     INTEGER                       :: NCALLS
         150             :     CHARACTER(LEN=255)            :: MSG, LOC
         151             :     LOGICAL                       :: SameGrid, verb
         152             :
         153             :     !=================================================================
         154             :     ! HCO_MESSY_REGRID begins here
         155             :     !=================================================================
         156             :
         157             :     ! For error handling
         158           0 :     LOC = 'HCO_MESSY_REGRID (HCOIO_MESSY_MOD.F90)'
         159           0 :     CALL HCO_ENTER ( HcoState%Config%Err, LOC, RC )
         160           0 :     IF ( RC /= HCO_SUCCESS ) THEN
         161           0 :         CALL HCO_ERROR( 'ERROR 0', RC, THISLOC=LOC )
         162           0 :         RETURN
         163             :     ENDIF
         164             :
         165             :     ! Init
         166           0 :     narr_src => NULL()
         167           0 :     narr_dst => NULL()
         168           0 :     axis_src => NULL()
         169           0 :     axis_dst => NULL()
         170           0 :     rg_type  => NULL()
         171           0 :     rcnt     => NULL()
         172           0 :     sovl     => NULL()
         173           0 :     dovl     => NULL()
         174           0 :     lon      => NULL()
         175           0 :     lat      => NULL()
         176           0 :     sigma    => NULL()
         177           0 :     ArrIn    => NULL()
         178           0 :     ArrOut   => NULL()
         179             :
         180             :     ! verbose?
         181           0 :     verb = HCO_IsVerb( HcoState%Config%Err )
         182             :
         183             :     ! Horizontal dimension of input data
         184           0 :     NXIN = SIZE(NcArr,1)
         185           0 :     NYIN = SIZE(NcArr,2)
         186             :
         187             :     ! Number of vertical levels of input data
         188           0 :     NZIN = SIZE(NcArr,3)
         189             :
         190             :     ! Number of time slices. All time slices will be regridded
         191             :     ! simultaneously.
         192           0 :     NTIME = SIZE(NcArr,4)
         193             :
         194             :     ! Error check: data must be 2D or 3D.
         195           0 :     IF ( Lct%Dct%Dta%SpaceDim /= 2 .AND. &
         196             :          Lct%Dct%Dta%SpaceDim /= 3         ) THEN
         197           0 :        MSG = 'Can only regrid 2D or 3D data: ' // TRIM(Lct%Dct%cName)
         198           0 :        CALL HCO_ERROR ( MSG, RC )
         199           0 :        RETURN
         200             :     ENDIF
         201             :
         202             :     !-----------------------------------------------------------------
         203             :     ! Shortcut if input field is already on output grid: directly
         204             :     ! pass data to Lct.
         205             :     ! NOTE: if the number of input levels matches the number of output
         206             :     ! levels (and the horizontal dimensions agree as well), it is
         207             :     ! assumed that they are on the same grid!
         208             :     !-----------------------------------------------------------------
         209             :
         210             :     ! Input grid = output grid?
         211           0 :     SameGrid = .FALSE.
         212             :
         213             :     ! Horizontal dimensions have to match
         214           0 :     IF ( (NXIN == HcoState%NX) .AND. (NYIN == HcoState%NY) ) THEN
         215             :
         216             :        ! Vertical dimensions have to match or be 1
         217           0 :        IF ( NZIN == 1 .OR. NZIN == HcoState%NZ ) THEN
         218             :
         219             :           ! Assume same grid
         220           0 :           SameGrid = .TRUE.
         221             :
         222             :           ! Check for same boundaries. Otherwise falsify SameGrid
         223             :           ! Lon ...
         224           0 :           IF ( MINVAL(LonEdge) /= MINVAL(HcoState%Grid%XEDGE%Val) .OR. &
         225             :                MAXVAL(LonEdge) /= MAXVAL(HcoState%Grid%XEDGE%Val) ) THEN
         226           0 :              SameGrid = .FALSE.
         227             :           ENDIF
         228             :           ! ... Lat ...
         229           0 :           IF ( MINVAL(LatEdge) /= MINVAL(HcoState%Grid%YEDGE%Val) .OR. &
         230             :                MAXVAL(LatEdge) /= MAXVAL(HcoState%Grid%YEDGE%Val) ) THEN
         231           0 :              SameGrid = .FALSE.
         232             :           ENDIF
         233             :           ! ... Lev
         234             :           IF ( NZIN > 1 ) THEN
         235             :              ! TODO: Eventually need to add level check here.
         236             :              ! For now, assume that input levels are equal to
         237             :              ! output level if they have the same dimension!
         238             :           ENDIF
         239             :        ENDIF
         240             :     ENDIF
         241             :
         242             :     !-----------------------------------------------------------------
         243             :     ! Define number of vertical levels on output grid
         244             :     !-----------------------------------------------------------------
         245             :
         246             :     ! Error check. If we are not on the same grid, LevEdge must be
         247             :     ! provided or IsModelLev must be set to true for 3D data.
         248           0 :     IF ( NZIN > 1 .AND. .NOT. SameGrid ) THEN
         249           0 :        IF ( .NOT. ASSOCIATED(LevEdge) .AND. .NOT. IsModelLev ) THEN
         250             :           MSG = 'Cannot regrid '//TRIM(Lct%Dct%cName)//'. Either level '//&
         251           0 :                 'edges must be provided or data must be on model levels.'
         252           0 :           CALL HCO_ERROR ( MSG, RC )
         253           0 :           RETURN
         254             :        ENDIF
         255             :     ENDIF
         256             :
         257           0 :     NZOUT = NZIN
         258           0 :     IF ( .NOT. SameGrid .AND. ASSOCIATED(LevEdge) ) THEN
         259             :
         260             :        ! Calculate sigma level for each grid point on the output grid.
         261             :        ! This is the pressure at location i,j,l normalized by surface
         262             :        ! pressure @ i,j: sigma(i,j,l) = p(i,j,l) / ps(i,j)
         263           0 :        ALLOCATE(sigout(HcoState%NX,HcoState%NY,HcoState%NZ+1),STAT=AS)
         264           0 :        IF ( AS/= 0 ) THEN
         265           0 :           CALL HCO_ERROR( 'Cannot allocate sigout', RC )
         266           0 :           RETURN
         267             :        ENDIF
         268           0 :        DO l = 1, HcoState%NZ+1
         269           0 :           sigout(:,:,l) = HcoState%Grid%PEDGE%Val(:,:,l) &
         270           0 :                         / HcoState%Grid%PEDGE%Val(:,:,1)
         271             :        ENDDO
         272             :
         273             :        ! Now find first level on output grid where all sigma levels are
         274             :        ! lower than the lowest sigma value on the input grid. This is
         275             :        ! the highest level that needs to be considered for regridding.
         276             :        IF ( ReduceVert ) THEN
         277             :           sigMin = MINVAL(LevEdge)
         278             :           DO l = 1, HcoState%NZ+1
         279             :              IF ( MINVAL(sigout(:,:,l)) < sigMin ) EXIT
         280             :           ENDDO
         281             :
         282             :           ! The output grid is at grid center, so use l-1. Must be at least one.
         283             :           NZOUT = max(1,l-1)
         284             :
         285             :        ! Use full vertical grid if vertical levels shall not be restricted
         286             :        ! to range of input data (default).
         287             :        ELSE
         288           0 :           NZOUT = HcoState%NZ
         289             :        ENDIF
         290             :     ENDIF
         291             :
         292             :     ! verbose mode
         293           0 :     IF ( verb ) THEN
         294           0 :        MSG = 'Do MESSy regridding: ' // TRIM(Lct%Dct%cName)
         295           0 :        CALL HCO_MSG(HcoState%Config%Err,MSG)
         296           0 :        WRITE(MSG,*) ' - SameGrid     ? ', SameGrid
         297           0 :        CALL HCO_MSG(HcoState%Config%Err,MSG)
         298           0 :        WRITE(MSG,*) ' - Model levels ? ', IsModelLev
         299           0 :        CALL HCO_MSG(HcoState%Config%Err,MSG)
         300             :     ENDIF
         301             :
         302             :     !-----------------------------------------------------------------
         303             :     ! Make sure output array is defined & allocated
         304             :     !-----------------------------------------------------------------
         305           0 :     IF ( Lct%Dct%Dta%SpaceDim == 2 ) THEN
         306             :        CALL FileData_ArrCheck( HcoState%Config, &
         307             :                                Lct%Dct%Dta, HcoState%NX, HcoState%NY, &
         308           0 :                                NTIME,       RC                         )
         309           0 :        IF ( RC /= HCO_SUCCESS ) THEN
         310           0 :            CALL HCO_ERROR( 'ERROR 1', RC, THISLOC=LOC )
         311           0 :            RETURN
         312             :        ENDIF
         313           0 :     ELSEIF ( Lct%Dct%Dta%SpaceDim == 3 ) THEN
         314             :        CALL FileData_ArrCheck( HcoState%Config, &
         315             :                                Lct%Dct%Dta, HcoState%NX, HcoState%NY, &
         316           0 :                                NZOUT,       NTIME,       RC            )
         317           0 :        IF ( RC /= HCO_SUCCESS ) THEN
         318           0 :            CALL HCO_ERROR( 'ERROR 2', RC, THISLOC=LOC )
         319           0 :            RETURN
         320             :        ENDIF
         321             :     ENDIF
         322             :
         323             :     !-----------------------------------------------------------------
         324             :     ! Do straight-forward mapping if input grid = output grid
         325             :     !-----------------------------------------------------------------
         326           0 :     IF ( SameGrid ) THEN
         327             :        MSG = 'Input grid seems to match output grid. ' // &
         328           0 :              'No regridding is performed: ' // TRIM(Lct%Dct%cName)
         329           0 :        CALL HCO_WARNING( HcoState%Config%Err, MSG, RC )
         330             :
         331             :        ! For every time slice...
         332           0 :        DO I = 1, NTIME
         333           0 :        DO L = 1, NZOUT
         334             :
         335             :           ! 3D data
         336           0 :           IF ( Lct%Dct%Dta%SpaceDim == 3 ) THEN
         337           0 :              Lct%Dct%Dta%V3(I)%Val(:,:,L) = NcArr(:,:,L,I)
         338             :
         339             :           ! 2D data
         340             :           ELSE
         341           0 :              Lct%Dct%Dta%V2(I)%Val(:,:) = NcArr(:,:,L,I)
         342             :           ENDIF
         343             :        ENDDO
         344             :        ENDDO
         345             :
         346             :        ! All done!
         347           0 :        CALL HCO_LEAVE ( HcoState%Config%Err, RC )
         348           0 :        RETURN
         349             :     ENDIF
         350             :
         351             :     !=================================================================
         352             :     ! MESSy regridding follows below
         353             :     !=================================================================
         354             :
         355             :     !-----------------------------------------------------------------
         356             :     ! Source grid description.
         357             :     ! This creates a MESSy axis object for the source grid.
         358             :     !-----------------------------------------------------------------
         359           0 :     lon   => LonEdge
         360           0 :     lat   => LatEdge
         361           0 :     sigma => LevEdge
         362             :
         363           0 :     CALL AXIS_CREATE( HcoState, lon, lat, sigma, axis_src, RC )
         364           0 :     IF ( RC /= HCO_SUCCESS ) THEN
         365           0 :         CALL HCO_ERROR( 'ERROR 3', RC, THISLOC=LOC )
         366           0 :         RETURN
         367             :     ENDIF
         368             :
         369             :     ! Free pointer
         370           0 :     lon   => NULL()
         371           0 :     lat   => NULL()
         372           0 :     sigma => NULL()
         373             :
         374             :     !-----------------------------------------------------------------
         375             :     ! Destination grid description.
         376             :     ! This creates a MESSy axis object for the target (=HEMCO) grid.
         377             :     !-----------------------------------------------------------------
         378             :
         379             :     ! Get horizontal grid directly from HEMCO state
         380           0 :     lon   => HcoState%Grid%XEDGE%Val(:,1)
         381           0 :     lat   => HcoState%Grid%YEDGE%Val(1,:)
         382           0 :     IF( ASSOCIATED(LevEdge) ) sigma => sigout(:,:,1:NZOUT+1)
         383             :
         384           0 :     CALL AXIS_CREATE( HcoState, lon, lat, sigma, axis_dst, RC )
         385           0 :     IF ( RC /= HCO_SUCCESS ) THEN
         386           0 :         CALL HCO_ERROR( 'ERROR 4', RC, THISLOC=LOC )
         387           0 :         RETURN
         388             :     ENDIF
         389             :
         390             :     ! Free pointer
         391           0 :     lon   => NULL()
         392           0 :     lat   => NULL()
         393           0 :     sigma => NULL()
         394           0 :     IF ( ALLOCATED(sigout) ) DEALLOCATE(sigout)
         395             :
         396             :     !-----------------------------------------------------------------
         397             :     ! Set all other regridding parameter
         398             :     !-----------------------------------------------------------------
         399             :
         400             :     ! rg_type denotes the regridding type for each array (i.e. time
         401             :     ! slice). Set to 'intensive quantity' for all concentrations (incl.
         402             :     ! unitless) data. Set to 'index distribution' for data marked as
         403             :     ! index data in the configuration file. This will remap discrete
         404             :     ! values without interpolation, i.e. each grid box on the new
         405             :     ! grid holds the value with most overlap in the original grid.
         406           0 :     ALLOCATE(rg_type(NTIME))
         407           0 :     IF ( HCO_IsIndexData(Lct%Dct%Dta%OrigUnit) ) THEN
         408           0 :        rg_type(:) = RG_IDX
         409           0 :        IF ( verb ) THEN
         410           0 :           MSG = ' - Remap as index data.'
         411           0 :           CALL HCO_MSG(HcoState%Config%Err,MSG)
         412             :        ENDIF
         413             :     ELSE
         414           0 :        rg_type(:) = RG_INT
         415           0 :        IF ( verb ) THEN
         416           0 :           MSG = ' - Remap as concentration data.'
         417           0 :           CALL HCO_MSG(HcoState%Config%Err,MSG)
         418             :        ENDIF
         419             :     ENDIF
         420             :
         421             :     !-----------------------------------------------------------------
         422             :     ! Number of times the regridding need to be performed. If input
         423             :     ! data is on model levels, the data is only horizontally regridded,
         424             :     ! e.g. the regridding routine is called for every horizontal level
         425             :     ! separately. The vertical interpolation is done afterwards using
         426             :     ! routine ModelLev_Interpolate.
         427             :     !-----------------------------------------------------------------
         428           0 :     NCALLS = 1
         429           0 :     IF ( IsModelLev .AND. .NOT. ASSOCIATED(LevEdge) .AND. NZIN > 1 ) THEN
         430           0 :        NCALLS = NZIN
         431           0 :        ALLOCATE(ArrOut(HcoState%NX,HcoState%NY,NZIN,NTIME),STAT=AS)
         432             :        IF ( AS /= 0 ) THEN
         433           0 :           CALL HCO_ERROR( 'Cannot allocate ArrOut', RC )
         434           0 :           RETURN
         435             :        ENDIF
         436           0 :        ArrOut = 0.0_sp
         437             :     ENDIF
         438             :
         439             :     ! Do for all level batches ...
         440           0 :     DO I = 1, NCALLS
         441             :
         442             :        ! ArrIn is the input array to be used
         443           0 :        IF ( NCALLS /= 1 ) THEN
         444           0 :           ArrIn => NcArr(:,:,I:I,:)
         445             :        ELSE
         446           0 :           ArrIn => NcArr
         447             :        ENDIF
         448             :
         449             :        !-----------------------------------------------------------------
         450             :        ! Map input array onto MESSy array. Different time slices are
         451             :        ! stored as individual vector elements of narr_src.
         452             :        !-----------------------------------------------------------------
         453           0 :        CALL HCO2MESSY( HcoState, ArrIn, narr_src, axis_src, RC )
         454           0 :        IF ( RC /= HCO_SUCCESS ) THEN
         455           0 :            CALL HCO_ERROR( 'ERROR 5', RC, THISLOC=LOC )
         456           0 :            RETURN
         457             :        ENDIF
         458             :
         459             :        !-----------------------------------------------------------------
         460             :        ! Do the regridding
         461             :        !-----------------------------------------------------------------
         462             :        CALL NREGRID(s=narr_src,      sax=axis_src, dax=axis_dst, d=narr_dst, &
         463           0 :                     rg_type=rg_type, sovl=sovl,    dovl=dovl,    rcnt=rcnt    )
         464             :
         465             :        !-----------------------------------------------------------------
         466             :        ! Map the destination array narr_dst onto the data vector in the
         467             :        ! HEMCO list container or onto the temporary array ArrOut.
         468             :        !-----------------------------------------------------------------
         469           0 :        CALL MESSY2HCO( HcoState, narr_dst, Lct, I, ArrOut, RC )
         470           0 :        IF ( RC /= HCO_SUCCESS ) THEN
         471           0 :            CALL HCO_ERROR( 'ERROR 6', RC, THISLOC=LOC )
         472           0 :            RETURN
         473             :        ENDIF
         474             :
         475             :        ! Cleanup
         476           0 :        ArrIn => NULL()
         477             :
         478             :     ENDDO !NCALLS
         479             :
         480             :     !-----------------------------------------------------------------
         481             :     ! If these are model levels, do vertical interpolation now
         482             :     !-----------------------------------------------------------------
         483           0 :     IF ( ASSOCIATED(ArrOut) ) THEN
         484           0 :        CALL ModelLev_Interpolate( HcoState, ArrOut, Lct, RC )
         485           0 :        IF ( RC /= HCO_SUCCESS ) THEN
         486           0 :            CALL HCO_ERROR( 'ERROR 7', RC, THISLOC=LOC )
         487           0 :            RETURN
         488             :        ENDIF
         489           0 :        DEALLOCATE(ArrOut)
         490             :     ENDIF
         491             :
         492             :     !-----------------------------------------------------------------
         493             :     ! Cleanup
         494             :     !-----------------------------------------------------------------
         495           0 :     DEALLOCATE( sovl, dovl, rcnt, STAT=AS)
         496           0 :     IF(AS/=0) THEN
         497           0 :        CALL HCO_ERROR('DEALLOCATION ERROR 1', RC )
         498           0 :        RETURN
         499             :     ENDIF
         500           0 :     NULLIFY(sovl, dovl, rcnt)
         501             :
         502           0 :     DO I=1, SIZE(narr_dst)
         503           0 :        CALL INIT_NARRAY(narr_dst(I))
         504             :     ENDDO
         505           0 :     DEALLOCATE(narr_dst, STAT=AS)
         506             :     IF(AS/=0) THEN
         507           0 :        CALL HCO_ERROR('DEALLOCATION ERROR 3', RC )
         508           0 :        RETURN
         509             :     ENDIF
         510             :     NULLIFY(narr_dst)
         511             :
         512           0 :     DO I=1, SIZE(narr_src)
         513           0 :        CALL INIT_NARRAY(narr_src(I))
         514             :     ENDDO
         515           0 :     DEALLOCATE(narr_src, STAT=AS)
         516             :     IF(AS/=0) THEN
         517           0 :        CALL HCO_ERROR('DEALLOCATION ERROR 2', RC )
         518           0 :        RETURN
         519             :     ENDIF
         520             :     NULLIFY(narr_src)
         521             :
         522           0 :     DEALLOCATE( rg_type )
         523             :     rg_type => NULL()
         524             :
         525           0 :     CALL AXIS_DELETE( axis_src, axis_dst, RC )
         526           0 :     IF ( RC /= HCO_SUCCESS ) THEN
         527           0 :         CALL HCO_ERROR( 'ERROR 8', RC, THISLOC=LOC )
         528           0 :         RETURN
         529             :     ENDIF
         530             :
         531             :     ! Return w/ success
         532           0 :     CALL HCO_LEAVE ( HcoState%Config%Err, RC )
         533             :
         534           0 :   END SUBROUTINE HCO_MESSY_REGRID
         535             : !EOC
         536             : !------------------------------------------------------------------------------
         537             : !                   Harmonized Emissions Component (HEMCO)                    !
         538             : !------------------------------------------------------------------------------
         539             : !BOP
         540             : !
         541             : ! !IROUTINE: Axis_Create
         542             : !
         543             : ! !DESCRIPTION: Subroutine AXIS\_CREATE creates a MESSy axis type
         544             : ! from the grid defined by mid points Lon, Lat, Lev. Lev must be in
         545             : ! (unitless) sigma coordinates: sigma(i,j,l) = p(i,j,l) / p\_surface(i,j)
         546             : !\\
         547             : !\\
         548             : ! !INTERFACE:
         549             : !
         550           0 :   SUBROUTINE AXIS_CREATE( HcoState, lon, lat, lev, ax, RC )
         551             : !
         552             : ! !USES:
         553             : !
         554             :   USE MESSY_NCREGRID_BASE,  ONLY : INIT_AXIS
         555             : !
         556             : ! !INPUT PARAMETERS:
         557             : !
         558             :     TYPE(HCO_State),  POINTER                 :: HcoState
         559             :     TYPE(ListCont),   POINTER                 :: Lct
         560             :     REAL(hp),         POINTER                 :: Lon(:)
         561             :     REAL(hp),         POINTER                 :: Lat(:)
         562             :     REAL(hp),         POINTER                 :: Lev(:,:,:)
         563             : !
         564             : ! !INPUT/OUTPUT PARAMETERS:
         565             : !
         566             :     TYPE(axis),       POINTER                 :: ax(:)
         567             :     INTEGER,          INTENT(INOUT)           :: RC
         568             : !
         569             : ! !REVISION HISTORY:
         570             : !  22 Jun 2014 - C. Keller - Initial version (from messy_ncregrid_geohyb.f90)
         571             : !  See https://github.com/geoschem/hemco for complete history
         572             : !EOP
         573             : !------------------------------------------------------------------------------
         574             : !BOC
         575             : !
         576             : ! !ROUTINE ARGUMENTS:
         577             : !
         578             :     INTEGER                :: fID, I, J, N, status
         579             :     INTEGER                :: LOW, UPP
         580             :     INTEGER                :: XLON, YLAT, XLEV, YLEV, ZLEV
         581             :     INTEGER                :: ndp, nlev, cnt
         582             :     INTEGER                :: vtype
         583             :     INTEGER                :: ndep_lon, ndep_lat
         584             :     CHARACTER(LEN=255)     :: MSG, LOC
         585             :
         586             :     !=================================================================
         587             :     ! AXIS_CREATE begins here
         588             :     !=================================================================
         589             :
         590             :     ! For error handling
         591           0 :     LOC = 'AXIS_CREATE (HCOIO_MESSY_MOD.F90)'
         592           0 :     CALL HCO_ENTER ( HcoState%Config%Err, LOC, RC )
         593           0 :     IF ( RC /= HCO_SUCCESS ) THEN
         594           0 :         CALL HCO_ERROR( 'ERROR 9', RC, THISLOC=LOC )
         595           0 :         RETURN
         596             :     ENDIF
         597             :
         598             :     ! ----------------------------------------------------------------
         599             :     ! Pass horizontal grid dimensions to local variables. For now,
         600             :     ! assume all grids to be regular, i.e. no curvilinear grids.
         601             :     ! In this case, lon and lat dimensions are 1D-vectors.
         602             :     ! ----------------------------------------------------------------
         603             :
         604             :     ! ndep_lon and ndep_lat are the axis number of lon and lat. Needed
         605             :     ! if we have to set dependencies for the vertical axis.
         606           0 :     ndep_lon = 0
         607           0 :     ndep_lat = 0
         608             :
         609             :     ! Count axes and get grid dimensions
         610             :     ! Last axis is always 'free' dimension. This is required for the
         611             :     ! regridding to work properly.
         612           0 :     N = 1
         613           0 :     IF ( ASSOCIATED(lon) ) N = N + 1
         614           0 :     IF ( ASSOCIATED(lat) ) N = N + 1
         615           0 :     IF ( ASSOCIATED(lev) ) N = N + 1
         616             :
         617             :     ! ALLOCATE AXIS
         618           0 :     ALLOCATE(ax(N), STAT=status)
         619           0 :     IF ( status /= 0 ) THEN
         620           0 :        CALL HCO_ERROR ( 'Cannot allocate axis', RC )
         621           0 :        RETURN
         622             :     ENDIF
         623           0 :     DO I=1, N
         624           0 :        CALL INIT_AXIS(ax(I))
         625             :     ENDDO
         626             :
         627             :     ! N is the current axis
         628           0 :     N = 0
         629             :
         630             :     ! ----------------------------------------------------------------
         631             :     ! Assign longitude: this is always the first dimension
         632             :     ! ----------------------------------------------------------------
         633           0 :     IF ( ASSOCIATED(lon) ) THEN
         634           0 :        N        = N + 1
         635           0 :        ax(N)%lm = .true.     ! LONGITUDE IS MODULO AXIS
         636             :
         637             :        ! Axis dimension
         638           0 :        XLON = SIZE(lon,1)
         639             :
         640             :        ! FOR NOW, ASSUME NO DEPENDENCIES. NEED TO EDIT HERE
         641             :        ! IF WE WANT TO USE CURVILINEAR GRIDS
         642           0 :        ax(N)%ndp    = 1          ! LONGITUDE IS ...
         643           0 :        ALLOCATE(ax(N)%dep(1), STAT=status)
         644             :        IF ( status/= 0 ) THEN
         645           0 :           CALL HCO_ERROR ( 'Cannot allocate lon dependencies', RC )
         646           0 :           RETURN
         647             :        ENDIF
         648           0 :        ax(N)%dep(1) = N          ! ... INDEPENDENT
         649           0 :        ndep_lon     = N
         650             :
         651           0 :        ax(N)%dat%n = 1          ! 1 dimension
         652           0 :        ALLOCATE(ax(N)%dat%dim(ax(N)%dat%n), STAT=status)
         653             :        IF ( status/= 0 ) THEN
         654           0 :           CALL HCO_ERROR( 'Cannot allocate lon dimensions', RC )
         655           0 :           RETURN
         656             :        ENDIF
         657           0 :        ax(N)%dat%dim(:) = XLON
         658             :
         659           0 :        ALLOCATE(ax(N)%dat%vd(XLON),STAT=status)
         660             :        IF ( status/= 0 ) THEN
         661           0 :           CALL HCO_ERROR( 'Cannot allocate lon axis', RC )
         662           0 :           RETURN
         663             :        ENDIF
         664           0 :        ax(N)%dat%vd(:) = lon
         665             :
         666             :     ENDIF !lon
         667             :
         668             :     ! ----------------------------------------------------------------
         669             :     ! Assign latitude: this is always the second dimension
         670             :     ! ----------------------------------------------------------------
         671           0 :     IF ( ASSOCIATED(lat) ) THEN
         672           0 :        N        = N + 1
         673           0 :        ax(N)%lm = .false.    ! LATITUDE IS NON-MODULO AXIS
         674             :
         675             :        ! Axis dimension
         676           0 :        YLAT = SIZE(lat,1)
         677             :
         678             :        ! FOR NOW, ASSUME NO DEPENDENCIES. NEED TO EDIT HERE
         679             :        ! IF WE WANT TO USE CURVILINEAR GRIDS
         680           0 :        ax(N)%ndp    = 1          ! LATITUDE IS ...
         681           0 :        ALLOCATE(ax(N)%dep(1), STAT=status)
         682             :        IF ( status/= 0 ) THEN
         683           0 :           CALL HCO_ERROR( 'Cannot allocate lat dependencies', RC )
         684           0 :           RETURN
         685             :        ENDIF
         686           0 :        ax(N)%dep(1) = N          ! ... INDEPENDENT
         687           0 :        ndep_lat     = N
         688             :
         689           0 :        ax(N)%dat%n = 1          ! 1 dimension
         690           0 :        ALLOCATE(ax(N)%dat%dim(ax(N)%dat%n), STAT=status)
         691             :        IF ( status/= 0 ) THEN
         692           0 :           CALL HCO_ERROR( 'Cannot allocate lat dimensions', RC )
         693           0 :           RETURN
         694             :        ENDIF
         695           0 :        ax(N)%dat%dim(:) = YLAT
         696             :
         697           0 :        ALLOCATE(ax(N)%dat%vd(YLAT),STAT=status)
         698             :        IF ( status/= 0 ) THEN
         699           0 :           CALL HCO_ERROR( 'Cannot allocate lat axis', RC )
         700           0 :           RETURN
         701             :        ENDIF
         702           0 :        ax(N)%dat%vd(:) = lat
         703             :
         704             :        ! TAKE INTO ACCOUNT SPHERICAL GEOMETRY ...
         705           0 :        ax(N)%dat%vd = COS( ( (ax(N)%dat%vd - 90.0_dp) / 180.0_dp ) * &
         706           0 :                            HcoState%Phys%PI )
         707             :
         708             :     ENDIF !lat
         709             :
         710             :     ! ----------------------------------------------------------------
         711             :     ! Assign vertical levels (if defined): this is the 3rd dimension.
         712             :     ! The vertical axis is assumed to be in sigma coordinates.
         713             :     ! ----------------------------------------------------------------
         714           0 :     IF ( ASSOCIATED(lev) ) THEN
         715             :
         716             :        ! -------------------------------------------------------------
         717             :        ! Initialize vertical axis. Set dependencies on other axis.
         718             :        ! Define dimensions in the following order: lev, lat, lon.
         719             :        ! (first dimension has to be the 'own' dimension!).
         720             :        ! -------------------------------------------------------------
         721           0 :        N        = N + 1
         722           0 :        ax(N)%lm = .false.     ! VERTICAL AXIS IS NON-MODULO AXIS
         723             :
         724             :        ! Axis dimension
         725           0 :        XLEV = SIZE(lev,1)
         726           0 :        YLEV = SIZE(lev,2)
         727           0 :        ZLEV = SIZE(lev,3)
         728             :
         729             :        ! Sanity check: if XLEV and YLEV are > 1, they must correspond
         730             :        ! to the lon/lat axis defined above
         731           0 :        IF ( XLEV > 1 .AND. XLEV /= (XLON-1) ) THEN
         732           0 :           CALL HCO_ERROR ( 'level lon has wrong dimension', RC )
         733           0 :           RETURN
         734             :        ENDIF
         735           0 :        IF ( YLEV > 1 .AND. YLEV /= (YLAT-1) ) THEN
         736           0 :           CALL HCO_ERROR ( 'level lat has wrong dimension', RC )
         737           0 :           RETURN
         738             :        ENDIF
         739             :
         740             :        ! Set dependencies. First dimension must be vertical axis!
         741           0 :        ndp = 1
         742           0 :        IF ( YLEV > 1 ) ndp = ndp + 1
         743           0 :        IF ( XLEV > 1 ) ndp = ndp + 1
         744             :
         745           0 :        ax(N)%ndp = ndp
         746           0 :        ALLOCATE(ax(N)%dep(ax(N)%ndp), STAT=status)
         747             :        IF ( status /= 0 ) THEN
         748           0 :           CALL HCO_ERROR( 'Cannot allocate lev dependencies', RC )
         749           0 :           RETURN
         750             :        ENDIF
         751             :
         752             :        ! The variable dat%n holds the number of axis that level depends
         753             :        ! upon, and dat%dim are the corresponding axis dimensions (lengths).
         754           0 :        ax(N)%dat%n = ndp
         755           0 :        ALLOCATE(ax(N)%dat%dim(ax(N)%dat%n), STAT=status)
         756             :        IF ( status/= 0 ) THEN
         757           0 :           CALL HCO_ERROR( 'Cannot allocate lat dimensions', RC )
         758           0 :           RETURN
         759             :        ENDIF
         760             :
         761             :        ! Set axis indices and lengths:
         762             :        ! - First dimension is vertical axis
         763           0 :        cnt                = 1
         764           0 :        ax(N)%dep(cnt)     = N
         765           0 :        ax(N)%dat%dim(cnt) = ZLEV
         766           0 :        nlev               = ZLEV  ! number of grid cells
         767             :        ! - Next dimension is latitude
         768           0 :        IF ( YLEV > 1 ) THEN
         769           0 :           cnt                = cnt + 1
         770           0 :           ax(N)%dep(cnt)     = ndep_lat
         771           0 :           ax(N)%dat%dim(cnt) = YLEV
         772           0 :           nlev               = nlev * YLEV
         773             :        ENDIF
         774             :        ! - Next dimension is longitude
         775           0 :        IF ( XLEV > 1 ) THEN
         776           0 :           cnt                = cnt + 1
         777           0 :           ax(N)%dep(cnt)     = ndep_lon
         778           0 :           ax(N)%dat%dim(cnt) = XLEV
         779           0 :           nlev               = nlev * XLEV
         780             :        ENDIF
         781             :
         782           0 :        ALLOCATE(ax(N)%dat%vd(nlev),STAT=status)
         783             :        IF ( status/= 0 ) THEN
         784           0 :           CALL HCO_ERROR( 'Cannot allocate lat axis', RC )
         785           0 :           RETURN
         786             :        ENDIF
         787             :
         788             :        ! -------------------------------------------------------------
         789             :        ! Pass vertical sigma coordinates to vector
         790             :        ! -------------------------------------------------------------
         791             :        LOW = 1
         792           0 :        UPP = 1
         793           0 :        DO J = 1, YLEV !lat
         794           0 :        DO I = 1, XLEV !lon
         795           0 :           UPP                   = LOW + ZLEV - 1 ! Upper index
         796           0 :           ax(N)%dat%vd(LOW:UPP) = lev(I,J,:)     ! Pass to vector
         797           0 :           LOW                   = UPP + 1        ! Next lower index
         798             :        ENDDO
         799             :        ENDDO
         800             :
         801             :     END IF !lev
         802             :
         803           0 :     CALL HCO_LEAVE ( HcoState%config%Err, RC )
         804             :
         805             :   END SUBROUTINE AXIS_CREATE
         806             : !EOC
         807             : !------------------------------------------------------------------------------
         808             : !                   Harmonized Emissions Component (HEMCO)                    !
         809             : !------------------------------------------------------------------------------
         810             : !BOP
         811             : !
         812             : ! !IROUTINE: Axis_Delete
         813             : !
         814             : ! !DESCRIPTION: Subroutine AXIS\_DELETE deletes the specified MESSy
         815             : ! axis.
         816             : !\\
         817             : !\\
         818             : ! !INTERFACE:
         819             : !
         820           0 :   SUBROUTINE AXIS_DELETE ( ax1, ax2, RC )
         821             : !
         822             : ! !USES:
         823             : !
         824             :   USE MESSY_NCREGRID_BASE,  ONLY : INIT_AXIS
         825             : !
         826             : ! !INPUT/OUTPUT PARAMETERS:
         827             : !
         828             :     TYPE(axis),       POINTER            :: ax1(:)
         829             :     TYPE(axis),       POINTER            :: ax2(:)
         830             :     INTEGER,          INTENT(INOUT)      :: RC
         831             : !
         832             : ! !REVISION HISTORY:
         833             : !  28 Aug 2013 - C. Keller - Initial version
         834             : !  See https://github.com/geoschem/hemco for complete history
         835             : !EOP
         836             : !------------------------------------------------------------------------------
         837             : !BOC
         838             : !
         839             : ! !ROUTINE ARGUMENTS:
         840             : !
         841             :     INTEGER     :: I, status
         842             :
         843             :     !=================================================================
         844             :     ! AXIS_DELETE begins here
         845             :     !=================================================================
         846             :
         847           0 :     DO I=1, SIZE(ax1)
         848           0 :        CALL INIT_AXIS(ax1(I))
         849           0 :        CALL INIT_AXIS(ax2(I))
         850             :     ENDDO
         851           0 :     DEALLOCATE(ax1, ax2, STAT=status)
         852           0 :     IF(status/=0) THEN
         853           0 :        CALL HCO_ERROR('AXIS DEALLOCATION ERROR', RC )
         854           0 :        RETURN
         855             :     ENDIF
         856           0 :     NULLIFY(ax1)
         857           0 :     NULLIFY(ax2)
         858             :
         859             :   END SUBROUTINE AXIS_DELETE
         860             : !EOC
         861             : !------------------------------------------------------------------------------
         862             : !                   Harmonized Emissions Component (HEMCO)                    !
         863             : !------------------------------------------------------------------------------
         864             : !BOP
         865             : !
         866             : ! !IROUTINE: Hco2Messy
         867             : !
         868             : ! !DESCRIPTION: Subroutine HCO2MESSY converts a HEMCO data array into a
         869             : ! messy array-structure.
         870             : !\\
         871             : !\\
         872             : ! !INTERFACE:
         873             : !
         874           0 :   SUBROUTINE HCO2MESSY( HcoState, InArr, narr, ax, RC )
         875             : !
         876             : ! !USES:
         877             : !
         878             : !
         879             : ! !INPUT/OUTPUT PARAMETERS:
         880             : !
         881             :     TYPE(HCO_State),  POINTER        :: HcoState
         882             :     REAL(sp),         POINTER        :: InArr(:,:,:,:)
         883             :     TYPE(narray),     POINTER        :: narr(:)
         884             :     TYPE(axis),       POINTER        :: ax(:)
         885             :     INTEGER,          INTENT(INOUT)  :: RC
         886             : !
         887             : ! !REVISION HISTORY:
         888             : !  27 Jun 2014 - C. Keller - Initial version
         889             : !  See https://github.com/geoschem/hemco for complete history
         890             : !EOP
         891             : !------------------------------------------------------------------------------
         892             : !BOC
         893             : !
         894             : ! !ROUTINE ARGUMENTS:
         895             : !
         896             :     INTEGER            :: NCELLS, NX, NY, NZ, NT
         897             :     INTEGER            :: TMP
         898             :     INTEGER            :: status
         899             :     INTEGER            :: J, L, T, LOW, UPP
         900             :     CHARACTER(LEN=255) :: MSG, LOC
         901             :
         902             :     !=================================================================
         903             :     ! HCO2MESSY begins here
         904             :     !=================================================================
         905             :
         906             :     ! For error handling
         907           0 :     LOC = 'HCO2MESSY (HCOIO_MESSY_MOD.F90)'
         908             :
         909             :     ! ----------------------------------------------------------------
         910             :     ! Number of grid cells
         911             :     ! ----------------------------------------------------------------
         912           0 :     NX = SIZE(InArr,1)
         913           0 :     NY = SIZE(InArr,2)
         914           0 :     NZ = SIZE(InArr,3)
         915           0 :     NT = SIZE(InArr,4)
         916           0 :     NCELLS = NX * NY * NZ
         917             :
         918             :     ! create
         919           0 :     ALLOCATE(narr(NT),STAT=status)
         920           0 :     IF(status/=0) THEN
         921           0 :        CALL HCO_ERROR( 'narr allocation error', RC, THISLOC=LOC )
         922           0 :        RETURN
         923             :     ENDIF
         924             :
         925             :     ! Set MESSy vector dimensions. Copy from source axis object.
         926             :     ! Array dimensions are always one less than axis interface dimensions!
         927             :     ! For last ('free') dimension, set dimension length to 1.
         928           0 :     DO T = 1, NT
         929           0 :        narr(T)%n = size(ax)
         930           0 :        ALLOCATE(narr(T)%dim(narr(T)%n),STAT=status)
         931             :        IF(status/=0) THEN
         932           0 :           CALL HCO_ERROR( 'Cannot allocate array dims', RC, THISLOC=LOC )
         933           0 :           RETURN
         934             :        ENDIF
         935             :
         936           0 :        DO J = 1, narr(T)%n
         937           0 :           IF ( ax(J)%dat%n > 0 ) THEN
         938           0 :              tmp = ax(J)%dat%dim(1)-1
         939             :           ELSE
         940             :              tmp = 1
         941             :           ENDIF
         942           0 :           narr(T)%dim(J) = tmp
         943             :        ENDDO
         944             :
         945           0 :        ALLOCATE(narr(T)%vd(NCELLS),STAT=status)
         946           0 :        IF(status/=0) THEN
         947           0 :           CALL HCO_ERROR( 'Cannot allocate array', RC, THISLOC=LOC )
         948           0 :           RETURN
         949             :        ENDIF
         950             :     ENDDO !T
         951             :
         952             :     ! ----------------------------------------------------------------
         953             :     ! Pass HEMCO data array in slices to MESSy data vector.
         954             :     ! The MESSy data vector is an 'unfolded' HEMCO data array, with
         955             :     ! the longitude dimension changing the fastest, then followed by
         956             :     ! latitude, and level. Hence, we pass the HEMCO array in slices
         957             :     ! along the longitude axis to the MESSy vector, starting with the
         958             :     ! slice at lat/lev position 1/1, followed by 2/1, 3/1, etc.
         959             :     ! For now, the MESSy vector is always of type double. In future,
         960             :     ! we may want to set it to the same type as the HEMCO array and
         961             :     ! use pointers to the HEMCO data slices!
         962             :     ! ----------------------------------------------------------------
         963             :
         964             :     ! Vector index counter
         965             :     LOW = 1
         966           0 :     UPP = 1
         967             :
         968             :     ! Loop over all higher dimensions. Don't loop over lon because we pass all
         969             :     ! lon-values in slices.
         970           0 :     DO L = 1, NZ  ! NZ is 1 for 2D arrays
         971           0 :     DO J = 1, NY
         972           0 :        UPP = LOW + NX - 1   ! Upper index of slice to be filled
         973           0 :        DO T = 1, NT
         974           0 :           narr(T)%vd(LOW:UPP) = InArr(:,J,L,T)   ! Pass this slice to vector
         975             :        ENDDO
         976           0 :        LOW = UPP + 1        ! Next slice begins right after this one
         977             :     ENDDO !J
         978             :     ENDDO !L
         979             :
         980             :     ! Return w/ success
         981           0 :     RC = HCO_SUCCESS
         982             :
         983             :   END SUBROUTINE HCO2MESSY
         984             : !EOC
         985             : !------------------------------------------------------------------------------
         986             : !                   Harmonized Emissions Component (HEMCO)                    !
         987             : !------------------------------------------------------------------------------
         988             : !BOP
         989             : !
         990             : ! !IROUTINE: Messy2Hco
         991             : !
         992             : ! !DESCRIPTION: Subroutine MESSY2HCO converts a MESSy array structure
         993             : ! into a HEMCO data array. This is basically the reverse function of
         994             : ! HCO2MESSY.
         995             : !\\
         996             : !\\
         997             : ! !INTERFACE:
         998             : !
         999           0 :   SUBROUTINE MESSY2HCO( HcoState, narr, Lct, LEV, Ptr4D, RC )
        1000             : !
        1001             : ! !USES:
        1002             : !
        1003             : !
        1004             : ! !INPUT/OUTPUT PARAMETERS:
        1005             : !
        1006             :     TYPE(HCO_State),  POINTER                :: HcoState
        1007             :     TYPE(narray),     POINTER                :: narr(:)
        1008             :     TYPE(ListCont),   POINTER                :: Lct
        1009             :     INTEGER,          INTENT(IN   )          :: LEV
        1010             :     REAL(sp),         POINTER                :: Ptr4D(:,:,:,:)
        1011             :     INTEGER,          INTENT(INOUT)          :: RC
        1012             : !
        1013             : ! !REVISION HISTORY:
        1014             : !  27 Jun 2014 - C. Keller - Initial version
        1015             : !  See https://github.com/geoschem/hemco for complete history
        1016             : !EOP
        1017             : !------------------------------------------------------------------------------
        1018             : !BOC
        1019             : !
        1020             : ! !ROUTINE ARGUMENTS:
        1021             : !
        1022             :     INTEGER            :: J, L, T
        1023             :     INTEGER            :: LOW, UPP
        1024             :     INTEGER            :: NX, NY, NZ, NT, Z1, Z2
        1025             :     CHARACTER(LEN=255) :: MSG
        1026             :
        1027             :     !=================================================================
        1028             :     ! MESSY2HCO begins here
        1029             :     !=================================================================
        1030             :
        1031             :     ! Grid dimensions
        1032           0 :     NX = HcoState%NX
        1033           0 :     NY = HcoState%NY
        1034           0 :     NT = SIZE(narr)
        1035             :
        1036             :     ! Vertical levels to be filled on output array. Only level LEV
        1037             :     ! is filled if NC is different than one!
        1038           0 :     IF ( Lct%Dct%Dta%SpaceDim == 2 ) THEN
        1039             :        Z1 = 1
        1040             :        Z2 = 1
        1041           0 :     ELSEIF ( Lct%Dct%Dta%SpaceDim == 3 ) THEN
        1042           0 :        IF ( ASSOCIATED(Ptr4D) ) THEN
        1043           0 :           Z1 = LEV
        1044           0 :           Z2 = LEV
        1045             :        ELSE
        1046           0 :           Z1 = 1
        1047           0 :           Z2 = SIZE(Lct%Dct%Dta%V3(1)%Val,3)
        1048             :        ENDIF
        1049             :     ENDIF
        1050             :
        1051             :     ! Vector index counter
        1052           0 :     LOW = 1
        1053           0 :     UPP = 1
        1054             :
        1055             :     ! If temporary array Ptr4D is provided, make sure that its
        1056             :     ! dimensions are correct
        1057           0 :     IF ( ASSOCIATED(Ptr4D) ) THEN
        1058           0 :        NZ = Z2 - Z1 + 1
        1059             :        IF ( ( SIZE(Ptr4D,1) /= NX ) .OR. &
        1060             :             ( SIZE(Ptr4D,2) /= NY ) .OR. &
        1061           0 :             ( SIZE(Ptr4D,3) /= NZ ) .OR. &
        1062             :             ( SIZE(Ptr4D,4) /= NT )       ) THEN
        1063           0 :           WRITE(MSG,*) 'Temporary pointer has wrong dimensions: ', &
        1064           0 :                        TRIM(Lct%Dct%cName), NX, NY, NZ, NT, SIZE(Ptr4D,1), SIZE(Ptr4D,2), SIZE(Ptr4D,3), SIZE(Ptr4D,4)
        1065             :           CALL HCO_ERROR ( MSG, RC, &
        1066           0 :                            THISLOC='MESSY2HCO (hcoio_messy_mod.F90)' )
        1067           0 :           RETURN
        1068             :        ENDIF
        1069             :     ENDIF
        1070             :
        1071             :     ! Loop over all higher dimensions
        1072           0 :     DO L = Z1, Z2
        1073           0 :     DO J = 1, NY
        1074             :
        1075             :        ! Upper index of slice to be filled
        1076           0 :        UPP = LOW + NX - 1
        1077             :
        1078             :        ! Do for every time slice
        1079           0 :        DO T = 1, NT
        1080             :
        1081             :           ! If provided, write into temporary array
        1082           0 :           IF ( ASSOCIATED(Ptr4D) ) THEN
        1083           0 :              Ptr4D(:,J,L,T) = narr(T)%vd(LOW:UPP)
        1084             :
        1085             :           ! If temporary array Ptr4D is not provided, write directly
        1086             :           ! into list container array
        1087             :           ELSE
        1088             :
        1089             :              ! 2D
        1090           0 :              IF ( Lct%Dct%Dta%SpaceDim == 2 ) THEN
        1091           0 :                 Lct%Dct%Dta%V2(T)%Val(:,J) = narr(T)%vd(LOW:UPP)
        1092             :              ! 3D
        1093             :              ELSE
        1094           0 :                 Lct%Dct%Dta%V3(T)%Val(:,J,L) = narr(T)%vd(LOW:UPP)
        1095             :              ENDIF
        1096             :
        1097             :           ENDIF
        1098             :        ENDDO !NT
        1099             :
        1100             :        ! Next slice begins right after this one
        1101           0 :        LOW = UPP + 1
        1102             :     ENDDO !J
        1103             :     ENDDO !L
        1104             :
        1105             :     ! Return w/ success
        1106           0 :     RC = HCO_SUCCESS
        1107             :
        1108             :   END SUBROUTINE MESSY2HCO
        1109             : !EOC
        1110             : END MODULE HCOIO_MESSY_MOD

Generated by: LCOV version 1.14