Index: branches/2012/dev_v3_4_STABLE_2012/NEMOGCM/NEMO/OPATAM_SRC/tstool_tam.F90
===================================================================
--- branches/2012/dev_v3_4_STABLE_2012/NEMOGCM/NEMO/OPATAM_SRC/tstool_tam.F90	(revision 5156)
+++ branches/2012/dev_v3_4_STABLE_2012/NEMOGCM/NEMO/OPATAM_SRC/tstool_tam.F90	(revision 5168)
@@ -14,9 +14,6 @@
       & e1u,                 &
       & e2u,                 &
-#if defined key_zco
       & e3t_0,               &
-#else
       & e3u,                 &
-#endif
       & umask,               &
       & mig,                 &
Index: branches/2012/dev_v3_4_STABLE_2012/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90
===================================================================
--- branches/2012/dev_v3_4_STABLE_2012/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90	(revision 5156)
+++ branches/2012/dev_v3_4_STABLE_2012/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90	(revision 5168)
@@ -9,5 +9,5 @@
    !!                 ! 2007-04  (A. Weaver)  Merge with OPAVAR/NEMOVAR
    !!   NEMO     3.3  ! 2010-05  (D. Lea)  Update to work with NEMO v3.2
-   !!             -   ! 2010-05  (D. Lea)  add calc_month_len routine based on day_init 
+   !!             -   ! 2010-05  (D. Lea)  add calc_month_len routine based on day_init
    !!            3.4  ! 2012-10  (A. Weaver and K. Mogensen) Fix for direct initialization
    !!----------------------------------------------------------------------
@@ -43,5 +43,6 @@
    IMPLICIT NONE
    PRIVATE
-   
+
+   PUBLIC asm_inc_rea_nam  !: Read namelist for ASM
    PUBLIC   asm_inc_init   !: Initialize the increment arrays and IAU weights
    PUBLIC   calc_date      !: Compute the calendar date YYYYMMDD on a given step
@@ -70,5 +71,5 @@
    REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE ::   u_bkg   , v_bkg      !: Background u- & v- velocity components
    REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE ::   t_bkginc, s_bkginc   !: Increment to the background T & S
-   REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE ::   u_bkginc, v_bkginc   !: Increment to the u- & v-components 
+   REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE ::   u_bkginc, v_bkginc   !: Increment to the u- & v-components
    REAL(wp), PUBLIC, DIMENSION(:)    , ALLOCATABLE ::   wgtiau               !: IAU weights for each time step
 #if defined key_asminc
@@ -78,9 +79,9 @@
    INTEGER , PUBLIC ::   nitbkg      !: Time step of the background state used in the Jb term
    INTEGER , PUBLIC ::   nitdin      !: Time step of the background state for direct initialization
-   INTEGER , PUBLIC ::   nitiaustr   !: Time step of the start of the IAU interval 
+   INTEGER , PUBLIC ::   nitiaustr   !: Time step of the start of the IAU interval
    INTEGER , PUBLIC ::   nitiaufin   !: Time step of the end of the IAU interval
-   ! 
+   !
    INTEGER , PUBLIC ::   niaufn      !: Type of IAU weighing function: = 0   Constant weighting
-   !                                 !: = 1   Linear hat-like, centred in middle of IAU interval 
+   !                                 !: = 1   Linear hat-like, centred in middle of IAU interval
    REAL(wp), PUBLIC ::   salfixmin   !: Ensure that the salinity is larger than this  value if (ln_salfix)
 
@@ -100,38 +101,5 @@
 CONTAINS
 
-   SUBROUTINE asm_inc_init
-      !!----------------------------------------------------------------------
-      !!                    ***  ROUTINE asm_inc_init  ***
-      !!          
-      !! ** Purpose : Initialize the assimilation increment and IAU weights.
-      !!
-      !! ** Method  : Initialize the assimilation increment and IAU weights.
-      !!
-      !! ** Action  : 
-      !!----------------------------------------------------------------------
-      !!
-      !!
-      INTEGER :: ji,jj,jk
-      INTEGER :: jt
-      INTEGER :: imid
-      INTEGER :: inum
-      INTEGER :: iiauper         ! Number of time steps in the IAU period
-      INTEGER :: icycper         ! Number of time steps in the cycle
-      INTEGER :: iitend_date     ! Date YYYYMMDD of final time step
-      INTEGER :: iitbkg_date     ! Date YYYYMMDD of background time step for Jb term
-      INTEGER :: iitdin_date     ! Date YYYYMMDD of background time step for DI
-      INTEGER :: iitiaustr_date  ! Date YYYYMMDD of IAU interval start time step
-      INTEGER :: iitiaufin_date  ! Date YYYYMMDD of IAU interval final time step
-
-      REAL(wp) :: znorm        ! Normalization factor for IAU weights
-      REAL(wp) :: ztotwgt      ! Value of time-integrated IAU weights 
-                               ! (should be equal to one)
-      REAL(wp) :: z_inc_dateb  ! Start date of interval on which increment is valid
-      REAL(wp) :: z_inc_datef  ! End date of interval on which increment is valid
-      REAL(wp) :: zdate_bkg    ! Date in background state file for DI
-      REAL(wp) :: zdate_inc    ! Time axis in increments file
-
-      REAL(wp), POINTER, DIMENSION(:,:) :: hdiv
-      !!
+   SUBROUTINE asm_inc_rea_nam
       NAMELIST/nam_asminc/ ln_bkgwri,                                      &
          &                 ln_trainc, ln_dyninc, ln_sshinc,                &
@@ -158,5 +126,5 @@
       salfixmin = -9999
       nitbkg    = 0
-      nitdin    = 0      
+      nitdin    = 0
       nitiaustr = 1
       nitiaufin = 150      ! = 10 days with ORCA2
@@ -187,4 +155,40 @@
          WRITE(numout,*) '      Minimum salinity after applying the increments           salfixmin = ', salfixmin
       ENDIF
+   END SUBROUTINE
+
+   SUBROUTINE asm_inc_init
+      !!----------------------------------------------------------------------
+      !!                    ***  ROUTINE asm_inc_init  ***
+      !!
+      !! ** Purpose : Initialize the assimilation increment and IAU weights.
+      !!
+      !! ** Method  : Initialize the assimilation increment and IAU weights.
+      !!
+      !! ** Action  :
+      !!----------------------------------------------------------------------
+      !!
+      !!
+      INTEGER :: ji,jj,jk
+      INTEGER :: jt
+      INTEGER :: imid
+      INTEGER :: inum
+      INTEGER :: iiauper         ! Number of time steps in the IAU period
+      INTEGER :: icycper         ! Number of time steps in the cycle
+      INTEGER :: iitend_date     ! Date YYYYMMDD of final time step
+      INTEGER :: iitbkg_date     ! Date YYYYMMDD of background time step for Jb term
+      INTEGER :: iitdin_date     ! Date YYYYMMDD of background time step for DI
+      INTEGER :: iitiaustr_date  ! Date YYYYMMDD of IAU interval start time step
+      INTEGER :: iitiaufin_date  ! Date YYYYMMDD of IAU interval final time step
+
+      REAL(wp) :: znorm        ! Normalization factor for IAU weights
+      REAL(wp) :: ztotwgt      ! Value of time-integrated IAU weights
+                               ! (should be equal to one)
+      REAL(wp) :: z_inc_dateb  ! Start date of interval on which increment is valid
+      REAL(wp) :: z_inc_datef  ! End date of interval on which increment is valid
+      REAL(wp) :: zdate_bkg    ! Date in background state file for DI
+      REAL(wp) :: zdate_inc    ! Time axis in increments file
+
+      REAL(wp), POINTER, DIMENSION(:,:) :: hdiv
+      !!
 
       nitbkg_r    = nitbkg    + nit000 - 1  ! Background time referenced to nit000
@@ -293,5 +297,5 @@
 
             !---------------------------------------------------------
-            ! Constant IAU forcing 
+            ! Constant IAU forcing
             !---------------------------------------------------------
 
@@ -303,5 +307,5 @@
 
             !---------------------------------------------------------
-            ! Linear hat-like, centred in middle of IAU interval 
+            ! Linear hat-like, centred in middle of IAU interval
             !---------------------------------------------------------
 
@@ -309,5 +313,5 @@
             znorm = 0.0
             IF ( MOD( iiauper, 2 ) == 0 ) THEN  ! Even number of time steps in IAU interval
-               imid = iiauper / 2 
+               imid = iiauper / 2
                DO jt = 1, imid
                   znorm = znorm + REAL( jt )
@@ -315,5 +319,5 @@
                znorm = 2.0 * znorm
             ELSE                               ! Odd number of time steps in IAU interval
-               imid = ( iiauper + 1 ) / 2        
+               imid = ( iiauper + 1 ) / 2
                DO jt = 1, imid - 1
                   znorm = znorm + REAL( jt )
@@ -342,11 +346,11 @@
              DO jt = 1, icycper
                 ztotwgt = ztotwgt + wgtiau(jt)
-                WRITE(numout,*) '         ', jt, '       ', wgtiau(jt) 
-             END DO   
+                WRITE(numout,*) '         ', jt, '       ', wgtiau(jt)
+             END DO
              WRITE(numout,*) '         ==================================='
              WRITE(numout,*) '         Time-integrated weight = ', ztotwgt
              WRITE(numout,*) '         ==================================='
           ENDIF
-         
+
       ENDIF
 
@@ -381,5 +385,5 @@
          CALL iom_open( c_asminc, inum )
 
-         CALL iom_get( inum, 'time', zdate_inc ) 
+         CALL iom_get( inum, 'time', zdate_inc )
 
          CALL iom_get( inum, 'z_inc_dateb', z_inc_dateb )
@@ -389,5 +393,5 @@
 
          IF(lwp) THEN
-            WRITE(numout,*) 
+            WRITE(numout,*)
             WRITE(numout,*) 'asm_inc_init : Assimilation increments valid ', &
                &            ' between dates ', NINT( z_inc_dateb ),' and ',  &
@@ -405,5 +409,5 @@
             &                ' not agree with Direct Initialization time' )
 
-         IF ( ln_trainc ) THEN   
+         IF ( ln_trainc ) THEN
             CALL iom_get( inum, jpdom_autoglo, 'bckint', t_bkginc, 1 )
             CALL iom_get( inum, jpdom_autoglo, 'bckins', s_bkginc, 1 )
@@ -417,7 +421,7 @@
          ENDIF
 
-         IF ( ln_dyninc ) THEN   
-            CALL iom_get( inum, jpdom_autoglo, 'bckinu', u_bkginc, 1 )              
-            CALL iom_get( inum, jpdom_autoglo, 'bckinv', v_bkginc, 1 )              
+         IF ( ln_dyninc ) THEN
+            CALL iom_get( inum, jpdom_autoglo, 'bckinu', u_bkginc, 1 )
+            CALL iom_get( inum, jpdom_autoglo, 'bckinv', v_bkginc, 1 )
             ! Apply the masks
             u_bkginc(:,:,:) = u_bkginc(:,:,:) * umask(:,:,:)
@@ -428,5 +432,5 @@
             WHERE( ABS( v_bkginc(:,:,:) ) > 1.0e+10 ) v_bkginc(:,:,:) = 0.0
          ENDIF
-        
+
          IF ( ln_sshinc ) THEN
             CALL iom_get( inum, jpdom_autoglo, 'bckineta', ssh_bkginc, 1 )
@@ -448,5 +452,5 @@
 
          CALL iom_close( inum )
- 
+
       ENDIF
 
@@ -457,5 +461,5 @@
       IF ( ln_dyninc .AND. nn_divdmp > 0 ) THEN
 
-         CALL wrk_alloc(jpi,jpj,hdiv) 
+         CALL wrk_alloc(jpi,jpj,hdiv)
 
          DO  jt = 1, nn_divdmp
@@ -482,8 +486,8 @@
                      u_bkginc(ji,jj,jk) = u_bkginc(ji,jj,jk) + 0.2_wp * ( e1t(ji+1,jj)*e2t(ji+1,jj) * hdiv(ji+1,jj)   &
                                                                         - e1t(ji  ,jj)*e2t(ji  ,jj) * hdiv(ji  ,jj) ) &
-                                                                      / e1u(ji,jj) * umask(ji,jj,jk) 
+                                                                      / e1u(ji,jj) * umask(ji,jj,jk)
                      v_bkginc(ji,jj,jk) = v_bkginc(ji,jj,jk) + 0.2_wp * ( e1t(ji,jj+1)*e2t(ji,jj+1) * hdiv(ji,jj+1)   &
                                                                         - e1t(ji,jj  )*e2t(ji,jj  ) * hdiv(ji,jj  ) ) &
-                                                                      / e2v(ji,jj) * vmask(ji,jj,jk) 
+                                                                      / e2v(ji,jj) * vmask(ji,jj,jk)
                   END DO
                END DO
@@ -493,5 +497,5 @@
          END DO
 
-         CALL wrk_dealloc(jpi,jpj,hdiv) 
+         CALL wrk_dealloc(jpi,jpj,hdiv)
 
       ENDIF
@@ -523,8 +527,8 @@
          CALL iom_open( c_asmdin, inum )
 
-         CALL iom_get( inum, 'rdastp', zdate_bkg ) 
-        
+         CALL iom_get( inum, 'rdastp', zdate_bkg )
+
          IF(lwp) THEN
-            WRITE(numout,*) 
+            WRITE(numout,*)
             WRITE(numout,*) 'asm_inc_init : Assimilation background state valid at : ', &
                &  NINT( zdate_bkg )
@@ -536,5 +540,5 @@
             &                ' not agree with Direct Initialization time' )
 
-         IF ( ln_trainc ) THEN   
+         IF ( ln_trainc ) THEN
             CALL iom_get( inum, jpdom_autoglo, 'tn', t_bkg )
             CALL iom_get( inum, jpdom_autoglo, 'sn', s_bkg )
@@ -543,5 +547,5 @@
          ENDIF
 
-         IF ( ln_dyninc ) THEN   
+         IF ( ln_dyninc ) THEN
             CALL iom_get( inum, jpdom_autoglo, 'un', u_bkg )
             CALL iom_get( inum, jpdom_autoglo, 'vn', v_bkg )
@@ -549,5 +553,5 @@
             v_bkg(:,:,:) = v_bkg(:,:,:) * vmask(:,:,:)
          ENDIF
-        
+
          IF ( ln_sshinc ) THEN
             CALL iom_get( inum, jpdom_autoglo, 'sshn', ssh_bkg )
@@ -565,10 +569,10 @@
       !!----------------------------------------------------------------------
       !!                    ***  ROUTINE calc_date  ***
-      !!          
+      !!
       !! ** Purpose : Compute the calendar date YYYYMMDD at a given time step.
       !!
       !! ** Method  : Compute the calendar date YYYYMMDD at a given time step.
       !!
-      !! ** Action  : 
+      !! ** Action  :
       !!----------------------------------------------------------------------
       INTEGER, INTENT(IN) :: kit000  ! Initial time step
@@ -595,5 +599,5 @@
       iyea0 =   kdate0 / 10000
       imon0 = ( kdate0 - ( iyea0 * 10000 ) ) / 100
-      iday0 =   kdate0 - ( iyea0 * 10000 ) - ( imon0 * 100 ) 
+      iday0 =   kdate0 - ( iyea0 * 10000 ) - ( imon0 * 100 )
 
       ! Check that kt >= kit000 - 1
@@ -610,5 +614,5 @@
       ! Compute the number of days from the initial date
       idaystp = INT( REAL( kt - kit000 ) * rdt / 86400. )
-   
+
       iday    = iday0
       imon    = imon0
@@ -628,5 +632,5 @@
             iyea = iyea + 1
             CALL calc_month_len( iyea, imonth_len )  ! update month lengths
-         ENDIF                 
+         ENDIF
          idaycnt = idaycnt + 1
       END DO
@@ -640,8 +644,8 @@
       !!----------------------------------------------------------------------
       !!                    ***  ROUTINE calc_month_len  ***
-      !!          
+      !!
       !! ** Purpose : Compute the number of days in a months given a year.
       !!
-      !! ** Method  : 
+      !! ** Method  :
       !!----------------------------------------------------------------------
       INTEGER, DIMENSION(12) ::   imonth_len    !: length in days of the months of the current year
@@ -650,5 +654,5 @@
       !
       ! length of the month of the current year (from nleapy, read in namelist)
-      IF ( nleapy < 2 ) THEN 
+      IF ( nleapy < 2 ) THEN
          imonth_len(:) = (/ 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31 /)
          IF ( nleapy == 1 ) THEN   ! we are using calendar with leap years
@@ -667,10 +671,10 @@
       !!----------------------------------------------------------------------
       !!                    ***  ROUTINE tra_asm_inc  ***
-      !!          
+      !!
       !! ** Purpose : Apply the tracer (T and S) assimilation increments
       !!
       !! ** Method  : Direct initialization or Incremental Analysis Updating
       !!
-      !! ** Action  : 
+      !! ** Action  :
       !!----------------------------------------------------------------------
       INTEGER, INTENT(IN) :: kt               ! Current time step
@@ -682,20 +686,20 @@
       !!----------------------------------------------------------------------
 
-      ! freezing point calculation taken from oc_fz_pt (but calculated for all depths) 
-      ! used to prevent the applied increments taking the temperature below the local freezing point 
-
-#if defined key_cice 
+      ! freezing point calculation taken from oc_fz_pt (but calculated for all depths)
+      ! used to prevent the applied increments taking the temperature below the local freezing point
+
+#if defined key_cice
         fzptnz(:,:,:) = -1.8_wp
-#else 
+#else
         DO jk = 1, jpk
            DO jj = 1, jpj
               DO ji = 1, jpk
-                 fzptnz (ji,jj,jk) = ( -0.0575_wp + 1.710523e-3_wp * SQRT( tsn(ji,jj,jk,jp_sal) )                   & 
-                                                  - 2.154996e-4_wp *       tsn(ji,jj,jk,jp_sal)   ) * tsn(ji,jj,jk,jp_sal)  & 
-                                                  - 7.53e-4_wp * fsdepw(ji,jj,jk)       ! (pressure in dbar) 
+                 fzptnz (ji,jj,jk) = ( -0.0575_wp + 1.710523e-3_wp * SQRT( tsn(ji,jj,jk,jp_sal) )                   &
+                                                  - 2.154996e-4_wp *       tsn(ji,jj,jk,jp_sal)   ) * tsn(ji,jj,jk,jp_sal)  &
+                                                  - 7.53e-4_wp * fsdepw(ji,jj,jk)       ! (pressure in dbar)
               END DO
            END DO
         END DO
-#endif 
+#endif
 
       IF ( ln_asmiau ) THEN
@@ -711,5 +715,5 @@
 
             IF(lwp) THEN
-               WRITE(numout,*) 
+               WRITE(numout,*)
                WRITE(numout,*) 'tra_asm_inc : Tracer IAU at time step = ', &
                   &  kt,' with IAU weight = ', wgtiau(it)
@@ -722,9 +726,9 @@
                   ! Do not apply negative increments if the temperature will fall below freezing
                   WHERE(t_bkginc(:,:,jk) > 0.0_wp .OR. &
-                     &   tsn(:,:,jk,jp_tem) + tsa(:,:,jk,jp_tem) + t_bkginc(:,:,jk) * wgtiau(it) > fzptnz(:,:,jk) ) 
-                     tsa(:,:,jk,jp_tem) = tsa(:,:,jk,jp_tem) + t_bkginc(:,:,jk) * zincwgt  
+                     &   t_bkg(:,:,jk) + tsa(:,:,jk,jp_tem) + t_bkginc(:,:,jk) * wgtiau(it) > fzptnz(:,:,jk) )
+                     tsa(:,:,jk,jp_tem) = tsa(:,:,jk,jp_tem) + t_bkginc(:,:,jk) * zincwgt
                   END WHERE
                ELSE
-                  tsa(:,:,jk,jp_tem) = tsa(:,:,jk,jp_tem) + t_bkginc(:,:,jk) * zincwgt  
+                  tsa(:,:,jk,jp_tem) = tsa(:,:,jk,jp_tem) + t_bkginc(:,:,jk) * zincwgt
                ENDIF
                IF (ln_salfix) THEN
@@ -732,5 +736,5 @@
                   ! minimum value salfixmin
                   WHERE(s_bkginc(:,:,jk) > 0.0_wp .OR. &
-                     &   tsn(:,:,jk,jp_sal) + tsa(:,:,jk,jp_sal) + s_bkginc(:,:,jk) * wgtiau(it) > salfixmin ) 
+                     &  s_bkg(:,:,jk) + tsa(:,:,jk,jp_sal) + s_bkginc(:,:,jk) * wgtiau(it) > salfixmin )
                      tsa(:,:,jk,jp_sal) = tsa(:,:,jk,jp_sal) + s_bkginc(:,:,jk) * zincwgt
                   END WHERE
@@ -753,5 +757,5 @@
          ! Direct Initialization
          !--------------------------------------------------------------------
-            
+
          IF ( kt == nitdin_r ) THEN
 
@@ -762,9 +766,9 @@
                ! Do not apply negative increments if the temperature will fall below freezing
                WHERE(t_bkginc(:,:,:) > 0.0_wp .OR. &
-                  &   tsn(:,:,:,jp_tem) + t_bkginc(:,:,:) > fzptnz(:,:,:) ) 
-                  tsn(:,:,:,jp_tem) = t_bkg(:,:,:) + t_bkginc(:,:,:)   
+                  &   tsn(:,:,:,jp_tem) + t_bkginc(:,:,:) > fzptnz(:,:,:) )
+                  tsn(:,:,:,jp_tem) = t_bkg(:,:,:) + t_bkginc(:,:,:)
                END WHERE
             ELSE
-               tsn(:,:,:,jp_tem) = t_bkg(:,:,:) + t_bkginc(:,:,:)   
+               tsn(:,:,:,jp_tem) = t_bkg(:,:,:) + t_bkginc(:,:,:)
             ENDIF
             IF (ln_salfix) THEN
@@ -772,9 +776,9 @@
                ! minimum value salfixmin
                WHERE(s_bkginc(:,:,:) > 0.0_wp .OR. &
-                  &   tsn(:,:,:,jp_sal) + s_bkginc(:,:,:) > salfixmin ) 
-                  tsn(:,:,:,jp_sal) = s_bkg(:,:,:) + s_bkginc(:,:,:)   
+                  &   tsn(:,:,:,jp_sal) + s_bkginc(:,:,:) > salfixmin )
+                  tsn(:,:,:,jp_sal) = s_bkg(:,:,:) + s_bkginc(:,:,:)
                END WHERE
             ELSE
-               tsn(:,:,:,jp_sal) = s_bkg(:,:,:) + s_bkginc(:,:,:)   
+               tsn(:,:,:,jp_sal) = s_bkg(:,:,:) + s_bkginc(:,:,:)
             ENDIF
 
@@ -782,5 +786,5 @@
 
             CALL eos( tsb, rhd, rhop )                ! Before potential and in situ densities
-         
+
             IF( ln_zps .AND. .NOT. lk_c1d ) &
                &  CALL zps_hde( nit000, jpts, tsb, &  ! Partial steps: before horizontal derivative
@@ -797,5 +801,5 @@
             DEALLOCATE( s_bkg    )
          ENDIF
-         !  
+         !
       ENDIF
       ! Perhaps the following call should be in step
@@ -808,10 +812,10 @@
       !!----------------------------------------------------------------------
       !!                    ***  ROUTINE dyn_asm_inc  ***
-      !!          
+      !!
       !! ** Purpose : Apply the dynamics (u and v) assimilation increments.
       !!
       !! ** Method  : Direct initialization or Incremental Analysis Updating.
       !!
-      !! ** Action  : 
+      !! ** Action  :
       !!----------------------------------------------------------------------
       INTEGER, INTENT(IN) :: kt   ! Current time step
@@ -834,5 +838,5 @@
 
             IF(lwp) THEN
-               WRITE(numout,*) 
+               WRITE(numout,*)
                WRITE(numout,*) 'dyn_asm_inc : Dynamics IAU at time step = ', &
                   &  kt,' with IAU weight = ', wgtiau(it)
@@ -845,5 +849,5 @@
                va(:,:,jk) = va(:,:,jk) + v_bkginc(:,:,jk) * zincwgt
             END DO
-           
+
             IF ( kt == nitiaufin_r ) THEN
                DEALLOCATE( u_bkginc )
@@ -853,10 +857,10 @@
          ENDIF
 
-      ELSEIF ( ln_asmdin ) THEN 
+      ELSEIF ( ln_asmdin ) THEN
 
          !--------------------------------------------------------------------
          ! Direct Initialization
          !--------------------------------------------------------------------
-         
+
          IF ( kt == nitdin_r ) THEN
 
@@ -865,9 +869,9 @@
             ! Initialize the now fields with the background + increment
             un(:,:,:) = u_bkg(:,:,:) + u_bkginc(:,:,:)
-            vn(:,:,:) = v_bkg(:,:,:) + v_bkginc(:,:,:)  
+            vn(:,:,:) = v_bkg(:,:,:) + v_bkginc(:,:,:)
 
             ub(:,:,:) = un(:,:,:)         ! Update before fields
             vb(:,:,:) = vn(:,:,:)
- 
+
             DEALLOCATE( u_bkg    )
             DEALLOCATE( v_bkg    )
@@ -884,10 +888,10 @@
       !!----------------------------------------------------------------------
       !!                    ***  ROUTINE ssh_asm_inc  ***
-      !!          
+      !!
       !! ** Purpose : Apply the sea surface height assimilation increment.
       !!
       !! ** Method  : Direct initialization or Incremental Analysis Updating.
       !!
-      !! ** Action  : 
+      !! ** Action  :
       !!----------------------------------------------------------------------
       INTEGER, INTENT(IN) :: kt   ! Current time step
@@ -910,5 +914,5 @@
 
             IF(lwp) THEN
-               WRITE(numout,*) 
+               WRITE(numout,*)
                WRITE(numout,*) 'ssh_asm_inc : SSH IAU at time step = ', &
                   &  kt,' with IAU weight = ', wgtiau(it)
@@ -938,8 +942,8 @@
 
             ! Initialize the now fields the background + increment
-            sshn(:,:) = ssh_bkg(:,:) + ssh_bkginc(:,:)  
+            sshn(:,:) = ssh_bkg(:,:) + ssh_bkginc(:,:)
 
             ! Update before fields
-            sshb(:,:) = sshn(:,:)         
+            sshb(:,:) = sshn(:,:)
 
             IF( lk_vvl ) THEN
@@ -961,10 +965,10 @@
       !!----------------------------------------------------------------------
       !!                    ***  ROUTINE seaice_asm_inc  ***
-      !!          
+      !!
       !! ** Purpose : Apply the sea ice assimilation increment.
       !!
       !! ** Method  : Direct initialization or Incremental Analysis Updating.
       !!
-      !! ** Action  : 
+      !! ** Action  :
       !!
       !! History :
@@ -998,9 +1002,9 @@
 
             it = kt - nit000 + 1
-            zincwgt = wgtiau(it)      ! IAU weight for the current time step 
+            zincwgt = wgtiau(it)      ! IAU weight for the current time step
             ! note this is not a tendency so should not be divided by rdt (as with the tracer and other increments)
 
             IF(lwp) THEN
-               WRITE(numout,*) 
+               WRITE(numout,*)
                WRITE(numout,*) 'seaice_asm_inc : sea ice conc IAU at time step = ', &
                   &  kt,' with IAU weight = ', wgtiau(it)
@@ -1021,6 +1025,6 @@
             ! Nudge sea ice depth to bring it up to a required minimum depth
 
-            WHERE( zseaicendg(:,:) > 0.0_wp .AND. hicif(:,:) < zhicifmin ) 
-               zhicifinc(:,:) = (zhicifmin - hicif(:,:)) * zincwgt    
+            WHERE( zseaicendg(:,:) > 0.0_wp .AND. hicif(:,:) < zhicifmin )
+               zhicifinc(:,:) = (zhicifmin - hicif(:,:)) * zincwgt
             ELSEWHERE
                zhicifinc(:,:) = 0.0_wp
@@ -1029,5 +1033,5 @@
 ! nudge ice depth
             hicif(:,:)=hicif(:,:) + zhicifinc(:,:)
-            phicif(:,:)=phicif(:,:) + zhicifinc(:,:)       
+            phicif(:,:)=phicif(:,:) + zhicifinc(:,:)
 
 ! seaice salinity balancing (to add)
@@ -1071,9 +1075,9 @@
             zofrld(:,:)=frld(:,:)
             zohicif(:,:)=hicif(:,:)
- 
+
             ! Initialize the now fields the background + increment
 
             frld(:,:) = MIN( MAX( frld(:,:) - seaice_bkginc(:,:), 0.0_wp), 1.0_wp)
-            pfrld(:,:) = frld(:,:) 
+            pfrld(:,:) = frld(:,:)
             fr_i(:,:) = 1.0_wp - frld(:,:)        ! adjust ice fraction
 
@@ -1082,6 +1086,6 @@
             ! Nudge sea ice depth to bring it up to a required minimum depth
 
-            WHERE( zseaicendg(:,:) > 0.0_wp .AND. hicif(:,:) < zhicifmin ) 
-               zhicifinc(:,:) = (zhicifmin - hicif(:,:)) * zincwgt    
+            WHERE( zseaicendg(:,:) > 0.0_wp .AND. hicif(:,:) < zhicifmin )
+               zhicifinc(:,:) = (zhicifmin - hicif(:,:)) * zincwgt
             ELSEWHERE
                zhicifinc(:,:) = 0.0_wp
@@ -1090,10 +1094,10 @@
 ! nudge ice depth
             hicif(:,:)=hicif(:,:) + zhicifinc(:,:)
-            phicif(:,:)=phicif(:,:)       
+            phicif(:,:)=phicif(:,:)
 
 ! seaice salinity balancing (to add)
-  
-#endif
- 
+
+#endif
+
 #if defined key_cice
 
@@ -1110,14 +1114,14 @@
 #if defined key_cice
 
-! Zero ice increment tendency into CICE 
+! Zero ice increment tendency into CICE
             ndaice_da(:,:) = 0.0_wp
 
 #endif
-         
+
          ENDIF
 
 !#if defined key_lim2 || defined key_cice
 !
-!            IF (ln_seaicebal ) THEN       
+!            IF (ln_seaicebal ) THEN
 !             !! balancing salinity increments
 !             !! simple case from limflx.F90 (doesn't include a mass flux)
@@ -1131,5 +1135,5 @@
 !
 !             DO jj = 1, jpj
-!               DO ji = 1, jpi 
+!               DO ji = 1, jpi
 !           ! calculate change in ice and snow mass per unit area
 !           ! positive values imply adding salt to the ocean (results from ice formation)
@@ -1142,10 +1146,10 @@
 !
 !           ! prevent small mld
-!           ! less than 10m can cause salinity instability 
+!           ! less than 10m can cause salinity instability
 !                 IF (mld < 10) mld=10
 !
-!           ! set to bottom of a level 
+!           ! set to bottom of a level
 !                 DO jk = jpk-1, 2, -1
-!                   IF ((mld > gdepw(ji,jj,jk)) .and. (mld < gdepw(ji,jj,jk+1))) THEN 
+!                   IF ((mld > gdepw(ji,jj,jk)) .and. (mld < gdepw(ji,jj,jk+1))) THEN
 !                     mld=gdepw(ji,jj,jk+1)
 !                     jkmax=jk
@@ -1154,5 +1158,5 @@
 !
 !            ! avoid applying salinity balancing in shallow water or on land
-!            ! 
+!            !
 !
 !            ! dsal_ocn (psu kg m^-2) / (kg m^-3 * m)
@@ -1165,9 +1169,9 @@
 !
 !           ! put increments in for levels in the mixed layer
-!           ! but prevent salinity below a threshold value 
-!
-!                   DO jk = 1, jkmax              
-!
-!                     IF (dsal_ocn > 0.0_wp .or. sb(ji,jj,jk)+dsal_ocn > sal_thresh) THEN 
+!           ! but prevent salinity below a threshold value
+!
+!                   DO jk = 1, jkmax
+!
+!                     IF (dsal_ocn > 0.0_wp .or. sb(ji,jj,jk)+dsal_ocn > sal_thresh) THEN
 !                           sb(ji,jj,jk) = sb(ji,jj,jk) + dsal_ocn
 !                           sn(ji,jj,jk) = sn(ji,jj,jk) + dsal_ocn
@@ -1180,7 +1184,7 @@
 !      !
 !      !! Adjust fsalt. A +ve fsalt means adding salt to ocean
-!      !!           fsalt(ji,jj) =  fsalt(ji,jj) + zpmess     ! adjust fsalt  
-!      !!               
-!      !!           emps(ji,jj) = emps(ji,jj) + zpmess        ! or adjust emps (see icestp1d) 
+!      !!           fsalt(ji,jj) =  fsalt(ji,jj) + zpmess     ! adjust fsalt
+!      !!
+!      !!           emps(ji,jj) = emps(ji,jj) + zpmess        ! or adjust emps (see icestp1d)
 !      !!                                                     ! E-P (kg m-2 s-2)
 !      !            emp(ji,jj) = emp(ji,jj) + zpmess          ! E-P (kg m-2 s-2)
Index: branches/2012/dev_v3_4_STABLE_2012/NEMOGCM/NEMO/OPA_SRC/DOM/istate.F90
===================================================================
--- branches/2012/dev_v3_4_STABLE_2012/NEMOGCM/NEMO/OPA_SRC/DOM/istate.F90	(revision 5156)
+++ branches/2012/dev_v3_4_STABLE_2012/NEMOGCM/NEMO/OPA_SRC/DOM/istate.F90	(revision 5168)
@@ -10,8 +10,8 @@
    !!            8.0  !  2001-09  (M. Levy, M. Ben Jelloul)  istate_uvg
    !!   NEMO     1.0  !  2003-08  (G. Madec, C. Talandier)  F90: Free form, modules + EEL R5
-   !!             -   !  2004-05  (A. Koch-Larrouy)  istate_gyre 
+   !!             -   !  2004-05  (A. Koch-Larrouy)  istate_gyre
    !!            2.0  !  2006-07  (S. Masson)  distributed restart using iom
    !!            3.3  !  2010-10  (C. Ethe) merge TRC-TRA
-   !!            3.4  !  2011-04  (G. Madec) Merge of dtatem and dtasal & suppression of tb,tn/sb,sn 
+   !!            3.4  !  2011-04  (G. Madec) Merge of dtatem and dtasal & suppression of tb,tn/sb,sn
    !!----------------------------------------------------------------------
 
@@ -24,6 +24,6 @@
    !!   istate_uvg    : initial velocity in geostropic balance
    !!----------------------------------------------------------------------
-   USE oce             ! ocean dynamics and active tracers 
-   USE dom_oce         ! ocean space and time domain 
+   USE oce             ! ocean dynamics and active tracers
+   USE dom_oce         ! ocean space and time domain
    USE daymod          ! calendar
    USE eosbn2          ! eq. of state, Brunt Vaisala frequency (eos     routine)
@@ -46,4 +46,5 @@
    USE wrk_nemo        ! Memory allocation
    USE timing          ! Timing
+   USE asminc
 
    IMPLICIT NONE
@@ -65,5 +66,5 @@
       !!----------------------------------------------------------------------
       !!                   ***  ROUTINE istate_init  ***
-      !! 
+      !!
       !! ** Purpose :   Initialization of the dynamics and tracer fields.
       !!----------------------------------------------------------------------
@@ -83,6 +84,6 @@
       rhd  (:,:,:  ) = 0.e0
       rhop (:,:,:  ) = 0.e0
-      rn2  (:,:,:  ) = 0.e0 
-      tsa  (:,:,:,:) = 0.e0    
+      rn2  (:,:,:  ) = 0.e0
+      tsa  (:,:,:,:) = 0.e0
 
       IF( ln_rstart ) THEN                    ! Restart from a file
@@ -100,8 +101,8 @@
          CALL day_init                           ! model calendar (using both namelist and restart infos)
          !                                       ! Initialization of ocean to zero
-         !   before fields      !       now fields     
+         !   before fields      !       now fields
          sshb (:,:)   = 0._wp   ;   sshn (:,:)   = 0._wp
          ub   (:,:,:) = 0._wp   ;   un   (:,:,:) = 0._wp
-         vb   (:,:,:) = 0._wp   ;   vn   (:,:,:) = 0._wp  
+         vb   (:,:,:) = 0._wp   ;   vn   (:,:,:) = 0._wp
          rotb (:,:,:) = 0._wp   ;   rotn (:,:,:) = 0._wp
          hdivb(:,:,:) = 0._wp   ;   hdivn(:,:,:) = 0._wp
@@ -109,5 +110,5 @@
          IF( cp_cfg == 'eel' ) THEN
             CALL istate_eel                      ! EEL   configuration : start from pre-defined U,V T-S fields
-         ELSEIF( cp_cfg == 'gyre' ) THEN         
+         ELSEIF( cp_cfg == 'gyre' ) THEN
             CALL istate_gyre                     ! GYRE  configuration : start from pre-defined T-S fields
          ELSEIF( ln_tsd_init      ) THEN         ! Initial T-S fields read in files
@@ -124,5 +125,5 @@
             &                                       rhd, gru , grv  )   ! of t,s,rd at ocean bottom
 #endif
-         !   
+         !
          ! - ML - sshn could be modified by istate_eel, so that initialization of fse3t_b is done here
          IF( lk_vvl ) THEN
@@ -133,5 +134,5 @@
          !                                       ! define e3u_b, e3v_b from e3t_b initialized in domzgr
          CALL dom_vvl_2( nit000, fse3u_b(:,:,:), fse3v_b(:,:,:) )
-         ! 
+         !
       ENDIF
       !
@@ -146,4 +147,8 @@
       ENDIF
       !
+      IF ( lk_asminc .AND. ln_asmdin ) THEN
+         neuler = 0            ! restart with an euler from the corrected background
+      ENDIF
+
       IF( nn_timing == 1 )  CALL timing_stop('istate_init')
       !
@@ -153,6 +158,6 @@
       !!---------------------------------------------------------------------
       !!                  ***  ROUTINE istate_t_s  ***
-      !!   
-      !! ** Purpose :   Intialization of the temperature field with an 
+      !!
+      !! ** Purpose :   Intialization of the temperature field with an
       !!      analytical profile or a file (i.e. in EEL configuration)
       !!
@@ -184,5 +189,5 @@
       !!----------------------------------------------------------------------
       !!                   ***  ROUTINE istate_eel  ***
-      !! 
+      !!
       !! ** Purpose :   Initialization of the dynamics and tracers for EEL R5
       !!      configuration (channel with or without a topographic bump)
@@ -195,5 +200,5 @@
       USE divcur     ! hor. divergence & rel. vorticity      (div_cur routine)
       USE iom
- 
+
       INTEGER  ::   inum              ! temporary logical unit
       INTEGER  ::   ji, jj, jk        ! dummy loop indices
@@ -207,5 +212,5 @@
       !!----------------------------------------------------------------------
 
-      SELECT CASE ( jp_cfg ) 
+      SELECT CASE ( jp_cfg )
          !                                              ! ====================
          CASE ( 5 )                                     ! EEL R5 configuration
@@ -244,8 +249,8 @@
             ! set the dynamics: U,V, hdiv, rot (and ssh if necessary)
             ! ----------------
-            ! Start EEL5 configuration with barotropic geostrophic velocities 
+            ! Start EEL5 configuration with barotropic geostrophic velocities
             ! according the sshb and sshn SSH imposed.
             ! we assume a uniform grid (hence the use of e1t(1,1) for delta_y)
-            ! we use the Coriolis frequency at mid-channel.   
+            ! we use the Coriolis frequency at mid-channel.
             ub(:,:,:) = zueel * umask(:,:,:)
             un(:,:,:) = ub(:,:,:)
@@ -254,5 +259,5 @@
             !
             DO jj = 1, jpjglo
-               zssh(:,jj) = - (FLOAT(jj)- FLOAT(jpjglo-1)/2.)*zueel*e1t(1,1)*zfcor/grav 
+               zssh(:,jj) = - (FLOAT(jj)- FLOAT(jpjglo-1)/2.)*zueel*e1t(1,1)*zfcor/grav
             END DO
             !
@@ -273,6 +278,6 @@
             sshn(:,:) = sshb(:,:)                   ! set now ssh to the before value
             !
-            IF( nn_rstssh /= 0 ) THEN  
-               nn_rstssh = 0                        ! hand-made initilization of ssh 
+            IF( nn_rstssh /= 0 ) THEN
+               nn_rstssh = 0                        ! hand-made initilization of ssh
                CALL ctl_warn( 'istate_eel: force nn_rstssh = 0' )
             ENDIF
@@ -326,5 +331,5 @@
       !!----------------------------------------------------------------------
       !!                   ***  ROUTINE istate_gyre  ***
-      !! 
+      !!
       !! ** Purpose :   Initialization of the dynamics and tracers for GYRE
       !!      configuration (double gyre with rotated domain)
@@ -351,6 +356,6 @@
                        &           * (-TANH( (500-fsdept(ji,jj,jk)) / 150 ) + 1) / 2               &
                        &       + (      15. * ( 1. - TANH( (fsdept(ji,jj,jk)-50.) / 1500.) )       &
-                       &                - 1.4 * TANH((fsdept(ji,jj,jk)-100.) / 100.)               &    
-                       &                + 7.  * (1500. - fsdept(ji,jj,jk)) / 1500.             )   & 
+                       &                - 1.4 * TANH((fsdept(ji,jj,jk)-100.) / 100.)               &
+                       &                + 7.  * (1500. - fsdept(ji,jj,jk)) / 1500.             )   &
                        &           * (-TANH( (fsdept(ji,jj,jk) - 500) / 150) + 1) / 2
                   tsn(ji,jj,jk,jp_tem) = tsn(ji,jj,jk,jp_tem) * tmask(ji,jj,jk)
@@ -363,5 +368,5 @@
                      &                + 0.2  * TANH( (fsdept(ji,jj,jk) - 35.  ) / 100. )       &
                      &                + 0.2  * TANH( (fsdept(ji,jj,jk) - 1000.) / 5000.)    )  &
-                     &              * (-TANH((fsdept(ji,jj,jk) - 500) / 150) + 1) / 2 
+                     &              * (-TANH((fsdept(ji,jj,jk) - 500) / 150) + 1) / 2
                   tsn(ji,jj,jk,jp_sal) = tsn(ji,jj,jk,jp_sal) * tmask(ji,jj,jk)
                   tsb(ji,jj,jk,jp_sal) = tsn(ji,jj,jk,jp_sal)
@@ -379,8 +384,8 @@
          ! ----------------------
          CALL iom_open ( 'data_tem', inum )
-         CALL iom_get ( inum, jpdom_data, 'votemper', tsn(:,:,:,jp_tem) ) 
+         CALL iom_get ( inum, jpdom_data, 'votemper', tsn(:,:,:,jp_tem) )
          CALL iom_close( inum )
 
-         tsn(:,:,:,jp_tem) = tsn(:,:,:,jp_tem) * tmask(:,:,:) 
+         tsn(:,:,:,jp_tem) = tsn(:,:,:,jp_tem) * tmask(:,:,:)
          tsb(:,:,:,jp_tem) = tsn(:,:,:,jp_tem)
 
@@ -388,8 +393,8 @@
          ! -------------------
          CALL iom_open ( 'data_sal', inum )
-         CALL iom_get ( inum, jpdom_data, 'vosaline', tsn(:,:,:,jp_sal) ) 
+         CALL iom_get ( inum, jpdom_data, 'vosaline', tsn(:,:,:,jp_sal) )
          CALL iom_close( inum )
 
-         tsn(:,:,:,jp_sal) = tsn(:,:,:,jp_sal) * tmask(:,:,:) 
+         tsn(:,:,:,jp_sal) = tsn(:,:,:,jp_sal) * tmask(:,:,:)
          tsb(:,:,:,jp_sal) = tsn(:,:,:,jp_sal)
 
@@ -412,5 +417,5 @@
       !! ** Purpose :   Compute the geostrophic velocities from (tn,sn) fields
       !!
-      !! ** Method  :   Using the hydrostatic hypothesis the now hydrostatic 
+      !! ** Method  :   Using the hydrostatic hypothesis the now hydrostatic
       !!      pressure is computed by integrating the in-situ density from the
       !!      surface to the bottom.
@@ -429,5 +434,5 @@
       CALL wrk_alloc( jpi, jpj, jpk, zprn)
       !
-      IF(lwp) WRITE(numout,*) 
+      IF(lwp) WRITE(numout,*)
       IF(lwp) WRITE(numout,*) 'istate_uvg : Start from Geostrophy'
       IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
@@ -437,5 +442,5 @@
 
       zalfg = 0.5 * grav * rau0
-      
+
       zprn(:,:,1) = zalfg * fse3w(:,:,1) * ( 1 + rhd(:,:,1) )       ! Surface value
 
@@ -443,5 +448,5 @@
          zprn(:,:,jk) = zprn(:,:,jk-1)   &
             &         + zalfg * fse3w(:,:,jk) * ( 2. + rhd(:,:,jk) + rhd(:,:,jk-1) )
-      END DO  
+      END DO
 
       ! Compute geostrophic balance
@@ -485,8 +490,8 @@
       CALL lbc_lnk( un, 'U', -1. )
       CALL lbc_lnk( vn, 'V', -1. )
-      
+
       ub(:,:,:) = un(:,:,:)
       vb(:,:,:) = vn(:,:,:)
-      
+
       ! WARNING !!!!!
       ! after initializing u and v, we need to calculate the initial streamfunction bsf.
@@ -516,5 +521,5 @@
       un(:,:,:) = ub(:,:,:)
       vn(:,:,:) = vb(:,:,:)
-       
+
       ! Compute the divergence and curl
 
Index: branches/2012/dev_v3_4_STABLE_2012/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90
===================================================================
--- branches/2012/dev_v3_4_STABLE_2012/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90	(revision 5156)
+++ branches/2012/dev_v3_4_STABLE_2012/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90	(revision 5168)
@@ -77,4 +77,5 @@
 #endif
    USE tamtrj          ! Output trajectory, needed for TAM
+   USE asminc
 
    IMPLICIT NONE
@@ -110,5 +111,5 @@
       !!              Madec, 2008, internal report, IPSL.
       !!----------------------------------------------------------------------
-      INTEGER ::   istp       ! time step index
+      INTEGER ::   istp, jk      ! time step index
       !!----------------------------------------------------------------------
       !
@@ -155,4 +156,17 @@
             CALL Agrif_Step( stp )           ! AGRIF: time stepping
 #else
+            !-------------------------------------------------------------!
+            ! This trick ensures a minimum consistency between mixing coef.
+            ! and other variables
+            IF( lk_zdftke .AND. lk_asminc .AND. ln_asmdin .AND. ( istp == nit000) ) THEN
+            CALL tke_avn                        ! recompute avt, avm,
+                                                ! avmu, avmv and dissl (approximation)
+               DO jk = nit000, nit000 + 2
+                  CALL zdf_tke( jk )   ;
+               END DO
+
+            END IF
+            !-------------------------------------------------------------!
+
             CALL stp( istp )                 ! standard time stepping
 #endif
@@ -322,5 +336,6 @@
                             CALL dyn_nept_init  ! simplified form of Neptune effect
                             CALL flush(numout)
-
+                            CALL asm_inc_rea_nam ! to fix properly neuler value
+                                                 ! in istate_init
                             CALL  istate_init   ! ocean initial state (Dynamics and tracers)
 
