Index: /branches/2017/dev_merge_2017/NEMOGCM/NEMO/LIM_SRC_3/icethd_zdf.F90
===================================================================
--- /branches/2017/dev_merge_2017/NEMOGCM/NEMO/LIM_SRC_3/icethd_zdf.F90	(revision 9118)
+++ /branches/2017/dev_merge_2017/NEMOGCM/NEMO/LIM_SRC_3/icethd_zdf.F90	(revision 9119)
@@ -119,6 +119,6 @@
       ENDIF
       !
-      IF ( ( ln_cndi_U64 .AND. ln_cndi_P07 ) .OR. ( .NOT.ln_cndi_U64 .AND. .NOT.ln_cndi_P07 ) ) THEN
-         CALL ctl_stop( 'ice_thd_zdf_init: choose one and only one formulation for thermal conduction (ln_cndi_U64 or ln_cndi_P07)' )
+      IF ( ( ln_cndi_U64 .AND. ln_cndi_P07 ) .OR. ( .NOT. ln_cndi_U64 .AND. .NOT. ln_cndi_P07 ) ) THEN
+         CALL ctl_stop( 'ice_thd_zdf_init: choose 1 and only 1 formulation for thermal conduction (ln_cndi_U64 or ln_cndi_P07)' )
       ENDIF
       !                             !== set the choice of ice vertical thermodynamic formulation ==!
Index: /branches/2017/dev_merge_2017/NEMOGCM/NEMO/LIM_SRC_3/icevar.F90
===================================================================
--- /branches/2017/dev_merge_2017/NEMOGCM/NEMO/LIM_SRC_3/icevar.F90	(revision 9118)
+++ /branches/2017/dev_merge_2017/NEMOGCM/NEMO/LIM_SRC_3/icevar.F90	(revision 9119)
@@ -236,5 +236,8 @@
       DO jk = 1, nlay_s
          WHERE( v_s(:,:,:) > epsi20 )        !--- icy area
-            t_s(:,:,jk,:) = MAX( -100._wp , MIN( r1_cpic * ( -r1_rhosn * (e_s(:,:,jk,:)/v_s(:,:,:)*zlay_s) + lfus ) , 0._wp ) ) + rt0
+            t_s(:,:,jk,:) = MAX( -100._wp                                                                          ,   &
+                 &               MIN( r1_cpic * ( -r1_rhosn * ( e_s(:,:,jk,:) / v_s(:,:,:) * zlay_s ) + lfus ) ,       &
+                 &                    0._wp                                                                      )   ) &
+                 &        + rt0
          ELSEWHERE                           !--- no ice
             t_s(:,:,jk,:) = rt0
Index: /branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90
===================================================================
--- /branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90	(revision 9118)
+++ /branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90	(revision 9119)
@@ -236,5 +236,6 @@
                            DO ji = 1, jpi
                               IF( mikt(ji,jj) <= jk .and. jk <=  nk_rnf(ji,jj)  ) THEN
-                                 e3t_b(ji,jj,jk) = e3t_b(ji,jj,jk) - zcoef *  (-rnf_b(ji,jj) + rnf(ji,jj))*(e3t_n(ji,jj,jk)/h_rnf(ji,jj)  )*tmask(ji,jj,jk)
+                                 e3t_b(ji,jj,jk) =   e3t_b(ji,jj,jk) - zcoef *  ( - rnf_b(ji,jj) + rnf(ji,jj) ) &
+                                      &          * ( e3t_n(ji,jj,jk) / h_rnf(ji,jj) ) * tmask(ji,jj,jk)
                               ENDIF
                            ENDDO
Index: /branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90
===================================================================
--- /branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90	(revision 9118)
+++ /branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90	(revision 9119)
@@ -1980,12 +1980,14 @@
 #endif
       ! outputs
-      IF( srcv(jpr_cal)%laction )    CALL iom_put('hflx_cal_cea' , - frcv(jpr_cal)%z3(:,:,1) * lfus                                  ) ! latent heat from calving
-      IF( srcv(jpr_icb)%laction )    CALL iom_put('hflx_icb_cea' , - frcv(jpr_icb)%z3(:,:,1) * lfus                                  ) ! latent heat from icebergs melting
-      IF( iom_use('hflx_snow_cea') ) CALL iom_put('hflx_snow_cea',  sprecip(:,:) * ( zcptsnw(:,:) - Lfus )                           ) ! heat flux from snow (cell average)
-      IF( iom_use('hflx_rain_cea') ) CALL iom_put('hflx_rain_cea',( tprecip(:,:) - sprecip(:,:) ) * zcptrain(:,:)                    ) ! heat flux from rain (cell average)
-      IF( iom_use('hflx_evap_cea') ) CALL iom_put('hflx_evap_cea',(frcv(jpr_tevp)%z3(:,:,1) - frcv(jpr_ievp)%z3(:,:,1) * picefr(:,:) & ! heat flux from from evap (cell average)
-         &                                                        ) * zcptn(:,:) * tmask(:,:,1) )
-      IF( iom_use('hflx_snow_ao_cea') ) CALL iom_put('hflx_snow_ao_cea',sprecip(:,:) * (zcptsnw(:,:) - Lfus) * (1._wp - zsnw(:,:))   ) ! heat flux from snow (over ocean)
-      IF( iom_use('hflx_snow_ai_cea') ) CALL iom_put('hflx_snow_ai_cea',sprecip(:,:) * (zcptsnw(:,:) - Lfus) *          zsnw(:,:)    ) ! heat flux from snow (over ice)
+      IF ( srcv(jpr_cal)%laction       ) CALL iom_put('hflx_cal_cea'    , - frcv(jpr_cal)%z3(:,:,1) * lfus )                       ! latent heat from calving
+      IF ( srcv(jpr_icb)%laction       ) CALL iom_put('hflx_icb_cea'    , - frcv(jpr_icb)%z3(:,:,1) * lfus )                       ! latent heat from icebergs melting
+      IF ( iom_use('hflx_rain_cea')    ) CALL iom_put('hflx_rain_cea'   , ( tprecip(:,:) - sprecip(:,:) ) * zcptrain(:,:) )        ! heat flux from rain (cell average)
+      IF ( iom_use('hflx_evap_cea')    ) CALL iom_put('hflx_evap_cea'   , ( frcv(jpr_tevp)%z3(:,:,1) - frcv(jpr_ievp)%z3(:,:,1) &
+           &                                                              * picefr(:,:) ) * zcptn(:,:) * tmask(:,:,1) )            ! heat flux from evap (cell average)
+      IF ( iom_use('hflx_snow_cea')    ) CALL iom_put('hflx_snow_cea'   , sprecip(:,:) * ( zcptsnw(:,:) - Lfus )   )               ! heat flux from snow (cell average)
+      IF ( iom_use('hflx_snow_ao_cea') ) CALL iom_put('hflx_snow_ao_cea', sprecip(:,:) * ( zcptsnw(:,:) - Lfus ) &
+           &                                                              * ( 1._wp - zsnw(:,:) )                  )               ! heat flux from snow (over ocean)
+      IF ( iom_use('hflx_snow_ai_cea') ) CALL iom_put('hflx_snow_ai_cea', sprecip(:,:) * ( zcptsnw(:,:) - Lfus ) &
+           &                                                              *           zsnw(:,:)                    )               ! heat flux from snow (over ice)
       ! note: hflx for runoff and iceshelf are done in sbcrnf and sbcisf resp.
       !
Index: /branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/TRA/trasbc.F90
===================================================================
--- /branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/TRA/trasbc.F90	(revision 9118)
+++ /branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/TRA/trasbc.F90	(revision 9119)
@@ -129,5 +129,6 @@
                   sbc_tsc(ji,jj,jp_tem) = r1_rau0_rcp * qns(ji,jj)   ! non solar heat flux
                ELSE IF ( sshn(ji,jj) + ht_0(ji,jj) >  rn_wdmin1 ) THEN
-                  sbc_tsc(ji,jj,jp_tem) = r1_rau0_rcp * qns(ji,jj) * (tanh(5._wp*( ( sshn(ji,jj) + ht_0(ji,jj) -  rn_wdmin1 )*r_rn_wdmin1)) )
+                  sbc_tsc(ji,jj,jp_tem) = r1_rau0_rcp * qns(ji,jj) &
+                       &                * tanh ( 5._wp * ( ( sshn(ji,jj) + ht_0(ji,jj) -  rn_wdmin1 ) * r_rn_wdmin1 ) )
                ELSE
                   sbc_tsc(ji,jj,jp_tem) = 0._wp
Index: /branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/TRD/trddyn.F90
===================================================================
--- /branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/TRD/trddyn.F90	(revision 9118)
+++ /branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/TRD/trddyn.F90	(revision 9119)
@@ -156,6 +156,8 @@
                                       ikbu = mbku(ji,jj)          ! deepest ocean u- & v-levels
                                           ikbv = mbkv(ji,jj)
-                                          z3dx(ji,jj,jk) = 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) )*un(ji,jj,ikbu)/e3u_n(ji,jj,ikbu)
-                                          z3dy(ji,jj,jk) = 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) )*vn(ji,jj,ikbv)/e3v_n(ji,jj,ikbv)
+                                          z3dx(ji,jj,jk) = 0.5 * ( rCdU_bot(ji+1,jj) + rCdU_bot(ji,jj) ) & 
+                                               &         * un(ji,jj,ikbu) / e3u_n(ji,jj,ikbu)
+                                          z3dy(ji,jj,jk) = 0.5 * ( rCdU_bot(ji,jj+1) + rCdU_bot(ji,jj) ) &
+                                               &         * vn(ji,jj,ikbv) / e3v_n(ji,jj,ikbv)
                                     END DO
                                  END DO
Index: /branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfosm.F90
===================================================================
--- /branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfosm.F90	(revision 9118)
+++ /branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfosm.F90	(revision 9119)
@@ -362,5 +362,6 @@
               ! It could represent the effects of the spread of wave directions
               ! around the mean wind. The effect of this adjustment needs to be tested.
-              zustke(ji,jj) = MAX (1.0 * ( zcos_wind(ji,jj) * ut0sd(ji,jj) + zsin_wind(ji,jj)  * vt0sd(ji,jj) ), zustar(ji,jj)/(0.45*0.45 ))
+              zustke(ji,jj) = MAX ( 1.0 * ( zcos_wind(ji,jj) * ut0sd(ji,jj ) + zsin_wind(ji,jj)  * vt0sd(ji,jj) ), &
+                   &                zustar(ji,jj) / ( 0.45 * 0.45 )                                                  )
               dstokes(ji,jj) = MAX(zfac * hsw(ji,jj)*hsw(ji,jj) / ( MAX(zustke(ji,jj)*wmp(ji,jj), 1.0e-7 ) ), 5.0e-1) !rn_osm_dstokes !
            END DO
@@ -454,6 +455,6 @@
                            &            + 0.135 * zla(ji,jj) * zwstrl(ji,jj)**3/hbl(ji,jj) )
 
-                      zvel_max = - ( 1.0 + 1.0 * ( zwstrl(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * rn_rdt / hbl(ji,jj) ) * zwb_ent(ji,jj) / &
-                           &       ( zwstrl(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird
+                      zvel_max =  - ( 1.0 + 1.0 * ( zwstrl(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * rn_rdt / hbl(ji,jj) ) &
+                           &   * zwb_ent(ji,jj) / ( zwstrl(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird
 ! Entrainment including component due to shear turbulence. Modified Langmuir component, but gives same result for La=0.3 For testing uncomment.
 !                      zwb_ent(ji,jj) = -( 2.0 * 0.2 * zwbav(ji,jj) &
@@ -471,7 +472,8 @@
                       IF ( zzdhdt < 0._wp ) THEN
                       ! For long timsteps factor in brackets slows the rapid collapse of the OSBL
-                          zpert = 2.0 * ( 1.0 + 2.0 * zwstrl(ji,jj) * rn_rdt / hbl(ji,jj) ) * zwstrl(ji,jj)**2 / hbl(ji,jj)
+                         zpert   = 2.0 * ( 1.0 + 2.0 * zwstrl(ji,jj) * rn_rdt / hbl(ji,jj) ) * zwstrl(ji,jj)**2 / hbl(ji,jj)
                       ELSE
-                          zpert = 2.0 * ( 1.0 + 2.0 * zwstrl(ji,jj) * rn_rdt / hbl(ji,jj) ) * zwstrl(ji,jj)**2 / hbl(ji,jj) + MAX( zdb_bl(ji,jj),0.0 )
+                         zpert   = 2.0 * ( 1.0 + 2.0 * zwstrl(ji,jj) * rn_rdt / hbl(ji,jj) ) * zwstrl(ji,jj)**2 / hbl(ji,jj) &
+                              &  + MAX( zdb_bl(ji,jj), 0.0 )
                       ENDIF
                       zzdhdt = 2.0 * zzdhdt / zpert
@@ -504,59 +506,61 @@
       DO jj = 2, jpjm1
          DO ji = 2, jpim1
-           IF ( ibld(ji,jj) - imld(ji,jj) > 1 ) THEN
+            IF ( ibld(ji,jj) - imld(ji,jj) > 1 ) THEN
 !
 ! If boundary layer changes by more than one level, need to check for stable layers between initial and final depths.
 !
-              zhbl_s = hbl(ji,jj)
-              jm = imld(ji,jj)
-              zthermal = rab_n(ji,jj,1,jp_tem)
-              zbeta = rab_n(ji,jj,1,jp_sal)
-              IF ( lconv(ji,jj) ) THEN
+               zhbl_s = hbl(ji,jj)
+               jm = imld(ji,jj)
+               zthermal = rab_n(ji,jj,1,jp_tem)
+               zbeta = rab_n(ji,jj,1,jp_sal)
+               IF ( lconv(ji,jj) ) THEN
 !unstable
-                 zvel_max = -( 1.0 + 1.0 * ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * rn_rdt / hbl(ji,jj) ) * zwb_ent(ji,jj) / &
-                      &      ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird
-
-                 DO jk = imld(ji,jj), ibld(ji,jj)
-                    zdb = MAX( grav * ( zthermal * ( zt_bl(ji,jj) - tsn(ji,jj,jm,jp_tem) ) &
-                         & - zbeta * ( zs_bl(ji,jj) - tsn(ji,jj,jm,jp_sal) ) ), 0.0 ) + zvel_max
-
-                    zhbl_s = zhbl_s + MIN( - zwb_ent(ji,jj) / zdb * rn_rdt / FLOAT(ibld(ji,jj)-imld(ji,jj) ), e3w_n(ji,jj,jk) )
-                    zhbl_s = MIN(zhbl_s, ht_n(ji,jj))
-
-                    IF ( zhbl_s >= gdepw_n(ji,jj,jm+1) ) jm = jm + 1
-                 END DO
-                 hbl(ji,jj) = zhbl_s
-                 ibld(ji,jj) = jm
-                 hbli(ji,jj) = hbl(ji,jj)
-             ELSE
+                  zvel_max =  - ( 1.0 + 1.0 * ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * rn_rdt / hbl(ji,jj) ) &
+                       &   * zwb_ent(ji,jj) / ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird
+
+                  DO jk = imld(ji,jj), ibld(ji,jj)
+                     zdb = MAX( grav * ( zthermal * ( zt_bl(ji,jj) - tsn(ji,jj,jm,jp_tem) ) &
+                          & - zbeta * ( zs_bl(ji,jj) - tsn(ji,jj,jm,jp_sal) ) ), 0.0 ) + zvel_max
+
+                     zhbl_s = zhbl_s + MIN( - zwb_ent(ji,jj) / zdb * rn_rdt / FLOAT(ibld(ji,jj)-imld(ji,jj) ), e3w_n(ji,jj,jk) )
+                     zhbl_s = MIN(zhbl_s, ht_n(ji,jj))
+
+                     IF ( zhbl_s >= gdepw_n(ji,jj,jm+1) ) jm = jm + 1
+                  END DO
+                  hbl(ji,jj) = zhbl_s
+                  ibld(ji,jj) = jm
+                  hbli(ji,jj) = hbl(ji,jj)
+               ELSE
 ! stable
-                 DO jk = imld(ji,jj), ibld(ji,jj)
-                    zdb = MAX( grav * ( zthermal * ( zt_bl(ji,jj) - tsn(ji,jj,jm,jp_tem) ) - zbeta * ( zs_bl(ji,jj) - tsn(ji,jj,jm,jp_sal) ) ), 0.0 ) + &
-             &       2.0 * zwstrl(ji,jj)**2 / zhbl_s
-
-                    zhbl_s = zhbl_s +  (&
-                         & 0.32 * ( hbli(ji,jj) / zhbl_s -1.0 ) * zwstrl(ji,jj)**3 / hbli(ji,jj) + &
-            &       ( ( 0.32 / 3.0 ) * EXP( -2.5 * ( hbli(ji,jj) / zhbl_s -1.0 )  ) - &
-            &       ( 0.32 / 3.0  - 0.0485 ) * EXP( -12.5 * ( hbli(ji,jj) / zhbl_s ) ) ) * zwstrl(ji,jj)**3 / hbli(ji,jj) ) / zdb * &
-            &       e3w_n(ji,jj,jk) / zdhdt(ji,jj)  ! ALMG to investigate whether need to include wn here
-
-                    zhbl_s = MIN(zhbl_s, ht_n(ji,jj))
-                    IF ( zhbl_s >= gdepw_n(ji,jj,jm) ) jm = jm + 1
-                 END DO
-                 hbl(ji,jj) = MAX(zhbl_s, gdepw_n(ji,jj,3) )
-                 ibld(ji,jj) = MAX(jm, 3 )
-                 IF ( hbl(ji,jj) > hbli(ji,jj) ) hbli(ji,jj) = hbl(ji,jj)
-             ENDIF   ! IF ( lconv )
-           ELSE
+                  DO jk = imld(ji,jj), ibld(ji,jj)
+                     zdb = MAX( grav * ( zthermal * ( zt_bl(ji,jj) - tsn(ji,jj,jm,jp_tem) )          &
+                          &               - zbeta * ( zs_bl(ji,jj) - tsn(ji,jj,jm,jp_sal) ) ), 0.0 ) &
+                          & + 2.0 * zwstrl(ji,jj)**2 / zhbl_s
+
+                     zhbl_s = zhbl_s +  (                                                                                &
+                          &                     0.32         *                         ( hbli(ji,jj) / zhbl_s -1.0 )     &
+                          &               * zwstrl(ji,jj)**3 / hbli(ji,jj)                                               &
+                          &               + ( ( 0.32 / 3.0 )           * EXP( -  2.5 * ( hbli(ji,jj) / zhbl_s -1.0 ) )   &
+                          &               -   ( 0.32 / 3.0  - 0.0485 ) * EXP( - 12.5 * ( hbli(ji,jj) / zhbl_s      ) ) ) &
+                          &          * zwstrl(ji,jj)**3 / hbli(ji,jj) ) / zdb * e3w_n(ji,jj,jk) / zdhdt(ji,jj)  ! ALMG to investigate whether need to include wn here
+
+                     zhbl_s = MIN(zhbl_s, ht_n(ji,jj))
+                     IF ( zhbl_s >= gdepw_n(ji,jj,jm) ) jm = jm + 1
+                  END DO
+                  hbl(ji,jj) = MAX(zhbl_s, gdepw_n(ji,jj,3) )
+                  ibld(ji,jj) = MAX(jm, 3 )
+                  IF ( hbl(ji,jj) > hbli(ji,jj) ) hbli(ji,jj) = hbl(ji,jj)
+               ENDIF   ! IF ( lconv )
+            ELSE
 ! change zero or one model level.
-             hbl(ji,jj) = zhbl_t(ji,jj)
-             IF ( lconv(ji,jj) ) THEN
-                 hbli(ji,jj) = hbl(ji,jj)
-             ELSE
-                 hbl(ji,jj) = MAX(hbl(ji,jj), gdepw_n(ji,jj,3) )
-                 IF ( hbl(ji,jj) > hbli(ji,jj) ) hbli(ji,jj) = hbl(ji,jj)
-             ENDIF
-           ENDIF
-           zhbl(ji,jj) = gdepw_n(ji,jj,ibld(ji,jj))
+               hbl(ji,jj) = zhbl_t(ji,jj)
+               IF ( lconv(ji,jj) ) THEN
+                  hbli(ji,jj) = hbl(ji,jj)
+               ELSE
+                  hbl(ji,jj) = MAX(hbl(ji,jj), gdepw_n(ji,jj,3) )
+                  IF ( hbl(ji,jj) > hbli(ji,jj) ) hbli(ji,jj) = hbl(ji,jj)
+               ENDIF
+            ENDIF
+            zhbl(ji,jj) = gdepw_n(ji,jj,ibld(ji,jj))
          END DO
       END DO
@@ -878,7 +882,8 @@
                     zznd_ml = gdepw_n(ji,jj,jk) / zhml(ji,jj)
                     !
-                    zdiffut(ji,jj,jk) = 0.8 * zdifml_sc(ji,jj) * zznd_ml * ( 1.0 - zbeta_d_sc(ji,jj) * zznd_ml )**1.5
+                    zdiffut(ji,jj,jk) = 0.8   * zdifml_sc(ji,jj) * zznd_ml * ( 1.0 - zbeta_d_sc(ji,jj) * zznd_ml    )**1.5
                     !
-                    zviscos(ji,jj,jk) = 0.375 * zvisml_sc(ji,jj) * zznd_ml * ( 1.0 - zbeta_v_sc(ji,jj) * zznd_ml ) * ( 1.0 - 0.5 * zznd_ml**2 )
+                    zviscos(ji,jj,jk) = 0.375 * zvisml_sc(ji,jj) * zznd_ml * ( 1.0 - zbeta_v_sc(ji,jj) * zznd_ml    ) &
+                         &            *                                      ( 1.0 -               0.5 * zznd_ml**2 )
                 END DO
                 ! pycnocline - if present linear profile
@@ -887,5 +892,5 @@
                        zznd_pyc = -( gdepw_n(ji,jj,jk) - zhml(ji,jj) ) / zdh(ji,jj)
                        !
-                       zdiffut(ji,jj,jk) = zdifpyc_sc(ji,jj)*( 1.0 + zznd_pyc )
+                       zdiffut(ji,jj,jk) = zdifpyc_sc(ji,jj) * ( 1.0 + zznd_pyc )
                        !
                        zviscos(ji,jj,jk) = zvispyc_sc(ji,jj) * ( 1.0 + zznd_pyc )
@@ -928,7 +933,7 @@
               DO jk = 2, imld(ji,jj)
                  zznd_d = gdepw_n(ji,jj,jk) / dstokes(ji,jj)
-                 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) +  1.35 * EXP ( -zznd_d ) * ( 1.0 - EXP ( -2.0 * zznd_d ) ) * zsc_wth_1(ji,jj)
+                 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 1.35 * EXP ( -zznd_d ) * ( 1.0 - EXP ( -2.0 * zznd_d ) ) * zsc_wth_1(ji,jj)
                  !
-                 ghams(ji,jj,jk) = ghams(ji,jj,jk) + 1.35 * EXP ( -zznd_d ) * ( 1.0 - EXP ( -2.0 * zznd_d ) ) * zsc_ws_1(ji,jj)
+                 ghams(ji,jj,jk) = ghams(ji,jj,jk) + 1.35 * EXP ( -zznd_d ) * ( 1.0 - EXP ( -2.0 * zznd_d ) ) *  zsc_ws_1(ji,jj)
               END DO ! end jk loop
             ELSE     ! else for if (lconv)
@@ -936,7 +941,9 @@
                DO jk = 2, ibld(ji,jj)
                   zznd_d=gdepw_n(ji,jj,jk) / dstokes(ji,jj)
-                  ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 1.5 * EXP ( -0.9 * zznd_d ) * ( 1.0 - EXP ( -4.0 * zznd_d ) ) * zsc_wth_1(ji,jj)
+                  ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 1.5 * EXP ( -0.9 * zznd_d ) &
+                       &          *                 ( 1.0 - EXP ( -4.0 * zznd_d ) ) * zsc_wth_1(ji,jj)
                   !
-                  ghams(ji,jj,jk) = ghams(ji,jj,jk) + 1.5 * EXP ( -0.9 * zznd_d ) * ( 1.0 - EXP ( -4.0 * zznd_d ) ) * zsc_ws_1(ji,jj)
+                  ghams(ji,jj,jk) = ghams(ji,jj,jk) + 1.5 * EXP ( -0.9 * zznd_d ) &
+                       &          *                 ( 1.0 - EXP ( -4.0 * zznd_d ) ) *  zsc_ws_1(ji,jj)
                END DO
             ENDIF               ! endif for check on lconv
@@ -961,7 +968,10 @@
                 DO jk = 2, imld(ji,jj)
                    zznd_d = gdepw_n(ji,jj,jk) / dstokes(ji,jj)
-                   ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + ( -0.05 * EXP ( -0.4 * zznd_d ) * zsc_uw_1(ji,jj) + 0.00125 * EXP ( -zznd_d ) * zsc_uw_2(ji,jj) ) * ( 1.0 - EXP ( -2.0 * zznd_d ) )
+                   ghamu(ji,jj,jk) = ghamu(ji,jj,jk) +      ( -0.05 * EXP ( -0.4 * zznd_d )   * zsc_uw_1(ji,jj)   &
+                        &          +                        0.00125 * EXP (      - zznd_d )   * zsc_uw_2(ji,jj) ) &
+                        &          *                          ( 1.0 - EXP ( -2.0 * zznd_d ) )
 !
-                   ghamv(ji,jj,jk) = ghamv(ji,jj,jk) - 0.65 * 0.15 * EXP ( -zznd_d ) * ( 1.0 - EXP ( -2.0 * zznd_d ) ) * zsc_vw_1(ji,jj)
+                   ghamv(ji,jj,jk) = ghamv(ji,jj,jk) - 0.65 *  0.15 * EXP (      - zznd_d )                       &
+                        &          *                          ( 1.0 - EXP ( -2.0 * zznd_d ) ) * zsc_vw_1(ji,jj)
                 END DO   ! end jk loop
              ELSE
@@ -969,5 +979,6 @@
                 DO jk = 2, ibld(ji,jj) ! corrected to ibld
                    zznd_d = gdepw_n(ji,jj,jk) / dstokes(ji,jj)
-                   ghamu(ji,jj,jk) = ghamu(ji,jj,jk) - 0.75 * 1.3 * EXP ( -0.5 * zznd_d ) * ( 1.0 - EXP ( -4.0 * zznd_d ) ) * zsc_uw_1(ji,jj)
+                   ghamu(ji,jj,jk) = ghamu(ji,jj,jk) - 0.75 *   1.3 * EXP ( -0.5 * zznd_d )                       &
+                        &                                   * ( 1.0 - EXP ( -4.0 * zznd_d ) ) * zsc_uw_1(ji,jj)
                    ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + 0._wp
                 END DO   ! end jk loop
@@ -980,5 +991,5 @@
        WHERE ( lconv )
           zsc_wth_1 = zwbav * zwth0 * ( 1.0 + EXP ( 0.2 * zhol ) ) / ( zvstr**3 + 0.5 * zwstrc**3 + epsln )
-          zsc_ws_1 = zwbav * zws0 * ( 1.0 + EXP ( 0.2 * zhol ) ) / ( zvstr**3 + 0.5 * zwstrc**3 + epsln )
+          zsc_ws_1  = zwbav * zws0  * ( 1.0 + EXP ( 0.2 * zhol ) ) / ( zvstr**3 + 0.5 * zwstrc**3 + epsln )
        ELSEWHERE
           zsc_wth_1 = 0._wp
@@ -991,16 +1002,18 @@
                 DO jk = 2, imld(ji,jj)
                    zznd_ml = gdepw_n(ji,jj,jk) / zhml(ji,jj)
-            ! calculate turbulent length scale
-                   zl_c = 0.9 * ( 1.0 - EXP ( -7.0 * ( zznd_ml - zznd_ml**3 / 3.0 ) ) ) * ( 1.0 - EXP ( -15.0 * ( 1.1 - zznd_ml ) ) )
-                   zl_l = 2.0 * ( 1.0 - EXP ( -2.0 * ( zznd_ml - zznd_ml**3 / 3.0 ) ) ) * (1.0 -  EXP ( -5.0 * ( 1.0 -zznd_ml ) ) ) * ( 1.0 + dstokes(ji,jj) / zhml (ji,jj) )
+                   ! calculate turbulent length scale
+                   zl_c = 0.9 * ( 1.0 - EXP ( - 7.0 * ( zznd_ml - zznd_ml**3 / 3.0 ) ) )                                           &
+                        &     * ( 1.0 - EXP ( -15.0 * (     1.1 - zznd_ml          ) ) )
+                   zl_l = 2.0 * ( 1.0 - EXP ( - 2.0 * ( zznd_ml - zznd_ml**3 / 3.0 ) ) )                                           &
+                        &     * ( 1.0 - EXP ( - 5.0 * (     1.0 - zznd_ml          ) ) ) * ( 1.0 + dstokes(ji,jj) / zhml (ji,jj) )
                    zl_eps = zl_l + ( zl_c - zl_l ) / ( 1.0 + EXP ( 3.0 * LOG10 ( - zhol(ji,jj) ) ) ) ** (3.0/2.0)
-            ! non-gradient buoyancy terms
+                   ! non-gradient buoyancy terms
                    ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 0.3 * 0.5 * zsc_wth_1(ji,jj) * zl_eps * zhml(ji,jj) / ( 0.15 + zznd_ml )
-                   ghams(ji,jj,jk) = ghams(ji,jj,jk) + 0.3 * 0.5 * zsc_ws_1(ji,jj) * zl_eps * zhml(ji,jj) / ( 0.15 + zznd_ml )
+                   ghams(ji,jj,jk) = ghams(ji,jj,jk) + 0.3 * 0.5 *  zsc_ws_1(ji,jj) * zl_eps * zhml(ji,jj) / ( 0.15 + zznd_ml )
                 END DO
              ELSE
                 DO jk = 2, ibld(ji,jj)
                    ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zsc_wth_1(ji,jj)
-                   ghams(ji,jj,jk) = ghams(ji,jj,jk) + zsc_ws_1(ji,jj)
+                   ghams(ji,jj,jk) = ghams(ji,jj,jk) +  zsc_ws_1(ji,jj)
                 END DO
              ENDIF
@@ -1011,5 +1024,5 @@
        WHERE ( lconv )
           zsc_uw_1 = -zwb0 * zustar**2 * zhml / ( zvstr**3 + 0.5 * zwstrc**3 + epsln )
-          zsc_uw_2 = zwb0 *  zustke * zhml / ( zvstr**3 + 0.5 * zwstrc**3 + epsln )**(2.0/3.0)
+          zsc_uw_2 =  zwb0 * zustke    * zhml / ( zvstr**3 + 0.5 * zwstrc**3 + epsln )**(2.0/3.0)
           zsc_vw_1 = 0._wp
        ELSEWHERE
@@ -1023,5 +1036,7 @@
                 DO jk = 2 , imld(ji,jj)
                    zznd_d = gdepw_n(ji,jj,jk) / dstokes(ji,jj)
-                   ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + 0.3 * 0.5 * ( zsc_uw_1(ji,jj) + 0.125 * EXP( -0.5 * zznd_d ) * ( 1.0 - EXP ( -0.5 * zznd_d) ) * zsc_uw_2(ji,jj) )
+                   ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + 0.3 * 0.5 * ( zsc_uw_1(ji,jj) +   0.125 * EXP( -0.5 * zznd_d )     &
+                        &                                                            * (   1.0 - EXP( -0.5 * zznd_d ) )   &
+                        &                                          * zsc_uw_2(ji,jj)                                    )
                    ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zsc_vw_1(ji,jj)
                 END DO  ! jk loop
@@ -1051,9 +1066,13 @@
                DO jk = 2, imld(ji,jj)
                   zznd_ml=gdepw_n(ji,jj,jk) / zhml(ji,jj)
-                  ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 0.3 * zsc_wth_1(ji,jj) * ( -2.0 + 2.75 * ( ( 1.0 + 0.6 * zznd_ml**4 ) - EXP ( -6.0 *zznd_ml ) ) )  &
-      &           * ( 1.0 -EXP ( -15.0 * ( 1.0 - zznd_ml ) ) )
+                  ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 0.3 * zsc_wth_1(ji,jj)                &
+                       &          * ( -2.0 + 2.75 * (       ( 1.0 + 0.6 * zznd_ml**4 )      &
+                       &                               - EXP(     - 6.0 * zznd_ml    ) ) )  &
+                       &          * ( 1.0 - EXP( - 15.0 * (         1.0 - zznd_ml    ) ) )
                   !
-                  ghams(ji,jj,jk) = ghams(ji,jj,jk) + 0.3 * zsc_ws_1(ji,jj) * ( -2.0 + 2.75 * ( ( 1.0 + 0.6 * zznd_ml**4 ) - EXP ( -6.0 *zznd_ml ) ) )  &
-      &           * ( 1.0 -EXP ( -15.0 * ( 1.0 - zznd_ml ) ) )
+                  ghams(ji,jj,jk) = ghams(ji,jj,jk) + 0.3 * zsc_ws_1(ji,jj)  &
+                       &          * ( -2.0 + 2.75 * (       ( 1.0 + 0.6 * zznd_ml**4 )      &
+                       &                               - EXP(     - 6.0 * zznd_ml    ) ) )  &
+                       &          * ( 1.0 - EXP ( -15.0 * (         1.0 - zznd_ml    ) ) )
                END DO
             ELSE
Index: /branches/2017/dev_merge_2017/NEMOGCM/NEMO/TOP_SRC/AGE/trcnam_age.F90
===================================================================
--- /branches/2017/dev_merge_2017/NEMOGCM/NEMO/TOP_SRC/AGE/trcnam_age.F90	(revision 9118)
+++ /branches/2017/dev_merge_2017/NEMOGCM/NEMO/TOP_SRC/AGE/trcnam_age.F90	(revision 9119)
@@ -61,6 +61,8 @@
          WRITE(numout,*) ' trc_nam_age: Read namage, namelist for Age passive tracer'
          WRITE(numout,*) ' ~~~~~~~'
-         WRITE(numout,*) '  depth over which age tracer reset to zero                              rn_age_depth      = ', rn_age_depth 
-         WRITE(numout,*) '  recip of relax. timescale (s) for age tracer shallower than age_depth  rn_age_kill_rate  = ', rn_age_kill_rate 
+         WRITE(numout,*) '  depth over which age tracer reset to zero                              rn_age_depth      = ', &
+              &          rn_age_depth 
+         WRITE(numout,*) '  recip of relax. timescale (s) for age tracer shallower than age_depth  rn_age_kill_rate  = ', &
+              &          rn_age_kill_rate 
          WRITE(numout,*) ''
       ENDIF
Index: /branches/2017/dev_merge_2017/NEMOGCM/NEMO/TOP_SRC/trcbc.F90
===================================================================
--- /branches/2017/dev_merge_2017/NEMOGCM/NEMO/TOP_SRC/trcbc.F90	(revision 9118)
+++ /branches/2017/dev_merge_2017/NEMOGCM/NEMO/TOP_SRC/trcbc.F90	(revision 9119)
@@ -198,5 +198,6 @@
          WRITE(numout,'(2a)') '   COASTAL BC data repository : ', TRIM(cn_dir_cbc)
          IF( .NOT.ln_rnf .OR. .NOT.ln_linssh ) ln_rnf_ctl = .FALSE.
-         IF( ln_rnf_ctl )  WRITE(numout,'(a)') ' -> Remove runoff dilution effect on tracers with absent river load (ln_rnf_ctl = .TRUE.)' 
+         IF( ln_rnf_ctl )  WRITE(numout,'(a)') &
+              &            ' -> Remove runoff dilution effect on tracers with absent river load (ln_rnf_ctl = .TRUE.)'
          WRITE(numout,*) ' '
          WRITE(numout,'(a,i3)') '   Total tracers to be initialized with OPEN BCs data:', nb_trcobc
@@ -205,12 +206,17 @@
             WRITE(numout,*) '   #trc        NAME        Boundary     Mult.Fact.   OBC Settings'
             DO jn = 1, ntrc
-               IF( ln_trc_obc(jn) )  WRITE(numout, 9001) jn, TRIM( sn_trcobc(jn)%clvar ), 'OBC', rn_trofac(jn), (trcdta_bdy(jn,ib)%cn_obc,ib=1,nb_bdy)
-               IF( .NOT. ln_trc_obc(jn) )  WRITE(numout, 9002) jn, 'Set data to IC and use default condition', (trcdta_bdy(jn,ib)%cn_obc,ib=1,nb_bdy)
+               IF (       ln_trc_obc(jn) )  WRITE(numout, 9001) jn, TRIM( sn_trcobc(jn)%clvar ), 'OBC', rn_trofac(jn), &
+                    &                                           (trcdta_bdy(jn,ib)%cn_obc,ib=1,nb_bdy)
+               IF ( .NOT. ln_trc_obc(jn) )  WRITE(numout, 9002) jn, 'Set data to IC and use default condition'       , &
+                    &                                           (trcdta_bdy(jn,ib)%cn_obc,ib=1,nb_bdy)
             ENDDO
             WRITE(numout,*) ' '
             DO ib = 1, nb_bdy
-                IF(nn_trcdmp_bdy(ib) .EQ. 0) WRITE(numout,9003) '   Boundary ',ib,' -> NO damping of tracers'
-                IF(nn_trcdmp_bdy(ib) .EQ. 1) WRITE(numout,9003) '   Boundary ',ib,' -> damping ONLY for tracers with external data provided'
-                IF(nn_trcdmp_bdy(ib) .EQ. 2) WRITE(numout,9003) '   Boundary ',ib,' -> damping of ALL tracers'
+                IF(nn_trcdmp_bdy(ib) .EQ. 0) WRITE(numout,9003) '   Boundary ', ib, &
+                     &                                          ' -> NO damping of tracers'
+                IF(nn_trcdmp_bdy(ib) .EQ. 1) WRITE(numout,9003) '   Boundary ', ib, &
+                     &                                          ' -> damping ONLY for tracers with external data provided'
+                IF(nn_trcdmp_bdy(ib) .EQ. 2) WRITE(numout,9003) '   Boundary ', ib, &
+                     &                                          ' -> damping of ALL tracers'
                 IF(nn_trcdmp_bdy(ib) .GT. 0) THEN
                    WRITE(numout,9003) '     USE damping parameters from nambdy for boundary ', ib,' : '
