SMS_Lm3_D_Mmpi-serial.1x1_brazil.I2000Clm50FatesCruRsGs.izumi_intel.clm-FatesColdHydro fails when t_soisno is initialized to 272K instead of 274K

[olyson]

Brief summary of bug

When t_soisno is initialized to 272K instead of 274K per PR #2355 , the test SMS_Lm3_D_Mmpi-serial.1x1_brazil.I2000Clm50FatesCruRsGs.izumi_intel.clm-FatesColdHydro fails, as described in PR #2371.

General bug information

The test fails at line 771 in FatesHydroWTFMod.F90:
psi = this%psi_sat*(th/this%th_sat)**(-this%beta)

where apparently th is zero.

th is passed in as h2osoi_liqvol_shell at line 578 in FatesPlantHydraulicsMod.F90:

cohort_hydr%psi_aroot(j) = csite_hydr%wrf_soil(j)%p%psi_from_th(csite_hydr%h2osoi_liqvol_shell(j,1))

where h2osoi_liqvol_shell is zero presumably because the soil is frozen.

CTSM version you are using: ctsm5.1.dev167-1-ge4dde90ef

Does this bug cause significantly incorrect results in the model's science? Model doesn't run

Configurations affected: Presumably any FATES hydro cold start would fail

Details of bug

This can be fixed by using a use_fates in TemperatureType.F90 such that the soil temperature is initialized to 272K if use_fates=.false., and 274K if use_fates=.true., but ideally it seems like fates should be able to cold start from an initial temperature of 272K. Particularly since if excess ice is on, the soil temperature is currently initialized to 5K below freezing, so presumably that configuration would fail as well.

Important output or errors that show the problem

Traceback is:

forrtl: error (73): floating divide by zero
Image              PC                Routine            Line        Source       
cesm.exe           00000000047FC6DB  Unknown               Unknown  Unknown
libpthread-2.17.s  00002B74A7C0D630  Unknown               Unknown  Unknown
cesm.exe           00000000047CE001  Unknown               Unknown  Unknown
cesm.exe           0000000003D270A5  fateshydrowtfmod_         771  FatesHydroWTFMod.F90
cesm.exe           0000000001DE07E5  fatesplanthydraul         578  FatesPlantHydraulicsMod.F90
cesm.exe           000000000178CF34  edcohortdynamicsm         264  EDCohortDynamicsMod.F90

[wwieder]

Thanks for looking into this @olyson. Interested to hear what @adrifoster, @rgknox, or @glemieux think here, but as a short term fix to let us move ahead with testing I think your simple fix for different initialized temperatures with fates vs. non-fates make sense.

[rgknox]

Thanks for pointing this out @olyson . My sense is that this is pointing out a vulnerability in the procedure: psi_from_th(). Initializing it with non frozen soils should work for that step, but it should just give us more problems later in the simulation.
It looks like we have provisions for saturated soils, but not completely frozen soils:

https://github.com/NGEET/fates/blob/main/biogeophys/FatesHydroWTFMod.F90#L768

I might take a look at how CLM calculates potential in completely frozen soils and use that as a guide.

[rgknox]

It is odd though, if its a div0, then the saturation water content is the zero, which should be a soil parameter.

[rgknox]

Noting that this is where we setup the saturated water content in a cold-start: https://github.com/NGEET/fates/blob/sci.1.72.0_api.33.0.0/biogeophys/FatesPlantHydraulicsMod.F90#L1575-L1581

[olyson]

I printed out the values involved in the equation. th was zero while th_sat was non-zero. However, that is then 0 raised to the negative beta power and I assume that the compiler is solving that as 1/0^(beta) and that's where the divide by zero comes in?

[rgknox]

That's really interesting, I did not expect that, but then again, it makes complete sense.

[olyson]

See NGEET/fates#1163 for corresponding FATES issue.