ecearth3/sources/nemo-3.6/NEMO/OPA_SRC/ASM/asminc.F90

MODULE asminc
   !!======================================================================
   !!                       ***  MODULE asminc  ***
   !! Assimilation increment : Apply an increment generated by data
   !!                          assimilation
   !!======================================================================
   !! History :       ! 2007-03  (M. Martin)  Met Office version
   !!                 ! 2007-04  (A. Weaver)  calc_date original code
   !!                 ! 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 
   !!            3.4  ! 2012-10  (A. Weaver and K. Mogensen) Fix for direct initialization
   !!----------------------------------------------------------------------

   !!----------------------------------------------------------------------
   !!   '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
   !!   seaice_asm_inc : Apply the seaice increment
   !!----------------------------------------------------------------------
   USE wrk_nemo         ! Memory Allocation
   USE par_oce          ! Ocean space and time domain variables
   USE dom_oce          ! Ocean space and time domain
   USE domvvl           ! domain: variable volume level
   USE oce              ! Dynamics and active tracers defined in memory
   USE ldfdyn_oce       ! ocean dynamics: lateral physics
   USE eosbn2           ! Equation of state - in situ and potential density
   USE zpshde           ! Partial step : Horizontal Derivative
   USE iom              ! Library to read input files
   USE asmpar           ! Parameters for the assmilation interface
   USE c1d              ! 1D initialization
   USE in_out_manager   ! I/O manager
   USE lib_mpp          ! MPP library
#if defined key_lim2
   USE ice_2            ! LIM2
#endif
   USE sbc_oce          ! Surface boundary condition variables.

   IMPLICIT NONE
   PRIVATE
   
   PUBLIC   asm_inc_init   !: Initialize the increment arrays and IAU weights
   PUBLIC   calc_date      !: Compute the calendar date YYYYMMDD on a given step
   PUBLIC   tra_asm_inc    !: Apply the tracer (T and S) increments
   PUBLIC   dyn_asm_inc    !: Apply the dynamic (u and v) increments
   PUBLIC   ssh_asm_inc    !: Apply the SSH increment
   PUBLIC   seaice_asm_inc !: Apply the seaice increment

#if defined key_asminc
    LOGICAL, PUBLIC, PARAMETER :: lk_asminc = .TRUE.   !: Logical switch for assimilation increment interface
#else
    LOGICAL, PUBLIC, PARAMETER :: lk_asminc = .FALSE.  !: No assimilation increments
#endif
   LOGICAL, PUBLIC :: ln_bkgwri = .FALSE.      !: No output of the background state fields
   LOGICAL, PUBLIC :: ln_asmiau = .FALSE.      !: No applying forcing with an assimilation increment
   LOGICAL, PUBLIC :: ln_asmdin = .FALSE.      !: No direct initialization
   LOGICAL, PUBLIC :: ln_trainc = .FALSE.      !: No tracer (T and S) assimilation increments
   LOGICAL, PUBLIC :: ln_dyninc = .FALSE.      !: No dynamics (u and v) assimilation increments
   LOGICAL, PUBLIC :: ln_sshinc = .FALSE.      !: No sea surface height assimilation increment
   LOGICAL, PUBLIC :: ln_seaiceinc             !: No sea ice concentration increment
   LOGICAL, PUBLIC :: ln_salfix = .FALSE.      !: Apply minimum salinity check
   LOGICAL, PUBLIC :: ln_temnofreeze = .FALSE. !: Don't allow the temperature to drop below freezing
   INTEGER, PUBLIC :: nn_divdmp                !: Apply divergence damping filter nn_divdmp times

   REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE ::   t_bkg   , s_bkg      !: Background temperature and salinity
   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 ::   wgtiau               !: IAU weights for each time step
#if defined key_asminc
   REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE ::   ssh_iau           !: IAU-weighted sea surface height increment
#endif
   !                                !!! time steps relative to the cycle interval [0,nitend-nit000-1]
   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 ::   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 
   REAL(wp), PUBLIC ::   salfixmin   !: Ensure that the salinity is larger than this  value if (ln_salfix)

   REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   ssh_bkg, ssh_bkginc   ! Background sea surface height and its increment
   REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   seaice_bkginc         ! Increment to the background sea ice conc
#if defined key_cice && defined key_asminc
   REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   ndaice_da             ! ice increment tendency into CICE
#endif

   !! * Substitutions
#  include "domzgr_substitute.h90"
#  include "ldfdyn_substitute.h90"
#  include "vectopt_loop_substitute.h90"
   !!----------------------------------------------------------------------
   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
   !! $Id: asminc.F90 5540 2015-07-02 15:11:23Z jchanut $
   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
   !!----------------------------------------------------------------------
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, jt  ! dummy loop indices
      INTEGER :: imid, inum      ! local integers
      INTEGER :: ios             ! Local integer output status for namelist read
      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   ! 2D workspace
      !!
      NAMELIST/nam_asminc/ ln_bkgwri,                                      &
         &                 ln_trainc, ln_dyninc, ln_sshinc,                &
         &                 ln_asmdin, ln_asmiau,                           &
         &                 nitbkg, nitdin, nitiaustr, nitiaufin, niaufn,   &
         &                 ln_salfix, salfixmin, nn_divdmp
      !!----------------------------------------------------------------------

      !-----------------------------------------------------------------------
      ! Read Namelist nam_asminc : assimilation increment interface
      !-----------------------------------------------------------------------
      ln_seaiceinc = .FALSE.
      ln_temnofreeze = .FALSE.

      REWIND( numnam_ref )              ! Namelist nam_asminc in reference namelist : Assimilation increment
      READ  ( numnam_ref, nam_asminc, IOSTAT = ios, ERR = 901)
901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_asminc in reference namelist', lwp )

      REWIND( numnam_cfg )              ! Namelist nam_asminc in configuration namelist : Assimilation increment
      READ  ( numnam_cfg, nam_asminc, IOSTAT = ios, ERR = 902 )
902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_asminc in configuration namelist', lwp )
      IF(lwm) WRITE ( numond, 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 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 applying sea ice increments        ln_seaiceinc = ', ln_seaiceinc
         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,*) '      Logical switch for ensuring that the sa > salfixmin      ln_salfix = ', 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

      CALL calc_date( nit000, nitend     , ndate0, iitend_date    )     ! Date of final time step
      CALL calc_date( nit000, nitbkg_r   , ndate0, iitbkg_date    )     ! Background time for Jb referenced to ndate0
      CALL calc_date( nit000, nitdin_r   , ndate0, iitdin_date    )     ! Background time for DI referenced to ndate0
      CALL calc_date( nit000, nitiaustr_r, ndate0, iitiaustr_date )     ! IAU start time referenced to ndate0
      CALL calc_date( nit000, nitiaufin_r, ndate0, iitiaufin_date )     ! IAU end time referenced to ndate0
      !
      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,*)
         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 ) .OR. ( ln_seaiceinc) )) &
         & CALL ctl_stop( ' One or more of ln_trainc, ln_dyninc, ln_sshinc and ln_seaiceinc is set to .true.', &
         &                ' but ln_asmdin and ln_asmiau are both set to .false. :', &
         &                ' Inconsistent options')

      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 ).AND.( .NOT. ln_seaiceinc ) &
         &                     )  &
         & CALL ctl_warn( ' ln_trainc, ln_dyninc, ln_sshinc and ln_seaiceinc 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)   )
      ALLOCATE( seaice_bkginc(jpi,jpj))
#if defined key_asminc
      ALLOCATE( ssh_iau(jpi,jpj)      )
#endif
#if defined key_cice && defined key_asminc
      ALLOCATE( ndaice_da(jpi,jpj)      )
#endif
      t_bkginc(:,:,:) = 0.0
      s_bkginc(:,:,:) = 0.0
      u_bkginc(:,:,:) = 0.0
      v_bkginc(:,:,:) = 0.0
      ssh_bkginc(:,:) = 0.0
      seaice_bkginc(:,:) = 0.0
#if defined key_asminc
      ssh_iau(:,:)    = 0.0
#endif
#if defined key_cice && defined key_asminc
      ndaice_da(:,:) = 0.0
#endif
      IF ( ( ln_trainc ).OR.( ln_dyninc ).OR.( ln_sshinc ).OR.( ln_seaiceinc ) ) 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

         IF ( ln_seaiceinc ) THEN
            CALL iom_get( inum, jpdom_autoglo, 'bckinseaice', seaice_bkginc, 1 )
            ! Apply the masks
            seaice_bkginc(:,:) = seaice_bkginc(:,:) * tmask(:,:,1)
            ! Set missing increments to 0.0 rather than 1e+20
            ! to allow for differences in masks
            WHERE( ABS( seaice_bkginc(:,:) ) > 1.0e+10 ) seaice_bkginc(:,:) = 0.0
         ENDIF

         CALL iom_close( inum )
 
      ENDIF

      !-----------------------------------------------------------------------
      ! Apply divergence damping filter
      !-----------------------------------------------------------------------

      IF ( ln_dyninc .AND. nn_divdmp > 0 ) THEN

         CALL wrk_alloc(jpi,jpj,hdiv) 

         DO  jt = 1, nn_divdmp

            DO jk = 1, jpkm1

               hdiv(:,:) = 0._wp

               DO jj = 2, jpjm1
                  DO ji = fs_2, fs_jpim1   ! vector opt.
                     hdiv(ji,jj) =   &
                        (  e2u(ji  ,jj  ) * fse3u(ji  ,jj  ,jk) * u_bkginc(ji  ,jj  ,jk)     &
                         - e2u(ji-1,jj  ) * fse3u(ji-1,jj  ,jk) * u_bkginc(ji-1,jj  ,jk)     &
                         + e1v(ji  ,jj  ) * fse3v(ji  ,jj  ,jk) * v_bkginc(ji  ,jj  ,jk)     &
                         - e1v(ji  ,jj-1) * fse3v(ji  ,jj-1,jk) * v_bkginc(ji  ,jj-1,jk)  )  &
                         / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
                  END DO
               END DO

               CALL lbc_lnk( hdiv, 'T', 1. )   ! lateral boundary cond. (no sign change)

               DO jj = 2, jpjm1
                  DO ji = fs_2, fs_jpim1   ! vector opt.
                     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) 
                     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) 
                  END DO
               END DO

            END DO

         END DO

         CALL wrk_dealloc(jpi,jpj,hdiv) 

      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, 'rdastp', 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  : 
      !!----------------------------------------------------------------------
      INTEGER, INTENT(IN) :: kit000  ! Initial time step
      INTEGER, INTENT(IN) :: kt      ! Current time step referenced to kit000
      INTEGER, INTENT(IN) :: kdate0  ! Initial date
      INTEGER, INTENT(OUT) :: kdate  ! Current date reference to kdate0
      !
      INTEGER :: iyea0    ! Initial year
      INTEGER :: imon0    ! Initial month
      INTEGER :: iday0    ! Initial day
      INTEGER :: iyea     ! Current year
      INTEGER :: imon     ! Current month
      INTEGER :: iday     ! Current day
      INTEGER :: idaystp  ! Number of days between initial and current date
      INTEGER :: 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  : 
      !!----------------------------------------------------------------------
      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  : 
      !!----------------------------------------------------------------------
      INTEGER, INTENT(IN) :: kt               ! Current time step
      !
      INTEGER :: ji,jj,jk
      INTEGER :: it
      REAL(wp) :: zincwgt  ! IAU weight for current time step
      REAL (wp), DIMENSION(jpi,jpj,jpk) :: fzptnz ! 3d freezing point values
      !!----------------------------------------------------------------------

      ! 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 

      DO jk = 1, jpkm1
        CALL eos_fzp( tsn(:,:,jk,jp_sal), fzptnz(:,:,jk), fsdept(:,:,jk) )
      END DO

      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
               IF (ln_temnofreeze) THEN
                  ! 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  
                  END WHERE
               ELSE
                  tsa(:,:,jk,jp_tem) = tsa(:,:,jk,jp_tem) + t_bkginc(:,:,jk) * zincwgt  
               ENDIF
               IF (ln_salfix) THEN
                  ! Do not apply negative increments if the salinity will fall below a specified
                  ! 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 ) 
                     tsa(:,:,jk,jp_sal) = tsa(:,:,jk,jp_sal) + s_bkginc(:,:,jk) * zincwgt
                  END WHERE
               ELSE
                  tsa(:,:,jk,jp_sal) = tsa(:,:,jk,jp_sal) + s_bkginc(:,:,jk) * zincwgt
               ENDIF
            END DO

         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
            IF (ln_temnofreeze) THEN
               ! 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(:,:,:)   
               END WHERE
            ELSE
               tsn(:,:,:,jp_tem) = t_bkg(:,:,:) + t_bkginc(:,:,:)   
            ENDIF
            IF (ln_salfix) THEN
               ! Do not apply negative increments if the salinity will fall below a specified
               ! minimum value salfixmin
               WHERE( s_bkginc(:,:,:) > 0.0_wp .OR. tsn(:,:,:,jp_sal) + s_bkginc(:,:,:) > salfixmin ) 
                  tsn(:,:,:,jp_sal) = s_bkg(:,:,:) + s_bkginc(:,:,:)   
               END WHERE
            ELSE
               tsn(:,:,:,jp_sal) = s_bkg(:,:,:) + s_bkginc(:,:,:)   
            ENDIF

            tsb(:,:,:,:) = tsn(:,:,:,:)                 ! Update before fields

            CALL eos( tsb, rhd, rhop, gdept_0(:,:,:) )  ! Before potential and in situ densities
!!gm  fabien
!            CALL eos( tsb, rhd, rhop )                ! Before potential and in situ densities
!!gm


            IF( ln_zps .AND. .NOT. lk_c1d .AND. .NOT. ln_isfcav)      &
               &  CALL zps_hde    ( kt, jpts, tsb, gtsu, gtsv,        &  ! Partial steps: before horizontal gradient
               &                              rhd, gru , grv          )  ! of t, s, rd at the last ocean level
            IF( ln_zps .AND. .NOT. lk_c1d .AND.       ln_isfcav)      &
               &  CALL zps_hde_isf( nit000, jpts, tsb, gtsu, gtsv,    &    ! Partial steps for top cell (ISF)
               &                                  rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv ,   &
               &                           gtui, gtvi, grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi    ) ! of t, s, rd at the last ocean level

#if defined key_zdfkpp
            CALL eos( tsn, rhd, fsdept_n(:,:,:) )                      ! Compute rhd
!!gm fabien            CALL eos( tsn, rhd )                      ! Compute rhd
#endif

            DEALLOCATE( t_bkginc )
            DEALLOCATE( s_bkginc )
            DEALLOCATE( t_bkg    )
            DEALLOCATE( s_bkg    )
         ENDIF
         !  
      ENDIF
      ! Perhaps the following call should be in step
      IF   ( ln_seaiceinc  )   CALL seaice_asm_inc ( kt )   ! apply sea ice concentration increment
      !
   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  : 
      !!----------------------------------------------------------------------
      INTEGER, INTENT(IN) :: kt   ! Current time step
      !
      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(:,:,:)
 
            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  : 
      !!----------------------------------------------------------------------
      INTEGER, INTENT(IN) :: kt   ! Current time step
      !
      INTEGER :: it
      INTEGER :: jk
      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
#if defined key_asminc
         ELSE IF( kt == nitiaufin_r+1 ) THEN
            !  
            ssh_iau(:,:) = 0._wp
            ! 
#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(:,:)  

            ! Update before fields
            sshb(:,:) = sshn(:,:)         

            IF( lk_vvl ) THEN
               DO jk = 1, jpk
                  fse3t_b(:,:,jk) = fse3t_n(:,:,jk)
               END DO
            ENDIF

            DEALLOCATE( ssh_bkg    )
            DEALLOCATE( ssh_bkginc )

         ENDIF
         !
      ENDIF
      !
   END SUBROUTINE ssh_asm_inc


   SUBROUTINE seaice_asm_inc( kt, kindic )
      !!----------------------------------------------------------------------
      !!                    ***  ROUTINE seaice_asm_inc  ***
      !!          
      !! ** Purpose : Apply the sea ice assimilation increment.
      !!
      !! ** Method  : Direct initialization or Incremental Analysis Updating.
      !!
      !! ** Action  : 
      !!
      !!----------------------------------------------------------------------
      IMPLICIT NONE
      !
      INTEGER, INTENT(in)           ::   kt   ! Current time step
      INTEGER, INTENT(in), OPTIONAL ::   kindic   ! flag for disabling the deallocation
      !
      INTEGER  ::   it
      REAL(wp) ::   zincwgt   ! IAU weight for current time step
#if defined key_lim2
      REAL(wp), DIMENSION(jpi,jpj) ::   zofrld, zohicif, zseaicendg, zhicifinc  ! LIM
      REAL(wp) ::   zhicifmin = 0.5_wp      ! ice minimum depth in metres
#endif
      !!----------------------------------------------------------------------

      IF ( ln_asmiau ) THEN

         !--------------------------------------------------------------------
         ! Incremental Analysis Updating
         !--------------------------------------------------------------------

         IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN

            it = kt - nit000 + 1
            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,*) 'seaice_asm_inc : sea ice conc IAU at time step = ', &
                  &  kt,' with IAU weight = ', wgtiau(it)
               WRITE(numout,*) '~~~~~~~~~~~~'
            ENDIF

            ! Sea-ice : LIM-3 case (to add)

#if defined key_lim2
            ! Sea-ice : LIM-2 case
            zofrld (:,:) = frld(:,:)
            zohicif(:,:) = hicif(:,:)
            !
            frld  = MIN( MAX( frld (:,:) - seaice_bkginc(:,:) * zincwgt, 0.0_wp), 1.0_wp)
            pfrld = MIN( MAX( pfrld(:,:) - seaice_bkginc(:,:) * zincwgt, 0.0_wp), 1.0_wp)
            fr_i(:,:) = 1.0_wp - frld(:,:)        ! adjust ice fraction
            !
            zseaicendg(:,:) = zofrld(:,:) - frld(:,:)   ! find out actual sea ice nudge applied
            !
            ! 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    
            ELSEWHERE
               zhicifinc(:,:) = 0.0_wp
            END WHERE
            !
            ! nudge ice depth
            hicif (:,:) = hicif (:,:) + zhicifinc(:,:)
            phicif(:,:) = phicif(:,:) + zhicifinc(:,:)       
            !
            ! seaice salinity balancing (to add)
#endif

#if defined key_cice && defined key_asminc
            ! Sea-ice : CICE case. Pass ice increment tendency into CICE
            ndaice_da(:,:) = seaice_bkginc(:,:) * zincwgt / rdt
#endif

            IF ( kt == nitiaufin_r ) THEN
               DEALLOCATE( seaice_bkginc )
            ENDIF

         ELSE

#if defined key_cice && defined key_asminc
            ! Sea-ice : CICE case. Zero ice increment tendency into CICE
            ndaice_da(:,:) = 0.0_wp
#endif

         ENDIF

      ELSEIF ( ln_asmdin ) THEN

         !--------------------------------------------------------------------
         ! Direct Initialization
         !--------------------------------------------------------------------

         IF ( kt == nitdin_r ) THEN

            neuler = 0                    ! Force Euler forward step

            ! Sea-ice : LIM-3 case (to add)

#if defined key_lim2
            ! Sea-ice : LIM-2 case.
            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(:,:) 
            fr_i (:,:) = 1.0_wp - frld(:,:)                ! adjust ice fraction
            zseaicendg(:,:) = zofrld(:,:) - frld(:,:)      ! find out actual sea ice nudge applied
            !
            ! 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    
            ELSEWHERE
               zhicifinc(:,:) = 0.0_wp
            END WHERE
            !
            ! nudge ice depth
            hicif (:,:) = hicif (:,:) + zhicifinc(:,:)
            phicif(:,:) = phicif(:,:)       
            !
            ! seaice salinity balancing (to add)
#endif
 
#if defined key_cice && defined key_asminc
            ! Sea-ice : CICE case. Pass ice increment tendency into CICE
           ndaice_da(:,:) = seaice_bkginc(:,:) / rdt
#endif
           IF ( .NOT. PRESENT(kindic) ) THEN
              DEALLOCATE( seaice_bkginc )
           END IF

         ELSE

#if defined key_cice && defined key_asminc
            ! Sea-ice : CICE case. Zero ice increment tendency into CICE 
            ndaice_da(:,:) = 0.0_wp
#endif
         
         ENDIF

!#if defined defined key_lim2 || defined key_cice
!
!            IF (ln_seaicebal ) THEN       
!             !! balancing salinity increments
!             !! simple case from limflx.F90 (doesn't include a mass flux)
!             !! assumption is that as ice concentration is reduced or increased
!             !! the snow and ice depths remain constant
!             !! note that snow is being created where ice concentration is being increased
!             !! - could be more sophisticated and
!             !! not do this (but would need to alter h_snow)
!
!             usave(:,:,:)=sb(:,:,:)   ! use array as a temporary store
!
!             DO jj = 1, jpj
!               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)
!           ! fwf : ice formation and melting
!
!                 zfons = ( -nfresh_da(ji,jj)*soce + nfsalt_da(ji,jj) )*rdt
!
!           ! change salinity down to mixed layer depth
!                 mld=hmld_kara(ji,jj)
!
!           ! prevent small mld
!           ! less than 10m can cause salinity instability 
!                 IF (mld < 10) mld=10
!
!           ! 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 
!                     mld=gdepw(ji,jj,jk+1)
!                     jkmax=jk
!                   ENDIF
!                 ENDDO
!
!            ! avoid applying salinity balancing in shallow water or on land
!            ! 
!
!            ! dsal_ocn (psu kg m^-2) / (kg m^-3 * m)
!
!                 dsal_ocn=0.0_wp
!                 sal_thresh=5.0_wp        ! minimum salinity threshold for salinity balancing
!
!                 if (tmask(ji,jj,1) > 0 .AND. tmask(ji,jj,jkmax) > 0 ) &
!                              dsal_ocn = zfons / (rhop(ji,jj,1) * mld)
!
!           ! 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 
!                           sb(ji,jj,jk) = sb(ji,jj,jk) + dsal_ocn
!                           sn(ji,jj,jk) = sn(ji,jj,jk) + dsal_ocn
!                     ENDIF
!
!                   ENDDO
!
!      !            !  salt exchanges at the ice/ocean interface
!      !            zpmess         = zfons / rdt_ice    ! rdt_ice is ice timestep
!      !
!      !! 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) 
!      !!                                                     ! E-P (kg m-2 s-2)
!      !            emp(ji,jj) = emp(ji,jj) + zpmess          ! E-P (kg m-2 s-2)
!               ENDDO !ji
!             ENDDO !jj!
!
!            ENDIF !ln_seaicebal
!
!#endif

      ENDIF

   END SUBROUTINE seaice_asm_inc
   
   !!======================================================================
END MODULE asminc