Index: NEMO/branches/UKMO/NEMO_4.0.4_CO9_shelf_climate/src/OCE/ZDF/zdfmxl.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.4_CO9_shelf_climate/src/OCE/ZDF/zdfmxl.F90	(revision 15382)
+++ NEMO/branches/UKMO/NEMO_4.0.4_CO9_shelf_climate/src/OCE/ZDF/zdfmxl.F90	(revision 15383)
@@ -2,10 +2,10 @@
    !!======================================================================
    !!                       ***  MODULE  zdfmxl  ***
-   !! Ocean physics: mixed layer depth 
+   !! Ocean physics: mixed layer depth
    !!======================================================================
    !! History :  1.0  ! 2003-08  (G. Madec)  original code
    !!            3.2  ! 2009-07  (S. Masson, G. Madec)  IOM + merge of DO-loop
-   !!            3.7  ! 2012-03  (G. Madec)  make public the density criteria for trdmxl 
-   !!             -   ! 2014-02  (F. Roquet)  mixed layer depth calculated using N2 instead of rhop 
+   !!            3.7  ! 2012-03  (G. Madec)  make public the density criteria for trdmxl
+   !!             -   ! 2014-02  (F. Roquet)  mixed layer depth calculated using N2 instead of rhop
    !!----------------------------------------------------------------------
    !!   zdf_mxl      : Compute the turbocline and mixed layer depths.
@@ -44,16 +44,16 @@
 
 
- 
-   ! Namelist variables for  namzdf_karaml 
+
+   ! Namelist variables for  namzdf_karaml
    LOGICAL   :: ln_kara            ! Logical switch for calculating kara mixed
                                      ! layer
    LOGICAL   :: ln_kara_write25h   ! Logical switch for 25 hour outputs
-   INTEGER   :: jpmld_type         ! mixed layer type             
-   REAL(wp)  :: ppz_ref            ! depth of initial T_ref 
-   REAL(wp)  :: ppdT_crit          ! Critical temp diff  
-   REAL(wp)  :: ppiso_frac         ! Fraction of ppdT_crit used 
-   
+   INTEGER   :: jpmld_type         ! mixed layer type
+   REAL(wp)  :: ppz_ref            ! depth of initial T_ref
+   REAL(wp)  :: ppdT_crit          ! Critical temp diff
+   REAL(wp)  :: ppiso_frac         ! Fraction of ppdT_crit used
+
    !Used for 25h mean
-   LOGICAL, PRIVATE :: kara_25h_init = .TRUE.   !Logical used to initalise 25h 
+   LOGICAL, PRIVATE :: kara_25h_init = .TRUE.   !Logical used to initalise 25h
                                                 !outputs. Necissary, because we need to
                                                 !initalise the kara_25h on the zeroth
@@ -69,8 +69,8 @@
 
    TYPE, PUBLIC :: MXL_ZINT   !: Structure for MLD defs
-      INTEGER   :: mld_type   ! mixed layer type     
+      INTEGER   :: mld_type   ! mixed layer type
       REAL(wp)  :: zref       ! depth of initial T_ref
       REAL(wp)  :: dT_crit    ! Critical temp diff
-      REAL(wp)  :: iso_frac   ! Fraction of rn_dT_crit 
+      REAL(wp)  :: iso_frac   ! Fraction of rn_dT_crit
    END TYPE MXL_ZINT
 
@@ -89,5 +89,5 @@
       IF( .NOT. ALLOCATED( nmln ) ) THEN
          ALLOCATE( nmln(jpi,jpj), hmld(jpi,jpj), hmlp(jpi,jpj), hmlpt(jpi,jpj), hmld_zint(jpi,jpj),     &
- 	&          htc_mld(jpi,jpj), ll_found(jpi,jpj), ll_belowml(jpi,jpj,jpk), STAT= zdf_mxl_alloc )         
+ 	&          htc_mld(jpi,jpj), ll_found(jpi,jpj), ll_belowml(jpi,jpj,jpk), STAT= zdf_mxl_alloc )
          !
          CALL mpp_sum ( 'zdfmxl', zdf_mxl_alloc )
@@ -101,9 +101,9 @@
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE zdfmxl  ***
-      !!                   
+      !!
       !! ** Purpose :   Compute the turbocline depth and the mixed layer depth
       !!              with density criteria.
       !!
-      !! ** Method  :   The mixed layer depth is the shallowest W depth with 
+      !! ** Method  :   The mixed layer depth is the shallowest W depth with
       !!      the density of the corresponding T point (just bellow) bellow a
       !!      given value defined locally as rho(10m) + rho_c
@@ -136,5 +136,5 @@
       zN2_c = grav * rho_c * r1_rau0   ! convert density criteria into N^2 criteria
       DO jk = nlb10, jpkm1
-         DO jj = 1, jpj                ! Mixed layer level: w-level 
+         DO jj = 1, jpj                ! Mixed layer level: w-level
             DO ji = 1, jpi
                ikt = mbkt(ji,jj)
@@ -147,8 +147,8 @@
       ! w-level of the turbocline and mixing layer (iom_use)
       imld(:,:) = mbkt(:,:) + 1        ! Initialization to the number of w ocean point
-      DO jk = jpkm1, nlb10, -1         ! from the bottom to nlb10 
+      DO jk = jpkm1, nlb10, -1         ! from the bottom to nlb10
          DO jj = 1, jpj
             DO ji = 1, jpi
-               IF( avt (ji,jj,jk) < avt_c * wmask(ji,jj,jk) )   imld(ji,jj) = jk      ! Turbocline 
+               IF( avt (ji,jj,jk) < avt_c * wmask(ji,jj,jk) )   imld(ji,jj) = jk      ! Turbocline
             END DO
          END DO
@@ -163,5 +163,5 @@
             iiki = imld(ji,jj)
             iikn = nmln(ji,jj)
-            hmld (ji,jj) = gdepw_n(ji,jj,iiki  ) * ssmask(ji,jj)    ! Turbocline depth 
+            hmld (ji,jj) = gdepw_n(ji,jj,iiki  ) * ssmask(ji,jj)    ! Turbocline depth
             hmlp (ji,jj) = gdepw_n(ji,jj,iikn  ) * ssmask(ji,jj)    ! Mixed layer depth
             hmlpt(ji,jj) = gdept_n(ji,jj,iikn-1) * ssmask(ji,jj)    ! depth of the last T-point inside the mixed layer
@@ -189,30 +189,30 @@
    END SUBROUTINE zdf_mxl
 
-   SUBROUTINE zdf_mxl_zint_mld( sf ) 
-      !!---------------------------------------------------------------------------------- 
-      !!                    ***  ROUTINE zdf_mxl_zint_mld  *** 
-      !                                                                        
-      !   Calculate vertically-interpolated mixed layer depth diagnostic. 
-      !            
+   SUBROUTINE zdf_mxl_zint_mld( sf )
+      !!----------------------------------------------------------------------------------
+      !!                    ***  ROUTINE zdf_mxl_zint_mld  ***
+      !
+      !   Calculate vertically-interpolated mixed layer depth diagnostic.
+      !
       !   This routine can calculate the mixed layer depth diagnostic suggested by
       !   Kara et al, 2000, JGR, 105, 16803, but is more general and can calculate
       !   vertically-interpolated mixed-layer depth diagnostics with other parameter
-      !   settings set in the namzdf_mldzint namelist.  
-      ! 
-      !   If mld_type=1 the mixed layer depth is calculated as the depth at which the  
-      !   density has increased by an amount equivalent to a temperature difference of  
-      !   0.8C at the surface. 
-      ! 
-      !   For other values of mld_type the mixed layer is calculated as the depth at  
-      !   which the temperature differs by 0.8C from the surface temperature.  
-      !                                                                        
-      !   David Acreman, Daley Calvert                                      
-      ! 
-      !!----------------------------------------------------------------------------------- 
+      !   settings set in the namzdf_mldzint namelist.
+      !
+      !   If mld_type=1 the mixed layer depth is calculated as the depth at which the
+      !   density has increased by an amount equivalent to a temperature difference of
+      !   0.8C at the surface.
+      !
+      !   For other values of mld_type the mixed layer is calculated as the depth at
+      !   which the temperature differs by 0.8C from the surface temperature.
+      !
+      !   David Acreman, Daley Calvert
+      !
+      !!-----------------------------------------------------------------------------------
 
       TYPE(MXL_ZINT), INTENT(in)  :: sf
 
       ! Diagnostic criteria
-      INTEGER   :: nn_mld_type   ! mixed layer type     
+      INTEGER   :: nn_mld_type   ! mixed layer type
       REAL(wp)  :: rn_zref       ! depth of initial T_ref
       REAL(wp)  :: rn_dT_crit    ! Critical temp diff
@@ -221,21 +221,21 @@
       ! Local variables
       REAL(wp), PARAMETER :: zepsilon = 1.e-30          ! local small value
-      INTEGER, DIMENSION(jpi,jpj) :: ikmt          ! number of active tracer levels 
-      INTEGER, DIMENSION(jpi,jpj) :: ik_ref        ! index of reference level 
-      INTEGER, DIMENSION(jpi,jpj) :: ik_iso        ! index of last uniform temp level 
-      REAL, DIMENSION(jpi,jpj,jpk)  :: zT            ! Temperature or density 
-      REAL, DIMENSION(jpi,jpj)    :: ppzdep        ! depth for use in calculating d(rho) 
-      REAL, DIMENSION(jpi,jpj)    :: zT_ref        ! reference temperature 
-      REAL    :: zT_b                                   ! base temperature 
-      REAL, DIMENSION(jpi,jpj,jpk)  :: zdTdz         ! gradient of zT 
-      REAL, DIMENSION(jpi,jpj,jpk)  :: zmoddT        ! Absolute temperature difference 
-      REAL    :: zdz                                    ! depth difference 
-      REAL    :: zdT                                    ! temperature difference 
-      REAL, DIMENSION(jpi,jpj)    :: zdelta_T      ! difference critereon 
-      REAL, DIMENSION(jpi,jpj)    :: zRHO1, zRHO2  ! Densities 
-      INTEGER :: ji, jj, jk                             ! loop counter 
-
-      !!------------------------------------------------------------------------------------- 
-      !  
+      INTEGER, DIMENSION(jpi,jpj) :: ikmt          ! number of active tracer levels
+      INTEGER, DIMENSION(jpi,jpj) :: ik_ref        ! index of reference level
+      INTEGER, DIMENSION(jpi,jpj) :: ik_iso        ! index of last uniform temp level
+      REAL, DIMENSION(jpi,jpj,jpk)  :: zT            ! Temperature or density
+      REAL, DIMENSION(jpi,jpj)    :: ppzdep        ! depth for use in calculating d(rho)
+      REAL, DIMENSION(jpi,jpj)    :: zT_ref        ! reference temperature
+      REAL    :: zT_b                                   ! base temperature
+      REAL, DIMENSION(jpi,jpj,jpk)  :: zdTdz         ! gradient of zT
+      REAL, DIMENSION(jpi,jpj,jpk)  :: zmoddT        ! Absolute temperature difference
+      REAL    :: zdz                                    ! depth difference
+      REAL    :: zdT                                    ! temperature difference
+      REAL, DIMENSION(jpi,jpj)    :: zdelta_T      ! difference critereon
+      REAL, DIMENSION(jpi,jpj)    :: zRHO1, zRHO2  ! Densities
+      INTEGER :: ji, jj, jk                             ! loop counter
+
+      !!-------------------------------------------------------------------------------------
+      !
       ! Unpack structure
       nn_mld_type = sf%mld_type
@@ -244,51 +244,51 @@
       rn_iso_frac = sf%iso_frac
 
-      ! Set the mixed layer depth criterion at each grid point 
+      ! Set the mixed layer depth criterion at each grid point
       IF( nn_mld_type == 0 ) THEN
          zdelta_T(:,:) = rn_dT_crit
          zT(:,:,:) = rhop(:,:,:)
       ELSE IF( nn_mld_type == 1 ) THEN
-         ppzdep(:,:)=0.0 
-         call eos ( tsn(:,:,1,:), ppzdep(:,:), zRHO1(:,:) ) 
-! Use zT temporarily as a copy of tsn with rn_dT_crit added to SST 
-! [assumes number of tracers less than number of vertical levels] 
-         zT(:,:,1:jpts)=tsn(:,:,1,1:jpts) 
-         zT(:,:,jp_tem)=zT(:,:,1)+rn_dT_crit 
-         CALL eos( zT(:,:,1:jpts), ppzdep(:,:), zRHO2(:,:) ) 
-         zdelta_T(:,:) = abs( zRHO1(:,:) - zRHO2(:,:) ) * rau0 
-         ! RHO from eos (2d version) doesn't calculate north or east halo: 
-         CALL lbc_lnk( 'zdfmxl', zdelta_T, 'T', 1. ) 
-         zT(:,:,:) = rhop(:,:,:) 
-      ELSE 
-         zdelta_T(:,:) = rn_dT_crit                      
-         zT(:,:,:) = tsn(:,:,:,jp_tem)                           
-      END IF 
-
-      ! Calculate the gradient of zT and absolute difference for use later 
-      DO jk = 1 ,jpk-2 
-         zdTdz(:,:,jk)  =    ( zT(:,:,jk+1) - zT(:,:,jk) ) / e3w_n(:,:,jk+1) 
-         zmoddT(:,:,jk) = abs( zT(:,:,jk+1) - zT(:,:,jk) ) 
-      END DO 
-
-      ! Find density/temperature at the reference level (Kara et al use 10m).          
-      ! ik_ref is the index of the box centre immediately above or at the reference level 
-      ! Find rn_zref in the array of model level depths and find the ref    
-      ! density/temperature by linear interpolation.                                   
-      DO jk = jpkm1, 2, -1 
-         WHERE ( gdept_n(:,:,jk) > rn_zref ) 
-           ik_ref(:,:) = jk - 1 
-           zT_ref(:,:) = zT(:,:,jk-1) + zdTdz(:,:,jk-1) * ( rn_zref - gdept_n(:,:,jk-1) ) 
-         END WHERE 
-      END DO 
-
-      ! If the first grid box centre is below the reference level then use the 
-      ! top model level to get zT_ref 
-      WHERE ( gdept_n(:,:,1) > rn_zref )  
-         zT_ref = zT(:,:,1) 
-         ik_ref = 1 
-      END WHERE 
-
-      ! The number of active tracer levels is 1 less than the number of active w levels 
-      ikmt(:,:) = mbkt(:,:) - 1 
+         ppzdep(:,:)=0.0
+         call eos ( tsn(:,:,1,:), ppzdep(:,:), zRHO1(:,:) )
+! Use zT temporarily as a copy of tsn with rn_dT_crit added to SST
+! [assumes number of tracers less than number of vertical levels]
+         zT(:,:,1:jpts)=tsn(:,:,1,1:jpts)
+         zT(:,:,jp_tem)=zT(:,:,1)+rn_dT_crit
+         CALL eos( zT(:,:,1:jpts), ppzdep(:,:), zRHO2(:,:) )
+         zdelta_T(:,:) = abs( zRHO1(:,:) - zRHO2(:,:) ) * rau0
+         ! RHO from eos (2d version) doesn't calculate north or east halo:
+         CALL lbc_lnk( 'zdfmxl', zdelta_T, 'T', 1. )
+         zT(:,:,:) = rhop(:,:,:)
+      ELSE
+         zdelta_T(:,:) = rn_dT_crit
+         zT(:,:,:) = tsn(:,:,:,jp_tem)
+      END IF
+
+      ! Calculate the gradient of zT and absolute difference for use later
+      DO jk = 1 ,jpk-2
+         zdTdz(:,:,jk)  =    ( zT(:,:,jk+1) - zT(:,:,jk) ) / e3w_n(:,:,jk+1)
+         zmoddT(:,:,jk) = abs( zT(:,:,jk+1) - zT(:,:,jk) )
+      END DO
+
+      ! Find density/temperature at the reference level (Kara et al use 10m).
+      ! ik_ref is the index of the box centre immediately above or at the reference level
+      ! Find rn_zref in the array of model level depths and find the ref
+      ! density/temperature by linear interpolation.
+      DO jk = jpkm1, 2, -1
+         WHERE ( gdept_n(:,:,jk) > rn_zref )
+           ik_ref(:,:) = jk - 1
+           zT_ref(:,:) = zT(:,:,jk-1) + zdTdz(:,:,jk-1) * ( rn_zref - gdept_n(:,:,jk-1) )
+         END WHERE
+      END DO
+
+      ! If the first grid box centre is below the reference level then use the
+      ! top model level to get zT_ref
+      WHERE ( gdept_n(:,:,1) > rn_zref )
+         zT_ref = zT(:,:,1)
+         ik_ref = 1
+      END WHERE
+
+      ! The number of active tracer levels is 1 less than the number of active w levels
+      ikmt(:,:) = mbkt(:,:) - 1
 
       ! Initialize / reset
@@ -296,71 +296,71 @@
 
       IF ( rn_iso_frac - zepsilon > 0. ) THEN
-         ! Search for a uniform density/temperature region where adjacent levels          
-         ! differ by less than rn_iso_frac * deltaT.                                      
-         ! ik_iso is the index of the last level in the uniform layer  
-         ! ll_found indicates whether the mixed layer depth can be found by interpolation 
-         ik_iso(:,:)   = ik_ref(:,:) 
-         DO jj = 1, nlcj 
-            DO ji = 1, nlci 
-!CDIR NOVECTOR 
-               DO jk = ik_ref(ji,jj), ikmt(ji,jj)-1 
-                  IF ( zmoddT(ji,jj,jk) > ( rn_iso_frac * zdelta_T(ji,jj) ) ) THEN 
-                     ik_iso(ji,jj)   = jk 
-                     ll_found(ji,jj) = ( zmoddT(ji,jj,jk) > zdelta_T(ji,jj) ) 
-                     EXIT 
-                  END IF 
-               END DO 
-            END DO 
-         END DO 
-
-         ! Use linear interpolation to find depth of mixed layer base where possible 
-         hmld_zint(:,:) = rn_zref 
-         DO jj = 1, jpj 
-            DO ji = 1, jpi 
-               IF (ll_found(ji,jj) .and. tmask(ji,jj,1) == 1.0) THEN 
-                  zdz =  abs( zdelta_T(ji,jj) / zdTdz(ji,jj,ik_iso(ji,jj)) ) 
-                  hmld_zint(ji,jj) = gdept_n(ji,jj,ik_iso(ji,jj)) + zdz 
-               END IF 
-            END DO 
-         END DO 
+         ! Search for a uniform density/temperature region where adjacent levels
+         ! differ by less than rn_iso_frac * deltaT.
+         ! ik_iso is the index of the last level in the uniform layer
+         ! ll_found indicates whether the mixed layer depth can be found by interpolation
+         ik_iso(:,:)   = ik_ref(:,:)
+         DO jj = 1, nlcj
+            DO ji = 1, nlci
+!CDIR NOVECTOR
+               DO jk = ik_ref(ji,jj), ikmt(ji,jj)-1
+                  IF ( zmoddT(ji,jj,jk) > ( rn_iso_frac * zdelta_T(ji,jj) ) ) THEN
+                     ik_iso(ji,jj)   = jk
+                     ll_found(ji,jj) = ( zmoddT(ji,jj,jk) > zdelta_T(ji,jj) )
+                     EXIT
+                  END IF
+               END DO
+            END DO
+         END DO
+
+         ! Use linear interpolation to find depth of mixed layer base where possible
+         hmld_zint(:,:) = rn_zref
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               IF (ll_found(ji,jj) .and. tmask(ji,jj,1) == 1.0) THEN
+                  zdz =  abs( zdelta_T(ji,jj) / zdTdz(ji,jj,ik_iso(ji,jj)) )
+                  hmld_zint(ji,jj) = gdept_n(ji,jj,ik_iso(ji,jj)) + zdz
+               END IF
+            END DO
+         END DO
       END IF
 
-      ! If ll_found = .false. then calculate MLD using difference of zdelta_T    
-      ! from the reference density/temperature 
- 
-! Prevent this section from working on land points 
-      WHERE ( tmask(:,:,1) /= 1.0 ) 
-         ll_found = .true. 
-      END WHERE 
- 
-      DO jk=1, jpk 
-         ll_belowml(:,:,jk) = abs( zT(:,:,jk) - zT_ref(:,:) ) >= zdelta_T(:,:)  
-      END DO 
- 
-! Set default value where interpolation cannot be used (ll_found=false)  
-      DO jj = 1, jpj 
-         DO ji = 1, jpi 
-            IF ( .not. ll_found(ji,jj) )  hmld_zint(ji,jj) = gdept_n(ji,jj,ikmt(ji,jj)) 
-         END DO 
-      END DO 
-
-      DO jj = 1, jpj 
-         DO ji = 1, jpi 
-!CDIR NOVECTOR 
-            DO jk = ik_ref(ji,jj)+1, ikmt(ji,jj) 
-               IF ( ll_found(ji,jj) ) EXIT 
-               IF ( ll_belowml(ji,jj,jk) ) THEN                
-                  zT_b = zT_ref(ji,jj) + zdelta_T(ji,jj) * SIGN(1.0, zdTdz(ji,jj,jk-1) ) 
-                  zdT  = zT_b - zT(ji,jj,jk-1)                                      
-                  zdz  = zdT / zdTdz(ji,jj,jk-1)                                       
-                  hmld_zint(ji,jj) = gdept_n(ji,jj,jk-1) + zdz 
-                  EXIT                                                   
-               END IF 
-            END DO 
-         END DO 
-      END DO 
-
-      hmld_zint(:,:) = hmld_zint(:,:)*tmask(:,:,1) 
-      !  
+      ! If ll_found = .false. then calculate MLD using difference of zdelta_T
+      ! from the reference density/temperature
+
+! Prevent this section from working on land points
+      WHERE ( tmask(:,:,1) /= 1.0 )
+         ll_found = .true.
+      END WHERE
+
+      DO jk=1, jpk
+         ll_belowml(:,:,jk) = abs( zT(:,:,jk) - zT_ref(:,:) ) >= zdelta_T(:,:)
+      END DO
+
+! Set default value where interpolation cannot be used (ll_found=false)
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            IF ( .not. ll_found(ji,jj) )  hmld_zint(ji,jj) = gdept_n(ji,jj,ikmt(ji,jj))
+         END DO
+      END DO
+
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+!CDIR NOVECTOR
+            DO jk = ik_ref(ji,jj)+1, ikmt(ji,jj)
+               IF ( ll_found(ji,jj) ) EXIT
+               IF ( ll_belowml(ji,jj,jk) ) THEN
+                  zT_b = zT_ref(ji,jj) + zdelta_T(ji,jj) * SIGN(1.0, zdTdz(ji,jj,jk-1) )
+                  zdT  = zT_b - zT(ji,jj,jk-1)
+                  zdz  = zdT / zdTdz(ji,jj,jk-1)
+                  hmld_zint(ji,jj) = gdept_n(ji,jj,jk-1) + zdz
+                  EXIT
+               END IF
+            END DO
+         END DO
+      END DO
+
+      hmld_zint(:,:) = hmld_zint(:,:)*tmask(:,:,1)
+      !
    END SUBROUTINE zdf_mxl_zint_mld
 
@@ -368,8 +368,8 @@
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE zdf_mxl_zint_htc  ***
-      !! 
-      !! ** Purpose :   
       !!
-      !! ** Method  :   
+      !! ** Purpose :
+      !!
+      !! ** Method  :
       !!----------------------------------------------------------------------
 
@@ -396,7 +396,7 @@
       zthick_0(:,:) = 0._wp
 
-      DO jk = 1, jpkm1  
+      DO jk = 1, jpkm1
          DO jj = 1, jpj
-            DO ji = 1, jpi                    
+            DO ji = 1, jpi
                zthick_0(ji,jj) = zthick_0(ji,jj) + e3t_n(ji,jj,jk)
                IF( zthick_0(ji,jj) < hmld_zint(ji,jj) )   ilevel(ji,jj) = jk
@@ -408,6 +408,6 @@
 
       ! Surface boundary condition
-      IF( ln_linssh ) THEN  ;   zthick(:,:) = sshn(:,:)   ;   htc_mld(:,:) = tsn(:,:,1,jp_tem) * sshn(:,:) * tmask(:,:,1)   
-      ELSE                  ;   zthick(:,:) = 0._wp       ;   htc_mld(:,:) = 0._wp                                   
+      IF( ln_linssh ) THEN  ;   zthick(:,:) = sshn(:,:)   ;   htc_mld(:,:) = tsn(:,:,1,jp_tem) * sshn(:,:) * tmask(:,:,1)
+      ELSE                  ;   zthick(:,:) = 0._wp       ;   htc_mld(:,:) = 0._wp
       ENDIF
 
@@ -431,5 +431,5 @@
       DO jj = 1, jpj
          DO ji = 1, jpi
-            htc_mld(ji,jj) = htc_mld(ji,jj) + tsn(ji,jj,ilevel(ji,jj)+1,jp_tem)  & 
+            htc_mld(ji,jj) = htc_mld(ji,jj) + tsn(ji,jj,ilevel(ji,jj)+1,jp_tem)  &
       &                      * MIN( e3t_n(ji,jj,ilevel(ji,jj)+1), zthick(ji,jj) ) * tmask(ji,jj,ilevel(ji,jj)+1)
          END DO
@@ -447,8 +447,8 @@
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE zdf_mxl_zint  ***
-      !! 
-      !! ** Purpose :   
       !!
-      !! ** Method  :   
+      !! ** Purpose :
+      !!
+      !! ** Method  :
       !!----------------------------------------------------------------------
 
@@ -468,11 +468,11 @@
 
       !!----------------------------------------------------------------------
-      
+
       IF( kt == nit000 ) THEN
-         REWIND( numnam_ref )              ! Namelist namzdf_mldzint in reference namelist 
+         REWIND( numnam_ref )              ! Namelist namzdf_mldzint in reference namelist
          READ  ( numnam_ref, namzdf_mldzint, IOSTAT = ios, ERR = 901)
 901      IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_mldzint in reference namelist' )
 
-         REWIND( numnam_cfg )              ! Namelist namzdf_mldzint in configuration namelist 
+         REWIND( numnam_cfg )              ! Namelist namzdf_mldzint in configuration namelist
          READ  ( numnam_cfg, namzdf_mldzint, IOSTAT = ios, ERR = 902 )
 902      IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_mldzint in configuration namelist' )
@@ -523,43 +523,44 @@
 
 
-   SUBROUTINE zdf_mxl_kara( kt ) 
-      !!---------------------------------------------------------------------------------- 
-      !!                    ***  ROUTINE zdf_mxl_kara  *** 
-      !                                                                        
-      !   Calculate mixed layer depth according to the definition of          
-      !   Kara et al, 2000, JGR, 105, 16803.  
-      ! 
-      !   If mld_type=1 the mixed layer depth is calculated as the depth at which the  
-      !   density has increased by an amount equivalent to a temperature difference of  
-      !   0.8C at the surface. 
-      ! 
-      !   For other values of mld_type the mixed layer is calculated as the depth at  
-      !   which the temperature differs by 0.8C from the surface temperature.  
-      !                                                                        
-      !   Original version: David Acreman                                      
-      ! 
+   SUBROUTINE zdf_mxl_kara( kt )
+      !!----------------------------------------------------------------------------------
+      !!                    ***  ROUTINE zdf_mxl_kara  ***
+      !
+      !   Calculate mixed layer depth according to the definition of
+      !   Kara et al, 2000, JGR, 105, 16803.
+      !
+      !   If mld_type=1 the mixed layer depth is calculated as the depth at which the
+      !   density has increased by an amount equivalent to a temperature difference of
+      !   0.8C at the surface.
+      !
+      !   For other values of mld_type the mixed layer is calculated as the depth at
+      !   which the temperature differs by 0.8C from the surface temperature.
+      !
+      !   Original version: David Acreman
+      !
       !!-----------------------------------------------------------------------------------
-     
-      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index 
- 
+
+      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
+
       NAMELIST/namzdf_karaml/ ln_kara,jpmld_type, ppz_ref, ppdT_crit, &
-      &                       ppiso_frac, ln_kara_write25h 
- 
-      ! Local variables                                                        
-      REAL, DIMENSION(jpi,jpk) :: ppzdep      ! depth for use in calculating d(rho) 
-      REAL(wp), DIMENSION(jpi,jpj,jpts) :: ztsn_2d  !Local version of tsn 
-      LOGICAL :: ll_found(jpi,jpj)              ! Is T_b to be found by interpolation ? 
-      LOGICAL :: ll_belowml(jpi,jpj,jpk)        ! Flag points below mixed layer when ll_found=F 
-      INTEGER :: ji, jj, jk                     ! loop counter 
-      INTEGER :: ik_ref(jpi,jpj)                ! index of reference level 
-      INTEGER :: ik_iso(jpi,jpj)                ! index of last uniform temp level 
-      REAL    :: zT(jpi,jpj,jpk)                ! Temperature or denisty 
-      REAL    :: zT_ref(jpi,jpj)                ! reference temperature 
-      REAL    :: zT_b                           ! base temperature 
-      REAL    :: zdTdz(jpi,jpj,jpk-2)           ! gradient of zT 
-      REAL    :: zmoddT(jpi,jpj,jpk-2)          ! Absolute temperature difference 
-      REAL    :: zdz                            ! depth difference 
-      REAL    :: zdT                            ! temperature difference 
-      REAL    :: zdelta_T(jpi,jpj)              ! difference critereon 
+      &                       ppiso_frac, ln_kara_write25h
+
+      ! Local variables
+      REAL, DIMENSION(jpi,jpk) :: ppzdep      ! depth for use in calculating d(rho)
+      REAL(wp), DIMENSION(jpi,jpj,jpts) :: ztsn_2d  !Local version of tsn
+      INTEGER, DIMENSION(jpi,jpj) :: ikmt       ! number of active tracer levels
+      LOGICAL :: ll_found(jpi,jpj)              ! Is T_b to be found by interpolation ?
+      LOGICAL :: ll_belowml(jpi,jpj,jpk)        ! Flag points below mixed layer when ll_found=F
+      INTEGER :: ji, jj, jk                     ! loop counter
+      INTEGER :: ik_ref(jpi,jpj)                ! index of reference level
+      INTEGER :: ik_iso(jpi,jpj)                ! index of last uniform temp level
+      REAL    :: zT(jpi,jpj,jpk)                ! Temperature or denisty
+      REAL    :: zT_ref(jpi,jpj)                ! reference temperature
+      REAL    :: zT_b                           ! base temperature
+      REAL    :: zdTdz(jpi,jpj,jpk-2)           ! gradient of zT
+      REAL    :: zmoddT(jpi,jpj,jpk-2)          ! Absolute temperature difference
+      REAL    :: zdz                            ! depth difference
+      REAL    :: zdT                            ! temperature difference
+      REAL    :: zdelta_T(jpi,jpj)              ! difference critereon
       REAL    :: zRHO1(jpi,jpj), zRHO2(jpi,jpj) ! Densities
       INTEGER, SAVE :: idt, i_steps
@@ -567,16 +568,16 @@
       INTEGER :: ios                 ! Local integer output status for namelist read
 
-     
- 
-      !!------------------------------------------------------------------------------------- 
- 
-      IF( kt == nit000 ) THEN 
-         ! Default values 
+
+
+      !!-------------------------------------------------------------------------------------
+
+      IF( kt == nit000 ) THEN
+         ! Default values
          ln_kara      = .FALSE.
-         ln_kara_write25h = .FALSE. 
-         jpmld_type   = 1     
-         ppz_ref      = 10.0 
-         ppdT_crit    = 0.2  
-         ppiso_frac   = 0.1   
+         ln_kara_write25h = .FALSE.
+         jpmld_type   = 1
+         ppz_ref      = 10.0
+         ppdT_crit    = 0.2
+         ppiso_frac   = 0.1
          ! Read namelist
          REWIND ( numnam_ref )
@@ -588,17 +589,17 @@
 902      IF( ios > 0 ) CALL ctl_nam ( ios , 'namzdf_karaml in configuration namelist' )
          IF(lwm) WRITE ( numond, namzdf_karaml )
- 
-
-         WRITE(numout,*) '===== Kara mixed layer =====' 
+
+
+         WRITE(numout,*) '===== Kara mixed layer ====='
          WRITE(numout,*) 'ln_kara = ',    ln_kara
-         WRITE(numout,*) 'jpmld_type = ', jpmld_type 
-         WRITE(numout,*) 'ppz_ref = ',    ppz_ref 
-         WRITE(numout,*) 'ppdT_crit = ',  ppdT_crit 
+         WRITE(numout,*) 'jpmld_type = ', jpmld_type
+         WRITE(numout,*) 'ppz_ref = ',    ppz_ref
+         WRITE(numout,*) 'ppdT_crit = ',  ppdT_crit
          WRITE(numout,*) 'ppiso_frac = ', ppiso_frac
          WRITE(numout,*) 'ln_kara_write25h = ', ln_kara_write25h
-         WRITE(numout,*) '============================' 
-      
+         WRITE(numout,*) '============================'
+
          IF ( .NOT.ln_kara ) THEN
-            WRITE(numout,*) "ln_kara not set; Kara mixed layer not calculated" 
+            WRITE(numout,*) "ln_kara not set; Kara mixed layer not calculated"
          ELSE
             IF (.NOT. ALLOCATED(hmld_kara) ) ALLOCATE( hmld_kara(jpi,jpj) )
@@ -615,136 +616,139 @@
 
               idt = INT(rdt)
-               IF( MOD( 3600,idt) == 0 ) THEN 
-                  i_steps = 3600 / idt  
-               ELSE 
+               IF( MOD( 3600,idt) == 0 ) THEN
+                  i_steps = 3600 / idt
+               ELSE
                   CALL ctl_stop('STOP', &
                   & 'zdf_mxl_kara: timestep must give MOD(3600,rdt)'// &
-                  & ' = 0 otherwise no hourly values are possible') 
-               ENDIF  
+                  & ' = 0 otherwise no hourly values are possible')
+               ENDIF
             ENDIF
          ENDIF
       ENDIF
-      
+
       IF ( ln_kara ) THEN
-         
+
          !set critical ln_kara
          ztsn_2d = tsn(:,:,1,:)
          ztsn_2d(:,:,jp_tem) = ztsn_2d(:,:,jp_tem) + ppdT_crit
-     
-         ! Set the mixed layer depth criterion at each grid point 
+
+         ! Set the mixed layer depth criterion at each grid point
          ppzdep = 0._wp
-         IF( jpmld_type == 1 ) THEN                                         
+         IF( jpmld_type == 1 ) THEN
             CALL eos ( tsn(:,:,1,:), &
-            &          ppzdep(:,:), zRHO1(:,:) ) 
+            &          ppzdep(:,:), zRHO1(:,:) )
             CALL eos ( ztsn_2d(:,:,:), &
-            &           ppzdep(:,:), zRHO2(:,:) ) 
-            zdelta_T(:,:) = abs( zRHO1(:,:) - zRHO2(:,:) ) * rau0 
-            ! RHO from eos (2d version) doesn't calculate north or east halo: 
-            CALL lbc_lnk( 'zdf_mxl_kara',zdelta_T, 'T', 1. ) 
-            zT(:,:,:) = rhop(:,:,:) 
-         ELSE 
-            zdelta_T(:,:) = ppdT_crit                      
-            zT(:,:,:) = tsn(:,:,:,jp_tem)                           
-         ENDIF 
-     
-         ! Calculate the gradient of zT and absolute difference for use later 
-         DO jk = 1 ,jpk-2 
-            zdTdz(:,:,jk)  =    ( zT(:,:,jk+1) - zT(:,:,jk) ) / e3w_n(:,:,jk+1) 
-            zmoddT(:,:,jk) = abs( zT(:,:,jk+1) - zT(:,:,jk) ) 
-         END DO 
-     
-         ! Find density/temperature at the reference level (Kara et al use 10m).          
-         ! ik_ref is the index of the box centre immediately above or at the reference level 
-         ! Find ppz_ref in the array of model level depths and find the ref    
-         ! density/temperature by linear interpolation.                                   
+            &           ppzdep(:,:), zRHO2(:,:) )
+            zdelta_T(:,:) = abs( zRHO1(:,:) - zRHO2(:,:) ) * rau0
+            ! RHO from eos (2d version) doesn't calculate north or east halo:
+            CALL lbc_lnk( 'zdf_mxl_kara',zdelta_T, 'T', 1. )
+            zT(:,:,:) = rhop(:,:,:)
+         ELSE
+            zdelta_T(:,:) = ppdT_crit
+            zT(:,:,:) = tsn(:,:,:,jp_tem)
+         ENDIF
+
+         ! Calculate the gradient of zT and absolute difference for use later
+         DO jk = 1 ,jpk-2
+            zdTdz(:,:,jk)  =    ( zT(:,:,jk+1) - zT(:,:,jk) ) / e3w_n(:,:,jk+1)
+            zmoddT(:,:,jk) = abs( zT(:,:,jk+1) - zT(:,:,jk) )
+         END DO
+
+         ! Find density/temperature at the reference level (Kara et al use 10m).
+         ! ik_ref is the index of the box centre immediately above or at the reference level
+         ! Find ppz_ref in the array of model level depths and find the ref
+         ! density/temperature by linear interpolation.
          ik_ref = -1
-         DO jk = jpkm1, 2, -1 
-            WHERE( gdept_n(:,:,jk) > ppz_ref ) 
-               ik_ref(:,:) = jk - 1 
+         DO jk = jpkm1, 2, -1
+            WHERE( gdept_n(:,:,jk) > ppz_ref )
+               ik_ref(:,:) = jk - 1
                zT_ref(:,:) = zT(:,:,jk-1) + &
-               &             zdTdz(:,:,jk-1) * ( ppz_ref - gdept_n(:,:,jk-1) ) 
-            ENDWHERE 
+               &             zdTdz(:,:,jk-1) * ( ppz_ref - gdept_n(:,:,jk-1) )
+            ENDWHERE
          END DO
          IF ( ANY( ik_ref  < 0 ) .OR. ANY( ik_ref  > jpkm1 ) ) THEN
             CALL ctl_stop( "STOP", &
-            & "zdf_mxl_kara: unable to find reference level for kara ML" ) 
+            & "zdf_mxl_kara: unable to find reference level for kara ML" )
          ELSE
-            ! If the first grid box centre is below the reference level then use the 
-            ! top model level to get zT_ref 
-            WHERE( gdept_n(:,:,1) > ppz_ref )  
-               zT_ref = zT(:,:,1) 
-               ik_ref = 1 
-            ENDWHERE 
-     
-            ! Search for a uniform density/temperature region where adjacent levels          
-            ! differ by less than ppiso_frac * deltaT.                                      
-            ! ik_iso is the index of the last level in the uniform layer  
-            ! ll_found indicates whether the mixed layer depth can be found by interpolation 
-            ik_iso(:,:)   = ik_ref(:,:) 
-            ll_found(:,:) = .false. 
-            DO jj = 1, nlcj 
-               DO ji = 1, nlci 
-                 !CDIR NOVECTOR 
-                  DO jk = ik_ref(ji,jj), mbkt(ji,jj)-1  !mbathy(ji,jj)-1 
-                     IF( zmoddT(ji,jj,jk) > ( ppiso_frac * zdelta_T(ji,jj) ) ) THEN 
-                        ik_iso(ji,jj)   = jk 
-                        ll_found(ji,jj) = ( zmoddT(ji,jj,jk) > zdelta_T(ji,jj) ) 
-                        EXIT 
-                     ENDIF 
-                  END DO 
-               END DO 
-            END DO 
-     
-            ! Use linear interpolation to find depth of mixed layer base where possible 
-            hmld_kara(:,:) = ppz_ref 
-            DO jj = 1, jpj 
-               DO ji = 1, jpi 
-                  IF( ll_found(ji,jj) .and. tmask(ji,jj,1) == 1.0 ) THEN 
-                     zdz =  abs( zdelta_T(ji,jj) / zdTdz(ji,jj,ik_iso(ji,jj)) ) 
-                     hmld_kara(ji,jj) = gdept_n(ji,jj,ik_iso(ji,jj)) + zdz 
-                  ENDIF 
-               END DO 
-            END DO 
-     
-            ! If ll_found = .false. then calculate MLD using difference of zdelta_T    
-            ! from the reference density/temperature 
-     
-            ! Prevent this section from working on land points 
-            WHERE( tmask(:,:,1) /= 1.0 ) 
-               ll_found = .true. 
-            ENDWHERE 
-     
-            DO jk = 1, jpk 
+            ! If the first grid box centre is below the reference level then use the
+            ! top model level to get zT_ref
+            WHERE( gdept_n(:,:,1) > ppz_ref )
+               zT_ref = zT(:,:,1)
+               ik_ref = 1
+            ENDWHERE
+
+            ! The number of active tracer levels is 1 less than the number of active w levels
+            ikmt(:,:) = mbkt(:,:) - 1
+
+            ! Search for a uniform density/temperature region where adjacent levels
+            ! differ by less than ppiso_frac * deltaT.
+            ! ik_iso is the index of the last level in the uniform layer
+            ! ll_found indicates whether the mixed layer depth can be found by interpolation
+            ik_iso(:,:)   = ik_ref(:,:)
+            ll_found(:,:) = .false.
+            DO jj = 1, nlcj
+               DO ji = 1, nlci
+                 !CDIR NOVECTOR
+                  DO jk = ik_ref(ji,jj),  ikmt(ji,jj)-1 ! mbkt(ji,jj)-1  !mbathy(ji,jj)-1
+                     IF( zmoddT(ji,jj,jk) > ( ppiso_frac * zdelta_T(ji,jj) ) ) THEN
+                        ik_iso(ji,jj)   = jk
+                        ll_found(ji,jj) = ( zmoddT(ji,jj,jk) > zdelta_T(ji,jj) )
+                        EXIT
+                     ENDIF
+                  END DO
+               END DO
+            END DO
+
+            ! Use linear interpolation to find depth of mixed layer base where possible
+            hmld_kara(:,:) = ppz_ref
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  IF( ll_found(ji,jj) .and. tmask(ji,jj,1) == 1.0 ) THEN
+                     zdz =  abs( zdelta_T(ji,jj) / zdTdz(ji,jj,ik_iso(ji,jj)) )
+                     hmld_kara(ji,jj) = gdept_n(ji,jj,ik_iso(ji,jj)) + zdz
+                  ENDIF
+               END DO
+            END DO
+
+            ! If ll_found = .false. then calculate MLD using difference of zdelta_T
+            ! from the reference density/temperature
+
+            ! Prevent this section from working on land points
+            WHERE( tmask(:,:,1) /= 1.0 )
+               ll_found = .true.
+            ENDWHERE
+
+            DO jk = 1, jpk
                ll_belowml(:,:,jk) = abs( zT(:,:,jk) - zT_ref(:,:) ) >= &
                & zdelta_T(:,:)
-            END DO 
-     
-            ! Set default value where interpolation cannot be used (ll_found=false)  
-            DO jj = 1, jpj 
-               DO ji = 1, jpi 
+            END DO
+
+            ! Set default value where interpolation cannot be used (ll_found=false)
+            DO jj = 1, jpj
+               DO ji = 1, jpi
                   IF( .NOT. ll_found(ji,jj) )  &
-                  &   hmld_kara(ji,jj) = gdept_n(ji,jj,mbkt(ji,jj))! mbathy(ji,jj)) 
-               END DO 
-            END DO 
-     
-            DO jj = 1, jpj 
-               DO ji = 1, jpi 
-                  !CDIR NOVECTOR 
-                  DO jk = ik_ref(ji,jj)+1, mbkt(ji,jj) !mbathy(ji,jj) 
-                     IF( ll_found(ji,jj) ) EXIT 
-                     IF( ll_belowml(ji,jj,jk) ) THEN                
+                  &   hmld_kara(ji,jj) = gdept_n(ji,jj,ikmt(ji,jj)) ! mbkt(ji,jj))! mbathy(ji,jj))
+               END DO
+            END DO
+
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  !CDIR NOVECTOR
+                  DO jk = ik_ref(ji,jj)+1, ikmt(ji,jj) !mbkt(ji,jj) !mbathy(ji,jj)
+                     IF( ll_found(ji,jj) ) EXIT
+                     IF( ll_belowml(ji,jj,jk) ) THEN
                         zT_b = zT_ref(ji,jj) + zdelta_T(ji,jj) * &
-                        &      SIGN(1.0, zdTdz(ji,jj,jk-1) ) 
-                        zdT  = zT_b - zT(ji,jj,jk-1)                                      
-                        zdz  = zdT / zdTdz(ji,jj,jk-1)                                       
-                        hmld_kara(ji,jj) = gdept_n(ji,jj,jk-1) + zdz 
-                        EXIT                                                   
-                     ENDIF 
-                  END DO 
-               END DO 
-            END DO 
-     
-            hmld_kara(:,:) = hmld_kara(:,:) * tmask(:,:,1) 
- 
+                        &      SIGN(1.0, zdTdz(ji,jj,jk-1) )
+                        zdT  = zT_b - zT(ji,jj,jk-1)
+                        zdz  = zdT / zdTdz(ji,jj,jk-1)
+                        hmld_kara(ji,jj) = gdept_n(ji,jj,jk-1) + zdz
+                        EXIT
+                     ENDIF
+                  END DO
+               END DO
+            END DO
+
+            hmld_kara(:,:) = hmld_kara(:,:) * tmask(:,:,1)
+
             IF(  ln_kara_write25h  ) THEN
                !Double IF required as i_steps not defined if ln_kara_write25h =
@@ -756,17 +760,17 @@
                ENDIF
             ENDIF
- 
-!#if defined key_iomput 
-            IF( kt /= nit000 ) THEN 
-               CALL iom_put( "mldkara"  , hmld_kara )   
+
+!#if defined key_iomput
+            IF( kt /= nit000 ) THEN
+               CALL iom_put( "mldkara"  , hmld_kara )
                IF( ( MOD( i_cnt_25h, 25) == 0 ) .AND.  ln_kara_write25h ) &
                   CALL iom_put( "kara25h"  , ( hmld_kara_25h / 25._wp ) )
             ENDIF
 !#endif
- 
+
          ENDIF
       ENDIF
-       
-   END SUBROUTINE zdf_mxl_kara 
+
+   END SUBROUTINE zdf_mxl_kara
 
    !!======================================================================
