subroutine subwp3n(AYR,AMN,IDY1,IDY2) C.....................................................Nov. 2012 C JPL results C.....................................................Oct.2012 C Update monthly Wp, Np, Pp files C C.....................................................Apr. 2012 C C IONEX maps C for glat=-87.5:2.5:87.5, glong=-180:5:175 (5112) C Correction of (W+max-W-max)*(1+cntstrm/5112) C.....................................................Feb. 2012 C C.....................................................Aug.2011 C C........................................................Apr.2008 C...... Count RANGE only when W+=>3, and/or W-<=-3 .............. C c........................................................Oct.2007 C cntstrm=total number of W+=3,W+=4, W-=-3,W-=-4 for the IONEX-map C c .......................................................Sep. 2007 C Range of TEC indices from 35 latitude results TEC.fYR......... C c from hourly IONEX maps c DIMENSION RESP(0:71),RESN(0:71),DHR(0:71),IM(12) +,cnts(71,3,0:23),xind(71,3,0:23) +,RES(0:71),RUP(0:71),RLO(0:71) +,clat(71,3,0:23),relat(71,3,0:23) C++ INTEGER*2 IYR,MN,IDY,idy1,idy2,iut CHARACTER(10) DD CHARACTER(10) TT CHARACTER(5) ZZ CHARACTER*66,infile,outfile,outfilep,outfilen !Ljuba поменяла 64 на 66 character*41 inf,ou !Ljuba character*41 ile,tfile !Ljuba CHARACTER*64 outinfo,title CHARACTER*133 wpast(0:365),wwpos(0:365),wwneg(0:365),txt CHARACTER*2 AYR,AMN,AUT,ADY +,aa(24),bb CHARACTER*6 YRMNDY CHARACTER*4 AYEAR,AMEAN DATA IM/31,28,31,30,31,30,31,31,30,31,30,31/ C Controlling Selected Latitudes 1=Yes, 0 = No: DATA aa/'00','01','02','03','04','05','06','07','08','09','10' +,'11','12','13','14','15','16','17','18','19','20','21','22','23'/ bb='UT' amean='mean' title='YYMMDD UT CNTNEG CNTPOS CNTTOT TOT% Wp' C CALL DATE_AND_TIME(DATE=DD,TIME=TT,ZONE=ZZ) !Tamara CALL DATE_AND_TIME(DD,TT,ZZ) !Ljuba TT(5:)=' ' WRITE(*,*) 'PC Date: Year,Month,Day = ',DD,'Time = ',TT C start AYEAR(1:2)='20' AYEAR(3:4)=ayr read(ayr,*) year iyr=int(year) read(amn,*) rmn mn=int(rmn) if(int(year/4.)*4.eq.iyr) THEN ndnr=366 im(2)=29 ELSE ndnr=365 im(2)=28 ENDIF 13 FORMAT(I2) CC CALCULATION OF DAYS OF YEAR C infile='d:\web\graf\dwc\YR\MN\ut\wcMNDYUT.jYR' !W(glat & glong map)!Tamara C infile(17:18)=ayr !Tamara C infile(20:21)=amn !Tamara C infile(36:37)=ayr !Tamara C infile(28:29)=amn !Tamara inf='/var/www/izmiran/ionosphere/weather/graf/'!W(glat & glong map)!Luba41 ile='dwc/YR/MN/ut/wcMNDYUT.jYR' !Ljuba 25 ile(5:6)=ayr !Ljuba YR ile(8:9)=amn !Ljuba MN ile(24:25)=ayr !Ljuba YR посл ile(16:17)=amn !Ljuba MN infile=inf//ile C outfile='d:\web\graf\datw\wptecYR.txt' !W-index from W-maps geogr.coord.!Tamara C outfile='d:\web\graf\datw\wptecYR.txt' ! !Tamara C outfile(23:24)=AYR !Tamara ou='/var/www/izmiran/ionosphere/weather/graf/' !Ljuba 41 tfile='datw/wptecYR.txt' !Ljuba tfile(11:12)=AYR !Ljuba YR outfile=ou//tfile !Ljuba outfilep=outfile outfilen=outfile C outfilep(18:18)='p' !Tamara C outfilen(18:18)='n' !Tamara outfilep(47:47)='p' !Ljuba wptecYR.txt 1 символ outfilen(47:47)='n' !Ljuba wptecYR.txt 1 символ C outinfo='wjinfoYR.txt' !Tamara C outinfo(7:8)=ayr !Tamara outinfo='/var/www/izmiran/ionosphere/weather/graf/wjinfoYR.txt' !Ljuba !!! outinfo(48:49)=ayr !Ljuba YR !!!!!!!!!!! OPEN(10,FILE=OUTINFO,ACCESS='APPEND') write(10,14) title 14 format(1X,A43) C OPEN(12,FILE=OUTFILE) OPEN(13,FILE=OUTFILEP) OPEN(14,FILE=OUTFILEN) 198 format(4X,A2,24(3X,A2),3X,A4) C lda1=ndoy(iyr,mn,idy1) ! day-of-year mndy=im(mn) 130 format(A133) if((mn.eq.1).and.(idy1.eq.1)) then write(12,198) bb,aa,amean write(13,198) bb,aa,amean write(14,198) bb,aa,amean write(*,198) bb,aa,amean goto 197 ! avoid read-re-write endif C C read-re-write outfile for days from 1 to lda1-1: ncnt=0 do 149 ipast=0,lda1-1 read(12,130,END=151,ERR=2) txt wpast(ipast)=txt YRMNDY=txt(1:6) 151 read(13,130,END=152,ERR=2) txt wwpos(ipast)=txt 152 read(14,130,END=153,ERR=2) txt wwneg(ipast)=txt if (yrmndy.eq.' ') goto 149 ncnt=ncnt+1 149 continue 153 close(unit=12) close(unit=13) close(unit=14) C OPEN(12,FILE=OUTFILE) OPEN(13,FILE=OUTFILEP) OPEN(14,FILE=OUTFILEN) C re-write past data: do ipast=0,ncnt-1 write(12,130) wpast(ipast) cc write(*,130) wpast(ipast) ! test write(13,130) wwpos(ipast) write(14,130) wwneg(ipast) enddo close(unit=12) close(unit=13) close(unit=14) CC 197 continue OPEN(12,FILE=OUTFILE,ACCESS='APPEND') OPEN(13,FILE=OUTFILEP,ACCESS='APPEND') OPEN(14,FILE=OUTFILEN,ACCESS='APPEND') C C day-to-day cycle C ======================================== DO 500 idy=idy1,idy2 lda=ndoy(iyr,mn,idy) ! day-of-year call blet2(idy,ady) C infile(30:31)=ady !Tamara infile(59:60)=ady !Ljuba DY C C cnts(71,3,0:23),xind(71,3,0:23) do i=0,23 ! ut=0,1,...,23 do k=1,3 do n=1,71 cnts(n,k,i)=0. xind(n,k,i)=0. clat(n,k,i)=0. !clat(71,3,0:23) relat(n,k,i)=0. enddo enddo enddo C UT do k=0,23 res(k)=0. rup(k)=0. rlo(k)=0. resp(k)=0. resn(k)=0. enddo C sum=0. sumn=0. sump=0. C cycle on UT: C ----------------------------------------------------- DO 801 iut=0,23 ! cycle on UT call blet2(iut,aut) C infile(32:33)=aut !Tamara infile(61:62)=aut !Ljuba UT OPEN(11,FILE=INFILE) c cnttot=0. !!!new for W-map cntneg=0. !!!new for W-map cntpos=0. !!!new for W-map C ===== Input datafile 1311 DO 1799 jj=1,71 ! glats= -87.5:2.5:87.5 c pause ' ' 1 READ (11,12,END=3,ERR=2) (DHR(k),k=0,71) ! W(longs=-180:5:175) 12 FORMAT(72(1X,F4.0)) C DO k=0,71 ! Cycle on Longi <<<<<<<< c== if (dhr(k).gt.90.) dhr(k)=0. C++ Count stormy indices W-: if (dhr(k).lt.-2) then dcnt=1. ! weight=1 for ind=-3, -4 cnts(jj,2,iut)=cnts(jj,2,iut)+dcnt !cnt(Wneg<-2) at lat jj cnts(jj,1,iut)=cnts(jj,1,iut)+dcnt !cnt(Wstr) at lat jj C Counts at each Lati (divided by 72 glongs): clat(jj,2,iut)=clat(jj,2,iut)+dcnt/72. !glat,Ind,UT_iut c!!! cntneg=clat(jj,2,iut) cntneg=cntneg+1. !!!new for W map clat(jj,1,iut)=clat(jj,1,iut)+dcnt/72. ! Range c!!! cnttot=clat(jj,1,iut) cnttot=cnttot+1. !!!new for W map endif C++ Count stormy indices W+: if (dhr(k).gt.2) then dcnt=1. ! weight=1 for ind=+3, +4 cnts(jj,3,iut)=cnts(jj,3,iut)+dcnt !cnt(Wpos>2) at lat jj cnts(jj,1,iut)=cnts(jj,1,iut)+dcnt !cnt(Wstr) at lat jj C Counts at each Lati: clat(jj,1,iut)=clat(jj,1,iut)+dcnt/72. ! Range c!!! cnttot=clat(jj,1,iut) cnttot=cnttot+1. !!!new for W map clat(jj,3,iut)=clat(jj,3,iut)+dcnt/72. ! W+ c!!! cntpos=clat(jj,1,iut) cntpos=cntpos+1. !!!new for W map endif C MAX index at 71 lati: if (rup(iut).lt.dhr(k)) rup(iut)=dhr(k) ! at current lat if (xind(jj,3,iut).lt.dhr(k)) xind(jj,3,iut)=dhr(k) !at lat=jj C MIN index at 71 lati: if (rlo(iut).gt.dhr(k)) rlo(iut)=dhr(k) ! at current lat if (xind(jj,2,iut).gt.dhr(k)) xind(jj,2,iut)=dhr(k) !at lat=jj enddo ! Cycle on Longi >>>>>>>>>>>>> C Range for given Lati, UT: amin=0. ! W<0. bmax=0. !W>0. do k=0,71 ! longs if ((dhr(k).lt.0.).and.(amin.gt.dhr(k))) amin=dhr(k) ! at lat=jj if ((dhr(k).gt.0.).and.(bmax.lt.dhr(k))) bmax=dhr(k) ! at lat=jj enddo cc wpos=bmax*(1.+clat(jj,3,iut)/35.) ! W+ /35????????????? wpos=bmax ! at lat=jj relat(jj,3,iut)=wpos !W+ ! at lat=jj cc wneg=amin*(1.+clat(jj,2,iut)/35.) ! W- /35????????????? wneg=amin ! at lat=jj relat(jj,2,iut)=wneg !W- ! at lat=jj if (clat(jj,1,iut).gt.0.) then ctemp wrange=(bmax-amin)*(1.+clat(jj,1,iut)) ! Range cc wrange=(bmax-amin)*(1.+clat(jj,1,iut)/35.) ! Range /35????????????? c!!! wrange=(bmax-amin)*(1.+clat(jj,1,iut)/10.) ! Range wrange=(bmax-amin) ! at lat=jj else wrange=amax1(bmax,abs(amin)) ! at lat=jj endif relat(jj,1,iut)=wrange ! W Range ! at lat=jj 1799 CONTINUE ! Cycle jj on glats (total W map) 3 CLOSE(unit=11) C---------------------------------------------------------- C average for 71 Lati => Basic Wp for each ut C cneg=0. cpos=0. ctot=0. dneg=0. dpos=0. dtot=0. lat1=1 !cnew lat2=71 !cnew alats=71. !cnew ! 71 glati -87.5:2.5:87.5 c<< grids=alats*72. ! alats*24 hrs UT grids=alats*72. ! alats*glongs (=0:71 glongs)?????????????/ do kk=lat1,lat2 ! dneg=dneg+relat(kk,2,iut) dpos=dpos+relat(kk,3,iut) dtot=dtot+relat(kk,1,iut) ctot=ctot+cnts(kk,1,iut) cneg=cneg+cnts(kk,2,iut) cpos=cpos+cnts(kk,3,iut) C ++++++++++++ enddo ! glati res(iut)=(1.+ctot/grids)*dtot/alats resn(iut)=(1.+cneg/grids)*dneg/alats resp(iut)=(1.+cpos/grids)*dpos/alats if (res(iut).ge.5.0) then totper=ctot/grids*100. write(*,73) AYR,AMN,ADY,iut,cneg,cpos,ctot,totper,res(iut) write(10,73) AYR,AMN,ADY,iut,cneg,cpos,ctot,totper,res(iut) endif c 801 continue ! end cycle on UT 73 format(1X,3A2,1X,I2,3(2X,F5.0),2(2X,F5.1)) C Daily mean: C sum=0. sump=0. sumn=0. do m=0,23 ! cycle on UT ====== 175 sum=sum+res(m) sump=sump+resp(m) sumn=sumn+resn(m) enddo !========= ct sum=sum/12. sum=sum/24. sump=sump/24. sumn=sumn/24. yrmndy='000000' yrmndy(1:2)=AYR yrmndy(3:4)=AMN yrmndy(5:6)=ADY write(12,72) yrmndy,(res(k),k=0,23),sum write(13,72) yrmndy,(resp(k),k=0,23),sump write(14,72) yrmndy,(resn(k),k=0,23),sumn write(*,72) yrmndy,(res(k),k=0,23),sum write(*,72) yrmndy,(resp(k),k=0,23),sump write(*,72) yrmndy,(resn(k),k=0,23),sumn C 500 CONTINUE ! next day C 72 format(A6,24(1X,F4.1),1X,F6.1) ct 75 continue ! <<<<<<<<<<<<<<<<<< CLOSE(unit=10) CLOSE(unit=12) CLOSE(unit=13) CLOSE(unit=14) goto 900 2 CONTINUE write(*,*) 'INPUT FILE IS NOT IN YOUR DIRECTORY ' 900 continue c pause ' ' RETURN END C ============================================================================ subroutine windex(xpar,xmed,delog,iw) C C Procedure to calculate the ionospheric weather W index C C Input: xpar = NeF2 - instantaneous value of the F2 layer peak electron density C or xpar = TECi - instantaneous value of GPS-derived total electron content C xmed = NmF2 - quiet reference (median) value of the F2 layer peak electron density C or xmed = TEC m - quiet reference (median) value of total electron content C C Output: ionospheric weather W index characterizing ionospheric and plasmaspheric activity: C W = 0, +1 or -1 (quiet state) C W = +2 or -2 (moderate disturbance) C W = +3 or -3 (moderate storm) C W = +4 or -4 (intense storm) C C Includes logarithmic scale indices -4,-3,-2,-1,0,1,2,3,4 C according to thresholds of C log10(xpar/xmed)=-0.999,-0.301,-0.155,-0.045,0,0.045,0.155,0.301,0.999 C C Ref.Gulyaeva, T.L.,, I. Stanislwska, and M. Tomasik. Ionospheric weather: Cloning missed C foF2 observations for derivation of variability index. Annales Geophysicae, 26, N.2, C 315-321, 2008. C Gulyaeva, T.L.,and I. Stanislwska, Derivation of a planetary ionospheric storm index. C Annales Geophysicae, 26, N.9, 2645-2648, 2008 . C DIMENSION TR(7) INTEGER IIN(7) data IIN/-3,-2,-1,0,1,2,3/ ! ID-index data TR/-.301,-.155,-.045,0.0 +,.045,.155,.301/ C Temp: c if (delog.ne.0.) goto 5 ! temp if delog is entered ! коммент внесен 08.12.2014 delog=alog10(XPAR)-alog10(XMED) 5 if (delog.gt.tr(7)) then iw=4 goto 2 endif if (delog.eq.0.) then iw=0 goto 2 endif if(delog.le.tr(1)) then iw=-4 ! W-index goto 2 endif if (delog.ge.tr(7)) then iw=4 ! W-index goto 2 endif C do 1 n=1,7 if (delog.lt.0.) then if ((delog.ge.tr(n)).and.(delog.lt.tr(n+1))) then iw=iin(n) ! W-index endif else if ((delog.gt.tr(n-1)).and.(delog.le.tr(n))) then iw=iin(n) ! W-index endif endif 1 continue C 2 continue C RETURN END C C C ===============================================================