*****************************************************************
*                                                               *
*****************************************************************      
      implicit double complex(c)
      IMPLICIT double precision (A-B,D-H,O-Z)
      parameter(lx = 600)
      parameter(ly = 600)
      parameter(ione = 1)
      parameter(itwo = 2)
      parameter(ithr = 3)
      parameter(zero = 0.0d0)
      parameter(rzero = 0.0d0)
      parameter(rnzero = 0.02d0)
      parameter(one = 1.0d0)
      parameter(rone = 1.0d0)
      parameter(two = 2.0d0)
      parameter(three = 3.0d0)
      parameter(four = 4.0d0)
      parameter(pi = 3.14159265d0)
      parameter(twopi = 6.2831853d0)
      parameter(hlf = 0.5d0)
      parameter(qrt = 0.25d0)
      parameter(itwo8 = 256)
      parameter(two8 = 256.0d0)
      dimension cpsi(lx+4,ly+4),cpsie(lx+4,ly+4), 
     &     ceaxe1(lx+4,ly+4),ceaxe2(lx+4,ly+4),
     &     ceaxo1(lx+4,ly+4),ceaxo2(lx+4,ly+4),
     &     ceaye1(lx+4,ly+4),ceaye2(lx+4,ly+4),
     &     ceayo1(lx+4,ly+4),ceayo2(lx+4,ly+4)
      czero = dcmplx(0.0d0,0.0d0)
      cone = dcmplx(1.0d0,0.0d0)
      ci = dcmplx(0.0d0,1.0d0)
      rt6 = sqrt(6.0d0)
      rt3 = sqrt(3.0d0)
      rt2 = sqrt(2.0d0)
      lxh = int(hlf*dble(lx))
      lyh = int(hlf*dble(ly))
      px = pi
      xmn = -12.0d0
      xmx = 12.0d0
      ymn = -12.0d0
      ymx = 12.0d0
      nfedti = 700
      nfedtd = 701
      nfedtp = 702
      nfeti = 399
      nfetd = 499
      nfetp = 599
      nfxv1 = 999
      nfxv2 = 1999
      nfgp1 = 2999
      open(unit = nfedti,status = "new")
      open(unit = nfedtd,status = "new")
      open(unit = nfedtp,status = "new")
      eps = (10.0d0)**(0.2d0)
      reps = one/(10.0d0)**(0.2d0)
       write(nfedti,"(a,a,x,2i6,1e12.5)") 
     &         '#', 'i', lx, ly, reps
       write(nfedtd,"(a,a,x,2i6,1e12.5)") 
     &         '#', 'd', lx, ly, reps
       write(nfedtp,"(a,a,x,2i6,1e12.5)") 
     &         '#', 'p', lx, ly, reps
c--------------------------------------------------------
      do 1986 ke = 1,8
      eps = reps*eps
      omega = 0.5d0
      px = pi
      py = zero
      x0 = zero
      y0 = px
c*-*-*-*-*-* lx-1 ---> lx, ly-1 ---> ly  *-*-*-*-*-*-*
      dx = (xmx-xmn)/(lx)
      dy = (ymx-ymn)/(ly)
      rdx = one/dx
      rdy = one/dy
      rdx2 = rdx*rdx
      rdy2 = rdy*rdy
c*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
      dt = dx*dx*eps
      bx = -hlf*rdx2
      by = -hlf*rdy2
      vx = hlf*rdx2
      vy = hlf*rdy2
      rcx = qrt*rdx
      rcy = qrt*rdy
c*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
      al = omega*dt
c*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
      n = int(two*pi/omega/dt)
      k1 = int(0.125d0*two*pi/omega/dt)
      ki = int(0.125d0*two*pi/omega/dt)
      n = n+1
c       write(6,"(2i6,2e12.5,2i9)" )
c     &         lx, ly, eps, dt,  n, k1
****************************************************************
        r1 = zero
        y = ymn
        do 23000 j = 3,ly+2 
          i1 = 1
          x = xmn
          do 23002 i = 3,lx+2
            rtmp1 = -hlf*((y-y0)*(y-y0)+(x-x0)*(x-x0))
            rtmp2 = px*x
            if (rtmp1.lt.-100.0d0) then
              rtmp1 = -100.0d0
            endif
            ctmp1 = dcmplx(rtmp1,rtmp2)
            cpsi(i,j) = (exp(ctmp1))
            tmp1 = dble(cpsi(i,j))
            tmp2 = dimag(cpsi(i,j))
            r1 = r1+(tmp1*tmp1+tmp2*tmp2)*dx*dy
            x = x+dx
23002     continue
          y = y+dy
23000   continue
        r1 = one/sqrt(r1)
        y = ymn
        do 23021 j = 3,ly+2 
          i1 = 1
          x = xmn
          do 23222 i = 3,lx+2
            cpsi(i,j) = cpsi(i,j)*r1
            x = x+dx
23222     continue
          y = y+dy
23021  continue
        do 450 j=3,ly+2
          cpsi(1,j)=cpsi(lx+1,j)
          cpsi(2,j)=cpsi(lx+2,j)
          cpsi(lx+3,j)=cpsi(3,j)
          cpsi(lx+4,j)=cpsi(4,j)
  450   continue
        do 550 i=3,lx+2
          cpsi(i,1)=cpsi(i,ly+1)
          cpsi(i,2)=cpsi(i,ly+2)
          cpsi(i,ly+3)=cpsi(i,3)
          cpsi(i,ly+4)=cpsi(i,4)
  550   continue
***************************************************************
c--------------< eaxe. >----------------------------------
       y = ymn
       do 23110 j = 3,ly+2
         x = xmn
         do 23108 i1 = 1,lxh
           i = itwo*i1+1
           ax0  = vpx(x,y)
           ay0  = vpy(x,y)
           ax1  = vpx(x+dx,y)
           ay1  = vpy(x,y+dy)
           v0 = vx+qrt*(vs(x,y)+hlf*(ax0*ax0+ay0*ay0))
           v1 = vx+qrt*(vs(x,y)+hlf*(ax1*ax1+ay0*ay0))
c
           a = hlf*(v0+v1)
           d = hlf*(v0-v1)
           rc01 = rcx*(ax0+ax1)
c
           b2c2 = bx*bx+rc01*rc01
           cbpic = dcmplx(bx,rc01)
           cbmic = dcmplx(bx,-rc01)
           rl = sqrt(b2c2+d*d)
           rlma = (rl-a)*dt
           rlpa = (rl+a)*dt
           rlmd =  rl-d
           rlmd2 = rlmd*rlmd
           rtworl = one/(two*rl)
           rtworld = rtworl/rlmd
c
           rcop = cos(rlpa)
           rsip = sin(rlpa)
           rcom = cos(rlma)
           rsim = sin(rlma)
           cep = dcmplx(rcom, rsim)
           cem = dcmplx(rcop,-rsip)
c
           ceaxe1(i,j) =   rtworld*(rlmd2*cep+b2c2*cem)
           ceaxe1(i+1,j) = -rtworl*(cbpic*cep-cbpic*cem)
           ceaxe2(i,j) =   -rtworl*(cbmic*cep-cbmic*cem)
           ceaxe2(i+1,j) = rtworld*(b2c2*cep+rlmd2*cem)
           x = x+two*dx
23108    continue
           ceaxe1(1,j) = ceaxe1(lx+1,j)
           ceaxe2(1,j) = ceaxe2(lx+1,j)
           ceaxe1(2,j) = ceaxe1(lx+2,j)
           ceaxe2(2,j) = ceaxe2(lx+2,j)
           ceaxe1(lx+3,j) = ceaxe1(3,j)
           ceaxe2(lx+3,j) = ceaxe2(3,j)
           ceaxe1(lx+4,j) = ceaxe1(4,j)
           ceaxe2(lx+4,j) = ceaxe2(4,j)
         y = y+dy
23110  continue
       do 23220 i = 3,lx+2
           ceaxe1(i,1) = ceaxe1(i,ly+1)
           ceaxe2(i,1) = ceaxe2(i,ly+1)
           ceaxe1(i,2) = ceaxe1(i,ly+2)
           ceaxe2(i,2) = ceaxe2(i,ly+2)
           ceaxe1(i,ly+3) = ceaxe1(i,3)
           ceaxe2(i,ly+3) = ceaxe2(i,3)
           ceaxe1(i,ly+4) = ceaxe1(i,4)
           ceaxe2(i,ly+4) = ceaxe2(i,4)
23220  continue
c--------------< eaxo. >----------------------------------
       y = ymn
       do 24110 j = 3,ly+2
         x = xmn-dx
         do 24108 i1 = 1,lxh+1
           i = itwo*i1
           ax0  = vpx(x,y)
           ay0  = vpy(x,y)
           ax1  = vpx(x+dx,y)
           ay1  = vpy(x,y+dy)
           v0 = vx+qrt*(vs(x,y)+hlf*(ax0*ax0+ay0*ay0))
           v1 = vx+qrt*(vs(x,y)+hlf*(ax1*ax1+ay0*ay0))
c
           a = hlf*(v0+v1)
           d = hlf*(v0-v1)
           rc01 = rcx*(ax0+ax1)
c
           b2c2 = bx*bx+rc01*rc01
           cbpic = dcmplx(bx,rc01)
           cbmic = dcmplx(bx,-rc01)
           rl = sqrt(b2c2+d*d)
           rlma = (rl-a)*dt
           rlpa = (rl+a)*dt
           rlmd =  rl-d
           rlmd2 = rlmd*rlmd
           rtworl = one/(two*rl)
           rtworld = rtworl/rlmd
c
           rcop = cos(rlpa)
           rsip = sin(rlpa)
           rcom = cos(rlma)
           rsim = sin(rlma)
           cep = dcmplx(rcom, rsim)
           cem = dcmplx(rcop,-rsip)
           ceaxo1(i,j) =   rtworld*(rlmd2*cep+b2c2*cem)
           ceaxo1(i+1,j) = -rtworl*(cbpic*cep-cbpic*cem)
           ceaxo2(i,j) =   -rtworl*(cbmic*cep-cbmic*cem)
           ceaxo2(i+1,j) = rtworld*(b2c2*cep+rlmd2*cem)
           x = x+two*dx
24108    continue
           ceaxo1(1,j) = ceaxo1(lx+1,j)
           ceaxo2(1,j) = ceaxo2(lx+1,j)
           ceaxo1(2,j) = ceaxo1(lx+2,j)
           ceaxo2(2,j) = ceaxo2(lx+2,j)
           ceaxo1(lx+3,j) = ceaxo1(3,j)
           ceaxo2(lx+3,j) = ceaxo2(3,j)
           ceaxo1(lx+4,j) = ceaxo1(4,j)
           ceaxo2(lx+4,j) = ceaxo2(4,j)
         y = y+dy
24110  continue
       do 24220 i = 3,lx+2
           ceaxo1(i,1) = ceaxo1(i,ly+1)
           ceaxo2(i,1) = ceaxo2(i,ly+1)
           ceaxo1(i,2) = ceaxo1(i,ly+2)
           ceaxo2(i,2) = ceaxo2(i,ly+2)
           ceaxo1(i,ly+3) = ceaxo1(i,3)
           ceaxo2(i,ly+3) = ceaxo2(i,3)
           ceaxo1(i,ly+4) = ceaxo1(i,4)
           ceaxo2(i,ly+4) = ceaxo2(i,4)
24220  continue
c--------------< eaye. >----------------------------------
       x = xmn
       do 25110 i = 3,lx+2
         y = ymn
         do 25108 j1 = 1,lyh+1
           j = itwo*j1+1
           ax0  = vpx(x,y)
           ay0  = vpy(x,y)
           ax1  = vpx(x+dx,y)
           ay1  = vpy(x,y+dy)
           v0 = vy+qrt*(vs(x,y)+hlf*(ax0*ax0+ay0*ay0))
           v1 = vy+qrt*(vs(x,y)+hlf*(ax0*ax0+ay1*ay1))
c
           a = hlf*(v0+v1)
           d = hlf*(v0-v1)
           rc01 = rcy*(ay0+ay1)
c
           b2c2 = by*by+rc01*rc01
           cbpic = dcmplx(by,rc01)
           cbmic = dcmplx(by,-rc01)
           rl = sqrt(b2c2+d*d)
           rlma = (rl-a)*dt
           rlpa = (rl+a)*dt
           rlmd =  rl-d
           rlmd2 = rlmd*rlmd
           rtworl = one/(two*rl)
           rtworld = rtworl/rlmd
c
           rcop = cos(rlpa)
           rsip = sin(rlpa)
           rcom = cos(rlma)
           rsim = sin(rlma)
           cep = dcmplx(rcom, rsim)
           cem = dcmplx(rcop,-rsip)
           ceaye1(i,j) =   rtworld*(rlmd2*cep+b2c2*cem)
           ceaye1(i,j+1) = -rtworl*(cbpic*cep-cbpic*cem)
           ceaye2(i,j) =   -rtworl*(cbmic*cep-cbmic*cem)
           ceaye2(i,j+1) = rtworld*(b2c2*cep+rlmd2*cem)
           y = y+two*dy
25108    continue
           ceaye1(i,1) = ceaye1(i,ly+1)
           ceaye2(i,1) = ceaye2(i,ly+1)
           ceaye1(i,2) = ceaye1(i,ly+2)
           ceaye2(i,2) = ceaye2(i,ly+2)
           ceaye1(i,ly+3) = ceaye1(i,3)
           ceaye2(i,ly+3) = ceaye2(i,3)
           ceaye1(i,ly+4) = ceaye1(i,4)
           ceaye2(i,ly+4) = ceaye2(i,4)
         x = x+dx
25110  continue
       do 25220 j = 3,ly+2
           ceaye1(1,j) = ceaye1(lx+1,j)
           ceaye2(1,j) = ceaye2(lx+1,j)
           ceaye1(2,j) = ceaye1(lx+2,j)
           ceaye2(2,j) = ceaye2(lx+2,j)
           ceaye1(lx+3,j) = ceaye1(3,j)
           ceaye2(lx+3,j) = ceaye2(3,j)
           ceaye1(lx+4,j) = ceaye1(4,j)
           ceaye2(lx+4,j) = ceaye2(4,j)
25220  continue
c--------------< eayo. >----------------------------------
       x = xmn
       do 26110 i = 3,lx+2
         y = ymn-dy
         do 26108 j1 = 1,lyh+1
           j = itwo*j1
           ax0  = vpx(x,y)
           ay0  = vpy(x,y)
           ax1  = vpx(x+dx,y)
           ay1  = vpy(x,y+dy)
           v0 = vy+qrt*(vs(x,y)+hlf*(ax0*ax0+ay0*ay0))
           v1 = vy+qrt*(vs(x,y)+hlf*(ax0*ax0+ay1*ay1))
c
           a = hlf*(v0+v1)
           d = hlf*(v0-v1)
           rc01 = rcy*(ay0+ay1)
c
           b2c2 = by*by+rc01*rc01
           cbpic = dcmplx(by,rc01)
           cbmic = dcmplx(by,-rc01)
           rl = sqrt(b2c2+d*d)
           rlma = (rl-a)*dt
           rlpa = (rl+a)*dt
           rlmd =  rl-d
           rlmd2 = rlmd*rlmd
           rtworl = one/(two*rl)
           rtworld = rtworl/rlmd
c
           rcop = cos(rlpa)
           rsip = sin(rlpa)
           rcom = cos(rlma)
           rsim = sin(rlma)
           cep = dcmplx(rcom, rsim)
           cem = dcmplx(rcop,-rsip)
           ceayo1(i,j) =   rtworld*(rlmd2*cep+b2c2*cem)
           ceayo1(i,j+1) = -rtworl*(cbpic*cep-cbpic*cem)
           ceayo2(i,j) =   -rtworl*(cbmic*cep-cbmic*cem)
           ceayo2(i,j+1) = rtworld*(b2c2*cep+rlmd2*cem)
           y = y+two*dy
26108    continue
           ceayo1(i,1) = ceayo1(i,ly+1)
           ceayo2(i,1) = ceayo2(i,ly+1)
           ceayo1(i,2) = ceayo1(i,ly+2)
           ceayo2(i,2) = ceayo2(i,ly+2)
           ceayo1(i,ly+3) = ceayo1(i,3)
           ceayo2(i,ly+3) = ceayo2(i,3)
           ceayo1(i,ly+4) = ceayo1(i,4)
           ceayo2(i,ly+4) = ceayo2(i,4)
         x = x+dx
26110  continue
       do 26220 j = 3,ly+2
           ceayo1(1,j) = ceayo1(lx+1,j)
           ceayo2(1,j) = ceayo2(lx+1,j)
           ceayo1(2,j) = ceayo1(lx+2,j)
           ceayo2(2,j) = ceayo2(lx+2,j)
           ceayo1(lx+3,j) = ceayo1(3,j)
           ceayo2(lx+3,j) = ceayo2(3,j)
           ceayo1(lx+4,j) = ceayo1(4,j)
           ceayo2(lx+4,j) = ceayo2(4,j)
26220  continue
***************************************************************
      nfeti = nfeti+1
      nfetd = nfetd+1
      nfetp = nfetp+1
      open(unit = nfeti,status = "new")
      open(unit = nfetd,status = "new")
      open(unit = nfetp,status = "new")
      t = zero
      write(nfeti,"(a,a,x,2i6,3e12.5,3i9)") 
     &             '#', 'i', lx, ly, eps, al, dt, n, k1,ki
      write(nfetd,"(a,a,x,2i6,3e12.5,3i9)") 
     &             '#', 'd', lx, ly, eps, al, dt, n, k1,ki
      write(nfetp,"(a,a,x,2i6,3e12.5,3i9)") 
     &             '#', 'p', lx, ly, eps, al, dt, n, k1,ki
c      write(6,*)   'dt = ', dt, ', n = ', n, ', k1 = ', k1
      do 53016 ka = 1,n 
        if(mod(ka,k1).eq.1) then
***************************************************************
        if(mod(ka,ki).eq.1) then
          nfxv1 = nfxv1+1
          nfxv2 = nfxv2+1
          nfgp1 = nfgp1+1
          open(unit = nfxv1,status = "new")
          write(nfxv1,"(a)") "P3"
          write(nfxv1,"(2i4)") lx, ly
          write(nfxv1,"(a)") "255"
          open(unit = nfxv2,status = "new")
          write(nfxv2,"(a)") "P3"
          write(nfxv2,"(2i4)") lx, ly
          write(nfxv2,"(a)") "255"
          open(unit = nfgp1,status = "new")
           write(nfgp1,"(a,x,2i6,4e12.5,1i6)") 
     &             '#', lx, ly, eps, al, dt, t, ka-1
        endif
**************************************************************
        cerrp = czero
        errd = zero
        s = zero
        sn = zero
        y = ymn
        do 23020 j = 3,ly+2
          x = xmn
          do 23022 i = 3,lx+2
            psix = dble(cpsi(i,j))
            psiy = dimag(cpsi(i,j))
c           --------------------------
            rtmp1 = -hlf*((y-y0)*(y-y0)+(x-x0)*(x-x0))
            rtmp2 = px*x
            if (rtmp1.lt.-100.0d0) then
              rtmp1 = -100.0d0
            endif
            ctmp1 = dcmplx(rtmp1,rtmp2)
            cpsic = r1*(exp(ctmp1))
            tmp1 = dble(cpsic)
            tmp2 = dimag(cpsic)
            psicx = dble(cpsic)
            psicy = dimag(cpsic)
c           -------------------------
            ra2 = psix*psix+psiy*psiy
            ra = sqrt(ra2)
            rac2 = psicx*psicx+psicy*psicy
            rac = sqrt(rac2)
            rae2 = (psix-psicx)*(psix-psicx)
     &            +(psiy-psicy)*(psiy-psicy)
c           ---------------------------
            s = s+rac2*dx*dy
            sn = sn+ra2*dx*dy
            cip = conjg(cpsi(i,j))*cpsic*dx*dy
            cerrp = cerrp+cip
            errd  =  errd+rae2*dx*dy
c           -------------------------
            if(mod(ka,ki).eq.1) then
              rtmp1 = hlf+sign(rone,ra2-rone)
     &              *(hlf-ra/(rone+ra2))
              rtmp2 = psix/(rt6*(rone+ra2))
              rtmp3 = psiy/(rt2*(rone+ra2))
              ir = int((rtmp1+two*rtmp2)*two8)
              ig = int((rtmp1-rtmp2+rtmp3)*two8)
              ib = int((rtmp1-rtmp2-rtmp3)*two8)
              ip = int(200.0d0*ra2)
              if (ip.gt.255) then
                ip = 255
              endif
              if ((mod(j,1).eq.0).and.(mod(i,1).eq.0)) then
                write(nfxv1,"(3i4)") ir,ig,ib
              endif
              rtmp1 = hlf+sign(rone,rac2-rone)
     &              *(hlf-rac/(rone+rac2))
              rtmp2 = psicx/(rt6*(rone+rac2))
              rtmp3 = psicy/(rt2*(rone+rac2))
              ir = int((rtmp1+two*rtmp2)*two8)
              ig = int((rtmp1-rtmp2+rtmp3)*two8)
              ib = int((rtmp1-rtmp2-rtmp3)*two8)
              ip = int(200.0d0*rac2)
              if (ip.gt.255) then
                ip = 255
              endif
              if ((mod(j,1).eq.0).and.(mod(i,1).eq.0)) then
                write(nfxv2,"(3i4)") ir,ig,ib
              endif
              if ((mod(j,1).eq.0).and.(mod(i,1).eq.0)) then
                write(nfgp1,"(f14.5)") ra2
              endif
            endif
            x = x+dx
23022     continue
          if(mod(ka,ki).eq.1) then
            write(nfgp1,*)
          endif
          y = y+dy
23020   continue
        if(mod(ka,ki).eq.1) then
          close(unit = nfxv1 , status = "keep")
          close(unit = nfxv2 , status = "keep")
          close(unit = nfgp1 , status = "keep")
          call system('gzip fort.????')
        endif
        tmp1 = dble(cerrp)
        tmp2 = dimag(cerrp)
        write(nfeti,"(2e16.7,1i8)") 
     &       t*omega, (s-sqrt(tmp1*tmp1+tmp2*tmp2))/s, ka-1
        write(nfetd,"(2e16.7,1i8)") 
     &       t*omega, sqrt(errd), ka-1
        write(nfetp,"(2e16.7,1i8)") 
     &       t*omega, dabs(atan(tmp2/tmp1)), ka-1
c        write(6,"(4e13.5)")     t, sn, s 
      endif
      if (ka.eq.n) then
        write(nfedti,"(3e16.7,1i8)") 
     &       omega*dt,(s-sqrt(tmp1*tmp1+tmp2*tmp2))/s, t, n
        write(nfedtd,"(3e16.7,1i8)") 
     &       omega*dt,sqrt(errd), t, n
        write(nfedtp,"(3e16.7,1i8)") 
     &       omega*dt,dabs(atan(tmp2/tmp1)), t, n
      endif
************************************************************
c-----------< eaxe. >--------------------------------------
      do 39020 j = 3,ly+2
       i = -1
       do 39030 icount = 1,lxh+2
        i = i+itwo
        ip = i+1
        cpsie(i,j) = ceaxe1(i,j)*cpsi(i,j)
     &              +ceaxe1(ip,j)*cpsi(ip,j)
        cpsie(ip,j) = ceaxe2(i,j)*cpsi(i,j)
     &               +ceaxe2(ip,j)*cpsi(ip,j)
39030  continue
39020 continue
c-----------< eaxo. >--------------------------------------
      do 39620 j = 3,ly+2
       i = 0
       do 39630 icount = 1,lxh+1
        i = i+itwo
        ip = i+1
        cpsi(i,j) = ceaxo1(i,j)*cpsie(i,j)
     &             +ceaxo1(ip,j)*cpsie(ip,j)
        cpsi(ip,j) = ceaxo2(i,j)*cpsie(i,j)
     &              +ceaxo2(ip,j)*cpsie(ip,j)
39630  continue
39620 continue
      do 40010 j = 3,ly+2
        cpsi(1,j) = cpsi(lx+1,j)
        cpsi(2,j) = cpsi(lx+2,j)
        cpsi(lx+3,j) = cpsi(3,j)
        cpsi(lx+4,j) = cpsi(4,j)
40010 continue 
      do 40020 i = 3,lx+2
        cpsi(i,1) = cpsi(i,ly+1)
        cpsi(i,2) = cpsi(i,ly+2)
        cpsi(i,ly+3) = cpsi(i,3)
        cpsi(i,ly+4) = cpsi(i,4)
40020 continue 
c-----------< eaye. >--------------------------------------
      do 99020 i = 3,lx+2
       j = -1
       do 99030 jcount = 1,lyh+2
        j = j+2
        jp = j+1
        cpsie(i,j) = ceaye1(i,j)*cpsi(i,j)
     &              +ceaye1(i,jp)*cpsi(i,jp)
        cpsie(i,jp) = ceaye2(i,j)*cpsi(i,j)
     &               +ceaye2(i,jp)*cpsi(i,jp)
99030  continue
99020 continue
c-----------< eayo. >--------------------------------------
      do 99620 i = 3,lx+2
       j = 0
       do 99630 jcount = 1,lyh+1
        j = j+2
        jp = j+1
        cpsi(i,j) = ceayo1(i,j)*cpsie(i,j)
     &             +ceayo1(i,jp)*cpsie(i,jp)
        cpsi(i,jp) = ceayo2(i,j)*cpsie(i,j)
     &              +ceayo2(i,jp)*cpsie(i,jp)
99630  continue
99620 continue
      do 90010 j = 3,ly+2
       cpsi(1,j) = cpsi(lx+1,j)
       cpsi(2,j) = cpsi(lx+2,j)
       cpsi(lx+3,j) = cpsi(3,j)
       cpsi(lx+4,j) = cpsi(4,j)
90010 continue 
      do 90020 i = 3,lx+2
        cpsi(i,1) = cpsi(i,ly+1)
        cpsi(i,2) = cpsi(i,ly+2)
        cpsi(i,ly+3) = cpsi(i,3)
        cpsi(i,ly+4) = cpsi(i,4)
90020 continue 
      t = t+dt
53016 continue
      close(unit = nfeti , status = "keep")
      close(unit = nfetd , status = "keep")
      close(unit = nfetp , status = "keep")
 1986 continue      
      close(unit = nfedti,status = "keep")
      close(unit = nfedtd,status = "keep")
      close(unit = nfedtp,status = "keep")
***********************************************************
      stop
      end
***********************************************************
      double precision function vpx(x,y)
      implicit double precision (a-h,o-z)
      pi = 3.14159265359d0
      hlf = 0.5d0
      b = 1.0d0
      vpx = -hlf*b*y
      return
      end
***********************************************************
      double precision function vpy(x,y)
      implicit double precision (a-h,o-z)
      pi = 3.14159265359d0
      hlf = 0.5d0
      b = 1.0d0
      vpy = hlf*b*x
      return
      end
***********************************************************
      double precision function vs(x,y)
      implicit double precision (a-h,o-z)
      zero = 0.0d0
      pi = 3.1415926535d0
      pr2 = 0.25d0
      vtmp = 10.0d0
      if (x*x+(y+pi)*(y+pi).le.pr2) then
        vs = zero
      else
        vs = zero
      endif
      return
      end