Index: /branches/UKMO/dev_r5107_mld_zint/NEMOGCM/CONFIG/SHARED/field_def.xml
===================================================================
--- /branches/UKMO/dev_r5107_mld_zint/NEMOGCM/CONFIG/SHARED/field_def.xml	(revision 6071)
+++ /branches/UKMO/dev_r5107_mld_zint/NEMOGCM/CONFIG/SHARED/field_def.xml	(revision 6072)
@@ -53,5 +53,14 @@
          <field id="mldr10_1max"  long_name="Max of Mixed Layer Depth (dsigma = 0.01 wrt 10m)"   field_ref="mldr10_1"   operation="maximum"                                                                          />
          <field id="mldr10_1min"  long_name="Min of Mixed Layer Depth (dsigma = 0.01 wrt 10m)"   field_ref="mldr10_1"   operation="minimum"                                                                          />
-         <field id="mldzint"      long_name="Mixed Layer Depth interpolated"                     standard_name="ocean_mixed_layer_thickness_defined_by_sigma_theta"                                unit="m"          />
+         <field id="mldzint_1"    long_name="Mixed Layer Depth interpolated"                     standard_name="ocean_mixed_layer_thickness"                                                       unit="m"          />
+         <field id="mldzint_2"     field_ref="mldzint_1" />
+         <field id="mldzint_3"     field_ref="mldzint_1" />
+         <field id="mldzint_4"     field_ref="mldzint_1" />
+         <field id="mldzint_5"     field_ref="mldzint_1" />
+         <field id="mldhtc_1"     long_name="Mixed Layer Depth integrated heat content"          standard_name="integral_of_sea_water_potential_temperature_wrt_depth_expressed_as_heat_content"   unit="J/m2"       />
+         <field id="mldhtc_2"      field_ref="mldhtc_1"  />
+         <field id="mldhtc_3"      field_ref="mldhtc_1"  />
+         <field id="mldhtc_4"      field_ref="mldhtc_1"  />
+         <field id="mldhtc_5"      field_ref="mldhtc_1"  />
          <field id="heatc"        long_name="Heat content vertically integrated"                 standard_name="integral_of_sea_water_potential_temperature_wrt_depth_expressed_as_heat_content"   unit="J/m2"       />
          <field id="saltc"        long_name="Salt content vertically integrated"                                                                                                                   unit="1e-3*kg/m2" />
Index: /branches/UKMO/dev_r5107_mld_zint/NEMOGCM/CONFIG/SHARED/namelist_ref
===================================================================
--- /branches/UKMO/dev_r5107_mld_zint/NEMOGCM/CONFIG/SHARED/namelist_ref	(revision 6071)
+++ /branches/UKMO/dev_r5107_mld_zint/NEMOGCM/CONFIG/SHARED/namelist_ref	(revision 6072)
@@ -1014,8 +1014,13 @@
 &namzdf_mldzint    !   Parameters for vertically-interpolated mixed-layer depth diagnostic
 !------------------------------------------------------------------------------------------
-         nn_mld_type   = 1      ! mixed layer type
-         rn_zref       = 10.0   ! depth of initial reference temperature
-         rn_dT_crit    = 0.2    ! critical temperature difference
-         rn_iso_frac   = 0.1    ! fraction of critical temperature difference used
+   nn_mld_diag = 0         !  Number of MLD diagnostics to use from below
+
+!              ! MLD criterion ! Reference ! Finite difference ! Gradient layer !
+!              ! type          ! depth     ! criterion         ! criterion      !
+   sn_mld1     =       1       ,    10.0   ,        0.2        ,       0.1
+   sn_mld2     =       ''      ,      ''   ,         ''        ,        ''
+   sn_mld3     =       ''      ,      ''   ,         ''        ,        ''
+   sn_mld4     =       ''      ,      ''   ,         ''        ,        ''
+   sn_mld5     =       ''      ,      ''   ,         ''        ,        ''
 /
 
Index: /branches/UKMO/dev_r5107_mld_zint/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfmxl.F90
===================================================================
--- /branches/UKMO/dev_r5107_mld_zint/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfmxl.F90	(revision 6071)
+++ /branches/UKMO/dev_r5107_mld_zint/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfmxl.F90	(revision 6072)
@@ -33,5 +33,6 @@
    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_zint  !: vertically-interpolated mixed layer depth   [m] 
+   REAL(wp), PUBLIC, ALLOCATABLE,       DIMENSION(:,:) ::   hmld_zint  !: vertically-interpolated mixed layer depth   [m] 
+   REAL(wp), PUBLIC, ALLOCATABLE,       DIMENSION(:,:) ::   htc_mld    ! Heat content of hmld_zint
    LOGICAL, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)    :: ll_found   ! Is T_b to be found by interpolation ? 
    LOGICAL, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  :: ll_belowml ! Flag points below mixed layer when ll_found=F
@@ -40,9 +41,10 @@
    REAL(wp)         ::   avt_c = 5.e-4_wp   ! Kz criterion for the turbocline depth
 
-   ! Namelist variables for  namzdf_mldzint
-   INTEGER          :: nn_mld_type         ! mixed layer type            
-   REAL(wp)         :: rn_zref            ! depth of initial T_ref
-   REAL(wp)         :: rn_dT_crit          ! Critical temp diff 
-   REAL(wp)         :: rn_iso_frac         ! Fraction of rn_dT_crit used 
+   TYPE, PUBLIC :: MXL_ZINT   !: Structure for MLD defs
+      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 used
+   END TYPE MXL_ZINT
 
    !! * Substitutions
@@ -62,4 +64,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 )
          !
@@ -149,8 +152,5 @@
       
       ! Vertically-interpolated mixed-layer depth diagnostic
-      IF( iom_use( "mldzint" ) ) THEN
-         CALL zdf_mxl_zint( kt )
-         CALL iom_put( "mldzint" , hmld_zint )
-      ENDIF
+      CALL zdf_mxl_zint( kt )
 
       IF(ln_ctl)   CALL prt_ctl( tab2d_1=REAL(nmln,wp), clinfo1=' nmln : ', tab2d_2=hmlp, clinfo2=' hmlp : ', ovlap=1 )
@@ -162,7 +162,7 @@
    END SUBROUTINE zdf_mxl
 
-   SUBROUTINE zdf_mxl_zint( kt ) 
+   SUBROUTINE zdf_mxl_zint_mld( sf ) 
       !!---------------------------------------------------------------------------------- 
-      !!                    ***  ROUTINE zdf_mxl_zint  *** 
+      !!                    ***  ROUTINE zdf_mxl_zint_mld  *** 
       !                                                                        
       !   Calculate vertically-interpolated mixed layer depth diagnostic. 
@@ -184,6 +184,12 @@
       !!----------------------------------------------------------------------------------- 
 
-      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
-      !
+      TYPE(MXL_ZINT), INTENT(in)  :: sf
+
+      ! Diagnostic criteria
+      INTEGER   :: nn_mld_type   ! mixed layer type     
+      REAL(wp)  :: rn_zref       ! depth of initial T_ref
+      REAL(wp)  :: rn_dT_crit    ! Critical temp diff
+      REAL(wp)  :: rn_iso_frac   ! Fraction of rn_dT_crit used
+
       ! Local variables
       INTEGER, POINTER, DIMENSION(:,:) :: ikmt          ! number of active tracer levels 
@@ -201,7 +207,4 @@
       REAL, POINTER, DIMENSION(:,:)    :: zRHO1, zRHO2  ! Densities 
       INTEGER :: ji, jj, jk                             ! loop counter 
-      INTEGER :: ios
-
-      NAMELIST/namzdf_mldzint/ nn_mld_type, rn_zref, rn_dT_crit, rn_iso_frac
 
       !!------------------------------------------------------------------------------------- 
@@ -210,25 +213,16 @@
       CALL wrk_alloc( jpi, jpj, ppzdep, zT_ref, zdelta_T, zRHO1, zRHO2 ) 
       CALL wrk_alloc( jpi, jpj, jpk, zT, zdTdz, zmoddT ) 
- 
-      IF( kt == nit000 ) THEN
-         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', lwp )
-
-         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', lwp )
-         IF(lwm) WRITE ( numond, namzdf_mldzint )
-
-         WRITE(numout,*) '===== Vertically-interpolated mixed layer ====='
-         WRITE(numout,*) 'nn_mld_type = ',nn_mld_type
-         WRITE(numout,*) 'rn_zref = ',rn_zref
-         WRITE(numout,*) 'rn_dT_crit = ',rn_dT_crit
-         WRITE(numout,*) 'rn_iso_frac = ',rn_iso_frac
-         WRITE(numout,*) '==============================================='
-      ENDIF
- 
+
+      ! Unpack structure
+      nn_mld_type = sf%mld_type
+      rn_zref     = sf%zref
+      rn_dT_crit  = sf%dT_crit
+      rn_iso_frac = sf%iso_frac
+
       ! Set the mixed layer depth criterion at each grid point 
-      IF (nn_mld_type == 1) THEN                                         
+      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(:,:) ) 
@@ -345,4 +339,160 @@
       CALL wrk_dealloc( jpi,jpj, jpk, zT, zdTdz, zmoddT ) 
       ! 
+   END SUBROUTINE zdf_mxl_zint_mld
+
+   SUBROUTINE zdf_mxl_zint_htc( kt )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE zdf_mxl_zint_htc  ***
+      !! 
+      !! ** Purpose :   
+      !!
+      !! ** Method  :   
+      !!----------------------------------------------------------------------
+
+      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
+
+      INTEGER :: ji, jj, jk
+      INTEGER :: ikmax
+      REAL(wp) :: zc, zcoef
+      !
+      INTEGER,  ALLOCATABLE, DIMENSION(:,:) ::   ilevel
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zthick_0, zthick
+
+      !!----------------------------------------------------------------------
+
+      IF( .NOT. ALLOCATED(ilevel) ) THEN
+         ALLOCATE( ilevel(jpi,jpj), zthick_0(jpi,jpj), &
+         &         zthick(jpi,jpj), STAT=ji )
+         IF( lk_mpp  )   CALL mpp_sum(ji)
+         IF( ji /= 0 )   CALL ctl_stop( 'STOP', 'zdf_mxl_zint_htc : unable to allocate arrays' )
+      ENDIF
+
+      ! Find last whole model T level above the MLD
+      ilevel(:,:)   = 0
+      zthick_0(:,:) = 0._wp
+
+      DO jk = 1, jpkm1  
+         DO jj = 1, jpj
+            DO ji = 1, jpi                    
+               zthick_0(ji,jj) = zthick_0(ji,jj) + fse3t(ji,jj,jk)
+               IF( zthick_0(ji,jj) < hmld_zint(ji,jj) )   ilevel(ji,jj) = jk
+            END DO
+         END DO
+         WRITE(numout,*) 'zthick_0(jk =',jk,') =',zthick_0(2,2)
+         WRITE(numout,*) 'fsdepw(jk+1 =',jk+1,') =',fsdepw(2,2,jk+1)
+      END DO
+
+      ! Surface boundary condition
+      IF( lk_vvl ) THEN   ;   zthick(:,:) = 0._wp       ;   htc_mld(:,:) = 0._wp                                   
+      ELSE                ;   zthick(:,:) = sshn(:,:)   ;   htc_mld(:,:) = tsn(:,:,1,jp_tem) * sshn(:,:) * tmask(:,:,1)   
+      ENDIF
+
+      ! Deepest whole T level above the MLD
+      ikmax = MIN( MAXVAL( ilevel(:,:) ), jpkm1 )
+
+      ! Integration down to last whole model T level
+      DO jk = 1, ikmax
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               zc = fse3t(ji,jj,jk) * REAL( MIN( MAX( 0, ilevel(ji,jj) - jk + 1 ) , 1  )  )    ! 0 below ilevel
+               zthick(ji,jj) = zthick(ji,jj) + zc
+               htc_mld(ji,jj) = htc_mld(ji,jj) + zc * tsn(ji,jj,jk,jp_tem) * tmask(ji,jj,jk)
+            END DO
+         END DO
+      END DO
+
+      ! Subsequent partial T level
+      zthick(:,:) = hmld_zint(:,:) - zthick(:,:)   !   remaining thickness to reach MLD
+
+      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)  & 
+      &                      * MIN( fse3t(ji,jj,ilevel(ji,jj)+1), zthick(ji,jj) ) * tmask(ji,jj,ilevel(ji,jj)+1)
+         END DO
+      END DO
+
+      WRITE(numout,*) 'htc_mld(after) =',htc_mld(2,2)
+
+      ! Convert to heat content
+      zcoef = rau0 * rcp
+      htc_mld(:,:) = zcoef * htc_mld(:,:)
+
+   END SUBROUTINE zdf_mxl_zint_htc
+
+   SUBROUTINE zdf_mxl_zint( kt )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE zdf_mxl_zint  ***
+      !! 
+      !! ** Purpose :   
+      !!
+      !! ** Method  :   
+      !!----------------------------------------------------------------------
+
+      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
+
+      INTEGER :: ios
+      INTEGER :: jn
+
+      INTEGER :: nn_mld_diag = 0    ! number of diagnostics
+
+      CHARACTER(len=1) :: cmld
+
+      TYPE(MXL_ZINT) :: sn_mld1, sn_mld2, sn_mld3, sn_mld4, sn_mld5
+      TYPE(MXL_ZINT), SAVE, DIMENSION(5) ::   mld_diags
+
+      NAMELIST/namzdf_mldzint/ nn_mld_diag, sn_mld1, sn_mld2, sn_mld3, sn_mld4, sn_mld5
+
+      !!----------------------------------------------------------------------
+      
+      IF( kt == nit000 ) THEN
+         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', lwp )
+
+         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', lwp )
+         IF(lwm) WRITE ( numond, namzdf_mldzint )
+
+         IF( nn_mld_diag > 5 )   CALL ctl_stop( 'STOP', 'zdf_mxl_ini: Specify no more than 5 MLD definitions' )
+
+         mld_diags(1) = sn_mld1
+         mld_diags(2) = sn_mld2
+         mld_diags(3) = sn_mld3
+         mld_diags(4) = sn_mld4
+         mld_diags(5) = sn_mld5
+
+         IF( nn_mld_diag > 0 ) THEN
+            WRITE(numout,*) '=============== Vertically-interpolated mixed layer ================'
+            WRITE(numout,*) '(Diagnostic number, nn_mld_type, rn_zref, rn_dT_crit, rn_iso_frac)'
+            DO jn = 1, nn_mld_diag
+               WRITE(numout,*) 'MLD criterion',jn,':'
+               WRITE(numout,*) '    nn_mld_type =', mld_diags(jn)%mld_type
+               WRITE(numout,*) '    rn_zref ='    , mld_diags(jn)%zref
+               WRITE(numout,*) '    rn_dT_crit =' , mld_diags(jn)%dT_crit
+               WRITE(numout,*) '    rn_iso_frac =', mld_diags(jn)%iso_frac
+            END DO
+            WRITE(numout,*) '===================================================================='
+         ENDIF
+      ENDIF
+
+      IF( nn_mld_diag > 0 ) THEN
+         DO jn = 1, nn_mld_diag
+            WRITE(cmld,'(I1)') jn
+            IF( iom_use( "mldzint_"//cmld ) .OR. iom_use( "mldhtc_"//cmld ) ) THEN
+               CALL zdf_mxl_zint_mld( mld_diags(jn) )
+
+               IF( iom_use( "mldzint_"//cmld ) ) THEN
+                  CALL iom_put( "mldzint_"//cmld, hmld_zint(:,:) )
+               ENDIF
+
+               IF( iom_use( "mldhtc_"//cmld ) )  THEN
+                  CALL zdf_mxl_zint_htc( kt )
+                  CALL iom_put( "mldhtc_"//cmld , htc_mld(:,:)   )
+               ENDIF
+            ENDIF
+         END DO
+      ENDIF
+
    END SUBROUTINE zdf_mxl_zint
 
