Index: trunk/NEMO/OPA_SRC/DTA/dtasal.F90
===================================================================
--- trunk/NEMO/OPA_SRC/DTA/dtasal.F90	(revision 15)
+++ trunk/NEMO/OPA_SRC/DTA/dtasal.F90	(revision 16)
@@ -23,7 +23,7 @@
    
    !! * Shared module variables
-   LOGICAL , PUBLIC, PARAMETER ::   lk_dtasal = .TRUE.    ! salinity data flag
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   &
-      s_dta       ! salinity data at given time-step
+   LOGICAL , PUBLIC, PARAMETER ::   lk_dtasal = .TRUE.    !: salinity data flag
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   &  !:
+      s_dta       !: salinity data at given time-step
 
    !! * Module variables
@@ -81,6 +81,7 @@
       INTEGER ::   ji, jj, jk, jl    ! dummy loop indicies
       INTEGER ::   &
-         imois, iman, ik, i15,  &  ! temporary integers
-         ipi, ipj, ipk, itime      !    "          "
+         imois, iman, ik, i15,       &  ! temporary integers
+         ipi, ipj, ipk, itime,       &  !    "          "
+         il0, il1, ii0, ii1, ij0, ij1   !    "          "
       INTEGER, DIMENSION(jpmois) ::   istep
       REAL(wp) ::   &
@@ -186,7 +187,8 @@
             !                                        !  ORCA_R2 configuration
             !                                        ! =======================
-         
-            DO jj = mj0(101), mj1(109)                      ! Reduced salinity in the Alboran Sea
-               DO ji = mi0(141), mi1(155)
+            ij0 = 101   ;   ij1 = 109
+            ii0 = 141   ;   ii1 = 155   
+            DO jj = mj0(ij0), mj1(ij1)                      ! Reduced salinity in the Alboran Sea
+               DO ji = mi0(ii0), mi1(ii1)
                   DO jk = 13, 13
                      saldta(ji,jj,jk,:) = saldta(ji,jj,jk,:) - 0.15
@@ -205,11 +207,18 @@
             IF( n_cla == 1 ) THEN 
                !                                         ! New salinity profile at Gibraltar
-               saldta( mi0(139):mi1(139) , mj0(101):mj1(101) , : , : ) =   &
-                  &                                    saldta( mi0(138):mi1(138) , mj0(101):mj1(101) , : , : )
-               saldta( mi0(139):mi1(139) , mj0(102):mj1(102) , : , : ) =   &
-                  &                                    saldta( mi0(138):mi1(138) , mj0(102):mj1(102) , : , : )
-               DO jl = mi0(138), mi1(138)                ! New temperature profile at Gibraltar
-                  DO jj = mj0(101), mj1(102)
-                     DO ji = mi0(139), mi1(139)
+               il0 = 138   ;   il1 = 138   
+               ij0 = 101   ;   ij1 = 101
+               ii0 = 139   ;   ii1 = 139   
+               saldta( mi0(ii0):mi1(ii1), mj0(ij0):mj1(ij1) , : , : ) =   &
+                  &                                    saldta( mi0(il0):mi1(il1) , mj0(ij0):mj1(ij1) , : , : )
+               ij0 = 101   ;   ij1 = 101
+               saldta( mi0(ii0):mi1(ii1), mj0(ij0):mj1(ij1) , : , : ) =   &
+                  &                                    saldta( mi0(il0):mi1(il1) , mj0(ij0):mj1(ij1) , : , : )
+               il0 = 138   ;   il1 = 138   
+               ij0 = 101   ;   ij1 = 102
+               ii0 = 139   ;   ii1 = 139   
+               DO jl = mi0(ii0), mi1(ii1)                ! New salinity profile at Gibraltar
+                  DO jj = mj0(ij0), mj1(ij1)
+                     DO ji = mi0(ii0), mi1(ii1)
                         saldta(ji,jj,:,:) = saldta(jl,jj,:,:)
                      END DO
@@ -217,12 +226,16 @@
                END DO
 
-               DO jl = mi0(164), mi1(164)                ! New salinity profile at Bab el Mandeb
-                  DO jj = mj0(88), mj1(88)
-                     DO ji = mi0(161), mi1(163)
+               il0 = 164   ;   il1 = 164   
+               ij0 =  88   ;   ij1 =  88
+               ii0 = 161   ;   ii1 = 163   
+               DO jl = mi0(ii0), mi1(ii1)                ! New salinity profile at Bab el Mandeb
+                  DO jj = mj0(ij0), mj1(ij1)
+                     DO ji = mi0(ii0), mi1(ii1)
                         saldta(ji,jj,:,:) = saldta(jl,jj,:,:)
                      END DO
                   END DO
-                  DO jj = mj0(87), mj1(87)
-                     DO ji = mi0(161), mi1(163)
+                  ij0 =  87   ;   ij1 =  87
+                  DO jj = mj0(ij0), mj1(ij1)
+                     DO ji = mi0(ii0), mi1(ii1)
                         saldta(ji,jj,:,:) = saldta(jl,jj,:,:)
                      END DO
@@ -280,8 +293,8 @@
    !!   Default option:                                    NO salinity data
    !!----------------------------------------------------------------------
-   LOGICAL, PUBLIC, PARAMETER ::   lk_dtasal = .FALSE.   ! salinity data flag
+   LOGICAL , PUBLIC, PARAMETER ::   lk_dtasal = .FALSE.   !: salinity data flag
 CONTAINS
    SUBROUTINE dta_sal( kt )        ! Empty routine
-      WRITE(*,*) kt
+      WRITE(*,*) 'dta_sal: You should not have seen this print! error?', kt
    END SUBROUTINE dta_sal
 #endif
Index: trunk/NEMO/OPA_SRC/DTA/dtasst.F90
===================================================================
--- trunk/NEMO/OPA_SRC/DTA/dtasst.F90	(revision 15)
+++ trunk/NEMO/OPA_SRC/DTA/dtasst.F90	(revision 16)
@@ -26,12 +26,12 @@
    !! * Shared module variables
 #if defined key_dtasst
-   LOGICAL , PUBLIC, PARAMETER ::   lk_dtasst = .TRUE.   ! sst data flag
+   LOGICAL , PUBLIC, PARAMETER ::   lk_dtasst = .TRUE.   !: sst data flag
 #else
-   LOGICAL , PUBLIC, PARAMETER ::   lk_dtasst = .TRUE.   ! sst data flag
+   LOGICAL , PUBLIC, PARAMETER ::   lk_dtasst = .FALSE.  !: sst data flag
 #endif
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   &
-      sst             ! surface temperature
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,2) ::   &
-      rclice          ! climatological ice index (0/1) (2 months)
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   &  !:
+      sst             !: surface temperature
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,2) ::   &  !:
+      rclice          !: climatological ice index (0/1) (2 months)
    !!----------------------------------------------------------------------
    !!   OPA 9.0 , IPSL-LODYC  (2003)
Index: trunk/NEMO/OPA_SRC/DTA/dtatem.F90
===================================================================
--- trunk/NEMO/OPA_SRC/DTA/dtatem.F90	(revision 15)
+++ trunk/NEMO/OPA_SRC/DTA/dtatem.F90	(revision 16)
@@ -23,7 +23,7 @@
 
    !! * Shared module variables
-   LOGICAL , PUBLIC, PARAMETER ::   lk_dtatem = .TRUE.   ! temperature data flag
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   &
-      t_dta             ! temperature data at given time-step
+   LOGICAL , PUBLIC, PARAMETER ::   lk_dtatem = .TRUE.   !: temperature data flag
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   &  !:
+      t_dta             !: temperature data at given time-step
 
    !! * Module variables
@@ -91,5 +91,6 @@
       INTEGER ::   &
          imois, iman, itime, ik ,    &  ! temporary integers
-         i15, ipi, ipj, ipk             !    "          "
+         i15, ipi, ipj, ipk,         &  !    "          "
+         il0, il1, ii0, ii1, ij0, ij1   !    "          "
 
       INTEGER, DIMENSION(jpmois) ::   istep
@@ -192,6 +193,8 @@
             !                                        ! ======================= 
 
-            DO jj = mj0(101), mj1(109)                   ! Reduced temperature at Alboran Sea
-               DO ji = mi0(141), mi1(155)
+            ij0 = 101   ;   ij1 = 109
+            ii0 = 141   ;   ii1 = 155
+            DO jj = mj0(ij0), mj1(ij1)                      ! Reduced temperature in the Alboran Sea
+               DO ji = mi0(ii0), mi1(ii1)
                   temdta(ji,jj, 13:13 ,:) = temdta(ji,jj, 13:13 ,:) - 0.20
                   temdta(ji,jj, 14:15 ,:) = temdta(ji,jj, 14:15 ,:) - 0.35
@@ -202,23 +205,32 @@
             IF( n_cla == 0 ) THEN 
                !                                         ! Reduced temperature at Red Sea
-               temdta( mi0(148):mi1(160) , mj0(87):mj1(96) ,  4:10 , : ) = 7.0 
-               temdta( mi0(148):mi1(160) , mj0(87):mj1(96) , 11:13 , : ) = 6.5 
-               temdta( mi0(148):mi1(160) , mj0(87):mj1(96) , 14:20 , : ) = 6.0
+               ij0 =  87   ;   ij1 =  96
+               ii0 = 148   ;   ii1 = 160
+               temdta( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ,  4:10 , : ) = 7.0 
+               temdta( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 11:13 , : ) = 6.5 
+               temdta( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 14:20 , : ) = 6.0
             ELSE
-               DO jl = mi0(138), mi1(138)                ! New temperature profile at Gibraltar
-                  DO jj = mj0(101), mj1(102)
-                     DO ji = mi0(139), mi1(139)
+               il0 = 138   ;   il1 = 138
+               ij0 = 101   ;   ij1 = 102
+               ii0 = 139   ;   ii1 = 139
+               DO jl = mi0(ii0), mi1(ii1)                ! New temperature profile at Gibraltar
+                  DO jj = mj0(ij0), mj1(ij1)
+                     DO ji = mi0(ii0), mi1(ii1)
                         temdta(ji,jj,:,:) = temdta(jl,jj,:,:)
                      END DO
                   END DO
                END DO
-               DO jl = mi0(164), mi1(164)                ! New temperature profile at Bab el Mandeb
-                  DO jj = mj0(88), mj1(88)
-                     DO ji = mi0(161), mi1(163)
+               il0 = 164   ;   il1 = 164
+               ij0 =  88   ;   ij1 =  88
+               ii0 = 161   ;   ii1 = 163
+               DO jl = mi0(ii0), mi1(ii1)                ! New temperature profile at Bab el Mandeb
+                  DO jj = mj0(ij0), mj1(ij1)
+                     DO ji = mi0(ii0), mi1(ii1)
                         temdta(ji,jj,:,:) = temdta(jl,jj,:,:)
                      END DO
                   END DO
-                  DO jj = mj0(87), mj1(87)
-                     DO ji = mi0(161), mi1(163)
+                  ij0 =  87   ;   ij1 =  87
+                  DO jj = mj0(ij0), mj1(ij1)
+                     DO ji = mi0(ii0), mi1(ii1)
                         temdta(ji,jj,:,:) = temdta(jl,jj,:,:)
                      END DO
@@ -274,8 +286,8 @@
    !!   Default case                           NO 3D temperature data field
    !!----------------------------------------------------------------------
-   LOGICAL , PUBLIC, PARAMETER ::   lk_dtatem = .FALSE.   ! temperature data flag
+   LOGICAL , PUBLIC, PARAMETER ::   lk_dtatem = .FALSE.   !: temperature data flag
 CONTAINS
    SUBROUTINE dta_tem( kt )        ! Empty routine
-      WRITE(*,*) kt
+      WRITE(*,*) 'dta_tem: You should not have seen this print! error?', kt
    END SUBROUTINE dta_tem
 #endif
Index: trunk/NEMO/OPA_SRC/FLO/flo4rk.F90
===================================================================
--- trunk/NEMO/OPA_SRC/FLO/flo4rk.F90	(revision 15)
+++ trunk/NEMO/OPA_SRC/FLO/flo4rk.F90	(revision 16)
@@ -15,4 +15,5 @@
    USE oce             ! ocean dynamics and tracers
    USE dom_oce         ! ocean space and time domain
+   USE in_out_manager  ! I/O manager
 
    IMPLICIT NONE
@@ -24,10 +25,10 @@
    !! * Module variables
    REAL(wp), DIMENSION (4) ::   &   ! RK4 and Lagrange interpolation
-      tcoef1 = /  1.0  ,  0.5  ,  0.5  ,  0.0  / ,  &  ! coeffients  for 
-      tcoef2 = /  0.0  ,  0.5  ,  0.5  ,  1.0  / ,  &  ! lagrangian interp.
-      scoef2 = /  1.0  ,  2.0  ,  2.0  ,  1.0  / ,  &  ! RK4 coefficients
-      rcoef  = /-1./6. , 1./2. ,-1./2. , 1./6. /       ! ???
+      tcoef1 = (/  1.0  ,  0.5  ,  0.5  ,  0.0  /) ,  &  ! coeffients  for 
+      tcoef2 = (/  0.0  ,  0.5  ,  0.5  ,  1.0  /) ,  &  ! lagrangian interp.
+      scoef2 = (/  1.0  ,  2.0  ,  2.0  ,  1.0  /) ,  &  ! RK4 coefficients
+      rcoef  = (/-1./6. , 1./2. ,-1./2. , 1./6. /)       ! ???
    REAL(wp), DIMENSION (3) ::   &
-      scoef1 = / .5, .5, 1. /       ! compute position with interpolated
+      scoef1 = (/ .5, .5, 1. /)       ! compute position with interpolated
    !!----------------------------------------------------------------------
    !!   OPA 9.0 , LODYC-IPSL  (2003)
Index: trunk/NEMO/OPA_SRC/FLO/flo_oce.F90
===================================================================
--- trunk/NEMO/OPA_SRC/FLO/flo_oce.F90	(revision 15)
+++ trunk/NEMO/OPA_SRC/FLO/flo_oce.F90	(revision 16)
@@ -11,5 +11,5 @@
    !!  OPA 9.0 , LODYC-IPSL (2003)
    !!----------------------------------------------------------------------
-#if defined key_floats
+#if   defined key_floats   ||   defined key_esopa
    !!----------------------------------------------------------------------
    !!   'key_floats'                                        drifting floats
@@ -20,12 +20,12 @@
    IMPLICIT NONE
 
-   LOGICAL, PUBLIC, PARAMETER ::   lk_floats = .TRUE.    ! float flag
+   LOGICAL, PUBLIC, PARAMETER ::   lk_floats = .TRUE.    !: float flag
 
    !! float parameters
    !! ----------------
    INTEGER, PARAMETER ::   &
-      jpnfl      = 23 ,          &  ! total number of floats during the run
-      jpnnewfl   =  0 ,          &  ! number of floats added in a new run
-      jpnrstarfl = jpnfl-jpnnewfl   ! number of floats for the restart
+      jpnfl     = 23 ,            &  ! total number of floats during the run
+      jpnnewflo =  0 ,            &  ! number of floats added in a new run
+      jpnrstflo = jpnfl-jpnnewflo    ! number of floats for the restart
 
    !! float variables
@@ -61,5 +61,5 @@
    !!   Default option :                                 NO drifting floats
    !!----------------------------------------------------------------------
-   LOGICAL, PUBLIC, PARAMETER ::   lk_floats = .FALSE.   ! float flag
+   LOGICAL, PUBLIC, PARAMETER ::   lk_floats = .FALSE.   !: float flag
 #endif
 
Index: trunk/NEMO/OPA_SRC/FLO/floats.F90
===================================================================
--- trunk/NEMO/OPA_SRC/FLO/floats.F90	(revision 15)
+++ trunk/NEMO/OPA_SRC/FLO/floats.F90	(revision 16)
@@ -67,7 +67,6 @@
       ENDIF
 
-# if defined key_mpp
-      CALL mppsync
-# endif
+      IF( lk_mpp )   CALL mppsync   ! synchronization of all the processor
+
 
       ! Writing and restart      
@@ -98,10 +97,10 @@
 
       !! * Local declarations
-      NAMELIST/namflo/ ln_rstarfl, nwritefl, nstockfl 
+      NAMELIST/namflo/ ln_rstflo, nwritefl, nstockfl 
       !!---------------------------------------------------------------------
       ! Namelist namflo : floats
       
       ! default values
-      ln_rstarfl  = .FALSE.
+      ln_rstflo  = .FALSE.
       nwritefl  = 150
       nstockfl  = 450
@@ -114,7 +113,7 @@
          WRITE(numout,*) ' '
          WRITE(numout,*) '         Namelist floats :'
-         WRITE(numout,*) '            restart                          ln_rstarfl = ', ln_rstarfl
-         WRITE(numout,*) '            frequency of float output file   nwritefl   = ', nwritefl
-         WRITE(numout,*) '            frequency of float restart file  nstockfl   = ', nstockfl
+         WRITE(numout,*) '            restart                          ln_rstflo = ', ln_rstflo
+         WRITE(numout,*) '            frequency of float output file   nwritefl  = ', nwritefl
+         WRITE(numout,*) '            frequency of float restart file  nstockfl  = ', nstockfl
          WRITE(numout,*) ' '
       ENDIF
@@ -128,5 +127,5 @@
 CONTAINS
    SUBROUTINE flo_stp( kt )          ! Empty routine
-      WRITE(*,*) kt
+      WRITE(*,*) 'flo_stp: You should not have seen this print! error?', kt
    END SUBROUTINE flo_stp
 #endif
Index: trunk/NEMO/OPA_SRC/FLO/floblk.F90
===================================================================
--- trunk/NEMO/OPA_SRC/FLO/floblk.F90	(revision 15)
+++ trunk/NEMO/OPA_SRC/FLO/floblk.F90	(revision 16)
@@ -17,4 +17,5 @@
    USE dom_oce         ! ocean space and time domain
    USE phycst          ! physical constants
+   USE in_out_manager  ! I/O manager
    USE lib_mpp         ! distribued memory computing library
 
@@ -105,5 +106,5 @@
       iloop = 0
 222   DO jfl = 1, jpnfl
-# if defined key_mpp
+# if   defined key_mpp_mpi   ||   defined key_mpp_shmem
          IF( (iil(jfl) >= (mig(nldi)-jpizoom+1)) .AND. (iil(jfl) <= (mig(nlei)-jpizoom+1)) .AND.   &
              (ijl(jfl) >= (mjg(nldj)-jpjzoom+1)) .AND. (ijl(jfl) <= (mjg(nlej)-jpjzoom+1)) ) THEN
@@ -320,5 +321,5 @@
             ! reinitialisation of the age of FLOAT
             zagefl(jfl) = zagenewfl(jfl)
-# if defined key_mpp
+# if   defined key_mpp_mpi   ||   defined key_mpp_shmem
          ELSE
             ! we put zgifl, zgjfl, zgkfl, zagefl
@@ -334,14 +335,10 @@
       
       ! synchronisation
-      ! sum of this arrays
-            
-# if defined key_mpp
-      CALL mpp_sum( zgifl , jpnfl )
-      CALL mpp_sum( zgjfl , jpnfl )
-      CALL mpp_sum( zgkfl , jpnfl )
-      CALL mpp_sum( zagefl, jpnfl )
-      CALL mpp_sum( iil   , jpnfl )
-      CALL mpp_sum( ijl   , jpnfl )
-# endif
+      IF( lk_mpp )   CALL mpp_sum( zgifl , jpnfl )   ! sums over the global domain
+      IF( lk_mpp )   CALL mpp_sum( zgjfl , jpnfl )
+      IF( lk_mpp )   CALL mpp_sum( zgkfl , jpnfl )
+      IF( lk_mpp )   CALL mpp_sum( zagefl, jpnfl )
+      IF( lk_mpp )   CALL mpp_sum( iil   , jpnfl )
+      IF( lk_mpp )   CALL mpp_sum( ijl   , jpnfl )
       
       ! in the case of open boundaries we need to test if the floats don't
Index: trunk/NEMO/OPA_SRC/FLO/flodom.F90
===================================================================
--- trunk/NEMO/OPA_SRC/FLO/flodom.F90	(revision 15)
+++ trunk/NEMO/OPA_SRC/FLO/flodom.F90	(revision 16)
@@ -13,7 +13,9 @@
    !!----------------------------------------------------------------------
    !! * Modules used
-   USE flo_oce         ! ocean drifting floats
    USE oce             ! ocean dynamics and tracers
    USE dom_oce         ! ocean space and time domain
+   USE flo_oce         ! ocean drifting floats
+!  USE floats
+   USE in_out_manager  ! I/O manager
    USE lib_mpp         ! distribued memory computing library
 
@@ -43,14 +45,14 @@
       !!----------------------------------------------------------------------      
       !! * Local declarations
-      LOGICAL   :: llinmesh
+      LOGICAL  :: llinmesh
       CHARACTER (len=21) ::  clname
-      INTEGER   :: ji, jj, jk               ! DO loop index on 3 directions
-      INTEGER   :: jfl, jfl1                ! number of floats   
-      INTEGER   :: inum = 11                ! logical unit for file read
-      INTEGER , DIMENSION ( jpnfl    )  ::   &
+      INTEGER  :: ji, jj, jk               ! DO loop index on 3 directions
+      INTEGER  :: jfl, jfl1                ! number of floats   
+      INTEGER  :: inum = 11                ! logical unit for file read
+      INTEGER, DIMENSION ( jpnfl    )  ::   &
          iimfl, ijmfl, ikmfl,    &          ! index mesh of floats
          idomfl,  ivtest, ihtest
-      REAL(wp)  :: zdxab,zdyad
-      REAL(wp) , DIMENSION ( jpnnewfl )  :: zgifl, zgjfl,  zgkfl
+      REAL(wp) :: zdxab, zdyad
+      REAL(wp), DIMENSION ( jpnnewflo+1 )  :: zgifl, zgjfl,  zgkfl
       !!---------------------------------------------------------------------
       
@@ -61,5 +63,5 @@
       IF(lwp) WRITE(numout,*) '           jpnfl = ',jpnfl
       
-      IF(ln_rstarfl) THEN
+      IF(ln_rstflo) THEN
          IF(lwp) WRITE(numout,*) '        float restart file read'
          
@@ -70,14 +72,14 @@
 
          ! read of the restart file
-         READ(inum) ( tpifl  (jfl), jfl=1, jpnrstarfl),   & 
-                        ( tpjfl  (jfl), jfl=1, jpnrstarfl),   &
-                        ( tpkfl  (jfl), jfl=1, jpnrstarfl),   &
-                        ( nisobfl(jfl), jfl=1, jpnrstarfl),   &
-                        ( ngrpfl (jfl), jfl=1, jpnrstarfl)    
+         READ(inum) ( tpifl  (jfl), jfl=1, jpnrstflo),   & 
+                        ( tpjfl  (jfl), jfl=1, jpnrstflo),   &
+                        ( tpkfl  (jfl), jfl=1, jpnrstflo),   &
+                        ( nisobfl(jfl), jfl=1, jpnrstflo),   &
+                        ( ngrpfl (jfl), jfl=1, jpnrstflo)    
          CLOSE(inum)
 
          ! if we want a  surface drift  ( like PROVOR floats )
          IF( ln_argo ) THEN
-            DO jfl = 1, jpnrstarfl
+            DO jfl = 1, jpnrstflo
                nisobfl(jfl) = 0
             END DO
@@ -87,10 +89,10 @@
          
          ! It is possible to add new floats.          
-         IF(lwp) WRITE(numout,*)' flo_dom:jpnfl jpnrstarfl ',jpnfl,jpnrstarfl
-         IF( jpnfl > jpnrstarfl ) THEN
+         IF(lwp) WRITE(numout,*)' flo_dom:jpnfl jpnrstflo ',jpnfl,jpnrstflo
+         IF( jpnfl > jpnrstflo ) THEN
             ! open the init file 
             clname='init_float'
             OPEN(inum,FILE=clname,FORM='FORMATTED')
-            DO jfl = jpnrstarfl+1, jpnfl
+            DO jfl = jpnrstflo+1, jpnfl
                READ(inum,*) flxx(jfl),flyy(jfl),flzz(jfl), nisobfl(jfl),ngrpfl(jfl),jfl1
             END DO
@@ -99,9 +101,9 @@
             
             ! Test to find the grid point coordonate with the geographical position            
-            DO jfl = jpnrstarfl+1, jpnfl
+            DO jfl = jpnrstflo+1, jpnfl
                ihtest(jfl) = 0
                ivtest(jfl) = 0
                ikmfl(jfl) = 0
-# if defined key_mpp
+# if   defined key_mpp_mpi   ||   defined key_mpp_shmem
                DO ji = MAX(nldi,2), nlei
                   DO jj = MAX(nldj,2), nlej
@@ -140,9 +142,9 @@
             
             ! A zero in the sum of the arrays "ihtest" and "ivtest"             
-# if defined key_mpp
-            CALL mpp_sum(ihtest,jpnfl,iwork)
-            CALL mpp_sum(ivtest,jpnfl,iwork)
+# if   defined key_mpp_mpi   ||   defined key_mpp_shmem
+            CALL mpp_sum(ihtest,jpnfl)
+            CALL mpp_sum(ivtest,jpnfl)
 # endif 
-            DO jfl = jpnrstarfl+1, jpnfl
+            DO jfl = jpnrstflo+1, jpnfl
                IF( (ihtest(jfl) > 1 ) .OR. ( ivtest(jfl) > 1) ) THEN
                   IF(lwp) WRITE(numout,*) 'THE FLOAT',jfl,' IS NOT IN ONLY ONE MESH'
@@ -156,5 +158,5 @@
             
             ! We compute the distance between the float and the face of the mesh            
-            DO jfl = jpnrstarfl+1, jpnfl               
+            DO jfl = jpnrstflo+1, jpnfl               
                ! Made only if the float is in the domain of the processor              
                IF( (iimfl(jfl) >= 0) .AND. (ijmfl(jfl) >= 0) ) THEN
@@ -182,31 +184,30 @@
                   ! Translation of this distances (in meter) in indexes
                   
-                  zgifl(jfl-jpnrstarfl)= (iimfl(jfl)-0.5) + zdxab/e1u(iimfl(jfl)-1,ijmfl(jfl)) + (mig(1)-jpizoom)
-                  zgjfl(jfl-jpnrstarfl)= (ijmfl(jfl)-0.5) + zdyad/e2v(iimfl(jfl),ijmfl(jfl)-1) + (mjg(1)-jpjzoom)
-                  zgkfl(jfl-jpnrstarfl) = (( fsdepw(ji,jj,ikmfl(jfl)+1) - flzz(jfl) )* ikmfl(jfl))  &
+                  zgifl(jfl-jpnrstflo)= (iimfl(jfl)-0.5) + zdxab/e1u(iimfl(jfl)-1,ijmfl(jfl)) + (mig(1)-jpizoom)
+                  zgjfl(jfl-jpnrstflo)= (ijmfl(jfl)-0.5) + zdyad/e2v(iimfl(jfl),ijmfl(jfl)-1) + (mjg(1)-jpjzoom)
+                  zgkfl(jfl-jpnrstflo) = (( fsdepw(ji,jj,ikmfl(jfl)+1) - flzz(jfl) )* ikmfl(jfl))  &
                                         / (  fsdepw(ji,jj,ikmfl(jfl)+1) - fsdepw(ji,jj,ikmfl(jfl) ) )   &
                                         + (( flzz(jfl)-fsdepw(ji,jj,ikmfl(jfl)) ) *(ikmfl(jfl)+1))   &
                                         / (  fsdepw(ji,jj,ikmfl(jfl)+1) - fsdepw(ji,jj,ikmfl(jfl)) ) 
                ELSE
-                  zgifl(jfl-jpnrstarfl) = 0.
-                  zgjfl(jfl-jpnrstarfl) = 0.
-                  zgkfl(jfl-jpnrstarfl) = 0.
+                  zgifl(jfl-jpnrstflo) = 0.
+                  zgjfl(jfl-jpnrstflo) = 0.
+                  zgkfl(jfl-jpnrstflo) = 0.
                ENDIF
             END DO
             
             ! The sum of all the arrays zgifl, zgjfl, zgkfl give 3 arrays with the positions of all the floats.
-# if defined key_mpp
-
-            CALL mpp_sum( zgjfl, jpnnewfl )
-            CALL mpp_sum( zgkfl, jpnnewfl )
-            IF(lwp) WRITE(numout,*) (zgifl(jfl),jfl=1,jpnnewfl)
-            IF(lwp) WRITE(numout,*) (zgjfl(jfl),jfl=1,jpnnewfl)
-            IF(lwp) WRITE(numout,*) (zgkfl(jfl),jfl=1,jpnnewfl) 
-# endif
+            IF( lk_mpp )   THEN
+               CALL mpp_sum( zgjfl, jpnnewflo )   ! sums over the global domain
+               CALL mpp_sum( zgkfl, jpnnewflo )
+               IF(lwp) WRITE(numout,*) (zgifl(jfl),jfl=1,jpnnewflo)
+               IF(lwp) WRITE(numout,*) (zgjfl(jfl),jfl=1,jpnnewflo)
+               IF(lwp) WRITE(numout,*) (zgkfl(jfl),jfl=1,jpnnewflo) 
+            ENDIF
            
-            DO jfl = jpnrstarfl+1, jpnfl
-               tpifl(jfl) = zgifl(jfl-jpnrstarfl)
-               tpjfl(jfl) = zgjfl(jfl-jpnrstarfl)
-               tpkfl(jfl) = zgkfl(jfl-jpnrstarfl)
+            DO jfl = jpnrstflo+1, jpnfl
+               tpifl(jfl) = zgifl(jfl-jpnrstflo)
+               tpjfl(jfl) = zgjfl(jfl-jpnrstflo)
+               tpkfl(jfl) = zgkfl(jfl-jpnrstflo)
             END DO
          ENDIF
@@ -234,5 +235,5 @@
             ivtest(jfl) = 0
             ikmfl(jfl) = 0
-# if defined key_mpp
+# if   defined key_mpp_mpi   ||   defined key_mpp_shmem
             DO ji = MAX(nldi,2), nlei
                DO jj = MAX(nldj,2), nlej
@@ -271,8 +272,6 @@
          
          ! A zero in the sum of the arrays "ihtest" and "ivtest"          
-# if defined key_mpp
-         CALL mpp_sum(ihtest,jpnfl,iwork)
-         CALL mpp_sum(ivtest,jpnfl,iwork)
-# endif
+         IF( lk_mpp )   CALL mpp_sum(ihtest,jpnfl)   ! sums over the global domain
+         IF( lk_mpp )   CALL mpp_sum(ivtest,jpnfl)
 
          DO jfl = 1, jpnfl
@@ -327,14 +326,12 @@
          
          ! The sum of all the arrays tpifl, tpjfl, tpkfl give 3 arrays with the positions of all the floats. 
-# if defined key_mpp
-         CALL mpp_sum( tpifl , jpnfl )
-         CALL mpp_sum( tpjfl , jpnfl )
-         CALL mpp_sum( tpkfl , jpnfl )
-         CALL mpp_sum( idomfl, jpnfl )
-# endif
+         IF( lk_mpp )   CALL mpp_sum( tpifl , jpnfl )   ! sums over the global domain
+         IF( lk_mpp )   CALL mpp_sum( tpjfl , jpnfl )
+         IF( lk_mpp )   CALL mpp_sum( tpkfl , jpnfl )
+         IF( lk_mpp )   CALL mpp_sum( idomfl, jpnfl )
       ENDIF
             
       ! Print the initial positions of the floats
-      IF( .NOT. ln_rstarfl ) THEN 
+      IF( .NOT. ln_rstflo ) THEN 
          ! WARNING : initial position not in the sea         
          DO jfl = 1, jpnfl
Index: trunk/NEMO/OPA_SRC/FLO/flowri.F90
===================================================================
--- trunk/NEMO/OPA_SRC/FLO/flowri.F90	(revision 15)
+++ trunk/NEMO/OPA_SRC/FLO/flowri.F90	(revision 16)
@@ -15,4 +15,5 @@
    USE dom_oce         ! ocean space and time domain
    USE lib_mpp         ! distribued memory computing library
+   USE daymod
    USE in_out_manager  ! I/O manager
 
@@ -52,4 +53,5 @@
 
       !! * Local declarations
+      CHARACTER (len=21) ::  clname
       INTEGER ::   inum = 11       ! temporary logical unit for restart file
       INTEGER  ::   &
@@ -62,7 +64,5 @@
       REAL(wp) :: zafl,zbfl,zcfl,zdtj
       REAL(wp) :: zxxu, zxxu_01,zxxu_10, zxxu_11
-      REAL(wp) , DIMENSION ( jpk  , jpnfl) :: ztemp, zsal
-
-      CHARACTER (len=21) ::  clname
+      REAL(wp), DIMENSION (jpk,jpnfl) :: ztemp, zsal
       !!---------------------------------------------------------------------
       
@@ -86,28 +86,75 @@
             IF(lwp) WRITE(numflo)cexper,no,irecflo,jpnfl,nwritefl
          ENDIF
-         zdtj = rdt/86400.      
+         zdtj = rdt / 86400.      !!bug   use of 86400 instead of the phycst parameter
 
          ! translation of index position in geographical position
 
-         DO jfl = 1, jpnfl
-            iafl  = INT (tpifl(jfl))
-            ibfl  = INT (tpjfl(jfl))
-            icfl  = INT (tpkfl(jfl))
-            iafln = NINT(tpifl(jfl))
-            ibfln = NINT(tpjfl(jfl))
-            ia1fl = iafl+1
-            ib1fl = ibfl+1
-            ic1fl = icfl+1
-            zafl  = tpifl(jfl) - FLOAT(iafl)
-            zbfl  = tpjfl(jfl) - FLOAT(ibfl)
-            zcfl  = tpkfl(jfl) - FLOAT(icfl)
-# if defined key_mpp
-            IF( (iafl >= (mig(nldi)-jpizoom+1)) .AND. (iafl <= (mig(nlei)-jpizoom+1)) .AND.   &
-              (  ibfl >= (mjg(nldj)-jpjzoom+1)) .AND. (ibfl <= (mjg(nlej)-jpjzoom+1)) ) THEN
-
-               ! local index
-
-               iafloc  = iafl -(mig(1)-jpizoom+1) + 1
-               ibfloc  = ibfl -(mjg(1)-jpjzoom+1) + 1
+         IF( lk_mpp ) THEN
+            DO jfl = 1, jpnfl
+               iafl  = INT ( tpifl(jfl) )
+               ibfl  = INT ( tpjfl(jfl) )
+               icfl  = INT ( tpkfl(jfl) )
+               iafln = NINT( tpifl(jfl) )
+               ibfln = NINT( tpjfl(jfl) )
+               ia1fl = iafl + 1
+               ib1fl = ibfl + 1
+               ic1fl = icfl + 1
+               zafl  = tpifl(jfl) - FLOAT( iafl )
+               zbfl  = tpjfl(jfl) - FLOAT( ibfl )
+               zcfl  = tpkfl(jfl) - FLOAT( icfl )
+               IF(   iafl >= mig(nldi)-jpizoom+1 .AND. iafl <= mig(nlei)-jpizoom+1 .AND.   &
+                  &  ibfl >= mjg(nldj)-jpjzoom+1 .AND. ibfl <= mjg(nlej)-jpjzoom+1       ) THEN
+
+                  ! local index
+
+                  iafloc  = iafl -(mig(1)-jpizoom+1) + 1
+                  ibfloc  = ibfl -(mjg(1)-jpjzoom+1) + 1
+                  ia1floc = iafloc + 1
+                  ib1floc = ibfloc + 1
+
+                  flyy(jfl) = (1.-zafl)*(1.-zbfl)*gphit(iafloc ,ibfloc ) + (1.-zafl) * zbfl * gphit(iafloc ,ib1floc)   &
+                     &      +     zafl *(1.-zbfl)*gphit(ia1floc,ibfloc ) +     zafl  * zbfl * gphit(ia1floc,ib1floc)
+                  flxx(jfl) = (1.-zafl)*(1.-zbfl)*glamt(iafloc ,ibfloc ) + (1.-zafl) * zbfl * glamt(iafloc ,ib1floc)   &
+                     &      +     zafl *(1.-zbfl)*glamt(ia1floc,ibfloc ) +     zafl  * zbfl * glamt(ia1floc,ib1floc)
+                  flzz(jfl) = (1.-zcfl)*fsdepw(iafloc,ibfloc,icfl ) + zcfl * fsdepw(iafloc,ibfloc,ic1fl)
+
+                  ! Change  by Alexandra Bozec et Jean-Philippe Boulanger
+                  ! We save  the instantaneous profile of T and S of the column     
+                  ! ztemp(jfl)=tn(iafloc,ibfloc,icfl)
+                  ! zsal(jfl)=sn(iafloc,ibfloc,icfl)
+                  ztemp(1:jpk,jfl) = tn(iafloc,ibfloc,1:jpk)
+                  zsal (1:jpk,jfl) = sn(iafloc,ibfloc,1:jpk)            
+               ELSE
+                  flxx(jfl) = 0.
+                  flyy(jfl) = 0.
+                  flzz(jfl) = 0.
+                  ztemp(1:jpk,jfl) = 0.
+                  zsal (1:jpk,jfl) = 0.
+               ENDIF
+            END DO
+
+            CALL mpp_sum( flxx, jpnfl )   ! sums over the global domain
+            CALL mpp_sum( flyy, jpnfl )
+            CALL mpp_sum( flzz, jpnfl )
+            ! these 2 lines have accendentaly been removed from ATL6-V8 run hence
+            ! giving 0 salinity and temperature on the float trajectory
+            CALL mpp_sum( ztemp, jpk*jpnfl )
+            CALL mpp_sum( zsal , jpk*jpnfl )
+
+         ELSE
+            DO jfl = 1, jpnfl
+               iafl  = INT (tpifl(jfl))
+               ibfl  = INT (tpjfl(jfl))
+               icfl  = INT (tpkfl(jfl))
+               iafln = NINT(tpifl(jfl))
+               ibfln = NINT(tpjfl(jfl))
+               ia1fl = iafl+1
+               ib1fl = ibfl+1
+               ic1fl = icfl+1
+               zafl  = tpifl(jfl) - FLOAT(iafl)
+               zbfl  = tpjfl(jfl) - FLOAT(ibfl)
+               zcfl  = tpkfl(jfl) - FLOAT(icfl)
+               iafloc  = iafl
+               ibfloc  = ibfl
                ia1floc = iafloc + 1
                ib1floc = ibfloc + 1
@@ -118,60 +165,25 @@
                          +     zafl *(1.-zbfl)*glamt(ia1floc,ibfloc ) +     zafl  * zbfl * glamt(ia1floc,ib1floc)
                flzz(jfl) = (1.-zcfl)*fsdepw(iafloc,ibfloc,icfl ) + zcfl * fsdepw(iafloc,ibfloc,ic1fl)
-
+               !ALEX
+               ! Astuce pour ne pas avoir des flotteurs qui se baladent sur IDL
+               zxxu_11 = glamt(iafloc ,ibfloc )
+               zxxu_10 = glamt(iafloc ,ib1floc)
+               zxxu_01 = glamt(ia1floc,ibfloc )
+               zxxu    = glamt(ia1floc,ib1floc)
+
+               IF( iafloc == 52 )  zxxu_10 = -181
+               IF( iafloc == 52 )  zxxu_11 = -181
+               flxx(jfl)=(1.-zafl)*(1.-zbfl)* zxxu_11 + (1.-zafl)*    zbfl * zxxu_10   &
+                        +    zafl *(1.-zbfl)* zxxu_01 +     zafl *    zbfl * zxxu
+               !ALEX         
                ! Change  by Alexandra Bozec et Jean-Philippe Boulanger
                ! We save  the instantaneous profile of T and S of the column     
-               ! ztemp(jfl)=tn(iafloc,ibfloc,icfl)
-               ! zsal(jfl)=sn(iafloc,ibfloc,icfl)
+               !     ztemp(jfl)=tn(iafloc,ibfloc,icfl)
+               !     zsal(jfl)=sn(iafloc,ibfloc,icfl)
                ztemp(1:jpk,jfl) = tn(iafloc,ibfloc,1:jpk)
-               zsal (1:jpk,jfl) = sn(iafloc,ibfloc,1:jpk)            
-            ELSE
-               flxx(jfl) = 0.
-               flyy(jfl) = 0.
-               flzz(jfl) = 0.
-               ztemp(1:jpk,jfl) = 0.
-               zsal (1:jpk,jfl) = 0.
-            ENDIF
-# else
-            iafloc  = iafl
-            ibfloc  = ibfl
-            ia1floc = iafloc + 1
-            ib1floc = ibfloc + 1
-            !
-            flyy(jfl) = (1.-zafl)*(1.-zbfl)*gphit(iafloc ,ibfloc ) + (1.-zafl) * zbfl * gphit(iafloc ,ib1floc)   &
-                      +     zafl *(1.-zbfl)*gphit(ia1floc,ibfloc ) +     zafl  * zbfl * gphit(ia1floc,ib1floc)
-            flxx(jfl) = (1.-zafl)*(1.-zbfl)*glamt(iafloc ,ibfloc ) + (1.-zafl) * zbfl * glamt(iafloc ,ib1floc)   &
-                      +     zafl *(1.-zbfl)*glamt(ia1floc,ibfloc ) +     zafl  * zbfl * glamt(ia1floc,ib1floc)
-            flzz(jfl) = (1.-zcfl)*fsdepw(iafloc,ibfloc,icfl ) + zcfl * fsdepw(iafloc,ibfloc,ic1fl)
-            !ALEX
-            ! Astuce pour ne pas avoir des flotteurs qui se baladent sur IDL
-            zxxu_11 = glamt(iafloc ,ibfloc )
-            zxxu_10 = glamt(iafloc ,ib1floc)
-            zxxu_01 = glamt(ia1floc,ibfloc )
-            zxxu    = glamt(ia1floc,ib1floc)
-
-            IF( iafloc == 52 )  zxxu_10 = -181
-            IF( iafloc == 52 )  zxxu_11 = -181
-            flxx(jfl)=(1.-zafl)*(1.-zbfl)* zxxu_11 + (1.-zafl)*    zbfl * zxxu_10   &
-                     +    zafl *(1.-zbfl)* zxxu_01 +     zafl *    zbfl * zxxu
-            !ALEX         
-            ! Change  by Alexandra Bozec et Jean-Philippe Boulanger
-            ! We save  the instantaneous profile of T and S of the column     
-            !     ztemp(jfl)=tn(iafloc,ibfloc,icfl)
-            !     zsal(jfl)=sn(iafloc,ibfloc,icfl)
-            ztemp(1:jpk,jfl) = tn(iafloc,ibfloc,1:jpk)
-            zsal (1:jpk,jfl) = sn(iafloc,ibfloc,1:jpk)
-# endif
-         END DO
-
-# if defined key_mpp
-         CALL mpp_sum( flxx, jpnfl )
-         CALL mpp_sum( flyy, jpnfl )
-         CALL mpp_sum( flzz, jpnfl )
-         ! these 2 lines have accendentaly been removed from ATL6-V8 run hence
-         ! giving 0 salinity and temperature on the float trajectory
-         CALL mpp_sum( ztemp, jpk*jpnfl )
-         CALL mpp_sum( zsal , jpk*jpnfl )
-
-# endif
+               zsal (1:jpk,jfl) = sn(iafloc,ibfloc,1:jpk)
+            END DO
+         ENDIF
+
          !
          WRITE(numflo) flxx,flyy,flzz,nisobfl,ngrpfl,ztemp,zsal, FLOAT(ndastp)
@@ -187,5 +199,5 @@
       !         iafln=NINT(tpifl(jfl))
       !         ibfln=NINT(tpjfl(jfl))
-      !# if defined key_mpp
+      !# if defined key_mpp_mpi   ||   defined key_mpp_shmem
       !        IF ( (iafl >= (mig(nldi)-jpizoom+1)) .AND.
       !     $       (iafl <= (mig(nlei)-jpizoom+1)) .AND.
@@ -206,5 +218,5 @@
       !         ztemp(jfl)=tn(iafloc,ibfloc,jk)
       !         zsal(jfl)=sn(iaflo!,ibfloc,jk)
-      !# if defined key_mpp
+      !# if defined key_mpp_mpi   ||   defined key_mpp_shmem
       !        ELSE
       !         ztemp(jfl) = 0.
@@ -214,8 +226,7 @@
       !! ... next float
       !        END DO
-      !#if defined key_mpp
-      !      CALL mpp_sum( ztemp, jpnfl )
-      !      CALL mpp_sum( zsal , jpnfl )
-      !# endif
+      !      IF( lk_mpp )   CALL mpp_sum( ztemp, jpnfl )
+      !      IF( lk_mpp )   CALL mpp_sum( zsal , jpnfl )
+      !
       !      IF (lwp) THEN 
       !         WRITE(numflo) ztemp, zsal
@@ -268,24 +279,24 @@
          ! Compute the number of trajectories for each processor
          !
-# if defined key_mpp
-         DO jfl = 1, jpnfl
-            IF( (INT(tpifl(jfl)) >= (mig(nldi)-jpizoom+1)) .AND.   &
-              (  INT(tpifl(jfl)) <= (mig(nlei)-jpizoom+1)) .AND.   &
-              (  INT(tpjfl(jfl)) >= (mjg(nldj)-jpjzoom+1)) .AND.   &
-              (  INT(tpjfl(jfl)) <= (mjg(nlej)-jpjzoom+1)) ) THEN
-               iproc(narea) = iproc(narea)+1
-            ENDIF
-         END DO
-         CALL mpp_sum( iproc, jpnij )
-         !
-         IF(lwp) THEN 
-            WRITE(numout,*) 'DATE',adatrj
-            DO jpn = 1, jpnij
-               IF( iproc(jpn) /= 0 ) THEN
-                  WRITE(numout,*)'PROCESSOR',jpn-1,'compute',iproc(jpn), 'trajectories.'
+         IF( lk_mpp ) THEN
+            DO jfl = 1, jpnfl
+               IF( (INT(tpifl(jfl)) >= (mig(nldi)-jpizoom+1)) .AND.   &
+                  &(INT(tpifl(jfl)) <= (mig(nlei)-jpizoom+1)) .AND.   &
+                  &(INT(tpjfl(jfl)) >= (mjg(nldj)-jpjzoom+1)) .AND.   &
+                  &(INT(tpjfl(jfl)) <= (mjg(nlej)-jpjzoom+1)) ) THEN
+                  iproc(narea) = iproc(narea)+1
                ENDIF
             END DO
-         ENDIF
-# endif
+            CALL mpp_sum( iproc, jpnij )
+            !
+            IF(lwp) THEN 
+               WRITE(numout,*) 'DATE',adatrj
+               DO jpn = 1, jpnij
+                  IF( iproc(jpn) /= 0 ) THEN
+                     WRITE(numout,*)'PROCESSOR',jpn-1,'compute',iproc(jpn), 'trajectories.'
+                  ENDIF
+               END DO
+            ENDIF
+         ENDIF
       ENDIF 
 
Index: trunk/NEMO/OPA_SRC/SOL/sol_oce.F90
===================================================================
--- trunk/NEMO/OPA_SRC/SOL/sol_oce.F90	(revision 15)
+++ trunk/NEMO/OPA_SRC/SOL/sol_oce.F90	(revision 16)
@@ -16,43 +16,48 @@
 
    IMPLICIT NONE
+   PRIVATE
 
-   !!----------------------------------------------------------------------
+   !!-----------------------------------
    !! elliptic solver: SOR, PCG or FETI
-   !! ---------------------------------------------------
-   INTEGER  ::              & !!! namsol   elliptic solver / island / free surface
-      nsolv =    1 ,        &  ! = 1/2/3 type of elliptic solver
-      nmax  =  800 ,        &  ! maximum of iterations for the solver
-      nmisl = 4000             ! maximum pcg iterations for island
+   !! ----------------------------------
+   INTEGER , PUBLIC ::      & !!: namsol   elliptic solver / island / free surface
+      nsolv =    1 ,        &  !: = 1/2/3 type of elliptic solver
+      nmax  =  800 ,        &  !: maximum of iterations for the solver
+      nmisl = 4000             !: maximum pcg iterations for island
      
-   REAL(wp) ::              & !!! namsol   elliptic solver / island / free surface
-      eps    = 1.e-6_wp ,   &  ! absolute precision of the solver
-      sor    = 1.76_wp  ,   &  ! optimal coefficient for sor solver
-      epsisl = 1.e-10_wp,   &  ! absolute precision on stream function solver
-      rnu    = 1.0_wp          ! strength of the additional force used in free surface
+   REAL(wp), PUBLIC ::      & !!: namsol   elliptic solver / island / free surface
+      eps    = 1.e-6_wp ,   &  !: absolute precision of the solver
+      sor    = 1.76_wp  ,   &  !: optimal coefficient for sor solver
+      epsisl = 1.e-10_wp,   &  !: absolute precision on stream function solver
+      rnu    = 1.0_wp          !: strength of the additional force used in free surface
 
-   INTEGER  ::   &
-      ncut,         &  ! indicator of solver convergence
-      niter            ! number of iteration done by the solver
+   CHARACTER(len=1), PUBLIC ::   &  !:
+      c_solver_pt = 'T'        !: nature of grid-points T (S) for free surface case
+      !                        !                       F (G) for rigid-lid case
 
-   REAL(wp) ::   &
-      epsr,         &  ! relative precision for SOR & PCG solvers
-      epsilo,       &  ! precision for the FETI solver
-      rnorme, res,  &  ! intermediate modulus, solver residu
-      alph,         &  ! coefficient  =(gcr,gcr)/(gcx,gccd)
-      beta,         &  ! coefficient  =(rn+1,rn+1)/(rn,rn)
-      radd,         &  ! coefficient  =(gccd,gcdes)
-      rr               ! coefficient  =(rn,rn)
+   INTEGER , PUBLIC ::   &  !:
+      ncut,         &  !: indicator of solver convergence
+      niter            !: number of iteration done by the solver
 
-   REAL(wp), DIMENSION(jpi,jpj,4) ::   &
-      gcp              ! barotropic matrix extra-diagonal elements
+   REAL(wp), PUBLIC ::   &  !:
+      epsr,         &  !: relative precision for SOR & PCG solvers
+      epsilo,       &  !: precision for the FETI solver
+      rnorme, res,  &  !: intermediate modulus, solver residu
+      alph,         &  !: coefficient  =(gcr,gcr)/(gcx,gccd)
+      beta,         &  !: coefficient  =(rn+1,rn+1)/(rn,rn)
+      radd,         &  !: coefficient  =(gccd,gcdes)
+      rr               !: coefficient  =(rn,rn)
 
-   REAL(wp), DIMENSION(jpi,jpj) ::   &
-      gcx, gcxb,    &  ! now, before solution of the elliptic equation
-      gcdprc,       &  ! inverse diagonal preconditioning matrix
-      gcdmat,       &  ! diagonal preconditioning matrix
-      gcb,          &  ! second member of the barotropic linear system
-      gcr,          &  ! residu =b-a.x
-      gcdes,        &  ! vector descente
-      gccd             ! vector such that ca.gccd=a.d (ca-1=gcdprc)
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,4) ::   &  !:
+      gcp              !: barotropic matrix extra-diagonal elements
+
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   &  !:
+      gcx, gcxb,    &  !: now, before solution of the elliptic equation
+      gcdprc,       &  !: inverse diagonal preconditioning matrix
+      gcdmat,       &  !: diagonal preconditioning matrix
+      gcb,          &  !: second member of the barotropic linear system
+      gcr,          &  !: residu =b-a.x
+      gcdes,        &  !: vector descente
+      gccd             !: vector such that ca.gccd=a.d (ca-1=gcdprc)
 
 #if defined key_feti
@@ -67,5 +72,5 @@
    !!      malistin()       : concatened list of interface nodes
 
-   INTEGER :: nim,nxm,   &
+   INTEGER, PUBLIC :: nim,nxm,   &
        malxm,malim,malxmax,malimax,   &
        nifmat,njfmat,nelem,npe,matopo,   &
@@ -85,7 +90,8 @@
        madwork
 
-   INTEGER :: mfet(jpi*jpj+2*jpi+2*jpj+51)
+   INTEGER, PUBLIC :: mfet(jpi*jpj+2*jpi+2*jpj+51)
 
-   REAL(wp) ::  wfeti(jpj*jpi*jpi+13*jpi*jpj+19*(jpi+jpj)   &
+   REAL(wp), PUBLIC ::   &  !:
+      wfeti(jpj*jpi*jpi+13*jpi*jpj+19*(jpi+jpj)   &
        +4*jpnij+33   &
        +2*(jpi+jpj)*(jpnij-jpni)*jpi   &
@@ -94,5 +100,5 @@
        +3*(jpnij-jpnj+jperio)*jpj) 
 
-   REAL(wp) ::   res2, rcompt
+   REAL(wp), PUBLIC ::   res2, rcompt
 
 #endif
Index: trunk/NEMO/OPA_SRC/SOL/solfet.F90
===================================================================
--- trunk/NEMO/OPA_SRC/SOL/solfet.F90	(revision 15)
+++ trunk/NEMO/OPA_SRC/SOL/solfet.F90	(revision 16)
@@ -33,5 +33,5 @@
       !!     Solve the ellipic equation for the barotropic stream function
       !!     system (default option) or the transport divergence system
-      !!     ("key_dynspg_fsc") using a Finite Elements Tearing & 
+      !!     (lk_dynspg_fsc=T) using a Finite Elements Tearing and 
       !!      Interconnecting (FETI) approach.
       !!     In the former case, the barotropic stream function trend has a
@@ -142,26 +142,5 @@
       CALL feti_vmov( noeuds, wfeti(miax), gcx )
 
-      ! boundary conditions   !!bug ???  check  arguments...
-#   if defined key_dynspg_fsc
-#      if defined key_mpp
-      !   Mpp: export boundary values to neighbouring processors
-      CALL lbc_lnk( gcx, 'S', 1. )
-#      else
-      !   mono- or macro-tasking: W-point, >0, 2D array, no slab
-      IF( nperio /= 0 ) THEN
-         CALL lbc_lnk( gcx, 'T', 1. )
-      ENDIF
-#      endif
-#   else
-#      if defined key_mpp
-      !   Mpp: export boundary values to neighbouring processors
-      CALL lbc_lnk( gcx, 'G', 1. )
-#      else
-      !   mono- or macro-tasking: W-point, >0, 2D array, no slab
-      IF( nperio /= 0 ) THEN
-         CALL lbc_lnk( gcx, 'F', 1. )
-      ENDIF
-#      endif
-#   endif
+      CALL lbc_lnk( gcx, c_solver_pt, 1. )   ! lateral boundary condition
 
    END SUBROUTINE sol_fet
Index: trunk/NEMO/OPA_SRC/SOL/solisl.F90
===================================================================
--- trunk/NEMO/OPA_SRC/SOL/solisl.F90	(revision 15)
+++ trunk/NEMO/OPA_SRC/SOL/solisl.F90	(revision 16)
@@ -36,6 +36,5 @@
 
    !! * Shared module variables
-   LOGICAL, PUBLIC ::   &
-      l_isl = .TRUE.          ! 'key_islands' flag
+   LOGICAL, PUBLIC, PARAMETER ::   l_isl = .TRUE.    !: 'key_islands' flag
 
    !! * module variable
@@ -157,58 +156,55 @@
          zwb(jpi,:) = 0.e0
       ENDIF
-# if defined key_mpp
-      ! Mpp: export boundary values to neighboring processors
-      CALL lbc_lnk( zwb, 'G', 1. )
-# endif
+      IF( lk_mpp )   CALL lbc_lnk( zwb, 'G', 1. )
 
 
       ! 1. Initialization for the search of island grid-points
       ! ------------------------------------------------------
-# if defined key_mpp
-
-      ! Mpp : The overlap region are not taken into account
-      ! (islands bondaries are searched over subdomain only)
-      iista =  1   + jpreci
-      iiend = nlci - jpreci
-      ijsta =  1   + jprecj
-      ijend = nlcj - jprecj
-      ijstm1=  1   + jprecj
-      ijenm1= nlcj - jprecj
-      IF( nbondi == -1 .OR. nbondi == 2 ) THEN
+
+      IF( lk_mpp ) THEN
+
+         ! Mpp : The overlap region are not taken into account
+         ! (islands bondaries are searched over subdomain only)
+         iista =  1   + jpreci
+         iiend = nlci - jpreci
+         ijsta =  1   + jprecj
+         ijend = nlcj - jprecj
+         ijstm1=  1   + jprecj
+         ijenm1= nlcj - jprecj
+         IF( nbondi == -1 .OR. nbondi == 2 ) THEN
+            iista  = 1
+         ENDIF
+         IF( nbondi ==  1 .OR. nbondi == 2 ) THEN
+            iiend  = nlci
+         ENDIF
+         IF( nbondj == -1 .OR. nbondj == 2 ) THEN
+            ijsta  = 1
+            ijstm1 = 2
+         ENDIF
+         IF( nbondj ==  1 .OR. nbondj == 2 ) THEN
+            ijend  = nlcj
+            ijenm1 = nlcj-1
+         ENDIF
+         IF( npolj == 3 .OR. npolj == 4 ) THEN
+            ijend  = nlcj-2
+            ijenm1 = nlcj-2
+         ENDIF 
+      ELSE
+         ! mono- or macro-tasking environnement: full domain scan
          iista  = 1
-      ENDIF
-      IF( nbondi ==  1 .OR. nbondi == 2 ) THEN
-         iiend  = nlci
-      ENDIF
-      IF( nbondj == -1 .OR. nbondj == 2 ) THEN
+         iiend  = jpi
          ijsta  = 1
          ijstm1 = 2
-      ENDIF
-      IF( nbondj ==  1 .OR. nbondj == 2 ) THEN
-         ijend  = nlcj
-         ijenm1 = nlcj-1
-      ENDIF
-      IF( npolj == 3 .OR. npolj == 4 ) THEN
-         ijend  = nlcj-2
-         ijenm1 = nlcj-2
-      ENDIF 
-# else
-
-      ! mono- or macro-tasking environnement: full domain scan
-      iista  = 1
-      iiend  = jpi
-      ijsta  = 1
-      ijstm1 = 2
-      IF( nperio == 3 .OR. nperio == 4 ) THEN
-         ijend  = jpj-2
-         ijenm1 = jpj-2
-      ELSEIF( nperio == 5 .OR. nperio == 6 ) THEN
-         ijend  = jpj-1
-         ijenm1 = jpj-1
-      ELSE
-         ijend  = jpj
-         ijenm1 = jpj-1
-      ENDIF
-# endif
+         IF( nperio == 3 .OR. nperio == 4 ) THEN
+            ijend  = jpj-2
+            ijenm1 = jpj-2
+         ELSEIF( nperio == 5 .OR. nperio == 6 ) THEN
+            ijend  = jpj-1
+            ijenm1 = jpj-1
+         ELSE
+            ijend  = jpj
+            ijenm1 = jpj-1
+         ENDIF
+      ENDIF
 
 
@@ -247,7 +243,6 @@
             inilt = inilt + indil(jj)
          END DO
-# if defined key_mpp
-         CALL mpp_sum( inilt )
-# endif
+         IF( lk_mpp )   CALL mpp_sum( inilt )   ! sum over the global domain
+
          IF( inilt == 0 ) THEN
             IF(lwp) THEN
@@ -255,5 +250,5 @@
                WRITE(numout,*) ' change parameter.h'
             ENDIF
-            STOP 'isldom'
+            STOP 'isldom'      !cr replace by nstop
          ENDIF
          
@@ -381,7 +376,5 @@
          ! Take account of redundant points
          
-# if defined key_mpp
-         CALL mpp_sum( ip )
-# endif
+         IF( lk_mpp )   CALL mpp_sum( ip )   ! sum over the global domain
          
          IF( ip > jpnisl ) THEN
@@ -391,5 +384,5 @@
                WRITE(numout,*) ' change parameter.h'
             ENDIF
-            STOP 'isldom'
+            STOP 'isldom'    !cr => nstop
          ENDIF
          
@@ -409,7 +402,5 @@
 
       inilt = isrchne( jpij, zwb(1,1), 1, 0. )
-# if defined key_mpp
-      CALL mpp_min( inilt )
-# endif
+      IF( lk_mpp )   CALL mpp_min( inilt )   ! min over the global domain
 
       IF( inilt /= jpij+1 ) THEN
@@ -426,5 +417,5 @@
       ! ----------------------------------------
 
-      CALL islpri
+      CALL isl_pri
 
 
@@ -432,5 +423,5 @@
       ! -------------------------------------------------------
 
-      CALL islpth
+      CALL isl_pth
 
    END SUBROUTINE isl_dom
@@ -466,11 +457,11 @@
          ipe = mnisl(3,jni)
          ipw = mnisl(4,jni)
-# if defined key_mpp
-         CALL mpp_sum( ip  )
-         CALL mpp_sum( ipn )
-         CALL mpp_sum( ips )
-         CALL mpp_sum( ipe )
-         CALL mpp_sum( ipw )
-# endif
+         IF( lk_mpp ) THEN
+            CALL mpp_sum( ip  )   ! sums over the global domain
+            CALL mpp_sum( ipn )
+            CALL mpp_sum( ips )
+            CALL mpp_sum( ipe )
+            CALL mpp_sum( ipw )
+         ENDIF
          IF(lwp) THEN
             WRITE(numout,9000) jni
@@ -484,6 +475,5 @@
       END DO
 
-      ! FORMAT 
-
+      ! FORMAT   !!cr => no more format
  9000 FORMAT(/, /, 'island number= ', i2 )
  9010 FORMAT(/, 'npil=',i4,' npn=',i3,' nps=',i3,' npe=',i3,' npw=',i3 )
@@ -514,6 +504,6 @@
       !!----------------------------------------------------------------------
       !! * Local declarations
-      INTEGER ::   ji, jj, jni, jii, jnp    ! dummy loop indices
-      INTEGER ::   ii, ij                   ! temporary integers
+      INTEGER ::   jni, jii, jnp    ! dummy loop indices
+      INTEGER ::   ii, ij           ! temporary integers
       !!----------------------------------------------------------------------
       
@@ -587,7 +577,5 @@
       REAL(wp), DIMENSION(jpi,jpj) ::   zlamt, zphit
       REAL(wp), DIMENSION(jpi,jpj,2) ::   zwx
-# if defined key_mpp
       REAL(wp), DIMENSION(jpisl*jpisl) ::   ztab
-# endif
       !!----------------------------------------------------------------------
 
@@ -674,14 +662,13 @@
          
       END DO
-# if defined key_mpp
-      DO jnj=1,jpisl
-         DO jni=1,jpisl
-            ztab(jni+(jnj-1)*jpisl) = aisl(jni,jnj)
-         END DO
-      END DO
-
-      CALL mpp_sum( ztab, jpisl*jpisl )
-!!    CALL mpp_sum( aisl, jpisl*jpisl )
-# endif
+      IF( lk_mpp ) THEN
+         DO jnj = 1, jpisl
+            DO jni = 1, jpisl
+               ztab(jni+(jnj-1)*jpisl) = aisl(jni,jnj)
+            END DO
+         END DO
+         CALL mpp_sum( ztab, jpisl*jpisl )   ! sum over the global domain
+!!       CALL mpp_sum( aisl, jpisl*jpisl )
+      ENDIF
 
       ! 1.3 Control print
@@ -775,7 +762,6 @@
       REAL(wp) ::   zep(jpisl), zlamt(jpi,jpj), zphit(jpi,jpj), zdept(1), zprec(4)
       REAL(wp) ::   zdate0, zdt
-# if defined key_mpp
+      REAL(wp) ::   t2p1(jpi,1,1)
       INTEGER  ::   iloc
-# endif
       !!----------------------------------------------------------------------
 
@@ -907,79 +893,78 @@
          ! Right hand side of the streamfunction equation
          
-# if defined key_mpp
-
-         ! north fold treatment
-         IF( npolj == 3 .OR. npolj == 5)   iloc=jpiglo-(nimpp-1+nimppt(nono+1)-1)
-         IF( npolj == 4 .OR. npolj == 6)   iloc=jpiglo-2*(nimpp-1)
-         t2p1(:,1,1) = 0.e0
-         ! north and south grid-points
-         DO jii = 1, 2
-            DO jnp = 1, mnisl(jii,jni)
-               ii = miisl(jnp,jii,jni)
-               ij = mjisl(jnp,jii,jni)
-               IF( ( npolj == 3 .OR. npolj == 4 ) .AND.   &
-                  ( ij == nlcj-1 .AND. jii == 1) ) THEN 
-                  iju=iloc-ii+1
-                  t2p1(iju,1,1) =  t2p1(iju,1,1) + hur(ii,ij) * e1u(ii,ij) / e2u(ii,ij) 
-               ELSEIF( ( npolj == 5 .OR. npolj == 6 ) .AND.   &
-                  ( ij == nlcj-1 .AND. jii == 1) ) THEN 
-                  iju=iloc-ii
-                  gcb(ii,ij) =  gcb(ii,ij) + hur(ii,ij) * e1u(ii,ij) / e2u(ii,ij) 
-                  t2p1(iju,1,1) =  t2p1(iju,1,1) + hur(ii,ij) * e1u(ii,ij) / e2u(ii,ij) 
-               ELSE  
-                  gcb(ii,ij-jii+1) =  gcb(ii,ij-jii+1) + hur(ii,ij) * e1u(ii,ij) / e2u(ii,ij) 
-               ENDIF
-            END DO
-         END DO
-         
-         ! east and west grid-points
-         
-         DO jii = 3, 4
-            DO jnp = 1, mnisl(jii,jni)
-               ii = miisl(jnp,jii,jni)
-               ij = mjisl(jnp,jii,jni)
-               gcb(ii-jii+3,ij) = gcb(ii-jii+3,ij) + hvr(ii,ij) * e2v(ii,ij) / e1v(ii,ij)
-            END DO
-         END DO
-         CALL mpplnks( gcb )
-
-# else
-
-         ! north and south grid-points
-         DO jii = 1, 2
-            DO jnp = 1, mnisl(jii,jni)
-               ii = miisl(jnp,jii,jni)
-               ij = mjisl(jnp,jii,jni)
-               IF( ( nperio == 3 .OR. nperio == 4 ) .AND.   &
-                  ( ij == jpj-1 .AND. jii == 1) ) THEN 
-                  gcb(jpi-ii+1,ij-1) = gcb(jpi-ii+1,ij-1) + hur(ii,ij) * e1u(ii,ij) / e2u(ii,ij) 
-               ELSEIF( ( nperio == 5 .OR. nperio == 6 ) .AND.   &
-                  ( ij == jpj-1 .AND. jii == 1) ) THEN 
-                  gcb(ii,ij) =  gcb(ii,ij) + hur(ii,ij) * e1u(ii,ij) / e2u(ii,ij)
-                  gcb(jpi-ii,ij) = gcb(jpi-ii,ij) + hur(ii,ij) * e1u(ii,ij) / e2u(ii,ij) 
-               ELSE  
-                  gcb(ii,ij-jii+1) =  gcb(ii,ij-jii+1) + hur(ii,ij) * e1u(ii,ij) / e2u(ii,ij)
-               ENDIF
-            END DO
-         END DO
-
-         ! east and west grid-points
-         DO jii = 3, 4
-            DO jnp = 1, mnisl(jii,jni)
-               ii = miisl(jnp,jii,jni)
-               ij = mjisl(jnp,jii,jni)
-               IF( bmask(ii-jii+3,ij) /= 0. ) THEN
+         IF( lk_mpp ) THEN
+
+            ! north fold treatment
+            IF( npolj == 3 .OR. npolj == 5)   iloc=jpiglo-(nimpp-1+nimppt(nono+1)-1)
+            IF( npolj == 4 .OR. npolj == 6)   iloc=jpiglo-2*(nimpp-1)
+            t2p1(:,1,1) = 0.e0
+            ! north and south grid-points
+            DO jii = 1, 2
+               DO jnp = 1, mnisl(jii,jni)
+                  ii = miisl(jnp,jii,jni)
+                  ij = mjisl(jnp,jii,jni)
+                  IF( ( npolj == 3 .OR. npolj == 4 ) .AND.   &
+                     ( ij == nlcj-1 .AND. jii == 1) ) THEN 
+                     iju=iloc-ii+1
+                     t2p1(iju,1,1) =  t2p1(iju,1,1) + hur(ii,ij) * e1u(ii,ij) / e2u(ii,ij) 
+                  ELSEIF( ( npolj == 5 .OR. npolj == 6 ) .AND.   &
+                     ( ij == nlcj-1 .AND. jii == 1) ) THEN 
+                     iju=iloc-ii
+                     gcb(ii,ij) =  gcb(ii,ij) + hur(ii,ij) * e1u(ii,ij) / e2u(ii,ij) 
+                     t2p1(iju,1,1) =  t2p1(iju,1,1) + hur(ii,ij) * e1u(ii,ij) / e2u(ii,ij) 
+                  ELSE  
+                     gcb(ii,ij-jii+1) =  gcb(ii,ij-jii+1) + hur(ii,ij) * e1u(ii,ij) / e2u(ii,ij) 
+                  ENDIF
+               END DO
+            END DO
+         
+            ! east and west grid-points
+         
+            DO jii = 3, 4
+               DO jnp = 1, mnisl(jii,jni)
+                  ii = miisl(jnp,jii,jni)
+                  ij = mjisl(jnp,jii,jni)
                   gcb(ii-jii+3,ij) = gcb(ii-jii+3,ij) + hvr(ii,ij) * e2v(ii,ij) / e1v(ii,ij)
-               ELSE
-
-                  ! east-west cyclic boundary conditions
-                  IF( ii-jii+3 == 1 ) THEN
-                     gcb(jpim1,ij) = gcb(jpim1,ij) + hvr(ii,ij) * e2v(ii,ij) / e1v(ii,ij)
+               END DO
+            END DO
+
+            IF( lk_mpp )   CALL mpplnks( gcb )   !!bug ? should use an lbclnk ? is it possible?
+
+         ELSE
+
+            ! north and south grid-points
+            DO jii = 1, 2
+               DO jnp = 1, mnisl(jii,jni)
+                  ii = miisl(jnp,jii,jni)
+                  ij = mjisl(jnp,jii,jni)
+                  IF( ( nperio == 3 .OR. nperio == 4 ) .AND.   &
+                     ( ij == jpj-1 .AND. jii == 1) ) THEN 
+                     gcb(jpi-ii+1,ij-1) = gcb(jpi-ii+1,ij-1) + hur(ii,ij) * e1u(ii,ij) / e2u(ii,ij) 
+                  ELSEIF( ( nperio == 5 .OR. nperio == 6 ) .AND.   &
+                     ( ij == jpj-1 .AND. jii == 1) ) THEN 
+                     gcb(ii,ij) =  gcb(ii,ij) + hur(ii,ij) * e1u(ii,ij) / e2u(ii,ij)
+                     gcb(jpi-ii,ij) = gcb(jpi-ii,ij) + hur(ii,ij) * e1u(ii,ij) / e2u(ii,ij) 
+                  ELSE  
+                     gcb(ii,ij-jii+1) =  gcb(ii,ij-jii+1) + hur(ii,ij) * e1u(ii,ij) / e2u(ii,ij)
                   ENDIF
-               ENDIF
-            END DO
-         END DO
-
-# endif
+               END DO
+            END DO
+
+            ! east and west grid-points
+            DO jii = 3, 4
+               DO jnp = 1, mnisl(jii,jni)
+                  ii = miisl(jnp,jii,jni)
+                  ij = mjisl(jnp,jii,jni)
+                  IF( bmask(ii-jii+3,ij) /= 0. ) THEN
+                     gcb(ii-jii+3,ij) = gcb(ii-jii+3,ij) + hvr(ii,ij) * e2v(ii,ij) / e1v(ii,ij)
+                  ELSE
+                     ! east-west cyclic boundary conditions
+                     IF( ii-jii+3 == 1 ) THEN
+                        gcb(jpim1,ij) = gcb(jpim1,ij) + hvr(ii,ij) * e2v(ii,ij) / e1v(ii,ij)
+                     ENDIF
+                  ENDIF
+               END DO
+            END DO
+         ENDIF
 
          ! Preconditioned right hand side and absolute precision
@@ -1011,7 +996,6 @@
                END DO
             END DO
-# if defined key_mpp
-            CALL mpp_sum( rnorme )
-# endif
+            IF( lk_mpp )   CALL mpp_sum( rnorme )   ! sum over the global domain
+
             IF(lwp) WRITE(numout,*) 'rnorme ', rnorme
             epsr = epsisl * epsisl * rnorme
@@ -1070,7 +1054,5 @@
             END DO
          ENDIF
-# if defined key_mpp
-         CALL lbc_lnk( bsfisl(:,:,jni), 'G', 1. )
-# endif
+         IF( lk_mpp )   CALL lbc_lnk( bsfisl(:,:,jni), 'G', 1. )   ! link at G-point
 
 
@@ -1212,8 +1194,6 @@
          END DO
       END DO
-#   if defined key_mpp
-      !  Mpp : global sum to obtain global dot from local ones
-      CALL mpp_sum( bisl, jpisl )
-#   endif
+      IF( lk_mpp )   CALL mpp_sum( bisl, jpisl )   ! sum over the global domain
+
       DO jni = 1, jpisl                     ! Island stream function trend
          visl(jni) = 0.e0
@@ -1270,7 +1250,6 @@
             zfact = 1.e-6 * bsfn(miisl(1,0,jni),mjisl(1,0,jni))
          ENDIF
-#      if defined key_mpp
-         CALL mpp_isl( zfact )
-#      endif
+         IF( lk_mpp )   CALL mpp_isl( zfact )
+
          IF(lwp) WRITE(numisp,9300) kt, jni, zfact, visl(jni)
          IF( MOD( kt, nwrite ) == 0 .OR. kindic < 0     &
@@ -1293,5 +1272,5 @@
    !!   Default option                                         Empty module
    !!----------------------------------------------------------------------
-   LOGICAL, PUBLIC ::   l_isl = .FALSE.    ! 'key_islands' flag
+   LOGICAL, PUBLIC, PARAMETER ::   l_isl = .FALSE.    !: 'key_islands' flag
 CONTAINS
    SUBROUTINE isl_dom                        ! Empty routine
@@ -1304,5 +1283,5 @@
    END SUBROUTINE isl_dyn_spg
    SUBROUTINE isl_stp_ctl( kt, kindic )      ! Empty routine
-      WRITE(*,*) kt, kindic                     ! no compilation warning
+      WRITE(*,*) 'isl_stp_ctl: You should not have seen this print! error?', kt, kindic
    END SUBROUTINE isl_stp_ctl
 #endif
Index: trunk/NEMO/OPA_SRC/SOL/solisl_fdir.h90
===================================================================
--- trunk/NEMO/OPA_SRC/SOL/solisl_fdir.h90	(revision 15)
+++ trunk/NEMO/OPA_SRC/SOL/solisl_fdir.h90	(revision 16)
@@ -109,10 +109,7 @@
 
       END DO
-#   if defined key_mpp
-      CALL mpp_sum( aisl, jpisl*jpisl )
-#   endif
+      IF( lk_mpp )   CALL mpp_sum( aisl, jpisl*jpisl )   ! sum over the global domain
 
       ! 1.3 Control print
-      
       IF(lwp) THEN
          WRITE(numout,*)
@@ -296,67 +293,65 @@
          ! 1.2 Right hand side of the stream FUNCTION equation
          
-#    if defined key_mpp
-         
-         ! north fold treatment
-         IF( npolj == 3 ) iloc = jpiglo -(nimpp-1+nimppt(nono+1)-1)
-         IF( npolj == 4 ) iloc = jpiglo - 2*(nimpp-1)
-         t2p1(:,1,1) = 0.
-         ! north and south grid-points
-         DO jii = 1, 2
-            DO jnp = 1, mnisl(jii,jni)
-               ii = miisl(jnp,jii,jni)
-               ij = mjisl(jnp,jii,jni)
-               IF( ( npolj == 3 .OR. npolj == 4 ) .AND. ( ij == nlcj-1 .AND. jii == 1) ) THEN 
-                  iju=iloc-ii+1
-                  t2p1(iju,1,1) =  t2p1(iju,1,1) + hur(ii,ij) * e1u(ii,ij) / e2u(ii,ij) 
-               ELSE  
-                  gcb(ii,ij-jii+1) =  gcb(ii,ij-jii+1) + hur(ii,ij) * e1u(ii,ij) / e2u(ii,ij) 
-               ENDIF
-            END DO
-         END DO
-         
-         ! east and west grid-points
-
-         DO jii = 3, 4
-            DO jnp = 1, mnisl(jii,jni)
-               ii = miisl(jnp,jii,jni)
-               ij = mjisl(jnp,jii,jni)
-               gcb(ii-jii+3,ij) = gcb(ii-jii+3,ij) + hvr(ii,ij) * e2v(ii,ij) / e1v(ii,ij)
-            END DO
-         END DO
-         CALL mpplnks( gcb )
-        
-#      else
-
-         ! north and south grid-points
-         DO jii = 1, 2
-            DO jnp = 1, mnisl(jii,jni)
-               ii = miisl(jnp,jii,jni)
-               ij = mjisl(jnp,jii,jni)
-               IF( ( nperio == 3 .OR. nperio == 4 ) .AND. ( ij == jpj-1 .AND. jii == 1) ) THEN 
-                  gcb(jpi-ii+1,ij-1) = gcb(jpi-ii+1,ij-1) + hur(ii,ij) * e1u(ii,ij) / e2u(ii,ij) 
-               ELSE  
-                  gcb(ii,ij-jii+1) =  gcb(ii,ij-jii+1) + hur(ii,ij) * e1u(ii,ij) / e2u(ii,ij)
-               ENDIF
-            END DO
-         END DO
-
-         ! east and west grid-points
-         DO jii = 3, 4
-            DO jnp = 1, mnisl(jii,jni)
-               ii = miisl(jnp,jii,jni)
-               ij = mjisl(jnp,jii,jni)
-               IF( bmask(ii-jii+3,ij) /= 0. ) THEN
+         IF( lk_mpp ) THEN
+            ! north fold treatment
+            IF( npolj == 3 ) iloc = jpiglo -(nimpp-1+nimppt(nono+1)-1)
+            IF( npolj == 4 ) iloc = jpiglo - 2*(nimpp-1)
+            t2p1(:,1,1) = 0.
+            ! north and south grid-points
+            DO jii = 1, 2
+               DO jnp = 1, mnisl(jii,jni)
+                  ii = miisl(jnp,jii,jni)
+                  ij = mjisl(jnp,jii,jni)
+                  IF( ( npolj == 3 .OR. npolj == 4 ) .AND. ( ij == nlcj-1 .AND. jii == 1) ) THEN 
+                     iju=iloc-ii+1
+                     t2p1(iju,1,1) =  t2p1(iju,1,1) + hur(ii,ij) * e1u(ii,ij) / e2u(ii,ij) 
+                  ELSE  
+                     gcb(ii,ij-jii+1) =  gcb(ii,ij-jii+1) + hur(ii,ij) * e1u(ii,ij) / e2u(ii,ij) 
+                  ENDIF
+               END DO
+            END DO
+         
+            ! east and west grid-points
+
+            DO jii = 3, 4
+               DO jnp = 1, mnisl(jii,jni)
+                  ii = miisl(jnp,jii,jni)
+                  ij = mjisl(jnp,jii,jni)
                   gcb(ii-jii+3,ij) = gcb(ii-jii+3,ij) + hvr(ii,ij) * e2v(ii,ij) / e1v(ii,ij)
-               ELSE
-                  ! east-west cyclic boundary conditions
-                  IF( ii-jii+3 == 1 ) THEN
-                     gcb(jpim1,ij) = gcb(jpim1,ij) + hvr(ii,ij) * e2v(ii,ij) / e1v(ii,ij)
+               END DO
+            END DO
+
+            IF( lk_mpp )   CALL mpplnks( gcb )   !!bug ? should use an lbclnk ? is it possible???
+
+         ELSE
+            ! north and south grid-points
+            DO jii = 1, 2
+               DO jnp = 1, mnisl(jii,jni)
+                  ii = miisl(jnp,jii,jni)
+                  ij = mjisl(jnp,jii,jni)
+                  IF( ( nperio == 3 .OR. nperio == 4 ) .AND. ( ij == jpj-1 .AND. jii == 1) ) THEN 
+                     gcb(jpi-ii+1,ij-1) = gcb(jpi-ii+1,ij-1) + hur(ii,ij) * e1u(ii,ij) / e2u(ii,ij) 
+                  ELSE  
+                     gcb(ii,ij-jii+1) =  gcb(ii,ij-jii+1) + hur(ii,ij) * e1u(ii,ij) / e2u(ii,ij)
                   ENDIF
-               ENDIF
-            END DO
-         END DO
-        
-#    endif
+               END DO
+            END DO
+
+            ! east and west grid-points
+            DO jii = 3, 4
+               DO jnp = 1, mnisl(jii,jni)
+                  ii = miisl(jnp,jii,jni)
+                  ij = mjisl(jnp,jii,jni)
+                  IF( bmask(ii-jii+3,ij) /= 0. ) THEN
+                     gcb(ii-jii+3,ij) = gcb(ii-jii+3,ij) + hvr(ii,ij) * e2v(ii,ij) / e1v(ii,ij)
+                  ELSE
+                     ! east-west cyclic boundary conditions
+                     IF( ii-jii+3 == 1 ) THEN
+                        gcb(jpim1,ij) = gcb(jpim1,ij) + hvr(ii,ij) * e2v(ii,ij) / e1v(ii,ij)
+                     ENDIF
+                  ENDIF
+               END DO
+            END DO
+         ENDIF
 
          ! 1.4 Preconditioned right hand side and absolute precision
@@ -388,7 +383,6 @@
                END DO
             END DO
-#if defined key_mpp
-            CALL mpp_sum( rnorme )
-#endif
+            IF( lk_mpp )   CALL mpp_sum( rnorme )   ! sum over the global domain
+
             IF(lwp) WRITE(numout,*) 'rnorme ', rnorme
             epsr = epsisl * epsisl * rnorme
@@ -451,7 +445,5 @@
             END DO
          ENDIF
-#if defined key_mpp
-         CALL lbc_lnk( bsfisl(:,:,jni), 'G', 1. )
-#endif
+         IF( lk_mpp )   CALL lbc_lnk( bsfisl(:,:,jni), 'G', 1. )   ! link at G-point
          
          
Index: trunk/NEMO/OPA_SRC/SOL/solmat.F90
===================================================================
--- trunk/NEMO/OPA_SRC/SOL/solmat.F90	(revision 15)
+++ trunk/NEMO/OPA_SRC/SOL/solmat.F90	(revision 16)
@@ -18,4 +18,6 @@
    USE obc_oce         ! ocean open boundary conditions
    USE lib_mpp         ! distributed memory computing
+   USE dynspg_rl
+   USE dynspg_fsc
 
    IMPLICIT NONE
@@ -38,8 +40,8 @@
       !!
       !! ** Method  :   The matrix depends on the type of free surface:
-      !!       * default option: rigid lid and bsf
+      !!       * lk_dynspg_rl=T: rigid lid formulation
       !!      The matrix is built for the barotropic stream function system.
       !!      a diagonal preconditioning matrix is also defined.
-      !!       * 'key_dynspg_fsc' defined: free surface
+      !!       * lk_dynspg_fsc=T: free surface formulation
       !!      The matrix is built for the divergence of the transport system
       !!      a diagonal preconditioning matrix is also defined.
@@ -67,8 +69,5 @@
       INTEGER ::   ii, ij, iiend, ijend      ! temporary integers
       REAL(wp) ::   zcoefs, zcoefw, zcoefe, zcoefn  ! temporary scalars
-      REAL(wp) ::   z2dt
-#if defined key_dynspg_fsc
-      REAL(wp) ::   zcoef
-#endif
+      REAL(wp) ::   z2dt, zcoef
       !!----------------------------------------------------------------------
 
@@ -83,4 +82,5 @@
       
       ! initialize to zero
+      zcoef = 0.e0
       gcp(:,:,1) = 0.e0
       gcp(:,:,2) = 0.e0
@@ -94,4 +94,5 @@
 
 #if defined key_dynspg_fsc && ! defined key_obc
+!!cr      IF( lk_dynspg_fsc .AND. .NOT.lk_obc ) THEN
 
       ! defined the coefficients for free surface elliptic system
@@ -99,5 +100,5 @@
       DO jj = 2, jpjm1
          DO ji = 2, jpim1
-            zcoef = z2dt * z2dt * g * rnu * bmask(ji,jj)
+            zcoef = z2dt * z2dt * grav * rnu * bmask(ji,jj)
             zcoefs = -zcoef * hv(ji  ,jj-1) * e1v(ji  ,jj-1) / e2v(ji  ,jj-1)    ! south coefficient
             zcoefw = -zcoef * hu(ji-1,jj  ) * e2u(ji-1,jj  ) / e1u(ji-1,jj  )    ! west coefficient
@@ -109,9 +110,10 @@
             gcp(ji,jj,4) = zcoefn
             gcdmat(ji,jj) = e1t(ji,jj) * e2t(ji,jj) * bmask(ji,jj)    &          ! diagonal coefficient
-                          - zcoefs -zcoefw -zcoefe -zcoefn
+               &          - zcoefs -zcoefw -zcoefe -zcoefn
          END DO
       END DO
       
 #  elif defined key_dynspg_fsc && defined key_obc
+!!cr      ELSEIF( lk_dynspg_fsc .AND. lk_obc ) THEN
 
       !   defined gcdmat in the case of open boundaries
@@ -119,5 +121,5 @@
       DO jj = 2, jpjm1
          DO ji = 2, jpim1
-            zcoef = z2dt * z2dt * g * rnu * bmask(ji,jj)
+            zcoef = z2dt * z2dt * grav * rnu * bmask(ji,jj)
             !  south coefficient
             IF( ( lpsouthobc ) .AND. ( jj == njs0p1 ) ) THEN
@@ -159,4 +161,5 @@
 
 #  else
+!!cr      ELSE
 
       !   defined the coefficients for bsf elliptic system
@@ -173,8 +176,8 @@
             gcp(ji,jj,4) = zcoefn
             gcdmat(ji,jj) = -zcoefs -zcoefw -zcoefe -zcoefn                             ! diagonal coefficient
-            
          END DO
       END DO
       
+!!cr  ENDIF
 #endif
 
@@ -194,4 +197,5 @@
       ! account for the existence of the south symmetric bassin.
       
+!!cr      IF( .NOT.lk_dynspg_fsc ) THEN
 #if ! defined key_dynspg_fsc
       IF( nperio == 2 ) THEN
@@ -203,4 +207,5 @@
          END DO
       ENDIF
+!!cr      ENDIF
 #endif
       
@@ -225,4 +230,5 @@
          gcp(:,:,3) = gcp(:,:,3) * gcdprc(:,:)
          gcp(:,:,4) = gcp(:,:,4) * gcdprc(:,:)
+         IF( nsolv == 2 )   gccd(:,:) = sor * gcp(:,:,2)
 
       ELSE
@@ -467,5 +473,5 @@
       nnitot = nni
 
-      CALL mpp_sum(nnitot,1,numit0ete)
+      CALL mpp_sum( nnitot, 1, numit0ete )
       CALL feti_creadr(malxm,malxmax,nxm,npe*npe,maae,'ae')
 
Index: trunk/NEMO/OPA_SRC/SOL/solpcg.F90
===================================================================
--- trunk/NEMO/OPA_SRC/SOL/solpcg.F90	(revision 15)
+++ trunk/NEMO/OPA_SRC/SOL/solpcg.F90	(revision 16)
@@ -14,4 +14,5 @@
    USE lib_mpp         ! distributed memory computing
    USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
+   USE in_out_manager  ! I/O manager
 
    IMPLICIT NONE
@@ -32,6 +33,6 @@
       !!                    
       !! ** Purpose :   Solve the ellipic equation for the barotropic stream 
-      !!      function system (default option) or the transport divergence 
-      !!      system ("key_dynspg_fsc") using a diagonal preconditionned
+      !!      function system (lk_dynspg_rl=T) or the transport divergence 
+      !!      system (lk_dynspg_fsc=T) using a diagonal preconditionned
       !!      conjugate gradient method.
       !!      In the former case, the barotropic stream function trend has a
@@ -93,42 +94,17 @@
          !                                             !================
 
-         !,,,,,,,,,,,,,,,,,,,,,,,,synchro,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
-         
-         IF( jn == 1 ) THEN
-            
-            ! 1.0 Initialization of the algorithm
-            ! -----------------------------------
-            
-#if defined key_dynspg_fsc
-#   if defined key_mpp
-            ! Mpp: export boundary values to neighbouring processors
-            CALL lbc_lnk( gcx, 'S', 1. )
-#   else
-            !   mono- or macro-tasking: W-point, >0, 2D array, no slab
-            CALL lbc_lnk( gcx, 'T', 1. )
-#   endif
-#else
-#   if defined key_mpp
-            ! ... Mpp: export boundary values to neighbouring processors
-            CALL lbc_lnk( gcx, 'G', 1. )
-#   else
-            !   ... mono- or macro-tasking: F-point, >0, 2D array, no slab
-            CALL lbc_lnk( gcx, 'F', 1. )
-#   endif
-#endif
+         IF( jn == 1 ) THEN           ! Initialization of the algorithm
 
-            !,,,,,,,,,,,,,,,,,,,,,,,,synchro ,,,,,,,,,,,,,,,,,,,,,,,
-            
+            CALL lbc_lnk( gcx, c_solver_pt, 1. )   ! lateral boundary condition
+
             ! gcr   = gcb-a.gcx
             ! gcdes = gsr
-            
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  zgcad = bmask(ji,jj)*(                         &
-                    gcb(ji,jj  ) -              gcx(ji  ,jj  )   &
-                                 - gcp(ji,jj,1)*gcx(ji  ,jj-1)   &
-                                 - gcp(ji,jj,2)*gcx(ji-1,jj  )   &
-                                 - gcp(ji,jj,3)*gcx(ji+1,jj  )   &
-                                 - gcp(ji,jj,4)*gcx(ji  ,jj+1)   )
+                  zgcad = bmask(ji,jj) * ( gcb(ji,jj  ) -                gcx(ji  ,jj  )   &
+                     &                                  - gcp(ji,jj,1) * gcx(ji  ,jj-1)   &
+                     &                                  - gcp(ji,jj,2) * gcx(ji-1,jj  )   &
+                     &                                  - gcp(ji,jj,3) * gcx(ji+1,jj  )   &
+                     &                                  - gcp(ji,jj,4) * gcx(ji  ,jj+1)   )
                   gcr  (ji,jj) = zgcad
                   gcdes(ji,jj) = zgcad
@@ -136,140 +112,74 @@
             END DO
             
-            !,,,,,,,,,,,,,,,,,,,,,,,,synchro ,,,,,,,,,,,,,,,,,,,,,,,
-            
             rnorme = SUM(  gcr(:,:) * gcdmat(:,:) * gcr(:,:)  )
-
-#if defined key_mpp
-            ! Mpp: sum over all the global domain
-            CALL mpp_sum( rnorme )
-#endif
+            IF( lk_mpp )   CALL mpp_sum( rnorme )   ! sum over the global domain
             rr = rnorme
 
-        ENDIF
-        !,,,,,,,,,,,,,,,,,,,,,,,,synchro ,,,,,,,,,,,,,,,,,,,,,,,
+         ENDIF
         
+         !                             ! Algorithm
         
-        ! 1.1 Algorithm
-        ! -------------
+         CALL lbc_lnk( gcdes, c_solver_pt, 1. )   ! lateral boundary condition
         
-        ! boundary condition on gcdes (only cyclic bc are required)
-#if defined key_dynspg_fsc
-#   if defined key_mpp
-        !   Mpp: export boundary values to neighbouring processors
-        CALL lbc_lnk( gcdes, 'S', 1. )
-#   else
-        !   mono- or macro-tasking: W-point, >0, 2D array, no slab
-        CALL lbc_lnk( gcdes, 'T', 1. )
-#   endif
-#else
-#   if defined key_mpp
-        !   Mpp: export boundary values to neighbouring processors
-        CALL lbc_lnk( gcdes, 'G', 1. )
-#   else
-        !   mono- or macro-tasking: F-point, >0, 2D array, no slab
-        CALL lbc_lnk( gcdes, 'F', 1. )
-#   endif
-#endif
+         ! ... gccd = matrix . gcdes
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               gccd(ji,jj) = bmask(ji,jj)*( gcdes(ji,jj)   &
+                  &        +gcp(ji,jj,1)*gcdes(ji,jj-1)+gcp(ji,jj,2)*gcdes(ji-1,jj)   &
+                  &        +gcp(ji,jj,4)*gcdes(ji,jj+1)+gcp(ji,jj,3)*gcdes(ji+1,jj)   )
+            END DO
+         END DO
+ 
+         ! alph = (gcr,gcr)/(gcdes,gccd)
+         radd = SUM(  gcdes(:,:) * gcdmat(:,:) * gccd(:,:)  )
+         IF( lk_mpp )   CALL mpp_sum( radd )   ! sum over the global domain
+         alph = rr / radd
+         
+         ! gcx = gcx + alph * gcdes
+         ! gcr = gcr - alph * gccd
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               gcx(ji,jj) = bmask(ji,jj) * ( gcx(ji,jj) + alph * gcdes(ji,jj) )
+               gcr(ji,jj) = bmask(ji,jj) * ( gcr(ji,jj) - alph * gccd (ji,jj) )
+            END DO
+         END DO
         
-        !,,,,,,,,,,,,,,,,,,,,,,,,synchro if macrotasking,,,,,,,,,,,,,,,,,,,,,,,
+         ! rnorme = (gcr,gcr)
+         rnorme = SUM(  gcr(:,:) * gcdmat(:,:) * gcr(:,:)  )
+         IF( lk_mpp )   CALL  mpp_sum( rnorme )   ! sum over the global domain
         
-        ! ... gccd = matrix . gcdes
-        DO jj = 2, jpjm1
-           DO ji = fs_2, fs_jpim1   ! vector opt.
-              gccd(ji,jj) = bmask(ji,jj)*(   &
-                 gcdes(ji,jj)   &
-                +gcp(ji,jj,1)*gcdes(ji,jj-1)+gcp(ji,jj,2)*gcdes(ji-1,jj)   &
-                +gcp(ji,jj,4)*gcdes(ji,jj+1)+gcp(ji,jj,3)*gcdes(ji+1,jj)   &
-                )
-           END DO
-        END DO
+         ! test of convergence
+         IF( rnorme < epsr .OR. jn == nmax ) THEN
+            res = SQRT( rnorme )
+            niter = jn
+            ncut = 999
+         ENDIF
+        
+         ! beta = (rk+1,rk+1)/(rk,rk)
+         beta = rnorme / rr
+         rr   = rnorme
 
-        !,,,,,,,,,,,,,,,,,,,,,,,,synchro if macrotasking,,,,,,,,,,,,,,,,,,,,,,,
+         ! indicator of non-convergence or explosion
+         IF( jn == nmax .OR. SQRT(epsr)/eps > 1.e+20 ) kindic = -2
+         IF( ncut == 999 ) GOTO 999
+
+         ! gcdes = gcr + beta * gcdes
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               gcdes(ji,jj) = bmask(ji,jj)*( gcr(ji,jj) + beta * gcdes(ji,jj) )
+            END DO
+         END DO
         
-        ! alph = (gcr,gcr)/(gcdes,gccd)
-
-        radd = SUM(  gcdes(:,:) * gcdmat(:,:) * gccd(:,:)  )
-
-#if defined key_mpp
-        ! Mpp: sum over all the global domain
-        CALL mpp_sum( radd )
-#endif
-        alph = rr / radd
-        
-        !,,,,,,,,,,,,,,,,,,,,,,,,synchro if macrotasking,,,,,,,,,,,,,,,,,,,,,,,
-        
-        ! gcx = gcx + alph * gcdes
-        ! gcr = gcr - alph * gccd
-        DO jj = 2, jpjm1
-           DO ji = fs_2, fs_jpim1   ! vector opt.
-              gcx(ji,jj) = bmask(ji,jj) * ( gcx(ji,jj) + alph * gcdes(ji,jj) )
-              gcr(ji,jj) = bmask(ji,jj) * ( gcr(ji,jj) - alph * gccd (ji,jj) )
-           END DO
-        END DO
-        
-        !,,,,,,,,,,,,,,,,,,,,,,,,synchro if macrotasking,,,,,,,,,,,,,,,,,,,,,,,
-        
-        ! rnorme = (gcr,gcr)
-
-        rnorme = SUM(  gcr(:,:) * gcdmat(:,:) * gcr(:,:)  )
-
-#if defined key_mpp
-        ! Mpp: sum over all the global domain
-        CALL  mpp_sum( rnorme )
-#endif
-        
-        ! test of convergence
-        IF ( rnorme < epsr .OR. jn == nmax ) THEN
-           res = SQRT( rnorme )
-           niter = jn
-           ncut = 999
-        ENDIF
-        
-        ! beta = (rk+1,rk+1)/(rk,rk)
-        beta = rnorme / rr
-        rr   = rnorme
-
-        !,,,,,,,,,,,,,,,,,,,,,,,,synchro if macrotasking,,,,,,,,,,,,,,,,,,,,,,,
-
-        ! indicator of non-convergence or explosion
-        IF( jn == nmax .OR. sqrt(epsr)/eps > 1.e+20 ) kindic = -2
-        IF( ncut == 999 ) GOTO 999
-
-        ! gcdes = gcr + beta * gcdes
-        DO jj = 2, jpjm1
-           DO ji = fs_2, fs_jpim1   ! vector opt.
-              gcdes(ji,jj) = bmask(ji,jj)*( gcr(ji,jj) + beta * gcdes(ji,jj) )
-           END DO
-        END DO
-        
-        !                                             !================
-     END DO                                           !    End Loop
-     !                                                !================
+         !                                             !================
+      END DO                                           !    End Loop
+      !                                                !================
      
-999  CONTINUE
+999   CONTINUE
      
      
-     !  2. Output in gcx with lateral b.c. applied
-     !  ------------------------------------------
+      ! Output in gcx with lateral b.c. applied
+      ! ---------------------------------------
      
-     ! boundary conditions   !!bug ???
-#if defined key_mpp
-     ! Mpp: export boundary values to neighbouring processors
-# if defined key_dynspg_fsc
-     CALL lbc_lnk( gcx, 'S', 1. )
-# else
-     CALL lbc_lnk( gcx, 'G', 1. )
-# endif
-#else
-     IF ( nperio /= 0 ) THEN
-# if defined key_dynspg_fsc
-        ! mono- or macro-tasking: W-point, >0, 2D array, no slab
-        CALL lbc_lnk( gcx, 'T', 1. )
-# else
-        ! mono- or macro-tasking: F-point, >0, 2D array, no slab
-        CALL lbc_lnk( gcx, 'F', 1. )
-# endif
-     ENDIF
-#endif
+      CALL lbc_lnk( gcx, c_solver_pt, 1. )
      
    END SUBROUTINE sol_pcg
Index: trunk/NEMO/OPA_SRC/SOL/solsor.F90
===================================================================
--- trunk/NEMO/OPA_SRC/SOL/solsor.F90	(revision 15)
+++ trunk/NEMO/OPA_SRC/SOL/solsor.F90	(revision 16)
@@ -22,4 +22,5 @@
    !! * Routine accessibility
    PUBLIC sol_sor              ! ???
+
    !!----------------------------------------------------------------------
    !!   OPA 9.0 , LODYC-IPSL  (2003)
@@ -28,11 +29,11 @@
 CONTAINS
       
-   SUBROUTINE sol_sor( kt, kindic )
+   SUBROUTINE sol_sor( kindic )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE sol_sor  ***
       !!                 
       !! ** Purpose :   Solve the ellipic equation for the barotropic stream 
-      !!      function system (default option) or the transport divergence 
-      !!      system (key_dynspg_fsc = T) using a successive-over-relaxation
+      !!      function system (lk_dynspg_rl=T) or the transport divergence 
+      !!      system (lk_dynspg_fsc=T) using a successive-over-relaxation
       !!      method.
       !!       In the former case, the barotropic stream function trend has a
@@ -59,5 +60,4 @@
       !!----------------------------------------------------------------------
       !! * Arguments
-      INTEGER, INTENT(  in   ) ::   kt       ! ocean time-step
       INTEGER, INTENT( inout ) ::   kindic   ! solver indicator, < 0 if the conver-
       !                                      ! gence is not reached: the model is
@@ -67,137 +67,46 @@
       !! * Local declarations
       INTEGER  ::   ji, jj, jn               ! dummy loop indices
-      REAL(wp) ::   zgwgt                    ! temporary scalar
       !!----------------------------------------------------------------------
       
-      
-      ! Iterative loop 
-      ! ==============
-      
-      IF( kt == nit000 )  gccd(:,:) = sor * gcp(:,:,2)
+      !                                                       ! ==============
+      DO jn = 1, nmax                                         ! Iterative loop 
+         !                                                    ! ==============
 
-
-      DO jn = 1, nmax
-
-         !,,,,,,,,,,,,,,,,,,,,,,,,,,,,,synchro,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
+         CALL lbc_lnk( gcx, c_solver_pt, 1. )   ! applied the lateral boubary conditions
          
-         ! boundary conditions (at each sor iteration) only cyclic b. c. are required
-#if defined key_dynspg_fsc
-# if defined key_mpp
-            ! Mpp: export boundary values to neighbouring processors
-            CALL lbc_lnk( gcx, 'S', 1. )   ! S=T with special staff ???
-# else
-            CALL lbc_lnk( gcx, 'T', 1. )
-# endif
-#else
-# if defined key_mpp
-            ! Mpp: export boundary values to neighbouring processors
-            CALL lbc_lnk( gcx, 'G', 1. )   ! G= F with special staff ???
-# else
-            CALL lbc_lnk( gcx, 'F', 1. )
-# endif
-#endif
-         
-         !,,,,,,,,,,,,,,,,,,,,,,,,,,,,,synchro,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
-         
-         ! 1. Residus
-         ! ----------
- 
+         ! Residus
+         ! -------
          DO jj = 2, jpjm1
             DO ji = 2, jpim1
-               gcr(ji,jj) =  gcb(ji,jj  ) -             gcx(ji  ,jj  )   &
-                                          -gcp(ji,jj,1)*gcx(ji  ,jj-1)   &
-                                          -gcp(ji,jj,2)*gcx(ji-1,jj  )   &
-                                          -gcp(ji,jj,3)*gcx(ji+1,jj  )   &
-                                          -gcp(ji,jj,4)*gcx(ji  ,jj+1)
+               gcr(ji,jj) =  gcb(ji,jj  ) -                gcx(ji  ,jj  )   &
+                                          - gcp(ji,jj,1) * gcx(ji  ,jj-1)   &
+                                          - gcp(ji,jj,2) * gcx(ji-1,jj  )   &
+                                          - gcp(ji,jj,3) * gcx(ji+1,jj  )   &
+                                          - gcp(ji,jj,4) * gcx(ji  ,jj+1)
             END DO
          END DO
+         CALL lbc_lnk( gcr, c_solver_pt, 1. )   ! applied the lateral boubary conditions
 
-         !,,,,,,,,,,,,,,,,,,,,,,,,,,,,,synchro,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
-         
-         ! 1.1 Boundary conditions (at each sor iteration) only cyclic b. c. are required
-#if defined key_dynspg_fsc
-#   if defined key_mpp
-            ! Mpp: export boundary values to neighbouring processors
-            CALL lbc_lnk( gcr, 'S', 1. )
-#   else
-            ! mono- or macro-tasking: W-point, >0, 2D array, no slab
-            CALL lbc_lnk( gcr, 'T', 1. )
-#   endif
-#else
-#   if defined key_mpp
-            ! Mpp: export boundary values to neighbouring processors
-            CALL lbc_lnk( gcr, 'G', 1. )
-#   else
-            ! mono- or macro-tasking: W-point, >0, 2D array, no slab
-            CALL lbc_lnk( gcr, 'F', 1. )
-#   endif
-#endif
-
-         ! 1.2 Successive over relaxation
-         
+         ! Successive over relaxation
          DO jj = 2, jpj
             DO ji = 1, jpi
-               gcr(ji,jj) = gcr(ji,jj) - sor*gcp(ji,jj,1)*gcr(ji,jj-1)
+               gcr(ji,jj) = gcr(ji,jj) - sor * gcp(ji,jj,1) * gcr(ji,jj-1)
             END DO
             DO ji = 2, jpi
-               gcr(ji,jj) = gcr(ji,jj) - sor*gcp(ji,jj,2)*gcr(ji-1,jj)
+               gcr(ji,jj) = gcr(ji,jj) - sor * gcp(ji,jj,2) * gcr(ji-1,jj)
             END DO
          END DO
          
-         !,,,,,,,,,,,,,,,,,,,,,,,,,,,,,synchro,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
-         
          ! gcx guess
-         
          DO jj = 2, jpjm1
             DO ji = 1, jpi
-               gcx(ji,jj)  = (gcx(ji,jj)+sor*gcr(ji,jj))*bmask(ji,jj)
+               gcx(ji,jj)  = ( gcx(ji,jj) + sor * gcr(ji,jj) ) * bmask(ji,jj)
             END DO
          END DO
-         
-         !,,,,,,,,,,,,,,,,,,,,,,,,,,,,,synchro,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
-         
-         ! boundary conditions (at each sor iteration) only cyclic b. c. are required
-#if defined key_dynspg_fsc
-#   if defined key_mpp
-         ! Mpp: export boundary values to neighbouring processors
-         CALL lbc_lnk( gcx, 'S', 1. )
-#   else
-         ! mono- or macro-tasking: W-point, >0, 2D array, no slab
-         CALL lbc_lnk( gcx, 'T', 1. )
-#   endif
-#else
-#   if defined key_mpp
-         ! Mpp: export boundary values to neighbouring processors
-         CALL lbc_lnk( gcx, 'G', 1. )
-#   else
-         ! mono- or macro-tasking: W-point, >0, 2D array, no slab
-         CALL lbc_lnk( gcx, 'F', 1. )
-#   endif
-#endif
-         
-         ! maximal residu (old exit test on the maximum value of residus)
-         ! 
-         ! imax = isamax( jpi*jpj, gcr, 1 )
-         
-         ! avoid an out of bound in no bounds compilation
-         
-         ! iimax1 = mod( imax, jpi )
-         ! ijmax1 = int( float(imax) / float(jpi)) + 1
-         ! resmax = abs( gcr(iimax1,ijmax1) )
+         CALL lbc_lnk( gcx, c_solver_pt, 1. )
          
          ! relative precision
-         
-         rnorme = 0.
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               zgwgt = gcdmat(ji,jj) * gcr(ji,jj)
-               rnorme= rnorme + gcr(ji,jj)*zgwgt
-            END DO
-         END DO
-         
-#if defined key_mpp
-         ! mpp sum over all the global domain
-         CALL  mpp_sum( rnorme )
-#endif
+         rnorme = SUM( gcr(:,:) * gcdmat(:,:) * gcr(:,:) )
+         IF( lk_mpp )   CALL mpp_sum( rnorme )   ! sum over the global domain
          
          ! test of convergence
@@ -216,45 +125,19 @@
          !****
          
-         !,,,,,,,,,,,,,,,,,,,,,,,,,,,,,synchro,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
-         
          ! indicator of non-convergence or explosion
-         
          IF( jn == nmax .OR. SQRT(epsr)/eps > 1.e+20 ) kindic = -2
          IF( ncut == 999 ) GOTO 999
          
-         
-         ! END of iterative loop
-         ! =====================
-         
-      END DO
-      
+         !                                                 ! =====================
+      END DO                                               ! END of iterative loop
+      !                                                    ! =====================
       
 999   CONTINUE
       
       
-      !  2. Output in gcx
-      !  -----------------
+      !  Output in gcx
+      !  -------------
       
-      ! boundary conditions (est-ce necessaire? je ne crois pas!!!!)
-      
-#if defined key_dynspg_fsc
-# if defined key_mpp
-      ! Mpp: export boundary values to neighbouring processors
-      CALL lbc_lnk( gcx, 'S', 1. )
-# else
-      IF( nperio /= 0 ) THEN
-         CALL lbc_lnk( gcx, 'T', 1. )
-      ENDIF
-# endif
-#else
-# if defined key_mpp
-      ! Mpp: export boundary values to neighbouring processors
-      CALL lbc_lnk( gcx, 'G', 1. )
-# else
-      IF( nperio /= 0 ) THEN
-         CALL lbc_lnk( gcx, 'F', 1. )
-      ENDIF
-# endif
-#endif
+      CALL lbc_lnk( gcx, c_solver_pt, 1. )    ! boundary conditions (est-ce necessaire? je ne crois pas!!!!)
       
    END SUBROUTINE sol_sor
Index: trunk/NEMO/OPA_SRC/SOL/solver.F90
===================================================================
--- trunk/NEMO/OPA_SRC/SOL/solver.F90	(revision 15)
+++ trunk/NEMO/OPA_SRC/SOL/solver.F90	(revision 16)
@@ -18,4 +18,7 @@
    USE in_out_manager  ! I/O manager
    USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
+   USE lib_mpp
+   USE dynspg_rl       
+   USE dynspg_fsc      
 
    IMPLICIT NONE
@@ -34,10 +37,10 @@
       !!       * default option: barotropic stream function system
       !!         and islands initialization (if l_isl=T)
-      !!       * key_dynspg_fsc = T : transport divergence system. No specific
+      !!       * lk_dynspg_fsc = T : transport divergence system. No specific
       !!         treatment of islands.
       !!      
       !! ** Method :
       !!       - Compute the local depth of the water column at u- and v-point
-      !!      (key_dynspg_fsc = T) or its inverse (key_dynspg_rl = T).
+      !!      (lk_dynspg_fsc = T) or its inverse (lk_dynspg_rl = T).
       !!      The local depth of the water column is computed by summing 
       !!      the vertical scale factors. For its inverse, the thickness of
@@ -56,7 +59,9 @@
       !!
       !! ** Action : - hur, hvr : masked inverse of the local depth at
-      !!                                u- and v-point. (key_dynspg_rl = T)
+      !!                                u- and v-point. (lk_dynspg_rl = T)
       !!             - hu, hv   : masked local depth at u- and v- points
-      !!                                (key_dynspg_fsc = T)
+      !!                                (lk_dynspg_fsc = T)
+      !!             - c_solver_pt : nature of the gridpoint at which the
+      !!                                solver is applied
       !! References :
       !!      Jensen, 1986: adv. phys. oceanogr. num. mod.,ed. o brien,87-110.
@@ -115,5 +120,5 @@
       ENDIF
 
-#if defined key_dynspg_fsc 
+      IF( lk_dynspg_fsc ) THEN
          IF(lwp) WRITE(numout,*)
          IF(lwp) WRITE(numout,*) '          *** free surface formulation'
@@ -123,14 +128,17 @@
             nstop = nstop + 1
          ENDIF
-#endif
-#if defined key_dynspg_rl
+      ELSEIF( lk_dynspg_rl ) THEN
          IF(lwp) WRITE(numout,*)
          IF(lwp) WRITE(numout,*) '          *** Rigid lid formulation'
-#endif
-#if defined key_dynspg_fsc && defined key_dynspg_rl
+      ELSE
+         IF(lwp) WRITE(numout,cform_err)
+         IF(lwp) WRITE(numout,*) '          Chose at least one surface pressure gradient calculation: free surface or rigid-lid'
+         nstop = nstop + 1
+      ENDIF
+      IF( lk_dynspg_fsc .AND. lk_dynspg_rl ) THEN
          IF(lwp) WRITE(numout,cform_err)
          IF(lwp) WRITE(numout,*) '          Chose between free surface or rigid-lid, not both'
          nstop = nstop + 1
-#endif
+      ENDIF
 
       SELECT CASE ( nsolv )
@@ -144,13 +152,13 @@
       CASE ( 3 )                ! FETI solver
          IF(lwp) WRITE(numout,*) '          use the FETI solver'
-#if ! defined key_mpp
-         IF(lwp) WRITE(numout,*) ' The FETI algorithm is used only with the key_mpp option'
-         nstop = nstop + 1
-#else
-         IF( jpnij == 1 ) THEN
-            IF(lwp) WRITE(numout,*) ' The FETI algorithm needs more than one processor'
+         IF( .NOT.lk_mpp ) THEN
+            IF(lwp) WRITE(numout,*) ' The FETI algorithm is used only with the key_mpp_... option'
             nstop = nstop + 1
-         ENDIF
-#endif
+         ELSE
+            IF( jpnij == 1 ) THEN
+               IF(lwp) WRITE(numout,*) ' The FETI algorithm needs more than one processor'
+               nstop = nstop + 1
+            ENDIF
+         ENDIF
          
       CASE DEFAULT
@@ -161,4 +169,21 @@
       END SELECT
 
+      ! Grid-point at which the solver is applied
+      ! -----------------------------------------
+
+      IF( lk_dynspg_rl ) THEN       ! rigid-lid
+         IF( lk_mpp ) THEN
+            c_solver_pt = 'G'   ! G= F with special staff ??? which one?
+         ELSE
+            c_solver_pt = 'F'
+         ENDIF
+      ELSE                          ! free surface T-point
+         IF( lk_mpp ) THEN
+            c_solver_pt = 'S'   ! S=T with special staff ??? which one?
+         ELSE
+            c_solver_pt = 'T'
+         ENDIF
+      ENDIF
+
 
       ! Construction of the elliptic system matrix
Index: trunk/NEMO/OPA_SRC/TRD/trddyn.F90
===================================================================
--- trunk/NEMO/OPA_SRC/TRD/trddyn.F90	(revision 15)
+++ trunk/NEMO/OPA_SRC/TRD/trddyn.F90	(revision 16)
@@ -34,5 +34,5 @@
 
    !! * Shared module vaiables
-   LOGICAL, PUBLIC, PARAMETER ::   lk_trddyn = .TRUE.   ! momentum trend flag
+   LOGICAL, PUBLIC, PARAMETER ::   lk_trddyn = .TRUE.    !: momentum trend flag
 
    !! * Substitutions
@@ -231,6 +231,6 @@
             DO ji = 1, jpi
                zhke(10) = zhke(10)   &
-               &   + ub(ji,jj,1) * tautrd(ji,jj,1) * e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,1)   &
-               &   + vb(ji,jj,1) * tautrd(ji,jj,2) * e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,1)
+                  &     + ub(ji,jj,1) * tautrd(ji,jj,1) * e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,1)   &
+                  &     + vb(ji,jj,1) * tautrd(ji,jj,2) * e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,1)
             END DO
          END DO
@@ -240,16 +240,16 @@
             DO jj = 1, jpj
                DO ji = 1, jpi
-                  zpeke    = zpeke + zkepe(ji,jj,jk) * g * fsdept(ji,jj,jk)   &
-                  &                        * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk)
-               END DO
-            END DO
-         END DO
-         
-# if defined key_mpp
-         CALL mpp_sum( zpeke )
-         CALL mpp_sum( zumo , 11 )
-         CALL mpp_sum( zvmo , 11 )
-         CALL mpp_sum( zhke , 10 )
-# endif
+                  zpeke    = zpeke + zkepe(ji,jj,jk) * grav * fsdept(ji,jj,jk)   &
+                     &                     * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk)
+               END DO
+            END DO
+         END DO
+         
+         IF( lk_mpp ) THEN
+            CALL mpp_sum( zpeke )
+            CALL mpp_sum( zumo , 11 )
+            CALL mpp_sum( zvmo , 11 )
+            CALL mpp_sum( zhke , 10 )
+         ENDIF
 
 
@@ -423,9 +423,7 @@
          END DO
       END DO
-# if defined key_mpp
-      CALL mpp_sum( tvols )
-      CALL mpp_sum( tvolu )
-      CALL mpp_sum( tvolv )
-# endif
+      IF( lk_mpp )   CALL mpp_sum( tvols )   ! sums over the global domain
+      IF( lk_mpp )   CALL mpp_sum( tvolu )
+      IF( lk_mpp )   CALL mpp_sum( tvolv )
 
       IF(lwp) THEN
@@ -446,8 +444,8 @@
    !!   Default option :                      NO mementum trend diagnostics
    !!----------------------------------------------------------------------
-   LOGICAL, PUBLIC, PARAMETER ::   lk_trddyn = .FALSE.   ! momentum trend flag
+   LOGICAL, PUBLIC, PARAMETER ::   lk_trddyn = .FALSE.   !: momentum trend flag
 CONTAINS
    SUBROUTINE trd_dyn( kt )        ! Empty routine
-      WRITE(*,*) kt
+      WRITE(*,*) 'trd_dyn: You should not have seen this print! error?', kt
    END SUBROUTINE trd_dyn
    SUBROUTINE trd_dyn_init         ! Empty routine
Index: trunk/NEMO/OPA_SRC/TRD/trdmld.F90
===================================================================
--- trunk/NEMO/OPA_SRC/TRD/trdmld.F90	(revision 15)
+++ trunk/NEMO/OPA_SRC/TRD/trdmld.F90	(revision 16)
@@ -37,5 +37,5 @@
 
    !! * Shared module variables
-   LOGICAL, PUBLIC, PARAMETER ::   lk_trdmld = .TRUE.   ! momentum trend flag
+   LOGICAL, PUBLIC, PARAMETER ::   lk_trdmld = .TRUE.    !: momentum trend flag
 
    !! * Module variables
@@ -609,8 +609,8 @@
    !!   Default option :                                       Empty module
    !!----------------------------------------------------------------------
-   LOGICAL, PUBLIC, PARAMETER ::   lk_trdmld = .FALSE.   ! momentum trend flag
+   LOGICAL, PUBLIC, PARAMETER ::   lk_trdmld = .FALSE.   !: momentum trend flag
 CONTAINS
    SUBROUTINE trd_mld( kt )        ! Empty routine
-      WRITE(*,*) kt
+      WRITE(*,*) 'trd_mld: You should not have seen this print! error?', kt
    END SUBROUTINE trd_mld
 #endif
Index: trunk/NEMO/OPA_SRC/TRD/trdtra.F90
===================================================================
--- trunk/NEMO/OPA_SRC/TRD/trdtra.F90	(revision 15)
+++ trunk/NEMO/OPA_SRC/TRD/trdtra.F90	(revision 16)
@@ -33,5 +33,5 @@
 
    !! * Shared module variables
-   LOGICAL, PUBLIC, PARAMETER ::   lk_trdtra = .TRUE.   ! momentum trend flag
+   LOGICAL, PUBLIC, PARAMETER ::   lk_trdtra = .TRUE.    !: momentum trend flag
 
    !! * Substitutions
@@ -95,7 +95,6 @@
             END DO
          END DO
-#if defined key_mpp
-         CALL mpp_sum( tvolt )
-#endif
+         IF( lk_mpp )   CALL mpp_sum( tvolt )   ! sum over the global domain
+
          IF(lwp) THEN
             WRITE(numout,*)
@@ -208,10 +207,10 @@
          END DO
          
-#if defined key_mpp
-         CALL mpp_sum( ztmo, 10 )
-         CALL mpp_sum( zsmo, 10 )
-         CALL mpp_sum( zt2 , 10 )
-         CALL mpp_sum( zs2 , 10 )
-#endif
+         IF( lk_mpp ) THEN
+            CALL mpp_sum( ztmo, 10 )   ! sums over the global domain
+            CALL mpp_sum( zsmo, 10 )
+            CALL mpp_sum( zt2 , 10 )
+            CALL mpp_sum( zs2 , 10 )
+         ENDIF
          
          ! 4. Print
@@ -358,10 +357,7 @@
          END DO
       END DO
-#if defined key_mpp
-      CALL mpp_sum( tvolt )
-#endif
-      IF(lwp) THEN
-         WRITE(numout,*) '          total ocean volume at T-point   tvolt = ',tvolt
-      ENDIF
+      IF( lk_mpp )   CALL mpp_sum( tvolt )   ! sum over the global domain
+
+      IF(lwp) WRITE(numout,*) '          total ocean volume at T-point   tvolt = ',tvolt
 
    END SUBROUTINE trd_tra_init
@@ -371,8 +367,8 @@
    !!   Default case :                                         Empty module
    !!----------------------------------------------------------------------
-   LOGICAL, PUBLIC, PARAMETER ::   lk_trdtra = .FALSE.   ! momentum trend flag
+   LOGICAL, PUBLIC, PARAMETER ::   lk_trdtra = .FALSE.   !: momentum trend flag
 CONTAINS
    SUBROUTINE trd_tra( kt )        ! Empty routine
-      WRITE(*,*) kt
+      WRITE(*,*) 'trd_tra: You should not have seen this print! error?', kt
    END SUBROUTINE trd_tra
    SUBROUTINE trd_tra_init         ! Empty routine
Index: trunk/NEMO/OPA_SRC/TRD/trdtra_oce.F90
===================================================================
--- trunk/NEMO/OPA_SRC/TRD/trdtra_oce.F90	(revision 15)
+++ trunk/NEMO/OPA_SRC/TRD/trdtra_oce.F90	(revision 16)
@@ -12,5 +12,5 @@
    PUBLIC
 
-   INTEGER  ::      & !!! namdia :  diagnostics on dynamics and/or tracer trends
+   INTEGER  ::      & !!: namdia :  diagnostics on dynamics and/or tracer trends
       ntrd  = 10 ,  &  !: time step frequency dynamics and tracers trends
       nctls =  0       !: control surface type for trends vertical integration
@@ -25,5 +25,5 @@
    !! Trends diagnostics parameters
    !!---------------------------------------------------------------------
-   INTEGER, PARAMETER ::            &
+   INTEGER, PARAMETER ::            &  !:
 # if defined key_traldf_eiv
       jptrdh = 4,   &  !: number of 3D horiz trends arrays
@@ -38,7 +38,7 @@
    !! Trends diagnostics variables
    !!---------------------------------------------------------------------
-   REAL(wp) ::   &
-      tvolt         ! volume of the whole ocean computed at t-points
-   REAL(wp), DIMENSION(jpi,jpj,jpk,7) ::   &
+   REAL(wp) ::   &  !:
+      tvolt         !: volume of the whole ocean computed at t-points
+   REAL(wp), DIMENSION(jpi,jpj,jpk,7) ::   &  !:
       ttrd             !: trends of the temperature tracer equations
       !                !  ttrd(,,,1) : horizontal advection
@@ -49,13 +49,13 @@
       !                !  ttrd(,,,6) : damping OR vertical EIV
       !                !  ttrd(,,,7) : penetrative solar radiation (T only)
-   REAL(wp), DIMENSION(jpi,jpj,jpk,6) ::   &
+   REAL(wp), DIMENSION(jpi,jpj,jpk,6) ::   &  !:
       strd             !: trends of the salinity tracer equations
       !                !  same as ttrd()
-   REAL(wp), DIMENSION(jpi,jpj,jpk,jptrdh) ::   &
+   REAL(wp), DIMENSION(jpi,jpj,jpk,jptrdh) ::   &  !:
       ttrdh, strdh     !: ttrdh(,,,1) : zonal advection
       !                !  ttrdh(,,,2) : meridional advection
       !                !  ttrdh(,,,3) : zonal EIV
       !                !  ttrdh(,,,4) : meridional EIV
-   REAL(wp), DIMENSION(jpi,jpj,2) ::   &
+   REAL(wp), DIMENSION(jpi,jpj,2) ::   &  !:
       flxtrd,       &  !: tracer forcing trends
       bbltrd           !: tracer bottom boundary layer trends
Index: trunk/NEMO/OPA_SRC/ZDF/zdf_oce.F90
===================================================================
--- trunk/NEMO/OPA_SRC/ZDF/zdf_oce.F90	(revision 15)
+++ trunk/NEMO/OPA_SRC/ZDF/zdf_oce.F90	(revision 16)
@@ -15,5 +15,5 @@
 
    !! * Share Module variables
-   LOGICAL, PARAMETER, PUBLIC ::    &
+   LOGICAL, PARAMETER, PUBLIC ::    &   !:
 #if defined key_zdfcst   ||   defined key_esopa
       lk_zdfcst        = .TRUE.         !: constant vertical mixing flag
@@ -21,9 +21,9 @@
       lk_zdfcst        = .FALSE.        !: constant vertical mixing flag
 #endif
-   LOGICAL, PUBLIC ::    &
+   LOGICAL, PUBLIC ::    &   !:
       ln_zdfevd        = .TRUE.  ,   &  !: convection: enhanced vertical diffusion flag
       ln_zdfnpc        = .FALSE.        !: convection: non-penetrative convection flag
 
-   LOGICAL, PUBLIC ::    &
+   LOGICAL, PUBLIC ::    &   !:
       l_trazdf_exp     = .FALSE. ,   &  !: ???
       l_trazdf_imp     = .FALSE. ,   &  !: 
@@ -32,19 +32,19 @@
       l_dynzdf_imp_tsk = .FALSE.        !:
 
-   INTEGER, PUBLIC ::    & !!! namzdf:  vertical diffusion
+   INTEGER, PUBLIC ::    & !!: namzdf:  vertical diffusion
       n_zdfexp = 3    ,  &  !: number of sub-time step (explicit time stepping)
       nevdm    = 1          !: =0/1 flag to apply enhanced avm or not
  
-   REAL(wp), PUBLIC ::   & !!! namzdf   vertical diffusion
+   REAL(wp), PUBLIC ::   & !!: namzdf   vertical diffusion
       avm0  = 1.e-4_wp,  &  !: vertical eddy viscosity (m2/s)
       avt0  = 1.e-5_wp,  &  !: vertical eddy diffusivity (m2/s)
       avevd = 1._wp         !: vertical eddy coeff. for enhanced vert. diff. (m2/s)
 
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   &
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   &   !:
       avmu,              &  !: vertical viscosity coeff. at uw-, vw-points
       avmv,              &  !: vertical viscosity coeff. at uw-, vw-points
       avt                   !: vertical diffusivity coeff. at w-point
  
-   REAL(wp), PUBLIC, DIMENSION(jpk) ::   &
+   REAL(wp), PUBLIC, DIMENSION(jpk) ::   &   !:
       avmb, avtb            !: background profile of avm and avt
  
Index: trunk/NEMO/OPA_SRC/ZDF/zdfddm.F90
===================================================================
--- trunk/NEMO/OPA_SRC/ZDF/zdfddm.F90	(revision 15)
+++ trunk/NEMO/OPA_SRC/ZDF/zdfddm.F90	(revision 16)
@@ -25,6 +25,6 @@
 
    !! * Shared module variables
-   LOGICAL, PUBLIC ::   lk_zdfddm = .TRUE.    !: double diffusive mixing flag
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   &
+   LOGICAL, PUBLIC, PARAMETER ::   lk_zdfddm = .TRUE.    !: double diffusive mixing flag
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   &   !:
       avs ,               &  !: salinity vertical diffusivity coeff. at w-point
       rrau                   !: heat/salt buoyancy flux ratio
@@ -242,8 +242,8 @@
    !!   Default option :          Dummy module          No double diffusion
    !!----------------------------------------------------------------------
-   LOGICAL, PUBLIC ::   lk_zdfddm = .FALSE.   !: double diffusion flag
+   LOGICAL, PUBLIC, PARAMETER ::   lk_zdfddm = .FALSE.   !: double diffusion flag
 CONTAINS
    SUBROUTINE zdf_ddm( kt )           ! Dummy routine
-      WRITE(*,*) kt                          ! avoid compil warning
+      WRITE(*,*) 'zdf_ddm: You should not have seen this print! error?', kt
    END SUBROUTINE zdf_ddm
 #endif
Index: trunk/NEMO/OPA_SRC/ZDF/zdfmxl.F90
===================================================================
--- trunk/NEMO/OPA_SRC/ZDF/zdfmxl.F90	(revision 15)
+++ trunk/NEMO/OPA_SRC/ZDF/zdfmxl.F90	(revision 16)
@@ -21,7 +21,7 @@
 
    !! * Shared module variables
-   INTEGER, PUBLIC, DIMENSION(jpi,jpj) ::   &
+   INTEGER, PUBLIC, DIMENSION(jpi,jpj) ::   &   !:
       nmln                  !: number of level in the mixed layer
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   &
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   &   !:
       hmld ,             &  !: mixing layer depth (turbocline) (m)
       hmlp ,             &  !: mixed layer depth  (rho=rho0+zdcrit) (m)
Index: trunk/NEMO/OPA_SRC/ZDF/zdfric.F90
===================================================================
--- trunk/NEMO/OPA_SRC/ZDF/zdfric.F90	(revision 15)
+++ trunk/NEMO/OPA_SRC/ZDF/zdfric.F90	(revision 16)
@@ -28,5 +28,5 @@
 
    !! * Shared module variables
-   LOGICAL, PUBLIC, PARAMETER ::   lk_zdfric = .TRUE.   !: Richardson vertical mixing flag
+   LOGICAL, PUBLIC, PARAMETER ::   lk_zdfric = .TRUE.    !: Richardson vertical mixing flag
 
    !! * Module variables
@@ -256,5 +256,5 @@
 CONTAINS
    SUBROUTINE zdf_ric( kt )        ! Dummy routine
-      WRITE(*,*) kt 
+      WRITE(*,*) 'zdf_ric: You should not have seen this print! error?', kt
    END SUBROUTINE zdf_ric
 #endif
Index: trunk/NEMO/OPA_SRC/ZDF/zdftke.F90
===================================================================
--- trunk/NEMO/OPA_SRC/ZDF/zdftke.F90	(revision 15)
+++ trunk/NEMO/OPA_SRC/ZDF/zdftke.F90	(revision 16)
@@ -28,9 +28,9 @@
 
    !! * Share Module variables
-   LOGICAL, PUBLIC, PARAMETER ::   lk_zdftke = .TRUE.   !: TKE vertical mixing flag
-   LOGICAL, PUBLIC ::         & !!! ** tke namelist (namtke) **
+   LOGICAL, PUBLIC, PARAMETER ::   lk_zdftke = .TRUE.    !: TKE vertical mixing flag
+   LOGICAL, PUBLIC ::         & !!: ** tke namelist (namtke) **
      ln_rstke = .FALSE.          !: =T restart with tke from a run without tke with 
      !                           !  a none zero initial value for en
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   &
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   &   !:
       en                         !: now turbulent kinetic energy
 
@@ -92,5 +92,5 @@
       !!         d(en)/dt = eboost eav (d(u)/dz)**2       ! shear production
       !!                  + d( efave eav d(en)/dz )/dz    ! diffusion of tke
-      !!                  + g/rau0 pdl eav d(rau)/dz      ! stratif. destruc.
+      !!                  + grav/rau0 pdl eav d(rau)/dz   ! stratif. destruc.
       !!                  - ediss / emxl en**(2/3)        ! dissipation
       !!      with the boundary conditions:
@@ -791,8 +791,8 @@
    !!   Dummy module :                                        NO TKE scheme
    !!----------------------------------------------------------------------
-   LOGICAL, PUBLIC, PARAMETER ::   lk_zdftke = .FALSE.   ! TKE flag
+   LOGICAL, PUBLIC, PARAMETER ::   lk_zdftke = .FALSE.   !: TKE flag
 CONTAINS
    SUBROUTINE zdf_tke( kt )          ! Empty routine
-      WRITE(*,*) kt                     ! no compilation warning
+      WRITE(*,*) 'zdf_tke: You should not have seen this print! error?', kt
    END SUBROUTINE zdf_tke
 #endif
