      implicit real * 8 (a-h,o-z)
      parameter(ldat=12345, npa=12345, mpa=npa*5,mdl=23333)
      common /data/ xx(ldat),yy(ldat),zz(ldat),zmg(ldat),nn
      common /triang/ adt(mdl),ldt(mdl,3),ndt(mdl,3),np,npext,ndelo
      common /coef/ aa(npa),bb(npa),cc(npa)
      common /approx/ xmg0
      dimension b(mpa),bv(npa),bv1(npa)
      dimension xi(50000),xj(50000),xoi(50000),bvxy(50000)
      dimension xl1(50000),xl2(50000)
      dimension xl(50000,60)
      dimension xd(3),yd(3),zd(3),zb(3)
      dimension b1(npa),b2(npa),b4(npa),b5(npa),b7(npa)
      dimension xx1(npa),yy1(npa),zz1(npa),zmg1(npa)
      dimension bb2(npa),bb4(npa),bb5(npa),bb7(npa)
      dimension bbv(npa),a11(npa)
      dimension aa1(npa),bb1(npa),cc1(npa)
      dimension xxb(npa),yyb(npa),ldtb(mdl,3)
      det(x1,x2,x3,y1,y2,y3)=x1*y2+x2*y3+x3*y1-x1*y3-x2*y1-x3*y2
      eps=1.d-10
      call redata(xx,yy,zz,zmg,tx,ty,tz,nn,xmin,ymind)
      txx=tx/cos((ymind+ty/2)*3.14159/180.)
      write(6,*) txx
c     open(8,file='work.configb')
c     read(8,*) np,npex,nd,tx,ty
c     read(8,*) (idum,xxb(ipt),yyb(ipt),d1,d2,d3,d4,d5,ipt=1,npex)
c     do i=1,npex
c     xxb(i)=xxb(i)*cos((20.d0+ty/2)*3.14159/180.)
c     end do
c     tx=tx*cos((20.d0+ty/2)*3.14159/180.)
c     do i=1,nd
c     read(8,*) idum,(ldtb(i,j),j=1,3),adum
c     enddo
c     close(8)
c     open(2,file='alljpM40nBvalues.out')
c     read(2,*)  w1,abic,mx1
c     read(2,*) (bv(i), i=1,mx1)
c     close(2)
c     open(2,file='work.bval')
c     read(2,*) nxy
c     read(2,7) (xoi(i),xj(i),bvxy(i),i=1,nxy)
      open(8,file='hist-etas-mk-forecast.conf')
      read(8,*) tx,ty,tz
      read(8,*) xxmin,yymin
      yminnode=yymin
       close(8)
      open(2,file='node.dat')
      i=1
  555 read(2,*,end=554) xoi(i),xj(i)
      i=i+1
      go to 555
  554 nxy=i-1
      write(6,*) nxy,xxmin,yymin
      do i=1,nxy
      xi(i)=(xoi(i)-xxmin)*cos((yymin+ty/2)*3.14159/180.)
      xj(i)=xj(i)-yymin
      end do
    7 format(2f7.2,d20.12)
      close(2)
      ll=nn
      open(2, file='hist-etas-mk.upda')
      read(2,*) w1,w2,fmin,mx
      read(2,*) a1,a2,a3,aa4,a5,a6,a7
      read(2,*) (b(i),i=1,mx)
      close(2)
      a4=aa4/2
      bmin=1.d10
      bmax=-1.d10
      do 210 i=1,mx/2
        b1(i)=a1*exp(b(i))
        b2(i)=a2*exp(b(i+mx/2))
        b4(i)=a4
        b5(i)=a5
        b7(i)=a7
        bv1(i)=exp(bv(i))
  210 continue
c     ymin=ymind
c     write(6,*) 'xmin',xmin,ymin
c     open(3,file='update7pa.data')
      open(3,file='work.para')
      n1=1
c 220 read(3,*,end=230) nn1,along,alatt,amg,xt,dd,ee,ff,
c    &           aklog1,bv2
  220 read(3,*,end=230)  amulog1,aklog1,
     &                  along,alatt,amg,xt,dd,ee,ff
      xx1(n1)=along-xxmin
      xx1(n1)=xx1(n1)*cos((yymin+ty/2)*3.14159/180.)
      yy1(n1)=alatt-yymin
      zz1(n1)=xt
      zmg1(n1)=amg-xmg0
      aa1(n1)=ee/dd/sqrt(1-ff**2)
      bb1(n1)=ff   /sqrt(1-ff**2)
      cc1(n1)=dd/ee/sqrt(1-ff**2)
      a11(n1)=a1*exp(amulog1)
      bb2(n1)=a2*exp(aklog1)
      bb4(n1)=a4
      bb5(n1)=a5
      bb7(n1)=a7
c     bbv(n1)=exp(bv2)
      n1=n1+1
      go to 220
  230 n1=n1-1
      xlmin=1.d10
      xlmax=-1.d10
c     open(4,file='img7pa.data')
      open(9,file='img1.data')
      iday1=1
      nday=1
      open(19,file='time.conf')
      read(19,*) tday1
c     read(19,*) nday
      nday=1
      close(19)
      write(6,*) iday1,nday
      do jday=iday1,iday1+nday-1
c     zzz=45290.5d0+jday
c 2011.4.11
      zzz=tz+tday1
      xlmin=1.d10
      xlmax=-1.d10
      do 100 i=1,nxy
      xxx=xi(i)
      yyy=xj(i)
c     bvi=exp(bvxy(i))
      do 320 ipt=1,ndelo
        xmin=1.d69
        ymin=1.d69
        xmax=-1.d69
        ymax=-1.d69
        do 330 k=1,3
          xd(k)=xx(ldt(ipt,k))
          yd(k)=yy(ldt(ipt,k))
          zd(k)=b1(ldt(ipt,k))
          if(xd(k).lt.xmin) xmin=xd(k)
          if(xd(k).gt.xmax) xmax=xd(k)
          if(yd(k).lt.ymin) ymin=yd(k)
          if(yd(k).gt.ymax) ymax=yd(k)
  330     continue
           if(xxx.lt.xmin) go to 340
           if(xxx.gt.xmax) go to 340
           if(yyy.lt.ymin) go to 340
           if(yyy.gt.ymax) go to 340
            denom=det(xd(1),xd(2),xd(3),yd(1),yd(2),yd(3))
      if(denom.eq.0.0) write(6,*) 'denominator=0.0 for ipt =', ipt
      if(denom.eq.0.0) go to 340
            ax1=det(xxx,xd(2),xd(3),yyy,yd(2),yd(3))/denom
            ax2=det(xd(1),xxx,xd(3),yd(1),yyy,yd(3))/denom
            ax3=det(xd(1),xd(2),xxx,yd(1),yd(2),yyy)/denom
      if(ax1.lt.eps.or. ax2.lt.eps.or. ax3.lt.eps) go to 340
            bx1=ax1*zd(1)+ax2*zd(2)+ax3*zd(3)
            go to 325
  340     continue
  320 continue
  325 continue
      rmdi=0.0
      do 120 l=1,nn+n1
      if(l.le.nn) then
      xxi=xx(l)
      yyi=yy(l)
      zzi=zz(l)
      zmgi=zmg(l)
      aai=aa(l)
      bbi=bb(l)
      cci=cc(l)
      b2i=b2(l)
      b4i=b4(l)
      b5i=b5(l)
      b7i=b7(l)
      else
      jj=l-nn
      xxi=xx1(jj)
      yyi=yy1(jj)
      zzi=zz1(jj)
c     if(zzi.gt.zzz-1.d0) go to 130
c 2011.4.11
      if(zzi.gt.zzz) go to 130
      zmgi=zmg1(jj)
      aai=aa1(jj)
      bbi=bb1(jj)
      cci=cc1(jj)
      b2i=bb2(jj)
      b4i=bb4(jj)
      b5i=bb5(jj)
      b7i=bb7(jj)
      end if
      xxij=(xxx-xxi)*exp(-b4i*zmgi)
      yyij=(yyy-yyi)*exp(-b4i*zmgi)
      space=(aai*xxij**2-2*bbi*xxij*yyij+
     &       cci*yyij**2+a6)**(-b7i)
      rmdi=rmdi+b2i/(zzz-zzi+a3)**b5i *space
  120 continue
  130 continue
      l1=l-1
c     write(6,*) ' l1',l1,jj,bx1,rmdi
      ramdai=bx1+rmdi
c     write(6,*) 'mu ',l,jj,bx1,a11(jj)
      if(ramdai.gt.0.9999) then
      write(6,*) j,i,ramdai,rmdi,b4i
c     write(6,*) xxx,xoi(j)+120.,yyy+20.,zzz
      endif
      sum=0
c     jm=51
c     do im=1,jm
c     bmag=4.d0+0.1d0*(im-1)
c     xl(i,im)=-bvi*(bmag-3.95)+log10(1.d0-10.d0**(-bvi*0.1))+
c    &        log10(ramdai)
c     if(xlmin.gt.xl(i,im)) xlmin=xl(i,im)
c     if(xlmax.lt.xl(i,im)) xlmax=xl(i,im)
c     sum=sum+10.d0**(xl(i,im))
c     end do
      xl1(i)=log10(ramdai)
      xl2(i)=log10(ramdai)
  110 continue
      if(mod(i,100).eq.0) write(6,*) i
  100 continue
  200 continue
    8 format(f6.1,f9.4,f9.4,f4.1,f10.5)
      write(6,*) 'max',jday,l1,l1-nn,xlmax
        sum=0
c       do i=1,nxy
c       yyy=xj(i)
c       xxx=xoi(i)
c       do im=1,jm
c       bmag=4.d0+0.1d0*(im-1)
c       sum=sum+10.**(xl(i,im)-log10(100.))
c       write(4,8) zzz-45291.d0,xxx,yyy,bmag,xl(i,im)-log10(100.)
c       end do
c       end do
        sumx=0
        do i=1,nxy
        yyy=xj(i)+yminnode
        xxx=xoi(i)
!       sumx=sumx+10.**(xl2(i)-log10(100.))
        sumx=sumx+10.**(xl2(i))
!       write(9,8) zzz-tz,xxx,yyy,bmag,xl2(i)-log10(100.) !#events/(0.1deg)^2/day
        write(9,88) zzz-tz,xxx,yyy,xl2(i)
   88 format(f11.5,f9.4,f9.4,f10.5)
   98 format(f6.1,f9.4,f9.4,f10.5)
        end do
        write(6,*) 'sum',sum,sumx
      end do
      close(4)
      close(9)
      stop
      end
      subroutine redata(xx,yy,zz,zmg,tx,ty,tz,nn,xmin,ymin)
      implicit real * 8 (a-h,o-z)
      parameter(npa=12345,nda=12345,mdl=23333,nr=020)
      parameter(npall=npa*5)
      common /coef/ aa(nda),bb(nda),cc(nda)
      common /triang/ adt(mdl),ldt(mdl,3),ndt(mdl,3),np,npext,ndelo
      common /delink/ nzu(npa),iu(nr,npa)
      common /xad/xd(npall),m
      common /jfun/rmax(nda),jxy(nda,nda),jnn(nda)
      common /approx/ xmg0
      dimension xx(1),yy(1),zz(1),zmg(1)
      character*8 fmt(10)
      character*80 maindata
      common /range/ tstart
      open(1,file='hist-etas-mk.conf')
      read(1,112) maindata
  112 format(a)
c     read(1,1) nfunct
    1 format(i10)
c     read(1,2) fmt
    2 format(10a8)
      read(1,*) tx,ty,tz,nn
    3 format(3f10.1,i10)
      read(1,*) xmin,ymin,xmg0,zmin,tstart,bi2
      write(6,112) maindata
      write(6,2) fmt
      write(6,3) tx,ty,tz,nn
      write(6,8) xmin,ymin,xmg0,zmin,tstart
    8 format(8f10.0)
      in=10
      open(in,file=maindata)
      read(in,*) np,npext,ndelo,tx,ty
    7 format(1h ,'data set  ',10a8)
      m=npext*5
      nnc=0
      do 20 i=1,npext
      read(in,*,end=25) idummy,xx(i),yy(i),zmg(i),zz(i),dd,ee,ff
      if(i.gt.np) go to 20
      aa(i)=ee/dd/sqrt(1-ff**2)
      bb(i)=ff   /sqrt(1-ff**2)
      cc(i)=dd/ee/sqrt(1-ff**2)
   11 format(i5,4f8.4)
   20 continue
      i=nn+1
   25 nn=i-1
      do 509 i=1,ndelo
      read(in,*) idummy,(ldt(i,j),j=1,3),adt(i)
  509 continue
      do 520 i = 1, npext
        read(in,*) idummy, nzu(i), (iu(j,i), j=1,nzu(i))
  520 continue
      ii=0
      do 10 i=1,nn
      ii=ii+1
       xx(ii)=xx(i)
       yy(ii)=yy(i)
       zmg(ii)=zmg(i)-xmg0
       aa(ii)=aa(i)
       bb(ii)=bb(i)
       cc(ii)=cc(i)
       if(zz(ii).gt.tstart) nnc=nnc+1
       xx(ii)=xx(ii)*cos((ymin+ty/2)*3.14159/180.)
   10 continue
      tx=tx*cos((ymin+ty/2)*3.14159/180.)
      do 40 i=1,ndelo
        adt(i)=adt(i)*cos((ymin+ty/2)*3.14159/180.)
   40 continue
      nn=ii
      write(6,4) in
      write(6,*) 'nfunct=',nfunct
      write(6,5) tx,ty,tz,xmin,ymin,xmg0,zmin,tstart,nn,nnc
    4 format(1h ,'input device',i10)
    5 format(1h0,'tx,ty,tz,xmin,ymin,xmg0,zmin,tsta'/8f10.3,
     &/' nn =',i5,' nnc =',i5)
    6 format(///1h ,i5,4f13.5,i5,4f13.5/(1h ,i5,4f13.5,i5,4f13.5))
      close(1)
      return
      end
