New GM runtime parameter to allow in/excluding dynamically computed K…

…appa contributions in Redi-tensor (MITgcm#769)

* new GM runtime parm. to allow in/excluding dynamic K in Redi-tensor

- eg. for Visbeck
- true by default (current behavior)

* rename GM_varKredi -> GM_isoFac_calcK and make it the runtime parm

* address 3 comments:

- add the new factor GM_isoFac_calcK to the Kwz-coefficient (was missing)
- update setting of GM_ExtraDiag the skew-flux form
- remove GM_advect

* remove leftover GM_useVarKforRedi

* Fix setting of GM_ExtraDiag

- move resetting of GM_useLeithQG to F (if viscC2LeithQG=0) from gmredi_check.F to gmredi_readparms.F
- fix the setting of GM_ExtraDiag (for background 3-D read from file, for QG-Leith, for Bates 3-D K)
  This does matter since, when incorrectly set to F, was missing 2 terms in isopycnal diffusion tensor
- add more check and stop for the case #undef GM_NON_UNITY_DIAGONAL

* refine test on 3-D file names

* cleaning

Remove Viscbeck code inside #ifndef GM_NON_UNITY_DIAGONAL since a check & stop
was added for this case.
Also add 2 lines of comments to justify keeping this block (from PR MITgcm#537)

* update ref. output after fixing bug

Setting of GM_ExtraDiag has been fixed (if using background 3D-K from file,
or QG-Leith or Bates 3-D K) and is now set to True in this test experiment.

* Also account for pkg/ctrl when setting GM_ExtraDiag

for now, just used CPP option ALLOW_KAPREDI_CONTROL and ALLOW_KAPGM_CONTROL

* Document pkg/gmredi addition (new param) and fix

---------

Co-authored-by: Jean-Michel Campin <[email protected]>
Co-authored-by: Jean-Michel Campin <[email protected]>

doc/tag-index

@@ -1,6 +1,13 @@
     Notes on tags used in MITgcmUV
     ==============================
 
+o pkg/gmredi:
+  - fix setting of GM_ExtraDiag that was not (always) turned on with 3D Kappa
+    read from file, with Bates-K3d or QG-Leith, therefore missing 2 terms in
+    isopycnal diffusion tensor; update ref. output of MLAdjust.QGLthGM test ;
+  - add more check and stop for the case #undef GM_NON_UNITY_DIAGONAL ;
+  - add runtime parameter GM_isoFac_calcK to allow to scale dynamically computed
+    Kappa contribution in the Redi-tensor (default = 1.).
 o pkg/autodiff:
   - adjust call to S/R DIAGNOSTICS_SWITCH_ONOFF in autodiff_inadmode_set_ad.F
     so that snap-shot Adjoint diagnostics are filled within the same time-step

pkg/gmredi/GMREDI.h

@@ -12,7 +12,8 @@ C--   COMMON /GM_PARAMS_L/ GM/Redi Logical-type parameters
 C     GM_AdvForm       :: use Advective Form (instead of Skew-Flux form)
 C     GM_AdvSeparate   :: do separately advection by Eulerian and Bolus velocity
 C     GM_useBVP        :: use Boundary-Value-Problem method for Bolus transport
-C     GM_useSubMeso    :: use parameterization of mixed layer (Sub-Mesoscale) eddies
+C     GM_useSubMeso    :: use parameterization of mixed layer (Sub-Mesoscale)
+C                         eddies
 C     GM_ExtraDiag     :: select extra diagnostics
 C     GM_InMomAsStress :: apply GM as a stress in momentum Eq.
 C     GM_MNC           ::
@@ -92,6 +93,8 @@ C     GM_Kmin_horiz    :: minimum horizontal diffusivity [m^2/s]
 C     GM_Small_Number  :: epsilon used in computing the slope
 C     GM_slopeSqCutoff :: slope^2 cut-off value
 C     GM_Scrit, GM_Sd  :: parameter for 'dm95' & 'ldd97' tapering fct
+C     GM_isoFac_calcK  :: add fraction of  dynamically computed variable K,
+C                         e.g. Visbeck, also to Redi tensor (default = 1.)
 C-    Transition layer thickness definition:
 C     GM_facTrL2dz   :: minimum Trans. Layer Thick. as a factor of local dz
 C     GM_facTrL2ML   :: maximum Trans. Layer Thick. as a factor of Mix-Layer Depth
@@ -130,6 +133,7 @@ C     GM_Bates_maxRenorm :: maximum value for the renormalisation factor
       _RL GM_Small_Number
       _RL GM_slopeSqCutoff
       _RL GM_Scrit, GM_Sd
+      _RL GM_isoFac_calcK
       _RL GM_facTrL2dz
       _RL GM_facTrL2ML
       _RL GM_maxTransLay
@@ -167,7 +171,7 @@ C     GM_Bates_maxRenorm :: maximum value for the renormalisation factor
      &                 GM_maxSlope,
      &                 GM_Kmin_horiz,
      &                 GM_Small_Number, GM_slopeSqCutoff,
-     &                 GM_Scrit, GM_Sd,
+     &                 GM_Scrit, GM_Sd, GM_isoFac_calcK,
      &                 GM_facTrL2dz, GM_facTrL2ML, GM_maxTransLay,
      &                 GM_BVP_cMin,
      &                 subMeso_Ceff, subMeso_invTau, subMeso_LfMin,
@@ -191,12 +195,11 @@ C--   COMMON /GM_DERIVED_PAR/ other GM/Redi parameters
 C     (derived from previous block and not directly user configured)
       _RL GM_rMaxSlope
       _RL GM_skewflx
-      _RL GM_advect
       _RL GM_BVP_rModeNumber
       _RL GM_BVP_cHat2Min
       COMMON /GM_DERIVED_PAR/
      &                   GM_rMaxSlope,
-     &                   GM_skewflx, GM_advect,
+     &                   GM_skewflx,
      &                   GM_BVP_rModeNumber, GM_BVP_cHat2Min
 
 C--   COMMON /GM_COEFFICIENTS/ GM/Redi scaling coefficients

pkg/gmredi/gmredi_calc_tensor.F

@@ -568,17 +568,17 @@ SUBROUTINE GMREDI_CALC_TENSOR(
      &           + GM_skewflx*bolus_K*GM_bolFac2d(i,j,bi,bj)
 #endif
 #ifdef GM_VISBECK_VARIABLE_K
-     &           + VisbeckK(i,j,bi,bj)*(1. _d 0 + GM_skewflx)
+     &           + VisbeckK(i,j,bi,bj)*(GM_isoFac_calcK + GM_skewflx)
 #endif
 #ifdef ALLOW_GM_LEITH_QG
      &           + op5*( GM_LeithQG_K(i,j,km1,bi,bj)
      &                 + GM_LeithQG_K(i,j,k,bi,bj) )
-     &                *(1. _d 0 + GM_skewflx)
+     &                *(GM_isoFac_calcK + GM_skewflx)
 #endif
 #if ((defined GM_BATES_K3D) && ! (defined GM_BATES_PASSIVE))
      &           + op5*( GM_BatesK3d(i,j,km1,bi,bj)
      &                 + GM_BatesK3d(i,j,k,bi,bj) )
-     &                *(1. _d 0 + GM_skewflx)
+     &                *(GM_isoFac_calcK + GM_skewflx)
 #endif
          Kwx(i,j,k,bi,bj)= Kgm_tmp*Kwx(i,j,k,bi,bj)
          Kwy(i,j,k,bi,bj)= Kgm_tmp*Kwy(i,j,k,bi,bj)
@@ -589,15 +589,17 @@ SUBROUTINE GMREDI_CALC_TENSOR(
          Kwz(i,j,k,bi,bj)= ( isopycK*GM_isoFac2d(i,j,bi,bj)
 #endif
 #ifdef GM_VISBECK_VARIABLE_K
-     &                     + VisbeckK(i,j,bi,bj)
+     &                     + VisbeckK(i,j,bi,bj)*GM_isoFac_calcK
 #endif
 #ifdef ALLOW_GM_LEITH_QG
      &                     + op5*( GM_LeithQG_K(i,j,km1,bi,bj)
      &                           + GM_LeithQG_K(i,j,k,bi,bj) )
+     &                                          *GM_isoFac_calcK
 #endif
 #if ((defined GM_BATES_K3D) && ! (defined GM_BATES_PASSIVE))
      &                     + op5*( GM_BatesK3d(i,j,km1,bi,bj)
      &                           + GM_BatesK3d(i,j,k,bi,bj) )
+     &                                          *GM_isoFac_calcK
 #endif
      &                     )*Kwz(i,j,k,bi,bj)
         ENDDO
@@ -755,12 +757,15 @@ SUBROUTINE GMREDI_CALC_TENSOR(
 #endif
 #ifdef GM_VISBECK_VARIABLE_K
      &     + op5*(VisbeckK(i-1,j,bi,bj)+VisbeckK(i,j,bi,bj))
+     &          *GM_isoFac_calcK
 #endif
 #ifdef ALLOW_GM_LEITH_QG
      &     + op5*(GM_LeithQG_K(i-1,j,k,bi,bj)+GM_LeithQG_K(i,j,k,bi,bj))
+     &          *GM_isoFac_calcK
 #endif
 #if ((defined GM_BATES_K3D) && ! (defined GM_BATES_PASSIVE))
      &     + op5*(GM_BatesK3d(i-1,j,k,bi,bj)+GM_BatesK3d(i,j,k,bi,bj))
+     &          *GM_isoFac_calcK
 #endif
      &     )*taperFct(i,j)
          ENDDO
@@ -802,15 +807,18 @@ SUBROUTINE GMREDI_CALC_TENSOR(
      &        *op5*(GM_bolFac2d(i-1,j,bi,bj)+GM_bolFac2d(i,j,bi,bj))
 #endif
 #ifdef GM_VISBECK_VARIABLE_K
-     &     + op5*(VisbeckK(i-1,j,bi,bj)+VisbeckK(i,j,bi,bj))*GM_advect
+     &     + op5*(VisbeckK(i-1,j,bi,bj)+VisbeckK(i,j,bi,bj))
+     &          *(GM_isoFac_calcK - GM_skewflx)
 #endif
 #ifdef ALLOW_GM_LEITH_QG
      &     + op5*( GM_LeithQG_K(i-1,j,k,bi,bj)
-     &           + GM_LeithQG_K(i,j,k,bi,bj) )*GM_advect
+     &           + GM_LeithQG_K(i,j,k,bi,bj) )
+     &          *(GM_isoFac_calcK - GM_skewflx)
 #endif
 #if ((defined GM_BATES_K3D) && ! (defined GM_BATES_PASSIVE))
      &     + op5*( GM_BatesK3d(i-1,j,k,bi,bj)
-     &           + GM_BatesK3d(i,j,k,bi,bj) )*GM_advect
+     &           + GM_BatesK3d(i,j,k,bi,bj) )
+     &          *(GM_isoFac_calcK - GM_skewflx)
 #endif
      &     )*SlopeX(i,j)*taperFct(i,j)
          ENDDO
@@ -840,13 +848,16 @@ SUBROUTINE GMREDI_CALC_TENSOR(
 #endif
 #ifdef GM_VISBECK_VARIABLE_K
      &      + op5*(VisbeckK(i-1,j,bi,bj)+VisbeckK(i,j,bi,bj))
+     &           *GM_isoFac_calcK
 #endif
 #ifdef ALLOW_GM_LEITH_QG
      &      + op5*( GM_LeithQG_K(i-1,j,k,bi,bj)
      &            + GM_LeithQG_K(i,j,k,bi,bj) )
+     &           *GM_isoFac_calcK
 #endif
 #if ((defined GM_BATES_K3D) && ! (defined GM_BATES_PASSIVE))
      &      + op5*(GM_BatesK3d(i-1,j,k,bi,bj)+GM_BatesK3d(i,j,k,bi,bj))
+     &           *GM_isoFac_calcK
 #endif
      &      )*SlopeX(i,j)*taperFct(i,j)
          ENDDO
@@ -982,12 +993,15 @@ SUBROUTINE GMREDI_CALC_TENSOR(
 #endif
 #ifdef GM_VISBECK_VARIABLE_K
      &     + op5*(VisbeckK(i,j-1,bi,bj)+VisbeckK(i,j,bi,bj))
+     &          *GM_isoFac_calcK
 #endif
 #ifdef ALLOW_GM_LEITH_QG
      &     + op5*(GM_LeithQG_K(i,j-1,k,bi,bj)+GM_LeithQG_K(i,j,k,bi,bj))
+     &          *GM_isoFac_calcK
 #endif
 #if ((defined GM_BATES_K3D) && ! (defined GM_BATES_PASSIVE))
      &     + op5*(GM_BatesK3d(i,j-1,k,bi,bj)+GM_BatesK3d(i,j,k,bi,bj))
+     &          *GM_isoFac_calcK
 #endif
      &     )*taperFct(i,j)
          ENDDO
@@ -1029,15 +1043,18 @@ SUBROUTINE GMREDI_CALC_TENSOR(
      &        *op5*(GM_bolFac2d(i,j-1,bi,bj)+GM_bolFac2d(i,j,bi,bj))
 #endif
 #ifdef GM_VISBECK_VARIABLE_K
-     &     + op5*(VisbeckK(i,j-1,bi,bj)+VisbeckK(i,j,bi,bj))*GM_advect
+     &     + op5*(VisbeckK(i,j-1,bi,bj)+VisbeckK(i,j,bi,bj))
+     &          *(GM_isoFac_calcK - GM_skewflx)
 #endif
 #ifdef ALLOW_GM_LEITH_QG
      &     + op5*( GM_LeithQG_K(i,j-1,k,bi,bj)
-     &           + GM_LeithQG_K(i,j,k,bi,bj) )*GM_advect
+     &           + GM_LeithQG_K(i,j,k,bi,bj) )
+     &          *(GM_isoFac_calcK - GM_skewflx)
 #endif
 #if ((defined GM_BATES_K3D) && ! (defined GM_BATES_PASSIVE))
      &     + op5*( GM_BatesK3d(i,j-1,k,bi,bj)
-     &           + GM_BatesK3d(i,j,k,bi,bj) )*GM_advect
+     &           + GM_BatesK3d(i,j,k,bi,bj) )
+     &          *(GM_isoFac_calcK - GM_skewflx)
 #endif
      &     )*SlopeY(i,j)*taperFct(i,j)
          ENDDO
@@ -1067,13 +1084,16 @@ SUBROUTINE GMREDI_CALC_TENSOR(
 #endif
 #ifdef GM_VISBECK_VARIABLE_K
      &      + op5*(VisbeckK(i,j-1,bi,bj)+VisbeckK(i,j,bi,bj))
+     &           *GM_isoFac_calcK
 #endif
 #ifdef ALLOW_GM_LEITH_QG
      &      + op5*( GM_LeithQG_K(i,j-1,k,bi,bj)
      &            + GM_LeithQG_K(i,j,k,bi,bj) )
+     &           *GM_isoFac_calcK
 #endif
 #if ((defined GM_BATES_K3D) && ! (defined GM_BATES_PASSIVE))
      &      + op5*(GM_BatesK3d(i,j-1,k,bi,bj)+GM_BatesK3d(i,j,k,bi,bj))
+     &           *GM_isoFac_calcK
 #endif
      &      )*SlopeY(i,j)*taperFct(i,j)
          ENDDO
@@ -1109,6 +1129,8 @@ SUBROUTINE GMREDI_CALC_TENSOR(
 #endif /* GM_NON_UNITY_DIAGONAL || GM_EXTRA_DIAGONAL */
 
 #ifndef GM_NON_UNITY_DIAGONAL
+C--   keep a simplified setting here (used to be inside gmredi_x/ytransport)
+C     for backeard compatibility + for trying to generate a simpler adjoint code
       DO k=1,Nr
        DO j=1-OLy+1,sNy+OLy-1
         DO i=1-OLx+1,sNx+OLx-1
@@ -1118,9 +1140,6 @@ SUBROUTINE GMREDI_CALC_TENSOR(
      &           + GM_inpK3dRedi(i,j,k,bi,bj) )
 #else
      &       GM_isopycK
-#endif
-#ifdef GM_VISBECK_VARIABLE_K
-     &     + op5*(VisbeckK(i-1,j,bi,bj)+VisbeckK(i,j,bi,bj))
 #endif
      &                       )
         ENDDO
@@ -1133,9 +1152,6 @@ SUBROUTINE GMREDI_CALC_TENSOR(
      &           + GM_inpK3dRedi(i,j,k,bi,bj) )
 #else
      &       GM_isopycK
-#endif
-#ifdef GM_VISBECK_VARIABLE_K
-     &     + op5*(VisbeckK(i,j-1,bi,bj)+VisbeckK(i,j,bi,bj))
 #endif
      &                       )
         ENDDO

pkg/gmredi/gmredi_check.F

@@ -44,6 +44,7 @@ SUBROUTINE GMREDI_CHECK( myThid )
 C     msgBuf :: Informational/error message buffer
       CHARACTER*(MAX_LEN_MBUF) msgBuf
       INTEGER iL, errCount
+      _RL     tmpVar
 CEOP
 #ifdef ALLOW_PTRACERS
       INTEGER iTr
@@ -57,15 +58,7 @@ SUBROUTINE GMREDI_CHECK( myThid )
        CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
      &                     SQUEEZE_RIGHT , myThid )
 
-       IF ( GM_useLeithQG .AND. viscC2LeithQG.EQ.zeroRL ) THEN
-        WRITE(msgBuf,'(2A,I3)') '** WARNING ** GMREDI_CHECK: ',
-     &   'switch OFF GM_useLeithQG since viscC2LeithQG = 0'
-        CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
-     &                      SQUEEZE_RIGHT, myThid )
-        GM_useLeithQG = .FALSE.
-       ENDIF
-
-C- print out some kee parameters :
+C- print out some key parameters :
        CALL WRITE_0D_L( GM_AdvForm, INDEX_NONE,
      &  'GM_AdvForm =', '     /* if FALSE => use SkewFlux Form */')
        CALL WRITE_0D_L( GM_InMomAsStress, INDEX_NONE,
@@ -76,12 +69,18 @@ SUBROUTINE GMREDI_CHECK( myThid )
      &  'GM_ExtraDiag =','   /* Tensor Extra Diag (line 1&2) non 0 */')
        CALL WRITE_0D_RL( GM_isopycK, INDEX_NONE, 'GM_isopycK =',
      &  '    /* Background Isopyc. Diffusivity [m^2/s] */')
+c      CALL WRITE_0D_RL( GM_background_K, INDEX_NONE,
+c    &  'GM_background_K =',
+c    &  ' /* Background GM (=Bolus) Diffusivity [m^2/s] */')
+       tmpVar = GM_background_K*( oneRL - GM_skewflx )
+       CALL WRITE_0D_RL( tmpVar, INDEX_NONE, 'GM_advec*K =',
+     &  '    /* Backg. GM-Advec(=Bolus) Diffusivity [m^2/s] */')
        CALL WRITE_0D_RL( GM_background_K*GM_skewflx, INDEX_NONE,
      &  'GM_skewflx*K =',
      &  '  /* Background GM_SkewFlx Diffusivity [m^2/s] */')
-       CALL WRITE_0D_RL( GM_background_K*GM_advect, INDEX_NONE,
-     &  'GM_advec*K =',
-     &  '    /* Backg. GM-Advec(=Bolus) Diffusivity [m^2/s]*/')
+       CALL WRITE_0D_RL( GM_isoFac_calcK, INDEX_NONE,
+     &  'GM_isoFac_calcK =',
+     &  ' /* Fraction of dynamic K added to Redi tensor */')
        CALL WRITE_0D_RL( GM_Kmin_horiz, INDEX_NONE, 'GM_Kmin_horiz =',
      &  ' /* Minimum Horizontal Diffusivity [m^2/s] */')
        CALL WRITE_0D_RL( GM_Visbeck_alpha, INDEX_NONE,
@@ -276,15 +275,30 @@ SUBROUTINE GMREDI_CHECK( myThid )
        CALL PRINT_ERROR( msgBuf, myThid )
        errCount = errCount + 1
       ENDIF
-      IF ( GM_useLeithQG ) THEN
-       WRITE(msgBuf,'(A)')
-     &   'GMREDI_CHECK: GM_useLeithQG is only implemented for cases'
-       CALL PRINT_ERROR( msgBuf, myThid )
+      IF ( GM_isoFac_calcK .NE. zeroRL .AND.
+     &     ( GM_Visbeck_alpha.NE.zeroRL .OR.
+     &       GM_useBatesK3d .OR. GM_useLeithQG ) ) THEN
        WRITE(msgBuf,'(A)')
-     &   'GMREDI_CHECK: with #define GM_NON_UNITY_DIAGONAL'
+     &   'GMREDI_CHECK: needs #define GM_NON_UNITY_DIAGONAL'
        CALL PRINT_ERROR( msgBuf, myThid )
+       IF ( GM_Visbeck_alpha.NE.zeroRL ) THEN
+        WRITE(msgBuf,'(A)')
+     &   'GMREDI_CHECK: to use Visbeck computed K'
+        CALL PRINT_ERROR( msgBuf, myThid )
+       ENDIF
+       IF ( GM_useBatesK3d .OR. GM_useLeithQG ) THEN
+        WRITE(msgBuf,'(A)')
+     &   'GMREDI_CHECK: to use GM_useBatesK3d or GM_useLeithQG'
+        CALL PRINT_ERROR( msgBuf, myThid )
+       ENDIF
        errCount = errCount + 1
       ENDIF
+      IF ( errCount .EQ. 0 ) THEN
+        WRITE(msgBuf,'(2A)') '** WARNING ** GMREDI_CHECK: ',
+     &       '#undef GM_NON_UNITY_DIAGONAL not recommended'
+        CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
+     &                      SQUEEZE_RIGHT, myThid )
+      ENDIF
 #endif
 
 #ifndef ALLOW_GM_LEITH_QG