Index: branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/CONFIG/SHARED/namelist_ref
===================================================================
--- branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/CONFIG/SHARED/namelist_ref	(revision 7162)
+++ branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/CONFIG/SHARED/namelist_ref	(revision 7163)
@@ -270,5 +270,5 @@
    nn_ice      = 2         !  =0 no ice boundary condition   ,
                            !  =1 use observed ice-cover      ,
-                           !  =2 ice-model used                         ("key_lim3", "key_lim2", or "key_cice")
+                           !  =2 to 4 :  ice-model used (LIM2, LIM3 or CICE)                         ("key_lim3", "key_lim2", or "key_cice")
    nn_ice_embd = 1         !  =0 levitating ice (no mass exchange, concentration/dilution effect)
                            !  =1 levitating ice with mass and salt exchange but no presure effect
Index: branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90
===================================================================
--- branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90	(revision 7162)
+++ branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90	(revision 7163)
@@ -110,5 +110,5 @@
 
       IF(lwp) WRITE(numout,*)
-      IF(lwp) WRITE(numout,*) 'lim_istate: sea-ice initialization '
+      IF(lwp) WRITE(numout,*) 'lim_istate : sea-ice initialization '
       IF(lwp) WRITE(numout,*) '~~~~~~~~~~ '
 
Index: branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90
===================================================================
--- branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90	(revision 7162)
+++ branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90	(revision 7163)
@@ -969,17 +969,18 @@
       IF (lwp) THEN                          ! control print
          WRITE(numout,*)
-         WRITE(numout,*)' lim_itd_me_init : ice parameters for mechanical ice redistribution '
-         WRITE(numout,*)' ~~~~~~~~~~~~~~~'
-         WRITE(numout,*)'   Fraction of shear energy contributing to ridging        rn_cs       = ', rn_cs 
-         WRITE(numout,*)'   Fraction of snow volume conserved during ridging        rn_fsnowrdg = ', rn_fsnowrdg 
-         WRITE(numout,*)'   Fraction of snow volume conserved during ridging        rn_fsnowrft = ', rn_fsnowrft 
-         WRITE(numout,*)'   Fraction of total ice coverage contributing to ridging  rn_gstar    = ', rn_gstar
-         WRITE(numout,*)'   Equivalent to G* for an exponential part function       rn_astar    = ', rn_astar
-         WRITE(numout,*)'   Quantity playing a role in max ridged ice thickness     rn_hstar    = ', rn_hstar
-         WRITE(numout,*)'   Rafting of ice sheets or not                            ln_rafting  = ', ln_rafting
-         WRITE(numout,*)'   Parmeter thickness (threshold between ridge-raft)       rn_hraft    = ', rn_hraft
-         WRITE(numout,*)'   Rafting hyperbolic tangent coefficient                  rn_craft    = ', rn_craft  
-         WRITE(numout,*)'   Initial porosity of ridges                              rn_por_rdg  = ', rn_por_rdg
-         WRITE(numout,*)'   Switch for part. function (0) linear (1) exponential    nn_partfun  = ', nn_partfun
+         WRITE(numout,*) ' lim_itd_me_init : ice parameters for mechanical ice redistribution '
+         WRITE(numout,*) ' ~~~~~~~~~~~~~~~'
+         WRITE(numout,*) '   Namelist namiceitdme :'
+         WRITE(numout,*) '      Fraction of shear energy contributing to ridging        rn_cs       = ', rn_cs 
+         WRITE(numout,*) '      Fraction of snow volume conserved during ridging        rn_fsnowrdg = ', rn_fsnowrdg 
+         WRITE(numout,*) '      Fraction of snow volume conserved during ridging        rn_fsnowrft = ', rn_fsnowrft 
+         WRITE(numout,*) '      Fraction of total ice coverage contributing to ridging  rn_gstar    = ', rn_gstar
+         WRITE(numout,*) '      Equivalent to G* for an exponential part function       rn_astar    = ', rn_astar
+         WRITE(numout,*) '      Quantity playing a role in max ridged ice thickness     rn_hstar    = ', rn_hstar
+         WRITE(numout,*) '      Rafting of ice sheets or not                            ln_rafting  = ', ln_rafting
+         WRITE(numout,*) '      Parmeter thickness (threshold between ridge-raft)       rn_hraft    = ', rn_hraft
+         WRITE(numout,*) '      Rafting hyperbolic tangent coefficient                  rn_craft    = ', rn_craft  
+         WRITE(numout,*) '      Initial porosity of ridges                              rn_por_rdg  = ', rn_por_rdg
+         WRITE(numout,*) '      Switch for part. function (0) linear (1) exponential    nn_partfun  = ', nn_partfun
       ENDIF
       !
Index: branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90
===================================================================
--- branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90	(revision 7162)
+++ branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90	(revision 7163)
@@ -621,6 +621,6 @@
       IF(lwp) THEN
          WRITE(numout,*)
-         WRITE(numout,*) 'lim_thd : Ice Thermodynamics'
-         WRITE(numout,*) '~~~~~~~'
+         WRITE(numout,*) 'lim_thd_init : Ice Thermodynamics initialization'
+         WRITE(numout,*) '~~~~~~~~~~~~'
       ENDIF
       !
@@ -634,12 +634,5 @@
       IF(lwm) WRITE ( numoni, namicethd )
       !
-      IF ( ( jpl > 1 ) .AND. ( nn_monocat == 1 ) ) THEN 
-         nn_monocat = 0
-         IF(lwp) WRITE(numout, *) '   nn_monocat must be 0 in multi-category case '
-      ENDIF
-
-      !
       IF(lwp) THEN                          ! control print
-         WRITE(numout,*)
          WRITE(numout,*)'   Namelist of ice parameters for ice thermodynamic computation '
          WRITE(numout,*)'      ice thick. for lateral accretion                        rn_hnewice   = ', rn_hnewice
@@ -659,4 +652,9 @@
          WRITE(numout,*)'      iterate the surface non-solar flux (T) or not (F)       ln_it_qnsice = ', ln_it_qnsice
       ENDIF
+      IF( jpl > 1 .AND. nn_monocat == 1 ) THEN 
+         nn_monocat = 0
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) '   nn_monocat forced to 0 as jpl>1, i.e. multi-category case is chosen'
+      ENDIF
       !
    END SUBROUTINE lim_thd_init
Index: branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/LIM_SRC_3/limthd_sal.F90
===================================================================
--- branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/LIM_SRC_3/limthd_sal.F90	(revision 7162)
+++ branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/LIM_SRC_3/limthd_sal.F90	(revision 7163)
@@ -144,12 +144,13 @@
          WRITE(numout,*) 'lim_thd_sal_init : Ice parameters for salinity '
          WRITE(numout,*) '~~~~~~~~~~~~~~~~'
-         WRITE(numout,*) '   switch for salinity nn_icesal        = ', nn_icesal
-         WRITE(numout,*) '   bulk salinity value if nn_icesal = 1 = ', rn_icesal
-         WRITE(numout,*) '   restoring salinity for GD            = ', rn_sal_gd
-         WRITE(numout,*) '   restoring time for GD                = ', rn_time_gd
-         WRITE(numout,*) '   restoring salinity for flushing      = ', rn_sal_fl
-         WRITE(numout,*) '   restoring time for flushing          = ', rn_time_fl
-         WRITE(numout,*) '   Maximum tolerated ice salinity       = ', rn_simax
-         WRITE(numout,*) '   Minimum tolerated ice salinity       = ', rn_simin
+         WRITE(numout,*) '   Namelist namicesal :'
+         WRITE(numout,*) '      switch for salinity nn_icesal        = ', nn_icesal
+         WRITE(numout,*) '      bulk salinity value if nn_icesal = 1 = ', rn_icesal
+         WRITE(numout,*) '      restoring salinity for GD            = ', rn_sal_gd
+         WRITE(numout,*) '      restoring time for GD                = ', rn_time_gd
+         WRITE(numout,*) '      restoring salinity for flushing      = ', rn_sal_fl
+         WRITE(numout,*) '      restoring time for flushing          = ', rn_time_fl
+         WRITE(numout,*) '      Maximum tolerated ice salinity       = ', rn_simax
+         WRITE(numout,*) '      Minimum tolerated ice salinity       = ', rn_simin
       ENDIF
       !
Index: branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/DIA/diahsb.F90
===================================================================
--- branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/DIA/diahsb.F90	(revision 7162)
+++ branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/DIA/diahsb.F90	(revision 7163)
@@ -230,98 +230,98 @@
 
    SUBROUTINE dia_hsb_rst( kt, cdrw )
-     !!---------------------------------------------------------------------
-     !!                   ***  ROUTINE limdia_rst  ***
-     !!                     
-     !! ** Purpose :   Read or write DIA file in restart file
-     !!
-     !! ** Method  :   use of IOM library
-     !!----------------------------------------------------------------------
-     INTEGER         , INTENT(in) ::   kt     ! ocean time-step
-     CHARACTER(len=*), INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
-     !
-     INTEGER ::   ji, jj, jk   ! dummy loop indices
-     INTEGER ::   id1          ! local integers
-     !!----------------------------------------------------------------------
-     !
-     IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise 
-        IF( ln_rstart ) THEN                   !* Read the restart file
-           !id1 = iom_varid( numror, 'frc_vol'  , ldstop = .FALSE. )
-           !
-           IF(lwp) WRITE(numout,*) '~~~~~~~'
-           IF(lwp) WRITE(numout,*) ' dia_hsb_rst at it= ', kt,' date= ', ndastp
-           IF(lwp) WRITE(numout,*) '~~~~~~~'
-           CALL iom_get( numror, 'frc_v', frc_v )
-           CALL iom_get( numror, 'frc_t', frc_t )
-           CALL iom_get( numror, 'frc_s', frc_s )
-           IF( ln_linssh ) THEN
-              CALL iom_get( numror, 'frc_wn_t', frc_wn_t )
-              CALL iom_get( numror, 'frc_wn_s', frc_wn_s )
-           ENDIF
-           CALL iom_get( numror, jpdom_autoglo, 'surf_ini', surf_ini ) ! ice sheet coupling
-           CALL iom_get( numror, jpdom_autoglo, 'ssh_ini', ssh_ini )
-           CALL iom_get( numror, jpdom_autoglo, 'e3t_ini', e3t_ini )
-           CALL iom_get( numror, jpdom_autoglo, 'hc_loc_ini', hc_loc_ini )
-           CALL iom_get( numror, jpdom_autoglo, 'sc_loc_ini', sc_loc_ini )
-           IF( ln_linssh ) THEN
-              CALL iom_get( numror, jpdom_autoglo, 'ssh_hc_loc_ini', ssh_hc_loc_ini )
-              CALL iom_get( numror, jpdom_autoglo, 'ssh_sc_loc_ini', ssh_sc_loc_ini )
-           ENDIF
-       ELSE
-          IF(lwp) WRITE(numout,*) '~~~~~~~'
-          IF(lwp) WRITE(numout,*) ' dia_hsb at initial state '
-          IF(lwp) WRITE(numout,*) '~~~~~~~'
-          surf_ini(:,:) = e1e2t(:,:) * tmask_i(:,:)         ! initial ocean surface
-          ssh_ini(:,:) = sshn(:,:)                          ! initial ssh
-          DO jk = 1, jpk
-             ! if ice sheet/oceqn coupling, need to mask ini variables here (mask could change at the next NEMO instance).
-             e3t_ini   (:,:,jk) = e3t_n(:,:,jk)                      * tmask(:,:,jk)  ! initial vertical scale factors
-             hc_loc_ini(:,:,jk) = tsn(:,:,jk,jp_tem) * e3t_n(:,:,jk) * tmask(:,:,jk)  ! initial heat content
-             sc_loc_ini(:,:,jk) = tsn(:,:,jk,jp_sal) * e3t_n(:,:,jk) * tmask(:,:,jk)  ! initial salt content
-          END DO
-          frc_v = 0._wp                                           ! volume       trend due to forcing
-          frc_t = 0._wp                                           ! heat content   -    -   -    -   
-          frc_s = 0._wp                                           ! salt content   -    -   -    -        
-          IF( ln_linssh ) THEN
-             IF ( ln_isfcav ) THEN
-                DO ji=1,jpi
-                   DO jj=1,jpj
-                      ssh_hc_loc_ini(ji,jj) = tsn(ji,jj,mikt(ji,jj),jp_tem) * sshn(ji,jj)   ! initial heat content in ssh
-                      ssh_sc_loc_ini(ji,jj) = tsn(ji,jj,mikt(ji,jj),jp_sal) * sshn(ji,jj)   ! initial salt content in ssh
-                   ENDDO
-                ENDDO
-             ELSE
-                ssh_hc_loc_ini(:,:) = tsn(:,:,1,jp_tem) * sshn(:,:)   ! initial heat content in ssh
-                ssh_sc_loc_ini(:,:) = tsn(:,:,1,jp_sal) * sshn(:,:)   ! initial salt content in ssh
-             END IF
-             frc_wn_t = 0._wp                                       ! initial heat content misfit due to free surface
-             frc_wn_s = 0._wp                                       ! initial salt content misfit due to free surface
-          ENDIF
-       ENDIF
-
-     ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
-        !                                   ! -------------------
-        IF(lwp) WRITE(numout,*) '~~~~~~~'
-        IF(lwp) WRITE(numout,*) ' dia_hsb_rst at it= ', kt,' date= ', ndastp
-        IF(lwp) WRITE(numout,*) '~~~~~~~'
-
-        CALL iom_rstput( kt, nitrst, numrow, 'frc_v'   , frc_v     )
-        CALL iom_rstput( kt, nitrst, numrow, 'frc_t'   , frc_t     )
-        CALL iom_rstput( kt, nitrst, numrow, 'frc_s'   , frc_s     )
-        IF( ln_linssh ) THEN
-           CALL iom_rstput( kt, nitrst, numrow, 'frc_wn_t', frc_wn_t )
-           CALL iom_rstput( kt, nitrst, numrow, 'frc_wn_s', frc_wn_s )
-        ENDIF
-        CALL iom_rstput( kt, nitrst, numrow, 'surf_ini', surf_ini )      ! ice sheet coupling
-        CALL iom_rstput( kt, nitrst, numrow, 'ssh_ini', ssh_ini )
-        CALL iom_rstput( kt, nitrst, numrow, 'e3t_ini', e3t_ini )
-        CALL iom_rstput( kt, nitrst, numrow, 'hc_loc_ini', hc_loc_ini )
-        CALL iom_rstput( kt, nitrst, numrow, 'sc_loc_ini', sc_loc_ini )
-        IF( ln_linssh ) THEN
-           CALL iom_rstput( kt, nitrst, numrow, 'ssh_hc_loc_ini', ssh_hc_loc_ini )
-           CALL iom_rstput( kt, nitrst, numrow, 'ssh_sc_loc_ini', ssh_sc_loc_ini )
-        ENDIF
-        !
-     ENDIF
-     !
+      !!---------------------------------------------------------------------
+      !!                   ***  ROUTINE limdia_rst  ***
+      !!                     
+      !! ** Purpose :   Read or write DIA file in restart file
+      !! 
+      !! ** Method  :   use of IOM library
+      !!----------------------------------------------------------------------
+      INTEGER         , INTENT(in) ::   kt     ! ocean time-step
+      CHARACTER(len=*), INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
+      !
+      INTEGER ::   ji, jj, jk   ! dummy loop indices
+      INTEGER ::   id1          ! local integers
+      !!----------------------------------------------------------------------
+      !
+      IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise 
+         IF( ln_rstart ) THEN                   !* Read the restart file
+            !id1 = iom_varid( numror, 'frc_vol'  , ldstop = .FALSE. )
+            !
+            IF(lwp) WRITE(numout,*)
+            IF(lwp) WRITE(numout,*) '   dia_hsb_rst : read restart at it= ', kt,' date= ', ndastp
+            IF(lwp) WRITE(numout,*) '   ~~~~~~~~~~~'
+            CALL iom_get( numror, 'frc_v', frc_v )
+            CALL iom_get( numror, 'frc_t', frc_t )
+            CALL iom_get( numror, 'frc_s', frc_s )
+            IF( ln_linssh ) THEN
+               CALL iom_get( numror, 'frc_wn_t', frc_wn_t )
+               CALL iom_get( numror, 'frc_wn_s', frc_wn_s )
+            ENDIF
+            CALL iom_get( numror, jpdom_autoglo, 'surf_ini', surf_ini ) ! ice sheet coupling
+            CALL iom_get( numror, jpdom_autoglo, 'ssh_ini', ssh_ini )
+            CALL iom_get( numror, jpdom_autoglo, 'e3t_ini', e3t_ini )
+            CALL iom_get( numror, jpdom_autoglo, 'hc_loc_ini', hc_loc_ini )
+            CALL iom_get( numror, jpdom_autoglo, 'sc_loc_ini', sc_loc_ini )
+            IF( ln_linssh ) THEN
+               CALL iom_get( numror, jpdom_autoglo, 'ssh_hc_loc_ini', ssh_hc_loc_ini )
+               CALL iom_get( numror, jpdom_autoglo, 'ssh_sc_loc_ini', ssh_sc_loc_ini )
+            ENDIF
+         ELSE
+            IF(lwp) WRITE(numout,*)
+            IF(lwp) WRITE(numout,*) '   dia_hsb_rst : no restart, set value at initial state '
+            IF(lwp) WRITE(numout,*) '   ~~~~~~~~~~~'
+            surf_ini(:,:) = e1e2t(:,:) * tmask_i(:,:)         ! initial ocean surface
+            ssh_ini(:,:) = sshn(:,:)                          ! initial ssh
+            DO jk = 1, jpk
+               ! if ice sheet/oceqn coupling, need to mask ini variables here (mask could change at the next NEMO instance).
+               e3t_ini   (:,:,jk) = e3t_n(:,:,jk)                      * tmask(:,:,jk)  ! initial vertical scale factors
+               hc_loc_ini(:,:,jk) = tsn(:,:,jk,jp_tem) * e3t_n(:,:,jk) * tmask(:,:,jk)  ! initial heat content
+               sc_loc_ini(:,:,jk) = tsn(:,:,jk,jp_sal) * e3t_n(:,:,jk) * tmask(:,:,jk)  ! initial salt content
+            END DO
+            frc_v = 0._wp                                           ! volume       trend due to forcing
+            frc_t = 0._wp                                           ! heat content   -    -   -    -   
+            frc_s = 0._wp                                           ! salt content   -    -   -    -        
+            IF( ln_linssh ) THEN
+               IF( ln_isfcav ) THEN
+                  DO ji=1,jpi
+                     DO jj=1,jpj
+                        ssh_hc_loc_ini(ji,jj) = tsn(ji,jj,mikt(ji,jj),jp_tem) * sshn(ji,jj)   ! initial heat content in ssh
+                        ssh_sc_loc_ini(ji,jj) = tsn(ji,jj,mikt(ji,jj),jp_sal) * sshn(ji,jj)   ! initial salt content in ssh
+                     END DO
+                  END DO
+               ELSE
+                  ssh_hc_loc_ini(:,:) = tsn(:,:,1,jp_tem) * sshn(:,:)   ! initial heat content in ssh
+                  ssh_sc_loc_ini(:,:) = tsn(:,:,1,jp_sal) * sshn(:,:)   ! initial salt content in ssh
+               END IF
+               frc_wn_t = 0._wp                                       ! initial heat content misfit due to free surface
+               frc_wn_s = 0._wp                                       ! initial salt content misfit due to free surface
+            ENDIF
+         ENDIF
+         !
+      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
+         !                                   ! -------------------
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) '   dia_hsb_rst : write restart at it= ', kt,' date= ', ndastp
+         IF(lwp) WRITE(numout,*) '   ~~~~~~~~~~~'
+
+         CALL iom_rstput( kt, nitrst, numrow, 'frc_v'   , frc_v     )
+         CALL iom_rstput( kt, nitrst, numrow, 'frc_t'   , frc_t     )
+         CALL iom_rstput( kt, nitrst, numrow, 'frc_s'   , frc_s     )
+         IF( ln_linssh ) THEN
+            CALL iom_rstput( kt, nitrst, numrow, 'frc_wn_t', frc_wn_t )
+            CALL iom_rstput( kt, nitrst, numrow, 'frc_wn_s', frc_wn_s )
+         ENDIF
+         CALL iom_rstput( kt, nitrst, numrow, 'surf_ini', surf_ini )      ! ice sheet coupling
+         CALL iom_rstput( kt, nitrst, numrow, 'ssh_ini', ssh_ini )
+         CALL iom_rstput( kt, nitrst, numrow, 'e3t_ini', e3t_ini )
+         CALL iom_rstput( kt, nitrst, numrow, 'hc_loc_ini', hc_loc_ini )
+         CALL iom_rstput( kt, nitrst, numrow, 'sc_loc_ini', sc_loc_ini )
+         IF( ln_linssh ) THEN
+            CALL iom_rstput( kt, nitrst, numrow, 'ssh_hc_loc_ini', ssh_hc_loc_ini )
+            CALL iom_rstput( kt, nitrst, numrow, 'ssh_sc_loc_ini', ssh_sc_loc_ini )
+         ENDIF
+         !
+      ENDIF
+      !
    END SUBROUTINE dia_hsb_rst
 
@@ -342,14 +342,8 @@
       INTEGER ::   ierror   ! local integer
       INTEGER ::   ios
-      !
+      !!
       NAMELIST/namhsb/ ln_diahsb
       !!----------------------------------------------------------------------
-
-      IF(lwp) THEN
-         WRITE(numout,*)
-         WRITE(numout,*) 'dia_hsb_init : check the heat and salt budgets'
-         WRITE(numout,*) '~~~~~~~~ '
-      ENDIF
-
+      !
       REWIND( numnam_ref )              ! Namelist namhsb in reference namelist
       READ  ( numnam_ref, namhsb, IOSTAT = ios, ERR = 901)
@@ -368,12 +362,7 @@
          WRITE(numout,*) '   Namelist namhsb : set hsb parameters'
          WRITE(numout,*) '      Switch for hsb diagnostic (T) or not (F)  ln_diahsb  = ', ln_diahsb
-         WRITE(numout,*)
       ENDIF
 
       IF( .NOT. ln_diahsb )   RETURN
-         !      IF( .NOT. lk_mpp_rep ) &
-         !        CALL ctl_stop (' Your global mpp_sum if performed in single precision - 64 bits -', &
-         !             &         ' whereas the global sum to be precise must be done in double precision ',&
-         !             &         ' please add key_mpp_rep')
 
       ! ------------------- !
@@ -383,10 +372,10 @@
          &      e3t_ini(jpi,jpj,jpk), surf(jpi,jpj),  ssh_ini(jpi,jpj), STAT=ierror  )
       IF( ierror > 0 ) THEN
-         CALL ctl_stop( 'dia_hsb: unable to allocate hc_loc_ini' )   ;   RETURN
+         CALL ctl_stop( 'dia_hsb_init: unable to allocate hc_loc_ini' )   ;   RETURN
       ENDIF
 
       IF( ln_linssh )   ALLOCATE( ssh_hc_loc_ini(jpi,jpj), ssh_sc_loc_ini(jpi,jpj),STAT=ierror )
       IF( ierror > 0 ) THEN
-         CALL ctl_stop( 'dia_hsb: unable to allocate hc_loc_ini' )   ;   RETURN
+         CALL ctl_stop( 'dia_hsb_init: unable to allocate hc_loc_ini' )   ;   RETURN
       ENDIF
 
@@ -394,10 +383,10 @@
       ! 2 - Time independant variables and file opening !
       ! ----------------------------------------------- !
-      IF(lwp) WRITE(numout,*) "dia_hsb: heat salt volume budgets activated"
-      IF(lwp) WRITE(numout,*) '~~~~~~~'
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*) "   heat salt volume budgets activated"
       surf(:,:) = e1t(:,:) * e2t(:,:) * tmask_i(:,:)      ! masked surface grid cell area
       surf_tot  = glob_sum( surf(:,:) )                                       ! total ocean surface area
 
-      IF( lk_bdy ) CALL ctl_warn( 'dia_hsb does not take open boundary fluxes into account' )         
+      IF( lk_bdy ) CALL ctl_warn( 'dia_hsb_init: heat/salt budget does not consider open boundary fluxes' )         
       !
       ! ---------------------------------- !
Index: branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/DOM/dtatsd.F90
===================================================================
--- branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/DOM/dtatsd.F90	(revision 7162)
+++ branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/DOM/dtatsd.F90	(revision 7163)
@@ -53,6 +53,6 @@
       LOGICAL, INTENT(in), OPTIONAL ::   ld_tradmp   ! force the initialization when tradp is used
       !
-      INTEGER ::   ierr0, ierr1, ierr2, ierr3   ! temporary integers
-      !
+      INTEGER ::   ios, ierr0, ierr1, ierr2, ierr3   ! local integers
+      !!
       CHARACTER(len=100)            ::   cn_dir          ! Root directory for location of ssr files
       TYPE(FLD_N), DIMENSION( jpts) ::   slf_i           ! array of namelist informations on the fields to read
@@ -60,5 +60,4 @@
       !!
       NAMELIST/namtsd/   ln_tsd_init, ln_tsd_tradmp, cn_dir, sn_tem, sn_sal
-      INTEGER  ::   ios
       !!----------------------------------------------------------------------
       !
@@ -117,5 +116,5 @@
          !                         ! fill sf_tsd with sn_tem & sn_sal and control print
          slf_i(jp_tem) = sn_tem   ;   slf_i(jp_sal) = sn_sal
-         CALL fld_fill( sf_tsd, slf_i, cn_dir, 'dta_tsd', 'Temperature & Salinity data', 'namtsd' )
+         CALL fld_fill( sf_tsd, slf_i, cn_dir, 'dta_tsd', 'Temperature & Salinity data', 'namtsd', no_print )
          !
       ENDIF
Index: branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/DOM/phycst.F90
===================================================================
--- branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/DOM/phycst.F90	(revision 7162)
+++ branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/DOM/phycst.F90	(revision 7163)
@@ -106,6 +106,6 @@
 
       IF(lwp) WRITE(numout,*)
-      IF(lwp) WRITE(numout,*) ' phy_cst : initialization of ocean parameters and constants'
-      IF(lwp) WRITE(numout,*) ' ~~~~~~~'
+      IF(lwp) WRITE(numout,*) 'phy_cst : initialization of ocean parameters and constants'
+      IF(lwp) WRITE(numout,*) '~~~~~~~'
 
       ! Ocean Parameters
Index: branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/DYN/dynadv.F90
===================================================================
--- branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/DYN/dynadv.F90	(revision 7162)
+++ branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/DYN/dynadv.F90	(revision 7163)
@@ -106,11 +106,11 @@
          WRITE(numout,*)
          WRITE(numout,*) 'dyn_adv_init : choice/control of the momentum advection scheme'
-         WRITE(numout,*) '~~~~~~~~~~~'
-         WRITE(numout,*) '       Namelist namdyn_adv : chose a advection formulation & scheme for momentum'
-         WRITE(numout,*) '          Vector/flux form (T/F)                           ln_dynadv_vec  = ', ln_dynadv_vec
-         WRITE(numout,*) '          = 0 standard scheme  ; =1 Hollingsworth scheme   nn_dynkeg      = ', nn_dynkeg
-         WRITE(numout,*) '          2nd order centred advection scheme               ln_dynadv_cen2 = ', ln_dynadv_cen2
-         WRITE(numout,*) '          3rd order UBS advection scheme                   ln_dynadv_ubs  = ', ln_dynadv_ubs
-         WRITE(numout,*) '          Sub timestepping of vertical advection           ln_dynzad_zts  = ', ln_dynzad_zts
+         WRITE(numout,*) '~~~~~~~~~~~~'
+         WRITE(numout,*) '   Namelist namdyn_adv : chose a advection formulation & scheme for momentum'
+         WRITE(numout,*) '      Vector/flux form (T/F)                           ln_dynadv_vec  = ', ln_dynadv_vec
+         WRITE(numout,*) '      = 0 standard scheme  ; =1 Hollingsworth scheme   nn_dynkeg      = ', nn_dynkeg
+         WRITE(numout,*) '      2nd order centred advection scheme               ln_dynadv_cen2 = ', ln_dynadv_cen2
+         WRITE(numout,*) '      3rd order UBS advection scheme                   ln_dynadv_ubs  = ', ln_dynadv_ubs
+         WRITE(numout,*) '      Sub timestepping of vertical advection           ln_dynzad_zts  = ', ln_dynzad_zts
       ENDIF
 
@@ -134,12 +134,12 @@
       IF(lwp) THEN                    ! Print the choice
          WRITE(numout,*)
-         IF( nadv ==  0 )   WRITE(numout,*) '         vector form : keg + zad + vor is used' 
-         IF( nadv ==  1 )   WRITE(numout,*) '         vector form : keg + zad_zts + vor is used'
+         IF( nadv ==  0 )   WRITE(numout,*) '      ===>>   vector form : keg + zad + vor is used' 
+         IF( nadv ==  1 )   WRITE(numout,*) '      ===>>   vector form : keg + zad_zts + vor is used'
          IF( nadv ==  0 .OR. nadv ==  1 ) THEN
-            IF( nn_dynkeg == nkeg_C2  )   WRITE(numout,*) 'with Centered standard keg scheme'
-            IF( nn_dynkeg == nkeg_HW  )   WRITE(numout,*) 'with Hollingsworth keg scheme'
+            IF( nn_dynkeg == nkeg_C2  )   WRITE(numout,*) '              with Centered standard keg scheme'
+            IF( nn_dynkeg == nkeg_HW  )   WRITE(numout,*) '              with Hollingsworth keg scheme'
          ENDIF
-         IF( nadv ==  2 )   WRITE(numout,*) '         flux form   : 2nd order scheme is used'
-         IF( nadv ==  3 )   WRITE(numout,*) '         flux form   : UBS       scheme is used'
+         IF( nadv ==  2 )   WRITE(numout,*) '      ===>>   flux form   : 2nd order scheme is used'
+         IF( nadv ==  3 )   WRITE(numout,*) '      ===>>   flux form   : UBS       scheme is used'
       ENDIF
       !
Index: branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf.F90
===================================================================
--- branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf.F90	(revision 7162)
+++ branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf.F90	(revision 7163)
@@ -110,11 +110,11 @@
          WRITE(numout,*)
          WRITE(numout,*) 'dyn_ldf_init : Choice of the lateral diffusive operator on dynamics'
-         WRITE(numout,*) '~~~~~~~~~~~'
-         WRITE(numout,*) '       Namelist nam_dynldf : set lateral mixing parameters (type, direction, coefficients)'
-         WRITE(numout,*) '          laplacian operator          ln_dynldf_lap = ', ln_dynldf_lap
-         WRITE(numout,*) '          bilaplacian operator        ln_dynldf_blp = ', ln_dynldf_blp
-         WRITE(numout,*) '          iso-level                   ln_dynldf_lev = ', ln_dynldf_lev
-         WRITE(numout,*) '          horizontal (geopotential)   ln_dynldf_hor = ', ln_dynldf_hor
-         WRITE(numout,*) '          iso-neutral                 ln_dynldf_iso = ', ln_dynldf_iso
+         WRITE(numout,*) '~~~~~~~~~~~~'
+         WRITE(numout,*) '   Namelist nam_dynldf : set lateral mixing parameters (type, direction, coefficients)'
+         WRITE(numout,*) '      laplacian operator          ln_dynldf_lap = ', ln_dynldf_lap
+         WRITE(numout,*) '      bilaplacian operator        ln_dynldf_blp = ', ln_dynldf_blp
+         WRITE(numout,*) '      iso-level                   ln_dynldf_lev = ', ln_dynldf_lev
+         WRITE(numout,*) '      horizontal (geopotential)   ln_dynldf_hor = ', ln_dynldf_hor
+         WRITE(numout,*) '      iso-neutral                 ln_dynldf_iso = ', ln_dynldf_iso
       ENDIF
       !                                   ! use of lateral operator or not
@@ -180,8 +180,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/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg.F90
===================================================================
--- branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg.F90	(revision 7162)
+++ branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg.F90	(revision 7163)
@@ -216,7 +216,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/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/DYN/dynvor.F90
===================================================================
--- branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/DYN/dynvor.F90	(revision 7162)
+++ branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/DYN/dynvor.F90	(revision 7163)
@@ -626,11 +626,11 @@
          WRITE(numout,*) 'dyn_vor_init : vorticity term : read namelist and control the consistency'
          WRITE(numout,*) '~~~~~~~~~~~~'
-         WRITE(numout,*) '        Namelist namdyn_vor : choice of the vorticity term scheme'
-         WRITE(numout,*) '           energy    conserving scheme                    ln_dynvor_ene = ', ln_dynvor_ene
-         WRITE(numout,*) '           enstrophy conserving scheme                    ln_dynvor_ens = ', ln_dynvor_ens
-         WRITE(numout,*) '           mixed enstrophy/energy conserving scheme       ln_dynvor_mix = ', ln_dynvor_mix
-         WRITE(numout,*) '           enstrophy and energy conserving scheme         ln_dynvor_een = ', ln_dynvor_een
-         WRITE(numout,*) '              e3f = averaging /4 (=0) or /sum(tmask) (=1)    nn_een_e3f = ', nn_een_e3f
-         WRITE(numout,*) '           masked (=T) or unmasked(=F) vorticity          ln_dynvor_msk = ', ln_dynvor_msk
+         WRITE(numout,*) '   Namelist namdyn_vor : choice of the vorticity term scheme'
+         WRITE(numout,*) '      energy    conserving scheme                    ln_dynvor_ene = ', ln_dynvor_ene
+         WRITE(numout,*) '      enstrophy conserving scheme                    ln_dynvor_ens = ', ln_dynvor_ens
+         WRITE(numout,*) '      mixed enstrophy/energy conserving scheme       ln_dynvor_mix = ', ln_dynvor_mix
+         WRITE(numout,*) '      enstrophy and energy conserving scheme         ln_dynvor_een = ', ln_dynvor_een
+         WRITE(numout,*) '         e3f = averaging /4 (=0) or /sum(tmask) (=1)    nn_een_e3f = ', nn_een_e3f
+         WRITE(numout,*) '      masked (=T) or unmasked(=F) vorticity          ln_dynvor_msk = ', ln_dynvor_msk
       ENDIF
 
@@ -639,5 +639,5 @@
       ! at angles with three ocean points and one land point
       IF(lwp) WRITE(numout,*)
-      IF(lwp) WRITE(numout,*) '           namlbc: change fmask value in the angles (T)   ln_vorlat = ', ln_vorlat
+      IF(lwp) WRITE(numout,*) '      change fmask value in the angles (T)           ln_vorlat = ', ln_vorlat
       IF( ln_vorlat .AND. ( ln_dynvor_ene .OR. ln_dynvor_ens .OR. ln_dynvor_mix ) ) THEN
          DO jk = 1, jpk
@@ -666,9 +666,9 @@
       ncor = np_COR
       IF( ln_dynadv_vec ) THEN     
-         IF(lwp) WRITE(numout,*) '         Vector form advection : vorticity = Coriolis + relative vorticity'
+         IF(lwp) WRITE(numout,*) '      ===>>   Vector form advection : vorticity = Coriolis + relative vorticity'
          nrvm = np_RVO        ! relative vorticity
          ntot = np_CRV        ! relative + planetary vorticity
       ELSE                        
-         IF(lwp) WRITE(numout,*) '         Flux form advection   : vorticity = Coriolis + metric term'
+         IF(lwp) WRITE(numout,*) '      ===>>   Flux form advection   : vorticity = Coriolis + metric term'
          nrvm = np_MET        ! metric term
          ntot = np_CME        ! Coriolis + metric term
@@ -677,8 +677,8 @@
       IF(lwp) THEN                   ! Print the choice
          WRITE(numout,*)
-         IF( nvor_scheme ==  np_ENE )   WRITE(numout,*) '         vorticity scheme ==>> energy conserving scheme'
-         IF( nvor_scheme ==  np_ENS )   WRITE(numout,*) '         vorticity scheme ==>> enstrophy conserving scheme'
-         IF( nvor_scheme ==  np_MIX )   WRITE(numout,*) '         vorticity scheme ==>> mixed enstrophy/energy conserving scheme'
-         IF( nvor_scheme ==  np_EEN )   WRITE(numout,*) '         vorticity scheme ==>> energy and enstrophy conserving scheme'
+         IF( nvor_scheme ==  np_ENE )   WRITE(numout,*) '      ===>>   energy conserving scheme'
+         IF( nvor_scheme ==  np_ENS )   WRITE(numout,*) '      ===>>   enstrophy conserving scheme'
+         IF( nvor_scheme ==  np_MIX )   WRITE(numout,*) '      ===>>   mixed enstrophy/energy conserving scheme'
+         IF( nvor_scheme ==  np_EEN )   WRITE(numout,*) '      ===>>   energy and enstrophy conserving scheme'
       ENDIF
       !
Index: branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf.F90
===================================================================
--- branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf.F90	(revision 7162)
+++ branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf.F90	(revision 7163)
@@ -119,6 +119,6 @@
          WRITE(numout,*) 'dyn_zdf_init : vertical dynamics physics scheme'
          WRITE(numout,*) '~~~~~~~~~~~'
-         IF( nzdf ==  0 )   WRITE(numout,*) '              Explicit time-splitting scheme'
-         IF( nzdf ==  1 )   WRITE(numout,*) '              Implicit (euler backward) scheme'
+         IF( nzdf ==  0 )   WRITE(numout,*) '      ===>>   Explicit time-splitting scheme'
+         IF( nzdf ==  1 )   WRITE(numout,*) '      ===>>   Implicit (euler backward) scheme'
       ENDIF
       !
Index: branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/DYN/wet_dry.F90
===================================================================
--- branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/DYN/wet_dry.F90	(revision 7162)
+++ branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/DYN/wet_dry.F90	(revision 7163)
@@ -1,3 +1,2 @@
-
 MODULE wet_dry
    !!==============================================================================
@@ -77,5 +76,5 @@
          WRITE(numout,*)
          WRITE(numout,*) 'wad_init : Wetting and drying initialization through namelist read'
-         WRITE(numout,*) '~~~~~~~ '
+         WRITE(numout,*) '~~~~~~~~'
          WRITE(numout,*) '   Namelist namwad'
          WRITE(numout,*) '      Logical activation                 ln_wd      = ', ln_wd
@@ -116,5 +115,4 @@
       REAL(wp), POINTER,  DIMENSION(:,:) ::   zflxu,  zflxv            ! local 2D workspace
       REAL(wp), POINTER,  DIMENSION(:,:) ::   zflxu1, zflxv1           ! local 2D workspace
-
       !!----------------------------------------------------------------------
       !
@@ -246,4 +244,5 @@
    END SUBROUTINE wad_lmt
 
+
    SUBROUTINE wad_lmt_bt( zflxu, zflxv, sshn_e, zssh_frc, rdtbt )
       !!----------------------------------------------------------------------
@@ -269,5 +268,4 @@
       REAL(wp), POINTER,  DIMENSION(:,:) ::   zflxp,  zflxn            ! local 2D workspace
       REAL(wp), POINTER,  DIMENSION(:,:) ::   zflxu1, zflxv1           ! local 2D workspace
-
       !!----------------------------------------------------------------------
       !
@@ -389,4 +387,7 @@
 
       IF( nn_timing == 1 )  CALL timing_stop('wad_lmt')
+      !
    END SUBROUTINE wad_lmt_bt
+
+   !!==============================================================================
 END MODULE wet_dry
Index: branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/ICB/icbini.F90
===================================================================
--- branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/ICB/icbini.F90	(revision 7162)
+++ branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/ICB/icbini.F90	(revision 7163)
@@ -372,7 +372,7 @@
       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,*) ' namelist namberg not read'
+         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.      
@@ -381,6 +381,6 @@
          IF(lwp) THEN
             WRITE(numout,*)
-            WRITE(numout,*) 'icbini :   Namelist namberg ln_icebergs = F , NO icebergs used'
-            WRITE(numout,*) '~~~~~~~~ '
+            WRITE(numout,*) 'icbini : Namelist namberg ln_icebergs = F , NO icebergs used'
+            WRITE(numout,*) '~~~~~~ '
          ENDIF
          RETURN
Index: branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/IOM/in_out_manager.F90
===================================================================
--- branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/IOM/in_out_manager.F90	(revision 7162)
+++ branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/IOM/in_out_manager.F90	(revision 7163)
@@ -18,6 +18,4 @@
    PUBLIC
 
- 
-   !
    !!----------------------------------------------------------------------
    !!                   namrun namelist parameters
@@ -95,17 +93,16 @@
    !!                    output monitoring
    !!----------------------------------------------------------------------
-   LOGICAL ::   ln_ctl       !: run control for debugging
-   INTEGER ::   nn_timing    !: run control for timing
-   INTEGER ::   nn_diacfl    !: flag whether to create CFL diagnostics
-   INTEGER ::   nn_print     !: level of print (0 no print)
-   INTEGER ::   nn_ictls     !: Start i indice for the SUM control
-   INTEGER ::   nn_ictle     !: End   i indice for the SUM control
-   INTEGER ::   nn_jctls     !: Start j indice for the SUM control
-   INTEGER ::   nn_jctle     !: End   j indice for the SUM control
-   INTEGER ::   nn_isplt     !: number of processors following i
-   INTEGER ::   nn_jsplt     !: number of processors following j
-   INTEGER ::   nn_bench     !: benchmark parameter (0/1)
-   INTEGER ::   nn_bit_cmp   =    0    !: bit reproducibility  (0/1)
-
+   LOGICAL ::   ln_ctl           !: run control for debugging
+   INTEGER ::   nn_timing        !: run control for timing
+   INTEGER ::   nn_diacfl        !: flag whether to create CFL diagnostics
+   INTEGER ::   nn_print         !: level of print (0 no print)
+   INTEGER ::   nn_ictls         !: Start i indice for the SUM control
+   INTEGER ::   nn_ictle         !: End   i indice for the SUM control
+   INTEGER ::   nn_jctls         !: Start j indice for the SUM control
+   INTEGER ::   nn_jctle         !: End   j indice for the SUM control
+   INTEGER ::   nn_isplt         !: number of processors following i
+   INTEGER ::   nn_jsplt         !: number of processors following j
+   INTEGER ::   nn_bench         !: benchmark parameter (0/1)
+   INTEGER ::   nn_bit_cmp = 0   !: bit reproducibility  (0/1)
    !                                          
    INTEGER ::   nprint, nictls, nictle, njctls, njctle, isplt, jsplt, nbench    !: OLD namelist names
@@ -138,4 +135,5 @@
    !!                          Run control  
    !!----------------------------------------------------------------------
+   INTEGER       ::   no_print = 0          !: optional argument of fld_fill (if present, suppress some control print)
    INTEGER       ::   nstop = 0             !: error flag (=number of reason for a premature stop run)
    INTEGER       ::   nwarn = 0             !: warning flag (=number of warning found during the run)
Index: branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/LBC/lib_mpp.F90
===================================================================
--- branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/LBC/lib_mpp.F90	(revision 7162)
+++ branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/LBC/lib_mpp.F90	(revision 7163)
@@ -4303,4 +4303,5 @@
             WRITE(kout,*)
          ENDIF
+         CALL FLUSH(kout) 
          STOP 'ctl_opn bad opening'
       ENDIF
Index: branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/LBC/mppini.F90
===================================================================
--- branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/LBC/mppini.F90	(revision 7162)
+++ branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/LBC/mppini.F90	(revision 7163)
@@ -65,5 +65,5 @@
          WRITE(numout,*)
          WRITE(numout,*) 'mpp_init(2) : NO massively parallel processing'
-         WRITE(numout,*) '~~~~~~~~~~~: '
+         WRITE(numout,*) '~~~~~~~~~~~ '
          WRITE(numout,*) '         nperio = ', nperio
          WRITE(numout,*) '         npolj  = ', npolj
@@ -265,10 +265,7 @@
       
       IF(lwp) WRITE(numout,*)
-      IF(lwp) WRITE(numout,*) ' mpp_init: defines mpp subdomains'
-      IF(lwp) WRITE(numout,*) ' ~~~~~~  ----------------------'
-      IF(lwp) WRITE(numout,*)
-      IF(lwp) WRITE(numout,*) 'iresti=',iresti,' irestj=',irestj
-      IF(lwp) WRITE(numout,*)
-      IF(lwp) WRITE(numout,*) 'jpni=',jpni,' jpnj=',jpnj
+      IF(lwp) WRITE(numout,*) '   defines mpp subdomains'
+      IF(lwp) WRITE(numout,*) '      jpni=', jpni, ' iresti=', iresti
+      IF(lwp) WRITE(numout,*) '      jpnj=', jpnj, ' irestj=', irestj
       zidom = nreci
       DO ji = 1, jpni
@@ -276,5 +273,5 @@
       END DO
       IF(lwp) WRITE(numout,*)
-      IF(lwp) WRITE(numout,*)' sum ilcit(i,1)=', zidom, ' jpiglo=', jpiglo
+      IF(lwp) WRITE(numout,*)'      sum ilcit(i,1)=', zidom, ' jpiglo=', jpiglo
 
       zjdom = nrecj
@@ -282,6 +279,5 @@
          zjdom = zjdom + ilcjt(1,jj) - nrecj
       END DO
-      IF(lwp) WRITE(numout,*)' sum ilcit(1,j)=', zjdom, ' jpjglo=', jpjglo
-      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*)'      sum ilcit(1,j)=', zjdom, ' jpjglo=', jpjglo
 
       IF(lwp) THEN
@@ -360,20 +356,19 @@
       njmpp  = njmppt(narea)  
 
-     ! Save processor layout in layout.dat file 
-       IF (lwp) THEN
-        CALL ctl_opn( inum, 'layout.dat', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, .FALSE., narea )
-        WRITE(inum,'(a)') '   jpnij     jpi     jpj     jpk  jpiglo  jpjglo'
-        WRITE(inum,'(6i8)') jpnij,jpi,jpj,jpk,jpiglo,jpjglo
-        WRITE(inum,'(a)') 'NAREA nlci nlcj nldi nldj nlei nlej nimpp njmpp'
-
-        DO  jn = 1, jpnij
-         WRITE(inum,'(9i5)') jn, nlcit(jn), nlcjt(jn), &
-                                      nldit(jn), nldjt(jn), &
-                                      nleit(jn), nlejt(jn), &
-                                      nimppt(jn), njmppt(jn)
-        END DO
-        CLOSE(inum)   
+      ! Save processor layout in layout.dat file 
+      IF(lwp) THEN
+         CALL ctl_opn( inum, 'layout.dat', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, .FALSE., narea )
+         WRITE(inum,'(a)') '   jpnij     jpi     jpj     jpk  jpiglo  jpjglo'
+         WRITE(inum,'(6i8)') jpnij,jpi,jpj,jpk,jpiglo,jpjglo
+         WRITE(inum,'(a)') 'NAREA nlci nlcj nldi nldj nlei nlej nimpp njmpp'
+         !
+         DO jn = 1, jpnij
+            WRITE(inum,'(9i5)') jn, nlcit(jn), nlcjt(jn), &
+               &                    nldit(jn), nldjt(jn), &
+               &                    nleit(jn), nlejt(jn), &
+               &                    nimppt(jn), njmppt(jn)
+         END DO
+         CLOSE(inum)   
       END IF
-
 
       ! w a r n i n g  narea (zone) /= nproc (processors)!
@@ -428,19 +423,16 @@
 
       IF(lwp) THEN
-         WRITE(numout,*) ' nproc  = ', nproc
-         WRITE(numout,*) ' nowe   = ', nowe  , ' noea   =  ', noea
-         WRITE(numout,*) ' nono   = ', nono  , ' noso   =  ', noso
-         WRITE(numout,*) ' nbondi = ', nbondi
-         WRITE(numout,*) ' nbondj = ', nbondj
-         WRITE(numout,*) ' npolj  = ', npolj
-         WRITE(numout,*) ' nperio = ', nperio
-         WRITE(numout,*) ' nlci   = ', nlci
-         WRITE(numout,*) ' nlcj   = ', nlcj
-         WRITE(numout,*) ' nimpp  = ', nimpp
-         WRITE(numout,*) ' njmpp  = ', njmpp
-         WRITE(numout,*) ' nreci  = ', nreci  , ' npse   = ', npse
-         WRITE(numout,*) ' nrecj  = ', nrecj  , ' npsw   = ', npsw
-         WRITE(numout,*) ' jpreci = ', jpreci , ' npne   = ', npne
-         WRITE(numout,*) ' jprecj = ', jprecj , ' npnw   = ', npnw
+         WRITE(numout,*) '      nproc  = ', nproc
+         WRITE(numout,*) '      nowe   = ', nowe  , '      noea   =  ', noea
+         WRITE(numout,*) '      nono   = ', nono  , '      noso   =  ', noso
+         WRITE(numout,*) '      nbondi = ', nbondi, '      nbondj = ', nbondj
+         WRITE(numout,*) '      npolj  = ', npolj
+         WRITE(numout,*) '      nperio = ', nperio
+         WRITE(numout,*) '      nlci   = ', nlci  , '      nlcj   = ', nlcj
+         WRITE(numout,*) '      nimpp  = ', nimpp , '      njmpp  = ', njmpp
+         WRITE(numout,*) '      nreci  = ', nreci , '      npse   = ', npse
+         WRITE(numout,*) '      nrecj  = ', nrecj , '      npsw   = ', npsw
+         WRITE(numout,*) '      jpreci = ', jpreci, '      npne   = ', npne
+         WRITE(numout,*) '      jprecj = ', jprecj, '      npnw   = ', npnw
          WRITE(numout,*)
       ENDIF
Index: branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/LDF/ldfc1d_c2d.F90
===================================================================
--- branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/LDF/ldfc1d_c2d.F90	(revision 7162)
+++ branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/LDF/ldfc1d_c2d.F90	(revision 7163)
@@ -59,4 +59,10 @@
       REAL(wp) ::   zw    , zdep2   !   -      -
       !!----------------------------------------------------------------------
+
+      IF(lwp) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) '   ldf_c1d : set a given profile to eddy diffusivity/viscosity coefficients'
+         WRITE(numout,*) '   ~~~~~~~'
+      ENDIF
 
       ! initialization of the profile
@@ -130,6 +136,7 @@
       !
       IF(lwp) THEN
-         WRITE(numout,*) 'ldf_c2d :     aht = rn_aht0 *  max(e1,e2)/e_equator     (  laplacian) '
-         WRITE(numout,*) '~~~~~~~       or  = rn_bht0 * [max(e1,e2)/e_equator]**3 (bilaplacian)'
+         WRITE(numout,*)
+         WRITE(numout,*) '   ldf_c2d :     aht = rn_aht0 *  max(e1,e2)/e_equator     (  laplacian) '
+         WRITE(numout,*) '   ~~~~~~~       or  = rn_bht0 * [max(e1,e2)/e_equator]**3 (bilaplacian)'
          WRITE(numout,*)
       ENDIF
Index: branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/SBC/fldread.F90
===================================================================
--- branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/SBC/fldread.F90	(revision 7162)
+++ branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/SBC/fldread.F90	(revision 7163)
@@ -4,12 +4,25 @@
    !! Ocean forcing:  read input field for surface boundary condition
    !!=====================================================================
-   !! History :  2.0  !  06-2006  (S. Masson, G. Madec)  Original code
-   !!                 !  05-2008  (S. Alderson)  Modified for Interpolation in memory from input grid to model grid
-   !!                 !  10-2013  (D. Delrosso, P. Oddo)  suppression of land point prior to interpolation
+   !! History :  2.0  !  2006-06  (S. Masson, G. Madec)  Original code
+   !!            3.0  !  2008-05  (S. Alderson)  Modified for Interpolation in memory from input grid to model grid
+   !!            3.4  !  2013-10  (D. Delrosso, P. Oddo)  suppression of land point prior to interpolation
    !!----------------------------------------------------------------------
 
    !!----------------------------------------------------------------------
-   !!   fld_read      : read input fields used for the computation of the
-   !!                   surface boundary condition
+   !!   fld_read      : read input fields used for the computation of the surface boundary condition
+   !!   fld_init      : initialization of field read
+   !!   fld_rec       : determined the record(s) to be read
+   !!   fld_get       : read the data
+   !!   fld_map       : read global data from file and map onto local data using a general mapping (use for open boundaries)
+   !!   fld_rot       : rotate the vector fields onto the local grid direction
+   !!   fld_clopn     : update the data file name and close/open the files
+   !!   fld_fill      : fill the data structure with the associated information read in namelist
+   !!   wgt_list      : manage the weights used for interpolation
+   !!   wgt_print     : print the list of known weights
+   !!   fld_weight    : create a WGT structure and fill in data from file, restructuring as required
+   !!   apply_seaoverland : fill land with ocean values
+   !!   seaoverland   : create shifted matrices for seaoverland application
+   !!   fld_interp    : apply weights to input gridded data to create data on model grid
+   !!   ksec_week     : function returning the first 3 letters of the first day of the weekly file
    !!----------------------------------------------------------------------
    USE oce            ! ocean dynamics and tracers
@@ -274,9 +287,9 @@
             IF( sd(jf)%ln_tint ) THEN              ! temporal interpolation
                IF(lwp .AND. kt - nit000 <= 100 ) THEN 
-                  clfmt = "('fld_read: var ', a, ' kt = ', i8, ' (', f9.4,' days), Y/M/D = ', i4.4,'/', i2.2,'/', i2.2," //   &
+                  clfmt = "('   fld_read: var ', a, ' kt = ', i8, ' (', f9.4,' days), Y/M/D = ', i4.4,'/', i2.2,'/', i2.2," //   &
                      &    "', records b/a: ', i6.4, '/', i6.4, ' (days ', f9.4,'/', f9.4, ')')"
                   WRITE(numout, clfmt)  TRIM( sd(jf)%clvar ), kt, REAL(isecsbc,wp)/rday, nyear, nmonth, nday,   &            
                      & sd(jf)%nrec_b(1), sd(jf)%nrec_a(1), REAL(sd(jf)%nrec_b(2),wp)/rday, REAL(sd(jf)%nrec_a(2),wp)/rday
-                  WRITE(numout, *) 'it_offset is : ',it_offset
+                  WRITE(numout, *) '      it_offset is : ',it_offset
                ENDIF
                ! temporal interpolation weights
@@ -286,5 +299,5 @@
             ELSE   ! nothing to do...
                IF(lwp .AND. kt - nit000 <= 100 ) THEN
-                  clfmt = "('fld_read: var ', a, ' kt = ', i8,' (', f9.4,' days), Y/M/D = ', i4.4,'/', i2.2,'/', i2.2," //   &
+                  clfmt = "('   fld_read: var ', a, ' kt = ', i8,' (', f9.4,' days), Y/M/D = ', i4.4,'/', i2.2,'/', i2.2," //   &
                      &    "', record: ', i6.4, ' (days ', f9.4, ' <-> ', f9.4, ')')"
                   WRITE(numout, clfmt) TRIM(sd(jf)%clvar), kt, REAL(isecsbc,wp)/rday, nyear, nmonth, nday,    &
@@ -407,5 +420,5 @@
          CALL fld_get( sdjf, map )         ! read before data in after arrays(as we will swap it later)
          !
-         clfmt = "('fld_init : time-interpolation for ', a, ' read previous record = ', i6, ' at time = ', f7.2, ' days')"
+         clfmt = "('   fld_init : time-interpolation for ', a, ' read previous record = ', i6, ' at time = ', f7.2, ' days')"
          IF(lwp) WRITE(numout, clfmt) TRIM(sdjf%clvar), sdjf%nrec_a(1), REAL(sdjf%nrec_a(2),wp)/rday
          !
@@ -791,5 +804,5 @@
       !!                    ***  ROUTINE fld_clopn  ***
       !!
-      !! ** Purpose :   update the file name and open the file
+      !! ** Purpose :   update the file name and close/open the files
       !!----------------------------------------------------------------------
       TYPE(FLD)        , INTENT(inout) ::   sdjf     ! input field related variables
@@ -882,18 +895,20 @@
 
 
-   SUBROUTINE fld_fill( sdf, sdf_n, cdir, cdcaller, cdtitle, cdnam )
+   SUBROUTINE fld_fill( sdf, sdf_n, cdir, cdcaller, cdtitle, cdnam, knoprint )
       !!---------------------------------------------------------------------
       !!                    ***  ROUTINE fld_fill  ***
       !!
-      !! ** Purpose :   fill sdf with sdf_n and control print
-      !!----------------------------------------------------------------------
-      TYPE(FLD)  , DIMENSION(:), INTENT(inout) ::   sdf        ! structure of input fields (file informations, fields read)
-      TYPE(FLD_N), DIMENSION(:), INTENT(in   ) ::   sdf_n      ! array of namelist information structures
-      CHARACTER(len=*)         , INTENT(in   ) ::   cdir       ! Root directory for location of flx files
-      CHARACTER(len=*)         , INTENT(in   ) ::   cdcaller   ! 
-      CHARACTER(len=*)         , INTENT(in   ) ::   cdtitle    ! 
-      CHARACTER(len=*)         , INTENT(in   ) ::   cdnam      ! 
-      !
-      INTEGER  ::   jf       ! dummy indices
+      !! ** Purpose :   fill the data structure (sdf) with the associated information 
+      !!              read in namelist (sdf_n) and control print
+      !!----------------------------------------------------------------------
+      TYPE(FLD)  , DIMENSION(:)          , INTENT(inout) ::   sdf        ! structure of input fields (file informations, fields read)
+      TYPE(FLD_N), DIMENSION(:)          , INTENT(in   ) ::   sdf_n      ! array of namelist information structures
+      CHARACTER(len=*)                   , INTENT(in   ) ::   cdir       ! Root directory for location of flx files
+      CHARACTER(len=*)                   , INTENT(in   ) ::   cdcaller   ! name of the calling routine
+      CHARACTER(len=*)                   , INTENT(in   ) ::   cdtitle    ! description of the calling routine 
+      CHARACTER(len=*)                   , INTENT(in   ) ::   cdnam      ! name of the namelist from which sdf_n comes
+      INTEGER                  , OPTIONAL, INTENT(in   ) ::   knoprint   ! no calling routine information printed
+      !
+      INTEGER  ::   jf   ! dummy indices
       !!---------------------------------------------------------------------
       !
@@ -922,10 +937,13 @@
       IF(lwp) THEN      ! control print
          WRITE(numout,*)
-         WRITE(numout,*) TRIM( cdcaller )//' : '//TRIM( cdtitle )
-         WRITE(numout,*) (/ ('~', jf = 1, LEN_TRIM( cdcaller ) ) /)
-         WRITE(numout,*) '   Namelist '//TRIM( cdnam )
+         IF( .NOT.PRESENT( knoprint) ) THEN
+            WRITE(numout,*) TRIM( cdcaller )//' : '//TRIM( cdtitle )
+            WRITE(numout,*) (/ ('~', jf = 1, LEN_TRIM( cdcaller ) ) /)
+         ENDIF
+         WRITE(numout,*) '   fld_fill : fill data structure with information from namelist '//TRIM( cdnam )
+         WRITE(numout,*) '   ~~~~~~~~'
          WRITE(numout,*) '      list of files and frequency (>0: in hours ; <0 in months)'
          DO jf = 1, SIZE(sdf)
-            WRITE(numout,*) '      root filename: '  , TRIM( sdf(jf)%clrootname ), '   variable name: ', TRIM( sdf(jf)%clvar)
+            WRITE(numout,*) '      root filename: '  , TRIM( sdf(jf)%clrootname ), '   variable name: ', TRIM( sdf(jf)%clvar )
             WRITE(numout,*) '         frequency: '      ,       sdf(jf)%nfreqh      ,   &
                &                  '   time interp: '    ,       sdf(jf)%ln_tint     ,   &
@@ -946,8 +964,7 @@
       !!                    ***  ROUTINE wgt_list  ***
       !!
-      !! ** Purpose :   search array of WGTs and find a weights file
-      !!                entry, or return a new one adding it to the end
-      !!                if it is a new entry, the weights data is read in and
-      !!                restructured (fld_weight)
+      !! ** Purpose :   search array of WGTs and find a weights file entry,
+      !!                or return a new one adding it to the end if new entry.
+      !!                the weights data is read in and restructured (fld_weight)
       !!----------------------------------------------------------------------
       TYPE( FLD ), INTENT(in   ) ::   sd        ! field with name of weights file
@@ -1018,6 +1035,6 @@
       !!                    ***  ROUTINE fld_weight  ***
       !!
-      !! ** Purpose :   create a new WGT structure and fill in data from  
-      !!                file, restructuring as required
+      !! ** Purpose :   create a new WGT structure and fill in data from file,
+      !!              restructuring as required
       !!----------------------------------------------------------------------
       TYPE( FLD ), INTENT(in) ::   sd   ! field with name of weights file
@@ -1162,5 +1179,5 @@
 
    SUBROUTINE apply_seaoverland( clmaskfile, zfieldo, jpi1_lsm, jpi2_lsm, jpj1_lsm,   &
-                          &      jpj2_lsm, itmpi, itmpj, itmpz, rec1_lsm, recn_lsm )
+      &                          jpj2_lsm, itmpi, itmpj, itmpz, rec1_lsm, recn_lsm )
       !!---------------------------------------------------------------------
       !!                    ***  ROUTINE apply_seaoverland  ***
@@ -1491,11 +1508,10 @@
       !!                    ***  FUNCTION kshift_week *** 
       !!
-      !! ** Purpose :  
-      !!---------------------------------------------------------------------
-      CHARACTER(len=*), INTENT(in)   ::   cdday   !3 first letters of the first day of the weekly file
-      !!
-      INTEGER                        ::   ksec_week  ! output variable
-      INTEGER                        ::   ijul       !temp variable
-      INTEGER                        ::   ishift     !temp variable
+      !! ** Purpose :   return the first 3 letters of the first day of the weekly file
+      !!---------------------------------------------------------------------
+      CHARACTER(len=*), INTENT(in)   ::   cdday   ! first 3 letters of the first day of the weekly file
+      !!
+      INTEGER                        ::   ksec_week      ! output variable
+      INTEGER                        ::   ijul, ishift   ! local integer
       CHARACTER(len=3),DIMENSION(7)  ::   cl_week 
       !!----------------------------------------------------------------------
Index: branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk.F90
===================================================================
--- branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk.F90	(revision 7162)
+++ branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk.F90	(revision 7163)
@@ -6,24 +6,22 @@
    !!                        SUCCESSOR OF "sbcblk_core"
    !!=====================================================================
-   !! History :  1.0  !  2004-08  (U. Schweckendiek)  Original code
-   !!            2.0  !  2005-04  (L. Brodeau, A.M. Treguier) additions:
-   !!                           -  new bulk routine for efficiency
-   !!                           -  WINDS ARE NOW ASSUMED TO BE AT T POINTS in input files
-   !!                           -  file names and file characteristics in namelist
-   !!                           -  Implement reading of 6-hourly fields
-   !!            3.0  !  2006-06  (G. Madec) sbc rewritting
-   !!             -   !  2006-12  (L. Brodeau) Original code for turb_core
+   !! History :  1.0  !  2004-08  (U. Schweckendiek)  Original CORE code
+   !!            2.0  !  2005-04  (L. Brodeau, A.M. Treguier)  improved CORE bulk and its user interface
+   !!            3.0  !  2006-06  (G. Madec)  sbc rewritting
+   !!             -   !  2006-12  (L. Brodeau)  Original code for turb_core
    !!            3.2  !  2009-04  (B. Lemaire)  Introduce iom_put
    !!            3.3  !  2010-10  (S. Masson)  add diurnal cycle
-   !!            3.4  !  2011-11  (C. Harris) Fill arrays required by CICE
-   !!            3.7  !  2014-06  (L. Brodeau) simplification and optimization of CORE bulk
-   !!            4.0  !  2016-06  (L. Brodeau) sbcblk_core becomes sbcblk and is not restricted to the CORE algorithm anymore
+   !!            3.4  !  2011-11  (C. Harris)  Fill arrays required by CICE
+   !!            3.7  !  2014-06  (L. Brodeau)  simplification and optimization of CORE bulk
+   !!            4.0  !  2016-06  (L. Brodeau)  sbcblk_core becomes sbcblk and is not restricted to the CORE algorithm anymore
    !!                                          ==> based on AeroBulk (http://aerobulk.sourceforge.net/)
+   !!            4.0  !  2016-10  (G. Madec)  introduce a sbc_blk_init routine
    !!----------------------------------------------------------------------
 
    !!----------------------------------------------------------------------
-   !!   sbc_blk       : bulk formulation as ocean surface boundary condition (forced mode, CORE bulk formulea)
+   !!   sbc_blk_init  : initialisation of the chosen bulk formulation as ocean surface boundary condition
+   !!   sbc_blk       : bulk formulation as ocean surface boundary condition
    !!   blk_oce       : computes momentum, heat and freshwater fluxes over ocean
-   !!   blk_ice       : computes momentum, heat and freshwater fluxes over ice
+   !!   blk_ice       : computes momentum, heat and freshwater fluxes over sea ice
    !!   rho_air       : density of (moist) air (depends on T_air, q_air and SLP
    !!   cp_air        : specific heat of (moist) air (depends spec. hum. q_air)
@@ -64,5 +62,6 @@
    PRIVATE
 
-   PUBLIC   sbc_blk       ! routine called in sbcmod module
+   PUBLIC   sbc_blk_init  ! called in sbcmod
+   PUBLIC   sbc_blk       ! called in sbcmod
 #if defined key_lim2 || defined key_lim3
    PUBLIC   blk_ice_tau   ! routine called in sbc_ice_lim module
@@ -82,8 +81,8 @@
    INTEGER , PARAMETER ::   jp_wndi = 1           ! index of 10m wind velocity (i-component) (m/s)    at T-point
    INTEGER , PARAMETER ::   jp_wndj = 2           ! index of 10m wind velocity (j-component) (m/s)    at T-point
-   INTEGER , PARAMETER ::   jp_humi = 3           ! index of specific humidity               ( % )
-   INTEGER , PARAMETER ::   jp_qsr  = 4           ! index of solar heat                      (W/m2)
-   INTEGER , PARAMETER ::   jp_qlw  = 5           ! index of Long wave                       (W/m2)
-   INTEGER , PARAMETER ::   jp_tair = 6           ! index of 10m air temperature             (Kelvin)
+   INTEGER , PARAMETER ::   jp_tair = 3           ! index of 10m air temperature             (Kelvin)
+   INTEGER , PARAMETER ::   jp_humi = 4           ! index of specific humidity               ( % )
+   INTEGER , PARAMETER ::   jp_qsr  = 5           ! index of solar heat                      (W/m2)
+   INTEGER , PARAMETER ::   jp_qlw  = 6           ! index of Long wave                       (W/m2)
    INTEGER , PARAMETER ::   jp_prec = 7           ! index of total precipitation (rain+snow) (Kg/m2/s)
    INTEGER , PARAMETER ::   jp_snow = 8           ! index of snow (solid prcipitation)       (kg/m2/s)
@@ -129,4 +128,112 @@
 CONTAINS
 
+   SUBROUTINE sbc_blk_init
+      !!---------------------------------------------------------------------
+      !!                    ***  ROUTINE sbc_blk_init  ***
+      !!
+      !! ** Purpose :   choose and initialize a bulk formulae formulation
+      !!
+      !! ** Method  : 
+      !!
+      !!      C A U T I O N : never mask the surface stress fields
+      !!                      the stress is assumed to be in the (i,j) mesh referential
+      !!
+      !! ** Action  :   
+      !!
+      !!----------------------------------------------------------------------
+      INTEGER  ::   ifpr, jfld            ! dummy loop indice and argument
+      INTEGER  ::   ios, ierror, ioptio   ! Local integer
+      !!
+      CHARACTER(len=100)            ::   cn_dir                ! Root directory for location of atmospheric forcing files
+      TYPE(FLD_N), DIMENSION(jpfld) ::   slf_i                 ! array of namelist informations on the fields to read
+      TYPE(FLD_N) ::   sn_wndi, sn_wndj, sn_humi, sn_qsr       ! informations about the fields to be read
+      TYPE(FLD_N) ::   sn_qlw , sn_tair, sn_prec, sn_snow      !       "                        "
+      TYPE(FLD_N) ::   sn_slp , sn_tdif                        !       "                        "
+      NAMELIST/namsbc_blk/ sn_wndi, sn_wndj, sn_humi, sn_qsr, sn_qlw ,                &   ! input fields
+         &                 sn_tair, sn_prec, sn_snow, sn_slp, sn_tdif,                &
+         &                 ln_NCAR, ln_COARE_3p0, ln_COARE_3p5, ln_ECMWF,             &   ! bulk algorithm
+         &                 cn_dir , ln_taudif, rn_zqt, rn_zu, rn_pfac, rn_efac, rn_vfac
+      !!---------------------------------------------------------------------
+      !
+      !                             !** read bulk namelist  
+      REWIND( numnam_ref )                !* Namelist namsbc_blk in reference namelist : bulk parameters
+      READ  ( numnam_ref, namsbc_blk, IOSTAT = ios, ERR = 901)
+901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_blk in reference namelist', lwp )
+      !
+      REWIND( numnam_cfg )                !* Namelist namsbc_blk in configuration namelist : bulk parameters
+      READ  ( numnam_cfg, namsbc_blk, IOSTAT = ios, ERR = 902 )
+902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_blk in configuration namelist', lwp )
+      !
+      IF(lwm) WRITE( numond, namsbc_blk )
+      !
+      !                             !** initialization of the chosen bulk formulae (+ check)
+      !                                   !* select the bulk chosen in the namelist and check the choice
+      ;                                                        ioptio = 0
+      IF( ln_NCAR      ) THEN   ;   nblk =  np_NCAR        ;   ioptio = ioptio + 1   ;   ENDIF
+      IF( ln_COARE_3p0 ) THEN   ;   nblk =  np_COARE_3p0   ;   ioptio = ioptio + 1   ;   ENDIF
+      IF( ln_COARE_3p5 ) THEN   ;   nblk =  np_COARE_3p5   ;   ioptio = ioptio + 1   ;   ENDIF
+      IF( ln_ECMWF     ) THEN   ;   nblk =  np_ECMWF       ;   ioptio = ioptio + 1   ;   ENDIF
+      !
+      IF( ioptio /= 1 )   CALL ctl_stop( 'sbc_blk_init: Choose one and only one bulk algorithm' )
+      !
+      IF( ln_dm2dc ) THEN                 !* check: diurnal cycle on Qsr
+         IF( sn_qsr%nfreqh /= 24 )   CALL ctl_stop( 'sbc_blk_init: ln_dm2dc=T only with daily short-wave input' )
+         IF( sn_qsr%ln_tint ) THEN 
+            CALL ctl_warn( 'sbc_blk_init: ln_dm2dc=T daily qsr time interpolation done by sbcdcy module',   &
+               &           '              ==> We force time interpolation = .false. for qsr' )
+            sn_qsr%ln_tint = .false.
+         ENDIF
+      ENDIF
+      !                                   !* set the bulk structure
+      !                                      !- store namelist information in an array
+      slf_i(jp_wndi) = sn_wndi   ;   slf_i(jp_wndj) = sn_wndj
+      slf_i(jp_qsr ) = sn_qsr    ;   slf_i(jp_qlw ) = sn_qlw
+      slf_i(jp_tair) = sn_tair   ;   slf_i(jp_humi) = sn_humi
+      slf_i(jp_prec) = sn_prec   ;   slf_i(jp_snow) = sn_snow
+      slf_i(jp_slp)  = sn_slp    ;   slf_i(jp_tdif) = sn_tdif
+      !
+      lhftau = ln_taudif                     !- add an extra field if HF stress is used
+      jfld = jpfld - COUNT( (/.NOT.lhftau/) )
+      !
+      !                                      !- allocate the bulk structure
+      ALLOCATE( sf(jfld), STAT=ierror )
+      IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_blk_init: unable to allocate sf structure' )
+      DO ifpr= 1, jfld
+         ALLOCATE( sf(ifpr)%fnow(jpi,jpj,1) )
+         IF( slf_i(ifpr)%ln_tint )   ALLOCATE( sf(ifpr)%fdta(jpi,jpj,1,2) )
+      END DO
+      !                                      !- fill the bulk structure with namelist informations
+      CALL fld_fill( sf, slf_i, cn_dir, 'sbc_blk_init', 'surface boundary condition -- bulk formulae', 'namsbc_blk' )
+      !
+      !           
+      IF(lwp) THEN                     !** Control print
+         !
+         WRITE(numout,*)                  !* namelist 
+         WRITE(numout,*) '   Namelist namsbc_blk (other than data information):'
+         WRITE(numout,*) '      "NCAR"      algorithm   (Large and Yeager 2008)     ln_NCAR      = ', ln_NCAR
+         WRITE(numout,*) '      "COARE 3.0" algorithm   (Fairall et al. 2003)       ln_COARE_3p0 = ', ln_COARE_3p0
+         WRITE(numout,*) '      "COARE 3.5" algorithm   (Edson et al. 2013)         ln_COARE_3p5 = ', ln_COARE_3p0
+         WRITE(numout,*) '      "ECMWF"     algorithm   (IFS cycle 31)              ln_ECMWF     = ', ln_ECMWF
+         WRITE(numout,*) '      add High freq.contribution to the stress module     ln_taudif    = ', ln_taudif
+         WRITE(numout,*) '      Air temperature and humidity reference height (m)   rn_zqt       = ', rn_zqt
+         WRITE(numout,*) '      Wind vector reference height (m)                    rn_zu        = ', rn_zu
+         WRITE(numout,*) '      factor applied on precipitation (total & snow)      rn_pfac      = ', rn_pfac
+         WRITE(numout,*) '      factor applied on evaporation                       rn_efac      = ', rn_efac
+         WRITE(numout,*) '      factor applied on ocean/ice velocity                rn_vfac      = ', rn_vfac
+         WRITE(numout,*) '         (form absolute (=0) to relative winds(=1))'
+         !
+         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)'
+         END SELECT
+         !
+      ENDIF
+      !
+   END SUBROUTINE sbc_blk_init
+
+
    SUBROUTINE sbc_blk( kt )
       !!---------------------------------------------------------------------
@@ -164,100 +271,8 @@
       !!----------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt   ! ocean time step
-      !
-      INTEGER  ::   ifpr, jfld            ! dummy loop indice and argument
-      INTEGER  ::   ios, ierror, ioptio   ! Local integer
-      !
-      CHARACTER(len=100)            ::   cn_dir                ! Root directory for location of atmospheric forcing files
-      TYPE(FLD_N), DIMENSION(jpfld) ::   slf_i                 ! array of namelist informations on the fields to read
-      TYPE(FLD_N) ::   sn_wndi, sn_wndj, sn_humi, sn_qsr       ! informations about the fields to be read
-      TYPE(FLD_N) ::   sn_qlw , sn_tair, sn_prec, sn_snow      !   "                                 "
-      TYPE(FLD_N) ::   sn_slp , sn_tdif                        !   "                                 "
-      NAMELIST/namsbc_blk/ sn_wndi, sn_wndj, sn_humi, sn_qsr, sn_qlw ,                &   ! input fields
-         &                 sn_tair, sn_prec, sn_snow, sn_slp, sn_tdif,                &
-         &                 ln_NCAR, ln_COARE_3p0, ln_COARE_3p5, ln_ECMWF,             &   ! bulk algorithm
-         &                 cn_dir , ln_taudif, rn_zqt, rn_zu, rn_pfac, rn_efac, rn_vfac
-      !!---------------------------------------------------------------------
-      !
-      !                                         ! ====================== !
-      IF( kt == nit000 ) THEN                   !  First call kt=nit000  !
-         !                                      ! ====================== !
-         !
-         REWIND( numnam_ref )              ! Namelist namsbc_blk in reference namelist : bulk parameters
-         READ  ( numnam_ref, namsbc_blk, IOSTAT = ios, ERR = 901)
-901      IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_blk in reference namelist', lwp )
-         !
-         REWIND( numnam_cfg )              ! Namelist namsbc_blk in configuration namelist : bulk parameters
-         READ  ( numnam_cfg, namsbc_blk, IOSTAT = ios, ERR = 902 )
-902      IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_blk in configuration namelist', lwp )
-         !
-         IF(lwm) WRITE( numond, namsbc_blk )
-         !
-         !                          ! Control of surface pressure gradient scheme options
-         ;                                                        ioptio = 0
-         IF( ln_NCAR      ) THEN   ;   nblk =  np_NCAR        ;   ioptio = ioptio + 1   ;   ENDIF
-         IF( ln_COARE_3p0 ) THEN   ;   nblk =  np_COARE_3p0   ;   ioptio = ioptio + 1   ;   ENDIF
-         IF( ln_COARE_3p5 ) THEN   ;   nblk =  np_COARE_3p5   ;   ioptio = ioptio + 1   ;   ENDIF
-         IF( ln_ECMWF     ) THEN   ;   nblk =  np_ECMWF       ;   ioptio = ioptio + 1   ;   ENDIF
-         !
-         IF( ioptio /= 1 )   CALL ctl_stop( 'sbc_blk: Choose one and only one bulk algorithm' )
-         !
-         !                                         ! check: do we plan to use ln_dm2dc with non-daily forcing?
-         IF( ln_dm2dc .AND. sn_qsr%nfreqh /= 24 )   &
-            &   CALL ctl_stop( 'sbc_blk: ln_dm2dc can be activated only with daily short-wave forcing' )
-         IF( ln_dm2dc .AND. sn_qsr%ln_tint ) THEN
-            CALL ctl_warn( 'sbc_blk: ln_dm2dc is taking care of the temporal interpolation of daily qsr',   &
-               &         '              ==> We force time interpolation = .false. for qsr' )
-            sn_qsr%ln_tint = .false.
-         ENDIF
-         !                                         ! store namelist information in an array
-         slf_i(jp_wndi) = sn_wndi   ;   slf_i(jp_wndj) = sn_wndj
-         slf_i(jp_qsr ) = sn_qsr    ;   slf_i(jp_qlw ) = sn_qlw
-         slf_i(jp_tair) = sn_tair   ;   slf_i(jp_humi) = sn_humi
-         slf_i(jp_prec) = sn_prec   ;   slf_i(jp_snow) = sn_snow
-         slf_i(jp_slp)  = sn_slp    ;   slf_i(jp_tdif) = sn_tdif
-         !
-         lhftau = ln_taudif                        ! do we use HF tau information?
-         jfld = jpfld - COUNT( (/.NOT. lhftau/) )
-         !
-         ALLOCATE( sf(jfld), STAT=ierror )         ! set sf structure
-         IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_blk: unable to allocate sf structure' )
-         DO ifpr= 1, jfld
-            ALLOCATE( sf(ifpr)%fnow(jpi,jpj,1) )
-            IF( slf_i(ifpr)%ln_tint )   ALLOCATE( sf(ifpr)%fdta(jpi,jpj,1,2) )
-         END DO
-         !                                         ! fill sf with slf_i and control print
-         CALL fld_fill( sf, slf_i, cn_dir, 'sbc_blk', 'surface boundary condition -- bulk formula', 'namsbc_blk' )
-         !
-         IF(lwp) THEN                              ! Control print (other namelist variable)
-            WRITE(numout,*)
-            WRITE(numout,*) '      "NCAR"      algorithm   (Large and Yeager 2008)     ln_NCAR      = ', ln_NCAR
-            WRITE(numout,*) '      "COARE 3.0" algorithm   (Fairall et al. 2003)       ln_COARE_3p0 = ', ln_COARE_3p0
-            WRITE(numout,*) '      "COARE 3.5" algorithm   (Edson et al. 2013)         ln_COARE_3p5 = ', ln_COARE_3p0
-            WRITE(numout,*) '      "ECMWF"     algorithm   (IFS cycle 31)              ln_ECMWF     = ', ln_ECMWF
-            WRITE(numout,*) '      add High freq.contribution to the stress module     ln_taudif    = ', ln_taudif
-            WRITE(numout,*) '      Air temperature and humidity reference height (m)   rn_zqt       = ', rn_zqt
-            WRITE(numout,*) '      Wind vector reference height (m)                    rn_zu        = ', rn_zu
-            WRITE(numout,*) '      factor applied on precipitation (total & snow)      rn_pfac      = ', rn_pfac
-            WRITE(numout,*) '      factor applied on evaporation                       rn_efac      = ', rn_efac
-            WRITE(numout,*) '      factor applied on ocean/ice velocity                rn_vfac      = ', rn_vfac
-            WRITE(numout,*) '         (form absolute (=0) to relative winds(=1))'
-         ENDIF
-         !
-         sfx(:,:) = 0._wp                          ! salt flux; zero unless ice is present (computed in limsbc(_2).F90)
-         !
-         IF(lwp) THEN                           ! Print the choice of bulk algorithm
-            WRITE(numout,*)
-            SELECT CASE ( nblk )
-            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
-            WRITE(numout,*)
-         ENDIF
-      ENDIF
-
+      !!---------------------------------------------------------------------
+      !
       CALL fld_read( kt, nn_fsbc, sf )             ! input fields provided at the current time-step
-
+      !
       !                                            ! compute the surface ocean fluxes using bulk formulea
       IF( MOD( kt - 1, nn_fsbc ) == 0 )   CALL blk_oce( kt, sf, sst_m, ssu_m, ssv_m )
@@ -265,6 +280,6 @@
 #if defined key_cice
       IF( MOD( kt - 1, nn_fsbc ) == 0 )   THEN
-         qlw_ice(:,:,1)   = sf(jp_qlw)%fnow(:,:,1)
-         qsr_ice(:,:,1)   = sf(jp_qsr)%fnow(:,:,1)
+         qlw_ice(:,:,1)   = sf(jp_qlw )%fnow(:,:,1)
+         qsr_ice(:,:,1)   = sf(jp_qsr )%fnow(:,:,1)
          tatm_ice(:,:)    = sf(jp_tair)%fnow(:,:,1)
          qatm_ice(:,:)    = sf(jp_humi)%fnow(:,:,1)
@@ -398,5 +413,5 @@
          &                                               Cd, Ch, Ce, zt_zu, zq_zu, zU_zu )
       CASE DEFAULT
-         CALL ctl_stop( 'STOP', 'sbc_blk: non-existing bulk formula selected' )
+         CALL ctl_stop( 'STOP', 'sbc_oce: non-existing bulk formula selected' )
       END SELECT
 
Index: branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90
===================================================================
--- branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90	(revision 7162)
+++ branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90	(revision 7163)
@@ -260,5 +260,5 @@
       IF(lwm) CALL ctl_opn( numoni, 'output.namelist.ice', 'UNKNOWN', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, 1 )
       !
-      CALL ice_run                     ! set some ice run parameters
+      CALL lim_run_init                ! set some ice run parameters
       !
       !                                ! Allocate the ice arrays
@@ -323,7 +323,7 @@
 
 
-   SUBROUTINE ice_run
+   SUBROUTINE lim_run_init
       !!-------------------------------------------------------------------
-      !!                  ***  ROUTINE ice_run ***
+      !!                  ***  ROUTINE lim_run_init ***
       !!
       !! ** Purpose :   Definition some run parameter for ice model
@@ -350,17 +350,18 @@
       IF(lwp) THEN                        ! control print
          WRITE(numout,*)
-         WRITE(numout,*) 'ice_run : ice share parameters for dynamics/advection/thermo of sea-ice'
-         WRITE(numout,*) ' ~~~~~~'
-         WRITE(numout,*) '   number of ice  categories                               = ', jpl
-         WRITE(numout,*) '   number of ice  layers                                   = ', nlay_i
-         WRITE(numout,*) '   number of snow layers                                   = ', nlay_s
-         WRITE(numout,*) '   switch for ice dynamics (1) or not (0)      ln_limdyn   = ', ln_limdyn
-         WRITE(numout,*) '   maximum ice concentration for NH                        = ', rn_amax_n
-         WRITE(numout,*) '   maximum ice concentration for SH                        = ', rn_amax_s
-         WRITE(numout,*) '   Diagnose heat/salt budget or not          ln_limdiahsb  = ', ln_limdiahsb
-         WRITE(numout,*) '   Output   heat/salt budget or not          ln_limdiaout  = ', ln_limdiaout
-         WRITE(numout,*) '   control prints in ocean.out for (i,j)=(iiceprt,jiceprt) = ', ln_icectl
-         WRITE(numout,*) '   i-index for control prints (ln_icectl=true)             = ', iiceprt
-         WRITE(numout,*) '   j-index for control prints (ln_icectl=true)             = ', jiceprt
+         WRITE(numout,*) 'lim_run_init : ice share parameters for dynamics/advection/thermo of sea-ice'
+         WRITE(numout,*) '~~~~~~~~~~~~'
+         WRITE(numout,*) '   Namelist namicerun'
+         WRITE(numout,*) '      number of ice  categories                               = ', jpl
+         WRITE(numout,*) '      number of ice  layers                                   = ', nlay_i
+         WRITE(numout,*) '      number of snow layers                                   = ', nlay_s
+         WRITE(numout,*) '      switch for ice dynamics (1) or not (0)      ln_limdyn   = ', ln_limdyn
+         WRITE(numout,*) '      maximum ice concentration for NH                        = ', rn_amax_n
+         WRITE(numout,*) '      maximum ice concentration for SH                        = ', rn_amax_s
+         WRITE(numout,*) '      Diagnose heat/salt budget or not          ln_limdiahsb  = ', ln_limdiahsb
+         WRITE(numout,*) '      Output   heat/salt budget or not          ln_limdiaout  = ', ln_limdiaout
+         WRITE(numout,*) '      control prints in ocean.out for (i,j)=(iiceprt,jiceprt) = ', ln_icectl
+         WRITE(numout,*) '      i-index for control prints (ln_icectl=true)             = ', iiceprt
+         WRITE(numout,*) '      j-index for control prints (ln_icectl=true)             = ', jiceprt
       ENDIF
       !
@@ -377,5 +378,5 @@
 #endif
       !
-   END SUBROUTINE ice_run
+   END SUBROUTINE lim_run_init
 
 
@@ -407,8 +408,9 @@
       IF(lwp) THEN                        ! control print
          WRITE(numout,*)
-         WRITE(numout,*) 'ice_itd : ice cat distribution'
-         WRITE(numout,*) ' ~~~~~~'
-         WRITE(numout,*) '   shape of ice categories distribution                     nn_catbnd = ', nn_catbnd
-         WRITE(numout,*) '   mean ice thickness in the domain (used if nn_catbnd=2)   rn_himean = ', rn_himean
+         WRITE(numout,*) 'lim_itd_init : Initialization of ice cat distribution '
+         WRITE(numout,*) '~~~~~~~~~~~~'
+         WRITE(numout,*) '   Namelist namiceitd'
+         WRITE(numout,*) '      shape of ice categories distribution                     nn_catbnd = ', nn_catbnd
+         WRITE(numout,*) '      mean ice thickness in the domain (used if nn_catbnd=2)   rn_himean = ', rn_himean
       ENDIF
       !
@@ -416,7 +418,4 @@
       !- Thickness categories boundaries
       !----------------------------------
-      IF(lwp) WRITE(numout,*)
-      IF(lwp) WRITE(numout,*) 'lim_itd_init : Initialization of ice cat distribution '
-      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
       !
       hi_max(:) = 0._wp
@@ -450,6 +449,7 @@
       hi_max(jpl) = 99._wp          ! set to a big value to ensure that all ice is thinner than hi_max(jpl)
       !
-      IF(lwp) WRITE(numout,*) ' Thickness category boundaries '
-      IF(lwp) WRITE(numout,*) ' hi_max ', hi_max(0:jpl)
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*) '      Thickness category boundaries '
+      IF(lwp) WRITE(numout,*) '         hi_max ', hi_max(0:jpl)
       !
    END SUBROUTINE lim_itd_init
@@ -588,6 +588,6 @@
       afx_dyn(:,:) = 0._wp   ;   afx_thd(:,:) = 0._wp
       !
-      diag_heat(:,:) = 0._wp ;   diag_smvi(:,:) = 0._wp ;
-      diag_vice(:,:) = 0._wp ;   diag_vsnw(:,:) = 0._wp ;
+      diag_heat(:,:) = 0._wp ;   diag_smvi(:,:) = 0._wp
+      diag_vice(:,:) = 0._wp ;   diag_vsnw(:,:) = 0._wp
       !
    END SUBROUTINE sbc_lim_diag0
Index: branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90
===================================================================
--- branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90	(revision 7162)
+++ branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90	(revision 7163)
@@ -36,5 +36,5 @@
    USE sbcice_lim_2   ! surface boundary condition: LIM 2.0 sea-ice model
    USE sbcice_cice    ! surface boundary condition: CICE    sea-ice model
-   USE sbccpl         ! surface boundary condition: coupled florulation
+   USE sbccpl         ! surface boundary condition: coupled formulation
    USE cpl_oasis3     ! OASIS routines for coupling
    USE sbcssr         ! surface boundary condition: sea surface restoring
@@ -83,16 +83,15 @@
       !!              - nsbc: type of sbc
       !!----------------------------------------------------------------------
-      INTEGER ::   icpt   ! local integer
-      !!
-      NAMELIST/namsbc/ nn_fsbc  , ln_ana   , ln_flx, ln_blk, ln_cpl   , ln_mixcpl,        &
-         &             nn_components      , nn_limflx  ,                                  &
-         &             ln_traqsr, ln_dm2dc ,                                              &
-         &             nn_ice   , nn_ice_embd,                                            &
-         &             ln_rnf   , ln_ssr   , ln_isf   , nn_fwb    , ln_apr_dyn,           &
-         &             ln_wave  ,                                                         &
+      INTEGER ::   ios, icpt                         ! local integer
+      LOGICAL ::   ll_purecpl, ll_opa, ll_not_nemo   ! local logical
+      !!
+      NAMELIST/namsbc/ nn_fsbc  ,                                                    &
+         &             ln_ana   , ln_flx   , ln_blk       ,                          &
+         &             ln_cpl   , ln_mixcpl, nn_components, nn_limflx,               &
+         &             nn_ice   , nn_ice_embd,                                       &
+         &             ln_traqsr, ln_dm2dc ,                                         &
+         &             ln_rnf   , nn_fwb   , ln_ssr       , ln_isf   , ln_apr_dyn,   &
+         &             ln_wave  ,                                                    &
          &             nn_lsm
-      INTEGER  ::   ios
-      INTEGER  ::   ierr, ierr0, ierr1, ierr2, ierr3, jpm
-      LOGICAL  ::   ll_purecpl
       !!----------------------------------------------------------------------
       !
@@ -103,14 +102,15 @@
       ENDIF
       !
-      REWIND( numnam_ref )              ! Namelist namsbc in reference namelist : Surface boundary
+      !                       !**  read Surface Module namelist
+      REWIND( numnam_ref )          !* Namelist namsbc in reference namelist : Surface boundary
       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
+      REWIND( numnam_cfg )          !* Namelist namsbc in configuration namelist : Parameters of the run
       READ  ( numnam_cfg, namsbc, IOSTAT = ios, ERR = 902 )
 902   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namsbc in configuration namelist', lwp )
       IF(lwm) WRITE( numond, namsbc )
       !
-      !                          ! overwrite namelist parameter using CPP key information
+      !                             !* overwrite namelist parameter using CPP key information
       IF( Agrif_Root() ) THEN                ! AGRIF zoom
          IF( lk_lim2 )   nn_ice      = 2
@@ -123,101 +123,135 @@
       ENDIF
       !
-      IF(lwp) THEN               ! Control print
-         WRITE(numout,*) '        Namelist namsbc (partly overwritten with CPP key setting)'
-         WRITE(numout,*) '           frequency update of sbc (and ice)             nn_fsbc     = ', nn_fsbc
-         WRITE(numout,*) '           Type of air-sea fluxes : '
-         WRITE(numout,*) '              analytical formulation                     ln_ana        = ', ln_ana
-         WRITE(numout,*) '              flux       formulation                     ln_flx        = ', ln_flx
-         WRITE(numout,*) '              bulk       formulation                     ln_blk        = ', ln_blk
-         WRITE(numout,*) '           Type of coupling (Ocean/Ice/Atmosphere) : '
-         WRITE(numout,*) '              ocean-atmosphere coupled formulation       ln_cpl        = ', ln_cpl
-         WRITE(numout,*) '              forced-coupled mixed formulation           ln_mixcpl     = ', ln_mixcpl
-         WRITE(numout,*) '              OASIS coupling (with atm or sas)           lk_oasis      = ', lk_oasis
-         WRITE(numout,*) '              components of your executable              nn_components = ', nn_components
-         WRITE(numout,*) '              Multicategory heat flux formulation (LIM3) nn_limflx     = ', nn_limflx
-         WRITE(numout,*) '           Sea-ice : '
-         WRITE(numout,*) '              ice management in the sbc (=0/1/2/3)       nn_ice        = ', nn_ice
-         WRITE(numout,*) '              ice-ocean embedded/levitating (=0/1/2)     nn_ice_embd   = ', nn_ice_embd
-         WRITE(numout,*) '           Misc. options of sbc : '
-         WRITE(numout,*) '              Light penetration in temperature Eq.       ln_traqsr     = ', ln_traqsr
-         WRITE(numout,*) '                 daily mean to diurnal cycle qsr            ln_dm2dc   = ', ln_dm2dc
-         WRITE(numout,*) '              Sea Surface Restoring on SST and/or SSS    ln_ssr        = ', ln_ssr
-         WRITE(numout,*) '              FreshWater Budget control  (=0/1/2)        nn_fwb        = ', nn_fwb
-         WRITE(numout,*) '              Patm gradient added in ocean & ice Eqs.    ln_apr_dyn    = ', ln_apr_dyn
-         WRITE(numout,*) '              runoff / runoff mouths                     ln_rnf        = ', ln_rnf
-         WRITE(numout,*) '              iceshelf formulation                       ln_isf        = ', ln_isf
-         WRITE(numout,*) '              closed sea (=0/1) (set in namdom)          nn_closea     = ', nn_closea
-         WRITE(numout,*) '              nb of iterations if land-sea-mask applied  nn_lsm        = ', nn_lsm
-         WRITE(numout,*) '              surface wave                               ln_wave       = ', ln_wave
-      ENDIF
-      !
-      IF(lwp) THEN
-         WRITE(numout,*)
-         SELECT CASE ( nn_limflx )        ! LIM3 Multi-category heat flux formulation
-         CASE ( -1 )   ;   WRITE(numout,*) '   LIM3: use per-category fluxes (nn_limflx = -1) '
-         CASE ( 0  )   ;   WRITE(numout,*) '   LIM3: use average per-category fluxes (nn_limflx = 0) '
-         CASE ( 1  )   ;   WRITE(numout,*) '   LIM3: use average then redistribute per-category fluxes (nn_limflx = 1) '
-         CASE ( 2  )   ;   WRITE(numout,*) '   LIM3: Redistribute a single flux over categories (nn_limflx = 2) '
+      IF(lwp) THEN                  !* Control print
+         WRITE(numout,*) '   Namelist namsbc (partly overwritten with CPP key setting)'
+         WRITE(numout,*) '      frequency update of sbc (and ice)             nn_fsbc       = ', nn_fsbc
+         WRITE(numout,*) '      Type of air-sea fluxes : '
+         WRITE(numout,*) '         analytical formulation                     ln_ana        = ', ln_ana
+         WRITE(numout,*) '         flux       formulation                     ln_flx        = ', ln_flx
+         WRITE(numout,*) '         bulk       formulation                     ln_blk        = ', ln_blk
+         WRITE(numout,*) '      Type of coupling (Ocean/Ice/Atmosphere) : '
+         WRITE(numout,*) '         ocean-atmosphere coupled formulation       ln_cpl        = ', ln_cpl
+         WRITE(numout,*) '         mixed forced-coupled     formulation       ln_mixcpl     = ', ln_mixcpl
+!!gm  lk_oasis is controlled by key_oasis3  ===>>>  It shoud be removed from the namelist 
+         WRITE(numout,*) '         OASIS coupling (with atm or sas)           lk_oasis      = ', lk_oasis
+         WRITE(numout,*) '         components of your executable              nn_components = ', nn_components
+         WRITE(numout,*) '         Multicategory heat flux formulation (LIM3) nn_limflx     = ', nn_limflx
+         WRITE(numout,*) '      Sea-ice : '
+         WRITE(numout,*) '         ice management in the sbc (=0/1/2/3)       nn_ice        = ', nn_ice
+         WRITE(numout,*) '         ice-ocean embedded/levitating (=0/1/2)     nn_ice_embd   = ', nn_ice_embd
+         WRITE(numout,*) '      Misc. options of sbc : '
+         WRITE(numout,*) '         Light penetration in temperature Eq.       ln_traqsr     = ', ln_traqsr
+         WRITE(numout,*) '            daily mean to diurnal cycle qsr            ln_dm2dc   = ', ln_dm2dc
+         WRITE(numout,*) '         Sea Surface Restoring on SST and/or SSS    ln_ssr        = ', ln_ssr
+         WRITE(numout,*) '         FreshWater Budget control  (=0/1/2)        nn_fwb        = ', nn_fwb
+         WRITE(numout,*) '         Patm gradient added in ocean & ice Eqs.    ln_apr_dyn    = ', ln_apr_dyn
+         WRITE(numout,*) '         runoff / runoff mouths                     ln_rnf        = ', ln_rnf
+         WRITE(numout,*) '         iceshelf formulation                       ln_isf        = ', ln_isf
+         WRITE(numout,*) '         closed sea (=0/1) (set in namdom)          nn_closea     = ', nn_closea
+         WRITE(numout,*) '         nb of iterations if land-sea-mask applied  nn_lsm        = ', nn_lsm
+         WRITE(numout,*) '         surface wave                               ln_wave       = ', ln_wave
+      ENDIF
+      !
+      !                       !**  check option consistency
+      !
+      IF(lwp) WRITE(numout,*)       !* Single / Multi - executable (NEMO / OPA+SAS) 
+      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' 
+      CASE( jp_iam_opa  )
+         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( .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' )
+      CASE DEFAULT
+         CALL ctl_stop( 'sbc_init : unsupported value for nn_components' )
+      END SELECT
+      !                             !* coupled options
+      IF( ln_cpl ) THEN
+         IF( .NOT. lk_oasis )   CALL ctl_stop( 'sbc_init : coupled mode with an atmosphere model (ln_cpl=T)',   &
+            &                                  '           required to defined key_oasis3' )
+      ENDIF
+      IF( ln_mixcpl ) THEN
+         IF( .NOT. lk_oasis )   CALL ctl_stop( 'sbc_init : mixed forced-coupled mode (ln_mixcpl=T) ',   &
+            &                                  '           required to defined key_oasis3' )
+         IF( .NOT.ln_cpl    )   CALL ctl_stop( 'sbc_init : mixed forced-coupled mode (ln_mixcpl=T) requires ln_cpl = T' )
+         IF( nn_components /= jp_iam_nemo )    &
+            &                   CALL ctl_stop( 'sbc_init : the mixed forced-coupled mode (ln_mixcpl=T) ',   &
+            &                                   '          not yet working with sas-opa coupling via oasis' )
+      ENDIF
+      !                             !* sea-ice
+      SELECT CASE( nn_ice )
+      CASE( 0 )                        !- no ice in the domain
+      CASE( 1 )                        !- Ice-cover climatology ("Ice-if" model)  
+      CASE( 2 )                        !- LIM2 ice model
+         IF( .NOT.( ln_blk .OR. ln_cpl ) )   CALL ctl_stop( 'sbc_init : LIM2 sea-ice model requires ln_blk or ln_cpl = T' )
+      CASE( 3 )                        !- LIM3 ice model
+         IF( .NOT.( ln_blk .OR. ln_cpl ) )   CALL ctl_stop( 'sbc_init : LIM3 sea-ice model requires ln_blk or ln_cpl = T' )
+         IF( nn_ice_embd == 0            )   CALL ctl_stop( 'sbc_init : LIM3 sea-ice models require nn_ice_embd = 1 or 2' )
+      CASE( 4 )                        !- CICE ice model
+         IF( .NOT.( ln_blk .OR. ln_cpl ) )   CALL ctl_stop( 'sbc_init : CICE sea-ice model requires ln_blk or ln_cpl = T' )
+         IF( nn_ice_embd == 0            )   CALL ctl_stop( 'sbc_init : CICE sea-ice models require nn_ice_embd = 1 or 2' )
+         IF( lk_agrif                    )   CALL ctl_stop( 'sbc_init : CICE sea-ice model not currently available with AGRIF' )
+      CASE DEFAULT                     !- not supported
+      END SELECT
+      !
+      IF( nn_ice == 3 ) THEN           !- LIM3 case: multi-category flux option
+         IF(lwp) WRITE(numout,*)
+         SELECT CASE( nn_limflx )         ! LIM3 Multi-category heat flux formulation
+         CASE ( -1 )   
+            IF(lwp) WRITE(numout,*) '   LIM3: use per-category fluxes (nn_limflx = -1) '
+            IF( ln_cpl )   CALL ctl_stop( 'sbc_init : the chosen nn_limflx for LIM3 in coupled mode must be 0 or 2' )
+         CASE ( 0  )   
+            IF(lwp) WRITE(numout,*) '   LIM3: use average per-category fluxes (nn_limflx = 0) '
+         CASE ( 1  )
+            IF(lwp) WRITE(numout,*) '   LIM3: use average then redistribute per-category fluxes (nn_limflx = 1) '
+            IF( ln_cpl )   CALL ctl_stop( 'sbc_init : the chosen nn_limflx for LIM3 in coupled mode must be 0 or 2' )
+         CASE ( 2  )
+            IF(lwp) WRITE(numout,*) '   LIM3: Redistribute a single flux over categories (nn_limflx = 2) '
+            IF( .NOT.ln_cpl )   CALL ctl_stop( 'sbc_init : the chosen nn_limflx for LIM3 in forced mode cannot be 2' )
+         CASE DEFAULT
+            CALL ctl_stop( 'sbcmod: LIM3 option, nn_limflx, should be between -1 and 2' )
          END SELECT
-      ENDIF
-      !
-      IF( nn_components /= jp_iam_nemo .AND. .NOT. lk_oasis )   &
-         &      CALL ctl_stop( 'sbc_init : OPA-SAS coupled via OASIS, but key_oasis3 disabled' )
-      IF( nn_components == jp_iam_opa .AND. ln_cpl )   &
-         &      CALL ctl_stop( 'sbc_init : OPA-SAS coupled via OASIS, but ln_cpl = T in OPA' )
-      IF( nn_components == jp_iam_opa .AND. ln_mixcpl )   &
-         &      CALL ctl_stop( 'sbc_init : OPA-SAS coupled via OASIS, but ln_mixcpl = T in OPA' )
-      IF( ln_cpl .AND. .NOT. lk_oasis )    &
-         &      CALL ctl_stop( 'sbc_init : OASIS-coupled atmosphere model, but key_oasis3 disabled' )
-      IF( ln_mixcpl .AND. .NOT. lk_oasis )    &
-         &      CALL ctl_stop( 'the forced-coupled mixed mode (ln_mixcpl) requires the cpp key key_oasis3' )
-      IF( ln_mixcpl .AND. .NOT. ln_cpl )    &
-         &      CALL ctl_stop( 'the forced-coupled mixed mode (ln_mixcpl) requires ln_cpl = T' )
-      IF( ln_mixcpl .AND. nn_components /= jp_iam_nemo )    &
-         &      CALL ctl_stop( 'the forced-coupled mixed mode (ln_mixcpl) is not yet working with sas-opa coupling via oasis' )
-
-      !                              ! allocate sbc arrays
+      ELSE                             ! other sea-ice model
+         IF( nn_limflx >= 0  )   CALL ctl_warn( 'sbc_init : multi-category flux option (nn_limflx) only available in LIM3' )
+      ENDIF
+      !
+      !                       !**  allocate and set required variables
+      !
+      !                             !* allocate sbc arrays
       IF( sbc_oce_alloc() /= 0 )   CALL ctl_stop( 'sbc_init : unable to allocate sbc_oce arrays' )
-
-      !                          ! Checks:
-      IF( .NOT. ln_isf ) THEN                      ! variable initialisation if no ice shelf
+      !
+      IF( .NOT.ln_isf ) THEN        !* No ice-shelf in the domain : allocate and set to zero
          IF( sbc_isf_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'sbc_init : unable to allocate sbc_isf arrays' )
-         fwfisf  (:,:)   = 0.0_wp ; fwfisf_b  (:,:)   = 0.0_wp
-         risf_tsc(:,:,:) = 0.0_wp ; risf_tsc_b(:,:,:) = 0.0_wp
+         fwfisf  (:,:)   = 0._wp   ;   risf_tsc  (:,:,:) = 0._wp
+         fwfisf_b(:,:)   = 0._wp   ;   risf_tsc_b(:,:,:) = 0._wp
       END IF
-      IF( nn_ice == 0 .AND. nn_components /= jp_iam_opa )   fr_i(:,:) = 0._wp    ! no ice in the domain, ice fraction is always zero
-
-      sfx(:,:) = 0._wp                             ! the salt flux due to freezing/melting will be computed (i.e. will be non-zero)
-      !                                            ! only if sea-ice is present
-
-      fmmflx(:,:) = 0._wp                          ! freezing-melting array initialisation
-
-      taum(:,:) = 0._wp                            ! Initialise taum for use in gls in case of reduced restart
-
-      !                                            ! restartability
-      IF( ( nn_ice == 2 .OR. nn_ice ==3 ) .AND. .NOT.( ln_blk .OR. ln_cpl ) )   &
-         &   CALL ctl_stop( 'LIM sea-ice model requires a bulk formulation or coupled configuration' )
-      IF( nn_ice == 4 .AND. .NOT.( ln_blk .OR. ln_cpl ) )   &
-         &   CALL ctl_stop( 'CICE sea-ice model requires ln_blk or ln_cpl' )
-      IF( nn_ice == 4 .AND. lk_agrif )   &
-         &   CALL ctl_stop( 'CICE sea-ice model not currently available with AGRIF' )
-      IF( ( nn_ice == 3 .OR. nn_ice == 4 ) .AND. nn_ice_embd == 0 )   &
-         &   CALL ctl_stop( 'LIM3 and CICE sea-ice models require nn_ice_embd = 1 or 2' )
-      IF( ( nn_ice /= 3 ) .AND. ( nn_limflx >= 0 ) )   &
-         &   WRITE(numout,*) 'The nn_limflx>=0 option has no effect if sea ice model is not LIM3'
-      IF( ( nn_ice == 3 ) .AND. ( ln_cpl ) .AND. ( ( nn_limflx == -1 ) .OR. ( nn_limflx == 1 ) ) )   &
-         &   CALL ctl_stop( 'The chosen nn_limflx for LIM3 in coupled mode must be 0 or 2' )
-      IF( ( nn_ice == 3 ) .AND. ( .NOT. ln_cpl ) .AND. ( nn_limflx == 2 ) )   &
-         &   CALL ctl_stop( 'The chosen nn_limflx for LIM3 in forced mode cannot be 2' )
-
-      IF( ln_dm2dc )   nday_qsr = -1   ! initialisation flag
-
-      IF( ln_dm2dc .AND. .NOT.( ln_flx .OR. ln_blk ) .AND. nn_components /= jp_iam_opa )   &
-         &   CALL ctl_stop( 'diurnal cycle into qsr field from daily values requires a flux or the bulk formulation' )
-
-      !                          ! Choice of the Surface Boudary Condition (set nsbc)
-      ll_purecpl = ln_cpl .AND. .NOT. ln_mixcpl
-      !
+      IF( nn_ice == 0 ) THEN        !* No sea-ice in the domain : ice fraction is always zero
+         IF( nn_components /= jp_iam_opa )   fr_i(:,:) = 0._wp    ! except for OPA in SAS-OPA coupled case
+      ENDIF
+      !
+      sfx   (:,:) = 0._wp           !* salt flux due to freezing/melting
+      fmmflx(:,:) = 0._wp           !* freezing minus melting flux
+
+      taum(:,:) = 0._wp             !* wind stress module (needed in GLS in case of reduced restart)
+
+
+      IF( ln_dm2dc ) THEN           !* daily mean to diurnal cycle
+         nday_qsr = -1   ! allow initialization at the 1st call
+         IF( .NOT.( ln_flx .OR. ln_blk ) .AND. nn_components /= jp_iam_opa )   &
+            &   CALL ctl_stop( 'qsr diurnal cycle from daily values requires a flux or bulk formulation' )
+      ENDIF
+
+      !                             !* Choice of the Surface Boudary Condition (set nsbc)
+      !
+      ll_purecpl  = ln_cpl .AND. .NOT.ln_mixcpl
+      ll_opa      = nn_components == jp_iam_opa
+      ll_not_nemo = nn_components /= jp_iam_nemo
       icpt = 0
+      !
       IF( ln_ana          ) THEN   ;   nsbc = jp_ana     ; icpt = icpt + 1   ;   ENDIF       ! analytical           formulation
       IF( ln_flx          ) THEN   ;   nsbc = jp_flx     ; icpt = icpt + 1   ;   ENDIF       ! flux                 formulation
@@ -225,28 +259,29 @@
       IF( ll_purecpl      ) THEN   ;   nsbc = jp_purecpl ; icpt = icpt + 1   ;   ENDIF       ! Pure Coupled         formulation
       IF( cp_cfg == 'gyre') THEN   ;   nsbc = jp_gyre                        ;   ENDIF       ! GYRE analytical      formulation
-      IF( nn_components == jp_iam_opa )   &
-         &                  THEN   ;   nsbc = jp_none    ; icpt = icpt + 1   ;   ENDIF       ! opa coupling via SAS module
-      !
-      IF( icpt /= 1 )    CALL ctl_stop( 'sbc_init: choose ONE and only ONE sbc option' )
-      !
-      IF(lwp) THEN
+      IF( ll_opa          ) THEN   ;   nsbc = jp_none    ; icpt = icpt + 1   ;   ENDIF       ! opa coupling via SAS module
+      !
+      IF( icpt /= 1 )    CALL ctl_stop( 'sbc_init : choose ONE and only ONE sbc option' )
+      !
+      IF(lwp) THEN                     !- print the choice of surface flux formulation
          WRITE(numout,*)
          SELECT CASE( nsbc )
-         CASE( jp_gyre    )   ;   WRITE(numout,*) '   GYRE analytical formulation'
-         CASE( jp_ana     )   ;   WRITE(numout,*) '   analytical 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_none    )   ;   WRITE(numout,*) '   OPA coupled to SAS via oasis'
-            IF( ln_mixcpl )       WRITE(numout,*) '       + forced-coupled mixed formulation'
+         CASE( jp_gyre    )   ;   WRITE(numout,*) '      ===>>   GYRE analytical formulation'
+         CASE( jp_ana     )   ;   WRITE(numout,*) '      ===>>   analytical 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'
          END SELECT
-         IF( nn_components/= jp_iam_nemo )  &
-            &                     WRITE(numout,*) '       + OASIS coupled SAS'
-      ENDIF
-      !
-      IF( lk_oasis )   CALL sbc_cpl_init (nn_ice)   ! OASIS initialisation. must be done before: (1) first time step
-      !                                             !                                            (2) the use of nn_fsbc
+         IF( ll_not_nemo  )       WRITE(numout,*) '                  + OASIS coupled SAS'
+      ENDIF
+      !
+      !                             !* OASIS initialization
+      !
+      IF( lk_oasis )   CALL sbc_cpl_init( nn_ice )   ! Must be done before: (1) first time step
+      !                                              !                      (2) the use of nn_fsbc
       !     nn_fsbc initialization if OPA-SAS coupling via OASIS
-      !     sas model time step has to be declared in OASIS (mandatory) -> nn_fsbc has to be modified accordingly
+      !     SAS time-step has to be declared in OASIS (mandatory) -> nn_fsbc has to be modified accordingly
       IF( nn_components /= jp_iam_nemo ) THEN
          IF( nn_components == jp_iam_opa )   nn_fsbc = cpl_freq('O_SFLX') / NINT(rdt)
@@ -260,7 +295,8 @@
       ENDIF
       !
+      !                             !* check consistency between model timeline and nn_fsbc
       IF( MOD( nitend - nit000 + 1, nn_fsbc) /= 0 .OR.   &
           MOD( nstock             , nn_fsbc) /= 0 ) THEN
-         WRITE(ctmp1,*) 'experiment length (', nitend - nit000 + 1, ') or nstock (', nstock,   &
+         WRITE(ctmp1,*) 'sbc_init : experiment length (', nitend - nit000 + 1, ') or nstock (', nstock,   &
             &           ' is NOT a multiple of nn_fsbc (', nn_fsbc, ')'
          CALL ctl_stop( ctmp1, 'Impossible to properly do model restart' )
@@ -268,18 +304,23 @@
       !
       IF( MOD( rday, REAL(nn_fsbc, wp) * rdt ) /= 0 )   &
-         &  CALL ctl_warn( 'nn_fsbc is NOT a multiple of the number of time steps in a day' )
-      !
-      IF( ln_dm2dc .AND. ( ( NINT(rday) / ( nn_fsbc * NINT(rdt) ) )  < 8 ) )   &
-         &   CALL ctl_warn( 'diurnal cycle for qsr: the sampling of the diurnal cycle is too small...' )
-      !
-                          CALL sbc_ssm_init               ! Sea-surface mean fields initialisation
-      !
-      IF( ln_ssr      )   CALL sbc_ssr_init               ! Sea-Surface Restoring initialisation
-      !
-                          CALL sbc_rnf_init               ! Runof initialisation
-      !
-      IF( nn_ice == 3 )   CALL sbc_lim_init               ! LIM3 initialisation
-      !
-      IF( nn_ice == 4 )   CALL cice_sbc_init( nsbc )      ! CICE initialisation
+         &  CALL ctl_warn( 'sbc_init : nn_fsbc is NOT a multiple of the number of time steps in a day' )
+      !
+      IF( ln_dm2dc .AND. NINT(rday) / ( nn_fsbc * NINT(rdt) ) < 8  )   &
+         &   CALL ctl_warn( 'sbc_init : diurnal cycle for qsr: the sampling of the diurnal cycle is too small...' )
+      !
+   
+      !                       !**  associated modules : initialization
+      !
+                          CALL sbc_ssm_init            ! Sea-surface mean fields initialization
+      !
+      IF( ln_blk      )   CALL sbc_blk_init            ! bulk formulae initialization
+
+      IF( ln_ssr      )   CALL sbc_ssr_init            ! Sea-Surface Restoring initialization
+      !
+                          CALL sbc_rnf_init            ! Runof initialization
+      !
+      IF( nn_ice == 3 )   CALL sbc_lim_init            ! LIM3 initialization
+      !
+      IF( nn_ice == 4 )   CALL cice_sbc_init( nsbc )   ! CICE initialization
       !
    END SUBROUTINE sbc_init
@@ -337,7 +378,7 @@
       SELECT CASE( nsbc )                                ! Compute ocean surface boundary condition
       !                                                  ! (i.e. utau,vtau, qns, qsr, emp, sfx)
-      CASE( jp_gyre    )   ;   CALL sbc_gyre  ( kt )                     ! analytical formulation : GYRE configuration
-      CASE( jp_ana     )   ;   CALL sbc_ana   ( kt )                     ! analytical formulation : uniform sbc
-      CASE( jp_flx     )   ;   CALL sbc_flx   ( kt )                     ! flux formulation
+      CASE( jp_gyre    )   ;   CALL sbc_gyre   ( kt )                    ! analytical formulation : GYRE configuration
+      CASE( jp_ana     )   ;   CALL sbc_ana    ( kt )                    ! analytical formulation : uniform sbc
+      CASE( jp_flx     )   ;   CALL sbc_flx    ( kt )                    ! flux formulation
       CASE( jp_blk     )
          IF( ll_sas    )       CALL sbc_cpl_rcv( kt, nn_fsbc, nn_ice )   ! OPA-SAS coupling: SAS receiving fields from OPA
Index: branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcrnf.F90
===================================================================
--- branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcrnf.F90	(revision 7162)
+++ branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcrnf.F90	(revision 7163)
@@ -279,6 +279,6 @@
       IF(lwp) THEN
          WRITE(numout,*)
-         WRITE(numout,*) 'sbc_rnf : runoff '
-         WRITE(numout,*) '~~~~~~~ '
+         WRITE(numout,*) 'sbc_rnf_init : runoff '
+         WRITE(numout,*) '~~~~~~~~~~~~ '
          WRITE(numout,*) '   Namelist namsbc_rnf'
          WRITE(numout,*) '      specific river mouths treatment            ln_rnf_mouth = ', ln_rnf_mouth
@@ -296,9 +296,9 @@
          IF(lwp) WRITE(numout,*) '          runoffs inflow read in a file'
          IF( ierror > 0 ) THEN
-            CALL ctl_stop( 'sbc_rnf: unable to allocate sf_rnf structure' )   ;   RETURN
+            CALL ctl_stop( 'sbc_rnf_init: unable to allocate sf_rnf structure' )   ;   RETURN
          ENDIF
          ALLOCATE( sf_rnf(1)%fnow(jpi,jpj,1)   )
          IF( sn_rnf%ln_tint ) ALLOCATE( sf_rnf(1)%fdta(jpi,jpj,1,2) )
-         CALL fld_fill( sf_rnf, (/ sn_rnf /), cn_dir, 'sbc_rnf_init', 'read runoffs data', 'namsbc_rnf' )
+         CALL fld_fill( sf_rnf, (/ sn_rnf /), cn_dir, 'sbc_rnf_init', 'read runoffs data', 'namsbc_rnf', no_print )
       ENDIF
       !
@@ -312,5 +312,5 @@
          ALLOCATE( sf_t_rnf(1)%fnow(jpi,jpj,1)   )
          IF( sn_t_rnf%ln_tint ) ALLOCATE( sf_t_rnf(1)%fdta(jpi,jpj,1,2) )
-         CALL fld_fill (sf_t_rnf, (/ sn_t_rnf /), cn_dir, 'sbc_rnf_init', 'read runoff temperature data', 'namsbc_rnf' )
+         CALL fld_fill (sf_t_rnf, (/ sn_t_rnf /), cn_dir, 'sbc_rnf_init', 'read runoff temperature data', 'namsbc_rnf', no_print )
       ENDIF
       !
@@ -324,5 +324,5 @@
          ALLOCATE( sf_s_rnf(1)%fnow(jpi,jpj,1)   )
          IF( sn_s_rnf%ln_tint ) ALLOCATE( sf_s_rnf(1)%fdta(jpi,jpj,1,2) )
-         CALL fld_fill (sf_s_rnf, (/ sn_s_rnf /), cn_dir, 'sbc_rnf_init', 'read runoff salinity data', 'namsbc_rnf' )
+         CALL fld_fill (sf_s_rnf, (/ sn_s_rnf /), cn_dir, 'sbc_rnf_init', 'read runoff salinity data', 'namsbc_rnf', no_print )
       ENDIF
       !
@@ -452,5 +452,5 @@
             DO WHILE( nkrnf /= jpkm1 .AND. gdepw_1d(nkrnf+1) < rn_hrnf )   ;   nkrnf = nkrnf + 1
             END DO
-            IF( ln_sco )   CALL ctl_warn( 'sbc_rnf: number of levels over which Kz is increased is computed for zco...' )
+            IF( ln_sco )   CALL ctl_warn( 'sbc_rnf_init: number of levels over which Kz is increased is computed for zco...' )
          ENDIF
          IF(lwp) WRITE(numout,*)
@@ -499,6 +499,6 @@
       !
       IF(lwp) WRITE(numout,*)
-      IF(lwp) WRITE(numout,*) 'rnf_mouth : river mouth mask'
-      IF(lwp) WRITE(numout,*) '~~~~~~~~~ '
+      IF(lwp) WRITE(numout,*) '   rnf_mouth : river mouth mask'
+      IF(lwp) WRITE(numout,*) '   ~~~~~~~~~ '
       !
       cl_rnfile = TRIM( cn_dir )//TRIM( sn_cnf%clname )
Index: branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcssm.F90
===================================================================
--- branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcssm.F90	(revision 7162)
+++ branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcssm.F90	(revision 7163)
@@ -221,6 +221,5 @@
             !
             IF( zf_sbc /= REAL( nn_fsbc, wp ) ) THEN      ! nn_fsbc has changed between 2 runs
-               IF(lwp) WRITE(numout,*) '~~~~~~~   restart with a change in the frequency of mean ',   &
-                  &                    'from ', zf_sbc, ' to ', nn_fsbc 
+               IF(lwp) WRITE(numout,*) '   restart with a change in the frequency of mean from ', zf_sbc, ' to ', nn_fsbc 
                zcoef = REAL( nn_fsbc - 1, wp ) / zf_sbc 
                ssu_m(:,:) = zcoef * ssu_m(:,:) 
@@ -232,5 +231,5 @@
                frq_m(:,:) = zcoef * frq_m(:,:)
             ELSE
-               IF(lwp) WRITE(numout,*) '~~~~~~~   mean fields read in the ocean restart file'
+               IF(lwp) WRITE(numout,*) '   mean fields read in the ocean restart file'
             ENDIF
          ENDIF
@@ -239,9 +238,9 @@
       IF( .NOT. l_ssm_mean ) THEN   ! default initialisation. needed by lim_istate
          !
-         IF(lwp) WRITE(numout,*) '          default initialisation of ss?_m arrays'
+         IF(lwp) WRITE(numout,*) '   default initialisation of ss._m arrays'
          ssu_m(:,:) = ub(:,:,1)
          ssv_m(:,:) = vb(:,:,1)
          IF( l_useCT )  THEN    ;   sst_m(:,:) = eos_pt_from_ct( tsn(:,:,1,jp_tem), tsn(:,:,1,jp_sal) )
-         ELSE                    ;   sst_m(:,:) = tsn(:,:,1,jp_tem)
+         ELSE                   ;   sst_m(:,:) = tsn(:,:,1,jp_tem)
          ENDIF
          sss_m(:,:) = tsn  (:,:,1,jp_sal)
Index: branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcssr.F90
===================================================================
--- branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcssr.F90	(revision 7162)
+++ branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcssr.F90	(revision 7163)
@@ -199,5 +199,5 @@
          !
          ! fill sf_sst with sn_sst and control print
-         CALL fld_fill( sf_sst, (/ sn_sst /), cn_dir, 'sbc_ssr', 'SST restoring term toward SST data', 'namsbc_ssr' )
+         CALL fld_fill( sf_sst, (/ sn_sst /), cn_dir, 'sbc_ssr', 'SST restoring term toward SST data', 'namsbc_ssr', no_print )
          IF( sf_sst(1)%ln_tint )   ALLOCATE( sf_sst(1)%fdta(jpi,jpj,1,2), STAT=ierror )
          IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_ssr: unable to allocate sf_sst data array' )
@@ -213,5 +213,5 @@
          !
          ! fill sf_sss with sn_sss and control print
-         CALL fld_fill( sf_sss, (/ sn_sss /), cn_dir, 'sbc_ssr', 'SSS restoring term toward SSS data', 'namsbc_ssr' )
+         CALL fld_fill( sf_sss, (/ sn_sss /), cn_dir, 'sbc_ssr', 'SSS restoring term toward SSS data', 'namsbc_ssr', no_print )
          IF( sf_sss(1)%ln_tint )   ALLOCATE( sf_sss(1)%fdta(jpi,jpj,1,2), STAT=ierror )
          IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_ssr: unable to allocate sf_sss data array' )
Index: branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/TRA/traadv.F90
===================================================================
--- branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/TRA/traadv.F90	(revision 7162)
+++ branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/TRA/traadv.F90	(revision 7163)
@@ -255,13 +255,13 @@
          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,   &
+         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,   &
+         CASE( np_FCT     )   ;   WRITE(numout,*) '      ===>>   FCT      scheme is used. Horizontal order: ', nn_fct_h,   &
             &                                                                      ' Vertical   order: ', nn_fct_v
-         CASE( np_FCT_zts )   ;   WRITE(numout,*) '         use 2nd order FCT with ', nn_fct_zts,'vertical sub-timestepping'
-         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_FCT_zts )   ;   WRITE(numout,*) '      ===>>   use 2nd order FCT with ', nn_fct_zts,'vertical sub-timestepping'
+         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/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_fct.F90
===================================================================
--- branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_fct.F90	(revision 7162)
+++ branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_fct.F90	(revision 7163)
@@ -728,8 +728,8 @@
                !
                IF( tmask(ji,jj,jk+1) == 0._wp) THEN   ! Switch to second order centered at bottom
-                  zwd (ji,jj,jk) = 1._wp
-                  zwi (ji,jj,jk) = 0._wp
-                  zws (ji,jj,jk) = 0._wp
-                  zwrm(ji,jj,jk) = 0.5 * ( pt_in(ji,jj,jk-1) + pt_in(ji,jj,jk) )    
+                  zwd (ji,jj,jk) = 2._wp
+                  zwi (ji,jj,jk) = 0._wp     ! 1  ici ou à l'autre
+                  zws (ji,jj,jk) = 0._wp     !    car 1 dans la diag inferieur au fond et superior en surf
+                  zwrm(ji,jj,jk) = 3._wp * pt_in(ji,jj,jk)  
                ENDIF
             END DO
Index: branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_mle.F90
===================================================================
--- branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_mle.F90	(revision 7162)
+++ branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_mle.F90	(revision 7163)
@@ -308,9 +308,9 @@
          WRITE(numout,*)
          IF( ln_mle ) THEN
-            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'
+            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/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/TRA/trabbc.F90
===================================================================
--- branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/TRA/trabbc.F90	(revision 7162)
+++ branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/TRA/trabbc.F90	(revision 7163)
@@ -176,5 +176,5 @@
             ! fill sf_chl with sn_chl and control print
             CALL fld_fill( sf_qgh, (/ sn_qgh /), cn_dir, 'tra_bbc_init',   &
-               &          'bottom temperature boundary condition', 'nambbc' )
+               &          'bottom temperature boundary condition', 'nambbc', no_print )
 
             CALL fld_read( nit000, 1, sf_qgh )                         ! Read qgh data
Index: branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/TRA/trabbl.F90
===================================================================
--- branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/TRA/trabbl.F90	(revision 7162)
+++ branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/TRA/trabbl.F90	(revision 7163)
@@ -519,9 +519,9 @@
          WRITE(numout,*) 'tra_bbl_init : bottom boundary layer initialisation'
          WRITE(numout,*) '~~~~~~~~~~~~'
-         WRITE(numout,*) '       Namelist nambbl : set bbl parameters'
-         WRITE(numout,*) '          diffusive bbl (=1)   or not (=0)    nn_bbl_ldf = ', nn_bbl_ldf
-         WRITE(numout,*) '          advective bbl (=1/2) or not (=0)    nn_bbl_adv = ', nn_bbl_adv
-         WRITE(numout,*) '          diffusive bbl coefficient           rn_ahtbbl  = ', rn_ahtbbl, ' m2/s'
-         WRITE(numout,*) '          advective bbl coefficient           rn_gambbl  = ', rn_gambbl, ' s'
+         WRITE(numout,*) '   Namelist nambbl : set bbl parameters'
+         WRITE(numout,*) '      diffusive bbl (=1)   or not (=0)    nn_bbl_ldf = ', nn_bbl_ldf
+         WRITE(numout,*) '      advective bbl (=1/2) or not (=0)    nn_bbl_adv = ', nn_bbl_adv
+         WRITE(numout,*) '      diffusive bbl coefficient           rn_ahtbbl  = ', rn_ahtbbl, ' m2/s'
+         WRITE(numout,*) '      advective bbl coefficient           rn_gambbl  = ', rn_gambbl, ' s'
       ENDIF
 
Index: branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/TRA/tradmp.F90
===================================================================
--- branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/TRA/tradmp.F90	(revision 7162)
+++ branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/TRA/tradmp.F90	(revision 7163)
@@ -192,5 +192,5 @@
          WRITE(numout,*)
          WRITE(numout,*) 'tra_dmp_init : T and S newtonian relaxation'
-         WRITE(numout,*) '~~~~~~~~~~~'
+         WRITE(numout,*) '~~~~~~~~~~~~'
          WRITE(numout,*) '   Namelist namtra_dmp : set relaxation parameters'
          WRITE(numout,*) '      Apply relaxation   or not       ln_tradmp = ', ln_tradmp
Index: branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/TRA/traldf.F90
===================================================================
--- branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/TRA/traldf.F90	(revision 7162)
+++ branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/TRA/traldf.F90	(revision 7163)
@@ -110,8 +110,7 @@
          WRITE(numout,*)
          WRITE(numout,*) 'tra_ldf_init : lateral tracer diffusive operator'
-         WRITE(numout,*) '~~~~~~~~~~~'
+         WRITE(numout,*) '~~~~~~~~~~~~'
          WRITE(numout,*) '   Namelist namtra_ldf: already read in ldftra module'
          WRITE(numout,*) '      see ldf_tra_init report for lateral mixing parameters'
-         WRITE(numout,*)
       ENDIF
       !                                   ! use of lateral operator or not
@@ -187,11 +186,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/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/TRA/traqsr.F90
===================================================================
--- branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/TRA/traqsr.F90	(revision 7162)
+++ branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/TRA/traqsr.F90	(revision 7163)
@@ -406,5 +406,5 @@
             !                                        ! fill sf_chl with sn_chl and control print
             CALL fld_fill( sf_chl, (/ sn_chl /), cn_dir, 'tra_qsr_init',   &
-               &           'Solar penetration function of read chlorophyll', 'namtra_qsr' )
+               &           'Solar penetration function of read chlorophyll', 'namtra_qsr' , no_print )
          ENDIF
          IF( nqsr == np_RGB ) THEN                 ! constant Chl
Index: branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/TRA/trazdf.F90
===================================================================
--- branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/TRA/trazdf.F90	(revision 7162)
+++ branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/TRA/trazdf.F90	(revision 7163)
@@ -141,6 +141,6 @@
          WRITE(numout,*) 'tra_zdf_init : vertical tracer physics scheme'
          WRITE(numout,*) '~~~~~~~~~~~'
-         IF( nzdf ==  0 )   WRITE(numout,*) '              Explicit time-splitting scheme'
-         IF( nzdf ==  1 )   WRITE(numout,*) '              Implicit (euler backward) scheme'
+         IF( nzdf ==  0 )   WRITE(numout,*) '      ===>>   Explicit time-splitting scheme'
+         IF( nzdf ==  1 )   WRITE(numout,*) '      ===>>   Implicit (euler backward) scheme'
       ENDIF
       !
Index: branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/TRD/trdini.F90
===================================================================
--- branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/TRD/trdini.F90	(revision 7162)
+++ branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/TRD/trdini.F90	(revision 7163)
@@ -57,6 +57,6 @@
       IF(lwp) THEN                  ! control print
          WRITE(numout,*)
-         WRITE(numout,*) ' trd_init : Momentum/Tracers trends'
-         WRITE(numout,*) ' ~~~~~~~~~~'
+         WRITE(numout,*) 'trd_init : Momentum/Tracers trends'
+         WRITE(numout,*) '~~~~~~~~'
          WRITE(numout,*) '   Namelist namtrd : set trends parameters'
          WRITE(numout,*) '      global domain averaged dyn & tra trends   ln_glo_trd  = ', ln_glo_trd
Index: branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfbfr.F90
===================================================================
--- branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfbfr.F90	(revision 7162)
+++ branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfbfr.F90	(revision 7163)
@@ -99,10 +99,4 @@
       !
       IF( nn_timing == 1 )  CALL timing_start('zdf_bfr')
-      !
-      IF( kt == nit000 .AND. lwp ) THEN
-         WRITE(numout,*)
-         WRITE(numout,*) 'zdf_bfr : Set bottom friction coefficient (non-linear case)'
-         WRITE(numout,*) '~~~~~~~~'
-      ENDIF
       !
       IF( nn_bfr == 2 ) THEN                 ! quadratic bottom friction only
@@ -259,5 +253,5 @@
       IF(lwp) WRITE(numout,*)
       IF(lwp) WRITE(numout,*) 'zdf_bfr_init : momentum bottom friction'
-      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~'
+      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
       IF(lwp) WRITE(numout,*) '   Namelist nam_bfr : set bottom friction parameters'
       !
@@ -266,8 +260,8 @@
       CASE( 0 )
          IF(lwp) WRITE(numout,*) '      free-slip '
-         bfrua(:,:) = 0.e0
-         bfrva(:,:) = 0.e0
-         tfrua(:,:) = 0.e0
-         tfrva(:,:) = 0.e0
+         bfrua(:,:) = 0._wp
+         bfrva(:,:) = 0._wp
+         tfrua(:,:) = 0._wp
+         tfrva(:,:) = 0._wp
          !
       CASE( 1 )
@@ -321,5 +315,5 @@
          IF(lwp) WRITE(numout,*) '      log formulation   ln_bfr2d = ', ln_loglayer
          IF(lwp) WRITE(numout,*) '      bottom roughness  rn_bfrz0 [m] = ', rn_bfrz0
-         IF( rn_bfrz0<=0.e0 ) THEN
+         IF( rn_bfrz0 <= 0._wp ) THEN
             WRITE(ctmp1,*) '      bottom roughness must be strictly positive'
             CALL ctl_stop( ctmp1 )
@@ -336,5 +330,5 @@
             IF(lwp) WRITE(numout,*) '      log formulation   ln_tfr2d     = ', ln_loglayer
             IF(lwp) WRITE(numout,*) '      top roughness     rn_tfrz0 [m] = ', rn_tfrz0
-            IF( rn_tfrz0<=0.e0 ) THEN
+            IF( rn_tfrz0 <= 0._wp ) THEN
                WRITE(ctmp1,*) '      top roughness must be strictly positive'
                CALL ctl_stop( ctmp1 )
Index: branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfini.F90
===================================================================
--- branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfini.F90	(revision 7162)
+++ branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfini.F90	(revision 7163)
@@ -66,5 +66,5 @@
       IF(lwp) THEN               !* Parameter print
          WRITE(numout,*)
-         WRITE(numout,*) 'zdf_init: vertical physics'
+         WRITE(numout,*) 'zdf_init : vertical physics'
          WRITE(numout,*) '~~~~~~~~'
          WRITE(numout,*) '   Namelist namzdf : set vertical mixing mixing parameters'
Index: branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90
===================================================================
--- branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90	(revision 7162)
+++ branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90	(revision 7163)
@@ -558,4 +558,23 @@
          zmxlm(:,:,1) = rn_mxl0
       ENDIF
+      
+!!gm  copy from GLS:
+!      ! Set surface roughness length
+!      SELECT CASE ( nn_z0_met )
+!      !
+!      CASE ( 0 )             ! Constant roughness          
+!         zhsro(:,:) = rn_hsro
+!      CASE ( 1 )             ! Standard Charnock formula for surface roughness
+!         zhsro(:,:) = MAX(rsbc_zs1 * ustars2(:,:), rn_hsro)
+!  with:           rsbc_zs1  = rn_charn/grav
+!      CASE ( 2 )             ! Roughness formulae according to Rascle et al., Ocean Modelling (2008)
+!         zdep(:,:)  = 30.*TANH(2.*0.3/(28.*SQRT(MAX(ustars2(:,:),rsmall))))             ! Wave age (eq. 10)
+!         zhsro(:,:) = MAX(rsbc_zs2 * ustars2(:,:) * zdep(:,:)**1.5, rn_hsro) ! zhsro = rn_frac_hs * Hsw (eq. 11)
+!      !
+!      END SELECT
+!!gm end
+
+
+
       !
       DO jk = 2, jpkm1              ! interior value : l=sqrt(2*e/n^2)
@@ -788,4 +807,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'
          rn_mxl0 = rmxl_min
Index: branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90
===================================================================
--- branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90	(revision 7162)
+++ branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90	(revision 7163)
@@ -519,6 +519,6 @@
       IF(lwp) THEN                  ! control print
          WRITE(numout,*)
-         WRITE(numout,*) 'nemo_ctl: Control prints & Benchmark'
-         WRITE(numout,*) '~~~~~~~ '
+         WRITE(numout,*) 'nemo_ctl : Control prints & Benchmark'
+         WRITE(numout,*) '~~~~~~~~'
          WRITE(numout,*) '   Namelist namctl'
          WRITE(numout,*) '      run control (for debugging)     ln_ctl     = ', ln_ctl
@@ -545,6 +545,6 @@
       IF(lwp) THEN                  ! control print
          WRITE(numout,*)
-         WRITE(numout,*) 'namcfg  : configuration initialization through namelist read'
-         WRITE(numout,*) '~~~~~~~ '
+         WRITE(numout,*) 'namcfg : configuration initialization through namelist read'
+         WRITE(numout,*) '~~~~~~ '
          WRITE(numout,*) '   Namelist namcfg'
          WRITE(numout,*) '      configuration name                               cp_cfg  = ', TRIM(cp_cfg)
Index: branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/trc_oce.F90
===================================================================
--- branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/trc_oce.F90	(revision 7162)
+++ branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/trc_oce.F90	(revision 7163)
@@ -110,6 +110,6 @@
       IF(lwp) THEN
          WRITE(numout,*)
-         WRITE(numout,*) 'trc_oce_rgb : Initialisation of the optical look-up table'
-         WRITE(numout,*) '~~~~~~~~~~~ '
+         WRITE(numout,*) '   trc_oce_rgb : Initialisation of the optical look-up table'
+         WRITE(numout,*) '   ~~~~~~~~~~~ '
       ENDIF
       !
