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

Jump to bottom

Generalize computation of moist specific heats #221

Merged
merged 4 commits into from
Sep 20, 2021
Merged

Generalize computation of moist specific heats #221

Changes from 2 commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
a61a8de
Generalize computation of moist specific heats
spencerkclark Sep 13, 2021
fc270f0
Consolidate common operations
spencerkclark Sep 13, 2021
39eba1f
Implement special case for nwat = 2
spencerkclark Sep 15, 2021
9eadfca
Merge branch 'master' into general-moist-specific-heat
spencerkclark Sep 16, 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
107 changes: 61 additions & 46 deletions FV3/gfsphysics/GFS_layer/GFS_physics_driver.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -5462,11 +5462,11 @@ subroutine GFS_physics_driver &
nwat = Statein%nwat

if (Statein%dycore_hydrostatic) then
call moist_cp_nwat6(Statein%qgrs(1:im,1:levs,1:nwat), Stateout%gq0(1:im,1:levs,1:nwat), &
call moist_cp(Statein%qgrs(1:im,1:levs,1:nwat), Stateout%gq0(1:im,1:levs,1:nwat), &
Statein%prsi(1:im,1:levs+1), im, levs, nwat, 1, Model%ntcw, Model%ntiw, &
Model%ntrw, Model%ntsw, Model%ntgl, specific_heat)
else
call moist_cv_nwat6(Statein%qgrs(1:im,1:levs,1:nwat), Stateout%gq0(1:im,1:levs,1:nwat), &
call moist_cv(Statein%qgrs(1:im,1:levs,1:nwat), Stateout%gq0(1:im,1:levs,1:nwat), &
Statein%prsi(1:im,1:levs+1), im, levs, nwat, 1, Model%ntcw, Model%ntiw, &
Model%ntrw, Model%ntsw, Model%ntgl, specific_heat)
endif
Expand Down Expand Up @@ -5677,7 +5677,38 @@ subroutine moist_bud2(im,ix,ix2,levs,me,kdt,grav,dtp,delp,rain, &

end subroutine moist_bud2

subroutine moist_cv_nwat6(initial_dynamics_q, physics_q, pressure_on_interfaces, &
subroutine compute_q_liquid(im, levs, nwat, ntcw, ntrw, new_dynamics_q, q_liquid)
spencerkclark marked this conversation as resolved.
Show resolved Hide resolved
integer, intent(in) :: im, levs, nwat, ntcw, ntrw
real(kind=kind_phys), dimension(1:im,1:levs,1:nwat), intent(in) :: new_dynamics_q
real(kind=kind_phys), dimension(1:im,1:levs), intent(out) :: q_liquid

q_liquid = 0.0
if (ntcw .gt. 0) then
q_liquid = q_liquid + new_dynamics_q(:,:,ntcw)
endif
if (ntrw .gt. 0) then
q_liquid = q_liquid + new_dynamics_q(:,:,ntrw)
endif
end subroutine compute_q_liquid

subroutine compute_q_ice(im, levs, nwat, ntiw, ntsw, ntgl, new_dynamics_q, q_ice)
integer, intent(in) :: im, levs, nwat, ntiw, ntsw, ntgl
real(kind=kind_phys), dimension(1:im,1:levs,1:nwat), intent(in) :: new_dynamics_q
real(kind=kind_phys), dimension(1:im,1:levs), intent(out) :: q_ice

q_ice = 0.0
if (ntiw .gt. 0) then
q_ice = q_ice + new_dynamics_q(:,:,ntiw)
endif
if (ntsw .gt. 0) then
q_ice = q_ice + new_dynamics_q(:,:,ntsw)
endif
if (ntgl .gt. 0) then
q_ice = q_ice + new_dynamics_q(:,:,ntgl)
endif
end subroutine compute_q_ice

subroutine moist_cv(initial_dynamics_q, physics_q, pressure_on_interfaces, &
spencerkclark marked this conversation as resolved.
Show resolved Hide resolved
im, levs, nwat, ntqv, ntcw, ntiw, ntrw, ntsw, ntgl, cvm)
integer, intent(in) :: im, levs, nwat, ntqv, ntcw, ntiw, ntrw, ntsw, ntgl
real(kind=kind_phys), dimension(1:im,1:levs,1:nwat), intent(in) :: initial_dynamics_q, physics_q
Expand All @@ -5695,35 +5726,27 @@ subroutine moist_cv_nwat6(initial_dynamics_q, physics_q, pressure_on_interfaces,
real(kind=kind_phys) :: c_liq = 4.1855e+3 ! Hard-coded in fv_mapz.F90
real(kind=kind_phys) :: c_ice = 1972.0 ! Hard-coded in fv_mapz.F90

! fv_mapz.moist_cv defines branches for using other moist tracer configurations.
! For simplicity we choose not to replicate that behavior here, since we have
! only run in one tracer configuration (nwat = 6) so far. We also do not implement
! the branch of code that is run if the compiler directive MULTI_GASES is defined.
! In those cases we default to using the specific heat at constant volume for dry
! air, and emit a warning.
! We do not currently implement the branch of code that is run if the compiler
! directive MULTI_GASES is defined. In those cases we default to using the
! specific heat at constant volume for dry air, and emit a warning.
#ifdef MULTI_GASES
call mpp_error (NOTE, 'GFS_physics_driver::moist_cv - moist_cv for tracer configuration not implemented; using default cv_air for t_dt diagnostics')
cvm = cv_air
#else
if (nwat /= 6) then
call mpp_error (NOTE, 'GFS_physics_driver::moist_cv - moist_cv for tracer configuration not implemented; using default cv_air for t_dt diagnostics')
cvm = cv_air
else
call physics_to_dycore_mass_fraction(initial_dynamics_q, physics_q, &
pressure_on_interfaces, im, levs, nwat, new_dynamics_q)
call physics_to_dycore_mass_fraction(initial_dynamics_q, physics_q, &
pressure_on_interfaces, im, levs, nwat, new_dynamics_q)

q_vapor = new_dynamics_q(:,:,ntqv)
q_liquid = new_dynamics_q(:,:,ntcw) + new_dynamics_q(:,:,ntrw)
q_ice = new_dynamics_q(:,:,ntiw) + new_dynamics_q(:,:,ntsw) + new_dynamics_q(:,:,ntgl)
q_dry_air = 1.0 - q_vapor - q_liquid - q_ice
q_vapor = new_dynamics_q(:,:,ntqv)
call compute_q_liquid(im, levs, nwat, ntcw, ntrw, new_dynamics_q, q_liquid)
call compute_q_ice(im, levs, nwat, ntiw, ntsw, ntgl, new_dynamics_q, q_ice)
q_dry_air = 1.0 - q_vapor - q_liquid - q_ice

! By definition now, the weights sum to 1.0.
cvm = cv_air * q_dry_air + cv_vap * q_vapor + c_liq * q_liquid + c_ice * q_ice
endif
! By definition now, the weights sum to 1.0.
cvm = cv_air * q_dry_air + cv_vap * q_vapor + c_liq * q_liquid + c_ice * q_ice
#endif
end subroutine moist_cv_nwat6
end subroutine moist_cv

subroutine moist_cp_nwat6(initial_dynamics_q, physics_q, pressure_on_interfaces, &
subroutine moist_cp(initial_dynamics_q, physics_q, pressure_on_interfaces, &
spencerkclark marked this conversation as resolved.
Show resolved Hide resolved
im, levs, nwat, ntqv, ntcw, ntiw, ntrw, ntsw, ntgl, cpm)
integer, intent(in) :: im, levs, nwat, ntqv, ntcw, ntiw, ntrw, ntsw, ntgl
real(kind=kind_phys), dimension(1:im,1:levs,1:nwat), intent(in) :: initial_dynamics_q, physics_q
Expand All @@ -5741,33 +5764,25 @@ subroutine moist_cp_nwat6(initial_dynamics_q, physics_q, pressure_on_interfaces,
real(kind=kind_phys) :: c_liq = 4.1855e+3 ! Hard-coded in fv_mapz.F90
real(kind=kind_phys) :: c_ice = 1972.0 ! Hard-coded in fv_mapz.F90

! fv_mapz.moist_cp defines branches for using other moist tracer configurations.
! For simplicity we choose not to replicate that behavior here, since we have
! only run in one tracer configuration (nwat = 6) so far. We also do not implement
! the branch of code that is run if the compiler directive MULTI_GASES is defined.
! In those cases we default to using the specific heat at constant volume for dry
! air, and emit a warning.
! We do not currently implement the branch of code that is run if the compiler
! directive MULTI_GASES is defined. In those cases we default to using the
! specific heat at constant volume for dry air, and emit a warning.
#ifdef MULTI_GASES
call mpp_error (NOTE, 'GFS_physics_driver::moist_cp - moist_cp for tracer configuration not implemented; using default cp_air for t_dt diagnostics')
cpm = cp_air
#else
if (nwat /= 6) then
call mpp_error (NOTE, 'GFS_physics_driver::moist_cp - moist_cp for tracer configuration not implemented; using default cp_air for t_dt diagnostics')
cpm = cp_air
else
call physics_to_dycore_mass_fraction(initial_dynamics_q, physics_q, &
pressure_on_interfaces, im, levs, nwat, new_dynamics_q)

q_vapor = new_dynamics_q(:,:,ntqv)
q_liquid = new_dynamics_q(:,:,ntcw) + new_dynamics_q(:,:,ntrw)
q_ice = new_dynamics_q(:,:,ntiw) + new_dynamics_q(:,:,ntsw) + new_dynamics_q(:,:,ntgl)
q_dry_air = 1.0 - q_vapor - q_liquid - q_ice

! By definition now, the weights sum to 1.0.
cpm = cp_air * q_dry_air + cp_vap * q_vapor + c_liq * q_liquid + c_ice * q_ice
endif
call physics_to_dycore_mass_fraction(initial_dynamics_q, physics_q, &
pressure_on_interfaces, im, levs, nwat, new_dynamics_q)

q_vapor = new_dynamics_q(:,:,ntqv)
call compute_q_liquid(im, levs, nwat, ntcw, ntrw, new_dynamics_q, q_liquid)
call compute_q_ice(im, levs, nwat, ntiw, ntsw, ntgl, new_dynamics_q, q_ice)
q_dry_air = 1.0 - q_vapor - q_liquid - q_ice

! By definition now, the weights sum to 1.0.
cpm = cp_air * q_dry_air + cp_vap * q_vapor + c_liq * q_liquid + c_ice * q_ice
#endif
end subroutine moist_cp_nwat6
end subroutine moist_cp

subroutine physics_to_dycore_mass_fraction(initial_dynamics_q, physics_q, &
pressure_on_interfaces, im, levs, nwat, new_dynamics_q)
Expand Down