#include #include module vertical_diffusion !--------------------------------------------------------------------------------- ! Module to compute vertical diffusion of momentum, moisture, trace constituents ! and temperature. Uses turbulence module to get eddy diffusivities, including boundary ! layer diffusivities and nonlocal transport terms. Computes and applies molecular ! diffusion, if model top is high enough (above ~90 km). ! ! calling sequence: ! ! vd_inti initializes vertical diffustion constants ! trbinti initializes turblence/pbl constants ! . ! . ! vd_intr interface for vertical diffusion and pbl scheme ! vdiff performs vert diff and pbl ! trbintr compute turbulent diffusivities and boundary layer cg terms ! vd_add_cg add boudary layer cg term to profiles ! vd_lu_decomp perform lu decomposition of vertical diffusion matrices ! vd_lu_qmolec update decomposition for molecular diffusivity of a constituent ! vd_lu_solve solve the euler backward problem for vertical diffusion ! !---------------------------Code history-------------------------------- ! Standardized: J. Rosinski, June 1992 ! Reviewed: P. Rasch, B. Boville, August 1992 ! Reviewed: P. Rasch, April 1996 ! Reviewed: B. Boville, April 1996 ! rewritten: B. Boville, May 2000 ! more changes B. Boville, May-Aug 2001 !--------------------------------------------------------------------------------- use shr_kind_mod, only: r8 => shr_kind_r8 use ppgrid, only: pcols, pver, pverp use constituents, only: pcnst, pnats, cnst_name, cnst_add, qmincg use constituents, only: cnst_get_ind, cnst_get_type_byind use pmgrid, only: masterproc use physconst, only: karman implicit none private ! Make default type private to the module save ! ! Public interfaces ! public vd_inti ! Initialization public vd_intr ! Full routine ! ! Private data ! real(r8), private :: cpair ! Specific heat of dry air real(r8), private :: gravit ! Acceleration due to gravity real(r8), private :: rair ! Gas const for dry air real(r8), private :: zvir ! rh2o/rair - 1 real(r8), parameter :: km_fac = 3.55E-7 ! molecular viscosity constant real(r8), parameter :: pr_num = 1. ! Prandtl number real(r8), parameter :: d0 = 1.52E20 ! diffusion factor m-1 s-1 molec sqrt(kg/kmol/K) ! Aerononmy, Part B, Banks and Kockarts (1973), p39 ! note text cites 1.52E18 cm-1 ... real(r8), parameter :: n_avog = 6.022E26 ! Avogadro's number (molec/kmol) real(r8), parameter :: kq_fac = 2.52E-13 ! molecular constituent diffusivity constant real(r8), parameter :: mw_dry = 29. ! molecular weight of dry air ****REMOVE THIS*** real(r8) :: mw_fac(pcnst+pnats) ! sqrt(1/M_q + 1/M_d) in constituent diffusivity ! turbulent mountain stress constants real(r8), parameter :: z0fac = 0.025 ! factor determining z_0 from orographic standard deviation real(r8), parameter :: z0max = 100. ! max value of z_0 for orography real(r8), parameter :: horomin= 10. ! min value of subgrid orographic height for mountain stress real(r8), parameter :: dv2min = 0.01 ! minimum shear squared real(r8), private :: oroconst ! converts from standard deviation to height integer, private :: ntop ! Top level to which vertical diffusion is applied (=1). integer, private :: nbot ! Bottom level to which vertical diffusion is applied (=pver). integer, private :: ntop_eddy ! Top level to which eddy vertical diffusion is applied. integer, private :: nbot_eddy ! Bottom level to which eddy vertical diffusion is applied (=pver). integer, private :: ntop_molec ! Top level to which molecular vertical diffusion is applied (=1). integer, private :: nbot_molec ! Bottom level to which molecular vertical diffusion is applied. character(len=8), private :: vdiffnam(pcnst+pnats) ! names of v-diff tendencies contains !=============================================================================== subroutine vd_inti(cpairx ,cpwvx ,gravx ,rairx, zvirx, & hypm ,vkx) !----------------------------------------------------------------------- ! Initialization of time independent fields for vertical diffusion. ! Calls initialization routine for boundary layer scheme. !----------------------------------------------------------------------- use history, only: addfld, add_default, phys_decomp use turbulence, only: trbinti use dycore, only: dycore_is !------------------------------Arguments-------------------------------- real(r8), intent(in) :: cpairx ! specific heat of dry air real(r8), intent(in) :: cpwvx ! spec. heat of water vapor at const. pressure real(r8), intent(in) :: gravx ! acceleration due to gravity real(r8), intent(in) :: rairx ! gas constant for dry air real(r8), intent(in) :: zvirx ! rh2o/rair - 1 real(r8), intent(in) :: hypm(pver) ! reference pressures at midpoints real(r8), intent(in) :: vkx ! Von Karman's constant !---------------------------Local workspace----------------------------- integer :: k ! vertical loop index !----------------------------------------------------------------------- ! Set physical constants for vertical diffusion and pbl: cpair = cpairx gravit = gravx rair = rairx zvir = zvirx ! Range of levels to which v-diff is applied. ntop_eddy = 1 ! no reason not to make this 1, if >1, must be <= nbot_molec nbot_eddy = pver ! should always be pver ntop_molec = 1 ! should always be 1 nbot_molec = 0 ! should be set below about 70 km ! Molecular diffusion turned on above ~60 km (50 Pa) if model top is above ~90 km (.1 Pa). ! Note that computing molecular diffusivities is a trivial expense, but constituent ! diffusivities depend on their molecular weights. Decomposing the diffusion matric ! for each constituent is a needless expense unless the diffusivity is significant. if (hypm(1) .lt. 0.1) then do k = 1, pver if (hypm(k) .lt. 50.) nbot_molec = k end do if (masterproc) then write (6,fmt='(a,i3,5x,a,i3)') 'NTOP_MOLEC =',ntop_molec, 'NBOT_MOLEC =',nbot_molec write (6,fmt='(a,i3,5x,a,i3)') 'NTOP_EDDY =',ntop_eddy, 'NBOT_EDDY =',nbot_eddy endif end if ! The vertical diffusion solver must operate over the full range of molecular and eddy diffusion ntop = 1 nbot = pver ! Molecular weight factor in constitutent diffusivity ! ***** FAKE THIS FOR NOW USING MOLECULAR WEIGHT OF DRY AIR FOR ALL TRACERS **** do k = 1, pcnst+pnats mw_fac(k) = d0 * mw_dry * sqrt(1./mw_dry + 1./mw_dry) / n_avog end do ! Initialize turbulence variables call trbinti(gravit, cpair, rair, zvir, ntop_eddy, nbot_eddy, hypm, vkx) ! Set names of diffused variable tendencies and declare them as history variables do k = 1, pcnst+pnats vdiffnam(k) = 'VD'//cnst_name(k) if(k==1)vdiffnam(k) = 'VD01' !**** compatibility with old code **** call addfld (vdiffnam(k),'kg/kg/s ',pver, 'A','Vertical diffusion of '//cnst_name(k),phys_decomp) end do ! Only tracer 1 (Q) is output by default call add_default (vdiffnam(1), 1, ' ') ! Turbulent mountain stress constant ! This is NOT tunable - just a convenient way to turn mountain stress off. if ( dycore_is ('LR') ) then oroconst = 0.0 ! Restore this to work with turbulent mountain stress ! oroconst = 4.0 else oroconst = 0.0 endif ! Turbulent mountain stress terms call addfld ('TAUTMSX' ,'N/m2 ',1, 'A','Zonal turbulent mountain surface stress' , phys_decomp) call addfld ('TAUTMSY' ,'N/m2 ',1, 'A','Meridional turbulent mountain surface stress', phys_decomp) ! call add_default ('TAUTMSX ', 1, ' ') ! call add_default ('TAUTMSY ', 1, ' ') return end subroutine vd_inti !=============================================================================== subroutine vd_intr( & ztodt ,state , & taux ,tauy ,shflx ,cflx ,pblh , & tpert ,qpert ,ustar ,obklen ,ptend , & cldn ,ocnfrac ,landfrac ,sgh ) !----------------------------------------------------------------------- ! interface routine for vertical diffusion and pbl scheme !----------------------------------------------------------------------- use physics_types, only: physics_state, physics_ptend use history, only: outfld !!$ use geopotential, only: geopotential_dse !------------------------------Arguments-------------------------------- real(r8), intent(in) :: taux(pcols) ! x surface stress (N/m2) real(r8), intent(in) :: tauy(pcols) ! y surface stress (N/m2) real(r8), intent(in) :: shflx(pcols) ! surface sensible heat flux (w/m2) real(r8), intent(in) :: cflx(pcols,pcnst+pnats)! surface constituent flux (kg/m2/s) real(r8), intent(in) :: ztodt ! 2 delta-t real(r8), intent(in) :: cldn(pcols,pver) ! new cloud fraction real(r8), intent(in) :: ocnfrac(pcols) ! Ocean fraction real(r8), intent(in) :: landfrac(pcols) ! Land fraction real(r8), intent(in) :: sgh(pcols) ! standard deviation of orography type(physics_state), intent(in) :: state ! Physics state variables type(physics_ptend), intent(inout) :: ptend ! indivdual parameterization tendencies real(r8), intent(out) :: pblh(pcols) ! planetary boundary layer height real(r8), intent(out) :: tpert(pcols) ! convective temperature excess real(r8), intent(out) :: qpert(pcols,pcnst+pnats) ! convective humidity and constituent excess real(r8), intent(out) :: ustar(pcols) ! surface friction velocity real(r8), intent(out) :: obklen(pcols) ! Obukhov length ! !---------------------------Local storage------------------------------- integer :: lchnk ! chunk identifier integer :: ncol ! number of atmospheric columns integer :: i,k,m ! longitude,level,constituent indices integer :: ixcldice, ixcldliq ! constituent indices for cloud liquid and ice water real(r8) :: dtk(pcols,pver) ! T tendency from KE dissipation real(r8) :: kvh(pcols,pverp) ! diffusion coefficient for heat real(r8) :: kvm(pcols,pverp) ! diffusion coefficient for momentum real(r8) :: cgs(pcols,pverp) ! counter-gradient star (cg/flux) real(r8) :: rztodt ! 1./ztodt real(r8) :: tautmsx(pcols) ! turbulent mountain surface stress (zonal) real(r8) :: tautmsy(pcols) ! turbulent mountain surface stress (meridional) !++pjr ! array used to calculate tend wrt dry pressure type(physics_ptend) ptendx ! indivdual parameterization tendencies real(r8) :: dtkx(pcols,pver) ! T tendency from KE dissipation !Trash variables for second call to vdiff real(r8) pblh_trash(pcols) ! planetary boundary layer height real(r8) tpert_trash(pcols) ! convective temperature excess real(r8) qpert_trash(pcols) ! convective humidity and constituent excess real(r8) ustar_trash(pcols) ! surface friction velocity real(r8) obklen_trash(pcols) ! Obukhov length real(r8) tautmsx_trash(pcols) ! turbulent mountain surface stress (zonal) real(r8) tautmsy_trash(pcols) ! turbulent mountain surface stress (meridional) logical :: dryflag ! = .true. if there are dry constituents !--pjr !----------------------------------------------------------------------- ! local constants rztodt = 1./ztodt lchnk = state%lchnk ncol = state%ncol ! Operate on copies of the input states, convert to tendencies at end. ! The input values of u,v are required for computing the kinetic energy dissipation. ptend%q(:ncol,:,:) = state%q(:ncol,:,:) ptend%s(:ncol,:) = state%s(:ncol,:) ptend%u(:ncol,:) = state%u(:ncol,:) ptend%v(:ncol,:) = state%v(:ncol,:) ! All variables are modified by vertical diffusion ptend%name = "vertical diffusion" ptend%lq(:) = .TRUE. ptend%ls = .TRUE. ptend%lu = .TRUE. ptend%lv = .TRUE. ! Call the real vertical diffusion code. call vdiff( & lchnk ,ncol , & ptend%u ,ptend%v ,state%t ,ptend%q ,ptend%s , & state%pmid ,state%pint ,state%lnpint,state%lnpmid,state%pdel , & state%rpdel ,state%zm ,state%zi ,ztodt ,taux , & tauy ,shflx ,cflx ,pblh ,ustar , & kvh ,kvm ,tpert ,qpert(1,1) ,cgs , & obklen ,state%exner ,dtk ,cldn ,ocnfrac , & landfrac ,sgh ,tautmsx ,tautmsy ) ! Convert the new profiles into vertical diffusion tendencies. ! Convert KE dissipative heat change into "temperature" tendency. do k = 1, pver do i = 1, ncol ptend%s(i,k) = (ptend%s(i,k) - state%s(i,k)) * rztodt ptend%u(i,k) = (ptend%u(i,k) - state%u(i,k)) * rztodt ptend%v(i,k) = (ptend%v(i,k) - state%v(i,k)) * rztodt end do do m = 1, pcnst+pnats do i = 1, ncol ptend%q(i,k,m) = (ptend%q(i,k,m) - state%q(i,k,m))*rztodt end do end do end do #ifdef PERGRO ! For pergro case, do not diffuse cloud water: replace with input values call cnst_get_ind('CLDLIQ', ixcldliq) call cnst_get_ind('CLDICE', ixcldice) ptend%lq(ixcldice) = .FALSE. ptend%lq(ixcldliq) = .FALSE. ptend%q(:ncol,:,ixcldice) = 0. ptend%q(:ncol,:,ixcldliq) = 0. #endif !++pjr ! ! kludge to allow calculation wrt dry pressures. This lets ! mixing ratios be wrt dry air rather than moist air ! it is expensive, but easy ! ! Operate on copies of the input states, convert to tendencies at end. ! The input values of u,v are required for computing the kinetic energy dissipation. ! only do this if there are some dry-type constituents dryflag = .false. do m = 1,pcnst+pnats if (cnst_get_type_byind(m).eq.'dry' ) dryflag = .true. end do ! Only do the vertical diffusion for dry if there are dry-type constituents if ( dryflag ) then ptendx%q(:ncol,:,:) = state%q(:ncol,:,:) ptendx%s(:ncol,:) = state%s(:ncol,:) ptendx%u(:ncol,:) = state%u(:ncol,:) ptendx%v(:ncol,:) = state%v(:ncol,:) ! All variables are modified by vertical diffusion ptendx%name = "vertical diffusion dry" ptendx%lq(:) = .TRUE. ptendx%ls = .TRUE. ptendx%lu = .TRUE. ptendx%lv = .TRUE. ! make temporary variables to hold fields that we don't want changed by (second) call to vdiff pblh_trash = pblh ustar_trash = ustar tpert_trash = tpert qpert_trash = qpert(1:pcols,1) obklen_trash = obklen tautmsx_trash = tautmsx tautmsy_trash = tautmsy ! Call the real vertical diffusion code. ! Note that dryflag (=.true.) is an optional argument that prevents vdiff ! (calling trbintd) from writing out trash fields a second time. call vdiff( & lchnk ,ncol , & ptendx%u ,ptendx%v ,state%t ,ptendx%q ,ptendx%s , & state%pmiddry ,state%pintdry ,state%lnpintdry ,state%lnpmiddry,state%pdeldry, & state%rpdeldry,state%zm ,state%zi ,ztodt ,taux , & tauy ,shflx ,cflx ,pblh_trash ,ustar_trash , & kvh ,kvm ,tpert_trash,qpert_trash ,cgs , & obklen_trash ,state%exner ,dtkx ,cldn ,ocnfrac , & landfrac ,sgh ,tautmsx_trash ,tautmsy_trash , dryflag) ! Convert the new profiles into vertical diffusion tendencies. ! Convert KE dissipative heat change into "temperature" tendency. do k = 1, pver do i = 1, ncol ptendx%s(i,k) = (ptendx%s(i,k) - state%s(i,k)) * rztodt ptendx%u(i,k) = (ptendx%u(i,k) - state%u(i,k)) * rztodt ptendx%v(i,k) = (ptendx%v(i,k) - state%v(i,k)) * rztodt end do do m = 1, pcnst+pnats do i = 1, ncol ptendx%q(i,k,m) = (ptendx%q(i,k,m) - state%q(i,k,m))*rztodt end do end do end do ! do m = 1,pcnst+pnats if (cnst_get_type_byind(m).eq.'dry' ) then ptend%q(:ncol,:,m) = ptendx%q(:ncol,:,m) endif end do endif ! if ( do_dry_vdiff ) !--pjr ! Save the vertical diffusion variables on the history file dtk(:ncol,:) = dtk(:ncol,:)/cpair ! normalize heating for history call outfld ('DTVKE ',dtk,pcols,lchnk) dtk(:ncol,:) = ptend%s(:ncol,:)/cpair ! normalize heating for history using dtk call outfld ('DTV ',dtk ,pcols,lchnk) call outfld ('DUV ',ptend%u,pcols,lchnk) call outfld ('DVV ',ptend%v,pcols,lchnk) do m = 1, pcnst+pnats call outfld(vdiffnam(m),ptend%q(1,1,m),pcols,lchnk) end do call outfld ('TAUTMSX', tautmsx, pcols, lchnk) call outfld ('TAUTMSY', tautmsy, pcols, lchnk) return end subroutine vd_intr !=============================================================================== subroutine vdiff (lchnk ,ncol , & u ,v ,t ,q ,dse , & pmid ,pint ,piln ,pmln ,pdel , & rpdel ,zm ,zi ,ztodt ,taux , & tauy ,shflx ,cflx ,pblh ,ustar , & kvh ,kvm ,tpert ,qpert ,cgs , & obklen ,exner ,dtk ,cldn ,ocnfrac , & landfrac ,sgh ,tautmsx ,tautmsy ,dryflag_in ) !----------------------------------------------------------------------- ! Driver routine to compute vertical diffusion of momentum, moisture, trace ! constituents and dry static energy. The new temperature is computed from ! the diffused dry static energy. ! Turbulent diffusivities and boundary layer nonlocal transport terms are ! obtained from the turbulence module. !----------------------------------------------------------------------- use turbulence, only: trbintr !------------------------------Arguments-------------------------------- integer, intent(in) :: lchnk ! chunk identifier integer, intent(in) :: ncol ! number of atmospheric columns real(r8), intent(in) :: exner(pcols,pver) ! (ps/pmid)**cappa real(r8), intent(in) :: pmid(pcols,pver) ! midpoint pressures real(r8), intent(in) :: pint(pcols,pverp) ! interface pressures real(r8), intent(in) :: piln (pcols,pverp) ! Log interface pressures real(r8), intent(in) :: pmln (pcols,pver) ! Log midpoint pressures real(r8), intent(in) :: pdel(pcols,pver) ! thickness between interfaces real(r8), intent(in) :: rpdel(pcols,pver) ! 1./pdel real(r8), intent(in) :: t(pcols,pver) ! temperature input real(r8), intent(in) :: ztodt ! 2 delta-t real(r8), intent(in) :: taux(pcols) ! x surface stress (N/m2) real(r8), intent(in) :: tauy(pcols) ! y surface stress (N/m2) real(r8), intent(in) :: shflx(pcols) ! surface sensible heat flux (W/m2) real(r8), intent(in) :: cflx(pcols,pcnst+pnats)! surface constituent flux (kg/m2/s) real(r8), intent(in) :: cldn(pcols,pver) ! new cloud fraction real(r8), intent(in) :: ocnfrac(pcols) ! Ocean fraction real(r8), intent(in) :: landfrac(pcols) ! Land fraction real(r8), intent(in) :: sgh(pcols) ! standard deviation of orography real(r8), intent(inout) :: u(pcols,pver) ! u wind real(r8), intent(inout) :: v(pcols,pver) ! v wind real(r8), intent(inout) :: q(pcols,pver,pcnst+pnats)! moisture and trace constituents real(r8), intent(inout) :: dse(pcols,pver) ! dry static energy [J/kg] real(r8), intent(in) :: zm(pcols,pver) ! midpoint geoptl height above sfc real(r8), intent(in) :: zi(pcols,pverp) ! interface geoptl height above sfc real(r8), intent(out) :: pblh(pcols) ! planetary boundary layer height real(r8), intent(out) :: ustar(pcols) ! surface friction velocity real(r8), intent(out) :: kvh(pcols,pverp) ! diffusivity for heat real(r8), intent(out) :: kvm(pcols,pverp) ! viscosity (diffusivity for momentum) real(r8), intent(out) :: tpert(pcols) ! convective temperature excess real(r8), intent(out) :: qpert(pcols) ! convective humidity excess real(r8), intent(out) :: cgs(pcols,pverp) ! counter-grad star (cg/flux) real(r8), intent(out) :: obklen(pcols) ! Obukhov length real(r8), intent(out) :: dtk(pcols,pver) ! T tendency from KE dissipation real(r8), intent(out) :: tautmsx(pcols) ! turbulent mountain surface stress (zonal) real(r8), intent(out) :: tautmsy(pcols) ! turbulent mountain surface stress (meridional) logical, intent(in), optional :: dryflag_in ! if vdiff is called just for dry tendencies ! !---------------------------Local workspace----------------------------- real(r8) :: tmpm(pcols,pver) ! potential temperature, ze term in tri-diag sol'n real(r8) :: ca(pcols,pver) ! -upper diag of tri-diag matrix real(r8) :: cc(pcols,pver) ! -lower diag of tri-diag matrix real(r8) :: dnom(pcols,pver) ! 1./(1. + ca(k) + cc(k) - cc(k)*ze(k-1)) real(r8) :: kqfs(pcols,pcnst+pnats) ! kinematic surf constituent flux (kg/m2/s) real(r8) :: kvq(pcols,pverp) ! diffusivity for constituents real(r8) :: kq_scal(pcols,pverp) ! kq_fac*sqrt(T)*m_d/rho for molecular diffusivity real(r8) :: mkvisc ! molecular kinematic viscosity c*tint**(2/3)/rho real(r8) :: tint ! interface temperature real(r8) :: tmp1(pcols) ! temporary storage real(r8) :: tmpi1(pcols,pverp) ! heat cg term [J/kg/m], or interface KE dissipation real(r8) :: tmpi2(pcols,pverp) ! rho or dt*(g*rho)**2/dp at interfaces real(r8) :: zero (pcols) ! zero array for surface heat exchange coefficients real(r8) :: ksrftms(pcols) ! surface "drag" coefficient for mountains real(r8) :: tautotx(pcols) ! total surface stress (zonal) real(r8) :: tautoty(pcols) ! total surface stress (meridional) !!$ real(r8) :: dampx,dampy ! damping constant for applying surface stress real(r8) :: dinp_u(pcols,pverp) ! vertical difference at interfaces, input u real(r8) :: dinp_v(pcols,pverp) ! vertical difference at interfaces, input v real(r8) :: dout_u ! vertical difference at interfaces, output u real(r8) :: dout_v ! vertical difference at interfaces, output v integer :: i,k,m ! longitude,level,constituent indices integer :: ktopbl(pcols) ! index of first midpoint inside pbl integer :: ktopblmn ! min value of ktopbl logical :: dryflag ! send to vdiff if called just for dry tendencies !----------------------------------------------------------------------- zero(:) = 0. ! ! Check if vdiff has just been called to calculate dry tendencies. If so, outfld call ! will not be made ! dryflag = .false. if ( present(dryflag_in) ) then dryflag = dryflag_in endif ! Compute the vertical differences of the input u,v for KE dissipation, interface k- ! Note, velocity=0 at surface, so then difference at the bottom interface is -u,v(pver) do i = 1, ncol dinp_u(i,1) = 0. dinp_v(i,1) = 0. dinp_u(i,pverp) = -u(i,pver) dinp_v(i,pverp) = -v(i,pver) end do do k = 2, pver do i = 1, ncol dinp_u(i,k) = u(i,k) - u(i,k-1) dinp_v(i,k) = v(i,k) - v(i,k-1) end do end do ! compute the turbulent mountain stress call trb_mtn_stress( lchnk , ncol , & u , v , t , pmid , exner , zm , sgh , & ksrftms, tautmsx, tautmsy ) ! add the surface stress and the mountain stress scaled by land fraction do i = 1, ncol ksrftms(i) = landfrac(i) * ksrftms(i) tautmsx(i) = landfrac(i) * tautmsx(i) tautmsy(i) = landfrac(i) * tautmsy(i) tautotx(i) = taux(i) + tautmsx(i) tautoty(i) = tauy(i) + tautmsy(i) end do ! Compute potential temperature tmpm(:ncol,:pver) = t(:ncol,:pver) * exner(:ncol,:pver) ! Get diffusivities and c-g terms from turbulence module call trbintr(lchnk ,ncol , & tmpm ,t ,q ,zm ,zi , & pmid ,u ,v ,tautotx ,tautoty , & shflx ,cflx ,obklen ,ustar ,pblh , & kvm ,kvh ,tmpi1 ,cgs ,kqfs , & tpert ,qpert ,ktopbl ,ktopblmn,cldn , & ocnfrac ,dryflag ) ! Compute rho at interfaces p/RT, Ti = (Tm_k + Tm_k-1)/2, interface k- do k = 2, pver do i = 1, ncol tmpi2(i,k) = pint(i,k) * 2. / (rair*(t(i,k) + t(i,k-1))) end do end do do i = 1, ncol tmpi2(i,pverp) = pint(i,pverp) / (rair*t(i,pver)) end do ! Copy turbulent diffusivity for constituents from heat diffusivity kvq(:ncol,:) = kvh(:ncol,:) ! Compute molecular kinematic viscosity, heat diffusivity and factor for constituent diffusivity ! For the moment, keep molecular diffusivities all the same kq_scal(:ncol,ntop_molec) = 0. do k = ntop_molec+1, nbot_molec do i = 1, ncol tint = 0.5*(t(i,k) + t(i,k-1)) mkvisc = km_fac * tint**(2./3.) / tmpi2(i,k) kvm(i,k) = kvm(i,k) + mkvisc kvh(i,k) = kvh(i,k) + mkvisc*pr_num kq_scal(i,k) = sqrt(tint) / tmpi2(i,k) end do end do ! Call gravity wave drag to get gw tendency and diffusivities ! ************ ADD GW CALL HERE ***************************** ! Add gw momentum forcing and KE dissipation ! ************ ADD CODE HERE ******************************** ! Add the nonlocal transport terms to dry static energy, specific humidity and ! other constituents in the boundary layer. call vd_add_cg( & lchnk ,ncol , & ktopblmn ,q ,dse ,tmpi2 ,kvh , & tmpi1 ,cgs ,kqfs ,rpdel ,ztodt ) ! Add the explicit surface fluxes to the lowest layer (do not include moutain stress) do i = 1, ncol tmp1(i) = ztodt*gravit*rpdel(i,pver) u(i,pver) = u(i,pver) + tmp1(i) * taux(i) v(i,pver) = v(i,pver) + tmp1(i) * tauy(i) !!$ dampx = - gravit * tautotx(i) * rpdel(i,pver) / u(i,pver) !!$ dampy = - gravit * tautoty(i) * rpdel(i,pver) / v(i,pver) !!$ print *, "dampx, dampy", lchnk,i, dampx, dampy !!$ u(i,pver) = u(i,pver) / (1. + ztodt*dampx) !!$ v(i,pver) = v(i,pver) / (1. + ztodt*dampx) dse(i,pver) = dse(i,pver) + tmp1(i) * shflx(i) end do do m = 1, pcnst+pnats do i = 1, ncol q(i,pver,m) = q(i,pver,m) + tmp1(i) * cflx(i,m) end do end do ! Calculate dt * (g*rho)^2/dp at interior interfaces, interface k- do k = 2, pver do i = 1, ncol tmpi2(i,k) = ztodt * (gravit*tmpi2(i,k))**2 / (pmid(i,k) - pmid(i,k-1)) end do end do ! Diffuse momentum call vd_lu_decomp( & lchnk ,ncol , & ksrftms ,kvm ,tmpi2 ,rpdel ,ztodt , & ca ,cc ,dnom ,tmpm ,ntop ,nbot ) call vd_lu_solve( & lchnk ,ncol , & u ,ca ,tmpm ,dnom ,ntop ,nbot ) call vd_lu_solve( & lchnk ,ncol , & v ,ca ,tmpm ,dnom ,ntop ,nbot ) ! Recalculate the mountain and total stresses do i = 1, ncol tautmsx(i) = - ksrftms(i) * u(i,pver) tautmsy(i) = - ksrftms(i) * v(i,pver) tautotx(i) = taux(i) + tautmsx(i) tautoty(i) = tauy(i) + tautmsy(i) end do ! Calculate kinetic energy dissipation ! 1. compute dissipation term at interfaces k = pverp do i = 1, ncol tmpi1(i,1) = 0. tmpi1(i,k) = 0.5 * ztodt * gravit * & ( (-u(i,k-1) + dinp_u(i,k))*tautotx(i) + (-v(i,k-1)+dinp_v(i,k))*tautoty(i) ) end do do k = 2, pver do i = 1, ncol dout_u = u(i,k) - u(i,k-1) dout_v = v(i,k) - v(i,k-1) tmpi1(i,k) = 0.25 * tmpi2(i,k) * kvm(i,k) * & (dout_u**2 + dout_v**2 + dout_u*dinp_u(i,k) + dout_v*dinp_v(i,k)) end do end do ! 2. Compute dissipation term at midpoints, add to dry static energy do k = 1, pver do i = 1, ncol dtk(i,k) = (tmpi1(i,k+1) + tmpi1(i,k)) * rpdel(i,k) dse(i,k) = dse(i,k) + dtk(i,k) end do end do ! Diffuse dry static energy call vd_lu_decomp( & lchnk ,ncol , & zero ,kvh ,tmpi2 ,rpdel ,ztodt , & ca ,cc ,dnom ,tmpm ,ntop ,nbot ) call vd_lu_solve( & lchnk ,ncol , & dse ,ca ,tmpm ,dnom ,ntop ,nbot ) ! Diffuse constituents. ! New decomposition required only if kvq != kvh (gw or molec diffusion) call vd_lu_decomp( & lchnk ,ncol , & zero ,kvq ,tmpi2 ,rpdel ,ztodt , & ca ,cc ,dnom ,tmpm ,ntop_eddy,nbot_eddy) do m = 1, pcnst+pnats ! update decomposition in molecular diffusion range if (ntop_molec < nbot_molec) then call vd_lu_qmolec( & ncol , & kvq ,kq_scal ,mw_fac(m) ,tmpi2 ,rpdel , & ztodt ,ca ,cc ,dnom ,tmpm , & ntop_molec ,nbot_molec ) end if call vd_lu_solve( & lchnk ,ncol , & q(1,1,m),ca ,tmpm ,dnom ,ntop ,nbot ) end do return end subroutine vdiff !============================================================================== subroutine vd_add_cg ( & lchnk ,ncol , & ktopblmn ,q ,dse ,rhoi ,kvh , & cgh ,cgs ,kqfs ,rpdel ,ztodt ) !----------------------------------------------------------------------- ! Add the "counter-gradient" term to the dry static energy and tracer fields ! in the boundary layer. ! Note, ktopblmn gives the minimum vertical index of the first midpoint ! within the boundary layer. ! This is really a boundary layer routine, but, since the pbl also supplies the ! diffusivities in the boundary layer, there is not a clear divsion between ! pbl and eddy vertical diffusion code. !------------------------------Arguments-------------------------------- integer, intent(in) :: lchnk ! chunk identifier integer, intent(in) :: ncol ! number of atmospheric columns integer, intent(in) :: ktopblmn ! min value of ktopbl (highest bl top) real(r8), intent(in) :: rhoi(pcols,pverp) ! density (p/RT) at interfaces real(r8), intent(in) :: kvh(pcols,pverp) ! coefficient for heat and tracers real(r8), intent(in) :: cgh(pcols,pverp) ! countergradient term for heat [J/kg/m] real(r8), intent(in) :: cgs(pcols,pverp) ! unscaled cg term for constituents real(r8), intent(in) :: kqfs(pcols,pcnst+pnats)! kinematic surf constituent flux (kg/m2/s) real(r8), intent(in) :: rpdel(pcols,pver) ! 1./pdel (thickness bet interfaces) real(r8), intent(in) :: ztodt ! 2 delta-t real(r8), intent(inout) :: q(pcols,pver,pcnst+pnats) ! moisture and trace constituent input real(r8), intent(inout) :: dse(pcols,pver) ! dry static energy !---------------------------Local workspace----------------------------- real(r8) :: qtm(pcols,pver) ! temporary trace constituent input logical lqtst(pcols) ! adjust vertical profiles integer i ! longitude index integer k ! vertical index integer m ! constituent index !----------------------------------------------------------------------- ! Add counter-gradient to input static energy profiles do k=ktopblmn-1,pver do i=1,ncol dse(i,k) = dse(i,k) + ztodt*rpdel(i,k)*gravit & *(rhoi(i,k+1)*kvh(i,k+1)*cgh(i,k+1) & - rhoi(i,k )*kvh(i,k )*cgh(i,k )) end do end do ! Add counter-gradient to input mixing ratio profiles. ! Check for neg q's in each constituent and put the original vertical ! profile back if a neg value is found. A neg value implies that the ! quasi-equilibrium conditions assumed for the countergradient term are ! strongly violated. do m = 1, pcnst+pnats qtm(:ncol,ktopblmn-1:pver) = q(:ncol,ktopblmn-1:pver,m) do k = ktopblmn-1, pver do i = 1, ncol q(i,k,m) = q(i,k,m) + ztodt*rpdel(i,k)*gravit*kqfs(i,m) & *(rhoi(i,k+1)*kvh(i,k+1)*cgs(i,k+1) & - rhoi(i,k )*kvh(i,k )*cgs(i,k )) end do end do lqtst(:ncol) = all(q(:ncol,ktopblmn-1:pver,m) >= qmincg(m), 2) do k = ktopblmn-1, pver q(:ncol,k,m) = merge (q(:ncol,k,m), qtm(:ncol,k), lqtst(:ncol)) end do end do return end subroutine vd_add_cg !============================================================================== subroutine vd_lu_decomp( & lchnk ,ncol , & ksrf ,kv ,tmpi ,rpdel ,ztodt , & ca ,cc ,dnom ,ze ,ntop , & nbot ) !----------------------------------------------------------------------- ! Determine superdiagonal (ca(k)) and subdiagonal (cc(k)) coeffs of the ! tridiagonal diffusion matrix. ! The diagonal elements (1+ca(k)+cc(k)) are not required by the solver. ! Also determine ze factor and denominator for ze and zf (see solver). !------------------------------Arguments-------------------------------- integer, intent(in) :: lchnk ! chunk identifier integer, intent(in) :: ncol ! number of atmospheric columns integer, intent(in) :: ntop ! top level to operate on integer, intent(in) :: nbot ! bottom level to operate on real(r8), intent(in) :: ksrf(pcols) ! surface "drag" coefficient real(r8), intent(in) :: kv(pcols,pverp) ! vertical diffusion coefficients real(r8), intent(in) :: tmpi(pcols,pverp) ! dt*(g/R)**2/dp*pi(k+1)/(.5*(tm(k+1)+tm(k))**2 real(r8), intent(in) :: rpdel(pcols,pver) ! 1./pdel (thickness bet interfaces) real(r8), intent(in) :: ztodt ! 2 delta-t real(r8), intent(out) :: ca(pcols,pver) ! -upper diagonal real(r8), intent(out) :: cc(pcols,pver) ! -lower diagonal real(r8), intent(out) :: dnom(pcols,pver) ! 1./(1. + ca(k) + cc(k) - cc(k)*ze(k-1)) ! 1./(b(k) - c(k)*e(k-1)) real(r8), intent(out) :: ze(pcols,pver) ! term in tri-diag. matrix system !---------------------------Local workspace----------------------------- integer :: i ! longitude index integer :: k ! vertical index !----------------------------------------------------------------------- ! Determine superdiagonal (ca(k)) and subdiagonal (cc(k)) coeffs of the ! tridiagonal diffusion matrix. The diagonal elements (cb=1+ca+cc) are ! a combination of ca and cc; they are not required by the solver. do k = nbot-1, ntop, -1 do i = 1, ncol ca(i,k ) = kv(i,k+1)*tmpi(i,k+1)*rpdel(i,k ) cc(i,k+1) = kv(i,k+1)*tmpi(i,k+1)*rpdel(i,k+1) end do end do ! The bottom element of the upper diagonal (ca) is zero (element not used). ! The subdiagonal (cc) is not needed in the solver. do i=1,ncol ca(i,nbot) = 0. end do ! Calculate e(k). This term is ! required in solution of tridiagonal matrix defined by implicit diffusion eqn. do i = 1, ncol dnom(i,nbot) = 1./(1. + cc(i,nbot) + ksrf(i)*ztodt*gravit*rpdel(i,nbot)) ze(i,nbot) = cc(i,nbot)*dnom(i,nbot) end do do k = nbot-1, ntop+1, -1 do i=1,ncol dnom(i,k) = 1./ (1. + ca(i,k) + cc(i,k) - ca(i,k)*ze(i,k+1)) ze(i,k) = cc(i,k)*dnom(i,k) end do end do do i=1,ncol dnom(i,ntop) = 1./ (1. + ca(i,ntop) - ca(i,ntop)*ze(i,ntop+1)) end do return end subroutine vd_lu_decomp !============================================================================== subroutine vd_lu_qmolec( & ncol , & kv ,kq_scal ,mw_facm ,tmpi ,rpdel , & ztodt ,ca ,cc ,dnom ,ze , & ntop ,nbot ) !----------------------------------------------------------------------- ! Add the molecular diffusivity to the turbulent and gw diffusivity for a ! consitutent. ! Update the superdiagonal (ca(k)) and subdiagonal (cc(k)) coeffs of the ! tridiagonal diffusion matrix, also ze and denominator. !------------------------------Arguments-------------------------------- integer, intent(in) :: ncol ! number of atmospheric columns integer, intent(in) :: ntop ! top level to operate on integer, intent(in) :: nbot ! bottom level to operate on real(r8), intent(in) :: kv(pcols,pverp) ! vertical diffusion coefficients real(r8), intent(in) :: kq_scal(pcols,pverp) ! kq_fac*sqrt(T)*m_d/rho for molecular diffusivity real(r8), intent(in) :: mw_facm ! sqrt(1/M_q + 1/M_d) for this constituent real(r8), intent(in) :: tmpi(pcols,pverp) ! dt*(g/R)**2/dp*pi(k+1)/(.5*(tm(k+1)+tm(k))**2 real(r8), intent(in) :: rpdel(pcols,pver) ! 1./pdel (thickness bet interfaces) real(r8), intent(in) :: ztodt ! 2 delta-t real(r8), intent(inout) :: ca(pcols,pver) ! -upper diagonal real(r8), intent(inout) :: cc(pcols,pver) ! -lower diagonal real(r8), intent(inout) :: dnom(pcols,pver) ! 1./(1. + ca(k) + cc(k) - cc(k)*ze(k-1)) ! 1./(b(k) - c(k)*e(k-1)) real(r8), intent(inout) :: ze(pcols,pver) ! term in tri-diag. matrix system !---------------------------Local workspace----------------------------- integer :: i ! longitude index integer :: k ! vertical index real(r8) :: kmq ! molecular diffusivity for constituent !----------------------------------------------------------------------- ! Determine superdiagonal (ca(k)) and subdiagonal (cc(k)) coeffs of the ! tridiagonal diffusion matrix. The diagonal elements (cb=1+ca+cc) are ! a combination of ca and cc; they are not required by the solver. do k = nbot-1, ntop, -1 do i = 1, ncol kmq = kq_scal(i,k+1) * mw_facm ca(i,k ) = (kv(i,k+1)+kmq) * tmpi(i,k+1) * rpdel(i,k ) cc(i,k+1) = (kv(i,k+1)+kmq) * tmpi(i,k+1) * rpdel(i,k+1) end do end do ! Calculate e(k). This term is ! required in solution of tridiagonal matrix defined by implicit diffusion eqn. do k = nbot, ntop+1, -1 do i=1,ncol dnom(i,k) = 1./ (1. + ca(i,k) + cc(i,k) - ca(i,k)*ze(i,k+1)) ze(i,k) = cc(i,k)*dnom(i,k) end do end do do i=1,ncol dnom(i,ntop) = 1./ (1. + ca(i,ntop) - ca(i,ntop)*ze(i,ntop+1)) end do return end subroutine vd_lu_qmolec !============================================================================== subroutine vd_lu_solve( & lchnk ,ncol , & q ,ca ,ze ,dnom , & ntop ,nbot ) !----------------------------------------------------------------------- ! Solve the implicit vertical diffusion equation with zero flux boundary conditions. ! Actual surface fluxes are explicit (rather than implicit) and are applied separately. ! Procedure for solution of the implicit equation follows Richtmyer and ! Morton (1967,pp 198-200). ! The equation solved is ! ! -ca(k)*q(k+1) + cb(k)*q(k) - cc(k)*q(k-1) = d(k), ! ! where d(k) is the input profile and q(k) is the output profile ! ! The solution has the form ! ! q(k) = ze(k)*q(k-1) + zf(k) ! ! ze(k) = cc(k) * dnom(k) ! ! zf(k) = [d(k) + ca(k)*zf(k+1)] * dnom(k) ! ! dnom(k) = 1/[cb(k) - ca(k)*ze(k+1)] = 1/[1 + ca(k) + cc(k) - ca(k)*ze(k+1)] ! Note that the same routine is used for temperature, momentum and tracers, ! and that input variables are replaced. !------------------------------Arguments-------------------------------- integer, intent(in) :: lchnk ! chunk identifier integer, intent(in) :: ncol ! number of atmospheric columns integer, intent(in) :: ntop ! top level to operate on integer, intent(in) :: nbot ! bottom level to operate on real(r8), intent(in) :: ca(pcols,pver) ! -upper diag coeff.of tri-diag matrix real(r8), intent(in) :: ze(pcols,pver) ! term in tri-diag solution real(r8), intent(in) :: dnom(pcols,pver) ! 1./(1. + ca(k) + cc(k) - ca(k)*ze(k+1)) real(r8), intent(inout) :: q(pcols,pver) ! constituent field !---------------------------Local workspace----------------------------- real(r8) :: zf(pcols,pver) ! term in tri-diag solution integer i,k ! longitude, vertical indices !----------------------------------------------------------------------- ! Calculate zf(k). Terms zf(k) and ze(k) are required in solution of ! tridiagonal matrix defined by implicit diffusion eqn. ! Note that only levels ntop through nbot need be solved for. do i = 1, ncol zf(i,nbot) = q(i,nbot)*dnom(i,nbot) end do do k = nbot-1, ntop, -1 do i=1,ncol zf(i,k) = (q(i,k) + ca(i,k)*zf(i,k+1))*dnom(i,k) end do end do ! Perform back substitution do i=1,ncol q(i,ntop) = zf(i,ntop) end do do k=ntop+1, nbot, +1 do i=1,ncol q(i,k) = zf(i,k) + ze(i,k)*q(i,k-1) end do end do return end subroutine vd_lu_solve !=============================================================================== subroutine trb_mtn_stress(lchnk, ncol , & u , v , t , pmid , exner , zm , sgh , & ksrf , taux , tauy ) !----------------------------------------------------------------------- ! Turbulent mountain stress parameterization ! ! Returns the surface stress associated with subgrid mountains ! For points where the orographic variance is small (including ocean), ! the returned stress is zero. !----------------------------------------------------------------------- !------------------------------Arguments-------------------------------- integer, intent(in) :: lchnk ! chunk identifier integer, intent(in) :: ncol ! number of atmospheric columns real(r8), intent(in) :: u(pcols,pver) ! midpoint zonal wind real(r8), intent(in) :: v(pcols,pver) ! midpoint meridional wind real(r8), intent(in) :: t(pcols,pver) ! midpoint temperatures real(r8), intent(in) :: pmid (pcols,pver) ! midpoint pressures real(r8), intent(in) :: exner(pcols,pver) ! exner function real(r8), intent(in) :: zm (pcols,pver) ! midpoint height real(r8), intent(in) :: sgh(pcols) ! standard deviation of orography real(r8), intent(out) :: ksrf(pcols) ! surface "drag" coefficient real(r8), intent(out) :: taux(pcols) ! turbulent mountain surface stress (zonal) real(r8), intent(out) :: tauy(pcols) ! turbulent mountain surface stress (meridional) !---------------------------Local storage------------------------------- integer :: i,k ! loop indexes integer :: kb,kt ! bottom and top of source region real(r8) :: horo ! orographic height real(r8) :: z0oro ! orographic z0 momentum real(r8) :: dv2 ! (delta v)**2 real(r8) :: ri ! richardson number real(r8) :: stabfri ! stability function of richardson number real(r8) :: rho ! density real(r8) :: cd ! drag coefficient real(r8) :: vmag ! velocity magnitude !--------------------------------------------------------------------------- do i = 1, ncol ! determine subgrid orgraphic height (trough to peak) horo = oroconst *sgh(i) ! no mountain stress if horo is too small if (horo < horomin) then ksrf(i) = 0. taux(i) = 0. tauy(i) = 0. else ! determine z0m for orography z0oro = min(z0fac * horo, z0max) ! calculate neutral drag coefficient k = pver cd = ( karman / log((zm(i,k) + z0oro)/z0oro) )**2 ! calculate the Richardson number over 1st 2 layers kt = pver-1 kb = pver dv2 = max((u(i,kt) - u(i,kb))**2 + (v(i,kt) - v(i,kb))**2, dv2min) ri = 2. * gravit * (t(i,kt)*exner(i,kt) - t(i,kb)*exner(i,kb)) * (zm(i,kt)-zm(i,kb)) & / ((t(i,kt) + t(i,kb)) * dv2) ! calculate the stability function and modify the neutral drag cofficient ! should probably follow Louis et al (1982) but for now just 1 for ri<0, 0 for ri>1, linear in 1-ri between 0 and 1 stabfri = max(0._r8,min(1._r8, 1. - ri)) cd = cd * stabfri ! compute density, velocity magnitude and stress using bottom level properties rho = pmid(i,k) / (rair * t(i,k)) vmag = sqrt(u(i,k)**2 + v(i,k)**2) ksrf(i) = rho * cd * vmag taux(i) = -ksrf(i) * u(i,k) tauy(i) = -ksrf(i) * v(i,k) end if end do return end subroutine trb_mtn_stress end module vertical_diffusion