Skip to content

Commit

Permalink
Merge branch 'fix_issue965'
Browse files Browse the repository at this point in the history 
Fix bug in snow aerosol numerics (causes occasional HUGE aerosol values)

Multiplications / divisions by dtime were missing in a few places, in
some code that ensures that aerosol mass fluxes between the snow layers
don't exceed the total mass. This bug was introduced in ctsm1.0.dev065,
when a term was changed from a mass to a flux (per unit time), but some
of the code wasn't changed correspondingly.

In one case this bug caused a crash, but in other cases it just causes
incorrect answers. The issue can be seen by examining the history
variables SNOBCMCL, SNOOCMCL and SNODSTMCL. In the case that led to the
crash, SNOBCMCL and SNOOCMCL became huge in one grid cell (~ 10^12) for
a long period before the crash.

Resolves ESCOMP#965
  • Loading branch information
@billsacks
billsacks committed Apr 7, 2020
2 parents 9e25889 + 16e4f58 commit 0dea14e
Show file tree
Hide file tree
Showing 3 changed files with 213 additions and 48 deletions.
148 changes: 148 additions & 0 deletions doc/ChangeLog
View file
Original file line number Diff line number Diff line change
@@ -1,4 +1,152 @@
===============================================================
Tag name: ctsm1.0.dev088
Originator(s): sacks (Bill Sacks)
Date: Tue Apr 7 10:43:03 MDT 2020
One-line Summary: Fix bug in snow aerosol numerics (causes occasional HUGE aerosol values)

Purpose of changes
------------------

Multiplications / divisions by dtime were missing in a few places, in
some code that ensures that aerosol mass fluxes between the snow layers
don't exceed the total mass. This bug was introduced in ctsm1.0.dev065,
when a term was changed from a mass to a flux (per unit time), but some
of the code wasn't changed correspondingly.

In one case this bug caused a crash, but in other cases it just causes
incorrect answers. The issue can be seen by examining the history
variables SNOBCMCL, SNOOCMCL and SNODSTMCL. In the case that led to the
crash, SNOBCMCL and SNOOCMCL became huge in one grid cell (~ 10^12) for
a long period before the crash.

Bugs fixed or introduced
------------------------

Issues fixed (include CTSM Issue #):
- Resolves ESCOMP/ctsm#965 (Bug in snow aerosol numerics can cause
aerosol masses to blow up, starting with ctsm1.0.dev065)

Significant changes to scientifically-supported configurations
--------------------------------------------------------------

Does this tag change answers significantly for any of the following physics configurations?
(Details of any changes will be given in the "Answer changes" section below.)

[Put an [X] in the box for any configuration with significant answer changes.]

[ ] clm5_0

[ ] ctsm5_0-nwp

[ ] clm4_5

Notes of particular relevance for users
---------------------------------------

Caveats for users (e.g., need to interpolate initial conditions): none

Changes to CTSM's user interface (e.g., new/renamed XML or namelist variables): none

Changes made to namelist defaults (e.g., changed parameter values): none

Changes to the datasets (e.g., parameter, surface or initial files): none

Substantial timing or memory changes: none (not checked, but none expected)

Notes of particular relevance for developers: (including Code reviews and testing)
---------------------------------------------
NOTE: Be sure to review the steps in README.CHECKLIST.master_tags as well as the coding style in the Developers Guide

Caveats for developers (e.g., code that is duplicated that requires double maintenance): none

Changes to tests or testing: none

Code reviewed by: Erik Kluzek, Keith Oleson


CTSM testing:

[PASS means all tests PASS and OK means tests PASS other than expected fails.]

build-namelist tests:

cheyenne - not run

tools-tests (test/tools):

cheyenne - not run

PTCLM testing (tools/shared/PTCLM/test):

cheyenne - not run

python testing (see instructions in python/README.md; document testing done):

(any machine) - not run

regular tests (aux_clm):

cheyenne ---- ok
izumi ------- pass

ok means tests pass, some answer changes as expected

If the tag used for baseline comparisons was NOT the previous tag, note that here:


Answer changes
--------------

Changes answers relative to baseline: YES

If a tag changes answers relative to baseline comparison the
following should be filled in (otherwise remove this section):

Summarize any changes to answers, i.e.,
- what code configurations: all
- what platforms/compilers: all
- nature of change (roundoff; larger than roundoff/same climate; new climate):
larger than roundoff/same climate

Answer changes arise due to a bug fix for a condition that only
appears infrequently. When it does appear, answer changes can
potentially be enormous (e.g., snow aerosol concentrations 20
orders of magnitude too large), though answer changes observed in
the test suite all appear relatively small.

Answer changes just appeared for the following tests:

FAIL ERP_D_P36x2_Ld3.f10_f10_musgs.I2000Clm45BgcCrop.cheyenne_gnu.clm-no_subgrid_fluxes BASELINE ctsm1.0.dev087: DIFF
FAIL SMS_Ld1.nldas2_rnldas2_mnldas2.I2000Ctsm50NwpSpNldasGs.cheyenne_gnu.clm-default BASELINE ctsm1.0.dev087: DIFF
FAIL SMS_Ld1.nldas2_rnldas2_mnldas2.I2000Ctsm50NwpSpNldasRsGs.cheyenne_gnu.clm-default BASELINE ctsm1.0.dev087: DIFF
FAIL ERI_D_Ld9.ne30_g16.I2000Clm50BgcCruGs.cheyenne_intel.clm-vrtlay BASELINE ctsm1.0.dev087: DIFF
FAIL ERI_D_Ld9.T31_g37.I2000Clm50Sp.cheyenne_intel.clm-SNICARFRC BASELINE ctsm1.0.dev087: DIFF
FAIL ERI_Ld9.f45_g37.I2000Clm50BgcCruGs.cheyenne_intel.clm-nofire BASELINE ctsm1.0.dev087: DIFF
FAIL ERI_N2_Ld9.f19_g17.I2000Clm50BgcCrop.cheyenne_intel.clm-default BASELINE ctsm1.0.dev087: DIFF
FAIL SMS_D_Lm6.f45_f45_mg37.I2000Clm50Fates.cheyenne_intel.clm-FatesColdDef BASELINE ctsm1.0.dev087: DIFF
FAIL SMS_D_Lm6_P144x1.f45_f45_mg37.I2000Clm50Fates.cheyenne_intel.clm-FatesColdDef BASELINE ctsm1.0.dev087: DIFF
FAIL SMS_Lm13.f19_g17.I2000Clm50BgcCrop.cheyenne_intel.clm-cropMonthOutput BASELINE ctsm1.0.dev087: DIFF

If bitwise differences were observed, how did you show they were no worse
than roundoff? N/A

If this tag changes climate describe the run(s) done to evaluate the new
climate (put details of the simulations in the experiment database)
- casename: N/A

URL for LMWG diagnostics output used to validate new climate: N/A


Detailed list of changes
------------------------

List any externals directories updated (cime, rtm, mosart, cism, fates, etc.): none

Pull Requests that document the changes (include PR ids):
https://github.com/ESCOMP/CTSM/pull/971

===============================================================
===============================================================
Tag name: ctsm1.0.dev087
Originator(s): sacks (Bill Sacks)
Date: Mon Apr 6 11:23:51 MDT 2020
Expand Down
1 change: 1 addition & 0 deletions doc/ChangeSum
View file
Original file line number Diff line number Diff line change
@@ -1,5 +1,6 @@
Tag Who Date Summary
============================================================================================================================
ctsm1.0.dev088 sacks 04/07/2020 Fix bug in snow aerosol numerics (causes occasional HUGE aerosol values)
ctsm1.0.dev087 sacks 04/06/2020 Change hard-coded soil layers in phenology to use a target depth
ctsm1.0.dev086 mvertens 03/19/2020 Updates to NUOPC cap and minor fixes to ncdio_pio
ctsm1.0.dev085 sacks 03/16/2020 Update cime
Expand Down
112 changes: 64 additions & 48 deletions src/biogeophys/SnowHydrologyMod.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1438,22 +1438,22 @@ subroutine CalcAndApplyAerosolFluxes(bounds, num_snowc, filter_snowc, &
integer :: fc, c
integer :: j
real(r8) :: mss_liqice ! mass of liquid+ice in a layer
real(r8) :: qin_bc_phi (bounds%begc:bounds%endc) ! flux of hydrophilic BC into layer [kg]
real(r8) :: qout_bc_phi (bounds%begc:bounds%endc) ! flux of hydrophilic BC out of layer [kg]
real(r8) :: qin_bc_pho (bounds%begc:bounds%endc) ! flux of hydrophobic BC into layer [kg]
real(r8) :: qout_bc_pho (bounds%begc:bounds%endc) ! flux of hydrophobic BC out of layer [kg]
real(r8) :: qin_oc_phi (bounds%begc:bounds%endc) ! flux of hydrophilic OC into layer [kg]
real(r8) :: qout_oc_phi (bounds%begc:bounds%endc) ! flux of hydrophilic OC out of layer [kg]
real(r8) :: qin_oc_pho (bounds%begc:bounds%endc) ! flux of hydrophobic OC into layer [kg]
real(r8) :: qout_oc_pho (bounds%begc:bounds%endc) ! flux of hydrophobic OC out of layer [kg]
real(r8) :: qin_dst1 (bounds%begc:bounds%endc) ! flux of dust species 1 into layer [kg]
real(r8) :: qout_dst1 (bounds%begc:bounds%endc) ! flux of dust species 1 out of layer [kg]
real(r8) :: qin_dst2 (bounds%begc:bounds%endc) ! flux of dust species 2 into layer [kg]
real(r8) :: qout_dst2 (bounds%begc:bounds%endc) ! flux of dust species 2 out of layer [kg]
real(r8) :: qin_dst3 (bounds%begc:bounds%endc) ! flux of dust species 3 into layer [kg]
real(r8) :: qout_dst3 (bounds%begc:bounds%endc) ! flux of dust species 3 out of layer [kg]
real(r8) :: qin_dst4 (bounds%begc:bounds%endc) ! flux of dust species 4 into layer [kg]
real(r8) :: qout_dst4 (bounds%begc:bounds%endc) ! flux of dust species 4 out of layer [kg]
real(r8) :: qin_bc_phi (bounds%begc:bounds%endc) ! flux of hydrophilic BC into layer [kg/s]
real(r8) :: qout_bc_phi (bounds%begc:bounds%endc) ! flux of hydrophilic BC out of layer [kg/s]
real(r8) :: qin_bc_pho (bounds%begc:bounds%endc) ! flux of hydrophobic BC into layer [kg/s]
real(r8) :: qout_bc_pho (bounds%begc:bounds%endc) ! flux of hydrophobic BC out of layer [kg/s]
real(r8) :: qin_oc_phi (bounds%begc:bounds%endc) ! flux of hydrophilic OC into layer [kg/s]
real(r8) :: qout_oc_phi (bounds%begc:bounds%endc) ! flux of hydrophilic OC out of layer [kg/s]
real(r8) :: qin_oc_pho (bounds%begc:bounds%endc) ! flux of hydrophobic OC into layer [kg/s]
real(r8) :: qout_oc_pho (bounds%begc:bounds%endc) ! flux of hydrophobic OC out of layer [kg/s]
real(r8) :: qin_dst1 (bounds%begc:bounds%endc) ! flux of dust species 1 into layer [kg/s]
real(r8) :: qout_dst1 (bounds%begc:bounds%endc) ! flux of dust species 1 out of layer [kg/s]
real(r8) :: qin_dst2 (bounds%begc:bounds%endc) ! flux of dust species 2 into layer [kg/s]
real(r8) :: qout_dst2 (bounds%begc:bounds%endc) ! flux of dust species 2 out of layer [kg/s]
real(r8) :: qin_dst3 (bounds%begc:bounds%endc) ! flux of dust species 3 into layer [kg/s]
real(r8) :: qout_dst3 (bounds%begc:bounds%endc) ! flux of dust species 3 out of layer [kg/s]
real(r8) :: qin_dst4 (bounds%begc:bounds%endc) ! flux of dust species 4 into layer [kg/s]
real(r8) :: qout_dst4 (bounds%begc:bounds%endc) ! flux of dust species 4 out of layer [kg/s]

character(len=*), parameter :: subname = 'CalcAndApplyAerosolFluxes'
!-----------------------------------------------------------------------
Expand Down Expand Up @@ -1519,73 +1519,89 @@ subroutine CalcAndApplyAerosolFluxes(bounds, num_snowc, filter_snowc, &
! BCPHI:
! 1. flux with meltwater:
qout_bc_phi(c) = qflx_snow_percolation(c,j)*scvng_fct_mlt_bcphi*(mss_bcphi(c,j)/mss_liqice)
if (qout_bc_phi(c) > mss_bcphi(c,j)) then
qout_bc_phi(c) = mss_bcphi(c,j)
endif
mss_bcphi(c,j) = mss_bcphi(c,j) - qout_bc_phi(c) * dtime
if (qout_bc_phi(c)*dtime > mss_bcphi(c,j)) then
qout_bc_phi(c) = mss_bcphi(c,j)/dtime
mss_bcphi(c,j) = 0._r8
else
mss_bcphi(c,j) = mss_bcphi(c,j) - qout_bc_phi(c) * dtime
end if
qin_bc_phi(c) = qout_bc_phi(c)

! BCPHO:
! 1. flux with meltwater:
qout_bc_pho(c) = qflx_snow_percolation(c,j)*scvng_fct_mlt_bcpho*(mss_bcpho(c,j)/mss_liqice)
if (qout_bc_pho(c) > mss_bcpho(c,j)) then
qout_bc_pho(c) = mss_bcpho(c,j)
endif
mss_bcpho(c,j) = mss_bcpho(c,j) - qout_bc_pho(c) * dtime
if (qout_bc_pho(c)*dtime > mss_bcpho(c,j)) then
qout_bc_pho(c) = mss_bcpho(c,j)/dtime
mss_bcpho(c,j) = 0._r8
else
mss_bcpho(c,j) = mss_bcpho(c,j) - qout_bc_pho(c) * dtime
end if
qin_bc_pho(c) = qout_bc_pho(c)

! OCPHI:
! 1. flux with meltwater:
qout_oc_phi(c) = qflx_snow_percolation(c,j)*scvng_fct_mlt_ocphi*(mss_ocphi(c,j)/mss_liqice)
if (qout_oc_phi(c) > mss_ocphi(c,j)) then
qout_oc_phi(c) = mss_ocphi(c,j)
endif
mss_ocphi(c,j) = mss_ocphi(c,j) - qout_oc_phi(c) * dtime
if (qout_oc_phi(c)*dtime > mss_ocphi(c,j)) then
qout_oc_phi(c) = mss_ocphi(c,j)/dtime
mss_ocphi(c,j) = 0._r8
else
mss_ocphi(c,j) = mss_ocphi(c,j) - qout_oc_phi(c) * dtime
end if
qin_oc_phi(c) = qout_oc_phi(c)

! OCPHO:
! 1. flux with meltwater:
qout_oc_pho(c) = qflx_snow_percolation(c,j)*scvng_fct_mlt_ocpho*(mss_ocpho(c,j)/mss_liqice)
if (qout_oc_pho(c) > mss_ocpho(c,j)) then
qout_oc_pho(c) = mss_ocpho(c,j)
endif
mss_ocpho(c,j) = mss_ocpho(c,j) - qout_oc_pho(c) * dtime
if (qout_oc_pho(c)*dtime > mss_ocpho(c,j)) then
qout_oc_pho(c) = mss_ocpho(c,j)/dtime
mss_ocpho(c,j) = 0._r8
else
mss_ocpho(c,j) = mss_ocpho(c,j) - qout_oc_pho(c) * dtime
end if
qin_oc_pho(c) = qout_oc_pho(c)

! DUST 1:
! 1. flux with meltwater:
qout_dst1(c) = qflx_snow_percolation(c,j)*scvng_fct_mlt_dst1*(mss_dst1(c,j)/mss_liqice)
if (qout_dst1(c) > mss_dst1(c,j)) then
qout_dst1(c) = mss_dst1(c,j)
endif
mss_dst1(c,j) = mss_dst1(c,j) - qout_dst1(c) * dtime
if (qout_dst1(c)*dtime > mss_dst1(c,j)) then
qout_dst1(c) = mss_dst1(c,j)/dtime
mss_dst1(c,j) = 0._r8
else
mss_dst1(c,j) = mss_dst1(c,j) - qout_dst1(c) * dtime
end if
qin_dst1(c) = qout_dst1(c)

! DUST 2:
! 1. flux with meltwater:
qout_dst2(c) = qflx_snow_percolation(c,j)*scvng_fct_mlt_dst2*(mss_dst2(c,j)/mss_liqice)
if (qout_dst2(c) > mss_dst2(c,j)) then
qout_dst2(c) = mss_dst2(c,j)
endif
mss_dst2(c,j) = mss_dst2(c,j) - qout_dst2(c) * dtime
if (qout_dst2(c)*dtime > mss_dst2(c,j)) then
qout_dst2(c) = mss_dst2(c,j)/dtime
mss_dst2(c,j) = 0._r8
else
mss_dst2(c,j) = mss_dst2(c,j) - qout_dst2(c) * dtime
end if
qin_dst2(c) = qout_dst2(c)

! DUST 3:
! 1. flux with meltwater:
qout_dst3(c) = qflx_snow_percolation(c,j)*scvng_fct_mlt_dst3*(mss_dst3(c,j)/mss_liqice)
if (qout_dst3(c) > mss_dst3(c,j)) then
qout_dst3(c) = mss_dst3(c,j)
endif
mss_dst3(c,j) = mss_dst3(c,j) - qout_dst3(c) * dtime
if (qout_dst3(c)*dtime > mss_dst3(c,j)) then
qout_dst3(c) = mss_dst3(c,j)/dtime
mss_dst3(c,j) = 0._r8
else
mss_dst3(c,j) = mss_dst3(c,j) - qout_dst3(c) * dtime
end if
qin_dst3(c) = qout_dst3(c)

! DUST 4:
! 1. flux with meltwater:
qout_dst4(c) = qflx_snow_percolation(c,j)*scvng_fct_mlt_dst4*(mss_dst4(c,j)/mss_liqice)
if (qout_dst4(c) > mss_dst4(c,j)) then
qout_dst4(c) = mss_dst4(c,j)
endif
mss_dst4(c,j) = mss_dst4(c,j) - qout_dst4(c) * dtime
if (qout_dst4(c)*dtime > mss_dst4(c,j)) then
qout_dst4(c) = mss_dst4(c,j)/dtime
mss_dst4(c,j) = 0._r8
else
mss_dst4(c,j) = mss_dst4(c,j) - qout_dst4(c) * dtime
end if
qin_dst4(c) = qout_dst4(c)

end if
Expand Down

0 comments on commit 0dea14e

Please sign in to comment.