Index: /trunk/NEMO/LIM_SRC/ice.F90
===================================================================
--- /trunk/NEMO/LIM_SRC/ice.F90	(revision 87)
+++ /trunk/NEMO/LIM_SRC/ice.F90	(revision 88)
@@ -31,10 +31,10 @@
       resl   = 5.0e-05,  &  !: maximum value for the residual of relaxation
       cw     = 5.0e-03,  &  !: drag coefficient for oceanic stress
-      angvg  = 0.0    ,  &  !: turning angle for oceanic stress
+      angvg  = 0.e0   ,  &  !: turning angle for oceanic stress
       pstar  = 1.0e+04,  &  !: first bulk-rheology parameter
-      c_rhg  = 20.0   ,  &  !: second bulk-rhelogy parameter
-      etamn  = 0.0e+07,  &  !: minimun value for viscosity
-      creepl = 2.0e-08,  &  !: creep limit
-      ecc    = 2.0    ,  &  !: eccentricity of the elliptical yield curve
+      c_rhg  = 20.e0  ,  &  !: second bulk-rhelogy parameter
+      etamn  = 0.e+07,   &  !: minimun value for viscosity
+      creepl = 2.e-08,   &  !: creep limit
+      ecc    = 2.e0   ,  &  !: eccentricity of the elliptical yield curve
       ahi0   = 350.e0       !: sea-ice hor. eddy diffusivity coeff. (m2/s)
 
@@ -98,11 +98,11 @@
       tio_u, tio_v       !: two components of the ice-ocean stress (N/m2)
 
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,nsmax) ::   &  !:
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpsmax) ::   &  !:
       scal0              !: ???
 
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,nlayersp1) ::   &  !:
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jplayersp1) ::   &  !:
       tbif          !: Temperature inside the ice/snow layer
 
-   REAL(wp), DIMENSION(jpi,jpj,0:kmax+1) ::    &  !:
+   REAL(wp), DIMENSION(jpi,jpj,0:jpkmax+1) ::    &  !:
       reslum        !: Relative absorption of solar radiation in each ocean level
 
Index: /trunk/NEMO/LIM_SRC/icestp.F90
===================================================================
--- /trunk/NEMO/LIM_SRC/icestp.F90	(revision 87)
+++ /trunk/NEMO/LIM_SRC/icestp.F90	(revision 88)
@@ -16,4 +16,5 @@
    USE flx_oce         ! forcings variables
    USE dom_ice
+   USE cpl_oce
    USE daymod
    USE phycst          ! Define parameters for the routines
@@ -76,13 +77,13 @@
 
       IF( kt == nit000 ) THEN
-# if defined key_coupled
-         IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) 'ice_stp : Louvain la Neuve Ice Model (LIM)'
-         IF(lwp) WRITE(numout,*) '~~~~~~~   coupled case'
-#else
-         IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) 'ice_stp : Louvain la Neuve Ice Model (LIM)' 
-         IF(lwp) WRITE(numout,*) '~~~~~~~   forced case using bulk formulea'
-#endif
+         IF( lk_cpl ) THEN
+            IF(lwp) WRITE(numout,*)
+            IF(lwp) WRITE(numout,*) 'ice_stp : Louvain la Neuve Ice Model (LIM)'
+            IF(lwp) WRITE(numout,*) '~~~~~~~   coupled case'
+         ELSE
+            IF(lwp) WRITE(numout,*)
+            IF(lwp) WRITE(numout,*) 'ice_stp : Louvain la Neuve Ice Model (LIM)' 
+            IF(lwp) WRITE(numout,*) '~~~~~~~   forced case using bulk formulea'
+         ENDIF
          !  Initialize fluxes fields
          gtaux(:,:) = 0.e0
@@ -183,5 +184,5 @@
          
 
-         IF( l_ctl .AND. lwp ) THEN         ! print mean trends (used for debugging)
+         IF( l_ctl ) THEN         ! print mean trends (used for debugging)
             WRITE(numout,*) 'Ice Forcings '
             WRITE(numout,*) ' qsr_oce  : ', SUM( qsr_oce (:,:) ), ' qsr_ice  : ', SUM( qsr_ice (:,:) )
@@ -202,5 +203,5 @@
          CALL lim_dyn                                                ! Ice dynamics !   ( rheology/dynamics )
          !                                                           !--------------!
-         IF( l_ctl .AND. lwp ) THEN
+         IF( l_ctl ) THEN
             WRITE(numout,*) ' hsnif  2 : ', SUM( hsnif  (:,:) ), ' hicnif   : ', SUM( hicif  (:,:) )
             WRITE(numout,*) ' frld   2 : ', SUM( frld   (:,:) ), ' sist     : ', SUM( sist   (:,:) )
@@ -211,5 +212,5 @@
          CALL lim_trp                                                ! Ice transport !  ( Advection/diffusion )
          !                                                           !---------------!
-         IF( l_ctl .AND. lwp ) THEN
+         IF( l_ctl ) THEN
             WRITE(numout,*) ' hsnif  3 : ', SUM( hsnif  (:,:) ), ' hicnif   : ', SUM( hicif  (:,:) )
             WRITE(numout,*) ' frld   3 : ', SUM( frld   (:,:) ), ' sist     : ', SUM( sist   (:,:) )
@@ -220,5 +221,5 @@
          CALL lim_thd                                                ! Ice thermodynamics !
          !                                                           !--------------------!
-         IF( l_ctl .AND. lwp ) THEN
+         IF( l_ctl ) THEN
             WRITE(numout,*) ' hsnif  4 : ', SUM( hsnif  (:,:) ), ' hicnif   : ', SUM( hicif  (:,:) )
             WRITE(numout,*) ' frld   4 : ', SUM( frld   (:,:) ), ' sist     : ', SUM( sist   (:,:) )
@@ -270,21 +271,10 @@
 #else
    !!----------------------------------------------------------------------
-   !!   Default option                                 NO LIM sea-ice model
-   !!----------------------------------------------------------------------
-   USE in_out_manager
-
+   !!   Default option           Dummy module          NO LIM sea-ice model
+   !!----------------------------------------------------------------------
 CONTAINS
-
-   SUBROUTINE ice_stp ( kt )            ! Empty routine
-      INTEGER, INTENT( in ) ::   kt     ! ocean time-step index
-
-      IF( kt == nit000 ) THEN
-         IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) 'No Sea Ice Model'
-         IF(lwp) WRITE(numout,*) '~~~~~~~'
-      ENDIF
-
+   SUBROUTINE ice_stp ( kt )     ! Dummy routine
+      WRITE(*,*) 'ice_stp: You should not have seen this print! error?', kt
    END SUBROUTINE ice_stp
-
 #endif
 
Index: /trunk/NEMO/LIM_SRC/limadv.F90
===================================================================
--- /trunk/NEMO/LIM_SRC/limadv.F90	(revision 87)
+++ /trunk/NEMO/LIM_SRC/limadv.F90	(revision 88)
@@ -1,9 +1,10 @@
 MODULE limadv 
-#if defined key_ice_lim
    !!======================================================================
    !!                       ***  MODULE limadv   ***
    !! LIM sea-ice model : sea-ice advection
    !!======================================================================
-
+#if defined key_ice_lim
+   !!----------------------------------------------------------------------
+   !!   'key_ice_lim'                                     LIM sea-ice model
    !!----------------------------------------------------------------------
    !!   lim_adv_x  : advection of sea ice on x axis
@@ -221,5 +222,5 @@
       CALL lbc_lnk( psxy, 'T', 1. )
 
-      IF( l_ctl .AND. lwp ) THEN
+      IF(l_ctl) THEN
          WRITE(numout,*) ' lim_adv_x: psm  ', SUM( psm  ), ' ps0  ', SUM( ps0  )
          WRITE(numout,*) ' lim_adv_x: psx  ', SUM( psx  ), ' psxx ', SUM( psxx )
@@ -421,5 +422,5 @@
       CALL lbc_lnk( psxy, 'T', 1. )
 
-      IF( l_ctl .AND. lwp ) THEN
+      IF(l_ctl) THEN
          WRITE(numout,*) ' lim_adv_y: psm  ', SUM( psm  ), ' ps0  ', SUM( ps0  )
          WRITE(numout,*) ' lim_adv_y: psx  ', SUM( psx  ), ' psxx ', SUM( psxx )
@@ -430,10 +431,8 @@
    END SUBROUTINE lim_adv_y
 
-   !!======================================================================
 #else
-   !!==============================================================================
-   !!                       ***  MODULE limadv   ***
-   !!                              No sea ice
-   !!==============================================================================
+   !!----------------------------------------------------------------------
+   !!   Default option            Dummy module         NO LIM sea-ice model
+   !!----------------------------------------------------------------------
 CONTAINS
    SUBROUTINE lim_adv_x         ! Empty routine
Index: /trunk/NEMO/LIM_SRC/limflx.F90
===================================================================
--- /trunk/NEMO/LIM_SRC/limflx.F90	(revision 87)
+++ /trunk/NEMO/LIM_SRC/limflx.F90	(revision 88)
@@ -9,4 +9,5 @@
    !!----------------------------------------------------------------------
    !!   lim_flx  : flux at the ice / ocean interface
+   !!----------------------------------------------------------------------
    !! * Modules used
    USE par_oce
@@ -26,9 +27,10 @@
    PUBLIC lim_flx       ! called by lim_step
 
-  !! * Module variables
-     REAL(wp)  ::            &  ! constant values
-         epsi16 = 1e-16   ,  &
-         rzero  = 0.0    ,  &
-         rone   = 1.0
+   !! * Module variables
+   REAL(wp)  ::            &  ! constant values
+      epsi16 = 1.e-16  ,  &
+      rzero  = 0.e0    ,  &
+      rone   = 1.e0
+
    !! * Substitutions
 #  include "vectopt_loop_substitute.h90"
@@ -41,5 +43,4 @@
       !!-------------------------------------------------------------------
       !!                ***  ROUTINE lim_flx ***
-      !!                
       !!  
       !! ** Purpose : Computes the mass and heat fluxes to the ocean
@@ -53,5 +54,4 @@
       !!              - fmass   : freshwater flux at sea ice/ocean interface
       !!
-      !!
       !! ** References :
       !!       H. Goosse et al. 1996, Bul. Soc. Roy. Sc. Liege, 65, 87-90
@@ -59,5 +59,4 @@
       !!         addition    : 02-07 (C. Ethe, G. Madec)
       !!---------------------------------------------------------------------
-      !! * Modules used
       !! * Local variables
       INTEGER ::   ji, jj         ! dummy loop indices
@@ -174,12 +173,8 @@
       !-----------------------------------------------!
 
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            ftaux (ji,jj) = - tio_u(ji,jj) * rau0   ! taux ( ice: N/m2/rau0, ocean: N/m2 )
-            ftauy (ji,jj) = - tio_v(ji,jj) * rau0   ! tauy ( ice: N/m2/rau0, ocean: N/m2 )                
-            freeze(ji,jj) = 1.0 - frld(ji,jj)       ! Sea ice cover            
-            tn_ice(ji,jj) = sist(ji,jj)             ! Ice surface temperature                      
-         END DO
-      END DO
+      ftaux (:,:) = - tio_u(:,:) * rau0   ! taux ( ice: N/m2/rau0, ocean: N/m2 )
+      ftauy (:,:) = - tio_v(:,:) * rau0   ! tauy ( ice: N/m2/rau0, ocean: N/m2 )                
+      freeze(:,:) = 1.0 - frld(:,:)       ! Sea ice cover            
+      tn_ice(:,:) = sist(:,:)             ! Ice surface temperature                      
 
 #if defined key_coupled            
@@ -194,12 +189,8 @@
       CALL flx_blk_albedo( zalb, zalcn, zalbp, zaldum )
 
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            alb_ice(ji,jj) =  0.5 * zalbp(ji,jj) + 0.5 * zalb (ji,jj)  ! Ice albedo                       
-         END DO
-      END DO
-#endif
-
-      IF( l_ctl .AND. lwp ) THEN
+      alb_ice(:,:) =  0.5 * zalbp(:,:) + 0.5 * zalb (:,:)   ! Ice albedo                       
+#endif
+
+      IF(l_ctl) THEN
          WRITE(numout,*) ' lim_flx  '
          WRITE(numout,*) ' fsolar ', SUM(fsolar), ' fnsolar', SUM( fnsolar )
@@ -208,5 +199,4 @@
          WRITE(numout,*) ' freeze ', SUM(freeze), ' tn_ice ', SUM(tn_ice)
       ENDIF 
- 
    
     END SUBROUTINE lim_flx
@@ -214,10 +204,11 @@
 #else
    !!----------------------------------------------------------------------
-   !!   Default option :        Empty module           NO LIM sea-ice model
+   !!   Default option :        Dummy module           NO LIM sea-ice model
    !!----------------------------------------------------------------------
 CONTAINS
-   SUBROUTINE lim_flx         ! Empty routine
+   SUBROUTINE lim_flx         ! Dummy routine
    END SUBROUTINE lim_flx
 #endif 
 
+   !!======================================================================
 END MODULE limflx
Index: /trunk/NEMO/LIM_SRC/limhdf.F90
===================================================================
--- /trunk/NEMO/LIM_SRC/limhdf.F90	(revision 87)
+++ /trunk/NEMO/LIM_SRC/limhdf.F90	(revision 88)
@@ -1,9 +1,10 @@
 MODULE limhdf
-#if defined key_ice_lim
    !!======================================================================
    !!                    ***  MODULE limhdf   ***
-   !! LIM diffusion ice model : sea-ice variables horizontal diffusion
+   !! LIM ice model : horizontal diffusion of sea-ice quantities
    !!======================================================================
-
+#if defined key_ice_lim
+   !!----------------------------------------------------------------------
+   !!   'key_ice_lim'                                     LIM sea-ice model
    !!----------------------------------------------------------------------
    !!   lim_hdf  : diffusion trend on sea-ice variable
@@ -96,5 +97,5 @@
             DO ji = fs_2 , fs_jpim1   ! vector opt.
                zfact(ji,jj) = ( e2u(ji,jj) + e2u(ji-1,jj  ) + e1v(ji,jj) + e1v(ji,jj-1) ) &
-                             / ( e1t(ji,jj) * e2t(ji,jj) )
+                  &          / ( e1t(ji,jj) * e2t(ji,jj) )
             END DO
          END DO
@@ -144,5 +145,5 @@
 
          ! convergence test
-         zconv = 0.0
+         zconv = 0.e0
          DO jj = 2, jpjm1
             DO ji = 2, jpim1
@@ -165,14 +166,13 @@
 
       ptab(:,:) = ptab(:,:)
-      IF( l_ctl .AND. lwp ) THEN
-         WRITE(numout,*) ' lim_hdf  : ', SUM( ptab-ptab0 ), ' zconv= ', zconv, ' iter= ', iter
-      ENDIF
+
+      IF(l_ctl)   WRITE(numout,*) ' lim_hdf  : ', SUM( ptab-ptab0 ), ' zconv= ', zconv, ' iter= ', iter
 
    END SUBROUTINE lim_hdf
+
 #else
-   !!======================================================================
-   !!                       ***  MODULE limhdf   ***
-   !!                          no sea ice model
-   !!======================================================================
+   !!----------------------------------------------------------------------
+   !!   Default option          Dummy module           NO LIM sea-ice model
+   !!----------------------------------------------------------------------
 CONTAINS
    SUBROUTINE lim_hdf         ! Empty routine
Index: /trunk/NEMO/LIM_SRC/limrst.F90
===================================================================
--- /trunk/NEMO/LIM_SRC/limrst.F90	(revision 87)
+++ /trunk/NEMO/LIM_SRC/limrst.F90	(revision 88)
@@ -84,7 +84,7 @@
       zinfo(2) = FLOAT( it0   )   ! iteration number
 
-      zsec     = 0.
+      zsec     = 0.e0
       itime    = 0
-      zdept(1) = 0.
+      zdept(1) = 0.e0
       zdt      = rdt_ice * nstock
 
Index: /trunk/NEMO/LIM_SRC/limthd_lac.F90
===================================================================
--- /trunk/NEMO/LIM_SRC/limthd_lac.F90	(revision 87)
+++ /trunk/NEMO/LIM_SRC/limthd_lac.F90	(revision 88)
@@ -22,10 +22,10 @@
    PUBLIC lim_thd_lac     ! called by lim_thd
 
-  !! * Module variables
-     REAL(wp)  ::            &  ! constant values
-         epsi20 = 1e-20   ,  &
-         epsi13 = 1e-13   ,  &
-         zzero  = 0.0     ,  &
-         zone   = 1.0
+   !! * Module variables
+   REAL(wp)  ::            &  ! constant values
+      epsi20 = 1.e-20  ,  &
+      epsi13 = 1.e-13  ,  &
+      zzero  = 0.e0    ,  &
+      zone   = 1.e0
    !!----------------------------------------------------------------------
    !!   LIM 2.0 , UCL-LODYC-IPSL  (2003)
@@ -33,5 +33,5 @@
 CONTAINS
     
-    SUBROUTINE lim_thd_lac( kideb, kiut )
+   SUBROUTINE lim_thd_lac( kideb, kiut )
       !!-------------------------------------------------------------------
       !!               ***   ROUTINE lim_thd_lac  ***
@@ -63,5 +63,4 @@
       !!      Fichefet T. and M. Maqueda 1997, J. Geo. Res., 102(C6), 
       !!                                                12609 -12646   
-      !!
       !! History :
       !!   1.0  !  01-04 (LIM)  original code
Index: /trunk/NEMO/LIM_SRC/limthd_zdf.F90
===================================================================
--- /trunk/NEMO/LIM_SRC/limthd_zdf.F90	(revision 87)
+++ /trunk/NEMO/LIM_SRC/limthd_zdf.F90	(revision 88)
@@ -1,11 +1,13 @@
 MODULE limthd_zdf
-#if defined key_ice_lim
    !!======================================================================
    !!                       ***  MODULE limthd_zdf ***
    !!                thermodynamic growth and decay of the ice 
    !!======================================================================
-
+#if defined key_ice_lim
+   !!----------------------------------------------------------------------
+   !!   'key_ice_lim'                                     LIM sea-ice model
    !!----------------------------------------------------------------------
    !!   lim_thd_zdf  : vertical accr./abl. and lateral ablation of sea ice
+   !!----------------------------------------------------------------------
    !! * Modules used
    USE par_oce          ! ocean parameters
@@ -25,8 +27,8 @@
    !! * Module variables
    REAL(wp)  ::           &  ! constant values
-      epsi20 = 1e-20   ,  &
-      epsi13 = 1e-13   ,  &
-      zzero  = 0.0     ,  &
-      zone   = 1.0
+      epsi20 = 1.e-20  ,  &
+      epsi13 = 1.e-13  ,  &
+      zzero  = 0.e0    ,  &
+      zone   = 1.e0
    !!----------------------------------------------------------------------
    !!   LIM 2.0 , UCL-LODYC-IPSL  (2003)
@@ -35,64 +37,64 @@
 
    SUBROUTINE lim_thd_zdf( kideb , kiut )
-       !!------------------------------------------------------------------
-       !!                ***  ROUTINE lim_thd_zdf  ***
-       !!              
-       !! ** Purpose : This routine determines the time evolution of snow 
-       !!      and sea-ice thicknesses, concentration and heat content 
-       !!      due to the vertical and lateral thermodynamic accretion-
-       !!      ablation processes. One only treats the case of lat. abl.
-       !!      For lateral accretion, see routine lim_lat_accr 
-       !! 
-       !! ** Method  : The representation of vertical growth and decay of 
-       !!      the sea-ice model is based upon the diffusion of heat 
-       !!      through the external and internal boundaries of a 
-       !!      three-layer system (two layers of ice and one layer and 
-       !!      one layer of snow, if present, on top of the ice).
-       !! 
-       !! ** Action  : - Calculation of some intermediates variables
-       !!              - Calculation of surface temperature
-       !!              - Calculation of available heat for surface ablation
-       !!              - Calculation of the changes in internal temperature 
-       !!                of the three-layer system, due to vertical diffusion 
-       !!                processes
-       !!              - Performs surface ablation and bottom accretion-ablation
-       !!              - Performs snow-ice formation
-       !!              - Performs lateral ablation
-       !!
-       !! References :
-       !!   Fichefet T. and M. Maqueda 1997, J. Geophys. Res., 102(C6), 12609-12646   
-       !!   Fichefet T. and M. Maqueda 1999, Clim. Dyn, 15(4), 251-268  
-       !!
-       !! History :
-       !!   original    : 01-04 (LIM)
-       !!   addition    : 02-08 (C. Ethe, G. Madec)
-       !!------------------------------------------------------------------
-       !! * Arguments
-       INTEGER , INTENT (in) ::  &
-          kideb ,  &  ! Start point on which the  the computation is applied
-          kiut        ! End point on which the  the computation is applied
-
-       !! * Local variables
-       INTEGER ::   ji       ! dummy loop indices
-
-       REAL(wp) , DIMENSION(jpij,2) ::  &
-          zqcmlt        ! energy due to surface( /1 ) and bottom melting( /2 )
-
-       REAL(wp), DIMENSION(jpij) ::  &
-          ztsmlt      &    ! snow/ice surface melting temperature
-          ,ztbif      &    ! int. temp. at the mid-point of the 1st layer of the snow/ice sys. 
-          ,zksn       &    ! effective conductivity of snow
-          ,zkic       &    ! effective conductivity of ice
-          ,zksndh     &    ! thermal cond. at the mid-point of the 1st layer of the snow/ice sys. 
-          , zfcsu     &    ! conductive heat flux at the surface of the snow/ice system 
-          , zfcsudt   &    ! = zfcsu * dt
-          , zi0       &    ! frac. of the net SW rad. which is not absorbed at the surface
-          , z1mi0     &    ! fraction of the net SW radiation absorbed at the surface
-          , zqmax     &    ! maximum energy stored in brine pockets
-          , zrcpdt    &    ! h_su*rho_su*cp_su/dt(h_su being the thick. of surf. layer)
-          , zts_old   &    ! previous surface temperature
-          , zidsn , z1midsn , zidsnic ! tempory variables
-
-       REAL(wp), DIMENSION(jpij) :: &
+      !!------------------------------------------------------------------
+      !!                ***  ROUTINE lim_thd_zdf  ***
+      !!              
+      !! ** Purpose : This routine determines the time evolution of snow 
+      !!      and sea-ice thicknesses, concentration and heat content 
+      !!      due to the vertical and lateral thermodynamic accretion-
+      !!      ablation processes. One only treats the case of lat. abl.
+      !!      For lateral accretion, see routine lim_lat_accr 
+      !! 
+      !! ** Method  : The representation of vertical growth and decay of 
+      !!      the sea-ice model is based upon the diffusion of heat 
+      !!      through the external and internal boundaries of a 
+      !!      three-layer system (two layers of ice and one layer and 
+      !!      one layer of snow, if present, on top of the ice).
+      !! 
+      !! ** Action  : - Calculation of some intermediates variables
+      !!              - Calculation of surface temperature
+      !!              - Calculation of available heat for surface ablation
+      !!              - Calculation of the changes in internal temperature 
+      !!                of the three-layer system, due to vertical diffusion 
+      !!                processes
+      !!              - Performs surface ablation and bottom accretion-ablation
+      !!              - Performs snow-ice formation
+      !!              - Performs lateral ablation
+      !!
+      !! References :
+      !!   Fichefet T. and M. Maqueda 1997, J. Geophys. Res., 102(C6), 12609-12646   
+      !!   Fichefet T. and M. Maqueda 1999, Clim. Dyn, 15(4), 251-268  
+      !!
+      !! History :
+      !!   original    : 01-04 (LIM)
+      !!   addition    : 02-08 (C. Ethe, G. Madec)
+      !!------------------------------------------------------------------
+      !! * Arguments
+      INTEGER , INTENT (in) ::  &
+         kideb ,  &  ! Start point on which the  the computation is applied
+         kiut        ! End point on which the  the computation is applied
+
+      !! * Local variables
+      INTEGER ::   ji       ! dummy loop indices
+
+      REAL(wp) , DIMENSION(jpij,2) ::  &
+         zqcmlt        ! energy due to surface( /1 ) and bottom melting( /2 )
+
+      REAL(wp), DIMENSION(jpij) ::  &
+         ztsmlt      &    ! snow/ice surface melting temperature
+         ,ztbif      &    ! int. temp. at the mid-point of the 1st layer of the snow/ice sys. 
+         ,zksn       &    ! effective conductivity of snow
+         ,zkic       &    ! effective conductivity of ice
+         ,zksndh     &    ! thermal cond. at the mid-point of the 1st layer of the snow/ice sys. 
+         , zfcsu     &    ! conductive heat flux at the surface of the snow/ice system 
+         , zfcsudt   &    ! = zfcsu * dt
+         , zi0       &    ! frac. of the net SW rad. which is not absorbed at the surface
+         , z1mi0     &    ! fraction of the net SW radiation absorbed at the surface
+         , zqmax     &    ! maximum energy stored in brine pockets
+         , zrcpdt    &    ! h_su*rho_su*cp_su/dt(h_su being the thick. of surf. layer)
+         , zts_old   &    ! previous surface temperature
+         , zidsn , z1midsn , zidsnic ! tempory variables
+
+      REAL(wp), DIMENSION(jpij) :: &
           zfnet       &  ! net heat flux at the top surface( incl. conductive heat flux)
           , zsprecip  &    ! snow accumulation
@@ -420,11 +422,11 @@
           zplediag(3) = 1 + 3.0 * sbeta * zkhic   
 
-          zsubdiag(1) = 0.0              
-          zsubdiag(2) = -1.0 * z1midsn(ji) * sbeta * zkhicint
-          zsubdiag(3) = -1.0 * sbeta * zkhic 
-
-          zsupdiag(1) =  -1.0 * z1midsn(ji) * sbeta * zkhsnint 
+          zsubdiag(1) =  0.e0              
+          zsubdiag(2) = -1.e0 * z1midsn(ji) * sbeta * zkhicint
+          zsubdiag(3) = -1.e0 * sbeta * zkhic 
+
+          zsupdiag(1) = -1.e0 * z1midsn(ji) * sbeta * zkhsnint 
           zsupdiag(2) = zsubdiag(3)
-          zsupdiag(3) = 0.0
+          zsupdiag(3) =  0.e0
           
           !     6.3. Fulfill the idependent term vector.
Index: /trunk/NEMO/LIM_SRC/limwri.F90
===================================================================
--- /trunk/NEMO/LIM_SRC/limwri.F90	(revision 87)
+++ /trunk/NEMO/LIM_SRC/limwri.F90	(revision 88)
@@ -32,22 +32,25 @@
 
    !! * Module variables
+   INTEGER, PARAMETER ::   &  !:
+      jpnoumax = 40             !: maximum number of variable for ice output
    INTEGER  ::                                &
       noumef                                     ! number of fields
-   REAL(wp)           , DIMENSION(noumax) ::  &
+   REAL(wp)           , DIMENSION(jpnoumax) ::  &
       cmulti ,                                &  ! multiplicative constant
       cadd                                       ! additive constant
-   CHARACTER(len = 35), DIMENSION(noumax) ::  &
+   CHARACTER(len = 35), DIMENSION(jpnoumax) ::  &
       titn                                       ! title of the field
-   CHARACTER(len = 8 ), DIMENSION(noumax) ::  &
+   CHARACTER(len = 8 ), DIMENSION(jpnoumax) ::  &
       nam                                        ! name of the field
-   CHARACTER(len = 8 ), DIMENSION(noumax) ::  &
+   CHARACTER(len = 8 ), DIMENSION(jpnoumax) ::  &
       uni                                        ! unit of the field
-   INTEGER            , DIMENSION(noumax) ::  &
+   INTEGER            , DIMENSION(jpnoumax) ::  &
       nc                                         ! switch for saving field ( = 1 ) or not ( = 0 )
 
    REAL(wp)  ::            &  ! constant values
-      epsi16 = 1e-16   ,  &
-      zzero  = 0.0     ,  &
-      zone   = 1.0
+      epsi16 = 1.e-16   ,  &
+      zzero  = 0.e0     ,  &
+      zone   = 1.e0
+   !!-------------------------------------------------------------------
 
 CONTAINS
@@ -70,5 +73,5 @@
          zindh,zinda,zindb,  &
          ztmu
-      REAL(wp), DIMENSION(jpi,jpj,noumax) :: &
+      REAL(wp), DIMENSION(jpi,jpj,jpnoumax) :: &
          zcmo
       REAL(wp), DIMENSION(jpi,jpj) ::  &
@@ -84,6 +87,4 @@
          ndex51  
       !!-------------------------------------------------------------------
- 
-
       
       IF ( numit == nstart ) THEN 
@@ -111,7 +112,7 @@
          CALL wheneq  ( jpij , tmask(:,:,1), 1, 1., ndex51, ndim)
          
-         DO jf = 1 , noumef
+         DO jf = 1, noumef
             IF ( nc(jf) == 1 ) THEN
-               CALL histdef( nice, nam(jf), titn(jf), uni(jf), jpi, jpj &
+               CALL histdef( nice, nam(jf), titn(jf), uni(jf), jpi, jpj   &
                   , nhorid, 1, 1, 1, -99, 32, clop, zsto, zout )
             ENDIF
@@ -133,5 +134,5 @@
       !-- calculs des valeurs instantanees
       
-      zcmo( 1:jpi, 1:jpj, 1:noumax ) = 0.0 
+      zcmo(:,:, 1:jpnoumax ) = 0.e0 
       DO jj = 2 , jpjm1
          DO ji = 2 , jpim1
@@ -229,5 +230,5 @@
          field_19
 
-      TYPE(FIELD) , DIMENSION(noumax) :: zfield
+      TYPE(FIELD) , DIMENSION(jpnoumax) :: zfield
 
       NAMELIST/namiceout/ noumef, &
Index: /trunk/NEMO/LIM_SRC/par_ice.F90
===================================================================
--- /trunk/NEMO/LIM_SRC/par_ice.F90	(revision 87)
+++ /trunk/NEMO/LIM_SRC/par_ice.F90	(revision 88)
@@ -7,23 +7,14 @@
    USE par_oce
 
-   INTEGER, PARAMETER ::   &  !:
-      kmax =   1    ,      &  !: ???
-      nsmax =  2              !: ???
+   IMPLICIT NONE
+   PUBLIC               ! allows par_oce and par_kind to be known in ice modules
 
-   !!-- Parameter providing the number of vertical ice layers
-   
-   INTEGER , PARAMETER ::      &  !:
-      nlayers   = 2 ,          &  !: ???
-      nlayersp1 = nlayers + 1     !: ???
+   INTEGER, PUBLIC, PARAMETER ::   &  !:
+      jpkmax =  1    ,      &  !: ???
+      jpsmax =  2              !: ???
 
-   ! maximum number of variable for output
-   INTEGER, PARAMETER ::   &  !:
-      noumax = 40             !: ???
-
-   ! Parameters for outputs to files "evolu" made by routine "informe" : ???
-   INTEGER, PARAMETER ::       &  !:
-      ninfmx = 100        ,    &  !: maximum number of key variables
-      nchinf = 5          ,    &  !: ???
-      nchsep = nchinf + 2         !: ???
+   INTEGER, PUBLIC, PARAMETER ::   &  !: 
+      jplayers   = 2 ,           &  !: number of vertical ice layers
+      jplayersp1 = jplayers + 1     !: ???
 
    !!======================================================================
Index: /trunk/NEMO/LIM_SRC/thd_ice.F90
===================================================================
--- /trunk/NEMO/LIM_SRC/thd_ice.F90	(revision 87)
+++ /trunk/NEMO/LIM_SRC/thd_ice.F90	(revision 88)
@@ -76,5 +76,5 @@
       dqla_ice_1d          !:    "                  "      dqla_ice
 
-   REAL(wp), PUBLIC, DIMENSION(jpij,nlayersp1) ::   &  !:
+   REAL(wp), PUBLIC, DIMENSION(jpij,jplayersp1) ::   &  !:
       tbif_1d              !: corresponding to the 2D var  tbif
 
Index: /trunk/NEMO/OPA_SRC/cpl_oce.F90
===================================================================
--- /trunk/NEMO/OPA_SRC/cpl_oce.F90	(revision 87)
+++ /trunk/NEMO/OPA_SRC/cpl_oce.F90	(revision 88)
@@ -234,8 +234,8 @@
       cpl_f_readflx,   &  !: coupler to ocean file name for flx.coupled
       cpl_f_readtau,   &  !: coupler to ocean file name for tau.coupled
-      cpl_f_writ   ,   &  !: ocean to coupler file name for stp_cmo
+      cpl_f_writ   ,   &  !: ocean to coupler file name for cpl_stp
       cpl_readflx  ,   &  !: coupler to ocean field name for flx.coupled
       cpl_readtau  ,   &  !: coupler to ocean field name for tau.coupled
-      cpl_writ            !: ocean to coupler field name for stp_cmo
+      cpl_writ            !: ocean to coupler field name for cpl_stp
 
   REAL(wp), DIMENSION(jpi,jpj) ::   &  !:
Index: /trunk/NEMO/OPA_SRC/geo2ocean.F90
===================================================================
--- /trunk/NEMO/OPA_SRC/geo2ocean.F90	(revision 87)
+++ /trunk/NEMO/OPA_SRC/geo2ocean.F90	(revision 88)
@@ -233,5 +233,5 @@
          pte, ptn, ptv
       REAL(wp), PARAMETER :: rpi = 3.141592653E0
-      REAL(wp), PARAMETER :: rad = rpi / 180.0E0
+      REAL(wp), PARAMETER :: rad = rpi / 180.e0
 
       !! * Local variables
Index: /trunk/NEMO/OPA_SRC/in_out_manager.F90
===================================================================
--- /trunk/NEMO/OPA_SRC/in_out_manager.F90	(revision 87)
+++ /trunk/NEMO/OPA_SRC/in_out_manager.F90	(revision 88)
@@ -64,11 +64,12 @@
       numwri     = 40 ,      &  !: logical unit for output write
       numisp     = 41 ,      &  !: logical unit for island statistics
+      numgap     = 45 ,      &  !: logical unit for differences diagnostic
       numwrs     = 46 ,      &  !: logical unit for output restart
       numtdt     = 62 ,      &  !: logical unit for data temperature
       numsdt     = 63 ,      &  !: logical unit for data salinity
+      numrnf     = 64 ,      &  !: logical unit for runoff data
       numwso     = 71 ,      &  !: logical unit for 2d output write
       numwvo     = 72 ,      &  !: logical unit for 3d output write
       numsst     = 65 ,      &  !: logical unit for surface temperature data
-      numgap     = 45 ,      &  !: logical unit for differences diagnostic
       numbol     = 67 ,      &  !: logical unit for "bol" diagnostics
       numptr     = 68 ,      &  !: logical unit for Poleward TRansports
Index: /trunk/NEMO/OPA_SRC/lbclnk.F90
===================================================================
--- /trunk/NEMO/OPA_SRC/lbclnk.F90	(revision 87)
+++ /trunk/NEMO/OPA_SRC/lbclnk.F90	(revision 88)
@@ -72,4 +72,6 @@
          !             !   =-1 , the sign is changed if north fold boundary
          !             !   = 1 , no sign change
+         !             !   = 0 , no sign change and > 0 required (use the inner
+         !             !         row/column if closed boundary)
 
       !! * Local declarations
@@ -89,19 +91,14 @@
 
          CASE ( 1 , 4 , 6 )                    ! * cyclic east-west
-
             pt3d( 1 ,:,jk) = pt3d(jpim1,:,jk)          ! all points
             pt3d(jpi,:,jk) = pt3d(  2  ,:,jk)
 
          CASE DEFAULT                          ! * closed
-
             SELECT CASE ( cd_type )
-
             CASE ( 'T' , 'U' , 'V' , 'W' )             ! T-, U-, V-, W-points
                pt3d( 1 ,:,jk) = 0.e0
                pt3d(jpi,:,jk) = 0.e0
-
             CASE ( 'F' )                               ! F-point
                pt3d(jpi,:,jk) = 0.e0
-
             END SELECT
 
@@ -115,16 +112,13 @@
 
             SELECT CASE ( cd_type )
-
             CASE ( 'T' , 'U' , 'W' )                   ! T-, U-, W-points
                pt3d(:, 1 ,jk) = pt3d(:,3,jk)
                pt3d(:,jpj,jk) = 0.e0
-
             CASE ( 'V' , 'F' )                         ! V-, F-points
                pt3d(:, 1 ,jk) = psgn * pt3d(:,2,jk)
                pt3d(:,jpj,jk) = 0.e0
-
             END SELECT
 
-         CASE ( 3 , 4 )                            ! *  North fold  T-point pivot
+         CASE ( 3 , 4 )                        ! *  North fold  T-point pivot
 
             pt3d( 1 ,jpj,jk) = 0.e0
@@ -132,5 +126,4 @@
 
             SELECT CASE ( cd_type )
-
             CASE ( 'T' , 'W' )                         ! T-, W-point
                DO ji = 2, jpi
@@ -143,5 +136,4 @@
                   pt3d(ji,jpjm1,jk) = psgn * pt3d(ijt,jpjm1,jk)
                END DO
-
             CASE ( 'U' )                               ! U-point
                DO ji = 1, jpi-1
@@ -154,5 +146,4 @@
                   pt3d(ji,jpjm1,jk) = psgn * pt3d(iju,jpjm1,jk)
                END DO
-
             CASE ( 'V' )                               ! V-point
                   DO ji = 2, jpi
@@ -162,5 +153,4 @@
                      pt3d(ji,jpj  ,jk) = psgn * pt3d(ijt,jpj-3,jk)
                   END DO
-
             CASE ( 'F' )                               ! F-point
                   DO ji = 1, jpi-1
@@ -169,5 +159,4 @@
                      pt3d(ji,jpj  ,jk) = pt3d(iju,jpj-3,jk)
                   END DO
-
             END SELECT
 
@@ -178,5 +167,4 @@
 
             SELECT CASE ( cd_type )
-
             CASE ( 'T' , 'W' )                         ! T-, W-point
                DO ji = 1, jpi
@@ -185,5 +173,4 @@
                   pt3d(ji,jpj,jk) = psgn * pt3d(ijt,jpj-1,jk)
                END DO
-
             CASE ( 'U' )                               ! U-point
                   DO ji = 1, jpi-1
@@ -192,5 +179,4 @@
                      pt3d(ji,jpj,jk) = psgn * pt3d(iju,jpj-1,jk)
                   END DO
-
             CASE ( 'V' )                               ! V-point
                   DO ji = 1, jpi
@@ -203,5 +189,4 @@
                      pt3d(ji,jpjm1,jk) = psgn * pt3d(ijt,jpjm1,jk)
                   END DO
-
             CASE ( 'F' )                               ! F-point
                   DO ji = 1, jpi-1
@@ -215,15 +200,12 @@
             END SELECT
 
-         CASE DEFAULT                           ! *  closed
+         CASE DEFAULT                          ! *  closed
 
             SELECT CASE ( cd_type )
-
             CASE ( 'T' , 'U' , 'V' , 'W' )             ! T-, U-, V-, W-points
                pt3d(:, 1 ,jk) = 0.e0
                pt3d(:,jpj,jk) = 0.e0
-
             CASE ( 'F' )                               ! F-point
                pt3d(:,jpj,jk) = 0.e0
-
             END SELECT
 
@@ -273,19 +255,17 @@
 
       CASE ( 1 , 4 , 6 )                       ! * cyclic east-west
-
          pt2d( 1 ,:) = pt2d(jpim1,:)
          pt2d(jpi,:) = pt2d(  2  ,:)
 
       CASE DEFAULT                             ! * closed 
-
          SELECT CASE ( cd_type )
-
          CASE ( 'T' , 'U' , 'V' , 'W' )                ! T-, U-, V-, W-points
             pt2d( 1 ,:) = 0.e0
             pt2d(jpi,:) = 0.e0
-
-         CASE ( 'F' , 'I' )                            ! F-point, ice U-V point
+         CASE ( 'F' )                                  ! F-point, ice U-V point
             pt2d(jpi,:) = 0.e0 
-
+         CASE ( 'I' )                                  ! F-point, ice U-V point
+            pt2d( 1 ,:) = 0.e0 
+            pt2d(jpi,:) = 0.e0 
          END SELECT
 
@@ -299,13 +279,10 @@
 
          SELECT CASE ( cd_type )
-
          CASE ( 'T' , 'U' , 'W' )                      ! T-, U-, W-points
             pt2d(:, 1 ) = pt2d(:,3)
             pt2d(:,jpj) = 0.e0
-
-         CASE ( 'V' , 'F' )                            ! V-, F-points, ice U-V point
+         CASE ( 'V' , 'F' , 'I' )                      ! V-, F-points, ice U-V point
             pt2d(:, 1 ) = psgn * pt2d(:,2)
             pt2d(:,jpj) = 0.e0
-
          END SELECT
 
@@ -419,16 +396,12 @@
 
          SELECT CASE ( cd_type )
-
          CASE ( 'T' , 'U' , 'V' , 'W' )                ! T-, U-, V-, W-points
             pt2d(:, 1 ) = 0.e0
             pt2d(:,jpj) = 0.e0
-
          CASE ( 'F' )                                  ! F-point
             pt2d(:,jpj) = 0.e0
-
          CASE ( 'I' )                                  ! ice U-V point
             pt2d(:, 1 ) = 0.e0
             pt2d(:,jpj) = 0.e0
-
          END SELECT
 
Index: /trunk/NEMO/OPA_SRC/lib_cray.f90
===================================================================
--- /trunk/NEMO/OPA_SRC/lib_cray.f90	(revision 87)
+++ /trunk/NEMO/OPA_SRC/lib_cray.f90	(revision 88)
@@ -9,5 +9,5 @@
 !  isrchne
 !---------------------------------------------------------
-        FUNCTION sdot( I, X, J, Y, K )
+     FUNCTION sdot( I, X, J, Y, K )
         DIMENSION X(1), Y(1)
         SDOT = 0.
@@ -15,7 +15,7 @@
         SDOT = SDOT + X(1+(N-1)*J) * Y(1+(N-1)*K)
         END DO
-        END FUNCTION sdot
+     END FUNCTION sdot
 !---------------------------------------------------------
-        SUBROUTINE wheneq ( i, x, j, t, ind, nn )
+     SUBROUTINE wheneq ( i, x, j, t, ind, nn )
         IMPLICIT NONE
 
@@ -35,24 +35,22 @@
         END DO 
 
-        END SUBROUTINE wheneq
+     END SUBROUTINE wheneq
 !---------------------------------------------------------
-        SUBROUTINE SAXPY(I,A,X,J,Y,K)
+     SUBROUTINE saxpy( I, A, X, J, Y, K )
         DIMENSION X(1),Y(1)
-        DO 1 N=1,I
-        Y(1+(N-1)*K)=A*X(1+(N-1)*J)+Y(1+(N-1)*K)
-  1     CONTINUE
-        RETURN
-        END
+        DO N = 1, I
+           Y(1+(N-1)*K)=A*X(1+(N-1)*J)+Y(1+(N-1)*K)
+        END DO
+     END SUBROUTINE saxpy
 !---------------------------------------------------------
-        FUNCTION ISRCHNE(K,X,I,B)
+     FUNCTION isrchne( K, X, I, B )
         DIMENSION X(1)
-        DO 1 N=1,K
-        IF(X(1+(N-1)*I) /= B)THEN
-         ISRCHNE=N
-        RETURN
-        ELSE
-          ISRCHNE=N+1
-        ENDIF
- 1      CONTINUE
-        RETURN
-        END
+        DO N = 1, K
+           IF( X(1+(N-1)*I) /= B ) THEN
+              ISRCHNE = N
+              RETURN
+           ELSE
+              ISRCHNE = N + 1
+           ENDIF
+        END DO
+     END FUNCTION isrchne
Index: /trunk/NEMO/OPA_SRC/lib_isml.f90
===================================================================
--- /trunk/NEMO/OPA_SRC/lib_isml.f90	(revision 87)
+++ /trunk/NEMO/OPA_SRC/lib_isml.f90	(revision 88)
@@ -14,7 +14,7 @@
 !
 !---------------------------------------------------------
-      SUBROUTINE linrg(kn,pa,klda,painv,kldainv)
-
-! compute inverse matrix
+   SUBROUTINE linrg(kn,pa,klda,painv,kldainv)
+
+      !! compute inverse matrix
 
       IMPLICIT NONE 
@@ -26,23 +26,23 @@
       INTEGER ji
 
-      IF (kn /= klda.or.kn /= kldainv) THEN 
+      IF( kn /= klda .OR. kn /= kldainv ) THEN 
           write(0,*)'change your parameters'
           STOP 
       ENDIF 
 
-      CALL vmov(kn*kn,pa,painv)
-
-      CALL gauss(kn,painv,iplin,zv)
-
-      zb(:,:) = 0.0   
-      DO ji=1,kn
-        zb(ji,ji)=1.
-        CALL desremopt(kn,painv,iplin,zb(1,ji),zb(1,ji),zv)
+      CALL vmov( kn*kn, pa, painv )
+
+      CALL gauss( kn, painv, iplin, zv )
+
+      zb(:,:) = 0.e0   
+      DO ji = 1, kn
+        zb(ji,ji) = 1.e0
+        CALL desremopt( kn, painv, iplin, zb(1,ji), zb(1,ji), zv )
       END DO 
-      CALL vmov(kn*kn,zb,painv)
-
-      END SUBROUTINE linrg
-!---------------------------------------------------------
-      SUBROUTINE gauss(kn,pa,kplin,pv)
+      CALL vmov( kn*kn, zb, painv )
+
+   END SUBROUTINE linrg
+!---------------------------------------------------------
+   SUBROUTINE gauss(kn,pa,kplin,pv)
 
       IMPLICIT NONE 
@@ -108,20 +108,19 @@
       END DO 
 
-      END SUBROUTINE gauss
-!---------------------------------------------------------
-        FUNCTION isamax( I, X )
-        DIMENSION X(I)
-        ISAMAX=0
-        XMIN=-1E+50
-        DO 1 N=1,I
-        IF(ABS(X(N)) > XMIN) THEN
-        XMIN=X(N)
-        ISAMAX=N
-        ENDIF
-  1     CONTINUE
-        RETURN
-        END
-!---------------------------------------------------------
-      SUBROUTINE vmov(kn,px,py)
+   END SUBROUTINE gauss
+!---------------------------------------------------------
+   FUNCTION isamax( I, X )
+      DIMENSION X(I)
+      ISAMAX = 0
+      XMIN = -1e+50
+      DO N = 1, I
+         IF(ABS(X(N)) > XMIN ) THEN
+            XMIN = X(N)
+            ISAMAX = N
+         ENDIF
+      END DO
+   END FUNCTION isamax
+!---------------------------------------------------------
+   SUBROUTINE vmov(kn,px,py)
 
       IMPLICIT NONE 
@@ -131,11 +130,10 @@
 
       DO ji=1,kn
-        py(ji)=px(ji)
+         py(ji)=px(ji)
       END DO 
 
-      RETURN 
-      END 
-!---------------------------------------------------------
-      subroutine desremopt(n,a,plin,y,x,v)
+   END SUBROUTINE vmov
+!---------------------------------------------------------
+   subroutine desremopt(n,a,plin,y,x,v)
       implicit none
       integer n,i,  j0
@@ -176,7 +174,7 @@
       end do
 
-      end SUBROUTINE desremopt
-!---------------------------------------------------------
-      SUBROUTINE DTRSV ( UPLO, TRANS, DIAG, N, A, LDA, X, INCX )
+   end SUBROUTINE desremopt
+!---------------------------------------------------------
+   SUBROUTINE DTRSV ( UPLO, TRANS, DIAG, N, A, LDA, X, INCX )
 !!    .. Scalar Arguments ..
       INTEGER            INCX, LDA, N
@@ -457,11 +455,7 @@
       END IF
 
-      RETURN
-
-!!    End of DTRSV .
-
-      END
-!---------------------------------------------------------
-      SUBROUTINE DGER  ( M, N, ALPHA, X, INCX, Y, INCY, A, LDA )
+   END SUBROUTINE DTRSV
+!---------------------------------------------------------
+   SUBROUTINE DGER  ( M, N, ALPHA, X, INCX, Y, INCY, A, LDA )
 !!    .. Scalar Arguments ..
 !     DOUBLE PRECISION   ALPHA
@@ -620,11 +614,7 @@
       END IF
 
-      RETURN
-
-!!    End of DGER  .
-
-      END
-!---------------------------------------------------------
-      SUBROUTINE XERBLA ( SRNAME, INFO )
+   END SUBROUTINE DGER
+!---------------------------------------------------------
+   SUBROUTINE XERBLA ( SRNAME, INFO )
 !!    ..    Scalar Arguments ..
       INTEGER            INFO
@@ -668,9 +658,7 @@
                ' had an illegal value' )
 
-!!    End of XERBLA.
-
-      END
+   END SUBROUTINE XERBLA
 !-----------------------------------------------------------
-      FUNCTION lsame( c1, c2 )
+   FUNCTION lsame( c1, c2 )
       logical lsame
       CHARACTER (len=*), INTENT(in) ::   c1, c2
@@ -680,4 +668,3 @@
           lsame=.FALSE.
       ENDIF
-
-      END FUNCTION lsame
+   END FUNCTION lsame
Index: /trunk/NEMO/OPA_SRC/mppini.F90
===================================================================
--- /trunk/NEMO/OPA_SRC/mppini.F90	(revision 87)
+++ /trunk/NEMO/OPA_SRC/mppini.F90	(revision 88)
@@ -9,4 +9,5 @@
    !!   mpp_init2      : Lay out the global domain over processors 
    !!                    with land processor elimination
+   !!   mpp_init_ioispl: IOIPSL initialization in mpp
    !!----------------------------------------------------------------------
    !! * Modules used
@@ -66,4 +67,5 @@
          WRITE(numout,*) '~~~~~~~~~~~: '
          WRITE(numout,*) '         nperio = ', nperio
+         WRITE(numout,*) '         npolj  = ', npolj
          WRITE(numout,*) '         nimpp  = ', nimpp
          WRITE(numout,*) '         njmpp  = ', njmpp
@@ -140,7 +142,7 @@
 
 #if defined key_mpp_shmem
-      IF(lwp)WRITE(numout,*)
-      IF(lwp)WRITE(numout,*) 'mpp_init : Message Passing PVM T3E + SHMEM'
-      IF(lwp)WRITE(numout,*) '~~~~~~~~'
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*) 'mpp_init : Message Passing PVM T3E + SHMEM'
+      IF(lwp) WRITE(numout,*) '~~~~~~~~'
 
       CALL mppshmem                           ! Initialisation of shmem array
@@ -148,7 +150,7 @@
 #endif
 #if defined key_mpp_mpi
-      IF(lwp)WRITE(numout,*)
-      IF(lwp)WRITE(numout,*) 'mpp_init : Message Passing MPI'
-      IF(lwp)WRITE(numout,*) '~~~~~~~~'
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*) 'mpp_init : Message Passing MPI'
+      IF(lwp) WRITE(numout,*) '~~~~~~~~'
 #endif
 
@@ -472,11 +474,72 @@
       END IF
 
+      ! Prepare NetCDF output file (if necessary)
+      CALL mpp_init_ioipsl
 
    END SUBROUTINE mpp_init
 
-
 #  include "mppini_2.h90"
 
+# if defined key_fdir   ||   defined key_dimgout
+   !!----------------------------------------------------------------------
+   !!   'key_fdir' OR 'key_dimgout'                  NO use of NetCDF files
+   !!----------------------------------------------------------------------
+   SUBROUTINE mpp_init_ioipsl       ! Dummy routine
+   END SUBROUTINE mpp_init_ioipsl  
+# else
+   SUBROUTINE mpp_init_ioipsl
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE mpp_init_ioipsl  ***
+      !!
+      !! ** Purpose :   
+      !!
+      !! ** Method  :   
+      !!
+      !! History :
+      !!   9.0  !  04-03  (G. Madec)  MPP-IOIPSL 
+      !!----------------------------------------------------------------------
+      USE ioipsl
+      INTEGER, DIMENSION(4) ::   &
+         iglo, iloc, iabsf, iabsl, ihals, ihale   ! ???
+      !!----------------------------------------------------------------------
+      iglo(1) = jpiglo
+      iglo(2) = jpjglo
+      iglo(3) = jpk
+      iglo(4) = 1
+      iloc(1) = nlci
+      iloc(2) = nlcj
+      iloc(3) = jpk
+      iloc(4) = 1
+      iabsf(1) = nimppt(narea)
+      iabsf(2) = njmppt(narea)
+      iabsf(3) = 1
+      iabsf(4) = 1
+      iabsl(:) = iabsf(:) + iloc(:) - 1
+      ihals(1) = jpreci
+      ihals(2) = jprecj
+      ihals(3) = 0
+      ihals(4) = 0
+      ihale(1) = jpreci
+      ihale(2) = jprecj
+      ihale(3) = 0
+      ihale(4) = 0
+      IF( nbondi == -1 .OR. nbondi == 2 )   ihals(1) = 0
+      IF( nbondi ==  1 .OR. nbondi == 2 )   ihale(1) = 0
+      IF( nbondj == -1 .OR. nbondj == 2 )   ihals(2) = 0
+      IF( nbondj ==  1 .OR. nbondj == 2 )   ihale(2) = 0
+      IF(lwp) THEN
+          WRITE(numout,*) 'mpp_init_ioipsl :   iloc  = ', iloc (1), iloc (2), iloc (3), iloc (4)
+          WRITE(numout,*) '~~~~~~~~~~~~~~~     iabsf = ', iabsf(1), iabsf(2), iabsf(3), iabsf(4)
+          WRITE(numout,*) '                    ihals = ', ihals(1), ihals(2), ihals(3), ihals(4)
+          WRITE(numout,*) '                    ihale = ', ihale(1), ihale(2), ihale(3), ihale(4)
+      ENDIF
+
+      CALL ioipsl_inimpp( jpnij, nproc, iglo, iloc, iabsf, iabsl, ihals, ihale )
+
+   END SUBROUTINE mpp_init_ioipsl  
+
+# endif
 #endif
+
    !!======================================================================
 END MODULE mppini
Index: /trunk/NEMO/OPA_SRC/oce.F90
===================================================================
--- /trunk/NEMO/OPA_SRC/oce.F90	(revision 87)
+++ /trunk/NEMO/OPA_SRC/oce.F90	(revision 88)
@@ -22,5 +22,6 @@
    !! --------------------------
    REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   &   !:
-      ! before !  now      !  after  !
+      ! before !  now      !  after  !      ! the after trends becomes the fields
+      ! fields !  fields   !  trends !      ! only in dyn(tra)_zdf and dyn(tra)_nxt
       ub       ,  un       ,  ua     ,   &  !: i-horizontal velocity (m/s)
       vb       ,  vn       ,  va     ,   &  !: j-horizontal velocity (m/s)
Index: /trunk/NEMO/OPA_SRC/opa.F90
===================================================================
--- /trunk/NEMO/OPA_SRC/opa.F90	(revision 87)
+++ /trunk/NEMO/OPA_SRC/opa.F90	(revision 88)
@@ -131,5 +131,4 @@
       IF(lwp) WRITE(numout,cform_aaa)       ! Flag AAAAAAA
 
-      CALL opa_flg
                                             ! Domain decomposition
       IF( jpni*jpnj == jpnij ) THEN
@@ -173,7 +172,11 @@
 #endif
 
+      !                                     ! Ocean scheme
+
+      CALL opa_flg                              ! Choice of algorithms
+
       !                                     ! Ocean physics
 
-      CALL tra_qsr_init                         ! Salor radiation penetration
+      CALL tra_qsr_init                         ! Solar radiation penetration
 
       CALL ldf_dyn_init                         ! Lateral ocean momentum physics
Index: /trunk/NEMO/OPA_SRC/phycst.F90
===================================================================
--- /trunk/NEMO/OPA_SRC/phycst.F90	(revision 87)
+++ /trunk/NEMO/OPA_SRC/phycst.F90	(revision 88)
@@ -91,5 +91,5 @@
       !!----------------------------------------------------------------------
       !! * Local variables
-      CHARACTER (len=64) ::   cform = "( 4(A13, I7) )" 
+      CHARACTER (len=64) ::   cform = "(A9, 3(A13, I7) )" 
       !!----------------------------------------------------------------------
 
@@ -97,42 +97,35 @@
       IF(lwp) WRITE(numout,*) ' phy_cst : initialization of ocean parameters and constants'
       IF(lwp) WRITE(numout,*) ' ~~~~~~~'
-      IF(lwp) WRITE(numout,*)
-
 
       ! Ocean Parameters
       ! ----------------
       IF(lwp) THEN
+         WRITE(numout,*) '       parameter file'
          WRITE(numout,*)
-         WRITE(numout,*) ' parameter file'
-         WRITE(numout,*) ' ~~~~~~~~~~~~~~'
+         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,*) '          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,*) ' dimension of model'
-         WRITE(numout,*) ' local domain                              Global domain           Data domain '
-         WRITE(numout,cform) '   jpi     : ', jpi, '   jpim1   : ', jpim1, '   jpiglo  : ', jpiglo, '   jpidta  : ', jpidta
-         WRITE(numout,cform) '   jpj     : ', jpj, '   jpjm1   : ', jpjm1, '   jpjglo  : ', jpjglo, '   jpjdta  : ', jpjdta
-         WRITE(numout,cform) '   jpk     : ', jpk, '   jpkm1   : ', jpkm1, '   jpk     : ', jpk   , '   jpkdta  : ', jpkdta
-         WRITE(numout,*)      '  jpij    : ', jpij
-         WRITE(numout,*)
-         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,*) '          lateral domain boundary condition type : jperio  = ', jperio
+         WRITE(numout,*) '          domain island (use in rigid-lid case)  : jpisl   = ', jpisl 
+         WRITE(numout,*) '                                                   jpnisl  = ', jpnisl
       ENDIF
-
 
       ! Define constants
       ! ----------------
       IF(lwp) WRITE(numout,*)
-      IF(lwp) WRITE(numout,*) ' constants'
-      IF(lwp) WRITE(numout,*) ' ~~~~~~~~~'
+      IF(lwp) WRITE(numout,*) '       constants'
 
       IF(lwp) WRITE(numout,*)
-      IF(lwp) WRITE(numout,*) ' mathematical constant rpi = ', rpi
+      IF(lwp) WRITE(numout,*) '          mathematical constant                 rpi = ', rpi
 
       rsiyea = 365.25 * rday * 2. * rpi / 6.283076
@@ -140,58 +133,58 @@
       omega  = 2. * rpi / rsiday 
       IF(lwp) WRITE(numout,*)
-      IF(lwp) WRITE(numout,*) ' day          rday   = ', rday,   ' s'
-      IF(lwp) WRITE(numout,*) ' sideral year rsiyea = ', rsiyea, ' s'
-      IF(lwp) WRITE(numout,*) ' sideral day  rsiday = ', rsiday, ' s'
-      IF(lwp) WRITE(numout,*) ' omega        omega  = ', omega,  ' s-1'
+      IF(lwp) WRITE(numout,*) '          day                                rday   = ', rday,   ' s'
+      IF(lwp) WRITE(numout,*) '          sideral year                       rsiyea = ', rsiyea, ' s'
+      IF(lwp) WRITE(numout,*) '          sideral day                        rsiday = ', rsiday, ' s'
+      IF(lwp) WRITE(numout,*) '          omega                              omega  = ', omega,  ' s-1'
 
       rjjss = rjjhh * rhhmm * rmmss
       IF(lwp) WRITE(numout,*)
-      IF(lwp) WRITE(numout,*) ' nb of months per year    raamo = ', raamo, ' months'
-      IF(lwp) WRITE(numout,*) ' nb of hours per day      rjjhh = ', rjjhh, ' hours'
-      IF(lwp) WRITE(numout,*) ' nb of minutes per hour   rhhmm = ', rhhmm, ' mn'
-      IF(lwp) WRITE(numout,*) ' nb of seconds per minute rmmss = ', rmmss, ' s'
-      IF(lwp) WRITE(numout,*) ' nb of seconds per day    rjjss = ', rjjss, ' s'
+      IF(lwp) WRITE(numout,*) '          nb of months per year               raamo = ', raamo, ' months'
+      IF(lwp) WRITE(numout,*) '          nb of hours per day                 rjjhh = ', rjjhh, ' hours'
+      IF(lwp) WRITE(numout,*) '          nb of minutes per hour              rhhmm = ', rhhmm, ' mn'
+      IF(lwp) WRITE(numout,*) '          nb of seconds per minute            rmmss = ', rmmss, ' s'
+      IF(lwp) WRITE(numout,*) '          nb of seconds per day               rjjss = ', rjjss, ' s'
 
       IF(lwp) WRITE(numout,*)
-      IF(lwp) WRITE(numout,*) ' earth radius ra   = ', ra, ' m'
-      IF(lwp) WRITE(numout,*) ' gravity      grav = ', grav , ' m/s2'
+      IF(lwp) WRITE(numout,*) '          earth radius                         ra   = ', ra, ' m'
+      IF(lwp) WRITE(numout,*) '          gravity                              grav = ', grav , ' m/s^2'
 
       IF(lwp) WRITE(numout,*)
-      IF(lwp) WRITE(numout,*) ' triple point of temperature      rtt      = ', rtt     , ' K'
-      IF(lwp) WRITE(numout,*) ' freezing point of water ( C)     rt0      = ', rt0     , ' K'
-      IF(lwp) WRITE(numout,*) ' melting point of snow            rt0_snow = ', rt0_snow, ' K'
-      IF(lwp) WRITE(numout,*) ' melting point of ice             rt0_ice  = ', rt0_ice , ' K'
+      IF(lwp) WRITE(numout,*) '          triple point of temperature      rtt      = ', rtt     , ' K'
+      IF(lwp) WRITE(numout,*) '          freezing point of water          rt0      = ', rt0     , ' K'
+      IF(lwp) WRITE(numout,*) '          melting point of snow            rt0_snow = ', rt0_snow, ' K'
+      IF(lwp) WRITE(numout,*) '          melting point of ice             rt0_ice  = ', rt0_ice , ' K'
 
-      ro0cpr = 1. / ( rau0 * rcp )    ! 
+      ro0cpr = 1. / ( rau0 * rcp )
       IF(lwp) WRITE(numout,*)
-      IF(lwp) WRITE(numout,*) ' volumic mass of pure water (kg/m3)  rauw   = ', rauw, ' kg/m3'
-      IF(lwp) WRITE(numout,*) ' volumic mass of reference  (kg/m3)  rau0   = ', rau0, ' kg/m3'
-      IF(lwp) WRITE(numout,*) ' ocean specific heat                 rcp    = ', rcp
-      IF(lwp) WRITE(numout,*) '               1. / ( rau0 * rcp ) = ro0cpr = ', ro0cpr
+      IF(lwp) WRITE(numout,*) '          volumic mass of pure water         rauw   = ', rauw, ' kg/m^3'
+      IF(lwp) WRITE(numout,*) '          volumic mass of reference          rau0   = ', rau0, ' kg/m^3'
+      IF(lwp) WRITE(numout,*) '          ocean specific heat                rcp    = ', rcp
+      IF(lwp) WRITE(numout,*) '                       1. / ( rau0 * rcp ) = ro0cpr = ', ro0cpr
 
       IF(lwp) THEN
          WRITE(numout,*)
-         WRITE(numout,*) ' thermal conductivity of the snow          = ', rcdsn   , ' J.s-1.m-1.K-1'
-         WRITE(numout,*) ' thermal conductivity of the ice           = ', rcdic   , ' J.s-1.m-1.K-1'
-         WRITE(numout,*) ' density times specific heat for snow      = ', rcpsn   , ' J.m-3.K-1' 
-         WRITE(numout,*) ' density times specific heat for ice       = ', rcpic   , ' J.m-3.K-1' 
-         WRITE(numout,*) ' volumetric latent heat fusion of sea ice  = ', xlic    , ' J.m-3' 
-         WRITE(numout,*) ' volumetric latent heat fusion of snow     = ', xlsn    , ' J.m-3' 
-         WRITE(numout,*) ' latent heat of sublimation of snow        = ', xsn     , ' J.kg-1' 
-         WRITE(numout,*) ' density of sea ice                        = ', rhoic   , ' kg.m-3'
-         WRITE(numout,*) ' density of snow                           = ', rhosn   , ' kg.m-3'
-         WRITE(numout,*) ' emissivity of snow or ice                 = ', emic  
-         WRITE(numout,*) ' salinity of ice                           = ', sice    , ' psu'
-         WRITE(numout,*) ' salinity of sea                           = ', soce    , ' psu'
-         WRITE(numout,*) ' latent heat of evaporation (water)        = ', cevap   , ' J.m-3' 
-         WRITE(numout,*) ' correction factor for solar radiation     = ', srgamma 
-         WRITE(numout,*) ' von Karman constant                       = ', vkarmn 
-         WRITE(numout,*) ' Stefan-Boltzmann constant                 = ', stefan  , ' J.s-1.m-2.K-4'
+         WRITE(numout,*) '          thermal conductivity of the snow          = ', rcdsn   , ' J/s/m/K'
+         WRITE(numout,*) '          thermal conductivity of the ice           = ', rcdic   , ' J/s/m/K'
+         WRITE(numout,*) '          density times specific heat for snow      = ', rcpsn   , ' J/m^3/K' 
+         WRITE(numout,*) '          density times specific heat for ice       = ', rcpic   , ' J/m^3/K'
+         WRITE(numout,*) '          volumetric latent heat fusion of sea ice  = ', xlic    , ' J/m' 
+         WRITE(numout,*) '          volumetric latent heat fusion of snow     = ', xlsn    , ' J/m' 
+         WRITE(numout,*) '          latent heat of sublimation of snow        = ', xsn     , ' J/kg' 
+         WRITE(numout,*) '          density of sea ice                        = ', rhoic   , ' kg/m^3'
+         WRITE(numout,*) '          density of snow                           = ', rhosn   , ' kg/m^3'
+         WRITE(numout,*) '          emissivity of snow or ice                 = ', emic  
+         WRITE(numout,*) '          salinity of ice                           = ', sice    , ' psu'
+         WRITE(numout,*) '          salinity of sea                           = ', soce    , ' psu'
+         WRITE(numout,*) '          latent heat of evaporation (water)        = ', cevap   , ' J/m^3' 
+         WRITE(numout,*) '          correction factor for solar radiation     = ', srgamma 
+         WRITE(numout,*) '          von Karman constant                       = ', vkarmn 
+         WRITE(numout,*) '          Stefan-Boltzmann constant                 = ', stefan  , ' J/s/m^2/K^4'
 
          WRITE(numout,*)
-         WRITE(numout,*) ' conversion: degre ==> radian rad = ', rad
+         WRITE(numout,*) '          conversion: degre ==> radian          rad = ', rad
 
          WRITE(numout,*)
-         WRITE(numout,*) ' smallest real computer value= ', rsmall
+         WRITE(numout,*) '          smallest real computer value       rsmall = ', rsmall
       ENDIF
 
