diff --git a/physics/samfdeepcnv.f b/physics/samfdeepcnv.f index a7ea0009b..2f34041c2 100644 --- a/physics/samfdeepcnv.f +++ b/physics/samfdeepcnv.f @@ -151,7 +151,7 @@ subroutine samfdeepcnv_run (im,km,itc,ntc,cliq,cp,cvap, & & dt2, dtmax, dtmin, & dxcrtas, dxcrtuf, & dv1h, dv2h, dv3h, - & dv1q, dv2q, dv3q, + & dv2q, & dz, dz1, e1, edtmax, & edtmaxl, edtmaxs, el2orc, elocp, & es, etah, @@ -243,7 +243,11 @@ subroutine samfdeepcnv_run (im,km,itc,ntc,cliq,cp,cvap, & & uo(im,km), vo(im,km), qeso(im,km), & ctr(im,km,ntr), ctro(im,km,ntr) ! for aerosol transport - real(kind=kind_phys) qaero(im,km,ntc) +! real(kind=kind_phys) qaero(im,km,ntc) +c variables for tracer wet deposition, + real(kind=kind_phys), dimension(im,km,ntc) :: chem_c, chem_pw, + & wet_dep +! ! for updraft velocity calculation real(kind=kind_phys) wu2(im,km), buo(im,km), drag(im,km) real(kind=kind_phys) wc(im), scaldfunc(im), sigmagfm(im) @@ -266,6 +270,14 @@ subroutine samfdeepcnv_run (im,km,itc,ntc,cliq,cp,cvap, & & pwo(im,km), pwdo(im,km), c0t(im,km), & tx1(im), sumx(im), cnvwt(im,km) &, rhbar(im) +! +! variables for Total Variation Diminishing (TVD) flux-limiter scheme +! on environmental subsidence and uplifting +! + real(kind=kind_phys) q_diff(im,0:km-1), e_diff(im,0:km-1,ntr), + & flxtvd(im,0:km-1) + real(kind=kind_phys) rrkp, phkp + real(kind=kind_phys) tsumn(im), tsump(im), rtnp(im) ! logical do_aerosols, totflg, cnvflg(im), asqecflg(im), flg(im) ! @@ -318,6 +330,10 @@ subroutine samfdeepcnv_run (im,km,itc,ntc,cliq,cp,cvap, & c c initialize arrays c + chem_c = 0. + chem_pw = 0. + wet_dep = 0. +! do i=1,im cnvflg(i) = .true. sfcpbl(i) = sfclfac * hpbl(i) @@ -1164,7 +1180,12 @@ subroutine samfdeepcnv_run (im,km,itc,ntc,cliq,cp,cvap, & enddo enddo if (.not.hwrf_samfdeep) then - do n = 1, ntr + if (do_aerosols) then + kk = itc -3 + else + kk = ntr + endif + do n = 1, kk do k = 2, km1 do i = 1, im if (cnvflg(i)) then @@ -1179,6 +1200,28 @@ subroutine samfdeepcnv_run (im,km,itc,ntc,cliq,cp,cvap, & enddo enddo enddo + if (do_aerosols) then + do n = 1, ntc + kk = n + itc -3 + do k = 2, km1 + do i = 1, im + if (cnvflg(i)) then + if(k > kb(i) .and. k < kmax(i)) then + dz = zi(i,k) - zi(i,k-1) + tem = 0.25 * (xlamue(i,k)+xlamue(i,k-1)) * dz + factor = 1. + tem + ecko(i,k,kk) = ((1. - tem) * ecko(i,k-1,kk) + tem * + & (ctro(i,k,kk) + ctro(i,k-1,kk))) / factor + chem_c(i,k,n) = fscav(n) * ecko(i,k,kk) + tem = chem_c(i,k,n) / (1. + c0t(i,k) * dz) + chem_pw(i,k,n) = c0t(i,k) * dz * tem * eta(i,k-1) + ecko(i,k,kk) = tem + ecko(i,k,kk) - chem_c(i,k,n) + endif + endif + enddo + enddo + enddo + endif endif c c taking account into convection inhibition due to existence of @@ -1749,8 +1792,8 @@ subroutine samfdeepcnv_run (im,km,itc,ntc,cliq,cp,cvap, & do i = 1, im if (cnvflg(i)) then if(k > kb(i) .and. k <= ktcon(i)) then - shear= sqrt((uo(i,k)-uo(i,k-1)) ** 2 - & + (vo(i,k)-vo(i,k-1)) ** 2) + shear = sqrt((uo(i,k)-uo(i,k-1)) ** 2 + & + (vo(i,k)-vo(i,k-1)) ** 2) vshear(i) = vshear(i) + shear endif endif @@ -1949,18 +1992,18 @@ subroutine samfdeepcnv_run (im,km,itc,ntc,cliq,cp,cvap, & do i = 1, im if (cnvflg(i) .and. k < jmin(i)) then gamma = el2orc * qeso(i,k) / to(i,k)**2 - dhh=hcdo(i,k) - dt=to(i,k) - dg=gamma - dh=heso(i,k) - dz=-1.*(zo(i,k+1)-zo(i,k)) + dhh = hcdo(i,k) + dt = to(i,k) + dg = gamma + dh = heso(i,k) + dz = zo(i,k) - zo(i,k+1) ! aa1(i)=aa1(i)+edto(i)*dz*etad(i,k) - aa1(i)=aa1(i)+edto(i)*dz - & *(grav/(cp*dt))*((dhh-dh)/(1.+dg)) - & *(1.+fv*cp*dg*dt/hvap) - val=0. + aa1(i) = aa1(i)+edto(i)*dz + & *(grav/(cp*dt))*((dhh-dh)/(1.+dg)) + & *(1.+fv*cp*dg*dt/hvap) + val = 0. ! aa1(i)=aa1(i)+edto(i)*dz*etad(i,k) - aa1(i)=aa1(i)+edto(i)*dz + aa1(i) = aa1(i)+edto(i)*dz & *grav*fv*max(val,(qeso(i,k)-qo(i,k))) endif enddo @@ -2007,14 +2050,11 @@ subroutine samfdeepcnv_run (im,km,itc,ntc,cliq,cp,cvap, & do i = 1, im if(cnvflg(i)) then dp = 1000. * del(i,1) - dellah(i,1) = edto(i) * etad(i,1) * (hcdo(i,1) - & - heo(i,1)) * grav / dp - dellaq(i,1) = edto(i) * etad(i,1) * (qrcdo(i,1) - & - qo(i,1)) * grav / dp - dellau(i,1) = edto(i) * etad(i,1) * (ucdo(i,1) - & - uo(i,1)) * grav / dp - dellav(i,1) = edto(i) * etad(i,1) * (vcdo(i,1) - & - vo(i,1)) * grav / dp + tem = edto(i) * etad(i,1) * grav / dp + dellah(i,1) = tem * (hcdo(i,1) - heo(i,1)) + dellaq(i,1) = tem * qrcdo(i,1) + dellau(i,1) = tem * (ucdo(i,1) - uo(i,1)) + dellav(i,1) = tem * (vcdo(i,1) - vo(i,1)) endif enddo if (.not.hwrf_samfdeep) then @@ -2022,8 +2062,8 @@ subroutine samfdeepcnv_run (im,km,itc,ntc,cliq,cp,cvap, & do i = 1, im if(cnvflg(i)) then dp = 1000. * del(i,1) - dellae(i,1,n) = edto(i) * etad(i,1) * (ecdo(i,1,n) - & - ctro(i,1,n)) * grav / dp + dellae(i,1,n) = edto(i) * etad(i,1) * ecdo(i,1,n) + & * grav / dp endif enddo enddo @@ -2044,9 +2084,7 @@ subroutine samfdeepcnv_run (im,km,itc,ntc,cliq,cp,cvap, & dv1h = heo(i,k) dv2h = .5 * (heo(i,k) + heo(i,k-1)) dv3h = heo(i,k-1) - dv1q = qo(i,k) dv2q = .5 * (qo(i,k) + qo(i,k-1)) - dv3q = qo(i,k-1) c tem = 0.5 * (xlamue(i,k)+xlamue(i,k-1)) tem1 = 0.5 * (xlamud(i,k)+xlamud(i,k-1)) @@ -2058,6 +2096,8 @@ subroutine samfdeepcnv_run (im,km,itc,ntc,cliq,cp,cvap, & ptem = xlamde ptem1 = xlamdd endif + + factor = grav / dp cj dellah(i,k) = dellah(i,k) + & ((aup*eta(i,k)-adw*edto(i)*etad(i,k))*dv1h @@ -2065,29 +2105,27 @@ subroutine samfdeepcnv_run (im,km,itc,ntc,cliq,cp,cvap, & & - (aup*tem*eta(i,k-1)+adw*edto(i)*ptem*etad(i,k))*dv2h*dz & + aup*tem1*eta(i,k-1)*.5*(hcko(i,k)+hcko(i,k-1))*dz & + adw*edto(i)*ptem1*etad(i,k)*.5*(hcdo(i,k)+hcdo(i,k-1))*dz - & ) *grav/dp + & ) * factor cj dellaq(i,k) = dellaq(i,k) + - & ((aup*eta(i,k)-adw*edto(i)*etad(i,k))*dv1q - & - (aup*eta(i,k-1)-adw*edto(i)*etad(i,k-1))*dv3q - & - (aup*tem*eta(i,k-1)+adw*edto(i)*ptem*etad(i,k))*dv2q*dz + & (- (aup*tem*eta(i,k-1)+adw*edto(i)*ptem*etad(i,k))*dv2q*dz & + aup*tem1*eta(i,k-1)*.5*(qrcko(i,k)+qcko(i,k-1))*dz & + adw*edto(i)*ptem1*etad(i,k)*.5*(qrcdo(i,k)+qcdo(i,k-1))*dz - & ) *grav/dp + & ) * factor cj tem1=eta(i,k)*(uo(i,k)-ucko(i,k)) tem2=eta(i,k-1)*(uo(i,k-1)-ucko(i,k-1)) ptem1=etad(i,k)*(uo(i,k)-ucdo(i,k)) ptem2=etad(i,k-1)*(uo(i,k-1)-ucdo(i,k-1)) dellau(i,k) = dellau(i,k) + - & (aup*(tem1-tem2)-adw*edto(i)*(ptem1-ptem2))*grav/dp + & (aup*(tem1-tem2)-adw*edto(i)*(ptem1-ptem2))*factor cj tem1=eta(i,k)*(vo(i,k)-vcko(i,k)) tem2=eta(i,k-1)*(vo(i,k-1)-vcko(i,k-1)) ptem1=etad(i,k)*(vo(i,k)-vcdo(i,k)) ptem2=etad(i,k-1)*(vo(i,k-1)-vcdo(i,k-1)) dellav(i,k) = dellav(i,k) + - & (aup*(tem1-tem2)-adw*edto(i)*(ptem1-ptem2))*grav/dp + & (aup*(tem1-tem2)-adw*edto(i)*(ptem1-ptem2))*factor cj endif enddo @@ -2103,10 +2141,10 @@ subroutine samfdeepcnv_run (im,km,itc,ntc,cliq,cp,cvap, & if(k > jmin(i)) adw = 0. dp = 1000. * del(i,k) cj - tem1=eta(i,k)*(ctro(i,k,n)-ecko(i,k,n)) - tem2=eta(i,k-1)*(ctro(i,k-1,n)-ecko(i,k-1,n)) - ptem1=etad(i,k)*(ctro(i,k,n)-ecdo(i,k,n)) - ptem2=etad(i,k-1)*(ctro(i,k-1,n)-ecdo(i,k-1,n)) + tem1 = -eta(i,k) * ecko(i,k,n) + tem2 = -eta(i,k-1) * ecko(i,k-1,n) + ptem1 = -etad(i,k) * ecdo(i,k,n) + ptem2 = -etad(i,k-1) * ecdo(i,k-1,n) dellae(i,k,n) = dellae(i,k,n) + & (aup*(tem1-tem2)-adw*edto(i)*(ptem1-ptem2))*grav/dp cj @@ -2122,21 +2160,15 @@ subroutine samfdeepcnv_run (im,km,itc,ntc,cliq,cp,cvap, & if(cnvflg(i)) then indx = ktcon(i) dp = 1000. * del(i,indx) - dv1h = heo(i,indx-1) - dellah(i,indx) = eta(i,indx-1) * - & (hcko(i,indx-1) - dv1h) * grav / dp - dv1q = qo(i,indx-1) - dellaq(i,indx) = eta(i,indx-1) * - & (qcko(i,indx-1) - dv1q) * grav / dp - dellau(i,indx) = eta(i,indx-1) * - & (ucko(i,indx-1) - uo(i,indx-1)) * grav / dp - dellav(i,indx) = eta(i,indx-1) * - & (vcko(i,indx-1) - vo(i,indx-1)) * grav / dp + tem = eta(i,indx-1) * grav / dp + dellah(i,indx) = tem * (hcko(i,indx-1) - heo(i,indx-1)) + dellaq(i,indx) = tem * qcko(i,indx-1) + dellau(i,indx) = tem * (ucko(i,indx-1) - uo(i,indx-1)) + dellav(i,indx) = tem * (vcko(i,indx-1) - vo(i,indx-1)) c c cloud water c - dellal(i,indx) = eta(i,indx-1) * - & qlko_ktcon(i) * grav / dp + dellal(i,indx) = tem * qlko_ktcon(i) endif enddo if (.not.hwrf_samfdeep) then @@ -2146,10 +2178,179 @@ subroutine samfdeepcnv_run (im,km,itc,ntc,cliq,cp,cvap, & indx = ktcon(i) dp = 1000. * del(i,indx) dellae(i,indx,n) = eta(i,indx-1) * - & (ecko(i,indx-1,n) - ctro(i,indx-1,n)) * grav / dp + & ecko(i,indx-1,n) * grav / dp + endif + enddo + enddo + endif +! +! compute change rates due to environmental subsidence & uplift +! using a positive definite TVD flux-limiter scheme +! +! for moisture +! + do k=1,km1 + do i=1,im + if(cnvflg(i) .and. k <= ktcon(i)) then + q_diff(i,k) = q1(i,k) - q1(i,k+1) + endif + enddo + enddo + do i=1,im + if(cnvflg(i)) then + if(q1(i,1) >= 0.) then + q_diff(i,0) = max(0.,2.*q1(i,1)-q1(i,2))- + & q1(i,1) + else + q_diff(i,0) = min(0.,2.*q1(i,1)-q1(i,2))- + & q1(i,1) + endif endif enddo +! + flxtvd = 0. + do k = 1, km1 + do i = 1, im + if(cnvflg(i) .and. k < ktcon(i)) then + tem = 0. + if(k >= kb(i)) tem = eta(i,k) + if(k <= jmin(i)) then + tem = tem - edto(i) * etad(i,k) + endif + if(tem > 0.) then + rrkp = 0. + if(abs(q_diff(i,k)) > 1.e-22) + & rrkp = q_diff(i,k+1) / q_diff(i,k) + phkp = (rrkp+abs(rrkp)) / (1.+abs(rrkp)) + tem1 = q1(i,k+1) + + & phkp*(qo(i,k)-q1(i,k+1)) + flxtvd(i,k) = tem * tem1 + elseif(tem < 0.) then + rrkp = 0. + if(abs(q_diff(i,k)) > 1.e-22) + & rrkp = q_diff(i,k-1) / q_diff(i,k) + phkp = (rrkp+abs(rrkp)) / (1.+abs(rrkp)) + tem1 = q1(i,k) + + & phkp*(qo(i,k)-q1(i,k)) + flxtvd(i,k) = tem * tem1 + else + tem1 = qo(i,k) + flxtvd(i,k) = 0. + endif +! +! subtract the double counting change rates at jmin+1 & kb beforehand +! + if(k == jmin(i)) then + dp = 1000. * del(i,k+1) + dellaq(i,k+1) = dellaq(i,k+1) - + & edto(i) * etad(i,k) * tem1 * grav/dp + endif + if(k == kb(i)) then + dp = 1000. * del(i,k) + dellaq(i,k) = dellaq(i,k) - + & eta(i,k) * tem1 * grav/dp + endif +! + endif + enddo + enddo +! + do k=1,km1 + do i=1,im + if(cnvflg(i) .and. k <= ktcon(i)) then + dp = 1000. * del(i,k) + dellaq(i,k) = dellaq(i,k) + + & (flxtvd(i,k) - flxtvd(i,k-1)) * grav/dp + endif + enddo + enddo +! +! for tracers including TKE & ozone +! + if (.not.hwrf_samfdeep) then +! + do n=1,ntr + do k=1,km1 + do i=1,im + if(cnvflg(i) .and. k <= ktcon(i)) then + e_diff(i,k,n) = ctr(i,k,n) - ctr(i,k+1,n) + endif + enddo + enddo + do i=1,im + if(cnvflg(i)) then + if(ctr(i,1,n) >= 0.) then + e_diff(i,0,n) = max(0.,2.*ctr(i,1,n)-ctr(i,2,n))- + & ctr(i,1,n) + else + e_diff(i,0,n) = min(0.,2.*ctr(i,1,n)-ctr(i,2,n))- + & ctr(i,1,n) + endif + endif + enddo + enddo +! + do n=1,ntr +! + flxtvd = 0. + do k = 1, km1 + do i = 1, im + if(cnvflg(i) .and. k < ktcon(i)) then + tem = 0. + if(k >= kb(i)) tem = eta(i,k) + if(k <= jmin(i)) then + tem = tem - edto(i) * etad(i,k) + endif + if(tem > 0.) then + rrkp = 0. + if(abs(e_diff(i,k,n)) > 1.e-22) + & rrkp = e_diff(i,k+1,n) / e_diff(i,k,n) + phkp = (rrkp+abs(rrkp)) / (1.+abs(rrkp)) + tem1 = ctr(i,k+1,n) + + & phkp*(ctro(i,k,n)-ctr(i,k+1,n)) + flxtvd(i,k) = tem * tem1 + elseif(tem < 0.) then + rrkp = 0. + if(abs(e_diff(i,k,n)) > 1.e-22) + & rrkp = e_diff(i,k-1,n) / e_diff(i,k,n) + phkp = (rrkp+abs(rrkp)) / (1.+abs(rrkp)) + tem1 = ctr(i,k,n) + + & phkp*(ctro(i,k,n)-ctr(i,k,n)) + flxtvd(i,k) = tem * tem1 + else + tem1 = ctro(i,k,n) + flxtvd(i,k) = 0. + endif +! +! subtract the double counting change rates at jmin+1 & kb beforehand +! + if(k == jmin(i)) then + dp = 1000. * del(i,k+1) + dellae(i,k+1,n) = dellae(i,k+1,n) - + & edto(i)*etad(i,k) * tem1 * grav/dp + endif + if(k == kb(i)) then + dp = 1000. * del(i,k) + dellae(i,k,n) = dellae(i,k,n) - + & eta(i,k) * tem1 * grav/dp + endif +! + endif + enddo + enddo +! + do k=1,km1 + do i=1,im + if(cnvflg(i) .and. k <= ktcon(i)) then + dp = 1000. * del(i,k) + dellae(i,k,n) = dellae(i,k,n) + + & (flxtvd(i,k) - flxtvd(i,k-1)) * grav/dp + endif + enddo + enddo +! enddo +! endif c c------- final changed variable per unit mass flux @@ -2431,19 +2632,19 @@ subroutine samfdeepcnv_run (im,km,itc,ntc,cliq,cp,cvap, & do i = 1, im if (asqecflg(i) .and. k < jmin(i)) then gamma = el2orc * qeso(i,k) / to(i,k)**2 - dhh=hcdo(i,k) - dt= to(i,k) - dg= gamma - dh= heso(i,k) - dz=-1.*(zo(i,k+1)-zo(i,k)) + dhh = hcdo(i,k) + dt = to(i,k) + dg = gamma + dh = heso(i,k) + dz = zo(i,k) - zo(i,k+1) ! xaa0(i)=xaa0(i)+edtx(i)*dz*etad(i,k) - xaa0(i)=xaa0(i)+edtx(i)*dz - & *(grav/(cp*dt))*((dhh-dh)/(1.+dg)) - & *(1.+fv*cp*dg*dt/hvap) - val=0. + xaa0(i) = xaa0(i)+edtx(i)*dz + & *(grav/(cp*dt))*((dhh-dh)/(1.+dg)) + & *(1.+fv*cp*dg*dt/hvap) + val = 0. ! xaa0(i)=xaa0(i)+edtx(i)*dz*etad(i,k) - xaa0(i)=xaa0(i)+edtx(i)*dz - & *grav*fv*max(val,(qeso(i,k)-qo(i,k))) + xaa0(i) = xaa0(i)+edtx(i)*dz + & *grav*fv*max(val,(qeso(i,k)-qo(i,k))) endif enddo enddo @@ -2680,13 +2881,13 @@ subroutine samfdeepcnv_run (im,km,itc,ntc,cliq,cp,cvap, & endif !> - Transport aerosols if present - - if (do_aerosols) - & call samfdeepcnv_aerosols(im, im, km, itc, ntc, ntr, delt, - & xlamde, xlamdd, cnvflg, jmin, kb, kmax, kd94, ktcon, fscav, - & edto, xlamd, xmb, c0t, eta, etad, zi, xlamue, xlamud, delp, - & qtr, qaero) - +! +! if (do_aerosols) +! & call samfdeepcnv_aerosols(im, im, km, itc, ntc, ntr, delt, +! & xlamde, xlamdd, cnvflg, jmin, kb, kmax, kd94, ktcon, fscav, +! & edto, xlamd, xmb, c0t, eta, etad, zi, xlamue, xlamud, delp, +! & qtr, qaero) +! c c restore to,qo,uo,vo to t1,q1,u1,v1 in case convection stops c @@ -2706,10 +2907,11 @@ subroutine samfdeepcnv_run (im,km,itc,ntc,cliq,cp,cvap, & enddo if (.not.hwrf_samfdeep) then do n = 1, ntr + kk = n+2 do k = 1, km do i = 1, im if (cnvflg(i) .and. k <= kmax(i)) then - ctro(i,k,n) = ctr(i,k,n) + ctro(i,k,n) = qtr(i,k,kk) endif enddo enddo @@ -2744,40 +2946,184 @@ subroutine samfdeepcnv_run (im,km,itc,ntc,cliq,cp,cvap, & do i = 1, im if (cnvflg(i) .and. k <= kmax(i)) then if(k <= ktcon(i)) then + tem2 = xmb(i) * dt2 dellat = (dellah(i,k) - hvap * dellaq(i,k)) / cp - t1(i,k) = t1(i,k) + dellat * xmb(i) * dt2 - q1(i,k) = q1(i,k) + dellaq(i,k) * xmb(i) * dt2 -! tem = 1./rcs(i) -! u1(i,k) = u1(i,k) + dellau(i,k) * xmb(i) * dt2 * tem -! v1(i,k) = v1(i,k) + dellav(i,k) * xmb(i) * dt2 * tem - u1(i,k) = u1(i,k) + dellau(i,k) * xmb(i) * dt2 - v1(i,k) = v1(i,k) + dellav(i,k) * xmb(i) * dt2 + t1(i,k) = t1(i,k) + tem2 * dellat + q1(i,k) = q1(i,k) + tem2 * dellaq(i,k) +! tem = tem2 / rcs(i) +! u1(i,k) = u1(i,k) + dellau(i,k) * tem +! v1(i,k) = v1(i,k) + dellav(i,k) * tem + u1(i,k) = u1(i,k) + tem2 * dellau(i,k) + v1(i,k) = v1(i,k) + tem2 * dellav(i,k) dp = 1000. * del(i,k) - delhbar(i) = delhbar(i) + dellah(i,k)*xmb(i)*dp/grav - delqbar(i) = delqbar(i) + dellaq(i,k)*xmb(i)*dp/grav - deltbar(i) = deltbar(i) + dellat*xmb(i)*dp/grav - delubar(i) = delubar(i) + dellau(i,k)*xmb(i)*dp/grav - delvbar(i) = delvbar(i) + dellav(i,k)*xmb(i)*dp/grav + tem = xmb(i) * dp / grav + delhbar(i) = delhbar(i) + tem * dellah(i,k) + delqbar(i) = delqbar(i) + tem * dellaq(i,k) + deltbar(i) = deltbar(i) + tem * dellat + delubar(i) = delubar(i) + tem * dellau(i,k) + delvbar(i) = delvbar(i) + tem * dellav(i,k) + endif + endif + enddo + enddo +! +! Negative moisture is set to zero after borrowing it from +! positive values within the mass-flux transport layers +! + do i = 1,im + tsumn(i) = 0. + tsump(i) = 0. + rtnp(i) = 1. + enddo + do k = 1,km + do i = 1,im + if(cnvflg(i) .and. k <= ktcon(i)) then + if(q1(i,k) < 0.) tsumn(i) = tsumn(i) + q1(i,k) + if(q1(i,k) > 0.) tsump(i) = tsump(i) + q1(i,k) + endif + enddo + enddo + do i = 1,im + if(cnvflg(i)) then + if(tsump(i) > 0. .and. tsumn(i) < 0.) then + if(tsump(i) > abs(tsumn(i))) then + rtnp(i) = tsumn(i) / tsump(i) + else + rtnp(i) = tsump(i) / tsumn(i) + endif + endif + endif + enddo + do k = 1,km + do i = 1,im + if(cnvflg(i) .and. k <= ktcon(i)) then + if(rtnp(i) < 0.) then + if(tsump(i) > abs(tsumn(i))) then + if(q1(i,k) < 0.) q1(i,k) = 0. + if(q1(i,k) > 0.) q1(i,k) = (1.+rtnp(i))*q1(i,k) + else + if(q1(i,k) < 0.) q1(i,k) = (1.+rtnp(i))*q1(i,k) + if(q1(i,k) > 0.) q1(i,k) = 0. + endif endif endif enddo enddo +! if (.not.hwrf_samfdeep) then do n = 1, ntr - kk = n+2 +! do k = 1, km do i = 1, im if (cnvflg(i) .and. k <= kmax(i)) then if(k <= ktcon(i)) then ctr(i,k,n) = ctr(i,k,n)+dellae(i,k,n)*xmb(i)*dt2 + dp = 1000. * del(i,k) delebar(i,n)=delebar(i,n)+dellae(i,k,n)*xmb(i)*dp/grav - qtr(i,k,kk) = ctr(i,k,n) endif endif enddo enddo +! +! Negative TKE, ozone, and aerosols are set to zero after borrowing them +! from positive values within the mass-flux transport layers +! + do i = 1,im + tsumn(i) = 0. + tsump(i) = 0. + rtnp(i) = 1. + enddo + do k = 1,km + do i = 1,im + if(cnvflg(i) .and. k <= ktcon(i)) then + if(ctr(i,k,n) < 0.) tsumn(i) = tsumn(i) + ctr(i,k,n) + if(ctr(i,k,n) > 0.) tsump(i) = tsump(i) + ctr(i,k,n) + endif + enddo + enddo + do i = 1,im + if(cnvflg(i)) then + if(tsump(i) > 0. .and. tsumn(i) < 0.) then + if(tsump(i) > abs(tsumn(i))) then + rtnp(i) = tsumn(i) / tsump(i) + else + rtnp(i) = tsump(i) / tsumn(i) + endif + endif + endif + enddo + do k = 1,km + do i = 1,im + if(cnvflg(i) .and. k <= ktcon(i)) then + if(rtnp(i) < 0.) then + if(tsump(i) > abs(tsumn(i))) then + if(ctr(i,k,n)<0.) ctr(i,k,n)=0. + if(ctr(i,k,n)>0.) ctr(i,k,n)=(1.+rtnp(i))*ctr(i,k,n) + else + if(ctr(i,k,n)<0.) ctr(i,k,n)=(1.+rtnp(i))*ctr(i,k,n) + if(ctr(i,k,n)>0.) ctr(i,k,n)=0. + endif + endif + endif + enddo + enddo +! + kk = n+2 + do k = 1, km + do i = 1, im + if(cnvflg(i) .and. k <= ktcon(i)) then + qtr(i,k,kk) = ctr(i,k,n) + endif + enddo + enddo +! enddo +! + if (do_aerosols) then +! + do n = 1, ntc +! +! convert wet deposition to total mass deposited over dt2 and dp + do k = 2, km1 + do i = 1, im + if (cnvflg(i)) then + if(k > kb(i) .and. k < ktcon(i)) then + dp = 1000. * del(i,k) + wet_dep(i,k,n) = chem_pw(i,k,n)*grav/dp + wet_dep(i,k,n) = wet_dep(i,k,n)*xmb(i)*dt2*dp + endif + endif + enddo + enddo +! + kk = n + itc - 1 + do k = 2, km1 + do i = 1, im + if (cnvflg(i)) then + if(k > kb(i) .and. k < ktcon(i)) then + dp = 1000. * del(i,k) + if (qtr(i,k,kk) < 0.) then +! borrow negative mass from wet deposition + tem = -qtr(i,k,kk)*dp + if(wet_dep(i,k,n) >= tem) then + wet_dep(i,k,n) = wet_dep(i,k,n) - tem + qtr(i,k,kk) = 0. + else + wet_dep(i,k,n) = 0. + qtr(i,k,kk) = qtr(i,k,kk)+wet_dep(i,k,n)/dp + endif + endif + endif + endif + enddo + enddo +! + enddo +! + endif +! endif +! !> - Recalculate saturation specific humidity using the updated temperature. do k = 1, km do i = 1, im @@ -2837,26 +3183,28 @@ subroutine samfdeepcnv_run (im,km,itc,ntc,cliq,cp,cvap, & qcond(i) = evef * (q1(i,k) - qeso(i,k)) & / (1. + el2orc * qeso(i,k) / t1(i,k)**2) dp = 1000. * del(i,k) + tem = grav / dp + tem1 = dp / grav if(rn(i) > 0. .and. qcond(i) < 0.) then qevap(i) = -qcond(i) * (1.-exp(-.32*sqrt(dt2*rn(i)))) - qevap(i) = min(qevap(i), rn(i)*1000.*grav/dp) - delq2(i) = delqev(i) + .001 * qevap(i) * dp / grav + qevap(i) = min(qevap(i), rn(i)*1000.*tem) + delq2(i) = delqev(i) + .001 * qevap(i) * tem1 endif if(rn(i) > 0. .and. qcond(i) < 0. .and. & delq2(i) > rntot(i)) then - qevap(i) = 1000.* grav * (rntot(i) - delqev(i)) / dp + qevap(i) = 1000.* tem * (rntot(i) - delqev(i)) flg(i) = .false. endif if(rn(i) > 0. .and. qevap(i) > 0.) then q1(i,k) = q1(i,k) + qevap(i) t1(i,k) = t1(i,k) - elocp * qevap(i) - rn(i) = rn(i) - .001 * qevap(i) * dp / grav + rn(i) = rn(i) - .001 * qevap(i) * tem1 deltv(i) = - elocp*qevap(i)/dt2 delq(i) = + qevap(i)/dt2 - delqev(i) = delqev(i) + .001*dp*qevap(i)/grav + delqev(i) = delqev(i) + .001 * qevap(i) * tem1 endif - delqbar(i) = delqbar(i) + delq(i)*dp/grav - deltbar(i) = deltbar(i) + deltv(i)*dp/grav + delqbar(i) = delqbar(i) + delq(i) * tem1 + deltbar(i) = deltbar(i) + deltv(i) * tem1 endif endif enddo @@ -2869,11 +3217,14 @@ subroutine samfdeepcnv_run (im,km,itc,ntc,cliq,cp,cvap, & enddo endif cj -! do i = 1, im +! do i = 1, 4 ! if(me == 31 .and. cnvflg(i)) then ! if(cnvflg(i)) then +! if(i==1) print*,'ntr:ntk= ',ntr,ntk ! print *, ' deep delhbar, delqbar, deltbar = ', ! & delhbar(i),hvap*delqbar(i),cp*deltbar(i) +! print *, ' deep delebar ozone = ',delebar(i,1) +! print *, ' deep delebar tke = ',delebar(i,2) ! print *, ' deep delubar, delvbar = ',delubar(i),delvbar(i) ! print *, ' precip =', hvap*rn(i)*1000./dt2 ! print*,'pdif= ',pfld(i,kbcon(i))-pfld(i,ktcon(i)) @@ -2975,27 +3326,38 @@ subroutine samfdeepcnv_run (im,km,itc,ntc,cliq,cp,cvap, & do i = 1, im if(cnvflg(i) .and. rn(i) <= 0.) then if (k <= kmax(i)) then - ctr(i,k,n)= ctro(i,k,n) - qtr(i,k,kk)= ctr(i,k,n) + qtr(i,k,kk)= ctro(i,k,n) endif endif enddo enddo enddo - -!> - Store aerosol concentrations if present if (do_aerosols) then do n = 1, ntc - kk = n + itc - 1 - do k = 1, km + do k = 2, km1 do i = 1, im - if(cnvflg(i) .and. rn(i) > 0.) then - if (k <= kmax(i)) qtr(i,k,kk) = qaero(i,k,n) + if(cnvflg(i) .and. rn(i) <= 0.) then + if (k <= ktcon(i)) then + wet_dep(i,k,n) = 0. + endif endif enddo enddo enddo - endif + endif +!> - Store aerosol concentrations if present +! if (do_aerosols) then +! do n = 1, ntc +! kk = n + itc - 1 +! do k = 1, km +! do i = 1, im +! if(cnvflg(i) .and. rn(i) > 0.) then +! if (k <= kmax(i)) qtr(i,k,kk) = qaero(i,k,n) +! endif +! enddo +! enddo +! enddo +! endif endif ! ! hchuang code change diff --git a/physics/samfshalcnv.f b/physics/samfshalcnv.f index c314809cc..846fb30c1 100644 --- a/physics/samfshalcnv.f +++ b/physics/samfshalcnv.f @@ -104,7 +104,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & & dq, dqsdp, dqsdt, dt, & dt2, dtmax, dtmin, dxcrt, & dv1h, dv2h, dv3h, - & dv1q, dv2q, dv3q, + & dv2q, & dz, dz1, e1, & el2orc, elocp, aafac, cm, & es, etah, h1, @@ -188,7 +188,12 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & & uo(im,km), vo(im,km), qeso(im,km), & ctr(im,km,ntr), ctro(im,km,ntr) ! for aerosol transport - real(kind=kind_phys) qaero(im,km,ntc) +! real(kind=kind_phys) qaero(im,km,ntc) +c variables for tracer wet deposition, + real(kind=kind_phys), dimension(im,km,ntc) :: chem_c, chem_pw, + & wet_dep + real(kind=kind_phys), parameter :: escav = 0.8 ! wet scavenging efficiency +! ! for updraft velocity calculation real(kind=kind_phys) wu2(im,km), buo(im,km), drag(im,km) real(kind=kind_phys) wc(im), scaldfunc(im), sigmagfm(im) @@ -208,6 +213,14 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & & zi(im,km), pwo(im,km), c0t(im,km), & sumx(im), tx1(im), cnvwt(im,km) &, rhbar(im) +! +! variables for Total Variation Diminishing (TVD) flux-limiter scheme +! on environmental subsidence and uplifting +! + real(kind=kind_phys) q_diff(im,0:km-1), e_diff(im,0:km-1,ntr), + & flxtvd(im,km-1) + real(kind=kind_phys) rrkp, phkp + real(kind=kind_phys) tsumn(im), tsump(im), rtnp(im) ! logical do_aerosols, totflg, cnvflg(im), flg(im) ! @@ -252,6 +265,11 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & c initialize arrays c !> - Initialize column-integrated and other single-value-per-column variable arrays. +! + chem_c = 0. + chem_pw = 0. + wet_dep = 0. +! if(hwrf_samfshal) then do i=1,im cnvflg(i) = .true. @@ -966,7 +984,12 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & enddo if (.not.hwrf_samfshal) then - do n = 1, ntr + if (do_aerosols) then + kk = itc -3 + else + kk = ntr + endif + do n = 1, kk do k = 2, km1 do i = 1, im if (cnvflg(i)) then @@ -981,6 +1004,28 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & enddo enddo enddo + if (do_aerosols) then + do n = 1, ntc + kk = n + itc -3 + do k = 2, km1 + do i = 1, im + if (cnvflg(i)) then + if(k > kb(i) .and. k < kmax(i)) then + dz = zi(i,k) - zi(i,k-1) + tem = 0.25 * (xlamue(i,k)+xlamue(i,k-1)) * dz + factor = 1. + tem + ecko(i,k,kk) = ((1. - tem) * ecko(i,k-1,kk) + tem * + & (ctro(i,k,kk) + ctro(i,k-1,kk))) / factor + chem_c(i,k,n) = escav * fscav(n) * ecko(i,k,kk) + tem = chem_c(i,k,n) / (1. + c0t(i,k) * dz) + chem_pw(i,k,n) = c0t(i,k) * dz * tem * eta(i,k-1) + ecko(i,k,kk) = tem + ecko(i,k,kk) - chem_c(i,k,n) + endif + endif + enddo + enddo + enddo + endif endif c c taking account into convection inhibition due to existence of @@ -1520,32 +1565,31 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & dv1h = heo(i,k) dv2h = .5 * (heo(i,k) + heo(i,k-1)) dv3h = heo(i,k-1) - dv1q = qo(i,k) dv2q = .5 * (qo(i,k) + qo(i,k-1)) - dv3q = qo(i,k-1) c tem = 0.5 * (xlamue(i,k)+xlamue(i,k-1)) tem1 = xlamud(i) + + factor = grav / dp cj dellah(i,k) = dellah(i,k) + & ( eta(i,k)*dv1h - eta(i,k-1)*dv3h & - tem*eta(i,k-1)*dv2h*dz & + tem1*eta(i,k-1)*.5*(hcko(i,k)+hcko(i,k-1))*dz - & ) *grav/dp + & ) * factor cj dellaq(i,k) = dellaq(i,k) + - & ( eta(i,k)*dv1q - eta(i,k-1)*dv3q - & - tem*eta(i,k-1)*dv2q*dz + & ( - tem*eta(i,k-1)*dv2q*dz & + tem1*eta(i,k-1)*.5*(qrcko(i,k)+qcko(i,k-1))*dz - & ) *grav/dp + & ) * factor cj tem1=eta(i,k)*(uo(i,k)-ucko(i,k)) tem2=eta(i,k-1)*(uo(i,k-1)-ucko(i,k-1)) - dellau(i,k) = dellau(i,k) + (tem1-tem2) * grav/dp + dellau(i,k) = dellau(i,k) + (tem1-tem2) * factor cj tem1=eta(i,k)*(vo(i,k)-vcko(i,k)) tem2=eta(i,k-1)*(vo(i,k-1)-vcko(i,k-1)) - dellav(i,k) = dellav(i,k) + (tem1-tem2) * grav/dp + dellav(i,k) = dellav(i,k) + (tem1-tem2) * factor cj endif endif @@ -1559,8 +1603,8 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & if(k > kb(i) .and. k < ktcon(i)) then dp = 1000. * del(i,k) cj - tem1=eta(i,k)*(ctro(i,k,n)-ecko(i,k,n)) - tem2=eta(i,k-1)*(ctro(i,k-1,n)-ecko(i,k-1,n)) + tem1 = -eta(i,k) * ecko(i,k,n) + tem2 = -eta(i,k-1) * ecko(i,k-1,n) dellae(i,k,n) = dellae(i,k,n) + (tem1-tem2) * grav/dp cj endif @@ -1576,21 +1620,15 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & if(cnvflg(i)) then indx = ktcon(i) dp = 1000. * del(i,indx) - dv1h = heo(i,indx-1) - dellah(i,indx) = eta(i,indx-1) * - & (hcko(i,indx-1) - dv1h) * grav / dp - dv1q = qo(i,indx-1) - dellaq(i,indx) = eta(i,indx-1) * - & (qcko(i,indx-1) - dv1q) * grav / dp - dellau(i,indx) = eta(i,indx-1) * - & (ucko(i,indx-1) - uo(i,indx-1)) * grav / dp - dellav(i,indx) = eta(i,indx-1) * - & (vcko(i,indx-1) - vo(i,indx-1)) * grav / dp + tem = eta(i,indx-1) * grav / dp + dellah(i,indx) = tem * (hcko(i,indx-1) - heo(i,indx-1)) + dellaq(i,indx) = tem * qcko(i,indx-1) + dellau(i,indx) = tem * (ucko(i,indx-1) - uo(i,indx-1)) + dellav(i,indx) = tem * (vcko(i,indx-1) - vo(i,indx-1)) c c cloud water c - dellal(i,indx) = eta(i,indx-1) * - & qlko_ktcon(i) * grav / dp + dellal(i,indx) = tem * qlko_ktcon(i) endif enddo if (.not.hwrf_samfshal) then @@ -1600,12 +1638,125 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & indx = ktcon(i) dp = 1000. * del(i,indx) dellae(i,indx,n) = eta(i,indx-1) * - & (ecko(i,indx-1,n) - ctro(i,indx-1,n)) * grav / dp + & ecko(i,indx-1,n) * grav / dp endif enddo enddo endif ! +! compute change rates due to environmental subsidence & uplift +! using a positive definite TVD flux-limiter scheme +! +! for moisture +! + do k=1,km1 + do i=1,im + if(cnvflg(i) .and. k <= ktcon(i)) then + q_diff(i,k) = q1(i,k) - q1(i,k+1) + endif + enddo + enddo + do i=1,im + if(cnvflg(i)) then + if(q1(i,1) >= 0.) then + q_diff(i,0) = max(0.,2.*q1(i,1)-q1(i,2))- + & q1(i,1) + else + q_diff(i,0) = min(0.,2.*q1(i,1)-q1(i,2))- + & q1(i,1) + endif + endif + enddo +! + flxtvd = 0. + do k = 1, km1 + do i = 1, im + if(cnvflg(i) .and. + & (k >= kb(i) .and. k < ktcon(i))) then + if(eta(i,k) > 0.) then + rrkp = 0. + if(abs(q_diff(i,k)) > 1.e-22) + & rrkp = q_diff(i,k+1) / q_diff(i,k) + phkp = (rrkp+abs(rrkp)) / (1.+abs(rrkp)) + tem1 = q1(i,k+1) + + & phkp*(qo(i,k)-q1(i,k+1)) + flxtvd(i,k) = eta(i,k) * tem1 + endif + endif + enddo + enddo +! + do k = 2, km1 + do i = 1, im + if(cnvflg(i) .and. + & (k > kb(i) .and. k <= ktcon(i))) then + dp = 1000. * del(i,k) + dellaq(i,k) = dellaq(i,k) + + & (flxtvd(i,k) - flxtvd(i,k-1)) * grav/dp + endif + enddo + enddo +! +! for tracers including TKE & ozone +! + if (.not.hwrf_samfshal) then +! + do n=1,ntr + do k=1,km1 + do i=1,im + if(cnvflg(i) .and. k <= ktcon(i)) then + e_diff(i,k,n) = ctr(i,k,n) - ctr(i,k+1,n) + endif + enddo + enddo + do i=1,im + if(cnvflg(i)) then + if(ctr(i,1,n) >= 0.) then + e_diff(i,0,n) = max(0.,2.*ctr(i,1,n)-ctr(i,2,n))- + & ctr(i,1,n) + else + e_diff(i,0,n) = min(0.,2.*ctr(i,1,n)-ctr(i,2,n))- + & ctr(i,1,n) + endif + endif + enddo + enddo +! + do n=1,ntr +! + flxtvd = 0. + do k= 1, km1 + do i = 1, im + if(cnvflg(i) .and. + & (k >= kb(i) .and. k < ktcon(i))) then + if(eta(i,k) > 0.) then + rrkp = 0. + if(abs(e_diff(i,k,n)) > 1.e-22) + & rrkp = e_diff(i,k+1,n) / e_diff(i,k,n) + phkp = (rrkp+abs(rrkp)) / (1.+abs(rrkp)) + tem1 = ctr(i,k+1,n) + + & phkp*(ctro(i,k,n)-ctr(i,k+1,n)) + flxtvd(i,k) = eta(i,k) * tem1 + endif + endif + enddo + enddo +! + do k = 2, km1 + do i = 1, im + if(cnvflg(i) .and. + & (k > kb(i) .and. k <= ktcon(i))) then + dp = 1000. * del(i,k) + dellae(i,k,n) = dellae(i,k,n) + + & (flxtvd(i,k) - flxtvd(i,k-1)) * grav/dp + endif + enddo + enddo +! + enddo +! + endif +! ! compute convective turn-over time ! !> - Following Bechtold et al. (2008) \cite bechtold_et_al_2008, calculate the convective turnover time using the mean updraft velocity (wc) and the cloud depth. It is also proportional to the grid size (gdx). @@ -1694,15 +1845,15 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & ! !> - Transport aerosols if present ! - if (.not.hwrf_samfshal) then - if (do_aerosols) - & call samfshalcnv_aerosols(im, im, km, itc, ntc, ntr, delt, -! & xlamde, xlamdd, cnvflg, jmin, kb, kmax, kbcon, ktcon, fscav, - & cnvflg, kb, kmax, kbcon, ktcon, fscav, -! & edto, xlamd, xmb, c0t, eta, etad, zi, xlamue, xlamud, delp, - & xmb, c0t, eta, zi, xlamue, xlamud, delp, - & qtr, qaero) - endif +! if (.not.hwrf_samfshal) then +! if (do_aerosols) +! & call samfshalcnv_aerosols(im, im, km, itc, ntc, ntr, delt, +!! & xlamde, xlamdd, cnvflg, jmin, kb, kmax, kbcon, ktcon, fscav, +! & cnvflg, kb, kmax, ktcon, fscav, +!! & edto, xlamd, xmb, c0t, eta, etad, zi, xlamue, xlamud, delp, +! & xmb, c0t, eta, zi, xlamue, xlamud, delp, +! & qtr, qaero) +! endif ! !> ## For the "feedback control", calculate updated values of the state variables by multiplying the cloud base mass flux and the tendencies calculated per unit cloud base mass flux from the static control. !! - Recalculate saturation specific humidity. @@ -1752,30 +1903,183 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & u1(i,k) = u1(i,k) + dellau(i,k) * xmb(i) * dt2 v1(i,k) = v1(i,k) + dellav(i,k) * xmb(i) * dt2 dp = 1000. * del(i,k) - delhbar(i) = delhbar(i) + dellah(i,k)*xmb(i)*dp/grav - delqbar(i) = delqbar(i) + dellaq(i,k)*xmb(i)*dp/grav - deltbar(i) = deltbar(i) + dellat*xmb(i)*dp/grav - delubar(i) = delubar(i) + dellau(i,k)*xmb(i)*dp/grav - delvbar(i) = delvbar(i) + dellav(i,k)*xmb(i)*dp/grav + tem = xmb(i) * dp / grav + delhbar(i) = delhbar(i) + tem * dellah(i,k) + delqbar(i) = delqbar(i) + tem * dellaq(i,k) + deltbar(i) = deltbar(i) + tem * dellat + delubar(i) = delubar(i) + tem * dellau(i,k) + delvbar(i) = delvbar(i) + tem * dellav(i,k) + endif + endif + enddo + enddo +! +! Negative moisture is set to zero after borrowing it from +! positive values within the mass-flux transport layers +! + do i = 1,im + tsumn(i) = 0. + tsump(i) = 0. + rtnp(i) = 1. + enddo + do k = 1,km + do i = 1,im + if (cnvflg(i)) then + if(k > kb(i) .and. k <= ktcon(i)) then + if(q1(i,k) < 0.) tsumn(i) = tsumn(i) + q1(i,k) + if(q1(i,k) > 0.) tsump(i) = tsump(i) + q1(i,k) endif endif enddo enddo + do i = 1,im + if(cnvflg(i)) then + if(tsump(i) > 0. .and. tsumn(i) < 0.) then + if(tsump(i) > abs(tsumn(i))) then + rtnp(i) = tsumn(i) / tsump(i) + else + rtnp(i) = tsump(i) / tsumn(i) + endif + endif + endif + enddo + do k = 1,km + do i = 1,im + if (cnvflg(i)) then + if(k > kb(i) .and. k <= ktcon(i)) then + if(rtnp(i) < 0.) then + if(tsump(i) > abs(tsumn(i))) then + if(q1(i,k) < 0.) q1(i,k)= 0. + if(q1(i,k) > 0.) q1(i,k)=(1.+rtnp(i))*q1(i,k) + else + if(q1(i,k) < 0.) q1(i,k)=(1.+rtnp(i))*q1(i,k) + if(q1(i,k) > 0.) q1(i,k)=0. + endif + endif + endif + endif + enddo + enddo +! if (.not.hwrf_samfshal) then +! do n = 1, ntr - kk = n+2 - do k = 1, km +! + do k = 1, km + do i = 1, im + if (cnvflg(i)) then + if(k > kb(i) .and. k <= ktcon(i)) then + ctr(i,k,n) = ctr(i,k,n)+dellae(i,k,n)*xmb(i)*dt2 + dp = 1000. * del(i,k) + delebar(i,n)=delebar(i,n)+dellae(i,k,n)*xmb(i)*dp/grav + endif + endif + enddo + enddo +! +! Negative TKE, ozone, and aerosols are set to zero after borrowing them +! from positive values within the mass-flux transport layers +! + do i = 1,im + tsumn(i) = 0. + tsump(i) = 0. + rtnp(i) = 1. + enddo + do k = 1,km + do i = 1,im + if (cnvflg(i)) then + if(k > kb(i) .and. k <= ktcon(i)) then + if(ctr(i,k,n) < 0.) tsumn(i) = tsumn(i) + ctr(i,k,n) + if(ctr(i,k,n) > 0.) tsump(i) = tsump(i) + ctr(i,k,n) + endif + endif + enddo + enddo + do i = 1,im + if(cnvflg(i)) then + if(tsump(i) > 0. .and. tsumn(i) < 0.) then + if(tsump(i) > abs(tsumn(i))) then + rtnp(i) = tsumn(i) / tsump(i) + else + rtnp(i) = tsump(i) / tsumn(i) + endif + endif + endif + enddo + do k = 1,km + do i = 1,im + if (cnvflg(i)) then + if(k > kb(i) .and. k <= ktcon(i)) then + if(rtnp(i) < 0.) then + if(tsump(i) > abs(tsumn(i))) then + if(ctr(i,k,n)<0.) ctr(i,k,n)=0. + if(ctr(i,k,n)>0.) ctr(i,k,n)=(1.+rtnp(i))*ctr(i,k,n) + else + if(ctr(i,k,n)<0.) ctr(i,k,n)=(1.+rtnp(i))*ctr(i,k,n) + if(ctr(i,k,n)>0.) ctr(i,k,n)=0. + endif + endif + endif + endif + enddo + enddo +! + kk = n+2 + do k = 1, km do i = 1, im - if (cnvflg(i) .and. k <= kmax(i)) then - if(k <= ktcon(i)) then - ctr(i,k,n) = ctr(i,k,n)+dellae(i,k,n)*xmb(i)*dt2 - delebar(i,n)=delebar(i,n)+dellae(i,k,n)*xmb(i)*dp/grav + if (cnvflg(i)) then + if(k > kb(i) .and. k <= ktcon(i)) then qtr(i,k,kk) = ctr(i,k,n) endif endif enddo + enddo +! enddo - enddo +! + if (do_aerosols) then +! + do n = 1, ntc +! +! convert wet deposition to total mass deposited over dt2 and dp + do k = 2, km1 + do i = 1, im + if (cnvflg(i)) then + if(k > kb(i) .and. k < ktcon(i)) then + dp = 1000. * del(i,k) + wet_dep(i,k,n) = chem_pw(i,k,n)*grav/dp + wet_dep(i,k,n) = wet_dep(i,k,n)*xmb(i)*dt2*dp + endif + endif + enddo + enddo +! + kk = n + itc - 1 + do k = 2, km1 + do i = 1, im + if (cnvflg(i)) then + if(k > kb(i) .and. k < ktcon(i)) then + dp = 1000. * del(i,k) + if (qtr(i,k,kk) < 0.) then +! borrow negative mass from wet deposition + tem = -qtr(i,k,kk)*dp + if(wet_dep(i,k,n) >= tem) then + wet_dep(i,k,n) = wet_dep(i,k,n) - tem + qtr(i,k,kk) = 0. + else + wet_dep(i,k,n) = 0. + qtr(i,k,kk) = qtr(i,k,kk)+wet_dep(i,k,n)/dp + endif + endif + endif + endif + enddo + enddo +! + enddo +! + endif +! endif ! !> - Recalculate saturation specific humidity using the updated temperature. @@ -1832,10 +2136,11 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & qcond(i) = evef * (q1(i,k) - qeso(i,k)) & / (1. + el2orc * qeso(i,k) / t1(i,k)**2) dp = 1000. * del(i,k) + factor = dp / grav if(rn(i) > 0. .and. qcond(i) < 0.) then qevap(i) = -qcond(i) * (1.-exp(-.32*sqrt(dt2*rn(i)))) qevap(i) = min(qevap(i), rn(i)*1000.*grav/dp) - delq2(i) = delqev(i) + .001 * qevap(i) * dp / grav + delq2(i) = delqev(i) + .001 * qevap(i) * factor endif if(rn(i) > 0. .and. qcond(i) < 0. .and. & delq2(i) > rntot(i)) then @@ -1843,7 +2148,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & flg(i) = .false. endif if(rn(i) > 0. .and. qevap(i) > 0.) then - tem = .001 * dp / grav + tem = .001 * factor tem1 = qevap(i) * tem if(tem1 > rn(i)) then qevap(i) = rn(i) / tem @@ -1855,10 +2160,10 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & t1(i,k) = t1(i,k) - elocp * qevap(i) deltv(i) = - elocp*qevap(i)/dt2 delq(i) = + qevap(i)/dt2 - delqev(i) = delqev(i) + .001*dp*qevap(i)/grav + delqev(i) = delqev(i) + tem * qevap(i) endif - delqbar(i) = delqbar(i) + delq(i)*dp/grav - deltbar(i) = deltbar(i) + deltv(i)*dp/grav + delqbar(i) = delqbar(i) + delq(i) * factor + deltbar(i) = deltbar(i) + deltv(i) * factor endif endif enddo @@ -1946,20 +2251,20 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & ! endif !> - Store aerosol concentrations if present - if (.not. hwrf_samfshal) then - if (do_aerosols) then - do n = 1, ntc - kk = n + itc - 1 - do k = 1, km - do i = 1, im - if(cnvflg(i) .and. rn(i) > 0.) then - if (k <= kmax(i)) qtr(i,k,kk) = qaero(i,k,n) - endif - enddo - enddo - enddo - endif - endif +! if (.not. hwrf_samfshal) then +! if (do_aerosols) then +! do n = 1, ntc +! kk = n + itc - 1 +! do k = 1, km +! do i = 1, im +! if(cnvflg(i) .and. rn(i) > 0.) then +! if (k <= kmax(i)) qtr(i,k,kk) = qaero(i,k,n) +! endif +! enddo +! enddo +! enddo +! endif +! endif ! ! hchuang code change ! diff --git a/physics/satmedmfvdifq.F b/physics/satmedmfvdifq.F index bb02ccba9..6e0c1bd80 100644 --- a/physics/satmedmfvdifq.F +++ b/physics/satmedmfvdifq.F @@ -67,7 +67,7 @@ end subroutine satmedmfvdifq_finalize !! (mfscuq.f). !! \section detail_satmedmfvidfq GFS satmedmfvdifq Detailed Algorithm !! @{ - subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntiw,ntke, & + subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw,ntiw,ntke, & & grav,rd,cp,rv,hvap,hfus,fv,eps,epsm1, & & dv,du,tdt,rtg,u1,v1,t1,q1,swh,hlw,xmu,garea,zvfun, & & psk,rbsoil,zorl,u10m,v10m,fm,fh, & @@ -76,7 +76,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntiw,ntke, & & dspheat,dusfc,dvsfc,dtsfc,dqsfc,hpbl,dkt,dku, & & kinver,xkzm_m,xkzm_h,xkzm_s,dspfac,bl_upfr,bl_dnfr, & & rlmx,elmx,sfc_rlm, & - & ntoz,ntqv,dtend,dtidx,index_of_temperature,index_of_x_wind, & + & ntqv,dtend,dtidx,index_of_temperature,index_of_x_wind, & & index_of_y_wind,index_of_process_pbl,gen_tend,ldiag3d, & & errmsg,errflg) ! @@ -86,7 +86,8 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntiw,ntke, & implicit none ! !---------------------------------------------------------------------- - integer, intent(in) :: im, km, ntrac, ntcw, ntiw, ntke, ntoz,ntqv + integer, intent(in) :: im, km, ntrac, ntcw, ntrw, ntiw, & + & ntke, ntqv integer, intent(in) :: sfc_rlm integer, intent(in) :: kinver(:) integer, intent(out) :: kpbl(:) @@ -134,7 +135,8 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntiw,ntke, & !*** !*** local variables !*** - integer i,is,k,kk,n,ndt,km1,kmpbl,kmscu,ntrac1,idtend + integer i,is,k,n,ndt,km1,kmpbl,kmscu,ntrac1,idtend + integer kps,kbx,kmx integer lcld(im),kcld(im),krad(im),mrad(im) integer kx1(im), kpblx(im) ! @@ -187,6 +189,15 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntiw,ntke, & real(kind=kind_phys) tcdo(im,km), qcdo(im,km,ntrac), & ucdo(im,km), vcdo(im,km), & buod(im,km), xmfd(im,km) +! +! variables for Total Variation Diminishing (TVD) flux-limiter scheme +! + real(kind=kind_phys) e_half(im,km-1), e_diff(im,0:km-1), + & q_half(im,km-1,ntrac-1), + & qh(im,km-1,ntrac-1), + & q_diff(im,0:km-1,ntrac-1) + real(kind=kind_phys) rrkp, phkp + real(kind=kind_phys) tsumn(im), tsump(im), rtnp(im) ! logical pblflg(im), sfcflg(im), flg(im) logical scuflg(im), pcnvflg(im) @@ -241,7 +252,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntiw,ntke, & parameter(qmin=1.e-8,qlmin=1.e-12,zfmin=1.e-8) parameter(aphi5=5.,aphi16=16.) parameter(elmfac=1.0,elefac=1.0,cql=100.) - parameter(dw2min=1.e-4,dkmax=1000.,xkgdx=3000.) + parameter(dw2min=1.e-4,dkmax=1000.,xkgdx=1000.) parameter(qlcr=3.5e-5,zstblmax=2500.) parameter(xkinv1=0.15,xkinv2=0.3) parameter(h1=0.33333333,hcrinv=250.) @@ -249,15 +260,15 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntiw,ntke, & parameter(ce0=0.4,cs0=0.2) parameter(rchck=1.5,ndt=20) - gravi=1.0/grav - g=grav - gocp=g/cp + gravi = 1.0 / grav + g = grav + gocp = g / cp ! cont=cp/g ! conq=hvap/g ! conw=1.0/g ! for del in pa !! parameter(cont=1000.*cp/g,conq=1000.*hvap/g,conw=1000./g) !kpa - elocp=hvap/cp - el2orc=hvap*hvap/(rv*cp) + elocp = hvap / cp + el2orc = hvap * hvap / (rv * cp) ! !************************************************************************ ! Initialize CCPP error handling variables @@ -450,11 +461,11 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntiw,ntke, & cfly(i,k) = 0. clwt = 1.0e-6 * (plyr(i,k)*0.001) if (qlx(i,k) > clwt) then - onemrh= max(1.e-10, 1.0-rhly(i,k)) - tem1 = min(max((onemrh*qstl(i,k))**0.49,0.0001),1.0) - tem1 = cql / tem1 - value = max(min( tem1*qlx(i,k), 50.0), 0.0) - tem2 = sqrt(sqrt(rhly(i,k))) + onemrh = max(1.e-10, 1.0-rhly(i,k)) + tem1 = min(max((onemrh*qstl(i,k))**0.49,0.0001),1.0) + tem1 = cql / tem1 + value = max(min( tem1*qlx(i,k), 50.0), 0.0) + tem2 = sqrt(sqrt(rhly(i,k))) cfly(i,k) = min(max(tem2*(1.0-exp(-value)), 0.0), 1.0) endif enddo @@ -815,14 +826,14 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntiw,ntke, & endif enddo enddo - do kk = 1, ntrac1 + do n = 1, ntrac1 do k = 1, km do i = 1, im if(pcnvflg(i)) then - qcko(i,k,kk) = q1(i,k,kk) + qcko(i,k,n) = q1(i,k,n) endif if(scuflg(i)) then - qcdo(i,k,kk) = q1(i,k,kk) + qcdo(i,k,n) = q1(i,k,n) endif enddo enddo @@ -1302,6 +1313,126 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntiw,ntke, & enddo enddo ! +!-------------------------------------------------------- +! compute variables for TVD flux-limiter scheme +! on environmental subsidence and uplifting +! + kps = max(kmpbl, kmscu) +! +! for moisture and tracers including hydrometeors +! + do n=1,ntrac1 + do k=1,kps + do i=1,im + qh(i,k,n) = 0.5 * (q1(i,k,n)+q1(i,k+1,n)) + enddo + enddo + do k=1,kps + do i=1,im + q_diff(i,k,n) = q1(i,k,n) - q1(i,k+1,n) + enddo + enddo + do i=1,im + if(q1(i,1,n) >= 0.) then + q_diff(i,0,n) = max(0.,2.*q1(i,1,n)-q1(i,2,n))- + & q1(i,1,n) + else + q_diff(i,0,n) = min(0.,2.*q1(i,1,n)-q1(i,2,n))- + & q1(i,1,n) + endif + enddo + enddo +! + do n = 1, ntrac1 +! + do k = 1, kps + do i = 1, im + kmx = max(kpbl(i), krad(i)) + q_half(i,k,n) = qh(i,k,n) + if((pcnvflg(i) .or. scuflg(i)) .and. k < kmx) then + tem = 0. + if(pcnvflg(i) .and. k < kpbl(i)) then + tem = xmf(i,k) + endif + if(scuflg(i) .and. + & (k >= mrad(i) .and. k < krad(i))) then + tem = tem - xmfd(i,k) + endif + if(tem > 0.) then + rrkp = 0. + if(abs(q_diff(i,k,n)) > 1.e-22) + & rrkp = q_diff(i,k+1,n) / q_diff(i,k,n) + phkp = (rrkp+abs(rrkp)) / (1.+abs(rrkp)) + q_half(i,k,n) = q1(i,k+1,n) + + & phkp*(qh(i,k,n)-q1(i,k+1,n)) + elseif (tem < 0.) then + rrkp = 0. + if(abs(q_diff(i,k,n)) > 1.e-22) + & rrkp = q_diff(i,k-1,n) / q_diff(i,k,n) + phkp = (rrkp+abs(rrkp)) / (1.+abs(rrkp)) + q_half(i,k,n) = q1(i,k,n) + + & phkp*(qh(i,k,n)-q1(i,k,n)) + endif + endif + enddo + enddo +! + enddo +! +! for tke +! + do k=1,kps + do i=1,im + tkeh(i,k) = 0.5 * (tke(i,k)+tke(i,k+1)) + enddo + enddo + do k=1,kps + do i=1,im + e_diff(i,k) = tke(i,k) - tke(i,k+1) + enddo + enddo + do i=1,im + if(tke(i,1) >= 0.) then + e_diff(i,0) = max(0.,2.*tke(i,1)-tke(i,2))- + & tke(i,1) + else + e_diff(i,0) = min(0.,2.*tke(i,1)-tke(i,2))- + & tke(i,1) + endif + enddo +! + do k = 1, kps + do i = 1, im + kmx = max(kpbl(i), krad(i)) + e_half(i,k) = tkeh(i,k) + if((pcnvflg(i) .or. scuflg(i)) .and. (k < kmx)) then + tem = 0. + if(pcnvflg(i) .and. k < kpbl(i)) then + tem = xmf(i,k) + endif + if(scuflg(i) .and. + & (k >= mrad(i) .and. k < krad(i))) then + tem = tem - xmfd(i,k) + endif + if(tem > 0.) then + rrkp = 0. + if(abs(e_diff(i,k)) > 1.e-22) + & rrkp = e_diff(i,k+1) / e_diff(i,k) + phkp = (rrkp+abs(rrkp)) / (1.+abs(rrkp)) + e_half(i,k) = tke(i,k+1) + + & phkp*(tkeh(i,k)-tke(i,k+1)) + elseif (tem < 0.) then + rrkp = 0. + if(abs(e_diff(i,k)) > 1.e-22) + & rrkp = e_diff(i,k-1) / e_diff(i,k) + phkp = (rrkp+abs(rrkp)) / (1.+abs(rrkp)) + e_half(i,k) = tke(i,k) + + & phkp*(tkeh(i,k)-tke(i,k)) + endif + endif + enddo + enddo +! !---------------------------------------------------------------------- !> - Compute tridiagonal matrix elements for turbulent kinetic energy ! @@ -1328,10 +1459,9 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntiw,ntke, & ptem = 0.5 * tem2 * xmf(i,k) ptem1 = dtodsd * ptem ptem2 = dtodsu * ptem - tem = tke(i,k) + tke(i,k+1) ptem = qcko(i,k,ntke) + qcko(i,k+1,ntke) - f1(i,k) = f1(i,k)-(ptem-tem)*ptem1 - f1(i,k+1) = tke(i,k+1)+(ptem-tem)*ptem2 + f1(i,k) = f1(i,k) - ptem * ptem1 + f1(i,k+1) = tke(i,k+1) + ptem * ptem2 else f1(i,k+1) = tke(i,k+1) endif @@ -1341,12 +1471,20 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntiw,ntke, & ptem = 0.5 * tem2 * xmfd(i,k) ptem1 = dtodsd * ptem ptem2 = dtodsu * ptem - tem = tke(i,k) + tke(i,k+1) ptem = qcdo(i,k,ntke) + qcdo(i,k+1,ntke) - f1(i,k) = f1(i,k) + (ptem - tem) * ptem1 - f1(i,k+1) = f1(i,k+1) - (ptem - tem) * ptem2 + f1(i,k) = f1(i,k) + ptem * ptem1 + f1(i,k+1) = f1(i,k+1) - ptem * ptem2 endif endif +! + kmx = max(kpbl(i), krad(i)) + if((pcnvflg(i) .or. scuflg(i)) .and. (k < kmx)) then + ptem = tem2 * (xmf(i,k) - xmfd(i,k)) + ptem1 = dtodsd * ptem + ptem2 = dtodsu * ptem + f1(i,k) = f1(i,k) + e_half(i,k) * ptem1 + f1(i,k+1) = f1(i,k+1) - e_half(i,k) * ptem2 + endif ! enddo enddo @@ -1354,6 +1492,109 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntiw,ntke, & !> - Call tridit() to solve tridiagonal problem for TKE c call tridit(im,km,1,al,ad,au,f1,au,f1) +! +! Negative TKE is set to zero after borrowing it from positive +! values within the mass-flux transport layers +! + do i = 1,im + tsumn(i) = 0. + tsump(i) = 0. + rtnp(i) = 1. + enddo + do k = 1,kps + do i = 1,im + if(pcnvflg(i) .and. scuflg(i)) then + kbx = 1 + kmx = max(kpbl(i), krad(i)) + elseif(pcnvflg(i) .and. .not. scuflg(i)) then + kbx = 1 + kmx = kpbl(i) + elseif(.not. pcnvflg(i) .and. scuflg(i)) then + kbx = mrad(i) + kmx = krad(i) + endif + if((pcnvflg(i) .or. scuflg(i)) .and. + & (k >= kbx .and. k <= kmx)) then + if(f1(i,k) < 0.) tsumn(i) = tsumn(i) + f1(i,k) + if(f1(i,k) > 0.) tsump(i) = tsump(i) + f1(i,k) + endif + enddo + enddo + do i = 1,im + if(pcnvflg(i) .or. scuflg(i)) then + if(tsump(i) > 0. .and. tsumn(i) < 0.) then + if(tsump(i) > abs(tsumn(i))) then + rtnp(i) = tsumn(i) / tsump(i) + else + rtnp(i) = tsump(i) / tsumn(i) + endif + endif + endif + enddo + do k = 1,kps + do i = 1,im + if(pcnvflg(i) .and. scuflg(i)) then + kbx = 1 + kmx = max(kpbl(i), krad(i)) + elseif(pcnvflg(i) .and. .not. scuflg(i)) then + kbx = 1 + kmx = kpbl(i) + elseif(.not. pcnvflg(i) .and. scuflg(i)) then + kbx = mrad(i) + kmx = krad(i) + endif + if((pcnvflg(i) .or. scuflg(i)) .and. + & (k >= kbx .and. k <= kmx)) then + if(rtnp(i) < 0.) then + if(tsump(i) > abs(tsumn(i))) then + if(f1(i,k) < 0.) f1(i,k) = 0. + if(f1(i,k) > 0.) f1(i,k) = (1.+rtnp(i))*f1(i,k) + else + if(f1(i,k) < 0.) f1(i,k) = (1.+rtnp(i))*f1(i,k) + if(f1(i,k) > 0.) f1(i,k) = 0. + endif + endif + endif + enddo + enddo +! +! To remove negative TKEs which were leaked out of the mass-flux transport layers +! by eddy diffusion or potential negative TKEs from the diffusion scheme, +! positive TKEs are borrowed again now from the entire layers +! + do i = 1,im + tsumn(i) = 0. + tsump(i) = 0. + rtnp(i) = 1. + enddo + do k = 1,km + do i = 1,im + if(f1(i,k) < 0.) tsumn(i) = tsumn(i) + f1(i,k) + if(f1(i,k) > 0.) tsump(i) = tsump(i) + f1(i,k) + enddo + enddo + do i = 1,im + if(tsump(i) > 0. .and. tsumn(i) < 0.) then + if(tsump(i) > abs(tsumn(i))) then + rtnp(i) = tsumn(i) / tsump(i) + else + rtnp(i) = tsump(i) / tsumn(i) + endif + endif + enddo + do k = 1,km + do i = 1,im + if(rtnp(i) < 0.) then + if(tsump(i) > abs(tsumn(i))) then + if(f1(i,k) < 0.) f1(i,k) = 0. + if(f1(i,k) > 0.) f1(i,k) = (1.+rtnp(i))*f1(i,k) + else + if(f1(i,k) < 0.) f1(i,k) = (1.+rtnp(i))*f1(i,k) + if(f1(i,k) > 0.) f1(i,k) = 0. + endif + endif + enddo + enddo c !> - Recover the tendency of tke c @@ -1380,10 +1621,10 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntiw,ntke, & f2(i,1) = q1(i,1,1) + dtdz1(i) * evap(i) enddo if(ntrac1 >= 2) then - do kk = 2, ntrac1 - is = (kk-1) * km + do n = 2, ntrac1 + is = (n-1) * km do i = 1, im - f2(i,1+is) = q1(i,1,kk) + f2(i,1+is) = q1(i,1,n) enddo enddo endif @@ -1411,10 +1652,9 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntiw,ntke, & ptem = tcko(i,k) + tcko(i,k+1) f1(i,k) = f1(i,k)+dtodsd*dsdzt-(ptem-tem)*ptem1 f1(i,k+1) = t1(i,k+1)-dtodsu*dsdzt+(ptem-tem)*ptem2 - tem = q1(i,k,1) + q1(i,k+1,1) ptem = qcko(i,k,1) + qcko(i,k+1,1) - f2(i,k) = f2(i,k) - (ptem - tem) * ptem1 - f2(i,k+1) = q1(i,k+1,1) + (ptem - tem) * ptem2 + f2(i,k) = f2(i,k) - ptem * ptem1 + f2(i,k+1) = q1(i,k+1,1) + ptem * ptem2 else f1(i,k) = f1(i,k)+dtodsd*dsdzt f1(i,k+1) = t1(i,k+1)-dtodsu*dsdzt @@ -1430,51 +1670,64 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntiw,ntke, & tem = t1(i,k) + t1(i,k+1) f1(i,k) = f1(i,k) + (ptem - tem) * ptem1 f1(i,k+1) = f1(i,k+1) - (ptem - tem) * ptem2 - tem = q1(i,k,1) + q1(i,k+1,1) ptem = qcdo(i,k,1) + qcdo(i,k+1,1) - f2(i,k) = f2(i,k) + (ptem - tem) * ptem1 - f2(i,k+1) = f2(i,k+1) - (ptem - tem) * ptem2 + f2(i,k) = f2(i,k) + ptem * ptem1 + f2(i,k+1) = f2(i,k+1) - ptem * ptem2 endif endif +! + kmx = max(kpbl(i), krad(i)) + if((pcnvflg(i) .or. scuflg(i)) .and. (k < kmx)) then + ptem = tem2 * (xmf(i,k) - xmfd(i,k)) + ptem1 = dtodsd * ptem + ptem2 = dtodsu * ptem + f2(i,k) = f2(i,k) + q_half(i,k,1) * ptem1 + f2(i,k+1) = f2(i,k+1) - q_half(i,k,1) * ptem2 + endif +! enddo enddo ! if(ntrac1 >= 2) then - do kk = 2, ntrac1 - is = (kk-1) * km + do n = 2, ntrac1 + is = (n-1) * km do k = 1, km1 do i = 1, im + dtodsd = dt2/del(i,k) + dtodsu = dt2/del(i,k+1) + dsig = prsl(i,k)-prsl(i,k+1) + tem2 = dsig * rdzt(i,k) +! if(pcnvflg(i) .and. k < kpbl(i)) then - dtodsd = dt2/del(i,k) - dtodsu = dt2/del(i,k+1) - dsig = prsl(i,k)-prsl(i,k+1) - tem = dsig * rdzt(i,k) - ptem = 0.5 * tem * xmf(i,k) + ptem = 0.5 * tem2 * xmf(i,k) ptem1 = dtodsd * ptem ptem2 = dtodsu * ptem - tem1 = qcko(i,k,kk) + qcko(i,k+1,kk) - tem2 = q1(i,k,kk) + q1(i,k+1,kk) - f2(i,k+is) = f2(i,k+is) - (tem1 - tem2) * ptem1 - f2(i,k+1+is)= q1(i,k+1,kk) + (tem1 - tem2) * ptem2 + ptem = qcko(i,k,n) + qcko(i,k+1,n) + f2(i,k+is) = f2(i,k+is) - ptem * ptem1 + f2(i,k+1+is)= q1(i,k+1,n) + ptem * ptem2 else - f2(i,k+1+is) = q1(i,k+1,kk) + f2(i,k+1+is) = q1(i,k+1,n) endif ! if(scuflg(i)) then if(k >= mrad(i) .and. k < krad(i)) then - dtodsd = dt2/del(i,k) - dtodsu = dt2/del(i,k+1) - dsig = prsl(i,k)-prsl(i,k+1) - tem = dsig * rdzt(i,k) - ptem = 0.5 * tem * xmfd(i,k) + ptem = 0.5 * tem2 * xmfd(i,k) ptem1 = dtodsd * ptem ptem2 = dtodsu * ptem - tem1 = qcdo(i,k,kk) + qcdo(i,k+1,kk) - tem2 = q1(i,k,kk) + q1(i,k+1,kk) - f2(i,k+is) = f2(i,k+is) + (tem1 - tem2) * ptem1 - f2(i,k+1+is)= f2(i,k+1+is) - (tem1 - tem2) * ptem2 + ptem = qcdo(i,k,n) + qcdo(i,k+1,n) + f2(i,k+is) = f2(i,k+is) + ptem * ptem1 + f2(i,k+1+is)= f2(i,k+1+is) - ptem * ptem2 endif endif +! + kmx = max(kpbl(i), krad(i)) + if((pcnvflg(i) .or. scuflg(i)) .and. (k < kmx)) then + ptem = tem2 * (xmf(i,k) - xmfd(i,k)) + ptem1 = dtodsd * ptem + ptem2 = dtodsu * ptem + f2(i,k+is) = f2(i,k+is) + q_half(i,k,n) * ptem1 + f2(i,k+1+is) = f2(i,k+1+is) - q_half(i,k,n) * ptem2 + endif ! enddo enddo @@ -1484,6 +1737,296 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntiw,ntke, & !> - Call tridin() to solve tridiagonal problem for heat and moisture c call tridin(im,km,ntrac1,al,ad,au,f1,f2,au,f1,f2) +! +! Negative moisture is set to zero after borrowing it from +! positive values within the mass-flux transport layers +! + do i = 1,im + tsumn(i) = 0. + tsump(i) = 0. + rtnp(i) = 1. + enddo + do k = 1,kps + do i = 1,im + if(pcnvflg(i) .and. scuflg(i)) then + kbx = 1 + kmx = max(kpbl(i), krad(i)) + elseif(pcnvflg(i) .and. .not. scuflg(i)) then + kbx = 1 + kmx = kpbl(i) + elseif(.not. pcnvflg(i) .and. scuflg(i)) then + kbx = mrad(i) + kmx = krad(i) + endif + if((pcnvflg(i) .or. scuflg(i)) .and. + & (k >= kbx .and. k <= kmx)) then + if(f2(i,k) < 0.) tsumn(i) = tsumn(i) + f2(i,k) + if(f2(i,k) > 0.) tsump(i) = tsump(i) + f2(i,k) + endif + enddo + enddo + do i = 1,im + if(pcnvflg(i) .or. scuflg(i)) then + if(tsump(i) > 0. .and. tsumn(i) < 0.) then + if(tsump(i) > abs(tsumn(i))) then + rtnp(i) = tsumn(i) / tsump(i) + else + rtnp(i) = tsump(i) / tsumn(i) + endif + endif + endif + enddo + do k = 1,kps + do i = 1,im + if(pcnvflg(i) .and. scuflg(i)) then + kbx = 1 + kmx = max(kpbl(i), krad(i)) + elseif(pcnvflg(i) .and. .not. scuflg(i)) then + kbx = 1 + kmx = kpbl(i) + elseif(.not. pcnvflg(i) .and. scuflg(i)) then + kbx = mrad(i) + kmx = krad(i) + endif + if((pcnvflg(i) .or. scuflg(i)) .and. + & (k >= kbx .and. k <= kmx)) then + if(rtnp(i) < 0.) then + if(tsump(i) > abs(tsumn(i))) then + if(f2(i,k) < 0.) f2(i,k) = 0. + if(f2(i,k) > 0.) f2(i,k) = (1.+rtnp(i))*f2(i,k) + else + if(f2(i,k) < 0.) f2(i,k) = (1.+rtnp(i))*f2(i,k) + if(f2(i,k) > 0.) f2(i,k) = 0. + endif + endif + endif + enddo + enddo +! +! To remove negative moistures which were leaked out of the mass-flux transport layers +! by eddy diffusion or potential negative moistures from the diffusion scheme +! especially due to downward surface latent heat flux during nighttime, +! positive moistures are borrowed again now from the entire layers +! + do i = 1,im + tsumn(i) = 0. + tsump(i) = 0. + rtnp(i) = 1. + enddo + do k = 1,km + do i = 1,im + if(f2(i,k) < 0.) tsumn(i) = tsumn(i) + f2(i,k) + if(f2(i,k) > 0.) tsump(i) = tsump(i) + f2(i,k) + enddo + enddo + do i = 1,im + if(tsump(i) > 0. .and. tsumn(i) < 0.) then + if(tsump(i) > abs(tsumn(i))) then + rtnp(i) = tsumn(i) / tsump(i) + else + rtnp(i) = tsump(i) / tsumn(i) + endif + endif + enddo + do k = 1,km + do i = 1,im + if(rtnp(i) < 0.) then + if(tsump(i) > abs(tsumn(i))) then + if(f2(i,k) < 0.) f2(i,k) = 0. + if(f2(i,k) > 0.) f2(i,k) = (1.+rtnp(i))*f2(i,k) + else + if(f2(i,k) < 0.) f2(i,k) = (1.+rtnp(i))*f2(i,k) + if(f2(i,k) > 0.) f2(i,k) = 0. + endif + endif + enddo + enddo +! +! Negative hydrometeors & tracers are set to zero after +! borrowing them from positive values within the mass-flux +! transport layers +! +! For the negative liquid water, first borrow water from vapor +! and then borrow it from the other layers if there is still +! negative water +! + if(ntrac1 >= 2) then + is = (ntcw-1) * km + do k = 1,kps + do i = 1,im + if(pcnvflg(i) .and. scuflg(i)) then + kbx = 1 + kmx = max(kpbl(i), krad(i)) + elseif(pcnvflg(i) .and. .not. scuflg(i)) then + kbx = 1 + kmx = kpbl(i) + elseif(.not. pcnvflg(i) .and. scuflg(i)) then + kbx = mrad(i) + kmx = krad(i) + endif + if((pcnvflg(i) .or. scuflg(i)) .and. + & (k >= kbx .and. k <= kmx)) then + if(f2(i,k+is) < 0.) then + tem = f2(i,k) + f2(i,k+is) + if(tem >= 0.0) then + f2(i,k) = tem + f1(i,k) = f1(i,k) - elocp * f2(i,k+is) + f2(i,k+is) = 0. + elseif (f2(i,k) > 0.0) then + f2(i,k+is) = tem + f1(i,k) = f1(i,k) + elocp * f2(i,k) + f2(i,k) = 0. + endif + endif + endif + enddo + enddo + endif +! +! For the negative rain water, first borrow water from vapor +! and then borrow it from the other layers if there is still +! negative water +! + if(ntrac1 >= 2 .and. ntrw > 0) then + is = (ntrw-1) * km + do k = 1,kps + do i = 1,im + if(pcnvflg(i) .and. scuflg(i)) then + kbx = 1 + kmx = max(kpbl(i), krad(i)) + elseif(pcnvflg(i) .and. .not. scuflg(i)) then + kbx = 1 + kmx = kpbl(i) + elseif(.not. pcnvflg(i) .and. scuflg(i)) then + kbx = mrad(i) + kmx = krad(i) + endif + if((pcnvflg(i) .or. scuflg(i)) .and. + & (k >= kbx .and. k <= kmx)) then + if(f2(i,k+is) < 0.) then + tem = f2(i,k) + f2(i,k+is) + if(tem >= 0.0) then + f2(i,k) = tem + f1(i,k) = f1(i,k) - elocp * f2(i,k+is) + f2(i,k+is) = 0. + elseif (f2(i,k) > 0.0) then + f2(i,k+is) = tem + f1(i,k) = f1(i,k) + elocp * f2(i,k) + f2(i,k) = 0. + endif + endif + endif + enddo + enddo + endif +! + if(ntrac1 >= 2) then + do n = 2, ntrac1 + is = (n-1) * km +! + do i = 1,im + tsumn(i) = 0. + tsump(i) = 0. + rtnp(i) = 1. + enddo + do k = 1,kps + do i = 1,im + if(pcnvflg(i) .and. scuflg(i)) then + kbx = 1 + kmx = max(kpbl(i), krad(i)) + elseif(pcnvflg(i) .and. .not. scuflg(i)) then + kbx = 1 + kmx = kpbl(i) + elseif(.not. pcnvflg(i) .and. scuflg(i)) then + kbx = mrad(i) + kmx = krad(i) + endif + if((pcnvflg(i) .or. scuflg(i)) .and. + & (k >= kbx .and. k <= kmx)) then + if(f2(i,k+is) < 0.) tsumn(i) = tsumn(i) + f2(i,k+is) + if(f2(i,k+is) > 0.) tsump(i) = tsump(i) + f2(i,k+is) + endif + enddo + enddo + do i = 1,im + if(pcnvflg(i) .or. scuflg(i)) then + if(tsump(i) > 0. .and. tsumn(i) < 0.) then + if(tsump(i) > abs(tsumn(i))) then + rtnp(i) = tsumn(i) / tsump(i) + else + rtnp(i) = tsump(i) / tsumn(i) + endif + endif + endif + enddo + do k = 1,kps + do i = 1,im + if(pcnvflg(i) .and. scuflg(i)) then + kbx = 1 + kmx = max(kpbl(i), krad(i)) + elseif(pcnvflg(i) .and. .not. scuflg(i)) then + kbx = 1 + kmx = kpbl(i) + elseif(.not. pcnvflg(i) .and. scuflg(i)) then + kbx = mrad(i) + kmx = krad(i) + endif + if((pcnvflg(i) .or. scuflg(i)) .and. + & (k >= kbx .and. k <= kmx)) then + if(rtnp(i) < 0.) then + if(tsump(i) > abs(tsumn(i))) then + if(f2(i,k+is)<0.) f2(i,k+is)=0. + if(f2(i,k+is)>0.) f2(i,k+is)=(1.+rtnp(i))*f2(i,k+is) + else + if(f2(i,k+is)<0.) f2(i,k+is)=(1.+rtnp(i))*f2(i,k+is) + if(f2(i,k+is)>0.) f2(i,k+is)=0. + endif + endif + endif + enddo + enddo +! +! To remove negative hydrometeors & tracers which were leaked out of the mass-flux transport layers +! by eddy diffusion or potential negative hydrometeors & tracers from the diffusion scheme +! especially due to downward surface fluxes during nighttime, +! positive hydrometeors & tracers are borrowed again now from the entire layers +! + do i = 1,im + tsumn(i) = 0. + tsump(i) = 0. + rtnp(i) = 1. + enddo + do k = 1,km + do i = 1,im + if(f2(i,k+is) < 0.) tsumn(i) = tsumn(i) + f2(i,k+is) + if(f2(i,k+is) > 0.) tsump(i) = tsump(i) + f2(i,k+is) + enddo + enddo + do i = 1,im + if(tsump(i) > 0. .and. tsumn(i) < 0.) then + if(tsump(i) > abs(tsumn(i))) then + rtnp(i) = tsumn(i) / tsump(i) + else + rtnp(i) = tsump(i) / tsumn(i) + endif + endif + enddo + do k = 1,km + do i = 1,im + if(rtnp(i) < 0.) then + if(tsump(i) > abs(tsumn(i))) then + if(f2(i,k+is)<0.) f2(i,k+is)=0. + if(f2(i,k+is)>0.) f2(i,k+is)=(1.+rtnp(i))*f2(i,k+is) + else + if(f2(i,k+is)<0.) f2(i,k+is)=(1.+rtnp(i))*f2(i,k+is) + if(f2(i,k+is)>0.) f2(i,k+is)=0. + endif + endif + enddo + enddo +! + enddo + endif c !> - Recover the tendencies of heat and moisture c @@ -1503,7 +2046,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntiw,ntke, & enddo ! if(ldiag3d .and. .not. gen_tend) then - idtend=dtidx(index_of_temperature,index_of_process_pbl) + idtend = dtidx(index_of_temperature,index_of_process_pbl) if(idtend>=1) then do k = 1,km do i = 1,im @@ -1513,7 +2056,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntiw,ntke, & enddo endif ! Send tendencies just for QV; other tracers are below. - idtend=dtidx(100+ntqv,index_of_process_pbl) + idtend = dtidx(100+ntqv,index_of_process_pbl) if(idtend>=1) then do k = 1,km do i = 1,im @@ -1525,25 +2068,25 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntiw,ntke, & endif ! if(ntrac1 >= 2) then - do kk = 2, ntrac1 - is = (kk-1) * km + do n = 2, ntrac1 + is = (n-1) * km do k = 1, km do i = 1, im - qtend = (f2(i,k+is)-q1(i,k,kk))*rdt - rtg(i,k,kk) = rtg(i,k,kk)+qtend + qtend = (f2(i,k+is)-q1(i,k,n))*rdt + rtg(i,k,n) = rtg(i,k,n)+qtend enddo enddo enddo if(ldiag3d .and. .not. gen_tend) then ! Send tendencies for all tracers that were selected. - do kk = 2, ntrac1 - is = (kk-1) * km - idtend = dtidx(kk+100,index_of_process_pbl) + do n = 2, ntrac1 + is = (n-1) * km + idtend = dtidx(n+100,index_of_process_pbl) if(idtend>=1) then - if(kk/=ntke) then + if(n/=ntke) then do k = 1, km do i = 1, im - qtend = (f2(i,k+is)-q1(i,k,kk))*rdt + qtend = (f2(i,k+is)-q1(i,k,n))*rdt dtend(i,k,idtend) = dtend(i,k,idtend)+qtend*delt enddo enddo @@ -1565,7 +2108,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntiw,ntke, & enddo enddo if(ldiag3d .and. .not. gen_tend) then - idtend=dtidx(index_of_temperature,index_of_process_pbl) + idtend = dtidx(index_of_temperature,index_of_process_pbl) if(idtend>=1) then do k = 1,km1 do i = 1,im @@ -1657,7 +2200,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntiw,ntke, & enddo ! if(ldiag3d .and. .not. gen_tend) then - idtend=dtidx(index_of_x_wind,index_of_process_pbl) + idtend = dtidx(index_of_x_wind,index_of_process_pbl) if(idtend>=1) then do k = 1,km do i = 1,im @@ -1667,7 +2210,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntiw,ntke, & enddo endif - idtend=dtidx(index_of_y_wind,index_of_process_pbl) + idtend = dtidx(index_of_y_wind,index_of_process_pbl) if(idtend>=1) then do k = 1,km do i = 1,im diff --git a/physics/satmedmfvdifq.meta b/physics/satmedmfvdifq.meta index 3bed902bb..153a81bd4 100644 --- a/physics/satmedmfvdifq.meta +++ b/physics/satmedmfvdifq.meta @@ -76,6 +76,13 @@ dimensions = () type = integer intent = in +[ntrw] + standard_name = index_for_rain_water_vertical_diffusion_tracer + long_name = tracer index for rain water in the vertically diffused tracer array + units = index + dimensions = () + type = integer + intent = in [ntiw] standard_name = index_for_ice_cloud_condensate_vertical_diffusion_tracer long_name = tracer index for ice water in the vertically diffused tracer array @@ -566,13 +573,6 @@ dimensions = () type = integer intent = in -[ntoz] - standard_name = index_of_ozone_mixing_ratio_in_tracer_concentration_array - long_name = tracer index for ozone mixing ratio - units = index - dimensions = () - type = integer - intent = in [ntqv] standard_name = index_of_specific_humidity_in_tracer_concentration_array long_name = tracer index for water vapor (specific humidity)