Index: branches/UKMO/dev_r6501_GO6_package_trunk/NEMOGCM/CONFIG/ORCA2_LIM3/EXP00/iodef.xml
===================================================================
--- branches/UKMO/dev_r6501_GO6_package_trunk/NEMOGCM/CONFIG/ORCA2_LIM3/EXP00/iodef.xml	(revision 6507)
+++ branches/UKMO/dev_r6501_GO6_package_trunk/NEMOGCM/CONFIG/ORCA2_LIM3/EXP00/iodef.xml	(revision 6825)
@@ -61,4 +61,6 @@
 	<file id="file2" name_suffix="_SBC" description="surface fluxes variables" > <!-- time step automaticaly defined based on nn_fsbc -->
 	  <field field_ref="empmr"        name="wfo"      />
+          <field field_ref="emp_oce"      name="emp_oce"  long_name="Evap minus Precip over ocean"              />
+          <field field_ref="emp_ice"      name="emp_ice"  long_name="Evap minus Precip over ice"                />
 	  <field field_ref="qsr_oce"      name="qsr_oce"  />
 	  <field field_ref="qns_oce"      name="qns_oce"  />
Index: branches/UKMO/dev_r6501_GO6_package_trunk/NEMOGCM/CONFIG/SHARED/field_def.xml
===================================================================
--- branches/UKMO/dev_r6501_GO6_package_trunk/NEMOGCM/CONFIG/SHARED/field_def.xml	(revision 6507)
+++ branches/UKMO/dev_r6501_GO6_package_trunk/NEMOGCM/CONFIG/SHARED/field_def.xml	(revision 6825)
@@ -61,4 +61,5 @@
          <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="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" />
@@ -204,4 +205,5 @@
 
          <!-- * variable related to ice shelf forcing * -->
+         <field id="berg_calve"   long_name="Iceberg calving"                               unit="kg/m2/s"  />
          <field id="fwfisf"       long_name="Ice shelf melting"                                            unit="Kg/m2/s"  />
          <field id="qisf"         long_name="Ice Shelf Heat Flux"                                          unit="W/m2"     />
@@ -501,4 +503,8 @@
         <field id="aht2d"        long_name="lateral eddy diffusivity"       standard_name="ocean_tracer_xy_laplacian_diffusivity"      unit="m2/s"   grid_ref="grid_W_2D" />
         <field id="aht2d_eiv"    long_name="EIV lateral eddy diffusivity"   standard_name="ocean_tracer_bolus_laplacian_diffusivity"   unit="m2/s"   grid_ref="grid_W_2D" />
+        <!-- variables added for nn_etau revision testing -->
+        <field id="htau"         long_name="e-decay length scale of nn_etau TKE"       unit="m"        grid_ref="grid_W_2D" />
+        <field id="tke"          long_name="Turbulent Kinetic Energy"                  unit="m2/s2"    />
+        <field id="tke_niw"      long_name="Turbulent Kinetic Energy due to near-inertial wave breaking"  unit="m2/s2" />
       </field_group>
         
@@ -507,4 +513,23 @@
 	      <field id="ahmf_2d"      long_name=" surface f-eddy viscosity coefficient"   unit="m2/s or m4/s" />
 	      <field id="ahmf_3d"      long_name=" 3D      f-eddy viscosity coefficient"   unit="m2/s or m4/s"                           grid_ref="grid_T_3D"/>
+
+      <!-- product fields -->
+      <field_group id="diaprod">
+	<field id="ut"           long_name="product_of_sea_water_x_velocity_and_potential_temperature"      unit="degC m/s"      grid_ref="grid_U_3D"   />
+        <field id="ut_e3u"       long_name="product_of_sea_water_x_velocity_and_potential_temperature * e3u"  unit="degC m2/s"   grid_ref="grid_U_3D" > ut * e3u </field >
+	<field id="us"           long_name="product_of_sea_water_x_velocity_and_salinity"                   unit="PSU m/s"       grid_ref="grid_U_3D"   />
+        <field id="us_e3u"       long_name="product_of_sea_water_x_velocity_and_salinity * e3u"             unit="PSU m2/s"      grid_ref="grid_U_3D" > us * e3u </field >
+	<field id="urhop"        long_name="product_of_sea_water_x_velocity_and_potential_density"          unit="(kg/m3).(m/s)" grid_ref="grid_U_3D"   />
+        <field id="urhop_e3u"    long_name="product_of_sea_water_x_velocity_and_potential_density * e3u"    unit="(kg/m3).(m2/s)"   grid_ref="grid_U_3D" > urhop * e3u </field >
+	<field id="vt"           long_name="product_of_sea_water_y_velocity_and_potential_temperature"      unit="degC m/s"      grid_ref="grid_V_3D"   />
+        <field id="vt_e3v"       long_name="product_of_sea_water_y_velocity_and_potential_temperature * e3v"  unit="degC m2/s"   grid_ref="grid_V_3D" > vt * e3v </field >
+	<field id="vs"           long_name="product_of_sea_water_y_velocity_and_salinity"                   unit="PSU m/s"       grid_ref="grid_V_3D"   />
+        <field id="vs_e3v"       long_name="product_of_sea_water_y_velocity_and_salinity * e3t"             unit="PSU m2/s"      grid_ref="grid_V_3D" > vs * e3v </field >
+	<field id="vrhop"        long_name="product_of_sea_water_y_velocity_and_potential_density"          unit="(kg/m3).(m/s)" grid_ref="grid_V_3D"   />
+        <field id="vrhop_e3v"    long_name="product_of_sea_water_y_velocity_and_potential_density * e3t"    unit="(kg/m3).(m2/s)"  grid_ref="grid_V_3D" > vrhop * e3v </field >
+	<field id="wt"           long_name="product_of_upward_sea_water_velocity_and_potential_temperature" unit="degC m/s"      grid_ref="grid_W_3D"   />
+	<field id="ws"           long_name="product_of_upward_sea_water_velocity_and_salinity"              unit="PSU m/s"       grid_ref="grid_W_3D"   />
+	<field id="wrhop"        long_name="product_of_upward_sea_water_velocity_and_potential_density"     unit="(kg/m3).(m/s)" grid_ref="grid_W_3D"   />
+      </field_group>
 
       <!-- scalar variables available with key_diaar5 -->
@@ -608,5 +633,7 @@
       <!-- variables available with iceberg trajectories -->
 
-      <field_group id="icbvar" domain_ref="grid_T"  > 
+      <field_group id="icbvar" grid_ref="grid_T_2D"  > 
+        <field id="berg_total_melt"    long_name="icb melt rate 2 of icebergs"                     unit="kg/m2/s"                    />
+        <field id="berg_total_heat_flux"  long_name="icb latent heat of melting of icebergs"          unit="kg/m2/s"                    />
         <field id="berg_melt"          long_name="icb melt rate of icebergs"                       unit="kg/m2/s"                    />
         <field id="berg_buoy_melt"     long_name="icb buoyancy component of iceberg melt rate"     unit="kg/m2/s"                    />
Index: branches/UKMO/dev_r6501_GO6_package_trunk/NEMOGCM/CONFIG/SHARED/namelist_ref
===================================================================
--- branches/UKMO/dev_r6501_GO6_package_trunk/NEMOGCM/CONFIG/SHARED/namelist_ref	(revision 6507)
+++ branches/UKMO/dev_r6501_GO6_package_trunk/NEMOGCM/CONFIG/SHARED/namelist_ref	(revision 6825)
@@ -395,4 +395,9 @@
    ln_usecplmask = .false. !  use a coupling mask file to merge data received from several models
    !                       !   -> file cplmask.nc with the float variable called cplmask (jpi,jpj,nn_cplmodel)
+   ln_coupled_iceshelf_fluxes = .false. ! If true use rate of change of mass of Greenland and Antarctic icesheets to set the 
+                                        ! combined magnitude of the iceberg calving and iceshelf melting freshwater fluxes.
+   rn_greenland_calving_fraction = 0.5  ! Set fraction of total freshwater flux for iceberg calving - remainder goes to iceshelf melting.
+   rn_antarctica_calving_fraction = 0.5 ! Set fraction of total freshwater flux for iceberg calving - remainder goes to iceshelf melting.
+   rn_iceshelf_fluxes_tolerance = 1e-6  ! Fractional threshold for detecting differences in icesheet masses (must be positive definite).
 /
 !-----------------------------------------------------------------------
Index: branches/UKMO/dev_r6501_GO6_package_trunk/NEMOGCM/NEMO/OPA_SRC/DIA/diaprod.F90
===================================================================
--- branches/UKMO/dev_r6501_GO6_package_trunk/NEMOGCM/NEMO/OPA_SRC/DIA/diaprod.F90	(revision 6825)
+++ branches/UKMO/dev_r6501_GO6_package_trunk/NEMOGCM/NEMO/OPA_SRC/DIA/diaprod.F90	(revision 6825)
@@ -0,0 +1,228 @@
+MODULE diaprod
+! Requires key_iom_put
+# if defined key_iomput
+   !!======================================================================
+   !!                     ***  MODULE  diaprod  ***
+   !! Ocean diagnostics :  write ocean product diagnostics
+   !!=====================================================================
+   !! History :  3.4  ! 2012  (D. Storkey)  Original code
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   dia_prod      : calculate and write out product diagnostics
+   !!----------------------------------------------------------------------
+   USE oce             ! ocean dynamics and tracers 
+   USE dom_oce         ! ocean space and time domain
+   USE domvvl          ! for thickness weighted diagnostics if key_vvl
+   USE eosbn2          ! equation of state  (eos call)
+   USE phycst          ! physical constants
+   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
+   USE in_out_manager  ! I/O manager
+   USE iom
+   USE ioipsl
+   USE lib_mpp         ! MPP library
+   USE timing          ! preformance summary
+   USE wrk_nemo        ! working array
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   dia_prod                 ! routines called by step.F90
+
+   !! * Substitutions
+#  include "vectopt_loop_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 3.4 , NEMO Consortium (2012)
+   !! $Id $
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE dia_prod( kt )
+      !!---------------------------------------------------------------------
+      !!                  ***  ROUTINE dia_prod  ***
+      !!                   
+      !! ** Purpose :   Write out product diagnostics (uT, vS etc.)
+      !!
+      !! ** Method  :  use iom_put
+      !!               Product diagnostics are not thickness-weighted in this routine.
+      !!               They should be thickness-weighted using XIOS if key_vvl is set. 
+      !!----------------------------------------------------------------------
+      !!
+      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
+      !!
+      INTEGER                      ::   ji, jj, jk              ! dummy loop indices
+      REAL(wp)                     ::   zztmp, zztmpx, zztmpy   ! 
+      !!
+      REAL(wp), POINTER, DIMENSION(:,:)   :: z2d      ! 2D workspace
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: z3d      ! 3D workspace
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: zrhop    ! potential density
+      !!----------------------------------------------------------------------
+      ! 
+      IF( nn_timing == 1 )   CALL timing_start('dia_prod')
+      ! 
+      CALL wrk_alloc( jpi , jpj      , z2d )
+      CALL wrk_alloc( jpi , jpj, jpk , z3d )
+      CALL wrk_alloc( jpi , jpj, jpk , zrhop )
+      !
+
+      IF( iom_use("urhop") .OR. iom_use("vrhop") .OR. iom_use("wrhop") &
+#if ! defined key_diaar5
+     &  .OR. iom_use("rhop") &
+#endif
+     & ) THEN
+         CALL eos( tsn, z3d, zrhop )                 ! now in situ and potential density
+         zrhop(:,:,:) = zrhop(:,:,:)-1000.e0         ! reference potential density to 1000 to avoid precision issues in rhop2 calculation
+         zrhop(:,:,jpk) = 0._wp
+#if ! defined key_diaar5
+         CALL iom_put( 'rhop', zrhop )
+#else
+         ! If key_diaar5 set then there is already an iom_put call to output rhop.
+         ! Really should be a standard diagnostics option?
+#endif
+      ENDIF
+
+      IF( iom_use("ut") ) THEN
+         z3d(:,:,:) = 0.e0 
+         DO jk = 1, jpkm1
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  z3d(ji,jj,jk) = un(ji,jj,jk) * 0.5 * ( tsn(ji,jj,jk,jp_tem) + tsn(ji+1,jj,jk,jp_tem) )
+               END DO
+            END DO
+         END DO
+         CALL iom_put( "ut", z3d )                  ! product of temperature and zonal velocity at U points
+      ENDIF
+
+      IF( iom_use("vt") ) THEN
+         z3d(:,:,:) = 0.e0 
+         DO jk = 1, jpkm1
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  z3d(ji,jj,jk) = vn(ji,jj,jk) * 0.5 * ( tsn(ji,jj,jk,jp_tem) + tsn(ji,jj+1,jk,jp_tem) )
+               END DO
+            END DO
+         END DO
+         CALL iom_put( "vt", z3d )                  ! product of temperature and meridional velocity at V points
+      ENDIF
+
+      IF( iom_use("wt") ) THEN
+         z3d(:,:,:) = 0.e0 
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               z3d(ji,jj,1) = wn(ji,jj,1) * tsn(ji,jj,1,jp_tem)
+            END DO
+         END DO
+         DO jk = 2, jpkm1
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  z3d(ji,jj,jk) = wn(ji,jj,jk) * 0.5 * ( tsn(ji,jj,jk-1,jp_tem) + tsn(ji,jj,jk,jp_tem) )
+               END DO
+            END DO
+         END DO
+         CALL iom_put( "wt", z3d )                  ! product of temperature and vertical velocity at W points
+      ENDIF
+
+      IF( iom_use("us") ) THEN
+         z3d(:,:,:) = 0.e0 
+         DO jk = 1, jpkm1
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  z3d(ji,jj,jk) = un(ji,jj,jk) * 0.5 * ( tsn(ji,jj,jk,jp_sal) + tsn(ji+1,jj,jk,jp_sal) )
+               END DO
+            END DO
+         END DO
+         CALL iom_put( "us", z3d )                  ! product of salinity and zonal velocity at U points
+      ENDIF
+
+      IF( iom_use("vs") ) THEN
+         z3d(:,:,:) = 0.e0 
+         DO jk = 1, jpkm1
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  z3d(ji,jj,jk) = vn(ji,jj,jk) * 0.5 * ( tsn(ji,jj,jk,jp_sal) + tsn(ji,jj+1,jk,jp_sal) )
+               END DO
+            END DO
+         END DO
+         CALL iom_put( "vs", z3d )                  ! product of salinity and meridional velocity at V points
+      ENDIF
+
+      IF( iom_use("ws") ) THEN
+         z3d(:,:,:) = 0.e0 
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               z3d(ji,jj,1) = wn(ji,jj,1) * tsn(ji,jj,1,jp_sal)
+            END DO
+         END DO
+         DO jk = 2, jpkm1
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  z3d(ji,jj,jk) = wn(ji,jj,jk) * 0.5 * ( tsn(ji,jj,jk-1,jp_sal) + tsn(ji,jj,jk,jp_sal) )
+               END DO
+            END DO
+         END DO
+         CALL iom_put( "ws", z3d )                  ! product of salinity and vertical velocity at W points
+      ENDIF
+
+      IF( iom_use("urhop") ) THEN
+         z3d(:,:,:) = 0.e0 
+         DO jk = 1, jpkm1
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  z3d(ji,jj,jk) = un(ji,jj,jk) * 0.5 * ( zrhop(ji,jj,jk) + zrhop(ji+1,jj,jk) )
+               END DO
+            END DO
+         END DO
+         CALL iom_put( "urhop", z3d )                  ! product of density and zonal velocity at U points
+      ENDIF
+
+      IF( iom_use("vrhop") ) THEN
+         z3d(:,:,:) = 0.e0 
+         DO jk = 1, jpkm1
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  z3d(ji,jj,jk) = vn(ji,jj,jk) * 0.5 * ( zrhop(ji,jj,jk) + zrhop(ji,jj+1,jk) )
+               END DO
+            END DO
+         END DO
+         CALL iom_put( "vrhop", z3d )                  ! product of density and meridional velocity at V points
+      ENDIF
+
+      IF( iom_use("wrhop") ) THEN
+         z3d(:,:,:) = 0.e0 
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               z3d(ji,jj,1) = wn(ji,jj,1) * zrhop(ji,jj,1)
+            END DO
+         END DO
+         DO jk = 2, jpkm1
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  z3d(ji,jj,jk) = wn(ji,jj,jk) * 0.5 * ( zrhop(ji,jj,jk-1) + zrhop(ji,jj,jk) )
+               END DO
+            END DO
+         END DO
+         CALL iom_put( "wrhop", z3d )                  ! product of density and vertical velocity at W points
+      ENDIF
+
+      !
+      CALL wrk_dealloc( jpi , jpj      , z2d )
+      CALL wrk_dealloc( jpi , jpj, jpk , z3d )
+      CALL wrk_dealloc( jpi , jpj, jpk , zrhop )
+      !
+      IF( nn_timing == 1 )   CALL timing_stop('dia_prod')
+      !
+   END SUBROUTINE dia_prod
+#else
+   !!----------------------------------------------------------------------
+   !!   Default option :                                         NO diaprod
+   !!----------------------------------------------------------------------
+   LOGICAL, PUBLIC, PARAMETER :: lk_diaprod = .FALSE.   ! coupled flag
+CONTAINS
+   SUBROUTINE dia_prod( kt )   ! Empty routine
+      INTEGER ::   kt
+      WRITE(*,*) 'dia_prod: You should not have seen this print! error?', kt
+   END SUBROUTINE dia_prod
+#endif
+   !!======================================================================
+END MODULE diaprod
Index: branches/UKMO/dev_r6501_GO6_package_trunk/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90
===================================================================
--- branches/UKMO/dev_r6501_GO6_package_trunk/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90	(revision 6507)
+++ branches/UKMO/dev_r6501_GO6_package_trunk/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90	(revision 6825)
@@ -234,5 +234,4 @@
       IF( lk_zdftke ) THEN   
          CALL iom_put( "tke"      , en                               )    ! TKE budget: Turbulent Kinetic Energy   
-         CALL iom_put( "tke_niw"  , e_niw                            )    ! TKE budget: Near-inertial waves   
       ENDIF 
       CALL iom_put( "avs" , fsavs(:,:,:)               )    ! S vert. eddy diff. coef. (useful only with key_zdfddm)
@@ -982,8 +981,4 @@
       CALL histwrite( id_i, "sozotaux", kt, utau             , jpi*jpj    , idex )    ! i-wind stress
       CALL histwrite( id_i, "sometauy", kt, vtau             , jpi*jpj    , idex )    ! j-wind stress
-      IF( lk_vvl ) THEN
-         CALL histwrite( id_i, "vovvldep", kt, gdept_n(:,:,:), jpi*jpj*jpk, idex )!  T-cell depth       
-         CALL histwrite( id_i, "vovvle3t", kt, e3t_n (:,:,:), jpi*jpj*jpk, idex )!  T-cell thickness  
-      END IF
 
       IF(  .NOT.ln_linssh  ) THEN             
Index: branches/UKMO/dev_r6501_GO6_package_trunk/NEMOGCM/NEMO/OPA_SRC/ICB/icbrst.F90
===================================================================
--- branches/UKMO/dev_r6501_GO6_package_trunk/NEMOGCM/NEMO/OPA_SRC/ICB/icbrst.F90	(revision 6507)
+++ branches/UKMO/dev_r6501_GO6_package_trunk/NEMOGCM/NEMO/OPA_SRC/ICB/icbrst.F90	(revision 6825)
@@ -170,4 +170,12 @@
       cl_path = TRIM(cn_ocerst_outdir)
       IF( cl_path(LEN_TRIM(cl_path):) /= '/' ) cl_path = TRIM(cl_path) // '/'
+      IF ( ln_rstdate ) THEN
+         CALL ju2ymds( fjulday + rdt / rday, iyear, imonth, iday, zsec )           
+         WRITE(clkt, '(i4.4,2i2.2)') iyear, imonth, iday
+      ELSE
+         IF( kt > 999999999 ) THEN   ;   WRITE(clkt, *       ) kt
+         ELSE                        ;   WRITE(clkt, '(i8.8)') kt
+         ENDIF
+      ENDIF
       IF( lk_mpp ) THEN
          WRITE(cl_filename,'(A,"_icebergs_",A,"_restart_",I4.4,".nc")') TRIM(cexper), TRIM(ADJUSTL(clkt)), narea-1
Index: branches/UKMO/dev_r6501_GO6_package_trunk/NEMOGCM/NEMO/OPA_SRC/IOM/restart.F90
===================================================================
--- branches/UKMO/dev_r6501_GO6_package_trunk/NEMOGCM/NEMO/OPA_SRC/IOM/restart.F90	(revision 6507)
+++ branches/UKMO/dev_r6501_GO6_package_trunk/NEMOGCM/NEMO/OPA_SRC/IOM/restart.F90	(revision 6825)
@@ -27,4 +27,6 @@
    USE in_out_manager  ! I/O manager
    USE iom             ! I/O module
+   USE ioipsl, ONLY : ju2ymds    ! for calendar
+   USE sbc_oce         ! for icesheet freshwater input variables
    USE diurnal_bulk
    
@@ -87,5 +89,5 @@
          IF( nitrst <= nitend .AND. nitrst > 0 ) THEN 
             IF ( ln_rstdate ) THEN
-               CALL ju2ymds( fjulday + rdttra(1) / rday, iyear, imonth, iday, zsec )           
+               CALL ju2ymds( fjulday + rdt / rday, iyear, imonth, iday, zsec )           
                WRITE(clkt, '(i4.4,2i2.2)') iyear, imonth, iday
             ELSE
Index: branches/UKMO/dev_r6501_GO6_package_trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_ice.F90
===================================================================
--- branches/UKMO/dev_r6501_GO6_package_trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_ice.F90	(revision 6507)
+++ branches/UKMO/dev_r6501_GO6_package_trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_ice.F90	(revision 6825)
@@ -108,5 +108,6 @@
    ! variables used in the coupled interface
    INTEGER , PUBLIC, PARAMETER ::   jpl = ncat
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   u_ice, v_ice          ! jpi, jpj
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   u_ice, v_ice       ! jpi, jpj
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   a_p, ht_p          ! Meltpond fraction and depth
 #endif
    
Index: branches/UKMO/dev_r6501_GO6_package_trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90
===================================================================
--- branches/UKMO/dev_r6501_GO6_package_trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90	(revision 6507)
+++ branches/UKMO/dev_r6501_GO6_package_trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90	(revision 6825)
@@ -208,5 +208,5 @@
          CALL wrk_alloc( jpi,jpj,jpk, ztfrz3d ) 
          DO jk=1,jpk
-             CALL eos_fzp( tsn(:,:,jk,jp_sal), ztfrz3d(:,:,jk), fsdept_n(:,:,jk) )
+             CALL eos_fzp( tsn(:,:,jk,jp_sal), ztfrz3d(:,:,jk), gdept_n(:,:,jk) )
          ENDDO
          tsn(:,:,:,jp_tem) = MAX( tsn(:,:,:,jp_tem), ztfrz3d )
@@ -224,5 +224,5 @@
 
       ! calculate surface freezing temperature and send to CICE
-      CALL  eos_fzp(sss_m(:,:), sstfrz(:,:), fsdept_n(:,:,1)) 
+      CALL  eos_fzp(sss_m(:,:), sstfrz(:,:), gdept_n(:,:,1)) 
       CALL nemo2cice(sstfrz,Tf, 'T', 1. )
 
Index: branches/UKMO/dev_r6501_GO6_package_trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90
===================================================================
--- branches/UKMO/dev_r6501_GO6_package_trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90	(revision 6507)
+++ branches/UKMO/dev_r6501_GO6_package_trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90	(revision 6825)
@@ -38,4 +38,5 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   qisf                  !: net heat flux from ice shelf      [W/m2]
    REAL(wp), PUBLIC ::   rn_hisf_tbl                 !: thickness of top boundary layer [m]
+   LOGICAL , PUBLIC ::   ln_divisf                   !: flag to correct divergence 
    INTEGER , PUBLIC ::   nn_isf                      !: flag to choose between explicit/param/specified  
    INTEGER , PUBLIC ::   nn_isfblk                   !: flag to choose the bulk formulation to compute the ice shelf melting
@@ -92,4 +93,7 @@
       INTEGER               :: ji, jj               ! loop index
       REAL(wp), DIMENSION (:,:), POINTER :: zt_frz, zdep ! freezing temperature (zt_frz) at depth (zdep) 
+      REAL(wp)                     ::   zgreenland_fwfisf_sum, zantarctica_fwfisf_sum
+      INTEGER                      ::   ji, jj, jk, ijkmin, inum, ierror
+      REAL(wp)                     ::   zt_frz, zpress
       !!---------------------------------------------------------------------
       !
@@ -231,7 +235,4 @@
          CALL eos_fzp( stbl(:,:), zt_frz(:,:), zdep(:,:) )
          
-         ! if apply only on the trend and not as a volume flux (rdivisf = 0), fwfisf have to be set to 0 now
-         fwfisf(:,:) = rdivisf * fwfisf(:,:)         
- 
          risf_tsc(:,:,jp_tem) = qisf(:,:) * r1_rau0_rcp - fwfisf(:,:) * zt_frz(:,:) * r1_rau0 !
          risf_tsc(:,:,jp_sal) = 0.0_wp
Index: branches/UKMO/dev_r6501_GO6_package_trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90
===================================================================
--- branches/UKMO/dev_r6501_GO6_package_trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90	(revision 6507)
+++ branches/UKMO/dev_r6501_GO6_package_trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90	(revision 6825)
@@ -189,5 +189,4 @@
          fwfisf  (:,:)   = 0.0_wp ; fwfisf_b  (:,:)   = 0.0_wp
          risf_tsc(:,:,:) = 0.0_wp ; risf_tsc_b(:,:,:) = 0.0_wp
-         rdivisf       = 0.0_wp
       END IF
       IF( nn_ice == 0 .AND. nn_components /= jp_iam_opa )   fr_i(:,:) = 0._wp    ! no ice in the domain, ice fraction is always zero
Index: branches/UKMO/dev_r6501_GO6_package_trunk/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfmxl.F90
===================================================================
--- branches/UKMO/dev_r6501_GO6_package_trunk/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfmxl.F90	(revision 6507)
+++ branches/UKMO/dev_r6501_GO6_package_trunk/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfmxl.F90	(revision 6825)
@@ -28,4 +28,5 @@
 
    PUBLIC   zdf_mxl       ! called by step.F90
+   PUBLIC   zdf_mxl_alloc ! Used in zdf_tke_init
 
    INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   nmln    !: number of level in the mixed layer (used by TOP)
@@ -39,4 +40,10 @@
    REAL(wp), PUBLIC ::   rho_c = 0.01_wp    !: density criterion for mixed layer depth
    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 
 
    !!----------------------------------------------------------------------
Index: branches/UKMO/dev_r6501_GO6_package_trunk/NEMOGCM/NEMO/OPA_SRC/step.F90
===================================================================
--- branches/UKMO/dev_r6501_GO6_package_trunk/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 6507)
+++ branches/UKMO/dev_r6501_GO6_package_trunk/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 6825)
@@ -236,4 +236,5 @@
       IF( lk_diaar5  )   CALL dia_ar5( kstp )         ! ar5 diag
       IF( lk_diaharm )   CALL dia_harm( kstp )        ! Tidal harmonic analysis
+                         CALL dia_prod( kstp )
                          CALL dia_wri( kstp )         ! ocean model: outputs
       !
