Index: branches/2011/dev_r2855_INGV2_3_blk_wave/DOC/TexFiles/Biblio/Biblio.bib
===================================================================
--- branches/2011/dev_r2855_INGV2_3_blk_wave/DOC/TexFiles/Biblio/Biblio.bib	(revision 2989)
+++ branches/2011/dev_r2855_INGV2_3_blk_wave/DOC/TexFiles/Biblio/Biblio.bib	(revision 2990)
@@ -355,4 +355,14 @@
 }
 
+@ARTICLE{Bignami_al_JGR95,
+  author = {F. Bignami and S. Marullo and R. Santoleri and M. E. Schiano},
+  title = {Longwave radiation budget in the Mediterranean Sea}, 
+  journal = JGR,
+  year = {1995},
+  volume = {100},  number = {C2},
+  pages = {2501--2514},
+  doi = {10.1029/94JC02496},
+}
+
 @ARTICLE{Blanke_al_JPO99,
   author = {B. Blanke and M. Arhan and G. Madec and S. Roche},
@@ -590,4 +600,13 @@
   doi = {10.1029/2002GL016473},
   url = {http://dx.doi.org/10.1029/2002GL016473}
+}
+
+@ARTICLE{Castellari_al_JMS1998,
+  author = {S> Castellari and  N. Pinardi and K. Leaman },
+  title = {A model study of air-sea interactions in the Mediterranean Sea.},
+  journal = JMS,
+  year = {1998},
+  volume = {18},
+  pages = {89--114}
 }
 
@@ -1247,4 +1266,13 @@
 }
 
+@ARTICLE{Hellerman_Rosenstein_JPO83,
+  author = {S. Hellerman and  M. Rosenstein },
+  title = {Normal monthly wind stress over the world ocean with error estimates},
+  journal = JPO,
+  year = {1983},
+  volume = {13},
+  pages = {1093--1104},
+}
+
 @ARTICLE{He_Ding_JSC01,
   author = {Y. He and C. H. Q. Ding},
@@ -1511,4 +1539,21 @@
   volume = {6},
   pages = {56--58}
+}
+@ARTICLE{Kondo1975,
+  author = {J. Kondo},
+  title = {Air-sea bulk transfer coefficients in diabatic conditions},
+  journal = {Boundary-Layer Meteorol},
+  year = {1975},
+  volume = {9},
+  pages = {91--112}
+}
+
+@ARTICLE{Lermusiaux2001,
+  author = {P. F. J. Lermusiaux},
+  title = {Evolving the subspace of three-dimensional miltiscale ocean variability: Massachusetts Bay},
+  journal = JMS,
+  year = {2001},
+  volume = {29},
+  pages = {385--422}
 }
 
@@ -1906,4 +1951,13 @@
   volume = {125},  number = {5},
   pages = {958--971}
+}
+
+@ARTICLE{Maggiore_al_PCE98,
+  author = {A. Maggiore and  M. Zavatarelli and M. G. Angelucci and N. Pinardi},
+  title = {Surface heat and water fluxes in the Adriatic Sea: seasonal and interannual variability},
+  journal = {Phys Chem Earth},
+  year = {1998},
+  volume = {23},
+  pages = {561--567}
 }
 
@@ -2115,4 +2169,13 @@
 }
 
+@ARTICLE{Oddo_al_OS09,
+  author = {P. Oddo and M. Adani and N. Pinardi and C. Fratianni and M. Tonani and D. Pettenuzzo},
+  title = {A nested Atlantic-Mediterranean Sea general circulation model for operational forecasting},
+  journal = OS,
+  year = {2009},
+  volume = {5},
+  pages = {1--13},
+}
+
 @PHDTHESIS{Olivier_PhD01,
   author = {F. Olivier},
@@ -2170,4 +2233,13 @@
 }
 
+@ARTICLE{Payne_JAS72,
+  author = {R. E. Payne},
+  title = {Albedo of the Sea Surface},
+  journal = JAS,
+  year = {1972},
+  volume = {29}
+  pages = {959--970}
+}
+
 @ARTICLE{Penduff_al_OM06,
   author = {T. Penduff and B. Barnier and J.-M. Molines and G. Madec},
@@ -2234,4 +2306,13 @@
   volume = {13},
   pages = {1154--1158}
+}
+
+@ARTICLE{Reed_JPO77,
+  author = {R. K. Reed},
+  title = {On estimating insolation over the ocean},
+  journal = JPO,
+  year = {1977},
+  volume = {1},
+  pages = {874--971}
 }
 
@@ -2502,4 +2583,13 @@
   volume = {8},
   pages = {175--201}
+}
+
+@ARTICLE{Tonani_al_OS08,
+  author = {M. Tonani and N. Pinardi and S. Dobricic and I. Pujol and C. Fratianni},
+  title = {A high-resolution free-surface model of the Mediterranean Sea},
+  journal = OS,
+  year = {2008},
+  volume = {4},
+  pages = {1--14}
 }
 
Index: branches/2011/dev_r2855_INGV2_3_blk_wave/DOC/TexFiles/Chapters/Chap_SBC.tex
===================================================================
--- branches/2011/dev_r2855_INGV2_3_blk_wave/DOC/TexFiles/Chapters/Chap_SBC.tex	(revision 2989)
+++ branches/2011/dev_r2855_INGV2_3_blk_wave/DOC/TexFiles/Chapters/Chap_SBC.tex	(revision 2990)
@@ -24,13 +24,13 @@
 \end{itemize}
 
-Four different ways to provide the first six fields to the ocean are available which 
+Five different ways to provide the first six fields to the ocean are available which 
 are controlled by namelist variables: an analytical formulation (\np{ln\_ana}~=~true), 
 a flux formulation (\np{ln\_flx}~=~true), a bulk formulae formulation (CORE 
-(\np{ln\_core}~=~true) or CLIO (\np{ln\_clio}~=~true) bulk formulae) and a coupled 
+(\np{ln\_core}~=~true), CLIO (\np{ln\_clio}~=~true) or MFS
+\footnote { Note that MFS bulk formulae compute fluxes only for the ocean component}
+(\np{ln\_ecmwf}~=~true) bulk formulae) and a coupled 
 formulation (exchanges with a atmospheric model via the OASIS coupler) 
 (\np{ln\_cpl}~=~true). When used, the atmospheric pressure forces both 
-ocean and ice dynamics (\np{ln\_apr\_dyn}~=~true)
-\footnote{The surface pressure field could be use in bulk formulae, nevertheless 
-none of the current bulk formulea (CLIO and CORE) uses the it.}. 
+ocean and ice dynamics (\np{ln\_apr\_dyn}~=~true).
 The frequency at which the six or seven fields have to be updated is the \np{nn\_fsbc} 
 namelist parameter. 
@@ -46,10 +46,11 @@
 (\np{nn\_ice}~=~0,1, 2 or 3); the addition of river runoffs as surface freshwater 
 fluxes or lateral inflow (\np{ln\_rnf}~=~true); the addition of a freshwater flux adjustment 
-in order to avoid a mean sea-level drift (\np{nn\_fwb}~=~0,~1~or~2); and the 
+in order to avoid a mean sea-level drift (\np{nn\_fwb}~=~0,~1~or~2); the 
 transformation of the solar radiation (if provided as daily mean) into a diurnal 
-cycle (\np{ln\_dm2dc}~=~true).
+cycle (\np{ln\_dm2dc}~=~true); and a neutral drag coefficient can be read from an external wave 
+model(\np{ln\_cdgw}~=~true). The latter option is possible only in case core or ecmwf bulk formulas are selected.
 
 In this chapter, we first discuss where the surface boundary condition appears in the
-model equations. Then we present the four ways of providing the surface boundary condition, 
+model equations. Then we present the five ways of providing the surface boundary condition, 
 followed by the description of the atmospheric pressure and the river runoff. 
 Next the scheme for interpolation on the fly is described.
@@ -480,6 +481,6 @@
 % Bulk formulation
 % ================================================================
-\section  [Bulk formulation (\textit{sbcblk\_core} or \textit{sbcblk\_clio}) ]
-		{Bulk formulation \small{(\mdl{sbcblk\_core} or \mdl{sbcblk\_clio} module)} }
+\section  [Bulk formulation (\textit{sbcblk\_core}, \textit{sbcblk\_clio} or \textit{sbcblk\_ecmwf}) ]
+		{Bulk formulation \small{(\mdl{sbcblk\_core} \mdl{sbcblk\_clio} \mdl{sbcblk\_ecmwf} modules)} }
 \label{SBC_blk}
 
@@ -487,10 +488,10 @@
 using bulk formulae and atmospheric fields and ocean (and ice) variables. 
 
-The atmospheric fields used depend on the bulk formulae used. Two bulk formulations 
-are available : the CORE and CLIO bulk formulea. The choice is made by setting to true
-one of the following namelist variable : \np{ln\_core} and \np{ln\_clio}.
-
-Note : in forced mode, when a sea-ice model is used, a bulk formulation have to be used. 
-Therefore the two bulk formulea provided include the computation of the fluxes over both 
+The atmospheric fields used depend on the bulk formulae used. Three bulk formulations 
+are available : the CORE, the CLIO and the MFS bulk formulea. The choice is made by setting to true
+one of the following namelist variable : \np{ln\_core} ; \np{ln\_clio} or  \np{ln\_ecmwf}.
+
+Note : in forced mode, when a sea-ice model is used, a bulk formulation (CLIO or CORE) have to be used. 
+Therefore the two bulk (CLIO and CORE) formulea include the computation of the fluxes over both 
 an ocean and an ice surface. 
 
@@ -583,4 +584,42 @@
 namelist (see \S\ref{SBC_fldread}). 
 
+% -------------------------------------------------------------------------------------------------------------
+%        ECMWF Bulk formulea
+% -------------------------------------------------------------------------------------------------------------
+\subsection    [MFS Bulk formulea (\np{ln\_ecmwf}=true)]
+		      {MFS Bulk formulea (\np{ln\_ecmwf}=true, \mdl{sbcblk\_ecmwf})}
+\label{SBC_blk_ecmwf}
+%------------------------------------------namsbc_ecmwf----------------------------------------------------
+\namdisplay{namsbc_ecmwf} 
+%----------------------------------------------------------------------------------------------------------
+
+The MFS (Mediterranean Forecasting System) bulk formulae have been developed by
+ \citet{Castellari_al_JMS1998}. 
+They have been designed to handle the ECMWF operational data and are currently 
+in use in the MFS operational system \citep{Tonani_al_OS08}, \citep{Oddo_al_OS09}.
+The wind stress computation uses a drag coefficient computed according to \citet{Hellerman_Rosenstein_JPO83}.
+The surface boundary condition for temperature involves the balance between surface solar radiation,
+net long-wave radiation, the latent and sensible heat fluxes.
+Solar radiation is dependent on cloud cover and is computed by means of
+an astronomical formula \citep{Reed_JPO77}. Albedo monthly values are from \citet{Payne_JAS72} 
+as means of the values at $40^{o}N$ and $30^{o}N$ for the Atlantic Ocean (hence the same latitudinal
+band of the Mediterranean Sea). The net long-wave radiation flux
+\citep{Bignami_al_JGR95} is a function of
+air temperature, sea-surface temperature, cloud cover and relative humidity.
+Sensible heat and latent heat fluxes are computed by classical
+bulk formulae parameterized according to \citet{Kondo1975}.
+Details on the bulk formulae used can be found in \citet{Maggiore_al_PCE98} and \citet{Castellari_al_JMS1998}.
+
+The required 7 input fields must be provided on the model Grid-T and  are:
+\begin{itemize}
+\item          Zonal Component of the 10m wind ($ms^{-1}$)  (\np{sn\_windi})
+\item          Meridional Component of the 10m wind ($ms^{-1}$)  (\np{sn\_windj})
+\item          Total Claud Cover (\%)  (\np{sn\_clc})
+\item          2m Air Temperature ($K$) (\np{sn\_tair})
+\item          2m Dew Point Temperature ($K$)  (\np{sn\_rhm})
+\item          Total Precipitation ${Kg} m^{-2} s^{-1}$ (\np{sn\_prec})
+\item          Mean Sea Level Pressure (${Pa}) (\np{sn\_msl})
+\end{itemize}
+% -------------------------------------------------------------------------------------------------------------
 % ================================================================
 % Coupled formulation
@@ -938,4 +977,17 @@
 \end{description}
 
+% -------------------------------------------------------------------------------------------------------------
+%        Neutral Drag Coefficient from external wave model
+% -------------------------------------------------------------------------------------------------------------
+\subsection   [Neutral drag coefficient from external wave model (\textit{sbcwave})]
+                        {Neutral drag coefficient from external wave model (\mdl{sbcwave})}
+\label{SBC_wave}
+%------------------------------------------namwave----------------------------------------------------
+\namdisplay{namsbc_wave}
+%-------------------------------------------------------------------------------------------------------------
+\begin{description}
+If (\np{ln\_cdgw}~=~true) in namsbc namelist is activated the \mdl{sbcwave} module which contains the routine \np{sbc\_wave}.This routine reads the namelist namsbc\_wave and the neutral drag coefficient. Then using the routine TURB\_CORE\_1Z or TURB\_CORE\_2Z the drag coefficient is computed according to stable/unstable conditions of the air-sea interface starting from the neutral drag coefficient.
+\end{description}
+
 % Griffies doc:
 % When running ocean-ice simulations, we are not explicitly representing land processes, such as rivers, catchment areas, snow accumulation, etc. However, to reduce model drift, it is important to balance the hydrological cycle in ocean-ice models. We thus need to prescribe some form of global normalization to the precipitation minus evaporation plus river runoff. The result of the normalization should be a global integrated zero net water input to the ocean-ice system over a chosen time scale. 
@@ -944,3 +996,2 @@
 
 
-
Index: branches/2011/dev_r2855_INGV2_3_blk_wave/DOC/TexFiles/Chapters/Chap_ZDF.tex
===================================================================
--- branches/2011/dev_r2855_INGV2_3_blk_wave/DOC/TexFiles/Chapters/Chap_ZDF.tex	(revision 2989)
+++ branches/2011/dev_r2855_INGV2_3_blk_wave/DOC/TexFiles/Chapters/Chap_ZDF.tex	(revision 2990)
@@ -100,4 +100,28 @@
 $a=5$ and $n=2$. The last three values can be modified by setting the 
 \np{rn\_avmri}, \np{rn\_alp} and \np{nn\_ric} namelist parameters, respectively.
+
+A simple mixing-layer model to transfer and dissipate the atmospheric
+ forcings (wind-stress and buoyancy fluxes) can be activated setting 
+the \np{ln\_mldw} =.true. in the namelist.
+
+In this case, the local depth of turbulent wind-mixing or "Ekman depth"
+ $h_{e}(x,y,t)$ is evaluated and the vertical eddy coefficients prescribed within this layer.
+
+This depth is assumed proportional to the "depth of frictional influence" that is limited by rotation:
+\begin{equation}
+         h_{e} = Ek \frac {u^{*}} {f_{0}}  	\\
+\end{equation}
+where, $Ek$ is an empirical parameter, $u^{*}$ is the friction velocity and $f_{0}$ is the Coriolis 
+parameter.
+
+In this similarity height relationship, the turbulent friction velocity:
+\begin{equation}
+         u^{*} = \sqrt \frac {|\tau|} {\rho_o}  	\\
+\end{equation}
+
+is computed from the wind stress vector $|\tau|$ and the reference dendity $ \rho_o$.
+The final $h_{e}$ is further constrained by the adjustable bounds \np{rn\_mldmin} and \np{rn\_mldmax}.
+Once $h_{e}$ is computed, the vertical eddy coefficients within $h_{e}$ are set to 
+the empirical values \np{rn\_wtmix} and \np{rn\_wvmix} \citep{Lermusiaux2001}.
 
 % -------------------------------------------------------------------------------------------------------------
Index: branches/2011/dev_r2855_INGV2_3_blk_wave/DOC/TexFiles/Namelist/namsbc
===================================================================
--- branches/2011/dev_r2855_INGV2_3_blk_wave/DOC/TexFiles/Namelist/namsbc	(revision 2989)
+++ branches/2011/dev_r2855_INGV2_3_blk_wave/DOC/TexFiles/Namelist/namsbc	(revision 2990)
@@ -8,4 +8,5 @@
    ln_blk_clio = .false.   !  CLIO bulk formulation                     (T => fill namsbc_clio) 
    ln_blk_core = .true.    !  CORE bulk formulation                     (T => fill namsbc_core) 
+   ln_blk_ecmwf= .false.   !  MFS bulk formulation                      (T => fill namsbc_ecmwf)
    ln_cpl      = .false.   !  Coupled formulation                       (T => fill namsbc_cpl )
    ln_apr_dyn  = .false.   !  Patm gradient added in ocean & ice Eqs.   (T => fill namsbc_apr )
@@ -20,3 +21,4 @@
                            !     =2 annual global mean of e-p-r set to zero
                            !     =3 global emp set to zero and spread out over erp area
+   ln_cdgw = .true.        ! true if neutral drag coefficient read from wave model
 /
Index: branches/2011/dev_r2855_INGV2_3_blk_wave/DOC/TexFiles/Namelist/namsbc_wave
===================================================================
--- branches/2011/dev_r2855_INGV2_3_blk_wave/DOC/TexFiles/Namelist/namsbc_wave	(revision 2990)
+++ branches/2011/dev_r2855_INGV2_3_blk_wave/DOC/TexFiles/Namelist/namsbc_wave	(revision 2990)
@@ -0,0 +1,9 @@
+!-----------------------------------------------------------------------
+&namsbc_wave   ! External fields from wave model
+!-----------------------------------------------------------------------
+!              !  file name  ! frequency (hours) ! variable  ! time interp. !  clim  ! 'yearly'/ ! weights  ! rotation !
+!              !             !  (if <0  months)  !   name    !   (logical)  !  (T/F) ! 'monthly' ! filename ! pairing  !
+   sn_cdg      =  'cdg_wave' ,        1          , 'drag_coeff' , .true.    , .false. , 'daily'  , ''       , ''
+!
+   cn_dir_cdg  = './'  !  root directory for the location of drag coefficient files
+/
Index: branches/2011/dev_r2855_INGV2_3_blk_wave/DOC/TexFiles/Namelist/namzdf_ric
===================================================================
--- branches/2011/dev_r2855_INGV2_3_blk_wave/DOC/TexFiles/Namelist/namzdf_ric	(revision 2989)
+++ branches/2011/dev_r2855_INGV2_3_blk_wave/DOC/TexFiles/Namelist/namzdf_ric	(revision 2990)
@@ -5,3 +5,9 @@
    rn_alp      =   5.      !  coefficient of the parameterization
    nn_ric      =   2       !  coefficient of the parameterization
+   rn_ekmfc    =   0.7     !  Factor in the Ekman depth Equation
+   rn_mldmin   =   1.0     !  minimum allowable mixed-layer depth estimate (m)
+   rn_mldmax   =1000.0     !  maximum allowable mixed-layer depth estimate (m)
+   rn_wtmix    =  10.0     !  vertical eddy viscosity coeff [m2/s] in the mixed-layer
+   rn_wvmix    =  10.0     !  vertical eddy diffusion coeff [m2/s] in the mixed-layer
+   ln_mldw     = .true.    !  Flag to use or not the mized layer depth param.
 /
Index: branches/2011/dev_r2855_INGV2_3_blk_wave/NEMOGCM/CONFIG/GYRE/EXP00/namelist
===================================================================
--- branches/2011/dev_r2855_INGV2_3_blk_wave/NEMOGCM/CONFIG/GYRE/EXP00/namelist	(revision 2989)
+++ branches/2011/dev_r2855_INGV2_3_blk_wave/NEMOGCM/CONFIG/GYRE/EXP00/namelist	(revision 2990)
@@ -138,4 +138,5 @@
    ln_blk_clio = .false.   !  CLIO bulk formulation                     (T => fill namsbc_clio) 
    ln_blk_core = .false.   !  CORE bulk formulation                     (T => fill namsbc_core) 
+   ln_blk_ecmwf= .false.   !  MFS bulk formulation                      (T => fill namsbc_ecmwf)
    ln_cpl      = .false.   !  Coupled formulation                       (T => fill namsbc_cpl )
    ln_apr_dyn  = .false.   !  Patm gradient added in ocean & ice Eqs.   (T => fill namsbc_apr )
@@ -150,4 +151,5 @@
                            !     =2 annual global mean of e-p-r set to zero
                            !     =3 global emp set to zero and spread out over erp area
+   ln_cdgw = .false.       !  Neutral drag coefficient read from wave model (T => fill namsbc_wave)
 /
 !-----------------------------------------------------------------------
@@ -208,4 +210,19 @@
    ln_taudif   = .false.   !  HF tau contribution: use "mean of stress module - module of the mean stress" data
    rn_pfac     = 1.        !  multiplicative factor for precipitation (total & snow)
+/
+!-----------------------------------------------------------------------
+&namsbc_ecmwf   !   namsbc_ecmwf  MFS bulk formulea
+!-----------------------------------------------------------------------
+!              !  file name  ! frequency (hours) ! variable  ! time interp. !  clim  ! 'yearly'/ ! weights  ! rotation !
+!              !             !  (if <0  months)  !   name    !   (logical)  !  (T/F) ! 'monthly' ! filename ! pairing  !
+   sn_wndi     =   'ecmwf'   ,        6          , 'u10'     ,    .true.    , .false. , 'daily'  ,'bicubic.nc' , ''
+   sn_wndj     =   'ecmwf'   ,        6          , 'v10'     ,    .true.    , .false. , 'daily'  ,'bicubic.nc' , ''
+   sn_clc      =   'ecmwf'   ,        6          , 'clc'     ,    .true.    , .false. , 'daily'  ,'bilinear.nc', ''
+   sn_msl      =   'ecmwf'   ,        6          , 'msl'     ,    .true.    , .false. , 'daily'  ,'bicubic.nc' , ''
+   sn_tair     =   'ecmwf'   ,        6          , 't2'      ,    .true.    , .false. , 'daily'  ,'bicubic.nc' , ''
+   sn_rhm      =   'ecmwf'   ,        6          , 'rh'      ,    .true.    , .false. , 'daily'  ,'bilinear.nc', ''
+   sn_prec     =   'ecmwf'   ,        6          , 'precip'  ,    .true.    , .true.  , 'daily'  ,'bicubic.nc' , ''
+
+   cn_dir      = './ECMWF/'      !  root directory for the location of the bulk files
 /
 !-----------------------------------------------------------------------
@@ -863,2 +880,11 @@
     salfixmin = -9999      !  Minimum salinity after applying the increments
 /
+!-----------------------------------------------------------------------
+&namsbc_wave   ! External fields from wave model
+!-----------------------------------------------------------------------
+!              !  file name  ! frequency (hours) ! variable  ! time interp. !  clim  ! 'yearly'/ ! weights  ! rotation !
+!              !             !  (if <0  months)  !   name    !   (logical)  !  (T/F) ! 'monthly' ! filename ! pairing  !
+   sn_cdg      =  'cdg_wave' ,        1          , 'drag_coeff' , .true.   , .false. , 'daily'  ,''         , ''
+!
+   cn_dir_cdg  = './'  !  root directory for the location of drag coefficient files
+/
Index: branches/2011/dev_r2855_INGV2_3_blk_wave/NEMOGCM/CONFIG/ORCA2_LIM/EXP00/namelist
===================================================================
--- branches/2011/dev_r2855_INGV2_3_blk_wave/NEMOGCM/CONFIG/ORCA2_LIM/EXP00/namelist	(revision 2989)
+++ branches/2011/dev_r2855_INGV2_3_blk_wave/NEMOGCM/CONFIG/ORCA2_LIM/EXP00/namelist	(revision 2990)
@@ -138,4 +138,5 @@
    ln_blk_clio = .false.   !  CLIO bulk formulation                     (T => fill namsbc_clio) 
    ln_blk_core = .true.    !  CORE bulk formulation                     (T => fill namsbc_core) 
+   ln_blk_ecmwf= .false.   !  MFS bulk formulation                      (T => fill namsbc_ecmwf)
    ln_cpl      = .false.   !  Coupled formulation                       (T => fill namsbc_cpl )
    ln_apr_dyn  = .false.   !  Patm gradient added in ocean & ice Eqs.   (T => fill namsbc_apr )
@@ -150,4 +151,5 @@
                            !     =2 annual global mean of e-p-r set to zero
                            !     =3 global emp set to zero and spread out over erp area
+   ln_cdgw = .false.       !  Neutral drag coefficient read from wave model (T => fill namsbc_wave)
 /
 !-----------------------------------------------------------------------
@@ -208,4 +210,19 @@
    ln_taudif   = .false.   !  HF tau contribution: use "mean of stress module - module of the mean stress" data
    rn_pfac     = 1.        !  multiplicative factor for precipitation (total & snow)
+/
+!-----------------------------------------------------------------------
+&namsbc_ecmwf   !   namsbc_ecmwf  MFS bulk formulea
+!-----------------------------------------------------------------------
+!              !  file name  ! frequency (hours) ! variable  ! time interp. !  clim  ! 'yearly'/ ! weights  ! rotation !
+!              !             !  (if <0  months)  !   name    !   (logical)  !  (T/F) ! 'monthly' ! filename ! pairing  !
+   sn_wndi     =   'ecmwf'   ,        6          , 'u10'     ,    .true.    , .false. , 'daily'  ,'bicubic.nc' , ''
+   sn_wndj     =   'ecmwf'   ,        6          , 'v10'     ,    .true.    , .false. , 'daily'  ,'bicubic.nc' , ''
+   sn_clc      =   'ecmwf'   ,        6          , 'clc'     ,    .true.    , .false. , 'daily'  ,'bilinear.nc', ''
+   sn_msl      =   'ecmwf'   ,        6          , 'msl'     ,    .true.    , .false. , 'daily'  ,'bicubic.nc' , ''
+   sn_tair     =   'ecmwf'   ,        6          , 't2'      ,    .true.    , .false. , 'daily'  ,'bicubic.nc' , ''
+   sn_rhm      =   'ecmwf'   ,        6          , 'rh'      ,    .true.    , .false. , 'daily'  ,'bilinear.nc', ''
+   sn_prec     =   'ecmwf'   ,        6          , 'precip'  ,    .true.    , .true.  , 'daily'  ,'bicubic.nc' , ''
+
+   cn_dir      = './ECMWF/'      !  root directory for the location of the bulk files
 /
 !-----------------------------------------------------------------------
@@ -863,2 +880,11 @@
     salfixmin = -9999      !  Minimum salinity after applying the increments
 /
+!-----------------------------------------------------------------------
+&namsbc_wave   ! External fields from wave model
+!-----------------------------------------------------------------------
+!              !  file name  ! frequency (hours) ! variable  ! time interp. !  clim  ! 'yearly'/ ! weights  ! rotation !
+!              !             !  (if <0  months)  !   name    !   (logical)  !  (T/F) ! 'monthly' ! filename ! pairing  !
+   sn_cdg      =  'cdg_wave' ,        1          , 'drag_coeff' , .true.   , .false. , 'daily'  ,''         , ''
+!
+   cn_dir_cdg  = './'  !  root directory for the location of drag coefficient files
+/
Index: branches/2011/dev_r2855_INGV2_3_blk_wave/NEMOGCM/CONFIG/ORCA2_OFF_PISCES/EXP00/namelist
===================================================================
--- branches/2011/dev_r2855_INGV2_3_blk_wave/NEMOGCM/CONFIG/ORCA2_OFF_PISCES/EXP00/namelist	(revision 2989)
+++ branches/2011/dev_r2855_INGV2_3_blk_wave/NEMOGCM/CONFIG/ORCA2_OFF_PISCES/EXP00/namelist	(revision 2990)
@@ -138,4 +138,5 @@
    ln_blk_clio = .false.   !  CLIO bulk formulation                     (T => fill namsbc_clio) 
    ln_blk_core = .true.    !  CORE bulk formulation                     (T => fill namsbc_core) 
+   ln_blk_ecmwf= .false.   !  MFS bulk formulation                      (T => fill namsbc_ecmwf)
    ln_cpl      = .false.   !  Coupled formulation                       (T => fill namsbc_cpl )
    ln_apr_dyn  = .false.   !  Patm gradient added in ocean & ice Eqs.   (T => fill namsbc_apr )
@@ -150,4 +151,5 @@
                            !     =2 annual global mean of e-p-r set to zero
                            !     =3 global emp set to zero and spread out over erp area
+   ln_cdgw = .false.       !  Neutral drag coefficient read from wave model (T => fill namsbc_wave )
 /
 !-----------------------------------------------------------------------
@@ -208,4 +210,19 @@
    ln_taudif   = .false.   !  HF tau contribution: use "mean of stress module - module of the mean stress" data ?
    rn_pfac     = 1.        !  multiplicative factor for precipitation (total & snow)
+/
+!-----------------------------------------------------------------------
+&namsbc_ecmwf   !   namsbc_ecmwf  MFS bulk formulea
+!-----------------------------------------------------------------------
+!              !  file name  ! frequency (hours) ! variable  ! time interp. !  clim  ! 'yearly'/ ! weights  ! rotation !
+!              !             !  (if <0  months)  !   name    !   (logical)  !  (T/F) ! 'monthly' ! filename ! pairing  !
+   sn_wndi     =   'ecmwf'   ,        6          , 'u10'     ,    .true.    , .false. , 'daily'  ,'bicubic.nc' , ''
+   sn_wndj     =   'ecmwf'   ,        6          , 'v10'     ,    .true.    , .false. , 'daily'  ,'bicubic.nc' , ''
+   sn_clc      =   'ecmwf'   ,        6          , 'clc'     ,    .true.    , .false. , 'daily'  ,'bilinear.nc', ''
+   sn_msl      =   'ecmwf'   ,        6          , 'msl'     ,    .true.    , .false. , 'daily'  ,'bicubic.nc' , ''
+   sn_tair     =   'ecmwf'   ,        6          , 't2'      ,    .true.    , .false. , 'daily'  ,'bicubic.nc' , ''
+   sn_rhm      =   'ecmwf'   ,        6          , 'rh'      ,    .true.    , .false. , 'daily'  ,'bilinear.nc', ''
+   sn_prec     =   'ecmwf'   ,        6          , 'precip'  ,    .true.    , .true.  , 'daily'  ,'bicubic.nc' , ''
+
+   cn_dir      = './ECMWF/'      !  root directory for the location of the bulk files
 /
 !-----------------------------------------------------------------------
@@ -878,2 +895,11 @@
     salfixmin = -9999      !  Minimum salinity after applying the increments
 /
+!-----------------------------------------------------------------------
+&namsbc_wave   ! External fields from wave model
+!-----------------------------------------------------------------------
+!              !  file name  ! frequency (hours) ! variable  ! time interp. !  clim  ! 'yearly'/ ! weights  ! rotation !
+!              !             !  (if <0  months)  !   name    !   (logical)  !  (T/F) ! 'monthly' ! filename ! pairing  !
+   sn_cdg      =  'cdg_wave' ,        1          , 'drag_coeff' , .true.   , .false. , 'daily'  ,''         , ''
+!
+   cn_dir_cdg  = './'  !  root directory for the location of drag coefficient files
+/
Index: branches/2011/dev_r2855_INGV2_3_blk_wave/NEMOGCM/CONFIG/POMME/EXP00/namelist
===================================================================
--- branches/2011/dev_r2855_INGV2_3_blk_wave/NEMOGCM/CONFIG/POMME/EXP00/namelist	(revision 2989)
+++ branches/2011/dev_r2855_INGV2_3_blk_wave/NEMOGCM/CONFIG/POMME/EXP00/namelist	(revision 2990)
@@ -138,4 +138,5 @@
    ln_blk_clio = .false.   !  CLIO bulk formulation                     (T => fill namsbc_clio) 
    ln_blk_core = .true.    !  CORE bulk formulation                     (T => fill namsbc_core) 
+   ln_blk_ecmwf= .false.   !  MFS bulk formulation                      (T => fill namsbc_ecmwf)
    ln_cpl      = .false.   !  Coupled formulation                       (T => fill namsbc_cpl )
    ln_apr_dyn  = .false.   !  Patm gradient added in ocean & ice Eqs.   (T => fill namsbc_apr )
@@ -150,4 +151,5 @@
                            !     =2 annual global mean of e-p-r set to zero
                            !     =3 global emp set to zero and spread out over erp area
+   ln_cdgw = .false.       !  Neutral drag coefficient read from wave model (T => fill namsbc_wave)
 /
 !-----------------------------------------------------------------------
@@ -208,4 +210,19 @@
    ln_taudif   = .false.   !  HF tau contribution: use "mean of stress module - module of the mean stress" data
    rn_pfac     = 1.        !  multiplicative factor for precipitation (total & snow)
+/
+!-----------------------------------------------------------------------
+&namsbc_ecmwf   !   namsbc_ecmwf  MFS bulk formulea
+!-----------------------------------------------------------------------
+!              !  file name  ! frequency (hours) ! variable  ! time interp. !  clim  ! 'yearly'/ ! weights  ! rotation !
+!              !             !  (if <0  months)  !   name    !   (logical)  !  (T/F) ! 'monthly' ! filename ! pairing  !
+   sn_wndi     =   'ecmwf'   ,        6          , 'u10'     ,    .true.    , .false. , 'daily'  ,'bicubic.nc' , ''
+   sn_wndj     =   'ecmwf'   ,        6          , 'v10'     ,    .true.    , .false. , 'daily'  ,'bicubic.nc' , ''
+   sn_clc      =   'ecmwf'   ,        6          , 'clc'     ,    .true.    , .false. , 'daily'  ,'bilinear.nc', ''
+   sn_msl      =   'ecmwf'   ,        6          , 'msl'     ,    .true.    , .false. , 'daily'  ,'bicubic.nc' , ''
+   sn_tair     =   'ecmwf'   ,        6          , 't2'      ,    .true.    , .false. , 'daily'  ,'bicubic.nc' , ''
+   sn_rhm      =   'ecmwf'   ,        6          , 'rh'      ,    .true.    , .false. , 'daily'  ,'bilinear.nc', ''
+   sn_prec     =   'ecmwf'   ,        6          , 'precip'  ,    .true.    , .true.  , 'daily'  ,'bicubic.nc' , ''
+
+   cn_dir      = './ECMWF/'      !  root directory for the location of the bulk files
 /
 !-----------------------------------------------------------------------
@@ -868,2 +885,11 @@
     salfixmin = -9999      !  Minimum salinity after applying the increments
 /
+!-----------------------------------------------------------------------
+&namsbc_wave   ! External fields from wave model
+!-----------------------------------------------------------------------
+!              !  file name  ! frequency (hours) ! variable  ! time interp. !  clim  ! 'yearly'/ ! weights  ! rotation !
+!              !             !  (if <0  months)  !   name    !   (logical)  !  (T/F) ! 'monthly' ! filename ! pairing  !
+   sn_cdg      =  'cdg_wave' ,        1          , 'drag_coeff' , .true.   , .false. , 'daily'  ,''         , ''
+!
+   cn_dir_cdg  = './'  !  root directory for the location of drag coefficient files
+/
Index: branches/2011/dev_r2855_INGV2_3_blk_wave/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90
===================================================================
--- branches/2011/dev_r2855_INGV2_3_blk_wave/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90	(revision 2989)
+++ branches/2011/dev_r2855_INGV2_3_blk_wave/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90	(revision 2990)
@@ -33,4 +33,5 @@
    LOGICAL , PUBLIC ::   ln_blk_clio = .FALSE.   !: CLIO bulk formulation
    LOGICAL , PUBLIC ::   ln_blk_core = .FALSE.   !: CORE bulk formulation
+   LOGICAL , PUBLIC ::   ln_blk_ecmwf= .FALSE.   !: MFS  bulk formulation
    LOGICAL , PUBLIC ::   ln_cpl      = .FALSE.   !: coupled   formulation (overwritten by key_sbc_coupled )
    LOGICAL , PUBLIC ::   ln_dm2dc    = .FALSE.   !: Daily mean to Diurnal Cycle short wave (qsr)
@@ -43,4 +44,5 @@
    !                                             !:  = 1 global mean of e-p-r set to zero at each nn_fsbc time step
    !                                             !:  = 2 annual global mean of e-p-r set to zero
+   LOGICAL , PUBLIC ::   ln_cdgw     = .FALSE.   !: true if neutral drag coefficient read from wave model
 
    !!----------------------------------------------------------------------
Index: branches/2011/dev_r2855_INGV2_3_blk_wave/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90
===================================================================
--- branches/2011/dev_r2855_INGV2_3_blk_wave/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90	(revision 2989)
+++ branches/2011/dev_r2855_INGV2_3_blk_wave/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90	(revision 2990)
@@ -34,4 +34,5 @@
    USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
    USE prtctl          ! Print control
+   USE sbcwave,ONLY :  cdn_wave !wave module 
 #if defined key_lim3
    USE sbc_ice         ! Surface boundary condition: ice fields
@@ -43,4 +44,5 @@
    PUBLIC   sbc_blk_core         ! routine called in sbcmod module
    PUBLIC   blk_ice_core         ! routine called in sbc_ice_lim module
+   PUBLIC   turb_core_2z         ! routine calles in sbcblk_ecmwf module
 
    INTEGER , PARAMETER ::   jpfld   = 9           ! maximum number of files to read 
@@ -682,5 +684,10 @@
       !! Neutral Drag Coefficient
       stab    = 0.5 + sign(0.5,dT)    ! stable : stab = 1 ; unstable : stab = 0 
-      Cd_n10  = 1E-3 * ( 2.7/dU10 + 0.142 + dU10/13.09 )    !   L & Y eq. (6a)
+      IF  ( ln_cdgw ) THEN
+        cdn_wave = cdn_wave - rsmall*(tmask(:,:,1)-1)
+        Cd_n10(:,:) =   cdn_wave
+      ELSE
+        Cd_n10  = 1E-3 * ( 2.7/dU10 + 0.142 + dU10/13.09 )    !   L & Y eq. (6a)
+      ENDIF
       sqrt_Cd_n10 = sqrt(Cd_n10)
       Ce_n10  = 1E-3 * ( 34.6 * sqrt_Cd_n10 )               !   L & Y eq. (6b)
@@ -705,18 +712,22 @@
          zpsi_m  = psi_m(zeta)
 
-         !! Shifting the wind speed to 10m and neutral stability :
-         U_n10 = dU10*1./(1. + sqrt_Cd_n10/kappa*(log(zu/10.) - zpsi_m)) !  L & Y eq. (9a)
-
-         !! Updating the neutral 10m transfer coefficients :
-         Cd_n10  = 1E-3 * (2.7/U_n10 + 0.142 + U_n10/13.09)              !  L & Y eq. (6a)
-         sqrt_Cd_n10 = sqrt(Cd_n10)
-         Ce_n10  = 1E-3 * (34.6 * sqrt_Cd_n10)                           !  L & Y eq. (6b)
-         stab    = 0.5 + sign(0.5,zeta)
-         Ch_n10  = 1E-3*sqrt_Cd_n10*(18.*stab + 32.7*(1-stab))           !  L & Y eq. (6c), (6d)
-
-         !! Shifting the neutral  10m transfer coefficients to ( zu , zeta ) :
-         !!
-         xct = 1. + sqrt_Cd_n10/kappa*(log(zu/10) - zpsi_m)
-         Cd  = Cd_n10/(xct*xct) ;  sqrt_Cd = sqrt(Cd)
+         IF  ( ln_cdgw ) THEN
+           sqrt_Cd=kappa/((kappa/sqrt_Cd_n10) - zpsi_m) ; Cd=sqrt_Cd*sqrt_Cd;
+         ELSE
+           !! Shifting the wind speed to 10m and neutral stability :
+           U_n10 = dU10*1./(1. + sqrt_Cd_n10/kappa*(log(zu/10.) - zpsi_m)) !  L & Y eq. (9a)
+
+           !! Updating the neutral 10m transfer coefficients :
+           Cd_n10  = 1E-3 * (2.7/U_n10 + 0.142 + U_n10/13.09)              !  L & Y eq. (6a)
+           sqrt_Cd_n10 = sqrt(Cd_n10)
+           Ce_n10  = 1E-3 * (34.6 * sqrt_Cd_n10)                           !  L & Y eq. (6b)
+           stab    = 0.5 + sign(0.5,zeta)
+           Ch_n10  = 1E-3*sqrt_Cd_n10*(18.*stab + 32.7*(1-stab))           !  L & Y eq. (6c), (6d)
+
+           !! Shifting the neutral  10m transfer coefficients to ( zu , zeta ) :
+           !!
+           xct = 1. + sqrt_Cd_n10/kappa*(log(zu/10) - zpsi_m)
+           Cd  = Cd_n10/(xct*xct) ;  sqrt_Cd = sqrt(Cd)
+         ENDIF
          !!
          xlogt = log(zu/10.) - zpsi_h
@@ -812,5 +823,10 @@
       !! Neutral Drag Coefficient :
       stab = 0.5 + sign(0.5,dT)                 ! stab = 1  if dT > 0  -> STABLE
-      Cd_n10  = 1E-3*( 2.7/dU10 + 0.142 + dU10/13.09 ) 
+      IF  ( ln_cdgw ) THEN
+        cdn_wave = cdn_wave - rsmall*(tmask(:,:,1)-1)
+        Cd_n10(:,:) =   cdn_wave
+      ELSE
+        Cd_n10  = 1E-3*( 2.7/dU10 + 0.142 + dU10/13.09 ) 
+      ENDIF
       sqrt_Cd_n10 = sqrt(Cd_n10)
       Ce_n10  = 1E-3*( 34.6 * sqrt_Cd_n10 )
@@ -853,18 +869,20 @@
          stab = 0.5 + sign(0.5,q_zu) ;  q_zu = stab*q_zu
          !!
-         !! Updating the neutral 10m transfer coefficients :
-         Cd_n10  = 1E-3 * (2.7/U_n10 + 0.142 + U_n10/13.09)    ! L & Y eq. (6a)
-         sqrt_Cd_n10 = sqrt(Cd_n10)
-         Ce_n10  = 1E-3 * (34.6 * sqrt_Cd_n10)                 ! L & Y eq. (6b)
-         stab    = 0.5 + sign(0.5,zeta_u)
-         Ch_n10  = 1E-3*sqrt_Cd_n10*(18.*stab + 32.7*(1-stab)) ! L & Y eq. (6c-6d)
-         !!
-         !!
-         !! Shifting the neutral 10m transfer coefficients to (zu,zeta_u) :
-!        xct = 1. + sqrt_Cd_n10/kappa*(log(zu/10.) - psi_m(zeta_u))
-         xct = 1. + sqrt_Cd_n10/kappa*(log(zu/10.) - zpsi_m)
-         Cd = Cd_n10/(xct*xct) ; sqrt_Cd = sqrt(Cd)
-         !!
-!        xlogt = log(zu/10.) - psi_h(zeta_u)
+         IF  ( ln_cdgw ) THEN
+            sqrt_Cd=kappa/((kappa/sqrt_Cd_n10) - zpsi_m) ; Cd=sqrt_Cd*sqrt_Cd;
+         ELSE
+           !! Updating the neutral 10m transfer coefficients :
+           Cd_n10  = 1E-3 * (2.7/U_n10 + 0.142 + U_n10/13.09)    ! L & Y eq. (6a)
+           sqrt_Cd_n10 = sqrt(Cd_n10)
+           Ce_n10  = 1E-3 * (34.6 * sqrt_Cd_n10)                 ! L & Y eq. (6b)
+           stab    = 0.5 + sign(0.5,zeta_u)
+           Ch_n10  = 1E-3*sqrt_Cd_n10*(18.*stab + 32.7*(1-stab)) ! L & Y eq. (6c-6d)
+           !!
+           !!
+           !! Shifting the neutral 10m transfer coefficients to (zu,zeta_u) :
+           xct = 1. + sqrt_Cd_n10/kappa*(log(zu/10.) - zpsi_m)
+           Cd = Cd_n10/(xct*xct) ; sqrt_Cd = sqrt(Cd)
+         ENDIF
+         !!
          xlogt = log(zu/10.) - zpsi_hu
          !!
Index: branches/2011/dev_r2855_INGV2_3_blk_wave/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_ecmwf.F90
===================================================================
--- branches/2011/dev_r2855_INGV2_3_blk_wave/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_ecmwf.F90	(revision 2990)
+++ branches/2011/dev_r2855_INGV2_3_blk_wave/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_ecmwf.F90	(revision 2990)
@@ -0,0 +1,659 @@
+MODULE sbcblk_ecmwf
+   !!======================================================================
+   !!                       ***  MODULE  sbcblk_ecmwf  ***
+   !! Ocean forcing:  momentum, heat and freshwater flux formulation
+   !!=====================================================================
+   !! History :  3.3  !   2010-05 (P. Oddo) Original Code
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   sbc_blk_ecmwf  : bulk formulation as ocean surface boundary condition
+   !!                   (forced mode, ecmwf bulk formulea)
+   !!   blk_oce_ecmwf  : ocean: computes momentum, heat and freshwater fluxes
+   !!----------------------------------------------------------------------
+   USE oce             ! ocean dynamics and tracers
+   USE dom_oce         ! ocean space and time domain
+   USE phycst          ! physical constants
+   USE fldread         ! read input fields
+   USE sbc_oce         ! Surface boundary condition: ocean fields
+   USE iom             ! I/O manager library
+   USE in_out_manager  ! I/O manager
+   USE lib_mpp         ! distribued memory computing library
+   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
+   USE prtctl          ! Print control
+   USE sbcwave,ONLY : cdn_wave !wave module
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   sbc_blk_ecmwf       ! routine called in sbcmod module
+      
+   INTEGER , PARAMETER ::   jpfld   = 7         ! maximum number of files to read 
+   INTEGER , PARAMETER ::   jp_wndi = 1         ! index of 10m wind velocity (i-component) (m/s) at T-point
+   INTEGER , PARAMETER ::   jp_wndj = 2         ! index of 10m wind velocity (j-component) (m/s) at T-point
+   INTEGER , PARAMETER ::   jp_clc  = 3         ! index of total cloud cover               ( % )
+   INTEGER , PARAMETER ::   jp_msl  = 4         ! index of mean sea level pressure         (Pa)
+   INTEGER , PARAMETER ::   jp_tair = 5         ! index of 10m air temperature             (Kelvin)
+   INTEGER , PARAMETER ::   jp_rhm  = 6         ! index of dew point temperature           (Kelvin)
+   INTEGER , PARAMETER ::   jp_prec = 7         ! index of total precipitation (rain+snow) (Kg/m2/s)
+   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf ! structure of input fields (file informations, fields read)
+         
+   !! * Substitutions
+#  include "domzgr_substitute.h90"
+#  include "vectopt_loop_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009) 
+   !! $Id: sbcblk_ecmwf.F90 1730 2009-11-16 14:34:19Z poddo $
+   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+
+CONTAINS
+
+   SUBROUTINE sbc_blk_ecmwf( kt )
+      !!---------------------------------------------------------------------
+      !!                    ***  ROUTINE sbc_blk_ecmwf  ***
+      !!                   
+      !! ** Purpose :   provide at each time step the surface ocean fluxes
+      !!      (momentum, heat, freshwater, runoff is added later in the code) 
+      !!
+      !! ** Method  : (1) READ Atmospheric data from NetCDF files:
+      !!      the 10m wind velocity (i-component) (m/s) at T-point
+      !!      the 10m wind velocity (j-component) (m/s) at T-point
+      !!      the 2m Dew point Temperature        (k)
+      !!      the Cloud COver                     (%)
+      !!      the 2m air temperature              (Kelvin)
+      !!      the Mean Sea Level Preesure         (hPa)
+      !!      the Climatological Precipitation    (kg/m2/s)
+      !!              (2) CALL blk_oce_ecmwf
+      !!
+      !!      Computes:
+      !!      Solar Radiation using Reed formula (1975, 1977)
+      !!      Net Long wave radiation using Bignami et al. (1995)
+      !!      Latent and Sensible heat using Kondo (1975)
+      !!      Drag coeff using Hllerman and Rosenstein (1983)
+      !!      C A U T I O N : never mask the surface stress fields
+      !!                      the stress is assumed to be in the mesh referential
+      !!                      i.e. the (i,j) referential
+      !!
+      !! ** Action  :   defined at each time-step at the air-sea interface
+      !!              - utau, vtau  i- and j-component of the wind stress
+      !!              - taum        wind stress module at T-point
+      !!              - wndm        10m wind module at T-point
+      !!              - qns, qsr    non-slor and solar heat flux
+      !!              - emp, emps   evaporation minus precipitation
+      !!----------------------------------------------------------------------
+      REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::  sh_now   ! specific humidity at T-point 
+      REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::  catm     ! Cover 
+      REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::  alonl    ! Lon 
+      REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::  alatl    ! Lat 
+      REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::  gsst     ! SST 
+      !!---------------------------------------------------------------------
+      !! Local fluxes variables
+      !!---------------------------------------------------------------------
+      REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::  qbw     ! Net Long wave 
+      REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::  ha      ! Sesnible 
+      REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::  elat    ! Latent 
+      REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::  evap    ! evaporation rate 
+
+      INTEGER, INTENT( in  ) ::   kt   ! ocean time step
+      !!
+      INTEGER  :: ierror                          ! return error code
+      INTEGER  :: ifpr     ! dummy loop indice
+      INTEGER  :: jj,ji    ! dummy loop arguments
+      REAL(wp) :: act_hour
+      !!--------------------------------------------------------------------
+      !! Variables for specific humidity computation
+      !!--------------------------------------------------------------------
+      REAL(wp) :: onsea,par1,par2
+      DATA onsea,par1,par2 / 0.98, 640380., -5107.4 /
+      !!                      par1 [Kg/m3], par2 [K]
+
+      CHARACTER(len=100) ::  cn_dir                           ! Root directory for location of ecmwf files
+      TYPE(FLD_N), DIMENSION(jpfld) ::   slf_i                ! array of namelist informations on the fields to read
+      TYPE(FLD_N) ::   sn_wndi, sn_wndj, sn_clc, sn_msl       ! informations about the fields to be read
+      TYPE(FLD_N) ::   sn_tair , sn_rhm, sn_prec              !   "                                 "
+      !!---------------------------------------------------------------------
+      NAMELIST/namsbc_ecmwf/ cn_dir ,                                        &
+         &                  sn_wndi , sn_wndj, sn_clc   , sn_msl ,           &
+         &                  sn_tair , sn_rhm , sn_prec 
+      !!---------------------------------------------------------------------
+
+      !                                         ! ====================== !
+      IF( kt == nit000 ) THEN                   !  First call kt=nit000  !
+         !                                      ! ====================== !
+      ALLOCATE( sh_now(jpi,jpj), catm(jpi,jpj), alonl(jpi,jpj), alatl(jpi,jpj),     &
+         &        gsst(jpi,jpj),  qbw(jpi,jpj),    ha(jpi,jpj),  elat(jpi,jpj),     &
+         &        evap(jpi,jpj), STAT=ierror )
+
+      IF( ierror /= 0 )   CALL ctl_warn('sbc_blk_ecmwf: failed to allocate arrays')
+         !
+
+         ! set file information (default values)
+         cn_dir = './'       ! directory in which the model is executed
+         !
+         ! (NB: frequency positive => hours, negative => months)
+         !            !    file     ! frequency !  variable  ! time intep !  clim   ! 'yearly' or ! weights  ! rotation   !
+         !            !    name     !  (hours)  !   name     !   (T/F)    !  (T/F)  !  'monthly'  ! filename ! pairs      !
+         sn_wndi = FLD_N( 'ecmwf'   ,    24     ,  'u10'     ,  .false.   , .false. ,   'daily'   , ''       , ''         )
+         sn_wndj = FLD_N( 'ecmwf'   ,    24     ,  'v10'     ,  .false.   , .false. ,   'daily'   , ''       , ''         )
+         sn_clc  = FLD_N( 'ecmwf'   ,    24     ,  'clc'     ,  .false.   , .false. ,   'daily'   , ''       , ''         )
+         sn_msl  = FLD_N( 'ecmwf'   ,    24     ,  'msl'     ,  .false.   , .false. ,   'daily'   , ''       , ''         )
+         sn_tair = FLD_N( 'ecmwf'   ,    24     ,  't2'      ,  .false.   , .false. ,   'daily'   , ''       , ''         )
+         sn_rhm  = FLD_N( 'ecmwf'   ,    24     ,  'rh'      ,  .false.   , .false. ,   'daily'   , ''       , ''         )
+         sn_prec = FLD_N( 'precip_cmap' ,  -1   ,  'precip'  ,  .true.    ,  .true. ,   'yearly'  , ''       , ''         )
+         !
+         REWIND( numnam )                    ! ... read in namlist namsbc_ecmwf
+         READ  ( numnam, namsbc_ecmwf )
+         !
+         ! store namelist information in an array
+         slf_i(jp_wndi) = sn_wndi   ;   slf_i(jp_wndj) = sn_wndj
+         slf_i(jp_clc ) = sn_clc    ;   slf_i(jp_msl ) = sn_msl
+         slf_i(jp_tair) = sn_tair   ;   slf_i(jp_rhm)  = sn_rhm
+         slf_i(jp_prec) = sn_prec   ;  
+         !
+         ALLOCATE( sf(jpfld), STAT=ierror )         ! set sf structure
+         IF( ierror > 0 ) THEN
+            CALL ctl_stop( 'sbc_blk_ecmwf: unable to allocate sf structure' )   ;   RETURN
+         ENDIF
+         DO ifpr= 1, jpfld
+            ALLOCATE( sf(ifpr)%fnow(jpi,jpj,1) )
+            IF( slf_i(ifpr)%ln_tint ) ALLOCATE( sf(ifpr)%fdta(jpi,jpj,1,2) )
+         END DO
+         ! fill sf with slf_i and control print
+         CALL fld_fill( sf, slf_i, cn_dir,'sbc_blk_ecmwf','bulk formulation for ocean SBC', 'namsbc_ecmwf' )
+         !
+      ENDIF
+
+      CALL fld_read( kt, nn_fsbc, sf )                   ! input fields provided at the current time-step
+
+      catm(:,:)   = 0.0    ! initializze cloud cover variable
+      sh_now(:,:) = 0.0    ! initializze specifif humidity variable
+
+         DO JJ=2,jpj-1           ! Here loop over 2 jpj-1 to avoid computation on ghost cells
+          DO JI=2,jpi-1          ! Here loop over 2 jpi-1 to avoid computation on ghost cells
+
+         ! Calculate Specific Humidity 
+         !-------------------------------------------------
+           sh_now(ji,jj) = (1/1.22) * onsea * par1 * EXP(par2/sf(jp_rhm)%fnow(ji,jj,1))
+
+         ! Normalize Clouds
+         !-------------------------------------------------
+           catm(ji,jj)   = max(0.0,min(1.0,sf(jp_clc)%fnow(ji,jj,1)*0.01))
+
+          END DO
+         END DO
+
+         ! wind module at 10m 
+         !--------------------------------------------------
+           wndm(:,:) = SQRT(  sf(jp_wndi)%fnow(:,:,1) * sf(jp_wndi)%fnow(:,:,1)   &
+              &             + sf(jp_wndj)%fnow(:,:,1) * sf(jp_wndj)%fnow(:,:,1)  )
+
+         ! Some conv for fluxes computation
+         !-------------------------------------------------
+         alonl(:,:) = glamt(:,:) * rad
+         alatl(:,:) = gphit(:,:) * rad
+         gsst(:,:)  = tn(:,:,1)  * tmask(:,:,1)
+
+      IF( MOD( kt - 1, nn_fsbc ) == 0 ) THEN
+
+         ! Force to zero the output of fluxes 
+         !-------------------------------------------------
+          qsr(:,:)  = 0.0 ; qbw(:,:)  = 0.0 ; 
+          ha(:,:)   = 0.0 ; elat(:,:) = 0.0 ; 
+          evap(:,:) = 0.0 ; utau(:,:) = 0.0 ; 
+          vtau(:,:) = 0.0
+
+      CALL lbc_lnk( sf(jp_wndi)%fnow(:,:,1), 'T', -1. )
+      CALL lbc_lnk( sf(jp_wndj)%fnow(:,:,1), 'T', -1. )
+
+      act_hour = (( nsec_year / rday ) - INT (nsec_year / rday)) * rjjhh
+
+                call fluxes_ecmwf(alatl,alonl,act_hour,                              &     ! input static
+                                  gsst(:,:),sf(jp_tair)%fnow(:,:,1),sh_now(:,:),     &     ! input dynamic
+                                  sf(jp_wndi)%fnow(:,:,1), sf(jp_wndj)%fnow(:,:,1),  &     ! input dynamic
+                                  sf(jp_msl)%fnow(:,:,1) , catm(:,:) ,               &     ! input dynamic
+                                  qsr,qbw,ha,elat,evap,utau,vtau)                          ! output
+
+         ! Shortwave radiation
+         !--------------------------------------------------
+           qsr(:,:) = qsr(:,:) * tmask(:,:,1)
+
+         ! total non solar heat flux over water
+         !--------------------------------------------------
+           qns(:,:) = -1. * ( qbw(:,:) + ha(:,:) + elat(:,:) )
+           qns(:,:) = qns(:,:)*tmask(:,:,1)
+
+         ! mask the wind module at 10m 
+         !--------------------------------------------------
+           wndm(:,:) = wndm(:,:) * tmask(:,:,1)
+
+         !   wind stress module (taum) into T-grid
+         !--------------------------------------------------
+           taum(:,:) = SQRT( utau(:,:) * utau(:,:) + vtau(:,:) * vtau(:,:) ) * tmask(:,:,1)
+
+         CALL lbc_lnk( taum, 'T', 1. )
+
+         ! Interpolate utau, vtau into the grid_V and grid_V
+         !-------------------------------------------------
+
+         DO jj = 1, jpjm1
+          DO ji = 1, fs_jpim1
+            utau(ji,jj) = 0.5 * ( 2. - umask(ji,jj,1) ) * ( utau(ji,jj) * tmask(ji,jj,1) &
+               &                                + utau(ji+1,jj) * tmask(ji+1,jj,1) )
+            vtau(ji,jj) = 0.5 * ( 2. - vmask(ji,jj,1) ) * ( vtau(ji,jj) * tmask(ji,jj,1) &
+               &                                + vtau(ji,jj+1) * tmask(ji,jj+1,1) )
+          END DO
+         END DO
+
+         CALL lbc_lnk( utau(:,:), 'U', -1. )
+         CALL lbc_lnk( vtau(:,:), 'V', -1. )
+
+         ! for basin budget and cooerence
+         !--------------------------------------------------
+!CDIR COLLAPSE
+           emp (:,:) = evap(:,:) - sf(jp_prec)%fnow(:,:,1) * tmask(:,:,1)
+!CDIR COLLAPSE
+           emps(:,:) = emp(:,:)
+
+         CALL iom_put( "qlw_oce",   qbw  )                 ! output downward longwave heat over the ocean
+         CALL iom_put( "qsb_oce", - ha   )                 ! output downward sensible heat over the ocean
+         CALL iom_put( "qla_oce", - elat )                 ! output downward latent   heat over the ocean
+         CALL iom_put( "qns_oce",   qns  )                 ! output downward non solar heat over the ocean
+
+      ENDIF
+
+   END SUBROUTINE sbc_blk_ecmwf
+   
+   SUBROUTINE fluxes_ecmwf(alat,alon,hour,                             &
+        sst,tnow,shnow,unow,vnow,mslnow,cldnow,qsw,qbw,ha,elat,        &
+        evap,taux,tauy)
+      !!----------------------------------------------------------------------
+      !!                    ***  ROUTINE fluxes_ecmwf  ***
+      !!
+      !! --- it provides SURFACE HEAT and MOMENTUM FLUXES in MKS :
+      !!
+      !!  1)   Water flux (WFLUX)                 [ watt/m*m ]
+      !!  2)   Short wave flux (QSW)              [ watt/m*m ] Reed 1977
+      !!  3)   Long wave flux backward (QBW)      [ watt/m*m ]
+      !!  4)   Latent heat of evaporation (ELAT)  [ watt/m*m ]
+      !!  5)   Sensible heat flux   (HA)          [ watt/m*m ]
+      !!  6)   Wind stress x-component   (TAUX)   [ newton/m*m ]
+      !!  7)   Wind stress y-component   (TAUY)   [ newton/m*m ]
+      !!
+      !!----------------------------------------------------------------------
+      !!
+
+      USE sbcblk_core, ONLY: turb_core_2z ! For wave coupling and Tair/rh from 2 to 10m
+      USE wrk_nemo, ONLY:  wrk_in_use, wrk_not_released
+      USE wrk_nemo, ONLY:  rspeed     => wrk_2d_1
+      USE wrk_nemo, ONLY:  sh10now    => wrk_2d_2
+      USE wrk_nemo, ONLY:  t10now     => wrk_2d_3
+      USE wrk_nemo, ONLY:  cdx        => wrk_2d_4   ! --- drag coeff.
+      USE wrk_nemo, ONLY:  ce         => wrk_2d_5   ! --- turbulent exchange coefficients
+      USE wrk_nemo, ONLY:  shms       => wrk_2d_6
+      USE wrk_nemo, ONLY:  rhom       => wrk_2d_7
+      USE wrk_nemo, ONLY:  sstk       => wrk_2d_8
+      USE wrk_nemo, ONLY:  ch         => wrk_2d_10
+      USE wrk_nemo, ONLY:  rel_windu  => wrk_2d_11
+      USE wrk_nemo, ONLY:  rel_windv  => wrk_2d_12
+
+      REAL(wp), INTENT(in   ) :: hour
+      REAL(wp), INTENT(in   ), DIMENSION (:,:) :: sst, unow, alat , alon
+      REAL(wp), INTENT(in   ), DIMENSION (:,:) :: vnow, cldnow, mslnow
+      REAL(wp), INTENT(out  ), DIMENSION (:,:) :: qsw, qbw, ha, elat
+      REAL(wp), INTENT(out  ), DIMENSION (:,:) :: evap,taux,tauy
+      REAL(wp), INTENT(inout), DIMENSION (:,:) :: tnow , shnow
+
+      INTEGER :: ji,jj 
+      REAL(wp)  :: wair, vtnow, ea, deltemp, s, stp , fh , fe
+      REAL(wp)  :: esre, cseep
+
+      !!----------------------------------------------------------------------
+      !!     coefficients ( in MKS )  :
+      !!----------------------------------------------------------------------
+
+      REAL(wp), PARAMETER ::  ps = 1013.    ! --- surface air pressure
+      REAL(wp), PARAMETER ::  expsi=0.622   ! --- expsi
+      REAL(wp), PARAMETER ::  rd=287.       ! --- dry air gas constant
+      REAL(wp), PARAMETER ::  cp=1005.      ! --- specific heat capacity
+      REAL(wp), PARAMETER ::  onsea=0.98    ! --- specific humidity factors
+      REAL(wp), PARAMETER ::  par1=640380.  ! [Kg/m3]
+      REAL(wp), PARAMETER ::  par2=-5107.4  ! [K]
+
+      !---------------------------------------------------------------------
+      !--- define Kondo parameters
+      !---------------------------------------------------------------------
+
+      REAL(wp), DIMENSION(5) :: a_h = (/0.0,0.927,1.15,1.17,1.652/)
+      REAL(wp), DIMENSION(5) :: a_e = (/0.0,0.969,1.18,1.196,1.68/)
+      REAL(wp), DIMENSION(5) :: b_h = (/1.185,0.0546,0.01,0.0075,-0.017/)
+      REAL(wp), DIMENSION(5) :: b_e = (/1.23,0.0521,0.01,0.008,-0.016/)
+      REAL(wp), DIMENSION(5) :: c_h = (/0.0,0.0,0.0,-0.00045,0.0/)
+      REAL(wp), DIMENSION(5) :: c_e = (/0.0,0.0,0.0,-0.0004,0.0/)
+      REAL(wp), DIMENSION(5) :: p_h = (/-0.157,1.0,1.0,1.0,1.0/)
+      REAL(wp), DIMENSION(5) :: p_e = (/-0.16,1.0,1.0,1.0,1.0/)
+      INTEGER :: kku                        !index varing with wind speed
+
+      ! Set-up access to workspace arrays
+      IF( wrk_in_use(2, 1,2,3,4,5,6,7,8,10,11,12) ) THEN
+         CALL ctl_stop('blk_ecmwf: requested workspace arrays unavailable')   ;   RETURN
+      END IF
+
+      !!----------------------------------------------------------------------
+      !! ------------------ (i)      short wave
+      !!----------------------------------------------------------------------
+
+       call qshort(hour,alat,alon,cldnow,qsw)
+
+      DO jj=1,jpj
+       DO ji=1,jpi
+          rel_windu(ji,jj) = unow(ji,jj) - 0.5_wp * ( ssu_m(ji-1,jj) + ssu_m(ji,jj) )
+          rel_windv(ji,jj) = vnow(ji,jj) - 0.5_wp * ( ssv_m(ji,jj-1) + ssv_m(ji,jj) )
+       END DO
+      END DO
+
+          rspeed(:,:)= SQRT(rel_windu(:,:)*rel_windu(:,:)   &
+         &                   + rel_windv(:,:)*rel_windv(:,:)) 
+
+          sstk(:,:) = sst(:,:) + rtt                          !- SST data in Kelvin degrees
+          shms(:,:) = (1/1.22)*onsea*par1*EXP(par2/sstk(:,:)) !- Saturation Specific Humidity
+
+      ! --- Transport temperature and humidity from 2m to 10m
+      !----------------------------------------------------------------------
+
+      t10now(:,:) = 0.0           ;   sh10now(:,:)= 0.0
+      ! Note that air temp is converted in potential temp
+      CALL turb_core_2z(2.,10.,sstk,tnow+2*0.0098,shms,shnow,rspeed,        &
+         &              Cdx,Ch,Ce,t10now,sh10now )
+      tnow(:,:)  = t10now(:,:)    ;    shnow(:,:) = sh10now(:,:)
+
+      !!----------------------------------------------------------------------
+      !! ------------------ (ii)    net long wave
+      !!----------------------------------------------------------------------
+
+      DO jj=1,jpj
+       DO ji=1,jpi
+
+          wair = shnow(ji,jj) / (1 - shnow(ji,jj))    ! mixing ratio of the air (Wallace and Hobbs)
+          vtnow = (tnow(ji,jj)*(expsi+wair))/(expsi*(1.+wair))   ! virtual temperature of air
+          rhom(ji,jj) = 100.*(ps/rd)/vtnow                       ! density of the moist air
+          ea   = (wair  / (wair  + 0.622 )) * mslnow(ji,jj)
+
+          qbw(ji,jj) = emic*stefan*( sstk(ji,jj)**4. )                    &
+                - ( stefan*( tnow(ji,jj)**4. ) * ( 0.653 + 0.00535*ea ) )  &
+                 * ( 1. + 0.1762*( cldnow(ji,jj)**2. ) )
+
+       END DO
+      END DO
+
+      DO jj=1,jpj
+       DO ji=1,jpi
+
+      !!----------------------------------------------------------------------
+      !! ------------------ (iii)   sensible heat
+      !!----------------------------------------------------------------------
+
+      !! --- calculates the term :      ( Ts - Ta )
+      !!----------------------------------------------------------------------
+      deltemp = sstk(ji,jj) - tnow (ji,jj)
+
+      !!----------------------------------------------------------------------
+      !! --- variable turbulent exchange coefficients ( from Kondo 1975 )
+      !! --- calculate the Neutral Transfer Coefficent using an empiric formula
+      !! --- by Kondo et al. Then it applies the diabatic approximation.
+      !!----------------------------------------------------------------------
+
+             s=deltemp/(wndm(ji,jj)**2.)   !! --- calculate S 
+
+             stp=s*abs(s)/(abs(s)+0.01)     !! --- calculate the Stability Parameter 
+
+      !!----------------------------------------------------------------------
+      !! --- for stable condition (sst-t_air < 0):
+      !!----------------------------------------------------------------------
+
+        IF(s.lt.0. .and. ((stp.gt.-3.3).and.(stp.lt.0.))) then
+                  fh=0.1+0.03*stp+0.9*exp(4.8*stp)
+                  fe=fh
+        ELSE IF(s.lt.0. .and. stp.le.-3.3) THEN
+                  fh=0.
+                  fe=fh
+        ELSE                                       ! --- for unstable condition 
+                fh=1.0+0.63*sqrt(stp)
+                fe=fh
+        ENDIF
+
+      !!----------------------------------------------------------------------
+      !! --- calculate the coefficient CH,CE,CD
+      !!----------------------------------------------------------------------
+
+        IF      (wndm(ji,jj)>=0. .AND. wndm(ji,jj)<=2.2)  THEN
+              kku=1
+        ELSE IF (wndm(ji,jj)>2.2 .AND. wndm(ji,jj)<=5.0)  THEN
+              kku=2
+        ELSE IF (wndm(ji,jj)>5.0 .AND. wndm(ji,jj)<=8.0)  THEN
+              kku=3
+        ELSE IF (wndm(ji,jj)>8.0 .AND. wndm(ji,jj)<=25.0) THEN
+              kku=4
+        ELSE IF (wndm(ji,jj)>25.0 ) THEN
+              kku=5
+        ENDIF
+
+        ch(ji,jj) = ( a_h(kku) + b_h(kku) * wndm(ji,jj) ** p_h(kku)      &
+                    + c_h(kku) * (wndm(ji,jj)-8 ) **2) * fh
+
+        ce(ji,jj) = ( a_e(kku) + b_e(kku) * wndm(ji,jj) ** p_e(kku)      &
+                    + c_e(kku) * (wndm(ji,jj)-8 ) **2) * fe
+
+        ch(ji,jj) = ch(ji,jj) / 1000.0
+        ce(ji,jj) = ce(ji,jj) / 1000.0
+
+        IF (wndm(ji,jj)<0.3) THEN
+
+            ch(ji,jj) = 1.3e-03 * fh
+            ce(ji,jj) = 1.5e-03 * fe
+
+        ELSE IF(wndm(ji,jj)>50.0) THEN
+
+            ch(ji,jj) = 1.25e-03 * fh
+            ce(ji,jj) = 1.30e-03 * fe
+
+        ENDIF
+
+      !!----------------------------------------------------------------------
+      !! --- calculates the SENSIBLE HEAT FLUX in MKS ( watt/m*m )
+      !!----------------------------------------------------------------------
+
+      HA(ji,jj) = rhom(ji,jj)*cp*ch(ji,jj)*wndm(ji,jj)*deltemp
+
+      !!----------------------------------------------------------------------
+      !! ------------------ (iv)  latent heat
+      !! --- calculates the LATENT HEAT FLUX  ( watt/m*m )
+      !! --- ELAT = L*rho*Ce*|V|*[qs(Ts)-qa(t2d)]
+      !!----------------------------------------------------------------------
+
+      esre  = shms(ji,jj) - shnow(ji,jj)   ! --- calculates the term : qs(Ta)-qa(t2d)
+
+      cseep = ce(ji,jj) * wndm(ji,jj) * esre     ! --- calculates the term : Ce*|V|*[qs(Ts)-qa(t2d)]
+
+      evap(ji,jj) = (cseep * rhom(ji,jj))  ! in [kg/m2/sec] !! --- calculates the EVAPORATION RATE [m/yr]
+
+      elat(ji,jj) = rhom(ji,jj) * cseep * heatlat(sst(ji,jj))
+
+      !!----------------------------------------------------------------------
+      !! --- calculates the Drag Coefficient
+      !!----------------------------------------------------------------------
+
+      !!----------------------------------------------------------------------
+      !! --- deltemp should be (Ts - Ta) in the formula estimating
+      !! --- drag coefficient
+      !!----------------------------------------------------------------------
+
+        IF( .NOT. ln_cdgw ) THEN
+        cdx(ji,jj) = cd_HR(wndm(ji,jj),deltemp)
+        ENDIF
+
+       END DO
+      END DO
+
+      !!----------------------------------------------------------------------
+      !! --- calculates the wind stresses in MKS ( newton/m*m )
+      !! ---            taux= rho*Cd*|V|u      tauy= rho*Cd*|V|v
+      !!----------------------------------------------------------------------
+
+        taux(:,:)= rhom(:,:) * cdx(:,:) * rspeed(:,:) * rel_windu(:,:)
+        tauy(:,:)= rhom(:,:) * cdx(:,:) * rspeed(:,:) * rel_windv(:,:)
+
+
+      IF( wrk_not_released(2, 1,2,3,4,5,6,7,8,10,11,12) )    &
+         CALL ctl_stop('fluxes_ecmwf: failed to release workspace arrays')
+
+   END SUBROUTINE fluxes_ecmwf
+
+      REAL(wp) FUNCTION cd_HR(speed,delt)
+      !!----------------------------------------------------------------------
+      !! --- calculates the Drag Coefficient as a function of the abs. value
+      !! --- of the wind velocity ( Hellermann and Rosenstein )
+      !!----------------------------------------------------------------------
+
+       REAL(wp), INTENT(in) :: speed,delt
+       REAL(wp), PARAMETER  :: a1=0.934e-3 , a2=0.788e-4, a3=0.868e-4,     &
+        &                a4=-0.616e-6, a5=-.120e-5, a6=-.214e-5
+
+        cd_HR = a1 + a2*speed + a3*delt + a4*speed*speed        &
+           + a5*delt*delt  + a6*speed*delt
+
+      END FUNCTION cd_HR
+
+      REAL(wp) function HEATLAT(t)
+      !!----------------------------------------------------------------------
+      !! --- calculates the Latent Heat of Vaporization ( J/kg ) as function of
+      !! --- the temperature ( Celsius degrees )
+      !! --- ( from A. Gill  pag. 607 )
+      !!
+      !! --- Constant Latent Heat of Vaporization ( Rosati,Miyakoda 1988 )
+      !!     L = 2.501e+6  (MKS)
+      !!----------------------------------------------------------------------
+
+      REAL(wp) , intent(in) :: t
+
+      heatlat = 2.5008e+6 -2.3e+3*t
+
+      END FUNCTION HEATLAT
+
+   SUBROUTINE qshort(hour,alat,alon,cldnow,qsw)
+      !!----------------------------------------------------------------------
+      !!                    ***  ROUTINE qshort  ***
+      !!
+      !! ** Purpose :   Compute Solar Radiation
+      !!
+      !! ** Method  :   Compute Solar Radiation according Astronomical
+      !!                formulae
+      !!
+      !! References :   Reed RK (1975) and Reed RK (1977)
+      !!
+      !! Note: alat,alon - (lat, lon)  in radians
+      !!----------------------------------------------------------------------
+        REAL(wp), INTENT (in) :: hour
+
+        REAL(wp), INTENT(in ), DIMENSION(:,:) :: alat,alon
+        REAL(wp), INTENT(in ), DIMENSION(:,:) :: cldnow
+        REAL(wp), INTENT(out), DIMENSION(:,:) :: qsw
+        REAL(wp), DIMENSION(12) :: alpham
+
+        REAL(wp), PARAMETER ::   eclips=23.439* (3.141592653589793_wp / 180._wp)
+        REAL(wp), PARAMETER ::   solar = 1350.
+        REAL(wp), PARAMETER ::   tau = 0.7
+        REAL(wp), PARAMETER ::   aozone = 0.09
+        REAL(wp), PARAMETER ::   yrdays = 360.
+        REAL(wp) :: days, th0,th02,th03, sundec, thsun, coszen, qatten
+        REAL(wp) :: qzer, qdir,qdiff,qtot,tjul,sunbet
+        REAL(wp) :: albedo
+        INTEGER :: jj, ji
+
+      !!----------------------------------------------------------------------
+      !! --- albedo monthly values from Payne (1972) as means of the values
+      !! --- at 40N and 30N for the Atlantic Ocean ( hence the same latitudinal
+      !! --- band of the Mediterranean Sea ) :
+      !!----------------------------------------------------------------------
+
+        data alpham /0.095,0.08,0.065,0.065,0.06,0.06,0.06,0.06,        &
+                    0.065,0.075,0.09,0.10/
+
+      !!----------------------------------------------------------------------
+      !!   days is the number of days elapsed until the day=nday_year
+      !!----------------------------------------------------------------------
+        days = nday_year -1.
+        th0  = 2.*rpi*days/yrdays
+        th02 = 2.*th0
+        th03 = 3.*th0
+
+      !! --- sun declination :
+      !!----------------------------------------------------------------------
+        sundec = 0.006918 - 0.399912*cos(th0) + 0.070257*sin(th0) -   &
+                          0.006758*cos(th02) + 0.000907*sin(th02) -   &
+                          0.002697*cos(th03) + 0.001480*sin(th03)
+
+      DO jj=1,jpj
+       DO ji=1,jpi
+
+      !! --- sun hour angle :
+      !!----------------------------------------------------------------------
+        thsun = (hour -12.)*15.*rad + alon(ji,jj)
+
+      !! --- cosine of the solar zenith angle :
+      !!----------------------------------------------------------------------
+        coszen =sin(alat(ji,jj))*sin(sundec)                 &
+                    +cos(alat(ji,jj))*cos(sundec)*cos(thsun)
+
+        IF(coszen .LE. 5.035D-04) THEN
+          coszen = 0.0
+          qatten = 0.0
+        ELSE
+          qatten = tau**(1./coszen)
+        END IF
+
+        qzer  = coszen * solar *                                 &
+                (1.0+1.67E-2*cos(rpi*2.*(days-3.0)/365.0))**2
+        qdir  = qzer * qatten
+        qdiff = ((1.-aozone)*qzer - qdir) * 0.5
+        qtot  =  qdir + qdiff
+
+       tjul = (days -81.)*rad
+
+      !! --- sin of the solar noon altitude in radians :
+      !!----------------------------------------------------------------------
+       sunbet=sin(alat(ji,jj))*sin(eclips*sin(tjul)) +   &
+              cos(alat(ji,jj))*cos(eclips*sin(tjul))
+
+      !! --- solar noon altitude in degrees :
+      !!----------------------------------------------------------------------
+
+      sunbet = asin(sunbet)/rad
+
+      !!----------------------------------------------------------------------
+      !! --- calculates the albedo according to Payne (1972)
+      !!----------------------------------------------------------------------
+
+       albedo = alpham(nmonth)
+
+      !!----------------------------------------------------------------------
+      !! --- ( radiation as from Reed(1977), Simpson and Paulson(1979) )
+      !! --- calculates SHORT WAVE FLUX ( watt/m*m )
+      !! --- ( Rosati,Miyakoda 1988 ; eq. 3.8 )
+      !!----------------------------------------------------------------------
+
+       IF(cldnow(ji,jj).LT.0.3) THEN
+        qsw(ji,jj) = qtot * (1.-albedo)
+       ELSE
+        qsw(ji,jj) = qtot*(1.-0.62*cldnow(ji,jj)              &
+                                + .0019*sunbet)*(1.-albedo)
+       ENDIF
+
+       END DO
+      END DO
+
+   END SUBROUTINE qshort
+
+   !!======================================================================
+
+END MODULE sbcblk_ecmwf
Index: branches/2011/dev_r2855_INGV2_3_blk_wave/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90
===================================================================
--- branches/2011/dev_r2855_INGV2_3_blk_wave/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90	(revision 2989)
+++ branches/2011/dev_r2855_INGV2_3_blk_wave/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90	(revision 2990)
@@ -29,4 +29,5 @@
    USE sbcblk_clio      ! surface boundary condition: bulk formulation : CLIO
    USE sbcblk_core      ! surface boundary condition: bulk formulation : CORE
+   USE sbcblk_ecmwf     ! surface boundary condition: bulk formulation : MFS
    USE sbcice_if        ! surface boundary condition: ice-if sea-ice model
    USE sbcice_lim       ! surface boundary condition: LIM 3.0 sea-ice model
@@ -46,4 +47,5 @@
    USE in_out_manager   ! I/O manager
    USE lib_mpp          ! MPP library
+   USE sbcwave          ! Wave module
 
    IMPLICIT NONE
@@ -78,5 +80,5 @@
       !!
       NAMELIST/namsbc/ nn_fsbc   , ln_ana , ln_flx  , ln_blk_clio, ln_blk_core, ln_cpl,   &
-         &             ln_apr_dyn, nn_ice , ln_dm2dc, ln_rnf     , ln_ssr     , nn_fwb
+         &             ln_blk_ecmwf, ln_apr_dyn, nn_ice , ln_dm2dc, ln_rnf, ln_ssr     , nn_fwb, ln_cdgw
       !!----------------------------------------------------------------------
 
@@ -107,5 +109,6 @@
          WRITE(numout,*) '              flux       formulation                     ln_flx      = ', ln_flx
          WRITE(numout,*) '              CLIO bulk  formulation                     ln_blk_clio = ', ln_blk_clio
-         WRITE(numout,*) '              CLIO bulk  formulation                     ln_blk_core = ', ln_blk_core
+         WRITE(numout,*) '              CORE bulk  formulation                     ln_blk_core = ', ln_blk_core
+         WRITE(numout,*) '              MFS  bulk  formulation                     ln_blk_ecmwf= ', ln_blk_ecmwf
          WRITE(numout,*) '              coupled    formulation (T if key_sbc_cpl)  ln_cpl      = ', ln_cpl
          WRITE(numout,*) '           Misc. options of sbc : '
@@ -154,4 +157,8 @@
       IF( ln_dm2dc .AND. ( ( NINT(rday) / ( nn_fsbc * NINT(rdt) ) )  < 8 ) )   &
          &   CALL ctl_warn( 'diurnal cycle for qsr: the sampling of the diurnal cycle is too small...' )
+
+       !drag coefficient read from wave model definable only with ecmwf bulk formulae and core 
+       IF(ln_cdgw .AND. .NOT.(ln_blk_ecmwf .OR. ln_blk_core) )              &
+          &   CALL ctl_stop( 'drag coefficient read from wave model definable only with ecmwf bulk formulae and core')
       
       !                          ! Choice of the Surface Boudary Condition (set nsbc)
@@ -161,4 +168,5 @@
       IF( ln_blk_clio     ) THEN   ;   nsbc =  3   ; icpt = icpt + 1   ;   ENDIF       ! CLIO bulk       formulation
       IF( ln_blk_core     ) THEN   ;   nsbc =  4   ; icpt = icpt + 1   ;   ENDIF       ! CORE bulk       formulation
+      IF( ln_blk_ecmwf    ) THEN   ;   nsbc =  6   ; icpt = icpt + 1   ;   ENDIF       ! MFS  bulk       formulation
       IF( ln_cpl          ) THEN   ;   nsbc =  5   ; icpt = icpt + 1   ;   ENDIF       ! Coupled         formulation
       IF( cp_cfg == 'gyre') THEN   ;   nsbc =  0                       ;   ENDIF       ! GYRE analytical formulation
@@ -181,4 +189,5 @@
          IF( nsbc ==  4 )   WRITE(numout,*) '              CORE bulk formulation'
          IF( nsbc ==  5 )   WRITE(numout,*) '              coupled formulation'
+         IF( nsbc ==  6 )   WRITE(numout,*) '              MFS Bulk formulation'
       ENDIF
       !
@@ -228,4 +237,5 @@
       !                                                  ! averaged over nf_sbc time-step
 
+      IF (ln_cdgw) CALL sbc_wave( kt )
                                                    !==  sbc formulation  ==!
                                                             
@@ -238,4 +248,5 @@
       CASE(  4 )   ;   CALL sbc_blk_core( kt )                    ! bulk formulation : CORE for the ocean
       CASE(  5 )   ;   CALL sbc_cpl_rcv ( kt, nn_fsbc, nn_ice )   ! coupled formulation
+      CASE(  6 )   ;   CALL sbc_blk_ecmwf( kt )                   ! bulk formulation : MFS for the ocean
       CASE( -1 )                                
                        CALL sbc_ana     ( kt )                    ! ESOPA, test ALL the formulations
Index: branches/2011/dev_r2855_INGV2_3_blk_wave/NEMOGCM/NEMO/OPA_SRC/SBC/sbcwave.F90
===================================================================
--- branches/2011/dev_r2855_INGV2_3_blk_wave/NEMOGCM/NEMO/OPA_SRC/SBC/sbcwave.F90	(revision 2990)
+++ branches/2011/dev_r2855_INGV2_3_blk_wave/NEMOGCM/NEMO/OPA_SRC/SBC/sbcwave.F90	(revision 2990)
@@ -0,0 +1,88 @@
+MODULE sbcwave
+   !!======================================================================
+   !!                       ***  MODULE  sbcwave  ***
+   !! Wave module 
+   !!======================================================================
+   !! History :  3.3.1  !   2011-09  (Adani M)  Original code
+   !!----------------------------------------------------------------------
+   USE iom             ! I/O manager library
+   USE in_out_manager  ! I/O manager
+   USE lib_mpp         ! distribued memory computing library
+   USE fldread	       ! read input fields
+   USE sbc_oce	       ! Surface boundary condition: ocean fields
+
+   
+   !!----------------------------------------------------------------------
+   !!   sbc_wave       : read drag coefficient from wave model in netcdf files 
+   !!----------------------------------------------------------------------
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   sbc_wave    ! routine called in sbc_blk_core or sbc_blk_ecmwf
+   
+   TYPE(FLD), ALLOCATABLE, DIMENSION(:)  :: sf_wave	  ! structure of input fields (file informations, fields read)
+   REAL(wp),PUBLIC,ALLOCATABLE,DIMENSION (:,:)       :: cdn_wave 
+
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 4.0 , NEMO Consortium (2011) 
+   !! $Id: $
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE sbc_wave( kt )
+      !!---------------------------------------------------------------------
+      !!                     ***  ROUTINE sbc_apr  ***
+      !!
+      !! ** Purpose :   read drag coefficient from wave model  in netcdf files.
+      !!
+      !! ** Method  : - Read namelist namsbc_wave
+      !!              - Read Cd_n10 fields in netcdf files 
+      !! ** action  :   
+      !!               
+      !!---------------------------------------------------------------------
+      INTEGER, INTENT( in  ) ::  kt   ! ocean time step
+      INTEGER                ::  ierror   ! return error code
+      CHARACTER(len=100)     ::  cn_dir_cdg                       ! Root directory for location of drag coefficient files
+      TYPE(FLD_N)            ::  sn_cdg                          ! informations about the fields to be read
+      !!---------------------------------------------------------------------
+      NAMELIST/namsbc_wave/  sn_cdg, cn_dir_cdg
+      !!---------------------------------------------------------------------
+
+      !!----------------------------------------------------------------------
+      !
+      !
+      !                                         ! -------------------- !
+      IF( kt == nit000 ) THEN                   ! First call kt=nit000 !
+         !                                      ! -------------------- !
+         !                                            !* set file information (default values)
+         ! ... default values (NB: frequency positive => hours, negative => months)
+         !              !   file   ! frequency !  variable  ! time intep !  clim   ! 'yearly' or ! weights  ! rotation !
+         !              !   name   !  (hours)  !   name     !   (T/F)    !  (T/F)  !  'monthly'  ! filename ! pairs    !
+         sn_cdg = FLD_N('cdg_wave'  ,    1     ,'drag_coeff',  .true.    , .false. ,   'daily'   , ''       , ''       )
+         cn_dir_cdg = './'          ! directory in which the Patm data are 
+         
+
+         REWIND( numnam )                             !* read in namlist namsbc_wave
+         READ  ( numnam, namsbc_wave ) 
+         !
+
+         ALLOCATE( sf_wave(1), STAT=ierror )           !* allocate and fill sf_wave with sn_cdg
+         IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave: unable to allocate sf_wave structure' )
+         !
+         CALL fld_fill( sf_wave, (/ sn_cdg /), cn_dir_cdg, 'sbc_wave', 'Wave module ', 'namsbc_wave' )
+                                ALLOCATE( sf_wave(1)%fnow(jpi,jpj,1)   )
+         IF( sn_cdg%ln_tint )   ALLOCATE( sf_wave(1)%fdta(jpi,jpj,1,2) )
+         ALLOCATE( cdn_wave(jpi,jpj) )
+         cdn_wave(:,:) = 0.0
+      ENDIF
+         !
+         !
+      CALL fld_read( kt, nn_fsbc, sf_wave )      !* read drag coefficient from external forcing
+      cdn_wave(:,:) = sf_wave(1)%fnow(:,:,1)
+
+   END SUBROUTINE sbc_wave
+      
+   !!======================================================================
+END MODULE sbcwave
