C ///////////////////////////////////////////////////////////////      
	 SUBROUTINE SUBDEVGXC(ist)
C ..........................................................June 2009
C Enter station number <ist> as formal parameter
C
C ........................................................Aug. 2007
C Move median calculation into separate subroutine MEDIAN................
C Apply all calculations to fnF2 reduced by XHI-function: ind,dev,dis,mpr
C then restore foF2 from fnF2 (and foF2med from fnF2med).................
C Accordingly avoid median calibration by ITU-R........................
C Include <status.txt> file for set of stations....................
C
C ........................................................Mar. 2007
C Cycle on stations/data
C
C Corrected                                              .......Dec. 2006
C Change Criterion for selecting .dis periods:
C  (1) average sum of positive indices >=3
C  (2) average sum of negative indices <=-3
C
C Data R12 are entered for given year RZS(12 months) from <rzyr.dat> file 

C Corrected median.using IRI-CCIR fc=foF2 trend prediction......Nov. 2006

C T.L.Gulyaeva .................................................Sep. 2006

C Array fc(0:12,0:24) includes fc for mon (Dec05...Dec06), UT=0,...23 hrs
C
C DEVFF2.FOR DEVIATIONS FOR EVENT>= 3 h : +-DEV% FROM DETRENDED MEDIAN :Nm=foF2^2~ZZ !!!
C IDCE FORMAT from 27 prec. days
C REM =============== Medians amd .mpr of foF2 => indices for NmF2 =======================
C Include logarithmic scale indices -4,-3,-2,-1,0,1,2,3
C Equivalent to log10(dNe/NmF2)=-0.999,-0.2,-0.1,-0.05,0
C                              0.05,0.1,0.2,0.999  
C Ref. Gulyaeva T.L. Daily assessment of ionosphere variability,
C      Acta Geod. Geophys. Hung., 37, Nos.3-4, 160-163, 2002.
C Add selection of D+,D- intervals>3h with devNm+>40%, devNm-<-40% for IDCE
C Add regression of daily y=foF2 with x=med27 for "detrending" med27 => quiet reg-med.
C Additional output in file ***.med include : 
C XREG(1)=xav,XREG(2)=yav,XREG(3)=coef.RA (slope),XREG(4) =Daily Standard Deviation, DSD
C Additional : XREG(5) = DSD-, XREG(6) = DSD+
C
      DIMENSION AV(12,24),SD(12,24),DEV(34,0:25),ICN(12,0:24)
	+,ICD(12,0:24),DX(24),XMED(0:12,0:24),IM(0:12),YY(0:24)
     +,XLEV(12,0:24)
	+,X1(31),X2(31),R1(12),R2(12),DA(34,0:25),XREG(0:12,6)
	+,GLATS(80),GLONS(80)
	dimension TR(7),rx(24),ry(24)
	integer*4 IRES(0:12,0:25),LB3(34,4),imed(0:24)
     +,LB1(0:90,5),LB2(0:12,10),LC(5)
	real RES(0:12,0:24),RESP(0:12,0:24)
	integer*2 IY,NN,ID,JYR,JMN,JDY,mon,JD3(34,3)
	+,JFS(80),JHS(80),JTS(80),ID1,ID2
CTemp	integer*2 IDEL(12,0:24),IIN(6)     ! ID-index
	integer*4 IDEL(12,0:24),iin(6)             ! delog
	CHARACTER*128 INFILE,OUTFILE1,outfile2,outfile3,outfile4,outfile5
	CHARACTER*4 TI(12),CTIME,CYEAR,YEAR,TTT
c	+,YEARR(80)
	CHARACTER*2 ss,outmn,outyr,outdy,ww,yr,sta,CMN,CDY,outdy0
	+,AMN,STS(80),AMNS(80),DY22
	CHARACTER*6 AYMNDY
	CHARACTER(8) DDDD
      CHARACTER(10) TT
      CHARACTER(5) ZZZ
	CHARACTER*1,IN(0:9)
	CHARACTER*13 STAT,STATS(80)
C DIM RE$(12),RN$(12)
      DATA TI/'JAN.','FEB.','MAR.','APR.','MAY ','JUN.','JUL.','AUG.'
	+,'SEP.','OCT.','NOV.','DEC.'/
      DATA IM/31,31,28,31,30,31,30,31,31,30,31,30,31/
	DATA IN/'0','1','2','3','4','5','6','7','8','9'/
	data IIN/-3,-2,-1,1,2,3/           ! ID-index
c	data TR/-.2,-.1,-.05,0.,.05,.1,.2/
C Relevant dev% [-50%, -30, -10, 0, 11, 43%]
	data TR/-.301,-.155,-.045,0.,.045,.155,.301/
	COMMON /BL1/STA,YEAR,YR,AMN,DY1,DY2,DD1,DD2,ID1,ID2,JF,JH,JT,DY22
  	COMMON /BL2/ STAT,GLAT,GLON,DDDD,TTT,ZZZ,R12
	COMMON /BL3/ STATS,STS,GLATS,GLONS,JFS,JHS,II,AMNS,JTS
	COMMON /BL5/ APROX,BPROX
 	COMMON   /CONST/UMR
C
	iic=ii
	 UMR=.0174533
      PHI=3.141593

	outdy0='00'
	CALL DATE_AND_TIME(DATE=DDDD,TIME=TT,ZONE=ZZZ)
	cyear=DDDD(1:4)
	cmn=DDDD(5:6)
	cdy=DDDD(7:8)
	ctime=TT(1:4)
C
  101 continue
C Start cycle on stations:
cc     	do 102 ist=1,80
	jf=jfs(ist)

 	if (jf.eq.0) then
	 goto 102
	endif
c	YEAR=YEARR(ist)
c	YR=YEAR(3:4)
	stat=stats(ist)
 	sta=sts(ist)
	amn=amns(ist)
	write(*,*) stat,sta
200	CONTINUE
 	AYMNDY(1:2)=YR
	AYMNDY(3:4)=AMN
	AYMNDY(5:6)='01'     ! Start from 1st day of current month
cd15	FORMAT(A6)
 
      infile='c:\web\idce\'
	      INFILE(13:14)=sta
	      INFILE(15:16)=yr
	      INFILE(17:)='f.br'
c      infile='c:\web\idce\eb0500.br'
1100	OUTFILE1='c:\web\graf\dat3\07\'
	OUTFILE1(18:19)=yr
	OUTFILE2='c:\web\graf\dat2\07\'
	OUTFILE2(18:19)=yr
	OUTFILE3='c:\web\graf\tempmpr\07mpr\'
	OUTFILE3(21:22)=yr
	OUTFILE4='c:\web\graf\dat5\'
	OUTFILE5='c:\web\graf\dat4\07\'
	OUTFILE5(18:19)=yr
C+++
	outyr=yr
	read(aymndy,*) datas
	iymnid=int(datas)

 	iy=iymnid/10000
	nn=(iymnid-iy*10000)/100
	mon=nn
	id=iymnid-iy*10000-nn*100
	if (mon.eq.0) then
	iend=2
	goto 6
	endif
  
	if (iy.lt.30) then
	iyyyy=2000+iy
	              else
	iyyyy=1900+iy
	endif

     	z1=iyyyy/4.0
      jz=int(z1)*4

      IF(jz.EQ.iyyyy) THEN
               IM(2)=29
					dnr=366.
        ELSE
                IM(2)=28
	  	dnr=365.
	       ENDIF

      DO I=1,128
      IF(INFILE(I:I).EQ.'.') THEN
      J=I
      EXIT
      ENDIF
      ENDDO
	ss=infile(j-6:j-5)
	ww=infile(j-2:j-1)      ! ww='00' refers to foF2// ww='99' refers to hmF2
C
	outmn=AMN
	OUTFILE1(21:22)=sta
  	OUTFILE2(21:22)=sta
	OUTFILE3(27:28)=sta
	OUTFILE4(18:19)=sta
	OUTFILE5(21:22)=sta

	OUTFILE1(23:24)=outyr
  	OUTFILE2(23:24)=outyr
	OUTFILE3(29:30)=outyr
	OUTFILE4(20:21)=outyr
	OUTFILE5(23:24)=outyr


	OUTFILE1(25:26)=outmn
 	OUTFILE2(25:26)=outmn
	OUTFILE3(31:32)=outmn
	OUTFILE4(22:23)=outmn
	OUTFILE5(25:26)=outmn


      OUTFILE1(27:)='d.txt'
	OUTFILE2(27:)='m.txt'
	OUTFILE3(33:)='.mpr'
	OUTFILE4(24:)='c.txt'
	OUTFILE5(27:)='i.txt'

      OPEN(11,FILE=INFILE)

      OPEN(12,FILE=OUTFILE1)
	OPEN(13,FILE=OUTFILE2)
	OPEN(14,FILE=OUTFILE3) 
	OPEN(15,FILE=OUTFILE4)
	OPEN(16,FILE=OUTFILE5)
		write(*,17) stat
 	write(12,17) stat
	write(13,17) stat
	write(14,17) stat
	write(15,17) stat
	write(16,17) stat
   17	format(A13)

c-	WRITE(*,18) cyear,cmn,cdy,ctime
 	WRITE(12,18) cyear,cmn,cdy,ctime
   18	format('Date of issue: ',A4,'/',A2,'/',A2,' Time:',A4)
		write(12,201) 
  201	format('YYMMDD|UT 0    1    2    3    4    5    6    7    8'
  	&'    9   10   11   12   13   14   15   16   17   18   19   20'
	&'   21   22   23'/'-------------------------------------------'
     &'-------------------------------------------------------------'
     &'------------------------')
 
	WRITE(13,18) cyear,cmn,cdy,ctime
	write(13,201) 
	WRITE(14,18) cyear,cmn,cdy,ctime
	WRITE(15,18) cyear,cmn,cdy,ctime
	write(15,202)
  202	format('_____ START___PEAK___END__DURATION')
	WRITE(16,18) cyear,cmn,cdy,ctime
	write(16,201) 

	icalls=0
	iend=0
 
	      DO I=1,34
                   DO K=0,25
	      DA(I,K)=0.
         DEV(I,k)=0.        
	 enddo
	      ENDDO
c
 	do k=0,23
	res(0,k)=0.
	resp(0,k)=0.
	xmed(0,k)=0.
	enddo
	lb1(0,2)=0
		do k=1,6
 	xreg(0,k)=0.
	enddo

C $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
1      DO J=1,12
      R1(J)=0
	R2(J)=0
      ENDDO 
	do j=1,90
	do k=1,5
	lb1(j,k)=0
	enddo
	enddo
	      do J=1,12
	ires(j,25)=0
	  do i=1,10
	  lb2(j,i)=0
	  enddo
      do K=0,24
      AV(J,K)=0.
      SD(J,K)=0.
	XMED(j,k)=0.
	XLEV(j,k)=0.
      ICN(J,K)=0 
	ICD(J,K)=0
	ires(j,k)=0
	res(j,k)=0.
	resp(j,k)=0.
	IDEL(j,k)=0
      enddo
      enddo
c IFL=0
      DO I=1,34
	do j=1,4
	LB3(i,j)=0
	enddo
CDD__       DO K=0,25
CDD__       DEV(I,K)=0.
CDD__	 enddo
      enddo
C 
  10	CONTINUE
C Parameter jjj=1,...7 defines number of days to analyse:
	jjj=1
C Reading foF2 data from styr00.br file:
	do 3 i=1,966
  59   READ (11,12,END=5,ERR=2) JYR,JMN,JDY,DX
   12	FORMAT(3I2,24(2X,F3.0))
c580 IF EOF(2) GOTO 670
c PRINT DA$(34)
	icalls=icalls+1
	JD3(34,1)=JYR
 	 JD3(34,2)=JMN
	  JD3(34,3)=JDY
C Add producing fnF2 reduced by XHI-function from foF2
	   do k=0,23
	ddd=DX(k+1)
	ut=float(k)
	 hour=ut+glon/15.
	if (hour.ge.24.) hour=hour-24.
	DEV(34,k)=ddd	! input foF2 ===> dev

	    DA(34,k)=ddd    ! fnF2 reduced by XHI-function
	   enddo
 	if((JD3(28,1).eq.iy).and.(JD3(28,2).eq.nn).and.(JD3(28,3).eq.id))
	* goto 4
C      if(jjj.gt.7) goto 4
c		if (jmn.eq.0) goto 4
	if(jjj.gt.1) goto 61
	if((jyr.eq.iy).and.(jmn.eq.nn).and.(jdy.eq.id)) then
	goto 61
	else
		goto 62
		endif
  61  jjj=jjj+1
 
c      Compare input data DA(28) equal to given yr,mn,dy:
  62	 do J=1,33
	jd3(j,1)=jd3(j+1,1)
	jd3(j,2)=jd3(j+1,2)
	jd3(j,3)=jd3(j+1,3)
	do k=0,23
	DA(j,k)=DA(j+1,k)
	DEV(j,k)=dev(j+1,k)
	enddo
	enddo
c	
   3  CONTINUE		! end cycle of input
C
   4	CONTINUE
	     if (jjj.lt.7) then
C     Move data up
	mov=7-jjj
	          do m=1,mov
	 do J=1,33
 	jd3(j,1)=jd3(j+1,1)
	jd3(j,2)=jd3(j+1,2)
	jd3(j,3)=jd3(j+1,3)
	do k=0,23
	DA(j,k)=DA(j+1,k)
	DEV(j,k)=dev(j+1,k)
	enddo
	enddo
	          enddo        ! End moving data up
	 	jjj=jjj-1
	endif

C ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      LDNR=34
      LDY=0
      KPN=24*LDNR 
	LDMN=1

C Daily mean, peak and rate-variability of foF2 => results.mpr:	
 71   L=0 
      LH1=0
	LT2=-1
	YY(0)=0.
 72   LDY=LDY+1 
	IF (LDY.gt.LDNR) GOTO 113

      LCNT=0 
	AVE=0. 
	PE=0. 
	DI=0.
  75  L=L+1
C
      DO K=0,23

c REM REPLACE z=z*z for NmF2 later when calculating results.dev, results.ind, results.dis

	Z=DA(LDY,K)
      IF (Z.gt.0.) then
	DEV(L,24)=DEV(L,24)+1.
	LH1=K
	else
	LCNT=LCNT+1
	endif
CDD__      DEV(L,K)=Z      ! NOT!!!!!!
      AVE=AVE+Z
      IF (PE.LT.Z) PE=Z
      YY(K+1)=Z
      ENDDO                       ! End diurnal cycle
      LCNN=24-LCNT 
	IF (LCNN.EQ.0) GOTO 110
 95   IF (YY(0).gt.0) GOTO 97 
      YY(0)=YY(LH1+1) 
      LT2=LH1-24
 97   AVE=AVE/LCNN 
      LT1=LT2 
	Y1=YY(0)
       DO 103 I=1,24
      IF (YY(I).EQ.0.) GOTO 103
       LT2=I-1
       DI=DI+ABS((YY(I)-Y1)/(LT2-LT1))
      LT1=LT2 
	Y1=YY(I)
  103  CONTINUE                   
      YY(0)=YY(LH1+1)
	LT2=LH1-24
      DD=DI/LCNN
c REM PRINT DY,AVE;PE;DD;CNT
C For NmF2: AVE=INT(AVE/100+.5) : PE=INT(PE/100+.5) : DI=INT(DD/10+.5)
C For foF2:
C      AVE=INT(AVE+.5) : PE=INT(PE+.5) : DI=INT(DD*10+.5)
 	DI=DD*10.
  110 LB3(L,1)=int(AVE +0.5)
      LB3(L,2)=int(PE+0.5)
	LB3(L,3)=int(DI+0.5)
	LB3(L,4)=LCNT
c PRINT LDY,AV;PE;DI;LCNT
  120 GOTO 72
  113 	CONTINUE

C ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c  MONTHLY MEAN & SD
C  VERIFICATION OF MEDIAN
  115  LDNR=27
c      DO 150 J=1,7
      DO 150 J=1,jjj
        DO 144 M=0,23          ! Hour-to-hour cycle: variable <m>
        MM=0
          DO 124 N=1,LDNR
        NJ=N+J-1
CDD__       X1(N)=DEV(NJ,M)
	X1(N)=DA(NJ,M)
CDD__       AV(J,M)=AV(J,M)+DEV(NJ,M)
	AV(J,M)=AV(J,M)+DA(NJ,M)
CDD__       IF (DEV(NJ,M).GT.0.) ICN(J,M)=ICN(J,M)+1
	 IF (DA(NJ,M).GT.0.) ICN(J,M)=ICN(J,M)+1
 124   CONTINUE
      ICN(J,24)=ICN(J,24)+ICN(J,M)
      IF (ICN(J,M).GT.0) AV(J,M)=AV(J,M)/ICN(J,M)
C  prepare array x2:
 	do nn=1,ldnr
	x2(nn)=0.
	enddo

 127  XX=0.
C  : MAX X1(N)
 128  DO 130 N=1,LDNR
      IF (XX.LT.X1(N)) XX=X1(N)
 130  CONTINUE

      IF (XX.EQ.0.) GOTO 137
      DO 135 N=1,LDNR
      IF (X1(N).LT.XX) GOTO 135
      MM=MM+1
	X2(MM)=XX 
	X1(N)=0.
 135  CONTINUE
      IF (MM.LT.LDNR) GOTO 127
 137  XMM=float(MM)
      XD1=XMM/2.
	ND1=INT(XD1)
	CD1=float(ND1) 
      ND2=ND1+1 
c//	IF ((XD1-CD1).LT.0.0001) THEN 
      IF (CD1.LT.XD1) THEN 
	XMD=X2(ND2)
	else
	XMD=(X2(ND1)+X2(ND2))/2.
	endif
      IF (MM.EQ.0) XMD=xmed(j-1,m)
      XMED(J,M)=XMD
C
 144  CONTINUE        ! cycle M
C 
C Include Daily Regression of foF2 on fcmed (calibrated median for 27 preceding days):
	ik=0
C Add count of DEV- and DEV+
	ikn=0
	ikp=0

	sx=0.
	sy=0.
	sx2=0.
	sxy=0.

	do 146 k=0,23
CDD__	zz=dev(j+27,k)
	zz=da(j+27,k)
	if (zz.eq.0.) goto 146
	ik=ik+1 
      ry(ik)=zz
	rx(ik)=xmed(j,k)
	sx=sx+rx(ik)
	sx2=sx2+rx(ik)*rx(ik)
	sy=sy+ry(ik)
	sxy=sxy+rx(ik)*ry(ik)
  146	continue                ! cycle K
	if (ik.lt.3) then
	xreg(j,1)=xreg(j-1,1)
	xreg(j,2)=xreg(j-1,2)
	xreg(j,3)=xreg(j-1,3)
	sx=xreg(j,1)
	sy=xreg(j,2)
	ra=xreg(j,3)
	goto 145
	             else
	sx=sx/ik
	sx2=sx2/ik
	sy=sy/ik
	sxy=sxy/ik
	 RA=(SXY-SX*SY)/(SX2-SX*SX)
	xreg(j,1)=sx
	xreg(j,2)=sy
	xreg(j,3)=RA
	            endif
  145	continue
C Include Daily Standard Deviation (DSD) from daily regression line:
C Additional : DSD- and DSD+
	 xreg(j,4)=0.	  ! DSD
	xreg(j,5)=0.	  ! DSD-
	xreg(j,6)=0.	  ! DSD+
	do i=0,23
c	if (ik.lt.3) then
c	xlev(j,i)=xmed(j,i)
c	 else
	xlev(j,i)=sy+ra*(xmed(j,i)-sx)
c	endif
CDD__	     if (dev(j+27,i).eq.0.) goto 149
	     if (da(j+27,i).eq.0.) goto 149
CDD__	xreg(j,4)=xreg(j,4)+(dev(j+27,i)-xlev(j,i))**2.	      ! DSD
	xreg(j,4)=xreg(j,4)+(da(j+27,i)-xlev(j,i))**2.	      ! DSD
CDD__	if (dev(j+27,i).lt.xlev(j,i)) then
	if (da(j+27,i).lt.xlev(j,i)) then
	 ikn=ikn+1
CDD__	 xreg(j,5)=xreg(j,5)+(dev(j+27,i)-xlev(j,i))**2.	  ! DSD-
	 xreg(j,5)=xreg(j,5)+(da(j+27,i)-xlev(j,i))**2.	  ! DSD-
	endif
CDD__	if (dev(j+27,i).gt.xlev(j,i)) then
	if (da(j+27,i).gt.xlev(j,i)) then
	ikp=ikp+1
CDD__	 xreg(j,6)=xreg(j,6)+(dev(j+27,i)-xlev(j,i))**2.	  ! DSD+
	 xreg(j,6)=xreg(j,6)+(da(j+27,i)-xlev(j,i))**2.	  ! DSD+
	endif
  149	continue
	enddo           ! cycle i
  	if (ik.gt.0) xreg(j,4)=sqrt(xreg(j,4)/ik)			  ! DSD
	if (ikn.gt.0) xreg(j,5)=-sqrt(xreg(j,5)/ikn)		  ! DSD-
	if (ikp.gt.0) xreg(j,6)=sqrt(xreg(j,6)/ikp)			  ! DSD+
C		 
 150	continue         ! cycle j
C !!!!!!! Not Use XLEV instead of XMED as quiet conditions foF2 !!!!!!!!!!!!!!!!!!
C
C Median results:
C      DO MJ=1,7
      DO MJ=1,jjj
C Limit output by results for given month:
	if (JD3(MJ+27,2).ne.mon) then 
	iend=1
      goto 1650
	endif

	do i=0,23
C Restore foF2med from proxy fnF2 reduced by XHI-function:
	ddd=xmed(mj,i)
	ut=float(i)
	 hour=ut+glon/15.
	if (hour.ge.24.) hour=hour-24.
C Old: med for 28th day:
	ndyout=ndoy(JD3(MJ+27,1),JD3(MJ+27,2),JD3(MJ+27,3))
C NEW!!!! Recalculate	med for 14th day before (centered):
	ndyout=ndyout-14
	if (ndyout.le.0) ndyout=365+ndyout
	xmed(mj,i)=ddd  ! Back from xhi-reduction but day=(28-14) !!!!
C	imed(i)=nint(xmed(mj,i)) ! 27-days median => output
		imed(i)=nint(ddd)
c	imed(i)=nint(xlev(mj,i))
	enddo
	nd=JD3(MJ+27,3)
	 n10=nd/10
      n1=nd-n10*10
	do i=0,9
	 if (n10.eq.i) then
	outdy(1:1)=IN(i)
	endif
	if (n1.eq.i) then
 	outdy(2:2)=IN(i)
	endif
	enddo
	if (outdy0.eq.outdy) goto 1650
	WRITE(13,158) outyr,outmn,outdy
c	*,(IMED(K),K=0,23),(xreg(mj,i),i=1,4)
	*,(IMED(K),K=0,23)
c	WRITE(*,158) JD3(MJ+27,1),JD3(MJ+27,2),JD3(MJ+27,3)
CC	WRITE(*,158) outyr,outmn,outdy
c	*,(IMED(K),K=0,23),(xreg(mj,i),i=1,4)
cdd	*,(IMED(K),K=0,23)
c 158  FORMAT(3I2,24(2X,I3),4(1X,F8.4))
 158  FORMAT(3A2,24(2X,I3))
	outdy0=outdy
	ENDDO
C ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 1650 continue
	outdy0='00'
C MPR- results:
CD     	WRITE(14,166) JD3(28,1),JD3(28,2),JD3(28,3)
CD	   	WRITE(*,166) JD3(28,1),JD3(28,2),JD3(28,3)
CD 166  FORMAT(1X,3I2)
c      DO MJ=1,7
      DO MJ=1,jjj
C Limit output by results for given month:
	if (JD3(MJ+27,2).ne.mon) goto 173
	 	imj0=JD3(mj+26,2)
  	imj1=JD3(mj+27,2)
	if (imj1.eq.imj0) goto 168
c//	if (imj1.gt.imj0) then
	if (jd3(mj+27,3).eq.1) then
	write (14,167) iyyyy,ti(imj1)
CC	write (*,167) iyyyy,ti(imj1)
  167	format(1X,I4,1X,A4)
	endif
 168  WRITE(14,169) jd3(mj+27,3),LB3(MJ+27,1),LB3(MJ+27,2),LB3(MJ+27,3)
     +,LB3(MJ+27,4),xreg(mj,4),xreg(mj,5),xreg(mj,6)
CC	WRITE(*,169) jd3(mj+27,3),LB3(MJ+27,1),LB3(MJ+27,2),LB3(MJ+27,3)
Cdd	+,LB3(MJ+27,4),xreg(mj,4),xreg(mj,5),xreg(mj,6)
 169   FORMAT(1X,I2,3(I4),1X,I2,3(1X,F6.2))
      ENDDO

C ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C ========= Percentage deviations for hourly NmF2 ================
 173  L=27
      LC=27
	MN=0      ! Prepare line of array [1:12] for 7 days analysed
c	LDNR=7
	LDNR=jjj
 174  LDY=0 
      LC=LC+LDNR
C-----------------------------------
C Cycle for jjj=1...7 days:
C-----------------------------------
 175  CONTINUE

C  Limit calculations to dev% output
 179     LDY=LDY+1 
       MN=MN+1 
	 IF (LDY.GT.LDNR) GOTO 301
 180  L=L+1
      DO 191 K=0,23
      DD=999.
	delog=999.
	DE=0.
      IF (DEV(L,K).EQ.0.) GOTO 190
C REM DE=DEV(L,K)-AV(MN,K)
C      DE=DEV(L,K)-XLEV(MN,K)
C ========== Put: foF2*foF2 for NmF2==================:
C--	XME=XLEV(MN,K)*XLEV(MN,K)
	XME=XMED(MN,K)*XMED(MN,K)
	DE=DEV(L,K)*DEV(L,K)-XME  ! dev=fob**2-fmed**2

C ***TWO OPTIONS: AV(MN,K) OR MED(MN,K)***Nm ~ ZZ !!!
C IF AV(MN,K)>0 THEN DD=INT(DE/AV(MN,K)*100+.5)
C--      IF (XLEV(MN,K).GT.0.) then
      IF (XMED(MN,K).GT.0.) then
C	DD=DE/XMED(MN,K)*100.	  	 
C ========== Put: foF2*foF2 for NmF2==================:
		DD=DE/XME*100.		  ! dev=(fob**2-fmed**2)/fmed**2*100%
C--	delog=2.*(alog10(dev(l,k))-alog10(XLEV(mn,k)))
	delog=2.*(alog10(dev(l,k))-alog10(XMED(mn,k)))
C-NO!	DD=delog                  ! dev from ratio of 2*(log10(fob/fmed))
	endif
 190     SD(MN,K)=SD(MN,K)+DE*DE
c	res(mn,k)=dd
	res(mn,k)=delog		   ! Log(Ne/Nm)
c	resp(mn,k)=de/xme*100. 
	resp(mn,k)=dd		   ! Percentage deviations of NmF2 from median
       if (DD.GT.0.) then
C	 IRES(MN,K)=nint(DD+0.5) ! ID-index
	 IRES(MN,K)=nint(DELOG*1000+0.5)  ! delog-index
	if (ires(mn,k).gt.9999) ires(mn,k)=9999
	                else
c	IRES(MN,K)=nint(DD-0.5)  ! ID-index
C	IRES(MN,K)=nint(DD*1000-0.5)  ! delog-index
	 IRES(MN,K)=nint(DELOG*1000-0.5)  ! delog-index
	if (ires(mn,k).le.-1000) ires(mn,k)=-999
	endif
	if (delog.eq.999.) then 
	idel(mn,k)=99			  ! ID-index
	ires(mn,k)=999
C	idel(mn,k)=999			  ! delog
	goto 191
	           endif
	if (delog.eq.0.) then
	IDEL(mn,k)=0			  ! ID-index=delog
	goto 191
	endif
      if(delog.lt.tr(1)) then
	idel(mn,k)=-4                 ! ID-index
C	idel(mn,k)=nint(delog*100+.5)  ! delog
	goto 191
	endif
	if (delog.gt.tr(7)) then
	idel(mn,k)=4				  ! ID-index
C	idel(mn,k)=nint(delog*100+.5)  ! delog
 	goto 191
	endif
	if (delog.eq.0.) then
	 idel(mn,k)=0				  ! ID-index
C	idel(mn,k)=int(delog*100+.5)  ! delog
  	goto 191
	endif
	do 189 n=1,6
	   if (delog.lt.0.) then
	        if ((delog.ge.tr(n)).and.(delog.lt.tr(n+1))) then
	        idel(mn,k)=iin(n)              ! ID-index
	         endif
	                 else
	        if ((delog.gt.tr(n)).and.(delog.le.tr(n+1))) then
 	        idel(mn,k)=iin(n)              ! ID-index
	        endif
	                    endif
  189	continue

 191  CONTINUE
      GOTO 175
C-----------------------------------
 301  CONTINUE
C  Deviations results
       DO I=1,LDNR
C Limit output by results for given month:
	if (JD3(i+27,2).ne.mon) goto 1324
	nd=JD3(I+27,3)
	  n10=nd/10
      n1=nd-n10*10
	do k=0,9
	 if (n10.eq.k) then
	outdy(1:1)=IN(k)
	endif
	if (n1.eq.k) then
 	outdy(2:2)=IN(k)
	endif
	enddo
	if (outdy0.eq.outdy) goto 1324
	WRITE(12,323) outyr,outmn,outdy
	+,(IRES(I,K),K=0,23)
CC	WRITE(*,323) outyr,outmn,outdy
Cdd	+,(IRES(I,K),K=0,23)
 323  FORMAT(3A2,24(1X,I4))
	outdy0=outdy
	ENDDO

C ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C Include selection of D+,D- intervals duration > 3 hrs
C
C Change criteria : average sum(-indices)<=-3 //
C Change criteria : average sum(+indices)>=+3 //
C
 1324 L=0 
	outdy0='00'
	LC=0
	LDY=0 
	  LC=LC+LDNR 
920	LDY=LDY+1 
	  IF (LDY.gt.LDNR) GOTO 1540
	L=L+1 
	KMI=0 
	KMA=0 
	idin=0
	indin=0
	idax=0
	indax=0
	LT1=0
	LT2=0
	LT3=0
	LT4=0
	LDT12=0
	LDT34=0
C
	 do 1060 k=0,23
c	idd=ires(l,k)   ! Log(Ne/Nm)
      pdd=resp(l,k)	! = percentage deviation
	idd=nint(pdd)	! dev%
	indd=idel(l,k)	 ! Indices
	IF (K.gt.0) THEN
c	 id1=ires(L,K-1)
C ///	 id1=nint(resp(L,K-1)) 
	id1=idel(L,K-1)	 ! Indices
	 ELSE 
c	 id1=ires(L-1,23)
C ///	 id1=nint(resp(L-1,23))
	 id1=idel(L-1,23)  ! Indices
	endif
	IF (K.LT.23) THEN
c	 id2=ires(L,K+1)
C ///	 id2=nint(resp(L,K+1)) 
	id2=idel(L,K+1)     ! Indices
	              ELSE
		IF (L.lt.LDNR) then
c	    id2=ires(L+1,0)
C ///	    id2=nint(resp(L+1,0))
	id2=idel(L+1,0)     ! Indices
	    else
	    id2=999
	    endif
	ENDIF
 1000	IF (idd.eq.999) GOTO 1060
	 IF((id1.eq.99).AND.(id2.eq.99)) GOTO 1060
      IF (idd.lt.idin) then 
	idin=idd
	rdin=pdd
	indin=indd
	endif
	IF (idd.gt.idax) then 
	idax=idd
	rdax=pdd
	indax=indd
	endif
      IF (idin.eq.idd) KMI=K
      IF (idax.eq.idd) KMA=K
 1060	continue
C	PRINT L,"din=";DIN,"dax=";DAX

C !!!!!!!!!!!! Different thresholds for D- and D+ !!!!!!!!!!!!!!!!!!!!!!
C	DURATION FOR -40% =<(DEV%NmF2)<= +40% is changed to:
C				 -30% =<(DEV%NmF2med-reg)<= +43%		 !
C                   -0.155=<alog10(devNm/Nmed-reg)<=0.155
C                  -16% < dev%fof2/fomed-reg < 20%       
C ///	LTHN=-30
		LTHN=-3	  !	Indices threshold

C	LTHN=-20	 ! According to index=-4 for delog<-.2
C Ole:	LTHN=-40
C ///	LTHP=43
	LTHP=3        !	Indices	threshold

C Ole:	LTHP=40
C      LTHP=20		 ! According to index=+4 for delog>+.2
	LTM=jd3(L+27,2)*10000+jd3(L+27,3)*100+KMI
	 LHM=jd3(L+27,2)*10000+jd3(L+27,3)*100+KMA
C Negative event:
	IF (indin.gt.LTHN) GOTO 1280
	K1=KMI  
	J1=L  
	ISM=indin 
	 KK=1
1140	K2=K1  
	J2=J1
	IF ((J1.eq.0).AND.(K1.eq.0)) GOTO 1200
	K1=K1-1  
	IF (K1.ge.0) GOTO 1180
	J1=J1-1 
	K1=K1+24
c1180	IF (IRES(J1,K1)*IRES(J2,K2).le.0) GOTO 1200
c \\\1180	IF (RESP(J1,K1)*RESP(J2,K2).le.0.) GOTO 1200
1180	IF (IDEL(J1,K1)*IDEL(J2,K2).le.0.) GOTO 1200
	KK=KK+1  
c	ISM=ISM+IRES(J1,K1)
C \\\ ISM=ISM+nint(RESP(J1,K1)) 
	ISM=ISM+idel(J1,K1)       ! Indices
	 IF (ISM/KK.le.LTHN) then
	 GOTO 1140
	else
c		ISM=ISM-IRES(J1,K1)
C \\\		ISM=ISM-nint(RESP(J1,K1))
 	 		ISM=ISM-idel(J1,K1)     ! Indices
 	 KK=KK-1
	endif
1200	LT1=jd3(J2+27,2)*10000+jd3(j2+27,3)*100+K2  
	LD1=J2  				 
	LTD1=K2
	K2=KMI  
	J2=L  
c	ISM=ISM-IRES(J1,K1) 
c	 KK=KK-1
1220	K1=K2  
	J1=J2
	K2=K2+1  
	IF (K2.lt.24) GOTO 1250
	J2=J2+1 
	K2=K2-24
c1250	IF (IRES(J1,K1)*IRES(J2,K2).le.0) GOTO 1270
C \\\1250	IF (RESP(J1,K1)*RESP(J2,K2).le.0.) GOTO 1270
 1250	IF (idel(J1,K1)*idel(J2,K2).le.0.) GOTO 1270      ! Indices
	KK=KK+1  
c	ISM=ISM+IRES(J2,K2) 
C \\\	ISM=ISM+nint(RESP(J2,K2))  
		ISM=ISM+idel(J2,K2)    ! Indices
	IF (ISM/KK.le.LTHN) GOTO 1220
1270	LT2=jd3(J1+27,2)*10000+jd3(j1+27,3)*100+K1
	LD2=J1  
	LTD2=K1
c Positive event:
1280	IF (INDAX.lt.LTHP) GOTO 1440
	K1=KMA  
	J1=L  
	ISM=INDAX 
	 KK=1
1300	K2=K1 
	 J2=J1
	IF ((J1.eq.0).AND.(K1.eq.0)) GOTO 1360
	K1=K1-1 
	 IF (K1.ge.0) GOTO 1340
	J1=J1-1 
	K1=K1+24
c1340	IF (IRES(J1,K1)*IRES(J2,K2).le.0) GOTO 1360
C \\\1340	IF (RESP(J1,K1)*RESP(J2,K2).le.0.) GOTO 1360
1340	IF (idel(J1,K1)*idel(J2,K2).le.0.) GOTO 1360
	KK=KK+1 
c	 ISM=ISM+IRES(J1,K1)
C \\\	 ISM=ISM+nint(RESP(J1,K1)) 
		 ISM=ISM+idel(J1,K1)      ! Indices
C	 IF ((ISM/KK.ge.LTHP).AND.(IRES(J1,K1).lt.990)) then
	 IF ((ISM/KK.ge.LTHP).AND.(RESP(J1,K1).lt.999.)) then
	 GOTO 1300
	else
c		  ISM=ISM-IRES(J1,K1) 
C \\\		  ISM=ISM-nint(RESP(J1,K1)) 
			  ISM=ISM-idel(J1,K1)      ! Indices
	   KK=KK-1
	   endif
1360  LT3=jd3(J2+27,2)*10000+jd3(j2+27,3)*100+K2  
       LD3=J2 
	  LTD3=K2
	  K2=KMA 
	  J2=L 
c	  ISM=ISM-IRES(J1,K1) 
c	   KK=KK-1
1380	K1=K2 
	 J1=J2
	K2=K2+1 
	 IF (K2.lt.24) GOTO 1410
	J2=J2+1 
	K2=K2-24
c1410 	  IF (IRES(J1,K1)*IRES(J2,K2).le.0) GOTO 1430
C \\\1410 	  IF (RESP(J1,K1)*RESP(J2,K2).le.0.) GOTO 1430
 1410 	  IF (idel(J1,K1)*idel(J2,K2).le.0.) GOTO 1430
	KK=KK+1 
C	 ISM=ISM+IRES(J2,K2) 
C \\\	 ISM=ISM+nint(RESP(J2,K2)) 
		 ISM=ISM+idel(J2,K2)	 ! Indices
c	  IF ((ISM/KK.ge.LTHP).AND.(IRES(J2,K2).lt.990)) GOTO 1380
	  IF ((ISM/KK.ge.LTHP).AND.(RESP(J2,K2).lt.999.)) GOTO 1380
1430  LT4=jd3(J1+27,2)*10000+jd3(j1+27,3)*100+K1
      LD4=J1 
	 LTD4=K1
1440	IF (INDIN.gt.LTHN) GOTO 1460
	 LDT12=(LD2-LD1)*24+LTD2-LTD1+1
1460	IF (INDAX.lt.LTHP) GOTO 1490
	LDT34=(LD4-LD3)*24+LTD4-LTD3+1
      
C	pause ' '
1490	CONTINUE
CC	write(*,*) INDIN,rdin,LT1,LTM,LT2,LDT12,L
CC	write(*,*) INDAX,rdax,LT3,LHM,LT4,LDT34
		II=L
c	LB2(L,1)=IDIN    		   ! I4
 	LB2(L,1)=nint(RDIN)		   ! I4  min D-
	LB2(L,2)=LT1			   ! I4	 hr1-
	LB2(L,3)=LTM 			   ! I4	 hrmin
	LB2(L,4)=LT2 			   ! I4	 hr2-
	LB2(L,5)=LDT12			   ! Duration D-
c      LB2(L,6)=IDAX   		   ! I4
      LB2(L,6)=nint(RDAX)   		   ! I4	 max D+
	LB2(L,7)=LT3               ! I4	 hr1+
	 LB2(L,8)=LHM 			   ! I4	 hrmax
	LB2(L,9)=LT4			   ! I4	 hr2+
	LB2(L,10)=LDT34			   ! I3	 Duration D+

	GOTO 920
C
1540	continue
C
1670  continue

C 

	jj=0
	do 1970 j=1,ii
	if (lb2(j,5).lt.3) goto 1920
c	if (jj.eq.0) goto 1919
c	if (lb2(jj,3).eq.lb2(jj-1,3)) goto 1920
1919		jj=jj+1
	do k=1,5
	LB1(JJ,K)=LB2(J,K)
	enddo
 1920	IF (LB2(J,10).lt.3) GOTO 1970
c		if (jj.eq.0) goto 1918
c	if (lb2(jj,8).eq.lb2(jj-1,8)) goto 1970
1918	jj=jj+1
	do k=6,10
	LB1(JJ,K-5)=LB2(J,K)
	enddo
 1970	continue
 	 ii=jj
1980	if (ii.eq.0) goto 1700
      if (ii.eq.1) goto 2250
      j=1
 1990	IF (LB1(J+1,2).eq.0) GOTO 2140
       IF ((j.gt.1).and.(LB1(j-1,2).eq.LB1(j-2,2))) then
	j=j-2
	GOTO 2140
	endif

c	PRINT LB1(J,2),LB1(J+1,2)
	IF (lb1(j+1,2).gt.lb1(j,4)) goto 2100
	IF((LB1(J+1,2).ge.LB1(J,2)).AND.(LB1(J+1,4).le.LB1(J,4))) 
	+ GOTO 2140

C Add control lb1(j+1,2).le.lb1(j,4) :

	if (lb1(j+1,2).le.lb1(j,4)) then 
	if (lb1(j+1,1)*lb1(j,1).lt.0) then
	 if (lb1(j+1,2).lt.lb1(j,2)) goto 2090
	else
	 goto 2140
	endif
	if (lb1(j+1,2).eq.lb1(j,4)) then
	lb1(j,5)=lb1(j,5)+lb1(j+1,5)-1
	lb1(j,4)=lb1(j+1,4)
	goto 2140
	                             else
	if (lb1(j,5).gt.lb1(j+1,5)) then
	goto 2140
	                            else
	goto 2100
		endif
	endif
	endif

	IF (LB1(J,2).lt.LB1(J+1,2)) goto 2100
 2090 	do k=1,5
	LC(K)=LB1(J+1,K)
	LB1(J+1,K)=LB1(J,K)
	LB1(J,K)=LC(K)
	enddo
C	PRINT J,LB1(J,1),LB1(J,2),LB1(J,3),LB1(J,4),LB1(J,5)
 2100	j=j+1
C
2110	IF (J.lt.II) GOTO 1990
	goto 2250
 2140	ii=ii-1
	i1=j+1
		IF (LB1(J,5).le.LB1(J+1,5)) THEN
 	 I1=J 
	 ELSE 
	 I1=J+1
	endif

	IF((LB1(J+1,2).eq.LB1(J,2)).AND.(LB1(J+1,4).eq.LB1(J,4))) then
	if (abs(lb1(j,1)).lt.abs(lb1(j+1,1))) then 
 		do k=1,5
	LB1(J,K)=LB1(J+1,K)
	enddo
	endif
						j=j+1
      i1=j
	endif
C
2150	DO i=i1,ii
	do k=1,5
	LB1(I,K)=LB1(I+1,K)
	enddo
	ENDDO
 	goto 1980

 2250	continue

C  Results for disturbed periods:
 	DO 290 J=1,II
	imon0=lb1(j-1,2)/10000
	imon1=lb1(j,2)/10000
C Limit output by results for given month:
	if (imon1.ne.mon) goto 1700

	if (imon1.gt.imon0) then
	write (15,270) iyyyy,ti(imon1)
CC	write (*,270) iyyyy,ti(imon1)
	endif
C--	ltemp=LB1(j,1)
C--	LB1(j,1)=LB1(J,1)/10
  270	format(1X,I4,1X,A4)
	write(15,280) (LB1(j,k),k=1,5)
CC	write(*,280) (LB1(j,k),k=1,5)
  280	Format(1X,I4,3(1X,I6),1X,I3)
C--	LB1(j,1)=ltemp
  290		continue
C ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 1700	continue
C
C  Logarithmic indices results
C 
	       DO I=1,LDNR
C Limit output by results for given month:
	if (JD3(i+27,2).ne.mon) goto 1325
	nd=JD3(i+27,3)
	  n10=nd/10
      n1=nd-n10*10
	do k=0,9
	 if (n10.eq.k) then
	outdy(1:1)=IN(k)
	endif
	if (n1.eq.k) then
 	outdy(2:2)=IN(k)
	endif
	enddo
	if (outdy0.eq.outdy) goto 1325
 	WRITE(16,324) outyr,outmn,outdy
	+,(IDEL(I,K),K=0,23)
c	WRITE(*,324) JD3(I+27,1),JD3(I+27,2),JD3(I+27,3)
c-		WRITE(*,324) outyr,outmn,outdy
c-	+,(IDEL(I,K),K=0,23)
C 324  FORMAT(1X,3I2,24(1X,I2))         ! ID-index
 324  FORMAT(3A2,24(1X,I4))         !
  	outdy0=outdy
	ENDDO

1325	 outdy0='00'
       iy=jd3(34,1)
	nn=jd3(34,2)
	id=jd3(34,3)+1
	if (id.gt.im(nn)) then
	id=id-im(nn)
	nn=nn+1
	if (nn.gt.12) then
	nn=1
	iy=iy+1
	endif
	 endif
C Move former data for one day forware:
  	 do J=1,33
	jd3(j,1)=jd3(j+1,1)
	jd3(j,2)=jd3(j+1,2)
	jd3(j,3)=jd3(j+1,3)
	  do k=0,23
	  DA(j,k)=DA(j+1,k)
	DEV(j,k)=DEV(j+1,k)
	  enddo
	enddo
	do i=1,10         ! move lb2 up 
	lb2(0,i)=lb2(jjj,i)
	 enddo
c
	do k=0,23
c	res(0,k)=res(7,k)
c	xmed(0,k)=xmed(7,k)
		res(0,k)=res(jjj,k)
			ires(0,k)=ires(jjj,k)
	resp(0,k)=resp(jjj,k)
	xmed(0,k)=xmed(jjj,k)
	enddo
	do k=1,6
c	xreg(0,k)=xreg(7,k)
		xreg(0,k)=xreg(jjj,k)
 	enddo
C Remember last month/day START:
	 lb1(0,2)=lb1(ii,2)
CTemporary:
C	goto 6
C
	if (iend.eq.0) then
      GOTO 1
	                else
	goto 6
	endif
c
    2 CONTINUE
      if (icalls.eq.0) then
	       write(*,*) 'INPUT FILE IS NOT IN YOUR DIRECTORY '
	else
	write(*,*) 'Error at icalls=',icalls
	   	pause ' '
      endif
    5	iend=1
	goto 4
    6 close(11)
     	close(12)
	close(13)
	close(14)
	close(15)
	close(16)
c-	 if((iy.gt.0).and.(iend.eq.1)) goto 1100 
c==	if (iend.eq.2) goto 102
	
c      if ((iend.eq.2).and.) goto 102

c==	goto 1100 
 102	continue
cc	close(unit=108)
	ii=iic    ! restore number of stations
	RETURN              
	END
