Quadratic solution error

Brief summary of bug

At year 1577 of a spinup with Princeton atm forcing data, I got a quadratic solution error (b^2 - 4ac is negative) that stopped the model

General bug information

CTSM version you are using: release-clm5.0.24

Does this bug cause significantly incorrect results in the model's science? Doubtful

Configurations affected: 1850_DATM%GSWP3v1_CLM50%BGC-CROP_SICE_SOCN_MOSART_CISM2%NOEVOLVE_SWAV --res f09_g17_gl4

Details of bug

Here is output from cesm log file:

1275: quadratic ERROR: Quadratic solution error: b^2 1275: -5.648060482703832E-003 1275: ERROR in quadraticMod.F90 at line 55 1275: 1275: 1275: 1275: 1275: 1275: 1275: ENDRUN: 1275: ERROR: quadratic ERROR: Quadratic solution error: 1275:Image PC Routine 1275:cesm.exe 00000000019B220D Unknown 1275:cesm.exe 0000000001131C32 shr_abort_modmp 1275:cesm.exe 00000000004FA47F abortutils_mp_end 1275:cesm.exe 0000000000B7F8ED quadraticmod_mp_q 1275:cesm.exe 00000000008AFD21 photosynthesismod 1275:cesm.exe 00000000008AAC14 photosynthesismod 1275:cesm.exe 00000000006FEF8D canopyfluxesmod_m 1275:cesm.exe 00000000004FFD22 clm_driver_mp_clm 1275:cesm.exe 00000000004EF2CB lnd_comp_mct_mp_l 1275:cesm.exe 0000000000424E84 component_modmp 1275:cesm.exe 000000000040A5AD cime_comp_modmp 1275:cesm.exe 0000000000424B2C MAIN__ 1275:cesm.exe 000000000040829E Unknown 1275:libc-2.19.so 00002AAAB08F9B25 __libc_start_main 1275:cesm.exe 00000000004081A9 Unknown 1275:MPT ERROR: Rank 1275(g:1275) is aborting with 1275: Process ID: 47548, Host: r4i7n17, Program: 1275: MPT Version: HPE MPT 2.19 12/07/18 05:31:15 - 4ac is negative = b^2 - 4ac is negative Line Source Unknown Unknown 114 shr_abort_mod.F90 50 abortutils.F90 56 quadraticMod.F90 4158 PhotosynthesisMod.F90 3661 PhotosynthesisMod.F90 852 CanopyFluxesMod.F90 543 clm_driver.F90 456 lnd_comp_mct.F90 728 component_mod.F90 2712 cime_comp_mod.F90 125 cime_driver.F90 Unknown Unknown Unknown Unknown Unknown Unknown error code 1001. /gpfs/fs1/scratch/oleson/clm50_release-clm5.0.24_1deg_PrincetonV2_1850pAD/bld/cesm.exe

Important details of your setup / configuration so we can reproduce the bug

Here is the case directory:

/glade/work/oleson/release-clm5.0.24/cime/scripts/clm50_release-clm5.0.24_1deg_PrincetonV2_1850pAD

I'll look into this once cheyenne is back up.

The situation here is that net shaded canopy photosynthesis is strongly negative (-9400) while net sunlit canopy photosynthesis is positive (10). We only exit the calculations in ci_func_PHS prior to the quadratic solution if both shaded and sunlit net photosynthesis are negative. The negative net shaded canopy photosynthesis is caused by very large shaded leaf respiration (lmr_z_sha), which in turn is caused by the fact that the leaf to canopy scaling factor is large (vcmaxcintsha = 6575). The exposed leaf area is very small (1e-13; phytoplankton-size leaves). That combined with a very small fsha (near zero) creates the large leaf to canopy scaling factor. One fix is that we could have if statements in ci_func_PHS checking to see if an_sha and an_sun are >= 0 before proceeding with the quadratic solution for gs_mol_sha and gs_mol_sun, which is only valid if an_sha and an_sun are >= 0 anyway. There are calculations for gs_mol_sha and gs_mol_sun above this if an_sha < 0 and an_sun < 0. This is similar to what is done at the end of the routine when calculating fvalsun and fvalsha. A simulation would be needed of course to see what answer changes were. Diffs shown below.

diff --git a/glade/work/oleson/release-clm5.0.24/src/biogeophys/PhotosynthesisMod.F90 b/PhotosynthesisMod.F90
index b5f2a6781..c4e8dd011 100644
--- a/glade/work/oleson/release-clm5.0.24/src/biogeophys/PhotosynthesisMod.F90
+++ b/PhotosynthesisMod.F90
@@ -4128,10 +4128,13 @@ contains
     ! With an <= 0, then gs_mol = bbb

     ! Sunlit
-    cs_sun = cair - 1.4_r8/gb_mol * an_sun(p,iv) * forc_pbot(c)
-    cs_sun = max(cs_sun,10.e-06_r8)
+    if (an_sun(p,iv) >= 0._r8) then
+      cs_sun = cair - 1.4_r8/gb_mol * an_sun(p,iv) * forc_pbot(c)
+      cs_sun = max(cs_sun,10.e-06_r8)
+    end if

     if ( stomatalcond_mtd == stomatalcond_mtd_medlyn2011 )then
+      if (an_sun(p,iv) >= 0._r8) then
        term = 1.6_r8 * an_sun(p,iv) / (cs_sun / forc_pbot(c) * 1.e06_r8)
        aquad = 1.0_r8
        bquad = -(2.0 * (medlynintercept(patch%itype(p))*1.e-06_r8 + term) + (medlynslope(patch%itype(p)) * term)**2 / &
@@ -4142,7 +4145,9 @@ contains

        call quadratic (aquad, bquad, cquad, r1, r2)
        gs_mol_sun = max(r1,r2) * 1.e06_r8
+      end if

+      if (an_sha(p,iv) >= 0._r8) then
        ! Shaded
        cs_sha = cair - 1.4_r8/gb_mol * an_sha(p,iv) * forc_pbot(c)
        cs_sha = max(cs_sha,10.e-06_r8)
@@ -4157,13 +4162,16 @@ contains

        call quadratic (aquad, bquad, cquad, r1, r2)
        gs_mol_sha = max(r1,r2)* 1.e06_r8
+      end if
     else if ( stomatalcond_mtd == stomatalcond_mtd_bb1987 )then
+      if (an_sun(p,iv) >= 0._r8) then
        aquad = cs_sun
        bquad = cs_sun*(gb_mol - max(bsun*bbb(p),1._r8)) - mbb(p)*an_sun(p,iv)*forc_pbot(c)
        cquad = -gb_mol*(cs_sun*max(bsun*bbb(p),1._r8) + mbb(p)*an_sun(p,iv)*forc_pbot(c)*rh_can)
        call quadratic (aquad, bquad, cquad, r1, r2)
        gs_mol_sun = max(r1,r2)
-
+      end if
+      if (an_sha(p,iv) >= 0._r8) then
        ! Shaded
        cs_sha = cair - 1.4_r8/gb_mol * an_sha(p,iv) * forc_pbot(c)
        cs_sha = max(cs_sha,10.e-06_r8)
@@ -4173,6 +4181,7 @@ contains
        cquad = -gb_mol*(cs_sha*max(bsha*bbb(p),1._r8) + mbb(p)*an_sha(p,iv)*forc_pbot(c)*rh_can)
        call quadratic (aquad, bquad, cquad, r1, r2)
        gs_mol_sha = max(r1,r2)
+      end if
     end if

     ! Derive new estimate for cisun,cisha

ESCOMP / CTSM