Merge branch 'bogensch/atm/SCMmods_for_v3' (PR #6254)

Single Column Model updates & fixes for v3

Modify and fix the single column model (SCM) to make it compatible with v3. Specifically:

1. The config_component for SCM_EAM has been changed to “scam_generic”
   (which should have never been changed).  When running cases with the standardized
   scripts prescribed aerosols are used, which is not compatible with the
   scm_generic_chemUCI-Linoz-mam5-vbs component currently specified.
   (Aside: a new prescribed aerosol file has been generated by running E3SM with mam5.
   This file has been uploaded to the E3SM data server and the E3SM SCM scripts
   have been updated to use this file).

2. Idealization flags have been added to the P3 microphysics to allow for:
    i. Runs with no precipitation.  Note that in the P3 implementation this only
       turns off liquid precipitation, which is used in several GCSS boundary
       layer cloud cases in the E3SM SCM library.
   ii. Prescribed droplet concentration.  This was already a feature in place
       in P3 but this PR connects this functionality to the EAM namelist.

3. Fixes to allow the ability to “replay” a column once again (back by semi-popular demand!).
   This was a feature that worked in E3SMv1 but has not been functional for sometime.  The fix in
   this PR allows users to generate output needed to make IOP forcing files to replay a single column.
   Note that additional post-processing is required to actually do an SCM run (in which the user
   needs to consult with the E3SM wiki page).

Further documentation and scripts have been updated on the E3SM SCM wiki to reflect changes made in this PR
(https://github.com/E3SM-Project/scmlib/wiki/E3SM-Single-Column-Model-Home).

[BFB]

components/eam/bld/namelist_files/namelist_defaults_eam.xml

@@ -917,6 +917,7 @@
 <cld_sed              microphys="mg2">  1.0D0       </cld_sed>
 <mg_prc_coeff_fix     microphys="mg2">  .false.     </mg_prc_coeff_fix>
 <micro_nccons         microphys="mg2">  100.D6      </micro_nccons>
+<micro_nccons         microphys="p3">   200.D6      </micro_nccons>
 <micro_nicons         microphys="mg2">  0.1D6       </micro_nicons>
 <micro_mincdnc        microphys="mg2">  -999.       </micro_mincdnc>
 <p3_mincdnc           microphys="p3">  -999.        </p3_mincdnc> 

components/eam/bld/namelist_files/namelist_definition.xml

@@ -4726,7 +4726,9 @@ Default: FALSE
 
 <entry id="precip_off" type="logical" category="scm"
        group="cam_inparm" valid_values="">
-Turn off microphysics computation.
+Turn off microphysics computation for MG2.  For P3 this disables the generation of
+liquid precipitation via autoconversion only, to be used for idealized simulations of
+warm phase boundary layer clouds.
 Default: FALSE
 </entry>                                                                        |
                                                                                 |

components/eam/bld/namelist_files/use_cases/scam_generic.xml

@@ -3,7 +3,6 @@
 
 <namelist_defaults>
 
-<ncdata             >atm/cam/inic/homme/cami_mam3_Linoz_ne4np4_L72_c160909.nc</ncdata>
 <soag_ext_file      >atm/cam/chem/trop_mozart_aero/emis/aces4bgc_nvsoa_soag_elev_2000_c160427.nc</soag_ext_file>
 
 <!-- SCM is configured with ne4 grid, but we want settings for an ne30 simulation -->

components/eam/cime_config/config_component.xml

@@ -136,7 +136,7 @@
       <value compset="20TR(?:SOI)?_EAM.*CHEMUCI.*LINOZ.*MAM5"  >20TR_eam_chemUCI-Linoz-mam5</value>
       <value compset="ARM95_EAM"                        >scam_arm95</value>
       <value compset="ARM97_EAM"                        >scm_arm97_chemUCI-Linoz-mam5-vbs</value>
-      <value compset="SCM_EAM"                          >scm_generic_chemUCI-Linoz-mam5-vbs</value>
+      <value compset="SCM_EAM"                          >scam_generic</value>
       <value compset="PIPD_EAM"                         >1850-PD_cam5</value>
       <value compset="_EAM%AQP_"                        >aquaplanet_EAMv1</value>
       <value compset="_EAM%RCE_"                        >RCEMIP_EAMv1</value>

components/eam/src/dynamics/se/se_iop_intr_mod.F90

@@ -56,7 +56,7 @@ subroutine iop_setinitial(elem)
   integer i, j, k, cix, ie, thelev
   integer inumliq, inumice, icldliq, icldice
 
-  if (.not. use_replay .and. get_nstep() .eq. 0 .and. par%dynproc) then
+  if (get_nstep() .eq. 0 .and. par%dynproc) then
     call cnst_get_ind('NUMLIQ', inumliq, abrtf=.false.)
     call cnst_get_ind('NUMICE', inumice, abrtf=.false.)
     call cnst_get_ind('CLDLIQ', icldliq)
@@ -215,8 +215,7 @@ subroutine iop_setfield(elem,iop_update_phase1)
   integer i, j, k, ie
 
   do ie=1,nelemd
-    if (have_ps .and. use_replay .and. .not. iop_update_phase1) elem(ie)%state%ps_v(:,:,:) = psobs
-    if (have_ps .and. .not. use_replay) elem(ie)%state%ps_v(:,:,:) = psobs
+    if (have_ps) elem(ie)%state%ps_v(:,:,:) = psobs
     do i=1, PLEV
       ! If DP CRM mode do NOT write over dycore vertical velocity
       if ((have_omega .and. iop_update_phase1) .and. .not. dp_crm) elem(ie)%derived%omega_p(:,:,i)=wfld(i)  !     set t to tobs at first

components/eam/src/dynamics/se/stepon.F90

@@ -526,9 +526,10 @@ subroutine stepon_run3(dtime, cam_out, phys_state, dyn_in, dyn_out)
    use hycoef,      only: hyam, hybm
    use dimensions_mod, only: nlev, nelemd, np, npsq
    use se_iop_intr_mod, only: iop_setfield, iop_setinitial
-   use dyn_comp, only: TimeLevel
+   use dyn_comp, only: TimeLevel, hvcoord
    use cam_history,     only: outfld   
    use cam_logfile, only: iulog
+   use element_ops,     only: get_temperature
    use mpishorthand
    real(r8), intent(in) :: dtime   ! Time-step
    real(r8) :: ftmp_temp(np,np,nlev,nelemd), ftmp_q(np,np,nlev,pcnst,nelemd)
@@ -542,6 +543,7 @@ subroutine stepon_run3(dtime, cam_out, phys_state, dyn_in, dyn_out)
    type (dyn_export_t), intent(inout) :: dyn_out ! Dynamics export container
    type (element_t), pointer :: elem(:)
    integer :: rc, i, j, k, p, ie, tl_f
+   real(r8) :: temperature(np,np,nlev)   ! Temperature from dynamics
 #if defined (E3SM_SCM_REPLAY)
    real(r8) :: forcing_temp(npsq,nlev), forcing_q(npsq,nlev,pcnst)
 #endif   
@@ -554,7 +556,10 @@ subroutine stepon_run3(dtime, cam_out, phys_state, dyn_in, dyn_out)
 
    ! Save ftmp stuff to get state before dynamics is called
    do ie=1,nelemd
-     ftmp_temp(:,:,:,ie) = dyn_in%elem(ie)%state%T(:,:,:,tl_f)
+
+     call get_temperature(dyn_in%elem(ie),temperature,hvcoord,tl_f)
+
+     ftmp_temp(:,:,:,ie) = temperature(:,:,:)
      ftmp_q(:,:,:,:,ie) = dyn_in%elem(ie)%state%Q(:,:,:,:)
    enddo
 
@@ -599,16 +604,20 @@ subroutine stepon_run3(dtime, cam_out, phys_state, dyn_in, dyn_out)
    tl_f = TimeLevel%n0
 
    do ie=1,nelemd
+
+     ! Get temperature from dynamics state
+     call get_temperature(dyn_in%elem(ie),temperature,hvcoord,tl_f)
+
      do k=1,nlev
        do j=1,np
          do i=1,np
 
            ! Note that this calculation will not provide b4b results with 
            !  an E3SM because the dynamics tendency is not computed in the exact
            !  same way as an E3SM run, introducing error with roundoff 
-           forcing_temp(i+(j-1)*np,k) = (dyn_in%elem(ie)%state%T(i,j,k,tl_f) - &
+           forcing_temp(i+(j-1)*np,k) = (temperature(i,j,k) - &
              ftmp_temp(i,j,k,ie))/dtime - dyn_in%elem(ie)%derived%FT(i,j,k)
-           out_temp(i+(j-1)*np,k) = dyn_in%elem(ie)%state%T(i,j,k,tl_f)
+           out_temp(i+(j-1)*np,k) = temperature(i,j,k)
            out_u(i+(j-1)*np,k) = dyn_in%elem(ie)%state%v(i,j,1,k,tl_f)
            out_v(i+(j-1)*np,k) = dyn_in%elem(ie)%state%v(i,j,2,k,tl_f)
            out_q(i+(j-1)*np,k) = dyn_in%elem(ie)%state%Q(i,j,k,1)

components/eam/src/physics/p3/eam/micro_p3.F90

@@ -56,7 +56,7 @@ module micro_p3
    use phys_control,  only: use_hetfrz_classnuc
 
    ! physical and mathematical constants
-   use micro_p3_utils, only: rho_1000mb,rho_600mb,ar,br,f1r,f2r,rho_h2o,kr,kc,aimm,mi0,nccnst,  &
+   use micro_p3_utils, only: rho_1000mb,rho_600mb,ar,br,f1r,f2r,rho_h2o,kr,kc,aimm,mi0, &
        eci,eri,bcn,cpw,cons1,cons3,cons4,cons5,cons6,cons7,         &
        inv_rho_h2o,inv_dropmass,qsmall,nsmall,cp,g,rd,rv,ep_2,inv_cp,   &
        thrd,sxth,piov6,rho_rimeMin,     &
@@ -441,15 +441,15 @@ SUBROUTINE p3_main_part1(kts, kte, kbot, ktop, kdir, do_predict_nc, do_prescribe
        t_atm, rho, inv_rho, qv_sat_l, qv_sat_i, qv_supersat_i, rhofacr, rhofaci, acn, qv, th_atm,           &
        qc, nc, qr, nr,                                                                                      &
        qi, ni, qm, bm, qc_incld, qr_incld, qi_incld, qm_incld,                                              &
-       nc_incld, nr_incld, ni_incld, bm_incld, is_nucleat_possible, is_hydromet_present)
+       nc_incld, nr_incld, ni_incld, bm_incld, is_nucleat_possible, is_hydromet_present, nccnst)
 
     implicit none
 
     ! args
 
     integer, intent(in) :: kts, kte, kbot, ktop, kdir
     logical(btype), intent(in) :: do_predict_nc
-    real(rtype), intent(in) :: dt
+    real(rtype), intent(in) :: dt, nccnst
 
     real(rtype), intent(in), dimension(kts:kte) :: pres, dpres, dz, nc_nuceat_tend, exner, inv_exner, &
          inv_cld_frac_l, inv_cld_frac_i, inv_cld_frac_r, latent_heat_vapor, latent_heat_sublim, latent_heat_fusion, nccn_prescribed
@@ -561,15 +561,15 @@ SUBROUTINE p3_main_part2(kts, kte, kbot, ktop, kdir, do_predict_nc, do_prescribe
        qm, bm, latent_heat_vapor, latent_heat_sublim, latent_heat_fusion, qc_incld, qr_incld, qi_incld, qm_incld, nc_incld, nr_incld, &
        ni_incld, bm_incld, mu_c, nu, lamc, cdist, cdist1, cdistr, mu_r, lamr, logn0r, qv2qi_depos_tend, precip_total_tend, &
        nevapr, qr_evap_tend, vap_liq_exchange, vap_ice_exchange, liq_ice_exchange, pratot, &
-       prctot, frzimm, frzcnt, frzdep, p3_tend_out, is_hydromet_present)
+       prctot, frzimm, frzcnt, frzdep, p3_tend_out, is_hydromet_present, do_precip_off, nccnst)
 
     implicit none
 
     ! args
 
     integer, intent(in) :: kts, kte, kbot, ktop, kdir
-    logical(btype), intent(in) :: do_predict_nc, do_prescribed_CCN
-    real(rtype), intent(in) :: dt, inv_dt
+    logical(btype), intent(in) :: do_predict_nc, do_prescribed_CCN, do_precip_off
+    real(rtype), intent(in) :: dt, inv_dt, nccnst
     real(rtype), intent(in) :: p3_autocon_coeff, p3_accret_coeff, p3_qc_autocon_expon, p3_nc_autocon_expon, p3_qc_accret_expon, &
          p3_wbf_coeff, p3_embryonic_rain_size, p3_max_mean_rain_size
 
@@ -875,7 +875,7 @@ SUBROUTINE p3_main_part2(kts, kte, kbot, ktop, kdir, do_predict_nc, do_prescribe
       ! NOTE: cloud_water_autoconversion must be called before droplet_self_collection
       call cloud_water_autoconversion(rho(k),qc_incld(k),nc_incld(k),inv_qc_relvar(k),&
            p3_autocon_coeff,p3_qc_autocon_expon,p3_nc_autocon_expon,p3_embryonic_rain_size,&
-           qc2qr_autoconv_tend,nc2nr_autoconv_tend,ncautr)
+           do_precip_off,qc2qr_autoconv_tend,nc2nr_autoconv_tend,ncautr)
 
       !............................
       ! self-collection of droplets
@@ -978,7 +978,7 @@ SUBROUTINE p3_main_part2(kts, kte, kbot, ktop, kdir, do_predict_nc, do_prescribe
       !-- warm-phase only processes:
       call update_prognostic_liquid(qc2qr_accret_tend, nc_accret_tend, qc2qr_autoconv_tend, nc2nr_autoconv_tend, ncautr, &
            nc_selfcollect_tend, qr2qv_evap_tend, nr_evap_tend, nr_selfcollect_tend,           &
-           do_predict_nc, do_prescribed_CCN, inv_rho(k), exner(k), latent_heat_vapor(k), dt,                     &
+           do_predict_nc, nccnst, do_prescribed_CCN, inv_rho(k), exner(k), latent_heat_vapor(k), dt,                     &
            th_atm(k), qv(k), qc(k), nc(k), qr(k), nr(k))
 
       !==
@@ -1259,7 +1259,7 @@ SUBROUTINE p3_main(qc,nc,qr,nr,th_atm,qv,dt,qi,qm,ni,bm,
        p3_wbf_coeff,p3_mincdnc,p3_max_mean_rain_size,p3_embryonic_rain_size,                                                                                             &
        dpres,exner,qv2qi_depos_tend,precip_total_tend,nevapr,qr_evap_tend,precip_liq_flux,precip_ice_flux,rflx,sflx,cflx,cld_frac_r,cld_frac_l,cld_frac_i,               &
        p3_tend_out,mu_c,lamc,liq_ice_exchange,vap_liq_exchange,                                                                                                          &
-       vap_ice_exchange,qv_prev,t_prev,col_location,diag_equiv_reflectivity,diag_ze_rain,diag_ze_ice                                                                     &
+       vap_ice_exchange,qv_prev,t_prev,col_location,do_precip_off,nccnst,diag_equiv_reflectivity,diag_ze_rain,diag_ze_ice                                                                     &
 #ifdef SCREAM_CONFIG_IS_CMAKE
        ,elapsed_s &
 #endif
@@ -1337,6 +1337,10 @@ SUBROUTINE p3_main(qc,nc,qr,nr,th_atm,qv,dt,qi,qm,ni,bm,
     ! INPUT for prescribed CCN option
     logical(btype), intent(in)                                  :: do_prescribed_CCN
 
+    ! INPUT for idealization options
+    logical(btype), intent(in)                                  :: do_precip_off
+    real(rtype), intent(in)                                     :: nccnst
+
     ! INPUT for p3 tuning parameters
     real(rtype), intent(in)                                     :: p3_autocon_coeff         ! autconversion coefficient
     real(rtype), intent(in)                                     :: p3_accret_coeff          ! accretion coefficient
@@ -1517,7 +1521,7 @@ SUBROUTINE p3_main(qc,nc,qr,nr,th_atm,qv,dt,qi,qm,ni,bm,
             rhofaci(i,:), acn(i,:), qv(i,:), th_atm(i,:), qc(i,:), nc(i,:), qr(i,:), nr(i,:), &
             qi(i,:), ni(i,:), qm(i,:), bm(i,:), qc_incld(i,:), qr_incld(i,:), &
             qi_incld(i,:), qm_incld(i,:), nc_incld(i,:), nr_incld(i,:), &
-            ni_incld(i,:), bm_incld(i,:), is_nucleat_possible, is_hydromet_present)
+            ni_incld(i,:), bm_incld(i,:), is_nucleat_possible, is_hydromet_present, nccnst)
 
       if (debug_ON) then
          tmparr1(i,:) = th_atm(i,:)*inv_exner(i,:)!(pres(i,:)*1.e-5)**(rd*inv_cp)
@@ -1541,7 +1545,8 @@ SUBROUTINE p3_main(qc,nc,qr,nr,th_atm,qv,dt,qi,qm,ni,bm,
             bm_incld(i,:), mu_c(i,:), nu(i,:), lamc(i,:), cdist(i,:), cdist1(i,:), &
             cdistr(i,:), mu_r(i,:), lamr(i,:), logn0r(i,:), qv2qi_depos_tend(i,:), precip_total_tend(i,:), &
             nevapr(i,:), qr_evap_tend(i,:), vap_liq_exchange(i,:), vap_ice_exchange(i,:), &
-            liq_ice_exchange(i,:), pratot(i,:), prctot(i,:), frzimm(i,:), frzcnt(i,:), frzdep(i,:), p3_tend_out(i,:,:), is_hydromet_present)
+            liq_ice_exchange(i,:), pratot(i,:), prctot(i,:), frzimm(i,:), frzcnt(i,:), frzdep(i,:), p3_tend_out(i,:,:), is_hydromet_present, &
+	    do_precip_off, nccnst)
 
 
        ! measure microphysics processes tendency output
@@ -2950,7 +2955,7 @@ end subroutine rain_self_collection
 
 subroutine cloud_water_autoconversion(rho,qc_incld,nc_incld,inv_qc_relvar,          &
    p3_autocon_coeff,p3_qc_autocon_expon,p3_nc_autocon_expon,p3_embryonic_rain_size, &
-   qc2qr_autoconv_tend,nc2nr_autoconv_tend,ncautr)
+   do_precip_off,qc2qr_autoconv_tend,nc2nr_autoconv_tend,ncautr)
 
    implicit none
 
@@ -2963,6 +2968,8 @@ subroutine cloud_water_autoconversion(rho,qc_incld,nc_incld,inv_qc_relvar,
    real(rtype), intent(in) :: p3_nc_autocon_expon
    real(rtype), intent(in) :: p3_embryonic_rain_size
 
+   logical(btype), intent(in) :: do_precip_off
+
    real(rtype), intent(out) :: qc2qr_autoconv_tend
    real(rtype), intent(out) :: nc2nr_autoconv_tend
    real(rtype), intent(out) :: ncautr
@@ -2981,8 +2988,9 @@ subroutine cloud_water_autoconversion(rho,qc_incld,nc_incld,inv_qc_relvar,
       ncautr = qc2qr_autoconv_tend*cons3*(1._rtype/bfb_pow(p3_embryonic_rain_size,3._rtype))
       nc2nr_autoconv_tend = qc2qr_autoconv_tend*nc_incld/qc_incld
 
-      if (qc2qr_autoconv_tend .eq.0._rtype) nc2nr_autoconv_tend = 0._rtype
-      if (nc2nr_autoconv_tend.eq.0._rtype) qc2qr_autoconv_tend  = 0._rtype
+      if (qc2qr_autoconv_tend .eq.0._rtype .or. do_precip_off) nc2nr_autoconv_tend = 0._rtype
+      if (nc2nr_autoconv_tend.eq.0._rtype .or. do_precip_off) qc2qr_autoconv_tend  = 0._rtype
+      if (do_precip_off) ncautr = 0._rtype
 
    endif qc_not_small
 
@@ -3502,7 +3510,7 @@ end subroutine update_prognostic_ice
 
 subroutine update_prognostic_liquid(qc2qr_accret_tend,nc_accret_tend,qc2qr_autoconv_tend,nc2nr_autoconv_tend, &
      ncautr,nc_selfcollect_tend, qr2qv_evap_tend,nr_evap_tend,nr_selfcollect_tend,         &
-    do_predict_nc, do_prescribed_CCN, inv_rho,exner,latent_heat_vapor,dt,                                      &
+    do_predict_nc, nccnst, do_prescribed_CCN, inv_rho,exner,latent_heat_vapor,dt,          &
     th_atm,qv,qc,nc,qr,nr)
 
    !-- warm-phase only processes:
@@ -3520,6 +3528,7 @@ subroutine update_prognostic_liquid(qc2qr_accret_tend,nc_accret_tend,qc2qr_autoc
 
 
    logical(btype), intent(in) :: do_predict_nc, do_prescribed_CCN
+   real(rtype), intent(in) :: nccnst
    real(rtype), intent(in) :: inv_rho
    real(rtype), intent(in) :: exner
    real(rtype), intent(in) :: latent_heat_vapor

components/eam/src/physics/p3/eam/micro_p3_interface.F90

@@ -42,6 +42,7 @@ module micro_p3_interface
   use ncdio_atm,       only: infld
   use ppgrid,         only: begchunk, endchunk, pcols, pver, pverp,psubcols
   use cam_history_support, only: add_hist_coord
+  use iop_data_mod,   only: precip_off
 
   implicit none
   save
@@ -131,8 +132,8 @@ module micro_p3_interface
       p3_wbf_coeff             = huge(1.0_rtype), &
       p3_mincdnc               = huge(1.0_rtype), & 
       p3_max_mean_rain_size    = huge(1.0_rtype), &
-      p3_embryonic_rain_size   = huge(1.0_rtype)
-
+      p3_embryonic_rain_size   = huge(1.0_rtype), &
+      micro_nccons             = huge(1.0_rtype)
 
    integer :: ncnst
 
@@ -165,7 +166,7 @@ subroutine micro_p3_readnl(nlfile)
        micro_p3_tableversion, micro_p3_lookup_dir, micro_aerosolactivation, micro_subgrid_cloud, &
        micro_tend_output, p3_autocon_coeff, p3_qc_autocon_expon, p3_nc_autocon_expon, p3_accret_coeff, &
        p3_qc_accret_expon, p3_wbf_coeff, p3_max_mean_rain_size, p3_embryonic_rain_size, &
-       do_prescribed_CCN, do_Cooper_inP3, p3_mincdnc 
+       do_prescribed_CCN, do_Cooper_inP3, p3_mincdnc, micro_nccons
 
   !-----------------------------------------------------------------------------
 
@@ -220,6 +221,7 @@ subroutine micro_p3_readnl(nlfile)
   call mpibcast(p3_embryonic_rain_size,  1 ,                         mpir8,   0, mpicom)
   call mpibcast(do_prescribed_CCN,       1,                          mpilog,  0, mpicom)
   call mpibcast(do_Cooper_inP3,          1,                          mpilog,  0, mpicom)
+  call mpibcast(micro_nccons,            1,                          mpir8,   0, mpicom)
 
 #endif
 
@@ -1363,6 +1365,8 @@ subroutine micro_p3_tend(state, ptend, dtime, pbuf)
          qv_prev(its:ite,kts:kte),         & ! IN  qv at end of prev p3_main call   kg kg-1
          t_prev(its:ite,kts:kte),          & ! IN  t at end of prev p3_main call    K
          col_location(its:ite,:3),         & ! IN column locations
+         precip_off,                       & ! IN Option to turn precip (liquid) off
+         micro_nccons,                     & ! IN Option for constant droplet concentration
          diag_equiv_reflectivity(its:ite,kts:kte), & !OUT equivalent reflectivity (rain + ice) [dBz]
          diag_ze_rain(its:ite,kts:kte),diag_ze_ice(its:ite,kts:kte)) !OUT equivalent reflectivity for rain and ice [dBz]
 

components/eam/src/physics/p3/eam/micro_p3_utils.F90

@@ -33,9 +33,6 @@ module micro_p3_utils
     ! maximum total ice concentration (sum of all categories)
     real(rtype), public, parameter :: max_total_ni = 500.e+3_rtype  ! (m)
 
-    ! droplet concentration (m-3)
-    real(rtype), public, parameter :: nccnst = 200.e+6_rtype
-
     ! parameters for Seifert and Beheng (2001) autoconversion/accretion
     real(rtype), public, parameter :: kc     = 9.44e+9_rtype
     real(rtype), public, parameter :: kr     = 5.78e+3_rtype