	program delpch1fun
c .......................................................Nov 2019
C Add data/model assimilation using chain power model(obser)
C
C Combine Pch+ 10 kinds of indices
C
C
C Select peaks of delPch=Pchmax(3d)-Pchmin(-3d) for del > 100 i.u.
C movind 3day window, 
C Select all Del >0: output => 'delpchpos.res'
C for 2011-2018
C
c
	DIMENSION ryears(303),del(303),ivsw(303),IAE(303),iaa(303)
	+,iap(303),fx(500),fy(500),ifx(500),ify(500)
     +,IDST(303),wp(303),iarn(303),iarp(303),iann(303),ianp(303)
C	+,pmin(303),pmax(303)
      INTEGER*2 lag(303,10)
	CHARACTER(10) DD
      CHARACTER(10) TT
      CHARACTER(5) ZZ
	CHARACTER*94,infile1,infile2,outfile1,outfile2,txt
	CHARACTER*2 AMD(303,2)
	CHARACTER*16 txt1(303)
C
	CALL DATE_AND_TIME(DATE=DD,TIME=TT,ZONE=ZZ)
      TT(5:)='      '
      WRITE(*,*) 'PC Date: Year,Month,Day = ',DD,'Time = ',TT
C
  791 FORMAT(A10,1X,I4)
C      outfile1='catpch2thr.res' !  only Pch
      outfile1='catpch2thd.res' !   Pch/ICME
	OPEN(12,FILE=OUTFILE1)		! data for 10 indices combined
C	outfile2='extpch2fun.res' !  only Pch
	outfile2='extpch2fud.res' !  Pch/ICME
	OPEN(13,FILE=OUTFILE2,access='APPEND')  ! data + model extra points
C
C Year    MN DD  Pchmin  Pchmax  delPch  d Vswmax  d AEmax  d aamax  d apmax  d Dstmin d Wpmax nn  d ICME
C 2011.211  3 19  235.94  446.21  210.27  4   516  4   587  4    45  3    18  3   -17  4  3.7    1
   36	format(1X,F8.3,2(1X,A2),A16,1X,F7.2,5(1X,I2,2X,I4),1X,I2,2X,F4.1)
C Year     MN DD  Pchmin  Pchmax  delPch +d   ARn +d   ARp +d   ANn +d   ANp   nn
C 2011.211  3 19  235.94  446.21  210.27  5   159  1   119  3   240  4   265    1
   37	format(39X,4(1X,I2,2X,I4))
   38	format(1X,F8.3,2(1X,A2),A16,1X,F7.2,5(1X,I2,2X,I4),1X,I2,2X,F4.1
     +,4(1X,I2,2X,I4),1X,I4)
   30	format(A94)

C
C+	infile1='catpchind5d-.txt' ! only Pch
CC      	infile1='extpchcat4n-.txt' ! only Pch 6-param
C+	infile1='catpchICME5d.txt'	 ! Pch/ICME
      	infile1='extpchcat4d.txt' ! Pch/ICME
		OPEN(10,FILE=infile1,ACTION='READ')
	read(10,30) txt !~ title line
	write(12,30) txt
C	infile2='catpcharcICME4.txt'
C      infile2='extpcharc4n-.txt' ! only Pch AR+AN
       infile2='extpcharc4d.txt' ! Pch/ICME AR+AN

	OPEN(11,FILE=infile2,ACTION='READ')
	read(11,30) txt !~ title line
	write(12,30) txt
C
 	icntot=0 ! count of dPch
	do 101 ii=1,303
	read(10,36,end=102) RYEARS(ii),AMD(ii,1),AMD(ii,2),txt1(ii)
	+,del(ii),lag(ii,1),ivsw(ii),lag(ii,2),IAE(ii),lag(ii,3),iaa(ii)
     +,lag(ii,4),iap(ii),lag(ii,5),IDST(ii),lag(ii,6),WP(ii)
      read(11,37,end=102)	lag(ii,7),iarn(ii),lag(ii,8),iarp(ii)
	+,lag(ii,9),iann(ii),lag(ii,10),ianp(ii)
      icntot=icntot+1
  101	CONTINUE
  102  CLOSE(unit=10)
      CLOSE(unit=11)
CREM	GOTO 103         ! TEMPORARY!!!!!!!!!!!!!!!!!!
C--------------------------------------------------------------------
C
      icnt=0
      do 100 ii=1,icntot 
	icnt=icnt+1
	write(*,38) RYEARS(ii),AMD(ii,1),AMD(ii,2),txt1(ii)
	+,del(ii),lag(ii,1),ivsw(ii),lag(ii,2),IAE(ii),lag(ii,3),iaa(ii)
     +,lag(ii,4),iap(ii),lag(ii,5),IDST(ii),lag(ii,6),WP(ii)
	+,lag(ii,7),iarn(ii),lag(ii,8),iarp(ii),lag(ii,9),iann(ii)
     +,lag(ii,10),ianp(ii),icnt
	write(12,38) RYEARS(ii),AMD(ii,1),AMD(ii,2),txt1(ii)
	+,del(ii),lag(ii,1),ivsw(ii),lag(ii,2),IAE(ii),lag(ii,3),iaa(ii)
     +,lag(ii,4),iap(ii),lag(ii,5),IDST(ii),lag(ii,6),WP(ii)
	+,lag(ii,7),iarn(ii),lag(ii,8),iarp(ii),lag(ii,9),iann(ii)
     +,lag(ii,10),ianp(ii),icnt
C
  100	CONTINUE
  103  CLOSE(unit=12)
C
C Extrapolation
C
CC Del Vsw(dPch)
  104 write(*,*) '% dPch  Vsw' 
	write(13,*) '% dPch  Vsw' 
      nn=icntot
      do kk=1,nn
	fx(kk)=del(kk)
	ify(kk)=ivsw(kk)
	enddo
C		  
      nbeg=nn+1
      IF (icntot.gt.200) THEN	 
C only PCH x= 90: 800 step 5	 n=142
      nend=nbeg+141
	x=85.
	step=5.
       a = 1.301
       b = 0.7853
       c = 479.4
                            ELSE
C PCH/ICME x=90: 500 step	20 n=21
      nend=nbeg+20
	x=70.
	step=20.
       a = 141.7
       b = 0.2562
	 c = 0.
	                    ENDIF
C
         do kk=nbeg,nend
  	x=x+step
	fx(kk)=x
	y = a*x**b+c
      ify(kk)=nint(y)
	enddo
C output
   40 format(1X,F7.2,2X,I4)
       do kk=1,nend
	write(*,40) fx(kk),ify(kk)
	write(13,40) fx(kk),ify(kk)
	 enddo
CC
CC Vsw(dPch) AE
  105 write(*,*) '% Vsw   AE' 
	write(13,*) '% Vsw   AE' 
      nn=icntot
      do kk=1,nn
	ifx(kk)=ivsw(kk)
	ify(kk)=iae(kk)
	enddo
C		  
      nbeg=nn+1
C		 
      IF (icntot.gt.200) THEN	 
C only PCH x= 350: 900 step 5	 n=110 Pch 
      nend=nbeg+109
	step=5.
	x=345.
       a = 12.76 
       b = 0.7244 
       c = -459.6 
                            ELSE
C PCH/ICME  x= 350: 800 step 20 n=22 PCH/ICME
      nend=nbeg+21
	step=20.
	x=330.
       a = 48.39
       b = 0.5753
       c = -829.9
	                    ENDIF
C
         do kk=nbeg,nend
  	x=x+step
	ifx(kk)=nint(x)
	y = a*x**b+c
      ify(kk)=nint(y)
	enddo
C output
   41 format(2(2X,I4))
       do kk=1,nend
	write(*,41) ifx(kk),ify(kk)
	write(13,41) ifx(kk),ify(kk)
	 enddo
CC
CC AE aa
  106 write(*,*) '% AE  aa' 
	write(13,*) '% AE  aa' 
      nn=icntot
      do kk=1,nn
	ifx(kk)=iae(kk)
	ify(kk)=iaa(kk)
	enddo
C		  
      nbeg=nn+1
      IF (icntot.gt.200) THEN	 
C only PCH x= 200: 1800 step 10	 n=160
      nend=nbeg+159
	x=190.
	step=10.
       a =   5.371e-05
       b =       2.015
       c =       23.41
                            ELSE
C PCH/ICME x=200: 1800 step 80    n=21
      nend=nbeg+20
	x=120.
	step=80.
       a =     0.00163
       b =       1.543
       c =       17.56
	                    ENDIF
C
      do kk=nbeg,nend
  	x=x+step
	ifx(kk)=nint(x)
	y = a*x**b+c
      ify(kk)=nint(y)
	enddo
C output
       do kk=1,nend
	write(*,41) ifx(kk),ify(kk)
	write(13,41) ifx(kk),ify(kk)
	 enddo
CC
CC aa  ap
  107 write(*,*) '% aa  ap' 
	write(13,*) '% aa  ap' 
      nn=icntot
      do kk=1,nn
	ifx(kk)=iaa(kk)
	ify(kk)=iap(kk)
	enddo
C		  
      nbeg=nn+1
      IF (icntot.gt.200) THEN	 
C only PCH x= 0: 350 step 5	 n=70
      nend=nbeg+69
	x=-5.
	step=5.
       a = 0.3032
       b = 1.143
       c = 0.
                            ELSE
C PCH/ICME x=0: 240 step 10    n=25
      nend=nbeg+24
	x=-5.
	step=10.
       a = 0.2927
       b = 1.208
	 c = 0.
	                    ENDIF
C
      do kk=nbeg,nend
  	x=x+step
	ifx(kk)=nint(x)
	y = a*x**b+c
      ify(kk)=nint(y)
	enddo
C output
       do kk=1,nend
	write(*,41) ifx(kk),ify(kk)
	write(13,41) ifx(kk),ify(kk)
	 enddo
CC
CC ap  Dst
  108 write(*,*) '% ap  Dst' 
	write(13,*) '% ap  Dst' 
      nn=icntot
      do kk=1,nn
	ifx(kk)=iap(kk)
	ify(kk)=idst(kk)
	enddo
C		  
      nbeg=nn+1
      IF (icntot.gt.200) THEN	 
C only PCH x= 0: 240 step 5	 n=49
      nend=nbeg+48
	x=0.
	step=5.
       a = -1.806 
       b =  0.8432
       c = 0.
                            ELSE
C PCH/ICME x=0: 180 step 10    n=19
      nend=nbeg+18
	x=-10.
	step=10.
       a = -2.013 
       b =  0.8437
	 c = 0.
	                    ENDIF
C
      do kk=nbeg,nend
  	x=x+step
	ifx(kk)=nint(x)
	y = a*x**b+c
      ify(kk)=nint(y)
	enddo
C output
       do kk=1,nend
	write(*,41) ifx(kk),ify(kk)
	write(13,41) ifx(kk),ify(kk)
	 enddo
CC
CC Dst  Wp
  109 write(*,*) '% Dst  Wp' 
	write(13,*) '% Dst  Wp' 
      nn=icntot
      do kk=1,nn
	ifx(kk)=-idst(kk)
	fy(kk)=wp(kk)
	enddo
C		  
      nbeg=nn+1
      IF (icntot.gt.200) THEN	 
C only PCH x= 0: 240 step 5	 n=49
      nend=nbeg+48
	x=0.
	step=5.
       a = 0.5668
       b = 0.4482
       c = 1.852 
                            ELSE
C PCH/ICME x=0: 240 step 10    n=24
      nend=nbeg+23
	x=-5.
	step=10.
       a = 0.1115
       b = 0.7507
       c = 3.271
	                    ENDIF
C
      do kk=nbeg,nend
  	x=x+step
	ifx(kk)=nint(x)
	y = a*x**b+c
      fy(kk)=y
	enddo
C output
   42 format(1X,I4,2X,F5.1)
       do kk=1,nend
	write(*,42) ifx(kk),fy(kk)
	write(13,42) ifx(kk),fy(kk)
	 enddo
CC
CC Wp  ARn
  110 write(*,*) '% Wp  ARn' 
	write(13,*) '% Wp  ARn' 
      nn=icntot
      do kk=1,nn
	fx(kk)=wp(kk)
	fy(kk)=iarn(kk)/10.
	enddo
C		  
      nbeg=nn+1
      IF (icntot.gt.200) THEN	 
C only PCH x= 2.5: 8.5 step 0.05	 n=120
      nend=nbeg+119
	x=2.45
	step=0.05
       a = 1.228
       b = 1.849
       c = 0.
                            ELSE
C PCH/ICME x=2.5: 8.5 step 0.3    n=20
      nend=nbeg+19
	x=2.2
	step=0.3
       a = 4.068
       b = 1.358
       c = 0.
	                    ENDIF
C
      do kk=nbeg,nend
  	x=x+step
	fx(kk)=x
	y = a*x**b+c
      fy(kk)=y
	enddo
C output
   43 format(2(2X,F6.1))
       do kk=1,nend
	write(*,43) fx(kk),fy(kk)
	write(13,43) fx(kk),fy(kk)
	 enddo
CC
CC Wp  ARp
  111 write(*,*) '% Wp  ARp' 
	write(13,*) '% Wp  ARp' 
      nn=icntot
      do kk=1,nn
	fx(kk)=wp(kk)
	fy(kk)=iarp(kk)/10.
	enddo
C		  
      nbeg=nn+1
      IF (icntot.gt.200) THEN	 
C only PCH x= 2.5: 8.5 step 0.05	 n=120
      nend=nbeg+119
	x=2.45
	step=0.05
       a = 2.631
       b = 1.564
       c = 0.
                            ELSE
C PCH/ICME x=2.5: 10.0 step 0.5    n=15
      nend=nbeg+14
	x=2.0
	step=0.5
       a = 6.979 
       b = 0.9323
       c = 0.
	                    ENDIF
C
      do kk=nbeg,nend
  	x=x+step
	fx(kk)=x
	y = a*x**b+c
      fy(kk)=y
	enddo
C output
       do kk=1,nend
	write(*,43) fx(kk),fy(kk)
	write(13,43) fx(kk),fy(kk)
	 enddo
CC
CC Wp  ANn
  112 write(*,*) '% Wp  ANn' 
	write(13,*) '% Wp  ANn' 
      nn=icntot
      do kk=1,nn
	fx(kk)=wp(kk)
	fy(kk)=iann(kk)/10.
	enddo
C		  
      nbeg=nn+1
      IF (icntot.gt.200) THEN	 
C only PCH x= 2.5: 8.5 step 0.05	 n=120
      nend=nbeg+119
	x=2.45
	step=0.05
       a = 3.356
       b = 1.398
       c = 0.
                            ELSE
C PCH/ICME x=2.5: 9.5 step 0.5    n=14
      nend=nbeg+13
	x=2.0
	step=0.5
       a =       4.156
       b =       1.442
	 c = 0.
	                    ENDIF
C
      do kk=nbeg,nend
  	x=x+step
	fx(kk)=x
	y = a*x**b+c
      fy(kk)=y
	enddo
C output
       do kk=1,nend
	write(*,43) fx(kk),fy(kk)
	write(13,43) fx(kk),fy(kk)
	 enddo
CC
CC Wp  ANp
  113 write(*,*) '% Wp  ANp' 
	write(13,*) '% Wp  ANp' 
      nn=icntot
      do kk=1,nn
	fx(kk)=wp(kk)
	fy(kk)=ianp(kk)/10.
	enddo
C		  
      nbeg=nn+1
      IF (icntot.gt.200) THEN	 
C only PCH x= 2.5: 8.5 step 0.05	 n=120
      nend=nbeg+119
	x=2.45
	step=0.05
       a = 3.603 
       b = 1.47  
       c = 0.
                            ELSE
C PCH/ICME x=2.5: 10.0 step 0.5    n=15
      nend=nbeg+14
	x=2.0
	step=0.5
       a = 16
       b = 0.7125
       c = -13.7
	                    ENDIF
C
      do kk=nbeg,nend
  	x=x+step
	fx(kk)=x
	y = a*x**b+c
      fy(kk)=y
	enddo
C output
       do kk=1,nend
	write(*,43) fx(kk),fy(kk)
	write(13,43) fx(kk),fy(kk)
	 enddo
CC
      close(unit=13)
  900 pause ' '
      STOP
      END

