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
   !!                 ! 2014-09  (D. Lea)  Local calc_date removed use routine from OBS
   !!                 ! 2015-11  (D. Lea)  Handle non-zero initial time of day
   !!----------------------------------------------------------------------

   !!----------------------------------------------------------------------
   !!   asm_inc_init   : Initialize the increment arrays and IAU weights
   !!   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
   !!   ssh_asm_div    : Apply divergence associated with SSH increment
   !!   seaice_asm_inc : Apply the seaice increment
   !!----------------------------------------------------------------------
   USE oce             ! Dynamics and active tracers defined in memory
   USE par_oce         ! Ocean space and time domain variables
   USE dom_oce         ! Ocean space and time domain
   USE domvvl          ! domain: variable volume level
   USE ldfdyn          ! lateral diffusion: eddy viscosity coefficients
   USE eosbn2          ! Equation of state - in situ and potential density
   USE zpshde          ! Partial step : Horizontal Derivative
   USE asmpar          ! Parameters for the assmilation interface
   USE asmbkg          !
   USE c1d             ! 1D initialization
   USE sbc_oce         ! Surface boundary condition variables.
   USE diaobs   , ONLY : calc_date     ! Compute the calendar date on a given step
#if defined key_si3
   USE ice      , ONLY : hm_i, at_i, at_i_b
#endif
   !
   USE in_out_manager  ! I/O manager
   USE iom             ! Library to read input files
   USE lib_mpp         ! MPP library

   IMPLICIT NONE
   PRIVATE

   PUBLIC   asm_inc_init   !: Initialize the increment arrays and IAU weights
   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   ssh_asm_div    !: Apply the SSH divergence
   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     !: Logical switch for writing out background state
   LOGICAL, PUBLIC :: ln_asmiau     !: Logical switch for Incremental Analysis Updating (IAU)
   LOGICAL, PUBLIC :: ln_asmdin     !: Logical switch for Direct Initialization (DI)
   LOGICAL, PUBLIC :: ln_trainc     !: Logical switch for applying tracer increments
   LOGICAL, PUBLIC :: ln_dyninc     !: Logical switch for applying velocity increments
   LOGICAL, PUBLIC :: ln_sshinc     !: Logical switch for applying SSH increments
   LOGICAL, PUBLIC :: ln_seaiceinc  !: Logical switch for applying Sea ice concentration increments
   LOGICAL, PUBLIC :: ln_salfix     !: Logical switch for ensuring that the sa > salfixmin 
   LOGICAL, PUBLIC :: ln_temnofreeze!: Don't allow the temperature to drop below freezing
   INTEGER, PUBLIC :: nn_divdmp     !: Number of iterations of divergence damping operator

   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 "do_loop_substitute.h90"
#  include "single_precision_substitute.h90"
#  include "domzgr_substitute.h90"
   !!----------------------------------------------------------------------
   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
   !! $Id: asminc.F90 15058 2021-06-25 09:15:15Z clem $
   !! Software governed by the CeCILL license (see ./LICENSE)
   !!----------------------------------------------------------------------
CONTAINS

   SUBROUTINE asm_inc_init( Kbb, Kmm, Krhs )
      !!----------------------------------------------------------------------
      !!                    ***  ROUTINE asm_inc_init  ***
      !!
      !! ** Purpose : Initialize the assimilation increment and IAU weights.
      !!
      !! ** Method  : Initialize the assimilation increment and IAU weights.
      !!
      !! ** Action  :
      !!----------------------------------------------------------------------
      INTEGER, INTENT(in) ::  Kbb, Kmm, Krhs  ! time level indices
      !
      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
      REAL(KIND=wp) :: ditend_date     ! Date YYYYMMDD.HHMMSS of final time step
      REAL(KIND=wp) :: ditbkg_date     ! Date YYYYMMDD.HHMMSS of background time step for Jb term
      REAL(KIND=wp) :: ditdin_date     ! Date YYYYMMDD.HHMMSS of background time step for DI
      REAL(KIND=wp) :: ditiaustr_date  ! Date YYYYMMDD.HHMMSS of IAU interval start time step
      REAL(KIND=wp) :: ditiaufin_date  ! Date YYYYMMDD.HHMMSS 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), ALLOCATABLE, DIMENSION(:,:) ::   zhdiv   ! 2D workspace
      !!
      NAMELIST/nam_asminc/ ln_bkgwri,                                      &
         &                 ln_trainc, ln_dyninc, ln_sshinc,                &
         &                 ln_asmdin, ln_asmiau,                           &
         &                 nitbkg, nitdin, nitiaustr, nitiaufin, niaufn,   &
         &                 ln_seaiceinc, ln_salfix, salfixmin,             &
         &                 ln_temnofreeze, nn_divdmp
      !!----------------------------------------------------------------------

      !-----------------------------------------------------------------------
      ! Read Namelist nam_asminc : assimilation increment interface
      !-----------------------------------------------------------------------
      READ  ( numnam_ref, nam_asminc, IOSTAT = ios, ERR = 901)
901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nam_asminc in reference namelist' )
      READ  ( numnam_cfg, nam_asminc, IOSTAT = ios, ERR = 902 )
902   IF( ios >  0 )   CALL ctl_nam ( ios , 'nam_asminc in configuration namelist' )
      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
         WRITE(numout,*) '      Do not apply negative increments if the T < freezing     ln_temnofreeze = ',ln_temnofreeze
      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( nitend     , ditend_date    )   ! Date of final time step
      CALL calc_date( nitbkg_r   , ditbkg_date    )   ! Background time for Jb referenced to ndate0
      CALL calc_date( nitdin_r   , ditdin_date    )   ! Background time for DI referenced to ndate0
      CALL calc_date( nitiaustr_r, ditiaustr_date )   ! IAU start time referenced to ndate0
      CALL calc_date( nitiaufin_r, ditiaufin_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,*) '       nn_time0       = ', nn_time0
         WRITE(numout,*) '       ditend_date    = ', ditend_date
         WRITE(numout,*) '       ditbkg_date    = ', ditbkg_date
         WRITE(numout,*) '       ditdin_date    = ', ditdin_date
         WRITE(numout,*) '       ditiaustr_date = ', ditiaustr_date
         WRITE(numout,*) '       ditiaufin_date = ', ditiaufin_date
      ENDIF


      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._wp
         !
         !                                !---------------------------------------------------------
         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._wp
            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) )   ;   t_bkginc     (:,:,:) = 0._wp
      ALLOCATE( s_bkginc     (jpi,jpj,jpk) )   ;   s_bkginc     (:,:,:) = 0._wp
      ALLOCATE( u_bkginc     (jpi,jpj,jpk) )   ;   u_bkginc     (:,:,:) = 0._wp
      ALLOCATE( v_bkginc     (jpi,jpj,jpk) )   ;   v_bkginc     (:,:,:) = 0._wp
      ALLOCATE( ssh_bkginc   (jpi,jpj)     )   ;   ssh_bkginc   (:,:)   = 0._wp
      ALLOCATE( seaice_bkginc(jpi,jpj)     )   ;   seaice_bkginc(:,:)   = 0._wp
#if defined key_asminc
      ALLOCATE( ssh_iau      (jpi,jpj)     )   ;   ssh_iau      (:,:)   = 0._wp
#endif
#if defined key_cice && defined key_asminc
      ALLOCATE( ndaice_da    (jpi,jpj)     )   ;   ndaice_da    (:,:)   = 0._wp
#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 ', z_inc_dateb,' and ', z_inc_datef
            WRITE(numout,*) '~~~~~~~~~~~~'
         ENDIF
         !
         IF ( ( z_inc_dateb < ndastp + nn_time0*0.0001_wp ) .OR.   &
            & ( z_inc_datef > ditend_date ) ) &
            &    CALL ctl_warn( ' Validity time of assimilation increments is ', &
            &                   ' outside the assimilation interval' )

         IF ( ( ln_asmdin ).AND.( zdate_inc /= ditdin_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_auto, 'bckint', t_bkginc, 1 )
            CALL iom_get( inum, jpdom_auto, '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_auto, 'bckinu', u_bkginc, 1, cd_type = 'U', psgn = -1._wp )
            CALL iom_get( inum, jpdom_auto, 'bckinv', v_bkginc, 1, cd_type = 'V', psgn = -1._wp )
            ! 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_auto, '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_auto, '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
      !
      !                                            !--------------------------------------
      IF ( ln_dyninc .AND. nn_divdmp > 0 ) THEN    ! Apply divergence damping filter
         !                                         !--------------------------------------
         ALLOCATE( zhdiv(jpi,jpj) )
         !
         DO jt = 1, nn_divdmp
            !
            DO jk = 1, jpkm1           ! zhdiv = e1e1 * div
               zhdiv(:,:) = 0._wp
               DO_2D( 0, 0, 0, 0 )
                  zhdiv(ji,jj) = (  e2u(ji  ,jj) * e3u(ji  ,jj,jk,Kmm) * u_bkginc(ji  ,jj,jk)    &
                     &            - e2u(ji-1,jj) * e3u(ji-1,jj,jk,Kmm) * u_bkginc(ji-1,jj,jk)    &
                     &            + e1v(ji,jj  ) * e3v(ji,jj  ,jk,Kmm) * v_bkginc(ji,jj  ,jk)    &
                     &            - e1v(ji,jj-1) * e3v(ji,jj-1,jk,Kmm) * v_bkginc(ji,jj-1,jk)  ) &
                     &            / e3t(ji,jj,jk,Kmm)
               END_2D
               CALL lbc_lnk( 'asminc', zhdiv, 'T', 1.0_wp )   ! lateral boundary cond. (no sign change)
               !
               DO_2D( 0, 0, 0, 0 )
                  u_bkginc(ji,jj,jk) = u_bkginc(ji,jj,jk)                         &
                     &               + 0.2_wp * ( zhdiv(ji+1,jj) - zhdiv(ji  ,jj) ) * r1_e1u(ji,jj) * umask(ji,jj,jk)
                  v_bkginc(ji,jj,jk) = v_bkginc(ji,jj,jk)                         &
                     &               + 0.2_wp * ( zhdiv(ji,jj+1) - zhdiv(ji,jj  ) ) * r1_e2v(ji,jj) * vmask(ji,jj,jk)
               END_2D
            END DO
            !
         END DO
         !
         DEALLOCATE( zhdiv )
         !
      ENDIF
      !
      !                             !-----------------------------------------------------
      IF ( ln_asmdin ) THEN         ! Allocate and initialize the background state arrays
         !                          !-----------------------------------------------------
         !
         ALLOCATE( t_bkg  (jpi,jpj,jpk) )   ;   t_bkg  (:,:,:) = 0._wp
         ALLOCATE( s_bkg  (jpi,jpj,jpk) )   ;   s_bkg  (:,:,:) = 0._wp
         ALLOCATE( u_bkg  (jpi,jpj,jpk) )   ;   u_bkg  (:,:,:) = 0._wp
         ALLOCATE( v_bkg  (jpi,jpj,jpk) )   ;   v_bkg  (:,:,:) = 0._wp
         ALLOCATE( ssh_bkg(jpi,jpj)     )   ;   ssh_bkg(:,:)   = 0._wp
         !
         !
         !--------------------------------------------------------------------
         ! 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,*) '   ==>>>  Assimilation background state valid at : ', zdate_bkg
            WRITE(numout,*)
         ENDIF
         !
         IF ( zdate_bkg /= ditdin_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_auto, 'tn', t_bkg )
            CALL iom_get( inum, jpdom_auto, 'sn', s_bkg )
            t_bkg(:,:,:) = t_bkg(:,:,:) * tmask(:,:,:)
            s_bkg(:,:,:) = s_bkg(:,:,:) * tmask(:,:,:)
         ENDIF
         !
         IF ( ln_dyninc ) THEN
            CALL iom_get( inum, jpdom_auto, 'un', u_bkg, cd_type = 'U', psgn = -1._wp )
            CALL iom_get( inum, jpdom_auto, 'vn', v_bkg, cd_type = 'V', psgn = -1._wp )
            u_bkg(:,:,:) = u_bkg(:,:,:) * umask(:,:,:)
            v_bkg(:,:,:) = v_bkg(:,:,:) * vmask(:,:,:)
         ENDIF
         !
         IF ( ln_sshinc ) THEN
            CALL iom_get( inum, jpdom_auto, 'sshn', ssh_bkg )
            ssh_bkg(:,:) = ssh_bkg(:,:) * tmask(:,:,1)
         ENDIF
         !
         CALL iom_close( inum )
         !
      ENDIF
      !
      IF(lwp) WRITE(numout,*) '   ==>>>   Euler time step switch is ', l_1st_euler
      !
      IF( lk_asminc ) THEN                            !==  data assimilation  ==!
         IF( ln_bkgwri )   CALL asm_bkg_wri( nit000 - 1, Kmm )      ! Output background fields
         IF( ln_asmdin ) THEN                                  ! Direct initialization
            IF( ln_trainc )   CALL tra_asm_inc( nit000 - 1, Kbb, Kmm, ts    , Krhs )      ! Tracers
            IF( ln_dyninc )   CALL dyn_asm_inc( nit000 - 1, Kbb, Kmm, uu, vv, Krhs )      ! Dynamics
            IF( ln_sshinc )   CALL ssh_asm_inc( nit000 - 1, Kbb, Kmm )                    ! SSH
         ENDIF
      ENDIF
      !
   END SUBROUTINE asm_inc_init


   SUBROUTINE tra_asm_inc( kt, Kbb, Kmm, pts, Krhs )
      !!----------------------------------------------------------------------
      !!                    ***  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                                  , INTENT(in   ) :: Kbb, Kmm, Krhs ! Time level indices
      REAL(dp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) :: pts            ! active tracers and RHS of tracer equation
      !
      INTEGER  :: ji, jj, jk
      INTEGER  :: it
      REAL(wp) :: zincwgt  ! IAU weight for current time step
      REAL(wp), DIMENSION(A2D(nn_hls),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
      IF( ln_temnofreeze ) THEN
         DO jk = 1, jpkm1
           CALL eos_fzp( CASTSP(pts(:,:,jk,jp_sal,Kmm)), fzptnz(:,:,jk), gdept(:,:,jk,Kmm) )
         END DO
      ENDIF
         !
         !                             !--------------------------------------
      IF ( ln_asmiau ) THEN            ! Incremental Analysis Updating
         !                             !--------------------------------------
         !
         IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN
            !
            it = kt - nit000 + 1
            zincwgt = wgtiau(it) / rn_Dt   ! IAU weight for the current time step
            !
            IF( .NOT. l_istiled .OR. ntile == 1 )  THEN                       ! Do only on the first tile
               IF(lwp) THEN
                  WRITE(numout,*)
                  WRITE(numout,*) 'tra_asm_inc : Tracer IAU at time step = ', kt,' with IAU weight = ', wgtiau(it)
                  WRITE(numout,*) '~~~~~~~~~~~~'
               ENDIF
            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(A2D(0),jk) > 0.0_wp .OR. &
                     &   pts(A2D(0),jk,jp_tem,Kmm) + pts(A2D(0),jk,jp_tem,Krhs) + t_bkginc(A2D(0),jk) * wgtiau(it) > fzptnz(:,:,jk) )
                     pts(A2D(0),jk,jp_tem,Krhs) = pts(A2D(0),jk,jp_tem,Krhs) + t_bkginc(A2D(0),jk) * zincwgt
                  END WHERE
               ELSE
                  DO_2D( 0, 0, 0, 0 )
                     pts(ji,jj,jk,jp_tem,Krhs) = pts(ji,jj,jk,jp_tem,Krhs) + t_bkginc(ji,jj,jk) * zincwgt
                  END_2D
               ENDIF
               IF (ln_salfix) THEN
                  ! Do not apply negative increments if the salinity will fall below a specified
                  ! minimum value salfixmin
                  WHERE(s_bkginc(A2D(0),jk) > 0.0_wp .OR. &
                     &   pts(A2D(0),jk,jp_sal,Kmm) + pts(A2D(0),jk,jp_sal,Krhs) + s_bkginc(A2D(0),jk) * wgtiau(it) > salfixmin )
                     pts(A2D(0),jk,jp_sal,Krhs) = pts(A2D(0),jk,jp_sal,Krhs) + s_bkginc(A2D(0),jk) * zincwgt
                  END WHERE
               ELSE
                  DO_2D( 0, 0, 0, 0 )
                     pts(ji,jj,jk,jp_sal,Krhs) = pts(ji,jj,jk,jp_sal,Krhs) + s_bkginc(ji,jj,jk) * zincwgt
                  END_2D
               ENDIF
            END DO
            !
         ENDIF
         !
         IF( .NOT. l_istiled .OR. ntile == nijtile )  THEN                ! Do only on the last tile
            IF ( kt == nitiaufin_r + 1  ) THEN   ! For bias crcn to work
               DEALLOCATE( t_bkginc )
               DEALLOCATE( s_bkginc )
            ENDIF
         ENDIF
         !                             !--------------------------------------
      ELSEIF ( ln_asmdin ) THEN        ! Direct Initialization
         !                             !--------------------------------------
         !
         IF ( kt == nitdin_r ) THEN
            !
            l_1st_euler = .TRUE.  ! 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. pts(:,:,:,jp_tem,Kmm) + t_bkginc(:,:,:) > fzptnz(:,:,:) )
                  pts(:,:,:,jp_tem,Kmm) = t_bkg(:,:,:) + t_bkginc(:,:,:)
               END WHERE
            ELSE
               pts(:,:,:,jp_tem,Kmm) = 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. pts(:,:,:,jp_sal,Kmm) + s_bkginc(:,:,:) > salfixmin )
                  pts(:,:,:,jp_sal,Kmm) = s_bkg(:,:,:) + s_bkginc(:,:,:)
               END WHERE
            ELSE
               pts(:,:,:,jp_sal,Kmm) = s_bkg(:,:,:) + s_bkginc(:,:,:)
            ENDIF

            pts(:,:,:,:,Kbb) = pts(:,:,:,:,Kmm)                 ! Update before fields

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

            IF( ln_zps .AND. .NOT. ln_c1d .AND. .NOT. ln_isfcav)           &
               &  CALL zps_hde    ( kt, jpts, pts(:,:,:,:,Kbb), gtsu, gtsv,        &  ! Partial steps: before horizontal gradient
               &                              rhd, gru , grv               )  ! of t, s, rd at the last ocean level
            IF( ln_zps .AND. .NOT. ln_c1d .AND.       ln_isfcav)                       &
               &  CALL zps_hde_isf( nit000, jpts, pts(:,:,:,:,Kbb), gtsu, gtsv, gtui, gtvi,    &  ! Partial steps for top cell (ISF)
               &                                  rhd, gru , grv , grui, grvi          )  ! of t, s, rd at the last ocean level

            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, Kbb, Kmm, puu, pvv, Krhs )
      !!----------------------------------------------------------------------
      !!                    ***  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             ! ocean time-step index
      INTEGER                             , INTENT( in )  ::  Kbb, Kmm, Krhs ! ocean time level indices
      REAL(dp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) ::  puu, pvv       ! ocean velocities and RHS of momentum equation
      !
      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) / rn_Dt   ! IAU weight for the current time step
            !
            IF( .NOT. l_istiled .OR. ntile == 1 )  THEN                       ! Do only on the first tile
               IF(lwp) THEN
                  WRITE(numout,*)
                  WRITE(numout,*) 'dyn_asm_inc : Dynamics IAU at time step = ', kt,' with IAU weight = ', wgtiau(it)
                  WRITE(numout,*) '~~~~~~~~~~~~'
               ENDIF
            ENDIF
            !
            ! Update the dynamic tendencies
            DO_3D( 0, 0, 0, 0, 1, jpkm1 )
               puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + u_bkginc(ji,jj,jk) * zincwgt
               pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + v_bkginc(ji,jj,jk) * zincwgt
            END_3D
            !
            IF( .NOT. l_istiled .OR. ntile == nijtile )  THEN                ! Do only on the last tile
               IF ( kt == nitiaufin_r ) THEN
                  DEALLOCATE( u_bkginc )
                  DEALLOCATE( v_bkginc )
               ENDIF
            ENDIF
            !
         ENDIF
         !                          !-----------------------------------------
      ELSEIF ( ln_asmdin ) THEN     ! Direct Initialization
         !                          !-----------------------------------------
         !
         IF ( kt == nitdin_r ) THEN
            !
            l_1st_euler = .TRUE.                    ! Force Euler forward step
            !
            ! Initialize the now fields with the background + increment
            puu(:,:,:,Kmm) = u_bkg(:,:,:) + u_bkginc(:,:,:)
            pvv(:,:,:,Kmm) = v_bkg(:,:,:) + v_bkginc(:,:,:)
            !
            puu(:,:,:,Kbb) = puu(:,:,:,Kmm)         ! Update before fields
            pvv(:,:,:,Kbb) = pvv(:,:,:,Kmm)
            !
            DEALLOCATE( u_bkg    )
            DEALLOCATE( v_bkg    )
            DEALLOCATE( u_bkginc )
            DEALLOCATE( v_bkginc )
         ENDIF
         !
      ENDIF
      !
   END SUBROUTINE dyn_asm_inc


   SUBROUTINE ssh_asm_inc( kt, Kbb, Kmm )
      !!----------------------------------------------------------------------
      !!                    ***  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, INTENT(IN) :: Kbb, Kmm   ! Current time step
      !
      INTEGER :: it
      INTEGER :: ji, jj, 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) / rn_Dt   ! IAU weight for the current time step
            !
            IF( .NOT. l_istiled .OR. ntile == 1 )  THEN                       ! Do only on the first tile
               IF(lwp) THEN
                  WRITE(numout,*)
                  WRITE(numout,*) 'ssh_asm_inc : SSH IAU at time step = ', &
                     &  kt,' with IAU weight = ', wgtiau(it)
                  WRITE(numout,*) '~~~~~~~~~~~~'
               ENDIF
            ENDIF
            !
            ! Save the tendency associated with the IAU weighted SSH increment
            ! (applied in dynspg.*)
#if defined key_asminc
            DO_2D_OVR( nn_hls, nn_hls, nn_hls, nn_hls )
               ssh_iau(ji,jj) = ssh_bkginc(ji,jj) * zincwgt
            END_2D
#endif
            !
         ELSE IF( kt == nitiaufin_r+1 ) THEN
            !
            ! test on ssh_bkginc needed as ssh_asm_inc is called twice by time step
            IF( .NOT. l_istiled .OR. ntile == nijtile )  THEN                ! Do only on the last tile
               IF (ALLOCATED(ssh_bkginc)) DEALLOCATE( ssh_bkginc )
            ENDIF
            !
#if defined key_asminc
            DO_2D_OVR( nn_hls, nn_hls, nn_hls, nn_hls )
               ssh_iau(ji,jj) = 0._wp
            END_2D
#endif
            !
         ENDIF
         !                          !-----------------------------------------
      ELSEIF ( ln_asmdin ) THEN     ! Direct Initialization
         !                          !-----------------------------------------
         !
         IF ( kt == nitdin_r ) THEN
            !
            l_1st_euler = .TRUE.                            ! Force Euler forward step
            !
            ssh(:,:,Kmm) = ssh_bkg(:,:) + ssh_bkginc(:,:)   ! Initialize the now fields the background + increment
            !
            ssh(:,:,Kbb) = ssh(:,:,Kmm)                        ! Update before fields
#if ! defined key_qco
            e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm)
#endif
!!gm why not e3u(:,:,:,Kbb), e3v(:,:,:,Kbb), gdept(:,:,:,Kbb) ????
            !
            DEALLOCATE( ssh_bkg    )
            DEALLOCATE( ssh_bkginc )
            !
         ENDIF
         !
      ENDIF
      !
   END SUBROUTINE ssh_asm_inc


   SUBROUTINE ssh_asm_div( kt, Kbb, Kmm, phdivn )
      !!----------------------------------------------------------------------
      !!                  ***  ROUTINE ssh_asm_div  ***
      !!
      !! ** Purpose :   ssh increment with z* is incorporated via a correction of the local divergence
      !!                across all the water column
      !!
      !! ** Method  :
      !!                CAUTION : sshiau is positive (inflow) decreasing the
      !!                          divergence and expressed in m/s
      !!
      !! ** Action  :   phdivn   decreased by the ssh increment
      !!----------------------------------------------------------------------
      INTEGER, INTENT(IN) :: kt                               ! ocean time-step index
      INTEGER, INTENT(IN) :: Kbb, Kmm                         ! time level indices
      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   phdivn   ! horizontal divergence
      !!
      INTEGER  ::   ji, jj, jk                                ! dummy loop index
      REAL(wp), DIMENSION(:,:)  , POINTER       ::   ztim     ! local array
      !!----------------------------------------------------------------------
      !
#if defined key_asminc
      CALL ssh_asm_inc( kt, Kbb, Kmm ) !==   (calculate increments)
      !
      IF( ln_linssh ) THEN
         DO_2D_OVR( nn_hls-1, nn_hls, nn_hls-1, nn_hls )
            phdivn(ji,jj,1) = phdivn(ji,jj,1) - ssh_iau(ji,jj) / e3t(ji,jj,1,Kmm) * tmask(ji,jj,1)
         END_2D
      ELSE
         ALLOCATE( ztim(A2D(nn_hls)) )
         DO_2D_OVR( nn_hls-1, nn_hls, nn_hls-1, nn_hls )
            ztim(ji,jj) = ssh_iau(ji,jj) / ( ht(ji,jj) + 1.0 - ssmask(ji,jj) )
            DO jk = 1, jpkm1
               phdivn(ji,jj,jk) = phdivn(ji,jj,jk) - ztim(ji,jj) * tmask(ji,jj,jk)
            END DO
         END_2D
         !
         DEALLOCATE(ztim)
      ENDIF
#endif
      !
   END SUBROUTINE ssh_asm_div


   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  :
      !!
      !!----------------------------------------------------------------------
      INTEGER, INTENT(in)           ::   kt       ! Current time step
      INTEGER, INTENT(in), OPTIONAL ::   kindic   ! flag for disabling the deallocation
      !
      INTEGER  ::   ji, jj
      INTEGER  ::   it
      REAL(wp) ::   zincwgt   ! IAU weight for current time step
#if defined key_si3
      REAL(wp), DIMENSION(A2D(nn_hls)) ::   zofrld, zohicif, zseaicendg, zhicifinc
      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 rn_Dt (as with the tracer and other increments)
            !
            IF( .NOT. l_istiled .OR. ntile == 1 )  THEN                       ! Do only on the first tile
               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
            ENDIF
            !
            ! Sea-ice : SI3 case
            !
#if defined key_si3
            DO_2D( 0, 0, 0, 0 )
               zofrld (ji,jj) = 1._wp - at_i(ji,jj)
               zohicif(ji,jj) = hm_i(ji,jj)
               !
               at_i  (ji,jj) = 1. - MIN( MAX( 1.-at_i  (ji,jj) - seaice_bkginc(ji,jj) * zincwgt, 0.0_wp), 1.0_wp)
               at_i_b(ji,jj) = 1. - MIN( MAX( 1.-at_i_b(ji,jj) - seaice_bkginc(ji,jj) * zincwgt, 0.0_wp), 1.0_wp)
               fr_i(ji,jj) = at_i(ji,jj)        ! adjust ice fraction
               !
               zseaicendg(ji,jj) = zofrld(ji,jj) - (1. - at_i(ji,jj))   ! find out actual sea ice nudge applied
            END_2D
            !
            ! Nudge sea ice depth to bring it up to a required minimum depth
            WHERE( zseaicendg(:,:) > 0.0_wp .AND. hm_i(A2D(0)) < zhicifmin )
               zhicifinc(:,:) = (zhicifmin - hm_i(A2D(0))) * zincwgt
            ELSEWHERE
               zhicifinc(:,:) = 0.0_wp
            END WHERE
            !
            ! nudge ice depth
            DO_2D( 0, 0, 0, 0 )
               hm_i (ji,jj) = hm_i (ji,jj) + zhicifinc(ji,jj)
            END_2D
            !
            ! seaice salinity balancing (to add)
#endif
            !
#if defined key_cice && defined key_asminc
            ! Sea-ice : CICE case. Pass ice increment tendency into CICE
            DO_2D( 0, 0, 0, 0 )
               ndaice_da(ji,jj) = seaice_bkginc(ji,jj) * zincwgt / rn_Dt
            END_2D
#endif
            !
            IF( .NOT. l_istiled .OR. ntile == nijtile )  THEN                ! Do only on the last tile
               IF ( kt == nitiaufin_r ) THEN
                  DEALLOCATE( seaice_bkginc )
               ENDIF
            ENDIF
            !
         ELSE
            !
#if defined key_cice && defined key_asminc
            DO_2D( 0, 0, 0, 0 )
               ndaice_da(ji,jj) = 0._wp        ! Sea-ice : CICE case. Zero ice increment tendency into CICE
            END_2D
#endif
            !
         ENDIF
         !                          !-----------------------------------------
      ELSEIF ( ln_asmdin ) THEN     ! Direct Initialization
         !                          !-----------------------------------------
         !
         IF ( kt == nitdin_r ) THEN
            !
            l_1st_euler = .TRUE.              ! Force Euler forward step
            !
            ! Sea-ice : SI3 case
            !
#if defined key_si3
            DO_2D( 0, 0, 0, 0 )
               zofrld (ji,jj) = 1._wp - at_i(ji,jj)
               zohicif(ji,jj) = hm_i(ji,jj)
               !
               ! Initialize the now fields the background + increment
               at_i(ji,jj) = 1. - MIN( MAX( 1.-at_i(ji,jj) - seaice_bkginc(ji,jj), 0.0_wp), 1.0_wp)
               at_i_b(ji,jj) = at_i(ji,jj)
               fr_i(ji,jj) = at_i(ji,jj)        ! adjust ice fraction
               !
               zseaicendg(ji,jj) = zofrld(ji,jj) - (1. - at_i(ji,jj))   ! find out actual sea ice nudge applied
            END_2D
            !
            ! Nudge sea ice depth to bring it up to a required minimum depth
            WHERE( zseaicendg(:,:) > 0.0_wp .AND. hm_i(A2D(0)) < zhicifmin )
               zhicifinc(:,:) = zhicifmin - hm_i(A2D(0))
            ELSEWHERE
               zhicifinc(:,:) = 0.0_wp
            END WHERE
            !
            ! nudge ice depth
            DO_2D( 0, 0, 0, 0 )
               hm_i(ji,jj) = hm_i (ji,jj) + zhicifinc(ji,jj)
            END_2D
            !
            ! seaice salinity balancing (to add)
#endif
            !
#if defined key_cice && defined key_asminc
            ! Sea-ice : CICE case. Pass ice increment tendency into CICE
            DO_2D( 0, 0, 0, 0 )
               ndaice_da(ji,jj) = seaice_bkginc(ji,jj) / rn_Dt
            END_2D
#endif
            IF( .NOT. l_istiled .OR. ntile == nijtile )  THEN                ! Do only on the last tile
               IF ( .NOT. PRESENT(kindic) ) THEN
                  DEALLOCATE( seaice_bkginc )
               END IF
            ENDIF
            !
         ELSE
            !
#if defined key_cice && defined key_asminc
            DO_2D( 0, 0, 0, 0 )
               ndaice_da(ji,jj) = 0._wp     ! Sea-ice : CICE case. Zero ice increment tendency into CICE
            END_2D
#endif
            !
         ENDIF

!#if defined defined key_si3 || 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) )*rn_Dt
!
!           ! 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,Kmm)) .and. (mld < gdepw(ji,jj,jk+1,Kmm))) THEN
!                     mld=gdepw(ji,jj,jk+1,Kmm)
!                     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
