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

Move number concentration bug fix, reorganization of interstitial code #15

Merged
merged 7 commits into from
Feb 18, 2020
Merged

Move number concentration bug fix, reorganization of interstitial code #15

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
1c6cad5
Number Concentrated code moved to interstitial code
hannahcbarnes Jan 6, 2020
16f6a02
Merge branch 'NCmove' of https://github.com/hannahcbarnes/ccpp-physic…
climbfuji Jan 22, 2020
8d5fe8c
physics/module_mp_thompson.F90: bugfix, remove threaded computation/r…
climbfuji Jan 23, 2020
d2f38dd
Reorganize interstitial code around convection, bugfixes for Hannah's…
climbfuji Feb 13, 2020
f0e8369
Merge branch 'gsd/develop' of https://github.com/noaa-gsd/ccpp-physic…
climbfuji Feb 13, 2020
c796456
Merge branch 'bugfix_threading_thompson' of https://github.com/climbf…
climbfuji Feb 13, 2020
cebdfa4
Minor cleanup of physics/GFS_suite_interstitial.F90
climbfuji Feb 17, 2020
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
117 changes: 53 additions & 64 deletions physics/GFS_suite_interstitial.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -463,13 +463,13 @@ end subroutine GFS_suite_interstitial_3_finalize
subroutine GFS_suite_interstitial_3_run (im, levs, nn, cscnv, &
satmedmf, trans_trac, do_shoc, ltaerosol, ntrac, ntcw, &
ntiw, ntclamt, ntrw, ntsw, ntrnc, ntsnc, ntgl, ntgnc, &
xlat, gq0, imp_physics, imp_physics_mg, &
xlat, gt0, gq0, imp_physics, imp_physics_mg, &
imp_physics_zhao_carr, imp_physics_zhao_carr_pdf, &
imp_physics_gfdl, imp_physics_thompson, &
imp_physics_wsm6, imp_physics_fer_hires, prsi, &
prsl, prslk, rhcbot,rhcpbl, rhctop, rhcmax, islmsk, &
work1, work2, kpbl, kinver,clw, rhc, save_qc, save_qi, &
errmsg, errflg)
save_tcp, errmsg, errflg)

use machine, only: kind_phys

Expand All @@ -487,11 +487,13 @@ subroutine GFS_suite_interstitial_3_run (im, levs, nn, cscnv, &
real(kind=kind_phys), dimension(im, levs), intent(in) :: prsl, prslk
real(kind=kind_phys), dimension(im, levs+1), intent(in) :: prsi
real(kind=kind_phys), dimension(im), intent(in) :: xlat
real(kind=kind_phys), dimension(im, levs), intent(in) :: gt0
real(kind=kind_phys), dimension(im, levs, ntrac), intent(in) :: gq0

real(kind=kind_phys), dimension(im, levs), intent(inout) :: rhc, save_qc
! save_qi is not allocated for Zhao-Carr MP
real(kind=kind_phys), dimension(:, :), intent(inout) :: save_qi
real(kind=kind_phys), dimension(:, :), intent(inout) :: save_tcp ! ONLY ALLOCATE FOR THOMPSON! TODO
real(kind=kind_phys), dimension(im, levs, nn), intent(inout) :: clw

character(len=*), intent(out) :: errmsg
Expand All @@ -512,33 +514,6 @@ subroutine GFS_suite_interstitial_3_run (im, levs, nn, cscnv, &
errmsg = ''
errflg = 0

!GF* The following section (initializing convective variables) is already executed in GFS_typedefs%interstitial_phys_reset
! do k=1,levs
! do i=1,im
! clw(i,k,1) = 0.0
! clw(i,k,2) = -999.9
! enddo
! enddo
! if (Model%imfdeepcnv >= 0 .or. Model%imfshalcnv > 0 .or. &
! (Model%npdf3d == 3 .and. Model%num_p3d == 4) .or. &
! (Model%npdf3d == 0 .and. Model%ncnvcld3d == 1) ) then
! do k=1,levs
! do i=1,im
! cnvc(i,k) = 0.0
! cnvw(i,k) = 0.0
! enddo
! enddo
! endif
! if(imp_physics == 8) then
! if(Model%ltaerosol) then
! ice00 (:,:) = 0.0
! liq0 (:,:) = 0.0
! else
! ice00 (:,:) = 0.0
! endif
! endif
!*GF

if (cscnv .or. satmedmf .or. trans_trac ) then
tracers = 2
do n=2,ntrac
Expand Down Expand Up @@ -596,6 +571,8 @@ subroutine GFS_suite_interstitial_3_run (im, levs, nn, cscnv, &
enddo
enddo
endif
else
rhc(:,:) = 1.0
endif

if (imp_physics == imp_physics_zhao_carr .or. imp_physics == imp_physics_zhao_carr_pdf) then ! zhao-carr microphysics
Expand All @@ -615,8 +592,9 @@ subroutine GFS_suite_interstitial_3_run (im, levs, nn, cscnv, &
elseif (imp_physics == imp_physics_thompson) then
do k=1,levs
do i=1,im
clw(i,k,1) = gq0(i,k,ntiw) ! ice
clw(i,k,2) = gq0(i,k,ntcw) ! water
clw(i,k,1) = gq0(i,k,ntiw) ! ice
clw(i,k,2) = gq0(i,k,ntcw) ! water
save_tcp(i,k) = gt0(i,k)
enddo
enddo
if(ltaerosol) then
Expand All @@ -632,15 +610,7 @@ subroutine GFS_suite_interstitial_3_run (im, levs, nn, cscnv, &
clw(i,k,2) = gq0(i,k,ntcw) ! water
enddo
enddo
else ! if_ntcw
!GF* never executed unless imp_physics = imp_physics_zhao_carr or imp_physics_zhao_carr_pdf
! do i=1,im
! psautco_l(i) = Model%psautco(1)*work1(i) + Model%psautco(2)*work2(i)
! prautco_l(i) = Model%prautco(1)*work1(i) + Model%prautco(2)*work2(i)
! enddo
!*GF
rhc(:,:) = 1.0
endif ! end if_ntcw
endif

end subroutine GFS_suite_interstitial_3_run

Expand All @@ -662,9 +632,10 @@ end subroutine GFS_suite_interstitial_4_finalize
subroutine GFS_suite_interstitial_4_run (im, levs, ltaerosol, cplchm, tracers_total, ntrac, ntcw, ntiw, ntclamt, &
ntrw, ntsw, ntrnc, ntsnc, ntgl, ntgnc, ntlnc, ntinc, nn, imp_physics, imp_physics_gfdl, imp_physics_thompson, &
imp_physics_zhao_carr, imp_physics_zhao_carr_pdf, dtf, save_qc, save_qi, con_pi, &
gq0, clw, dqdti, imfdeepcnv, imfdeepcnv_gf, errmsg, errflg)
gq0, clw, prsl, save_tcp, con_rd, nwfa, spechum, dqdti, imfdeepcnv, imfdeepcnv_gf, errmsg, errflg)

use machine, only: kind_phys
use module_mp_thompson_make_number_concentrations, only: make_IceNumber, make_DropletNumber

implicit none

Expand All @@ -683,6 +654,11 @@ subroutine GFS_suite_interstitial_4_run (im, levs, ltaerosol, cplchm, tracers_to

real(kind=kind_phys), dimension(im,levs,ntrac), intent(inout) :: gq0
real(kind=kind_phys), dimension(im,levs,nn), intent(inout) :: clw
real(kind=kind_phys), dimension(im,levs), intent(in) :: prsl
real(kind=kind_phys), intent(in) :: con_rd
real(kind=kind_phys), dimension(:,:), intent(in) :: nwfa, save_tcp
real(kind=kind_phys), dimension(im,levs), intent(in) :: spechum

! dqdti may not be allocated
real(kind=kind_phys), dimension(:,:), intent(inout) :: dqdti

Expand All @@ -693,10 +669,12 @@ subroutine GFS_suite_interstitial_4_run (im, levs, ltaerosol, cplchm, tracers_to
! local variables
integer :: i,k,n,tracers

real(kind=kind_phys) :: liqm, icem

liqm = 4./3.*con_pi*1.e-12
icem = 4./3.*con_pi*3.2768*1.e-14*890.
real(kind=kind_phys), dimension(im,levs) :: rho_dryair
real(kind=kind_phys), dimension(im,levs) :: qv_mp !< kg kg-1 (dry mixing ratio)
real(kind=kind_phys), dimension(im,levs) :: qc_mp !< kg kg-1 (dry mixing ratio)
real(kind=kind_phys), dimension(im,levs) :: qi_mp !< kg kg-1 (dry mixing ratio)
real(kind=kind_phys), dimension(im,levs) :: nc_mp !< kg-1 (dry mixing ratio)
real(kind=kind_phys), dimension(im,levs) :: ni_mp !< kg-1 (dry mixing ratio)

! Initialize CCPP error handling variables
errmsg = ''
Expand Down Expand Up @@ -729,32 +707,42 @@ subroutine GFS_suite_interstitial_4_run (im, levs, ltaerosol, cplchm, tracers_to
imp_physics == imp_physics_zhao_carr_pdf .or. &
imp_physics == imp_physics_gfdl) then
gq0(1:im,:,ntcw) = clw(1:im,:,1) + clw(1:im,:,2)

elseif (ntiw > 0) then
do k=1,levs
do i=1,im
gq0(i,k,ntiw) = clw(i,k,1) ! ice
gq0(i,k,ntcw) = clw(i,k,2) ! water
enddo
enddo
! if (imp_physics == imp_physics_thompson) then
if (imp_physics == imp_physics_thompson .and. imfdeepcnv /= imfdeepcnv_gf) then
if (ltaerosol) then
do k=1,levs
do i=1,im
gq0(i,k,ntlnc) = gq0(i,k,ntlnc) &
+ max(0.0, (clw(i,k,2)-save_qc(i,k))) / liqm
gq0(i,k,ntinc) = gq0(i,k,ntinc) &
+ max(0.0, (clw(i,k,1)-save_qi(i,k))) / icem
enddo
enddo
else
do k=1,levs
do i=1,im
gq0(i,k,ntinc) = gq0(i,k,ntinc) &
+ max(0.0, (clw(i,k,1)-save_qi(i,k))) / icem
enddo
enddo
endif

if (imp_physics == imp_physics_thompson .and. (ntlnc>0 .or. ntinc>0)) then
do k=1,levs
do i=1,im
!> - Density of air in kg m-3
rho_dryair(i,k) = prsl(i,k)/(con_rd*save_tcp(i,k))
!> - Convert specific humidity to dry mixing ratio
qv_mp(i,k) = spechum(i,k)/(1.0_kind_phys-spechum(i,k))
if (ntlnc>0) then
!> - Convert moist mixing ratio to dry mixing ratio
qc_mp(i,k) = save_qc(i,k)/(1.0_kind_phys-spechum(i,k))
!> - Convert number concentration from moist to dry
nc_mp(i,k) = gq0(i,k,ntlnc)/(1.0_kind_phys-spechum(i,k))
nc_mp(i,k) = nc_mp(i,k) + max(0.0, make_DropletNumber(qc_mp(i,k) * rho_dryair(i,k), nwfa(i,k)) * (1.0/rho_dryair(i,k)))
!> - Convert number concentrations from dry to moist
gq0(i,k,ntlnc) = nc_mp(i,k)/(1.0_kind_phys+qv_mp(i,k))
endif
if (ntinc>0) then
!> - Convert moist mixing ratio to dry mixing ratio
qi_mp(i,k) = save_qi(i,k)/(1.0_kind_phys-spechum(i,k))
!> - Convert number concentration from moist to dry
ni_mp(i,k) = gq0(i,k,ntinc)/(1.0_kind_phys-spechum(i,k))
ni_mp(i,k) = ni_mp(i,k) + max(0.0, make_IceNumber(qi_mp(i,k) * rho_dryair(i,k), save_tcp(i,k)) * (1.0/rho_dryair(i,k)))
!> - Convert number concentrations from dry to moist
gq0(i,k,ntinc) = ni_mp(i,k)/(1.0_kind_phys+qv_mp(i,k))
endif
enddo
enddo
endif

else
Expand All @@ -764,6 +752,7 @@ subroutine GFS_suite_interstitial_4_run (im, levs, ltaerosol, cplchm, tracers_to
enddo
enddo
endif ! end if_ntiw

else
do k=1,levs
do i=1,im
Expand Down
63 changes: 63 additions & 0 deletions physics/GFS_suite_interstitial.meta
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1218,6 +1218,15 @@
kind = kind_phys
intent = in
optional = F
[gt0]
standard_name = air_temperature_updated_by_physics
long_name = temperature updated by physics
units = K
dimensions = (horizontal_dimension,vertical_dimension)
type = real
kind = kind_phys
intent = in
optional = F
[gq0]
standard_name = tracer_concentration_updated_by_physics
long_name = tracer concentration updated by physics
Expand Down Expand Up @@ -1432,6 +1441,15 @@
kind = kind_phys
intent = inout
optional = F
[save_tcp]
standard_name = air_temperature_save_from_cumulus_paramterization
long_name = air temperature after cumulus parameterization
units = K
dimensions = (horizontal_dimension,vertical_dimension)
type = real
kind = kind_phys
intent = inout
optional = F
[errmsg]
standard_name = ccpp_error_message
long_name = error message for error handling in CCPP
Expand Down Expand Up @@ -1692,6 +1710,51 @@
kind = kind_phys
intent = inout
optional = F
[prsl]
standard_name = air_pressure
long_name = mean layer pressure
units = Pa
dimensions = (horizontal_dimension,vertical_dimension)
type = real
kind = kind_phys
intent = in
optional = F
[save_tcp]
standard_name = air_temperature_save_from_cumulus_paramterization
long_name = air temperature after cumulus parameterization
units = K
dimensions = (horizontal_dimension,vertical_dimension)
type = real
kind = kind_phys
intent = in
optional = F
[con_rd]
standard_name = gas_constant_dry_air
long_name = ideal gas constant for dry air
units = J kg-1 K-1
dimensions = ()
type = real
kind = kind_phys
intent = in
optional = F
[nwfa]
standard_name = water_friendly_aerosol_number_concentration
long_name = number concentration of water-friendly aerosols
units = kg-1
dimensions = (horizontal_dimension,vertical_dimension)
type = real
kind = kind_phys
intent = in
optional = F
[spechum]
standard_name = water_vapor_specific_humidity
long_name = water vapor specific humidity
units = kg kg-1
dimensions = (horizontal_dimension,vertical_dimension)
type = real
kind = kind_phys
intent = inout
optional = F
[dqdti]
standard_name = instantaneous_water_vapor_specific_humidity_tendency_due_to_convection
long_name = instantaneous moisture tendency due to convection
Expand Down
26 changes: 0 additions & 26 deletions physics/cu_gf_driver.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -9,7 +9,6 @@ module cu_gf_driver
use machine , only: kind_phys
use cu_gf_deep, only: cu_gf_deep_run,neg_check,autoconv,aeroevap,fct1d3
use cu_gf_sh , only: cu_gf_sh_run
use module_mp_thompson_make_number_concentrations, only: make_IceNumber, make_DropletNumber

implicit none

Expand Down Expand Up @@ -74,7 +73,6 @@ subroutine cu_gf_driver_run(ntracer,garea,im,ix,km,dt,cactiv, &
us,vs,t2di,w,qv2di_spechum,p2di,psuri, &
hbot,htop,kcnv,xland,hfx2,qfx2,cliw,clcw, &
pbl,ud_mf,dd_mf,dt_mf,cnvw_moist,cnvc,imfshalcnv, &
nwfa,con_rd,gq0,ntinc,ntlnc,imp_physics,imp_physics_thompson, &
errmsg,errflg)
!-------------------------------------------------------------
implicit none
Expand Down Expand Up @@ -126,12 +124,6 @@ subroutine cu_gf_driver_run(ntracer,garea,im,ix,km,dt,cactiv, &
real(kind=kind_phys), dimension( im ),intent(in) :: garea
real(kind=kind_phys), intent(in ) :: dt

! additional variables for number concentrations
real(kind=kind_phys), intent(in) :: nwfa(1:im,1:km)
real(kind=kind_phys), intent(in) :: con_rd
real(kind=kind_phys), dimension(im,km,ntracer), intent(inout) :: gq0
integer, intent(in) :: imp_physics,imp_physics_thompson,ntlnc,ntinc

integer, intent(in ) :: imfshalcnv
character(len=*), intent(out) :: errmsg
integer, intent(out) :: errflg
Expand Down Expand Up @@ -826,26 +818,8 @@ subroutine cu_gf_driver_run(ntracer,garea,im,ix,km,dt,cactiv, &
cliw(i,k) = max(0.,cliw(i,k) + tem)
endif

!
!> calculate cloud water and cloud ice number concentrations
!
rho_dryar(i,k) = p2di(i,k)/(con_rd*t(i,k)) ! Density of dry air in kg m-3
if (imp_physics == imp_physics_thompson) then
if ((tem*tem1)>1.e-5) then
gq0(i,k,ntinc) = max(0., gq0(i,k,ntinc) + &
make_IceNumber(tem*tem1*rho_dryar(i,k), t(i,k)) * &
(1/rho_dryar(i,k)))
end if
if ((tem*(1-tem1))>1.e-5) then
gq0(i,k,ntlnc) = max(0., gq0(i,k,ntlnc) + &
make_DropletNumber(tem*(1-tem1)*rho_dryar(i,k), nwfa(i,k)) &
* (1/rho_dryar(i,k)))
end if
end if

enddo


gdc(i,1,10)=forcing(i,1)
gdc(i,2,10)=forcing(i,2)
gdc(i,3,10)=forcing(i,3)
Expand Down
Loading