enthalpy flux implemented with physics buffers

src/control/camsrfexch.F90

@@ -27,7 +27,7 @@ module camsrfexch
   public atm2hub_deallocate
   public hub2atm_deallocate
   public cam_export
-
+  public get_prec_vars
   ! Public data types
   public cam_out_t                  ! Data from atmosphere
   public cam_in_t                   ! Merged surface data
@@ -415,7 +415,8 @@ subroutine cam_export(state,cam_in,cam_out,pbuf)
    use physics_buffer,   only: pbuf_get_index, pbuf_get_field, physics_buffer_desc, pbuf_set_field
    use rad_constituents, only: rad_cnst_get_gas
    use cam_control_mod,  only: simple_phys
-   use air_composition,  only: hliq_idx, hice_idx, fliq_idx, fice_idx, compute_enthalpy_flux
+   use air_composition,  only: hliq_idx, hice_idx, fliq_idx, fice_idx
+   use air_composition,  only: compute_enthalpy_flux, num_enthalpy_vars
    implicit none
 
    ! Input arguments
@@ -438,9 +439,12 @@ subroutine cam_export(state,cam_in,cam_out,pbuf)
 
    real(r8), pointer :: o3_ptr(:,:), srf_o3_ptr(:)
    real(r8), pointer :: lightning_ptr(:)
-
-   real(r8), dimension(:,:), pointer :: enthalpy_prec_bc, enthalpy_prec_ac
-   real(r8), dimension(pcols)        :: fliq_tot, fice_tot
+   !
+   ! enthalpy variables (if applicable)
+   !
+   real(r8), dimension(:,:), pointer            :: enthalpy_prec_ac
+   real(r8), dimension(pcols)                   :: fliq_tot, fice_tot
+   real(r8), dimension(pcols,num_enthalpy_vars) :: enthalpy_prec_bc
 
    character(len=*), parameter :: sub = 'cam_export'
    !-----------------------------------------------------------------------
@@ -475,13 +479,13 @@ subroutine cam_export(state,cam_in,cam_out,pbuf)
       !
       ! Idem for ice
       !
-      enthalpy_prec_bc(:ncol,fice_idx) = fice_tot(:ncol)-enthalpy_prec_ac(:ncol,fice_idx)
-      enthalpy_prec_bc(:ncol,fliq_idx) = fliq_tot(:ncol)-enthalpy_prec_ac(:ncol,fliq_idx)
+      enthalpy_prec_bc(:ncol,fice_idx) = fice_tot(:ncol) -enthalpy_prec_ac(:ncol,fice_idx)
+      enthalpy_prec_bc(:ncol,fliq_idx) = fliq_tot(:ncol) -enthalpy_prec_ac(:ncol,fliq_idx)
       !
       ! compute precipitation enthalpy fluxes from tphysbc
       !
       enthalpy_prec_bc(:ncol,hice_idx) =  -enthalpy_prec_bc(:ncol,fice_idx)*cpice*state%T(:ncol,pver)
-      enthalpy_prec_bc(:ncol,hliq_idx) =  -enthalpy_prec_bc(:ncol,fliq_idx)*cpice*state%T(:ncol,pver)
+      enthalpy_prec_bc(:ncol,hliq_idx) =  -enthalpy_prec_bc(:ncol,fliq_idx)*cpliq*state%T(:ncol,pver)
       call pbuf_set_field(pbuf, enthalpy_prec_bc_idx, enthalpy_prec_bc)
       !
       ! compute evaporation enthalpy flux

src/physics/cam/cam_diagnostics.F90

@@ -1,3 +1,4 @@
+#define debug_entalpy
 module cam_diagnostics
 
 !---------------------------------------------------------------------------------
@@ -394,6 +395,18 @@ subroutine diag_init_dry(pbuf2d)
     call addfld( 'CPAIRV', (/ 'lev' /), 'I', 'J/K/kg', 'Variable specific heat cap air' )
     call addfld( 'RAIRV', (/ 'lev' /), 'I', 'J/K/kg', 'Variable dry air gas constant' )
 
+#ifdef debug_entalpy
+    call addfld('enth_prec_ac_hice', horiz_only, 'I', 'W/m2', '' )
+    call addfld('enth_prec_ac_hliq', horiz_only, 'I', 'W/m2', '' )
+    call addfld('enth_prec_bc_hice', horiz_only, 'I', 'W/m2', '' )
+    call addfld('enth_prec_bc_hliq', horiz_only, 'I', 'W/m2', '' )
+    call addfld('enth_prec_ac_fice', horiz_only, 'I', 'W/m2', '' )
+    call addfld('enth_prec_ac_fliq', horiz_only, 'I', 'W/m2', '' )
+    call addfld('enth_prec_bc_fice', horiz_only, 'I', 'W/m2', '' )
+    call addfld('enth_prec_bc_fliq', horiz_only, 'I', 'W/m2', '' )
+    call addfld('enth_evap_hevap', horiz_only, 'I', 'W/m2', '' )
+#endif
+
     if (thermo_budget_history) then
        !
        ! energy diagnostics addflds for vars_stage combinations plus e_m_snapshots

src/physics/cam/check_energy.F90

@@ -1,4 +1,4 @@
-
+#define debug_entalpy
 module check_energy
 
 !---------------------------------------------------------------------------------
@@ -52,6 +52,7 @@ module check_energy
 
   public :: tot_energy_phys ! calculate and output total energy and axial angular momentum diagnostics
 
+  public :: enthalpy_adjustment
 ! Private module data
 
   logical  :: print_energy_errors = .false.
@@ -964,5 +965,73 @@ subroutine tot_energy_phys(state, outfld_name_suffix,vc)
 
   end subroutine tot_energy_phys
 
-
+  subroutine enthalpy_adjustment(ncol, lchnk, state, cam_in, cam_out, pbuf, ztodt, itim_old,&
+       qini,totliqini,toticeini,tend)
+    use camsrfexch,      only: cam_out_t, cam_in_t, get_prec_vars
+    use physics_buffer,  only: pbuf_get_index, physics_buffer_desc, pbuf_set_field, pbuf_get_field
+    use cam_abortutils,  only: endrun
+    use air_composition, only: hliq_idx, hice_idx, fliq_idx, fice_idx, num_enthalpy_vars
+    use physconst,       only: cpliq, cpice
+#ifdef debug_entalpy
+    use cam_history,     only: outfld
+#endif
+    integer,             intent(in)    :: ncol, lchnk
+    type(physics_state), intent(inout) :: state
+    type(cam_in_t),      intent(inout) :: cam_in
+    type(cam_out_t),     intent(inout) :: cam_out
+    type(physics_buffer_desc), pointer :: pbuf(:)
+    real(r8),            intent(in)    :: ztodt
+    integer,             intent(in)    :: itim_old
+
+    real(r8), dimension(pcols,pver), intent(in) :: qini, totliqini, toticeini
+    type(physics_tend )    , intent(inout) :: tend!xxx
+
+    integer:: enthalpy_prec_bc_idx, enthalpy_prec_ac_idx, enthalpy_evap_idx
+    real(r8), dimension(:,:), pointer            :: enthalpy_prec_bc
+    real(r8), dimension(pcols,num_enthalpy_vars) :: enthalpy_prec_ac
+    real(r8), dimension(:)  , pointer            :: enthalpy_evap
+    real(r8), dimension(pcols)                   :: fliq_tot, fice_tot
+
+    integer :: i
+
+    enthalpy_prec_bc_idx = pbuf_get_index('ENTHALPY_PREC_BC', errcode=i)
+    enthalpy_prec_ac_idx = pbuf_get_index('ENTHALPY_PREC_AC', errcode=i)
+    enthalpy_evap_idx    = pbuf_get_index('ENTHALPY_EVAP'   , errcode=i)
+    if (enthalpy_prec_bc_idx==0.or.enthalpy_prec_ac_idx==0.or.enthalpy_evap_idx==0) then
+      call endrun("pbufs for enthalpy flux not allocated")
+    end if
+    call pbuf_get_field(pbuf, enthalpy_prec_bc_idx, enthalpy_prec_bc)
+    call pbuf_get_field(pbuf, enthalpy_evap_idx   , enthalpy_evap   )
+    !
+    !------------------------------------------------------------------
+    !
+    ! compute precipitation fluxes and set associated physics buffers
+    !
+    !------------------------------------------------------------------
+    !
+    call get_prec_vars(ncol,pbuf,fliq=fliq_tot,fice=fice_tot)
+    !
+    ! fliq_tot holds liquid precipitation from tphysbc and tphysac; idem for ice
+    !
+    enthalpy_prec_ac(:ncol,fice_idx) = fice_tot(:ncol)-enthalpy_prec_bc(:ncol,fice_idx)
+    enthalpy_prec_ac(:ncol,fliq_idx) = fliq_tot(:ncol)-enthalpy_prec_bc(:ncol,fliq_idx)
+    !
+    ! compute precipitation enthalpy fluxes from tphysbc
+    !
+    enthalpy_prec_ac(:ncol,hice_idx) =  -enthalpy_prec_ac(:ncol,fice_idx)*cpice*state%T(:ncol,pver)
+    enthalpy_prec_ac(:ncol,hliq_idx) =  -enthalpy_prec_ac(:ncol,fliq_idx)*cpliq*state%T(:ncol,pver)
+
+#ifdef debug_entalpy
+    call outfld("enth_prec_ac_hice"  , enthalpy_prec_ac(:,hice_idx)     , pcols   ,lchnk   )
+    call outfld("enth_prec_ac_hliq"  , enthalpy_prec_ac(:,hliq_idx)     , pcols   ,lchnk   )
+    call outfld("enth_prec_bc_hice"  , enthalpy_prec_bc(:,hice_idx)     , pcols   ,lchnk   )
+    call outfld("enth_prec_bc_hliq"  , enthalpy_prec_bc(:,hliq_idx)     , pcols   ,lchnk   )
+    call outfld("enth_prec_ac_fice"  , enthalpy_prec_ac(:,fice_idx)     , pcols   ,lchnk   )
+    call outfld("enth_prec_ac_fliq"  , enthalpy_prec_ac(:,fliq_idx)     , pcols   ,lchnk   )
+    call outfld("enth_prec_bc_fice"  , enthalpy_prec_bc(:,fice_idx)     , pcols   ,lchnk   )
+    call outfld("enth_prec_bc_fliq"  , enthalpy_prec_bc(:,fliq_idx)     , pcols   ,lchnk   )
+    call outfld("enth_evap_hevap"    , enthalpy_evap   (:)              , pcols   ,lchnk   )    
+#endif
+    call pbuf_set_field(pbuf, enthalpy_prec_ac_idx, enthalpy_prec_ac)
+  end subroutine enthalpy_adjustment
 end module check_energy

src/physics/cam_dev/physpkg.F90

@@ -1372,7 +1372,7 @@ subroutine tphysac (ztodt,   cam_in,  &
     use aoa_tracers,        only: aoa_tracers_timestep_tend
     use physconst,          only: rhoh2o
     use aero_model,         only: aero_model_drydep
-    use check_energy,       only: check_energy_chng, tot_energy_phys
+    use check_energy,       only: check_energy_chng, tot_energy_phys, enthalpy_adjustment
     use check_energy,       only: check_tracers_data, check_tracers_init, check_tracers_chng
     use time_manager,       only: get_nstep
     use cam_abortutils,     only: endrun
@@ -2367,6 +2367,8 @@ subroutine tphysac (ztodt,   cam_in,  &
           call endrun("Explicit enthalpy flux functionality only supported for SE dycore")
        end if
        !Thomas: we will add call here to subroutine that does the simple dme bflx etc.
+       call enthalpy_adjustment(ncol,lchnk,state,cam_in,cam_out,pbuf,ztodt,itim_old,&
+         qini,totliqini,toticeini,tend)
     else
        !-------------- Energy budget checks vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
        ! Save total energy for global fixer in next timestep
@@ -2511,6 +2513,8 @@ subroutine tphysbc (ztodt, state,  &
     use constituents,    only: qmin
     use air_composition, only: thermodynamic_active_species_liq_num,thermodynamic_active_species_liq_idx
     use air_composition, only: thermodynamic_active_species_ice_num,thermodynamic_active_species_ice_idx
+    use air_composition, only: compute_enthalpy_flux, num_enthalpy_vars
+    use physics_buffer,  only: pbuf_set_field
     use convect_deep,    only: convect_deep_tend
     use time_manager,    only: is_first_step, get_nstep
     use convect_diagnostics,only: convect_diagnostics_calc
@@ -2604,6 +2608,8 @@ subroutine tphysbc (ztodt, state,  &
     real(r8),pointer :: prec_sed(:)     ! total precip from cloud sedimentation
     real(r8),pointer :: snow_sed(:)     ! snow from cloud ice sedimentation
 
+    real(r8) :: enthalpy_prec_ac(pcols,num_enthalpy_vars)
+
     ! energy checking variables
     real(r8) :: zero(pcols)                    ! array of zeros
     real(r8) :: zero_sc(pcols*psubcols)        ! array of zeros
@@ -2867,6 +2873,18 @@ subroutine tphysbc (ztodt, state,  &
       snow_sed = 0._r8
       prec_str = 0._r8
       snow_str = 0._r8
+      !
+      ! In first time-step tphysac variables need to be zero'd out 
+      !
+      if (compute_enthalpy_flux) then
+        ifld = pbuf_get_index('ENTHALPY_PREC_AC', errcode=i)
+        enthalpy_prec_ac = 0._r8
+        if (ifld>0) then
+          call pbuf_set_field(pbuf, ifld, enthalpy_prec_ac)
+        else
+          call endrun('tphysbc: pbuf ENTHALPY_PREC_AC not allocated')
+        end if
+      end if
 
       if (is_subcol_on()) then
         prec_str_sc = 0._r8

src/utils/air_composition.F90

@@ -25,7 +25,7 @@ module air_composition
    ! get_mbarv: molecular weight of dry air
    public :: get_mbarv
 
-   logical, public :: compute_enthalpy_flux=.false.
+   logical, public :: compute_enthalpy_flux=.true.
    !
    ! for book keeping of enthalpy variables in physics buffer
    !