Negative SIF value from the CLM SourceMods

XueliHuo commented 5 months ago

:bug: Your bug may already be reported! Please search on the issue tracker before creating a new issue.

Describe the bug

This is not a bug in DART. This is a bug in the SourceMods of CLM SIF module. Just would like to mention here in case the SIF SourceMods will be included in the future DART versions.

I ran CLM with the SIF SourceMods. It was an 80-ensemble free run without assimilating any observations. The simulation period is from 2003 to 2010. SIF should always be positive. However, a negative SIF value occurred in April 2005.

Error Message

CLM didn't blow up. Instead, it continued running until the end of the simulation period, though the negative SIF doesn't make sense.

Which model(s) are you working with?

CLM with the SIF module developed by Li et al 2022 (https://doi.org/10.1029/2021MS002747 https://github.com/ESCOMP/CTSM/issues/1672) in SourceModes. The version of CLM is release-cesm2.2.01

Screenshots

The negative SIF occurred in the PFT (with vegetation type as 10) within the grid cell indexed with (9,11) (latidx, lonidx) at one time step on the simulation date 2005-04-28 in ensemble member 37.

Version of DART

Which version of DART are you using? You can find the version using git describe --tags
v10.7.0-111-g70e6af803

Have you modified the DART code?

No
If your code changes are available on GitHub, please provide the repository.

Build information

Please describe:

The machine you are running on (e.g. windows laptop, NSF NCAR supercomputer Derecho). CHPC university of utah.
The compiler you are using (e.g. gnu, intel).
Intel

XueliHuo commented 5 months ago

Diving into the issue, we found that it was caused by a very large ftin (78642838.3848912) calculated in the TwoStream subroutine in the SurfaceAlbedoMod.F90 in the SourceMods for SIF. Here is a snapshot of the code:

           albd(p,ib) = h1/sigma + h2 + h3
            ftid(p,ib) = h4*s2/sigma + h5*s1 + h6/s1
            ! for narid viewing                       
     **ftin(p,ib) = vh4*s3/vsigma + vh5*s4 + vh6/s4**
            ! singularity                       
     if (abs(sigma) < 0.00000001_r8) then
        albd(p,ib) = -h1/avmu(p)/avmu(p)*(1._r8/2._r8/twostext(p)) + h2 + h3
               ftid(p,ib) = 1._r8/c1*(-1._r8/2._r8/twostext(p)*h1/avmu(p)/avmu(p)*(p3*(elai(p)+esai(p))*CI(patch%itype(p))+p4/2._r8/twostext(p))-d)*s2 + h5*s1 + h6/s1
        ftin(p,ib) = 1._r8/vc1*(-1._r8/2._r8/vtwostext(p)*vh1/avmu(p)/avmu(p)*(vp3*(elai(p)+esai(p))*CI(patch%itype(p))+vp4/2._r8/vtwostext(p))-d)*s3 + vh5*s4 + vh6/s4
     end if

It provides the solutions for albd and ftid when sigma is very small and doesn't provide the solution for ftin when vsigma is very small. Confirming with the developer Li, we modified the if statements in the singularity section as follows to solve this issue:

   ! singularity                       
     if (abs(sigma) < 0.00000001_r8) then
        albd(p,ib) = -h1/avmu(p)/avmu(p)*(1._r8/2._r8/twostext(p)) + h2 + h3
               ftid(p,ib) = 1._r8/c1*(-1._r8/2._r8/twostext(p)*h1/avmu(p)/avmu(p)*(p3*(elai(p)+esai(p))*CI(patch%itype(p))+p4/2._r8/twostext(p))-d)*s2 + h5*s1 + h6/s1
     end if

      if (abs(sigma) < 0.00000001_r8) then
             ftin(p,ib) = 1._r8/vc1*(-1._r8/2._r8/vtwostext(p)*vh1/avmu(p)/avmu(p)*(vp3*(elai(p)+esai(p))*CI(patch%itype(p))+vp4/2._r8/vtwostext(p))-d)*s3 + vh5*s4 + vh6/s4
     end if

This modification was tested and it prevented the negative SIF value from occurring on the date. I am currently running the 80-ensemble free run with this bug fix for the entire simulation period to make sure that the negative SIF value doesn't occur. I will confirm this once the free run is done.

XueliHuo commented 5 months ago

The free run with this bug fix for the entire simulation period from 1998-01-01 to 2011-01-01 is done. There are no negative SIF values anymore, and all SIF values are positive, proving the bug fix works.

hkershaw-brown commented 5 months ago

Thanks for recording this in an issue and all the work diagnosing this @XueliHuo. Super helpful for anyone working with SIFs.

braczka commented 5 months ago

Quick note to myself -- that the current cesm2_2 SourceMods included with DART, do not have these two fixes implemented yet for this bug (ftin, albedo fix), and also separate bug (degree of light saturation, #664 ). I will generate PR to update these changes in CLM.

NCAR / DART