Index: /NEMO/trunk/src/ICE/icecor.F90
===================================================================
--- /NEMO/trunk/src/ICE/icecor.F90	(revision 13496)
+++ /NEMO/trunk/src/ICE/icecor.F90	(revision 13497)
@@ -105,5 +105,5 @@
       !                             !-----------------------------------------------------
       IF( kn == 2 ) THEN            !  Ice drift case: Corrections to avoid wrong values !
-         DO_2D( 0, 0, 0, 0 )
+         DO_2D( 0, 0, 0, 0 )        !-----------------------------------------------------
             IF ( at_i(ji,jj) == 0._wp ) THEN    ! what to do if there is no ice
                IF ( at_i(ji+1,jj) == 0._wp )   u_ice(ji  ,jj) = 0._wp   ! right side
Index: /NEMO/trunk/src/ICE/icedyn_adv_pra.F90
===================================================================
--- /NEMO/trunk/src/ICE/icedyn_adv_pra.F90	(revision 13496)
+++ /NEMO/trunk/src/ICE/icedyn_adv_pra.F90	(revision 13497)
@@ -383,5 +383,5 @@
 
          !  Calculate fluxes and moments between boxes i<-->i+1              
-         DO_2D( 0, 0, 1, 1 )
+         DO_2D( 0, 0, 1, 1 )                   !  Flux from i to i+1 WHEN u GT 0
             zbet(ji,jj)  =  MAX( 0._wp, SIGN( 1._wp, put(ji,jj) ) )
             zalf         =  MAX( 0._wp, put(ji,jj) ) * pdt / psm(ji,jj,jl)
@@ -408,5 +408,5 @@
          END_2D
 
-         DO_2D( 0, 0, 1, 0 )
+         DO_2D( 0, 0, 1, 0 )                   !  Flux from i+1 to i when u LT 0.
             zalf          = MAX( 0._wp, -put(ji,jj) ) * pdt / psm(ji+1,jj,jl) 
             zalg  (ji,jj) = zalf
@@ -427,5 +427,5 @@
          END_2D
 
-         DO_2D( 0, 0, 0, 0 )
+         DO_2D( 0, 0, 0, 0 )                   !  Readjust moments remaining in the box.
             zbt  =       zbet(ji-1,jj)
             zbt1 = 1.0 - zbet(ji-1,jj)
@@ -441,5 +441,5 @@
 
          !   Put the temporary moments into appropriate neighboring boxes.    
-         DO_2D( 0, 0, 0, 0 )
+         DO_2D( 0, 0, 0, 0 )                   !   Flux from i to i+1 IF u GT 0.
             zbt  =       zbet(ji-1,jj)
             zbt1 = 1.0 - zbet(ji-1,jj)
@@ -461,5 +461,5 @@
          END_2D
 
-         DO_2D( 0, 0, 0, 0 )
+         DO_2D( 0, 0, 0, 0 )                   !  Flux from i+1 to i IF u LT 0.
             zbt  =       zbet(ji,jj)
             zbt1 = 1.0 - zbet(ji,jj)
@@ -540,5 +540,5 @@
  
          !  Calculate fluxes and moments between boxes j<-->j+1              
-         DO_2D( 1, 1, 0, 0 )
+         DO_2D( 1, 1, 0, 0 )                   !  Flux from j to j+1 WHEN v GT 0
             zbet(ji,jj)  =  MAX( 0._wp, SIGN( 1._wp, pvt(ji,jj) ) )
             zalf         =  MAX( 0._wp, pvt(ji,jj) ) * pdt / psm(ji,jj,jl)
@@ -565,5 +565,5 @@
          END_2D
          !
-         DO_2D( 1, 0, 0, 0 )
+         DO_2D( 1, 0, 0, 0 )                   !  Flux from j+1 to j when v LT 0.
             zalf          = MAX( 0._wp, -pvt(ji,jj) ) * pdt / psm(ji,jj+1,jl) 
             zalg  (ji,jj) = zalf
@@ -599,5 +599,5 @@
 
          !   Put the temporary moments into appropriate neighboring boxes.    
-         DO_2D( 0, 0, 0, 0 )
+         DO_2D( 0, 0, 0, 0 )                   !  Flux from j to j+1 IF v GT 0.
             zbt  =       zbet(ji,jj-1)
             zbt1 = 1.0 - zbet(ji,jj-1)
@@ -620,5 +620,5 @@
          END_2D
 
-         DO_2D( 0, 0, 0, 0 )
+         DO_2D( 0, 0, 0, 0 )                   !  Flux from j+1 to j IF v LT 0.
             zbt  =       zbet(ji,jj)
             zbt1 = 1.0 - zbet(ji,jj)
Index: /NEMO/trunk/src/ICE/icedyn_adv_umx.F90
===================================================================
--- /NEMO/trunk/src/ICE/icedyn_adv_umx.F90	(revision 13496)
+++ /NEMO/trunk/src/ICE/icedyn_adv_umx.F90	(revision 13497)
@@ -1006,8 +1006,8 @@
       !                                                     !--  Laplacian in j-direction  --!
       DO jl = 1, jpl
-         DO_2D( 1, 0, 0, 0 )
+         DO_2D( 1, 0, 0, 0 )         ! First derivative (gradient)
             ztv1(ji,jj,jl) = ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) * r1_e2v(ji,jj) * vmask(ji,jj,1)
          END_2D
-         DO_2D( 0, 0, 0, 0 )
+         DO_2D( 0, 0, 0, 0 )         ! Second derivative (Laplacian)
             ztv2(ji,jj,jl) = ( ztv1(ji,jj,jl) - ztv1(ji,jj-1,jl) ) * r1_e2t(ji,jj)
          END_2D
@@ -1017,8 +1017,8 @@
       !                                                     !--  BiLaplacian in j-direction  --!
       DO jl = 1, jpl
-         DO_2D( 1, 0, 0, 0 )
+         DO_2D( 1, 0, 0, 0 )         ! First derivative
             ztv3(ji,jj,jl) = ( ztv2(ji,jj+1,jl) - ztv2(ji,jj,jl) ) * r1_e2v(ji,jj) * vmask(ji,jj,1)
          END_2D
-         DO_2D( 0, 0, 0, 0 )
+         DO_2D( 0, 0, 0, 0 )         ! Second derivative
             ztv4(ji,jj,jl) = ( ztv3(ji,jj,jl) - ztv3(ji,jj-1,jl) ) * r1_e2t(ji,jj)
          END_2D
Index: /NEMO/trunk/src/ICE/icedyn_rhg_evp.F90
===================================================================
--- /NEMO/trunk/src/ICE/icedyn_rhg_evp.F90	(revision 13496)
+++ /NEMO/trunk/src/ICE/icedyn_rhg_evp.F90	(revision 13497)
@@ -370,5 +370,5 @@
          CALL lbc_lnk( 'icedyn_rhg_evp', zds, 'F', 1.0_wp )
 
-         DO_2D( 0, 1, 0, 1 )
+         DO_2D( 0, 1, 0, 1 )   ! loop to jpi,jpj to avoid making a communication for zs1,zs2,zs12 ! no vector loop
 
             ! shear**2 at T points (doc eq. A16)
@@ -704,5 +704,5 @@
       END_2D
       
-      DO_2D( 0, 0, 0, 0 )
+      DO_2D( 0, 0, 0, 0 )   ! no vector loop
          
          ! tension**2 at T points
Index: /NEMO/trunk/src/ICE/iceupdate.F90
===================================================================
--- /NEMO/trunk/src/ICE/iceupdate.F90	(revision 13496)
+++ /NEMO/trunk/src/ICE/iceupdate.F90	(revision 13497)
@@ -332,5 +332,5 @@
       !
       IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN     !==  Ice time-step only  ==!   (i.e. surface module time-step)
-         DO_2D( 0, 0, 0, 0 )
+         DO_2D( 0, 0, 0, 0 )                          !* update the modulus of stress at ocean surface (T-point)
             !                                               ! 2*(U_ice-U_oce) at T-point
             zu_t = u_ice(ji,jj) + u_ice(ji-1,jj) - u_oce(ji,jj) - u_oce(ji-1,jj)   
@@ -358,5 +358,5 @@
       ENDIF
       !
-      DO_2D( 0, 0, 0, 0 )
+      DO_2D( 0, 0, 0, 0 )                             !* update the stress WITHOUT an ice-ocean rotation angle
          ! ice area at u and v-points 
          zat_u  = ( at_i(ji,jj) * tmask(ji,jj,1) + at_i (ji+1,jj    ) * tmask(ji+1,jj  ,1) )  &
Index: /NEMO/trunk/src/OCE/C1D/dtauvd.F90
===================================================================
--- /NEMO/trunk/src/OCE/C1D/dtauvd.F90	(revision 13496)
+++ /NEMO/trunk/src/OCE/C1D/dtauvd.F90	(revision 13497)
@@ -158,5 +158,5 @@
          ENDIF
          !
-         DO_2D( 1, 1, 1, 1 )
+         DO_2D( 1, 1, 1, 1 )           ! vertical interpolation of U & V current:
             DO jk = 1, jpk
                zl = gdept(ji,jj,jk,Kmm)
Index: /NEMO/trunk/src/OCE/DIA/diaar5.F90
===================================================================
--- /NEMO/trunk/src/OCE/DIA/diaar5.F90	(revision 13496)
+++ /NEMO/trunk/src/OCE/DIA/diaar5.F90	(revision 13497)
@@ -144,10 +144,8 @@
          IF( ln_linssh ) THEN
             IF( ln_isfcav ) THEN
-               DO ji = 1, jpi
-                  DO jj = 1, jpj
-                     iks = mikt(ji,jj)
-                     zbotpres(ji,jj) = zbotpres(ji,jj) + ssh(ji,jj,Kmm) * zrhd(ji,jj,iks) + riceload(ji,jj)
-                  END DO
-               END DO
+               DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+                  iks = mikt(ji,jj)
+                  zbotpres(ji,jj) = zbotpres(ji,jj) + ssh(ji,jj,Kmm) * zrhd(ji,jj,iks) + riceload(ji,jj)
+               END_2D
             ELSE
                zbotpres(:,:) = zbotpres(:,:) + ssh(:,:,Kmm) * zrhd(:,:,1)
@@ -385,5 +383,5 @@
          zvol0 (:,:) = 0._wp
          thick0(:,:) = 0._wp
-         DO_3D( 1, 1, 1, 1, 1, jpkm1 )
+         DO_3D( 1, 1, 1, 1, 1, jpkm1 )   ! interpolation of salinity at the last ocean level (i.e. the partial step)
             idep = tmask(ji,jj,jk) * e3t_0(ji,jj,jk)
             zvol0 (ji,jj) = zvol0 (ji,jj) +  idep * e1e2t(ji,jj)
@@ -403,5 +401,5 @@
             sn0(:,:,:) = sn0(:,:,:) * tmask(:,:,:)
             IF( ln_zps ) THEN               ! z-coord. partial steps
-               DO_2D( 1, 1, 1, 1 )
+               DO_2D( 1, 1, 1, 1 )          ! interpolation of salinity at the last ocean level (i.e. the partial step)
                   ik = mbkt(ji,jj)
                   IF( ik > 1 ) THEN
Index: /NEMO/trunk/src/OCE/DIA/diacfl.F90
===================================================================
--- /NEMO/trunk/src/OCE/DIA/diacfl.F90	(revision 13496)
+++ /NEMO/trunk/src/OCE/DIA/diacfl.F90	(revision 13497)
@@ -66,5 +66,5 @@
       llmsk(:,Nje1: jpj,:) = .FALSE.
       !
-      DO_3D( 0, 0, 0, 0, 1, jpk )
+      DO_3D( 0, 0, 0, 0, 1, jpk )      ! calculate Courant numbers
          zCu_cfl(ji,jj,jk) = ABS( uu(ji,jj,jk,Kmm) ) * rDt / e1u  (ji,jj)      ! for i-direction
          zCv_cfl(ji,jj,jk) = ABS( vv(ji,jj,jk,Kmm) ) * rDt / e2v  (ji,jj)      ! for j-direction
Index: /NEMO/trunk/src/OCE/DIA/diahth.F90
===================================================================
--- /NEMO/trunk/src/OCE/DIA/diahth.F90	(revision 13496)
+++ /NEMO/trunk/src/OCE/DIA/diahth.F90	(revision 13497)
@@ -170,5 +170,5 @@
             ! MLD: rho = rho(1) + zrho1                                     !
             ! ------------------------------------------------------------- !
-            DO_3DS( 1, 1, 1, 1, jpkm1, 2, -1 )
+            DO_3DS( 1, 1, 1, 1, jpkm1, 2, -1 )   ! loop from bottom to 2
                !
                zzdep = gdepw(ji,jj,jk,Kmm)
@@ -207,5 +207,5 @@
             ! depth of temperature inversion                                !
             ! ------------------------------------------------------------- !
-            DO_3DS( 1, 1, 1, 1, jpkm1, nlb10, -1 )
+            DO_3DS( 1, 1, 1, 1, jpkm1, nlb10, -1 )   ! loop from bottom to nlb10
                !
                zzdep = gdepw(ji,jj,jk,Kmm) * tmask(ji,jj,1)
@@ -305,5 +305,5 @@
       ! --------------------------------------- !
       iktem(:,:) = 1
-      DO_3D( 1, 1, 1, 1, 1, jpkm1 )
+      DO_3D( 1, 1, 1, 1, 1, jpkm1 )   ! beware temperature is not always decreasing with depth => loop from top to bottom
          zztmp = ts(ji,jj,jk,jp_tem,Kmm)
          IF( zztmp >= ptem )   iktem(ji,jj) = jk
Index: /NEMO/trunk/src/OCE/DIA/diawri.F90
===================================================================
--- /NEMO/trunk/src/OCE/DIA/diawri.F90	(revision 13496)
+++ /NEMO/trunk/src/OCE/DIA/diawri.F90	(revision 13497)
@@ -270,5 +270,5 @@
          
       IF ( iom_use("sstgrad") .OR. iom_use("sstgrad2") ) THEN
-         DO_2D( 0, 0, 0, 0 )
+         DO_2D( 0, 0, 0, 0 )                                 ! sst gradient
             zztmp  = ts(ji,jj,1,jp_tem,Kmm)
             zztmpx = ( ts(ji+1,jj,1,jp_tem,Kmm) - zztmp ) * r1_e1u(ji,jj) + ( zztmp - ts(ji-1,jj  ,1,jp_tem,Kmm) ) * r1_e1u(ji-1,jj)
Index: /NEMO/trunk/src/OCE/DOM/domvvl.F90
===================================================================
--- /NEMO/trunk/src/OCE/DOM/domvvl.F90	(revision 13496)
+++ /NEMO/trunk/src/OCE/DOM/domvvl.F90	(revision 13497)
@@ -202,5 +202,5 @@
       gdept(:,:,1,Kbb) = 0.5_wp * e3w(:,:,1,Kbb)
       gdepw(:,:,1,Kbb) = 0.0_wp
-      DO_3D( 1, 1, 1, 1, 2, jpk )
+      DO_3D( 1, 1, 1, 1, 2, jpk )                     ! vertical sum
          !    zcoef = tmask - wmask    ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt
          !                             ! 1 everywhere from mbkt to mikt + 1 or 1 (if no isf)
@@ -420,5 +420,5 @@
          zwu(:,:) = 0._wp
          zwv(:,:) = 0._wp
-         DO_3D( 1, 0, 1, 0, 1, jpkm1 )
+         DO_3D( 1, 0, 1, 0, 1, jpkm1 )   ! a - first derivative: diffusive fluxes
             un_td(ji,jj,jk) = rn_ahe3 * umask(ji,jj,jk) * e2_e1u(ji,jj)           &
                &            * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji+1,jj  ,jk) )
@@ -428,14 +428,14 @@
             zwv(ji,jj) = zwv(ji,jj) + vn_td(ji,jj,jk)
          END_3D
-         DO_2D( 1, 1, 1, 1 )
+         DO_2D( 1, 1, 1, 1 )             ! b - correction for last oceanic u-v points
             un_td(ji,jj,mbku(ji,jj)) = un_td(ji,jj,mbku(ji,jj)) - zwu(ji,jj)
             vn_td(ji,jj,mbkv(ji,jj)) = vn_td(ji,jj,mbkv(ji,jj)) - zwv(ji,jj)
          END_2D
-         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
+         DO_3D( 0, 0, 0, 0, 1, jpkm1 )   ! c - second derivative: divergence of diffusive fluxes
             tilde_e3t_a(ji,jj,jk) = tilde_e3t_a(ji,jj,jk) + (   un_td(ji-1,jj  ,jk) - un_td(ji,jj,jk)    &
                &                                          +     vn_td(ji  ,jj-1,jk) - vn_td(ji,jj,jk)    &
                &                                            ) * r1_e1e2t(ji,jj)
          END_3D
-         !                       ! d - thickness diffusion transport: boundary conditions
+         !                               ! d - thickness diffusion transport: boundary conditions
          !                             (stored for tracer advction and continuity equation)
          CALL lbc_lnk_multi( 'domvvl', un_td , 'U' , -1._wp, vn_td , 'V' , -1._wp)
Index: /NEMO/trunk/src/OCE/DOM/dtatsd.F90
===================================================================
--- /NEMO/trunk/src/OCE/DOM/dtatsd.F90	(revision 13496)
+++ /NEMO/trunk/src/OCE/DOM/dtatsd.F90	(revision 13497)
@@ -186,5 +186,5 @@
          ENDIF
          !
-         DO_2D( 1, 1, 1, 1 )
+         DO_2D( 1, 1, 1, 1 )                  ! vertical interpolation of T & S
             DO jk = 1, jpk                        ! determines the intepolated T-S profiles at each (i,j) points
                zl = gdept_0(ji,jj,jk)
Index: /NEMO/trunk/src/OCE/DYN/divhor.F90
===================================================================
--- /NEMO/trunk/src/OCE/DYN/divhor.F90	(revision 13496)
+++ /NEMO/trunk/src/OCE/DYN/divhor.F90	(revision 13497)
@@ -77,5 +77,5 @@
       ENDIF
       !
-      DO_3D( 0, 0, 0, 0, 1, jpkm1 )
+      DO_3D( 0, 0, 0, 0, 1, jpkm1 )                                    !==  Horizontal divergence  ==!
          hdiv(ji,jj,jk) = (  e2u(ji  ,jj) * e3u(ji  ,jj,jk,Kmm) * uu(ji  ,jj,jk,Kmm)      &
             &               - e2u(ji-1,jj) * e3u(ji-1,jj,jk,Kmm) * uu(ji-1,jj,jk,Kmm)      &
Index: /NEMO/trunk/src/OCE/DYN/dynadv_cen2.F90
===================================================================
--- /NEMO/trunk/src/OCE/DYN/dynadv_cen2.F90	(revision 13496)
+++ /NEMO/trunk/src/OCE/DYN/dynadv_cen2.F90	(revision 13497)
@@ -72,5 +72,5 @@
          zfu(:,:,jk) = 0.25_wp * e2u(:,:) * e3u(:,:,jk,Kmm) * puu(:,:,jk,Kmm)
          zfv(:,:,jk) = 0.25_wp * e1v(:,:) * e3v(:,:,jk,Kmm) * pvv(:,:,jk,Kmm)
-         DO_2D( 1, 0, 1, 0 )
+         DO_2D( 1, 0, 1, 0 )              ! horizontal momentum fluxes (at T- and F-point)
             zfu_t(ji+1,jj  ,jk) = ( zfu(ji,jj,jk) + zfu(ji+1,jj,jk) ) * ( puu(ji,jj,jk,Kmm) + puu(ji+1,jj  ,jk,Kmm) )
             zfv_f(ji  ,jj  ,jk) = ( zfv(ji,jj,jk) + zfv(ji+1,jj,jk) ) * ( puu(ji,jj,jk,Kmm) + puu(ji  ,jj+1,jk,Kmm) )
@@ -78,5 +78,5 @@
             zfv_t(ji  ,jj+1,jk) = ( zfv(ji,jj,jk) + zfv(ji,jj+1,jk) ) * ( pvv(ji,jj,jk,Kmm) + pvv(ji  ,jj+1,jk,Kmm) )
          END_2D
-         DO_2D( 0, 0, 0, 0 )
+         DO_2D( 0, 0, 0, 0 )              ! divergence of horizontal momentum fluxes
             puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) - (  zfu_t(ji+1,jj,jk) - zfu_t(ji,jj  ,jk)    &
                &                           + zfv_f(ji  ,jj,jk) - zfv_f(ji,jj-1,jk)  ) * r1_e1e2u(ji,jj)   &
@@ -98,5 +98,5 @@
       !                             !==  Vertical advection  ==!
       !
-      DO_2D( 0, 0, 0, 0 )
+      DO_2D( 0, 0, 0, 0 )                 ! surface/bottom advective fluxes set to zero
          zfu_uw(ji,jj,jpk) = 0._wp   ;   zfv_vw(ji,jj,jpk) = 0._wp
          zfu_uw(ji,jj, 1 ) = 0._wp   ;   zfv_vw(ji,jj, 1 ) = 0._wp
@@ -109,5 +109,5 @@
       ENDIF
       DO jk = 2, jpkm1                    ! interior advective fluxes
-         DO_2D( 0, 1, 0, 1 )
+         DO_2D( 0, 1, 0, 1 )                  ! 1/4 * Vertical transport
             zfw(ji,jj,jk) = 0.25_wp * e1e2t(ji,jj) * ww(ji,jj,jk)
          END_2D
@@ -117,5 +117,5 @@
          END_2D
       END DO
-      DO_3D( 0, 0, 0, 0, 1, jpkm1 )
+      DO_3D( 0, 0, 0, 0, 1, jpkm1 )       ! divergence of vertical momentum flux divergence
          puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) - ( zfu_uw(ji,jj,jk) - zfu_uw(ji,jj,jk+1) ) * r1_e1e2u(ji,jj)   &
             &                                      / e3u(ji,jj,jk,Kmm)
Index: /NEMO/trunk/src/OCE/DYN/dynadv_ubs.F90
===================================================================
--- /NEMO/trunk/src/OCE/DYN/dynadv_ubs.F90	(revision 13496)
+++ /NEMO/trunk/src/OCE/DYN/dynadv_ubs.F90	(revision 13497)
@@ -108,5 +108,5 @@
          zfv(:,:,jk) = e1v(:,:) * e3v(:,:,jk,Kmm) * pvv(:,:,jk,Kmm)
          !            
-         DO_2D( 0, 0, 0, 0 )
+         DO_2D( 0, 0, 0, 0 )                       ! laplacian
             zlu_uu(ji,jj,jk,1) = ( puu (ji+1,jj  ,jk,Kbb) - 2.*puu (ji,jj,jk,Kbb) + puu (ji-1,jj  ,jk,Kbb) ) * umask(ji,jj,jk)
             zlv_vv(ji,jj,jk,1) = ( pvv (ji  ,jj+1,jk,Kbb) - 2.*pvv (ji,jj,jk,Kbb) + pvv (ji  ,jj-1,jk,Kbb) ) * vmask(ji,jj,jk)
@@ -136,5 +136,5 @@
          zfv(:,:,jk) = 0.25_wp * e1v(:,:) * e3v(:,:,jk,Kmm) * pvv(:,:,jk,Kmm)
          !
-         DO_2D( 1, 0, 1, 0 )
+         DO_2D( 1, 0, 1, 0 )                       ! horizontal momentum fluxes at T- and F-point
             zui = ( puu(ji,jj,jk,Kmm) + puu(ji+1,jj  ,jk,Kmm) )
             zvj = ( pvv(ji,jj,jk,Kmm) + pvv(ji  ,jj+1,jk,Kmm) )
@@ -168,5 +168,5 @@
                &                * ( pvv(ji,jj,jk,Kmm) + pvv(ji+1,jj  ,jk,Kmm) - gamma1 * zl_v )
          END_2D
-         DO_2D( 0, 0, 0, 0 )
+         DO_2D( 0, 0, 0, 0 )                       ! divergence of horizontal momentum fluxes
             puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) - (  zfu_t(ji+1,jj,jk) - zfu_t(ji,jj  ,jk)    &
                &                           + zfv_f(ji  ,jj,jk) - zfv_f(ji,jj-1,jk)  ) * r1_e1e2u(ji,jj)   &
@@ -187,5 +187,5 @@
       !                                      !  Vertical advection  !
       !                                      ! ==================== !
-      DO_2D( 0, 0, 0, 0 )
+      DO_2D( 0, 0, 0, 0 )                          ! surface/bottom advective fluxes set to zero
          zfu_uw(ji,jj,jpk) = 0._wp
          zfv_vw(ji,jj,jpk) = 0._wp
@@ -208,5 +208,5 @@
          END_2D
       END DO
-      DO_3D( 0, 0, 0, 0, 1, jpkm1 )
+      DO_3D( 0, 0, 0, 0, 1, jpkm1 )             ! divergence of vertical momentum flux divergence
          puu(ji,jj,jk,Krhs) =  puu(ji,jj,jk,Krhs) - ( zfu_uw(ji,jj,jk) - zfu_uw(ji,jj,jk+1) ) * r1_e1e2u(ji,jj)   &
             &                                       / e3u(ji,jj,jk,Kmm)
Index: /NEMO/trunk/src/OCE/DYN/dynkeg.F90
===================================================================
--- /NEMO/trunk/src/OCE/DYN/dynkeg.F90	(revision 13496)
+++ /NEMO/trunk/src/OCE/DYN/dynkeg.F90	(revision 13497)
@@ -125,5 +125,5 @@
       END SELECT 
       !
-      DO_3D( 0, 0, 0, 0, 1, jpkm1 )
+      DO_3D( 0, 0, 0, 0, 1, jpkm1 )       !==  grad( KE ) added to the general momentum trends  ==!
          puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) - ( zhke(ji+1,jj  ,jk) - zhke(ji,jj,jk) ) / e1u(ji,jj)
          pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) - ( zhke(ji  ,jj+1,jk) - zhke(ji,jj,jk) ) / e2v(ji,jj)
Index: /NEMO/trunk/src/OCE/DYN/dynldf_iso.F90
===================================================================
--- /NEMO/trunk/src/OCE/DYN/dynldf_iso.F90	(revision 13496)
+++ /NEMO/trunk/src/OCE/DYN/dynldf_iso.F90	(revision 13497)
@@ -128,5 +128,5 @@
       IF( ln_dynldf_hor .AND. ln_traldf_iso ) THEN
          !
-         DO_3D( 0, 0, 0, 0, 1, jpk )
+         DO_3D( 0, 0, 0, 0, 1, jpk )      ! set the slopes of iso-level
             uslp (ji,jj,jk) = - ( gdept(ji+1,jj,jk,Kbb) - gdept(ji ,jj ,jk,Kbb) ) * r1_e1u(ji,jj) * umask(ji,jj,jk)
             vslp (ji,jj,jk) = - ( gdept(ji,jj+1,jk,Kbb) - gdept(ji ,jj ,jk,Kbb) ) * r1_e2v(ji,jj) * vmask(ji,jj,jk)
@@ -268,5 +268,5 @@
          ! Second derivative (divergence) and add to the general trend
          ! -----------------------------------------------------------
-         DO_2D( 0, 0, 0, 0 )
+         DO_2D( 0, 0, 0, 0 )      !!gm Question vectop possible??? !!bug
             puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + (  ziut(ji+1,jj) - ziut(ji,jj  )    &
                &                           + zjuf(ji  ,jj) - zjuf(ji,jj-1)  ) * r1_e1e2u(ji,jj)   &
Index: /NEMO/trunk/src/OCE/DYN/dynldf_lap_blp.F90
===================================================================
--- /NEMO/trunk/src/OCE/DYN/dynldf_lap_blp.F90	(revision 13496)
+++ /NEMO/trunk/src/OCE/DYN/dynldf_lap_blp.F90	(revision 13497)
@@ -84,5 +84,5 @@
          END_2D
          !
-         DO_2D( 0, 0, 0, 0 )
+         DO_2D( 0, 0, 0, 0 )                       ! - curl( curl) + grad( div )
             pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zsign * umask(ji,jj,jk) * (    &    ! * by umask is mandatory for dyn_ldf_blp use
                &              - ( zcur(ji  ,jj) - zcur(ji,jj-1) ) * r1_e2u(ji,jj) / e3u(ji,jj,jk,Kmm)   &
Index: /NEMO/trunk/src/OCE/DYN/dynspg.F90
===================================================================
--- /NEMO/trunk/src/OCE/DYN/dynspg.F90	(revision 13496)
+++ /NEMO/trunk/src/OCE/DYN/dynspg.F90	(revision 13497)
@@ -102,5 +102,5 @@
          IF( ln_apr_dyn .AND. .NOT.ln_dynspg_ts ) THEN   !==  Atmospheric pressure gradient (added later in time-split case) ==!
             zg_2 = grav * 0.5
-            DO_2D( 0, 0, 0, 0 )
+            DO_2D( 0, 0, 0, 0 )                       ! gradient of Patm using inverse barometer ssh
                spgu(ji,jj) = spgu(ji,jj) + zg_2 * (  ssh_ib (ji+1,jj) - ssh_ib (ji,jj)    &
                   &                                + ssh_ibb(ji+1,jj) - ssh_ibb(ji,jj)  ) * r1_e1u(ji,jj)
@@ -117,5 +117,5 @@
             CALL upd_tide(zt0step, Kmm)
             !
-            DO_2D( 0, 0, 0, 0 )
+            DO_2D( 0, 0, 0, 0 )                      ! add tide potential forcing
                spgu(ji,jj) = spgu(ji,jj) + grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) * r1_e1u(ji,jj)
                spgv(ji,jj) = spgv(ji,jj) + grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) * r1_e2v(ji,jj)
@@ -124,5 +124,5 @@
             IF (ln_scal_load) THEN
                zld = rn_scal_load * grav
-               DO_2D( 0, 0, 0, 0 )
+               DO_2D( 0, 0, 0, 0 )                   ! add scalar approximation for load potential
                   spgu(ji,jj) = spgu(ji,jj) + zld * ( pssh(ji+1,jj,Kmm) - pssh(ji,jj,Kmm) ) * r1_e1u(ji,jj)
                   spgv(ji,jj) = spgv(ji,jj) + zld * ( pssh(ji,jj+1,Kmm) - pssh(ji,jj,Kmm) ) * r1_e2v(ji,jj)
@@ -143,5 +143,5 @@
          ENDIF
          !
-         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
+         DO_3D( 0, 0, 0, 0, 1, jpkm1 )       !== Add all terms to the general trend
             puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + spgu(ji,jj)
             pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + spgv(ji,jj)
Index: /NEMO/trunk/src/OCE/DYN/dynspg_exp.F90
===================================================================
--- /NEMO/trunk/src/OCE/DYN/dynspg_exp.F90	(revision 13496)
+++ /NEMO/trunk/src/OCE/DYN/dynspg_exp.F90	(revision 13497)
@@ -74,10 +74,10 @@
       IF( ln_linssh ) THEN          !* linear free surface : add the surface pressure gradient trend
          !
-         DO_2D( 0, 0, 0, 0 )
+         DO_2D( 0, 0, 0, 0 )                 ! now surface pressure gradient
             spgu(ji,jj) = - grav * ( ssh(ji+1,jj,Kmm) - ssh(ji,jj,Kmm) ) * r1_e1u(ji,jj)
             spgv(ji,jj) = - grav * ( ssh(ji,jj+1,Kmm) - ssh(ji,jj,Kmm) ) * r1_e2v(ji,jj)
          END_2D
          !
-         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
+         DO_3D( 0, 0, 0, 0, 1, jpkm1 )       ! Add it to the general trend
             puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + spgu(ji,jj)
             pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + spgv(ji,jj)
Index: /NEMO/trunk/src/OCE/DYN/dynspg_ts.F90
===================================================================
--- /NEMO/trunk/src/OCE/DYN/dynspg_ts.F90	(revision 13496)
+++ /NEMO/trunk/src/OCE/DYN/dynspg_ts.F90	(revision 13497)
@@ -279,5 +279,5 @@
       ENDIF
       !
-      DO_2D( 0, 0, 0, 0 )
+      DO_2D( 0, 0, 0, 0 )                          ! Remove coriolis term (and possibly spg) from barotropic trend
           zu_frc(ji,jj) = zu_frc(ji,jj) - zu_trd(ji,jj) * ssumask(ji,jj)
           zv_frc(ji,jj) = zv_frc(ji,jj) - zv_trd(ji,jj) * ssvmask(ji,jj)
@@ -475,10 +475,10 @@
             !
             !                          ! ocean u- and v-depth at mid-step   (separate DO-loops remove the need of a lbc_lnk)
-            DO_2D( 1, 1, 1, 0 )
+            DO_2D( 1, 1, 1, 0 )   ! not jpi-column
                zhup2_e(ji,jj) = hu_0(ji,jj) + r1_2 * r1_e1e2u(ji,jj)                        &
                     &                              * (  e1e2t(ji  ,jj) * zsshp2_e(ji  ,jj)  &
                     &                                 + e1e2t(ji+1,jj) * zsshp2_e(ji+1,jj)  ) * ssumask(ji,jj)
             END_2D
-            DO_2D( 1, 0, 1, 1 )
+            DO_2D( 1, 0, 1, 1 )   ! not jpj-row
                zhvp2_e(ji,jj) = hv_0(ji,jj) + r1_2 * r1_e1e2v(ji,jj)                        &
                     &                              * (  e1e2t(ji,jj  ) * zsshp2_e(ji,jj  )  &
@@ -1308,5 +1308,5 @@
       !!----------------------------------------------------------------------
       !
-      DO_2D( 1, 1, 1, 0 )
+      DO_2D( 1, 1, 1, 0 )   ! not jpi-column
          IF ( phU(ji,jj) > 0._wp ) THEN   ;   pUmsk(ji,jj) = pTmsk(ji  ,jj) 
          ELSE                             ;   pUmsk(ji,jj) = pTmsk(ji+1,jj)  
@@ -1316,5 +1316,5 @@
       END_2D
       !
-      DO_2D( 1, 0, 1, 1 )
+      DO_2D( 1, 0, 1, 1 )   ! not jpj-row
          IF ( phV(ji,jj) > 0._wp ) THEN   ;   pVmsk(ji,jj) = pTmsk(ji,jj  )
          ELSE                             ;   pVmsk(ji,jj) = pTmsk(ji,jj+1)  
Index: /NEMO/trunk/src/OCE/DYN/dynvor.F90
===================================================================
--- /NEMO/trunk/src/OCE/DYN/dynvor.F90	(revision 13496)
+++ /NEMO/trunk/src/OCE/DYN/dynvor.F90	(revision 13497)
@@ -246,5 +246,5 @@
       CASE ( np_CRV )                           !* Coriolis + relative vorticity
          DO jk = 1, jpkm1                                 ! Horizontal slab
-            DO_2D( 1, 0, 1, 0 )
+            DO_2D( 1, 0, 1, 0 )                          ! relative vorticity
                zwz(ji,jj,jk) = (   e2v(ji+1,jj) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk)   &
                   &              - e1u(ji,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)   ) * r1_e1e2f(ji,jj)
Index: /NEMO/trunk/src/OCE/DYN/dynzad.F90
===================================================================
--- /NEMO/trunk/src/OCE/DYN/dynzad.F90	(revision 13496)
+++ /NEMO/trunk/src/OCE/DYN/dynzad.F90	(revision 13497)
@@ -71,5 +71,5 @@
       ENDIF
 
-      IF( l_trddyn )   THEN         ! Save puu(:,:,:,Krhs) and pvv(:,:,:,Krhs) trends
+      IF( l_trddyn )   THEN           ! Save puu(:,:,:,Krhs) and pvv(:,:,:,Krhs) trends
          ALLOCATE( ztrdu(jpi,jpj,jpk) , ztrdv(jpi,jpj,jpk) ) 
          ztrdu(:,:,:) = puu(:,:,:,Krhs) 
@@ -77,9 +77,9 @@
       ENDIF
       
-      DO jk = 2, jpkm1              ! Vertical momentum advection at level w and u- and v- vertical
-         DO_2D( 0, 1, 0, 1 )
+      DO jk = 2, jpkm1                ! Vertical momentum advection at level w and u- and v- vertical
+         DO_2D( 0, 1, 0, 1 )              ! vertical fluxes
             zww(ji,jj) = 0.25_wp * e1e2t(ji,jj) * ww(ji,jj,jk)
          END_2D
-         DO_2D( 0, 0, 0, 0 )
+         DO_2D( 0, 0, 0, 0 )              ! vertical momentum advection at w-point
             zwuw(ji,jj,jk) = ( zww(ji+1,jj  ) + zww(ji,jj) ) * ( puu(ji,jj,jk-1,Kmm) - puu(ji,jj,jk,Kmm) )
             zwvw(ji,jj,jk) = ( zww(ji  ,jj+1) + zww(ji,jj) ) * ( pvv(ji,jj,jk-1,Kmm) - pvv(ji,jj,jk,Kmm) )
@@ -95,5 +95,5 @@
       END_2D
       !
-      DO_3D( 0, 0, 0, 0, 1, jpkm1 )
+      DO_3D( 0, 0, 0, 0, 1, jpkm1 )   ! Vertical momentum advection at u- and v-points
          puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) - ( zwuw(ji,jj,jk) + zwuw(ji,jj,jk+1) ) * r1_e1e2u(ji,jj)   &
             &                                      / e3u(ji,jj,jk,Kmm)
@@ -102,5 +102,5 @@
       END_3D
 
-      IF( l_trddyn ) THEN           ! save the vertical advection trends for diagnostic
+      IF( l_trddyn ) THEN             ! save the vertical advection trends for diagnostic
          ztrdu(:,:,:) = puu(:,:,:,Krhs) - ztrdu(:,:,:)
          ztrdv(:,:,:) = pvv(:,:,:,Krhs) - ztrdv(:,:,:)
@@ -108,5 +108,5 @@
          DEALLOCATE( ztrdu, ztrdv ) 
       ENDIF
-      !                             ! Control print
+      !                               ! Control print
       IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab3d_1=puu(:,:,:,Krhs), clinfo1=' zad  - Ua: ', mask1=umask,   &
          &                                  tab3d_2=pvv(:,:,:,Krhs), clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
Index: /NEMO/trunk/src/OCE/DYN/dynzdf.F90
===================================================================
--- /NEMO/trunk/src/OCE/DYN/dynzdf.F90	(revision 13496)
+++ /NEMO/trunk/src/OCE/DYN/dynzdf.F90	(revision 13497)
@@ -131,5 +131,5 @@
             pvv(ji,jj,jk,Kaa) = ( pvv(ji,jj,jk,Kaa) - vv_b(ji,jj,Kaa) ) * vmask(ji,jj,jk)
          END_3D
-         DO_2D( 0, 0, 0, 0 )
+         DO_2D( 0, 0, 0, 0 )      ! Add bottom/top stress due to barotropic component only
             iku = mbku(ji,jj)         ! ocean bottom level at u- and v-points 
             ikv = mbkv(ji,jj)         ! (deepest ocean u- and v-points)
@@ -190,5 +190,5 @@
             END_3D
          END SELECT
-         DO_2D( 0, 0, 0, 0 )
+         DO_2D( 0, 0, 0, 0 )     !* Surface boundary conditions
             zwi(ji,jj,1) = 0._wp
             ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,1,Kmm)    &
@@ -227,5 +227,5 @@
             END_3D
          END SELECT
-         DO_2D( 0, 0, 0, 0 )
+         DO_2D( 0, 0, 0, 0 )     !* Surface boundary conditions
             zwi(ji,jj,1) = 0._wp
             zwd(ji,jj,1) = 1._wp - zws(ji,jj,1)
@@ -366,5 +366,5 @@
             END_3D
          END SELECT
-         DO_2D( 0, 0, 0, 0 )
+         DO_2D( 0, 0, 0, 0 )        !* Surface boundary conditions
             zwi(ji,jj,1) = 0._wp
             zwd(ji,jj,1) = 1._wp - zws(ji,jj,1)
Index: /NEMO/trunk/src/OCE/DYN/sshwzv.F90
===================================================================
--- /NEMO/trunk/src/OCE/DYN/sshwzv.F90	(revision 13496)
+++ /NEMO/trunk/src/OCE/DYN/sshwzv.F90	(revision 13497)
@@ -203,5 +203,5 @@
       ELSE                                            !==  Quasi-Eulerian vertical coordinate  ==!   ('key_qco')
          !                                            !==========================================!
-         DO jk = jpkm1, 1, -1                       ! integrate from the bottom the hor. divergence
+         DO jk = jpkm1, 1, -1                               ! integrate from the bottom the hor. divergence
             pww(:,:,jk) = pww(:,:,jk+1) - (  e3t(:,:,jk,Kmm) * hdiv(:,:,jk)                 &
                &                            + r1_Dt * (  e3t(:,:,jk,Kaa)        &
@@ -393,5 +393,5 @@
       !
       IF( MAXVAL( Cu_adv(:,:,:) ) > Cu_min ) THEN       ! Quick check if any breaches anywhere
-         DO_3DS( 1, 1, 1, 1, jpkm1, 2, -1 )
+         DO_3DS( 1, 1, 1, 1, jpkm1, 2, -1 )             ! or scan Courant criterion and partition ! w where necessary
             !
             zCu = MAX( Cu_adv(ji,jj,jk) , Cu_adv(ji,jj,jk-1) )
Index: /NEMO/trunk/src/OCE/DYN/wet_dry.F90
===================================================================
--- /NEMO/trunk/src/OCE/DYN/wet_dry.F90	(revision 13496)
+++ /NEMO/trunk/src/OCE/DYN/wet_dry.F90	(revision 13497)
@@ -307,5 +307,5 @@
       zwdlmtv(:,:) = 1._wp
       !
-      DO_2D( 0, 1, 0, 1 )
+      DO_2D( 0, 1, 0, 1 )      ! Horizontal Flux in u and v direction
          !
          IF( tmask(ji, jj, 1 ) < 0.5_wp) CYCLE   ! we don't care about land cells
Index: /NEMO/trunk/src/OCE/LDF/ldfc1d_c2d.F90
===================================================================
--- /NEMO/trunk/src/OCE/LDF/ldfc1d_c2d.F90	(revision 13496)
+++ /NEMO/trunk/src/OCE/LDF/ldfc1d_c2d.F90	(revision 13497)
@@ -80,5 +80,5 @@
             pah1(:,:,jk) = pahs1(:,:) * (  zratio + zc * ( 1._wp + TANH( - ( gdept_0(:,:,jk) - zh ) * zw) )  )
          END DO
-         DO_3DS( 1, 0, 1, 0, jpkm1, 1, -1 )
+         DO_3DS( 1, 0, 1, 0, jpkm1, 1, -1 )  ! pah2 at F-point (zdep2 is an approximation in zps-coord.)
             zdep2 = (  gdept_0(ji,jj+1,jk) + gdept_0(ji+1,jj+1,jk)   &
                &     + gdept_0(ji,jj  ,jk) + gdept_0(ji+1,jj  ,jk)  ) * r1_4
Index: /NEMO/trunk/src/OCE/LDF/ldfdyn.F90
===================================================================
--- /NEMO/trunk/src/OCE/LDF/ldfdyn.F90	(revision 13496)
+++ /NEMO/trunk/src/OCE/LDF/ldfdyn.F90	(revision 13497)
@@ -311,5 +311,5 @@
             IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'ldf_dyn_init: failed to allocate Smagorinsky arrays')
             !
-            DO_2D( 1, 1, 1, 1 )
+            DO_2D( 1, 1, 1, 1 )        ! Set local gridscale values
                esqt(ji,jj) = ( 2._wp * e1e2t(ji,jj) / ( e1t(ji,jj) + e2t(ji,jj) ) )**2 
                esqf(ji,jj) = ( 2._wp * e1e2f(ji,jj) / ( e1f(ji,jj) + e2f(ji,jj) ) )**2 
@@ -434,5 +434,5 @@
             DO jk = 1, jpkm1
               !
-               DO_2D( 0, 0, 0, 0 )
+               DO_2D( 0, 0, 0, 0 )                                   ! T-point value
                   !
                   zu2pv2_ij    = uu(ji  ,jj  ,jk,Kbb) * uu(ji  ,jj  ,jk,Kbb) + vv(ji  ,jj  ,jk,Kbb) * vv(ji  ,jj  ,jk,Kbb)
@@ -448,5 +448,5 @@
                END_2D
                !
-               DO_2D( 1, 0, 1, 0 )
+               DO_2D( 1, 0, 1, 0 )                                   ! F-point value
                   !
                   zu2pv2_ij_p1 = uu(ji  ,jj+1,jk, kbb) * uu(ji  ,jj+1,jk, kbb) + vv(ji+1,jj  ,jk, kbb) * vv(ji+1,jj  ,jk, kbb)
Index: /NEMO/trunk/src/OCE/LDF/ldfslp.F90
===================================================================
--- /NEMO/trunk/src/OCE/LDF/ldfslp.F90	(revision 13496)
+++ /NEMO/trunk/src/OCE/LDF/ldfslp.F90	(revision 13497)
@@ -128,5 +128,5 @@
       IF( ln_timing )   CALL timing_start('ldf_slp')
       !
-      zeps   =  1.e-20_wp        !==   Local constant initialization   ==!
+      zeps   =  1.e-20_wp           !==   Local constant initialization   ==!
       z1_16  =  1.0_wp / 16._wp
       zm1_g  = -1.0_wp / grav
@@ -137,5 +137,5 @@
       zwz(:,:,:) = 0._wp
       !
-      DO_3D( 1, 0, 1, 0, 1, jpk )
+      DO_3D( 1, 0, 1, 0, 1, jpk )   !==   i- & j-gradient of density   ==!
          zgru(ji,jj,jk) = umask(ji,jj,jk) * ( prd(ji+1,jj  ,jk) - prd(ji,jj,jk) )
          zgrv(ji,jj,jk) = vmask(ji,jj,jk) * ( prd(ji  ,jj+1,jk) - prd(ji,jj,jk) )
@@ -154,5 +154,5 @@
       ENDIF
       !
-      zdzr(:,:,1) = 0._wp        !==   Local vertical density gradient at T-point   == !   (evaluated from N^2)
+      zdzr(:,:,1) = 0._wp           !==   Local vertical density gradient at T-point   == !   (evaluated from N^2)
       DO jk = 2, jpkm1
          !                                ! zdzr = d/dz(prd)= - ( prd ) / grav * mk(pn2) -- at t point
@@ -165,5 +165,5 @@
       END DO
       !
-      !                          !==   Slopes just below the mixed layer   ==!
+      !                             !==   Slopes just below the mixed layer   ==!
       CALL ldf_slp_mxl( prd, pn2, zgru, zgrv, zdzr, Kmm )        ! output: uslpml, vslpml, wslpiml, wslpjml
 
@@ -186,5 +186,5 @@
       END IF
 
-      DO_3D( 0, 0, 0, 0, 2, jpkm1 )
+      DO_3D( 0, 0, 0, 0, 2, jpkm1 )        !* Slopes at u and v points
          !                                      ! horizontal and vertical density gradient at u- and v-points
          zau = zgru(ji,jj,jk) * r1_e1u(ji,jj)
@@ -231,7 +231,7 @@
       CALL lbc_lnk_multi( 'ldfslp', zwz, 'U', -1.0_wp,  zww, 'V', -1.0_wp )      ! lateral boundary conditions
       !
-      !                                            !* horizontal Shapiro filter
+      !                                    !* horizontal Shapiro filter
       DO jk = 2, jpkm1
-         DO_2D( 0, 0, 0, 0 )
+         DO_2D( 0, 0, 0, 0 )                                 ! rows jj=2 and =jpjm1 only
             uslp(ji,jj,jk) = z1_16 * (        zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk)      &
                &                       +      zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk)      &
@@ -245,5 +245,5 @@
                &                       + 4.*  zww(ji,jj    ,jk)                       )
          END_2D
-         DO jj = 3, jpj-2                               ! other rows
+         DO jj = 3, jpj-2                                    ! other rows
             DO ji = 2, jpim1   ! vector opt.
                uslp(ji,jj,jk) = z1_16 * (        zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk)      &
@@ -259,5 +259,5 @@
             END DO
          END DO
-         !                                        !* decrease along coastal boundaries
+         !                                 !* decrease along coastal boundaries
          DO_2D( 0, 0, 0, 0 )
             uslp(ji,jj,jk) = uslp(ji,jj,jk) * ( umask(ji,jj+1,jk) + umask(ji,jj-1,jk  ) ) * 0.5_wp   &
@@ -307,5 +307,5 @@
       !                                           !* horizontal Shapiro filter
       DO jk = 2, jpkm1
-         DO_2D( 0, 0, 0, 0 )
+         DO_2D( 0, 0, 0, 0 )                             ! rows jj=2 and =jpjm1 only
             zcofw = wmask(ji,jj,jk) * z1_16
             wslpi(ji,jj,jk) = (         zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk)     &
@@ -401,5 +401,5 @@
          !
          ip = jl   ;   jp = jl                ! guaranteed nonzero gradients ( absolute value larger than repsln)
-         DO_3D( 1, 0, 1, 0, 1, jpkm1 )
+         DO_3D( 1, 0, 1, 0, 1, jpkm1 )        ! done each pair of triad ! NB: not masked ==>  a minimum value is set
             zdit = ( ts(ji+1,jj,jk,jp_tem,Kbb) - ts(ji,jj,jk,jp_tem,Kbb) )    ! i-gradient of T & S at u-point
             zdis = ( ts(ji+1,jj,jk,jp_sal,Kbb) - ts(ji,jj,jk,jp_sal,Kbb) )
@@ -427,6 +427,6 @@
 
       DO kp = 0, 1                            !==  unmasked before density i- j-, k-gradients  ==!
-         DO_3D( 1, 1, 1, 1, 1, jpkm1 )
-            IF( jk+kp > 1 ) THEN        ! k-gradient of T & S a jk+kp
+         DO_3D( 1, 1, 1, 1, 1, jpkm1 )        ! done each pair of triad ! NB: not masked ==>  a minimum value is set
+            IF( jk+kp > 1 ) THEN              ! k-gradient of T & S a jk+kp
                zdkt = ( ts(ji,jj,jk+kp-1,jp_tem,Kbb) - ts(ji,jj,jk+kp,jp_tem,Kbb) )
                zdks = ( ts(ji,jj,jk+kp-1,jp_sal,Kbb) - ts(ji,jj,jk+kp,jp_sal,Kbb) )
@@ -442,5 +442,5 @@
       END DO
       !
-      DO_2D( 1, 1, 1, 1 )
+      DO_2D( 1, 1, 1, 1 )                     !==  Reciprocal depth of the w-point below ML base  ==!
          jk = MIN( nmln(ji,jj), mbkt(ji,jj) ) + 1     ! MIN in case ML depth is the ocean depth
          z1_mlbw(ji,jj) = 1._wp / gdepw(ji,jj,jk,Kmm)
@@ -628,5 +628,5 @@
       !
       !                                            !==   surface mixed layer mask   !
-      DO_3D( 1, 1, 1, 1, 1, jpk )
+      DO_3D( 1, 1, 1, 1, 1, jpk )                  ! =1 inside the mixed layer, =0 otherwise
          ik = nmln(ji,jj) - 1
          IF( jk <= ik ) THEN   ;   omlmask(ji,jj,jk) = 1._wp
Index: /NEMO/trunk/src/OCE/LDF/ldftra.F90
===================================================================
--- /NEMO/trunk/src/OCE/LDF/ldftra.F90	(revision 13496)
+++ /NEMO/trunk/src/OCE/LDF/ldftra.F90	(revision 13497)
@@ -694,5 +694,5 @@
       CALL lbc_lnk( 'ldftra', zaeiw(:,:), 'W', 1.0_wp )       ! lateral boundary condition
       !               
-      DO_2D( 0, 0, 0, 0 )
+      DO_2D( 0, 0, 0, 0 )                       !== aei at u- and v-points  ==!
          paeiu(ji,jj,1) = 0.5_wp * ( zaeiw(ji,jj) + zaeiw(ji+1,jj  ) ) * umask(ji,jj,1)
          paeiv(ji,jj,1) = 0.5_wp * ( zaeiw(ji,jj) + zaeiw(ji  ,jj+1) ) * vmask(ji,jj,1)
@@ -813,5 +813,5 @@
       CALL iom_put( "voce_eiv", zw3d )
       !
-      DO_3D( 0, 0, 0, 0, 1, jpkm1 )
+      DO_3D( 0, 0, 0, 0, 1, jpkm1 )                            ! e1 e2 w_eiv = dk[psix] + dk[psix]
          zw3d(ji,jj,jk) = (  psi_vw(ji,jj,jk) - psi_vw(ji  ,jj-1,jk)    &
             &              + psi_uw(ji,jj,jk) - psi_uw(ji-1,jj  ,jk)  ) / e1e2t(ji,jj)
Index: /NEMO/trunk/src/OCE/SBC/sbccpl.F90
===================================================================
--- /NEMO/trunk/src/OCE/SBC/sbccpl.F90	(revision 13496)
+++ /NEMO/trunk/src/OCE/SBC/sbccpl.F90	(revision 13497)
@@ -1195,5 +1195,5 @@
             !                              
             IF( srcv(jpr_otx1)%clgrid == 'T' ) THEN
-               DO_2D( 0, 0, 0, 0 )
+               DO_2D( 0, 0, 0, 0 )                                        ! T ==> (U,V)
                   frcv(jpr_otx1)%z3(ji,jj,1) = 0.5 * ( frcv(jpr_otx1)%z3(ji+1,jj  ,1) + frcv(jpr_otx1)%z3(ji,jj,1) )
                   frcv(jpr_oty1)%z3(ji,jj,1) = 0.5 * ( frcv(jpr_oty1)%z3(ji  ,jj+1,1) + frcv(jpr_oty1)%z3(ji,jj,1) )
@@ -1586,5 +1586,5 @@
             p_tauj(:,:) = frcv(jpr_ity1)%z3(:,:,1)
          CASE( 'T' )
-            DO_2D( 0, 0, 0, 0 )
+            DO_2D( 0, 0, 0, 0 )                    ! T ==> (U,V)
                ! take care of the land-sea mask to avoid "pollution" of coastal stress. p[uv]taui used in frazil and  rheology 
                zztmp1 = 0.5_wp * ( 2. - umask(ji,jj,1) ) * MAX( tmask(ji,jj,1),tmask(ji+1,jj  ,1) )
Index: /NEMO/trunk/src/OCE/SBC/sbcflx.F90
===================================================================
--- /NEMO/trunk/src/OCE/SBC/sbcflx.F90	(revision 13496)
+++ /NEMO/trunk/src/OCE/SBC/sbcflx.F90	(revision 13497)
@@ -133,5 +133,5 @@
             END_2D
          ENDIF
-         DO_2D( 0, 0, 0, 0 )
+         DO_2D( 0, 0, 0, 0 )                                      ! set the ocean fluxes from read fields
             utau(ji,jj) =   sf(jp_utau)%fnow(ji,jj,1)                              * umask(ji,jj,1)
             vtau(ji,jj) =   sf(jp_vtau)%fnow(ji,jj,1)                              * vmask(ji,jj,1)
Index: /NEMO/trunk/src/OCE/SBC/sbcrnf.F90
===================================================================
--- /NEMO/trunk/src/OCE/SBC/sbcrnf.F90	(revision 13496)
+++ /NEMO/trunk/src/OCE/SBC/sbcrnf.F90	(revision 13497)
@@ -215,7 +215,7 @@
             END_2D
          ELSE                    !* variable volume case
-            DO_2D( 1, 1, 1, 1 )
+            DO_2D( 1, 1, 1, 1 )              ! update the depth over which runoffs are distributed
                h_rnf(ji,jj) = 0._wp
-               DO jk = 1, nk_rnf(ji,jj)                           ! recalculates h_rnf to be the depth in metres
+               DO jk = 1, nk_rnf(ji,jj)                             ! recalculates h_rnf to be the depth in metres
                   h_rnf(ji,jj) = h_rnf(ji,jj) + e3t(ji,jj,jk,Kmm)   ! to the bottom of the relevant grid box
                END DO
@@ -374,5 +374,5 @@
             ENDIF
          END_2D
-         DO_2D( 1, 1, 1, 1 )
+         DO_2D( 1, 1, 1, 1 )                           ! set the associated depth
             h_rnf(ji,jj) = 0._wp
             DO jk = 1, nk_rnf(ji,jj)
@@ -404,5 +404,5 @@
          WHERE( zrnfcl(:,:,1) > 0._wp )  h_rnf(:,:) = zacoef * zrnfcl(:,:,1)   ! compute depth for all runoffs
          !
-         DO_2D( 1, 1, 1, 1 )
+         DO_2D( 1, 1, 1, 1 )                ! take in account min depth of ocean rn_hmin
             IF( zrnfcl(ji,jj,1) > 0._wp ) THEN
                jk = mbkt(ji,jj)
@@ -423,5 +423,5 @@
          END_2D
          !
-         DO_2D( 1, 1, 1, 1 )
+         DO_2D( 1, 1, 1, 1 )                          ! set the associated depth
             h_rnf(ji,jj) = 0._wp
             DO jk = 1, nk_rnf(ji,jj)
Index: /NEMO/trunk/src/OCE/SBC/sbcwave.F90
===================================================================
--- /NEMO/trunk/src/OCE/SBC/sbcwave.F90	(revision 13496)
+++ /NEMO/trunk/src/OCE/SBC/sbcwave.F90	(revision 13497)
@@ -121,5 +121,5 @@
             zk_t(ji,jj) = ABS( tsd2d(ji,jj) ) / MAX( ABS( 5.97_wp*ztransp ), 0.0000001_wp )
          END_2D
-         DO_2D( 1, 0, 1, 0 )
+         DO_2D( 1, 0, 1, 0 )          ! exp. wave number & Stokes drift velocity at u- & v-points
             zk_u(ji,jj) = 0.5_wp * ( zk_t(ji,jj) + zk_t(ji+1,jj) )
             zk_v(ji,jj) = 0.5_wp * ( zk_t(ji,jj) + zk_t(ji,jj+1) )
@@ -164,5 +164,5 @@
          zsqrtpi = SQRT(rpi)
          z_two_thirds = 2.0_wp / 3.0_wp
-         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
+         DO_3D( 0, 0, 0, 0, 1, jpkm1 )       ! exp. wave number & Stokes drift velocity at u- & v-points
             zbot_u = ( gdepw(ji,jj,jk+1,Kmm) + gdepw(ji+1,jj,jk+1,Kmm) )  ! 2 * bottom depth
             zbot_v = ( gdepw(ji,jj,jk+1,Kmm) + gdepw(ji,jj+1,jk+1,Kmm) )  ! 2 * bottom depth
@@ -204,5 +204,5 @@
       !                       !==  vertical Stokes Drift 3D velocity  ==!
       !
-      DO_3D( 0, 1, 0, 1, 1, jpkm1 )
+      DO_3D( 0, 1, 0, 1, 1, jpkm1 )    ! Horizontal e3*divergence
          ze3divh(ji,jj,jk) = (  e2u(ji  ,jj) * e3u(ji  ,jj,jk,Kmm) * usd(ji  ,jj,jk)    &
             &                 - e2u(ji-1,jj) * e3u(ji-1,jj,jk,Kmm) * usd(ji-1,jj,jk)    &
Index: /NEMO/trunk/src/OCE/TRA/eosbn2.F90
===================================================================
--- /NEMO/trunk/src/OCE/TRA/eosbn2.F90	(revision 13496)
+++ /NEMO/trunk/src/OCE/TRA/eosbn2.F90	(revision 13497)
@@ -873,5 +873,5 @@
       IF( ln_timing )   CALL timing_start('bn2')
       !
-      DO_3D( 1, 1, 1, 1, 2, jpkm1 )
+      DO_3D( 1, 1, 1, 1, 2, jpkm1 )      ! interior points only (2=< jk =< jpkm1 ); surface and bottom value set to zero one for all in istate.F90
          zrw =   ( gdepw(ji,jj,jk  ,Kmm) - gdept(ji,jj,jk,Kmm) )   &
             &  / ( gdept(ji,jj,jk-1,Kmm) - gdept(ji,jj,jk,Kmm) ) 
Index: /NEMO/trunk/src/OCE/TRA/traadv_cen.F90
===================================================================
--- /NEMO/trunk/src/OCE/TRA/traadv_cen.F90	(revision 13496)
+++ /NEMO/trunk/src/OCE/TRA/traadv_cen.F90	(revision 13497)
@@ -112,5 +112,5 @@
             ztu(:,:,jpk) = 0._wp                   ! Bottom value : flux set to zero
             ztv(:,:,jpk) = 0._wp
-            DO_3D( 0, 0, 0, 0, 1, jpkm1 )
+            DO_3D( 0, 0, 0, 0, 1, jpkm1 )          ! masked gradient
                ztu(ji,jj,jk) = ( pt(ji+1,jj  ,jk,jn,Kmm) - pt(ji,jj,jk,jn,Kmm) ) * umask(ji,jj,jk)
                ztv(ji,jj,jk) = ( pt(ji  ,jj+1,jk,jn,Kmm) - pt(ji,jj,jk,jn,Kmm) ) * vmask(ji,jj,jk)
@@ -118,5 +118,5 @@
             CALL lbc_lnk_multi( 'traadv_cen', ztu, 'U', -1.0_wp , ztv, 'V', -1.0_wp )   ! Lateral boundary cond.
             !
-            DO_3D( 0, 0, 0, 0, 1, jpkm1 )
+            DO_3D( 0, 0, 0, 0, 1, jpkm1 )           ! Horizontal advective fluxes
                zC2t_u = pt(ji,jj,jk,jn,Kmm) + pt(ji+1,jj  ,jk,jn,Kmm)   ! C2 interpolation of T at u- & v-points (x2)
                zC2t_v = pt(ji,jj,jk,jn,Kmm) + pt(ji  ,jj+1,jk,jn,Kmm)
@@ -159,5 +159,5 @@
          ENDIF
          !               
-         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
+         DO_3D( 0, 0, 0, 0, 1, jpkm1 )   !--  Divergence of advective fluxes  --!
             pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs)    &
                &             - (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )    &
@@ -166,5 +166,5 @@
                &                * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
          END_3D
-         !                             ! trend diagnostics
+         !                               ! trend diagnostics
          IF( l_trd ) THEN
             CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_xad, zwx, pU, pt(:,:,:,jn,Kmm) )
Index: /NEMO/trunk/src/OCE/TRA/traadv_fct.F90
===================================================================
--- /NEMO/trunk/src/OCE/TRA/traadv_fct.F90	(revision 13496)
+++ /NEMO/trunk/src/OCE/TRA/traadv_fct.F90	(revision 13497)
@@ -160,16 +160,16 @@
             zwy(ji,jj,jk) = 0.5 * ( zfp_vj * pt(ji,jj,jk,jn,Kbb) + zfm_vj * pt(ji  ,jj+1,jk,jn,Kbb) )
          END_3D
-         !                    !* upstream tracer flux in the k direction *!
-         DO_3D( 1, 1, 1, 1, 2, jpkm1 )
+         !                               !* upstream tracer flux in the k direction *!
+         DO_3D( 1, 1, 1, 1, 2, jpkm1 )      ! Interior value ( multiplied by wmask)
             zfp_wk = pW(ji,jj,jk) + ABS( pW(ji,jj,jk) )
             zfm_wk = pW(ji,jj,jk) - ABS( pW(ji,jj,jk) )
             zwz(ji,jj,jk) = 0.5 * ( zfp_wk * pt(ji,jj,jk,jn,Kbb) + zfm_wk * pt(ji,jj,jk-1,jn,Kbb) ) * wmask(ji,jj,jk)
          END_3D
-         IF( ln_linssh ) THEN    ! top ocean value (only in linear free surface as zwz has been w-masked)
-            IF( ln_isfcav ) THEN             ! top of the ice-shelf cavities and at the ocean surface
+         IF( ln_linssh ) THEN               ! top ocean value (only in linear free surface as zwz has been w-masked)
+            IF( ln_isfcav ) THEN                        ! top of the ice-shelf cavities and at the ocean surface
                DO_2D( 1, 1, 1, 1 )
                   zwz(ji,jj, mikt(ji,jj) ) = pW(ji,jj,mikt(ji,jj)) * pt(ji,jj,mikt(ji,jj),jn,Kbb)   ! linear free surface 
                END_2D
-            ELSE                             ! no cavities: only at the ocean surface
+            ELSE                                        ! no cavities: only at the ocean surface
                DO_2D( 1, 1, 1, 1 )
                   zwz(ji,jj,1) = pW(ji,jj,1) * pt(ji,jj,1,jn,Kbb)
@@ -178,10 +178,10 @@
          ENDIF
          !               
-         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
-            !                             ! total intermediate advective trends
+         DO_3D( 0, 0, 0, 0, 1, jpkm1 )   !* trend and after field with monotonic scheme
+            !                               ! total intermediate advective trends
             ztra = - (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   &
                &      + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  )   &
                &      + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1) ) * r1_e1e2t(ji,jj)
-            !                             ! update and guess with monotonic sheme
+            !                               ! update and guess with monotonic sheme
             pt(ji,jj,jk,jn,Krhs) =                   pt(ji,jj,jk,jn,Krhs) +       ztra   &
                &                                  / e3t(ji,jj,jk,Kmm ) * tmask(ji,jj,jk)
@@ -194,5 +194,5 @@
             !
             ztw(:,:,1) = 0._wp ; ztw(:,:,jpk) = 0._wp ;
-            DO_3D( 0, 0, 0, 0, 2, jpkm1 )
+            DO_3D( 0, 0, 0, 0, 2, jpkm1 )       ! Interior value ( multiplied by wmask)
                zfp_wk = wi(ji,jj,jk) + ABS( wi(ji,jj,jk) )
                zfm_wk = wi(ji,jj,jk) - ABS( wi(ji,jj,jk) )
@@ -227,9 +227,9 @@
             zltv(:,:,jpk) = 0._wp
             DO jk = 1, jpkm1                 ! Laplacian
-               DO_2D( 1, 0, 1, 0 )
+               DO_2D( 1, 0, 1, 0 )                 ! 1st derivative (gradient)
                   ztu(ji,jj,jk) = ( pt(ji+1,jj  ,jk,jn,Kmm) - pt(ji,jj,jk,jn,Kmm) ) * umask(ji,jj,jk)
                   ztv(ji,jj,jk) = ( pt(ji  ,jj+1,jk,jn,Kmm) - pt(ji,jj,jk,jn,Kmm) ) * vmask(ji,jj,jk)
                END_2D
-               DO_2D( 0, 0, 0, 0 )
+               DO_2D( 0, 0, 0, 0 )                 ! 2nd derivative * 1/ 6
                   zltu(ji,jj,jk) = (  ztu(ji,jj,jk) + ztu(ji-1,jj,jk)  ) * r1_6
                   zltv(ji,jj,jk) = (  ztv(ji,jj,jk) + ztv(ji,jj-1,jk)  ) * r1_6
@@ -238,8 +238,8 @@
             CALL lbc_lnk_multi( 'traadv_fct', zltu, 'T', 1.0_wp , zltv, 'T', 1.0_wp )   ! Lateral boundary cond. (unchanged sgn)
             !
-            DO_3D( 1, 0, 1, 0, 1, jpkm1 )
+            DO_3D( 1, 0, 1, 0, 1, jpkm1 )    ! Horizontal advective fluxes
                zC2t_u = pt(ji,jj,jk,jn,Kmm) + pt(ji+1,jj  ,jk,jn,Kmm)   ! 2 x C2 interpolation of T at u- & v-points
                zC2t_v = pt(ji,jj,jk,jn,Kmm) + pt(ji  ,jj+1,jk,jn,Kmm)
-               !                                                  ! C4 minus upstream advective fluxes 
+               !                                                        ! C4 minus upstream advective fluxes 
                zwx(ji,jj,jk) =  0.5_wp * pU(ji,jj,jk) * ( zC2t_u + zltu(ji,jj,jk) - zltu(ji+1,jj,jk) ) - zwx(ji,jj,jk)
                zwy(ji,jj,jk) =  0.5_wp * pV(ji,jj,jk) * ( zC2t_v + zltv(ji,jj,jk) - zltv(ji,jj+1,jk) ) - zwy(ji,jj,jk)
@@ -249,5 +249,5 @@
             ztu(:,:,jpk) = 0._wp             ! Bottom value : flux set to zero
             ztv(:,:,jpk) = 0._wp
-            DO_3D( 1, 0, 1, 0, 1, jpkm1 )
+            DO_3D( 1, 0, 1, 0, 1, jpkm1 )    ! 1st derivative (gradient)
                ztu(ji,jj,jk) = ( pt(ji+1,jj  ,jk,jn,Kmm) - pt(ji,jj,jk,jn,Kmm) ) * umask(ji,jj,jk)
                ztv(ji,jj,jk) = ( pt(ji  ,jj+1,jk,jn,Kmm) - pt(ji,jj,jk,jn,Kmm) ) * vmask(ji,jj,jk)
@@ -255,5 +255,5 @@
             CALL lbc_lnk_multi( 'traadv_fct', ztu, 'U', -1.0_wp , ztv, 'V', -1.0_wp )   ! Lateral boundary cond. (unchanged sgn)
             !
-            DO_3D( 0, 0, 0, 0, 1, jpkm1 )
+            DO_3D( 0, 0, 0, 0, 1, jpkm1 )    ! Horizontal advective fluxes
                zC2t_u = pt(ji,jj,jk,jn,Kmm) + pt(ji+1,jj  ,jk,jn,Kmm)   ! 2 x C2 interpolation of T at u- & v-points (x2)
                zC2t_v = pt(ji,jj,jk,jn,Kmm) + pt(ji  ,jj+1,jk,jn,Kmm)
@@ -288,6 +288,6 @@
          !         
          IF ( ll_zAimp ) THEN
-            DO_3D( 0, 0, 0, 0, 1, jpkm1 )
-               !                             ! total intermediate advective trends
+            DO_3D( 0, 0, 0, 0, 1, jpkm1 )    !* trend and after field with monotonic scheme
+               !                                                ! total intermediate advective trends
                ztra = - (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   &
                   &      + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  )   &
@@ -298,5 +298,5 @@
             CALL tridia_solver( zwdia, zwsup, zwinf, ztw, ztw , 0 )
             !
-            DO_3D( 0, 0, 0, 0, 2, jpkm1 )
+            DO_3D( 0, 0, 0, 0, 2, jpkm1 )       ! Interior value ( multiplied by wmask)
                zfp_wk = wi(ji,jj,jk) + ABS( wi(ji,jj,jk) )
                zfm_wk = wi(ji,jj,jk) - ABS( wi(ji,jj,jk) )
@@ -324,5 +324,5 @@
             !
             ztw(:,:,1) = 0._wp ; ztw(:,:,jpk) = 0._wp
-            DO_3D( 0, 0, 0, 0, 2, jpkm1 )
+            DO_3D( 0, 0, 0, 0, 2, jpkm1 )      ! Interior value ( multiplied by wmask)
                zfp_wk = wi(ji,jj,jk) + ABS( wi(ji,jj,jk) )
                zfm_wk = wi(ji,jj,jk) - ABS( wi(ji,jj,jk) )
@@ -454,6 +454,6 @@
          pbb(ji,jj,jk) = pbb(ji,jj,jk) * ( zcv * zav + ( 1._wp - zcv) * zbv )
 
-! monotonic flux in the k direction, i.e. pcc
-! -------------------------------------------
+      ! monotonic flux in the k direction, i.e. pcc
+      ! -------------------------------------------
          za = MIN( 1., zbetdo(ji,jj,jk+1), zbetup(ji,jj,jk) )
          zb = MIN( 1., zbetup(ji,jj,jk+1), zbetdo(ji,jj,jk) )
@@ -481,5 +481,5 @@
       !!----------------------------------------------------------------------
       
-      DO_3D( 1, 1, 1, 1, 3, jpkm1 )
+      DO_3D( 1, 1, 1, 1, 3, jpkm1 )       !==  build the three diagonal matrix  ==!
          zwd (ji,jj,jk) = 4._wp
          zwi (ji,jj,jk) = 1._wp
@@ -495,5 +495,5 @@
       END_3D
       !
-      jk = 2                                          ! Switch to second order centered at top
+      jk = 2                                    ! Switch to second order centered at top
       DO_2D( 1, 1, 1, 1 )
          zwd (ji,jj,jk) = 1._wp
@@ -504,5 +504,5 @@
       !
       !                       !==  tridiagonal solve  ==!
-      DO_2D( 1, 1, 1, 1 )
+      DO_2D( 1, 1, 1, 1 )           ! first recurrence
          zwt(ji,jj,2) = zwd(ji,jj,2)
       END_2D
@@ -511,5 +511,5 @@
       END_3D
       !
-      DO_2D( 1, 1, 1, 1 )
+      DO_2D( 1, 1, 1, 1 )           ! second recurrence:    Zk = Yk - Ik / Tk-1  Zk-1
          pt_out(ji,jj,2) = zwrm(ji,jj,2)
       END_2D
@@ -518,5 +518,5 @@
       END_3D
 
-      DO_2D( 1, 1, 1, 1 )
+      DO_2D( 1, 1, 1, 1 )           ! third recurrence: Xk = (Zk - Sk Xk+1 ) / Tk
          pt_out(ji,jj,jpkm1) = pt_out(ji,jj,jpkm1) / zwt(ji,jj,jpkm1)
       END_2D
@@ -546,5 +546,5 @@
       !                      !==  build the three diagonal matrix & the RHS  ==!
       !
-      DO_3D( 0, 0, 0, 0, 3, jpkm1 )
+      DO_3D( 0, 0, 0, 0, 3, jpkm1 )    ! interior (from jk=3 to jpk-1)
          zwd (ji,jj,jk) = 3._wp * wmask(ji,jj,jk) + 1._wp                 !       diagonal
          zwi (ji,jj,jk) =         wmask(ji,jj,jk)                         ! lower diagonal
@@ -565,5 +565,5 @@
       END IF
       !
-      DO_2D( 0, 0, 0, 0 )
+      DO_2D( 0, 0, 0, 0 )              ! 2nd order centered at top & bottom
          ikt = mikt(ji,jj) + 1            ! w-point below the 1st  wet point
          ikb = MAX(mbkt(ji,jj), 2)        !     -   above the last wet point
@@ -582,5 +582,5 @@
       !                       !==  tridiagonal solver  ==!
       !
-      DO_2D( 0, 0, 0, 0 )
+      DO_2D( 0, 0, 0, 0 )           !* 1st recurrence:   Tk = Dk - Ik Sk-1 / Tk-1
          zwt(ji,jj,2) = zwd(ji,jj,2)
       END_2D
@@ -589,5 +589,5 @@
       END_3D
       !
-      DO_2D( 0, 0, 0, 0 )
+      DO_2D( 0, 0, 0, 0 )           !* 2nd recurrence:    Zk = Yk - Ik / Tk-1  Zk-1
          pt_out(ji,jj,2) = zwrm(ji,jj,2)
       END_2D
@@ -596,5 +596,5 @@
       END_3D
 
-      DO_2D( 0, 0, 0, 0 )
+      DO_2D( 0, 0, 0, 0 )           !* 3d recurrence:    Xk = (Zk - Sk Xk+1 ) / Tk
          pt_out(ji,jj,jpkm1) = pt_out(ji,jj,jpkm1) / zwt(ji,jj,jpkm1)
       END_2D
@@ -638,5 +638,5 @@
       kstart =  1  + klev
       !
-      DO_2D( 0, 0, 0, 0 )
+      DO_2D( 0, 0, 0, 0 )                         !* 1st recurrence:   Tk = Dk - Ik Sk-1 / Tk-1
          zwt(ji,jj,kstart) = pD(ji,jj,kstart)
       END_2D
@@ -645,5 +645,5 @@
       END_3D
       !
-      DO_2D( 0, 0, 0, 0 )
+      DO_2D( 0, 0, 0, 0 )                        !* 2nd recurrence:    Zk = Yk - Ik / Tk-1  Zk-1
          pt_out(ji,jj,kstart) = pRHS(ji,jj,kstart)
       END_2D
@@ -652,5 +652,5 @@
       END_3D
 
-      DO_2D( 0, 0, 0, 0 )
+      DO_2D( 0, 0, 0, 0 )                       !* 3d recurrence:    Xk = (Zk - Sk Xk+1 ) / Tk
          pt_out(ji,jj,jpkm1) = pt_out(ji,jj,jpkm1) / zwt(ji,jj,jpkm1)
       END_2D
Index: /NEMO/trunk/src/OCE/TRA/traadv_mus.F90
===================================================================
--- /NEMO/trunk/src/OCE/TRA/traadv_mus.F90	(revision 13496)
+++ /NEMO/trunk/src/OCE/TRA/traadv_mus.F90	(revision 13497)
@@ -148,5 +148,5 @@
          END_3D
          !
-         DO_3D( 0, 1, 0, 1, 1, jpkm1 )
+         DO_3D( 0, 1, 0, 1, 1, jpkm1 )    !-- Slopes limitation
             zslpx(ji,jj,jk) = SIGN( 1.0_wp, zslpx(ji,jj,jk) ) * MIN(    ABS( zslpx(ji  ,jj,jk) ),   &
                &                                                     2.*ABS( zwx  (ji-1,jj,jk) ),   &
@@ -157,5 +157,5 @@
          END_3D
          !
-         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
+         DO_3D( 0, 0, 0, 0, 1, jpkm1 )    !-- MUSCL horizontal advective fluxes
             ! MUSCL fluxes
             z0u = SIGN( 0.5_wp, pU(ji,jj,jk) )
@@ -175,5 +175,5 @@
          CALL lbc_lnk_multi( 'traadv_mus', zwx, 'U', -1.0_wp , zwy, 'V', -1.0_wp )   ! lateral boundary conditions   (changed sign)
          !
-         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
+         DO_3D( 0, 0, 0, 0, 1, jpkm1 )    !-- Tracer advective trend
             pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - ( zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )       &
             &                                     + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  ) )     &
@@ -204,10 +204,10 @@
                &            * (  0.25 + SIGN( 0.25_wp, zwx(ji,jj,jk) * zwx(ji,jj,jk+1) )  )
          END_3D
-         DO_3D( 1, 1, 1, 1, 2, jpkm1 )
+         DO_3D( 1, 1, 1, 1, 2, jpkm1 )    !-- Slopes limitation
             zslpx(ji,jj,jk) = SIGN( 1.0_wp, zslpx(ji,jj,jk) ) * MIN(    ABS( zslpx(ji,jj,jk  ) ),   &
                &                                                     2.*ABS( zwx  (ji,jj,jk+1) ),   &
                &                                                     2.*ABS( zwx  (ji,jj,jk  ) )  )
          END_3D
-         DO_3D( 0, 0, 0, 0, 1, jpk-2 )
+         DO_3D( 0, 0, 0, 0, 1, jpk-2 )    !-- vertical advective flux
             z0w = SIGN( 0.5_wp, pW(ji,jj,jk+1) )
             zalpha = 0.5 + z0w
@@ -227,5 +227,5 @@
          ENDIF
          !
-         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
+         DO_3D( 0, 0, 0, 0, 1, jpkm1 )     !-- vertical advective trend
             pt(ji,jj,jk,jn,Krhs) =  pt(ji,jj,jk,jn,Krhs) - ( zwx(ji,jj,jk) - zwx(ji,jj,jk+1) )   &
                &                                      * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
Index: /NEMO/trunk/src/OCE/TRA/traadv_qck.F90
===================================================================
--- /NEMO/trunk/src/OCE/TRA/traadv_qck.F90	(revision 13496)
+++ /NEMO/trunk/src/OCE/TRA/traadv_qck.F90	(revision 13497)
@@ -142,5 +142,5 @@
          !
 !!gm why not using a SHIFT instruction...
-         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
+         DO_3D( 0, 0, 0, 0, 1, jpkm1 )     !--- Computation of the ustream and downstream value of the tracer and the mask
             zfc(ji,jj,jk) = pt(ji-1,jj,jk,jn,Kbb)        ! Upstream   in the x-direction for the tracer
             zfd(ji,jj,jk) = pt(ji+1,jj,jk,jn,Kbb)        ! Downstream in the x-direction for the tracer
@@ -327,5 +327,5 @@
          !                                                       ! ===========
          !
-         DO_3D( 0, 0, 0, 0, 2, jpkm1 )
+         DO_3D( 0, 0, 0, 0, 2, jpkm1 )       !* Interior point   (w-masked 2nd order centered flux)
             zwz(ji,jj,jk) = 0.5 * pW(ji,jj,jk) * ( pt(ji,jj,jk-1,jn,Kmm) + pt(ji,jj,jk,jn,Kmm) ) * wmask(ji,jj,jk)
          END_3D
@@ -340,5 +340,5 @@
          ENDIF
          !
-         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
+         DO_3D( 0, 0, 0, 0, 1, jpkm1 )   !==  Tracer flux divergence added to the general trend  ==!
             pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - ( zwz(ji,jj,jk) - zwz(ji,jj,jk+1) )   &
                &                                * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
Index: /NEMO/trunk/src/OCE/TRA/traadv_ubs.F90
===================================================================
--- /NEMO/trunk/src/OCE/TRA/traadv_ubs.F90	(revision 13496)
+++ /NEMO/trunk/src/OCE/TRA/traadv_ubs.F90	(revision 13497)
@@ -124,6 +124,6 @@
          !                                                       ! ===========
          !                                              
-         DO jk = 1, jpkm1        !==  horizontal laplacian of before tracer ==!
-            DO_2D( 1, 0, 1, 0 )
+         DO jk = 1, jpkm1                !==  horizontal laplacian of before tracer ==!
+            DO_2D( 1, 0, 1, 0 )                   ! First derivative (masked gradient)
                zeeu = e2_e1u(ji,jj) * e3u(ji,jj,jk,Kmm) * umask(ji,jj,jk)
                zeev = e1_e2v(ji,jj) * e3v(ji,jj,jk,Kmm) * vmask(ji,jj,jk)
@@ -131,5 +131,5 @@
                ztv(ji,jj,jk) = zeev * ( pt(ji  ,jj+1,jk,jn,Kbb) - pt(ji,jj,jk,jn,Kbb) )
             END_2D
-            DO_2D( 0, 0, 0, 0 )
+            DO_2D( 0, 0, 0, 0 )                   ! Second derivative (divergence)
                zcoef = 1._wp / ( 6._wp * e3t(ji,jj,jk,Kmm) )
                zltu(ji,jj,jk) = (  ztu(ji,jj,jk) - ztu(ji-1,jj,jk)  ) * zcoef
@@ -140,6 +140,6 @@
          CALL lbc_lnk( 'traadv_ubs', zltu, 'T', 1.0_wp )   ;    CALL lbc_lnk( 'traadv_ubs', zltv, 'T', 1.0_wp )   ! Lateral boundary cond. (unchanged sgn)
          !    
-         DO_3D( 1, 0, 1, 0, 1, jpkm1 )
-            zfp_ui = pU(ji,jj,jk) + ABS( pU(ji,jj,jk) )      ! upstream transport (x2)
+         DO_3D( 1, 0, 1, 0, 1, jpkm1 )   !==  Horizontal advective fluxes  ==!     (UBS)
+            zfp_ui = pU(ji,jj,jk) + ABS( pU(ji,jj,jk) )        ! upstream transport (x2)
             zfm_ui = pU(ji,jj,jk) - ABS( pU(ji,jj,jk) )
             zfp_vj = pV(ji,jj,jk) + ABS( pV(ji,jj,jk) )
@@ -166,5 +166,5 @@
          !
          zltu(:,:,:) = pt(:,:,:,jn,Krhs) - zltu(:,:,:)    ! Horizontal advective trend used in vertical 2nd order FCT case
-         !                                            ! and/or in trend diagnostic (l_trd=T) 
+         !                                                ! and/or in trend diagnostic (l_trd=T) 
          !                
          IF( l_trd ) THEN                  ! trend diagnostics
@@ -187,5 +187,5 @@
             IF( l_trd )   zltv(:,:,:) = pt(:,:,:,jn,Krhs)          ! store pt(:,:,:,:,Krhs) if trend diag.
             !
-            !                          !*  upstream advection with initial mass fluxes & intermediate update  ==!
+            !                               !*  upstream advection with initial mass fluxes & intermediate update  ==!
             DO_3D( 1, 1, 1, 1, 2, jpkm1 )
                zfp_wk = pW(ji,jj,jk) + ABS( pW(ji,jj,jk) )
@@ -193,15 +193,15 @@
                ztw(ji,jj,jk) = 0.5_wp * (  zfp_wk * pt(ji,jj,jk,jn,Kbb) + zfm_wk * pt(ji,jj,jk-1,jn,Kbb)  ) * wmask(ji,jj,jk)
             END_3D
-            IF( ln_linssh ) THEN             ! top ocean value (only in linear free surface as ztw has been w-masked)
-               IF( ln_isfcav ) THEN                ! top of the ice-shelf cavities and at the ocean surface
+            IF( ln_linssh ) THEN                ! top ocean value (only in linear free surface as ztw has been w-masked)
+               IF( ln_isfcav ) THEN                   ! top of the ice-shelf cavities and at the ocean surface
                   DO_2D( 1, 1, 1, 1 )
                      ztw(ji,jj, mikt(ji,jj) ) = pW(ji,jj,mikt(ji,jj)) * pt(ji,jj,mikt(ji,jj),jn,Kbb)   ! linear free surface 
                   END_2D
-               ELSE                                ! no cavities: only at the ocean surface
+               ELSE                                   ! no cavities: only at the ocean surface
                   ztw(:,:,1) = pW(:,:,1) * pt(:,:,1,jn,Kbb)
                ENDIF
             ENDIF
             !
-            DO_3D( 0, 0, 0, 0, 1, jpkm1 )
+            DO_3D( 0, 0, 0, 0, 1, jpkm1 )   !* trend and after field with monotonic scheme
                ztak = - ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) )    &
                   &     * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
@@ -230,11 +230,11 @@
          END SELECT
          !
-         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
+         DO_3D( 0, 0, 0, 0, 1, jpkm1 )   !  final trend with corrected fluxes
             pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) )    &
                &                                        * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
          END_3D
          !
-         IF( l_trd )  THEN       ! vertical advective trend diagnostics
-            DO_3D( 0, 0, 0, 0, 1, jpkm1 )
+         IF( l_trd )  THEN               ! vertical advective trend diagnostics
+            DO_3D( 0, 0, 0, 0, 1, jpkm1 )                 ! (compute -w.dk[ptn]= -dk[w.ptn] + ptn.dk[w])
                zltv(ji,jj,jk) = pt(ji,jj,jk,jn,Krhs) - zltv(ji,jj,jk)                          &
                   &           + pt(ji,jj,jk,jn,Kmm) * (  pW(ji,jj,jk) - pW(ji,jj,jk+1)  )   &
Index: /NEMO/trunk/src/OCE/TRA/trabbl.F90
===================================================================
--- /NEMO/trunk/src/OCE/TRA/trabbl.F90	(revision 13496)
+++ /NEMO/trunk/src/OCE/TRA/trabbl.F90	(revision 13497)
@@ -197,5 +197,5 @@
          END_2D
          !               
-         DO_2D( 0, 0, 0, 0 )
+         DO_2D( 0, 0, 0, 0 )                               ! Compute the trend
             ik = mbkt(ji,jj)                            ! bottom T-level index
             pt_rhs(ji,jj,ik,jn) = pt_rhs(ji,jj,ik,jn)                                                  &
@@ -358,5 +358,5 @@
       IF( nn_bbl_ldf == 1 ) THEN          !   diffusive bbl   !
          !                                !-------------------!
-         DO_2D( 1, 0, 1, 0 )
+         DO_2D( 1, 0, 1, 0 )                   ! (criteria for non zero flux: grad(rho).grad(h) < 0 )
             !                                                   ! i-direction
             za = zab(ji+1,jj,jp_tem) + zab(ji,jj,jp_tem)              ! 2*(alpha,beta) at u-point
@@ -388,5 +388,5 @@
          !
          CASE( 1 )                                   != use of upper velocity
-            DO_2D( 1, 0, 1, 0 )
+            DO_2D( 1, 0, 1, 0 )                              ! criteria: grad(rho).grad(h)<0  and grad(rho).grad(h)<0
                !                                                  ! i-direction
                za = zab(ji+1,jj,jp_tem) + zab(ji,jj,jp_tem)               ! 2*(alpha,beta) at u-point
@@ -417,5 +417,5 @@
          CASE( 2 )                                 != bbl velocity = F( delta rho )
             zgbbl = grav * rn_gambbl
-            DO_2D( 1, 0, 1, 0 )
+            DO_2D( 1, 0, 1, 0 )                         ! criteria: rho_up > rho_down
                !                                                  ! i-direction
                ! down-slope T-point i/k-index (deep)  &   up-slope T-point i/k-index (shelf)
@@ -509,5 +509,5 @@
       !
       !                             !* vertical index of  "deep" bottom u- and v-points
-      DO_2D( 1, 0, 1, 0 )
+      DO_2D( 1, 0, 1, 0 )                 ! (the "shelf" bottom k-indices are mbku and mbkv)
          mbku_d(ji,jj) = MAX(  mbkt(ji+1,jj  ) , mbkt(ji,jj)  )   ! >= 1 as mbkt=1 over land
          mbkv_d(ji,jj) = MAX(  mbkt(ji  ,jj+1) , mbkt(ji,jj)  )
@@ -530,5 +530,5 @@
       END_2D
       !
-      DO_2D( 1, 0, 1, 0 )
+      DO_2D( 1, 0, 1, 0 )           !* bbl thickness at u- (v-) point; minimum of top & bottom e3u_0 (e3v_0)
          e3u_bbl_0(ji,jj) = MIN( e3u_0(ji,jj,mbkt(ji+1,jj  )), e3u_0(ji,jj,mbkt(ji,jj)) )
          e3v_bbl_0(ji,jj) = MIN( e3v_0(ji,jj,mbkt(ji  ,jj+1)), e3v_0(ji,jj,mbkt(ji,jj)) )
Index: /NEMO/trunk/src/OCE/TRA/traldf_iso.F90
===================================================================
--- /NEMO/trunk/src/OCE/TRA/traldf_iso.F90	(revision 13496)
+++ /NEMO/trunk/src/OCE/TRA/traldf_iso.F90	(revision 13497)
@@ -205,5 +205,5 @@
          END_3D
          IF( ln_zps ) THEN      ! botton and surface ocean correction of the horizontal gradient
-            DO_2D( 1, 0, 1, 0 )
+            DO_2D( 1, 0, 1, 0 )           ! bottom correction (partial bottom cell)
                zdit(ji,jj,mbku(ji,jj)) = pgu(ji,jj,jn)          
                zdjt(ji,jj,mbkv(ji,jj)) = pgv(ji,jj,jn)
@@ -229,5 +229,5 @@
             ELSE                 ;   zdkt(:,:) = ( pt(:,:,jk-1,jn) - pt(:,:,jk,jn) ) * wmask(:,:,jk)
             ENDIF
-            DO_2D( 1, 0, 1, 0 )
+            DO_2D( 1, 0, 1, 0 )           !==  Horizontal fluxes
                zabe1 = pahu(ji,jj,jk) * e2_e1u(ji,jj) * e3u(ji,jj,jk,Kmm)
                zabe2 = pahv(ji,jj,jk) * e1_e2v(ji,jj) * e3v(ji,jj,jk,Kmm)
@@ -250,5 +250,5 @@
             END_2D
             !
-            DO_2D( 0, 0, 0, 0 )
+            DO_2D( 0, 0, 0, 0 )           !== horizontal divergence and add to pta
                pt_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn)    &
                   &       + zsign * (  zftu(ji,jj,jk) - zftu(ji-1,jj,jk) + zftv(ji,jj,jk) - zftv(ji,jj-1,jk)  )   &
@@ -266,5 +266,5 @@
          ztfw(:,:, 1 ) = 0._wp      ;      ztfw(:,:,jpk) = 0._wp
          
-         DO_3D( 0, 0, 0, 0, 2, jpkm1 )
+         DO_3D( 0, 0, 0, 0, 2, jpkm1 )    ! interior (2=<jk=<jpk-1)
             !
             zmsku = wmask(ji,jj,jk) / MAX(   umask(ji  ,jj,jk-1) + umask(ji-1,jj,jk)          &
@@ -311,5 +311,5 @@
          ENDIF
          !         
-         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
+         DO_3D( 0, 0, 0, 0, 1, jpkm1 )    !==  Divergence of vertical fluxes added to pta  ==!
             pt_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn) + zsign * (  ztfw (ji,jj,jk) - ztfw(ji,jj,jk+1)  ) * r1_e1e2t(ji,jj)   &
                &                                             / e3t(ji,jj,jk,Kmm)
Index: /NEMO/trunk/src/OCE/TRA/traldf_lap_blp.F90
===================================================================
--- /NEMO/trunk/src/OCE/TRA/traldf_lap_blp.F90	(revision 13496)
+++ /NEMO/trunk/src/OCE/TRA/traldf_lap_blp.F90	(revision 13497)
@@ -108,14 +108,14 @@
          !                          ! =========== !    
          !                               
-         DO_3D( 1, 0, 1, 0, 1, jpkm1 )
+         DO_3D( 1, 0, 1, 0, 1, jpkm1 )            !== First derivative (gradient)  ==!
             ztu(ji,jj,jk) = zaheeu(ji,jj,jk) * ( pt(ji+1,jj  ,jk,jn) - pt(ji,jj,jk,jn) )
             ztv(ji,jj,jk) = zaheev(ji,jj,jk) * ( pt(ji  ,jj+1,jk,jn) - pt(ji,jj,jk,jn) )
          END_3D
-         IF( ln_zps ) THEN                ! set gradient at bottom/top ocean level
-            DO_2D( 1, 0, 1, 0 )
+         IF( ln_zps ) THEN                             ! set gradient at bottom/top ocean level
+            DO_2D( 1, 0, 1, 0 )                              ! bottom
                ztu(ji,jj,mbku(ji,jj)) = zaheeu(ji,jj,mbku(ji,jj)) * pgu(ji,jj,jn)
                ztv(ji,jj,mbkv(ji,jj)) = zaheev(ji,jj,mbkv(ji,jj)) * pgv(ji,jj,jn)
             END_2D
-            IF( ln_isfcav ) THEN                ! top in ocean cavities only
+            IF( ln_isfcav ) THEN                             ! top in ocean cavities only
                DO_2D( 1, 0, 1, 0 )
                   IF( miku(ji,jj) > 1 )   ztu(ji,jj,miku(ji,jj)) = zaheeu(ji,jj,miku(ji,jj)) * pgui(ji,jj,jn) 
@@ -125,5 +125,5 @@
          ENDIF
          !
-         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
+         DO_3D( 0, 0, 0, 0, 1, jpkm1 )            !== Second derivative (divergence) added to the general tracer trends  ==!
             pt_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn) + (  ztu(ji,jj,jk) - ztu(ji-1,jj,jk)     &
                &                                      +    ztv(ji,jj,jk) - ztv(ji,jj-1,jk) )   &
Index: /NEMO/trunk/src/OCE/TRA/traldf_triad.F90
===================================================================
--- /NEMO/trunk/src/OCE/TRA/traldf_triad.F90	(revision 13496)
+++ /NEMO/trunk/src/OCE/TRA/traldf_triad.F90	(revision 13497)
@@ -211,10 +211,10 @@
          zftv(:,:,:) = 0._wp
          !
-         DO_3D( 1, 0, 1, 0, 1, jpkm1 )
+         DO_3D( 1, 0, 1, 0, 1, jpkm1 )    !==  before lateral T & S gradients at T-level jk  ==!
             zdit(ji,jj,jk) = ( pt(ji+1,jj  ,jk,jn) - pt(ji,jj,jk,jn) ) * umask(ji,jj,jk)
             zdjt(ji,jj,jk) = ( pt(ji  ,jj+1,jk,jn) - pt(ji,jj,jk,jn) ) * vmask(ji,jj,jk)
          END_3D
          IF( ln_zps .AND. l_grad_zps ) THEN    ! partial steps: correction at top/bottom ocean level
-            DO_2D( 1, 0, 1, 0 )
+            DO_2D( 1, 0, 1, 0 )                    ! bottom level
                zdit(ji,jj,mbku(ji,jj)) = pgu(ji,jj,jn)
                zdjt(ji,jj,mbkv(ji,jj)) = pgv(ji,jj,jn)
@@ -361,5 +361,5 @@
          ENDIF
          !
-         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
+         DO_3D( 0, 0, 0, 0, 1, jpkm1 )      !==  Divergence of vertical fluxes added to pta  ==!
             pt_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn)    &
             &                                  + zsign * (  ztfw(ji,jj,jk+1) - ztfw(ji,jj,jk)  )   &
Index: /NEMO/trunk/src/OCE/TRA/tramle.F90
===================================================================
--- /NEMO/trunk/src/OCE/TRA/tramle.F90	(revision 13496)
+++ /NEMO/trunk/src/OCE/TRA/tramle.F90	(revision 13497)
@@ -100,5 +100,5 @@
       inml_mle(:,:) = mbkt(:,:) + 1                    ! init. to number of ocean w-level (T-level + 1)
       IF ( nla10 > 0 ) THEN                            ! avoid case where first level is thicker than 10m
-         DO_3DS( 1, 1, 1, 1, jpkm1, nlb10, -1 )
+         DO_3DS( 1, 1, 1, 1, jpkm1, nlb10, -1 )        ! from the bottom to nlb10 (10m)
             IF( rhop(ji,jj,jk) > rhop(ji,jj,nla10) + rn_rho_c_mle )   inml_mle(ji,jj) = jk      ! Mixed layer
          END_3D
@@ -110,5 +110,5 @@
       zbm (:,:) = 0._wp
       zn2 (:,:) = 0._wp
-      DO_3D( 1, 1, 1, 1, 1, ikmax )
+      DO_3D( 1, 1, 1, 1, 1, ikmax )                    ! MLD and mean buoyancy and N2 over the mixed layer
          zc = e3t(ji,jj,jk,Kmm) * REAL( MIN( MAX( 0, inml_mle(ji,jj)-jk ) , 1  )  )    ! zc being 0 outside the ML t-points
          zmld(ji,jj) = zmld(ji,jj) + zc
@@ -182,5 +182,5 @@
       zpsi_vw(:,:,:) = 0._wp
       !
-      DO_3D( 1, 0, 1, 0, 2, ikmax )
+      DO_3D( 1, 0, 1, 0, 2, ikmax )                ! start from 2 : surface value = 0
          zcuw = 1._wp - ( gdepw(ji+1,jj,jk,Kmm) + gdepw(ji,jj,jk,Kmm) ) * zhu(ji,jj)
          zcvw = 1._wp - ( gdepw(ji,jj+1,jk,Kmm) + gdepw(ji,jj,jk,Kmm) ) * zhv(ji,jj)
@@ -196,5 +196,5 @@
       !                                      !==  transport increased by the MLE induced transport ==!
       DO jk = 1, ikmax
-         DO_2D( 1, 0, 1, 0 )
+         DO_2D( 1, 0, 1, 0 )                      ! CAUTION pu,pv must be defined at row/column i=1 / j=1
             pu(ji,jj,jk) = pu(ji,jj,jk) + ( zpsi_uw(ji,jj,jk) - zpsi_uw(ji,jj,jk+1) )
             pv(ji,jj,jk) = pv(ji,jj,jk) + ( zpsi_vw(ji,jj,jk) - zpsi_vw(ji,jj,jk+1) )
@@ -283,5 +283,5 @@
             IF( ierr /= 0 )   CALL ctl_stop( 'tra_adv_mle_init: failed to allocate arrays' )
             z1_t2 = 1._wp / ( rn_time * rn_time )
-            DO_2D( 0, 1, 0, 1 )
+            DO_2D( 0, 1, 0, 1 )                      ! "coriolis+ time^-1" at u- & v-points
                zfu = ( ff_f(ji,jj) + ff_f(ji,jj-1) ) * 0.5_wp
                zfv = ( ff_f(ji,jj) + ff_f(ji-1,jj) ) * 0.5_wp
Index: /NEMO/trunk/src/OCE/TRA/tranpc.F90
===================================================================
--- /NEMO/trunk/src/OCE/TRA/tranpc.F90	(revision 13496)
+++ /NEMO/trunk/src/OCE/TRA/tranpc.F90	(revision 13497)
@@ -103,5 +103,5 @@
          inpcc = 0
          !
-         DO_2D( 0, 0, 0, 0 )
+         DO_2D( 0, 0, 0, 0 )                                ! interior column only
             !
             IF( tmask(ji,jj,2) == 1 ) THEN      ! At least 2 ocean points
Index: /NEMO/trunk/src/OCE/TRA/traqsr.F90
===================================================================
--- /NEMO/trunk/src/OCE/TRA/traqsr.F90	(revision 13496)
+++ /NEMO/trunk/src/OCE/TRA/traqsr.F90	(revision 13497)
@@ -231,5 +231,5 @@
          END_2D
          !
-         !* interior equi-partition in R-G-B depending on vertical profile of Chl
+         !                                    !* interior equi-partition in R-G-B depending on vertical profile of Chl
          DO_3D( 0, 0, 0, 0, 2, nksr + 1 )
             ze3t = e3t(ji,jj,jk-1,Kmm)
@@ -246,5 +246,5 @@
          END_3D
          !
-         DO_3D( 0, 0, 0, 0, 1, nksr )
+         DO_3D( 0, 0, 0, 0, 1, nksr )          !* now qsr induced heat content
             qsr_hc(ji,jj,jk) = r1_rho0_rcp * ( ztmp3d(ji,jj,jk) - ztmp3d(ji,jj,jk+1) )
          END_3D
@@ -256,5 +256,5 @@
          zz0 =        rn_abs   * r1_rho0_rcp      ! surface equi-partition in 2-bands
          zz1 = ( 1. - rn_abs ) * r1_rho0_rcp
-         DO_3D( 0, 0, 0, 0, 1, nksr )
+         DO_3D( 0, 0, 0, 0, 1, nksr )             ! solar heat absorbed at T-point in the top 400m 
             zc0 = zz0 * EXP( -gdepw(ji,jj,jk  ,Kmm)*xsi0r ) + zz1 * EXP( -gdepw(ji,jj,jk  ,Kmm)*xsi1r )
             zc1 = zz0 * EXP( -gdepw(ji,jj,jk+1,Kmm)*xsi0r ) + zz1 * EXP( -gdepw(ji,jj,jk+1,Kmm)*xsi1r )
@@ -264,4 +264,6 @@
       END SELECT
       !
+      !                          !-----------------------------!
+      !                          !  update to the temp. trend  !
       !                          !-----------------------------!
       DO_3D( 0, 0, 0, 0, 1, nksr )
Index: /NEMO/trunk/src/OCE/TRA/trasbc.F90
===================================================================
--- /NEMO/trunk/src/OCE/TRA/trasbc.F90	(revision 13496)
+++ /NEMO/trunk/src/OCE/TRA/trasbc.F90	(revision 13497)
@@ -129,8 +129,8 @@
       END_2D
       IF( ln_linssh ) THEN                !* linear free surface  
-         DO_2D( 0, 1, 0, 0 )
+         DO_2D( 0, 1, 0, 0 )                    !==>> add concentration/dilution effect due to constant volume cell
             sbc_tsc(ji,jj,jp_tem) = sbc_tsc(ji,jj,jp_tem) + r1_rho0 * emp(ji,jj) * pts(ji,jj,1,jp_tem,Kmm)
             sbc_tsc(ji,jj,jp_sal) = sbc_tsc(ji,jj,jp_sal) + r1_rho0 * emp(ji,jj) * pts(ji,jj,1,jp_sal,Kmm)
-         END_2D
+         END_2D                                 !==>> output c./d. term
          IF( iom_use('emp_x_sst') )   CALL iom_put( "emp_x_sst", emp (:,:) * pts(:,:,1,jp_tem,Kmm) )
          IF( iom_use('emp_x_sss') )   CALL iom_put( "emp_x_sss", emp (:,:) * pts(:,:,1,jp_sal,Kmm) )
Index: /NEMO/trunk/src/OCE/TRA/trazdf.F90
===================================================================
--- /NEMO/trunk/src/OCE/TRA/trazdf.F90	(revision 13496)
+++ /NEMO/trunk/src/OCE/TRA/trazdf.F90	(revision 13497)
@@ -208,5 +208,5 @@
             !   used as a work space array: its value is modified.
             !
-            DO_2D( 0, 0, 0, 0 )
+            DO_2D( 0, 0, 0, 0 )      !* 1st recurrence:   Tk = Dk - Ik Sk-1 / Tk-1   (increasing k) ! done one for all passive tracers (so included in the IF instruction)
                zwt(ji,jj,1) = zwd(ji,jj,1)
             END_2D
@@ -217,5 +217,5 @@
          ENDIF 
          !         
-         DO_2D( 0, 0, 0, 0 )
+         DO_2D( 0, 0, 0, 0 )         !* 2nd recurrence:    Zk = Yk - Ik / Tk-1  Zk-1
             pt(ji,jj,1,jn,Kaa) =        e3t(ji,jj,1,Kbb) * pt(ji,jj,1,jn,Kbb)    &
                &               + p2dt * e3t(ji,jj,1,Kmm) * pt(ji,jj,1,jn,Krhs)
@@ -227,5 +227,5 @@
          END_3D
          !
-         DO_2D( 0, 0, 0, 0 )
+         DO_2D( 0, 0, 0, 0 )         !* 3d recurrence:    Xk = (Zk - Sk Xk+1 ) / Tk   (result is the after tracer)
             pt(ji,jj,jpkm1,jn,Kaa) = pt(ji,jj,jpkm1,jn,Kaa) / zwt(ji,jj,jpkm1) * tmask(ji,jj,jpkm1)
          END_2D
Index: /NEMO/trunk/src/OCE/TRA/zpshde.F90
===================================================================
--- /NEMO/trunk/src/OCE/TRA/zpshde.F90	(revision 13496)
+++ /NEMO/trunk/src/OCE/TRA/zpshde.F90	(revision 13497)
@@ -167,5 +167,5 @@
          CALL eos( ztj, zhj, zrj )        ! at the partial step depth output in  zri, zrj 
          !
-         DO_2D( 1, 0, 1, 0 )
+         DO_2D( 1, 0, 1, 0 )              ! Gradient of density at the last level
             iku = mbku(ji,jj)
             ikv = mbkv(ji,jj)
@@ -329,5 +329,5 @@
          CALL eos( ztj, zhj, zrj )
 
-         DO_2D( 1, 0, 1, 0 )
+         DO_2D( 1, 0, 1, 0 )            ! Gradient of density at the last level
             iku = mbku(ji,jj)
             ikv = mbkv(ji,jj)
@@ -420,5 +420,5 @@
          CALL eos( ztj, zhj, zrj )        ! at the partial step depth output in  zri, zrj 
          !
-         DO_2D( 1, 0, 1, 0 )
+         DO_2D( 1, 0, 1, 0 )              ! Gradient of density at the last level
             iku = miku(ji,jj) 
             ikv = mikv(ji,jj) 
Index: /NEMO/trunk/src/OCE/TRD/trddyn.F90
===================================================================
--- /NEMO/trunk/src/OCE/TRD/trddyn.F90	(revision 13496)
+++ /NEMO/trunk/src/OCE/TRD/trddyn.F90	(revision 13497)
@@ -124,5 +124,5 @@
                               z3dx(:,:,:) = 0._wp                  ! U.dxU & V.dyV (approximation)
                               z3dy(:,:,:) = 0._wp
-                              DO_3D( 0, 0, 0, 0, 1, jpkm1 )
+                              DO_3D( 0, 0, 0, 0, 1, jpkm1 )   ! no mask as un,vn are masked
                                  z3dx(ji,jj,jk) = uu(ji,jj,jk,Kmm) * ( uu(ji+1,jj,jk,Kmm) - uu(ji-1,jj,jk,Kmm) ) / ( 2._wp * e1u(ji,jj) )
                                  z3dy(ji,jj,jk) = vv(ji,jj,jk,Kmm) * ( vv(ji,jj+1,jk,Kmm) - vv(ji,jj-1,jk,Kmm) ) / ( 2._wp * e2v(ji,jj) )
Index: /NEMO/trunk/src/OCE/TRD/trdglo.F90
===================================================================
--- /NEMO/trunk/src/OCE/TRD/trdglo.F90	(revision 13496)
+++ /NEMO/trunk/src/OCE/TRD/trdglo.F90	(revision 13497)
@@ -86,5 +86,5 @@
       	!
       	CASE( 'TRA' )          !==  Tracers (T & S)  ==!
-      	   DO_3D( 1, 1, 1, 1, 1, jpkm1 )
+      	   DO_3D( 1, 1, 1, 1, 1, jpkm1 )   ! global sum of mask volume trend and trend*T (including interior mask)
 	            zvm = e1e2t(ji,jj) * e3t(ji,jj,jk,Kmm) * tmask(ji,jj,jk) * tmask_i(ji,jj)
 	            zvt = ptrdx(ji,jj,jk) * zvm
@@ -218,5 +218,5 @@
          END_3D
          
-         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
+         DO_3D( 0, 0, 0, 0, 1, jpkm1 )    ! Density flux divergence at t-point
             zkepe(ji,jj,jk) = - (  zkz(ji,jj,jk) - zkz(ji  ,jj  ,jk+1)               &
                &                 + zkx(ji,jj,jk) - zkx(ji-1,jj  ,jk  )               &
Index: /NEMO/trunk/src/OCE/TRD/trdmxl.F90
===================================================================
--- /NEMO/trunk/src/OCE/TRD/trdmxl.F90	(revision 13496)
+++ /NEMO/trunk/src/OCE/TRD/trdmxl.F90	(revision 13497)
@@ -120,5 +120,5 @@
          !
          wkx(:,:,:) = 0._wp         !==  now ML weights for vertical averaging  ==!
-         DO_3D( 1, 1, 1, 1, 1, jpktrd )
+         DO_3D( 1, 1, 1, 1, 1, jpktrd )  ! initialize wkx with vertical scale factor in mixed-layer
             IF( jk - kmxln(ji,jj) < 0 )   THEN
                wkx(ji,jj,jk) = e3t(ji,jj,jk,Kmm) * tmask(ji,jj,jk)
Index: /NEMO/trunk/src/OCE/TRD/trdtra.F90
===================================================================
--- /NEMO/trunk/src/OCE/TRD/trdtra.F90	(revision 13496)
+++ /NEMO/trunk/src/OCE/TRD/trdtra.F90	(revision 13497)
@@ -210,5 +210,5 @@
       !!----------------------------------------------------------------------
       !
-      SELECT CASE( cdir )      ! shift depending on the direction
+      SELECT CASE( cdir )             ! shift depending on the direction
       CASE( 'X' )   ;   ii = 1   ;   ij = 0   ;   ik = 0      ! i-trend
       CASE( 'Y' )   ;   ii = 0   ;   ij = 1   ;   ik = 0      ! j-trend
@@ -216,10 +216,10 @@
       END SELECT
       !
-      !                        ! set to zero uncomputed values
+      !                               ! set to zero uncomputed values
       ptrd(jpi,:,:) = 0._wp   ;   ptrd(1,:,:) = 0._wp
       ptrd(:,jpj,:) = 0._wp   ;   ptrd(:,1,:) = 0._wp
       ptrd(:,:,jpk) = 0._wp
       !
-      DO_3D( 0, 0, 0, 0, 1, jpkm1 )
+      DO_3D( 0, 0, 0, 0, 1, jpkm1 )   ! advective trend
          ptrd(ji,jj,jk) = - (     pf (ji,jj,jk) - pf (ji-ii,jj-ij,jk-ik)                        &
            &                  - ( pu(ji,jj,jk) - pu(ji-ii,jj-ij,jk-ik) ) * pt(ji,jj,jk)  )   &
Index: /NEMO/trunk/src/OCE/TRD/trdvor.F90
===================================================================
--- /NEMO/trunk/src/OCE/TRD/trdvor.F90	(revision 13496)
+++ /NEMO/trunk/src/OCE/TRD/trdvor.F90	(revision 13497)
@@ -103,7 +103,7 @@
       CASE( jpdyn_zad )   ;   CALL trd_vor_zint( putrd, pvtrd, jpvor_zad, Kmm )   ! Vertical Advection 
       CASE( jpdyn_spg )   ;   CALL trd_vor_zint( putrd, pvtrd, jpvor_spg, Kmm )   ! Surface Pressure Grad. 
-      CASE( jpdyn_zdf )                                                      ! Vertical Diffusion 
+      CASE( jpdyn_zdf )                                                           ! Vertical Diffusion 
          ztswu(:,:) = 0.e0   ;   ztswv(:,:) = 0.e0
-         DO_2D( 0, 0, 0, 0 )
+         DO_2D( 0, 0, 0, 0 )                                                               ! wind stress trends
             ztswu(ji,jj) = 0.5 * ( utau_b(ji,jj) + utau(ji,jj) ) / ( e3u(ji,jj,1,Kmm) * rho0 )
             ztswv(ji,jj) = 0.5 * ( vtau_b(ji,jj) + vtau(ji,jj) ) / ( e3v(ji,jj,1,Kmm) * rho0 )
Index: /NEMO/trunk/src/OCE/USR/usrdef_istate.F90
===================================================================
--- /NEMO/trunk/src/OCE/USR/usrdef_istate.F90	(revision 13496)
+++ /NEMO/trunk/src/OCE/USR/usrdef_istate.F90	(revision 13497)
@@ -57,9 +57,9 @@
       IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~   Ocean at rest, with an horizontally uniform T and S profiles'
       !
-      pu  (:,:,:) = 0._wp        ! ocean at rest
+      pu  (:,:,:) = 0._wp           ! ocean at rest
       pv  (:,:,:) = 0._wp
       pssh(:,:)   = 0._wp
       !
-      DO_3D( 1, 1, 1, 1, 1, jpk )
+      DO_3D( 1, 1, 1, 1, 1, jpk )   ! horizontally uniform T & S profiles
          pts(ji,jj,jk,jp_tem) =  (  (  16. - 12. * TANH( (pdept(ji,jj,jk) - 400) / 700 ) )   &
               &           * (-TANH( (500. - pdept(ji,jj,jk)) / 150. ) + 1.) / 2.             &
Index: /NEMO/trunk/src/OCE/ZDF/zdfddm.F90
===================================================================
--- /NEMO/trunk/src/OCE/ZDF/zdfddm.F90	(revision 13496)
+++ /NEMO/trunk/src/OCE/ZDF/zdfddm.F90	(revision 13497)
@@ -94,5 +94,5 @@
 !!gm                            and many acces in memory
          
-         DO_2D( 1, 1, 1, 1 )
+         DO_2D( 1, 1, 1, 1 )           !==  R=zrau = (alpha / beta) (dk[t] / dk[s])  ==!
             zrw =   ( gdepw(ji,jj,jk  ,Kmm) - gdept(ji,jj,jk,Kmm) )   &
 !!gm please, use e3w at Kmm below 
@@ -110,5 +110,5 @@
          END_2D
 
-         DO_2D( 1, 1, 1, 1 )
+         DO_2D( 1, 1, 1, 1 )           !==  indicators  ==!
             ! stability indicator: msks=1 if rn2>0; 0 elsewhere
             IF( rn2(ji,jj,jk) + 1.e-12  <= 0. ) THEN   ;   zmsks(ji,jj) = 0._wp
Index: /NEMO/trunk/src/OCE/ZDF/zdfdrg.F90
===================================================================
--- /NEMO/trunk/src/OCE/ZDF/zdfdrg.F90	(revision 13496)
+++ /NEMO/trunk/src/OCE/ZDF/zdfdrg.F90	(revision 13497)
@@ -431,5 +431,5 @@
             l_log_not_linssh = .FALSE.    !- don't update Cd at each time step
             !
-            DO_2D( 1, 1, 1, 1 )
+            DO_2D( 1, 1, 1, 1 )              ! pCd0 = mask (and boosted) logarithmic drag coef.
                zzz =  0.5_wp * e3t_0(ji,jj,k_mk(ji,jj))
                zcd = (  vkarmn / LOG( zzz / rn_z0 )  )**2
Index: /NEMO/trunk/src/OCE/ZDF/zdfgls.F90
===================================================================
--- /NEMO/trunk/src/OCE/ZDF/zdfgls.F90	(revision 13496)
+++ /NEMO/trunk/src/OCE/ZDF/zdfgls.F90	(revision 13497)
@@ -179,5 +179,5 @@
       
       ! Compute surface, top and bottom friction at T-points
-      DO_2D( 0, 0, 0, 0 )
+      DO_2D( 0, 0, 0, 0 )          !==  surface ocean friction
          ustar2_surf(ji,jj) = r1_rho0 * taum(ji,jj) * tmask(ji,jj,1)   ! surface friction
       END_2D
@@ -185,6 +185,6 @@
       !!gm Rq we may add here r_ke0(_top/_bot) ?  ==>> think about that...
       !    
-      IF( .NOT.ln_drg_OFF ) THEN   !== top/bottom friction   (explicit before friction)
-         DO_2D( 0, 0, 0, 0 )       ! bottom friction (explicit before friction)
+      IF( .NOT.ln_drg_OFF ) THEN     !== top/bottom friction   (explicit before friction)
+         DO_2D( 0, 0, 0, 0 )         ! bottom friction (explicit before friction)
             zmsku = ( 2._wp - umask(ji-1,jj,mbkt(ji,jj)) * umask(ji,jj,mbkt(ji,jj)) )
             zmskv = ( 2._wp - vmask(ji,jj-1,mbkt(ji,jj)) * vmask(ji,jj,mbkt(ji,jj)) )     ! (CAUTION: CdU<0)
@@ -193,5 +193,5 @@
          END_2D
          IF( ln_isfcav ) THEN
-            DO_2D( 0, 0, 0, 0 )    ! top friction
+            DO_2D( 0, 0, 0, 0 )      ! top friction
                zmsku = ( 2. - umask(ji-1,jj,mikt(ji,jj)) * umask(ji,jj,mikt(ji,jj)) )
                zmskv = ( 2. - vmask(ji,jj-1,mikt(ji,jj)) * vmask(ji,jj,mikt(ji,jj)) )     ! (CAUTION: CdU<0)
@@ -220,5 +220,5 @@
       zhsro(:,:) = ( (1._wp-zice_fra(:,:)) * zhsro(:,:) + zice_fra(:,:) * rn_hsri )*tmask(:,:,1)  + (1._wp - tmask(:,:,1))*rn_hsro
       !
-      DO_3D( 1, 0, 1, 0, 2, jpkm1 )
+      DO_3D( 1, 0, 1, 0, 2, jpkm1 )  !==  Compute dissipation rate  ==!
          eps(ji,jj,jk)  = rc03 * en(ji,jj,jk) * SQRT( en(ji,jj,jk) ) / hmxl_n(ji,jj,jk)
       END_3D
@@ -416,11 +416,11 @@
       ! ----------------------------------------------------------
       !
-      DO_3D( 0, 0, 0, 0, 2, jpkm1 )
+      DO_3D( 0, 0, 0, 0, 2, jpkm1 )                ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1
          zdiag(ji,jj,jk) = zdiag(ji,jj,jk) - zd_lw(ji,jj,jk) * zd_up(ji,jj,jk-1) / zdiag(ji,jj,jk-1)
       END_3D
-      DO_3D( 0, 0, 0, 0, 2, jpk )
+      DO_3D( 0, 0, 0, 0, 2, jpk )                  ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1
          zd_lw(ji,jj,jk) = en(ji,jj,jk) - zd_lw(ji,jj,jk) / zdiag(ji,jj,jk-1) * zd_lw(ji,jj,jk-1)
       END_3D
-      DO_3DS( 0, 0, 0, 0, jpk-1, 2, -1 )
+      DO_3DS( 0, 0, 0, 0, jpk-1, 2, -1 )           ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk
          en(ji,jj,jk) = ( zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) * en(ji,jj,jk+1) ) / zdiag(ji,jj,jk)
       END_3D
@@ -610,11 +610,11 @@
       ! ----------------
       !
-      DO_3D( 0, 0, 0, 0, 2, jpkm1 )
+      DO_3D( 0, 0, 0, 0, 2, jpkm1 )                ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1
          zdiag(ji,jj,jk) = zdiag(ji,jj,jk) - zd_lw(ji,jj,jk) * zd_up(ji,jj,jk-1) / zdiag(ji,jj,jk-1)
       END_3D
-      DO_3D( 0, 0, 0, 0, 2, jpk )
+      DO_3D( 0, 0, 0, 0, 2, jpk )                  ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1
          zd_lw(ji,jj,jk) = psi(ji,jj,jk) - zd_lw(ji,jj,jk) / zdiag(ji,jj,jk-1) * zd_lw(ji,jj,jk-1)
       END_3D
-      DO_3DS( 0, 0, 0, 0, jpk-1, 2, -1 )
+      DO_3DS( 0, 0, 0, 0, jpk-1, 2, -1 )           ! Third recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk
          psi(ji,jj,jk) = ( zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) * psi(ji,jj,jk+1) ) / zdiag(ji,jj,jk)
       END_3D
@@ -652,5 +652,5 @@
       ! Limit dissipation rate under stable stratification
       ! --------------------------------------------------
-      DO_3D( 0, 0, 0, 0, 1, jpkm1 )
+      DO_3D( 0, 0, 0, 0, 1, jpkm1 )   ! Note that this set boundary conditions on hmxl_n at the same time
          ! limitation
          eps   (ji,jj,jk)  = MAX( eps(ji,jj,jk), rn_epsmin )
@@ -717,5 +717,5 @@
       ! default value, in case jpk > mbkt(ji,jj)+1. Not needed but avoid a bug when looking for undefined values (-fpe0)
       zstm(:,:,jpk) = 0.  
-      DO_2D( 0, 0, 0, 0 )
+      DO_2D( 0, 0, 0, 0 )             ! update bottom with good values
          zstm(ji,jj,mbkt(ji,jj)+1) = zstm(ji,jj,mbkt(ji,jj))
       END_2D
Index: /NEMO/trunk/src/OCE/ZDF/zdfiwm.F90
===================================================================
--- /NEMO/trunk/src/OCE/ZDF/zdfiwm.F90	(revision 13496)
+++ /NEMO/trunk/src/OCE/ZDF/zdfiwm.F90	(revision 13497)
@@ -164,5 +164,5 @@
       !                       !* Critical slope mixing: distribute energy over the time-varying ocean depth,
       !                                                 using an exponential decay from the seafloor.
-      DO_2D( 0, 0, 0, 0 )
+      DO_2D( 0, 0, 0, 0 )             ! part independent of the level
          zhdep(ji,jj) = gdepw_0(ji,jj,mbkt(ji,jj)+1)       ! depth of the ocean
          zfact(ji,jj) = rho0 * (  1._wp - EXP( -zhdep(ji,jj) / hcri_iwm(ji,jj) )  )
@@ -170,5 +170,5 @@
       END_2D
 !!gm gde3w ==>>>  check for ssh taken into account.... seem OK gde3w_n=gdept(:,:,:,Kmm) - ssh(:,:,Kmm)
-      DO_3D( 0, 0, 0, 0, 2, jpkm1 )
+      DO_3D( 0, 0, 0, 0, 2, jpkm1 )   ! complete with the level-dependent part
          IF ( zfact(ji,jj) == 0._wp .OR. wmask(ji,jj,jk) == 0._wp ) THEN   ! optimization
             zemx_iwm(ji,jj,jk) = 0._wp
@@ -293,6 +293,6 @@
       END_3D
       !
-      IF( ln_mevar ) THEN              ! Variable mixing efficiency case : modify zav_wave in the
-         DO_3D( 0, 0, 0, 0, 2, jpkm1 )
+      IF( ln_mevar ) THEN                ! Variable mixing efficiency case : modify zav_wave in the
+         DO_3D( 0, 0, 0, 0, 2, jpkm1 )   ! energetic (Reb > 480) and buoyancy-controlled (Reb <10.224 ) regimes
             IF( zReb(ji,jj,jk) > 480.00_wp ) THEN
                zav_wave(ji,jj,jk) = 3.6515_wp * znu_w(ji,jj,jk) * SQRT( zReb(ji,jj,jk) )
@@ -303,5 +303,5 @@
       ENDIF
       !
-      DO_3D( 0, 0, 0, 0, 2, jpkm1 )          ! Bound diffusivity by molecular value and 100 cm2/s
+      DO_3D( 0, 0, 0, 0, 2, jpkm1 )      ! Bound diffusivity by molecular value and 100 cm2/s
          zav_wave(ji,jj,jk) = MIN(  MAX( 1.4e-7_wp, zav_wave(ji,jj,jk) ), 1.e-2_wp  ) * wmask(ji,jj,jk)
       END_3D
@@ -330,7 +330,7 @@
       !                          ! ----------------------- !
       !      
-      IF( ln_tsdiff ) THEN          !* Option for differential mixing of salinity and temperature
+      IF( ln_tsdiff ) THEN                !* Option for differential mixing of salinity and temperature
          ztmp1 = 0.505_wp + 0.495_wp * TANH( 0.92_wp * ( LOG10( 1.e-20_wp ) - 0.60_wp ) )
-         DO_3D( 0, 0, 0, 0, 2, jpkm1 )
+         DO_3D( 0, 0, 0, 0, 2, jpkm1 )       ! Calculate S/T diffusivity ratio as a function of Reb
             ztmp2 = zReb(ji,jj,jk) * 5._wp * r1_6
             IF ( ztmp2 > 1.e-20_wp .AND. wmask(ji,jj,jk) == 1._wp ) THEN
@@ -347,5 +347,5 @@
          END_3D
          !
-      ELSE                          !* update momentum & tracer diffusivity with wave-driven mixing
+      ELSE                                !* update momentum & tracer diffusivity with wave-driven mixing
          DO_3D( 0, 0, 0, 0, 2, jpkm1 )
             p_avs(ji,jj,jk) = p_avs(ji,jj,jk) + zav_wave(ji,jj,jk)
@@ -355,9 +355,9 @@
       ENDIF
 
-      !                             !* output internal wave-driven mixing coefficient
+      !                                   !* output internal wave-driven mixing coefficient
       CALL iom_put( "av_wave", zav_wave )
-                                    !* output useful diagnostics: Kz*N^2 , 
+                                          !* output useful diagnostics: Kz*N^2 , 
 !!gm Kz*N2 should take into account the ratio avs/avt if it is used.... (see diaar5)
-                                    !  vertical integral of rho0 * Kz * N^2 , energy density (zemx_iwm)
+                                          !  vertical integral of rho0 * Kz * N^2 , energy density (zemx_iwm)
       IF( iom_use("bflx_iwm") .OR. iom_use("pcmap_iwm") ) THEN
          ALLOCATE( z2d(jpi,jpj) , z3d(jpi,jpj,jpk) )
Index: /NEMO/trunk/src/OCE/ZDF/zdfmxl.F90
===================================================================
--- /NEMO/trunk/src/OCE/ZDF/zdfmxl.F90	(revision 13496)
+++ /NEMO/trunk/src/OCE/ZDF/zdfmxl.F90	(revision 13497)
@@ -96,8 +96,8 @@
       !
       ! w-level of the mixing and mixed layers
-      nmln(:,:)  = nlb10               ! Initialization to the number of w ocean point
-      hmlp(:,:)  = 0._wp               ! here hmlp used as a dummy variable, integrating vertically N^2
-      zN2_c = grav * rho_c * r1_rho0   ! convert density criteria into N^2 criteria
-      DO_3D( 1, 1, 1, 1, nlb10, jpkm1 )
+      nmln(:,:)  = nlb10                  ! Initialization to the number of w ocean point
+      hmlp(:,:)  = 0._wp                  ! here hmlp used as a dummy variable, integrating vertically N^2
+      zN2_c = grav * rho_c * r1_rho0      ! convert density criteria into N^2 criteria
+      DO_3D( 1, 1, 1, 1, nlb10, jpkm1 )   ! Mixed layer level: w-level
          ikt = mbkt(ji,jj)
          hmlp(ji,jj) =   &
@@ -107,6 +107,6 @@
       !
       ! w-level of the turbocline and mixing layer (iom_use)
-      imld(:,:) = mbkt(:,:) + 1        ! Initialization to the number of w ocean point
-      DO_3DS( 1, 1, 1, 1, jpkm1, nlb10, -1 )
+      imld(:,:) = mbkt(:,:) + 1                ! Initialization to the number of w ocean point
+      DO_3DS( 1, 1, 1, 1, jpkm1, nlb10, -1 )   ! from the bottom to nlb10
          IF( avt (ji,jj,jk) < avt_c * wmask(ji,jj,jk) )   imld(ji,jj) = jk      ! Turbocline 
       END_3D
Index: /NEMO/trunk/src/OCE/ZDF/zdfosm.F90
===================================================================
--- /NEMO/trunk/src/OCE/ZDF/zdfosm.F90	(revision 13496)
+++ /NEMO/trunk/src/OCE/ZDF/zdfosm.F90	(revision 13497)
@@ -1184,5 +1184,5 @@
 ! KPP-style Ri# mixing
        IF( ln_kpprimix) THEN
-          DO_3D( 1, 0, 1, 0, 2, jpkm1 )
+          DO_3D( 1, 0, 1, 0, 2, jpkm1 )      !* Shear production at uw- and vw-points (energy conserving form)
              z3du(ji,jj,jk) = 0.5 * (  uu(ji,jj,jk-1,Kmm) -  uu(ji  ,jj,jk,Kmm) )   &
                   &                 * (  uu(ji,jj,jk-1,Kbb) -  uu(ji  ,jj,jk,Kbb) ) * wumask(ji,jj,jk) &
@@ -1516,5 +1516,5 @@
      !
      hbl(:,:)  = 0._wp              ! here hbl used as a dummy variable, integrating vertically N^2
-     DO_3D( 1, 1, 1, 1, 1, jpkm1 )
+     DO_3D( 1, 1, 1, 1, 1, jpkm1 )  ! Mixed layer level: w-level
         ikt = mbkt(ji,jj)
         hbl(ji,jj) = hbl(ji,jj) + MAX( rn2(ji,jj,jk) , 0._wp ) * e3w(ji,jj,jk,Kmm)
@@ -1629,5 +1629,5 @@
       !code saving tracer trends removed, replace with trdmxl_oce
 
-      DO_3D( 0, 0, 0, 0, 1, jpkm1 )
+      DO_3D( 0, 0, 0, 0, 1, jpkm1 )       ! add non-local u and v fluxes
          puu(ji,jj,jk,Krhs) =  puu(ji,jj,jk,Krhs)                      &
             &                 - (  ghamu(ji,jj,jk  )  &
Index: /NEMO/trunk/src/OCE/ZDF/zdfric.F90
===================================================================
--- /NEMO/trunk/src/OCE/ZDF/zdfric.F90	(revision 13496)
+++ /NEMO/trunk/src/OCE/ZDF/zdfric.F90	(revision 13497)
@@ -160,5 +160,5 @@
       !
       !                       !==  avm and avt = F(Richardson number)  ==!
-      DO_3D( 1, 0, 1, 0, 2, jpkm1 )
+      DO_3D( 1, 0, 1, 0, 2, jpkm1 )       ! coefficient = F(richardson number) (avm-weighted Ri)
          zcfRi = 1._wp / (  1._wp + rn_alp * MAX(  0._wp , avm(ji,jj,jk) * rn2(ji,jj,jk) / ( p_sh2(ji,jj,jk) + 1.e-20 ) )  )
          zav   = rn_avmri * zcfRi**nn_ric
@@ -173,10 +173,10 @@
       IF( ln_mldw ) THEN      !==  set a minimum value in the Ekman layer  ==!
          !
-         DO_2D( 0, 0, 0, 0 )
+         DO_2D( 0, 0, 0, 0 )             !* Ekman depth
             zustar = SQRT( taum(ji,jj) * r1_rho0 )
             zhek   = rn_ekmfc * zustar / ( ABS( ff_t(ji,jj) ) + rsmall )   ! Ekman depth
             zh_ekm(ji,jj) = MAX(  rn_mldmin , MIN( zhek , rn_mldmax )  )   ! set allowed range
          END_2D
-         DO_3D( 0, 0, 0, 0, 2, jpkm1 )
+         DO_3D( 0, 0, 0, 0, 2, jpkm1 )   !* minimum mixing coeff. within the Ekman layer
             IF( gdept(ji,jj,jk,Kmm) < zh_ekm(ji,jj) ) THEN
                p_avm(ji,jj,jk) = MAX(  p_avm(ji,jj,jk), rn_wvmix  ) * wmask(ji,jj,jk)
Index: /NEMO/trunk/src/OCE/ZDF/zdfsh2.F90
===================================================================
--- /NEMO/trunk/src/OCE/ZDF/zdfsh2.F90	(revision 13496)
+++ /NEMO/trunk/src/OCE/ZDF/zdfsh2.F90	(revision 13497)
@@ -60,5 +60,5 @@
       !
       DO jk = 2, jpkm1
-         DO_2D( 1, 0, 1, 0 )
+         DO_2D( 1, 0, 1, 0 )     !* 2 x shear production at uw- and vw-points (energy conserving form)
             zsh2u(ji,jj) = ( p_avm(ji+1,jj,jk) + p_avm(ji,jj,jk) ) &
                &         * (   uu(ji,jj,jk-1,Kmm) -   uu(ji,jj,jk,Kmm) ) &
@@ -72,5 +72,5 @@
                &         * wvmask(ji,jj,jk)
          END_2D
-         DO_2D( 0, 0, 0, 0 )
+         DO_2D( 0, 0, 0, 0 )     !* shear production at w-point ! coast mask: =2 at the coast ; =1 otherwise (NB: wmask useless as zsh2 are masked)
             p_sh2(ji,jj,jk) = 0.25 * (   ( zsh2u(ji-1,jj) + zsh2u(ji,jj) ) * ( 2. - umask(ji-1,jj,jk) * umask(ji,jj,jk) )   &
                &                       + ( zsh2v(ji,jj-1) + zsh2v(ji,jj) ) * ( 2. - vmask(ji,jj-1,jk) * vmask(ji,jj,jk) )   )
Index: /NEMO/trunk/src/OCE/ZDF/zdftke.F90
===================================================================
--- /NEMO/trunk/src/OCE/ZDF/zdftke.F90	(revision 13496)
+++ /NEMO/trunk/src/OCE/ZDF/zdftke.F90	(revision 13497)
@@ -238,5 +238,5 @@
       !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
       !
-      DO_2D( 0, 0, 0, 0 )
+      DO_2D( 0, 0, 0, 0 )         ! en(1)   = rn_ebb taum / rau0  (min value rn_emin0)
 !! clem: this should be the right formulation but it makes the model unstable unless drags are calculated implicitly
 !!       one way around would be to increase zbbirau 
@@ -325,5 +325,5 @@
       !                     ! zdiag : diagonal zd_up : upper diagonal zd_lw : lower diagonal
       !
-      IF( nn_pdl == 1 ) THEN      !* Prandtl number = F( Ri )
+      IF( nn_pdl == 1 ) THEN          !* Prandtl number = F( Ri )
          DO_3D( 0, 0, 0, 0, 2, jpkm1 )
             !                             ! local Richardson number
@@ -338,5 +338,5 @@
       ENDIF
       !         
-      DO_3D( 0, 0, 0, 0, 2, jpkm1 )
+      DO_3D( 0, 0, 0, 0, 2, jpkm1 )   !* Matrix and right hand side in en
          zcof   = zfact1 * tmask(ji,jj,jk)
          !                                   ! A minimum of 2.e-5 m2/s is imposed on TKE vertical
@@ -358,8 +358,8 @@
       END_3D
       !                          !* Matrix inversion from level 2 (tke prescribed at level 1)
-      DO_3D( 0, 0, 0, 0, 3, jpkm1 )
+      DO_3D( 0, 0, 0, 0, 3, jpkm1 )                ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1
          zdiag(ji,jj,jk) = zdiag(ji,jj,jk) - zd_lw(ji,jj,jk) * zd_up(ji,jj,jk-1) / zdiag(ji,jj,jk-1)
       END_3D
-      DO_2D( 0, 0, 0, 0 )
+      DO_2D( 0, 0, 0, 0 )                          ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1
          zd_lw(ji,jj,2) = en(ji,jj,2) - zd_lw(ji,jj,2) * en(ji,jj,1)    ! Surface boudary conditions on tke
       END_2D
@@ -367,5 +367,5 @@
          zd_lw(ji,jj,jk) = en(ji,jj,jk) - zd_lw(ji,jj,jk) / zdiag(ji,jj,jk-1) *zd_lw(ji,jj,jk-1)
       END_3D
-      DO_2D( 0, 0, 0, 0 )
+      DO_2D( 0, 0, 0, 0 )                          ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk
          en(ji,jj,jpkm1) = zd_lw(ji,jj,jpkm1) / zdiag(ji,jj,jpkm1)
       END_2D
@@ -373,5 +373,5 @@
          en(ji,jj,jk) = ( zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) * en(ji,jj,jk+1) ) / zdiag(ji,jj,jk)
       END_3D
-      DO_3D( 0, 0, 0, 0, 2, jpkm1 )
+      DO_3D( 0, 0, 0, 0, 2, jpkm1 )                ! set the minimum value of tke
          en(ji,jj,jk) = MAX( en(ji,jj,jk), rn_emin ) * wmask(ji,jj,jk)
       END_3D
@@ -396,5 +396,5 @@
          END_2D
       ELSEIF( nn_etau == 3 ) THEN       !* penetration belox the mixed layer (HF variability)
-         DO_3D( 0, 0, 0, 0, 2, jpkm1 )
+         DO_3D( 0, 0, 0, 0, 2, jpkm1 )         ! nn_eice=0 : ON below sea-ice ; nn_eice>0 : partly OFF
             ztx2 = utau(ji-1,jj  ) + utau(ji,jj)
             zty2 = vtau(ji  ,jj-1) + vtau(ji,jj)
@@ -470,5 +470,5 @@
          zraug = vkarmn * 2.e5_wp / ( rho0 * grav )
 #if ! defined key_si3 && ! defined key_cice
-         DO_2D( 0, 0, 0, 0 )
+         DO_2D( 0, 0, 0, 0 )                  ! No sea-ice
             zmxlm(ji,jj,1) =  zraug * taum(ji,jj) * tmask(ji,jj,1)
          END_2D
@@ -481,5 +481,5 @@
             END_2D
             !
-         CASE( 1 )                           ! scaling with constant sea-ice thickness
+         CASE( 1 )                      ! scaling with constant sea-ice thickness
             DO_2D( 0, 0, 0, 0 )
                zmxlm(ji,jj,1) =  ( ( 1._wp - fr_i(ji,jj) ) * zraug * taum(ji,jj) + &
@@ -487,5 +487,5 @@
             END_2D
             !
-         CASE( 2 )                                 ! scaling with mean sea-ice thickness
+         CASE( 2 )                      ! scaling with mean sea-ice thickness
             DO_2D( 0, 0, 0, 0 )
 #if defined key_si3
@@ -499,5 +499,5 @@
             END_2D
             !
-         CASE( 3 )                                 ! scaling with max sea-ice thickness
+         CASE( 3 )                      ! scaling with max sea-ice thickness
             DO_2D( 0, 0, 0, 0 )
                zmaxice = MAXVAL( h_i(ji,jj,:) )
@@ -551,9 +551,9 @@
          !
       CASE ( 2 )           ! |dk[xml]| bounded by e3t :
-         DO_3D( 0, 0, 0, 0, 2, jpkm1 )
+         DO_3D( 0, 0, 0, 0, 2, jpkm1 )        ! from the surface to the bottom :
             zmxlm(ji,jj,jk) =   &
                &    MIN( zmxlm(ji,jj,jk-1) + e3t(ji,jj,jk-1,Kmm), zmxlm(ji,jj,jk) )
          END_3D
-         DO_3DS( 0, 0, 0, 0, jpkm1, 2, -1 )
+         DO_3DS( 0, 0, 0, 0, jpkm1, 2, -1 )   ! from the bottom to the surface :
             zemxl = MIN( zmxlm(ji,jj,jk+1) + e3t(ji,jj,jk+1,Kmm), zmxlm(ji,jj,jk) )
             zmxlm(ji,jj,jk) = zemxl
@@ -562,9 +562,9 @@
          !
       CASE ( 3 )           ! lup and ldown, |dk[xml]| bounded by e3t :
-         DO_3D( 0, 0, 0, 0, 2, jpkm1 )
+         DO_3D( 0, 0, 0, 0, 2, jpkm1 )        ! from the surface to the bottom : lup
             zmxld(ji,jj,jk) =    &
                &    MIN( zmxld(ji,jj,jk-1) + e3t(ji,jj,jk-1,Kmm), zmxlm(ji,jj,jk) )
          END_3D
-         DO_3DS( 0, 0, 0, 0, jpkm1, 2, -1 )
+         DO_3DS( 0, 0, 0, 0, jpkm1, 2, -1 )   ! from the bottom to the surface : ldown
             zmxlm(ji,jj,jk) =   &
                &    MIN( zmxlm(ji,jj,jk+1) + e3t(ji,jj,jk+1,Kmm), zmxlm(ji,jj,jk) )
@@ -582,5 +582,5 @@
       !                     !  Vertical eddy viscosity and diffusivity  (avm and avt)
       !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-      DO_3D( 0, 0, 0, 0, 1, jpkm1 )
+      DO_3D( 0, 0, 0, 0, 1, jpkm1 )   !* vertical eddy viscosity & diffivity at w-points
          zsqen = SQRT( en(ji,jj,jk) )
          zav   = rn_ediff * zmxlm(ji,jj,jk) * zsqen
@@ -591,5 +591,5 @@
       !
       !
-      IF( nn_pdl == 1 ) THEN      !* Prandtl number case: update avt
+      IF( nn_pdl == 1 ) THEN          !* Prandtl number case: update avt
          DO_3D( 0, 0, 0, 0, 2, jpkm1 )
             p_avt(ji,jj,jk)   = MAX( apdlr(ji,jj,jk) * p_avt(ji,jj,jk), avtb_2d(ji,jj) * avtb(jk) ) * wmask(ji,jj,jk)
Index: /NEMO/trunk/src/OFF/dtadyn.F90
===================================================================
--- /NEMO/trunk/src/OFF/dtadyn.F90	(revision 13496)
+++ /NEMO/trunk/src/OFF/dtadyn.F90	(revision 13497)
@@ -412,6 +412,5 @@
             ENDIF
          END_2D
-!!st pourquoi on n'utilise pas le gde3w ici plutôt que de faire une boucle ? 
-         DO_2D( 1, 1, 1, 1 )
+         DO_2D( 1, 1, 1, 1 )                           ! set the associated depth
             h_rnf(ji,jj) = 0._wp
             DO jk = 1, nk_rnf(ji,jj)
@@ -688,6 +687,5 @@
       !!----------------------------------------------------------------------
       !
-!!st code dupliqué même remarque que plus haut pourquoi ne pas utiliser gdepw ?
-      DO_2D( 1, 1, 1, 1 )
+      DO_2D( 1, 1, 1, 1 )               ! update the depth over which runoffs are distributed
          h_rnf(ji,jj) = 0._wp
          DO jk = 1, nk_rnf(ji,jj)                           ! recalculates h_rnf to be the depth in metres
Index: /NEMO/trunk/src/TOP/C14/trcatm_c14.F90
===================================================================
--- /NEMO/trunk/src/TOP/C14/trcatm_c14.F90	(revision 13496)
+++ /NEMO/trunk/src/TOP/C14/trcatm_c14.F90	(revision 13497)
@@ -120,5 +120,5 @@
             IF( ierr3 /= 0 )   CALL ctl_stop( 'STOP', 'trc_atm_c14_ini: unable to allocate fareaz' )
       !
-            DO_2D( 1, 1, 1, 1 )
+            DO_2D( 1, 1, 1, 1 )                 ! from C14b package
               IF( gphit(ji,jj) >= yn40 ) THEN
                  fareaz(ji,jj,1) = 0.
Index: /NEMO/trunk/src/TOP/CFC/trcsms_cfc.F90
===================================================================
--- /NEMO/trunk/src/TOP/CFC/trcsms_cfc.F90	(revision 13496)
+++ /NEMO/trunk/src/TOP/CFC/trcsms_cfc.F90	(revision 13497)
@@ -126,6 +126,6 @@
          
          !                                                         !------------!
-         DO_2D( 1, 1, 1, 1 )
- 
+         DO_2D( 1, 1, 1, 1 )                                       !  i-j loop  !
+            !                                                      !------------!
             ! space interpolation
             zpp_cfc  =       xphem(ji,jj)   * zpatm(1,jl)   &
Index: /NEMO/trunk/src/TOP/PISCES/P2Z/p2zopt.F90
===================================================================
--- /NEMO/trunk/src/TOP/PISCES/P2Z/p2zopt.F90	(revision 13496)
+++ /NEMO/trunk/src/TOP/PISCES/P2Z/p2zopt.F90	(revision 13497)
@@ -95,5 +95,5 @@
       !                                          ! Photosynthetically Available Radiation (PAR)
       zcoef = 12 * redf / rcchl / rpig           ! --------------------------------------
-      DO_3D( 1, 1, 1, 1, 2, jpk )
+      DO_3D( 1, 1, 1, 1, 2, jpk )                     ! local par at w-levels
          zpig = LOG(  MAX( TINY(0.), tr(ji,jj,jk-1,jpphy,Kmm) ) * zcoef  )
          zkr  = xkr0 + xkrp * EXP( xlr * zpig )
@@ -102,5 +102,5 @@
          zparg(ji,jj,jk) = zparg(ji,jj,jk-1) * EXP( -zkg * e3t(ji,jj,jk-1,Kmm) )
       END_3D
-      DO_3D( 1, 1, 1, 1, 1, jpkm1 )
+      DO_3D( 1, 1, 1, 1, 1, jpkm1 )                   ! mean par at t-levels
          zpig = LOG(  MAX( TINY(0.), tr(ji,jj,jk,jpphy,Kmm) ) * zcoef  )
          zkr  = xkr0 + xkrp * EXP( xlr * zpig )
@@ -114,5 +114,5 @@
       !                                          ! --------------
       neln(:,:) = 1                                   ! euphotic layer level
-      DO_3D( 1, 1, 1, 1, 1, jpkm1 )
+      DO_3D( 1, 1, 1, 1, 1, jpkm1 )                   ! (i.e. 1rst T-level strictly below EL bottom)
         IF( etot(ji,jj,jk) >= zpar100(ji,jj) )   neln(ji,jj) = jk + 1 
       END_3D
Index: /NEMO/trunk/src/TOP/TRP/trdmxl_trc.F90
===================================================================
--- /NEMO/trunk/src/TOP/TRP/trdmxl_trc.F90	(revision 13496)
+++ /NEMO/trunk/src/TOP/TRP/trdmxl_trc.F90	(revision 13497)
@@ -148,5 +148,5 @@
          ! ... Weights for vertical averaging
          wkx_trc(:,:,:) = 0.e0
-         DO_3D( 1, 1, 1, 1, 1, jpktrd_trc )
+         DO_3D( 1, 1, 1, 1, 1, jpktrd_trc )                       ! initialize wkx_trc with vertical scale factor in mixed-layer
             IF( jk - nmld_trc(ji,jj) < 0 )   wkx_trc(ji,jj,jk) = e3t(ji,jj,jk,Kmm) * tmask(ji,jj,jk)
          END_3D
Index: /NEMO/trunk/src/TOP/trcdta.F90
===================================================================
--- /NEMO/trunk/src/TOP/trcdta.F90	(revision 13496)
+++ /NEMO/trunk/src/TOP/trcdta.F90	(revision 13497)
@@ -199,5 +199,5 @@
                WRITE(numout,*) 'trc_dta: interpolates passive tracer data onto the s- or mixed s-z-coordinate mesh'
             ENDIF
-            DO_2D( 1, 1, 1, 1 )
+            DO_2D( 1, 1, 1, 1 )                 ! vertical interpolation of T & S
                DO jk = 1, jpk                        ! determines the intepolated T-S profiles at each (i,j) points
                   zl = gdept(ji,jj,jk,Kmm)
