      program cntpch3d
C...............................................Jan. 2019
C
C cnt delpch >0 delpch(i)=Pchmax(3d_i)-Pchmin(-3d_i-1)
C to produce histogram of delpch >0
C
      dimension pch(3),pch_pre(3),pchs(3214),iyrs(3214,3)
      character*50 infile,outfile
      character*2 AMN,ADY

	infile='pchdata.dat'
	outfile='delpch3d.txt'	 ! all min(-3d),max(+3d)+- delpch date:2+day
	OPEN(10,file=infile,action='READ')
	OPEN(11,file=outfile)
	nn=0
C
      DO ii=1,3214
	read(10,30,end=51) iyr,iMN,iDY,pchdy
	iyrs(ii,1)=iyr
      iyrs(ii,2)=imn
	iyrs(ii,3)=idy
	pchs(ii)=pchdy
	icnt=ii
	ENDDO
C
   51   pcmax_pre=0.0
	nn=0
      DO 44 ii=1,icnt
      pch_pre(1)=pchs(ii)
      pch_pre(2)=pchs(ii+1)
	pch_pre(3)=pchs(ii+2)
      pch(1)=pchs(ii+3)
      pch(2)=pchs(ii+4)
	pch(3)=pchs(ii+5)
	 pcmin=AMIN1(pch_pre(1),pch_pre(2),pch_pre(3))
	pcmax=AMAX1(pch(1),pch(2),pch(3))
	do k=1,3
	if (pch(k).eq.pcmax) then
	nk=k+ii+2
	exit
	endif
	enddo
	del=pcmax-pcmin
	iyr=iyrs(nk,1)
	imn=iyrs(nk,2)
	idy=iyrs(nk,3)
      lda=ndoy(iyr,imn,idy)
	ryr=float(iyr)
		if(int(ryr/4.)*4.eq.iyr) THEN 
	dnr=366.
	                          ELSE
	dnr=365.
	ENDIF
      delyr=1./dnr
	ryear=2000.+ryr+(lda-1)*delyr
	call blet2(imn,AMN)
      call blet2(idy,ADY)
	if (pcmax.ne.pcmax_pre) then
	nn=nn+1
	write(*,31) ryear,AMN,ADY,pcmin,pcmax,del,nn
	write(11,31) ryear,AMN,ADY,pcmin,pcmax,del,nn
	pcmax_pre=pcmax
	endif
   44 CONTINUE
   50 close(unit=10)
       close(unit=11)
   30 format(1X,3I2,1X,F7.2)
   31 format(1X,F8.3,2(1X,A2),3(1X,F7.2),1X,I4)
      pause ' '
	STOP
	END