QG Leith viscosity implementation

Configurations Global
, , ,
1

Hi all,

I am currently attempting to implement the QG Leith viscous coefficient in nemo v4.0.4. The choice of v4.0.4 is because there is an ORCA025 setup ready to go, which is not yet the case for v4.2. I am also new to the nemo ocean model, but am familiar with the MITgcm. The QG Leith has been implemented in MITgcm, so I intend to lean heavily on this.

In the images below is a latex document detailing QG Leith and some of the nemo basics. There are also some questions regarding the calculation of some terms and how the lateral boundary condition function works.

QG_leith_pg1
QG_leith_pg1501×782 24.8 KB
QG_Leith_pg2
QG_Leith_pg2505×683 24.5 KB

In addition, I have made a start on coding up this scheme, however, I will update this thread with my progress at a later date.

In the mean time, I would be grateful for any tips on this work, in particular on the calculation of some terms (vorticity, divergence), as well as the lateral boundary conditions function.

Thanks

Bachman et al. (2017):
https://agupubs.onlinelibrary.wiley.com/doi/full/10.1002/2016JC012265

NEMO v4.0.4
https://forge.ipsl.jussieu.fr/nemo/svn/NEMO/releases/r4.0/r4.0.4/

2

There is a first set or remarks to start…

In your equation 2:

  • For a gradient at F point, you have
    • the zonal component is at V point: ( q2(ji,jj) - q2(ji-1,jj) ) / e1v(ji,jj)
    • the meridional component is at U point: ( q2(ji,jj) - q2(ji,jj-1) ) / e2u(ji,jj)
  • For a gradient at T point, you have
    • the zonal component is at U point: ( b(ji+1,jj) - b(ji,jj) ) / e1u(ji,jj)
    • the meridional component is at V point: ( b(ji,jj+1) - b(ji,jj) ) / e2v(ji,jj)

→ you cannot directly add these 2 gradient as their zonal and meridional components are not located at the same point. I guess you should first move them to the T point with a simple average: 0.5*(u(ji,jj)+u(ji-1,jj)) or 0.5*(v(ji,jj)+v(ji,jj-1))

Question 1:
dynvor is computing other things that the current curl. You need only the relative vorticity part.
I think you should recompute it (using un and not ub) with this :

       DO jk = 1, jpkm1                                 ! Horizontal slab
          DO jj = 1, jpjm1
             DO ji = 1, fs_jpim1   ! vector opt.
                zwz(ji,jj) = (  e2v(ji+1,jj  ) * un(ji+1,jj  ,jk) - e2v(ji,jj) * vn(ji,jj,jk)    &
                   &          - e1u(ji  ,jj+1) * un(ji  ,jj+1,jk) + e1u(ji,jj) * un(ji,jj,jk)  ) * r1_e1e2f(ji,jj)
             END DO
          END DO
       END DO

Question 2:
in ldfdyn, dshe is not put through the lateral boundary condition function because it is not needed!
dshesq is computed from 1 to jpi-1 and 1 to jpj-1 and is not use in jpi and jpj in the following loop (when you compute aht or ahmf)

3

Thanks for these initial comments. Below I have made some progress on implementing the QG Leith code.

The first step is to calculate buoyancy at T points from points jk=1:jpkm1 and then putting this on W points, which is the same as the buoyancy frequency, pn2. Buoyancy is b = -g(\rho - \rho_0)/\rho_0 . Since we take the horizontal gradient of buoyancy, we can use the density anomaly.

DO jk = 1, jpkm1
   DO jj = 1, jpjm1
      DO ji = 1, jpim1
         IF( jk < 2 )
            ! buoyancy at surface on w-grid !
            zbu(ji,jj,jk) = - grav * prd(ji,jj,jk)
         ELSE
            ! buoyancy below surface on w-grid !
            zbuup = - grav * prd(ji,jj,jk-1)
            zbulw = - grav * prd(ji,jj,jk  )
            zbu(ji,jj,jk) = 0.5_wp * ( zbuup + zbulw )
         ENDIF
      END DO
   END DO
END DO

Then, we find (f/N^2)b and take its vertical derivative to form zst

DO jk = 1, jpkm1
   DO jj = 1, jpjm1
      DO ji = 1, jpim1
         IF( jk < jpkm1 )
            zker1 = ( ff_t(ji,jj) / MAX( pn2(ji,jj,jk  ), zqglep ) ) * zbu(ji,jj,jk  )
            zker2 = ( ff_t(ji,jj) / MAX( pn2(ji,jj,jk+1), zqglep ) ) * zbu(ji,jj,jk+1)
            zst(ji,jj,jk) = r1_e3w(ji,jj,jk) * ( zker1 - zker2 ) * tmask(ji,jj,jk)
         ELSE
            zker1 = ( ff_t(ji,jj) / MAX( pn2(ji,jj,jk  ), zqglep ) ) * zbu(ji,jj,jk  )
            zker2 = 0._wp ! no value set at ocean bottom jk = jpk as it lies outside domain !
            zst(ji,jj,jk) = r1_e3w(ji,jj,jk) * ( zker1 - zker2 ) * tmask(ji,jj,jk)
         ENDIF   
      END DO
   END DO
END DO

Since the bottom T point is outside the domain, I have taken zker2=0._wp at jpk. In zker we find the max of pn2 and a small number to avoid division by zero.

These routines calculate the term \partial_z (f/N^2)b.

4

Following on, we then calculate the absolute vorticity (Coriolis + relative vorticity) on F point, which is actually q_2 in the initial equation, Eq (1).

DO jk = 1, jpkm1                                 ! Horizontal slab
   DO jj = 1, jpjm1
      DO ji = 1, fs_jpim1   ! vector opt.
         zwz(ji,jj,jk) = ff_f(ji,jj) + (  e2v(ji+1,jj  ) * vn(ji+1,jj  ,jk) - e2v(ji,jj) * vn(ji,jj,jk)            &
            &          - e1u(ji  ,jj+1) * un(ji  ,jj+1,jk) + e1u(ji,jj) * un(ji,jj,jk)  ) * r1_e1e2f(ji,jj)
      END DO
   END DO
END DO

We then find the square of the Rossby and Froude numbers, and limit the stretching term (f/N^2)b if needed based on these dimensionless numbers. This is on the F point

DO jk = 1, jpkm1
   DO jj = 1, jpjm1
      DO ji = 1, jpim1
         zusq = r1_2 * ( ( un(ji,jj,jk)**2 + un(ji,jj+1,jk)**2 ) + ( vn(ji,jj,jk)**2 + vn(ji+1,jj,jk)**2 ) )
         zrosq = ( zusq / ( MAX( ff_f(ji,jj), zQGLep ) * ( eqsf(ji,jj)**2 ) ) ) * fmask(ji,jj,jk)
         znsq = r1_2 * ( r1_4 * ( pn2(ji,jj,jk) + pn2(ji+1,jj,jk) + pn2(ji,jj+1,jk) + pn2(ji+1,jj+1,jk) )          &
            &  + r1_4 * ( pn2(ji,jj,jk+1) + pn2(ji+1,jj,jk+1) + pn2(ji,jj+1,jk+1) + pn2(ji+1,jj+1,jk+1) ) )         ! pn2 on f-point !
         zfrsq = ( ( zusq * PI * PI ) / MAX( ( SQRT( znsq ) * e3f(ji,jj) )**2, zQGLep ) ) * fmask(ji,jj,jk)
         ! limit the stretching term !
         zsthld = min( ABS( r1_4 * ( zst(ji,jj,jk) + zst(ji+1,jj,jk) + zst(ji,jj+1,jk) + zst(ji+1,jj+1,jk) ) ),    &
            &  ABS( ( zwz(ji,jj,jk) * zfrsq ) /                                                                    &
            &  ( zrosq + zfrsq**2 + zQGLep ) ) )
         zst(ji,jj,jk) = SIGN( zsthld, zst(ji,jj,jk) ) ! possibly define key_nosignedzero? !
      END DO
   END DO
END DO

It doesn’t seem like using SIGN will provide any issues as zst is unlikely to be zero.

5

With absolute vorticity q_2 = f+\zeta and stretching \partial_z (f/N^2)b, their combined horizontal gradients and magnitude are found on the T point

!== calculate gradients of vorticity and stretching, then square of magnitude (t-point) ==!
DO jk = 1, jpkm1
   DO jj = 2, jpjm1
      DO ji = 2, jpim1
         zvsdx = r1_2 * ( r1_e1f(ji,jj) * ( zst(ji,jj-1,jk) + zwz(ji,jj-1,jk) -                &
            &              zst(ji-1,jj-1,jk) - zwz(ji-1,jj-1,jk) )                             &
            &                 + r1_e1f(ji,jj) * ( zst(ji,jj,jk) + zwz(ji,jj,jk) -              &
            &                     zst(ji-1,jj,jk) - zwz(ji-1,jj,jk) ) )
         zvsdy = r1_2 * ( r1_e2f(ji,jj) * ( zst(ji-1,jj,jk) + zwz(ji-1,jj,jk) -                &
            &              zst(ji-1,jj-1,jk) - zwz(ji-1,jj-1,jk) )                             &
            &                 + r1_e2f(ji,jj) * ( zst(ji,jj,jk) + zwz(ji,jj,jk) -              &
            &                     zst(ji,jj-1,jk) - zwz(ji,jj-1,jk) ) ) 
         dvsmagsq(ji,jj,jk) = ( zvsdx**2 + zvsdy**2 ) * tmask(ji,jj,jk)
      END DO
   END DO
END DO
!
CALL lbc_lnk_multi( 'ldfdyn', dvsmagsq, 'T', 1. )

and dvsmagsq is applied through the lateral boundary condition.

Then, gradients of divergence and its magnitude are found on the F point. Horizontal divergence is calculated in the module divhor.F90 as hdivn

!== calculate gradients of divergence, then square of magnitude (f-point) ==!
DO jk = 1, jpkm1
   DO jj = 1, jpjm1
      DO ji = 1, jpim1
         zdivdx = r1_2 * ( r1_e1t(ji,jj) * ( hdivn(ji+1,jj+1,jk) - hdivn(ji,jj+1,jk) )               &
            &  + r1_e1t(ji,jj) * ( hdivn(ji+1,jj,jk) - hdivn(ji,jj,jk) ) )
         zdivdy = r1_2 * ( r1_e2t(ji,jj) * ( hdivn(ji+1,jj+1,jk) - hdivn(ji+1,jj,jk) )               &
            &  + r1_e2t(ji,jj) * ( hdivn(ji,jj+1,jk) - hdivn(ji,jj,jk) ) )
         ddivmagsq(ji,jj,jk) = ( zdivdx**2 + zdivdy**2 ) * fmask(ji,jj,jk)
      END DO
   END DO
END DO

We now have nearly everything to calculate the viscous coefficient A_{qg}, but first we find (\Lambda/\pi)^6 (n.b. this is later square rooted)

zcmqgl  = (rn_cqgc/rpi)**6         			! (C_qg/pi)^6

The viscous coefficients are found

DO jk = 1, jpkm1	         !== QG Leith viscosity coefficient on T-point ==!
   DO jj = 2, jpjm1
      DO ji = fs_2, fs_jpim1 ! vector opt.
         zsqQG = r1_4 * ( ddivmag(ji,jj,jk) + ddivmag(ji-1,jj,jk) + ddivmag(ji,jj-1,jk) +     &
            &  ddivmag(ji-1,jj-1,jk) ) + dvsmag(ji,jj,jk)
         amht(ji,jj,jk) = SQRT( zcmqgl * esqt(ji,jj)**3 * zsqQG )
      END DO
   END DO
END DO
!
DO jk = 1, jpkm1            !== QG Leith viscosity coefficient on F-point ==!
   DO jj = 1, jpjm1
      DO ji = 1, fs_jpim1 ! vector opt.
         zsqQG = r1_4 * ( dvsmag(ji,jj,jk) + dvsmag(ji+1,jj,jk) + dvsmag(ji,jj+1,jk) +     &
            &  dvsmag(ji+1,jj+1,jk) ) + ddivmag(ji,jj,jk)
         amhf(ji,jj,jk) = SQRT( zcmqgl * esqf(ji,jj)**3 * zsqQG )
      END DO
   END DO
END DO

And then amht and amhf are passed through lbc_lnk().

I believe this to be almost complete. As always, I would appreciate any further tips if possible e.g. optimisation. For example, is it worth calculating Ro and Fr numbers every n time steps?

7

Hi,

I wasn’t sure whether to create a new post for this issue, or continue here on this thread. I’ll continue here for now.

Having successfully compiled the nemo model with updated ldfdyn.F90 and step.F90, I have attempted to run a simulation using the QG Leith code. However, the model crashes and I receive the following errors in testing.output:
*** glibc detected *** ./nemo: free(): invalid pointer
I think the above error implies the code is trying to access an index that does not exist. I can’t see anything immediately obvious with the below routine,

CASE( 33 )       !==  time varying 3D field  ==!   = F(QG PV gradient, divergence, and gridscale) (QG Leith)
!
IF( ln_dynldf_lap ) THEN        ! laplacian operator
   !
   ! allocate local variables !
   ALLOCATE( zbu(jpi,jpj,jpk) , zst(jpi,jpj,jpk) , zwz(jpi,jpj,jpk) )
   !
   zcmqgl = (rn_cqgc/rpi)**6         			! (C_qg/pi)^6
   zqglep1 = 1.e-12_wp
   zqglep2 = 1.e-24_wp
   !== begin calculation of stretching term d/dz[(f/(N**2))*b] ==!
   ! find buoyancy and interpolate onto w-grid !
   DO jk = 1, jpkm1
      DO jj = 1, jpjm1
         DO ji = 1, jpim1
            IF( jk < 2 ) THEN
               ! buoyancy at surface on w-grid !
               zbu(ji,jj,jk) = - grav * prd(ji,jj,jk)
            ELSE
               ! buoyancy below surface on w-grid !
               zbuup = - grav * prd(ji,jj,jk-1)
               zbulw = - grav * prd(ji,jj,jk  )
               zbu(ji,jj,jk) = 0.5_wp * ( zbuup + zbulw )
            ENDIF
         END DO
      END DO
   END DO
   !
   !== finish calculation and take vertical derivative (t-point) ==!
   DO jk = 1, jpkm1
      DO jj = 1, jpjm1
         DO ji = 1, jpim1
            IF( jk < jpkm1 ) THEN
               zker1 = ( ff_t(ji,jj) / MAX( pn2(ji,jj,jk  ), zqglep1 ) ) * zbu(ji,jj,jk  )
               zker2 = ( ff_t(ji,jj) / MAX( pn2(ji,jj,jk+1), zqglep1 ) ) * zbu(ji,jj,jk+1)
               zst(ji,jj,jk) = (( zker1 - zker2 )  / e3w_n(ji,jj,jk) ) * tmask(ji,jj,jk)
            ELSE
               zker1 = ( ff_t(ji,jj) / MAX( pn2(ji,jj,jk  ), zqglep1 ) ) * zbu(ji,jj,jk  )
               zker2 = 0._wp ! no value set at ocean bottom jk = jpk as it lies outside domain !
               zst(ji,jj,jk) =  (( zker1 - zker2  ) / e3w_n(ji,jj,jk) ) * tmask(ji,jj,jk)
            ENDIF   
         END DO
      END DO
   END DO
   !
   !== calculate absolute vorticity (f+zeta) on f-point ==!
   DO jk = 1, jpkm1                                 ! Horizontal slab
      DO jj = 1, jpjm1
         DO ji = 1, fs_jpim1   ! vector opt.
            zwz(ji,jj,jk) = ff_f(ji,jj) + (  e2v(ji+1,jj  ) * vn(ji+1,jj  ,jk) - e2v(ji,jj) * vn(ji,jj,jk)            &
               &          - e1u(ji  ,jj+1) * un(ji  ,jj+1,jk) + e1u(ji,jj) * un(ji,jj,jk)  ) * r1_e1e2f(ji,jj)
         END DO
      END DO
   END DO
   !
   !== calculate square of Rossby and Froude numbers, then limit stretching term if needed (f-point) ==!
   DO jk = 1, jpkm1
      DO jj = 1, jpjm1
         DO ji = 1, jpim1
            zusq = r1_2 * ( ( un(ji,jj,jk)**2 + un(ji,jj+1,jk)**2 ) + ( vn(ji,jj,jk)**2 + vn(ji+1,jj,jk)**2 ) )
            zrosq = ( zusq / ( MAX( ff_f(ji,jj)**2, zqglep2 ) * ( esqf(ji,jj)**2 ) ) ) * fmask(ji,jj,jk)
            znsq = r1_2 * ( r1_4 * ( pn2(ji,jj,jk) + pn2(ji+1,jj,jk) + pn2(ji,jj+1,jk) + pn2(ji+1,jj+1,jk) )          &
               &  + r1_4 * ( pn2(ji,jj,jk+1) + pn2(ji+1,jj,jk+1) + pn2(ji,jj+1,jk+1) + pn2(ji+1,jj+1,jk+1) ) )         ! pn2 on f-point !
            zfrsq = ( ( zusq * rpi * rpi ) / MAX( ( znsq * e3f_n(ji,jj,jk)**2 ), zqglep2 ) ) * fmask(ji,jj,jk)
            ! limit the stretching term !
            zsthld = min( ABS( r1_4 * ( zst(ji,jj,jk) + zst(ji+1,jj,jk) + zst(ji,jj+1,jk) + zst(ji+1,jj+1,jk) ) ),    &
               &  ABS( ( zwz(ji,jj,jk) * zfrsq ) /                                                                    &
               &  ( zrosq + zfrsq**2 + zqglep2 ) ) )
            zst(ji,jj,jk) = SIGN( zsthld, zst(ji,jj,jk) ) ! possibly define key_nosignedzero? !
         END DO
      END DO
   END DO
   !
   !== calculate gradients of vorticity and stretching, then square of magnitude (t-point) ==!
   DO jk = 1, jpkm1
      DO jj = 2, jpjm1
         DO ji = 2, jpim1
            zvsdx = r1_2 * ( r1_e1f(ji,jj-1) * ( zst(ji,jj-1,jk) + zwz(ji,jj-1,jk) -                &
               &              zst(ji-1,jj-1,jk) - zwz(ji-1,jj-1,jk) )                             &
               &                 + r1_e1f(ji,jj) * ( zst(ji,jj,jk) + zwz(ji,jj,jk) -              &
               &                     zst(ji-1,jj,jk) - zwz(ji-1,jj,jk) ) )
            zvsdy = r1_2 * ( r1_e2f(ji-1,jj) * ( zst(ji-1,jj,jk) + zwz(ji-1,jj,jk) -                &
               &              zst(ji-1,jj-1,jk) - zwz(ji-1,jj-1,jk) )                             &
               &                 + r1_e2f(ji,jj) * ( zst(ji,jj,jk) + zwz(ji,jj,jk) -              &
               &                     zst(ji,jj-1,jk) - zwz(ji,jj-1,jk) ) ) 
            dvsmagsq(ji,jj,jk) = ( zvsdx**2 + zvsdy**2 ) * tmask(ji,jj,jk)
         END DO
      END DO
   END DO
   !
   CALL lbc_lnk_multi( 'ldfdyn', dvsmagsq, 'T', 1. )
   !
   !== calculate gradients of divergence, then square of magnitude (f-point) ==!
   DO jk = 1, jpkm1
      DO jj = 1, jpjm1
         DO ji = 1, jpim1
            zdivdx = r1_2 * ( r1_e1t(ji,jj+1) * ( hdivn(ji+1,jj+1,jk) - hdivn(ji,jj+1,jk) )               &
               &  + r1_e1t(ji,jj) * ( hdivn(ji+1,jj,jk) - hdivn(ji,jj,jk) ) )
            zdivdy = r1_2 * ( r1_e2t(ji+1,jj) * ( hdivn(ji+1,jj+1,jk) - hdivn(ji+1,jj,jk) )               &
               &  + r1_e2t(ji,jj) * ( hdivn(ji,jj+1,jk) - hdivn(ji,jj,jk) ) )
            ddivmagsq(ji,jj,jk) = ( zdivdx**2 + zdivdy**2 ) * fmask(ji,jj,jk)
         END DO
      END DO
   END DO
   !
   DO jk = 1, jpkm1	         !== QG Leith viscosity coefficient on T-point ==!
      DO jj = 2, jpjm1
         DO ji = fs_2, fs_jpim1 ! vector opt.
            zsqqg = r1_4 * ( ddivmagsq(ji,jj,jk) + ddivmagsq(ji-1,jj,jk) + ddivmagsq(ji,jj-1,jk) +     &
               &  ddivmagsq(ji-1,jj-1,jk) ) + dvsmagsq(ji,jj,jk)
            ahmt(ji,jj,jk) = SQRT( zcmqgl * esqt(ji,jj)**3 * zsqqg )
         END DO
      END DO
   END DO
   !
   DO jk = 1, jpkm1            !== QG Leith viscosity coefficient on F-point ==!
      DO jj = 1, jpjm1
         DO ji = 1, fs_jpim1 ! vector opt.
            zsqqg = r1_4 * ( dvsmagsq(ji,jj,jk) + dvsmagsq(ji+1,jj,jk) + dvsmagsq(ji,jj+1,jk) +     &
               &  dvsmagsq(ji+1,jj+1,jk) ) + ddivmagsq(ji,jj,jk)
            ahmf(ji,jj,jk) = SQRT( zcmqgl * esqf(ji,jj)**3 * zsqqg )
         END DO
      END DO
   END DO
   !
   DEALLOCATE( zbu, zst, zwz )
   !
ENDIF
!
CALL lbc_lnk_multi( 'ldfdyn', ahmt, 'T', 1.,  ahmf, 'F', 1. )
!
END SELECT

Part of me wonders whether I should replace un/vn with ub/vb and then calculate hdiv inside the routine using ub/vb.

Anyway, has anyone encountered an error like this before? And how did you fix it?

Cheers

8

Could there error come from a variable that hasn’t been (properly) allocated, or deallocated in the wrong place?

9

The allocation may have been the error. However, earlier in the routine (which I haven’t shown) I had incorrectly entered indices into the allocation of a variable, so ddivmagsq(ji,jj,jk) instead of ddivmagsq(jpi,jpj,jpk). Making these changes, the model now runs with no error… but the viscosity output is not physical. I will go over the maths later with a fine toothpick.