! +-======-+ ! Copyright (c) 2003-2018 United States Government as represented by ! the Admistrator of the National Aeronautics and Space Administration. ! All Rights Reserved. ! ! THIS OPEN SOURCE AGREEMENT ("AGREEMENT") DEFINES THE RIGHTS OF USE, ! REPRODUCTION, DISTRIBUTION, MODIFICATION AND REDISTRIBUTION OF CERTAIN ! COMPUTER SOFTWARE ORIGINALLY RELEASED BY THE UNITED STATES GOVERNMENT AS ! REPRESENTED BY THE GOVERNMENT AGENCY LISTED BELOW ("GOVERNMENT AGENCY"). ! THE UNITED STATES GOVERNMENT, AS REPRESENTED BY GOVERNMENT AGENCY, IS AN ! INTENDED THIRD-PARTY BENEFICIARY OF ALL SUBSEQUENT DISTRIBUTIONS OR ! REDISTRIBUTIONS OF THE SUBJECT SOFTWARE. ANYONE WHO USES, REPRODUCES, ! DISTRIBUTES, MODIFIES OR REDISTRIBUTES THE SUBJECT SOFTWARE, AS DEFINED ! HEREIN, OR ANY PART THEREOF, IS, BY THAT ACTION, ACCEPTING IN FULL THE ! RESPONSIBILITIES AND OBLIGATIONS CONTAINED IN THIS AGREEMENT. ! ! Government Agency: National Aeronautics and Space Administration ! Government Agency Original Software Designation: GSC-15354-1 ! Government Agency Original Software Title: GEOS-5 GCM Modeling Software ! User Registration Requested. Please Visit http://opensource.gsfc.nasa.gov ! Government Agency Point of Contact for Original Software: ! Dale Hithon, SRA Assistant, (301) 286-2691 ! ! +-======-+ module shoc ! Implementation of the Simplified High Order Closure (SHOC) scheme ! of Bogenschutz and Krueger (2013), J. Adv. Model. Earth Syst, 5, 195-211, ! doi: 10.1002/jame.200118. (further referred to as BK13) ! in a single column form suitable for use in a GCM physics package. ! Alex Belochitski, heavily based on the code of Peter Bogenschutz. ! S Moorthi - optimization, cleanup, improve and customize for gsm ! ! N Arnold - implemented in GEOS ! use MAPL_ConstantsMod, only: ggr => MAPL_GRAV, & cp => MAPL_CP, & rgas => MAPL_RGAS, & rv => MAPL_RVAP, & lcond => MAPL_ALHL, & ! lfus => MAPL_ALHS, & pi => MAPL_PI, & MAPL_H2OMW, MAPL_AIRMW use MAPL_Mod, only: MAPL_UNDEF use MAPL_SatVaporMod, only: MAPL_EQsat implicit none private public run_shoc contains subroutine run_shoc( nx, ny, nzm, nz, dtn, & ! in prsl_inv, phii_inv, phil_inv, u_inv, v_inv, & ! in omega_inv, hflx, evap, & ! in tabs_inv, qwv_inv, qi_inv, qc_inv, qpi_inv, & ! inout qpl_inv, cld_sgs_inv, & ! inout tke_inv, prnum, tkh_inv, wthv_sec_inv, & ! inout tkesbdiss_inv, tkesbbuoy_inv, & ! out tkesbshear_inv,tkesbtrans_inv, & ! out smixt_inv,brunt_inv,shear_inv, & ! out LAM, TSCL, VON, CKVAL, CEFAC, CESFAC, & ! tuning param in THLTUN, QWTUN, QWTHLTUN, DO_TRANS ) real, parameter :: lfus = 1.0 real, parameter :: lsub = lcond+lfus, fac_cond = lcond/cp, fac_fus = lfus/cp, & fac_sub = lsub/cp, ggri = 1.0/ggr, kapa = rgas/cp, & gocp = ggr/cp, rog = rgas*ggri, sqrt2 = sqrt(2.0), & sqrtpii = 1.0/sqrt(pi+pi), epsterm = rgas/rv, twoby3 = 2.0/3.0, & onebeps = 1.0/epsterm, twoby15 = 2.0 / 15.0, & onebrvcp= 1.0/(rv*cp), skew_facw=1.2, skew_fact=1.0,& tkef1=0.5, tkef2=1.0-tkef1, tkhmax=500.0, epsv=MAPL_H2OMW/MAPL_AIRMW ! skew_facw=1.2, skew_fact=0.5 ! onebeps = 1.0/epsterm, twoby15 = 2.0 / 15.0, skew_facw=1.2 ! orig integer, intent(in) :: nx ! Number of points in the physics window in the x integer, intent(in) :: ny ! and y directions integer, intent(in) :: nzm ! Number of vertical layers integer, intent(in) :: nz ! Number of layer interfaces (= nzm + 1) integer, intent(in) :: DO_TRANS ! TKE transport term on/off real, intent(in) :: dtn ! Physics time step, s real, intent(in) :: hflx(nx) real, intent(in) :: evap(nx) real, intent(in) :: prsl_inv (nx,ny,nzm) ! mean layer presure real, intent(in) :: phii_inv (nx,ny,nz ) ! interface geopotential height real, intent(in) :: phil_inv (nx,ny,nzm) ! layer geopotential height real, intent(in) :: u_inv (nx,ny,nzm) ! u-wind, m/s real, intent(in) :: v_inv (nx,ny,nzm) ! v-wind, m/s real, intent(in) :: omega_inv(nx,ny,nzm) ! omega, Pa/s real, intent(inout) :: tabs_inv (nx,ny,nzm) ! temperature, K real, intent(inout) :: qwv_inv (nx,ny,nzm) ! water vapor mixing ratio, kg/kg real, intent(inout) :: qc_inv (nx,ny,nzm) ! cloud water mixing ratio, kg/kg real, intent(inout) :: qi_inv (nx,ny,nzm) ! cloud ice mixing ratio, kg/kg real, intent(inout) :: qpl_inv (nx,ny,nzm) ! rain mixing ratio, kg/kg real, intent(inout) :: qpi_inv (nx,ny,nzm) ! snow mixing ratio, kg/kg real, intent(inout) :: cld_sgs_inv(nx,ny,nzm) ! sgs cloud fraction real, intent(inout) :: tke_inv (nx,ny,nzm) ! turbulent kinetic energy. m**2/s**2 real, intent(inout) :: tkh_inv (nx,ny,nzm) ! eddy diffusivity real, intent(inout) :: prnum (nx,ny,nzm) ! turbulent Prandtl number real, intent(inout) :: wthv_sec_inv (nx,ny,nzm) ! Buoyancy flux, K*m/s real, intent( out) :: tkesbdiss_inv(nx,ny,nzm) ! dissipation real, intent( out) :: tkesbbuoy_inv(nx,ny,nzm) ! buoyancy production real, intent( out) :: tkesbshear_inv(nx,ny,nzm) ! shear production real, intent( out) :: tkesbtrans_inv(nx,ny,nzm) ! tke transport real, intent( out) :: smixt_inv(nx,ny,nzm) ! dissipation length scale real, intent( out) :: brunt_inv(nx,ny,nzm) ! Brunt vaisala frequency real, intent( out) :: shear_inv(nx,ny,nzm) ! Brunt vaisala frequency real, intent(in ) :: LAM, TSCL, VON, CKVAL, CEFAC, & ! tuning param CESFAC, THLTUN, QWTUN, QWTHLTUN ! SHOC tunable parameters real :: lambda real, parameter :: min_tke = 1e-6 ! Minumum TKE value, m**2/s**2 real, parameter :: max_tke = 100.0 ! Maximum TKE value, m**2/s**2 ! real, parameter :: max_tke = 5. ! Maximum TKE value, m**2/s**2 ! Maximum turbulent eddy length scale, m real, parameter :: max_eddy_length_scale = 2000. ! Maximum "return-to-isotropy" time scale, s real, parameter :: max_eddy_dissipation_time_scale = 2000. ! Constants for the TKE dissipation term based on Deardorff (1980) real, parameter :: pt19=0.19, pt51=0.51 real, parameter :: Cs = 0.15 real :: Ck, Ce, Ces ! real, parameter :: Ce = Ck**3/(0.7*Cs**4) ! real, parameter :: Ce = Ck**3/(0.7*Cs**4) * 2.2 ! real, parameter :: Ce = Ck**3/(0.7*Cs**4) * 3.0 , Ces = Ce ! real, parameter :: Ce = Ck**3/(0.7*Cs**4) * 2.5 , Ces = Ce * 3.0 / 2.5 ! real, parameter :: Ces = Ce/0.7*3.0 ! real, parameter :: Ce = Ck**3/(0.7*Cs**4), Ces = Ce*3.0/0.7 ! Commented Moor ! real, parameter :: vonk=0.35 ! Von Karman constant real :: vonk, tscale real, parameter :: w_tol_sqd = 4.0e-04 ! Min vlaue of second moment of w real, parameter :: w_thresh = 0.0, thresh = 0.0 ! These parameters are a tie-in with a microphysical scheme ! Double check their values for the Zhao-Carr scheme. real, parameter :: tbgmin = 233.16 ! Minimum temperature for cloud water., K ! real, parameter :: tbgmin = 258.16 ! Minimum temperature for cloud water., K (ZC) ! real, parameter :: tbgmin = 253.16 ! Minimum temperature for cloud water., K real, parameter :: tbgmax = 273.16 ! Maximum temperature for cloud ice, K real, parameter :: a_bg = 1.0/(tbgmax-tbgmin) ! ! Parameters to tune the second order moments- No tuning is performed currently real :: thl2tune, qw2tune, qwthl2tune real, parameter :: thl_tol = 1.e-2, rt_tol = 1.e-4, basetemp = 300.0 integer, parameter :: nitr=6 ! Local variables. Note that pressure is in millibars in the SHOC code. real zl (nx,ny,nzm) ! height of the pressure levels above surface, m real zi (nx,ny,nz) ! height of the interface levels, m real adzl (nx,ny,nzm) ! layer thickness i.e. zi(k+1)-zi(k) - defined at levels real adzi (nx,ny,nz) ! level thickness i.e. zl(k)-zl(k-1) - defined at interface real hl (nx,ny,nzm) ! liquid/ice water static energy , K real qv (nx,ny,nzm) ! water vapor, kg/kg real qcl (nx,ny,nzm) ! liquid water (condensate), kg/kg real qci (nx,ny,nzm) ! ice water (condensate), kg/kg real w (nx,ny,nzm) ! z-wind, m/s real bet (nx,ny,nzm) ! ggr/tv0 real gamaz (nx,ny,nzm) ! ggr/cp*z real phii (nx,ny,nz) real phil (nx,ny,nzm) real prsl (nx,ny,nzm) real u (nx,ny,nzm) real v (nx,ny,nzm) real omega (nx,ny,nzm) real tabs (nx,ny,nzm) real qwv (nx,ny,nzm) real qc (nx,ny,nzm) real qi (nx,ny,nzm) real qpl (nx,ny,nzm) real qpi (nx,ny,nzm) real cld_sgs (nx,ny,nzm) real tke (nx,ny,nzm) real tkh (nx,ny,nzm) real wthv_sec(nx,ny,nzm) real tkesbdiss(nx,ny,nzm) real tkesbbuoy(nx,ny,nzm) real tkesbshear(nx,ny,nzm) real tkesbtrans(nx,ny,nzm) ! Moments of the trivariate double Gaussian PDF for the SGS total water mixing ratio ! SGS liquid/ice static energy, and vertical velocity real qw_sec (nx,ny,nzm) ! Second moment total water mixing ratio, kg^2/kg^2 real thl_sec (nx,ny,nzm) ! Second moment liquid/ice static energy, K^2 real qwthl_sec(nx,ny,nzm) ! Covariance tot. wat. mix. ratio and static energy, K*kg/kg real wqw_sec (nx,ny,nzm) ! Turbulent flux of tot. wat. mix., kg/kg*m/s real wthl_sec (nx,ny,nzm) ! Turbulent flux of liquid/ice static energy, K*m/s real w_sec (nx,ny,nzm) ! Second moment of vertical velocity, m**2/s**2 real w3 (nx,ny,nzm) ! Third moment of vertical velocity, m**3/s**3 real wqp_sec (nx,ny,nzm) ! Turbulent flux of precipitation, kg/kg*m/s ! Eddy length formulation real smixt (nx,ny,nzm) ! Turbulent length scale, m real isotropy (nx,ny,nzm) ! "Return-to-isotropy" eddy dissipation time scale, s ! real isotropy_debug (nx,ny,nzm) ! Return to isotropy scale, s without artificial limits real brunt (nx,ny,nzm) ! Moist Brunt-Vaisalla frequency, s^-1 real conv_vel2(nx,ny,nzm) ! Convective velocity scale cubed, m^3/s^3 ! Output of SHOC real diag_frac, diag_qn, diag_qi, diag_ql ! real diag_frac(nx,ny,nzm) ! SGS cloud fraction ! real diag_qn (nx,ny,nzm) ! SGS cloud+ice condensate, kg/kg ! real diag_qi (nx,ny,nzm) ! SGS ice condensate, kg/kg ! real diag_ql (nx,ny,nzm) ! SGS liquid condensate, kg/kg ! Horizontally averaged variables ! real conv_vel(nzm) ! Convective velocity scale cubed, m^3/s^3 real wqlsb (nzm) ! liquid water flux, kg/kg/ m/s real wqisb (nzm) ! ice flux, kg/kg m/s real thlv (nzm) ! Grid-scale level-average virtual potential temperature ! Local variables !, tkesbdiss, tkesbbuoy_debug & ! tkebuoy_sgs, total_water, tscale1_debug, brunt2 real, dimension(nx,ny,nzm) :: total_water, brunt2, def2, thv real, dimension(nx,ny) :: denom, numer, l_inf, cldarr real lstarn, depth, omn, betdz, bbb, term, qsatt, dqsat, & thedz, conv_var, tkes, skew_w, skew_qw, aterm, w1_1, w1_2, w2_1, & w2_2, w3var, thl1_1, thl1_2, thl2_1, thl2_2, qw1_1, qw1_2, qw2_1, & qw2_2, ql1, ql2, w_ql1, w_ql2, & r_qwthl_1, r_wqw_1, r_wthl_1, testvar, s1, s2, std_s1, std_s2, C1, C2, & thl_first, qw_first, w_first, Tl1_1, Tl1_2, betatest, pval, pkap, & w2thl, w2qw,w2ql, w2ql_1, w2ql_2, & thec, thlsec, qwsec, qwthlsec, wqwsec, wthlsec, thestd,dum, & cqt1, cthl1, cqt2, cthl2, qn1, qn2, qi1, qi2, omn1, omn2, & basetemp2, beta1, beta2, qs1, qs2, & esval1_1, esval2_1, esval1_2, esval2_2, om1, om2, & lstarn1, lstarn2, sqrtw2, sqrtthl, sqrtqt, & sqrtstd1, sqrtstd2, tsign, tvar, sqrtw2t, wqls, wqis, & sqrtqw2_1, sqrtqw2_2, sqrtthl2_1, sqrtthl2_2, sm, prespot, & corrtest1, corrtest2, wrk, wrk1, wrk2, wrk3, onema, tkeavg, dtqw, dtqi integer i,j,k,km1,ku,kd,ka,kb,kinv,strt,fnsh,cnvl ! set parameter values lambda = LAM ! used in return-to-isotropy timescale Ck = CKVAL ! Coeff in the eddy diffusivity - TKE relationship, see Eq. 7 in BK13 Ce = CEFAC*Ck**3/Cs**4 ! diss ~= Ce * sqrt(tke) Ces = CESFAC*Ce vonk = VON ! Von Karman constant Moorthi - as in GFS tscale = TSCL ! time scale set based off of similarity results of BK13, s thl2tune = THLTUN qw2tune = QWTUN qwthl2tune = QWTHLTUN wthl_sec = 0. ! Map GEOS variables to those of SHOC do k=1,nz do j=1,ny do i=1,nx kinv = nz-k+1 zi(i,j,kinv) = phii_inv(i,j,k) enddo enddo enddo do k=1,nzm do j=1,ny do i=1,nx kinv = nzm-k+1 zl(i,j,kinv) = phil_inv(i,j,k) tkh(i,j,kinv) = tkh_inv(i,j,k) prsl(i,j,kinv) = prsl_inv(i,j,k) u(i,j,kinv) = u_inv(i,j,k) v(i,j,kinv) = v_inv(i,j,k) omega(i,j,kinv) = omega_inv(i,j,k) tabs(i,j,kinv) = tabs_inv(i,j,k) qwv(i,j,kinv) = qwv_inv(i,j,k) qc(i,j,kinv) = qc_inv(i,j,k) qi(i,j,kinv) = qi_inv(i,j,k) cld_sgs(i,j,kinv) = cld_sgs_inv(i,j,k) tke(i,j,kinv) = tke_inv(i,j,k) wthv_sec(i,j,kinv) = wthv_sec_inv(i,j,k) enddo enddo enddo ! print *,'phil_inv=',phil_inv(50,50,:) ! print *,'phii_inv(1:4)=',phii_inv(20,20,1:4) ! print *,'phii_inv(nz-3:nz)=',phii_inv(20,20,nz-3:nz) ! print *,'zi(1:4)=',zi(20,20,1:4) ! print *,'zi(nz-3:nz)=',zi(20,20,nz-3:nz) ! print *,'zl=',zl(50,50,:) ! print *,'zi=',zi(50,50,:) ! print *,'prsl=',prsl(50,50,:) ! print *,'tabs=',tabs(50,50,:) ! print *,'u=',u(50,50,:) ! print *,'old tke(1,1,:)=',tke(1,1,:) ! ! move water from vapor to condensate if the condensate is negative ! do k=1,nzm do j=1,ny do i=1,nx if (qc(i,j,k) < 0.0) then wrk = qwv(i,j,k) + qc(i,j,k) if (wrk >= 0.0) then qwv(i,j,k) = wrk tabs(i,j,k) = tabs(i,j,k) - fac_cond * qc(i,j,k) ! add heat used for condensation qc(i,j,k) = 0.0 else ! if total water 0 or negative, and negative... qc(i,j,k) = 0.0 tabs(i,j,k) = tabs(i,j,k) + fac_cond * qwv(i,j,k) qwv(i,j,k) = 0.0 endif endif if (qi(i,j,k) < 0.0) then wrk = qwv(i,j,k) + qi(i,j,k) if (wrk >= 0.0) then qwv(i,j,k) = wrk tabs(i,j,k) = tabs(i,j,k) - fac_sub * qi(i,j,k) qi(i,j,k) = 0.0 else qi(i,j,k) = 0.0 tabs(i,j,k) = tabs(i,j,k) + fac_sub * qwv(i,j,k) qwv(i,j,k) = 0.0 endif endif enddo enddo enddo do k=1,nzm do j=1,ny do i=1,nx ! zl(i,j,k) = phil(i,j,k) !* ggri wrk = 1.0 / prsl(i,j,k) qv(i,j,k) = max(qwv(i,j,k), 0.0) thv(i,j,k) = tabs(i,j,k) * (1.0+epsv*qv(i,j,k)) w(i,j,k) = - rog * omega(i,j,k) * thv(i,j,k) * wrk qcl(i,j,k) = max(qc(i,j,k), 0.0) qci(i,j,k) = max(qi(i,j,k), 0.0) ! qpl(i,j,k) = 0.0 ! comment or remove when using with prognostic rain/snow qpi(i,j,k) = 0.0 ! comment or remove when using with prognostic rain/snow wqp_sec(i,j,k) = 0.0 ! Turbulent flux of precipiation ! total_water(i,j,k) = qcl(i,j,k) + qci(i,j,k) + qv(i,j,k) prespot = (100000.0*wrk) ** kapa ! Exner function bet(i,j,k) = ggr/(tabs(i,j,k)*prespot) ! Moorthi thv(i,j,k) = thv(i,j,k)*prespot ! Moorthi ! ! Lapse rate * height = reference temperature gamaz(i,j,k) = gocp * zl(i,j,k) ! Liquid/ice water static energy - ! Note the the units are degrees K hl(i,j,k) = tabs(i,j,k) + gamaz(i,j,k) - fac_cond*(qcl(i,j,k)+qpl(i,j,k)) & - fac_fus *(qci(i,j,k)+qpi(i,j,k)) w3(i,j,k) = 0.0 enddo enddo enddo ! print *,'wthv_sec/thv=',wthv_sec(20,20,:)/thv(20,20,:) ! Define vertical grid increments for later use in the vertical differentiation ! print *,'nx=',nx ! print *,'ny=',ny ! print *,'zl=',zl(50,50,:) ! print *,'zi=',zi(50,50,:) do k=2,nzm km1 = k - 1 do j=1,ny do i=1,nx adzi(i,j,k) = (zl(i,j,k) - zl(i,j,km1)) ! print *,'zl(i,j,k)=',zl(i,j,k) ! print *,'zl(i,j,km1)=',zl(i,j,km1) adzl(i,j,km1) = (zi(i,j,k) - zi(i,j,km1)) ! from 1 to nzm-1 enddo enddo enddo do j=1,ny do i=1,nx adzi(i,j,1) = (zl(i,j,1)-zi(i,j,1)) ! unused in the code adzi(i,j,nz) = zi(i,j,nz)-zl(i,j,nzm) ! at the top - probably unused adzl(i,j,nzm) = adzi(i,j,nzm) ! wthl_sec(i,j,1) = hflx(i)/cp wqw_sec(i,j,1) = evap(i) enddo enddo ! print *,'adzi=',adzi(20,20,:) ! print *,'adzl=',adzl(20,20,:) call tke_shoc() ! Integrate prognostic TKE equation forward in time ! print *,'new tke(1,1,:)=',tke(1,1,:) ! diagnose second order moments of the subgrid PDF following ! Redelsperger J.L., and G. Sommeria, 1986, JAS, 43, 2619-2635 sans the use of stabilty ! weighting functions - Result is in global variables w_sec, thl_sec, qw_sec, and qwthl_sec ! call diag_moments(total_water,tke,tkh) ! Second moment of vertical velocity. ! Note that Eq 6 in BK13 gives a different expression that is dependent on ! vertical gradient of grid scale vertical velocity do k=1,nzm ku = k+1 kd = k-1 ka = ku kb = k if (k == 1) then kd = k kb = ka elseif (k == nzm) then ku = k ka = kb endif do j=1,ny do i=1,nx if (tke(i,j,k) > 0.0) then ! JDC Modify ! wrk = 0.5*(tkh(i,j,ka)+tkh(i,j,kb))*(w(i,j,ku) - w(i,j,kd)) & ! / (sqrt(tke(i,j,k)) * (zl(i,j,ku) - zl(i,j,kd))) wrk = 0.0 w_sec(i,j,k) = max(twoby3 * tke(i,j,k) - twoby15 * wrk, 0.0) else w_sec(i,j,k) = 0.0 endif enddo enddo enddo do k=2,nzm km1 = k-1 do j=1,ny do i=1,nx ! Use backward difference in the vertical, use averaged values of "return-to-isotropy" ! time scale and diffusion coefficient wrk1 = 1.0 / adzi(i,j,k) ! adzi(k) = (zl(k)-zl(km1)) wrk3 = tkh(i,j,k) * wrk1 sm = 0.5*(isotropy(i,j,k)+isotropy(i,j,km1))*wrk1*wrk3 ! Tau*Kh/dz^2 ! SGS vertical flux liquid/ice water static energy. Eq 1 in BK13 wrk1 = hl(i,j,k) - hl(i,j,km1) wthl_sec(i,j,k) = - wrk3 * wrk1 ! SGS vertical flux of total water. Eq 2 in BK13 wrk2 = total_water(i,j,k) - total_water(i,j,km1) wqw_sec(i,j,k) = - wrk3 * wrk2 ! Second moment of liquid/ice water static energy. Eq 4 in BK13 thl_sec(i,j,k) = thl2tune * sm * wrk1 * wrk1 ! Second moment of total water mixing ratio. Eq 3 in BK13 qw_sec(i,j,k) = qw2tune * sm * wrk2 * wrk2 ! Covariance of total water mixing ratio and liquid/ice water static energy. ! Eq 5 in BK13 qwthl_sec(i,j,k) = qwthl2tune * sm * wrk1 * wrk2 enddo ! i loop enddo ! j loop enddo ! k loop ! These would be at the surface - do we need them? do j=1,ny do i=1,nx ! wthl_sec(i,j,1) = wthl_sec(i,j,2) ! wqw_sec(i,j,1) = wqw_sec(i,j,2) thl_sec(i,j,1) = thl_sec(i,j,2) qw_sec(i,j,1) = qw_sec(i,j,2) qwthl_sec(i,j,1) = qwthl_sec(i,j,2) enddo enddo ! Diagnose the third moment of SGS vertical velocity call canuto() ! Recover parameters of the subgrid PDF using diagnosed moments ! and calculate SGS cloudiness, condensation and it's effects on temeperature ! and moisture variables call assumed_pdf() ! do k=1,nz ! do j=1,ny ! do i=1,nx ! kinv = nz-k+1 ! tkh(i,j,kinv) = tkh(i,j,k) ! enddo ! enddo ! enddo do k=1,nzm do j=1,ny do i=1,nx kinv = nzm-k+1 tabs_inv(i,j,kinv) = tabs(i,j,k) tkh_inv(i,j,kinv) = tkh(i,j,k) qwv_inv(i,j,kinv) = qwv(i,j,k) qc_inv(i,j,kinv) = qc(i,j,k) qi_inv(i,j,kinv) = qi(i,j,k) cld_sgs_inv(i,j,kinv) = cld_sgs(i,j,k) tke_inv(i,j,kinv) = tke(i,j,k) wthv_sec_inv(i,j,kinv) = wthv_sec(i,j,k) tkesbdiss_inv(i,j,kinv) = tkesbdiss(i,j,k) tkesbbuoy_inv(i,j,kinv) = tkesbbuoy(i,j,k) tkesbtrans_inv(i,j,kinv) = tkesbtrans(i,j,k) tkesbshear_inv(i,j,kinv) = tkesbshear(i,j,k) smixt_inv(i,j,kinv) = smixt(i,j,k) brunt_inv(i,j,kinv) = brunt2(i,j,k) shear_inv(i,j,kinv) = def2(i,j,k) enddo enddo enddo contains subroutine tke_shoc() ! This subroutine solves the TKE equation, ! Heavily based on SAM's tke_full.f90 by Marat Khairoutdinov real grd,betdz,Cek,Cee,lstarn, lstarp, bbb, omn, omp,qsatt,dqsat, smix, & buoy_sgs,ratio,a_prod_sh,a_prod_bu,a_diss,a_prod_bu_debug, buoy_sgs_debug, & tscale1, wrk, wrk1, wtke, wtk2, rdtn integer i,j,k,ku,kd,itr rdtn = 1.0 / dtn call tke_shear_prod(def2) ! Calculate shear production of TKE ! print *,'def2=',sqrt(def2(20,20,:)) ! Ensure values of TKE are reasonable do k=1,nzm do j=1,ny do i=1,nx tke(i,j,k) = max(min_tke,tke(i,j,k)) tkesbdiss(i,j,k) = 0. tkesbshear(i,j,k) = 0. tkesbbuoy(i,j,k) = 0. tkesbtrans(i,j,k) = 0. enddo enddo enddo call eddy_length() ! Find turbulent mixing length call check_eddy() ! Make sure it's reasonable ! temporary estimate of transport within convective layers ! - need to adjust a_prod_tr to conserve tke ! - need to make pressure weighted, rather than height ! - add option to extend 'convection' across single stable layer if (DO_TRANS /= 0) then cnvl = 0 do i=1,nx do j=1,ny do k=1,nzm if (wthv_sec(i,j,k).ge.0.01 .and. cnvl.eq.0) then cnvl = 1 strt = k end if if (cnvl.eq.1 .and. wthv_sec(i,j,k).lt.0.01) then fnsh = k-1 if (strt 0.0 .and. numer(i,j) > 0.0) then l_inf(i,j) = 0.1 * (numer(i,j)/denom(i,j)) else l_inf(i,j) = 100. endif enddo enddo !Calculate length scale outside of cloud, Eq. 10 in BK13 (Eq. 4.12 in Pete's dissertation) do k=1,nzm kb = k-1 kc = k+1 do j=1,ny do i=1,nx ! vars module variable bet (=ggr/tv0) ; grid module variable adzi if (k == 1) then kb = 1 kc = 2 thedz = adzi(i,j,kc) elseif (k == nzm) then kb = nzm-1 kc = nzm thedz = adzi(i,j,k) else thedz = (adzi(i,j,kc)+adzi(i,j,k)) ! = (z(k+1)-z(k-1)) endif betdz = bet(i,j,k) / thedz tkes = sqrt(tke(i,j,k)) ! Compute local Brunt-Vaisalla frequency wrk = qcl(i,j,k) + qci(i,j,k) if (wrk > 0.0) then ! If in the cloud ! Find the in-cloud Brunt-Vaisalla frequency omn = qcl(i,j,k) / (wrk+1.e-20) ! Ratio of liquid water to total water ! Latent heat of phase transformation based on relative water phase content ! fac_cond = lcond/cp, fac_fus = lfus/cp lstarn = fac_cond + (1.-omn)*fac_fus ! Saturation mixing ratio over water/ice wrt temp based on relative water phase content ! qsatt = omn * qsatw(tabs(i,j,k),prsl(i,j,k)) & ! + (1.-omn) * qsati(tabs(i,j,k),prsl(i,j,k)) qsatt = omn * MAPL_EQsat(tabs(i,j,k),prsl(i,j,k),dtqw) & + (1.-omn) * MAPL_EQsat(tabs(i,j,k),prsl(i,j,k),dtqi,OverIce=.TRUE.) ! Derivative of saturation mixing ratio over water/ice wrt temp. based on relative water phase content ! dqsat = omn * dtqsatw(tabs(i,j,k),prsl(i,j,k)) & ! + (1.-omn) * dtqsati(tabs(i,j,k),prsl(i,j,k)) dqsat = omn * dtqw + (1.-omn) * dtqi ! liquid/ice moist static energy static energy divided by cp? bbb = (1. + epsv*qsatt-wrk-qpl(i,j,k)-qpi(i,j,k) & + 1.61*tabs(i,j,k)*dqsat) / (1.+lstarn*dqsat) ! Calculate Brunt-Vaisalla frequency using centered differences in the vertical brunt(i,j,k) = betdz*(bbb*(hl(i,j,kc)-hl(i,j,kb)) & + (bbb*lstarn - (1.+lstarn*dqsat)*tabs(i,j,k)) & * (total_water(i,j,kc)-total_water(i,j,kb)) & + (bbb*fac_cond - (1.+fac_cond*dqsat)*tabs(i,j,k))*(qpl(i,j,kc)-qpl(i,j,kb)) & + (bbb*fac_sub - (1.+fac_sub*dqsat)*tabs(i,j,k))*(qpi(i,j,kc)-qpi(i,j,kb)) ) else ! outside of cloud ! Find outside-of-cloud Brunt-Vaisalla frequency ! Only unsaturated air, rain and snow contribute to virt. pot. temp. ! liquid/ice moist static energy divided by cp? bbb = 1. + epsv*qv(i,j,k) - qpl(i,j,k) - qpi(i,j,k) brunt(i,j,k) = betdz*( bbb*(hl(i,j,kc)-hl(i,j,kb)) & + epsv*tabs(i,j,k)*(total_water(i,j,kc)-total_water(i,j,kb)) & + (bbb*fac_cond-tabs(i,j,k))*(qpl(i,j,kc)-qpl(i,j,kb)) & + (bbb*fac_sub -tabs(i,j,k))*(qpi(i,j,kc)-qpi(i,j,kb)) ) end if ! Reduction of mixing length in the stable regions (where B.-V. freq. > 0) is required. ! Here we find regions of Brunt-Vaisalla freq. > 0 for later use. if (brunt(i,j,k) >= 1e-5) then brunt2(i,j,k) = brunt(i,j,k) else brunt2(i,j,k) = 1e-5 endif ! Calculate turbulent length scale in the boundary layer. ! See Eq. 10 in BK13 (Eq. 4.12 in Pete's dissertation) ! Keep the length scale adequately small near the surface following Blackadar (1984) ! Note that this is not documented in BK13 and was added later for SP-CAM runs ! if (k == 1) then ! term = 600.*tkes ! smixt(i,j,k) = term + (0.4*zl(i,j,k)-term)*exp(-zl(i,j,k)*0.01) ! else ! tscale is the eddy turnover time scale in the boundary layer and is ! an empirically derived constant if (tkes > 0.0 .and. l_inf(i,j) > 0.0) then wrk1 = 1.0 / (tscale*tkes*vonk*zl(i,j,k)) wrk2 = 1.0 / (tscale*tkes*l_inf(i,j)) ! wrk1 = 1.0 / (vonk*zl(i,j,k)) ! wrk2 = 1.0 / (l_inf(i,j)) wrk1 = 1.0 / (wrk1 + wrk2 + 0.01 * brunt2(i,j,k) / tke(i,j,k)) ! smixt(i,j,k) = min(max_eddy_length_scale, 2.8284*sqrt(wrk1)/0.3) smixt(i,j,k) = min(max_eddy_length_scale, sqrt(wrk1)/0.3) endif ! overwrite with simple Blackadar ! wrk1 = 1./200. ! wrk2 = 1./(vonk*zl(i,j,k)) ! wrk1 = 1./(wrk1+wrk2) ! smixt(i,j,k) = min(max_eddy_length_scale, wrk1) ! end if ! not k=1 end do end do end do ! Now find the in-cloud turbulence length scale ! See Eq. 13 in BK13 (Eq. 4.18 in Pete's disseration) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Remove after coupling to subgrid PDF. !wthv_sec = -300/ggr*brunt*tk !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! determine cubed convective velocity scale (conv_vel2) inside the cloud ! call conv_scale() ! inlining the relevant code do j=1,ny do i=1,nx conv_vel2(i,j,1) = 0. ! Convective velocity scale cubed enddo enddo ! Integrate velocity scale in the vertical do k=2,nzm do j=1,ny do i=1,nx conv_vel2(i,j,k) = conv_vel2(i,j,k-1) & + 2.5*adzi(i,j,k)*bet(i,j,k)*wthv_sec(i,j,k) enddo enddo enddo do j=1,ny do i=1,nx if (cldarr(i,j) == 1) then ! If there's a cloud in this column kl = 0 ku = 0 do k=2,nzm-3 ! Look for the cloud base in this column ! thresh (=0) is a variable local to eddy_length(). Should be a module constant. wrk = qcl(i,j,k) + qci(i,j,k) if (wrk > thresh .and. kl == 0) then kl = k endif ! Look for the cloud top in this column if (wrk > thresh .and. qcl(i,j,k+1)+qci(i,j,k+1) <= thresh) then ku = k ! conv_vel2 (Cubed convective velocity scale) is calculated in conv_scale() ! Use the value of conv_vel2 at the top of the cloud. conv_var = conv_vel2(i,j,k)**(1./3.) endif ! Compute the mixing length scale for the cloud layer that we just found if (kl > 0 .and. ku > 0 .and. ku-kl > 1) then if (conv_var > 0) then ! If convective vertical velocity scale > 0 depth = (zl(i,j,ku)-zl(i,j,kl)) + adzl(i,j,kl) do kk=kl,ku ! in-cloud turbulence length scale, Eq. 13 in BK13 (Eq. 4.18) ! wrk = conv_var/(depth*sqrt(tke(i,j,kk))) ! wrk = wrk * wrk + 0.01*brunt2(i,j,kk)/tke(i,j,kk) ! JDC Modify (based on Eq 13 in BK13) wrk = conv_var/(depth*depth*sqrt(tke(i,j,kk))) wrk = wrk + 0.01*brunt2(i,j,kk)/tke(i,j,kk) smixt(i,j,kk) = min(max_eddy_length_scale, (1.0/0.3)*sqrt(1.0/wrk)) enddo endif ! If convective vertical velocity scale > 0 kl = 0. ku = 0. endif ! if inside the cloud layer enddo ! k=2,nzm-3 endif ! if in the cloudy column enddo ! i=1,nx enddo ! j=1,ny end subroutine eddy_length subroutine conv_scale() ! This subroutine calculates the cubed convective velocity scale needed ! for the definition of the length scale in clouds ! See Eq. 16 in BK13 (Eq. 4.21 in Pete's dissertation) integer i, j, k !!!!!!!!! !! A bug in formulation of conv_vel ! Obtain it by averaging conv_vel2 in the horizontal !!!!!!!!!! ! conv_vel(1)=0. ! Horizontally averaged convective velocity scale cubed do j=1,ny do i=1,nx conv_vel2(i,j,1) = 0. ! Convective velocity scale cubed enddo enddo ! Integrate velocity scale in the vertical do k=2,nzm ! conv_vel(k)=conv_vel(k-1) do j=1,ny do i=1,nx !********************************************************************** !Do not include grid-scale contribution to convective velocity scale in GCM applications ! conv_vel(k)=conv_vel(k-1)+2.5*adzi(k)*bet(k)*(tvwle(k)+tvws(k)) ! conv_vel(k)=conv_vel(k)+2.5*adzi(i,j,k)*bet(i,j,k)*(tvws(k)) !Do not include grid-scale contribution to convective velocity scale in GCM applications ! conv_vel2(i,j,k)=conv_vel2(i,j,k-1)+2.5*adzi(k)*bet(k)*(tvwle(k)+wthv_sec(i,j,k)) !********************************************************************** conv_vel2(i,j,k) = conv_vel2(i,j,k-1) & + 2.5*adzi(i,j,k)*bet(i,j,k)*wthv_sec(i,j,k) enddo enddo enddo end subroutine conv_scale subroutine check_eddy() ! This subroutine checks eddy length values integer i, j, k, kb, ks, zend real wrk ! real zstart, zthresh, qthresh ! Temporary kludge for marine stratocumulus under very strong inversions at coarse resolution ! Placement until some explicity PBL top is put in ! Not used. ! zthresh = 100. ! qthresh = -6.0 do k=1,nzm if (k == nzm) then kb = k else kb = k+1 endif do j=1,ny do i=1,nx wrk = 0.1*adzl(i,j,k) ! Minimum 0.1 of local dz smixt(i,j,k) = max(wrk, min(max_eddy_length_scale,smixt(i,j,k))) ! If chracteristic grid dimension in the horizontal< 1000m, set lengthscale to ! be not larger that that. ! if (sqrt(dx*dy) .le. 1000.) smixt(i,j,k)=min(sqrt(dx*dy),smixt(i,j,k)) if (qcl(i,j,kb) == 0 .and. qcl(i,j,k) > 0 .and. brunt(i,j,k) > 1.e-4) then !If just above the cloud top and atmosphere is stable, set to 0.1 of local dz smixt(i,j,k) = wrk endif end do ! i end do ! j end do ! k end subroutine check_eddy subroutine canuto() ! Subroutine impements an analytic expression for the third moment of SGS vertical velocity ! based on Canuto et at, 2001, JAS, 58, 1169-1172 (further referred to as C01) ! This allows to avoid having a prognostic equation for the third moment. ! Result is returned in a global variable w3 defined at the interface levels. ! Local variables integer i, j, k, kb, kc real bet2, f0, f1, f2, f3, f4, f5, iso, isosqr, & omega0, omega1, omega2, X0, Y0, X1, Y1, AA0, AA1, buoy_sgs2, & thedz, thedz2, cond, wrk, wrk1, wrk2, wrk3, avew ! ! See Eq. 7 in C01 (B.7 in Pete's dissertation) real, parameter :: c=7.0, a0=0.52/(c*c*(c-2.)), a1=0.87/(c*c), & a2=0.5/c, a3=0.6/(c*(c-2.)), a4=2.4/(3.*c+5.), & a5=0.6/(c*(3.*c+5)) !Moorthi a5=0.6/(c*(3.+5.*c)) ! do k=1,nzm do k=2,nzm kb = k-1 kc = k+1 do j=1,ny do i=1,nx if(k == 1) then kb = 1 kc = 2 thedz = adzl(i,j,kc) thedz2 = thedz elseif(k == nzm) then kb = nzm-1 kc = nzm thedz = adzl(i,j,k) thedz2 = thedz else ! thedz = adzl(i,j,k) ! thedz2 = adzl(i,j,kc)+adzl(i,j,k) thedz = adzl(i,j,k) ! Moorthi jul08 thedz2 = adzl(i,j,k)+adzl(i,j,kb) ! Moorthi jul08 endif ! if (abs(thedz).lt.1e-20) print *,'thedz=',thedz ! if (abs(thedz2).lt.1e-20) print *,'thedz2=',thedz2 if (abs(thedz).le.1e-10) thedz = sign(1e-10,thedz) if (abs(thedz).eq.1e-10) print *,'thedz' if (abs(thedz2).le.1e-10) thedz2 = sign(1e-10,thedz2) if (abs(thedz2).eq.1e-10) print *,'thedz2' thedz = 1. / thedz thedz2 = 1. / thedz2 iso = 0.5*(isotropy(i,j,k)+isotropy(i,j,kb)) isosqr = iso*iso ! Two-level average of "return-to-isotropy" time scale squared buoy_sgs2 = isosqr*0.5*(brunt(i,j,k)+brunt(i,j,kb)) bet2 = 0.5*(bet(i,j,k)+bet(i,j,kb)) !Two-level average of BV frequency squared ! Compute functions f0-f5, see Eq, 8 in C01 (B.8 in Pete's dissertation) avew = 0.5*(w_sec(i,j,k)+w_sec(i,j,kb)) if (abs(avew).ge.1e10) avew = sign(1e10,avew) if (abs(avew).eq.1e10) print *,'avew' cond = 1.2*sqrt(max(1.0e-16,2.*avew*avew*avew)) wrk1 = bet2*iso wrk2 = thedz2*wrk1*wrk1*iso wrk3 = thl_sec(i,j,kc) - thl_sec(i,j,kb) f0 = wrk2 * wrk1 * wthl_sec(i,j,k) * wrk3 wrk = wthl_sec(i,j,kc) - wthl_sec(i,j,kb) f1 = wrk2 * (wrk*wthl_sec(i,j,k) + 0.5*avew*wrk3) wrk1 = bet2*isosqr ! JDC Modify: different from SHOC and Canuto et al 2001 ! f2 = thedz*wrk1*wthl_sec(i,j,k)*(w_sec(i,j,k)-w_sec(i,j,kb)) & ! + (thedz2+thedz2)*bet(i,j,k)*isosqr*wrk f2 = thedz*wrk1*wthl_sec(i,j,k)*(w_sec(i,j,k)-w_sec(i,j,kb)) & + (thedz2+thedz2)*bet(i,j,k)*isosqr*avew*wrk ! JDC modify: different from SHOC and Canuto et al 2001 ! f3 = thedz2*wrk1*wrk + thedz*bet2*isosqr*(wthl_sec(i,j,k)*(tke(i,j,k)-tke(i,j,kb))) f3 = thedz2*wrk1*avew*wrk + thedz*wrk1*(wthl_sec(i,j,k)*(tke(i,j,k)-tke(i,j,kb))) wrk1 = thedz*iso*avew f4 = wrk1*(w_sec(i,j,k)-w_sec(i,j,kb) + tke(i,j,k)-tke(i,j,kb)) f5 = wrk1*(w_sec(i,j,k)-w_sec(i,j,kb)) ! Compute the "omega" terms, see Eq. 6 in C01 (B.6 in Pete's dissertation) dum = 1.-a5*buoy_sgs2 if (abs(dum).le.1e-16) dum = sign(1e-16,dum) if (abs(dum).eq.1e-16) print *,'1.-a5*buoy_sgs2' omega0 = a4 / dum omega1 = omega0 / (c+c) omega2 = omega1*f3+(5./4.)*omega0*f4 ! Compute the X0, Y0, X1, Y1 terms, see Eq. 5 a-b in C01 (B.5 in Pete's dissertation) dum = 1.-(a1+a3)*buoy_sgs2 if (abs(dum).le.1e-16) dum = sign(1e-16,dum) if (abs(dum).eq.1e-16) print *,'1.-(a1+a3)*buoy_sgs2' wrk1 = 1.0 / dum dum = 1.-a3*buoy_sgs2 if (abs(dum).le.1e-16) dum = sign(1e-16,dum) if (abs(dum).eq.1e-16) print *,'1.-a3*buoy_sgs2' wrk2 = 1.0 / dum X0 = wrk1 * (a2*buoy_sgs2*(1.-a3*buoy_sgs2)) Y0 = wrk2 * (2.*a2*buoy_sgs2*X0) X1 = wrk1 * (a0*f0+a1*f1+a2*(1.-a3*buoy_sgs2)*f2) Y1 = wrk2 * (2.*a2*(buoy_sgs2*X1+(a0/a1)*f0+f1)) ! Compute the A0, A1 terms, see Eq. 5d in C01 (B.5 in Pete's dissertation) AA0 = omega0*X0 + omega1*Y0 AA1 = omega0*X1 + omega1*Y1 + omega2 ! Finally, we have the third moment of w, see Eq. 4c in C01 (B.4 in Pete's dissertation) ! cond is an estimate of third moment from second oment - If the third moment is larger ! than the estimate - limit w3. dum = c-1.2*X0+AA0 if (abs(dum).le.1e-16) dum = sign(1e-16,dum) if (abs(dum).eq.1e-16) print *,'c-1.2*X0+AA0=',dum w3(i,j,k) = max(-cond, min(cond, (AA1-1.2*X1-1.5*f5)/dum)) ! Implemetation of the C01 approach in this subroutine is nearly complete ! (the missing part are Eqs. 5c and 5e which are very simple) ! therefore it's easy to diagnose other third order moments obtained in C01 using this code. end do end do end do do j=1,ny do i=1,nx w3(i,j,1) = w3(i,j,2) enddo enddo end subroutine canuto subroutine assumed_pdf() ! Compute SGS buoyancy flux, SGS cloud fraction, and SGS condensation ! using assumed analytic double-gaussian PDF for SGS vertical velocity, ! moisture, and liquid/ice water static energy, based on the ! general approach of Larson et al 2002, JAS, 59, 3519-3539, ! and Golaz et al 2002, JAS, 59, 3540-3551 ! References in the comments in this code are given to ! the Appendix A of Pete Bogenschutz's dissertation. ! Local variables integer i,j,k,ku,kd real wrk, wrk1, wrk2, wrk3, wrk4, bastoeps ! bastoeps = basetemp / epsterm ! Initialize for statistics do k=1,nzm wqlsb(k) = 0. wqisb(k) = 0. enddo DO k=1,nzm kd = k ku = k + 1 if (k == nzm) ku = k DO j=1,ny DO i=1,nx ! Initialize cloud variables to zero diag_qn = 0.0 diag_frac = 0.0 diag_ql = 0.0 diag_qi = 0.0 pval = prsl(i,j,k) pkap = (pval/100000.0) ** kapa ! Read in liquid/ice static energy, total water mixing ratio, ! and vertical velocity to variables PDF needs thl_first = hl(i,j,k) qw_first = total_water(i,j,k) ! w_first = 0.5*(w(i,j,kd)+w(i,j,ku)) w_first = w(i,j,k) ! GET ALL INPUT VARIABLES ON THE SAME GRID ! Points to be computed with relation to thermo point ! Read in points that need to be averaged w3var = 0.5*(w3(i,j,kd)+w3(i,j,ku)) thlsec = max(0., 0.5*(thl_sec(i,j,kd)+thl_sec(i,j,ku)) ) qwsec = max(0., 0.5*(qw_sec(i,j,kd)+qw_sec(i,j,ku)) ) qwthlsec = 0.5 * (qwthl_sec(i,j,kd) + qwthl_sec(i,j,ku)) wqwsec = 0.5 * (wqw_sec(i,j,kd) + wqw_sec(i,j,ku)) wthlsec = 0.5 * (wthl_sec(i,j,kd) + wthl_sec(i,j,ku)) ! w3var = w3(i,j,k) ! thlsec = max(0.,thl_sec(i,j,k)) ! qwsec = max(0.,qw_sec(i,j,k)) ! qwthlsec = qwthl_sec(i,j,k) ! wqwsec = wqw_sec(i,j,k) ! wthlsec = wthl_sec(i,j,k) ! Compute square roots of some variables so we don't have to do it again if (w_sec(i,j,k) > 0.0) then sqrtw2 = sqrt(w_sec(i,j,k)) else sqrtw2 = 0.0 endif if (thlsec > 0.0) then sqrtthl = sqrt(thlsec) else sqrtthl = 0.0 endif if (qwsec > 0.0) then sqrtqt = sqrt(qwsec) else sqrtqt = 0.0 endif ! Find parameters of the double Gaussian PDF of vertical velocity ! Skewness of vertical velocity ! Skew_w = w3var / w_sec(i,j,k)**(3./2.) ! Skew_w = w3var / (sqrtw2*sqrtw2*sqrtw2) ! Moorthi IF (w_sec(i,j,k) <= w_tol_sqd) THEN ! If variance of w is too small then ! PDF is a sum of two delta functions Skew_w = 0. w1_1 = w_first w1_2 = w_first w2_1 = 0. w2_2 = 0. aterm = 0.5 onema = 0.5 ELSE Skew_w = w3var / (sqrtw2*sqrtw2*sqrtw2) ! Moorthi ! Proportionality coefficients between widths of each vertical velocity ! gaussian and the sqrt of the second moment of w w2_1 = 0.4 w2_2 = 0.4 ! Compute realtive weight of the first PDF "plume" ! See Eq A4 in Pete's dissertaion - Ensure 0.01 < a < 0.99 wrk = 1.0 - w2_1 aterm = max(0.01,min(0.5*(1.-Skew_w*sqrt(1./(4.*wrk*wrk*wrk+Skew_w*Skew_w))),0.99)) onema = 1.0 - aterm sqrtw2t = sqrt(wrk) ! Eq. A.5-A.6 wrk = sqrt(onema/aterm) w1_1 = sqrtw2t * wrk w1_2 = - sqrtw2t / wrk w2_1 = w2_1 * w_sec(i,j,k) w2_2 = w2_2 * w_sec(i,j,k) ENDIF ! Find parameters of the PDF of liquid/ice static energy IF (thlsec <= thl_tol*thl_tol .or. abs(w1_2-w1_1) <= w_thresh) THEN thl1_1 = thl_first thl1_2 = thl_first thl2_1 = 0. thl2_2 = 0. sqrtthl2_1 = 0. sqrtthl2_2 = 0. ELSE corrtest1 = max(-1.0,min(1.0,wthlsec/(sqrtw2*sqrtthl))) thl1_1 = -corrtest1 / w1_2 ! A.7 thl1_2 = -corrtest1 / w1_1 ! A.8 wrk1 = thl1_1 * thl1_1 wrk2 = thl1_2 * thl1_2 wrk3 = 1.0 - aterm*wrk1 - onema*wrk2 wrk4 = -skew_fact*Skew_w - aterm*wrk1*thl1_1 - onema*wrk2*thl1_2 ! testing - Moorthi ! wrk4 = - aterm*wrk1*thl1_1 - onema*wrk2*thl1_2 wrk = 3. * (thl1_2-thl1_1) if (wrk /= 0.0) then thl2_1 = thlsec * min(100.,max(0.,( 3.*thl1_2*wrk3-wrk4)/(aterm*wrk))) ! A.10 thl2_2 = thlsec * min(100.,max(0.,(-3.*thl1_1*wrk3+wrk4)/(onema*wrk))) ! A.11 else thl2_1 = 0.0 thl2_2 = 0.0 endif ! thl1_1 = thl1_1*sqrtthl + thl_first thl1_2 = thl1_2*sqrtthl + thl_first sqrtthl2_1 = sqrt(thl2_1) sqrtthl2_2 = sqrt(thl2_2) ENDIF ! FIND PARAMETERS FOR TOTAL WATER MIXING RATIO IF (qwsec <= rt_tol*rt_tol .or. abs(w1_2-w1_1) <= w_thresh) THEN qw1_1 = qw_first qw1_2 = qw_first qw2_1 = 0. qw2_2 = 0. sqrtqw2_1 = 0. sqrtqw2_2 = 0. ELSE corrtest2 = max(-1.0,min(1.0,wqwsec/(sqrtw2*sqrtqt))) qw1_1 = - corrtest2 / w1_2 ! A.7 qw1_2 = - corrtest2 / w1_1 ! A.8 tsign = abs(qw1_2-qw1_1) Skew_qw = skew_facw*Skew_w ! IF (tsign > 0.4) THEN ! Skew_qw = skew_facw*Skew_w ! ELSE IF (tsign <= 0.2) THEN ! Skew_qw = 0. ! ELSE ! Skew_qw = (skew_facw/0.2) * Skew_w * (tsign-0.2) ! ENDIF wrk1 = qw1_1 * qw1_1 wrk2 = qw1_2 * qw1_2 wrk3 = 1. - aterm*wrk1 - onema*wrk2 wrk4 = Skew_qw - aterm*wrk1*qw1_1 - onema*wrk2*qw1_2 wrk = 3. * (qw1_2-qw1_1) if (wrk /= 0.0) then qw2_1 = qwsec * min(100.,max(0.,( 3.*qw1_2*wrk3-wrk4)/(aterm*wrk))) ! A.10 qw2_2 = qwsec * min(100.,max(0.,(-3.*qw1_1*wrk3+wrk4)/(onema*wrk))) ! A.11 else qw2_1 = 0.0 qw2_2 = 0.0 endif ! qw1_1 = qw1_1*sqrtqt + qw_first qw1_2 = qw1_2*sqrtqt + qw_first sqrtqw2_1 = sqrt(qw2_1) sqrtqw2_2 = sqrt(qw2_2) ENDIF ! CONVERT FROM TILDA VARIABLES TO "REAL" VARIABLES w1_1 = w1_1*sqrtw2 + w_first w1_2 = w1_2*sqrtw2 + w_first ! FIND WITHIN-PLUME CORRELATIONS testvar = aterm*sqrtqw2_1*sqrtthl2_1 + onema*sqrtqw2_2*sqrtthl2_2 IF (testvar == 0) THEN r_qwthl_1 = 0. ELSE r_qwthl_1 = max(-1.0,min(1.0,(qwthlsec-aterm*(qw1_1-qw_first)*(thl1_1-thl_first)-onema*(qw1_2-qw_first)*(thl1_2-thl_first))/testvar)) ! A.12 ENDIF ! BEGIN TO COMPUTE CLOUD PROPERTY STATISTICS Tl1_1 = thl1_1 - gamaz(i,j,k) + fac_cond*qpl(i,j,k) + fac_sub*qpi(i,j,k) Tl1_2 = thl1_2 - gamaz(i,j,k) + fac_cond*qpl(i,j,k) + fac_sub*qpi(i,j,k) ! Now compute qs esval1_1 = 0. esval1_2 = 0. esval2_1 = 0. esval2_2 = 0. om1 = 1. om2 = 1. ! Partition based on temperature for the first plume IF (Tl1_1 >= tbgmax) THEN ! esval1_1 = fpvs(Tl1_1) !! esval1_1 = fpvsl(Tl1_1) esval1_1 = MAPL_EQsat(Tl1_1) ! esval1_1 = esatw(Tl1_1) lstarn1 = lcond ELSE IF (Tl1_1 < tbgmin) THEN ! esval1_1 = fpvs(Tl1_1) !! esval1_1 = fpvsi(Tl1_1) esval1_1 = MAPL_EQsat(Tl1_1,OverIce=.TRUE.) ! esval1_1 = esati(Tl1_1) lstarn1 = lsub ELSE ! esval1_1 = fpvs(Tl1_1) ! esval2_1 = fpvs(Tl1_1) !! esval1_1 = fpvsl(Tl1_1) !! esval2_1 = fpvsi(Tl1_1) esval1_1 = MAPL_EQsat(Tl1_1) esval2_1 = MAPL_EQsat(Tl1_1,OverIce=.TRUE.) ! esval1_1 = esatw(Tl1_1) ! esval2_1 = esati(Tl1_1) om1 = max(0.,min(1.,a_bg*(Tl1_1-tbgmin))) lstarn1 = lcond + (1.-om1)*lfus ENDIF qs1 = om1 * (0.622*esval1_1/max(esval1_1,pval-0.378*esval1_1)) & + (1.-om1) * (0.622*esval2_1/max(esval2_1,pval-0.378*esval2_1)) ! qs1 = om1 * (0.622*esval1_1/max(esval1_1,pval-esval1_1)) & ! + (1.-om1) * (0.622*esval2_1/max(esval2_1,pval-esval2_1)) ! beta1 = (rgas/rv)*(lstarn1/(rgas*Tl1_1))*(lstarn1/(cp*Tl1_1)) ! JDC Modify: Moorthi used beta2, likely a bug. switch to beta1 ! beta2 = (lstarn1*lstarn1*onebrvcp) / (Tl1_1*Tl1_1) ! A.18 beta1 = (lstarn1*lstarn1*onebrvcp) / (Tl1_1*Tl1_1) ! Are the two plumes equal? If so then set qs and beta ! in each column to each other to save computation IF (Tl1_1 == Tl1_2) THEN qs2 = qs1 beta2 = beta1 ELSE IF (Tl1_2 < tbgmin) THEN ! esval1_2 = fpvs(Tl1_2) !! esval1_2 = fpvsi(Tl1_2) esval1_2 = MAPL_EQsat(Tl1_2,OverIce=.TRUE.) ! esval1_2 = esati(Tl1_2) lstarn2 = lsub ELSE IF (Tl1_2 >= tbgmax) THEN ! esval1_2 = fpvs(Tl1_2) !! esval1_2 = fpvsl(Tl1_2) esval1_2 = MAPL_EQsat(Tl1_2) ! esval1_2 = esatw(Tl1_2) lstarn2 = lcond ELSE ! esval1_2 = fpvs(Tl1_2) ! esval2_2 = fpvs(Tl1_2) !! esval1_2 = fpvsl(Tl1_2) !! esval2_2 = fpvsi(Tl1_2) esval1_2 = MAPL_EQsat(Tl1_2) esval2_2 = MAPL_EQsat(Tl1_2,OverIce=.TRUE.) ! esval1_2 = esatw(Tl1_2) ! esval2_2 = esati(Tl1_2) om2 = max(0.,min(1.,a_bg*(Tl1_2-tbgmin))) lstarn2 = lcond + (1.-om2)*lfus ENDIF qs2 = om2 * (0.622*esval1_2/max(esval1_2,pval-0.378*esval1_2)) & + (1.-om2) * (0.622*esval2_2/max(esval2_2,pval-0.378*esval2_2)) ! qs2 = om2 * (0.622*esval1_2/max(esval1_2,pval-esval1_2)) & ! + (1.-om2) * (0.622*esval2_2/max(esval2_2,pval-esval2_2)) ! beta2 = (rgas/rv)*(lstarn2/(rgas*Tl1_2))*(lstarn2/(cp*Tl1_2)) ! A.18 beta2 = (lstarn2*lstarn2*onebrvcp) / (Tl1_2*Tl1_2) ! A.18 ENDIF ! qs1 = qs1 * rhc(i,j,k) ! qs2 = qs2 * rhc(i,j,k) ! Now compute cloud stuff - compute s term cqt1 = 1.0 / (1.0+beta1*qs1) ! A.19 wrk = (1.0+beta1*qw1_1) * cqt1 s1 = qw1_1 - qs1* wrk ! A.17 cthl1 = cqt1*wrk*(cp/lcond)*beta1*qs1*pkap ! A.20 wrk1 = cthl1 * cthl1 wrk2 = cqt1 * cqt1 std_s1 = sqrt(max(0.,wrk1*thl2_1+wrk2*qw2_1-2.*cthl1*sqrtthl2_1*cqt1*sqrtqw2_1*r_qwthl_1)) qn1 = 0. C1 = 0. IF (std_s1 /= 0) THEN wrk = s1 / (std_s1*sqrt2) C1 = 0.5*(1.+erf(wrk)) ! A.15 IF (C1 /= 0) qn1 = s1*C1 + (std_s1*sqrtpii)*exp(-wrk*wrk) ! A.16 ELSEIF (s1 > 0) THEN C1 = 1.0 qn1 = s1 ENDIF ! now compute non-precipitating cloud condensate ! If two plumes exactly equal, then just set many of these ! variables to themselves to save on computation. IF (qw1_1 == qw1_2 .and. thl2_1 == thl2_2 .and. qs1 == qs2) THEN s2 = s1 cthl2 = cthl1 cqt2 = cqt1 std_s2 = std_s1 C2 = C1 qn2 = qn1 ELSE cqt2 = 1.0 / (1.0+beta2*qs2) wrk = (1.0+beta2*qw1_2) * cqt2 s2 = qw1_2 - qs2*wrk cthl2 = wrk*cqt2*(cp/lcond)*beta2*qs2*pkap wrk1 = cthl2 * cthl2 wrk2 = cqt2 * cqt2 std_s2 = sqrt(max(0.,wrk1*thl2_2+wrk2*qw2_2-2.*cthl2*sqrtthl2_2*cqt2*sqrtqw2_2*r_qwthl_1)) qn2 = 0. C2 = 0. IF (std_s2 /= 0) THEN wrk = s2 / (std_s2*sqrt2) C2 = 0.5*(1.+erf(wrk)) IF (C2 /= 0) qn2 = s2*C2 + (std_s2*sqrtpii)*exp(-wrk*wrk) ELSEIF (s2 > 0) THEN C2 = 1.0 qn2 = s2 ENDIF ENDIF ! finally, compute the SGS cloud fraction diag_frac = aterm*C1 + onema*C2 om1 = max(0.,min(1.,(Tl1_1-tbgmin)*a_bg)) om2 = max(0.,min(1.,(Tl1_2-tbgmin)*a_bg)) qn1 = min(qn1,qw1_1) qn2 = min(qn2,qw1_2) ql1 = qn1*om1 ql2 = qn2*om2 qi1 = qn1 - ql1 qi2 = qn2 - ql2 diag_qn = min(max(0.0, aterm*qn1 + onema*qn2), total_water(i,j,k)) diag_ql = min(max(0.0, aterm*ql1 + onema*ql2), diag_qn) diag_qi = diag_qn - diag_ql ! Update temperature variable based on diagnosed cloud properties om1 = max(0.,min(1.,(tabs(i,j,k)-tbgmin)*a_bg)) lstarn1 = lcond + (1.-om1)*lfus tabs(i,j,k) = hl(i,j,k) - gamaz(i,j,k) + fac_cond*(diag_ql+qpl(i,j,k)) & + fac_sub *(diag_qi+qpi(i,j,k)) & + tkesbdiss(i,j,k) * (dtn/cp) ! tke dissipative heating ! Update moisture fields qc(i,j,k) = diag_ql qi(i,j,k) = diag_qi qwv(i,j,k) = total_water(i,j,k) - diag_qn cld_sgs(i,j,k) = diag_frac ! Compute the liquid water flux wqls = aterm * ((w1_1-w_first)*ql1) + onema * ((w1_2-w_first)*ql2) wqis = aterm * ((w1_1-w_first)*qi1) + onema * ((w1_2-w_first)*qi2) ! Compute statistics for the fluxes so we don't have to save these variables wqlsb(k) = wqlsb(k) + wqls wqisb(k) = wqisb(k) + wqis ! diagnostic buoyancy flux. Includes effects from liquid water, ice ! condensate, liquid & ice precipitation ! wrk = epsv * basetemp wrk = epsv * thv(i,j,k) bastoeps = onebeps * thv(i,j,k) wthv_sec(i,j,k) = wthlsec + wrk*wqwsec & + (fac_cond-bastoeps)*wqls & + (fac_sub-bastoeps) *wqis & + ((lstarn1/cp)-thv(i,j,k))*0.5*(wqp_sec(i,j,kd)+wqp_sec(i,j,ku)) ! wthv_sec(i,j,k) = wthlsec + wrk*wqwsec & ! + (fac_cond-bastoeps)*wqls & ! + (fac_sub-bastoeps)*wqis & ! + ((lstarn1/cp)-basetemp)*0.5*(wqp_sec(i,j,kd)+wqp_sec(i,j,ku)) ENDDO ENDDO ENDDO end subroutine assumed_pdf ! Saturation vapor pressure and mixing ratio subroutines ! Based on Flatau et al (1992), J. App. Met., 31, 1507-1513 ! Code by Marat Khairoutdinov real function esatw(t) !!! returns e in mb real t ! temperature (K) real a0,a1,a2,a3,a4,a5,a6,a7,a8 data a0,a1,a2,a3,a4,a5,a6,a7,a8 / & 6.11239921, 0.443987641, 0.142986287e-1, & 0.264847430e-3, 0.302950461e-5, 0.206739458e-7, & 0.640689451e-10, -0.952447341e-13,-0.976195544e-15/ real dt dt = max(-80.,t-273.16) esatw = a0 + dt*(a1+dt*(a2+dt*(a3+dt*(a4+dt*(a5+dt*(a6+dt*(a7+a8*dt))))))) end function esatw real function qsatw(t,p) ! implicit none real t ! temperature (K) real p ! pressure (Pa) real esat ! esat = fpvs(t) !! esat = fpvsl(t) esat = MAPL_EQsat(t) qsatw = 0.622 * esat/max(esat,p-0.378*esat) ! esat = esatw(t) ! qsatw = 0.622 * esat/max(esat,p-esat) end function qsatw real function esati(t) !!! returns e in mb real t ! temperature (K) real a0,a1,a2,a3,a4,a5,a6,a7,a8 data a0,a1,a2,a3,a4,a5,a6,a7,a8 / & 6.11147274, 0.503160820, 0.188439774e-1, & 0.420895665e-3, 0.615021634e-5, 0.602588177e-7, & 0.385852041e-9, 0.146898966e-11, 0.252751365e-14/ real dt ! real esatw if(t > 273.15) then esati = esatw(t) else if(t.gt.185.) then dt = t-273.16 esati = a0 + dt*(a1+dt*(a2+dt*(a3+dt*(a4+dt*(a5+dt*(a6+dt*(a7+a8*dt))))))) else ! use some additional interpolation below 184K dt = max(-100.,t-273.16) esati = 0.00763685 + dt*(0.000151069+dt*7.48215e-07) end if end function esati real function qsati(t,p) real t ! temperature (K) real p ! pressure (Pa) real esat !,esati ! esat = fpvs(t) !! esat = fpvsi(t) esat = MAPL_EQsat(t,OverIce=.TRUE.) qsati = 0.622 * esat/max(esat,p-0.378*esat) ! esat = esati(t) ! qsati = 0.622 * esat/max(esat,p-esat) end function qsati real function dtesatw(t) real t ! temperature (K) real a0,a1,a2,a3,a4,a5,a6,a7,a8 data a0,a1,a2,a3,a4,a5,a6,a7,a8 / & 0.443956472, 0.285976452e-1, 0.794747212e-3, & 0.121167162e-4, 0.103167413e-6, 0.385208005e-9, & -0.604119582e-12, -0.792933209e-14, -0.599634321e-17/ real dt dt = max(-80.,t-273.16) dtesatw = a0 + dt* (a1+dt*(a2+dt*(a3+dt*(a4+dt*(a5+dt*(a6+dt*(a7+a8*dt))))))) end function dtesatw real function dtqsatw(t,p) real t ! temperature (K) real p ! pressure (Pa) ! real dtesatw dtqsatw = 100.0*0.622*dtesatw(t)/p end function dtqsatw real function dtesati(t) real t ! temperature (K) real a0,a1,a2,a3,a4,a5,a6,a7,a8 data a0,a1,a2,a3,a4,a5,a6,a7,a8 / & 0.503223089, 0.377174432e-1, 0.126710138e-2, & 0.249065913e-4, 0.312668753e-6, 0.255653718e-8, & 0.132073448e-10, 0.390204672e-13, 0.497275778e-16/ real dt ! real dtesatw if(t > 273.15) then dtesati = dtesatw(t) else if(t > 185.) then dt = t-273.16 dtesati = a0 + dt*(a1+dt*(a2+dt*(a3+dt*(a4+dt*(a5+dt*(a6+dt*(a7+a8*dt))))))) else ! use additional interpolation below 185K dt = max(-100.,t-273.16) dtesati = 0.0013186 + dt*(2.60269e-05+dt*1.28676e-07) end if end function dtesati real function dtqsati(t,p) real t ! temperature (K) real p ! pressure (Pa) ! real dtesati dtqsati = 100.0*0.622*dtesati(t)/p end function dtqsati end subroutine run_shoc end module shoc