!################################################################# !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 minum value for u* to trap division by zero (MK). ! !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 ! --- in/out ------------------------------------ private public :: Diffusion_Init, Diffusion_Done public :: Calc_Kzz, Read_Diffusion, Write_Diffusion 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 Meteo , only : Set use Meteo , only : pu_dat, pv_dat use Meteo , only : temper_dat, humid_dat, gph_dat use Meteo , only : oro_dat, lsmask_dat use Meteo , only : sr_ecm_dat, sr_ols_dat use Meteo , only : u10m_dat, v10m_dat use Meteo , only : slhf_dat, sshf_dat use Meteo , only : ewss_dat, nsss_dat use GO , only : TrcFile, Init, Done, ReadRc use toolbox , only : escape_tm ! --- 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 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: isr,jsr,ier,jer 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 meteo , only: spm_dat use meteo , only: pu_dat, pv_dat use meteo , only: temper_dat, humid_dat, gph_dat use meteodata , only : m_dat #ifdef MPI use mpi_const , only: myid, root_t, ierr, com_trac, my_real, lmloc #endif ! --- in/out ---------------------------- ! --- 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(:,:,:) ! --- begin ---------------------------------- if (okdebug) write(*,*) 'start of calc_kzz' #ifdef MPI if(myid == root_t) then ! reading /calculations by root! #endif do region = 1, nregions allocate(ustar_loc(region)%surf(im(region),jm(region))) allocate(sr_mix(region)%surf(im(region),jm(region))) 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 do region = 1,nregions ! local grid sizes imr = im(region) ; jmr = jm(region) ; lmr = lm(region) 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(1:imr,1:jmr,1:lmr) ) allocate ( phlb_adj(imr,jmr,lmr+1)) allocate ( p_adj(imr,jmr)) am => wind_dat(region)%am_t bm => wind_dat(region)%bm_t cm => wind_dat(region)%cm_t dxyp => region_dat(region)%dxyp ntimes = (nread/(ndyn))*tref(region) m_adj = m_dat(region)%data(1:imr,1:jmr,1:lmr) do i=1,ntimes m_adj = m_adj + revert*am(0:imr-1,1:jmr ,1:lmr ) & !east - revert*am(1:imr ,1:jmr ,1:lmr ) & !west + revert*bm(1:imr ,1:jmr ,1:lmr ) & !south - revert*bm(1:imr ,2:jmr+1,1:lmr ) & !north + revert*cm(1:imr ,1:jmr ,0:lmr-1) & !lower - revert*cm(1:imr ,1:jmr ,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=jsr(region), jer(region) do i=isr(region), ier(region) 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=jsr(region), jer(region) do i=isr(region), ier(region) 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(1:imr,1:jmr,1:lmr+1) ) do l = 1, lmr+1 phlb_t(:,:,l) = at(l) + bt(l)*spm_dat(region)%data(1:imr,1:jmr,1) end do allocate( m_t(-1:imr+2,-1:jmr+2,1:lmr) ) m_t = 0.0 do l = 1, lmr do j = 1, jmr m_t(1:imr,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) #ifndef without_convection #ifndef without_diffusion ! calculate kzz call bldiff( region, phlb, gph, t, q, m, pu, pv, dkg ) #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 #ifdef MPI end if !myid == root_t #endif #ifdef MPI #ifndef without_convection #ifndef without_diffusion do region = 1, nregions nsend = im(region)*jm(region)*lmax_conv call mpi_bcast(conv_dat(region)%dkg ,nsend, my_real, root_t, com_trac, ierr) end do #endif #endif if ( myid == root_t ) then #endif if (okdebug) print *, 'vertical diffusion kzz (dkg) broadcasted' do region = 1,nregions deallocate(sr_mix(region)%surf) deallocate(ustar_loc(region)%surf) end do #ifdef MPI endif #endif 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). ! !------------------------ subroutine dd_calc_ustar(region) use binas, only : grav, vKarman use dims, only : im, jm use dims, only : okdebug use meteo, only : lsmask_dat, sr_ecm_dat, sr_ols_dat use meteo, only : u10m_dat, v10m_dat use meteo, only : slhf_dat, sshf_dat use meteo, only : ewss_dat, nsss_dat !--- in/out ------------------------------------ integer,intent(in) :: region ! --- 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 real, allocatable :: field(:,:) ! --- begin ---------------------------------- allocate(field(im(region),jm(region))) do j=1,jm(region) do i=1,im(region) 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 sr_land=max(1e-2,sr_ols_dat(region)%data(i,j,1)) ! occurs at Islands, etc. 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) 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 field = ustar_loc(region)%surf if (okdebug) call dd_field_statistics('ustar_loc',field,im(region),jm(region)) !if (okdebug) call dumpfield(0,'ustar_loc'//c_time//'.hdf',field,'ustar_loc') deallocate(field) 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, gph, t, q, m, pu, pv, dkg ) 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 use omp_partools,only : my_omp_domain ! --- in/out ---------------------------------------------- integer,intent(in) :: region 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(jm(region)) :: 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 ! Openmp parameters integer :: js, je, jss, jee ! --- begin ------------------------------------------ if (okdebug) write(*,*) 'start of bldiff' allocate(wkvf(im(region),jm(region),lm(region))) allocate(zdup2(im(region),jm(region),lm(region))) allocate(zrinub(im(region),jm(region),lm(region))) allocate(wthm(im(region),jm(region),lm(region))) allocate(zthvk(im(region),jm(region),lm(region))) allocate(pf(im(region),jm(region),lm(region))) allocate(whm(im(region),jm(region),lm(region))) allocate(uwind(im(region),jm(region),lm(region))) allocate(vwind(im(region),jm(region),lm(region))) allocate(kvh(im(region),jm(region),lm(region))) allocate(ustr(im(region),jm(region))) allocate(wheat(im(region),jm(region))) allocate(wqflx(im(region),jm(region))) allocate(wobkl(im(region),jm(region))) allocate(wheatv(im(region),jm(region))) allocate(wtseff(im(region),jm(region))) pblh => conv_dat(region)%blh ! mid-layer pressure levels ! potential temperature field ! virtual potential temperature field ! !$OMP PARALLEL & !$OMP default (none) & !$OMP shared (jsr, jer, isr, ier, region, pf, phlb, t, q, wthm, & !$OMP zthvk, whm, gph) & !$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=jss, jee do i=isr(region), ier(region) 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 !$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 ! !$OMP PARALLEL & !$OMP default (none) & !$OMP shared (region, wheat, sshf_dat, gph, phlb, slhf_dat, & !$OMP wqflx, ustr, ustar_loc) & !$OMP private (i, j, js, je) call my_omp_domain (1, jm(region), js, je) do j = js, je do i = 1, im(region) 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 !$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 = 1,jm(region) lat(j) = ybeg(region) + 0.5 * yres + yres * (j-1) enddo dxx(:) = ae * xres * gtor * cos(lat(:)*gtor) dyy = ae * yres * gtor !$OMP PARALLEL & !$OMP default (none) & !$OMP shared (jsr, jer, isr, ier, region, dxx, dyy, pu, pv, m, & !$OMP uwind, vwind) & !$OMP private (i, j, l, jss, jee) call my_omp_domain (jsr(region), jer(region), jss, jee) do l=1,lm(region) do j=jss,jee do i=isr(region), ier(region) uwind(i,j,l) = dxx(j)*(pu(i,j,l) + pu(i-1,j,l))*0.5 / m(i,j,l) ! cmk Bug repared, 06/2005 yres --> dyy vwind(i,j,l) = dyy* (pv(i,j,l) + pv(i,j+1,l))*0.5 / m(i,j,l) enddo enddo enddo !$OMP END PARALLEL ! ! Initialize gradient Richardson number and free tropospheric k values ! if (okdebug) write(*,*) ' Richardson number' !$OMP PARALLEL & !$OMP default (none) & !$OMP shared (jsr, jer, isr, ier, region, gph, zdup2, uwind, & !$OMP vwind, whm, phlb, pf, zthvk, zrinub, wkvf) & !$OMP private (i, j, l, jss, jee, jq, zlambdac, cml2, zdz, arg, & !$OMP zthva, zsstab, zfunst, zfstabh, zkvn) call my_omp_domain (jsr(region), jer(region), jss, jee) do l=1,lm(region)-1 do j=jss,jee do i=isr(region), ier(region) ! ! 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 !$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' !$OMP PARALLEL & !$OMP default (none) & !$OMP shared (jsr, jer, isr, ier, region, wtseff, wthm, q, & !$OMP wheatv, wheat, wqflx, wobkl, ustr, pblh, uwind, vwind, & !$OMP whm ) & !$OMP private (i, j, l, jss, jee, jwork, wrino, vvk, tkv, dtkv, & !$OMP zrino, jq, fmt, wsc, zexcess, thvparc ) call my_omp_domain (jsr(region), jer(region), jss, jee) do j=jss,jee do i=isr(region),ier(region) ! 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 !$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 ! !$OMP PARALLEL & !$OMP default (none) & !$OMP shared (jsr, jer, isr, ier, region, pblh, gph, wobkl, & !$OMP wheatv, whm, zrinub, zdup2, kvhmin, wkvf, ustr, & !$OMP wtseff, kvh, zthvk) & !$OMP private (i, j, l, jss, jee, jq, z, zh, zl, zzh, jq1, jq2, jq3, & !$OMP jq4, jq5, jq6, jq7, jqq, cml2, zfstabh, zpblk, zkvh, & !$OMP zslask, fmt, wsc, zwstr, obukov, term1, term2, term3, & !$OMP pr, thvgrad, kvhentr) call my_omp_domain (jsr(region), jer(region), jss, jee) do l=1,lm(region)-1 do j=jss,jee do i=isr(region), ier(region) 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 !$OMP END PARALLEl ! calculate vertical diffusion !$OMP PARALLEL & !$OMP default (none) & !$OMP shared (jsr, jer, isr, ier, region, dkg, kvh, m, gph) & !$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=jss,jee do i=isr(region), ier(region) 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(:,jss:jee,lmax_conv) = 0.0 !$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(*,*) 'start of bldiff' end subroutine bldiff #endif #endif ! end subroutine calc_kzz ! *** !----------------------------------------------- ! ! Purpose: ! ------- ! this subroutine calculate the min,mean,max of a real field ! ! External ! -------- ! none ! ! Reference ! --------- ! none !------------------------------------ subroutine dd_field_statistics(name,field,ix,jx) ! --- in/out ------------------------------- character(len=*),intent(in) :: name integer,intent(in) :: ix,jx real,dimension(ix,jx),intent(in) :: field ! --- local -------------------------------- integer :: i,j,ntel_non_zero real :: maxf,minf,meanf,mean_non_zero ! --- begin ------------------------------- maxf=-1e20 minf=1e12 meanf=0. mean_non_zero=0. ntel_non_zero=0 do i=1,ix do j=1,jx meanf=meanf+field(i,j) maxf=max(maxf,field(i,j)) minf=min(minf,field(i,j)) if (field(i,j).ne.0) then ntel_non_zero=ntel_non_zero+1 mean_non_zero=mean_non_zero+field(i,j) endif enddo enddo meanf=meanf/ix/jx if (ntel_non_zero.gt.0) mean_non_zero= mean_non_zero/ntel_non_zero write(6,'(a10,4(a3,1x,1pe10.3))') name,'max',maxf,'min',minf,'mean',meanf,'mn0',mean_non_zero end subroutine dd_field_statistics ! **** !------------------------------------ ! ! Purpose: ! ------- ! this subroutine calculate the min,mean,max of a field of int8 ! ! External ! -------- ! none ! ! Reference ! --------- ! none !------------------------------------ subroutine dd_field_statistics_i8(name,field,ix,jx) ! --- in/out ---------------------------- character(len=*),intent(in) :: name integer,intent(in) :: ix,jx integer(kind=1),dimension(ix,jx),intent(in) :: field ! --- local ---------------------------- integer :: i,j,ntel_non_zero real :: maxf,minf,meanf,mean_non_zero ! --- begin --------------------------- maxf=-1e20 minf=1e12 meanf=0. mean_non_zero=0. ntel_non_zero=0 do i=1,ix do j=1,jx meanf=meanf+field(i,j) maxf=max(maxf,real(field(i,j))) minf=min(minf,real(field(i,j))) if (field(i,j).ne.0) then ntel_non_zero=ntel_non_zero+1 mean_non_zero=mean_non_zero+field(i,j) endif enddo enddo meanf=meanf/ix/jx if (ntel_non_zero.gt.0) mean_non_zero= mean_non_zero/ntel_non_zero write(6,'(a10,4(a3,1x,1pe10.3))') name,'max',maxf,'min',minf,'mean',meanf,'mn0',mean_non_zero end subroutine dd_field_statistics_i8 ! *** ! ! 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 ! ! --- modules ------------------------------ use dims, only : itau, revert, region_name, lm use global_data, only : conv_dat use datetime, only : tau2date use file_hdf, only : Init, ReadData, Done, THdfFile, TSds use toolbox, only : escape_tm ! --- in/out ---------------------------------------------- integer, intent(out) :: status ! --- const ----------------------------------------------- character(len=*), parameter :: rname = mname//', read_diffusion' ! --- local ----------------------------------------------- character(len=1000) :: fname integer :: region integer, dimension(6) :: idater logical :: exist type(THdfFile) :: hdf type(TSds) :: sds ! --- 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 ! 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 return end if ! Open file ! This is a bit too verbose, leading to long output files. ! write(*,'(a,": reading ", a)') rname, trim(fname) call Init( hdf, fname, 'read', status ) IF_NOTOK_RETURN(status=1) ! Read dkg call Init( sds, hdf, 'dkg', status ) IF_NOTOK_RETURN(status=1) call ReadData( sds, conv_dat(region)%dkg ,status) call Done( sds, status ) IF_NOTOK_RETURN(status=1) ! Read blh call Init( sds, hdf, 'blh', status ) IF_NOTOK_RETURN(status=1) call ReadData( sds, conv_dat(region)%blh , status) call Done( sds, status ) IF_NOTOK_RETURN(status=1) ! Close file call Done( hdf, status ) IF_NOTOK_RETURN(status=1) 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 ! ! --- modules ------------------------------ use dims, only : itau, revert, region_name, lm use global_data, only : conv_dat use datetime, only : tau2date use file_hdf, only : Init, WriteData, Done, THdfFile, TSds use toolbox , only : escape_tm use Partools, only : myid, root ! --- in/out ---------------------------------------------- integer, intent(out) :: status ! --- const ----------------------------------------------- character(len=*), parameter :: rname = mname//', write_diffusion' ! --- local ----------------------------------------------- character(len=1000) :: fname,targetdir integer :: region integer, dimension(6) :: idater type(THdfFile) :: hdf type(TSds) :: sds ! --- 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 if (myid /= root) return !WP! only root PE writes fields do region = 1, nregions ! 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) ! Create file write(*,'(a,": creating ", a)') rname, trim(fname) call Init( hdf, fname, 'create', status ) IF_NOTOK_RETURN(status=1) ! Write dkg call Init( sds, hdf, 'dkg', shape(conv_dat(region)%dkg), 'real(8)', status ) IF_NOTOK_RETURN(status=1) call WriteData( sds, conv_dat(region)%dkg ,status) call Done( sds, status ) IF_NOTOK_RETURN(status=1) ! Write blh call Init( sds, hdf, 'blh', shape(conv_dat(region)%blh), 'real(8)', status ) IF_NOTOK_RETURN(status=1) call WriteData( sds, conv_dat(region)%blh,status ) call Done( sds, status ) IF_NOTOK_RETURN(status=1) ! Close file call Done( hdf, status ) IF_NOTOK_RETURN(status=1) end do status = 0 end subroutine write_diffusion !===================================================================================================== !===================================================================================================== end module diffusion