!
!  Module and data types:
!
!    use grid_type_hyb
!
!    type(TLevelInfo)    ::  levi
!    type(TLevelInfo)    ::  levi_parent
!
!    integer              ::  status
!
!  Initialize in one of the following 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
!       '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 )
!
!  Copy 3D array llX into ll, 
!  How levels are combined is specified by a key:
!   'bottom'     :  use bottom value (most close to the ground)
!   'top'        :  use top value (most close to the model top)
!   'sum'        :  sum values
!   'aver'       :  average of all levels
!   'mass-aver'  : mass weighted average
!  Eventually reverse.
!
!      call FillLevels( levi, ps, ll, leviX, llX, combine_key, status )
!
! 
!

module grid_type_hyb

  implicit none
  
  ! --- in/out --------------------------
  
  private
  
  public  ::  TLevelInfo
  public  ::  Init, Done
  public  ::  Check
  public  ::  Compare
  public  ::  HPressure, FPressure
  public  ::  FillLevels
  public  ::  PrintInfo
 
  ! --- const ---------------------------------------
  
  character(len=*), parameter  ::  mname = 'grid_type_hyb'
  
  
  ! --- types ---------------------------
  
  type TLevelInfo
    ! name used for messages:
    character(len=32)                   ::  name
    ! number of levels
    integer                             ::  nlev
    ! hybride half level coeff; indices 0, 1, ..., nlev
    real, pointer                       ::  a(:), b(:)
    ! hybride full level coeff; indices 1, ..., nlev
    real, pointer                       ::  fa(:), fb(:)
    ! half level coeff per full level (1:2,1:nlev)
    real, pointer                       ::  fa_bounds(:,:), fb_bounds(:,:)
    ! standard pressures and mass (surface pressure 1e5)
    real, pointer                       ::  p0(:), fp0(:), m0(:)
    ! upwards od downwards ?
    character(len=1)                    ::  updo
    ! hybride increments; indices 1, ..., nlev
    real, pointer                       ::  da(:), db(:)
    ! -- parent levels
    logical                             ::  selection
    integer, pointer                    ::  hlevs(:)
    integer, pointer                    ::  flevs(:,:)
    integer                             ::  nlev_parent
    real, pointer                       ::  a_parent(:), b_parent(:)
  end type TLevelInfo
  
  
  ! --- const ------------------------------------------------------------------
  
  ! *** 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                 /)   

  ! *** 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 /)


  ! *** 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
  
  
contains



  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:  PrintInfo
  !
  ! !DESCRIPTION: Basic print method for hybrid level objects.
  !\\
  !\\
  ! !INTERFACE:
  !
  subroutine PrintInfo( LevelInfo )
    !
    ! !USES:
    !
    use GO,                 only : gol, goPr
    !
    ! !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
    
    !    ! hybride increments; indices 1, ..., nlev
    !    real, pointer                       ::  da(:), db(:)
    !    ! -- parent levels
    !    logical                             ::  selection
    !    integer, pointer                    ::  hlevs(:)
    !    integer, pointer                    ::  flevs(:,:)
    !    integer                             ::  nlev_parent
    !    real, pointer                       ::  a_parent(:), b_parent(:)
    
    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
      write (*,'("ERROR in ",a)') rname; 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.
    nullify( levi%hlevs )
    nullify( levi%flevs )
    nullify( levi%a_parent )
    nullify( levi%b_parent )
    
    ! ok
    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
    
    ! --- const --------------------------------------
    
    character(len=*), parameter  ::  rname = mname//'/levi_Init_levi'
    
    ! --- local --------------------------------------
    
    integer             ::  l
    
    ! --- begin ---------------------------------------
    
    ! copy fields:
    call Init( levi, levi2%nlev, levi2%a, levi2%b, status, name=trim(levi2%name) )
    if (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; return; end if

    ! ok
    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 ( 'ec19' )
        call Init( levi, 19, a_ec19, b_ec19, status, name=trim(key) )
        if (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; return; end if

      case ( 'ec31' )
        call Init( levi, 31, a_ec31, b_ec31, status, name=trim(key) )
        if (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; return; end if

      case ( 'ec40' )
        call Init( levi, 40, a_ec40, b_ec40, status, name=trim(key) )
        if (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; return; end if

      case ( 'ec60' )
        call Init( levi, 60, a_ec60, b_ec60, status, name=trim(key) )
        if (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; return; end if

      case ( 'ec62' )
        call Init( levi, 62, a_ec62, b_ec62, status, name=trim(key) )
        if (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; return; end if

      case ( 'ec91' )
        call Init( levi, 91, a_ec91, b_ec91, status, name=trim(key) )
        if (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; return; end if

      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 (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; return; end if
        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 (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; return; end if
        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 (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; return; end if
        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 (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; return; end if
        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 (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; return; end if
        deallocate( a )
        deallocate( b )

      case ( 'nc28' )
        call Init( levi, 28, a_nc28, b_nc28, status, name=trim(key) )
        if (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; return; end if

      case ( 'msc71' )
        call Init( levi, 71, a_msc71, b_msc71, status, name=trim(key) )
        if (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; return; end if

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

    ! ok
    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 ---------------------------------------
    
    ! check length of array:
    if ( (size(hlev_select) < 2) .or. (size(hlev_select) > levi_parent%nlev+1) ) then
      write (*,'("ERROR - strange length of array with selected levels:")')
      write (*,'("ERROR -   expected : ",i4,", ..,",i4)') 2, levi_parent%nlev+1
      write (*,'("ERROR -   found    : ",i4)') size(hlev_select)
      write (*,'("ERROR in ",a)') rname; status=1; return
    end if
    
    ! check range of values:
    if ( (minval(hlev_select) /= 0) .or. (maxval(hlev_select) /= levi_parent%nlev) ) then
      write (*,'("ERROR - invalid range of selected levels:")')
      write (*,'("ERROR -   expected : ",i4,", ..,",i4)') 0, levi_parent%nlev
      write (*,'("ERROR -   found    : ",i4,", ..,",i4)') minval(hlev_select), maxval(hlev_select)
      write (*,'("ERROR in ",a)') rname; 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 (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; return; end if

    ! store selection info:    
    levi%selection = .true.
    ! structure will contain selection info:
    !  o levi_parent half levels 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
    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) )
    levi%a_parent = levi_parent%a
    levi%b_parent = levi_parent%b
    
    ! ok
    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 )
    end if
    
    ! ok
    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
      write (*,'("ERROR in ",a)') rname; 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.")')
      write (*,'("ERROR in ",a)') rname; 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)
      write (*,'("ERROR in ",a)') rname; 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.")')
      write (*,'("ERROR in ",a)') rname; 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)
      write (*,'("ERROR in ",a)') rname; status=1; return
    end if
    
    ! ok
    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 (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; return; end if
    
    if ( size(ll,3) /= levi%nlev ) then
      write (*,'("ERROR - size of 3D grid does not match with level info")')
      write (*,'("ERROR -   size(ll,3) : ",i4)') size(ll,3)
      write (*,'("ERROR -   levi%nlev  : ",i4)') levi%nlev
      write (*,'("ERROR in ",a)') rname; status=1; return
    end if
    
    ! ok
    status = 0
    
  end subroutine levi_Check_3d


  
  ! ##############################################################
  ! ###
  ! ### combining levels
  ! ###
  ! ##############################################################
  
  !
  ! How could levi be formed from leviX ?
  !
  !  o levi could be combined out of leviX levels (thus leviX is a superset)
  !  o levels (superset?) in leviX might be a reversed version of levi levels
  !


  subroutine levi_Compare( levi, leviX, verbose, fillkey, reverse, status )
  
    ! --- in/out ----------------------------------

    type(TLevelInfo), intent(in)         ::  levi
    type(TLevelInfo), intent(in)         ::  leviX
    integer, intent(in)                  ::  verbose
    character(len=10), intent(out)       ::  fillkey
    logical, intent(out)                 ::  reverse
    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
    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

    ! by 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
      
    ! selection ?
    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
      
    end if
    
    ! info on matching ?
    !if ( verbose > 0 ) then
    !  write (*,'("DEBUG - try to match levels:")')
    !  write (*,'("DEBUG -   requested :")')
    !  write (*,'("DEBUG -        l        a              b        ")')
    !  write (*,'("DEBUG -     ----  --------------  --------------")')
    !  do l = 0, levi%nlev
    !    write (*,'("DEBUG -     ",i4,2f16.8)') l, levi%a(l), levi%b(l)
    !  end do
    !  if ( levi%selection ) then
    !    write (*,'("DEBUG -   parent :")')
    !    write (*,'("DEBUG -        l        a              b        ")')
    !    write (*,'("DEBUG -     ----  --------------  --------------")')
    !    do l = 0, levi%nlev_parent
    !      write (*,'("DEBUG -     ",i4,2f16.8)') l, levi%a_parent(l), levi%b_parent(l)
    !    end do
    !  end if
    !  write (*,'("DEBUG -   found :")')
    !  write (*,'("DEBUG -        l        a              b        ")')
    !  write (*,'("DEBUG -     ----  --------------  --------------")')
    !  do l = 0, leviX%nlev
    !    write (*,'("DEBUG -     ",i4,2f16.8)') l, leviX%a(l), leviX%b(l)
    !  end do
    !  write (*,'("DEBUG -   fillkey : ",a)') trim(fillkey)
    !  write (*,'("DEBUG -   reverse : ",l2)') reverse
    !  write (*,'("DEBUG in ",a)') rname; status = 1 ; return
    !end if

    ! done    
    deallocate( Xa_rev, Xb_rev )
    
    ! ok
    status = 0
    
  end subroutine levi_Compare
  
  
  ! =====================================================================
  ! ===
  ! === pressure levels
  ! ===
  ! =====================================================================
  
  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)
      write (*,'("ERROR in ",a)') rname; status=1; return
    end if
    
    ! half level pressure
    ph = levi%a + levi%b * sp
    
    ! ok
    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 ---------------------------------
    
    ! check ...
    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)
      write (*,'("ERROR in ",a)') rname; status=1; return
    end if
    
    ! half level pressure
    do i = 1, size(sp)
      ph(i,:) = levi%a + levi%b * sp(i)
    end do
    
    ! ok
    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 ---------------------------------
    
    ! check ...
    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)
      write (*,'("ERROR in ",a)') rname; 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
    
    ! ok
    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 ---------------------------------
    
    ! check ...
    if ( size(pf) /= levi%nlev ) then
      write (*,'("ERROR - shape of output grid does not match definition:")')
      write (*,'("ERROR -   full levels  : ",i4)') levi%nlev
      write (*,'("ERROR -   pf shape     : ",i4)') shape(pf)
      write (*,'("ERROR in ",a)') rname; status=1; return
    end if
    
    ! half level pressure
    pf = levi%fa + levi%fb * sp
    
    ! ok
    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 ---------------------------------
    
    ! check ...
    if ( any( shape(pf) /= (/shape(sp),levi%nlev/) ) ) then
      write (*,'("ERROR - shape of output grid does not match definition:")')
      write (*,'("ERROR -   sp shape     : ",i6)') shape(sp)
      write (*,'("ERROR -   full levels  : ",i4)') levi%nlev
      write (*,'("ERROR -   pf shape     : ",i6,i4)') shape(pf)
      write (*,'("ERROR in ",a)') rname; status=1; return
    end if
    
    ! half level pressure
    do i = 1, size(sp)
      pf(i,:) = levi%fa + levi%fb * sp(i)
    end do
    
    ! ok
    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 ---------------------------------
    
    ! check ...
    if ( any( shape(pf) /= (/shape(sp),levi%nlev/) ) ) then
      write (*,'("ERROR - shape of output grid does not match definition:")')
      write (*,'("ERROR -   sp shape     : ",2i6)') shape(sp)
      write (*,'("ERROR -   full levels  : ",i4)') levi%nlev
      write (*,'("ERROR -   pf shape     : ",2i6,i4)') shape(pf)
      write (*,'("ERROR in ",a)') rname; status=1; return
    end if
    
    ! half 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
    
    ! ok
    status = 0
    
  end subroutine levi_FPressure_3d

  
   
  ! =====================================================================
  ! ===
  ! === 3D conversion
  ! ===
  ! =====================================================================
   
  
  !
  ! Copy llX (defined by leviX) into ll (defined by levi),
  ! eventually by combining levels.
  ! How levels are combined is specified by a key:
  !  'bottom'     :  use bottom value (most close to the ground)
  !  'top'        :  use top value (most close to the model top)
  !  'sum'        :  sum values
  !  'aver'       :  average of all levels
  !  'mass-aver'  :  mass weighted average
  ! Eventually reverse.
  ! Surface pressure field is used to compute massed in case of 'mass-aver' combination.
  !

  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 -------------------------------
    
    ! for safety ...
    status = 1
    
    ! correct number of levels ?
    if ( size(ll) /= levi%nlev ) then
      write (*,'("ERROR - number of levels in output grid does not match definition:")')
      write (*,'("ERROR -   ll  levels  : ",i3)') size(ll)
      write (*,'("ERROR -   levi nlev   : ",i3)') levi%nlev
      write (*,'("ERROR in ",a)') rname; status=1; return
    end if
    
    ! necessary to combine or reverse ?
    verbose = 1
    call Compare( levi, leviX, verbose, fillkey, reverse, status )
    if (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; return; end if

    ! fill given key:
    select case ( fillkey )

      !
      ! ** copy levels
      !
          
      case ( 'copy' )

        ! loop over all TM levels
        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

      !
      ! ** combine levels
      !
          
      case ( 'combine' )
      
        ! loop over all TM levels
        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
                write (*,'("ERROR - combine key `",a,"` not supported.")') trim(combine_key)
                write (*,'("ERROR in ",a)') rname; status=1; return
            end select
          end do
        end do

    end select
    
    ! ok
    status = 0

  end subroutine levi_FillLevels_1d
  
  
  ! ***
  

!  subroutine levi_FillLevels_2d( 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_2d'
!    
!    ! --- local ------------------------------
!    
!    integer            ::  verbose
!    character(len=10)  ::  fillkey
!    logical            ::  reverse
!
!    integer            ::  k, l
!    integer            ::  le, nle, le1, le2
!    real               ::  Xfp0, Xfp0_save
!    
!    ! --- begin -------------------------------
!    
!    ! for safety ...
!    status = 1
!    
!    ! correct number of levels ?
!    if ( size(ll,2) /= levi%nlev ) then
!      write (*,'("ERROR - number of levels in output grid does not match definition:")')
!      write (*,'("ERROR -   ll  levels  : ",i3)') size(ll,2)
!      write (*,'("ERROR -   levi nlev   : ",i3)') levi%nlev
!      write (*,'("ERROR in ",a)') rname; status=1; return
!    end if
!    
!    ! same grid sizes of ll and llX ?
!    if ( (size(ll,1) /= size(llX,1)) ) then
!      write (*,'("ERROR - horizontal sizes do not match:")')
!      write (*,'("ERROR -   ll     : ",i3)') size(ll ,1)
!      write (*,'("ERROR -   llX    : ",i3)') size(llX,1)
!      write (*,'("ERROR in ",a)') rname; status=1; return
!    end if
!
!    ! correct size of surface pressure field ?
!    if ( combine_key == 'mass-aver' ) then
!      if ( (size(ps,1) /= size(ll,1)) ) then
!        write (*,'("ERROR - horizontal sizes do not match:")')
!        write (*,'("ERROR -   ps     : ",i3)') size(ps ,1)
!        write (*,'("ERROR -   ll     : ",i3)') size(ll ,1)
!        write (*,'("ERROR in ",a)') rname; status=1; return
!      end if
!    end if
!
!    ! necessary to combine or reverse ?
!    verbose = 1
!    call Compare( levi, leviX, verbose, fillkey, reverse, status )
!    if (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; return; end if
!
!    ! fill given key:
!    select case ( fillkey )
!
!      !
!      ! ** copy levels
!      !
!          
!      case ( 'copy' )
!
!        ! loop over all TM levels
!        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
!
!      !
!      ! ** combine levels
!      !
!          
!      case ( 'combine' )
!      
!        ! loop over all TM levels
!        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
!                write (*,'("ERROR - combine key `",a,"` not supported.")') trim(combine_key)
!                write (*,'("ERROR in ",a)') rname; status=1; return
!            end select
!          end do
!        end do
!
!      !
!      ! ** error ..
!      !
!          
!      case default
!        write (*,'("ERROR - fill key `",a,"` not supported.")') trim(fillkey)
!        write (*,'("ERROR in ",a)') rname; status=1; return
!    end select
!    
!    ! ok
!    status = 0
!
!  end subroutine levi_FillLevels_2d
  

  ! ***
  

  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 -------------------------------
  
    ! for safety ...
    status = 1
    
    ! 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)
      write (*,'("ERROR in ",a)') rname; 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)
      write (*,'("ERROR in ",a)') rname; 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)
      write (*,'("ERROR in ",a)') rname; status=1; return
    end if

    ! necessary to combine or reverse ?
    verbose = 0
    call Compare( levi, leviX, verbose, fillkey, reverse, status )
    if (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; return; end if

    ! 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
        write (*,'("ERROR in ",a)') rname; status=1; return
      end if
    end if
    
    ! levels or half-levels ?
    select case ( nw )
    
      !
      ! === levels ===================================
      !
      
      case ( 'n' )

        ! fill given key:
        select case ( fillkey )

          !
          ! ** copy levels
          !

          case ( 'copy' )

            ! loop over all TM levels
            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

          !
          ! ** combine levels
          !

          case ( 'combine' )

            ! loop over all TM levels
            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
                    write (*,'("ERROR - combine key not supported:")')
                    write (*,'("ERROR -   combine key : ",a)') trim(combine_key)
                    write (*,'("ERROR -   fill key    : ",a)') trim(fillkey)
                    write (*,'("ERROR -   nw key      : ",a)') trim(nw)
                    write (*,'("ERROR in ",a)') rname; status=1; return
                end select
              end do
            end do

          !
          ! ** interpol
          !

          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 : ",es10.4)') leviX%a
                      write (*,'("ERROR -   leviX%b : ",es10.4)') leviX%b
                      write (*,'("ERROR in ",a)') rname; 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 (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; return; end if
                    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..

                  ! loop over all levels:
                  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 (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; return; end if
                  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..

                  ! loop over all levels:
                  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 (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; return; end if
                  end do

                end do

              case default
                write (*,'("ERROR - combine key not supported:")')
                write (*,'("ERROR -   combine key : ",a)') trim(combine_key)
                write (*,'("ERROR -   fill key    : ",a)') trim(fillkey)
                write (*,'("ERROR -   nw key      : ",a)') trim(nw)
                write (*,'("ERROR in ",a)') rname; status=1; return
            end select

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

          !
          ! ** error ..
          !

          case default
            write (*,'("ERROR - fill key not supported:")')
            write (*,'("ERROR -   fill key    : ",a)') trim(fillkey)
            write (*,'("ERROR -   nw key      : ",a)') trim(nw)
            write (*,'("ERROR in ",a)') rname; status=1; return
        end select
        
      !
      ! === half levels ========================================
      !

      case ( 'w' )

        ! fill given key:
        select case ( fillkey )

          !
          ! ** copy levels
          !

          case ( 'copy' )

            ! loop over all TM levels
            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

          !
          ! ** combine levels
          !

          case ( 'combine' )
          
            ll = 0.0

            ! loop over target half levels
            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

          !
          ! ** interpol
          !

          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

                  ! loop over all half levels:
                  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 (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; return; end if
                  end do

                end do

              case default
                write (*,'("ERROR - combine key not supported:")')
                write (*,'("ERROR -   combine key : ",a)') trim(combine_key)
                write (*,'("ERROR -   fill key    : ",a)') trim(fillkey)
                write (*,'("ERROR -   nw key      : ",a)') trim(nw)
                write (*,'("ERROR in ",a)') rname; status=1; return
            end select

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

          !
          ! ** error ..
          !

          case default
            write (*,'("ERROR - fill key not supported:")')
            write (*,'("ERROR -   fill key    : ",a)') trim(fillkey)
            write (*,'("ERROR -   nw key      : ",a)') trim(nw)
            write (*,'("ERROR in ",a)') rname; status=1; return
        end select
        
      !
      ! === error ========================================
      !

      case default

        write (*,'("ERROR - nw key `",a,"` not supported.")') trim(nw)
        write (*,'("ERROR in ",a)') rname; status=1; return

    end select

    ! ok
    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 -------------------------------
  
    ! for safety ...
    status = 1
    
    ! correct number of levels ?
    if ( ((nw == 'n') .and. (size(ll,3) /= levi%nlev  )) .or. &
         ((nw == 'w') .and. (size(ll,3) /= 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,3)
      write (*,'("ERROR in ",a)') rname; status=1; return
    end if
    
    ! correct number of levels ?
    if ( ((nw == 'n') .and. (size(llX,3) /= leviX%nlev  )) .or. &
         ((nw == 'w') .and. (size(llX,3) /= 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,3)
      write (*,'("ERROR in ",a)') rname; 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 (*,'("ERROR - horizontal sizes do not match:")')
      write (*,'("ERROR -   ll     : ",i3," x ",i3)') size(ll,1), size(ll,2)
      write (*,'("ERROR -   llX    : ",i3," x ",i3)') size(llX,1), size(llX,2)
      write (*,'("ERROR in ",a)') rname; status=1; return
    end if

    ! necessary to combine or reverse ?
    verbose = 0
    call Compare( levi, leviX, verbose, fillkey, reverse, status )
    if (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; return; end if

    ! 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 (*,'("ERROR - horizontal sizes do not match:")')
        write (*,'("ERROR -   ps          : ",i3," x ",i3)') size(ps ,1), size(ps ,2)
        write (*,'("ERROR -   ll          : ",i3," x ",i3)') size(ll ,1), size(ll ,2)
        write (*,'("ERROR -   combine_key : ",a)') combine_key
        write (*,'("ERROR -   fillkey     : ",a)') fillkey
        write (*,'("ERROR -   reverse     : ",l1)') reverse
        write (*,'("ERROR in ",a)') rname; status=1; return
      end if
    end if
    
    ! levels or half-levels ?
    select case ( nw )
    
      !
      ! === levels ===================================
      !
      
      case ( 'n' )

        ! fill given key:
        select case ( fillkey )

          !
          ! ** copy levels
          !

          case ( 'copy' )

            ! loop over all TM levels
            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

          !
          ! ** combine levels
          !

          case ( 'combine' )

            ! loop over all TM levels
            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')
                    !PLS : assumes that height-ave is akind to so-called 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
                    write (*,'("ERROR - combine key not supported:")')
                    write (*,'("ERROR -   combine key : ",a)') trim(combine_key)
                    write (*,'("ERROR -   fill key    : ",a)') trim(fillkey)
                    write (*,'("ERROR -   nw key      : ",a)') trim(nw)
                    write (*,'("ERROR in ",a)') rname; status=1; return
                end select
              end do
            end do

          !
          ! ** interpol
          !

          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 (*,'("ERROR - from IntervalSum (combine key `sum`)")')
                        write (*,'("ERROR -   i, j    : ",2i4   )') i, j
                        write (*,'("ERROR -   ps      : ", f16.4)') ps(i,j)
                        write (*,'("ERROR -   leviX%a : ",es10.4)') leviX%a
                        write (*,'("ERROR -   leviX%b : ",es10.4)') leviX%b
                        write (*,'("ERROR in ",a)') rname; 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' )

                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 (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; return; end if
                      ll(i,j,l) = ll(i,j,l) / (phalf(l)-phalf(l-1))*pfac
                    end do

                  end do
                end do

             case ( 'height-ave' )

                !PLS : assumes that height-ave is akind to so-called mass-aver.
                
                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 (status/=0) then;
                            write (*,'("ERROR in ",a)') rname; status=1; return;
                         end if
                         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..

                    ! loop over all levels:
                    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 (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; return; end if
                    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..

                    ! loop over all levels:
                    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 (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; return; end if
                    end do

                  end do
                end do

              case default
                write (*,'("ERROR - combine key not supported:")')
                write (*,'("ERROR -   combine key : ",a)') trim(combine_key)
                write (*,'("ERROR -   fill key    : ",a)') trim(fillkey)
                write (*,'("ERROR -   nw key      : ",a)') trim(nw)
                write (*,'("ERROR in ",a)') rname; status=1; return
            end select

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

          !
          ! ** error ..
          !

          case default
            write (*,'("ERROR - fill key not supported:")')
            write (*,'("ERROR -   fill key    : ",a)') trim(fillkey)
            write (*,'("ERROR -   nw key      : ",a)') trim(nw)
            write (*,'("ERROR in ",a)') rname; status=1; return
        end select
        
      !
      ! === half levels ========================================
      !

      case ( 'w' )

        ! fill given key:
        select case ( fillkey )

          !
          ! ** copy levels
          !

          case ( 'copy' )

            ! loop over all TM levels
            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

          !
          ! ** combine levels
          !

          case ( 'combine' )
          
            ll = 0.0

            ! loop over target half levels
            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

          !
          ! ** interpol
          !

          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

                    ! loop over all half levels:
                    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 (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; return; end if
                    end do

                  end do
                end do

              case default
                write (*,'("ERROR - combine key not supported:")')
                write (*,'("ERROR -   combine key : ",a)') trim(combine_key)
                write (*,'("ERROR -   fill key    : ",a)') trim(fillkey)
                write (*,'("ERROR -   nw key      : ",a)') trim(nw)
                write (*,'("ERROR in ",a)') rname; status=1; return
            end select

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

          !
          ! ** error ..
          !

          case default
            write (*,'("ERROR - fill key not supported:")')
            write (*,'("ERROR -   fill key    : ",a)') trim(fillkey)
            write (*,'("ERROR -   nw key      : ",a)') trim(nw)
            write (*,'("ERROR in ",a)') rname; status=1; return
        end select
        
      !
      ! === error ========================================
      !

      case default

        write (*,'("ERROR - nw key `",a,"` not supported.")') trim(nw)
        write (*,'("ERROR in ",a)') rname; status=1; return

    end select

    ! ok
    status = 0

  end subroutine levi_FillLevels_3d
  
  
  
end module grid_type_hyb
