Index: /branches/UKMO/dev_r5518_bdy_sponge_temp/NEMOGCM/NEMO/OPA_SRC/BDY/bdy_oce.F90
===================================================================
--- /branches/UKMO/dev_r5518_bdy_sponge_temp/NEMOGCM/NEMO/OPA_SRC/BDY/bdy_oce.F90	(revision 5877)
+++ /branches/UKMO/dev_r5518_bdy_sponge_temp/NEMOGCM/NEMO/OPA_SRC/BDY/bdy_oce.F90	(revision 5878)
@@ -83,4 +83,5 @@
    LOGICAL                    ::   ln_mask_file             !: =T read bdymask from file
    LOGICAL                    ::   ln_vol                   !: =T volume correction             
+   LOGICAL, DIMENSION(jp_bdy) ::   ln_sponge                !: =T use sponge layer 
    !
    INTEGER                    ::   nb_bdy                   !: number of open boundary sets
@@ -88,4 +89,5 @@
    INTEGER                    ::   nn_volctl                !: = 0 the total volume will have the variability of the surface Flux E-P 
    !                                                        !  = 1 the volume will be constant during all the integration.
+   REAL(wp)                   ::   rn_sponge                !: multiplier of diffusion for sponge layer
    CHARACTER(len=20), DIMENSION(jp_bdy) ::   cn_dyn2d       ! Choice of boundary condition for barotropic variables (U,V,SSH)
    INTEGER, DIMENSION(jp_bdy)           ::   nn_dyn2d_dta   !: = 0 use the initial state as bdy dta ; 
@@ -118,7 +120,7 @@
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), TARGET ::   bdyumask   !: Mask defining computational domain at U-points
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), TARGET ::   bdyvmask   !: Mask defining computational domain at V-points
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), TARGET ::   sponge_factor !: Multiplier for diffusion for sponge layer
 
    REAL(wp)                                    ::   bdysurftot !: Lateral surface of unstructured open boundary
-
    !!----------------------------------------------------------------------
    !! open boundary data variables
@@ -148,5 +150,5 @@
       !
       ALLOCATE( bdytmask(jpi,jpj) , bdyumask(jpi,jpj), bdyvmask(jpi,jpj),     &  
-         &      STAT=bdy_oce_alloc )
+         &      sponge_factor(jpi,jpj), STAT=bdy_oce_alloc )
       !
       ! Initialize masks 
Index: /branches/UKMO/dev_r5518_bdy_sponge_temp/NEMOGCM/NEMO/OPA_SRC/BDY/bdyini.F90
===================================================================
--- /branches/UKMO/dev_r5518_bdy_sponge_temp/NEMOGCM/NEMO/OPA_SRC/BDY/bdyini.F90	(revision 5877)
+++ /branches/UKMO/dev_r5518_bdy_sponge_temp/NEMOGCM/NEMO/OPA_SRC/BDY/bdyini.F90	(revision 5878)
@@ -102,5 +102,6 @@
          &             cn_ice_lim, nn_ice_lim_dta,                           &
          &             rn_ice_tem, rn_ice_sal, rn_ice_age,                 &
-         &             ln_vol, nn_volctl, nn_rimwidth
+         &             ln_vol, nn_volctl, nn_rimwidth,                     &
+         &             ln_sponge, rn_sponge
       !!
       NAMELIST/nambdy_index/ ctypebdy, nbdyind, nbdybeg, nbdyend
@@ -365,4 +366,12 @@
         IF(lwp) WRITE(numout,*) '      Width of relaxation zone = ', nn_rimwidth(ib_bdy)
         IF(lwp) WRITE(numout,*)
+        IF( ln_sponge(ib_bdy) ) THEN                     ! check sponge layer choice
+             IF(lwp) WRITE(numout,*) 'Sponge layer applied at open boundaries'
+             IF(lwp) WRITE(numout,*) 'Multiplier for diffusion in sponge layer : ', rn_sponge
+             IF(lwp) WRITE(numout,*)
+           ELSE
+             IF(lwp) WRITE(numout,*) 'No volume correction applied at open boundaries'
+             IF(lwp) WRITE(numout,*)
+           ENDIF
 
       ENDDO
@@ -383,4 +392,6 @@
         ENDIF
      ENDIF
+
+     sponge_factor(:,:) = 1.0
 
       ! -------------------------------------------------
@@ -1103,4 +1114,18 @@
             END DO
          END DO 
+      ! Compute multiplier for diffusion for sponge layer
+           ! -------------------------------------------------
+           IF( ln_sponge(ib_bdy) ) THEN
+              igrd=1
+              DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd)
+                 nbi => idx_bdy(ib_bdy)%nbi(ib,igrd)
+                 nbj => idx_bdy(ib_bdy)%nbj(ib,igrd)
+                 nbr => idx_bdy(ib_bdy)%nbr(ib,igrd)
+                 sponge_factor(nbi,nbj) = 1.0 + (rn_sponge-1.0) * ( 1.- TANH( FLOAT( nbr - 1 ) *0.5 ) )
+                 ! 0.5 factor says how fast tanh goes to 1 - reduce this to have
+                 ! effect on more of the rimwidth
+              END DO
+           ENDIF
+
 
       ENDDO
Index: /branches/UKMO/dev_r5518_bdy_sponge_temp/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_iso.F90
===================================================================
--- /branches/UKMO/dev_r5518_bdy_sponge_temp/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_iso.F90	(revision 5877)
+++ /branches/UKMO/dev_r5518_bdy_sponge_temp/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_iso.F90	(revision 5878)
@@ -30,4 +30,7 @@
    USE wrk_nemo        ! Memory Allocation
    USE timing          ! Timing
+#if defined key_bdy 
+   USE bdy_oce         ! needed for extra diffusion in rim 
+#endif
 
    IMPLICIT NONE
@@ -114,4 +117,5 @@
       REAL(wp) ::   zcoef0, zcoef3, zcoef4, zmkt, zmkf               !   -      -
       REAL(wp) ::   zuav, zvav, zuwslpi, zuwslpj, zvwslpi, zvwslpj   !   -      -
+      REAL(wp), DIMENSION(jpi,jpj) :: zfactor  ! multiplier for diffusion
       !
       REAL(wp), POINTER, DIMENSION(:,:) :: ziut, zjuf, zjvt, zivf, zdku, zdk1u, zdkv, zdk1v
@@ -129,4 +133,14 @@
          IF( dyn_ldf_iso_alloc() /= 0 )   CALL ctl_stop('STOP', 'dyn_ldf_iso: failed to allocate arrays')
       ENDIF
+#if defined key_bdy
+      zfactor(:,:) = sponge_factor(:,:)
+#else
+      zfactor(:,:) = 1.0
+#endif
+#if defined key_bdy
+      zfactor(:,:) = sponge_factor(:,:)
+#else
+      zfactor(:,:) = 1.0
+#endif
 
       ! s-coordinate: Iso-level diffusion on tracer, but geopotential level diffusion on momentum
@@ -184,10 +198,10 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, jpi   ! vector opt.
-                  zabe1 = (fsahmt(ji,jj,jk)+ahmb0) * e2t(ji,jj) * MIN( fse3u(ji,jj,jk), fse3u(ji-1,jj,jk) ) / e1t(ji,jj)
+                  zabe1 = zfactor(ji,jj) *(fsahmt(ji,jj,jk)+ahmb0) * e2t(ji,jj) * MIN( fse3u(ji,jj,jk), fse3u(ji-1,jj,jk) ) / e1t(ji,jj)
 
                   zmskt = 1./MAX(  umask(ji-1,jj,jk  )+umask(ji,jj,jk+1)   &
                      &           + umask(ji-1,jj,jk+1)+umask(ji,jj,jk  ), 1. )
 
-                  zcof1 = - aht0 * e2t(ji,jj) * zmskt * 0.5  * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) )
+                  zcof1 = - zfactor(ji,jj) * aht0 * e2t(ji,jj) * zmskt * 0.5  * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) )
    
                   ziut(ji,jj) = (  zabe1 * ( ub(ji,jj,jk) - ub(ji-1,jj,jk) )   &
@@ -199,10 +213,10 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, jpi   ! vector opt.
-                  zabe1 = (fsahmt(ji,jj,jk)+ahmb0) * e2t(ji,jj) * fse3t(ji,jj,jk) / e1t(ji,jj)
+                  zabe1 = zfactor(ji,jj) * (fsahmt(ji,jj,jk)+ahmb0) * e2t(ji,jj) * fse3t(ji,jj,jk) / e1t(ji,jj)
 
                   zmskt = 1./MAX(  umask(ji-1,jj,jk  )+umask(ji,jj,jk+1)   &
                      &           + umask(ji-1,jj,jk+1)+umask(ji,jj,jk  ), 1. )
 
-                  zcof1 = - aht0 * e2t(ji,jj) * zmskt * 0.5  * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) )
+                  zcof1 = - zfactor(ji,jj) * aht0 * e2t(ji,jj) * zmskt * 0.5  * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) )
 
                   ziut(ji,jj) = (  zabe1 * ( ub(ji,jj,jk) - ub(ji-1,jj,jk) )   &
@@ -216,10 +230,10 @@
          DO jj = 1, jpjm1
             DO ji = 1, fs_jpim1   ! vector opt.
-               zabe2 = ( fsahmf(ji,jj,jk) + ahmb0 ) * e1f(ji,jj) * fse3f(ji,jj,jk) / e2f(ji,jj)
+               zabe2 = zfactor(ji,jj) * ( fsahmf(ji,jj,jk) + ahmb0 ) * e1f(ji,jj) * fse3f(ji,jj,jk) / e2f(ji,jj)
 
                zmskf = 1./MAX(  umask(ji,jj+1,jk  )+umask(ji,jj,jk+1)   &
                   &           + umask(ji,jj+1,jk+1)+umask(ji,jj,jk  ), 1. )
 
-               zcof2 = - aht0 * e1f(ji,jj) * zmskf * 0.5  * ( vslp(ji+1,jj,jk) + vslp(ji,jj,jk) )
+               zcof2 = - zfactor(ji,jj) * aht0 * e1f(ji,jj) * zmskf * 0.5  * ( vslp(ji+1,jj,jk) + vslp(ji,jj,jk) )
 
                zjuf(ji,jj) = (  zabe2 * ( ub(ji,jj+1,jk) - ub(ji,jj,jk) )   &
@@ -237,10 +251,10 @@
          DO jj = 2, jpjm1
             DO ji = 1, fs_jpim1   ! vector opt.
-               zabe1 = ( fsahmf(ji,jj,jk) + ahmb0 ) * e2f(ji,jj) * fse3f(ji,jj,jk) / e1f(ji,jj)
+               zabe1 = zfactor(ji,jj) * ( fsahmf(ji,jj,jk) + ahmb0 ) * e2f(ji,jj) * fse3f(ji,jj,jk) / e1f(ji,jj)
 
                zmskf = 1./MAX(  vmask(ji+1,jj,jk  )+vmask(ji,jj,jk+1)   &
                   &           + vmask(ji+1,jj,jk+1)+vmask(ji,jj,jk  ), 1. )
 
-               zcof1 = - aht0 * e2f(ji,jj) * zmskf * 0.5 * ( uslp(ji,jj+1,jk) + uslp(ji,jj,jk) )
+               zcof1 = - zfactor(ji,jj) * aht0 * e2f(ji,jj) * zmskf * 0.5 * ( uslp(ji,jj+1,jk) + uslp(ji,jj,jk) )
 
                zivf(ji,jj) = (  zabe1 * ( vb(ji+1,jj,jk) - vb(ji,jj,jk) )   &
@@ -254,10 +268,10 @@
             DO jj = 2, jpj
                DO ji = 1, fs_jpim1   ! vector opt.
-                  zabe2 = (fsahmt(ji,jj,jk)+ahmb0) * e1t(ji,jj) * MIN( fse3v(ji,jj,jk), fse3v(ji,jj-1,jk) ) / e2t(ji,jj)
+                  zabe2 = zfactor(ji,jj) * (fsahmt(ji,jj,jk)+ahmb0) * e1t(ji,jj) * MIN( fse3v(ji,jj,jk), fse3v(ji,jj-1,jk) ) / e2t(ji,jj)
 
                   zmskt = 1./MAX(  vmask(ji,jj-1,jk  )+vmask(ji,jj,jk+1)   &
                      &           + vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk  ), 1. )
 
-                  zcof2 = - aht0 * e1t(ji,jj) * zmskt * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) )
+                  zcof2 = - zfactor(ji,jj) * aht0 * e1t(ji,jj) * zmskt * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) )
 
                   zjvt(ji,jj) = (  zabe2 * ( vb(ji,jj,jk) - vb(ji,jj-1,jk) )   &
@@ -269,10 +283,10 @@
             DO jj = 2, jpj
                DO ji = 1, fs_jpim1   ! vector opt.
-                  zabe2 = (fsahmt(ji,jj,jk)+ahmb0) * e1t(ji,jj) * fse3t(ji,jj,jk) / e2t(ji,jj)
+                  zabe2 = zfactor(ji,jj) * (fsahmt(ji,jj,jk)+ahmb0) * e1t(ji,jj) * fse3t(ji,jj,jk) / e2t(ji,jj)
 
                   zmskt = 1./MAX(  vmask(ji,jj-1,jk  )+vmask(ji,jj,jk+1)   &
                      &           + vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk  ), 1. )
 
-                  zcof2 = - aht0 * e1t(ji,jj) * zmskt * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) )
+                  zcof2 = - zfactor(ji,jj) * aht0 * e1t(ji,jj) * zmskt * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) )
 
                   zjvt(ji,jj) = (  zabe2 * ( vb(ji,jj,jk) - vb(ji,jj-1,jk) )   &
Index: /branches/UKMO/dev_r5518_bdy_sponge_temp/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_iso.F90
===================================================================
--- /branches/UKMO/dev_r5518_bdy_sponge_temp/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_iso.F90	(revision 5877)
+++ /branches/UKMO/dev_r5518_bdy_sponge_temp/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_iso.F90	(revision 5878)
@@ -32,4 +32,7 @@
    USE wrk_nemo        ! Memory Allocation
    USE timing          ! Timing
+#if defined key_bdy 
+   USE bdy_oce 
+#endif 
 
    IMPLICIT NONE
@@ -110,4 +113,5 @@
       REAL(wp), POINTER, DIMENSION(:,:  ) ::  z2d
       REAL(wp), POINTER, DIMENSION(:,:,:) ::  zdkt, zdk1t, zdit, zdjt, ztfw 
+      REAL(wp), DIMENSION(jpi,jpj)     ::   zfactor  ! multiplier for diffusion
       !!----------------------------------------------------------------------
       !
@@ -124,4 +128,9 @@
       ENDIF
       !
+#if defined key_bdy
+      zfactor(:,:) = sponge_factor(:,:)
+#else
+      zfactor(:,:) = 1.0
+#endif
       !                                                          ! ===========
       DO jn = 1, kjpt                                            ! tracer loop
@@ -200,6 +209,6 @@
             DO jj = 1 , jpjm1
                DO ji = 1, fs_jpim1   ! vector opt.
-                  zabe1 = ( fsahtu(ji,jj,jk) + pahtb0 ) * re2u_e1u(ji,jj) * fse3u_n(ji,jj,jk)
-                  zabe2 = ( fsahtv(ji,jj,jk) + pahtb0 ) * re1v_e2v(ji,jj) * fse3v_n(ji,jj,jk)
+                  zabe1 = zfactor(ji,jj) * ( fsahtu(ji,jj,jk) + pahtb0 ) * re2u_e1u(ji,jj) * fse3u_n(ji,jj,jk)
+                  zabe2 = zfactor(ji,jj) * ( fsahtv(ji,jj,jk) + pahtb0 ) * re1v_e2v(ji,jj) * fse3v_n(ji,jj,jk)
                   !
                   zmsku = 1. / MAX(  tmask(ji+1,jj,jk  ) + tmask(ji,jj,jk+1)   &
Index: /branches/UKMO/dev_r5518_bdy_sponge_temp/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_lap.F90
===================================================================
--- /branches/UKMO/dev_r5518_bdy_sponge_temp/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_lap.F90	(revision 5877)
+++ /branches/UKMO/dev_r5518_bdy_sponge_temp/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_lap.F90	(revision 5878)
@@ -26,4 +26,7 @@
    USE lib_mpp         ! MPP library
    USE timing          ! Timing
+#if defined key_bdy 
+   USE bdy_oce 
+#endif
 
    IMPLICIT NONE
@@ -77,4 +80,5 @@
       INTEGER  ::   iku, ikv, ierr       ! local integers
       REAL(wp) ::   zabe1, zabe2, zbtr   ! local scalars
+      REAL(wp), DIMENSION(jpi,jpj)     ::   zfactor  ! multiplier for diffusion
       !!----------------------------------------------------------------------
       !
@@ -86,4 +90,9 @@
          IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ '
       ENDIF
+#if defined key_bdy
+     zfactor(:,:) = sponge_factor(:,:)
+#else
+     zfactor(:,:) = 1.0
+#endif
 
       !                                                          ! =========== !
@@ -96,6 +105,6 @@
             DO jj = 1, jpjm1
                DO ji = 1, fs_jpim1   ! vector opt.
-                  zabe1 = fsahtu(ji,jj,jk) * umask(ji,jj,jk) * re2u_e1u(ji,jj) * fse3u_n(ji,jj,jk)
-                  zabe2 = fsahtv(ji,jj,jk) * vmask(ji,jj,jk) * re1v_e2v(ji,jj) * fse3v_n(ji,jj,jk)
+                  zabe1 = zfactor(ji,jj) * fsahtu(ji,jj,jk) * umask(ji,jj,jk) * re2u_e1u(ji,jj) * fse3u_n(ji,jj,jk)
+                  zabe2 = zfactor(ji,jj) * fsahtv(ji,jj,jk) * vmask(ji,jj,jk) * re1v_e2v(ji,jj) * fse3v_n(ji,jj,jk)
                   ztu(ji,jj,jk) = zabe1 * ( ptb(ji+1,jj  ,jk,jn) - ptb(ji,jj,jk,jn) )
                   ztv(ji,jj,jk) = zabe2 * ( ptb(ji  ,jj+1,jk,jn) - ptb(ji,jj,jk,jn) )
@@ -109,9 +118,9 @@
                      ikv = mbkv(ji,jj)
                      IF( iku == jk ) THEN
-                        zabe1 = fsahtu(ji,jj,iku) * umask(ji,jj,iku) * re2u_e1u(ji,jj) * fse3u_n(ji,jj,iku)
+                        zabe1 = zfactor(ji,jj) * fsahtu(ji,jj,iku) * umask(ji,jj,iku) * re2u_e1u(ji,jj) * fse3u_n(ji,jj,iku)
                         ztu(ji,jj,jk) = zabe1 * pgu(ji,jj,jn)
                      ENDIF
                      IF( ikv == jk ) THEN
-                        zabe2 = fsahtv(ji,jj,ikv) * vmask(ji,jj,ikv) * re1v_e2v(ji,jj) * fse3v_n(ji,jj,ikv)
+                        zabe2 = zfactor(ji,jj) * fsahtv(ji,jj,ikv) * vmask(ji,jj,ikv) * re1v_e2v(ji,jj) * fse3v_n(ji,jj,ikv)
                         ztv(ji,jj,jk) = zabe2 * pgv(ji,jj,jn)
                      ENDIF
@@ -128,9 +137,9 @@
                      ikv = mikv(ji,jj) 
                      IF( iku == MAX(2,jk) ) THEN 
-                        zabe1 = fsahtu(ji,jj,iku) * umask(ji,jj,iku) * re2u_e1u(ji,jj) * fse3u_n(ji,jj,iku) 
+                        zabe1 = zfactor(ji,jj) * fsahtu(ji,jj,iku) * umask(ji,jj,iku) * re2u_e1u(ji,jj) * fse3u_n(ji,jj,iku) 
                         ztu(ji,jj,jk) = zabe1 * pgui(ji,jj,jn) 
                      ENDIF 
                      IF( ikv == MAX(2,jk) ) THEN 
-                        zabe2 = fsahtv(ji,jj,ikv) * vmask(ji,jj,ikv) * re1v_e2v(ji,jj) * fse3v_n(ji,jj,ikv) 
+                        zabe2 = zfactor(ji,jj) * fsahtv(ji,jj,ikv) * vmask(ji,jj,ikv) * re1v_e2v(ji,jj) * fse3v_n(ji,jj,ikv) 
                         ztv(ji,jj,jk) = zabe2 * pgvi(ji,jj,jn) 
                      END IF 
