Index: branches/dev_003_CPL/NEMO/LIM_SRC_3/limtrp.F90
===================================================================
--- branches/dev_003_CPL/NEMO/LIM_SRC_3/limtrp.F90	(revision 893)
+++ branches/dev_003_CPL/NEMO/LIM_SRC_3/limtrp.F90	(revision 990)
@@ -57,5 +57,5 @@
 CONTAINS
 
-   SUBROUTINE lim_trp
+   SUBROUTINE lim_trp( kt ) 
       !!-------------------------------------------------------------------
       !!                   ***  ROUTINE lim_trp ***
@@ -75,7 +75,8 @@
       !!   3.0  !  05-11 (M. Vancoppenolle)   Multi-layer sea ice, salinity variations
       !!---------------------------------------------------------------------
+      INTEGER, INTENT(in) ::   kt     ! number of iteration
       !! * Local Variables
       INTEGER  ::   ji, jj, jk, jl, layer, &  ! dummy loop indices
-                    initad           ! number of sub-timestep for the advection
+         initad           ! number of sub-timestep for the advection
       INTEGER  ::   ji_maxu, ji_maxv, jj_maxu, jj_maxv
 
@@ -102,5 +103,5 @@
          zs0sm , zs0oi
 
-! MHE Multilayer heat content
+      ! MHE Multilayer heat content
       REAL(wp), DIMENSION(jpi,jpj,jkmax,jpl)  ::   &  ! temporary workspace
          zs0e
@@ -111,18 +112,19 @@
 
       zsm(:,:) = area(:,:)
-      
+
       IF( ln_limdyn ) THEN
-         WRITE(numout,*)
-         WRITE(numout,*) ' lim_trp : Ice Advection'
-         WRITE(numout,*) ' ~~~~~~~'
-
-!-----------------------------------------------------------------------------!
-! 1) CFL Test                                                             
-!-----------------------------------------------------------------------------!
+         IF( kt == nit000 .AND. lwp ) THEN
+            WRITE(numout,*) ' lim_trp : Ice Advection'
+            WRITE(numout,*) ' ~~~~~~~'
+         ENDIF
+
+         !-----------------------------------------------------------------------------!
+         ! 1) CFL Test                                                             
+         !-----------------------------------------------------------------------------!
 
          !------------------------------------------
          ! ice velocities at ocean U- and V-points 
          !------------------------------------------
-         
+
          ! zvbord factor between 1 and 2 to take into account slip or no-slip boundary conditions.        
          zvbord = 1.0 + ( 1.0 - bound )
@@ -166,9 +168,9 @@
 
          IF ( zcfl > 0.5 .AND. lwp ) &
-         WRITE(numout,*) 'lim_trp : violation of cfl criterion the ',nday,'th day, cfl = ',zcfl
-
-!-----------------------------------------------------------------------------!
-! 2) Computation of transported fields                                        
-!-----------------------------------------------------------------------------!
+            WRITE(numout,*) 'lim_trp : violation of cfl criterion the ',nday,'th day, cfl = ',zcfl
+
+         !-----------------------------------------------------------------------------!
+         ! 2) Computation of transported fields                                        
+         !-----------------------------------------------------------------------------!
 
          !------------------------------------------------------
@@ -185,7 +187,7 @@
             zs0oi (:,:,jl) =  oa_i (:,:,jl)  * area(:,:)    ! Age content
 
-         !----------------------------------
-         ! 1.2) Ice and snow heat contents
-         !----------------------------------
+            !----------------------------------
+            ! 1.2) Ice and snow heat contents
+            !----------------------------------
 
             zs0c0 (:,:,jl)     = e_s(:,:,1,jl)              ! Snow heat cont.
@@ -195,73 +197,73 @@
          END DO
 
-!-----------------------------------------------------------------------------!
-! 3) Advection of Ice fields                                              
-!-----------------------------------------------------------------------------!
+         !-----------------------------------------------------------------------------!
+         ! 3) Advection of Ice fields                                              
+         !-----------------------------------------------------------------------------!
 
          ! If ice drift field is too fast, use an appropriate time step for advection.         
          initad = 1 + INT( MAX( rzero, SIGN( rone, zcfl-0.5 ) ) )
          zusnit = 1.0 / REAL( initad ) 
-         
+
          IF ( MOD( nday , 2 ) == 0) THEN
             DO jk = 1,initad
                !--- ice open water area
                CALL lim_adv_x( zusnit, zui_u, rone , zsm, zs0ow(:,:) , sxopw(:,:) , & 
-                                                          sxxopw(:,:), syopw(:,:) , & 
-                                                          syyopw(:,:), sxopw(:,:) )
+                  sxxopw(:,:), syopw(:,:) , & 
+                  syyopw(:,:), sxopw(:,:) )
                CALL lim_adv_y( zusnit, zvi_v, rzero, zsm, zs0ow(:,:) , sxopw (:,:) , &
-                                                          sxxopw(:,:), syopw (:,:) , & 
-                                                          syyopw(:,:), sxyopw(:,:) )
+                  sxxopw(:,:), syopw (:,:) , & 
+                  syyopw(:,:), sxyopw(:,:) )
                DO jl = 1, jpl
                   !--- ice volume  ---
                   CALL lim_adv_x( zusnit, zui_u, rone , zsm, zs0ice(:,:,jl) , sxice (:,:,jl) , & 
-                                                             sxxice(:,:,jl) , syice (:,:,jl) , & 
-                                                             syyice(:,:,jl) , sxyice(:,:,jl) )
+                     sxxice(:,:,jl) , syice (:,:,jl) , & 
+                     syyice(:,:,jl) , sxyice(:,:,jl) )
                   CALL lim_adv_y( zusnit, zvi_v, rzero, zsm, zs0ice(:,:,jl) , sxice (:,:,jl) , &
-                                                             sxxice(:,:,jl) , syice (:,:,jl) , & 
-                                                             syyice(:,:,jl) , sxyice(:,:,jl) )
+                     sxxice(:,:,jl) , syice (:,:,jl) , & 
+                     syyice(:,:,jl) , sxyice(:,:,jl) )
                   !--- snow volume  ---
                   CALL lim_adv_x( zusnit, zui_u, rone , zsm, zs0sn (:,:,jl) , sxsn  (:,:,jl) , &
-                                                             sxxsn (:,:,jl) , sysn  (:,:,jl) , &
-                                                             syysn (:,:,jl) , sxysn (:,:,jl) )
+                     sxxsn (:,:,jl) , sysn  (:,:,jl) , &
+                     syysn (:,:,jl) , sxysn (:,:,jl) )
                   CALL lim_adv_y( zusnit, zvi_v, rzero, zsm, zs0sn (:,:,jl) , sxsn  (:,:,jl) , &
-                                                             sxxsn (:,:,jl) , sysn  (:,:,jl) , &
-                                                             syysn (:,:,jl) , sxysn (:,:,jl) )
+                     sxxsn (:,:,jl) , sysn  (:,:,jl) , &
+                     syysn (:,:,jl) , sxysn (:,:,jl) )
                   !--- ice salinity ---
                   CALL lim_adv_x( zusnit, zui_u, rone , zsm, zs0sm (:,:,jl) , sxsal (:,:,jl) , &
-                                                             sxxsal(:,:,jl) , sysal (:,:,jl) , &
-                                                             syysal(:,:,jl) , sxysal(:,:,jl)  )
+                     sxxsal(:,:,jl) , sysal (:,:,jl) , &
+                     syysal(:,:,jl) , sxysal(:,:,jl)  )
                   CALL lim_adv_y( zusnit, zvi_v, rzero, zsm, zs0sm (:,:,jl) , sxsal (:,:,jl) , &
-                                                             sxxsal(:,:,jl) , sysal (:,:,jl) , &
-                                                             syysal(:,:,jl) , sxysal(:,:,jl)  )
+                     sxxsal(:,:,jl) , sysal (:,:,jl) , &
+                     syysal(:,:,jl) , sxysal(:,:,jl)  )
                   !--- ice age      ---     
                   CALL lim_adv_x( zusnit, zui_u, rone , zsm, zs0oi (:,:,jl) , sxage (:,:,jl) , &
-                                                             sxxage(:,:,jl) , syage (:,:,jl) , &
-                                                             syyage(:,:,jl) , sxyage(:,:,jl)  )
+                     sxxage(:,:,jl) , syage (:,:,jl) , &
+                     syyage(:,:,jl) , sxyage(:,:,jl)  )
                   CALL lim_adv_y( zusnit, zvi_v, rzero, zsm, zs0oi (:,:,jl) , sxage (:,:,jl) , &
-                                                             sxxage(:,:,jl) , syage (:,:,jl) , &
-                                                             syyage(:,:,jl) , sxyage(:,:,jl)  )
+                     sxxage(:,:,jl) , syage (:,:,jl) , &
+                     syyage(:,:,jl) , sxyage(:,:,jl)  )
                   !--- ice concentrations ---
                   CALL lim_adv_x( zusnit, zui_u, rone , zsm, zs0a  (:,:,jl) , sxa   (:,:,jl) , &
-                                                             sxxa  (:,:,jl) , sya   (:,:,jl) , & 
-                                                             syya  (:,:,jl) , sxya  (:,:,jl)  )
+                     sxxa  (:,:,jl) , sya   (:,:,jl) , & 
+                     syya  (:,:,jl) , sxya  (:,:,jl)  )
                   CALL lim_adv_y( zusnit, zvi_v, rzero, zsm, zs0a  (:,:,jl) , sxa   (:,:,jl) , & 
-                                                             sxxa  (:,:,jl) , sya   (:,:,jl) , & 
-                                                             syya  (:,:,jl) , sxya  (:,:,jl)  )
+                     sxxa  (:,:,jl) , sya   (:,:,jl) , & 
+                     syya  (:,:,jl) , sxya  (:,:,jl)  )
                   !--- ice / snow thermal energetic contents ---
                   CALL lim_adv_x( zusnit, zui_u, rone , zsm, zs0c0 (:,:,jl) , sxc0  (:,:,jl) , & 
-                                                             sxxc0 (:,:,jl) , syc0  (:,:,jl) , &
-                                                             syyc0 (:,:,jl) , sxyc0 (:,:,jl)  )
+                     sxxc0 (:,:,jl) , syc0  (:,:,jl) , &
+                     syyc0 (:,:,jl) , sxyc0 (:,:,jl)  )
                   CALL lim_adv_y( zusnit, zvi_v, rzero, zsm, zs0c0 (:,:,jl) , sxc0  (:,:,jl) , &
-                                                             sxxc0 (:,:,jl) , syc0  (:,:,jl) , &
-                                                             syyc0 (:,:,jl) , sxyc0 (:,:,jl)  )
+                     sxxc0 (:,:,jl) , syc0  (:,:,jl) , &
+                     syyc0 (:,:,jl) , sxyc0 (:,:,jl)  )
                   DO layer = 1, nlay_i
                      CALL lim_adv_x( zusnit, zui_u, rone , zsm, &
-                                                             zs0e(:,:,layer,jl) , sxe (:,:,layer,jl) , & 
-                                                             sxxe(:,:,layer,jl) , sye (:,:,layer,jl) , &
-                                                             syye(:,:,layer,jl) , sxye(:,:,layer,jl) )
+                        zs0e(:,:,layer,jl) , sxe (:,:,layer,jl) , & 
+                        sxxe(:,:,layer,jl) , sye (:,:,layer,jl) , &
+                        syye(:,:,layer,jl) , sxye(:,:,layer,jl) )
                      CALL lim_adv_y( zusnit, zvi_v, rzero, zsm, & 
-                                                             zs0e(:,:,layer,jl) , sxe (:,:,layer,jl) , & 
-                                                             sxxe(:,:,layer,jl) , sye (:,:,layer,jl) , &
-                                                             syye(:,:,layer,jl) , sxye(:,:,layer,jl) )
+                        zs0e(:,:,layer,jl) , sxe (:,:,layer,jl) , & 
+                        sxxe(:,:,layer,jl) , sye (:,:,layer,jl) , &
+                        syye(:,:,layer,jl) , sxye(:,:,layer,jl) )
                   END DO
                END DO
@@ -271,59 +273,59 @@
                !--- ice volume  ---
                CALL lim_adv_y( zusnit, zvi_v, rone , zsm, zs0ow (:,:) , sxopw (:,:) , &
-                                                          sxxopw(:,:) , syopw (:,:) , & 
-                                                          syyopw(:,:) , sxyopw(:,:) )
+                  sxxopw(:,:) , syopw (:,:) , & 
+                  syyopw(:,:) , sxyopw(:,:) )
                CALL lim_adv_x( zusnit, zui_u, rzero, zsm, zs0ow (:,:) , sxopw (:,:) , & 
-                                                          sxxopw(:,:) , syopw (:,:) , &
-                                                          syyopw(:,:) , sxyopw(:,:) )
+                  sxxopw(:,:) , syopw (:,:) , &
+                  syyopw(:,:) , sxyopw(:,:) )
                DO jl = 1, jpl
                   !--- ice volume  ---
                   CALL lim_adv_y( zusnit, zvi_v, rone , zsm, zs0ice(:,:,jl) , sxice (:,:,jl) , &
-                                                             sxxice(:,:,jl) , syice (:,:,jl) , & 
-                                                             syyice(:,:,jl) , sxyice(:,:,jl) )
+                     sxxice(:,:,jl) , syice (:,:,jl) , & 
+                     syyice(:,:,jl) , sxyice(:,:,jl) )
                   CALL lim_adv_x( zusnit, zui_u, rzero, zsm, zs0ice(:,:,jl) , sxice (:,:,jl) , & 
-                                                             sxxice(:,:,jl) , syice (:,:,jl) , &
-                                                             syyice(:,:,jl) , sxyice(:,:,jl) )
+                     sxxice(:,:,jl) , syice (:,:,jl) , &
+                     syyice(:,:,jl) , sxyice(:,:,jl) )
                   !--- snow volume  ---
                   CALL lim_adv_y( zusnit, zvi_v, rone , zsm, zs0sn (:,:,jl) , sxsn  (:,:,jl) , & 
-                                                             sxxsn (:,:,jl) , sysn  (:,:,jl) , & 
-                                                             syysn (:,:,jl) , sxysn (:,:,jl) )
+                     sxxsn (:,:,jl) , sysn  (:,:,jl) , & 
+                     syysn (:,:,jl) , sxysn (:,:,jl) )
                   CALL lim_adv_x( zusnit, zui_u, rzero, zsm, zs0sn (:,:,jl) , sxsn  (:,:,jl) , & 
-                                                             sxxsn (:,:,jl) , sysn  (:,:,jl) , & 
-                                                             syysn (:,:,jl) , sxysn (:,:,jl) )
+                     sxxsn (:,:,jl) , sysn  (:,:,jl) , & 
+                     syysn (:,:,jl) , sxysn (:,:,jl) )
                   !--- ice salinity ---
                   CALL lim_adv_y( zusnit, zvi_v, rone , zsm, zs0sm (:,:,jl) , sxsal (:,:,jl) , &
-                                                             sxxsal(:,:,jl) , sysal (:,:,jl) , &
-                                                             syysal(:,:,jl) , sxysal(:,:,jl) )
+                     sxxsal(:,:,jl) , sysal (:,:,jl) , &
+                     syysal(:,:,jl) , sxysal(:,:,jl) )
                   CALL lim_adv_x( zusnit, zui_u, rzero, zsm, zs0sm (:,:,jl) , sxsal (:,:,jl) , &
-                                                             sxxsal(:,:,jl) , sysal (:,:,jl) , &
-                                                             syysal(:,:,jl) , sxysal(:,:,jl) )
+                     sxxsal(:,:,jl) , sysal (:,:,jl) , &
+                     syysal(:,:,jl) , sxysal(:,:,jl) )
                   !--- ice age      ---
                   CALL lim_adv_y( zusnit, zvi_v, rone , zsm, zs0oi (:,:,jl) , sxage (:,:,jl) , &
-                                                             sxxage(:,:,jl) , syage (:,:,jl) , & 
-                                                             syyage(:,:,jl) , sxyage(:,:,jl)  )
+                     sxxage(:,:,jl) , syage (:,:,jl) , & 
+                     syyage(:,:,jl) , sxyage(:,:,jl)  )
                   CALL lim_adv_x( zusnit, zui_u, rzero, zsm, zs0oi (:,:,jl) , sxage (:,:,jl) , &
-                                                             sxxage(:,:,jl) , syage (:,:,jl) , &
-                                                             syyage(:,:,jl) , sxyage(:,:,jl)   )
+                     sxxage(:,:,jl) , syage (:,:,jl) , &
+                     syyage(:,:,jl) , sxyage(:,:,jl)   )
                   !--- ice concentration ---
                   CALL lim_adv_y( zusnit, zvi_v, rone , zsm, zs0a  (:,:,jl) , sxa   (:,:,jl) , & 
-                                                             sxxa  (:,:,jl) , sya   (:,:,jl) , & 
-                                                             syya  (:,:,jl) , sxya  (:,:,jl)  )
+                     sxxa  (:,:,jl) , sya   (:,:,jl) , & 
+                     syya  (:,:,jl) , sxya  (:,:,jl)  )
                   CALL lim_adv_x( zusnit, zui_u, rzero, zsm, zs0a  (:,:,jl) , sxa   (:,:,jl) , & 
-                                                             sxxa  (:,:,jl) , sya   (:,:,jl) , & 
-                                                             syya  (:,:,jl) , sxya  (:,:,jl)  )
+                     sxxa  (:,:,jl) , sya   (:,:,jl) , & 
+                     syya  (:,:,jl) , sxya  (:,:,jl)  )
                   !--- ice / snow thermal energetic contents ---
                   CALL lim_adv_y( zusnit, zvi_v, rone , zsm, zs0c0 (:,:,jl) , sxc0  (:,:,jl) , & 
-                                                             sxxc0 (:,:,jl) , syc0  (:,:,jl) , &
-                                                             syyc0 (:,:,jl) , sxyc0 (:,:,jl)  )
+                     sxxc0 (:,:,jl) , syc0  (:,:,jl) , &
+                     syyc0 (:,:,jl) , sxyc0 (:,:,jl)  )
                   CALL lim_adv_x( zusnit, zui_u, rzero, zsm, zs0c0 (:,:,jl) , sxc0  (:,:,jl) , &
-                                                             sxxc0 (:,:,jl) , syc0  (:,:,jl) , &
-                                                             syyc0 (:,:,jl) , sxyc0 (:,:,jl)  )
+                     sxxc0 (:,:,jl) , syc0  (:,:,jl) , &
+                     syyc0 (:,:,jl) , sxyc0 (:,:,jl)  )
                   DO layer = 1, nlay_i
                      CALL lim_adv_y( zusnit, zvi_v, rone , zsm, zs0e(:,:,layer,jl) , &
-                                     sxe (:,:,layer,jl) , sxxe (:,:,layer,jl) , sye (:,:,layer,jl) , &
-                                     syye (:,:,layer,jl), sxye (:,:,layer,jl) )
+                        sxe (:,:,layer,jl) , sxxe (:,:,layer,jl) , sye (:,:,layer,jl) , &
+                        syye (:,:,layer,jl), sxye (:,:,layer,jl) )
                      CALL lim_adv_x( zusnit, zui_u, rzero, zsm, zs0e(:,:,layer,jl) , &
-                                     sxe (:,:,layer,jl) , sxxe (:,:,layer,jl) , sye (:,:,layer,jl) , &
-                                     syye (:,:,layer,jl), sxye (:,:,layer,jl)  )
+                        sxe (:,:,layer,jl) , sxxe (:,:,layer,jl) , sye (:,:,layer,jl) , &
+                        syye (:,:,layer,jl), sxye (:,:,layer,jl)  )
                   END DO
 
@@ -349,11 +351,11 @@
          END DO
 
-!------------------------------------------------------------------------------!
-! 4) Diffusion of Ice fields                  
-!------------------------------------------------------------------------------!
-
-      !------------------------------------
-      ! 4.1) diffusion of open water area
-      !------------------------------------
+         !------------------------------------------------------------------------------!
+         ! 4) Diffusion of Ice fields                  
+         !------------------------------------------------------------------------------!
+
+         !------------------------------------
+         ! 4.1) diffusion of open water area
+         !------------------------------------
 
          ! Compute total ice fraction
@@ -364,5 +366,5 @@
                   zs0at (ji,jj) = zs0at(ji,jj) + zs0a(ji,jj,jl) !
                END DO
-            END DO 
+            END DO
          END DO
 
@@ -380,10 +382,10 @@
          CALL lim_hdf( zs0ow (:,:) )
 
-      !----------------------------------------
-      ! 4.2) Diffusion of other ice variables
-      !----------------------------------------
-         DO jl = 1, jpl
-
-         ! Masked eddy diffusivity coefficient at ocean U- and V-points
+         !----------------------------------------
+         ! 4.2) Diffusion of other ice variables
+         !----------------------------------------
+         DO jl = 1, jpl
+
+            ! Masked eddy diffusivity coefficient at ocean U- and V-points
             DO jj = 1, jpjm1          ! NB: has not to be defined on jpj line and jpi row
                DO ji = 1 , fs_jpim1   ! vector opt.
@@ -406,7 +408,7 @@
          END DO !jl
 
-      !-----------------------------------------
-      ! 4.3) Remultiply everything by ice area
-      !-----------------------------------------
+         !-----------------------------------------
+         ! 4.3) Remultiply everything by ice area
+         !-----------------------------------------
          zs0ow(:,:) = MAX(rzero, zs0ow(:,:) * area(:,:) )
          DO jl = 1, jpl
@@ -422,11 +424,11 @@
          END DO ! jl
 
-!------------------------------------------------------------------------------!
-! 5) Update and limit ice properties after transport                           
-!------------------------------------------------------------------------------!
-
-      !--------------------------------------------------
-      ! 5.1) Recover mean values over the grid squares.
-      !--------------------------------------------------
+         !------------------------------------------------------------------------------!
+         ! 5) Update and limit ice properties after transport                           
+         !------------------------------------------------------------------------------!
+
+         !--------------------------------------------------
+         ! 5.1) Recover mean values over the grid squares.
+         !--------------------------------------------------
 
          DO jl = 1, jpl
@@ -435,5 +437,5 @@
                   DO ji = 1, jpi
                      zs0e (ji,jj,jk,jl) =  &
-                     MAX( rzero, zs0e (ji,jj,jk,jl) / area(ji,jj) )
+                        MAX( rzero, zs0e (ji,jj,jk,jl) / area(ji,jj) )
                   END DO
                END DO
@@ -446,5 +448,5 @@
             END DO
          END DO
-         
+
          zs0at(:,:) = 0.0
          DO jl = 1, jpl
@@ -462,7 +464,7 @@
          END DO
 
-      !---------------------------------------------------------
-      ! 5.2) Snow thickness, Ice thickness, Ice concentrations
-      !---------------------------------------------------------
+         !---------------------------------------------------------
+         ! 5.2) Snow thickness, Ice thickness, Ice concentrations
+         !---------------------------------------------------------
 
          DO jj = 1, jpj
@@ -501,7 +503,7 @@
          END DO
 
-      !----------------------
-      ! 5.3) Ice properties
-      !----------------------
+         !----------------------
+         ! 5.3) Ice properties
+         !----------------------
 
          zbigval         =  1.0d+13
@@ -521,12 +523,12 @@
                   ! Ice salinity and age
                   zsal            = MAX( MIN( (rhoic-rhosn)/rhoic*sss_m(ji,jj)  , &
-                                            zusvoic * zs0sm(ji,jj,jl) ), s_i_min ) * &
-                                            v_i(ji,jj,jl)
+                     zusvoic * zs0sm(ji,jj,jl) ), s_i_min ) * &
+                     v_i(ji,jj,jl)
                   IF ( ( num_sal .EQ. 2 ) .OR. ( num_sal .EQ. 4 ) ) & 
                      smv_i(ji,jj,jl) = zindic*zsal + (1.0-zindic)*0.0
 
                   zage            = MAX( MIN( zbigval, zs0oi(ji,jj,jl) / & 
-                                              MAX( a_i(ji,jj,jl), epsi16 )  ), 0.0 ) * &
-                                              a_i(ji,jj,jl)
+                     MAX( a_i(ji,jj,jl), epsi16 )  ), 0.0 ) * &
+                     a_i(ji,jj,jl)
                   oa_i (ji,jj,jl)  = zindic*zage 
 
@@ -616,5 +618,5 @@
          WRITE(numout,*) 
       ENDIF
-            
+
    END SUBROUTINE lim_trp_init
 
