bugs on nSubsteps

[weiyuan-jiang]

It was reported and confirmed by me that nSubSteps becomes zero here ( both main and develop branch)

GOCART/Process_Library/GOCART2G_Process.F90

Lines 1579 to 1580 in c393e57

	nSubSteps = ceiling(cdt / dt_cfl)
	dt = cdt/nSubSteps

It turns out dt_cfl is a big negative value when it crashes. We can add a protection to set nSubSteps to zero. But I am wondering if the bug is from the other place that makes dt_cfl negative

@amdasilva

[RatkoVasic-NOAA]

@vbuchard @weiyuan-jiang
Even after fix #314 , there is one more case in subroutine SettlingSolver when model crashes with division by zero. It happens in boundary i,j point (1,1), and probably is not used in the model, so only debug mode can catch that.

What happens is that VSETTLE (or VS) is very small negative number (~ -10^-7)
TAU becomes very small negative number (~ -10^-7)
DT_CFL becomes very large negative number (~ -10^7)
since DT_CFL is not greater than CDT:   if (dt_cfl > cdt) or  if (-10^7 > 720) then
we get nSubSteps = ceiling(cdt / dt_cfl) --- gets close (rounded) to zero
and dt = cdt/nSubSteps --- getting division with zero

Can you please consider changing logic to include this extreme condition when division with zero appears.
Thanks!

[RatkoVasic-NOAA]

Here is suggestion, it worked for me:

          dp_  = delp(i,j,:)
          tau_ = tau(i,j,:)

          dt_cfl  = 1 / maxval(tau_)

          dt_cfl  = abs(dt_cfl)

Just this one line: dt_cfl = abs(dt_cfl)

[weiyuan-jiang]

If this line dt_cfl = abs(dt_cfl) is added, it probably would jump to (dt_cfl >cdt) and set nSubSteps=1. The original fix will skip the loop do iit=1, nSubSteps. I am not sure which one is desirable. @vbuchard

[weiyuan-jiang]

After tracing the problem, I think the real issue is in dz. According to the comments in the code, the variable hghte may be from (0,km) or (1:km+1). The lines like

GOCART/Process_Library/GOCART2G_Process.F90

Line 1173 in c393e57

hsurf => hghte(i1:i2,j1:j2,km)

are suspicious if the index range is (1, km+1). It would lead to dz(:,:, km) = 0 from the lines like

GOCART/Process_Library/GOCART2G_Process.F90

Line 1188 in c393e57

dz(:,:,k) = hghte(:,:,k+dk)-hsurf(:,:)

The subroutines here with hghte are called by the other subroutines which are out of control here. So we better take a closer look at this issue. I remembered some one complained dz became zero and the model crashed, but could not find it now. @RatkoVasic-NOAA @vbuchard

[weiyuan-jiang]

Also, Should this block always be called to guarantee vs > 0 when crashing?

GOCART/Process_Library/GOCART2G_Process.F90

Lines 1570 to 1574 in 71e75ca

	if(present(correctionMaring)) then
	if (correctionMaring) then
	vsettle = max(1.0e-9, vsettle - v_upwardMaring)
	endif
	endif

[weiyuan-jiang]

Going a little bit deeper, I found some calls without the correction ( and the problem happens when there is no correction) :

GOCART/ESMF/GOCART2G_GridComp/CA2G_GridComp/CA2G_GridCompMod.F90

Line 1035 in 837737f

call Chem_SettlingSimple (self%km, self%klid, self%diag_Mie, n, self%cdt, MAPL_GRAV, &

GOCART/ESMF/GOCART2G_GridComp/NI2G_GridComp/NI2G_GridCompMod.F90

Line 791 in 837737f

call Chem_SettlingSimple (self%km, self%klid, self%diag_Mie, 1, self%cdt, MAPL_GRAV, &

and call with correction

GOCART/ESMF/GOCART2G_GridComp/DU2G_GridComp/DU2G_GridCompMod.F90

Line 974 in 837737f

call Chem_SettlingSimple (self%km, self%klid, self%diag_Mie, n, self%cdt, MAPL_GRAV, &

I will leave this to someone who knows about the issue. It is beyond the scope of my knowledge.

[vbuchard]

Tagging @pcolarco and @amdasilva for visibility about this issue

[weiyuan-jiang]

In general, I agree with @RatkoVasic-NOAA 's change that would move one time step forward. Having said that, I also notice a few subroutines with name "Orig" skip the time forwarding. For example

GOCART/Process_Library/GOCART2G_Process.F90

Line 2091 in c393e57

nSubSteps = 0

[lipan-NOAA]

@weiyuan-jiang Would you like to submit a PR for this issue?

[weiyuan-jiang]

I think nSubsteps would not be 0 according to the logic . But since @RatkoVasic-NOAA has a better solution, I open the PR #323 here.

[mathomp4]

@RatkoVasic-NOAA @lipan-NOAA We have released a GOCART v2.4.1 which has these fixes from UFS:

• Update for ESMF8.8.0 and MAPL2.53.0 #310
• UFS Wet deposition and fengsha updates #305
• avoid dividing nSubsteps when it is zero #314

NOTE: It does not have #323 because we want to test that here before merging it in.