Merge branch 'beharrop/mmf/add_macv2-sp_rebase' (E3SM-Project#5905)

Adds in the MACv2-SP aerosol scheme for a simplified form of anthropogenic aerosols
that combines with the prescribed natural aerosol. The current implementation is only
for the E3SM-MMF. These changes are controlled via atm namelist options, the
default setting is off for the MACv2 scheme.

[stealth]

components/eam/bld/namelist_files/namelist_definition.xml

@@ -5528,6 +5528,25 @@ Rinaldi et al. (JGR, 2013).
 Default:1
 </entry>
 
+<entry id="prescribed_macv2_datapath" type="char*256" input_pathname="abs" category="cam_chem"
+       group="prescribed_macv2_nl" valid_values="" >
+Full pathname of the directory that contains the files specified in
+<varname>prescribed_macv2_filelist</varname>.
+Default: set by build-namelist.
+</entry>
+
+<entry id="prescribed_macv2_file" type="char*256" input_pathname="rel:prescribed_macv2_datapath" category="cam_chem"
+       group="prescribed_macv2_nl" valid_values="" >
+Filename of dataset for prescribed macv2 anthropogenic aerosol forcing.
+Default: set by build-namelist.
+</entry>
+
+<entry id="do_macv2sp" type="logical"  category="cam_chem"
+       group="prescribed_macv2_nl" valid_values="">
+If set to .true., use MACv2-SP scheme for anthropogenic aerosol effects.
+Default: .false.
+</entry>
+
 <entry id="prescribed_ghg_datapath" type="char*256" input_pathname="abs" category="cam_chem"
        group="prescribed_ghg_nl" valid_values="" >
 Full pathname of the directory that contains the files specified in

components/eam/src/control/runtime_opts.F90

@@ -279,6 +279,9 @@ subroutine read_namelist(single_column_in, scmlon_in, scmlat_in, scm_multcols_in
 #endif
    use radiation,           only: radiation_readnl
    use conditional_diag,    only: cnd_diag_readnl
+#if defined(MMF_SAMXX) || defined(MMF_PAM)
+   use prescribed_macv2,    only: prescribed_macv2_readnl
+#endif
 
 !---------------------------Arguments-----------------------------------
 
@@ -542,6 +545,10 @@ subroutine read_namelist(single_column_in, scmlon_in, scmlat_in, scm_multcols_in
    ! Read radiation namelist
    call radiation_readnl(nlfilename, dtime_in=dtime)
 
+#if defined(MMF_SAMXX) || defined(MMF_PAM)
+   call prescribed_macv2_readnl(nlfilename)
+#endif
+
    ! Print cam_inparm input variables to standard output
    if (masterproc) then
       write(iulog,*)' ------------------------------------------'

components/eam/src/physics/crm/physpkg.F90

@@ -549,6 +549,7 @@ subroutine phys_init( phys_state, phys_tend, pbuf2d, cam_out )
   use modal_aero_wateruptake,only: modal_aero_wateruptake_init
   use nucleate_ice_cam,      only: nucleate_ice_cam_init
   use hetfrz_classnuc_cam,   only: hetfrz_classnuc_cam_init
+  use prescribed_macv2,      only: macv2_rad_props_init
   !-----------------------------------------------------------------------------
   ! Input/output arguments
   !-----------------------------------------------------------------------------
@@ -630,6 +631,9 @@ subroutine phys_init( phys_state, phys_tend, pbuf2d, cam_out )
   ! Initialize ocean data
   if (has_mam_mom) call init_ocean_data()
 
+  ! Initialize MACv2-SP aerosols
+  call macv2_rad_props_init()
+
   ! co2 cycle            
   if (co2_transport()) call co2_init()
   call co2_diags_init(phys_state)

components/eam/src/physics/crm/rrtmgp/radiation.F90

@@ -25,7 +25,7 @@ module radiation
       get_sw_spectral_boundaries, &
       rrtmg_to_rrtmgp_swbands
    use cam_history_support, only: add_hist_coord
-   use physconst, only: cpair, cappa, stebol
+   use physconst, only: cpair, cappa, gravit
 
    ! RRTMGP gas optics object to store coefficient information. This is imported
    ! here so that we can make the k_dist objects module data and only load them
@@ -901,6 +901,8 @@ subroutine radiation_init(state)
                      sampling_seq='rad_lwsw', flag_xyfill=.true.)
       endif
 
+
+
    end subroutine radiation_init
 
    subroutine radiation_final()
@@ -1103,7 +1105,7 @@ subroutine radiation_tend(state_in,ptend,    pbuf,          cam_out, cam_in,  &
 
       ! For running CFMIP Observation Simulator Package (COSP)
       use cospsimulator_intr, only: docosp, cospsimulator_intr_run, cosp_nradsteps
-
+      
       ! ---------------------------------------------------------------------------
       ! Arguments
       ! ---------------------------------------------------------------------------
@@ -1274,6 +1276,8 @@ subroutine radiation_tend(state_in,ptend,    pbuf,          cam_out, cam_in,  &
       ! Loop variables
       integer :: icol, ilay, iday
 
+
+
       !----------------------------------------------------------------------
 
       ! Copy state so we can use CAM routines with arrays replaced with data
@@ -1371,13 +1375,16 @@ subroutine radiation_tend(state_in,ptend,    pbuf,          cam_out, cam_in,  &
                aer_tau_bnd_lw = 0
                if (do_aerosol_rad) then
                   if (radiation_do('sw')) then
+
                      call t_startf('rad_aerosol_optics_sw')
                      call set_aerosol_optics_sw( &
                           icall, dt, state, pbuf, night_indices(1:nnight), is_cmip6_volc, &
                           aer_tau_bnd_sw, aer_ssa_bnd_sw, aer_asm_bnd_sw,  &
                           clear_rh=clear_rh)
                      ! Now reorder bands to be consistent with RRTMGP
+
                      ! TODO: fix the input files themselves!
+                     ! Note that MACv2 may need changing too if the input files are reordered!
                      do icol = 1,size(aer_tau_bnd_sw,1)
                         do ilay = 1,size(aer_tau_bnd_sw,2)
                            aer_tau_bnd_sw(icol,ilay,:) = reordered(aer_tau_bnd_sw(icol,ilay,:), rrtmg_to_rrtmgp_swbands)
@@ -1386,6 +1393,9 @@ subroutine radiation_tend(state_in,ptend,    pbuf,          cam_out, cam_in,  &
                         end do
                      end do
                      call t_stopf('rad_aerosol_optics_sw')
+
+
+
                   end if ! radiation_do('sw')
                   if (radiation_do('lw')) then
                      call t_startf('rad_aerosol_optics_lw')

components/eam/src/physics/rrtmgp/cam_optics.F90

@@ -450,7 +450,11 @@ subroutine set_aerosol_optics_sw(icall, dt, state, pbuf, &
       use physics_types, only: physics_state
       use physics_buffer, only: physics_buffer_desc
       use aer_rad_props, only: aer_rad_props_sw
-      use radconstants, only: nswbands
+      use radconstants, only: nswbands, idx_sw_diag
+#if defined(MMF_SAMXX) || defined(MMF_PAM)
+      use prescribed_macv2, only: do_macv2sp
+#endif
+      use cam_history, only: outfld
       integer, intent(in) :: icall
       real(r8), intent(in):: dt
       type(physics_state), intent(in) :: state
@@ -469,13 +473,24 @@ subroutine set_aerosol_optics_sw(icall, dt, state, pbuf, &
       !
       ! NOTE: dimension ordering is different than for cloud optics!
       real(r8), dimension(pcols,0:pver,nswbands) :: tau, tau_w, tau_w_g, tau_w_f
+      ! NOTE: MACv2 array only has pver levels, so initialize differently
+      real(r8), dimension(pcols,1:pver,nswbands) :: mac_tau, mac_tau_w, &
+           mac_tau_w_g, mac_tau_w_f
+      real(r8) :: itau
+      real(r8) :: itauw(2)
+      real(r8) :: itaug(2)
 
       integer :: ncol
       integer :: icol, ilay
+      integer :: isw              ! additional loop indices
 
       ! Everyone needs a name
       character(len=*), parameter :: subroutine_name = 'set_aerosol_optics_sw'
 
+#if !defined(MMF_SAMXX) && !defined(MMF_PAM)
+      logical :: do_macv2sp = .false.
+#endif
+
       ncol = state%ncol
 
       ! Get aerosol absorption optical depth from CAM routine
@@ -487,6 +502,59 @@ subroutine set_aerosol_optics_sw(icall, dt, state, pbuf, &
            count(night_indices > 0), night_indices, is_cmip6_volc, &
            tau, tau_w, tau_w_g, tau_w_f, clear_rh=clear_rh)
 
+      ! Combine MACv2-SP aerosol optical properties with natural aerosols
+      ! note that AEROD_v is already written to history file in aer_rad_props_sw
+      ! so the MACv2SP effect is not included
+
+      if ( (icall == 0) .AND. (do_macv2sp) ) then
+         mac_tau     = 0._r8
+         mac_tau_w   = 0._r8
+         mac_tau_w_g = 0._r8
+         mac_tau_w_f = 0._r8
+
+         call outfld('NAT_TAU',     tau(    :,1:pver,idx_sw_diag), pcols, state%lchnk)
+         call outfld('NAT_TAU_W',   tau_w(  :,1:pver,idx_sw_diag), pcols, state%lchnk)
+         call outfld('NAT_TAU_W_G', tau_w_g(:,1:pver,idx_sw_diag), pcols, state%lchnk)
+         call outfld('NAT_TAU_W_F', tau_w_f(:,1:pver,idx_sw_diag), pcols, state%lchnk)
+
+         call set_macv2_aerosol_optics(state, pbuf, mac_tau, &
+              mac_tau_w, mac_tau_w_g, mac_tau_w_f)
+
+         call outfld('MAC_TAU',     mac_tau(    :,1:pver,idx_sw_diag), pcols, state%lchnk)
+         call outfld('MAC_TAU_W',   mac_tau_w(  :,1:pver,idx_sw_diag), pcols, state%lchnk)
+         call outfld('MAC_TAU_W_G', mac_tau_w_g(:,1:pver,idx_sw_diag), pcols, state%lchnk)
+         call outfld('MAC_TAU_W_F', mac_tau_w_f(:,1:pver,idx_sw_diag), pcols, state%lchnk)
+
+         do isw = 1, nswbands
+            do ilay = 1, pver
+               do icol = 1, ncol
+                  ! A quick reference for how to combine optical properties
+                  ! t = tau; w = SSA; g = asymmetry parameter
+                  ! t = t1 + t2
+                  ! w = (w1*t1 + w2*t2) / (t1 + t2)
+                  ! g = (g1*w1*t1 + g2*w2*t2) / (w1*t1 + w2*t2)
+                  !
+                  ! So if tau_w_g = (g1*w1*t1 + g2*tw*t2)
+                  !       tau_w   = (w1*t1 + w2*t2)
+                  !       tau     = (t1 + t2)
+                  ! Then, g = tau_w_g / tau_w, w = tau_w / tau, and t = tau, as desired
+                  tau_w_f(icol,ilay,isw) = tau_w_f(icol,ilay,isw) + mac_tau_w_f(icol,ilay,isw)
+                  tau_w_g(icol,ilay,isw) = tau_w_g(icol,ilay,isw) + mac_tau_w_g(icol,ilay,isw)
+                  tau_w(  icol,ilay,isw) = tau_w(  icol,ilay,isw) + mac_tau_w(  icol,ilay,isw)
+                  tau(    icol,ilay,isw) = tau(    icol,ilay,isw) + mac_tau(    icol,ilay,isw)
+
+               end do
+            end do
+         end do
+
+         !these 3 variables has pver+1 levels
+         call outfld('AER_TAU',     tau(    :,1:pver,idx_sw_diag), pcols, state%lchnk)
+         call outfld('AER_TAU_W',   tau_w(  :,1:pver,idx_sw_diag), pcols, state%lchnk)
+         call outfld('AER_TAU_W_G', tau_w_g(:,1:pver,idx_sw_diag), pcols, state%lchnk)
+         call outfld('AER_TAU_W_F', tau_w_f(:,1:pver,idx_sw_diag), pcols, state%lchnk)
+
+      end if
+
       ! Extract quantities from products
       do icol = 1,ncol
          ! Copy cloud optical depth over directly
@@ -536,6 +604,99 @@ subroutine set_aerosol_optics_lw(icall, dt, state, pbuf, is_cmip6_volc, tau)
 
    end subroutine set_aerosol_optics_lw
 
+   !----------------------------------------------------------------------------
+
+   subroutine set_macv2_aerosol_optics(state, pbuf, &
+                                       tau, tau_w, tau_w_g, tau_w_f)
+
+      use ppgrid,           only: pcols, pver
+      use physics_types,    only: physics_state
+      use physics_buffer,   only: physics_buffer_desc
+      use aer_rad_props,    only: aer_rad_props_sw
+      use radconstants,     only: nswbands, wavenum_sw_lower, wavenum_sw_upper
+#if defined(MMF_SAMXX) || defined(MMF_PAM)
+      use prescribed_macv2, only: do_macv2sp, sp_aop_profile, swbandnum
+#endif
+      use time_manager,     only: get_curr_date, get_curr_calday
+      use physconst,        only: gravit
+      use cam_history,      only: addfld, outfld
+      use phys_grid,        only: get_rlat_all_p, get_rlon_all_p
+
+      type(physics_state), intent(in)         :: state
+      type(physics_buffer_desc), pointer      :: pbuf(:)
+      real(r8), intent(out), dimension(:,:,:) :: tau, tau_w, tau_w_g, tau_w_f
+
+      ! add variables for MACv2-SP
+      integer  :: yr, mon, day     ! year, month, and day components of date
+      integer  :: ncsec            ! current time of day [seconds]
+      real(r8) :: calday           ! current calendar day
+      real(r8) :: clat(pcols)      ! current latitudes(radians)
+      real(r8) :: clon(pcols)      ! current longitudes(radians)
+      integer  :: isw              ! additional loop indices
+      real(r8) :: lambda           !SW wavelengh input to MACV2
+      real(r8) :: year_fr
+      real(r8) :: aod_prof(pcols,pver,nswbands)  ! profile of aerosol optical depth
+      real(r8) :: ssa_prof(pcols,pver,nswbands)  ! profile of single scattering albedo
+      real(r8) :: asy_prof(pcols,pver,nswbands)  ! profile of asymmetry parameter
+      integer  :: ncol
+      integer  :: icol, ilay
+
+      ! Subroutine name for error messages
+      character(len=*), parameter :: subroutine_name = 'set_macv2_aerosol_optics'
+
+#if !defined(MMF_SAMXX) && !defined(MMF_PAM)
+      logical :: do_macv2sp = .false.
+      character(len=5) :: swbandnum(nswbands) =(/'_sw01','_sw02','_sw03','_sw04','_sw05','_sw06','_sw07','_sw08','_sw09','_sw10','_sw11','_sw12','_sw13','_sw14'/)
+#endif
+
+      if ( .not. do_macv2sp ) return
+
+      ! calculate MACv2-SP aerosol direct effects 
+      ! get year fraction for MACv2-SP's prescribed annual cycle and year-to-year variability
+
+      ncol    = state%ncol
+
+      call get_curr_date(yr, mon, day, ncsec)
+      calday  = get_curr_calday()
+      year_fr = yr + calday/365.0_r8
+
+      call get_rlat_all_p(state%lchnk, ncol, clat(1:ncol))
+      call get_rlon_all_p(state%lchnk, ncol, clon(1:ncol))
+
+      do isw = 1, nswbands
+
+      !get the bin-center wavelength from the wave number bin  (in the units of nm)
+      !this is an input to the sp_aop_profile subroutine of MACv2-SP
+         lambda = 10.0_r8**7*( wavenum_sw_lower(isw)**(-1) + wavenum_sw_upper(isw)**(-1) )/2.0_r8 
+
+#if defined(MMF_SAMXX) || defined(MMF_PAM)
+         call sp_aop_profile (ncol, lambda, state%phis/gravit, clon, clat, year_fr, state%zm, &
+              aod_prof(:,:,isw), ssa_prof(:,:,isw), asy_prof(:,:,isw), state%lchnk, isw)
+#endif
+      ! state%phis/gravit for orography in [m]
+
+      ! output diagnostic variables for MACv2-SP, only for the mid-visible wavelength
+         call outfld('MACv2_aod'//swbandnum(isw),  aod_prof(:,:,isw),  pcols, state%lchnk)
+         call outfld('MACv2_ssa'//swbandnum(isw),  ssa_prof(:,:,isw),  pcols, state%lchnk)
+         call outfld('MACv2_asy'//swbandnum(isw),  asy_prof(:,:,isw),  pcols, state%lchnk)
+
+      ! prepare arrays to be combined with EAM's aerosol optical parameters 
+         ! e.g, single-scattering albedo multiplied by optical depth
+         do icol = 1, ncol
+            do ilay = 1, pver
+               tau(    icol, ilay, isw) = aod_prof(icol, ilay, isw)
+               tau_w(  icol, ilay, isw) = aod_prof(icol, ilay, isw) * ssa_prof(icol, ilay, isw)
+               tau_w_g(icol, ilay, isw) = tau_w(   icol, ilay, isw) * asy_prof(icol, ilay, isw)
+               tau_w_f(icol, ilay, isw) = tau_w_g( icol, ilay, isw) * asy_prof(icol, ilay, isw)
+            end do
+         end do
+
+      end do ! isw = 1, nswbands
+
+
+
+   end subroutine set_macv2_aerosol_optics
+
    !----------------------------------------------------------------------------
    !----------------------------------------------------------------------------
 

components/eam/src/physics/rrtmgp/radconstants.F90

@@ -40,10 +40,10 @@ module radconstants
 ! properties are output before bands are reordered to the expected RRTMGP order.
 ! Optical properties should be reordered to RRTMGP order in the RRTMGP driver
 ! interface codes.
-real(r8), parameter :: wavenum_sw_lower(nswbands) = & ! in cm^-1
+real(r8), parameter, public :: wavenum_sw_lower(nswbands) = & ! in cm^-1
   (/2600._r8, 3250._r8, 4000._r8, 4650._r8, 5150._r8, 6150._r8, 7700._r8, &
     8050._r8,12850._r8,16000._r8,22650._r8,29000._r8,38000._r8,  820._r8/)
-real(r8), parameter :: wavenum_sw_upper(nswbands) = & ! in cm^-1
+real(r8), parameter, public :: wavenum_sw_upper(nswbands) = & ! in cm^-1
   (/3250._r8, 4000._r8, 4650._r8, 5150._r8, 6150._r8, 7700._r8, 8050._r8, &
    12850._r8,16000._r8,22650._r8,29000._r8,38000._r8,50000._r8, 2600._r8/)
 real(r8), parameter :: wavenum_lw_lower(nlwbands) = &! Longwave spectral band limits (cm-1)