Index: branches/dev_1784_ASM/NEMO/OPA_SRC/ASM/asminc.F90
===================================================================
--- branches/dev_1784_ASM/NEMO/OPA_SRC/ASM/asminc.F90	(revision 2002)
+++ branches/dev_1784_ASM/NEMO/OPA_SRC/ASM/asminc.F90	(revision 2002)
@@ -0,0 +1,976 @@
+MODULE asminc
+   !!======================================================================
+   !!                       ***  MODULE asminc  ***
+   !! Assimilation increment : Apply an increment generated by data
+   !!                          assimilation
+   !!======================================================================
+
+   !!----------------------------------------------------------------------
+   !!   'key_asminc' : Switch on the assimilation increment interface
+   !!----------------------------------------------------------------------
+   !!   asm_inc_init : Initialize the increment arrays and IAU weights
+   !!   calc_date    : Compute the calendar date YYYYMMDD on a given step
+   !!   tra_asm_inc  : Apply the tracer (T and S) increments
+   !!   dyn_asm_inc  : Apply the dynamic (u and v) increments
+   !!   ssh_asm_inc  : Apply the SSH increment
+   !!----------------------------------------------------------------------
+   !! * Modules used   
+   USE par_kind, ONLY : &       ! Precision variables
+      & wp
+   USE in_out_manager, ONLY : & ! I/O manager 
+      & lwp,      &
+      & numnam,   &
+      & numout,   &
+      & ctl_warn, &
+      & ctl_stop, &
+      & nit000,   &
+      & nstop,    &
+      & ln_rstart
+   USE par_oce, ONLY : & ! Ocean space and time domain variables
+      & jpi,  & 
+      & jpj,  &
+      & jpk,  &
+      & jpkm1
+   USE dom_oce, ONLY : & ! Ocean space and time domain
+      & rdt,    &
+      & n_cla,  &
+      & neuler, &
+      & ln_zps, &
+      & tmask,  &
+      & umask,  &
+      & vmask
+   USE c1d, ONLY : &    ! 1D initialization
+      & lk_c1d
+   USE oce, ONLY : &      ! Dynamics and active tracers defined in memory
+      & ub, un, ua,    &
+      & vb, vn, va,    &
+      & tb, tn, ta,    &
+      & sb, sn, sa,    &
+      & sshb, sshn,    &
+      & rhd, rhop,     &
+      & rotb, rotn,    &
+      & hdivb, hdivn,  &
+      & gtu, gsu, gru, &
+      & gtv, gsv, grv 
+   USE divcur, ONLY : &   ! Horizontal divergence and relative vorticity
+      & div_cur
+   USE cla_div, ONLY : &  ! Specific update of the horizontal divergence
+      & div_cla           ! (specific to ORCA_R2)
+   USE eosbn2, ONLY : &   ! Equation of state - in situ and potential density
+      & eos
+   USE zpshde, ONLY : &   ! Partial step : Horizontal Derivative
+      & zps_hde
+!   USE phycst, ONLY : &   ! Calendar parameters
+!      & rjjss  
+   USE iom                ! Library to read input files
+   USE asmpar             ! Parameters for the assmilation interface
+   USE dom_oce, ONLY : &
+      & ndastp
+!   USE daymod, ONLY : &
+!      & nmonth_len         ! length of month in days
+!      & nbiss, &
+!      & nobis, &
+ 
+
+   IMPLICIT NONE
+
+   !! * Routine accessibility
+   PRIVATE
+   PUBLIC asm_inc_init, & !: Initialize the increment arrays and IAU weights
+      &   calc_date,    & !: Compute the calendar date YYYYMMDD on a given step
+      &   tra_asm_inc,  & !: Apply the tracer (T and S) increments
+      &   dyn_asm_inc,  & !: Apply the dynamic (u and v) increments
+      &   ssh_asm_inc     !: Apply the SSH increment
+
+   !! * Private Module variables
+   REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: &
+      & t_bkg,    &       !: Background temperature 
+      & s_bkg,    &       !: Background salinity
+      & u_bkg,    &       !: Background u-component velocity
+      & v_bkg,    &       !: Background v-component velocity
+      & t_bkginc, &       !: Increment to the background temperature
+      & s_bkginc, &       !: Increment to the background salinity
+      & u_bkginc, &       !: Increment to the u-component velocity 
+      & v_bkginc          !: Increment to the v-component velocity 
+         
+   REAL(wp), PRIVATE, DIMENSION(:,:), ALLOCATABLE :: &
+      & ssh_bkg,   &      !: Background sea surface height
+      & ssh_bkginc        !: Increment to the background sea surface height
+
+   REAL(wp), PUBLIC, DIMENSION(:), ALLOCATABLE :: &
+      & wgtiau            !: IAU weights for each time step
+
+   !! * Shared Module variables
+   LOGICAL, PUBLIC, PARAMETER :: &
+#if defined key_asminc
+      & lk_asminc = .TRUE.   !: Logical switch for assimilation increment interface
+#else
+      & lk_asminc = .FALSE.  !: No assimilation increments
+#endif
+
+   LOGICAL, PUBLIC :: & 
+      & ln_bkgwri = .FALSE., & !: No output of the background state fields
+      & ln_trjwri = .FALSE., & !: No output of the state trajectory fields
+      & ln_asmiau = .FALSE., & !: No applying forcing with an assimilation increment
+      & ln_asmdin = .FALSE., & !: No direct initialization
+      & ln_trainc = .FALSE., & !: No tracer (T and S) assimilation increments
+      & ln_dyninc = .FALSE., & !: No dynamics (u and v) assimilation increments
+      & ln_sshinc = .FALSE., & !: No sea surface height assimilation increment
+      & ln_salfix = .FALSE.    !: Apply minimum salinity check
+
+   REAL, PUBLIC :: &
+      & salfixmin              !: Ensure that the salinity is larger than 
+   !                           !: this  value if (ln_salfix)
+#if defined key_asminc
+   REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: &
+      & ssh_iau           !: IAU-weighted sea surface height increment
+#endif
+   INTEGER, PUBLIC :: &
+      & nitbkg,    &   !: Time step of the background state used in the Jb term
+                       !: (relative to the cycle interval [0,nitend-nit000-1])
+      & nitdin,    &   !: Time step of the background state for direct initialization
+                       !: (relative to the cycle interval [0,nitend-nit000-1])
+      & nitiaustr, &   !: Time step of the start of the IAU interval 
+                       !: (relative to the cycle interval [0,nitend-nit000-1])
+                       !: 
+      & nitiaufin, &   !: Time step of the end of the IAU interval 
+                       !: (relative to the cycle interval [0,nitend-nit000-1])
+      & niaufn         !: Type of IAU weighing function
+                       !:   0 = Constant weighting
+                       !:   1 = Linear hat-like, centred in middle of IAU interval 
+
+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  : 
+      !!
+      !! History :
+      !!        !  07-03  (M. Martin) Met Office version
+      !!        !  07-04  (A. Weaver) Merge with OPAVAR/NEMOVAR
+      !!----------------------------------------------------------------------
+
+      IMPLICIT NONE
+
+      !! * Modules used
+      NAMELIST/nam_asminc/ ln_bkgwri, ln_trjwri, &
+         &                 ln_trainc, ln_dyninc, ln_sshinc, &
+         &                 ln_asmdin, ln_asmiau, &
+         &                 nitbkg, nitdin, nitiaustr, nitiaufin, niaufn, &
+         &                 nittrjfrq, ln_salfix, salfixmin
+
+      !! * Local declarations
+      INTEGER :: &
+         & jt
+      INTEGER :: &
+         & imid,           &
+         & inum,           &
+         & iiauper,        & ! Number of time steps in the IAU period
+         & icycper,        & ! Number of time steps in the cycle
+         & iitend_date,    & ! Date YYYYMMDD of final time step
+         & iitbkg_date,    & ! Date YYYYMMDD of background time step for Jb term
+         & iitdin_date,    & ! Date YYYYMMDD of background time step for DI
+         & iitiaustr_date, & ! Date YYYYMMDD of IAU interval start time step
+         & iitiaufin_date    ! Date YYYYMMDD of IAU interval final time step
+      REAL(wp) :: &
+         & znorm,       & ! Normalization factor for IAU weights
+         & ztotwgt,     & ! Value of time-integrated IAU weights 
+                          ! (should be equal to one)
+         & z_inc_dateb, & ! Start date of interval on which increment is valid
+         & z_inc_datef, & ! End date of interval on which increment is valid
+         & zdate_bkg      ! Date in background state file for DI
+      REAL(wp) :: &
+         & zdate_inc      ! Time axis in increments file
+
+      !-----------------------------------------------------------------------
+      ! Read Namelist nam_asminc : assimilation increment interface
+      !-----------------------------------------------------------------------
+
+      ! Set default values
+      ln_bkgwri = .FALSE.
+      ln_trjwri = .FALSE.
+      ln_trainc = .FALSE.
+      ln_dyninc = .FALSE.
+      ln_sshinc = .FALSE.
+      ln_asmdin = .FALSE.
+      ln_asmiau = .TRUE.
+      ln_salfix = .FALSE.
+      salfixmin = -9999
+      nitbkg    = 0
+      nitdin    = 0      
+      nitiaustr = 1
+      nitiaufin = 150      ! = 10 days with ORCA2
+      niaufn    = 0
+      nittrjfrq = 1
+
+      REWIND ( numnam )
+      READ   ( numnam, nam_asminc )
+
+      ! Control print
+      IF(lwp) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) 'asm_inc_init : Assimilation increment initialization :'
+         WRITE(numout,*) '~~~~~~~~~~~~'
+         WRITE(numout,*) '          Namelist namasm : set assimilation increment parameters'
+         WRITE(numout,*) '             Logical switch for writing out background state         ', &
+            &            ' ln_bkgwri = ', ln_bkgwri
+         WRITE(numout,*) '             Logical switch for writing out state trajectory         ', &
+            &            ' ln_trjwri = ', ln_trjwri
+         WRITE(numout,*) '             Logical switch for applying tracer increments           ', &
+            &            ' ln_trainc = ', ln_trainc
+         WRITE(numout,*) '             Logical switch for applying velocity increments         ', &
+            &            ' ln_dyninc = ', ln_dyninc
+         WRITE(numout,*) '             Logical switch for applying SSH increments              ', &
+            &            ' ln_sshinc = ', ln_sshinc
+         WRITE(numout,*) '             Logical switch for Direct Initialization (DI)           ', &
+            &            ' ln_asmdin = ', ln_asmdin
+         WRITE(numout,*) '             Logical switch for Incremental Analysis Updating (IAU)  ', &
+            &            ' ln_asmiau = ', ln_asmiau
+         WRITE(numout,*) '             Timestep of background in [0,nitend-nit000-1]           ', &
+            &            '    nitbkg = ', nitbkg
+         WRITE(numout,*) '             Timestep of background for DI in [0,nitend-nit000-1]    ', &
+            &            '    nitdin = ', nitdin
+         WRITE(numout,*) '             Timestep of start of IAU interval in [0,nitend-nit000-1]', &
+            &            ' nitiaustr = ', nitiaustr
+         WRITE(numout,*) '             Timestep of end of IAU interval in [0,nitend-nit000-1]  ', &
+            &            ' nitiaufin = ', nitiaufin
+         WRITE(numout,*) '             Type of IAU weighting function                          ', &
+            &            '    niaufn = ', niaufn
+         WRITE(numout,*) '             Frequency of trajectory output for 4D-VAR               ', &
+            &            ' nittrjfrq = ', nittrjfrq
+         WRITE(numout,*) '             Logical switch for ensuring that the sa > salfixmin     ', &
+            &            ' ln_salinc = ', ln_salfix
+         WRITE(numout,*) '             Minimum salinity after applying the increments          ', &
+            &            ' salfixmin = ', salfixmin
+      ENDIF
+
+      nitbkg_r    = nitbkg    + nit000 - 1  ! Background time referenced to nit000
+      nitdin_r    = nitdin    + nit000 - 1  ! Background time for DI referenced to nit000
+      nitiaustr_r = nitiaustr + nit000 - 1  ! Start of IAU interval referenced to nit000
+      nitiaufin_r = nitiaufin + nit000 - 1  ! End of IAU interval referenced to nit000
+
+      iiauper = nitiaufin_r - nitiaustr_r + 1  ! IAU interval length
+      icycper = nitend      - nit000      + 1  ! Cycle interval length
+
+      ! Date of final time step
+      CALL calc_date( nit000, nitend, ndate0, iitend_date )
+
+      ! Background time for Jb referenced to ndate0
+      CALL calc_date( nit000, nitbkg_r, ndate0, iitbkg_date )
+
+      ! Background time for DI referenced to ndate0
+      CALL calc_date( nit000, nitdin_r, ndate0, iitdin_date )
+
+      ! IAU start time referenced to ndate0
+      CALL calc_date( nit000, nitiaustr_r, ndate0, iitiaustr_date )
+
+      ! IAU end time referenced to ndate0
+      CALL calc_date( nit000, nitiaufin_r, ndate0, iitiaufin_date )
+
+      IF(lwp) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) '   Time steps referenced to current cycle:'
+         WRITE(numout,*) '       iitrst      = ', nit000 - 1
+         WRITE(numout,*) '       nit000      = ', nit000
+         WRITE(numout,*) '       nitend      = ', nitend
+         WRITE(numout,*) '       nitbkg_r    = ', nitbkg_r
+         WRITE(numout,*) '       nitdin_r    = ', nitdin_r
+         WRITE(numout,*) '       nitiaustr_r = ', nitiaustr_r
+         WRITE(numout,*) '       nitiaufin_r = ', nitiaufin_r
+         WRITE(numout,*) '       nittrjfrq   = ', nittrjfrq
+         WRITE(numout,*)
+         WRITE(numout,*) '   Dates referenced to current cycle:'
+         WRITE(numout,*) '       ndastp         = ', ndastp
+         WRITE(numout,*) '       ndate0         = ', ndate0
+         WRITE(numout,*) '       iitend_date    = ', iitend_date
+         WRITE(numout,*) '       iitbkg_date    = ', iitbkg_date
+         WRITE(numout,*) '       iitdin_date    = ', iitdin_date
+         WRITE(numout,*) '       iitiaustr_date = ', iitiaustr_date
+         WRITE(numout,*) '       iitiaufin_date = ', iitiaufin_date
+      ENDIF
+
+      IF ( nacc /= 0 ) &
+         & CALL ctl_stop( ' nacc /= 0 and key_asminc :',  &
+         &                ' Assimilation increments have only been implemented', &
+         &                ' for synchronous time stepping' )
+
+      IF ( ( ln_asmdin ).AND.( ln_asmiau ) )   &
+         & CALL ctl_stop( ' ln_asmdin and ln_asmiau :', &
+         &                ' Choose Direct Initialization OR Incremental Analysis Updating')
+
+      IF (      ( ( .NOT. ln_asmdin ).AND.( .NOT. ln_asmiau ) ) &
+           .AND.( ( ln_trainc ).OR.( ln_dyninc ).OR.( ln_sshinc ) ) ) &
+         & CALL ctl_stop( ' One or more of ln_trainc, ln_dyninc and ln_sshinc is set to .true.', &
+         &                ' but ln_asmdin and ln_asmiau are both set to .false. :', &
+         &                ' Inconsistent options')
+
+      IF ( ( ln_bkgwri ).AND.( ( ln_asmdin ).OR.( ln_asmiau ) ) )  &
+         & CALL ctl_stop( ' ln_bkgwri and either ln_asmdin or ln_asmiau are set to .true.:', &
+         &                ' The background state must be written before applying the increments')
+
+      IF ( ( niaufn /= 0 ).AND.( niaufn /= 1 ) ) &
+         & CALL ctl_stop( ' niaufn /= 0 or niaufn /=1 :',  &
+         &                ' Type IAU weighting function is invalid')
+
+      IF ( ( .NOT. ln_trainc ).AND.( .NOT. ln_dyninc ).AND.( .NOT. ln_sshinc ) &
+         &                     )  &
+         & CALL ctl_warn( ' ln_trainc, ln_dyninc and ln_sshinc are set to .false. :', &
+         &                ' The assimilation increments are not applied')
+
+      IF ( ( ln_asmiau ).AND.( nitiaustr == nitiaufin ) ) &
+         & CALL ctl_stop( ' nitiaustr = nitiaufin :',  &
+         &                ' IAU interval is of zero length')
+
+      IF ( ( ln_asmiau ).AND.( ( nitiaustr_r < nit000 ).OR.( nitiaufin_r > nitend ) ) ) &
+         & CALL ctl_stop( ' nitiaustr or nitiaufin :',  &
+         &                ' IAU starting or final time step is outside the cycle interval', &
+         &                 ' Valid range nit000 to nitend')
+
+      IF ( ( nitbkg_r < nit000 - 1 ).OR.( nitbkg_r > nitend ) ) &
+         & CALL ctl_stop( ' nitbkg :',  &
+         &                ' Background time step is outside the cycle interval')
+
+      IF ( ( nitdin_r < nit000 - 1 ).OR.( nitdin_r > nitend ) ) &
+         & CALL ctl_stop( ' nitdin :',  &
+         &                ' Background time step for Direct Initialization is outside', &
+         &                ' the cycle interval')
+
+      IF ( nstop > 0 ) RETURN       ! if there are any errors then go no further
+
+      !--------------------------------------------------------------------
+      ! Initialize the Incremental Analysis Updating weighting function
+      !--------------------------------------------------------------------
+
+      IF ( ln_asmiau ) THEN
+
+         ALLOCATE( wgtiau( icycper ) )
+
+         wgtiau(:) = 0.0
+
+         IF ( niaufn == 0 ) THEN
+
+            !---------------------------------------------------------
+            ! Constant IAU forcing 
+            !---------------------------------------------------------
+
+            DO jt = 1, iiauper
+               wgtiau(jt+nitiaustr-1) = 1.0 / REAL( iiauper )
+            END DO
+
+         ELSEIF ( niaufn == 1 ) THEN
+
+            !---------------------------------------------------------
+            ! Linear hat-like, centred in middle of IAU interval 
+            !---------------------------------------------------------
+
+            ! Compute the normalization factor
+            znorm = 0.0
+            IF ( MOD( iiauper, 2 ) == 0 ) THEN  ! Even number of time steps in IAU interval
+               imid = iiauper / 2 
+               DO jt = 1, imid
+                  znorm = znorm + REAL( jt )
+               END DO
+               znorm = 2.0 * znorm
+            ELSE                               ! Odd number of time steps in IAU interval
+               imid = ( iiauper + 1 ) / 2        
+               DO jt = 1, imid - 1
+                  znorm = znorm + REAL( jt )
+               END DO
+               znorm = 2.0 * znorm + REAL( imid )
+            ENDIF
+            znorm = 1.0 / znorm
+
+            DO jt = 1, imid - 1
+               wgtiau(jt+nitiaustr-1) = REAL( jt ) * znorm
+            END DO
+            DO jt = imid, iiauper
+               wgtiau(jt+nitiaustr-1) = REAL( iiauper - jt + 1 ) * znorm
+            END DO
+
+         ENDIF
+
+         ! Test that the integral of the weights over the weighting interval equals 1
+          IF(lwp) THEN
+             WRITE(numout,*)
+             WRITE(numout,*) 'asm_inc_init : IAU weights'
+             WRITE(numout,*) '~~~~~~~~~~~~'
+             WRITE(numout,*) '             time step         IAU  weight'
+             WRITE(numout,*) '             =========     ====================='
+             ztotwgt = 0.0
+             DO jt = 1, icycper
+                ztotwgt = ztotwgt + wgtiau(jt)
+                WRITE(numout,*) '         ', jt, '       ', wgtiau(jt) 
+             END DO   
+             WRITE(numout,*) '         ==================================='
+             WRITE(numout,*) '         Time-integrated weight = ', ztotwgt
+             WRITE(numout,*) '         ==================================='
+          ENDIF
+         
+      ENDIF
+
+      !--------------------------------------------------------------------
+      ! Allocate and initialize the increment arrays
+      !--------------------------------------------------------------------
+
+      ALLOCATE( t_bkginc(jpi,jpj,jpk) )
+      ALLOCATE( s_bkginc(jpi,jpj,jpk) )
+      ALLOCATE( u_bkginc(jpi,jpj,jpk) )
+      ALLOCATE( v_bkginc(jpi,jpj,jpk) )
+      ALLOCATE( ssh_bkginc(jpi,jpj)   )
+#if defined key_asminc
+      ALLOCATE( ssh_iau(jpi,jpj)      )
+#endif
+      t_bkginc(:,:,:) = 0.0
+      s_bkginc(:,:,:) = 0.0
+      u_bkginc(:,:,:) = 0.0
+      v_bkginc(:,:,:) = 0.0
+      ssh_bkginc(:,:) = 0.0
+#if defined key_asminc
+      ssh_iau(:,:)    = 0.0
+#endif
+      IF ( ( ln_trainc ).OR.( ln_dyninc ).OR.( ln_sshinc ) ) THEN
+
+         !--------------------------------------------------------------------
+         ! Read the increments from file
+         !--------------------------------------------------------------------
+
+         CALL iom_open( c_asminc, inum )
+
+         CALL iom_get( inum, 'time', zdate_inc ) 
+
+         CALL iom_get( inum, 'z_inc_dateb', z_inc_dateb )
+         CALL iom_get( inum, 'z_inc_datef', z_inc_datef )
+         z_inc_dateb = zdate_inc
+         z_inc_datef = zdate_inc
+
+         IF(lwp) THEN
+            WRITE(numout,*) 
+            WRITE(numout,*) 'asm_inc_init : Assimilation increments valid ', &
+               &            ' between dates ', NINT( z_inc_dateb ),' and ',  &
+               &            NINT( z_inc_datef )
+            WRITE(numout,*) '~~~~~~~~~~~~'
+         ENDIF
+
+         IF (     ( NINT( z_inc_dateb ) < ndastp      ) &
+            & .OR.( NINT( z_inc_datef ) > iitend_date ) ) &
+            & CALL ctl_warn( ' Validity time of assimilation increments is ', &
+            &                ' outside the assimilation interval' )
+
+         IF ( ( ln_asmdin ).AND.( NINT( zdate_inc ) /= iitdin_date ) ) &
+            & CALL ctl_warn( ' Validity time of assimilation increments does ', &
+            &                ' not agree with Direct Initialization time' )
+
+         IF ( ln_trainc ) THEN   
+            CALL iom_get( inum, jpdom_autoglo, 'bckint', t_bkginc, 1 )
+            CALL iom_get( inum, jpdom_autoglo, 'bckins', s_bkginc, 1 )
+            ! Apply the masks
+            t_bkginc(:,:,:) = t_bkginc(:,:,:) * tmask(:,:,:)
+            s_bkginc(:,:,:) = s_bkginc(:,:,:) * tmask(:,:,:)
+            ! Set missing increments to 0.0 rather than 1e+20
+            ! to allow for differences in masks
+            WHERE( ABS( t_bkginc(:,:,:) ) > 1.0e+10 ) t_bkginc(:,:,:) = 0.0
+            WHERE( ABS( s_bkginc(:,:,:) ) > 1.0e+10 ) s_bkginc(:,:,:) = 0.0
+         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 )              
+            ! Apply the masks
+            u_bkginc(:,:,:) = u_bkginc(:,:,:) * umask(:,:,:)
+            v_bkginc(:,:,:) = v_bkginc(:,:,:) * vmask(:,:,:)
+            ! Set missing increments to 0.0 rather than 1e+20
+            ! to allow for differences in masks
+            WHERE( ABS( u_bkginc(:,:,:) ) > 1.0e+10 ) u_bkginc(:,:,:) = 0.0
+            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 )
+            ! Apply the masks
+            ssh_bkginc(:,:) = ssh_bkginc(:,:) * tmask(:,:,1)
+            ! Set missing increments to 0.0 rather than 1e+20
+            ! to allow for differences in masks
+            WHERE( ABS( ssh_bkginc(:,:) ) > 1.0e+10 ) ssh_bkginc(:,:) = 0.0
+         ENDIF
+
+         CALL iom_close( inum )
+ 
+      ENDIF
+
+      !-----------------------------------------------------------------------
+      ! Allocate and initialize the background state arrays
+      !-----------------------------------------------------------------------
+
+      IF ( ln_asmdin ) THEN
+
+         ALLOCATE( t_bkg(jpi,jpj,jpk) )
+         ALLOCATE( s_bkg(jpi,jpj,jpk) )
+         ALLOCATE( u_bkg(jpi,jpj,jpk) )
+         ALLOCATE( v_bkg(jpi,jpj,jpk) )
+         ALLOCATE( ssh_bkg(jpi,jpj)   )
+
+         t_bkg(:,:,:) = 0.0
+         s_bkg(:,:,:) = 0.0
+         u_bkg(:,:,:) = 0.0
+         v_bkg(:,:,:) = 0.0
+         ssh_bkg(:,:) = 0.0
+
+         !--------------------------------------------------------------------
+         ! Read from file the background state at analysis time
+         !--------------------------------------------------------------------
+
+         CALL iom_open( c_asmdin, inum )
+
+         CALL iom_get( inum, 'zdate', zdate_bkg ) 
+        
+         IF(lwp) THEN
+            WRITE(numout,*) 
+            WRITE(numout,*) 'asm_inc_init : Assimilation background state valid at : ', &
+               &  NINT( zdate_bkg )
+            WRITE(numout,*) '~~~~~~~~~~~~'
+         ENDIF
+
+         IF ( NINT( zdate_bkg ) /= iitdin_date ) &
+            & CALL ctl_warn( ' Validity time of assimilation background state does', &
+            &                ' not agree with Direct Initialization time' )
+
+         IF ( ln_trainc ) THEN   
+            CALL iom_get( inum, jpdom_autoglo, 'tn', t_bkg )
+            CALL iom_get( inum, jpdom_autoglo, 'sn', s_bkg )
+            t_bkg(:,:,:) = t_bkg(:,:,:) * tmask(:,:,:)
+            s_bkg(:,:,:) = s_bkg(:,:,:) * tmask(:,:,:)
+         ENDIF
+
+         IF ( ln_dyninc ) THEN   
+            CALL iom_get( inum, jpdom_autoglo, 'un', u_bkg )
+            CALL iom_get( inum, jpdom_autoglo, 'vn', v_bkg )
+            u_bkg(:,:,:) = u_bkg(:,:,:) * umask(:,:,:)
+            v_bkg(:,:,:) = v_bkg(:,:,:) * vmask(:,:,:)
+         ENDIF
+        
+         IF ( ln_sshinc ) THEN
+            CALL iom_get( inum, jpdom_autoglo, 'sshn', ssh_bkg )
+            ssh_bkg(:,:) = ssh_bkg(:,:) * tmask(:,:,1)
+         ENDIF
+
+         CALL iom_close( inum )
+
+      ENDIF
+
+   END SUBROUTINE asm_inc_init
+
+   SUBROUTINE calc_date( kit000, kt, kdate0, kdate )
+      !!----------------------------------------------------------------------
+      !!                    ***  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  : 
+      !!
+      !! History :
+      !!        !  07-04  (A. Weaver) 
+      !!        !  10-05  (D. Lea)        Update to work with NEMO vn3.2
+      !!----------------------------------------------------------------------
+
+      IMPLICIT NONE
+
+      !! * Arguments
+      INTEGER, INTENT(IN) :: &
+         & kit000, &  ! Initial time step
+         & kt,     &  ! Current time step referenced to kit000
+         & kdate0     ! Initial date
+      INTEGER, INTENT(OUT) :: &
+         & kdate      ! Current date reference to kdate0
+
+      !! * Local declarations
+      INTEGER :: &
+         & iyea0,   & ! Initial year
+         & imon0,   & ! Initial month
+         & iday0,   & ! Initial day
+         & iyea,    & ! Current year
+         & imon,    & ! Current month
+         & iday,    & ! Current day
+         & idaystp, & ! Number of days between initial and current date
+         & idaycnt    ! Day counter
+
+      INTEGER, DIMENSION(12) ::   imonth_len    !: length in days of the months of the current year
+
+      !-----------------------------------------------------------------------
+      ! Compute the calendar date YYYYMMDD
+      !-----------------------------------------------------------------------
+
+      ! Initial date
+      iyea0 =   kdate0 / 10000
+      imon0 = ( kdate0 - ( iyea0 * 10000 ) ) / 100
+      iday0 =   kdate0 - ( iyea0 * 10000 ) - ( imon0 * 100 ) 
+
+      ! Check that kt >= kit000 - 1
+      IF ( kt < kit000 - 1 ) CALL ctl_stop( ' kt must be >= kit000 - 1')
+
+      ! If kt = kit000 - 1 then set the date to the restart date
+      IF ( kt == kit000 - 1 ) THEN
+
+         kdate = ndastp
+         RETURN
+
+      ENDIF
+
+      ! Compute the number of days from the initial date
+      idaystp = INT( REAL( kt - kit000 ) * rdt / 86400. )
+   
+      iday    = iday0
+      imon    = imon0
+      iyea    = iyea0
+      idaycnt = 0
+
+      CALL calc_month_len( iyea, imonth_len )
+
+      DO WHILE ( idaycnt < idaystp )
+         iday = iday + 1
+         IF ( iday > imonth_len(imon) )  THEN
+            iday = 1
+            imon = imon + 1
+         ENDIF
+         IF ( imon > 12 ) THEN
+            imon = 1
+            iyea = iyea + 1
+            CALL calc_month_len( iyea, imonth_len )  ! update month lengths
+         ENDIF                 
+         idaycnt = idaycnt + 1
+      END DO
+
+      kdate = iyea * 10000 + imon * 100 + iday
+
+   END SUBROUTINE
+
+   SUBROUTINE calc_month_len( iyear, imonth_len )
+      !!----------------------------------------------------------------------
+      !!                    ***  ROUTINE calc_month_len  ***
+      !!          
+      !! ** Purpose : Compute the number of days in a months given a year.
+      !!
+      !! ** Method  : 
+      !!
+      !! ** Action  : 
+      !!
+      !! History :
+      !!        !  10-05  (D. Lea)   New routine based on day_init 
+      !!----------------------------------------------------------------------
+
+      INTEGER, DIMENSION(12) ::   imonth_len    !: length in days of the months of the current year
+      INTEGER :: iyear         !: year
+
+      ! length of the month of the current year (from nleapy, read in namelist)
+      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
+            IF ( MOD(iyear, 4) == 0 .AND. ( MOD(iyear, 400) == 0 .OR. MOD(iyear, 100) /= 0 ) ) THEN
+               imonth_len(2) = 29
+            ENDIF
+         ENDIF
+      ELSE
+         imonth_len(:) = nleapy   ! all months with nleapy days per year
+      ENDIF
+
+   END SUBROUTINE
+
+   SUBROUTINE tra_asm_inc( kt )
+      !!----------------------------------------------------------------------
+      !!                    ***  ROUTINE tra_asm_inc  ***
+      !!          
+      !! ** Purpose : Apply the tracer (T and S) assimilation increments
+      !!
+      !! ** Method  : Direct initialization or Incremental Analysis Updating
+      !!
+      !! ** Action  : 
+      !!
+      !! History :
+      !!        !  07-03  (M. Martin) Met Office version
+      !!        !  07-04  (A. Weaver) Merge with OPAVAR/NEMOVAR
+      !!----------------------------------------------------------------------
+
+      IMPLICIT NONE
+
+      !! * Arguments
+      INTEGER, INTENT(IN) :: &
+         & kt               ! Current time step
+
+      !! * Local declarations
+      INTEGER :: &
+         & ji,jj,jk
+      INTEGER :: &
+         & it
+      REAL(wp) :: &
+         & zincwgt  ! IAU weight for current time step
+
+      IF ( ln_asmiau ) THEN
+
+         !--------------------------------------------------------------------
+         ! Incremental Analysis Updating
+         !--------------------------------------------------------------------
+
+         IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN
+
+            it = kt - nit000 + 1
+            zincwgt = wgtiau(it) / rdt   ! IAU weight for the current time step
+
+            IF(lwp) THEN
+               WRITE(numout,*) 
+               WRITE(numout,*) 'tra_asm_inc : Tracer IAU at time step = ', &
+                  &  kt,' with IAU weight = ', wgtiau(it)
+               WRITE(numout,*) '~~~~~~~~~~~~'
+            ENDIF
+
+            ! Update the tracer tendencies
+            DO jk = 1, jpkm1
+               ta(:,:,jk) = ta(:,:,jk) + t_bkginc(:,:,jk) * zincwgt  
+               sa(:,:,jk) = sa(:,:,jk) + s_bkginc(:,:,jk) * zincwgt
+            END DO
+
+            ! Salinity fix
+            IF (ln_salfix) THEN
+               DO jk = 1, jpkm1
+                  DO jj = 1, jpj
+                     DO ji= 1, jpi
+                        sa(ji,jj,jk) = MAX( sa(ji,jj,jk), salfixmin )
+                     END DO
+                  ENDDO
+               END DO
+            ENDIF
+
+         ENDIF
+
+         IF ( kt == nitiaufin_r + 1  ) THEN   ! For bias crcn to work
+            DEALLOCATE( t_bkginc )
+            DEALLOCATE( s_bkginc )
+         ENDIF
+
+
+      ELSEIF ( ln_asmdin ) THEN
+
+         !--------------------------------------------------------------------
+         ! Direct Initialization
+         !--------------------------------------------------------------------
+            
+         IF ( kt == nitdin_r ) THEN
+
+            neuler = 0  ! Force Euler forward step
+
+            ! Initialize the now fields with the background + increment
+            tn(:,:,:) = t_bkg(:,:,:) + t_bkginc(:,:,:)   
+            sn(:,:,:) = s_bkg(:,:,:) + s_bkginc(:,:,:)   
+
+            ! Optional salinity fix
+            IF (ln_salfix) THEN
+               DO jk = 1, jpkm1
+                  DO jj = 1, jpj
+                     DO ji= 1, jpi
+                        sn(ji,jj,jk) = MAX( sn(ji,jj,jk), salfixmin )
+                     END DO
+                  ENDDO
+               END DO
+            ENDIF
+
+            tb(:,:,:) = tn(:,:,:)                        ! Update before fields
+            sb(:,:,:) = sn(:,:,:)
+
+            CALL eos( tb, sb, rhd, rhop )                ! Before potential and in situ densities
+         
+            IF( ln_zps .AND. .NOT. lk_c1d ) &
+               &  CALL zps_hde( nit000, tb, sb, rhd,  &  ! Partial steps: before horizontal derivative
+               &                gtu, gsu, gru,        &  ! of T, S, rd at the bottom ocean level
+               &                gtv, gsv, grv )
+
+            DEALLOCATE( t_bkginc )
+            DEALLOCATE( s_bkginc )
+            DEALLOCATE( t_bkg    )
+            DEALLOCATE( s_bkg    )
+         ENDIF
+         
+      ENDIF
+
+   END SUBROUTINE tra_asm_inc
+
+   SUBROUTINE dyn_asm_inc( kt )
+      !!----------------------------------------------------------------------
+      !!                    ***  ROUTINE dyn_asm_inc  ***
+      !!          
+      !! ** Purpose : Apply the dynamics (u and v) assimilation increments.
+      !!
+      !! ** Method  : Direct initialization or Incremental Analysis Updating.
+      !!
+      !! ** Action  : 
+      !!
+      !! History :
+      !!        !  07-03  (M. Martin) Met Office version
+      !!        !  07-04  (A. Weaver) Merge with OPAVAR/NEMOVAR
+      !!----------------------------------------------------------------------
+
+      IMPLICIT NONE
+
+      !! * Arguments
+      INTEGER, INTENT(IN) :: &
+         & kt   ! Current time step
+
+      !! * Local declarations
+      INTEGER :: &
+         & jk
+      INTEGER :: &
+         & it
+      REAL(wp) :: &
+         & zincwgt  ! IAU weight for current time step
+
+      IF ( ln_asmiau ) THEN
+
+         !--------------------------------------------------------------------
+         ! Incremental Analysis Updating
+         !--------------------------------------------------------------------
+
+         IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN
+
+            it = kt - nit000 + 1
+            zincwgt = wgtiau(it) / rdt   ! IAU weight for the current time step
+
+            IF(lwp) THEN
+               WRITE(numout,*) 
+               WRITE(numout,*) 'dyn_asm_inc : Dynamics IAU at time step = ', &
+                  &  kt,' with IAU weight = ', wgtiau(it)
+               WRITE(numout,*) '~~~~~~~~~~~~'
+            ENDIF
+
+            ! Update the dynamic tendencies
+            DO jk = 1, jpkm1
+               ua(:,:,jk) = ua(:,:,jk) + u_bkginc(:,:,jk) * zincwgt
+               va(:,:,jk) = va(:,:,jk) + v_bkginc(:,:,jk) * zincwgt
+            END DO
+           
+            IF ( kt == nitiaufin_r ) THEN
+               DEALLOCATE( u_bkginc )
+               DEALLOCATE( v_bkginc )
+            ENDIF
+
+         ENDIF
+
+      ELSEIF ( ln_asmdin ) THEN 
+
+         !--------------------------------------------------------------------
+         ! Direct Initialization
+         !--------------------------------------------------------------------
+         
+         IF ( kt == nitdin_r ) THEN
+
+            neuler = 0                    ! Force Euler forward step
+
+            ! Initialize the now fields with the background + increment
+            un(:,:,:) = u_bkg(:,:,:) + u_bkginc(:,:,:)
+            vn(:,:,:) = v_bkg(:,:,:) + v_bkginc(:,:,:)  
+
+            ub(:,:,:) = un(:,:,:)         ! Update before fields
+            vb(:,:,:) = vn(:,:,:)
+ 
+            CALL div_cur( kt )                  ! Compute divergence and curl for now fields
+            IF( n_cla == 1 ) CALL div_cla( kt ) ! Cross Land Advection (Update Hor. divergence)
+
+            rotb (:,:,:) = rotn (:,:,:)   ! Update before fields
+            hdivb(:,:,:) = hdivn(:,:,:)
+
+            DEALLOCATE( u_bkg    )
+            DEALLOCATE( v_bkg    )
+            DEALLOCATE( u_bkginc )
+            DEALLOCATE( v_bkginc )
+
+         ENDIF
+
+      ENDIF
+
+   END SUBROUTINE dyn_asm_inc
+
+   SUBROUTINE ssh_asm_inc( kt )
+      !!----------------------------------------------------------------------
+      !!                    ***  ROUTINE ssh_asm_inc  ***
+      !!          
+      !! ** Purpose : Apply the sea surface height assimilation increment.
+      !!
+      !! ** Method  : Direct initialization or Incremental Analysis Updating.
+      !!
+      !! ** Action  : 
+      !!
+      !! History :
+      !!        !  07-03  (M. Martin) Met Office version
+      !!        !  07-04  (A. Weaver) Merge with OPAVAR/NEMOVAR
+      !!----------------------------------------------------------------------
+
+      IMPLICIT NONE
+
+      !! * Arguments
+      INTEGER, INTENT(IN) :: &
+         & kt   ! Current time step
+
+      !! * Local declarations
+      INTEGER :: &
+         & it
+      REAL(wp) :: &
+         & zincwgt  ! IAU weight for current time step
+
+      IF ( ln_asmiau ) THEN
+
+         !--------------------------------------------------------------------
+         ! Incremental Analysis Updating
+         !--------------------------------------------------------------------
+
+         IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN
+
+            it = kt - nit000 + 1
+            zincwgt = wgtiau(it) / rdt   ! IAU weight for the current time step
+
+            IF(lwp) THEN
+               WRITE(numout,*) 
+               WRITE(numout,*) 'ssh_asm_inc : SSH IAU at time step = ', &
+                  &  kt,' with IAU weight = ', wgtiau(it)
+               WRITE(numout,*) '~~~~~~~~~~~~'
+            ENDIF
+
+            ! Save the tendency associated with the IAU weighted SSH increment
+            ! (applied in dynspg.*)
+#if defined key_asminc
+            ssh_iau(:,:) = ssh_bkginc(:,:) * zincwgt
+#endif
+            IF ( kt == nitiaufin_r ) THEN
+               DEALLOCATE( ssh_bkginc )
+            ENDIF
+
+         ENDIF
+
+      ELSEIF ( ln_asmdin ) THEN
+
+         !--------------------------------------------------------------------
+         ! Direct Initialization
+         !--------------------------------------------------------------------
+
+         IF ( kt == nitdin_r ) THEN
+
+            neuler = 0                    ! Force Euler forward step
+
+            ! Initialize the now fields the background + increment
+            sshn(:,:) = ssh_bkg(:,:) + ssh_bkginc(:,:)  
+
+            sshb(:,:) = sshn(:,:)         ! Update before fields
+
+            DEALLOCATE( ssh_bkg    )
+            DEALLOCATE( ssh_bkginc )
+
+         ENDIF
+
+      ENDIF
+
+   END SUBROUTINE ssh_asm_inc
+
+END MODULE asminc
Index: branches/dev_1784_ASM/NEMO/OPA_SRC/ASM/asmpar.F90
===================================================================
--- branches/dev_1784_ASM/NEMO/OPA_SRC/ASM/asmpar.F90	(revision 2002)
+++ branches/dev_1784_ASM/NEMO/OPA_SRC/ASM/asmpar.F90	(revision 2002)
@@ -0,0 +1,33 @@
+MODULE asmpar
+   !!======================================================================
+   !!                       ***  MODULE asmpar  ***
+   !! Assimilation increment : Parameters for assimilation interface
+   !!======================================================================
+
+   IMPLICIT NONE
+
+   !! * Routine accessibility
+   PRIVATE
+
+   !! * Shared Modules variables
+   CHARACTER (LEN=40), PUBLIC, PARAMETER :: &
+      & c_asmbkg = 'assim_background_state_Jb',  & !: Filename for storing the 
+                                                   !: background state for use 
+                                                   !: in the Jb term
+      & c_asmdin = 'assim_background_state_DI',  & !: Filename for storing the 
+                                                   !: background state for direct 
+                                                   !: initialization
+      & c_asmtrj = 'assim_trj',                  & !: Filename for storing the 
+                                                   !: reference trajectory
+      & c_asminc = 'assim_background_increments'   !: Filename for storing the 
+                                                   !: increments to the background
+                                                   !: state
+
+   INTEGER, PUBLIC :: &          
+      & nitbkg_r,    &    !: Background time step referenced to nit000
+      & nitdin_r,    &    !: Direct Initialization time step referenced to nit000
+      & nitiaustr_r, &    !: IAU starting time step referenced to nit000
+      & nitiaufin_r, &    !: IAU final time step referenced to nit000
+      & nittrjfrq         !: Frequency of trajectory output for 4D-VAR
+
+END MODULE asmpar
Index: branches/dev_1784_ASM/NEMO/OPA_SRC/ASM/asmtrj.F90
===================================================================
--- branches/dev_1784_ASM/NEMO/OPA_SRC/ASM/asmtrj.F90	(revision 2002)
+++ branches/dev_1784_ASM/NEMO/OPA_SRC/ASM/asmtrj.F90	(revision 2002)
@@ -0,0 +1,293 @@
+MODULE asmtrj
+   !!======================================================================
+   !!                       ***  MODULE asmtrj  ***
+   !! Assimilation trajectory interface: Write to file the background state
+   !!                                    and the model state trajectory
+   !!======================================================================
+
+   !!----------------------------------------------------------------------
+   !!   'key_asminc' : Switch on the assimilation increment interface
+   !!----------------------------------------------------------------------
+   !!   asm_bkg_wri  : Write out the background state
+   !!   asm_trj_wri  : Write out the model state trajectory (used with 4D-Var)
+   !!----------------------------------------------------------------------
+   !! * Modules used   
+   USE oce                ! Dynamics and active tracers defined in memory
+   USE sbc_oce            ! Ocean surface boundary conditions
+   USE zdf_oce            ! Vertical mixing variables
+   USE zdfddm             ! Double diffusion mixing parameterization
+   USE ldftra_oce         ! Lateral tracer mixing coefficient defined in memory
+   USE ldfslp             ! Slopes of neutral surfaces
+   USE tradmp             ! Tracer damping
+
+#if defined key_zdftke
+   USE zdftke             ! TKE vertical physics
+#endif
+   USE eosbn2             ! Equation of state (eos_bn2 routine)
+   USE zdfmxl             ! Mixed layer depth
+   USE sol_oce, ONLY : &  ! Solver variables defined in memory
+      & gcx
+   USE in_out_manager, ONLY : &  ! I/O manager
+      & lwp,     &
+      & numout
+   USE dom_oce, ONLY : &
+      & ndastp
+   USE iom                 ! I/O module
+   USE asmpar              ! Parameters for the assmilation interface
+   USE zdfmxl, ONLY : &   ! mixed layer depth
+      & hmlp
+#if defined key_traldf_c2d
+   USE ldfeiv          ! eddy induced velocity coef.      (ldf_eiv routine)
+#endif
+
+   IMPLICIT NONE
+
+   !! * Routine accessibility
+   PRIVATE
+   PUBLIC asm_bkg_wri, &  !: Write out the background state
+      &   asm_trj_wri     !: Write out the background state
+
+CONTAINS
+
+   SUBROUTINE asm_bkg_wri( kt )
+      !!-----------------------------------------------------------------------
+      !!
+      !!                  ***  ROUTINE asm_bkg_wri ***
+      !!
+      !! ** Purpose : Write to file the background state for later use in the
+      !!              inner loop of data assimilation or for direct initialization
+      !!              in the outer loop.
+      !!
+      !! ** Method  : Write out the background state for use in the Jb term
+      !!              in the cost function and for use with direct initialization
+      !!              at analysis time.
+      !!
+      !! ** Action  :
+      !!                   
+      !! References : 
+      !!
+      !! History :
+      !!        ! 07-03 (M. Martin) Met. Office version
+      !!        ! 07-03 (K. Mogensen) Adapt to NEMOVAR and use IOM instead of IOIPSL
+      !!        ! 07-04 (A. Weaver) Name change (formally asmbkg.F90). Distinguish
+      !!                            background states in Jb term and at analysis time.
+      !!                            Include state trajectory routine (currently empty)
+      !!        ! 07-07 (A. Weaver) Add tke_rst and flt_rst for case nitbkg=0 
+      !!-----------------------------------------------------------------------
+
+      !! * Arguments
+      INTEGER, INTENT( IN ) :: &
+         & kt               ! Current time-step
+
+      !! * Local declarations
+      CHARACTER (LEN=50) :: &
+         & cl_asmbkg, &
+         & cl_asmdin
+      LOGICAL :: &
+         & llok             ! Check if file exists
+      INTEGER :: &
+         & inum             ! File unit number
+      REAL(wp) :: &
+         & zdate            ! Date
+
+      !--------------------------------------------------------------------
+      ! Write out background at time step nitbkg_r or nitdin_r
+      !--------------------------------------------------------------------
+
+      IF ( kt == nitbkg_r ) THEN
+
+         WRITE(cl_asmbkg, FMT='(A,".nc")' ) TRIM( c_asmbkg )
+         cl_asmbkg = TRIM( cl_asmbkg )
+
+         INQUIRE( FILE = cl_asmbkg, EXIST = llok )
+
+         IF( .NOT. llok ) THEN
+            IF(lwp) WRITE(numout,*) ' Setting up assimilation background file '// &
+               &                    TRIM( c_asmbkg )
+
+            ! Define the output file        
+            CALL iom_open( c_asmbkg, inum, ldwrt = .TRUE., kiolib = jprstlib)
+
+            ! Treat special case when nitbkg = 0
+            IF ( nitbkg_r == nit000 - 1 ) THEN
+               zdate = REAL( ndastp )
+#if defined key_zdftke
+               IF(lwp) WRITE(numout,*) ' Reading TKE (en) from restart...'
+               ! Read turbulent kinetic energy ( en )
+               CALL tke_rst( nit000, 'READ' )
+#endif
+            ELSE
+               zdate = REAL( ndastp )
+            ENDIF
+
+            ! Write the information
+            CALL iom_rstput( kt, nitbkg_r, inum, 'rdastp' , zdate   )
+            CALL iom_rstput( kt, nitbkg_r, inum, 'un'     , un      )
+            CALL iom_rstput( kt, nitbkg_r, inum, 'vn'     , vn      )
+            CALL iom_rstput( kt, nitbkg_r, inum, 'tn'     , tn      )
+            CALL iom_rstput( kt, nitbkg_r, inum, 'sn'     , sn      )
+            CALL iom_rstput( kt, nitbkg_r, inum, 'sshn'   , sshn    )
+#if defined key_zdftke
+            CALL iom_rstput( kt, nitbkg_r, inum, 'en'     , en      )
+#endif
+            CALL iom_rstput( kt, nitbkg_r, inum, 'gcx'    , gcx     )
+
+            CALL iom_close( inum )
+
+         ENDIF
+
+      ENDIF
+
+      IF ( kt == nitdin_r ) THEN
+
+         WRITE(cl_asmdin, FMT='(A,".nc")' ) TRIM( c_asmdin )
+         cl_asmdin = TRIM( cl_asmdin )
+
+         INQUIRE( FILE = cl_asmdin, EXIST = llok )
+
+         IF( .NOT. llok ) THEN
+            IF(lwp) WRITE(numout,*) ' Setting up assimilation background file '// &
+               &                    TRIM( c_asmdin )
+
+            ! Define the output file        
+            CALL iom_open( c_asmdin, inum, ldwrt = .TRUE., kiolib = jprstlib)
+
+            ! Treat special case when nitbkg = 0
+            IF ( nitdin_r == nit000 - 1) THEN
+               zdate = REAL( ndastp )
+            ELSE
+               zdate = REAL( ndastp )
+            ENDIF
+
+            ! Write the information
+            CALL iom_rstput( kt, nitdin_r, inum, 'rdastp' , zdate   )
+            CALL iom_rstput( kt, nitdin_r, inum, 'un'     , un      )
+            CALL iom_rstput( kt, nitdin_r, inum, 'vn'     , vn      )
+            CALL iom_rstput( kt, nitdin_r, inum, 'tn'     , tn      )
+            CALL iom_rstput( kt, nitdin_r, inum, 'sn'     , sn      )
+            CALL iom_rstput( kt, nitdin_r, inum, 'sshn'   , sshn    )
+
+            CALL iom_close( inum )
+
+         ENDIF
+
+      ENDIF
+                                  
+   END SUBROUTINE asm_bkg_wri
+
+   SUBROUTINE asm_trj_wri( kt )
+      !!-----------------------------------------------------------------------
+      !!
+      !!                  ***  ROUTINE asm_trj_wri ***
+      !!
+      !! ** Purpose : Write to file the model state trajectory for use with
+      !!              4D-Var.
+      !!
+      !! ** Method  : 
+      !!
+      !! ** Action  :
+      !!                   
+      !! References : 
+      !!
+      !! History :
+      !!        ! 07-04 (A. Weaver) 
+      !!        ! 09-03 (F. Vigilant) Add hmlp (zdfmxl) for no tracer nmldp=2 
+      !!        ! 09-06 (F. Vigilant) special case when kt=nit000-1
+      !!        ! 09-07 (F. Vigilant) add computation of eiv at restart
+      !!-----------------------------------------------------------------------
+
+      !! * Arguments
+      INTEGER, INTENT( IN ) :: &
+         & kt                    ! Current time-step
+
+      !! * Local declarations
+      INTEGER :: &
+         & inum                  ! File unit number
+      INTEGER :: &
+         & it
+      CHARACTER (LEN=50) :: &
+         & cl_asmtrj
+      REAL(wp) :: &
+         & zdate            ! Date
+
+      !------------------------------------------------------------------------
+      ! Write a single file for each trajectory time step
+      !------------------------------------------------------------------------
+      IF ( ( MOD( kt - nit000 + 1, nittrjfrq ) == 0 ) .OR. &
+         & ( kt == nitend ) ) THEN
+         
+         ! Treat special case when kt = nit000-1
+         IF ( kt == nit000 - 1 ) THEN
+#if defined key_zdftke
+            IF(lwp) WRITE(numout,*) ' Computing  zdf_tke coeff. form restart...'
+            ! Compute the vertical eddy viscosity and diffusivity coefficients
+            CALL zdf_tke( nit000 )
+#endif
+#if defined key_zdfddm
+            IF(lwp) WRITE(numout,*) ' Computing zdf_ddm coeff. from restart...'
+            ! Compute the vertical eddy viscosity and diffusivity coefficients (salt effect)
+            CALL zdf_ddm( nit000 )
+#endif
+            IF(lwp) WRITE(numout,*) ' Computing zdf_mxl coeff. from restart...'
+            ! Compute the turbocline depth and the mixed layer depth
+            CALL zdf_mxl( nit000 ) 
+#if defined key_ldfslp
+            IF(lwp) WRITE(numout,*) ' Compute the slopes of neutral surface...'
+            CALL bn2( tb, sb, rn2 ) 
+            CALL ldf_slp( nit000, rhd, rn2 )
+#endif
+#if defined key_traldf_c2d
+            IF(lwp) WRITE(numout,*) ' Computing ldf_eiv coeff. from restart...'
+            ! Compute eddy induced velocity coefficient
+            IF( lk_traldf_eiv )   CALL ldf_eiv( nit000 )
+#endif
+          ENDIF
+
+
+         it = kt - nit000 + 1
+
+         ! Define the output file        
+         WRITE(cl_asmtrj, FMT='(A,A,I5.5)' ) TRIM( c_asmtrj ), '_', it
+         cl_asmtrj = TRIM( cl_asmtrj )
+         CALL iom_open( cl_asmtrj, inum, ldwrt = .TRUE., kiolib = jprstlib)
+
+         ! Output trajectory fields
+         CALL iom_rstput( it, it, inum, 'emp'   , emp    )
+         CALL iom_rstput( it, it, inum, 'emps'  , emps   )
+         CALL iom_rstput( it, it, inum, 'un'    , un     )
+         CALL iom_rstput( it, it, inum, 'vn'    , vn     )
+         CALL iom_rstput( it, it, inum, 'tn'    , tn     )
+         CALL iom_rstput( it, it, inum, 'sn'    , sn     )
+         CALL iom_rstput( it, it, inum, 'avmu'  , avmu   )
+         CALL iom_rstput( it, it, inum, 'avmv'  , avmv   )
+         CALL iom_rstput( it, it, inum, 'avt'   , avt    )
+#if defined key_ldfslp
+         CALL iom_rstput( it, it, inum, 'uslp'  , uslp   )
+         CALL iom_rstput( it, it, inum, 'vslp'  , vslp   )
+         CALL iom_rstput( it, it, inum, 'wslpi' , wslpi  )
+         CALL iom_rstput( it, it, inum, 'wslpj' , wslpj  )
+#endif
+#if defined key_zdfddm
+         CALL iom_rstput( it, it, inum, 'avs'   , avs    )
+#endif
+         CALL iom_rstput( it, it, inum, 'ta'    , ta     )
+         CALL iom_rstput( it, it, inum, 'sa'    , sa     )
+         CALL iom_rstput( it, it, inum, 'tb'    , tb     )
+         CALL iom_rstput( it, it, inum, 'sb'    , sb     )
+#if defined key_tradmp
+         CALL iom_rstput( it, it, inum, 'strdmp', strdmp )
+         CALL iom_rstput( it, it, inum, 'hmlp'  , hmlp   )
+#endif
+         CALL iom_rstput( it, it, inum, 'aeiu'  , aeiu   )
+         CALL iom_rstput( it, it, inum, 'aeiv'  , aeiv   )
+         CALL iom_rstput( it, it, inum, 'aeiw'  , aeiw   )
+
+         CALL iom_close( inum )
+         
+      ENDIF
+
+   END SUBROUTINE asm_trj_wri
+
+END MODULE asmtrj
+   
+   
Index: branches/dev_1784_ASM/NEMO/OPA_SRC/DYN/sshwzv.F90
===================================================================
--- branches/dev_1784_ASM/NEMO/OPA_SRC/DYN/sshwzv.F90	(revision 2001)
+++ branches/dev_1784_ASM/NEMO/OPA_SRC/DYN/sshwzv.F90	(revision 2002)
@@ -5,4 +5,5 @@
    !!==============================================================================
    !! History :  3.1  !  2009-02  (G. Madec, M. Leclair)  Original code
+   !!            3.3  !  2010-05  (K. Mogensen, A. Weaver, M. Martin, D. Lea) Assimilation interface
    !!----------------------------------------------------------------------
 
@@ -27,4 +28,5 @@
    USE diaar5, ONLY :   lk_diaar5
    USE iom
+   USE asminc          ! Assimilation increment
 
    IMPLICIT NONE
@@ -187,4 +189,10 @@
       ENDIF
 
+! Include the IAU weighted SSH increment
+      IF( ( lk_asminc ).AND.( ln_sshinc ).AND.( ln_asmiau ) ) THEN
+         CALL ssh_asm_inc( kt )
+         ssha(:,:) = ssha(:,:) + z2dt * ssh_iau(:,:)
+      ENDIF
+
       !                                           !------------------------------!
       !                                           !     Now Vertical Velocity    !
Index: branches/dev_1784_ASM/NEMO/OPA_SRC/TRA/tradmp.F90
===================================================================
--- branches/dev_1784_ASM/NEMO/OPA_SRC/TRA/tradmp.F90	(revision 2001)
+++ branches/dev_1784_ASM/NEMO/OPA_SRC/TRA/tradmp.F90	(revision 2002)
@@ -47,4 +47,6 @@
    LOGICAL, PUBLIC            ::   lk_tradmp = .TRUE.     !: internal damping flag
 #endif
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   strdmp   !: damping salinity trend (psu/s)
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   ttrdmp   !: damping temperature trend (Centigrade/s)
    REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   resto    !: restoring coeff. on T and S (s-1)
    
@@ -91,4 +93,5 @@
       INTEGER, INTENT(in) ::   kt   ! ocean time-step index
       !!
+      REAL(wp) ::   zta, zsa       ! temporary scalars
       INTEGER ::   ji, jj, jk   ! dummy loop indices
       !!----------------------------------------------------------------------
@@ -107,6 +110,12 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  ta(ji,jj,jk) = ta(ji,jj,jk) + resto(ji,jj,jk) * ( t_dta(ji,jj,jk) - tb(ji,jj,jk) )
-                  sa(ji,jj,jk) = sa(ji,jj,jk) + resto(ji,jj,jk) * ( s_dta(ji,jj,jk) - sb(ji,jj,jk) )
+                  zta = resto(ji,jj,jk) * ( t_dta(ji,jj,jk) - tb(ji,jj,jk) )
+                  zsa = resto(ji,jj,jk) * ( s_dta(ji,jj,jk) - sb(ji,jj,jk) )
+                  ! add the trends to the general tracer trends
+                  ta(ji,jj,jk) = ta(ji,jj,jk) + zta
+                  sa(ji,jj,jk) = sa(ji,jj,jk) + zsa
+                  ! save the salinity trend (used in asmtrj)
+                  strdmp(ji,jj,jk) = zsa
+                  ttrdmp(ji,jj,jk) = zta
                END DO
             END DO
@@ -118,7 +127,16 @@
                DO ji = fs_2, fs_jpim1   ! vector opt.
                   IF( avt(ji,jj,jk) <= 5.e-4 ) THEN
-                     ta(ji,jj,jk) = ta(ji,jj,jk) + resto(ji,jj,jk) * ( t_dta(ji,jj,jk) - tb(ji,jj,jk) )
-                     sa(ji,jj,jk) = sa(ji,jj,jk) + resto(ji,jj,jk) * ( s_dta(ji,jj,jk) - sb(ji,jj,jk) )
+                     zta = resto(ji,jj,jk) * ( t_dta(ji,jj,jk) - tb(ji,jj,jk) )
+                     zsa = resto(ji,jj,jk) * ( s_dta(ji,jj,jk) - sb(ji,jj,jk) )
+                  ELSE
+                     zta = 0.e0
+                     zsa = 0.e0   
                   ENDIF
+                  ! add the trends to the general tracer trends
+                  ta(ji,jj,jk) = ta(ji,jj,jk) + zta
+                  sa(ji,jj,jk) = sa(ji,jj,jk) + zsa
+                  ! save the salinity trend (used in asmtrj)
+                  strdmp(ji,jj,jk) = zsa
+                  ttrdmp(ji,jj,jk) = zta
                END DO
             END DO
@@ -130,7 +148,16 @@
                DO ji = fs_2, fs_jpim1   ! vector opt.
                   IF( fsdept(ji,jj,jk) >= hmlp (ji,jj) ) THEN
-                     ta(ji,jj,jk) = ta(ji,jj,jk) + resto(ji,jj,jk) * ( t_dta(ji,jj,jk) - tb(ji,jj,jk) )
-                     sa(ji,jj,jk) = sa(ji,jj,jk) + resto(ji,jj,jk) * ( s_dta(ji,jj,jk) - sb(ji,jj,jk) )
+                     zta = resto(ji,jj,jk) * ( t_dta(ji,jj,jk) - tb(ji,jj,jk) )
+                     zsa = resto(ji,jj,jk) * ( s_dta(ji,jj,jk) - sb(ji,jj,jk) )
+                  ELSE
+                     zta = 0.e0
+                     zsa = 0.e0   
                   ENDIF
+                  ! add the trends to the general tracer trends
+                  ta(ji,jj,jk) = ta(ji,jj,jk) + zta
+                  sa(ji,jj,jk) = sa(ji,jj,jk) + zsa
+                  ! save the salinity trend (used in asmtrj)
+                  strdmp(ji,jj,jk) = zsa
+                  ttrdmp(ji,jj,jk) = zta
                END DO
             END DO
@@ -199,4 +226,6 @@
          &   CALL ctl_stop( 'no temperature and/or salinity data define key_dtatem and key_dtasal' )
 
+      strdmp(:,:,:) = 0.e0       ! internal damping salinity trend (used in asmtrj)
+      ttrdmp(:,:,:) = 0.e0
       !                          ! Damping coefficients initialization
       IF( lzoom ) THEN   ;   CALL dtacof_zoom
Index: branches/dev_1784_ASM/NEMO/OPA_SRC/ZDF/zdftke.F90
===================================================================
--- branches/dev_1784_ASM/NEMO/OPA_SRC/ZDF/zdftke.F90	(revision 2001)
+++ branches/dev_1784_ASM/NEMO/OPA_SRC/ZDF/zdftke.F90	(revision 2002)
@@ -87,5 +87,5 @@
 
    REAL(wp), DIMENSION(jpi,jpj)     ::   htau      ! depth of tke penetration (nn_htau)
-   REAL(wp), DIMENSION(jpi,jpj,jpk) ::   en        ! now turbulent kinetic energy   [m2/s2]
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   en        ! now turbulent kinetic energy   [m2/s2]
    REAL(wp), DIMENSION(jpi,jpj,jpk) ::   dissl     ! now mixing lenght of dissipation
 
Index: branches/dev_1784_ASM/NEMO/OPA_SRC/opa.F90
===================================================================
--- branches/dev_1784_ASM/NEMO/OPA_SRC/opa.F90	(revision 2001)
+++ branches/dev_1784_ASM/NEMO/OPA_SRC/opa.F90	(revision 2002)
@@ -25,4 +25,5 @@
    !!             -   ! 2007-07  (J. Chanut, A. Sellar) Unstructured open boundaries (BDY)
    !!            3.2  ! 2009-08  (S. Masson)  open/write in the listing file in mpp
+   !!            3.3  ! 2010-05  (K. Mogensen, A. Weaver, M. Martin, D. Lea) Assimilation interface 
    !!----------------------------------------------------------------------
 
@@ -53,4 +54,7 @@
    USE phycst          ! physical constant                  (par_cst routine)
    USE trdmod          ! momentum/tracers trends       (trd_mod_init routine)
+   USE asminc          ! assimilation increments       (asm_inc_init routine)
+   USE asmtrj          ! writing out state trajectory
+   USE sshwzv          ! vertical velocity used in asm
    USE diaptr          ! poleward transports           (dia_ptr_init routine)
    USE step            ! OPA time-stepping                  (stp     routine)
@@ -133,4 +137,17 @@
          END DO
       ELSE                              !==  3D ocean with  ==!
+          IF( lk_asminc ) THEN
+             IF( ln_bkgwri ) CALL asm_bkg_wri( nit000 - 1 )    ! Output background fields
+             IF( ln_trjwri ) CALL asm_trj_wri( nit000 - 1 )    ! Output trajectory fields
+             IF( ln_asmdin ) THEN                        ! Direct initialization
+                IF( ln_trainc ) CALL tra_asm_inc( nit000 - 1 )    ! Tracers
+                IF( ln_dyninc ) THEN 
+                   CALL dyn_asm_inc( nit000 - 1 )    ! Dynamics
+                   IF ( ln_asmdin ) CALL ssh_wzv ( nit000 - 1 )      ! update vertical velocity 
+                ENDIF
+                IF( ln_sshinc ) CALL ssh_asm_inc( nit000 - 1 )    ! SSH
+             ENDIF
+          ENDIF
+        
          DO WHILE ( istp <= nitend .AND. nstop == 0 )
 #if defined key_agrif
@@ -271,4 +288,8 @@
       CALL dia_ptr_init                         ! Poleward TRansports initialization
       CALL trd_mod_init                         ! Mixed-layer/Vorticity/Integral constraints trends
+
+      IF( lk_asminc ) CALL asm_inc_init     ! Initialize assimilation increments
+      IF(lwp) WRITE(numout,*)'Euler time step switch is ', neuler
+
       !
    END SUBROUTINE opa_init
Index: branches/dev_1784_ASM/NEMO/OPA_SRC/step.F90
===================================================================
--- branches/dev_1784_ASM/NEMO/OPA_SRC/step.F90	(revision 2001)
+++ branches/dev_1784_ASM/NEMO/OPA_SRC/step.F90	(revision 2002)
@@ -21,4 +21,5 @@
    !!            3.2  !  2009-02  (G. Madec, R. Benshila)  reintroduicing z*-coordinate
    !!             -   !  2009-06  (S. Masson, G. Madec)  TKE restart compatible with key_cpl
+   !!            3.3  !  2010-05  (K. Mogensen, A. Weaver, M. Martin, D. Lea) Assimilation interface
    !!----------------------------------------------------------------------
 
@@ -113,4 +114,5 @@
    USE floats          ! floats computation               (flo_stp routine)
 
+   USE asminc          ! assimilation increments    (tra_asm_inc, dyn_asm_inc routines)
    USE stpctl          ! time stepping control            (stp_ctl routine)
    USE restart         ! ocean restart                    (rst_wri routine)
@@ -262,4 +264,6 @@
                              sa(:,:,:) = 0.e0
 
+      IF(  ln_asmiau .AND. &
+         & ln_trainc     )   CALL tra_asm_inc( kstp )       ! apply tracer assimilation increment
                              CALL tra_sbc    ( kstp )       ! surface boundary condition
       IF( ln_traqsr      )   CALL tra_qsr    ( kstp )       ! penetrative solar radiation qsr
@@ -300,4 +304,6 @@
                                va(:,:,:) = 0.e0
 
+      IF(  ln_asmiau .AND. &
+         & ln_dyninc     )     CALL dyn_asm_inc( kstp )       ! apply dynamics assimilation increment
                                CALL dyn_adv( kstp )         ! advection (vector or flux form)
                                CALL dyn_vor( kstp )         ! vorticity term including Coriolis
