Index: /branches/dev_r1821_mixed_ldfdyn/NEMO/OPA_SRC/DYN/dynldf.F90
===================================================================
--- /branches/dev_r1821_mixed_ldfdyn/NEMO/OPA_SRC/DYN/dynldf.F90	(revision 1838)
+++ /branches/dev_r1821_mixed_ldfdyn/NEMO/OPA_SRC/DYN/dynldf.F90	(revision 1839)
@@ -69,4 +69,10 @@
       CASE ( 2 )    ;   CALL dyn_ldf_bilap  ( kt )      ! iso-level bilaplacian
       CASE ( 3 )    ;   CALL dyn_ldf_bilapg ( kt )      ! s-coord. horizontal bilaplacian
+      CASE ( 4 )                                        ! iso-level laplacian + bilaplacian
+         CALL dyn_ldf_lap    ( kt )
+         CALL dyn_ldf_bilap  ( kt )
+      CASE ( 5 )                                        ! rotated laplacian + bilaplacian (s-coord)
+         CALL dyn_ldf_iso    ( kt )
+         CALL dyn_ldf_bilapg ( kt )
       !
       CASE ( -1 )                                       ! esopa: test all possibility with control print
@@ -83,4 +89,11 @@
                         CALL prt_ctl( tab3d_1=ua, clinfo1=' ldf3 - Ua: ', mask1=umask,   &
             &                         tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
+      !
+      CASE ( -2 )                                       ! neither laplacian nor bilaplacian schemes used
+         IF( kt == nit000 ) THEN
+            IF(lwp) WRITE(numout,*)
+            IF(lwp) WRITE(numout,*) 'dyn_ldf : no lateral diffusion on momentum setup'
+            IF(lwp) WRITE(numout,*) '~~~~~~~ '
+         ENDIF
       END SELECT
 
@@ -124,5 +137,5 @@
       IF( ln_dynldf_lap   )   ioptio = ioptio + 1
       IF( ln_dynldf_bilap )   ioptio = ioptio + 1
-      IF( ioptio /= 1 ) CALL ctl_stop( '          use ONE of the 2 lap/bilap operator type on dynamics' )
+      IF( ioptio <  1 ) CALL ctl_warn( '          neither laplacian nor bilaplacian operator set for dynamics' )
       ioptio = 0
       IF( ln_dynldf_level )   ioptio = ioptio + 1
@@ -144,5 +157,5 @@
             IF ( ln_dynldf_iso   )   nldf = 1      ! isoneutral (   rotation)
          ENDIF
-         IF ( ln_sco ) THEN             ! z-coordinate
+         IF ( ln_sco ) THEN             ! s-coordinate
             IF ( ln_dynldf_level )   nldf = 0      ! iso-level  (no rotation)
             IF ( ln_dynldf_hor   )   nldf = 1      ! horizontal (   rotation)
@@ -162,5 +175,5 @@
             IF ( ln_dynldf_iso   )   ierr = 2      ! isoneutral (   rotation)
          ENDIF
-         IF ( ln_sco ) THEN             ! z-coordinate
+         IF ( ln_sco ) THEN             ! s-coordinate
             IF ( ln_dynldf_level )   nldf = 2      ! iso-level  (no rotation)
             IF ( ln_dynldf_hor   )   nldf = 3      ! horizontal (   rotation)
@@ -169,4 +182,22 @@
       ENDIF
       
+      IF( ln_dynldf_lap .AND. ln_dynldf_bilap ) THEN  ! mixed laplacian and bilaplacian operators
+         IF ( ln_zco ) THEN                ! z-coordinate
+            IF ( ln_dynldf_level )   nldf = 4      ! iso-level  (no rotation)
+            IF ( ln_dynldf_hor   )   nldf = 4      ! horizontal (no rotation)
+            IF ( ln_dynldf_iso   )   ierr = 2      ! isoneutral (   rotation)
+         ENDIF
+         IF ( ln_zps ) THEN             ! z-coordinate
+            IF ( ln_dynldf_level )   ierr = 1      ! iso-level not allowed 
+            IF ( ln_dynldf_hor   )   nldf = 4      ! horizontal (no rotation)
+            IF ( ln_dynldf_iso   )   ierr = 2      ! isoneutral (   rotation)
+         ENDIF
+         IF ( ln_sco ) THEN             ! s-coordinate
+            IF ( ln_dynldf_level )   nldf = 4      ! iso-level  (no rotation)
+            IF ( ln_dynldf_hor   )   nldf = 5      ! horizontal (   rotation)
+            IF ( ln_dynldf_iso   )   ierr = 2      ! isoneutral (   rotation)
+         ENDIF
+      ENDIF
+
       IF( lk_esopa )                 nldf = -1     ! esopa test
 
@@ -179,9 +210,12 @@
       IF(lwp) THEN
          WRITE(numout,*)
+         IF( nldf == -2 )   WRITE(numout,*) '              neither laplacian nor bilaplacian schemes used'
          IF( nldf == -1 )   WRITE(numout,*) '              ESOPA test All scheme used'
          IF( nldf ==  0 )   WRITE(numout,*) '              laplacian operator'
-         IF( nldf ==  1 )   WRITE(numout,*) '              Rotated laplacian operator'
+         IF( nldf ==  1 )   WRITE(numout,*) '              rotated laplacian operator'
          IF( nldf ==  2 )   WRITE(numout,*) '              bilaplacian operator'
-         IF( nldf ==  3 )   WRITE(numout,*) '              Rotated bilaplacian'
+         IF( nldf ==  3 )   WRITE(numout,*) '              rotated bilaplacian'
+         IF( nldf ==  4 )   WRITE(numout,*) '              laplacian and bilaplacian operators'
+         IF( nldf ==  5 )   WRITE(numout,*) '              rotated laplacian and bilaplacian operators'
       ENDIF
       !
Index: /branches/dev_r1821_mixed_ldfdyn/NEMO/OPA_SRC/DYN/dynldf_lap.F90
===================================================================
--- /branches/dev_r1821_mixed_ldfdyn/NEMO/OPA_SRC/DYN/dynldf_lap.F90	(revision 1838)
+++ /branches/dev_r1821_mixed_ldfdyn/NEMO/OPA_SRC/DYN/dynldf_lap.F90	(revision 1839)
@@ -78,5 +78,5 @@
       IF( kt == nit000 ) THEN
          IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) 'dyn_ldf : iso-level harmonic (laplacien) operator'
+         IF(lwp) WRITE(numout,*) 'dyn_ldf : iso-level harmonic (laplacian) operator'
          IF(lwp) WRITE(numout,*) '~~~~~~~ '
       ENDIF
Index: /branches/dev_r1821_mixed_ldfdyn/NEMO/OPA_SRC/LDF/ldfdyn.F90
===================================================================
--- /branches/dev_r1821_mixed_ldfdyn/NEMO/OPA_SRC/LDF/ldfdyn.F90	(revision 1838)
+++ /branches/dev_r1821_mixed_ldfdyn/NEMO/OPA_SRC/LDF/ldfdyn.F90	(revision 1839)
@@ -67,5 +67,5 @@
       NAMELIST/namdyn_ldf/ ln_dynldf_lap  , ln_dynldf_bilap,                  &
          &                 ln_dynldf_level, ln_dynldf_hor  , ln_dynldf_iso,   &
-         &                 rn_ahm_0       , rn_ahmb_0
+         &                 rn_ahm_0       , rn_ahmb_0      , rn_ahm_0_blp
       !!----------------------------------------------------------------------
 
@@ -78,15 +78,18 @@
          WRITE(numout,*) '~~~~~~~'
          WRITE(numout,*) '   Namelist nam_dynldf : set lateral mixing parameters'
-         WRITE(numout,*) '      laplacian operator          ln_dynldf_lap   = ', ln_dynldf_lap
-         WRITE(numout,*) '      bilaplacian operator        ln_dynldf_bilap = ', ln_dynldf_bilap
-         WRITE(numout,*) '      iso-level                   ln_dynldf_level = ', ln_dynldf_level
-         WRITE(numout,*) '      horizontal (geopotential)   ln_dynldf_hor   = ', ln_dynldf_hor
-         WRITE(numout,*) '      iso-neutral                 ln_dynldf_iso   = ', ln_dynldf_iso
-         WRITE(numout,*) '      horizontal eddy viscosity   rn_ahm_0        = ', rn_ahm_0
-         WRITE(numout,*) '      background viscosity        rn_ahmb_0       = ', rn_ahmb_0
-      ENDIF
-
-      ahm0  = rn_ahm_0                  ! OLD namelist variables defined from DOCTOR namelist variables
-      ahmb0 = rn_ahmb_0
+         WRITE(numout,*) '      laplacian operator                      ln_dynldf_lap   = ', ln_dynldf_lap
+         WRITE(numout,*) '      bilaplacian operator                    ln_dynldf_bilap = ', ln_dynldf_bilap
+         WRITE(numout,*) '      iso-level                               ln_dynldf_level = ', ln_dynldf_level
+         WRITE(numout,*) '      horizontal (geopotential)               ln_dynldf_hor   = ', ln_dynldf_hor
+         WRITE(numout,*) '      iso-neutral                             ln_dynldf_iso   = ', ln_dynldf_iso
+         WRITE(numout,*) '      horizontal laplacian eddy viscosity     rn_ahm_0        = ', rn_ahm_0
+         WRITE(numout,*) '      background viscosity                    rn_ahmb_0       = ', rn_ahmb_0
+         WRITE(numout,*) '      horizontal bilaplacian eddy viscosity   rn_ahm_0        = ', rn_ahm_0
+
+      ENDIF
+
+      ahm0     = rn_ahm_0                  ! OLD namelist variables defined from DOCTOR namelist variables
+      ahmb0    = rn_ahmb_0
+      ahm0_blp = rn_ahm_0_blp
 
       ! ... check of lateral diffusive operator on tracers
@@ -117,5 +120,7 @@
       IF( ln_dynldf_bilap ) THEN
          IF(lwp) WRITE(numout,*) '   biharmonic momentum diffusion'
-         IF( ahm0 > 0 .AND. .NOT. lk_esopa )   CALL ctl_stop( 'The horizontal viscosity coef. ahm0 must be negative' )
+         IF( ahm0_blp == 0.0 ) ahm0_blp = ahm0       ! Old namelist method: bilap specified with ahm0
+         IF( .NOT. ln_dynldf_lap ) ahm0 = ahm0_blp   ! Allow spatially varying coefs, which use ahm0 as input
+         IF( ahm0_blp > 0 .AND. .NOT. lk_esopa )   CALL ctl_stop( 'The horizontal viscosity coef. ahm0 must be negative' )
       ELSE
          IF(lwp) WRITE(numout,*) '   harmonic momentum diff. (default)'
Index: /branches/dev_r1821_mixed_ldfdyn/NEMO/OPA_SRC/LDF/ldfdyn_c1d.h90
===================================================================
--- /branches/dev_r1821_mixed_ldfdyn/NEMO/OPA_SRC/LDF/ldfdyn_c1d.h90	(revision 1838)
+++ /branches/dev_r1821_mixed_ldfdyn/NEMO/OPA_SRC/LDF/ldfdyn_c1d.h90	(revision 1839)
@@ -17,5 +17,6 @@
       !! ** Method  :   1D eddy viscosity coefficients ( depth )
       !!      ahm3, ahm4 never used
-      !!      biharmonic or harmonic operator : ahm1=ahm2 defined at T-level
+      !!      harmonic   operator : ahm1 defined at T-level
+      !!      biharmonic operator : ahm2 defined at T-level
       !!      isopycnal or geopotential harmonic operator
       !!                                      : ahm1 defined at T-level
@@ -27,5 +28,6 @@
 
       !! * Local variables
-      REAL(wp) ::   zdam, zwam, zm00, zm01, zmhf, zmhs
+      REAL(wp) ::    zdam,  zwam,  zm00,  zm01,  zmhf,  zmhs
+      REAL(wp) ::   zdam2, zwam2, zm200, zm201, zmh2f, zmh2s
       REAL(wp) ::   zahmf, zahms
       !!----------------------------------------------------------------------
@@ -36,5 +38,5 @@
       IF(lwp) WRITE(numout,*)
 
-      ! Set ahm1=ahm2   (always at t-level)
+      ! Set ahm1 for laplacian     (always at t-level)
       ! =============
       ! (USER: modify ahm1 following your desiderata)
@@ -53,8 +55,26 @@
       zmhs = zahms-zmhf * zm00
 
-      ! set ahm1=ahm2 at T-level
+      ! Set ahm2 for bilaplacian   (always at t-level)
+      ! =============
+      ! (USER: modify ahm2 following your desiderata)
+
+      ! initialization of the profile
+      ! ahms, ahmf: surface and bottom values
+      zahm2s = ahm0_blp
+      zahm2f = ahm0_blp/4.
+      ! zdam, zwam: depth of the inflection pt and width of inflection
+      zdam2  = -300.
+      zwam2  =  300.
+      ! computation coefficients
+      zm200 = TANH( (0-zdam2)/zwam2 )
+      zm201 = TANH( (-fsdept(1,1,jpk)-zdam2)/zwam2 )
+      zmh2f = (zahm2s-zahm2f)/(zm200-zm201)
+      zmh2s = zahm2s-zmh2f * zm00
+
+
+      ! set ahm1 and ahm2 at T-level
       DO jk = 1, jpk
-         ahm1(jk) = zmhs + zmhf * TANH( (-fsdept(1,1,jk)-zdam) / zwam )
-         ahm2(jk) = ahm1(jk)
+         ahm1(jk) = zmhs  + zmhf  * TANH( (-fsdept(1,1,jk)-zdam ) / zwam  )
+         ahm2(jk) = zmh2s + zmh2f * TANH( (-fsdept(1,1,jk)-zdam2) / zwam2 )
       END DO
 
@@ -62,9 +82,16 @@
       IF(lwp .AND. ld_print ) THEN
          WRITE(numout,*)
-         WRITE(numout,*) '         ahm profile : '
+         WRITE(numout,*) '         ahm profile (laplacian): '
          WRITE(numout,*)
          WRITE(numout,9100)
          DO jk = 1, jpk
             WRITE(numout,9110) jk, ahm1(jk), fsdept(1,1,jk)
+         END DO
+         WRITE(numout,*)
+         WRITE(numout,*) '         ahm profile (bilaplacian): '
+         WRITE(numout,*)
+         WRITE(numout,9100)
+         DO jk = 1, jpk
+            WRITE(numout,9110) jk, ahm2(jk), fsdept(1,1,jk)
          END DO
       ENDIF
Index: /branches/dev_r1821_mixed_ldfdyn/NEMO/OPA_SRC/LDF/ldfdyn_c2d.h90
===================================================================
--- /branches/dev_r1821_mixed_ldfdyn/NEMO/OPA_SRC/LDF/ldfdyn_c2d.h90	(revision 1838)
+++ /branches/dev_r1821_mixed_ldfdyn/NEMO/OPA_SRC/LDF/ldfdyn_c2d.h90	(revision 1839)
@@ -27,7 +27,7 @@
       !!           iso-model level : ahm3, ahm4 not used
       !!
-      !!       biharmonic operator : ahm1 is defined at u-point
-      !!                             ahm2 is defined at v-point
-      !!                           : ahm3, ahm4 not used
+      !!       biharmonic operator : ahm3 is defined at u-point
+      !!                             ahm4 is defined at v-point
+      !!                           : ahm1, ahm2 not used
       !!
       !!----------------------------------------------------------------------
@@ -102,5 +102,5 @@
          IF(lwp) WRITE(numout,*) '              maximum grid-spacing = ', zd_max, ' maximum value for ahm = ', ahm0
 
-         za00 = ahm0 / ( zd_max * zd_max * zd_max )
+         za00 = ahm0_blp / ( zd_max * zd_max * zd_max )
          DO jj = 1, jpj
             DO ji = 1, jpi
Index: /branches/dev_r1821_mixed_ldfdyn/NEMO/OPA_SRC/LDF/ldfdyn_c3d.h90
===================================================================
--- /branches/dev_r1821_mixed_ldfdyn/NEMO/OPA_SRC/LDF/ldfdyn_c3d.h90	(revision 1838)
+++ /branches/dev_r1821_mixed_ldfdyn/NEMO/OPA_SRC/LDF/ldfdyn_c3d.h90	(revision 1839)
@@ -122,5 +122,5 @@
          IF(lwp) WRITE(numout,*) '              maximum grid-spacing = ', zd_max, ' maximum value for ahm = ', ahm0
 
-         za00 = ahm0 / ( zd_max * zd_max * zd_max )
+         za00 = ahm0_blp / ( zd_max * zd_max * zd_max )
          DO jj = 1, jpj
             DO ji = 1, jpi
Index: /branches/dev_r1821_mixed_ldfdyn/NEMO/OPA_SRC/LDF/ldfdyn_oce.F90
===================================================================
--- /branches/dev_r1821_mixed_ldfdyn/NEMO/OPA_SRC/LDF/ldfdyn_oce.F90	(revision 1838)
+++ /branches/dev_r1821_mixed_ldfdyn/NEMO/OPA_SRC/LDF/ldfdyn_oce.F90	(revision 1839)
@@ -17,8 +17,9 @@
    LOGICAL , PUBLIC ::   ln_dynldf_hor   = .TRUE.      !: horizontal (geopotential) direction
    LOGICAL , PUBLIC ::   ln_dynldf_iso   = .FALSE.     !: iso-neutral direction
-   REAL(wp), PUBLIC ::   rn_ahm_0        = 40000._wp   !: lateral eddy viscosity (m2/s)
-   REAL(wp), PUBLIC ::   rn_ahmb_0       =     0._wp   !: lateral background eddy viscosity (m2/s)
+   REAL(wp), PUBLIC ::   rn_ahm_0        = 40000._wp   !: lateral laplacian eddy viscosity (m2/s)
+   REAL(wp), PUBLIC ::   rn_ahmb_0       =     0._wp   !: lateral laplacian background eddy viscosity (m2/s)
+   REAL(wp), PUBLIC ::   rn_ahm_0_blp    =     0._wp   !: lateral bilaplacian eddy viscosity (m4/s)
 
-   REAL(wp), PUBLIC ::   ahm0, ahmb0                 ! OLD namelist names
+   REAL(wp), PUBLIC ::   ahm0, ahmb0, ahm0_blp       ! OLD namelist names
 
 #if defined key_dynldf_c3d
Index: /branches/dev_r1821_mixed_ldfdyn/NEMO/OPA_SRC/LDF/ldfdyn_substitute.h90
===================================================================
--- /branches/dev_r1821_mixed_ldfdyn/NEMO/OPA_SRC/LDF/ldfdyn_substitute.h90	(revision 1838)
+++ /branches/dev_r1821_mixed_ldfdyn/NEMO/OPA_SRC/LDF/ldfdyn_substitute.h90	(revision 1839)
@@ -10,4 +10,8 @@
    !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 
    !!----------------------------------------------------------------------
+   !!
+   !! fsahmt, fsahmf - used for laplaian operator only
+   !! fsahmu, fsahmv - used for bilaplacian operator only
+   !!
 #if defined key_dynldf_c3d
 !   ' key_dynldf_c3d' :                  3D coefficient
@@ -26,11 +30,11 @@
 #  define   fsahmt(i,j,k)    ahm1(k)
 #  define   fsahmf(i,j,k)    ahm1(k)
-#  define   fsahmu(i,j,k)    ahm1(k)
-#  define   fsahmv(i,j,k)    ahm1(k)
+#  define   fsahmu(i,j,k)    ahm2(k)
+#  define   fsahmv(i,j,k)    ahm2(k)
 #else
 !   default option :               Constant coefficient
 #  define   fsahmt(i,j,k)   ahm0
 #  define   fsahmf(i,j,k)   ahm0
-#  define   fsahmu(i,j,k)   ahm0
-#  define   fsahmv(i,j,k)   ahm0
+#  define   fsahmu(i,j,k)   ahm0_blp
+#  define   fsahmv(i,j,k)   ahm0_blp
 #endif
