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

Updates to advanced snow physics implementation #449

Merged
merged 5 commits into from
Aug 18, 2023
Merged

Updates to advanced snow physics implementation #449

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
0bccb30
snwgrain bug fix; use snwredist and snwgrain instead of tr_snow in co…
eclare108213 Jul 28, 2023
2e12190
Use snow tracers only for snwredist or snwgrain cases that require th…
eclare108213 Aug 3, 2023
e4e04d4
initialize snow density tracer as rhosnew, as in E3SM
eclare108213 Aug 4, 2023
d6d8d45
do not set rhos_cmpn when the tracer might be undefined in the driver
eclare108213 Aug 4, 2023
8fe89cd
assume all condensation becomes solid ice, for snow mass tracers
eclare108213 Aug 18, 2023
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
6 changes: 3 additions & 3 deletions columnphysics/icepack_flux.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -9,10 +9,10 @@
module icepack_flux

use icepack_kinds
use icepack_parameters, only: c1, emissivity
use icepack_parameters, only: c1, emissivity, snwgrain
use icepack_warnings, only: warnstr, icepack_warnings_add
use icepack_warnings, only: icepack_warnings_setabort, icepack_warnings_aborted
use icepack_tracers, only: tr_iso, tr_snow
use icepack_tracers, only: tr_iso

implicit none
private
Expand Down Expand Up @@ -217,7 +217,7 @@ subroutine merge_fluxes (aicen, &
meltt = meltt + melttn * aicen
meltb = meltb + meltbn * aicen
melts = melts + meltsn * aicen
if (tr_snow) then
if (snwgrain) then
meltsliq = meltsliq + meltsliqn * aicen
endif
if (present(dsnow)) then
Expand Down
14 changes: 9 additions & 5 deletions columnphysics/icepack_itd.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -28,13 +28,13 @@ module icepack_itd
use icepack_kinds
use icepack_parameters, only: c0, c1, c2, c3, c15, c25, c100, p1, p01, p001, p5, puny
use icepack_parameters, only: Lfresh, rhos, ice_ref_salinity, hs_min, cp_ice, Tocnfrz, rhoi
use icepack_parameters, only: rhosi, sk_l, hs_ssl, min_salin, rsnw_fall
use icepack_parameters, only: rhosi, sk_l, hs_ssl, min_salin, rsnw_fall, rhosnew
use icepack_tracers, only: nt_Tsfc, nt_qice, nt_qsno, nt_aero, nt_isosno, nt_isoice
use icepack_tracers, only: nt_apnd, nt_hpnd, nt_fbri, tr_brine, nt_bgc_S, bio_index
use icepack_tracers, only: n_iso, tr_iso, tr_snow, nt_smice, nt_rsnw, nt_rhos, nt_sice
use icepack_tracers, only: n_iso, tr_iso, nt_smice, nt_rsnw, nt_rhos, nt_sice
use icepack_tracers, only: icepack_compute_tracers
use icepack_parameters, only: solve_zsal, skl_bgc, z_tracers
use icepack_parameters, only: kcatbound, kitd, saltflux_option
use icepack_parameters, only: kcatbound, kitd, saltflux_option, snwgrain, snwredist
use icepack_therm_shared, only: Tmin, hi_min
use icepack_warnings, only: warnstr, icepack_warnings_add
use icepack_warnings, only: icepack_warnings_setabort, icepack_warnings_aborted
Expand Down Expand Up @@ -1227,9 +1227,13 @@ subroutine zap_small_areas (dt, ntrcr, &
enddo
endif
if (tr_brine) trcrn(nt_fbri,n) = c1
if (tr_snow) then
if (snwredist(1:3) == 'ITD') then
do k = 1, nslyr
trcrn(nt_rhos +k-1,n) = rhosnew
enddo
endif
if (snwgrain) then
do k = 1, nslyr
trcrn(nt_rhos +k-1,n) = rhos
trcrn(nt_smice+k-1,n) = rhos
trcrn(nt_rsnw +k-1,n) = rsnw_fall
enddo
Expand Down
4 changes: 1 addition & 3 deletions columnphysics/icepack_snow.F90
View file
Open in desktop
Copy link
Contributor

@njeffery njeffery Aug 7, 2023 •
edited by eclare108213
Loading

Choose a reason for hiding this comment

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

@eclare108213 : This expression (L1245 and 1246) forces the ice mass tracer to maintain it's current ice density by adding/subtracting mass from the snow liquid tracer.

1240  if (hsn > puny) then    ! add snow with enthalpy zqsn(1)                                                            
1241        dhs = econ / (zqsn(1) - rhos*Lvap) ! econ < 0, dhs > 0  
1242         mass  = massice(1) + massliq(1)                                                                                  
1243         massi = c0                                                                                                       
1244         if (dzs(1) > puny) massi = c1 + dhs/dzs(1)
1245         massice(1) = massice(1) * massi                                                                                  
1246        massliq(1) = max(c0, mass + rhos*dhs - massice(1)) ! conserve new total mass

It's only going to be round off level changes at this point because the ice mass is always near rho, but it's not what we want. It should just be

massice(1) = massice(1) + dhs * rhos

[edited formatting only]

Copy link
Contributor Author

Choose a reason for hiding this comment

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

To clarify, @njeffery, you are suggesting that we should only change massice(1) and not change massliq(1)? Having thought through this, I think you are right. It's pretty confusing, though, so I'm recording my thought process here for later reference:
Condensation could be done either way, with all of the condensation becoming solid ice or dividing it between ice and liquid. I think the massi code above divides it evenly based on what the tracers were to begin with. But the latent heat values would probably be different in that case. Is that the issue here? This calculation is using Lvaporization, which is between liquid and gas, but if we're condensing from gas directly to ice, shouldn't it be Lsublimation? Answer: No, because the enthalpy definition that we use in the model includes the heat needed to convert between solid and liquid phases (and to bring the liquid to 0 C). Here, we are only changing the ice and liquid mass tracers, which shouldn't affect the solution at this point in the code -- the liquid mass may be used later for melt ponds, which would change the solution in those configurations.

Copy link
Contributor

Choose a reason for hiding this comment

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

@eclare108213 : yes. Just massice. New snow formation has massice = rhos and this is most consistent with both our energy and mass exchanges. The points you're making above about vapor to liquid and vapor to ice are also worth thinking more about, especially if we ever start changing snow density in the thermodynamics. But we okay for conservation. You'll see differences through the melt ponds and the snow grain radius.

Original file line number Diff line number Diff line change
Expand Up @@ -309,14 +309,12 @@ subroutine icepack_step_snow(dt, nilyr, &
! Initialize effective snow density (compaction) for new snow
!-----------------------------------------------------------------

if (trim(snwredist) /= 'none') then
if (snwredist(1:3) == 'ITD') then
do n = 1, ncat
do k = 1, nslyr
if (rhos_cmpn(k,n) < rhosmin) rhos_cmpn(k,n) = rhosnew
enddo
enddo
else
rhos_cmpn(:,:) = rhos
endif

!-----------------------------------------------------------------
Expand Down
8 changes: 4 additions & 4 deletions columnphysics/icepack_therm_itd.F90
View file
Open in desktop
Copy link
Contributor

Choose a reason for hiding this comment

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

@eclare108213 : For the snwredist options, it seems we're only supporting ITDrdg and none. Are bulk, 30percent and ITDsd defunct?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

We're supporting ITDrdg, bulk and none. The old 30percent is available - it's set using 'bulk' with a (new) namelist parameter that can be any percentage. ITDsd is defunct.

Copy link
Contributor

Choose a reason for hiding this comment

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

Good to know. I'll have to look into how the bulk option would work in mpas-seaice. Does 'bulk' use the effective snow density?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

I think bulk works the same way as 30percent did before, just generalized the percentage. So it only rearranges the snow without changing the density. And yes the density tracer is still available in that case - it should remain constant, although I haven't looked at that yet after this latest round of changes.

Copy link
Contributor

Choose a reason for hiding this comment

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

I'm concerned that the way it's implemented, we are requiring the density tracer for the bulk option. Otherwise it will fail. And it doesn't do anything.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

It doesn't fail in Icepack (or CICE). Can MPAS-SI not implement it that way? With it, analysis scripts looking for it (for comparison with ITDrdg) won't fail, and its values can be verified to be constant.
Have you tried running this icepack version in your E3SM tests? I think it's coded so that if the driver doesn't have the tracer, icepack will still handle it gracefully. But I'm not sure about that. If it doesn't work, then I can change it.

Copy link
Contributor

Choose a reason for hiding this comment

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

I haven't tested (recently) snwredist = bulk and snow effective density off, but I will. But tracers are passed to icepack in the tracer array which in mpas is only made up of active tracers. The index for nonactive tracers is set to 0. So in something like

do k = nt_rhos, nt_rhos+nslyr-1
trcrn(k,n) = max(trcrn(k,n)-rhosmin, c0)

I suspect F90 won't like that nt_rhos = 0, and the rest of the code will overwrite something valuable.

I could be wrong. I've almost got your changes implemented. I'll try running with different options.

Original file line number Diff line number Diff line change
Expand Up @@ -24,7 +24,7 @@ module icepack_therm_itd
use icepack_parameters, only: p001, p1, p333, p5, p666, puny, bignum
use icepack_parameters, only: rhos, rhoi, Lfresh, ice_ref_salinity
use icepack_parameters, only: phi_init, dsin0_frazil, hs_ssl, salt_loss
use icepack_parameters, only: rhosi, conserv_check, rhosmin
use icepack_parameters, only: rhosi, conserv_check, rhosmin, snwredist
use icepack_parameters, only: kitd, ktherm
use icepack_parameters, only: z_tracers, solve_zsal, hfrazilmin
use icepack_parameters, only: saltflux_option
Expand All @@ -35,7 +35,7 @@ module icepack_therm_itd
use icepack_tracers, only: nt_apnd, nt_hpnd, nt_aero, nt_isosno, nt_isoice
use icepack_tracers, only: nt_Tsfc, nt_iage, nt_FY, nt_fsd, nt_rhos, nt_sice
use icepack_tracers, only: nt_alvl, nt_vlvl
use icepack_tracers, only: tr_pond_lvl, tr_pond_topo, tr_snow
use icepack_tracers, only: tr_pond_lvl, tr_pond_topo
use icepack_tracers, only: tr_iage, tr_FY, tr_lvl, tr_aero, tr_iso, tr_brine, tr_fsd
use icepack_tracers, only: n_aero, n_iso
use icepack_tracers, only: bio_index
Expand Down Expand Up @@ -569,7 +569,7 @@ subroutine linear_itd (ncat, hin_max, &
enddo
enddo
! maintain rhos_cmp positive definiteness
if (tr_snow) then
if (snwredist(1:3) == 'ITD') then
do n = 1, ncat
do k = nt_rhos, nt_rhos+nslyr-1
trcrn(k,n) = max(trcrn(k,n)-rhosmin, c0)
Expand All @@ -596,7 +596,7 @@ subroutine linear_itd (ncat, hin_max, &
enddo
enddo
! maintain rhos_cmp positive definiteness
if (tr_snow) then
if (snwredist(1:3) == 'ITD') then
do n = 1, ncat
do k = nt_rhos, nt_rhos+nslyr-1
trcrn(k,n) = trcrn(k,n) + rhosmin
Expand Down
2 changes: 1 addition & 1 deletion columnphysics/icepack_therm_mushy.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3307,7 +3307,7 @@ subroutine flood_ice(hsn, hin, &
! density of newly formed snow-ice
rho_snowice = phi_snowice * rho_ocn + (c1 - phi_snowice) * rhoi
endif ! freeboard_density > c0
! endif ! tr_snow
! endif ! snwgrain

if (freeboard_density > c0) then ! ice is flooded

Expand Down
68 changes: 35 additions & 33 deletions columnphysics/icepack_therm_vertical.F90
View file
Open in desktop
Copy link
Contributor

Choose a reason for hiding this comment

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

The changes 141, 267 and 268 would be a problem with other ktherm options and the snow grain tracer active. Since the internal assignments of massice and massliq are wrapped in a ktherm == 2 statement, they would be sent as zeros to thickness_changes.

Copy link
Contributor

Choose a reason for hiding this comment

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

My first look didn't catch this. Probably the simplest fix is to just pull the assignment statements out of the ktherm == 2 if-else and move to after both the temperature_changes routines.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Do we need to support advanced snow physics on BL99 thermo? We could abort if that combination is attempted. Or we could initialize massice and massliq according to the BL99 assumptions (use rhos for massice and 0 for massliq?).

Copy link
Contributor

Choose a reason for hiding this comment

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

I think just moving it out of the if statement is easiest. It has it's own if snwgrain ... wrapper. I'm certainly not in favor of supporting a massive number of config combinations, but the snow stuff doesn't really care about mushy or BL99. That's not to say I've tested it with BL99...

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Do you mean line 333ff?

Copy link
Contributor

Choose a reason for hiding this comment

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

Yes

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Okay, I agree that should be changed.

Original file line number Diff line number Diff line change
Expand Up @@ -136,7 +136,9 @@ subroutine thermo_vertical (nilyr, nslyr, &
zSin , & ! internal ice layer salinities
rsnw , & ! snow grain radius (10^-6 m)
smice , & ! ice mass tracer in snow (kg/m^3)
smliq , & ! liquid water mass tracer in snow (kg/m^3)
smliq ! liquid water mass tracer in snow (kg/m^3)

real (kind=dbl_kind), dimension (:), intent(out) :: &
massice , & ! ice mass in snow (kg/m^2)
massliq ! liquid water mass in snow (kg/m^2)

Expand Down Expand Up @@ -262,6 +264,8 @@ subroutine thermo_vertical (nilyr, nslyr, &
zTsn(:) = c0
zTin(:) = c0
meltsliq= c0
massice(:) = c0
massliq(:) = c0

if (calc_Tsfc) then
fsensn = c0
Expand Down Expand Up @@ -325,12 +329,6 @@ subroutine thermo_vertical (nilyr, nslyr, &
smice, smliq)
if (icepack_warnings_aborted(subname)) return

! reinitialize mass in case of snow-ice formation
if (snwgrain) then
massice(:) = smice(:) * hslyr
massliq(:) = smliq(:) * hslyr
endif

else ! ktherm

call temperature_changes(dt, &
Expand All @@ -353,6 +351,12 @@ subroutine thermo_vertical (nilyr, nslyr, &

endif ! ktherm

! mass of ice and liquid water in snow
if (snwgrain) then
massice(:) = smice(:) * hslyr
massliq(:) = smliq(:) * hslyr
endif

! intermediate energy for error check

einter = c0
Expand Down Expand Up @@ -1235,11 +1239,8 @@ subroutine thickness_changes (nilyr, nslyr, &
if (hsn > puny) then ! add snow with enthalpy zqsn(1)
dhs = econ / (zqsn(1) - rhos*Lvap) ! econ < 0, dhs > 0

mass = massice(1) + massliq(1)
massi = c0
if (dzs(1) > puny) massi = c1 + dhs/dzs(1)
massice(1) = massice(1) * massi
massliq(1) = max(c0, mass + rhos*dhs - massice(1)) ! conserve new total mass
! assume all condensation becomes ice (no liquid)
massice(1) = massice(1) + dhs*rhos

dzs(1) = dzs(1) + dhs
evapn = evapn + dhs*rhos
Expand Down Expand Up @@ -2357,8 +2358,9 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, &
n ! category index

real (kind=dbl_kind) :: &
rnslyr , & ! 1 / nslyr
worka, workb ! temporary variables
worka , & ! temporary variables
workb , &
workc

! 2D coupler variables (computed for each category, then aggregated)
real (kind=dbl_kind) :: &
Expand Down Expand Up @@ -2465,6 +2467,10 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, &
l_meltsliq = c0
l_meltsliqn = c0

! solid and liquid components of snow mass
massicen(:,:) = c0
massliqn(:,:) = c0

!-----------------------------------------------------------------
! Initialize rate of snow loss to leads
!-----------------------------------------------------------------
Expand All @@ -2489,21 +2495,21 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, &
fsnow = fsnow*(c1-worka)
endif ! snwredist

!echmod - remove all of this code - non-BFB because of values carried in tracer arrays when the snow physics options are not active
!-----------------------------------------------------------------
! solid and liquid components of snow mass
!-----------------------------------------------------------------

massicen(:,:) = c0
massliqn(:,:) = c0
if (snwgrain) then
rnslyr = c1 / real(nslyr, dbl_kind)
do n = 1, ncat
do k = 1, nslyr
massicen(k,n) = smicen(k,n) * vsnon(n) * rnslyr ! kg/m^2
massliqn(k,n) = smliqn(k,n) * vsnon(n) * rnslyr
enddo
enddo
endif
! massicen(:,:) = c0
! massliqn(:,:) = c0
! if (snwgrain) then
! rnslyr = c1 / real(nslyr, dbl_kind)
! do n = 1, ncat
! do k = 1, nslyr
! massicen(k,n) = smicen(k,n) * vsnon(n) * rnslyr ! kg/m^2
! massliqn(k,n) = smliqn(k,n) * vsnon(n) * rnslyr
! enddo
! enddo
! endif

!-----------------------------------------------------------------
! Update the neutral drag coefficients to account for form drag
Expand Down Expand Up @@ -2904,14 +2910,10 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, &
if (snwgrain) then
do n = 1, ncat
if (vsnon(n) > puny) then
workc = real(nslyr, dbl_kind) * aicen(n) / vsnon(n)
do k = 1, nslyr
smicen(k,n) = massicen(k,n) / (vsnon(n) * rnslyr)
smliqn(k,n) = massliqn(k,n) / (vsnon(n) * rnslyr)
worka = smicen(k,n) + smliqn(k,n)
if (worka > puny) then
smicen(k,n) = rhos * smicen(k,n) / worka
smliqn(k,n) = rhos * smliqn(k,n) / worka
endif
smicen(k,n) = massicen(k,n) * workc
smliqn(k,n) = massliqn(k,n) * workc
enddo
else ! reset to default values
do k = 1, nslyr
Expand Down
23 changes: 14 additions & 9 deletions columnphysics/icepack_tracers.F90
View file
Open in desktop
Copy link
Contributor

Choose a reason for hiding this comment

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

@eclare108213 : line 1294 in icepack_tracers.F90 should also be changed to rhosnew.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

done

Copy link
Contributor

Choose a reason for hiding this comment

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

This else part of this statement in icepack_snow.F90 is a problem. rhos_cmpn can't be used at all if snow effective density is false otherwise it will overwrite ice area.
Lines 318 to 320
else
rhos_cmpn(:,:) = rhos
endif

Copy link
Contributor

Choose a reason for hiding this comment

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

And the first part of the if statement
if (trim(snwredist) /= 'none') then
should be changed to
if(snwredist(1:3) = 'ITD') then

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Agreed. Will it work to just remove the else part? This is part of the reason for defining the all of the tracers when tr_snow=T.

Copy link
Contributor

Choose a reason for hiding this comment

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

yes, deleting L318 and 319 works

Original file line number Diff line number Diff line change
Expand Up @@ -7,7 +7,8 @@
module icepack_tracers

use icepack_kinds
use icepack_parameters, only: c0, c1, puny, Tocnfrz, rhos, rsnw_fall
use icepack_parameters, only: c0, c1, puny, Tocnfrz, rhos, rhosnew, rsnw_fall
use icepack_parameters, only: snwredist, snwgrain
use icepack_warnings, only: warnstr, icepack_warnings_add
use icepack_warnings, only: icepack_warnings_setabort, icepack_warnings_aborted

Expand Down Expand Up @@ -105,7 +106,7 @@ module icepack_tracers
tr_pond = .false., & ! if .true., use melt pond tracer
tr_pond_lvl = .false., & ! if .true., use level-ice pond tracer
tr_pond_topo = .false., & ! if .true., use explicit topography-based ponds
tr_snow = .false., & ! if .true., use snow metamorphosis tracers
tr_snow = .false., & ! if .true., use snow redistribution or metamorphosis tracers
tr_iso = .false., & ! if .true., use isotope tracers
tr_aero = .false., & ! if .true., use aerosol tracers
tr_brine = .false., & ! if .true., brine height differs from ice thickness
Expand Down Expand Up @@ -219,7 +220,7 @@ subroutine icepack_init_tracer_flags(&
tr_pond_in , & ! if .true., use melt pond tracer
tr_pond_lvl_in , & ! if .true., use level-ice pond tracer
tr_pond_topo_in , & ! if .true., use explicit topography-based ponds
tr_snow_in , & ! if .true., use snow metamorphosis tracers
tr_snow_in , & ! if .true., use snow redistribution or metamorphosis tracers
tr_fsd_in , & ! if .true., use floe size distribution tracers
tr_iso_in , & ! if .true., use isotope tracers
tr_aero_in , & ! if .true., use aerosol tracers
Expand Down Expand Up @@ -286,7 +287,7 @@ subroutine icepack_query_tracer_flags(&
tr_pond_out , & ! if .true., use melt pond tracer
tr_pond_lvl_out , & ! if .true., use level-ice pond tracer
tr_pond_topo_out , & ! if .true., use explicit topography-based ponds
tr_snow_out , & ! if .true., use snow metamorphosis tracers
tr_snow_out , & ! if .true., use snow redistribution or metamorphosis tracers
tr_fsd_out , & ! if .true., use floe size distribution
tr_iso_out , & ! if .true., use isotope tracers
tr_aero_out , & ! if .true., use aerosol tracers
Expand Down Expand Up @@ -1288,11 +1289,15 @@ subroutine icepack_compute_tracers (ntrcr, trcr_depend, &
enddo

if (vicen <= c0 .and. tr_brine) trcrn(nt_fbri) = c1
if (vsnon <= c0 .and. tr_snow) then
trcrn(nt_rsnw :nt_rsnw +nslyr-1) = rsnw_fall
trcrn(nt_smice:nt_smice+nslyr-1) = rhos
trcrn(nt_rhos :nt_rhos +nslyr-1) = rhos
endif
if (vsnon <= c0) then
if (snwredist(1:3) == 'ITD') then
trcrn(nt_rhos :nt_rhos +nslyr-1) = rhosnew
endif
if (snwgrain) then
trcrn(nt_rsnw :nt_rsnw +nslyr-1) = rsnw_fall
trcrn(nt_smice:nt_smice+nslyr-1) = rhos
endif
endif ! vsnon <= 0

end subroutine icepack_compute_tracers

Expand Down
9 changes: 7 additions & 2 deletions configuration/driver/icedrv_init.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -504,7 +504,7 @@ subroutine input_data
shortwave = 'dEdd'
endif

if (snwredist(1:4) /= 'none' .and. .not. tr_snow) then
if (snwredist(1:3) == 'ITD' .and. .not. tr_snow) then
write (nu_diag,*) 'WARNING: snwredist on but tr_snow=F'
call icedrv_system_abort(file=__FILE__,line=__LINE__)
endif
Expand Down Expand Up @@ -567,11 +567,16 @@ subroutine input_data
shortwave = 'dEdd'
endif

if (tr_snow .and. trim(shortwave) /= 'dEdd') then
if (snwgrain .and. trim(shortwave) /= 'dEdd') then
write (nu_diag,*) 'WARNING: snow grain radius activated but'
write (nu_diag,*) 'WARNING: dEdd shortwave is not.'
endif

if (snwredist(1:4) /= 'none' .and. trim(shortwave) /= 'dEdd') then
write (nu_diag,*) 'WARNING: snow redistribution activated but'
write (nu_diag,*) 'WARNING: dEdd shortwave is not.'
endif

rfracmin = min(max(rfracmin,c0),c1)
rfracmax = min(max(rfracmax,c0),c1)

Expand Down
13 changes: 7 additions & 6 deletions configuration/driver/icedrv_step.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -155,7 +155,7 @@ subroutine step_therm1 (dt)
nt_aero, nt_isosno, nt_isoice, nt_rsnw, nt_smice, nt_smliq

logical (kind=log_kind) :: &
tr_iage, tr_FY, tr_aero, tr_iso, calc_Tsfc, tr_snow
tr_iage, tr_FY, tr_aero, tr_iso, calc_Tsfc, snwgrain

real (kind=dbl_kind), dimension(n_aero,2,ncat) :: &
aerosno, aeroice ! kg/m^2
Expand All @@ -177,6 +177,7 @@ subroutine step_therm1 (dt)

call icepack_query_parameters(puny_out=puny)
call icepack_query_parameters(calc_Tsfc_out=calc_Tsfc)
call icepack_query_parameters(snwgrain_out=snwgrain)
call icepack_warnings_flush(nu_diag)
if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, &
file=__FILE__,line= __LINE__)
Expand All @@ -189,7 +190,7 @@ subroutine step_therm1 (dt)

call icepack_query_tracer_flags( &
tr_iage_out=tr_iage, tr_FY_out=tr_FY, &
tr_aero_out=tr_aero, tr_iso_out=tr_iso, tr_snow_out=tr_snow)
tr_aero_out=tr_aero, tr_iso_out=tr_iso)
call icepack_warnings_flush(nu_diag)
if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, &
file=__FILE__,line= __LINE__)
Expand Down Expand Up @@ -261,15 +262,15 @@ subroutine step_therm1 (dt)
enddo
endif ! tr_iso

if (tr_snow) then
if (snwgrain) then
do n = 1, ncat
do k = 1, nslyr
rsnwn (k,n) = trcrn(i,nt_rsnw +k-1,n)
smicen(k,n) = trcrn(i,nt_smice+k-1,n)
smliqn(k,n) = trcrn(i,nt_smliq+k-1,n)
enddo
enddo
endif ! tr_snow
endif ! snwgrain

call icepack_step_therm1(dt=dt, ncat=ncat, nilyr=nilyr, nslyr=nslyr, &
aicen_init = aicen_init(i,:), &
Expand Down Expand Up @@ -400,15 +401,15 @@ subroutine step_therm1 (dt)
enddo
endif ! tr_iso

if (tr_snow) then
if (snwgrain) then
do n = 1, ncat
do k = 1, nslyr
trcrn(i,nt_rsnw +k-1,n) = rsnwn (k,n)
trcrn(i,nt_smice+k-1,n) = smicen(k,n)
trcrn(i,nt_smliq+k-1,n) = smliqn(k,n)
enddo
enddo
endif ! tr_snow
endif ! snwgrain

enddo ! i
call icepack_warnings_flush(nu_diag)
Expand Down
Loading