ln_traldf_hor option crashes with zps coordinate

[hadcv]

Context

When using ln_traldf_hor = .true. with the zps coordinate, NEMO crashes with:

  ===>>> : E R R O R
 
          ===========
 
   stp_ctl: |ssh| > 20 m  or  |U| > 10 m/s  or  S <= 0  or  S >= 100  or  NaN en
 counter in the tests
 
 kt 3 |ssh| max   -Infinity at i j    43   3    found in 36 MPI tasks, spread ou
 t among ranks  0 to 35
 kt 3 |U|   max   -Infinity at i j k  43   3  1 found in 36 MPI tasks, spread ou
 t among ranks  0 to 35
 kt 3 Sal   min    Infinity at i j k  43   3  1 found in 36 MPI tasks, spread ou
 t among ranks  0 to 35
 kt 3 Sal   max   -Infinity at i j k  43   3  1 found in 36 MPI tasks, spread ou
 t among ranks  0 to 35

This affects both the trunk and r4.0-HEAD.

Analysis

​tra_ldf_lap requires gtsu/gtsv/gtui/gtvi to be defined (pgu/pgv/pgui/pgvi below):

147            IF( ln_zps ) THEN                             ! set gradient at bottom/top ocean level
148               DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 )                              ! bottom
149                  ztu(ji,jj,mbku(ji,jj)) = zaheeu(ji,jj,mbku(ji,jj)) * pgu(ji,jj,jn)
150                  ztv(ji,jj,mbkv(ji,jj)) = zaheev(ji,jj,mbkv(ji,jj)) * pgv(ji,jj,jn)
151               END_2D
152               IF( ln_isfcav ) THEN                             ! top in ocean cavities only
153                  DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 )
154                     IF( miku(ji,jj) > 1 )   ztu(ji,jj,miku(ji,jj)) = zaheeu(ji,jj,miku(ji,jj)) * pgui(ji,jj,jn)
155                     IF( mikv(ji,jj) > 1 )   ztv(ji,jj,mikv(ji,jj)) = zaheev(ji,jj,mikv(ji,jj)) * pgvi(ji,jj,jn)
156                  END_2D
157               ENDIF
158            ENDIF

However, in ​step.F90 these are only defined when l_ldfslp = .true.:

173         !  LATERAL  PHYSICS
174         !
175         IF( l_ldfslp ) THEN                             ! slope of lateral mixing
176                            CALL eos( ts(:,:,:,:,Nbb), rhd, gdept_0(:,:,:) )               ! before in situ density
177   
178            IF( ln_zps .AND. .NOT. ln_isfcav)                                    &
179               &            CALL zps_hde    ( kstp, Nnn, jpts, ts(:,:,:,:,Nbb), gtsu, gtsv,  &  ! Partial steps: before horizontal gradient
180               &                                          rhd, gru , grv    )       ! of t, s, rd at the last ocean level
181   
182            IF( ln_zps .AND.       ln_isfcav)                                                &
183               &            CALL zps_hde_isf( kstp, Nnn, jpts, ts(:,:,:,:,Nbb), gtsu, gtsv, gtui, gtvi,  &  ! Partial steps for top cell (ISF)
184               &                                          rhd, gru , grv , grui, grvi   )       ! of t, s, rd at the first ocean level
185            IF( ln_traldf_triad ) THEN
186                            CALL ldf_slp_triad( kstp, Nbb, Nnn )             ! before slope for triad operator
187            ELSE
188                            CALL ldf_slp     ( kstp, rhd, rn2b, Nbb, Nnn )   ! before slope for standard operator
189            ENDIF
190         ENDIF

which ​only occurs when:

250         IF(  nldf_tra == np_lap_i .OR. nldf_tra == np_lap_it .OR. &
251            & nldf_tra == np_blp_i .OR. nldf_tra == np_blp_it  )   l_ldfslp = .TRUE.    ! slope of neutral surfaces required

For ln_traldf_hor = .true. with the zps coordinate, this condition is not met (nldf_tra = np_lap/np_blp) and gtsu/gtsv/gtui/gtvi are therefore not defined.

Fix

Forcing the computation of gtsu/gtsv/gtui/gtvi resolves the issue:

 !  LATERAL  PHYSICS
      !
-     IF( l_ldfslp ) THEN                             ! slope of lateral mixing
-                        CALL eos( ts(:,:,:,:,Nbb), rhd, gdept_0(:,:,:) )               ! before in situ density
+        IF( ln_zps .OR. l_ldfslp ) CALL eos( ts(:,:,:,:,Nbb), rhd, gdept_0(:,:,:) )    ! before in situ density

         IF( ln_zps .AND. .NOT. ln_isfcav)                                    &
            &            CALL zps_hde    ( kstp, Nnn, jpts, ts(:,:,:,:,Nbb), gtsu, gtsv,  &  ! Partial steps: before horizontal gradient
            &                                          rhd, gru , grv    )       ! of t, s, rd at the last ocean level

         IF( ln_zps .AND.       ln_isfcav)                                                &
            &            CALL zps_hde_isf( kstp, Nnn, jpts, ts(:,:,:,:,Nbb), gtsu, gtsv, gtui, gtvi,  &  ! Partial steps for top cell (ISF)
            &                                          rhd, gru , grv , grui, grvi   )       ! of t, s, rd at the first ocean level
+     IF( l_ldfslp ) THEN                             ! slope of lateral mixing
         IF( ln_traldf_triad ) THEN
                         CALL ldf_slp_triad( kstp, Nbb, Nnn )             ! before slope for triad operator
         ELSE
                         CALL ldf_slp     ( kstp, rhd, rn2b, Nbb, Nnn )   ! before slope for standard operator
         ENDIF
      ENDIF

An alternative solution might be to include np_lap/np_blp in the conditions for which l_ldfslp = .true.. However tra_ldf_lap does not define akz/ah_wslp2, which are required by tra_zdf when l_ldfslp = .true..

[smasson]

I report Gurvan's comments:

• OK for the fix. This bug is in the code for decades...
• BTW, in realistic configuration, using an horizontal diffusion with a zps is a strange idea (which explain why we never faced this bug).

[hadcv]

In 15005:

Error: Failed to load processor CommitTicketReference

No macro or processor named 'CommitTicketReference' found