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