trabbl not conserving tracers at the north-fold

[acc]

Context

v3.6_STABLE (and probably trunk, though not tested for trunk since trend diagnostics not yet implemented there).
Analyses of trend diagnostics and evidence from passive tracers has highlighted a conservation issue with the diffusive bbl code (nn_bbl_ldf=1). This first came to light with eORCA1 studies but is reproducible with ORCA2_LIM3

Analysis

The cause can be traced to north-fold points with zero bathymetry gradients across the fold. Consistent calculation of the diffusive coefficients relies on density gradients and bathymetry gradients taking opposite signs either side of the pivot, I.e. in routine bbl:

          zgdrho = (  za * ( zts(ji,jj+1,jp_tem) - zts(ji,jj,jp_tem))    &
             &      - zb * ( zts(ji,jj+1,jp_sal) - zts(ji,jj,jp_sal))  ) * vmask(ji,jj,1)
          !
          zsign = SIGN(  0.5, -zgdrho * REAL( mgrhv(ji,jj) )  )     !  sign of ( j-gradient * j-slope )
          ahv_bbl(ji,jj) = ( 0.5 - zsign ) * ahv_bbl_0(ji,jj)

but mgrhv (for example) is set in tra_bbl_init by:

                                       !* sign of grad(H) at u- and v-points
      mgrhu(jpi,:) = 0   ;   mgrhu(:,jpj) = 0   ;   mgrhv(jpi,:) = 0   ; mgrhv(:,jpj) = 0
      DO jj = 1, jpjm1
         DO ji = 1, jpim1
            mgrhu(ji,jj) = INT(  SIGN( 1.e0, gdept_0(ji+1,jj,mbkt(ji+1,jj)) - gdept_0(ji,jj,mbkt(ji,jj)) )  )
            mgrhv(ji,jj) = INT(  SIGN( 1.e0, gdept_0(ji,jj+1,mbkt(ji,jj+1)) - gdept_0(ji,jj,mbkt(ji,jj)) )  )
         END DO
      END DO

and when the gdept_0 values are equal the result is always positive. zgdrho, however, will always have opposite signs either side of the pivot. Thus ahv_bbl is non-zero one side of the pivot and zero on the other; leading to the non-conservation.

For example consider this T-pivot case:

    +---+                                  +---+
   | | C |  <------------\                  | A |
   v +---+                \                 +---+
     | B |-----------------o----------------| B |
     +---+                  \               +---+ ^
     | A |                   \------------> | C | |
     +---+                                  +---+

The cell marked B on the lhs is an inner point and its ahv_bbl value must match that calculated for the rhs inner point marked C. Any density gradient will be reversed by the folding of tracer values as indicated but the sign of grad(H) is unaltered in the cases of grad(H)=0.

Fix

The bbl should not be active when grad(H)=0 so the easiest solution is to force a zero coefficient in these cases. This can be done by setting the sign of the gradient to zero when the gradient is zero. This in turn sets zsign to 0.5 providing the sign function adheres to the f90 standard and sets sign of zero to be positive. Use key_nosignedzero if in doubt.

This code replacement in tra_bbl_init (trabbl.F90) achieves this:

                                       !* sign of grad(H) at u- and v-points; zero if grad(H) = 0
      mgrhu(:,:) = 0   ;   mgrhv(:,:) = 0
      DO jj = 1, jpjm1
         DO ji = 1, jpim1
            IF( gdept_0(ji+1,jj,mbkt(ji+1,jj)) - gdept_0(ji,jj,mbkt(ji,jj)) /= 0._wp ) THEN
               mgrhu(ji,jj) = INT(  SIGN( 1.e0, gdept_0(ji+1,jj,mbkt(ji+1,jj)) - gdept_0(ji,jj,mbkt(ji,jj)) )  )
            ENDIF
            !
            IF( gdept_0(ji,jj+1,mbkt(ji,jj+1)) - gdept_0(ji,jj,mbkt(ji,jj)) /= 0._wp ) THEN
               mgrhv(ji,jj) = INT(  SIGN( 1.e0, gdept_0(ji,jj+1,mbkt(ji,jj+1)) - gdept_0(ji,jj,mbkt(ji,jj)) )  )
            ENDIF
         END DO
      END DO

It is inelegant but only performed once at start-up. This has been verified to fix the issue in ORCA2_LIM3 with:

• diffusive bbl without any advective bbl i.e. nn_bbl_ldf=1 nn_bbl_adv=0
• diffusive bbl + advective bbl option 1 i.e. nn_bbl_ldf=1 nn_bbl_adv=1
• diffusive bbl + advective bbl option 2 i.e. nn_bbl_ldf=1 nn_bbl_adv=2

It also has been verified to fix the issue in the Met O. GO6 configuration at 1 deg res (ORCA1-CICE) being used as the UK submision to CMIP6, for

• diffusive bbl without any advective bbl i.e. nn_bbl_ldf=1 nn_bbl_adv=0.

Other options have not yet been checked.

[acc]

Bug fix applied at revision 8507:

trabbl.F90

@@ -547 +547 @@
       zmbk(:,:) = REAL( mbku_d(:,:), wp )   ;   CALL lbc_lnk(zmbk,'U',1.)   ;   mbku_d(:,:) = MAX( INT( zmbk(:,:) ), 1 )
       zmbk(:,:) = REAL( mbkv_d(:,:), wp )   ;   CALL lbc_lnk(zmbk,'V',1.)   ;   mbkv_d(:,:) = MAX( INT( zmbk(:,:) ), 1 )
 
-                                        !* sign of grad(H) at u- and v-points
-      mgrhu(jpi,:) = 0   ;   mgrhu(:,jpj) = 0   ;   mgrhv(jpi,:) = 0   ;   mgrhv(:,jpj) = 0
+                                        !* sign of grad(H) at u- and v-points; zero if grad(H) = 0
+      mgrhu(:,:) = 0   ;   mgrhv(:,:) = 0
       DO jj = 1, jpjm1
          DO ji = 1, jpim1
-            mgrhu(ji,jj) = INT(  SIGN( 1.e0, gdept_0(ji+1,jj,mbkt(ji+1,jj)) - gdept_0(ji,jj,mbkt(ji,jj)) )  )
-            mgrhv(ji,jj) = INT(  SIGN( 1.e0, gdept_0(ji,jj+1,mbkt(ji,jj+1)) - gdept_0(ji,jj,mbkt(ji,jj)) )  )
+            IF( gdept_0(ji+1,jj,mbkt(ji+1,jj)) - gdept_0(ji,jj,mbkt(ji,jj)) /= 0._wp ) THEN
+               mgrhu(ji,jj) = INT(  SIGN( 1.e0, gdept_0(ji+1,jj,mbkt(ji+1,jj)) - gdept_0(ji,jj,mbkt(ji,jj)) )  )
+            ENDIF
+            !
+            IF( gdept_0(ji,jj+1,mbkt(ji,jj+1)) - gdept_0(ji,jj,mbkt(ji,jj)) /= 0._wp ) THEN
+               mgrhv(ji,jj) = INT(  SIGN( 1.e0, gdept_0(ji,jj+1,mbkt(ji,jj+1)) - gdept_0(ji,jj,mbkt(ji,jj)) )  )
+            ENDIF
          END DO
       END DO
 

[acc]

And committed to trunk at revision r8509