      program IMFpchmax
C...............................................Jan. 2019
C
C Fix days of Pch peaks 2011-2018 at IMF sector boundaries <+> towards the Earth
C Similar for Wp storm peaks
C Svalgaard IMF_SB list infilesb = 'IMF-sector.txt'
C
      dimension wps(840),iyrs(840,3),years(840)
	character*50 infile,outfile,infilesb
	character*2 s1,s2,s1_nex,s2_nex ! <+> towards Earth, <-> towards Sun
C
C+	infile='pchmax.txt'
C+	outfile='imfpchmax.txt'	 ! Pch peaks at SB <+>
	infile='wpstmax1.txt'
	outfile='imfwpmax.txt'	 ! Wp storm peaks at SB <+>
c
	infilesb = 'imfsectorpch.txt'
C
	OPEN(12,file=infilesb,action='READ')
C
	OPEN(10,file=infile,action='READ')
C
	OPEN(11,file=outfile,access='APPEND')
C
   28 format(A18)
   29 format(2(1X,A2),2X,I4,2(1X,I2))
   30 format(21X,I2,3(1X,I2),1X,F4.1)
   31 format(3X,I2,5X,I2,2(1X,I2),1X,F4.1)
C   32 format(1X,F8.3,2(1X,I2),1X,F7.2,18X,I4)	! Pch
   32 format(1X,F8.3,2(1X,I2),4X,F4.1,1X,I4)	! Wpmax
   33 format(1X,F8.3,2(1X,I2),1X,F7.2,2(1X,I4))
C
      ii=0
      DO 49 i=1,632
c	read(10,*) txt
	read(10,32,end=50) xyear,imn,idy,pch,ks
	iyear=int(xyear)
C--	if (iyear.lt.1997) goto 49
	ii=ii+1
	years(ii)=xyear
	iyrs(ii,1)=iyear-iyear/100*100
	iyrs(ii,2)=imn
	iyrs(ii,3)=idy
	wps(ii)=pch
   49	continue
   50 close(unit=10)
      iev=ii
C		  
	nm=0  ! number of Pch peak dates
	nres=0 ! number of Pch(<+>) peak dates
C
C Read extra lines IMF_SB
c    1 	read(12,29,end=45) txt
    1	read(12,29,end=45) s1,s2,iyearsb,imnsb,idysb
c      s1=txt(1:1)
c	s2=txt(3:3)
      if (iyearsb.lt.2011) goto 1		! extra lines
	IF (s2.eq.'+1') THEN
      read(12,29) s1_nex,s2_nex,iyearsb_nex,imnsb_nex,idysb_nex
	   if (s1_nex.eq.'+1') then
	   goto 5
	   else
         backspace(12)
	   goto 1
	   endif
	ELSE
	goto 1
	ENDIF
C End pair lines for IMF_SB <+> sector
   5	zyearsb=float(iyearsb)
		if(int(zyearsb/4.)*4.eq.iyearsb) THEN 
	dnrsb=366.
	                          ELSE
	dnrsb=365.
	ENDIF
	delyrsb=1./dnrsb
	iyrcb=iyearsb-iyearsb/100*100
	ldasb=ndoy(iyrsb,imnsb,idysb)
	ryearsb=zyearsb+delyrsb*(ldasb-1)
CC
	zyearsb_nex=float(iyearsb_nex)
		if(int(zyearsb_nex/4.)*4.eq.iyearsb_nex) THEN 
	dnrsb_nex=366.
	                          ELSE
	dnrsb_nex=365.
	ENDIF
	delyrsb_nex=1./dnrsb_nex
	iyrcb_nex=iyearsb_nex-iyearsb_nex/100*100
	ldasb_nex=ndoy(iyrsb_nex,imnsb_nex,idysb_nex)
	ryearsb_nex=zyearsb_nex+delyrsb_nex*(ldasb_nex-1)
C  
    2	continue
C
C
C prepare output:
C
C
   43	nm=nm+1        ! Pch event
      if (nm.gt.iev) goto 44
      ryear=years(nm)
      nyr=iyrs(nm,1)
	nmn=iyrs(nm,2)
	ndy=iyrs(nm,3)
	pch=wps(nm)
C
	if ((ryear.ge.ryearsb).and.(ryear.lt.ryearsb_nex)) then  !  >>>
      nres=nres+1
C+	nout=nm ! Pch
	nout=nm+492 ! Wp max
	write(*,33)  ryear,nmn,ndy,pch,nres,nout
	write(11,33) ryear,nmn,ndy,pch,nres,nout
      goto 43
      	else											  !
	if (ryear.lt.ryearsb) then
	goto 43
	else
             backspace(12)
	nm=nm-1
	goto 1
	endif
	endif		
C
   44 CONTINUE
   45  close(unit=11)
      close(unit=12)
      pause ' '
	STOP
	END