Skip to content

Commit

Permalink
Fixes conservation problem for non-water species
Browse files Browse the repository at this point in the history 
All non-water tracers are kept as dry mixing ratios, but CLUBB and the
gravity wave parameterizations assume wet mixing ratios, and break
conservation when they are called. This fix converts dry mixing ratios
to wet mixing ratios prior to calling the routines that need wet mixing
ratios. The code fixes follow similar fixes from CESM2, with additional
modifications for E3SM-specific components.

Fixes #2704

[Non-BFB] - Non Bit-For-Bit

See confluence for a more detailed description about these tags.
  • Loading branch information
@beharrop
beharrop committed Jun 4, 2019
1 parent 25c9436 commit 69c6cda
Show file tree
Hide file tree
Showing 7 changed files with 118 additions and 37 deletions.
19 changes: 19 additions & 0 deletions components/cam/src/control/cam_history.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4064,6 +4064,13 @@ subroutine h_define (t, restart)
'h_define: cannot define units for '//trim(fname_tmp))
end if

str = tape(t)%hlist(f)%field%mixing_ratio
if (len_trim(str) > 0) then
ierr=pio_put_att (tape(t)%File, varid, 'mixing_ratio', trim(str))
call cam_pio_handle_error(ierr, &
'h_define: cannot define mixing_ratio for '//trim(fname_tmp))
end if

str = tape(t)%hlist(f)%field%long_name
ierr=pio_put_att (tape(t)%File, varid, 'long_name', trim(str))
call cam_pio_handle_error(ierr, &
Expand Down Expand Up @@ -4813,6 +4820,7 @@ subroutine addfld_nd(fname, dimnames, avgflag, units, long_name, &
!-----------------------------------------------------------------------
use cam_history_support, only: fillvalue, hist_coord_find_levels
use cam_grid_support, only: cam_grid_id
use constituents, only: pcnst, cnst_get_ind, cnst_get_type_byind

!
! Arguments
Expand All @@ -4836,9 +4844,11 @@ subroutine addfld_nd(fname, dimnames, avgflag, units, long_name, &
!
character(len=max_fieldname_len) :: fname_tmp ! local copy of fname
character(len=128) :: errormsg
character(len=3) :: mixing_ratio
type(master_entry), pointer :: listentry

integer :: dimcnt
integer :: idx

if (htapes_defined) then
call endrun ('ADDFLD: Attempt to add field '//trim(fname)//' after history files set')
Expand Down Expand Up @@ -4871,6 +4881,14 @@ subroutine addfld_nd(fname, dimnames, avgflag, units, long_name, &
call endrun ('ADDFLD: '//fname//' already on list')
end if

! If the field is an advected constituent determine whether its concentration
! is based on dry or wet air.
call cnst_get_ind(fname_tmp, idx, abort=.false.)
mixing_ratio = ''
if (idx > 0) then
mixing_ratio = cnst_get_type_byind(idx)
end if

!
! Add field to Master Field List arrays fieldn and iflds
!
Expand All @@ -4879,6 +4897,7 @@ subroutine addfld_nd(fname, dimnames, avgflag, units, long_name, &
listentry%field%long_name = long_name
listentry%field%numlev = 1 ! Will change if lev or ilev in shape
listentry%field%units = units
listentry%field%mixing_ratio = mixing_ratio
listentry%field%meridional_complement = -1
listentry%field%zonal_complement = -1
listentry%htapeindx(:) = -1
Expand Down
1 change: 1 addition & 0 deletions components/cam/src/control/cam_history_support.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -121,6 +121,7 @@ module cam_history_support
character(len=max_fieldname_len) :: name ! field name
character(len=max_chars) :: long_name ! long name
character(len=max_chars) :: units ! units
character(len=3) :: mixing_ratio ! 'dry' or 'wet'
character(len=max_chars) :: sampling_seq ! sampling sequence - if not every timestep, how often field is sampled
! (i.e., how often "outfld" is called): every other; only during LW/SW
! radiation calcs; etc.
Expand Down
38 changes: 33 additions & 5 deletions components/cam/src/physics/cam/clubb_intr.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -903,15 +903,15 @@ subroutine clubb_tend_cam( &

use physics_types, only: physics_state, physics_ptend, &
physics_state_copy, physics_ptend_init, &
physics_ptend_sum
physics_ptend_sum, set_dry_to_wet

use physics_update_mod, only: physics_update

use physics_buffer, only: pbuf_get_index, pbuf_old_tim_idx, pbuf_get_field, &
pbuf_set_field, physics_buffer_desc

use ppgrid, only: pver, pverp, pcols
use constituents, only: cnst_get_ind
use constituents, only: cnst_get_ind, cnst_type
use camsrfexch, only: cam_in_t
use ref_pres, only: top_lev => trop_cloud_top_lev
use time_manager, only: is_first_step
Expand Down Expand Up @@ -1257,6 +1257,9 @@ subroutine clubb_tend_cam( &

call physics_state_copy(state,state1)

! constituents are all treated as wet mmr by clubb
call set_dry_to_wet(state1)

if (micro_do_icesupersat) then
naai_idx = pbuf_get_index('NAAI')
call pbuf_get_field(pbuf, naai_idx, naai)
Expand Down Expand Up @@ -2202,6 +2205,18 @@ subroutine clubb_tend_cam( &
call physics_ptend_sum(ptend_loc,ptend_all,ncol)
call physics_update(state1,ptend_loc,hdtime)

! ptend_all now has all accumulated tendencies. Convert the tendencies for the
! dry constituents to dry air basis.
do ixind = 1, pcnst
if (lq(ixind) .and. cnst_type(ixind).eq.'dry') then
do k = 1, pver
do i = 1, ncol
ptend_all%q(i,k,ixind) = ptend_all%q(i,k,ixind)*state1%pdel(i,k)/state1%pdeldry(i,k)
end do
end do
end if
end do

! ------------------------------------------------- !
! Diagnose relative cloud water variance !
! ------------------------------------------------- !
Expand Down Expand Up @@ -2548,10 +2563,12 @@ subroutine clubb_surface (state, ptend, ztodt, cam_in, ustar, obklen)
! None
!-------------------------------------------------------------------------------

use physics_types, only: physics_state, physics_ptend, physics_ptend_init
use physics_types, only: physics_state, physics_ptend, &
physics_ptend_init, &
set_dry_to_wet, set_wet_to_dry
use physconst, only: gravit, zvir, latvap
use ppgrid, only: pver, pcols
use constituents, only: pcnst, cnst_get_ind
use constituents, only: pcnst, cnst_get_ind, cnst_type
use camsrfexch, only: cam_in_t

implicit none
Expand All @@ -2560,7 +2577,7 @@ subroutine clubb_surface (state, ptend, ztodt, cam_in, ustar, obklen)
! Input Auguments !
! --------------- !

type(physics_state), intent(in) :: state ! Physics state variables
type(physics_state), intent(inout) :: state ! Physics state variables
type(cam_in_t), intent(in) :: cam_in

real(r8), intent(in) :: ztodt ! 2 delta-t [ s ]
Expand Down Expand Up @@ -2604,6 +2621,8 @@ subroutine clubb_surface (state, ptend, ztodt, cam_in, ustar, obklen)
! Main Computation Begins !
! ----------------------- !

! Assume 'wet' mixing ratios in surface diffusion code.
call set_dry_to_wet(state)

call cnst_get_ind('Q',ixq)

Expand Down Expand Up @@ -2635,6 +2654,15 @@ subroutine clubb_surface (state, ptend, ztodt, cam_in, ustar, obklen)
enddo

ptend%q(:ncol,:pver,:) = (ptend%q(:ncol,:pver,:) - state%q(:ncol,:pver,:)) * rztodt

! Convert tendencies of dry constituents to dry basis.
do m = 1,pcnst
if (cnst_type(m).eq.'dry') then
ptend%q(:ncol,:pver,m) = ptend%q(:ncol,:pver,m)*state%pdel(:ncol,:pver)/state%pdeldry(:ncol,:pver)
endif
end do
! convert wet mmr back to dry before conservation check
call set_wet_to_dry(state)

return

Expand Down
2 changes: 1 addition & 1 deletion components/cam/src/physics/cam/constituents.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -51,7 +51,7 @@ module constituents
real(r8), public :: cnst_cp (pcnst) ! specific heat at constant pressure (J/kg/K)
real(r8), public :: cnst_cv (pcnst) ! specific heat at constant volume (J/kg/K)
real(r8), public :: cnst_mw (pcnst) ! molecular weight (kg/kmole)
character*3, public :: cnst_type(pcnst) ! wet or dry mixing ratio
character*3, public, protected :: cnst_type(pcnst) ! wet or dry mixing ratio
character*5, public :: cnst_molec(pcnst) ! major or minor species molecular diffusion
real(r8), public :: cnst_rgas(pcnst) ! gas constant ()
real(r8), public :: qmin (pcnst) ! minimum permitted constituent concentration (kg/kg)
Expand Down
62 changes: 42 additions & 20 deletions components/cam/src/physics/cam/gw_drag.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -558,6 +558,8 @@ subroutine gw_tend(state, sgh, pbuf, dt, ptend, cam_in)
!-----------------------------------------------------------------------
! Interface for multiple gravity wave drag parameterization.
!-----------------------------------------------------------------------
use physics_types, only: physics_state_copy, set_dry_to_wet
use constituents, only: cnst_type
use physics_buffer, only: physics_buffer_desc, pbuf_get_field
use camsrfexch, only: cam_in_t
! Location-dependent cpair
Expand All @@ -578,10 +580,13 @@ subroutine gw_tend(state, sgh, pbuf, dt, ptend, cam_in)
type(cam_in_t), intent(in) :: cam_in

!---------------------------Local storage-------------------------------

type(physics_state) :: state1 ! Local copy of state variable

integer :: lchnk ! chunk identifier
integer :: ncol ! number of atmospheric columns

integer :: k ! loop index
integer :: i, k ! loop indices

real(r8) :: ttgw(state%ncol,pver) ! temperature tendency
real(r8) :: utgw(state%ncol,pver) ! zonal wind tendency
Expand Down Expand Up @@ -658,23 +663,29 @@ subroutine gw_tend(state, sgh, pbuf, dt, ptend, cam_in)

!------------------------------------------------------------------------

lchnk = state%lchnk
ncol = state%ncol

dse = state%s(:ncol,:)
t = state%t(:ncol,:)
u = state%u(:ncol,:)
v = state%v(:ncol,:)
q = state%q(:ncol,:,:)
pmid = state%pmid(:ncol,:)
pint = state%pint(:ncol,:)
piln = state%lnpint(:ncol,:)
dpm = state%pdel(:ncol,:)
rdpm = state%rpdel(:ncol,:)
zm = state%zm(:ncol,:)
! Make local copy of input state.
call physics_state_copy(state, state1)

! constituents are all treated as wet mmr
call set_dry_to_wet(state1)

lchnk = state1%lchnk
ncol = state1%ncol

dse = state1%s(:ncol,:)
t = state1%t(:ncol,:)
u = state1%u(:ncol,:)
v = state1%v(:ncol,:)
q = state1%q(:ncol,:,:)
pmid = state1%pmid(:ncol,:)
pint = state1%pint(:ncol,:)
piln = state1%lnpint(:ncol,:)
dpm = state1%pdel(:ncol,:)
rdpm = state1%rpdel(:ncol,:)
zm = state1%zm(:ncol,:)

lq = .true.
call physics_ptend_init(ptend, state%psetcols, "Grav_wave_drag", &
call physics_ptend_init(ptend, state1%psetcols, "Grav_wave_drag", &
ls=.true., lu=.true., lv=.true., lq=lq)

! Profiles of background state variables
Expand Down Expand Up @@ -722,13 +733,13 @@ subroutine gw_tend(state, sgh, pbuf, dt, ptend, cam_in)
call pbuf_get_field(pbuf, ttend_dp_idx, ttend_dp)

! Determine wave sources for Beres04 scheme
call gw_beres_src(ncol, pgwv, state%lat(:ncol), u, v, ttend_dp, &
call gw_beres_src(ncol, pgwv, state1%lat(:ncol), u, v, ttend_dp, &
zm, src_level, tend_level, tau, ubm, ubi, xv, yv, c, &
hdepth, maxq0)

! Solve for the drag profile with Beres source spectrum.
call gw_drag_prof(ncol, pgwv, src_level, tend_level, .false., dt, &
state%lat(:ncol), t, ti, pmid, pint, dpm, rdpm, &
state1%lat(:ncol), t, ti, pmid, pint, dpm, rdpm, &
piln, rhoi, nm, ni, ubm, ubi, xv, yv, &
effgw_beres, c, kvtt, q, dse, tau, utgw, vtgw, &
ttgw, qtgw, taucd, egwdffi, gwut, dttdf, dttke)
Expand Down Expand Up @@ -785,7 +796,7 @@ subroutine gw_tend(state, sgh, pbuf, dt, ptend, cam_in)

! Solve for the drag profile with C&M source spectrum.
call gw_drag_prof(ncol, pgwv, src_level, tend_level, .true., dt, &
state%lat(:ncol), t, ti, pmid, pint, dpm, rdpm, &
state1%lat(:ncol), t, ti, pmid, pint, dpm, rdpm, &
piln, rhoi, nm, ni, ubm, ubi, xv, yv, &
effgw_cm, c, kvtt, q, dse, tau, utgw, vtgw, &
ttgw, qtgw, taucd, egwdffi, gwut, dttdf, dttke)
Expand Down Expand Up @@ -838,7 +849,7 @@ subroutine gw_tend(state, sgh, pbuf, dt, ptend, cam_in)

! Solve for the drag profile with orographic sources.
call gw_drag_prof(ncol, 0, src_level, tend_level, .false., dt, &
state%lat(:ncol), t, ti, pmid, pint, dpm, rdpm, &
state1%lat(:ncol), t, ti, pmid, pint, dpm, rdpm, &
piln, rhoi, nm, ni, ubm, ubi, xv, yv, &
effgw_oro, c, kvtt, q, dse, tau, utgw, vtgw, &
ttgw, qtgw, taucd, egwdffi, gwut(:,:,0:0), dttdf, dttke)
Expand Down Expand Up @@ -878,6 +889,17 @@ subroutine gw_tend(state, sgh, pbuf, dt, ptend, cam_in)

end if

! Convert the tendencies for the dry constituents to dry air basis.
do m = 1, pcnst
if (cnst_type(m).eq.'dry') then
do k = 1, pver
do i = 1, ncol
ptend%q(i,k,m) = ptend%q(i,k,m)*state1%pdel(i,k)/state1%pdeldry(i,k)
end do
end do
end if
end do

! Write total temperature tendency to history file
call outfld ('TTGW', ptend%s/cpairv(:,:,lchnk), pcols, lchnk)

Expand Down
4 changes: 2 additions & 2 deletions components/cam/src/physics/cam/physics_types.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -504,7 +504,7 @@ subroutine physics_state_check(state, name)
shr_infnan_posinf, shr_infnan_neginf
use shr_assert_mod, only: shr_assert, shr_assert_in_domain
use physconst, only: pi
use constituents, only: pcnst, qmin
use constituents, only: pcnst

!------------------------------Arguments--------------------------------
! State to check.
Expand Down Expand Up @@ -677,7 +677,7 @@ subroutine physics_state_check(state, name)

! 3-D variables
do m = 1,pcnst
call shr_assert_in_domain(state%q(:ncol,:,m), lt=posinf_r8, ge=qmin(m), &
call shr_assert_in_domain(state%q(:ncol,:,m), lt=posinf_r8, gt=neginf_r8, &
varname="state%q ("//trim(cnst_name(m))//")", msg=msg)
end do

Expand Down
29 changes: 20 additions & 9 deletions components/cam/src/physics/cam/vertical_diffusion.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -443,11 +443,8 @@ subroutine vertical_diffusion_init(pbuf2d)
enddo
end if

if( cnst_get_type_byind(k) .eq. 'wet' ) then
if( vdiff_select( fieldlist_wet, 'q', k ) .ne. '' ) call endrun( vdiff_select( fieldlist_wet, 'q', k ) )
else
if( vdiff_select( fieldlist_dry, 'q', k ) .ne. '' ) call endrun( vdiff_select( fieldlist_dry, 'q', k ) )
endif
! Convert all constituents to wet before doing diffusion.
if( vdiff_select( fieldlist_wet, 'q', k ) .ne. '' ) call endrun( vdiff_select( fieldlist_wet, 'q', k ) )

! ----------------------------------------------- !
! Select constituents for molecular diffusion !
Expand Down Expand Up @@ -650,14 +647,15 @@ subroutine vertical_diffusion_tend( &
!---------------------------------------------------- !
use physics_buffer, only : physics_buffer_desc, pbuf_get_field, pbuf_set_field
use physics_types, only : physics_state, physics_ptend, physics_ptend_init
use physics_types, only : set_dry_to_wet, set_wet_to_dry
use cam_history, only : outfld

use trb_mtn_stress, only : compute_tms
use eddy_diff, only : compute_eddy_diff
use hb_diff, only : compute_hb_diff
use wv_saturation, only : qsat
use molec_diff, only : compute_molec_diff, vd_lu_qdecomp
use constituents, only : qmincg, qmin
use constituents, only : qmincg, qmin, cnst_type
use diffusion_solver, only : compute_vdiff, any, operator(.not.)
use physconst, only : cpairv, rairv !Needed for calculation of upward H flux
use time_manager, only : get_nstep
Expand All @@ -669,7 +667,7 @@ subroutine vertical_diffusion_tend( &
! Input Arguments !
! --------------- !

type(physics_state), intent(in) :: state ! Physics state variables
type(physics_state), intent(inout) :: state ! Physics state variables

real(r8), intent(in) :: taux(pcols) ! x surface stress [ N/m2 ]
real(r8), intent(in) :: tauy(pcols) ! y surface stress [ N/m2 ]
Expand Down Expand Up @@ -821,6 +819,9 @@ subroutine vertical_diffusion_tend( &
! Main Computation Begins !
! ----------------------- !

! Assume 'wet' mixing ratios in diffusion code.
call set_dry_to_wet(state)

rztodt = 1._r8 / ztodt
lchnk = state%lchnk
ncol = state%ncol
Expand Down Expand Up @@ -1117,7 +1118,7 @@ subroutine vertical_diffusion_tend( &
if (prog_modal_aero) then

! Modal aerosol species not diffused, so just add the explicit surface fluxes to the
! lowest layer
! lowest layer. **NOTE** This code assumes wet mmr.

tmp1(:ncol) = ztodt * gravit * state%rpdel(:ncol,pver)
do m = 1, pmam_ncnst
Expand Down Expand Up @@ -1217,6 +1218,16 @@ subroutine vertical_diffusion_tend( &
ptend%u(:ncol,:) = ( u_tmp(:ncol,:) - state%u(:ncol,:) ) * rztodt
ptend%v(:ncol,:) = ( v_tmp(:ncol,:) - state%v(:ncol,:) ) * rztodt
ptend%q(:ncol,:pver,:) = ( q_tmp(:ncol,:pver,:) - state%q(:ncol,:pver,:) ) * rztodt

! Convert tendencies of dry constituents to dry basis.
do m = 1,pcnst
if (cnst_type(m).eq.'dry') then
ptend%q(:ncol,:pver,m) = ptend%q(:ncol,:pver,m)*state%pdel(:ncol,:pver)/state%pdeldry(:ncol,:pver)
endif
end do
! convert wet mmr back to dry before conservation check
call set_wet_to_dry(state)

slten(:ncol,:) = ( sl(:ncol,:) - sl_prePBL(:ncol,:) ) * rztodt
qtten(:ncol,:) = ( qt(:ncol,:) - qt_prePBL(:ncol,:) ) * rztodt

Expand Down Expand Up @@ -1311,7 +1322,7 @@ subroutine vertical_diffusion_tend( &
'MASSCHECK vert diff : nstep,lon,lat,mass1,mass2,sum3,sflx,rel-diff : ', &
trim(cnst_name(m)), ' : ', nstep, state%lon(i)*180._r8/pi, state%lat(i)*180._r8/pi, &
sum1, sum2, sum3, sflx, abs(sum2-sum1)/sum1
call endrun('vertical_diffusion_tend : mass not conserved' )
!xxx call endrun('vertical_diffusion_tend : mass not conserved' )
endif
endif
enddo col_loop
Expand Down

0 comments on commit 69c6cda

Please sign in to comment.