Index: trunk/NEMO/LIM_SRC/icestp.F90
===================================================================
--- trunk/NEMO/LIM_SRC/icestp.F90	(revision 257)
+++ trunk/NEMO/LIM_SRC/icestp.F90	(revision 258)
@@ -31,4 +31,5 @@
    USE limwri
    USE limrst
+   USE prtctl          ! Print control
 
    IMPLICIT NONE
@@ -186,15 +187,15 @@
          
 
-         IF(l_ctl) THEN         ! print mean trends (used for debugging)
-            WRITE(numout,*) 'Ice Forcings '
-            WRITE(numout,*) ' qsr_oce  : ', SUM( qsr_oce (2:nictl,2:njctl) ), ' qsr_ice  : ', SUM( qsr_ice (2:nictl,2:njctl) )
-            WRITE(numout,*) ' qnsr_oce : ', SUM( qnsr_oce(2:nictl,2:njctl) ), ' qnsr_ice : ', SUM( qnsr_ice(2:nictl,2:njctl) )
-            WRITE(numout,*) ' evap     : ', SUM( evap    (2:nictl,2:njctl) )
-            WRITE(numout,*) ' precip   : ', SUM( tprecip(2:nictl,2:njctl)  ), ' Snow     : ', SUM( sprecip (2:nictl,2:njctl) )
-            WRITE(numout,*) ' u-stress : ', SUM( gtaux  (2:nictl,2:njctl)  ), ' v-stress : ', SUM( gtauy  (2:nictl,2:njctl) )
-            WRITE(numout,*) ' sst      : ', SUM( sst_io (2:nictl,2:njctl)  ), ' sss      : ', SUM( sss_io (2:nictl,2:njctl) )
-            WRITE(numout,*) ' u_io     : ', SUM( u_io   (2:nictl,2:njctl)  ), ' v_io     : ', SUM( v_io   (2:nictl,2:njctl) )
-            WRITE(numout,*) ' hsnif  1 : ', SUM( hsnif  (2:nictl,2:njctl)  ), ' hicnif   : ', SUM( hicif  (2:nictl,2:njctl) )
-            WRITE(numout,*) ' frld   1 : ', SUM( frld   (2:nictl,2:njctl)  ), ' sist     : ', SUM( sist   (2:nictl,2:njctl) )
+         IF(ln_ctl) THEN         ! print mean trends (used for debugging)
+            CALL prt_ctl_info('Ice Forcings ')
+            CALL prt_ctl(tab2d_1=qsr_oce ,clinfo1=' qsr_oce  : ', tab2d_2=qsr_ice , clinfo2=' qsr_ice   : ')
+            CALL prt_ctl(tab2d_1=qnsr_oce,clinfo1=' qnsr_oce : ', tab2d_2=qnsr_ice, clinfo2=' qnsr_ice  : ')
+            CALL prt_ctl(tab2d_1=evap    ,clinfo1=' evap     : ')
+            CALL prt_ctl(tab2d_1=tprecip ,clinfo1=' precip   : ', tab2d_2=sprecip , clinfo2=' Snow      : ')
+            CALL prt_ctl(tab2d_1=gtaux   ,clinfo1=' u-stress : ', tab2d_2=gtauy   , clinfo2=' v-stress  : ')
+            CALL prt_ctl(tab2d_1=sst_io  ,clinfo1=' sst      : ', tab2d_2=sss_io  , clinfo2=' sss       : ')
+            CALL prt_ctl(tab2d_1=u_io    ,clinfo1=' u_io     : ', tab2d_2=v_io    , clinfo2=' v_io      : ')
+            CALL prt_ctl(tab2d_1=hsnif   ,clinfo1=' hsnif  1 : ', tab2d_2=hicif   , clinfo2=' hicif     : ')
+            CALL prt_ctl(tab2d_1=frld    ,clinfo1=' frld   1 : ', tab2d_2=sist    , clinfo2=' sist      : ')
          ENDIF
 
@@ -205,7 +206,7 @@
          CALL lim_dyn                                                ! Ice dynamics !   ( rheology/dynamics )
          !                                                           !--------------!
-         IF(l_ctl) THEN
-            WRITE(numout,*) ' hsnif  2 : ', SUM( hsnif  (2:nictl,2:njctl) ), ' hicnif   : ', SUM( hicif  (2:nictl,2:njctl) )
-            WRITE(numout,*) ' frld   2 : ', SUM( frld   (2:nictl,2:njctl) ), ' sist     : ', SUM( sist   (2:nictl,2:njctl) )
+         IF(ln_ctl) THEN
+            CALL prt_ctl(tab2d_1=hsnif ,clinfo1=' hsnif  2 : ', tab2d_2=hicif , clinfo2=' hicif     : ')
+            CALL prt_ctl(tab2d_1=frld  ,clinfo1=' frld   2 : ', tab2d_2=sist  , clinfo2=' sist      : ')
          ENDIF
 
@@ -214,7 +215,7 @@
          CALL lim_trp                                                ! Ice transport !  ( Advection/diffusion )
          !                                                           !---------------!
-         IF(l_ctl) THEN
-            WRITE(numout,*) ' hsnif  3 : ', SUM( hsnif  (2:nictl,2:njctl) ), ' hicnif   : ', SUM( hicif  (2:nictl,2:njctl) )
-            WRITE(numout,*) ' frld   3 : ', SUM( frld   (2:nictl,2:njctl) ), ' sist     : ', SUM( sist   (2:nictl,2:njctl) )
+         IF(ln_ctl) THEN
+            CALL prt_ctl(tab2d_1=hsnif ,clinfo1=' hsnif  3 : ', tab2d_2=hicif , clinfo2=' hicif     : ')
+            CALL prt_ctl(tab2d_1=frld  ,clinfo1=' frld   3 : ', tab2d_2=sist  , clinfo2=' sist      : ')
          ENDIF
 
@@ -223,7 +224,7 @@
          CALL lim_thd                                                ! Ice thermodynamics !
          !                                                           !--------------------!
-         IF(l_ctl) THEN
-            WRITE(numout,*) ' hsnif  4 : ', SUM( hsnif  (2:nictl,2:njctl) ), ' hicnif   : ', SUM( hicif  (2:nictl,2:njctl) )
-            WRITE(numout,*) ' frld   4 : ', SUM( frld   (2:nictl,2:njctl) ), ' sist     : ', SUM( sist   (2:nictl,2:njctl) )
+         IF(ln_ctl) THEN
+            CALL prt_ctl(tab2d_1=hsnif ,clinfo1=' hsnif  4 : ', tab2d_2=hicif , clinfo2=' hicif     : ')
+            CALL prt_ctl(tab2d_1=frld  ,clinfo1=' frld   4 : ', tab2d_2=sist  , clinfo2=' sist      : ')
          ENDIF
 
Index: trunk/NEMO/LIM_SRC/limadv.F90
===================================================================
--- trunk/NEMO/LIM_SRC/limadv.F90	(revision 257)
+++ trunk/NEMO/LIM_SRC/limadv.F90	(revision 258)
@@ -17,4 +17,5 @@
    USE in_out_manager  ! I/O manager
    USE lbclnk
+   USE prtctl          ! Print control
 
    IMPLICIT NONE
@@ -224,9 +225,9 @@
       CALL lbc_lnk( psxy, 'T', 1. )
 
-      IF(l_ctl) THEN
-         WRITE(numout,*) ' lim_adv_x: psm  ', SUM( psm (2:nictl,2:njctl) ), ' ps0  ', SUM( ps0 (2:nictl,2:njctl) )
-         WRITE(numout,*) ' lim_adv_x: psx  ', SUM( psx (2:nictl,2:njctl) ), ' psxx ', SUM( psxx(2:nictl,2:njctl) )
-         WRITE(numout,*) ' lim_adv_x: psy  ', SUM( psy (2:nictl,2:njctl) ), ' psyy ', SUM( psyy(2:nictl,2:njctl) )
-         WRITE(numout,*) ' lim_adv_x: psxy ', SUM( psxy(2:nictl,2:njctl) )
+      IF(ln_ctl)   THEN
+         CALL prt_ctl(tab2d_1=psm  , clinfo1=' lim_adv_x: psm  :', tab2d_2=ps0 , clinfo2=' ps0  : ')
+         CALL prt_ctl(tab2d_1=psx  , clinfo1=' lim_adv_x: psx  :', tab2d_2=psxx, clinfo2=' psxx : ')
+         CALL prt_ctl(tab2d_1=psy  , clinfo1=' lim_adv_x: psy  :', tab2d_2=psyy, clinfo2=' psyy : ')
+         CALL prt_ctl(tab2d_1=psxy , clinfo1=' lim_adv_x: psxy :')
       ENDIF
 
@@ -424,9 +425,9 @@
       CALL lbc_lnk( psxy, 'T', 1. )
 
-      IF(l_ctl) THEN
-         WRITE(numout,*) ' lim_adv_y: psm  ', SUM( psm (2:nictl,2:njctl) ), ' ps0  ', SUM( ps0 (2:nictl,2:njctl) )
-         WRITE(numout,*) ' lim_adv_y: psx  ', SUM( psx (2:nictl,2:njctl) ), ' psxx ', SUM( psxx(2:nictl,2:njctl) )
-         WRITE(numout,*) ' lim_adv_y: psy  ', SUM( psy (2:nictl,2:njctl) ), ' psyy ', SUM( psyy(2:nictl,2:njctl) )
-         WRITE(numout,*) ' lim_adv_y: psxy ', SUM( psxy(2:nictl,2:njctl) )
+      IF(ln_ctl) THEN
+         CALL prt_ctl(tab2d_1=psm  , clinfo1=' lim_adv_y: psm  :', tab2d_2=ps0 , clinfo2=' ps0  : ')
+         CALL prt_ctl(tab2d_1=psx  , clinfo1=' lim_adv_y: psx  :', tab2d_2=psxx, clinfo2=' psxx : ')
+         CALL prt_ctl(tab2d_1=psy  , clinfo1=' lim_adv_y: psy  :', tab2d_2=psyy, clinfo2=' psyy : ')
+         CALL prt_ctl(tab2d_1=psxy , clinfo1=' lim_adv_y: psxy :')
       ENDIF
 
Index: trunk/NEMO/LIM_SRC/limdyn.F90
===================================================================
--- trunk/NEMO/LIM_SRC/limdyn.F90	(revision 257)
+++ trunk/NEMO/LIM_SRC/limdyn.F90	(revision 258)
@@ -23,4 +23,5 @@
    USE lbclnk
    USE lib_mpp
+   USE prtctl          ! Print control
 
    IMPLICIT NONE
@@ -92,5 +93,7 @@
             i_j1 = 1   
             i_jpj = jpj
-            IF(l_ctl)   WRITE(numout,*) 'lim_dyn :    i_j1 = ', i_j1, ' ij_jpj = ', i_jpj
+            IF(ln_ctl)    THEN
+               CALL prt_ctl_info('lim_dyn  :    i_j1 = ', ivar1=i_j1, clinfo2=' ij_jpj = ', ivar2=i_jpj)
+            ENDIF
             CALL lim_rhg( i_j1, i_jpj )
 
@@ -113,5 +116,5 @@
                END DO
                i_j1 = MAX( 1, i_j1-1 )
-               IF(l_ctl)   WRITE(numout,*) 'lim_dyn : NH i_j1 = ', i_j1, ' ij_jpj = ', i_jpj
+               IF(ln_ctl)   WRITE(numout,*) 'lim_dyn : NH i_j1 = ', i_j1, ' ij_jpj = ', i_jpj
     
                CALL lim_rhg( i_j1, i_jpj )
@@ -124,5 +127,5 @@
                END DO
                i_jpj = MIN( jpj, i_jpj+2 )
-               IF(l_ctl)   WRITE(numout,*) 'lim_dyn : SH i_j1 = ', i_j1, ' ij_jpj = ', i_jpj
+               IF(ln_ctl)   WRITE(numout,*) 'lim_dyn : SH i_j1 = ', i_j1, ' ij_jpj = ', i_jpj
     
                CALL lim_rhg( i_j1, i_jpj )
@@ -143,5 +146,5 @@
                i_jpj = MIN( jpj, i_jpj+2)
     
-               IF(l_ctl)   WRITE(numout,*) 'lim_dyn : one hemisphere: i_j1 = ', i_j1, ' ij_jpj = ', i_jpj
+               IF(ln_ctl)   WRITE(numout,*) 'lim_dyn : one hemisphere: i_j1 = ', i_j1, ' ij_jpj = ', i_jpj
     
                CALL lim_rhg( i_j1, i_jpj )
@@ -154,7 +157,7 @@
          v_ice(:,1) = 0.e0
 
-         IF(l_ctl) THEN 
-            WRITE(numout,*) ' lim_dyn  : u_oce ', SUM( u_oce(2:nictl,2:njctl) ), ' v_oce ', SUM( v_oce(2:nictl,2:njctl) )
-            WRITE(numout,*) ' lim_dyn  : u_ice ', SUM( u_ice(2:nictl,2:njctl) ), ' v_ice ', SUM( v_ice(2:nictl,2:njctl) )
+         IF(ln_ctl)   THEN 
+            CALL prt_ctl(tab2d_1=u_oce , clinfo1=' lim_dyn  : u_oce :', tab2d_2=v_oce , clinfo2=' v_oce :')
+            CALL prt_ctl(tab2d_1=u_ice , clinfo1=' lim_dyn  : u_ice :', tab2d_2=v_ice , clinfo2=' v_ice :')
          ENDIF
          
@@ -244,7 +247,7 @@
       CALL lbc_lnk( tio_v, 'I', -1. )   ! I-point (i.e. ice U-V point)
 
-      IF(l_ctl) THEN 
-         WRITE(numout,*) ' lim_dyn  : tio_u ', SUM( tio_u(2:nictl,2:njctl) ), ' tio_v ', SUM( tio_v(2:nictl,2:njctl) )
-         WRITE(numout,*) ' lim_dyn  : ust2s ', SUM( ust2s(2:nictl,2:njctl) )
+      IF(ln_ctl) THEN 
+            CALL prt_ctl(tab2d_1=tio_u , clinfo1=' lim_dyn  : tio_u :', tab2d_2=tio_v , clinfo2=' tio_v :')
+            CALL prt_ctl(tab2d_1=ust2s , clinfo1=' lim_dyn  : ust2s :')
       ENDIF
 
Index: trunk/NEMO/LIM_SRC/limflx.F90
===================================================================
--- trunk/NEMO/LIM_SRC/limflx.F90	(revision 257)
+++ trunk/NEMO/LIM_SRC/limflx.F90	(revision 258)
@@ -11,14 +11,15 @@
    !!----------------------------------------------------------------------
    !! * Modules used
-   USE par_oce
-   USE phycst
-   USE ocfzpt
-   USE ice_oce
-   USE flx_oce
-   USE ice
-   USE flxblk
-   USE lbclnk
-   USE in_out_manager
-   USE albedo
+   USE par_oce          ! ocean parameters
+   USE phycst           ! physical constants
+   USE ocfzpt           ! surface ocean freezing point
+   USE ice_oce          ! sea-ice variable
+   USE flx_oce          ! sea-ice/ocean forcings variables
+   USE ice              ! LIM sea-ice variables
+   USE flxblk           ! bulk formulea
+   USE lbclnk           ! ocean lateral boundary condition
+   USE in_out_manager   ! I/O manager
+   USE albedo           ! albedo parameters
+   USE prtctl           ! Print control
 
    IMPLICIT NONE
@@ -195,10 +196,9 @@
 #endif
 
-      IF(l_ctl) THEN
-         WRITE(numout,*) ' lim_flx  '
-         WRITE(numout,*) ' fsolar ', SUM( fsolar(2:nictl,2:njctl) ), ' fnsolar', SUM( fnsolar(2:nictl,2:njctl) )
-         WRITE(numout,*) ' fmass  ', SUM( fmass (2:nictl,2:njctl) ), ' fsalt  ', SUM( fsalt  (2:nictl,2:njctl) )
-         WRITE(numout,*) ' ftaux  ', SUM( ftaux (2:nictl,2:njctl) ), ' ftauy  ', SUM( ftauy  (2:nictl,2:njctl) )
-         WRITE(numout,*) ' freeze ', SUM( freeze(2:nictl,2:njctl) ), ' tn_ice ', SUM( tn_ice (2:nictl,2:njctl) )
+      IF(ln_ctl) THEN
+         CALL prt_ctl(tab2d_1=fsolar, clinfo1=' lim_flx: fsolar : ', tab2d_2=fnsolar, clinfo2=' fnsolar : ')
+         CALL prt_ctl(tab2d_1=fmass , clinfo1=' lim_flx: fmass  : ', tab2d_2=fsalt  , clinfo2=' fsalt   : ')
+         CALL prt_ctl(tab2d_1=ftaux , clinfo1=' lim_flx: ftaux  : ', tab2d_2=ftauy  , clinfo2=' ftauy   : ')
+         CALL prt_ctl(tab2d_1=freeze, clinfo1=' lim_flx: freeze : ', tab2d_2=tn_ice , clinfo2=' tn_ice  : ')
       ENDIF 
    
Index: trunk/NEMO/LIM_SRC/limhdf.F90
===================================================================
--- trunk/NEMO/LIM_SRC/limhdf.F90	(revision 257)
+++ trunk/NEMO/LIM_SRC/limhdf.F90	(revision 258)
@@ -18,4 +18,5 @@
    USE lbclnk
    USE lib_mpp
+   USE prtctl          ! Print control
 
    IMPLICIT NONE
@@ -67,4 +68,5 @@
       INTEGER ::  &
          its, iter            ! temporary integers
+      CHARACTER (len=55) :: charout
       REAL(wp) ::  &
          zalfa, zrlxint, zconv, zeps   ! temporary scalars
@@ -168,7 +170,9 @@
 
       ptab(:,:) = ptab(:,:)
-
-      IF(l_ctl)   WRITE(numout,*) ' lim_hdf  : ', SUM( ptab(2:nictl,2:njctl)-ptab0(2:nictl,2:njctl) ), &
-      &                                        ' zconv= ', zconv, ' iter= ', iter
+      IF(ln_ctl)   THEN
+         zrlx(:,:) = ptab(:,:) - ptab0(:,:)
+         WRITE(charout,FMT="(' lim_hdf  : zconv =',D23.16, ' iter =',I4,2X)") zconv, iter
+         CALL prt_ctl(tab2d_1=zrlx, clinfo1=charout)
+      ENDIF
 
    END SUBROUTINE lim_hdf
Index: trunk/NEMO/LIM_SRC/limrhg.F90
===================================================================
--- trunk/NEMO/LIM_SRC/limrhg.F90	(revision 257)
+++ trunk/NEMO/LIM_SRC/limrhg.F90	(revision 258)
@@ -19,4 +19,5 @@
    USE lib_mpp
    USE in_out_manager  ! I/O manager
+   USE prtctl          ! Print control
 
    IMPLICIT NONE
@@ -64,4 +65,6 @@
          iim1, ijm1, iip1 , ijp1   , & ! temporary integers
          iter, jter                    !    "          "
+
+      CHARACTER (len=50) :: charout
 
       REAL(wp) :: &
@@ -632,7 +635,8 @@
       !                                                   ! ==================== !
 
-      IF(l_ctl) THEN
-         WRITE(numout,*) ' lim_rhg  : res= ', zresm, 'iter= ', jter,' u_ice ', SUM( u_ice(2:nictl,2:njctl) ) , &
-         &                                                          ' v_ice ', SUM( v_ice(2:nictl,2:njctl) )
+      IF(ln_ctl) THEN
+         WRITE(charout,FMT="('lim_rhg  : res =',D23.16, ' iter =',I4)") zresm, jter
+         CALL prt_ctl_info(charout)
+         CALL prt_ctl(tab2d_1=u_ice, clinfo1=' lim_rhg  : u_ice :', tab2d_2=v_ice, clinfo2=' v_ice :')
       ENDIF
 
Index: trunk/NEMO/LIM_SRC/limthd.F90
===================================================================
--- trunk/NEMO/LIM_SRC/limthd.F90	(revision 257)
+++ trunk/NEMO/LIM_SRC/limthd.F90	(revision 258)
@@ -25,4 +25,5 @@
    USE limthd_lac
    USE limtab
+   USE prtctl          ! Print control
       
    IMPLICIT NONE
@@ -78,4 +79,5 @@
          nbpb  ,               &   ! nb of icy pts for thermo. cal.
          nbpac                     ! nb of pts for lateral accretion 
+      CHARACTER (len=22) :: charout
       REAL(wp) ::  &
          zfric_umin = 5e-03 ,  &   ! lower bound for the friction velocity
@@ -92,4 +94,6 @@
          zhicifp            ,  &   ! ice thickness for outputs
          zqlbsbq                   ! link with lead energy budget qldif
+      REAL(wp), DIMENSION(jpi,jpj,2) :: &
+         zmsk                      ! working array
       !!-------------------------------------------------------------------
 
@@ -112,4 +116,5 @@
       rdmicif(:,:) = 0.e0   ! variation of ice mass per unit area
       hicifp (:,:) = 0.e0   ! daily thermodynamic ice production. 
+      zmsk (:,:,:) = 0.e0
 
       DO jj = 1, jpj
@@ -118,6 +123,6 @@
          END DO
       END DO
-      IF(l_ctl)   WRITE(numout,*) 'lim_thd  : ', SUM( hsnif(2:nictl,2:njctl) ) , ' hsnif'
-      
+
+      IF(ln_ctl)   CALL prt_ctl(tab2d_1=hsnif     , clinfo1=' lim_thd: hsnif   : ')
       
       !-----------------------------------!
@@ -150,10 +155,11 @@
          END DO
       END DO
-      IF(l_ctl) THEN
-         WRITE(numout,*) 'lim_thd: hicif : ', SUM( hicif (2:nictl,2:njctl) )
-         WRITE(numout,*) 'lim_thd: hsnif : ', SUM( hsnif (2:nictl,2:njctl) )
-         WRITE(numout,*) 'lim_thd: dmgwi : ', SUM( dmgwi (2:nictl,2:njctl) )
-         WRITE(numout,*) 'lim_thd: qstoif: ', SUM( qstoif(2:nictl,2:njctl) )
-         WRITE(numout,*) 'lim_thd: frld  : ', SUM( frld  (2:nictl,2:njctl) )
+
+      IF(ln_ctl) THEN
+         CALL prt_ctl(tab2d_1=hicif  , clinfo1=' lim_thd: hicif   : ')
+         CALL prt_ctl(tab2d_1=hsnif  , clinfo1=' lim_thd: hsnif   : ')
+         CALL prt_ctl(tab2d_1=dmgwi  , clinfo1=' lim_thd: dmgwi   : ')
+         CALL prt_ctl(tab2d_1=qstoif , clinfo1=' lim_thd: qstoif  : ')
+         CALL prt_ctl(tab2d_1=frld   , clinfo1=' lim_thd: frld    : ')
       ENDIF
 
@@ -221,12 +227,15 @@
          END DO
       END DO
-      IF(l_ctl) THEN
-         WRITE(numout,*) 'lim_thd: pfrld ' , SUM( pfrld  (2:nictl,2:njctl) ), ' thcm    ', SUM( thcm   (2:nictl,2:njctl) )
-         WRITE(numout,*) 'lim_thd: fdtcn ' , SUM( fdtcn  (2:nictl,2:njctl) ), ' qdtcn   ', SUM( qdtcn  (2:nictl,2:njctl) )
-         WRITE(numout,*) 'lim_thd: qldif ' , SUM( qldif  (2:nictl,2:njctl) ), ' zqlbsbq ', SUM( zqlbsbq(2:nictl,2:njctl) )
-         WRITE(numout,*) 'lim_thd: qcmif ' , SUM( qcmif  (2:nictl,2:njctl) ), ' fbif    ', SUM( fbif   (2:nictl,2:njctl) )
-         WRITE(numout,*) 'lim_thd: qcmif  ', SUM( qcmif  (2:nictl,2:njctl)*tms(2:nictl,2:njctl) )
-         WRITE(numout,*) 'lim_thd: zhicifp', SUM( zhicifp(2:nictl,2:njctl) )
-         WRITE(numout,*) 'limthd : nbpb = ', nbpb
+
+      IF(ln_ctl) THEN
+         CALL prt_ctl(tab2d_1=pfrld, clinfo1=' lim_thd: pfrld   : ', tab2d_2=thcm   , clinfo2='  thcm    : ')
+         CALL prt_ctl(tab2d_1=fdtcn, clinfo1=' lim_thd: fdtcn   : ', tab2d_2=qdtcn  , clinfo2='  qdtcn   : ')
+         CALL prt_ctl(tab2d_1=qldif, clinfo1=' lim_thd: qldif   : ', tab2d_2=zqlbsbq, clinfo2='  zqlbsbq : ')
+         CALL prt_ctl(tab2d_1=qcmif, clinfo1=' lim_thd: qcmif   : ', tab2d_2=fbif   , clinfo2='  fbif    : ')
+         zmsk(:,:,1) = tms(:,:)
+         CALL prt_ctl(tab2d_1=qcmif  , clinfo1=' lim_thd: qcmif   : ', mask1=zmsk)
+         CALL prt_ctl(tab2d_1=zhicifp, clinfo1=' lim_thd: zhicifp : ')
+         WRITE(charout, FMT="('lim_thd: nbpb = ',I4)") nbpb
+         CALL prt_ctl_info(charout)
       ENDIF
       
@@ -327,7 +336,8 @@
       END DO
       
-      IF(l_ctl) THEN
-         WRITE(numout,*) 'lim_thd : phicif ', SUM( phicif(2:nictl,2:njctl) ), ' hicif ', SUM( hicif(2:nictl,2:njctl) )
-         WRITE(numout,*) 'lim_thd : nbpac = ', nbpac
+      IF(ln_ctl) THEN
+         CALL prt_ctl(tab2d_1=phicif, clinfo1=' lim_thd: phicif  : ', tab2d_2=hicif, clinfo2=' hicif : ')
+         WRITE(charout, FMT="('lim_thd: nbpac = ',I4)") nbpac
+         CALL prt_ctl_info(charout)
       ENDIF
 
@@ -381,13 +391,15 @@
       END DO
 
-      IF(l_ctl) THEN
-         WRITE(numout,*) ' lim_thd  end  '
-         WRITE(numout,*) '  hicif ', SUM( hicif (2:nictl,2:njctl)   ), '  hsnif ', SUM( hsnif (2:nictl,2:njctl)   )
-         WRITE(numout,*) '  frld  ', SUM( frld  (2:nictl,2:njctl)   ), '  hicifp', SUM( hicifp(2:nictl,2:njctl)   )
-         WRITE(numout,*) '  phicif', SUM( phicif(2:nictl,2:njctl)   ), '  pfrld ', SUM( pfrld (2:nictl,2:njctl)   )
-         WRITE(numout,*) '  sist  ', SUM( sist  (2:nictl,2:njctl)   ), '  tbif 1', SUM( tbif  (2:nictl,2:njctl,1) )
-         WRITE(numout,*) '  tbif 2', SUM( tbif  (2:nictl,2:njctl,2) ), '  tbif 3', SUM( tbif  (2:nictl,2:njctl,3) )
-         WRITE(numout,*) '  fdtcn ', SUM( fdtcn (2:nictl,2:njctl)   ), '  qdtcn ', SUM( qdtcn (2:nictl,2:njctl)   )
-         WRITE(numout,*) '  qstoif', SUM( qstoif(2:nictl,2:njctl)   ), '  fsbbq ', SUM( fsbbq (2:nictl,2:njctl)   )
+      IF(ln_ctl) THEN
+         CALL prt_ctl_info(' lim_thd  end  ')
+         CALL prt_ctl(tab2d_1=hicif , clinfo1=' lim_thd: hicif   : ', tab2d_2=hsnif , clinfo2=' hsnif  : ')
+         CALL prt_ctl(tab2d_1=frld  , clinfo1=' lim_thd: frld    : ', tab2d_2=hicifp, clinfo2=' hicifp : ')
+         CALL prt_ctl(tab2d_1=phicif, clinfo1=' lim_thd: phicif  : ', tab2d_2=pfrld , clinfo2=' pfrld  : ')
+         CALL prt_ctl(tab2d_1=sist  , clinfo1=' lim_thd: sist    : ')
+         CALL prt_ctl(tab2d_1=tbif(:,:,1), clinfo1=' lim_thd: tbif 1  : ')
+         CALL prt_ctl(tab2d_1=tbif(:,:,2), clinfo1=' lim_thd: tbif 2  : ')
+         CALL prt_ctl(tab2d_1=tbif(:,:,3), clinfo1=' lim_thd: tbif 3  : ')
+         CALL prt_ctl(tab2d_1=fdtcn , clinfo1=' lim_thd: fdtcn   : ', tab2d_2=qdtcn , clinfo2=' qdtcn  : ')
+         CALL prt_ctl(tab2d_1=qstoif, clinfo1=' lim_thd: qstoif  : ', tab2d_2=fsbbq , clinfo2=' fsbbq  : ')
       ENDIF
 
Index: trunk/NEMO/LIM_SRC/limwri.F90
===================================================================
--- trunk/NEMO/LIM_SRC/limwri.F90	(revision 257)
+++ trunk/NEMO/LIM_SRC/limwri.F90	(revision 258)
@@ -136,5 +136,5 @@
       !-----------------------------------------------------------------------
 
-      IF(lwp .AND. L_ctl) THEN
+      IF(lwp) THEN
          WRITE(numout,*)
          WRITE(numout,*) 'lim_wri : write ice outputs in NetCDF files at time : ', nyear, nmonth, nday, numit
Index: trunk/NEMO/OPA_SRC/DOM/dom_oce.F90
===================================================================
--- trunk/NEMO/OPA_SRC/DOM/dom_oce.F90	(revision 257)
+++ trunk/NEMO/OPA_SRC/DOM/dom_oce.F90	(revision 258)
@@ -66,4 +66,5 @@
    INTEGER, PUBLIC, DIMENSION(jpnij) ::   &  !:
       nimppt, njmppt,  &  !: i-, j-indexes for each processor
+      ibonit, ibonjt,  &  !: i-, j- processor neighbour existence
       nlcit, nlcjt,    &  !: dimensions of every subdomain
       nldit, nldjt,    &  !: first, last indoor index for each i-domain
Index: trunk/NEMO/OPA_SRC/DOM/domain.F90
===================================================================
--- trunk/NEMO/OPA_SRC/DOM/domain.F90	(revision 257)
+++ trunk/NEMO/OPA_SRC/DOM/domain.F90	(revision 258)
@@ -156,8 +156,9 @@
       !! * Modules used
       USE ioipsl
-      NAMELIST/namrun/ no    , cexper   , ln_rstart , nrstdt , nit000,         &
-         &             nitend, ndate0   , nleapy   , ninist , nstock,          &
-         &             nprint, nwrite   , nrunoff  , ln_ctl , nictl , njctl,   &
-         &             nbench
+      NAMELIST/namrun/ no    , cexper   , ln_rstart , nrstdt , nit000,          &
+         &             nitend, ndate0   , nleapy   , ninist , nstock,           &
+         &             nprint, nwrite   , nrunoff  , ln_ctl , nictls, nictle,   &
+         &             njctls, njctle   , nbench   , isplt  , jsplt
+
       NAMELIST/namdom/ ntopo , e3zps_min, e3zps_rat, ngrid  , nmsh  ,   &
          &             nacc  , atfp     , rdt      , rdtmin , rdtmax,   &
@@ -192,13 +193,76 @@
          WRITE(numout,*) '           runoff option                   nrunoff   = ', nrunoff
          WRITE(numout,*) '           run control (for debugging)     ln_ctl    = ', ln_ctl
-         WRITE(numout,*) '           Max i indice for SUM control    nictl     = ', nictl
-         WRITE(numout,*) '           Max j indice for SUM control    njctl     = ', njctl
+         WRITE(numout,*) '           Start i indice for SUM control  nictls    = ', nictls
+         WRITE(numout,*) '           End i indice for SUM control    nictle    = ', nictle
+         WRITE(numout,*) '           Start j indice for SUM control  njctls    = ', njctls
+         WRITE(numout,*) '           End j indice for SUM control    njctle    = ', njctle
+         WRITE(numout,*) '           number of proc. following i     isplt     = ', isplt
+         WRITE(numout,*) '           number of proc. following j     jsplt     = ', jsplt
          WRITE(numout,*) '           benchmark parameter (0/1)       nbench    = ', nbench
       ENDIF
 
-      l_ctl = ln_ctl .AND. lwp       ! trend control print on the 1st processor only
-
       ndastp = ndate0                ! Assign initial date to current date
 
+! ... Control the sub-domain area indices for the print control
+      IF(ln_ctl)   THEN
+         IF( lk_mpp ) THEN
+            ! the domain is forced to the real splitted domain in MPI
+            isplt = jpni ; jsplt = jpnj ; ijsplt = jpni*jpnj
+         ELSE
+            IF( isplt == 1 .AND. jsplt == 1  ) THEN
+               IF(lwp) WRITE(numout,cform_war)
+               IF(lwp) WRITE(numout,*)'          - isplt & jsplt are equal to 1'
+               IF(lwp) WRITE(numout,*)'          - the print control will be done over the whole domain'
+               IF(lwp) WRITE(numout,*)
+            ENDIF
+
+            ! compute the total number of processors ijsplt
+            ijsplt = isplt*jsplt
+         ENDIF
+
+         IF(lwp) WRITE(numout,*)'          - The total number of processors over which the'
+         IF(lwp) WRITE(numout,*)'            print control will be done is ijsplt : ', ijsplt
+
+         ! Control the indices used for the SUM control
+         IF( nictls+nictle+njctls+njctle == 0 )   THEN
+            ! the print control is done over the default area
+            lsp_area = .FALSE.
+         ELSE
+            ! the print control is done over a specific  area
+            lsp_area = .TRUE.
+            IF( nictls < 1 .OR. nictls > jpiglo )   THEN
+               IF(lwp) WRITE(numout,cform_war)
+               IF(lwp) WRITE(numout,*)'          - nictls must be 1<=nictls>=jpiglo, it is forced to 1'
+               IF(lwp) WRITE(numout,*)
+               nwarn = nwarn + 1
+               nictls = 1
+            ENDIF
+
+            IF( nictle < 1 .OR. nictle > jpiglo )   THEN
+               IF(lwp) WRITE(numout,cform_war)
+               IF(lwp) WRITE(numout,*)'          - nictle must be 1<=nictle>=jpiglo, it is forced to jpiglo'
+               IF(lwp) WRITE(numout,*)
+               nwarn = nwarn + 1
+               nictle = jpjglo
+            ENDIF
+
+            IF( njctls < 1 .OR. njctls > jpjglo )   THEN
+               IF(lwp) WRITE(numout,cform_war)
+               IF(lwp) WRITE(numout,*)'          - njctls must be 1<=njctls>=jpjglo, it is forced to 1'
+               IF(lwp) WRITE(numout,*)
+               nwarn = nwarn + 1
+               njctls = 1
+            ENDIF
+
+            IF( njctle < 1 .OR. njctle > jpjglo )   THEN
+               IF(lwp) WRITE(numout,cform_war)
+               IF(lwp) WRITE(numout,*)'          - njctle must be 1<=njctle>= jpjglo, it is forced to jpjglo'
+               IF(lwp) WRITE(numout,*)
+               nwarn = nwarn + 1
+               njctle = jpjglo
+            ENDIF
+
+         ENDIF          ! IF( nictls+nictle+njctls+njctle == 0 )
+       ENDIF            ! IF(ln_ctl)
 
 ! ... Control of output frequency
@@ -247,25 +311,4 @@
          WRITE(numout,*) '           nb of seconds per year   raass = ', raass, ' s'
          WRITE(numout,*) '           nb of seconds per month  rmoss = ', rmoss, ' s'
-      ENDIF
-
-! ... Control the Max i and j indices used for the SUM control (i.e. when ln_ctl=.true.)
-      IF ( nictl > jpim1 ) THEN
-          IF(lwp) THEN
-             WRITE(numout,cform_war)
-             WRITE(numout,*) '           nictl = ', nictl, ' must be <= to jpim1 '
-             WRITE(numout,*) '           nictl forced to be equal to jpim1 '
-          ENDIF
-          nwarn = nwarn + 1
-          nictl = jpim1
-      ENDIF
-
-      IF ( njctl > jpjm1 ) THEN
-          IF(lwp) THEN
-             WRITE(numout,cform_war)
-             WRITE(numout,*) '           njctl = ', njctl, ' must be <= to jpjm1 '
-             WRITE(numout,*) '           njctl forced to be equal to jpjm1 '
-          ENDIF
-          nwarn = nwarn + 1
-          njctl = jpjm1
       ENDIF
 
Index: trunk/NEMO/OPA_SRC/DYN/dynhpg.F90
===================================================================
--- trunk/NEMO/OPA_SRC/DYN/dynhpg.F90	(revision 257)
+++ trunk/NEMO/OPA_SRC/DYN/dynhpg.F90	(revision 258)
@@ -21,4 +21,5 @@
    USE trdmod          ! ocean dynamics trends 
    USE trdmod_oce      ! ocean variables trends
+   USE prtctl          ! Print control
 
    IMPLICIT NONE
@@ -172,9 +173,7 @@
       ENDIF
 
-      IF(l_ctl) THEN         ! print sum trends (used for debugging)
-         zuap = SUM( ua(2:nictl,2:njctl,1:jpkm1) * umask(2:nictl,2:njctl,1:jpkm1) )
-         zvap = SUM( va(2:nictl,2:njctl,1:jpkm1) * vmask(2:nictl,2:njctl,1:jpkm1) )
-         WRITE(numout,*) ' hpg  - Ua: ', zuap-u_ctl, ' Va: ', zvap-v_ctl
-         u_ctl = zuap   ;   v_ctl = zvap
+      IF(ln_ctl) THEN         ! print sum trends (used for debugging)
+         CALL prt_ctl(tab3d_1=ua, clinfo1=' hpg  - Ua: ', mask1=umask, &
+            &         tab3d_2=va, clinfo2=' Va: ', mask2=vmask, clinfo3='dyn')
       ENDIF
 
@@ -325,9 +324,7 @@
       ENDIF
 
-      IF(l_ctl) THEN         ! print sum trends (used for debugging)
-         zuap = SUM( ua(2:nictl,2:njctl,1:jpkm1) * umask(2:nictl,2:njctl,1:jpkm1) )
-         zvap = SUM( va(2:nictl,2:njctl,1:jpkm1) * vmask(2:nictl,2:njctl,1:jpkm1) )
-         WRITE(numout,*) ' hpg  - Ua: ', zuap-u_ctl, ' Va: ', zvap-v_ctl
-         u_ctl = zuap   ;   v_ctl = zvap
+      IF(ln_ctl) THEN         ! print sum trends (used for debugging)
+         CALL prt_ctl(tab3d_1=ua, clinfo1=' hpg  - Ua: ', mask1=umask, &
+            &         tab3d_2=va, clinfo2=' Va: ', mask2=vmask, clinfo3='dyn')
       ENDIF
 
@@ -443,9 +440,7 @@
       ENDIF
 
-      IF(l_ctl) THEN         ! print sum trends (used for debugging)
-         zuap = SUM( ua(2:nictl,2:njctl,1:jpkm1) * umask(2:nictl,2:njctl,1:jpkm1) )
-         zvap = SUM( va(2:nictl,2:njctl,1:jpkm1) * vmask(2:nictl,2:njctl,1:jpkm1) )
-         WRITE(numout,*) ' hpg  - Ua: ', zuap-u_ctl, ' Va: ', zvap-v_ctl
-         u_ctl = zuap   ;   v_ctl = zvap
+      IF(ln_ctl) THEN         ! print sum trends (used for debugging)
+         CALL prt_ctl(tab3d_1=ua, clinfo1=' hpg  - Ua: ', mask1=umask, &
+            &         tab3d_2=va, clinfo2=' Va: ', mask2=vmask, clinfo3='dyn')
       ENDIF
 
Index: trunk/NEMO/OPA_SRC/DYN/dynhpg_atsk.F90
===================================================================
--- trunk/NEMO/OPA_SRC/DYN/dynhpg_atsk.F90	(revision 257)
+++ trunk/NEMO/OPA_SRC/DYN/dynhpg_atsk.F90	(revision 258)
@@ -21,4 +21,5 @@
    USE trdmod          ! ocean dynamics trends 
    USE trdmod_oce      ! ocean variables trends
+   USE prtctl          ! Print control
 
    IMPLICIT NONE
@@ -162,9 +163,7 @@
       ENDIF
 
-      IF(l_ctl) THEN         ! print sum trends (used for debugging)
-         zuap = SUM( ua(2:nictl,2:njctl,1:jpkm1) * umask(2:nictl,2:njctl,1:jpkm1) )
-         zvap = SUM( va(2:nictl,2:njctl,1:jpkm1) * vmask(2:nictl,2:njctl,1:jpkm1) )
-         WRITE(numout,*) ' hpg  - Ua: ', zuap-u_ctl, ' Va: ', zvap-v_ctl
-         u_ctl = zuap   ;   v_ctl = zvap
+      IF(ln_ctl) THEN         ! print sum trends (used for debugging)
+         CALL prt_ctl(tab3d_1=ua, clinfo1=' hpg  - Ua: ', mask1=umask, &
+            &         tab3d_2=va, clinfo2=' Va: ', mask2=vmask, clinfo3='dyn')
       ENDIF
 
@@ -306,9 +305,7 @@
       ENDIF
 
-      IF(l_ctl) THEN         ! print sum trends (used for debugging)
-         zuap = SUM( ua(2:nictl,2:njctl,1:jpkm1) * umask(2:nictl,2:njctl,1:jpkm1) )
-         zvap = SUM( va(2:nictl,2:njctl,1:jpkm1) * vmask(2:nictl,2:njctl,1:jpkm1) )
-         WRITE(numout,*) ' hpg  - Ua: ', zuap-u_ctl, ' Va: ', zvap-v_ctl
-         u_ctl = zuap   ;   v_ctl = zvap
+      IF(ln_ctl) THEN         ! print sum trends (used for debugging)
+         CALL prt_ctl(tab3d_1=ua, clinfo1=' hpg  - Ua: ', mask1=umask, &
+            &         tab3d_2=va, clinfo2=' Va: ', mask2=vmask, clinfo3='dyn')
       ENDIF   
 
@@ -429,9 +426,7 @@
       ENDIF
 
-      IF(l_ctl) THEN         ! print sum trends (used for debugging)
-         zuap = SUM( ua(2:nictl,2:njctl,1:jpkm1) * umask(2:nictl,2:njctl,1:jpkm1) )
-         zvap = SUM( va(2:nictl,2:njctl,1:jpkm1) * vmask(2:nictl,2:njctl,1:jpkm1) )
-         WRITE(numout,*) ' hpg  - Ua: ', zuap-u_ctl, ' Va: ', zvap-v_ctl
-         u_ctl = zuap   ;   v_ctl = zvap
+      IF(ln_ctl) THEN         ! print sum trends (used for debugging)
+         CALL prt_ctl(tab3d_1=ua, clinfo1=' hpg  - Ua: ', mask1=umask, &
+            &         tab3d_2=va, clinfo2=' Va: ', mask2=vmask, clinfo3='dyn')
       ENDIF
 
Index: trunk/NEMO/OPA_SRC/DYN/dynkeg.F90
===================================================================
--- trunk/NEMO/OPA_SRC/DYN/dynkeg.F90	(revision 257)
+++ trunk/NEMO/OPA_SRC/DYN/dynkeg.F90	(revision 258)
@@ -14,4 +14,5 @@
    USE trdmod          ! ocean dynamics trends 
    USE trdmod_oce      ! ocean variables trends
+   USE prtctl          ! Print control
 
    IMPLICIT NONE
@@ -119,9 +120,7 @@
       ENDIF
 
-      IF(l_ctl) THEN         ! print sum trends (used for debugging)
-         zua = SUM( ua(2:nictl,2:njctl,1:jpkm1) * umask(2:nictl,2:njctl,1:jpkm1) )
-         zva = SUM( va(2:nictl,2:njctl,1:jpkm1) * vmask(2:nictl,2:njctl,1:jpkm1) )
-         WRITE(numout,*) ' keg  - Ua: ', zua-u_ctl, ' Va: ', zva-v_ctl
-         u_ctl = zua   ;   v_ctl = zva
+      IF(ln_ctl) THEN         ! print sum trends (used for debugging)
+         CALL prt_ctl(tab3d_1=ua, clinfo1=' keg  - Ua: ', mask1=umask, &
+            &         tab3d_2=va, clinfo2=' Va: ', mask2=vmask, clinfo3='dyn')
       ENDIF
 
Index: trunk/NEMO/OPA_SRC/DYN/dynldf_bilap.F90
===================================================================
--- trunk/NEMO/OPA_SRC/DYN/dynldf_bilap.F90	(revision 257)
+++ trunk/NEMO/OPA_SRC/DYN/dynldf_bilap.F90	(revision 258)
@@ -17,4 +17,5 @@
    USE trdmod_oce      ! ocean variables trends
    USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
+   USE prtctl          ! Print control
 
    IMPLICIT NONE
@@ -237,9 +238,7 @@
       ENDIF
 
-      IF(l_ctl) THEN         ! print sum trends (used for debugging)
-         zua = SUM( ua(2:nictl,2:njctl,1:jpkm1) * umask(2:nictl,2:njctl,1:jpkm1) )
-         zva = SUM( va(2:nictl,2:njctl,1:jpkm1) * vmask(2:nictl,2:njctl,1:jpkm1) )
-         WRITE(numout,*) ' ldf  - Ua: ', zua-u_ctl, ' Va: ', zva-v_ctl
-         u_ctl = zua   ;   v_ctl = zva
+      IF(ln_ctl) THEN         ! print sum trends (used for debugging)
+         CALL prt_ctl(tab3d_1=ua, clinfo1=' ldf  - Ua: ', mask1=umask, &
+            &         tab3d_2=va, clinfo2=' Va: ', mask2=vmask, clinfo3='dyn')
       ENDIF
 
Index: trunk/NEMO/OPA_SRC/DYN/dynldf_bilapg.F90
===================================================================
--- trunk/NEMO/OPA_SRC/DYN/dynldf_bilapg.F90	(revision 257)
+++ trunk/NEMO/OPA_SRC/DYN/dynldf_bilapg.F90	(revision 258)
@@ -22,4 +22,5 @@
    USE ldfslp          ! iso-neutral slopes available
    USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
+   USE prtctl          ! Print control
 
    IMPLICIT NONE
@@ -81,5 +82,4 @@
       !! * Local declarations
       INTEGER ::   ji, jj, jk                 ! dummy loop indices
-      REAL(wp) ::   zua, zva                  ! temporary scalars
       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   &
          zwk1, zwk2                ! work array used for rotated biharmonic
@@ -138,9 +138,7 @@
       ENDIF
 
-      IF(l_ctl) THEN         ! print sum trends (used for debugging)
-         zua = SUM( ua(2:nictl,2:njctl,1:jpkm1) * umask(2:nictl,2:njctl,1:jpkm1) )
-         zva = SUM( va(2:nictl,2:njctl,1:jpkm1) * vmask(2:nictl,2:njctl,1:jpkm1) )
-         WRITE(numout,*) ' ldf  - Ua: ', zua-u_ctl, ' Va: ', zva-v_ctl
-         u_ctl = zua   ;   v_ctl = zva
+      IF(ln_ctl) THEN         ! print sum trends (used for debugging)
+         CALL prt_ctl(tab3d_1=ua, clinfo1=' ldf  - Ua: ', mask1=umask, &
+            &         tab3d_2=va, clinfo2=' Va: ', mask2=vmask, clinfo3='dyn')
       ENDIF
 
Index: trunk/NEMO/OPA_SRC/DYN/dynldf_iso.F90
===================================================================
--- trunk/NEMO/OPA_SRC/DYN/dynldf_iso.F90	(revision 257)
+++ trunk/NEMO/OPA_SRC/DYN/dynldf_iso.F90	(revision 258)
@@ -22,4 +22,5 @@
    USE ldfslp          ! iso-neutral slopes 
    USE in_out_manager  ! I/O manager
+   USE prtctl          ! Print control
 
    IMPLICIT NONE
@@ -98,5 +99,5 @@
          zabe1, zabe2, zcof1, zcof2,   &  ! temporary scalars
          zmskt, zmskf, zbu, zbv,       &
-         zuah, zvah, zua, zva
+         zuah, zvah
       REAL(wp), DIMENSION(jpi,jpj) ::   &
          ziut, zjuf, zjvt, zivf,        & ! temporary workspace
@@ -259,9 +260,7 @@
       ENDIF
 
-      IF(l_ctl) THEN         ! print sum trends (used for debugging)
-         zua = SUM( ua(2:nictl,2:njctl,1:jpkm1) * umask(2:nictl,2:njctl,1:jpkm1) )
-         zva = SUM( va(2:nictl,2:njctl,1:jpkm1) * vmask(2:nictl,2:njctl,1:jpkm1) )
-         WRITE(numout,*) ' ldf  - Ua: ', zua-u_ctl, ' Va: ', zva-v_ctl
-         u_ctl = zua   ;   v_ctl = zva
+      IF(ln_ctl) THEN         ! print sum trends (used for debugging)
+         CALL prt_ctl(tab3d_1=ua, clinfo1=' ldf  - Ua: ', mask1=umask, &
+            &         tab3d_2=va, clinfo2=' Va: ', mask2=vmask, clinfo3='dyn')
       ENDIF
 
Index: trunk/NEMO/OPA_SRC/DYN/dynldf_lap.F90
===================================================================
--- trunk/NEMO/OPA_SRC/DYN/dynldf_lap.F90	(revision 257)
+++ trunk/NEMO/OPA_SRC/DYN/dynldf_lap.F90	(revision 258)
@@ -18,4 +18,5 @@
    USE trdmod_oce      ! ocean variables trends
    USE ldfslp          ! iso-neutral slopes 
+   USE prtctl          ! Print control
 
    IMPLICIT NONE
@@ -136,9 +137,7 @@
       ENDIF
 
-      IF(l_ctl) THEN         ! print sum trends (used for debugging)
-         zua = SUM( ua(2:nictl,2:njctl,1:jpkm1) * umask(2:nictl,2:njctl,1:jpkm1) )
-         zva = SUM( va(2:nictl,2:njctl,1:jpkm1) * vmask(2:nictl,2:njctl,1:jpkm1) )
-         WRITE(numout,*) ' ldf  - Ua: ', zua-u_ctl, ' Va: ', zva-v_ctl
-         u_ctl = zua   ;   v_ctl = zva
+      IF(ln_ctl) THEN         ! print sum trends (used for debugging)
+         CALL prt_ctl(tab3d_1=ua, clinfo1=' ldf  - Ua: ', mask1=umask, &
+            &         tab3d_2=va, clinfo2=' Va: ', mask2=vmask, clinfo3='dyn')
       ENDIF
 
Index: trunk/NEMO/OPA_SRC/DYN/dynnxt.F90
===================================================================
--- trunk/NEMO/OPA_SRC/DYN/dynnxt.F90	(revision 257)
+++ trunk/NEMO/OPA_SRC/DYN/dynnxt.F90	(revision 258)
@@ -14,4 +14,5 @@
    USE obcdyn          ! open boundary condition for momentum (obc_dyn routine)
    USE lbclnk          ! lateral boundary condition (or mpp link)
+   USE prtctl          ! Print control
 
    IMPLICIT NONE
@@ -135,6 +136,8 @@
       !                                                ! ===============
 
-      IF(l_ctl)   WRITE(numout,*) ' nxt  - Un: ', SUM(un(2:nictl,2:njctl,1:jpkm1)*umask(2:nictl,2:njctl,1:jpkm1)), &
-      &                                  ' Vn: ', SUM(vn(2:nictl,2:njctl,1:jpkm1)*vmask(2:nictl,2:njctl,1:jpkm1))
+      IF(ln_ctl)   THEN
+         CALL prt_ctl(tab3d_1=un, clinfo1=' nxt  - Un: ', mask1=umask, &
+            &         tab3d_2=vn, clinfo2=' Vn: ', mask2=vmask)
+      ENDIF
 
    END SUBROUTINE dyn_nxt
Index: trunk/NEMO/OPA_SRC/DYN/dynspg_fsc.F90
===================================================================
--- trunk/NEMO/OPA_SRC/DYN/dynspg_fsc.F90	(revision 257)
+++ trunk/NEMO/OPA_SRC/DYN/dynspg_fsc.F90	(revision 258)
@@ -33,4 +33,5 @@
    USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
    USE cla_dynspg      ! cross land advection
+   USE prtctl          ! Print control
 
    IMPLICIT NONE
@@ -337,7 +338,7 @@
       ENDIF
 
-      IF(l_ctl) THEN         ! print sum trends (used for debugging)
-         WRITE(numout,*) ' spg  - Ua: ', SUM( ua(2:nictl,2:njctl,1:jpkm1)*umask(2:nictl,2:njctl,1:jpkm1) ),   &
-            &                   ' Va: ', SUM( va(2:nictl,2:njctl,1:jpkm1)*vmask(2:nictl,2:njctl,1:jpkm1) )
+      IF(ln_ctl) THEN         ! print sum trends (used for debugging)
+         CALL prt_ctl(tab3d_1=ua, clinfo1=' spg  - Ua: ', mask1=umask, &
+            &         tab3d_2=va, clinfo2=' Va: ', mask2=vmask)
       ENDIF
 
@@ -370,6 +371,6 @@
 
 
-      IF(l_ctl) THEN         ! print sum trends (used for debugging)
-         WRITE(numout,*) ' spg - ssh:', SUM( sshn(2:nictl,2:njctl)*tmask(2:nictl,2:njctl,1) )
+      IF(ln_ctl) THEN         ! print sum trends (used for debugging)
+         CALL prt_ctl(tab2d_1=sshn, clinfo1=' spg  - ssh: ', mask1=tmask)
       ENDIF
 
Index: trunk/NEMO/OPA_SRC/DYN/dynspg_fsc_atsk.F90
===================================================================
--- trunk/NEMO/OPA_SRC/DYN/dynspg_fsc_atsk.F90	(revision 257)
+++ trunk/NEMO/OPA_SRC/DYN/dynspg_fsc_atsk.F90	(revision 258)
@@ -37,4 +37,5 @@
    USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
    USE cla_dynspg      ! cross land advection
+   USE prtctl          ! Print control
 
    IMPLICIT NONE
@@ -354,7 +355,7 @@
          END DO
 
-         IF(l_ctl) THEN         ! print sum trends (used for debugging)
-            WRITE(numout,*) ' spg  - Ua: ', SUM( ua(2:nictl,2:njctl,1:jpkm1)*umask(2:nictl,2:njctl,1:jpkm1) ),   &
-               &                   ' Va: ', SUM( va(2:nictl,2:njctl,1:jpkm1)*vmask(2:nictl,2:njctl,1:jpkm1) )
+         IF(ln_ctl) THEN         ! print sum trends (used for debugging)
+            CALL prt_ctl(tab3d_1=ua, clinfo1=' spg  - Ua: ', mask1=umask, &
+               &         tab3d_2=va, clinfo2=' Va: ', mask2=vmask)
          ENDIF
 
@@ -398,6 +399,6 @@
       CALL lbc_lnk( sshn, 'T', 1. )
 
-      IF(l_ctl) THEN         ! print sum trends (used for debugging)
-         WRITE(numout,*) ' spg - ssh:', SUM( sshn(2:nictl,2:njctl)*tmask(2:nictl,2:njctl,1) )
+      IF(ln_ctl) THEN         ! print sum trends (used for debugging)
+         CALL prt_ctl(tab2d_1=sshn, clinfo1=' spg  - ssh: ', mask1=tmask)
       ENDIF
 
Index: trunk/NEMO/OPA_SRC/DYN/dynvor.F90
===================================================================
--- trunk/NEMO/OPA_SRC/DYN/dynvor.F90	(revision 257)
+++ trunk/NEMO/OPA_SRC/DYN/dynvor.F90	(revision 258)
@@ -20,4 +20,5 @@
    USE trdmod_oce      ! ocean variables trends
    USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
+   USE prtctl          ! Print control
 
    IMPLICIT NONE
@@ -188,9 +189,7 @@
       ENDIF
 
-      IF(l_ctl) THEN         ! print sum trends (used for debugging)
-         zua = SUM( ua(2:nictl,2:njctl,1:jpkm1) * umask(2:nictl,2:njctl,1:jpkm1) )
-         zva = SUM( va(2:nictl,2:njctl,1:jpkm1) * vmask(2:nictl,2:njctl,1:jpkm1) )
-         WRITE(numout,*) ' vor  - Ua: ', zua-u_ctl, ' Va: ', zva-v_ctl
-         u_ctl = zua   ;   v_ctl = zva
+      IF(ln_ctl) THEN         ! print sum trends (used for debugging)
+         CALL prt_ctl(tab3d_1=ua, clinfo1=' vor  - Ua: ', mask1=umask, &
+            &         tab3d_2=va, clinfo2=' Va: ', mask2=vmask, clinfo3='dyn')
       ENDIF
 
@@ -342,9 +341,7 @@
       ENDIF
 
-      IF(l_ctl) THEN         ! print sum trends (used for debugging)
-         zua = SUM( ua(2:nictl,2:njctl,1:jpkm1) * umask(2:nictl,2:njctl,1:jpkm1) )
-         zva = SUM( va(2:nictl,2:njctl,1:jpkm1) * vmask(2:nictl,2:njctl,1:jpkm1) )
-         WRITE(numout,*) ' vor  - Ua: ', zua-u_ctl, ' Va: ', zva-v_ctl
-         u_ctl = zua   ;   v_ctl = zva
+      IF(ln_ctl) THEN         ! print sum trends (used for debugging)
+         CALL prt_ctl(tab3d_1=ua, clinfo1=' vor  - Ua: ', mask1=umask, &
+            &         tab3d_2=va, clinfo2=' Va: ', mask2=vmask, clinfo3='dyn')
       ENDIF
 
@@ -499,9 +496,7 @@
       ENDIF
 
-      IF(l_ctl) THEN         ! print sum trends (used for debugging)
-         zua = SUM( ua(2:nictl,2:njctl,1:jpkm1) * umask(2:nictl,2:njctl,1:jpkm1) )
-         zva = SUM( va(2:nictl,2:njctl,1:jpkm1) * vmask(2:nictl,2:njctl,1:jpkm1) )
-         WRITE(numout,*) ' vor  - Ua: ', zua-u_ctl, ' Va: ', zva-v_ctl
-         u_ctl = zua   ;   v_ctl = zva
+      IF(ln_ctl) THEN         ! print sum trends (used for debugging)
+         CALL prt_ctl(tab3d_1=ua, clinfo1=' vor  - Ua: ', mask1=umask, &
+            &         tab3d_2=va, clinfo2=' Va: ', mask2=vmask, clinfo3='dyn')
       ENDIF
 
@@ -662,9 +657,7 @@
       ENDIF
 
-      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,*) ' vor  een - Ua: ', zua-u_ctl, ' Va: ', zva-v_ctl
-         u_ctl = zua   ;   v_ctl = zva
+      IF(ln_ctl) THEN         ! print sum trends (used for debugging)
+         CALL prt_ctl(tab3d_1=ua, clinfo1=' vor  - Ua: ', mask1=umask, &
+            &         tab3d_2=va, clinfo2=' Va: ', mask2=vmask, clinfo3='dyn')
       ENDIF
 
Index: trunk/NEMO/OPA_SRC/DYN/dynzad.F90
===================================================================
--- trunk/NEMO/OPA_SRC/DYN/dynzad.F90	(revision 257)
+++ trunk/NEMO/OPA_SRC/DYN/dynzad.F90	(revision 258)
@@ -15,4 +15,5 @@
    USE trdmod_oce      ! ocean variables trends
    USE flxrnf          ! ocean runoffs
+   USE prtctl          ! Print control
 
    IMPLICIT NONE
@@ -142,9 +143,7 @@
       ENDIF
 
-      IF(l_ctl) THEN         ! print sum trends (used for debugging)
-         zua = SUM( ua(2:nictl,2:njctl,1:jpkm1) * umask(2:nictl,2:njctl,1:jpkm1) )
-         zva = SUM( va(2:nictl,2:njctl,1:jpkm1) * vmask(2:nictl,2:njctl,1:jpkm1) )
-         WRITE(numout,*) ' zad  - Ua: ', zua-u_ctl, ' Va: ', zva-v_ctl
-         u_ctl = zua   ;   v_ctl = zva
+      IF(ln_ctl) THEN         ! print sum trends (used for debugging)
+         CALL prt_ctl(tab3d_1=ua, clinfo1=' zad  - Ua: ', mask1=umask, &
+            &         tab3d_2=va, clinfo2=' Va: ', mask2=vmask, clinfo3='dyn')
       ENDIF
 
@@ -255,9 +254,7 @@
       ENDIF
 
-      IF(l_ctl) THEN         ! print sum trends (used for debugging)
-         zua = SUM( ua(2:nictl,2:njctl,1:jpkm1) * umask(2:nictl,2:njctl,1:jpkm1) )
-         zva = SUM( va(2:nictl,2:njctl,1:jpkm1) * vmask(2:nictl,2:njctl,1:jpkm1) )
-         WRITE(numout,*) ' zad  - Ua: ', zua-u_ctl, ' Va: ', zva-v_ctl
-         u_ctl = zua   ;   v_ctl = zva
+      IF(ln_ctl) THEN         ! print sum trends (used for debugging)
+         CALL prt_ctl(tab3d_1=ua, clinfo1=' zad  - Ua: ', mask1=umask, &
+            &         tab3d_2=va, clinfo2=' Va: ', mask2=vmask, clinfo3='dyn')
       ENDIF
 
Index: trunk/NEMO/OPA_SRC/DYN/dynzdf_exp.F90
===================================================================
--- trunk/NEMO/OPA_SRC/DYN/dynzdf_exp.F90	(revision 257)
+++ trunk/NEMO/OPA_SRC/DYN/dynzdf_exp.F90	(revision 258)
@@ -18,4 +18,5 @@
    USE trdmod          ! ocean dynamics trends 
    USE trdmod_oce      ! ocean variables trends
+   USE prtctl          ! Print control
 
    IMPLICIT NONE
@@ -197,9 +198,7 @@
       ENDIF
 
-      IF(l_ctl) THEN         ! print sum trends (used for debugging)
-         zua = SUM( ua(2:nictl,2:njctl,1:jpkm1) * umask(2:nictl,2:njctl,1:jpkm1) )
-         zva = SUM( va(2:nictl,2:njctl,1:jpkm1) * vmask(2:nictl,2:njctl,1:jpkm1) )
-         WRITE(numout,*) ' zdf  - Ua: ', zua-u_ctl, ' Va: ', zva-v_ctl
-         u_ctl = zua   ;   v_ctl = zva
+      IF(ln_ctl) THEN         ! print sum trends (used for debugging)
+         CALL prt_ctl(tab3d_1=ua, clinfo1=' zdf  - Ua: ', mask1=umask, &
+            &         tab3d_2=va, clinfo2=' Va: ', mask2=vmask, clinfo3='dyn')
       ENDIF
 
Index: trunk/NEMO/OPA_SRC/DYN/dynzdf_imp.F90
===================================================================
--- trunk/NEMO/OPA_SRC/DYN/dynzdf_imp.F90	(revision 257)
+++ trunk/NEMO/OPA_SRC/DYN/dynzdf_imp.F90	(revision 258)
@@ -22,4 +22,5 @@
    USE trdmod          ! ocean dynamics trends 
    USE trdmod_oce      ! ocean variables trends
+   USE prtctl          ! Print control
 
    IMPLICIT NONE
@@ -79,5 +80,5 @@
          ikbu, ikbum1, ikbv, ikbvm1      ! temporary integers
       REAL(wp) ::   &
-         zrau0r, z2dt, zua, zva,      &  ! temporary scalars
+         zrau0r, z2dt,                &  ! temporary scalars
          z2dtf, zcoef, zzws, zrhs        !    "         "
       REAL(wp), DIMENSION(jpi,jpj) ::   &
@@ -355,9 +356,7 @@
       ENDIF
 
-      IF(l_ctl) THEN         ! print sum trends (used for debugging)
-         zua = SUM( ua(2:nictl,2:njctl,1:jpkm1) * umask(2:nictl,2:njctl,1:jpkm1) )
-         zva = SUM( va(2:nictl,2:njctl,1:jpkm1) * vmask(2:nictl,2:njctl,1:jpkm1) )
-         WRITE(numout,*) ' zdf  - Ua: ', zua-u_ctl, ' Va: ', zva-v_ctl
-         u_ctl = zua   ;   v_ctl = zva
+      IF(ln_ctl) THEN         ! print sum trends (used for debugging)
+         CALL prt_ctl(tab3d_1=ua, clinfo1=' zdf  - Ua: ', mask1=umask, &
+            &         tab3d_2=va, clinfo2=' Va: ', mask2=vmask, clinfo3='dyn')
       ENDIF
 
Index: trunk/NEMO/OPA_SRC/DYN/dynzdf_imp_atsk.F90
===================================================================
--- trunk/NEMO/OPA_SRC/DYN/dynzdf_imp_atsk.F90	(revision 257)
+++ trunk/NEMO/OPA_SRC/DYN/dynzdf_imp_atsk.F90	(revision 258)
@@ -19,4 +19,5 @@
    USE trdmod          ! ocean dynamics trends 
    USE trdmod_oce      ! ocean variables trends
+   USE prtctl          ! Print control
 
    IMPLICIT NONE
@@ -74,5 +75,5 @@
          ikbu, ikbum1, ikbv, ikbvm1      !    "         "      
       REAL(wp) ::   &
-         zrau0r, z2dt, zua, zva,       & !temporary scalars
+         zrau0r, z2dt,                 & !temporary scalars
          z2dtf, zcoef, zzws
       REAL(wp), DIMENSION(jpi,jpk) ::  &
@@ -338,9 +339,7 @@
       ENDIF
 
-      IF(l_ctl) THEN         ! print sum trends (used for debugging)
-         zua = SUM( ua(2:nictl,2:njctl,1:jpkm1) * umask(2:nictl,2:njctl,1:jpkm1) )
-         zva = SUM( va(2:nictl,2:njctl,1:jpkm1) * vmask(2:nictl,2:njctl,1:jpkm1) )
-         WRITE(numout,*) ' zdf  - Ua: ', zua-u_ctl, ' Va: ', zva-v_ctl
-         u_ctl = zua   ;   v_ctl = zva
+      IF(ln_ctl) THEN         ! print sum trends (used for debugging)
+         CALL prt_ctl(tab3d_1=ua, clinfo1=' zdf  - Ua: ', mask1=umask, &
+            &         tab3d_2=va, clinfo2=' Va: ', mask2=vmask, clinfo3='dyn')
       ENDIF
 
Index: trunk/NEMO/OPA_SRC/DYN/dynzdf_iso.F90
===================================================================
--- trunk/NEMO/OPA_SRC/DYN/dynzdf_iso.F90	(revision 257)
+++ trunk/NEMO/OPA_SRC/DYN/dynzdf_iso.F90	(revision 258)
@@ -21,4 +21,5 @@
    USE trdmod          ! ocean dynamics trends 
    USE trdmod_oce      ! ocean variables trends
+   USE prtctl          ! Print control
 
    IMPLICIT NONE
@@ -110,5 +111,5 @@
       REAL(wp) ::   &
          zrau0r, z2dt,                       & ! temporary scalars
-         z2dtf, zua, zva, zcoef, zzws
+         z2dtf, zcoef, zzws
       REAL(wp) ::   &
          zcoef0, zcoef3, zcoef4, zbu, zbv, zmkt, zmkf,   &
@@ -406,9 +407,7 @@
       ENDIF
 
-      IF(l_ctl) THEN         ! print sum trends (used for debugging)
-         zua = SUM( ua(2:nictl,2:njctl,1:jpkm1) * umask(2:nictl,2:njctl,1:jpkm1) )
-         zva = SUM( va(2:nictl,2:njctl,1:jpkm1) * vmask(2:nictl,2:njctl,1:jpkm1) )
-         WRITE(numout,*) ' zdf  - Ua: ', zua-u_ctl, ' Va: ', zva-v_ctl
-         u_ctl = zua   ;   v_ctl = zva
+      IF(ln_ctl) THEN         ! print sum trends (used for debugging)
+         CALL prt_ctl(tab3d_1=ua, clinfo1=' zdf  - Ua: ', mask1=umask, &
+            &         tab3d_2=va, clinfo2=' Va: ', mask2=vmask, clinfo3='dyn')
       ENDIF
 
Index: trunk/NEMO/OPA_SRC/DYN/wzvmod.F90
===================================================================
--- trunk/NEMO/OPA_SRC/DYN/wzvmod.F90	(revision 257)
+++ trunk/NEMO/OPA_SRC/DYN/wzvmod.F90	(revision 258)
@@ -12,4 +12,5 @@
    USE dom_oce         ! ocean space and time domain variables 
    USE in_out_manager  ! I/O manager
+   USE prtctl          ! Print control
 
    IMPLICIT NONE
@@ -80,7 +81,5 @@
       !                                                ! ===============
 
-      IF(l_ctl) THEN                ! print mean trends (used for debugging)
-         WRITE(numout,*) ' w**2 -   : ', SUM( wn(2:nictl,2:njctl,1:jpkm1) * wn   (2:nictl,2:njctl,1:jpkm1) )
-      ENDIF
+      IF(ln_ctl)   CALL prt_ctl(tab3d_1=wn, clinfo1=' w**2 -   : ', mask1=wn)
 
    END SUBROUTINE wzv
@@ -131,7 +130,5 @@
       END DO
 
-      IF(l_ctl) THEN                ! print mean trends (used for debugging)
-         WRITE(numout,*) ' w**2 -   : ', SUM( wn(2:nictl,2:njctl,1:jpkm1) * wn   (2:nictl,2:njctl,1:jpkm1) )
-      ENDIF
+      IF(ln_ctl)   CALL prt_ctl(tab3d_1=wn, clinfo1=' w**2 -   : ', mask1=wn)
 
    END SUBROUTINE wzv
Index: trunk/NEMO/OPA_SRC/LDF/ldfeiv.F90
===================================================================
--- trunk/NEMO/OPA_SRC/LDF/ldfeiv.F90	(revision 257)
+++ trunk/NEMO/OPA_SRC/LDF/ldfeiv.F90	(revision 258)
@@ -22,4 +22,5 @@
    USE in_out_manager  ! I/O manager
    USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
+   USE prtctl          ! Print control
 
    IMPLICIT NONE
@@ -162,7 +163,9 @@
       CALL lbc_lnk( aeiv, 'V', 1. )
 
-      IF(l_ctl)   WRITE(numout,*) ' eiv  - u: ', SUM( aeiu(1:nictl+1,1:njctl+1) ), ' v: ', SUM( aeiv(1:nictl+1,1:njctl+1) )
-
-
+      IF(ln_ctl)   THEN
+         CALL prt_ctl(tab2d_1=aeiu, clinfo1=' eiv  - u: ', ovlap=1)
+         CALL prt_ctl(tab2d_1=aeiv, clinfo1=' eiv  - v: ', ovlap=1)
+      ENDIF
+      
       ! ORCA R05: add a space variation on aht (=aeiv except at the equator and river mouth)
       IF( cp_cfg == "orca" .AND. jp_cfg == 05 ) THEN
@@ -178,8 +181,8 @@
             END DO
          END DO
-         IF(l_ctl) THEN
-            WRITE(numout,*)' aht  - u: ', SUM( ahtu(1:nictl+1,1:njctl+1) )
-            WRITE(numout,*)'      - v: ', SUM( ahtv(1:nictl+1,1:njctl+1) )
-            WRITE(numout,*)'      - w: ', SUM( ahtw(1:nictl+1,1:njctl+1) )
+         IF(ln_ctl) THEN
+            CALL prt_ctl(tab2d_1=ahtu, clinfo1=' aht  - u: ', ovlap=1)
+            CALL prt_ctl(tab2d_1=ahtv, clinfo1=' aht  - v: ', ovlap=1)
+            CALL prt_ctl(tab2d_1=ahtw, clinfo1=' aht  - w: ', ovlap=1)
          ENDIF
       ENDIF
@@ -327,5 +330,8 @@
       CALL lbc_lnk( aeiv, 'V', 1. )
 
-      IF(l_ctl)   WRITE(numout,*) ' eiv  - u: ', SUM( aeiu(1:nictl+1,1:njctl+1) ), ' v: ', SUM( aeiv(1:nictl+1,1:njctl+1) )
+      IF(ln_ctl)   THEN
+         CALL prt_ctl(tab2d_1=aeiu, clinfo1=' eiv  - u: ', ovlap=1)
+         CALL prt_ctl(tab2d_1=aeiv, clinfo1=' eiv  - v: ', ovlap=1)
+      ENDIF
 
       ! ORCA R05: add a space variation on aht (=aeiv except at the equator and river mouth)
@@ -342,8 +348,8 @@
             END DO
          END DO
-         IF(l_ctl) THEN
-            WRITE(numout,*) ' aht  - u: ', SUM( ahtu(1:nictl+1,1:njctl+1) )
-            WRITE(numout,*) '      - v: ', SUM( ahtv(1:nictl+1,1:njctl+1) )
-            WRITE(numout,*) '      - w: ', SUM( ahtw(1:nictl+1,1:njctl+1) )
+         IF(ln_ctl) THEN
+            CALL prt_ctl(tab2d_1=ahtu, clinfo1=' aht  - u: ', ovlap=1)
+            CALL prt_ctl(tab2d_1=ahtv, clinfo1=' aht  - v: ', ovlap=1)
+            CALL prt_ctl(tab2d_1=ahtw, clinfo1=' aht  - w: ', ovlap=1)
          ENDIF
       ENDIF
Index: trunk/NEMO/OPA_SRC/LDF/ldfslp.F90
===================================================================
--- trunk/NEMO/OPA_SRC/LDF/ldfslp.F90	(revision 257)
+++ trunk/NEMO/OPA_SRC/LDF/ldfslp.F90	(revision 258)
@@ -21,4 +21,5 @@
    USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
    USE in_out_manager  ! I/O manager
+   USE prtctl          ! Print control
 
    IMPLICIT NONE
@@ -480,7 +481,7 @@
       CALL lbc_lnk( wslpj, 'W', -1. )
 
-      IF(l_ctl) THEN
-         WRITE(numout,*) ' slp  - u : ', SUM( uslp (1:nictl,1:njctl,:) ), ' v : ', SUM( vslp (1:nictl,1:njctl,:) )
-         WRITE(numout,*) '      - wi: ', SUM( wslpi(1:nictl,1:njctl,:) ), ' wj: ', SUM( wslpj(1:nictl,1:njctl,:) )
+      IF(ln_ctl) THEN
+         CALL prt_ctl(tab3d_1=uslp , clinfo1=' slp  - u : ', tab3d_2=vslp,  clinfo2=' v : ', kdim=jpk)
+         CALL prt_ctl(tab3d_1=wslpi, clinfo1=' slp  - wi: ', tab3d_2=wslpj, clinfo2=' wj: ', kdim=jpk)
       ENDIF
 
Index: trunk/NEMO/OPA_SRC/SBC/bulk.F90
===================================================================
--- trunk/NEMO/OPA_SRC/SBC/bulk.F90	(revision 257)
+++ trunk/NEMO/OPA_SRC/SBC/bulk.F90	(revision 258)
@@ -21,4 +21,5 @@
    USE phycst          ! physical constants
    USE in_out_manager  ! I/O manager
+   USE prtctl          ! Print control
 
    IMPLICIT NONE
@@ -97,18 +98,17 @@
 
 # if ! defined key_ice_lim
-         IF(l_ctl) THEN         ! print mean trends (used for debugging)
-            WRITE(numout,*) ' Forcings '
-            WRITE(numout,*) ' qsr_oce  : ', SUM( qsr_oce (1:nictl+1,1:njctl+1) * tmask(1:nictl+1,1:njctl+1,1) )
-            WRITE(numout,*) ' qsr_ice  : ', SUM( qsr_ice (1:nictl+1,1:njctl+1) * tmask(1:nictl+1,1:njctl+1,1) )
-            WRITE(numout,*) ' qnsr_oce : ', SUM( qnsr_oce(1:nictl+1,1:njctl+1) * tmask(1:nictl+1,1:njctl+1,1) )
-            WRITE(numout,*) ' qnsr_ice : ', SUM( qnsr_ice(1:nictl+1,1:njctl+1) * tmask(1:nictl+1,1:njctl+1,1) )
-            WRITE(numout,*) ' evap     : ', SUM( evap    (1:nictl+1,1:njctl+1) * tmask(1:nictl+1,1:njctl+1,1) )
-            WRITE(numout,*) ' precip   : ', SUM( tprecip (1:nictl+1,1:njctl+1) * tmask(1:nictl+1,1:njctl+1,1) ) / rday
-            WRITE(numout,*) ' Snow     : ', SUM( sprecip (1:nictl+1,1:njctl+1) * tmask(1:nictl+1,1:njctl+1,1) ) / rday
-            WRITE(numout,*) ' u-stress : ', SUM( taux    (1:nictl+1,1:njctl+1) * umask(1:nictl+1,1:njctl+1,1) )
-            WRITE(numout,*) ' v-stress : ', SUM( tauy    (1:nictl+1,1:njctl+1) * vmask(1:nictl+1,1:njctl+1,1) )
-            WRITE(numout,*) ' sst      : ', SUM( zsst    (1:nictl+1,1:njctl+1) * tmask(1:nictl+1,1:njctl+1,1) )
-            WRITE(numout,*) ' sss      : ', SUM( zsss    (1:nictl+1,1:njctl+1) * tmask(1:nictl+1,1:njctl+1,1) )
-            WRITE(numout,*)
+         IF(ln_ctl) THEN         ! print mean trends (used for debugging)
+            CALL prt_ctl_info(' Forcings ')
+            CALL prt_ctl(tab2d_1=qsr_oce , clinfo1=' qsr_oce   : ', mask1=tmask, ovlap=1)
+            CALL prt_ctl(tab2d_1=qsr_ice , clinfo1=' qsr_ice   : ', mask1=tmask, ovlap=1)
+            CALL prt_ctl(tab2d_1=qnsr_oce, clinfo1=' qnsr_oce  : ', mask1=tmask, ovlap=1)
+            CALL prt_ctl(tab2d_1=qnsr_ice, clinfo1=' qnsr_ice  : ', mask1=tmask, ovlap=1)
+            CALL prt_ctl(tab2d_1=evap    , clinfo1=' evap      : ', mask1=tmask, ovlap=1)
+            CALL prt_ctl(tab2d_1=tprecip , clinfo1=' precip    : ', mask1=tmask, ovlap=1)
+            CALL prt_ctl(tab2d_1=sprecip , clinfo1=' Snow      : ', mask1=tmask, ovlap=1)
+            CALL prt_ctl(tab2d_1=taux    , clinfo1=' u-stress  : ', mask1=tmask, ovlap=1)
+            CALL prt_ctl(tab2d_1=tauy    , clinfo1=' v-stress  : ', mask1=tmask, ovlap=1)
+            CALL prt_ctl(tab2d_1=zsst    , clinfo1=' sst       : ', mask1=tmask, ovlap=1)
+            CALL prt_ctl(tab2d_1=zsss    , clinfo1=' sss       : ', mask1=tmask, ovlap=1)
          ENDIF
 # endif   
Index: trunk/NEMO/OPA_SRC/SBC/flxblk.F90
===================================================================
--- trunk/NEMO/OPA_SRC/SBC/flxblk.F90	(revision 257)
+++ trunk/NEMO/OPA_SRC/SBC/flxblk.F90	(revision 258)
@@ -24,4 +24,5 @@
    USE lbclnk
    USE albedo
+   USE prtctl          ! Print control
 
    IMPLICIT NONE
@@ -319,5 +320,7 @@
 
       iday   = INT( zxday )
-      IF(l_ctl)   WRITE(numout,*) ' declin : iday ', iday, ' nfbulk= ', nfbulk
+
+      IF(ln_ctl) CALL prt_ctl_info('declin : iday ', ivar1=iday, clinfo2=' nfbulk= ', ivar2=nfbulk)
+
       !   computation of the solar declination, his sine and his cosine
       CALL flx_blk_declin( indaet, iday, zdecl )
Index: trunk/NEMO/OPA_SRC/SBC/ocesbc.F90
===================================================================
--- trunk/NEMO/OPA_SRC/SBC/ocesbc.F90	(revision 257)
+++ trunk/NEMO/OPA_SRC/SBC/ocesbc.F90	(revision 258)
@@ -13,18 +13,19 @@
    USE dom_oce        ! ocean space domain variables
    USE cpl_oce        ! coupled ocean-atmosphere variables
-   USE ice_oce
-   USE blk_oce
-   USE flx_oce
+   USE ice_oce        ! sea-ice variable
+   USE blk_oce        ! bulk variables
+   USE flx_oce        ! sea-ice/ocean forcings variables
    USE phycst         ! Define parameters for the routines
-   USE taumod
-   USE flxmod
-   USE flxrnf
+   USE taumod         ! surface stress forcing
+   USE flxmod         ! thermohaline fluxes
+   USE flxrnf         ! runoffs forcing
    USE tradmp         ! damping salinity trend
-   USE dtatem
-   USE dtasal
-   USE ocfzpt
-   USE lbclnk
-   USE lib_mpp
+   USE dtatem         ! ocean temperature data
+   USE dtasal         ! ocean salinity data
+   USE ocfzpt         ! surface ocean freezing point
+   USE lbclnk         ! ocean lateral boundary condition
+   USE lib_mpp        ! distribued memory computing library
    USE in_out_manager ! I/O manager
+   USE prtctl         ! Print control
 
    IMPLICIT NONE
@@ -301,5 +302,5 @@
       !!----------------------------------------------------------------------
       !! * Modules used
-      USE cpl_oce 
+      USE cpl_oce       ! coupled ocean-atmosphere variables
 
       !! * Arguments
@@ -400,13 +401,13 @@
       !!----------------------------------------------------------------------
       !! * Modules used
-      USE daymod
+      USE daymod                 ! calendar
 #if ! defined key_dtasst
-      USE dtasst, ONLY : rclice
+      USE dtasst, ONLY : rclice  ! sea surface temperature data
 #endif
 #if defined key_flx_bulk_monthly || defined key_flx_bulk_daily
-      USE blk_oce
+      USE blk_oce                ! bulk variables
 #endif
 #if defined key_flx_forced_daily
-      USE flx_oce
+      USE flx_oce                ! sea-ice/ocean forcings variables
 #endif
 
@@ -769,4 +770,5 @@
       REAL(wp), DIMENSION(jpi,jpj)  :: zsss, zfreeze
       REAL(wp) ::   zerp, ztrp, zsrp
+      CHARACTER (len=71) :: charout
 #if defined key_dynspg_fsc
       REAL(wp) ::   zwei
@@ -852,5 +854,9 @@
       IF( lk_mpp )   CALL mpp_sum( aplus  )   ! sums over the global domain
       IF( lk_mpp )   CALL mpp_sum( aminus )
-      IF(l_ctl)   WRITE(numout,*) ' oce_sbc_dmp : a+ = ', aplus, ' a- = ', aminus
+
+      IF(ln_ctl)   THEN
+         WRITE(charout,FMT="('oce_sbc_dmp : a+ = ',D23.16, ' a- = ',D23.16)") aplus, aminus
+         CALL prt_ctl_info(charout)
+      ENDIF
 
       zadefi = MIN( aplus, aminus )
Index: trunk/NEMO/OPA_SRC/TRA/traadv_cen2.F90
===================================================================
--- trunk/NEMO/OPA_SRC/TRA/traadv_cen2.F90	(revision 257)
+++ trunk/NEMO/OPA_SRC/TRA/traadv_cen2.F90	(revision 258)
@@ -24,4 +24,5 @@
    USE dynspg_fsc      ! 
    USE dynspg_fsc_atsk ! 
+   USE prtctl          ! Print control
 
    IMPLICIT NONE
@@ -299,9 +300,7 @@
       ENDIF
 
-      IF(l_ctl) THEN
-         zta = SUM( ta(2:nictl,2:njctl,1:jpkm1) * tmask(2:nictl,2:njctl,1:jpkm1) )
-         zsa = SUM( sa(2:nictl,2:njctl,1:jpkm1) * tmask(2:nictl,2:njctl,1:jpkm1) )
-         WRITE(numout,*) ' had  - Ta: ', zta-t_ctl, ' Sa: ', zsa-s_ctl
-         t_ctl = zta   ;   s_ctl = zsa
+      IF(ln_ctl)   THEN
+          CALL prt_ctl(tab3d_1=ta, clinfo1=' centered2 had  - Ta: ', mask1=tmask, &
+             &         tab3d_2=sa, clinfo2=' Sa: ', mask2=tmask, clinfo3='tra')
       ENDIF
 
@@ -400,9 +399,7 @@
       ENDIF
 
-      IF(l_ctl) THEN
-         zta = SUM( ta(2:nictl,2:njctl,1:jpkm1) * tmask(2:nictl,2:njctl,1:jpkm1) )
-         zsa = SUM( sa(2:nictl,2:njctl,1:jpkm1) * tmask(2:nictl,2:njctl,1:jpkm1) )
-         WRITE(numout,*) ' zad  - Ta: ', zta-t_ctl, ' Sa: ', zsa-s_ctl, ' centered2'
-         t_ctl = zta   ;   s_ctl = zsa
+      IF(ln_ctl)   THEN
+          CALL prt_ctl(tab3d_1=ta, clinfo1=' centered2 zad  - Ta: ', mask1=tmask, &
+             &         tab3d_2=sa, clinfo2=' Sa: ', mask2=tmask, clinfo3='tra')
       ENDIF
 
Index: trunk/NEMO/OPA_SRC/TRA/traadv_cen2_atsk.h90
===================================================================
--- trunk/NEMO/OPA_SRC/TRA/traadv_cen2_atsk.h90	(revision 257)
+++ trunk/NEMO/OPA_SRC/TRA/traadv_cen2_atsk.h90	(revision 258)
@@ -252,9 +252,7 @@
       ENDIF
 
-      IF(l_ctl) THEN
-         zta = SUM( ta(2:nictl,2:njctl,1:jpkm1) * tmask(2:nictl,2:njctl,1:jpkm1) )
-         zsa = SUM( sa(2:nictl,2:njctl,1:jpkm1) * tmask(2:nictl,2:njctl,1:jpkm1) )
-         WRITE(numout,*) ' had  - Ta: ', zta-t_ctl, ' Sa: ', zsa-s_ctl
-         t_ctl = zta   ;   s_ctl = zsa
+      IF(ln_ctl)   THEN
+          CALL prt_ctl(tab3d_1=ta, clinfo1=' centered2 had  - Ta: ', mask1=tmask, &
+             &         tab3d_2=sa, clinfo2=' Sa: ', mask2=tmask, clinfo3='tra')
       ENDIF
 
@@ -350,9 +348,7 @@
       ENDIF
 
-      IF(l_ctl) THEN
-         zta = SUM( ta(2:nictl,2:njctl,1:jpkm1) * tmask(2:nictl,2:njctl,1:jpkm1) )
-         zsa = SUM( sa(2:nictl,2:njctl,1:jpkm1) * tmask(2:nictl,2:njctl,1:jpkm1) )
-         WRITE(numout,*) ' zad  - Ta: ', zta-t_ctl, ' Sa: ', zsa-s_ctl, ' centered2 autotsk'
-         t_ctl = zta   ;   s_ctl = zsa
+      IF(ln_ctl)   THEN
+          CALL prt_ctl(tab3d_1=ta, clinfo1=' centered2 zad  - Ta: ', mask1=tmask, &
+             &         tab3d_2=sa, clinfo2=' Sa: ', mask2=tmask, clinfo3='tra')
       ENDIF
 
Index: trunk/NEMO/OPA_SRC/TRA/traadv_muscl.F90
===================================================================
--- trunk/NEMO/OPA_SRC/TRA/traadv_muscl.F90	(revision 257)
+++ trunk/NEMO/OPA_SRC/TRA/traadv_muscl.F90	(revision 258)
@@ -21,4 +21,5 @@
    USE lbclnk          ! ocean lateral boundary condition (or mpp link) 
    USE diaptr          ! poleward transport diagnostics
+   USE prtctl          ! Print control
 
    IMPLICIT NONE
@@ -279,9 +280,7 @@
       ENDIF
 
-      IF(l_ctl) THEN
-         zta = SUM( ta(2:nictl,2:njctl,1:jpkm1) * tmask(2:nictl,2:njctl,1:jpkm1) )
-         zsa = SUM( sa(2:nictl,2:njctl,1:jpkm1) * tmask(2:nictl,2:njctl,1:jpkm1) )
-         WRITE(numout,*) ' had  - Ta: ', zta-t_ctl, ' Sa: ', zsa-s_ctl, ' muscl'
-         t_ctl = zta   ;   s_ctl = zsa
+      IF(ln_ctl) THEN
+         CALL prt_ctl(tab3d_1=ta, clinfo1=' muscl had  - Ta: ', mask1=tmask ,  &
+            &         tab3d_2=sa, clinfo2=' Sa: ', mask2=tmask, clinfo3='tra' )
       ENDIF
 
@@ -412,9 +411,7 @@
       ENDIF
 
-      IF(l_ctl) THEN
-         zta = SUM( ta(2:nictl,2:njctl,1:jpkm1) * tmask(2:nictl,2:njctl,1:jpkm1) )
-         zsa = SUM( sa(2:nictl,2:njctl,1:jpkm1) * tmask(2:nictl,2:njctl,1:jpkm1) )
-         WRITE(numout,*) ' zad  - Ta: ', zta-t_ctl, ' Sa: ', zsa-s_ctl, ' muscl'
-         t_ctl = zta   ;   s_ctl = zsa
+      IF(ln_ctl) THEN
+         CALL prt_ctl(tab3d_1=ta, clinfo1=' muscl zad  - Ta: ', mask1=tmask ,  &
+            &         tab3d_2=sa, clinfo2=' Sa: ', mask2=tmask, clinfo3='tra')
       ENDIF
 
Index: trunk/NEMO/OPA_SRC/TRA/traadv_muscl2.F90
===================================================================
--- trunk/NEMO/OPA_SRC/TRA/traadv_muscl2.F90	(revision 257)
+++ trunk/NEMO/OPA_SRC/TRA/traadv_muscl2.F90	(revision 258)
@@ -21,4 +21,5 @@
    USE lbclnk          ! ocean lateral boundary condition (or mpp link) 
    USE diaptr          ! poleward transport diagnostics
+   USE prtctl          ! Print control
 
    IMPLICIT NONE
@@ -340,9 +341,7 @@
       ENDIF
 
-      IF(l_ctl) THEN
-         zta = SUM( ta(2:nictl,2:njctl,1:jpkm1) * tmask(2:nictl,2:njctl,1:jpkm1) )
-         zsa = SUM( sa(2:nictl,2:njctl,1:jpkm1) * tmask(2:nictl,2:njctl,1:jpkm1) )
-         WRITE(numout,*) ' had  - Ta: ', zta-t_ctl, ' Sa: ', zsa-s_ctl, ' muscl'
-         t_ctl = zta   ;   s_ctl = zsa
+      IF(ln_ctl) THEN
+         CALL prt_ctl(tab3d_1=ta, clinfo1=' muscl2 had  - Ta: ', mask1=tmask, &
+            &         tab3d_2=sa, clinfo2=' Sa: ', mask2=tmask, clinfo3='tra')
       ENDIF
 
@@ -486,9 +485,8 @@
       ENDIF
 
-      IF(l_ctl) THEN
-         zta = SUM( ta(2:nictl,2:njctl,1:jpkm1) * tmask(2:nictl,2:njctl,1:jpkm1) )
-         zsa = SUM( sa(2:nictl,2:njctl,1:jpkm1) * tmask(2:nictl,2:njctl,1:jpkm1) )
-         WRITE(numout,*) ' zad  - Ta: ', zta-t_ctl, ' Sa: ', zsa-s_ctl, ' muscl2'
-         t_ctl = zta   ;   s_ctl = zsa
+      IF(ln_ctl) THEN
+         CALL prt_ctl(tab3d_1=ta, clinfo1=' muscl2 zad  - Ta: ', mask1=tmask, &
+            &         tab3d_2=sa, clinfo2=' Sa: ', mask2=tmask, clinfo3='tra')
+
       ENDIF
 
Index: trunk/NEMO/OPA_SRC/TRA/traadv_tvd.F90
===================================================================
--- trunk/NEMO/OPA_SRC/TRA/traadv_tvd.F90	(revision 257)
+++ trunk/NEMO/OPA_SRC/TRA/traadv_tvd.F90	(revision 258)
@@ -23,4 +23,5 @@
    USE lbclnk          ! ocean lateral boundary condition (or mpp link) 
    USE diaptr          ! poleward transport diagnostics
+   USE prtctl          ! Print control
 
 
@@ -86,5 +87,5 @@
       !! * Local declarations
       INTEGER  ::   ji, jj, jk              ! dummy loop indices
-      REAL(wp) ::   zta, zsa,            &  ! temporary scalar
+      REAL(wp) ::                        &  ! temporary scalar
          ztai, ztaj, ztak,               &  !    "         "   
          zsai, zsaj, zsak                   !    "         "   
@@ -343,9 +344,7 @@
       ENDIF
 
-      IF(l_ctl) THEN
-         zta = SUM( ta(2:nictl,2:njctl,1:jpkm1) * tmask(2:nictl,2:njctl,1:jpkm1) )
-         zsa = SUM( sa(2:nictl,2:njctl,1:jpkm1) * tmask(2:nictl,2:njctl,1:jpkm1) )
-         WRITE(numout,*) ' zad  - Ta: ', zta-t_ctl, ' Sa: ', zsa-s_ctl, ' tvd'
-         t_ctl = zta   ;   s_ctl = zsa
+      IF(ln_ctl) THEN
+         CALL prt_ctl(tab3d_1=ta, clinfo1=' tvd adv  - Ta: ', mask1=tmask, &
+            &         tab3d_2=sa, clinfo2=' Sa: ', mask2=tmask, clinfo3='tra')
       ENDIF
 
Index: trunk/NEMO/OPA_SRC/TRA/trabbc.F90
===================================================================
--- trunk/NEMO/OPA_SRC/TRA/trabbc.F90	(revision 257)
+++ trunk/NEMO/OPA_SRC/TRA/trabbc.F90	(revision 258)
@@ -16,4 +16,5 @@
    USE phycst          ! physical constants
    USE in_out_manager  ! I/O manager
+   USE prtctl          ! Print control
 
    IMPLICIT NONE
@@ -83,5 +84,4 @@
       INTEGER ::   ji, jj              ! dummy loop indices
 #endif
-      REAL(wp) ::   zta                ! temporary scalar
       !!----------------------------------------------------------------------
 
@@ -107,8 +107,6 @@
 #endif
 
-      IF(l_ctl) THEN
-         zta = SUM( ta(2:nictl,2:njctl,1:jpkm1) * tmask(2:nictl,2:njctl,1:jpkm1) )
-         WRITE(numout,*) ' bbc  - Ta: ', zta-t_ctl
-         t_ctl = zta
+      IF(ln_ctl) THEN
+         CALL prt_ctl(tab3d_1=ta, clinfo1=' bbc  - Ta: ', mask1=tmask, clinfo3='tra-ta')
       ENDIF
 
Index: trunk/NEMO/OPA_SRC/TRA/trabbl.F90
===================================================================
--- trunk/NEMO/OPA_SRC/TRA/trabbl.F90	(revision 257)
+++ trunk/NEMO/OPA_SRC/TRA/trabbl.F90	(revision 258)
@@ -20,4 +20,5 @@
    USE trdmod_oce      ! ocean variables trends
    USE in_out_manager  ! I/O manager
+   USE prtctl          ! Print control
 
    IMPLICIT NONE
@@ -370,9 +371,7 @@
       ENDIF
 
-      IF(l_ctl) THEN
-         zta = SUM( ta(2:nictl,2:njctl,1:jpkm1) * tmask(2:nictl,2:njctl,1:jpkm1) )
-         zsa = SUM( sa(2:nictl,2:njctl,1:jpkm1) * tmask(2:nictl,2:njctl,1:jpkm1) )
-         WRITE(numout,*) ' bbl  - Ta: ', zta-t_ctl, ' Sa: ', zsa-s_ctl
-         t_ctl = zta   ;   s_ctl = zsa
+      IF(ln_ctl) THEN
+         CALL prt_ctl(tab3d_1=ta, clinfo1=' bbl  - Ta: ', mask1=tmask, &
+            &         tab3d_2=sa, clinfo2=' Sa: ', mask2=tmask, clinfo3='tra')
       ENDIF
 
Index: trunk/NEMO/OPA_SRC/TRA/trabbl_adv.h90
===================================================================
--- trunk/NEMO/OPA_SRC/TRA/trabbl_adv.h90	(revision 257)
+++ trunk/NEMO/OPA_SRC/TRA/trabbl_adv.h90	(revision 258)
@@ -451,9 +451,7 @@
       ENDIF
 
-      IF(l_ctl) THEN
-         zta = SUM( ta(2:nictl,2:njctl,1:jpkm1) * tmask(2:nictl,2:njctl,1:jpkm1) )
-         zsa = SUM( sa(2:nictl,2:njctl,1:jpkm1) * tmask(2:nictl,2:njctl,1:jpkm1) )
-         WRITE(numout,*) ' bbl  - Ta: ', zta-t_ctl, ' Sa: ', zsa-s_ctl
-         t_ctl = zta   ;   s_ctl = zsa
+      IF(ln_ctl) THEN
+         CALL prt_ctl(tab3d_1=ta, clinfo1=' bbl  - Ta: ', mask1=tmask, &
+            &         tab3d_2=sa, clinfo2=' Sa: ', mask2=tmask, clinfo3='tra')
       ENDIF
 
Index: trunk/NEMO/OPA_SRC/TRA/tradmp.F90
===================================================================
--- trunk/NEMO/OPA_SRC/TRA/tradmp.F90	(revision 257)
+++ trunk/NEMO/OPA_SRC/TRA/tradmp.F90	(revision 258)
@@ -26,4 +26,5 @@
    USE zdfmxl          ! mixed layer depth
    USE lib_mpp         ! distribued memory computing
+   USE prtctl          ! Print control
 
    IMPLICIT NONE
@@ -184,9 +185,7 @@
       ENDIF
 
-      IF(l_ctl) THEN         ! print mean trends (used for debugging)
-         zta = SUM( ta(2:nictl,2:njctl,1:jpkm1) * tmask(2:nictl,2:njctl,1:jpkm1) )
-         zsa = SUM( sa(2:nictl,2:njctl,1:jpkm1) * tmask(2:nictl,2:njctl,1:jpkm1) )
-         WRITE(numout,*) ' dmp  - Ta: ', zta-t_ctl, ' Sa: ', zsa-s_ctl
-         t_ctl = zta   ;   s_ctl = zsa
+      IF(ln_ctl) THEN         ! print mean trends (used for debugging)
+         CALL prt_ctl(tab3d_1=ta, clinfo1=' dmp  - Ta: ', mask1=tmask, &
+            &         tab3d_2=sa, clinfo2=' Sa: ', mask2=tmask, clinfo3='tra')
       ENDIF
 
Index: trunk/NEMO/OPA_SRC/TRA/traldf_bilap.F90
===================================================================
--- trunk/NEMO/OPA_SRC/TRA/traldf_bilap.F90	(revision 257)
+++ trunk/NEMO/OPA_SRC/TRA/traldf_bilap.F90	(revision 258)
@@ -19,4 +19,5 @@
    USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
    USE diaptr          ! poleward transport diagnostics
+   USE prtctl          ! Print control
 
    IMPLICIT NONE
@@ -231,9 +232,7 @@
       ENDIF
 
-      IF(l_ctl) THEN         ! print mean trends (used for debugging)
-         zta = SUM( ta(2:nictl,2:njctl,1:jpkm1) * tmask(2:nictl,2:njctl,1:jpkm1) )
-         zsa = SUM( sa(2:nictl,2:njctl,1:jpkm1) * tmask(2:nictl,2:njctl,1:jpkm1) )
-         WRITE(numout,*) ' ldf  - Ta: ', zta-t_ctl, ' Sa: ', zsa-s_ctl
-         t_ctl = zta   ;   s_ctl = zsa
+      IF(ln_ctl) THEN         ! print mean trends (used for debugging)
+         CALL prt_ctl(tab3d_1=ta, clinfo1=' ldf  - Ta: ', mask1=tmask, &
+            &         tab3d_2=sa, clinfo2=' Sa: ', mask2=tmask, clinfo3='tra')
       ENDIF
 
Index: trunk/NEMO/OPA_SRC/TRA/traldf_bilapg.F90
===================================================================
--- trunk/NEMO/OPA_SRC/TRA/traldf_bilapg.F90	(revision 257)
+++ trunk/NEMO/OPA_SRC/TRA/traldf_bilapg.F90	(revision 258)
@@ -22,4 +22,5 @@
    USE lbclnk          ! ocean lateral boundary condition (or mpp link)
    USE diaptr          ! poleward transport diagnostics 
+   USE prtctl          ! Print control
 
    IMPLICIT NONE
@@ -80,5 +81,4 @@
       !! * Local declarations
       INTEGER ::   ji, jj, jk                 ! dummy loop indices
-      REAL(wp) ::   zta, zsa                  ! temporary scalars 
       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   &
          wk3, wk4                ! work array used for rotated biharmonic
@@ -134,9 +134,7 @@
       ENDIF
 
-      IF(l_ctl) THEN         ! print mean trends (used for debugging)
-         zta = SUM( ta(2:nictl,2:njctl,1:jpkm1) * tmask(2:nictl,2:njctl,1:jpkm1) )
-         zsa = SUM( sa(2:nictl,2:njctl,1:jpkm1) * tmask(2:nictl,2:njctl,1:jpkm1) )
-         WRITE(numout,*) ' ldf  - Ta: ', zta-t_ctl, ' Sa: ', zsa-s_ctl
-         t_ctl = zta   ;   s_ctl = zsa
+      IF(ln_ctl) THEN         ! print mean trends (used for debugging)
+         CALL prt_ctl(tab3d_1=ta, clinfo1=' ldf  - Ta: ', mask1=tmask, &
+            &         tab3d_2=sa, clinfo2=' Sa: ', mask2=tmask, clinfo3='tra')
       ENDIF
 
Index: trunk/NEMO/OPA_SRC/TRA/traldf_iso.F90
===================================================================
--- trunk/NEMO/OPA_SRC/TRA/traldf_iso.F90	(revision 257)
+++ trunk/NEMO/OPA_SRC/TRA/traldf_iso.F90	(revision 258)
@@ -21,4 +21,5 @@
    USE ldfslp          ! iso-neutral slopes
    USE diaptr          ! poleward transport diagnostics
+   USE prtctl          ! Print control
 
    IMPLICIT NONE
@@ -254,9 +255,7 @@
       ENDIF
 
-      IF(l_ctl) THEN         ! print mean trends (used for debugging)
-         zta = SUM( ta(2:nictl,2:njctl,1:jpkm1) * tmask(2:nictl,2:njctl,1:jpkm1) )
-         zsa = SUM( sa(2:nictl,2:njctl,1:jpkm1) * tmask(2:nictl,2:njctl,1:jpkm1) )
-         WRITE(numout,*) ' ldf  - Ta: ', zta-t_ctl, ' Sa: ', zsa-s_ctl
-         t_ctl = zta   ;   s_ctl = zsa
+      IF(ln_ctl) THEN         ! print mean trends (used for debugging)
+         CALL prt_ctl(tab3d_1=ta, clinfo1=' ldf  - Ta: ', mask1=tmask, &
+            &         tab3d_2=sa, clinfo2=' Sa: ', mask2=tmask, clinfo3='tra')
       ENDIF
 
Index: trunk/NEMO/OPA_SRC/TRA/traldf_iso_zps.F90
===================================================================
--- trunk/NEMO/OPA_SRC/TRA/traldf_iso_zps.F90	(revision 257)
+++ trunk/NEMO/OPA_SRC/TRA/traldf_iso_zps.F90	(revision 258)
@@ -21,4 +21,5 @@
    USE ldfslp          ! iso-neutral slopes
    USE diaptr          ! poleward transport diagnostics
+   USE prtctl          ! Print control
 
 
@@ -279,9 +280,7 @@
       ENDIF
 
-      IF(l_ctl) THEN         ! print mean trends (used for debugging)
-         zta = SUM( ta(2:nictl,2:njctl,1:jpkm1) * tmask(2:nictl,2:njctl,1:jpkm1) )
-         zsa = SUM( sa(2:nictl,2:njctl,1:jpkm1) * tmask(2:nictl,2:njctl,1:jpkm1) )
-         WRITE(numout,*) ' ldf  - Ta: ', zta-t_ctl, ' Sa: ', zsa-s_ctl
-         t_ctl = zta   ;   s_ctl = zsa
+      IF(ln_ctl) THEN         ! print mean trends (used for debugging)
+         CALL prt_ctl(tab3d_1=ta, clinfo1=' ldf  - Ta: ', mask1=tmask, &
+            &         tab3d_2=sa, clinfo2=' Sa: ', mask2=tmask, clinfo3='tra')
       ENDIF
 
Index: trunk/NEMO/OPA_SRC/TRA/traldf_lap.F90
===================================================================
--- trunk/NEMO/OPA_SRC/TRA/traldf_lap.F90	(revision 257)
+++ trunk/NEMO/OPA_SRC/TRA/traldf_lap.F90	(revision 258)
@@ -17,4 +17,5 @@
    USE in_out_manager  ! I/O manager
    USE diaptr          ! poleward transport diagnostics
+   USE prtctl          ! Print control
 
 
@@ -157,9 +158,7 @@
       ENDIF
 
-      IF(l_ctl) THEN         ! print mean trends (used for debugging)
-         zta = SUM( ta(2:nictl,2:njctl,1:jpkm1) * tmask(2:nictl,2:njctl,1:jpkm1) )
-         zsa = SUM( sa(2:nictl,2:njctl,1:jpkm1) * tmask(2:nictl,2:njctl,1:jpkm1) )
-         WRITE(numout,*) ' ldf  - Ta: ', zta-t_ctl, ' Sa: ', zsa-s_ctl
-         t_ctl = zta   ;   s_ctl = zsa
+      IF(ln_ctl) THEN         ! print mean trends (used for debugging)
+         CALL prt_ctl(tab3d_1=ta, clinfo1=' ldf  - Ta: ', mask1=tmask, &
+            &         tab3d_2=sa, clinfo2=' Sa: ', mask2=tmask, clinfo3='tra')
       ENDIF
 
Index: trunk/NEMO/OPA_SRC/TRA/tranxt.F90
===================================================================
--- trunk/NEMO/OPA_SRC/TRA/tranxt.F90	(revision 257)
+++ trunk/NEMO/OPA_SRC/TRA/tranxt.F90	(revision 258)
@@ -15,4 +15,5 @@
    USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
    USE obctra          ! open boundary condition (obc_tra routine)
+   USE prtctl          ! Print control
 
    IMPLICIT NONE
@@ -165,7 +166,7 @@
       !                                                ! ===============
 
-      IF(l_ctl) THEN         ! print mean field (used for debugging)
-         WRITE(numout,*) ' nxt  - Tn: ', SUM(tn(2:nictl,2:njctl,1:jpkm1)*tmask(2:nictl,2:njctl,1:jpkm1)), &
-         &                      ' Sn: ', SUM(sn(2:nictl,2:njctl,1:jpkm1)*tmask(2:nictl,2:njctl,1:jpkm1)) 
+      IF(ln_ctl) THEN         ! print mean field (used for debugging)
+         CALL prt_ctl(tab3d_1=tn, clinfo1=' nxt  - Tn: ', mask1=tmask, &
+            &         tab3d_2=sn, clinfo2=' Sn: ', mask2=tmask)
       ENDIF
 
Index: trunk/NEMO/OPA_SRC/TRA/traqsr.F90
===================================================================
--- trunk/NEMO/OPA_SRC/TRA/traqsr.F90	(revision 257)
+++ trunk/NEMO/OPA_SRC/TRA/traqsr.F90	(revision 258)
@@ -15,9 +15,8 @@
    USE trdmod_oce      ! ocean variables trends
    USE in_out_manager  ! I/O manager
-
    USE trc_oce         ! share SMS/Ocean variables
-
    USE ocesbc          ! thermohaline fluxes
    USE phycst          ! physical constants
+   USE prtctl          ! Print control
 
    IMPLICIT NONE
@@ -235,8 +234,6 @@
 
 
-      IF(l_ctl) THEN         ! print mean trends (used for debugging)
-         zta = SUM( ta(2:nictl,2:njctl,1:jpkm1) * tmask(2:nictl,2:njctl,1:jpkm1) )
-         WRITE(numout,*) ' qsr  - Ta: ', zta-t_ctl
-         t_ctl = zta 
+      IF(ln_ctl) THEN         ! print mean trends (used for debugging)
+         CALL prt_ctl(tab3d_1=ta, clinfo1=' qsr  - Ta: ', mask1=tmask, clinfo3='tra-ta')
       ENDIF
 
Index: trunk/NEMO/OPA_SRC/TRA/trasbc.F90
===================================================================
--- trunk/NEMO/OPA_SRC/TRA/trasbc.F90	(revision 257)
+++ trunk/NEMO/OPA_SRC/TRA/trasbc.F90	(revision 258)
@@ -17,4 +17,5 @@
    USE trdmod_oce      ! ocean variables trends
    USE in_out_manager  ! I/O manager
+   USE prtctl          ! Print control
 
    IMPLICIT NONE
@@ -138,9 +139,7 @@
       ENDIF
       
-      IF(l_ctl) THEN         ! print mean trends (used for debugging)
-         zta = SUM( ta(2:nictl,2:njctl,1:jpkm1) * tmask(2:nictl,2:njctl,1:jpkm1) )
-         zsa = SUM( sa(2:nictl,2:njctl,1:jpkm1) * tmask(2:nictl,2:njctl,1:jpkm1) )
-         WRITE(numout,*) ' sbc  - Ta: ', zta, ' Sa: ', zsa
-         t_ctl = zta   ;   s_ctl = zsa
+      IF(ln_ctl) THEN         ! print mean trends (used for debugging)
+         CALL prt_ctl(tab3d_1=ta, clinfo1=' sbc  - Ta: ', mask1=tmask, &
+            &         tab3d_2=sa, clinfo2=' Sa: ', mask2=tmask, clinfo3='tra')
       ENDIF
 
Index: trunk/NEMO/OPA_SRC/TRA/trazdf_exp.F90
===================================================================
--- trunk/NEMO/OPA_SRC/TRA/trazdf_exp.F90	(revision 257)
+++ trunk/NEMO/OPA_SRC/TRA/trazdf_exp.F90	(revision 258)
@@ -18,4 +18,5 @@
    USE zdfddm          ! ocean vertical physics: double diffusion
    USE in_out_manager  ! I/O manager
+   USE prtctl          ! Print control
 
    IMPLICIT NONE
@@ -180,9 +181,7 @@
       ENDIF
 
-      IF(l_ctl) THEN         ! print mean trends (used for debugging)
-         zta = SUM( ta(2:nictl,2:njctl,1:jpkm1) * tmask(2:nictl,2:njctl,1:jpkm1) )
-         zsa = SUM( sa(2:nictl,2:njctl,1:jpkm1) * tmask(2:nictl,2:njctl,1:jpkm1) )
-         WRITE(numout,*) ' zdf  - Ta: ', zta-t_ctl, ' Sa: ', zsa-s_ctl
-         t_ctl = zta   ;   s_ctl = zsa
+      IF(ln_ctl) THEN         ! print mean trends (used for debugging)
+         CALL prt_ctl(tab3d_1=ta, clinfo1=' zdf  - Ta: ', mask1=tmask, &
+            &         tab3d_2=sa, clinfo2=' Sa: ', mask2=tmask, clinfo3='tra')
       ENDIF
 
Index: trunk/NEMO/OPA_SRC/TRA/trazdf_imp.F90
===================================================================
--- trunk/NEMO/OPA_SRC/TRA/trazdf_imp.F90	(revision 257)
+++ trunk/NEMO/OPA_SRC/TRA/trazdf_imp.F90	(revision 258)
@@ -19,4 +19,5 @@
    USE trdmod_oce      ! ocean variables trends
    USE in_out_manager  ! I/O manager
+   USE prtctl          ! Print control
 
 
@@ -84,5 +85,4 @@
       INTEGER ::   ji, jj, jk                 ! dummy loop indices
       INTEGER ::   ikst, ikenm2, ikstp1
-      REAL(wp)  ::   zta, zsa                 ! temporary scalars
       REAL(wp), DIMENSION(jpi,jpk) ::   &
          zwd, zws, zwi,          &  ! ???
@@ -243,9 +243,7 @@
       ENDIF
 
-      IF(l_ctl) THEN         ! print mean trends (used for debugging)
-         zta = SUM( ta(2:nictl,2:njctl,1:jpkm1) * tmask(2:nictl,2:njctl,1:jpkm1) )
-         zsa = SUM( sa(2:nictl,2:njctl,1:jpkm1) * tmask(2:nictl,2:njctl,1:jpkm1) )
-         WRITE(numout,*) ' zdf  - Ta: ', zta-t_ctl, ' Sa: ', zsa-s_ctl
-         t_ctl = zta   ;   s_ctl = zsa
+      IF(ln_ctl) THEN         ! print mean trends (used for debugging)
+         CALL prt_ctl(tab3d_1=ta, clinfo1=' zdf  - Ta: ', mask1=tmask, &
+            &         tab3d_2=sa, clinfo2=' Sa: ', mask2=tmask, clinfo3='tra')
       ENDIF
 
Index: trunk/NEMO/OPA_SRC/TRA/trazdf_iso.F90
===================================================================
--- trunk/NEMO/OPA_SRC/TRA/trazdf_iso.F90	(revision 257)
+++ trunk/NEMO/OPA_SRC/TRA/trazdf_iso.F90	(revision 258)
@@ -24,4 +24,5 @@
    USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
    USE zdfkpp          ! KPP parameterisation
+   USE prtctl          ! Print control
 
    IMPLICIT NONE
@@ -408,9 +409,7 @@
       ENDIF
 
-      IF(l_ctl) THEN         ! print mean trends (used for debugging)
-         zta = SUM( ta(2:nictl,2:njctl,1:jpkm1) * tmask(2:nictl,2:njctl,1:jpkm1) )
-         zsa = SUM( sa(2:nictl,2:njctl,1:jpkm1) * tmask(2:nictl,2:njctl,1:jpkm1) )
-         WRITE(numout,*) ' zdf 1- Ta: ', zta-t_ctl, ' Sa: ', zsa-s_ctl
-         t_ctl = zta   ;   s_ctl = zsa
+      IF(ln_ctl) THEN         ! print mean trends (used for debugging)
+         CALL prt_ctl(tab3d_1=ta, clinfo1=' zdf 1- Ta: ', mask1=tmask, &
+            &         tab3d_2=sa, clinfo2=' Sa: ', mask2=tmask, clinfo3='tra')
       ENDIF
 
@@ -578,8 +577,7 @@
       ENDIF
 
-      IF(l_ctl) THEN         ! print mean trends (used for debugging)
-         zta = SUM( ta(2:nictl,2:njctl,1:jpkm1) * tmask(2:nictl,2:njctl,1:jpkm1) )
-         zsa = SUM( sa(2:nictl,2:njctl,1:jpkm1) * tmask(2:nictl,2:njctl,1:jpkm1) )
-         WRITE(numout,*) ' zdf 2- Ta: ', zta, ' Sa: ', zsa
+      IF(ln_ctl) THEN         ! print mean trends (used for debugging)
+         CALL prt_ctl(tab3d_1=ta, clinfo1=' zdf 2- Ta: ', mask1=tmask, &
+            &         tab3d_2=sa, clinfo2=' Sa: ', mask2=tmask)
       ENDIF
 
Index: trunk/NEMO/OPA_SRC/TRA/trazdf_iso_vopt.F90
===================================================================
--- trunk/NEMO/OPA_SRC/TRA/trazdf_iso_vopt.F90	(revision 257)
+++ trunk/NEMO/OPA_SRC/TRA/trazdf_iso_vopt.F90	(revision 258)
@@ -27,4 +27,5 @@
    USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
    USE zdfkpp          ! KPP parameterisation
+   USE prtctl          ! Print control
 
    IMPLICIT NONE
@@ -70,6 +71,4 @@
       !! * Arguments
       INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
-      !! * Local variables
-      REAL(wp) :: zta, zsa               ! temporary scalars
 
       IF( kt == nit000 ) THEN
@@ -91,9 +90,7 @@
       CALL tra_zdf_iso( kt )
 
-      IF(l_ctl) THEN         ! print mean trends (used for debugging)
-         zta = SUM( ta(2:nictl,2:njctl,1:jpkm1) * tmask(2:nictl,2:njctl,1:jpkm1) )
-         zsa = SUM( sa(2:nictl,2:njctl,1:jpkm1) * tmask(2:nictl,2:njctl,1:jpkm1) )
-         WRITE(numout,*) ' zdf 1- Ta: ', zta-t_ctl, ' Sa: ', zsa-s_ctl
-         t_ctl = zta   ;   s_ctl = zsa
+      IF(ln_ctl) THEN         ! print mean trends (used for debugging)
+         CALL prt_ctl(tab3d_1=ta, clinfo1=' zdf 1- Ta: ', mask1=tmask, &
+            &         tab3d_2=sa, clinfo2=' Sa: ', mask2=tmask, clinfo3='tra')
       ENDIF
 
@@ -103,8 +100,7 @@
       CALL tra_zdf_zdf( kt )
 
-      IF(l_ctl) THEN         ! print mean trends (used for debugging)
-         zta = SUM( ta(2:nictl,2:njctl,1:jpkm1) * tmask(2:nictl,2:njctl,1:jpkm1) )
-         zsa = SUM( sa(2:nictl,2:njctl,1:jpkm1) * tmask(2:nictl,2:njctl,1:jpkm1) )
-         WRITE(numout,*) ' zdf 2- Ta: ', zta, ' Sa: ', zsa
+      IF(ln_ctl) THEN         ! print mean trends (used for debugging)
+         CALL prt_ctl(tab3d_1=ta, clinfo1=' zdf 2- Ta: ', mask1=tmask, &
+            &         tab3d_2=sa, clinfo2=' Sa: ', mask2=tmask)
       ENDIF
 
Index: trunk/NEMO/OPA_SRC/ZDF/zdfbfr.F90
===================================================================
--- trunk/NEMO/OPA_SRC/ZDF/zdfbfr.F90	(revision 257)
+++ trunk/NEMO/OPA_SRC/ZDF/zdfbfr.F90	(revision 258)
@@ -17,4 +17,5 @@
    USE in_out_manager  ! I/O manager
    USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
+   USE prtctl          ! Print control
 
    IMPLICIT NONE
@@ -170,5 +171,7 @@
       CALL lbc_lnk( avmv, 'V', 1. )
 
-      IF(l_ctl)   WRITE(numout,*) ' bfr  u : ', SUM( avmu(1:nictl+1,1:njctl+1,:) ), ' v : ', SUM( avmv(1:nictl+1,1:njctl+1,:) )
+      IF(ln_ctl)   THEN
+         CALL prt_ctl(tab3d_1=avmu, clinfo1=' bfr  - u: ', tab3d_2=avmv, clinfo2=' v: ', ovlap=1, kdim=jpk)
+      ENDIF
 
    END SUBROUTINE zdf_bfr
Index: trunk/NEMO/OPA_SRC/ZDF/zdfddm.F90
===================================================================
--- trunk/NEMO/OPA_SRC/ZDF/zdfddm.F90	(revision 257)
+++ trunk/NEMO/OPA_SRC/ZDF/zdfddm.F90	(revision 258)
@@ -17,4 +17,5 @@
    USE in_out_manager  ! I/O manager
    USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
+   USE prtctl          ! Print control
 
    IMPLICIT NONE
@@ -201,7 +202,7 @@
       CALL lbc_lnk( avmv, 'V', 1. )
 
-      IF(l_ctl) THEN
-         WRITE(numout,*) ' ddm  t : ', SUM( avt (1:nictl+1,1:njctl+1,:) ), ' s : ', SUM( avs (1:nictl+1,1:njctl+1,:) )
-         WRITE(numout,*) '      u : ', SUM( avmu(1:nictl+1,1:njctl+1,:) ), ' v : ', SUM( avmv(1:nictl+1,1:njctl+1,:) )
+      IF(ln_ctl) THEN
+         CALL prt_ctl(tab3d_1=avt , clinfo1=' ddm  - t: ', tab3d_2=avs , clinfo2=' s: ', ovlap=1, kdim=jpk)
+         CALL prt_ctl(tab3d_1=avmu, clinfo1=' ddm  - u: ', tab3d_2=avmv, clinfo2=' v: ', ovlap=1, kdim=jpk)
       ENDIF
       
Index: trunk/NEMO/OPA_SRC/ZDF/zdfkpp.F90
===================================================================
--- trunk/NEMO/OPA_SRC/ZDF/zdfkpp.F90	(revision 257)
+++ trunk/NEMO/OPA_SRC/ZDF/zdfkpp.F90	(revision 258)
@@ -23,4 +23,5 @@
    USE ocesbc          ! thermohaline fluxes
    USE zdfddm          ! double diffusion mixing
+   USE prtctl          ! Print control
 
    IMPLICIT NONE
@@ -145,6 +146,6 @@
    !!----------------------------------------------------------------------
    !!   OPA 9.0 , LOCEAN-IPSL   (2005)
-   !! $Header$
-   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
+   !! $Header$ 
+   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 
    !!----------------------------------------------------------------------
 
@@ -1247,4 +1248,15 @@
          CALL lbc_lnk( avmu(:,:,:), 'U', 1. )   ;    CALL lbc_lnk( avmv(:,:,:), 'V', 1. )  
  
+         IF(ln_ctl) THEN
+#if defined key_zdfddm
+            CALL prt_ctl(tab3d_1=avt  , clinfo1=' kpp  - t: ', tab3d_2=avs , clinfo2=' s: ', &
+               &         ovlap=1, kdim=jpk)
+#else
+            CALL prt_ctl(tab3d_1=avt  , clinfo1=' kpp  - t: ', ovlap=1, kdim=jpk)
+#endif
+            CALL prt_ctl(tab3d_1=avmu  , clinfo1='       u: ', tab3d_2=avmv , clinfo2=' v: ', &
+               &         ovlap=1, kdim=jpk)
+         ENDIF
+
    END SUBROUTINE zdf_kpp
 
@@ -1268,5 +1280,4 @@
       !!     9.0  ! 05-01 (C. Ethe) F90 : free form
       !!----------------------------------------------------------------------
-      !! * Module used
       !! * local declarations
 
Index: trunk/NEMO/OPA_SRC/ZDF/zdfmxl.F90
===================================================================
--- trunk/NEMO/OPA_SRC/ZDF/zdfmxl.F90	(revision 257)
+++ trunk/NEMO/OPA_SRC/ZDF/zdfmxl.F90	(revision 258)
@@ -13,4 +13,5 @@
    USE zdf_oce         ! ocean vertical physics
    USE in_out_manager  ! I/O manager
+   USE prtctl          ! Print control
 
    IMPLICIT NONE
@@ -134,5 +135,8 @@
       !                                                ! ===============
 
-      IF(l_ctl)   WRITE(numout,*) ' nmln : ', SUM( nmln(1:nictl+1,1:njctl+1) ), ' hmld: ', SUM( hmld(1:nictl+1,1:njctl+1) )
+      IF(ln_ctl)   THEN
+         CALL prt_ctl(tab2d_1=FLOAT(nmln), clinfo1=' nmln : ', tab2d_2=hmld, clinfo2=' hmld : ', ovlap=1)
+      ENDIF
+      
 
    END SUBROUTINE zdf_mxl
@@ -235,5 +239,7 @@
       END DO
 
-      IF(l_ctl)   WRITE(numout,*) ' nmln : ', SUM( nmln(1:nictl+1,1:njctl+1) ), ' hmld: ', SUM( hmld(1:nictl+1,1:njctl+1) )
+      IF(ln_ctl)   THEN
+         CALL prt_ctl(tab2d_1=FLOAT(nmln), clinfo1=' nmln : ', tab2d_2=hmld, clinfo2=' hmld : ', ovlap=1)
+      ENDIF
 
    END SUBROUTINE zdf_mxl
Index: trunk/NEMO/OPA_SRC/ZDF/zdftke.F90
===================================================================
--- trunk/NEMO/OPA_SRC/ZDF/zdftke.F90	(revision 257)
+++ trunk/NEMO/OPA_SRC/ZDF/zdftke.F90	(revision 258)
@@ -20,4 +20,5 @@
    USE phycst          ! physical constants
    USE taumod          ! surface stress
+   USE prtctl          ! Print control
 
    IMPLICIT NONE
@@ -614,7 +615,7 @@
       CALL lbc_lnk( avt, 'W', 1. )
 
-      IF(l_ctl) THEN
-         WRITE(numout,*) ' tke  e : ', SUM( en  (1:nictl+1,1:njctl+1,:) ), ' t : ', SUM( avt (1:nictl+1,1:njctl+1,:) )
-         WRITE(numout,*) '      u : ', SUM( avmu(1:nictl+1,1:njctl+1,:) ), ' v : ', SUM( avmv(1:nictl+1,1:njctl+1,:) )
+      IF(ln_ctl) THEN
+         CALL prt_ctl(tab3d_1=en  , clinfo1=' tke  - e: ', tab3d_2=avt , clinfo2=' t: ', ovlap=1, kdim=jpk)
+         CALL prt_ctl(tab3d_1=avmu, clinfo1=' tke  - u: ', tab3d_2=avmv, clinfo2=' v: ', ovlap=1, kdim=jpk)
       ENDIF
 
Index: trunk/NEMO/OPA_SRC/ZDF/zdftke_atsk.h90
===================================================================
--- trunk/NEMO/OPA_SRC/ZDF/zdftke_atsk.h90	(revision 257)
+++ trunk/NEMO/OPA_SRC/ZDF/zdftke_atsk.h90	(revision 258)
@@ -478,7 +478,7 @@
       CALL lbc_lnk( avt, 'W', 1. )
 
-      IF(l_ctl) THEN
-         WRITE(numout,*) ' tke  e : ', SUM( en  (1:nictl+1,1:njctl+1,:) ), ' t : ', SUM( avt (1:nictl+1,1:njctl+1,:) )
-         WRITE(numout,*) '      u : ', SUM( avmu(1:nictl+1,1:njctl+1,:) ), ' v : ', SUM( avmv(1:nictl+1,1:njctl+1,:) )
+      IF(ln_ctl) THEN
+         CALL prt_ctl(tab3d_1=en  , clinfo1=' tke  - e: ', tab3d_2=avt , clinfo2=' t: ', ovlap=1, kdim=jpk)
+         CALL prt_ctl(tab3d_1=avmu, clinfo1=' tke  - u: ', tab3d_2=avmv, clinfo2=' v: ', ovlap=1, kdim=jpk)
       ENDIF
 
Index: trunk/NEMO/OPA_SRC/daymod.F90
===================================================================
--- trunk/NEMO/OPA_SRC/daymod.F90	(revision 257)
+++ trunk/NEMO/OPA_SRC/daymod.F90	(revision 258)
@@ -12,4 +12,5 @@
    USE phycst          ! physical constants
    USE in_out_manager  ! I/O manager
+   USE prtctl          ! Print control
 
    IMPLICIT NONE
@@ -69,4 +70,5 @@
       INTEGER  ::   iend, iday0, iday1   ! temporary integers
       REAL(wp) :: zadatrjn, zadatrjb     ! adatrj at timestep kt-1 and kt-2 
+      CHARACTER (len=25) :: charout
       !!----------------------------------------------------------------------
 
@@ -153,5 +155,8 @@
       ENDIF
 
-      IF(l_ctl) WRITE(numout,*)' kt =', kt, 'd/m/y =', nday, nmonth, nyear
+      IF(ln_ctl) THEN
+         WRITE(charout,FMT="('kt =', I4,'  d/m/y =',I2,I2,I4)") kt, nday, nmonth, nyear
+         CALL prt_ctl_info(charout)
+      ENDIF
 
    END SUBROUTINE day
Index: trunk/NEMO/OPA_SRC/eosbn2.F90
===================================================================
--- trunk/NEMO/OPA_SRC/eosbn2.F90	(revision 257)
+++ trunk/NEMO/OPA_SRC/eosbn2.F90	(revision 258)
@@ -20,4 +20,5 @@
    USE in_out_manager  ! I/O manager
    USE zdfddm          ! vertical physics: double diffusion
+   USE prtctl          ! Print control
 
    IMPLICIT NONE
@@ -224,5 +225,7 @@
       END SELECT
 
-      IF(l_ctl)   WRITE(numout,*) ' eos  : ', SUM( prd(1:nictl+1,1:njctl+1,:) )
+      IF(ln_ctl)   THEN
+         CALL prt_ctl(tab3d_1=prd, clinfo1=' eos  : ', ovlap=1, kdim=jpk)
+      ENDIF
 
    END SUBROUTINE eos_insitu
@@ -410,5 +413,7 @@
       END SELECT
 
-      IF(l_ctl)   WRITE(numout,*) ' eos-p: ', SUM( prd(1:nictl+1,1:njctl+1,:) ), ' pot : ', SUM( prhop(1:nictl+1,1:njctl+1,:) )
+      IF(ln_ctl)   THEN
+         CALL prt_ctl(tab3d_1=prd, clinfo1=' eos-p: ', tab3d_2=prhop, clinfo2=' pot : ', ovlap=1, kdim=jpk)
+      ENDIF
 
    END SUBROUTINE eos_insitu_pot
@@ -587,5 +592,5 @@
       END SELECT
 
-      IF(l_ctl)   WRITE(numout,*) ' eos2d: ', SUM( prd(1:nictl,1:njctl) )
+      IF(ln_ctl)   CALL prt_ctl(tab2d_1=prd, clinfo1=' eos2d: ')
 
    END SUBROUTINE eos_insitu_2d
@@ -762,8 +767,10 @@
       END SELECT
 
-      IF(l_ctl)   WRITE(numout,*) ' bn2  : ', SUM( pn2 (1:nictl+1,1:njctl+1,:) )
+      IF(ln_ctl)   THEN
+         CALL prt_ctl(tab3d_1=pn2, clinfo1=' bn2  : ', ovlap=1, kdim=jpk)
 #if defined key_zdfddm
-      IF(l_ctl)   WRITE(numout,*) ' rrau : ', SUM( rrau(1:nictl+1,1:njctl+1,:) )
+         CALL prt_ctl(tab3d_1=rrau, clinfo1=' rrau : ', ovlap=1, kdim=jpk)
 #endif
+      ENDIF
 
    END SUBROUTINE eos_bn2
Index: trunk/NEMO/OPA_SRC/in_out_manager.F90
===================================================================
--- trunk/NEMO/OPA_SRC/in_out_manager.F90	(revision 257)
+++ trunk/NEMO/OPA_SRC/in_out_manager.F90	(revision 258)
@@ -52,5 +52,5 @@
    LOGICAL ::   &               !:
       lwp                ,   &  !: boolean : true on the 1st processor only
-      l_ctl                     !: = ln_ctl.AND.lwp (print control on the 1st proc)
+      lsp_area = .TRUE.         !: to make a control print over a specific area
 
    INTEGER ::                &  !:
@@ -58,6 +58,12 @@
       nprint =    0 ,        &  !: level of print (0 no print)
       nwrite =   10 ,        &  !: restart file frequency
-      nictl  = jpim1,        &  !: max i indice for the SUM control
-      njctl  = jpjm1            !: max j indice for the SUM control
+      nictls =    0 ,        &  !: Start i indice for the SUM control
+      nictle =    0 ,        &  !: End   i indice for the SUM control
+      njctls =    0 ,        &  !: Start j indice for the SUM control
+      njctle =    0 ,        &  !: End   j indice for the SUM control
+      isplt  =    1 ,        &  !: number of processors following i
+      jsplt  =    1 ,        &  !: number of processors following j
+      ijsplt =    1             !: nb of local domain = nb of processors
+
    !!----------------------------------------------------------------------
    !! logical units
Index: trunk/NEMO/OPA_SRC/mppini.F90
===================================================================
--- trunk/NEMO/OPA_SRC/mppini.F90	(revision 257)
+++ trunk/NEMO/OPA_SRC/mppini.F90	(revision 258)
@@ -136,6 +136,4 @@
          inum                               ! temporary logical unit
 
-      INTEGER, DIMENSION(jpnij) ::       &
-         ibonit, ibonjt                     ! temporary workspace
       INTEGER, DIMENSION(jpni,jpnj) ::   &
          iimppt, ijmppt, ilcit, ilcjt       ! temporary workspace
Index: trunk/NEMO/OPA_SRC/opa.F90
===================================================================
--- trunk/NEMO/OPA_SRC/opa.F90	(revision 257)
+++ trunk/NEMO/OPA_SRC/opa.F90	(revision 258)
@@ -46,4 +46,5 @@
 
    USE step            ! OPA time-stepping                  (stp     routine)
+   USE prtctl          ! Print control                 (prt_ctl_init routine)
    USE ini1d           ! re-initialization of u-v mask for the 1D configuration
    USE dyncor1d        ! Coriolis factor at T-point
@@ -126,17 +127,4 @@
       lwp   = narea == 1
 
-      IF( lk_mpp )   THEN
-!! #############################  MPI debug  ###########################  !!
-!!  To allow each processor to write its own mpp.output_XXX file, decomment 
-!!  the 5 following lines, else the output will be done only by the first
-!!  processor in the ocean.output file:
-!!  this one        numout = 80
-!!  this one        numout = numout + narea
-!!  this one        lwp    = .TRUE.
-!!  this one        WRITE(file_out,FMT="('mpp.output_',I3.3)") narea-1
-!!  this one        OPEN( UNIT=numout, FILE=TRIM(file_out),FORM='FORMATTED' )
-!! #############################  MPI debug  ###########################  !!
-      ENDIF
-
       IF(lwp) THEN
          WRITE(numout,*)
@@ -168,4 +156,6 @@
       
       CALL dom_init                         ! Domain
+
+      IF( ln_ctl )      CALL prt_ctl_init   ! Print control
 
       IF( lk_cfg_1d )   CALL fcorio_1d      ! redefine Coriolis at T-point
Index: trunk/NEMO/OPA_SRC/prtctl.f90
===================================================================
--- trunk/NEMO/OPA_SRC/prtctl.f90	(revision 258)
+++ trunk/NEMO/OPA_SRC/prtctl.f90	(revision 258)
@@ -0,0 +1,563 @@
+MODULE prtctl
+   !!==============================================================================
+   !!                       ***  MODULE prtctl   ***
+   !! Ocean system   : print all SUM trends for each processor domain
+   !!==============================================================================
+   USE dom_oce          ! ocean space and time domain variables
+   USE in_out_manager   ! I/O manager
+   USE lib_mpp          ! distributed memory computing
+
+   IMPLICIT NONE
+   PRIVATE
+
+   !! * Module declaration
+   INTEGER, DIMENSION(:), ALLOCATABLE, SAVE ::   &  !:
+      nlditl , nldjtl ,   &  !: first, last indoor index for each i-domain
+      nleitl , nlejtl ,   &  !: first, last indoor index for each j-domain
+      nimpptl, njmpptl,   &  !: i-, j-indexes for each processor
+      nlcitl , nlcjtl ,   &  !: dimensions of every subdomain
+      ibonitl, ibonjtl
+
+   REAL(wp), DIMENSION(:), ALLOCATABLE ::   &  !:
+      t_ctll , s_ctll ,   &  !: previous trend values
+      u_ctll , v_ctll
+
+   !! * Routine accessibility
+   PUBLIC prt_ctl         ! called by all subroutines
+   PUBLIC prt_ctl_info    ! called by all subroutines
+   PUBLIC prt_ctl_init    ! called by opa.F90
+   !!----------------------------------------------------------------------
+   !!   OPA 9.0 , LOCEAN-IPSL (2005) 
+   !! $Header$ 
+   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 
+   !!----------------------------------------------------------------------
+
+
+CONTAINS
+
+   SUBROUTINE prt_ctl (tab2d_1, tab3d_1, mask1, clinfo1, tab2d_2, tab3d_2, mask2, clinfo2, ovlap, kdim, clinfo3)
+      !!----------------------------------------------------------------------
+      !!                     ***  ROUTINE prt_ctl  ***
+      !!
+      !! ** Purpose : - print sum control of 2D or 3D arrays over the same area 
+      !!                in mono and mpp case. This way can be usefull when
+      !!                debugging a new parametrization in mono or mpp. 
+      !!
+      !! ** Method  : 2 possibilities exist when setting the ln_ctl parameter to
+      !!                .true. in the ocean namelist:
+      !!              - to debug a MPI run .vs. a mono-processor one; 
+      !!                the control print will be done over each sub-domain.
+      !!                The nictl[se] and njctl[se] parameters in the namelist must 
+      !!                be set to zero and [ij]splt to the corresponding splitted
+      !!                domain in MPI along respectively i-, j- directions.
+      !!              - to debug a mono-processor run over the whole domain/a specific area; 
+      !!                in the first case the nictl[se] and njctl[se] parameters must be set
+      !!                to zero else to the indices of the area to be controled. In both cases
+      !!                isplt and jsplt must be set to 1.
+      !!              - All arguments of the above calling sequence are optional so their
+      !!                name must be explicitly typed if used. For instance if the 3D
+      !!                array tn(:,:,:) must be passed through the prt_ctl subroutine, 
+      !!                it must looks like: CALL prt_ctl(tab3d_1=tn).
+      !!
+      !!                    tab2d_1 : first 2D array
+      !!                    tab3d_1 : first 3D array
+      !!                    mask1   : mask (3D) to apply to the tab[23]d_1 array
+      !!                    clinfo1 : information about the tab[23]d_1 array
+      !!                    tab2d_2 : second 2D array
+      !!                    tab3d_2 : second 3D array
+      !!                    mask2   : mask (3D) to apply to the tab[23]d_2 array
+      !!                    clinfo2 : information about the tab[23]d_2 array
+      !!                    ovlap   : overlap value
+      !!                    kdim    : k- direction for 3D arrays 
+      !!                    clinfo3 : additional information 
+      !!
+      !! History :
+      !!   9.0  !  05-07  (C. Talandier) original code
+      !!----------------------------------------------------------------------
+      !! * Arguments
+      REAL(wp), DIMENSION(:,:)  , INTENT(in), OPTIONAL :: tab2d_1
+      REAL(wp), DIMENSION(:,:,:), INTENT(in), OPTIONAL :: tab3d_1
+      REAL(wp), DIMENSION(:,:,:), INTENT(in), OPTIONAL :: mask1
+      CHARACTER (len=*), INTENT(in), OPTIONAL :: clinfo1
+      REAL(wp), DIMENSION(:,:)  , INTENT(in), OPTIONAL :: tab2d_2
+      REAL(wp), DIMENSION(:,:,:), INTENT(in), OPTIONAL :: tab3d_2
+      REAL(wp), DIMENSION(:,:,:), INTENT(in), OPTIONAL :: mask2
+      CHARACTER (len=*), INTENT(in), OPTIONAL :: clinfo2
+      INTEGER, INTENT(in), OPTIONAL :: ovlap
+      INTEGER, INTENT(in), OPTIONAL :: kdim
+      CHARACTER (len=*), INTENT(in), OPTIONAL :: clinfo3
+
+      !! * Local declarations
+      INTEGER :: overlap, numid, jn, sind, eind, kdir
+      CHARACTER (len=15) :: cl2
+      REAL(wp) :: zsum1, zsum2, zvctl1, zvctl2
+      REAL(wp), DIMENSION(jpi,jpj) :: ztab2d_1, ztab2d_2
+      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zmask1, zmask2, ztab3d_1, ztab3d_2
+      !!----------------------------------------------------------------------
+
+      ! Arrays, scalars initialization 
+      overlap   = 0
+      kdir      = jpkm1
+      cl2       = ''
+      zsum1     = 0.e0
+      zsum2     = 0.e0
+      zvctl1    = 0.e0
+      zvctl2    = 0.e0
+      ztab2d_1(:,:)   = 0.e0
+      ztab2d_2(:,:)   = 0.e0
+      ztab3d_1(:,:,:) = 0.e0
+      ztab3d_2(:,:,:) = 0.e0
+      zmask1  (:,:,:) = 1.e0
+      zmask2  (:,:,:) = 1.e0
+
+      ! Control of optional arguments
+      IF( PRESENT(clinfo2) )  cl2            = clinfo2
+      IF( PRESENT(ovlap)   )  overlap        = ovlap
+      IF( PRESENT(kdim)    )  kdir           = kdim
+      IF( PRESENT(tab2d_1) )  ztab2d_1(:,:)  = tab2d_1(:,:)
+      IF( PRESENT(tab2d_2) )  ztab2d_2(:,:)  = tab2d_2(:,:)
+      IF( PRESENT(tab3d_1) )  ztab3d_1(:,:,:)= tab3d_1(:,:,:)
+      IF( PRESENT(tab3d_2) )  ztab3d_2(:,:,:)= tab3d_2(:,:,:)
+      IF( PRESENT(mask1)   )  zmask1  (:,:,:)= mask1  (:,:,:)
+      IF( PRESENT(mask2)   )  zmask2  (:,:,:)= mask2  (:,:,:)
+
+      IF( lk_mpp )   THEN
+         ! processor number
+         sind = narea
+         eind = narea
+      ELSE
+         ! processors total number
+         sind = 1
+         eind = ijsplt
+      ENDIF
+
+      ! Loop over each sub-domain, i.e. the total number of processors ijsplt
+      DO jn = sind, eind
+
+         numid = 80 + jn
+
+         ! Set indices for the SUM control
+         IF( .NOT. lsp_area ) THEN
+            IF (lk_mpp )   THEN
+               nictls = MAX( 1, nlditl(jn) - overlap )
+               nictle = nleitl(jn) + overlap * MIN( 1, nlcitl(jn) - nleitl(jn)) 
+               njctls = MAX( 1, nldjtl(jn) - overlap )
+               njctle = nlejtl(jn) + overlap * MIN( 1, nlcjtl(jn) - nlejtl(jn))
+               ! Do not take into account the bound of the domain
+               IF( ibonitl(jn) == -1 ) nictls = MAX(2, nictls)
+               IF( ibonitl(jn) ==  1 ) nictle = MIN(nictle, nleitl(jn) - 1)
+               IF( ibonjtl(jn) == -1 ) njctls = MAX(2, njctls)
+               IF( ibonjtl(jn) ==  1 ) njctle = MIN(njctle, nlejtl(jn) - 1)
+            ELSE
+               nictls = MAX( 1, nimpptl(jn) + nlditl(jn) - 1 - overlap )
+               nictle = nimpptl(jn) + nleitl(jn) - 1 + overlap * MIN( 1, nlcitl(jn) - nleitl(jn) ) 
+               njctls = MAX( 1, njmpptl(jn) + nldjtl(jn) - 1 - overlap )
+               njctle = njmpptl(jn) + nlejtl(jn) - 1 + overlap * MIN( 1, nlcjtl(jn) - nlejtl(jn) ) 
+               ! Do not take into account the bound of the domain
+               IF( ibonitl(jn) == -1 ) nictls = MAX(2, nictls)
+               IF( ibonjtl(jn) == -1 ) njctls = MAX(2, njctls)
+               IF( ibonitl(jn) ==  1 ) nictle = MIN(nictle, nimpptl(jn) + nleitl(jn) - 2)
+               IF( ibonjtl(jn) ==  1 ) njctle = MIN(njctle, njmpptl(jn) + nlejtl(jn) - 2)
+            ENDIF
+         ENDIF
+
+         IF ( clinfo3 == 'tra' )  THEN
+             zvctl1 = t_ctll(jn)
+             zvctl2 = s_ctll(jn)
+         ELSEIF ( clinfo3 == 'dyn' )   THEN
+             zvctl1 = u_ctll(jn)
+             zvctl2 = v_ctll(jn)
+         ENDIF
+
+         ! Compute the sum control
+         ! 2D arrays
+         IF( PRESENT(tab2d_1) )   THEN
+            zsum1 = SUM( ztab2d_1(nictls:nictle,njctls:njctle)*zmask1(nictls:nictle,njctls:njctle,1) )
+            zsum2 = SUM( ztab2d_2(nictls:nictle,njctls:njctle)*zmask2(nictls:nictle,njctls:njctle,1) )
+         ENDIF
+
+         ! 3D arrays
+         IF( PRESENT(tab3d_1) )   THEN
+            zsum1 = SUM( ztab3d_1(nictls:nictle,njctls:njctle,1:kdir)*zmask1(nictls:nictle,njctls:njctle,1:kdir) )
+            zsum2 = SUM( ztab3d_2(nictls:nictle,njctls:njctle,1:kdir)*zmask2(nictls:nictle,njctls:njctle,1:kdir) )
+         ENDIF
+
+         ! Print the result
+         IF( PRESENT(clinfo2) .OR. PRESENT(tab2d_2) .OR. PRESENT(tab3d_2) )   THEN
+            WRITE(numid,FMT='(a,D23.16,3x,a,D23.16)')clinfo1, zsum1-zvctl1, cl2, zsum2-zvctl2
+             IF( clinfo3 == 'tra-ta' )  THEN
+                t_ctll(jn) = zsum1
+             ENDIF 
+             IF( clinfo3 == 'tra' )  THEN
+                t_ctll(jn) = zsum1
+                s_ctll(jn) = zsum2
+             ENDIF 
+             IF( clinfo3 == 'dyn' )   THEN
+                u_ctll(jn) = zsum1
+                v_ctll(jn) = zsum2 
+             ENDIF
+         ELSE
+            WRITE(numid,FMT='(a,D23.16)')clinfo1, zsum1
+         ENDIF
+
+      ENDDO
+
+   END SUBROUTINE prt_ctl
+
+
+   SUBROUTINE prt_ctl_info (clinfo1, ivar1, clinfo2, ivar2)
+      !!----------------------------------------------------------------------
+      !!                     ***  ROUTINE prt_ctl_info  ***
+      !!
+      !! ** Purpose : - print information without any computation
+      !!
+      !! ** Action  : - input arguments
+      !!                    clinfo1 : information about the ivar1
+      !!                    ivar1   : value to print
+      !!                    clinfo2 : information about the ivar2
+      !!                    ivar2   : value to print
+      !!
+      !! History :
+      !!   9.0  !  05-07  (C. Talandier) original code
+      !!----------------------------------------------------------------------
+      !! * Arguments
+      CHARACTER (len=*), INTENT(in) ::   clinfo1
+      INTEGER, INTENT(in), OPTIONAL ::   ivar1
+      CHARACTER (len=*), INTENT(in), OPTIONAL ::   clinfo2
+      INTEGER, INTENT(in), OPTIONAL ::   ivar2
+
+      !! * Local declarations
+      INTEGER ::  numid, jn, sind, eind
+      !!----------------------------------------------------------------------
+
+      IF( lk_mpp )   THEN
+         ! processor number
+         sind = narea
+         eind = narea
+      ELSE
+         ! total number of processors
+         sind = 1
+         eind = ijsplt
+      ENDIF
+
+      ! Loop over each sub-domain, i.e. number of processors ijsplt
+      DO jn = sind, eind
+
+         numid = 80 + jn
+
+         IF( PRESENT(ivar1) .AND. PRESENT(clinfo2) .AND. PRESENT(ivar2) )   THEN
+            WRITE(numid,*)clinfo1, ivar1, clinfo2, ivar2
+         ELSEIF ( PRESENT(ivar1) .AND. PRESENT(clinfo2) .AND. .NOT. PRESENT(ivar2) )   THEN
+            WRITE(numid,*)clinfo1, ivar1, clinfo2
+         ELSEIF ( PRESENT(ivar1) .AND. .NOT. PRESENT(clinfo2) .AND. PRESENT(ivar2) )   THEN
+            WRITE(numid,*)clinfo1, ivar1, ivar2
+         ELSEIF ( PRESENT(ivar1) .AND. .NOT. PRESENT(clinfo2) .AND. .NOT. PRESENT(ivar2) )   THEN
+            WRITE(numid,*)clinfo1, ivar1
+         ELSE
+            WRITE(numid,*)clinfo1
+         ENDIF
+
+      ENDDO
+
+
+      END SUBROUTINE prt_ctl_info
+
+   SUBROUTINE prt_ctl_init
+      !!----------------------------------------------------------------------
+      !!                     ***  ROUTINE prt_ctl_init  ***
+      !!
+      !! ** Purpose :   open ASCII files & compute indices
+      !!
+      !! History :
+      !!   9.0  !  05-07  (C. Talandier) original code
+      !!----------------------------------------------------------------------
+      !! * Local declarations
+      INTEGER ::   jn, numid, sind, eind
+      CHARACTER (len=28) :: clfile_out
+      CHARACTER (len=23) :: clb_name
+      CHARACTER (len=19) :: cl_run
+      !!----------------------------------------------------------------------
+
+      ! Allocate arrays
+      ALLOCATE(nlditl (ijsplt))
+      ALLOCATE(nldjtl (ijsplt))
+      ALLOCATE(nleitl (ijsplt))
+      ALLOCATE(nlejtl (ijsplt))
+      ALLOCATE(nimpptl(ijsplt))
+      ALLOCATE(njmpptl(ijsplt))
+      ALLOCATE(nlcitl (ijsplt))
+      ALLOCATE(nlcjtl (ijsplt))
+      ALLOCATE(t_ctll (ijsplt))
+      ALLOCATE(s_ctll (ijsplt))
+      ALLOCATE(u_ctll (ijsplt))
+      ALLOCATE(v_ctll (ijsplt))
+      ALLOCATE(ibonitl(ijsplt))
+      ALLOCATE(ibonjtl(ijsplt))
+
+      ! Initialization 
+      t_ctll (:)=0.e0
+      s_ctll (:)=0.e0
+      u_ctll (:)=0.e0
+      v_ctll (:)=0.e0
+
+      IF( lk_mpp ) THEN
+         sind = narea
+         eind = narea
+         clb_name = "('mpp.output_',I3.3)"
+         cl_run = 'MULTI processor run'
+         ! use indices for each area computed by mpp_init subroutine
+         nlditl(:) = nldit(:) 
+         nleitl(:) = nleit(:) 
+         nldjtl(:) = nldjt(:) 
+         nlejtl(:) = nlejt(:) 
+         !
+         nimpptl(:) = nimppt(:)
+         njmpptl(:) = njmppt(:)
+         !
+         nlcitl(:) = nlcit(:)
+         nlcjtl(:) = nlcjt(:)
+         !
+         ibonitl(:) = ibonit(:)
+         ibonjtl(:) = ibonjt(:)
+      ELSE
+         sind = 1
+         eind = ijsplt
+         clb_name = "('mono.output_',I3.3)"
+         cl_run = 'MONO processor run '
+         ! compute indices for each area as done in mpp_init subroutine
+         CALL sub_dom
+      ENDIF
+
+      DO jn = sind, eind
+         numid = 80 + jn
+         WRITE(clfile_out,FMT=clb_name) jn-1
+         OPEN ( UNIT=numid, FILE=TRIM(clfile_out),FORM='FORMATTED' )
+         WRITE(numid,*)
+         WRITE(numid,*) '                 L O D Y C - I P S L'
+         WRITE(numid,*) '                     O P A model'
+         WRITE(numid,*) '            Ocean General Circulation Model'
+         WRITE(numid,*) '               version OPA 9.0  (2005) '
+         WRITE(numid,*)
+         WRITE(numid,*) '                   PROC number: ', jn
+         WRITE(numid,*)
+         WRITE(numid,FMT="(19x,a20)")cl_run
+
+         ! Print the SUM control indices
+         IF( .NOT. lsp_area )   THEN
+            IF ( lk_mpp )   THEN
+               nictls = nlditl(jn) 
+               nictle = nleitl(jn)
+               njctls = nldjtl(jn)
+               njctle = nlejtl(jn)
+            ELSE
+               nictls = nimpptl(jn) + nlditl(jn) - 1
+               nictle = nimpptl(jn) + nleitl(jn) - 1
+               njctls = njmpptl(jn) + nldjtl(jn) - 1
+               njctle = njmpptl(jn) + nlejtl(jn) - 1
+            ENDIF
+         ENDIF
+         WRITE(numid,*) 
+         WRITE(numid,*) 'prt_ctl :  Sum control indices'
+         WRITE(numid,*) '~~~~~~~'
+         WRITE(numid,*)
+         WRITE(numid,9000)'                                nlej   = ', nlejtl(jn), '              '
+         WRITE(numid,9000)'                  ------------- njctle = ', njctle, ' -------------'
+         WRITE(numid,9001)'                  |                                       |'
+         WRITE(numid,9001)'                  |                                       |'
+         WRITE(numid,9001)'                  |                                       |'
+         WRITE(numid,9002)'           nictls = ', nictls,  '                           nictle = ', nictle
+         WRITE(numid,9002)'           nldi   = ', nlditl(jn),  '                           nlei   = ', nleitl(jn)
+         WRITE(numid,9001)'                  |                                       |'
+         WRITE(numid,9001)'                  |                                       |'
+         WRITE(numid,9001)'                  |                                       |'
+         WRITE(numid,9004)'  njmpp  = ',njmpptl(jn),'   ------------- njctls = ', njctls, ' -------------'
+         WRITE(numid,9003)'           nimpp  = ', nimpptl(jn), '        nldj   = ', nldjtl(jn), '              '
+         WRITE(numid,*)
+         WRITE(numid,*)
+
+9000     FORMAT(a41,i4.4,a14)
+9001     FORMAT(a59)
+9002     FORMAT(a20,i4.4,a36,i3.3)
+9003     FORMAT(a20,i4.4,a17,i4.4)
+9004     FORMAT(a11,i4.4,a26,i4.4,a14)
+      ENDDO
+
+   END SUBROUTINE prt_ctl_init
+
+
+   SUBROUTINE sub_dom
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE sub_dom  ***
+      !!                    
+      !! ** Purpose :   Lay out the global domain over processors. 
+      !!                CAUTION: 
+      !!                This part has been extracted from the mpp_init
+      !!                subroutine and names of variables/arrays have been 
+      !!                slightly changed to avoid confusion but the computation
+      !!                is exactly the same. Any modification about indices of
+      !!                each sub-domain in the mppini.F90 module should be reported 
+      !!                here.
+      !!
+      !! ** Method  :   Global domain is distributed in smaller local domains.
+      !!                Periodic condition is a function of the local domain position
+      !!                (global boundary or neighbouring domain) and of the global
+      !!                periodic
+      !!                Type :         jperio global periodic condition
+      !!                               nperio local  periodic condition
+      !!
+      !! ** Action  : - set domain parameters
+      !!                    nimpp     : longitudinal index 
+      !!                    njmpp     : latitudinal  index
+      !!                    nperio    : lateral condition type 
+      !!                    narea     : number for local area
+      !!                    nlcil      : first dimension
+      !!                    nlcjl      : second dimension
+      !!                    nbondil    : mark for "east-west local boundary"
+      !!                    nbondjl    : mark for "north-south local boundary"
+      !!
+      !! History :
+      !!        !  94-11  (M. Guyon)  Original code
+      !!        !  95-04  (J. Escobar, M. Imbard)
+      !!        !  98-02  (M. Guyon)  FETI method
+      !!        !  98-05  (M. Imbard, J. Escobar, L. Colombet )  SHMEM and MPI versions
+      !!   8.5  !  02-08  (G. Madec)  F90 : free form
+      !!----------------------------------------------------------------------
+      !! * Local variables
+      INTEGER ::   ji, jj, jn               ! dummy loop indices
+      INTEGER ::   &
+         ii, ij, ifreq, il1, il2,        &  ! temporary integers
+         irestil, irestjl,               &  !    "          "
+         ijpi  , ijpj, nlcil,            &  ! temporary logical unit
+         nlcjl , nbondil, nbondjl,       &
+         nrecil, nrecjl, nldil, nleil, nldjl, nlejl
+
+      INTEGER, DIMENSION(:,:), ALLOCATABLE ::   &
+         iimpptl, ijmpptl, ilcitl, ilcjtl       ! temporary workspace
+      REAL(wp) ::   zidom, zjdom            ! temporary scalars
+      !!----------------------------------------------------------------------
+
+      !  1. Dimension arrays for subdomains
+      ! -----------------------------------
+      !  Computation of local domain sizes ilcitl() ilcjtl()
+      !  These dimensions depend on global sizes isplt,jsplt and jpiglo,jpjglo
+      !  The subdomains are squares leeser than or equal to the global
+      !  dimensions divided by the number of processors minus the overlap
+      !  array (cf. par_oce.F90).
+
+      ijpi = ( jpiglo-2*jpreci + (isplt-1) ) / isplt + 2*jpreci
+      ijpj = ( jpjglo-2*jprecj + (jsplt-1) ) / jsplt + 2*jprecj
+
+      ALLOCATE(ilcitl (isplt,jsplt))
+      ALLOCATE(ilcjtl (isplt,jsplt))
+
+      nrecil  = 2 * jpreci
+      nrecjl  = 2 * jprecj
+      irestil = MOD( jpiglo - nrecil , isplt )
+      irestjl = MOD( jpjglo - nrecjl , jsplt )
+
+      IF(  irestil == 0 )   irestil = isplt
+      DO jj = 1, jsplt
+         DO ji = 1, irestil
+            ilcitl(ji,jj) = ijpi
+         END DO
+         DO ji = irestil+1, isplt
+            ilcitl(ji,jj) = ijpi -1
+         END DO
+      END DO
+      
+      IF( irestjl == 0 )   irestjl = jsplt
+      DO ji = 1, isplt
+         DO jj = 1, irestjl
+            ilcjtl(ji,jj) = ijpj
+         END DO
+         DO jj = irestjl+1, jsplt
+            ilcjtl(ji,jj) = ijpj -1
+         END DO
+      END DO
+      
+      zidom = nrecil
+      DO ji = 1, isplt
+         zidom = zidom + ilcitl(ji,1) - nrecil
+      END DO
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*)' sum ilcitl(i,1) = ', zidom, ' jpiglo = ', jpiglo
+      
+      zjdom = nrecjl
+      DO jj = 1, jsplt
+         zjdom = zjdom + ilcjtl(1,jj) - nrecjl
+      END DO
+      IF(lwp) WRITE(numout,*)' sum ilcitl(1,j) = ', zjdom, ' jpjglo = ', jpjglo
+      IF(lwp) WRITE(numout,*)
+      
+
+      !  2. Index arrays for subdomains
+      ! -------------------------------
+
+      ALLOCATE(iimpptl(isplt,jsplt))
+      ALLOCATE(ijmpptl(isplt,jsplt))
+      
+      iimpptl(:,:) = 1
+      ijmpptl(:,:) = 1
+      
+      IF( isplt > 1 ) THEN
+         DO jj = 1, jsplt
+            DO ji = 2, isplt
+               iimpptl(ji,jj) = iimpptl(ji-1,jj) + ilcitl(ji-1,jj) - nrecil
+            END DO
+         END DO
+      ENDIF
+
+      IF( jsplt > 1 ) THEN
+         DO jj = 2, jsplt
+            DO ji = 1, isplt
+               ijmpptl(ji,jj) = ijmpptl(ji,jj-1)+ilcjtl(ji,jj-1)-nrecjl
+            END DO
+         END DO
+      ENDIF
+      
+      ! 3. Subdomain description
+      ! ------------------------
+
+      DO jn = 1, ijsplt
+         ii = 1 + MOD( jn-1, isplt )
+         ij = 1 + (jn-1) / isplt
+         nimpptl(jn) = iimpptl(ii,ij)
+         njmpptl(jn) = ijmpptl(ii,ij)
+         nlcitl (jn) = ilcitl (ii,ij)     
+         nlcil       = nlcitl (jn)     
+         nlcjtl (jn) = ilcjtl (ii,ij)     
+         nlcjl       = nlcjtl (jn)
+         nbondjl = -1                                    ! general case
+         IF( jn   >  isplt          )   nbondjl = 0      ! first row of processor
+         IF( jn   >  (jsplt-1)*isplt )  nbondjl = 1     ! last  row of processor
+         IF( jsplt == 1             )   nbondjl = 2      ! one processor only in j-direction
+         ibonjtl(jn) = nbondjl
+         
+         nbondil = 0                                     ! 
+         IF( MOD( jn, isplt ) == 1 )   nbondil = -1      !
+         IF( MOD( jn, isplt ) == 0 )   nbondil =  1      !
+         IF( isplt            == 1 )   nbondil =  2      ! one processor only in i-direction
+         ibonitl(jn) = nbondil
+         
+         nldil =  1   + jpreci
+         nleil = nlcil - jpreci
+         IF( nbondil == -1 .OR. nbondil == 2 )   nldil = 1
+         IF( nbondil ==  1 .OR. nbondil == 2 )   nleil = nlcil
+         nldjl =  1   + jprecj
+         nlejl = nlcjl - jprecj
+         IF( nbondjl == -1 .OR. nbondjl == 2 )   nldjl = 1
+         IF( nbondjl ==  1 .OR. nbondjl == 2 )   nlejl = nlcjl
+         nlditl(jn) = nldil
+         nleitl(jn) = nleil
+         nldjtl(jn) = nldjl
+         nlejtl(jn) = nlejl
+      END DO
+
+      DEALLOCATE(iimpptl)
+      DEALLOCATE(ijmpptl)
+      DEALLOCATE(ilcitl)
+      DEALLOCATE(ilcjtl)
+
+   END SUBROUTINE sub_dom
+
+END MODULE prtctl
Index: trunk/NEMO/OPA_SRC/step.F90
===================================================================
--- trunk/NEMO/OPA_SRC/step.F90	(revision 257)
+++ trunk/NEMO/OPA_SRC/step.F90	(revision 258)
@@ -117,4 +117,5 @@
    USE restart         ! ocean restart                    (rst_wri routine)
    USE cpl             ! exchanges in coupled mode        (cpl_stp routine)
+   USE prtctl          ! Print control                    (prt_ctl routine)
 
    IMPLICIT NONE
@@ -212,14 +213,14 @@
       ENDIF
 
-      IF(l_ctl) THEN         ! print mean trends (used for debugging)
-         WRITE(numout,*) ' emp  -   : ', SUM( emp   (1:nictl+1,1:njctl+1)   * tmask(1:nictl+1,1:njctl+1,1) )
-         WRITE(numout,*) ' emps -   : ', SUM( emps  (1:nictl+1,1:njctl+1)   * tmask(1:nictl+1,1:njctl+1,1) )
-         WRITE(numout,*) ' qt   -   : ', SUM( qt    (1:nictl+1,1:njctl+1)   * tmask(1:nictl+1,1:njctl+1,1) )
-         WRITE(numout,*) ' qsr  -   : ', SUM( qsr   (1:nictl+1,1:njctl+1)   * tmask(1:nictl+1,1:njctl+1,1) )
-         WRITE(numout,*) ' runoff   : ', SUM( runoff(1:nictl+1,1:njctl+1)   * tmask(1:nictl+1,1:njctl+1,1) )
-         WRITE(numout,*) ' tmask    : ', SUM( tmask (1:nictl+1,1:njctl+1,:) )
-         WRITE(numout,*) ' sst  -   : ', SUM( tn    (1:nictl+1,1:njctl+1,1) * tmask(1:nictl+1,1:njctl+1,1) )
-         WRITE(numout,*) ' sss  -   : ', SUM( sn    (1:nictl+1,1:njctl+1,1) * tmask(1:nictl+1,1:njctl+1,1) )
-         WRITE(numout,*) ' tau  - x : ', SUM( taux  (1:nictl+1,1:njctl+1) ), ' - y : ', SUM( tauy(1:nictl+1,1:njctl+1) )
+      IF(ln_ctl) THEN         ! print mean trends (used for debugging)
+         CALL prt_ctl(tab2d_1=emp    , clinfo1=' emp  -   : ', mask1=tmask, ovlap=1)
+         CALL prt_ctl(tab2d_1=emps   , clinfo1=' emps -   : ', mask1=tmask, ovlap=1)
+         CALL prt_ctl(tab2d_1=qt     , clinfo1=' qt   -   : ', mask1=tmask, ovlap=1)
+         CALL prt_ctl(tab2d_1=qsr    , clinfo1=' qsr  -   : ', mask1=tmask, ovlap=1)
+         CALL prt_ctl(tab2d_1=runoff , clinfo1=' runoff   : ', mask1=tmask, ovlap=1)
+         CALL prt_ctl(tab3d_1=tmask  , clinfo1=' tmask    : ', mask1=tmask, ovlap=1, kdim=jpk)
+         CALL prt_ctl(tab3d_1=tn     , clinfo1=' sst  -   : ', mask1=tmask, ovlap=1, kdim=1)
+         CALL prt_ctl(tab3d_1=sn     , clinfo1=' sss  -   : ', mask1=tmask, ovlap=1, kdim=1)
+         CALL prt_ctl(tab2d_1=taux   , clinfo1=' tau  - x : ', tab2d_2=tauy, clinfo2='      - y : ', ovlap=1)
       ENDIF
 
