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 12775)
+++ NEMO/branches/UKMO/NEMO_4.0.1_biharmonic_GM/src/OCE/DYN/dynldf.F90	(revision 12776)
@@ -67,8 +67,9 @@
       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
       !
       END SELECT
+
+      IF( ln_dynldf_bgm )      CALL dyn_ldf_bgm( kt, ub, vb, ua, va)  ! bi-Laplacian GM parameterisation
 
       IF( l_trddyn ) THEN                        ! save the horizontal diffusive trends for further diagnostics
@@ -106,7 +107,7 @@
          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
+         IF( ln_dynldf_bgm )     WRITE(numout,*) '   ==>>>   biharmonic Gent-McWilliams operator turned on'
       ENDIF
       !
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 12775)
+++ NEMO/branches/UKMO/NEMO_4.0.1_biharmonic_GM/src/OCE/LDF/ldfdyn.F90	(revision 12776)
@@ -61,5 +61,4 @@
    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.
@@ -180,7 +179,6 @@
       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( (ioptio < 1).OR.(ioptio > 6) )   CALL ctl_stop( 'dyn_ldf_init: use at least ONE of the 4 operator options (NONE/lap/blp) but not (lap+blp)' )
       !
 
@@ -230,26 +228,6 @@
          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' )
          !
 
@@ -265,5 +243,4 @@
          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,*)
@@ -277,9 +254,6 @@
       !
       IF( ln_dynldf_OFF ) THEN
-         !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, bhm are not allocated'
-         !
-         RETURN
+         !
+         IF(lwp) WRITE(numout,*) '   ==>>>   No viscous operator selected. ahmt, ahmf are not allocated'
          !
       ELSE                                   !==  a lateral diffusion operator is used  ==!
@@ -292,9 +266,7 @@
          ahmt(:,:,:) = 0._wp                       ! init to 0 needed 
          ahmf(:,:,:) = 0._wp
-         bhm(:,:,:)  = 0._wp                       ! init to 0 needed 
          !                                        ! 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
@@ -372,142 +344,5 @@
             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
-                     bhm(ji,jj,:) = bhm0
-                  ENDDO
-               ENDDO
-
-               !     Surface BC : set to zero
-               bhm(:,:,1)   = 0._wp
-               !     Flat bottom BC : set to zero
-               bhm(:,:,jpk) = 0._wp 
-               !     Variable bathymetry case: diffusive fluxes masked in dyn_ldf_bgm
-
-               !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
-                     ikw = mbkt(ji,jj) !bottom last wet T-level
-                     DO jk = 2, ikw
-                        lhscrit=bhm0*rdt/(e3t_n(ji,jj,jk)**2*e1t(ji,jj)**2)
-                        IF(lhscrit .gt. lhscritmax) lhscritmax=lhscrit
-                     ENDDO
-                  ENDDO
-               ENDDO
-
-               IF ( lhscritmax .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
-                           bhm(ji,jj,jk) = bhm0 * gdepw_1d(jk)/rn_bgmzcrit
-                        ENDDO
-                     ENDDO
-
-                  ELSE
-
-                     ! constant at depth rn_bgmzcrit and larger
-                     DO jj = 2, jpjm1
-                        DO ji = 2, jpim1
-                           bhm(ji,jj,jk) = bhm0
-                        ENDDO
-                     ENDDO
-
-                  ENDIF
-
-               ENDDO
-
-               !     Surface BC : set to zero
-               bhm(:,:,1)   = 0._wp 
-               !     Flat bottom BC : set to zero
-               bhm(:,:,jpk) = 0._wp 
-               !     Variable bathymetry case: diffusive fluxes masked in dyn_ldf_bgm
-
-               !--------------------------------
-            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
-                        bhm(ji,jj,jk) = bhm0*e1t(ji,jj)**2*e3w_n(ji,jj,jk)**2
-                     ENDDO
-                  ENDDO
-               ENDDO
-
-               !     Surface BC : set to zero
-               bhm(:,:,1)   = 0._wp 
-               !     Flat bottom BC : set to zero
-               bhm(:,:,jpk) = 0._wp 
-               !     Variable bathymetry case: diffusive fluxes masked in dyn_ldf_bgm
-
-               !--------------------------------
-
-            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 
@@ -518,13 +353,155 @@
                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)
-               bhm(:,:,1:jpkm1) =       bhm(:,:,1:jpkm1)   * wmask(:,:,1:jpkm1)
-               CALL lbc_lnk( 'ldf_dyn_init', bhm , 'W', 1. ) 
-               !
             ENDIF
          ENDIF
          !
       ENDIF
+      !
+!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
+         IF(lwp) WRITE(numout,*) '  ==>>> Initialising bi-Laplacian GM eddy viscosity coefficient.'
+         ALLOCATE( bhm(jpi,jpj,jpk) , STAT=ierr )
+         bhm(:,:,:)  = 0._wp                       ! init to 0 needed 
+         bhm0 = rn_bhm_b       ;             ;   cl_Units = ' m4/s'   ! bi-Laplacian GM
+         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
+                  bhm(ji,jj,:) = bhm0
+               ENDDO
+            ENDDO
+
+            !     Surface BC : set to zero
+            bhm(:,:,1)   = 0._wp
+            !     Flat bottom BC : set to zero
+            bhm(:,:,jpk) = 0._wp 
+            !     Variable bathymetry case: diffusive fluxes masked in dyn_ldf_bgm
+
+            !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
+                  ikw = mbkt(ji,jj) !bottom last wet T-level
+                  DO jk = 2, ikw
+                     lhscrit=bhm0*rdt/(e3t_n(ji,jj,jk)**2*e1t(ji,jj)**2)
+                     IF(lhscrit .gt. lhscritmax) lhscritmax=lhscrit
+                  ENDDO
+               ENDDO
+            ENDDO
+
+            IF ( lhscritmax .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
+                        bhm(ji,jj,jk) = bhm0 * gdepw_1d(jk)/rn_bgmzcrit
+                     ENDDO
+                  ENDDO
+
+               ELSE
+
+                  ! constant at depth rn_bgmzcrit and larger
+                  DO jj = 2, jpjm1
+                     DO ji = 2, jpim1
+                        bhm(ji,jj,jk) = bhm0
+                     ENDDO
+                  ENDDO
+
+               ENDIF
+
+            ENDDO
+
+            !     Surface BC : set to zero
+            bhm(:,:,1)   = 0._wp 
+            !     Flat bottom BC : set to zero
+            bhm(:,:,jpk) = 0._wp 
+            !     Variable bathymetry case: diffusive fluxes masked in dyn_ldf_bgm
+
+            !--------------------------------
+         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
+                     bhm(ji,jj,jk) = bhm0*e1t(ji,jj)**2*e3w_n(ji,jj,jk)**2
+                  ENDDO
+               ENDDO
+            ENDDO
+
+            !     Surface BC : set to zero
+            bhm(:,:,1)   = 0._wp 
+            !     Flat bottom BC : set to zero
+            bhm(:,:,jpk) = 0._wp 
+            !     Variable bathymetry case: diffusive fluxes masked in dyn_ldf_bgm
+
+            !--------------------------------
+
+         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
+         !
+         ! mask and lbc_lnk
+         bhm(:,:,1:jpkm1) =       bhm(:,:,1:jpkm1)   * wmask(:,:,1:jpkm1)
+         CALL lbc_lnk( 'ldf_dyn_init', bhm , 'W', 1. ) 
+         !
+      ENDIF ! ln_dynldf_bgm
       !
     END SUBROUTINE ldf_dyn_init
