Update CICE to run with eclare108213/Icepack branch snicar

- including eclare108213/Icepack#13, Sept 11, 2022
- Passes full CICE test suite on cheyenne with 3 compilers except alt04 changes
  answers for all compilers and all tests.  CICE #fea412a55f was baseline.
- Icepack submodule still points to standard version on main, need to be
  swapped manually to appropriate development version.

cicecore/cicedynB/general/ice_forcing_bgc.F90

@@ -17,13 +17,13 @@ module ice_forcing_bgc
       use ice_calendar, only: dt, istep, msec, mday, mmonth
       use ice_fileunits, only: nu_diag
       use ice_arrays_column, only: restore_bgc, &
-         bgc_data_dir, fe_data_type, optics_file, optics_file_fieldname
+         bgc_data_dir, fe_data_type
       use ice_constants, only: c0, p1
       use ice_constants, only: field_loc_center, field_type_scalar
       use ice_exit, only: abort_ice
       use ice_forcing, only: bgc_data_type
       use icepack_intfc, only: icepack_warnings_flush, icepack_warnings_aborted
-      use icepack_intfc, only: icepack_nspint, icepack_max_aero, &
+      use icepack_intfc, only: icepack_nspint_3bd, icepack_max_aero, &
           icepack_max_algae, icepack_max_doc, icepack_max_dic
       use icepack_intfc, only: icepack_query_tracer_flags, &
           icepack_query_parameters, icepack_query_parameters, &
@@ -32,8 +32,7 @@ module ice_forcing_bgc
       implicit none
       private
       public :: get_forcing_bgc, get_atm_bgc, fzaero_data, alloc_forcing_bgc, &
-                init_bgc_data, faero_data, faero_default, faero_optics, &
-                fiso_default
+                init_bgc_data, faero_data, faero_default, fiso_default
 
       integer (kind=int_kind) :: &
          bgcrecnum = 0   ! old record number (save between steps)
@@ -840,6 +839,7 @@ subroutine init_bgc_data (fed1,fep1)
 
       end subroutine init_bgc_data
 
+#if (1 == 0)
 !=======================================================================
 !
 ! Aerosol optical properties for bulk and modal aerosol formulation
@@ -854,14 +854,6 @@ subroutine faero_optics
       use ice_broadcast, only: broadcast_array
       use ice_read_write, only: ice_open_nc, ice_close_nc
       use ice_communicate, only: my_task, master_task
-      use ice_arrays_column, only: &
-         kaer_tab, & ! aerosol mass extinction cross section (m2/kg)
-         waer_tab, & ! aerosol single scatter albedo (fraction)
-         gaer_tab, & ! aerosol asymmetry parameter (cos(theta))
-         kaer_bc_tab, & ! BC mass extinction cross section (m2/kg)
-         waer_bc_tab, & ! BC single scatter albedo (fraction)
-         gaer_bc_tab, & ! BC aerosol asymmetry parameter (cos(theta))
-         bcenh          ! BC absorption enhancement factor
 
 #ifdef USE_NETCDF
       use netcdf
@@ -895,23 +887,23 @@ subroutine faero_optics
            2665.85867,   2256.71027,    820.36024, &
             840.78295,   1028.24656,   1163.03298, &
             387.51211,    414.68808,    450.29814/), &
-            (/icepack_nspint,icepack_max_aero/))
+            (/icepack_nspint_3bd,icepack_max_aero/))
       waer_tab = reshape((/ &      ! aerosol single scatter albedo (fraction)
               0.29003,      0.17349,      0.06613, &
               0.51731,      0.41609,      0.21324, &
               0.84467,      0.94216,      0.95666, &
               0.97764,      0.99402,      0.98552, &
               0.94146,      0.98527,      0.99093, &
               0.90034,      0.96543,      0.97678/), &
-              (/icepack_nspint,icepack_max_aero/))
+              (/icepack_nspint_3bd,icepack_max_aero/))
       gaer_tab = reshape((/ &      ! aerosol asymmetry parameter (cos(theta))
               0.35445,      0.19838,      0.08857, &
               0.52581,      0.32384,      0.14970, &
               0.83162,      0.78306,      0.74375, &
               0.68861,      0.70836,      0.54171, &
               0.70239,      0.66115,      0.71983, &
               0.78734,      0.73580,      0.64411/), &
-              (/icepack_nspint,icepack_max_aero/))
+              (/icepack_nspint_3bd,icepack_max_aero/))
 
       ! this data is used in MODAL AEROSOL treatment in dEdd radiation
       kaer_bc_tab = reshape((/ &      ! aerosol mass extinction cross section (m2/kg)
@@ -925,7 +917,7 @@ subroutine faero_optics
               3934.17604,   4020.20799, 3543.27199, &
               3461.20656,   3587.80962, 3289.98060, &
               3083.03396,   3226.27231, 3052.91441/), &
-              (/icepack_nspint,10/))
+              (/icepack_nspint_3bd,10/))
 
       waer_bc_tab = reshape((/ &      ! aerosol single scatter albedo (fraction)
               0.26107,      0.15861,    0.06535, &
@@ -938,7 +930,7 @@ subroutine faero_optics
               0.49440,      0.46328,    0.42008, &
               0.50131,      0.47070,    0.43128, &
               0.50736,      0.47704,    0.44056/), &
-              (/icepack_nspint,10/))
+              (/icepack_nspint_3bd,10/))
 
       gaer_bc_tab = reshape((/ &      ! aerosol asymmetry parameter (cos(theta))
               0.28328,      0.19644,      0.10498, &
@@ -951,7 +943,7 @@ subroutine faero_optics
               0.75064,      0.64959,      0.47551, &
               0.76663,      0.67461,      0.50415, &
               0.77926,      0.69561,      0.52981/),&
-              (/icepack_nspint,10/))
+              (/icepack_nspint_3bd,10/))
 
       bcenh(:,:,:)     = c0
 
@@ -999,7 +991,7 @@ subroutine faero_optics
       endif      ! modal_aero
 
       end subroutine faero_optics
-
+#endif
 !=======================================================================
 
       end module ice_forcing_bgc

cicecore/cicedynB/general/ice_init.F90

@@ -158,7 +158,7 @@ subroutine input_data
 
       character (len=char_len) :: shortwave, albedo_type, conduct, fbot_xfer_type, &
         tfrz_option, frzpnd, atmbndy, wave_spec_type, snwredist, snw_aging_table, &
-        capping_method
+        capping_method, snw_ssp_table
 
       logical (kind=log_kind) :: calc_Tsfc, formdrag, highfreq, calc_strair, wave_spec, &
         sw_redist, calc_dragio, use_smliq_pnd, snwgrain
@@ -255,7 +255,7 @@ subroutine input_data
         Cf,             Pstar,          Cstar,          Ktens
 
       namelist /shortwave_nml/ &
-        shortwave,      albedo_type,                                    &
+        shortwave,      albedo_type,     snw_ssp_table,                 &
         albicev,        albicei,         albsnowv,      albsnowi,       &
         ahmax,          R_ice,           R_pnd,         R_snw,          &
         sw_redist,      sw_frac,         sw_dtemp,                      &
@@ -437,6 +437,7 @@ subroutine input_data
       advection  = 'remap'     ! incremental remapping transport scheme
       conserv_check = .false.  ! tracer conservation check
       shortwave = 'ccsm3'      ! 'ccsm3' or 'dEdd' (delta-Eddington)
+      snw_ssp_table = 'test'   ! 'test' or 'snicar' dEdd_snicar_ad table data
       albedo_type = 'ccsm3'    ! 'ccsm3' or 'constant'
 #ifdef UNDEPRECATE_0LAYER
       ktherm = 1               ! -1 = OFF, 0 = 0-layer, 1 = BL99, 2 = mushy thermo
@@ -902,6 +903,7 @@ subroutine input_data
       call broadcast_scalar(advection,            master_task)
       call broadcast_scalar(conserv_check,        master_task)
       call broadcast_scalar(shortwave,            master_task)
+      call broadcast_scalar(snw_ssp_table,        master_task)
       call broadcast_scalar(albedo_type,          master_task)
       call broadcast_scalar(ktherm,               master_task)
       call broadcast_scalar(coriolis,             master_task)
@@ -1188,7 +1190,7 @@ subroutine input_data
                write(nu_diag,*) subname//' ERROR: invalid seabed stress method'
                write(nu_diag,*) subname//' ERROR: seabed_stress_method should be LKD or probabilistic'
             endif
-            abort_list = trim(abort_list)//":34"
+            abort_list = trim(abort_list)//":48"
          endif
       endif
 
@@ -1283,10 +1285,10 @@ subroutine input_data
          abort_list = trim(abort_list)//":7"
       endif
 
-      if (trim(shortwave) /= 'dEdd' .and. tr_pond .and. calc_tsfc) then
+      if (shortwave(1:4) /= 'dEdd' .and. tr_pond .and. calc_tsfc) then
          if (my_task == master_task) then
             write(nu_diag,*) subname//' ERROR: tr_pond=T, calc_tsfc=T, invalid shortwave'
-            write(nu_diag,*) subname//' ERROR:   Must use shortwave=dEdd'
+            write(nu_diag,*) subname//' ERROR:   Must use shortwave=dEdd or dEdd_snicar_ad'
          endif
          abort_list = trim(abort_list)//":8"
       endif
@@ -1399,19 +1401,20 @@ subroutine input_data
          abort_list = trim(abort_list)//":36"
       endif
 
-      if (trim(shortwave) /= 'dEdd' .and. tr_aero) then
+      if (shortwave(1:4) /= 'dEdd' .and. tr_aero) then
          if (my_task == master_task) then
             write(nu_diag,*) subname//' ERROR: tr_aero=T, invalid shortwave'
-            write(nu_diag,*) subname//' ERROR:   Must use shortwave=dEdd'
+            write(nu_diag,*) subname//' ERROR:   Must use shortwave=dEdd or dEdd_snicar_ad'
          endif
          abort_list = trim(abort_list)//":10"
       endif
 
-      if (trim(shortwave) /= 'dEdd' .and. snwgrain) then
+      if (shortwave(1:4) /= 'dEdd' .and. snwgrain) then
          if (my_task == master_task) then
-            write (nu_diag,*) 'WARNING: snow grain radius activated but'
-            write (nu_diag,*) 'WARNING: dEdd shortwave is not.'
+            write (nu_diag,*) subname//' ERROR: snow grain radius is activated'
+            write (nu_diag,*) subname//' ERROR:   Must use shortwave=dEdd or dEdd_snicar_ad'
          endif
+         abort_list = trim(abort_list)//":29"
       endif
 
       if ((rfracmin < -puny .or. rfracmin > c1+puny) .or. &
@@ -1689,7 +1692,7 @@ subroutine input_data
          write(nu_diag,1020) ' nilyr            = ', nilyr, ' : number of ice layers (equal thickness)'
          write(nu_diag,1020) ' nslyr            = ', nslyr, ' : number of snow layers (equal thickness)'
          write(nu_diag,1020) ' nblyr            = ', nblyr, ' : number of bio layers (equal thickness)'
-         if (trim(shortwave) == 'dEdd') &
+         if (shortwave(1:4) == 'dEdd') &
          write(nu_diag,*) 'dEdd interior and sfc scattering layers are used in both ice, snow (unequal)'
          write(nu_diag,1020) ' ncat             = ', ncat,  ' : number of ice categories'
          if (kcatbound == 0) then
@@ -1937,19 +1940,24 @@ subroutine input_data
             write(nu_diag,*) '--------------------------------'
             if (trim(shortwave) == 'dEdd') then
                tmpstr2 = ' : delta-Eddington multiple-scattering method'
+            elseif (trim(shortwave) == 'dEdd_snicar_ad') then
+               tmpstr2 = ' : delta-Eddington multiple-scattering method with SNICAR AD'
             elseif (trim(shortwave) == 'ccsm3') then
                tmpstr2 = ' : NCAR CCSM3 distribution method'
             else
                tmpstr2 = ' : unknown value'
             endif
             write(nu_diag,1030) ' shortwave        = ', trim(shortwave),trim(tmpstr2)
-            if (trim(shortwave) == 'dEdd') then
+            if (shortwave(1:4) == 'dEdd') then
                write(nu_diag,1002) ' R_ice            = ', R_ice,' : tuning parameter for sea ice albedo'
                write(nu_diag,1002) ' R_pnd            = ', R_pnd,' : tuning parameter for ponded sea ice albedo'
                write(nu_diag,1002) ' R_snw            = ', R_snw,' : tuning parameter for snow broadband albedo'
                write(nu_diag,1002) ' dT_mlt           = ', dT_mlt,' : change in temperature per change in snow grain radius'
                write(nu_diag,1002) ' rsnw_mlt         = ', rsnw_mlt,' : maximum melting snow grain radius'
                write(nu_diag,1002) ' kalg             = ', kalg,' : absorption coefficient for algae'
+             if (trim(shortwave) == 'dEdd_snicar_ad') then
+               write(nu_diag,1030) ' snw_ssp_table    = ', trim(snw_ssp_table)
+             endif
             else
                if (trim(albedo_type) == 'ccsm3') then
                   tmpstr2 = ' : NCAR CCSM3 albedos'
@@ -2120,7 +2128,7 @@ subroutine input_data
             write(nu_diag,*) 'Using default dEdd melt pond scheme for testing only'
          endif
 
-         if (trim(shortwave) == 'dEdd') then
+         if (shortwave(1:4) == 'dEdd') then
             write(nu_diag,1002) ' hs0              = ', hs0,' : snow depth of transition to bare sea ice'
          endif
 
@@ -2415,7 +2423,7 @@ subroutine input_data
 
       call icepack_init_parameters(ustar_min_in=ustar_min, albicev_in=albicev, albicei_in=albicei, &
          albsnowv_in=albsnowv, albsnowi_in=albsnowi, natmiter_in=natmiter, atmiter_conv_in=atmiter_conv, &
-         emissivity_in=emissivity, &
+         emissivity_in=emissivity, snw_ssp_table_in=snw_ssp_table, &
          ahmax_in=ahmax, shortwave_in=shortwave, albedo_type_in=albedo_type, R_ice_in=R_ice, R_pnd_in=R_pnd, &
          R_snw_in=R_snw, dT_mlt_in=dT_mlt, rsnw_mlt_in=rsnw_mlt, &
          kstrength_in=kstrength, krdg_partic_in=krdg_partic, krdg_redist_in=krdg_redist, mu_rdg_in=mu_rdg, &

cicecore/cicedynB/general/ice_step_mod.F90

@@ -146,8 +146,7 @@ subroutine prep_radiation (iblk)
             alidr_init(i,j,iblk) = alidr_ai(i,j,iblk)
             alidf_init(i,j,iblk) = alidf_ai(i,j,iblk)
 
-            call icepack_prep_radiation (ncat=ncat, nilyr=nilyr, nslyr=nslyr,                 &
-                        scale_factor=scale_factor(i,j,iblk),                                  &
+            call icepack_prep_radiation (scale_factor=scale_factor(i,j,iblk),                 &
                         aice     = aice    (i,j,    iblk), aicen    = aicen   (i,j,  :,iblk), &
                         swvdr    = swvdr   (i,j,    iblk), swvdf    = swvdf   (i,j,    iblk), &
                         swidr    = swidr   (i,j,    iblk), swidf    = swidf   (i,j,    iblk), &
@@ -1223,8 +1222,7 @@ subroutine step_radiation (dt, iblk)
           fswthrun, fswthrun_vdr, fswthrun_vdf, fswthrun_idr, fswthrun_idf, &
           albicen, albsnon, albpndn, &
           alvdrn, alidrn, alvdfn, alidfn, apeffn, trcrn_sw, snowfracn, &
-          kaer_tab, waer_tab, gaer_tab, kaer_bc_tab, waer_bc_tab, &
-          gaer_bc_tab, bcenh, swgrid, igrid
+          swgrid, igrid
       use ice_blocks, only: block, get_block
       use ice_calendar, only: calendar_type, days_per_year, nextsw_cday, yday, msec
       use ice_domain, only: blocks_ice
@@ -1333,9 +1331,7 @@ subroutine step_radiation (dt, iblk)
 
          if (tmask(i,j,iblk)) then
 
-            call icepack_step_radiation (dt=dt,   ncat=ncat,                  &
-                         nblyr=nblyr, nilyr=nilyr, nslyr=nslyr,               &
-                         dEdd_algae=dEdd_algae,                               &
+            call icepack_step_radiation (dt=dt,                               &
                          swgrid=swgrid(:),        igrid=igrid(:),             &
                          fbri=fbri(:),                                        &
                          aicen=aicen(i,j,        :,iblk),                     &
@@ -1355,11 +1351,6 @@ subroutine step_radiation (dt, iblk)
                          days_per_year=days_per_year,                         &
                          nextsw_cday=nextsw_cday, yday=yday,                  &
                          sec=msec,                                             &
-                         kaer_tab=kaer_tab, kaer_bc_tab=kaer_bc_tab(:,:),     &
-                         waer_tab=waer_tab, waer_bc_tab=waer_bc_tab(:,:),     &
-                         gaer_tab=gaer_tab, gaer_bc_tab=gaer_bc_tab(:,:),     &
-                         bcenh=bcenh(:,:,:),                                  &
-                         modal_aero=modal_aero,                               &
                          swvdr    =swvdr    (i,j    ,iblk), swvdf   =swvdf   (i,j    ,iblk), &
                          swidr    =swidr    (i,j    ,iblk), swidf   =swidf   (i,j    ,iblk), &
                          coszen   =coszen   (i,j    ,iblk), fsnow   =fsnow   (i,j    ,iblk), &

cicecore/drivers/direct/hadgem3/CICE_InitMod.F90

@@ -18,7 +18,7 @@ module CICE_InitMod
       use icepack_intfc, only: icepack_aggregate
       use icepack_intfc, only: icepack_init_itd, icepack_init_itd_hist
       use icepack_intfc, only: icepack_init_fsd_bounds, icepack_init_wave
-      use icepack_intfc, only: icepack_configure
+      use icepack_intfc, only: icepack_configure, icepack_init_radiation
       use icepack_intfc, only: icepack_warnings_flush, icepack_warnings_aborted
       use icepack_intfc, only: icepack_query_parameters, icepack_query_tracer_flags, &
           icepack_query_tracer_indices, icepack_query_tracer_sizes
@@ -78,7 +78,7 @@ subroutine cice_init
       use ice_forcing, only: init_forcing_ocn, init_forcing_atmo, &
           get_forcing_atmo, get_forcing_ocn, alloc_forcing, get_wave_spec
       use ice_forcing_bgc, only: get_forcing_bgc, get_atm_bgc, &
-          faero_data, faero_default, faero_optics, alloc_forcing_bgc
+          faero_data, faero_default, alloc_forcing_bgc
       use ice_grid, only: init_grid1, init_grid2, alloc_grid
       use ice_history, only: init_hist, accum_hist
       use ice_restart_shared, only: restart, runid, runtype
@@ -170,15 +170,13 @@ subroutine cice_init
       call init_history_therm   ! initialize thermo history variables
       call init_history_dyn     ! initialize dynamic history variables
       call calc_timesteps       ! update timestep counter if not using npt_unit="1"
+      call icepack_init_radiation ! initialize icepack shortwave tables
 
       call icepack_query_tracer_flags(tr_aero_out=tr_aero, tr_zaero_out=tr_zaero)
       call icepack_warnings_flush(nu_diag)
       if (icepack_warnings_aborted()) call abort_ice(subname, &
           file=__FILE__,line= __LINE__)
 
-      if (tr_aero .or. tr_zaero) call faero_optics !initialize aerosol optical
-                                                   !property tables
-
       ! Initialize shortwave components using swdn from previous timestep
       ! if restarting. These components will be scaled to current forcing
       ! in prep_radiation.

cicecore/drivers/direct/nemo_concepts/CICE_InitMod.F90

@@ -18,7 +18,7 @@ module CICE_InitMod
       use icepack_intfc, only: icepack_aggregate
       use icepack_intfc, only: icepack_init_itd, icepack_init_itd_hist
       use icepack_intfc, only: icepack_init_fsd_bounds, icepack_init_wave
-      use icepack_intfc, only: icepack_configure
+      use icepack_intfc, only: icepack_configure, icepack_init_radiation
       use icepack_intfc, only: icepack_warnings_flush, icepack_warnings_aborted
       use icepack_intfc, only: icepack_query_parameters, icepack_query_tracer_flags, &
           icepack_query_tracer_indices, icepack_query_tracer_sizes
@@ -78,7 +78,7 @@ subroutine cice_init
       use ice_forcing, only: init_forcing_ocn, init_forcing_atmo, &
           get_forcing_atmo, get_forcing_ocn, alloc_forcing, get_wave_spec
       use ice_forcing_bgc, only: get_forcing_bgc, get_atm_bgc, &
-          faero_data, faero_default, faero_optics, alloc_forcing_bgc
+          faero_data, faero_default, alloc_forcing_bgc
       use ice_grid, only: init_grid1, init_grid2, alloc_grid
       use ice_history, only: init_hist, accum_hist
       use ice_restart_shared, only: restart, runid, runtype
@@ -170,15 +170,13 @@ subroutine cice_init
       call init_history_therm   ! initialize thermo history variables
       call init_history_dyn     ! initialize dynamic history variables
       call calc_timesteps       ! update timestep counter if not using npt_unit="1"
+      call icepack_init_radiation ! initialize icepack shortwave tables
 
       call icepack_query_tracer_flags(tr_aero_out=tr_aero, tr_zaero_out=tr_zaero)
       call icepack_warnings_flush(nu_diag)
       if (icepack_warnings_aborted()) call abort_ice(subname, &
           file=__FILE__,line= __LINE__)
 
-      if (tr_aero .or. tr_zaero) call faero_optics !initialize aerosol optical
-                                                   !property tables
-
       ! Initialize shortwave components using swdn from previous timestep
       ! if restarting. These components will be scaled to current forcing
       ! in prep_radiation.