Index: /trunk/NEMO/OPA_SRC/DOM/closea.F90
===================================================================
--- /trunk/NEMO/OPA_SRC/DOM/closea.F90	(revision 702)
+++ /trunk/NEMO/OPA_SRC/DOM/closea.F90	(revision 703)
@@ -2,17 +2,22 @@
    !!======================================================================
    !!                       ***  MODULE  closea  ***
-   !! Closed Seas  : 
+   !! Closed Seas  : specific treatments associated with closed seas
    !!======================================================================
+   !! History :   8.2  !  00-05  (O. Marti)  Original code
+   !!             8.5  !  02-06  (E. Durand, G. Madec)  F90
+   !!             9.0  !  06-07  (G. Madec)  add clo_rnf, clo_ups, clo_bat
+   !!----------------------------------------------------------------------
 
    !!----------------------------------------------------------------------
    !!   dom_clo    : modification of the ocean domain for closed seas cases
-   !!   flx_clo    : Special handling of closed seas
-   !!----------------------------------------------------------------------
-   !! * Modules used
+   !!   sbc_clo    : Special handling of closed seas
+   !!   clo_rnf    : set close sea outflows as river mouths (see sbcrnf)
+   !!   clo_ups    : set mixed centered/upstream scheme in closed sea (see traadv_cen2)
+   !!   clo_bat    : set to zero a field over closed sea (see domzrg)
+   !!----------------------------------------------------------------------
    USE oce             ! dynamics and tracers
    USE dom_oce         ! ocean space and time domain
    USE in_out_manager  ! I/O manager
-   USE ocesbc          ! ocean surface boundary conditions (fluxes)
-   USE flxrnf          ! runoffs
+   USE sbc_oce         ! ocean surface boundary conditions
    USE lib_mpp         ! distributed memory computing library
    USE lbclnk          ! ???
@@ -21,32 +26,29 @@
    PRIVATE
 
-   !! * Accessibility
-   PUBLIC dom_clo      ! routine called by dom_init
-   PUBLIC flx_clo      ! routine called by step
-
-   !! * Share module variables
-   INTEGER, PUBLIC, PARAMETER ::   &  !:
-      jpncs   = 4               !: number of closed sea
-   INTEGER, PUBLIC ::          & !!: namclo : closed seas and lakes
-      nclosea =  0                !: = 0 no closed sea or lake
-      !                           !  = 1 closed sea or lake in the domain
-   INTEGER, PUBLIC, DIMENSION (jpncs) ::   &  !:
-      ncstt,           &  !: Type of closed sea
-      ncsi1, ncsj1,    &  !: closed sea limits                                                                 
-      ncsi2, ncsj2,    &  !: 
-      ncsnr               !: number of point where run-off pours
-   INTEGER, PUBLIC, DIMENSION (jpncs,4) ::   &
-      ncsir, ncsjr        !: Location of run-off
-
-   !! * Module variable
-   REAL(wp), DIMENSION (jpncs+1) ::   &
-      surf               ! closed sea surface
+   PUBLIC dom_clo      ! routine called by domain module
+   PUBLIC sbc_clo      ! routine called by step module
+   PUBLIC clo_rnf      ! routine called by sbcrnf module
+   PUBLIC clo_ups      ! routine called in traadv_cen2(_jki) module
+   PUBLIC clo_bat      ! routine called in domzgr module
+
+   !!* Namelist namclo : closed seas and lakes
+   INTEGER, PUBLIC                     ::   nclosea =  0     !: = 0 no closed sea or lake
+      !                                                      !  = 1 closed sea or lake in the domain
+      
+   INTEGER, PUBLIC, PARAMETER          ::   jpncs   = 4      !: number of closed sea
+   INTEGER, PUBLIC, DIMENSION(jpncs)   ::   ncstt            !: Type of closed sea
+   INTEGER, PUBLIC, DIMENSION(jpncs)   ::   ncsi1, ncsj1     !: south-west closed sea limits (i,j)
+   INTEGER, PUBLIC, DIMENSION(jpncs)   ::   ncsi2, ncsj2     !: north-east closed sea limits (i,j)
+   INTEGER, PUBLIC, DIMENSION(jpncs)   ::   ncsnr            !: number of point where run-off pours
+   INTEGER, PUBLIC, DIMENSION(jpncs,4) ::   ncsir, ncsjr     !: Location of runoff
+
+   REAL(wp), DIMENSION (jpncs+1)       ::   surf             ! closed sea surface
 
    !! * Substitutions
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
-   !!  OPA 9.0 , LOCEAN-IPSL (2005) 
+   !!  OPA 9.0 , LOCEAN-IPSL (2006) 
    !! $Id$
-   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 
+   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
 
@@ -60,19 +62,14 @@
       !!
       !! ** Method  :   if a closed sea is located only in a model grid point
-      !!      just the thermodynamic processes are applied.
-      !!
-      !! ** Action :   ncsi1(), ncsj1() : south-west closed sea limits (i,j)
-      !!               ncsi2(), ncsj2() : north-east Closed sea limits (i,j)
-      !!               ncsir(), ncsjr() : Location of runoff
-      !!               ncsnr            : number of point where run-off pours
-      !!               ncstt            : Type of closed sea
-      !!                                  =0 spread over the world ocean
-      !!                                  =2 put at location runoff
-      !!
-      !! History :
-      !!        !  01-04  (E. Durand)  Original code
-      !!   8.5  !  02-06  (G. Madec)  F90: Free form and module
-      !!----------------------------------------------------------------------
-      !! * Local variables
+      !!                just the thermodynamic processes are applied.
+      !!
+      !! ** Action  :   ncsi1(), ncsj1() : south-west closed sea limits (i,j)
+      !!                ncsi2(), ncsj2() : north-east Closed sea limits (i,j)
+      !!                ncsir(), ncsjr() : Location of runoff
+      !!                ncsnr            : number of point where run-off pours
+      !!                ncstt            : Type of closed sea
+      !!                                   =0 spread over the world ocean
+      !!                                   =2 put at location runoff
+      !!----------------------------------------------------------------------
       INTEGER ::   jc            ! dummy loop indices
       !!----------------------------------------------------------------------
@@ -90,10 +87,9 @@
 
       IF( cp_cfg == "orca" ) THEN
-   
+         !
          SELECT CASE ( jp_cfg )
          !                                           ! =======================
          CASE ( 2 )                                  !  ORCA_R2 configuration
             !                                        ! =======================
-
             !                                            ! Caspian Sea
             ncsnr(1)   =   1  ;  ncstt(1)   =   0           ! spread over the globe
@@ -116,9 +112,7 @@
             ncsi2(4)   =   6  ;  ncsj2(4)   = 112
             ncsir(4,1) = 171  ;  ncsjr(4,1) = 106 
-
             !                                        ! =======================
          CASE ( 4 )                                  !  ORCA_R4 configuration
             !                                        ! =======================
-
             !                                            ! Caspian Sea
             ncsnr(1)   =  1  ;  ncstt(1)   =  0  
@@ -144,5 +138,4 @@
             ncsi2(4)   = 76  ;  ncsj2(4)   = 61
             ncsir(4,1) = 84  ;  ncsjr(4,1) = 59 
-
             !                                        ! =======================
          CASE ( 025 )                                ! ORCA_R025 configuration
@@ -157,7 +150,7 @@
             ncsi2(2)   = 1304 ; ncsj2(2)   = 747
             ncsir(2,1) = 1    ; ncsjr(2,1) = 1
-
+            !
          END SELECT
-
+         !
       ENDIF
 
@@ -171,12 +164,11 @@
          ncsj2(jc)   = mj1( ncsj2(jc) )  
       END DO
-         
-
+      !
    END SUBROUTINE dom_clo
 
 
-   SUBROUTINE flx_clo( kt )
-      !!---------------------------------------------------------------------
-      !!                  ***  ROUTINE flx_clo  ***
+   SUBROUTINE sbc_clo( kt )
+      !!---------------------------------------------------------------------
+      !!                  ***  ROUTINE sbc_clo  ***
       !!                    
       !! ** Purpose :   Special handling of closed seas
@@ -186,25 +178,18 @@
       !!      put as run-off in open ocean.
       !!
-      !! ** Action :
-      !!
-      !! History :
-      !!   8.2  !  00-05  (O. Marti)  Original code
-      !!   8.5  !  02-07  (G. Madec)  Free form, F90
-      !!----------------------------------------------------------------------
-      !! * Arguments
-      INTEGER, INTENT (in) :: kt
-
-      !! * Local declarations
-      REAL(wp), DIMENSION (jpncs) :: zemp
-      INTEGER  :: ji, jj, jc, jn
-      REAL(wp) :: zze2
-      !!----------------------------------------------------------------------
-
-      ! 1 - Initialisation
-      ! ------------------
-
-      IF( kt == nit000 ) THEN 
+      !! ** Action  :   emp, emps   updated surface freshwater fluxes at kt
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT(in) ::   kt   ! ocean model time step
+      !
+      INTEGER                     ::   ji, jj, jc, jn   ! dummy loop indices
+      REAL(wp)                    ::   zze2
+      REAL(wp), DIMENSION (jpncs) ::   zemp
+      !!----------------------------------------------------------------------
+      !
+      !                                                   !------------------!
+      IF( kt == nit000 ) THEN                             !  Initialisation  !
+         !                                                !------------------!
          IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*)'flx_clo : closed seas '
+         IF(lwp) WRITE(numout,*)'sbc_clo : closed seas '
          IF(lwp) WRITE(numout,*)'~~~~~~~'
 
@@ -216,18 +201,7 @@
             DO jj = ncsj1(jc), ncsj2(jc)
                DO ji = ncsi1(jc), ncsi2(jc)
-                  ! surface of closed seas
-                  surf(jc) = surf(jc) + e1t(ji,jj)*e2t(ji,jj)*tmask_i(ji,jj)
-                  ! upstream in closed seas
-                  upsadv(ji,jj) = 0.5
+                  surf(jc) = surf(jc) + e1t(ji,jj) * e2t(ji,jj) * tmask_i(ji,jj)      ! surface of closed seas
                END DO 
             END DO 
-            ! upstream at closed sea outflow
-            IF( ncstt(jc) >= 1 ) THEN 
-                DO jn = 1, 4
-                  ji = mi0( ncsir(jc,jn) )
-                  jj = mj0( ncsjr(jc,jn) )
-                  upsrnfh(ji,jj) = MAX( upsrnfh(ji,jj), 1.0 )
-                END DO 
-            ENDIF 
          END DO 
          IF( lk_mpp )   CALL mpp_sum ( surf, jpncs+1 )       ! mpp: sum over all the global domain
@@ -235,6 +209,5 @@
          IF(lwp) WRITE(numout,*)'     Closed sea surfaces'
          DO jc = 1, jpncs
-            IF(lwp) WRITE(numout,FMT='(1I3,4I4,5X,F16.2)')    &
-                jc, ncsi1(jc), ncsi2(jc), ncsj1(jc), ncsj2(jc), surf(jc)
+            IF(lwp)WRITE(numout,FMT='(1I3,4I4,5X,F16.2)') jc, ncsi1(jc), ncsi2(jc), ncsj1(jc), ncsj2(jc), surf(jc)
          END DO
 
@@ -243,11 +216,9 @@
             surf(jpncs+1) = surf(jpncs+1) - surf(jc)
          END DO           
- 
+         !
       ENDIF
-
-      ! 2 - Computation
-      ! ---------------
-      zemp = 0.e0
-
+      !                                                   !--------------------!
+      !                                                   !  update emp, emps  !
+      zemp = 0.e0                                         !--------------------!
       DO jc = 1, jpncs
          DO jj = ncsj1(jc), ncsj2(jc)
@@ -257,5 +228,5 @@
          END DO 
       END DO
-      IF( lk_mpp )   CALL mpp_sum ( zemp , jpncs )       ! mpp: sum over all the global domain
+      IF( lk_mpp )   CALL mpp_sum ( zemp(:) , jpncs )       ! mpp: sum over all the global domain
 
       IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN      ! Black Sea case for ORCA_R2 configuration
@@ -266,5 +237,5 @@
 
       DO jc = 1, jpncs
-
+         !
          IF( ncstt(jc) == 0 ) THEN 
             ! water/evap excess is shared by all open ocean
@@ -303,5 +274,5 @@
             ENDIF 
          ENDIF 
-
+         !
          DO jj = ncsj1(jc), ncsj2(jc)
             DO ji = ncsi1(jc), ncsi2(jc)
@@ -310,14 +281,99 @@
             END DO  
          END DO 
-
+         !
       END DO 
-
-
-      ! 5. Boundary condition on emp and emps
-      ! -------------------------------------
+      !
       CALL lbc_lnk( emp , 'T', 1. )
       CALL lbc_lnk( emps, 'T', 1. )
-
-   END SUBROUTINE flx_clo
+      !
+   END SUBROUTINE sbc_clo
+   
+   
+   SUBROUTINE clo_rnf( p_rnfmsk )
+      !!---------------------------------------------------------------------
+      !!                  ***  ROUTINE sbc_rnf  ***
+      !!                    
+      !! ** Purpose :   allow the treatment of closed sea outflow grid-points
+      !!                to be the same as river mouth grid-points
+      !!
+      !! ** Method  :   set to 1 the runoff mask (mskrnf, see sbcrnf module)
+      !!                at the closed sea outflow grid-point.
+      !!
+      !! ** Action  :   update (p_)mskrnf (set 1 at closed sea outflow)
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   p_rnfmsk   ! river runoff mask (rnfmsk array)
+      !
+      INTEGER  ::   jc, jn      ! dummy loop indices
+      INTEGER  ::   ii, ij      ! temporary integer
+      !!----------------------------------------------------------------------
+      !
+      DO jc = 1, jpncs
+         IF( ncstt(jc) >= 1 ) THEN            ! runoff mask set to 1 at closed sea outflows
+             DO jn = 1, 4
+               ii = mi0( ncsir(jc,jn) )
+               ij = mj0( ncsjr(jc,jn) )
+               p_rnfmsk(ii,ij) = MAX( p_rnfmsk(ii,ij), 1.0 )
+            END DO 
+         ENDIF 
+      END DO 
+      !
+   END SUBROUTINE clo_rnf
+
+   
+   SUBROUTINE clo_ups( p_upsmsk )
+      !!---------------------------------------------------------------------
+      !!                  ***  ROUTINE sbc_rnf  ***
+      !!                    
+      !! ** Purpose :   allow the treatment of closed sea outflow grid-points
+      !!                to be the same as river mouth grid-points
+      !!
+      !! ** Method  :   set to 0.5 the upstream mask (upsmsk, see traadv_cen2 
+      !!                module) over the closed seas.
+      !!
+      !! ** Action  :   update (p_)upsmsk (set 0.5 over closed seas)
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   p_upsmsk   ! upstream mask (upsmsk array)
+      !
+      INTEGER  ::   jc, ji, jj      ! dummy loop indices
+      !!----------------------------------------------------------------------
+      !
+      DO jc = 1, jpncs
+         DO jj = ncsj1(jc), ncsj2(jc)
+            DO ji = ncsi1(jc), ncsi2(jc)
+               p_upsmsk(ji,jj) = 0.5            ! mixed upstream/centered scheme over closed seas
+            END DO 
+         END DO 
+       END DO 
+       !
+   END SUBROUTINE clo_ups
+   
+      
+   SUBROUTINE clo_bat( pbat, kbat )
+      !!---------------------------------------------------------------------
+      !!                  ***  ROUTINE clo_bat  ***
+      !!                    
+      !! ** Purpose :   suppress closed sea from the domain
+      !!
+      !! ** Method  :   set to 0 the meter and level bathymetry (given in 
+      !!                arguments) over the closed seas.
+      !!
+      !! ** Action  :   set pbat=0 and kbat=0 over closed seas
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   pbat   ! bathymetry in meters (bathy array)
+      INTEGER , DIMENSION(jpi,jpj), INTENT(inout) ::   kbat   ! bathymetry in levels (mbathy array)
+      !
+      INTEGER  ::   jc, ji, jj      ! dummy loop indices
+      !!----------------------------------------------------------------------
+      !
+      DO jc = 1, jpncs
+         DO jj = ncsj1(jc), ncsj2(jc)
+            DO ji = ncsi1(jc), ncsi2(jc)
+               pbat(ji,jj) = 0.e0   
+               kbat(ji,jj) = 0   
+            END DO 
+         END DO 
+       END DO 
+       !
+   END SUBROUTINE clo_bat
 
    !!======================================================================
Index: /trunk/NEMO/OPA_SRC/SBC/albedo.F90
===================================================================
--- /trunk/NEMO/OPA_SRC/SBC/albedo.F90	(revision 702)
+++ /trunk/NEMO/OPA_SRC/SBC/albedo.F90	(revision 703)
@@ -4,49 +4,46 @@
    !! Ocean forcing:  bulk thermohaline forcing of the ocean (or ice)
    !!=====================================================================
-   !!----------------------------------------------------------------------
-   !!   flx_blk_albedo : albedo for ocean and ice (clear and overcast skies)
-   !!----------------------------------------------------------------------
-   !! * Modules used
+   !! History :  8.0  !  01-04  (LIM 1.0)
+   !!            8.5  !  03-07  (C. Ethe, G. Madec)  Optimization (old name:shine)
+   !!            9.0  !  04-11  (C. Talandier)  add albedo_init
+   !!            9.0  !  06-08  (G. Madec)  cleaning for surface module
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   blk_albedo  : albedo for ocean and ice (clear and overcast skies)
+   !!   albedo_init : initialisation
+   !!----------------------------------------------------------------------
    USE oce             ! ocean dynamics and tracers
-   USE dom_oce         ! ocean space and time domain
-   USE cpl_oce         ! ???
    USE phycst          ! physical constants
-   USE daymod
-   USE blk_oce         ! bulk variables
-   USE flx_oce         ! forcings variables
-   USE ocfzpt          ! ???
    USE in_out_manager
-   USE lbclnk
 
    IMPLICIT NONE
    PRIVATE
 
-   !! * Accessibility
-   PUBLIC flx_blk_albedo ! routine called by limflx.F90 in coupled
-                         ! and in flxblk.F90 in forced
-   !! * Module variables
-   INTEGER  ::             &  !: nameos : ocean physical parameters
-      albd_init = 0           !: control flag for initialization
-
-   REAL(wp)  ::            &  ! constant values
-      zzero   = 0.e0    ,  &
-      zone    = 1.0
-
-   !! * constants for albedo computation (flx_blk_albedo)
+   PUBLIC   blk_albedo   ! routine called by sbcice_lim module
+
+   INTEGER  ::   albd_init = 0    !: control flag for initialization
+
+   REAL(wp) ::   zzero   = 0.e0   ! constant values
+   REAL(wp) ::   zone    = 1.e0   !    "       "
+
+   REAL(wp) ::   c1     = 0.05    ! constants values
+   REAL(wp) ::   c2     = 0.10    !    "        "
+   REAL(wp) ::   cmue   = 0.40    !  cosine of local solar altitude
+
+   !!* namelist namalb
    REAL(wp) ::   &
-      c1     = 0.05  ,     &   ! constants values
-      c2     = 0.10  ,     &
+      cgren  = 0.06  ,     &   !  correction of the snow or ice albedo to take into account
+      !                        !  effects of cloudiness (Grenfell & Perovich, 1984)
       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
-
-   !!----------------------------------------------------------------------
-   !!   OPA 9.0 , LOCEAN-IPSL (2005) 
+      alphc  = 0.65 
+   NAMELIST/namalb/ cgren, albice, alphd, alphdi, alphc
+
+   !!----------------------------------------------------------------------
+   !!   OPA 9.0 , LOCEAN-IPSL (2006) 
    !! $Id$
-   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 
+   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
 
@@ -58,37 +55,26 @@
    !!----------------------------------------------------------------------
 
-   SUBROUTINE flx_blk_albedo( palb , palcn , palbp , palcnp )
-      !!----------------------------------------------------------------------
-      !!               ***  ROUTINE flx_blk_albedo  ***
+   SUBROUTINE blk_albedo( palb , palcn , palbp , palcnp )
+      !!----------------------------------------------------------------------
+      !!               ***  ROUTINE blk_albedo  ***
       !!          
       !! ** Purpose :   Computation of the albedo of the snow/ice system 
-      !!      as well as the ocean one
+      !!                as well as the ocean one
       !!       
       !! ** Method  : - Computation of the albedo of snow or ice (choose the 
-      !!      rignt one by a large number of tests
+      !!                rignt one by a large number of tests
       !!              - Computation of the albedo of the ocean
       !!
-      !! References :
-      !!      Shine and Hendersson-Sellers 1985, JGR, 90(D1), 2243-2250.
-      !!
-      !! History :
-      !!  8.0   !  01-04  (LIM 1.0)
-      !!  8.5   !  03-07  (C. Ethe, G. Madec)  Optimization (old name:shine)
-      !!----------------------------------------------------------------------
-      !! * Modules used
-      USE ice                   ! ???
-
-      !! * Arguments
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(out) ::  &
-         palb         ,     &    !  albedo of ice under overcast sky
-         palcn        ,     &    !  albedo of ocean under overcast sky
-         palbp        ,     &    !  albedo of ice under clear sky 
-         palcnp                  !  albedo of ocean under clear sky
-
-      !! * Local variables
-      INTEGER ::    &
-         ji, jj                   ! dummy loop indices
-      REAL(wp) ::   & 
-         zmue14         ,     &   !  zmue**1.4
+      !! References :   Shine and Hendersson-Sellers 1985, JGR, 90(D1), 2243-2250.
+      !!----------------------------------------------------------------------
+      USE ice             ! ???
+      !!
+      REAL(wp), INTENT(out), DIMENSION(jpi,jpj) ::   palb     ! albedo of ice under overcast sky
+      REAL(wp), INTENT(out), DIMENSION(jpi,jpj) ::   palcn    ! albedo of ocean under overcast sky
+      REAL(wp), INTENT(out), DIMENSION(jpi,jpj) ::   palbp    ! albedo of ice under clear sky 
+      REAL(wp), INTENT(out), DIMENSION(jpi,jpj) ::   palcnp   ! albedo of ocean under clear sky
+      !!
+      INTEGER  ::   ji, jj                   ! dummy loop indices
+      REAL(wp) ::   zcoef,    &   ! temporary scalar
          zalbpsnm       ,     &   !  albedo of ice under clear sky when snow is melting
          zalbpsnf       ,     &   !  albedo of ice under clear sky when snow is freezing
@@ -99,17 +85,14 @@
          zihsc1         ,     &   !  = 1 hsn <= c1 ; = 0 hsn > c1
          zihsc2                   !  = 1 hsn >= c2 ; = 0 hsn < c2
-      REAL(wp), DIMENSION(jpi,jpj) ::  &
-         zalbfz         ,     &   !  ( = alphdi for freezing ice ; = albice for melting ice )
-         zficeth                  !  function of ice thickness
-      LOGICAL , DIMENSION(jpi,jpj) ::  &
-         llmask
+      LOGICAL , DIMENSION(jpi,jpj) ::   llmask    ! 
+      REAL(wp), DIMENSION(jpi,jpj) ::   zalbfz    ! ( = alphdi for freezing ice ; = albice for melting ice )
+      REAL(wp), DIMENSION(jpi,jpj) ::   zficeth   ! function of ice thickness
       !!---------------------------------------------------------------------
       
-      ! initialization 
-      IF( albd_init == 0 )   CALL albedo_init
-
-      !-------------------------                                                             
+      IF( albd_init == 0 )   CALL albedo_init      ! initialization 
+
+      !---------------------------
       !  Computation of  zficeth
-      !-------------------------- 
+      !---------------------------
       
       llmask = (hsnif == 0.e0) .AND. ( sist >= rt0_ice )
@@ -175,12 +158,9 @@
       !-------------------------- -----------------                                                          
       
-      !  Parameterization of Briegled and Ramanathan, 1982 
-      zmue14      = zmue**1.4                                       
-      palcnp(:,:) = 0.05 / ( 1.1 * zmue14 + 0.15 )                
-      
-      !  Parameterization of Kondratyev, 1969 and Payne, 1972
-      palcn(:,:)  = 0.06                                                 
-      
-   END SUBROUTINE flx_blk_albedo
+      zcoef = 0.05 / ( 1.1 * cmue**1.4 + 0.15 )        ! Parameterization of Briegled and Ramanathan, 1982 
+      palcnp(:,:) = zcoef
+      palcn(:,:)  = 0.06                               ! Parameterization of Kondratyev, 1969 and Payne, 1972
+      !
+   END SUBROUTINE blk_albedo
 
 # else
@@ -189,45 +169,29 @@
    !!----------------------------------------------------------------------
 
-   SUBROUTINE flx_blk_albedo( palb , palcn , palbp , palcnp )
-      !!----------------------------------------------------------------------
-      !!               ***  ROUTINE flx_blk_albedo  ***
+   SUBROUTINE blk_albedo( palb , palcn , palbp , palcnp )
+      !!----------------------------------------------------------------------
+      !!               ***  ROUTINE blk_albedo  ***
       !! 
-      !! ** Purpose :   Computation of the albedo of the snow/ice system
-      !!      as well as the ocean one
-      !!
-      !! ** Method  :   Computation of the albedo of snow or ice (choose the
-      !!      wright one by a large number of tests Computation of the albedo
-      !!      of the ocean
-      !!
-      !! History :
-      !!  8.0   !  01-04  (LIM 1.0)
-      !!  8.5   !  03-07  (C. Ethe, G. Madec)  Optimization (old name:shine)
-      !!----------------------------------------------------------------------
-      !! * Arguments
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(out) ::  &
-         palb         ,     &    !  albedo of ice under overcast sky
-         palcn        ,     &    !  albedo of ocean under overcast sky
-         palbp        ,     &    !  albedo of ice under clear sky
-         palcnp                  !  albedo of ocean under clear sky
-
-      REAL(wp) ::   &
-         zmue14                 !  zmue**1.4
-      !!----------------------------------------------------------------------
-
-      !--------------------------------------------
-      !    Computation of the albedo of the ocean
-      !-------------------------- -----------------
-
-      !  Parameterization of Briegled and Ramanathan, 1982
-      zmue14      = zmue**1.4
-      palcnp(:,:) = 0.05 / ( 1.1 * zmue14 + 0.15 )
-
-      !  Parameterization of Kondratyev, 1969 and Payne, 1972
-      palcn(:,:)  = 0.06
-
-      palb (:,:)  = palcn(:,:)
-      palbp(:,:)  = palcnp(:,:)
-
-   END SUBROUTINE flx_blk_albedo
+      !! ** Purpose :   Computation of the albedo of the ocean
+      !!
+      !! ** Method  :   ....
+      !!----------------------------------------------------------------------
+      REAL(wp), INTENT(out), DIMENSION(jpi,jpj) ::   palb     ! albedo of ice under overcast sky
+      REAL(wp), INTENT(out), DIMENSION(jpi,jpj) ::   palcn    ! albedo of ocean under overcast sky
+      REAL(wp), INTENT(out), DIMENSION(jpi,jpj) ::   palbp    ! albedo of ice under clear sky
+      REAL(wp), INTENT(out), DIMENSION(jpi,jpj) ::   palcnp   ! albedo of ocean under clear sky
+      !!
+      REAL(wp) ::   zcoef    ! temporary scalar
+      !!----------------------------------------------------------------------
+      !
+      zcoef = 0.05 / ( 1.1 * cmue**1.4 + 0.15 )
+
+      palcnp(:,:) = zcoef           ! Parameterization of Briegled and Ramanathan, 1982
+      palcn(:,:)  = 0.06            ! Parameterization of Kondratyev, 1969 and Payne, 1972
+
+      palb (:,:)  = zcoef           ! ice overcast  albedo set to oceanvalue
+      palbp(:,:)  = 0.06            ! ice clear sky albedo set to oceanvalue
+      !
+   END SUBROUTINE blk_albedo
 
 #endif
@@ -240,39 +204,24 @@
       !!
       !! ** Method  :   Read the namelist namalb
-      !!
-      !! ** Action  :  
-      !!
-      !!
-      !! History :
-      !!   9.0  !  04-11  (C. Talandier)  Original code
-      !!----------------------------------------------------------------------
-      NAMELIST/namalb/ cgren, albice, alphd, alphdi, alphc
-      !!----------------------------------------------------------------------
-      !!  OPA 9.0, LODYC-IPSL (2004)
-      !!----------------------------------------------------------------------
-
-      ! set the initialization flag to 1
-      albd_init = 1           ! indicate that the initialization has been done
-
-      ! Read Namelist namalb : albedo parameters
-      REWIND( numnam )
+      !!----------------------------------------------------------------------
+      !
+      albd_init = 1              ! set the initialization flag to 1 (done)
+
+      REWIND( numnam )           ! Read Namelist namalb : albedo parameters
       READ  ( numnam, namalb )
 
-      ! Control print
-      IF(lwp) THEN
+      IF(lwp) THEN               ! Control print
          WRITE(numout,*)
-         WRITE(numout,*) 'albedo_init : albedo '
+         WRITE(numout,*) 'albedo_init : set albedo parameters from namelist namalb'
          WRITE(numout,*) '~~~~~~~~~~~'
-         WRITE(numout,*) '          Namelist namalb : set albedo parameters'
-         WRITE(numout,*)
-         WRITE(numout,*) '             correction of the snow or ice albedo to take into account cgren = ', cgren
-         WRITE(numout,*) '             albedo of melting ice in the arctic and antarctic        albice = ', albice
-         WRITE(numout,*) '             coefficients for linear                                   alphd = ', alphd
-         WRITE(numout,*) '             interpolation used to compute albedo                     alphdi = ', alphdi
-         WRITE(numout,*) '             between two extremes values (Pyane, 1972)                 alphc = ', alphc
-         WRITE(numout,*)
+         WRITE(numout,*) '             correction for snow and ice albedo                    cgren  = ', cgren
+         WRITE(numout,*) '             albedo of melting ice in the arctic and antarctic     albice = ', albice
+         WRITE(numout,*) '             coefficients for linear                               alphd  = ', alphd
+         WRITE(numout,*) '             interpolation used to compute albedo                  alphdi = ', alphdi
+         WRITE(numout,*) '             between two extremes values (Pyane, 1972)             alphc  = ', alphc
       ENDIF
-
+      !
    END SUBROUTINE albedo_init
+
    !!======================================================================
 END MODULE albedo
Index: /trunk/NEMO/OPA_SRC/daymod.F90
===================================================================
--- /trunk/NEMO/OPA_SRC/daymod.F90	(revision 702)
+++ /trunk/NEMO/OPA_SRC/daymod.F90	(revision 703)
@@ -4,9 +4,16 @@
    !! Ocean        :  calendar 
    !!=====================================================================
+   !! History :        !  94-09  (M. Pontaud M. Imbard)  Original code
+   !!                  !  97-03  (O. Marti)
+   !!                  !  97-05  (G. Madec) 
+   !!                  !  97-08  (M. Imbard)
+   !!             9.0  !  03-09  (G. Madec)  F90 + nyear, nmonth, nday
+   !!                  !  04-01  (A.M. Treguier) new calculation based on adatrj
+   !!                  !  06-08  (G. Madec)  surface module major update
+   !!----------------------------------------------------------------------      
 
    !!----------------------------------------------------------------------
    !!   day        : calendar
    !!----------------------------------------------------------------------
-   !! * Modules used
    USE dom_oce         ! ocean space and time domain
    USE phycst          ! physical constants
@@ -17,23 +24,33 @@
    PRIVATE
 
-   !! * Routine accessibility
    PUBLIC day        ! called by step.F90
 
-   !! * Shared module variables
-   INTEGER , PUBLIC ::   &  !:
-      nyear     ,   &  !: current year
-      nmonth    ,   &  !: current month
-      nday      ,   &  !: current day of the month
-      nday_year ,   &  !: curent day counted from jan 1st of the current year
-      ndastp           !: time step date in year/month/day aammjj
-   REAL(wp), PUBLIC ::   &  !:
-       adatrj   ,   &  !: number of elapsed days since the begining of the run
-       adatrj0         !: value of adatrj at nit000-1 (before the present run).
-       !               !  it is the accumulated duration of previous runs
-       !               !  that may have been run with different time steps.
+   INTEGER , PUBLIC ::   nyear       !: current year
+   INTEGER , PUBLIC ::   nmonth      !: current month
+   INTEGER , PUBLIC ::   nday        !: current day of the month
+   INTEGER , PUBLIC ::   nday_year   !: current day counted from jan 1st of the current year
+   REAL(wp), PUBLIC ::   rsec_year   !: current time step counted in second since 00h jan 1st of the current year
+   REAL(wp), PUBLIC ::   rsec_month  !: current time step counted in second since 00h 1st day of the current month
+   REAL(wp), PUBLIC ::   rsec_day    !: current time step counted in second since 00h         of the current day
+   INTEGER , PUBLIC ::   ndastp      !: time step date in year/month/day aammjj
+
+!!gm supprimer adatrj et adatrj0 ==> remplacer par rsecday.....
+   REAL(wp), PUBLIC ::   adatrj      !: number of elapsed days since the begining of the run
+   REAL(wp), PUBLIC ::   adatrj0     !: value of adatrj at nit000-1 (before the present run).
+   !                                 !  it is the accumulated duration of previous runs
+   !                                 !  that may have been run with different time steps.
+   INTEGER , PUBLIC, DIMENSION(0:13) ::   nmonth_len   !: length of the current year
+
+   INTEGER, PUBLIC, DIMENSION(12) ::   nbiss = (/ 31, 29, 31, 30, 31, 30,    &  !: number of days per month
+      &                                           31, 31, 30, 31, 30, 31 /)     !: (leap-year)
+   INTEGER, PUBLIC, DIMENSION(12) ::   nobis = (/ 31, 28, 31, 30, 31, 30,    &  !: number of days per month
+      &                                           31, 31, 30, 31, 30, 31 /)     !: (365 days a year)
+
+   REAL(wp), PUBLIC, DIMENSION(0:14) ::   rmonth_half(0:14)
+
    !!----------------------------------------------------------------------
-   !!  OPA 9.0 , LOCEAN-IPSL (2005) 
+   !!  OPA 9.0 , LOCEAN-IPSL (2006) 
    !! $Id$
-   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
+   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
 
@@ -54,21 +71,10 @@
       !!              - ndastp    : =nyear*10000+nmonth*100+nday
       !!              - adatrj    : date in days since the beginning of the run
-      !!
-      !! History :
-      !!        !  94-09  (M. Pontaud M. Imbard)  Original code
-      !!        !  97-03  (O. Marti)
-      !!        !  97-05  (G. Madec) 
-      !!        !  97-08  (M. Imbard)
-      !!   9.0  !  03-09  (G. Madec)  F90 + nyear, nmonth, nday
-      !!        !  04-01  (A.M. Treguier) new calculation based on adatrj
+      !!              - rsec_year : current time of the year (in second since 00h, jan 1st)
       !!----------------------------------------------------------------------      
-      !! * Arguments
-      INTEGER, INTENT( in ) ::   kt      ! ocean time-step indices
-
-      !! * Local declarations
-      INTEGER  ::   js                   ! dummy loop indice
-      INTEGER  ::   iend, iday0, iday1   ! temporary integers
-      REAL(wp) :: zadatrjn, zadatrjb     ! adatrj at timestep kt-1 and kt-2 
-      CHARACTER (len=25) :: charout
+      INTEGER, INTENT(in) ::   kt        ! ocean time-step indices
+      !
+      INTEGER  ::   js, jm               ! dummy loop indice
+      CHARACTER (len=25) ::   charout
       !!----------------------------------------------------------------------
 
@@ -77,85 +83,110 @@
       !-----------------------------------------------------------------
 
-      IF( kt == nit000 ) THEN
-
+      !                        ! ---------------- !
+      IF( kt == -1 ) THEN      !  Initialisation  !
+         !                     ! ---------------- !
+         !
          IF( .NOT.ln_rstart )   adatrj0 = 0.e0      ! adatrj0 initialized in rst_read when restart 
 
-         adatrj  = adatrj0
+         ! set the calandar from adatrj0 and ndastp (read in restart file and namelist)
+         adatrj  =   adatrj0      !???? bug.... toujours rest   !!gm
          nyear   =   ndastp / 10000
          nmonth  = ( ndastp - (nyear * 10000) ) / 100
          nday    =   ndastp - (nyear * 10000) - ( nmonth * 100 ) 
 
-         ! Calculates nday_year, day since january 1st (useful to read  daily forcing fields)
+         ! length of the month of the current year (from nleapy, read in namelist)
+         nmonth_len(0) = nbiss(12)   ;   nmonth_len(13) = nbiss(1)
+         SELECT CASE( nleapy )
+         CASE( 1  )   
+            IF( MOD( nyear, 4 ) == 0 ) THEN
+               ;          nmonth_len(1:12) = nbiss(:)      ! 366 days per year (leap year)
+            ELSE
+               ;          nmonth_len(1:12) = nobis(:)      ! 365 days per year
+            ENDIF
+         CASE( 0  )   ;   nmonth_len(1:12) = nobis(:)      ! 365 days per year
+         CASE( 2: )   ;   nmonth_len(1:13) = nleapy        ! 12*nleapy days per year
+         END SELECT
+
+         ! half month in second since the bigining of the year
+         rmonth_half(0) = - 0.5 * rday * REAL( nmonth_len( 0 ) )
+         DO jm = 1, 12
+            rmonth_half(jm) = rmonth_half(jm-1) + 0.5 * rday * REAL( nmonth_len(jm-1) + nmonth_len(jm) )
+         END DO
+         rmonth_half(13) = rmonth_half( 1 ) + 365. * rday
+         rmonth_half(14) = rmonth_half( 2 ) + 365. * rday
+
+         ! day since january 1st (useful to read  daily forcing fields)
          nday_year =  nday
-         !                               ! accumulates days of previous months of this year
-         DO js = 1, nmonth-1
-            IF( nleapy == 1 .AND. MOD( nyear, 4 ) == 0 ) THEN
-               nday_year = nday_year + nbiss(js)
-            ELSE
-               nday_year = nday_year + nobis(js)
-            ENDIF
+         DO js = 1, nmonth - 1             ! accumulates days of previous months of this year
+            nday_year = nday_year + nmonth_len(js)
          END DO
 
-      ENDIF
+         ! number of seconds since...
+         rsec_year  = REAL( nday_year - 1 ) * rday - rdttra(1)      ! 00h 1st day of the current year
+         rsec_day   = REAL( nday      - 1 ) * rday - rdttra(1)      ! 00h 1st day of the current month
+         rsec_month =                              - rdttra(1)      ! 00h         of the current day
 
-      ! I.  calculates adatrj, zadatrjn, zadatrjb.
-      ! ------------------------------------------------------------------
+         ! control print
+         IF(lwp) WRITE(numout,*)' ==============>> time-step =', kt, ' Initial DATE= ',   &
+               &                   nyear, '/', nmonth, '/', nday, '  rsec_day:', rsec_day
 
-      adatrj    = adatrj0 + ( kt - nit000 + 1 ) * rdttra(1) / rday
-      zadatrjn  = adatrj0 + ( kt - nit000     ) * rdttra(1) / rday
-      zadatrjb  = adatrj0 + ( kt - nit000 - 1 ) * rdttra(1) / rday
+         !
+         !                     ! -------------------------------- ! 
+      ELSE                     !  Model calendar at time-step kt  !
+         !                     ! -------------------------------- ! 
 
+         rsec_year  = rsec_year  + rdttra(1)                 ! New time-step
+         rsec_month = rsec_month + rdttra(1)                 ! New time-step
+         rsec_day   = rsec_day   + rdttra(1)                 ! New time-step
 
-      ! II.  increment the date.  The date corresponds to 'now' variables (kt-1),
-      !      which is the time step of forcing fields. 
-      !      Do not do this at nit000  unless nrstdt= 2
-      !      In that case ndastp (read in restart) was for step nit000-2
-      ! -------------------------------------------------------------------
+         adatrj    = adatrj0 + ( kt - nit000 + 1 ) * rdttra(1) / rday
 
-      iday0 = INT( zadatrjb )
-      iday1 = INT( zadatrjn )
+         IF( rsec_day >= rday ) THEN
+            !
+            rsec_day  = 0.e0                               ! NEW day
+            nday      = nday + 1
+            nday_year = nday_year + 1
+            !
+            IF( nday == nmonth_len(nmonth) + 1 ) THEN      ! NEW month
+               nday  = 1
+               rsec_month = 0.e0   
+               nmonth = nmonth + 1
+               IF( nmonth == 13 ) THEN                     ! NEW year
+                  nyear     = nyear + 1
+                  nmonth    = 1
+                  nday_year = 1
+                  rsec_year = 0.e0
+                  !                                        ! update the length of the month
+                  IF( nleapy == 1 ) THEN                   ! of the current year (if necessary)
+                     IF( MOD( nyear, 4 ) == 0 ) THEN
+                        nmonth_len(1:12) = nbiss(:)              ! 366 days per year (leap year)
+                     ELSE
+                        nmonth_len(1:12) = nobis(:)              ! 365 days per year
+                     ENDIF
+                     ! half month in second since the bigining of the year
+                     rmonth_half(0) = - 0.5 * rday * REAL( nmonth_len( 0 ) )
+                     DO jm = 1, 12
+                        rmonth_half(jm) = rmonth_half(jm-1) + 0.5 * rday * REAL( nmonth_len(jm-1) + nmonth_len(jm) )
+                     END DO
+                     rmonth_half(13) = rmonth_half( 1 ) + 365. * rday
+                     rmonth_half(14) = rmonth_half( 2 ) + 365. * rday
+                  ENDIF
+               ENDIF
+            ENDIF
 
-      IF( iday1 - iday0 >= 1 .AND. ( kt /= nit000 .OR. nrstdt == 2 ) ) THEN
+            !
+            ndastp = nyear * 10000 + nmonth * 100 + nday   ! NEW date
+            !
+!           IF(lwp) WRITE(numout,'(a,i8,a,i4,a,i2,a,i2,a,i3)') '======>> time-step =', kt,   &
+!              &   '      New day, DATE= ', nyear, '/', nmonth, '/', nday, '      nday_year = ', nday_year
+!           IF(lwp) WRITE(numout,'(a,F9.0,a,F9.0,a,F9.0)') '         rsec_year = ', rsec_year,   &
+!              &   '   rsec_month = ', rsec_month, '   rsec_day = ', rsec_day
+         ENDIF
 
-         ! increase calendar
-         nyear  =   ndastp / 10000
-         nmonth = ( ndastp - (nyear * 10000) ) / 100
-         nday   =   ndastp - (nyear * 10000) - ( nmonth * 100 ) 
-         nday = nday + 1
-         IF( nleapy == 1 .AND. MOD( nyear, 4 ) == 0 ) THEN
-            iend = nbiss(nmonth)
-         ELSEIF( nleapy > 1 ) THEN 
-            iend = nleapy
-         ELSE 
-            iend = nobis(nmonth)
+         IF(ln_ctl) THEN
+            WRITE(charout,FMT="('kt =', I4,'  d/m/y =',I2,I2,I4)") kt, nday, nmonth, nyear
+            CALL prt_ctl_info(charout)
          ENDIF
-         IF( nday == iend + 1 ) THEN
-            nday  = 1
-            nmonth = nmonth + 1
-            IF( nmonth == 13 ) THEN
-               nmonth  = 1
-               nyear = nyear + 1
-            ENDIF
-         ENDIF
-         ndastp = nyear * 10000 + nmonth * 100 + nday
-
-         ! Calculates nday_year, day since january 1st (useful to read  daily forcing fields)
-         nday_year =  nday
-         !                                ! accumulates days of previous months of this year
-         DO js = 1, nmonth-1
-            IF( nleapy == 1 .AND. MOD( nyear, 4 ) == 0 ) THEN
-               nday_year = nday_year + nbiss(js)
-            ELSE
-               nday_year = nday_year + nobis(js)
-            ENDIF
-         END DO
-
-         IF(lwp) WRITE(numout,*)' ==============>> time-step =', kt, ' New day, DATE= ',   &
-            &                   nyear, '/', nmonth, '/', nday, 'nday_year:', nday_year
-      ENDIF
-
-      IF(ln_ctl) THEN
-         WRITE(charout,FMT="('kt =', I4,'  d/m/y =',I2,I2,I4)") kt, nday, nmonth, nyear
-         CALL prt_ctl_info(charout, itime=kt)
+         !
       ENDIF
 
Index: /trunk/NEMO/OPA_SRC/eosbn2.F90
===================================================================
--- /trunk/NEMO/OPA_SRC/eosbn2.F90	(revision 702)
+++ /trunk/NEMO/OPA_SRC/eosbn2.F90	(revision 703)
@@ -5,4 +5,17 @@
    !!                                               - Brunt-Vaisala frequency 
    !!==============================================================================
+   !! History :       !  89-03  (O. Marti)  Original code
+   !!            6.0  !  94-07  (G. Madec, M. Imbard)  add bn2
+   !!            6.0  !  94-08  (G. Madec)  Add Jackett & McDougall eos
+   !!            7.0  !  96-01  (G. Madec)  statement function for e3
+   !!            8.1  !  97-07  (G. Madec)  introduction of neos, OPA8.1
+   !!            8.1  !  97-07  (G. Madec)  density instead of volumic mass
+   !!                 !  99-02  (G. Madec, N. Grima) semi-implicit pressure gradient
+   !!                 !  01-09  (M. Ben Jelloul)  bugfix onlinear eos
+   !!            8.5  !  02-10  (G. Madec)  add eos_init
+   !!            8.5  !  02-11  (G. Madec, A. Bozec)  partial step, eos_insitu_2d
+   !!            9.0  !  03-08  (G. Madec)  F90, free form
+   !!            9.0  !  06-08  (G. Madec)  add tfreez function
+   !!----------------------------------------------------------------------
 
    !!----------------------------------------------------------------------
@@ -13,7 +26,7 @@
    !!   eos_insitu_2d  : Compute the in situ density for 2d fields
    !!   eos_bn2        : Compute the Brunt-Vaisala frequency
+   !!   tfreez         : Compute the surface freezing temperature
    !!   eos_init       : set eos parameters (namelist)
    !!----------------------------------------------------------------------
-   !! * Modules used
    USE dom_oce         ! ocean space and time domain
    USE phycst          ! physical constants
@@ -33,30 +46,28 @@
    END INTERFACE 
 
-   !! * Routine accessibility
-   PUBLIC eos        ! called by step.F90, inidtr.F90, tranpc.F90 and intgrd.F90
-   PUBLIC bn2        ! called by step.F90
-   PUBLIC eos_init   ! called by step.F90
-
-   !! * Share module variables
-   INTEGER , PUBLIC ::   &  !: nameos : ocean physical parameters
-      neos      = 0,     &  !: = 0/1/2 type of eq. of state and Brunt-Vaisala frequ.
-      neos_init = 0         !: control flag for initialization
-
-   REAL(wp), PUBLIC ::   &  !: nameos : ocean physical parameters
-      ralpha = 2.0e-4,   &  !: thermal expension coeff. (linear equation of state)
-      rbeta  = 7.7e-4       !: saline  expension coeff. (linear equation of state)
+   PUBLIC   eos        ! called by step, istate, tranpc and zpsgrd modules
+   PUBLIC   bn2        ! called by step module
+   PUBLIC   tfreez     ! called by sbcice_... modules
+
+   !!* Namelist (nameos)
+   INTEGER , PUBLIC ::   neos   = 0        !: = 0/1/2 type of eq. of state and Brunt-Vaisala frequ.
+   REAL(wp), PUBLIC ::   ralpha = 2.0e-4   !: thermal expension coeff. (linear equation of state)
+   REAL(wp), PUBLIC ::   rbeta  = 7.7e-4   !: saline  expension coeff. (linear equation of state)
+   NAMELIST/nameos/ neos, ralpha, rbeta
    
+   INTEGER ::   neos_init = 0         !: control flag for initialization
+
    !! * Substitutions
 #  include "domzgr_substitute.h90"
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
-   !!   OPA 9.0 , LOCEAN-IPSL (2005) 
+   !!   OPA 9.0 , LOCEAN-IPSL (2006) 
    !! $Id$
-   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 
+   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
 
 CONTAINS
 
-   SUBROUTINE eos_insitu ( ptem, psal, prd )
+   SUBROUTINE eos_insitu( ptem, psal, prd )
       !!----------------------------------------------------------------------
       !!                   ***  ROUTINE eos_insitu  ***
@@ -92,24 +103,10 @@
       !! ** Action  :   compute prd , the in situ density (no units)
       !!
-      !! References :
-      !!      Jackett, D.R., and T.J. McDougall. J. Atmos. Ocean. Tech., 1994
-      !!
-      !! History :
-      !!        !  89-03 (o. Marti)  Original code
-      !!        ! 94-08 (G. Madec)
-      !!        !  96-01 (G. Madec) statement function for e3
-      !!        !  97-07 (G. Madec) introduction of neos, OPA8.1
-      !!        !  97-07 (G. Madec) density instead of volumic mass
-      !!        !  99-02 (G. Madec, N. Grima) semi-implicit pressure gradient
-      !!        !  01-09 (M. Ben Jelloul) bugfix   
-      !!----------------------------------------------------------------------
-      !! * Arguments
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( in ) ::   &
-         ptem,                 &  ! potential temperature
-         psal                     ! salinity
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( out ) ::   &
-         prd                      ! potential density (surface referenced)
-
-      !! * Local declarations
+      !! References :   Jackett and McDougall, J. Atmos. Ocean. Tech., 1994
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   ptem   ! potential temperature  [Celcius]
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   psal   ! salinity               [psu]
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(  out) ::   prd    ! in situ density 
+      !!
       INTEGER ::  ji, jj, jk      ! dummy loop indices
       REAL(wp) ::   &          
@@ -119,20 +116,15 @@
          zd , zc , zaw, za ,   &  !    "         "
          zb1, za1, zkw, zk0       !    "         "
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   &
-         zws                      ! temporary workspace
-      !!----------------------------------------------------------------------
-
-
-      ! initialization (in not already done)
-      IF( neos_init == 0 ) CALL eos_init
-
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zws   ! temporary workspace
+      !!----------------------------------------------------------------------
+
+      IF( neos_init == 0 ) CALL eos_init      ! initialization (in not already done)
 
       SELECT CASE ( neos )
-
+      !
       CASE ( 0 )               ! Jackett and McDougall (1994) formulation
-
+         !
 !CDIR NOVERRCHK
          zws(:,:,:) = SQRT( ABS( psal(:,:,:) ) )
-
          !                                                ! ===============
          DO jk = 1, jpkm1                                 ! Horizontal slab
@@ -181,8 +173,7 @@
          END DO                                           !   End of slab
          !                                                ! ===============
-
-
+         !
       CASE ( 1 )               ! Linear formulation function of temperature only
-
+         !
          !                                                ! ===============
          DO jk = 1, jpkm1                                 ! Horizontal slab
@@ -199,8 +190,7 @@
          END DO                                           !   End of slab
          !                                                ! ===============
-
-
+         !
       CASE ( 2 )               ! Linear formulation function of temperature and salinity
-
+         !
          !                                                ! ===============
          DO jk = 1, jpkm1                                 ! Horizontal slab
@@ -217,20 +207,18 @@
          END DO                                           !   End of slab
          !                                                ! ===============
-
+         !
       CASE DEFAULT
-
+         !
          WRITE(ctmp1,*) '          bad flag value for neos = ', neos
          CALL ctl_stop( ctmp1 )
-
+         !
       END SELECT
-
-      IF(ln_ctl)   THEN
-         CALL prt_ctl(tab3d_1=prd, clinfo1=' eos  : ', ovlap=1, kdim=jpk)
-      ENDIF
-
+      !
+      IF(ln_ctl)   CALL prt_ctl(tab3d_1=prd, clinfo1=' eos  : ', ovlap=1, kdim=jpk)
+      !
    END SUBROUTINE eos_insitu
 
 
-   SUBROUTINE eos_insitu_pot ( ptem, psal, prd, prhop)
+   SUBROUTINE eos_insitu_pot( ptem, psal, prd, prhop )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE eos_insitu_pot  ***
@@ -275,45 +263,27 @@
       !!              - prhop, the potential volumic mass (Kg/m3)
       !!
-      !! References :
-      !!      Jackett, D.R., and T.J. McDougall. J. Atmos. Ocean. Tech., 1994
-      !!      Brown, J. A. and K. A. Campana. Mon. Weather Rev., 1978
-      !!
-      !! History :
-      !!   4.0  !  89-03  (O. Marti)
-      !!        !  94-08  (G. Madec)
-      !!        !  96-01  (G. Madec) statement function for e3
-      !!        !  97-07  (G. Madec) introduction of neos, OPA8.1
-      !!        !  97-07  (G. Madec) density instead of volumic mass
-      !!        !  99-02  (G. Madec, N. Grima) semi-implicit pressure gradient
-      !!        !  01-09  (M. Ben Jelloul) bugfix   
-      !!   9.0  !  03-08  (G. Madec)  F90, free form
-      !!----------------------------------------------------------------------
-      !! * Arguments
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( in ) ::   &
-         ptem,   &  ! potential temperature
-         psal       ! salinity
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( out ) ::   &
-         prd,    &  ! potential density (surface referenced)
-         prhop      ! potential density (surface referenced)
-
-      !! * Local declarations
+      !! References :   Jackett and McDougall, J. Atmos. Ocean. Tech., 1994
+      !!                Brown and Campana, Mon. Weather Rev., 1978
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   ptem   ! potential temperature  [Celcius]
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   psal   ! salinity               [psu]
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(  out) ::   prd    ! in situ density 
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(  out) ::   prhop  ! potential density (surface referenced)
+
       INTEGER  ::  ji, jj, jk                ! dummy loop indices
       REAL(wp) ::   &             ! temporary scalars
          zt, zs, zh, zsr, zr1, zr2, zr3, zr4, zrhop, ze, zbw,   &
          zb, zd, zc, zaw, za, zb1, za1, zkw, zk0
-      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zws
-      !!----------------------------------------------------------------------
-
-      ! initialization (in not already done)
-      IF( neos_init == 0 ) CALL eos_init
-
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zws
+      !!----------------------------------------------------------------------
+
+      IF( neos_init == 0 ) CALL eos_init      ! initialization (in not already done)
 
       SELECT CASE ( neos )
-
+      !
       CASE ( 0 )               ! Jackett and McDougall (1994) formulation
-
+         !
 !CDIR NOVERRCHK
          zws(:,:,:) = SQRT( ABS( psal(:,:,:) ) )
-
          !                                                ! ===============
          DO jk = 1, jpkm1                                 ! Horizontal slab
@@ -326,5 +296,4 @@
                   zh = fsdept(ji,jj,jk)
                   ! square root salinity
-!!Edmee           zsr= SQRT( ABS( zs ) )
                   zsr= zws(ji,jj,jk)
                   ! compute volumic mass pure water at atm pressure
@@ -366,8 +335,7 @@
          END DO                                           !   End of slab
          !                                                ! ===============
-
-
+         !
       CASE ( 1 )               ! Linear formulation function of temperature only
-
+         !
          !                                                ! ===============
          DO jk = 1, jpkm1                                 ! Horizontal slab
@@ -385,8 +353,7 @@
          END DO                                           !   End of slab
          !                                                ! ===============
-
-
+         !
       CASE ( 2 )               ! Linear formulation function of temperature and salinity
-
+         !
          !                                                ! ===============
          DO jk = 1, jpkm1                                 ! Horizontal slab
@@ -404,19 +371,18 @@
          END DO                                           !   End of slab
          !                                                ! ===============
-
+         !
       CASE DEFAULT
-
+         !
          WRITE(ctmp1,*) '          bad flag value for neos = ', neos
          CALL ctl_stop( ctmp1 )
-
+         !
       END SELECT
-
-      IF(ln_ctl)   THEN
-         CALL prt_ctl(tab3d_1=prd, clinfo1=' eos-p: ', tab3d_2=prhop, clinfo2=' pot : ', ovlap=1, kdim=jpk)
-      ENDIF
-
+      !
+      IF(ln_ctl)   CALL prt_ctl( tab3d_1=prd, clinfo1=' eos-p: ', tab3d_2=prhop, clinfo2=' pot : ', ovlap=1, kdim=jpk )
+      !
    END SUBROUTINE eos_insitu_pot
 
-   SUBROUTINE eos_insitu_2d ( ptem, psal, pdep, prd )
+
+   SUBROUTINE eos_insitu_2d( ptem, psal, pdep, prd )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE eos_insitu_2d  ***
@@ -452,37 +418,27 @@
       !! ** Action  : - prd , the in situ density (no units)
       !!
-      !! References :
-      !!      Jackett, D.R., and T.J. McDougall. J. Atmos. Ocean. Tech., 1994
-      !!
-      !! History :
-      !!   8.5  !  02-11  (G. Madec, A. Bozec)  partial step
-      !!----------------------------------------------------------------------
-      !! * Arguments
-      REAL(wp), DIMENSION(jpi,jpj), INTENT( in ) ::   &
-         ptem,                           &  ! potential temperature
-         psal,                           &  ! salinity
-         pdep                               ! depth
-      REAL(wp), DIMENSION(jpi,jpj), INTENT( out ) ::   &
-         prd                                ! potential density (surface referenced)
-
-      !! * Local declarations
-      INTEGER ::  ji, jj                    ! dummy loop indices
+      !! References :   Jackett and McDougall, J. Atmos. Ocean. Tech., 1994
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   ptem   ! potential temperature  [Celcius]
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   psal   ! salinity               [psu]
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pdep   ! depth                  [m]
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::   prd    ! in situ density 
+      !!
+      INTEGER  ::  ji, jj                    ! dummy loop indices
       REAL(wp) ::   &             ! temporary scalars
          zt, zs, zh, zsr, zr1, zr2, zr3, zr4, zrhop, ze, zbw,   &
          zb, zd, zc, zaw, za, zb1, za1, zkw, zk0,               &
          zmask
-      REAL(wp), DIMENSION(jpi,jpj) :: zws
-      !!----------------------------------------------------------------------
-
-
-      ! initialization (in not already done)
-      IF( neos_init == 0 ) CALL eos_init
+      REAL(wp), DIMENSION(jpi,jpj) ::   zws
+      !!----------------------------------------------------------------------
+
+      IF( neos_init == 0 ) CALL eos_init      ! initialization (in not already done)
 
       prd(:,:) = 0.e0
 
       SELECT CASE ( neos )
-
+      !
       CASE ( 0 )               ! Jackett and McDougall (1994) formulation
-
+      !
 !CDIR NOVERRCHK
          DO jj = 1, jpjm1
@@ -496,5 +452,4 @@
             END DO
          END DO
-
          !                                                ! ===============
          DO jj = 1, jpjm1                                 ! Horizontal slab
@@ -505,5 +460,4 @@
             DO ji = 1, fs_jpim1   ! vector opt.
 #endif
-
                zmask = tmask(ji,jj,1)      ! land/sea bottom mask = surf. mask
 
@@ -543,13 +497,12 @@
                ! masked in situ density anomaly
                prd(ji,jj) = ( zrhop / (  1.0 - zh / ( zk0 - zh * ( za - zh * zb ) )  ) - rau0 )   &
-                          / rau0 * zmask
-            END DO
-            !                                             ! ===============
-         END DO                                           !   End of slab
-         !                                                ! ===============
-
-
+                  &       / rau0 * zmask
+            END DO
+            !                                             ! ===============
+         END DO                                           !   End of slab
+         !                                                ! ===============
+         !
       CASE ( 1 )               ! Linear formulation function of temperature only
-
+         !
          !                                                ! ===============
          DO jj = 1, jpjm1                                 ! Horizontal slab
@@ -565,8 +518,7 @@
          END DO                                           !   End of slab
          !                                                ! ===============
-
-
+         !
       CASE ( 2 )               ! Linear formulation function of temperature and salinity
-
+         !
          !                                                ! ===============
          DO jj = 1, jpjm1                                 ! Horizontal slab
@@ -582,14 +534,14 @@
          END DO                                           !   End of slab
          !                                                ! ===============
-
+         !
       CASE DEFAULT
-
+         !
          WRITE(ctmp1,*) '          bad flag value for neos = ', neos
          CALL ctl_stop( ctmp1 )
-
+         !
       END SELECT
 
-      IF(ln_ctl)   CALL prt_ctl(tab2d_1=prd, clinfo1=' eos2d: ')
-
+      IF(ln_ctl)   CALL prt_ctl( tab2d_1=prd, clinfo1=' eos2d: ' )
+      !
    END SUBROUTINE eos_insitu_2d
 
@@ -623,31 +575,17 @@
       !! ** Action  : - pn2 : the brunt-vaisala frequency
       !!
-      !! References :
-      !!      McDougall, T. J., J. Phys. Oceanogr., 17, 1950-1964, 1987.
-      !!
-      !! History :
-      !!   6.0  !  94-07  (G. Madec, M. Imbard)  Original code
-      !!   8.0  !  97-07  (G. Madec) introduction of statement functions
-      !!   8.5  !  02-07  (G. Madec) Free form, F90
-      !!   8.5  !  02-08  (G. Madec) introduction of arguments
-      !!----------------------------------------------------------------------
-      !! * Arguments
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( in ) ::   &
-         ptem,                           &  ! potential temperature
-         psal                               ! salinity
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( out ) ::   &
-         pn2                               ! Brunt-Vaisala frequency
-
-      !! * Local declarations
+      !! References :   McDougall, J. Phys. Oceanogr., 17, 1950-1964, 1987.
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   ptem   ! potential temperature   [Celcius]
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   psal   ! salinity                [psu]
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(  out) ::   pn2    ! Brunt-Vaisala frequency [s-1]
+
       INTEGER  ::   ji, jj, jk   ! dummy loop indices
-      REAL(wp) ::   &
-         zgde3w, zt, zs, zh,  &  ! temporary scalars 
-         zalbet, zbeta           !    "         "
+      REAL(wp) ::   zgde3w, zt, zs, zh,  &  ! temporary scalars 
+         &          zalbet, zbeta           !    "         "
 #if defined key_zdfddm
       REAL(wp) ::   zds          ! temporary scalars
 #endif
       !!----------------------------------------------------------------------
-      !!  OPA8.5, LODYC-IPSL (2002)
-      !!----------------------------------------------------------------------
 
       ! pn2 : first and last levels
@@ -660,7 +598,7 @@
 
       SELECT CASE ( neos )
-
+      !
       CASE ( 0 )               ! Jackett and McDougall (1994) formulation
-
+         !
          !                                                ! ===============
          DO jk = 2, jpkm1                                 ! Horizontal slab
@@ -712,8 +650,7 @@
          END DO                                           !   End of slab
          !                                                ! ===============
-
-
+         !
       CASE ( 1 )               ! Linear formulation function of temperature only
-
+         !
          !                                                ! ===============
          DO jk = 2, jpkm1                                 ! Horizontal slab
@@ -728,8 +665,7 @@
          END DO                                           !   End of slab
          !                                                ! ===============
-
-
+         !
       CASE ( 2 )               ! Linear formulation function of temperature and salinity
-
+         !
          !                                                ! ===============
          DO jk = 2, jpkm1                                 ! Horizontal slab
@@ -756,22 +692,40 @@
          END DO                                           !   End of slab
          !                                                ! ===============
-
+         !
       CASE DEFAULT
-
+         !
          WRITE(ctmp1,*) '          bad flag value for neos = ', neos
          CALL ctl_stop( ctmp1 )
-
+         !
       END SELECT
 
-      IF(ln_ctl)   THEN
-         CALL prt_ctl(tab3d_1=pn2, clinfo1=' bn2  : ', ovlap=1, kdim=jpk)
+      IF(ln_ctl)   CALL prt_ctl(tab3d_1=pn2, clinfo1=' bn2  : ', ovlap=1, kdim=jpk)
 #if defined key_zdfddm
-         CALL prt_ctl(tab3d_1=rrau, clinfo1=' rrau : ', ovlap=1, kdim=jpk)
+      IF(ln_ctl)   CALL prt_ctl(tab3d_1=rrau, clinfo1=' rrau : ', ovlap=1, kdim=jpk)
 #endif
-      ENDIF
-
+      !
    END SUBROUTINE eos_bn2
 
 
+   FUNCTION tfreez( psal ) RESULT( ptf )
+      !!----------------------------------------------------------------------
+      !!                 ***  ROUTINE eos_init  ***
+      !!
+      !! ** Purpose :   Compute the sea surface freezing temperature [Celcius]
+      !!
+      !! ** Method  :   UNESCO freezing point at the surface (pressure = 0???)
+      !!       freezing point [Celcius]=(-.0575+1.710523e-3*sqrt(abs(s))-2.154996e-4*s)*s-7.53e-4*p
+      !!       checkvalue: tf= -2.588567 Celsius for s=40.0psu, p=500. decibars
+      !!
+      !! Reference  :   UNESCO tech. papers in the marine science no. 28. 1978
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   psal   ! salinity             [psu]
+      REAL(wp), DIMENSION(jpi,jpj)                ::   ptf    ! freezing temperature [Celcius]
+      !!----------------------------------------------------------------------
+      ptf(:,:) = ( - 0.0575 + 1.710523e-3 * SQRT( psal(:,:) )   &
+         &                  - 2.154996e-4 *       psal(:,:)   ) * psal(:,:)
+   END FUNCTION tfreez
+
+
    SUBROUTINE eos_init
       !!----------------------------------------------------------------------
@@ -780,23 +734,10 @@
       !! ** Purpose :   initializations for the equation of state
       !!
-      !! ** Method  :   Read the namelist nameos
-      !!
-      !! ** Action  :   blahblah....
-      !!
-      !! History :
-      !!   8.5  !  02-10  (G. Madec)  Original code
-      !!----------------------------------------------------------------------
-      NAMELIST/nameos/ neos, ralpha, rbeta
-      !!----------------------------------------------------------------------
-      !!  OPA 8.5, LODYC-IPSL (2002)
-      !!----------------------------------------------------------------------
-
-      ! set the initialization flag to 1
-      neos_init = 1           ! indicate that the initialization has been done
-
-      ! namelist nameos : ocean physical parameters
-
-      ! Read Namelist nameos : equation of state
-      REWIND( numnam )
+      !! ** Method  :   Read the namelist nameos and control the parameters
+      !!----------------------------------------------------------------------
+
+      neos_init = 1               ! indicate that the initialization has been done
+
+      REWIND( numnam )            ! Read Namelist nameos : equation of state
       READ  ( numnam, nameos )
 
@@ -807,9 +748,7 @@
          WRITE(numout,*) '~~~~~~~~'
          WRITE(numout,*) '          Namelist nameos : set eos parameters'
-         WRITE(numout,*)
          WRITE(numout,*) '             flag for eq. of state and N^2  neos   = ', neos
          WRITE(numout,*) '             thermal exp. coef. (linear)    ralpha = ', ralpha
          WRITE(numout,*) '             saline  exp. coef. (linear)    rbeta  = ', rbeta
-         WRITE(numout,*)
       ENDIF
 
@@ -817,23 +756,21 @@
 
       CASE ( 0 )               ! Jackett and McDougall (1994) formulation
-
+         IF(lwp) WRITE(numout,*)
          IF(lwp) WRITE(numout,*) '          use of Jackett & McDougall (1994) equation of state and'
          IF(lwp) WRITE(numout,*) '                 McDougall (1987) Brunt-Vaisala frequency'
-
+         !
       CASE ( 1 )               ! Linear formulation function of temperature only
-
+         IF(lwp) WRITE(numout,*)
          IF(lwp) WRITE(numout,*) '          use of linear eos rho(T) = rau0 * ( 1.0285 - ralpha * T )'
          IF( lk_zdfddm ) CALL ctl_stop( '          double diffusive mixing parameterization requires',   &
               &                         ' that T and S are used as state variables' )
-
+         !
       CASE ( 2 )               ! Linear formulation function of temperature and salinity
-
+         IF(lwp) WRITE(numout,*)
          IF(lwp) WRITE(numout,*) '          use of linear eos rho(T,S) = rau0 * ( rbeta * S - ralpha * T )'
-
-      CASE DEFAULT
-
+         !
+      CASE DEFAULT             ! E R R O R in neos 
          WRITE(ctmp1,*) '          bad flag value for neos = ', neos
          CALL ctl_stop( ctmp1 )
-
       END SELECT
 
Index: /trunk/NEMO/OPA_SRC/phycst.F90
===================================================================
--- /trunk/NEMO/OPA_SRC/phycst.F90	(revision 702)
+++ /trunk/NEMO/OPA_SRC/phycst.F90	(revision 703)
@@ -4,5 +4,14 @@
    !!     Definition of of both ocean and ice parameters used in the code
    !!=====================================================================
-   !! * Modules used
+   !! History :        !  90-10  (C. Levy - G. Madec)  Original code
+   !!                  !  91-11  (G. Madec)
+   !!                  !  91-12  (M. Imbard)
+   !!             8.5  !  02-08  (G. Madec, C. Ethe)  F90, add ice constants
+   !!             9.0  !  06-08  (G. Madec)  style 
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   phy_cst  : define and print physical constant and domain parameters
+   !!----------------------------------------------------------------------
    USE par_oce          ! ocean parameters
    USE in_out_manager   ! I/O manager
@@ -11,18 +20,9 @@
    PRIVATE
 
-   !! * Routine accessibility
-   PUBLIC phy_cst          ! routine called by inipar.F90
+   PUBLIC   phy_cst     ! routine called by inipar.F90
 
-   !! * Shared module variables
-   INTEGER, PUBLIC, DIMENSION(12) ::   &  !:
-      nbiss = (/ 31, 29, 31, 30, 31, 30,      &  !: number of days per month
-         &       31, 31, 30, 31, 30, 31 /) ,  &  !  (leap-year)
-      nobis = (/ 31, 28, 31, 30, 31, 30,      &  !: number of days per month
-         &       31, 31, 30, 31, 30, 31 /)       !  (365 days a year)
-   
-   REAL(wp), PUBLIC ::                        &  !:
-      rpi = 3.141592653589793_wp           ,  &  !: pi
-      rad = 3.141592653589793_wp / 180._wp ,  &  !: conversion from degre into radian
-      rsmall = 0.5 * EPSILON( 1. )               !: smallest real computer value
+   REAL(wp), PUBLIC ::   rpi = 3.141592653589793_wp             !: pi
+   REAL(wp), PUBLIC ::   rad = 3.141592653589793_wp / 180._wp   !: conversion from degre into radian
+   REAL(wp), PUBLIC ::   rsmall = 0.5 * EPSILON( 1. )           !: smallest real computer value
    
    REAL(wp), PUBLIC ::          & !:
@@ -61,6 +61,6 @@
       xlic    = 300.33e+6_wp  ,   &  !: volumetric latent heat fusion of ice
       xsn     =   2.8e+6      ,   &  !: latent heat of sublimation of snow
-      rhoic   = 900._wp       ,   &  !: density of sea ice (kg/m3)
-      rhosn   = 330._wp       ,   &  !: density of snow (kg/m3)
+      rhoic   = 900._wp       ,   &  !: volumic mass of sea ice (kg/m3)
+      rhosn   = 330._wp       ,   &  !: volumic mass of snow (kg/m3)
       emic    =   0.97_wp     ,   &  !: emissivity of snow or ice
       sice    =   6.0_wp      ,   &  !: salinity of ice (psu)
@@ -70,9 +70,10 @@
       vkarmn  =   0.4_wp      ,   &  !: von Karman constant
       stefan  =   5.67e-8_wp         !: Stefan-Boltzmann constant 
-      !!----------------------------------------------------------------------
-      !!  OPA 9.0 , LOCEAN-IPSL (2005) 
-      !! $Id$
-      !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
-      !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!  OPA 9.0 , LOCEAN-IPSL (2006) 
+   !! $Id$
+   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
    
 CONTAINS
@@ -83,15 +84,6 @@
       !!
       !! ** Purpose :   Print model parameters and set and print the constants
-      !!
-      !! ** Method  :   no
-      !!
-      !! History :
-      !!        !  90-10  (C. Levy - G. Madec)  Original code
-      !!        !  91-11  (G. Madec)
-      !!        !  91-12  (M. Imbard)
-      !!   8.5  !  02-08  (G. Madec, C. Ethe)  F90, add ice constants 
       !!----------------------------------------------------------------------
-      !! * Local variables
-      CHARACTER (len=64) ::   cform = "(A9, 3(A13, I7) )" 
+      CHARACTER (len=64) ::   cform = "(A12, 3(A13, I7) )" 
       !!----------------------------------------------------------------------
 
@@ -103,22 +95,17 @@
       ! ----------------
       IF(lwp) THEN
-         WRITE(numout,*) '       parameter file'
-         WRITE(numout,*)
+         WRITE(numout,*) '       Domain info'
          WRITE(numout,*) '          dimension of model'
-         WRITE(numout,*) '              Local domain      Global domain       Data domain '
-         WRITE(numout,cform) '         ','   jpi     : ', jpi, '   jpiglo  : ', jpiglo, '   jpidta  : ', jpidta
-         WRITE(numout,cform) '         ','   jpj     : ', jpj, '   jpjglo  : ', jpjglo, '   jpjdta  : ', jpjdta
-         WRITE(numout,cform) '         ','   jpk     : ', jpk, '   jpk     : ', jpk   , '   jpkdta  : ', jpkdta
-         WRITE(numout,*)      '        ','   jpij    : ', jpij
-         WRITE(numout,*)
+         WRITE(numout,*) '                 Local domain      Global domain       Data domain '
+         WRITE(numout,cform) '            ','   jpi     : ', jpi, '   jpiglo  : ', jpiglo, '   jpidta  : ', jpidta
+         WRITE(numout,cform) '            ','   jpj     : ', jpj, '   jpjglo  : ', jpjglo, '   jpjdta  : ', jpjdta
+         WRITE(numout,cform) '            ','   jpk     : ', jpk, '   jpk     : ', jpk   , '   jpkdta  : ', jpkdta
+         WRITE(numout,*)      '           ','   jpij    : ', jpij
          WRITE(numout,*) '          mpp local domain info (mpp)'
          WRITE(numout,*) '             jpni    : ', jpni, '   jpreci  : ', jpreci
          WRITE(numout,*) '             jpnj    : ', jpnj, '   jprecj  : ', jprecj
          WRITE(numout,*) '             jpnij   : ', jpnij
-
-         WRITE(numout,*)
          WRITE(numout,*) '          lateral domain boundary condition type : jperio  = ', jperio
-         WRITE(numout,*) '          domain island (use in rigid-lid case)  : jpisl   = ', jpisl 
-         WRITE(numout,*) '                                                   jpnisl  = ', jpnisl
+         WRITE(numout,*) '          domain island (use in rigid-lid case)  : jpisl   = ', jpisl, '   jpnisl  = ', jpnisl
       ENDIF
 
@@ -126,5 +113,5 @@
       ! ----------------
       IF(lwp) WRITE(numout,*)
-      IF(lwp) WRITE(numout,*) '       constants'
+      IF(lwp) WRITE(numout,*) '       Constants'
 
       IF(lwp) WRITE(numout,*)
Index: /trunk/NEMO/OPA_SRC/step.F90
===================================================================
--- /trunk/NEMO/OPA_SRC/step.F90	(revision 702)
+++ /trunk/NEMO/OPA_SRC/step.F90	(revision 703)
@@ -4,6 +4,5 @@
    !! Time-stepping    : manager of the ocean, tracer and ice time stepping
    !!======================================================================
-   !! History :        !  91-03  ()  Original code
-   !!                  !  91-11  (G. Madec)
+   !! History :        !  91-03  (G. Madec)  Original code
    !!                  !  92-06  (M. Imbard)  add a first output record
    !!                  !  96-04  (G. Madec)  introduction of dynspg
@@ -20,4 +19,7 @@
    !!             " "  !  06-01  (L. Debreu, C. Mazauric)  Agrif implementation
    !!             " "  !  06-07  (S. Masson)  restart using iom
+   !!             " "  !  06-08  (G. Madec)  surface module 
+   !!----------------------------------------------------------------------
+
    !!----------------------------------------------------------------------
    !!   stp            : OPA system time-stepping
@@ -30,5 +32,5 @@
    USE cpl_oce         ! coupled ocean-atmosphere variables
    USE in_out_manager  ! I/O manager
-   USE iom
+   USE iom             !
    USE lbclnk
 
@@ -37,12 +39,6 @@
    USE dtatem          ! ocean temperature data           (dta_tem routine)
    USE dtasal          ! ocean salinity    data           (dta_sal routine)
-   USE dtasst          ! ocean sea surface temperature    (dta_sst routine)
-   USE dtasss          ! ocean sea surface salinity       (dta_sss routine)
-   USE taumod          ! surface stress                   (tau     routine)
-   USE flxmod          ! thermohaline fluxes              (flx     routine)
-   USE ocesbc          ! thermohaline fluxes              (oce_sbc routine)
-   USE flxrnf          ! runoffs                          (flx_rnf routine)
-   USE flxfwb          ! freshwater budget correction     (flx_fwb routine)
-   USE closea          ! closed sea freshwater budget     (flx_clo routine)
+   USE sbcmod          ! surface boundary condition       (sbc     routine)
+   USE sbcrnf          ! surface boundary condition: runoff variables
    USE ocfzpt          ! surface ocean freezing point    (oc_fz_pt routine)
 
@@ -161,12 +157,10 @@
       !!              -8- Outputs and diagnostics
       !!----------------------------------------------------------------------
-      !! * Arguments
 #if defined key_agrif   
-      INTEGER               :: kstp   ! ocean time-step index
+      INTEGER             ::   kstp   ! ocean time-step index
 #else
-      INTEGER, INTENT( in ) :: kstp   ! ocean time-step index
+      INTEGER, INTENT(in) ::   kstp   ! ocean time-step index
 #endif      
-
-      !! * local declarations
+      INTEGER ::   jk       ! dummy loop indice
       INTEGER ::   indic    ! error indicator if < 0
       !! ---------------------------------------------------------------------
@@ -179,5 +173,7 @@
       indic = 1                    ! reset to no error condition
 
+!!gm: attention n'est plus ds le step de gm
       adatrj = adatrj + rdt/86400._wp
+!!gm: attention n'est plus ds le step de gm
 
       CALL day( kstp )             ! Calendar
@@ -186,50 +182,18 @@
 
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-      ! Update data, open boundaries and Forcings
+      ! Update data, open boundaries, surface boundary condition (including sea-ice)
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 
       IF( lk_dtatem  )   CALL dta_tem( kstp )         ! update 3D temperature data
-
-      IF( lk_dtasal  )   CALL dta_sal( kstp )         ! Salinity data
-
-      IF( lk_dtasst  )   CALL dta_sst( kstp )         ! Sea Surface Temperature data
-
-      IF( lk_dtasss  )   CALL dta_sss( kstp )         ! Sea Surface Salinity data
+      IF( lk_dtasal  )   CALL dta_sal( kstp )         ! update 3D salinity data
+
+                         CALL sbc    ( kstp )         ! Sea Boundary Condition (including sea-ice)
 
       IF( lk_obc     )   CALL obc_dta( kstp )         ! update dynamic and tracer data at open boundaries
-
       IF( lk_obc     )   CALL obc_rad( kstp )         ! compute phase velocities at open boundaries
 
-      IF( .NOT. lk_core )    CALL tau( kstp )         ! wind stress
-
-                         CALL flx_rnf( kstp )         ! runoff data
-
-                         CALL flx( kstp )             ! heat and freshwater fluxes
-
-      IF( lk_ice_lim )   CALL ice_stp( kstp )         ! sea-ice model (Update stress & fluxes)
-
-                         CALL oce_sbc( kstp )         ! ocean surface boudaries
-
-      IF( ln_fwb     )   CALL flx_fwb( kstp )         ! freshwater budget
-
-      IF( nclosea == 1 ) CALL flx_clo( kstp )         ! closed sea in the domain (update freshwater fluxes)
-
-      IF( kstp == nit000 ) THEN 
-         IF( ninist == 1 ) THEN                       ! Output the initial state and forcings
-            CALL dia_wri_state( 'output.init' )
-         ENDIF
-      ENDIF
-
-      IF(ln_ctl) THEN         ! print mean trends (used for debugging)
-         CALL prt_ctl(tab2d_1=emp    , clinfo1=' emp  -   : ', mask1=tmask, ovlap=1)
-         CALL prt_ctl(tab2d_1=emps   , clinfo1=' emps -   : ', mask1=tmask, ovlap=1)
-         CALL prt_ctl(tab2d_1=qt     , clinfo1=' qt   -   : ', mask1=tmask, ovlap=1)
-         CALL prt_ctl(tab2d_1=qsr    , clinfo1=' qsr  -   : ', mask1=tmask, ovlap=1)
-         CALL prt_ctl(tab2d_1=runoff , clinfo1=' runoff   : ', mask1=tmask, ovlap=1)
-         CALL prt_ctl(tab3d_1=tmask  , clinfo1=' tmask    : ', mask1=tmask, ovlap=1, kdim=jpk)
-         CALL prt_ctl(tab3d_1=tn     , clinfo1=' sst  -   : ', mask1=tmask, ovlap=1, kdim=1)
-         CALL prt_ctl(tab3d_1=sn     , clinfo1=' sss  -   : ', mask1=tmask, ovlap=1, kdim=1)
-         CALL prt_ctl(tab2d_1=taux   , clinfo1=' tau  - x : ', mask1=umask, &
-            &         tab2d_2=tauy   , clinfo2='      - y : ', mask2=vmask,ovlap=1)
+      IF( ninist == 1 ) THEN                       ! Output the initial state and forcings
+                        CALL dia_wri_state( 'output.init' )
+                        ninist = 0
       ENDIF
 
@@ -244,5 +208,5 @@
       !-----------------------------------------------------------------------
 
-                       CALL bn2( tb, sb, rn2 )              ! before Brunt-Vaisala frequency
+                        CALL bn2( tb, sb, rn2 )              ! before Brunt-Vaisala frequency
       
       !                                                     ! Vertical eddy viscosity and diffusivity coefficients
@@ -261,16 +225,11 @@
       ENDIF
 
-      IF( cp_cfg == "orca" ) THEN                           ! ORCA: Reduce vertical mixing in some specific areas
-         SELECT CASE ( jp_cfg )
-            CASE ( 05 )                         ! ORCA R2 configuration
-               avt  (:,:,2) = avt  (:,:,2) + 1.e-3 * upsrnfh(:,:)   ! increase diffusivity of rivers mouths
-            CASE ( 025 )                         ! ORCA R025 configuration
-               avt  (:,:,2) = avt  (:,:,2) + 2.e-3 * upsrnfh(:,:)   ! increase diffusivity of rivers mouths
-         END SELECT
+      IF( nn_runoff /=0 ) THEN                              ! increase diffusivity at rivers mouths
+         DO jk = 2, nkrnf   ;   avt(:,:,jk) = avt(:,:,jk) + rn_avt_rnf * rnfmsk(:,:)   ;   END DO
       ENDIF
 
       IF( ln_zdfevd )   CALL zdf_evd( kstp )                 ! enhanced vertical eddy diffusivity
 
-      IF( lk_zdfddm .AND. .NOT. lk_zdfkpp)   &
+      IF( lk_zdfddm .AND. .NOT. lk_zdfkpp )   &
          &              CALL zdf_ddm( kstp )                 ! double diffusive mixing
 
@@ -285,12 +244,9 @@
       ! N.B. ua, va, ta, sa arrays are used as workspace in this section
       !-----------------------------------------------------------------------
-
       IF( lk_ldfslp     )   CALL ldf_slp( kstp, rhd, rn2 )       ! before slope of the lateral mixing
-
 #if defined key_traldf_c2d
       IF( lk_traldf_eiv )   CALL ldf_eiv( kstp )                 ! eddy induced velocity coefficient
 #endif
 
-
 #if defined key_passivetrc
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
@@ -299,9 +255,6 @@
       ! N.B. ua, va, ta, sa arrays are used as workspace in this section
       !-----------------------------------------------------------------------
-
                              CALL trc_stp( kstp, indic )            ! time-stepping
-
-#endif
-
+#endif
 
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
@@ -310,25 +263,16 @@
       ! N.B. ua, va arrays are used as workspace in this section
       !-----------------------------------------------------------------------
-
                              ta(:,:,:) = 0.e0               ! set tracer trends to zero
                              sa(:,:,:) = 0.e0
 
                              CALL tra_sbc    ( kstp )       ! surface boundary condition
-
       IF( ln_traqsr      )   CALL tra_qsr    ( kstp )       ! penetrative solar radiation qsr
-
       IF( lk_trabbc      )   CALL tra_bbc    ( kstp )       ! bottom heat flux
-
       IF( lk_trabbl_dif  )   CALL tra_bbl_dif( kstp )       ! diffusive bottom boundary layer scheme
       IF( lk_trabbl_adv  )   CALL tra_bbl_adv( kstp )       ! advective (and/or diffusive) bottom boundary layer scheme
-
       IF( lk_tradmp      )   CALL tra_dmp    ( kstp )       ! internal damping trends
-
                              CALL tra_adv    ( kstp )       ! horizontal & vertical advection
-
       IF( n_cla == 1     )   CALL tra_cla    ( kstp )       ! Cross Land Advection (Update Hor. advection)
-
       IF( lk_zdfkpp )        CALL tra_kpp    ( kstp )       ! KPP non-local tracer fluxes
-
                              CALL tra_ldf    ( kstp )       ! lateral mixing
 #if defined key_agrif
@@ -336,7 +280,5 @@
 #endif
                              CALL tra_zdf    ( kstp )       ! vertical mixing
-
                              CALL tra_nxt( kstp )           ! tracer fields at next time step
-
       IF( ln_zdfnpc      )   CALL tra_npc( kstp )           ! update the new (t,s) fields by non
       !                                                     ! penetrative convective adjustment
@@ -359,13 +301,9 @@
       ! N.B. ta, sa arrays are used as workspace in this section 
       !-----------------------------------------------------------------------
-
-
                                ua(:,:,:) = 0.e0               ! set dynamics trends to zero
                                va(:,:,:) = 0.e0
 
                                CALL dyn_adv( kstp )           ! advection (vector or flux form)
-
                                CALL dyn_vor( kstp )           ! vorticity term including Coriolis
-
                                CALL dyn_ldf( kstp )           ! lateral mixing
 #if defined key_agrif
@@ -373,18 +311,11 @@
 #endif
                                CALL dyn_hpg( kstp )           ! horizontal gradient of Hydrostatic pressure
-
                                CALL dyn_zdf( kstp )           ! vertical diffusion
-
       IF( lk_dynspg_rl ) THEN
          IF( lk_obc    )       CALL obc_spg( kstp )           ! surface pressure gradient at open boundaries
       ENDIF
-                       indic=0
-!i bug lbc sur emp
-      CALL lbc_lnk( emp, 'T', 1. )
-!i
+                               indic=0
                                CALL dyn_spg( kstp, indic )    ! surface pressure gradient
-
                                CALL dyn_nxt( kstp )           ! lateral velocity at next time step
-
       IF( lk_vvl )             CALL dom_vvl                   ! vertical mesh at next time step
 
@@ -395,24 +326,14 @@
       ! N.B. ua, va, ta, sa arrays are used as workspace in this section
       !-----------------------------------------------------------------------
-
                        CALL oc_fz_pt                        ! ocean surface freezing temperature
-
                        CALL div_cur( kstp )                 ! Horizontal divergence & Relative vorticity
-
       IF( n_cla == 1 ) CALL div_cla( kstp )                 ! Cross Land Advection (Update Hor. divergence)
-
                        CALL wzv( kstp )                     ! Vertical velocity
 
-
-
-      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-      ! Control, and restarts
-      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-      ! N.B. ua, va, ta, sa arrays are used as workspace in this section
-      !-----------------------------------------------------------------------
-      !                                                     ! Time loop: control and print
+      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+      ! Control and restarts
+      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
                                  CALL stp_ctl( kstp, indic )
       IF( indic < 0          )   CALL ctl_stop( 'step: indic < 0' )
-
       IF( kstp == nit000     )   CALL iom_close( numror )             ! close input  ocean restart file
       IF( lrst_oce           )   CALL rst_write  ( kstp )             ! write output ocean restart file
@@ -425,5 +346,5 @@
       !-----------------------------------------------------------------------
 
-      IF ( nstop == 0 ) THEN                                ! Diagnostics
+      IF( nstop == 0 ) THEN                                 ! Diagnostics:
          IF( lk_floats  )   CALL flo_stp( kstp )                 ! drifting Floats
          IF( lk_trddyn  )   CALL trd_dwr( kstp )                 ! trends: dynamics 
@@ -437,7 +358,6 @@
          IF( lk_diafwb  )   CALL dia_fwb( kstp )                 ! Fresh water budget diagnostics
          IF( ln_diaptr  )   CALL dia_ptr( kstp )                 ! Poleward TRansports diagnostics
-
-         !                                                 ! Outputs
-                            CALL dia_wri    ( kstp, indic )      ! ocean model: outputs
+         !                                                 ! outputs
+                            CALL dia_wri( kstp, indic )          ! ocean model: outputs
       ENDIF
 
