1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677787980818283848586878889909192939495969798991001011021031041051061071081091101111121131141151161171181191201211221231241251261271281291301311321331341351361371381391401411421431441451461471481491501511521531541551561571581591601611621631641651661671681691701711721731741751761771781791801811821831841851861871881891901911921931941951961971981992002012022032042052062072082092102112122132142152162172182192202212222232242252262272282292302312322332342352362372382392402412422432442452462472482492502512522532542552562572582592602612622632642652662672682692702712722732742752762772782792802812822832842852862872882892902912922932942952962972982993003013023033043053063073083093103113123133143153163173183193203213223233243253263273283293303313323333343353363373383393403413423433443453463473483493503513523533543553563573583593603613623633643653663673683693703713723733743753763773783793803813823833843853863873883893903913923933943953963973983994004014024034044054064074084094104114124134144154164174184194204214224234244254264274284294304314324334344354364374384394404414424434444454464474484494504514524534544554564574584594604614624634644654664674684694704714724734744754764774784794804814824834844854864874884894904914924934944954964974984995005015025035045055065075085095105115125135145155165175185195205215225235245255265275285295305315325335345355365375385395405415425435445455465475485495505515525535545555565575585595605615625635645655665675685695705715725735745755765775785795805815825835845855865875885895905915925935945955965975985996006016026036046056066076086096106116126136146156166176186196206216226236246256266276286296306316326336346356366376386396406416426436446456466476486496506516526536546556566576586596606616626636646656666676686696706716726736746756766776786796806816826836846856866876886896906916926936946956966976986997007017027037047057067077087097107117127137147157167177187197207217227237247257267277287297307317327337347357367377387397407417427437447457467477487497507517527537547557567577587597607617627637647657667677687697707717727737747757767777787797807817827837847857867877887897907917927937947957967977987998008018028038048058068078088098108118128138148158168178188198208218228238248258268278288298308318328338348358368378388398408418428438448458468478488498508518528538548558568578588598608618628638648658668678688698708718728738748758768778788798808818828838848858868878888898908918928938948958968978988999009019029039049059069079089099109119129139149159169179189199209219229239249259269279289299309319329339349359369379389399409419429439449459469479489499509519529539549559569579589599609619629639649659669679689699709719729739749759769779789799809819829839849859869879889899909919929939949959969979989991000100110021003100410051006100710081009101010111012101310141015101610171018101910201021102210231024102510261027102810291030103110321033103410351036103710381039104010411042104310441045104610471048104910501051105210531054105510561057105810591060106110621063106410651066106710681069107010711072107310741075107610771078107910801081108210831084108510861087108810891090109110921093109410951096109710981099110011011102110311041105110611071108110911101111111211131114111511161117111811191120112111221123112411251126112711281129113011311132113311341135113611371138113911401141114211431144114511461147114811491150115111521153115411551156115711581159116011611162116311641165116611671168116911701171117211731174117511761177117811791180118111821183118411851186118711881189119011911192119311941195119611971198119912001201120212031204120512061207120812091210121112121213121412151216121712181219122012211222122312241225122612271228122912301231123212331234123512361237123812391240124112421243124412451246124712481249125012511252125312541255125612571258125912601261126212631264126512661267126812691270127112721273127412751276127712781279128012811282128312841285128612871288128912901291129212931294129512961297129812991300130113021303130413051306130713081309131013111312131313141315131613171318131913201321132213231324132513261327132813291330133113321333133413351336133713381339134013411342134313441345134613471348134913501351135213531354135513561357135813591360136113621363136413651366136713681369137013711372137313741375137613771378137913801381138213831384138513861387138813891390139113921393139413951396139713981399140014011402140314041405140614071408140914101411141214131414141514161417141814191420142114221423142414251426142714281429143014311432143314341435143614371438143914401441144214431444144514461447144814491450145114521453145414551456145714581459146014611462146314641465146614671468146914701471147214731474147514761477147814791480148114821483148414851486148714881489149014911492149314941495149614971498149915001501150215031504150515061507150815091510151115121513151415151516151715181519152015211522152315241525152615271528152915301531153215331534153515361537153815391540154115421543154415451546154715481549155015511552155315541555155615571558155915601561156215631564156515661567156815691570157115721573157415751576157715781579158015811582158315841585158615871588158915901591159215931594159515961597159815991600160116021603160416051606160716081609161016111612161316141615161616171618161916201621162216231624162516261627162816291630163116321633163416351636163716381639164016411642164316441645164616471648164916501651165216531654165516561657165816591660166116621663166416651666166716681669167016711672167316741675167616771678167916801681168216831684168516861687168816891690169116921693169416951696169716981699170017011702170317041705170617071708170917101711171217131714171517161717171817191720172117221723172417251726172717281729173017311732173317341735173617371738173917401741174217431744174517461747174817491750175117521753175417551756175717581759176017611762176317641765176617671768176917701771177217731774177517761777177817791780178117821783178417851786178717881789179017911792179317941795179617971798179918001801180218031804180518061807180818091810181118121813181418151816181718181819182018211822182318241825182618271828182918301831183218331834183518361837183818391840184118421843184418451846184718481849185018511852185318541855185618571858185918601861186218631864186518661867186818691870187118721873187418751876187718781879188018811882188318841885188618871888188918901891189218931894189518961897189818991900190119021903190419051906190719081909191019111912191319141915191619171918191919201921192219231924192519261927192819291930193119321933193419351936193719381939194019411942194319441945194619471948194919501951195219531954195519561957195819591960196119621963196419651966196719681969197019711972197319741975197619771978197919801981198219831984198519861987198819891990199119921993199419951996199719981999200020012002200320042005200620072008200920102011201220132014201520162017201820192020202120222023202420252026202720282029203020312032203320342035203620372038203920402041204220432044204520462047204820492050205120522053205420552056205720582059206020612062206320642065206620672068206920702071207220732074207520762077207820792080208120822083208420852086208720882089209020912092209320942095209620972098209921002101210221032104210521062107210821092110211121122113211421152116211721182119212021212122212321242125212621272128212921302131213221332134213521362137213821392140214121422143214421452146214721482149215021512152215321542155215621572158215921602161216221632164216521662167216821692170217121722173217421752176217721782179218021812182218321842185218621872188218921902191219221932194219521962197219821992200220122022203220422052206220722082209221022112212221322142215221622172218221922202221222222232224222522262227222822292230223122322233223422352236223722382239224022412242224322442245224622472248224922502251225222532254225522562257225822592260226122622263226422652266226722682269227022712272227322742275227622772278227922802281228222832284228522862287228822892290229122922293229422952296229722982299230023012302230323042305230623072308230923102311231223132314231523162317231823192320232123222323232423252326232723282329233023312332233323342335233623372338233923402341234223432344234523462347234823492350235123522353235423552356235723582359236023612362236323642365236623672368236923702371237223732374237523762377237823792380238123822383238423852386238723882389239023912392239323942395239623972398239924002401240224032404240524062407240824092410241124122413241424152416241724182419242024212422242324242425242624272428242924302431243224332434243524362437243824392440244124422443244424452446244724482449245024512452245324542455245624572458245924602461246224632464246524662467246824692470247124722473247424752476247724782479248024812482248324842485248624872488248924902491249224932494249524962497249824992500250125022503250425052506250725082509251025112512251325142515251625172518251925202521252225232524252525262527252825292530253125322533253425352536253725382539254025412542254325442545254625472548254925502551255225532554255525562557255825592560256125622563256425652566256725682569257025712572257325742575257625772578257925802581258225832584258525862587258825892590259125922593259425952596259725982599260026012602260326042605260626072608260926102611261226132614261526162617261826192620262126222623262426252626262726282629263026312632263326342635263626372638263926402641264226432644264526462647264826492650265126522653265426552656265726582659266026612662266326642665266626672668266926702671267226732674267526762677267826792680268126822683268426852686268726882689269026912692269326942695269626972698269927002701270227032704270527062707270827092710271127122713271427152716271727182719272027212722272327242725272627272728272927302731273227332734273527362737273827392740274127422743274427452746274727482749275027512752275327542755275627572758275927602761276227632764276527662767276827692770277127722773277427752776277727782779278027812782278327842785278627872788278927902791279227932794279527962797279827992800280128022803280428052806280728082809281028112812281328142815281628172818281928202821282228232824282528262827282828292830283128322833283428352836283728382839284028412842284328442845284628472848284928502851285228532854285528562857285828592860286128622863286428652866286728682869287028712872287328742875287628772878287928802881288228832884288528862887288828892890289128922893289428952896289728982899290029012902290329042905290629072908290929102911291229132914291529162917291829192920292129222923292429252926292729282929293029312932293329342935293629372938293929402941294229432944294529462947294829492950295129522953295429552956295729582959296029612962296329642965296629672968296929702971297229732974297529762977297829792980298129822983298429852986298729882989299029912992299329942995299629972998299930003001300230033004300530063007300830093010301130123013301430153016301730183019302030213022302330243025302630273028302930303031303230333034303530363037303830393040304130423043304430453046304730483049305030513052305330543055305630573058305930603061306230633064306530663067306830693070307130723073307430753076307730783079308030813082308330843085308630873088308930903091309230933094309530963097309830993100310131023103310431053106310731083109311031113112311331143115311631173118311931203121312231233124312531263127312831293130313131323133313431353136313731383139314031413142314331443145314631473148314931503151315231533154315531563157315831593160316131623163316431653166316731683169317031713172317331743175317631773178317931803181318231833184318531863187318831893190319131923193319431953196319731983199320032013202320332043205320632073208320932103211321232133214321532163217321832193220322132223223322432253226322732283229323032313232323332343235323632373238323932403241324232433244324532463247324832493250325132523253325432553256325732583259326032613262326332643265326632673268326932703271327232733274327532763277327832793280328132823283328432853286328732883289329032913292329332943295329632973298329933003301330233033304330533063307330833093310331133123313331433153316331733183319332033213322332333243325332633273328332933303331333233333334333533363337333833393340334133423343334433453346334733483349335033513352335333543355335633573358335933603361336233633364336533663367336833693370337133723373337433753376337733783379338033813382338333843385338633873388338933903391339233933394339533963397 |
- !
- ! may 2002, Arjo Segers
- !
- ! 23 Oct 2013 - Ph. Le Sager - added '==' operator to compare grids,
- ! and '=' assignment to copy data instead of pointers,
- ! and initialization of TllGridInfo pointers (requires F95)
- !
- ! NOTES: (0) "lli" is commonly used in the code, and stands for LonLatInfo
- ! (1) if you declare an object of TllGridInfo type in a module, it must
- ! have the 'save' attribute (unless it is a pointer or allocatable array)
- ! (2) you must call DONE on all lli that have been initialized
- ! (either through a call to INIT or as a copy of another lli) to
- ! avoid memory leak.
- ! (3) because of the associated status of the pointers of a
- ! TllGridInfo type, an argument object of that type should never
- ! have the intent(out) attribute, but the intent(inout) instead
- !
- ! STAND ALONE VERSION FOR LAT/LON GRIDS
- !
- ! Apply factor over region.
- !
- ! use grid, only : TRegion
- !
- ! type(TllRegion) :: llreg
- !
- ! ! define region by either calling INIT or COPY (i.e. "=")
- ! call Init( llreg, west_deg, east_deg, south_deg, north_deg, status )
- !
- ! ! apply factor to data set within box, or to the complement;
- ! ! if a grid cell only partly covers the region, the factor
- ! ! is applied according to the area ratio;
- ! ! x is at least 2d, if it is 3d the factor is applied for all levels:
- ! call Region_Apply_Factor( lli, x, llreg, fac, status [,complement=.false.] )
- !
- ! ! done
- ! call Done( llreg, status )
- !
- module grid_type_ll
-
- implicit none
-
- ! --- in/out --------------------------------
-
- private
-
- public :: TllGridInfo
-
- public :: Init, Done
- public :: Check
- public :: operator(==), assignment(=)
-
- public :: AreaOper
- public :: InterpolFractions
- public :: Interpol
-
- public :: BalanceMassFluxes
- public :: CheckMassBalance
-
- public :: GetRefinement
- public :: Match
-
- public :: FillGrid
-
- public :: EquivLat
-
- ! ~~ region
-
- public :: TllRegion
- public :: Region_Apply_Factor
-
-
- ! --- const ---------------------------------
-
- character(len=*), parameter :: mname = 'grid_type_ll'
-
- ! --- types ---------------------------------
-
- ! *** location, size, etc
-
- type TllGridInfo
- character(len=32) :: name
- ! * radius
- real :: R ! m
- ! * spacing
- real :: dlon_deg, dlat_deg ! degrees
- real :: dlon, dlat ! radians
- ! * size
- integer :: im, nlon
- integer :: jm, nlat
- ! * coordinates of gridpoint (cell center)
- ! indices 1, 2, ...
- real, pointer :: lon_deg(:) => null() ! degrees
- real, pointer :: lat_deg(:) => null() ! degrees
- real, pointer :: lon(:) => null() ! rad
- real, pointer :: lat(:) => null() ! rad
- ! * boundaries in a rank-2 array:
- real, pointer :: lon_bounds_deg(:,:) => null() ! degrees
- real, pointer :: lat_bounds_deg(:,:) => null() ! degrees
- real, pointer :: lon_bounds(:,:) => null() ! rad
- real, pointer :: lat_bounds(:,:) => null() ! rad
- ! * coordinates of boundaries for cell around grid point;
- ! indices 0, 1, 2, ...
- real, pointer :: blon_deg(:) => null() ! degrees
- real, pointer :: blat_deg(:) => null() ! degrees
- real, pointer :: blon(:) => null() ! rad
- real, pointer :: blat(:) => null() ! rad
- ! * area for cell in certain row
- real, pointer :: area(:) => null() ! rad^2
- real, pointer :: area_m2(:) => null() ! m^2
- ! * cell length in m
- real, pointer :: dx(:) => null() ! m
- real, pointer :: bdx(:) => null() ! m
- real :: dy, bdy ! m
- end type TllGridInfo
-
-
- ! ~~ region
-
-
- type TllRegion
- ! region boundaries in degrees:
- real :: west_deg, east_deg, south_deg, north_deg
- ! idem in radians:
- real :: west, east, south, north
- end type TllRegion
- ! --- interfaces ----------------------------
-
- interface Init
- module procedure llgridinfo_Init
- module procedure llreg_Init
- end interface
-
- interface Done
- module procedure llgridinfo_Done
- module procedure llreg_Done
- end interface
-
- interface Check
- module procedure llgrid_Check_i
- module procedure llgrid_Check_r
- end interface
-
- interface AreaOper
- module procedure llgrid_AreaOper_2d
- module procedure llgrid_AreaOper_3d
- end interface
-
- interface InterpolFractions
- module procedure llgrid_InterpolFractions
- end interface
-
- interface Interpol
- module procedure llgrid_Eval_2d
- module procedure llgrid_Eval_3d
- end interface
-
- interface Match
- module procedure llgrid_Match
- end interface
- interface operator(==)
- module procedure llgrid_equal_llgrid
- end interface
- interface assignment(=)
- module procedure copy_llgrid
- end interface
-
-
- interface BalanceMassFluxes
- module procedure BalanceMassFluxes_sm
- module procedure BalanceMassFluxes_m
- end interface
- interface EquivLat
- module procedure llgrid_EquivLat
- module procedure llgrid_EquivLat_sort
- end interface
-
- ! ~~ region
-
- interface Region_Apply_Factor
- module procedure llreg_Region_Apply_Factor_2d
- module procedure llreg_Region_Apply_Factor_3d
- end interface
-
- contains
- ! ========================================================
- ! ===
- ! === Init, Done
- ! ===
- ! ========================================================
- !
- ! blat(j) +-------+
- ! | | |
- ! lat(j) |---o---|
- ! | | |
- ! +-------+
- ! lon(i)
- ! blon(i)
- !
- !--------------------------------------------------------------------------
- ! TM5 !
- !--------------------------------------------------------------------------
- !BOP
- !
- ! !IROUTINE: LLGRIDINFO_INIT
- !
- ! !DESCRIPTION: Initialize Longitude-Latitude Grid object
- !\\
- !\\
- ! !INTERFACE:
- !
- SUBROUTINE LLGridInfo_Init( lli, west_deg, dlon_deg, im, &
- south_deg, dlat_deg, jm, status, name )
- !
- ! !USES:
- !
- use Grid_Tools, only : ll_area
- use grid_tools, only : deg2rad, ae
- !
- ! !OUTPUT PARAMETERS:
- !
- type(TllGridInfo), intent(inout) :: lli
- integer, intent(out) :: status
- !
- ! !INPUT PARAMETERS:
- !
- real, intent(in) :: west_deg ! longitude center of 1st grid cell
- real, intent(in) :: dlon_deg ! spacing in degree
- integer, intent(in) :: im
- real, intent(in) :: south_deg ! latitude center of 1st grid cell
- real, intent(in) :: dlat_deg ! spacing in degree
- integer, intent(in) :: jm
- character(len=*), intent(in), optional :: name
- !
- ! !REVISION HISTORY:
- ! 1 Apr 2011 - P. Le Sager - some doc
- !
- ! !REMARKS: the 1st grid cell is the one at lower left (i.e. the further
- ! west and further south cell)
- !
- !EOP
- !------------------------------------------------------------------------
- !BOC
- integer :: i, j
-
- ! --- begin ---------------------------------
- if ( present(name) ) then
- lli%name = name
- else
- lli%name = 'anonymous'
- end if
- ! *** radius
-
- lli%R = ae
-
- ! *** grid spacing
-
- lli%dlon_deg = dlon_deg
- lli%dlon = dlon_deg*deg2rad
- lli%dlat_deg = dlat_deg
- lli%dlat = dlat_deg*deg2rad
-
- ! *** grid range
- lli%im = im
- lli%nlon = im
- lli%jm = jm
- lli%nlat = jm
- ! *** grid points
-
- ! * coor of grid points
-
- ! east-west
- if ( associated( lli%lon_deg) ) deallocate(lli%lon_deg)
- allocate( lli%lon_deg(im) )
- do i = 1, im
- lli%lon_deg(i) = west_deg + (i-1)*dlon_deg
- end do
- if ( associated( lli%lon) ) deallocate(lli%lon)
- allocate( lli%lon(im) )
- lli%lon = lli%lon_deg * deg2rad
- ! south-north
- if ( associated( lli%lat_deg) ) deallocate(lli%lat_deg)
- allocate( lli%lat_deg(jm) )
- do j = 1, jm
- lli%lat_deg(j) = south_deg + (j-1)*dlat_deg
- end do
- if ( associated( lli%lat) ) deallocate(lli%lat)
- allocate( lli%lat(jm) )
- lli%lat = lli%lat_deg * deg2rad
- ! *** cells with grid point in center;
- ! grid point at pole is at top of triangle !
- ! * bounds
- ! west-east
- if ( associated( lli%blon_deg) ) deallocate(lli%blon_deg)
- allocate( lli%blon_deg(0:im) )
- do i = 0, im
- lli%blon_deg(i) = west_deg + (i-0.5)*dlon_deg
- end do
- if ( associated( lli%blon) ) deallocate(lli%blon)
- allocate( lli%blon(0:im) )
- lli%blon = lli%blon_deg * deg2rad
- ! south-north
- if ( associated( lli%blat_deg) ) deallocate(lli%blat_deg)
- allocate( lli%blat_deg(0:jm) )
- do j = 0, jm
- lli%blat_deg(j) = south_deg + (j-0.5)*dlat_deg
- end do
- if ( lli%blat_deg(0) < -90.0 ) lli%blat_deg(0) = -90.0
- if ( lli%blat_deg(0) > 90.0 ) lli%blat_deg(0) = 90.0
- if ( lli%blat_deg(jm) < -90.0 ) lli%blat_deg(jm) = -90.0
- if ( lli%blat_deg(jm) > 90.0 ) lli%blat_deg(jm) = 90.0
- if ( associated( lli%blat) ) deallocate(lli%blat)
- allocate( lli%blat(0:jm) )
- lli%blat = lli%blat_deg * deg2rad
- ! * bounds in a rank-2 array:
- ! west-east
- if ( associated( lli%lon_bounds_deg) ) deallocate(lli%lon_bounds_deg)
- allocate( lli%lon_bounds_deg(2,im) )
- lli%lon_bounds_deg(1,:) = lli%blon_deg(0:im-1)
- lli%lon_bounds_deg(2,:) = lli%blon_deg(1:im)
- if ( associated( lli%lon_bounds) ) deallocate(lli%lon_bounds)
- allocate( lli%lon_bounds(2,im) )
- lli%lon_bounds = lli%lon_bounds_deg * deg2rad
- ! south-north
- if ( associated( lli%lat_bounds_deg) ) deallocate(lli%lat_bounds_deg)
- allocate( lli%lat_bounds_deg(2,jm) )
- lli%lat_bounds_deg(1,:) = lli%blat_deg(0:jm-1)
- lli%lat_bounds_deg(2,:) = lli%blat_deg(1:jm)
- if ( associated( lli%lat_bounds) ) deallocate(lli%lat_bounds)
- allocate( lli%lat_bounds(2,jm) )
- lli%lat_bounds = lli%lat_bounds_deg * deg2rad
- ! * area of cell in lat band
- ! rad^2 :
- if ( associated( lli%area) ) deallocate(lli%area)
- allocate( lli%area(jm) )
- do j = 1, jm
- lli%area(j) = ll_area( 0.0, lli%dlon, lli%blat(j-1), lli%blat(j) )
- end do
- ! m^2 :
- if ( associated( lli%area_m2) ) deallocate(lli%area_m2)
- allocate( lli%area_m2(jm) )
- lli%area_m2 = lli%area * lli%R**2
- ! * length in m
- ! east-west, mid latitude of cell
- if ( associated( lli%dx) ) deallocate(lli%dx)
- allocate( lli%dx(jm) )
- lli%dx = lli%dlon * lli%R * cos(lli%lat)
- ! east-west, boundaries
- if ( associated( lli%bdx) ) deallocate(lli%bdx)
- allocate( lli%bdx(0:jm) )
- lli%bdx = lli%dlon * lli%R * cos(lli%blat)
- ! north-south, the same for each longitude
- lli%dy = lli%dlat * lli%R
- lli%bdy = lli%dlat * lli%R
- ! ok
- status = 0
- END SUBROUTINE LLGRIDINFO_INIT
- !EOC
- !--------------------------------------------------------------------------
- ! TM5 !
- !--------------------------------------------------------------------------
- !BOP
- !
- ! !IROUTINE: COPY_LLGRID
- !
- ! !DESCRIPTION: make a copy of one LLI
- !\\
- !\\
- ! !INTERFACE:
- !
- SUBROUTINE COPY_LLGRID(tgt_grid, src_grid)
- !
- ! !INPUT/OUTPUT PARAMETERS:
- !
- type(TllGridInfo), intent(inout) :: tgt_grid
- !
- ! !INPUT PARAMETERS:
- !
- type(TllGridInfo), intent(in) :: src_grid
- !
- ! !REVISION HISTORY:
- ! 1 Nov 2013 - Ph. Le Sager -
- !
- ! !REMARKS:
- !
- !EOP
- !------------------------------------------------------------------------
- !BOC
- character(len=*), parameter :: rname = mname//'/copy_llgrid'
- integer :: status
- if (.not.(associated( src_grid%lon_deg ))) then
- write(*,*) " : WARNING : source InfoGrid is not initialized"
- write(*,*) " : no data to copy. Returning."
- return
- endif
- call llgridinfo_Init( tgt_grid, src_grid%lon_deg(1), src_grid%dlon_deg, src_grid%im, &
- src_grid%lat_deg(1), src_grid%dlat_deg, src_grid%jm, status, src_grid%name )
- if (status/=0) write (*,'("ERROR in ",a)') rname
- END SUBROUTINE COPY_LLGRID
- !EOC
-
- !--------------------------------------------------------------------------
- ! TM5 !
- !--------------------------------------------------------------------------
- !BOP
- !
- ! !FUNCTION: LLGRID_EQUAL_LLGRID
- !
- ! !DESCRIPTION: Compare two LLIs. Returns .true. if both have been
- ! initialized and describe the same grid.
- !\\
- !\\
- ! !INTERFACE:
- !
- LOGICAL FUNCTION LLGRID_EQUAL_LLGRID( lliA , lliB )
- !
- ! !INPUT PARAMETERS:
- !
- type(TllGridInfo), intent(in) :: lliA, lliB
- !
- ! !REVISION HISTORY:
- !
- ! 23 Oct 2013 - Ph. Le Sager - v0
- !
- ! !REMARKS:
- ! (1) If both grids are not initialized then .false. is return.
- !
- !EOP
- !------------------------------------------------------------------------
- !BOC
- character(len=*), parameter :: name = mname//'/LLGRID_EQUAL_LLGRID'
- ! -- begin
- ! This would be problematic if using F90 instead of F95
- if ((.not.(associated( lliA%lon_deg ))) .or. (.not.(associated( lliB%lon_deg )))) then
- ! write(gol,*)'WARNING : LL ggrid A not initialized'; call goPr
- llgrid_equal_llgrid = .false.
- return
- end if
-
- ! Just compare the inputs of llgridinfo_Init to ensure equivalence
- llgrid_equal_llgrid = &
- ( lliA%dlon_deg == lliB%dlon_deg ).and. &
- ( lliA%nlon == lliB%nlon ).and. &
- ( lliA%dlat_deg == lliB%dlat_deg ).and. &
- ( lliA%nlat == lliB%nlat ).and. &
- ( lliA%lon_deg(1) == lliB%lon_deg(1) ).and. &
- ( lliA%lat_deg(1) == lliB%lat_deg(1) )
- END FUNCTION LLGRID_EQUAL_LLGRID
- !EOC
-
- ! ===
-
- subroutine llgridinfo_Done( lli, status )
-
- ! --- in/out ---------------------------------
-
- type(TllGridInfo), intent(inout) :: lli
- integer, intent(out) :: status
-
- ! --- begin ---------------------------------
-
- if (associated( lli%lon_deg )) deallocate( lli%lon_deg )
- if (associated( lli%lat_deg )) deallocate( lli%lat_deg )
- if (associated( lli%lon )) deallocate( lli%lon )
- if (associated( lli%lat )) deallocate( lli%lat )
- if (associated( lli%blon_deg )) deallocate( lli%blon_deg )
- if (associated( lli%blat_deg )) deallocate( lli%blat_deg )
- if (associated( lli%blon )) deallocate( lli%blon )
- if (associated( lli%blat )) deallocate( lli%blat )
- if (associated( lli%lon_bounds_deg )) deallocate( lli%lon_bounds_deg )
- if (associated( lli%lat_bounds_deg )) deallocate( lli%lat_bounds_deg )
- if (associated( lli%lon_bounds )) deallocate( lli%lon_bounds )
- if (associated( lli%lat_bounds )) deallocate( lli%lat_bounds )
- if (associated( lli%area )) deallocate( lli%area )
- if (associated( lli%area_m2 )) deallocate( lli%area_m2 )
-
- if (associated( lli%dx )) deallocate( lli%dx )
- if (associated( lli%bdx )) deallocate( lli%bdx )
-
- status = 0
-
- end subroutine llgridinfo_Done
-
- ! =============================================================
-
-
- subroutine llgrid_Check_i( lli, nuv, ll, status )
-
- ! --- in/out ----------------------------------
-
- type(TllGridInfo), intent(in) :: lli
- character(len=*), intent(in) :: nuv
- integer, intent(in) :: ll(:,:)
- integer, intent(out) :: status
-
- ! --- const -----------------------------
-
- character(len=*), parameter :: rname = mname//'/llgrid_Check_i'
-
- ! --- begin ----------------------------------
-
- ! check shape of target grid:
- if ( ((nuv == 'n') .and. ((size(ll,1) /= lli%nlon ) .or. (size(ll,2) /= lli%nlat ))) .or. &
- ((nuv == 'u') .and. ((size(ll,1) /= lli%nlon+1) .or. (size(ll,2) /= lli%nlat ))) .or. &
- ((nuv == 'v') .and. ((size(ll,1) /= lli%nlon ) .or. (size(ll,2) /= lli%nlat+1))) ) then
- write (*,'("ERROR - target array does not match with grid definition:")')
- write (*,'("ERROR - lli : ",i3," x ",i3)') lli%nlon, lli%nlat
- write (*,'("ERROR - nuv : ",a )') nuv
- write (*,'("ERROR - ll : ",i3," x ",i3)') shape(ll)
- write (*,'("ERROR in ",a)') rname; status=1; return
- end if
-
- ! ok
- status = 0
- end subroutine llgrid_Check_i
- ! ***
-
-
- subroutine llgrid_Check_r( lli, nuv, ll, status )
-
- ! --- in/out ----------------------------------
-
- type(TllGridInfo), intent(in) :: lli
- character(len=*), intent(in) :: nuv
- real, intent(in) :: ll(:,:)
- integer, intent(out) :: status
-
- ! --- const -----------------------------
-
- character(len=*), parameter :: rname = mname//'/llgrid_Check_r'
-
- ! --- begin ----------------------------------
-
- ! check shape of target grid:
- if ( ((nuv == 'n') .and. ((size(ll,1) /= lli%nlon ) .or. (size(ll,2) /= lli%nlat ))) .or. &
- ((nuv == 'u') .and. ((size(ll,1) /= lli%nlon+1) .or. (size(ll,2) /= lli%nlat ))) .or. &
- ((nuv == 'v') .and. ((size(ll,1) /= lli%nlon ) .or. (size(ll,2) /= lli%nlat+1))) ) then
- write (*,'("ERROR - target array does not match with grid definition:")')
- write (*,'("ERROR - lli : ",i3," x ",i3)') lli%nlon, lli%nlat
- write (*,'("ERROR - nuv : ",a )') nuv
- write (*,'("ERROR - ll : ",i3," x ",i3)') shape(ll)
- write (*,'("ERROR in ",a)') rname; status=1; return
- end if
-
- ! ok
- status = 0
- end subroutine llgrid_Check_r
- ! =====================================================
-
- ! call AreaOper( lli, ll, '/' | '*' | '=', 'rad2' | 'm2' )
-
- subroutine llgrid_AreaOper_2d( lli, ll, oper, unit, status )
-
- ! --- in/out ----------------------------------
-
- type(TllGridInfo), intent(in) :: lli
- real, intent(inout) :: ll(:,:)
- character(len=*), intent(in) :: unit, oper
- integer, intent(out) :: status
-
- ! --- const -----------------------------
-
- character(len=*), parameter :: rname = mname//'/llgrid_AreaOper_2d'
-
- ! --- local --------------------------------
-
- integer :: j
- real :: cell_area
-
- ! --- begin ----------------------------------
-
- ! check ...
- if ( size(ll,2) /= lli%nlat ) then
- write (*,'("ERROR - unexpected size of ll grid:")')
- write (*,'("ERROR - shape(ll) : ",i4," x ",i4)') shape(ll)
- write (*,'("ERROR - lli%nlat : ",i4)') lli%nlat
- write (*,'("ERROR in ",a)') rname; status=1; return
- end if
- ! loop over latitudes:
- do j = 1, lli%nlat
-
- ! select correct area for cells in this row:
- select case ( unit )
- case ( 'rad2' )
- cell_area = lli%area(j)
- case ( 'm2' )
- cell_area = lli%area_m2(j)
- case default
- write (*,'("ERROR - unknown unit : ",a)') trim(unit)
- write (*,'("ERROR in ",a)') rname; status=1; return
- end select
-
- ! assign/mult/div by cell area:
- select case ( oper )
- case ( '=' )
- ll(:,j) = cell_area
- case ( '/' )
- ll(:,j) = ll(:,j) / cell_area
- case ( '*' )
- ll(:,j) = ll(:,j) * cell_area
- case default
- write (*,'("ERROR - unknown operation : ",a)') trim(oper)
- write (*,'("ERROR in ",a)') rname; status=1; return
- end select
- end do
-
- ! ok
- status = 0
-
- end subroutine llgrid_AreaOper_2d
-
-
- ! *
-
-
- subroutine llgrid_AreaOper_3d( lli, ll, oper, unit, status )
-
- ! --- in/out ----------------------------------
-
- type(TllGridInfo), intent(in) :: lli
- real, intent(inout) :: ll(:,:,:)
- character(len=*), intent(in) :: oper
- character(len=*), intent(in) :: unit
- integer, intent(out) :: status
-
- ! --- const -----------------------------
-
- character(len=*), parameter :: rname = mname//'/llgrid_AreaOper_3d'
-
- ! --- local --------------------------------
-
- integer :: l
-
- ! --- begin ----------------------------------
-
- ! loop over layers
- do l = 1, size(ll,3)
- ! apply 2d operator:
- call AreaOper( lli, ll(:,:,l), oper, unit, status )
- if (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; return; end if
- end do
-
- ! ok
- status = 0
-
- end subroutine llgrid_AreaOper_3d
-
-
- ! =====================================
-
- !
- ! Interpolation to (lon,lat) in deg.
- !
-
- subroutine llgrid_InterpolFractions( lli, lon, lat, ii, jj, ff )
-
- ! --- in/out ---------------------------
-
- type(TllGridInfo), intent(in) :: lli
- real, intent(in) :: lon, lat ! deg
- integer, intent(out) :: ii(4)
- integer, intent(out) :: jj(4)
- real, intent(out) :: ff(4)
-
- ! --- local -----------------------------
-
- real :: lonX, latX
- real :: ir
- integer :: i1, i2
- real :: i1f, i2f
- real :: jr
- integer :: j1, j2
- real :: j1f, j2f
-
- ! --- begin -----------------------------
-
- ! bring lon in [-180,180.0]
- lonX = modulo(lon,360.0)
- if ( lonX > 180.0 ) lonX = lonX - 360.0
-
- ! check lat
- latX = lat
- if ( latX < -90.0 .or. latX > 90.0 ) then
- write (*,'("ERROR - invalid lat (deg) :",f12.4)') latX
- write (*,'("ERROR in ",a)') 'llgrid_InterpolFractions'; stop
- end if
-
- !
- ! 1 2
- ! 3 4
- !
-
- ! i fractions ; circular
- !
- ! *
- ! [----+----][----+----] .. [----+----]
- ! ir 1.0 2.0 120.0
- !
- ir = ( lonX - lli%lon_deg(1) ) / lli%dlon_deg + 1.0
- i1 = floor(ir)
- i2 = i1 + 1
- !
- i2f = ( ir - i1 ) / ( i2 - i1 )
- i1f = 1.0 - i2f
- !
- if ( i1 < 1 ) i1 = lli%nlon
- if ( i2 > lli%nlon ) i2 = 1
-
- ! j fractions ; constant in half cells at poles
- jr = ( latX - lli%lat_deg(1) ) / lli%dlat_deg + 1.0
- j1 = floor(jr)
- j2 = j1 + 1
- !
- if ( j1 < 1 ) then
- j2f = ( jr - 0.5 ) / ( j2 - 0.5 )
- j1f = 1.0 - j2f
- else if ( j2 > lli%nlat ) then
- j2f = ( jr - j1 ) / ( lli%nlat+0.5 - j1 )
- j1f = 1.0 - j2f
- else
- j2f = ( jr - j1 ) / ( j2 - j1 )
- j1f = 1.0 - j2f
- end if
-
- ! fill output
- !
- ! j2 3 4
- ! j1 1 2
- !
- ! i1 i2
- !
- ii = (/ i1, i2, i1, i2 /)
- jj = (/ j1, j1, j2, j2 /)
- ff = (/ i1f*j1f, i2f*j1f, i1f*j2f, i2f*j2f /)
-
- ! write (*,'(" ",2i5)') ii(3), ii(4)
- ! write (*,'(i4," ",f4.2," ",f4.2," ",i4)') jj(3),ff(3),ff(4),jj(4)
- ! write (*,'(i4," ",f4.2," ",f4.2," ",i4)') jj(1),ff(1),ff(2),jj(2)
- ! write (*,'(" ",2i5)') ii(1), ii(2)
-
- end subroutine llgrid_InterpolFractions
-
-
- ! ***
-
-
- subroutine llgrid_Eval_2d( lli, ll, lon, lat, res )
-
- ! --- in/out ---------------------------
-
- type(TllGridInfo), intent(in) :: lli
- real, intent(in) :: ll(:,:)
- real, intent(in) :: lon, lat ! deg
- real, intent(out) :: res
-
- integer :: status
-
- ! --- const -----------------------------
-
- character(len=*), parameter :: rname = mname//'/llgrid_Eval_2d'
-
- ! --- local -----------------------------
-
- integer :: ii(4)
- integer :: jj(4)
- real :: ff(4)
-
- integer :: k
- real :: value
-
- ! --- begin -----------------------------
-
- ! check ...
- call Check( lli, 'n', ll, status )
- if (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; stop; end if
-
- ! indices and fractions
- call InterpolFractions( lli, lon, lat, ii, jj, ff )
-
- ! init zero
- res = 0.0
-
- ! add contributions:
- do k = 1, 4
-
- ! handle poles
- if ( jj(k) < 1 ) then
- value = sum(ll(:,1))/lli%nlon
- else if ( jj(k) > lli%nlat ) then
- value = sum(ll(:,lli%nlat))/lli%nlon
- else
- value = ll(ii(k),jj(k))
- end if
-
- ! add fraction
- res = res + value * ff(k)
-
- end do
-
- end subroutine llgrid_Eval_2d
-
-
- ! ***
- subroutine llgrid_Eval_3d( lli, ll, lon, lat, res )
-
- ! --- in/out ---------------------------
-
- type(TllGridInfo), intent(in) :: lli
- real, intent(in) :: ll(:,:,:)
- real, intent(in) :: lon, lat ! deg
- real, intent(out) :: res(size(ll,3))
-
- integer :: status
-
- ! --- const -----------------------------
-
- character(len=*), parameter :: rname = mname//'/llgrid_Eval_3d'
-
- ! --- local -----------------------------
-
- integer :: ii(4)
- integer :: jj(4)
- real :: ff(4)
-
- integer :: k
- real :: value(size(ll,3))
-
- ! --- begin -----------------------------
-
- ! check ...
- call Check( lli, 'n', ll(:,:,1), status )
- if (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; stop; end if
-
- ! indices and fractions
- call InterpolFractions( lli, lon, lat, ii, jj, ff )
-
- ! init zero
- res = 0.0
-
- ! add contributions:
- do k = 1, 4
-
- ! handle poles
- if ( jj(k) < 1 ) then
- value = sum(ll(:,1,:),1) / lli%nlon
- else if ( jj(k) > lli%nlat ) then
- value = sum(ll(:,lli%nlat,:),1) / lli%nlon
- else
- value = ll(ii(k),jj(k),:)
- end if
-
- ! add fraction
- res = res + value * ff(k)
-
- end do
-
- end subroutine llgrid_Eval_3d
-
-
- ! =================================================================
- ! ===
- ! === match fine grid with coarse grid
- ! ===
- ! =================================================================
-
-
- subroutine GetRefinement( cgi, fgi, &
- refine_i, refine_j, &
- cg_i1, cg_i2, cg_j1, cg_j2, status )
-
- ! --- in/out ------------------------------
-
- type(TllGridInfo), intent(in) :: cgi
- type(TllGridInfo), intent(in) :: fgi
- integer, intent(out) :: refine_i, refine_j
- integer, intent(out) :: cg_i1, cg_i2, cg_j1, cg_j2
- integer, intent(out) :: status
-
- ! --- const -----------------------------
-
- character(len=*), parameter :: rname = mname//'/GetRefinement'
-
- ! --- local -----------------------------
-
- integer :: i, j
-
- ! --- begin -----------------------------
-
- ! *** determine refinement
-
- refine_i = nint( cgi%dlon_deg / fgi%dlon_deg )
- refine_j = nint( cgi%dlat_deg / fgi%dlat_deg )
-
-
- ! *** position of fine grid within coarse grid:
-
- ! search column in coarse grid with same west bound as fine grid:
- cg_i1 = 0
- do i = 1, cgi%im
- if ( cgi%blon_deg(i-1) == fgi%blon_deg(0) ) then
- cg_i1 = i
- exit
- end if
- end do
- if ( cg_i1 < 1 ) then
- write (*,'("ERROR - could not match west bound of fine grid:")')
- write (*,'("ERROR - cgi%blon : ",f12.4)') cgi%blon_deg
- write (*,'("ERROR - target : ",f12.4)') fgi%blon_deg(0)
- write (*,'("ERROR in ",a)') rname; status=1; return
- end if
-
- ! search column in coarse grid with same east bound as fine grid:
- cg_i2 = 0
- do i = 1, cgi%im
- if ( cgi%blon_deg(i) == fgi%blon_deg(fgi%im) ) then
- cg_i2 = i
- exit
- end if
- end do
- if ( cg_i2 < 1 ) then
- write (*,'("ERROR - could not match east bound of fine grid")')
- write (*,'("ERROR - cgi%blon : ",f12.4)') cgi%blon_deg
- write (*,'("ERROR - target : ",f12.4)') fgi%blon_deg(fgi%im)
- write (*,'("ERROR in ",a)') rname; status=1; return
- end if
-
- ! check ...
- if ( (cg_i2-cg_i1+1)*refine_i /= fgi%im ) then
- write (*,'("ERROR - i refinement not ok:")')
- write (*,'("ERROR - coarse cells : ",f12.4)') cg_i1, cg_i2
- write (*,'("ERROR - refinement : ",f12.4)') refine_i
- write (*,'("ERROR - fine cells : ",f12.4)') fgi%im
- write (*,'("ERROR in ",a)') rname; status=1; return
- end if
-
- ! search row in coarse grid with same south bound as fine grid:
- cg_j1 = 0
- do j = 1, cgi%jm
- if ( cgi%blat_deg(j-1) == fgi%blat_deg(0) ) then
- cg_j1 = j
- exit
- end if
- end do
- if ( cg_j1 < 1 ) then
- write (*,'("ERROR - could not match south bound of fine grid")')
- write (*,'("ERROR in ",a)') rname; status=1; return
- end if
-
- ! search row in coarse grid with same south bound as fine grid:
- cg_j2 = 0
- do j = 1, cgi%jm
- if ( cgi%blat_deg(j) == fgi%blat_deg(fgi%jm) ) then
- cg_j2 = j
- exit
- end if
- end do
- if ( cg_j2 < 1 ) then
- write (*,'("ERROR - could not match north bound of fine grid")')
- write (*,'("ERROR in ",a)') rname; status=1; return
- end if
-
- ! check ...
- if ( (cg_j2-cg_j1+1)*refine_j /= fgi%jm ) then
- write (*,'("ERROR - j refinement not ok:")')
- write (*,'("ERROR - coarse cells : ",2i5)') cg_j1, cg_j2
- write (*,'("ERROR - refinement : ",i5)') refine_j
- write (*,'("ERROR - fine cells : ",i5)') fgi%jm
- write (*,'("ERROR in ",a)') rname; status=1; return
- end if
-
- ! ok
- status = 0
-
- end subroutine GetRefinement
-
-
- ! ***
-
- ! Relate fine grid and coarse grid:
- !
- ! call Match_cell( action, key, cgi, cg, fgi, fg )
- !
- ! 1. The coarse grid covers the same area as the fine grid:
- !
- ! action:
- ! 'combine' : fill coarse grid by combining cells of fine grid
- ! key:
- ! 'sum' : sum values in fine grid
- ! 'aver' : aver values in fine grid
- ! 'area-aver' : idem, weighted with area
- !
- ! 2. The fine grid is a subset of the coarse grid:
- !
- ! action:
- ! 'subset' : fill fine grid as a subset of coarse grid
- ! key:
- ! not used
- !
- ! 3. The fine grid is a zooming area of the coarse grid:
- !
- ! action:
- ! 'match' : adjust values in fine grid to match coarse grid
- ! key:
- ! 'sum' : sum values in fine grid
- ! 'aver' : aver values in fine grid
- ! 'area-aver' : idem, weighted with area
- !
-
- subroutine llgrid_Match( key, nuv, pgi, pg, tgi, tg, status )
-
- ! --- in/out ------------------------------
-
- ! character(len=*), intent(in) :: action
- character(len=*), intent(in) :: key
- character(len=1), intent(in) :: nuv
- type(TllGridInfo), intent(in) :: pgi
- real, intent(in) :: pg(pgi%im,pgi%jm)
- type(TllGridInfo), intent(in) :: tgi
- real, intent(inout) :: tg(tgi%im,tgi%jm)
- integer, intent(out) :: status
-
- ! --- const ---------------------------------
- character(len=*), parameter :: rname = mname//'/llgrid_Match'
-
- ! --- begin -----------------------------
-
- ! call nuv specific routine
- select case ( nuv )
- case ( 'n' )
- call Match_cell( key, pgi, pg, tgi, tg, status )
- if (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; return; end if
- case ( 'u' )
- call Match_u( key, pgi, pg, tgi, tg, status )
- if (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; return; end if
- case ( 'v' )
- call Match_v( key, pgi, pg, tgi, tg, status )
- if (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; return; end if
- case default
- write (*,'("ERROR - unsupported nuv `",a,"`")') nuv
- write (*,'("ERROR in ",a)') rname; status=1; return
- end select
-
- ! ok
- status = 0
-
- end subroutine llgrid_Match
-
- ! ***
- subroutine Match_cell( key, pgi, pg, tgi, tg, status )
-
- ! --- in/out ------------------------------
-
- ! character(len=*), intent(in) :: action
- character(len=*), intent(in) :: key
- type(TllGridInfo), intent(in) :: pgi
- real, intent(in) :: pg(pgi%im,pgi%jm)
- type(TllGridInfo), intent(in) :: tgi
- real, intent(inout) :: tg(tgi%im,tgi%jm)
- integer, intent(out) :: status
-
- ! --- const ---------------------------------
- character(len=*), parameter :: rname = mname//'/Match_cell'
-
- ! --- local -----------------------------
-
- integer :: refine_i, refine_j
- integer :: cg_i1, cg_i2, cg_j1, cg_j2
- integer :: ci, cj
- integer :: fg_i1, fg_i2, fg_j1, fg_j2
- integer :: fi, fj
- real :: fsum
-
- ! --- begin -----------------------------
-
- ! select case ( action )
- !
- ! !
- ! ! match fine grid with coarse cells:
- ! !
- ! case ( 'match' )
- ! determine refinement
- call GetRefinement( pgi, tgi, refine_i, refine_j, cg_i1, cg_i2, cg_j1, cg_j2, status )
- if (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; return; end if
- ! loop over cells in coarse grid covering fine grid:
- do cj = cg_j1, cg_j2
- do ci = cg_i1, cg_i2
- fg_i1 = (ci-cg_i1)*refine_i + 1 ; fg_i2 = fg_i1-1 + refine_i
- fg_j1 = (cj-cg_j1)*refine_j + 1 ; fg_j2 = fg_j1-1 + refine_j
- ! sum over cells in fine grid:
- select case ( key )
- case ( 'sum' )
- fsum = sum( tg(fg_i1:fg_i2,fg_j1:fg_j2) )
- ! distribute difference equally over all cells in fine grid:
- ! (/6,4/) + ( 14 - 10 )/ 2 = (/8,6/)
- tg(fg_i1:fg_i2,fg_j1:fg_j2) = tg(fg_i1:fg_i2,fg_j1:fg_j2) + (pg(ci,cj)-fsum)/(refine_j*refine_i)
- ! cmk corrected: divide by (refine_i * refine_j)
- case ( 'aver' )
- fsum = sum( tg(fg_i1:fg_i2,fg_j1:fg_j2) )/(refine_i*refine_j)
- ! add difference in averages to all cells in fine grid:
- ! (/6,4/) + ( 7 - 5 ) = (/8,6/)
- tg(fg_i1:fg_i2,fg_j1:fg_j2) = tg(fg_i1:fg_i2,fg_j1:fg_j2) + (pg(ci,cj)-fsum)
- case ( 'area-aver' )
- fsum = 0.0
- do fj = fg_j1, fg_j2
- fsum = fsum + sum(tg(fg_i1:fg_i2,fj))*tgi%area(fj)
- end do
- fsum = fsum / pgi%area(cj)
- ! add difference in averages to all cells in fine grid:
- ! (/6,4/) + ( 7 - 5 ) = (/8,6/)
- tg(fg_i1:fg_i2,fg_j1:fg_j2) = tg(fg_i1:fg_i2,fg_j1:fg_j2) + (pg(ci,cj)-fsum)
- case default
- write (*,'("ERROR - unsupported key for match action:")')
- ! write (*,'("ERROR - action : ",a)') trim(action)
- write (*,'("ERROR - key : ",a)') trim(key)
- write (*,'("ERROR in ",a)') rname; status=1; return
- end select
- ! ! match fine grid with coarse cell:
- ! if ( fsum == 0.0 ) then
- ! if ( abs(pg(ci,cj)) > 1.0e-5 ) then
- ! write (*,'("ERROR - zero sum over coarse cell (",i3,",",i3,")")') ci, cj
- ! write (*,'("ERROR - coarse cell : ",es12.4)') pg(ci,cj)
- ! write (*,'("ERROR - refine : ",2i6)') refine_i, refine_j
- ! write (*,'("ERROR - fine grid : ",es12.4)') tg(fg_i1:fg_i2,fg_j1:fg_j2)
- ! write (*,'("ERROR in ",a)') rname; status=1; return
- ! end if
- ! else
- ! ! scale towards value of coarse grid:
- ! tg(fg_i1:fg_i2,fg_j1:fg_j2) = tg(fg_i1:fg_i2,fg_j1:fg_j2) * pg(ci,cj)/fsum
- ! end if
- end do
- end do
-
- ! !
- ! ! fine grid is subset of coarse grid:
- ! !
- ! case ( 'subset' )
- !
- ! ! determine refinement
- ! call GetRefinement( pgi, tgi, refine_i, refine_j, cg_i1, cg_i2, cg_j1, cg_j2 )
- !
- ! if ( (refine_i /= 1) .or. (refine_j /= 1) ) then
- ! write (*,'("ERROR - for this action the fine grid should a subset of a coarse grid:")')
- ! write (*,'("ERROR - action : ",a)') trim(action)
- ! write (*,'("ERROR - refinement : ",2i6)') refine_i, refine_j
- ! write (*,'("ERROR in ",a)') rname; status=1; return
- ! end if
- !
- ! ! loop over cells in coarse grid covering fine grid:
- ! do cj = cg_j1, cg_j2
- ! do ci = cg_i1, cg_i2
- !
- ! fg_i1 = (ci-cg_i1) + 1
- ! fg_j1 = (cj-cg_j1) + 1
- !
- ! ! fine grid is subset of coarse grid
- ! tg(fg_i1,fg_j1) = pg(ci,cj)
- !
- ! end do
- ! end do
- !
- ! !
- ! ! collect cells in fine grid to coarse grid
- ! !
- ! case ( 'combine' )
- !
- ! ! determine refinement
- ! ! NOTE: parent grid is fine now, target is coarse !
- ! call GetRefinement( tgi, pgi, refine_i, refine_j, cg_i1, cg_i2, cg_j1, cg_j2 )
- !
- ! if ( (cg_i1 /= 1) .or. (cg_i2 /= tgi%im) .or. &
- ! (cg_j1 /= 1) .or. (cg_j2 /= tgi%jm) ) then
- ! write (*,'("ERROR - for this action the fine grid should cover the complete coarse grid:")')
- ! write (*,'("ERROR - action : ",a)') trim(action)
- ! write (*,'("ERROR - covered : [",i4,",",i4,"] x [",i4,",","]")') cg_i1,cg_i2,cg_j1,cg_j2
- ! write (*,'("ERROR - of : ",2i5)') tgi%im, tgi%jm
- ! write (*,'("ERROR in ",a)') rname; status=1; return
- ! end if
- !
- ! ! loop over cells in coarse grid covering fine grid:
- ! do cj = cg_j1, cg_j2
- ! do ci = cg_i1, cg_i2
- !
- ! fg_i1 = (ci-cg_i1)*refine_i + 1 ; fg_i2 = fg_i1-1 + refine_i
- ! fg_j1 = (cj-cg_j1)*refine_j + 1 ; fg_j2 = fg_j1-1 + refine_j
- !
- ! ! sum over cells in fine grid:
- ! select case ( key )
- ! case ( 'sum' )
- ! fsum = sum( pg(fg_i1:fg_i2,fg_j1:fg_j2) )
- ! case ( 'aver' )
- ! fsum = sum( pg(fg_i1:fg_i2,fg_j1:fg_j2) )/(refine_i*refine_j)
- ! case ( 'area-aver' )
- ! fsum = 0.0
- ! do fj = fg_j1, fg_j2
- ! fsum = fsum + sum(pg(fg_i1:fg_i2,fj))*pgi%area(fj)
- ! end do
- ! fsum = fsum / tgi%area(cj)
- ! case default
- ! write (*,'("ERROR - unsupported key for match action:")')
- ! write (*,'("ERROR - action : ",a)') trim(action)
- ! write (*,'("ERROR - key : ",a)') trim(key)
- ! write (*,'("ERROR in ",a)') rname; status=1; return
- ! end select
- !
- ! ! collect cells in fine parent grid to target coarse grid
- ! tg(ci,cj) = fsum
- !
- ! end do
- ! end do
- !
- ! case default
- ! write (*,'("ERROR - unknown action `",a,"`")') trim(action)
- ! write (*,'("ERROR in ",a)') rname; status=1; return
- ! end select
-
- ! ok
- status = 0
-
- end subroutine Match_cell
-
- ! ***
-
- ! flux through east/west boundaries
- subroutine Match_u( key, pgi, pg, tgi, tg, status )
-
- ! --- in/out ------------------------------
-
- ! character(len=*), intent(in) :: action
- character(len=*), intent(in) :: key
- type(TllGridInfo), intent(in) :: pgi
- real, intent(in) :: pg(0:pgi%im,pgi%jm)
- type(TllGridInfo), intent(in) :: tgi
- real, intent(inout) :: tg(0:tgi%im,tgi%jm)
- integer, intent(out) :: status
-
- ! --- const ---------------------------------
- character(len=*), parameter :: rname = mname//'/Match_u'
-
- ! --- local -----------------------------
-
- integer :: refine_i, refine_j
- integer :: cg_i1, cg_i2, cg_j1, cg_j2
- integer :: ci, cj
- integer :: fg_i, fg_j1, fg_j2
- integer :: fi, fj
- real :: fsum
-
- ! --- begin -----------------------------
-
- ! select case ( action )
- !
- ! !
- ! ! match fine grid with coarse cells:
- ! !
- ! case ( 'match' )
- ! determine refinement
- call GetRefinement( pgi, tgi, refine_i, refine_j, cg_i1, cg_i2, cg_j1, cg_j2, status )
- if (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; return; end if
- ! loop over cells in coarse grid covering fine grid:
- do cj = cg_j1, cg_j2
- do ci = cg_i1-1, cg_i2
- fg_i = (ci-(cg_i1-1))*refine_i
- fg_j1 = (cj-cg_j1)*refine_j + 1 ; fg_j2 = fg_j1-1 + refine_j
- ! sum over cells in fine grid:
- select case ( key )
- case ( 'sum' )
- fsum = sum( tg(fg_i,fg_j1:fg_j2) )
- ! distribute difference equally over all cells in fine grid:
- ! (/6,4/) + ( 14 - 10 )/ 2 = (/8,6/)
- tg(fg_i,fg_j1:fg_j2) = tg(fg_i,fg_j1:fg_j2) + (pg(ci,cj)-fsum)/refine_j
- case ( 'aver' )
- fsum = sum( tg(fg_i,fg_j1:fg_j2) )/(refine_j)
- ! add difference in averages to all cells in fine grid:
- ! (/6,4/) + ( 7 - 5 ) = (/8,6/)
- tg(fg_i,fg_j1:fg_j2) = tg(fg_i,fg_j1:fg_j2) + (pg(ci,cj)-fsum)
- case default
- write (*,'("ERROR - unsupported key for match action:")')
- ! write (*,'("ERROR - action : ",a)') trim(action)
- write (*,'("ERROR - key : ",a)') trim(key)
- write (*,'("ERROR in ",a)') rname; status=1; return
- end select
- ! ! match fine grid with coarse cell:
- ! if ( fsum == 0.0 ) then
- ! if ( pg(ci,cj) /= 0.0 ) then
- ! write (*,'("ERROR - zero sum over coarse cell (",i3,",",i3,")")') ci, cj
- ! write (*,'("ERROR - coarse cell : ",es12.4)') pg(ci,cj)
- ! write (*,'("ERROR - refine : ",2i6)') refine_i, refine_j
- ! write (*,'("ERROR - target grid : ",es12.4)') tg(fg_i,fg_j1:fg_j2)
- ! write (*,'("ERROR in ",a)') rname; status=1; return
- ! end if
- ! else
- ! ! scale towards value of coarse grid:
- ! tg(fg_i,fg_j1:fg_j2) = tg(fg_i,fg_j1:fg_j2)/fsum * pg(ci,cj)
- ! end if
- end do
- end do
-
- ! !
- ! ! fine grid is subset of coarse grid:
- ! !
- ! case ( 'subset' )
- !
- ! ! determine refinement
- ! call GetRefinement( pgi, tgi, refine_i, refine_j, cg_i1, cg_i2, cg_j1, cg_j2 )
- !
- ! if ( (refine_i /= 1) .or. (refine_j /= 1) ) then
- ! write (*,'("ERROR - for this action the fine grid should a subset of a coarse grid:")')
- ! write (*,'("ERROR - action : ",a)') trim(action)
- ! write (*,'("ERROR - refinement : ",2i6)') refine_i, refine_j
- ! write (*,'("ERROR in ",a)') rname; status=1; return
- ! end if
- !
- ! ! loop over cells in coarse grid covering fine grid:
- ! do cj = cg_j1, cg_j2
- ! do ci = cg_i1-1, cg_i2
- !
- ! fg_i = (ci-(cg_i1-1))
- ! fg_j1 = (cj-cg_j1) + 1
- !
- ! ! fine grid is subset of coarse grid
- ! tg(fg_i,fg_j1) = pg(ci,cj)
- !
- ! end do
- ! end do
- !
- ! !
- ! ! collect cells in fine grid to coarse grid
- ! !
- ! case ( 'combine' )
- !
- ! ! determine refinement
- ! ! NOTE: parent grid is fine, targer is coarse !
- ! call GetRefinement( tgi, pgi, refine_i, refine_j, cg_i1, cg_i2, cg_j1, cg_j2 )
- !
- ! if ( (cg_i1 /= 1) .or. (cg_i2 /= tgi%im) .or. &
- ! (cg_j1 /= 1) .or. (cg_j2 /= tgi%jm) ) then
- ! write (*,'("ERROR - for this action the fine grid should cover the complete coarse grid:")')
- ! write (*,'("ERROR - action : ",a)') trim(action)
- ! write (*,'("ERROR - covered : [",i4,",",i4,"] x [",i4,",","]")') cg_i1,cg_i2,cg_j1,cg_j2
- ! write (*,'("ERROR - of : ",2i5)') tgi%im, tgi%jm
- ! write (*,'("ERROR in ",a)') rname; status=1; return
- ! end if
- !
- ! ! loop over cells in coarse grid covering fine grid:
- ! do cj = cg_j1, cg_j2
- ! do ci = cg_i1-1, cg_i2
- !
- ! fg_i = (ci-(cg_i1-1))*refine_i
- ! fg_j1 = (cj-cg_j1)*refine_j + 1 ; fg_j2 = fg_j1-1 + refine_j
- !
- ! ! sum over cells in fine grid:
- ! select case ( key )
- ! case ( 'sum' )
- ! fsum = sum( pg(fg_i,fg_j1:fg_j2) )
- ! case ( 'aver' )
- ! fsum = sum( pg(fg_i,fg_j1:fg_j2) )/(refine_j)
- ! case default
- ! write (*,'("ERROR - unsupported key for match action:")')
- ! write (*,'("ERROR - action : ",a)') trim(action)
- ! write (*,'("ERROR - key : ",a)') trim(key)
- ! write (*,'("ERROR in ",a)') rname; status=1; return
- ! end select
- !
- ! ! collect cells in fine parent grid to target coarse grid
- ! tg(ci,cj) = fsum
- !
- ! end do
- ! end do
- !
- ! case default
- ! write (*,'("ERROR - unknown action `",a,"`")') trim(action)
- ! write (*,'("ERROR in ",a)') rname; status=1; return
- ! end select
-
- ! ok
- status = 0
-
- end subroutine Match_u
-
- ! ***
-
- ! flux through north/south boundaries
- subroutine Match_v( key, pgi, pg, tgi, tg, status )
-
- ! --- in/out ------------------------------
-
- ! character(len=*), intent(in) :: action
- character(len=*), intent(in) :: key
- type(TllGridInfo), intent(in) :: pgi
- real, intent(in) :: pg(pgi%im,0:pgi%jm)
- type(TllGridInfo), intent(in) :: tgi
- real, intent(inout) :: tg(tgi%im,0:tgi%jm)
- integer, intent(out) :: status
-
- ! --- const ---------------------------------
- character(len=*), parameter :: rname = mname//'/Match_v'
-
- ! --- local -----------------------------
-
- integer :: refine_i, refine_j
- integer :: cg_i1, cg_i2, cg_j1, cg_j2
- integer :: ci, cj
- integer :: fg_i1, fg_i2, fg_j
- integer :: fi, fj
- real :: fsum
-
- ! --- begin -----------------------------
-
- ! select case ( action )
- !
- ! !
- ! ! match fine grid with coarse cells:
- ! !
- ! case ( 'match' )
- ! determine refinement
- call GetRefinement( pgi, tgi, refine_i, refine_j, cg_i1, cg_i2, cg_j1, cg_j2, status )
- if (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; return; end if
- ! loop over cells in coarse grid covering fine grid:
- do cj = cg_j1-1, cg_j2
- do ci = cg_i1, cg_i2
- fg_i1 = (ci-cg_i1)*refine_i + 1 ; fg_i2 = fg_i1-1 + refine_i
- fg_j = (cj-(cg_j1-1))*refine_j
- ! sum over cells in fine grid:
- select case ( key )
- case ( 'sum' )
- fsum = sum( tg(fg_i1:fg_i2,fg_j) )
- ! distribute difference equally over all cells in fine grid:
- ! (/6,4/) + ( 14 - 10 )/ 2 = (/8,6/)
- tg(fg_i1:fg_i2,fg_j) = tg(fg_i1:fg_i2,fg_j) + (pg(ci,cj)-fsum)/refine_i
- case ( 'aver' )
- fsum = sum( tg(fg_i1:fg_i2,fg_j) )/(refine_i)
- ! add difference in averages to all cells in fine grid:
- ! (/6,4/) + ( 7 - 5 ) = (/8,6/)
- tg(fg_i1:fg_i2,fg_j) = tg(fg_i1:fg_i2,fg_j) + (pg(ci,cj)-fsum)
- case default
- write (*,'("ERROR - unsupported key for match action:")')
- ! write (*,'("ERROR - action : ",a)') trim(action)
- write (*,'("ERROR - key : ",a)') trim(key)
- write (*,'("ERROR in ",a)') rname; status=1; return
- end select
- ! ! match fine grid with coarse cell:
- ! if ( fsum == 0.0 ) then
- ! if ( pg(ci,cj) /= 0.0 ) then
- ! write (*,'("ERROR - zero sum over coarse cell (",i3,",",i3,")")') ci, cj
- ! write (*,'("ERROR - coarse cell : ",es12.4)') pg(ci,cj)
- ! write (*,'("ERROR - refine : ",2i6)') refine_i, refine_j
- ! write (*,'("ERROR - fine grid : ",es12.4)') tg(fg_i1:fg_i2,fg_j)
- ! write (*,'("ERROR in ",a)') rname; status=1; return
- ! end if
- ! else
- ! ! scale towards value of coarse grid:
- ! tg(fg_i1:fg_i2,fg_j) = tg(fg_i1:fg_i2,fg_j)/fsum * pg(ci,cj)
- ! end if
- end do
- end do
-
- ! !
- ! ! fine grid is subset of coarse grid:
- ! !
- ! case ( 'subset' )
- !
- ! ! determine refinement
- ! call GetRefinement( pgi, tgi, refine_i, refine_j, cg_i1, cg_i2, cg_j1, cg_j2 )
- !
- ! if ( (refine_i /= 1) .or. (refine_j /= 1) ) then
- ! write (*,'("ERROR - for this action the fine grid should a subset of a coarse grid:")')
- ! write (*,'("ERROR - action : ",a)') trim(action)
- ! write (*,'("ERROR - refinement : ",2i6)') refine_i, refine_j
- ! write (*,'("ERROR in ",a)') rname; status=1; return
- ! end if
- !
- ! ! loop over cells in coarse grid covering fine grid:
- ! do cj = cg_j1-1, cg_j2
- ! do ci = cg_i1, cg_i2
- !
- ! fg_i1 = (ci-cg_i1) + 1
- ! fg_j = (cj-(cg_j1-1))
- !
- ! ! fine grid is subset of coarse grid
- ! tg(fg_i1,fg_j) = pg(ci,cj)
- !
- ! end do
- ! end do
- !
- ! !
- ! ! collect cells in fine grid to coarse grid
- ! !
- ! case ( 'combine' )
- !
- ! ! determine refinement
- ! ! NOTE: parent grid is fine, target is coarse !
- ! call GetRefinement( tgi, pgi, refine_i, refine_j, cg_i1, cg_i2, cg_j1, cg_j2 )
- !
- ! if ( (cg_i1 /= 1) .or. (cg_i2 /= tgi%im) .or. &
- ! (cg_j1 /= 1) .or. (cg_j2 /= tgi%jm) ) then
- ! write (*,'("ERROR - for this action the fine grid should cover the complete coarse grid:")')
- ! write (*,'("ERROR - action : ",a)') trim(action)
- ! write (*,'("ERROR - covered : [",i4,",",i4,"] x [",i4,",","]")') cg_i1,cg_i2,cg_j1,cg_j2
- ! write (*,'("ERROR - of : ",2i5)') tgi%im, tgi%jm
- ! write (*,'("ERROR in ",a)') rname; status=1; return
- ! end if
- !
- ! ! loop over cells in coarse grid covering fine grid:
- ! do cj = cg_j1-1, cg_j2
- ! do ci = cg_i1, cg_i2
- !
- ! fg_i1 = (ci-cg_i1)*refine_i + 1 ; fg_i2 = fg_i1-1 + refine_i
- ! fg_j = (cj-(cg_j1-1))*refine_j
- !
- ! ! sum over cells in fine grid:
- ! select case ( key )
- ! case ( 'sum' )
- ! fsum = sum( pg(fg_i1:fg_i2,fg_j) )
- ! case ( 'aver' )
- ! fsum = sum( pg(fg_i1:fg_i2,fg_j) )/(refine_i)
- ! case default
- ! write (*,'("ERROR - unsupported key for match action:")')
- ! write (*,'("ERROR - action : ",a)') trim(action)
- ! write (*,'("ERROR - key : ",a)') trim(key)
- ! write (*,'("ERROR in ",a)') rname; status=1; return
- ! end select
- !
- ! ! collect cells in fine parent grid to target coarse grid
- ! tg(ci,cj) = fsum
- !
- ! end do
- ! end do
- !
- ! case default
- ! write (*,'("ERROR - unknown action `",a,"`")') trim(action)
- ! write (*,'("ERROR in ",a)') rname; status=1; return
- ! end select
-
- ! ok
- status = 0
-
- end subroutine Match_v
-
-
- ! ========================================================
- ! ===
- ! === fill grid from other grid
- ! ===
- ! ========================================================
-
- !
- ! NUV
- !
- ! Key to identify data positions:
- ! 'n' : value valid for cell (center) ll(1:nlon ,1:nlat )
- ! 'u' : value valid for east/west boundaries ll(1:nlon+1,1:nlat )
- ! 'v' : value valid for north/south boundaries ll(1:nlon ,1:nlat+1)
- !
- ! ROUTINES
- !
- ! call FillGrid( lli, nuv, ll, lliX, nuvX, llX, combkey, status [,llX_w] )
- !
- ! Fill ll (defined by lli,nuv) with values from llX (defined by lliX,nuvX)
- !
- ! Coverage of lli by lliX :
- ! o lliX is larger than or equal to lli -> all cells in ll changed
- ! o lliX is smaller than lli -> only part of ll is changed
- !
- ! Create new ll from llX:
- ! o llX is superset -> copy values from llX into ll
- ! o llX is fine -> fill ll by combining cells in llX
- ! (average/sum/etc given the combine key)
- !
- ! Combine keys:
- !
- ! 'sum' : sum_i llX_i
- !
- ! 'aver' : sum_i llX_i / sum_i i
- !
- ! 'area-aver' : sum_i llX_i A_i / sum_i A_i
- !
- ! 'weight' : sum_i llX_i llX_w_i / sum_i llX_w_i
- ! (only for nuv='n')
- !
- ! Return status:
- ! 0 : ok
- ! -1 : only part of ll is filled
- !
- !
- ! AdjustGrid NOT IMPLEMENTED YET
- !
- ! Adjust ll given llX:
- ! o llX is coarse -> adjust ll such that average/sum/.. of lli matches
- !
- ! 15 Dec 2010 - P. Le Sager - added fix for 'mass-aver' (always
- ! understood as area-averaged)
- !
- subroutine FillGrid( lli, nuv, ll, lliX, nuvX, llX, combkey, status, llX_w )
- use dims, only : okdebug
- use GO, only : gol, goPr
-
- ! --- in/out --------------------------------
- type(TllGridInfo), intent(in) :: lli
- character(len=*), intent(in) :: nuv
- real, intent(out) :: ll(:,:)
- type(TllGridInfo), intent(in) :: lliX
- character(len=*), intent(in) :: nuvX
- real, intent(in) :: llX(:,:)
- character(len=*), intent(in) :: combkey
- integer, intent(out) :: status
- real, intent(in), optional :: llX_w(:,:)
- ! --- const ---------------------------------
- character(len=*), parameter :: rname = mname//'/FillGrid'
-
- ! --- local ---------------------------------
- character(len=10) :: action
- integer :: di, dj
- integer :: i1, i2, j1, j2
- integer :: js, je
- integer :: i, j
- integer :: i1X, i2X, j1X, j2X
- integer :: iX, jX
- integer :: diX, djX, nX
- real :: res, resw
- integer :: ia, ib, ja, jb
- integer :: iaX, ibX, jaX, jbX
- real :: llX_ab
-
- real, allocatable :: wwX(:,:)
- logical :: wwdiv
- logical :: only_part_of_ll
- ! --- begin ---------------------------------
-
- ! check input ...
- if ( nuv /= nuvX ) then
- write (*,'("ERROR - nuv keys should be equal:")') combkey
- write (*,'("ERROR - nuv : `",a,"`")') nuv
- write (*,'("ERROR - nuvX : `",a,"`")') nuvX
- write (*,'("ERROR in ",a)') rname; status=1; return
- end if
- ! determine how lliX is related to lli:
- ! i1, i2, j1, j2 : cell ranges in lli covered by cells of lliX
- ! i1X, i2X, j1X, j2X : cell ranges in lliX covering cells of lliX
- ! action : 'copy', 'combine'
- call Relate( lli , i1 , i2 , j1 , j2 , &
- lliX, i1X, i2X, j1X, j2X, &
- action, status )
- select case ( status )
- case ( -1 )
- only_part_of_ll = .true.
- case ( 0 )
- only_part_of_ll = .false.
- case default
- write (*,'("ERROR in ",a)') rname; return
- end select
- ! IF(okdebug)THEN
- ! WRITE(gol,'(a," : action = ",a)') rname, action; CALL goPr
- ! ENDIF
-
- ! what to do with cells in llX?
- select case ( action )
- !
- ! * copy
- !
- ! 'copy' : lli and lliX define same resolution;
- ! the cells from llX area [i1X,i2X] x [j1X,j2X]
- ! should be copied into ll area [i1,i2] x [j1,j2]
- case ( 'copy' )
- select case ( nuv )
-
- case ( 'n' )
- ! loop over (selection of) cells of target grid lli:
- !$OMP PARALLEL &
- !$OMP default (none) &
- !$OMP shared (i1, i2, j1, j2, i1x, j1x, llx, ll) &
- !$OMP private (i, j, js, je, ix, jx)
- do j = j1, j2
- do i = i1, i2
- ! source cell in llX:
- iX = i1X + i-i1
- jX = j1X + j-j1
- ! copy cell:
- ll(i,j) = llX(iX,jX)
- end do
- end do
- !$OMP END PARALLEL
- case ( 'u' )
- ! loop over (selection of) cells of target grid lli:
- do j = j1, j2
- do i = i1, i2+1
- ! source cell in llX:
- iX = i1X + i-i1
- jX = j1X + j-j1
- ! copy cell:
- ll(i,j) = llX(iX,jX)
- end do
- end do
- case ( 'v' )
- ! loop over (selection of) cells of target grid lli:
- do j = j1, j2+1
- do i = i1, i2
- ! source cell in llX:
- iX = i1X + i-i1
- jX = j1X + j-j1
- ! copy cell:
- ll(i,j) = llX(iX,jX)
- end do
- end do
- case default
- write (*,'("ERROR - unsupported nuv `",a,"`")') nuv
- write (*,'("ERROR - action : `",a,"`")') action
- write (*,'("ERROR in ",a)') rname; status=1; return
- end select
- !
- ! * combine
- !
- ! 'combine' : lliX defines a fine resolution;
- ! the cells from llX area [i1X,i2X] x [j1X,j2X]
- ! should be combined and copied into ll area [i1,i2] x [j1,j2]
- case ( 'combine' )
-
- ! resolution of fine cells (lliX) in cell of lli,
- ! for example 3 x 2 :
- diX = (i2X-i1X+1) / (i2-i1+1)
- djX = (j2X-j1X+1) / (j2-j1+1)
- nX = diX * djX
- select case ( nuv )
-
- case ( 'n' )
- !
- ! set weight:
- !
- ! 'weight' : (sum_i f_i w_i) / (sum_i w_i)
- !
- ! 'sum' : (sum_i f_i) : w = 1.0
- ! 'aver' : (sum_i f_i) / (sum_i) : w = 1.0, wwdiv
- ! 'area-aver' : (sum_i f_i A_i) / (sum_i A_i) : w = A_i, wwdiv
- ! same size as input grid:
- allocate( wwX(size(llX,1),size(llX,2)) )
- ! weight provided as argument ?
- if ( combkey == 'weight' ) then
- if ( .not. present(llX_w) ) then
- write (*,'("ERROR - combkey `weight` but llX_w not present ...")')
- write (*,'("ERROR in ",a)') rname; status=1; return
- end if
- if ( any( shape(llX_w) /= shape(llX) ) ) then
- write (*,'("ERROR - weight should have shape of input grid:")')
- write (*,'("ERROR - shape(llX_w) : `",2i6)') shape(llX_w)
- write (*,'("ERROR - shape(llX) : `",2i6)') shape(llX)
- write (*,'("ERROR in ",a)') rname; status=1; return
- end if
- wwX = llX_w
- wwdiv = .true.
- else
- if ( present(llX_w) ) then
- write (*,'("ERROR - llX_w pressent but no combkey `weight` ...")')
- write (*,'("ERROR in ",a)') rname; status=1; return
- end if
- ! fill weight given combkey:
- select case ( combkey )
- case ( 'sum' )
- wwX = 1.0
- wwdiv = .false.
- case ( 'aver' )
- wwX = 1.0
- wwdiv = .true.
- case ( 'area-aver','mass-aver' )
- call AreaOper( lliX, wwX, '=', 'm2', status )
- if (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; return; end if
- wwdiv = .true.
- case default
- write (*,'("ERROR - unsupported combkey `",a,"`")') combkey
- write (*,'("ERROR - nuv : `",a,"`")') nuv
- write (*,'("ERROR - action : `",a,"`")') action
- write (*,'("ERROR in ",a)') rname; status=1; return
- end select
- end if
- ! loop over (selection of) cells of target grid lli:
- !$OMP PARALLEL &
- !$OMP default (none) &
- !$OMP shared (i1, i2, j1, j2, i1x, j1x, djx, dix, llx, wwx, wwdiv, &
- !$OMP ll) &
- !$OMP private (i, j, js, je, res, resw, ix, jx)
- do j = j1, j2
- do i = i1, i2
- ! start with zero result:
- res = 0.0
- resw = 0.0
- ! loop over source cells in llX:
- do jX = j1X + (j-j1)*djX, j1X + (j-j1+1)*djX-1
- do iX = i1X + (i-i1)*diX, i1X + (i-i1+1)*diX-1
- res = res + llX(iX,jX) * wwX(iX,jX)
- resw = resw + wwX(iX,jX)
- end do
- end do
- ! fill result:
- if ( wwdiv ) then
- ll(i,j) = res / resw
- else
- ll(i,j) = res
- end if
- end do
- end do
- !$OMP END PARALLEL
-
- ! clear
- deallocate( wwX )
- case ( 'u' )
- ! loop over (selection of) cells of target grid lli:
- do j = j1, j2
- do i = i1, i2+1
- ! start with zero result:
- res = 0.0
- ! loop over points on east/west boundary in llX:
- iX = i1X + (i-i1)*diX
- do jX = j1X + (j-j1)*djX, j1X + (j-j1+1)*djX-1
- ! add contribution of this llX cell:
- select case ( combkey )
- case ( 'sum' )
- res = res + llX(iX,jX)
- case ( 'aver', 'area-aver','mass-aver' )
- res = res + llX(iX,jX)/real(djX)
- case default
- write (*,'("ERROR - unsupported combkey `",a,"`")') combkey
- write (*,'("ERROR - nuv : `",a,"`")') nuv
- write (*,'("ERROR - action : `",a,"`")') action
- write (*,'("ERROR in ",a)') rname; status=1; return
- end select
- end do
- ! fill result:
- ll(i,j) = res
- end do
- end do
-
- case ( 'v' )
- ! loop over (selection of) cells of target grid lli:
- do j = j1, j2+1
- do i = i1, i2
- ! start with zero result:
- res = 0.0
- ! loop over points on north/south boundary in llX:
- jX = j1X + (j-j1)*djX
- do iX = i1X + (i-i1)*diX, i1X + (i-i1+1)*diX-1
- ! add contribution of this llX cell:
- select case ( combkey )
- case ( 'sum' )
- res = res + llX(iX,jX)
- case ( 'aver', 'area-aver','mass-aver' )
- res = res + llX(iX,jX)/real(diX)
- case default
- write (*,'("ERROR - unsupported combkey `",a,"`")') combkey
- write (*,'("ERROR - nuv : `",a,"`")') nuv
- write (*,'("ERROR - action : `",a,"`")') action
- write (*,'("ERROR in ",a)') rname; status=1; return
- end select
- end do
- ! fill result:
- ll(i,j) = res
- end do
- end do
-
- case default
- write (*,'("ERROR - unsupported nuv `",a,"`")') nuv
- write (*,'("ERROR - action : `",a,"`")') action
- write (*,'("ERROR in ",a)') rname; status=1; return
- end select
- !
- ! * distribute
- !
- ! lliX defines a coarse resolution;
- ! the cells from llX area [i1X,i2X] x [j1X,j2X]
- ! should be sampled onto ll area [i1,i2] x [j1,j2]
- !
- ! Note: it is not possible to produce weighted distributions
- ! to have different values within a coarse grid,
- ! for example based on area, since what to do with a cell with
- ! zero area ?
-
- case ( 'distribute' )
-
- ! resolution of fine cells (lli) in cell of lliX,
- ! for example 3 x 2 :
- di = (i2-i1+1) / (i2X-i1X+1)
- dj = (j2-j1+1) / (j2X-j1X+1)
- select case ( nuv )
-
- case ( 'n' )
- ! loop over (selection of) coarse cells in coarse source grid lliX:
- do jX = j1X, j2X
- do iX = i1X, i2X
-
- ! loop over cells in fine target grid covered by coarse cell:
- do j = j1 + (jX-j1X)*dj, j1 + (jX-j1X)*dj + dj-1
- do i = i1 + (iX-i1X)*di, i1 + (iX-i1X)*di + di-1
-
- ! copy or take fraction of coarse value:
- select case ( combkey )
- case ( 'sum' )
- ! coarse cell is sum of finer; take fraction:
- !ll(i,j) = llX(iX,jX) / real(di*dj)
- ! coarse cell is sum of finer; take area fractions:
- ll(i,j) = llX(iX,jX) * lli%area(j) / lliX%area(jX)
- case ( 'aver', 'area-aver', 'weight','mass-aver')
- ! coarse cell is aver of finer; take all the same:
- ll(i,j) = llX(iX,jX)
- case default
- write (*,'("ERROR - unsupported combkey `",a,"`")') combkey
- write (*,'("ERROR - nuv : `",a,"`")') nuv
- write (*,'("ERROR - action : `",a,"`")') action
- write (*,'("ERROR in ",a)') rname; status=1; return
- end select
- end do
- end do
-
- end do
- end do
-
- case ( 'u' )
-
- !
- ! coarse iaX ibX
- ! | | |
- ! | | | | | | | | |
- ! fine ia ib
- ! i
- !
- ! loop over (selection of) u boundaries of target grid lli:
- do j = j1, j2
- do i = i1, i2+1
-
- ! coarse cell in j direction:
- jX = j1X + floor(real(j-j1)/real(dj))
- ! left and right bound in fine grid that are covered by coarse bound:
- ia = i1 + floor(real(i-i1)/real(di))*di
- ib = i1 + ceiling(real(i-i1)/real(di))*di
-
- ! corresponding coarse boundaries:
- iaX = i1X + (ia-i1)/di
- ibX = i1X + (ib-i1)/di
-
- ! interpolation in i direction of surounding lliX boundaries:
- if ( iaX == ibX ) then
- llX_ab = llX(iaX,jX)
- else
- llX_ab = llX(iaX,jX) * real(ib-i)/real(di) + &
- llX(ibX,jX) * real(i-ia)/real(di)
- end if
- ! copy or take fraction of coarse value:
- select case ( combkey )
- case ( 'sum' )
- ! coarse cell is sum of finer; take fraction:
- ll(i,j) = llX_ab / real(dj)
- case ( 'aver', 'area-aver','mass-aver' )
- ! coarse cell is aver of finer; take all the same:
- ll(i,j) = llX_ab
- case default
- write (*,'("ERROR - unsupported combkey `",a,"`")') combkey
- write (*,'("ERROR - nuv : `",a,"`")') nuv
- write (*,'("ERROR - action : `",a,"`")') action
- write (*,'("ERROR in ",a)') rname; status=1; return
- end select
- end do
- end do
-
- case ( 'v' )
-
- !
- ! coarse fine
- !
- ! --
- ! jaX -- -- ja
- ! -- j
- ! --
- ! jbX -- -- jb
- ! --
- !
- ! loop over (selection of) u boundaries of target grid lli:
- do i = i1, i2
- do j = j1, j2+1
-
- ! coarse cell in i direction:
- iX = i1X + floor(real(i-i1)/real(di))
- ! left and right bound in fine grid that are covered by coarse bound:
- ja = j1 + floor(real(j-j1)/real(dj))*dj
- jb = j1 + ceiling(real(j-j1)/real(dj))*dj
-
- ! corresponding coarse boundaries:
- jaX = j1X + (ja-j1)/dj
- jbX = j1X + (jb-j1)/dj
-
- ! interpolation in j direction of surounding lliX boundaries:
- if ( jaX == jbX ) then
- llX_ab = llX(iX,jaX)
- else
- llX_ab = llX(iX,jaX) * real(jb-j)/real(dj) + &
- llX(iX,jbX) * real(j-ja)/real(dj)
- end if
- ! copy or take fraction of coarse value:
- select case ( combkey )
- case ( 'sum' )
- ! coarse cell is sum of finer; take fraction:
- ll(i,j) = llX_ab / real(di)
- case ( 'aver', 'area-aver','mass-aver' )
- ! coarse cell is aver of finer; take all the same:
- ll(i,j) = llX_ab
- case default
- write (*,'("ERROR - unsupported combkey `",a,"`")') combkey
- write (*,'("ERROR - nuv : `",a,"`")') nuv
- write (*,'("ERROR - action : `",a,"`")') action
- write (*,'("ERROR in ",a)') rname; status=1; return
- end select
- end do
- end do
-
- case default
- write (*,'("ERROR - unsupported nuv `",a,"`")') nuv
- write (*,'("ERROR - action : `",a,"`")') action
- write (*,'("ERROR in ",a)') rname; status=1; return
- end select
- !
- ! * error ....
- !
-
- case default
- write (*,'("ERROR - unsupported action `",a,"`")') action
- write (*,'("ERROR - lliX lon : ",f7.2,f6.2,i4)') lliX%lon_deg(1), lliX%dlon_deg, lliX%nlon
- write (*,'("ERROR - lat : ",f7.2,f6.2,i4)') lliX%lat_deg(1), lliX%dlat_deg, lliX%nlat
- write (*,'("ERROR - lli lon : ",f7.2,f6.2,i4)') lli%lon_deg(1), lli%dlon_deg, lli%nlon
- write (*,'("ERROR - lat : ",f7.2,f6.2,i4)') lli%lat_deg(1), lli%dlat_deg, lli%nlat
- write (*,'("ERROR in ",a)') rname; status=1; return
- end select
- ! ok
- if ( only_part_of_ll ) then
- status = -1
- else
- status = 0
- end if
- end subroutine FillGrid
-
- !
- ! determine how lliX is related to lli:
- ! i1, i2, j1, j2 : cell ranges in lli covered by cells of lliX
- ! i1X, i2X, j1X, j2X : cell ranges in lliX covering cells of lliX
- ! action : 'copy', 'combine'
- !
- subroutine Relate( lli , i1 , i2 , j1 , j2 , &
- lliX, i1X, i2X, j1X, j2X, &
- action, status )
- ! --- in/out --------------------------------
- type(TllGridInfo), intent(in) :: lli
- integer, intent(out) :: i1, i2, j1, j2
- type(TllGridInfo), intent(in) :: lliX
- integer, intent(out) :: i1X, i2X, j1X, j2X
- character(len=10), intent(out) :: action
- integer, intent(out) :: status
- ! --- const ---------------------------------
- character(len=*), parameter :: name = mname//'/Relate'
-
- ! --- local ---------------------------------
- real :: dlon, dlonX
- integer :: nlon, nlonX
- real :: dlat, dlatX
- integer :: nlat, nlatX
- real :: west, westX
- real :: east, eastX
- real :: south, southX
- real :: north, northX
- real :: r, rX
-
- logical :: only_part_of_ll
-
- ! --- begin ---------------------------------
-
- ! shorthands ...
- dlon = lli%dlon_deg ; dlonX = lliX%dlon_deg
- nlon = lli%nlon ; nlonX = lliX%nlon
- dlat = lli%dlat_deg ; dlatX = lliX%dlat_deg
- nlat = lli%nlat ; nlatX = lliX%nlat
- west = lli%blon_deg(0) ; westX = lliX%blon_deg(0)
- east = lli%blon_deg(nlon) ; eastX = lliX%blon_deg(nlonX)
- south = lli%blat_deg(0) ; southX = lliX%blat_deg(0)
- north = lli%blat_deg(nlat) ; northX = lliX%blat_deg(nlatX)
- ! same spacing ?
- if ( (dlonX == dlon) .and. (dlatX == dlat) ) then
- ! cells from lliX should be copied to lli
- action = 'copy'
- else if ( (dlonX <= dlon) .and. (dlatX <= dlat) ) then
- ! cells from lliX should be combined:
- action = 'combine'
- else if ( (dlonX >= dlon) .and. (dlatX >= dlat) ) then
- ! distribute coarse cells of lliX over fine cells of lli:
- action = 'distribute'
- else
- write (*,'("ERROR - do not know how to relate grids:")')
- write (*,'("ERROR - lli dlon,dlat :",2f10.4)') dlon , dlat
- write (*,'("ERROR - lliX dlon,dlat :",2f10.4)') dlonX, dlatX
- write (*,'("ERROR in ",a)') name; status=1; return
- end if
-
- ! assume by default that all is filled
- only_part_of_ll = .false.
-
- ! west boundary
- r = abs(west - westX) / dlon
- rX = abs(west - westX) / dlonX
- ! relate lliX to lli:
- if ( (westX <= west) .and. (rX == nint(rX)) .and. (nint(rX)+1 <= nlonX) ) then
- ! at this side, all of lli is covered by lliX:
- i1 = 1
- i1X = nint(rX) + 1
- else if ( (westX > west) .and. (r == nint(r)) .and. (nint(r)+1 <= nlon) ) then
- ! at this side, only a part of lli is covered by lliX:
- i1 = nint(r) + 1
- i1X = 1
- only_part_of_ll = .true.
- else
- write (*,'("ERROR - do not know how to relate west bounds:")')
- write (*,'("ERROR - lli bound, spacing, number :",2f10.4,i6)') west , dlon , nlon
- write (*,'("ERROR - lliX bound, spacing, number :",2f10.4,i6)') westX, dlonX, nlonX
- write (*,'("ERROR in ",a)') name; status=1; return
- end if
-
- ! east boundary
- r = abs(east - eastX) / dlon
- rX = abs(east - eastX) / dlonX
- ! relate lliX to lli:
- if ( (eastX >= east) .and. (rX == nint(rX)) .and. (nlonX-nint(rX) >= 1) ) then
- ! at this side, all of lli is covered by lliX:
- i2 = nlon
- i2X = nlonX - nint(rX)
- else if ( (eastX < east) .and. (r == nint(r)) .and. (nlon-nint(r) >= 1) ) then
- ! at this side, only a part of lli is covered by lliX:
- i2 = nlon - nint(r)
- i2X = nlonX
- only_part_of_ll = .true.
- else
- write (*,'("ERROR - do not know how to relate east bounds:")')
- write (*,'("ERROR - lli bound, spacing, number :",2f10.4,i6)') east , dlon , nlon
- write (*,'("ERROR - lliX bound, spacing, number :",2f10.4,i6)') eastX, dlonX, nlonX
- write (*,'("ERROR in ",a)') name; status=1; return
- end if
-
- ! south boundary
- r = abs(south - southX) / dlat
- rX = abs(south - southX) / dlatX
- ! relate lliX to lli:
- if ( (southX <= south) .and. (rX == nint(rX)) .and. (nint(rX)+1 <= nlatX) ) then
- ! at this side, all of lli is covered by lliX:
- j1 = 1
- j1X = nint(rX) + 1
- else if ( (southX > south) .and. (r == nint(r)) .and. (nint(r)+1 <= nlat) ) then
- ! at this side, only a part of lli is covered by lliX:
- j1 = nint(r) + 1
- j1X = 1
- only_part_of_ll = .true.
- else
- write (*,'("ERROR - do not know how to relate south bounds:")')
- write (*,'("ERROR - lli bound, spacing, number :",2f10.4,i6)') south , dlat , nlat
- write (*,'("ERROR - lliX bound, spacing, number :",2f10.4,i6)') southX, dlatX, nlatX
- write (*,'("ERROR in ",a)') name; status=1; return
- end if
-
- ! north boundary
- r = abs(north - northX) / dlat
- rX = abs(north - northX) / dlatX
- ! relate lliX to lli:
- if ( (northX >= north) .and. (rX == nint(rX)) .and. (nlatX-nint(rX) >= 1) ) then
- ! at this side, all of lli is covered by lliX:
- j2 = nlat
- j2X = nlatX - nint(rX)
- else if ( (northX < north) .and. (r == nint(r)) .and. (nlat-nint(r) >= 1) ) then
- ! at this side, only a part of lli is covered by lliX:
- j2 = nlat - nint(r)
- j2X = nlatX
- only_part_of_ll = .true.
- else
- write (*,'("ERROR - do not know how to relate north bounds:")')
- write (*,'("ERROR - lli bound, spacing, number :",2f10.4,i6)') north , dlat , nlat
- write (*,'("ERROR - lliX bound, spacing, number :",2f10.4,i6)') northX, dlatX, nlatX
- write (*,'("ERROR in ",a)') name; status=1; return
- end if
- ! ok
- if ( only_part_of_ll ) then
- status = -1
- else
- status = 0
- end if
- end subroutine Relate
-
-
- ! ========================================================
- ! ===
- ! === grid data operations
- ! ===
- ! ========================================================
-
- ! subroutine Divergence_ll( lli, u, v, div )
- !
- ! ! --- in/out -------------------------------
- !
- ! type(TllGridInfo), intent(in) :: lli
- ! type(TllGrid), intent(in) :: u, v
- ! type(TllGrid), intent(out) :: div
- !
- ! ! --- local --------------------------------
- !
- ! integer :: i, j
- ! integer :: im1, ip1
- !
- ! real :: wuim1(0:lme)
- ! real :: wuip1(0:lme)
- ! real :: wvjm1(0:lme)
- ! real :: wvjp1(0:lme)
- !
- ! ! --- begin ---------------------------------
- !
- ! if ( lli%i1 == lli%im ) then
- ! stop 'at least grid of 2 columns'
- ! end if
- !
- ! if ( lli%j1 == lli%jm ) then
- ! stop 'at least grid of 2 rows'
- ! end if
- !
- ! stop 'circular ?'
- !
- ! do j = lli%j1, lli%jm
- !
- ! do i = lli%i1, lli%im
- !
- ! if ( j==lli%j1 .or. j==lli%jm ) then
- !
- ! div%d(i,j) = 0.0
- !
- ! else
- !
- ! if ( i==lli%i1 ) then
- ! im1 = lli%im
- ! ip1 = i + 1
- ! else if ( i==lli%im ) then
- ! im1 = lli%im
- ! ip1 = i + 1
- ! else
- ! im1 = i - 1
- ! ip1 = i + 1
- ! end if
- !
- ! wuim1 = wu(im1,j,:)
- ! wuip1 = wu(ip1,j,:)
- !
- ! wvjm1 = wv(i,j-1,:)
- ! wvjp1 = wv(i,j+1,:)
- !
- ! dphi2 = eclon(2) - eclon(0)
- ! dlat2 = eclat(j+1) - eclat(j-1)
- !
- ! qam(i,j,:) = ( wuip1 - wuim1 ) / (ceclat(j)*dphi2*ae) + &
- ! ( wvjp1 - wvjm1 ) / (dlat2*ae)
- !
- ! end if
- !
- ! end do
- ! end do
- !
- ! end subroutine divergence_uv
-
-
- ! ==================================================================
- !
- ! mass flux balance
- !
- ! ==================================================================
-
- !
- ! (copied from TMPP;
- ! code should be cleaned ..)
- !
- !
- ! lli : fine grid definition
- ! cgi : coarse grid covering fine grid.
- !
- ! On input, fine grid fields pu/pv/pw should match with
- ! their coarse grid equivalents.
- !
- subroutine BalanceMassFluxes_sm( lli, lme, pu, pv, ww, dmdt, cgi, b_ec, status )
-
- ! --- in/out -------------------------------------
- type(TllGridInfo), intent(in) :: lli
- integer, intent(in) :: lme
- real,intent(inout) :: pu(0:lli%im,lli%jm,lme)
- real,intent(inout) :: pv(lli%im,0:lli%jm,lme)
- real, intent(in) :: ww(lli%im,lli%jm,lme)
- real, intent(in) :: dmdt(lli%im,lli%jm)
- type(TllGridInfo), intent(in) :: cgi
- real, intent(in) :: b_ec(0:lme)
- integer, intent(out) :: status
-
- ! --- const -----------------------------
-
- character(len=*), parameter :: rname = mname//'/BalanceMassFluxes_sm'
-
- ! --- local ----------------------------------
- integer :: refine_i, refine_j
- integer :: cg_i1, cg_i2, cg_j1, cg_j2
- integer :: ci, cj
- integer :: fg_i1, fg_i2, fg_j1, fg_j2
- integer :: fi, fj
- integer :: i,j,l
- integer :: ifail
- real :: massc(lli%im,lli%jm)
- real :: dm(lli%im,lli%jm)
- real, allocatable :: cqu(:,:), cqv(:,:)
- ! --- begin -----------------------------------
-
- !print *, '<pascha_zoom>'
-
- !print *, ' correct layers 1 to',lme,' with (eta dp/deta) and dps/dt '
- ! determine refinement
- call GetRefinement( cgi, lli, refine_i, refine_j, cg_i1, cg_i2, cg_j1, cg_j2, status )
- if (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; return; end if
-
- ! correction fluxes:
- allocate( cqu(0:refine_i,refine_j) )
- allocate( cqv(refine_i,0:refine_j) )
-
- do l = 1, lme
- ! total horizontal massconvergence
- massc = 0.0
- massc = massc + sum(pu(1:lli%im ,:,1:l),3) ! leaving through east
- massc = massc - sum(pu(0:lli%im-1,:,1:l),3) ! enter through west
- massc = massc + sum(pv(:,1:lli%jm ,1:l),3) ! leaving through north
- massc = massc - sum(pv(:,0:lli%jm-1,1:l),3) ! enter through south
- ! set difference
- dm = (-massc) - ww(:,:,l) - b_ec(l)*dmdt
- ! loop over cells in coarse grid covering fine grid:
- do cj = cg_j1, cg_j2
- do ci = cg_i1, cg_i2
- fg_i1 = (ci-cg_i1)*refine_i + 1 ; fg_i2 = fg_i1-1 + refine_i
- fg_j1 = (cj-cg_j1)*refine_j + 1 ; fg_j2 = fg_j1-1 + refine_j
- ! do FFT
- call SolvePoissonEq_zoom( refine_i, refine_j, &
- dm(fg_i1:fg_i2,fg_j1:fg_j2), &
- cqu, cqv, status )
- if (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; return; end if
- ! add corrections to original flux:
- pu(fg_i1-1:fg_i2,fg_j1 :fg_j2,l) = pu(fg_i1-1:fg_i2,fg_j1 :fg_j2,l) + cqu
- pv(fg_i1 :fg_i2,fg_j1-1:fg_j2,l) = pv(fg_i1 :fg_i2,fg_j1-1:fg_j2,l) + cqv
- end do
- end do
- end do
-
- ! done
- deallocate( cqu )
- deallocate( cqv )
-
- ! ok
- status = 0
-
- end subroutine BalanceMassFluxes_sm
-
-
- ! ***
-
- !
- ! pw : vertical flux in direction of increasing level
- !
- subroutine BalanceMassFluxes_m( lli, pu, pv, pw, dm_dt, cgi, dt_sec, status )
- ! --- in/out -------------------------------------
- type(TllGridInfo), intent(in) :: lli
- real,intent(inout) :: pu(:,:,:)
- real,intent(inout) :: pv(:,:,:)
- real, intent(in) :: pw(:,:,:)
- real, intent(in) :: dm_dt(:,:,:)
- type(TllGridInfo), intent(in) :: cgi
- real, intent(in) :: dt_sec
- integer, intent(out) :: status
-
- ! --- const -----------------------------
-
- character(len=*), parameter :: rname = mname//'/BalanceMassFluxes_m'
-
- ! --- local ----------------------------------
- integer :: refine_i, refine_j
- integer :: cg_i1, cg_i2, cg_j1, cg_j2
- integer :: ci, cj
- integer :: fg_i1, fg_i2, fg_j1, fg_j2
- integer :: fi, fj
- integer :: i,j,l
- integer :: nlev
- real :: massc(lli%im,lli%jm)
- real :: dm(lli%im,lli%jm)
- real, allocatable :: cqu(:,:), cqv(:,:)
- ! --- begin -----------------------------------
- ! number of levels
- nlev = size(pu,3)
- ! check ...
- if ( (size(pu ,1) /= lli%nlon+1) .or. (size(pu ,2) /= lli%nlat ) .or. &
- (size(pv ,1) /= lli%nlon ) .or. (size(pv ,2) /= lli%nlat+1) .or. &
- (size(pw ,1) /= lli%nlon ) .or. (size(pw ,2) /= lli%nlat ) .or. &
- (size(dm_dt,1) /= lli%nlon ) .or. (size(dm_dt,2) /= lli%nlat ) .or. &
- (size(pv ,3) /= nlev) .or. (size(pw,3) /= nlev+1) .or. &
- (size(dm_dt,3) /= nlev) ) then
- write (*,'("ERROR - input arrays do not match with grid definition or levels:")')
- write (*,'("ERROR - lli : ",i3," x ",i3)') lli%nlon, lli%nlat
- write (*,'("ERROR - pu : ",i3," x ",i3," x ",i3)') shape(pu)
- write (*,'("ERROR - pv : ",i3," x ",i3," x ",i3)') shape(pv)
- write (*,'("ERROR - pw : ",i3," x ",i3," x ",i3)') shape(pw)
- write (*,'("ERROR - dm_dt : ",i3," x ",i3," x ",i3)') shape(dm_dt)
- write (*,'("ERROR in ",a)') rname; status=1; return
- end if
- ! determine refinement:
- ! refine_i, refine_j
- ! cg_i1, cg_i2, cg_j1, cg_j2 : cells in coarse grid cgi covering fine grid lli:
- call GetRefinement( cgi, lli, refine_i, refine_j, cg_i1, cg_i2, cg_j1, cg_j2, status )
- if (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; return; end if
-
- !
- ! determine how lliX (finer?) is related to lli:
- ! i1, i2, j1, j2 : cell ranges in lli covered by cells of lliX
- ! i1X, i2X, j1X, j2X : cell ranges in lliX covering cells of lliX
- ! action : 'copy', 'combine' to map lliX to lli
- !
- !
- !subroutine Relate( lli , i1 , i2 , j1 , j2 , &
- ! lliX, i1X, i2X, j1X, j2X, &
- ! action, status )
- ! correction fluxes:
- allocate( cqu(0:refine_i,refine_j) )
- allocate( cqv(refine_i,0:refine_j) )
- ! loop over levels
- do l = 1, nlev
- ! total horizontal mass convergence
- massc = 0.0
- massc = massc + pu(1:lli%im ,:,l) ! enter through west
- massc = massc - pu(2:lli%im+1,:,l) ! leaving through east
- massc = massc + pv(:,1:lli%jm ,l) ! enter through south
- massc = massc - pv(:,2:lli%jm+1,l) ! leaving through north
- massc = massc + pw(:,:,l ) ! enter through bottom
- massc = massc - pw(:,:,l+1) ! leaving through top
- ! set difference
- dm = massc - dm_dt(:,:,l)
- ! loop over cells in coarse grid covering fine grid:
- do cj = cg_j1, cg_j2
- do ci = cg_i1, cg_i2
- fg_i1 = (ci-cg_i1)*refine_i + 1 ; fg_i2 = fg_i1-1 + refine_i
- fg_j1 = (cj-cg_j1)*refine_j + 1 ; fg_j2 = fg_j1-1 + refine_j
- ! do FFT
- call SolvePoissonEq_zoom( refine_i, refine_j, &
- dm(fg_i1:fg_i2,fg_j1:fg_j2), &
- cqu, cqv, status )
- if (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; return; end if
- ! add corrections to original flux:
- pu(fg_i1:fg_i2+1,fg_j1:fg_j2 ,l) = pu(fg_i1:fg_i2+1,fg_j1:fg_j2 ,l) + cqu
- pv(fg_i1:fg_i2 ,fg_j1:fg_j2+1,l) = pv(fg_i1:fg_i2 ,fg_j1:fg_j2+1,l) + cqv
-
- end do ! loop over coarse cells i
- end do ! loop over coarse cells j
- end do ! loop over levels
-
- ! done
- deallocate( cqu )
- deallocate( cqv )
-
- ! ok
- status = 0
-
- end subroutine BalanceMassFluxes_m
-
- ! ***
-
- subroutine CheckMassBalance( lli, pu, pv, sp1, sp2, dt, &
- rms_max, mav_max, status )
-
- use Binas, only : grav
-
- ! --- in/out -------------------------------------
- type(TllGridInfo), intent(in) :: lli
- real,intent(inout) :: pu(:,:,:)
- real,intent(inout) :: pv(:,:,:)
- real, intent(in) :: sp1(:,:)
- real, intent(in) :: sp2(:,:)
- real, intent(in) :: dt
- real, intent(in) :: rms_max
- real, intent(in) :: mav_max
- integer, intent(out) :: status
-
- ! --- const -----------------------------
-
- character(len=*), parameter :: rname = mname//'/CheckMassBalance'
-
- ! --- local ----------------------------------
- real :: dmuv_dt(lli%im,lli%jm)
- real :: dmsp_dt(lli%im,lli%jm)
- real :: dm_dt_diff(lli%im,lli%jm)
- real :: rms_diff
- real :: mav_diff
- ! --- begin -----------------------------------
- ! check ...
- if ( (size(pu ,1) /= lli%nlon+1) .or. (size(pu ,2) /= lli%nlat ) .or. &
- (size(pv ,1) /= lli%nlon ) .or. (size(pv ,2) /= lli%nlat+1) .or. &
- (size(sp1,1) /= lli%nlon ) .or. (size(sp1,2) /= lli%nlat ) .or. &
- (size(sp2,1) /= lli%nlon ) .or. (size(sp2,2) /= lli%nlat ) .or. &
- (size(pu,3) /= size(pv,3)) ) then
- write (*,'("ERROR - input arrays do not match with grid definition or levels:")')
- write (*,'("ERROR - lli : ",i3," x ",i3)') lli%nlon, lli%nlat
- write (*,'("ERROR - pu : ",i3," x ",i3," x ",i3)') shape(pu)
- write (*,'("ERROR - pv : ",i3," x ",i3," x ",i3)') shape(pv)
- write (*,'("ERROR - sp1 : ",i3," x ",i3)') shape(sp1)
- write (*,'("ERROR - sp2 : ",i3," x ",i3)') shape(sp2)
- write (*,'("ERROR in ",a)') rname; status=1; return
- end if
- ! total horizontal mass convergence (kg/s)
- dmuv_dt = 0.0
- dmuv_dt = dmuv_dt + sum(pu(1:lli%nlon ,:,:),3) ! enter through west
- dmuv_dt = dmuv_dt - sum(pu(2:lli%nlon+1,:,:),3) ! leaving through east
- dmuv_dt = dmuv_dt + sum(pv(:,1:lli%nlat ,:),3) ! enter through south
- dmuv_dt = dmuv_dt - sum(pv(:,2:lli%nlat+1,:),3) ! leaving through north
- ! convert from kg/s to Pa/s
- dmuv_dt = dmuv_dt * grav ! Pa/s/m2
- call AreaOper( lli, dmuv_dt, '/', 'm2', status ) ! Pa/s
- if (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; return; end if
-
- ! same according to surface pressures:
- dmsp_dt = ( sp2 - sp1 ) / dt
-
- ! imbalance relative to average surface pressure:
- dm_dt_diff = ( dmuv_dt - dmsp_dt ) / ( (sp1+sp2)/2.0 )
-
- ! statistics
- rms_diff = sqrt( sum(dm_dt_diff**2) / size(dm_dt_diff) )
- mav_diff = maxval( abs( dm_dt_diff ) )
-
- ! check ...
- if ( (rms_diff > rms_max) .or. (mav_diff > mav_max) ) then
- write (*,'("ERROR - found relative surface pressure imbalance;")')
- write (*,'("ERROR - difference [dmuv/dt]/[(sp1+sp2)/2] and [dsp/dt]/[(sp1+sp2)/2] :")')
- write (*,'("ERROR - rms : ",es12.4," (",es12.4,")")') rms_diff, rms_max
- write (*,'("ERROR - max abs. : ",es12.4," (",es12.4,")")') mav_diff, mav_max
- write (*,'("ERROR in ",a)') rname; status=1; return
- end if
- ! ok
- status = 0
-
- end subroutine CheckMassBalance
-
- ! ***
-
-
- ! Return flux field (u,v) such that netto flux in a cell matches m:
- !
- ! du dv
- ! -- + -- = M
- ! dx dy
- !
- ! Fluxes at grid boundaries set to zero.
- !
- ! IN/OUT
- ! input:
- ! m(1:im,1:jm) : target netto flux
- ! output:
- ! u(0:im,1:jm) : u flux through cell boundaries
- ! v(1:im,0:jm) : v vlux
- !
- !
- ! +.......+.......+.......+.......+
- ! : : : : :
- ! : o : o : o : o :
- ! : : : : :
- ! jm +.......+---v---+---v---+---v---+---v---+.......+
- ! : | | | | | :
- ! : o u * u * u * u * u o :
- ! : | | | | | :
- ! : +.......+---v---+---v---+---v---+---v---+.......+
- ! : | | | | | :
- ! : o u * u * u * u * u o :
- ! : | | | | | :
- ! 1 +.......+---v---+---v---+---v---+---v---+.......+
- ! : | | | | | :
- ! : o u * u * u * u * u o :
- ! : | | | | | :
- ! j = 0 +.......+---v---+---v---+---v---+---v---+.......+
- ! : : : : :
- ! : o : o : o : o :
- ! : : : : :
- ! +.......+.......+.......+.......+
- !
- ! i= 0 1 2 .. im
- !
- ! * = m, Psi
- ! o = Psi (periodic)
- !
- !
- ! ALGORITHM
- !
- ! 1. Solve the Poisson equation:
- !
- ! d^2 d^2
- ! ( ---- + ---- ) Psi = M
- ! dx^2 dy^2
- !
- ! and return F = (u,v) = ( dPsi/dx, dPsi/dy )
- ! The differential equations is actually a difference equation:
- ! sum of fluxes is prescribed, we assume that the fluxes are differences
- ! of a discrete potential Psi.
- !
- ! 2. The sollution is periodic in i and j.
- ! To obtain the zero boundary conditions, substract the boundaries
- ! from each row and column.
- !
-
- subroutine SolvePoissonEq_zoom( im, jm, m, u, v, status )
- use Binas, only : pi
- use grid_singleton, only : fftkind, fft
- ! --- in/out -------------------------------------
- integer, intent(in) :: im, jm
- real, intent(in) :: m(im,jm)
- real, intent(out) :: u(0:im,jm)
- real, intent(out) :: v(im,0:jm)
- integer, intent(out) :: status
-
- ! --- const -----------------------------
-
- character(len=*), parameter :: rname = mname//'/SolvePoissonEq_zoom'
-
- ! --- local ----------------------------------
- integer :: i,j
- integer :: js, je
- real :: fac(im,jm)
- complex(fftkind) :: A(im,jm)
- complex(fftkind) :: X(im,jm)
- real :: psi(im,jm)
- real :: row(im), col(jm)
- ! --- begin -----------------------------------
- ! precalculate cos/sin array fac
- !$OMP PARALLEL &
- !$OMP default (none) &
- !$OMP shared (im, jm, fac) &
- !$OMP private (i, j, js, je)
- do j = 1, jm
- do i=1,im
- fac(i,j) = 2.0*(cos(2*pi*(i-1)/im)+cos(2*pi*(j-1)/jm)-2.)
- end do
- end do
- !$OMP END PARALLEL
- fac(1,1) = 1.0 !to avoid division by zero
- ! do FFT
- A = cmplx( m, 0.0 )
- X = fft( A, stat=status )
- if ( status /= 0 ) then
- write (*,'("ERROR - from first fft; stat = ",i6)') status
- write (*,'("ERROR in ",a)') rname; status=1; return
- end if
- ! compute fourier coefficients of correction potential
- A = cmplx( real(X)/fac, -aimag(X)/fac )
- A(1,1) = (0.0,0.0)
- ! get correction potential
- X = fft( A, stat=status )
- if ( status /= 0 ) then
- write (*,'("ERROR - from second fft; stat = ",i6)') status
- write (*,'("ERROR in ",a)') rname; status=1; return
- end if
- psi = real( X )
- ! correction flux in lon direction:
- do j = 1, jm
- !u(:,j) = (/ psi(1,j)-psi(im,j), psi(2:im,j)-psi(1:im-1,j), psi(1,j)-psi(im,j) /)
- u(0 ,j) = psi(1 ,j)-psi(im ,j)
- u(1:im-1,j) = psi(2:im,j)-psi(1:im-1,j)
- u(im ,j) = psi(1 ,j)-psi(im ,j)
- end do
- ! substract flux over east/west boundary:
- col = u(0,:)
- do i = 0, im
- u(i,:) = u(i,:) - col
- end do
- ! correction flux in lat direction;
- do i = 1, im
- !v(i,:) = (/ psi(i,1)-psi(i,jm), psi(i,2:jm)-psi(i,1:jm-1), psi(i,1)-psi(i,jm) /)
- v(i,0 ) = psi(i,1 ) - psi(i,jm )
- v(i,1:jm-1) = psi(i,2:jm) - psi(i,1:jm-1)
- v(i,jm ) = psi(i,1 ) - psi(i,jm )
- end do
- ! substract flux over north/west boundary:
- row = v(:,0)
- do j = 0, jm
- v(:,j) = v(:,j) - row
- end do
-
- ! Correction for flux over poles.
- ! The equation is solved with a 2D Fourrier transform, which
- ! assumes that the grid is part of an into infinity periodical
- ! extended grid; the solution is thus periodic from east into west
- ! but also from southpole into north pole (donut configuration).
- ! To obtain a zero flux over the poles, first the orginal flux over
- ! the poles is computed:
- ! vpole = psi(:,1) - psi(:,jm)
- ! This flux is substracted from each of the lat fluxes:
- ! v(:,j) := v(:,j) - vpole
- ! In this way, the netto mass convergence is maintained,
- ! but (u,v) is not a potential flow anymore!
- ! We will not bother about this small 'error', since we are dealing
- ! with a correction to a much larger error.
- end subroutine SolvePoissonEq_zoom
-
- ! =================================================================
- ! ===
- ! === equivalent latitudes
- ! ===
- ! =================================================================
- !
- ! eqvlatb(1:jm+1) : bounds of eqv.lat. bins (/-90.0,..,90.0/)
- !
- ! inds(im,jm) : cell (i,j) is in eqvlatb( inds(i,j), inds(i,j)+1 )
- !
- subroutine llgrid_EquivLat( lli, ff, eqvlatb, inds )
-
- use Binas, only : deg2rad, pi
- use Num, only : Interval, Interp_Lin
- use grid_tools, only : ll_area
-
- ! --- in/out -----------------------------
-
- type(TllGridInfo), intent(in) :: lli
- real, intent(in) :: ff(:,:)
- real, intent(out) :: eqvlatb(:)
- integer, intent(out) :: inds(:,:)
-
- integer :: status
-
- ! --- const ------------------------------
-
- character(len=*), parameter :: rname = mname//'/llgrid_EquivLat'
-
- ! --- local ------------------------------
-
- real :: fmin, fmax, f
-
- integer :: nbins
- real, allocatable :: bin_inds(:)
- real, allocatable :: bin_fval(:)
- real, allocatable :: bin_area(:)
- real, allocatable :: bin_areacum(:)
-
- real :: indf
-
- real :: area_tot
- real :: area_glob, area_south
-
- integer :: ilast
- integer :: i, j
-
- ! --- begin ------------------------------
-
- !if ( (maxval(lli%blon_deg)-minval(lli%blon_deg) < 360.0) .or. &
- ! (maxval(lli%blat_deg)-minval(lli%blat_deg) < 180.0) ) then
- ! write (*,'("ERROR - equivalent lats for global grids only yet")')
- ! write (*,'("ERROR - lons : ",f8.2,"-",f8.2)') minval(lli%blon_deg), maxval(lli%blon_deg)
- ! write (*,'("ERROR - lats : ",f8.2,"-",f8.2)') minval(lli%blat_deg), maxval(lli%blat_deg)
- ! write (*,'("ERROR in ",a)') name; stop
- !end if
-
- call Check( lli, 'n', ff, status )
- if (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; stop; end if
- call Check( lli, 'n', inds, status )
- if (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; stop; end if
-
- ! for the moment ...
- nbins = lli%jm
- if ( size(eqvlatb) /= nbins+1 ) then
- write (*,'("ERROR - routine returns jm+1 bounds of jm eqv.lat bins:")')
- write (*,'("ERROR - size(eqvlatb) : ",i6)') size(eqvlatb)
- write (*,'("ERROR - nbins+1 (=jm+1) : ",i6)') nbins+1
- write (*,'("ERROR in ",a)') rname; stop
- end if
-
- ! number of bins, field values, total area
- allocate( bin_inds(nbins) )
- allocate( bin_fval(nbins) )
- allocate( bin_area(nbins) )
- allocate( bin_areacum(nbins) )
-
- ! linear ax of values in field:
- fmin = minval(ff)
- fmax = maxval(ff)
- do j = 1, nbins
- bin_inds(j) = j*1.0
- bin_fval(j) = fmin*(nbins-j)/real(nbins-1) + fmax*(j-1)/real(nbins-1)
- end do
-
- ! add weights for each cell to corresponding bin:
- ilast = 0
- bin_area = 0.0
- do i = 1, lli%im
- do j = 1, lli%jm
- ! search index in bins:
- call Interp_Lin( bin_fval, bin_inds, ff(i,j), indf, ilast, status )
- if (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; stop; end if
- ! store index:
- inds(i,j) = nint(indf)
- ! add area to corresponding bin:
- bin_area(inds(i,j)) = bin_area(inds(i,j)) + lli%area(j) ! rad^2
- end do
- end do
-
- ! total area:
- area_tot = lli%im * sum(lli%area)
-
- ! check ...
- if ( (sum(bin_area)-area_tot)/area_tot > 0.01 ) then
- write (*,'("ERROR - total area and areas in bins do not match:")')
- write (*,'("ERROR - total area : ",es12.2)') area_tot
- write (*,'("ERROR - collected area : ",es12.2)') sum(bin_area)
- write (*,'("ERROR in ",a)') rname; stop
- end if
-
- ! cumulative bin area
- bin_areacum(1) = bin_area(1)
- do j = 2, nbins
- bin_areacum(j) = bin_areacum(j-1) + bin_area(j)
- end do
- !
- ! blat(nlat) ------------------------------- 90
- ! | | \
- ! |----------| |
- ! |//////////| |
- ! |----------| | glob
- ! | | \ |
- ! | | | south |
- ! | | / /
- ! blat(nlat) ------------------------------- -90
- ! blon(0) blon(nlon)
-
- ! global area and southern part (all in rad!)
- area_glob = ll_area( 0.0, lli%blon(lli%nlon)-lli%blon(0), -0.5*pi, 0.5*pi ) ! rad^2
- area_south = ll_area( 0.0, lli%blon(lli%nlon)-lli%blon(0), -0.5*pi, lli%blat(0) ) ! rad^2
- ! convert cumulative ax to (part of) [-1,1]
- bin_areacum = (area_south + bin_areacum ) / area_glob ! [ 0,1]
- bin_areacum = -1.0 + 2.0*bin_areacum ! [-1,1]
- where ( bin_areacum < -1.0 ) bin_areacum = -1.0
- where ( bin_areacum > 1.0 ) bin_areacum = 1.0
-
- ! convert to latitude:
- eqvlatb(1) = lli%blat_deg(0) ! deg
- do j = 1, nbins
- eqvlatb(j+1) = asin( bin_areacum(j) ) / deg2rad
- end do
- !print *, 0, area_south/area_glob, eqvlatb(1)
- !do j=1,nbins
- ! print *, j, bin_areacum(j), eqvlatb(j+1)
- !end do
- !stop
-
- ! done
- deallocate( bin_inds )
- deallocate( bin_fval )
- deallocate( bin_area )
- deallocate( bin_areacum )
-
- end subroutine llgrid_EquivLat
-
-
-
- !
- ! eqvlat1, eqvlat2 : lower and upper bound of eqv.lat.
- !
- subroutine llgrid_EquivLat_sort( lli, ff, eqvlatb1, eqvlatb2, status )
-
- use Binas, only : deg2rad, pi
- use Num, only : Sort
- use grid_tools, only : ll_area
-
- ! --- in/out -----------------------------
-
- type(TllGridInfo), intent(in) :: lli
- real, intent(in) :: ff(:,:)
- real, intent(out) :: eqvlatb1(:,:)
- real, intent(out) :: eqvlatb2(:,:)
- integer, intent(out) :: status
-
- ! --- const ------------------------------
-
- character(len=*), parameter :: rname = mname//'/llgrid_EquivLat_sort'
-
- ! --- local ------------------------------
-
- real :: ffsort(lli%nlon,lli%nlat)
- integer :: ijsort(lli%nlon,lli%nlat,2)
- real :: areacum
- real :: bin_areacum(lli%nlon,lli%nlat,2)
-
- real :: area_tot
- real :: area_glob, area_south
-
- integer :: i, j
- integer :: is, js
-
- ! --- begin ------------------------------
- ! check input
- call Check( lli, 'n', ff, status )
- if (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; return; end if
-
- ! sort input array
- call Sort( ff, ffsort, ijsort, status )
- if (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; return; end if
-
- ! loop over sorted values to compute cumulative areas:
- areacum = 0.0
- do js = 1, lli%nlat
- do is = 1, lli%nlon
-
- ! original indices
- i = ijsort(is,js,1)
- j = ijsort(is,js,2)
-
- ! store lower cumulative area:
- bin_areacum(i,j,1) = areacum ! rad^2
-
- ! add area of this cell:
- areacum = areacum + lli%area(j) ! rad^2
-
- ! store upper cumulative area:
- bin_areacum(i,j,2) = areacum ! rad^2
-
- end do
- end do
-
- ! total area:
- area_tot = lli%nlon * sum(lli%area)
-
- ! check ...
- if ( (bin_areacum(lli%nlon,lli%nlat,2)-area_tot)/area_tot > 0.01 ) then
- write (*,'("ERROR - total area and areas in bins do not match:")')
- write (*,'("ERROR - total area : ",es12.2)') area_tot
- write (*,'("ERROR - collected area : ",es12.2)') bin_areacum(lli%nlon,lli%nlat,2)
- write (*,'("ERROR in ",a)') rname; status=1; return
- end if
-
- !
- ! blat(nlat) ------------------------------- 90
- ! | | \
- ! |----------| |
- ! |//////////| |
- ! |----------| | glob
- ! | | \ |
- ! | | | south |
- ! | | / /
- ! blat(nlat) ------------------------------- -90
- ! blon(0) blon(nlon)
-
- ! global area and southern part (all in rad!)
- area_glob = ll_area( 0.0, lli%blon(lli%nlon)-lli%blon(0), -0.5*pi, 0.5*pi ) ! rad^2
- area_south = ll_area( 0.0, lli%blon(lli%nlon)-lli%blon(0), -0.5*pi, lli%blat(0) ) ! rad^2
- ! convert cumulative ax to (part of) [-1,1]
- bin_areacum = (area_south + bin_areacum ) / area_glob ! [ 0,1]
- bin_areacum = -1.0 + 2.0*bin_areacum ! [-1,1]
- ! adhoc truncations ...
- where ( bin_areacum < -1.0 ) bin_areacum = -1.0
- where ( bin_areacum > 1.0 ) bin_areacum = 1.0
-
- ! convert to latitude:
- eqvlatb1 = asin( bin_areacum(:,:,1) ) / deg2rad
- eqvlatb2 = asin( bin_areacum(:,:,2) ) / deg2rad
- ! ok
- status = 0
-
- end subroutine llgrid_EquivLat_sort
-
-
- ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
- ! ~~~
- ! ~~~ region box
- ! ~~~
- ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-
- subroutine llreg_Init( llreg, west_deg, east_deg, south_deg, north_deg, status )
-
- use Binas, only : deg2rad
-
- ! --- in/out --------------------------------
-
- type(TllRegion), intent(inout) :: llreg
- real, intent(in) :: west_deg, east_deg, south_deg, north_deg
- integer, intent(out) :: status
-
- ! --- const --------------------------------
-
- character(len=*), parameter :: rname = 'llreg_Init'
-
- ! --- begin --------------------------------
-
- ! store region boundaries:
- llreg%west_deg = west_deg
- llreg%east_deg = east_deg
- llreg%south_deg = south_deg
- llreg%north_deg = north_deg
-
- ! idem in radians:
- llreg%west = west_deg * deg2rad
- llreg%east = east_deg * deg2rad
- llreg%south = south_deg * deg2rad
- llreg%north = north_deg * deg2rad
-
- ! ok
- status = 0
-
- end subroutine llreg_Init
-
- ! ***
- subroutine llreg_Done( llreg, status )
-
- ! --- in/out --------------------------------
-
- type(TllRegion), intent(inout) :: llreg
- integer, intent(out) :: status
-
- ! --- const --------------------------------
-
- character(len=*), parameter :: rname = 'llreg_Done'
-
- ! --- begin --------------------------------
-
- ! nothing to be done
-
- ! ok
- status = 0
-
- end subroutine llreg_Done
-
-
- ! ***
-
-
- subroutine llreg_Region_Apply_Factor_2d( lli, x, llreg, fac, status, complement )
-
- use grid_tools, only : ll_area_frac
-
- ! --- in/out ---------------------------
-
- type(TllGridInfo), intent(in) :: lli
- real, intent(inout) :: x(:,:)
- type(TllRegion), intent(in) :: llreg
- real, intent(in) :: fac
- integer, intent(out) :: status
- logical, intent(in), optional :: complement
-
- ! --- const --------------------------------
-
- character(len=*), parameter :: rname = 'llreg_Region_Apply_Factor_2d'
-
- ! --- local --------------------------------
-
- logical :: do_complement
- real :: frac
- integer :: i, j
-
- ! --- begin --------------------------------
-
- ! apply to interior of region or to complement ?
- do_complement = .false.
- if ( present(complement) ) do_complement = complement
-
- ! check ...
- call Check( lli, 'n', x, status )
- if (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; return; end if
- ! loop over all cells in grid:
- do j = 1, lli%nlat
- do i = 1, lli%nlon
-
- ! fraction of grid cell covered by region;
- ! input to routine in radians:
- frac = ll_area_frac( lli%blon(i-1), lli%blon(i), lli%blat(j-1), lli%blat(j), &
- llreg%west , llreg%east , llreg%south , llreg%north )
-
- ! apply to interior or to complement ?
- if ( do_complement ) then
- ! apply factor if cell is (partly) outside region box:
- if ( frac < 1.0 ) then
- x(i,j) = fac * ( 1.0 - frac ) * x(i,j)
- end if
- else
- ! apply factor if cell is (partly) inside region box:
- if ( frac > 0.0 ) then
- x(i,j) = fac * frac * x(i,j)
- end if
- end if
-
- end do ! i
- end do ! j
-
- ! ok
- status = 0
-
- end subroutine llreg_Region_Apply_Factor_2d
-
-
-
- ! ***
-
-
- subroutine llreg_Region_Apply_Factor_3d( lli, x, llreg, fac, status, complement )
-
- use grid_tools, only : ll_area_frac
-
- ! --- in/out ---------------------------
-
- type(TllGridInfo), intent(in) :: lli
- real, intent(inout) :: x(:,:,:)
- type(TllRegion), intent(in) :: llreg
- real, intent(in) :: fac
- integer, intent(out) :: status
- logical, intent(in), optional :: complement
-
- ! --- const --------------------------------
-
- character(len=*), parameter :: rname = 'llreg_Region_Apply_Factor_3d'
-
- ! --- local --------------------------------
-
- logical :: do_complement
- real :: frac
- integer :: i, j
-
- ! --- begin --------------------------------
-
- ! apply to interior of region or to complement ?
- do_complement = .false.
- if ( present(complement) ) do_complement = complement
-
- ! check ...
- call Check( lli, 'n', x(:,:,1), status )
- if (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; return; end if
-
- ! loop over all cells in grid:
- do j = 1, lli%nlat
- do i = 1, lli%nlon
-
- ! fraction of grid cell covered by region;
- ! input to routine in radians:
- frac = ll_area_frac( lli%blon(i-1), lli%blon(i), lli%blat(j-1), lli%blat(j), &
- llreg%west , llreg%east , llreg%south , llreg%north )
-
- ! apply to interior or to complement ?
- if ( do_complement ) then
- ! apply factor if cell is (partly) outside region box:
- if ( frac < 1.0 ) then
- x(i,j,:) = fac * ( 1.0 - frac ) * x(i,j,:)
- end if
- else
- ! apply factor if cell is (partly) inside region box:
- if ( frac > 0.0 ) then
- x(i,j,:) = fac * frac * x(i,j,:)
- end if
- end if
-
- end do ! i
- end do ! j
-
- ! ok
- status = 0
-
- end subroutine llreg_Region_Apply_Factor_3d
-
-
-
- end module grid_type_ll
-
|