bug in Smagorinsky formulation at 4.0

[davestorkey]

BE CAREFUL !!! Due to dynamic behaviour of this ticket creation page, it is highly recommend to set first all other fields before writing the ticket description below.
If you have lost your draft after an unwanted reload, you can click on the link 'Restore Form' in the contextual menu upper right to recover it.
Remove these lines after reading.

Context

The bilaplacian Smagorinsky formulation at 4.0 gives much smaller viscosity values than the bilaplacian Smagorinsky formulation at 3.6 for same "C" value and the same velocity field in tests with ORCA025 and GYRE. There is about an order of magnitude difference.

Analysis

The main issue appears to be a missing factor of 2 in the definition of esqt and esqf in ldf_dyn_init:

      DO jj = 1, jpj             ! Set local gridscale values
        DO ji = 1, jpi
           esqt(ji,jj) = ( e1e2t(ji,jj) / ( e1t(ji,jj) + e2t(ji,jj) ) )**2 
           esqf(ji,jj) = ( e1e2f(ji,jj) / ( e1f(ji,jj) + e2f(ji,jj) ) )**2      
        END DO
      END DO

instead should be:

      DO jj = 1, jpj             ! Set local gridscale values
        DO ji = 1, jpi
           esqt(ji,jj) = ( 2._wp * e1e2t(ji,jj) / ( e1t(ji,jj) + e2t(ji,jj) ) )**2 
           esqf(ji,jj) = ( 2._wp * e1e2f(ji,jj) / ( e1f(ji,jj) + e2f(ji,jj) ) )**2      
        END DO
      END DO

Then for the special case e1t = e2t = dx we have esqt = dx2 and esqt corresponds to Delta2 in the Griffies and Hallberg (2000) MWR paper. This length scale is raised to the power of 4 in the definition of the bilaplacian viscosity, so a missing factor of 2 would give a bilaplacian viscosity 16x smaller at 4.0 than at 3.6 which is roughly what I see in my tests.

Fix

I have tested this change and the values at 3.6 and 4.0 are now the same order of magnitude but it seems that the 4.0 values are still systematically a bit smaller than the 3.6 values.