Merge pull request E3SM-Project#13 from apcraig/aerclean

Update implementation of g(k), w0(k) in icepack_shortwave, clean up options

columnphysics/icepack_parameters.F90

@@ -214,7 +214,7 @@ module icepack_parameters
          awtidf = 0.36218_dbl_kind   ! near IR, diffuse
 
       character (len=char_len), public :: &
-         shortwave   = 'dEdd', & ! shortwave method, 'ccsm3' or 'dEdd' or 'dEdd_snicar'
+         shortwave   = 'dEdd', & ! shortwave method, 'ccsm3' or 'dEdd' or 'dEdd_snicar_ad'
          albedo_type = 'ccsm3'   ! albedo parameterization, 'ccsm3' or 'constant'
                                  ! shortwave='dEdd' overrides this parameter
 
@@ -226,9 +226,7 @@ module icepack_parameters
          sw_frac      = 0.9_dbl_kind    , & ! Fraction of internal shortwave moved to surface
          sw_dtemp     = 0.02_dbl_kind       ! temperature difference from melting
 
-      ! Parameters for dEdd_snicar
-      logical (kind=log_kind), public :: &
-         use_snicar = .false.     ! .true. use 5-band SNICAR-AD approach
+      ! Parameters for dEdd_snicar_ad
       character (len=char_len), public :: &
          snw_ssp_table = 'test'   ! lookup table: 'snicar' or 'test'
 
@@ -462,7 +460,7 @@ subroutine icepack_init_parameters(   &
          update_ocn_f_in, ustar_min_in, a_rapid_mode_in, &
          Rac_rapid_mode_in, aspect_rapid_mode_in, &
          dSdt_slow_mode_in, phi_c_slow_mode_in, &
-         phi_i_mushy_in, shortwave_in, use_snicar_in, albedo_type_in, albsnowi_in, &
+         phi_i_mushy_in, shortwave_in, albedo_type_in, albsnowi_in, &
          albicev_in, albicei_in, albsnowv_in, &
          ahmax_in, R_ice_in, R_pnd_in, R_snw_in, dT_mlt_in, rsnw_mlt_in, &
          kalg_in, kstrength_in, krdg_partic_in, krdg_redist_in, mu_rdg_in, &
@@ -612,7 +610,7 @@ subroutine icepack_init_parameters(   &
          awtidf_in        ! near IR, diffuse
 
       character (len=*), intent(in), optional :: &
-         shortwave_in, & ! shortwave method, 'ccsm3' or 'dEdd' or 'dEdd_snicar'
+         shortwave_in, & ! shortwave method, 'ccsm3' or 'dEdd' or 'dEdd_snicar_ad'
          albedo_type_in  ! albedo parameterization, 'ccsm3' or 'constant'
                          ! shortwave='dEdd' overrides this parameter
 
@@ -635,7 +633,6 @@ subroutine icepack_init_parameters(   &
          kalg_in         ! algae absorption coefficient for 0.5 m thick layer
 
       logical (kind=log_kind), intent(in), optional :: &
-         use_snicar_in,& ! snicar adjustments to dEdd radiation for snow
          sw_redist_in    ! redistribute shortwave
 
       real (kind=dbl_kind), intent(in), optional :: &
@@ -939,7 +936,6 @@ subroutine icepack_init_parameters(   &
       if (present(phi_c_slow_mode_in)   ) phi_c_slow_mode  = phi_c_slow_mode_in
       if (present(phi_i_mushy_in)       ) phi_i_mushy      = phi_i_mushy_in
       if (present(shortwave_in)         ) shortwave        = shortwave_in
-      if (present(use_snicar_in)        ) use_snicar       = use_snicar_in
       if (present(albedo_type_in)       ) albedo_type      = albedo_type_in
       if (present(albicev_in)           ) albicev          = albicev_in
       if (present(albicei_in)           ) albicei          = albicei_in
@@ -1178,7 +1174,7 @@ subroutine icepack_query_parameters(   &
          Lfresh_out, cprho_out, Cp_out, ustar_min_out, a_rapid_mode_out, &
          ktherm_out, conduct_out, fbot_xfer_type_out, calc_Tsfc_out, dts_b_out, &
          Rac_rapid_mode_out, aspect_rapid_mode_out, dSdt_slow_mode_out, &
-         phi_c_slow_mode_out, phi_i_mushy_out, shortwave_out, use_snicar_out, &
+         phi_c_slow_mode_out, phi_i_mushy_out, shortwave_out, &
          albedo_type_out, albicev_out, albicei_out, albsnowv_out, &
          albsnowi_out, ahmax_out, R_ice_out, R_pnd_out, R_snw_out, dT_mlt_out, &
          rsnw_mlt_out, dEdd_algae_out, &
@@ -1338,7 +1334,7 @@ subroutine icepack_query_parameters(   &
          awtidf_out        ! near IR, diffuse
 
       character (len=*), intent(out), optional :: &
-         shortwave_out, & ! shortwave method, 'ccsm3' or 'dEdd' or 'dEdd_snicar'
+         shortwave_out, & ! shortwave method, 'ccsm3' or 'dEdd' or 'dEdd_snicar_ad'
          albedo_type_out  ! albedo parameterization, 'ccsm3' or 'constant'
                              ! shortwave='dEdd' overrides this parameter
 
@@ -1361,7 +1357,6 @@ subroutine icepack_query_parameters(   &
          kalg_out         ! algae absorption coefficient for 0.5 m thick layer
 
       logical (kind=log_kind), intent(out), optional :: &
-         use_snicar_out,& ! snicar adjustments to dEdd radiation for snow
          sw_redist_out    ! redistribute shortwave
 
       real (kind=dbl_kind), intent(out), optional :: &
@@ -1697,7 +1692,6 @@ subroutine icepack_query_parameters(   &
       if (present(phi_c_slow_mode_out)   ) phi_c_slow_mode_out = phi_c_slow_mode
       if (present(phi_i_mushy_out)       ) phi_i_mushy_out  = phi_i_mushy
       if (present(shortwave_out)         ) shortwave_out    = shortwave
-      if (present(use_snicar_out)        ) use_snicar_out   = use_snicar
       if (present(albedo_type_out)       ) albedo_type_out  = albedo_type
       if (present(albicev_out)           ) albicev_out      = albicev
       if (present(albicei_out)           ) albicei_out      = albicei
@@ -1908,7 +1902,6 @@ subroutine icepack_write_parameters(iounit)
         write(iounit,*) "  phi_c_slow_mode   = ", phi_c_slow_mode
         write(iounit,*) "  phi_i_mushy       = ", phi_i_mushy
         write(iounit,*) "  shortwave     = ", shortwave
-        write(iounit,*) "  use_snicar    = ", use_snicar
         write(iounit,*) "  albedo_type   = ", albedo_type
         write(iounit,*) "  albicev       = ", albicev
         write(iounit,*) "  albicei       = ", albicei

columnphysics/icepack_shortwave.F90

@@ -54,7 +54,7 @@ module icepack_shortwave
       use icepack_parameters, only: z_tracers, skl_bgc, calc_tsfc, shortwave, kalg
       use icepack_parameters, only: R_ice, R_pnd, R_snw, dT_mlt, rsnw_mlt, hs0, hs1, hp1
       use icepack_parameters, only: pndaspect, albedo_type, albicev, albicei, albsnowv, albsnowi, ahmax
-      use icepack_parameters, only: snw_ssp_table, use_snicar, modal_aero
+      use icepack_parameters, only: snw_ssp_table, modal_aero
       use icepack_parameters, only: dEdd_algae
 
       use icepack_tracers,    only: ncat, nilyr, nslyr, nblyr
@@ -173,7 +173,7 @@ subroutine icepack_init_radiation()
          if (icepack_warnings_aborted(subname)) return
       endif
 
-      if (use_snicar) then
+      if (trim(shortwave) == 'dEdd_snicar_ad') then
          call icepack_shortwave_init_dEdd5band()
          if (icepack_warnings_aborted(subname)) return
 
@@ -1507,7 +1507,7 @@ subroutine shortwave_dEdd  (coszen,                &
                ! calculate snow covered sea ice
 
                srftyp = 1
-               if (use_snicar) then
+               if (trim(shortwave) == 'dEdd_snicar_ad') then
                 call compute_dEdd_5bd(klev, klevp,                        &
                         zbio,                                             &
                         fnidr,     coszen,                                &
@@ -2361,12 +2361,10 @@ subroutine compute_dEdd_3bd( &
             ! aerosol in snow
             if (tr_zaero .and. dEdd_algae) then
                do k = 0,nslyr
-                  gzaer(ns,k) = gzaer(ns,k)/(wzaer(ns,k)+puny)
-                  wzaer(ns,k) = wzaer(ns,k)/(tzaer(ns,k)+puny)
-                  g  (k) = (g (k)*w0(k)*tau(k) + gzaer(ns,k)*wzaer(ns,k)*tzaer(ns,k)) &
-                         /       (w0(k)*tau(k) +             wzaer(ns,k)*tzaer(ns,k))
-                  w0 (k) = (      w0(k)*tau(k) +             wzaer(ns,k)*tzaer(ns,k)) &
-                         /       (      tau(k) +                         tzaer(ns,k))
+                  g(k)   = (g(k)*w0(k)*tau(k) + gzaer(ns,k)) / &
+                                (w0(k)*tau(k) + wzaer(ns,k))
+                  w0(k)  =      (w0(k)*tau(k) + wzaer(ns,k)) / &
+                                      (tau(k) + tzaer(ns,k))
                   tau(k) = tau(k) + tzaer(ns,k)
                enddo
             elseif (tr_aero) then
@@ -2412,6 +2410,8 @@ subroutine compute_dEdd_3bd( &
                gaer = gaer/(waer+puny)
                waer = waer/(taer+puny)
 
+! tcraig, why does the above section exist if taer=waer=gaer=0 below
+
                do k=1,nslyr
                   taer = c0
                   waer = c0
@@ -2461,12 +2461,10 @@ subroutine compute_dEdd_3bd( &
                             waer_3bd(ns,(1+(na-1)/4))*gaer_3bd(ns,(1+(na-1)/4))
                   endif       ! modal_aero
                   enddo       ! na
-                  gaer = gaer/(waer+puny)
-                  waer = waer/(taer+puny)
-                  g(k)   = (g(k)*w0(k)*tau(k) + gaer*waer*taer) / &
-                           (w0(k)*tau(k) + waer*taer)
-                  w0(k)  = (w0(k)*tau(k) + waer*taer) / &
-                           (tau(k) + taer)
+                  g(k)   = (g(k)*w0(k)*tau(k) + gaer) / &
+                                (w0(k)*tau(k) + waer)
+                  w0(k)  =      (w0(k)*tau(k) + waer) / &
+                                      (tau(k) + taer)
                   tau(k) = tau(k) + taer
                enddo       ! k
             endif     ! tr_aero
@@ -2522,12 +2520,10 @@ subroutine compute_dEdd_3bd( &
             ! aerosol in sea ice
             if (tr_zaero .and. dEdd_algae) then
                do k = kii, klev
-                  gzaer(ns,k) = gzaer(ns,k)/(wzaer(ns,k)+puny)
-                  wzaer(ns,k) = wzaer(ns,k)/(tzaer(ns,k)+puny)
-                  g(k)   = (g(k)*w0(k)*tau(k) + gzaer(ns,k)*wzaer(ns,k)*tzaer(ns,k)) &
-                         /      (w0(k)*tau(k) +             wzaer(ns,k)*tzaer(ns,k))
-                  w0(k)  = (w0(k)*tau(k) + wzaer(ns,k)*tzaer(ns,k)) &
-                         /       (tau(k) +             tzaer(ns,k))
+                  g(k)   = (g(k)*w0(k)*tau(k) + gzaer(ns,k)) / &
+                                (w0(k)*tau(k) + wzaer(ns,k))
+                  w0(k)  =      (w0(k)*tau(k) + wzaer(ns,k)) / &
+                                      (tau(k) + tzaer(ns,k))
                   tau(k) = tau(k) + tzaer(ns,k)
                enddo
             elseif (tr_aero) then
@@ -2578,12 +2574,10 @@ subroutine compute_dEdd_3bd( &
                 endif     ! modal_aero
                enddo      ! na
 
-               gaer = gaer/(waer+puny)
-               waer = waer/(taer+puny)
-               g(k)   = (g(k)*w0(k)*tau(k) + gaer*waer*taer) / &
-                    (w0(k)*tau(k) + waer*taer)
-               w0(k)  = (w0(k)*tau(k) + waer*taer) / &
-                    (tau(k) + taer)
+               g(k)   = (g(k)*w0(k)*tau(k) + gaer) / &
+                             (w0(k)*tau(k) + waer)
+               w0(k)  =      (w0(k)*tau(k) + waer) / &
+                                   (tau(k) + taer)
                tau(k) = tau(k) + taer
                do k = kii+1, klev
                   taer = c0
@@ -2631,12 +2625,10 @@ subroutine compute_dEdd_3bd( &
                           waer_3bd(ns,(1+(na-1)/4))*gaer_3bd(ns,(1+(na-1)/4))
                   endif       ! modal_aero
                   enddo       ! na
-                  gaer = gaer/(waer+puny)
-                  waer = waer/(taer+puny)
-                  g(k)   = (g(k)*w0(k)*tau(k) + gaer*waer*taer) / &
-                           (w0(k)*tau(k) + waer*taer)
-                  w0(k)  = (w0(k)*tau(k) + waer*taer) / &
-                           (tau(k) + taer)
+                  g(k)   = (g(k)*w0(k)*tau(k) + gaer) / &
+                                (w0(k)*tau(k) + waer)
+                  w0(k)  =      (w0(k)*tau(k) + waer) / &
+                                      (tau(k) + taer)
                   tau(k) = tau(k) + taer
                enddo ! k
             endif ! tr_aero
@@ -4958,9 +4950,9 @@ subroutine compute_dEdd_5bd (klev,  klevp, &
           if (tr_zaero .and. dEdd_algae) then
             do k = 0,nslyr
                g(k)   = (g(k)*w0(k)*tau(k) + gzaer_5bd(ns,k)) / &
-                               (w0(k)*tau(k) + wzaer_5bd(ns,k))
-               w0(k)  = (w0(k)*tau(k) + wzaer_5bd(ns,k)) / &
-                                (tau(k) + tzaer_5bd(ns,k))
+                             (w0(k)*tau(k) + wzaer_5bd(ns,k))
+               w0(k)  =      (w0(k)*tau(k) + wzaer_5bd(ns,k)) / &
+                                   (tau(k) + tzaer_5bd(ns,k))
                tau(k) = tau(k) + tzaer_5bd(ns,k)
             enddo
           elseif (tr_aero) then
@@ -5019,6 +5011,8 @@ subroutine compute_dEdd_5bd (klev,  klevp, &
                gaer = gaer/(waer+puny)
                waer = waer/(taer+puny)
 
+! tcraig, again why does the above exist if taer=waer=gaer=0 below
+
                do k=1,nslyr
                   taer = c0
                   waer = c0
@@ -5072,12 +5066,10 @@ subroutine compute_dEdd_5bd (klev,  klevp, &
                   endif       ! modal_aero
 !mgf--
                   enddo       ! na
-                  gaer = gaer/(waer+puny)
-                  waer = waer/(taer+puny)
-                  g(k)   = (g(k)*w0(k)*tau(k) + gaer*waer*taer) / &
-                           (w0(k)*tau(k) + waer*taer)
-                  w0(k)  = (w0(k)*tau(k) + waer*taer) / &
-                           (tau(k) + taer)
+                  g(k)   = (g(k)*w0(k)*tau(k) + gaer) / &
+                                (w0(k)*tau(k) + waer)
+                  w0(k)  =      (w0(k)*tau(k) + waer) / &
+                                      (tau(k) + taer)
                   tau(k) = tau(k) + taer
                enddo       ! k
             endif     ! tr_aero
@@ -5123,9 +5115,9 @@ subroutine compute_dEdd_5bd (klev,  klevp, &
                if (tr_zaero .and. dEdd_algae) then
                   do k = kii, klev
                      g(k)   = (g(k)*w0(k)*tau(k) + gzaer_5bd(ns,k)) / &
-                                     (w0(k)*tau(k) + wzaer_5bd(ns,k))
-                     w0(k)  = (w0(k)*tau(k) + wzaer_5bd(ns,k)) / &
-                                      (tau(k) + tzaer_5bd(ns,k))
+                                   (w0(k)*tau(k) + wzaer_5bd(ns,k))
+                     w0(k)  =      (w0(k)*tau(k) + wzaer_5bd(ns,k)) / &
+                                         (tau(k) + tzaer_5bd(ns,k))
                      tau(k) = tau(k) + tzaer_5bd(ns,k)
                   enddo
                elseif (tr_aero) then
@@ -5181,12 +5173,10 @@ subroutine compute_dEdd_5bd (klev,  klevp, &
    !mgf--
                   enddo      ! na
 
-                  gaer = gaer/(waer+puny)
-                  waer = waer/(taer+puny)
-                  g(k)   = (g(k)*w0(k)*tau(k) + gaer*waer*taer) / &
-                       (w0(k)*tau(k) + waer*taer)
-                  w0(k)  = (w0(k)*tau(k) + waer*taer) / &
-                       (tau(k) + taer)
+                  g(k)   = (g(k)*w0(k)*tau(k) + gaer) / &
+                                (w0(k)*tau(k) + waer)
+                  w0(k)  =      (w0(k)*tau(k) + waer) / &
+                                      (tau(k) + taer)
                   tau(k) = tau(k) + taer
                   do k = kii+1, klev
                      taer = c0
@@ -5241,12 +5231,10 @@ subroutine compute_dEdd_5bd (klev,  klevp, &
                      endif       ! modal_aero
    !mgf--
                      enddo       ! na
-                     gaer = gaer/(waer+puny)
-                     waer = waer/(taer+puny)
-                     g(k)   = (g(k)*w0(k)*tau(k) + gaer*waer*taer) / &
-                              (w0(k)*tau(k) + waer*taer)
-                     w0(k)  = (w0(k)*tau(k) + waer*taer) / &
-                              (tau(k) + taer)
+                     g(k)   = (g(k)*w0(k)*tau(k) + gaer) / &
+                                   (w0(k)*tau(k) + waer)
+                     w0(k)  =      (w0(k)*tau(k) + waer) / &
+                                         (tau(k) + taer)
                      tau(k) = tau(k) + taer
                   enddo ! k
                endif ! tr_aero

columnphysics/icepack_shortwave_data.F90

@@ -5,7 +5,6 @@ module icepack_shortwave_data
       use icepack_warnings, only: warnstr, icepack_warnings_add
       use icepack_warnings, only: icepack_warnings_setabort, icepack_warnings_aborted
       use icepack_parameters, only: c0
-      use icepack_parameters, only: use_snicar, snw_ssp_table
       use icepack_tracers, only: nmodal1, nmodal2, max_aero
 
       implicit none

configuration/driver/icedrv_init.F90

@@ -715,7 +715,7 @@ subroutine input_data
          write(nu_diag,1000) ' ahmax                     = ', ahmax
          endif
 
-         if (trim(shortwave) == 'dEdd_snicar') then
+         if (trim(shortwave) == 'dEdd_snicar_ad') then
          write(nu_diag,1030) ' snw_ssp_table             = ', trim(snw_ssp_table)
          endif
 

configuration/driver/icedrv_init_column.F90

@@ -120,7 +120,6 @@ subroutine init_shortwave
 
       logical (kind=log_kind) :: &
          l_print_point, & ! flag to print designated grid point diagnostics
-         use_snicar,    & ! use 5-band SNICAR radiation scheme for snow
          dEdd_algae,    & ! BGC - radiation interactions
          snwgrain         ! use variable snow grain size
 
@@ -221,9 +220,8 @@ subroutine init_shortwave
             if (icepack_warnings_aborted()) &
                call icedrv_system_abort(i, istep1, subname, __FILE__, __LINE__)
 
-            if (trim(shortwave) == 'dEdd_snicar') then
-               use_snicar = .true. ! 5-band SNICAR scheme for snow cover
-               call icepack_init_parameters(use_snicar_in=use_snicar, snw_ssp_table_in=snw_ssp_table)
+            if (trim(shortwave) == 'dEdd_snicar_ad') then
+               call icepack_init_parameters(snw_ssp_table_in=snw_ssp_table)
                call icepack_warnings_flush(nu_diag)
                if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, &
                    file=__FILE__,line= __LINE__)

configuration/scripts/options/set_nml.snicar

@@ -1,2 +1,2 @@
-  shortwave          = 'dEdd_snicar'
+  shortwave          = 'dEdd_snicar_ad'
   snw_ssp_table      = 'snicar'

configuration/scripts/options/set_nml.snicartest

@@ -1,2 +1,2 @@
-  shortwave       = 'dEdd_snicar'
+  shortwave       = 'dEdd_snicar_ad'
   snw_ssp_table   = 'test'