Index: branches/UKMO/dev_r5518_convadj/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_iso.F90
===================================================================
--- branches/UKMO/dev_r5518_convadj/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_iso.F90	(revision 6836)
+++ branches/UKMO/dev_r5518_convadj/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_iso.F90	(revision 6837)
@@ -32,4 +32,6 @@
    USE wrk_nemo        ! Memory Allocation
    USE timing          ! Timing
+   USE zdfconvadj
+   USE zdfmxl, ONLY : nmln
 
    IMPLICIT NONE
@@ -103,4 +105,5 @@
       REAL(wp)                             , INTENT(in   ) ::   pahtb0     ! background diffusion coef
       !
+      INTEGER  :: zconvadjw, zconvadju, zconvadjv 
       INTEGER  ::  ji, jj, jk, jn   ! dummy loop indices
       INTEGER  ::  ikt
@@ -212,8 +215,15 @@
                   zcof2 = - fsahtv(ji,jj,jk) * e1v(ji,jj) * vslp(ji,jj,jk) * zmskv
                   !
-                  zftu(ji,jj,jk ) = (  zabe1 * zdit(ji,jj,jk)   &
+                  IF(ln_convadj .AND. (jk .GT. nmln(ji,jj)) ) THEN
+                     zconvadju = 1 - MAX( convadj_mask(ji,jj), convadj_mask(ji+1,jj) )
+                     zconvadjv = 1 - MAX( convadj_mask(ji,jj), convadj_mask(ji,jj+1) )
+                  ELSE
+                     zconvadju = 1
+                     zconvadjv = 1
+                  END IF
+                  zftu(ji,jj,jk ) = zconvadju * (  zabe1 * zdit(ji,jj,jk)   &
                      &              + zcof1 * (  zdkt (ji+1,jj,jk) + zdk1t(ji,jj,jk)      &
                      &                         + zdk1t(ji+1,jj,jk) + zdkt (ji,jj,jk)  )  ) * umask(ji,jj,jk)
-                  zftv(ji,jj,jk) = (  zabe2 * zdjt(ji,jj,jk)   &
+                  zftv(ji,jj,jk) = zconvadjv * (  zabe2 * zdjt(ji,jj,jk)   &
                      &              + zcof2 * (  zdkt (ji,jj+1,jk) + zdk1t(ji,jj,jk)      &
                      &                         + zdk1t(ji,jj+1,jk) + zdkt (ji,jj,jk)  )  ) * vmask(ji,jj,jk)                  
@@ -296,6 +306,11 @@
                      &            + vmask(ji,jj-1,jk-1) + vmask(ji,jj  ,jk), 1.  )
                   !
-                  zcoef3 = zcoef0 * e2t(ji,jj) * zmsku * wslpi (ji,jj,jk)
-                  zcoef4 = zcoef0 * e1t(ji,jj) * zmskv * wslpj (ji,jj,jk)
+                  IF(ln_convadj .AND. ( jk .GT. nmln(ji,jj)) ) THEN
+                     zconvadjw = 1 - convadj_mask(ji,jj)
+                  ELSE
+                     zconvadjw = 1
+                  ENDIF
+                  zcoef3 = zcoef0 * e2t(ji,jj) * zmsku * wslpi (ji,jj,jk) * zconvadjw
+                  zcoef4 = zcoef0 * e1t(ji,jj) * zmskv * wslpj (ji,jj,jk) * zconvadjw
                   !
                   ztfw(ji,jj,jk) = zcoef3 * (   zdit(ji  ,jj  ,jk-1) + zdit(ji-1,jj  ,jk)      &
Index: branches/UKMO/dev_r5518_convadj/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfconvadj.F90
===================================================================
--- branches/UKMO/dev_r5518_convadj/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfconvadj.F90	(revision 6837)
+++ branches/UKMO/dev_r5518_convadj/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfconvadj.F90	(revision 6837)
@@ -0,0 +1,308 @@
+MODULE zdfconvadj
+   !!======================================================================
+   !!                       ***  MODULE  zdfconvadj  ***
+   !! Perform a convective adjustment of an unstable water column at overflow
+   !! grid points.
+   !!======================================================================
+   !! History :  3.4  !  2013-08  (T. Graham)  Original code adapted from UKMO HadCM3 code
+   !!
+   !!----------------------------------------------------------------------
+   USE oce             ! ocean dynamics and tracers variables
+   USE dom_oce         ! ocean space and time domain variables
+   USE zdf_oce         ! ocean vertical physics variables
+   USE in_out_manager  ! I/O manager
+   USE iom             ! for iom_put
+   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
+   USE timing          ! Timing
+   USE wrk_nemo        ! Wrk arrays
+   USE zdfmxl, ONLY : nmln
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   zdf_convadj         ! called by step.F90
+   PUBLIC   zdf_convadj_init    ! called by nemogcm.f90
+
+   ! 
+   ! Namelist zdf_conadj
+   INTEGER nn_con                      ! Number of times to convect
+   CHARACTER(300) cn_convadj_mask      ! NetCDF mask file indicating points
+   LOGICAL , PUBLIC ::   ln_zdfconvadj ! convadj flag - used in step and tra_ldfiso
+   
+   REAL(wp) , PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:) :: convadj_mask !: mask of points where adjustment applied
+   
+
+   !! * Substitutions
+#  include "domzgr_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 4.0 , NEMO Consortium (2011)
+   !! $Id$
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+
+CONTAINS
+
+   SUBROUTINE zdf_convadj( kt )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE zdf_convadj  ***
+      !!                 
+      !! ** Purpose: Perform a convective adjustment of an unstable water column
+      !!                
+      !! ** Method: Given an unstable profile, the dense water will be moved
+      !!                down the water column to a depth at which it is stable, and the
+      !!                intervening boxes shuffled upwards, rather than using convection
+      !!                to stabilise the water column.
+      !!                Note that this requires a "cliff" to be dug out in bathymetry at
+      !!                overflow gridpoints.
+      !!                NB: This scheme modifies tsn in place to remove the
+      !!                    instabilities in these regions before other physics is
+      !!                    called.
+      !!
+      !! ** Action  :   T and S updated in static instability cases
+      !!----------------------------------------------------------------------
+      USE par_oce, ONLY: jp_tem, jp_sal, jpts
+
+      INTEGER, INTENT( in ) ::   kt   ! ocean time-step indexocean time step
+      INTEGER  ::   ji, jj, jk, jn, icon           ! dummy loop indices  
+      INTEGER  ::   itop, icona, iconb, ilevel
+      REAL(wp) , POINTER, DIMENSION(:) :: rhoup, rholo, trasum, tra
+      REAL(wp) ::   dztdif, tramix, rho, rho_insitu
+
+      IF( nn_timing == 1 )  CALL timing_start('zdf_convadj')
+      !
+      IF( kt == nit000 ) THEN
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) 'zdf_convadj : Convective adjustment at overflow points'
+         IF(lwp) WRITE(numout,*) '~~~~~~~ '
+         IF(lwp) WRITE(numout,*)
+      ENDIF
+
+      !Allocation of arrays
+      CALL wrk_alloc( jpkm1, rhoup )
+      CALL wrk_alloc( jpkm1, rholo )
+      CALL wrk_alloc( jpts, tra )
+      CALL wrk_alloc( jpts, trasum )
+      
+      DO ji = 1,jpi !Do we need to include halos here?
+        DO jj = 1,jpj
+          IF (convadj_mask(ji,jj) .EQ. 1) THEN
+            DO icon = 1,nn_con ! Loop to repeat action nn_con times
+              rhoup(:)=0 
+              rholo(:)=0
+              itop=nmln(ji,jj)+1 ! Only apply below mixed layer
+              DO WHILE (itop .lt. mbkt(ji,jj))
+                ! Loop over depth and calculate potential density at W points (not surface)
+                ! *Do we need to worry about the before density or just the now?
+                DO jk = nmln(ji,jj),mbkt(ji,jj)-1
+                  rhoup(jk) = sigmai( tsn(ji,jj,jk,jp_tem), tsn(ji,jj,jk,jp_sal), fsdepw(ji,jj,jk) )
+                  rholo(jk) = sigmai( tsn(ji,jj,jk+1,jp_tem), tsn(ji,jj,jk+1,jp_sal), fsdepw(ji,jj,jk) )
+                END DO
+  
+                ! Loop over depth from bottom to find shallowest 
+                ! depth where column is unstable
+                icona = 0
+                DO jk = mbkt(ji,jj),itop,-1
+                  IF ( rhoup(jk) > rholo(jk) ) THEN
+                    icona = jk
+                    iconb = jk+1
+                  ENDIF
+                END DO
+
+                ! Loop over depth to find depth where this water is stable.
+                ! *Need to deal with key_vvl here.
+                IF (icona .NE. 0) THEN
+                  jk=iconb+1
+                  rho=1
+                  rho_insitu=0
+                  ilevel=iconb
+                  DO WHILE ((jk .LE. mbkt(ji,jj)) .AND. (rho .GT. rho_insitu))
+                    ilevel=jk-1
+                    rho = sigmai( tsn(ji,jj,icona,jp_tem), tsn(ji,jj,icona,jp_sal), fsdept(ji,jj,jk) )
+                    rho_insitu = sigmai( tsn(ji,jj,jk,jp_tem), tsn(ji,jj,jk,jp_sal), fsdept(ji,jj,jk) )
+                    jk = jk+1
+                  END DO
+
+                  !Found the right level so replace water in lower box
+                  dztdif = fse3t(ji,jj,ilevel)-fse3t(ji,jj,icona)
+                  DO jn = 1, jpts ! Loop over tracers
+                    tra(jn) = tsn(ji,jj,ilevel,jn)
+                    trasum(jn) = tsn(ji,jj,icona,jn)*fse3t(ji,jj,icona) + tsn(ji,jj,ilevel,jn)*dztdif
+                    tramix = trasum(jn) / fse3t(ji,jj,ilevel)
+                    tsn(ji,jj,ilevel,jn) = tramix
+                  END DO
+
+
+                  !Shuffle up the intervening boxes
+                  DO jn = 1,jpts
+                    DO jk = ilevel-1,icona,-1
+                      dztdif = fse3t(ji,jj,jk)-fse3t(ji,jj,icona)
+                      trasum(jn) = tra(jn) * fse3t(ji,jj,icona) + tsn(ji,jj,jk,jn)*dztdif
+                      tramix = trasum(jn)/fse3t(ji,jj,jk)
+                      tra(jn) = tsn(ji,jj,jk,jn)
+                      tsn(ji,jj,jk,jn) = tramix
+                    END DO
+                  END DO
+
+                  itop = icona+1
+                ELSE
+                  itop=jpk
+                ENDIF
+              END DO ! End of WHILE loop
+            END DO !End loop over nn_con
+          END IF
+        END DO
+      END DO 
+
+
+      CALL wrk_dealloc(jpkm1, rhoup, rholo)
+      CALL wrk_dealloc(jpts, tra, trasum)
+   END SUBROUTINE zdf_convadj
+
+
+   SUBROUTINE zdf_convadj_init
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE zdf_convadj  ***
+      !! ** Purpose: Initialise zdf_convadj module
+      !!
+      !! ** Method: Read namelist options
+      !!            Read mask from file
+      !!----------------------------------------------------------------------
+      USE iom
+
+      INTEGER               ::   imask
+
+      NAMELIST/namzdf_convadj/ nn_con, cn_convadj_mask, ln_convadj
+
+      REWIND ( numnam_ref )               !* Read Namelist namzdf_tke : Turbulente Kinetic Energy
+      READ   ( numnam_ref, namzdf_convadj )
+
+
+      REWIND ( numnam_cfg )               !* Read Namelist namzdf_tke : Turbulente Kinetic Energy
+      READ   ( numnam_cfg, namzdf_convadj )
+
+      IF(lwp) THEN                    !* Control print
+         WRITE(numout,*)
+         WRITE(numout,*) 'zdf_convadj_init : Convective adjustment scheme - initialisation'
+         WRITE(numout,*) '~~~~~~~~~~~~'
+         WRITE(numout,*) '   Namelist namzdf_convadj : set convadj parameters'
+         WRITE(numout,*) '      Use convective adjustment (cliff) scheme    ln_convadj= ', ln_convadj
+         WRITE(numout,*) '      Number of times to convect                  nn_con    = ', nn_con
+         WRITE(numout,*) '      NetCDF mask file of points to convect       cn_convadj_mask  = ',cn_convadj_mask
+      ENDIF
+
+      !Read in mask
+      IF( ln_convadj) THEN
+         ALLOCATE( convadj_mask(jpi,jpj) )
+         CALL iom_open ( cn_convadj_mask, imask )
+         CALL iom_get  ( imask, jpdom_autoglo, 'convadj_mask', convadj_mask)
+         CALL iom_close( imask )
+      ENDIF
+
+   END SUBROUTINE zdf_convadj_init
+
+   FUNCTION sigmai ( ptem, psal, pref )
+       !! --------------------------------------------------------------------
+       !! ** Purpose :   Compute the  density referenced to pref (ratio rho/rau0) 
+       !!       from potential temperature and
+       !!      salinity fields using an equation of state defined through the
+       !!     namelist parameter neos.
+       !!
+       !! ** Method  :
+       !!       Jackett and McDougall (1994) equation of state.
+       !!         the in situ density is computed directly as a function of
+       !!         potential temperature relative to the surface (the opa t
+       !!         variable), salt and pressure (assuming no pressure variation
+       !!         along geopotential surfaces, i.e. the pressure p in decibars
+       !!         is approximated by the depth in meters.
+       !!              prd(t,s,p) = ( rho(t,s,p) - rau0 ) / rau0
+       !!              rhop(t,s)  = rho(t,s,0)
+       !!         with pressure                      p        decibars
+       !!              potential temperature         t        deg celsius
+       !!              salinity                      s        psu
+       !!              reference volumic mass        rau0     kg/m**3
+       !!              in situ volumic mass          rho      kg/m**3
+       !!              in situ density anomalie      prd      no units
+       !! --------------------------------------------------------------------
+ 
+       !! * Arguments
+       REAL(wp), INTENT(in) :: ptem, psal
+       REAL(wp), INTENT(in) :: pref      !: reference pressure (meters or db)
+       REAL(wp)             :: sigmai
+ 
+       REAL(kind=8),PARAMETER :: dpr4=4.8314d-4,dpd=-2.042967d-2 , dprau0 = 1000;
+
+       !! * local variables
+       INTEGER :: ji,jj 
+       REAL(wp) ::  dlrs
+       REAL(wp) ::  dlt, dls      
+       REAL(wp) :: dla,dla1,dlaw,dlb,dlb1,dlbw,dlc,dle,dlk0,dlkw 
+       REAL(wp) :: dlrhop, dlr1,dlr2,dlr3, dlref
+      
+       dlref = pref
+       sigmai = 0.d0
+       dlrs = SQRT( ABS( psal ) )
+
+! Convert T and S to double precision.
+        dlt=DBLE(ptem)
+        dls=DBLE(psal)
+
+
+! Compute the volumic mass of pure water at atmospheric pressure.
+        dlr1=((((6.536332d-9*dlt-1.120083d-6)&
+                *dlt+1.001685d-4)&
+               *dlt-9.095290d-3)&
+              *dlt+6.793952d-2)&
+             *dlt+999.842594d0
+
+! Compute the seawater volumic mass at atmospheric pressure.
+        dlr2=(((5.3875d-9*dlt-8.2467d-7)&
+               *dlt+7.6438d-5)&
+              *dlt-4.0899d-3)&
+             *dlt+0.824493d0
+
+        dlr3=(-1.6546d-6*dlt+1.0227d-4)&
+             *dlt-5.72466d-3
+
+! Compute the potential volumic mass (referenced to the surface).
+        dlrhop=(dpr4*dls+dlr3*dlrs+dlr2)*dls+dlr1
+
+! Compute the compression terms.
+        dle=(-3.508914d-8*dlt-1.248266d-8)&
+            *dlt-2.595994d-6
+
+        dlbw=(1.296821d-6*dlt-5.782165d-9)&
+             *dlt+1.045941d-4
+
+        dlb=dlbw+dle*dls
+
+        dlc=(-7.267926d-5*dlt+2.598241d-3 )&
+            *dlt+0.1571896d0
+
+        dlaw=((5.939910d-6*dlt+2.512549d-3)&
+              *dlt-0.1028859d0)&
+             *dlt-4.721788d0
+
+        dla=(dpd*dlrs+dlc)*dls+dlaw
+
+        dlb1=(-0.1909078d0*dlt+7.390729d0)&
+             *dlt-55.87545d0
+
+        dla1=((2.326469d-3*dlt+1.553190d0)&
+              *dlt-65.00517d0)&
+             *dlt+1044.077d0
+
+        dlkw=(((-1.361629d-4*dlt-1.852732d-2)&
+               *dlt-30.41638d0)&
+              *dlt+2098.925d0)&
+             *dlt+190925.6d0
+
+        dlk0=(dlb1*dlrs+dla1)*dls+dlkw
+
+! Compute the potential density anomaly.
+        sigmai=dlrhop/(1.0d0-dlref/(dlk0-dlref*(dla-dlref*dlb)))&
+                       -dprau0
+ 
+    END FUNCTION sigmai
+
+END MODULE zdfconvadj
+
Index: branches/UKMO/dev_r5518_convadj/NEMOGCM/NEMO/OPA_SRC/step.F90
===================================================================
--- branches/UKMO/dev_r5518_convadj/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 6836)
+++ branches/UKMO/dev_r5518_convadj/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 6837)
@@ -98,5 +98,5 @@
 
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-      ! Update data, open boundaries, surface boundary condition (including sea-ice)
+      ! Update data, open boundaries, surface boundary condition (including sea-ice), convective adjustment (cliff)
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
       IF( lk_tide    )   CALL sbc_tide( kstp )
@@ -105,4 +105,8 @@
                          CALL bdy_dta ( kstp, time_offset=+1 )   ! update dynamic & tracer data at open boundaries
       ENDIF
+
+      IF( ln_convadj )   CALL zdf_convadj( kstp )     ! Called here so that tsn is updated before rest of physics is called
+
+
                          CALL sbc    ( kstp )         ! Sea Boundary Condition (including sea-ice)
                                                       ! clem: moved here for bdy ice purpose
Index: branches/UKMO/dev_r5518_convadj/NEMOGCM/NEMO/OPA_SRC/step_oce.F90
===================================================================
--- branches/UKMO/dev_r5518_convadj/NEMOGCM/NEMO/OPA_SRC/step_oce.F90	(revision 6836)
+++ branches/UKMO/dev_r5518_convadj/NEMOGCM/NEMO/OPA_SRC/step_oce.F90	(revision 6837)
@@ -84,4 +84,5 @@
    USE zdfric           ! Richardson vertical mixing       (zdf_ric routine)
    USE zdfmxl           ! Mixed-layer depth                (zdf_mxl routine)
+   USE zdfconvadj       ! Convective adjustment            (zdf_convadj routine)
 
    USE zpshde           ! partial step: hor. derivative     (zps_hde routine)
