Context

Attempts to derive a SETTE test for SAS configurations by comparing restarts in restart ability tests has revealed a minor issue with sbc_tsc_b values differing between continuous runs and restarted runs covering the same period. Differences are found in halo regions which have no impact on results but invalidate the proposed new test.

Analysis

The cause has been traced to the initialisation carried out in trasbc.F90. Currently the code is (in shorter, pre-OMP form):

      IF( kt == nit000 ) THEN             !* 1st time-step
         IF( ln_rstart .AND.    &               ! Restart: read in restart file
              & iom_varid( numror, 'sbc_hc_b', ldstop = .FALSE. ) > 0 ) THEN
            IF(lwp) WRITE(numout,*) '          nit000-1 sbc tracer content field read in the restart file'
            zfact = 0.5_wp
            CALL iom_get( numror, jpdom_autoglo, 'sbc_hc_b', sbc_tsc_b(:,:,jp_tem) )   ! before heat content sbc trend
            CALL iom_get( numror, jpdom_autoglo, 'sbc_sc_b', sbc_tsc_b(:,:,jp_sal) )   ! before salt content sbc trend
         ELSE                                   ! No restart or restart not found: Euler forward time stepping
            zfact = 1._wp
            sbc_tsc(:,:,:) = 0._wp
            sbc_tsc_b(:,:,:) = 0._wp
         ENDIF
      ELSE                                !* other time-steps: swap of forcing fields
         zfact = 0.5_wp
         sbc_tsc_b(:,:,:) = sbc_tsc(:,:,:)
      ENDIF
      !                             !==  Now sbc tracer content fields  ==!
      DO jj = 2, jpj
         DO ji = fs_2, fs_jpim1   ! vector opt.
            sbc_tsc(ji,jj,jp_tem) = r1_rau0_rcp * qns(ji,jj)   ! non solar heat flux
            sbc_tsc(ji,jj,jp_sal) = r1_rau0     * sfx(ji,jj)   ! salt flux due to freezing/melting
         END DO
      END DO

Note that on a restart, sbc_tsc is not initialised and only inner points are subsequently set before the:

         sbc_tsc_b(:,:,:) = sbc_tsc(:,:,:)

instruction on the next time step. Thus the unused haloes can (and do) contain arbitrary values unless unset values are automatically set to zero by the compiler.

Recommendation

The simplest solution is to introduce a:

            sbc_tsc(:,:,:) = 0._wp

immediately before the first call to iom_get.

Of course for OMP operability what actually needs to be inserted is:

            DO jn = 1, jpts
!$OMP PARALLEL DO schedule(static) private(jj, ji)
               DO jj = 1, jpj
                  DO ji = 1, jpi
                     sbc_tsc(ji,jj,jn) = 0._wp  ! needed just to ensure haloes are consistent across restarts
                  END DO
               END DO
            END DO

Context

With the move to full dynamic allocation there is no need to retain a uniform jpi and jpj across all processing regions. Allowing these to vary will permit simplification by removing the possibility of 'ghost' rows or columns and therefore remove the need for the confusing number of start and end markers such as nldi and nlci. This work will also simplify the introduction of wider haloes and will be used as an opportunity to reduce the size of the lib_mpp.F90 code base by introducing generic functions in .h90 files that can be included multiple times with different preprocessor settings. This work builds on the developments completed in branches/2017/dev_r7832_HPC08_lbclnk_3rd_dim and #1880

Implementation

See https://forge.ipsl.jussieu.fr/nemo/wiki/2017WP/ROBUST-08_AndrewC-MPP_no_ghost for details and associated development branch: branches/2017/dev_r8126_ROBUST08_no_ghost

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.