diff --git a/components/cam/bld/namelist_files/namelist_definition.xml b/components/cam/bld/namelist_files/namelist_definition.xml index 0a114a82b61b..fbd2548d4eb1 100644 --- a/components/cam/bld/namelist_files/namelist_definition.xml +++ b/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. + +IEFLX fixer options +Default: 0 + + + +Whether or not to fix ieflx +Default: .false. + + 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) !!................................................................... diff --git a/components/cam/src/physics/cam/phys_control.F90 b/components/cam/src/physics/cam/phys_control.F90 index 3f4cfdcb9762..9accf1f67738 100644 --- a/components/cam/src/physics/cam/phys_control.F90 +++ b/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 diff --git a/components/cam/src/physics/cam/physpkg.F90 b/components/cam/src/physics/cam/physpkg.F90 index 478a208f9c7b..10b440a34b18 100644 --- a/components/cam/src/physics/cam/physpkg.F90 +++ b/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 !