C /////////////////////////////////////////////////////////////// SUBROUTINE SUBDEVGXc(ist) C...........................................................Dec 2023 C C ..........................................................June 2009 C Enter station number 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 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 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\graf\temp\' INFILE(18:19)=sta INFILE(20:21)=yr INFILE(22:)='f.br' c infile='c:\web\graf\temp\stYRf.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 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 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=+.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 forward: 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