1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677787980818283848586878889909192939495969798991001011021031041051061071081091101111121131141151161171181191201211221231241251261271281291301311321331341351361371381391401411421431441451461471481491501511521531541551561571581591601611621631641651661671681691701711721731741751761771781791801811821831841851861871881891901911921931941951961971981992002012022032042052062072082092102112122132142152162172182192202212222232242252262272282292302312322332342352362372382392402412422432442452462472482492502512522532542552562572582592602612622632642652662672682692702712722732742752762772782792802812822832842852862872882892902912922932942952962972982993003013023033043053063073083093103113123133143153163173183193203213223233243253263273283293303313323333343353363373383393403413423433443453463473483493503513523533543553563573583593603613623633643653663673683693703713723733743753763773783793803813823833843853863873883893903913923933943953963973983994004014024034044054064074084094104114124134144154164174184194204214224234244254264274284294304314324334344354364374384394404414424434444454464474484494504514524534544554564574584594604614624634644654664674684694704714724734744754764774784794804814824834844854864874884894904914924934944954964974984995005015025035045055065075085095105115125135145155165175185195205215225235245255265275285295305315325335345355365375385395405415425435445455465475485495505515525535545555565575585595605615625635645655665675685695705715725735745755765775785795805815825835845855865875885895905915925935945955965975985996006016026036046056066076086096106116126136146156166176186196206216226236246256266276286296306316326336346356366376386396406416426436446456466476486496506516526536546556566576586596606616626636646656666676686696706716726736746756766776786796806816826836846856866876886896906916926936946956966976986997007017027037047057067077087097107117127137147157167177187197207217227237247257267277287297307317327337347357367377387397407417427437447457467477487497507517527537547557567577587597607617627637647657667677687697707717727737747757767777787797807817827837847857867877887897907917927937947957967977987998008018028038048058068078088098108118128138148158168178188198208218228238248258268278288298308318328338348358368378388398408418428438448458468478488498508518528538548558568578588598608618628638648658668678688698708718728738748758768778788798808818828838848858868878888898908918928938948958968978988999009019029039049059069079089099109119129139149159169179189199209219229239249259269279289299309319329339349359369379389399409419429439449459469479489499509519529539549559569579589599609619629639649659669679689699709719729739749759769779789799809819829839849859869879889899909919929939949959969979989991000100110021003100410051006100710081009101010111012101310141015101610171018101910201021102210231024102510261027102810291030103110321033103410351036103710381039104010411042104310441045104610471048104910501051105210531054105510561057105810591060106110621063106410651066106710681069107010711072107310741075107610771078107910801081108210831084108510861087108810891090109110921093109410951096109710981099110011011102110311041105110611071108110911101111111211131114111511161117111811191120112111221123112411251126112711281129113011311132113311341135113611371138113911401141114211431144114511461147114811491150115111521153115411551156115711581159116011611162116311641165116611671168116911701171117211731174117511761177117811791180118111821183118411851186118711881189119011911192119311941195119611971198119912001201120212031204120512061207120812091210121112121213121412151216121712181219122012211222122312241225122612271228122912301231123212331234123512361237123812391240124112421243124412451246124712481249125012511252125312541255125612571258125912601261126212631264126512661267126812691270127112721273127412751276127712781279128012811282128312841285128612871288128912901291129212931294129512961297129812991300130113021303130413051306130713081309131013111312131313141315131613171318131913201321132213231324132513261327132813291330133113321333133413351336133713381339134013411342134313441345134613471348134913501351135213531354135513561357135813591360136113621363136413651366136713681369137013711372137313741375137613771378137913801381138213831384138513861387138813891390139113921393139413951396139713981399140014011402140314041405140614071408140914101411141214131414141514161417141814191420142114221423142414251426142714281429143014311432143314341435143614371438143914401441144214431444144514461447144814491450145114521453145414551456145714581459146014611462146314641465146614671468146914701471147214731474147514761477147814791480148114821483 |
- !#################################################################
- !BOP
- !
- ! !MODULE: DIFFUSION
- !
- ! !DESCRIPTION:
- !
- ! Diffusion routines.
- !
- !
- ! !REVISION HISTORY: -----------------------------
- !
- ! ????-?? Bram Bregman.
- ! Adjusted for TM5 by
- !
- ! 2005-04 ??
- ! Adjoint version
- !
- ! 2010-08 Arjo Segers
- ! Debugged computation of u* over sea:
- ! abs. surface stress should be devided by air density.
- ! Set minimum value for u* to trap division by zero (MK).
- !
- ! 2012-01 Philippe Le Sager
- ! Adapted for Lat-Lon MPI domain decomposition
- !
- !EOP
- !#################################################################
- !
- ! global input fields:
- ! phlb(:,:,lm+1) ! half level pressure [Pa]
- ! gph(:,:,lm+1) ! height of half-levels [m] !CMK note gph starts NOW at 1!
- ! t ! temperature [K]
- ! q ! specific humidity [kg/kg]
- ! uwind ! mean wind W-E [m/s]
- ! vwind ! mean wind S-N [m/s]
- ! sshf ! surface sensible heat flux [W / m^2]
- ! slhf ! surface latent heat flux [W / m^2]
- ! ustar (ustar_loc) ! velocity scale
- !
- ! global output fields:
- ! kvh ! eddy diffusivity for heat at top [m^2/s]
- ! pblh ! boundary layer height[m]
- ! dkg ! vertical diffusion coefficient [ks s-1]
- !
- !### macro's #####################################################
- !
- #define TRACEBACK write (gol,'("in ",a," (",a,", line",i5,")")') rname, __FILE__, __LINE__; call goErr
- #define IF_NOTOK_RETURN(action) if (status/=0) then; TRACEBACK; action; return; end if
- #define IF_ERROR_RETURN(action) if (status> 0) then; TRACEBACK; action; return; end if
- !
- #include "tm5.inc"
- !
- !#################################################################
- MODULE DIFFUSION
- use GO, only : gol, goPr, goErr
- use dims, only : nregions
- use global_types, only : emis_data
- implicit none
- private
- ! methods
- public :: Diffusion_Init, Diffusion_Done
- public :: Calc_Kzz, Read_Diffusion, Write_Diffusion
- ! flag
- public :: diffusion_write
-
- ! --- const --------------------------------------
- character(len=*), parameter :: mname = 'diffusion'
- logical :: diffusion_write = .false.
- character(len=2000) :: diffusion_dir = ''
- ! --- local ------------------------------------
- type(emis_data), dimension(nregions), target :: ustar_loc, sr_mix
- CONTAINS
- ! ================================================================
-
-
- subroutine Diffusion_Init( status )
-
- use global_data, only : rcfile
- use MeteoData, only : Set
- use MeteoData, only : pu_dat, pv_dat
- use MeteoData, only : temper_dat, humid_dat, gph_dat
- use MeteoData, only : oro_dat, lsmask_dat
- use MeteoData, only : sr_ecm_dat, sr_ols_dat
- use MeteoData, only : u10m_dat, v10m_dat
- use MeteoData, only : slhf_dat, sshf_dat
- use MeteoData, only : ewss_dat, nsss_dat
- use GO, only : TrcFile, Init, Done, ReadRc
- ! --- in/out --------------------------------
-
- integer, intent(out) :: status
-
- ! --- const ------------------------------
- character(len=*), parameter :: rname = mname//'/Diffusion_Init'
-
- ! --- local -------------------------------------
-
- integer :: region
- type(TrcFile) :: rcF
-
- ! --- begin --------------------------------
-
- ! loop over all (zoom) regions:
- do region = 1, nregions
- ! meteo used by diffusion
- call Set( pu_dat(region), status, used=.true. )
- call Set( pv_dat(region), status, used=.true. )
- call Set( temper_dat(region), status, used=.true. )
- call Set( humid_dat(region), status, used=.true. )
- call Set( gph_dat(region), status, used=.true. )
- call Set( oro_dat(region), status, used=.true. )
- call Set( lsmask_dat(region), status, used=.true. )
- call Set( sr_ecm_dat(region), status, used=.true. )
- call Set( sr_ols_dat(region), status, used=.true. )
- call Set( u10m_dat(region), status, used=.true. )
- call Set( v10m_dat(region), status, used=.true. )
- call Set( sshf_dat(region), status, used=.true. )
- call Set( slhf_dat(region), status, used=.true. )
- call Set( ewss_dat(region), status, used=.true. )
- call Set( nsss_dat(region), status, used=.true. )
- end do
-
- call Init( rcF, rcfile , status)
- IF_NOTOK_RETURN(status=1)
- call ReadRc(rcF, 'diffusion.write', diffusion_write, status, default=.false.)
- IF_ERROR_RETURN(status=1)
- if (diffusion_write) then
- call ReadRc(rcF, 'diffusion.dir', diffusion_dir, status)
- IF_NOTOK_RETURN(status=1)
- endif
-
- call Done(rcF, status)
- IF_NOTOK_RETURN(status=1)
-
- ! ok
- status = 0
-
- end subroutine Diffusion_Init
-
-
- ! ***
-
-
- subroutine Diffusion_Done( status )
-
- ! --- in/out --------------------------------
-
- integer, intent(out) :: status
-
- ! --- const ------------------------------
- character(len=*), parameter :: rname = mname//'/Diffusion_Done'
- ! --- begin --------------------------------
-
- ! nothing to be done
- ! ok
- status = 0
-
- end subroutine Diffusion_Done
-
-
-
- ! ================================================================
-
- !------------------------------------
- !
- ! Purpose:
- ! -------
- ! this subroutine reads and prepares all fields needed for the calculation of kzz
- !
- ! External
- ! --------
- ! dd_get_3_hourly_surface_fields
- ! dd_calc_ustar_raero_rb
- ! dd_coarsen_fields
- !
- ! Reference
- ! ---------
- ! Ganzeveld and Lelieveld (1996) and references therein.
- ! Adjusted for TM5, Bram Bregman, August 2003
- !
- !
- subroutine calc_kzz( status )
- use binas, only : grav
- use dims, only : im, jm, lmax_conv, idate, lm, okdebug
- use dims, only : nread, ndyn, tref, at, bt
- use dims, only : revert
- use global_data, only : mass_dat, wind_dat, region_dat
- #ifndef without_convection
- #ifndef without_diffusion
- use global_data, only : conv_dat
- #endif
- #endif
- use meteoData, only : spm_dat, pu_dat, pv_dat, m_dat
- use meteoData, only : temper_dat, humid_dat, gph_dat
- use tm5_distgrid, only : dgrid, Get_DistGrid, update_halo, dist_arr_stat
- integer, intent(out) :: status
-
- ! --- local ------------------------------
- character(len=10) :: c_time
- real,dimension(:,:,:),pointer :: phlb,gph,t,q,m,pu,pv
- real,dimension(:,:,:),pointer :: am,bm,cm ! fluxes in kg/timestep(region)
- #ifndef without_convection
- #ifndef without_diffusion
- real,dimension(:,:,:),pointer :: dkg !vert diff (ks/s)
- #endif
- #endif
- real,dimension(:),pointer :: dxyp
- integer,parameter :: avg_field=1, nlon360=360, nlat180=180
- integer :: i,region, nsend, ntimes
- real, dimension(:,:,:),allocatable,target :: m_adj, phlb_adj
- real, dimension(:,:),allocatable :: p_adj
- integer :: imr, jmr, lmr, j,l
- real, allocatable, target :: phlb_t(:,:,:)
- real, allocatable, target :: m_t(:,:,:)
-
- integer :: i1, i2, j1, j2
- character(len=*), parameter :: rname = mname//'/Calc_kzz'
- ! --- begin ----------------------------------
- if (okdebug) write(*,*) 'start of calc_kzz'
- ! Allocate work arrays, and fill ghost cells
- DO region = 1, nregions
- CALL GET_DISTGRID( DGRID(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 )
- ALLOCATE(ustar_loc(region)%surf( i1:i2, j1:j2))
- ALLOCATE( sr_mix(region)%surf( i1:i2, j1:j2))
- CALL UPDATE_HALO( dgrid(region), wind_dat(region)%am, wind_dat(region)%halo, status)
- IF_NOTOK_RETURN(status=1)
-
- CALL UPDATE_HALO( dgrid(region), wind_dat(region)%bm, wind_dat(region)%halo, status)
- IF_NOTOK_RETURN(status=1)
- CALL UPDATE_HALO( dgrid(region), wind_dat(region)%cm, wind_dat(region)%halo, status)
- IF_NOTOK_RETURN(status=1)
- CALL UPDATE_HALO( dgrid(region), pu_dat(region)%data, pu_dat(region)%halo, status)
- IF_NOTOK_RETURN(status=1)
- CALL UPDATE_HALO( dgrid(region), pv_dat(region)%data, pv_dat(region)%halo, status)
- IF_NOTOK_RETURN(status=1)
-
- ENDDO
- ! -- for output of time header on debug fields
- write(c_time,'(i4,3(i2))') idate(1:4) !CMK corrected
- do i=1,10
- if (c_time(i:i)==' ') c_time(i:i)='0'
- enddo
- !
- ! Surface data sets may have the following characteristics
- ! 1. instanteneous. The time written in the file is valid from t-1.5 to t+1.5
- ! 2. Accumulated. The time written in the file is valid from t to t+3
- ! 3. Daily averaged. The time on the file is always 00 hours, and valid until t+24
- !
- ! *** It is essential to understand this timing when reading and applying these sets.
- !
- !
- ! fd2mk it has to be decided to 'save' fields like vgrat_low and daily average surface fields
- ! for later use, or to read it every 3 hours time step
- REG: do region = 1,nregions
- ! local grid sizes
- lmr=lm(region)
- call Get_DistGrid( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 )
- if (revert == -1) then
- !mkadj
- !mass and pressure are at t=t+dt in adjoint
- !these should be brought back to t to get
- !exactly the same diffusion as in forward run
- allocate ( m_adj (i1:i2, j1:j2, lmr ) )
- allocate ( phlb_adj(i1:i2, j1:j2, lmr+1) )
- allocate ( p_adj (i1:i2, j1:j2 ) )
-
- am => wind_dat(region)%am
- bm => wind_dat(region)%bm
- cm => wind_dat(region)%cm
- dxyp => region_dat(region)%dxyp
-
- ntimes = (nread/(ndyn))*tref(region)
- m_adj = m_dat(region)%data(i1:i2, j1:j2, 1:lmr)
-
- do i=1,ntimes
- m_adj = m_adj + revert*am(i1-1:i2-1, j1 :j2 , 1:lmr ) & !east
- - revert*am(i1 :i2 , j1 :j2 , 1:lmr ) & !west
- + revert*bm(i1 :i2 , j1 :j2 , 1:lmr ) & !south
- - revert*bm(i1 :i2 , j1+1:j2+1, 1:lmr ) & !north
- + revert*cm(i1 :i2 , j1 :j2 , 0:lmr-1) & !lower
- - revert*cm(i1 :i2 , j1 :j2 , 1:lmr ) !upper
- enddo
- p_adj(:,:) = 0.0
-
- do l=1,lmr
- ! note that on the EDGES masses are coarse>
- ! diffusion iis applied only on core zoom>
- do j = j1, j2
- do i = i1, i2
- p_adj(i,j) = p_adj(i,j) + m_adj(i,j,l)*grav/dxyp(j)
- end do
- end do
- end do
- do l=1,lmr+1
- do j = j1, j2
- do i = i1, i2
- phlb_adj(i,j,l) = at(l)+bt(l)*p_adj(i,j)
- end do
- end do
- end do
- deallocate(p_adj)
- nullify(am)
- nullify(bm)
- nullify(cm)
- nullify(dxyp)
- phlb => phlb_adj
- m => m_adj
- else
- ! Bug fix 2006-09-21 (AJS)
- ! The two half level pressure arrays phbl_t and phlb_k
- ! are sometimes not updated together when running in parallel
- ! with zoom regions. To avoid problems, pressure and mass
- ! are computed here given the surface pressure in the middle
- ! of the current time interval, stored in spm_dat.
- !--
- !phlb => mass_dat(region)%phlb_t
- !m => mass_dat(region)%m_t
- !--
- allocate( phlb_t(i1:i2, j1:j2, 1:lmr+1) )
- do l = 1, lmr+1
- phlb_t(:,:,l) = at(l) + bt(l)*spm_dat(region)%data(i1:i2, j1:j2, 1)
- end do
- allocate( m_t(i1-2:i2+2, j1-2:j2+2, 1:lmr) )
- m_t = 0.0
- do l = 1, lmr
- do j = j1, j2
- m_t(i1:i2,j,l) = ( phlb_t(:,j,l) - phlb_t(:,j,l+1) )*region_dat(region)%dxyp(j)/grav
- end do
- end do
- phlb => phlb_t
- m => m_t
- !--
- endif
- pu => pu_dat(region)%data
- pv => pv_dat(region)%data
- gph => gph_dat(region)%data
- t => temper_dat(region)%data
- q => humid_dat(region)%data
- #ifndef without_convection
- #ifndef without_diffusion
- dkg => conv_dat(region)%dkg
- #endif
- #endif
- ! ustar_loc
- call dd_calc_ustar(region, i1, i2, j1, j2)
- if (okdebug) then
- ! call dd_field_statistics(region, 'ustar_loc', ustar_loc(region)%surf, i1, j1, status)
- call dist_arr_stat(dgrid(region), 'ustar_loc', ustar_loc(region)%surf, 0, status)
- end if
-
- #ifndef without_convection
- #ifndef without_diffusion
- ! calculate kzz
- call bldiff( region, phlb, m, i1, i2, j1, j2 )
- #endif
- #endif
- nullify(phlb)
- nullify(m)
- !--
- deallocate( phlb_t )
- deallocate( m_t )
- !--
- nullify(pu)
- nullify(pv)
- #ifndef without_convection
- #ifndef without_diffusion
- nullify(dkg)
- #endif
- #endif
- nullify(t)
- nullify(q)
- nullify(gph)
- if(revert == -1) then
- deallocate(m_adj)
- deallocate(phlb_adj)
- endif
- end do REG
- ! Done
- do region = 1,nregions
- deallocate(sr_mix(region)%surf)
- deallocate(ustar_loc(region)%surf)
- end do
- if (okdebug) write(*,*) 'end of calc_kzz'
- end subroutine calc_kzz
-
- !------------------------------------------------------------------------------
- ! TM5 !
- !------------------------------------------------------------------------------
- !BOP
- !
- ! !IROUTINE: dd_calc_ustar
- !
- ! !DESCRIPTION:
- !
- ! Calculate friction velocity (ustar)
- !
- ! Uses the log normal wind profile information stored by ECMWF in 10 meter wind
- ! to estimate a local ustar over land
- ! uses Charnock equation over sea to estimate ustar
- ! aerodynamic resistance is calculated from heat fluxes and ustar using a constant reference height
- ! sub laminar resistance from ustar
- !
- ! Reference:
- ! o M.Z. Jacobson, "Fundamentals of Atmospheric Modeling", sections 8.3-8.4
- !
- ! !REVISION HISTORY:
- !
- ! 2003 ??; Ge Verver (personal communication, 2003)
- !
- ! 2010-08 Arjo Segers
- ! Debugged computation of u* over sea:
- ! abs. surface stress should be devided by air density.
- ! Set minum value for u* to trap division by zero (MK).
- !
- ! 2012-01 P. Le Sager
- ! adapted for lon-lat MPI domain decomposition
- !------------------------
-
- subroutine dd_calc_ustar(region, i1, i2, j1, j2)
- use binas, only : grav, vKarman
- use dims, only : okdebug
- use meteoData, only : lsmask_dat, sr_ecm_dat, sr_ols_dat
- use meteoData, only : u10m_dat, v10m_dat
- use meteoData, only : slhf_dat, sshf_dat
- use meteoData, only : ewss_dat, nsss_dat
- !--- in/out ------------------------------------
- integer,intent(in) :: region, i1, i2, j1, j2
-
- ! --- const -------------------------------------
- ! Garret relation:
- real,parameter :: alfa_garret = 0.11
- ! air density at seal level and 293 K :
- real, parameter :: rho_air = 1.2 ! kg/m3
- ! Charnock relation:
- real,parameter :: alfa_charnock = 0.018
- ! kinematic viscosity of air at about 300 K :
- real,parameter :: nu_air = 1.5e-5 ! m2/s
-
- real, parameter :: Href=30. ! constant reference height for calculations of raero
- ! some constants specific to calculation of ustar and raero
- real,parameter :: rhoCp = 1231.0
- real,parameter :: rhoLv = 3013000.0
-
- real, parameter :: rz0 =2.0
- !--- local -------------------------------------
- integer :: i,j
- real :: tau_z_abs
- real :: ustar_sea,ustar_land,xland,sr_sea,sr_land,sr_help,wind10m
- ! --- begin ----------------------------------
- do j=j1, j2
- do i=i1, i2
- xland=lsmask_dat(region)%data(i,j,1)/100. ! 0-1
- ! SEA:
- !
- ! From M.Z. Jacobson, "Fundamentals of Atmospheric Modeling",
- ! section "8.3 Friction velocity" :
- ! The friction velocity :
- ! u* = sqrt( |tau_z|/rho_a ) [m/s]
- ! where:
- ! tau_z = (ewss,nsss) surface stress [N/m2=kg/m/s2]
- ! rho_air air density [kg/m3] ; about 1.2 at sea level and 293 K
- !
- tau_z_abs = sqrt(ewss_dat(region)%data(i,j,1)**2+&
- nsss_dat(region)%data(i,j,1)**2 ) ! N/m2
- ustar_sea = sqrt(tau_z_abs/rho_air) ! m/s
- !
- ! limitation to avoid division by zero:
- ustar_sea = max( 0.001, ustar_sea ) ! m/s ; minium of 1 mm/s
- !
- ! From M.Z. Jacobson, "Fundamentals of Atmospheric Modeling",
- ! section "8.4 Surface roughness lengths" :
- ! "For smooth surfaces, such as over a smooth ocean with low wind speeds (Garret):
- ! z0m = 0.11 nu_a / u*
- ! where nu_a is the kinematic viscosity of air.
- ! Over a rough ocean with high wind speeds:
- ! z0m = alpha_c (u*)**2 / grav
- ! which is the 'Charnock relation',
- ! where alpha_c ~ 0.016 is the 'Charnock constant'."
- !
- ! Here the sum of these two parameterizations is used,
- ! where the first part is dominant for u*<0.1 and the second for u*>0.1 ;
- ! minimum value for u* prevents division by zero:
- !
- sr_sea = alfa_garret * nu_air / ustar_sea + &
- alfa_charnock * ustar_sea**2 / grav
- !
- ! LAND
- ! calculate the 'local' surface roughness for momentum that is consistent with 10 m wind
- ! and the Olsson data base, we assume that the windspeed at 75 m is independent of
- ! surface roughness
- !
- wind10m = sqrt( u10m_dat(region)%data(i,j,1)**2 + v10m_dat(region)%data(i,j,1)**2 )
- if ( xland > 0.0 ) then
- !>>> TvN
- ! Over land, the friction velocity ustar
- ! is calculated from the 10 m wind speed, u10.
- ! In IFS u10 is a diagnostic variable, calculated to be
- ! compatible with "open-terrain" wind speed observations.
- ! Over rough or inhomogeneous terrain (z0M > 0.03 m),
- ! it is calculated using an aerodynamic roughness
- ! length typical for open terrain with grassland (0.03 m)
- ! (see IFS cy31r1 documentation, p. 46).
- ! In the expressions applied in TM5,
- ! the stability profile functions Psi_M have disappeared,
- ! and only the logarithmic terms are kept.
- ! Apart from this, the expression for ustar was incorrect:
- ! 10./sr_land should be 75./sr_land.
- ! This has now been corrected.
- ! Moreover, over islands and near coast lines,
- ! zeros can occur in sr_ols_dat,
- ! which have to be removed.
- ! However, a lower bound of 1e-2 m = 1 cm
- ! seems too high for smooth surfaces,
- ! like deserts and ice caps.
- ! The minimum positive value in sr_ols_dat
- ! is 2.5e-4 m = 0.025 cm,
- ! which is obtained in parts of Antarctica.
- ! This is likely the result of regridding of
- ! a field with minimum value of 0.1 cm
- ! from 0.5x0.5 to 1x1 degrees.
- ! Please note that with the introduction of CY31R1,
- ! the orographic contribution to the aerodynamic roughness length
- ! in IFS has been replaced by an explicit specification of stress
- ! on model levels due to turbulent orographic form drag,
- ! and the climatological variable previously used
- ! has been replaced by a prognostic variable.
- ! In TM5, the climatological variable is still used,
- ! but it would be better to use the prognostic variable instead.
- ! Note also that the same calculation
- ! is done in dry_deposition.F90.
- !sr_land=max(1e-2,sr_ols_dat(region)%data(i,j,1)) ! occurs at Islands, etc.
- sr_land=max(1e-3,sr_ols_dat(region)%data(i,j,1))
- sr_help=min(sr_ecm_dat(region)%data(i,j,1),0.03)
- !fd ustar_land=vKarman*wind10m(i,j)/alog(10./sr_help) !ustar consistent with 'clipped' large scale roughness
- !ustar_land=vKarman*wind10m/&
- ! alog(10./sr_help)*alog(75./sr_help)/alog(10./sr_land)
- ustar_land=vKarman*wind10m/&
- alog(10./sr_help)*alog(75./sr_help)/alog(75./sr_land)
- ! <<< TvN
- else
- sr_land = 0.0
- ustar_land = 0.0
- endif
- ! interpolate between land and sea for mixed pixels:
- ustar_loc(region)%surf(i,j) = xland*ustar_land+(1-xland)*ustar_sea
- sr_mix (region)%surf(i,j) = xland* sr_land+(1-xland)* sr_sea
- end do !i
- end do !j
- if (okdebug) write(*,*) 'end calc_ustar'
-
- end subroutine dd_calc_ustar
-
-
- ! ***
-
-
- #ifndef without_convection
- #ifndef without_diffusion
- !------------------------------------------------------------------
- ! Procedure:
- ! (1) calc_Kv: Diffusion coefficients full atmosphere from shear
- ! (2) pblhght: Determine the boundary layer height
- ! (3) difcoef: Replace diffusion coefficients in the boundary layer
- !------------------------------------------------------------------
- ! subroutine adjusted for TM5, Bram Bregman, August 2003.
- !-------------------------IO---------------------------------------
-
- subroutine bldiff( region, phlb, m, i1, i2, j1, j2 )
- use binas, only: ae, cp_air, Rgas, grav, Lvap
- use binas, only: xmair
- use binas, only: vkarman
- use dims, only: im, jm, lm, lmax_conv
- use dims, only: okdebug
- use dims, only: gtor, dx, dy, ybeg, xref, yref
- use dims, only: isr, ier, jsr, jer
- use toolbox, only: dumpfield
- use global_data, only : conv_dat
- use meteodata , only : slhf_dat, sshf_dat, pu_dat, pv_dat, gph_dat, temper_dat, humid_dat
- ! --- in/out ----------------------------------------------
- integer,intent(in) :: region, i1, i2, j1, j2
- real,dimension(:,:,:),pointer :: phlb, gph, t, q, m, pu, pv, dkg
- ! --- const ----------------------------------
- ! minimum value diffusion coefficient
- real, parameter :: ckmin = 1.e-15
- real, parameter :: sffrac= 0.1
- real, parameter :: onet = 1./3.
- real, parameter :: rair = Rgas*1.e3/xmair
- real, parameter :: ccpq = 1869.46/cp_air -1.
- real, parameter :: epsilo = rair/461.51
- real, parameter :: apzero = 100000. ! reference pressure (usually 1000 hPa)
- real, parameter :: ccon = sffrac*vkarman
- real, parameter :: betam=15.,betah=15.
- real, parameter :: binm = betam*sffrac
- real, parameter :: binh = betah*sffrac
- real, parameter :: fr_excess = 1.
- real, parameter :: fak = fr_excess * 8.5
- real, parameter :: fakn = fr_excess * 7.2
- real, parameter :: zkappa = rair / cp_air
- real, parameter :: zcrdq = 1.0/epsilo - 1.0
- real, parameter :: tiny = 1.e-9 ! bound wind**2 to avoid dividing by 0
- real, parameter :: zacb = 100. !factor in ustr-shearproduction term
- real, parameter :: ricr = 0.3 ! critical richardson number
- real, parameter :: ssfrac=0.1
- ! --- local ------------------------------------
- real, dimension(:,:,:), allocatable:: wkvf,& ! outer layer Kv (at the top of each layer)
- zdup2,& ! vertical shear squared
- zrinub,& ! Richardson number
- wthm,& ! potential temperature
- zthvk,& ! virtual potential temperature
- pf,& ! full level pressure
- whm, & ! height of layer centers
- uwind,vwind,& ! wind velocities [m s-1]
- kvh !
-
- real, dimension(:,:), allocatable :: ustr,& ! velocity scale
- wheat,& ! surface sensible heat flux [K m/s]
- wqflx,& ! surface latent heat flux [kg m/ (kg s)]
- wobkl,& ! Monin Obukhov length
- wheatv,& ! surface virtual heat flux [K m/s]
- wtseff ! theta_v at lowest level [K]
- integer :: i, j, l
- real :: intfac !cmk BUG was integer!
- ! variables for the calculation of free atmosphere vertical diffusion coefficients wkvf
- real :: cml2,zlambdac,arg,zdz,zthva,zfunst, zfstabh,zsstab,zkvn,dyy
- real, dimension(j1:j2) :: dxx, lat
- integer :: jq! jqif function
- ! variabales for the calculation of the boundary layer height
- real :: wrino ! Richardson number
- real :: thvparc ! parcel virtual potential temperature
- real :: fmt,wsc,tkv,zrino,vvk,dtkv,zexcess
- integer :: jwork
- real,dimension(:,:),pointer :: pblh ! planetary boundary layer height [m]
- ! variables to calculate difusion coefficient in the boundary layer
- !real :: wsc,cml2
- real :: z,zwstr,zkvh,kvhmin,zm,zp,zslask,zsl1, &
- zh,zzh,zl,pr,zpblk,zfstabh1,zfstabh2,&
- term1,term2,term3,obukov
- integer :: jq1,jq2,jq3,jq4,jq5,jck,jqq, jq6, jq7
- real :: yres
- real :: xres
- real :: thvgrad, kvhentr
- ! --- begin ------------------------------------------
- if (okdebug) write(*,*) 'start of bldiff'
-
- allocate(wkvf ( i1:i2, j1:j2, lm(region)))
- allocate(zdup2 ( i1:i2, j1:j2, lm(region)))
- allocate(zrinub( i1:i2, j1:j2, lm(region)))
- allocate(wthm ( i1:i2, j1:j2, lm(region)))
- allocate(zthvk ( i1:i2, j1:j2, lm(region)))
- allocate(pf ( i1:i2, j1:j2, lm(region)))
- allocate(whm ( i1:i2, j1:j2, lm(region)))
- allocate(uwind ( i1:i2, j1:j2, lm(region)))
- allocate(vwind ( i1:i2, j1:j2, lm(region)))
- allocate(kvh ( i1:i2, j1:j2, lm(region)))
- allocate(ustr ( i1:i2, j1:j2 ))
- allocate(wheat ( i1:i2, j1:j2 ))
- allocate(wqflx ( i1:i2, j1:j2 ))
- allocate(wobkl ( i1:i2, j1:j2 ))
- allocate(wheatv( i1:i2, j1:j2 ))
- allocate(wtseff( i1:i2, j1:j2 ))
- pblh => conv_dat(region)%blh
- pu => pu_dat(region)%data
- pv => pv_dat(region)%data
- gph => gph_dat(region)%data
- t => temper_dat(region)%data
- q => humid_dat(region)%data
- dkg => conv_dat(region)%dkg
-
- ! mid-layer pressure levels
- ! potential temperature field
- ! virtual potential temperature field
- !
- !PLS $OMP PARALLEL &
- !PLS $OMP default (none) &
- !PLS $OMP shared (lm, jsr, jer, isr, ier, region, pf, phlb, t, q, wthm, &
- !PLS $OMP zthvk, whm, gph) &
- !PLS $OMP private (i, j, l, jss, jee, intfac)
- ! call my_omp_domain (jsr(region), jer(region), jss, jee)
- do l=1,lm(region)
- do j=j1,j2
- do i=i1,i2
- pf(i,j,l) = (phlb(i,j,l) + phlb(i,j,l+1) )*0.5
- wthm(i,j,l) = t(i,j,l) * ( apzero/pf(i,j,l) )**zkappa
- zthvk(i,j,l) = wthm(i,j,l) * (1.0 + zcrdq*q(i,j,l))
- ! mid-layer heights (note: gph array is for levels 1:lm+1)
- intfac = (phlb(i,j,l)-pf(i,j,l)) / (phlb(i,j,l)-phlb(i,j,l+1))
- ! cmk: intfac is now real!
- whm(i,j,l) = (gph(i,j,l)-gph(i,j,1)) * (1.0-intfac) + (gph(i,j,l+1) -gph(i,j,1)) * intfac
- ! cmk bug corrected:
- ! cmk old: whm(i,j,l) = (gph(i,j,l+1)-gph(i,j,1)) * (1.0-intfac) + (gph(i,j,l+1) -gph(i,j,1)) * intfac
- end do
- end do
- end do
- !PLS $OMP END PARALLEL
- ! convert units of sensible and latent heat fluxes
- ! sshf: W/m2 / (J/kg/K) / (kg/m3) = J/s/m2 /J kg K /kg m3 = K m/s
- ! slhf: W/m2 / (J/kg) / (kg/m3) = J/s/m2 /J kg /kg m3 = m/s
- ! ustr is local velocity scale, NOT from ECMWF surface stresses
- !
- !PLS $OMP PARALLEL &
- !PLS $OMP default (none) &
- !PLS $OMP shared (jm, im, region, wheat, sshf_dat, gph, phlb, slhf_dat, &
- !PLS $OMP wqflx, ustr, ustar_loc) &
- !PLS $OMP private (i, j, js, je)
- ! call my_omp_domain (1, jm(region), js, je)
- do j = j1,j2
- do i = i1,i2
- wheat(i,j) = -sshf_dat(region)%data(i,j,1) * (gph(i,j,2) - &
- gph(i,j,1)) * grav / (phlb(i,j,1) - &
- phlb(i,j,2)) / cp_air
- wqflx(i,j) = -slhf_dat(region)%data(i,j,1) * (gph(i,j,2) - &
- gph(i,j,1)) * grav / (phlb(i,j,1) - &
- phlb(i,j,2)) / lvap
- ustr(i,j) = max (ustar_loc(region)%surf(i,j), 0.01)
- end do
- end do
- !PLS $OMP END PARALLEL
-
- !
- !-----------------------------------------------------------------
- ! compute the free atmosphere vertical diffusion coefficients wkvf
- ! height-dependent asymptotic mixing length implemented
- ! Holtslag and Boville, 1993, J. Climate, 6, 1825-1842.
- !------------------------------------------------------------------
-
- ! first calculate wind velocities [m s-1] from the mass fluxes
-
- yres = dy/yref(region)
- xres = dx/xref(region)
- do j = j1,j2
- lat(j) = ybeg(region) + 0.5 * yres + yres * (j-1)
- enddo
- dxx(j1:j2) = ae * xres * gtor * cos(lat(j1:j2)*gtor)
- dyy = ae * yres * gtor
- !PLS $OMP PARALLEL &
- !PLS $OMP default (none) &
- !PLS $OMP shared (lm, jsr, jer, isr, ier, region, dxx, dyy, pu, pv, m, &
- !PLS $OMP uwind, vwind) &
- !PLS $OMP private (i, j, l, jss, jee)
- ! call my_omp_domain (jsr(region), jer(region), jss, jee)
- do l=1,lm(region)
- do j=j1,j2
- do i=i1,i2
- uwind(i,j,l) = dxx(j)*(pu(i,j,l) + pu(i-1,j,l))*0.5 / m(i,j,l)
- vwind(i,j,l) = dyy *(pv(i,j,l) + pv(i,j+1,l))*0.5 / m(i,j,l)
- enddo
- enddo
- enddo
- !PLS $OMP END PARALLEL
- !
- ! Initialize gradient Richardson number and free tropospheric k values
- !
- if (okdebug) write(*,*) ' Richardson number'
- !PLS $OMP PARALLEL &
- !PLS $OMP default (none) &
- !PLS $OMP shared (lm, jsr, jer, isr, ier, region, gph, zdup2, uwind, &
- !PLS $OMP vwind, whm, phlb, pf, zthvk, zrinub, wkvf) &
- !PLS $OMP private (i, j, l, jss, jee, jq, zlambdac, cml2, zdz, arg, &
- !PLS $OMP zthva, zsstab, zfunst, zfstabh, zkvn)
- ! call my_omp_domain (jsr(region), jer(region), jss, jee)
- do l=1,lm(region)-1
- do j=j1,j2
- do i=i1,i2
- !
- ! if (delta gph < 1000) then jq=1
- jq = jqif( 1000.0, (gph(i,j,l+1)-gph(i,j,1)) )
- zlambdac=jq*300.+(1-jq)*(30. + 270. * exp (1.-(gph(i,j,l+1)-gph(i,j,1))/1000.))
- cml2=(1./( 1./zlambdac + 1./(vkarman*(gph(i,j,l+1)-gph(i,j,1))) ))**2.
- !
- ! vertical shear squared,
- ! min value of (delta v)**2 prevents zero shear
- !
- zdup2(i,j,l) = (uwind(i,j,l+1)-uwind(i,j,l))**2 + &
- (vwind(i,j,l+1)-vwind(i,j,l))**2
- zdup2(i,j,l) = max(zdup2(i,j,l),1.e-10)
- zdz = whm(i,j,l+1) - whm(i,j,l)
- zdup2(i,j,l) = zdup2(i,j,l)/(zdz**2)
- !
- ! static stability (use virtual potential temperature)
- ! dv now calculated at interface height (lineair in pressure)
- !
- arg=(log(phlb(i,j,l+1))-log(pf(i,j,l)))/(log(pf(i,j,l+1))-log(pf(i,j,l)))
- zthva = zthvk(i,j,l) + arg * (zthvk(i,j,l+1)-zthvk(i,j,l))
- zsstab = grav/zthva*(zthvk(i,j,l+1) - zthvk(i,j,l))/zdz
- !
- ! gradient Richardson number
- !
- zrinub(i,j,l) = zsstab/zdup2(i,j,l)
- !
- ! stability functions
- !
- zfunst = max(1. - 18.*zrinub(i,j,l),0.)
- zfunst = sqrt(zfunst)
- !
- zfstabh = 1.0 / (1.0 + 10.0*zrinub(i,j,l)*(1.0+8.0*zrinub(i,j,l)))
- !
- ! neutral diffusion coefficient
- !
- zkvn = cml2 * sqrt(zdup2(i,j,l))
- !
- ! correction with stability functions
- !
- jq = jqif( zrinub(i,j,l),0.0 ) ! thus: if (zrinub>0) then jq=1
- wkvf(i,j,l)=jq*(max(ckmin,zkvn*zfstabh))+(1-jq)*(max(ckmin,zkvn*zfunst))
- end do
- end do
- end do
- !PLS $OMP END PARALLEL
- !
- ! Determine the boundary layer height
- !--------------------------------------------------------------------------------
- ! based on: Vogelezang and Holtslag, 1996, Bound. Layer Meteorol., 81, 245-269.
- !--------------------------------------------------------------------------------
- if (okdebug) write(*,*) ' boundary layer height'
- !PLS $OMP PARALLEL &
- !PLS $OMP default (none) &
- !PLS $OMP shared (lm, jsr, jer, isr, ier, region, wtseff, wthm, q, &
- !PLS $OMP wheatv, wheat, wqflx, wobkl, ustr, pblh, uwind, vwind, &
- !PLS $OMP whm ) &
- !PLS $OMP private (i, j, l, jss, jee, jwork, wrino, vvk, tkv, dtkv, &
- !PLS $OMP zrino, jq, fmt, wsc, zexcess, thvparc )
- ! call my_omp_domain (jsr(region), jer(region), jss, jee)
- do j=j1,j2
- do i=i1,i2
- ! compute bottom level virtual temperature and heat flux and Monin-Obukhov Length
- wtseff(i,j) = wthm(i,j,1)*(1. + zcrdq*q(i,j,1))
- wheatv(i,j) = wheat(i,j) + zcrdq*wthm(i,j,1)*wqflx(i,j)
- wobkl(i,j) = -wtseff(i,j)*ustr(i,j)**3 / (grav*vkarman*(wheatv(i,j)+sign(1.e-10,wheatv(i,j))) )
- ! initialize
- pblh(i,j) = 0.0
- jwork = 1
- wrino = 0.0
- do l=1,lm(region)
- vvk = (uwind(i,j,l)-uwind(i,j,1))**2 + (vwind(i,j,l)-vwind(i,j,1))**2 + &
- zacb*ustr(i,j)*ustr(i,j)
- vvk = max(vvk,tiny)
- tkv = wthm(i,j,l)*(1. + zcrdq*q(i,j,l))
- dtkv = tkv-wtseff(i,j)
- !
- ! Bulk Richardson number
- !
- zrino = grav*dtkv * (whm(i,j,l)-whm(i,j,1)) / (wtseff(i,j)*vvk)
- zrino = zrino+sign(1.e-10,zrino) ! prevent zrino becoming zero
- zrino = jwork*zrino
- jq = jqif(ricr,zrino) ! thus: if (zrino < ricr) then jq=1
- !
- ! calculate pblh from linear interpolation in Ri(bulk)
- !
- if ( l == 1 ) then
- pblh(i,j) = jq*pblh(i,j) + (1-jq) * &
- ( (ricr-wrino)/(zrino-wrino)*(whm(i,j,l) ) )
- else
- pblh(i,j) = jq*pblh(i,j) + (1-jq) * &
- ( whm(i,j,l-1) + (ricr-wrino)/(zrino-wrino)*(whm(i,j,l)-whm(i,j,l-1)) )
- end if
- !
- ! first time zrino > ricr we set jwork to zero, thus all further times zrino=0
- ! and pblh does not change anymore
- !
- jwork = jq*jwork
- !
- ! if pblh is found already avoid dividing by zero next level by a faked value
- ! of wrino (jwork=0 already)
- !
- wrino = zrino+(1-jwork)*0.1
- end do
- pblh(i,j) = (1-jwork)*pblh(i,j) + whm(i,j,lm(region))*jwork
- !
- ! second calculation of pblh including excess surface temperature using
- ! a convective velocity scale wsc (wm, not equal to w*!!!), using result
- ! of the first calculation of pblh
- !
- jq = jqif(wheatv(i,j),0.) ! thus: if (wheatv(i,j) > 0.) then jq=1
- jwork = jq
- fmt = ( jq*(1.0 - binm*pblh(i,j)/wobkl(i,j)) + (1.0-jq) )**onet
- wsc = ustr(i,j)*fmt
- zexcess = wheatv(i,j)*fak/wsc
- !
- ! improve pblh estimate under convective conditions using
- ! convective temperature excess (zexcess)
- !
- thvparc = wtseff(i,j)+jq*zexcess
- !
- jwork = 1
- wrino = 0.0
- !
- do l=2,lm(region)
- !
- ! Bulk Richardson number:
- ! if stable corrected with extra shear term
- ! if unstable now also corrected for temperature excess
- !
- vvk = (uwind(i,j,l)-uwind(i,j,1))**2 + (vwind(i,j,l)-vwind(i,j,1))**2 + &
- zacb*ustr(i,j)*ustr(i,j)
- vvk = max(vvk,tiny)
- tkv = wthm(i,j,l)*(1. + zcrdq*q(i,j,l))
- zrino = grav*(tkv-thvparc) * (whm(i,j,l)-whm(i,j,1))/(vvk*wtseff(i,j))
- zrino = zrino+sign(1.e-10,zrino) ! prevent zrino becoming zero
- zrino = zrino*jwork
- jq = jqif(ricr,zrino) ! thus: if (zrino < ricr) then jq=0
- !
- ! calculate pblh from linear interpolation in Ri(bulk)
- !
- pblh(i,j) = jq*pblh(i,j)+(1-jq)* (whm(i,j,l-1)+(ricr-wrino)/ &
- (zrino-wrino)*(whm(i,j,l)-whm(i,j,l-1)))
- !
- ! first time zrino > ricr we set jwork to zero, thus all further times zrino=0
- ! and pblh does not change anymore
- !
- jwork = jq*jwork
- !
- ! if pblh is found already avoid dividing by zero next level by a faked value
- ! of wrino (jwork=0 already)
- !
- wrino = zrino+(1-jwork)*0.1
- end do
- pblh(i,j) = (1-jwork)*pblh(i,j) + jwork*whm(i,j,lm(region))
- pblh(i,j) = max(pblh(i,j), 100.0) ! set minimum value for pblh
- end do
- end do
- !PLS $OMP END PARALLEL
-
-
- !
- ! Diffusion coefficients in the surface layer
- !
- !--------------------------------------------------------------------------------
- ! the revised LTG scheme (Beljaars and Viterbo, 1998)
- ! Modification - April 1, 1999: prevent zpblk to become too small in very stable conditions
- !--------------------------------------------------------------------------------
-
- if (okdebug) write(*,*) ' diffusion coeff'
- ! initialize output array with minimum value
- kvhmin = 0.1
- !
- ! Loop to calculate the diffusivity
- !
- !PLS $OMP PARALLEL &
- !PLS $OMP default (none) &
- !PLS $OMP shared (lm, jsr, jer, isr, ier, region, pblh, gph, wobkl, &
- !PLS $OMP wheatv, whm, zrinub, zdup2, kvhmin, wkvf, ustr, &
- !PLS $OMP wtseff, kvh, zthvk) &
- !PLS $OMP private (i, j, l, jss, jee, jq, z, zh, zl, zzh, jq1, jq2, jq3, &
- !PLS $OMP jq4, jq5, jq6, jq7, jqq, cml2, zfstabh, zpblk, zkvh, &
- !PLS $OMP zslask, fmt, wsc, zwstr, obukov, term1, term2, term3, &
- !PLS $OMP pr, thvgrad, kvhentr)
- ! call my_omp_domain (jsr(region), jer(region), jss, jee)
- do l=1,lm(region)-1
- do j=j1,j2
- do i=i1,i2
- jq = jqif(pblh(i,j),(gph(i,j,l+1)-gph(i,j,1))) ! if top of level hh < pblh jq=1 inside BL
- z = gph(i,j,l+1)-gph(i,j,1)
- zh = jq*z/pblh(i,j) ! only defined in BL (jq=1)
- ! zl must not be zero
- zl = jq*z/wobkl(i,j)+(1-jq) ! z/L outside BL = 1
- zzh = jq*(1.-zh)**2 ! (1-z/h)^2 only in BL
- jq1 = jqif(wheatv(i,j),0.) ! jq1 is 1 in unstable case
- jq2 = jq*(1-jq1) ! jq2 is 1 in stable BL
- jq3 = jq*jq1 ! jq3 is 1 in unstable BL
- jq6 = jqif(pblh(i,j), whm(i,j,l))
- jq7 = jqif(whm(i,j,l+1), pblh(i,j))
- jq6 = jq1*jq6*jq7 ! jq6 is 1 at at the kvh-level that is located between the
- ! mid-levels that surround pblh
- jq1 = jqif(sffrac,zh) ! jq1 is 1 in surface layer
- jq4 = jq3*jq1 ! jq4 is 1 in unstable surface layer
- jq5 = jq3*(1-jq1) ! jq5 is 1 outside unstable surface layer
- !
- ! calculate coefficients for momentum, the values for heat are obtained by dividing by Prantl number
- !
- ! stable and neutral:
- !
- cml2 = (1./( 1./(vkarman*z) + 1./450. ))**2.
- jqq = jqif(zrinub(i,j,l),0.)
- !
- ! stability function for momentum
- !
- zfstabh = jqq/(1.+10.*zrinub(i,j,l)*sqrt(1.+jqq*zrinub(i,j,l))) + (1-jqq)
- zpblk = cml2 * sqrt(zdup2(i,j,l))*zfstabh
- zkvh = jq2*max(zpblk,kvhmin)+(1-jq2)*wkvf(i,j,l) ! take free troposp. values outside BL
- !
- ! unstable case in surface layer
- !
- zslask = 1.-betah*zl
- zslask = jq4*zslask+(1-jq4)
- zpblk = (1-jq4)*zkvh+jq4*(ustr(i,j)*vkarman*z*zzh*zslask**(1./3.))
- !
- ! unstable case above surface layer
- !
- ! evaluate stability function for momentum at height z = 0.1*pblh (top of surface layer)
- !
- zslask = 1.-ssfrac*betah*pblh(i,j)/wobkl(i,j)
- zslask = jq5*zslask+(1-jq5)
- fmt = zslask**(1./3.)
- !
- ! use surface layer definition for w_m
- wsc = ustr(i,j)*fmt
- zpblk = (1-jq5)*zpblk + jq5*wsc*vkarman*z*zzh
- !
- ! Determine Prandt Number
- !
- ! Pr-number is assumed to be constant throughout the CBL, i.e. in surface and mixed layer
- ! NOTE: this is different from the original formulation
- ! zwstr not used if wheatv<0.
- !
- zwstr = (abs(wheatv(i,j))*grav*pblh(i,j)/wtseff(i,j))**(1./3.) ! bh
- obukov = jq3*wobkl(i,j)-(1-jq3)
- term1 = (1.-ssfrac*betah*pblh(i,j)/obukov)**(1./3.)
- term2 = sqrt(1.-ssfrac*betah*pblh(i,j)/obukov)
- term3 = ccon*fakn*zwstr/(ustr(i,j)*(1.-ssfrac*betah*pblh(i,j)/obukov)**(1./3.))
- pr = jq3*(term1/term2+term3)+(1-jq3)
- !
- ! NOTE that kvh applies to the top of the level
- !
- kvh(i,j,l) = jq3*max(zpblk/pr,kvhmin)+(1-jq3)*zkvh
- ! if in the entrainment zone, override with prescribed entrainment
- thvgrad = (zthvk(i,j,l+1)-zthvk(i,j,l))/(whm(i,j,l+1)-whm(i,j,l))
- kvhentr = 0.2*wheatv(i,j)/thvgrad
- kvh(i,j,l) = jq6*kvhentr + (1-jq6)*kvh(i,j,l)
- end do
- end do
- end do
- !PLS $OMP END PARALLEl
- ! calculate vertical diffusion
- !PLS $OMP PARALLEL &
- !PLS $OMP default (none) &
- !PLS $OMP shared (jsr, jer, isr, ier, region, dkg, kvh, m, gph) &
- !PLS $OMP private (i, j, l, jss, jee)
- ! call my_omp_domain (jsr(region), jer(region), jss, jee)
- do l=1,lmax_conv-1
- do j=j1,j2
- do i=i1, i2
- dkg(i,j,l)=max(0.,kvh(i,j,l))*2.*(m(i,j,l+1)+m(i,j,l))/(gph(i,j,l+2)-gph(i,j,l))**2
- ! CMK dec 2002 ----> gph changed to 1--->lm+1 boundaries
- enddo
- enddo
- enddo
- dkg(i1:i2,j1:j2,lmax_conv) = 0.0
- !PLS $OMP END PARALLEL
- deallocate(wkvf)
- deallocate(zdup2)
- deallocate(zrinub)
- deallocate(wthm)
- deallocate(zthvk)
- deallocate(pf)
- deallocate(whm)
- deallocate(ustr)
- deallocate(wheat)
- deallocate(wqflx)
- deallocate(wobkl)
- deallocate(wheatv)
- deallocate(wtseff)
- deallocate(uwind)
- deallocate(vwind)
- deallocate(kvh)
- nullify(pblh)
- if (okdebug) write(*,*) 'end of bldiff'
- end subroutine bldiff
- #endif /* diffusion */
- #endif /* convection */
- !
- ! define vectorizable 'if' function
- ! if xxxz>yyyz jqif=1 else jqif=0
- !
- INTEGER FUNCTION JQIF( xxxz, yyyz )
- ! --- in/out ---------------------------
- real, intent(in) :: xxxz, yyyz
- ! --- begin ----------------------------
- !
- !jqif = nint( 0.5 - sign( 0.5, -dim(xxxz,yyyz)) )
- !
- jqif = nint( 0.5 + sign( 0.5, xxxz-yyyz ) )
- end function jqif
- !=====================================================================================================
- !=====================================================================================================
-
- subroutine read_diffusion( status )
- !
- ! Read vertical diffusion from file and store in conv_dat%dkg
- !
- ! Output:
- ! o status 0 = ok
- ! -1 = file not available for at least one region
- ! 1 = read error
- !
- ! 19 Dec 2012 - P. Le Sager - TM5-MP version, using MDF
- ! --- modules ------------------------------
- use dims, only : itau, revert, region_name, lm, okdebug
- use global_data, only : conv_dat
- use datetime, only : tau2date
- use MDF, only : MDF_OPEN, MDF_HDF4, MDF_Inq_VarID, MDF_Get_Var, MDF_Close, MDF_READ
- use Partools, only : isRoot
- use meteodata, only : global_lli
- use tm5_distgrid, only : dgrid, scatter
- ! --- in/out ----------------------------------------------
- integer, intent(out) :: status
- ! --- const -----------------------------------------------
- character(len=*), parameter :: rname = mname//', read_diffusion'
- ! --- local -----------------------------------------------
- character(len=1000) :: fname
- integer :: region, fid, varid, imr, jmr, sz(3)
- integer, dimension(6) :: idater
- logical :: exist
- real, allocatable :: global_arr(:,:,:)
- ! --- begin -----------------------------------------------
- if ( revert == 1 ) then
- call tau2date( itau, idater )
- else
- ! read field from 3 hours before
- call tau2date( itau - 3600*3, idater )
- end if
- do region = 1, nregions
- ! entire region grid size
- imr = global_lli(region)%nlon
- jmr = global_lli(region)%nlat
- ! for gathering
- sz = shape(conv_dat(region)%dkg) ! to get 3rd dimension
- if(isRoot)then
- allocate(global_arr(imr,jmr,sz(3)))
- else
- allocate(global_arr(1,1,1))
- end if
- if (isRoot) then
-
- ! Create file name
- write(fname,'(2a,i4.4,a,i2.2,a,i2.2,3a,i2,a,i4,5i2.2,a)') trim(diffusion_dir), &
- '/',idater(1),'/',idater(2),'/dkg/',idater(3),'/dkg_', &
- trim(region_name(region)), '_z', lm(region), '_', idater, '.hdf'
- inquire( file=fname, exist=exist )
-
- if ( .not. exist ) then
- write(*,'(a,": dkg file not found for region",i2)') rname, region
- status = -1
- TRACEBACK; return
- end if
- ! Open file
- ! This is a bit too verbose, leading to long output files.
- ! write(*,'(a,": reading ", a)') rname, trim(fname)
-
- call MDF_OPEN(TRIM(fname), MDF_HDF4, MDF_READ, fid, status )
- IF_NOTOK_RETURN(status=1)
- ! Read dkg
- CALL MDF_Inq_VarID( fid, 'dkg', varid, status )
- IF_ERROR_RETURN(status=1)
- CALL MDF_Get_Var( fid, varid, global_arr, status )
- IF_NOTOK_RETURN(status=1)
- end if
- call scatter( dgrid(region), conv_dat(region)%dkg, global_arr, 0, status)
- IF_NOTOK_RETURN(status=1)
-
- if(isRoot) then
- ! Read blh
- CALL MDF_Inq_VarID( fid, 'blh', varid, status )
- IF_ERROR_RETURN(status=1)
- CALL MDF_Get_Var( fid, varid, global_arr(:,:,1), status )
- IF_NOTOK_RETURN(status=1)
- end if
- call scatter( dgrid(region), conv_dat(region)%blh, global_arr(:,:,1), 0, status)
- IF_NOTOK_RETURN(status=1)
-
- ! Close file
- if (isRoot) then
- call MDF_Close( fid, status )
- IF_NOTOK_RETURN(status=1)
- end if
- deallocate(global_arr)
-
- end do
- status = 0
- end subroutine read_diffusion
- !=====================================================================================================
- !=====================================================================================================
- subroutine write_diffusion( status )
- !
- ! Write vertical diffusion to files (for all regions)
- !
- ! Output:
- ! o status 0 = ok, else not ok
- !
- ! 19 Dec 2012 - P. Le Sager - made a TM5-MP version , that works using MDF
-
- ! --- modules ------------------------------
- use dims, only : itau, revert, region_name, lm, okdebug
- use global_data, only : conv_dat
- use datetime, only : tau2date
- use MDF, only : MDF_CREATE, MDF_HDF4, MDF_Def_Dim, MDF_Def_Var, MDF_Close
- use MDF, only : MDF_DOUBLE, MDF_WRITE, MDF_Put_Var, MDF_EndDef
- use Partools, only : isRoot
- use meteodata, only : global_lli
- use tm5_distgrid, only : dgrid, gather
- ! --- in/out ----------------------------------------------
- integer, intent(out) :: status
- ! --- const -----------------------------------------------
- character(len=*), parameter :: rname = mname//'write_diffusion'
- ! --- local -----------------------------------------------
- character(len=1000) :: fname, targetdir
- integer :: region, fid, dkgid, blhid, dimlon, dimlat, dimlev, imr, jmr, sz(3)
- integer, dimension(6) :: idater
- real, allocatable :: global_arr(:,:,:)
- ! --- begin -----------------------------------------------
- if ( revert == 1 ) then
- call tau2date( itau, idater )
- else
- ! read field from 3 hours before
- call tau2date( itau - 3600*3, idater )
- end if
- status=0
- do region = 1, nregions
- ! entire region grid size
- imr = global_lli(region)%nlon
- jmr = global_lli(region)%nlat
-
- ! for gathering
- sz = shape(conv_dat(region)%dkg) ! to get 3rd dimension
-
- if(isRoot)then
- allocate(global_arr(imr,jmr,sz(3)))
- else
- allocate(global_arr(1,1,1))
- end if
- ! Create and define
- if(isRoot) then
-
- ! Create file name
- write(fname,'(2a,i4.4,a,i2.2,a,i2.2,3a,i2,a,i4,5i2.2,a)') trim(diffusion_dir), &
- '/',idater(1),'/',idater(2),'/dkg/',idater(3),'/dkg_', &
- trim(region_name(region)), '_z', lm(region), '_', idater, '.hdf'
- ! Create directory if necessary
- write(targetdir,'(2a,i4.4,a,i2.2,a,i2.2)') trim(diffusion_dir), &
- '/',idater(1),'/',idater(2),'/dkg/',idater(3)
- call system('mkdir -p ' // targetdir)
- if(okdebug) then
- write(gol,'(a,": creating ", a)') rname, trim(fname); call goPr
- end if
- call MDF_CREATE(TRIM(fname), MDF_HDF4, MDF_WRITE, fid, status )
- IF_NOTOK_RETURN(status=1)
-
- call MDF_Def_Dim( fid, 'lon', imr, dimlon, status )
- call MDF_Def_Dim( fid, 'lat', jmr, dimlat, status )
- call MDF_Def_Dim( fid, 'levels', sz(3), dimlev, status )
-
- call MDF_Def_Var( fid, 'dkg', MDF_DOUBLE, (/dimlon,dimlat,dimlev/), dkgid, status)
- call MDF_Def_Var( fid, 'blh', MDF_DOUBLE, (/dimlon,dimlat/), blhid, status)
-
- call MDF_EndDef(fid, status)
- end if
- ! gather and write dkg/blh
- call gather( dgrid(region), conv_dat(region)%dkg, global_arr, 0, status)
- IF_NOTOK_RETURN(status=1)
- if (isRoot) then
- call MDF_Put_Var( fid, dkgid, global_arr, status)
- IF_NOTOK_RETURN(status=1)
- end if
- call gather( dgrid(region), conv_dat(region)%blh, global_arr(:,:,1), 0, status)
- IF_NOTOK_RETURN(status=1)
- if(isRoot) then
- call MDF_Put_Var( fid, blhid, global_arr(:,:,1), status)
- IF_NOTOK_RETURN(status=1)
- ! Close file
- call MDF_CLOSE( fid, status )
- IF_NOTOK_RETURN(status=1)
- end if
- deallocate(global_arr)
-
- end do
- status = 0
- end subroutine write_diffusion
- !=====================================================================================================
- !=====================================================================================================
- END MODULE DIFFUSION
|