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.

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

Merge gsd/develop into master #216

Merged
merged 4 commits into from
Mar 1, 2019
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
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
145 changes: 98 additions & 47 deletions physics/cu_gf_deep.F90
Original file line number Diff line number Diff line change
Expand Up @@ -6,9 +6,9 @@ module cu_gf_deep
real(kind=kind_phys), parameter::r_v=461.
real(kind=kind_phys), parameter :: tcrit=258.
! tuning constant for cloudwater/ice detrainment
real(kind=kind_phys), parameter:: c1=.002 ! .0005
real(kind=kind_phys), parameter:: c1= 0.003 !.002 ! .0005
! parameter to turn on or off evaporation of rainwater as done in sas
integer, parameter :: irainevap=1
integer, parameter :: irainevap=0
! max allowed fractional coverage (frh_thresh)
real(kind=kind_phys), parameter::frh_thresh = .9
! rh threshold. if fractional coverage ~ frh_thres, do not use cupa any further
Expand Down Expand Up @@ -362,6 +362,8 @@ subroutine cu_gf_deep_run( &
el2orc=xlv*xlv/(r_v*cp)
evfact=.2
evfactl=.2
!evfact=.0 ! for 4F5f
!evfactl=.4


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Expand All @@ -370,11 +372,11 @@ subroutine cu_gf_deep_run( &
!
! ecmwf
pgcon=0.
lambau(:)=5.
if(imid.eq.1)lambau(:)=5.
lambau(:)=2.0
if(imid.eq.1)lambau(:)=2.0
! here random must be between -1 and 1
if(nranflag == 1)then
lambau(:)=4.5+rand_mom(:)
lambau(:)=1.5+rand_mom(:)
endif
! sas
! lambau=0.
Expand Down Expand Up @@ -445,7 +447,7 @@ subroutine cu_gf_deep_run( &
entr_rate(i)=7.e-5 - min(20.,float(csum(i))) * 3.e-6
if(xland1(i) == 0)entr_rate(i)=7.e-5
if(imid.eq.1)entr_rate(i)=3.e-4
if(imid.eq.1)c1d(i,:)=c1 ! comment to test warm bias 08/14/17
! if(imid.eq.1)c1d(i,:)=c1 ! comment to test warm bias 08/14/17
radius=.2/entr_rate(i)
frh=min(1.,3.14*radius*radius/dx(i)/dx(i))
if(frh > frh_thresh)then
Expand Down Expand Up @@ -756,10 +758,18 @@ subroutine cu_gf_deep_run( &
!
! calculate mass entrainment and detrainment
!
if(imid.eq.1)then
call get_lateral_massflux(itf,ktf, its,ite, kts,kte &
,ierr,ktop,zo_cup,zuo,cd,entr_rate_2d &
,up_massentro, up_massdetro ,up_massentr, up_massdetr &
,'mid',kbcon,k22,up_massentru,up_massdetru,lambau)
else
call get_lateral_massflux(itf,ktf, its,ite, kts,kte &
,ierr,ktop,zo_cup,zuo,cd,entr_rate_2d &
,up_massentro, up_massdetro ,up_massentr, up_massdetr &
,'deep',kbcon,k22,up_massentru,up_massdetru,lambau)
endif


!
! note: ktop here already includes overshooting, ktopdby is without
Expand Down Expand Up @@ -892,10 +902,27 @@ subroutine cu_gf_deep_run( &
do i=its,itf
if(ierr(i) /= 0) cycle
! do k=kbcon(i)+1,ktop(i)-1
do k=jmin(i)+1,ktop(i)-1
c1d(i,k)=c1
enddo
!c do k=jmin(i)+1,ktop(i)-1
!c c1d(i,k)=c1
!c enddo
!if(imid.eq.1)c1d(i,:)=0.
! do k=kts,ktop(i)
! if(po(i,k).gt.700.)then
! c1d(i,k)=0.
! elseif(po(i,k).gt.600.)then
! c1d(i,k)=0.001
! elseif(po(i,k).gt.500.)then
! c1d(i,k)=0.002
! elseif(po(i,k).gt.400.)then
! c1d(i,k)=0.003
! elseif(po(i,k).gt.300.)then
! c1d(i,k)=0.004
! elseif(po(i,k).gt.200.)then
! c1d(i,k)=0.005
! endif
! enddo
! if(imid.eq.1)c1d(i,:)=0.003

do k=ktop(i)+1,ktf
hco(i,k)=heso_cup(i,k)
dbyo(i,k)=0.
Expand Down Expand Up @@ -1161,31 +1188,31 @@ subroutine cu_gf_deep_run( &
endif
if(zdo(i,ki+1).gt.0.)cdd(i,ki)= dd_massdetro(i,ki)/(dzo*zdo(i,ki+1))
enddo
! cbeg=po_cup(i,kbcon(i)) !850.
! cend=min(po_cup(i,ktop(i)),400.)
! cbeg=800. !po_cup(i,kbcon(i)) !850.
! cend=min(po_cup(i,ktop(i)),200.)
! cmid=.5*(cbeg+cend) !600.
! const_b=c1/((cmid*cmid-cbeg*cbeg)*(cbeg-cend)/(cend*cend-cbeg*cbeg)+cmid-cbeg)
! const_a=const_b*(cbeg-cend)/(cend*cend-cbeg*cbeg)
! const_c=-const_a*cbeg*cbeg-const_b*cbeg
! do k=kbcon(i)+1,ktop(i)-1
! c1d(i,k)=const_a*po_cup(i,k)*po_cup(i,k)+const_b*po_cup(i,k)+const_c
! c1d(i,k)=max(0.,c1d(i,k))
! c1d(i,k)=c1
!! c1d(i,k)=c1
! enddo
! if(imid.eq.1)c1d(i,:)=0.
! do k=1,jmin(i)
! c1d(i,k)=0.
! enddo
! c1d(i,jmin(i)-2)=c1/40.
! if(imid.eq.1)c1d(i,jmin(i)-2)=c1/20.
! do k=jmin(i)-1,ktop(i)
! dz=zo_cup(i,ktop(i))-zo_cup(i,jmin(i))
! c1d(i,k)=c1d(i,k-1)+c1*(zo_cup(i,k+1)-zo_cup(i,k))/dz
! c1d(i,k)=max(0.,c1d(i,k))
! c1d(i,k)=min(.002,c1d(i,k))
! enddo


!! if(imid.eq.1)c1d(i,:)=0.
!! do k=1,jmin(i)
!! c1d(i,k)=0.
!! enddo
!! c1d(i,jmin(i)-2)=c1/40.
!! if(imid.eq.1)c1d(i,jmin(i)-2)=c1/20.
!! do k=jmin(i)-1,ktop(i)
!! dz=zo_cup(i,ktop(i))-zo_cup(i,jmin(i))
!! c1d(i,k)=c1d(i,k-1)+c1*(zo_cup(i,k+1)-zo_cup(i,k))/dz
!! c1d(i,k)=max(0.,c1d(i,k))
!! c1d(i,k)=min(.002,c1d(i,k))
!! enddo
!
!
! downdraft moist static energy + moisture budget
do k=2,jmin(i)+1
dd_massentru(i,k-1)=dd_massentro(i,k-1)+lambau(i)*dd_massdetro(i,k-1)
Expand Down Expand Up @@ -1621,16 +1648,16 @@ subroutine cu_gf_deep_run( &
!-- take out cloud liquid water for detrainment
detup=up_massdetro(i,k)
dz=zo_cup(i,k)-zo_cup(i,k-1)
! if(k.lt.ktop(i) .and. k.ge.jmin(i)) then
! if(k.lt.ktop(i) .and. c1d(i,k).gt.0) then
!! if(k.lt.ktop(i) .and. k.ge.jmin(i)) then
!! if(k.lt.ktop(i) .and. c1d(i,k).gt.0) then
if(k.lt.ktop(i)) then
dellaqc(i,k) = zuo(i,k)*c1d(i,k)*qrco(i,k)*dz/dp*g
else
dellaqc(i,k)= detup*0.5*(qrco(i,k+1)+qrco(i,k)) *g/dp
endif
! if(imid.eq.1) dellaqc(i,k)= detup*0.5*(qrco(i,k+1)+qrco(i,k)) *g/dp
if(k.eq.ktop(i))dellaqc(i,k)= detup*0.5*(qrco(i,k+1)+qrco(i,k)) *g/dp
!---
!! if(imid.eq.1) dellaqc(i,k)= detup*0.5*(qrco(i,k+1)+qrco(i,k)) *g/dp
! if(k.eq.ktop(i))dellaqc(i,k)= detup*0.5*(qrco(i,k+1)+qrco(i,k)) *g/dp
! !---
g_rain= 0.5*(pwo (i,k)+pwo (i,k+1))*g/dp
e_dn = -0.5*(pwdo(i,k)+pwdo(i,k+1))*g/dp*edto(i) ! pwdo < 0 and e_dn must > 0
!-- condensation source term = detrained + flux divergence of
Expand Down Expand Up @@ -1939,7 +1966,7 @@ subroutine cu_gf_deep_run( &
do k = ktop(i), 1, -1
rain = pwo(i,k) + edto(i) * pwdo(i,k)
rn(i) = rn(i) + rain * xmb(i) * .001 * dtime
if(po(i,k).gt.700.)then
!if(po(i,k).gt.700.)then
if(flg(i))then
q1=qo(i,k)+(outq(i,k))*dtime
t1=tn(i,k)+(outt(i,k))*dtime
Expand All @@ -1964,7 +1991,7 @@ subroutine cu_gf_deep_run( &
pre(i)=max(pre(i),0.)
delqev(i) = delqev(i) + .001*dp*qevap(i)/g
endif
endif ! 700mb
!endif ! 700mb
endif
enddo
! pre(i)=1000.*rn(i)/dtime
Expand Down Expand Up @@ -3626,7 +3653,7 @@ subroutine cup_up_moisture(name,ierr,z_cup,qc,qrc,pw,pwav, &
iprop,iall,i,k
integer :: start_level(its:ite)
real(kind=kind_phys) :: &
prop_ave,qrcb_h,bdsp,dp,rhoc,qrch,qaver, &
prop_ave,qrcb_h,bdsp,dp,rhoc,qrch,qaver,clwdet, &
c0,dz,berryc0,q1,berryc
real(kind=kind_phys) :: &
denom, c0t
Expand All @@ -3636,6 +3663,7 @@ subroutine cup_up_moisture(name,ierr,z_cup,qc,qrc,pw,pwav, &
prop_b(kts:kte)=0
iall=0
c0=.002
clwdet=100.
bdsp=bdispm
!
!--- no precip for small clouds
Expand Down Expand Up @@ -3813,7 +3841,14 @@ subroutine cup_up_moisture(name,ierr,z_cup,qc,qrc,pw,pwav, &
pw(i,k)=(qc(i,k)-qrch)*zu(i,k)
if(pw(i,k).lt.0.)pw(i,k)=0.
else
! create clw detrainment profile that depends on mass detrainment and
! in-cloud clw/ice
!
c1d(i,k)=clwdet*up_massdetr(i,k-1)*qrc(i,k-1)
qrc(i,k)=(qc(i,k)-qrch)/(1.+(c1d(i,k)+c0t)*dz)
if(qrc(i,k).lt.0.)then ! hli new test 02/12/19
qrc(i,k)=0.
endif
pw(i,k)=c0t*dz*qrc(i,k)*zu(i,k)
!-----srf-08aug2017-----begin
! pw(i,k)=(c1d(i,k)+c0)*dz*max(0.,qrc(i,k) -qrc_crit)! units kg[rain]/kg[air]
Expand Down Expand Up @@ -3919,7 +3954,7 @@ subroutine rates_up_pdf(rand_vmas,ipr,name,ktop,ierr,p_cup,entr_rate_2d,hkbo,heo
integer, dimension (its:ite) :: start_level
!
zustart=.1
dbythresh= 1. !.0.95 ! 0.85, 0.6
dbythresh= 0.8 !.0.95 ! 0.85, 0.6
if(name == 'shallow' .or. name == 'mid') dbythresh=1.
dby(:)=0.

Expand Down Expand Up @@ -3969,6 +4004,7 @@ subroutine rates_up_pdf(rand_vmas,ipr,name,ktop,ierr,p_cup,entr_rate_2d,hkbo,heo
kfinalzu=ktf-2
ktop(i)=kfinalzu
412 continue
kklev=min(kklev+3,ktop(i)-2)
!
! at least overshoot by one level
!
Expand Down Expand Up @@ -4017,9 +4053,10 @@ subroutine get_zu_zd_pdf_fim(kklev,p,rand_vmas,zubeg,ipr,xland,zuh2,draft,ierr,k
implicit none
! real(kind=kind_phys), parameter :: beta_deep=1.3,g_beta_deep=0.8974707
! real(kind=kind_phys), parameter :: beta_deep=1.2,g_beta_deep=0.8974707
real(kind=kind_phys), parameter :: beta_sh=2.5,g_beta_sh=1.329340
! real(kind=kind_phys), parameter :: beta_mid=1.3,g_beta_mid=0.8974707
real(kind=kind_phys), parameter :: beta_mid=2.2,g_beta_mid=0.8974707
! real(kind=kind_phys), parameter :: beta_sh=2.5,g_beta_sh=1.329340
real(kind=kind_phys), parameter :: beta_sh=2.2,g_beta_sh=0.8974707
real(kind=kind_phys), parameter :: beta_mid=1.3,g_beta_mid=0.8974707
! real(kind=kind_phys), parameter :: beta_mid=1.8,g_beta_mid=0.8974707
real(kind=kind_phys), parameter :: beta_dd=4.0,g_beta_dd=6.
integer, intent(in) ::ipr,xland,kb,kklev,kt,kts,kte,ktf,kpbli,csum,pmin_lev
real(kind=kind_phys), intent(in) ::max_mass,zubeg
Expand Down Expand Up @@ -4064,6 +4101,7 @@ subroutine get_zu_zd_pdf_fim(kklev,p,rand_vmas,zubeg,ipr,xland,zuh2,draft,ierr,k
lev_start=min(.9,.1+csum*.013)
kb_adj=max(kb,2)
tunning=max(p(kklev+1),.5*(p(kpbli)+p(kt)))
tunning=p(kklev)
! tunning=p(kklev+1) !p(kpbli+1) !p(kklev) !p(kt)+(p(kpbli)-p(kt))*lev_start
! tunning=.5*(p(kb_adj)+p(kt)) !p(kpbli+1) !p(kklev) !p(kt)+(p(kpbli)-p(kt))*lev_start
trash=-p(kt)+p(kb_adj)
Expand Down Expand Up @@ -4175,7 +4213,7 @@ subroutine get_zu_zd_pdf_fim(kklev,p,rand_vmas,zubeg,ipr,xland,zuh2,draft,ierr,k
kb_adj=max(kb,2)
tunning=.5*(p(kt)+p(kpbli)) !p(kt)+(p(kb_adj)-p(kt))*.9 !*.33
!new nov18
tunning=p(kpbli) !p(kpbli+1) !p(kklev) !p(kt)+(p(kpbli)-p(kt))*lev_start
! tunning=p(kpbli) !p(kpbli+1) !p(kklev) !p(kt)+(p(kpbli)-p(kt))*lev_start
!end new
tunning =min(0.95, (tunning-p(kb_adj))/(p(kt)-p(kb_adj))) !=.6
tunning =max(0.02, tunning)
Expand Down Expand Up @@ -4517,7 +4555,8 @@ subroutine get_lateral_massflux(itf,ktf, its,ite, kts,kte
character *(*), intent (in) :: draft
integer, intent(in):: itf,ktf, its,ite, kts,kte
integer, intent(in) , dimension(its:ite) :: ierr,ktop,kbcon,k22
real(kind=kind_phys), intent(in), optional , dimension(its:ite):: lambau
!real(kind=kind_phys), intent(in), optional , dimension(its:ite):: lambau
real(kind=kind_phys), intent(inout), optional , dimension(its:ite):: lambau
real(kind=kind_phys), intent(in) , dimension(its:ite,kts:kte) :: zo_cup,zuo
real(kind=kind_phys), intent(inout), dimension(its:ite,kts:kte) :: cd,entr_rate_2d
real(kind=kind_phys), intent( out), dimension(its:ite,kts:kte) :: up_massentro, up_massdetro &
Expand Down Expand Up @@ -4585,13 +4624,25 @@ subroutine get_lateral_massflux(itf,ktf, its,ite, kts,kte
up_massentr (i,k-1)=up_massentro(i,k-1)
up_massdetr (i,k-1)=up_massdetro(i,k-1)
enddo
if(present(up_massentru) .and. present(up_massdetru))then
turn=maxloc(zuo(i,:),1)
do k=2,turn
up_massentru(i,k-1)=up_massentro(i,k-1)+.1*lambau(i)*up_massentro(i,k-1)
up_massdetru(i,k-1)=up_massdetro(i,k-1)+.1*lambau(i)*up_massentro(i,k-1)
if(present(up_massentru) .and. present(up_massdetru) .and. draft == 'deep')then
!turn=maxloc(zuo(i,:),1)
!do k=2,turn
! up_massentru(i,k-1)=up_massentro(i,k-1)+.1*lambau(i)*up_massentro(i,k-1)
! up_massdetru(i,k-1)=up_massdetro(i,k-1)+.1*lambau(i)*up_massentro(i,k-1)
!enddo
!do k=turn+1,ktf-1
do k=2,ktf-1
up_massentru(i,k-1)=up_massentro(i,k-1)+lambau(i)*up_massdetro(i,k-1)
up_massdetru(i,k-1)=up_massdetro(i,k-1)+lambau(i)*up_massdetro(i,k-1)
enddo
else if(present(up_massentru) .and. present(up_massdetru) .and. draft == 'shallow')then
do k=2,ktf-1
up_massentru(i,k-1)=up_massentro(i,k-1)+lambau(i)*up_massdetro(i,k-1)
up_massdetru(i,k-1)=up_massdetro(i,k-1)+lambau(i)*up_massdetro(i,k-1)
enddo
do k=turn+1,ktf-1
else if(present(up_massentru) .and. present(up_massdetru) .and. draft == 'mid')then
lambau(i)=0.
do k=2,ktf-1
up_massentru(i,k-1)=up_massentro(i,k-1)+lambau(i)*up_massdetro(i,k-1)
up_massdetru(i,k-1)=up_massdetro(i,k-1)+lambau(i)*up_massdetro(i,k-1)
enddo
Expand Down Expand Up @@ -4791,7 +4842,7 @@ subroutine get_cloud_top(name,ktop,ierr,p_cup,entr_rate_2d,hkbo,heo,heso_cup,z_c
integer, dimension (its:ite) :: start_level
integer,parameter :: find_ktop_option = 1 !0=original, 1=new

dbythresh=1. !0.95 ! the range of this parameter is 0-1, higher => lower
dbythresh=0.8 !0.95 ! the range of this parameter is 0-1, higher => lower
! overshoting (cheque aa0 calculation)
! rainfall is too sensible this parameter
! for now, keep =1.
Expand Down
13 changes: 6 additions & 7 deletions physics/cu_gf_driver.F90
Original file line number Diff line number Diff line change
Expand Up @@ -374,8 +374,7 @@ subroutine cu_gf_driver_run(tottracer,ntrac,garea,im,ix,km,dt,cactiv, &
ierrs(:)=0
cuten(:)=0.
cutenm(:)=0.
cutens(:)=1.
if(ishallow_g3.eq.0)cutens(:)=0.
cutens(:)=0.
ierrc(:)=" "

kbcon(:)=0
Expand Down Expand Up @@ -531,21 +530,21 @@ subroutine cu_gf_driver_run(tottracer,ntrac,garea,im,ix,km,dt,cactiv, &
!> if ishallow_g3=1, call shallow: cup_gf_sh()
!
! print*,'hli bf shallow t2d',t2d
call cu_gf_sh_run ( &
call cu_gf_sh_run (us,vs, &
! input variables, must be supplied
zo,t2d,q2d,ter11,tshall,qshall,p2d,psur,dhdt,kpbli, &
rhoi,hfx,qfx,xlandi,ichoice_s,tcrit,dt, &
! input variables. ierr should be initialized to zero or larger than zero for
! turning off shallow convection for grid points
zus,xmbs,kbcons,ktops,k22s,ierrs,ierrcs, &
! output tendencies
outts,outqs,outqcs,cnvwt,prets,cupclws, &
outts,outqs,outqcs,outus,outvs,cnvwt,prets,cupclws, &
! dimesnional variables
itf,ktf,its,ite, kts,kte,ipr,tropics)


do i=its,itf
if(xmbs(i).le.0.)cutens(i)=0.
if(xmbs(i).gt.0.)cutens(i)=1.
enddo
call neg_check('shallow',ipn,dt,qcheck,outqs,outts,outus,outvs, &
outqcs,prets,its,ite,kts,kte,itf,ktf,ktops)
Expand Down Expand Up @@ -779,8 +778,8 @@ subroutine cu_gf_driver_run(tottracer,ntrac,garea,im,ix,km,dt,cactiv, &
t(i,k)=t(i,k)+dt*(cutens(i)*outts(i,k)+cutenm(i)*outtm(i,k)+outt(i,k)*cuten(i))
qv(i,k)=max(1.e-16,qv(i,k)+dt*(cutens(i)*outqs(i,k)+cutenm(i)*outqm(i,k)+outq(i,k)*cuten(i)))
gdc(i,k,7)=sqrt(us(i,k)**2 +vs(i,k)**2)
us(i,k)=us(i,k)+outu(i,k)*cuten(i)*dt +outum(i,k)*cutenm(i)*dt
vs(i,k)=vs(i,k)+outv(i,k)*cuten(i)*dt +outvm(i,k)*cutenm(i)*dt
us(i,k)=us(i,k)+outu(i,k)*cuten(i)*dt +outum(i,k)*cutenm(i)*dt +outus(i,k)*cutens(i)*dt
vs(i,k)=vs(i,k)+outv(i,k)*cuten(i)*dt +outvm(i,k)*cutenm(i)*dt +outvs(i,k)*cutens(i)*dt

!hj 10/11/2016: don't need gdc and gdc2 yet for gsm.
!hli 08/18/2017: couple gdc to radiation
Expand Down
Loading