Anisotropic tracer diffusion

[francesco.tuccia]

Hello everyone,
I need to implement a particular tracer diffusion scheme, with a fully anisotropic and location dependent diffusion tensor. Right now I coded it with a “brute force” approach, coding it from scratch following exactly what the equations suggests. However, when I look at the nemo implementations of the diffusion operators, they are far more complex so I assume there are more details to take care of. I need some guidance, to know if there is already some routine in Nemo performing such an operation that I can exploit or easily modify, and what are the biggest topics about implementation of diffusion in numerical schemes (if someone have ever written a paper about the different approaches and schemes).

Thank you in advance
Francesco

[jlesommer]

Hi @francesco.tuccia : quite a while ago Nicolas Jourdain and I spent some time implementing anisotropic diffusion dependent on velocity gradients. We may be able to dig this old code out if that’s useful. As least this could give ideas as to the different discretization options.

[francesco.tuccia]

Dear @jlesommer
First of all thank you for the answer. If it is ok for you and your collaborator I would be glad to take a look on that method to have an idea on where to start.
Thank you in advance

[francesco.tuccia]

Dear all,
I’m still developing this method and potentially expanding it to be a diffusion operator for momentum as well. Right now I am following some guidelines: the scheme must be conservative, so the divergence (of the diffusive flux) is implemented in such a way that its integral over the domain is actually equal to the flux balance at the boundary. The inner gradient, characterising the diffusive flux, is implemented as a forward derivative. The components of the diffusion tensor furthermore are located at specific points: diagonal components are at T-points and extra diagonal components at “mixed points”. This choice is arbitrary, but can be explained briefly taking into account what the divergence of a tensor represents: a velocity. This can be justified both dimensionally and physically. So all the components of the diffusion tensor are located in such a way that a (column-wise) divergence moves everything to the corresponding velocity point. For the horizontal diffusion operator, it is all summarised in the attached picture.

Every suggestion and contribution (especially regarding the order of the interpolation and/or their minimisation in number) is more than welcome.

I will keep updating this discussion while I finalise this work
Francesco

[francesco.tuccia]

Using the following convention for the naming of variables:

ztuu = a_{11}\partial_{x}\theta
ztuv = a_{12}\partial_{y}\theta
ztvu = a_{21}\partial_{x}\theta
ztvv = a_{22}\partial_{y}\theta

so that the diffusive flux is

q_{1} = a_{11}\partial_{x}\theta + a_{12}\partial_{y}\theta = ztuu + ztuv
q_{2} = a_{21}\partial_{x}\theta + a_{22}\partial_{y}\theta = ztvu + ztvv,

the code to perform the diffusion of the tracer that I have developed is the following:

      DO jk = 1, jpkm1              !==  Interpolation of diffusion tensor  ==!
         DO jj = 2, jpjm1
            DO ji =  2, jpim1
               int_var11(ji,jj,jk) = ( a_ten(ji+1,jj  ,jk,ia11) + a_ten(ji  ,jj  ,jk,ia11) ) * 0.5_wp
               int_var12(ji,jj,jk) = ( a_ten(ji  ,jj  ,jk,ia12) + a_ten(ji-1,jj  ,jk,ia12) ) * 0.5_wp
               int_var21(ji,jj,jk) = ( a_ten(ji  ,jj  ,jk,ia12) + a_ten(ji  ,jj-1,jk,ia12) ) * 0.5_wp
               int_var22(ji,jj,jk) = ( a_ten(ji  ,jj+1,jk,ia22) + a_ten(ji  ,jj  ,jk,ia22) ) * 0.5_wp
            END DO
         END DO
      END DO  
      !
      ! Lateral boundary condition transfer across nodes
      !
      CALL lbc_lnk_multi( 'new_trahhdiff', int_var11 , 'U', 1.,   &
                        &                  int_var12 , 'V', 1.,   &
                        &                  int_var21 , 'U', 1.,   &
                        &                  int_var22 , 'V', 1.    )
      !
      !
      zw0 = 0._wp
      zw1 = 0._wp
      zw2 = 0._wp
      zw3 = 0._wp
      !
      !                             ! =========== !
      DO jn = 1, jpts               ! tracer loop !
         !                          ! =========== !   
         ztuu = 0._wp 
         ztuv = 0._wp 
         ztvu = 0._wp 
         ztvv = 0._wp 
         !                               
         DO jk = 1, jpkm1              !==  First derivative (gradient)  ==!
            DO jj = 1, jpjm1
               DO ji = 1, jpim1
                  ! Derivative of the tracer
                  ztuu(ji,jj,jk) = r1_e1u(ji,jj) * ( ptb(ji+1,jj  ,jk,jn) - ptb(ji,jj,jk,jn) )
                  ztuv(ji,jj,jk) = r1_e1v(ji,jj) * ( ptb(ji  ,jj+1,jk,jn) - ptb(ji,jj,jk,jn) ) 
                  ztvu(ji,jj,jk) = r1_e1u(ji,jj) * ( ptb(ji+1,jj  ,jk,jn) - ptb(ji,jj,jk,jn) ) 
                  ztvv(ji,jj,jk) = r1_e1v(ji,jj) * ( ptb(ji  ,jj+1,jk,jn) - ptb(ji,jj,jk,jn) )
                  ! Flux
                  ztuu(ji,jj,jk) = e2u(ji,jj) * e3u_n(ji,jj,jk) * int_var11(ji,jj,jk) * ztuu(ji,jj,jk)
                  ztuv(ji,jj,jk) = e1v(ji,jj) * e3v_n(ji,jj,jk) * int_var12(ji,jj,jk) * ztuv(ji,jj,jk)
                  ztvu(ji,jj,jk) = e2u(ji,jj) * e3u_n(ji,jj,jk) * int_var21(ji,jj,jk) * ztvu(ji,jj,jk)
                  ztvv(ji,jj,jk) = e1v(ji,jj) * e3v_n(ji,jj,jk) * int_var22(ji,jj,jk) * ztvv(ji,jj,jk)
                  ! Masking
                  ztuu(ji,jj,jk) = ztuu(ji,jj,jk) * umask(ji,jj,jk) 
                  ztuv(ji,jj,jk) = ztuv(ji,jj,jk) * vmask(ji,jj,jk) 
                  ztvu(ji,jj,jk) = ztvu(ji,jj,jk) * umask(ji,jj,jk) 
                  ztvv(ji,jj,jk) = ztvv(ji,jj,jk) * vmask(ji,jj,jk) 
               END DO
            END DO
         END DO  
         !
         CALL lbc_lnk_multi( 'new_trahhdiff', ztuv, 'F', 1.,   &
                        &                     ztvu, 'F', 1.    )
         !
         DO jk = 1, jpkm1              !==  Second derivative (divergence) added to the general tracer trends  ==!
            DO jj = 2, jpjm1
               DO ji = 2, jpim1
                  !
                  ! Diagonal terms of the diffusion
                  !
                  zw0(ji,jj,jk,jn) = (  ztuu(ji,jj,jk) - ztuu(ji-1,jj  ,jk)     &
                     &               +  ztvv(ji,jj,jk) - ztvv(ji  ,jj-1,jk) )   &
                     &               * tmask(ji,jj,jk) / ( e1e2t(ji,jj) * e3t_n(ji,jj,jk) )
                  !
                  ! Interpolation of extradiagonal term a_12\partial_y\theta
                  !
                  zw1(ji,jj,jk,jn) = 0.25_wp * (  ztuv(ji  ,jj  ,jk) + ztuv(ji+1,jj  ,jk)     &
                     &                         +  ztuv(ji+1,jj-1,jk) + ztuv(ji  ,jj-1,jk) )   &
                     &                         * umask(ji  ,jj  ,jk)
                  !
                  ! Interpolation of extradiagonal term a_21\partial_x\theta
                  !
                  zw2(ji,jj,jk,jn) = 0.25_wp * (  ztvu(ji-1,jj+1,jk) + ztvu(ji  ,jj+1,jk)     &
                     &                         +  ztvu(ji-1,jj  ,jk) + ztvu(ji  ,jj  ,jk) )   &
                     &                         * vmask(ji  ,jj  ,jk)
                  !
               END DO
            END DO
         END DO  
         !
         CALL lbc_lnk_multi( 'new_trahhdiff', zw0(:,:,:,jn), 'T', 1.,   &
                        &                     zw1(:,:,:,jn), 'U', 1.,   &
                        &                     zw2(:,:,:,jn), 'V', 1.    )
         !
         !
         DO jk = 1, jpkm1              !==  Second derivative (divergence) added to the general tracer trends  ==!
            DO jj = 2, jpjm1
               DO ji = 2, jpim1
                  !
                  ! Diagonal terms (0.25 only due to interpolation) 
                  !
                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + zw0(ji,jj,jk,jn) 
                  !
                  ! Extra-diagonal terms
                  !
                  zw3(ji,jj,jk,jn) =  (  zw1(ji,jj,jk,jn) - zw1(ji-1,jj  ,jk,jn)     &
                     &                +  zw2(ji,jj,jk,jn) - zw2(ji  ,jj-1,jk,jn) )   &
                     &               * tmask(ji,jj,jk) / ( e1e2t(ji,jj) * e3t_n(ji,jj,jk) )
                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + zw3(ji,jj,jk,jn)
               END DO
            END DO
         END DO  
         !                          ! ==================
      END DO                        ! end of tracer loop
      !                             ! ==================

I hope that the code is sufficiently self-documented, readable and understandable. Whoever wants to give me suggestions is welcome to do so. I’ll keep updating this post with the other components.
Francesco

[francesco.tuccia]

Using the same conventions explained in my previous post, this is the first draft for the horizontal anisotropic diffusion of momentum. The subroutine implementing the procedure requires the velocity component pcn and the “kind” of velocity kcmp as input.

      SELECT CASE( kcmp )

         !> @par               U Velocity
         CASE ( np_ucmp ) !For U velocity
 
            zwuu = 0._wp 
            zwuv = 0._wp 
            zwvu = 0._wp 
            zwvv = 0._wp 
            !                               
            DO jk = 1, jpkm1              !==  First derivative (gradient)  ==!
               DO jj = 2, jpjm1
                  DO ji = 2, jpim1
                     ! Derivative of the velocity
                     zwuu(ji,jj,jk) = r1_e1t(ji,jj) * ( pcn(ji  ,jj  ,jk) - pcn(ji-1,jj,jk) ) 
                     zwuv(ji,jj,jk) = r1_e1f(ji,jj) * ( pcn(ji  ,jj+1,jk) - pcn(ji  ,jj,jk) )
                     zwvu(ji,jj,jk) = r1_e1t(ji,jj) * ( pcn(ji  ,jj  ,jk) - pcn(ji-1,jj,jk) )
                     zwvv(ji,jj,jk) = r1_e1f(ji,jj) * ( pcn(ji  ,jj+1,jk) - pcn(ji  ,jj,jk) ) 
                     ! Diffusive velocity in x directions
                     zwuu(ji,jj,jk) = a_ten(ji,jj,jk,ia11) * zwuu(ji,jj,jk) * e2t(ji,jj) * e3t_n(ji,jj,jk) 
                     zwuv(ji,jj,jk) = a_ten(ji,jj,jk,ia12) * zwuv(ji,jj,jk)
                     ! Masking
                     zwuu(ji,jj,jk) = zwuu(ji,jj,jk) * tmask(ji,jj,jk) 
                     zwuv(ji,jj,jk) = zwuv(ji,jj,jk) * fmask(ji,jj,jk) 
                     zwvu(ji,jj,jk) = zwvu(ji,jj,jk) * tmask(ji,jj,jk) 
                     zwvv(ji,jj,jk) = zwvv(ji,jj,jk) * fmask(ji,jj,jk) 
                  END DO
               END DO
            END DO             
            !
            CALL lbc_lnk_multi( 'new_hhdiff', zwuu, 'T', 1.,   &
                           &                  zwuv, 'F', 1.,   &
                           &                  zwvu, 'T', 1.,   &
                           &                  zwvv, 'F', 1.    )
            !
            DO jk = 1, jpkm1              !==  First derivative (gradient)  ==!
               DO jj = 1, jpjm1
                  DO ji = 1, jpim1
                     !
                     ! Interpolation of zwuv from f-point to T-point
                     !
                     zw0(ji,jj,jk) = 0.25_wp * ( zwuv(ji-1,jj  ,jk) + zwuv(ji  ,jj  ,jk)     &
                     &                         + zwuv(ji-1,jj-1,jk) + zwuv(ji  ,jj-1,jk) )   &
                     &                       *  tmask(ji  ,jj  ,jk) * e2t(ji,jj) * e3t_n(ji,jj,jk)
                     !
                     ! Interpolation of \partial_{x}u from T-point to f-point and multiplication by a_{12} 
                     !
                     zw1(ji,jj,jk) = 0.25_wp * ( zwvu(ji  ,jj+1,jk) + zwvu(ji+1,jj+1,jk)     &
                     &                         + zwvu(ji  ,jj  ,jk) + zwvu(ji+1,jj  ,jk) )   &
                     &                         *  e1f(ji,jj)         * e3f_n(ji,jj,jk)       &
                     &                      * var_ten(ji,jj,jk,ia12) * fmask(ji,jj,jk) 
                     !
                     ! Interpolation of a_{22} from T-point to f-point and multiplication by \partial_{y}u
                     !
                     zw2(ji,jj,jk) = 0.25_wp * ( a_ten(ji  ,jj+1,jk,ia22) + a_ten(ji+1,jj+1,jk,ia22)     &
                        &                      + a_ten(ji  ,jj  ,jk,ia22) + a_ten(ji+1,jj  ,jk,ia22) )   &
                        &                      *     e1f(ji,jj)    * e3f_n(ji,jj,jk)       &
                        &                         * zwvv(ji,jj,jk) * fmask(ji,jj,jk)
                     !
                  END DO
               END DO
            END DO  
            !
            CALL lbc_lnk_multi( 'new_hhdiff', zw0, 'T', 1.,   &
                           &                  zw1, 'F', 1.,   &
                           &                  zw2, 'F', 1.    )
            !
            !
            !
            DO jk = 1, jpkm1              !==  Second derivative (divergence) added to the general tracer trends  ==!
               DO jj = 2, jpjm1
                  DO ji = 2, jpim1
                     pca(ji,jj,jk) = pca(ji,jj,jk) +                                                &
                        &                            (  zwuu(ji+1,jj  ,jk) - zwuu(ji  ,jj  ,jk)     &
                        &                            +   zw0(ji+1,jj  ,jk) -  zw0(ji  ,jj  ,jk)     &
                        &                            +   zw1(ji  ,jj  ,jk) -  zw1(ji-1,jj  ,jk)     &
                        &                            +   zw2(ji  ,jj  ,jk) -  zw2(ji-1,jj  ,jk) )   &
                        &                            * umask(ji,jj,jk) / ( e1e2u(ji,jj) * e3u_n(ji,jj,jk) )
                     !
                  END DO
               END DO
            END DO 
 
         !> @par               V Velocity
         CASE ( np_vcmp ) !For V velocity 
           
            zwuu = 0._wp 
            zwuv = 0._wp 
            zwvu = 0._wp 
            zwvv = 0._wp 
            !                               
            DO jk = 1, jpkm1              !==  First derivative (gradient)  ==!
               DO jj = 1, jpjm1
                  DO ji = 1, jpim1
                     ! Derivative of the velocity
                     zwuu(ji,jj,jk) = r1_e1f(ji,jj) * ( pcn(ji+1,jj,jk) - pcn(ji,jj  ,jk) ) 
                     zwuv(ji,jj,jk) = r1_e1t(ji,jj) * ( pcn(ji  ,jj,jk) - pcn(ji,jj-1,jk) )
                     zwvu(ji,jj,jk) = r1_e1f(ji,jj) * ( pcn(ji+1,jj,jk) - pcn(ji,jj  ,jk) ) 
                     zwvv(ji,jj,jk) = r1_e1t(ji,jj) * ( pcn(ji  ,jj,jk) - pcn(ji,jj-1,jk) ) 
                     ! Diffusive velocity in x directions
                     zwvu(ji,jj,jk) = a_ten(ji,jj,jk,ia12) * zwvu(ji,jj,jk) 
                     zwvv(ji,jj,jk) = a_ten(ji,jj,jk,ia22) * zwvv(ji,jj,jk) * e1f(ji,jj) * e3f_n(ji,jj,jk)
                     ! Masking
                     zwuu(ji,jj,jk) = zwuu(ji,jj,jk) * fmask(ji,jj,jk) 
                     zwuv(ji,jj,jk) = zwuv(ji,jj,jk) * tmask(ji,jj,jk) 
                     zwvu(ji,jj,jk) = zwvu(ji,jj,jk) * fmask(ji,jj,jk) 
                     zwvv(ji,jj,jk) = zwvv(ji,jj,jk) * tmask(ji,jj,jk) 
                  END DO
               END DO
            END DO             
            !
            CALL lbc_lnk_multi( 'new_hhdiff', zwuu, 'F', 1.,   &
                           &                  zwuv, 'T', 1.,   &
                           &                  zwvu, 'F', 1.,   &
                           &                  zwvv, 'T', 1.    )
            !
            DO jk = 1, jpkm1              !==  First derivative (gradient)  ==!
               DO jj = 2, jpjm1
                  DO ji = 2, jpim1
                     !
                     ! Interpolation of a_{11} from T-point to f-point and multiplication by \partial_{x}v
                     !
                     zw0(ji,jj,jk) = 0.25_wp * ( a_ten(ji  ,jj+1,jk,ia11) + a_ten(ji+1,jj+1,jk,ia11)     &
                     &                         + a_ten(ji  ,jj  ,jk,ia11) + a_ten(ji+1,jj  ,jk,ia11) )   &
                     &                         *     e1f(ji,jj)    * e3f_n(ji,jj,jk)       &
                     &                         *    zwuu(ji,jj,jk) * fmask(ji,jj,jk)
                     !
                     ! Interpolation of \partial_{y}v from T-point to f-point and multiplication by a_{12} 
                     !
                     zw1(ji,jj,jk) = 0.25_wp * ( zwuv(ji  ,jj+1,jk) + zwuv(ji+1,jj+1,jk)     &
                     &                         + zwuv(ji  ,jj  ,jk) + zwuv(ji+1,jj  ,jk) )   &
                     &                         *  e1f(ji,jj)         * e3f_n(ji,jj,jk)       &
                     &                      * a_ten(ji,jj,jk,ia12) * fmask(ji,jj,jk) 
                     !
                     ! Interpolation of zwvu from f-point to T-point
                     !
                     zw2(ji,jj,jk) = 0.25_wp * ( zwvu(ji-1,jj  ,jk) + zwvu(ji  ,jj  ,jk)     &
                     &                         + zwvu(ji-1,jj-1,jk) + zwvu(ji  ,jj-1,jk) )   &
                     &                       *  tmask(ji  ,jj  ,jk) * e2t(ji,jj) * e3t_n(ji,jj,jk)
                     !
                  END DO
               END DO
            END DO  
            !
            CALL lbc_lnk_multi( 'new_hhdiff', zw0, 'F', 1.,   &
                           &                  zw1, 'F', 1.,   &
                           &                  zw2, 'T', 1.    )
            !
            !
            !
            DO jk = 1, jpkm1              !==  Second derivative (divergence) added to the general tracer trends  ==!
               DO jj = 2, jpjm1
                  DO ji = 2, jpim1
                     pca(ji,jj,jk) = pca(ji,jj,jk) +                                                &
                        &                            (   zw0(ji  ,jj  ,jk) -  zw0(ji-1,jj  ,jk)     &
                        &                            +   zw1(ji  ,jj  ,jk) -  zw1(ji-1,jj  ,jk)     &
                        &                            +   zw2(ji+1,jj  ,jk) -  zw2(ji  ,jj  ,jk)     &
                        &                            +  zwvv(ji+1,jj  ,jk) - zwvv(ji  ,jj  ,jk) )   &
                        &                            * vmask(ji,jj,jk) / ( e1e2v(ji,jj) * e3v_n(ji,jj,jk) )
                     !
                  END DO
               END DO
            END DO 
 
         CASE DEFAULT          
           CALL ctl_stop('STOP','new_hhdiff: wrong value for kcmp'  )
      END SELECT

On a side note, all these codes are not optimised, right now I am more concerned about having a functioning implementation rather than a fast and optimised implementation.
Stay tuned for purely vertical diffusion
Francesco

[francesco.tuccia]

Vertical diffusion of tracer:

      int_var33 = 0._wp

      DO jk = 2, jpk               !==  Interpolation of diffusion tensor  ==!  
         DO jj = 1, jpj
            DO ji = 1, jpi
               int_var33(ji,jj,jk) = ( a_ten(ji,jj,jk  ,ia33)               & ! average()
                                   & + a_ten(ji,jj,jk-1,ia33) ) * 0.5_wp
            END DO
         END DO
      END DO  

      DO jn = 1, jpts              !==  loop over the tracers  ==!
         !
         ztww = 0._wp   
         ! Surface boundary fluxes imposed through a specific subroutine,
         ! that for now just imposes ztww(:,:,1) = 0._wp
         CALL new_tradflux( kt, jn, ptb(:,:,1,jn), ztww(:,:,1)  )                 
         !
         DO jk = 2, jpk             !==  First derivative (gradient)  ==!
            DO jj = 1, jpj
               DO ji = 1, jpi
                  ztww(ji,jj,jk) = int_var33(ji,jj,jk) * ( ptb(ji,jj,jk-1,jn) - ptb(ji,jj,jk,jn) )   &
                      &                                / ( e3w_n(ji,jj,jk) )
               END DO
            END DO
         END DO
         !
         ztww(:,:,jpk) = 0._wp      ! Avoid bottom boudary fluxes
         
         DO jk = 1, jpkm1           !==  Second derivative (divergence) added to the general tracer trends  ==!
            DO jj = 1, jpj
               DO ji = 1, jpi
                  !
                  ! Vertical term of the diffusion
                  !
                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + (  ztww(ji,jj,jk) - ztww(ji,jj,jk+1) ) / ( e3t_n(ji,jj,jk) )
                  !
               END DO
            END DO
         END DO

      END DO

[francesco.tuccia]

And vertical diffusion of momentum

      ! Fluxes
      zwww = 0._wp

      SELECT CASE( kcmp )

         !> @par               U Velocity
         CASE ( np_ucmp ) !For U velocity

            ! construction of the vertical flux by interpolation of the variance tensor in 
            ! z and x and multiplication by the derivative of velocity in z
            DO jk = 2, jpkm1
               DO jj = 1, jpjm1
                  DO ji = 1, jpi   
                     zwww(ji,jj,jk) = ( a_ten(ji+1,jj,jk-1,ia33) + a_ten(ji  ,jj,jk-1,ia33)   +           &
                     &                  a_ten(ji+1,jj,jk  ,ia33) + a_ten(ji  ,jj,jk  ,ia33) ) * 0.25_wp * &
                     &                (     pcn(ji  ,jj,jk-1)          - pcn(ji  ,jj,jk  ) * wumask(ji,jj,jk) / e3uw_n(ji,jj,jk)
                  END DO
               END DO
            END DO
            !
            zwww(:,:,jpk) = 0._wp      ! Avoid bottom boudary fluxes
            !
            !                                ! ================
            DO jk = 1, jpkm1                 ! Horizontal slab
               !                             ! ================
               DO jj = 1, jpj
                  DO ji = 1, jpim1
                     ! Averaging and update of the velocity component
                     ! This is a diffusion component added to the RIGHT-HAND SIDE, so positive
                     pca(ji,jj,jk) = pca(ji,jj,jk) + umask(ji,jj,jk) *  &
                     &             ( zwww(ji,jj,jk-1) - zwww(ji,jj,jk) ) / e3uw_n(ji,jj,jk)
                  END DO  
               END DO
               !                             ! ================
            END DO                           !   End of slab
            !                                ! ================


         !> @par               V Velocity
         CASE ( np_vcmp ) !For V velocity 

            ! construction of the vertical flux by interpolation of the variance tensor in 
            ! z and y and multiplication by the derivative of velocity in z
            DO jk = 2, jpkm1
               DO jj = 1, jpj
                  DO ji = 1, jpim1 
                     zwww(ji,jj,jk) = ( a_ten(ji,jj+1,jk,ia33) + a_ten(ji,jj  ,jk,ia33) +               &
                     &                  a_ten(ji,jj+1,jk,ia33) + a_ten(ji,jj  ,jk,ia33) ) * 0.25_wp *   &
                     &                (     pcn(ji,jj  ,jk-1)        - pcn(ji,jj  ,jk  ) * wvmask(ji,jj,jk) / e3vw_n(ji,jj,jk)
                  END DO
               END DO
            END DO
            !
            zwww(:,:,jpk) = 0._wp      ! Avoid bottom boudary fluxes
            !
            !                                ! ================
            DO jk = 2, jpkm1                 ! Horizontal slab
               !                             ! ================
               DO jj = 1, jpjm1
                  DO ji = 1, jpi
                     ! Averaging and update of the velocity component
                     ! This is a diffusion component added to the RIGHT-HAND SIDE, so positive
                     pca(ji,jj,jk) = pca(ji,jj,jk) + vmask(ji,jj,jk) *  &
                     &             ( zwww(ji,jj,jk-1) - zwww(ji,jj,jk) ) / e3uw_n(ji,jj,jk)
                  END DO  
               END DO
               !                             ! ================
            END DO                           !   End of slab
            !                                ! ================
            ! 

         CASE DEFAULT          
           CALL ctl_stop('STOP','new_zzdiff: wrong value for kcmp'  )
      END SELECT

These last two codes are the most obscure. Going through the NEMO documentation I’ve noticed that vertical diffusion is imposed implicitly, but I can’t find more about the schemes. Can someone give me some hint?