Index: trunk/NEMO/OPA_SRC/DYN/dynhpg.F90
===================================================================
--- trunk/NEMO/OPA_SRC/DYN/dynhpg.F90	(revision 456)
+++ trunk/NEMO/OPA_SRC/DYN/dynhpg.F90	(revision 474)
@@ -203,9 +203,6 @@
       IF( ln_hpg_djc )   ioptio = ioptio + 1
       IF( ln_hpg_rot )   ioptio = ioptio + 1
-      IF ( ioptio > 1 ) THEN
-          IF(lwp) WRITE(numout,cform_err)
-          IF(lwp) WRITE(numout,*) ' several hydrostatic pressure gradient options used'
-          nstop = nstop + 1
-      ENDIF
+      IF ( ioptio > 1 )   &
+           &   CALL ctl_stop( ' several hydrostatic pressure gradient options used' )
 
       IF( lk_dynhpg_jki ) THEN
Index: trunk/NEMO/OPA_SRC/DYN/dynldf.F90
===================================================================
--- trunk/NEMO/OPA_SRC/DYN/dynldf.F90	(revision 456)
+++ trunk/NEMO/OPA_SRC/DYN/dynldf.F90	(revision 474)
@@ -152,18 +152,10 @@
       IF( ln_dynldf_lap   )   ioptio = ioptio + 1
       IF( ln_dynldf_bilap )   ioptio = ioptio + 1
-      IF( ioptio /= 1 )   THEN
-          IF(lwp) WRITE(numout,cform_err)
-          IF(lwp) WRITE(numout,*) '          use ONE of the 2 lap/bilap operator type on dynamics'
-          nstop = nstop + 1
-      ENDIF
+      IF( ioptio /= 1 ) CALL ctl_stop( '          use ONE of the 2 lap/bilap operator type on dynamics' )
       ioptio = 0
       IF( ln_dynldf_level )   ioptio = ioptio + 1
       IF( ln_dynldf_hor   )   ioptio = ioptio + 1
       IF( ln_dynldf_iso   )   ioptio = ioptio + 1
-      IF( ioptio /= 1 ) THEN
-         IF(lwp) WRITE(numout,cform_err)
-         IF(lwp) WRITE(numout,*) '          use only ONE direction (level/hor/iso)'
-         nstop = nstop + 1
-      ENDIF
+      IF( ioptio /= 1 ) CALL ctl_stop( '          use only ONE direction (level/hor/iso)' )
 
       ! defined the type of lateral diffusion from ln_dynldf_... logicals
@@ -205,20 +197,11 @@
       ENDIF
 
-      IF( ierr == 1 ) THEN
-         IF(lwp) WRITE(numout,cform_err)
-         IF(lwp) WRITE(numout,*) ' iso-level in z-coordinate - partial step, not allowed'
-         nstop = nstop + 1
-      ENDIF
-      IF( ierr == 2 ) THEN
-         IF(lwp) WRITE(numout,cform_err)
-         IF(lwp) WRITE(numout,*) ' isoneutral bilaplacian operator does not exist'
-         nstop = nstop + 1
-      ENDIF
+      IF( ierr == 1 )   &
+           &   CALL ctl_stop( ' iso-level in z-coordinate - partial step, not allowed' )
+      IF( ierr == 2 )   &
+           &   CALL ctl_stop( ' isoneutral bilaplacian operator does not exist' )
       IF( nldf == 1 .OR. nldf == 3 ) THEN      ! rotation
-         IF( .NOT.lk_ldfslp ) THEN
-            IF(lwp) WRITE(numout,cform_err)
-            IF(lwp) WRITE(numout,*) '          the rotation of the diffusive tensor require key_ldfslp'
-            nstop = nstop + 1
-         ENDIF
+         IF( .NOT.lk_ldfslp )   &
+           &   CALL ctl_stop(  '          the rotation of the diffusive tensor require key_ldfslp' )
       ENDIF
 
Index: trunk/NEMO/OPA_SRC/DYN/dynldf_bilap.F90
===================================================================
--- trunk/NEMO/OPA_SRC/DYN/dynldf_bilap.F90	(revision 456)
+++ trunk/NEMO/OPA_SRC/DYN/dynldf_bilap.F90	(revision 474)
@@ -86,5 +86,7 @@
       REAL(wp) ::   zua, zva, zbt, ze2u, ze2v ! temporary scalar
       REAL(wp), DIMENSION(jpi,jpj) ::   &
-         zuf, zut, zlu, zlv, zcu, zcv         ! temporary workspace
+         zcu, zcv                             ! temporary workspace
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   &
+         zuf, zut, zlu, zlv                   ! temporary workspace
       !!----------------------------------------------------------------------
       !!  OPA 8.5, LODYC-IPSL (2002)
@@ -96,8 +98,14 @@
          IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~'
       ENDIF
-      zuf(:,:) = 0.e0
-      zut(:,:) = 0.e0
-      zlu(:,:) = 0.e0
-      zlv(:,:) = 0.e0
+
+!!bug gm this should be enough
+!!$      zuf(:,:,jpk) = 0.e0
+!!$      zut(:,:,jpk) = 0.e0
+!!$      zlu(:,:,jpk) = 0.e0
+!!$      zlv(:,:,jpk) = 0.e0
+      zuf(:,:,:) = 0.e0
+      zut(:,:,:) = 0.e0
+      zlu(:,:,:) = 0.e0
+      zlv(:,:,:) = 0.e0
 
       !                                                ! ===============
@@ -108,11 +116,11 @@
 
          IF( ln_sco .OR. ln_zps ) THEN   ! s-coordinate or z-coordinate with partial steps
-            zuf(:,:) = rotb(:,:,jk) * fse3f(:,:,jk)
+            zuf(:,:,jk) = rotb(:,:,jk) * fse3f(:,:,jk)
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  zlu(ji,jj) = - ( zuf(ji,jj) - zuf(ji,jj-1) ) / ( e2u(ji,jj) * fse3u(ji,jj,jk) )   &
+                  zlu(ji,jj,jk) = - ( zuf(ji,jj,jk) - zuf(ji,jj-1,jk) ) / ( e2u(ji,jj) * fse3u(ji,jj,jk) )   &
                      &         + ( hdivb(ji+1,jj,jk) - hdivb(ji,jj,jk) ) / e1u(ji,jj)
    
-                  zlv(ji,jj) = + ( zuf(ji,jj) - zuf(ji-1,jj) ) / ( e1v(ji,jj) * fse3v(ji,jj,jk) )   &
+                  zlv(ji,jj,jk) = + ( zuf(ji,jj,jk) - zuf(ji-1,jj,jk) ) / ( e1v(ji,jj) * fse3v(ji,jj,jk) )   &
                      &         + ( hdivb(ji,jj+1,jk) - hdivb(ji,jj,jk) ) / e2v(ji,jj)
                END DO
@@ -121,33 +129,36 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  zlu(ji,jj) = - ( rotb (ji  ,jj,jk) - rotb (ji,jj-1,jk) ) / e2u(ji,jj)   &
+                  zlu(ji,jj,jk) = - ( rotb (ji  ,jj,jk) - rotb (ji,jj-1,jk) ) / e2u(ji,jj)   &
                      &         + ( hdivb(ji+1,jj,jk) - hdivb(ji,jj  ,jk) ) / e1u(ji,jj)
    
-                  zlv(ji,jj) = + ( rotb (ji,jj  ,jk) - rotb (ji-1,jj,jk) ) / e1v(ji,jj)   &
+                  zlv(ji,jj,jk) = + ( rotb (ji,jj  ,jk) - rotb (ji-1,jj,jk) ) / e1v(ji,jj)   &
                      &         + ( hdivb(ji,jj+1,jk) - hdivb(ji  ,jj,jk) ) / e2v(ji,jj)
                END DO  
             END DO  
          ENDIF
-
-         ! Boundary conditions on the laplacian  (zlu,zlv)
-         CALL lbc_lnk( zlu, 'U', -1. )
-         CALL lbc_lnk( zlv, 'V', -1. )
-         
-         
+      ENDDO
+
+      ! Boundary conditions on the laplacian  (zlu,zlv)
+      CALL lbc_lnk( zlu, 'U', -1. )
+      CALL lbc_lnk( zlv, 'V', -1. )
+         
+         
+      DO jk = 1, jpkm1
+   
          ! Third derivative
          ! ----------------
          
          ! Multiply by the eddy viscosity coef. (at u- and v-points)
-         zlu(:,:) = zlu(:,:) * fsahmu(:,:,jk)
-         zlv(:,:) = zlv(:,:) * fsahmv(:,:,jk)
+         zlu(:,:,jk) = zlu(:,:,jk) * fsahmu(:,:,jk)
+         zlv(:,:,jk) = zlv(:,:,jk) * fsahmv(:,:,jk)
          
          ! Contravariant "laplacian"
-         zcu(:,:) = e1u(:,:) * zlu(:,:)
-         zcv(:,:) = e2v(:,:) * zlv(:,:)
+         zcu(:,:) = e1u(:,:) * zlu(:,:,jk)
+         zcv(:,:) = e2v(:,:) * zlv(:,:,jk)
          
          ! Laplacian curl ( * e3f if s-coordinates or z-coordinate with partial steps)
          DO jj = 1, jpjm1
             DO ji = 1, fs_jpim1   ! vector opt.
-               zuf(ji,jj) = fmask(ji,jj,jk) * (  zcv(ji+1,jj  ) - zcv(ji,jj)      &
+               zuf(ji,jj,jk) = fmask(ji,jj,jk) * (  zcv(ji+1,jj  ) - zcv(ji,jj)      &
                   &                            - zcu(ji  ,jj+1) + zcu(ji,jj)  )   &
 #if defined key_zco
@@ -163,9 +174,9 @@
             DO ji = 1, fs_jpim1   ! vector opt.
 #if defined key_zco
-               zlu(ji,jj) = e2u(ji,jj) * zlu(ji,jj)
-               zlv(ji,jj) = e1v(ji,jj) * zlv(ji,jj)
-#else
-               zlu(ji,jj) = e2u(ji,jj) * fse3u(ji,jj,jk) * zlu(ji,jj)
-               zlv(ji,jj) = e1v(ji,jj) * fse3v(ji,jj,jk) * zlv(ji,jj)
+               zlu(ji,jj,jk) = e2u(ji,jj) * zlu(ji,jj,jk)
+               zlv(ji,jj,jk) = e1v(ji,jj) * zlv(ji,jj,jk)
+#else
+               zlu(ji,jj,jk) = e2u(ji,jj) * fse3u(ji,jj,jk) * zlu(ji,jj,jk)
+               zlv(ji,jj,jk) = e1v(ji,jj) * fse3v(ji,jj,jk) * zlv(ji,jj,jk)
 #endif
             END DO
@@ -180,15 +191,18 @@
                zbt = e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk)
 #endif
-               zut(ji,jj) = (  zlu(ji,jj) - zlu(ji-1,jj  )   &
-                  &          + zlv(ji,jj) - zlv(ji  ,jj-1)  ) / zbt
+               zut(ji,jj,jk) = (  zlu(ji,jj,jk) - zlu(ji-1,jj  ,jk)   &
+                  &             + zlv(ji,jj,jk) - zlv(ji  ,jj-1,jk) ) / zbt
             END DO
          END DO
+      END DO
 
 
       ! boundary conditions on the laplacian curl and div (zuf,zut)
+!!bug gm no need to do this 2 following lbc...
       CALL lbc_lnk( zuf, 'F', 1. )
       CALL lbc_lnk( zut, 'T', 1. )
 
-         
+      DO jk = 1, jpkm1      
+   
          ! Bilaplacian
          ! -----------
@@ -204,9 +218,9 @@
 #endif
                ! horizontal biharmonic diffusive trends
-               zua = - ( zuf(ji  ,jj) - zuf(ji,jj-1) ) / ze2u   &
-                  &  + ( zut(ji+1,jj) - zut(ji,jj  ) ) / e1u(ji,jj)
-
-               zva = + ( zuf(ji,jj  ) - zuf(ji-1,jj) ) / ze2v   &
-                  &  + ( zut(ji,jj+1) - zut(ji  ,jj) ) / e2v(ji,jj)
+               zua = - ( zuf(ji  ,jj,jk) - zuf(ji,jj-1,jk) ) / ze2u   &
+                  &  + ( zut(ji+1,jj,jk) - zut(ji,jj  ,jk) ) / e1u(ji,jj)
+
+               zva = + ( zuf(ji,jj  ,jk) - zuf(ji-1,jj,jk) ) / ze2v   &
+                  &  + ( zut(ji,jj+1,jk) - zut(ji  ,jj,jk) ) / e2v(ji,jj)
                ! add it to the general momentum trends
                ua(ji,jj,jk) = ua(ji,jj,jk) + zua
Index: trunk/NEMO/OPA_SRC/DYN/dynspg.F90
===================================================================
--- trunk/NEMO/OPA_SRC/DYN/dynspg.F90	(revision 456)
+++ trunk/NEMO/OPA_SRC/DYN/dynspg.F90	(revision 474)
@@ -168,9 +168,6 @@
       IF(lk_dynspg_rl )   ioptio = ioptio + 1
 
-      IF( ( ioptio > 1 .AND. .NOT. lk_esopa ) .OR. ioptio == 0 ) THEN
-         IF(lwp) WRITE(numout,cform_err)
-         IF(lwp) WRITE(numout,*) ' Choose only one surface pressure gradient scheme with a key cpp'
-         nstop = nstop + 1
-      ENDIF
+      IF( ( ioptio > 1 .AND. .NOT. lk_esopa ) .OR. ioptio == 0 )   &
+           &   CALL ctl_stop( ' Choose only one surface pressure gradient scheme with a key cpp' )
 
       IF( lk_esopa     )   nspg = -1
@@ -199,9 +196,6 @@
       ! --------------------------
       IF( lk_dynspg_ts ) THEN
-         IF( MOD( rdt , rdtbt ) /= 0. ) THEN
-            IF(lwp) WRITE(numout,cform_err)
-            IF(lwp) WRITE(numout,*) ' The barotropic timestep must be an integer divisor of the baroclinic timestep'
-            nstop = nstop + 1
-         ENDIF
+         IF( MOD( rdt , rdtbt ) /= 0. )   &
+           &   CALL ctl_stop( ' The barotropic timestep must be an integer divisor of the baroclinic timestep' )
       ENDIF
 
Index: trunk/NEMO/OPA_SRC/DYN/dynspg_flt.F90
===================================================================
--- trunk/NEMO/OPA_SRC/DYN/dynspg_flt.F90	(revision 456)
+++ trunk/NEMO/OPA_SRC/DYN/dynspg_flt.F90	(revision 474)
@@ -286,8 +286,6 @@
             CALL sol_sor_e( kindic )
          ELSE                          ! e r r o r in nsolv namelist parameter
-            IF(lwp) WRITE(numout,cform_err)
-            IF(lwp) WRITE(numout,*) ' dyn_spg_flt : e r r o r, nsolv = 1, 2 ,3 or 4'
-            IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~                not = ', nsolv
-            nstop = nstop + 1
+            WRITE(ctmp1,*) ' ~~~~~~~~~~~                not = ', nsolv
+            CALL ctl_stop( ' dyn_spg_flt : e r r o r, nsolv = 1, 2 ,3 or 4', ctmp1 )
          ENDIF
       ENDIF
Index: trunk/NEMO/OPA_SRC/DYN/dynspg_flt_jki.F90
===================================================================
--- trunk/NEMO/OPA_SRC/DYN/dynspg_flt_jki.F90	(revision 456)
+++ trunk/NEMO/OPA_SRC/DYN/dynspg_flt_jki.F90	(revision 474)
@@ -300,9 +300,6 @@
             CALL sol_sor_e( kindic )
          ELSE                          ! e r r o r in nsolv namelist parameter
-            IF(lwp) WRITE(numout,cform_err)
-            IF(lwp) WRITE(numout,*) ' dyn_spg_flt_jki : e r r o r, nsolv = 1, 2, 3 or 4'
-            IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~~~~~                not = ', nsolv
-            nstop = nstop + 1
-         ENDIF
+            WRITE(ctmp1,*) ' ~~~~~~~~~~~~~~~~                not = ', nsolv
+            CALL ctl_stop( ' dyn_spg_flt_jki : e r r o r, nsolv = 1, 2, 3 or 4', ctmp1 )
       ENDIF
 
Index: trunk/NEMO/OPA_SRC/DYN/dynspg_rl.F90
===================================================================
--- trunk/NEMO/OPA_SRC/DYN/dynspg_rl.F90	(revision 456)
+++ trunk/NEMO/OPA_SRC/DYN/dynspg_rl.F90	(revision 474)
@@ -229,8 +229,6 @@
             CALL sol_sor_e( kindic )
          CASE DEFAULT                  ! e r r o r in nsolv namelist parameter
-            IF(lwp) WRITE(numout,cform_err)
-            IF(lwp) WRITE(numout,*) ' dyn_spg_rl : e r r o r, nsolv = 1, 2 ,3 or 4'
-            IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~                not = ', nsolv
-            nstop = nstop + 1
+            WRITE(ctmp1,*) ' ~~~~~~~~~~                not = ', nsolv
+            CALL ctl_stop( ' dyn_spg_rl : e r r o r, nsolv = 1, 2 ,3 or 4', ctmp1 )
          END SELECT
       ENDIF
Index: trunk/NEMO/OPA_SRC/DYN/dynvor.F90
===================================================================
--- trunk/NEMO/OPA_SRC/DYN/dynvor.F90	(revision 456)
+++ trunk/NEMO/OPA_SRC/DYN/dynvor.F90	(revision 474)
@@ -677,9 +677,5 @@
          ioptio = ioptio + 1
       ENDIF
-      IF( ioptio /= 1 .AND. .NOT. lk_esopa ) THEN
-         WRITE(numout,cform_err)
-         IF(lwp) WRITE(numout,*) ' use ONE and ONLY one vorticity scheme'
-         nstop = nstop + 1
-      ENDIF
+      IF( ioptio /= 1 .AND. .NOT. lk_esopa ) CALL ctl_stop( ' use ONE and ONLY one vorticity scheme' )
       IF( lk_esopa ) THEN
          nvor = -1
