C 
C     AUTHOR:             ALFIO BORZI'
C                         INSTITUT FÜR MATHEMATIK
C                         UND WISSENSCHAFTLICHES RECHNEN
C                         KARL-FRANZENS-UNIVRSITÄT GRAZ
C                         HEINRICHSTR. 36, A-8010 GRAZ
C                         AUSTRIA
C                         E-MAIL: ALFIO.BORZI@UNI-GRAZ.AT
C
C     DATE:         2004
C     PROJECT:      FAS MULTIGRID FOR REACTION-DIFFUSION OPTIMALITY SYSTEMS 
C                   OPEN-LOOP CONTROL AND RECEDING HORIZON TECHNIQUES
C     LANGUAGE:     FORTRAN 77
C     CONSTRAINTS:  NONE - PUBLIC DOMAIN SOFTWARE   
C
C     FOR THE PRESENTATION OF THIS CODE AND FURTHER REFERENCES SEE:
C     ALFIO BORZI'
C     Space-time multigrid methods for solving unsteady 
C     optimal control problems
C
C     appear in the proceedings of the
C     
C     Second Sandia Workshop on PDE-Constrained Optimization:
C     Toward Real-time and Online PDE-constrained Optimization
C	May 19-21, 2004
C	Bishop's Lodge, Santa Fe, New Mexico
C    
C     Workshop organizing committee:
C
C     Larry Biegler, Omar Ghattas, Matthias Heinkenschloss, 
C     David Keyes, Bart van Bloemen Waanders.
C        
C
       program santafe

       implicit real*8 (a-h,o-z)
c
c description: min J(u,f),  e(u,f)=0.
c
c constraint: 
c Parabolic (2+1)-dimensional FAS Multigrid solver (with FDM)
c Dirichlet b.c. + initial condition, 
c
c  Optimal control problem:
c
C  -Ut + GNL(U) + sigma * DELTA(U) = F    
C
C  min J(U,F):= alfa ||U-Zd||^2/2 + beta ||U(T)-ZT||^2/2 + rnu ||F||^2/2 
C  (distributed control)
C  
C  with Receding horizon techniques .....
C  
c 
c
c running indeces for space and time: i->x, j->y, n->t. 
c
c                 q(nst(i,j,le)+n)=var(x,y,t)
c                 q(nst(i,j,le+ne)+n)=rhs(var)
c
c   le=1,ne - le index of the equation; ne number of equations
c
c   The storage required is
c                               d   d
c                       2*ne*n*2 /(2 - 1)
c
c   where
c                ne   is the number of equations
c                n       the number of finest grid points
c                d       dimensions of the domain
c
c
c   Residuals are transferred by full weighting
c
c      meaning of some flags
c
c      iwrite      1     dump the solution (for debugging purpose)
c                  2     don't dump the solution
c      nvw         1     v cycle
c                  2     w cycle
c
c      x0 - x1     domain of integration in x
c      y0 - y1     domain of integration in y
c      t0 - t1     domain of integration in t
c      nx0         # of meshes of the coarsest grid in
c                  x direction
c      ny0         # of meshes of the coarsest grid in
c                  y direction
c      nt0         # of meshes of the coarsest grid in
c                  t (time) direction
c      m           # of level used
c      ncycle      # of cycle for each MG cycle
c      nrf         # of sweeps before CGC
c      nrc         # of sweeps after CGC
c      isor = 1    performs max relax sweeps starting from level
c                  'ksor'
c      itstl= 1    TS-CGS relax; itstl= 2 TL-CGS relax 
c      nwin        number of time windows of size t1-t0 for receding horizon
c
      common/qarr/q(50000000)
      common/pntrs/ist(129,129),isf(129,129,8),
     X             ism(129,129),isc(129,129,8)
      common/sor/isor,nmax,ksor,itstl
      common/grid/x0,x1,y0,y1,t0,t1
      common/twindow/t00
      common/tol/tol
      common/albe/sigma
      common/cweight/rnu,alfa,beta
      common/horiz/nwin
c
      external fsol,fqin,frhs,fqb
      save
c
      qdim=50000000
c
      open(unit=14,status='unknown',file='out.dat')
      open(unit=17,status='unknown',file='midsol.dat',access='append')
      open(unit=20,status='unknown',file='santafe.dat')
c
      pi=acos(-1.0d+0)
      t00=0.0d0
c
      ne=2
      read(20,*) x0,x1,y0,y1,t0,t1
      read(20,*) nx0,ny0,nt0,m
      read(20,*) ncycle,nrc,nrf,ncrst,nvw
      read(20,*) tol,wmax
      read(20,*) isor,nmax,ksor,itstl
      read(20,*) sigma
      read(20,*) rnu,alfa,beta
      read(20,*) nwin
      close(20)
      t00=t0      
c
      write(14,80)
      write(14,99)
      write(14,100) nx0,ny0,nt0
      write(14,105) m,ncycle
      write(14,106) nrf
      write(14,107) nrc
      write(14,110) x0,x1,y0,y1,t0,t1
c
      call cpu_time(time1)
c      
      call multig(nx0,ny0,nt0,x0,y0,t0,x1,y1,t1,m,qdim,wmax,fsol,
     *            frhs,fqin,fqb,ne,ncycle,nrc,nrf,ncrst,nvw)
      call cpu_time(time2)
      write(*,*)' Total CPU time (s) ',time2-time1
      write(14,175) time2-time1
c
   80 format(19x,44('*'),/,19x,'*',42x,'*',/,19x,'*  ','Multigrid sol',
     *      'ution of Poisson equation','  *',/,19x,'*',42x,'*',/,
     *      19x,44('*'),///)
   99 format('   coarsest grid:')
  100 format('   # of x-intervals',i5,
     *           '  # of y-intervals',i5,'  # of t-intervals',i5)
  105 format('   # of levels used =',i3,' # of mgm-cycles per level ='
     *       ,i3)
  106 format('   # of sweeps before coarse grid correction =',i3)
  107 format('   # of sweeps after coarse grid correction =',i3)
  110 format('   domain of integration :   x0=',f9.3,'  x1=',f9.3,
     *       /,29x,'y0=',f9.3,'  y1=',f9.3,
     *       /,29x,'t0=',f9.3,'  t1=',f9.3)
  175 format('   total cpu time = ',f12.2,' secs')
      stop
      end
c------------------------------------------------
      real*8 function gnl(u)
      implicit real*8 (a-h,o-z)
c
c---- the nonlinear function
c
      gnl=exp(u)
      return
      end
c------------------------------------------------
      real*8 function d1gnl(u)
      implicit real*8 (a-h,o-z)
c
c---- the 1st derivative of the nonlinear function
c
      d1gnl=exp(u)
      return
      end
c------------------------------------------------
      real*8 function d2gnl(u)
      implicit real*8 (a-h,o-z)
c
c---- the 2nd derivative of the nonlinear function
c
      d2gnl=exp(u)
      return
      end
c------------------------------------------------
      real*8 function zf(x,y,t)
      implicit real*8 (a-h,o-z)
      common/cweight/rnu,alfa,beta
      pi=acos(-1.0d+0)
      pi2=pi*pi
      pi4=pi2*pi2
c
c---- the objective function
c
      zf=(1.+t)*(x-x**2)*(y-y**2)*cos(4.0*pi*t)
      return
      end
c------------------------------------------------
      real*8 function fsol(x,y,t,i)
      implicit real*8 (a-h,o-z)
      common/cweight/rnu,alfa,beta
c
c---- the analytic solution (?)
c
      pi=acos(-1.0d+0)
      pi2=pi*pi
      if(i.eq.1) then
      fsol=(1.+t)*(x-x**2)*(y-y**2)*cos(4.0*pi*t)
      else
      fsol=0.0
      endif
      return
      end
c------------------------------------------------
      function fqb(x,y,t,i)
      implicit real*8 (a-h,o-z)
      if(i.eq.1) then
         fqb=0.0
      else
         fqb=0.0
      endif
      return
      end
c------------------------------------------------
      real*8 function fqin(x,y,t,i)
      implicit real*8 (a-h,o-z)
      external fsol
      if(i.eq.1) then
         fqin=fsol(x,y,t,1)
      else
         fqin=fsol(x,y,t,2)
      endif
      return
      end
c------------------------------------------------
      function frhs(x,y,t,i)
      implicit real*8 (a-h,o-z)
      pi=acos(-1.0d+0)

      if(i.eq.1) then
        frhs=0.0
      else
        frhs=0.0
      endif
      return
      end

c------------------------------------------------
      subroutine crsres(k,ne)
      implicit real*8 (a-h,o-z)
      common/qarr/q(50000000)
      common/pntrs/ist(129,129),is(129,129,8),
     X             isp(129,129),isf(129,129,8)
      common/albe/sigma
      common/cweight/rnu,alfa,beta
      external zf,gnl,d1gnl
c
      n2=2*ne
      do 2 l=1,n2
      call key(k,l,ist,imax,jmax,nmax,x0,y0,t0,hx,hy,ht)
      do 1 i=1,imax
      do 1 j=1,jmax
    1 is(i,j,l)=ist(i,j)
    2 continue
c
      imax1=imax-1
      jmax1=jmax-1
      nmax1=nmax-1
c
      coe1=sigma*ht/(hx*hx)
      coe2=sigma*ht/(hy*hy)
      coe3=1.0
      coe=2.0d0*(coe1+coe2)+coe3
c
c* state equation
      le=1
      do 20 i=2,imax1
      do 20 j=2,jmax1
      iom=is(i-1,j,le)
      io=is(i,j,le)
      iop=is(i+1,j,le)
      jom=is(i,j-1,le)
      jop=is(i,j+1,le)
      ifo=is(i,j,ne+le)
c
      x=(i-1)*hx
      y=(j-1)*hy

      do 40 n=2,nmax
      nm=n-1
      t=t0+(n-1)*ht

      ax=q(iom+n)-2.0d0*q(io+n)+q(iop+n)
      ay=q(jom+n)-2.0d0*q(io+n)+q(jop+n)
      at=q(io+nm)-q(io+n)
      a1=q(ifo+n)+(coe1*ax+coe2*ay+coe3*at)
      a1=a1+ht*gnl(q(io+n))-ht*q(is(i,j,2)+n)/rnu
      q(ifo+n)=a1
   40 continue
   20 continue
c
c* adjoint equation
      le=2
      do 30 i=2,imax1
      do 30 j=2,jmax1
      iom=is(i-1,j,le)
      io=is(i,j,le)
      iop=is(i+1,j,le)
      jom=is(i,j-1,le)
      jop=is(i,j+1,le)
      ifo=is(i,j,ne+le)

      x=(i-1)*hx
      y=(j-1)*hy
      do 30 na=2,nmax1
      n=nmax-na+1
      nm=n+1
      t=t0+(n-1)*ht

      ax=q(iom+n)-2.0d0*q(io+n)+q(iop+n)
      ay=q(jom+n)-2.0d0*q(io+n)+q(jop+n)
      at=q(io+nm)-q(io+n)
      a1=q(ifo+n)+(coe1*ax+coe2*ay+coe3*at)
      yvar=q(is(i,j,1)+n)
      a1=a1+ht*d1gnl(yvar)*q(io+n)+ht*alfa*(yvar-zf(x,y,t))
      q(ifo+n)=a1
   30 continue
      return
      end

c------------------------------------------------
      subroutine fas(fs,lev,mm,wu,ne,
     *               ncycle,nrc,nrf,ncrst,nvw)
      implicit real*8 (a-h,o-z)
c
c     fix algorithm
c
      dimension error(260),work(260)
      dimension err(12),ivw(10)
      common/qarr/q(50000000)
      common/pntrs/ist(129,129),isf(129,129,8),
     X             ism(129,129),isc(129,129,8)
      common/tol/tol
      common/cweight/rnu,alfa,beta
c
      external fs
c
      do 1 n=1,12
    1 err(n)=0.0d0
      write(14,170)
c
      m=mm
      lwork=lev
      do 5 k=1,m
    5 ivw(k)=0
      errmx=1.0d0
c
c* cycle
      do 100 nc=1,ncycle
      k=lev
      errold=errmx
c
   10 do 20 ir=1,nrf
      call relax(k,err,errm,ne)
c
      wu=wu+4.**(k-lwork)
      write(14,190) k,wu,(err(i),i=1,ne)
   20 continue
      if(lev.eq.1) goto 80
c
c  coarse grid correction
c  ----------------------
c                          H   H   h   h h     H  H h
c                         r = I  (f - L u ) + L (I u )
c                              h                  h
c
      call rescal(k,k-1,ne)
      call putu(k,k-1,ne)
      call crsres(k-1,ne)
c
      ivw(k)=ivw(k)+1
      k=k-1
      if(k.ne.1) goto 10
c
c     make now ncrst sweeps on the coarsest grid
c
      do 42 n=1,ncrst
      call relax(k,err,errm,ne)
      write(14,190) k,wu,(err(i),i=1,ne)
   42 continue
      if(ivw(lev).ne.1) goto 44

   44 write(14,190) k,wu,(err(i),i=1,ne)
c
   45 k=k+1
c
c   transfer to finer grid and addition
c   -----------------------------------
c
c                            h   h   h  H   h H
c                           u = u + I (u - I u )
c                                    H      H
      call subtrt(k-1,k,ne)
      call intadd(k-1,k,ne)
c
      do 60 ir=1,nrc
      call relax(k,err,errm,ne)
      wu=wu+4.0d0**(k-lwork)
      write(14,190) k,wu,(err(i),i=1,ne)
   60 continue
c
      if(k.ne.2) goto 65
      if(ivw(k).ne.1) goto 10
      ivw(k)=0
      if(k.eq.lev) goto 80
      goto 45
   65 if(k.eq.lev) goto 80
      if(ivw(k).ne.nvw) goto 10
      ivw(k)=0
      goto 45
c
   80 work(nc)=wu
c
      error(nc)=err(1)+err(2)
      write(14,200) nc
      if(k.eq.lev) errmx=(err(1)+err(2))
      write(14,130) nc, errmx/errold
      write(*,130) nc, errmx/errold
      if(errmx.lt.tol) goto 150
  100 continue
c
  150 if(ncycle.le.2) return
      ncf=nc
      do 110 nc=2,ncf
      tw=work(nc)-work(nc-1)
      ta=work(nc)-work(1)
      ew=error(nc)/error(nc-1)
      if(ew.lt.1e-20)ew=1e-20
      ea=error(nc)/error(1)
      if(ea.lt.1e-20)ea=1e-20
      efw=ew**(1./tw)
      efa=ea**(1./ta)
      write(14,210) nc,ew,ea
      write(*,210) nc,ew,ea      
  110 continue
c
  130 format(10x,'nc =',i4,15x,'obs. rho =',f8.4)
  170 format(' level',2x,'work',18x,'residual norms (i,b)')
  190 format(i2,2x,f6.2,2x,2(1pe8.1))
  200 format(1x,24(1h-),'end cycle no.',i3)
  210 format(14x,'cycle no.',i3,5x,' cycle eff.',f6.3,5x,'acc.eff.'
     1,f6.3)
      return
      end
c------------------------------------------------
      subroutine grdfn(k,l,imax,jmax,nmax,x0a,y0a,t0a,hhx,hhy,hht)
      implicit real*8 (a-h,o-z)
      common/grd/nst(20,10),imx(20),jmx(20),nmx(20),
     *           hx(20),hy(20),ht(20),x0,y0,t0
      common/iq/iq
      nst(k,l)=iq
      imx(k)=imax
      jmx(k)=jmax
      nmx(k)=nmax
      x0=x0a
      y0=y0a
      t0=t0a
      hx(k)=hhx
      hy(k)=hhy
      ht(k)=hht
      iq=iq+imax*jmax*nmax
      return
      end

c-----------------------------------------------------------------
      subroutine intadd(km,k,ne)
      implicit real*8 (a-h,o-z)
c
      common/qarr/q(50000000)
      common/pntrs/ist(129,129),isf(129,129,8),
     X             ism(129,129),isc(129,129,8)

      do 30 le=1,ne
      call key(k,le,ist,iif,jjf,nnf,x0,y0,t0,hxf,hyf,htf)
      call key(km,le,ism,iic,jjc,nnc,x0,y0,t0,hxc,hyc,htc)
c
c no time coarsening
c
c
c side i=1
c
      ic=1
      if=2*ic-1
      do 55 jc=2,jjc
      jf=2*jc-1
      ifq=ist(if,jf)
      jfqm=ist(if,jf-1)
      icq=ism(ic,jc)
      jcqm=ism(ic,jc-1)
      do 55 nc=2,nnc
      nf=nc
c
      a=.50d0*(q(icq+nc)+q(jcqm+nc))
      q(ifq+nf)=q(ifq+nf)+q(icq+nc)
      q(jfqm+nf)=q(jfqm+nf)+a
   55 continue
c
c edge i=1,j=1
      ic=1
      if=1
      jc=1 
      jf=1
c corner i=1,j=1,n=1
      ifq=ist(if,jf)
      icq=ism(ic,jc)
      q(ifq+1)=q(ifq+1)+q(icq+1)
c
      do 56 nc=2,nnc
      nf=nc
c
      q(ifq+nf)=q(ifq+nf)+q(icq+nc)
   56 continue
c side j=1
      jc=1
      jf=2*jc-1
      do 510 ic=2,iic
      if=2*ic-1
      ifq=ist(if,jf)
      ifqm=ist(if-1,jf)
      icq=ism(ic,jc)
      icqm=ism(ic-1,jc)
      do 510 nc=2,nnc
      nf=nc
c
      am=.50d0*(q(icqm+nc)+q(icq+nc))
      q(ifq+nf)=q(ifq+nf)+q(icq+nc)
      q(ifqm+nf)=q(ifqm+nf)+am
  510 continue
c
c edge j=1,n=1
      jc=1
      jf=1
      nc=1
      nf=1
      do 511 ic=2,iic
      if=2*ic-1
      ifq=ist(if,jf)
      ifqm=ist(if-1,jf)
      icq=ism(ic,jc)
      icqm=ism(ic-1,jc)
c
      am=.50d0*(q(icqm+nc)+q(icq+nc))
      q(ifq+nf)=q(ifq+nf)+q(icq+nc)
      q(ifqm+nf)=q(ifqm+nf)+am
  511 continue
c side n=1
      nc=1
      nf=1
      do 515 ic=2,iic
      if=2*ic-1
      do 515 jc=2,jjc
      jf=2*jc-1
      ifq=ist(if,jf)
      ifqm=ist(if-1,jf)
      jfqm=ist(if,jf-1)
      jfqmm=ist(if-1,jf-2)
      ijfqm=ist(if-1,jf-1)
      icq=ism(ic,jc)
      icqm=ism(ic-1,jc)
      jcqm=ism(ic,jc-1)
      ijcqm=ism(ic-1,jc-1)
c
      a=.50d0*(q(icq+nc)+q(jcqm+nc))
      q(ifq+nf)=q(ifq+nf)+q(icq+nc)
      q(jfqm+nf)=q(jfqm+nf)+a

      a=.50d0*(q(icqm+nc)+q(icq+nc))
      am=.250d0*(q(icq+nc)+q(icqm+nc)
     X          +q(jcqm+nc)+q(ijcqm+nc))
      q(ifqm+nf)=q(ifqm+nf)+a
      q(ijfqm+nf)=q(ijfqm+nf)+am
  515 continue
c edge n=1,i=1
      nc=1
      nf=1
      ic=1 
      if=1
      do 516 jc=2,jjc
      jf=2*jc-1
      ifq=ist(if,jf)
      jfqm=ist(if,jf-1)
      icq=ism(ic,jc)
      jcqm=ism(ic,jc-1)
      a=.50d0*(q(icq+nc)+q(jcqm+nc))
      q(ifq+nf)=q(ifq+nf)+q(icq+nc)
      q(jfqm+nf)=q(jfqm+nf)+a
  516 continue
c internal points and outer sides
      do 520 ic=2,iic
      if=2*ic-1
      do 520 jc=2,jjc
      jf=2*jc-1
      ifq=ist(if,jf)
      ifqm=ist(if-1,jf)
      jfqm=ist(if,jf-1)
      jfqmm=ist(if-1,jf-2)
      ijfqm=ist(if-1,jf-1)
      icq=ism(ic,jc)
      icqm=ism(ic-1,jc)
      jcqm=ism(ic,jc-1)
      ijcqm=ism(ic-1,jc-1)
      do 520 nc=2,nnc
      nf=nc
c
      a=.50d0*(q(icq+nc)+q(jcqm+nc))
      q(ifq+nf)=q(ifq+nf)+q(icq+nc)
      q(jfqm+nf)=q(jfqm+nf)+a

      a=.50d0*(q(icqm+nc)+q(icq+nc))
      am=.250d0*(q(icqm+nc)+q(icq+nc)+q(jcqm+nc)+q(ijcqm+nc))
      q(ifqm+nf)=q(ifqm+nf)+a
      q(ijfqm+nf)=q(ijfqm+nf)+am
  520 continue

   30 continue
      return
      end
c------------------------------------------------
      subroutine key(k,l,ist,imax,jmax,nmax,x0a,y0a,t0a,hhx,hhy,hht)
      implicit real*8 (a-h,o-z)
      common/grd/nst(20,10),imx(20),jmx(20),nmx(20),
     *           hx(20),hy(20),ht(20),x0,y0,t0
      common/twindow/t00
      dimension ist(129,129)
      imax=imx(k)
      jmax=jmx(k)
      nmax=nmx(k)
      is=nst(k,l)-nmax-1
      do 1 i=1,imax
         do 1 j=1,jmax
            is=is+nmax
            ist(i,j)=is
    1 continue
      x0a=x0
      y0a=y0
      t0a=t00
      hhx=hx(k)
      hhy=hy(k)
      hht=ht(k)
      return
      end
c------------------------------------------------
      subroutine multig(nx0,ny0,nt0,x0,y0,t0,x1,y1,t1,mm,
     * qdim,wmax,fs,fr,fi,fb,ne,ncycle,nrc,nrf,ncrst,nvw)
      implicit real*8 (a-h,o-z)
      dimension err(12)
      common/qarr/q(50000000)
      common/pntrs/ist(129,129),isf(129,129,8),
     X             ism(129,129),isc(129,129,8)
      common/sor/isor,nmax,ksor,itstl
      common/albe/sigma
      common/cweight/rnu,alfa,beta
      common/iqco/iqco
      common/finest/m
      common/tol/tol
      common/iq/iq
      common/twindow/t00
      common/horiz/nwin
      external fs,fr,fi,fb
c
      iq=1
      iqco=1
c
      wu=0.0d0
      m=mm
c
      hx0=(x1-x0)/nx0
      hy0=(y1-y0)/ny0
      ht0=(t1-t0)/nt0
      np=(nx0+1)*(ny0+1)*(nt0+1)
      n2=2*ne
      do 10 k=1,m
      k2=2**(k-1)
      k4=1
      npx=nx0*k2+1
      npy=ny0*k2+1
      npt=nt0*k4+1
      hx=hx0/k2
      hy=hy0/k2
      ht=ht0/k4
      write(*,*) 'k nx nt',k,npx,npt
c
      do 5 l=1,n2
      call grdfn(k,l,npx,npy,npt,x0,y0,t0,hx,hy,ht)
    5 continue
   10 continue
c
      if(iq.gt.qdim) write(*,*)' iq qdim ',iq,qdim
      if(iq.gt.qdim) stop'mem'
c
c     nwin=number of time windows with size dtwin
c
      dtwin=t1-t0
      ery=0.0
      erp=0.0
      erz=0.0
      do 100 iw=1,nwin
      write(*,*) '      time window: ', iw,dtwin
      t0=(iw-1)*dtwin
      t00=t0
      t1=t0+dtwin
      write(*,*) '      from ',t0,' to ',t1
      do 1 n=1,12
    1 err(n)=0.0d0
c
      do 12 k=1,m
      call putf(k,fr,ne)
   12 continue
c
      do 19 k=1,m
   19 call putb(k,ne,iw)
c
      if(isor.eq.0) goto 25
c
c  solve by sor
c
      wusor=0.
      errmx=1.0d0
      write(14,170)
      do 24 l=ksor,m
      write(14,*) 'Level ',l
      do 23 nrel=1,nmax
      errold=errmx
      call relax(l,err,errm,ne)
      errmx=(err(1)+err(2))
      wusor=wusor+1.
      write(14,130) errmx/errold
      write(*,130) errmx/errold
      res=(err(1)+err(2))
      write(14,190) wusor, (err(i),i=1,ne)
   23 continue
   24 continue
      return
c
c   multigrid
c
   25 continue
      k=m
      call fas(fs,k,m,wu,ne,ncycle,nrc,nrf,ncrst,nvw)
      call store(m,ne,iw,dtwin,erry,errp,errz)
      ery=ery+erry*erry
      erp=erp+errp*errp
      erz=erz+errz*errz
  100 continue
      write(*,*) '   ' 
      write(*,*) ' Receding horizon: Tracking   ' 
      write(*,*) ' Total time ', dtwin*nwin, ' No. windows ',nwin
      write(*,*) ' Total tracking error ', sqrt(erz)
c
  130 format(45x,'obs. rho =',f8.4)
  170 format('level',4x,'work',6x,'residual norms(i,b)')
  190 format('WU   ',f8.4,5x,2(1pe10.3))
      return
      end

c------------------------------------------------
      subroutine putb(k,ne,iwin)
      implicit real*8 (a-h,o-z)
      common/qarr/q(50000000)
      common/pntrs/ist(129,129),isf(129,129,8),
     X             ism(129,129),isc(129,129,8)
      external fsol
c
      do 30 l=1,ne
      call key(k,l,ist,ii,jj,nn,x0,y0,t0,hx,hy,ht)

c* lateral boundary
      i=1
      x=x0
      do 2 j=1,jj
      io=ist(i,j)
      y=y0+(j-1)*hy
      do 2 n=1,nn
      t=t0+(n-1)*ht
      q(io+n)=0.0
    2 continue
c
      i=ii
      x=x0+(i-1)*hx
      do 4 j=1,jj
      io=ist(i,j)
      y=y0+(j-1)*hy
      do 4 n=1,nn
      t=t0+(n-1)*ht
      q(io+n)=0.0
    4 continue
c
      j=1
      y=y0
      do 6 i=1,ii
      x=x0+(i-1)*hx
      io=ist(i,j)
      do 6 n=1,nn
      t=t0+(n-1)*ht
      q(io+n)=0.0
    6 continue
c
      j=jj
      y=y0+(j-1)*hy
      do 8 i=1,ii
      x=x0+(i-1)*hx
      io=ist(i,j)
      do 8 n=1,nn
      t=t0+(n-1)*ht
      q(io+n)=0.0
    8 continue
c* i.c.
      ii1=ii-1
      jj1=jj-1
c* initial condition for the state eq.
      if(l.eq.1) then
      if(iwin.eq.1) then
      n=1
      t=t0
      do 10 i=2,ii1
      x=x0+(i-1)*hx
      do 10 j=2,jj1
      y=y0+(j-1)*hy
      io=ist(i,j)
      q(io+n)=fsol(x,y,t,1)
   10 continue
      else 
      n=1
      t=t0
      do 11 i=2,ii1
      x=x0+(i-1)*hx
      do 11 j=2,jj1
      y=y0+(j-1)*hy
      io=ist(i,j)
      q(io+n)=q(io+nn)
   11 continue
      endif
      endif
c
c* final condition for the adjoint eq.
      if(l.eq.2) then
      n=nn
      t=t0+(nn-1)*ht
      do 12 i=2,ii1
      x=x0+(i-1)*hx
      do 12 j=2,jj1
      y=y0+(j-1)*hy
      io=ist(i,j)
cc      q(io+n)=0.0
   12 continue
      if(iwin.gt.1) then
      n=1
      do 14 i=2,ii1
      do 14 j=2,jj1
      io=ist(i,j)
      q(io+n)=q(io+nn)
   14 continue
      endif
      endif
   30 continue
c
c initialize
c      
      do 13 l=1,ne
      call key(k,l,ist,imax,jmax,nmax,x0,y0,t0,hx,hy,ht)
      do 1 i=1,imax
      x=x0+(i-1)*hx
      do 1 j=1,jmax
      io=ist(i,j)
      y=y0+(j-1)*hy
      do 1 n=2,nmax
      t=t0+(n-1)*ht
      q(io+n)=0.0
    1 continue
   13 continue
      return
      end


c------------------------------------------------
      subroutine putf(m,fr,ne)
      implicit real*8 (a-h,o-z)
      common/qarr/q(50000000)
      common/pntrs/ist(129,129),isf(129,129,8),
     X             ism(129,129),isc(129,129,8)
      common/albe/sigma
      external fr
      k=m
c
      do 20 l=1,ne
      call key(k,l+ne,ist,imax,jmax,nmax,x0,y0,t0,hx,hy,ht)
      hs=ht
c* state rhs
      if(l.eq.1) then
      do 3 i=1,imax
      x=x0+(i-1)*hx
      do 3 j=1,jmax
      io=ist(i,j)
      y=y0+(j-1)*hy
      do 3 n=1,nmax
      t=t0+(n-1)*ht
      q(io+n)=fr(x,y,t,l)*hs
    3 continue
      else
c* adjoint rhs
      do 6 i=1,imax
      x=x0+(i-1)*hx
      do 6 j=1,jmax
      io=ist(i,j)
      y=y0+(j-1)*hy
      do 6 na=1,nmax
      n=nmax-na+1
      t=t0+(n-1)*ht
      q(io+n)=fr(x,y,t,l)*hs
    6 continue
      endif
   20 continue
      return
      end
c------------------------------------------------
      subroutine putu(k,km,ne)
      implicit real*8 (a-h,o-z)
      common/qarr/q(50000000)
      common/pntrs/ist(129,129),isf(129,129,8),
     X             ism(129,129),isc(129,129,8)
      do 20 l=1,ne
      call key(k,l,ist,iif,jjf,nnf,x0,y0,t0,hxf,hyf,htf)
      call key(km,l,ism,iic,jjc,nnc,x0,y0,t0,hxc,hyc,htc)
      do 2 i=1,iif
      do 2 j=1,jjf
    2 isf(i,j,l)=ist(i,j)
      do 3 i=1,iic
      do 3 j=1,jjc
    3 isc(i,j,l)=ism(i,j)
      iic1=iic-1
      jjc1=jjc-1
      nnc1=nnc-1
c
      do 9 ic=1,iic
      if=2*ic-1
      do 9 jc=1,jjc
      jf=2*jc-1
      icq=isc(ic,jc,l)
      ifq=isf(if,jf,l)
      do 9 nc=1,nnc
      nf=nc
      q(icq+nc)=q(ifq+nf)
    9 continue
   20 continue
      return
      end

c------------------------------------------------
      subroutine relax(k,err,errmax,ne)
      implicit real*8 (a-h,o-z)
      dimension err(12)
      common/pntrs/ist(129,129),isf(129,129,8),
     X             ism(129,129),isc(129,129,8)
      common/sor/isor,nmax,ksor,itstl
      common/tol/tol
c
c*    TS/TL - Gauss-Seidel
      if(itstl.eq.1) call relaxts(k,err,errmax,ne)
      if(itstl.eq.2) call relaxtl(k,err,errmax,ne)
c
      write(*,100) k, (err(i),i=1,ne)
  100 format('err(i,e)(',i2,')=',d13.6,4x,d13.6)
      return
      end



c------------------------------------------------
      subroutine resca1(k,ne,le)
      implicit real*8 (a-h,o-z)
      common/qarr/q(50000000)
      common/pntrs/ist(129,129),isf(129,129,8),
     X             ism(129,129),isc(129,129,8)
      common/res/resf(129,129,830)
      common/albe/sigma
      common/cweight/rnu,alfa,beta
      external gnl
c
      do 3 l=1,ne
      call key(k,l,ist,iif,jjf,nnf,x0,y0,t0,hxf,hyf,htf)
      do 1 i=1,iif
      do 1 j=1,jjf
    1 isf(i,j,l)=ist(i,j)

      call key(k,l+ne,ist,iif,jjf,nnf,x0,y0,t0,hxf,hyf,htf)
      do 3 i=1,iif
      do 3 j=1,jjf
    3 isf(i,j,l+ne)=ist(i,j)
c
c  compute residuals on the fine grid
c
      coe1=sigma*htf/(hxf*hxf)
      coe2=sigma*htf/(hyf*hyf)
      coe3=1.0
      coe=2.0d0*(coe1+coe2)+coe3

      iif1=iif-1
      jjf1=jjf-1
      nnf1=nnf-1
c
      do 10 if=2,iif1
      do 10 jf=2,jjf1
      ifom=isf(if-1,jf,le)
      ifo=isf(if,jf,le)
      ifop=isf(if+1,jf,le)
      jfom=isf(if,jf-1,le)
      jfop=isf(if,jf+1,le)
      iffo=isf(if,jf,ne+le)

      ifo2=isf(if,jf,2)
      x=(if-1)*hxf
      y=(jf-1)*hyf
c
      do 10 nf=2,nnf
      nfm=nf-1
      t=t0+(nf-1)*htf
      ax=q(ifom+nf)-2.0d0*q(ifo+nf)+q(ifop+nf)
      ay=q(jfom+nf)-2.0d0*q(ifo+nf)+q(jfop+nf)
      at=q(ifo+nfm)-q(ifo+nf)
      a1=q(iffo+nf)-(coe1*ax+coe2*ay+coe3*at)
      a1=a1-htf*gnl(q(ifo+nf))+htf*q(ifo2+nf)/rnu
      resf(if,jf,nf)=a1
   10 continue
      return
      end
c------------------------------------------------
      subroutine resca2(k,ne,le)
      implicit real*8 (a-h,o-z)
      common/qarr/q(50000000)
      common/pntrs/ist(129,129),isf(129,129,8),
     X             ism(129,129),isc(129,129,8)
      common/res/resf(129,129,830)
      common/albe/sigma
      common/cweight/rnu,alfa,beta
      external zf,d1gnl
c
      do 3 l=1,ne
      call key(k,l,ist,iif,jjf,nnf,x0,y0,t0,hxf,hyf,htf)
      do 1 i=1,iif
      do 1 j=1,jjf
    1 isf(i,j,l)=ist(i,j)

      call key(k,l+ne,ist,iif,jjf,nnf,x0,y0,t0,hxf,hyf,htf)
      do 3 i=1,iif
      do 3 j=1,jjf
    3 isf(i,j,l+ne)=ist(i,j)
c
c  compute residuals on the fine grid
c
      coe1=sigma*htf/(hxf*hxf)
      coe2=sigma*htf/(hyf*hyf)
      coe3=1.0
      coe=2.0d0*(coe1+coe2)+coe3

      iif1=iif-1
      jjf1=jjf-1
      nnf1=nnf-1
c
      do 10 if=2,iif1
      do 10 jf=2,jjf1
      ifom=isf(if-1,jf,le)
      ifo=isf(if,jf,le)
      ifop=isf(if+1,jf,le)
      jfom=isf(if,jf-1,le)
      jfop=isf(if,jf+1,le)
      iffo=isf(if,jf,ne+le)

      ifo1=isf(if,jf,1)
      x=(if-1)*hxf
      y=(jf-1)*hyf
c
      do 10 nfa=2,nnf1
      nf=nnf-nfa+1
      nfm=nf+1

      t=t0+(nf-1)*htf

      ax=q(ifom+nf)-2.0d0*q(ifo+nf)+q(ifop+nf)
      ay=q(jfom+nf)-2.0d0*q(ifo+nf)+q(jfop+nf)
      at=q(ifo+nfm)-q(ifo+nf)
      a1=q(iffo+nf)-(coe1*ax+coe2*ay+coe3*at)
      yvar=q(ifo1+nf)
      a1=a1-htf*d1gnl(yvar)*q(ifo+nf)-htf*alfa*(yvar-zf(x,y,t))
      resf(if,jf,nf)=a1
   10 continue
      return
      end
c------------------------------------------------
      subroutine rescal(k,km,ne)
      implicit real*8 (a-h,o-z)
      common/res/resf(129,129,830)
      common/pntrs/ist(129,129),isf(129,129,8),
     X             ism(129,129),isc(129,129,8)
c
      do 30 le=1,ne
      call key(k,le,ist,iif,jjf,nnf,x0,y0,t0,hxf,hyf,htf)
      call key(km,le,ism,iic,jjc,nnc,x0,y0,t0,hxc,hyc,htc)
      do 1 if=1,iif
      do 1 jf=1,jjf
      do 1 nf=1,nnf
    1 resf(if,jf,nf)=0.0d0
      if(le.eq.1) then
         call resca1(k,ne,le)
         call trare1(k,km,ne,le)
      else
         call resca2(k,ne,le)
         call trare2(k,km,ne,le)
      endif
   30 continue
      return
      end
c------------------------------------------------
      subroutine store(k,ne,iwin,dtwin,erry,errp,errz)
c
c  store the solution, compute sol. error and stat. residual
c
      implicit real*8 (a-h,o-z)
      common/qarr/q(50000000)
      common/pntrs/ist(129,129),is(129,129,8),
     X             ism(129,129),isc(129,129,8)
      common/albe/sigma
      common/cweight/rnu,alfa,beta
      common/horiz/nwin
      dimension err(12)

      external fsol,zf,gnl,d1gnl
c
      n2=2*ne
      do 2 l=1,n2
      call key(k,l,ist,ii,jj,nn,x0,y0,t0,hx,hy,ht)
      do 1 i=1,ii
      do 1 j=1,jj
    1 is(i,j,l)=ist(i,j)
    2 continue
      i1=ii-1
      j1=jj-1
      n1=nn-1

      write(*,*) 
      write(*,*) ' ii=    ',ii
      write(*,*) ' jj=    ',jj
      write(*,*) ' nn=    ',nn
      write(*,*) ' h =    ',hx,hy,ht

c
      coe1=sigma*ht/(hx*hx)
      coe2=sigma*ht/(hy*hy)
      coe3=1.0
      coe=2.0d0*(coe1+coe2)+coe3
      write(*,*) ' sigma * gamma = ', coe1
c
c---- static residual
c
c* state
      le=1
      err(le)=0.0d0
      erry   =0.0d0
      errz   =0.0d0
      errt   =0.0d0
c
      do 20 n=2,n1
      t=t0+(n-1)*ht
      nm=n-1
      do 20 i=2,i1
      do 20 j=2,j1
      x=(i-1)*hx
      y=(j-1)*hy

      iom=is(i-1,j,le)
      io =is(i,j,le)
      iop=is(i+1,j,le)
      jom=is(i,j-1,le)
      jop=is(i,j+1,le)
      ifo=is(i,j,ne+le)

      io2 =is(i,j,2)
c
      ax=q(iom+n)-2.0d0*q(io+n)+q(iop+n)
      ay=q(jom+n)-2.0d0*q(io+n)+q(jop+n)
      at=q(io+nm)-q(io+n)
      a1=q(ifo+n)-(coe1*ax+coe2*ay+coe3*at)
      pvar=q(io2+n)
      a1=a1-ht*gnl(q(io+n))+ht*pvar/rnu
      err(le)=err(le)+a1**2

      erry=erry+(fsol(x,y,t,le)-q(io+n))**2
      errz=errz+(zf(x,y,t)-q(io+n))**2
   20 continue
      err(le)=sqrt(err(le)*hx*hy/ht)
      erry=sqrt(erry*hx*hy*ht)
      errz=sqrt(errz*hx*hy*ht)
      write(*,*)' eq.n and static residue ',le,err(le) 
      write(*,*)' eq.n and solution error ',le,erry
      write(*,*)

c* adjoint
      le=2
      err(le)=0.0d0
      errp   =0.0d0
c
      do 40 n=2,n1
      t=t0+(n-1)*ht
      nm=n+1
      do 40 i=2,i1
      do 40 j=2,j1
      x=(i-1)*hx
      y=(j-1)*hy

      iom=is(i-1,j,le)
      io =is(i,j,le)
      iop=is(i+1,j,le)
      jom=is(i,j-1,le)
      jop=is(i,j+1,le)
      ifo=is(i,j,ne+le)

      io1 =is(i,j,1)
c
      ax=q(iom+n)-2.0d0*q(io+n)+q(iop+n)
      ay=q(jom+n)-2.0d0*q(io+n)+q(jop+n)
      at=(q(io+nm)-q(io+n))

      a1=q(ifo+n)-(coe1*ax+coe2*ay+coe3*at)
      yvar=q(io1+n)
      a1=a1-ht*d1gnl(yvar)*q(io+n)-ht*alfa*(yvar-zf(x,y,t))
      err(le)=err(le)+a1**2
      errp=errp+(fsol(x,y,t,le)-q(io+n))**2
   40 continue
      err(le)=sqrt(err(le)*hx*hy/ht)
      errp=sqrt(errp*hx*hy*ht)
      write(*,*)' eq.n and static residue ',le,err(le) 
      write(*,*)' eq.n and solution error ',le,errp

c
c  terminal 
c
      erest   =0.0d0
      n=nn
      t=t0+(n-1)*ht
      do 50 i=2,i1
      do 50 j=2,j1
      x=(i-1)*hx
      y=(j-1)*hy

      io1 =is(i,j,1)
      io2 =is(i,j,2)

      errt  = errt+(zf(x,y,t)-q(io1+n))**2
      erest = erest+(q(io2+n)-beta*(zf(x,y,t)-q(io1+n)))**2
   50 continue
      errt=sqrt(errt*hx*hy)
      erest=sqrt(erest*hx*hy)
      write(*,*)'  terminal eq. res. ',erest
c
      write(*,*)'  tracking error ',errz
      write(*,*)'  terminal error ',errt
c
c     store for matlab
c  
cc      open(unit=17,status='unknown',file='midsol.dat',access='append')   

C  PRINT THE SOLUTION VALUE AT A MESH POINT FOR ALL TIMES
C
      i=(ii+1)/2
      j=(jj+1)/2
      x=(i-1)*hx
      y=(j-1)*hy
      do 95 n=1,nn-1
      t=t0+(n-1)*ht
      cont=q(is(i,j,2)+n)/rnu
 95   write(17,133) t,q(is(i,j,1)+n),zf(x,y,t),cont
 133  format(1x,4(1x,f16.8))
      return
      end
c------------------------------------------------
      subroutine subtrt(kc,kf,ne)
      implicit real*8 (a-h,o-z)
      common/qarr/q(50000000)
      common/pntrs/ist(129,129),isf(129,129,8),
     X             ism(129,129),isc(129,129,8)
      do 20 l=1,ne
      call key(kf,l,ist,iif,jjf,nnf,x0,y0,t0,hxf,hyf,htf)
      call key(kc,l,ism,iic,jjc,nnc,x0,y0,t0,hxc,hyc,htc)
      do 2 i=1,iif
      do 2 j=1,jjf
    2 isf(i,j,l)=ist(i,j)
      do 3 i=1,iic
      do 3 j=1,jjc
    3 isc(i,j,l)=ism(i,j)
      iic1=iic-1
      jjc1=jjc-1
      nnc1=nnc-1

c* no time coarsening 
      do 17 ic=1,iic
      if=2*ic-1
      do 17 jc=1,jjc
      jf=2*jc-1
      ifq=isf(if,jf,l)
      icq=isc(ic,jc,l)
      do 17 nc=1,nnc
      nf=nc
      q(icq+nc)=q(icq+nc)-q(ifq+nf)
   17 continue
   20 continue
      return
      end

c------------------------------------------------------------
      subroutine trare1(k,km,ne,le)
      implicit real*8 (a-h,o-z)
      common/qarr/q(50000000)
      common/pntrs/ist(129,129),isf(129,129,8),
     X             ism(129,129),isc(129,129,8)
      common/res/resf(129,129,830)
      common/albe/sigma
c
      call key(km,le+ne,ism,iic,jjc,nnc,x0,y0,t0,hxc,hyc,htc)
      do 3 i=1,iic
      do 3 j=1,jjc
    3 isc(i,j,le+ne)=ism(i,j)
c
c   half-full weighting of the internal residuals 
c
      iic1=iic-1
      jjc1=jjc-1
      nnc1=nnc-1
c* no time coarsening
      do 75 ic=2,iic1
      if=2*ic-1
      ifm=if-1
      ifp=if+1
      do 75 jc=2,jjc1
      jf=2*jc-1
      jfm=jf-1
      jfp=jf+1
      icf=isc(ic,jc,le+ne)
      do 75 nc=2,nnc
      nf=nc
      nfm=nf-1
      nfp=nf+1
C Half-weighting      
      sq=resf(if-1,jf,nf)+resf(if+1,jf,nf)
     X  +resf(if,jf-1,nf)+resf(if,jf+1,nf)
      q(icf+nc)=0.5*resf(if,jf,nf)+0.125*sq
   75 continue
      return
      end
c------------------------------------------------------------
      subroutine trare2(k,km,ne,le)
      implicit real*8 (a-h,o-z)
      common/qarr/q(50000000)
      common/pntrs/ist(129,129),isf(129,129,8),
     X             ism(129,129),isc(129,129,8)
      common/res/resf(129,129,830)
      common/albe/sigma


      call key(km,le+ne,ism,iic,jjc,nnc,x0,y0,t0,hxc,hyc,htc)
      do 3 i=1,iic
      do 3 j=1,jjc
    3 isc(i,j,le+ne)=ism(i,j)
c
c   half-full weighting of the internal residuals 
c
      iic1=iic-1
      jjc1=jjc-1
      nnc1=nnc-1
c* no time coarsening
      do 75 ic=2,iic1
      if=2*ic-1
      ifm=if-1
      ifp=if+1
      do 75 jc=2,jjc1
      jf=2*jc-1
      jfm=jf-1
      jfp=jf+1
      icf=isc(ic,jc,le+ne)
      do 75 nca=2,nnc
      nc=nnc-nca+1
      nf=nc
      nfm=nf+1
      nfp=nf-1
C Half-weighting      
      sq=resf(if-1,jf,nf)+resf(if+1,jf,nf)
     X  +resf(if,jf-1,nf)+resf(if,jf+1,nf)
      q(icf+nc)=0.5*resf(if,jf,nf)+0.125*sq
   75 continue
      return
      end

c------------------------------------------------------------

      subroutine intr2(km,k,ne)
      implicit real*8 (a-h,o-z)
c
      common/qarr/q(50000000)
      common/pntrs/ist(129,129),isf(129,129,8),
     X             ism(129,129),isc(129,129,8)
      do 30 le=1,ne
      call key(k,le,ist,iif,jjf,nnf,x0,y0,t0,hxf,hyf,htf)
      call key(km,le,ism,iic,jjc,nnc,x0,y0,t0,hxc,hyc,htc)
c
c no time coarsening
c
c
c side i=1
c
      ic=1
      if=2*ic-1
      do 55 jc=2,jjc
      jf=2*jc-1
      ifq=ist(if,jf)
      jfqm=ist(if,jf-1)
      icq=ism(ic,jc)
      jcqm=ism(ic,jc-1)
      do 55 nc=2,nnc
      nf=nc
c
      a=.50d0*(q(icq+nc)+q(jcqm+nc))
      q(ifq+nf)=q(icq+nc)
      q(jfqm+nf)=a
   55 continue
c
c edge i=1,j=1
      ic=1
      if=1
      jc=1 
      jf=1
c corner i=1,j=1,n=1
      ifq=ist(if,jf)
      icq=ism(ic,jc)
      q(ifq+1)=q(icq+1)
c
      do 56 nc=2,nnc
      nf=nc
c
      q(ifq+nf)=q(icq+nc)
   56 continue
c side j=1
      jc=1
      jf=2*jc-1
      do 510 ic=2,iic
      if=2*ic-1
      ifq=ist(if,jf)
      ifqm=ist(if-1,jf)
      icq=ism(ic,jc)
      icqm=ism(ic-1,jc)
      do 510 nc=2,nnc
      nf=nc
c
      am=.50d0*(q(icqm+nc)+q(icq+nc))
      q(ifq+nf)=q(icq+nc)
      q(ifqm+nf)=am
  510 continue
c
c edge j=1,n=1
      jc=1
      jf=1
      nc=1
      nf=1
      do 511 ic=2,iic
      if=2*ic-1
      ifq=ist(if,jf)
      ifqm=ist(if-1,jf)
      icq=ism(ic,jc)
      icqm=ism(ic-1,jc)
c
      am=.50d0*(q(icqm+nc)+q(icq+nc))
      q(ifq+nf)=q(icq+nc)
      q(ifqm+nf)=am
  511 continue
c side n=1
      nc=1
      nf=1
      do 515 ic=2,iic
      if=2*ic-1
      do 515 jc=2,jjc
      jf=2*jc-1
      ifq=ist(if,jf)
      ifqm=ist(if-1,jf)
      jfqm=ist(if,jf-1)
      jfqmm=ist(if-1,jf-2)
      ijfqm=ist(if-1,jf-1)
      icq=ism(ic,jc)
      icqm=ism(ic-1,jc)
      jcqm=ism(ic,jc-1)
      ijcqm=ism(ic-1,jc-1)
c
      a=.50d0*(q(icq+nc)+q(jcqm+nc))
      q(ifq+nf)=q(icq+nc)
      q(jfqm+nf)=a

      a=.50d0*(q(icqm+nc)+q(icq+nc))
      am=.250d0*(q(icq+nc)+q(icqm+nc)
     X          +q(jcqm+nc)+q(ijcqm+nc))
      q(ifqm+nf)=a
      q(ijfqm+nf)=am
  515 continue
c edge n=1,i=1
      nc=1
      nf=1
      ic=1 
      if=1
      do 516 jc=2,jjc
      jf=2*jc-1
      ifq=ist(if,jf)
      jfqm=ist(if,jf-1)
      icq=ism(ic,jc)
      jcqm=ism(ic,jc-1)
      a=.50d0*(q(icq+nc)+q(jcqm+nc))
      q(ifq+nf)=q(icq+nc)
      q(jfqm+nf)=a
  516 continue
c internal points and outer sides
      do 520 ic=2,iic
      if=2*ic-1
      do 520 jc=2,jjc
      jf=2*jc-1
      ifq=ist(if,jf)
      ifqm=ist(if-1,jf)
      jfqm=ist(if,jf-1)
      jfqmm=ist(if-1,jf-2)
      ijfqm=ist(if-1,jf-1)
      icq=ism(ic,jc)
      icqm=ism(ic-1,jc)
      jcqm=ism(ic,jc-1)
      ijcqm=ism(ic-1,jc-1)
      do 520 nc=2,nnc
      nf=nc
c
      a=.50d0*(q(icq+nc)+q(jcqm+nc))
      q(ifq+nf)=q(icq+nc)
      q(jfqm+nf)=a

      a=.50d0*(q(icqm+nc)+q(icq+nc))
      am=.250d0*(q(icqm+nc)+q(icq+nc)+q(jcqm+nc)+q(ijcqm+nc))
      q(ifqm+nf)=a
      q(ijfqm+nf)=am
  520 continue
   30 continue
      return
      end

c------------------------------------------------------------

      subroutine relaxts(k,err,errm,ne)
c---- time-splitted collective Gauss-Seidel iteration
      implicit real*8 (a-h,o-z)
      dimension err(12)
      common/qarr/q(50000000)
      common/pntrs/ist(129,129),iss(129,129,8),
     X             ism(129,129),isc(129,129,8)
      common/albe/sigma
      common/cweight/rnu,alfa,beta
      common/horiz/nwin

      external zf,gnl,d1gnl,d2gnl
c
      n2=2*ne
      do 2 l=1,n2
      call key(k,l,ist,imax,jmax,nmax,x0,y0,t0,hx,hy,ht)
      do 1 i=1,imax
      do 1 j=1,jmax
    1 iss(i,j,l)=ist(i,j)
    2 continue


      i1=imax-1
      j1=jmax-1
      n1=nmax-1
      ht2=ht*ht
      h2xy=hx*hy

      dt=ht
      h2=h2xy
c
      coe1=sigma*ht/(hx*hx)
      coe2=sigma*ht/(hy*hy)
      coe3=1.0
      coe=2.0d0*(coe1+coe2)+coe3
c
c* state & adjoint
      errm=0.0d0
      errq=0.0
      erru=0.0

      le1=1
      le2=2

      err(le1)=0.0 
      err(le2)=0.0

      do 40 i=2,i1
      do 40 j=2,j1

      iom1=iss(i-1,j,le1)
      io1=iss(i,j,le1)
      iop1=iss(i+1,j,le1)
      jom1=iss(i,j-1,le1)
      jop1=iss(i,j+1,le1)
      ifo1=iss(i,j,ne+le1)

      iom2=iss(i-1,j,le2)
      io2=iss(i,j,le2)
      iop2=iss(i+1,j,le2)
      jom2=iss(i,j-1,le2)
      jop2=iss(i,j+1,le2)
      ifo2=iss(i,j,ne+le2)

      x=(i-1)*hx
      y=(j-1)*hy

      do 20 n=2,nmax-1
      nm=n-1
      np=n+1
      t=t0+(n-1)*ht

      na=nmax-n+1
      nam=na-1
      nap=na+1
      ta=t0+(na-1)*ht
c
      ax=q(iom1+n)-2.0d0*q(io1+n)+q(iop1+n)
      ay=q(jom1+n)-2.0d0*q(io1+n)+q(jop1+n)
      at=q(io1+nm)-q(io1+n)
      a1=q(ifo1+n)-(coe1*ax+coe2*ay+coe3*at)
      pvar=q(io2+n)
      a1=a1-ht*gnl(q(io1+n))+ht*pvar/rnu

      axa=q(iom1+na)-2.0d0*q(io1+na)+q(iop1+na)
      aya=q(jom1+na)-2.0d0*q(io1+na)+q(jop1+na)
      ata=q(io1+nam)-q(io1+na)
      a1a=q(ifo1+na)-(coe1*axa+coe2*aya+coe3*ata)
      pvara=q(io2+na)
      a1a=a1a-ht*gnl(q(io1+na))+ht*pvara/rnu
c
      bx=q(iom2+n)-2.0d0*q(io2+n)+q(iop2+n)
      by=q(jom2+n)-2.0d0*q(io2+n)+q(jop2+n)
      bt=q(io2+np)-q(io2+n)
      b1=q(ifo2+n)-(coe1*bx+coe2*by+coe3*bt)
      yvar=q(io1+n)
      b1=b1-ht*d1gnl(yvar)*q(io2+n)-ht*alfa*(yvar-zf(x,y,t))

      bxa=q(iom2+na)-2.0d0*q(io2+na)+q(iop2+na)
      bya=q(jom2+na)-2.0d0*q(io2+na)+q(jop2+na)
      bta=q(io2+nap)-q(io2+na)
      b1a=q(ifo2+na)-(coe1*bxa+coe2*bya+coe3*bta)
      yvara=q(io1+na)
      b1a=b1a-ht*d1gnl(yvara)*q(io2+na)-ht*alfa*(yvara-zf(x,y,ta))
c
      rt1=-(1.0 + 4.0*coe1)+ht*d1gnl(q(io1+n))
      rt2=(ht*ht/rnu)*(alfa+d2gnl(q(io1+n))*q(io2+n))
      rden  = rt1*rt1+rt2
      rnum11= rt1*a1
      rnum12= (ht/rnu)*b1
      rnum1 = rnum11+rnum12
c
      rt1a=-(1.0 + 4.0*coe1)+ht*d1gnl(q(io1+na))
      rt2a=(ht*ht/rnu)*(alfa+d2gnl(q(io1+na))*q(io2+na))
      rdena =rt1a*rt1a+rt2a
      rnum21=-ht*(alfa+d2gnl(q(io1+na))*q(io2+na))*a1a
      rnum22= rt1a*b1a
      rnum2 = rnum21+rnum22

      qx=q(io1+n)+rnum1/rden
      ux=q(io2+na)+rnum2/rdena

      erq=q(io1+n)-qx
      eru=q(io2+na)-ux

      errq=errq+erq*erq
      erru=erru+eru*eru

      q(io1+n)=qx
      q(io2+na)=ux
   20 continue
c
c  terminal condition p(T) = beta*(y(T)-z_T)
c
      n=nmax
      t=t0+(n-1)*ht
      ax=q(iom1+n)-2.0d0*q(io1+n)+q(iop1+n)
      ay=q(jom1+n)-2.0d0*q(io1+n)+q(jop1+n)
      at=q(io1+nm)-q(io1+n)
      a1=q(ifo1+n)-(coe1*ax+coe2*ay+coe3*at)
      a1=a1-ht*gnl(q(io1+n))+ht*q(io2+n)/rnu

      b1=q(io2+n) - beta*(q(io1+n)-zf(x,y,t))

      rden=beta*dt+(1.0+4.0*coe1)*rnu-ht*rnu*d1gnl(q(io1+n))

      qx=q(io1+n)-(rnu*a1-ht*b1)/rden
      ux=q(io2+n)-(beta*rnu*a1
     X           +rnu*(1.0+4.0*coe1-ht*d1gnl(q(io1+n)))*b1)/rden
      q(io1+n)=qx
      q(io2+n)=ux
   40 continue

      errm=0.0
      err(le1)=sqrt(errq*hx*hy*ht) 
      err(le2)=sqrt(erru*hx*hy*ht)/rnu
      errm=err(le1)+err(le2)
      return
      end

c------------------------------------------------------------

      subroutine relaxtl(k,err,errm,ne)
c---- t-line collective Gauss-Seidel iteration    
      implicit real*8 (a-h,o-z)

      dimension err(12)
      common/qarr/q(50000000)
      common/pntrs/ist(129,129),iss(129,129,8),
     X             ism(129,129),isc(129,129,8)
      common/albe/sigma
      common/cweight/rnu,alfa,beta
      common/horiz/nwin
      dimension aa2(2,2,333),bb2(2,2,333),cc2(2,2,333),
     *          ff(2,333),xx2(2,333)

      external zf,gnl,d1gnl,d2gnl
c
      n2=2*ne
      do 2 l=1,n2
      call key(k,l,ist,imax,jmax,nmax,x0,y0,t0,hx,hy,ht)
      do 1 i=1,imax
      do 1 j=1,jmax
    1 iss(i,j,l)=ist(i,j)
    2 continue

      i1=imax-1
      j1=jmax-1
      n1=nmax-1
      ht2=ht*ht
      h2xy=hx*hy

      dt=ht
      h2=h2xy
c
      coe1=sigma*ht/(hx*hx)
      coe2=sigma*ht/(hy*hy)
      coe3=1.0
      coe=2.0d0*(coe1+coe2)+coe3
c
      do 5 n=1,nmax
      do 5 i=1,2
      ff(i,n)=0.
      xx2(i,n)=0.
      do 5 j=1,2
      aa2(i,j,n)=0.
      bb2(i,j,n)=0.
      cc2(i,j,n)=0.
    5 continue
c
c* state & adjoint
c
      le1=1
      le2=2

      erry=0.0d0
      errp=0.0d0

      err(le1)=0.0 
      err(le2)=0.0
c
      do 40 i=2,i1
      do 40 j=2,j1

      x=(i-1)*hx
      y=(j-1)*hy

      io1=iss(i,j,le1)
      io2=iss(i,j,le2)

      im=i-1
      if(i.eq.1) im=i+1
      iom1=iss(im,j,le1)
      iom2=iss(im,j,le2)


      ip=i+1
      if(i.eq.imax) ip=i-1
      iop1=iss(ip,j,le1)
      iop2=iss(ip,j,le2)

      jm=j-1
      if(j.eq.1) jm=j+1
      jom1=iss(i,jm,le1)
      jom2=iss(i,jm,le2)

      jp=j+1
      if(j.eq.jmax) jp=j-1
      jop1=iss(i,jp,le1)
      jop2=iss(i,jp,le2)

      ifo1=iss(i,j,ne+le1)
      ifo2=iss(i,j,ne+le2)
c
c in (0,T)
c
      n=2
      nm=n-1
      np=n+1
      t=t0+(n-1)*dt
cc y-rhs
      ax=q(iom1+n)-2.0d0*q(io1+n)+q(iop1+n)
      ay=q(jom1+n)-2.0d0*q(io1+n)+q(jop1+n)
      at=q(io1+nm)-q(io1+n)
      a1=q(ifo1+n)-(coe1*ax+coe2*ay+coe3*at)
      a1=a1-ht*gnl(q(io1+n))+ht*q(io2+n)/rnu
cc p-rhs
      bx=q(iom2+n)-2.0d0*q(io2+n)+q(iop2+n)
      by=q(jom2+n)-2.0d0*q(io2+n)+q(jop2+n)
      bt=q(io2+np)-q(io2+n)
      b1=q(ifo2+n)-(coe1*bx+coe2*by+coe3*bt)
      b1=b1-ht*d1gnl(q(io1+n))*q(io2+n)-ht*alfa*(q(io1+n)-zf(x,y,t))
cc
      aa2(1,1,n) = -(1.+4.*coe1) + ht*d1gnl(q(io1+n))
      aa2(2,2,n) = -(1.+4.*coe1) + ht*d1gnl(q(io1+n))
      aa2(1,2,n) = -dt/rnu
      aa2(2,1,n) = ht*(alfa + d2gnl(q(io1+n))*q(io2+n))
cc      
      cc2(1,1,n) =  0.0
      cc2(2,2,n) =  1.0
      cc2(1,2,n) =  0.0
      cc2(2,1,n) =  0.0
cc
      ff(1,n) = a1 
      ff(2,n) = b1 
cc
      do 30 n=3,n1
      nm=n-1
      np=n+1
      t=t0+(n-1)*dt
cc y-rhs
      ax=q(iom1+n)-2.0d0*q(io1+n)+q(iop1+n)
      ay=q(jom1+n)-2.0d0*q(io1+n)+q(jop1+n)
      at=q(io1+nm)-q(io1+n)
      a1=q(ifo1+n)-(coe1*ax+coe2*ay+coe3*at)
      a1=a1-ht*gnl(q(io1+n))+ht*q(io2+n)/rnu
cc p-rhs
      bx=q(iom2+n)-2.0d0*q(io2+n)+q(iop2+n)
      by=q(jom2+n)-2.0d0*q(io2+n)+q(jop2+n)
      bt=q(io2+np)-q(io2+n)
      b1=q(ifo2+n)-(coe1*bx+coe2*by+coe3*bt)
      b1=b1-ht*d1gnl(q(io1+n))*q(io2+n)-dt*alfa*(q(io1+n)-zf(x,y,t)) 
cc      
      aa2(1,1,n) = -(1.+4.*coe1) + ht*d1gnl(q(io1+n))
      aa2(2,2,n) = -(1.+4.*coe1) + ht*d1gnl(q(io1+n))
      aa2(1,2,n) = -dt/rnu
      aa2(2,1,n) =  ht*(alfa + d2gnl(q(io1+n))*q(io2+n))
cc      
      bb2(1,1,n) = 1.0
      bb2(2,2,n) = 0.0
      bb2(1,2,n) = 0.0
      bb2(2,1,n) = 0.0
cc      
      cc2(1,1,n) =  0.0
      cc2(2,2,n) =  1.0
      cc2(1,2,n) =  0.0
      cc2(2,1,n) =  0.0

      ff(1,n) = a1 
      ff(2,n) = b1      
   30 continue
c
c  terminal condition p(T)=beta*(y(T)-z_T)      
c
      n=nmax
      nm=n-1
      t=t0+(n-1)*dt
cc y-rhs
      ax=q(iom1+n)-2.0d0*q(io1+n)+q(iop1+n)
      ay=q(jom1+n)-2.0d0*q(io1+n)+q(jop1+n)
      at=q(io1+nm)-q(io1+n)
      a1=q(ifo1+n)-(coe1*ax+coe2*ay+coe3*at)
      a1=a1-ht*gnl(q(io1+n))+ht*q(io2+n)/rnu
cc p-rhs
      b1=q(io2+n) - beta*(q(io1+n)-zf(x,y,t))
cc      
      aa2(1,1,n) = -(1.+4.*coe1) + ht*d1gnl(q(io1+n))
      aa2(2,2,n) = -1.0
      aa2(1,2,n) = -dt/rnu
      aa2(2,1,n) =  beta
cc      
      bb2(1,1,n) = 1.0
      bb2(2,2,n) = 0.0
      bb2(1,2,n) = 0.0
      bb2(2,1,n) = 0.0

      ff(1,n) =  a1
      ff(2,n) =  b1
c
c     block tridiagonal solve
c

      call tridb4(aa2,bb2,cc2,ff,xx2,2,nmax)
      do 20 n=2,nmax
      ery=xx2(1,n)
      erry=erry+ery*ery
      erp=xx2(2,n)
      errp=errp+erp*erp

      q(io1+n) = q(io1+n) + ery
      q(io2+n) = q(io2+n) + erp
   20 continue

   40 continue

      errm=0.0
      err(le1)=sqrt(erry*hx*hy*ht) 
      err(le2)=sqrt(errp*hx*hy*ht)/rnu

      errm=err(le1)+err(le2)
      return
      end

c----------------------------------------------------

      subroutine tridb4(aa,bb,cc,df,x,ni,nf)
      implicit real*8 (a-h,o-z)
      parameter (np=2)
c
c---- solution of ( np x np )-block tridiagonal systems
c     b a c
      dimension aa(np,np,333),bb(np,np,333),cc(np,np,333),
     *          ga(np,np,333),af(np,np),y(np,333),x(np,333),
     *          df(np,333)
      dimension ax(np,np),bx(np,np)

      nip=ni+1
      k=ni
c
c af = aa
      do 5 i=1,np
      do 5 j=1,np
      af(i,j)=aa(i,j,k)
    5 continue
c
c solve af * y = df
c call gaussj(a,n,np,b,m,mp)
      do 50 i=1,np
      bx(i,1)=df(i,k)
      do 50 j=1,np
      ax(i,j)=af(i,j)
   50 continue
      call gaussj(ax,np,np,bx,1,np)
      do 55 i=1,np
      y(i,k)=bx(i,1)
   55 continue

      do 10 k=nip,nf
      k1=k-1
c
c solve af * ga = cc by columns for ga
c call gaussj(a,n,np,b,m,mp)
      do 60 i=1,np
      do 60 j=1,np
      bx(i,j)=cc(i,j,k1)
      ax(i,j)=af(i,j)
   60 continue
      call gaussj(ax,np,np,bx,np,np)
      do 65 i=1,np
      do 65 j=1,np
      ga(i,j,k) = bx(i,j)
   65 continue
c      
c  af = aa - bb * ga
      do 80 i=1,np
      do 80 j=1,np
      af(i,j)=aa(i,j,k)
   80 continue
      do 85 i=1,np
      do 85 j=1,np
      do 85 l=1,np
      af(i,j)=af(i,j)-bb(i,l,k)*ga(l,j,k)
   85 continue
c
c solve af * y = df - b * y
      do 90 i=1,np
      do 90 l=1,np
      df(i,k) = df(i,k)-bb(i,l,k)*y(l,k1)
   90 continue
c
c call gaussj(a,n,np,b,m,mp)
      do 70 i=1,np
      bx(i,1)=df(i,k)
      do 70 j=1,np
      ax(i,j)=af(i,j)
   70 continue
      call gaussj(ax,np,np,bx,1,np)
      do 75 i=1,np
      y(i,k)=bx(i,1)
   75 continue
   10 continue
c
c---- back substitution
      do 2 i=1,np
      x(i,nf)=y(i,nf)
    2 continue
      do 20 k=nf-1,1,-1 
      k1=k+1
      do 4 i=1,np
      x(i,k)=y(i,k)
      do 4 l=1,np
      x(i,k)=x(i,k)-ga(i,l,k1)*x(l,k1)
    4 continue
   20 continue      
      return
      end
c-----------------------------------------------------------------

      SUBROUTINE gaussj(a,n,np,b,m,mp)
      implicit real*8 (a-h,o-z)

      INTEGER m,mp,n,np,NMAX
      PARAMETER (NMAX=1000)
      INTEGER i,icol,irow,j,k,l,ll,indxc(NMAX),indxr(NMAX),ipiv(NMAX)
c      REAL a(np,np),b(np,mp)
      DIMENSION a(np,np),b(np,mp)

c      REAL big,dum,pivinv

      do 11 j=1,n
        ipiv(j)=0
11    continue
      do 22 i=1,n
        big=0.
        do 13 j=1,n
          if(ipiv(j).ne.1)then
            do 12 k=1,n
              if (ipiv(k).eq.0) then
                if (abs(a(j,k)).ge.big)then
                  big=abs(a(j,k))
                  irow=j
                  icol=k
                endif

              else if (ipiv(k).gt.1) then
                pause 'singular matrix in gaussj'
              endif
12          continue
          endif
13      continue
        ipiv(icol)=ipiv(icol)+1
        if (irow.ne.icol) then
          do 14 l=1,n
            dum=a(irow,l)
            a(irow,l)=a(icol,l)
            a(icol,l)=dum
14        continue
          do 15 l=1,m
            dum=b(irow,l)
            b(irow,l)=b(icol,l)
            b(icol,l)=dum
15        continue
        endif
        indxr(i)=irow
        indxc(i)=icol
        if (abs(a(icol,icol)).lt.1.0e-13) 
     X      pause 'singular matrix in gaussj'
      

        pivinv=1./a(icol,icol)
        a(icol,icol)=1.
        do 16 l=1,n
          a(icol,l)=a(icol,l)*pivinv
16      continue
        do 17 l=1,m
          b(icol,l)=b(icol,l)*pivinv
17      continue
        do 21 ll=1,n
          if(ll.ne.icol)then
            dum=a(ll,icol)
            a(ll,icol)=0.
            do 18 l=1,n
              a(ll,l)=a(ll,l)-a(icol,l)*dum
18          continue
            do 19 l=1,m
              b(ll,l)=b(ll,l)-b(icol,l)*dum
19          continue
          endif
21      continue
22    continue
      do 24 l=n,1,-1
        if(indxr(l).ne.indxc(l))then

          do 23 k=1,n
            dum=a(k,indxr(l))
            a(k,indxr(l))=a(k,indxc(l))
            a(k,indxc(l))=dum
23        continue
        endif
24    continue
      return
      END