Index: /branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90
===================================================================
--- /branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90	(revision 9189)
+++ /branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90	(revision 9190)
@@ -284,5 +284,5 @@
       IF(lwp) THEN
          WRITE(numout,*)
-         WRITE(numout,*) 'dom_nam  : domain initialization through namelist read'
+         WRITE(numout,*) 'dom_nam : domain initialization through namelist read'
          WRITE(numout,*) '~~~~~~~ '
       ENDIF
@@ -297,5 +297,5 @@
       !
       IF(lwp) THEN                  ! control print
-         WRITE(numout,*) '   Namelist : namrun'
+         WRITE(numout,*) '   Namelist : namrun   ---   run parameters'
          WRITE(numout,*) '      job number                      nn_no           = ', nn_no
          WRITE(numout,*) '      experiment name for output      cn_exp          = ', TRIM( cn_exp           )
@@ -359,14 +359,15 @@
       IF( Agrif_Root() ) THEN
 #endif
+      IF(lwp) WRITE(numout,*)
       SELECT CASE ( nleapy )        ! Choose calendar for IOIPSL
       CASE (  1 ) 
          CALL ioconf_calendar('gregorian')
-         IF(lwp) WRITE(numout,*) '   The IOIPSL calendar is "gregorian", i.e. leap year'
+         IF(lwp) WRITE(numout,*) '   ==>>>   The IOIPSL calendar is "gregorian", i.e. leap year'
       CASE (  0 )
          CALL ioconf_calendar('noleap')
-         IF(lwp) WRITE(numout,*) '   The IOIPSL calendar is "noleap", i.e. no leap year'
+         IF(lwp) WRITE(numout,*) '   ==>>>   The IOIPSL calendar is "noleap", i.e. no leap year'
       CASE ( 30 )
          CALL ioconf_calendar('360d')
-         IF(lwp) WRITE(numout,*) '   The IOIPSL calendar is "360d", i.e. 360 days in a year'
+         IF(lwp) WRITE(numout,*) '   ==>>>   The IOIPSL calendar is "360d", i.e. 360 days in a year'
       END SELECT
 #if defined key_agrif
@@ -512,7 +513,7 @@
          CALL iom_get( inum, 'ORCA_index', zorca_res )   ;   kk_cfg = INT( zorca_res )
          !
-         WRITE(ldtxt(ii),*) '       '                                                     ;   ii = ii+1
-         WRITE(ldtxt(ii),*) '       ==>>>   ORCA configuration '                          ;   ii = ii+1
-         WRITE(ldtxt(ii),*) '       '                                                     ;   ii = ii+1
+         WRITE(ldtxt(ii),*) '   .'                                                     ;   ii = ii+1
+         WRITE(ldtxt(ii),*) '   ==>>>   ORCA configuration '                           ;   ii = ii+1
+         WRITE(ldtxt(ii),*) '   .'                                                     ;   ii = ii+1
          !
       ELSE                                !- cd_cfg & k_cfg are not used
@@ -524,6 +525,6 @@
             CALL iom_getatt( inum, 'cn_cfg', cd_cfg )  ! returns   !  if not found
             CALL iom_getatt( inum, 'nn_cfg', kk_cfg )  ! returns -999 if not found
-            IF( TRIM(cd_cfg) .EQ. '!') cd_cfg = 'UNKNOWN'
-            IF( kk_cfg .EQ. -999     ) kk_cfg = -9999999
+            IF( TRIM(cd_cfg) == '!') cd_cfg = 'UNKNOWN'
+            IF( kk_cfg == -999     ) kk_cfg = -9999999
          ENDIF
          !
@@ -536,9 +537,9 @@
       CALL iom_close( inum )
       !
-      WRITE(ldtxt(ii),*) '   cn_cfg = ', TRIM(cd_cfg), '   nn_cfg = ', kk_cfg             ;   ii = ii+1
-      WRITE(ldtxt(ii),*) '   jpiglo = ', kpi                                              ;   ii = ii+1
-      WRITE(ldtxt(ii),*) '   jpjglo = ', kpj                                              ;   ii = ii+1
-      WRITE(ldtxt(ii),*) '   jpkglo = ', kpk                                              ;   ii = ii+1
-      WRITE(ldtxt(ii),*) '   type of global domain lateral boundary   jperio = ', kperio  ;   ii = ii+1
+      WRITE(ldtxt(ii),*) '      cn_cfg = ', TRIM(cd_cfg), '   nn_cfg = ', kk_cfg             ;   ii = ii+1
+      WRITE(ldtxt(ii),*) '      jpiglo = ', kpi                                              ;   ii = ii+1
+      WRITE(ldtxt(ii),*) '      jpjglo = ', kpj                                              ;   ii = ii+1
+      WRITE(ldtxt(ii),*) '      jpkglo = ', kpk                                              ;   ii = ii+1
+      WRITE(ldtxt(ii),*) '      type of global domain lateral boundary   jperio = ', kperio  ;   ii = ii+1
       !        
    END SUBROUTINE domain_cfg
Index: /branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DOM/domhgr.F90
===================================================================
--- /branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DOM/domhgr.F90	(revision 9189)
+++ /branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DOM/domhgr.F90	(revision 9190)
@@ -183,6 +183,6 @@
       IF(lwp) THEN
          WRITE(numout,*)
-         WRITE(numout,*) 'hgr_read : read the horizontal coordinates in mesh_mask'
-         WRITE(numout,*) '~~~~~~~~      jpiglo = ', jpiglo, ' jpjglo = ', jpjglo, ' jpk = ', jpk
+         WRITE(numout,*) '   hgr_read : read the horizontal coordinates in mesh_mask'
+         WRITE(numout,*) '   ~~~~~~~~      jpiglo = ', jpiglo, ' jpjglo = ', jpjglo, ' jpk = ', jpk
       ENDIF
       !
Index: /branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90
===================================================================
--- /branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90	(revision 9189)
+++ /branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90	(revision 9190)
@@ -1017,5 +1017,5 @@
       !
 #if defined key_agrif
-      IF( (.NOT.Agrif_Root()).AND.(.NOT.ln_vvl_zstar) )CALL ctl_stop( 'AGRIF is implemented with zstar coordinate only' )
+      IF( (.NOT.Agrif_Root()).AND.(.NOT.ln_vvl_zstar) )   CALL ctl_stop( 'AGRIF is implemented with zstar coordinate only' )
 #endif
       !
Index: /branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DYN/dynadv.F90
===================================================================
--- /branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DYN/dynadv.F90	(revision 9189)
+++ /branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DYN/dynadv.F90	(revision 9190)
@@ -131,10 +131,10 @@
          WRITE(numout,*)
          SELECT CASE( n_dynadv )
-         CASE( np_LIN_dyn )   ;   WRITE(numout,*) '      ===>>   linear dynamics : no momentum advection used'
-         CASE( np_VEC_c2  )   ;   WRITE(numout,*) '      ===>>   vector form : keg + zad + vor is used' 
+         CASE( np_LIN_dyn )   ;   WRITE(numout,*) '   ==>>>   linear dynamics : no momentum advection used'
+         CASE( np_VEC_c2  )   ;   WRITE(numout,*) '   ==>>>   vector form : keg + zad + vor is used' 
             IF( nn_dynkeg == nkeg_C2  )   WRITE(numout,*) '              with Centered standard keg scheme'
             IF( nn_dynkeg == nkeg_HW  )   WRITE(numout,*) '              with Hollingsworth keg scheme'
-         CASE( np_FLX_c2  )   ;   WRITE(numout,*) '      ===>>   flux form   : 2nd order scheme is used'
-         CASE( np_FLX_ubs )   ;   WRITE(numout,*) '      ===>>   flux form   : UBS       scheme is used'
+         CASE( np_FLX_c2  )   ;   WRITE(numout,*) '   ==>>>   flux form   : 2nd order scheme is used'
+         CASE( np_FLX_ubs )   ;   WRITE(numout,*) '   ==>>>   flux form   : UBS       scheme is used'
          END SELECT
       ENDIF
Index: /branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf.F90
===================================================================
--- /branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf.F90	(revision 9189)
+++ /branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf.F90	(revision 9190)
@@ -181,8 +181,8 @@
       IF(lwp) THEN
          WRITE(numout,*)
-         IF( nldf == np_no_ldf )   WRITE(numout,*) '      ===>>   NO lateral viscosity'
-         IF( nldf == np_lap    )   WRITE(numout,*) '      ===>>   iso-level laplacian operator'
-         IF( nldf == np_lap_i  )   WRITE(numout,*) '      ===>>   rotated laplacian operator with iso-level background'
-         IF( nldf == np_blp    )   WRITE(numout,*) '      ===>>   iso-level bi-laplacian operator'
+         IF( nldf == np_no_ldf )   WRITE(numout,*) '   ==>>>   NO lateral viscosity'
+         IF( nldf == np_lap    )   WRITE(numout,*) '   ==>>>   iso-level laplacian operator'
+         IF( nldf == np_lap_i  )   WRITE(numout,*) '   ==>>>   rotated laplacian operator with iso-level background'
+         IF( nldf == np_blp    )   WRITE(numout,*) '   ==>>>   iso-level bi-laplacian operator'
       ENDIF
       !
Index: /branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg.F90
===================================================================
--- /branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg.F90	(revision 9189)
+++ /branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg.F90	(revision 9190)
@@ -226,7 +226,7 @@
       IF(lwp) THEN
          WRITE(numout,*)
-         IF( nspg == np_EXP )   WRITE(numout,*) '      ===>>   explicit free surface'
-         IF( nspg == np_TS  )   WRITE(numout,*) '      ===>>   free surface with time splitting scheme'
-         IF( nspg == np_NO  )   WRITE(numout,*) '      ===>>   No surface surface pressure gradient trend in momentum Eqs.'
+         IF( nspg == np_EXP )   WRITE(numout,*) '   ==>>>   explicit free surface'
+         IF( nspg == np_TS  )   WRITE(numout,*) '   ==>>>   free surface with time splitting scheme'
+         IF( nspg == np_NO  )   WRITE(numout,*) '   ==>>>   No surface surface pressure gradient trend in momentum Eqs.'
       ENDIF
       !
Index: /branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DYN/dynvor.F90
===================================================================
--- /branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DYN/dynvor.F90	(revision 9189)
+++ /branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DYN/dynvor.F90	(revision 9190)
@@ -611,13 +611,13 @@
       SELECT CASE( n_dynadv )
       CASE( np_LIN_dyn )
-         IF(lwp) WRITE(numout,*) '      ===>>   linear dynamics : total vorticity = Coriolis'
+         IF(lwp) WRITE(numout,*) '   ==>>>   linear dynamics : total vorticity = Coriolis'
          nrvm = np_COR        ! planetary vorticity
          ntot = np_COR        !     -         -
       CASE( np_VEC_c2  )
-         IF(lwp) WRITE(numout,*) '      ===>>   vector form dynamics : total vorticity = Coriolis + relative vorticity' 
+         IF(lwp) WRITE(numout,*) '   ==>>>   vector form dynamics : total vorticity = Coriolis + relative vorticity' 
          nrvm = np_RVO        ! relative vorticity
          ntot = np_CRV        ! relative + planetary vorticity         
       CASE( np_FLX_c2 , np_FLX_ubs  )
-         IF(lwp) WRITE(numout,*) '      ===>>   flux form dynamics : total vorticity = Coriolis + metric term'
+         IF(lwp) WRITE(numout,*) '   ==>>>   flux form dynamics : total vorticity = Coriolis + metric term'
          nrvm = np_MET        ! metric term
          ntot = np_CME        ! Coriolis + metric term
@@ -627,8 +627,8 @@
          WRITE(numout,*)
          SELECT CASE( nvor_scheme )
-         CASE( np_ENE )   ;   WRITE(numout,*) '      ===>>   energy conserving scheme'
-         CASE( np_ENS )   ;   WRITE(numout,*) '      ===>>   enstrophy conserving scheme'
-         CASE( np_MIX )   ;   WRITE(numout,*) '      ===>>   mixed enstrophy/energy conserving scheme'
-         CASE( np_EEN )   ;   WRITE(numout,*) '      ===>>   energy and enstrophy conserving scheme'
+         CASE( np_ENE )   ;   WRITE(numout,*) '   ==>>>   energy conserving scheme'
+         CASE( np_ENS )   ;   WRITE(numout,*) '   ==>>>   enstrophy conserving scheme'
+         CASE( np_MIX )   ;   WRITE(numout,*) '   ==>>>   mixed enstrophy/energy conserving scheme'
+         CASE( np_EEN )   ;   WRITE(numout,*) '   ==>>>   energy and enstrophy conserving scheme'
          END SELECT         
       ENDIF
Index: /branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/ICB/icb_oce.F90
===================================================================
--- /branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/ICB/icb_oce.F90	(revision 9189)
+++ /branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/ICB/icb_oce.F90	(revision 9190)
@@ -1,4 +1,3 @@
 MODULE icb_oce
-
    !!======================================================================
    !!                       ***  MODULE  icb_oce  ***
Index: /branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/ICB/icbclv.F90
===================================================================
--- /branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/ICB/icbclv.F90	(revision 9189)
+++ /branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/ICB/icbclv.F90	(revision 9190)
@@ -1,4 +1,3 @@
 MODULE icbclv
-
    !!======================================================================
    !!                       ***  MODULE  icbclv  ***
@@ -11,4 +10,5 @@
    !!            -    !  2011-05  (Alderson)       budgets into separate module
    !!----------------------------------------------------------------------
+
    !!----------------------------------------------------------------------
    !!   icb_clv_flx   : transfer input flux of ice into iceberg classes
@@ -45,10 +45,10 @@
       !!
       !!----------------------------------------------------------------------
-      INTEGER, INTENT(in)             :: kt
+      INTEGER, INTENT(in) ::   kt
       !
-      REAL(wp)                        :: zcalving_used, zdist, zfact
-      INTEGER                         :: jn, ji, jj                    ! loop counters
-      INTEGER                         :: imx                           ! temporary integer for max berg class
-      LOGICAL, SAVE                   :: ll_first_call = .TRUE.
+      REAL(wp)      ::   zcalving_used, zdist, zfact
+      INTEGER       ::   jn, ji, jj                    ! loop counters
+      INTEGER       ::   imx                           ! temporary integer for max berg class
+      LOGICAL, SAVE ::   ll_first_call = .TRUE.
       !!----------------------------------------------------------------------
       !
@@ -70,8 +70,8 @@
          DO jj = 2, jpjm1
             DO ji = 2, jpim1
-               IF( berg_grid%calving(ji,jj) /= 0._wp )                                  &    ! Need units of J
-                  berg_grid%stored_heat(ji,jj) = SUM( berg_grid%stored_ice(ji,jj,:) ) *         &  ! initial stored ice in kg
-                                         berg_grid%calving_hflx(ji,jj) * e1e2t(ji,jj) /   &  ! J/s/m2 x m^2 = J/s
-                                         berg_grid%calving(ji,jj)                            ! /calving in kg/s
+               IF( berg_grid%calving(ji,jj) /= 0._wp )                                          &    ! Need units of J
+                  berg_grid%stored_heat(ji,jj) = SUM( berg_grid%stored_ice(ji,jj,:) ) *         &    ! initial stored ice in kg
+                     &                   berg_grid%calving_hflx(ji,jj) * e1e2t(ji,jj) / berg_grid%calving(ji,jj)   ! J/s/m2 x m^2 
+                     !                                                                                             ! = J/s/calving in kg/s
             END DO
          END DO
@@ -80,11 +80,11 @@
       ! assume that all calving flux must be distributed even if distribution array does not sum
       ! to one - this may not be what is intended, but it's what you've got
-      DO jj = 1,jpj
-         DO ji = 1,jpi
+      DO jj = 1, jpj
+         DO ji = 1, jpi
             imx = berg_grid%maxclass(ji,jj)
             zdist = SUM( rn_distribution(1:nclasses) ) / SUM( rn_distribution(1:imx) )
             DO jn = 1, imx
-               berg_grid%stored_ice(ji,jj,jn) = berg_grid%stored_ice(ji,jj,jn) + &
-                                          berg_dt * berg_grid%calving(ji,jj) * rn_distribution(jn) * zdist
+               berg_grid%stored_ice(ji,jj,jn) = berg_grid%stored_ice(ji,jj,jn)     &
+                  &                           + berg_dt * berg_grid%calving(ji,jj) * rn_distribution(jn) * zdist
             END DO
          END DO
@@ -98,4 +98,5 @@
       !
    END SUBROUTINE icb_clv_flx
+
 
    SUBROUTINE icb_clv()
@@ -171,5 +172,5 @@
       END DO
       !
-      DO jn = 1,nclasses
+      DO jn = 1, nclasses
          CALL lbc_lnk( berg_grid%stored_ice(:,:,jn), 'T', 1._wp )
       END DO
Index: /branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/ICB/icbdia.F90
===================================================================
--- /branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/ICB/icbdia.F90	(revision 9189)
+++ /branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/ICB/icbdia.F90	(revision 9190)
@@ -1,4 +1,3 @@
 MODULE icbdia
-
    !!======================================================================
    !!                       ***  MODULE  icbdia  ***
@@ -13,6 +12,20 @@
    !!            -  !                            from the right points in the code
    !!----------------------------------------------------------------------
+ 
    !!----------------------------------------------------------------------
-   !! icb_dia_init : initialise iceberg budgeting
+   !!   icb_dia_init  : initialise iceberg budgeting
+   !!   icb_dia       : global iceberg diagnostics
+   !!   icb_dia_step  : reset at the beginning of each timestep
+   !!   icb_dia_put   : output (via iom_put) iceberg fields
+   !!   icb_dia_calve : 
+   !!   icb_dia_income: 
+   !!   icb_dia_size  : 
+   !!   icb_dia_speed : 
+   !!   icb_dia_melt  : 
+   !!   report_state  : 
+   !!   report_consistant : 
+   !!   report_budget : 
+   !!   report_istate : 
+   !!   report_ibudget: 
    !!----------------------------------------------------------------------
    USE par_oce        ! ocean parameters
@@ -85,5 +98,5 @@
       !!----------------------------------------------------------------------
       !
-      IF( .NOT. ln_bergdia ) RETURN
+      IF( .NOT.ln_bergdia )   RETURN
 
       ALLOCATE( berg_melt    (jpi,jpj)   )           ;   berg_melt   (:,:)   = 0._wp
@@ -136,6 +149,6 @@
 
       floating_mass_start       = icb_utl_mass( first_berg )
-      bergs_mass_start          = icb_utl_mass( first_berg, justbergs=.true. )
-      bits_mass_start           = icb_utl_mass( first_berg, justbits=.true. )
+      bergs_mass_start          = icb_utl_mass( first_berg, justbergs=.TRUE. )
+      bits_mass_start           = icb_utl_mass( first_berg, justbits =.TRUE. )
       IF( lk_mpp ) THEN
          ALLOCATE( rsumbuf(23) )          ; rsumbuf(:) = 0._wp
@@ -159,11 +172,11 @@
       !! for this we pack variables into buffer so we only need one mpp_sum
       !!----------------------------------------------------------------------
-      LOGICAL, INTENT(in) ::   ld_budge
-      !
-      INTEGER             ::   ik
-      REAL(wp)            ::   zunused_calving, ztmpsum, zgrdd_berg_mass, zgrdd_bits_mass
-      !!----------------------------------------------------------------------
-      !
-      IF( .NOT. ln_bergdia )   RETURN
+      LOGICAL, INTENT(in) ::   ld_budge   !
+      !
+      INTEGER ::   ik
+      REAL(wp)::   zunused_calving, ztmpsum, zgrdd_berg_mass, zgrdd_bits_mass
+      !!----------------------------------------------------------------------
+      !
+      IF( .NOT.ln_bergdia )   RETURN
 
       zunused_calving      = SUM( berg_grid%calving(:,:) )
@@ -181,17 +194,17 @@
       ztmpsum              = SUM( berg_grid%calving_hflx(:,:) * e1e2t(:,:) * tmask_i(:,:) )
       calving_out_heat_net = calving_out_heat_net + ztmpsum * berg_dt   ! Units of J
-
+      !
       IF( ld_budge ) THEN
          stored_end        = SUM( berg_grid%stored_ice(:,:,:) )
          stored_heat_end   = SUM( berg_grid%stored_heat(:,:) )
          floating_mass_end = icb_utl_mass( first_berg )
-         bergs_mass_end    = icb_utl_mass( first_berg,justbergs=.true. )
-         bits_mass_end     = icb_utl_mass( first_berg,justbits=.true. )
+         bergs_mass_end    = icb_utl_mass( first_berg,justbergs=.TRUE. )
+         bits_mass_end     = icb_utl_mass( first_berg,justbits =.TRUE. )
          floating_heat_end = icb_utl_heat( first_berg )
-
+         !
          nbergs_end        = icb_utl_count()
          zgrdd_berg_mass   = SUM( berg_mass(:,:)*e1e2t(:,:)*tmask_i(:,:) )
          zgrdd_bits_mass   = SUM( bits_mass(:,:)*e1e2t(:,:)*tmask_i(:,:) )
-
+         !
          IF( lk_mpp ) THEN
             rsumbuf( 1) = stored_end
@@ -218,7 +231,7 @@
             rsumbuf(22) = zgrdd_berg_mass
             rsumbuf(23) = zgrdd_bits_mass
-
+            !
             CALL mpp_sum( rsumbuf(1:23), 23)
-
+            !
             stored_end                = rsumbuf( 1)
             stored_heat_end           = rsumbuf( 2)
@@ -244,5 +257,5 @@
             zgrdd_berg_mass           = rsumbuf(22)
             zgrdd_bits_mass           = rsumbuf(23)
-
+            !
             nsumbuf(1) = nbergs_end
             nsumbuf(2) = nbergs_calved
@@ -252,7 +265,6 @@
                nsumbuf(4+ik) = nbergs_calved_by_class(ik)
             END DO
-
             CALL mpp_sum( nsumbuf(1:nclasses+4), nclasses+4 )
-
+            !
             nbergs_end        = nsumbuf(1)
             nbergs_calved     = nsumbuf(2)
@@ -261,57 +273,57 @@
             DO ik = 1,nclasses
                nbergs_calved_by_class(ik)= nsumbuf(4+ik)
-            ENDDO
-
+            END DO
+            !
          ENDIF
-
-         CALL report_state( 'stored ice','kg','',stored_start,'',stored_end,'')
-         CALL report_state( 'floating','kg','',floating_mass_start,'',floating_mass_end,'',nbergs_end)
-         CALL report_state( 'icebergs','kg','',bergs_mass_start,'',bergs_mass_end,'')
-         CALL report_state( 'bits','kg','',bits_mass_start,'',bits_mass_end,'')
-         CALL report_istate( 'berg #','',nbergs_start,'',nbergs_end,'')
+         !
+         CALL report_state  ( 'stored ice','kg','',stored_start,'',stored_end,'')
+         CALL report_state  ( 'floating','kg','',floating_mass_start,'',floating_mass_end,'',nbergs_end )
+         CALL report_state  ( 'icebergs','kg','',bergs_mass_start,'',bergs_mass_end,'')
+         CALL report_state  ( 'bits','kg','',bits_mass_start,'',bits_mass_end,'')
+         CALL report_istate ( 'berg #','',nbergs_start,'',nbergs_end,'')
          CALL report_ibudget( 'berg #','calved',nbergs_calved, &
-                                       'melted',nbergs_melted, &
-                                       '#',nbergs_start,nbergs_end)
+            &                          'melted',nbergs_melted, &
+            &                          '#',nbergs_start,nbergs_end)
          CALL report_budget( 'stored mass','kg','calving used',calving_used_net, &
-                                           'bergs',calving_to_bergs_net, &
-                                           'stored mass',stored_start,stored_end)
+            &                              'bergs',calving_to_bergs_net, &
+            &                              'stored mass',stored_start,stored_end)
          CALL report_budget( 'floating mass','kg','calving used',calving_to_bergs_net, &
-                                             'bergs',melt_net, &
-                                             'stored mass',floating_mass_start,floating_mass_end)
+            &                                'bergs',melt_net, &
+            &                                'stored mass',floating_mass_start,floating_mass_end)
          CALL report_budget( 'berg mass','kg','calving',calving_to_bergs_net, &
-                                              'melt+eros',berg_melt_net, &
-                                              'berg mass',bergs_mass_start,bergs_mass_end)
+            &                                 'melt+eros',berg_melt_net, &
+            &                                 'berg mass',bergs_mass_start,bergs_mass_end)
          CALL report_budget( 'bits mass','kg','eros used',bits_src_net, &
-                                              'bergs',bits_melt_net, &
-                                              'stored mass',bits_mass_start,bits_mass_end)
+            &                                 'bergs',bits_melt_net, &
+            &                                 'stored mass',bits_mass_start,bits_mass_end)
          CALL report_budget( 'net mass','kg','recvd',calving_rcv_net, &
-                                             'rtrnd',calving_ret_net, &
-                                             'net mass',stored_start+floating_mass_start, &
-                                                        stored_end+floating_mass_end)
+            &                                'rtrnd',calving_ret_net, &
+            &                                'net mass',stored_start+floating_mass_start, &
+            &                                           stored_end+floating_mass_end)
          CALL report_consistant( 'iceberg mass','kg','gridded',zgrdd_berg_mass,'bergs',bergs_mass_end)
          CALL report_consistant( 'bits mass','kg','gridded',zgrdd_bits_mass,'bits',bits_mass_end)
          CALL report_state( 'net heat','J','',stored_heat_start+floating_heat_start,'', &
-                                              stored_heat_end+floating_heat_end,'')
+            &                                 stored_heat_end+floating_heat_end,'')
          CALL report_state( 'stored heat','J','',stored_heat_start,'',stored_heat_end,'')
          CALL report_state( 'floating heat','J','',floating_heat_start,'',floating_heat_end,'')
          CALL report_budget( 'net heat','J','net heat',calving_src_heat_net, &
-                                            'net heat',calving_out_heat_net, &
-                                            'net heat',stored_heat_start+floating_heat_start, &
-                                                       stored_heat_end+floating_heat_end)
+            &                               'net heat',calving_out_heat_net, &
+            &                               'net heat',stored_heat_start+floating_heat_start, &
+            &                                          stored_heat_end+floating_heat_end)
          CALL report_budget( 'stored heat','J','calving used',calving_src_heat_used_net, &
-                                               'bergs',heat_to_bergs_net, &
-                                               'net heat',stored_heat_start,stored_heat_end)
+            &                                  'bergs',heat_to_bergs_net, &
+            &                                  'net heat',stored_heat_start,stored_heat_end)
          CALL report_budget( 'flting heat','J','calved',heat_to_bergs_net, &
-                                               'melt',heat_to_ocean_net, &
-                                               'net heat',floating_heat_start,floating_heat_end)
+            &                                  'melt',heat_to_ocean_net, &
+            &                                  'net heat',floating_heat_start,floating_heat_end)
          IF (nn_verbose_level >= 1) THEN
             CALL report_consistant( 'top interface','kg','from SIS',calving_src_net, &
-                                    'received',calving_rcv_net)
+               &                    'received',calving_rcv_net)
             CALL report_consistant( 'bot interface','kg','sent',calving_out_net, &
-                                    'returned',calving_ret_net)
+               &                    'returned',calving_ret_net)
          ENDIF
          WRITE( numicb, '("calved by class = ",i6,20(",",i6))') (nbergs_calved_by_class(ik),ik=1,nclasses)
-         IF ( nspeeding_tickets > 0 ) WRITE( numicb, '("speeding tickets issued = ",i6)') nspeeding_tickets
-
+         IF( nspeeding_tickets > 0 )   WRITE( numicb, '("speeding tickets issued = ",i6)') nspeeding_tickets
+         !
          nbergs_start              = nbergs_end
          stored_start              = stored_end
@@ -350,15 +362,15 @@
       !!----------------------------------------------------------------------
       !
-      IF( .NOT. ln_bergdia ) RETURN
-      berg_melt    (:,:)   = 0._wp
-      buoy_melt    (:,:)   = 0._wp
-      eros_melt    (:,:)   = 0._wp
-      conv_melt    (:,:)   = 0._wp
-      bits_src     (:,:)   = 0._wp
-      bits_melt    (:,:)   = 0._wp
-      bits_mass    (:,:)   = 0._wp
-      berg_mass    (:,:)   = 0._wp
-      virtual_area (:,:)   = 0._wp
-      real_calving (:,:,:) = 0._wp
+      IF( .NOT.ln_bergdia )   RETURN
+      berg_melt   (:,:)   = 0._wp
+      buoy_melt   (:,:)   = 0._wp
+      eros_melt   (:,:)   = 0._wp
+      conv_melt   (:,:)   = 0._wp
+      bits_src    (:,:)   = 0._wp
+      bits_melt   (:,:)   = 0._wp
+      bits_mass   (:,:)   = 0._wp
+      berg_mass   (:,:)   = 0._wp
+      virtual_area(:,:)   = 0._wp
+      real_calving(:,:,:) = 0._wp
       !
    END SUBROUTINE icb_dia_step
@@ -369,5 +381,5 @@
       !!----------------------------------------------------------------------
       !
-      IF( .NOT. ln_bergdia )   RETURN            !!gm useless iom will control whether it is output or not
+      IF( .NOT.ln_bergdia )   RETURN            !!gm useless iom will control whether it is output or not
       !
       CALL iom_put( "berg_melt"        , berg_melt   (:,:)   )   ! Melt rate of icebergs                     [kg/m2/s]
@@ -388,5 +400,5 @@
       !!----------------------------------------------------------------------
       !!----------------------------------------------------------------------
-      INTEGER,  INTENT(in)  ::   ki, kj, kn
+      INTEGER , INTENT(in)  ::   ki, kj, kn
       REAL(wp), INTENT(in)  ::   pcalved
       REAL(wp), INTENT(in)  ::   pheated
@@ -411,5 +423,5 @@
       !!----------------------------------------------------------------------
       !
-      IF( .NOT. ln_bergdia ) RETURN
+      IF( .NOT.ln_bergdia )   RETURN
       !
       IF( kt == nit000 ) THEN
@@ -437,9 +449,9 @@
       !!----------------------------------------------------------------------
       !!----------------------------------------------------------------------
-      INTEGER,  INTENT(in) :: ki, kj
-      REAL(wp), INTENT(in) :: pWn, pLn, pAbits, pmass_scale, pMnew, pnMbits, pz1_e1e2
-      !!----------------------------------------------------------------------
-      !
-      IF( .NOT. ln_bergdia ) RETURN
+      INTEGER , INTENT(in) ::   ki, kj
+      REAL(wp), INTENT(in) ::   pWn, pLn, pAbits, pmass_scale, pMnew, pnMbits, pz1_e1e2
+      !!----------------------------------------------------------------------
+      !
+      IF( .NOT.ln_bergdia )   RETURN
       virtual_area(ki,kj) = virtual_area(ki,kj) + ( pWn * pLn + pAbits ) * pmass_scale      ! m^2
       berg_mass(ki,kj)    = berg_mass(ki,kj) + pMnew * pz1_e1e2                             ! kg/m2
@@ -453,5 +465,5 @@
       !!----------------------------------------------------------------------
       !
-      IF( .NOT. ln_bergdia ) RETURN
+      IF( .NOT.ln_bergdia )   RETURN
       nspeeding_tickets = nspeeding_tickets + 1
       !
@@ -459,7 +471,7 @@
 
 
-   SUBROUTINE icb_dia_melt(ki, kj, pmnew, pheat, pmass_scale,   &
-      &                   pdM, pdMbitsE, pdMbitsM, pdMb, pdMe,   &
-      &                   pdMv, pz1_dt_e1e2 )
+   SUBROUTINE icb_dia_melt(ki, kj, pmnew, pheat, pmass_scale,     &
+      &                    pdM, pdMbitsE, pdMbitsM, pdMb, pdMe,   &
+      &                    pdMv, pz1_dt_e1e2 )
       !!----------------------------------------------------------------------
       !!----------------------------------------------------------------------
@@ -469,6 +481,6 @@
       !!----------------------------------------------------------------------
       !
-      IF( .NOT. ln_bergdia ) RETURN
-
+      IF( .NOT.ln_bergdia )   RETURN
+      !
       berg_melt (ki,kj) = berg_melt (ki,kj) + pdM      * pz1_dt_e1e2   ! kg/m2/s
       bits_src  (ki,kj) = bits_src  (ki,kj) + pdMbitsE * pz1_dt_e1e2   ! mass flux into bergy bitskg/m2/s
@@ -492,17 +504,18 @@
       !!----------------------------------------------------------------------
       !
-      IF ( PRESENT(kbergs) ) THEN
+      IF( PRESENT(kbergs) ) THEN
          WRITE(numicb,100) cd_budgetstr // ' state:',                                    &
-                           cd_startstr  // ' start',  pstartval,         cd_budgetunits, &
-                           cd_endstr    // ' end',    pendval,           cd_budgetunits, &
-                           'Delta '     // cd_delstr, pendval-pstartval, cd_budgetunits, &
-                           '# of bergs', kbergs
+            &              cd_startstr  // ' start',  pstartval,         cd_budgetunits, &
+            &              cd_endstr    // ' end',    pendval,           cd_budgetunits, &
+            &              'Delta '     // cd_delstr, pendval-pstartval, cd_budgetunits, &
+            &              '# of bergs', kbergs
       ELSE
          WRITE(numicb,100) cd_budgetstr // ' state:',                                   &
-                           cd_startstr  // ' start', pstartval,         cd_budgetunits, &
-                           cd_endstr    // ' end',   pendval,           cd_budgetunits, &
-                           cd_delstr    // 'Delta',  pendval-pstartval, cd_budgetunits
+            &              cd_startstr  // ' start', pstartval,         cd_budgetunits, &
+            &              cd_endstr    // ' end',   pendval,           cd_budgetunits, &
+            &              cd_delstr    // 'Delta',  pendval-pstartval, cd_budgetunits
       ENDIF
-      100 FORMAT(a19,3(a18,"=",es14.7,x,a2,:,","),a12,i8)
+100   FORMAT(a19,3(a18,"=",es14.7,x,a2,:,","),a12,i8)
+      !
    END SUBROUTINE report_state
 
@@ -516,8 +529,9 @@
       !
       WRITE(numicb,200) cd_budgetstr // ' check:',                 &
-                        cd_startstr,    pstartval, cd_budgetunits, &
-                        cd_endstr,      pendval,   cd_budgetunits, &
-                        'error',        (pendval-pstartval)/((pendval+pstartval)+1e-30), 'nd'
-      200 FORMAT(a19,10(a18,"=",es14.7,x,a2,:,","))
+         &              cd_startstr,    pstartval, cd_budgetunits, &
+         &              cd_endstr,      pendval,   cd_budgetunits, &
+         &              'error',        (pendval-pstartval)/((pendval+pstartval)+1e-30), 'nd'
+200   FORMAT(a19,10(a18,"=",es14.7,x,a2,:,","))
+      !
    END SUBROUTINE report_consistant
 
@@ -530,10 +544,10 @@
       REAL(wp),      INTENT(in) :: pinval, poutval, pstartval, pendval
       !
-      REAL(wp)                  :: zval
+      REAL(wp) ::   zval
       !!----------------------------------------------------------------------
       !
       zval = ( ( pendval - pstartval ) - ( pinval - poutval ) ) /   &
-         &   MAX( 1.e-30, MAX( abs( pendval - pstartval ) , ABS( pinval - poutval ) ) )
-
+         &   MAX( 1.e-30, MAX( ABS( pendval - pstartval ) , ABS( pinval - poutval ) ) )
+         !
       WRITE(numicb,200) cd_budgetstr // ' budget:', &
          &              cd_instr     // ' in',      pinval,         cd_budgetunits, &
@@ -549,6 +563,7 @@
       !!----------------------------------------------------------------------
       !!----------------------------------------------------------------------
-      CHARACTER*(*), INTENT(in) :: cd_budgetstr, cd_startstr, cd_endstr, cd_delstr
-      INTEGER,       INTENT(in) :: pstartval, pendval
+      CHARACTER*(*), INTENT(in) ::   cd_budgetstr, cd_startstr, cd_endstr, cd_delstr
+      INTEGER      , INTENT(in) ::   pstartval, pendval
+      !!----------------------------------------------------------------------
       !
       WRITE(numicb,100) cd_budgetstr // ' state:',           &
@@ -570,9 +585,10 @@
       !
       WRITE(numicb,200) cd_budgetstr // ' budget:', &
-                        cd_instr     // ' in',      pinval, &
-                        cd_outstr    // ' out',     poutval, &
-                        'Delta '     // cd_delstr,  pinval-poutval, &
-                        'error',                    ( ( pendval - pstartval ) - ( pinval - poutval ) )
-      200 FORMAT(a19,10(a18,"=",i14,x,:,","))
+         &              cd_instr     // ' in',      pinval, &
+         &              cd_outstr    // ' out',     poutval, &
+         &              'Delta '     // cd_delstr,  pinval-poutval, &
+         &              'error',                    ( ( pendval - pstartval ) - ( pinval - poutval ) )
+200   FORMAT(a19,10(a18,"=",i14,x,:,","))
+      !
    END SUBROUTINE report_ibudget
 
Index: /branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/ICB/icbdyn.F90
===================================================================
--- /branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/ICB/icbdyn.F90	(revision 9189)
+++ /branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/ICB/icbdyn.F90	(revision 9190)
@@ -1,15 +1,13 @@
 MODULE icbdyn
-
    !!======================================================================
    !!                       ***  MODULE  icbdyn  ***
    !! Iceberg:  time stepping routine for iceberg tracking
    !!======================================================================
-   !! History : 3.3.1 !  2010-01  (Martin&Adcroft) Original code
-   !!            -    !  2011-03  (Madec)          Part conversion to NEMO form
-   !!            -    !                            Removal of mapping from another grid
-   !!            -    !  2011-04  (Alderson)       Split into separate modules
-   !!            -    !  2011-05  (Alderson)       Replace broken grounding routine
-   !!            -    !                            with one of Gurvan's suggestions (just like
-   !!            -    !                            the broken one)
+   !! History :  3.3  !  2010-01  (Martin&Adcroft)  Original code
+   !!             -   !  2011-03  (Madec)  Part conversion to NEMO form
+   !!             -   !                    Removal of mapping from another grid
+   !!             -   !  2011-04  (Alderson)  Split into separate modules
+   !!             -   !  2011-05  (Alderson)  Replace broken grounding routine with one of
+   !!             -   !                       Gurvan's suggestions (just like the broken one)
    !!----------------------------------------------------------------------
    USE par_oce        ! NEMO parameters
@@ -41,16 +39,17 @@
       !! ** Method  : - See Martin & Adcroft, Ocean Modelling 34, 2010
       !!----------------------------------------------------------------------
-      REAL(wp)                        ::   zuvel1 , zvvel1 , zu1, zv1, zax1, zay1, zxi1 , zyj1
-      REAL(wp)                        ::   zuvel2 , zvvel2 , zu2, zv2, zax2, zay2, zxi2 , zyj2
-      REAL(wp)                        ::   zuvel3 , zvvel3 , zu3, zv3, zax3, zay3, zxi3 , zyj3
-      REAL(wp)                        ::   zuvel4 , zvvel4 , zu4, zv4, zax4, zay4, zxi4 , zyj4
-      REAL(wp)                        ::   zuvel_n, zvvel_n, zxi_n   , zyj_n
-      REAL(wp)                        ::   zdt, zdt_2, zdt_6, ze1, ze2
-      LOGICAL                         ::   ll_bounced
-      TYPE(iceberg), POINTER          ::   berg
-      TYPE(point)  , POINTER          ::   pt
-      INTEGER                         ::   kt
-      !!----------------------------------------------------------------------
-
+      INTEGER, INTENT(in) ::   kt   !
+      !
+      LOGICAL  ::   ll_bounced
+      REAL(wp) ::   zuvel1 , zvvel1 , zu1, zv1, zax1, zay1, zxi1 , zyj1
+      REAL(wp) ::   zuvel2 , zvvel2 , zu2, zv2, zax2, zay2, zxi2 , zyj2
+      REAL(wp) ::   zuvel3 , zvvel3 , zu3, zv3, zax3, zay3, zxi3 , zyj3
+      REAL(wp) ::   zuvel4 , zvvel4 , zu4, zv4, zax4, zay4, zxi4 , zyj4
+      REAL(wp) ::   zuvel_n, zvvel_n, zxi_n   , zyj_n
+      REAL(wp) ::   zdt, zdt_2, zdt_6, ze1, ze2
+      TYPE(iceberg), POINTER ::   berg
+      TYPE(point)  , POINTER ::   pt
+      !!----------------------------------------------------------------------
+      !
       ! 4th order Runge-Kutta to solve:   d/dt X = V,  d/dt V = A
       !                    with I.C.'s:   X=X1 and V=V1
@@ -75,5 +74,5 @@
          pt => berg%current_point
 
-         ll_bounced = .false.
+         ll_bounced = .FALSE.
 
 
@@ -99,5 +98,5 @@
          !
          CALL icb_ground( zxi2, zxi1, zu1,   &
-         &                zyj2, zyj1, zv1, ll_bounced )
+            &             zyj2, zyj1, zv1, ll_bounced )
 
          !                                         !**   A2 = A(X2,V2)
@@ -115,5 +114,5 @@
          !
          CALL icb_ground( zxi3, zxi1, zu3,   &
-         &                zyj3, zyj1, zv3, ll_bounced )
+            &                zyj3, zyj1, zv3, ll_bounced )
 
          !                                         !**   A3 = A(X3,V3)
@@ -131,5 +130,5 @@
 
          CALL icb_ground( zxi4, zxi1, zu4,   &
-         &                zyj4, zyj1, zv4, ll_bounced )
+            &             zyj4, zyj1, zv4, ll_bounced )
 
          !                                         !**   A4 = A(X4,V4)
@@ -150,5 +149,5 @@
 
          CALL icb_ground( zxi_n, zxi1, zuvel_n,   &
-         &                      zyj_n, zyj1, zvvel_n, ll_bounced )
+            &             zyj_n, zyj1, zvvel_n, ll_bounced )
 
          pt%uvel = zuvel_n                        !** save in berg structure
@@ -169,5 +168,5 @@
 
    SUBROUTINE icb_ground( pi, pi0, pu,   &
-      &                         pj, pj0, pv, ld_bounced )
+      &                   pj, pj0, pv, ld_bounced )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE icb_ground  ***
@@ -216,18 +215,18 @@
       ibounce_method = 2
       SELECT CASE ( ibounce_method )
-         CASE ( 1 )
-            pi = pi0
-            pj = pj0
-            pu = 0._wp
-            pv = 0._wp
-         CASE ( 2 )
-            IF( ii0 /= ii ) THEN
-               pi = pi0                   ! return back to the initial position
-               pu = 0._wp                 ! zeroing of velocity in the direction of the grounding
-            ENDIF
-            IF( ij0 /= ij ) THEN
-               pj = pj0                   ! return back to the initial position
-               pv = 0._wp                 ! zeroing of velocity in the direction of the grounding
-            ENDIF
+      CASE ( 1 )
+         pi = pi0
+         pj = pj0
+         pu = 0._wp
+         pv = 0._wp
+      CASE ( 2 )
+         IF( ii0 /= ii ) THEN
+            pi = pi0                   ! return back to the initial position
+            pu = 0._wp                 ! zeroing of velocity in the direction of the grounding
+         ENDIF
+         IF( ij0 /= ij ) THEN
+            pj = pj0                   ! return back to the initial position
+            pv = 0._wp                 ! zeroing of velocity in the direction of the grounding
+         ENDIF
       END SELECT
       !
@@ -259,5 +258,6 @@
       !
       INTEGER  ::   itloop
-      REAL(wp) ::   zuo, zvo, zui, zvi, zua, zva, zuwave, zvwave, zssh_x, zssh_y, zsst, zcn, zhi
+      REAL(wp) ::   zuo, zui, zua, zuwave, zssh_x, zsst, zcn, zhi
+      REAL(wp) ::   zvo, zvi, zva, zvwave, zssh_y
       REAL(wp) ::   zff, zT, zD, zW, zL, zM, zF
       REAL(wp) ::   zdrag_ocn, zdrag_atm, zdrag_ice, zwave_rad
@@ -339,5 +339,5 @@
             zaxe = zaxe - zdrag_ocn*(puvel -zuo) - zdrag_atm*(puvel -zua) -zdrag_ice*(puvel -zui)
             zaye = zaye - zdrag_ocn*(pvvel -zvo) - zdrag_atm*(pvvel -zva) -zdrag_ice*(pvvel -zvi)
-         endif
+         ENDIF
 
          ! Solve for implicit accelerations
@@ -349,7 +349,7 @@
             pax     = zdetA * ( zA11*zaxe + zA12*zaye )
             pay     = zdetA * ( zA11*zaye - zA12*zaxe )
-         else
+         ELSE
             pax = zaxe   ;   pay = zaye
-         endif
+         ENDIF
 
          zuveln = puvel0 + pdt*pax
@@ -362,5 +362,5 @@
          IF( zspeed > 0._wp ) THEN
             zloc_dx = MIN( pe1, pe2 )                          ! minimum grid spacing
-            zspeed_new = zloc_dx / pdt * rn_speed_limit     ! Speed limit as a factor of dx / dt
+            zspeed_new = zloc_dx / pdt * rn_speed_limit        ! Speed limit as a factor of dx / dt
             IF( zspeed_new < zspeed ) THEN
                zuveln = zuveln * ( zspeed_new / zspeed )        ! Scale velocity to reduce speed
Index: /branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/ICB/icbini.F90
===================================================================
--- /branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/ICB/icbini.F90	(revision 9189)
+++ /branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/ICB/icbini.F90	(revision 9190)
@@ -361,5 +361,22 @@
       !!----------------------------------------------------------------------
 
-#if !defined key_agrif
+#if defined key_agrif
+      IF(lwp) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) 'icb_nam : AGRIF is not compatible with namelist namberg :  '
+         WRITE(numout,*) '~~~~~~~   definition of rn_initial_mass(nclasses) with nclasses as PARAMETER '
+         WRITE(numout,*)
+         WRITE(numout,*) '   ==>>>   force  NO icebergs used. The namelist namberg is not read'
+      ENDIF
+      ln_icebergs = .false.      
+      RETURN
+#else
+      IF(lwp) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) 'icb_nam : iceberg initialization through namberg namelist read'
+         WRITE(numout,*) '~~~~~~~~ '
+      ENDIF
+#endif   
+      !                             !==  read namelist  ==!
       REWIND( numnam_ref )              ! Namelist namberg in reference namelist : Iceberg parameters
       READ  ( numnam_ref, namberg, IOSTAT = ios, ERR = 901)
@@ -369,27 +386,19 @@
 902   IF( ios >  0 ) CALL ctl_nam ( ios , 'namberg in configuration namelist', lwp )
       IF(lwm) WRITE ( numond, namberg )
-#else
-      IF(lwp) THEN
-         WRITE(numout,*)
-         WRITE(numout,*) 'icbini : AGRIF is not compatible with namelist namberg :  '
-         WRITE(numout,*) '~~~~~~   definition of rn_initial_mass(nclasses) with nclasses as PARAMETER '
-         WRITE(numout,*) '         ==>>>   force  NO icebergs used. The namelist namberg is not read'
-      ENDIF
-      ln_icebergs = .false.      
-#endif   
-      IF( .NOT. ln_icebergs ) THEN   ! no icebergs
-         IF(lwp) THEN
-            WRITE(numout,*)
-            WRITE(numout,*) 'icbini : Namelist namberg ln_icebergs = F , NO icebergs used'
-            WRITE(numout,*) '~~~~~~ '
-         ENDIF
+      !
+      IF(lwp) WRITE(numout,*)
+      IF( ln_icebergs ) THEN
+         IF(lwp) WRITE(numout,*) '   ==>>>   icebergs are used'
+      ELSE
+         IF(lwp) WRITE(numout,*) '   ==>>>   No icebergs used'
          RETURN
       ENDIF
-
-
-!     IF( lk_lim3 .AND. ln_icebergs ) THEN
-!        CALL ctl_stop( 'icb_nam: the use of ICB with LIM3 not allowed. ice thickness missing in ICB' )
-!     ENDIF
-
+      !
+      IF( nn_test_icebergs > nclasses ) THEN
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) '   ==>>>   Resetting of nn_test_icebergs to ', nclasses
+         nn_test_icebergs = nclasses
+      ENDIF
+      !
       IF(lwp) THEN                  ! control print
          WRITE(numout,*)
@@ -399,18 +408,18 @@
          WRITE(numout,*) '   Period between sampling of position for trajectory storage   nn_sample_rate = ', nn_sample_rate
          WRITE(numout,*) '   Mass thresholds between iceberg classes (kg)                 rn_initial_mass     ='
-         DO jn=1,nclasses
-            WRITE(numout,'(a,f15.2)') '                                                                ',rn_initial_mass(jn)
+         DO jn = 1, nclasses
+            WRITE(numout,'(a,f15.2)') '                                                                ', rn_initial_mass(jn)
          ENDDO
          WRITE(numout,*) '   Fraction of calving to apply to this class (non-dim)         rn_distribution     ='
          DO jn = 1, nclasses
-            WRITE(numout,'(a,f10.4)') '                                                                ',rn_distribution(jn)
+            WRITE(numout,'(a,f10.4)') '                                                                ', rn_distribution(jn)
          END DO
          WRITE(numout,*) '   Ratio between effective and real iceberg mass (non-dim)      rn_mass_scaling     = '
          DO jn = 1, nclasses
-            WRITE(numout,'(a,f10.2)') '                                                                ',rn_mass_scaling(jn)
+            WRITE(numout,'(a,f10.2)') '                                                                ', rn_mass_scaling(jn)
          END DO
          WRITE(numout,*) '   Total thickness of newly calved bergs (m)                    rn_initial_thickness = '
          DO jn = 1, nclasses
-            WRITE(numout,'(a,f10.2)') '                                                                ',rn_initial_thickness(jn)
+            WRITE(numout,'(a,f10.2)') '                                                                ', rn_initial_thickness(jn)
          END DO
          WRITE(numout,*) '   Timesteps between verbose messages                           nn_verbose_write    = ', nn_verbose_write
@@ -435,9 +444,4 @@
       ENDIF
       !
-      IF( nn_test_icebergs > nclasses ) THEN
-         IF(lwp) WRITE(numout,*) '      ==>>>   Resetting of nn_test_icebergs to ', nclasses
-         nn_test_icebergs = nclasses
-      ENDIF
-
       ! ensure that the sum of berg input distribution is equal to one
       zfact = SUM( rn_distribution )
Index: /branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/ICB/icblbc.F90
===================================================================
--- /branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/ICB/icblbc.F90	(revision 9189)
+++ /branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/ICB/icblbc.F90	(revision 9190)
@@ -1,16 +1,15 @@
 MODULE icblbc
-
    !!======================================================================
    !!                       ***  MODULE  icblbc  ***
    !! Ocean physics:  routines to handle boundary exchanges for icebergs
    !!======================================================================
-   !! History : 3.3.1 !  2010-01  (Martin&Adcroft) Original code
-   !!            -    !  2011-03  (Madec)          Part conversion to NEMO form
-   !!            -    !                            Removal of mapping from another grid
-   !!            -    !  2011-04  (Alderson)       Split into separate modules
-   !!            -    !  2011-05  (Alderson)       MPP exchanges written based on lib_mpp
-   !!            -    !  2011-05  (Alderson)       MPP and single processor boundary
-   !!            -    !                            conditions added
+   !! History :  3.3  !  2010-01  (Martin&Adcroft) Original code
+   !!             -   !  2011-03  (Madec)          Part conversion to NEMO form
+   !!             -   !                            Removal of mapping from another grid
+   !!             -   !  2011-04  (Alderson)       Split into separate modules
+   !!             -   !  2011-05  (Alderson)       MPP exchanges written based on lib_mpp
+   !!             -   !  2011-05  (Alderson)       MPP and single processor boundary conditions added
    !!----------------------------------------------------------------------
+
    !!----------------------------------------------------------------------
    !!   icb_lbc       : -  Pass icebergs across cyclic boundaries
@@ -27,5 +26,4 @@
    !!                         nicbfldpts  - packed i,j point in exchanging processor
    !!----------------------------------------------------------------------
-
    USE par_oce                             ! ocean parameters
    USE dom_oce                             ! ocean domain
@@ -45,18 +43,18 @@
 
    TYPE, PUBLIC :: buffer
-      INTEGER :: size=0
-      REAL(wp), DIMENSION(:,:), POINTER :: data
+      INTEGER :: size = 0
+      REAL(wp), DIMENSION(:,:), POINTER ::   data
    END TYPE buffer
 
-   TYPE(buffer), POINTER           ::   obuffer_n=>NULL() , ibuffer_n=>NULL()
-   TYPE(buffer), POINTER           ::   obuffer_s=>NULL() , ibuffer_s=>NULL()
-   TYPE(buffer), POINTER           ::   obuffer_e=>NULL() , ibuffer_e=>NULL()
-   TYPE(buffer), POINTER           ::   obuffer_w=>NULL() , ibuffer_w=>NULL()
+   TYPE(buffer), POINTER       ::   obuffer_n=>NULL() , ibuffer_n=>NULL()
+   TYPE(buffer), POINTER       ::   obuffer_s=>NULL() , ibuffer_s=>NULL()
+   TYPE(buffer), POINTER       ::   obuffer_e=>NULL() , ibuffer_e=>NULL()
+   TYPE(buffer), POINTER       ::   obuffer_w=>NULL() , ibuffer_w=>NULL()
 
    ! north fold exchange buffers
-   TYPE(buffer), POINTER           ::   obuffer_f=>NULL() , ibuffer_f=>NULL()
-
-   INTEGER, PARAMETER, PRIVATE     ::   jp_delta_buf = 25             ! Size by which to increment buffers
-   INTEGER, PARAMETER, PRIVATE     ::   jp_buffer_width = 15+nkounts  ! items to store for each berg
+   TYPE(buffer), POINTER       ::   obuffer_f=>NULL() , ibuffer_f=>NULL()
+
+   INTEGER, PARAMETER, PRIVATE ::   jp_delta_buf = 25             ! Size by which to increment buffers
+   INTEGER, PARAMETER, PRIVATE ::   jp_buffer_width = 15+nkounts  ! items to store for each berg
 
 #endif
@@ -926,5 +924,4 @@
       WRITE(numout,*) 'icb_lbc_mpp: You should not have seen this message!!'
    END SUBROUTINE icb_lbc_mpp
-
 #endif
 
Index: /branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/ICB/icbrst.F90
===================================================================
--- /branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/ICB/icbrst.F90	(revision 9189)
+++ /branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/ICB/icbrst.F90	(revision 9190)
@@ -1,4 +1,3 @@
 MODULE icbrst
-
    !!======================================================================
    !!                       ***  MODULE  icbrst  ***
@@ -15,4 +14,5 @@
    !!                                              read single restart files
    !!----------------------------------------------------------------------
+
    !!----------------------------------------------------------------------
    !!   icb_rst_read    : read restart file
@@ -110,12 +110,11 @@
                CALL iom_get( ncid, 'mass_of_bits' , localpt%mass_of_bits  , ktime=jn )
                CALL iom_get( ncid, 'heat_density' , localpt%heat_density  , ktime=jn )
-
                !
                CALL icb_utl_add( localberg, localpt )
-
+               !
             ENDIF
-
+            !
          END DO
-
+         !
       ENDIF 
 
@@ -144,5 +143,5 @@
       CALL iom_close( ncid )
       !
-      IF( lwp .and. nn_verbose_level >= 0)  WRITE(numout,'(a)') 'icebergs, read_restart_bergs: completed'
+      IF( lwp .AND. nn_verbose_level >= 0)  WRITE(numout,'(a)') 'icebergs, read_restart_bergs: completed'
       !
    END SUBROUTINE icb_rst_read
@@ -361,3 +360,4 @@
    END SUBROUTINE icb_rst_write
    !
+   !!======================================================================
 END MODULE icbrst
Index: /branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/ICB/icbstp.F90
===================================================================
--- /branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/ICB/icbstp.F90	(revision 9189)
+++ /branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/ICB/icbstp.F90	(revision 9190)
@@ -171,4 +171,4 @@
    END SUBROUTINE icb_end
 
-   !!-------------------------------------------------------------------------
+   !!======================================================================
 END MODULE icbstp
Index: /branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/ICB/icbthm.F90
===================================================================
--- /branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/ICB/icbthm.F90	(revision 9189)
+++ /branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/ICB/icbthm.F90	(revision 9190)
@@ -1,4 +1,3 @@
 MODULE icbthm
-
    !!======================================================================
    !!                       ***  MODULE  icbthm  ***
@@ -31,5 +30,9 @@
    PUBLIC   icb_thm ! routine called in icbstp.F90 module
 
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 3.3 , NEMO Consortium (2011)
    !! $Id$
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
 CONTAINS
 
@@ -55,5 +58,5 @@
       !
       z1_rday = 1._wp / rday
-      
+      !
       ! we're either going to ignore berg fresh water melt flux and associated heat
       ! or we pass it into the ocean, so at this point we set them both to zero,
@@ -63,13 +66,13 @@
       berg_grid%floating_melt(:,:) = 0._wp
       berg_grid%calving_hflx(:,:)  = 0._wp
-    
+      !
       this => first_berg
-      DO WHILE( associated(this) )
+      DO WHILE( ASSOCIATED(this) )
          !
          pt => this%current_point
          nknberg = this%number(1)
-         CALL icb_utl_interp( pt%xi, pt%e1, pt%uo, pt%ui, pt%ua, pt%ssh_x, &
-         &                    pt%yj, pt%e2, pt%vo, pt%vi, pt%va, pt%ssh_y, &
-         &                 pt%sst, pt%cn, pt%hi, zff )
+         CALL icb_utl_interp( pt%xi, pt%e1, pt%uo, pt%ui, pt%ua, pt%ssh_x,   &
+            &                 pt%yj, pt%e2, pt%vo, pt%vi, pt%va, pt%ssh_y,   &
+            &                 pt%sst, pt%cn, pt%hi, zff )
          !
          zSST = pt%sst
@@ -98,51 +101,51 @@
          zMv = MAX( 7.62e-3*zSST+1.29e-3*(zSST**2)            , 0._wp ) * z1_rday   ! Buoyant convection at sides (eqn M.A10)
          zMb = MAX( 0.58*(zdvo**0.8)*(zSST+4.0)/(zL**0.2)      , 0._wp ) * z1_rday   ! Basal turbulent melting     (eqn M.A7 )
-         zMe = MAX( 1./12.*(zSST+2.)*zSs*(1+cos(rpi*(zIC**3))) , 0._wp ) * z1_rday   ! Wave erosion                (eqn M.A8 )
+         zMe = MAX( 1./12.*(zSST+2.)*zSs*(1+COS(rpi*(zIC**3))) , 0._wp ) * z1_rday   ! Wave erosion                (eqn M.A8 )
 
          IF( ln_operator_splitting ) THEN      ! Operator split update of volume/mass
             zTn    = MAX( zT - zMb*zdt , 0._wp )         ! new total thickness (m)
-            znVol  = zTn * zW * zL                        ! new volume (m^3)
-            zMnew1 = (znVol/zVol) * zM                    ! new mass (kg)
+            znVol  = zTn * zW * zL                       ! new volume (m^3)
+            zMnew1 = (znVol/zVol) * zM                   ! new mass (kg)
             zdMb   = zM - zMnew1                         ! mass lost to basal melting (>0) (kg)
             !
             zLn    = MAX( zL - zMv*zdt , 0._wp )         ! new length (m)
             zWn    = MAX( zW - zMv*zdt , 0._wp )         ! new width (m)
-            znVol  = zTn * zWn * zLn                      ! new volume (m^3)
-            zMnew2 = (znVol/zVol) * zM                    ! new mass (kg)
+            znVol  = zTn * zWn * zLn                     ! new volume (m^3)
+            zMnew2 = (znVol/zVol) * zM                   ! new mass (kg)
             zdMv   = zMnew1 - zMnew2                     ! mass lost to buoyant convection (>0) (kg)
             !
             zLn    = MAX( zLn - zMe*zdt , 0._wp )        ! new length (m)
             zWn    = MAX( zWn - zMe*zdt , 0._wp )        ! new width (m)
-            znVol  = zTn * zWn * zLn                      ! new volume (m^3)
-            zMnew  = ( znVol / zVol ) * zM                ! new mass (kg)
+            znVol  = zTn * zWn * zLn                     ! new volume (m^3)
+            zMnew  = ( znVol / zVol ) * zM               ! new mass (kg)
             zdMe   = zMnew2 - zMnew                      ! mass lost to erosion (>0) (kg)
             zdM    = zM - zMnew                          ! mass lost to all erosion and melting (>0) (kg)
             !
          ELSE                                         ! Update dimensions of berg
-            zLn = MAX( zL -(zMv+zMe)*zdt ,0._wp )         ! (m)
-            zWn = MAX( zW -(zMv+zMe)*zdt ,0._wp )         ! (m)
+            zLn = MAX( zL -(zMv+zMe)*zdt ,0._wp )        ! (m)
+            zWn = MAX( zW -(zMv+zMe)*zdt ,0._wp )        ! (m)
             zTn = MAX( zT - zMb    *zdt ,0._wp )         ! (m)
             ! Update volume and mass of berg
-            znVol = zTn*zWn*zLn                           ! (m^3)
-            zMnew = (znVol/zVol)*zM                       ! (kg)
+            znVol = zTn*zWn*zLn                          ! (m^3)
+            zMnew = (znVol/zVol)*zM                      ! (kg)
             zdM   = zM - zMnew                           ! (kg)
-            zdMb = (zM/zVol) * (zW*   zL ) *zMb*zdt         ! approx. mass loss to basal melting (kg)
-            zdMe = (zM/zVol) * (zT*(zW+zL)) *zMe*zdt         ! approx. mass lost to erosion (kg)
-            zdMv = (zM/zVol) * (zT*(zW+zL)) *zMv*zdt         ! approx. mass loss to buoyant convection (kg)
-         ENDIF
-
-         IF( rn_bits_erosion_fraction > 0._wp ) THEN      ! Bergy bits
+            zdMb = (zM/zVol) * (zW*   zL ) *zMb*zdt      ! approx. mass loss to basal melting (kg)
+            zdMe = (zM/zVol) * (zT*(zW+zL)) *zMe*zdt     ! approx. mass lost to erosion (kg)
+            zdMv = (zM/zVol) * (zT*(zW+zL)) *zMv*zdt     ! approx. mass loss to buoyant convection (kg)
+         ENDIF
+
+         IF( rn_bits_erosion_fraction > 0._wp ) THEN     ! Bergy bits
             !
             zMbits   = pt%mass_of_bits                                               ! mass of bergy bits (kg)
-            zdMbitsE = rn_bits_erosion_fraction * zdMe                        ! change in mass of bits (kg)
-            znMbits  = zMbits + zdMbitsE                                               ! add new bergy bits to mass (kg)
-            zLbits   = MIN( zL, zW, zT, 40._wp )                                        ! assume bergy bits are smallest dimension or 40 meters
-            zAbits   = ( zMbits / rn_rho_bergs ) / zLbits                           ! Effective bottom area (assuming T=Lbits)
-            zMbb     = MAX( 0.58*(zdvo**0.8)*(zSST+2.0)/(zLbits**0.2), 0.) * z1_rday    ! Basal turbulent melting (for bits)
-            zMbb     = rn_rho_bergs * zAbits * zMbb                                 ! in kg/s
-            zdMbitsM = MIN( zMbb*zdt , znMbits )                                       ! bergy bits mass lost to melting (kg)
-            znMbits  = znMbits-zdMbitsM                                                ! remove mass lost to bergy bits melt
+            zdMbitsE = rn_bits_erosion_fraction * zdMe                               ! change in mass of bits (kg)
+            znMbits  = zMbits + zdMbitsE                                             ! add new bergy bits to mass (kg)
+            zLbits   = MIN( zL, zW, zT, 40._wp )                                     ! assume bergy bits are smallest dimension or 40 meters
+            zAbits   = ( zMbits / rn_rho_bergs ) / zLbits                            ! Effective bottom area (assuming T=Lbits)
+            zMbb     = MAX( 0.58*(zdvo**0.8)*(zSST+2.0)/(zLbits**0.2), 0.) * z1_rday ! Basal turbulent melting (for bits)
+            zMbb     = rn_rho_bergs * zAbits * zMbb                                  ! in kg/s
+            zdMbitsM = MIN( zMbb*zdt , znMbits )                                     ! bergy bits mass lost to melting (kg)
+            znMbits  = znMbits-zdMbitsM                                              ! remove mass lost to bergy bits melt
             IF( zMnew == 0._wp ) THEN                                                ! if parent berg has completely melted then
-               zdMbitsM = zdMbitsM + znMbits                                           ! instantly melt all the bergy bits
+               zdMbitsM = zdMbitsM + znMbits                                         ! instantly melt all the bergy bits
                znMbits  = 0._wp
             ENDIF
@@ -163,6 +166,6 @@
             berg_grid%calving_hflx (ii,ij) = berg_grid%calving_hflx (ii,ij) + zheat    * z1_e1e2    ! W/m2
             CALL icb_dia_melt( ii, ij, zMnew, zheat, this%mass_scaling,       &
-            &                          zdM, zdMbitsE, zdMbitsM, zdMb, zdMe,   &
-            &                          zdMv, z1_dt_e1e2 )
+               &                       zdM, zdMbitsE, zdMbitsM, zdMb, zdMe,   &
+               &                       zdMv, z1_dt_e1e2 )
          ELSE
             WRITE(numout,*) 'icb_thm: berg ',this%number(:),' appears to have grounded  at ',narea,ii,ij
@@ -178,5 +181,5 @@
             zTn = zWn
             zWn = zT
-         endif
+         ENDIF
 
          ! Store the new state of iceberg (with L>W)
@@ -184,6 +187,6 @@
          pt%mass_of_bits = znMbits
          pt%thickness    = zTn
-         pt%width        = min(zWn,zLn)
-         pt%length       = max(zWn,zLn)
+         pt%width        = MIN( zWn , zLn )
+         pt%length       = MAX( zWn , zLn )
 
          next=>this%next
@@ -197,5 +200,5 @@
             z1_e1e2 = r1_e1e2t(ii,ij) * this%mass_scaling
             CALL icb_dia_size( ii, ij, zWn, zLn, zAbits,   &
-            &                  this%mass_scaling, zMnew, znMbits, z1_e1e2)
+               &               this%mass_scaling, zMnew, znMbits, z1_e1e2 )
          ENDIF
          !
@@ -203,5 +206,5 @@
          !
       END DO
-      
+
       ! now use melt and associated heat flux in ocean (or not)
       !
Index: /branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/ICB/icbtrj.F90
===================================================================
--- /branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/ICB/icbtrj.F90	(revision 9189)
+++ /branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/ICB/icbtrj.F90	(revision 9190)
@@ -1,15 +1,18 @@
 MODULE icbtrj
-
    !!======================================================================
    !!                       ***  MODULE  icbtrj  ***
    !! Ocean physics:  trajectory I/O routines
    !!======================================================================
-   !! History : 3.3.1 !  2010-01  (Martin&Adcroft) Original code
-   !!            -    !  2011-03  (Madec)          Part conversion to NEMO form
-   !!            -    !                            Removal of mapping from another grid
-   !!            -    !  2011-05  (Alderson)       New module to handle trajectory output
-   !!----------------------------------------------------------------------
-   !!----------------------------------------------------------------------
-   !!   icb_trj_init          :
+   !! History :  3.3  !  2010-01  (Martin&Adcroft) Original code
+   !!             -   !  2011-03  (Madec)          Part conversion to NEMO form
+   !!             -   !                            Removal of mapping from another grid
+   !!             -   !  2011-05  (Alderson)       New module to handle trajectory output
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   icb_trj_init  :
+   !!   icb_trj_write :
+   !!   icb_trj_sync  :
+   !!   icb_trj_end   :
    !!----------------------------------------------------------------------
    USE par_oce        ! NEMO parameters
@@ -49,6 +52,4 @@
    !!----------------------------------------------------------------------
 CONTAINS
-
-   !!-------------------------------------------------------------------------
 
    SUBROUTINE icb_trj_init( ktend )
@@ -252,5 +253,4 @@
    END SUBROUTINE icb_trj_write
 
-   !!-------------------------------------------------------------------------
 
    SUBROUTINE icb_trj_sync()
@@ -260,5 +260,5 @@
       !! ** Purpose :   
       !!----------------------------------------------------------------------
-      INTEGER                               :: iret
+      INTEGER ::   iret
       !!----------------------------------------------------------------------
       ! flush to file
@@ -270,6 +270,6 @@
 
    SUBROUTINE icb_trj_end()
-      ! Local variables
-      INTEGER                               :: iret
+      !!----------------------------------------------------------------------
+      INTEGER ::   iret
       !!----------------------------------------------------------------------
       ! Finish up
@@ -279,5 +279,4 @@
    END SUBROUTINE icb_trj_end
 
-   !!-------------------------------------------------------------------------
-
+   !!======================================================================
 END MODULE icbtrj
Index: /branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/ICB/icbutl.F90
===================================================================
--- /branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/ICB/icbutl.F90	(revision 9189)
+++ /branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/ICB/icbutl.F90	(revision 9190)
@@ -9,4 +9,5 @@
    !!            -    !  2011-04  (Alderson)       Split into separate modules
    !!----------------------------------------------------------------------
+
    !!----------------------------------------------------------------------
    !!   icb_utl_interp   :
@@ -48,6 +49,5 @@
    !! $Id$
    !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
-   !!-------------------------------------------------------------------------
-
+   !!----------------------------------------------------------------------
 CONTAINS
 
@@ -65,12 +65,12 @@
       ! and ssh which is used to calculate gradients
 
-      uo_e(:,:) = 0._wp ;   uo_e(1:jpi, 1:jpj) = ssu_m(:,:) * umask(:,:,1)
-      vo_e(:,:) = 0._wp ;   vo_e(1:jpi, 1:jpj) = ssv_m(:,:) * vmask(:,:,1)
-      ff_e(:,:) = 0._wp ;   ff_e(1:jpi, 1:jpj) = ff_f (:,:) 
-      tt_e(:,:) = 0._wp ;   tt_e(1:jpi, 1:jpj) = sst_m(:,:)
-      fr_e(:,:) = 0._wp ;   fr_e(1:jpi, 1:jpj) = fr_i (:,:)
-      ua_e(:,:) = 0._wp ;   ua_e(1:jpi, 1:jpj) = utau (:,:) * umask(:,:,1) ! maybe mask useless because mask applied in sbcblk
-      va_e(:,:) = 0._wp ;   va_e(1:jpi, 1:jpj) = vtau (:,:) * vmask(:,:,1) ! maybe mask useless because mask applied in sbcblk
-
+      uo_e(:,:) = 0._wp   ;   uo_e(1:jpi,1:jpj) = ssu_m(:,:) * umask(:,:,1)
+      vo_e(:,:) = 0._wp   ;   vo_e(1:jpi,1:jpj) = ssv_m(:,:) * vmask(:,:,1)
+      ff_e(:,:) = 0._wp   ;   ff_e(1:jpi,1:jpj) = ff_f (:,:) 
+      tt_e(:,:) = 0._wp   ;   tt_e(1:jpi,1:jpj) = sst_m(:,:)
+      fr_e(:,:) = 0._wp   ;   fr_e(1:jpi,1:jpj) = fr_i (:,:)
+      ua_e(:,:) = 0._wp   ;   ua_e(1:jpi,1:jpj) = utau (:,:) * umask(:,:,1) ! maybe mask useless because mask applied in sbcblk
+      va_e(:,:) = 0._wp   ;   va_e(1:jpi,1:jpj) = vtau (:,:) * vmask(:,:,1) ! maybe mask useless because mask applied in sbcblk
+      !
       CALL lbc_lnk_icb( uo_e, 'U', -1._wp, 1, 1 )
       CALL lbc_lnk_icb( vo_e, 'V', -1._wp, 1, 1 )
@@ -84,7 +84,8 @@
       ui_e(:,:) = 0._wp ;   ui_e(1:jpi, 1:jpj) = u_ice(:,:)
       vi_e(:,:) = 0._wp ;   vi_e(1:jpi, 1:jpj) = v_ice(:,:)
-      CALL lbc_lnk_icb(hicth, 'T', +1._wp, 1, 1 )
-      CALL lbc_lnk_icb( ui_e, 'U', -1._wp, 1, 1 )
-      CALL lbc_lnk_icb( vi_e, 'V', -1._wp, 1, 1 )
+      !
+      CALL lbc_lnk_icb( hicth, 'T', +1._wp, 1, 1 )
+      CALL lbc_lnk_icb( ui_e , 'U', -1._wp, 1, 1 )
+      CALL lbc_lnk_icb( vi_e , 'V', -1._wp, 1, 1 )
 #endif
 
@@ -149,7 +150,7 @@
 
 #if defined key_lim3
-      pui = icb_utl_bilin_h( ui_e, pi, pj, 'U' )              ! sea-ice velocities
-      pvi = icb_utl_bilin_h( vi_e, pi, pj, 'V' )
-      phi = icb_utl_bilin_h(hicth, pi, pj, 'T' )              ! ice thickness
+      pui = icb_utl_bilin_h( ui_e , pi, pj, 'U' )              ! sea-ice velocities
+      pvi = icb_utl_bilin_h( vi_e , pi, pj, 'V' )
+      phi = icb_utl_bilin_h( hicth, pi, pj, 'T' )              ! ice thickness
 #else
       pui = 0._wp
@@ -160,7 +161,7 @@
       ! Estimate SSH gradient in i- and j-direction (centred evaluation)
       pssh_i = ( icb_utl_bilin_h( ssh_e, pi+0.1_wp, pj, 'T' ) -   &
-          &      icb_utl_bilin_h( ssh_e, pi-0.1_wp, pj, 'T' )  ) / ( 0.2_wp * pe1 )
+         &       icb_utl_bilin_h( ssh_e, pi-0.1_wp, pj, 'T' )  ) / ( 0.2_wp * pe1 )
       pssh_j = ( icb_utl_bilin_h( ssh_e, pi, pj+0.1_wp, 'T' ) -   &
-          &      icb_utl_bilin_h( ssh_e, pi, pj-0.1_wp, 'T' )  ) / ( 0.2_wp * pe2 )
+         &       icb_utl_bilin_h( ssh_e, pi, pj-0.1_wp, 'T' )  ) / ( 0.2_wp * pe2 )
       !
    END SUBROUTINE icb_utl_interp
@@ -181,4 +182,68 @@
       REAL(wp)                            , INTENT(in) ::   pi, pj    ! targeted coordinates in (i,j) referential
       CHARACTER(len=1)                    , INTENT(in) ::   cd_type   ! type of pfld array grid-points: = T , U , V or F points
+      !
+      INTEGER  ::   ii, ij   ! local integer
+      REAL(wp) ::   zi, zj   ! local real
+      !!----------------------------------------------------------------------
+      !
+      SELECT CASE ( cd_type )
+      CASE ( 'T' )
+         ! note that here there is no +0.5 added
+         ! since we're looking for four T points containing quadrant we're in of 
+         ! current T cell
+         ii = MAX(1, INT( pi     ))
+         ij = MAX(1, INT( pj     ))    ! T-point
+         zi = pi - REAL(ii,wp)
+         zj = pj - REAL(ij,wp)
+      CASE ( 'U' )
+         ii = MAX(1, INT( pi-0.5 ))
+         ij = MAX(1, INT( pj     ))    ! U-point
+         zi = pi - 0.5 - REAL(ii,wp)
+         zj = pj - REAL(ij,wp)
+      CASE ( 'V' )
+         ii = MAX(1, INT( pi     ))
+         ij = MAX(1, INT( pj-0.5 ))    ! V-point
+         zi = pi - REAL(ii,wp)
+         zj = pj - 0.5 - REAL(ij,wp)
+      CASE ( 'F' )
+         ii = MAX(1, INT( pi-0.5 ))
+         ij = MAX(1, INT( pj-0.5 ))    ! F-point
+         zi = pi - 0.5 - REAL(ii,wp)
+         zj = pj - 0.5 - REAL(ij,wp)
+      END SELECT
+      !
+      ! find position in this processor. Prevent near edge problems (see #1389)
+      !
+      IF    ( ii < mig( 1 ) ) THEN   ;   ii = 1
+      ELSEIF( ii > mig(jpi) ) THEN   ;   ii = jpi
+      ELSE                           ;   ii = mi1(ii)
+      ENDIF
+      IF    ( ij < mjg( 1 ) ) THEN   ;   ij = 1
+      ELSEIF( ij > mjg(jpj) ) THEN   ;   ij = jpj
+      ELSE                           ;   ij  = mj1(ij)
+      ENDIF
+      !
+      IF( ii == jpi )   ii = ii-1      
+      IF( ij == jpj )   ij = ij-1
+      !
+      icb_utl_bilin_h = ( pfld(ii,ij  ) * (1.-zi) + pfld(ii+1,ij  ) * zi ) * (1.-zj)   &
+         &            + ( pfld(ii,ij+1) * (1.-zi) + pfld(ii+1,ij+1) * zi ) *     zj
+      !
+   END FUNCTION icb_utl_bilin_h
+
+
+   REAL(wp) FUNCTION icb_utl_bilin( pfld, pi, pj, cd_type )
+      !!----------------------------------------------------------------------
+      !!                  ***  FUNCTION icb_utl_bilin  ***
+      !!
+      !! ** Purpose :   bilinear interpolation at berg location depending on the grid-point type
+      !!
+      !!       !!gm  CAUTION an optional argument should be added to handle
+      !!             the slip/no-slip conditions  ==>>> to be done later
+      !!
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pfld      ! field to be interpolated
+      REAL(wp)                    , INTENT(in) ::   pi, pj    ! targeted coordinates in (i,j) referential
+      CHARACTER(len=1)            , INTENT(in) ::   cd_type   ! type of pfld array grid-points: = T , U , V or F points
       !
       INTEGER  ::   ii, ij   ! local integer
@@ -213,96 +278,16 @@
       !
       ! find position in this processor. Prevent near edge problems (see #1389)
-
-      if (ii.lt.mig(1)) then
-        ii = 1
-      else if (ii.gt.mig(jpi)) then
-        ii = jpi
-      else
-        ii  = mi1( ii  )
-      end if
-
-      if (ij.lt.mjg(1)) then
-        ij = 1
-      else if (ij.gt.mjg(jpj)) then
-        ij = jpj
-      else
-        ij  = mj1( ij  )
-      end if
-
-      if (ij.eq.jpj) ij=ij-1
-      if (ii.eq.jpi) ii=ii-1      
-
-      !
-      icb_utl_bilin_h = ( pfld(ii,ij  ) * (1.-zi) + pfld(ii+1,ij  ) * zi ) * (1.-zj)   &
-         &            + ( pfld(ii,ij+1) * (1.-zi) + pfld(ii+1,ij+1) * zi ) *     zj
-      !
-   END FUNCTION icb_utl_bilin_h
-
-
-   REAL(wp) FUNCTION icb_utl_bilin( pfld, pi, pj, cd_type )
-      !!----------------------------------------------------------------------
-      !!                  ***  FUNCTION icb_utl_bilin  ***
-      !!
-      !! ** Purpose :   bilinear interpolation at berg location depending on the grid-point type
-      !!
-      !!       !!gm  CAUTION an optional argument should be added to handle
-      !!             the slip/no-slip conditions  ==>>> to be done later
-      !!
-      !!----------------------------------------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pfld      ! field to be interpolated
-      REAL(wp)                    , INTENT(in) ::   pi, pj    ! targeted coordinates in (i,j) referential
-      CHARACTER(len=1)            , INTENT(in) ::   cd_type   ! type of pfld array grid-points: = T , U , V or F points
-      !
-      INTEGER  ::   ii, ij   ! local integer
-      REAL(wp) ::   zi, zj   ! local real
-      !!----------------------------------------------------------------------
-      !
-      SELECT CASE ( cd_type )
-         CASE ( 'T' )
-            ! note that here there is no +0.5 added
-            ! since we're looking for four T points containing quadrant we're in of 
-            ! current T cell
-            ii = MAX(1, INT( pi     ))
-            ij = MAX(1, INT( pj     ))    ! T-point
-            zi = pi - REAL(ii,wp)
-            zj = pj - REAL(ij,wp)
-         CASE ( 'U' )
-            ii = MAX(1, INT( pi-0.5 ))
-            ij = MAX(1, INT( pj     ))    ! U-point
-            zi = pi - 0.5 - REAL(ii,wp)
-            zj = pj - REAL(ij,wp)
-         CASE ( 'V' )
-            ii = MAX(1, INT( pi     ))
-            ij = MAX(1, INT( pj-0.5 ))    ! V-point
-            zi = pi - REAL(ii,wp)
-            zj = pj - 0.5 - REAL(ij,wp)
-         CASE ( 'F' )
-            ii = MAX(1, INT( pi-0.5 ))
-            ij = MAX(1, INT( pj-0.5 ))    ! F-point
-            zi = pi - 0.5 - REAL(ii,wp)
-            zj = pj - 0.5 - REAL(ij,wp)
-      END SELECT
-      !
-      ! find position in this processor. Prevent near edge problems (see #1389)
-
-      if (ii.lt.mig(1)) then
-        ii = 1
-      else if (ii.gt.mig(jpi)) then
-        ii = jpi
-      else
-        ii  = mi1( ii  )
-      end if
-
-      if (ij.lt.mjg(1)) then
-        ij = 1
-      else if (ij.gt.mjg(jpj)) then
-        ij = jpj
-      else
-        ij  = mj1( ij  )
-      end if
-
-      if (ij.eq.jpj) ij=ij-1
-      if (ii.eq.jpi) ii=ii-1
-
+      IF    ( ii < mig( 1 ) ) THEN   ;   ii = 1
+      ELSEIF( ii > mig(jpi) ) THEN   ;   ii = jpi
+      ELSE                           ;   ii = mi1(ii)
+      ENDIF
+      IF    ( ij < mjg( 1 ) ) THEN   ;   ij = 1
+      ELSEIF( ij > mjg(jpj) ) THEN   ;   ij = jpj
+      ELSE                           ;   ij  = mj1(ij)
+      ENDIF
+      !
+      IF( ii == jpi )   ii = ii-1      
+      IF( ij == jpj )   ij = ij-1
+      !
       icb_utl_bilin = ( pfld(ii,ij  ) * (1.-zi) + pfld(ii+1,ij  ) * zi ) * (1.-zj)   &
          &          + ( pfld(ii,ij+1) * (1.-zi) + pfld(ii+1,ij+1) * zi ) *     zj
@@ -340,24 +325,16 @@
       !
       ! find position in this processor. Prevent near edge problems (see #1389)
-
-      if (ii.lt.mig(1)) then
-        ii = 1
-      else if (ii.gt.mig(jpi)) then
-        ii = jpi
-      else
-        ii  = mi1( ii  )
-      end if
-
-      if (ij.lt.mjg(1)) then
-        ij = 1
-      else if (ij.gt.mjg(jpj)) then
-        ij = jpj
-      else
-        ij  = mj1( ij  )
-      end if
-
-      if (ij.eq.jpj) ij=ij-1
-      if (ii.eq.jpi) ii=ii-1
-
+      IF    ( ii < mig( 1 ) ) THEN   ;   ii = 1
+      ELSEIF( ii > mig(jpi) ) THEN   ;   ii = jpi
+      ELSE                           ;   ii = mi1(ii)
+      ENDIF
+      IF    ( ij < mjg( 1 ) ) THEN   ;   ij = 1
+      ELSEIF( ij > mjg(jpj) ) THEN   ;   ij = jpj
+      ELSE                           ;   ij  = mj1(ij)
+      ENDIF
+      !
+      IF( ii == jpi )   ii = ii-1      
+      IF( ij == jpj )   ij = ij-1
+      !
       z4(1) = pfld(ii  ,ij  )
       z4(2) = pfld(ii+1,ij  )
@@ -408,24 +385,16 @@
 
       ! find position in this processor. Prevent near edge problems (see #1389)
-
-      if (ii.lt.mig(1)) then
-        ii = 1
-      else if (ii.gt.mig(jpi)) then
-        ii = jpi
-      else
-        ii  = mi1( ii  )
-      end if
-
-      if (ij.lt.mjg(1)) then
-        ij = 1
-      else if (ij.gt.mjg(jpj)) then
-        ij = jpj
-      else
-        ij  = mj1( ij  )
-      end if
-
-      if (ij.eq.jpj) ij=ij-1
-      if (ii.eq.jpi) ii=ii-1
-
+      IF    ( ii < mig( 1 ) ) THEN   ;   ii = 1
+      ELSEIF( ii > mig(jpi) ) THEN   ;   ii = jpi
+      ELSE                           ;   ii = mi1(ii)
+      ENDIF
+      IF    ( ij < mjg( 1 ) ) THEN   ;   ij = 1
+      ELSEIF( ij > mjg(jpj) ) THEN   ;   ij = jpj
+      ELSE                           ;   ij  = mj1(ij)
+      ENDIF
+      !
+      IF( ii == jpi )   ii = ii-1      
+      IF( ij == jpj )   ij = ij-1
+      !
       IF(    0.0_wp <= zi .AND. zi < 0.5_wp   ) THEN
          IF( 0.0_wp <= zj .AND. zj < 0.5_wp        )   THEN        !  NE quadrant
Index: /branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/LBC/lbc_nfd_ext_generic.h90
===================================================================
--- /branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/LBC/lbc_nfd_ext_generic.h90	(revision 9189)
+++ /branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/LBC/lbc_nfd_ext_generic.h90	(revision 9190)
@@ -11,4 +11,5 @@
 
    SUBROUTINE ROUTINE_NFD( ptab, cd_nat, psgn, kextj )
+      !!----------------------------------------------------------------------
       ARRAY_TYPE(:,:,:,:,:)                             ! array or pointer of arrays on which the boundary condition is applied
       CHARACTER(len=1) , INTENT(in   ) ::   NAT_IN(:)   ! nature of array grid-points
Index: /branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/LBC/lbc_nfd_nogather_generic.h90
===================================================================
--- /branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/LBC/lbc_nfd_nogather_generic.h90	(revision 9189)
+++ /branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/LBC/lbc_nfd_nogather_generic.h90	(revision 9190)
@@ -55,5 +55,4 @@
 
    SUBROUTINE ROUTINE_NFD( ptab, ptab2, cd_nat, psgn, kfld )
-      !
       !!----------------------------------------------------------------------
       !!
@@ -85,4 +84,5 @@
       !
       DO jf = 1, ipf                      ! Loop over the number of arrays to be processed
+         !
          SELECT CASE ( npolj )
          !
@@ -90,17 +90,18 @@
             !
             SELECT CASE ( NAT_IN(jf) )
-               CASE ( 'T' , 'W' )                         ! T-, W-point
-                  IF ( nimpp /= 1 ) THEN   ;   startloop = 1
-                  ELSE                     ;   startloop = 2
-               ENDIF
-                  !
+            !
+            CASE ( 'T' , 'W' )                         ! T-, W-point
+               IF ( nimpp /= 1 ) THEN   ;   startloop = 1
+               ELSE                     ;   startloop = 2
+               ENDIF
+               !
                DO ji = startloop, nlci
-                 ijt=jpiglo - ji - nimpp - nfiimpp(isendto(1),jpnj) + 4
-                        ARRAY_IN(ji,ijpj,:,:,jf) = SGN_IN(jf) * ARRAY2_IN(ijt,ijpj-2,:,:,jf)
-               END DO
-               IF(nimpp .eq. 1) THEN
-                 ARRAY_IN(1,ijpj,:,:,jf) = SGN_IN(jf) * ARRAY_IN(3,ijpj-2,:,:,jf)
-               ENDIF
-   
+                  ijt = jpiglo - ji - nimpp - nfiimpp(isendto(1),jpnj) + 4
+                  ARRAY_IN(ji,ijpj,:,:,jf) = SGN_IN(jf) * ARRAY2_IN(ijt,ijpj-2,:,:,jf)
+               END DO
+               IF( nimpp == 1 ) THEN
+                  ARRAY_IN(1,ijpj,:,:,jf) = SGN_IN(jf) * ARRAY_IN(3,ijpj-2,:,:,jf)
+               ENDIF
+               !
                IF( nimpp >= jpiglo/2+1 ) THEN
                   startloop = 1
@@ -110,15 +111,15 @@
                   startloop = nlci + 1
                ENDIF
-               IF(startloop <= nlci) THEN
-               DO ji = startloop, nlci
-                  ijt = jpiglo - ji - nimpp - nfiimpp(isendto(1),jpnj) + 4
-                  jia = ji + nimpp - 1
-                  ijta = jpiglo - jia + 2
-                  IF( ijta >= startloop+nimpp-1 .AND. ijta < jia ) THEN
-                             ARRAY_IN(ji,ijpjm1,:,:,jf) = SGN_IN(jf) * ARRAY_IN(ijta-nimpp+1,ijpjm1,:,:,jf)
-                  ELSE
-                             ARRAY_IN(ji,ijpjm1,:,:,jf) = SGN_IN(jf) * ARRAY2_IN(ijt,ijpjm1,:,:,jf)
-                  ENDIF
-               END DO
+               IF( startloop <= nlci ) THEN
+                  DO ji = startloop, nlci
+                     ijt  = jpiglo - ji - nimpp - nfiimpp(isendto(1),jpnj) + 4
+                     jia  = ji + nimpp - 1
+                     ijta = jpiglo - jia + 2
+                     IF( ijta >= startloop+nimpp-1 .AND. ijta < jia ) THEN
+                        ARRAY_IN(ji,ijpjm1,:,:,jf) = SGN_IN(jf) * ARRAY_IN(ijta-nimpp+1,ijpjm1,:,:,jf)
+                     ELSE
+                        ARRAY_IN(ji,ijpjm1,:,:,jf) = SGN_IN(jf) * ARRAY2_IN(ijt,ijpjm1,:,:,jf)
+                     ENDIF
+                  END DO
                ENDIF
                !
@@ -289,19 +290,19 @@
                !
             CASE ( 'I' )                                  ! ice U-V point (I-point)
-                  IF( nimpp /= 1 ) THEN
-                     startloop = 1
-                  ELSE
-                     startloop = 2
-                  ENDIF
-                  IF( nimpp + nlci - 1 /= jpiglo ) THEN
-                     endloop = nlci
-                  ELSE
-                     endloop = nlci - 1
-                  ENDIF
-                  DO ji = startloop , endloop
-                     ijt = jpiglo - ji - nimpp - nfiimpp(isendto(1),jpnj) + 4
-                     ARRAY_IN(ji,ijpj,:,:,jf) = 0.5 * (ARRAY_IN(ji,ijpjm1,:,:,jf) + SGN_IN(jf) * ARRAY2_IN(ijt,ijpjm1,:,:,jf))
-                  END DO
-                  !
+               IF( nimpp /= 1 ) THEN
+                  startloop = 1
+               ELSE
+                  startloop = 2
+               ENDIF
+               IF( nimpp + nlci - 1 /= jpiglo ) THEN
+                  endloop = nlci
+               ELSE
+                  endloop = nlci - 1
+               ENDIF
+               DO ji = startloop , endloop
+                  ijt = jpiglo - ji - nimpp - nfiimpp(isendto(1),jpnj) + 4
+                  ARRAY_IN(ji,ijpj,:,:,jf) = 0.5 * (ARRAY_IN(ji,ijpjm1,:,:,jf) + SGN_IN(jf) * ARRAY2_IN(ijt,ijpjm1,:,:,jf))
+               END DO
+               !
             END SELECT
             !
@@ -310,11 +311,11 @@
             SELECT CASE ( NAT_IN(jf))
             CASE ( 'T' , 'U' , 'V' , 'W' )                 ! T-, U-, V-, W-points
-                  ARRAY_IN(:, 1  ,:,:,jf) = 0._wp
-                  ARRAY_IN(:,ijpj,:,:,jf) = 0._wp
+               ARRAY_IN(:, 1  ,:,:,jf) = 0._wp
+               ARRAY_IN(:,ijpj,:,:,jf) = 0._wp
             CASE ( 'F' )                                   ! F-point
-                  ARRAY_IN(:,ijpj,:,:,jf) = 0._wp
+               ARRAY_IN(:,ijpj,:,:,jf) = 0._wp
             CASE ( 'I' )                                   ! ice U-V point
-                  ARRAY_IN(:, 1  ,:,:,jf) = 0._wp
-                  ARRAY_IN(:,ijpj,:,:,jf) = 0._wp
+               ARRAY_IN(:, 1  ,:,:,jf) = 0._wp
+               ARRAY_IN(:,ijpj,:,:,jf) = 0._wp
             END SELECT
             !
Index: /branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/LBC/mppini.F90
===================================================================
--- /branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/LBC/mppini.F90	(revision 9189)
+++ /branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/LBC/mppini.F90	(revision 9190)
@@ -128,6 +128,7 @@
       INTEGER ::   iiea, ijea, iiwe, ijwe     !   -       -
       INTEGER ::   iresti, irestj, iproc      !   -       -
-      INTEGER, DIMENSION(jpnij) ::   iin, ii_nono, ii_noea   ! 1D workspace
-      INTEGER, DIMENSION(jpnij) ::   ijn, ii_noso, ii_nowe   !  -     -
+      REAL(wp)::   zidom, zjdom               ! local scalars
+      INTEGER, DIMENSION(jpnij)     ::   iin, ii_nono, ii_noea   ! 1D workspace
+      INTEGER, DIMENSION(jpnij)     ::   ijn, ii_noso, ii_nowe   !  -     -
       INTEGER, DIMENSION(jpni,jpnj) ::   iimppt, ilci, ibondi, ipproc   ! 2D workspace
       INTEGER, DIMENSION(jpni,jpnj) ::   ijmppt, ilcj, ibondj, ipolj    !  -     -
@@ -135,5 +136,4 @@
       INTEGER, DIMENSION(jpni,jpnj) ::   ilej, ildj, ioso, iowe         !  -     -
       INTEGER, DIMENSION(jpiglo,jpjglo) ::   imask   ! 2D golbal domain workspace
-      REAL(wp) ::   zidom, zjdom   ! local scalars
       !!----------------------------------------------------------------------
       !
@@ -532,10 +532,10 @@
       REAL(wp), DIMENSION(jpiglo,jpjglo) ::   zbot, zbdy          ! global workspace
       REAL(wp) ::   zidom , zjdom          ! local scalars
-      NAMELIST/nambdy/ ln_bdy, nb_bdy, ln_coords_file, cn_coords_file,         &
+      NAMELIST/nambdy/ ln_bdy, nb_bdy, ln_coords_file, cn_coords_file,           &
            &             ln_mask_file, cn_mask_file, cn_dyn2d, nn_dyn2d_dta,     &
            &             cn_dyn3d, nn_dyn3d_dta, cn_tra, nn_tra_dta,             &  
            &             ln_tra_dmp, ln_dyn3d_dmp, rn_time_dmp, rn_time_dmp_out, &
-           &             cn_ice_lim, nn_ice_lim_dta,                           &
-           &             rn_ice_tem, rn_ice_sal, rn_ice_age,                 &
+           &             cn_ice_lim, nn_ice_lim_dta,                             &
+           &             rn_ice_tem, rn_ice_sal, rn_ice_age,                     &
            &             ln_vol, nn_volctl, nn_rimwidth, nb_jpk_bdy
       !!----------------------------------------------------------------------
Index: /branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk.F90
===================================================================
--- /branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk.F90	(revision 9189)
+++ /branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk.F90	(revision 9190)
@@ -273,8 +273,8 @@
          WRITE(numout,*)
          SELECT CASE( nblk )              !* Print the choice of bulk algorithm
-         CASE( np_NCAR      )   ;   WRITE(numout,*) '      ===>>   "NCAR" algorithm        (Large and Yeager 2008)'
-         CASE( np_COARE_3p0 )   ;   WRITE(numout,*) '      ===>>   "COARE 3.0" algorithm   (Fairall et al. 2003)'
-         CASE( np_COARE_3p5 )   ;   WRITE(numout,*) '      ===>>   "COARE 3.5" algorithm   (Edson et al. 2013)'
-         CASE( np_ECMWF     )   ;   WRITE(numout,*) '      ===>>   "ECMWF" algorithm       (IFS cycle 31)'
+         CASE( np_NCAR      )   ;   WRITE(numout,*) '   ==>>>   "NCAR" algorithm        (Large and Yeager 2008)'
+         CASE( np_COARE_3p0 )   ;   WRITE(numout,*) '   ==>>>   "COARE 3.0" algorithm   (Fairall et al. 2003)'
+         CASE( np_COARE_3p5 )   ;   WRITE(numout,*) '   ==>>>   "COARE 3.5" algorithm   (Edson et al. 2013)'
+         CASE( np_ECMWF     )   ;   WRITE(numout,*) '   ==>>>   "ECMWF" algorithm       (IFS cycle 31)'
          END SELECT
          !
Index: /branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90
===================================================================
--- /branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90	(revision 9189)
+++ /branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90	(revision 9190)
@@ -110,5 +110,4 @@
       READ  ( numnam_ref, namsbc, IOSTAT = ios, ERR = 901)
 901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namsbc in reference namelist', lwp )
-      !
       REWIND( numnam_cfg )          !* Namelist namsbc in configuration namelist : Parameters of the run
       READ  ( numnam_cfg, namsbc, IOSTAT = ios, ERR = 902 )
@@ -188,12 +187,12 @@
       SELECT CASE( nn_components )
       CASE( jp_iam_nemo )
-         IF(lwp) WRITE(numout,*) '   NEMO configured as a single executable (i.e. including both OPA and Surface module'
+         IF(lwp) WRITE(numout,*) '   ==>>>   NEMO configured as a single executable (i.e. including both OPA and Surface module)'
       CASE( jp_iam_opa  )
-         IF(lwp) WRITE(numout,*) '   Multi executable configuration. Here, OPA component'
+         IF(lwp) WRITE(numout,*) '   ==>>>   Multi executable configuration. Here, OPA component'
          IF( .NOT.lk_oasis )   CALL ctl_stop( 'sbc_init : OPA-SAS coupled via OASIS, but key_oasis3 disabled' )
          IF( ln_cpl        )   CALL ctl_stop( 'sbc_init : OPA-SAS coupled via OASIS, but ln_cpl = T in OPA'   )
          IF( ln_mixcpl     )   CALL ctl_stop( 'sbc_init : OPA-SAS coupled via OASIS, but ln_mixcpl = T in OPA' )
       CASE( jp_iam_sas  )
-         IF(lwp) WRITE(numout,*) '   Multi executable configuration. Here, SAS component'
+         IF(lwp) WRITE(numout,*) '   ==>>>   Multi executable configuration. Here, SAS component'
          IF( .NOT.lk_oasis )   CALL ctl_stop( 'sbc_init : OPA-SAS coupled via OASIS, but key_oasis3 disabled' )
          IF( ln_mixcpl     )   CALL ctl_stop( 'sbc_init : OPA-SAS coupled via OASIS, but ln_mixcpl = T in OPA' )
@@ -272,13 +271,13 @@
          WRITE(numout,*)
          SELECT CASE( nsbc )
-         CASE( jp_usr     )   ;   WRITE(numout,*) '   user defined formulation'
-         CASE( jp_flx     )   ;   WRITE(numout,*) '      ===>>   flux formulation'
-         CASE( jp_blk     )   ;   WRITE(numout,*) '      ===>>   bulk formulation'
-         CASE( jp_purecpl )   ;   WRITE(numout,*) '      ===>>   pure coupled formulation'
+         CASE( jp_usr     )   ;   WRITE(numout,*) '   ==>>>   user defined forcing formulation'
+         CASE( jp_flx     )   ;   WRITE(numout,*) '   ==>>>   flux formulation'
+         CASE( jp_blk     )   ;   WRITE(numout,*) '   ==>>>   bulk formulation'
+         CASE( jp_purecpl )   ;   WRITE(numout,*) '   ==>>>   pure coupled formulation'
 !!gm abusive use of jp_none ??   ===>>> need to be check and changed by adding a jp_sas parameter
-         CASE( jp_none    )   ;   WRITE(numout,*) '      ===>>   OPA coupled to SAS via oasis'
-            IF( ln_mixcpl )       WRITE(numout,*) '                  + forced-coupled mixed formulation'
+         CASE( jp_none    )   ;   WRITE(numout,*) '   ==>>>   OPA coupled to SAS via oasis'
+            IF( ln_mixcpl )       WRITE(numout,*) '               + forced-coupled mixed formulation'
          END SELECT
-         IF( ll_not_nemo  )       WRITE(numout,*) '                  + OASIS coupled SAS'
+         IF( ll_not_nemo  )       WRITE(numout,*) '               + OASIS coupled SAS'
       ENDIF
       !
Index: /branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/TRA/traadv.F90
===================================================================
--- /branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/TRA/traadv.F90	(revision 9189)
+++ /branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/TRA/traadv.F90	(revision 9190)
@@ -256,12 +256,12 @@
          WRITE(numout,*)
          SELECT CASE ( nadv )
-         CASE( np_NO_adv  )   ;   WRITE(numout,*) '      ===>>   NO T-S advection'
-         CASE( np_CEN     )   ;   WRITE(numout,*) '      ===>>   CEN      scheme is used. Horizontal order: ', nn_cen_h,   &
-            &                                                                     ' Vertical   order: ', nn_cen_v
-         CASE( np_FCT     )   ;   WRITE(numout,*) '      ===>>   FCT      scheme is used. Horizontal order: ', nn_fct_h,   &
-            &                                                                      ' Vertical   order: ', nn_fct_v
-         CASE( np_MUS     )   ;   WRITE(numout,*) '      ===>>   MUSCL    scheme is used'
-         CASE( np_UBS     )   ;   WRITE(numout,*) '      ===>>   UBS      scheme is used'
-         CASE( np_QCK     )   ;   WRITE(numout,*) '      ===>>   QUICKEST scheme is used'
+         CASE( np_NO_adv  )   ;   WRITE(numout,*) '   ==>>>   NO T-S advection'
+         CASE( np_CEN     )   ;   WRITE(numout,*) '   ==>>>   CEN      scheme is used. Horizontal order: ', nn_cen_h,   &
+            &                                                                        ' Vertical   order: ', nn_cen_v
+         CASE( np_FCT     )   ;   WRITE(numout,*) '   ==>>>   FCT      scheme is used. Horizontal order: ', nn_fct_h,   &
+            &                                                                        ' Vertical   order: ', nn_fct_v
+         CASE( np_MUS     )   ;   WRITE(numout,*) '   ==>>>   MUSCL    scheme is used'
+         CASE( np_UBS     )   ;   WRITE(numout,*) '   ==>>>   UBS      scheme is used'
+         CASE( np_QCK     )   ;   WRITE(numout,*) '   ==>>>   QUICKEST scheme is used'
          END SELECT
       ENDIF
Index: /branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_mle.F90
===================================================================
--- /branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_mle.F90	(revision 9189)
+++ /branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_mle.F90	(revision 9190)
@@ -292,9 +292,9 @@
          WRITE(numout,*)
          IF( ln_mle ) THEN
-            WRITE(numout,*) '      ===>>   Mixed Layer Eddy induced transport added to tracer advection'
+            WRITE(numout,*) '   ==>>>   Mixed Layer Eddy induced transport added to tracer advection'
             IF( nn_mle == 0 )   WRITE(numout,*) '              Fox-Kemper et al 2010 formulation'
             IF( nn_mle == 1 )   WRITE(numout,*) '              New formulation'
          ELSE
-            WRITE(numout,*) '      ===>>   Mixed Layer Eddy parametrisation NOT used'
+            WRITE(numout,*) '   ==>>>   Mixed Layer Eddy parametrisation NOT used'
          ENDIF
       ENDIF
Index: /branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/TRA/trabbc.F90
===================================================================
--- /branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/TRA/trabbc.F90	(revision 9189)
+++ /branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/TRA/trabbc.F90	(revision 9190)
@@ -160,9 +160,9 @@
          !
          CASE ( 1 )                          !* constant flux
-            IF(lwp) WRITE(numout,*) '      ===>>   constant heat flux  =   ', rn_geoflx_cst
+            IF(lwp) WRITE(numout,*) '   ==>>>   constant heat flux  =   ', rn_geoflx_cst
             qgh_trd0(:,:) = r1_rau0_rcp * rn_geoflx_cst
             !
          CASE ( 2 )                          !* variable geothermal heat flux : read the geothermal fluxes in mW/m2
-            IF(lwp) WRITE(numout,*) '      ===>>   variable geothermal heat flux'
+            IF(lwp) WRITE(numout,*) '   ==>>>   variable geothermal heat flux'
             !
             ALLOCATE( sf_qgh(1), STAT=ierror )
@@ -186,5 +186,5 @@
          !
       ELSE
-         IF(lwp) WRITE(numout,*) '      ===>>   no geothermal heat flux'
+         IF(lwp) WRITE(numout,*) '   ==>>>   no geothermal heat flux'
       ENDIF
       !
Index: /branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/TRA/traldf.F90
===================================================================
--- /branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/TRA/traldf.F90	(revision 9189)
+++ /branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/TRA/traldf.F90	(revision 9190)
@@ -185,11 +185,11 @@
          WRITE(numout,*)
          SELECT CASE( nldf )
-         CASE( np_no_ldf )   ;   WRITE(numout,*) '      ===>>   NO lateral diffusion'
-         CASE( np_lap    )   ;   WRITE(numout,*) '      ===>>   laplacian iso-level operator'
-         CASE( np_lap_i  )   ;   WRITE(numout,*) '      ===>>   Rotated laplacian operator (standard)'
-         CASE( np_lap_it )   ;   WRITE(numout,*) '      ===>>   Rotated laplacian operator (triad)'
-         CASE( np_blp    )   ;   WRITE(numout,*) '      ===>>   bilaplacian iso-level operator'
-         CASE( np_blp_i  )   ;   WRITE(numout,*) '      ===>>   Rotated bilaplacian operator (standard)'
-         CASE( np_blp_it )   ;   WRITE(numout,*) '      ===>>   Rotated bilaplacian operator (triad)'
+         CASE( np_no_ldf )   ;   WRITE(numout,*) '   ==>>>   NO lateral diffusion'
+         CASE( np_lap    )   ;   WRITE(numout,*) '   ==>>>   laplacian iso-level operator'
+         CASE( np_lap_i  )   ;   WRITE(numout,*) '   ==>>>   Rotated laplacian operator (standard)'
+         CASE( np_lap_it )   ;   WRITE(numout,*) '   ==>>>   Rotated laplacian operator (triad)'
+         CASE( np_blp    )   ;   WRITE(numout,*) '   ==>>>   bilaplacian iso-level operator'
+         CASE( np_blp_i  )   ;   WRITE(numout,*) '   ==>>>   Rotated bilaplacian operator (standard)'
+         CASE( np_blp_it )   ;   WRITE(numout,*) '   ==>>>   Rotated bilaplacian operator (triad)'
          END SELECT
       ENDIF
Index: /branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfdrg.F90
===================================================================
--- /branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfdrg.F90	(revision 9189)
+++ /branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfdrg.F90	(revision 9190)
@@ -289,6 +289,6 @@
       IF( ln_boost ) THEN           !* regional boost:   boost factor = 1 + regional boost
          IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) '   ==>>   use a regional boost read in ', TRIM(cl_file), ' file'
-         IF(lwp) WRITE(numout,*) '          using enhancement factor of ', rn_boost
+         IF(lwp) WRITE(numout,*) '   ==>>>   use a regional boost read in ', TRIM(cl_file), ' file'
+         IF(lwp) WRITE(numout,*) '           using enhancement factor of ', rn_boost
          ! cl_varname is a coefficient in [0,1] giving where to apply the regional boost
          CALL iom_open ( TRIM(cl_file), inum )
@@ -309,5 +309,5 @@
       CASE( np_NONE )            !==  No top/bottom friction  ==!   (pCdU = 0)
          IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) '   ==>>   ',TRIM(cd_topbot),' free-slip, friction set to zero'
+         IF(lwp) WRITE(numout,*) '   ==>>>   ',TRIM(cd_topbot),' free-slip, friction set to zero'
          !
          l_zdfdrg = .FALSE.         ! no time variation of the drag: set it one for all
@@ -318,5 +318,5 @@
       CASE( np_lin )             !==  linear friction  ==!   (pCdU = Cd0 * Uc0)
          IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) '   ==>>   linear ',TRIM(cd_topbot),' friction (constant coef = Cd0*Uc0 = ', rn_Cd0*rn_Uc0, ')'
+         IF(lwp) WRITE(numout,*) '   ==>>>   linear ',TRIM(cd_topbot),' friction (constant coef = Cd0*Uc0 = ', rn_Cd0*rn_Uc0, ')'
          !
          l_zdfdrg = .FALSE.         ! no time variation of the Cd*|U| : set it one for all
@@ -327,7 +327,7 @@
       CASE( np_non_lin )         !== non-linear friction  ==!   (pCd0 = Cd0 )
          IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) '   ==>>   quadratic ',TRIM(cd_topbot),' friction (propotional to module of the velocity)'
-         IF(lwp) WRITE(numout,*) '   with   Cd0 = ', rn_Cd0, ', and',   &
-            &                             ' a background velocity module of (rn_ke0)^1/2 = ', SQRT(rn_ke0), 'm/s)'
+         IF(lwp) WRITE(numout,*) '   ==>>>   quadratic ',TRIM(cd_topbot),' friction (propotional to module of the velocity)'
+         IF(lwp) WRITE(numout,*) '   with    a drag coefficient Cd0 = ', rn_Cd0, ', and'
+         IF(lwp) WRITE(numout,*) '           a background velocity module of (rn_ke0)^1/2 = ', SQRT(rn_ke0), 'm/s)'
          !
          l_zdfdrg = .TRUE.          !* Cd*|U| updated at each time-step (it depends on ocean velocity)
@@ -338,5 +338,5 @@
       CASE( np_loglayer )       !== logarithmic layer formulation of friction  ==!   (CdU = (vkarman log(z/z0))^2 |U| )
          IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) '   ==>>   quadratic ',TRIM(cd_topbot),' drag (propotional to module of the velocity)'
+         IF(lwp) WRITE(numout,*) '   ==>>>   quadratic ',TRIM(cd_topbot),' drag (propotional to module of the velocity)'
          IF(lwp) WRITE(numout,*) '   with   a logarithmic Cd0 formulation Cd0 = ( vkarman log(z/z0) )^2 ,'
          IF(lwp) WRITE(numout,*) '          a background velocity module of (rn_ke0)^1/2 = ', SQRT(pke0), 'm/s), '
Index: /branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfosm.F90
===================================================================
--- /branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfosm.F90	(revision 9189)
+++ /branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfosm.F90	(revision 9190)
@@ -164,5 +164,5 @@
 
       REAL(wp) ::   ztx, zty, zflageos, zstabl, zbuofdep,zucube      !
-      REAL(wp) ::   zalbet, zbeta, zthermal, zatt1                   !
+      REAL(wp) ::   zbeta, zthermal                                  !
       REAL(wp) ::   zehat, zeta, zhrib, zsig, zscale, zwst, zws, zwm ! Velocity scales
       REAL(wp) ::   zwsun, zwmun, zcons, zconm, zwcons, zwconm       !
Index: /branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90
===================================================================
--- /branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90	(revision 9189)
+++ /branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90	(revision 9190)
@@ -298,5 +298,6 @@
                   !                                           ! TKE Langmuir circulation source term
                   en(ji,jj,jk) = en(ji,jj,jk) + rdt * MAX(0.,1._wp - 4.*fr_i(ji,jj) ) * ( zwlc * zwlc * zwlc )   &
-                     &                              / zhlc(ji,jj) * wmask(ji,jj,jk) * tmask(ji,jj,1)
+                     &                              / zhlc(ji,jj) * wmask(ji,jj,jk)
+!!gm                     &                          / zhlc(ji,jj) * wmask(ji,jj,jk) * tmask(ji,jj,1)
                END DO
             END DO
@@ -400,5 +401,7 @@
                DO ji = fs_2, fs_jpim1   ! vector opt.
                   en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -pdepw(ji,jj,jk) / htau(ji,jj) )   &
-                     &                                 * MAX(0.,1._wp - 4.*fr_i(ji,jj) )  * wmask(ji,jj,jk) * tmask(ji,jj,1)
+                     &                                 * MAX(0.,1._wp - 4.*fr_i(ji,jj) )  * wmask(ji,jj,jk)
+!!gm                     &                             * MAX(0.,1._wp - 4.*fr_i(ji,jj) )  * wmask(ji,jj,jk) * tmask(ji,jj,1)
+!!gm    multiplication by surface tmask useless (already includes in MAX( 0, 1-4*fr_i )
                END DO
             END DO
@@ -409,5 +412,6 @@
                jk = nmln(ji,jj)
                en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -pdepw(ji,jj,jk) / htau(ji,jj) )   &
-                  &                                 * MAX(0.,1._wp - 4.*fr_i(ji,jj) )  * wmask(ji,jj,jk) * tmask(ji,jj,1)
+                  &                                 * MAX(0.,1._wp - 4.*fr_i(ji,jj) )  * wmask(ji,jj,jk)
+!!gm                  &                             * MAX(0.,1._wp - 4.*fr_i(ji,jj) )  * wmask(ji,jj,jk) * tmask(ji,jj,1)
             END DO
          END DO
@@ -422,5 +426,6 @@
                   zdif = rhftau_scl * MAX( 0._wp, zdif + rhftau_add )  ! apply some modifications...
                   en(ji,jj,jk) = en(ji,jj,jk) + zbbrau * zdif * EXP( -pdepw(ji,jj,jk) / htau(ji,jj) )   &
-                     &                        * MAX(0.,1._wp - 4.*fr_i(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1)
+                     &                        * MAX(0.,1._wp - 4.*fr_i(ji,jj) ) * wmask(ji,jj,jk)
+!!gm                     &                    * MAX(0.,1._wp - 4.*fr_i(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1)
                END DO
             END DO
@@ -685,5 +690,5 @@
          ENDIF
          WRITE(numout,*)
-         WRITE(numout,*) '   ==>>   critical Richardson nb with your parameters  ri_cri = ', ri_cri
+         WRITE(numout,*) '   ==>>>   critical Richardson nb with your parameters  ri_cri = ', ri_cri
          WRITE(numout,*)
       ENDIF
@@ -692,8 +697,8 @@
          rn_emin  = 1.e-10_wp             ! specific values of rn_emin & rmxl_min are used
          rmxl_min = 1.e-03_wp             ! associated avt minimum = molecular salt diffusivity (10^-9 m2/s)
-         IF(lwp) WRITE(numout,*) '   ==>>   Internal wave-driven mixing case:   force   rn_emin = 1.e-10 and rmxl_min = 1.e-3'
+         IF(lwp) WRITE(numout,*) '   ==>>>   Internal wave-driven mixing case:   force   rn_emin = 1.e-10 and rmxl_min = 1.e-3'
       ELSE                          ! standard case : associated avt minimum = molecular viscosity (10^-6 m2/s)
          rmxl_min = 1.e-6_wp / ( rn_ediff * SQRT( rn_emin ) )    ! resulting minimum length to recover molecular viscosity
-         IF(lwp) WRITE(numout,*) '   ==>>   minimum mixing length with your parameters rmxl_min = ', rmxl_min
+         IF(lwp) WRITE(numout,*) '   ==>>>   minimum mixing length with your parameters rmxl_min = ', rmxl_min
       ENDIF
       !
@@ -709,5 +714,5 @@
       IF( ln_mxl0 ) THEN
          IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) '   ==>>   use a surface mixing length = F(stress) :   set rn_mxl0 = rmxl_min'
+         IF(lwp) WRITE(numout,*) '   ==>>>   use a surface mixing length = F(stress) :   set rn_mxl0 = rmxl_min'
          rn_mxl0 = rmxl_min
       ENDIF
@@ -764,5 +769,5 @@
             ELSE                                          ! start TKE from rest
                IF(lwp) WRITE(numout,*)
-               IF(lwp) WRITE(numout,*) '   ==>>   previous run without TKE scheme, set en to background values'
+               IF(lwp) WRITE(numout,*) '   ==>>>   previous run without TKE scheme, set en to background values'
                en   (:,:,:) = rn_emin * wmask(:,:,:)
                dissl(:,:,:) = 1.e-12_wp
@@ -771,5 +776,5 @@
          ELSE                                   !* Start from rest
             IF(lwp) WRITE(numout,*)
-            IF(lwp) WRITE(numout,*) '   ==>>   start from rest: set en to the background value'
+            IF(lwp) WRITE(numout,*) '   ==>>>   start from rest: set en to the background value'
             en   (:,:,:) = rn_emin * wmask(:,:,:)
             dissl(:,:,:) = 1.e-12_wp
