Index: /branches/NERC/dev_r5589_marine_glacier_termini/NEMOGCM/CONFIG/SHARED/field_def.xml
===================================================================
--- /branches/NERC/dev_r5589_marine_glacier_termini/NEMOGCM/CONFIG/SHARED/field_def.xml	(revision 5604)
+++ /branches/NERC/dev_r5589_marine_glacier_termini/NEMOGCM/CONFIG/SHARED/field_def.xml	(revision 5605)
@@ -193,5 +193,11 @@
 
          <!-- * variable related to ice shelf forcing * -->
-         <field id="fwfisf"       long_name="Ice shelf melting"                             unit="kg/m2/s"  />
+         <field id="nktop"        long_name="top plumes level"                              unit="lev"      />
+         <field id="nksta"        long_name="sta plumes level"                              unit="lev"      />
+         <field id="nkovt"        long_name="ovt plumes level"                              unit="lev"      />
+         <field id="icwovt"       long_name="plumes ovt"                                    unit="m3/s"     />
+         <field id="qsubg"        long_name="subglacial melting"                            unit="m3/s"     />
+         <field id="qmelt"        long_name="ice wall melting"                              unit="m3/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"     />
          <field id="isfgammat"    long_name="transfert coefficient for isf (temperature)"   unit="m/s"      />
Index: /branches/NERC/dev_r5589_marine_glacier_termini/NEMOGCM/NEMO/OPA_SRC/BDY/bdyvol.F90
===================================================================
--- /branches/NERC/dev_r5589_marine_glacier_termini/NEMOGCM/NEMO/OPA_SRC/BDY/bdyvol.F90	(revision 5604)
+++ /branches/NERC/dev_r5589_marine_glacier_termini/NEMOGCM/NEMO/OPA_SRC/BDY/bdyvol.F90	(revision 5605)
@@ -24,4 +24,5 @@
    USE in_out_manager  ! I/O manager
    USE sbc_oce         ! ocean surface boundary conditions
+   USE traicw
 
    IMPLICIT NONE
@@ -91,5 +92,5 @@
       ! Calculate the cumulate surface Flux z_cflxemp (m3/s) over all the domain
       ! -----------------------------------------------------------------------
-      z_cflxemp = SUM ( ( emp(:,:)-rnf(:,:)+rdivisf*fwfisf(:,:) ) * bdytmask(:,:) * e1t(:,:) * e2t(:,:) ) / rau0
+      z_cflxemp = SUM ( ( emp(:,:)-rnf(:,:)+rdivisf*fwfisf(:,:) -qtot(:,:)) * bdytmask(:,:) * e1t(:,:) * e2t(:,:) ) / rau0
       IF( lk_mpp )   CALL mpp_sum( z_cflxemp )     ! sum over the global domain
 
Index: /branches/NERC/dev_r5589_marine_glacier_termini/NEMOGCM/NEMO/OPA_SRC/DOM/domhgr.F90
===================================================================
--- /branches/NERC/dev_r5589_marine_glacier_termini/NEMOGCM/NEMO/OPA_SRC/DOM/domhgr.F90	(revision 5604)
+++ /branches/NERC/dev_r5589_marine_glacier_termini/NEMOGCM/NEMO/OPA_SRC/DOM/domhgr.F90	(revision 5605)
@@ -318,5 +318,6 @@
          !                          ==========
          glam0 = 0.e0
-         gphi0 = - ppe2_m * 1.e-3
+         !gphi0 = - ppe2_m * 1.e-3  ! need to be restore before submission into the trunk
+         gphi0 = ppgphi0 - ppe2_m * 1.e-3
          
 #if defined key_agrif 
@@ -649,31 +650,4 @@
       CALL iom_close( inum )
       
-! need to be define for the extended grid south of -80S
-! some point are undefined but you need to have e1 and e2 .NE. 0
-      WHERE (e1t==0.0_wp)
-         e1t=1.0e2
-      END WHERE
-      WHERE (e1v==0.0_wp)
-         e1v=1.0e2
-      END WHERE
-      WHERE (e1u==0.0_wp)
-         e1u=1.0e2
-      END WHERE
-      WHERE (e1f==0.0_wp)
-         e1f=1.0e2
-      END WHERE
-      WHERE (e2t==0.0_wp)
-         e2t=1.0e2
-      END WHERE
-      WHERE (e2v==0.0_wp)
-         e2v=1.0e2
-      END WHERE
-      WHERE (e2u==0.0_wp)
-         e2u=1.0e2
-      END WHERE
-      WHERE (e2f==0.0_wp)
-         e2f=1.0e2
-      END WHERE
-       
     END SUBROUTINE hgr_read
     
Index: /branches/NERC/dev_r5589_marine_glacier_termini/NEMOGCM/NEMO/OPA_SRC/DYN/divcur.F90
===================================================================
--- /branches/NERC/dev_r5589_marine_glacier_termini/NEMOGCM/NEMO/OPA_SRC/DYN/divcur.F90	(revision 5604)
+++ /branches/NERC/dev_r5589_marine_glacier_termini/NEMOGCM/NEMO/OPA_SRC/DYN/divcur.F90	(revision 5605)
@@ -29,4 +29,5 @@
    USE sbcrnf          ! river runoff 
    USE sbcisf          ! ice shelf 
+   USE traicw          ! calving face
    USE cla             ! cross land advection             (cla_div routine)
    USE in_out_manager  ! I/O manager
@@ -230,4 +231,5 @@
       IF( ln_rnf      )   CALL sbc_rnf_div( hdivn )          ! runoffs   (update hdivn field)
       IF( ln_divisf .AND. (nn_isf /= 0) )   CALL sbc_isf_div( hdivn )          ! ice shelf (update hdivn field)
+      IF( ln_traicw   )   CALL div_icw( hdivn )          ! ice shelf (update hdivn field)
       IF( nn_cla == 1 )   CALL cla_div    ( kt )             ! Cross Land Advection (Update Hor. divergence)
       
@@ -329,4 +331,5 @@
       IF( ln_rnf      )   CALL sbc_rnf_div( hdivn )                            ! runoffs (update hdivn field)
       IF( ln_divisf .AND. (nn_isf .GT. 0) )   CALL sbc_isf_div( hdivn )          ! ice shelf (update hdivn field)
+      IF( ln_traicw   )   CALL div_icw( hdivn )          ! ice shelf (update hdivn field)
       IF( nn_cla == 1 )   CALL cla_div    ( kt )             ! Cross Land Advection (update hdivn field)
       !
Index: /branches/NERC/dev_r5589_marine_glacier_termini/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90
===================================================================
--- /branches/NERC/dev_r5589_marine_glacier_termini/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90	(revision 5604)
+++ /branches/NERC/dev_r5589_marine_glacier_termini/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90	(revision 5605)
@@ -23,4 +23,5 @@
    USE sbc_oce         ! surface boundary condition: ocean
    USE sbcisf          ! ice shelf variable (fwfisf)
+   USE traicw          ! ice shelf variable (fwfisf)
    USE dynspg_oce      ! surface pressure gradient variables
    USE phycst          ! physical constants
@@ -457,5 +458,6 @@
       ELSE
          zssh_frc(:,:) = zraur * z1_2 * (  emp(:,:) + emp_b(:,:) - rnf(:,:) - rnf_b(:,:)   &
-                &                        + rdivisf * ( fwfisf(:,:) + fwfisf_b(:,:) )       )
+                &                        + rdivisf * ( fwfisf(:,:) + fwfisf_b(:,:) )       &
+                                         - qtot(:,:) - qtot_b(:,:)                         )
       ENDIF
 #if defined key_asminc
Index: /branches/NERC/dev_r5589_marine_glacier_termini/NEMOGCM/NEMO/OPA_SRC/SBC/sbcfwb.F90
===================================================================
--- /branches/NERC/dev_r5589_marine_glacier_termini/NEMOGCM/NEMO/OPA_SRC/SBC/sbcfwb.F90	(revision 5604)
+++ /branches/NERC/dev_r5589_marine_glacier_termini/NEMOGCM/NEMO/OPA_SRC/SBC/sbcfwb.F90	(revision 5605)
@@ -28,4 +28,5 @@
    USE lbclnk          ! ocean lateral boundary conditions
    USE lib_fortran
+   USE traicw
 
    IMPLICIT NONE
@@ -108,5 +109,5 @@
          !
          IF( MOD( kt-1, kn_fsbc ) == 0 ) THEN
-            z_fwf = glob_sum( e1e2t(:,:) * ( emp(:,:) - rnf(:,:) + rdivisf * fwfisf(:,:) -  snwice_fmass(:,:) ) ) / area   ! sum over the global domain
+            z_fwf = glob_sum( e1e2t(:,:) * ( emp(:,:) - rnf(:,:) + rdivisf * fwfisf(:,:) -qtot(:,:) - snwice_fmass(:,:) ) ) / area   ! sum over the global domain
             zcoef = z_fwf * rcp
             emp(:,:) = emp(:,:) - z_fwf              * tmask(:,:,1)
Index: /branches/NERC/dev_r5589_marine_glacier_termini/NEMOGCM/NEMO/OPA_SRC/SBC/sbcrnf.F90
===================================================================
--- /branches/NERC/dev_r5589_marine_glacier_termini/NEMOGCM/NEMO/OPA_SRC/SBC/sbcrnf.F90	(revision 5604)
+++ /branches/NERC/dev_r5589_marine_glacier_termini/NEMOGCM/NEMO/OPA_SRC/SBC/sbcrnf.F90	(revision 5605)
@@ -364,4 +364,5 @@
                   END DO
                   nk_rnf(ji,jj) = jk
+               ELSEIF( h_rnf(ji,jj) == 0._wp    ) THEN   ;  nk_rnf(ji,jj) = 1 ! (need to be remove before sumition into the trunk)
                ELSEIF( h_rnf(ji,jj) == -1._wp   ) THEN   ;  nk_rnf(ji,jj) = 1
                ELSEIF( h_rnf(ji,jj) == -999._wp ) THEN   ;  nk_rnf(ji,jj) = mbkt(ji,jj)
Index: /branches/NERC/dev_r5589_marine_glacier_termini/NEMOGCM/NEMO/OPA_SRC/TRA/eosbn2.F90
===================================================================
--- /branches/NERC/dev_r5589_marine_glacier_termini/NEMOGCM/NEMO/OPA_SRC/TRA/eosbn2.F90	(revision 5604)
+++ /branches/NERC/dev_r5589_marine_glacier_termini/NEMOGCM/NEMO/OPA_SRC/TRA/eosbn2.F90	(revision 5605)
@@ -56,5 +56,5 @@
    !                   !! * Interface
    INTERFACE eos
-      MODULE PROCEDURE eos_insitu, eos_insitu_pot, eos_insitu_2d
+      MODULE PROCEDURE eos_insitu, eos_insitu_pot, eos_insitu_2d, eos_insitu_0d
    END INTERFACE
    !
@@ -558,4 +558,80 @@
    END SUBROUTINE eos_insitu_2d
 
+   SUBROUTINE eos_insitu_0d( pts, pdep, prd )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE eos_insitu_2d  ***
+      !!
+      !! ** Purpose :   Compute the in situ density (ratio rho/rau0) from
+      !!      potential temperature and salinity using an equation of state
+      !!      defined through the namelist parameter nn_eos. * 2D field case
+      !!
+      !! ** Action  : - prd , the in situ density (no units) (unmasked)
+      !!
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION( 2 ), INTENT(in   ) ::   pts   ! 1 : potential temperature  [Celcius]
+      !                                                           ! 2 : salinity               [psu]
+      REAL(wp), INTENT(in   ) ::   pdep  ! depth                      [m]
+      REAL(wp), INTENT(  out) ::   prd   ! in situ density
+      !
+      INTEGER  ::   ji, jj, jk                ! dummy loop indices
+      REAL(wp) ::   zt , zh , zs              ! local scalars
+      REAL(wp) ::   zn , zn0, zn1, zn2, zn3   !   -      -
+      !!----------------------------------------------------------------------
+      !
+      IF( nn_timing == 1 )   CALL timing_start('eos2d')
+      !
+      prd = 0._wp
+      !
+      SELECT CASE( nn_eos )
+      !
+      CASE( -1, 0 )                !==  polynomial TEOS-10 / EOS-80 ==!
+         !
+               !
+               zh  = pdep * r1_Z0                                  ! depth
+               zt  = pts (jp_tem) * r1_T0                           ! temperature
+               zs  = SQRT( ABS( pts(jp_sal) + rdeltaS ) * r1_S0 )   ! square root salinity
+               !
+               zn3 = EOS013*zt   &
+                  &   + EOS103*zs+EOS003
+                  !
+               zn2 = (EOS022*zt   &
+                  &   + EOS112*zs+EOS012)*zt   &
+                  &   + (EOS202*zs+EOS102)*zs+EOS002
+
+               zn1 = (((EOS041*zt   &
+                  &   + EOS131*zs+EOS031)*zt   &
+                  &   + (EOS221*zs+EOS121)*zs+EOS021)*zt   &
+                  &   + ((EOS311*zs+EOS211)*zs+EOS111)*zs+EOS011)*zt   &
+                  &   + (((EOS401*zs+EOS301)*zs+EOS201)*zs+EOS101)*zs+EOS001
+                  !
+               zn0 = (((((EOS060*zt   &
+                  &   + EOS150*zs+EOS050)*zt   &
+                  &   + (EOS240*zs+EOS140)*zs+EOS040)*zt   &
+                  &   + ((EOS330*zs+EOS230)*zs+EOS130)*zs+EOS030)*zt   &
+                  &   + (((EOS420*zs+EOS320)*zs+EOS220)*zs+EOS120)*zs+EOS020)*zt   &
+                  &   + ((((EOS510*zs+EOS410)*zs+EOS310)*zs+EOS210)*zs+EOS110)*zs+EOS010)*zt   &
+                  &   + (((((EOS600*zs+EOS500)*zs+EOS400)*zs+EOS300)*zs+EOS200)*zs+EOS100)*zs+EOS000
+                  !
+               zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0
+               !
+               prd = zn * r1_rau0 - 1._wp               ! unmasked in situ density anomaly
+               !
+      CASE( 1 )                !==  simplified EOS  ==!
+         !
+               zt    = pts  (jp_tem)  - 10._wp
+               zs    = pts  (jp_sal)  - 35._wp
+               zh    = pdep                          ! depth at the partial step level
+               !
+               zn =  - rn_a0 * ( 1._wp + 0.5_wp*rn_lambda1*zt + rn_mu1*zh ) * zt   &
+                  &  + rn_b0 * ( 1._wp - 0.5_wp*rn_lambda2*zs - rn_mu2*zh ) * zs   &
+                  &  - rn_nu * zt * zs
+                  !
+               prd = zn * r1_rau0               ! unmasked in situ density anomaly
+               !
+      END SELECT
+      !
+      IF( nn_timing == 1 )   CALL timing_stop('eos0d')
+      !
+   END SUBROUTINE eos_insitu_0d
 
    SUBROUTINE rab_3d( pts, pab )
@@ -921,5 +997,5 @@
                pn2(ji,jj,jk) = grav * (  zaw * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) )     &
                   &                    - zbw * ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) )  )  &
-                  &            / fse3w(ji,jj,jk) * tmask(ji,jj,jk)
+                  &            / fse3w(ji,jj,jk) * wmask(ji,jj,jk)
             END DO
          END DO
Index: /branches/NERC/dev_r5589_marine_glacier_termini/NEMOGCM/NEMO/OPA_SRC/TRA/tradmp.F90
===================================================================
--- /branches/NERC/dev_r5589_marine_glacier_termini/NEMOGCM/NEMO/OPA_SRC/TRA/tradmp.F90	(revision 5604)
+++ /branches/NERC/dev_r5589_marine_glacier_termini/NEMOGCM/NEMO/OPA_SRC/TRA/tradmp.F90	(revision 5605)
@@ -114,4 +114,11 @@
       !
       CASE( 0 )                   !==  newtonian damping throughout the water column  ==!
+          resto = 0.0_wp                              ! need to be restore before submissionto the trunk
+          !WHERE (un(3:11,3:7,:) .GE. 1.e-2)
+          !   resto(3:11,3:7,:) = 1._wp/1440._wp
+          !END WHERE
+          WHERE (un(12,3:7,:) .GE. 1.e-2)
+             resto(12,3:7,:) = 1._wp/1440._wp
+          END WHERE
          DO jk = 1, jpkm1
             DO jj = 2, jpjm1
@@ -240,7 +247,12 @@
 
          !Read in mask from file
-         CALL iom_open ( cn_resto, imask)
-         CALL iom_get  ( imask, jpdom_autoglo, 'resto', resto)
-         CALL iom_close( imask )
+! need to be restore before submission to the trunk
+!         CALL iom_open ( cn_resto, imask)
+!         CALL iom_get  ( imask, jpdom_autoglo, 'resto', resto)
+!         CALL iom_close( imask )
+          resto = 0.0_wp
+          WHERE (un(2:10,2,:) .LE. 1.e-3)
+             resto(2:10,2,:) = 1._wp/86400._wp
+          END WHERE
        ENDIF
 
Index: /branches/NERC/dev_r5589_marine_glacier_termini/NEMOGCM/NEMO/OPA_SRC/TRA/traicw.F90
===================================================================
--- /branches/NERC/dev_r5589_marine_glacier_termini/NEMOGCM/NEMO/OPA_SRC/TRA/traicw.F90	(revision 5605)
+++ /branches/NERC/dev_r5589_marine_glacier_termini/NEMOGCM/NEMO/OPA_SRC/TRA/traicw.F90	(revision 5605)
@@ -0,0 +1,663 @@
+MODULE traicw
+   !!==============================================================================
+   !!                       ***  MODULE  trasbc  ***
+   !! Ocean active tracers:  surface boundary condition
+   !!==============================================================================
+   !! History : 02/2015  P. Mathiot : original code
+   !!----------------------------------------------------------------------
+   !!----------------------------------------------------------------------
+   !!   tra_icw      : parametrisation of ice wall processes
+   !!----------------------------------------------------------------------
+   USE oce             ! ocean dynamics and active tracers
+   USE dom_oce         ! ocean space domain variables
+   USE in_out_manager  ! I/O manager
+   USE fldread         ! read input field at current time step
+   USE iom
+   USE wrk_nemo        ! Memory Allocation
+   USE timing          ! Timing
+   USE eosbn2
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   tra_icw    ! routine called by step.F90
+   PUBLIC   tra_icw_alloc ! routine call in sbcmod module
+   PUBLIC   tra_icw_init  ! (PUBLIC for TAM)
+   PUBLIC   div_icw   ! routine called in sshwzv module
+
+   LOGICAL, PUBLIC  :: ln_plumes
+   LOGICAL, PUBLIC  :: ln_ovt
+   LOGICAL, PUBLIC  :: ln_traicw
+
+   CHARACTER(len=100), PUBLIC ::   cn_dir          !: Root directory for location of icw files
+
+   TYPE(FLD_N)       , PUBLIC ::   sn_qmelt
+   TYPE(FLD_N)       , PUBLIC ::   sn_qsubg
+   TYPE(FLD_N)       , PUBLIC ::   sn_ovt
+   TYPE(FLD_N)       , PUBLIC ::   sn_ztop
+   TYPE(FLD_N)       , PUBLIC ::   sn_zovt
+   TYPE(FLD_N)       , PUBLIC ::   sn_icwmsk
+
+   REAL(wp) :: rn_rztop, rn_rzsta, rn_rzovt
+
+   TYPE(FLD), PUBLIC, ALLOCATABLE, DIMENSION(:) ::   sf_qmelt
+   TYPE(FLD), PUBLIC, ALLOCATABLE, DIMENSION(:) ::   sf_qsubg
+   TYPE(FLD), PUBLIC, ALLOCATABLE, DIMENSION(:) ::   sf_ovt
+   TYPE(FLD), PUBLIC, ALLOCATABLE, DIMENSION(:) ::   sf_ztop
+   TYPE(FLD), PUBLIC, ALLOCATABLE, DIMENSION(:) ::   sf_zovt
+   TYPE(FLD), PUBLIC, ALLOCATABLE, DIMENSION(:) ::   sf_icwmsk
+
+   INTEGER , PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:) :: nk_top, nk_ovt, nk_bot, nk_sta
+   REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:) :: h_neg, h_pos, h_sta
+   REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:) :: rztop, rzovt, rzsta
+   REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:) :: rovt 
+   REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:) :: qtot, qtot_b 
+   REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:) :: ricwmsk
+   REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:) :: hdivicw_b, hdivicw 
+   REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:,:) :: ricw_tsc_b, ricw_tsc
+
+   REAL(wp), PUBLIC, SAVE ::   lfusicw= 0.334e6_wp    ! phycst ?
+   REAL(wp), PUBLIC, SAVE ::   E0, K0, r1_E0K0 ! entrainment
+   REAL(wp), PUBLIC, SAVE ::   FPa, FPb, FPc   ! Freezing point
+   REAL(wp), PUBLIC, SAVE ::   GamT, GamTS0, zeropt0
+   REAL(wp), PUBLIC, SAVE ::   ricwpol1, ricwpol2, ricwpol3, ricwpol4 , ricwpol5 , ricwpol6, &
+   &                           ricwpol7, ricwpol8, ricwpol9, ricwpol10, ricwpol11              ! poly coef
+   !! * Substitutions
+#  include "domzgr_substitute.h90"
+#  include "vectopt_loop_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
+   !! $Id: trasbc.F90 4726 2014-07-23 16:27:21Z mathiot $
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   INTEGER FUNCTION tra_icw_alloc()
+      !!----------------------------------------------------------------------
+      !!                ***  ROUTINE tra_icw_alloc  ***
+      !!----------------------------------------------------------------------
+      ALLOCATE( h_neg  (jpi,jpj), h_pos  (jpi,jpj), h_sta (jpi,jpj), &
+         &      nk_top (jpi,jpj), nk_ovt (jpi,jpj), nk_bot(jpi,jpj), nk_sta(jpi,jpj),                       &
+         &      rovt   (jpi,jpj), qtot   (jpi,jpj), qtot_b(jpi,jpj),                                        & 
+         &      rztop  (jpi,jpj), rzovt  (jpi,jpj), rzsta (jpi,jpj),                                        &
+         &      ricwmsk(jpi,jpj), hdivicw(jpi,jpj,jpk), hdivicw_b(jpi,jpj,jpk),                             &
+         &      ricw_tsc_b(jpi,jpj,jpk,2), ricw_tsc (jpi,jpj,jpk,2), STAT=tra_icw_alloc )
+         !
+      IF( lk_mpp            )   CALL mpp_sum ( tra_icw_alloc )
+      IF( tra_icw_alloc > 0 )   CALL ctl_warn('tra_icw_alloc: allocation of arrays failed')
+   END FUNCTION tra_icw_alloc
+
+
+   SUBROUTINE tra_icw_init
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE tra_icw_init  ***
+      !!
+      !! ** Purpose :   Initialisation of parametrisation data
+      !!
+      !! ** Method  : - read the runoff namtra_icw namelist
+      !!
+      !! ** Action  : - read parameters
+      !!----------------------------------------------------------------------------------
+      CHARACTER(len=32) ::   rn_icw_file   ! runoff file name
+      INTEGER           ::   ios           ! Local integer output status for namelist read
+      INTEGER           ::   ierror,inum   ! temporary integer
+      !
+      NAMELIST/namtra_icw/ cn_dir  , ln_traicw, ln_plumes, ln_ovt ,                        &
+        &                  sn_qmelt, sn_qsubg , sn_ovt   , sn_ztop, sn_zovt, sn_icwmsk, rn_rztop, rn_rzovt, rn_rzsta
+      !!----------------------------------------------------------------------------------
+      !
+      !                                   ! ============
+      !                                   !   Namelist
+      !                                   ! ============
+      !
+      REWIND( numnam_ref )              ! Namelist namtra_icw in reference namelist : Runoffs 
+      READ  ( numnam_ref, namtra_icw, IOSTAT = ios, ERR = 901)
+901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtra_icw in reference namelist', lwp )
+
+      REWIND( numnam_cfg )              ! Namelist namtra_icw in configuration namelist : Runoffs
+      READ  ( numnam_cfg, namtra_icw, IOSTAT = ios, ERR = 902 )
+902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtra_icw in configuration namelist', lwp )
+      IF(lwm) WRITE ( numond, namtra_icw )
+      !
+      !                                         ! Control print
+      IF(lwp) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) 'tra_icw : parametrisation of melting along calving face '
+         WRITE(numout,*) '~~~~~~~ '
+         WRITE(numout,*) ' use plumes model poly   = ', ln_plumes
+         WRITE(numout,*) ' use prescription of ovt = ', ln_ovt
+      ENDIF
+      !
+      !                                         !==  allocate runoff arrays
+      IF( tra_icw_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'tra_icw_alloc : unable to allocate arrays' )
+      IF (ln_traicw) THEN
+
+        IF( .NOT. ln_plumes ) THEN           !* set up file structure
+           IF(lwp) WRITE(numout,*)
+           IF(lwp) WRITE(numout,*) '        icw parametrisation variables read in a file'
+           ALLOCATE( sf_qmelt(1), STAT=ierror )
+           IF( ierror > 0 ) THEN
+              CALL ctl_stop( 'tra_icw_init: unable to allocate sf_qmelt structure' )   ;   RETURN
+           ENDIF
+           ALLOCATE( sf_qmelt(1)%fnow(jpi,jpj,1)   )
+           IF( sn_qmelt%ln_tint )ALLOCATE( sf_qmelt(1)%fdta(jpi,jpj,1,2) )
+           !                                        ! fill sf_chl with sn_chl and control print
+           CALL fld_fill( sf_qmelt, (/ sn_qmelt /), cn_dir, 'tra_icw_init',   &
+              &                                         'ice wall parametrisation function of read qmelt', 'namtra_icw' )
+               !
+           ALLOCATE( sf_ovt(1), STAT=ierror )
+           IF( ierror > 0 ) THEN
+              CALL ctl_stop( 'tra_icw_init: unable to allocate sf_ovt structure' )   ;   RETURN
+           ENDIF
+           ALLOCATE( sf_ovt(1)%fnow(jpi,jpj,1)   )
+           IF( sn_ovt%ln_tint )ALLOCATE( sf_ovt(1)%fdta(jpi,jpj,1,2) )
+           !                                        ! fill sf_chl with sn_chl and control print
+           CALL fld_fill( sf_ovt, (/ sn_ovt /), cn_dir, 'tra_icw_init',   &
+              &                                         'ice wall parametrisation function of read ovt', 'namtra_icw' )
+
+           ALLOCATE( sf_ztop(1), STAT=ierror )
+           IF( ierror > 0 ) THEN
+              CALL ctl_stop( 'tra_icw_init: unable to allocate sf_ztop structure' )   ;   RETURN
+           ENDIF
+           ALLOCATE( sf_ztop(1)%fnow(jpi,jpj,1)   )
+           IF( sn_ztop%ln_tint )ALLOCATE( sf_ztop(1)%fdta(jpi,jpj,1,2) )
+           !                                        ! fill sf_chl with sn_chl and control print
+           CALL fld_fill( sf_ztop, (/ sn_ztop /), cn_dir, 'tra_icw_init',   &
+              &                                         'ice wall parametrisation function of read ztop', 'namtra_icw' )
+        
+           ALLOCATE( sf_zovt(1), STAT=ierror )
+           IF( ierror > 0 ) THEN
+              CALL ctl_stop( 'tra_icw_init: unable to allocate sf_zovt structure' )   ;   RETURN
+           ENDIF
+           ALLOCATE( sf_zovt(1)%fnow(jpi,jpj,1)   )
+           IF( sn_zovt%ln_tint )ALLOCATE( sf_zovt(1)%fdta(jpi,jpj,1,2) )
+           !                                        ! fill sf_chl with sn_chl and control print
+           CALL fld_fill( sf_zovt, (/ sn_zovt /), cn_dir, 'tra_icw_init',   &
+              &                                         'ice wall parametrisation function of read zovt', 'namtra_icw' )
+        END IF
+
+        ! set up mask
+        ALLOCATE( sf_icwmsk(1), STAT=ierror )
+        IF( ierror > 0 ) THEN
+           CALL ctl_stop( 'tra_icw_init: unable to allocate sf_icwmsk structure' )   ;   RETURN
+        ENDIF
+        ALLOCATE( sf_icwmsk(1)%fnow(jpi,jpj,1)   )
+        IF( sn_icwmsk%ln_tint )ALLOCATE( sf_icwmsk(1)%fdta(jpi,jpj,1,2) )
+                                                ! fill sf_chl with sn_chl and control print
+        CALL fld_fill( sf_icwmsk, (/ sn_icwmsk /), cn_dir, 'tra_icw_init',   &
+        &                                         'ice wall parametrisation function of read icwmsk', 'namtra_icw' )
+         
+        IF ( .NOT. ln_plumes ) THEN
+           ! read depth of the top of the plumes
+           rn_icw_file = TRIM( cn_dir )//TRIM( sn_ztop%clname )
+           CALL iom_open ( rn_icw_file, inum )                       
+           CALL iom_get  ( inum, jpdom_data, sn_ztop%clvar, rztop ) 
+           CALL iom_close( inum )                                   
+           ! read depth of the top of the plumes
+           rn_icw_file = TRIM( cn_dir )//TRIM( sn_zovt%clname )
+           CALL iom_open ( rn_icw_file, inum )                       
+           CALL iom_get  ( inum, jpdom_data, sn_zovt%clvar, rzovt ) 
+           CALL iom_close( inum )                                   
+        ELSE
+           !%melt_param_pierre
+!           E0  =  3.60e-2 ! Entrainment coefficient
+!           K0  =  1.00e-3 !
+!           FPa = -5.73e-2 ! Seawater freezing point slope
+!           FPb =  8.32e-2 ! Seawater freezing point offset
+!           FPc =  7.61e-4 ! Depth dependance of freezing point
+!           GamT   = 3.479e-4 !
+!           GamTS0 = 1.897e-4 !
+!           zeropt0 = 0.56    !
+!           ricwpol1 =  2.3109067e5/10._wp ! poly coef 1
+!           ricwpol2 = -6.2399669e5/ 9._wp ! poly coef 2
+!           ricwpol3 =  7.1371620e5/ 8._wp ! poly coef 3
+!           ricwpol4 = -4.5038638e5/ 7._wp ! poly coef 4
+!           ricwpol5 =  1.7122183e5/ 6._wp ! poly coef 5
+!           ricwpol6 = -0.4025004e5/ 5._wp ! poly coef 6
+!           ricwpol7 =  0.0580962e5/ 4._wp ! poly coef 7
+!           ricwpol8 = -0.0050814e5/ 3._wp ! poly coef 8
+!           ricwpol9 =  0.0002680e5/ 2._wp ! poly coef 9
+!           ricwpol10=  0.0000005e5/ 1._wp ! poly coef 10
+!           ricwpol11=  0.0
+           !%melt_param_pierre_plumes_lab
+!           E0  =  7.20e-2 ! Entrainment coefficient
+!           K0  =  2.50e-3 !
+!           FPa = -5.73e-2 ! Seawater freezing point slope
+!           FPb =  8.32e-2 ! Seawater freezing point offset
+!           FPc =  7.61e-4 ! Depth dependance of freezing point
+!           GamT   = 1.100e-3 !
+!           GamTS0 = 6.000e-4 !
+!           zeropt0 = 0.56    !
+!           ricwpol1 =  1.0412536e6/10._wp ! poly coef 1
+!           ricwpol2 = -2.8114553e6/ 9._wp ! poly coef 2
+!           ricwpol3 =  3.2153640e6/ 8._wp ! poly coef 3
+!           ricwpol4 = -2.0287589e6/ 7._wp ! poly coef 4
+!           ricwpol5 =  0.7711177e6/ 6._wp ! poly coef 5
+!           ricwpol6 = -0.1812180e6/ 5._wp ! poly coef 6
+!           ricwpol7 =  0.0261438e6/ 4._wp ! poly coef 7
+!           ricwpol8 = -0.0022843e6/ 3._wp ! poly coef 8
+!           ricwpol9 =  0.0001203e6/ 2._wp ! poly coef 9
+!           ricwpol10=  0.0000002e6/ 1._wp ! poly coef 10
+!           ricwpol11=  0.0
+           !%melt_param_pierre_plumes_lab
+           E0  =  3.60e-2 ! Entrainment coefficient
+           K0  =  2.50e-3 !
+           FPa = -5.73e-2 ! Seawater freezing point slope
+           FPb =  8.32e-2 ! Seawater freezing point offset
+           FPc =  7.61e-4 ! Depth dependance of freezing point
+           GamT   = 1.100e-3 !
+           GamTS0 = 6.000e-4 !
+           zeropt0 = 0.56    !
+           ricwpol1 =  0.7373182e6/10._wp ! poly coef 1
+           ricwpol2 = -1.9907329e6/ 9._wp ! poly coef 2
+           ricwpol3 =  2.2766835e6/ 8._wp ! poly coef 3
+           ricwpol4 = -1.4364516e6/ 7._wp ! poly coef 4
+           ricwpol5 =  0.5459653e6/ 6._wp ! poly coef 5
+           ricwpol6 = -0.1282985e6/ 5._wp ! poly coef 6
+           ricwpol7 =  0.0185074e6/ 4._wp ! poly coef 7
+           ricwpol8 = -0.0016168e6/ 3._wp ! poly coef 8
+           ricwpol9 =  0.0000851e6/ 2._wp ! poly coef 9
+           ricwpol10=  0.0000001e6/ 1._wp ! poly coef 10
+           ricwpol11=  0.0
+           r1_E0K0 =  1.0/(E0+K0)**(1./3.)
+        END IF
+        
+        ! define subglacial file variable
+        ALLOCATE( sf_qsubg(1), STAT=ierror )
+        IF( ierror > 0 ) THEN
+           CALL ctl_stop( 'tra_icw_init: unable to allocate sf_qsubg structure' )   ;   RETURN
+        ENDIF
+        ALLOCATE( sf_qsubg(1)%fnow(jpi,jpj,1)   )
+        IF( sn_qsubg%ln_tint )ALLOCATE( sf_qsubg(1)%fdta(jpi,jpj,1,2) )
+        !                                        ! fill sf_chl with sn_chl and control print
+        CALL fld_fill( sf_qsubg, (/ sn_qsubg /), cn_dir, 'tra_icw_init',   &
+           &                                         'ice wall parametrisation function of read qmelt', 'namtra_icw' )
+
+        ! read mask
+        rn_icw_file = TRIM( cn_dir )//TRIM( sn_icwmsk%clname )
+        CALL iom_open ( rn_icw_file, inum )                       
+        CALL iom_get  ( inum, jpdom_data, sn_icwmsk%clvar, ricwmsk ) 
+        CALL iom_close( inum )                                   
+
+      END IF
+
+      ! initialisation
+      nk_top(:,:)=jpk     ; nk_ovt(:,:)=jpk     ; nk_bot(:,:)=jpk ; nk_sta(:,:)=jpk
+!      rztop(:,:)=0.0      ; rzovt(:,:)=0.0     ; rzsta(:,:)=0.0 
+      qtot (:,:) = 0.0_wp ; rovt (:,:) = 0.0_wp ;  qtot_b(:,:)=0.0
+      ricw_tsc(:,:,:,:) = 0.0_wp ; ricw_tsc_b(:,:,:,:) = 0.0_wp
+      h_pos(:,:) = 1.0_wp ; h_neg(:,:) = 1.0_wp ; h_sta(:,:) = 1.0_wp
+      hdivicw(:,:,:) = 0.0_wp
+
+   END SUBROUTINE tra_icw_init
+
+   SUBROUTINE tra_icw ( kt )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE tra_icw  ***
+      !!
+      !! ** Purpose :   correct T/S
+      !!
+      !! ** Method  :   apply "negative runoff" at the bottom and positive runoff at the top
+      !!
+      !! ** Action  :   apply the heat and salt flux when we impose the overturning
+      !!                along glacier termini face
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT(in) ::   kt          ! ocean time step
+      !
+      INTEGER  ::   ji, jj, jk             ! dummy loop indices
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(:,:)  , POINTER       ::   ztfrz_melt, ztfrz_sub   ! freezing point used for temperature correction
+      REAL(wp), DIMENSION(:,:)  , POINTER       ::   qmelt , qsubg
+      REAL(wp), DIMENSION(:,:)  , POINTER       ::   zglDe, zglDUf, zTw, zSw, zRw
+      REAL(wp), DIMENSION(:,:)  , POINTER       ::   zdeldfp, zGamTS, zxscale, zmscale, zmfscale
+      REAL(wp), DIMENSION(:,:)  , POINTER       ::   zmeltgl, zfluxgl, zxscalegl
+      REAL(wp), DIMENSION(:,:,:), POINTER       ::   zmeltfluxres, zmf, zPf, zMelt
+      REAL(wp), DIMENSION(:,:,:), POINTER       ::   zricw_tsc_sum
+      REAL(wp), DIMENSION(2)                    ::   zTSw
+      REAL(wp) :: zpress, zmeltfres, zrespts
+
+      !
+      CALL wrk_alloc( jpi,jpj,   ztfrz_melt, ztfrz_sub)
+      CALL wrk_alloc( jpi,jpj,   qmelt , qsubg )
+      CALL wrk_alloc( jpi,jpj,2, zricw_tsc_sum ) 
+      CALL wrk_alloc( jpi,jpj, zglDe, zglDUf, zTw, zSw, zRw)
+      CALL wrk_alloc( jpi,jpj, zdeldfp, zGamTS, zxscale, zmscale, zmfscale )
+      CALL wrk_alloc( jpi,jpj, zmeltgl, zfluxgl, zxscalegl)
+      CALL wrk_alloc( jpi,jpj,jpk, zmeltfluxres, zmf, zPf, zMelt) 
+     
+      IF( kt /= nit000 ) THEN                      !          Swap of forcing fields          !
+         !                                         ! ---------------------------------------- !
+         ricw_tsc_b(:,:,:,:) = ricw_tsc(:,:,:,:)               ! where before fields are set at the end of the routine
+         qtot_b(:,:) = qtot(:,:)               ! where before fields are set at the end of the routine
+         ricw_tsc(:,:,:,:) = 0.0_wp
+      ENDIF
+ 
+
+      ! read subglacial
+      CALL fld_read ( kt, 1, sf_qsubg  )   
+      qsubg(:,:) = sf_qsubg(1)%fnow(:,:,1)
+
+      nk_top(:,:) = 1     ; nk_ovt(:,:) = 1 ; nk_bot(:,:) = mbkt(:,:)
+
+      IF ( .NOT. ln_plumes ) THEN
+         ! all is define in a file define in the init process
+         ! read ovt
+         CALL fld_read ( kt, 1, sf_ovt   )   
+         rovt(:,:)  = sf_ovt(1)%fnow(:,:,1)
+         ! read melting along the face
+         CALL fld_read ( kt, 1, sf_qmelt  )   
+         qmelt(:,:) = sf_qmelt(1)%fnow(:,:,1)
+         ! compute nk_pos and nk_ovt to force the overturning
+         WHERE (ricwmsk == 1)
+            rztop= rn_rztop
+            rzsta= rn_rzsta
+            rzovt= rn_rzovt
+         END WHERE
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               jk = 1
+               DO WHILE ( jk .LE. mbkt(ji,jj) .AND. fsdept_n(ji,jj,jk+1) < rztop(ji,jj) ) ;  jk = jk + 1 ;  END DO
+               nk_top(ji,jj) = jk
+            END DO
+         END DO
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               jk = 1
+               DO WHILE ( jk .LE. mbkt(ji,jj) .AND. fsdept_n(ji,jj,jk+1) < rzovt(ji,jj) ) ;  jk = jk + 1 ;  END DO
+               nk_ovt(ji,jj) = MAX(jk,nk_top(ji,jj)+1)   ! at least 1 level beneath the top outflow level (traicw_div issue)
+            END DO
+         END DO
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               jk = 1
+               DO WHILE ( jk .LE. mbkt(ji,jj) .AND. fsdept_n(ji,jj,jk+1) < rzsta(ji,jj) ) ;  jk = jk + 1 ;  END DO
+               nk_sta(ji,jj) = jk
+            END DO
+         END DO
+      ELSE
+      !----------------------------------------------------------------------%
+      ! Define physical constants and parameters (as defined within plume model).
+      !De   (:)   = fsdepw_n(ji,jj,:)                ! depth [m]
+         DO jj=1,jpj
+            DO ji=1,jpi
+               zglDe (ji,jj) = -fsdepw_n(ji,jj,nk_bot(ji,jj)+1)  ! Grounding line [m]
+               zglDUf(ji,jj) = qsubg(ji,jj) / e2t (ji,jj)        ! subglacial runoff [m2/s] e2t should be in the input file
+               zTw   (ji,jj) = tsn(ji,jj,nk_bot(ji,jj),1) 
+               zSw   (ji,jj) = tsn(ji,jj,nk_bot(ji,jj),2) 
+            END DO
+         END DO
+
+         !----------------------------------------------------------------------%
+         ! Define basic melting plume scales.
+         zdeldfp(:,:)  = (zTw(:,:) - FPa*zSw(:,:)-FPb-FPc*zglDe(:,:))/FPc;
+   
+         zGamTS (:,:)  = GamT*(0.545+3.5e-5*zdeldfp(:,:)*E0/(GamTS0+E0))
+         zxscale(:,:)  = zeropt0*(zGamTS(:,:)+E0)/(zeropt0*zGamTS(:,:)+E0)
+         zmscale(:,:)  = 3.16e-7*SQRT(1._wp/(K0+E0))*SQRT(zGamTS(:,:)/(zGamTS(:,:)+E0))    &
+         &              * (E0/(zGamTS(:,:)+E0))*((FPc*zdeldfp(:,:))**2.0)
+         zmfscale(:,:) = FPc*zdeldfp(:,:)/(FPc*zdeldfp(:,:)+84.3)*zGamTS(:,:)/(zGamTS(:,:)+E0)
+   
+         !----------------------------------------------------------------------%
+         ! Define scales for freshwater discharge.
+         zmeltgl  (:,:) = 7.9e-3*r1_E0K0*(zGamTS(:,:)*E0/(zGamTS(:,:)+E0))  &
+         &               * (zglDUf(:,:)**(1./3.))*(FPc*zdeldfp(:,:))/1.8_wp
+         zfluxgl  (:,:) = 0.64*E0*r1_E0K0*(zglDUf(:,:)**(1.0/3.0))*zdeldfp(:,:)/zxscale(:,:)/1.2_wp
+         zxscalegl(:,:) = zglDUf(:,:)/(zmeltgl(:,:)+zmscale(:,:));
+
+         !----------------------------------------------------------------------%
+         DO jj=1,jpj
+            DO ji=1,jpi
+               DO jk=1,nk_bot(ji,jj) !useless, if 3d ovt needed have to uncomment it
+         ! Calculate scaled meltwater flux.
+                  zrespts = (-fsdepw_n(ji,jj,jk)-zglDe(ji,jj))*zxscale(ji,jj)/zdeldfp(ji,jj)
+                  zmeltfres = 1._wp-zrespts
+                  zmeltfluxres(ji,jj,jk) =   ricwpol1 * zrespts ** 10._wp  &
+                  &                        + ricwpol2 * zrespts ** 9._wp   &
+                  &                        + ricwpol3 * zrespts ** 8._wp   &
+                  &                        + ricwpol4 * zrespts ** 7._wp   &
+                  &                        + ricwpol5 * zrespts ** 6._wp   &
+                  &                        + ricwpol6 * zrespts ** 5._wp   &
+                  &                        + ricwpol7 * zrespts ** 4._wp   &
+                  &                        + ricwpol8 * zrespts ** 3._wp   &
+                  &                        + ricwpol9 * zrespts ** 2._wp   &
+                  &                        + ricwpol10* zrespts ** 1._wp   &
+                  &                        + ricwpol11* zrespts ** 0._wp
+   
+         ! Derive physical solution
+                  zmf(ji,jj,jk) = zmeltfluxres(ji,jj,jk)*zmscale(ji,jj)*zdeldfp(ji,jj)/zxscale(ji,jj)
+                  zPf(ji,jj,jk) = zmf(ji,jj,jk)*(1._wp-zmeltfres*zmfscale(ji,jj))/(zmeltfres*zmfscale(ji,jj))                      &
+                  &               + zfluxgl(ji,jj)*zrespts
+                  zMelt(ji,jj,jk) = (  zmeltfluxres(ji,jj,jk)*(zmscale(ji,jj)+zmeltgl(ji,jj))                                      &
+                  &                  + zrespts*zmeltgl(ji,jj)*(zxscalegl(ji,jj)*(zdeldfp(ji,jj)+zglDe(ji,jj))/zdeldfp(ji,jj)**2.)) &
+                  &                 / (1._wp+(zxscalegl(ji,jj)*(zdeldfp(ji,jj)+zglDe(ji,jj))/zdeldfp(ji,jj)**2.0))                 &
+                  &                 * zdeldfp(ji,jj)/zxscale(ji,jj) + zglDUf(ji,jj);
+               END DO
+            END DO
+         END DO
+
+         ! initialise stabilisation level and max ovt level
+         nk_sta(:,:) = nk_bot(:,:)
+         nk_ovt(:,:) = nk_bot(:,:)
+         zRw(:,:) = -99.9_wp
+
+         ! compute nk_ovt
+         DO jj=2,jpj
+            DO ji=2,jpi
+
+               ! Calculate freezing temperature
+               qsubg(ji,jj) = zglDUf(ji,jj) * ricwmsk(ji,jj) * rau0 / e2t(ji,jj)
+               zpress = grav*rau0*fsdept(ji,jj,nk_bot(ji,jj))*1.e-04
+               CALL eos_fzp(0.0_wp, ztfrz_sub(ji,jj), zpress)  ! subglacial runoff is fresh water
+
+   
+               DO WHILE (zRw(ji,jj) .LE. rhd(ji,jj,nk_sta(ji,jj)) .AND. (ricwmsk(ji,jj) == 1) .AND. nk_sta(ji,jj) .GE. 1)
+   
+                  ! update nk_sta
+                  nk_sta(ji,jj) = nk_sta(ji,jj)-1
+
+                  ! compute nk_ovt
+                  jk = nk_sta(ji,jj)
+                  DO WHILE ( jk .LT. mbkt(ji,jj) .AND.                                                                    &
+                  &          fsdept_n(ji,jj,jk) < fsdepw(ji,jj,nk_sta(ji,jj)) + MAX(100._wp,fse3t(ji,jj,nk_sta(ji,jj))))
+                    jk = jk + 1
+                  END DO
+                  nk_ovt(ji,jj) = jk
+                  
+                  ! compute overshoot level
+                  jk = 1
+                  DO WHILE ( jk .LT. mbkt(ji,jj) .AND.                                                                    &
+                  &          fsdept_n(ji,jj,jk) < MAX(1._wp, fsdepw(ji,jj,nk_sta(ji,jj)) - 100._wp))
+                    jk = jk + 1
+                  END DO
+                  nk_top(ji,jj) = jk
+
+                  rovt (ji,jj) = zPf(ji,jj,jk) * rau0 / e2t(ji,jj)
+                  qmelt(ji,jj) = (zMelt(ji,jj,jk) - zglDUf(ji,jj)) * rau0 / e2t(ji,jj)
+                  qtot (ji,jj) =qmelt(ji,jj) + qsubg(ji,jj)
+
+                  ! Calculate freezing temperature
+                  zpress = grav*rau0*fsdept(ji,jj,nk_ovt(ji,jj))*1.e-04
+                  CALL eos_fzp(tsn(ji,jj,jk,jp_sal), ztfrz_melt(ji,jj), zpress)
+
+                  ! compute heat and salt content
+                  zTw(ji,jj) = SUM(tsn(ji,jj,jk:nk_bot(ji,jj),1) * fse3t(ji,jj,jk:nk_bot(ji,jj)))/SUM(fse3t(ji,jj,jk:nk_bot(ji,jj)))
+                  zTw(ji,jj) = ( ( rovt(ji,jj) * zTw(ji,jj) + qmelt(ji,jj) * ztfrz_melt(ji,jj) + qsubg(ji,jj) * ztfrz_sub(ji,jj) ) &
+                     &           - qmelt(ji,jj) * lfusicw * r1_rcp ) / ( rovt(ji,jj) + qtot(ji,jj))
+                  zSw(ji,jj) = SUM(tsn(ji,jj,jk:nk_bot(ji,jj),2) * fse3t(ji,jj,jk:nk_bot(ji,jj)))/SUM(fse3t(ji,jj,jk:nk_bot(ji,jj)))
+                  zSw(ji,jj) = ( rovt(ji,jj) * zSw(ji,jj) ) / ( rovt(ji,jj) + qtot(ji,jj) )
+                  zTSw(1)=zTw(ji,jj)
+                  zTSw(2)=zSw(ji,jj)
+                  CALL eos(zTSw(:), fsdept(ji,jj,nk_sta(ji,jj)), zRw(ji,jj)) ! density at nk_sta to test position
+
+               END DO
+
+               ! compute kn_sta
+               nk_sta(ji,jj) = nk_sta(ji,jj)+1
+
+               ! compute nk_ovt
+               jk = nk_sta(ji,jj)
+               DO WHILE ( jk .LE. mbkt(ji,jj) .AND.                                                                    &
+               &          fsdept_n(ji,jj,jk) < fsdepw(ji,jj,nk_sta(ji,jj)) + MAX(100._wp,fse3t(ji,jj,nk_sta(ji,jj))))
+                 jk = jk + 1
+               END DO
+               nk_ovt(ji,jj) = jk
+
+               ! compute overshoot level
+               jk = 1
+               DO WHILE ( jk .LT. mbkt(ji,jj) .AND.                                                                    &
+               &          fsdept_n(ji,jj,jk) < MAX(1._wp, fsdepw(ji,jj,nk_sta(ji,jj)) - 100._wp))
+                 jk = jk + 1
+               END DO
+               nk_top(ji,jj) = jk
+
+               ! compute melt and ovt
+               qmelt(ji,jj) = (zMelt(ji,jj,nk_top(ji,jj)) - zglDUf(ji,jj)) * e2t(ji,jj)
+               qsubg(ji,jj) =                               zglDUf(ji,jj)  * e2t(ji,jj)
+               rovt (ji,jj) = - zPf(ji,jj,nk_top(ji,jj)) * e2t(ji,jj)
+
+            END DO
+         END DO
+!----------------------------------------------------------------------%
+      END IF
+         
+      ! convert ovt (m3/s) in ovt(kg/m2/s)
+      ! rovt negatif
+      rovt (:,:)=rovt (:,:) * rau0 / e1e2t(:,:) * ricwmsk(:,:)
+      qmelt(:,:)=qmelt(:,:) * rau0 / e1e2t(:,:) * ricwmsk(:,:)
+      qsubg(:,:)=qsubg(:,:) * rau0 / e1e2t(:,:) * ricwmsk(:,:)
+      ! compute total fresh water into the system
+      qtot (:,:) = qmelt(:,:) + qsubg(:,:)
+
+      ! useful to lock the loop in case of no ovt parametrisation 
+      IF ( ln_ovt .OR. ln_plumes ) THEN
+      ELSE
+         rovt(:,:)  = 0.0_wp
+         nk_ovt(:,:)=nk_bot(:,:)+1
+      END IF
+
+      ! compute h_pos  and h_neg  to force the overturning
+!      nk_sta = nk_top
+      DO ji=1,jpi
+         DO jj=1,jpj
+            h_pos (ji,jj) = MAX(1._wp,SUM(fse3t_n(ji,jj,nk_top(ji,jj):nk_ovt(ji,jj)-1)))
+            h_neg (ji,jj) = MAX(1._wp,SUM(fse3t_n(ji,jj,nk_top(ji,jj):nk_bot(ji,jj)  )))
+            h_sta (ji,jj) = MAX(1._wp,SUM(fse3t_n(ji,jj,nk_sta(ji,jj):nk_ovt(ji,jj)-1)))
+         END DO
+      END DO
+
+      ! compute heat and salt flux due to the ovt (negative part)
+      ricw_tsc(:,:,:,:) = 0.0_wp
+      DO ji=1,jpi
+         DO jj=1,jpj
+            DO jk=nk_top(ji,jj),nk_bot(ji,jj)
+               ricw_tsc(ji,jj,jk,1) = rovt(ji,jj) * tsn(ji,jj,jk,1) * r1_rau0 / h_neg (ji,jj) * ricwmsk(ji,jj)
+               ricw_tsc(ji,jj,jk,2) = rovt(ji,jj) * tsn(ji,jj,jk,2) * r1_rau0 / h_neg (ji,jj) * ricwmsk(ji,jj)
+            END DO
+         END DO
+      END DO
+
+      ! compute total heat and salt for negative part
+      zricw_tsc_sum(:,:,1) = SUM(ricw_tsc(:,:,:,1) * fse3t_n(:,:,:),3) 
+      zricw_tsc_sum(:,:,2) = SUM(ricw_tsc(:,:,:,2) * fse3t_n(:,:,:),3)
+     
+      ! compute heat and salt flux due to the ovt (positive part)
+      DO jj=1,jpj
+         DO ji=1,jpi
+            ! Calculate freezing temperature
+            zpress = grav*rau0*fsdept(ji,jj,nk_bot(ji,jj))*1.e-04
+            CALL eos_fzp(0.0_wp, ztfrz_sub(ji,jj), zpress)  ! subglacial runoff is fresh water
+
+            DO jk=nk_sta(ji,jj),nk_ovt(ji,jj) - 1
+               ! Calculate freezing temperature
+               zpress = grav*rau0*fsdept(ji,jj,jk)*1.e-04
+               CALL eos_fzp(tsn(ji,jj,jk,jp_sal), ztfrz_melt(ji,jj), zpress)
+
+               ! compute heat and salt content
+               ricw_tsc(ji,jj,jk,1) = ricw_tsc(ji,jj,jk,1) + ( - zricw_tsc_sum(ji,jj,1)                                    &
+                  &                    + ( qmelt(ji,jj) * ztfrz_melt(ji,jj) + qsubg(ji,jj) * ztfrz_sub(ji,jj) ) * r1_rau0  &
+                  &                    - qmelt(ji,jj) * lfusicw * r1_rau0_rcp                           )                  & 
+                  &                  / h_sta(ji,jj) * ricwmsk(ji,jj)
+               ricw_tsc(ji,jj,jk,2) = ricw_tsc(ji,jj,jk,2) + (- zricw_tsc_sum(ji,jj,2) + qtot(ji,jj) * 0.0_wp * r1_rau0        )  &
+                  &                  / h_sta(ji,jj) * ricwmsk(ji,jj)
+            END DO
+         END DO
+      END DO
+
+      IF(lwp .AND. kt == 1) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) 'tra_icw : parametrisation of melting along calving face '
+         WRITE(numout,*) '~~~~~~~ '
+         WRITE(numout,*) ' rzsta, nk_sta, h_sta        = ', MAXVAL(rzsta*ricwmsk), MAXVAL(nk_sta*ricwmsk), MAXVAL(h_sta*ricwmsk)
+         WRITE(numout,*) ' rzovt, nk_ovt, h_ovt        = ', MAXVAL(rzovt*ricwmsk), MAXVAL(nk_ovt*ricwmsk), MAXVAL(h_neg*ricwmsk)
+      ENDIF
+
+      IF( kt == nit000 ) THEN                          !   set the forcing field at nit000 - 1    !
+         !                                             ! ---------------------------------------- !
+     ! need to do something for restart before submission into the trunk
+     !    IF( ln_rstart .AND.    &                               !* Restart: read in restart file
+     !       & iom_varid( numror, 'rnf_b', ldstop = .FALSE. ) > 0 ) THEN
+     !       IF(lwp) WRITE(numout,*) '          nit000-1 runoff forcing fields red in the restart file'
+     !       CALL iom_get( numror, jpdom_autoglo, 'rnf_b', rnf_b )     ! before runoff
+     !       CALL iom_get( numror, jpdom_autoglo, 'rnf_hc_b', rnf_tsc_b(:,:,jp_tem) )   ! before heat content of runoff
+     !       CALL iom_get( numror, jpdom_autoglo, 'rnf_sc_b', rnf_tsc_b(:,:,jp_sal) )   ! before salinity content of runoff
+     !    ELSE                                                   !* no restart: set from nit000 values
+     !       IF(lwp) WRITE(numout,*) '          nit000-1 runoff forcing fields set to nit000'
+             ricw_tsc_b(:,:,:,:) = ricw_tsc(:,:,:,:)
+             qtot_b(:,:) = qtot(:,:)
+     !    ENDIF
+      ENDIF
+      
+      ! output
+      CALL iom_put('qmelt',  qmelt * r1_rau0 * e1e2t(:,:) )
+      CALL iom_put('qsubg',  qsubg * r1_rau0 * e1e2t(:,:) )
+      CALL iom_put('icwovt', rovt * r1_rau0 * e1e2t(:,:) )
+      CALL iom_put('nktop',  FLOAT(nk_top(:,:)) )
+      CALL iom_put('nksta',  FLOAT(nk_sta(:,:)) )
+      CALL iom_put('nkovt',  FLOAT(nk_ovt(:,:)) )
+
+      ! deallocate array
+      CALL wrk_dealloc( jpi,jpj,   ztfrz_melt, ztfrz_sub )
+      CALL wrk_dealloc( jpi,jpj,   qmelt , qsubg )
+      CALL wrk_dealloc( jpi,jpj,2, zricw_tsc_sum ) 
+      CALL wrk_dealloc( jpi,jpj, zglDe, zglDUf, zTw, zSw, zRw)
+      CALL wrk_dealloc( jpi,jpj, zdeldfp, zGamTS, zxscale, zmscale, zmfscale )
+      CALL wrk_dealloc( jpi,jpj, zmeltgl, zfluxgl, zxscalegl)
+      CALL wrk_dealloc( jpi,jpj,jpk, zmeltfluxres, zmf, zPf, zMelt) 
+   END SUBROUTINE tra_icw
+
+   SUBROUTINE div_icw( phdivn )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE div_icw  ***
+      !!
+      !! ** Purpose :   update the horizontal divergence with the ovt inflow/outflow
+      !!
+      !! ** Method  :
+      !!                CAUTION : icw contribution is positive (inflow) decreasing the
+      !!                          divergence and expressed in m/s
+      !!
+      !! ** Action  :   phdivn   decreased by the runoff inflow
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   phdivn   ! horizontal divergence
+      !!
+      INTEGER  ::   ji, jj, jk   ! dummy loop indices
+      !!----------------------------------------------------------------------
+      !
+      hdivicw_b(:,:,:) = hdivicw(:,:,:)
+      hdivicw(:,:,:) = 0.0_wp
+      DO jj = 1, jpj                   ! update the depth over which runoffs are distributed
+         DO ji = 1, jpi
+            DO jk = nk_top(ji,jj), nk_bot(ji,jj)                  
+               hdivicw(ji,jj,jk) = hdivicw(ji,jj,jk) -     rovt(ji,jj)                 * r1_rau0 / h_neg(ji,jj)
+            END DO
+            DO jk = nk_sta(ji,jj), nk_ovt(ji,jj)-1                    
+               hdivicw(ji,jj,jk) = hdivicw(ji,jj,jk) - ( - rovt(ji,jj) + qtot(ji,jj) ) * r1_rau0 / h_sta(ji,jj)
+            END DO
+         END DO
+      END DO
+      !
+      phdivn(:,:,:) = phdivn(:,:,:) + 0.5 * ( hdivicw(:,:,:) + hdivicw_b(:,:,:) )
+      !
+   END SUBROUTINE div_icw
+      
+END MODULE traicw
+
Index: /branches/NERC/dev_r5589_marine_glacier_termini/NEMOGCM/NEMO/OPA_SRC/TRA/trasbc.F90
===================================================================
--- /branches/NERC/dev_r5589_marine_glacier_termini/NEMOGCM/NEMO/OPA_SRC/TRA/trasbc.F90	(revision 5604)
+++ /branches/NERC/dev_r5589_marine_glacier_termini/NEMOGCM/NEMO/OPA_SRC/TRA/trasbc.F90	(revision 5605)
@@ -26,4 +26,5 @@
    USE trdtra          ! trends manager: tracers 
    !
+   USE traicw          ! ice wall parametrization
    USE in_out_manager  ! I/O manager
    USE prtctl          ! Print control
@@ -287,4 +288,11 @@
          END DO  
       ENDIF
+
+
+      ! compute trend
+      IF ( ln_traicw ) THEN
+         tsa(:,:,:,1) = tsa(:,:,:,1) + 0.5_wp * ( ricw_tsc(:,:,:,1) + ricw_tsc_b(:,:,:,1) )
+         tsa(:,:,:,2) = tsa(:,:,:,2) + 0.5_wp * ( ricw_tsc(:,:,:,2) + ricw_tsc_b(:,:,:,2) )
+      END IF
  
       IF( l_trdtra )   THEN                      ! send trends for further diagnostics
Index: /branches/NERC/dev_r5589_marine_glacier_termini/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90
===================================================================
--- /branches/NERC/dev_r5589_marine_glacier_termini/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90	(revision 5604)
+++ /branches/NERC/dev_r5589_marine_glacier_termini/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90	(revision 5605)
@@ -422,4 +422,5 @@
       !                                     ! Active tracers
                             CALL tra_qsr_init   ! penetrative solar radiation qsr
+                            CALL tra_icw_init   ! penetrative solar radiation qsr
                             CALL tra_bbc_init   ! bottom heat flux
       IF( lk_trabbl     )   CALL tra_bbl_init   ! advective (and/or diffusive) bottom boundary layer scheme
Index: /branches/NERC/dev_r5589_marine_glacier_termini/NEMOGCM/NEMO/OPA_SRC/step.F90
===================================================================
--- /branches/NERC/dev_r5589_marine_glacier_termini/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 5604)
+++ /branches/NERC/dev_r5589_marine_glacier_termini/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 5605)
@@ -120,4 +120,6 @@
                          CALL bn2    ( tsb, rab_b, rn2b ) ! before Brunt-Vaisala frequency
                          CALL bn2    ( tsn, rab_n, rn2  ) ! now    Brunt-Vaisala frequency
+      IF ( ln_traicw  )  CALL tra_icw( kstp )             ! name need to be changed the best way is a 3d runoff routine
+                                                          ! it should simplified this, sbcisf, sbciscpl, sbcrnf and trasbc
       !
       !  VERTICAL PHYSICS
Index: /branches/NERC/dev_r5589_marine_glacier_termini/NEMOGCM/NEMO/OPA_SRC/step_oce.F90
===================================================================
--- /branches/NERC/dev_r5589_marine_glacier_termini/NEMOGCM/NEMO/OPA_SRC/step_oce.F90	(revision 5604)
+++ /branches/NERC/dev_r5589_marine_glacier_termini/NEMOGCM/NEMO/OPA_SRC/step_oce.F90	(revision 5605)
@@ -30,4 +30,5 @@
 
    USE traqsr           ! solar radiation penetration      (tra_qsr routine)
+   USE traicw           ! ice wall melting parametrisation (tra_icw routine)
    USE trasbc           ! surface boundary condition       (tra_sbc routine)
    USE trabbc           ! bottom boundary condition        (tra_bbc routine)
