Index: trunk/NEMO/OPA_SRC/DIA/diafwb.F90
===================================================================
--- trunk/NEMO/OPA_SRC/DIA/diafwb.F90	(revision 83)
+++ trunk/NEMO/OPA_SRC/DIA/diafwb.F90	(revision 84)
@@ -109,6 +109,6 @@
       IF( lk_mpp )   CALL mpp_sum( a_rnf    )       ! sum over the global domain
 
-      IF( aminus /= 0.0 ) a_aminus = a_aminus + ( MIN( aplus, aminus ) / aminus )
-      IF( aplus  /= 0.0 ) a_aplus  = a_aplus  + ( MIN( aplus, aminus ) / aplus  )
+      IF( aminus /= 0.e0 ) a_aminus = a_aminus + ( MIN( aplus, aminus ) / aminus )
+      IF( aplus  /= 0.e0 ) a_aplus  = a_aplus  + ( MIN( aplus, aminus ) / aplus  )
 
       IF( kt == nitend ) THEN
@@ -205,5 +205,5 @@
          zu = un(ji,jj,jk) * fse3t(ji,jj,jk) * e2u(ji,jj)
 
-         IF( un(ji,jj,jk) > 0.0 ) THEN 
+         IF( un(ji,jj,jk) > 0.e0 ) THEN 
             zflxi(1) = zflxi(1) +    zu
             ztemi(1) = ztemi(1) + zt*zu
@@ -245,5 +245,5 @@
          zu = un(ji,jj,jk) * fse3t(ji,jj,jk) * e2u(ji,jj)
          
-         IF( un(ji,jj,jk) > 0.0 ) THEN 
+         IF( un(ji,jj,jk) > 0.e0 ) THEN 
             zflxi(2) = zflxi(2) +    zu
             ztemi(2) = ztemi(2) + zt*zu
@@ -285,5 +285,5 @@
          zu = un(ji,jj,jk) * fse3t(ji,jj,jk) * e2u(ji,jj)
          
-         IF( un(ji,jj,jk) > 0.0 ) THEN 
+         IF( un(ji,jj,jk) > 0.e0 ) THEN 
             zflxi(3) = zflxi(3) +    zu
             ztemi(3) = ztemi(3) + zt*zu
@@ -325,5 +325,5 @@
          zu = un(ji,jj,jk) * fse3t(ji,jj,jk) * e2u(ji,jj)
          
-         IF( un(ji,jj,jk) > 0.0 ) THEN 
+         IF( un(ji,jj,jk) > 0.e0 ) THEN 
             zflxi(4) = zflxi(4) +    zu
             ztemi(4) = ztemi(4) + zt*zu
@@ -338,5 +338,5 @@
       ! Sum at each time-step
       DO jt = 1, 4 
-         IF( zflxi(jt) /= 0.0 .AND. zflxo(jt) /= 0.0 ) THEN
+         IF( zflxi(jt) /= 0.e0 .AND. zflxo(jt) /= 0.e0 ) THEN
             a_flxi(jt) = a_flxi(jt) + zflxi(jt)
             a_temi(jt) = a_temi(jt) + ztemi(jt)/zflxi(jt)
@@ -350,10 +350,10 @@
       IF( kt == nitend ) THEN
          DO jt = 1, 4 
-            a_flxi(jt) = a_flxi(jt)/((nitend-nit000+1)*1.e6)
-            a_temi(jt) = a_temi(jt)/( nitend-nit000+1)
-            a_sali(jt) = a_sali(jt)/( nitend-nit000+1)
-            a_flxo(jt) = a_flxo(jt)/((nitend-nit000+1)*1.e6)
-            a_temo(jt) = a_temo(jt)/( nitend-nit000+1)
-            a_salo(jt) = a_salo(jt)/( nitend-nit000+1)
+            a_flxi(jt) = a_flxi(jt) / ( FLOAT( nitend - nit000 + 1 ) * 1.e6 )
+            a_temi(jt) = a_temi(jt) /   FLOAT( nitend - nit000 + 1 )
+            a_sali(jt) = a_sali(jt) /   FLOAT( nitend - nit000 + 1 )
+            a_flxo(jt) = a_flxo(jt) / ( FLOAT( nitend - nit000 + 1 ) * 1.e6 )
+            a_temo(jt) = a_temo(jt) /   FLOAT( nitend - nit000 + 1 )
+            a_salo(jt) = a_salo(jt) /   FLOAT( nitend - nit000 + 1 )
          END DO
       ENDIF
@@ -368,7 +368,7 @@
          WRITE(111,*)
          WRITE(111,*)    'Net freshwater budget '
-         WRITE(111,9010) '  emp    = ',a_emp,   ' m3 =', a_emp   /((nitend-nit000+1)*rdttra(1)) * 1.e-6,' Sv'
-         WRITE(111,9010) '  precip = ',a_precip,' m3 =', a_precip/((nitend-nit000+1)*rdttra(1)) * 1.e-6,' Sv'
-         WRITE(111,9010) '  a_rnf   = ',a_rnf,   ' m3 =', a_rnf   /((nitend-nit000+1)*rdttra(1)) * 1.e-6,' Sv'
+         WRITE(111,9010) '  emp    = ',a_emp,   ' m3 =', a_emp   /(FLOAT(nitend-nit000+1)*rdttra(1)) * 1.e-6,' Sv'
+         WRITE(111,9010) '  precip = ',a_precip,' m3 =', a_precip/(FLOAT(nitend-nit000+1)*rdttra(1)) * 1.e-6,' Sv'
+         WRITE(111,9010) '  a_rnf  = ',a_rnf,   ' m3 =', a_rnf   /(FLOAT(nitend-nit000+1)*rdttra(1)) * 1.e-6,' Sv'
          WRITE(111,*)
          WRITE(111,9010) '  zarea =',zarea
@@ -378,5 +378,5 @@
          WRITE(111,9010) '  at nit000 = ',a_sshb        ,' m3 '
          WRITE(111,9010) '  at nitend = ',a_sshn        ,' m3 '
-         WRITE(111,9010) '  diff      = ',(a_sshn-a_sshb),' m3 =', (a_sshn-a_sshb)/((nitend-nit000+1)*rdt) * 1.e-6,' Sv'
+         WRITE(111,9010) '  diff      = ',(a_sshn-a_sshb),' m3 =', (a_sshn-a_sshb)/(FLOAT(nitend-nit000+1)*rdt) * 1.e-6,' Sv'
          WRITE(111,9020) '  mean sea level elevation    =', a_sshn/zarea,' m'
          WRITE(111,*)
Index: trunk/NEMO/OPA_SRC/DIA/diawri.F90
===================================================================
--- trunk/NEMO/OPA_SRC/DIA/diawri.F90	(revision 83)
+++ trunk/NEMO/OPA_SRC/DIA/diawri.F90	(revision 84)
@@ -364,6 +364,4 @@
             CALL histdef( nid_W,"voddmavs","Salt Vertical Eddy Diffusivity"    , "m2/s"   ,   &  ! avs
                &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
-            CALL histdef( nid_W,"voddmrra","Heat/Salt buoyancy Ratio"          , "-"      ,   &  ! rrau
-               &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
          ENDIF
          !                                                                                      !!! nid_W : 2D
Index: trunk/NEMO/OPA_SRC/DOM/domwri.F90
===================================================================
--- trunk/NEMO/OPA_SRC/DOM/domwri.F90	(revision 83)
+++ trunk/NEMO/OPA_SRC/DOM/domwri.F90	(revision 84)
@@ -6,5 +6,5 @@
 
    !!----------------------------------------------------------------------
-   !!   dom_wri        : create mesh and mask file (s)
+   !!   dom_wri        : create and write mesh and mask file(s)
    !!                    nmsh = 1  :   mesh_mask file
    !!                         = 2  :   mesh and mask file
Index: trunk/NEMO/OPA_SRC/DTA/dtasst.F90
===================================================================
--- trunk/NEMO/OPA_SRC/DTA/dtasst.F90	(revision 83)
+++ trunk/NEMO/OPA_SRC/DTA/dtasst.F90	(revision 84)
@@ -96,4 +96,5 @@
          IF(lwp) WRITE(numout,*) 'dta_sst : DAILY sea surface temperature data'
          IF(lwp) WRITE(numout,*) '~~~~~~~   read in file: ', clname
+         sst(:,:) = 0.e0   ! required for extra halos in mpp
       ENDIF
 
Index: trunk/NEMO/OPA_SRC/DYN/dynhpg.F90
===================================================================
--- trunk/NEMO/OPA_SRC/DYN/dynhpg.F90	(revision 83)
+++ trunk/NEMO/OPA_SRC/DYN/dynhpg.F90	(revision 84)
@@ -162,5 +162,5 @@
       END DO
 
-      IF( l_ctl .AND. lwp ) THEN         ! print sum trends (used for debugging)
+      IF(l_ctl) THEN         ! print sum trends (used for debugging)
          zuap = SUM( ua(2:jpim1,2:jpjm1,1:jpkm1) * umask(2:jpim1,2:jpjm1,1:jpkm1) )
          zvap = SUM( va(2:jpim1,2:jpjm1,1:jpkm1) * vmask(2:jpim1,2:jpjm1,1:jpkm1) )
@@ -317,5 +317,5 @@
       END DO
 
-      IF( l_ctl .AND. lwp ) THEN         ! print sum trends (used for debugging)
+      IF(l_ctl) THEN         ! print sum trends (used for debugging)
          zuap = SUM( ua(2:jpim1,2:jpjm1,1:jpkm1) * umask(2:jpim1,2:jpjm1,1:jpkm1) )
          zvap = SUM( va(2:jpim1,2:jpjm1,1:jpkm1) * vmask(2:jpim1,2:jpjm1,1:jpkm1) )
@@ -428,5 +428,5 @@
       END DO
 
-      IF( l_ctl .AND. lwp ) THEN         ! print sum trends (used for debugging)
+      IF(l_ctl) THEN         ! print sum trends (used for debugging)
          zuap = SUM( ua(2:jpim1,2:jpjm1,1:jpkm1) * umask(2:jpim1,2:jpjm1,1:jpkm1) )
          zvap = SUM( va(2:jpim1,2:jpjm1,1:jpkm1) * vmask(2:jpim1,2:jpjm1,1:jpkm1) )
Index: trunk/NEMO/OPA_SRC/DYN/dynhpg_atsk.F90
===================================================================
--- trunk/NEMO/OPA_SRC/DYN/dynhpg_atsk.F90	(revision 83)
+++ trunk/NEMO/OPA_SRC/DYN/dynhpg_atsk.F90	(revision 84)
@@ -150,5 +150,5 @@
       !                                                ! ===============
 
-      IF( l_ctl .AND. lwp ) THEN         ! print sum trends (used for debugging)
+      IF(l_ctl) THEN         ! print sum trends (used for debugging)
          zuap = SUM( ua(2:jpim1,2:jpjm1,1:jpkm1) * umask(2:jpim1,2:jpjm1,1:jpkm1) )
          zvap = SUM( va(2:jpim1,2:jpjm1,1:jpkm1) * vmask(2:jpim1,2:jpjm1,1:jpkm1) )
@@ -295,4 +295,12 @@
       END DO                                           !   End of slab
       !                                                ! ===============
+
+      IF(l_ctl) THEN         ! print sum trends (used for debugging)
+         zuap = SUM( ua(2:jpim1,2:jpjm1,1:jpkm1) * umask(2:jpim1,2:jpjm1,1:jpkm1) )
+         zvap = SUM( va(2:jpim1,2:jpjm1,1:jpkm1) * vmask(2:jpim1,2:jpjm1,1:jpkm1) )
+         WRITE(numout,*) ' hpg  - Ua: ', zuap-u_ctl, ' Va: ', zvap-v_ctl
+         u_ctl = zuap   ;   v_ctl = zvap
+      ENDIF   
+
    END SUBROUTINE dyn_hpg_atsk
 
@@ -402,5 +410,5 @@
       !                                                ! ===============
 
-      IF( l_ctl .AND. lwp ) THEN         ! print sum trends (used for debugging)
+      IF(l_ctl) THEN         ! print sum trends (used for debugging)
          WRITE(numout,*) ' hpg  - Ua: ', SUM(ua*umask), ' Va: ', SUM(va*vmask)
          zuap = SUM( ua(2:jpim1,2:jpjm1,1:jpkm1) * umask(2:jpim1,2:jpjm1,1:jpkm1) )
Index: trunk/NEMO/OPA_SRC/DYN/dynkeg.F90
===================================================================
--- trunk/NEMO/OPA_SRC/DYN/dynkeg.F90	(revision 83)
+++ trunk/NEMO/OPA_SRC/DYN/dynkeg.F90	(revision 84)
@@ -119,5 +119,5 @@
       !                                                ! ===============
 
-      IF( l_ctl .AND. lwp ) THEN         ! print sum trends (used for debugging)
+      IF(l_ctl) THEN         ! print sum trends (used for debugging)
          zua = SUM( ua(2:jpim1,2:jpjm1,1:jpkm1) * umask(2:jpim1,2:jpjm1,1:jpkm1) )
          zva = SUM( va(2:jpim1,2:jpjm1,1:jpkm1) * vmask(2:jpim1,2:jpjm1,1:jpkm1) )
Index: trunk/NEMO/OPA_SRC/DYN/dynldf_bilap.F90
===================================================================
--- trunk/NEMO/OPA_SRC/DYN/dynldf_bilap.F90	(revision 83)
+++ trunk/NEMO/OPA_SRC/DYN/dynldf_bilap.F90	(revision 84)
@@ -219,4 +219,12 @@
       END DO                                           !   End of slab
       !                                                ! ===============
+
+      IF(l_ctl) THEN         ! print sum trends (used for debugging)
+         zua = SUM( ua(2:jpim1,2:jpjm1,1:jpkm1) * umask(2:jpim1,2:jpjm1,1:jpkm1) )
+         zva = SUM( va(2:jpim1,2:jpjm1,1:jpkm1) * vmask(2:jpim1,2:jpjm1,1:jpkm1) )
+         WRITE(numout,*) ' ldf  - Ua: ', zua-u_ctl, ' Va: ', zva-v_ctl
+         u_ctl = zua   ;   v_ctl = zva
+      ENDIF
+
    END SUBROUTINE dyn_ldf_bilap
 
Index: trunk/NEMO/OPA_SRC/DYN/dynldf_bilapg.F90
===================================================================
--- trunk/NEMO/OPA_SRC/DYN/dynldf_bilapg.F90	(revision 83)
+++ trunk/NEMO/OPA_SRC/DYN/dynldf_bilapg.F90	(revision 84)
@@ -73,4 +73,5 @@
       !! * Local declarations
       INTEGER ::   ji, jj, jk                 ! dummy loop indices
+      REAL(wp) ::   zua, zva                  ! temporary scalars
       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   &
          wk1, wk2,            &  ! work array used for rotated biharmonic
@@ -126,4 +127,11 @@
       END DO                                           !   End of slab
       !                                                ! ===============
+
+      IF(l_ctl) THEN         ! print sum trends (used for debugging)
+         zua = SUM( ua(2:jpim1,2:jpjm1,1:jpkm1) * umask(2:jpim1,2:jpjm1,1:jpkm1) )
+         zva = SUM( va(2:jpim1,2:jpjm1,1:jpkm1) * vmask(2:jpim1,2:jpjm1,1:jpkm1) )
+         WRITE(numout,*) ' ldf  - Ua: ', zua-u_ctl, ' Va: ', zva-v_ctl
+         u_ctl = zua   ;   v_ctl = zva
+      ENDIF
 
    END SUBROUTINE dyn_ldf_bilapg
Index: trunk/NEMO/OPA_SRC/DYN/dynldf_lap.F90
===================================================================
--- trunk/NEMO/OPA_SRC/DYN/dynldf_lap.F90	(revision 83)
+++ trunk/NEMO/OPA_SRC/DYN/dynldf_lap.F90	(revision 84)
@@ -118,5 +118,5 @@
       !                                                ! ===============
 
-      IF( l_ctl .AND. lwp ) THEN         ! print sum trends (used for debugging)
+      IF(l_ctl) THEN         ! print sum trends (used for debugging)
          zua = SUM( ua(2:jpim1,2:jpjm1,1:jpkm1) * umask(2:jpim1,2:jpjm1,1:jpkm1) )
          zva = SUM( va(2:jpim1,2:jpjm1,1:jpkm1) * vmask(2:jpim1,2:jpjm1,1:jpkm1) )
Index: trunk/NEMO/OPA_SRC/DYN/dynnxt.F90
===================================================================
--- trunk/NEMO/OPA_SRC/DYN/dynnxt.F90	(revision 83)
+++ trunk/NEMO/OPA_SRC/DYN/dynnxt.F90	(revision 84)
@@ -133,7 +133,5 @@
       !                                                ! ===============
 
-      IF( l_ctl .AND. lwp ) THEN         ! print sum fields (used for debugging)
-         WRITE(numout,*) ' nxt  - Un: ', SUM(un*umask), ' Vn: ', SUM(vn*vmask)
-      ENDIF
+      IF(l_ctl)   WRITE(numout,*) ' nxt  - Un: ', SUM(un*umask), ' Vn: ', SUM(vn*vmask)
 
    END SUBROUTINE dyn_nxt
Index: trunk/NEMO/OPA_SRC/DYN/dynspg_fsc.F90
===================================================================
--- trunk/NEMO/OPA_SRC/DYN/dynspg_fsc.F90	(revision 83)
+++ trunk/NEMO/OPA_SRC/DYN/dynspg_fsc.F90	(revision 84)
@@ -347,5 +347,5 @@
       END DO
 
-      IF( l_ctl .AND. lwp ) THEN         ! print sum trends (used for debugging)
+      IF(l_ctl) THEN         ! print sum trends (used for debugging)
          WRITE(numout,*) ' spg  - Ua: ', SUM( ua(2:jpim1,2:jpjm1,1:jpkm1)*umask(2:jpim1,2:jpjm1,1:jpkm1) ),   &
             &                   ' Va: ', SUM( va(2:jpim1,2:jpjm1,1:jpkm1)*vmask(2:jpim1,2:jpjm1,1:jpkm1) )
Index: trunk/NEMO/OPA_SRC/DYN/dynspg_fsc_atsk.F90
===================================================================
--- trunk/NEMO/OPA_SRC/DYN/dynspg_fsc_atsk.F90	(revision 83)
+++ trunk/NEMO/OPA_SRC/DYN/dynspg_fsc_atsk.F90	(revision 84)
@@ -18,9 +18,10 @@
    USE oce             ! ocean dynamics and tracers 
    USE dom_oce         ! ocean space and time domain
-   USE trdtra_oce     ! ocean active tracer trend 
-   USE trddyn_oce     ! ocean dynamics trend
+   USE zdf_oce         ! ocean vertical physics
+   USE trdtra_oce      ! ocean active tracer trend 
+   USE trddyn_oce      ! ocean dynamics trend
    USE in_out_manager  ! I/O manager
    USE phycst          ! physical constant
-   USE ocesbc          ! Ocen Surface Boundary condition
+   USE ocesbc          ! Ocean Surface Boundary condition
    USE flxrnf          ! ???
    USE sol_oce         ! ocean elliptic solver
@@ -30,5 +31,4 @@
    USE obc_oce         ! Lateral open boundary condition
    USE obcdyn          ! open boudary condition
-   USE obcdyn          !    "              "
    USE obcvol          !    "              "
    USE lib_mpp         ! ???
@@ -362,4 +362,9 @@
          END DO
 
+         IF(l_ctl) THEN         ! print sum trends (used for debugging)
+            WRITE(numout,*) ' spg  - Ua: ', SUM( ua(2:jpim1,2:jpjm1,1:jpkm1)*umask(2:jpim1,2:jpjm1,1:jpkm1) ),   &
+               &                   ' Va: ', SUM( va(2:jpim1,2:jpjm1,1:jpkm1)*vmask(2:jpim1,2:jpjm1,1:jpkm1) )
+         ENDIF
+
          ! 8. Sea surface elevation time stepping
          ! --------------------------------------
Index: trunk/NEMO/OPA_SRC/DYN/dynzad.F90
===================================================================
--- trunk/NEMO/OPA_SRC/DYN/dynzad.F90	(revision 83)
+++ trunk/NEMO/OPA_SRC/DYN/dynzad.F90	(revision 84)
@@ -126,5 +126,5 @@
       !                                                ! ===============
 
-      IF( l_ctl .AND. lwp ) THEN         ! print sum trends (used for debugging)
+      IF(l_ctl) THEN         ! print sum trends (used for debugging)
          zua = SUM( ua(2:jpim1,2:jpjm1,1:jpkm1) * umask(2:jpim1,2:jpjm1,1:jpkm1) )
          zva = SUM( va(2:jpim1,2:jpjm1,1:jpkm1) * vmask(2:jpim1,2:jpjm1,1:jpkm1) )
@@ -227,5 +227,5 @@
       END DO
 
-      IF( l_ctl .AND. lwp ) THEN         ! print sum trends (used for debugging)
+      IF(l_ctl) THEN         ! print sum trends (used for debugging)
          zua = SUM( ua(2:jpim1,2:jpjm1,1:jpkm1) * umask(2:jpim1,2:jpjm1,1:jpkm1) )
          zva = SUM( va(2:jpim1,2:jpjm1,1:jpkm1) * vmask(2:jpim1,2:jpjm1,1:jpkm1) )
Index: trunk/NEMO/OPA_SRC/DYN/dynzdf_imp.F90
===================================================================
--- trunk/NEMO/OPA_SRC/DYN/dynzdf_imp.F90	(revision 83)
+++ trunk/NEMO/OPA_SRC/DYN/dynzdf_imp.F90	(revision 84)
@@ -342,5 +342,5 @@
       END DO
 
-      IF( l_ctl .AND. lwp ) THEN         ! print sum trends (used for debugging)
+      IF(l_ctl) THEN         ! print sum trends (used for debugging)
          zua = SUM( ua(2:jpim1,2:jpjm1,1:jpkm1) * umask(2:jpim1,2:jpjm1,1:jpkm1) )
          zva = SUM( va(2:jpim1,2:jpjm1,1:jpkm1) * vmask(2:jpim1,2:jpjm1,1:jpkm1) )
Index: trunk/NEMO/OPA_SRC/DYN/dynzdf_imp_atsk.F90
===================================================================
--- trunk/NEMO/OPA_SRC/DYN/dynzdf_imp_atsk.F90	(revision 83)
+++ trunk/NEMO/OPA_SRC/DYN/dynzdf_imp_atsk.F90	(revision 84)
@@ -312,4 +312,12 @@
       END DO                                           !   End of slab
       !                                                ! ===============
+
+      IF(l_ctl) THEN         ! print sum trends (used for debugging)
+         zua = SUM( ua(2:jpim1,2:jpjm1,1:jpkm1) * umask(2:jpim1,2:jpjm1,1:jpkm1) )
+         zva = SUM( va(2:jpim1,2:jpjm1,1:jpkm1) * vmask(2:jpim1,2:jpjm1,1:jpkm1) )
+         WRITE(numout,*) ' zdf  - Ua: ', zua-u_ctl, ' Va: ', zva-v_ctl
+         u_ctl = zua   ;   v_ctl = zva
+      ENDIF
+
    END SUBROUTINE dyn_zdf_imp_tsk
 
Index: trunk/NEMO/OPA_SRC/FLO/floblk.F90
===================================================================
--- trunk/NEMO/OPA_SRC/FLO/floblk.F90	(revision 83)
+++ trunk/NEMO/OPA_SRC/FLO/floblk.F90	(revision 84)
@@ -348,5 +348,5 @@
 # if defined key_obc
       DO jfl = 1, jpnfl
-         IF( lpeastobc ) THEN
+         IF( lp_obc_east ) THEN
             IF( jped <=  zgjfl(jfl) .AND. zgjfl(jfl) <= jpef .AND. nieob-1 <=  zgifl(jfl) ) THEN
                zgifl (jfl) = INT(zgifl(jfl)) + 0.5
@@ -355,5 +355,5 @@
             END IF
          END IF
-         IF( lpwestobc ) THEN
+         IF( lp_obc_west ) THEN
             IF( jpwd <= zgjfl(jfl) .AND. zgjfl(jfl) <= jpwf .AND. niwob >=  zgifl(jfl) ) THEN
                zgifl (jfl) = INT(zgifl(jfl)) + 0.5
@@ -362,5 +362,5 @@
             END IF
          END IF
-         IF( lpnorthobc ) THEN
+         IF( lp_obc_north ) THEN
             IF( jpnd <=  zgifl(jfl) .AND. zgifl(jfl) <= jpnf .AND. njnob-1 >=  zgjfl(jfl) ) THEN
                zgifl (jfl) = INT(zgifl(jfl)) + 0.5
@@ -369,5 +369,5 @@
             END IF
          END IF
-         IF( lpsouthobc ) THEN
+         IF( lp_obc_south ) THEN
             IF( jpsd <=  zgifl(jfl) .AND. zgifl(jfl) <= jpsf .AND.  njsob >= zgjfl(jfl) ) THEN
                zgifl (jfl) = INT(zgifl(jfl)) + 0.5
@@ -408,5 +408,4 @@
       ENDIF
 
-      RETURN
    END SUBROUTINE flo_blk
 
Index: trunk/NEMO/OPA_SRC/SBC/bulk.F90
===================================================================
--- trunk/NEMO/OPA_SRC/SBC/bulk.F90	(revision 83)
+++ trunk/NEMO/OPA_SRC/SBC/bulk.F90	(revision 84)
@@ -95,5 +95,5 @@
 
 # if ! defined key_ice_lim
-         IF( l_ctl .AND. lwp ) THEN         ! print mean trends (used for debugging)
+         IF(l_ctl) THEN         ! print mean trends (used for debugging)
             WRITE(numout,*) ' Forcings '
             WRITE(numout,*) ' qsr_oce  : ', SUM( qsr_oce (:,:) * tmask(:,:,1) )
Index: trunk/NEMO/OPA_SRC/SBC/flxblk.F90
===================================================================
--- trunk/NEMO/OPA_SRC/SBC/flxblk.F90	(revision 83)
+++ trunk/NEMO/OPA_SRC/SBC/flxblk.F90	(revision 84)
@@ -33,55 +33,55 @@
 
    !! * Module variables
-      INTEGER, PARAMETER  ::   &
-         jpintsr = 24          ! number of time step between sunrise and sunset
-         !                     ! uses for heat flux computation
-      LOGICAL ::   &
-         lbulk_init = .TRUE.   ! flag, bulk initialization done or not)
-
-      REAL(wp), DIMENSION(jpi,jpj) ::   &
-         stauc            ,  &   ! cloud optical depth 
-         sbudyko   
-
-      !! * constants for bulk computation (flx_blk)
-      REAL(wp), DIMENSION(19)  ::  &
-         budyko                  ! BUDYKO's coefficient
-      ! BUDYKO's coefficient (cloudiness effect on LW radiation):
-      DATA budyko / 1.00, 0.98, 0.95, 0.92, 0.89, 0.86, 0.83, 0.80, 0.78, 0.75,  &
-                    0.72, 0.69, 0.67, 0.64, 0.61, 0.58, 0.56, 0.53, 0.50 /
-      REAL(wp), DIMENSION(20)  :: &
-         tauco                  ! cloud optical depth coefficient
-      ! Cloud optical depth coefficient
-      DATA tauco / 6.6, 6.6, 7.0, 7.2, 7.1, 6.8, 6.5, 6.6, 7.1, 7.6,   &
-                   6.6, 6.1, 5.6, 5.5, 5.8, 5.8, 5.6, 5.6, 5.6, 5.6 /
-      REAL(wp)  ::            &  ! constant values
-         zeps    = 1e-20   ,  &
-         zeps0   = 1e-13   ,  &
-         zeps1   = 1e-06   ,  &
-         zzero   = 0.0     ,  &
-         zone    = 1.0
-
-      !! * constants for albedo computation (flx_blk_albedo)
-      REAL(wp) ::   &
-         c1     = 0.05  ,     &   ! constants values
-         c2     = 0.1   ,     &
-         albice = 0.50  ,     &   !  albedo of melting ice in the arctic and antarctic (Shine & Hendersson-Sellers)
-         cgren  = 0.06  ,     &   !  correction of the snow or ice albedo to take into account
-                                  !  effects of cloudiness (Grenfell & Perovich, 1984)
-         alphd  = 0.80  ,     &   !  coefficients for linear interpolation used to compute
-         alphdi = 0.72  ,     &   !  albedo between two extremes values (Pyane, 1972)
-         alphc  = 0.65  ,     &
-         zmue   = 0.4             !  cosine of local solar altitude
-
-      !! * constants for solar declinaison computation (flx_blk_declin)
-      REAL(wp) ::                &
-         a0  =  0.39507671   ,   &  ! coefficients
-         a1  = 22.85684301   ,   &
-         a2  = -0.38637317   ,   &
-         a3  =  0.15096535   ,   &
-         a4  = -0.00961411   ,   &
-         b1  = -4.29692073   ,   &
-         b2  =  0.05702074   ,   &
-         b3  = -0.09028607   ,   &
-         b4  =  0.00592797
+   INTEGER, PARAMETER  ::   &
+      jpintsr = 24          ! number of time step between sunrise and sunset
+      !                     ! uses for heat flux computation
+   LOGICAL ::   &
+      lbulk_init = .TRUE.   ! flag, bulk initialization done or not)
+
+   REAL(wp), DIMENSION(jpi,jpj) ::   &
+      stauc            ,  &   ! cloud optical depth 
+      sbudyko   
+
+   !! * constants for bulk computation (flx_blk)
+   REAL(wp), DIMENSION(19)  ::  &
+      budyko                  ! BUDYKO's coefficient
+   ! BUDYKO's coefficient (cloudiness effect on LW radiation):
+   DATA budyko / 1.00, 0.98, 0.95, 0.92, 0.89, 0.86, 0.83, 0.80, 0.78, 0.75,  &
+      &          0.72, 0.69, 0.67, 0.64, 0.61, 0.58, 0.56, 0.53, 0.50 /
+   REAL(wp), DIMENSION(20)  :: &
+      tauco                  ! cloud optical depth coefficient
+   ! Cloud optical depth coefficient
+   DATA tauco / 6.6, 6.6, 7.0, 7.2, 7.1, 6.8, 6.5, 6.6, 7.1, 7.6,   &
+      &         6.6, 6.1, 5.6, 5.5, 5.8, 5.8, 5.6, 5.6, 5.6, 5.6 /
+   REAL(wp)  ::            &  ! constant values
+      zeps    = 1.e-20  ,  &
+      zeps0   = 1.e-13  ,  &
+      zeps1   = 1.e-06  ,  &
+      zzero   = 0.e0    ,  &
+      zone    = 1.0
+
+   !! * constants for albedo computation (flx_blk_albedo)
+   REAL(wp) ::   &
+      c1     = 0.05  ,     &   ! constants values
+      c2     = 0.10  ,     &
+      albice = 0.50  ,     &   !  albedo of melting ice in the arctic and antarctic (Shine & Hendersson-Sellers)
+      cgren  = 0.06  ,     &   !  correction of the snow or ice albedo to take into account
+                               !  effects of cloudiness (Grenfell & Perovich, 1984)
+      alphd  = 0.80  ,     &   !  coefficients for linear interpolation used to compute
+      alphdi = 0.72  ,     &   !  albedo between two extremes values (Pyane, 1972)
+      alphc  = 0.65  ,     &
+      zmue   = 0.40            !  cosine of local solar altitude
+
+   !! * constants for solar declinaison computation (flx_blk_declin)
+   REAL(wp) ::                &
+      a0  =  0.39507671   ,   &  ! coefficients
+      a1  = 22.85684301   ,   &
+      a2  = -0.38637317   ,   &
+      a3  =  0.15096535   ,   &
+      a4  = -0.00961411   ,   &
+      b1  = -4.29692073   ,   &
+      b2  =  0.05702074   ,   &
+      b3  = -0.09028607   ,   &
+      b4  =  0.00592797
    !!----------------------------------------------------------------------
    !!   OPA 9.0 , LODYC-IPSL  (2003)
@@ -237,21 +237,21 @@
          DO jj = 1, jpj  
             DO ji = 1 , jpi
-            zalat          = ( 90.0 - ABS( gphit(ji,jj) ) ) / 5.0
-            zclat          = ( 95.0 -      gphit(ji,jj)   ) / 10.0
-            indxb          = 1 + INT( zalat ) 
-            !  correction factor to account for the effect of clouds 
-            sbudyko(ji,jj) = budyko(indxb)  
-            indxc          = 1 + INT( zclat )  
-            zdl            = zclat - INT( zclat ) 
-            zdr            = 1.0 - zdl
-            stauc(ji,jj)   = zdr * tauco( indxc ) + zdl * tauco( indxc + 1 ) 
+               zalat          = ( 90.e0 - ABS( gphit(ji,jj) ) ) /  5.e0
+               zclat          = ( 95.e0 -      gphit(ji,jj)   ) / 10.e0
+               indxb          = 1 + INT( zalat ) 
+               !  correction factor to account for the effect of clouds 
+               sbudyko(ji,jj) = budyko(indxb)  
+               indxc          = 1 + INT( zclat )  
+               zdl            = zclat - INT( zclat ) 
+               zdr            = 1.0 - zdl
+               stauc(ji,jj)   = zdr * tauco( indxc ) + zdl * tauco( indxc + 1 ) 
             END DO
          END DO
          IF( nleapy == 1 ) THEN
-            yearday = 366.0
+            yearday = 366.e0
          ELSE IF( nleapy == 0 ) THEN
-            yearday = 365.0
+            yearday = 365.e0
          ELSEIF( nleapy == 30) THEN
-            yearday = 360.0
+            yearday = 360.e0
          ENDIF
          lbulk_init = .FALSE.
@@ -265,21 +265,21 @@
       zqsb_ice(:,:) = 0.e0
 
-      zpis2       = rpi / 2.              ! pi / 2
-      z2pi        = 2. * rpi              ! 2 * pi 
+      zpis2       = rpi / 2.
+      z2pi        = 2. * rpi
 
  !CDIR NOVERRCHK
-     DO jj = 1, jpj
+      DO jj = 1, jpj
  !CDIR NOVERRCHK
-        DO ji = 1, jpi
-
-      ztatm (ji,jj) = 273.15 + tatm  (ji,jj)  !  air temperature in Kelvins 
-      zcatm1(ji,jj) = 1.0    - catm  (ji,jj)  !  fractional cloud cover
-      zfrld (ji,jj) = 1.0    - freeze(ji,jj)  !  fractional sea ice cover
-      zpatm(ji,jj)  = 101000.               !  pressure 
-      
-      !  Computation of air density, obtained from the equation of state for dry air. 
-      zrhoa(ji,jj) = zpatm(ji,jj) / ( 287.04 * ztatm(ji,jj) )
-      
-      !  zes : Saturation water vapour
+         DO ji = 1, jpi
+
+            ztatm (ji,jj) = 273.15 + tatm  (ji,jj)  !  air temperature in Kelvins 
+            zcatm1(ji,jj) = 1.0    - catm  (ji,jj)  !  fractional cloud cover
+            zfrld (ji,jj) = 1.0    - freeze(ji,jj)  !  fractional sea ice cover
+            zpatm(ji,jj)  = 101000.               !  pressure 
+      
+            !  Computation of air density, obtained from the equation of state for dry air. 
+            zrhoa(ji,jj) = zpatm(ji,jj) / ( 287.04 * ztatm(ji,jj) )
+      
+            !  zes : Saturation water vapour
             ztamr = ztatm(ji,jj) - rtt
             zmt1  = SIGN( 17.269, ztamr )
@@ -289,16 +289,16 @@
                &                      / ( ztatm(ji,jj) - 35.86  + MAX( zzero, zmt3 ) ) )
 
-      !  zev : vapour pressure  (hatm is relative humidity)  
-      zev(ji,jj)   = hatm(ji,jj) * zes(ji,jj) 
-      !  square-root of vapour pressure
+            !  zev : vapour pressure  (hatm is relative humidity)  
+            zev(ji,jj)   = hatm(ji,jj) * zes(ji,jj) 
+            !  square-root of vapour pressure
 !CDIR NOVERRCHK
-      zevsqr(ji,jj) = SQRT( zev(ji,jj) * 0.01 )
-      !  zqapb  : specific humidity 
-      zqatm(ji,jj) = 0.622 * zev(ji,jj) / ( zpatm(ji,jj) - 0.378 * zev(ji,jj) )
-
-
-      !----------------------------------------------------
-      !   Computation of snow precipitation (Ledley, 1985) |
-      !----------------------------------------------------
+            zevsqr(ji,jj) = SQRT( zev(ji,jj) * 0.01 )
+            !  zqapb  : specific humidity 
+            zqatm(ji,jj) = 0.622 * zev(ji,jj) / ( zpatm(ji,jj) - 0.378 * zev(ji,jj) )
+
+
+            !----------------------------------------------------
+            !   Computation of snow precipitation (Ledley, 1985) |
+            !----------------------------------------------------
 
             zmt1  =   253.0 - ztatm(ji,jj)
@@ -333,5 +333,5 @@
 
       iday   = INT( zxday )
-      IF( l_ctl .AND. lwp ) WRITE(numout,*) ' declin : iday ', iday, ' nfbulk= ', nfbulk
+      IF(l_ctl)   WRITE(numout,*) ' declin : iday ', iday, ' nfbulk= ', nfbulk
       !   computation of the solar declination, his sine and his cosine
       CALL flx_blk_declin( indaet, iday, zdecl )
@@ -351,18 +351,18 @@
          DO ji = 1, jpi
             !  product of sine of latitude and sine of solar declination
-            zps(ji,jj)     = SIN( gphit(ji,jj) * rad ) * zsdecl
+            zps     (ji,jj) = SIN( gphit(ji,jj) * rad ) * zsdecl
             !  product of cosine of latitude and cosine of solar declination
-            zpc(ji,jj)     = COS( gphit(ji,jj) * rad ) * zcdecl
+            zpc     (ji,jj) = COS( gphit(ji,jj) * rad ) * zcdecl
             !  computation of the both local time of sunrise and sunset
-            zlsrise(ji,jj) =  ACOS( - SIGN( zone, zps(ji,jj) ) * MIN( zone, SIGN( zone, zps(ji,jj) )  &
+            zlsrise (ji,jj) = ACOS( - SIGN( zone, zps(ji,jj) ) * MIN( zone, SIGN( zone, zps(ji,jj) )  &
                &                     * ( zps(ji,jj) / zpc(ji,jj) ) ) ) 
-            zlsset(ji,jj)  = - zlsrise(ji,jj)
+            zlsset  (ji,jj) = - zlsrise(ji,jj)
             !  dividing the solar day into jpintsr segments of length zdlha
-            zdlha(ji,jj)   =  ( zlsrise(ji,jj) - zlsset(ji,jj) ) / REAL( jpintsr )
+            zdlha   (ji,jj) = ( zlsrise(ji,jj) - zlsset(ji,jj) ) / REAL( jpintsr )
             !  computation of the local noon solar altitude
-            zlmunoon(ji,jj)=  ASIN ( ( zps(ji,jj) + zpc(ji,jj) ) ) / rad
+            zlmunoon(ji,jj) = ASIN ( ( zps(ji,jj) + zpc(ji,jj) ) ) / rad
             
             !  cloud correction taken from Reed (1977) (imposed lower than 1)
-            zcldcor(ji,jj) = MIN( zone, ( 1 - 0.62 * catm(ji,jj) + 0.0019 * zlmunoon(ji,jj) ) )
+            zcldcor (ji,jj) = MIN( zone, ( 1.e0 - 0.62 * catm(ji,jj) + 0.0019 * zlmunoon(ji,jj) ) )
          END DO
       END DO
@@ -380,5 +380,5 @@
             DO ji = 1, jpi
                !  local hour angle
-               zlha (ji,jj,jt)= COS ( zlsrise(ji,jj) - zcoef * zdlha(ji,jj) )
+               zlha (ji,jj,jt) = COS ( zlsrise(ji,jj) - zcoef * zdlha(ji,jj) )
 
                ! cosine of local solar altitude
@@ -417,26 +417,26 @@
       !--------------------------------------------------------------------
 
-      zalbocsd(:,:) = 0.  
-      zqsro   (:,:) = 0.
-      zqsrics (:,:) = 0.
-      zqsrios (:,:) = 0.
+      zalbocsd(:,:) = 0.e0
+      zqsro   (:,:) = 0.e0
+      zqsrics (:,:) = 0.e0
+      zqsrios (:,:) = 0.e0
 
       DO jt = 1, jpintsr 
 #   if defined key_vectopt_loop && ! defined key_autotasking
          DO ji = 1, jpij  
-            zalbocsd(ji,1)  = zalbocsd(ji,1) + zdlha(ji,1) * zalbocs(ji,1,jt)   &
-               &             / MAX( 2.0 * zlsrise(ji,1) , zeps0 )
-            zqsro   (ji,1)  =   zqsro(ji,1) + zsqsro  (ji,1,jt)
-            zqsrics (ji,1)  = zqsrics(ji,1) + zsqsrics(ji,1,jt)
-            zqsrios (ji,1)  = zqsrios(ji,1) + zsqsrios(ji,1,jt)
+            zalbocsd(ji,1) = zalbocsd(ji,1) + zdlha   (ji,1) * zalbocs(ji,1,jt)   &
+               &                                             / MAX( 2.0 * zlsrise(ji,1) , zeps0 )
+            zqsro   (ji,1) = zqsro   (ji,1) + zsqsro  (ji,1,jt)
+            zqsrics (ji,1) = zqsrics (ji,1) + zsqsrics(ji,1,jt)
+            zqsrios (ji,1) = zqsrios (ji,1) + zsqsrios(ji,1,jt)
          END DO
 #  else
          DO jj = 1, jpj
             DO ji = 1, jpi  
-               zalbocsd(ji,jj)  = zalbocsd(ji,jj) + zdlha(ji,jj) * zalbocs(ji,jj,jt)   &
-                  &             / MAX( 2.0 * zlsrise(ji,jj) , zeps0 )
-               zqsro  (ji,jj)   =   zqsro(ji,jj) + zsqsro  (ji,jj,jt)
-               zqsrics(ji,jj)   = zqsrics(ji,jj) + zsqsrics(ji,jj,jt)
-               zqsrios(ji,jj)   = zqsrios(ji,jj) + zsqsrios(ji,jj,jt)
+               zalbocsd(ji,jj) = zalbocsd(ji,jj) + zdlha(ji,jj) * zalbocs(ji,jj,jt)   &
+                  &                                              / MAX( 2.0 * zlsrise(ji,jj) , zeps0 )
+               zqsro  (ji,jj)  = zqsro   (ji,jj) + zsqsro  (ji,jj,jt)
+               zqsrics(ji,jj)  = zqsrics (ji,jj) + zsqsrics(ji,jj,jt)
+               zqsrios(ji,jj)  = zqsrios (ji,jj) + zsqsrios(ji,jj,jt)
             END DO
          END DO
@@ -447,7 +447,7 @@
          DO ji = 1, jpi 
 
-      !------------------------------------------- 
-      !  Computation of shortwave radiation.
-      !-------------------------------------------
+            !------------------------------------------- 
+            !  Computation of shortwave radiation.
+            !-------------------------------------------
 
             ! the solar heat flux absorbed at ocean and snow/ice surfaces
@@ -478,7 +478,7 @@
             fr2_i0(ji,jj) = 0.82  * zcatm1(ji,jj) + 0.65 * catm(ji,jj)
 
-      !---------------------------------------------------------------------------
-      !   Computation of long-wave radiation  ( Berliand 1952 ; all latitudes )
-      !---------------------------------------------------------------------------
+            !---------------------------------------------------------------------------
+            !   Computation of long-wave radiation  ( Berliand 1952 ; all latitudes )
+            !---------------------------------------------------------------------------
 
             ! tempory variables
@@ -633,5 +633,5 @@
             !  latent heat flux 
             zqla_ice(ji,jj) = zrhovacshi * ( zqsati        - zqatm(ji,jj) )
-             qla_ice(ji,jj) = zqla_ice(ji,jj)
+            qla_ice (ji,jj) = zqla_ice(ji,jj)
               
             !  Computation of sensitivity of non solar fluxes (dQ/dT)
@@ -750,5 +750,5 @@
       !-------------------------- 
       
-      llmask = (hsnif == 0.0) .AND. ( sist >= rt0_ice )
+      llmask = (hsnif == 0.e0) .AND. ( sist >= rt0_ice )
       WHERE ( llmask )   !  ice free of snow and melts
          zalbfz = albice
Index: trunk/NEMO/OPA_SRC/SBC/flxfwb.F90
===================================================================
--- trunk/NEMO/OPA_SRC/SBC/flxfwb.F90	(revision 83)
+++ trunk/NEMO/OPA_SRC/SBC/flxfwb.F90	(revision 84)
@@ -130,6 +130,6 @@
       IF( lk_mpp )   CALL  mpp_sum( a_rnf    )   ! sum over the global domain
 
-      IF( aminus /= 0.0 ) a_aminus = a_aminus + ( MIN( aplus, aminus ) / aminus )
-      IF( aplus  /= 0.0 ) a_aplus  = a_aplus  + ( MIN( aplus, aminus ) / aplus  )
+      IF( aminus /= 0.e0 ) a_aminus = a_aminus + ( MIN( aplus, aminus ) / aminus )
+      IF( aplus  /= 0.e0 ) a_aplus  = a_aplus  + ( MIN( aplus, aminus ) / aplus  )
 
 
Index: trunk/NEMO/OPA_SRC/SBC/ocesbc.F90
===================================================================
--- trunk/NEMO/OPA_SRC/SBC/ocesbc.F90	(revision 83)
+++ trunk/NEMO/OPA_SRC/SBC/ocesbc.F90	(revision 84)
@@ -6,5 +6,6 @@
 
    !!----------------------------------------------------------------------
-   !!    oce_sbc : initialization and namelist read
+   !!   oce_sbc     : ???
+   !!   oce_sbc_dmp : ???
    !!----------------------------------------------------------------------
    !! * Modules used
@@ -34,8 +35,8 @@
 
    !! * Shared module variables
-   REAL(wp), PUBLIC ::   &
+   REAL(wp), PUBLIC ::   &  !:
       aplus, aminus,     &  !:
       empold = 0.e0         !: current year freshwater budget correction
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   &
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   &  !:
       qt  ,         &  !: total surface heat flux (w/m2)
       q   ,         &  !: surface heat flux (w/m2)
@@ -158,10 +159,10 @@
 
          ! Re-initialization of fluxes
-         sst_io(:,:) = 0.0
-         sss_io(:,:) = 0.0
-         u_io  (:,:) = 0.0
-         v_io  (:,:) = 0.0
-         gtaux (:,:) = 0.
-         gtauy (:,:) = 0.
+         sst_io(:,:) = 0.e0
+         sss_io(:,:) = 0.e0
+         u_io  (:,:) = 0.e0
+         v_io  (:,:) = 0.e0
+         gtaux (:,:) = 0.e0
+         gtauy (:,:) = 0.e0
 
       ENDIF
@@ -257,8 +258,8 @@
 
          ! Re-initialization of fluxes
-         sst_io(:,:) = 0.0
-         sss_io(:,:) = 0.0
-         u_io  (:,:) = 0.0
-         v_io  (:,:) = 0.0
+         sst_io(:,:) = 0.e0
+         sss_io(:,:) = 0.e0
+         u_io  (:,:) = 0.e0
+         v_io  (:,:) = 0.e0
 
       ENDIF
@@ -651,5 +652,5 @@
       DO jj = 1, jpj
          DO ji = 1, jpi
-            freezn(ji,jj) = MAX(0., SIGN(1., freeze(ji,jj)-rsmall))
+            freezn(ji,jj) = MAX(0., SIGN(1., freeze(ji,jj)-rsmall) )
          END DO
       END DO
@@ -683,5 +684,6 @@
       END DO
       ! volume flux associated to internal damping to climatology
-      dmp(:,:) = zstrdmp(:,:) * rauw / ( zsss(:,:) + rsmall )
+!!ibu dmp(:,:) = zstrdmp(:,:) * rauw / ( zsss(:,:) + rsmall )
+      dmp(:,:) = zstrdmp(:,:) * rauw / ( zsss(:,:) + 1.e-20 )
 # else
       dmp(:,:) = 0.e0            ! No internal damping
@@ -698,5 +700,6 @@
             zerp = ( 1. - 2.*upsrnfh(ji,jj) ) * zsrp   &
                & * ( zsss(ji,jj) - s_dta(ji,jj,1) )     &
-               & / ( zsss(ji,jj) + rsmall        )
+               & / ( zsss(ji,jj) + 1.e-20        )
+!ib            & / ( zsss(ji,jj) + rsmall        )
             
             zerp = MIN( zerp, zplus  )
@@ -719,8 +722,8 @@
       IF( lk_mpp )   CALL mpp_sum( aplus  )   ! sums over the global domain
       IF( lk_mpp )   CALL mpp_sum( aminus )
-      IF( l_ctl .AND. lwp ) WRITE(numout,*) ' oce_sbc_dmp : a+ = ', aplus, ' a- = ', aminus
+      IF(l_ctl)   WRITE(numout,*) ' oce_sbc_dmp : a+ = ', aplus, ' a- = ', aminus
 
       zadefi = MIN( aplus, aminus )
-      IF( zadefi == 0.0 ) THEN 
+      IF( zadefi == 0.e0 ) THEN 
          erp(:,:) = 0.e0
       ELSE
@@ -732,5 +735,6 @@
       erp(:,:) = ( 1. - zfreeze(:,:) ) * zsrp    &   ! surface restoring term
          &     * ( zsss(:,:) - s_dta(:,:,1) )     &
-         &     / ( zsss(:,:) + rsmall      )
+         &     / ( zsss(:,:) + 1.e-20      )
+!ib      &     / ( zsss(:,:) + rsmall      )
 #endif
 
