Program dsdcal

           parameter(nx=1000000,na=130)
           parameter(ny=100)
           
           Real pippsd(nx),corcal(nx)
           Real psdoriginal(nx,na)
           Real psdsave(ny),psds(nx,ny)

           Integer kyear,jday(nx),ka(nx),ma(nx)
           Integer jd(nx),kr(nx),mn(nx)
           Integer ipsd(nx,ny),isay(ny)
   
           double precision  dmid(na)

           Character*80 infile,infile2
           Character*80 infile0,infile1
           Character*10,y1,y2,y3,y4,y5,y6,y7,y8,y9

           double precision dmm,a1,a2

           open(10,file='filepsdadjustment',status='unknown')
           open(20,file='pip002_diameter_table_2023.txt',
     *     status='unknown')

           read(10,11)infile,infile2
c          read(10,11)infile0,infile1
 11        format(a35,1x,a30)

           open(12,file=infile,status='unknown')
           open(11,file=infile2,status='unknown')
c          open(14,file=infile0,status='unknown')
c          open(13,file=infile1,status='unknown')
  
           k2=1
           do i=1,nx
           read(12,*,err=998,end=998)iy,jd(k2),kr(k2),mn(k2),
     *     ktot,(ipsd(k2,k),k=1,ny)
           k2=k2+1
 998       enddo
c          do i=1,nx
c          read(14,*,err=997,end=997)iy,jd(k2),kr(k2),mn(k2),
c    *     ktot,(ipsd(k2,k),k=1,ny)
c          k2=k2+1
c997       enddo
c          k2=k2-1

           print*,k2

           do k=1,ny
           psdsave(k)=0.
           isay(k)=0
           enddo

           do i=1,9
           read(11,*)
c          read(13,*)
           enddo
           read(11,*)y1,y2,y3,y4,y5,y6,y7,y8,(dmid(k),k=1,na)
c          read(13,*)

           k1=1
           ksay=0
           do i=1,nx
           read(11,*,err=999,end=999)kyear,jday(k1),mon,iday,
     *     ka(k1),ma(k1),itot,y9,(psdoriginal(k1,k),k=1,na)
           do j=1,k2
           if(jday(k1).eq.jd(j) .and. ka(k1).eq.kr(j) .and.
     *     ma(k1).eq.mn(j))then
           ksay=ksay+1
           do k=1,41
           if(ipsd(j,k).eq.1 .and. psdoriginal(k1,k).gt.0.)then
           isay(k)=isay(k)+1
           psdsave(k)=psdoriginal(k1,k)
           psds(isay(k),k)=psdoriginal(k1,k)
           endif
           enddo
           do k=42,ny
           if(psdoriginal(k1,k).gt.0.)then
           isay(k)=isay(k)+1
           psdsave(k)=psdoriginal(k1,k)
           psds(isay(k),k)=psdoriginal(k1,k)
           endif
           enddo
           endif
           enddo
           k1=k1+1
 999       enddo
c          do i=1,nx
c          read(13,*,err=996,end=996)kyear,jday(k1),mon,iday,
c    *     ka(k1),ma(k1),itot,y9,(psdoriginal(k1,k),k=1,na)
c          do j=1,k2
c          if(jday(k1).eq.jd(j) .and. ka(k1).eq.kr(j) .and.
c    *     ma(k1).eq.mn(j))then
c          ksay=ksay+1
c          do k=1,41
c          if(ipsd(j,k).eq.1 .and. psdoriginal(k1,k).gt.0.)then
c          isay(k)=isay(k)+1
c          psdsave(k)=psdoriginal(k1,k)
c          psds(isay(k),k)=psdoriginal(k1,k)
c          endif
c          enddo
c          do k=42,ny
c          if(psdoriginal(k1,k).gt.0.)then
c          isay(k)=isay(k)+1
c          psdsave(k)=psdoriginal(k1,k)
c          psds(isay(k),k)=psdoriginal(k1,k)
c          endif
c          enddo
c          endif
c          enddo
c          k1=k1+1
c996       enddo
c          k1=k1-1

           print*,k1,k2,ksay

           do k=1,ny
           is=isay(k)
           write(16,30)dmid(k),psdsave(k),(psds(kk,k),kk=1,is)
 30        format(f10.1,1000f15.3)
           enddo


           frame=380.*60
           a1=0.0d0
           a2=0.2d0
           do k=1,ny
           dmm=dmid(k)
           dof=0.117/2.*dmm
           fov=1.e-6*(64.-dmm)*(48.-dmm)
           fov2=1.e-6*(32.-dmm)*(24.-dmm)
           vol=frame*dof*fov
           vol1=frame*dof*fov2
           psd=1./(vol*0.2)
           psd1=1./(vol1*0.2)
           cor=psdsave(k)/psd
           write(20,33)dmm,psdsave(k),psd,cor
 33        format(f8.1,4f15.5)
           a1=a1+0.2d0
           a2=a2+0.2d0
           enddo

           stop
           end