Index: /NEMO/trunk/src/OCE/SBC/sbcblk.F90
===================================================================
--- /NEMO/trunk/src/OCE/SBC/sbcblk.F90	(revision 12614)
+++ /NEMO/trunk/src/OCE/SBC/sbcblk.F90	(revision 12615)
@@ -639,14 +639,4 @@
       END IF
 
-      !!      CALL iom_put( "Cd_oce", zcd_oce)  ! output value of pure ocean-atm. transfer coef.
-      !!      CALL iom_put( "Ch_oce", zch_oce)  ! output value of pure ocean-atm. transfer coef.
-
-      IF( ABS(rn_zu - rn_zqt) < 0.1_wp ) THEN
-         !! If zu == zt, then ensuring once for all that:
-         t_zu(:,:) = ztpot(:,:)
-         q_zu(:,:) = zqair(:,:)
-      ENDIF
-
-
       !  Turbulent fluxes over ocean  => BULK_FORMULA @ sbcblk_phy.F90
       ! -------------------------------------------------------------
@@ -663,8 +653,8 @@
       ELSE                      !==  BLK formulation  ==!   turbulent fluxes computation
          CALL BULK_FORMULA( rn_zu, ptsk(:,:), pssq(:,:), t_zu(:,:), q_zu(:,:), &
-            &               zcd_oce(:,:), zch_oce(:,:), zce_oce(:,:),         &
-            &               wndm(:,:), zU_zu(:,:), pslp(:,:),                 &
-            &               taum(:,:), psen(:,:), zqla(:,:),                  &
-            &               pEvap=pevp(:,:), prhoa=rhoa(:,:) )
+            &               zcd_oce(:,:), zch_oce(:,:), zce_oce(:,:),          &
+            &               wndm(:,:), zU_zu(:,:), pslp(:,:),                  &
+            &               taum(:,:), psen(:,:), zqla(:,:),                   &
+            &               pEvap=pevp(:,:), prhoa=rhoa(:,:), pfact_evap=rn_efac )
 
          zqla(:,:) = zqla(:,:) * tmask(:,:,1)
@@ -1046,5 +1036,5 @@
       evap_ice (:,:,:) = rn_efac * qla_ice (:,:,:) * z1_rLsub    ! sublimation
       devap_ice(:,:,:) = rn_efac * dqla_ice(:,:,:) * z1_rLsub    ! d(sublimation)/dT
-      zevap    (:,:)   = rn_efac * ( emp(:,:) + tprecip(:,:) )   ! evaporation over ocean
+      zevap    (:,:)   = rn_efac * ( emp(:,:) + tprecip(:,:) )   ! evaporation over ocean  !LB: should we remove rn_efac here???
 
       ! --- evaporation minus precipitation --- !
Index: /NEMO/trunk/src/OCE/SBC/sbcblk_algo_coare3p0.F90
===================================================================
--- /NEMO/trunk/src/OCE/SBC/sbcblk_algo_coare3p0.F90	(revision 12614)
+++ /NEMO/trunk/src/OCE/SBC/sbcblk_algo_coare3p0.F90	(revision 12615)
@@ -194,6 +194,5 @@
       IF( kt == nit000 ) CALL SBCBLK_ALGO_COARE3P0_INIT(l_use_cs, l_use_wl)
 
-      l_zt_equal_zu = .FALSE.
-      IF( ABS(zu - zt) < 0.01_wp )   l_zt_equal_zu = .TRUE.    ! testing "zu == zt" is risky with double precision
+      l_zt_equal_zu = ( ABS(zu - zt) < 0.01_wp ) ! testing "zu == zt" is risky with double precision
       IF( .NOT. l_zt_equal_zu )  ALLOCATE( zeta_t(jpi,jpj) )
 
@@ -396,15 +395,15 @@
       !
       DO_2D_11_11
-         !
-         zw = pwnd(ji,jj)   ! wind speed
-         !
-         ! Charnock's constant, increases with the wind :
-         zgt10 = 0.5 + SIGN(0.5_wp,(zw - 10))  ! If zw<10. --> 0, else --> 1
-         zgt18 = 0.5 + SIGN(0.5_wp,(zw - 18.)) ! If zw<18. --> 0, else --> 1
-         !
-         alfa_charn_3p0(ji,jj) =  (1. - zgt10)*0.011    &    ! wind is lower than 10 m/s
-            &     + zgt10*((1. - zgt18)*(0.011 + (0.018 - 0.011) &
-            &      *(zw - 10.)/(18. - 10.)) + zgt18*( 0.018 ) )    ! Hare et al. (1999)
-         !
+      !
+      zw = pwnd(ji,jj)   ! wind speed
+      !
+      ! Charnock's constant, increases with the wind :
+      zgt10 = 0.5 + SIGN(0.5_wp,(zw - 10))  ! If zw<10. --> 0, else --> 1
+      zgt18 = 0.5 + SIGN(0.5_wp,(zw - 18.)) ! If zw<18. --> 0, else --> 1
+      !
+      alfa_charn_3p0(ji,jj) =  (1. - zgt10)*0.011    &    ! wind is lower than 10 m/s
+         &     + zgt10*((1. - zgt18)*(0.011 + (0.018 - 0.011) &
+         &      *(zw - 10.)/(18. - 10.)) + zgt18*( 0.018 ) )    ! Hare et al. (1999)
+      !
       END_2D
       !
@@ -432,26 +431,26 @@
       !
       DO_2D_11_11
-         !
-         zta = pzeta(ji,jj)
-         !
-         zphi_m = ABS(1. - 15.*zta)**.25    !!Kansas unstable
-         !
-         zpsi_k = 2.*LOG((1. + zphi_m)/2.) + LOG((1. + zphi_m*zphi_m)/2.)   &
-            & - 2.*ATAN(zphi_m) + 0.5*rpi
-         !
-         zphi_c = ABS(1. - 10.15*zta)**.3333                   !!Convective
-         !
-         zpsi_c = 1.5*LOG((1. + zphi_c + zphi_c*zphi_c)/3.) &
-            &     - 1.7320508*ATAN((1. + 2.*zphi_c)/1.7320508) + 1.813799447
-         !
-         zf = zta*zta
-         zf = zf/(1. + zf)
-         zc = MIN(50._wp, 0.35_wp*zta)
-         zstab = 0.5 + SIGN(0.5_wp, zta)
-         !
-         psi_m_coare(ji,jj) = (1. - zstab) * ( (1. - zf)*zpsi_k + zf*zpsi_c ) & ! (zta < 0)
-            &                -   zstab     * ( 1. + 1.*zta     &                ! (zta > 0)
-            &                         + 0.6667*(zta - 14.28)/EXP(zc) + 8.525 )   !     "
-         !
+      !
+      zta = pzeta(ji,jj)
+      !
+      zphi_m = ABS(1. - 15.*zta)**.25    !!Kansas unstable
+      !
+      zpsi_k = 2.*LOG((1. + zphi_m)/2.) + LOG((1. + zphi_m*zphi_m)/2.)   &
+         & - 2.*ATAN(zphi_m) + 0.5*rpi
+      !
+      zphi_c = ABS(1. - 10.15*zta)**.3333                   !!Convective
+      !
+      zpsi_c = 1.5*LOG((1. + zphi_c + zphi_c*zphi_c)/3.) &
+         &     - 1.7320508*ATAN((1. + 2.*zphi_c)/1.7320508) + 1.813799447
+      !
+      zf = zta*zta
+      zf = zf/(1. + zf)
+      zc = MIN(50._wp, 0.35_wp*zta)
+      zstab = 0.5 + SIGN(0.5_wp, zta)
+      !
+      psi_m_coare(ji,jj) = (1. - zstab) * ( (1. - zf)*zpsi_k + zf*zpsi_c ) & ! (zta < 0)
+         &                -   zstab     * ( 1. + 1.*zta     &                ! (zta > 0)
+         &                         + 0.6667*(zta - 14.28)/EXP(zc) + 8.525 )   !     "
+      !
       END_2D
       !
@@ -483,25 +482,25 @@
       !
       DO_2D_11_11
-         !
-         zta = pzeta(ji,jj)
-         !
-         zphi_h = (ABS(1. - 15.*zta))**.5  !! Kansas unstable   (zphi_h = zphi_m**2 when unstable, zphi_m when stable)
-         !
-         zpsi_k = 2.*LOG((1. + zphi_h)/2.)
-         !
-         zphi_c = (ABS(1. - 34.15*zta))**.3333   !! Convective
-         !
-         zpsi_c = 1.5*LOG((1. + zphi_c + zphi_c*zphi_c)/3.) &
-            &    -1.7320508*ATAN((1. + 2.*zphi_c)/1.7320508) + 1.813799447
-         !
-         zf = zta*zta
-         zf = zf/(1. + zf)
-         zc = MIN(50._wp,0.35_wp*zta)
-         zstab = 0.5 + SIGN(0.5_wp, zta)
-         !
-         psi_h_coare(ji,jj) = (1. - zstab) * ( (1. - zf)*zpsi_k + zf*zpsi_c ) &
-            &                -   zstab     * ( (ABS(1. + 2.*zta/3.))**1.5     &
-            &                           + .6667*(zta - 14.28)/EXP(zc) + 8.525 )
-         !
+      !
+      zta = pzeta(ji,jj)
+      !
+      zphi_h = (ABS(1. - 15.*zta))**.5  !! Kansas unstable   (zphi_h = zphi_m**2 when unstable, zphi_m when stable)
+      !
+      zpsi_k = 2.*LOG((1. + zphi_h)/2.)
+      !
+      zphi_c = (ABS(1. - 34.15*zta))**.3333   !! Convective
+      !
+      zpsi_c = 1.5*LOG((1. + zphi_c + zphi_c*zphi_c)/3.) &
+         &    -1.7320508*ATAN((1. + 2.*zphi_c)/1.7320508) + 1.813799447
+      !
+      zf = zta*zta
+      zf = zf/(1. + zf)
+      zc = MIN(50._wp,0.35_wp*zta)
+      zstab = 0.5 + SIGN(0.5_wp, zta)
+      !
+      psi_h_coare(ji,jj) = (1. - zstab) * ( (1. - zf)*zpsi_k + zf*zpsi_c ) &
+         &                -   zstab     * ( (ABS(1. + 2.*zta/3.))**1.5     &
+         &                           + .6667*(zta - 14.28)/EXP(zc) + 8.525 )
+      !
       END_2D
       !
Index: /NEMO/trunk/src/OCE/SBC/sbcblk_algo_coare3p6.F90
===================================================================
--- /NEMO/trunk/src/OCE/SBC/sbcblk_algo_coare3p6.F90	(revision 12614)
+++ /NEMO/trunk/src/OCE/SBC/sbcblk_algo_coare3p6.F90	(revision 12615)
@@ -194,6 +194,5 @@
       IF( kt == nit000 ) CALL SBCBLK_ALGO_COARE3P6_INIT(l_use_cs, l_use_wl)
 
-      l_zt_equal_zu = .FALSE.
-      IF( ABS(zu - zt) < 0.01_wp )   l_zt_equal_zu = .TRUE.    ! testing "zu == zt" is risky with double precision
+      l_zt_equal_zu = ( ABS(zu - zt) < 0.01_wp ) ! testing "zu == zt" is risky with double precision
       IF( .NOT. l_zt_equal_zu )  ALLOCATE( zeta_t(jpi,jpj) )
 
@@ -432,26 +431,26 @@
       !
       DO_2D_11_11
-         !
-         zta = pzeta(ji,jj)
-         !
-         zphi_m = ABS(1. - 15.*zta)**.25    !!Kansas unstable
-         !
-         zpsi_k = 2.*LOG((1. + zphi_m)/2.) + LOG((1. + zphi_m*zphi_m)/2.)   &
-            & - 2.*ATAN(zphi_m) + 0.5*rpi
-         !
-         zphi_c = ABS(1. - 10.15*zta)**.3333                   !!Convective
-         !
-         zpsi_c = 1.5*LOG((1. + zphi_c + zphi_c*zphi_c)/3.) &
-            &     - 1.7320508*ATAN((1. + 2.*zphi_c)/1.7320508) + 1.813799447
-         !
-         zf = zta*zta
-         zf = zf/(1. + zf)
-         zc = MIN(50._wp, 0.35_wp*zta)
-         zstab = 0.5 + SIGN(0.5_wp, zta)
-         !
-         psi_m_coare(ji,jj) = (1. - zstab) * ( (1. - zf)*zpsi_k + zf*zpsi_c ) & ! (zta < 0)
-            &                -   zstab     * ( 1. + 1.*zta     &                ! (zta > 0)
-            &                         + 0.6667*(zta - 14.28)/EXP(zc) + 8.525 )   !     "
-         !
+      !
+      zta = pzeta(ji,jj)
+      !
+      zphi_m = ABS(1. - 15.*zta)**.25    !!Kansas unstable
+      !
+      zpsi_k = 2.*LOG((1. + zphi_m)/2.) + LOG((1. + zphi_m*zphi_m)/2.)   &
+         & - 2.*ATAN(zphi_m) + 0.5*rpi
+      !
+      zphi_c = ABS(1. - 10.15*zta)**.3333                   !!Convective
+      !
+      zpsi_c = 1.5*LOG((1. + zphi_c + zphi_c*zphi_c)/3.) &
+         &     - 1.7320508*ATAN((1. + 2.*zphi_c)/1.7320508) + 1.813799447
+      !
+      zf = zta*zta
+      zf = zf/(1. + zf)
+      zc = MIN(50._wp, 0.35_wp*zta)
+      zstab = 0.5 + SIGN(0.5_wp, zta)
+      !
+      psi_m_coare(ji,jj) = (1. - zstab) * ( (1. - zf)*zpsi_k + zf*zpsi_c ) & ! (zta < 0)
+         &                -   zstab     * ( 1. + 1.*zta     &                ! (zta > 0)
+         &                         + 0.6667*(zta - 14.28)/EXP(zc) + 8.525 )   !     "
+      !
       END_2D
       !
@@ -483,25 +482,25 @@
       !
       DO_2D_11_11
-         !
-         zta = pzeta(ji,jj)
-         !
-         zphi_h = (ABS(1. - 15.*zta))**.5  !! Kansas unstable   (zphi_h = zphi_m**2 when unstable, zphi_m when stable)
-         !
-         zpsi_k = 2.*LOG((1. + zphi_h)/2.)
-         !
-         zphi_c = (ABS(1. - 34.15*zta))**.3333   !! Convective
-         !
-         zpsi_c = 1.5*LOG((1. + zphi_c + zphi_c*zphi_c)/3.) &
-            &    -1.7320508*ATAN((1. + 2.*zphi_c)/1.7320508) + 1.813799447
-         !
-         zf = zta*zta
-         zf = zf/(1. + zf)
-         zc = MIN(50._wp,0.35_wp*zta)
-         zstab = 0.5 + SIGN(0.5_wp, zta)
-         !
-         psi_h_coare(ji,jj) = (1. - zstab) * ( (1. - zf)*zpsi_k + zf*zpsi_c ) &
-            &                -   zstab     * ( (ABS(1. + 2.*zta/3.))**1.5     &
-            &                           + .6667*(zta - 14.28)/EXP(zc) + 8.525 )
-         !
+      !
+      zta = pzeta(ji,jj)
+      !
+      zphi_h = (ABS(1. - 15.*zta))**.5  !! Kansas unstable   (zphi_h = zphi_m**2 when unstable, zphi_m when stable)
+      !
+      zpsi_k = 2.*LOG((1. + zphi_h)/2.)
+      !
+      zphi_c = (ABS(1. - 34.15*zta))**.3333   !! Convective
+      !
+      zpsi_c = 1.5*LOG((1. + zphi_c + zphi_c*zphi_c)/3.) &
+         &    -1.7320508*ATAN((1. + 2.*zphi_c)/1.7320508) + 1.813799447
+      !
+      zf = zta*zta
+      zf = zf/(1. + zf)
+      zc = MIN(50._wp,0.35_wp*zta)
+      zstab = 0.5 + SIGN(0.5_wp, zta)
+      !
+      psi_h_coare(ji,jj) = (1. - zstab) * ( (1. - zf)*zpsi_k + zf*zpsi_c ) &
+         &                -   zstab     * ( (ABS(1. + 2.*zta/3.))**1.5     &
+         &                           + .6667*(zta - 14.28)/EXP(zc) + 8.525 )
+      !
       END_2D
       !
Index: /NEMO/trunk/src/OCE/SBC/sbcblk_algo_ecmwf.F90
===================================================================
--- /NEMO/trunk/src/OCE/SBC/sbcblk_algo_ecmwf.F90	(revision 12614)
+++ /NEMO/trunk/src/OCE/SBC/sbcblk_algo_ecmwf.F90	(revision 12615)
@@ -98,5 +98,5 @@
       &                      Qsw, rad_lw, slp, pdT_cs,                                & ! optionals for cool-skin (and warm-layer)
       &                      pdT_wl, pHz_wl )                                           ! optionals for warm-layer only
-      !!----------------------------------------------------------------------
+      !!----------------------------------------------------------------------------------
       !!                      ***  ROUTINE  turb_ecmwf  ***
       !!
@@ -184,7 +184,7 @@
       LOGICAL :: l_zt_equal_zu = .FALSE.      ! if q and t are given at same height as U
       !
-      REAL(wp), DIMENSION(jpi,jpj) ::  u_star, t_star, q_star
-      REAL(wp), DIMENSION(jpi,jpj) :: dt_zu, dq_zu     
-      REAL(wp), DIMENSION(jpi,jpj) :: znu_a !: Nu_air, Viscosity of air
+      REAL(wp), DIMENSION(jpi,jpj) :: u_star, t_star, q_star
+      REAL(wp), DIMENSION(jpi,jpj) :: dt_zu, dq_zu
+      REAL(wp), DIMENSION(jpi,jpj) :: znu_a         !: Nu_air, Viscosity of air
       REAL(wp), DIMENSION(jpi,jpj) :: Linv  !: 1/L (inverse of Monin Obukhov length...
       REAL(wp), DIMENSION(jpi,jpj) :: z0, z0t, z0q
@@ -196,9 +196,7 @@
       CHARACTER(len=40), PARAMETER :: crtnm = 'turb_ecmwf@sbcblk_algo_ecmwf.F90'
       !!----------------------------------------------------------------------------------
-
       IF( kt == nit000 ) CALL SBCBLK_ALGO_ECMWF_INIT(l_use_cs, l_use_wl)
 
-      l_zt_equal_zu = .FALSE.
-      IF( ABS(zu - zt) < 0.01_wp )   l_zt_equal_zu = .TRUE.    ! testing "zu == zt" is risky with double precision
+      l_zt_equal_zu = ( ABS(zu - zt) < 0.01_wp ) ! testing "zu == zt" is risky with double precision
 
       !! Initializations for cool skin and warm layer:
@@ -413,26 +411,26 @@
       !!----------------------------------------------------------------------------------
       DO_2D_11_11
-         !
-         zzeta = MIN( pzeta(ji,jj) , 5._wp ) !! Very stable conditions (L positif and big!):
-         !
-         ! Unstable (Paulson 1970):
-         !   eq.3.20, Chap.3, p.33, IFS doc - Cy31r1
-         zx = SQRT(ABS(1._wp - 16._wp*zzeta))
-         ztmp = 1._wp + SQRT(zx)
-         ztmp = ztmp*ztmp
-         psi_unst = LOG( 0.125_wp*ztmp*(1._wp + zx) )   &
-            &       -2._wp*ATAN( SQRT(zx) ) + 0.5_wp*rpi
-         !
-         ! Unstable:
-         ! eq.3.22, Chap.3, p.33, IFS doc - Cy31r1
-         psi_stab = -2._wp/3._wp*(zzeta - 5._wp/0.35_wp)*EXP(-0.35_wp*zzeta) &
-            &       - zzeta - 2._wp/3._wp*5._wp/0.35_wp
-         !
-         ! Combining:
-         stab = 0.5_wp + SIGN(0.5_wp, zzeta) ! zzeta > 0 => stab = 1
-         !
-         psi_m_ecmwf(ji,jj) = (1._wp - stab) * psi_unst & ! (zzeta < 0) Unstable
-            &                +      stab  * psi_stab      ! (zzeta > 0) Stable
-         !
+      !
+      zzeta = MIN( pzeta(ji,jj) , 5._wp ) !! Very stable conditions (L positif and big!):
+      !
+      ! Unstable (Paulson 1970):
+      !   eq.3.20, Chap.3, p.33, IFS doc - Cy31r1
+      zx = SQRT(ABS(1._wp - 16._wp*zzeta))
+      ztmp = 1._wp + SQRT(zx)
+      ztmp = ztmp*ztmp
+      psi_unst = LOG( 0.125_wp*ztmp*(1._wp + zx) )   &
+         &       -2._wp*ATAN( SQRT(zx) ) + 0.5_wp*rpi
+      !
+      ! Unstable:
+      ! eq.3.22, Chap.3, p.33, IFS doc - Cy31r1
+      psi_stab = -2._wp/3._wp*(zzeta - 5._wp/0.35_wp)*EXP(-0.35_wp*zzeta) &
+         &       - zzeta - 2._wp/3._wp*5._wp/0.35_wp
+      !
+      ! Combining:
+      stab = 0.5_wp + SIGN(0.5_wp, zzeta) ! zzeta > 0 => stab = 1
+      !
+      psi_m_ecmwf(ji,jj) = (1._wp - stab) * psi_unst & ! (zzeta < 0) Unstable
+         &                +      stab  * psi_stab      ! (zzeta > 0) Stable
+      !
       END_2D
    END FUNCTION psi_m_ecmwf
@@ -458,23 +456,23 @@
       !
       DO_2D_11_11
-         !
-         zzeta = MIN(pzeta(ji,jj) , 5._wp)   ! Very stable conditions (L positif and big!):
-         !
-         zx  = ABS(1._wp - 16._wp*zzeta)**.25        ! this is actually (1/phi_m)**2  !!!
-         !                                     ! eq.3.19, Chap.3, p.33, IFS doc - Cy31r1
-         ! Unstable (Paulson 1970) :
-         psi_unst = 2._wp*LOG(0.5_wp*(1._wp + zx*zx))   ! eq.3.20, Chap.3, p.33, IFS doc - Cy31r1
-         !
-         ! Stable:
-         psi_stab = -2._wp/3._wp*(zzeta - 5._wp/0.35_wp)*EXP(-0.35_wp*zzeta) & ! eq.3.22, Chap.3, p.33, IFS doc - Cy31r1
-            &       - ABS(1._wp + 2._wp/3._wp*zzeta)**1.5_wp - 2._wp/3._wp*5._wp/0.35_wp + 1._wp
-         ! LB: added ABS() to avoid NaN values when unstable, which contaminates the unstable solution...
-         !
-         stab = 0.5_wp + SIGN(0.5_wp, zzeta) ! zzeta > 0 => stab = 1
-         !
-         !
-         psi_h_ecmwf(ji,jj) = (1._wp - stab) * psi_unst &   ! (zzeta < 0) Unstable
-            &                +    stab    * psi_stab        ! (zzeta > 0) Stable
-         !
+      !
+      zzeta = MIN(pzeta(ji,jj) , 5._wp)   ! Very stable conditions (L positif and big!):
+      !
+      zx  = ABS(1._wp - 16._wp*zzeta)**.25        ! this is actually (1/phi_m)**2  !!!
+      !                                     ! eq.3.19, Chap.3, p.33, IFS doc - Cy31r1
+      ! Unstable (Paulson 1970) :
+      psi_unst = 2._wp*LOG(0.5_wp*(1._wp + zx*zx))   ! eq.3.20, Chap.3, p.33, IFS doc - Cy31r1
+      !
+      ! Stable:
+      psi_stab = -2._wp/3._wp*(zzeta - 5._wp/0.35_wp)*EXP(-0.35_wp*zzeta) & ! eq.3.22, Chap.3, p.33, IFS doc - Cy31r1
+         &       - ABS(1._wp + 2._wp/3._wp*zzeta)**1.5_wp - 2._wp/3._wp*5._wp/0.35_wp + 1._wp
+      ! LB: added ABS() to avoid NaN values when unstable, which contaminates the unstable solution...
+      !
+      stab = 0.5_wp + SIGN(0.5_wp, zzeta) ! zzeta > 0 => stab = 1
+      !
+      !
+      psi_h_ecmwf(ji,jj) = (1._wp - stab) * psi_unst &   ! (zzeta < 0) Unstable
+         &                +    stab    * psi_stab        ! (zzeta > 0) Stable
+      !
       END_2D
    END FUNCTION psi_h_ecmwf
Index: /NEMO/trunk/src/OCE/SBC/sbcblk_algo_ncar.F90
===================================================================
--- /NEMO/trunk/src/OCE/SBC/sbcblk_algo_ncar.F90	(revision 12614)
+++ /NEMO/trunk/src/OCE/SBC/sbcblk_algo_ncar.F90	(revision 12615)
@@ -112,7 +112,5 @@
       REAL(wp), DIMENSION(jpi,jpj) ::   stab          ! stability test integer
       !!----------------------------------------------------------------------------------
-      !
-      l_zt_equal_zu = .FALSE.
-      IF( ABS(zu - zt) < 0.01_wp )   l_zt_equal_zu = .TRUE.    ! testing "zu == zt" is risky with double precision
+      l_zt_equal_zu = ( ABS(zu - zt) < 0.01_wp ) ! testing "zu == zt" is risky with double precision
 
       U_blk = MAX( 0.5_wp , U_zu )   !  relative wind speed at zu (normally 10m), we don't want to fall under 0.5 m/s
@@ -143,6 +141,7 @@
       ENDIF
 
-      !! Initializing values at z_u with z_t values:
-      t_zu = t_zt   ;   q_zu = q_zt
+      !! First guess of temperature and humidity at height zu:
+      t_zu = MAX( t_zt ,  180._wp )   ! who knows what's given on masked-continental regions...
+      q_zu = MAX( q_zt , 1.e-6_wp )   !               "
 
       !! ITERATION BLOCK
Index: /NEMO/trunk/src/OCE/SBC/sbcblk_phy.F90
===================================================================
--- /NEMO/trunk/src/OCE/SBC/sbcblk_phy.F90	(revision 12614)
+++ /NEMO/trunk/src/OCE/SBC/sbcblk_phy.F90	(revision 12615)
@@ -520,9 +520,9 @@
          zCe = zz0*pqst(ji,jj)/zdq
 
-         CALL BULK_FORMULA( pzu, pTs(ji,jj), pqs(ji,jj), pTa(ji,jj), pqa(ji,jj), &
-            &              zCd, zCh, zCe,                                        &
-            &              pwnd(ji,jj), pUb(ji,jj), pslp(ji,jj),                 &
-            &              pTau(ji,jj), zQsen, zQlat )
-
+         CALL BULK_FORMULA_SCLR( pzu, pTs(ji,jj), pqs(ji,jj), pTa(ji,jj), pqa(ji,jj), &
+            &                    zCd, zCh, zCe,                                       &
+            &                    pwnd(ji,jj), pUb(ji,jj), pslp(ji,jj),                &
+            &                    pTau(ji,jj), zQsen, zQlat )
+         
          zTs2  = pTs(ji,jj)*pTs(ji,jj)
          zQlw  = emiss_w*(prlw(ji,jj) - stefan*zTs2*zTs2) ! Net longwave flux
@@ -535,8 +535,63 @@
 
 
-   SUBROUTINE BULK_FORMULA_VCTR( pzu, pTs, pqs, pTa, pqa,  &
-      &                          pCd, pCh, pCe,            &
-      &                          pwnd, pUb, pslp,          &
-      &                          pTau, pQsen, pQlat,  pEvap, prhoa )
+   SUBROUTINE BULK_FORMULA_SCLR( pzu, pTs, pqs, pTa, pqa, &
+      &                          pCd, pCh, pCe,           &
+      &                          pwnd, pUb, pslp,         &
+      &                          pTau, pQsen, pQlat,      &
+      &                          pEvap, prhoa, pfact_evap )
+      !!----------------------------------------------------------------------------------
+      REAL(wp),                     INTENT(in)  :: pzu  ! height above the sea-level where all this takes place (normally 10m)
+      REAL(wp), INTENT(in)  :: pTs  ! water temperature at the air-sea interface [K]
+      REAL(wp), INTENT(in)  :: pqs  ! satur. spec. hum. at T=pTs   [kg/kg]
+      REAL(wp), INTENT(in)  :: pTa  ! potential air temperature at z=pzu [K]
+      REAL(wp), INTENT(in)  :: pqa  ! specific humidity at z=pzu [kg/kg]
+      REAL(wp), INTENT(in)  :: pCd
+      REAL(wp), INTENT(in)  :: pCh
+      REAL(wp), INTENT(in)  :: pCe
+      REAL(wp), INTENT(in)  :: pwnd ! wind speed module at z=pzu [m/s]
+      REAL(wp), INTENT(in)  :: pUb  ! bulk wind speed at z=pzu (inc. pot. effect of gustiness etc) [m/s]
+      REAL(wp), INTENT(in)  :: pslp ! sea-level atmospheric pressure [Pa]
+      !!
+      REAL(wp), INTENT(out) :: pTau  ! module of the wind stress [N/m^2]
+      REAL(wp), INTENT(out) :: pQsen !  [W/m^2]
+      REAL(wp), INTENT(out) :: pQlat !  [W/m^2]
+      !!
+      REAL(wp), INTENT(out), OPTIONAL :: pEvap ! Evaporation [kg/m^2/s]
+      REAL(wp), INTENT(out), OPTIONAL :: prhoa ! Air density at z=pzu [kg/m^3]
+      REAL(wp), INTENT(in) , OPTIONAL :: pfact_evap  ! ABOMINATION: corrective factor for evaporation (doing this against my will! /laurent)
+      !!
+      REAL(wp) :: ztaa, zgamma, zrho, zUrho, zevap, zfact_evap
+      INTEGER  :: jq
+      !!----------------------------------------------------------------------------------
+      zfact_evap = 1._wp
+      IF( PRESENT(pfact_evap) ) zfact_evap = pfact_evap
+      
+      !! Need ztaa, absolute temperature at pzu (formula to estimate rho_air needs absolute temperature, not the potential temperature "pTa")
+      ztaa = pTa ! first guess...
+      DO jq = 1, 4
+         zgamma = gamma_moist( 0.5*(ztaa+pTs) , pqa )  !LOLO: why not "0.5*(pqs+pqa)" rather then "pqa" ???
+         ztaa = pTa - zgamma*pzu   ! Absolute temp. is slightly colder...
+      END DO
+      zrho = rho_air(ztaa, pqa, pslp)
+      zrho = rho_air(ztaa, pqa, pslp-zrho*grav*pzu) ! taking into account that we are pzu m above the sea level where SLP is given!
+
+      zUrho = pUb*MAX(zrho, 1._wp)     ! rho*U10
+
+      pTau = zUrho * pCd * pwnd ! Wind stress module
+
+      zevap = zUrho * pCe * (pqa - pqs)
+      pQsen = zUrho * pCh * (pTa - pTs) * cp_air(pqa)
+      pQlat = L_vap(pTs) * zevap
+
+      IF( PRESENT(pEvap) ) pEvap = - zfact_evap * zevap
+      IF( PRESENT(prhoa) ) prhoa = zrho
+
+   END SUBROUTINE BULK_FORMULA_SCLR
+
+   SUBROUTINE BULK_FORMULA_VCTR( pzu, pTs, pqs, pTa, pqa, &
+      &                          pCd, pCh, pCe,           &
+      &                          pwnd, pUb, pslp,         &
+      &                          pTau, pQsen, pQlat,      & 
+      &                          pEvap, prhoa, pfact_evap )      
       !!----------------------------------------------------------------------------------
       REAL(wp),                     INTENT(in)  :: pzu  ! height above the sea-level where all this takes place (normally 10m)
@@ -558,85 +613,25 @@
       REAL(wp), DIMENSION(jpi,jpj), INTENT(out), OPTIONAL :: pEvap ! Evaporation [kg/m^2/s]
       REAL(wp), DIMENSION(jpi,jpj), INTENT(out), OPTIONAL :: prhoa ! Air density at z=pzu [kg/m^3]
-      !!
-      REAL(wp) :: ztaa, zgamma, zrho, zUrho, zevap
-      INTEGER  :: ji, jj, jq     ! dummy loop indices
-      !!----------------------------------------------------------------------------------
-      DO_2D_11_11
-
-         !! Need ztaa, absolute temperature at pzu (formula to estimate rho_air needs absolute temperature, not the potential temperature "pTa")
-         ztaa = pTa(ji,jj) ! first guess...
-         DO jq = 1, 4
-            zgamma = gamma_moist( 0.5*(ztaa+pTs(ji,jj)) , pqa(ji,jj) )
-            ztaa = pTa(ji,jj) - zgamma*pzu   ! Absolute temp. is slightly colder...
-         END DO
-         zrho = rho_air(ztaa, pqa(ji,jj), pslp(ji,jj))
-         zrho = rho_air(ztaa, pqa(ji,jj), pslp(ji,jj)-zrho*grav*pzu) ! taking into account that we are pzu m above the sea level where SLP is given!
-
-         zUrho = pUb(ji,jj)*MAX(zrho, 1._wp)     ! rho*U10
-
-         pTau(ji,jj) = zUrho * pCd(ji,jj) * pwnd(ji,jj) ! Wind stress module
-
-         zevap        = zUrho * pCe(ji,jj) * (pqa(ji,jj) - pqs(ji,jj))
-         pQsen(ji,jj) = zUrho * pCh(ji,jj) * (pTa(ji,jj) - pTs(ji,jj)) * cp_air(pqa(ji,jj))
-         pQlat(ji,jj) = L_vap(pTs(ji,jj)) * zevap
-
-         IF( PRESENT(pEvap) ) pEvap(ji,jj) = - zevap
+      REAL(wp),                     INTENT(in) , OPTIONAL :: pfact_evap  ! ABOMINATION: corrective factor for evaporation (doing this against my will! /laurent)
+      !!
+      REAL(wp) :: ztaa, zgamma, zrho, zUrho, zevap, zfact_evap
+      INTEGER  :: ji, jj
+      !!----------------------------------------------------------------------------------
+      zfact_evap = 1._wp
+      IF( PRESENT(pfact_evap) ) zfact_evap = pfact_evap
+
+      DO_2D_11_11
+
+         CALL BULK_FORMULA_SCLR( pzu, pTs(ji,jj), pqs(ji,jj), pTa(ji,jj), pqa(ji,jj), &
+            &                    pCd(ji,jj), pCh(ji,jj), pCe(ji,jj),                  &
+            &                    pwnd(ji,jj), pUb(ji,jj), pslp(ji,jj),                &
+            &                    pTau(ji,jj), pQsen(ji,jj), pQlat(ji,jj),             &
+            &                    pEvap=zevap, prhoa=zrho, pfact_evap=zfact_evap       )
+
+         IF( PRESENT(pEvap) ) pEvap(ji,jj) = zevap
          IF( PRESENT(prhoa) ) prhoa(ji,jj) = zrho
-
+   
       END_2D
    END SUBROUTINE BULK_FORMULA_VCTR
-
-
-   SUBROUTINE BULK_FORMULA_SCLR( pzu, pTs, pqs, pTa, pqa, &
-      &                          pCd, pCh, pCe,           &
-      &                          pwnd, pUb, pslp,         &
-      &                          pTau, pQsen, pQlat,  pEvap, prhoa )
-      !!----------------------------------------------------------------------------------
-      REAL(wp),                     INTENT(in)  :: pzu  ! height above the sea-level where all this takes place (normally 10m)
-      REAL(wp), INTENT(in)  :: pTs  ! water temperature at the air-sea interface [K]
-      REAL(wp), INTENT(in)  :: pqs  ! satur. spec. hum. at T=pTs   [kg/kg]
-      REAL(wp), INTENT(in)  :: pTa  ! potential air temperature at z=pzu [K]
-      REAL(wp), INTENT(in)  :: pqa  ! specific humidity at z=pzu [kg/kg]
-      REAL(wp), INTENT(in)  :: pCd
-      REAL(wp), INTENT(in)  :: pCh
-      REAL(wp), INTENT(in)  :: pCe
-      REAL(wp), INTENT(in)  :: pwnd ! wind speed module at z=pzu [m/s]
-      REAL(wp), INTENT(in)  :: pUb  ! bulk wind speed at z=pzu (inc. pot. effect of gustiness etc) [m/s]
-      REAL(wp), INTENT(in)  :: pslp ! sea-level atmospheric pressure [Pa]
-      !!
-      REAL(wp), INTENT(out) :: pTau  ! module of the wind stress [N/m^2]
-      REAL(wp), INTENT(out) :: pQsen !  [W/m^2]
-      REAL(wp), INTENT(out) :: pQlat !  [W/m^2]
-      !!
-      REAL(wp), INTENT(out), OPTIONAL :: pEvap ! Evaporation [kg/m^2/s]
-      REAL(wp), INTENT(out), OPTIONAL :: prhoa ! Air density at z=pzu [kg/m^3]
-      !!
-      REAL(wp) :: ztaa, zgamma, zrho, zUrho, zevap
-      INTEGER  :: jq
-      !!----------------------------------------------------------------------------------
-
-      !! Need ztaa, absolute temperature at pzu (formula to estimate rho_air needs absolute temperature, not the potential temperature "pTa")
-      ztaa = pTa ! first guess...
-      DO jq = 1, 4
-         zgamma = gamma_moist( 0.5*(ztaa+pTs) , pqa )
-         ztaa = pTa - zgamma*pzu   ! Absolute temp. is slightly colder...
-      END DO
-      zrho = rho_air(ztaa, pqa, pslp)
-      zrho = rho_air(ztaa, pqa, pslp-zrho*grav*pzu) ! taking into account that we are pzu m above the sea level where SLP is given!
-
-      zUrho = pUb*MAX(zrho, 1._wp)     ! rho*U10
-
-      pTau = zUrho * pCd * pwnd ! Wind stress module
-
-      zevap = zUrho * pCe * (pqa - pqs)
-      pQsen = zUrho * pCh * (pTa - pTs) * cp_air(pqa)
-      pQlat = L_vap(pTs) * zevap
-
-      IF( PRESENT(pEvap) ) pEvap = - zevap
-      IF( PRESENT(prhoa) ) prhoa = zrho
-
-   END SUBROUTINE BULK_FORMULA_SCLR
-
-
 
 
Index: /NEMO/trunk/src/SAS/sbcssm.F90
===================================================================
--- /NEMO/trunk/src/SAS/sbcssm.F90	(revision 12614)
+++ /NEMO/trunk/src/SAS/sbcssm.F90	(revision 12615)
@@ -26,5 +26,5 @@
    USE lib_mpp        ! distributed memory computing library
    USE prtctl         ! print control
-   USE fldread        ! read input fields 
+   USE fldread        ! read input fields
    USE timing         ! Timing
 
@@ -38,5 +38,5 @@
    LOGICAL            ::   ln_3d_uve     ! specify whether input velocity data is 3D
    LOGICAL            ::   ln_read_frq   ! specify whether we must read frq or not
-   
+
    LOGICAL            ::   l_sasread     ! Ice intilisation: =T read a file ; =F anaytical initilaistion
    LOGICAL            ::   l_initdone = .false.
@@ -69,10 +69,10 @@
       !!               for an off-line simulation using surface processes only
       !!
-      !! ** Method : calculates the position of data 
+      !! ** Method : calculates the position of data
       !!             - interpolates data if needed
       !!----------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt   ! ocean time-step index
       INTEGER, INTENT(in) ::   Kbb, Kmm   ! ocean time level indices
-                          ! (not needed for SAS but needed to keep a consistent interface in sbcmod.F90)
+      ! (not needed for SAS but needed to keep a consistent interface in sbcmod.F90)
       !
       INTEGER  ::   ji, jj     ! dummy loop indices
@@ -82,25 +82,25 @@
       !
       IF( ln_timing )   CALL timing_start( 'sbc_ssm')
-     
+
       IF ( l_sasread ) THEN
          IF( nfld_3d > 0 ) CALL fld_read( kt, 1, sf_ssm_3d )      !==   read data at kt time step   ==!
          IF( nfld_2d > 0 ) CALL fld_read( kt, 1, sf_ssm_2d )      !==   read data at kt time step   ==!
-         ! 
+         !
          IF( ln_3d_uve ) THEN
             IF( .NOT. ln_linssh ) THEN
-               e3t_m(:,:) = sf_ssm_3d(jf_e3t)%fnow(:,:,1) * tmask(:,:,1) ! vertical scale factor 
+               e3t_m(:,:) = sf_ssm_3d(jf_e3t)%fnow(:,:,1) * tmask(:,:,1) ! vertical scale factor
             ELSE
                e3t_m(:,:) = e3t_0(:,:,1)                                 ! vertical scale factor
             ENDIF
             ssu_m(:,:) = sf_ssm_3d(jf_usp)%fnow(:,:,1) * umask(:,:,1)    ! u-velocity
-            ssv_m(:,:) = sf_ssm_3d(jf_vsp)%fnow(:,:,1) * vmask(:,:,1)    ! v-velocity 
+            ssv_m(:,:) = sf_ssm_3d(jf_vsp)%fnow(:,:,1) * vmask(:,:,1)    ! v-velocity
          ELSE
             IF( .NOT. ln_linssh ) THEN
-               e3t_m(:,:) = sf_ssm_2d(jf_e3t)%fnow(:,:,1) * tmask(:,:,1) ! vertical scale factor 
+               e3t_m(:,:) = sf_ssm_2d(jf_e3t)%fnow(:,:,1) * tmask(:,:,1) ! vertical scale factor
             ELSE
                e3t_m(:,:) = e3t_0(:,:,1)                                 ! vertical scale factor
             ENDIF
             ssu_m(:,:) = sf_ssm_2d(jf_usp)%fnow(:,:,1) * umask(:,:,1)    ! u-velocity
-            ssv_m(:,:) = sf_ssm_2d(jf_vsp)%fnow(:,:,1) * vmask(:,:,1)    ! v-velocity 
+            ssv_m(:,:) = sf_ssm_2d(jf_vsp)%fnow(:,:,1) * vmask(:,:,1)    ! v-velocity
          ENDIF
          !
@@ -123,5 +123,5 @@
          ssh  (:,:,Kmm) = 0._wp                              !              - -
       ENDIF
-      
+
       IF ( nn_ice == 1 ) THEN
          ts(:,:,1,jp_tem,Kmm) = sst_m(:,:)
@@ -132,5 +132,5 @@
       uu (:,:,1,Kbb) = ssu_m(:,:)
       vv (:,:,1,Kbb) = ssv_m(:,:)
- 
+
       IF(sn_cfctl%l_prtctl) THEN            ! print control
          CALL prt_ctl(tab2d_1=sst_m, clinfo1=' sst_m   - : ', mask1=tmask   )
@@ -162,8 +162,8 @@
       !!                  ***  ROUTINE sbc_ssm_init  ***
       !!
-      !! ** Purpose :   Initialisation of sea surface mean data     
-      !!----------------------------------------------------------------------
-      INTEGER, INTENT(in) ::   Kbb, Kmm   ! ocean time level indices 
-                          ! (not needed for SAS but needed to keep a consistent interface in sbcmod.F90)
+      !! ** Purpose :   Initialisation of sea surface mean data
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT(in) ::   Kbb, Kmm   ! ocean time level indices
+      ! (not needed for SAS but needed to keep a consistent interface in sbcmod.F90)
       INTEGER  :: ierr, ierr0, ierr1, ierr2, ierr3   ! return error code
       INTEGER  :: ifpr                               ! dummy loop indice
@@ -195,8 +195,8 @@
 902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namsbc_sas in configuration namelist' )
       IF(lwm) WRITE ( numond, namsbc_sas )
-      !           
+      !
       IF(lwp) THEN                              ! Control print
          WRITE(numout,*) '   Namelist namsbc_sas'
-         WRITE(numout,*) '      Initialisation using an input file                                 l_sasread   = ', l_sasread 
+         WRITE(numout,*) '      Initialisation using an input file                                 l_sasread   = ', l_sasread
          WRITE(numout,*) '      Are we supplying a 3D u,v and e3 field                             ln_3d_uve   = ', ln_3d_uve
          WRITE(numout,*) '      Are we reading frq (fraction of qsr absorbed in the 1st T level)   ln_read_frq = ', ln_read_frq
@@ -226,8 +226,8 @@
          ln_closea = .false.
       ENDIF
-      
-      !                  
+
+      !
       IF( l_sasread ) THEN                       ! store namelist information in an array
-         ! 
+         !
          !! following code is a bit messy, but distinguishes between when u,v are 3d arrays and
          !! when we have other 3d arrays that we need to read in
@@ -275,5 +275,5 @@
          ENDIF
          !
-         ierr1 = 0    ! default definition if slf_?d(ifpr)%ln_tint = .false. 
+         ierr1 = 0    ! default definition if slf_?d(ifpr)%ln_tint = .false.
          IF( nfld_3d > 0 ) THEN
             ALLOCATE( sf_ssm_3d(nfld_3d), STAT=ierr )         ! set sf structure
@@ -282,5 +282,5 @@
             ENDIF
             DO ifpr = 1, nfld_3d
-                                            ALLOCATE( sf_ssm_3d(ifpr)%fnow(jpi,jpj,jpk)    , STAT=ierr0 )
+               ALLOCATE( sf_ssm_3d(ifpr)%fnow(jpi,jpj,jpk)    , STAT=ierr0 )
                IF( slf_3d(ifpr)%ln_tint )   ALLOCATE( sf_ssm_3d(ifpr)%fdta(jpi,jpj,jpk,2)  , STAT=ierr1 )
                IF( ierr0 + ierr1 > 0 ) THEN
@@ -298,5 +298,5 @@
             ENDIF
             DO ifpr = 1, nfld_2d
-                                            ALLOCATE( sf_ssm_2d(ifpr)%fnow(jpi,jpj,1)    , STAT=ierr0 )
+               ALLOCATE( sf_ssm_2d(ifpr)%fnow(jpi,jpj,1)    , STAT=ierr0 )
                IF( slf_2d(ifpr)%ln_tint )   ALLOCATE( sf_ssm_2d(ifpr)%fdta(jpi,jpj,1,2)  , STAT=ierr1 )
                IF( ierr0 + ierr1 > 0 ) THEN
Index: /NEMO/trunk/tests/STATION_ASF/EXPREF/launch_sasf.sh
===================================================================
--- /NEMO/trunk/tests/STATION_ASF/EXPREF/launch_sasf.sh	(revision 12614)
+++ /NEMO/trunk/tests/STATION_ASF/EXPREF/launch_sasf.sh	(revision 12615)
@@ -2,5 +2,13 @@
 
 # NEMO directory where to fetch compiled STATION_ASF nemo.exe + setup:
-NEMO_DIR="${HOME}/NEMO/NEMOvdev_r11085_ASINTER-05_Brodeau_Advanced_Bulk"
+NEMO_DIR=`pwd | sed -e "s|/tests/STATION_ASF/EXPREF||g"`
+
+echo "Using NEMO_DIR=${NEMO_DIR}"
+
+# what directory inside "tests" actually contains the compiled test-case?
+TC_DIR="STATION_ASF2"
+
+# => so the executable to use is:
+NEMO_EXE="${NEMO_DIR}/tests/${TC_DIR}/BLD/bin/nemo.exe"
 
 # Directory where to run the simulation:
@@ -24,5 +32,5 @@
 mkdir -p ${WORK_DIR}
 
-NEMO_EXE="${NEMO_DIR}/tests/STATION_ASF/BLD/bin/nemo.exe"
+
 if [ ! -f ${NEMO_EXE} ]; then echo " Mhhh, no compiled nemo.exe found into ${NEMO_DIR}/tests/STATION_ASF/BLD/bin !"; exit; fi
 
@@ -40,5 +48,5 @@
 rsync -avP ${FORC_DIR}/Station_PAPA_50N-145W*.nc ${WORK_DIR}/
 
-for CASE in "ECMWF-noskin" "COARE3p6-noskin" "ECMWF" "COARE3p6" "NCAR"; do
+for CASE in "ECMWF" "COARE3p6" "NCAR" "ECMWF-noskin" "COARE3p6-noskin"; do
 
     echo ; echo
@@ -56,5 +64,5 @@
     echo
     echo "Launching NEMO !"
-    ./nemo.exe 1> out_nemo.out 2>err_nemo.err
+    ./nemo.exe 1>out_nemo.out 2>err_nemo.err
     echo "Done!"
     echo
Index: /NEMO/trunk/tests/STATION_ASF/EXPREF/namelist_coare3p6-noskin_cfg
===================================================================
--- /NEMO/trunk/tests/STATION_ASF/EXPREF/namelist_coare3p6-noskin_cfg	(revision 12614)
+++ /NEMO/trunk/tests/STATION_ASF/EXPREF/namelist_coare3p6-noskin_cfg	(revision 12615)
@@ -33,6 +33,6 @@
    nn_time0    =       0   !  initial time of day in hhmm
    nn_leapy    =       0   !  Leap year calendar (1) or not (0)
-   ln_rstart   =  .false.   !  start from rest (F) or from a restart file (T)
-      nn_euler    =    1      !  = 0 : start with forward time step if ln_rstart=T
+   ln_rstart   = .false.   !  start from rest (F) or from a restart file (T)
+      ln_1st_euler = .false.  !  =T force a start with forward time step (ln_rstart=T)
       nn_rstctl   = 2      !  restart control ==> activated only if ln_rstart=T
       !                          !    = 0 nn_date0 read in namelist ; nn_it000 : read in namelist
Index: /NEMO/trunk/tests/STATION_ASF/EXPREF/namelist_coare3p6_cfg
===================================================================
--- /NEMO/trunk/tests/STATION_ASF/EXPREF/namelist_coare3p6_cfg	(revision 12614)
+++ /NEMO/trunk/tests/STATION_ASF/EXPREF/namelist_coare3p6_cfg	(revision 12615)
@@ -33,6 +33,6 @@
    nn_time0    =       0   !  initial time of day in hhmm
    nn_leapy    =       0   !  Leap year calendar (1) or not (0)
-   ln_rstart   =  .false.   !  start from rest (F) or from a restart file (T)
-      nn_euler    =    1      !  = 0 : start with forward time step if ln_rstart=T
+   ln_rstart   = .false.   !  start from rest (F) or from a restart file (T)
+      ln_1st_euler = .false.  !  =T force a start with forward time step (ln_rstart=T)
       nn_rstctl   = 2      !  restart control ==> activated only if ln_rstart=T
       !                          !    = 0 nn_date0 read in namelist ; nn_it000 : read in namelist
Index: /NEMO/trunk/tests/STATION_ASF/EXPREF/namelist_ecmwf-noskin_cfg
===================================================================
--- /NEMO/trunk/tests/STATION_ASF/EXPREF/namelist_ecmwf-noskin_cfg	(revision 12614)
+++ /NEMO/trunk/tests/STATION_ASF/EXPREF/namelist_ecmwf-noskin_cfg	(revision 12615)
@@ -33,6 +33,6 @@
    nn_time0    =       0   !  initial time of day in hhmm
    nn_leapy    =       0   !  Leap year calendar (1) or not (0)
-   ln_rstart   =  .false.   !  start from rest (F) or from a restart file (T)
-      nn_euler    =    1      !  = 0 : start with forward time step if ln_rstart=T
+   ln_rstart   = .false.   !  start from rest (F) or from a restart file (T)
+      ln_1st_euler = .false.  !  =T force a start with forward time step (ln_rstart=T)
       nn_rstctl   = 2      !  restart control ==> activated only if ln_rstart=T
       !                          !    = 0 nn_date0 read in namelist ; nn_it000 : read in namelist
Index: /NEMO/trunk/tests/STATION_ASF/EXPREF/namelist_ecmwf_cfg
===================================================================
--- /NEMO/trunk/tests/STATION_ASF/EXPREF/namelist_ecmwf_cfg	(revision 12614)
+++ /NEMO/trunk/tests/STATION_ASF/EXPREF/namelist_ecmwf_cfg	(revision 12615)
@@ -33,6 +33,6 @@
    nn_time0    =       0   !  initial time of day in hhmm
    nn_leapy    =       0   !  Leap year calendar (1) or not (0)
-   ln_rstart   =  .false.   !  start from rest (F) or from a restart file (T)
-      nn_euler    =    1      !  = 0 : start with forward time step if ln_rstart=T
+   ln_rstart   = .false.   !  start from rest (F) or from a restart file (T)
+      ln_1st_euler = .false.  !  =T force a start with forward time step (ln_rstart=T)
       nn_rstctl   = 2      !  restart control ==> activated only if ln_rstart=T
       !                          !    = 0 nn_date0 read in namelist ; nn_it000 : read in namelist
Index: /NEMO/trunk/tests/STATION_ASF/EXPREF/namelist_ncar_cfg
===================================================================
--- /NEMO/trunk/tests/STATION_ASF/EXPREF/namelist_ncar_cfg	(revision 12614)
+++ /NEMO/trunk/tests/STATION_ASF/EXPREF/namelist_ncar_cfg	(revision 12615)
@@ -33,6 +33,6 @@
    nn_time0    =       0   !  initial time of day in hhmm
    nn_leapy    =       0   !  Leap year calendar (1) or not (0)
-   ln_rstart   =  .false.   !  start from rest (F) or from a restart file (T)
-      nn_euler    =    1      !  = 0 : start with forward time step if ln_rstart=T
+   ln_rstart   = .false.   !  start from rest (F) or from a restart file (T)
+      ln_1st_euler = .false.  !  =T force a start with forward time step (ln_rstart=T)
       nn_rstctl   = 2      !  restart control ==> activated only if ln_rstart=T
       !                          !    = 0 nn_date0 read in namelist ; nn_it000 : read in namelist
Index: /NEMO/trunk/tests/STATION_ASF/MY_SRC/diawri.F90
===================================================================
--- /NEMO/trunk/tests/STATION_ASF/MY_SRC/diawri.F90	(revision 12614)
+++ /NEMO/trunk/tests/STATION_ASF/MY_SRC/diawri.F90	(revision 12615)
@@ -35,4 +35,5 @@
    USE iom            !
    USE ioipsl         !
+
 #if defined key_si3
    USE ice
@@ -56,6 +57,6 @@
 
    !!----------------------------------------------------------------------
-   !! NEMO/SAS 4.0 , NEMO Consortium (2018)
-   !! $Id: diawri.F90 10425 2018-12-19 21:54:16Z smasson $
+   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
+   !! $Id: diawri.F90 12493 2020-03-02 07:56:31Z smasson $
    !! Software governed by the CeCILL license (see ./LICENSE)
    !!----------------------------------------------------------------------
@@ -114,11 +115,16 @@
       INTEGER, DIMENSION(2) :: ierr
       !!----------------------------------------------------------------------
-      ierr = 0
-      ALLOCATE( ndex_hT(jpi*jpj) , ndex_T(jpi*jpj*jpk) ,     &
-         &      ndex_hU(jpi*jpj) , ndex_U(jpi*jpj*jpk) ,     &
-         &      ndex_hV(jpi*jpj) , ndex_V(jpi*jpj*jpk) , STAT=ierr(1) )
-      !
-      dia_wri_alloc = MAXVAL(ierr)
-      CALL mpp_sum( 'diawri', dia_wri_alloc )
+      IF( nn_write == -1 ) THEN
+         dia_wri_alloc = 0
+      ELSE
+         ierr = 0
+         ALLOCATE( ndex_hT(jpi*jpj) , ndex_T(jpi*jpj*jpk) ,     &
+            &      ndex_hU(jpi*jpj) , ndex_U(jpi*jpj*jpk) ,     &
+            &      ndex_hV(jpi*jpj) , ndex_V(jpi*jpj*jpk) , STAT=ierr(1) )
+         !
+         dia_wri_alloc = MAXVAL(ierr)
+         CALL mpp_sum( 'diawri', dia_wri_alloc )
+         !
+      ENDIF
       !
    END FUNCTION dia_wri_alloc
@@ -374,5 +380,5 @@
       CALL iom_rstput( 0, 0, inum, 'vozocrtx', uu(:,:,:,Kmm)                )    ! now i-velocity
       CALL iom_rstput( 0, 0, inum, 'vomecrty', vv(:,:,:,Kmm)                )    ! now j-velocity
-         CALL iom_rstput( 0, 0, inum, 'vovecrtz', ww             )    ! now k-velocity
+      CALL iom_rstput( 0, 0, inum, 'vovecrtz', ww             )    ! now k-velocity
       CALL iom_rstput( 0, 0, inum, 'sowaflup', emp - rnf         )    ! freshwater budget
       CALL iom_rstput( 0, 0, inum, 'sohefldo', qsr + qns         )    ! total heat flux
Index: /NEMO/trunk/tests/STATION_ASF/MY_SRC/nemogcm.F90
===================================================================
--- /NEMO/trunk/tests/STATION_ASF/MY_SRC/nemogcm.F90	(revision 12614)
+++ /NEMO/trunk/tests/STATION_ASF/MY_SRC/nemogcm.F90	(revision 12615)
@@ -2,5 +2,5 @@
    !!======================================================================
    !!                       ***  MODULE nemogcm   ***
-   !! StandAlone Surface module : surface fluxes
+   !!                      STATION_ASF (SAS meets C1D)
    !!======================================================================
    !! History :  3.6  ! 2011-11  (S. Alderson, G. Madec) original code
@@ -19,16 +19,14 @@
    !!----------------------------------------------------------------------
    USE step_oce       ! module used in the ocean time stepping module (step.F90)
-   USE sbc_oce        ! surface boundary condition: ocean #LB: rm?
    USE phycst         ! physical constant                  (par_cst routine)
    USE domain         ! domain initialization   (dom_init & dom_cfg routines)
    USE closea         ! treatment of closed seas (for ln_closea)
    USE usrdef_nam     ! user defined configuration
+   USE istate         ! initial state setting          (istate_init routine)
    USE step, ONLY : Nbb, Nnn, Naa, Nrhs ! time level indices
    USE daymod         ! calendar
    USE restart        ! open  restart file
-   !LB:USE step           ! NEMO time-stepping                 (stp     routine)
    USE c1d            ! 1D configuration
    USE step_c1d       ! Time stepping loop for the 1D configuration
-   USE sbcssm         !
    !
    USE lib_mpp        ! distributed memory computing
@@ -49,5 +47,5 @@
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
-   !! $Id: nemogcm.F90 11536 2019-09-11 13:54:18Z smasson $
+   !! $Id: nemogcm.F90 12489 2020-02-28 15:55:11Z davestorkey $
    !! Software governed by the CeCILL license (see ./LICENSE)
    !!----------------------------------------------------------------------
@@ -84,4 +82,7 @@
       !                            !==   time stepping   ==!
       !                            !-----------------------!
+      !
+      !                                               !== set the model time-step  ==!
+      !
       istp = nit000
       !
@@ -106,8 +107,7 @@
       !
 #if defined key_iomput
-      CALL xios_finalize  ! end mpp communications with xios
+                                    CALL xios_finalize  ! end mpp communications with xios
 #else
-      IF( lk_mpp   ) THEN   ;   CALL mppstop      ! end mpp communications
-      ENDIF
+      IF( lk_mpp )                  CALL mppstop      ! end mpp communications
 #endif
       !
@@ -161,9 +161,9 @@
       IF( lwm )   CALL ctl_opn(     numout,        'ocean.output', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, -1, .FALSE. )
       ! open reference and configuration namelist files
-      CALL load_nml( numnam_ref,        'namelist_ref',                                           -1, lwm )
-      CALL load_nml( numnam_cfg,        'namelist_cfg',                                           -1, lwm )
+                  CALL load_nml( numnam_ref,        'namelist_ref',                                           -1, lwm )
+                  CALL load_nml( numnam_cfg,        'namelist_cfg',                                           -1, lwm )
       IF( lwm )   CALL ctl_opn(     numond, 'output.namelist.dyn', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, -1, .FALSE. )
       ! open /dev/null file to be able to supress output write easily
-      CALL ctl_opn(     numnul,           '/dev/null', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, -1, .FALSE. )
+                  CALL ctl_opn(     numnul,           '/dev/null', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, -1, .FALSE. )
       !
       !                             !--------------------!
@@ -235,5 +235,5 @@
 903   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namcfg in reference namelist' )
       READ  ( numnam_cfg, namcfg, IOSTAT = ios, ERR = 904 )
-904   IF( ios >  0 )   CALL ctl_nam ( ios , 'namcfg in configuration namelist' )
+904   IF( ios >  0 )   CALL ctl_nam ( ios , 'namcfg in configuration namelist' )   
       !
       IF( ln_read_cfg ) THEN            ! Read sizes in domain configuration file
@@ -266,26 +266,16 @@
       IF( ln_timing    )   CALL timing_start( 'nemo_init')
       !
-      CALL     phy_cst         ! Physical constants
-      CALL     eos_init        ! Equation of state
+                           CALL     phy_cst         ! Physical constants
+                           CALL     eos_init        ! Equation of state
       IF( lk_c1d       )   CALL     c1d_init        ! 1D column configuration
-      CALL     dom_init( Nbb, Nnn, Naa, "OPA") ! Domain
+                           CALL     dom_init( Nbb, Nnn, Naa, "OPA") ! Domain
       IF( sn_cfctl%l_prtctl )   &
          &                 CALL prt_ctl_init        ! Print control
-
-      IF( ln_rstart ) THEN                    ! Restart from a file                                                                                 
-         !                                    ! -------------------                                                                                 
-         CALL rst_read( Nbb, Nnn )            ! Read the restart file                                                                               
-         CALL day_init                        ! model calendar (using both namelist and restart infos)                                              
-         !                                                                                                                                          
-      ELSE                                    ! Start from rest                                                                                     
-         !                                    ! ---------------                                                                                     
-         numror = 0                           ! define numror = 0 -> no restart file to read                                                        
-         neuler = 0                           ! Set time-step indicator at nit000 (euler forward)                                                   
-         CALL day_init                        ! model calendar (using both namelist and restart infos)                                              
-      ENDIF
-      !
-
-      !                                      ! external forcing
-      CALL     sbc_init( Nbb, Nnn, Naa )    ! surface boundary conditions (including sea-ice)
+      !
+      
+                           CALL  istate_init( Nbb, Nnn, Naa )    ! ocean initial state (Dynamics and tracers)
+
+      !                                      ! external forcing 
+                           CALL     sbc_init( Nbb, Nnn, Naa )    ! surface boundary conditions (including sea-ice)
 
       !
@@ -321,8 +311,8 @@
          WRITE(numout,*) '                              sn_cfctl%l_prttrc  = ', sn_cfctl%l_prttrc
          WRITE(numout,*) '                              sn_cfctl%l_oasout  = ', sn_cfctl%l_oasout
-         WRITE(numout,*) '                              sn_cfctl%procmin   = ', sn_cfctl%procmin
-         WRITE(numout,*) '                              sn_cfctl%procmax   = ', sn_cfctl%procmax
-         WRITE(numout,*) '                              sn_cfctl%procincr  = ', sn_cfctl%procincr
-         WRITE(numout,*) '                              sn_cfctl%ptimincr  = ', sn_cfctl%ptimincr
+         WRITE(numout,*) '                              sn_cfctl%procmin   = ', sn_cfctl%procmin  
+         WRITE(numout,*) '                              sn_cfctl%procmax   = ', sn_cfctl%procmax  
+         WRITE(numout,*) '                              sn_cfctl%procincr  = ', sn_cfctl%procincr 
+         WRITE(numout,*) '                              sn_cfctl%ptimincr  = ', sn_cfctl%ptimincr 
          WRITE(numout,*) '      level of print                  nn_print   = ', nn_print
          WRITE(numout,*) '      Start i indice for SUM control  nn_ictls   = ', nn_ictls
@@ -439,5 +429,5 @@
       !!----------------------------------------------------------------------
       !
-      ierr =        oce_alloc    ()    ! ocean
+      ierr =        oce_alloc    ()    ! ocean 
       ierr = ierr + dia_wri_alloc()
       ierr = ierr + dom_oce_alloc()    ! ocean domain
@@ -448,5 +438,5 @@
    END SUBROUTINE nemo_alloc
 
-
+   
    SUBROUTINE nemo_set_cfctl(sn_cfctl, setto, for_all )
       !!----------------------------------------------------------------------
@@ -479,2 +469,3 @@
    !!======================================================================
 END MODULE nemogcm
+
Index: /NEMO/trunk/tests/STATION_ASF/MY_SRC/sbcssm.F90
===================================================================
--- /NEMO/trunk/tests/STATION_ASF/MY_SRC/sbcssm.F90	(revision 12614)
+++ /NEMO/trunk/tests/STATION_ASF/MY_SRC/sbcssm.F90	(revision 12615)
@@ -54,5 +54,5 @@
    !!----------------------------------------------------------------------
    !! NEMO/SAS 4.0 , NEMO Consortium (2018)
-   !! $Id: sbcssm.F90 10068 2018-08-28 14:09:04Z nicolasmartin $
+   !! $Id: sbcssm.F90 12377 2020-02-12 14:39:06Z acc $
    !! Software governed by the CeCILL license (see ./LICENSE)
    !!----------------------------------------------------------------------
Index: /NEMO/trunk/tests/STATION_ASF/MY_SRC/step_c1d.F90
===================================================================
--- /NEMO/trunk/tests/STATION_ASF/MY_SRC/step_c1d.F90	(revision 12614)
+++ /NEMO/trunk/tests/STATION_ASF/MY_SRC/step_c1d.F90	(revision 12615)
@@ -26,5 +26,5 @@
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
-   !! $Id: step_c1d.F90 10068 2018-08-28 14:09:04Z nicolasmartin $
+   !! $Id: step_c1d.F90 12377 2020-02-12 14:39:06Z acc $
    !! Software governed by the CeCILL license (see ./LICENSE)
    !!----------------------------------------------------------------------
