Bug in diaar5.F90 (v3.6_STABLE & trunk)

[gm]

Context

key_diaar5 defined (diagnostics used for CMIP6)
When double diffusion parameterization is used (key_zdfddm), the diagnostic of the work done against stratification by vertical mixing ( iom_use( 'tnpeo' ) ) is wrong
In addition there is a out-of-bound : the jk-do loop should start at jk=2, not 1 ! In case of lk_zdfddm=F we can also start with jk=2 since act(:,:,1) is always zero.

This concerns the v3.6_STABLE and trunk versions.

Analysis

The portion of diaar5.F90 involved in the calculation is:

      IF( iom_use( 'tnpeo' )) THEN    
      ! Work done against stratification by vertical mixing
      ! Exclude points where rn2 is negative as convection kicks in here and
      ! work is not being done against stratification
          CALL wrk_alloc( jpi, jpj, zpe )
          zpe(:,:) = 0._wp
          IF( lk_zdfddm ) THEN
             DO ji=1,jpi
                DO jj=1,jpj
                   DO jk=1,jpk
                      zrw =   ( gdepw_n(ji,jj,jk  ) - gdept_n(ji,jj,jk) )   &
                         &  / ( gdept_n(ji,jj,jk-1) - gdept_n(ji,jj,jk) )
                      !
                      zaw = rab_n(ji,jj,jk,jp_tem) * (1. - zrw) + rab_n(ji,jj,jk-1,jp_tem)* zrw
                      zbw = rab_n(ji,jj,jk,jp_sal) * (1. - zrw) + rab_n(ji,jj,jk-1,jp_sal)* zrw
                      !
                      zpe(ji, jj) = zpe(ji, jj) - MIN(0._wp, rn2(ji,jj,jk)) * &
                           &       grav * (  avt(ji,jj,jk) * zaw * (tsn(ji,jj,jk-1,jp_tem) - tsn(ji,jj,jk,jp_tem) )  &
                           &               - fsavs(ji,jj,jk) * zbw * (tsn(ji,jj,jk-1,jp_sal) - tsn(ji,jj,jk,jp_sal) ) )

                   ENDDO
                ENDDO
             ENDDO
          ELSE
             DO ji = 1, jpi
                DO jj = 1, jpj
                   DO jk = 1, jpk
                       zpe(ji,jj) = zpe(ji,jj) + avt(ji, jj, jk) * MIN(0._wp,rn2(ji, jj, jk)) * rau0 * e3w_n(ji, jj, jk)
                   ENDDO
                ENDDO
             ENDDO
          ENDIF
          CALL lbc_lnk( zpe, 'T', 1._wp)         
          CALL iom_put( 'tnpeo', zpe )
          CALL wrk_dealloc( jpi, jpj, zpe )
      ENDIF

When lk_zdfddm = F, the calculation is OK and corresponds to a vertical sum of rau0 * avt * N²
When lk_zdfddm = T, the loop involved both a multiplication by rn2 (more precisely MIN( 0, rn2 ) ) AND by the vertical diffusive density flux taking into account avt and avs.
It is almost as if the vertical sum of rau0 * avt * N⁴ were computed (avt is not significantly different from avs in most of the world ocean...)

Fix

A solution is to replace the lk_zdfddm = True case by:

         IF( ln_zdfddm ) THEN
            DO jk = 2, jpk
               DO jj = 1, jpj
                  DO ji = 1, jpi
                     IF( rn2(ji,jj,jk) > 0._wp ) THEN
                        zrw =   ( gdepw_n(ji,jj,jk  ) - gdept_n(ji,jj,jk) )   &
                           &  / ( gdept_n(ji,jj,jk-1) - gdept_n(ji,jj,jk) )
                        !
                        zaw = rab_n(ji,jj,jk,jp_tem) * (1. - zrw) + rab_n(ji,jj,jk-1,jp_tem)* zrw
                        zbw = rab_n(ji,jj,jk,jp_sal) * (1. - zrw) + rab_n(ji,jj,jk-1,jp_sal)* zrw
                        !
                        zpe(ji, jj) = zpe(ji, jj)            &
                           &        -  grav * (    avt(ji,jj,jk) * zaw * (tsn(ji,jj,jk-1,jp_tem) - tsn(ji,jj,jk,jp_tem) )  &
                           &                   - fsavs(ji,jj,jk) * zbw * (tsn(ji,jj,jk-1,jp_sal) - tsn(ji,jj,jk,jp_sal) ) )
                     ENDIF
                  END DO
               END DO
             END DO
          ELSE
            DO jk = 1, jpk
               DO ji = 1, jpi
                  DO jj = 1, jpj
                     zpe(ji,jj) = zpe(ji,jj) + avt(ji, jj, jk) * MIN(0._wp,rn2(ji, jj, jk)) * rau0 * e3w_n(ji, jj, jk)
                  END DO
               END DO
            END DO
         ENDIF

NB: I also invert the-do-loop order (i.e. from i-j-k to k-j-i) as it is faster due to the model 3D array structure

[timgraham]

Fixed at revision r8077 of nemo_v3_6_STABLE and r8078 of the trunk

[frrh]

In 10870:

Error: Failed to load processor CommitTicketReference

No macro or processor named 'CommitTicketReference' found