Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Update implementation of g(k), w0(k) in icepack_shortwave #13

Merged
merged 3 commits into from
Sep 18, 2022
Merged

Update implementation of g(k), w0(k) in icepack_shortwave #13

Show file tree
Hide file tree
Changes from 2 commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
44d377b
Update implementation of g(k), w0(k) in icepack_shortwave
apcraig Sep 11, 2022
6f8f21b
Get rid of use_snicar flag, redundant with shortwave setting
apcraig Sep 12, 2022
9178ede
cleanup
apcraig Sep 18, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
19 changes: 6 additions & 13 deletions columnphysics/icepack_parameters.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand All @@ -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'

Expand Down Expand Up @@ -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, &
Expand Down Expand Up @@ -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

Expand All @@ -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 :: &
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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, &
Expand Down Expand Up @@ -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

Expand All @@ -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 :: &
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down
102 changes: 45 additions & 57 deletions columnphysics/icepack_shortwave.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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

Expand Down Expand Up @@ -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, &
Expand Down Expand Up @@ -2362,12 +2362,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))
Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

is it possible for the denominators to be zero?

Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@njeffery
We are reducing the calculations here (and several other places) as was done in MPAS-SI -- it's much cleaner. The original has a bunch of puny's in the (unnecessary) denominators. My question is whether it's possible for the denominators to ever be zero in the new code, as long as tr_zaero is turned on?

tau(k) = tau(k) + tzaer(ns,k)
enddo
elseif (tr_aero) then
Expand Down Expand Up @@ -2413,6 +2411,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
Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Good question! It's possible that this code is left over from my trying to combine the 3-band and 5-band routines. If this module is diff'ed with the original icepack_shortwave.F90 in the latest release, does one of the duplicate code sections show up as an addition?

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The same coding does seem to also be in the Icepack main. See https://github.com/CICE-Consortium/Icepack/blob/main/columnphysics/icepack_shortwave.F90#L2657 which has exactly the same situation. So I think it's also in the current code. One thing to remember is that I don't think we've been testing the modal aerosols at all because they don't work in CICE and the Icepack driver doesn't seem to have all the required options to do those tests either. So we have no baseline to test against and haven't been exercising that part of the code.

Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It also appears in MPAS-SI in both compute_dEdd (3-band) and compute_dEdd_5bd: https://github.com/E3SM-Project/E3SM/blob/851fbfbbbdc13915de772519944bc44b4644dfb7/components/mpas-seaice/src/column/ice_shortwave.F90#L4854
@njeffery do you know what's going on here?

Copy link
Owner

@eclare108213 eclare108213 Sep 12, 2022 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Apparently this issue entered the code in 2015, in https://github.com/E3SM-Project/CICE-archive/commit/6e8e276fce7eb35338d21d708397cad62ae59d1d# (or maybe another branch feeding this one). !mgf seems to be Mark Flanner; Hailong Wang is also referenced in the commit message -- this code seems to have come via HiLAT.
There is a bit of code that was dropped, which used taer etc to calculate g(k), w0(k) and tau(k) for k=0 (lines 2199ff), before taer etc are reset to 0. k=0 is the surface scattering layer of snow. [Much of the rest of the compute_dEdd routine calculates g, w0 and tau for the other layers in both snow and ice, which are then used in the solution_dEdd routine.] The missing code is only for tr_aero = true (not tr_zaero). The g, w0 and tau vectors are computed near the end of the tr_aero conditional block, but only for k >= 1, never for k=0.
This looks like a bug to me. k=0 values appear to be initialized near the top of the "begin spectral loop", regardless of whether aerosol tracers are turned on, which might mask the impacts of a bug.

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I agree, it looks like setting of g(k), w0(k), and tau(k) is missing after that code block. But I don't really know for sure. As we move forward, we should take extra time to add test cases and validate the varied combination of options associated with aerosols, snicar, and bgc. Especially as we add an updated implementation of bgc.


do k=1,nslyr
taer = c0
waer = c0
Expand Down Expand Up @@ -2462,12 +2462,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
Expand Down Expand Up @@ -2523,12 +2521,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
Expand Down Expand Up @@ -2579,12 +2575,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
Expand Down Expand Up @@ -2632,12 +2626,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
Expand Down Expand Up @@ -4969,9 +4961,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
Expand Down Expand Up @@ -5030,6 +5022,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
Expand Down Expand Up @@ -5083,12 +5077,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
Expand Down Expand Up @@ -5134,9 +5126,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
Expand Down Expand Up @@ -5192,12 +5184,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
Expand Down Expand Up @@ -5252,12 +5242,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
Expand Down
2 changes: 1 addition & 1 deletion columnphysics/icepack_shortwave_data.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -5,7 +5,7 @@ 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_parameters, only: snw_ssp_table
use icepack_tracers, only: nmodal1, nmodal2, max_aero

implicit none
Expand Down
2 changes: 1 addition & 1 deletion configuration/driver/icedrv_init.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down
6 changes: 2 additions & 4 deletions configuration/driver/icedrv_init_column.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down Expand Up @@ -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__)
Expand Down
2 changes: 1 addition & 1 deletion configuration/scripts/options/set_nml.snicar
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,2 +1,2 @@
shortwave = 'dEdd_snicar'
shortwave = 'dEdd_snicar_ad'
snw_ssp_table = 'snicar'
2 changes: 1 addition & 1 deletion configuration/scripts/options/set_nml.snicartest
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,2 +1,2 @@
shortwave = 'dEdd_snicar'
shortwave = 'dEdd_snicar_ad'
snw_ssp_table = 'test'