Index: /NEMO/branches/UKMO/NEMO_4.0.1_biharmonic_GM/cfgs/SHARED/namelist_ref
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.1_biharmonic_GM/cfgs/SHARED/namelist_ref	(revision 12534)
+++ /NEMO/branches/UKMO/NEMO_4.0.1_biharmonic_GM/cfgs/SHARED/namelist_ref	(revision 12535)
@@ -924,4 +924,5 @@
    ln_dynldf_lap = .false.     !    laplacian operator
    ln_dynldf_blp = .false.     !  bilaplacian operator
+   ln_dynldf_bgm = .false.     !  bilaplician Gent-McWilliams 
    !                       !  Direction of action  :
    ln_dynldf_lev = .false.     !  iso-level
@@ -948,4 +949,10 @@
       !                       !  iso-neutral laplacian operator (ln_dynldf_iso=T) :
       rn_ahm_b      = 0.0         !  background eddy viscosity  [m2/s]
+   nn_bhm_ijk_t  = 0           !  space/time variation of bilaplacian GM coefficient :
+      !                             !  = 11  uniform coefficient in interior, zero top and bottom
+      !                             !  = 12  linear tapering with dz to surface, zero at top and bottom
+      !                             !  = 13  dx^2.dz^2 scaling, zero at top and bottom
+      rn_bhm_b      = 0.0         ! bilaplacian GM coefficient
+      rn_bgmzcrit   = 0.0         ! critical depth for depth variation of bilaplacian GM coefficient
 /
 !-----------------------------------------------------------------------
Index: /NEMO/branches/UKMO/NEMO_4.0.1_biharmonic_GM/src/OCE/DIA/diaar5.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.1_biharmonic_GM/src/OCE/DIA/diaar5.F90	(revision 12534)
+++ /NEMO/branches/UKMO/NEMO_4.0.1_biharmonic_GM/src/OCE/DIA/diaar5.F90	(revision 12535)
@@ -42,5 +42,5 @@
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
-   !! $Id$
+   !! $Id: diaar5.F90 10425 2018-12-19 21:54:16Z smasson $
    !! Software governed by the CeCILL license (see ./LICENSE)
    !!----------------------------------------------------------------------
@@ -136,5 +136,6 @@
          !                                         ! steric sea surface height
          CALL eos( tsn, zrhd, zrhop, gdept_n(:,:,:) )                 ! now in situ and potential density
-         zrhop(:,:,jpk) = 0._wp
+         !CW: commenting out the following line, as seems weird to zero the array before output & only way to get rhop output
+         !zrhop(:,:,jpk) = 0._wp
          CALL iom_put( 'rhop', zrhop )
          !
Index: /NEMO/branches/UKMO/NEMO_4.0.1_biharmonic_GM/src/OCE/DYN/dynldf.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.1_biharmonic_GM/src/OCE/DYN/dynldf.F90	(revision 12534)
+++ /NEMO/branches/UKMO/NEMO_4.0.1_biharmonic_GM/src/OCE/DYN/dynldf.F90	(revision 12535)
@@ -7,4 +7,5 @@
    !!            3.7  ! 2014-01  (F. Lemarie, G. Madec)  restructuration/simplification of ahm specification,
    !!                 !                                  add velocity dependent coefficient and optional read in file
+   !!            4.0  ! 2020-02  (C. Wilson, ...) add bi-Laplacian GM implementation via momentum     
    !!----------------------------------------------------------------------
 
@@ -38,5 +39,5 @@
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
-   !! $Id$
+   !! $Id: dynldf.F90 10068 2018-08-28 14:09:04Z nicolasmartin $
    !! Software governed by the CeCILL license (see ./LICENSE)
    !!----------------------------------------------------------------------
@@ -66,4 +67,5 @@
       CASE ( np_lap   )    ;   CALL dyn_ldf_lap( kt, ub, vb, ua, va, 1 )      ! iso-level    laplacian
       CASE ( np_lap_i )    ;   CALL dyn_ldf_iso( kt )                         ! rotated      laplacian
+      CASE ( np_bgm )  ;       CALL dyn_ldf_bgm( kt, ub, vb, ua, va)  ! bi-Laplacian GM parameterisation
       CASE ( np_blp   )    ;   CALL dyn_ldf_blp( kt, ub, vb, ua, va    )      ! iso-level bi-laplacian
       !
@@ -104,4 +106,5 @@
          CASE( np_lap    )   ;   WRITE(numout,*) '   ==>>>   iso-level laplacian operator'
          CASE( np_lap_i  )   ;   WRITE(numout,*) '   ==>>>   rotated laplacian operator with iso-level background'
+         CASE( np_bgm )      ;   WRITE(numout,*) '   ==>>>   bi-Laplacian GM parameterisation via momentum equation'
          CASE( np_blp    )   ;   WRITE(numout,*) '   ==>>>   iso-level bi-laplacian operator'
          END SELECT
Index: /NEMO/branches/UKMO/NEMO_4.0.1_biharmonic_GM/src/OCE/DYN/dynldf_lap_blp.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.1_biharmonic_GM/src/OCE/DYN/dynldf_lap_blp.F90	(revision 12534)
+++ /NEMO/branches/UKMO/NEMO_4.0.1_biharmonic_GM/src/OCE/DYN/dynldf_lap_blp.F90	(revision 12535)
@@ -5,4 +5,5 @@
    !!======================================================================
    !! History : 3.7  ! 2014-01  (G. Madec, S. Masson)  Original code, re-entrant laplacian
+   !!           4.0  ! 2020-02  (C. Wilson, ...) add bhm coefficient for bi-Laplacian GM implementation via momentum     
    !!----------------------------------------------------------------------
 
@@ -25,4 +26,6 @@
    PUBLIC dyn_ldf_lap  ! called by dynldf.F90
    PUBLIC dyn_ldf_blp  ! called by dynldf.F90
+   PUBLIC dyn_ldf_bgm  ! called by dynldf.F90
+   PRIVATE dyn_ldf_lap_no_ahm !needed by dyn_ldf_bgm
 
    !! * Substitutions
@@ -30,5 +33,5 @@
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
-   !! $Id$
+   !! $Id: dynldf_lap_blp.F90 10425 2018-12-19 21:54:16Z smasson $ 
    !! Software governed by the CeCILL license (see ./LICENSE)
    !!----------------------------------------------------------------------
@@ -104,4 +107,168 @@
    END SUBROUTINE dyn_ldf_lap
 
+!CW: new subroutine for the Laplacian of velocity, copied from dyn_ldf_lap above, but without eddy viscosity ahmf, ahmt
+
+   SUBROUTINE dyn_ldf_lap_no_ahm( kt, pub, pvb, pulap, pvlap, kpass )
+      !!----------------------------------------------------------------------
+      !!                     ***  ROUTINE dyn_ldf_lap_no_ahm  ***
+      !!                       
+      !! ** Purpose :   Companion subroutine to dyn_ldf_bgm to calculate 
+      !!      the before horizontal momentum Laplacian 
+      !!      trend and add it to the general trend of momentum equation.  
+      !!      Note: mimics dyn_ldf_lap but without eddy viscosity ahmf, ahmt, 
+      !!      because biharmonic GM eddy diffusivity is applied in dyn_bgm.
+      !!
+      !! ** Method  :   The Laplacian operator apply on horizontal velocity is 
+      !!      writen as :   grad_h( ahmt div_h(U )) - curl_h( ahmf curl_z(U) ) 
+      !!
+      !! ** Action : - pua, pva increased by the harmonic operator applied on pub, pvb.
+      !!----------------------------------------------------------------------
+      INTEGER                         , INTENT(in   ) ::   kt         ! ocean time-step index
+      INTEGER                         , INTENT(in   ) ::   kpass      ! =1/2 first or second passage
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pub, pvb   ! before velocity  [m/s]
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(out) ::   pulap, pvlap ! velocity trend   [m/s2]
+      !
+      INTEGER  ::   ji, jj, jk   ! dummy loop indices
+      REAL(wp) ::   zsign        ! local scalars
+      REAL(wp) ::   zua, zva     ! local scalars
+      REAL(wp), DIMENSION(jpi,jpj) ::   zcur, zdiv
+      !!----------------------------------------------------------------------
+      !
+      IF( kt == nit000 .AND. lwp ) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) 'dyn_ldf_lap_no_ahm : companion iso-level harmonic (laplacian) operator to dyn_bgm, pass=', kpass
+         WRITE(numout,*) '~~~~~~~ '
+      ENDIF
+      !
+      IF( kpass == 1 ) THEN   ;   zsign =  1._wp      ! bilaplacian operator require a minus sign
+      ELSE                    ;   zsign = -1._wp      !  (eddy viscosity coef. >0)
+      ENDIF
+      !
+      !                                                ! ===============
+      DO jk = 1, jpkm1                                 ! Horizontal slab
+         !                                             ! ===============
+         DO jj = 2, jpj
+            DO ji = fs_2, jpi   ! vector opt.
+               !                                      ! ahm * e3 * curl  (computed from 1 to jpim1/jpjm1)
+!!gm open question here : e3f  at before or now ?    probably now...
+
+!!CW: note that the removed ahmf had already been multiplied by fmask, so we multiply by fmask explicitly here,
+!!    with the i,j,k of fmask aligning with those of ahmf, as defined in ldfdyn.F90.
+               zcur(ji-1,jj-1) = fmask(ji-1,jj-1,jk) * e3f_n(ji-1,jj-1,jk) * r1_e1e2f(ji-1,jj-1)      &
+                  &     * (  e2v(ji  ,jj-1) * pvb(ji  ,jj-1,jk) - e2v(ji-1,jj-1) * pvb(ji-1,jj-1,jk)  &
+                  &        - e1u(ji-1,jj  ) * pub(ji-1,jj  ,jk) + e1u(ji-1,jj-1) * pub(ji-1,jj-1,jk)  )
+               !                                      ! ahm * div        (computed from 2 to jpi/jpj)
+!!CW: note that the removed ahmt had already been multiplied by tmask, so we multiply by tmask explicitly here,
+!!    with the i,j,k of tmask aligning with those of ahmt, as defined in ldfdyn.F90
+               zdiv(ji,jj)     = tmask(ji,jj,jk) * r1_e1e2t(ji,jj) / e3t_b(ji,jj,jk)                                         &
+                  &     * (  e2u(ji,jj)*e3u_b(ji,jj,jk) * pub(ji,jj,jk) - e2u(ji-1,jj)*e3u_b(ji-1,jj,jk) * pub(ji-1,jj,jk)  &
+                  &        + e1v(ji,jj)*e3v_b(ji,jj,jk) * pvb(ji,jj,jk) - e1v(ji,jj-1)*e3v_b(ji,jj-1,jk) * pvb(ji,jj-1,jk)  )
+            END DO  
+         END DO  
+         !
+!CW: multiply pulap, pvlap by umask, vmask
+         DO jj = 2, jpjm1                             ! - curl( curl) + grad( div )
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               pulap(ji,jj,jk) = umask(ji,jj,jk) * zsign * (                                            &
+                  &              - ( zcur(ji  ,jj) - zcur(ji,jj-1) ) * r1_e2u(ji,jj) / e3u_n(ji,jj,jk)   &
+                  &              + ( zdiv(ji+1,jj) - zdiv(ji,jj  ) ) * r1_e1u(ji,jj)                     )
+                  !
+               pvlap(ji,jj,jk) = vmask(ji,jj,jk) * zsign * (                                             &
+                  &                ( zcur(ji,jj  ) - zcur(ji-1,jj) ) * r1_e1v(ji,jj) / e3v_n(ji,jj,jk)   &
+                  &              + ( zdiv(ji,jj+1) - zdiv(ji  ,jj) ) * r1_e2v(ji,jj)                     )
+            END DO
+         END DO
+         !                                             ! ===============
+      END DO                                           !   End of slab
+      !                                                ! ===============
+      !
+   END SUBROUTINE dyn_ldf_lap_no_ahm
+
+
+
+
+!CW: new subroutine for biharmonic GM, copied from SUBROUTINE dyn_ldf_blp below
+
+   SUBROUTINE dyn_ldf_bgm( kt, pub, pvb, pua, pva )
+      !!----------------------------------------------------------------------
+      !!                 ***  ROUTINE dyn_bgm  ***
+      !!                    
+      !! ** Purpose :   Compute the before lateral momentum trend due to the bi-Laplacian GM parameterisation 
+      !!              and add it to the general trend of momentum equation.
+      !!
+      !! ** Method  :   The bi-Laplacian implementation of GM is via a -d/dz(diffusivity x d/dz(Laplacian of velocity))
+      !!      operator applied at the 'now' time level.  The existing code for the Laplacian contains the 'before' time also in zdiv.
+      !!      It is computed by a call to dyn_ldf_lap routine and vertical differentiation applied twice.
+      !!
+      !! ** Action :   pua, pva increased with the before bi-Laplacian GM momentum trend calculated from pub, pvb.
+      !!----------------------------------------------------------------------
+      INTEGER                         , INTENT(in   ) ::   kt         ! ocean time-step index
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pub, pvb   ! before velocity fields
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pua, pva   ! momentum trend
+      !   
+      INTEGER  ::   iku, ikv                         ! local integers
+      !CW
+      INTEGER  ::   ji, jj, jk   ! dummy loop indices
+      !
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zulap, zvlap   ! laplacian at u- and v-point
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zulapdz, zvlapdz   ! -1*bhm * d/dz(del^2 u) at u- and v-point
+      !!----------------------------------------------------------------------
+      !
+      IF( kt == nit000 )  THEN
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) 'dyn_bgm : bi-Laplacian GM operator via momentum '
+         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
+      ENDIF
+      !
+  
+      CALL dyn_ldf_lap_no_ahm( kt, pub, pvb, zulap, zvlap, 1 )   
+
+      ! ===============
+      !CW: calculate -bhm * d/dz(del^2 u)  
+      DO jk = 2, jpkm1 
+         DO jj = 2, jpjm1
+            DO ji = fs_2, jpim1   ! vector opt.
+                zulapdz(ji,jj,jk) = -bhmu(ji,jj,jk)*(zulap(ji,jj,jk-1) - zulap(ji,jj,jk)) / e3uw_n(ji,jj,jk)
+                zvlapdz(ji,jj,jk) = -bhmv(ji,jj,jk)*(zvlap(ji,jj,jk-1) - zvlap(ji,jj,jk)) / e3vw_n(ji,jj,jk)
+            ENDDO
+         ENDDO
+      ENDDO
+
+
+!CW: set boundary conditions: d/dz(del^2 u) = 0 at top and bottom, so that eddy-induced velocity, w*=0
+!     Surface
+      zulapdz(:,:,1)   = 0._wp   ;   zvlapdz(:,:,1)   = 0._wp 
+!     Flat bottom case
+      zulapdz(:,:,jpk) = 0._wp   ;   zvlapdz(:,:,jpk) = 0._wp
+!     Variable bathymetry case, including z-partial steps, as in dynhpg.F90, subroutine hpg_zps
+        DO jj = 2, jpjm1
+            DO ji = 2, jpim1
+                iku = mbku(ji,jj)+1
+                ikv = mbkv(ji,jj)+1
+                zulapdz(:,:,iku) = 0._wp
+                zvlapdz(:,:,ikv) = 0._wp
+            ENDDO
+        ENDDO
+
+
+
+
+!!    calculate d/dz(-bhm * d/dz(del^2 u))
+                                                       ! ===============
+      DO jk = 1, jpkm1                                 ! Horizontal slab
+         !                                             ! ===============
+         DO jj = 2, jpjm1
+            DO ji = fs_2, jpim1   ! vector opt.
+            pua(ji,jj,jk) = pua(ji,jj,jk) + (zulapdz(ji,jj,jk) - zulapdz(ji,jj,jk+1)) / e3u_n(ji,jj,jk)
+            pva(ji,jj,jk) = pva(ji,jj,jk) + (zvlapdz(ji,jj,jk) - zvlapdz(ji,jj,jk+1)) / e3v_n(ji,jj,jk)
+            ENDDO
+         ENDDO
+      ENDDO
+
+      ! -----
+
+   
+   END SUBROUTINE dyn_ldf_bgm
+   
 
    SUBROUTINE dyn_ldf_blp( kt, pub, pvb, pua, pva )
Index: /NEMO/branches/UKMO/NEMO_4.0.1_biharmonic_GM/src/OCE/LDF/ldfdyn.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.1_biharmonic_GM/src/OCE/LDF/ldfdyn.F90	(revision 12534)
+++ /NEMO/branches/UKMO/NEMO_4.0.1_biharmonic_GM/src/OCE/LDF/ldfdyn.F90	(revision 12535)
@@ -8,4 +8,5 @@
    !!            3.7  ! 2014-01  (F. Lemarie, G. Madec)  restructuration/simplification of ahm specification,
    !!                 !                                  add velocity dependent coefficient and optional read in file
+   !!            4.0  ! 2020-02  (C. Wilson, ...) add bhm coefficient for bi-Laplacian GM implementation via momentum     
    !!----------------------------------------------------------------------
 
@@ -15,5 +16,5 @@
    !!----------------------------------------------------------------------
    USE oce             ! ocean dynamics and tracers   
-   USE dom_oce         ! ocean space and time domain 
+   USE dom_oce         ! ocean space and time domain
    USE phycst          ! physical constants
    USE ldfslp          ! lateral diffusion: slopes of mixing orientation
@@ -35,4 +36,5 @@
    LOGICAL , PUBLIC ::   ln_dynldf_OFF   !: No operator (i.e. no explicit diffusion)
    LOGICAL , PUBLIC ::   ln_dynldf_lap   !: laplacian operator
+   LOGICAL , PUBLIC ::   ln_dynldf_bgm   !: bi-Laplacian GM implementation via momentum
    LOGICAL , PUBLIC ::   ln_dynldf_blp   !: bilaplacian operator
    LOGICAL , PUBLIC ::   ln_dynldf_lev   !: iso-level direction
@@ -42,4 +44,5 @@
    !                                        !  time invariant coefficients:  aht = 1/2  Ud*Ld   (lap case) 
    !                                           !                             bht = 1/12 Ud*Ld^3 (blp case)
+   INTEGER , PUBLIC ::   nn_bhm_ijk_t    !: choice of time & space variations of the lateral bi-Laplacian GM coef.
    REAL(wp), PUBLIC ::   rn_Uv                 !: lateral viscous velocity  [m/s]
    REAL(wp), PUBLIC ::   rn_Lv                 !: lateral viscous length    [m]
@@ -50,5 +53,6 @@
    !                                        ! iso-neutral laplacian (ln_dynldf_lap=ln_dynldf_iso=T)
    REAL(wp), PUBLIC ::   rn_ahm_b              !: lateral laplacian background eddy viscosity  [m2/s]
-
+   REAL(wp), PUBLIC ::   rn_bhm_b              !: lateral bi-Laplacian GM background eddy viscosity  [m4/s]
+   REAL(wp), PUBLIC ::   rn_bgmzcrit           !: critical depth (>=0) for linear tapering of bi-Lap. GM coef. to zero at surface [m] 
    !                                    !!* Parameter to control the type of lateral viscous operator
    INTEGER, PARAMETER, PUBLIC ::   np_ERROR  =-10                       !: error in setting the operator
@@ -57,9 +61,10 @@
    INTEGER, PARAMETER, PUBLIC ::   np_lap    = 10   ,   np_blp    = 20  !: iso-level operator
    INTEGER, PARAMETER, PUBLIC ::   np_lap_i  = 11                       !: iso-neutral or geopotential operator
-   !
+   INTEGER, PARAMETER, PUBLIC ::   np_bgm    = 30                       !: iso-level operator, bi-Laplacian GM
    INTEGER           , PUBLIC ::   nldf_dyn         !: type of lateral diffusion used defined from ln_dynldf_... (namlist logicals)
    LOGICAL           , PUBLIC ::   l_ldfdyn_time    !: flag for time variation of the lateral eddy viscosity coef.
 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ahmt, ahmf   !: eddy viscosity coef. at T- and F-points [m2/s or m4/s]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   bhmu, bhmv   !: eddy viscosity coef. for bi-Laplacian GM at u- and v-points [m4/s]
    REAL(wp),         ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   dtensq       !: horizontal tension squared         (Smagorinsky only)
    REAL(wp),         ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   dshesq       !: horizontal shearing strain squared (Smagorinsky only)
@@ -76,5 +81,5 @@
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
-   !! $Id$
+   !! $Id: ldfdyn.F90 11653 2019-10-04 12:34:02Z davestorkey $ 
    !! Software governed by the CeCILL license (see ./LICENSE)
    !!----------------------------------------------------------------------
@@ -91,4 +96,5 @@
       !!             ln_dynldf_lap = T     laplacian operator
       !!             ln_dynldf_blp = T   bilaplacian operator
+      !!             ln_dynldf_bgm = T   bi-Laplacian GM operator
       !!              - the parameter nn_ahm_ijk_t:
       !!    nn_ahm_ijk_t  =  0 => = constant
@@ -103,13 +109,22 @@
       !!                                                           (   L^2|D|      laplacian operator
       !!                                                           or  L^4|D|/8  bilaplacian operator )
+      !!                - the parameter nn_bhm_ijk_t:
+      !!    nn_bhm_ijk_t  = 11 => = F(z) :     = constant, with BCs set to zero at top and bottom
+      !!    nn_bhm_ijk_t  = 12 => = F(z) :     = constant in interior, linear taper to zero at surface,
+      !!                                         for depths < rn_bgmzcrit, with BCs set to zero at top and bottom
+      !!    nn_bhm_ijk_t  = 13 => = F(z) :     = bhm0 X dx^2 X dz^2, with BCs set to zero at top and bottom
       !!----------------------------------------------------------------------
       INTEGER  ::   ji, jj, jk                     ! dummy loop indices
       INTEGER  ::   ioptio, ierr, inum, ios, inn   ! local integer
+      INTEGER  ::   iku, ikv, ikt                  ! local integer
       REAL(wp) ::   zah0, zah_max, zUfac           ! local scalar
+      REAL(wp) ::   bhm0, lhscrit, lhscritmax      ! local scalar
       CHARACTER(len=5) ::   cl_Units               ! units (m2/s or m4/s)
       !!
       NAMELIST/namdyn_ldf/ ln_dynldf_OFF, ln_dynldf_lap, ln_dynldf_blp,   &   ! type of operator
+         &                 ln_dynldf_bgm,                                 &   ! type of operator
          &                 ln_dynldf_lev, ln_dynldf_hor, ln_dynldf_iso,   &   ! acting direction of the operator
          &                 nn_ahm_ijk_t , rn_Uv    , rn_Lv,   rn_ahm_b,   &   ! lateral eddy coefficient
+         &                 nn_bhm_ijk_t, rn_bhm_b, rn_bgmzcrit,           &   ! lateral bi-Lap. GM coefficient
          &                 rn_csmc      , rn_minfac    , rn_maxfac            ! Smagorinsky settings
       !!----------------------------------------------------------------------
@@ -134,4 +149,6 @@
          WRITE(numout,*) '         laplacian operator                   ln_dynldf_lap = ', ln_dynldf_lap
          WRITE(numout,*) '         bilaplacian operator                 ln_dynldf_blp = ', ln_dynldf_blp
+         WRITE(numout,*) '         bi-Laplacian GM operator             ln_dynldf_bgm = ', ln_dynldf_bgm
+         WRITE(numout,*) '         critical depth for taper (if used)   rn_bgmzcrit   = ', rn_bgmzcrit
          !
          WRITE(numout,*) '      direction of action :'
@@ -145,5 +162,5 @@
          WRITE(numout,*) '         lateral viscous length    (if cst)      rn_Lv      = ', rn_Lv, ' m'
          WRITE(numout,*) '         background viscosity (iso-lap case)     rn_ahm_b   = ', rn_ahm_b, ' m2/s'
-         !
+         WRITE(numout,*) '         bi-Laplacian GM background viscosity    rn_bhm_b   = ', rn_bhm_b, ' m4/s'
          WRITE(numout,*) '      Smagorinsky settings (nn_ahm_ijk_t  = 32) :'
          WRITE(numout,*) '         Smagorinsky coefficient              rn_csmc       = ', rn_csmc
@@ -159,9 +176,13 @@
       nldf_dyn = np_ERROR
       ioptio = 0
+      !CW exclude combination of lap+blp (as in existing code), but allow other combinations
       IF( ln_dynldf_OFF ) THEN   ;   nldf_dyn = np_no_ldf   ;   ioptio = ioptio + 1   ;   ENDIF
-      IF( ln_dynldf_lap ) THEN   ;                              ioptio = ioptio + 1   ;   ENDIF
-      IF( ln_dynldf_blp ) THEN   ;                              ioptio = ioptio + 1   ;   ENDIF
-      IF( ioptio /= 1   )   CALL ctl_stop( 'dyn_ldf_init: use ONE of the 3 operator options (NONE/lap/blp)' )
-      !
+      IF( ln_dynldf_lap ) THEN   ;                              ioptio = ioptio + 5   ;   ENDIF
+      IF( ln_dynldf_blp ) THEN   ;                              ioptio = ioptio + 5   ;   ENDIF
+      IF( ln_dynldf_bgm ) THEN   ;                              ioptio = ioptio + 1   ;   ENDIF
+      !
+      IF( (ioptio < 1).OR.(ioptio > 6) )   CALL ctl_stop( 'dyn_ldf_init: use at least ONE of the 4 operator options (NONE/lap/blp/bgm) but not (lap+blp)' )
+      !
+
       IF(.NOT.ln_dynldf_OFF ) THEN     !==  direction ==>> type of operator  ==!
          ioptio = 0
@@ -209,6 +230,28 @@
          ENDIF
          !
+         IF( ln_dynldf_bgm ) THEN         ! bi-Laplacian GM operator
+            IF( ln_zco ) THEN                   ! z-coordinate
+               IF( ln_dynldf_lev )   nldf_dyn = np_bgm   ! iso-level = horizontal (no rotation)
+               IF( ln_dynldf_hor )   nldf_dyn = np_bgm   ! iso-level = horizontal (no rotation)
+               IF( ln_dynldf_iso )   ierr = 3            ! iso-neutral            (   rotation)
+            ENDIF
+            IF( ln_zps ) THEN                   ! z-coordinate with partial step
+               IF( ln_dynldf_lev )   nldf_dyn = np_bgm   ! iso-level              (no rotation)
+               IF( ln_dynldf_hor )   nldf_dyn = np_bgm   ! iso-level              (no rotation)
+               IF( ln_dynldf_iso )   ierr = 3            ! iso-neutral            (   rotation)
+            ENDIF
+            IF( ln_sco ) THEN                   ! s-coordinate
+               IF( ln_dynldf_lev )   nldf_dyn = np_bgm   ! iso-level              (no rotation)
+               IF( ln_dynldf_hor )   ierr = 3            ! horizontal             (   rotation)
+               IF( ln_dynldf_iso )   ierr = 3            ! iso-neutral            (   rotation)
+            ENDIF
+         ENDIF
+         !         
+
          IF( ierr == 2 )   CALL ctl_stop( 'rotated bi-laplacian operator does not exist' )
          !
+         IF( ierr == 3 )   CALL ctl_stop( 'rotated bi-Laplacian GM operator does not exist' )
+         !
+
          IF( nldf_dyn == np_lap_i )   l_ldfslp = .TRUE.  ! rotation require the computation of the slopes
          !
@@ -222,4 +265,5 @@
          CASE( np_lap_i  )   ;   WRITE(numout,*) '   ==>>>   rotated laplacian operator with iso-level background'
          CASE( np_blp    )   ;   WRITE(numout,*) '   ==>>>   iso-level bi-laplacian operator'
+         CASE( np_bgm    )   ;   WRITE(numout,*) '   ==>>>   iso-level bi-Laplacian GM operator'
          END SELECT
          WRITE(numout,*)
@@ -233,5 +277,8 @@
       !
       IF( ln_dynldf_OFF ) THEN
-         IF(lwp) WRITE(numout,*) '   ==>>>   No viscous operator selected. ahmt and ahmf are not allocated'
+         !CW
+         !IF(lwp) WRITE(numout,*) '   ==>>>   No viscous operator selected. ahmt and ahmf are not allocated'
+         IF(lwp) WRITE(numout,*) '   ==>>>   No viscous operator selected. ahmt, ahmf, bhmu, bhmv are not allocated'
+         !
          RETURN
          !
@@ -239,13 +286,16 @@
          !
          !                                         ! allocate the ahm arrays
-         ALLOCATE( ahmt(jpi,jpj,jpk) , ahmf(jpi,jpj,jpk) , STAT=ierr )
+         ALLOCATE( ahmt(jpi,jpj,jpk) , ahmf(jpi,jpj,jpk) , bhmu(jpi,jpj,jpk) , bhmv(jpi,jpj,jpk) , STAT=ierr )
+         !
          IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'ldf_dyn_init: failed to allocate arrays')
          !
          ahmt(:,:,:) = 0._wp                       ! init to 0 needed 
          ahmf(:,:,:) = 0._wp
-         !
-         !                                         ! value of lap/blp eddy mixing coef.
+         bhmu(:,:,:) = 0._wp                       ! init to 0 needed 
+         bhmv(:,:,:) = 0._wp
+         !                                        ! value of lap/blp eddy mixing coef.
          IF(     ln_dynldf_lap ) THEN   ;   zUfac = r1_2 *rn_Uv   ;   inn = 1   ;   cl_Units = ' m2/s'   !   laplacian
          ELSEIF( ln_dynldf_blp ) THEN   ;   zUfac = r1_12*rn_Uv   ;   inn = 3   ;   cl_Units = ' m4/s'   ! bilaplacian
+         ELSEIF( ln_dynldf_bgm ) THEN   ;   bhm0 = rn_bhm_b       ;             ;   cl_Units = ' m4/s'   ! bi-Laplacian GM
          ENDIF
          zah0    = zUfac *    rn_Lv**inn              ! mixing coefficient
@@ -265,5 +315,5 @@
             ahmf(:,:,1) = zah0
             CALL ldf_c1d( 'DYN', ahmt(:,:,1), ahmf(:,:,1), ahmt, ahmf )
-            !
+            !       
          CASE ( -20 )      !== fixed horizontal shape read in file  ==!
             IF(lwp) WRITE(numout,*) '   ==>>>   eddy viscosity = F(i,j) read in eddy_viscosity.nc file'
@@ -323,4 +373,175 @@
             CALL ctl_stop('ldf_dyn_init: wrong choice for nn_ahm_ijk_t, the type of space-time variation of ahm')
          END SELECT
+!CW      Add separate structure options for bi-Laplacian GM, to allow it to be used in combination with other types of diffusion, above.
+
+         IF( ln_dynldf_bgm ) THEN
+            SELECT CASE (nn_bhm_ijk_t)
+               !CW Define bi-Laplacian GM diffusivity and test the related computational stability criterion   
+
+            CASE(  11  )      !==  fixed profile: constant except for zero at top and bottom, so that eddy-induced velocity, w*=0  ==!
+               IF(lwp) WRITE(numout,*) '   ==>>>   bi-Laplacian GM eddy viscosity = constant in interior, zero at top and bottom'
+               IF(lwp) WRITE(numout,*) '           interior viscous coef. = constant = ', bhm0, cl_Units
+
+               !     First set to constant in interior, all levels
+               DO jj = 2, jpjm1
+                  DO ji = 2, jpim1
+                     bhmu(ji,jj,:) = bhm0
+                     bhmv(ji,jj,:) = bhm0
+                  ENDDO
+               ENDDO
+
+               !     Surface BC : set to zero
+               bhmu(:,:,1)   = 0._wp   ;   bhmv(:,:,1)   = 0._wp
+               !     Flat bottom BC : set to zero
+               bhmu(:,:,jpk) = 0._wp   ;   bhmv(:,:,jpk) = 0._wp
+               !     Variable bathymetry case, including z-partial steps, as in dynhpg.F90, subroutine hpg_zps
+               !     and exactly mirroring top, bottom BC d/dz(del^2 u) = 0 for BGM calculation in dynldf_lap_blp
+               DO jj = 2, jpjm1
+                  DO ji = 2, jpim1
+                     iku = mbku(ji,jj)+1
+                     ikv = mbkv(ji,jj)+1
+                     bhmu(:,:,iku) = 0._wp
+                     bhmv(:,:,ikv) = 0._wp
+                  ENDDO
+               ENDDO
+
+               !CW Test criterion for two grid-length mode for BGM scheme with leapfrog timestepping
+               ! \frac{{\kappa \Delta t}{(\Delta z)^2 (\Delta x)^2}< 1/32 
+               ! test at T-point - may need to revise this choice!!!!
+
+               ! Find domain maximum value of LHS, then test inequality.
+
+               ! initialise
+               lhscritmax=0._wp
+
+               DO jj = 2, jpjm1
+                  DO ji = 2, jpim1
+                     ikt = mbkt(ji,jj) !bottom last wet T-level
+                     DO jk = 2, ikt
+                        lhscrit=bhm0*rdt/(e3t_n(ji,jj,jk)**2*e1t(ji,jj)**2)
+                        IF(lhscrit .gt. lhscritmax) lhscritmax=lhscrit
+                     ENDDO
+                  ENDDO
+               ENDDO
+
+               IF ( lhscrit .ge. 1._wp/32._wp) THEN
+                  IF(lwp) THEN
+                     WRITE(numout,*)
+                     WRITE(numout,*) ' WARNING : Bi-Laplacian GM eddy viscosity is not',              &
+                          &            ' consistent with the model time step and grid setup: ', &
+                          &            ' LHS=',lhscrit, & 
+                          &            'Should be < 0.03125 in order to avoid computational instability (Mike Bell, Dave Storkey)'
+                     WRITE(numout,*) ' ========='
+                     WRITE(numout,*)
+                  ENDIF
+                  ! CW: warn, don't stop, as we may wish to violate this theoretical limit.
+                  ! CALL ctl_stop( 'ldf_dyn_init', 'Inconsistent Bi-Lap. GM diffusivity')
+               ELSE
+                  IF(lwp) THEN
+                     WRITE(numout,*)
+                     WRITE(numout,*) ' INFORMATION : Bi-Laplacian GM eddy viscosity is ',              &
+                          &            ' consistent with the model time step and grid setup: ', &
+                          &            ' LHS=',lhscrit, & 
+                          &            'Should be < 0.03125 in order to avoid computational instability (Mike Bell, Dave Storkey)'
+                     WRITE(numout,*) ' ========='
+                     WRITE(numout,*)
+                  ENDIF
+               ENDIF
+               !-----------------------------
+            CASE(  12  )      !==  fixed profile: constant in interior, zero at bottom and linearly-tapering to zero at top, over depth shallower than rn_bgmzcrit  ==!
+               IF(lwp) WRITE(numout,*) '   ==>>>   bi-Laplacian GM eddy viscosity = constant in interior, zero at bottom and'
+               IF(lwp) WRITE(numout,*) '           linearly-tapering to zero at top, over depth shallower than ', rn_bgmzcrit, ' metres.'
+               IF(lwp) WRITE(numout,*) '           Interior viscous coef. = constant = ', bhm0, cl_Units
+
+               !     Impose linear taper from interior value to zero at zero depth, for depths < rn_bgmzcrit
+               !     N.B. Test criterion for two grid-length mode for BGM scheme with leapfrog timestepping not implemented for this case yet.
+
+               DO jk = 1, jpk
+
+                  IF (gdept_1d(jk) .lt. rn_bgmzcrit) THEN
+
+                     DO jj = 2, jpjm1
+                        DO ji = 2, jpim1
+                           bhmu(ji,jj,jk) = bhm0 * gdept_1d(jk)/rn_bgmzcrit
+                           bhmv(ji,jj,jk) = bhm0 * gdept_1d(jk)/rn_bgmzcrit
+                        ENDDO
+                     ENDDO
+
+                  ELSE
+
+                     ! constant at depth rn_bgmzcrit and larger
+                     DO jj = 2, jpjm1
+                        DO ji = 2, jpim1
+                           bhmu(ji,jj,jk) = bhm0
+                           bhmv(ji,jj,jk) = bhm0
+                        ENDDO
+                     ENDDO
+
+                  ENDIF
+
+               ENDDO
+
+
+               !     Surface BC : set to zero
+               bhmu(:,:,1)   = 0._wp   ;   bhmv(:,:,1)   = 0._wp
+               !     Flat bottom BC : set to zero
+               bhmu(:,:,jpk) = 0._wp   ;   bhmv(:,:,jpk) = 0._wp
+               !     Variable bathymetry case, including z-partial steps, as in dynhpg.F90, subroutine hpg_zps
+               !     and exactly mirroring top, bottom BC d/dz(del^2 u) = 0 for BGM calculation in dynldf_lap_blp
+               DO jj = 2, jpjm1
+                  DO ji = 2, jpim1
+                     iku = mbku(ji,jj)+1
+                     ikv = mbkv(ji,jj)+1
+                     bhmu(:,:,iku) = 0._wp
+                     bhmv(:,:,ikv) = 0._wp
+                  ENDDO
+               ENDDO
+
+               !--------------------------------
+            CASE(  13  )      !==  fixed profile: steady profile of the form bhm0 X dx^2 X dz^2 in interior, zero at bottom and top  ==!
+
+               cl_Units = ' 1/s'
+
+               IF(lwp) WRITE(numout,*) '   ==>>>   bi-Laplacian GM eddy viscosity : steady profile of the form'
+               IF(lwp) WRITE(numout,*) '           bhm0 X dx^2 X dz^2 in interior, zero boundary conds. at bottom and top.'
+               IF(lwp) WRITE(numout,*) '           In this case, bhm0 is the constant of proportionality,'
+               IF(lwp) WRITE(numout,*) '           dimensioned as [diffusivity]/L^4, i.e. bhm0=', bhm0, cl_Units
+
+               cl_Units = ' m4/s'
+
+               !     N.B. Test criterion for two grid-length mode for BGM scheme with leapfrog timestepping not implemented for this case yet.
+
+               DO jk = 1, jpk
+                  ! constant at depth rn_bgmzcrit and larger
+                  DO jj = 2, jpjm1
+                     DO ji = 2, jpim1
+                        bhmu(ji,jj,jk) = bhm0*e1u(ji,jj)**2*e3t_n(ji,jj,jk)**2
+                        bhmv(ji,jj,jk) = bhm0*e1v(ji,jj)**2*e3t_n(ji,jj,jk)**2
+                     ENDDO
+                  ENDDO
+               ENDDO
+
+
+               !     Surface BC : set to zero
+               bhmu(:,:,1)   = 0._wp   ;   bhmv(:,:,1)   = 0._wp
+               !     Flat bottom BC : set to zero
+               bhmu(:,:,jpk) = 0._wp   ;   bhmv(:,:,jpk) = 0._wp
+               !     Variable bathymetry case, including z-partial steps, as in dynhpg.F90, subroutine hpg_zps
+               !     and exactly mirroring top, bottom BC d/dz(del^2 u) = 0 for BGM calculation in dynldf_lap_blp
+               DO jj = 2, jpjm1
+                  DO ji = 2, jpim1
+                     iku = mbku(ji,jj)+1
+                     ikv = mbkv(ji,jj)+1
+                     bhmu(:,:,iku) = 0._wp
+                     bhmv(:,:,ikv) = 0._wp
+                  ENDDO
+               ENDDO
+
+               !--------------------------------
+
+            CASE DEFAULT
+               CALL ctl_stop('ldf_dyn_init: wrong choice for nn_bhm_ijk_t, the type of space-time variation of bhm')
+            END SELECT
+         ENDIF ! ln_dynldf_bgm
          !
          IF( .NOT.l_ldfdyn_time ) THEN             !* No time variation 
@@ -331,4 +552,9 @@
                ahmt(:,:,1:jpkm1) = SQRT( ahmt(:,:,1:jpkm1) ) * tmask(:,:,1:jpkm1)
                ahmf(:,:,1:jpkm1) = SQRT( ahmf(:,:,1:jpkm1) ) * fmask(:,:,1:jpkm1)
+               !CW
+            ELSEIF( ln_dynldf_bgm ) THEN                 !bi-Laplacian GM operator (mask only)
+               bhmu(:,:,1:jpkm1) =       bhmu(:,:,1:jpkm1)   * umask(:,:,1:jpkm1)
+               bhmv(:,:,1:jpkm1) =       bhmv(:,:,1:jpkm1)   * vmask(:,:,1:jpkm1)
+               !
             ENDIF
          ENDIF
@@ -336,8 +562,8 @@
       ENDIF
       !
-   END SUBROUTINE ldf_dyn_init
-
-
-   SUBROUTINE ldf_dyn( kt )
+    END SUBROUTINE ldf_dyn_init
+
+
+    SUBROUTINE ldf_dyn( kt )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE ldf_dyn  ***
@@ -366,5 +592,5 @@
       !
       SELECT CASE(  nn_ahm_ijk_t  )       !== Eddy vicosity coefficients ==!
-      !
+         !
       CASE(  31  )       !==  time varying 3D field  ==!   = F( local velocity )
          !
@@ -426,5 +652,5 @@
                   DO ji = 2, jpim1
                      zdb =   ( ub(ji,jj,jk) * r1_e2u(ji,jj) - ub(ji-1,jj,jk) * r1_e2u(ji-1,jj) ) * r1_e1t(ji,jj) * e2t(ji,jj)  &
-                        &  - ( vb(ji,jj,jk) * r1_e1v(ji,jj) - vb(ji,jj-1,jk) * r1_e1v(ji,jj-1) ) * r1_e2t(ji,jj) * e1t(ji,jj)
+                          &  - ( vb(ji,jj,jk) * r1_e1v(ji,jj) - vb(ji,jj-1,jk) * r1_e1v(ji,jj-1) ) * r1_e2t(ji,jj) * e1t(ji,jj)
                      dtensq(ji,jj,jk) = zdb * zdb * tmask(ji,jj,jk)
                   END DO
@@ -434,5 +660,5 @@
                   DO ji = 1, jpim1
                      zdb =   ( ub(ji,jj+1,jk) * r1_e1u(ji,jj+1) - ub(ji,jj,jk) * r1_e1u(ji,jj) ) * r1_e2f(ji,jj) * e1f(ji,jj)  &
-                        &  + ( vb(ji+1,jj,jk) * r1_e2v(ji+1,jj) - vb(ji,jj,jk) * r1_e2v(ji,jj) ) * r1_e1f(ji,jj) * e2f(ji,jj)
+                          &  + ( vb(ji+1,jj,jk) * r1_e2v(ji+1,jj) - vb(ji,jj,jk) * r1_e2v(ji,jj) ) * r1_e1f(ji,jj) * e2f(ji,jj)
                      dshesq(ji,jj,jk) = zdb * zdb * fmask(ji,jj,jk)
                   END DO
@@ -444,5 +670,5 @@
             !
             DO jk = 1, jpkm1
-              !
+               !
                DO jj = 2, jpjm1                                ! T-point value
                   DO ji = fs_2, fs_jpim1
@@ -453,6 +679,6 @@
                      zdelta         = zcmsmag * esqt(ji,jj)                                        ! L^2 * (C_smag/pi)^2
                      ahmt(ji,jj,jk) = zdelta * SQRT(          dtensq(ji  ,jj,jk) +                         &
-                        &                            r1_4 * ( dshesq(ji  ,jj,jk) + dshesq(ji  ,jj-1,jk) +  &
-                        &                                     dshesq(ji-1,jj,jk) + dshesq(ji-1,jj-1,jk) ) )
+                          &                            r1_4 * ( dshesq(ji  ,jj,jk) + dshesq(ji  ,jj-1,jk) +  &
+                          &                                     dshesq(ji-1,jj,jk) + dshesq(ji-1,jj-1,jk) ) )
                      ahmt(ji,jj,jk) = MAX( ahmt(ji,jj,jk), SQRT( (zu2pv2_ij + zu2pv2_ij_m1) * zdelta * zstabf_lo ) ) ! Impose lower limit == minfac  * |U|L/2
                      ahmt(ji,jj,jk) = MIN( ahmt(ji,jj,jk),                                    zdelta * zstabf_up )   ! Impose upper limit == maxfac  * L^2/(4*2dt)
@@ -469,6 +695,6 @@
                      zdelta         = zcmsmag * esqf(ji,jj)                                        ! L^2 * (C_smag/pi)^2
                      ahmf(ji,jj,jk) = zdelta * SQRT(          dshesq(ji  ,jj,jk) +                         &
-                        &                            r1_4 * ( dtensq(ji  ,jj,jk) + dtensq(ji  ,jj+1,jk) +  &
-                        &                                     dtensq(ji+1,jj,jk) + dtensq(ji+1,jj+1,jk) ) )
+                          &                            r1_4 * ( dtensq(ji  ,jj,jk) + dtensq(ji  ,jj+1,jk) +  &
+                          &                                     dtensq(ji+1,jj,jk) + dtensq(ji+1,jj+1,jk) ) )
                      ahmf(ji,jj,jk) = MAX( ahmf(ji,jj,jk), SQRT( (zu2pv2_ij + zu2pv2_ij_p1) * zdelta * zstabf_lo ) ) ! Impose lower limit == minfac  * |U|L/2
                      ahmf(ji,jj,jk) = MIN( ahmf(ji,jj,jk),                                    zdelta * zstabf_up )   ! Impose upper limit == maxfac  * L^2/(4*2dt)
