Neutral flat simulation requires bldT > 0.

A neutral simulation that has the floor covered with one large facet (see 001 in branch bss/example-simulations) crashes when the parameter bldT in namoptions section ENERGYBALANCE is not set. (default value is 0.). The simulation runs for values > 0.

Error stack: forrtl: error (73): floating divide by zero Image PC Routine Line Source u-dales 000000000078E7CF Unknown Unknown Unknown libpthread-2.17.s 00002B2CCE195370 Unknown Unknown Unknown u-dales 00000000007055BE modthermodynamics 419 modthermodynamics.f90 u-dales 0000000000701E61 modthermodynamics 313 modthermodynamics.f90 u-dales 00000000006F6DDA modthermodynamics 69 modthermodynamics.f90 u-dales 000000000070A6A2 MAIN__ 145 program.f90

Division by zero seems to come from this line: presh(k) = presh(k-1)rdocp - grav(pref0*rdocp)dzf(k-1) / (cpthvf(k-1))

@ivosuter any ideas how to remove bldT requirement for neutral simulations?

We are getting a division by zero because the slab averaged virtual potential temp at the level of the flat single block: thvf(kb) = 0.

I think this happens because thermodynamics is called in modstartup for cold start simulations (subroutine readinitfiles) and this is before modibm will overwrites the zero values within these kind of blocks (it will set it equal to kb+1 to avoid subgrid diffusion of temp). I assume these calls of thermodynamics should have a switch around them regardless:ltempeq?

In modmpi there is a switch for when a whole layer is filled with buildings - this avoids multiplying by the slab averaged masking matrix to avoid getting zero values just like this. However, if all of the values that are summed are zero regardless of this then we still get a zero.

Options I can think of:

• Change default bldT to something other than 0.
• Edit the code in slabave_xy to also check if it outputs a zero when summing and use the value from kb+1 if this is the case.

To note:

• We should run this example case without any blocks at all to test subroutine bottom anyway.
• Should we put the calls to thermodynamics in modstartup behind a ltempeq switch?

If the ltempeq switch is appropriate and it fixes the bug, then I would say we implement that. Definitely not changing the default of bldT, no parameter from the ENERGYBALANCE section should be used when lEB = .false. (that's why we moved nfcts to WALLS).

I'm not sure about changing slabave_xy to use kb+1 if all kb are zero, sounds like this could too easily lead to unintended results...

Note that bldT is also used as the initial building temperature in non-energybalance simulations, and therefore adapting the default value should be considered either way, see discussion in https://github.com/uDALES/u-dales/issues/39#issuecomment-624111119