      subroutine oiassi(j)
c
c=======================================================================
c                                                                    ===
c  Intermittent Optimal Interpolation, PE module:                    ===
c                                                                    ===
c  This scheme consists of a chaining number of assimilation cycles  ===
c  between model forecast and observation fields.  All error fields  ===
c  used in this  assimilation scheme  are assumed to be  normalized  ===
c  (by their variance) mean square errors.                           ===
c                                                                    ===
c  Calls:     NCDID, NCOPN, NCPOPT, NCVGT, NCVGT1, NCVID             ===
c             (NetCDF library)                                       ===
c             ASSIWGHT, EXITUS, LNBLK, OIBNDRY                       ===
c                                                                    ===
c  WARNING:   Character argument to NetCDF routines NCDID and NCVID  ===
c             is (upper/lower) case sensitive.                       ===
c                                                                    ===
c=======================================================================
c
c-----------------------------------------------------------------------
c  Define global data.
c-----------------------------------------------------------------------
c
#include <cdefs.h>
#include <param.h>
#include <pconst.h>
#include <fullwd.h>
#include <scalar.h>
#include <onedim.h>
#include <fields.h>
#include <workspa.h>
#include <options.h>
#include <iounits.h>
#include <rhomean.h>
#include <oiopts.h>
#include <ioi.h>
#ifdef bndy_rlx
#  include <bndyrlx.h>
#endif
#include <netcdf.inc>
#include <pe_netcdf.h>
#include <vertslabs.h>
c
c-----------------------------------------------------------------------
c  Define local data.
c-----------------------------------------------------------------------
c
      integer i,ifld,imax,j,jmax,k,kmax,lenstr,m,n,ntrc
#ifdef bndy_rlx
     *        ,ip1
#endif
      integer lnblk
      integer count(5),start(5)
      FLOAT
     *      a(imt,km),b(imt,km),fld(imtkm),tlapse
     *      ,umean(imt),vmean(imt)
#ifdef bndy_rlx
     *      ,diag1,diag2,ubwk,vbwk
#endif
      character*44 errnam,fldnam
      character*80 varnam
c
c=======================================================================
c  Begin executable code.
c=======================================================================
c
      if(ioi.eq.0) then
c
c-----------------------------------------------------------------------
c  On first call, read in assimilation parameters and switches from
c  assimilation script file unit IOIPAR.
c-----------------------------------------------------------------------
c
        ioi=1
        noi=0
        ittoi=0
c
c  Read in number of assimilation cycles (NOI), type observation error
c  (IOBSERR): [0] uniform [1] nonuniform, write out switch (IOIWRT),
c  and correlation (COR).
c
        open(ioipar,file=aparnam,status='old',form='formatted',err=20)
c
        read(ioipar,*,err=20) noi,iobserr,ioiwrt,cor,tsat
        if(noi.gt.mxoi) then
          write(stdout,900) ', MXOI = ',mxoi
          call exitus('OIASSI')
        endif
c
c  Read in assimilation times (in days) and switches for variables
c  to assimilate: [0] no, [1] yes.
c
        read(ioipar,*,err=20) (iaread(n),toisec(n),iap(n),iauv(n),
     *                         (iats(n,m),m=1,nt),n=1,noi)
c
c  Convert assimilation times to seconds.
c
        tsat=tsat*day2sec
        do 10 n=1,noi
          toisec(n)=toisec(n)*day2sec
  10    continue
c
c  Read in weight coefficients (from 0 to 1) needed to compute forecast
c  errors from provided observation errors. It is assumed that, at an
c  assimilation cycle, the model forecast can be relied by a percetage
c  derived from these weights.
c
        read(ioipar,*,err=20) ((obswgt(n,ifld),ifld=1,nt+3),n=1,noi)
c
c  If needed, read in uniform (per level) observation errors.
c
        if(iobserr.eq.0) then
          read(ioipar,*,err=20) (errp(n),n=1,noi)
          read(ioipar,*,err=20) ((erruv(n,k),k=1,km),n=1,noi)
          read(ioipar,*,err=20) (((errts(n,k,m),k=1,km),n=1,noi),m=1,nt)
        endif
        close(ioipar)
        goto 30
  20    continue
        write(stdout,919) aparnam
  30    continue
c
c  Write out summary of assimilation parameters to standard output.
c
        write(stdout,901)
        do 40 n=1,noi
          write(stdout,902) toisec(n)*sec2day,iaread(n),iap(n),iauv(n),
     *                      (iats(n,m),m=1,nt)
  40    continue
        write(stdout,903)
        do 50 n=1,noi
          write(stdout,904) toisec(n)*sec2day,obswgt(n,1),obswgt(n,2),
     *                      (obswgt(n,i),i=4,nt+3)
  50    continue
        if(iobserr.eq.0) then
          write(stdout,905)
          do 60 n=1,noi
            write(stdout,906) toisec(n)*sec2day
          do 60 k=1,km
            if(k.eq.1) then
              write(stdout,907) k,errp(n),erruv(n,k),errts(n,k,1)
            else
              write(stdout,908) k,erruv(n,k),errts(n,k,1)
            endif
  60      continue
        else
          write(stdout,909)
        endif
c
c=======================================================================
c  On first call, open input NetCDF file containing the observation
c  fields to assimilate.  Get NetCDF IDs.
c=======================================================================
c
c  Open assimilation NetCDF file.
c
        call ncpopt(ncverbos)
        lenstr=lnblk(assname,len(assname))
        ncassid=ncopn(assname(1:lenstr),ncnowrit,rcode)
        if(rcode.ne.0) then
          write(stdout,910) assname(1:lenstr)
          call exitus('OIASSI')
        endif
c
c  Read in and check main dimension parameters.
c
        varid=ncvid(ncassid,'imt',rcode)
        if(rcode.eq.0) then
          call ncvgt1(ncassid,varid,1,imax,rcode)
          if(imax.ne.imt) then
            write(stdout,911) 'IMT, IMAX: ',imt,imax
            call exitus('OIASSI')
          endif
        else
          write(stdout,912) 'imt'
          call exitus('OIASSI')
        endif
        varid=ncvid(ncassid,'jmt',rcode)
        if(rcode.eq.0) then
          call ncvgt1(ncassid,varid,1,jmax,rcode)
          if(imax.ne.imt) then
            write(stdout,911) 'JMT, JMAX: ',jmt,jmax
            call exitus('OIASSI')
          endif
        else
          write(stdout,912) 'jmt'
          call exitus('OIASSI')
        endif
        varid=ncvid(ncassid,'km',rcode)
        if(rcode.eq.0) then
          call ncvgt1(ncassid,varid,1,kmax,rcode)
          if(imax.ne.imt) then
            write(stdout,911) 'KM, KMAX: ',km,kmax
            call exitus('OIASSI')
          endif
        else
          write(stdout,912) 'km'
          call exitus('OIASSI')
        endif
        varid=ncvid(ncassid,'nt',rcode)
        if(rcode.eq.0) then
          call ncvgt1(ncassid,varid,1,ntrc,rcode)
          if(imax.ne.imt) then
            write(stdout,911) 'NT, NTRC: ',nt,ntrc
            call exitus('OIASSI')
          endif
        else
          write(stdout,912) 'nt'
          call exitus('OIASSI')
        endif
c
c  Inquire time dimension and get observations time ID and initialize
c  observations time counter.
c
        varid=ncdid(ncassid,'time',rcode)
        if(rcode.eq.0) then
          call ncdinq(ncassid,varid,varnam,tassmax,rcode)
        else
          write(stdout,913) 'time'
          call exitus('OIASSI')
        endif
        tassid=ncvid(ncassid,'time',rcode)
        if(rcode.ne.0) then
          write(stdout,912) 'time'
          call exitus('OIASSI')
        endif
        tassindx=0
c
c-----------------------------------------------------------------------
c  Get variables ID from NetCDF assimilation file.
c-----------------------------------------------------------------------
c
c  Transport streamfunction.
c
        pobsid=ncvid(ncassid,'pbar',rcode)
        if(rcode.ne.0) then
          write(stdout,912) 'pbar'
          call exitus('OIASSI')
        endif
        if(iobserr.ne.0) then
          perrid=ncvid(ncassid,'pbarerr',rcode)
          if(rcode.ne.0) then
            write(stdout,912) 'pbarerr'
            call exitus('OIASSI')
          endif
        endif
c
c  Tracers.
c
        do 70 m=1,nt
          fldnam=tname(1,m)
          lenstr=lnblk(fldnam,len(fldnam))
          tobsid(m)=ncvid(ncassid,fldnam(1:lenstr),rcode)
          if(rcode.ne.0) then
            write(stdout,912) fldnam(1:lenstr)
            call exitus('OIASSI')
          endif
          if(iobserr.ne.0) then
            errnam=fldnam(1:lenstr)//'err'
            lenstr=lnblk(errnam,len(errnam))
            terrid(m)=ncvid(ncassid,errnam(1:lenstr),rcode)
            if(rcode.ne.0) then
              write(stdout,912) errnam
              call exitus('OIASSI')
            endif
          endif
  70    continue
c
c  Internal mode velocity.
c
        vobsid=ncvid(ncassid,'vclin',rcode)
        if(rcode.ne.0) then
          write(stdout,912) 'vclin'
          call exitus('OIASSI')
        endif
        if(iobserr.ne.0) then
          verrid=ncvid(ncassid,'vclinerr',rcode)
          if(rcode.ne.0) then
            write(stdout,912) 'vclinerr'
            call exitus('OIASSI')
          endif
        endif
c
c  Error field ranges.
c
        if(iobserr.ne.0) then
          eminid=ncvid(ncassid,'oerrmin',rcode)
          if(rcode.ne.0) then
            write(stdout,912) 'oerrmin'
            call exitus('OIASSI')
          endif
          emaxid=ncvid(ncassid,'oerrmax',rcode)
          if(rcode.ne.0) then
            write(stdout,912) 'oerrmax'
            call exitus('OIASSI')
          endif
        endif
c
c-----------------------------------------------------------------------
c  On first call, Start OI assimilation cycle counter, ICOI.
c-----------------------------------------------------------------------
c
        icoi=1
      endif
c
c=======================================================================
c  End introductory section.
c=======================================================================
c
c  Activate switches on an assimilation timestep.
c
      if(ioi.eq.1) then
        if(toisec(icoi).le.ttsec) then
          ioi=2
          ittoi=1
          ioiuv=iauv(icoi)
          do 80 m=1,nt
            ioits(m)=iats(icoi,m)
  80      continue
          ioip=iap(icoi)
        endif
      endif
c
c=======================================================================
c  Begin row-by-row assimilation.
c=======================================================================
c
      if((ioi.eq.2).and.(ittoi.eq.2)) then
c
c-----------------------------------------------------------------------
c  Read in observation fields to assimilate.  If applicable, read in
c  nonuniform observation errors.
c-----------------------------------------------------------------------
c
c  Read in time since initialization and set-up time counter.
c
        if((iaread(icoi).eq.1).and.(j.eq.1)) then
          tassindx=tassindx+1
          if(tassindx.gt.tassmax) then
            write(stdout,914) 'TASSINDX, TASSMAX = ',tassindx,tassmax
            call exitus('OIASSI')
          endif
          start(1)=tassindx
          count(1)=1
          call ncvgt(ncassid,tassid,start,count,ttoi,rcode)
          if(rcode.ne.0) then
            write(stdout,915) 'time',tassindx
            call exitus('OIASSI')
          endif
        endif
c
c  Read in internal mode velocity components.
c
        start(2)=1
        count(2)=km
        start(3)=1
        count(3)=imt
        start(4)=j
        count(4)=1
        start(5)=tassindx
        count(5)=1
        if(ioiuv.ne.0) then
          start(1)=xindx
          count(1)=1
          call ncvgt(ncassid,vobsid,start,count,fld,rcode)
          if(rcode.ne.0) then
            write(stdout,916) 'vclin x-component',j
            call exitus('OIASSI')
          endif
          call flip(uobs,imt,km,fld,km,imt)
          start(1)=yindx
          count(1)=1
          call ncvgt(ncassid,vobsid,start,count,fld,rcode)
          if(rcode.ne.0) then
            write(stdout,916) 'vclin y-component',j
            call exitus('OIASSI')
          endif
          call flip(vobs,imt,km,fld,km,imt)
          if(iobserr.ne.0) then
            start(1)=xindx
            count(1)=1
            call ncvgt(ncassid,verrid,start,count,fld,rcode)
            if(rcode.ne.0) then
              write(stdout,916) 'vclinerr x-component',j
              call exitus('OIASSI')
            endif
            call flip(uobserr,imt,km,fld,km,imt)
            call ncvgt(ncassid,verrid,start,count,fld,rcode)
            if(rcode.ne.0) then
              write(stdout,916) 'vclinerr y-component',j
              call exitus('OIASSI')
            endif
            call flip(vobserr,imt,km,fld,km,imt)
          endif
        endif
c
c  Read in tracers (temperature, salinity, and other tracers if any).
c
        do 90 m=1,nt
          if(ioits(m).ne.0) then
            call ncvgt(ncassid,tobsid(m),start(2),count(2),fld,rcode)
            if(rcode.ne.0) then
              write(stdout,917) 'trc',j,m
              call exitus('OIASSI')
            endif
            call flip(tobs(1,1,m),imt,km,fld,km,imt)
            if(iobserr.ne.0) then
              call ncvgt(ncassid,terrid(m),start(2),count(2),fld,rcode)
              if(rcode.ne.0) then
                write(stdout,917) 'trcerr',j,m
                call exitus('OIASSI')
              endif
              call flip(tobserr(1,1,m),imt,km,fld,km,imt)
            endif
          endif
  90    continue
c
c  Scale salinity to gain accuracy.
c
        if(ioits(2).ne.0) then
          do 100 k=1,km
          do 100 i=1,imt
            tobs(i,k,2)=tobs(i,k,2)-smean
 100      continue
        endif
c
c  Read in transport streamfunction.
c
        if(ioip.ne.0) then
          call ncvgt(ncassid,pobsid,start(3),count(3),pobs,rcode)
          if(rcode.ne.0) then
            write(stdout,916) 'pbar',j
            call exitus('OIASSI')
          endif
          if(iobserr.ne.0) then
            call ncvgt(ncassid,perrid,start(3),count(3),pobserr,rcode)
            if(rcode.ne.0) then
              write(stdout,916) 'pbarerr',j
              call exitus('OIASSI')
            endif
          endif
        endif
c
c  Set or read in minimum and maximum observation error per field,
c  OERRMIN and OERRMAX.
c
        if((iobserr.eq.0).and.(j.eq.1)) then
          do 120 k=1,km
            oerrmin(1,k)=errp(icoi)
            oerrmax(1,k)=errp(icoi)
            oerrmin(2,k)=erruv(icoi,k)
            oerrmax(2,k)=erruv(icoi,k)
            oerrmin(3,k)=erruv(icoi,k)
            oerrmax(3,k)=erruv(icoi,k)
            do 110 m=1,nt
              oerrmin(m+3,k)=errts(icoi,k,m)
 110        continue
 120      continue
        elseif((iobserr.ne.0).and.(j.eq.1)) then
          start(1)=1
          count(1)=nt+3
          start(2)=1
          count(2)=km
          start(3)=tassindx
          count(3)=1
          call ncvgt(ncassid,eminid,start,count,oerrmin,rcode)
          if(rcode.ne.0) then
            write(stdout,918) 'oerrmin'
            call exitus('OIASSI')
          endif
          call ncvgt(ncassid,emaxid,start,count,oerrmax,rcode)
          if(rcode.ne.0) then
            write(stdout,918) 'oerrmax'
            call exitus('OIASSI')
          endif
        endif
c
c  Compute elapsed time since last asssimilation cycle.
c
        if(icoi.eq.1) then
          tlapse=toisec(icoi)
        else
          tlapse=toisec(icoi)-toisec(icoi-1)
        endif
c
c-----------------------------------------------------------------------
c  Assimilate transport stream function.
c-----------------------------------------------------------------------
c
        if(ioip.ne.0) then
          call assiwght(a,b,j,1,tlapse,1)
          do 130 i=1,imt
            p(i,j)=a(i,1)*pobs(i)+b(i,1)*p(i,j)
 130      continue
        endif
c
c-----------------------------------------------------------------------
c  Assimilate internal mode velocities.
c-----------------------------------------------------------------------
c
        if(ioiuv.ne.0) then
          if (iopt(5).eq.1) write(stdout,920) itt,j
          call assiwght(a,b,j,km,tlapse,2)
          do 140 k=1,km
          do 140 i=1,imt
            u(i,k)=a(i,k)*uobs(i,k)+b(i,k)*u(i,k)
 140      continue
          call assiwght(a,b,j,km,tlapse,3)
          do 150 k=1,km
          do 150 i=1,imt
            v(i,k)=a(i,k)*vobs(i,k)+b(i,k)*v(i,k)
 150      continue
c
c-----------------------------------------------------------------------
c  Determine the incorrect vertical means of the new velocities.
c-----------------------------------------------------------------------
c
          do 160 i=1,imt
             umean(i)=c0
             vmean(i)=c0
 160      continue
          do 170 k=1,km
          do 170 i=1,imt
             umean(i)=umean(i)+u(i,k)*dzvqz(i,k,0)
             vmean(i)=vmean(i)+v(i,k)*dzvqz(i,k,0)
 170      continue
          do 180 i=1,imt
             umean(i)=umean(i)*hvav(i,j)
             vmean(i)=vmean(i)*hvav(i,j)
 180      continue
c
c-----------------------------------------------------------------------
c  Subtract incorrect vertical mean to get internal mode.
c-----------------------------------------------------------------------
c
          do 190 k=1,km
          do 190 i=1,imt
             u(i,k)=(u(i,k)-umean(i))
             v(i,k)=(v(i,k)-vmean(i))
 190      continue

        endif
c
c-----------------------------------------------------------------------
c  Assimilate tracers.
c-----------------------------------------------------------------------
c
        do 210 m=1,nt
          if(ioits(m).ne.0) then
            if (iopt(5).eq.1) write(stdout,921)
            call assiwght(a,b,j,km,tlapse,m+3)
            do 200 k=1,km
            do 200 i=1,imt
              t(i,k,m)=a(i,k)*tobs(i,k,m)+b(i,k)*t(i,k,m)
 200        continue
          endif
 210    continue
c
c-----------------------------------------------------------------------
c  Update the boundary conditions based upon the assimilation.
c-----------------------------------------------------------------------
c
        call oibndry(j)
c
#ifdef bndy_rlx
c-----------------------------------------------------------------------
c  Update boundary relaxation arrays based upon the assimilation.
c-----------------------------------------------------------------------
c
      do 240 i=1,imt
c
         if (j.gt.2) then
c
c           -- Compute barotropic component for previous row.
c
            ip1  = min( i+1, imt)
            diag1=p (ip1,j)-p (i  ,j-1)
            diag2=p (i  ,j)-p (ip1,j-1)
            ubwk =-(diag1+diag2)*dyu2r(j-1)*hv(i,j-1)
            vbwk = (diag1-diag2)*dxu2r(i  )*hv(i,j-1)*csr(j-1)
c
         end if
c
         do 230 k=1,km
c
            u_0(i,j,k)=u(i,k)
            v_0(i,j,k)=v(i,k)
c
            if (j.gt.1) then
               u_0(i,j-1,k)=u_0(i,j-1,k)+ubwk
               v_0(i,j-1,k)=v_0(i,j-1,k)+vbwk
            end if
c
            if (j.eq.jmt) then
               u_0(i,j,k)=u_0(i,j,k)+ubwk
               v_0(i,j,k)=v_0(i,j,k)+vbwk
            end if
c
            do 220 m=1,nt
               t_0(i,j,k,m)=t(i,k,m)
 220        continue
c
 230     continue
c
 240  continue
#endif
c-----------------------------------------------------------------------
c  If last row, prepare for next assimilation cycle.
c-----------------------------------------------------------------------
c
        if(j.eq.jmt) then
          icoi=icoi+1
          if(icoi.gt.noi) then
            ioi=-1
          else
            ioi=1
          endif
        endif
      endif
#ifdef sunflush
c
c  Flush output buffers.
c
      call flush(stdout)
#endif
c
 900  format(/' OIASSI - small dimension parameter ',a,i4)
 901  format(/' OIASSI - Assimilation switches: [0] no, [1] yes'//
     *       15x,'DAY',5x,'READ',6x,'PSI',6x,'U,V',6x,'TRACERS'/)
 902  format(10x,f8.4,3(7x,i2),7x,6i2)
 903  format(/' OIASSI - Observation weights used to compute',
     *        ' forecast errors'//12x,'DAY',7x,'PSI',6x,'U,V',6x,
     *        'TRACERS'/)
 904  format(9x,f8.4,4(1x,f8.4))
 905  format(/' OIASSI - Uniform observation errors')
 906  format(/4x,' Time = ',f8.4//12x,'LEVEL',8x,'PSI',12x,'U,V',12x,
     *       'TRACERS'/)
 907  format(15x,i2,3(3x,1pe12.5))
 908  format(15x,i2,15x,2(3x,1pe12.5))
 909  format(/' OIASSI:  Expecting nonuniform observation errors')
 910  format(/' OIASSI - unable to open assimilation NetCDF file: ',a)
 911  format(/' OIASSI - inconsistent dimension parameters, ',a,2i4)
 912  format(/' OIASSI - cannot find variable: ',a,2x,
     *        ' in assimilation NetCDF file.')
 913  format(/' OIASSI - cannot find dimension: ',a,2x,
     *        ' in assimilation NetCDF file.')
 914  format(/' OIASSI - time dimension index exceeded, ',a,2x,2i4)
 915  format(/' OIASSI - error while reading variable: ',a,2x,
     *        ' at index = ',i4)
 916  format(/' OIASSI - error while reading variable: ',a,2x,
     *        ' at row J = ',i5)
 917  format(/' OIASSI - error while reading variable: ',a,2x,
     *        ' at row J = ',i5,2x,' and TRACER = ',i2)
 918  format(/' OIASSI - error while reading variable: ',a)
 919  format(/' OIASSI - error while reading assimilation parameters',
     *        ' from input file = ',a)
 920  format('Assimilating internal mode velocities (itt,j):',2(1x,i7))
 921  format('Assimilating tracers')
c
      return
      end
