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