1st simple version with explicit enthalpy flux representation

src/physics/cam/cam_diagnostics.F90

@@ -1,4 +1,3 @@
-#define debug_entalpy
 module cam_diagnostics
 
 !---------------------------------------------------------------------------------
@@ -395,7 +394,21 @@ 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
+    !
+    ! enthalpy variables
+    !
+    call addfld('STEND_HFIX', (/ 'lev' /), 'I', 'W/m2', 'Heating tendency from variable latent heat enthalpy fixer' )
+    call addfld('dEdt_enth_fix', horiz_only, 'I', 'W/m2', 'Column integrated dEdt from variable latent heat enthalpy fixer' )
+    call addfld('dEdt_efix_physics', horiz_only, 'I', 'W/m2', 'Column integrated physics energy fixer dEdt from enthalpy fixer' )
+    call addfld('residual', horiz_only, 'I', 'W/m2', '' )
+    call addfld('enth_fix_fct_bc', (/ 'lev' /), 'I', 'W/m2', '' )
+    call addfld('enth_fix_fct_ac', (/ 'lev' /), 'I', 'W/m2', '' )
+    call addfld('enth_fix_fct_bc_tot', horiz_only, 'I', 'W/m2', '' )
+    call addfld('enth_fix_fct_ac_tot', horiz_only, 'I', 'W/m2', '' )
+
+    call addfld('enthalpy_heating_fix_bc', horiz_only, 'I', 'W/m2', '' )
+    call addfld('enthalpy_heating_fix_ac', horiz_only, 'I', 'W/m2', '' )
+
     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', '' )
@@ -405,7 +418,21 @@ subroutine diag_init_dry(pbuf2d)
     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
+
+    call addfld('te_tnd', horiz_only, 'I', 'W/m2', 'Total column integrated energy tendency from CAM physics' )
+    call addfld('te_lat', horiz_only, 'I', 'W/m2', 'Total column integrated constant latent heat tendency from CAM physics')
+    call addfld('cnst_lat_heat_srf'       , horiz_only, 'I', 'W/m2', '' )!xxx diags will remove
+    call addfld('cpice_srf'    , horiz_only, 'I', 'W/m2', '' )!xxx diags will remove
+    call addfld('ls_srf'       , horiz_only, 'I', 'W/m2', '' )!xxx diags will remove
+    call addfld('lf_srf'       , horiz_only, 'I', 'W/m2', '' )!xxx diags will remove
+    call addfld('dEdt_dme'     , horiz_only, 'I', 'W/m2', 'Column integrated dEdt from water update')
+    call addfld('dEdt_physics' , horiz_only, 'I', 'W/m2', '' )!xxx diags will remove
+    call addfld('dEdt_cpdycore' , horiz_only, 'I', 'W/m2', 'Column integrated dEdt from updating cp')
+    !
+    !fincl1  = 'enth_prec_ac_hice','enth_prec_ac_hliq','enth_prec_bc_hice','enth_prec_bc_hliq','enth_prec_ac_fice',!xxx
+    !          'enth_prec_ac_fliq','enth_prec_bc_fice','enth_prec_bc_fliq','enth_evap_hevap',!xxx
+    !          'heating_cpice','te_tnd','te_lat','cnst_lat_heat_srf','ls_srf','lf_srf','dEdt_dme','dEdt_physics','dEdt_cpdycore'!xxx
+    !
 
     if (thermo_budget_history) then
        !

src/physics/cam/clubb_intr.F90

@@ -401,8 +401,8 @@ module clubb_intr
     ice_supersat_idx, & ! ice cloud fraction for SILHS
     rcm_idx, &          ! Cloud water mixing ratio for SILHS
     ztodt_idx,&         ! physics timestep for SILHS
-    clubbtop_idx        ! level index for CLUBB top
-
+    clubbtop_idx,     & ! level index for CLUBB top
+    rcmtend_clubb_idx
   ! Indices for microphysical covariance tendencies
   integer :: &
     rtp2_mc_zt_idx,   &
@@ -587,6 +587,7 @@ subroutine clubb_register_cam( )
     call pbuf_add_field('pdf_zm_var_w_2', 'global', dtype_r8, (/pcols,pverp,dyn_time_lvls/), pdf_zm_varnce_w_2_idx)
     call pbuf_add_field('pdf_zm_mixt_frac',  'global', dtype_r8, (/pcols,pverp,dyn_time_lvls/), pdf_zm_mixt_frac_idx)
 
+    call pbuf_add_field('rcmtend_clubb', 'global', dtype_r8, (/pcols, pver/), rcmtend_clubb_idx)
 #endif
 
   end subroutine clubb_register_cam
@@ -3907,8 +3908,11 @@ subroutine clubb_tend_cam( state,   ptend_all,   pbuf,     hdtime, &
     call outfld( 'RVMTEND_CLUBB', temp2d, pcols, lchnk)
 
     temp2d(:ncol,:pver) = ptend_loc%q(:ncol,:pver,ixcldliq)*state1%pdeldry(:ncol,:pver)/state1%pdel(:ncol,:pver)
+    call pbuf_set_field(pbuf, rcmtend_clubb_idx,temp2d)
     call outfld( 'RCMTEND_CLUBB', temp2d, pcols, lchnk)
 
+
+
     temp2d(:ncol,:pver) = ptend_loc%q(:ncol,:pver,ixcldice)*state1%pdeldry(:ncol,:pver)/state1%pdel(:ncol,:pver)
     call outfld( 'RIMTEND_CLUBB', temp2d, pcols, lchnk)
 

src/physics/cam/physics_types.F90

@@ -126,6 +126,7 @@ module physics_types
      real(r8), dimension(:),  allocatable        :: flx_net
      real(r8), dimension(:),  allocatable        :: &
           te_tnd,  &! cumulative boundary flux of total energy
+          te_lat,  &! cumulative latent heat fluxes (assuming constant latent heats)
           tw_tnd    ! cumulative boundary flux of total water
   end type physics_tend
 
@@ -1432,6 +1433,7 @@ subroutine physics_tend_init(tend)
     tend%dvdt    = 0._r8
     tend%flx_net = 0._r8
     tend%te_tnd  = 0._r8
+    tend%te_lat  = 0._r8
     tend%tw_tnd  = 0._r8
 
 end subroutine physics_tend_init
@@ -1847,6 +1849,9 @@ subroutine physics_tend_alloc(tend,psetcols)
   allocate(tend%te_tnd(psetcols), stat=ierr)
   if ( ierr /= 0 ) call endrun('physics_tend_alloc error: allocation error for tend%te_tnd')
 
+  allocate(tend%te_lat(psetcols), stat=ierr)
+  if ( ierr /= 0 ) call endrun('physics_tend_alloc error: allocation error for tend%te_lat')
+
   allocate(tend%tw_tnd(psetcols), stat=ierr)
   if ( ierr /= 0 ) call endrun('physics_tend_alloc error: allocation error for tend%tw_tnd')
 
@@ -1855,6 +1860,7 @@ subroutine physics_tend_alloc(tend,psetcols)
   tend%dvdt(:,:) = inf
   tend%flx_net(:) = inf
   tend%te_tnd(:) = inf
+  tend%te_lat(:) = inf
   tend%tw_tnd(:) = inf
 
 end subroutine physics_tend_alloc
@@ -1883,6 +1889,9 @@ subroutine physics_tend_dealloc(tend)
   deallocate(tend%te_tnd, stat=ierr)
   if ( ierr /= 0 ) call endrun('physics_tend_dealloc error: deallocation error for tend%te_tnd')
 
+  deallocate(tend%te_lat, stat=ierr)
+  if ( ierr /= 0 ) call endrun('physics_tend_dealloc error: deallocation error for tend%te_lat')
+
   deallocate(tend%tw_tnd, stat=ierr)
   if ( ierr /= 0 ) call endrun('physics_tend_dealloc error: deallocation error for tend%tw_tnd')
 end subroutine physics_tend_dealloc

src/physics/cam/zm_conv_intr.F90

@@ -41,6 +41,7 @@ module zm_conv_intr
    public zmconv_ke, zmconv_ke_lnd,  zmconv_org  ! needed by convect_shallow
 
    integer ::& ! indices for fields in the physics buffer
+      zmdt_idx,       &
       zm_mu_idx,      &
       zm_eu_idx,      &
       zm_du_idx,      &
@@ -62,7 +63,8 @@ module zm_conv_intr
       dnifzm_idx,    &     ! detrained convective cloud ice num concen.
       prec_dp_idx,   &
       snow_dp_idx,   &
-      mconzm_idx           ! convective mass flux
+      mconzm_idx,    &     ! convective mass flux
+      rprd_idx
 
    real(r8), parameter :: unset_r8 = huge(1.0_r8)
    real(r8) :: zmconv_c0_lnd = unset_r8
@@ -109,6 +111,7 @@ subroutine zm_conv_register
 
   integer idx
 
+   call pbuf_add_field('ZMDT' , 'physpkg', dtype_r8, (/pcols,pver/), zmdt_idx)
    call pbuf_add_field('ZM_MU', 'physpkg', dtype_r8, (/pcols,pver/), zm_mu_idx)
    call pbuf_add_field('ZM_EU', 'physpkg', dtype_r8, (/pcols,pver/), zm_eu_idx)
    call pbuf_add_field('ZM_DU', 'physpkg', dtype_r8, (/pcols,pver/), zm_du_idx)
@@ -160,6 +163,7 @@ subroutine zm_conv_register
       call cnst_add('ZM_ORG',0._r8,0._r8,0._r8,ixorg,longname='organization parameter')
    endif
 
+   call pbuf_add_field('rprd','physpkg',dtype_r8,(/pcols,pver/),rprd_idx)
 end subroutine zm_conv_register
 
 !=========================================================================================
@@ -382,6 +386,7 @@ subroutine zm_conv_tend(pblh    ,mcon    ,cme     , &
 
    use time_manager,  only: get_nstep, is_first_step
    use physics_buffer, only : pbuf_get_field, physics_buffer_desc, pbuf_old_tim_idx
+   use physics_buffer, only : pbuf_set_field
    use constituents,  only: pcnst, cnst_get_ind, cnst_is_convtran1
    use check_energy,  only: check_energy_chng
    use physconst,     only: gravit, latice, latvap, tmelt, cpwv, cpliq, rh2o
@@ -647,6 +652,9 @@ subroutine zm_conv_tend(pblh    ,mcon    ,cme     , &
 
    ftem(:ncol,:pver) = ptend_loc%s(:ncol,:pver)/cpair
    call outfld('ZMDT    ',ftem           ,pcols   ,lchnk   )
+
+   call pbuf_set_field(pbuf, zmdt_idx, ptend_loc%s(:,:pver))
+
    call outfld('ZMDQ    ',ptend_loc%q(1,1,1) ,pcols   ,lchnk   )
    call t_stopf ('zm_convr_run')
 
@@ -736,6 +744,7 @@ subroutine zm_conv_tend(pblh    ,mcon    ,cme     , &
    call outfld('CMFMC_DP   ',mcon ,  pcols   ,lchnk   )
    call outfld('PRECCDZM   ',prec,  pcols   ,lchnk   )
 
+   call pbuf_set_field(pbuf, rprd_idx, rprd)
 
    call t_stopf ('zm_conv_evap_run')
 

src/physics/cam_dev/micro_pumas_cam.F90

@@ -171,7 +171,9 @@ module micro_pumas_cam
    ls_flxsnw_idx,      &
    relvar_idx,         &
    cmeliq_idx,         &
-   accre_enhan_idx
+   accre_enhan_idx,    &
+   mpdt_idx,           &
+   mpdice_idx
 
 ! Fields for UNICON
 integer :: &
@@ -628,6 +630,9 @@ subroutine micro_pumas_cam_register
               longname='Grid box averaged graupel/hail number', is_convtran1=.true.)
    end if
 
+   call pbuf_add_field('MPDT',     'physpkg', dtype_r8, (/pcols,pver/), mpdt_idx)
+   call pbuf_add_field('MPDICE',   'physpkg', dtype_r8, (/pcols,pver/), mpdice_idx)
+
    ! Request physics buffer space for fields that persist across timesteps.
 
    call pbuf_add_field('CLDO','global',dtype_r8,(/pcols,pver,dyn_time_lvls/), cldo_idx)
@@ -2429,7 +2434,8 @@ subroutine micro_pumas_cam_tend(state, ptend, dtime, pbuf)
 
       call physics_ptend_init(ptend_loc, psetcols, "micro_pumas", &
                               ls=.true., lq=lq)
-
+      call pbuf_set_field(pbuf, mpdt_idx, tlat)
+      call pbuf_set_field(pbuf, mpdice_idx, qiten)
       ! Set local tendency.
       ptend_loc%s(:ncol,top_lev:) = tlat(:ncol,top_lev:)
       ptend_loc%q(:ncol,top_lev:,ixq) = qvlat(:ncol,top_lev:)

src/physics/cam_dev/physpkg.F90

@@ -2366,9 +2366,8 @@ subroutine tphysac (ztodt,   cam_in,  &
        if (.not.dycore_is('SE')) then
           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)
+       call enthalpy_adjustment(ncol,lchnk,state,cam_in,pbuf,ztodt,itim_old,&
+            qini,totliqini,toticeini,tend)
     else
        !-------------- Energy budget checks vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
        ! Save total energy for global fixer in next timestep
@@ -2466,7 +2465,6 @@ subroutine tphysac (ztodt,   cam_in,  &
           call endrun ('TPHYSAC error: in aquaplanet mode, but grid contains non-ocean point')
        endif
     endif
-
     call diag_phys_tend_writeout (state, pbuf,  tend, ztodt, qini, cldliqini, cldiceini)
 
     call clybry_fam_set( ncol, lchnk, map2chm, state%q, pbuf )
@@ -2513,7 +2511,7 @@ 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 air_composition, only: compute_enthalpy_flux, num_enthalpy_vars, cp_or_cv_dycore
     use physics_buffer,  only: pbuf_set_field
     use convect_deep,    only: convect_deep_tend
     use time_manager,    only: is_first_step, get_nstep
@@ -2532,7 +2530,9 @@ subroutine tphysbc (ztodt, state,  &
     use cam_snapshot,    only: cam_snapshot_all_outfld_tphysbc
     use cam_snapshot_common, only: cam_snapshot_ptend_outfld
     use dyn_tests_utils, only: vc_dycore
-
+    use air_composition, only: te_init,cpairv,compute_enthalpy_flux!xxx to be removed
+    use dyn_tests_utils, only: vc_dycore!xxx to be removed
+    use cam_thermo,      only: get_hydrostatic_energy!xxx to be removed
     ! Arguments
 
     real(r8), intent(in) :: ztodt                          ! 2 delta t (model time increment)
@@ -2608,8 +2608,6 @@ 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
@@ -2620,9 +2618,7 @@ subroutine tphysbc (ztodt, state,  &
     real(r8) :: flx_heat(pcols)
     type(check_tracers_data):: tracerint             ! energy integrals and cummulative boundary fluxes
     real(r8) :: zero_tracers(pcols,pcnst)
-
     logical   :: lq(pcnst)
-
     !-----------------------------------------------------------------------
 
     call t_startf('bc_init')
@@ -2691,9 +2687,39 @@ subroutine tphysbc (ztodt, state,  &
 
     call t_stopf('bc_init')
 
+    call cnst_get_ind('Q', ixq)
+    call cnst_get_ind('CLDLIQ', ixcldliq)
+    call cnst_get_ind('CLDICE', ixcldice)
+    qini     (:ncol,:pver) = state%q(:ncol,:pver,     ixq)
+    cldliqini(:ncol,:pver) = state%q(:ncol,:pver,ixcldliq)
+    cldiceini(:ncol,:pver) = state%q(:ncol,:pver,ixcldice)
+
+    totliqini(:ncol,:pver) = 0.0_r8
+    do m_cnst=1,thermodynamic_active_species_liq_num
+      m = thermodynamic_active_species_liq_idx(m_cnst)
+      totliqini(:ncol,:pver) = totliqini(:ncol,:pver)+state%q(:ncol,:pver,m)
+    end do
+    toticeini(:ncol,:pver) = 0.0_r8
+    do m_cnst=1,thermodynamic_active_species_ice_num
+      m = thermodynamic_active_species_ice_idx(m_cnst)
+      toticeini(:ncol,:pver) = toticeini(:ncol,:pver)+state%q(:ncol,:pver,m)
+    end do
+
+    !
+    ! compute energy variables for state at the beginning of physics - xxx to be remove
+    !
+    if (compute_enthalpy_flux) then
+      call get_hydrostatic_energy(state%q(1:ncol,1:pver,1:pcnst),.true.,          &
+           state%pdel(1:ncol,1:pver), cp_or_cv_dycore(:ncol,:,lchnk),             &
+           state%u(1:ncol,1:pver), state%v(1:ncol,1:pver), state%T(1:ncol,1:pver),&
+           vc_dycore, ptop=state%pintdry(1:ncol,1), phis = state%phis(1:ncol),    &
+           te = te_init(:ncol,1,lchnk), se=te_init(:ncol,2,lchnk), po=te_init(:ncol,3,lchnk), ke=te_init(:ncol,4,lchnk))
+    endif
+
     !===================================================
     ! Global mean total energy fixer
     !===================================================
+
     call t_startf('energy_fixer')
 
     call tot_energy_phys(state, 'phBF')
@@ -2711,25 +2737,6 @@ subroutine tphysbc (ztodt, state,  &
     ! Save state for convective tendency calculations.
     call diag_conv_tend_ini(state, pbuf)
 
-    call cnst_get_ind('Q', ixq)
-    call cnst_get_ind('CLDLIQ', ixcldliq)
-    call cnst_get_ind('CLDICE', ixcldice)
-    qini     (:ncol,:pver) = state%q(:ncol,:pver,     ixq)
-    cldliqini(:ncol,:pver) = state%q(:ncol,:pver,ixcldliq)
-    cldiceini(:ncol,:pver) = state%q(:ncol,:pver,ixcldice)
-
-    totliqini(:ncol,:pver) = 0.0_r8
-    do m_cnst=1,thermodynamic_active_species_liq_num
-      m = thermodynamic_active_species_liq_idx(m_cnst)
-      totliqini(:ncol,:pver) = totliqini(:ncol,:pver)+state%q(:ncol,:pver,m)
-    end do
-    toticeini(:ncol,:pver) = 0.0_r8
-    do m_cnst=1,thermodynamic_active_species_ice_num
-      m = thermodynamic_active_species_ice_idx(m_cnst)
-      toticeini(:ncol,:pver) = toticeini(:ncol,:pver)+state%q(:ncol,:pver,m)
-    end do
-
-
     call outfld('TEOUT', teout       , pcols, lchnk   )
     call outfld('TEINP', state%te_ini(:,dyn_te_idx), pcols, lchnk   )
     call outfld('TEFIX', state%te_cur(:,dyn_te_idx), pcols, lchnk   )
@@ -2874,16 +2881,11 @@ subroutine tphysbc (ztodt, state,  &
       prec_str = 0._r8
       snow_str = 0._r8
       !
-      ! In first time-step tphysac variables need to be zero'd out 
+      ! 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
+        if (ifld>0) call pbuf_set_field(pbuf, ifld, 0._r8)
       end if
 
       if (is_subcol_on()) then

src/utils/air_composition.F90

@@ -118,6 +118,7 @@ module air_composition
    ! cp_or_cv_dycore:  enthalpy or internal energy scaling factor for 
    !                   energy consistency
    real(r8), public, protected, allocatable :: cp_or_cv_dycore(:,:,:)
+   real(r8), public, allocatable :: te_init(:,:,:)!xxx to be removed
    !
    ! Interfaces for public routines
    interface get_cp_dry
@@ -354,7 +355,7 @@ subroutine air_composition_init()
       if (ierr /= 0) then
          call endrun(errstr//"cp_or_cv_dycore")
       end if
-
+      allocate(te_init(pcols,4,begchunk:endchunk), stat=ierr)!xxx to be removed
       thermodynamic_active_species_idx        = -HUGE(1)
       thermodynamic_active_species_idx_dycore = -HUGE(1)
       thermodynamic_active_species_cp         = 0.0_r8