ecearth3/sources/tm5mp/base/grid_type_hyb.F90

!------------------------
!  MODULE AND DATA TYPES
!------------------------
!
!    USE  grid_type_hyb
!    type(TLevelInfo) :: levi
!
!------------------------
!  INITIALIZATION
!------------------------
!    Three possible ways:
!    
!   o by half level coeff  a (Pa)  and   b (0-1)
!  
!      call Init( levi, 60, (/a0,..,a60/), (/b0,..,b60/), status )
!  
!   o by a character key; yet supported:
!       'ec19'   : echam 19 levels
!       'ec31'   : ecmwf 31 levels
!       'ec4'    : ecmwf  4 levels ( 4 levels extracted from 91 levels in IFS : TM5 sees a  4-levels IFS)
!       'ec10'   : ecmwf 10 levels (10 levels extracted from 91 levels in IFS : TM5 sees a 10-levels IFS)
!       'ec34'   : ecmwf 34 levels (34 levels extracted from 91 levels in IFS : TM5 sees a 34-levels IFS)
!       'ec40'   : ecmwf 40 levels
!       'ec60'   : ecmwf 60 levels
!       'ec62'   : ecmwf 62 levels
!       'ec91'   : ecmwf 91 levels
!       'tm31'   : tm 31 levels (reversed ec31)
!       'tm40'   : tm 40 levels (reversed ec40)
!       'tm60'   : tm 60 levels (reversed ec60)
!       'tm62'   : tm 62 levels (reversed ec62)
!       'tm91'   : tm 91 levels (reversed ec91)
!       
!      call Init( levi, 'ec60', status )
!  
!   o define a selection of half levels:
!
!      call Init( levi_parent, 'ec60', status )
!      call Init( levi, levi_parent, (/0,2,..,58,60/), status )
!
!--------------------------
!  VERTICAL REGRIDDING (I)
!--------------------------
!
! To regrid an array (InArr) defined on one grid (InGrid) into an array
! (OutArr) onto another vertical grid (OutGrid), use FillLevels.
!
! Input array can be 1D (column), 2D (latitudinal or longitudinal slice), or
! 3D. Data can be on full or half levels, except for 1D (only full levels, no
! interpolation available: easy to code though. FIXME).
! 
! The regridding is done according to 3 more inputs:
! 
!   (1) a combination key, used if 'combination of' or 'interpolation between'
!       input is needed:
!  
!       'bottom'     :  use closest-to-the-ground neighbor
!       'top'        :  use closest-to-the-model-top neighbor
!       'sum'        :  sum input values (use for mass [extensive])
!       'aver'       :  unweighted average across levels
!       'mass-aver'  :  air-mass weighted average (use for mixing ratio [intensive])
!
!   (2) a location key (nw), that determines if the data are pinned to the
!       edge or the center of the boxes:
! 
!       'n' : full levels, a/k/a center of box
!       'w' : half-levels, a/k/a box edges
!
!   (3) surface pressure (SP), which is one dimension less than InArr
!
!
! Call is:
! 
!   FillLevels( OutGrid,     SP, OutArr, InGrid, InArr, combine_key, status )  ! for 1D
!   FillLevels( OutGrid, nw, SP, OutArr, InGrid, InArr, combine_key, status )  ! for 2D and 3D
!  
! Note only two output : OutArr (regridded data) and Status (0/1=success/failure)
!   
!--------------------------
!  VERTICAL REGRIDDING (II)
!--------------------------
! There is a remapping routine specifically for cases when one of the levels set
! is a subset of the other, be it target or source grid.
!
!  call FillLevelsParents(levi, nw, ArrayIn, combine_key, ArrayOut, status, sp)
!
! - Works both ways:  parent-to-subset  and  subset-to-parent. Direction is determined by the arrays size.
! - Surface pressure is optional, since not always needed. If needed and not passed, run stops w/ a message.
! - 'nw' and 'combine_key' are the same as for FillLevels, but not all cases have been implemented yet.
! - Arrays follow the order ('u' or 'd') of the levels they are attached to.
! 
!------------------------
!  CHECKING
!------------------------
! You can check that a LevI is initalized:
!
!    call check(levi, status)
!
! or that a 3D array conforms to the LevI (3rd dim is same as nb of levels):
!
!    call check(levi, Arr, status)
!
! or just examine the content of a LevI:
!
!    call PrintInfo( LevI )
!
!------------------------
!  PRESSURE
!------------------------
! You can get the pressure at each half/full level of LevI from surface
! pressure sp:
!
!   o call H/FPressure( levi, sp, pressure, status)
!
! Output pressure can be 1-,2- or 3D. Then SP has one dimension less than pressure.
!
!------------------------
!  TODO - ISSUE
!------------------------
! When remapping:
!
! It is assumed that 'HEIGHT-AVE' is akin to 'AIR-MASS-AVER'. This is not
! true, and is a crude approximation at best (going from coarse to fine grid,
! results should be ok, but not when going from fine-to-coarse grid). Problem
! is to get grid box height on the non-TM5 levels. See comment in
! levi_FillParent_3d.
!
! In TM5, it is used for UDDR and DDDR. 
!-

#define TRACEBACK write (gol,'("in ",a," (",a,", line",i5,")")') rname, __FILE__, __LINE__; call goErr
#define IF_NOTOK_RETURN(action) if (status/=0) then; TRACEBACK; action; return; end if
#define IF_ERROR_RETURN(action) if (status> 0) then; TRACEBACK; action; return; end if

MODULE GRID_TYPE_HYB

  use go, only : gol, goPr, goErr
  
  implicit none
  private
  
  public :: TLevelInfo
  public :: Init, Done ! life cycle of a TLevelInfo instance
  public :: Check
  public :: Compare
  public :: HPressure  ! returns pressure [Pa] at each half level of LevI
  public :: FPressure  ! returns pressure [Pa] at each full level of LevI
  public :: FillLevels ! remapping utility b/w two LevelInfo
  public :: FillLevelsParents ! remapping utility b/w LevelInfo and its parents
  public :: PrintInfo
 
  ! --- const ---------------------------------------
  
  character(len=*), parameter  ::  mname = 'grid_type_hyb'
  
  ! --- types ---------------------------
  
  type TLevelInfo
     character(len=32) ::  name                           ! name used for messages
     integer           ::  nlev                           ! number of levels
     real,    pointer  ::  a(:), b(:)                     ! hybride half level coeff; indices 0, 1, ..., nlev
     real,    pointer  ::  fa(:), fb(:)                   ! hybride full level coeff; indices 1, ..., nlev
     real,    pointer  ::  fa_bounds(:,:), fb_bounds(:,:) ! half level coeff per full level (1:2,1:nlev)
     real,    pointer  ::  p0(:), fp0(:), m0(:)           ! standard pressures and mass (surface pressure 1e5)
     character(len=1)  ::  updo                           ! upwards or downwards ?
     real,    pointer  ::  da(:), db(:)                   ! hybride increments; indices 1, ..., nlev
     ! -- parent levels if any
     logical           ::  selection   ! has parent?
     integer, pointer  ::  hlevs(:)    ! indices of parent half levels used in defining child
     integer, pointer  ::  flevs(:,:)  ! indice range of parent full levels covered by child level (nlev,2) 
     integer           ::  nlev_parent ! number of levels
     character(len=1)  ::  updo_parent ! upwards or downwards?
     real,    pointer  ::  a_parent(:), b_parent(:)   ! hybride half level coeff
     real,    pointer  ::  da_parent(:), db_parent(:) ! hybride increments of parents; indices 1, ..., nlev
  end type TLevelInfo
  
  ! --- const ------------------------------------------------------------------
  
  ! *** ec4 **************************************************************
  
  ! Hybrid coordinate coefficients at half level interfaces,
  ! specifying 4 vertical ECMWF levels. These are subset of 91 levels, used in EC-Earth.
  ! Coefficient a is given in Pa, b in [0,1] .
  ! Should match what is in levels/ml4/const_ec_v__ml4.F90
  
  real, parameter  :: a_ec4(0:4) = (/   &
       0.000000, &
       15638.053711, &
       8356.252930, &
       6.575628, &
       0.000000 /)

  real, parameter  :: b_ec4(0:4) = (/   &
       0.000000, &
       0.012508, &
       0.698224, &
       0.994204, &
       1.000000 /)

  ! *** ec10 **************************************************************
  
  ! Hybrid coordinate coefficients at half level interfaces,
  ! specifying 10 vertical ECMWF levels. These are subset of 91 levels, used in EC-Earth.
  ! Coefficient a is given in Pa, b in [0,1] .
  ! Should match what is in levels/ml10/const_ec_v__ml10.F90
  
  real, parameter  :: a_ec10(0:10) = (/   &
       0.000000, &
       450.685791, &
       3358.425781, &
       7341.469727, &
       15638.053711, &
       20319.011719, &
       14665.645508, &
       3010.146973, &
       336.772369, &
       6.575628, &
       0.000000 /)

  real, parameter  :: b_ec10(0:10) = (/   &
       0.000000, &
       0.000000, &
       0.000000, &
       0.000000, &
       0.012508, &
       0.176091, &
       0.475016, &
       0.875518, &
       0.973466, &
       0.994204, &
       1.000000 /)


  ! *** ec19 **************************************************************
  
  ! Hybrid coordinate coefficients at half level interfaces,
  ! specifying 19 vertical ECHAM levels.
  ! Coefficient a is given in Pa, b in [0,1] .

  real, parameter :: a_ec19(0:19) = (/ &
             0.000   ,  2000.000   ,  4000.000   ,  6046.111   , &
          8267.928   , 10609.510   , 12851.100   , 14698.500   , &
         15861.130   , 16116.240   , 15356.920   , 13621.460   , &
         11101.560   ,  8127.144   ,  5125.142   ,  2549.969   , &
           783.195   ,     0.000   ,     0.000   ,     0.000        /)

  real, parameter :: b_ec19(0:19) = (/ &
             0.000000000, 0.00000000, 0.00000000, 0.0003389933, &
             0.003357187, 0.01307004, 0.03407715, 0.0706498300, &
             0.125916700, 0.20119540, 0.29551960, 0.4054092000, &
             0.524932200, 0.64610790, 0.75969840, 0.8564376000, &
             0.928746900, 0.97298520, 0.99228150, 1.0000000000    /)

  ! *** ec31 **************************************************************
  
  ! Hybrid coordinate coefficients at half level interfaces,
  ! specifying 31  vertical ECMWF levels.
  ! Coefficient a is given in Pa, b in [0,1] .

  real, parameter :: a_ec31(0:31) = (/ &
        0.0         ,  2000.0         ,  4000.0         ,  6000.0         , &     
     8000.0         ,  9976.135361    , 11820.539617    , 13431.393926    , &     
    14736.356909    , 15689.207458    , 16266.610500    , 16465.005734    , &     
    16297.619332    , 15791.598604    , 14985.269630    , 13925.517858    , &     
    12665.291662    , 11261.228878    ,  9771.406290    ,  8253.212096    , &     
     6761.341326    ,  5345.914240    ,  4050.717678    ,  2911.569385    , &     
     1954.805296    ,  1195.889791    ,   638.148911    ,   271.626545    , &     
       72.063577    ,     0.0         ,     0.0         ,     0.0             /)        

  real, parameter :: b_ec31(0:31) = (/ &
        0.0         , 0.0             , 0.0             , 0.0             , &
        0.0         , 0.0003908582    , 0.0029197006    , 0.0091941320    , &
        0.0203191555, 0.0369748598    , 0.0594876397    , 0.0878949492    , &
        0.1220035886, 0.1614415235    , 0.2057032385    , 0.2541886223    , &
        0.3062353873, 0.3611450218    , 0.4182022749    , 0.4766881754    , &
        0.5358865832, 0.5950842740    , 0.6535645569    , 0.7105944258    , &
        0.7654052430, 0.8171669567    , 0.8649558510    , 0.9077158297    , &
        0.9442132326, 0.9729851852    , 0.9922814815    , 1.0                 /)   

  ! *** ec34 **************************************************************
  
  ! Hybrid coordinate coefficients at half level interfaces,
  ! specifying 34 vertical ECMWF levels. These are subset of 91 levels, used in EC-Earth.
  ! Coefficient a is given in Pa, b in [0,1] .
  
  real, parameter  :: a_ec34(0:34) = (/   &
       0.000000,                                                              &  !  0
      21.413612,                                    76.167656,                &  !  5
                     204.637451,                                  450.685791, &  ! 10
                                  857.945801,                                 &  ! 15
    1463.163940,                                  2292.155518,                &  ! 20
                    3358.425781,                                 4663.776367, &  ! 25
                                  6199.839355,                   7341.469727, &  ! 30
                    8564.624023,                  9873.560547,                &  ! 35
    11262.484375,                 12713.897461,                 14192.009766, &  ! 40
                   15638.053711,                 16990.623047,                &  ! 45
    18191.029297,                 19184.544922,                 19919.796875, &  ! 50
                   20348.916016,                                20319.011719, &  ! 55
                                  19348.775391,                               &  ! 60
    17385.595703,                                14665.645508,                &  ! 65
                   11543.166992,                                 8356.252930, &  ! 70
                                  5422.802734,                                &  ! 75
     3010.146973,                                 1297.656128,                &  ! 80
                     336.772369,                                    6.575628, &  ! 85
                       0.000000 /)                                               ! 90

  real, parameter  :: b_ec34(0:34) = (/   &
       0.000000,                                                             & !  0 
       0.000000,                                    0.000000,                & !  5 
                      0.000000,                                    0.000000, & ! 10 
                                     0.000000,                               & ! 15 
       0.000000,                                    0.000000,                & ! 20 
                      0.000000,                                    0.000000, & ! 25 
                                     0.000000,                     0.000000, & ! 30 
                      0.000055,                     0.000279,                & ! 35 
       0.001000,                     0.002765,                     0.006322, & ! 40 
                      0.012508,                     0.022189,                & ! 45 
       0.036227,                     0.055474,                     0.080777, & ! 50 
                      0.112979,                                    0.176091, & ! 55 
                                     0.259554,                               & ! 60 
       0.362203,                                    0.475016,                & ! 65 
                      0.589317,                                    0.698224, & ! 70 
                                     0.795385,                               & ! 75 
       0.875518,                                    0.935157,                & ! 80 
                      0.973466,                                    0.994204, & ! 85 
                      1.000000 /)                                              ! 90 

  ! *** ec40 **************************************************************

  ! Hybrid coordinate coefficients at half level interfaces,
  ! specifying 40 vertical ECMWF levels.
  ! Coefficient a is given in [Pa] .

  real, parameter  :: a_ec40(0:40) = (/   &
             0.000000,  2000.000000,  4000.000000,  6000.000000, &
          8000.000000,  9988.882813, 11914.523438, 13722.941406, &
         15369.730469, 16819.476563, 18045.183594, 19027.695313, &
         19755.109375, 20222.205078, 20429.863281, 20384.480469, &
         20097.402344, 19584.330078, 18864.750000, 17961.357422, &
         16899.468750, 15706.447266, 14411.124023, 13043.218750, &
         11632.758789, 10209.500977,  8802.356445,  7438.803223, &
          6144.314941,  4941.778320,  3850.913330,  2887.696533, &
          2063.779785,  1385.912598,   855.361755,   467.333588, &
           210.393890,    65.889244,     7.367743,     0.000000, &
             0.000000/)

  real, parameter  :: b_ec40(0:40) = (/   &
             0.0000000000, 0.0000000000, 0.0000000000, 0.0000000000,  &
             0.0000000000, 0.0001971156, 0.0015112918, 0.0048841573,  &
             0.0110761747, 0.0206778906, 0.0341211632, 0.0516904071,  &
             0.0735338330, 0.0996747017, 0.1300225258, 0.1643843055,  &
             0.2024759650, 0.2439331412, 0.2883229256, 0.3351548910,  &
             0.3838921785, 0.4339629412, 0.4847716093, 0.5357099175,  &
             0.5861684084, 0.6355474591, 0.6832686067, 0.7287858129,  &
             0.7715966105, 0.8112534285, 0.8473749161, 0.8796569109,  &
             0.9078838825, 0.9319403172, 0.9518215060, 0.9676452279,  &
             0.9796627164, 0.9882701039, 0.9940194488, 0.9976301193,  &
             1.0000000000 /)
             
  ! *** ec60 **************************************************************
  
  ! Hybrid coordinate coefficients at half level interfaces,
  ! specifying 60  vertical ECMWF levels.
  ! Coefficient a is given in Pa, b in [0,1] .

  real, parameter  :: a_ec60(0:60) = (/   &
             0.000000,    20.000000,    38.425343,    63.647804, &
            95.636963,   134.483307,   180.584351,   234.779053, &
           298.495789,   373.971924,   464.618134,   575.651001, &
           713.218079,   883.660522,  1094.834717,  1356.474609, &
          1680.640259,  2082.273926,  2579.888672,  3196.421631, &
          3960.291504,  4906.708496,  6018.019531,  7306.631348, &
          8765.053711, 10376.126953, 12077.446289, 13775.325195, &
         15379.805664, 16819.474609, 18045.183594, 19027.695313, &
         19755.109375, 20222.205078, 20429.863281, 20384.480469, &
         20097.402344, 19584.330078, 18864.750000, 17961.357422, &
         16899.468750, 15706.447266, 14411.124023, 13043.218750, &
         11632.758789, 10209.500977,  8802.356445,  7438.803223, &
          6144.314941,  4941.778320,  3850.913330,  2887.696533, &
          2063.779785,  1385.912598,   855.361755,   467.333588, &
           210.393890,    65.889244,     7.367743,     0.000000, &
             0.000000/)

  real, parameter  :: b_ec60(0:60) = (/   &
             0.00000000, 0.00000000, 0.00000000, 0.00000000,  &
             0.00000000, 0.00000000, 0.00000000, 0.00000000,  &
             0.00000000, 0.00000000, 0.00000000, 0.00000000,  &
             0.00000000, 0.00000000, 0.00000000, 0.00000000,  &
             0.00000000, 0.00000000, 0.00000000, 0.00000000,  &
             0.00000000, 0.00000000, 0.00000000, 0.00000000,  &
             0.00007582, 0.00046139, 0.00181516, 0.00508112,  &
             0.01114291, 0.02067788, 0.03412116, 0.05169041,  &
             0.07353383, 0.09967469, 0.13002251, 0.16438432,  &
             0.20247590, 0.24393314, 0.28832296, 0.33515489,  &
             0.38389215, 0.43396294, 0.48477158, 0.53570992,  &
             0.58616841, 0.63554746, 0.68326861, 0.72878581,  &
             0.77159661, 0.81125343, 0.84737492, 0.87965691,  &
             0.90788388, 0.93194032, 0.95182151, 0.96764523,  &
             0.97966272, 0.98827010, 0.99401945, 0.99763012,  &
             1.00000000 /)

  ! *** ec62 **************************************************************

  real, parameter  :: a_ec62(0:62) = (/   &
                                                                       0.000000, &
         988.835876,   1977.676270,   2966.516602,   3955.356934,   4944.197266, &
        5933.037598,   6921.870117,   7909.441406,   8890.707031,   9860.528320, &
       10807.783203,  11722.749023,  12595.006836,  13419.463867,  14192.009766, &
       14922.685547,  15638.053711,  16329.560547,  16990.623047,  17613.281250, &
       18191.029297,  18716.968750,  19184.544922,  19587.513672,  19919.796875, &
       20175.394531,  20348.916016,  20434.158203,  20426.218750,  20319.011719, &
       20107.031250,  19785.357422,  19348.775391,  18798.822266,  18141.296875, &
       17385.595703,  16544.585938,  15633.566406,  14665.645508,  13653.219727, &
       12608.383789,  11543.166992,  10471.310547,   9405.222656,   8356.252930, &
        7335.164551,   6353.920898,   5422.802734,   4550.215820,   3743.464355, &
        3010.146973,   2356.202637,   1784.854614,   1297.656128,    895.193542, &
         576.314148,    336.772369,    162.043427,     54.208336,      6.575628, &
           0.003160,      0.000000 /)

  real, parameter  :: b_ec62(0:62) = (/   &
                                                                       0.000000, &
           0.000000,      0.000000,      0.000000,      0.000000,      0.000000, &
           0.000000,      0.000000,      0.000013,      0.000087,      0.000275, &
           0.000685,      0.001415,      0.002565,      0.004187,      0.006322, &
           0.009035,      0.012508,      0.016860,      0.022189,      0.028610, &
           0.036227,      0.045146,      0.055474,      0.067316,      0.080777, &
           0.095964,      0.112979,      0.131935,      0.152934,      0.176091, &
           0.201520,      0.229315,      0.259554,      0.291993,      0.326329, &
           0.362203,      0.399205,      0.436906,      0.475016,      0.513280, &
           0.551458,      0.589317,      0.626559,      0.662934,      0.698224, &
           0.732224,      0.764679,      0.795385,      0.824185,      0.850950, &
           0.875518,      0.897767,      0.917651,      0.935157,      0.950274, &
           0.963007,      0.973466,      0.982238,      0.989153,      0.994204, &
           0.997630,      1.000000 /)
             
  ! *** ec91 **************************************************************
  
  ! Hybrid coordinate coefficients at half level interfaces,
  ! specifying 91 vertical ECMWF levels.
  ! Coefficient a is given in Pa, b in [0,1] .

  real, parameter  :: a_ec91(0:91) = (/   &
           0.000000,      2.000040,      3.980832,      7.387186,     12.908319, &
          21.413612,     33.952858,     51.746601,     76.167656,    108.715561, &
         150.986023,    204.637451,    271.356506,    352.824493,    450.685791, &
         566.519226,    701.813354,    857.945801,   1036.166504,   1237.585449, &
        1463.163940,   1713.709595,   1989.874390,   2292.155518,   2620.898438, &
        2976.302246,   3358.425781,   3767.196045,   4202.416504,   4663.776367, &
        5150.859863,   5663.156250,   6199.839355,   6759.727051,   7341.469727, &
        7942.926270,   8564.624023,   9208.305664,   9873.560547,  10558.881836, &
       11262.484375,  11982.662109,  12713.897461,  13453.225586,  14192.009766, &
       14922.685547,  15638.053711,  16329.560547,  16990.623047,  17613.281250, &
       18191.029297,  18716.968750,  19184.544922,  19587.513672,  19919.796875, &
       20175.394531,  20348.916016,  20434.158203,  20426.218750,  20319.011719, &
       20107.031250,  19785.357422,  19348.775391,  18798.822266,  18141.296875, &
       17385.595703,  16544.585938,  15633.566406,  14665.645508,  13653.219727, &
       12608.383789,  11543.166992,  10471.310547,   9405.222656,   8356.252930, &
        7335.164551,   6353.920898,   5422.802734,   4550.215820,   3743.464355, &
        3010.146973,   2356.202637,   1784.854614,   1297.656128,    895.193542, &
         576.314148,    336.772369,    162.043427,     54.208336,      6.575628, &
           0.003160,      0.000000 /)

  real, parameter  :: b_ec91(0:91) = (/   &
           0.000000,      0.000000,      0.000000,      0.000000,      0.000000, &
           0.000000,      0.000000,      0.000000,      0.000000,      0.000000, &
           0.000000,      0.000000,      0.000000,      0.000000,      0.000000, &
           0.000000,      0.000000,      0.000000,      0.000000,      0.000000, &
           0.000000,      0.000000,      0.000000,      0.000000,      0.000000, &
           0.000000,      0.000000,      0.000000,      0.000000,      0.000000, &
           0.000000,      0.000000,      0.000000,      0.000000,      0.000000, &
           0.000014,      0.000055,      0.000131,      0.000279,      0.000548, &
           0.001000,      0.001701,      0.002765,      0.004267,      0.006322, &
           0.009035,      0.012508,      0.016860,      0.022189,      0.028610, &
           0.036227,      0.045146,      0.055474,      0.067316,      0.080777, &
           0.095964,      0.112979,      0.131935,      0.152934,      0.176091, &
           0.201520,      0.229315,      0.259554,      0.291993,      0.326329, &
           0.362203,      0.399205,      0.436906,      0.475016,      0.513280, &
           0.551458,      0.589317,      0.626559,      0.662934,      0.698224, &
           0.732224,      0.764679,      0.795385,      0.824185,      0.850950, &
           0.875518,      0.897767,      0.917651,      0.935157,      0.950274, &
           0.963007,      0.973466,      0.982238,      0.989153,      0.994204, &
           0.997630,      1.000000 /)

  ! *** ec137 **************************************************************
  
  ! Hybrid coordinate coefficients at half level interfaces,
  ! specifying 137 vertical ECMWF levels.
  ! Coefficient a is given in Pa, b in [0,1] .
  
  real, parameter  :: a_ec137(0:137) = (/   &
       0.000000,      2.000365,      3.102241,      4.666084,      6.827977,      &
       9.746966,      13.605424,     18.608931,     24.985718,     32.985710,     &
       42.879242,     54.955463,     69.520576,     86.895882,     107.415741,    &
       131.425507,    159.279404,    191.338562,    227.968948,    269.539581,    &
       316.420746,    368.982361,    427.592499,    492.616028,    564.413452,    &
       643.339905,    729.744141,    823.967834,    926.344910,    1037.201172,   &
       1156.853638,   1285.610352,   1423.770142,   1571.622925,   1729.448975,   &
       1897.519287,   2076.095947,   2265.431641,   2465.770508,   2677.348145,   &
       2900.391357,   3135.119385,   3381.743652,   3640.468262,   3911.490479,   &
       4194.930664,   4490.817383,   4799.149414,   5119.895020,   5452.990723,   &
       5798.344727,   6156.074219,   6526.946777,   6911.870605,   7311.869141,   &
       7727.412109,   8159.354004,   8608.525391,   9076.400391,   9562.682617,   &
       10065.978516,  10584.631836,  11116.662109,  11660.067383,  12211.547852,  &
       12766.873047,  13324.668945,  13881.331055,  14432.139648,  14975.615234,  &
       15508.256836,  16026.115234,  16527.322266,  17008.789063,  17467.613281,  &
       17901.621094,  18308.433594,  18685.718750,  19031.289063,  19343.511719,  &
       19620.042969,  19859.390625,  20059.931641,  20219.664063,  20337.863281,  &
       20412.308594,  20442.078125,  20425.718750,  20361.816406,  20249.511719,  &
       20087.085938,  19874.025391,  19608.572266,  19290.226563,  18917.460938,  &
       18489.707031,  18006.925781,  17471.839844,  16888.687500,  16262.046875,  &
       15596.695313,  14898.453125,  14173.324219,  13427.769531,  12668.257813,  &
       11901.339844,  11133.304688,  10370.175781,  9617.515625,   8880.453125,   &
       8163.375000,   7470.343750,   6804.421875,   6168.531250,   5564.382813,   &
       4993.796875,   4457.375000,   3955.960938,   3489.234375,   3057.265625,   &
       2659.140625,   2294.242188,   1961.500000,   1659.476563,   1387.546875,   &
       1143.250000,   926.507813,    734.992188,    568.062500,    424.414063,    &
       302.476563,    202.484375,    122.101563,    62.781250,     22.835938,     &
       3.757813,      0.000000,      0.000000 /)
       
  real, parameter  :: b_ec137(0:137) = (/   &
       0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
       0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
       0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
       0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
       0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
       0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
       0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
       0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
       0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
       0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
       0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
       0.000007, 0.000024, 0.000059, 0.000112, 0.000199, &
       0.000340, 0.000562, 0.000890, 0.001353, 0.001992, &
       0.002857, 0.003971, 0.005378, 0.007133, 0.009261, &
       0.011806, 0.014816, 0.018318, 0.022355, 0.026964, &
       0.032176, 0.038026, 0.044548, 0.051773, 0.059728, &
       0.068448, 0.077958, 0.088286, 0.099462, 0.111505, &
       0.124448, 0.138313, 0.153125, 0.168910, 0.185689, &
       0.203491, 0.222333, 0.242244, 0.263242, 0.285354, &
       0.308598, 0.332939, 0.358254, 0.384363, 0.411125, &
       0.438391, 0.466003, 0.493800, 0.521619, 0.549301, &
       0.576692, 0.603648, 0.630036, 0.655736, 0.680643, &
       0.704669, 0.727739, 0.749797, 0.770798, 0.790717, &
       0.809536, 0.827256, 0.843881, 0.859432, 0.873929, &
       0.887408, 0.899900, 0.911448, 0.922096, 0.931881, &
       0.940860, 0.949064, 0.956550, 0.963352, 0.969513, &
       0.975078, 0.980072, 0.984542, 0.988500, 0.991984, &
       0.995003, 0.997630, 1.000000 /)

  
       
    ! *** nc28 **************************************************************
  
  ! Hybrid coordinate coefficients at half level interfaces,
  ! specifying the 28 sigma levels of the NCEP reanalysis (surf -> top).
  ! Coefficient a is given in Pa, b in [0,1] .

  real, parameter  :: a_nc28(0:28) = 0.0

  real, parameter  :: b_nc28(0:28) = (/   &
              1.0000  ,  0.9900  ,   0.9742  ,   0.9546 ,  &
              0.9304  ,  0.9014  ,   0.8662  ,   0.8254 ,  &
              0.7774  ,  0.7242  ,   0.6644  ,   0.6014 ,  &
              0.5348  ,  0.4686  ,   0.4028  ,   0.3412 ,  &
              0.2838  ,  0.2326  ,   0.1876  ,   0.1488 ,  &
              0.1164  ,  0.0892  ,   0.0672  ,   0.0488 ,  &
              0.0348  ,  0.0228  ,   0.0138  ,   0.0064 ,  &
              0.0000 /)

  ! *** cmam **************************************************************
  
  real, parameter  :: a_msc71(0:71) = (/   &
  0.0000000E+00, 0.9697381E-01, 0.1359287E+00, 0.1903444E+00, 0.2661948E+00, 0.3729337E+00, &
  0.5219892E+00, 0.7302384E+00, 0.1023965E+01, 0.1429865E+01, 0.1996367E+01, 0.2795496E+01, &
  0.3910756E+01, 0.5469573E+01, 0.7652655E+01, 0.1064398E+02, 0.1476761E+02, 0.2035368E+02, &
  0.2782417E+02, 0.3784165E+02, 0.5115196E+02, 0.6867558E+02, 0.9154660E+02, 0.1215257E+03, &
  0.1601752E+03, 0.2096810E+03, 0.2730070E+03, 0.3534054E+03, 0.4539333E+03, 0.5791536E+03, &
  0.7337209E+03, 0.9020061E+03, 0.1151447E+04, 0.1438797E+04, 0.1756787E+04, 0.2144013E+04, &
  0.2592025E+04, 0.3131841E+04, 0.3718082E+04, 0.4406751E+04, 0.5184419E+04, 0.6002543E+04, &
  0.6901672E+04, 0.7882596E+04, 0.8892869E+04, 0.9923001E+04, 0.1096051E+05, 0.1196177E+05, &
  0.1289793E+05, 0.1371837E+05, 0.1436708E+05, 0.1479985E+05, 0.1496916E+05, 0.1481930E+05, &
  0.1430370E+05, 0.1346313E+05, 0.1242185E+05, 0.1128402E+05, 0.1013239E+05, 0.9018101E+04, &
  0.7976487E+04, 0.7027693E+04, 0.6182180E+04, 0.5440920E+04, 0.4799065E+04, 0.4202355E+04, &
  0.3590442E+04, 0.2963440E+04, 0.2321461E+04, 0.1664614E+04, 0.1115094E+04, 0.0000000E+00  /)
                 
  real, parameter  :: b_msc71(0:71) = (/    &
  0.0000000E+00, 0.2585571E-09, 0.7040398E-09, 0.1535705E-08, 0.3013138E-08, 0.5588928E-08, &
  0.9977033E-08, 0.1738750E-07, 0.2995330E-07, 0.5067781E-07, 0.8520696E-07, 0.1431463E-06, &
  0.2392853E-06, 0.3989986E-06, 0.6646806E-06, 0.1095722E-05, 0.1800122E-05, 0.2924593E-05, &
  0.4696273E-05, 0.7484670E-05, 0.1182433E-04, 0.1850007E-04, 0.2865571E-04, 0.4416054E-04, &
  0.6735843E-04, 0.1018461E-03, 0.1529120E-03, 0.2279373E-03, 0.3362286E-03, 0.4919699E-03, &
  0.7133741E-03, 0.9869111E-03, 0.1466181E-02, 0.2084517E-02, 0.2893929E-02, 0.4006976E-02, &
  0.5507053E-02, 0.7581515E-02, 0.1018475E-01, 0.1375098E-01, 0.1841276E-01, 0.2415572E-01, &
  0.3156661E-01, 0.4113112E-01, 0.5287339E-01, 0.6738057E-01, 0.8526936E-01, 0.1069703E+00, &
  0.1332617E+00, 0.1646431E+00, 0.2011737E+00, 0.2437835E+00, 0.2934350E+00, 0.3496911E+00, &
  0.4130113E+00, 0.4790453E+00, 0.5416321E+00, 0.5992498E+00, 0.6510819E+00, 0.6971171E+00, &
  0.7375002E+00, 0.7725258E+00, 0.8025842E+00, 0.8281592E+00, 0.8497921E+00, 0.8694991E+00, &
  0.8893561E+00, 0.9093621E+00, 0.9295158E+00, 0.9498163E+00, 0.9665902E+00, 0.1000000E+01  /)

  ! --- interface -----------------------
  
  interface Init
    module procedure levi_Init
    module procedure levi_Init_levi
    module procedure levi_Init_key
    module procedure levi_Init_select
  end interface
  
  interface Done
    module procedure levi_Done
  end interface
  
  interface Check
    module procedure levi_Check
    module procedure levi_Check_3d
  end interface
  
  interface Compare
    module procedure levi_Compare
  end interface
  
  interface HPressure
    module procedure levi_HPressure_1d
    module procedure levi_HPressure_2d
    module procedure levi_HPressure_3d
  end interface
  
  interface FPressure
    module procedure levi_FPressure_1d
    module procedure levi_FPressure_2d
    module procedure levi_FPressure_3d
  end interface
  
  interface FillLevels
    module procedure levi_FillLevels_1d
    module procedure levi_FillLevels_2d
    module procedure levi_FillLevels_3d
  end interface
  
  interface FillLevelsParents
!FIXME    module procedure levi_FillParent_1d
!FIXME    module procedure levi_FillParent_2d
     module procedure levi_FillParent_3d
  end interface FillLevelsParents

CONTAINS

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    PRINTINFO
  !
  ! !DESCRIPTION: Basic print method for hybrid level objects.
  !\\
  !\\
  ! !INTERFACE:
  !
  subroutine PrintInfo( LevelInfo )
    !
    ! !INPUT PARAMETERS:
    !
    type(TLevelInfo), intent(in) ::  LevelInfo
    !
    ! !REVISION HISTORY: 
    !   15 Dec 2010 - P. Le Sager - written
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    write(gol,*)'============Hyb. Lev. Info ==================='; call goPr
    write(gol,'("name            :", a  )') LevelInfo%name; call goPr
    write(gol,'("up/down         :", a  )') LevelInfo%updo; call goPr
    write(gol,'("number of levels:", i3 )') LevelInfo%nlev; call goPr
    if (associated(LevelInfo%a)) then
       write(gol,'("Edge a0..aN:",F12.3,"...",F12.3)') LevelInfo%a(0),LevelInfo%a(LevelInfo%nlev); call goPr
    else
       write(gol,*)"Edge A coeff not associated yet"; call goPr
    endif
    if (associated(LevelInfo%b)) then
       write(gol,'("Edge b0..bN:",F12.3,"...",F12.3)') LevelInfo%b(0),LevelInfo%b(LevelInfo%nlev); call goPr
    else
       write(gol,*)"Edge B coeff not associated yet"; call goPr
    endif
    if (associated(LevelInfo%fa)) then
       write(gol,'("Center a1..aN:",F12.3,"...",F12.3)') LevelInfo%fa(1),LevelInfo%fa(LevelInfo%nlev); call goPr
    else
       write(gol,*)" Center A coeff not associated yet"; call goPr
    endif
    if (associated(LevelInfo%fb)) then
       write(gol,'("Center b1..bN:",F12.3,"...",F12.3)') LevelInfo%fb(1),LevelInfo%fb(LevelInfo%nlev); call goPr
    else
       write(gol,*)"Center B coeff not associated yet"; call goPr
    endif


    if (associated(LevelInfo%p0)) then
       write(gol,'("Edge p0..pN:",F12.3,"...",F12.3)') LevelInfo%p0(0),LevelInfo%p0(LevelInfo%nlev); call goPr
    else
       write(gol,*)"Edge Press not associated yet"; call goPr
    endif
    if (associated(LevelInfo%fp0)) then
       write(gol,'("Center p1..pN:",F12.3,"...",F12.3)') LevelInfo%fp0(1),LevelInfo%fp0(LevelInfo%nlev); call goPr
    else
       write(gol,*)" Center Press coeff not associated yet"; call goPr
    endif
    if (associated(LevelInfo%m0)) then
       write(gol,'("Box mass m1..mN:",F12.3,"...",F12.3)') LevelInfo%m0(1),LevelInfo%m0(LevelInfo%nlev); call goPr
    else
       write(gol,*)"Mass not associated yet"; call goPr
    endif
        
    ! -- parent levels
    write(gol,'("has a parent    :",L1)') LevelInfo%selection
    if (LevelInfo%selection) then
       write(gol,'("       with ",i3, " levels")') LevelInfo%nlev_parent; call goPr
    endif
    
    write(gol,*)'=============================================='; call goPr

  end subroutine PrintInfo
  !EOC

  ! =========================================================

  subroutine levi_Init( levi, nlev, a, b, status, name )

    use Binas, only : p0, grav
  
    ! --- in/out ---------------------------------------
  
    type(TLevelInfo), intent(out) ::  levi
    integer,          intent(in)  ::  nlev
    real,             intent(in)  ::  a(:)
    real,             intent(in)  ::  b(:)
    integer,          intent(out) ::  status
    
    character(len=*), intent(in), optional  ::  name
    
    ! --- const --------------------------------------
    
    character(len=*), parameter  ::  rname = mname//'/levi_Init'
    
    ! --- local --------------------------------------
    
    integer             ::  l
    
    ! --- begin ---------------------------------------
    
    ! check param
    if ( (size(a) /= nlev+1) .or. (size(b) /= nlev+1) ) then
      write (*,'("ERROR - invalid length of half level coeff:")')
      write (*,'("ERROR -   nlev     : ",i4)') nlev
      write (*,'("ERROR -   size(a)  : ",i4," (",i4,")")') size(a), nlev+1
      write (*,'("ERROR -   size(b)  : ",i4," (",i4,")")') size(b), nlev+1
      TRACEBACK; status=1; return
    end if
    
    ! store name ?
    if ( present(name) ) then
      levi%name = name
    else
      levi%name = 'levs'
    end if
    
    ! number of levels
    levi%nlev = nlev
    
    ! store half level hybride coeff
    allocate( levi%a(0:nlev) )
    allocate( levi%b(0:nlev) )
    levi%a = a
    levi%b = b
    
    ! hybride coeff increments:
    allocate( levi%da(nlev) )
    allocate( levi%db(nlev) )
    levi%da = levi%a(1:nlev) - levi%a(0:nlev-1)
    levi%db = levi%b(1:nlev) - levi%b(0:nlev-1)
    
    ! full level hybride coeff:
    allocate( levi%fa(nlev) )
    allocate( levi%fb(nlev) )
    levi%fa = ( levi%a(0:nlev-1) + levi%a(1:nlev) )/2.0
    levi%fb = ( levi%b(0:nlev-1) + levi%b(1:nlev) )/2.0
    
    ! half level coeff per full level:
    allocate( levi%fa_bounds(1:2,1:nlev) )
    allocate( levi%fb_bounds(1:2,1:nlev) )
    levi%fa_bounds(1,:) = levi%a(0:nlev-1)
    levi%fa_bounds(2,:) = levi%a(1:nlev  )
    levi%fb_bounds(1,:) = levi%b(0:nlev-1)
    levi%fb_bounds(2,:) = levi%b(1:nlev  )
    
    ! fill standard pressures:
    allocate( levi%p0(0:nlev) )
    levi%p0(0:nlev) = levi%a(0:nlev) + levi%b(0:nlev) * p0     ! Pa
    allocate( levi%fp0(1:nlev) )
    levi%fp0 = levi%fa + levi%fb * p0     ! Pa
    allocate( levi%m0(1:nlev) )
    levi%m0 = abs( levi%da + levi%db * p0 )/grav     ! kg air / m2
    
    ! upwards (decreasing pressure) or downwards (increasing pressure)
    if ( levi%p0(0) > levi%p0(nlev) ) then
      levi%updo = 'u'
    else
      levi%updo = 'd'
    end if
    
    ! no selection info
    levi%selection = .false.
    levi%updo_parent='-'
    nullify( levi%hlevs )
    nullify( levi%flevs )
    nullify( levi%a_parent )
    nullify( levi%b_parent )
    nullify( levi%da_parent )
    nullify( levi%db_parent )
    
    status = 0
      
  end subroutine levi_Init

  ! ***
    
  subroutine levi_Init_levi( levi, levi2, status )
  
    ! --- in/out ---------------------------------------
  
    type(TLevelInfo), intent(out) ::  levi
    type(TLevelInfo), intent(in)  ::  levi2
    integer,          intent(out) ::  status
    
    ! --- local --------------------------------------
    
    character(len=*), parameter  ::  rname = mname//'/levi_Init_levi'
    
    ! --- begin ---------------------------------------
    
    ! copy fields:
    call Init( levi, levi2%nlev, levi2%a, levi2%b, status, name=trim(levi2%name) )
    IF_NOTOK_RETURN(status=1)

    status = 0
      
  end subroutine levi_Init_levi

  ! ***
    
  subroutine levi_Init_key( levi, key, status )
  
    ! --- in/out ---------------------------------------
  
    type(TLevelInfo), intent(out) ::  levi
    character(len=*), intent(in)  ::  key
    integer, intent(out)          ::  status
    
    ! --- const --------------------------------------
    
    character(len=*), parameter  ::  rname = mname//'/levi_Init_key'
    
    ! --- local ---------------------------------------
    
    real, allocatable   ::  a(:), b(:)
    integer             ::  l

    ! --- begin ---------------------------------------

    select case ( key )

      case ( 'ec4' )
        call Init( levi, 4, a_ec4, b_ec4, status, name=trim(key) )
        IF_NOTOK_RETURN(status=1)

      case ( 'ec10' )
        call Init( levi, 10, a_ec10, b_ec10, status, name=trim(key) )
        IF_NOTOK_RETURN(status=1)

      case ( 'ec19' )
        call Init( levi, 19, a_ec19, b_ec19, status, name=trim(key) )
        IF_NOTOK_RETURN(status=1)

      case ( 'ec31' )
        call Init( levi, 31, a_ec31, b_ec31, status, name=trim(key) )
        IF_NOTOK_RETURN(status=1)

      case ( 'ec34' )
        call Init( levi, 34, a_ec34, b_ec34, status, name=trim(key) )
        IF_NOTOK_RETURN(status=1)

      case ( 'ec40' )
        call Init( levi, 40, a_ec40, b_ec40, status, name=trim(key) )
        IF_NOTOK_RETURN(status=1)

      case ( 'ec60' )
        call Init( levi, 60, a_ec60, b_ec60, status, name=trim(key) )
        IF_NOTOK_RETURN(status=1)

      case ( 'ec62' )
        call Init( levi, 62, a_ec62, b_ec62, status, name=trim(key) )
        IF_NOTOK_RETURN(status=1)

      case ( 'ec91' )
        call Init( levi, 91, a_ec91, b_ec91, status, name=trim(key) )
        IF_NOTOK_RETURN(status=1)

     case ( 'ec137' )
        call Init( levi, 137, a_ec137, b_ec137, status, name=trim(key) )
        IF_NOTOK_RETURN(status=1)
        
      case ( 'tm4' )
        allocate( a(0:4) )
        allocate( b(0:4) )
        do l = 0, 4
          a(l) = a_ec4(4-l)
          b(l) = b_ec4(4-l)
        end do
        call Init( levi, 4, a, b, status, name=trim(key) )
        IF_NOTOK_RETURN(status=1)
        deallocate( a )
        deallocate( b )

      case ( 'tm10' )
        allocate( a(0:10) )
        allocate( b(0:10) )
        do l = 0, 10
          a(l) = a_ec10(10-l)
          b(l) = b_ec10(10-l)
        end do
        call Init( levi, 10, a, b, status, name=trim(key) )
        IF_NOTOK_RETURN(status=1)
        deallocate( a )
        deallocate( b )

     case ( 'tm31' )
        allocate( a(0:31) )
        allocate( b(0:31) )
        do l = 0, 31
          a(l) = a_ec31(31-l)
          b(l) = b_ec31(31-l)
        end do
        call Init( levi, 31, a, b, status, name=trim(key) )
        IF_NOTOK_RETURN(status=1)
        deallocate( a )
        deallocate( b )

      case ( 'tm34' )
        allocate( a(0:34) )
        allocate( b(0:34) )
        do l = 0, 34
           a(l) = a_ec34(34-l)
           b(l) = b_ec34(34-l)
        end do
        call Init( levi, 34, a, b, status, name=trim(key) )
        IF_NOTOK_RETURN(status=1)
        deallocate( a )
        deallocate( b )

      case ( 'tm40' )
        allocate( a(0:40) )
        allocate( b(0:40) )
        do l = 0, 40
          a(l) = a_ec40(40-l)
          b(l) = b_ec40(40-l)
        end do
        call Init( levi, 40, a, b, status, name=trim(key) )
        IF_NOTOK_RETURN(status=1)
        deallocate( a )
        deallocate( b )

      case ( 'tm60' )
        allocate( a(0:60) )
        allocate( b(0:60) )
        do l = 0, 60
          a(l) = a_ec60(60-l)
          b(l) = b_ec60(60-l)
        end do
        call Init( levi, 60, a, b, status, name=trim(key) )
        IF_NOTOK_RETURN(status=1)
        deallocate( a )
        deallocate( b )

      case ( 'tm62' )
        allocate( a(0:62) )
        allocate( b(0:62) )
        do l = 0, 62
          a(l) = a_ec62(62-l)
          b(l) = b_ec62(62-l)
        end do
        call Init( levi, 62, a, b, status, name=trim(key) )
        IF_NOTOK_RETURN(status=1)
        deallocate( a )
        deallocate( b )

      case ( 'tm91' )
        allocate( a(0:91) )
        allocate( b(0:91) )
        do l = 0, 91
          a(l) = a_ec91(91-l)
          b(l) = b_ec91(91-l)
        end do
        call Init( levi, 91, a, b, status, name=trim(key) )
        IF_NOTOK_RETURN(status=1)
        deallocate( a )
        deallocate( b )

     case ( 'tm137' )
        allocate( a(0:137) )
        allocate( b(0:137) )
        do l = 0, 137
           a(l) = a_ec137(137-l)
           b(l) = b_ec137(137-l)
        end do
        call Init( levi, 137, a, b, status, name=trim(key) )
        IF_NOTOK_RETURN(status=1)
        deallocate( a )
        deallocate( b )

      case ( 'nc28' )
        call Init( levi, 28, a_nc28, b_nc28, status, name=trim(key) )
        IF_NOTOK_RETURN(status=1)

      case ( 'msc71' )
        call Init( levi, 71, a_msc71, b_msc71, status, name=trim(key) )
        IF_NOTOK_RETURN(status=1)

      case default
        write (*,'("ERROR - unknown key `",a,"`")') trim(key)
        TRACEBACK; status=1; return
    end select

    status = 0
      
  end subroutine levi_Init_key

  ! ***
    
  subroutine levi_Init_select( levi, levi_parent, hlev_select, status, name )
  
    ! --- in/out ---------------------------------------
  
    type(TLevelInfo), intent(out) :: levi
    type(TLevelInfo), intent(in)  :: levi_parent
    integer,          intent(in)  :: hlev_select(:)
    integer,          intent(out) :: status
    
    character(len=*), intent(in), optional  ::  name
    
    ! --- const --------------------------------------
    
    character(len=*), parameter  ::  rname = mname//'/levi_Init_select'
    
    ! --- local --------------------------------------
    
    integer             ::  nlev, l
    
    ! --- begin ---------------------------------------
    
    if ( (size(hlev_select) < 2) .or. (size(hlev_select) > levi_parent%nlev+1) ) then
      write(gol,'(" Strange length of array with selected levels:")')     ; call goErr
      write(gol,'("   expected : ",i4,", ..,",i4)') 2, levi_parent%nlev+1 ; call goErr
      write(gol,'("   found    : ",i4)') size(hlev_select)                ; call goErr
      TRACEBACK; status=1; return
    end if
    
    ! check range of values:
    if ( (minval(hlev_select) /= 0) .or. (maxval(hlev_select) /= levi_parent%nlev) ) then
      write(gol,'(" Invalid range of selected levels:")')                                    ; call goErr
      write(gol,'("   expected : ",i4,", ..,",i4)') 0, levi_parent%nlev                      ; call goErr
      write(gol,'("   found    : ",i4,", ..,",i4)') minval(hlev_select), maxval(hlev_select) ; call goErr
      TRACEBACK; status=1; return
    end if
    
    ! set number of full levels
    nlev = size(hlev_select) - 1
    
    ! copy coeff
    call Init( levi, nlev, levi_parent%a(hlev_select), levi_parent%b(hlev_select), status, name=name )
    IF_NOTOK_RETURN(status=1)

    ! store selection info
    levi%selection = .true.

    !  o levi_parent half levels indices corresponding to levi half levels
    allocate( levi%hlevs(0:nlev) )
    levi%hlevs = hlev_select
    
    !  o range of levi_parent full levels covered by a levi level. Note that
    !  bounds from hlev_select are (1:nlev+1).
    allocate( levi%flevs(nlev,2) )
    do l = 1, nlev
      levi%flevs(l,1) = hlev_select(l+1) + 1
      levi%flevs(l,2) = hlev_select(l)
   end do
   
    !  o original hybride coeff:
    levi%nlev_parent = levi_parent%nlev
    allocate( levi%a_parent(0:levi%nlev_parent) )
    allocate( levi%b_parent(0:levi%nlev_parent) )
    allocate( levi%da_parent(levi%nlev_parent) )
    allocate( levi%db_parent(levi%nlev_parent) )
    levi%a_parent = levi_parent%a
    levi%b_parent = levi_parent%b
    levi%da_parent = levi_parent%da
    levi%db_parent = levi_parent%db
    
    ! upwards (decreasing pressure) or downwards (increasing pressure)
    levi%updo_parent = levi_parent%updo

    status = 0
      
  end subroutine levi_Init_select

  ! ***
  
  subroutine levi_Done( levi, status )
  
    ! --- in/out ---------------------------------------
  
    type(TLevelInfo), intent(inout)     ::  levi
    integer, intent(out)                ::  status
    
    ! --- begin ---------------------------------------
    
    ! deallocate storage for hybride coeff
    deallocate( levi%a )
    deallocate( levi%b )
    
    ! deallocate storage for full level hybride coeff
    deallocate( levi%fa )
    deallocate( levi%fb )
    
    ! deallocate storage for half level coeff per full level:
    deallocate( levi%fa_bounds )
    deallocate( levi%fb_bounds )
    
    ! deallocate storage for standard pressures and mass:
    deallocate( levi%p0  )
    deallocate( levi%fp0 )
    deallocate( levi%m0  )

    ! deallocate storage for increment coeff
    deallocate( levi%da )
    deallocate( levi%db )
    
    ! selection info ?
    if ( levi%selection ) then
      deallocate( levi%hlevs )
      deallocate( levi%flevs )
      deallocate( levi%a_parent )
      deallocate( levi%b_parent )
      deallocate( levi%da_parent )
      deallocate( levi%db_parent )
    end if
    
    status = 0
      
  end subroutine levi_Done

  ! ***  
  
  subroutine levi_Check( levi, status )
  
    ! --- in/out ------------------------------
    
    type(TLevelInfo), intent(in)  ::  levi
    integer,          intent(out) ::  status
    
    ! --- const -------------------------------
    
    character(len=*), parameter  ::  rname = mname//'/levi_Check'
    
    ! --- begin -------------------------------
    
    if ( levi%nlev < 1 ) then
      write (*,'("ERROR - level info contains strange number of levels:")')
      write (*,'("ERROR -   levi%nlev  : ",i4)') levi%nlev
      TRACEBACK; status=1; return
    end if
    
    if ( (.not. associated(levi%a)) .or. (.not. associated(levi%b)) ) then
      write (*,'("ERROR - hybride coeffs in level info not allocated properly.")')
      TRACEBACK; status=1; return
    end if
    
    if ( (size(levi%a)/=levi%nlev+1) .or. (size(levi%b)/=levi%nlev+1) ) then
      write (*,'("ERROR - hybride coeffs in level info wrong size:")')
      write (*,'("ERROR -   nlev   : ",i4)') levi%nlev
      write (*,'("ERROR -   size a : ",i4)') size(levi%a)
      write (*,'("ERROR -   size b : ",i4)') size(levi%b)
      TRACEBACK; status=1; return
    end if
    
    if ( (.not. associated(levi%fa)) .or. (.not. associated(levi%fb)) ) then
      write (*,'("ERROR - f hybride coeffs in level info not allocated properly.")')
      TRACEBACK; status=1; return
    end if
    
    if ( (size(levi%fa)/=levi%nlev) .or. (size(levi%fb)/=levi%nlev) ) then
      write (*,'("ERROR - f hybride coeffs in level info wrong size:")')
      write (*,'("ERROR -   nlev    : ",i4)') levi%nlev
      write (*,'("ERROR -   size fa : ",i4)') size(levi%fa)
      write (*,'("ERROR -   size fb : ",i4)') size(levi%fb)
      TRACEBACK; status=1; return
    end if
    
    status = 0
    
  end subroutine levi_Check

  ! ***  
  
  subroutine levi_Check_3d( levi, ll, status )
  
    ! --- in/out ------------------------------
    
    type(TLevelInfo), intent(in)  ::  levi
    real,             intent(in)  ::  ll(:,:,:)
    integer,          intent(out) ::  status
    
    ! --- const -------------------------------
    
    character(len=*), parameter  ::  rname = mname//'/levi_Check_3d'
    
    ! --- begin -------------------------------
    
    call Check( levi, status )
    IF_NOTOK_RETURN(status=1)
    
    if ( size(ll,3) /= levi%nlev ) then
      write(gol,'(" Size of 3D grid does not match with level info")'); call goErr 
      write(gol,'("   size(ll,3) : ",i4)') size(ll,3)                 ; call goErr 
      write(gol,'("   levi%nlev  : ",i4)') levi%nlev                  ; call goErr 
      TRACEBACK; status=1; return
    end if
    
    status = 0
    
  end subroutine levi_Check_3d
  
  ! ##############################################################
  ! ###
  ! ### combining levels
  ! ###
  ! ##############################################################  
  !
  ! How levi and leviX compare?
  ! - exact match             : copy
  ! - levi is subset of leviX : combine
  ! - other case              : interpol
  !
  ! Also returns if levi and leviX are going in the same direction or not.
  !
  ! Note that case of "source grid is a subset of target grid" is not handled
  ! specifically, and thus is treated as 'interpol'. See FillLevelsParent
  ! routine, which takes care of that case.
  ! 
  subroutine levi_Compare( levi, leviX, verbose, fillkey, reverse, status )
  
    ! --- in/out ----------------------------------

    type(TLevelInfo),  intent(in)  ::  levi    ! target grid
    type(TLevelInfo),  intent(in)  ::  leviX   ! source grid
    integer,           intent(in)  ::  verbose
    character(len=10), intent(out) ::  fillkey ! copy, combine, interpol
    logical,           intent(out) ::  reverse ! same direction
    integer,           intent(out) ::  status
    
    ! --- const --------------------------------------
    
    character(len=*), parameter  ::  rname = mname//'/levi_Compare'
    
    ! --- local -----------------------------------
    
    integer                ::  l
    real, allocatable      ::  Xa_rev(:), Xb_rev(:)
    
    ! --- begin ------------------------------------
    
    ! reverse hybride coeff of source grid
    allocate( Xa_rev(0:leviX%nlev), Xb_rev(0:leviX%nlev) )
    do l = 0, leviX%nlev
      Xa_rev(l) = leviX%a(leviX%nlev-l)
      Xb_rev(l) = leviX%b(leviX%nlev-l)
    end do

    ! default
    fillkey = 'interpol'
    reverse  = levi%updo /= leviX%updo

    ! same number of levels ?
    if ( leviX%nlev == levi%nlev ) then

      ! exact match or reverse ...
      if      ( all(abs(leviX%a-levi%a)<0.1) .and. all(abs(leviX%b-levi%b)<0.1) ) then
        fillkey = 'copy'
        reverse = .false.
      else if ( all(abs(Xa_rev -levi%a)<0.1) .and. all(abs(Xb_rev -levi%b)<0.1) ) then
        fillkey = 'copy'
        reverse = .true.
      end if
      
    ! is target grid a subset of source grid? 
    else if ( levi%selection .and. (leviX%nlev == levi%nlev_parent) ) then

      ! exact match or reverse ...
      if      ( all(abs(leviX%a-levi%a_parent)<0.1) .and. all(abs(leviX%b-levi%b_parent)<0.1) ) then
        fillkey = 'combine'
        reverse = .false.
      else if ( all(abs( Xa_rev-levi%a_parent)<0.1) .and. all(abs( Xb_rev-levi%b_parent)<0.1) ) then
        fillkey = 'combine'
        reverse  = .true.
      end if

      !PLS The case of "source grid a subset of target grid" is not specifically
      !    covered (ie it uses 'interpol'). Could be added here, but decided
      !    for the alternate FillLevelsParents routine, which should be a bit faster.
       
    end if
    
    deallocate( Xa_rev, Xb_rev )
    
    status = 0
    
  end subroutine levi_Compare
  
  ! =====================================================================
  ! ===
  ! === pressure levels - returns pressure [Pa] at each half level of LevI
  ! ===
  ! =====================================================================
  
  subroutine levi_HPressure_1d( levi, sp, ph, status )
  
    ! --- in/out -------------------------------
    
    type(TLevelInfo), intent(in)  ::  levi
    real,             intent(in)  ::  sp        ! Pa
    real,             intent(out) ::  ph(:)     ! Pa
    integer,          intent(out) ::  status
    
    ! --- const --------------------------------------
    
    character(len=*), parameter  ::  rname = mname//'/levi_HPressure_1d'
    
    ! --- begin ---------------------------------
    
    ! check ...
    if ( size(ph) /= levi%nlev+1 ) then
      write (*,'("ERROR - shape of output grid does not match definition:")')
      write (*,'("ERROR -   half levels  : ",i4)') levi%nlev+1
      write (*,'("ERROR -   ph shape     : ",i4)') shape(ph)
      TRACEBACK; status=1; return
    end if
    
    ! half level pressure
    ph = levi%a + levi%b * sp
    
    status = 0
    
  end subroutine levi_HPressure_1d

  ! *
    
  subroutine levi_HPressure_2d( levi, sp, ph, status )
  
    ! --- in/out -------------------------------
    
    type(TLevelInfo), intent(in)  ::  levi
    real,             intent(in)  ::  sp(:)       ! Pa
    real,             intent(out) ::  ph(:,:)     ! Pa
    integer,          intent(out) ::  status
    
    ! --- const --------------------------------------
    
    character(len=*), parameter  ::  rname = mname//'/levi_HPressure_2d'
    
    ! --- local ---------------------------------
    
    integer     ::  i
    
    ! --- begin ---------------------------------
    
    if ( any( shape(ph) /= (/shape(sp),levi%nlev+1/) ) ) then
      write (*,'("ERROR - shape of output grid does not match definition:")')
      write (*,'("ERROR -   sp shape     : ",i6)') shape(sp)
      write (*,'("ERROR -   half levels  : ",i4)') levi%nlev+1
      write (*,'("ERROR -   ph shape     : ",i6,i4)') shape(ph)
      TRACEBACK; status=1; return
    end if
    
    ! half level pressure
    do i = 1, size(sp)
      ph(i,:) = levi%a + levi%b * sp(i)
    end do
    
    status = 0
    
  end subroutine levi_HPressure_2d

  ! *
    
  subroutine levi_HPressure_3d( levi, sp, ph, status )
  
    ! --- in/out -------------------------------
    
    type(TLevelInfo), intent(in)  ::  levi
    real,             intent(in)  ::  sp(:,:)       ! Pa
    real,             intent(out) ::  ph(:,:,:)     ! Pa
    integer,          intent(out) ::  status
    
    ! --- const --------------------------------------
    
    character(len=*), parameter  ::  rname = mname//'/levi_HPressure_3d'
    
    ! --- local ---------------------------------
    
    integer     ::  i, j
    
    ! --- begin ---------------------------------
    
    if ( any( shape(ph) /= (/shape(sp),levi%nlev+1/) ) ) then
      write (*,'("ERROR - shape of output grid does not match definition:")')
      write (*,'("ERROR -   sp shape     : ",2i6)') shape(sp)
      write (*,'("ERROR -   half levels  : ",i4)') levi%nlev+1
      write (*,'("ERROR -   ph shape     : ",2i6,i4)') shape(ph)
      TRACEBACK; status=1; return
    end if
    
    ! half level pressure
    do i = 1, size(sp,1)
      do j = 1, size(sp,2)
        ph(i,j,:) = levi%a + levi%b * sp(i,j)
      end do
    end do
    
    status = 0
    
  end subroutine levi_HPressure_3d
  
  ! ***

  subroutine levi_FPressure_1d( levi, sp, pf, status )
  
    ! --- in/out -------------------------------
    
    type(TLevelInfo), intent(in)  ::  levi
    real,             intent(in)  ::  sp        ! Pa
    real,             intent(out) ::  pf(:)     ! Pa
    integer,          intent(out) ::  status
    
    ! --- const --------------------------------------
    
    character(len=*), parameter  ::  rname = mname//'/levi_FPressure_1d'
    
    ! --- begin ---------------------------------
    
    if ( size(pf) /= levi%nlev ) then
      write(gol,'("ERROR - shape of output grid does not match definition:")'); call goErr 
      write(gol,'("ERROR -   full levels  : ",i4)') levi%nlev                 ; call goErr 
      write(gol,'("ERROR -   pf shape     : ",i4)') shape(pf)                 ; call goErr 
      TRACEBACK; status=1; return
    end if
    
    ! full level pressure
    pf = levi%fa + levi%fb * sp
    
    status = 0
    
  end subroutine levi_FPressure_1d

  ! *
    
  subroutine levi_FPressure_2d( levi, sp, pf, status )
  
    ! --- in/out -------------------------------
    
    type(TLevelInfo), intent(in)  :: levi
    real,             intent(in)  :: sp(:)       ! Pa
    real,             intent(out) :: pf(:,:)     ! Pa
    integer,          intent(out) :: status
    
    ! --- const --------------------------------------
    
    character(len=*), parameter  ::  rname = mname//'/levi_FPressure_2d'
    
    ! --- local ---------------------------------
    
    integer     ::  i
    
    ! --- begin ---------------------------------
    
    if ( any( shape(pf) /= (/shape(sp),levi%nlev/) ) ) then
      write(gol,'("ERROR - shape of output grid does not match definition:")'); call goErr 
      write(gol,'("ERROR -   sp shape     : ",i6)') shape(sp)                 ; call goErr 
      write(gol,'("ERROR -   full levels  : ",i4)') levi%nlev                 ; call goErr 
      write(gol,'("ERROR -   pf shape     : ",i6,i4)') shape(pf)              ; call goErr 
      TRACEBACK; status=1; return
    end if
    
    ! full level pressure
    do i = 1, size(sp)
      pf(i,:) = levi%fa + levi%fb * sp(i)
    end do
    
    status = 0
    
  end subroutine levi_FPressure_2d

  ! *
    
  subroutine levi_FPressure_3d( levi, sp, pf, status )
  
    ! --- in/out -------------------------------
    
    type(TLevelInfo), intent(in)  :: levi
    real,             intent(in)  :: sp(:,:)       ! Pa
    real,             intent(out) :: pf(:,:,:)     ! Pa
    integer,          intent(out) :: status
    
    ! --- const --------------------------------------
    
    character(len=*), parameter  :: rname = mname//'/levi_FPressure_3d'
    
    ! --- local ---------------------------------
    
    integer     ::  i, j
    
    ! --- begin ---------------------------------
    
    if ( any( shape(pf) /= (/shape(sp),levi%nlev/) ) ) then
      write(gol,'(" Shape of output grid does not match definition:")'); call goErr 
      write(gol,'("   sp shape     : ",2i6)') shape(sp)                ; call goErr 
      write(gol,'("   full levels  : ",i4)') levi%nlev                 ; call goErr 
      write(gol,'("   pf shape     : ",2i6,i4)') shape(pf)             ; call goErr 
      TRACEBACK; status=1; return
    end if
    
    ! full level pressure
    do i = 1, size(sp,1)
      do j = 1, size(sp,2)
        pf(i,j,:) = levi%fa + levi%fb * sp(i,j)
      end do
    end do
    
    status = 0
    
  end subroutine levi_FPressure_3d

  ! =====================================================================
  ! ===
  ! ===   REMAPPING
  ! ===
  ! =====================================================================   
  !
  ! Interpolate llX (defined on leviX) to ll (defined on levi).
  ! 
  ! How levels are combined/interpolated is specified by a key:
  ! 
  !     'bottom'     :  use closest-to-the-ground neighbor
  !     'top'        :  use closest-to-the-model-top neighbor
  !     'sum'        :  sum input values (use for mass [extensive])
  !     'aver'       :  unweighted average across levels
  !     'mass-aver'  :  air-mass weighted average (use for mixing ratio [intensive])
  !  
  ! Surface pressure field is used to compute mass in case of 'mass-aver' combination.
  !
  subroutine FillLevels_errmess(fillkey,nw,combine_key)

    character(len=*), intent(in)            :: fillkey
    character(len=*), intent(in)            :: nw
    character(len=*), intent(in), optional  :: combine_key
    
    if (present(combine_key)) then
       write (gol,'(" Combine key not supported:")')           ; call goErr
       write (gol,'("   combine key : ",a)') trim(combine_key) ; call goErr
       write (gol,'("   fill key    : ",a)') trim(fillkey)     ; call goErr
       write (gol,'("   nw key      : ",a)') trim(nw)          ; call goErr
    end if

    write (gol,'(" Fill key not supported:")')          ; call goErr
    write (gol,'("   fill key    : ",a)') trim(fillkey) ; call goErr
    write (gol,'("   nw key      : ",a)') trim(nw)      ; call goErr

  end subroutine FillLevels_errmess

  ! ***

  subroutine levi_FillLevels_1d( levi, ps, ll, leviX, llX, combine_key, status )
  
    ! --- in/out ----------------------------------

    type(TLevelInfo), intent(in)  :: levi
    real,             intent(in)  :: ps        ! Pa
    real,             intent(out) :: ll(:)
    type(TLevelInfo), intent(in)  :: leviX
    real,             intent(in)  :: llX(:)
    character(len=*), intent(in)  :: combine_key
    integer,          intent(out) :: status
    
    ! --- const --------------------------------------
    
    character(len=*), parameter  ::  rname = mname//'/levi_FillLevels_1d'
    
    ! --- local ------------------------------
    
    integer            ::  verbose
    character(len=10)  ::  fillkey
    logical            ::  reverse

    integer            ::  k, l
    integer            ::  le, nle, le1, le2
    real               ::  Xfp0, Xfp0_save
    
    ! --- begin -------------------------------
    
    ! correct number of levels ?
    if ( size(ll) /= levi%nlev ) then
      write(gol,'(" number of levels in output grid does not match definition:")') ; call goErr 
      write(gol,'("   ll  levels  : ",i3)') size(ll)                               ; call goErr 
      write(gol,'("   levi nlev   : ",i3)') levi%nlev                              ; call goErr 
      write(gol,'(" **NOTE** NO HALF-LEVELS case in 1D vertical remap!")')         ; call goErr 
      TRACEBACK; status=1; return
    end if
    
    ! necessary to combine or reverse ?
    verbose = 1
    call Compare( levi, leviX, verbose, fillkey, reverse, status )
    IF_NOTOK_RETURN(status=1)

    status = 1
    select case ( fillkey )

    case ( 'copy' )

       do l = 1, levi%nlev
          ! index of corresponding level in field 'llX' read from file:
          k = l
          if ( reverse ) k = levi%nlev + 1 - l
          ! copy level:
          ll(l) = llX(k)
       end do

    case ( 'combine' )
      
       ll = 0.0
       do l = 1, levi%nlev
          ! loop over range of parent levels:
          le1 = levi%flevs(l,1)
          le2 = levi%flevs(l,2)
          nle = le2 - le1 + 1
          Xfp0_save = -1.0
          do le = le1, le2
            ! index of corresponding level in field 'llX' read from file:
            k = le
            if ( reverse ) k = levi%nlev_parent + 1 - le
            ! standard pressure for level in llX :
            Xfp0 = leviX%fp0(k)
            ! based on combine key, add contribution of this level:
            select case ( combine_key )
            case ( 'sum' )
               ll(l) = ll(l) + llX(k)
            case ( 'aver' )
               ll(l) = ll(l) + llX(k)/(nle*1.0)
            case ( 'mass-aver' )
               !
               !  m = dp * A / g = |da+db*ps| * A/g
               !  ( X1*mX1 + X2*mX2 + .. ) / m  
               !  ~  ( X1*|daX1+dbX1*ps| + X2*|daX2+dbX2*ps| + .. ) / |da+db*ps|
               !
               ll(l) = ll(l) + llX(k)*abs(leviX%da(k)+leviX%db(k)*ps) / &
                                      abs( levi%da(l)+ levi%db(l)*ps)
            case ( 'bottom' )
               if ( (Xfp0_save < 0.0) .or. (Xfp0 >= Xfp0_save) ) then
                  ll(l) = llX(k)
                  Xfp0_save = Xfp0
               end if
            case ( 'top' )
               if ( (Xfp0_save < 0.0) .or. (Xfp0 <= Xfp0_save) ) then
                  ll(l) = llX(k)
                  Xfp0_save = Xfp0
               end if
            case default
               call FillLevels_errmess(fillkey,'n',combine_key)
               IF_NOTOK_RETURN(status=1)
            end select
         end do
      end do
      
   case default
      call FillLevels_errmess(fillkey,'n')
      IF_NOTOK_RETURN(status=1)
   end select
    
   status = 0

 end subroutine levi_FillLevels_1d

  ! ***

 subroutine levi_FillLevels_2d( levi, nw, ps, ll, &
                                 leviX, llX, combine_key, status )
  
    use Num, only : IntervalSum, Interp_Lin

    ! --- in/out ----------------------------------

    type(TLevelInfo), intent(in)  :: levi
    character(len=*), intent(in)  :: nw
    real,             intent(in)  :: ps(:)        ! Pa
    real,             intent(out) :: ll(:,:)
    type(TLevelInfo), intent(in)  :: leviX
    real,             intent(in)  :: llX(:,:)
    character(len=*), intent(in)  :: combine_key
    integer,          intent(out) :: status
    
    ! --- const --------------------------------------
    
    character(len=*), parameter  :: rname = mname//'/levi_FillLevels_2d'
    
    ! --- local ------------------------------
    
    integer            ::  verbose
    character(len=10)  ::  fillkey
    logical            ::  reverse

    integer            ::  j, k, l
    integer            ::  le, nle, le1, le2
    real               ::  Xfp0, Xfp0_save
    
    real, allocatable  ::  phalf(:)
    real, allocatable  ::  phalfX(:)
    real, allocatable  ::  dpX(:)
    real               ::  pfac
    integer            ::  ilast
    
    ! --- begin -------------------------------
  
    ! correct number of levels ?
    if ( ((nw == 'n') .and. (size(ll,2) /= levi%nlev  )) .or. &
         ((nw == 'w') .and. (size(ll,2) /= levi%nlev+1)) ) then
      write (*,'("ERROR - number of levels in output grid does not match definition:")')
      write (*,'("ERROR -   levi nlev   : ",i3)') levi%nlev
      write (*,'("ERROR -   nw          : ",a )') nw
      write (*,'("ERROR -   ll  levels  : ",i3)') size(ll,2)
      TRACEBACK; status=1; return
    end if
    
    ! correct number of levels ?
    if ( ((nw == 'n') .and. (size(llX,2) /= leviX%nlev  )) .or. &
         ((nw == 'w') .and. (size(llX,2) /= leviX%nlev+1)) ) then
      write (*,'("ERROR - number of levels in input grid does not match definition:")')
      write (*,'("ERROR -   leviX nlev   : ",i3)') leviX%nlev
      write (*,'("ERROR -   nw           : ",a )') nw
      write (*,'("ERROR -   llX  levels  : ",i3)') size(llX,2)
      TRACEBACK; status=1; return
    end if
    
    ! same horizontal grid sizes of ll and llX ?
    if ( (size(ll,1) /= size(llX,1)) ) then
      write (*,'("ERROR - horizontal size does not match:")')
      write (*,'("ERROR -   ll     : ",i3)') size(ll,1)
      write (*,'("ERROR -   llX    : ",i3)') size(llX,1)
      TRACEBACK; status=1; return
    end if

    ! necessary to combine or reverse ?
    verbose = 0
    call Compare( levi, leviX, verbose, fillkey, reverse, status )
    IF_NOTOK_RETURN(status=1)

    ! correct size of surface pressure field ?
    if ( (combine_key == 'mass-aver') .or. (fillkey == 'interpol') ) then
      if ( (size(ps) /= size(ll,1)) ) then
        write (*,'("ERROR - horizontal sizes do not match:")')
        write (*,'("ERROR -   ps          : ",i3)') size(ps)
        write (*,'("ERROR -   ll          : ",i3)') size(ll,1)
        write (*,'("ERROR -   combine_key : ",a)') combine_key
        write (*,'("ERROR -   fillkey     : ",a)') fillkey
        write (*,'("ERROR -   reverse     : ",l1)') reverse
        TRACEBACK; status=1; return
      end if
    end if
    
    ! levels or half-levels ?
    select case ( nw )
    case ( 'n' )
       !
       ! === levels ===================================
       !
       select case ( fillkey )
       case ( 'copy' )

          if ( reverse ) then
             do l = 1, levi%nlev
                ! copy corresponding level in field 'llX' read from file:
                ll(:,l) = llX(:,levi%nlev+1-l)
             end do
          else
             ll = llX
          end if

       case ( 'combine' )

          ll = 0.0
          do l = 1, levi%nlev
             ! loop over range of parent levels:
             le1 = levi%flevs(l,1)
             le2 = levi%flevs(l,2)
             nle = le2 - le1 + 1
             Xfp0_save = -1.0
             do le = le1, le2
                ! index of corresponding level in field 'llX' read from file:
                k = le
                if ( reverse ) k = levi%nlev_parent + 1 - le
                ! standard pressure for level in llX :
                Xfp0 = leviX%fp0(k)
                ! based on combine key, add contribution of this level:
                select case ( combine_key )
                case ( 'sum' )
                   ll(:,l) = ll(:,l) + llX(:,k)
                case ( 'aver' )
                   ll(:,l) = ll(:,l) + llX(:,k)/(nle*1.0)
                case ( 'mass-aver' )
                   !
                   !  m = dp * A / g = |da+db*ps| * A/g
                   !  ( X1*mX1 + X2*mX2 + .. ) / m  
                   !  ~  ( X1*|daX1+dbX1*ps| + X2*|daX2+dbX2*ps| + .. ) / |da+db*ps|
                   !
                   ll(:,l) = ll(:,l) + llX(:,k)*abs(leviX%da(k)+leviX%db(k)*ps) / &
                        abs( levi%da(l)+ levi%db(l)*ps)
                case ( 'bottom' )
                   if ( (Xfp0_save < 0.0) .or. (Xfp0 >= Xfp0_save) ) then
                      ll(:,l) = llX(:,k)
                      Xfp0_save = Xfp0
                   end if
                case ( 'top' )
                   if ( (Xfp0_save < 0.0) .or. (Xfp0 <= Xfp0_save) ) then
                      ll(:,l) = llX(:,k)
                      Xfp0_save = Xfp0
                   end if
                case default
                   call FillLevels_errmess(fillkey,nw,combine_key)
                   TRACEBACK; status=1; return
                end select
             end do
          end do

       case ( 'interpol' )

          allocate( phalf (0: levi%nlev) )
          allocate( phalfX(0:leviX%nlev) )
          allocate(    dpX(  leviX%nlev) )

          ! phalfX should be increasing
          pfac = 1.0
          if ( leviX%updo == 'u' ) pfac = -1.0

          select case ( combine_key )
             
          case ( 'sum' )

             do j = 1, size(ll,1)

                phalf  =  levi%a  +  levi%b  * ps(j)    ! Pa
                phalfX = leviX%a  + leviX%b  * ps(j)    ! Pa

                ilast = 1
                do l = 1, levi%nlev
                   ! take partly sums of f(i)*m(i) within [phalf(l-1),phalf(l)]
                   ! NOTE: if phalf(l) < phalf(l-1), the result of IntervalSum is negative
                   call IntervalSum( phalfX*pfac, llX(j,:), &
                        phalf(l-1)*pfac, phalf(l)*pfac, &
                        ll(j,l), ilast, status )
                   if (status/=0) then; 
                      write (*,'("ERROR - from IntervalSum (combine key `sum`)")')
                      write (*,'("ERROR -    j    : ",2i4   )')  j
                      write (*,'("ERROR -   ps      : ", f16.4)') ps(j)
                      write (*,'("ERROR -   leviX%a : ",es11.4)') leviX%a
                      write (*,'("ERROR -   leviX%b : ",es11.4)') leviX%b
                      TRACEBACK; status=1; return
                   end if
                   ll(j,l) = ll(j,l) * sign(1.0,(phalf(l)-phalf(l-1))*pfac)
                end do

             end do

          case ( 'mass-aver' )

             do j = 1, size(ll,1)

                phalf  =  levi%a  +  levi%b  * ps(j)    ! Pa
                phalfX = leviX%a  + leviX%b  * ps(j)    ! Pa
                dpX = abs(leviX%da + leviX%db * ps(j))    ! Pa

                ilast = 1
                do l = 1, levi%nlev
                   ! take partly sums of f(i)*m(i) within [phalf(l-1),phalf(l)]
                   ! NOTE: if phalf(l) < phalf(l-1), the result of IntervalSum is negative
                   call IntervalSum( phalfX*pfac, llX(j,:)*dpX, &
                                     phalf(l-1)*pfac, phalf(l)*pfac, &
                                     ll(j,l), ilast, status )
                   IF_NOTOK_RETURN(status=1)
                   ll(j,l) = ll(j,l) / (phalf(l)-phalf(l-1))*pfac
                end do

             end do

          case ( 'top' )
             
             do j = 1, size(ll,1)

                phalf  =  levi%a  +  levi%b  * ps(j)    ! Pa, 0..
                phalfX = leviX%a  + leviX%b  * ps(j)    ! Pa, 0..

                do l = 1, levi%nlev
                   ! interpolate to top half level pressure:
                    call Interp_Lin( phalfX(1:leviX%nlev)*pfac, llX(j,:), &
                                     phalf (l)           *pfac,  ll(j,l), ilast, status )
                    IF_NOTOK_RETURN(status=1)
                end do

             end do

          case ( 'bottom' )

             do j = 1, size(ll,1)

                phalf  =  levi%a  +  levi%b  * ps(j)    ! Pa, 0..
                phalfX = leviX%a  + leviX%b  * ps(j)    ! Pa, 0..

                do l = 1, levi%nlev
                   ! interpolate to bottom half level pressure:
                   call Interp_Lin( phalfX(0:leviX%nlev-1)*pfac, llX(j,:), &
                                    phalf (l-1)           *pfac,  ll(j,l), ilast, status )
                   IF_NOTOK_RETURN(status=1)
                end do

             end do

          case default
             call FillLevels_errmess(fillkey,nw,combine_key)
             TRACEBACK; status=1; return
          end select

          deallocate( phalf  )
          deallocate( phalfX )
          deallocate(    dpX )

       case default
          call FillLevels_errmess(fillkey,nw)
          TRACEBACK; status=1; return
       end select
        

    case ( 'w' )
       !
       ! === half levels ========================================
       !
       select case ( fillkey )

       case ( 'copy' )

          if ( reverse ) then
             do l = 1, levi%nlev+1
                ll(:,l) = llX(:,levi%nlev+2-l)
             end do
          else
             ll = llX
          end if

       case ( 'combine' )
          
          ll = 0.0

          do l = 0, levi%nlev
             ! note: k=0,..,nlev
             k = levi%hlevs(l)
             if ( reverse ) k = levi%nlev_parent - k
             ! copy:
             ll(:,l+1) = llX(:,k+1)
          end do

       case ( 'interpol' )
          
          allocate( phalf (0: levi%nlev) )
          allocate( phalfX(0:leviX%nlev) )
          allocate(    dpX(  leviX%nlev) )

          ! phalfX should be increasing
          pfac = 1.0
          if ( leviX%updo == 'u' ) pfac = -1.0

          select case ( combine_key )

          case ( 'top', 'bottom' )

             do j = 1, size(ll,1)

                phalf  =  levi%a  +  levi%b  * ps(j)    ! Pa
                phalfX = leviX%a  + leviX%b  * ps(j)    ! Pa

                do l = 0, levi%nlev
                   ! interpolate to half level pressure:
                   call Interp_Lin( phalfX   *pfac, llX(j,:  ), &
                                    phalf (l)*pfac,  ll(j,l+1), ilast, status )
                   IF_NOTOK_RETURN(status=1)
                end do

             end do

          case default
             call FillLevels_errmess(fillkey,nw,combine_key)
             TRACEBACK; status=1; return
          end select

          deallocate( phalf  )
          deallocate( phalfX )
          deallocate(    dpX )

       case default
          call FillLevels_errmess(fillkey,nw)
          TRACEBACK; status=1; return
       end select
        
    case default
       write (*,'("ERROR - nw key `",a,"` not supported.")') trim(nw)
       TRACEBACK; status=1; return
    end select

    status = 0

  end subroutine levi_FillLevels_2d
  
  ! ***

  subroutine levi_FillLevels_3d( levi, nw, ps, ll, &
                                 leviX, llX, combine_key, status )
  
    use Num, only : IntervalSum, Interp_Lin

    ! --- in/out ----------------------------------

    type(TLevelInfo), intent(in)  ::  levi
    character(len=*), intent(in)  ::  nw
    real,             intent(in)  ::  ps(:,:)        ! Pa
    real,             intent(out) ::  ll(:,:,:)
    type(TLevelInfo), intent(in)  ::  leviX
    real,             intent(in)  ::  llX(:,:,:)
    character(len=*), intent(in)  ::  combine_key
    integer,          intent(out) ::  status
    
    ! --- const --------------------------------------
    
    character(len=*), parameter  ::  rname = mname//'/levi_FillLevels_3d'
    
    ! --- local ------------------------------
    
    integer            ::  verbose
    character(len=10)  ::  fillkey
    logical            ::  reverse

    integer            ::  i, j, k, l
    integer            ::  le, nle, le1, le2
    real               ::  Xfp0, Xfp0_save
    
    real, allocatable  ::  phalf(:)
    real, allocatable  ::  phalfX(:)
    real, allocatable  ::  dpX(:)
    real               ::  pfac
    integer            ::  ilast
    
    ! --- begin -------------------------------
  
    ! correct number of levels OUT?
    if ( ((nw == 'n') .and. (size(ll,3) /= levi%nlev  )) .or. &
         ((nw == 'w') .and. (size(ll,3) /= levi%nlev+1)) ) then
      write(gol,'(" number of levels in output grid does not match definition:")') ; call goErr 
      write(gol,'("   levi nlev   : ",i3)') levi%nlev                              ; call goErr 
      write(gol,'("   nw          : ",a )') nw                                     ; call goErr 
      write(gol,'("   ll  levels  : ",i3)') size(ll,3)                             ; call goErr 
      TRACEBACK; status=1; return
    end if
    
    ! correct number of levels IN?
    if ( ((nw == 'n') .and. (size(llX,3) /= leviX%nlev  )) .or. &
         ((nw == 'w') .and. (size(llX,3) /= leviX%nlev+1)) ) then
      write(gol,'(" number of levels in input grid does not match definition:")') ; call goErr 
      write(gol,'("   leviX nlev   : ",i3)') leviX%nlev                           ; call goErr 
      write(gol,'("   nw           : ",a )') nw                                   ; call goErr 
      write(gol,'("   llX  levels  : ",i3)') size(llX,3)                          ; call goErr 
      TRACEBACK; status=1; return
    end if
    
    ! same horizontal grid sizes of ll and llX ?
    if ( (size(ll,1) /= size(llX,1)) .or. (size(ll,2) /= size(llX,2)) ) then
      write(gol,'(" horizontal sizes do not match:")')                  ; call goErr 
      write(gol,'("   ll     : ",i3," x ",i3)') size(ll,1), size(ll,2)  ; call goErr 
      write(gol,'("   llX    : ",i3," x ",i3)') size(llX,1), size(llX,2); call goErr 
      TRACEBACK; status=1; return
    end if

    ! necessary to combine or reverse ?
    verbose = 0
    call Compare( levi, leviX, verbose, fillkey, reverse, status )
    IF_NOTOK_RETURN(status=1)

    ! correct size of surface pressure field ?
    if ( (combine_key == 'mass-aver') .or. (fillkey == 'interpol') ) then
      if ( (size(ps,1) /= size(ll,1)) .or. (size(ps,2) /= size(ll,2)) ) then
        write(gol,'(" horizontal sizes do not match:")')                        ; call goErr 
        write(gol,'("   ps          : ",i3," x ",i3)') size(ps ,1), size(ps ,2) ; call goErr 
        write(gol,'("   ll          : ",i3," x ",i3)') size(ll ,1), size(ll ,2) ; call goErr 
        write(gol,'("   combine_key : ",a)') combine_key                        ; call goErr 
        write(gol,'("   fillkey     : ",a)') fillkey                            ; call goErr
        write(gol,'("   reverse     : ",l1)') reverse                           ; call goErr
        TRACEBACK; status=1; return
      end if
    end if
    
    select case ( nw )
    case ( 'n' )
       !
       ! === full levels ===================================
       !
       select case ( fillkey )
       case ( 'copy' )

          if ( reverse ) then
             do l = 1, levi%nlev
                ll(:,:,l) = llX(:,:,levi%nlev+1-l)
             end do
          else
             ll = llX
          end if

       case ( 'combine' )

          ll = 0.0
          do l = 1, levi%nlev
             ! loop over range of parent levels:
             le1 = levi%flevs(l,1)
             le2 = levi%flevs(l,2)
             nle = le2 - le1 + 1
             Xfp0_save = -1.0
             do le = le1, le2
                ! index of corresponding level in field 'llX' read from file:
                k = le
                if ( reverse ) k = levi%nlev_parent + 1 - le
                ! standard pressure for level in llX :
                Xfp0 = leviX%fp0(k)

                ! based on combine key, add contribution of this level:
                select case ( combine_key )
                case ( 'sum' )
                   ll(:,:,l) = ll(:,:,l) + llX(:,:,k)
                case ( 'aver' )
                   ll(:,:,l) = ll(:,:,l) + llX(:,:,k)/(nle*1.0)
                case ( 'mass-aver', 'height-ave')
                   !
                   !  m = dp * A / g = |da+db*ps| * A/g
                   !  ( X1*mX1 + X2*mX2 + .. ) / m  
                   !  ~  ( X1*|daX1+dbX1*ps| + X2*|daX2+dbX2*ps| + .. ) / |da+db*ps|
                   !
                   ll(:,:,l) = ll(:,:,l) + llX(:,:,k)*abs(leviX%da(k)+leviX%db(k)*ps) / &
                        abs( levi%da(l)+ levi%db(l)*ps)
                case ( 'bottom' )
                   if ( (Xfp0_save < 0.0) .or. (Xfp0 >= Xfp0_save) ) then
                      ll(:,:,l) = llX(:,:,k)
                      Xfp0_save = Xfp0
                   end if
                case ( 'top' )
                   if ( (Xfp0_save < 0.0) .or. (Xfp0 <= Xfp0_save) ) then
                      ll(:,:,l) = llX(:,:,k)
                      Xfp0_save = Xfp0
                   end if
                case default
                   call FillLevels_errmess(fillkey,nw,combine_key)
                   TRACEBACK; status=1; return
                end select
             end do
          end do

       case ( 'interpol' )

          allocate( phalf (0: levi%nlev) )
          allocate( phalfX(0:leviX%nlev) )
          allocate(    dpX(  leviX%nlev) )
          
          ! phalfX should be increasing
          pfac = 1.0
          if ( leviX%updo == 'u' ) pfac = -1.0

          select case ( combine_key )             
          case ( 'sum' )

             do j = 1, size(ll,2)
             do i = 1, size(ll,1)

                phalf  =  levi%a  +  levi%b  * ps(i,j)    ! Pa
                phalfX = leviX%a  + leviX%b  * ps(i,j)    ! Pa

                ilast = 1
                do l = 1, levi%nlev
                   ! take partly sums of f(i)*m(i) within [phalf(l-1),phalf(l)]
                   ! NOTE: if phalf(l) < phalf(l-1), the result of IntervalSum is negative
                   call IntervalSum( phalfX*pfac, llX(i,j,:), &
                                     phalf(l-1)*pfac, phalf(l)*pfac, &
                                     ll(i,j,l), ilast, status )
                   if (status/=0) then; 
                      write(gol,'(" from IntervalSum (combine key `sum`)")') ; call goErr 
                      write(gol,'("   i, j    : ",2i4   )') i, j             ; call goErr 
                      write(gol,'("   ps      : ", f16.4)') ps(i,j)          ; call goErr 
                      write(gol,'("   leviX%a : ",es11.4)') leviX%a          ; call goErr 
                      write(gol,'("   leviX%b : ",es11.4)') leviX%b          ; call goErr
                      TRACEBACK; status=1; return
                   end if
                   ll(i,j,l) = ll(i,j,l) * sign(1.0,(phalf(l)-phalf(l-1))*pfac)
                end do

             end do
             end do

          case ( 'mass-aver', 'height-ave' )

             do j = 1, size(ll,2)
             do i = 1, size(ll,1)

                phalf  =  levi%a  +  levi%b  * ps(i,j)    ! Pa
                phalfX = leviX%a  + leviX%b  * ps(i,j)    ! Pa
                dpX = abs(leviX%da + leviX%db * ps(i,j))    ! Pa

                ilast = 1
                do l = 1, levi%nlev
                   ! take partly sums of f(i)*m(i) within [phalf(l-1),phalf(l)]
                   ! NOTE: if phalf(l) < phalf(l-1), the result of IntervalSum is negative
                   call IntervalSum( phalfX*pfac, llX(i,j,:)*dpX, &
                                     phalf(l-1)*pfac, phalf(l)*pfac, &
                                     ll(i,j,l), ilast, status )
                   IF_NOTOK_RETURN(status=1)
                   ll(i,j,l) = ll(i,j,l) / (phalf(l)-phalf(l-1))*pfac
                end do

             end do
             end do

          case ( 'top' )

             do j = 1, size(ll,2)
             do i = 1, size(ll,1)

                phalf  =  levi%a  +  levi%b  * ps(i,j)    ! Pa, 0..
                phalfX = leviX%a  + leviX%b  * ps(i,j)    ! Pa, 0..

                do l = 1, levi%nlev
                   ! interpolate to top half level pressure:
                   call Interp_Lin( phalfX(1:leviX%nlev)*pfac, llX(i,j,:), &
                                    phalf (l)           *pfac,  ll(i,j,l), ilast, status )
                   IF_NOTOK_RETURN(status=1)
                end do

             end do
             end do

          case ( 'bottom' )

             do j = 1, size(ll,2)
             do i = 1, size(ll,1)

                phalf  =  levi%a  +  levi%b  * ps(i,j)    ! Pa, 0..
                phalfX = leviX%a  + leviX%b  * ps(i,j)    ! Pa, 0..

                do l = 1, levi%nlev
                   ! interpolate to bottom half level pressure:
                   call Interp_Lin( phalfX(0:leviX%nlev-1)*pfac, llX(i,j,:), &
                                    phalf (l-1)           *pfac,  ll(i,j,l), ilast, status )
                   IF_NOTOK_RETURN(status=1)
                end do

             end do
             end do

          case default
             call FillLevels_errmess(fillkey,nw,combine_key)
             TRACEBACK; status=1; return
          end select

          deallocate( phalf  )
          deallocate( phalfX )
          deallocate(    dpX )

       case default
          call FillLevels_errmess(fillkey,nw)
          TRACEBACK; status=1; return
       end select
        
    case ( 'w' )
       !
       ! === half levels ========================================
       !
       select case ( fillkey )
       case ( 'copy' )

          if ( reverse ) then
             do l = 1, levi%nlev+1
                ! copy corresponding level in field 'llX' read from file:
                ll(:,:,l) = llX(:,:,levi%nlev+2-l)
             end do
          else
             ll = llX
          end if

       case ( 'combine' )

          ll = 0.0

          do l = 0, levi%nlev
             ! index of corresponding level in field 'llX' read from file:
             ! note: k=0,..,nlev
             k = levi%hlevs(l)
             if ( reverse ) k = levi%nlev_parent - k
             ! copy:
             ll(:,:,l+1) = llX(:,:,k+1)
          end do

       case ( 'interpol' )

          allocate( phalf (0: levi%nlev) )
          allocate( phalfX(0:leviX%nlev) )
          allocate(    dpX(  leviX%nlev) )

          ! phalfX should be increasing
          pfac = 1.0
          if ( leviX%updo == 'u' ) pfac = -1.0

          select case ( combine_key )

          case ( 'top', 'bottom' )

             do j = 1, size(ll,2)
             do i = 1, size(ll,1)

                phalf  =  levi%a  +  levi%b  * ps(i,j)    ! Pa
                phalfX = leviX%a  + leviX%b  * ps(i,j)    ! Pa

                do l = 0, levi%nlev
                   ! interpolate to half level pressure:
                   call Interp_Lin( phalfX   *pfac, llX(i,j,:  ), &
                                    phalf (l)*pfac,  ll(i,j,l+1), ilast, status )
                   IF_NOTOK_RETURN(status=1)
                end do
             end do
             end do

          case default
             call FillLevels_errmess(fillkey,nw,combine_key)
             TRACEBACK; status=1; return
          end select

          deallocate( phalf  )
          deallocate( phalfX )
          deallocate(    dpX )

       case default
          call FillLevels_errmess(fillkey,nw)
          TRACEBACK; status=1; return
       end select
        
    case default
       write (*,'("ERROR - nw key `",a,"` not supported.")') trim(nw)
       TRACEBACK; status=1; return
    end select

    status = 0

  end subroutine levi_FillLevels_3d

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    LEVI_FILLPARENT_3D
  !
  ! !DESCRIPTION: remap one array into another one, their 3rd dimension being
  !                along either LevInfo levels or its parent levels. 
  !\\
  !\\
  ! !INTERFACE:
  !
  subroutine levi_FillParent_3d( levi, nw, ArrIn, combine_key, ArrOut, status, ps)
    !
    ! !INPUT PARAMETERS:
    !
    type(TLevelInfo), intent(in) :: levi
    character(len=*), intent(in) :: nw
    real,             intent(in) :: ArrIn(:,:,:)
    character(len=*), intent(in) :: combine_key
    !
    real,   optional, intent(in) :: ps(:,:)        ! Pa
    !
    ! !OUTPUT PARAMETERS:
    !
    real,            intent(out) :: ArrOut(:,:,:)
    integer,         intent(out) :: status
    !
    ! !REMARKS:
    ! (1) developed to account for child-to-parent cases, which were handled always as 'interpol' in FillLevels.
    ! (2) by design, we are always in 'combine' filling mode. No 'interpolation' or 'copy'. No need to call levi_compare.
    ! (3) Arrays follow the order ('u' or 'd') of the levels they are attached to.
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    character(len=*), parameter  ::  rname = mname//'/levi_FillParent_3d'

    logical            ::  reverse, to_parent, to_child

    integer            ::  i, j, k, l
    integer            ::  le, nle, le1, le2
    real               ::  Xfp0, Xfp0_save

    real, allocatable  ::  phalf(:)
    real, allocatable  ::  phalfX(:)
    real, allocatable  ::  dpX(:)
    real               ::  pfac
    integer            ::  ilast

    ! --- begin -------------------------------

    if (.not.levi%selection) then
       write(gol,*) "Cannot use FillLevelsParents: no parent available in levI"; call goErr
       TRACEBACK; status=1; return
    endif

    select case ( nw )
    case ('n')
      to_parent = (size(ArrOut,3) == levi%nlev_parent).and.(size(ArrIn,3) == levi%nlev )
      to_child  = (size(ArrOut,3) == levi%nlev)       .and.(size(ArrIn,3) == levi%nlev_parent)
    case ('w')
      write(gol,*) "FillLevelsParents: not implemented yet for half-levels"; call goErr
      TRACEBACK; status=1; return
    case default
      write(gol,*) "FillLevelsParents: input nw must be 'n' and not "//trim(nw); call goErr
      TRACEBACK; status=1; return
    end select

    if (.not.(to_parent.or.to_child)) then
       write(gol,*) "Cannot use FillLevelsParents: array sizes do not match number of levels"; call goErr
       write(gol,'("   ArrOut nlevels: ",i3)') size(ArrOut,3)   ; call goErr
       write(gol,'("   ArrIn  nlevels: ",i3)') size(ArrIn,3)    ; call goErr
       write(gol,'("   child  nlevels: ",i3)') levi%nlev        ; call goErr
       write(gol,'("   parent nlevels: ",i3)') levi%nlev_parent ; call goErr
       TRACEBACK; status=1; return
    endif

    ! same horizontal grid sizes of ArrOut and ArrIn ?
    if ( (size(ArrOut,1) /= size(ArrIn,1)) .or. (size(ArrOut,2) /= size(ArrIn,2)) ) then
       write(gol,'(" Horizontal sizes do not match:")') ; call goErr
       write(gol,'("   ArrOut : ",i3," x ",i3)') size(ArrOut,1), size(ArrOut,2) ; call goErr 
       write(gol,'("   ArrIn  : ",i3," x ",i3)') size(ArrIn,1),  size(ArrIn,2)  ; call goErr 
       TRACEBACK; status=1; return
    end if

    ! check cases when surface pressure needed
    if ( ((combine_key == 'mass-aver') .and. to_child  ) .or. &
         ((combine_key == 'sum')       .and. to_parent ) ) then
       if (.not. present(ps)) then
          write(gol,*)" surface pressure required"; call goErr
          TRACEBACK; status=1; return
       endif
       if ( (size(ps,1) /= size(ArrOut,1)) .or. (size(ps,2) /= size(ArrOut,2)) ) then
          write(gol,'(" Horizontal sizes do not match:")')                           ; call goErr
          write(gol,'("   ps     : ",i3," x ",i3)') size(ps ,1),     size(ps ,2)     ; call goErr
          write(gol,'("   ArrOut : ",i3," x ",i3)') size(ArrOut ,1), size(ArrOut ,2) ; call goErr
          TRACEBACK; status=1; return
       end if
    end if

    ! reverse = levi%updo /= levi%updo_parent ! NOT USED YET
    ArrOut = 0.
    
    select case ( nw )
    case ( 'n' )
       !
       ! === full levels ===================================
       !
       select case ( combine_key )
       case ( 'mass-aver' )     ! Intensive w/r/t pressure/air_mass (eg, mass-mixing-ratio)
       
          if(to_parent) then
             do l = 1, levi%nlev
                le1 = levi%flevs(l,1)
                le2 = levi%flevs(l,2)
                do le = le1, le2
                   ArrOut(:,:,le)=ArrIn(:,:,l)
                enddo
             enddo
          else if(to_child) then  ! same as combine/mass-aver from FillLevels
             do l = 1, levi%nlev
                le1 = levi%flevs(l,1)
                le2 = levi%flevs(l,2)                
                do le = le1, le2
                   ArrOut(:,:,l) = ArrOut(:,:,l) + ArrIn(:,:,le)*abs(levi%da_parent(le)+levi%db_parent(le)*ps)
                enddo
                ArrOut(:,:,l) = ArrOut(:,:,l) / abs( levi%da(l)+ levi%db(l)*ps )
             enddo
          endif
       
       case ( 'sum' )           ! Extensive (eg, mass)
          
          if(to_parent) then
             do l = 1, levi%nlev
                le1 = levi%flevs(l,1)
                le2 = levi%flevs(l,2)
                do le = le1, le2
                   ArrOut(:,:,le) = ArrIn(:,:,l)*abs(levi%da_parent(le)+levi%db_parent(le)*ps) / &
                        abs( levi%da(l)+ levi%db(l)*ps )
                enddo
             enddo
          else if(to_child) then ! same as combine/sum from FillLevels
             do l = 1, levi%nlev
                le1 = levi%flevs(l,1)
                le2 = levi%flevs(l,2)                
                ArrOut(:,:,l) = sum(ArrIn(:,:,le1:le2),dim=3)
             enddo
          endif

       case ( 'aver' )          ! Nominal (dimensionless, eg, level number)

          if(to_parent) then
             do l = 1, levi%nlev
                le1 = levi%flevs(l,1)
                le2 = levi%flevs(l,2)
                do le = le1, le2
                   ArrOut(:,:,le)=ArrIn(:,:,l)
                enddo
             enddo
          else if(to_child) then ! same as combine/aver from FillLevels
             do l = 1, levi%nlev
                le1 = levi%flevs(l,1)
                le2 = levi%flevs(l,2)
                ArrOut(:,:,l) = sum(ArrIn(:,:,le1:le2),dim=3)/((le2 - le1 + 1)*1.0)
             enddo
          endif

       case ( 'height-ave')     ! Intensive w/r/t height

          if(to_parent) then
             do l = 1, levi%nlev
                le1 = levi%flevs(l,1)
                le2 = levi%flevs(l,2)
                do le = le1, le2
                   ArrOut(:,:,le)=ArrIn(:,:,l)
                enddo
             enddo
          else if(to_child) then
             ! Layer Height is available from gph for child levels only. See compute_gph and phys_gph.F90.
             ! FIXME - dH (delta Height on the parent grid is not known. Requires temperature and humidity.
             ! do l = 1, levi%nlev
             !    le1 = levi%flevs(l,1)
             !    le2 = levi%flevs(l,2)
             !    do le = le1, le2
             !       ArrOut(:,:,l)=ArrOut(:,:,l)+ArrIn(:,:,le)*dH(le)
             !    enddo
             !    ArrOut(:,:,l) = ArrOut(:,:,l)  /  ( gph(:,:,l+1)-gph(:,:,l) )
             ! enddo
             call FillLevels_errmess('combine',nw,combine_key)
             TRACEBACK; status=1; return
          end if
          
       case default
          call FillLevels_errmess('combine',nw,combine_key)
          TRACEBACK; status=1; return
       end select       
    
    ! case ( 'w' )
    !    !
    !    ! === half levels ========================================
    !    !
    case default
       write(gol,'("  nw key `",a,"` not supported.")') trim(nw) ; call goErr
       TRACEBACK; status=1; return
    end select
    
    status = 0

  END SUBROUTINE LEVI_FILLPARENT_3D
  !EOC

END MODULE GRID_TYPE_HYB