Merge branch 'kaizhangpnl/atm/check_energy_bugfix_AND_ieflx' (PR #1303)

This PR implements the IEFLX fixer
(https://acme-climate.atlassian.net/wiki/display/ATM/Simulation+with+the+IEFLX+fixer)

These modifications add to those in PR #1301.

The following variable will be added to the master list if ieflx_opt > 0
SHFLXFIX : SHFLX plus global mean IEFLX (sensible heat flux used in the atmosphere model)

To switch on the IEFLX fixer, set
ieflx_opt = 2

By default, the IEFLX fixer is switched off.

[BFB]

* kaizhangpnl/atm/check_energy_bugfix_AND_ieflx:
  Implement the IEFLX fixer

components/cam/bld/namelist_files/namelist_definition.xml

@@ -1116,6 +1116,18 @@ Whether or not to printout fixer diagnostic message.
 Default: set by build-namelist.
 </entry>
 
+<entry id="ieflx_opt" type="integer"  category="fixers"
+       group="phys_ctl_nl" valid_values="">
+IEFLX fixer options
+Default: 0
+</entry>
+
+<entry id="l_ieflx_fix" type="logical" category="fixers"
+       group="phys_ctl_nl" valid_values="" >
+Whether or not to fix ieflx
+Default: .false.
+</entry>
+
 <!-- Gravity Wave Drag -->
 
 <entry id="use_gw_oro" type="logical" category="gw_drag"

components/cam/src/physics/cam/cam_diagnostics.F90

@@ -1500,6 +1500,7 @@ subroutine diag_surf (cam_in, cam_out, ps, trefmxav, trefmnav )
    use time_manager,     only: is_end_curr_day
    use co2_cycle,        only: c_i, co2_transport
    use constituents,     only: sflxnam
+   use phys_control,     only: ieflx_opt, l_ieflx_fix
 
 !-----------------------------------------------------------------------
 !
@@ -1526,6 +1527,8 @@ subroutine diag_surf (cam_in, cam_out, ps, trefmxav, trefmnav )
     lchnk = cam_in%lchnk
     ncol  = cam_in%ncol
 
+    !! if l_ieflx_fix = true, SHFLX will be output in physpkg.F90 
+
     call outfld('SHFLX',    cam_in%shf,       pcols, lchnk)
     call outfld('LHFLX',    cam_in%lhf,       pcols, lchnk)
     call outfld('QFLX',     cam_in%cflx(1,1), pcols, lchnk)

components/cam/src/physics/cam/check_energy.F90

@@ -34,6 +34,7 @@ module check_energy
   use time_manager,    only: is_first_step
   use cam_logfile,     only: iulog
   use cam_abortutils,  only: endrun 
+  use phys_control,    only: ieflx_opt !!l_ieflx_fix
 
   implicit none
   private
@@ -59,6 +60,9 @@ module check_energy
   public :: check_prect             ! output prect at certain locations for water conservation check  
   public :: check_water             ! output water path at certain locations for water conservation check  
 
+  public :: ieflx_gmean             ! calculate global mean of ieflx 
+  public :: check_ieflx_fix         ! add ieflx to sensible heat flux 
+
 
 ! Private module data
 
@@ -70,9 +74,11 @@ module check_energy
   real(r8) :: psurf_glob           ! global mean surface pressure
   real(r8) :: ptopb_glob           ! global mean top boundary pressure
   real(r8) :: heat_glob            ! global mean heating rate
+  real(r8) :: ieflx_glob           ! global mean implied internal energy flux 
 
 ! Physics buffer indices
 
+  integer  :: ieflx_idx  = 0       ! teout index in physics buffer 
   integer  :: teout_idx  = 0       ! teout index in physics buffer 
   integer  :: dtcore_idx = 0       ! dtcore index in physics buffer 
 
@@ -137,11 +143,20 @@ subroutine check_energy_register()
 
     call pbuf_add_field('TEOUT', 'global',dtype_r8 , (/pcols,dyn_time_lvls/),      teout_idx)
     call pbuf_add_field('DTCORE','global',dtype_r8,  (/pcols,pver,dyn_time_lvls/),dtcore_idx)
+
     if(is_subcol_on()) then
       call pbuf_register_subcol('TEOUT', 'phys_register', teout_idx)
+      call pbuf_register_subcol('IEFLX', 'phys_register', ieflx_idx)
       call pbuf_register_subcol('DTCORE', 'phys_register', dtcore_idx)
     end if
 
+    if(ieflx_opt>0) then 
+       call pbuf_add_field('IEFLX', 'global',dtype_r8 , (/pcols,dyn_time_lvls/),      ieflx_idx)
+       if(is_subcol_on()) then
+         call pbuf_register_subcol('IEFLX', 'phys_register', ieflx_idx)
+       end if
+    end if 
+
   end subroutine check_energy_register
 
 !===============================================================================
@@ -220,6 +235,12 @@ subroutine check_energy_init()
        call add_default ('DTCORE   '  , history_budget_histfile_num, ' ')
     end if
 
+    if(ieflx_opt>0) then 
+       call addfld('IEFLX',    horiz_only, 'A', 'W/m2', 'Implied internal energy flux')
+       call addfld('SHFLXFIX', horiz_only, 'A', 'W/m2', 'SHFLX after adding IEFLX')
+       call add_default ('SHFLXFIX', 1, ' ') 
+    end if 
+
   end subroutine check_energy_init
 
 !===============================================================================
@@ -244,6 +265,8 @@ subroutine check_energy_timestep_init(state, tend, pbuf, col_type)
     real(r8) :: wl(state%ncol)                     ! vertical integral of water (liquid)
     real(r8) :: wi(state%ncol)                     ! vertical integral of water (ice)
 
+    real(r8) :: ieflx(pcols)                     ! vertical integral of kinetic energy
+
     real(r8),allocatable :: cpairv_loc(:,:,:)
 
     integer lchnk                                  ! chunk identifier
@@ -285,6 +308,9 @@ subroutine check_energy_timestep_init(state, tend, pbuf, col_type)
     wi = 0._r8
     wr = 0._r8
     ws = 0._r8
+
+    ieflx(1:ncol) = 0._r8
+
     do k = 1, pver
        do i = 1, ncol
           ke(i) = ke(i) + 0.5_r8*(state%u(i,k)**2 + state%v(i,k)**2)*state%pdel(i,k)/gravit
@@ -336,6 +362,9 @@ subroutine check_energy_timestep_init(state, tend, pbuf, col_type)
 ! initialize physics buffer
     if (is_first_step()) then
        call pbuf_set_field(pbuf, teout_idx, state%te_ini, col_type=col_type)
+       if(ieflx_opt>0) then 
+          call pbuf_set_field(pbuf, ieflx_idx, ieflx,        col_type=col_type)
+       end if
     end if
 
     deallocate(cpairv_loc)
@@ -434,6 +463,7 @@ subroutine check_energy_chng(state, tend, name, nstep, ztodt,        &
     wi = 0._r8
     wr = 0._r8
     ws = 0._r8
+
     do k = 1, pver
        do i = 1, ncol
           ke(i) = ke(i) + 0.5_r8*(state%u(i,k)**2 + state%v(i,k)**2)*state%pdel(i,k)/gravit
@@ -615,6 +645,123 @@ end subroutine check_energy_gmean
 
 !===============================================================================
 
+subroutine ieflx_gmean(state, tend, pbuf2d, cam_in, cam_out, nstep)
+
+!!...................................................................
+!! Calculate global mean of the implied internal energy flux 
+!! 
+!! Author: Kai Zhang 
+!!...................................................................
+
+    use camsrfexch,       only: cam_out_t, cam_in_t
+    use physics_buffer,   only: physics_buffer_desc, pbuf_get_field, pbuf_get_chunk, pbuf_set_field 
+    use cam_history,      only: outfld
+    use phys_control,     only: ieflx_opt
+
+    ! Compute global mean qflx
+
+    integer , intent(in) :: nstep        ! current timestep number
+    type(physics_state), intent(in   ), dimension(begchunk:endchunk) :: state
+    type(physics_tend ), intent(in   ), dimension(begchunk:endchunk) :: tend
+    type(cam_in_t),                     dimension(begchunk:endchunk) :: cam_in
+    type(cam_out_t),                    dimension(begchunk:endchunk) :: cam_out
+    type(physics_buffer_desc),          pointer :: pbuf2d(:,:)
+
+    real(r8), parameter :: cpsw = 3.996e3  !! cp of sea water used in MPAS [J/kg/K]
+    real(r8), parameter :: rhow = 1.e3     !! density of water [kg/m3]
+
+    integer :: ncol                      ! number of active columns
+    integer :: lchnk                     ! chunk index
+
+    real(r8), pointer :: ieflx(:) 
+
+    real(r8) :: qflx(pcols,begchunk:endchunk) !qflx [kg/m2/s]
+    real(r8) :: rain(pcols,begchunk:endchunk) !rain [m/s] 
+    real(r8) :: snow(pcols,begchunk:endchunk) !snow [m/s] 
+    real(r8) :: ienet(pcols,begchunk:endchunk) !ieflx net [W/m2] or [J/m2/s]
+
+!- 
+    ieflx_glob = 0._r8
+
+    qflx = 0._r8 
+    rain = 0._r8 
+    snow = 0._r8 
+    ienet = 0._r8 
+
+!DIR$ CONCURRENT
+    do lchnk = begchunk, endchunk
+
+       ncol = state(lchnk)%ncol
+       qflx(:ncol,lchnk) = cam_in(lchnk)%cflx(:ncol,1)
+       snow(:ncol,lchnk) = cam_out(lchnk)%precsc(:ncol) + cam_out(lchnk)%precsl(:ncol)
+       rain(:ncol,lchnk) = cam_out(lchnk)%precc(:ncol)  + cam_out(lchnk)%precl(:ncol) - snow(:ncol,lchnk) 
+
+       !! the calculation below (rhow*) converts the unit of precipitation from m/s to kg/m2/s 
+
+       select case (ieflx_opt) 
+
+       case(1) 
+          ienet(:ncol,lchnk) = cpsw * qflx(:ncol,lchnk) * cam_in(lchnk)%ts(:ncol) - & 
+                               cpsw * rhow * ( rain(:ncol,lchnk) + snow(:ncol,lchnk) ) * cam_out(lchnk)%tbot(:ncol)
+       case(2) 
+          ienet(:ncol,lchnk) = cpsw * qflx(:ncol,lchnk) * cam_in(lchnk)%ts(:ncol) - & 
+                               cpsw * rhow * ( rain(:ncol,lchnk) + snow(:ncol,lchnk) ) * cam_in(lchnk)%ts(:ncol)
+       case default 
+          call endrun('*** incorrect ieflx_opt ***')
+       end select 
+
+       !! put it to pbuf for more comprehensive treatment in the future 
+
+       call pbuf_get_field(pbuf_get_chunk(pbuf2d,lchnk),ieflx_idx, ieflx)
+
+       ieflx(:ncol) = ienet(:ncol,lchnk) 
+
+       call outfld('IEFLX', ieflx(:ncol), pcols, lchnk)
+
+    end do
+
+    call gmean(ienet, ieflx_glob)
+
+!!!    if (begchunk .le. endchunk) then
+!!!       if (masterproc) then
+!!!          write(iulog,'(1x,a12,1x,i8,4(1x,e25.17))') "nstep, ieflx, ieup, iedn ", nstep, ieflx_glob, ieup_glob, iedn_glob 
+!!!       end if
+!!!    end if
+
+  end subroutine ieflx_gmean
+
+
+!===============================================================================
+  subroutine check_ieflx_fix(lchnk, ncol, nstep, shflx)
+
+!!
+!! Add implied internal energy flux to the sensible heat flux 
+!!
+!! Called by typhsac 
+!! 
+
+    use cam_history,       only: outfld
+
+    integer, intent(in   ) :: nstep          ! time step number
+    integer, intent(in   ) :: lchnk  
+    integer, intent(in   ) :: ncol
+    real(r8),intent(inout) :: shflx(pcols) 
+
+    integer :: i
+
+    if(nstep>1) then 
+       do i = 1, ncol
+          shflx(i) = shflx(i) + ieflx_glob 
+       end do
+    end if 
+
+    call outfld('SHFLXFIX', shflx, pcols, lchnk)
+
+    return
+  end subroutine check_ieflx_fix
+
+!===============================================================================
+
 subroutine qflx_gmean(state, tend, cam_in, dtime, nstep)
 
 !!...................................................................

components/cam/src/physics/cam/phys_control.F90

@@ -105,6 +105,9 @@ module phys_control
 logical, public, protected :: use_qqflx_fixer      = .false.     ! switch on water vapor fixer to compensate changes in qflx
 logical, public, protected :: print_fixer_message  = .false.     ! switch on error message printout in log file
 
+integer, public, protected :: ieflx_opt = 0
+logical, public, protected :: l_ieflx_fix = .false.
+
 ! Macro/micro-physics co-substeps
 integer           :: cld_macmic_num_steps = 1
 
@@ -171,6 +174,8 @@ subroutine phys_ctl_readnl(nlfile)
       history_eddy, history_budget,  history_budget_histfile_num, history_waccm, &
       conv_water_in_rad, history_clubb, do_clubb_sgs, do_tms, state_debug_checks, &
       use_mass_borrower, & 
+      l_ieflx_fix, &
+      ieflx_opt, & 
       use_qqflx_fixer, & 
       print_fixer_message, & 
       use_hetfrz_classnuc, use_gw_oro, use_gw_front, use_gw_convect, &
@@ -224,6 +229,8 @@ subroutine phys_ctl_readnl(nlfile)
    call mpibcast(conv_water_in_rad,               1 , mpiint,  0, mpicom)
    call mpibcast(do_tms,                          1 , mpilog,  0, mpicom)
    call mpibcast(use_mass_borrower,               1 , mpilog,  0, mpicom)
+   call mpibcast(l_ieflx_fix,             1 , mpilog,  0, mpicom)
+   call mpibcast(ieflx_opt,               1 , mpiint,  0, mpicom)
    call mpibcast(use_qqflx_fixer,                 1 , mpilog,  0, mpicom)
    call mpibcast(print_fixer_message,             1 , mpilog,  0, mpicom)
    call mpibcast(micro_do_icesupersat,            1 , mpilog,  0, mpicom)
@@ -390,9 +397,10 @@ subroutine phys_getopts(deep_scheme_out, shallow_scheme_out, eddy_scheme_out, mi
                         radiation_scheme_out, use_subcol_microp_out, atm_dep_flux_out, &
                         history_amwg_out, history_vdiag_out, history_aerosol_out, history_aero_optics_out, history_eddy_out, &
                         history_budget_out, history_budget_histfile_num_out, history_waccm_out, &
-                        history_clubb_out, conv_water_in_rad_out, cam_chempkg_out, prog_modal_aero_out, macrop_scheme_out, &
+                        history_clubb_out, ieflx_opt_out, conv_water_in_rad_out, cam_chempkg_out, prog_modal_aero_out, macrop_scheme_out, &
                         do_clubb_sgs_out, do_tms_out, state_debug_checks_out, &
                         use_mass_borrower_out, & 
+                        l_ieflx_fix_out, & 
                         use_qqflx_fixer_out, & 
                         print_fixer_message_out, & 
                         cld_macmic_num_steps_out, micro_do_icesupersat_out, &
@@ -431,11 +439,13 @@ subroutine phys_getopts(deep_scheme_out, shallow_scheme_out, eddy_scheme_out, mi
    logical,           intent(out), optional :: history_clubb_out
    logical,           intent(out), optional :: do_clubb_sgs_out
    logical,           intent(out), optional :: micro_do_icesupersat_out
+   integer,           intent(out), optional :: ieflx_opt_out
    integer,           intent(out), optional :: conv_water_in_rad_out
    character(len=32), intent(out), optional :: cam_chempkg_out
    logical,           intent(out), optional :: prog_modal_aero_out
    logical,           intent(out), optional :: do_tms_out
    logical,           intent(out), optional :: use_mass_borrower_out
+   logical,           intent(out), optional :: l_ieflx_fix_out
    logical,           intent(out), optional :: use_qqflx_fixer_out
    logical,           intent(out), optional :: print_fixer_message_out
    logical,           intent(out), optional :: state_debug_checks_out
@@ -492,10 +502,12 @@ subroutine phys_getopts(deep_scheme_out, shallow_scheme_out, eddy_scheme_out, mi
    if ( present(do_clubb_sgs_out        ) ) do_clubb_sgs_out         = do_clubb_sgs
    if ( present(micro_do_icesupersat_out )) micro_do_icesupersat_out = micro_do_icesupersat
    if ( present(conv_water_in_rad_out   ) ) conv_water_in_rad_out    = conv_water_in_rad
+   if ( present(ieflx_opt_out   ) ) ieflx_opt_out    = ieflx_opt 
    if ( present(cam_chempkg_out         ) ) cam_chempkg_out          = cam_chempkg
    if ( present(prog_modal_aero_out     ) ) prog_modal_aero_out      = prog_modal_aero
    if ( present(do_tms_out              ) ) do_tms_out               = do_tms
    if ( present(use_mass_borrower_out   ) ) use_mass_borrower_out    = use_mass_borrower
+   if ( present(l_ieflx_fix_out   ) ) l_ieflx_fix_out    = l_ieflx_fix
    if ( present(use_qqflx_fixer_out     ) ) use_qqflx_fixer_out      = use_qqflx_fixer
    if ( present(print_fixer_message_out ) ) print_fixer_message_out  = print_fixer_message
    if ( present(state_debug_checks_out  ) ) state_debug_checks_out   = state_debug_checks

components/cam/src/physics/cam/physpkg.F90

@@ -909,6 +909,7 @@ subroutine phys_run1(phys_state, ztodt, phys_tend, pbuf2d,  cam_in, cam_out)
 #if ( defined OFFLINE_DYN )
      use metdata,       only: get_met_srf1
 #endif
+
     !
     ! Input arguments
     !
@@ -1115,6 +1116,7 @@ subroutine phys_run2(phys_state, ztodt, phys_tend, pbuf2d,  cam_out, &
     ! Purpose: 
     ! Second part of atmospheric physics package after updating of surface models
     ! 
+    ! Modified by Kai Zhang 2017-03: add IEFLX fixer treatment 
     !-----------------------------------------------------------------------
     use physics_buffer,         only: physics_buffer_desc, pbuf_get_chunk, pbuf_deallocate, pbuf_update_tim_idx
     use mo_lightning,   only: lightning_no_prod
@@ -1127,6 +1129,10 @@ subroutine phys_run2(phys_state, ztodt, phys_tend, pbuf2d,  cam_out, &
 #if ( defined OFFLINE_DYN )
     use metdata,        only: get_met_srf2
 #endif
+    use time_manager,   only: get_nstep
+    use check_energy,   only: ieflx_gmean 
+    use check_energy,   only: check_ieflx_fix 
+    use phys_control,   only: ieflx_opt !!l_ieflx_fix
     !
     ! Input arguments
     !
@@ -1146,6 +1152,7 @@ subroutine phys_run2(phys_state, ztodt, phys_tend, pbuf2d,  cam_out, &
     !
     integer :: c                                 ! chunk index
     integer :: ncol                              ! number of columns
+    integer :: nstep                             ! current timestep number
 #if (! defined SPMD)
     integer  :: mpicom = 0
 #endif
@@ -1177,11 +1184,29 @@ subroutine phys_run2(phys_state, ztodt, phys_tend, pbuf2d,  cam_out, &
     call t_startf ('ac_physics')
     !call t_adj_detailf(+1)
 
+    nstep = get_nstep()
+
+
+    !! calculate the global mean ieflx 
+
+    if(ieflx_opt>0) then
+       call ieflx_gmean(phys_state, phys_tend, pbuf2d, cam_in, cam_out, nstep)
+    end if
+
 !$OMP PARALLEL DO PRIVATE (C, NCOL, phys_buffer_chunk)
 
     do c=begchunk,endchunk
        ncol = get_ncols_p(c)
        phys_buffer_chunk => pbuf_get_chunk(pbuf2d, c)
+
+       !! 
+       !! add the implied internal energy flux to sensible heat flux
+       !! 
+
+       if(ieflx_opt>0) then
+          call check_ieflx_fix(c, ncol, nstep, cam_in(c)%shf)
+       end if
+
        !
        ! surface diagnostics for history files
        !