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 9106)
+++ 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
