Index: /branches/UKMO/dev_r5518_25hr_mean_assim_bkg/NEMOGCM/NEMO/OPA_SRC/ASM/asmbkg.F90
===================================================================
--- /branches/UKMO/dev_r5518_25hr_mean_assim_bkg/NEMOGCM/NEMO/OPA_SRC/ASM/asmbkg.F90	(revision 6760)
+++ /branches/UKMO/dev_r5518_25hr_mean_assim_bkg/NEMOGCM/NEMO/OPA_SRC/ASM/asmbkg.F90	(revision 6761)
@@ -120,4 +120,5 @@
 #endif
             CALL iom_rstput( kt, nitbkg_r, inum, 'gcx'    , gcx               )
+            CALL iom_rstput( kt, nitbkg_r, inum, 'avt'    , avt               )
             !
             CALL iom_close( inum )
Index: /branches/UKMO/dev_r5518_25hr_mean_assim_bkg/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90
===================================================================
--- /branches/UKMO/dev_r5518_25hr_mean_assim_bkg/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90	(revision 6760)
+++ /branches/UKMO/dev_r5518_25hr_mean_assim_bkg/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90	(revision 6761)
@@ -40,4 +40,16 @@
 #endif
    USE sbc_oce          ! Surface boundary condition variables.
+   USE zdfmxl, ONLY :  &  
+   &  hmld_tref,       &   
+#if defined key_karaml
+   &  hmld_kara,       &
+   &  ln_kara,         &
+#endif   
+   &  hmld,            & 
+   &  hmlp,            &
+   &  hmlpt
+#if defined key_bdy 
+   USE bdy_oce, ONLY: bdytmask  
+#endif  
 
    IMPLICIT NONE
@@ -88,4 +100,11 @@
    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   seaice_bkginc         ! Increment to the background sea ice conc
 
+   INTEGER :: mld_choice        = 4   !: choice of mld criteria to use for physics assimilation
+                                      !: 1) hmld      - Turbocline/mixing depth                           [W points]
+                                      !: 2) hmlp      - Density criterion (0.01 kg/m^3 change from 10m)   [W points]
+                                      !: 3) hmld_kara - Kara MLD                                          [Interpolated]
+                                      !: 4) hmld_tref - Temperature criterion (0.2 K change from surface) [T points]
+
+
    !! * Substitutions
 #  include "domzgr_substitute.h90"
@@ -119,4 +138,5 @@
       INTEGER :: iitiaustr_date  ! Date YYYYMMDD of IAU interval start time step
       INTEGER :: iitiaufin_date  ! Date YYYYMMDD of IAU interval final time step
+      INTEGER :: isurfstat       ! Local integer for status of reading surft variable
       !
       REAL(wp) :: znorm        ! Normalization factor for IAU weights
@@ -127,5 +147,14 @@
       REAL(wp) :: zdate_inc    ! Time axis in increments file
       !
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: & 
+          &       t_bkginc_2d  ! file for reading in 2D  
+      !                        ! temperature increments 
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: & 
+          &       z_mld     ! Mixed layer depth 
+          
       REAL(wp), POINTER, DIMENSION(:,:) ::   hdiv   ! 2D workspace
+      !
+      LOGICAL :: lk_surft      ! Logical: T => Increments file contains surft variable 
+                               !               so only apply surft increments.
       !!
       NAMELIST/nam_asminc/ ln_bkgwri,                                      &
@@ -133,5 +162,5 @@
          &                 ln_asmdin, ln_asmiau,                           &
          &                 nitbkg, nitdin, nitiaustr, nitiaufin, niaufn,   &
-         &                 ln_salfix, salfixmin, nn_divdmp
+         &                 ln_salfix, salfixmin, nn_divdmp, mld_choice
       !!----------------------------------------------------------------------
 
@@ -139,6 +168,21 @@
       ! Read Namelist nam_asminc : assimilation increment interface
       !-----------------------------------------------------------------------
+
+      ! Set default values
+      ln_bkgwri = .FALSE.
+      ln_trainc = .FALSE.
+      ln_dyninc = .FALSE.
+      ln_sshinc = .FALSE.
       ln_seaiceinc = .FALSE.
+      ln_asmdin = .FALSE.
+      ln_asmiau = .TRUE.
+      ln_salfix = .FALSE.
       ln_temnofreeze = .FALSE.
+      salfixmin = -9999
+      nitbkg    = 0
+      nitdin    = 0      
+      nitiaustr = 1
+      nitiaufin = 150
+      niaufn    = 0
 
       REWIND( numnam_ref )              ! Namelist nam_asminc in reference namelist : Assimilation increment
@@ -171,4 +215,5 @@
          WRITE(numout,*) '      Logical switch for ensuring that the sa > salfixmin      ln_salfix = ', ln_salfix
          WRITE(numout,*) '      Minimum salinity after applying the increments           salfixmin = ', salfixmin
+         WRITE(numout,*) '      Choice of MLD for physics assimilation                  mld_choice = ', mld_choice
       ENDIF
 
@@ -327,20 +372,26 @@
       !--------------------------------------------------------------------
 
-      ALLOCATE( t_bkginc(jpi,jpj,jpk) )
-      ALLOCATE( s_bkginc(jpi,jpj,jpk) )
-      ALLOCATE( u_bkginc(jpi,jpj,jpk) )
-      ALLOCATE( v_bkginc(jpi,jpj,jpk) )
-      ALLOCATE( ssh_bkginc(jpi,jpj)   )
-      ALLOCATE( seaice_bkginc(jpi,jpj))
+      IF ( ln_trainc ) THEN
+         ALLOCATE( t_bkginc(jpi,jpj,jpk) )
+         ALLOCATE( s_bkginc(jpi,jpj,jpk) )
+         t_bkginc(:,:,:) = 0.0
+         s_bkginc(:,:,:) = 0.0
+      ENDIF
+      IF ( ln_dyninc ) THEN 
+         ALLOCATE( u_bkginc(jpi,jpj,jpk) )
+         ALLOCATE( v_bkginc(jpi,jpj,jpk) )
+         u_bkginc(:,:,:) = 0.0
+         v_bkginc(:,:,:) = 0.0
+      ENDIF
+      IF ( ln_sshinc ) THEN
+         ALLOCATE( ssh_bkginc(jpi,jpj)   )
+         ssh_bkginc(:,:) = 0.0
+      ENDIF
+      IF ( ln_seaiceinc ) THEN 
+         ALLOCATE( seaice_bkginc(jpi,jpj))
+         seaice_bkginc(:,:) = 0.0
+      ENDIF
 #if defined key_asminc
       ALLOCATE( ssh_iau(jpi,jpj)      )
-#endif
-      t_bkginc(:,:,:) = 0.0
-      s_bkginc(:,:,:) = 0.0
-      u_bkginc(:,:,:) = 0.0
-      v_bkginc(:,:,:) = 0.0
-      ssh_bkginc(:,:) = 0.0
-      seaice_bkginc(:,:) = 0.0
-#if defined key_asminc
       ssh_iau(:,:)    = 0.0
 #endif
@@ -378,13 +429,72 @@
 
          IF ( ln_trainc ) THEN   
-            CALL iom_get( inum, jpdom_autoglo, 'bckint', t_bkginc, 1 )
-            CALL iom_get( inum, jpdom_autoglo, 'bckins', s_bkginc, 1 )
-            ! Apply the masks
-            t_bkginc(:,:,:) = t_bkginc(:,:,:) * tmask(:,:,:)
-            s_bkginc(:,:,:) = s_bkginc(:,:,:) * tmask(:,:,:)
-            ! Set missing increments to 0.0 rather than 1e+20
-            ! to allow for differences in masks
-            WHERE( ABS( t_bkginc(:,:,:) ) > 1.0e+10 ) t_bkginc(:,:,:) = 0.0
-            WHERE( ABS( s_bkginc(:,:,:) ) > 1.0e+10 ) s_bkginc(:,:,:) = 0.0
+            
+            !Test if the increments file contains the surft variable.
+            isurfstat = iom_varid( inum, 'bckinsurft', ldstop = .FALSE. )
+            IF ( isurfstat == -1 ) THEN
+               lk_surft = .FALSE.
+            ELSE
+               lk_surft = .TRUE.
+               CALL ctl_warn( ' Applying 2D temperature increment to bottom of ML: ', &
+            &                 ' bckinsurft found in increments file.' )
+            ENDIF             
+
+            IF (lk_surft) THEN 
+                
+               ALLOCATE(z_mld(jpi,jpj)) 
+               SELECT CASE(mld_choice) 
+               CASE(1) 
+                  z_mld = hmld 
+               CASE(2) 
+                  z_mld = hmlp 
+               CASE(3) 
+#if defined key_karaml
+                  IF ( ln_kara ) THEN
+                     z_mld = hmld_kara
+                  ELSE
+                     CALL ctl_stop("Kara mixed layer not calculated as ln_kara=.false.")
+                  ENDIF
+#else
+                  CALL ctl_stop("Kara mixed layer not defined in current version of NEMO")  ! JW: Safety feature, should be removed
+                                                                                            ! once the Kara mixed layer is available 
+#endif
+               CASE(4) 
+                  z_mld = hmld_tref 
+               END SELECT       
+                      
+               ALLOCATE( t_bkginc_2d(jpi,jpj) ) 
+               CALL iom_get( inum, jpdom_autoglo, 'bckinsurft', t_bkginc_2d, 1) 
+#if defined key_bdy                
+               DO jk = 1,jpkm1 
+                  WHERE( z_mld(:,:) > fsdepw(:,:,jk) ) 
+                     t_bkginc(:,:,jk) = t_bkginc_2d(:,:) * 0.5 * &
+                     &       ( 1 + cos( (fsdept(:,:,jk)/z_mld(:,:) ) * rpi ) )
+                     
+                     t_bkginc(:,:,jk) = t_bkginc(:,:,jk) * bdytmask(:,:) 
+                  ELSEWHERE 
+                     t_bkginc(:,:,jk) = 0. 
+                  ENDWHERE 
+               ENDDO 
+#else 
+               t_bkginc(:,:,:) = 0. 
+#endif                
+               s_bkginc(:,:,:) = 0. 
+                
+               DEALLOCATE(z_mld, t_bkginc_2d) 
+            
+            ELSE 
+               
+               CALL iom_get( inum, jpdom_autoglo, 'bckint', t_bkginc, 1 )
+               CALL iom_get( inum, jpdom_autoglo, 'bckins', s_bkginc, 1 )
+               ! Apply the masks
+               t_bkginc(:,:,:) = t_bkginc(:,:,:) * tmask(:,:,:)
+               s_bkginc(:,:,:) = s_bkginc(:,:,:) * tmask(:,:,:)
+               ! Set missing increments to 0.0 rather than 1e+20
+               ! to allow for differences in masks
+               WHERE( ABS( t_bkginc(:,:,:) ) > 1.0e+10 ) t_bkginc(:,:,:) = 0.0
+               WHERE( ABS( s_bkginc(:,:,:) ) > 1.0e+10 ) s_bkginc(:,:,:) = 0.0
+         
+            ENDIF
+         
          ENDIF
 
Index: /branches/UKMO/dev_r5518_25hr_mean_assim_bkg/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90
===================================================================
--- /branches/UKMO/dev_r5518_25hr_mean_assim_bkg/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90	(revision 6760)
+++ /branches/UKMO/dev_r5518_25hr_mean_assim_bkg/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90	(revision 6761)
@@ -74,5 +74,5 @@
       INTEGER, INTENT(in) ::   kt                      ! time step
       ! 
-      INTEGER             ::   jk                      ! dummy loop indice
+      INTEGER             ::   ji, jj, jk              ! dummy loop indices
       REAL(wp)            ::   z2dt, z1_rau0           ! local scalars
       !!----------------------------------------------------------------------
@@ -94,4 +94,19 @@
       z2dt = 2._wp * rdt                              ! set time step size (Euler/Leapfrog)
       IF( neuler == 0 .AND. kt == nit000 )   z2dt = rdt
+
+#if defined key_asminc
+      !                                                ! Include the IAU weighted SSH increment
+      IF( lk_asminc .AND. ln_sshinc .AND. ln_asmiau ) THEN
+         CALL ssh_asm_inc( kt )
+#if defined key_vvl
+! Don't directly adjust ssh but change hdivn at all levels instead
+! In trasbc also add in the heat and salt content associated with these changes at each level  
+        DO jk = 1, jpkm1                                 
+                 hdivn(:,:,jk) = hdivn(:,:,jk) - ( ssh_iau(:,:) / ( ht_0(:,:) + 1.0 - ssmask(:,:) ) ) * ( e3t_0(:,:,jk) / fse3t_n(:,:,jk) ) * tmask(:,:,jk) 
+        END DO
+        CALL lbc_lnk( hdivn, 'T', 1. ) ! Not sure that's necessary
+      ENDIF
+#endif
+#endif
 
       !                                           !------------------------------!
@@ -123,12 +138,4 @@
 #endif
 
-#if defined key_asminc
-      !                                                ! Include the IAU weighted SSH increment
-      IF( lk_asminc .AND. ln_sshinc .AND. ln_asmiau ) THEN
-         CALL ssh_asm_inc( kt )
-         ssha(:,:) = ssha(:,:) + z2dt * ssh_iau(:,:)
-      ENDIF
-#endif
-
       !                                           !------------------------------!
       !                                           !           outputs            !
Index: /branches/UKMO/dev_r5518_25hr_mean_assim_bkg/NEMOGCM/NEMO/OPA_SRC/TRA/trasbc.F90
===================================================================
--- /branches/UKMO/dev_r5518_25hr_mean_assim_bkg/NEMOGCM/NEMO/OPA_SRC/TRA/trasbc.F90	(revision 6760)
+++ /branches/UKMO/dev_r5518_25hr_mean_assim_bkg/NEMOGCM/NEMO/OPA_SRC/TRA/trasbc.F90	(revision 6761)
@@ -33,4 +33,7 @@
    USE timing          ! Timing
    USE eosbn2
+#if defined key_asminc   
+   USE asminc          ! Assimilation increment
+#endif
 
    IMPLICIT NONE
@@ -278,4 +281,23 @@
          END DO  
       ENDIF
+
+#if defined key_asminc
+! WARNING: THIS MAY WELL NOT BE REQUIRED - WE DON'T WANT TO CHANGE T&S BUT THIS MAY COMPENSATE ANOTHER TERM...
+! Rate of change in e3t for each level is ssh_iau*e3t_0/ht_0
+! Contribution to tsa should be rate of change in level / per m of ocean? (hence the division by fse3t_n)
+      IF( ln_sshinc ) THEN         ! input of heat and salt due to assimilation
+         DO jj = 2, jpj 
+            DO ji = fs_2, fs_jpim1
+               zdep = ssh_iau(ji,jj) / ( ht_0(ji,jj) + 1.0 - ssmask(ji, jj) )
+               DO jk = 1, jpkm1
+                  tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem)   &
+                                        &            + tsn(ji,jj,jk,jp_tem) * zdep * ( e3t_0(ji,jj,jk) / fse3t_n(ji,jj,jk) )
+                  tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal)   &
+                                        &            + tsn(ji,jj,jk,jp_sal) * zdep * ( e3t_0(ji,jj,jk) / fse3t_n(ji,jj,jk) )
+               END DO
+            END DO  
+         END DO  
+      ENDIF
+#endif
  
       IF( l_trdtra )   THEN                      ! send trends for further diagnostics
Index: /branches/UKMO/dev_r5518_25hr_mean_assim_bkg/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfmxl.F90
===================================================================
--- /branches/UKMO/dev_r5518_25hr_mean_assim_bkg/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfmxl.F90	(revision 6760)
+++ /branches/UKMO/dev_r5518_25hr_mean_assim_bkg/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfmxl.F90	(revision 6761)
@@ -27,4 +27,5 @@
 
    PUBLIC   zdf_mxl       ! called by step.F90
+   PUBLIC   zdf_mxl_tref  ! called by asminc.F90
    PUBLIC   zdf_mxl_alloc ! Used in zdf_tke_init
 
@@ -33,4 +34,5 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hmlp    !: mixed layer depth  (rho=rho0+zdcrit) [m]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hmlpt   !: mixed layer depth at t-points        [m]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hmld_tref  !: mixed layer depth at t-points - temperature criterion [m]
 
    REAL(wp), PUBLIC ::   rho_c = 0.01_wp    !: density criterion for mixed layer depth
@@ -52,5 +54,6 @@
       zdf_mxl_alloc = 0      ! set to zero if no array to be allocated
       IF( .NOT. ALLOCATED( nmln ) ) THEN
-         ALLOCATE( nmln(jpi,jpj), hmld(jpi,jpj), hmlp(jpi,jpj), hmlpt(jpi,jpj), STAT= zdf_mxl_alloc )
+         ALLOCATE( nmln(jpi,jpj), hmld(jpi,jpj), hmlp(jpi,jpj), hmlpt(jpi,jpj), &
+         &                           hmld_tref(jpi,jpj), STAT= zdf_mxl_alloc )
          !
          IF( lk_mpp             )   CALL mpp_sum ( zdf_mxl_alloc )
@@ -146,4 +149,48 @@
    END SUBROUTINE zdf_mxl
 
+
+   SUBROUTINE zdf_mxl_tref()
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE zdf_mxl_tref  ***
+      !!                   
+      !! ** Purpose :   Compute the mixed layer depth with temperature criteria.
+      !!
+      !! ** Method  :   The temperature-defined mixed layer depth is required
+      !!                   when assimilating SST in a 2D analysis. 
+      !!
+      !! ** Action  :   hmld_tref
+      !!----------------------------------------------------------------------
+      !
+      INTEGER  ::   ji, jj, jk   ! dummy loop indices
+      REAL(wp) ::   t_ref               ! Reference temperature  
+      REAL(wp) ::   temp_c = 0.2        ! temperature criterion for mixed layer depth  
+      !!----------------------------------------------------------------------
+      !
+      ! Initialise array
+      IF( zdf_mxl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'zdf_mxl_tref : unable to allocate arrays' )
+      
+      !For the AMM model assimiation uses a temperature based mixed layer depth  
+      !This is defined here  
+      DO jj = 1, jpj  
+         DO ji = 1, jpi  
+           hmld_tref(ji,jj)=fsdept(ji,jj,1  )   
+           IF(ssmask(ji,jj) > 0.)THEN  
+             t_ref=tsn(ji,jj,1,jp_tem) 
+             DO jk=2,jpk  
+               IF(ssmask(ji,jj)==0.)THEN  
+                  hmld_tref(ji,jj)=fsdept(ji,jj,jk )  
+                  EXIT  
+               ELSEIF( ABS(tsn(ji,jj,jk,jp_tem)-t_ref) < temp_c)THEN  
+                  hmld_tref(ji,jj)=fsdept(ji,jj,jk )  
+               ELSE  
+                  EXIT  
+               ENDIF  
+             ENDDO  
+           ENDIF  
+         ENDDO  
+      ENDDO
+
+   END SUBROUTINE zdf_mxl_tref
+
    !!======================================================================
 END MODULE zdfmxl
Index: /branches/UKMO/dev_r5518_25hr_mean_assim_bkg/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90
===================================================================
--- /branches/UKMO/dev_r5518_25hr_mean_assim_bkg/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90	(revision 6760)
+++ /branches/UKMO/dev_r5518_25hr_mean_assim_bkg/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90	(revision 6761)
@@ -473,5 +473,11 @@
 
       !                                     ! Assimilation increments
-      IF( lk_asminc     )   CALL asm_inc_init   ! Initialize assimilation increments
+      IF( lk_asminc ) THEN 
+#if defined key_shelf 
+         CALL  zdf_mxl_tref()     ! Initialization of hmld_tref
+#endif 
+         CALL asm_inc_init     ! Initialize assimilation increments 
+      ENDIF
+           
       IF(lwp) WRITE(numout,*) 'Euler time step switch is ', neuler
       !
