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