!   This code was generated by rnpl, a numerical programming language 
!   copyright (c) 1994-1998 by Robert L. Marsa and Matthew W. Choptuik 
c-----------------------------------------------------------------------
c
c     M  E  M  O  R  Y        M  A  N  A  G  E  M  E  N  T
c
c                    R  O  U  T  I  N  E  S
c
c-----------------------------------------------------------------------
c     These routines, in conjunction with the file 'mm.inc' reproduced
c     here:
c
c        integer       maxblk
c        parameter     ( maxblk = 1 000 )
c
c        integer       mmmloc(0:maxblk),
c    *                 mmsize(0:maxblk),
c    *                 mmnext(0:maxblk)
c        logical       mfree(0:maxblk)
c
c        common        / commm / mmmloc, mmsize, mmnext, mfree
c
c     constitute a FORTRAN implementation of the ``First-fit''
c     memory management scheme described in 
c 
c     D. Knuth
c     The Art of Computer Programming
c     Vol. 1, Fundamental Algorithms
c     Section 2.5: Dynamic Storage Allocation
c
c     Invoking routines are responsible for declaring actual storage 
c     array (arena) from which blocks are allocated/deallocated.
c    
c     Example arena declaration:
c
c        integer      memsiz
c        parameter  ( memsiz = 1 000 000 )
c        real*8       q(memsiz)
c
c-----------------------------------------------------------------------
c
c     Sample usage (given above declarations)
c
c     Initialize memory management structure ...
c
c        integer    mmaloc,   mmdeal
c        integer    ptr1,     ptr2
c 
c        call mmini(memsiz)
c
c     Allocate two blocks of storage ...
c
c        ptr1 = mmaloc(1000)
c        if( ptr1 .lt. 0 ) then
c           write(*,*) 'Unable to allocate block of size 1000'
c           stop
c        end if
c        ptr2 = mmaloc(1000)
c        if( ptr2 .lt. 0 ) then
c           write(*,*) 'Unable to allocate block of size 1000'
c           stop
c        end if
c
c     ... do something with the blocks ...
c
c        call subroutine_which_wants_2_vectors(q(ptr1),q(ptr2),1000)
c
c     ... and free them. Note that block size must be supplied to 
c     deallocating routine.
c
c        ptr1 = mmdeal(ptr1,1000)
c        ptr2 = mmdeal(ptr2,1000)
c
c
c-----------------------------------------------------------------------

c=======================================================================
c
c  subroutine  mmini(memsiz)
c
c  I: memsiz             integer        Arena size
c
c  O: None
c
c
c     Initializes memory management structure.
c
c=======================================================================
 
      subroutine mmini(memsiz)
 
         implicit      none
 
C-----------------------------------------------------------------------
C
C     M  E  M  O  R  Y      M  A  N  A  G  E  M  E  N  T
C
C                  S  T  R  U  C  T  U  R  E
C
C-----------------------------------------------------------------------
C
      integer       maxblk
      parameter     ( maxblk = 1 000 )
C
      integer       mmmloc(0:maxblk),
     *              mmsize(0:maxblk),
     *              mmnext(0:maxblk)
      logical       mfree(0:maxblk)
C
      common        / commm / mmmloc, mmsize, mmnext, mfree
C
C----------------------------------------------------------------------
 
         integer       memsiz
 
         call l_ivls(mmmloc(0),0,maxblk+1)
         call l_ivls(mmsize(0),0,maxblk+1)
         call l_ivls(mmnext(0),0,maxblk+1)
         mmnext(0) = 1
         mmmloc(1) = 1
         mmsize(1) = memsiz
 
         return
 
      end
 
c=======================================================================
c
c  integer function mmfree()
c
c  I: None
c
c  O: None
c
c     Returns first available location in list of memory block nodes,
c     -1 if none available.
c
c=======================================================================
 
      integer function mmfree()
 
         implicit     none
 
C-----------------------------------------------------------------------
C
C     M  E  M  O  R  Y      M  A  N  A  G  E  M  E  N  T
C
C                  S  T  R  U  C  T  U  R  E
C
C-----------------------------------------------------------------------
C
      integer       maxblk
      parameter     ( maxblk = 1 000 )
C
      integer       mmmloc(0:maxblk),
     *              mmsize(0:maxblk),
     *              mmnext(0:maxblk)
      logical       mfree(0:maxblk)
C
      common        / commm / mmmloc, mmsize, mmnext, mfree
C
C----------------------------------------------------------------------
 
         integer      p,     q
 
         do mmfree = 0 , maxblk
            mfree(mmfree) = .true.
         end do
         q = 0
 200     continue
            p = mmnext(q)
            if( p .eq. 0 ) go to 300
            mfree(p) = .false.
            q = p
         go to 200
 300     continue
         do mmfree = 1 , maxblk
            if( mfree(mmfree) ) return
         end do
         mmfree = -1
 
         return
 
      end
 
c=======================================================================
c
c  integer function mmaloc(reqsiz)
c
c  I: reqsiz             integer        Requested block size.
c
c  O: None
c
c
c     Allocates memory block of size REQSIZ, returning pointer to first
c     location, or -1 if unable to perform allocation.
c
c=======================================================================
 
      integer function mmaloc(reqsiz)
 
         implicit     none
 
C-----------------------------------------------------------------------
C
C     M  E  M  O  R  Y      M  A  N  A  G  E  M  E  N  T
C
C                  S  T  R  U  C  T  U  R  E
C
C-----------------------------------------------------------------------
C
      integer       maxblk
      parameter     ( maxblk = 1 000 )
C
      integer       mmmloc(0:maxblk),
     *              mmsize(0:maxblk),
     *              mmnext(0:maxblk)
      logical       mfree(0:maxblk)
C
      common        / commm / mmmloc, mmsize, mmnext, mfree
C
C----------------------------------------------------------------------
 
         integer      k, p, q, reqsiz
 
         q = 0
 100     continue
            p = mmnext(q)
            if( p .eq. 0 ) then
               mmaloc = -1
               return
            end if
            if( mmsize(p) .ge. reqsiz ) then
               k = mmsize(p) - reqsiz
               if( k .eq. 0 ) then
                  mmnext(q) = mmnext(p)
               else
                  mmsize(p) = k
               end if
               mmaloc = mmmloc(p) + k
               return
            end if
            q = mmnext(q)
         go to 100
 
         end
 
 
c=======================================================================
c
c  integer function mmdeal(p0loc,p0size)
c
c  I: p0loc              integer        Pointer to block to be freed.
c     p0size             integer        Size of block to be freed.
c
c  O: None
c
c
c        Returns block of storage of length P0SIZE, which begins at
c        location P0LOC to free storage list.
c
c=======================================================================
 
      integer function mmdeal(p0loc,p0size)
 
         implicit    none
 
C-----------------------------------------------------------------------
C
C     M  E  M  O  R  Y      M  A  N  A  G  E  M  E  N  T
C
C                  S  T  R  U  C  T  U  R  E
C
C-----------------------------------------------------------------------
C
      integer       maxblk
      parameter     ( maxblk = 1 000 )
C
      integer       mmmloc(0:maxblk),
     *              mmsize(0:maxblk),
     *              mmnext(0:maxblk)
      logical       mfree(0:maxblk)
C
      common        / commm / mmmloc, mmsize, mmnext, mfree
C
C----------------------------------------------------------------------
 
         integer     mmfree
         integer     p, p0, p0loc, p0size, q
 
         p0 = mmfree()
         if( p0 .lt. 0 ) then
            mmdeal = p0
            return
         else
            mmmloc(p0) = p0loc
            mmsize(p0) = p0size
         end if
         q = 0
 100     continue
            p = mmnext(q)
            if( mmmloc(p) .gt. mmmloc(p0)  .or.  p .eq. 0 ) then
               if( mmmloc(p0) + mmsize(p0) .eq. mmmloc(p) ) then
                  mmnext(p0) = mmnext(p)
                  mmsize(p0) = mmsize(p0) + mmsize(p)
               else
                  mmnext(p0) = p
               end if
               if( mmmloc(q) + mmsize(q) .eq. mmmloc(p0) ) then
                  mmnext(q) = mmnext(p0)
                  mmsize(q) = mmsize(q) + mmsize(p0)
               else
                  mmnext(q) = p0
               end if
               mmdeal = 0
               return
            end if
            q = p
         go to 100
 
      end
 
c=======================================================================
c
c  subroutine mmdump(unit)
c
c  I: unit               integer        Logical unit for dump.
c
c  O: None
c
c
c     Dumps map of free list on UNIT.
c
c=======================================================================
 
      subroutine mmdump(unit)
 
         implicit none
 
C-----------------------------------------------------------------------
C
C     M  E  M  O  R  Y      M  A  N  A  G  E  M  E  N  T
C
C                  S  T  R  U  C  T  U  R  E
C
C-----------------------------------------------------------------------
C
      integer       maxblk
      parameter     ( maxblk = 1 000 )
C
      integer       mmmloc(0:maxblk),
     *              mmsize(0:maxblk),
     *              mmnext(0:maxblk)
      logical       mfree(0:maxblk)
C
      common        / commm / mmmloc, mmsize, mmnext, mfree
C
C----------------------------------------------------------------------
 
         integer      n, p, unit
 
         write(unit,100)
 100     FORMAT(//T10,' <<<     Free Block Structure.     >>>'//
     *            T11,'Block',  T29,'Locations'/)
         p = 0
         n = 0
 200     continue
            p = mmnext(p)
            if( p .eq. 0 ) then
               return
            else
               n = n + 1
               write(unit,300) n, mmmloc(p), mmmloc(p) + mmsize(p) - 1
 300           FORMAT(T12,I3,T25,I7,' ... ',I8)
            end if
         go to 200
 
      end

!-----------------------------------------------------------------------
!     Front ends to latinfo for extracting:
!
!     getlatid(ilat):    Lattice id of ilat'th lattice
!     getlatcs(ilat):    Coordinate system of ilat'th lattice (0..)
!     getrank (ilat):    Rank of ilat'th lat
!     getpshape(ilat):   Pointer into PTRS for shape of lat
!                           ptrs(getpshape(ilat) - 1 + j) := shape(j)
!     getshape(ilat,j):  j'th dim (element of shape) of ilat'th lat
!                           ptrs(getpshape(ilat) - 1 + j) := shape(j)
!     getpbounds(ilat):  Pointer into PTRS for bounds of lat
!                           ptrs(getpbounds(ilat) - 1 + j) := bounds(j)
!     getbounds(ilat,j): j'th component of bounds
!     getpcoord(ilat,j)  Pointer into Q for j'th coords. of ilat'th lat
!     getsize(ilat):     Size of lattice (product reduction of shape)
!-----------------------------------------------------------------------
      integer function getnlat()
         implicit      none
         integer       latinfo
         getnlat = latinfo(0,0,0)
         return
      end

      integer function getlatcs(ilat)
         implicit      none
         integer       latinfo
         integer       ilat
         getlatcs = latinfo(ilat,0,9)
         return
      end

      integer function getlatid(ilat)
         implicit      none
         integer       latinfo
         integer       ilat
         getlatid = latinfo(ilat,0,1)
         return
      end

      integer function getrank(ilat)
         implicit      none
         integer       latinfo
         integer       ilat
         getrank  = latinfo(ilat,0,2)
         return
      end

      integer function getpshape(ilat)
         implicit      none
         integer       latinfo
         integer       ilat
         getpshape = latinfo(ilat,0,3)
         return
      end

      integer function getshape(ilat,j)
         implicit      none
         integer       latinfo
         integer       ilat,    j
         getshape = latinfo(ilat,j,4)
         return
      end

      integer function getpbounds(ilat)
         implicit      none
         integer       latinfo
         integer       ilat
         getpbounds = latinfo(ilat,0,7)
         return
      end

      integer function getbounds(ilat,j)
         implicit      none
         integer       latinfo
         integer       ilat,    j
         getbounds = latinfo(ilat,j,8)
         return
      end
      
      integer function getpcoord(ilat,j)
         implicit      none
         integer       latinfo
         integer       ilat,    j
         getpcoord = latinfo(ilat,j,5)
         return
      end

      integer function getsize(ilat)
         implicit      none
         integer       latinfo
         integer       ilat
         getsize = latinfo(ilat,0,6)
         return
      end

!-----------------------------------------------------------------------
!     Main decoding routine
!
!     Requires 'l_ivprod' for product reduction of integer vector.
!
!     option = 0:   Number of lattices
!     option = 1:   Lattice id
!     option = 2:   Rank of ilat'th lattice
!     option = 3:   Pointer into ptrs for shape of ilat'th lattice
!     option = 4:   j'th dim of ilat'th lattice
!     option = 5:   Pointer into q for coordinates of j'th dim of
!                   ilat'th lattice
!     option = 6:   Size of lattice (product reduction of shape)
!     option = 7:   Pointer into ptrs for bounds of ilat'th lattice
!     option = 8:   j'th bound of ilat'th lattice
!     option = 9:   Coordinate system
!-----------------------------------------------------------------------
      integer function latinfo(ilat,j,option)
         implicit      none

         include      'sys_param.inc'
         include      'gfuni0.inc'

         integer       l_ivprod

         integer       ilat,       j,       option, cs
         integer       nlat,       base,    latid,  rank
         integer  shape, bounds, cptr
         
         nlat = (nptrs - 2*ngfcn) / sizelatcb
         if( option .eq. 0 ) then
            latinfo = nlat
            return
         end if

         if( 1 .le. ilat  .and.  ilat .le. nlat ) then
            base    = 2*ngfcn + (ilat - 1)*sizelatcb + 1
            latid   = ptrs(base)
            cs      = ptrs(base+1)
            rank    = ptrs(base+2)
            shape   = base + 3
            bounds  = shape + maxrank
            cptr    = bounds + 2*maxrank
            if(      option .eq. 1 ) then
               latinfo = latid
            else if( option .eq. 2 ) then
               latinfo = rank
            else if( option .eq. 3 ) then
               latinfo = shape
            else if( option .eq. 4 ) then
               if( 1 .le. j  .and.  j .le. rank ) then
                  latinfo = ptrs(shape + j - 1)
               else
                  write(0,1000) j, rank, ilat
1000              format(' latinfo: j = ',i2,' out of bounds for ',
     &                   'rank ',i2,' lattice ',i2)
                  latinfo = -1
               end if
            else if( option .eq. 5 ) then
               if( 1 .le. j  .and.  j .le. rank ) then
                  latinfo = ptrs(cptr + j-1)
               else
                  write(0,1000) j, rank, ilat
                  latinfo = -1
               end if
            else if( option .eq. 6 ) then
                  latinfo = l_ivprod(ptrs(shape),rank)
            else if( option .eq. 7 ) then
                  latinfo = bounds
            else if( option .eq. 8 ) then
               if( 1 .le. j  .and.  j .le. 2*rank ) then
                  latinfo = ptrs(bounds + j - 1)
               else
                  write(0,1000) j, rank, ilat
                  latinfo = -1
               end if
            else if( option .eq. 9 ) then
                  latinfo = cs
            end if
         else
            write(0,1100) ilat, nlat
1100        format(' latinfo: ilat = ',i2,' out of bounds 1 - ',i2)
            latinfo = -1
         end if

         return

      end

!----------------------------------------------------------------------
!  init_coord_difs
!----------------------------------------------------------------------
      subroutine init_coord_difs()
        implicit none

        include 'sys_param.inc'
        include 'gfuni0.inc'
        include 'other_glbs.inc'
        include 'globals.inc'

10      format (a,a,i5,a)
11      format (a,i5,a)
        if(level .gt. 0 .and. mod(Nx0,2) .eq. 1) then
          write(*,11)'WARNING: Nx0=',Nx0,' is odd and level != 0'
          write(*,11)'WARNING: Adjusting Nx0 to ',Nx0-1,
     &      '.  If this is part of a convergence'
          write(*,10)'WARNING: test, rerun level 0 calculation with'
     &      ,' Nx0=',Nx0-1,'.'
          Nx0 = Nx0 - 1
        end if
        Nx = Nx0 * 2**level + 1
        if(level .gt. 0 .and. mod(Ny0,2) .eq. 1) then
          write(*,11)'WARNING: Ny0=',Ny0,' is odd and level != 0'
          write(*,11)'WARNING: Adjusting Ny0 to ',Ny0-1,
     &      '.  If this is part of a convergence'
          write(*,10)'WARNING: test, rerun level 0 calculation with'
     &      ,' Ny0=',Ny0-1,'.'
          Ny0 = Ny0 - 1
        end if
        Ny = Ny0 * 2**level + 1
        dx=(xmax - xmin)/(Nx - 1)
        dy=(ymax - ymin)/(Ny - 1)
        dt=lambda*sqrt((dx*dx+dy*dy)/2)

        return
      end

!----------------------------------------------------------------------
!  init_lats
!----------------------------------------------------------------------
      subroutine init_lats()
        implicit none

        include 'sys_param.inc'
        include 'gfuni0.inc'
        include 'other_glbs.inc'
        include 'globals.inc'

        integer i,mmaloc

        cname(1)='x|y'
        ncnames = 6
        ngfcn = 5
        nptrs = sizelatcb * 1 + 2 * ngfcn
        ptrs(11)=1
        ptrs(12)=1
        ptrs(13)=2
        ptrs(14)=Nx-1+1
        ptrs(15)=Ny-1+1
        ptrs(16)=1
        ptrs(17)=1
        ptrs(18)=1
        ptrs(19)=1
        ptrs(20)=Nx
        ptrs(21)=1
        ptrs(22)=Ny
        ptrs(23)=1
        ptrs(24)=1
        ptrs(25)=1
        ptrs(26)=1
        ptrs(27)=1
        ptrs(28)=1
        x=mmaloc(getsize(1))
        if(x .eq. -1) then
          write(*,*)'ERROR: unable to allocate coordinate'
          write(*,*)'increase memsiz and recompile'
          stop
        end if
        call rdvramp(q(x),getshape(1,1),xmin,dx)
        ptrs(29)=x
        y=x+getshape(1,1)
        call rdvramp(q(y),getshape(1,2),ymin,dy)
        ptrs(30)=y
        return
      end

!----------------------------------------------------------------------
!  init_gfuncs
!----------------------------------------------------------------------
      subroutine init_gfuncs()
        implicit none

        integer  mmaloc
        character*(128)  temp
        character*256  catsqz
        include 'sys_param.inc'
        include 'gfuni0.inc'
        include 'other_glbs.inc'
        include 'globals.inc'

901     format (a,i1)
902     format (a,i1,a)
        ptrs(ngfcn + 1) = 2*ngfcn + 0*sizelatcb + 1
        ptrs(1) = mmaloc(getsize(1))
        if(ptrs(1) .eq. -1) then
          write(*,*)'ERROR: unable to allocate grid function'
          write(*,*)'increase memsiz and recompile'
          stop
        end if
        ptrs(ngfcn + 2) = 2*ngfcn + 0*sizelatcb + 1
        ptrs(2) = mmaloc(getsize(1))
        if(ptrs(2) .eq. -1) then
          write(*,*)'ERROR: unable to allocate grid function'
          write(*,*)'increase memsiz and recompile'
          stop
        end if
        write(temp,901) 'u_',level
        attrfname(1) = catsqz(tag,temp)
        gf_ln(1)=2
        gf_st(1)=0
        ptrs(ngfcn + 3) = 2*ngfcn + 0*sizelatcb + 1
        ptrs(3) = mmaloc(getsize(1))
        if(ptrs(3) .eq. -1) then
          write(*,*)'ERROR: unable to allocate grid function'
          write(*,*)'increase memsiz and recompile'
          stop
        end if
        write(temp,901) 'rhs_',level
        attrfname(2) = catsqz(tag,temp)
        gf_ln(2)=1
        gf_st(2)=2
        ptrs(ngfcn + 4) = 2*ngfcn + 0*sizelatcb + 1
        ptrs(4) = mmaloc(getsize(1))
        if(ptrs(4) .eq. -1) then
          write(*,*)'ERROR: unable to allocate grid function'
          write(*,*)'increase memsiz and recompile'
          stop
        end if
        ptrs(ngfcn + 5) = 2*ngfcn + 0*sizelatcb + 1
        ptrs(5) = mmaloc(getsize(1))
        if(ptrs(5) .eq. -1) then
          write(*,*)'ERROR: unable to allocate grid function'
          write(*,*)'increase memsiz and recompile'
          stop
        end if
        write(temp,901) 'ures_',level
        attrfname(3) = catsqz(tag,temp)
        gf_ln(3)=2
        gf_st(3)=3
        return
      end

!----------------------------------------------------------------------
!  swap_levels
!----------------------------------------------------------------------
      subroutine swap_levels()
        implicit none

        include 'sys_param.inc'
        include 'gfuni0.inc'
        include 'globals.inc'

        integer tmp

        tmp=ptrs(gf_st(1)+2)
        ptrs(gf_st(1)+2)=ptrs(gf_st(1)+1)
        ptrs(gf_st(1)+1)=tmp
        tmp=ptrs(gf_st(3)+2)
        ptrs(gf_st(3)+2)=ptrs(gf_st(3)+1)
        ptrs(gf_st(3)+1)=tmp
        return
      end

!----------------------------------------------------------------------
!  initial_guess
!----------------------------------------------------------------------
      subroutine initial_guess()
        implicit none

        include 'sys_param.inc'
        include 'gfuni0.inc'
        include 'other_glbs.inc'
        include 'globals.inc'

        integer i,j,k
        call rdvcpy(q(ptrs(gf_st(1)+1)), q(ptrs(gf_st(1)+2)),getsize(1))
        call rdvcpy(q(ptrs(gf_st(3)+1)), q(ptrs(gf_st(3)+2)),getsize(1))
        return
      end

!----------------------------------------------------------------------
!  calc_resid
!----------------------------------------------------------------------
      real*8 function calc_resid()
        implicit none

        include 'sys_param.inc'
        include 'gfuni0.inc'
        include 'globals.inc'

        real*8   u,l2norm
        calc_resid=u
        return
      end

!----------------------------------------------------------------------
!  take_step
!----------------------------------------------------------------------
      subroutine take_step(t)
        implicit none

        include 'sys_param.inc'
        include 'gfuni0.inc'
        include 'other_glbs.inc'
        include 'globals.inc'

        integer mmaloc, mmdeal
        real*8  t

        call update0(q(ptrs(gf_st(1)+1)),q(ptrs(gf_st(1)+2)),
     &    q(ptrs(gf_st(2)+1)),ptrs(getpshape(1)+0),ptrs(getpshape(1)+1),
     &    t,q(x),q(y),dt,dx,dy,mgeps,mgmaxcyc)
        call update1(q(ptrs(gf_st(1)+1)),q(ptrs(gf_st(1)+2)),
     &    q(ptrs(gf_st(3)+2)),ptrs(getpshape(1)+0),ptrs(getpshape(1)+1),
     &    dt,dx,dy)
        return
      end

      integer function update_gfuncs(t)
        implicit none
        include 'sys_param.inc'
        include 'gfuni0.inc'
        include 'globals.inc'
        
        real*8   t
        
        call take_step(t)
        call swap_levels
        update_gfuncs=1
        return
      end
      
      

!----------------------------------------------------------------------
!  init_params_attribs
!----------------------------------------------------------------------
      subroutine init_params_attribs()
        implicit none

        include 'sys_param.inc'
        include 'gfuni0.inc'
        include 'other_glbs.inc'
        include 'globals.inc'

        set_amp=0
        amp=(1D0)
        set_camp=0
        camp=(0D0)
        set_epsiter=0
        epsiter=(1D-05)
        set_epsiterid=0
        epsiterid=(1D-05)
        set_fout=0
        fout=0
        set_in_file=0
        set_iter=0
        iter=100
        set_lambda=0
        lambda=(0.5D0)
        set_level=0
        level=0
        set_maxstep=0
        maxstep=50
        set_maxstepid=0
        maxstepid=50
        set_mgeps=0
        mgeps=(1D-09)
        set_mgmaxcyc=0
        mgmaxcyc=50
        set_mpi=0
        mpi=(3.141592653589793D0)
        set_Nx0=0
        Nx0=3
        set_Ny0=0
        Ny0=3
        set_out_file=0
        set_output=0
        output(1)=1
        output(2)=-1
        output(3)=-1
        output(4)=1
        output(5)=0
        set_s_step=0
        s_step=0
        set_ser=0
        ser=0
        set_start_t=0
        start_t=(0D0)
        set_tag=0
        tag=' '
        set_trace=0
        trace(1)=1
        trace(2)=-1
        trace(3)=-1
        trace(4)=1
        trace(5)=0
        set_xc=0
        xc=(0.5D0)
        set_xmax=0
        xmax=(1D0)
        set_xmin=0
        xmin=(0D0)
        set_xwid=0
        xwid=(0.05D0)
        set_yc=0
        yc=(0.5D0)
        set_ymax=0
        ymax=(1D0)
        set_ymin=0
        ymin=(0D0)
        set_ywid=0
        ywid=(0.05D0)
        set_out_gf=0
        out_gf(1)=1
        out_gf(2)=0
        out_gf(3)=1
        return
      end

!----------------------------------------------------------------------
!  read_parameters
!----------------------------------------------------------------------
      subroutine read_parameters(p_file)
        implicit none

        character*(*) p_file
        integer ret
        integer get_int_param, get_real_param
        integer get_str_param, get_ivec_param
        include 'sys_param.inc'
        include 'gfuni0.inc'
        include 'other_glbs.inc'
        include 'globals.inc'

      if(set_amp.ne.1) then
        ret=get_real_param(p_file,'amp',amp,1)
        if(ret.eq.1) then
          set_amp=1
        end if
      end if
      if(set_camp.ne.1) then
        ret=get_real_param(p_file,'camp',camp,1)
        if(ret.eq.1) then
          set_camp=1
        end if
      end if
      if(set_epsiter.ne.1) then
        ret=get_real_param(p_file,'epsiter',epsiter,1)
        if(ret.eq.1) then
          set_epsiter=1
        end if
      end if
      if(set_epsiterid.ne.1) then
        ret=get_real_param(p_file,'epsiterid',epsiterid,1)
        if(ret.eq.1) then
          set_epsiterid=1
        end if
      end if
      if(set_fout.ne.1) then
        ret=get_int_param(p_file,'fout',fout,1)
        if(ret.eq.1) then
          set_fout=1
        end if
      end if
      if(set_in_file.ne.1) then
        ret=get_str_param(p_file,'in_file',in_file,1)
        if(ret.eq.1) then
          set_in_file=1
        end if
      end if
      if(set_iter.ne.1) then
        ret=get_int_param(p_file,'iter',iter,1)
        if(ret.eq.1) then
          set_iter=1
        end if
      end if
      if(set_lambda.ne.1) then
        ret=get_real_param(p_file,'lambda',lambda,1)
        if(ret.eq.1) then
          set_lambda=1
        end if
      end if
      if(set_level.ne.1) then
        ret=get_int_param(p_file,'level',level,1)
        if(ret.eq.1) then
          set_level=1
        end if
      end if
      if(set_maxstep.ne.1) then
        ret=get_int_param(p_file,'maxstep',maxstep,1)
        if(ret.eq.1) then
          set_maxstep=1
        end if
      end if
      if(set_maxstepid.ne.1) then
        ret=get_int_param(p_file,'maxstepid',maxstepid,1)
        if(ret.eq.1) then
          set_maxstepid=1
        end if
      end if
      if(set_mgeps.ne.1) then
        ret=get_real_param(p_file,'mgeps',mgeps,1)
        if(ret.eq.1) then
          set_mgeps=1
        end if
      end if
      if(set_mgmaxcyc.ne.1) then
        ret=get_int_param(p_file,'mgmaxcyc',mgmaxcyc,1)
        if(ret.eq.1) then
          set_mgmaxcyc=1
        end if
      end if
      if(set_mpi.ne.1) then
        ret=get_real_param(p_file,'mpi',mpi,1)
        if(ret.eq.1) then
          set_mpi=1
        end if
      end if
      if(set_Nx0.ne.1) then
        ret=get_int_param(p_file,'Nx0',Nx0,1)
        if(ret.eq.1) then
          set_Nx0=1
        end if
      end if
      if(set_Ny0.ne.1) then
        ret=get_int_param(p_file,'Ny0',Ny0,1)
        if(ret.eq.1) then
          set_Ny0=1
        end if
      end if
      if(set_out_file.ne.1) then
        ret=get_str_param(p_file,'out_file',out_file,1)
        if(ret.eq.1) then
          set_out_file=1
        end if
      end if
      if(set_output.ne.1) then
        ret=get_ivec_param(p_file,'output',output,41)
        if(ret.eq.1) then
          set_output=1
        end if
      end if
      if(set_s_step.ne.1) then
        ret=get_int_param(p_file,'s_step',s_step,1)
        if(ret.eq.1) then
          set_s_step=1
        end if
      end if
      if(set_ser.ne.1) then
        ret=get_int_param(p_file,'ser',ser,1)
        if(ret.eq.1) then
          set_ser=1
        end if
      end if
      if(set_start_t.ne.1) then
        ret=get_real_param(p_file,'start_t',start_t,1)
        if(ret.eq.1) then
          set_start_t=1
        end if
      end if
      if(set_tag.ne.1) then
        ret=get_str_param(p_file,'tag',tag,1)
        if(ret.eq.1) then
          set_tag=1
        end if
      end if
      if(set_trace.ne.1) then
        ret=get_ivec_param(p_file,'trace',trace,41)
        if(ret.eq.1) then
          set_trace=1
        end if
      end if
      if(set_xc.ne.1) then
        ret=get_real_param(p_file,'xc',xc,1)
        if(ret.eq.1) then
          set_xc=1
        end if
      end if
      if(set_xmax.ne.1) then
        ret=get_real_param(p_file,'xmax',xmax,1)
        if(ret.eq.1) then
          set_xmax=1
        end if
      end if
      if(set_xmin.ne.1) then
        ret=get_real_param(p_file,'xmin',xmin,1)
        if(ret.eq.1) then
          set_xmin=1
        end if
      end if
      if(set_xwid.ne.1) then
        ret=get_real_param(p_file,'xwid',xwid,1)
        if(ret.eq.1) then
          set_xwid=1
        end if
      end if
      if(set_yc.ne.1) then
        ret=get_real_param(p_file,'yc',yc,1)
        if(ret.eq.1) then
          set_yc=1
        end if
      end if
      if(set_ymax.ne.1) then
        ret=get_real_param(p_file,'ymax',ymax,1)
        if(ret.eq.1) then
          set_ymax=1
        end if
      end if
      if(set_ymin.ne.1) then
        ret=get_real_param(p_file,'ymin',ymin,1)
        if(ret.eq.1) then
          set_ymin=1
        end if
      end if
      if(set_ywid.ne.1) then
        ret=get_real_param(p_file,'ywid',ywid,1)
        if(ret.eq.1) then
          set_ywid=1
        end if
      end if
        return
      end

!----------------------------------------------------------------------
!  sread_parameters
!----------------------------------------------------------------------
      subroutine sread_parameters(p_str)
        implicit none

        character*(*) p_str
        integer ret
        integer sget_int_param, sget_real_param
        integer sget_str_param, sget_ivec_param
        include 'sys_param.inc'
        include 'gfuni0.inc'
        include 'other_glbs.inc'
        include 'globals.inc'

        if(set_amp.ne.1) then
          ret=sget_real_param(p_str,'amp',amp,1)
          if(ret.eq.1) then
            set_amp=1
          end if
        end if
        if(set_camp.ne.1) then
          ret=sget_real_param(p_str,'camp',camp,1)
          if(ret.eq.1) then
            set_camp=1
          end if
        end if
        if(set_epsiter.ne.1) then
          ret=sget_real_param(p_str,'epsiter',epsiter,1)
          if(ret.eq.1) then
            set_epsiter=1
          end if
        end if
        if(set_epsiterid.ne.1) then
          ret=sget_real_param(p_str,'epsiterid',epsiterid,1)
          if(ret.eq.1) then
            set_epsiterid=1
          end if
        end if
        if(set_fout.ne.1) then
          ret=sget_int_param(p_str,'fout',fout,1)
          if(ret.eq.1) then
            set_fout=1
          end if
        end if
        if(set_in_file.ne.1) then
          ret=sget_str_param(p_str,'in_file',in_file,1)
          if(ret.eq.1) then
            set_in_file=1
          end if
        end if
        if(set_iter.ne.1) then
          ret=sget_int_param(p_str,'iter',iter,1)
          if(ret.eq.1) then
            set_iter=1
          end if
        end if
        if(set_lambda.ne.1) then
          ret=sget_real_param(p_str,'lambda',lambda,1)
          if(ret.eq.1) then
            set_lambda=1
          end if
        end if
        if(set_level.ne.1) then
          ret=sget_int_param(p_str,'level',level,1)
          if(ret.eq.1) then
            set_level=1
          end if
        end if
        if(set_maxstep.ne.1) then
          ret=sget_int_param(p_str,'maxstep',maxstep,1)
          if(ret.eq.1) then
            set_maxstep=1
          end if
        end if
        if(set_maxstepid.ne.1) then
          ret=sget_int_param(p_str,'maxstepid',maxstepid,1)
          if(ret.eq.1) then
            set_maxstepid=1
          end if
        end if
        if(set_mgeps.ne.1) then
          ret=sget_real_param(p_str,'mgeps',mgeps,1)
          if(ret.eq.1) then
            set_mgeps=1
          end if
        end if
        if(set_mgmaxcyc.ne.1) then
          ret=sget_int_param(p_str,'mgmaxcyc',mgmaxcyc,1)
          if(ret.eq.1) then
            set_mgmaxcyc=1
          end if
        end if
        if(set_mpi.ne.1) then
          ret=sget_real_param(p_str,'mpi',mpi,1)
          if(ret.eq.1) then
            set_mpi=1
          end if
        end if
        if(set_Nx0.ne.1) then
          ret=sget_int_param(p_str,'Nx0',Nx0,1)
          if(ret.eq.1) then
            set_Nx0=1
          end if
        end if
        if(set_Ny0.ne.1) then
          ret=sget_int_param(p_str,'Ny0',Ny0,1)
          if(ret.eq.1) then
            set_Ny0=1
          end if
        end if
        if(set_out_file.ne.1) then
          ret=sget_str_param(p_str,'out_file',out_file,1)
          if(ret.eq.1) then
            set_out_file=1
          end if
        end if
        if(set_output.ne.1) then
          ret=sget_ivec_param(p_str,'output',output,41)
          if(ret.eq.1) then
            set_output=1
          end if
        end if
        if(set_s_step.ne.1) then
          ret=sget_int_param(p_str,'s_step',s_step,1)
          if(ret.eq.1) then
            set_s_step=1
          end if
        end if
        if(set_ser.ne.1) then
          ret=sget_int_param(p_str,'ser',ser,1)
          if(ret.eq.1) then
            set_ser=1
          end if
        end if
        if(set_start_t.ne.1) then
          ret=sget_real_param(p_str,'start_t',start_t,1)
          if(ret.eq.1) then
            set_start_t=1
          end if
        end if
        if(set_tag.ne.1) then
          ret=sget_str_param(p_str,'tag',tag,1)
          if(ret.eq.1) then
            set_tag=1
          end if
        end if
        if(set_trace.ne.1) then
          ret=sget_ivec_param(p_str,'trace',trace,41)
          if(ret.eq.1) then
            set_trace=1
          end if
        end if
        if(set_xc.ne.1) then
          ret=sget_real_param(p_str,'xc',xc,1)
          if(ret.eq.1) then
            set_xc=1
          end if
        end if
        if(set_xmax.ne.1) then
          ret=sget_real_param(p_str,'xmax',xmax,1)
          if(ret.eq.1) then
            set_xmax=1
          end if
        end if
        if(set_xmin.ne.1) then
          ret=sget_real_param(p_str,'xmin',xmin,1)
          if(ret.eq.1) then
            set_xmin=1
          end if
        end if
        if(set_xwid.ne.1) then
          ret=sget_real_param(p_str,'xwid',xwid,1)
          if(ret.eq.1) then
            set_xwid=1
          end if
        end if
        if(set_yc.ne.1) then
          ret=sget_real_param(p_str,'yc',yc,1)
          if(ret.eq.1) then
            set_yc=1
          end if
        end if
        if(set_ymax.ne.1) then
          ret=sget_real_param(p_str,'ymax',ymax,1)
          if(ret.eq.1) then
            set_ymax=1
          end if
        end if
        if(set_ymin.ne.1) then
          ret=sget_real_param(p_str,'ymin',ymin,1)
          if(ret.eq.1) then
            set_ymin=1
          end if
        end if
        if(set_ywid.ne.1) then
          ret=sget_real_param(p_str,'ywid',ywid,1)
          if(ret.eq.1) then
            set_ywid=1
          end if
        end if
        return
      end

!----------------------------------------------------------------------
!  read_attributes
!----------------------------------------------------------------------
      subroutine read_attributes(p_file)
        implicit none

        character*(*) p_file
        integer ret
        integer get_int_param, get_real_param
        integer get_str_param, get_ivec_param
        include 'sys_param.inc'
        include 'gfuni0.inc'
        include 'other_glbs.inc'
        include 'globals.inc'

        if(set_out_gf.ne.1) then
          ret=get_int_param(p_file,'out_gf',out_gf,3)
          if(ret.eq.1) then
            set_out_gf=1
          end if
        end if
        return
      end

!----------------------------------------------------------------------
!  check_params_attribs
!----------------------------------------------------------------------
      integer function check_params_attribs()
        implicit none

        integer all_ok
        include 'sys_param.inc'
        include 'gfuni0.inc'
        include 'other_glbs.inc'
        include 'globals.inc'

        all_ok=1
        if(set_amp.eq.0 .and. (1).eq.0) then
          write(*,*) 'ERROR: parameter amp has not been set.'
          all_ok=0
        else if(set_amp.eq.0 .and. (1).eq.1) then
          write(*,*) 'WARNING: using default for parameter amp.'
        end if
        if(set_camp.eq.0 .and. (1).eq.0) then
          write(*,*) 'ERROR: parameter camp has not been set.'
          all_ok=0
        else if(set_camp.eq.0 .and. (1).eq.1) then
          write(*,*) 'WARNING: using default for parameter camp.'
        end if
        if(set_epsiter.eq.0 .and. (1).eq.0) then
          write(*,*) 'ERROR: parameter epsiter has not been set.'
          all_ok=0
        else if(set_epsiter.eq.0 .and. (1).eq.1) then
          write(*,*) 'WARNING: using default for parameter epsiter.'
        end if
        if(set_epsiterid.eq.0 .and. (1).eq.0) then
          write(*,*) 'ERROR: parameter epsiterid has not been set.'
          all_ok=0
        else if(set_epsiterid.eq.0 .and. (1).eq.1) then
          write(*,*) 'WARNING: using default for parameter epsiterid.'
        end if
        if(set_fout.eq.0 .and. (1).eq.0) then
          write(*,*) 'ERROR: parameter fout has not been set.'
          all_ok=0
        else if(set_fout.eq.0 .and. (1).eq.1) then
          write(*,*) 'WARNING: using default for parameter fout.'
        end if
        if(set_in_file.eq.0 .and. (0).eq.0) then
          write(*,*) 'ERROR: parameter in_file has not been set.'
          all_ok=0
        else if(set_in_file.eq.0 .and. (0).eq.1) then
          write(*,*) 'WARNING: using default for parameter in_file.'
        end if
        if(set_iter.eq.0 .and. (1).eq.0) then
          write(*,*) 'ERROR: parameter iter has not been set.'
          all_ok=0
        else if(set_iter.eq.0 .and. (1).eq.1) then
          write(*,*) 'WARNING: using default for parameter iter.'
        end if
        if(set_lambda.eq.0 .and. (1).eq.0) then
          write(*,*) 'ERROR: parameter lambda has not been set.'
          all_ok=0
        else if(set_lambda.eq.0 .and. (1).eq.1) then
          write(*,*) 'WARNING: using default for parameter lambda.'
        end if
        if(set_level.eq.0 .and. (1).eq.0) then
          write(*,*) 'ERROR: parameter level has not been set.'
          all_ok=0
        else if(set_level.eq.0 .and. (1).eq.1) then
          write(*,*) 'WARNING: using default for parameter level.'
        end if
        if(set_maxstep.eq.0 .and. (1).eq.0) then
          write(*,*) 'ERROR: parameter maxstep has not been set.'
          all_ok=0
        else if(set_maxstep.eq.0 .and. (1).eq.1) then
          write(*,*) 'WARNING: using default for parameter maxstep.'
        end if
        if(set_maxstepid.eq.0 .and. (1).eq.0) then
          write(*,*) 'ERROR: parameter maxstepid has not been set.'
          all_ok=0
        else if(set_maxstepid.eq.0 .and. (1).eq.1) then
          write(*,*) 'WARNING: using default for parameter maxstepid.'
        end if
        if(set_mgeps.eq.0 .and. (1).eq.0) then
          write(*,*) 'ERROR: parameter mgeps has not been set.'
          all_ok=0
        else if(set_mgeps.eq.0 .and. (1).eq.1) then
          write(*,*) 'WARNING: using default for parameter mgeps.'
        end if
        if(set_mgmaxcyc.eq.0 .and. (1).eq.0) then
          write(*,*) 'ERROR: parameter mgmaxcyc has not been set.'
          all_ok=0
        else if(set_mgmaxcyc.eq.0 .and. (1).eq.1) then
          write(*,*) 'WARNING: using default for parameter mgmaxcyc.'
        end if
        if(set_mpi.eq.0 .and. (1).eq.0) then
          write(*,*) 'ERROR: parameter mpi has not been set.'
          all_ok=0
        else if(set_mpi.eq.0 .and. (1).eq.1) then
          write(*,*) 'WARNING: using default for parameter mpi.'
        end if
        if(set_Nx0.eq.0 .and. (1).eq.0) then
          write(*,*) 'ERROR: parameter Nx0 has not been set.'
          all_ok=0
        else if(set_Nx0.eq.0 .and. (1).eq.1) then
          write(*,*) 'WARNING: using default for parameter Nx0.'
        end if
        if(set_Ny0.eq.0 .and. (1).eq.0) then
          write(*,*) 'ERROR: parameter Ny0 has not been set.'
          all_ok=0
        else if(set_Ny0.eq.0 .and. (1).eq.1) then
          write(*,*) 'WARNING: using default for parameter Ny0.'
        end if
        if(set_out_file.eq.0 .and. (0).eq.0) then
          write(*,*) 'ERROR: parameter out_file has not been set.'
          all_ok=0
        else if(set_out_file.eq.0 .and. (0).eq.1) then
          write(*,*) 'WARNING: using default for parameter out_file.'
        end if
        if(set_output.eq.0 .and. (1).eq.0) then
          write(*,*) 'ERROR: parameter output has not been set.'
          all_ok=0
        else if(set_output.eq.0 .and. (1).eq.1) then
          write(*,*) 'WARNING: using default for parameter output.'
        end if
        if(set_s_step.eq.0 .and. (1).eq.0) then
          write(*,*) 'ERROR: parameter s_step has not been set.'
          all_ok=0
        else if(set_s_step.eq.0 .and. (1).eq.1) then
          write(*,*) 'WARNING: using default for parameter s_step.'
        end if
        if(set_ser.eq.0 .and. (1).eq.0) then
          write(*,*) 'ERROR: parameter ser has not been set.'
          all_ok=0
        else if(set_ser.eq.0 .and. (1).eq.1) then
          write(*,*) 'WARNING: using default for parameter ser.'
        end if
        if(set_start_t.eq.0 .and. (1).eq.0) then
          write(*,*) 'ERROR: parameter start_t has not been set.'
          all_ok=0
        else if(set_start_t.eq.0 .and. (1).eq.1) then
          write(*,*) 'WARNING: using default for parameter start_t.'
        end if
        if(set_tag.eq.0 .and. (1).eq.0) then
          write(*,*) 'ERROR: parameter tag has not been set.'
          all_ok=0
        else if(set_tag.eq.0 .and. (1).eq.1) then
          write(*,*) 'WARNING: using default for parameter tag.'
        end if
        if(set_trace.eq.0 .and. (1).eq.0) then
          write(*,*) 'ERROR: parameter trace has not been set.'
          all_ok=0
        else if(set_trace.eq.0 .and. (1).eq.1) then
          write(*,*) 'WARNING: using default for parameter trace.'
        end if
        if(set_xc.eq.0 .and. (1).eq.0) then
          write(*,*) 'ERROR: parameter xc has not been set.'
          all_ok=0
        else if(set_xc.eq.0 .and. (1).eq.1) then
          write(*,*) 'WARNING: using default for parameter xc.'
        end if
        if(set_xmax.eq.0 .and. (1).eq.0) then
          write(*,*) 'ERROR: parameter xmax has not been set.'
          all_ok=0
        else if(set_xmax.eq.0 .and. (1).eq.1) then
          write(*,*) 'WARNING: using default for parameter xmax.'
        end if
        if(set_xmin.eq.0 .and. (1).eq.0) then
          write(*,*) 'ERROR: parameter xmin has not been set.'
          all_ok=0
        else if(set_xmin.eq.0 .and. (1).eq.1) then
          write(*,*) 'WARNING: using default for parameter xmin.'
        end if
        if(set_xwid.eq.0 .and. (1).eq.0) then
          write(*,*) 'ERROR: parameter xwid has not been set.'
          all_ok=0
        else if(set_xwid.eq.0 .and. (1).eq.1) then
          write(*,*) 'WARNING: using default for parameter xwid.'
        end if
        if(set_yc.eq.0 .and. (1).eq.0) then
          write(*,*) 'ERROR: parameter yc has not been set.'
          all_ok=0
        else if(set_yc.eq.0 .and. (1).eq.1) then
          write(*,*) 'WARNING: using default for parameter yc.'
        end if
        if(set_ymax.eq.0 .and. (1).eq.0) then
          write(*,*) 'ERROR: parameter ymax has not been set.'
          all_ok=0
        else if(set_ymax.eq.0 .and. (1).eq.1) then
          write(*,*) 'WARNING: using default for parameter ymax.'
        end if
        if(set_ymin.eq.0 .and. (1).eq.0) then
          write(*,*) 'ERROR: parameter ymin has not been set.'
          all_ok=0
        else if(set_ymin.eq.0 .and. (1).eq.1) then
          write(*,*) 'WARNING: using default for parameter ymin.'
        end if
        if(set_ywid.eq.0 .and. (1).eq.0) then
          write(*,*) 'ERROR: parameter ywid has not been set.'
          all_ok=0
        else if(set_ywid.eq.0 .and. (1).eq.1) then
          write(*,*) 'WARNING: using default for parameter ywid.'
        end if
        if(set_out_gf.eq.0 .and. (1).eq.0) then
          write(*,*) 'ERROR: attribute out_gf has not been set.'
          all_ok=0
        else if(set_out_gf.eq.0 .and. (1).eq.1) then
          write(*,*) 'WARNING: using default for attribute out_gf.'
        end if
        check_params_attribs=all_ok
        return
      end

!----------------------------------------------------------------------
!  read_params_attribs
!----------------------------------------------------------------------
      subroutine read_params_attribs()
        implicit none

        include 'sys_param.inc'
        include 'gfuni0.inc'
        include 'globals.inc'
        include 'other_glbs.inc'
        integer res, gft_read_id_int
        integer gft_read_id_float, gft_read_id_str

        if(set_amp.eq.0) then
          res=gft_read_id_float('amp',amp,1)
        end if
        if(set_camp.eq.0) then
          res=gft_read_id_float('camp',camp,1)
        end if
        if(set_epsiter.eq.0) then
          res=gft_read_id_float('epsiter',epsiter,1)
        end if
        if(set_epsiterid.eq.0) then
          res=gft_read_id_float('epsiterid',epsiterid,1)
        end if
        if(set_fout.eq.0) then
          res=gft_read_id_int('fout',fout,1)
        end if
        if(set_in_file.eq.0) then
          res=gft_read_id_str('in_file',in_file,1)
        end if
        if(set_iter.eq.0) then
          res=gft_read_id_int('iter',iter,1)
        end if
        if(set_lambda.eq.0) then
          res=gft_read_id_float('lambda',lambda,1)
        end if
        if(set_level.eq.0) then
          res=gft_read_id_int('level',level,1)
        end if
        if(set_maxstep.eq.0) then
          res=gft_read_id_int('maxstep',maxstep,1)
        end if
        if(set_maxstepid.eq.0) then
          res=gft_read_id_int('maxstepid',maxstepid,1)
        end if
        if(set_mgeps.eq.0) then
          res=gft_read_id_float('mgeps',mgeps,1)
        end if
        if(set_mgmaxcyc.eq.0) then
          res=gft_read_id_int('mgmaxcyc',mgmaxcyc,1)
        end if
        if(set_mpi.eq.0) then
          res=gft_read_id_float('mpi',mpi,1)
        end if
        if(set_Nx0.eq.0) then
          res=gft_read_id_int('Nx0',Nx0,1)
        end if
        if(set_Ny0.eq.0) then
          res=gft_read_id_int('Ny0',Ny0,1)
        end if
        if(set_out_file.eq.0) then
          res=gft_read_id_str('out_file',out_file,1)
        end if
        if(set_output.eq.0) then
          res=gft_read_id_int('output',output,5)
        end if
        if(set_s_step.eq.0) then
          res=gft_read_id_int('s_step',s_step,1)
        end if
        if(set_ser.eq.0) then
          res=gft_read_id_int('ser',ser,1)
        end if
        if(set_start_t.eq.0) then
          res=gft_read_id_float('start_t',start_t,1)
        end if
        if(set_tag.eq.0) then
          res=gft_read_id_str('tag',tag,1)
        end if
        if(set_trace.eq.0) then
          res=gft_read_id_int('trace',trace,5)
        end if
        if(set_xc.eq.0) then
          res=gft_read_id_float('xc',xc,1)
        end if
        if(set_xmax.eq.0) then
          res=gft_read_id_float('xmax',xmax,1)
        end if
        if(set_xmin.eq.0) then
          res=gft_read_id_float('xmin',xmin,1)
        end if
        if(set_xwid.eq.0) then
          res=gft_read_id_float('xwid',xwid,1)
        end if
        if(set_yc.eq.0) then
          res=gft_read_id_float('yc',yc,1)
        end if
        if(set_ymax.eq.0) then
          res=gft_read_id_float('ymax',ymax,1)
        end if
        if(set_ymin.eq.0) then
          res=gft_read_id_float('ymin',ymin,1)
        end if
        if(set_ywid.eq.0) then
          res=gft_read_id_float('ywid',ywid,1)
        end if
        return
      end

!----------------------------------------------------------------------
!  read_state
!----------------------------------------------------------------------
      subroutine read_state(argc,argv)
        implicit none

        include 'sys_param.inc'
        include 'gfuni0.inc'
        include 'globals.inc'

        integer argc,i
        character*64 argv(argc)
        character*256 command
        integer getu,ret
        integer fp, indlnb, gft_read_id_gf

        fp=getu()
        open(unit=fp,file=in_file(1:indlnb(in_file)),
     &    status = 'old', form = 'unformatted', iostat = ret)
        if(ret .ne. 0) then
          write (*,*) 'Can''t open ',in_file
          write (*,*) 'Calling initial data generator.'
          command='diff2dcn_init'
          do i=1,argc
            command=command(1:indlnb(command))//' '//argv(i)
          end do
          call system(command)
          open(unit=fp,file=in_file(1:indlnb(in_file)),
     &      status = 'old', form = 'unformatted', iostat = ret)
          if(ret .ne. 0) then
            write (*,*) 'Can''t open ',in_file
            write (*,*) 'Assuming updates will do initialization'
            return
          else
            close(fp)
          end if
        else
          close(fp)
        end if
        call gft_set_single(in_file(1:indlnb(in_file)))
        ret=gft_read_id_gf('u[0]',ptrs(getpshape(1)),
     &    getrank(1),q(ptrs(gf_st(1)+2)))
        ret=gft_read_id_gf('rhs[0]',ptrs(getpshape(1)),
     &    getrank(1),q(ptrs(gf_st(2)+1)))
        ret=gft_read_id_gf('ures[0]',ptrs(getpshape(1)),
     &    getrank(1),q(ptrs(gf_st(3)+2)))
        call read_params_attribs()
        call gft_set_multi()
        return
      end

!----------------------------------------------------------------------
!  write_params_attribs
!----------------------------------------------------------------------
      subroutine write_params_attribs()
        implicit none

        include 'sys_param.inc'
        include 'gfuni0.inc'
        include 'globals.inc'
        include 'other_glbs.inc'
        integer res, gft_write_id_int
        integer gft_write_id_float, gft_write_id_str

        res=gft_write_id_float('amp',amp,1)
        res=gft_write_id_float('camp',camp,1)
        res=gft_write_id_float('epsiter',epsiter,1)
        res=gft_write_id_float('epsiterid',epsiterid,1)
        res=gft_write_id_int('fout',fout,1)
        res=gft_write_id_str('in_file',in_file,1)
        res=gft_write_id_int('iter',iter,1)
        res=gft_write_id_float('lambda',lambda,1)
        res=gft_write_id_int('level',level,1)
        res=gft_write_id_int('maxstep',maxstep,1)
        res=gft_write_id_int('maxstepid',maxstepid,1)
        res=gft_write_id_float('mgeps',mgeps,1)
        res=gft_write_id_int('mgmaxcyc',mgmaxcyc,1)
        res=gft_write_id_float('mpi',mpi,1)
        res=gft_write_id_int('Nx0',Nx0,1)
        res=gft_write_id_int('Ny0',Ny0,1)
        res=gft_write_id_str('out_file',out_file,1)
        res=gft_write_id_int('output',output,5)
        res=gft_write_id_int('s_step',s_step,1)
        res=gft_write_id_int('ser',ser,1)
        res=gft_write_id_float('start_t',start_t,1)
        res=gft_write_id_str('tag',tag,1)
        res=gft_write_id_int('trace',trace,5)
        res=gft_write_id_float('xc',xc,1)
        res=gft_write_id_float('xmax',xmax,1)
        res=gft_write_id_float('xmin',xmin,1)
        res=gft_write_id_float('xwid',xwid,1)
        res=gft_write_id_float('yc',yc,1)
        res=gft_write_id_float('ymax',ymax,1)
        res=gft_write_id_float('ymin',ymin,1)
        res=gft_write_id_float('ywid',ywid,1)
        return
      end

!----------------------------------------------------------------------
!  dump_state
!----------------------------------------------------------------------
      subroutine dump_state(t, st, f_name)
        implicit none

        include 'sys_param.inc'
        include 'gfuni0.inc'
        include 'globals.inc'
        real*8  t
        integer st, res
        character*64 f_name
        integer gft_write_id_gf

        call gft_set_single(f_name)
        res=gft_write_id_gf('u[0]',ptrs(getpshape(1)),
     &    getrank(1),q(ptrs(gf_st(1)+2)))
        res=gft_write_id_gf('rhs[0]',ptrs(getpshape(1)),
     &    getrank(1),q(ptrs(gf_st(2)+1)))
        res=gft_write_id_gf('ures[0]',ptrs(getpshape(1)),
     &    getrank(1),q(ptrs(gf_st(3)+2)))
        start_t=t
        s_step=st
        call write_params_attribs()
        call gft_set_multi()
        return
      end
!----------------------------------------------------------------------
!  handler
!----------------------------------------------------------------------
      subroutine handler
        implicit none

        include 'sys_param.inc'
        include 'gfuni0.inc'
        include 'globals.inc'

        character*(5) temp
        integer quit,ch,i,encode
        integer rs
        quit=0
99      format (a)
98      format (i3,a)
100     continue
          write(*,*) '1.  Change output frequency'
          write(*,*) '2.  View out_gf'
          write(*,*) '3.  Change out_gf'
          write(*,*) '4.  Resume'
          write(*,*) '5.  Exit'
          write(*,99) 'Enter choice: '
          read(*,*,err=100) ch
          if(ch .eq. 1) then
            write(*,99) 'Enter new frequency 1/'
            read(*,*) rmod
          else if(ch .eq. 2) then
            do i=1,3
              write(*,98) i,'. '//attrfname(i)
              if(out_gf(i) .eq. 1) then
                write(*,*) ' ON'
              else
                write(*,*) ' OFF'
              end if
            end do
          else if(ch .eq. 3) then
            write(*,99) 'Enter grid function number: '
            read(*,*) i
            if(i .gt. 0 .and. i .le. 3) then
              out_gf(i)=mod(out_gf(i)+1,2)
              if(out_gf(i) .eq. 0) then
                call gft_close(attrfname(i))
              end if
            end if
          else if(ch .eq. 4) then
            quit=1
          else if(ch .eq. 5) then
            quit=1
            rnpldone=1
          end if
        if(quit .ne. 1) then
          goto 100
        end if
        nout_gf=encode(lout_gf,out_gf,3)
        return
      end

!----------------------------------------------------------------------
!  output_gfuncs
!----------------------------------------------------------------------
      subroutine output_gfuncs(o_fout, o_ser, t, o_out_gf, fname,u_n,rhs
     &  ,ures_n,g1_rank,g1_shape,g1_cnames,g1_crds)
        implicit none

        include 'sys_param.inc'
        include 'gfuni0.inc'
        include 'globals.inc'
        integer o_fout, o_ser
        real*8  t
        integer o_out_gf(3)
        character*64 fname(3)
        real*8 u_n
        real*8 rhs
        real*8 ures_n
        integer g1_rank
        integer g1_shape(*)
        character*64 g1_cnames
        real*8 g1_crds
        integer res
        integer gft_out_full,rvsxynt

        if(o_out_gf(1).eq.1) then
          if(o_fout.eq.1) then
            res=gft_out_full(fname(1),t,g1_shape,g1_cnames,g1_rank,
     &        g1_crds,u_n)
          end if
        end if
        if(o_out_gf(2).eq.1) then
          if(o_fout.eq.1) then
            res=gft_out_full(fname(2),t,g1_shape,g1_cnames,g1_rank,
     &        g1_crds,rhs)
          end if
        end if
        if(o_out_gf(3).eq.1) then
          if(o_fout.eq.1) then
            res=gft_out_full(fname(3),t,g1_shape,g1_cnames,g1_rank,
     &        g1_crds,ures_n)
          end if
        end if
        return
      end
!----------------------------------------------------------------------
!  cleanup
!----------------------------------------------------------------------
      subroutine cleanup
        implicit none

        include 'sys_param.inc'
        include 'gfuni0.inc'
        include 'globals.inc'

        call gft_close_all()
        return
      end

!-----------------------------------------------------------------------
!
!   This is the f77 main
!    
!-----------------------------------------------------------------------

      program     main

      implicit       none

      include       'sys_param.inc'
      include       'gfuni0.inc'
      integer argc
      character*64 argv(64)

      include       'globals.inc'
      include       'other_glbs.inc'

      integer      iargc, getu, indlnb,do_out
      integer      got_pf,argerr,do_ivec
      integer      IVEL
      parameter  ( IVEL = 4 )

      character*64 param_file,optarg
      integer      fp,i,j,steps,rc,ot,max,ret
      real*8       t
      integer      update_gfuncs, query_interrupt, check_params_attribs
      
      argc=iargc()
      do i=1,argc
        call getarg(i,argv(i))
      end do
      call mmini(memsiz)
      call init_params_attribs()
      i=1
      argerr=0
      got_pf=0
100   continue
      if(i.le.argc) then
        optarg=argv(i)
        if(optarg .eq. '-p') then
          i=i+1
          optarg=argv(i)
          call sread_parameters(optarg)
        else if(optarg(1:1).eq.'-') then
          argerr=1
        else
          param_file=optarg
          got_pf=1
        end if
        i=i+1
        goto 100
      end if
      if(argerr .eq. 1) then
        write(*,*) 'Usage: diff2dcn'
        write(*,*) '     [ -p "parameter:=value" ]'
        write(*,*) '     [ parameter_file ]'
        stop
      end if
      if(got_pf .ne. 1) then
101     continue 
        write(*,*) 'Enter name of parameter file: '
        read(*,1010,end=101,err=101) param_file
1010    format(a)  
      end if
      fp=getu()
      open(unit=fp,file=param_file(1:indlnb(param_file)),
     &     status = 'old', form = 'unformatted', iostat = rc)
      if(rc .ne. 0) then
        write(*,*) 'Unable to open initial data file '//
     &             param_file(1:indlnb(param_file))
        stop 
      end if
      close(fp)
      
      call read_parameters(param_file(1:indlnb(param_file)))
      call read_attributes(param_file(1:indlnb(param_file)))
      call read_attributes('.rnpl.attributes')
      call init_coord_difs()
      call init_lats()
      call init_gfuncs()
      call read_state(argc,argv)
      ret=check_params_attribs()
      if(ret.eq.0) then
        write(*,*) 'diff2dcn: Unable to continue.'
        stop
      end if

      iter = iter * 2.0**level

!------------------------------------------------------------------------
!   set up signal handler 
!------------------------------------------------------------------------

      call set_interrupt()
 
!------------------------------------------------------------------------
!   start main loop
!------------------------------------------------------------------------
      rmod = -1
      rnpldone = 0
      t = start_t
      i = s_step
      if(set_epsiterid.eq.0) then
        epsiterid=epsiter
      end if
      if(set_trace.eq.0) then
        call rivcpy(trace,output,output(1)*IVEL+1)
      end if
!   fix ivec indicies
      call fixup_ivec(0,iter,level,output)
      call fixup_ivec(0,iter,level,trace)

      if(i.eq.0) then
        if(do_ivec(i,iter,trace).eq.1) then
          write(*,*) 'Starting evolution.  step: ',i,' at t=',t
        end if
        if(((rmod.ne.-1).and.(mod(i,rmod).eq.0)).or.
     &     (do_ivec(i,iter,output).eq.1)) then
          call output_gfuncs(fout,ser,t,out_gf,attrfname,
     &      q(ptrs(gf_st(1)+2)),q(ptrs(gf_st(2)+1)),q(ptrs(gf_st(3)+2)),
     &      getrank(1),ptrs(getpshape(1)),cname(getlatcs(1)),
     &      q(getpcoord(1,1)))
        end if
      end if
      j=i+1
      do i=j,iter
        steps=update_gfuncs(t)
        t = t+dt
        if(do_ivec(i,iter,trace).eq.1) then
          write(*,*) 'step: ',i,' t=',t, ' steps=',steps
        end if
        if(((rmod.ne.-1).and.(mod(i,rmod).eq.0)).or.
     &     (do_ivec(i,iter,output).eq.1)) then
          call output_gfuncs(fout,ser,t,out_gf,attrfname,
     &      q(ptrs(gf_st(1)+2)),q(ptrs(gf_st(2)+1)),q(ptrs(gf_st(3)+2)),
     &      getrank(1),ptrs(getpshape(1)),cname(getlatcs(1)),
     &      q(getpcoord(1,1)))
        end if
        if(query_interrupt().ne.0) then
          call handler()
          call reset_interrupt()
        end if
        if(rnpldone .eq. 1) then
          goto 9999
        end if
      end do
9999  continue
      call dump_state(t,i,out_file)
      call cleanup()
      stop

      end