Merge branch 'huiwanpnnl/atm/v3atm_cflx_new_imp' into NGD_v3atm (PR #17)

A revised handling of cflx as a stealth feature

This PR adds an optional revision to the numerical process coupling
at the tphysbc/tphysac level to couple the surface emissions of tracers
more tightly with the turbulent transport. This is done by revising where
in the time integration loop the surface fluxes in cam_in%cflx(:,2:) are
applied to update the tracer mixing ratios.

The revision was motivated by the known issues of overly short dust lifetime
in EAMv1 and the very strong sensitivity of the dust lifetime to vertical
resolution, see Feng et al., 2022 for the v1 results.

The subroutine clubb_surface in clubb_intr.F90 was split into two parts:

The calculations of ustar and the Obukhov length remain in clubb_surface.
The calculation of tracer mixing ratio tendencies using cam_in%clfx was
moved to a new subroutine cflx_tend in a new module in cflx.F90
A new namelist variable cflx_cpl_opt was added.

cflx_cpl_opt = 1 (default) gives the original coupling in v2.
There, the subroutine cflx_tend is called immediately after clubb_surface
followed by a call physics_update(...). This gives the original process
coupling in EAMv2 and the results are BFB.

cflx_cpl_opt = 2 activates the revised coupling.
The call of cflx_tend and the corresponding call physics_update(...) are
moved to right before the cloud macro-microphysics sub-cycles in tphysbc,
so that the tracers emitted to the lowest model layer can be transported
to upper layers by CLUBB or dropmixnuc.

Note that cflx_cpl_opt = 2 does not affect water vapor. Neither does it
affect tracers whose portion of cam_in%cflx is zero when clubb_surface
is called.

[BFB] with default cflx_cpl_opt = 1

components/eam/bld/namelist_files/namelist_definition.xml

@@ -2644,6 +2644,12 @@ Do radiative transfer and heating/cooling calculation.
 Default: TRUE
 </entry>
 
+<entry id="cflx_cpl_opt" type="integer" category="process_coupling"
+       group="phys_ctl_nl" valid_values="" >
+Numerical scheme for coupling surface tracer fluxes with the rest of model.
+Default: 1 
+</entry>
+
 <!-- Physics sub-column switches -->
 
 <entry id="use_subcol_microp" type="logical" category="conv"

components/eam/src/cpl/atm_import_export.F90

@@ -52,7 +52,19 @@ subroutine atm_import( x2a, cam_in, restart_init )
        ! NOTE:overwrite_flds is .FALSE. for the first restart
        ! time step making cflx(:,1)=0.0 for the first restart time step.
        ! cflx(:,1) should not be zeroed out, start the second index of cflx from 2.
-       cam_in(c)%cflx(:,2:) = 0._r8 
+
+       ! +++ Update from 2022-09 +++
+       ! For some of the new process coupling options in EAM, some of the constituents'
+       ! cam_in%cflx are used not in tphysac but in the tphysbc call of the next time step. 
+       ! This means for an exact restart, we also need to write out cam_in%cflx(:,2:)
+       ! and then read them back in. Because the present subroutine is called after 
+       ! the cflx variables are read in in the subroutine read_restart_physics 
+       ! in physics/cam/restart_physics.F90, we need to move the following line 
+       ! to that read_restart_physics to avoid incorrectly zeroing out the needed values.
+       !
+       !cam_in(c)%cflx(:,2:) = 0._r8 
+       !
+       ! === Update from 2022-09 ===
 
        do i =1,ncols                                                               
           if (overwrite_flds) then

components/eam/src/physics/cam/cflx.F90

@@ -0,0 +1,100 @@
+module cflx
+
+  !-----------------------------------------------------------------------------------------------------
+  ! History:
+  !  Separated from module clubb_intr, subroutine clubb_surface by Hui Wan (PNNL), 2022
+  !-----------------------------------------------------------------------------------------------------
+
+  use shr_kind_mod,  only: r8=>shr_kind_r8
+
+  implicit none
+  public
+
+contains
+
+    subroutine cflx_tend (state, cam_in, ztodt, ptend)
+
+    use physics_types,          only: physics_state, physics_ptend, &
+                                      physics_ptend_init, &
+                                      set_dry_to_wet, set_wet_to_dry
+    use physconst,              only: gravit
+    use ppgrid,                 only: pver, pcols
+    use constituents,           only: pcnst, cnst_get_ind, cnst_type
+    use co2_cycle,              only: co2_cycle_set_cnst_type
+    use camsrfexch,             only: cam_in_t
+
+    implicit none
+
+    ! Input Auguments
+
+    type(physics_state), intent(inout)  :: state                ! Physics state variables
+    type(cam_in_t),      intent(in)     :: cam_in               ! contains surface fluxes of constituents
+    real(r8),            intent(in)     :: ztodt                ! 2 delta-t        [ s ]
+
+    ! Output Auguments
+
+    type(physics_ptend), intent(out)    :: ptend                ! Individual parameterization tendencies
+
+    ! Local Variables
+
+    integer :: i                                                ! indicees
+    integer :: ncol                                             ! # of atmospheric columns
+
+    real(r8) :: tmp1(pcols)
+    real(r8) :: rztodt                                          ! 1./ztodt
+    integer  :: m
+
+    logical  :: lq(pcnst)
+
+    character(len=3), dimension(pcnst) :: cnst_type_loc         ! local override option for constituents cnst_type
+
+
+    ncol = state%ncol
+
+    !-------------------------------------------------------
+    ! Assume 'wet' mixing ratios in surface diffusion code.
+    ! don't convert co2 tracers to wet mixing ratios
+
+    cnst_type_loc(:) = cnst_type(:)
+    call co2_cycle_set_cnst_type(cnst_type_loc, 'wet')
+    call set_dry_to_wet(state, cnst_type_loc)
+
+    !-------------------------------------------------------
+    ! Initialize ptend
+
+    lq(:) = .TRUE.
+    call physics_ptend_init(ptend, state%psetcols, 'clubb_srf', lq=lq)
+
+    !-------------------------------------------------------
+    ! Calculate tracer mixing ratio tendencies from cflx
+
+    rztodt                 = 1._r8/ztodt
+    ptend%q(:ncol,:pver,:) = state%q(:ncol,:pver,:)
+    tmp1(:ncol)            = ztodt * gravit * state%rpdel(:ncol,pver)
+
+    do m = 2, pcnst
+      ptend%q(:ncol,pver,m) = ptend%q(:ncol,pver,m) + tmp1(:ncol) * cam_in%cflx(:ncol,m)
+    enddo
+
+    ptend%q(:ncol,:pver,:) = (ptend%q(:ncol,:pver,:) - state%q(:ncol,:pver,:)) * rztodt
+
+    ! Convert tendencies of dry constituents to dry basis.
+    do m = 1,pcnst
+       if (cnst_type(m).eq.'dry') then
+          ptend%q(:ncol,:pver,m) = ptend%q(:ncol,:pver,m)*state%pdel(:ncol,:pver)/state%pdeldry(:ncol,:pver)
+       endif
+    end do
+
+    !-------------------------------------------------------
+    ! convert wet mmr back to dry before conservation check
+    ! avoid converting co2 tracers again
+
+    cnst_type_loc(:) = cnst_type(:)
+    call co2_cycle_set_cnst_type(cnst_type_loc, 'wet')
+    call set_wet_to_dry(state, cnst_type_loc)
+
+    return
+
+  end subroutine cflx_tend
+
+end module

components/eam/src/physics/cam/clubb_intr.F90

@@ -3055,10 +3055,10 @@ end subroutine clubb_tend_cam
   !                                                                                 !
   ! =============================================================================== !
 
-    subroutine clubb_surface (state, ptend, ztodt, cam_in, ustar, obklen)
+    subroutine clubb_surface (state, cam_in, ustar, obklen)
 
 !-------------------------------------------------------------------------------
-! Description: Provide the obukov length and the surface friction velocity
+! Description: Provide the obukhov length and the surface friction velocity
 !              for the dry deposition code in routine tphysac.  Since University
 !              of Washington Moist Turbulence (UWMT) scheme is not called when
 !              CLUBB is turned on the obukov length and ustar are never initialized
@@ -3071,13 +3071,10 @@ subroutine clubb_surface (state, ptend, ztodt, cam_in, ustar, obklen)
 !   None
 !-------------------------------------------------------------------------------
 
-    use physics_types,          only: physics_state, physics_ptend, &
-                                      physics_ptend_init, &
-                                      set_dry_to_wet, set_wet_to_dry
-    use physconst,              only: gravit, zvir, latvap
+    use physics_types,          only: physics_state
+    use physconst,              only: zvir
     use ppgrid,                 only: pver, pcols
-    use constituents,           only: pcnst, cnst_get_ind, cnst_type
-    use co2_cycle,              only: co2_cycle_set_cnst_type
+    use constituents,           only: cnst_get_ind
     use camsrfexch,             only: cam_in_t
 
     implicit none
@@ -3089,13 +3086,10 @@ subroutine clubb_surface (state, ptend, ztodt, cam_in, ustar, obklen)
     type(physics_state), intent(inout)  :: state                ! Physics state variables
     type(cam_in_t),      intent(in)     :: cam_in
 
-    real(r8),            intent(in)     :: ztodt                ! 2 delta-t        [ s ]
-
     ! ---------------- !
     ! Output Auguments !
     ! ---------------- !
 
-    type(physics_ptend), intent(out)    :: ptend                ! Individual parameterization tendencies
     real(r8),            intent(out)    :: obklen(pcols)        ! Obukhov length [ m ]
     real(r8),            intent(out)    :: ustar(pcols)         ! Surface friction velocity [ m/s ]
 
@@ -3113,15 +3107,9 @@ subroutine clubb_surface (state, ptend, ztodt, cam_in, ustar, obklen)
     real(r8) :: kinheat                                         ! kinematic surface heat flux
     real(r8) :: kinwat                                          ! kinematic surface vapor flux
     real(r8) :: kbfs                                            ! kinematic surface buoyancy flux
-    real(r8) :: tmp1(pcols)
-    real(r8) :: rztodt                                          ! 1./ztodt
-    integer  :: m
     integer  :: ixq,ixcldliq !PMA fix for thv
     real(r8) :: rrho                                            ! Inverse air density
 
-    logical  :: lq(pcnst)
-
-    character(len=3), dimension(pcnst) :: cnst_type_loc         ! local override option for constituents cnst_type
 
 #endif
     obklen(pcols) = 0.0_r8
@@ -3131,21 +3119,11 @@ subroutine clubb_surface (state, ptend, ztodt, cam_in, ustar, obklen)
     ! ----------------------- !
     ! Main Computation Begins !
     ! ----------------------- !
-
-    ! Assume 'wet' mixing ratios in surface diffusion code.
-    ! don't convert co2 tracers to wet mixing ratios
-    cnst_type_loc(:) = cnst_type(:)
-    call co2_cycle_set_cnst_type(cnst_type_loc, 'wet')
-    call set_dry_to_wet(state, cnst_type_loc)
-
     call cnst_get_ind('Q',ixq)
     if (use_sgv) then
        call cnst_get_ind('CLDLIQ',ixcldliq)
     endif
 
-    lq(:) = .TRUE.
-    call physics_ptend_init(ptend, state%psetcols, 'clubb_srf', lq=lq)
-
     ncol = state%ncol
 
     ! Compute the surface friction velocity and obukov length
@@ -3167,28 +3145,6 @@ subroutine clubb_surface (state, ptend, ztodt, cam_in, ustar, obklen)
                         kinheat, kinwat, kbfs, obklen(i) )
     enddo
 
-    rztodt                 = 1._r8/ztodt
-    ptend%q(:ncol,:pver,:) = state%q(:ncol,:pver,:)
-    tmp1(:ncol)            = ztodt * gravit * state%rpdel(:ncol,pver)
-
-    do m = 2, pcnst
-      ptend%q(:ncol,pver,m) = ptend%q(:ncol,pver,m) + tmp1(:ncol) * cam_in%cflx(:ncol,m)
-    enddo
-
-    ptend%q(:ncol,:pver,:) = (ptend%q(:ncol,:pver,:) - state%q(:ncol,:pver,:)) * rztodt
-
-    ! Convert tendencies of dry constituents to dry basis.
-    do m = 1,pcnst
-       if (cnst_type(m).eq.'dry') then
-          ptend%q(:ncol,:pver,m) = ptend%q(:ncol,:pver,m)*state%pdel(:ncol,:pver)/state%pdeldry(:ncol,:pver)
-       endif
-    end do
-    ! convert wet mmr back to dry before conservation check
-    ! avoid converting co2 tracers again
-    cnst_type_loc(:) = cnst_type(:)
-    call co2_cycle_set_cnst_type(cnst_type_loc, 'wet')
-    call set_wet_to_dry(state, cnst_type_loc)
-
     return
 
 #endif

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

@@ -189,6 +189,11 @@ module phys_control
 logical :: l_st_mic        = .true.
 logical :: l_rad           = .true.
 
+! Numerical schemes for process coupling
+
+integer :: cflx_cpl_opt = 1  ! When to apply surface tracer fluxes (not including water vapor).
+                             ! The default for aerosols is to do this 
+                             ! after tphysac:clubb_surface and before aerosol dry removal.
 
 !======================================================================= 
 contains
@@ -233,6 +238,7 @@ subroutine phys_ctl_readnl(nlfile)
       fix_g1_err_ndrop, ssalt_tuning, resus_fix, convproc_do_aer, &
       convproc_do_gas, convproc_method_activate, liqcf_fix, regen_fix, demott_ice_nuc, pergro_mods, pergro_test_active, &
       mam_amicphys_optaa, n_so4_monolayers_pcage,micro_mg_accre_enhan_fac, &
+      cflx_cpl_opt, &
       l_tracer_aero, l_vdiff, l_rayleigh, l_gw_drag, l_ac_energy_chk, &
       l_bc_energy_fix, l_dry_adj, l_st_mac, l_st_mic, l_rad, prc_coef1,prc_exp,prc_exp1,cld_sed,mg_prc_coeff_fix, &
       rrtmg_temp_fix
@@ -336,6 +342,7 @@ subroutine phys_ctl_readnl(nlfile)
    call mpibcast(demott_ice_nuc,                  1 , mpilog,  0, mpicom)
    call mpibcast(pergro_mods,                     1 , mpilog,  0, mpicom)
    call mpibcast(pergro_test_active,              1 , mpilog,  0, mpicom)
+   call mpibcast(cflx_cpl_opt,                    1 , mpiint,  0, mpicom)
    call mpibcast(l_tracer_aero,                   1 , mpilog,  0, mpicom)
    call mpibcast(l_vdiff,                         1 , mpilog,  0, mpicom)
    call mpibcast(l_rayleigh,                      1 , mpilog,  0, mpicom)
@@ -430,6 +437,12 @@ subroutine phys_ctl_readnl(nlfile)
       end if
    end if
 
+   ! Check settings for process coupling
+   if (masterproc) write(iulog,*) '**** phys_ctl_readnl: cflx_cpl_opt = ', cflx_cpl_opt
+   if (.not.( (cflx_cpl_opt==1).or.(cflx_cpl_opt==2) )) then
+      call endrun('phys_ctl_readnl: unsupported value of cflx_cpl_opt')
+   end if
+
    ! prog_modal_aero determines whether prognostic modal aerosols are present in the run.
    prog_modal_aero = (     cam_chempkg_is('trop_mam3') &
                       .or. cam_chempkg_is('trop_mam4') &
@@ -519,6 +532,7 @@ subroutine phys_getopts(deep_scheme_out, shallow_scheme_out, eddy_scheme_out, &
                         fix_g1_err_ndrop_out, ssalt_tuning_out,resus_fix_out,convproc_do_aer_out,  &
                         convproc_do_gas_out, convproc_method_activate_out, mam_amicphys_optaa_out, n_so4_monolayers_pcage_out, &
                         micro_mg_accre_enhan_fac_out, liqcf_fix_out, regen_fix_out,demott_ice_nuc_out, pergro_mods_out, pergro_test_active_out &
+                       ,cflx_cpl_opt_out &
                        ,l_tracer_aero_out, l_vdiff_out, l_rayleigh_out, l_gw_drag_out, l_ac_energy_chk_out  &
                        ,l_bc_energy_fix_out, l_dry_adj_out, l_st_mac_out, l_st_mic_out, l_rad_out  &
                        ,prc_coef1_out,prc_exp_out,prc_exp1_out, cld_sed_out,mg_prc_coeff_fix_out,rrtmg_temp_fix_out)
@@ -608,7 +622,7 @@ subroutine phys_getopts(deep_scheme_out, shallow_scheme_out, eddy_scheme_out, &
    logical,           intent(out), optional :: pergro_mods_out     
    logical,           intent(out), optional :: pergro_test_active_out     
 
-
+   integer,           intent(out), optional :: cflx_cpl_opt_out
    logical,           intent(out), optional :: l_tracer_aero_out
    logical,           intent(out), optional :: l_vdiff_out
    logical,           intent(out), optional :: l_rayleigh_out
@@ -706,6 +720,7 @@ subroutine phys_getopts(deep_scheme_out, shallow_scheme_out, eddy_scheme_out, &
    if ( present(pergro_mods_out         ) ) pergro_mods_out          = pergro_mods
    if ( present(pergro_test_active_out  ) ) pergro_test_active_out   = pergro_test_active
 
+   if ( present(cflx_cpl_opt_out        ) ) cflx_cpl_opt_out      = cflx_cpl_opt
    if ( present(l_tracer_aero_out       ) ) l_tracer_aero_out     = l_tracer_aero
    if ( present(l_vdiff_out             ) ) l_vdiff_out           = l_vdiff
    if ( present(l_rayleigh_out          ) ) l_rayleigh_out        = l_rayleigh