#if defined(ROW_LAND)
#define SEA_P .true.
#define SEA_U .true.
#define SEA_V .true.
#elif defined(ROW_ALLSEA)
#define SEA_P allip(j).or.ip(i,j).ne.0
#define SEA_U alliu(j).or.iu(i,j).ne.0
#define SEA_V alliv(j).or.iv(i,j).ne.0
#else
#define SEA_P ip(i,j).ne.0
#define SEA_U iu(i,j).ne.0
#define SEA_V iv(i,j).ne.0
#endif
      module mod_incupd
      use mod_xc  ! HYCOM communication interface
!
      implicit none
!
! --- HYCOM incremental updating (for data assimilation).
! --- Stochastic anomaly forcing is also handled here for convenence.
!
      integer, save, public  :: &
       incflg,    & ! incremental update flag (0=no,-/+1=yes,-/+2=full-velocity)
                    !   -ve to read in difference archive for increment
       incstp,    & ! no. timesteps for full update (1=full insertion)
       incupf,    & ! number of days of incremental updating input
                    !   -ve to write a restart after each update completes
       incice       ! direct insertion of sea ice concentration flag
!
      logical, save, private ::  &
       lrestart  ! write a restart after each update completes
!
      integer, save, private ::  &
       ncount,    & ! increment time step counter
       ncountd      ! increment day counter
!
      real*8, save, private ::  &
       dtimeu    ! next days increment field
!
      real,    allocatable, dimension(:,:), &
               save, private :: &
        ubinc,    & !  ubaro increment
        vbinc,    & !  vbaro increment
       covinc,    & ! covice increment for direct insertion
       thkinc       ! thkice increment for direct insertion

!
      real,    allocatable, dimension(:,:,:), &
               save, private :: &
        tinc,     & !     t increment
        sinc,     & !     s increment
       dpinc,     & !    dp increment
        uinc,     & !     u increment
        vinc        !     v increment

      contains

      subroutine incupd_init(dtime0)
!
      real*8 dtime0
!
! --- subroutine used to calculate increment field for the incremental updating
! --- version: dec 2005
!
      integer i,j,k
      logical lopen
!
! --- allocate arrays
!
      allocate(  tinc(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy,kdm), &
                 sinc(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy,kdm), &
                dpinc(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy,kdm), &
                 uinc(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy,kdm), &
                 vinc(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy,kdm), &
                ubinc(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy), &
                vbinc(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy), &
               covinc(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy), &
               thkinc(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy)     )
      call mem_stat_add( 5*(idm+2*nbdy)*(jdm+2*nbdy)*kdm )
      call mem_stat_add( 4*(idm+2*nbdy)*(jdm+2*nbdy) )
!
! --- set counter to zero
!
      ncount  = 0
      ncountd = 0
      dtimeu  = 1.d-6
!
! --- restart flag
!
      lrestart = incupf.lt.0
      incupf   = abs(incupf)
!
! --- read the target fields, and initialize the "inc" arrays.
!
      if     (incflg.gt.0) then
        call incupd_read_full(dtime0)
      else  !(incflg.le.0)
        call incupd_read_diff(dtime0)
      endif
!
      return
      end subroutine incupd_init

      subroutine incupd_rd(dtime0)
!
      real*8 dtime0
!
! --- subroutine used to calculate increment field for the incremental updating
! --- version: dec 2005
!
      integer i,j,k
      logical lopen
!
      if     (ncountd.gt.incupf) then
        if     (ncountd.eq.incupf+1) then
!         should never get here, see ncountd=incupf+99 in incupd_read_*
          if (mnproc.eq.1) then
          write(lp,*) '... ended updating fields with increments ....'
          write(lp,*) 'ncountd= ',ncountd
          write(lp,*)
          endif !1st tile
          call xcsync(flush_lp)
        endif
        return
      endif
!
! --- read the target fields, and initialize the "inc" arrays.
!
      if     (incflg.gt.0) then
        call incupd_read_full(dtime0)
      else  !(incflg.le.0)
        call incupd_read_diff(dtime0)
      endif
!
      return
      end subroutine incupd_rd

      subroutine incupd(n, restrt)
      use mod_cb_arrays  ! HYCOM saved arrays
      use mod_pipe       ! HYCOM debugging interface
!
      logical restrt
      integer n
!
!**********
!*
! 1)  update hycom variables with increments.
!
! 2)  parameters:
!
!     output:
!      incremental updated model variables
!
! 4)  Ole Martin Smedstad (PSI), December 2005
!
!**********
!
      logical, parameter :: lpipe_incupd=.false.  !extra checking
!
      character utxt*12,vtxt*12
      integer   i,j,k
      real      q,utotij,vtotij
!
# include "stmt_fns.h"
!
      if     (ncountd.gt.incupf) then
!-------if (mnproc.eq.1) then
!-------write(lp,*) '... ended updating fields with increments .....'
!-------write(lp,*) 'ncountd= ',ncountd
!-------write(lp,*)
!-------endif !1st tile
!-------call xcsync(flush_lp)
        return
      endif
!
! --- update counter
!
      if     (incstp.ne.1) then
        ncount=ncount+1
        if     (ncount.eq.1 .or. ncount.eq.incstp+1) then
          q = 0.5  !half increment on 1st and incstp+1th step
        else
          q = 1.0
        endif !ncount
      else
        q = 1.0
      endif
!
      if     (ncount.gt.incstp+1) then
        if     (ncount.eq.incstp+2) then
          if (mnproc.eq.1) then
          write(lp,*) '... ended updating fields with increments ...'
          write(lp,*) 'ncount= ',ncount
          write(lp,*)
          endif !1st tile
          call xcsync(flush_lp)
        endif !ncount==incstp+2
        return
      endif !ncount>incstp+1
!
! --- ncount <= incstp+1
!
      if (mnproc.eq.1) then
      write(lp,*)
      if     (abs(incflg).eq.1) then
        write(lp,'(2a)') 'update fields with increments, ', &
                         'but not ubavg and vbavg'
      else   !abs(incflg).eq.2
        write(lp,'(2a)') 'update fields with increments, ', &
                         'including ubavg and vbavg'
      endif !incflg
      write(lp,*) '..........ncount,q= ',ncount,q
      endif !1st tile
      call xcsync(flush_lp)
!
! --- incremental update of dp (dpu, dpv).
!
!$OMP PARALLEL DO PRIVATE(j,k,i) &
!$OMP          SCHEDULE(STATIC,jblk)
      do j=1,jj
        do i=1,ii
          if (SEA_P) then
            do k=1,kk-1
! ---         dp must be non-negative
              dp(i,j,k,n) = max( dp(i,j,k,n) + q*dpinc(i,j,k), 0.0 )
! ---          p must be at or above the bottom
               p(i,j,k+1) = min( p(i,j,k) + dp(i,j,k,n), &
                                 p(i,j,kk+1) )
              dp(i,j,k,n) =      p(i,j,k+1) - p(i,j,k)
            enddo !k
! ---       layer kk always touches the bottom
            dp(i,j,kk,n) = p(i,j,kk+1) - p(i,j,kk)
          endif !ip
        enddo !i
      enddo !j
!$OMP END PARALLEL DO
!
      call dpudpv(dpu(1-nbdy,1-nbdy,1,n), &
                  dpv(1-nbdy,1-nbdy,1,n), &
                  p,depthu,depthv, 0,0)
!
      if     (lpipe_incupd) then
        do k= 1,kk
          write (utxt,'(a9,i3)') 'up dpu k=',k
          write (vtxt,'(a9,i3)') 'up dpv k=',k
          call pipe_compare_sym2(dpu(1-nbdy,1-nbdy,1,n),iu,utxt, &
                                 dpv(1-nbdy,1-nbdy,1,n),iv,vtxt)
        enddo !k
      endif !lpipe_incupd
!
! --- incremental update of the other fields.
! --- salinity from updated th&S.
! --- rebalance u and v via utotij and vtotij.
!
!$OMP PARALLEL DO PRIVATE(j,i,k,utotij,vtotij) &
!$OMP          SCHEDULE(STATIC,jblk)
      do j=1,jj
        do i=1,ii
          if (SEA_P) then
            do k=1,kk
              if     (tinc(i,j,k).ne.0.0 .or. &
                      sinc(i,j,k).ne.0.0     ) then
                temp(i,j,k,n) = temp(i,j,k,n) + q*tinc(i,j,k)
                saln(i,j,k,n) = saln(i,j,k,n) + q*sinc(i,j,k)
                th3d(i,j,k,n) = sig(temp(i,j,k,n),saln(i,j,k,n))-thbase
              endif !non-zero increment
            enddo ! k
          endif !ip
          if (SEA_U) then
            utotij = 0.0
            do k=1,kk
              u(i,j,k,n) = u(i,j,k,n) + q*uinc(i,j,k)
              utotij = utotij + u(i,j,k,n)*dpu(i,j,k,n)
            enddo ! k
            utotij=utotij/depthu(i,j)
            do k=1,kk
              u(i,j,k,n) = u(i,j,k,n) - utotij
            enddo ! k
            if     (abs(incflg).eq.2) then !update ubavg
              ubavg(i,j,n) = ubavg(i,j,n) + q*ubinc(i,j)
!             ubavg(i,j,n) = ubavg(i,j,n) + q*ubinc(i,j) + utotij
            endif !incflg==2
          endif !iu
          if (SEA_V) then
            vtotij = 0.0
            do k=1,kk
              v(i,j,k,n) = v(i,j,k,n) + q*vinc(i,j,k)
              vtotij = vtotij + v(i,j,k,n)*dpv(i,j,k,n)
            enddo ! k
            vtotij=vtotij/depthv(i,j)
            do k=1,kk
              v(i,j,k,n) = v(i,j,k,n) - vtotij
            enddo ! k
            if     (abs(incflg).eq.2) then !update vbavg
              vbavg(i,j,n) = vbavg(i,j,n) + q*vbinc(i,j)
!             vbavg(i,j,n) = vbavg(i,j,n) + q*vbinc(i,j) + vtotij
            endif !incflg==2
          endif !iv
        enddo !i
      enddo ! j
!$OMP END PARALLEL DO
!
      if     (iniflg.lt.0 .and. incice.eq.1 .and. ncount.eq.incstp) then
        if (mnproc.eq.1) then
        write(lp,*) '... direct insertion of covice and thkice'
        write(lp,*)
        endif !1st tile
        call xcsync(flush_lp)
!$OMP   PARALLEL DO PRIVATE(j,k,i) &
!$OMP            SCHEDULE(STATIC,jblk)
        do j=1,jj
          do i=1,ii
            if (SEA_P) then
              covice(i,j)=max(0.0,min(1.0,covice(i,j)+covinc(i,j)))
              thkice(i,j)=max(0.0,        thkice(i,j)+thkinc(i,j) )
            endif !ip
          enddo !li
        enddo !j
!$OMP   END PARALLEL DO
      endif !covice and thkice
!
      if (mnproc.eq.1) then
       write(lp,*) 'finished incupdate',ncount
       write(lp,*)
      endif !1st tile
      call xcsync(flush_lp)
!
      if     (lrestart .and. ncount.eq.incstp) then
        restrt = .true.
!
        if (mnproc.eq.1) then
         write(lp,*) 'incupd: set the flag to write a restart'
         write(lp,*)
        endif !1st tile
        call xcsync(flush_lp)
      endif
      return
      end subroutine incupd

      subroutine incupd_read_full(dtime)
      use mod_cb_arrays  ! HYCOM saved arrays
      use mod_za         ! HYCOM I/O interface
      use mod_pipe       ! HYCOM debugging interface
!
      real*8    dtime
!
! --- input 3-d HYCOM fields (from an archive file) on model day dtime.
! --- directly insert the input covice and thkice (if they exist).
! --- calculate the increment between the input and the initial state.
!
! --- filenames incup/incupd.iyear_iday_ihour.[ab].
! --- I/O and array I/O unit 925 used here, but not reserved.
!
      logical, parameter :: ldebug_incupd_read=.false. !usually .false.
!
      character flnm*24, cline*80, cvarin*6, cfield*8
      character ptxt*12,utxt*12,vtxt*12
      integer   i,idmtst,ios,j,jdmtst,k,l,layer,nskip
      integer   iyear,iday,ihour
      logical   nodens
      real      tincstp
!     real      sumdpi
!
      integer   nstep0
      real*8    dtime0
!
# include "stmt_fns.h"
!
      call forday(dtime, yrflag, iyear,iday,ihour)
!
      write(flnm,'("incup/incupd.",i4.4,"_",i3.3,"_",i2.2)') &
                                 iyear,iday,ihour
!
      if(dtime.ge.dtimeu) then
!
      ncountd=ncountd+1
      ncount=0
!
      if     (ncountd.gt.incupf) then
        if     (ncountd.eq.incupf+1) then
          if (mnproc.eq.1) then
          write(lp,*) '... ended updating fields with increments ...'
          write(lp,*) 'ncountd= ',ncountd
          write(lp,*)
          endif !1st tile
          call xcsync(flush_lp)
          ncountd=incupf+99 !turn off "ended" printout
        endif
        return
      endif !ncountd>incupf
!
      if (mnproc.eq.1) then
      write(lp,*) 'read (full) incremental updating ...'
      write(lp,*) 'ncountd ...',ncountd
      write (lp,*) 'incupd_read: ',flnm
      write (lp,*) '       time: ',dtime
      write (lp,*) 'iyear,iday,ihour: ',iyear,iday,ihour
      endif !1st tile
      call xcsync(flush_lp)
!
      call zaiopf(flnm//'.a','old', 925)
      if     (mnproc.eq.1) then  ! .b file from 1st tile only
        open (unit=uoff+925,file=flnm//'.b',form='formatted', &
              status='old',action='read')
!
        read(uoff+925,'(a)') cline
        read(uoff+925,'(a)') cline
        read(uoff+925,'(a)') cline
        read(uoff+925,'(a)') cline
!
        read(uoff+925,'(a)') cline  !'iversn'
        read(uoff+925,'(a)') cline  !'iexpt '
        read(uoff+925,'(a)') cline  !'yrflag'
      endif !1st tile
!
      call zagetc(cline,ios, uoff+925)
      read(cline,*) idmtst,cvarin
!     if     (mnproc.eq.1) then
!     write(lp,*) cvarin,' = ',idmtst
!     endif !1st tile
      if (cvarin.ne.'idm   ') then
        if     (mnproc.eq.1) then
        write(lp,*)
        write(lp,*) 'error in incupd_read_full - input ',cvarin, &
                              ' but should be idm   '
        write(lp,*)
        endif !1st tile
        call xcstop('(incupd_read_full)')
               stop '(incupd_read_full)'
      endif
      call zagetc(cline,ios, uoff+925)
      read(cline,*) jdmtst,cvarin
!     if     (mnproc.eq.1) then
!     write(lp,*) cvarin,' = ',jdmtst
!     endif !1st tile
      if (cvarin.ne.'jdm   ') then
        if     (mnproc.eq.1) then
        write(lp,*)
        write(lp,*) 'error in incupd_read_full - input ',cvarin, &
                              ' but should be jdm   '
        write(lp,*)
        endif !1st tile
        call xcstop('(incupd_read_full)')
               stop '(incupd_read_full)'
      endif
!
      if (idmtst.ne.itdm .or. jdmtst.ne.jtdm) then
        if     (mnproc.eq.1) then
        write(lp,*)
        write(lp,*) 'error in incupd_read_full - input idm,jdm', &
                              ' not consistent with parameters'
        write(lp,*) 'idm,jdm = ',itdm,  jtdm,  '  (dimensions.h)'
        write(lp,*) 'idm,jdm = ',idmtst,jdmtst,'  (input)'
        write(lp,*)
        endif !1st tile
        call xcstop('(incupd_read_full)')
               stop '(incupd_read_full)'
      endif
!
      if     (mnproc.eq.1) then  ! .b file from 1st tile only
        read (uoff+925,*)
      endif
!
! --- skip (most) surface fields.
!
      call zaiosk(925)
      call zagetc(cline,ios, uoff+925)
      i = index(cline,'=')
      read(cline(i+1:),*) nstep0,dtime0,layer
      if     (mnproc.eq.1) then
        write(lp,*) 'dtime0= ',dtime0
      endif
      if (dtime0.ne.dtime) then
        if     (mnproc.eq.1) then
        write(lp,*)
        write(lp,*) 'error in incupd_read_full - input ',dtime0, &
                            ' but dtime should be ',dtime
        write(lp,*)
        endif !1st tile
        call xcstop('(incupd_read_full)')
               stop '(incupd_read_full)'
      endif
      nodens = layer.ne.0  !new or original archive type
      if     (nodens .and. layer.ne.sigver) then
        if     (mnproc.eq.1) then
        write(lp,*)
        write(lp,*) 'error in incupd_read_full - input ',layer, &
                           ' sigver but should be ',sigver
        write(lp,*)
        endif !1st tile
        call xcstop('(incupd_read_full)')
               stop '(incupd_read_full)'
      endif
!
! assumes that there is a new incremental updating file once a day
! for "incupf" days, see blkdat.input
!
      dtimeu=dtime0+1.d0
!
      if     (mnproc.eq.1) then
        write(lp,*)
        write(lp,*) 'dtime, dtime0, dtimeu = ',dtime, &
                     dtime0, dtimeu
        write(lp,*)
      endif !1st tile
      call xcsync(flush_lp)
!
      if     (nodens) then
        do i= 2,6
          if     (mnproc.eq.1) then  ! .b file from 1st tile only
            read (uoff+925,*)
          endif
          call zaiosk(925)
        enddo
      else
        do i= 2,11
          if     (mnproc.eq.1) then  ! .b file from 1st tile only
            read (uoff+925,*)
          endif
          call zaiosk(925)
        enddo
      endif
!
      call rd_archive(ubinc, cfield,layer, 925)  !u_btrop or covice or mix_dpth or kemix
      if     (cfield.eq.'mix_dpth' .or. cfield.eq.'kemix') then
! ---   archive contains 'steric  '
        call rd_archive(ubinc, cfield,layer, 925)  !u_btrop or covice
      endif
      if     (mnproc.eq.1) then
      write(lp,'(2a)') "surface: ",cfield
      endif
      call xcsync(flush_lp)
      if     (cfield.eq.'covice  ') then
!
! ---   directly insert covice and thkice.
!
        call rd_archive(util5, cfield,layer, 925)  !thkice
        if     (mnproc.eq.1) then
        write(lp,'(2a)') "surface: ",cfield
        endif
        call xcsync(flush_lp)
!$OMP   PARALLEL DO PRIVATE(j,k,i) &
!$OMP            SCHEDULE(STATIC,jblk)
        do j=1,jj
          do i=1,ii
            if (SEA_P) then
              covice(i,j)=ubinc(i,j)
              thkice(i,j)=util5(i,j)
            endif !ip
          enddo !li
        enddo !j
!$OMP   END PARALLEL DO
        call zaiosk(925)  !temice
        if     (mnproc.eq.1) then  ! .b file from 1st tile only
          read (uoff+925,*)
        endif
        call rd_archive(ubinc, cfield,layer, 925)
        if     (mnproc.eq.1) then
        write(lp,'(2a)') "surface: ",cfield
        endif
        call xcsync(flush_lp)
        incice =  1  !have     input covice, don't direct insert si_c
      else
        incice = -1  !have not input covice, might direct insert si_c
      endif
      call rd_archive(vbinc, cfield,layer, 925)
      if     (mnproc.eq.1) then
      write(lp,'(2a)') "surface: ",cfield
      endif
      call xcsync(flush_lp)
!
           if     (mnproc.eq.1) then
           write (lp,*) 'start 3-D archive file read'
           endif
           call xcsync(flush_lp)
!
! --- 3-d fields.
!
      nskip = 0
      do k=1,kk
        call rd_archive(uinc(1-nbdy,1-nbdy,k), cfield,layer, 925)
        if     (cfield.ne.'u-vel.  ' .and. k.ne.2) then
          if     (mnproc.eq.1) then
          write(lp,'(/ a / a,a /)') cfield, &
                 'error in incupd_read_full - expected ','u-vel.  '
          endif !1st tile
          call xcstop('(incupd_read_full)')
                 stop '(incupd_read_full)'
        elseif (cfield.ne.'u-vel.  ') then !k==2
!
! ---     count "tracer" fields (to be skipped)
!
          if     (mnproc.eq.1) then
          write(lp,'(2a)') "counting tracers: ",cfield
          endif
          do nskip= 2,99
            call rd_archive(uinc(1-nbdy,1-nbdy,k), cfield,layer, 925)
            if     (mnproc.eq.1) then
            write(lp,'(2a)') "counting tracers: ",cfield
            endif
            if     (cfield.eq.'u-vel.  ') then
              exit
            endif
          enddo !nskip
          nskip = nskip - 1
          write(lp,'(a,i3)') "nskip =",nskip
        endif
        call rd_archive(vinc(1-nbdy,1-nbdy,k), cfield,layer, 925)
        if     (cfield.ne.'v-vel.  ') then
          if     (mnproc.eq.1) then
          write(lp,'(/ a / a,a /)') cfield, &
                 'error in incupd_read_full - expected ','v-vel.  '
          endif !1st tile
          call xcstop('(incupd_read_full)')
                 stop '(incupd_read_full)'
        endif
!          if     (mnproc.eq.1) then
!            write (lp,*) 'read v-vel archive file'
!          endif
        call rd_archive(dpinc(1-nbdy,1-nbdy,k), cfield,layer, 925)
        if     (cfield.ne.'thknss  ') then
          if     (mnproc.eq.1) then
          write(lp,'(/ a / a,a /)') cfield, &
                   'error in incupd_read_full - expected ','thknss  '
          endif !1st tile
          call xcstop('(incupd_read_full)')
                 stop '(incupd_read_full)'
        endif
!          if     (mnproc.eq.1) then
!            write (lp,*) 'read dpinc archive file'
!          endif
        call rd_archive(tinc(1-nbdy,1-nbdy,k), cfield,layer, 925)
        if     (cfield.ne.'temp    ') then
          if     (mnproc.eq.1) then
          write(lp,'(/ a / a,a /)') cfield, &
                 'error in incupd_read_full - expected ','temp    '
          endif !1st tile
          call xcstop('(incupd_read_full)')
                 stop '(incupd_read_full)'
        endif
        call rd_archive(sinc(1-nbdy,1-nbdy,k), cfield,layer, 925)
        if     (cfield.ne.'salin   ') then
          if     (mnproc.eq.1) then
          write(lp,'(/ a / a,a /)') cfield, &
                 'error in incupd_read_full - expected ','salin   '
          endif !1st tile
          call xcstop('(incupd_read_full)')
                 stop '(incupd_read_full)'
        endif
        if     (.not. nodens) then
! ---     skip density
          if     (mnproc.eq.1) then  ! .b file from 1st tile only
            read (uoff+925,*)
          endif
          call zaiosk(925)
        endif !nodens:else
!
! ---   skip (nskip) tracers
!
        do l= 1,nskip
          if     (mnproc.eq.1) then  ! .b file from 1st tile only
            read (uoff+925,*)
          endif
          call zaiosk(925)
        enddo !l
      enddo !k
!
      call xctilr(dpinc,1,kk, 1,1, halo_ps)
!
      if     (mnproc.eq.1) then  ! .b file from 1st tile only
      close( unit=uoff+925)
      endif
      call zaiocl(925)
!
! --- calculate increments
! --- the "inc" reads, above, are full HYCOM fields (not increments yet).
!
      if(incstp.eq.1) then
        tincstp=1.0
        if     (mnproc.eq.1) then
        write(lp,*)
        write(lp,*) 'tincstp=1.0        ',tincstp,incstp
        endif
      else
! ---   allow for LeapFrog and Asselin filter
        if     (mod(incstp,2).eq.1) then
          incstp = incstp-1  !must be even
        endif
        tincstp=(2.0/sqrt(1.0-ra2fac))/real(incstp)
        if     (mnproc.eq.1) then
        write(lp,*)
        write(lp,'(a,f6.4,a,f12.7,i5)') &
          'tincstp=',2.0/sqrt(1.0-ra2fac),'/incstp ', &
           tincstp,incstp
        endif
      endif !incstp
!
      if (mnproc.eq.1) then
      write(lp,*)
      write(lp,*) 'calculate t,s,u,v and dp increments'
      endif !1st tile
      call xcsync(flush_lp)
!
!     if     (iutest.gt.0 .and. jutest.gt.0) then
!       write(lp,*) '*',' iutest= ',iutest+i0,' jutest= ',jutest+j0,' *'
!       write(lp,*) '*********** dpinc input ************'
!               sumdpi=0.0
!             write(lp,'(a)')
!    &                'k,dp1,dp2,dpinc='
!               do k= 1,kk
!                sumdpi=sumdpi+dpinc(iutest,jutest,k)
!                   write(lp,'(a,i3,3f20.5)')
!    &                'k= ',
!    &                 k,dp(iutest,jutest,k,1)*qonem,
!    &                 dp(iutest,jutest,k,2)*qonem,
!    &                 dpinc(iutest,jutest,k)*qonem
!                   call flush(lp)
!               enddo !k
!            write(lp,*) 'sumdpi= ', sumdpi*qonem
!            call flush(lp)
!     endif
!
!$OMP PARALLEL DO PRIVATE(j,k,i) &
!$OMP          SCHEDULE(STATIC,jblk)
      do j=1,jj
          do i=1,ii
            if (SEA_U) then
            ubinc(i,j)=(ubinc(i,j) - ubavg(i,j,1))*tincstp
            do k=1,kk
!             use an approximate 2*dpu
              if     (dpinc(i,j,k)+dpinc(i-1,j,k).gt.2.0*onemm) then
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!! MODIFICACAO FEITA PELO VLADIMIR/LUANA PARA NAO ATUALIZAR !!!
!!!!!!!       COM O INCUP AS VELOCIDADES BAROCLINICAS        !!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!                uinc(i,j,k)=(uinc(i,j,k) - u(i,j,k,1))*tincstp
                uinc(i,j,k)=0.0
              else
                uinc(i,j,k)=0.0  !thin target layer
              endif
            enddo !k
          endif !iu
          if (SEA_V) then
            vbinc(i,j)=(vbinc(i,j) - vbavg(i,j,1))*tincstp
            do k=1,kk
!             use an approximate 2*dpv
              if     (dpinc(i,j,k)+dpinc(i,j-1,k).gt.2.0*onemm) then
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!! MODIFICACAO FEITA PELO VLADIMIR/LUANA PARA NAO ATUALIZAR !!!
!!!!!!!       COM O INCUP AS VELOCIDADES BAROCLINICAS        !!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!                vinc(i,j,k)=(vinc(i,j,k) - v(i,j,k,1))*tincstp
                vinc(i,j,k)=0.0
              else  
                vinc(i,j,k)=0.0  !thin target layer
              endif
            enddo !k
          endif !iv
          if (SEA_P) then
            do k=1,kk
              if     (dpinc(i,j,k).gt.onemm) then
                sinc(i,j,k)=(sinc(i,j,k) - saln(i,j,k,1))*tincstp
                tinc(i,j,k)=(tinc(i,j,k) - temp(i,j,k,1))*tincstp
              else
                tinc(i,j,k)=0.0  !thin target layer
                sinc(i,j,k)=0.0  !thin target layer
              endif
              dpinc(i,j,k)=(dpinc(i,j,k) - dp(i,j,k,1))*tincstp
            enddo !k
          endif !ip
        enddo !li
      enddo !j
!$OMP END PARALLEL DO
!
      call xctilr(dpinc,1,kk, 1,1, halo_ps)
!
      if     (ldebug_incupd_read) then
         call pipe_compare_sym2(ubinc,iu,'incupd:ubinc', &
                                vbinc,iv,'incupd:vbinc')
         do k= 1,kk
           write (utxt,'(a9,i3)') '  uinc k=',k
           write (vtxt,'(a9,i3)') '  vinc k=',k
           call pipe_compare_sym2(uinc(1-nbdy,1-nbdy,k),iu,utxt, &
                                  vinc(1-nbdy,1-nbdy,k),iv,vtxt)
           write (ptxt,'(a9,i3)') ' dpinc k=',k
           call pipe_compare_sym1(dpinc(1-nbdy,1-nbdy,k),ip,ptxt)
           write (ptxt,'(a9,i3)') '  tinc k=',k
           call pipe_compare_sym1( tinc(1-nbdy,1-nbdy,k),ip,ptxt)
           write (ptxt,'(a9,i3)') '  sinc k=',k
           call pipe_compare_sym1( sinc(1-nbdy,1-nbdy,k),ip,ptxt)
         enddo !k
       endif !ldebug_incupd_read
!
!     if     (iutest.gt.0 .and. jutest.gt.0) then
!       write(lp,*) '*',' iutest= ',iutest+i0,' jutest= ',jutest+j0,' *'
!       write(lp,*) '*********** dpinc out ************'
!             write(lp,'(a)')
!    &                'k,dp1,dp2,dpinc='
!               sumdpi=0.0
!               do k= 1,kk
!                sumdpi=sumdpi+dpinc(iutest,jutest,k)
!                   write(lp,'(a,i3,3f20.5)')
!    &                'k= ',
!    &                 k,dp(iutest,jutest,k,1)*qonem,
!    &                 dp(iutest,jutest,k,2)*qonem,
!    &                 dpinc(iutest,jutest,k)*qonem
!                   call flush(lp)
!               enddo !k
!            write(lp,*) 'inc sumdpi= ', sumdpi*qonem
!            call flush(lp)
!     endif
!
      if (mnproc.eq.1) then
       write(lp,*) '... finished reading incupd',dtime,dtime0
      endif !1st tile
      call xcsync(flush_lp)
!
      endif ! dtime
!
      return
      end subroutine incupd_read_full

      subroutine incupd_read_diff(dtime)
      use mod_cb_arrays  ! HYCOM saved arrays
      use mod_za         ! HYCOM I/O interface
      use mod_pipe       ! HYCOM debugging interface
!
      real*8    dtime
!
! --- input 3-d HYCOM increment fields (from a difference archive file)
! --- on model day dtime or dtime + insertion period.
! --- read in covice and thkice increments, for later direct insertion.
!
! --- filenames incup/incupd.iyear_iday_ihour.[ab].
! --- I/O and array I/O unit 925 used here, but not reserved.
!
      logical, parameter :: ldebug_incupd_read=.false. !usually .false.
!
      character flnm*24, cline*80, cvarin*6, cfield*8
      character ptxt*12,utxt*12,vtxt*12
      integer   i,idmtst,ios,j,jdmtst,k,l,layer,nskip
      integer   iyear,iday,ihour
      real      tincstp
!     real      sumdpi
!
      logical   lexist
      integer   nstep0
      real*8    dtimei,dtime0
!
# include "stmt_fns.h"
!
      dtimei = dtime
!
      call forday(dtimei, yrflag, iyear,iday,ihour)
!
      write(flnm,'("incup/incupd.",i4.4,"_",i3.3,"_",i2.2)') &
                                 iyear,iday,ihour
!
      inquire(file=flnm//'.a',exist=lexist)
!
! --- is input at end of insertion period
!
      if     (.not. lexist) then
        dtimei = dtime + incstp*(baclin/86400.d0)
!
        call forday(dtimei, yrflag, iyear,iday,ihour)
!
        write(flnm,'("incup/incupd.",i4.4,"_",i3.3,"_",i2.2)') &
                                   iyear,iday,ihour
      endif
!
      if(dtime.ge.dtimeu) then
!
      ncountd=ncountd+1
      ncount=0
!
      if     (ncountd.gt.incupf) then
        if     (ncountd.eq.incupf+1) then
          if (mnproc.eq.1) then
          write(lp,*) '... ended updating fields with increments ...'
          write(lp,*) 'ncountd= ',ncountd
          write(lp,*)
          endif !1st tile
          call xcsync(flush_lp)
          ncountd=incupf+99 !turn off "ended" printout
        endif
        return
      endif !ncountd>incupf
!
      if (mnproc.eq.1) then
      write(lp,*) 'read (diff) incremental updating ...'
      write(lp,*) 'ncountd ...',ncountd
      write (lp,*) 'incupd_read: ',flnm
      if     (dtime.eq.dtimei) then
        write (lp,*) '       time: ',dtime
      else
        write (lp,*) '       time: ',dtime
        write (lp,*) '    time_in: ',dtimei
      endif
      write (lp,*) 'iyear,iday,ihour: ',iyear,iday,ihour
      endif !1st tile
      call xcsync(flush_lp)
!
      call zaiopf(flnm//'.a','old', 925)
      if     (mnproc.eq.1) then  ! .b file from 1st tile only
        open (unit=uoff+925,file=flnm//'.b',form='formatted', &
              status='old',action='read')
!
        read(uoff+925,'(a)') cline
        read(uoff+925,'(a)') cline
        read(uoff+925,'(a)') cline
        read(uoff+925,'(a)') cline
!
        read(uoff+925,'(a)') cline  !'iversn'
        read(uoff+925,'(a)') cline  !'jexpt '
        read(uoff+925,'(a)') cline  !'iexpt '
        read(uoff+925,'(a)') cline  !'yrflag'
      endif !1st tile
!
      call zagetc(cline,ios, uoff+925)
      read(cline,*) idmtst,cvarin
!     if     (mnproc.eq.1) then
!     write(lp,*) cvarin,' = ',idmtst
!     endif !1st tile
      if (cvarin.ne.'idm   ') then
        if     (mnproc.eq.1) then
        write(lp,*)
        write(lp,*) 'error in incupd_read_diff - input ',cvarin, &
                              ' but should be idm   '
        write(lp,*)
        endif !1st tile
        call xcstop('(incupd_read_diff)')
               stop '(incupd_read_diff)'
      endif
      call zagetc(cline,ios, uoff+925)
      read(cline,*) jdmtst,cvarin
!     if     (mnproc.eq.1) then
!     write(lp,*) cvarin,' = ',jdmtst
!     endif !1st tile
      if (cvarin.ne.'jdm   ') then
        if     (mnproc.eq.1) then
        write(lp,*)
        write(lp,*) 'error in incupd_read_diff - input ',cvarin, &
                              ' but should be jdm   '
        write(lp,*)
        endif !1st tile
        call xcstop('(incupd_read_diff)')
               stop '(incupd_read_diff)'
      endif
!
      if (idmtst.ne.itdm .or. jdmtst.ne.jtdm) then
        if     (mnproc.eq.1) then
        write(lp,*)
        write(lp,*) 'error in incupd_read_diff - input idm,jdm', &
                              ' not consistent with parameters'
        write(lp,*) 'idm,jdm = ',itdm,  jtdm,  '  (dimensions.h)'
        write(lp,*) 'idm,jdm = ',idmtst,jdmtst,'  (input)'
        write(lp,*)
        endif !1st tile
        call xcstop('(incupd_read_diff)')
               stop '(incupd_read_diff)'
      endif
!
! --- check for a difference archive
!
      call zagetc(cline,ios, uoff+925)
      i = index(cline,'diff day')
      if     (i.eq.0) then
        if     (mnproc.eq.1) then
        write(lp,*)
        write(lp,*) 'error in incupd_read_diff - ', &
                              'not a difference archive'
        write(lp,*) 'header = '
        write(lp,*)  trim(cline)
        write(lp,*)
        endif !1st tile
        call xcstop('(incupd_read_diff)')
               stop '(incupd_read_diff)'
      endif
!
! --- skip (most) surface fields.
!
      call zaiosk(925)
      call zagetc(cline,ios, uoff+925)
      i = index(cline,'=')
      read(cline(i+1:),*) nstep0,dtime0,layer
      if     (mnproc.eq.1) then
        write(lp,*) 'dtime0= ',dtime0
      endif
      if (dtime0.ne.dtimei) then
        if     (mnproc.eq.1) then
        write(lp,*)
        write(lp,*) 'error in incupd_read_diff - input ',dtime0, &
                                ' but dtimei should be ',dtimei
        write(lp,*)
        endif !1st tile
        call xcstop('(incupd_read_diff)')
               stop '(incupd_read_diff)'
      endif
!
! assumes that there is a new incremental updating file once a day
! for "incupf" days, see blkdat.input
!
      dtimeu=dtime0+1.d0
!
      if     (mnproc.eq.1) then
        write(lp,*)
        write(lp,*) 'dtime, dtime0, dtimeu = ',dtime, &
                     dtime0, dtimeu
        write(lp,*)
      endif !1st tile
      call xcsync(flush_lp)
!
      do i= 2,12
        if     (mnproc.eq.1) then  ! .b file from 1st tile only
          read (uoff+925,*)
        endif
        call zaiosk(925)
      enddo
!
      call rd_archive(ubinc, cfield,layer, 925)  !u_btrop or covice or kemix
      if     (cfield.eq.'kemix') then
! ---   archive contains 'steric  '
        call rd_archive(ubinc, cfield,layer, 925)  !u_btrop or covice
      endif
      if     (mnproc.eq.1) then
      write(lp,'(2a)') "surface: ",cfield
      endif
      call xcsync(flush_lp)
      if     (cfield.eq.'covice  ') then
        covinc(:,:) = ubinc(:,:)
        call rd_archive(thkinc, cfield,layer, 925)  !thkice
        if     (mnproc.eq.1) then
        write(lp,'(2a)') "surface: ",cfield
        endif
        call xcsync(flush_lp)
        call zaiosk(925)  !temice
        if     (mnproc.eq.1) then  ! .b file from 1st tile only
          read (uoff+925,*)
        endif
        call rd_archive(ubinc, cfield,layer, 925)
        if     (mnproc.eq.1) then
        write(lp,'(2a)') "surface: ",cfield
        endif
        call xcsync(flush_lp)
        incice =  1  !have     input covice, don't direct insert si_c
      else
        incice = -1  !have not input covice, might direct insert si_c
      endif
      call rd_archive(vbinc, cfield,layer, 925)
      if     (mnproc.eq.1) then
      write(lp,'(2a)') "surface: ",cfield
      endif
      call xcsync(flush_lp)
! --- kebtrop
      if     (mnproc.eq.1) then  ! .b file from 1st tile only
        read (uoff+925,*)
      endif
      call zaiosk(925)
!
           if     (mnproc.eq.1) then
           write (lp,*) 'start 3-D archive file read'
           endif
           call xcsync(flush_lp)
!
! --- 3-d fields.
!
      nskip = 0
      do k=1,kk
        call rd_archive(uinc(1-nbdy,1-nbdy,k), cfield,layer, 925)
        if     (cfield.ne.'u-vel.  ' .and. k.ne.2) then
          if     (mnproc.eq.1) then
          write(lp,'(/ a / a,a /)') cfield, &
                 'error in incupd_read_diff - expected ','u-vel.  '
          endif !1st tile
          call xcstop('(incupd_read_diff)')
                 stop '(incupd_read_diff)'
        elseif (cfield.ne.'u-vel.  ') then !k==2
!
! ---     count "tracer" fields (to be skipped)
!
          if     (mnproc.eq.1) then
          write(lp,'(2a)') "counting tracers: ",cfield
          endif
          do nskip= 2,99
            call rd_archive(uinc(1-nbdy,1-nbdy,k), cfield,layer, 925)
            if     (mnproc.eq.1) then
            write(lp,'(2a)') "counting tracers: ",cfield
            endif
            if     (cfield.eq.'u-vel.  ') then
              exit
            endif
          enddo !nskip
          nskip = nskip - 1
          write(lp,'(a,i3)') "nskip =",nskip
        endif
        call rd_archive(vinc(1-nbdy,1-nbdy,k), cfield,layer, 925)
        if     (cfield.ne.'v-vel.  ') then
          if     (mnproc.eq.1) then
          write(lp,'(/ a / a,a /)') cfield, &
                 'error in incupd_read_diff - expected ','v-vel.  '
          endif !1st tile
          call xcstop('(incupd_read_diff)')
                 stop '(incupd_read_diff)'
        endif
!          if     (mnproc.eq.1) then
!            write (lp,*) 'read v-vel archive file'
!          endif
! ---   skip k.e.
        if     (mnproc.eq.1) then  ! .b file from 1st tile only
          read (uoff+925,*)
        endif
        call zaiosk(925)
! ---   skip mnthknss
        if     (mnproc.eq.1) then  ! .b file from 1st tile only
          read (uoff+925,*)
        endif
        call zaiosk(925)
!
        call rd_archive(dpinc(1-nbdy,1-nbdy,k), cfield,layer, 925)
        if     (cfield.ne.'thknss  ') then
          if     (mnproc.eq.1) then
          write(lp,'(/ a / a,a /)') cfield, &
                   'error in incupd_read_diff - expected ','thknss  '
          endif !1st tile
          call xcstop('(incupd_read_diff)')
                 stop '(incupd_read_diff)'
        endif
!          if     (mnproc.eq.1) then
!            write (lp,*) 'read dpinc archive file'
!          endif
        call rd_archive(tinc(1-nbdy,1-nbdy,k), cfield,layer, 925)
        if     (cfield.ne.'temp    ') then
          if     (mnproc.eq.1) then
          write(lp,'(/ a / a,a /)') cfield, &
                 'error in incupd_read_diff - expected ','temp    '
          endif !1st tile
          call xcstop('(incupd_read_diff)')
                 stop '(incupd_read_diff)'
        endif
        call rd_archive(sinc(1-nbdy,1-nbdy,k), cfield,layer, 925)
        if     (cfield.ne.'salin   ') then
          if     (mnproc.eq.1) then
          write(lp,'(/ a / a,a /)') cfield, &
                 'error in incupd_read_diff - expected ','salin   '
          endif !1st tile
          call xcstop('(incupd_read_diff)')
                 stop '(incupd_read_diff)'
        endif
! ---   skip density
        if     (mnproc.eq.1) then  ! .b file from 1st tile only
          read (uoff+925,*)
        endif
        call zaiosk(925)
!
! ---   skip (nskip) tracers
!
        do l= 1,nskip
          if     (mnproc.eq.1) then  ! .b file from 1st tile only
            read (uoff+925,*)
          endif
          call zaiosk(925)
        enddo !l
      enddo !k
!
      call xctilr(dpinc,1,kk, 1,1, halo_ps)
!
      if     (mnproc.eq.1) then  ! .b file from 1st tile only
      close( unit=uoff+925)
      endif
      call zaiocl(925)
!
! --- calculate increments
! --- the "inc" reads, above, are HYCOM diff fields.
!
      if(incstp.eq.1) then
        tincstp=1.0
        if     (mnproc.eq.1) then
        write(lp,*)
        write(lp,*) 'tincstp=1.0        ',tincstp,incstp
        endif
      else
! ---   allow for LeapFrog and Asselin filter
        if     (mod(incstp,2).eq.1) then
          incstp = incstp-1  !must be even
        endif
        tincstp=(2.0/sqrt(1.0-ra2fac))/real(incstp)
        if     (mnproc.eq.1) then
        write(lp,*)
        write(lp,'(a,f6.4,a,f12.7,i5)') &
          'tincstp=',2.0/sqrt(1.0-ra2fac),'/incstp ', &
           tincstp,incstp
        endif
      endif !incstp
!
      if (mnproc.eq.1) then
      write(lp,*)
      write(lp,*) 'calculate t,s,u,v and dp increments'
      endif !1st tile
      call xcsync(flush_lp)
!
!     if     (iutest.gt.0 .and. jutest.gt.0) then
!       write(lp,*) '*',' iutest= ',iutest+i0,' jutest= ',jutest+j0,' *'
!       write(lp,*) '*********** dpinc input ************'
!               sumdpi=0.0
!             write(lp,'(a)')
!    &                'k,dp1,dp2,dpinc='
!               do k= 1,kk
!                sumdpi=sumdpi+dpinc(iutest,jutest,k)
!                   write(lp,'(a,i3,3f20.5)')
!    &                'k= ',
!    &                 k,dp(iutest,jutest,k,1)*qonem,
!    &                 dp(iutest,jutest,k,2)*qonem,
!    &                 dpinc(iutest,jutest,k)*qonem
!                   call flush(lp)
!               enddo !k
!            write(lp,*) 'sumdpi= ', sumdpi*qonem
!            call flush(lp)
!     endif
!
! --- use background thickness fields for thinness detection
!
!$OMP PARALLEL DO PRIVATE(j,k,i) &
!$OMP          SCHEDULE(STATIC,jblk)
      do j=1,jj
          do i=1,ii
            if (SEA_U) then
            ubinc(i,j)=ubinc(i,j)*tincstp
            do k=1,kk
              if     (dpu(i,j,k,1).gt.onemm) then
                uinc(i,j,k)=uinc(i,j,k)*tincstp
              else
                uinc(i,j,k)=0.0  !thin target layer
              endif
            enddo !k
          endif !iu
          if (SEA_V) then
            vbinc(i,j)=vbinc(i,j)*tincstp
            do k=1,kk
              if     (dpv(i,j,k,1).gt.onemm) then
                vinc(i,j,k)=vinc(i,j,k)*tincstp
              else  
                vinc(i,j,k)=0.0  !thin target layer
              endif
            enddo !k
          endif !iv
          if (SEA_P) then
            do k=1,kk
              if     (dp(i,j,k,1).gt.onemm) then
                sinc(i,j,k)=sinc(i,j,k)*tincstp
                tinc(i,j,k)=tinc(i,j,k)*tincstp
              else
                tinc(i,j,k)=0.0  !thin target layer
                sinc(i,j,k)=0.0  !thin target layer
              endif
              dpinc(i,j,k)=dpinc(i,j,k)*tincstp
            enddo !k
          endif !ip
        enddo !li
      enddo !j
!$OMP END PARALLEL DO
!
      call xctilr(dpinc,1,kk, 1,1, halo_ps)
!
      if     (ldebug_incupd_read) then
         call pipe_compare_sym2(ubinc,iu,'incupd:ubinc', &
                                vbinc,iv,'incupd:vbinc')
         do k= 1,kk
           write (utxt,'(a9,i3)') '  uinc k=',k
           write (vtxt,'(a9,i3)') '  vinc k=',k
           call pipe_compare_sym2(uinc(1-nbdy,1-nbdy,k),iu,utxt, &
                                  vinc(1-nbdy,1-nbdy,k),iv,vtxt)
           write (ptxt,'(a9,i3)') ' dpinc k=',k
           call pipe_compare_sym1(dpinc(1-nbdy,1-nbdy,k),ip,ptxt)
           write (ptxt,'(a9,i3)') '  tinc k=',k
           call pipe_compare_sym1( tinc(1-nbdy,1-nbdy,k),ip,ptxt)
           write (ptxt,'(a9,i3)') '  sinc k=',k
           call pipe_compare_sym1( sinc(1-nbdy,1-nbdy,k),ip,ptxt)
         enddo !k
       endif !ldebug_incupd_read
!
!     if     (iutest.gt.0 .and. jutest.gt.0) then
!       write(lp,*) '*',' iutest= ',iutest+i0,' jutest= ',jutest+j0,' *'
!       write(lp,*) '*********** dpinc out ************'
!             write(lp,'(a)')
!    &                'k,dp1,dp2,dpinc='
!               sumdpi=0.0
!               do k= 1,kk
!                sumdpi=sumdpi+dpinc(iutest,jutest,k)
!                   write(lp,'(a,i3,3f20.5)')
!    &                'k= ',
!    &                 k,dp(iutest,jutest,k,1)*qonem,
!    &                 dp(iutest,jutest,k,2)*qonem,
!    &                 dpinc(iutest,jutest,k)*qonem
!                   call flush(lp)
!               enddo !k
!            write(lp,*) 'inc sumdpi= ', sumdpi*qonem
!            call flush(lp)
!     endif
!
      if (mnproc.eq.1) then
       write(lp,*) '... finished reading incupd',dtime,dtime0
      endif !1st tile
      call xcsync(flush_lp)
!
      endif ! dtime
!
      return
      end subroutine incupd_read_diff

      subroutine incupd_si_c(dtime)
      use mod_cb_arrays  ! HYCOM saved arrays
      use mod_za         ! HYCOM I/O interface
      use mod_pipe       ! HYCOM debugging interface
!
      real*8    dtime
!
! --- directly insert si_c into covice and thkice.
!
      integer   i,j,l
!
      if     (incice.eq.-1) then
        incice = 0  !will have directly inserted si_c
!
! ---   directly insert covice and thkice.
!
!$OMP   PARALLEL DO PRIVATE(j,i) &
!$OMP            SCHEDULE(STATIC,jblk)
        do j=1-nbdy,jj+nbdy
          do i=1,ii
            if (SEA_P) then
              covice(i,j)=si_c(i,j)
              thkice(i,j)=covice(i,j)*hicemn
            endif !ip
          enddo !i
        enddo !j
!$OMP   END PARALLEL DO
!
        if (mnproc.eq.1) then
         write(lp,*) '... finished inserting si_c',dtime
        endif !1st tile
        call xcsync(flush_lp)
      endif !incice
!
      return
      end subroutine incupd_si_c
 
      subroutine stfupd(n)
      use mod_cb_arrays  ! HYCOM saved arrays
      use mod_pipe       ! HYCOM debugging interface
!
      integer n
!
!**********
!*
! 1)  update hycom variables with stochastic anomalies.
!
! 2)  parameters:
!
!     output:
!      updated model variables
!
! 4)  Alan Wallcaft, NRL, September 2017.
!     Based on incupd.
!
!**********
!
      logical, parameter :: lpipe_stfupd=.false.  !extra checking
!
      integer   i,j,k
      real      q,q_s,q_t,zij,utotij,vtotij
!
# include "stmt_fns.h"
!
      if     (stfflg.eq.0) then
        return
      endif
!
! --- incremental update of T&S&vel.
! --- rebalance u and v via utotij and vtotij.
!
!$OMP PARALLEL DO PRIVATE(j,i,k,utotij,vtotij) &
!$OMP          SCHEDULE(STATIC,jblk)
      do j=1,jj
        do i=1,ii
          if (SEA_P) then
            zij = -0.5*dp(i,j,1,n)*qonem  !layer 1 correction
            do k=1,kk
              zij = zij + 0.5*qonem*(dp(i,j,max(k-1,1),n)+dp(i,j,k,n))  !layer k center
              q_t = (delt1/3600.0)*exp(-zij/stfrdt)
              q_s = (delt1/3600.0)*exp(-zij/stfrds)
              temp(i,j,k,n) = temp(i,j,k,n) + q_t*(stoc_t(i,j,l0)*w0+ &
                                                   stoc_t(i,j,l1)*w1 )
              saln(i,j,k,n) = saln(i,j,k,n) + q_s*(stoc_s(i,j,l0)*w0+ &
                                                   stoc_s(i,j,l1)*w1 )
              th3d(i,j,k,n) = sig(temp(i,j,k,n),saln(i,j,k,n))-thbase
            enddo ! k
          endif !ip
          if     (stfflg.eq.2) then
          if (SEA_U) then
            utotij = 0.0
            zij = -0.5*dpu(i,j,1,n)*qonem  !layer 1 correction
            do k=1,kk
              zij = zij + 0.5*qonem*(dpu(i,j,max(k-1,1),n)+dpu(i,j,k,n))  !layer k center
              q   = (delt1/3600.0)*exp(-zij/stfrdv)
              u(i,j,k,n) = u(i,j,k,n) + q*(stoc_u(i,j,l0)*w0+ &
                                           stoc_u(i,j,l1)*w1 )
            enddo ! k
            utotij=utotij/depthu(i,j)
            do k=1,kk
              u(i,j,k,n) =   u(i,j,k,n) - utotij
            enddo ! k
            ubavg(i,j,n) = ubavg(i,j,n) + utotij
          endif !iu
          if (SEA_V) then
            vtotij = 0.0
            zij = -0.5*dpv(i,j,1,n)*qonem  !layer 1 correction
            do k=1,kk
              zij = zij + 0.5*qonem*(dpv(i,j,max(k-1,1),n)+dpv(i,j,k,n))  !layer k center
              q   = (delt1/3600.0)*exp(-zij/stfrdv)
              v(i,j,k,n) = v(i,j,k,n) + q*(stoc_v(i,j,l0)*w0+ &
                                           stoc_v(i,j,l1)*w1 )
              vtotij = vtotij + v(i,j,k,n)*dpv(i,j,k,n)
            enddo ! k
            vtotij=vtotij/depthv(i,j)
            do k=1,kk
              v(i,j,k,n) =   v(i,j,k,n) - vtotij
            enddo ! k
            vbavg(i,j,n) = vbavg(i,j,n) + vtotij
          endif !iv
          endif !stfflg==2
        enddo !i
      enddo ! j
!$OMP END PARALLEL DO
!
      return
      end subroutine stfupd
!
      end module mod_incupd
!
!
!> Revision history:
!>
!> Feb  2006 - 1st module version
!> May  2006 - changed to read multiple increment files
!> Jul  2011 - thin layer is now 1mm (no longer 1m)
!> Jul  2011 - replace thinc with sinc
!> Nov  2012 - bugfix: added xctilr(dpinc to update halo
!> Dec  2012 - cleaned up printing to .log file
!> Apr  2012 - added incice and incupd_si_c
!> May  2014 - use land/sea masks (e.g. ip) to skip land
!> Nov  2015 - added iniflg<0 for reading a difference archive
!> Dec  2015 - allow for Asselin filter when calculating tincstp
!> Dec  2015 - half increment for 1st and incstp+1th time step
!> Dec  2015 - difference archive can be at end of insertion period
!> Jul  2017 - bugfix for steric in archive
!> Aug  2017 - incupf -ve to write a restart after each update completes
!> Sep  2017 - added subroutine stfupd
