utrd_zdf and vtrd_zdf are velocities minus momentum trends

[chrenkl]

In the subroutine for vertical diffusion using a implicit time-stepping (dyn_zdf_imp in dynzdf_imp.F90), the momentum trends ua and va are integrated in time and become the after velocities (lines 122-136). Thus, in lines 90-95 of dynzdf.F90, when the momentum trend diagnostics are supposed to be computed, momentum trends (ztrdu, ztrdv) are subtracted from velocities (ua, va).

Therefore,

90   IF( l_trddyn )   THEN                      ! save the vertical diffusive trends for further diagnostics
91      ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
92      ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
93      CALL trd_dyn( ztrdu, ztrdv, jpdyn_zdf, kt )
94      CALL wrk_dealloc( jpi, jpj, jpk, ztrdu, ztrdv ) 
95   ENDIF

should be changed to

90   IF( l_trddyn )   THEN                      ! save the vertical diffusive trends for further diagnostics
91      ztrdu(:,:,:) = ( ua(:,:,:) - ub(:,:,:) ) / r2dt - ztrdu(:,:,:)
92      ztrdv(:,:,:) = ( va(:,:,:) - vb(:,:,:) ) / r2dt - ztrdu(:,:,:)
93      CALL trd_dyn( ztrdu, ztrdv, jpdyn_zdf, kt )
94      CALL wrk_dealloc( jpi, jpj, jpk, ztrdu, ztrdv ) 
95   ENDIF

as it is done in dynspg.F90. I have not checked whether this is also the case for vertical diffusion with explicit time-stepping.

[chrenkl]

It seems like my initial solution did not quite work. I have since modified the code (see attachment). This might not be the most straightforward solution, but it seems to work.

[chrenkl]

Modified dynzdf.F90

[nicolasmartin]

The status of this ticket is not clear.

2 commits (r6131 & r6132) have been linked to it dating back to the end of 2015 but the ticket had not been closed. Then the ZDF routines have been vastly re-factored in 2017 before merging with the trunk (see #1883).
Anyway the initial fix is already applied in the new routine:

IncludeSource(/branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf.F90, rev=8178, start=384, end=389)?

But a decision has to be made regarding the last fix suggested by Christoph (extract of dynzdf.F90 attached)

   SUBROUTINE dyn_zdf( kt )
      !!----------------------------------------------------------------------
      !!                  ***  ROUTINE dyn_zdf  ***
      !!
      !! ** Purpose :   compute the vertical ocean dynamics physics.
      !!---------------------------------------------------------------------
      !!
      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
      !
      REAL(wp), POINTER, DIMENSION(:,:,:) ::  ztrdu, ztrdv
! CR 20161028 --- begin ---
      INTEGER                             ::  jk   ! dummy loop indices
      REAL(wp)                            ::  z1_2dt 
      REAL(wp), POINTER, DIMENSION(:,:)   ::  zue, zve
      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zua, zva 
! CR 20161028 --- end ---
      !!---------------------------------------------------------------------
      !
      IF( nn_timing == 1 )  CALL timing_start('dyn_zdf')
! CR 20161028 --- begin ---
!       CALL wrk_alloc( jpi,jpj,jpk, zua, zva )
!       CALL wrk_alloc( jpi,jpj, zue, zve )
! CR 20161028 --- end ---
      !
      !                                          ! set time step
      IF( neuler == 0 .AND. kt == nit000     ) THEN   ;   r2dt =      rdt   ! = rdtra (restart with Euler time stepping)
      ELSEIF(               kt <= nit000 + 1 ) THEN   ;   r2dt = 2. * rdt   ! = 2 rdttra (leapfrog)
      ENDIF

      IF( l_trddyn )   THEN                      ! temporary save of ta and sa trends
         CALL wrk_alloc( jpi, jpj, jpk, ztrdu, ztrdv ) 
         ztrdu(:,:,:) = ua(:,:,:)
         ztrdv(:,:,:) = va(:,:,:)
      ENDIF

      SELECT CASE ( nzdf )                       ! compute lateral mixing trend and add it to the general trend
      !
      CASE ( 0 )   ;   CALL dyn_zdf_exp( kt, r2dt )      ! explicit scheme
      CASE ( 1 )   ;   CALL dyn_zdf_imp( kt, r2dt )      ! implicit scheme
      !
      CASE ( -1 )                                        ! esopa: test all possibility with control print
                       CALL dyn_zdf_exp( kt, r2dt )
                       CALL prt_ctl( tab3d_1=ua, clinfo1=' zdf0 - Ua: ', mask1=umask,               &
                          &          tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
                       CALL dyn_zdf_imp( kt, r2dt )
                       CALL prt_ctl( tab3d_1=ua, clinfo1=' zdf1 - Ua: ', mask1=umask,               &
                          &          tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
      END SELECT

      IF( l_trddyn )   THEN                      ! save the vertical diffusive trends for further diagnostics
! CR 20161028         ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
! CR 20161028         ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
! CR 20161028 --- begin ---
      ! In dyn_zdf_imp, ua & va become velocities and are not trends anymore. Thus, we have to back-calculate
      ! in order to obtain the trend due to vertical diffusion.
      ! Before, replace barotropic component: (copied from dynnxt.F90))
      ! Ensure below that barotropic velocities match time splitting estimate
      ! Compute actual transport and replace it with ts estimate at "after" time step
         CALL wrk_alloc( jpi,jpj,jpk, zua, zva )
         CALL wrk_alloc( jpi,jpj, zue, zve )
         zue(:,:) = fse3u_a(:,:,1) * ua(:,:,1) * umask(:,:,1)
         zve(:,:) = fse3v_a(:,:,1) * va(:,:,1) * vmask(:,:,1)
         DO jk = 2, jpkm1
            zue(:,:) = zue(:,:) + fse3u_a(:,:,jk) * ua(:,:,jk) * umask(:,:,jk)
            zve(:,:) = zve(:,:) + fse3v_a(:,:,jk) * va(:,:,jk) * vmask(:,:,jk)
         END DO
         DO jk = 1, jpkm1
            ua(:,:,jk) = ( ua(:,:,jk) - zue(:,:) * hur_a(:,:) + ua_b(:,:) ) * umask(:,:,jk)
            va(:,:,jk) = ( va(:,:,jk) - zve(:,:) * hvr_a(:,:) + va_b(:,:) ) * vmask(:,:,jk)
         END DO
         !
         z1_2dt = 1._wp / (2. * rdt)        ! Euler or leap-frog time step
         IF( neuler == 0 .AND. kt == nit000 )   z1_2dt = 1._wp / rdt     
         !
         zua(:,:,:) = ( ua(:,:,:) - ub(:,:,:) ) * z1_2dt - ztrdu
         zva(:,:,:) = ( va(:,:,:) - vb(:,:,:) ) * z1_2dt - ztrdv
         CALL trd_dyn( zua, zva, jpdyn_zdf, kt )
         CALL wrk_dealloc( jpi,jpj,jpk, zua, zva )
         CALL wrk_dealloc( jpi,jpj, zue, zve )
! CR 20161028 --- end ---
! CR 20161031         CALL trd_dyn( ztrdu, ztrdv, jpdyn_zdf, kt )
         CALL wrk_dealloc( jpi, jpj, jpk, ztrdu, ztrdv ) 
      ENDIF
      !                                          ! print mean trends (used for debugging)
      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' zdf  - Ua: ', mask1=umask,               &
            &                    tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
      !
      IF( nn_timing == 1 )  CALL timing_stop('dyn_zdf')
      !
   END SUBROUTINE dyn_zdf

[Robin_Waldman]

I confirm that with time splitting (key_dynspg_ts) and implicit vertical diffusion (ln_zdfexp = .false.), the last fix proposed by Christoph R. allows to close the momentum budget with very high accuracy (error <0.1% on the daily vertically-integrated momentum trend).

However, the update of the barotropic trend should not be performed twice (in dynzdf.F90 and in dynnxt.F90), and the conditions "# if defined key_dynspg_ts" and "IF( l_trddyn )" should be properly added to dyn_zdf routine. See the attachment for an updated version of dynzdf.F90 and dynnxt.F90.