Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Implement positive-definite TVD methods for tracer transport in EDMF-TKE PBL schemes and SAMF convective schemes #754

Merged
merged 10 commits into from
Nov 9, 2021
Merged

Implement positive-definite TVD methods for tracer transport in EDMF-TKE PBL schemes and SAMF convective schemes #754

Show file tree
Hide file tree
Changes from 2 commits
Commits
Show all changes
10 commits
Select commit Hold shift + click to select a range
ed3768b
Implement positive-definite Total Variation Diminishing (TVD)
rmontuoro Oct 5, 2021
c1f0706
Remove unused input argument for ozone tracer index.
rmontuoro Oct 6, 2021
2758810
Use proper index for rain water tracer in vertically-diffused tracer …
rmontuoro Oct 6, 2021
9a7131f
Adjust negative mixing ratios for rain water only if such tracer is a…
rmontuoro Oct 6, 2021
7d080b8
Improve test for rain water tracer.
rmontuoro Oct 6, 2021
975afe7
Merge branch 'main' into feature/tvd
rmontuoro Nov 5, 2021
0014964
Remove remaining optional arguments from metadata.
rmontuoro Nov 5, 2021
1455ddd
Adding spaces to improve code readability as suggested by S. Moorthi.
rmontuoro Nov 5, 2021
6bc90ec
Precompute a few quantities in SAMF deep & shallow convective schemes
rmontuoro Nov 8, 2021
3c00f02
Remove multiplication by -1.
rmontuoro Nov 8, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
140 changes: 69 additions & 71 deletions physics/samfdeepcnv.f
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1210,12 +1210,12 @@ subroutine samfdeepcnv_run (im,km,itc,ntc,cliq,cp,cvap, &
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)
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
Expand Down Expand Up @@ -1792,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
Expand Down Expand Up @@ -1992,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 = -1.*(zo(i,k+1)-zo(i,k))
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Not trying to be picky, but we really don't need this multiplication by "-1".

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

it could be just a "-" or you could reverse the two items inside "( )"

Copy link
Contributor Author

@rmontuoro rmontuoro Nov 8, 2021 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Agree! However, this should not add compute time (the compiler should recognize this).

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@SMoorthi-emc - I've implemented your suggested changes.

! 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
Expand Down Expand Up @@ -2050,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)
& * 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
Expand Down Expand Up @@ -2099,34 +2096,36 @@ 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
& - (aup*eta(i,k-1)-adw*edto(i)*etad(i,k-1))*dv3h
& - (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*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
Expand Down Expand Up @@ -2161,20 +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
dellaq(i,indx) = eta(i,indx-1) *
& qcko(i,indx-1) * 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
Expand Down Expand Up @@ -2638,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 = -1.*(zo(i,k+1)-zo(i,k))
! 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
Expand Down Expand Up @@ -2952,20 +2946,22 @@ 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
Expand Down Expand Up @@ -3187,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
Expand Down
61 changes: 30 additions & 31 deletions physics/samfshalcnv.f
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1014,12 +1014,12 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, &
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)
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
Expand Down Expand Up @@ -1569,25 +1569,27 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, &
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) +
& ( - 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
Expand Down Expand Up @@ -1618,20 +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
dellaq(i,indx) = eta(i,indx-1) *
& qcko(i,indx-1) * 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
Expand Down Expand Up @@ -1906,11 +1903,12 @@ 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
Expand Down Expand Up @@ -2138,18 +2136,19 @@ 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
qevap(i) = 1000.* grav * (rntot(i) - delqev(i)) / dp
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
Expand All @@ -2161,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
Expand Down
Loading