Index: tags/nemo_dev_x3/NEMO/OPA_SRC/SBC/ocesbc_cpl_ice.h90
===================================================================
--- branches/nemo/NEMO/OPA_SRC/SBC/ocesbc_cpl_ice.h90	(revision 4)
+++ 	(revision )
@@ -1,111 +1,0 @@
-   !!----------------------------------------------------------------------
-   !!                       ***  ocesbc_cpl_ice.h90  ***
-   !!----------------------------------------------------------------------
-
-   SUBROUTINE oce_sbc(kt)
-      !!---------------------------------------------------------------------
-      !!                 ***  ROUTINE oce_sbc  ***
-      !!                    
-      !! ** Purpose :   Ocean surface boundaries conditions with 
-      !!        Louvain la Neuve Sea Ice Model in coupled mode
-      !!
-      !! History :
-      !!   1.0  !  00-10  (O. Marti)  Original code
-      !!   2.0  !  02-12  (G. Madec)  F90: Free form and module
-      !!----------------------------------------------------------------------
-      !! * Local declarations
-      INTEGER ::   ji, jj                   ! dummy loop indices
-      REAL(wp) ::   ztx, ztaux, zty, ztauy, ztrp, zsrp
-      REAL(wp) ::   ztdta, ztgel, zqrp
-
-      !!----------------------------------------------------------------------
-      !!  OPA 8.5, LODYC-IPSL (2002)
-      !!----------------------------------------------------------------------
- 
-      !!----------------------------------------------------------------------        
-      ! 1. initialization to zero at kt = nit000
-      ! ---------------------------------------
-      
-      IF( kt == nit000 ) THEN     
-         qsr   (:,:) = 0.e0
-         freeze(:,:) = 0.e0
-         qt    (:,:) = 0.e0
-         q     (:,:) = 0.e0
-         qrp   (:,:) = 0.e0
-         emp   (:,:) = 0.e0
-         emps  (:,:) = 0.e0
-         erp   (:,:) = 0.e0
-#if defined key_dynspg_fsc 
-         dmp   (:,:) = 0.e0
-#endif
-      ENDIF
-
-      IF( MOD( kt-1, nfice ) == 0 ) THEN 
-
-         CALL oce_sbc_dmp   ! Computation of internal and evaporation damping terms       
-
-         ! Surface heat flux (W/m2)
-         ! -----------------------
-
-         ! restoring heat flux
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               ztgel = fzptn(ji,jj)
-#if defined key_dtasst
-               ztdta = MAX( sst(ji,jj),    ztgel )
-#else
-               ztdta = MAX( t_dta(ji,jj,1), ztgel )
-#endif
-               zqrp = ztrp * ( tb(ji,jj,1) - ztdta )
-
-               qrp(ji,jj) = (1.0-freeze(ji,jj) ) * zqrp
-            END DO
-         END DO
-
-         ! non solar heat flux + solar flux + restoring
-         q  (:,:) = fnsolar(:,:) + fsolar(:,:) + qrp(:,:)
-         qt (:,:) = q(:,:)
-
-         ! solar flux
-         qsr(:,:) = fsolar(:,:)
-
-#if defined key_dynspg_fsc  
-         ! total concentration/dilution effect (use on SSS)
-         emps(:,:) = fmass(:,:) + fsalt(:,:) + runoff(:,:) + erp(:,:)
-
-         ! total volume flux (use on sea-surface height)
-         emp (:,:) = fmass(:,:)  -  dmp(:,:) + runoff(:,:) + erp(:,:)
-
-#else
-         ! Rigid-lid (emp=emps=E-P-R+Erp) 
-         ! freshwater flux
-         emps(:,:) = fmass(:,:) + fsalt(:,:) + runoff(:,:) + erp(:,:)
-         emp (:,:) = emps(:,:)
-
-#endif
-
-         DO jj = 1, jpjm1
-            DO ji = 1, fs_jpim1   ! vertor opt.
-               ztx   = 0.5 * ( freeze(ji+1,jj) + freeze(ji+1,jj+1) )
-               ztaux = 0.5 * ( ftaux (ji+1,jj) + ftaux (ji+1,jj+1) )
-               taux(ji,jj) = (1.0-ztx) * taux(ji,jj) + ztx * ztaux
-
-               zty   = 0.5 * ( freeze(ji,jj+1) + freeze(ji+1,jj+1) )
-               ztauy = 0.5 * ( ftauy (ji,jj+1) + ftauy (ji+1,jj+1) )
-               tauy(ji,jj) = (1.0-zty) * tauy(ji,jj) + zty * ztauy
-            END DO
-         END DO
-         CALL lbc_lnk( taux, 'U', -1. )
-         CALL lbc_lnk( tauy, 'V', -1. )    
-
-         ! Re-initialization of fluxes
-         sst_io(:,:) = 0.0
-         sss_io(:,:) = 0.0
-         u_io  (:,:) = 0.0
-         v_io  (:,:) = 0.0
-         gtaux (:,:) = 0.
-         gtauy (:,:) = 0.
-
-      ENDIF
-
-   END SUBROUTINE oce_sbc
Index: tags/nemo_dev_x3/NEMO/OPA_SRC/SBC/ocesbc_cpl_noice.h90
===================================================================
--- branches/nemo/NEMO/OPA_SRC/SBC/ocesbc_cpl_noice.h90	(revision 4)
+++ 	(revision )
@@ -1,90 +1,0 @@
-   !!----------------------------------------------------------------------
-   !!                       ***  ocesbc_cpl_noice.h90  ***
-   !!----------------------------------------------------------------------
-
-   SUBROUTINE oce_sbc(kt)
-      !!---------------------------------------------------------------------
-      !!                 ***  ROUTINE oce_sbc  ***
-      !!                    
-      !! ** Purpose :   Ocean surface boundaries conditions in 
-      !!                coupled ocean/atmosphere case without sea-ice
-      !!
-      !! History :
-      !!   1.0  !  00-10  (O. Marti)  Original code
-      !!   2.0  !  02-12  (G. Madec)  F90: Free form and module
-      !!----------------------------------------------------------------------
-      !! * Modules used
-      USE cpl_oce 
-
-      !! * arguments
-      INTEGER, INTENT( in  ) ::   kt ! ocean time step
-
-      !! * Local declarations
-      INTEGER  ::   ji, jj, jf
-
-      REAL(wp) ::   ztrp
-      REAL(wp) ::   ztgel,zice,zhemis,zqrp,zqri,zq,zqi,zerp,ze,zei,zro
-
-      ! Compute fluxes
-      ! --------------
-
-      ! constant initialization
-      ztrp = -40.   ! restoring term for temperature (w/m2/k)
-
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-
-            ztgel = fzptn(ji,jj)   ! local freezing temperature
-
-            ! opa model ice freeze()
-
-            zice = tmask(ji,jj,1)
-            IF( tn(ji,jj,1) >=  ztgel ) zice = 0.
-            freeze(ji,jj) = zice
-
-            ! hemisphere indicator (=1 north, =-1 south)
-            
-            zhemis = float(isign(1, mjg(jj)-(jpjglo/2+1)))
-            
-            ! a) net downward radiative flux qsr()
-            ! - AGCM qsrc if no ice
-            ! - zero under ice (zice=1)
-
-            qsr(ji,jj) = (1.-zice)*qsrc(ji,jj)*tmask(ji,jj,1)
-
-            ! b) heat flux damping term qrp()
-            ! - no damping if no  ice      (zice=0) 
-            ! - gamma*min(0,t-tgel) if ice (zice=1)
-
-            zqrp = 0.
-            zqri = ztrp*MIN( 0., tb(ji,jj,1)-ztgel )
-            qrp(ji,jj) = ( ( 1. - zice ) * zqrp + zice * zqri ) * tmask(ji,jj,1)
-
-
-            ! c) net downward heat flux q() = q0 + qrp()
-            ! for q0
-            ! - AGCM qc if no  ice (zice=0)
-            ! - -2 watt/m2 (arctic) or -4 watt/m2 (antarctic) if ice (zice=1)
-            zq  = qc(ji,jj)
-            zqi = -3. + zhemis
-            qt(ji,jj) = ( (1.-zice) * zq + zice * zqi ) * tmask(ji,jj,1) + qrp(ji,jj)
-            q (ji,jj) = qt(ji,jj)
-            
-            ! d) water flux damping term erp()
-            ! - no damping
-            zerp = 0.
-            erp(ji,jj) = zerp
-            
-            ! e) net upward water flux e() = eo + runoff() + erp()
-            ! for e0
-            ! - AGCM if no ice (zice=0)
-            ! - 1.mm/day if climatological and opa ice (zice=1)
-            ze  = ec(ji,jj)
-            zei = 1./rday
-            zro = runoff(ji,jj)
-            emp(ji,jj) = ( ( 1. - zice ) *  ze + zice * zei + zro ) * tmask(ji,jj,1) + erp(ji,jj)
-            
-         END DO
-      END DO
-
-   END SUBROUTINE oce_sbc
Index: tags/nemo_dev_x3/NEMO/OPA_SRC/SBC/ocesbc_forced_ice.h90
===================================================================
--- branches/nemo/NEMO/OPA_SRC/SBC/ocesbc_forced_ice.h90	(revision 4)
+++ 	(revision )
@@ -1,101 +1,0 @@
-   !!----------------------------------------------------------------------
-   !!                     ***  ocesbc_forced_ice.h90  ***
-   !!----------------------------------------------------------------------
-
-   SUBROUTINE oce_sbc (kt)
-      !!---------------------------------------------------------------------
-      !!                   ***    ROUTINE oce_sbc ***
-      !!                    
-      !! ** Purpose : Ocean surface boundary conditions with 
-      !!        Louvain la Neuve Sea Ice Model in forced mode
-      !!        using 'clio' bulk formulea
-      !!
-      !! History :
-      !!   1.0  !  99-11  (M. Imbard)  Original code
-      !!        !  01-03  (D. Ludicone, E. Durand, G. Madec) free surf.
-      !!   2.0  !  02-09  (G. Madec, C. Ethe)  F90: Free form and module
-      !!----------------------------------------------------------------------
-      !! * arguments
-      INTEGER, INTENT( in  ) ::   kt   ! ocean time step
-
-      !! * Local declarations
-      INTEGER ::   ji, jj                   ! dummy loop indices
-      REAL ztx, ztaux, zty, ztauy
-
-      !!----------------------------------------------------------------------
-      !!  OPA 8.5, LODYC-IPSL (2002)
-      !!----------------------------------------------------------------------
- 
-
-      ! 1. initialization to zero at kt = nit000
-      ! ---------------------------------------
-      
-      IF( kt == nit000 ) THEN     
-         qsr    (:,:) = 0.e0
-         qt     (:,:) = 0.e0
-         q      (:,:) = 0.e0
-         qrp    (:,:) = 0.e0
-         emp    (:,:) = 0.e0
-         emps   (:,:) = 0.e0
-         erp    (:,:) = 0.e0
-#if defined key_dynspg_fsc 
-         dmp    (:,:) = 0.e0
-#endif
-      ENDIF
-
-      IF( MOD(kt-1,nfice) == 0 ) THEN
-
-         ! Computation of internal and evaporation damping terms       
-         CALL oce_sbc_dmp
-
-         ! Surface Ocean fluxes
-         ! ====================
-
-         ! Surface heat flux (W/m2)
-         ! -----------------
-
-         q   (:,:) = fnsolar(:,:) + fsolar(:,:)     ! non solar heat flux + solar flux
-         qt  (:,:) = q(:,:)
-         qsr (:,:) = fsolar(:,:)                     ! solar flux
-
-#if defined key_dynspg_fsc     
-         ! total concentration/dilution effect (use on SSS)
-         emps(:,:) = fmass(:,:) + fsalt(:,:) + runoff(:,:) + erp(:,:) + empold
-
-         ! total volume flux (use on sea-surface height)
-         emp (:,:) = fmass(:,:) -   dmp(:,:) + runoff(:,:) + erp(:,:) + empold      
-#else
-         ! Rigid-lid (emp=emps=E-P-R+Erp)
-         emps(:,:) = fmass(:,:) + fsalt(:,:) + runoff(:,:) + erp(:,:)     ! freshwater flux
-         emp (:,:) = emps(:,:)
-
-#endif
-
-         ! Surface stress
-         ! --------------
-
-         ! update the stress beloww sea-ice area
-         DO jj = 1, jpjm1
-            DO ji = 1, fs_jpim1   ! vertor opt.
-               ztx         = MAX( freezn(ji,jj), freezn(ji,jj+1) )   ! ice/ocean indicator at U- and V-points
-               zty         = MAX( freezn(ji,jj), freezn(ji+1,jj) )
-               ztaux       = 0.5 *( ftaux(ji+1,jj) + ftaux(ji+1,jj+1) ) ! ice-ocean stress at U- and V-points
-               ztauy       = 0.5 *( ftauy(ji,jj+1) + ftauy(ji+1,jj+1) )
-               taux(ji,jj) = (1.-ztx) * taux(ji,jj) + ztx * ztaux    ! stress at the ocean surface
-               tauy(ji,jj) = (1.-zty) * tauy(ji,jj) + zty * ztauy
-            END DO
-         END DO
-
-         ! boundary condition on the stress (taux,tauy)
-         CALL lbc_lnk( taux, 'U', -1. )
-         CALL lbc_lnk( tauy, 'V', -1. )
-
-         ! Re-initialization of fluxes
-         sst_io(:,:) = 0.0
-         sss_io(:,:) = 0.0
-         u_io  (:,:) = 0.0
-         v_io  (:,:) = 0.0
-
-      ENDIF
-
-   END SUBROUTINE oce_sbc
Index: tags/nemo_dev_x3/NEMO/OPA_SRC/SBC/ocesbc_forced_noice.h90
===================================================================
--- branches/nemo/NEMO/OPA_SRC/SBC/ocesbc_forced_noice.h90	(revision 4)
+++ 	(revision )
@@ -1,185 +1,0 @@
-   !!----------------------------------------------------------------------
-   !!                     ***  ocesbc_forced_noice.h90  ***
-   !!----------------------------------------------------------------------
-
-   SUBROUTINE oce_sbc(kt)
-      !!---------------------------------------------------------------------
-      !!                   ***    ROUTINE oce_sbc ***
-      !!                    
-      !! ** Purpose : Ocean surface boundary conditions 
-      !!        in forced mode using flux formulation or bulk formulation
-      !!       
-      !!
-      !! History :
-      !!   1.0  !  99-11  (M. Imbard)  Original code
-      !!        !  01-03  (D. Ludicone, E. Durand, G. Madec) free surf.
-      !!   2.0  !  02-09  (G. Madec, C. Ethe)  F90: Free form and module
-      !!----------------------------------------------------------------------
-      !! * Modules used
-      USE daymod
-#if ! defined key_dtasst
-      USE dtasst, ONLY : rclice
-#endif
-
-#if defined key_flx_bulk_monthly || defined key_flx_bulk_daily
-      USE blk_oce
-#endif
-
-#if defined key_flx_forced_daily
-      USE flx_oce
-#endif
-      !! * arguments
-      INTEGER, INTENT( in  ) ::   kt   ! ocean time step
-
-      !! * local declarations
-      INTEGER ::   ji, jj        ! dummy loop arguments
-      INTEGER ::   i15, ifreq             !      
-      REAL(wp) ::  zxy
-      REAL(wp) ::  zsice, zqri, zqrp, ztdta, zqrj
-      REAL(wp) ::  zq, zqi, zhemis, ztrp
-      REAL(wp), DIMENSION(jpi,jpj) :: zeri, zerps, ziclim
-      REAL(wp), DIMENSION(jpi,jpj) :: zqt, zqsr, zemp  
-
-      !!----------------------------------------------------------------------
-      !!  OPA 8.5, LODYC-IPSL (2002)
-      !!----------------------------------------------------------------------
- 
-      ! 1. initialization to zero at kt = nit000
-      ! ---------------------------------------
-      
-      IF( kt == nit000 ) THEN     
-         qsr    (:,:) = 0.e0
-         freeze (:,:) = 0.e0
-         qt     (:,:) = 0.e0
-         q      (:,:) = 0.e0
-         qrp    (:,:) = 0.e0
-         emp    (:,:) = 0.e0
-         emps   (:,:) = 0.e0
-         erp    (:,:) = 0.e0
-#if defined key_dynspg_fsc 
-         dmp    (:,:) = 0.e0
-#endif
-      ENDIF
-
-#if defined key_flx_bulk_monthly || defined key_flx_bulk_daily
-      ifreq      = nfbulk
-      zqt (:,:)  = qsr_oce(:,:) + qnsr_oce(:,:)
-      zqsr(:,:)  = qsr_oce(:,:)
-      zemp(:,:)  = evap(:,:) - tprecip(:,:)
-#endif 
-
-#if defined key_flx_forced_daily
-      ifreq      = 1
-      zqt (:,:)  = p_qt (:,:)
-      zqsr(:,:)  = p_qsr(:,:)
-      zemp(:,:)  = p_emp(:,:)
-#endif 
-
-      IF( MOD( kt-1, ifreq) == 0 ) THEN
-         ! Computation of internal and evaporation damping terms       
-         CALL oce_sbc_dmp
-
-         ztrp = -40.             ! restoring terme for temperature (w/m2/k)   
-         zsice = - 0.04 / 0.8    ! ratio of isohaline compressibility over isotherme compressibility 
-                                 ! ( d rho / dt ) / ( d rho / ds )      ( s = 34, t = -1.8 )
-         ! Flux computation
-         DO jj = 1, jpj
-            DO ji = 1, jpi      
-               ! climatological ice 
-#if defined key_dtasst 
-               ziclim(ji,jj) = FLOAT( NINT( rclice(ji,jj,1) ) )
-#else
-               ! tested only with key key_dtasst (A. Lazar 07/2001)
-               ! this loop in CASE of interpolation of monthly rclice
-               i15           = INT( 2.* FLOAT(nday) / (FLOAT( nobis(nmonth) ) + 0.5) )
-               zxy           = FLOAT(nday) / FLOAT(nobis(nmonth)) + 0.5 - i15
-               ziclim(ji,jj) = FLOAT( NINT( (1-zxy) * rclice(ji,jj,1) + zxy  * rclice(ji,jj,2) ) )
-#endif
-
-               ! avoid surfreezing point            
-               tn(ji,jj,1) = MAX( tn(ji,jj,1), fzptn(ji,jj) )
-
-               ! hemisphere indicator (=1 north, =-1 south)           
-               zhemis = FLOAT( isign(1, mjg(jj) - (jpjdta/2+1) ) )
-
-               ! restoring temperature (ztdta >= to local freezing temperature)            
-#if defined key_dtasst
-               ztdta = MAX( sst(ji,jj),    fzptn(ji,jj) )
-#else
-               ztdta = MAX( t_dta(ji,jj,1), fzptn(ji,jj) )
-#endif
-
-               ! a) net downward radiative flux qsr()           
-               qsr(ji,jj) = zqsr(ji,jj) * tmask(ji,jj,1)
-
-               ! b) heat flux damping term qrp()
-               ! - gamma*(t-tlevitus) if no  climatological ice (ziclim=0)
-               ! - gamma*(t-(tgel-1.))  if climatological ice and no opa ice   (ziclim=1 zicopa=0)
-               ! - gamma*min(0,t-tgel) if climatological and opa ice (ziclim=1 zicopa=1)
-
-               zqrp = ztrp * ( tb(ji,jj,1) - ztdta )
-               zqri = ztrp * ( tb(ji,jj,1) - ( fzptn(ji,jj) - 1.) )
-               zqrj = ztrp * MIN( 0., tb(ji,jj,1) - fzptn(ji,jj) )
-               qrp(ji,jj) = ( (1. - ziclim(ji,jj)) * zqrp   &
-                  + ziclim(ji,jj)  * ( ( 1 - freeze(ji,jj)) * zqri   &
-                  + freeze(ji,jj)  * zqrj ) ) * tmask(ji,jj,1)
-
-               ! c) net downward heat flux q() = q0 + qrp()
-               ! for q0
-               ! - ECMWF fluxes if no climatological ice      (ziclim=0)
-               ! - qrp if climatological ice and no opa ice   (ziclim=1 zicopa=0)
-               ! - -2 watt/m2 (arctic) or -4 watt/m2 (antarctic) if climatological and opa ice 
-               !                                              (ziclim=1 zicopa=1)
-               zq  = zqt(ji,jj)
-               zqi = -3. + zhemis
-               qt (ji,jj) = ( (1.-ziclim(ji,jj)) * zq   &
-                  +ziclim(ji,jj)  * freeze(ji,jj) * zqi )   &
-                  * tmask(ji,jj,1)   &
-                  + qrp(ji,jj)
-               q  (ji,jj) = qt (ji,jj)
-
-            END DO
-         END DO
-
-#if defined key_dynspg_fsc 
-         ! Free-surface
-
-         ! Water flux for zero buoyancy flux if no opa ice and ice clim
-         zeri(:,:) = -zsice * qrp(:,:) * ro0cpr * rauw / 34.0
-         zerps(:,:) = ziclim(:,:) * ( (1-freeze(:,:)) * zeri(:,:) )
-
-         ! Contribution to sea level:
-         ! net upward water flux emp() = e-p + runoff() + erp() + dmp + empold
-         emp (:,:) = zemp(:,:)     &   ! e-p data
-                   + runoff(:,:)   &   ! runoff data
-                   + erp(:,:)      &   ! restoring term to SSS data
-                   + dmp(:,:)      &   ! freshwater flux associated with internal damping
-                   + empold            ! domain averaged annual mean correction
-
-         ! Contribution to salinity:
-         ! net upward water flux emps() = e-p + runoff() + erp() + zerps + empold
-         emps(:,:) = zemp(:,:)     &
-                   + runoff(:,:)   &
-                   + erp(:,:)      &
-                   + zerps(:,:)    &
-                   + empold 
-#else
-         ! Rigid-lid (emp=emps=E-P-R+Erp)
-         ! freshwater flux
-         zeri(:,:)  = -zsice * qrp(:,:) * ro0cpr * rauw / 34.0
-         zerps(:,:) = ziclim(:,:) * ( (1-freeze(:,:)) * zeri(:,:) )
-         emps (:,:) = zemp(:,:)     &
-                    + runoff(:,:)   &
-                    + erp(:,:)      &
-                    + zerps(:,:)
-         emp (:,:) = emps(:,:)
-#endif  
-
-
-         ! Boundary condition on emp for free surface option
-         ! -------------------------------------------------
-         CALL lbc_lnk( emp, 'T', 1. )
-      
-      ENDIF
-
-   END SUBROUTINE oce_sbc
Index: tags/nemo_dev_x3/NEMO/OPA_SRC/SBC/blk_oce.F90
===================================================================
--- branches/nemo/NEMO/OPA_SRC/SBC/blk_oce.F90	(revision 4)
+++ tags/nemo_dev_x3/NEMO/OPA_SRC/SBC/blk_oce.F90	(revision 105)
@@ -17,5 +17,6 @@
 
    IMPLICIT NONE
-   LOGICAL, PUBLIC ::   l_bulk = .TRUE.   
+
+   LOGICAL, PUBLIC ::   l_bulk = .TRUE.   !: 
    
    !!----------------------------------------------------------------------
@@ -23,27 +24,27 @@
    !!----------------------------------------------------------------------
 
-   REAL(wp), DIMENSION(jpi,jpj) ::   & 
-      watm     ,      &   ! precipitation
-      tatm     ,      &   ! atmospheric temperature
-      hatm     ,      &   ! relative humidity
-      vatm     ,      &   ! wind speed
-      catm                ! percent of cloud cover
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   &  !:
+      watm     ,      &  !: precipitation
+      tatm     ,      &  !: atmospheric temperature
+      hatm     ,      &  !: relative humidity
+      vatm     ,      &  !: wind speed
+      catm               !: percent of cloud cover
 
-   REAL(wp), DIMENSION(jpi,jpj)    ::   &
-      gsst     ,      &  ! SST mean on nfbulk ocean time step
-      gsss               ! SSS mean on nfbulk ocean time step
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj)    ::   &  !:
+      gsst     ,      &  !: SST mean on nfbulk ocean time step
+      gsss               !: SSS mean on nfbulk ocean time step
 
    REAL(wp) ::        &
-      yearday  ,      &  ! number of days per year
-      rdtbs2             ! bulk time step divide by 2
+      yearday  ,      &  !: number of days per year
+      rdtbs2             !: bulk time step divide by 2
 #else
    !!----------------------------------------------------------------------
    !!   Default option                                         Empty module
    !!----------------------------------------------------------------------
-   LOGICAL, PUBLIC ::   l_bulk = .FALSE.  
+   LOGICAL, PUBLIC ::   l_bulk = .FALSE.  !:
 #endif
    
-   INTEGER ::         & !!! namdom : space/time domain (namlist)
-      nfbulk =  5         ! bulk computation frequency 
+   INTEGER ::         & !!: namdom : space/time domain (namlist)
+      nfbulk =  5        !: bulk computation frequency 
    !!----------------------------------------------------------------------
 END MODULE blk_oce
Index: tags/nemo_dev_x3/NEMO/OPA_SRC/SBC/bulk.F90
===================================================================
--- branches/nemo/NEMO/OPA_SRC/SBC/bulk.F90	(revision 4)
+++ tags/nemo_dev_x3/NEMO/OPA_SRC/SBC/bulk.F90	(revision 105)
@@ -3,5 +3,8 @@
    !!                           ***  bulk  ***
    !!======================================================================
-
+#if defined key_flx_bulk_monthly || defined key_flx_bulk_daily
+   !!----------------------------------------------------------------------
+   !!   'key_flx_bulk_monthly'                        monthly bulk formulea
+   !!   'key_flx_bulk_daily'                          monthly bulk formulea
    !!----------------------------------------------------------------------
    !!   bulk          : computation of fluxes using bulk formulation
@@ -14,4 +17,5 @@
    USE flxblk          ! bulk formulae
    USE blk_oce         ! bulk variable 
+   USE flx_oce
    USE taumod
    USE phycst          ! physical constants
@@ -91,5 +95,5 @@
 
 # if ! defined key_ice_lim
-         IF( l_ctl .AND. lwp ) THEN         ! print mean trends (used for debugging)
+         IF(l_ctl) THEN         ! print mean trends (used for debugging)
             WRITE(numout,*) ' Forcings '
             WRITE(numout,*) ' qsr_oce  : ', SUM( qsr_oce (:,:) * tmask(:,:,1) )
@@ -110,4 +114,14 @@
  
    END SUBROUTINE blk
+
+#else
+   !!----------------------------------------------------------------------
+   !!   Dummy module :                                     NO bulk formulea
+   !!----------------------------------------------------------------------
+CONTAINS
+   SUBROUTINE blk( kt )          ! Dummy routine
+      WRITE(*,*) 'blk: You should not see this print! error? ', kt
+   END SUBROUTINE blk
+#endif
  
    !!======================================================================
Index: tags/nemo_dev_x3/NEMO/OPA_SRC/SBC/flx_bulk_daily.h90
===================================================================
--- branches/nemo/NEMO/OPA_SRC/SBC/flx_bulk_daily.h90	(revision 4)
+++ tags/nemo_dev_x3/NEMO/OPA_SRC/SBC/flx_bulk_daily.h90	(revision 105)
@@ -59,4 +59,5 @@
       USE blk_oce         ! bulk variable
       USE bulk            ! bulk module
+
       !! * arguments
       INTEGER, INTENT( in  ) ::   kt ! ocean time step
@@ -64,20 +65,13 @@
       !! * Local declarations      
       INTEGER , PARAMETER :: jpday  = 365, jpmois = 12
-      INTEGER                    :: ji,jj
-      INTEGER                    :: i15,iday, idy
-      INTEGER , DIMENSION(jpday) :: istep_n
-      INTEGER , DIMENSION(jpmois):: istep_c, istep_x
-      INTEGER                    :: itime
-      REAL(wp) ::  zdate0
-      REAL(wp), DIMENSION(jpi,jpj) ::   &
-         zsst,       &  ! ???
-         zlon    ,   &  ! ???
-         zlat           ! ???
-
-      REAL(wp), DIMENSION(jpk) ::   &
-         zlev           ! ???
-      REAL(wp) ::   zxy
-      INTEGER  :: ipi,ipj,ipk
-      INTEGER  :: iman,imois,imois2
+      INTEGER  ::   i15,iday, idy
+      INTEGER  ::   ipi,ipj,ipk
+      INTEGER  ::   iman,imois,imois2
+      INTEGER, DIMENSION(jpday) :: istep_n
+      INTEGER, DIMENSION(jpmois):: istep_c, istep_x
+      INTEGER  ::   itime
+      REAL(wp) ::   zdate0, zxy
+      REAL(wp), DIMENSION(jpi,jpj) ::   zlon, zlat   ! ???
+      REAL(wp), DIMENSION(jpk) ::   zlev           ! ???
       CHARACTER(len=45)  ::  &
          clname_n = 'NCEP_1d_0101_0112_tair.nc',   &
Index: tags/nemo_dev_x3/NEMO/OPA_SRC/SBC/flx_forced_daily.h90
===================================================================
--- branches/nemo/NEMO/OPA_SRC/SBC/flx_forced_daily.h90	(revision 4)
+++ tags/nemo_dev_x3/NEMO/OPA_SRC/SBC/flx_forced_daily.h90	(revision 105)
@@ -70,8 +70,9 @@
 
 
-#if defined key_mpp
-      == >> forced compile error:
-      flx_forced_daily : Reading Netcdf file in mpp Not yet implemented
-#endif
+      IF( lk_mpp ) THEN
+         IF(lwp) WRITE(numout,cform_err)
+         IF(lwp) WRITE(numout,*) 'flx_forced_daily:  netcdf & mpp not yet possible'
+         nstop = nstop + 1
+      ENDIF
 
       ! Initialization
Index: tags/nemo_dev_x3/NEMO/OPA_SRC/SBC/flx_oce.F90
===================================================================
--- branches/nemo/NEMO/OPA_SRC/SBC/flx_oce.F90	(revision 4)
+++ tags/nemo_dev_x3/NEMO/OPA_SRC/SBC/flx_oce.F90	(revision 105)
@@ -13,4 +13,6 @@
 
    IMPLICIT NONE
+   PRIVATE
+   
 
    !!----------------------------------------------------------------------
@@ -21,29 +23,29 @@
    !! 'key_flx_forced_daily'
    !!----------------------------------------------------------------------
-   REAL(wp), DIMENSION(jpi,jpj,3) ::   &
-      p_qt ,        &   ! total heat flux ( solar + non solar)
-      p_qsr,        &   ! solar heat flux
-      p_emp             ! evaporation minus precipitation            
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   &  !:
+      p_qt ,        &   !: total heat flux ( solar + non solar)
+      p_qsr,        &   !: solar heat flux
+      p_emp             !: evaporation minus precipitation            
 
 #elif defined key_ice_lim || defined key_flx_bulk_monthly || defined key_flx_bulk_daily
-   REAL(wp), DIMENSION(jpi,jpj)    ::   &
-      qsr_ice  ,      &  ! solar flux over ice
-      qsr_oce  ,      &  ! solar flux over ocean
-      qnsr_oce ,      &  ! total non solar heat flux (Longwave downward radiation) over ocean 
-      qnsr_ice ,      &  ! total non solar heat flux (Longwave downward radiation) over ice
-      tprecip  ,      &  ! total precipitation ( or liquid precip minus evaporation in coupled mode)
-      sprecip  ,      &  ! solid (snow) precipitation
-      dqns_ice ,      &  ! total non solar sensibility over ice (LW+SEN+LA)
-      tn_ice   ,      &  ! ice surface temperature
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj)    ::   &  !:
+      qsr_ice  ,      &  !: solar flux over ice
+      qsr_oce  ,      &  !: solar flux over ocean
+      qnsr_oce ,      &  !: total non solar heat flux (Longwave downward radiation) over ocean 
+      qnsr_ice ,      &  !: total non solar heat flux (Longwave downward radiation) over ice
+      tprecip  ,      &  !: total precipitation ( or liquid precip minus evaporation in coupled mode)
+      sprecip  ,      &  !: solid (snow) precipitation
+      dqns_ice ,      &  !: total non solar sensibility over ice (LW+SEN+LA)
+      tn_ice   ,      &  !: ice surface temperature
 #if ! defined key_coupled
-      fr1_i0   ,      &  ! 1st part of the fraction of sol. rad.  which penetrate inside the ice cover
-      fr2_i0   ,      &  ! 2nd part of the fraction of sol. rad.  which penetrate inside the ice cover 
-      qla_ice  ,      &  ! latent flux over ice  
-      dqla_ice ,      &  ! latent sensibility over ice
-      evap               ! evaporation over ocean
+      fr1_i0   ,      &  !: 1st part of the fraction of sol. rad.  which penetrate inside the ice cover
+      fr2_i0   ,      &  !: 2nd part of the fraction of sol. rad.  which penetrate inside the ice cover 
+      qla_ice  ,      &  !: latent flux over ice  
+      dqla_ice ,      &  !: latent sensibility over ice
+      evap               !: evaporation over ocean
 #else
-      rrunoff  ,      &  ! runoff
-      calving  ,         ! calving
-      alb_ice            ! albedo of ice      
+      rrunoff  ,      &  !: runoff
+      calving  ,      &  !: calving
+      alb_ice            !: albedo of ice      
 #endif
 
Index: tags/nemo_dev_x3/NEMO/OPA_SRC/SBC/flxblk.F90
===================================================================
--- branches/nemo/NEMO/OPA_SRC/SBC/flxblk.F90	(revision 4)
+++ tags/nemo_dev_x3/NEMO/OPA_SRC/SBC/flxblk.F90	(revision 105)
@@ -16,4 +16,5 @@
    USE oce             ! ocean dynamics and tracers
    USE dom_oce         ! ocean space and time domain
+   USE cpl_oce         ! ???
    USE phycst          ! physical constants
    USE daymod
@@ -32,55 +33,55 @@
 
    !! * Module variables
-      INTEGER, PARAMETER  ::   &
-         jpintsr = 24          ! number of time step between sunrise and sunset
-         !                     ! uses for heat flux computation
-      LOGICAL ::   &
-         lbulk_init = .TRUE.   ! flag, bulk initialization done or not)
-
-      REAL(wp), DIMENSION(jpi,jpj) ::   &
-         stauc            ,  &   ! cloud optical depth 
-         sbudyko   
-
-      !! * constants for bulk computation (flx_blk)
-      REAL(wp), DIMENSION(19)  ::  &
-         budyko                  ! BUDYKO's coefficient
-      ! BUDYKO's coefficient (cloudiness effect on LW radiation):
-      DATA budyko / 1.00, 0.98, 0.95, 0.92, 0.89, 0.86, 0.83, 0.80, 0.78, 0.75,  &
-                    0.72, 0.69, 0.67, 0.64, 0.61, 0.58, 0.56, 0.53, 0.50 /
-      REAL(wp), DIMENSION(20)  :: &
-         tauco                  ! cloud optical depth coefficient
-      ! Cloud optical depth coefficient
-      DATA tauco / 6.6, 6.6, 7.0, 7.2, 7.1, 6.8, 6.5, 6.6, 7.1, 7.6,   &
-                   6.6, 6.1, 5.6, 5.5, 5.8, 5.8, 5.6, 5.6, 5.6, 5.6 /
-      REAL(wp)  ::            &  ! constant values
-         zeps    = 1e-20   ,  &
-         zeps0   = 1e-13   ,  &
-         zeps1   = 1e-06   ,  &
-         zzero   = 0.0     ,  &
-         zone    = 1.0
-
-      !! * constants for albedo computation (flx_blk_albedo)
-      REAL(wp) ::   &
-         c1     = 0.05  ,     &   ! constants values
-         c2     = 0.1   ,     &
-         albice = 0.50  ,     &   !  albedo of melting ice in the arctic and antarctic (Shine & Hendersson-Sellers)
-         cgren  = 0.06  ,     &   !  correction of the snow or ice albedo to take into account
-                                  !  effects of cloudiness (Grenfell & Perovich, 1984)
-         alphd  = 0.80  ,     &   !  coefficients for linear interpolation used to compute
-         alphdi = 0.72  ,     &   !  albedo between two extremes values (Pyane, 1972)
-         alphc  = 0.65  ,     &
-         zmue   = 0.4             !  cosine of local solar altitude
-
-      !! * constants for solar declinaison computation (flx_blk_declin)
-      REAL(wp) ::                &
-         a0  =  0.39507671   ,   &  ! coefficients
-         a1  = 22.85684301   ,   &
-         a2  = -0.38637317   ,   &
-         a3  =  0.15096535   ,   &
-         a4  = -0.00961411   ,   &
-         b1  = -4.29692073   ,   &
-         b2  =  0.05702074   ,   &
-         b3  = -0.09028607   ,   &
-         b4  =  0.00592797
+   INTEGER, PARAMETER  ::   &
+      jpintsr = 24          ! number of time step between sunrise and sunset
+      !                     ! uses for heat flux computation
+   LOGICAL ::   &
+      lbulk_init = .TRUE.   ! flag, bulk initialization done or not)
+
+   REAL(wp), DIMENSION(jpi,jpj) ::   &
+      stauc            ,  &   ! cloud optical depth 
+      sbudyko   
+
+   !! * constants for bulk computation (flx_blk)
+   REAL(wp), DIMENSION(19)  ::  &
+      budyko                  ! BUDYKO's coefficient
+   ! BUDYKO's coefficient (cloudiness effect on LW radiation):
+   DATA budyko / 1.00, 0.98, 0.95, 0.92, 0.89, 0.86, 0.83, 0.80, 0.78, 0.75,  &
+      &          0.72, 0.69, 0.67, 0.64, 0.61, 0.58, 0.56, 0.53, 0.50 /
+   REAL(wp), DIMENSION(20)  :: &
+      tauco                  ! cloud optical depth coefficient
+   ! Cloud optical depth coefficient
+   DATA tauco / 6.6, 6.6, 7.0, 7.2, 7.1, 6.8, 6.5, 6.6, 7.1, 7.6,   &
+      &         6.6, 6.1, 5.6, 5.5, 5.8, 5.8, 5.6, 5.6, 5.6, 5.6 /
+   REAL(wp)  ::            &  ! constant values
+      zeps    = 1.e-20  ,  &
+      zeps0   = 1.e-13  ,  &
+      zeps1   = 1.e-06  ,  &
+      zzero   = 0.e0    ,  &
+      zone    = 1.0
+
+   !! * constants for albedo computation (flx_blk_albedo)
+   REAL(wp) ::   &
+      c1     = 0.05  ,     &   ! constants values
+      c2     = 0.10  ,     &
+      albice = 0.50  ,     &   !  albedo of melting ice in the arctic and antarctic (Shine & Hendersson-Sellers)
+      cgren  = 0.06  ,     &   !  correction of the snow or ice albedo to take into account
+                               !  effects of cloudiness (Grenfell & Perovich, 1984)
+      alphd  = 0.80  ,     &   !  coefficients for linear interpolation used to compute
+      alphdi = 0.72  ,     &   !  albedo between two extremes values (Pyane, 1972)
+      alphc  = 0.65  ,     &
+      zmue   = 0.40            !  cosine of local solar altitude
+
+   !! * constants for solar declinaison computation (flx_blk_declin)
+   REAL(wp) ::                &
+      a0  =  0.39507671   ,   &  ! coefficients
+      a1  = 22.85684301   ,   &
+      a2  = -0.38637317   ,   &
+      a3  =  0.15096535   ,   &
+      a4  = -0.00961411   ,   &
+      b1  = -4.29692073   ,   &
+      b2  =  0.05702074   ,   &
+      b3  = -0.09028607   ,   &
+      b4  =  0.00592797
    !!----------------------------------------------------------------------
    !!   OPA 9.0 , LODYC-IPSL  (2003)
@@ -236,21 +237,21 @@
          DO jj = 1, jpj  
             DO ji = 1 , jpi
-            zalat          = ( 90.0 - ABS( gphit(ji,jj) ) ) / 5.0
-            zclat          = ( 95.0 -      gphit(ji,jj)   ) / 10.0
-            indxb          = 1 + INT( zalat ) 
-            !  correction factor to account for the effect of clouds 
-            sbudyko(ji,jj) = budyko(indxb)  
-            indxc          = 1 + INT( zclat )  
-            zdl            = zclat - INT( zclat ) 
-            zdr            = 1.0 - zdl
-            stauc(ji,jj)   = zdr * tauco( indxc ) + zdl * tauco( indxc + 1 ) 
+               zalat          = ( 90.e0 - ABS( gphit(ji,jj) ) ) /  5.e0
+               zclat          = ( 95.e0 -      gphit(ji,jj)   ) / 10.e0
+               indxb          = 1 + INT( zalat ) 
+               !  correction factor to account for the effect of clouds 
+               sbudyko(ji,jj) = budyko(indxb)  
+               indxc          = 1 + INT( zclat )  
+               zdl            = zclat - INT( zclat ) 
+               zdr            = 1.0 - zdl
+               stauc(ji,jj)   = zdr * tauco( indxc ) + zdl * tauco( indxc + 1 ) 
             END DO
          END DO
          IF( nleapy == 1 ) THEN
-            yearday = 366.0
+            yearday = 366.e0
          ELSE IF( nleapy == 0 ) THEN
-            yearday = 365.0
+            yearday = 365.e0
          ELSEIF( nleapy == 30) THEN
-            yearday = 360.0
+            yearday = 360.e0
          ENDIF
          lbulk_init = .FALSE.
@@ -264,21 +265,21 @@
       zqsb_ice(:,:) = 0.e0
 
-      zpis2       = rpi / 2.              ! pi / 2
-      z2pi        = 2. * rpi              ! 2 * pi 
+      zpis2       = rpi / 2.
+      z2pi        = 2. * rpi
 
  !CDIR NOVERRCHK
-     DO jj = 1, jpj
+      DO jj = 1, jpj
  !CDIR NOVERRCHK
-        DO ji = 1, jpi
-
-      ztatm (ji,jj) = 273.15 + tatm  (ji,jj)  !  air temperature in Kelvins 
-      zcatm1(ji,jj) = 1.0    - catm  (ji,jj)  !  fractional cloud cover
-      zfrld (ji,jj) = 1.0    - freeze(ji,jj)  !  fractional sea ice cover
-      zpatm(ji,jj)  = 101000.               !  pressure 
-      
-      !  Computation of air density, obtained from the equation of state for dry air. 
-      zrhoa(ji,jj) = zpatm(ji,jj) / ( 287.04 * ztatm(ji,jj) )
-      
-      !  zes : Saturation water vapour
+         DO ji = 1, jpi
+
+            ztatm (ji,jj) = 273.15 + tatm  (ji,jj)  !  air temperature in Kelvins 
+            zcatm1(ji,jj) = 1.0    - catm  (ji,jj)  !  fractional cloud cover
+            zfrld (ji,jj) = 1.0    - freeze(ji,jj)  !  fractional sea ice cover
+            zpatm(ji,jj)  = 101000.               !  pressure 
+      
+            !  Computation of air density, obtained from the equation of state for dry air. 
+            zrhoa(ji,jj) = zpatm(ji,jj) / ( 287.04 * ztatm(ji,jj) )
+      
+            !  zes : Saturation water vapour
             ztamr = ztatm(ji,jj) - rtt
             zmt1  = SIGN( 17.269, ztamr )
@@ -288,16 +289,16 @@
                &                      / ( ztatm(ji,jj) - 35.86  + MAX( zzero, zmt3 ) ) )
 
-      !  zev : vapour pressure  (hatm is relative humidity)  
-      zev(ji,jj)   = hatm(ji,jj) * zes(ji,jj) 
-      !  square-root of vapour pressure
+            !  zev : vapour pressure  (hatm is relative humidity)  
+            zev(ji,jj)   = hatm(ji,jj) * zes(ji,jj) 
+            !  square-root of vapour pressure
 !CDIR NOVERRCHK
-      zevsqr(ji,jj) = SQRT( zev(ji,jj) * 0.01 )
-      !  zqapb  : specific humidity 
-      zqatm(ji,jj) = 0.622 * zev(ji,jj) / ( zpatm(ji,jj) - 0.378 * zev(ji,jj) )
-
-
-      !----------------------------------------------------
-      !   Computation of snow precipitation (Ledley, 1985) |
-      !----------------------------------------------------
+            zevsqr(ji,jj) = SQRT( zev(ji,jj) * 0.01 )
+            !  zqapb  : specific humidity 
+            zqatm(ji,jj) = 0.622 * zev(ji,jj) / ( zpatm(ji,jj) - 0.378 * zev(ji,jj) )
+
+
+            !----------------------------------------------------
+            !   Computation of snow precipitation (Ledley, 1985) |
+            !----------------------------------------------------
 
             zmt1  =   253.0 - ztatm(ji,jj)
@@ -332,5 +333,5 @@
 
       iday   = INT( zxday )
-      IF( l_ctl .AND. lwp ) WRITE(numout,*) ' declin : iday ', iday, ' nfbulk= ', nfbulk
+      IF(l_ctl)   WRITE(numout,*) ' declin : iday ', iday, ' nfbulk= ', nfbulk
       !   computation of the solar declination, his sine and his cosine
       CALL flx_blk_declin( indaet, iday, zdecl )
@@ -350,18 +351,18 @@
          DO ji = 1, jpi
             !  product of sine of latitude and sine of solar declination
-            zps(ji,jj)     = SIN( gphit(ji,jj) * rad ) * zsdecl
+            zps     (ji,jj) = SIN( gphit(ji,jj) * rad ) * zsdecl
             !  product of cosine of latitude and cosine of solar declination
-            zpc(ji,jj)     = COS( gphit(ji,jj) * rad ) * zcdecl
+            zpc     (ji,jj) = COS( gphit(ji,jj) * rad ) * zcdecl
             !  computation of the both local time of sunrise and sunset
-            zlsrise(ji,jj) =  ACOS( - SIGN( zone, zps(ji,jj) ) * MIN( zone, SIGN( zone, zps(ji,jj) )  &
+            zlsrise (ji,jj) = ACOS( - SIGN( zone, zps(ji,jj) ) * MIN( zone, SIGN( zone, zps(ji,jj) )  &
                &                     * ( zps(ji,jj) / zpc(ji,jj) ) ) ) 
-            zlsset(ji,jj)  = - zlsrise(ji,jj)
+            zlsset  (ji,jj) = - zlsrise(ji,jj)
             !  dividing the solar day into jpintsr segments of length zdlha
-            zdlha(ji,jj)   =  ( zlsrise(ji,jj) - zlsset(ji,jj) ) / REAL( jpintsr )
+            zdlha   (ji,jj) = ( zlsrise(ji,jj) - zlsset(ji,jj) ) / REAL( jpintsr )
             !  computation of the local noon solar altitude
-            zlmunoon(ji,jj)=  ASIN ( ( zps(ji,jj) + zpc(ji,jj) ) ) / rad
+            zlmunoon(ji,jj) = ASIN ( ( zps(ji,jj) + zpc(ji,jj) ) ) / rad
             
             !  cloud correction taken from Reed (1977) (imposed lower than 1)
-            zcldcor(ji,jj) = MIN( zone, ( 1 - 0.62 * catm(ji,jj) + 0.0019 * zlmunoon(ji,jj) ) )
+            zcldcor (ji,jj) = MIN( zone, ( 1.e0 - 0.62 * catm(ji,jj) + 0.0019 * zlmunoon(ji,jj) ) )
          END DO
       END DO
@@ -379,5 +380,5 @@
             DO ji = 1, jpi
                !  local hour angle
-               zlha (ji,jj,jt)= COS ( zlsrise(ji,jj) - zcoef * zdlha(ji,jj) )
+               zlha (ji,jj,jt) = COS ( zlsrise(ji,jj) - zcoef * zdlha(ji,jj) )
 
                ! cosine of local solar altitude
@@ -416,26 +417,26 @@
       !--------------------------------------------------------------------
 
-      zalbocsd(:,:) = 0.  
-      zqsro   (:,:) = 0.
-      zqsrics (:,:) = 0.
-      zqsrios (:,:) = 0.
+      zalbocsd(:,:) = 0.e0
+      zqsro   (:,:) = 0.e0
+      zqsrics (:,:) = 0.e0
+      zqsrios (:,:) = 0.e0
 
       DO jt = 1, jpintsr 
 #   if defined key_vectopt_loop && ! defined key_autotasking
          DO ji = 1, jpij  
-            zalbocsd(ji,1)  = zalbocsd(ji,1) + zdlha(ji,1) * zalbocs(ji,1,jt)   &
-               &             / MAX( 2.0 * zlsrise(ji,1) , zeps0 )
-            zqsro   (ji,1)  =   zqsro(ji,1) + zsqsro  (ji,1,jt)
-            zqsrics (ji,1)  = zqsrics(ji,1) + zsqsrics(ji,1,jt)
-            zqsrios (ji,1)  = zqsrios(ji,1) + zsqsrios(ji,1,jt)
+            zalbocsd(ji,1) = zalbocsd(ji,1) + zdlha   (ji,1) * zalbocs(ji,1,jt)   &
+               &                                             / MAX( 2.0 * zlsrise(ji,1) , zeps0 )
+            zqsro   (ji,1) = zqsro   (ji,1) + zsqsro  (ji,1,jt)
+            zqsrics (ji,1) = zqsrics (ji,1) + zsqsrics(ji,1,jt)
+            zqsrios (ji,1) = zqsrios (ji,1) + zsqsrios(ji,1,jt)
          END DO
 #  else
          DO jj = 1, jpj
             DO ji = 1, jpi  
-               zalbocsd(ji,jj)  = zalbocsd(ji,jj) + zdlha(ji,jj) * zalbocs(ji,jj,jt)   &
-                  &             / MAX( 2.0 * zlsrise(ji,jj) , zeps0 )
-               zqsro  (ji,jj)   =   zqsro(ji,jj) + zsqsro  (ji,jj,jt)
-               zqsrics(ji,jj)   = zqsrics(ji,jj) + zsqsrics(ji,jj,jt)
-               zqsrios(ji,jj)   = zqsrios(ji,jj) + zsqsrios(ji,jj,jt)
+               zalbocsd(ji,jj) = zalbocsd(ji,jj) + zdlha(ji,jj) * zalbocs(ji,jj,jt)   &
+                  &                                              / MAX( 2.0 * zlsrise(ji,jj) , zeps0 )
+               zqsro  (ji,jj)  = zqsro   (ji,jj) + zsqsro  (ji,jj,jt)
+               zqsrics(ji,jj)  = zqsrics (ji,jj) + zsqsrics(ji,jj,jt)
+               zqsrios(ji,jj)  = zqsrios (ji,jj) + zsqsrios(ji,jj,jt)
             END DO
          END DO
@@ -446,7 +447,7 @@
          DO ji = 1, jpi 
 
-      !------------------------------------------- 
-      !  Computation of shortwave radiation.
-      !-------------------------------------------
+            !------------------------------------------- 
+            !  Computation of shortwave radiation.
+            !-------------------------------------------
 
             ! the solar heat flux absorbed at ocean and snow/ice surfaces
@@ -477,7 +478,7 @@
             fr2_i0(ji,jj) = 0.82  * zcatm1(ji,jj) + 0.65 * catm(ji,jj)
 
-      !---------------------------------------------------------------------------
-      !   Computation of long-wave radiation  ( Berliand 1952 ; all latitudes )
-      !---------------------------------------------------------------------------
+            !---------------------------------------------------------------------------
+            !   Computation of long-wave radiation  ( Berliand 1952 ; all latitudes )
+            !---------------------------------------------------------------------------
 
             ! tempory variables
@@ -555,5 +556,5 @@
             zpsil          = zpsih
             
-            zvatmg         = MAX( 0.032 * 1.5e-3 * vatm(ji,jj) * vatm(ji,jj) / g, zeps )
+            zvatmg         = MAX( 0.032 * 1.5e-3 * vatm(ji,jj) * vatm(ji,jj) / grav, zeps )
 !i
 !!          if( zvatmg == 0.e0 ) then
@@ -632,5 +633,5 @@
             !  latent heat flux 
             zqla_ice(ji,jj) = zrhovacshi * ( zqsati        - zqatm(ji,jj) )
-             qla_ice(ji,jj) = zqla_ice(ji,jj)
+            qla_ice (ji,jj) = zqla_ice(ji,jj)
               
             !  Computation of sensitivity of non solar fluxes (dQ/dT)
@@ -749,5 +750,5 @@
       !-------------------------- 
       
-      llmask = (hsnif == 0.0) .AND. ( sist >= rt0_ice )
+      llmask = (hsnif == 0.e0) .AND. ( sist >= rt0_ice )
       WHERE ( llmask )   !  ice free of snow and melts
          zalbfz = albice
Index: tags/nemo_dev_x3/NEMO/OPA_SRC/SBC/flxfwb.F90
===================================================================
--- branches/nemo/NEMO/OPA_SRC/SBC/flxfwb.F90	(revision 4)
+++ tags/nemo_dev_x3/NEMO/OPA_SRC/SBC/flxfwb.F90	(revision 105)
@@ -6,4 +6,6 @@
 #if defined key_dynspg_fsc
    !!----------------------------------------------------------------------
+   !!   'key_dynspg_fsc'               No free surface with constant volume
+   !!----------------------------------------------------------------------
    !!   flx_fwb      : freshwater budget for global ocean configurations
    !!                  in free surface and forced mode
@@ -13,4 +15,5 @@
    USE oce             ! ocean dynamics and tracers
    USE dom_oce         ! ocean space and time domain
+   USE cpl_oce         ! coupled atmosphere/ocean
    USE phycst          ! physical constants
    USE in_out_manager  ! I/O manager
@@ -30,9 +33,9 @@
 
    !! * Share module variables
-   LOGICAL, PUBLIC ::    &  !!! * namelist *
-      ln_fwb = .TRUE.        ! Flag to activate the fwb computation
-   REAL(wp), PUBLIC ::   &
-      a_fwb_b  ,  &  ! annual domain averaged freshwater budget
-      a_fwb          ! for 2 year before (_b) and before year.
+   LOGICAL, PUBLIC ::    & !!: * namelist *
+      ln_fwb = .TRUE.       !: Flag to activate the fwb computation
+   REAL(wp), PUBLIC ::   &  !:
+      a_fwb_b  ,         &  !: annual domain averaged freshwater budget
+      a_fwb                 !: for 2 year before (_b) and before year.
 
    !! * Module variables
@@ -103,4 +106,5 @@
          e1e2_i(:,:) = e1t(:,:) * e2t(:,:) * tmask_i(:,:)
          a_sshb = SUM( e1e2_i(:,:) * sshn(:,:) )
+         IF( lk_mpp )   CALL  mpp_sum( a_sshb   )   ! sum over all the global domain
          DO jk = 1, jpkm1
             DO jj = 2, jpjm1
@@ -111,9 +115,5 @@
             END DO
          END DO
-#if defined key_mpp
-      ! Mpp: sum over all the global domain
-      CALL  mpp_sum( a_sshb   )        
-      CALL  mpp_sum( a_sal000 )        
-#endif
+         IF( lk_mpp )   CALL  mpp_sum( a_sal000 )   ! sum over the global domain
 
       ENDIF
@@ -122,38 +122,29 @@
       ! --------------------------------------====----------
       a_emp    = SUM( e1e2_i(:,:) * emp   (:,:) )
+      IF( lk_mpp )   CALL  mpp_sum( a_emp    )   ! sum over the global domain
 #if defined key_flx_bulk_monthly || defined key_flx_bulk_daily
       a_precip = SUM( e1e2_i(:,:) * watm  (:,:) )
+      IF( lk_mpp )   CALL  mpp_sum( a_precip )   ! sum over the global domain
 #endif
       a_rnf    = SUM( e1e2_i(:,:) * runoff(:,:) )
-
-      IF( aminus /= 0.0 ) a_aminus = a_aminus + ( MIN( aplus, aminus ) / aminus )
-      IF( aplus  /= 0.0 ) a_aplus  = a_aplus  + ( MIN( aplus, aminus ) / aplus  )
-
-#if defined key_mpp
-      ! Mpp: sum over all the global domain
-      CALL  mpp_sum( a_emp    )   
-#if defined key_flx_bulk_monthly || defined key_flx_bulk_daily
-      CALL  mpp_sum( a_precip )   
-#endif
-      CALL  mpp_sum( a_rnf    )   
-#endif
+      IF( lk_mpp )   CALL  mpp_sum( a_rnf    )   ! sum over the global domain
+
+      IF( aminus /= 0.e0 ) a_aminus = a_aminus + ( MIN( aplus, aminus ) / aminus )
+      IF( aplus  /= 0.e0 ) a_aplus  = a_aplus  + ( MIN( aplus, aminus ) / aplus  )
 
 
       ! Update empold if new year start
-      ikty = 365 * 86400 /rdttra(1)
+      ikty = 365 * 86400 / rdttra(1)    !!bug  use of 365 days leap year or 360d year !!!!!!!
       IF( MOD( kt, ikty ) == 0 ) THEN
-            zarea    = SUM( e1e2_i(:,:)             )
-            a_fwb_b = a_fwb
-            a_fwb   = SUM( e1e2_i(:,:) * sshn(:,:) )
-#if defined key_mpp
-      ! Mpp: sum over all the global domain
-      CALL  mpp_sum( zarea    )        
-      CALL  mpp_sum( a_fwb    )        
-#endif
-            a_fwb   = a_fwb * 1.e+3 / ( zarea * 86400. * 365. )    ! convert in Kg/m3/s = mm/s
-
-            empold = 2. * a_fwb - a_fwb_b   ! current year freshwater budget correction
-            !                               ! estimate from the 2 previous years budget
-
+         zarea    = SUM( e1e2_i(:,:)            )
+         IF( lk_mpp )   CALL  mpp_sum( zarea    )   ! sum over the global domain
+         a_fwb_b = a_fwb
+         a_fwb   = SUM( e1e2_i(:,:) * sshn(:,:) )
+         IF( lk_mpp )   CALL  mpp_sum( a_fwb    )   ! sum over the global domain
+
+         a_fwb   = a_fwb * 1.e+3 / ( zarea * 86400. * 365. )    ! convert in Kg/m3/s = mm/s
+         !                                                      !!bug 365d year 
+         empold = 2. * a_fwb - a_fwb_b   ! current year freshwater budget correction
+         !                               ! estimate from the 2 previous years budget
       ENDIF
 
@@ -176,10 +167,9 @@
             END DO
          END DO
-#if defined key_mpp
-      ! Mpp: sum over all the global domain
-      CALL  mpp_sum( zarea    )        
-      CALL  mpp_sum( a_sshend )        
-      CALL  mpp_sum( a_salend )        
-#endif
+         IF( lk_mpp ) THEN
+            CALL  mpp_sum( zarea    )   ! sums over all the global domain
+            CALL  mpp_sum( a_sshend )        
+            CALL  mpp_sum( a_salend )        
+         ENDIF
 
          a_aminus = a_aminus / ( nitend - nit000 + 1 )
@@ -281,12 +271,12 @@
       IF(lwp) WRITE '               Rigid-lid option, fwb correction is useless, but valid'
 #endif
-#if defined key_coupled
-      IF(lwp) WRITE '               Coupled option, fwb correction is a flux correction ! '
-      IF(lwp) WRITE '               ln_fwb = .FALSE. is recommanded'
-#endif
-
-      !                                       ! ==============================
+      IF( lk_cpl ) THEN
+         IF(lwp) WRITE(numout,*) '               Coupled option, fwb correction is a flux correction ! '
+         IF(lwp) WRITE(numout,*) '               ln_fwb = .FALSE. is recommanded'
+      ENDIF
+
+      !                                        ! ==============================
       IF( ln_fwb ) THEN                        !  Freshwater budget correction 
-         !                                    ! ==============================
+         !                                     ! ==============================
          ! Read the corrective factor on precipitations (empold)
          INQUIRE( FILE=clname, EXIST=llbon )
@@ -317,6 +307,11 @@
 
 #else
+   !!----------------------------------------------------------------------
+   !!   Default case :                       
+   !!----------------------------------------------------------------------
+   LOGICAL, PUBLIC ::   ln_fwb = .FALSE.   !: no fwb forced
 CONTAINS
-   SUBROUTINE flx_fwb
+   SUBROUTINE flx_fwb( kt )                ! dummy routine
+      WRITE(*,*) 'flx_fwb: You should not have seen this print! error?', kt
    END SUBROUTINE flx_fwb
    SUBROUTINE flx_fwb_init
Index: tags/nemo_dev_x3/NEMO/OPA_SRC/SBC/flxrnf.F90
===================================================================
--- branches/nemo/NEMO/OPA_SRC/SBC/flxrnf.F90	(revision 4)
+++ tags/nemo_dev_x3/NEMO/OPA_SRC/SBC/flxrnf.F90	(revision 105)
@@ -27,20 +27,21 @@
    PRIVATE
 
+   !! * Routine accessibility
+   PUBLIC flx_rnf          ! routine call in step module
+
    !! * Shared module variables
-   PUBLIC flx_rnf         ! routine call by step.F90
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   &
-      runoff,          &  ! monthly runoff (kg/m2/s)
-      upsadv,          &  ! mixed adv scheme in straits vicinity (hori.)
-      upsrnfh             ! mixed adv scheme in runoffs vicinity (hori.)
-   REAL(wp), PUBLIC, DIMENSION(jpk) ::   &
-      upsrnfz             ! mixed adv scheme in runoffs vicinity (vert.)
-   INTEGER, PUBLIC ::   &
-      numrof  = 48 ,    &  ! logical unit for runoff data
-      nrunoff =  0 ,    &  ! runoff option (namelist)
-      nrnf1, nrnf2         ! first and second record used
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   &  !:
+      runoff,           &  !: monthly runoff (kg/m2/s)
+      upsadv,           &  !: mixed adv scheme in straits vicinity (hori.)
+      upsrnfh              !: mixed adv scheme in runoffs vicinity (hori.)
+   REAL(wp), PUBLIC, DIMENSION(jpk) ::   &  !:
+      upsrnfz              !: mixed adv scheme in runoffs vicinity (vert.)
+   INTEGER, PUBLIC ::   &  !:
+      nrunoff =  0 ,    &  !: runoff option (namelist)
+      nrnf1, nrnf2         !: first and second record used
 
    !! * Module variable
-   REAL(wp), DIMENSION(jpi,jpj,2) ::   &
-      rnfdta              ! monthly runoff data array (kg/m2/s)
+   REAL(wp), DIMENSION(jpi,jpj,2) ::   &  !:
+      rnfdta               !: monthly runoff data array (kg/m2/s)
    !!----------------------------------------------------------------------
    !!   OPA 9.0 , LODYC-IPSL  (2003)
@@ -104,7 +105,8 @@
       INTEGER  ::   ji, jj                 ! dummy loop indices
       INTEGER  ::   ipi, ipj, ipk          ! temporary integers
+      INTEGER  ::   ii0, ii1, ij0, ij1     !    "          "
       INTEGER, DIMENSION(jpmois) ::     &
          istep                             ! temporary workspace
-      REAL(wp) ::   zdate0                 ! temporary scalars
+      REAL(wp) ::   zdate0, zdt            ! temporary scalars
       REAL(wp), DIMENSION(jpk) ::       &
          zlev                              ! temporary workspace
@@ -135,5 +137,5 @@
          CASE DEFAULT
             IF(lwp) WRITE(numout,cform_err)
-            IF(lwp) WRITE(numout,*)' Error nrunoff = ', nrunoff, ' /= 0, 1 or 2'
+            IF(lwp) WRITE(numout,*) ' Error nrunoff = ', nrunoff, ' /= 0, 1 or 2'
             nstop = nstop + 1
 
@@ -175,4 +177,5 @@
          ipj = jpjglo
          ipk = jpk
+         zdt = rdt
          
          ! Open file
@@ -180,7 +183,7 @@
          IF( kt == nit000 ) THEN
             CALL flinopen( clname, mig(1), nlci, mjg(1), nlcj,    &
-                           .false., ipi, ipj, ipk, zlon,        &
-                           zlat, zlev, jpmois, istep, zdate0,   &
-                           rdt, numrof )
+               &           .false., ipi, ipj, ipk, zlon,        &
+               &           zlat, zlev, jpmois, istep, zdate0,   &
+               &           zdt, numrnf )
             !   Title, dimensions and tests
 # if ! defined key_coupled
@@ -191,8 +194,8 @@
                nstop = nstop + 1
             ENDIF
-            IF(lwp) WRITE(numout,*) iman, istep, zdate0, rdt, numrof
-            IF(lwp) WRITE(numout,*) 'numrof=', numrof
+            IF(lwp) WRITE(numout,*) iman, istep, zdate0, rdt, numrnf
+            IF(lwp) WRITE(numout,*) 'numrnf=', numrnf
             IF(lwp) WRITE(numout,*) 'jpmois=', jpmois
-            IF(lwp) WRITE(numout,*) 'rdt=', rdt
+            IF(lwp) WRITE(numout,*) 'zdt=', zdt
 # endif
             IF(ipi /= jpidta .AND. ipj /= jpjdta .AND. ipk /= 1) THEN
@@ -237,9 +240,11 @@
             
             ! Read monthly runoff data in kg/m2/s
-
-            CALL flinget(numrof,'sorunoff',jpidta,jpjdta,1,jpmois   &
-                  ,nrnf1,nrnf1,mig(1),nlci,mjg(1),nlcj,rnfdta(1:nlci,1:nlcj,1))
-            CALL flinget(numrof,'sorunoff',jpidta,jpjdta,1,jpmois   &
-                  ,nrnf2,nrnf2,mig(1),nlci,mjg(1),nlcj,rnfdta(1:nlci,1:nlcj,2))
+!ibug
+            IF( kt == nit000 )   rnfdta(:,:,:) = 0.e0
+!ibug
+            CALL flinget( numrnf, 'sorunoff', jpidta, jpjdta, 1, jpmois   &
+               &        , nrnf1, nrnf1, mig(1), nlci, mjg(1), nlcj, rnfdta(1:nlci,1:nlcj,1) )
+            CALL flinget( numrnf, 'sorunoff', jpidta, jpjdta, 1, jpmois   &
+               &        , nrnf2, nrnf2, mig(1), nlci, mjg(1), nlcj, rnfdta(1:nlci,1:nlcj,2) )
 
             IF(lwp) WRITE(numout,*)
@@ -254,12 +259,10 @@
          zxy = FLOAT( nday + idmeom - idbd * i15 ) / idbd
 
-         runoff(:,:) = -( (1.-zxy)*rnfdta(:,:,1) + zxy*rnfdta(:,:,2) )
+         runoff(:,:) = -( ( 1.e0 - zxy ) * rnfdta(:,:,1) + zxy * rnfdta(:,:,2) )
 
          ! Runoff reduction
          DO jj = 1, jpj
             DO ji = 1, jpi
-               IF( gphit(ji,jj) > 40 .AND. gphit(ji,jj) < 65 ) THEN
-                  runoff(ji,jj) = 0.85 * runoff(ji,jj)
-               ENDIF
+               IF( gphit(ji,jj) > 40 .AND. gphit(ji,jj) < 65 )   runoff(ji,jj) = 0.85 * runoff(ji,jj)
             END DO
          END DO
@@ -281,7 +284,10 @@
          !  coefr * upstream + (1- coefr) centered
          !  coefr must be between 0 and 1.
-
-         CALL flinget(numrof,'socoefr',jpidta,jpjdta,1,jpmois,nrnf1   &
-                     ,nrnf1,mig(1),nlci,mjg(1),nlcj,zcoefr(1:nlci,1:nlcj))
+!ibug
+         zcoefr(:,:) = 0.e0
+!ibug
+
+         CALL flinget( numrnf, 'socoefr', jpidta, jpjdta, 1, jpmois, nrnf1,   &
+            &          nrnf1, mig(1), nlci, mjg(1), nlcj, zcoefr(1:nlci,1:nlcj) )
 
          IF(lwp) WRITE(numout,*)
@@ -298,9 +304,9 @@
          
          IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN
-            ! ORCA_R2 configuration :
-            ! upstream scheme in the Sound Strait
-            upsrnfh( mi0(144):mi1(144) , mj0(116):mj1(116) ) = 0.25
-            upsrnfh( mi0(145):mi1(147) , mj0(116):mj1(116) ) = 0.50
-            upsrnfh( mi0(148):mi1(148) , mj0(116):mj1(116) ) = 0.25
+            ! ORCA_R2 configuration : upstream scheme in the Sound Strait
+            ij0 = 116   ;   ij1 = 116
+            ii0 = 144   ;   ii1 = 144   ;   upsrnfh( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.25
+            ii0 = 145   ;   ii1 = 147   ;   upsrnfh( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.50
+            ii0 = 148   ;   ii1 = 148   ;   upsrnfh( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.25
          ENDIF
 
@@ -318,6 +324,7 @@
                !                                     ! =======================
 
-               !                                                       ! Gibraltar Strait
-               upsadv( mi0(70):mi1(71) , mj0(52):mj1(53) ) = 0.50
+               !                                          ! Gibraltar Strait
+               ii0 =  70   ;   ii1 =  71
+               ij0 =  52   ;   ij1 =  53   ;   upsadv( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.50
          
                !                                     ! =======================
@@ -326,19 +333,26 @@
 
                !                                          ! Gibraltar Strait
-               upsadv( mi0(138):mi1(138) , mj0(102):mj1(102) ) = 0.20
-               upsadv( mi0(139):mi1(139) , mj0(102):mj1(102) ) = 0.40
-               upsadv( mi0(140):mi1(140) , mj0(102):mj1(102) ) = 0.50
-               upsadv( mi0(141):mi1(141) , mj0(101):mj1(102) ) = 0.50
+               ij0 = 102   ;   ij1 = 102
+               ii0 = 138   ;   ii1 = 138   ;   upsadv( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.20
+               ii0 = 139   ;   ii1 = 139   ;   upsadv( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.40
+               ii0 = 140   ;   ii1 = 140   ;   upsadv( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.50
+               ij0 = 101   ;   ij1 = 102
+               ii0 = 141   ;   ii1 = 141   ;   upsadv( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.50
 
                !                                          ! Bab el Mandeb Strait
-               upsadv( mi0(164):mi1(164) , mj0( 87):mj1( 88) ) = 0.10
-               upsadv( mi0(163):mi1(163) , mj0( 88):mj1( 88) ) = 0.25
-               upsadv( mi0(162):mi1(162) , mj0( 88):mj1( 88) ) = 0.40
-               upsadv( mi0(160):mi1(161) , mj0( 88):mj1( 88) ) = 0.50
-               upsadv( mi0(158):mi1(160) , mj0( 89):mj1( 89) ) = 0.25
-               upsadv( mi0(160):mi1(160) , mj0( 90):mj1( 90) ) = 0.25
+               ij0 =  87   ;   ij1 =  88
+               ii0 = 164   ;   ii1 = 164   ;   upsadv( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.10
+               ij0 =  88   ;   ij1 =  88
+               ii0 = 163   ;   ii1 = 163   ;   upsadv( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.25
+               ii0 = 162   ;   ii1 = 162   ;   upsadv( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.40
+               ii0 = 160   ;   ii1 = 161   ;   upsadv( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.50
+               ij0 =  89   ;   ij1 =  89
+               ii0 = 158   ;   ii1 = 160   ;   upsadv( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.25
+               ij0 =  90   ;   ij1 =  90
+               ii0 = 160   ;   ii1 = 160   ;   upsadv( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.25
 
                !                                          ! Sound Strait
-               upsadv( mi0(145):mi1(147) , mj0(116):mj1(116) ) = 0.50
+               ij0 = 116   ;   ij1 = 116
+               ii0 = 145   ;   ii1 = 147   ;   upsadv( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.50
          
             END SELECT 
@@ -351,5 +365,5 @@
       ! --------------------
 
-      IF( kt == nitend .AND. nrunoff >= 1 )   CALL flinclo( numrof )
+      IF( kt == nitend .AND. nrunoff >= 1 )   CALL flinclo( numrnf )
 
    END SUBROUTINE flx_rnf
Index: tags/nemo_dev_x3/NEMO/OPA_SRC/SBC/flxrnf_ORCA_R05.h90
===================================================================
--- branches/nemo/NEMO/OPA_SRC/SBC/flxrnf_ORCA_R05.h90	(revision 4)
+++ tags/nemo_dev_x3/NEMO/OPA_SRC/SBC/flxrnf_ORCA_R05.h90	(revision 105)
@@ -13,20 +13,22 @@
    PRIVATE
 
+   !! * Routine accessibility
+   PUBLIC flx_rnf         ! routine call by step.F90
+
    !! * Shared module variables
-   PUBLIC flx_rnf         ! routine call by step.F90
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   &
-      runoff,          &  ! monthly runoff (kg/m2/s)
-      upsadv,          &  ! mixed adv scheme in straits vicinity (hori.)
-      upsrnfh             ! mixed adv scheme in runoffs vicinity (hori.)
-   REAL(wp), PUBLIC, DIMENSION(jpk) ::   &
-      upsrnfz             ! mixed adv scheme in runoffs vicinity (vert.)
-   INTEGER, PUBLIC ::   &
-      numrof  = 48 ,    &  ! logical unit for runoff data
-      nrunoff =  0 ,    &  ! runoff option (namelist)
-      nrnf1, nrnf2         ! first and second record used
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   &  !:
+      runoff,            &  !: monthly runoff (kg/m2/s)
+      upsadv,            &  !: mixed adv scheme in straits vicinity (hori.)
+      upsrnfh               !: mixed adv scheme in runoffs vicinity (hori.)
+   REAL(wp), PUBLIC, DIMENSION(jpk) ::   &  !:
+      upsrnfz               !: mixed adv scheme in runoffs vicinity (vert.)
+   INTEGER , PUBLIC ::   &  !:
+      numrof  = 48 ,     &  !: logical unit for runoff data
+      nrunoff =  0 ,     &  !: runoff option (namelist)
+      nrnf1, nrnf2          !: first and second record used
 
    !! * Module variable
-   REAL(wp), DIMENSION(jpi,jpj,2) ::   &
-      rnfdta              ! monthly runoff data array (kg/m2/s)
+   REAL(wp), DIMENSION(jpi,jpj,2) ::   &  !:
+      rnfdta                !: monthly runoff data array (kg/m2/s)
    !!----------------------------------------------------------------------
    !!   OPA 9.0 , LODYC-IPSL  (2003)
@@ -88,14 +90,12 @@
       !!----------------------------------------------------------------------
       !! * Local declarations
-      INTEGER  :: jpmois
-      INTEGER  ::   &
-         jpriv,    & ! jpriv : maximum river number
-         jpcoef      ! jpcoef: maximum point number
-      INTEGER, PARAMETER ::   jpmois=12, jpriv=120, jpcoef=20
+      INTEGER, PARAMETER ::   &
+         jpmois=12,   &  ! number of month in the year !
+         jpriv=120,   &  ! maximum number of rivers
+         jpcoef=20       ! maximum number of gridpoints for mouth rivers
 
       INTEGER  ::   ji, jj, jl, jm, jn, jr, jc
       INTEGER  ::   i15, imois, iman
       INTEGER  ::   inbp, idbd, idmeom
-      INTEGER  ::   iipt0, iipt1, ijpt0, ijpt1, ii, ij
       INTEGER  ::   iirnf(jpcoef,jpriv),ijrnf(jpcoef,jpriv)
       INTEGER  ::   inb(jpriv)
@@ -1040,8 +1040,6 @@
                      END DO
                      ! ji-1,jj-2
-                     ijpt0 = MAX( 1, MIN( ij-2,jpj  ) )
-                     ijpt1 = MAX( 0, MIN( ij-2,jpj-1) )
-                     DO jj = ijpt0, ijpt1
-                        DO ji = iipt0, iipt1
+                     DO jj = mj0(ijrnf(jl,jr) - 2), mj1(ijrnf(jl,jr) - 2)
+                        DO ji = mi0(iirnf(jl,jr) - 1), mi1(iirnf(jl,jr) - 1)
                            DO jn = 1, 3 
                               zcoefr(ji,jj,jn) = MAX( zrup3, zcoefr(ji,jj,jn) )
@@ -1050,8 +1048,6 @@
                      END DO
                      ! ji,jj-2
-                     iipt0 = MAX( 1, MIN( ii,jpi  ) )
-                     iipt1 = MAX( 0, MIN( ii,jpi-1) )
-                     DO jj = ijpt0, ijpt1
-                        DO ji = iipt0, iipt1
+                     DO jj = mj0(ijrnf(jl,jr) - 2), mj1(ijrnf(jl,jr) - 2)
+                        DO ji = mi0(iirnf(jl,jr)    ), mi1(iirnf(jl,jr)    )
                            DO jn = 1, 3 
                               zcoefr(ji,jj,jn) = MAX( zrup3, zcoefr(ji,jj,jn) )
@@ -1060,8 +1056,6 @@
                      END DO
                      ! ji,jj+2
-                     ijpt0 = MAX( 1, MIN( ij+2,jpj  ) )
-                     ijpt1 = MAX( 0, MIN( ij+2,jpj-1) )
-                     DO jj = ijpt0, ijpt1
-                        DO ji = iipt0, iipt1
+                     DO jj = mj0(ijrnf(jl,jr) + 2), mj1(ijrnf(jl,jr) + 2)
+                        DO ji = mi0(iirnf(jl,jr)    ), mi1(iirnf(jl,jr)    )
                            DO jn = 1, 3
                               zcoefr(ji,jj,jn) = MAX( zrup3, zcoefr(ji,jj,jn) )
@@ -1070,8 +1064,6 @@
                      END DO
                      ! ji+1,jj-2
-                     iipt0 = MAX( 1, MIN( ii+1,jpi  ) )
-                     iipt1 = MAX( 0, MIN( ii+1,jpi-1) )
-                     DO jj = ijpt0, ijpt1
-                        DO ji = iipt0, iipt1
+                     DO jj = mj0(ijrnf(jl,jr) - 2), mj1(ijrnf(jl,jr) - 2)
+                        DO ji = mi0(iirnf(jl,jr) + 1), mi1(iirnf(jl,jr) + 1)
                            DO jn = 1, 3 
                               zcoefr(ji,jj,jn) = MAX( zrup3, zcoefr(ji,jj,jn) )
@@ -1080,8 +1072,6 @@
                      END DO
                      ! ji+1,jj+2
-                     ijpt0 = MAX( 1, MIN( ij+2,jpj  ) )
-                     ijpt1 = MAX( 0, MIN( ij+2,jpj-1) )
-                     DO jj = ijpt0, ijpt1
-                        DO ji = iipt0, iipt1
+                     DO jj = mj0(ijrnf(jl,jr) + 2), mj1(ijrnf(jl,jr) + 2)
+                        DO ji = mi0(iirnf(jl,jr) + 1), mi1(iirnf(jl,jr) + 1)
                            DO jn = 1, 3 
                               zcoefr(ji,jj,jn) = MAX( zrup3, zcoefr(ji,jj,jn) )
@@ -1090,8 +1080,6 @@
                      END DO
                      ! ji+2,jj-2
-                     iipt0 = MAX( 1, MIN( ii+2,jpi  ) )
-                     iipt1 = MAX( 0, MIN( ii+2,jpi-1) )
-                     DO jj = ijpt0, ijpt1
-                        DO ji = iipt0, iipt1
+                     DO jj = mj0(ijrnf(jl,jr) - 2), mj1(ijrnf(jl,jr) - 2)
+                        DO ji = mi0(iirnf(jl,jr) + 2), mi1(iirnf(jl,jr) + 2)
                            DO jn = 1, 3 
                               zcoefr(ji,jj,jn) = MAX( zrup3, zcoefr(ji,jj,jn) )
@@ -1100,8 +1088,6 @@
                      END DO
                      ! ji+2,jj+2
-                     ijpt0 = MAX( 1, MIN( ij+2,jpj  ) )
-                     ijpt1 = MAX( 0, MIN( ij+2,jpj-1) )
-                     DO jj = ijpt0, ijpt1
-                        DO ji = iipt0, iipt1
+                     DO jj = mj0(ijrnf(jl,jr) + 2), mj1(ijrnf(jl,jr) + 2)
+                        DO ji = mi0(iirnf(jl,jr) + 1), mi1(iirnf(jl,jr) + 1)
                            DO jn = 1, 3  
                               zcoefr(ji,jj,jn) = MAX( zrup3, zcoefr(ji,jj,jn) )
@@ -1110,8 +1096,6 @@
                      END DO
                      ! ji+2,jj-1
-                     ijpt0 = MAX( 1, MIN( ij-1,jpj  ) )
-                     ijpt1 = MAX( 0, MIN( ij-1,jpj-1) )
-                     DO jj = ijpt0, ijpt1
-                        DO ji = iipt0, iipt1
+                     DO jj = mj0(ijrnf(jl,jr) - 1), mj1(ijrnf(jl,jr) - 1)
+                        DO ji = mi0(iirnf(jl,jr) + 2), mi1(iirnf(jl,jr) + 2)
                            DO jn = 1, 3
                               zcoefr(ji,jj,jn) = MAX( zrup3, zcoefr(ji,jj,jn) )
@@ -1120,8 +1104,6 @@
                      END DO
                      ! ji+2,jj+1
-                     ijpt0 = MAX( 1, MIN( ij+1,jpj  ) )
-                     ijpt1 = MAX( 0, MIN( ij+1,jpj-1) )
-                     DO jj = ijpt0, ijpt1
-                        DO ji = iipt0, iipt1
+                     DO jj = mj0(ijrnf(jl,jr) + 1), mj1(ijrnf(jl,jr) + 1)
+                        DO ji = mi0(iirnf(jl,jr) + 2), mi1(iirnf(jl,jr) + 2)
                            DO jn = 1, 3 
                               zcoefr(ji,jj,jn) = MAX( zrup3, zcoefr(ji,jj,jn) )
@@ -1130,8 +1112,6 @@
                      END DO
                      ! ji+2,jj
-                     ijpt0 = MAX( 1, MIN( ij,jpj  ) )
-                     ijpt1 = MAX( 0, MIN( ij,jpj-1) )
-                     DO jj = ijpt0, ijpt1
-                        DO ji = iipt0, iipt1
+                     DO jj = mj0(ijrnf(jl,jr)    ), mj1(ijrnf(jl,jr)    )
+                        DO ji = mi0(iirnf(jl,jr) + 2), mi1(iirnf(jl,jr) + 2)
                            DO jn = 1, 3
                               zcoefr(ji,jj,jn) = MAX( zrup3, zcoefr(ji,jj,jn) )
@@ -1140,8 +1120,6 @@
                      END DO
                      ! ji-2,jj
-                     iipt0 = MAX( 1, MIN( ii-2,jpi  ) )
-                     iipt1 = MAX( 0, MIN( ii-2,jpi-1) )
-                     DO jj = ijpt0, ijpt1
-                        DO ji = iipt0, iipt1
+                     DO jj = mj0(ijrnf(jl,jr)    ), mj1(ijrnf(jl,jr)    )
+                        DO ji = mi0(iirnf(jl,jr) - 2), mi1(iirnf(jl,jr) - 2)
                            DO jn = 1, 3        
                               zcoefr(ji,jj,jn) = MAX( zrup3, zcoefr(ji,jj,jn) )
@@ -1150,8 +1128,6 @@
                      END DO
                      ! ji-2,jj-2
-                     ijpt0 = MAX( 1, MIN( ij-2,jpj  ) )
-                     ijpt1 = MAX( 0, MIN( ij-2,jpj-1) )
-                     DO jj = ijpt0, ijpt1
-                        DO ji = iipt0, iipt1
+                     DO jj = mj0(ijrnf(jl,jr) - 2), mj1(ijrnf(jl,jr) - 2)
+                        DO ji = mi0(iirnf(jl,jr) - 2), mi1(iirnf(jl,jr) - 2)
                            DO jn = 1, 3 
                               zcoefr(ji,jj,jn) = MAX( zrup3, zcoefr(ji,jj,jn) )
@@ -1160,8 +1136,6 @@
                      END DO
                      ! ji-2,jj+2
-                     ijpt0 = MAX( 1, MIN( ij+2,jpj  ) )
-                     ijpt1 = MAX( 0, MIN( ij+2,jpj-1) )
-                     DO jj = ijpt0, ijpt1
-                        DO ji = iipt0, iipt1
+                     DO jj = mj0(ijrnf(jl,jr) + 2), mj1(ijrnf(jl,jr) + 2)
+                        DO ji = mi0(iirnf(jl,jr) - 2), mi1(iirnf(jl,jr) - 2)
                            DO jn = 1, 3
                               zcoefr(ji,jj,jn) = MAX( zrup3, zcoefr(ji,jj,jn) )
@@ -1170,8 +1144,6 @@
                      END DO
                      ! ji-2,jj-1
-                     ijpt0 = MAX( 1, MIN( ij-1,jpj  ) )
-                     ijpt1 = MAX( 0, MIN( ij-1,jpj-1) )
-                     DO jj = ijpt0, ijpt1
-                        DO ji = iipt0, iipt1
+                     DO jj = mj0(ijrnf(jl,jr) - 1), mj1(ijrnf(jl,jr) - 1)
+                        DO ji = mi0(iirnf(jl,jr) - 2), mi1(iirnf(jl,jr) - 2)
                            DO jn = 1, 3 
                               zcoefr(ji,jj,jn) = MAX( zrup3, zcoefr(ji,jj,jn) )
@@ -1180,8 +1152,6 @@
                      END DO
                      ! ji-2,jj+1
-                     ijpt0 = MAX( 1, MIN( ij+1,jpj  ) )
-                     ijpt1 = MAX( 0, MIN( ij+1,jpj-1) )
-                     DO jj = ijpt0, ijpt1
-                        DO ji = iipt0, iipt1
+                     DO jj = mj0(ijrnf(jl,jr) + 1), mj1(ijrnf(jl,jr) + 1)
+                        DO ji = mi0(iirnf(jl,jr) - 2), mi1(iirnf(jl,jr) - 2)
                            DO jn = 1, 3 
                               zcoefr(ji,jj,jn) = MAX( zrup3, zcoefr(ji,jj,jn) )
Index: tags/nemo_dev_x3/NEMO/OPA_SRC/SBC/ocesbc.F90
===================================================================
--- branches/nemo/NEMO/OPA_SRC/SBC/ocesbc.F90	(revision 4)
+++ tags/nemo_dev_x3/NEMO/OPA_SRC/SBC/ocesbc.F90	(revision 105)
@@ -6,5 +6,6 @@
 
    !!----------------------------------------------------------------------
-   !!    oce_sbc : initialization and namelist read
+   !!   oce_sbc     : ???
+   !!   oce_sbc_dmp : ???
    !!----------------------------------------------------------------------
    !! * Modules used
@@ -14,4 +15,5 @@
    USE ice_oce
    USE blk_oce
+   USE flx_oce
    USE phycst         ! Define parameters for the routines
    USE taumod
@@ -23,4 +25,5 @@
    USE ocfzpt
    USE lbclnk
+   USE lib_mpp
    USE in_out_manager ! I/O manager
 
@@ -32,18 +35,18 @@
 
    !! * Shared module variables
-   REAL(wp), PUBLIC ::   &
-      aplus, aminus,     &
-      empold = 0.e0  ! current year freshwater budget correction
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   &
-      qt  ,         &  ! total surface heat flux (w/m2)
-      q   ,         &  ! surface heat flux (w/m2)
-      qsr ,         &  ! solar radiation (w/m2)
-      emp ,         &  ! evaporation minus precipitation (kg/m2/s = mm/s)
-      emps,         &  ! evaporation - precipitation (free surface)
-      qrp ,         &  ! heat flux damping (w/m2)
-      erp              ! evaporation damping (kg/m2/s = mm/s)
+   REAL(wp), PUBLIC ::   &  !:
+      aplus, aminus,     &  !:
+      empold = 0.e0         !: current year freshwater budget correction
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   &  !:
+      qt  ,         &  !: total surface heat flux (w/m2)
+      q   ,         &  !: surface heat flux (w/m2)
+      qsr ,         &  !: solar radiation (w/m2)
+      emp ,         &  !: evaporation minus precipitation (kg/m2/s = mm/s)
+      emps,         &  !: evaporation - precipitation (free surface)
+      qrp ,         &  !: heat flux damping (w/m2)
+      erp              !: evaporation damping (kg/m2/s = mm/s)
 #if defined key_dynspg_fsc
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   &
-      dmp              ! internal dampind term
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   &  !:
+      dmp              !: internal dampind term
 #endif
 
@@ -156,10 +159,10 @@
 
          ! Re-initialization of fluxes
-         sst_io(:,:) = 0.0
-         sss_io(:,:) = 0.0
-         u_io  (:,:) = 0.0
-         v_io  (:,:) = 0.0
-         gtaux (:,:) = 0.
-         gtauy (:,:) = 0.
+         sst_io(:,:) = 0.e0
+         sss_io(:,:) = 0.e0
+         u_io  (:,:) = 0.e0
+         v_io  (:,:) = 0.e0
+         gtaux (:,:) = 0.e0
+         gtauy (:,:) = 0.e0
 
       ENDIF
@@ -188,6 +191,6 @@
 
       !! * Local declarations
-      INTEGER ::   ji, jj                   ! dummy loop indices
-      REAL ztx, ztaux, zty, ztauy
+      INTEGER  ::   ji, jj                   ! dummy loop indices
+      REAL(wp) ::   ztx, ztaux, zty, ztauy
       !!----------------------------------------------------------------------
 
@@ -255,8 +258,8 @@
 
          ! Re-initialization of fluxes
-         sst_io(:,:) = 0.0
-         sss_io(:,:) = 0.0
-         u_io  (:,:) = 0.0
-         v_io  (:,:) = 0.0
+         sst_io(:,:) = 0.e0
+         sss_io(:,:) = 0.e0
+         u_io  (:,:) = 0.e0
+         v_io  (:,:) = 0.e0
 
       ENDIF
@@ -560,5 +563,5 @@
       !!              
       !! ** Purpose :   provide the thermohaline fluxes (heat and freshwater)
-      !!      to the ocean at each time step.
+      !!                to the ocean at each time step.
       !!
       !! ** Method  :   Constant surface fluxes (read in namelist (namflx))
@@ -569,7 +572,9 @@
       !!        !  91-03  ()  Original code
       !!   8.5  !  02-09  (G. Madec)  F90: Free form and module
+      !!   9.0  !  04-05  (A. Koch-Larrouy) Add Gyre configuration 
       !!----------------------------------------------------------------------
       !! * Modules used
       USE flxrnf                       ! ocean runoffs
+      USE daymod, ONLY : nyear         ! calendar
 
       !! * arguments
@@ -581,42 +586,163 @@
          qsr0 = 0.e0,               &  ! solar heat flux
          emp0 = 0.e0                   ! net freshwater flux
+      REAL(wp) ::   ztrp, zemp_S, zemp_N, zemp_sais, zTstar, zcos_sais, zconv
+      REAL(wp) ::           &
+         zsumemp,           &          ! tampon used for the emp sum
+         zsurf,             &          ! tampon used for the domain sum
+         ztime,             &          ! time in hour
+         ztimemax, ztimemin            ! 21th june, and 21th december if date0 = 1st january
+      REAL(wp), DIMENSION(jpi,jpj) :: t_star
+      INTEGER  ::   ji, jj,  &         ! dummy loop indices
+         js                            ! indice for months
+      INTEGER  ::           &
+         zyear0,            &          ! initial year
+         zmonth0,           &          ! initial month
+         zday0,             &          ! initial day
+         zday_year0,        &          ! initial day since january 1st
+         zdaymax
 
       NAMELIST/namflx/ q0, qsr0, emp0
       !!---------------------------------------------------------------------
 
-      IF( kt == nit000 ) THEN
-
-         ! Read Namelist namflx : surface thermohaline fluxes
-         ! --------------------
-         REWIND ( numnam )
-         READ   ( numnam, namflx )
-
-         IF(lwp) WRITE(numout,*)' '
-         IF(lwp) WRITE(numout,*)' ocesbc  : Constant surface fluxes read in namelist'
-         IF(lwp) WRITE(numout,*)' ~~~~~~~ '
-         IF(lwp) WRITE(numout,*)'           Namelist namflx: set the constant flux values'
-         IF(lwp) WRITE(numout,*)'              net heat flux          q0   = ', q0  , ' W/m2'
-         IF(lwp) WRITE(numout,*)'              solar heat flux        qsr0 = ', qsr0, ' W/m2'
-         IF(lwp) WRITE(numout,*)'              net heat flux          emp0 = ', emp0, ' W/m2'
-
-         qt    (:,:) = q0
-         qsr   (:,:) = qsr0
-         q     (:,:) = q0
-         emp   (:,:) = emp0
-         emps  (:,:) = emp0
-         qrp   (:,:) = 0.e0
-         erp   (:,:) = 0.e0
-
-         runoff(:,:) = 0.e0
-     ENDIF
-
-  END SUBROUTINE oce_sbc
+      !same temperature, E-P as in HAZELEGER 2000
+
+      IF( cp_cfg == 'gyre' ) THEN
+
+          zyear0  = ndate0 / 10000
+          zmonth0 = ( ndate0 - zyear0 * 10000 ) / 100
+          zday0   = ndate0 - zyear0 * 10000 - zmonth0 * 100
+          !Calculates nday_year, day since january 1st
+          zday_year0 = zday0
+          !accumulates days of previous months of this year
+
+          DO js = 1, zmonth0
+            IF(nleapy > 1) THEN
+                zday_year0 = zday_year0 + nleapy
+            ELSE
+                IF( MOD(zyear0, 4 ) == 0 ) THEN
+                    zday_year0 = zday_year0 + nbiss(js)
+                ELSE
+                    zday_year0 = zday_year0 + nobis(js)
+                ENDIF
+            ENDIF
+          END DO
+          ! day (in hours) since january the 1st
+          ztime = FLOAT( kt ) * rdt / (rmmss * rhhmm)  &  ! incrementation in hour
+             &     - (nyear - 1) * rjjhh * raajj       &  !  - nber of hours the precedent years
+             &     + zday_year0 / 24                      ! nber of hours initial date
+          ! day 21th counted since the 1st January
+          zdaymax = 21                                    ! 21th day of the month
+          DO js = 1, 5                                    ! count each day  until end May
+            IF(nleapy > 1) THEN
+                zdaymax = zdaymax + nleapy
+            ELSE
+                IF( MOD(zyear0, 4 ) == 0 ) THEN
+                    zdaymax = zdaymax + nbiss(js)
+                ELSE
+                    zdaymax = zdaymax + nobis(js)
+                ENDIF
+            ENDIF
+          END DO
+          ! 21th june in hours
+          ztimemax = zdaymax * 24
+          ! 21th december day in hours
+          ztimemin = ztimemax + rjjhh * raajj / 2         ! rjjhh * raajj / 4 = 1 seasonal cycle in hours
+          ! amplitudes
+          zemp_S   = 0.7       ! intensity of COS in the South
+          zemp_N   = 0.8       ! intensity of COS in the North
+          zemp_sais= 0.1
+          zTstar   = 28.3      ! intemsity from 28.3 à -5°
+          zcos_sais = COS( (ztime - ztimemax) / (ztimemin - ztimemax) * rpi )         ! 1/2 period between 21th June and 21th December
+          ztrp= - 40.          ! retroaction term (W/m2/K)
+          zconv = 3.16e-5      ! convert 1m/yr->3.16e-5mm/s
+          DO jj = 1, jpj
+            DO ji = 1, jpi
+              t_star (ji,jj) = zTstar * ( 1 + 1. / 50. * zcos_sais )   &                  ! domain from 15° to 50°
+                 &             * COS( rpi * (gphit(ji,jj) - 5.)        &                  ! between 27 and 28 °C at 15N, -3 
+                 &                        / (53.5 * ( 1 + 11 / 53.5 * zcos_sais ) * 2.) ) ! and 13°C at 50N 53.5 + or - 11 
+                 !                                                                        ! = 1/4 period :
+                 !                                                                        ! 64.5 in summer, 42.5 in winter
+              qt  (ji,jj) = ztrp * ( tb(ji,jj,1) - t_star(ji,jj) )
+              IF( gphit(ji,jj) >= 14.845 .AND. 37.2 >= gphit(ji,jj)) THEN
+                 emp  (ji,jj) =   zemp_S * zconv   &
+                    &         * SIN( rpi / 2 * (gphit(ji,jj) - 37.2) / (24.6 - 37.2) )  & ! zero at 37.8°, max at 24.6°
+                    &         * ( 1 - zemp_sais / zemp_S * zcos_sais)
+                 emps (ji,jj) = emp  (ji,jj)
+              ELSE
+                 emp (ji,jj) =  - zemp_N * zconv   &
+                    &         * SIN( rpi / 2 * (gphit(ji,jj) - 37.2) / (46.8 - 37.2) )  & ! zero at 37.8°, max at 46.8°
+                    &         * ( 1 - zemp_sais / zemp_N * zcos_sais )
+                  emps (ji,jj) = emp  (ji,jj)
+              ENDIF
+              qsr (ji,jj) =  230 * COS( 3.1415 * ( gphit(ji,jj) - 23.5 * zcos_sais ) / ( 0.9 * 180 ) )   ! 23.5° : tropics
+            END DO
+          END DO
+          ! compute the emp flux such as its integration on the whole domain and at each time be zero
+          zsumemp = 0.
+          zsurf = 0.
+          DO jj = 1, jpj
+            DO ji = 1, jpi
+              zsumemp = zsumemp + emp(ji, jj) * tmask(ji, jj,  1)
+              zsurf =  zsurf +  tmask(ji, jj,  1)
+            END DO
+          END DO
+
+          IF( lk_mpp )   CALL mpp_sum( zsumemp  )       ! sum over the global domain
+          IF( lk_mpp )   CALL mpp_sum( zsurf    )       ! sum over the global domain
+
+          IF( nbench /= 0 ) THEN
+             ! Benchmark GYRE configuration (to allow the bit to bit comparison between Mpp/Mono case)
+             zsumemp = 0.e0
+          ELSE
+             ! Default GYRE configuration
+             zsumemp = zsumemp / zsurf
+          ENDIF
+          DO jj = 1, jpj
+            DO ji = 1, jpi
+              emp(ji, jj)=  emp(ji, jj) - zsumemp * tmask(ji, jj,  1)
+            END DO
+          END DO
+
+      ELSE
+
+         IF( kt == nit000 ) THEN
+
+            ! Read Namelist namflx : surface thermohaline fluxes
+            ! --------------------
+            REWIND ( numnam )
+            READ   ( numnam, namflx )
+   
+            IF(lwp) WRITE(numout,*)' '
+            IF(lwp) WRITE(numout,*)' ocesbc  : Constant surface fluxes read in namelist'
+            IF(lwp) WRITE(numout,*)' ~~~~~~~ '
+            IF(lwp) WRITE(numout,*)'           Namelist namflx: set the constant flux values'
+            IF(lwp) WRITE(numout,*)'              net heat flux          q0   = ', q0  , ' W/m2'
+            IF(lwp) WRITE(numout,*)'              solar heat flux        qsr0 = ', qsr0, ' W/m2'
+            IF(lwp) WRITE(numout,*)'              net heat flux          emp0 = ', emp0, ' W/m2'
+
+            qt    (:,:) = q0
+            qsr   (:,:) = qsr0
+            q     (:,:) = q0
+            emp   (:,:) = emp0
+            emps  (:,:) = emp0
+            qrp   (:,:) = 0.e0
+            erp   (:,:) = 0.e0
+   
+            runoff(:,:) = 0.e0
+         ENDIF
+      ENDIF
+
+   END SUBROUTINE oce_sbc
 
 # endif
 #endif
 
+#if defined key_dtasal
+   !!----------------------------------------------------------------------
+   !!   'key_dtasal'                                          salinity data
+   !!----------------------------------------------------------------------
    SUBROUTINE oce_sbc_dmp
       !!---------------------------------------------------------------------
-      !!                   ***    ROUTINE oce_sbc_dmp ***
+      !!                   ***  ROUTINE oce_sbc_dmp  ***
       !!                    
       !! ** Purpose : Computation of internal and evaporation damping terms 
@@ -624,15 +750,14 @@
       !!
       !! History :
-      !!   1.0  !  99-11  (M. Imbard)  Original code
-      !!        !  01-03  (D. Ludicone, E. Durand, G. Madec) free surf.
-      !!   2.0  !  02-09  (G. Madec, C. Ethe)  F90: Free form and module
+      !!   9.0  !  04-01  (G. Madec, C. Ethe)  Original code
       !!----------------------------------------------------------------------
       !! * Local declarations
       INTEGER ::   ji, jj                   ! dummy loop indices
       REAL(wp), DIMENSION(jpi,jpj)  :: zsss, zfreeze
+      REAL(wp) ::   zerp, ztrp, zsrp
 #if defined key_dynspg_fsc
-      REAL zwei, zerp, ztrp, zsrp
-      REAL zerpplus(jpi,jpj), zerpminus(jpi,jpj)
-      REAL zplus, zminus, zadefi
+      REAL(wp) ::   zwei
+      REAL(wp) ::   zerpplus(jpi,jpj), zerpminus(jpi,jpj)
+      REAL(wp) ::   zplus, zminus, zadefi
 # if defined key_tradmp
       INTEGER jk
@@ -646,5 +771,5 @@
       DO jj = 1, jpj
          DO ji = 1, jpi
-            freezn(ji,jj) = MAX(0., SIGN(1., freeze(ji,jj)-rsmall))
+            freezn(ji,jj) = MAX(0., SIGN(1., freeze(ji,jj)-rsmall) )
          END DO
       END DO
@@ -678,5 +803,6 @@
       END DO
       ! volume flux associated to internal damping to climatology
-      dmp(:,:) = zstrdmp(:,:) * rauw / ( zsss(:,:) + rsmall )
+!!ibu dmp(:,:) = zstrdmp(:,:) * rauw / ( zsss(:,:) + rsmall )
+      dmp(:,:) = zstrdmp(:,:) * rauw / ( zsss(:,:) + 1.e-20 )
 # else
       dmp(:,:) = 0.e0            ! No internal damping
@@ -693,5 +819,6 @@
             zerp = ( 1. - 2.*upsrnfh(ji,jj) ) * zsrp   &
                & * ( zsss(ji,jj) - s_dta(ji,jj,1) )     &
-               & / ( zsss(ji,jj) + rsmall        )
+               & / ( zsss(ji,jj) + 1.e-20        )
+!ib            & / ( zsss(ji,jj) + rsmall        )
             
             zerp = MIN( zerp, zplus  )
@@ -712,12 +839,10 @@
          END DO
       END DO
-# if defined key_mpp
-      CALL mpp_sum( aplus  )   ! mpp: sum over all the global domain
-      CALL mpp_sum( aminus )
-# endif
-      IF( l_ctl .AND. lwp ) WRITE(numout,*) ' oce_sbc_dmp : a+ = ', aplus, ' a- = ', aminus
+      IF( lk_mpp )   CALL mpp_sum( aplus  )   ! sums over the global domain
+      IF( lk_mpp )   CALL mpp_sum( aminus )
+      IF(l_ctl)   WRITE(numout,*) ' oce_sbc_dmp : a+ = ', aplus, ' a- = ', aminus
 
       zadefi = MIN( aplus, aminus )
-      IF( zadefi == 0.0 ) THEN 
+      IF( zadefi == 0.e0 ) THEN 
          erp(:,:) = 0.e0
       ELSE
@@ -729,8 +854,18 @@
       erp(:,:) = ( 1. - zfreeze(:,:) ) * zsrp    &   ! surface restoring term
          &     * ( zsss(:,:) - s_dta(:,:,1) )     &
-         &     / ( zsss(:,:) + rsmall      )
+         &     / ( zsss(:,:) + 1.e-20      )
+!ib      &     / ( zsss(:,:) + rsmall      )
 #endif
 
    END SUBROUTINE oce_sbc_dmp
+
+#else
+   !!----------------------------------------------------------------------
+   !!   Dummy routine                                      NO salinity data
+   !!----------------------------------------------------------------------
+   SUBROUTINE oce_sbc_dmp         ! Dummy routine
+      WRITE(*,*) 'oce_sbc_dmp: you should not have seen that print! error?'
+   END SUBROUTINE oce_sbc_dmp
+#endif
 
    !!======================================================================
Index: tags/nemo_dev_x3/NEMO/OPA_SRC/SBC/taumod.F90
===================================================================
--- branches/nemo/NEMO/OPA_SRC/SBC/taumod.F90	(revision 4)
+++ tags/nemo_dev_x3/NEMO/OPA_SRC/SBC/taumod.F90	(revision 105)
@@ -13,4 +13,5 @@
    USE in_out_manager  ! I/O manager
    USE daymod          ! calendar
+   USE lbclnk          ! 
 
    IMPLICIT NONE
@@ -22,7 +23,7 @@
    !! * Share modules variables
    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   &
-      taux, tauy,      &  ! surface stress components in (i,j) referential
-      tauxg, tauyg        ! surface stress components in geographical
-      !                   ! referential (used in output)
+      taux, tauy,      &  !: surface stress components in (i,j) referential
+      tauxg, tauyg        !: surface stress components in geographical
+      !                   !  referential (used in output)
    !!----------------------------------------------------------------------
    !!   OPA 9.0 , LODYC-IPSL  (2003)
@@ -105,39 +106,137 @@
       !!----------------------------------------------------------------------
       !! * Arguments
-      INTEGER, INTENT( in  ) ::   kt   ! ocean time step
+      INTEGER, INTENT( in  ) ::   kt    ! ocean time step
+      REAL(wp) ::   ztau, ztau_sais, &  ! wind intensity and of the seasonal cycle
+         ztime,                      &  ! time in hour
+         ztimemax, ztimemin,         &  ! 21th June, and 21th decem. if date0 = 1st january
+         ztaun                          ! intensity 
+      INTEGER  ::   ji, jj, &           ! dummy loop indices
+         js                             ! indice for months
+      INTEGER  ::           &
+         zyear0,            &           ! initial year
+         zmonth0,           &           ! initial month
+         zday0,             &           ! initial day
+         zday_year0,        &           ! initial day since january 1st
+         zdaymax
 
       !! * Local declarations
-      REAL(wp) ::   zfacto             ! 
+      REAL(wp) ::   zfacto              ! 
 
       NAMELIST/namtau/ ntau000, tau0x, tau0y
       !!---------------------------------------------------------------------
 
-      IF( kt == nit000 ) THEN
-
-         ! Read Namelist namtau : surface wind stress
-         ! --------------------
-         REWIND ( numnam )
-         READ   ( numnam, namtau )
-
-         IF(lwp) WRITE(numout,*)' '
-         IF(lwp) WRITE(numout,*)' tau     : Constant surface wind stress read in namelist'
-         IF(lwp) WRITE(numout,*)' ~~~~~~~ '
-         IF(lwp) WRITE(numout,*)'           Namelist namtau: set the constant stress values'
-         IF(lwp) WRITE(numout,*)'              spin up of the stress  ntau000 = ', ntau000, ' time-steps'
-         IF(lwp) WRITE(numout,*)'              constant i-stress      tau0x   = ', tau0x  , ' N/m2'
-         IF(lwp) WRITE(numout,*)'              constant j-stress      tau0y   = ', tau0y  , ' N/m2'
-
-         ntau000 = MAX( ntau000, 1 )   ! must be >= 1
-
-      ENDIF
-
-      ! Increase the surface stress to its nominal value in ntau000 time-step
-      
-      IF( kt <= ntau000 ) THEN
-         zfacto = 0.5 * (  1. - COS( rpi * FLOAT( kt ) / FLOAT( ntau000 ) )  )
-         taux (:,:) = zfacto * tau0x
-         tauy (:,:) = zfacto * tau0y
-         tauxg(:,:) = zfacto * tau0x
-         tauyg(:,:) = zfacto * tau0y
+      IF( cp_cfg == 'gyre' ) THEN
+
+         ! same wind as in Wico
+         !test date0 : ndate0 = 010203
+         zyear0  = ndate0 / 10000
+         zmonth0 = ( ndate0 - zyear0 * 10000 ) / 100
+         zday0   = ndate0 - zyear0 * 10000 - zmonth0 * 100
+         !Calculates nday_year, day since january 1st
+         zday_year0 = zday0
+         !accumulates days of previous months of this year
+
+         DO js = 1, zmonth0
+            IF(nleapy > 1) THEN
+               zday_year0 = zday_year0 + nleapy
+            ELSE
+               IF( MOD(zyear0, 4 ) == 0 ) THEN
+                  zday_year0 = zday_year0 + nbiss(js)
+               ELSE
+                  zday_year0 = zday_year0 + nobis(js)
+               ENDIF
+            ENDIF
+         END DO
+   
+         ! day (in hours) since january the 1st
+         ztime = FLOAT( kt ) * rdt / (rmmss * rhhmm)  &  ! incrementation in hour
+            &     - (nyear - 1) * rjjhh * raajj       &  !  - nber of hours the precedent years
+            &     + zday_year0 / 24                      ! nber of hours initial date
+         ! day 21th counted since the 1st January
+         zdaymax = 21                                    ! 21th day of the month
+         DO js = 1, 5                                    ! count each day  until end May
+           IF(nleapy > 1) THEN
+               zdaymax = zdaymax + nleapy
+           ELSE
+               IF( MOD(zyear0, 4 ) == 0 ) THEN
+                   zdaymax = zdaymax + nbiss(js)
+               ELSE
+                   zdaymax = zdaymax + nobis(js)
+               ENDIF
+           ENDIF
+         END DO
+        ! 21th june in hours
+         ztimemax = zdaymax * 24
+         ! 21th december day in hours
+         ztimemin = ztimemax + rjjhh * raajj / 2  ! rjjhh * raajj / 4 = 1 seasonal cycle in hours
+   
+         ztau = 0.105 / SQRT(2.)            ! mean intensity at 0.105 ; srqt(2) because projected with 45° angle
+         ztau_sais = 0.015                  ! seasonal oscillation intensity
+         ztaun = ztau - ztau_sais * COS( (ztime - ztimemax) / (ztimemin - ztimemax) * rpi )
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+              ! domain from 15° to 50° and 1/2 period along 14° so 5/4 of half period with seasonal cycle
+              taux (ji,jj) = - ztaun * SIN( rpi * (gphiu(ji,jj) - 15.) / (29.-15.) )
+              tauy (ji,jj) =   ztaun * SIN( rpi * (gphiv(ji,jj) - 15.) / (29.-15.) )
+            END DO
+         END DO
+   
+         IF( kt == nit000 ) THEN
+            IF(lwp) WRITE(numout,*)' tau     : Constant surface wind stress read in namelist'
+            IF(lwp) WRITE(numout,*)' ~~~~~~~ '
+            IF(lwp) WRITE(numout,*)'nyear      = ', nyear
+            IF(lwp) WRITE(numout,*)'nmonth     = ', nmonth
+            IF(lwp) WRITE(numout,*)'nday       = ', nday
+            IF(lwp) WRITE(numout,*)'nday_year  = ',nday_year
+            IF(lwp) WRITE(numout,*)'ndastp     = ',ndastp
+            IF(lwp) WRITE(numout,*)'adatrj     = ',adatrj
+            IF(lwp) WRITE(numout,*)'ztime      = ',ztime
+            IF(lwp) WRITE(numout,*)'zdaymax    = ',zdaymax
+   
+            IF(lwp) WRITE(numout,*)'ztimemax   = ',ztimemax
+            IF(lwp) WRITE(numout,*)'ztimemin   = ',ztimemin
+            IF(lwp) WRITE(numout,*)'zyear0     = ', zyear0
+            IF(lwp) WRITE(numout,*)'zmonth0    = ', zmonth0
+            IF(lwp) WRITE(numout,*)'zday0      = ', zday0
+            IF(lwp) WRITE(numout,*)'zday_year0 = ',zday_year0
+            IF(lwp) WRITE(numout,*)'nobis(2)', nobis(2)
+            IF(lwp) WRITE(numout,*)'nobis(5)', nobis(5)
+            IF(lwp) WRITE(numout,*)'nobis(6)', nobis(6)
+            IF(lwp) WRITE(numout,*)'nobis(1)', nobis(1)
+            IF(lwp) WRITE(numout,*)'nobis(zmonth0 -1)', nobis(zmonth0 - 1)
+            IF(lwp) WRITE(numout,*)'raajj  = ', raajj
+         ENDIF
+   
+      ELSE
+
+         IF( kt == nit000 ) THEN
+   
+            ! Read Namelist namtau : surface wind stress
+            ! --------------------
+            REWIND ( numnam )
+            READ   ( numnam, namtau )
+   
+            IF(lwp) WRITE(numout,*)' '
+            IF(lwp) WRITE(numout,*)' tau     : Constant surface wind stress read in namelist'
+            IF(lwp) WRITE(numout,*)' ~~~~~~~ '
+            IF(lwp) WRITE(numout,*)'           Namelist namtau: set the constant stress values'
+            IF(lwp) WRITE(numout,*)'              spin up of the stress  ntau000 = ', ntau000, ' time-steps'
+            IF(lwp) WRITE(numout,*)'              constant i-stress      tau0x   = ', tau0x  , ' N/m2'
+            IF(lwp) WRITE(numout,*)'              constant j-stress      tau0y   = ', tau0y  , ' N/m2'
+   
+            ntau000 = MAX( ntau000, 1 )   ! must be >= 1
+   
+         ENDIF
+   
+         ! Increase the surface stress to its nominal value in ntau000 time-step
+         
+         IF( kt <= ntau000 ) THEN
+            zfacto = 0.5 * (  1. - COS( rpi * FLOAT( kt ) / FLOAT( ntau000 ) )  )
+            taux (:,:) = zfacto * tau0x
+            tauy (:,:) = zfacto * tau0y
+            tauxg(:,:) = zfacto * tau0x
+            tauyg(:,:) = zfacto * tau0y
+         ENDIF
+
       ENDIF
       
