Index: branches/2013/dev_LOCEAN_CMCC_INGV_MERC_UKMO_2013/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90
===================================================================
--- branches/2013/dev_LOCEAN_CMCC_INGV_MERC_UKMO_2013/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90	(revision 4244)
+++ branches/2013/dev_LOCEAN_CMCC_INGV_MERC_UKMO_2013/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90	(revision 4245)
@@ -79,4 +79,7 @@
    REAL(wp) ::   rn_efac     ! multiplication factor for evaporation (clem)
    REAL(wp) ::   rn_vfac     ! multiplication factor for ice/ocean velocity in the calculation of wind stress (clem)
+   LOGICAL  ::   ln_bulk2z   ! logical flag for case where z(q,t) and z(u) are specified in the namelist
+   REAL(wp) ::   rn_zqt      ! z(q,t) : height of humidity and temperature measurements
+   REAL(wp) ::   rn_zu       ! z(u)   : height of wind measurements
 
    !! * Substitutions
@@ -134,5 +137,6 @@
       NAMELIST/namsbc_core/ cn_dir , ln_2m  , ln_taudif, rn_pfac, rn_efac, rn_vfac,  &
          &                  sn_wndi, sn_wndj, sn_humi  , sn_qsr ,           &
-         &                  sn_qlw , sn_tair, sn_prec  , sn_snow, sn_tdif
+         &                  sn_qlw , sn_tair, sn_prec  , sn_snow,           &
+         &                  sn_tdif, rn_zqt , ln_bulk2z, rn_zu
       !!---------------------------------------------------------------------
 
@@ -316,4 +320,18 @@
             &                      Cd    , Ch              , Ce  ,   &
             &                      zt_zu , zq_zu                   )
+      ELSE IF( ln_bulk2z ) THEN
+         !! If the height of the air temp./spec. hum. and wind are to be specified by hand :
+         IF( rn_zqt == rn_zu ) THEN
+            !! If air temp. and spec. hum. are at the same height as wind :
+            CALL TURB_CORE_1Z( rn_zu, zst   , sf(jp_tair)%fnow(:,:,1),       &
+               &                      zqsatw, sf(jp_humi)%fnow(:,:,1), wndm, &
+               &                      Cd    , Ch                     , Ce  )
+         ELSE
+            !! If air temp. and spec. hum. are at a different height to wind :
+            CALL TURB_CORE_2Z(rn_zqt, rn_zu , zst   , sf(jp_tair)%fnow,         &
+               &                              zqsatw, sf(jp_humi)%fnow, wndm,   &
+               &                              Cd    , Ch              , Ce  ,   &
+               &                              zt_zu , zq_zu                 )
+         ENDIF
       ELSE
          !! If air temp. and spec. hum. are given at same height than wind (10m) :
@@ -360,6 +378,6 @@
       !  Turbulent fluxes over ocean
       ! -----------------------------
-      IF( ln_2m ) THEN
-         ! Values of temp. and hum. adjusted to 10m must be used instead of 2m values
+      IF( ln_2m .OR. ( ln_bulk2z .AND. rn_zqt /= rn_zu ) ) THEN
+         ! Values of temp. and hum. adjusted to height of wind must be used
          zevap(:,:) = rn_efac * MAX( 0.e0, rhoa    *Ce(:,:)*( zqsatw(:,:) - zq_zu(:,:) ) * wndm(:,:) )   ! Evaporation
          zqsb (:,:) =                      rhoa*cpa*Ch(:,:)*( zst   (:,:) - zt_zu(:,:) ) * wndm(:,:)     ! Sensible Heat
@@ -789,6 +807,6 @@
       !!      Momentum, Latent and sensible heat exchange coefficients
       !!      Caution: this procedure should only be used in cases when air
-      !!      temperature (T_air) and air specific humidity (q_air) are at 2m
-      !!      whereas wind (dU) is at 10m.
+      !!      temperature (T_air) and air specific humidity (q_air) are at a
+      !!      different height to wind (dU).
       !!
       !! References :   Large & Yeager, 2004 : ???
@@ -808,5 +826,5 @@
 
       INTEGER :: j_itt
-      INTEGER , PARAMETER :: nb_itt = 3              ! number of itterations
+      INTEGER , PARAMETER :: nb_itt = 5              ! number of itterations
       REAL(wp), PARAMETER ::   grav   = 9.8          ! gravity                       
       REAL(wp), PARAMETER ::   kappa  = 0.4          ! von Karman's constant
@@ -905,5 +923,5 @@
            !!
            !! Shifting the neutral 10m transfer coefficients to (zu,zeta_u) :
-           xct = 1. + sqrt_Cd_n10/kappa*(log(zu/10.) - zpsi_m)
+           xct = 1. + sqrt_Cd_n10/kappa*(log(zu/10.) - zpsi_m)   ! L & Y eq. (10a)
            Cd = Cd_n10/(xct*xct) ; sqrt_Cd = sqrt(Cd)
          ENDIF
@@ -911,8 +929,8 @@
          xlogt = log(zu/10.) - zpsi_hu
          !!
-         xct = 1. + Ch_n10*xlogt/kappa/sqrt_Cd_n10
+         xct = 1. + Ch_n10*xlogt/kappa/sqrt_Cd_n10               ! L & Y eq. (10b)
          Ch  = Ch_n10*sqrt_Cd/sqrt_Cd_n10/xct
          !!
-         xct = 1. + Ce_n10*xlogt/kappa/sqrt_Cd_n10
+         xct = 1. + Ce_n10*xlogt/kappa/sqrt_Cd_n10               ! L & Y eq. (10c)
          Ce  = Ce_n10*sqrt_Cd/sqrt_Cd_n10/xct
          !!
