variable partially uninitialized in traadv_tvd_zts

[mathiot]

Context

On the MetO machine, traadv_tvd_zts is very unstable. NaN appeared sometimes during the restart time step with the sub time stepping activated.

Analysis

In traadv_tvd_zts, ztrs(:,:,:,1) is inialized to ptb(:,:,:,jn). However, ztrs(:,:,1,2:3) are not initialised or computed and are used to compute the interior value of zwz (l. 496).

            DO jk = 2, jpkm1          ! Interior value
               DO jj = 2, jpjm1
                  DO ji = fs_2, fs_jpim1
                     zwz(ji,jj,jk) = 0.5_wp * pwn(ji,jj,jk) * ( ztrs(ji,jj,jk,jtn) + ztrs(ji,jj,jk-1,jtn) )
                     IF( jtaken == 0 ) zwzts(ji,jj,jk) = zwzts(ji,jj,jk) + zwz(ji,jj,jk)*zts(jk)           ! Accumulate time-weighted vertcal flux
                  END DO
               END DO
            END DO

Fix

As fix, I suggest, in traadv_tvd_zts, set ztrs(:,:,1,:) = 0.0. Maybe initialisation of the surface value to ptb(:,:,1,jn) for jtb, jtn and jta is better? other solution ? Could someone comment on it?

[acc]

A couple of comments:

1. This routine was modified by bugfix #1603 which I think benign but unnecessary since ztrs is used entirely within a loop over tracers.

2. Setting ztrs(:,:,:,1) should be sufficient since the leap-frog pointers are shuffled such that only set entries are used. For jl =1, only entry '1' appears on the RHS and entry '2' is set. For jl=2, entries '1' and '2' appear on the RHS and '3' is set. After that all should be defined.

If you are encountering undefined values in ztrs then they must be upstream in ptb itself unless I'm missing something more subtle?

[mathiot]

1. : I agree with you, only ztrs(:,:,:,jl) is used with the index jl=jtb/jtn/jta used between 1 and 3 whatever is the tracer.

2. : I agree ”for jl=1, only entry '1' appears on the RHS and entry '2' is set” but only for jk=2,jpkm1 (l.506, loop "interior value",“total advective trend”). For jl=2, ztrs(ji,jj,1,2) is used to compute zwz(ji,jj,2) and not set in the case jl=1 or in the intialisation.

I suggest to change the initialisation:

         ! antidiffusive flux on k
         zwz(:,:,1) = 0._wp        ! Surface value
         zwz_sav(:,:,:) = zwz(:,:,:)
         !
         ztrs(:,:,:,1) = ptb(:,:,:,jn)
         zwzts(:,:,:) = 0._wp

by:

         ! antidiffusive flux on k
         zwz(:,:,1) = 0._wp        ! Surface value
         zwz_sav(:,:,:) = zwz(:,:,:)
         !
         ztrs(:,:,:,1) = ptb(:,:,:,jn)
         ztrs(:,:,1,2) = ptb(:,:,1,jn)
         ztrs(:,:,1,3) = ptb(:,:,1,jn)
         zwzts(:,:,:) = 0._wp

[acc]

Yes, I see your point about the surface values not being initialised for the second and third entries of ztrs(:,:,:,n). The proposed fix is valid and necessary.

[mathiot]

Fix in the trunk at commit r6771 and in the 3.6 at commit r6770.