trasbc with key_vvl

[jamesharle]

Around ln 160 in trasbc.F90 the surface boundary condition content trend (sbc_tsc) is populated having been previously allocated, but not intialised (from a cold start). It's populated for i=2:end-1 and j=2:end, which struck me as a little odd (why end-1 in i and end in j). I first noticed this as the uninitialised row and columns occasionally (depending on compile options etc) contain NaNs? and blow up my model run. I first thought that these uninitialised cells shouldn't matter as they're in the halo, which when needed are updated with an exchange with neighbouring processors, thus removing the uninitialised values. However if I've tracked this correctly the sbc_tsc array is used to update tsb (aka ptb) in tranxt.F90 (around ln 322 and 328 - subroutine tra_nxt_vvl). This is looped over the full i and j so my NaN jumps into tsb (ptb). Now when tra_adv and tra_adv_tvd are called the NaN from tsb (ptb) makes it's way into the halo of the upstream flux zwx & zwy (~ln 128 traadv_tvd.F90) . In calculating the total advective trend (~ln 157) the NaN jumps to an interior ocean point as ztra at j=2 relies on zwy(j) - zwy(j-1) - this is subsequently propagated into tsa (aka pta) and tsn (aka ptn) so in 3 time steps the model blows up. Now I can solve this by intialising sbc_tsc(:,:) = 0 around ln 147 in the trasbc.F90 routine or by looping over i=1:end j=1:end at ln 160 in trasbc.F90 - both give slightly different answers. Looking at the code I would have thought one would loop over the full i and j for sbc_tsc to make the halo calculation in the total advective trend consistent with the
neighbouring processor - as by setting sbc_tsc to zero in the first row and first and last columns would remove the correction to the surface tracer field prior to the flux calculation and be inconsistent with the neighbouring processor interior points - e.g the last interior row of the processor to the south would have the surface flux correction and therefore there will be an albeit small difference in the horizontal flux calculation (there isn't, as far as I can see, lbc_lnk prior to this - hence the propagation of the NaN into the interior).

[charris]

I can confirm this is still a problem with the 3.4 beta release and agree that the simplest way to fix it would appear to be changing the loop indices as James suggests. To get complete restartability I also find I need to add lbclnk calls on tsb in tra_nxt (in tra_nxt_vvl or after the call to it). I haven't had chance to investigate why this is required and whether this is showing there might be some other problems so would appreciate someone else looking at the code before I just commit my changes.

Chris

[rblod]

I've checked on a idealized configuration with key_vvl. Same problem ocured. It seems to me the correct fix is to change the loop index as suggested by James.Initialization to 0 leed to incorrect results.
Thus a similar problem seems to appear in traqsr (ln_traqsr and key_vvl) line 248 where qsr_hc is computed only in the interior domain and leads to incorrect value of before temperature (extra term addded in tra_nxt_vvl). Doing the concerned loop on the whole domain in both trasbc and traqsr allow me to have full reproducibility (standard GYRE with key_vvl added). I guess doing the full loop in traqsr is equivalent do add lbc on tsb in tra_nxt_vvl as suggested by Chris.
The only reason I would see to compute these termes on the interior domain only would be if the forcing fields were not read/defined on the halo. I have to confess my total ignorance on this, Sebastien explained me thousands of time but I systematically forget.
Finally It would be appreciable at some point to have some feedback from people who introduced this astonishing feature of forcing at dt+0.5.

Rachid

[charris]

Thanks Rachid. I confirm that the loop change in traqsr has the same effect as my lbclnk calls in tranxt. I'll commit the changes in traqsr and trasbc early next week unless someone has done it first (or someone disagrees before then).

[gm]

The deep reason of this "astonish feature" is described in Leclair & Madec (Ocean Modelling 2010).
I check and agree there is a Pb ! Sorry about that !

The reason why these termes on the interior domain only is that historically all tendency term routines perform the calculation in the interior only: at the beginning of tranxt, a single lbc is applied on tsa and every thing works fine. Except that with the modified leap-frog/time filter, the forcing term are also used in the time filter (in order to remove their effect on the filter and so exactly conserve heat and salt content).

So the solution is to change the loop index in both trasbc and traqsr when computing (sbc_tsc sbc_tsc _b and qsr_hc, qsr_hc_b). In traqsr, only the 2 bands case is concerned, since it is already done in the RGB case.
I prefer this to the addition of a extra lbc on before t & s (tsb) in tranxt, something that also solve the issue.

Gurvan

[charris]

Loop index changes made in traqsr and trasbc (see [3275]). If anyone thinks equivalent changes are required in the non-vvl case then please add them, but for the moment this ticket will be closed.