Index: branches/UKMO/r6232_collate_bgc_diagnostics/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90
===================================================================
--- branches/UKMO/r6232_collate_bgc_diagnostics/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90	(revision 11132)
+++ branches/UKMO/r6232_collate_bgc_diagnostics/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90	(revision 11134)
@@ -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
@@ -57,4 +69,5 @@
 #endif
    LOGICAL, PUBLIC :: ln_bkgwri = .FALSE.      !: No output of the background state fields
+   LOGICAL, PUBLIC :: ln_avgbkg = .FALSE.      !: No output of the mean background state fields
    LOGICAL, PUBLIC :: ln_asmiau = .FALSE.      !: No applying forcing with an assimilation increment
    LOGICAL, PUBLIC :: ln_asmdin = .FALSE.      !: No direct initialization
@@ -80,4 +93,5 @@
    INTEGER , PUBLIC ::   nitiaustr   !: Time step of the start of the IAU interval 
    INTEGER , PUBLIC ::   nitiaufin   !: Time step of the end of the IAU interval
+   INTEGER , PUBLIC ::   nitavgbkg   !: Number of timesteps to average assim bkg [0,nitavgbkg]
    ! 
    INTEGER , PUBLIC ::   niaufn      !: Type of IAU weighing function: = 0   Constant weighting
@@ -87,4 +101,11 @@
    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   ssh_bkg, ssh_bkginc   ! Background sea surface height and its increment
    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
@@ -119,4 +140,6 @@
       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
+      INTEGER :: iitavgbkg_date  ! Date YYYYMMDD of end of assim bkg averaging period
       !
       REAL(wp) :: znorm        ! Normalization factor for IAU weights
@@ -127,11 +150,20 @@
       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
-      !!
-      NAMELIST/nam_asminc/ ln_bkgwri,                                      &
+      !
+      LOGICAL :: lk_surft      ! Logical: T => Increments file contains surft variable 
+                               !               so only apply surft increments.
+      !!
+      NAMELIST/nam_asminc/ ln_bkgwri, ln_avgbkg,                           &
          &                 ln_trainc, ln_dyninc, ln_sshinc,                &
          &                 ln_asmdin, ln_asmiau,                           &
          &                 nitbkg, nitdin, nitiaustr, nitiaufin, niaufn,   &
-         &                 ln_salfix, salfixmin, nn_divdmp
+         &                 ln_salfix, salfixmin, nn_divdmp, nitavgbkg, mld_choice
       !!----------------------------------------------------------------------
 
@@ -139,6 +171,23 @@
       ! Read Namelist nam_asminc : assimilation increment interface
       !-----------------------------------------------------------------------
+
+      ! Set default values
+      ln_bkgwri = .FALSE.
+      ln_avgbkg = .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
+      nitavgbkg = 1
 
       REWIND( numnam_ref )              ! Namelist nam_asminc in reference namelist : Assimilation increment
@@ -158,4 +207,5 @@
          WRITE(numout,*) '   Namelist namasm : set assimilation increment parameters'
          WRITE(numout,*) '      Logical switch for writing out background state          ln_bkgwri = ', ln_bkgwri
+         WRITE(numout,*) '      Logical switch for writing mean background state         ln_avgbkg = ', ln_avgbkg
          WRITE(numout,*) '      Logical switch for applying tracer increments            ln_trainc = ', ln_trainc
          WRITE(numout,*) '      Logical switch for applying velocity increments          ln_dyninc = ', ln_dyninc
@@ -168,7 +218,9 @@
          WRITE(numout,*) '      Timestep of start of IAU interval in [0,nitend-nit000-1] nitiaustr = ', nitiaustr
          WRITE(numout,*) '      Timestep of end of IAU interval in [0,nitend-nit000-1]   nitiaufin = ', nitiaufin
+         WRITE(numout,*) '      Number of timesteps to average assim bkg [0,nitavgbkg]   nitavgbkg = ', nitavgbkg
          WRITE(numout,*) '      Type of IAU weighting function                           niaufn    = ', niaufn
          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
 
@@ -177,4 +229,5 @@
       nitiaustr_r = nitiaustr + nit000 - 1  ! Start of IAU interval referenced to nit000
       nitiaufin_r = nitiaufin + nit000 - 1  ! End of IAU interval referenced to nit000
+      nitavgbkg_r = nitavgbkg + nit000 - 1  ! Averaging period referenced to nit000
 
       iiauper = nitiaufin_r - nitiaustr_r + 1  ! IAU interval length
@@ -186,4 +239,5 @@
       CALL calc_date( nit000, nitiaustr_r, ndate0, iitiaustr_date )     ! IAU start time referenced to ndate0
       CALL calc_date( nit000, nitiaufin_r, ndate0, iitiaufin_date )     ! IAU end time referenced to ndate0
+      CALL calc_date( nit000, nitavgbkg_r, ndate0, iitavgbkg_date )     ! End of assim bkg averaging period referenced to ndate0
       !
       IF(lwp) THEN
@@ -197,4 +251,5 @@
          WRITE(numout,*) '       nitiaustr_r = ', nitiaustr_r
          WRITE(numout,*) '       nitiaufin_r = ', nitiaufin_r
+         WRITE(numout,*) '       nitavgbkg_r = ', nitavgbkg_r
          WRITE(numout,*)
          WRITE(numout,*) '   Dates referenced to current cycle:'
@@ -206,4 +261,5 @@
          WRITE(numout,*) '       iitiaustr_date = ', iitiaustr_date
          WRITE(numout,*) '       iitiaufin_date = ', iitiaufin_date
+         WRITE(numout,*) '       iitavgbkg_date = ', iitavgbkg_date
       ENDIF
 
@@ -248,4 +304,9 @@
          & CALL ctl_stop( ' nitdin :',  &
          &                ' Background time step for Direct Initialization is outside', &
+         &                ' the cycle interval')
+
+      IF ( nitavgbkg_r > nitend ) &
+         & CALL ctl_stop( ' nitavgbkg_r :',  &
+         &                ' Assim bkg averaging period is outside', &
          &                ' the cycle interval')
 
@@ -327,20 +388,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 +445,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
 
@@ -892,4 +1018,8 @@
             ENDIF
 
+#if defined key_asminc
+         ELSE
+            ssh_iau(:,:) = 0._wp
+#endif
          ENDIF
 
