C.......................................................May 2014 C Plas-cor using Geomagnetic Local Time, GMLT C Feb. 2014 C C Count global CCIR-based median afc(UT) and ahc*UT) from CCIR maps fcF2 and hcF2(from M3000 map) C for given UT subroutine subsmi4(AYEAR,AMN,ADY,IMEDF,IMEDH,IRIT) C C DIMENSION im(12),av(0:72) +,IMEDF(71,0:72,0:23),IMEDH(71,0:72,0:23) ! results of median maps +,resmodip(71,0:72) ! modip data for IONEX map +,IRIT(71,0:72,0:23) ! TEC-NeQ-P dimension R12y(0:12),F2(13,76,2),FM3(9,49,2) +,F2n(13,76,2),FM3n(9,49,2) !current and prev/next month integer month CHARACTER*48 R12file !Ljuba character*2 ayr,amn,ady character*4 AYEAR C++ integer*2 iyr,imn,idy,iut DATA IM/31,28,31,30,31,30,31,31,30,31,30,31/ common/BLOKM/resmodip COMMON /bplas/r12i,glat,ylon,ut,ldaynr,heitop,eltop C c PHI=ATAN(1.0)*4. UMR=PHI/180. XFOF2=0. XHMF2=0. XHI=0.0 c C ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ c AYR=AYEAR(3:4) read(ayear,*) xyear iyear=int(xyear) read(ayr,*) yr read(amn,*) rmn imn=int(rmn) month=int(rmn) C iyr=int(yr) if(int(yr/4.)*4.eq.iyr) THEN im(2)=29 ELSE IM(2)=28 ENDIF ndnr=im(imn) C read(ady,*) rdy idy=int(rdy) ldaynr=ndoy(iyr,imn,idy) C mnmidy=im(imn)/2 if (idy.lt.mnmidy+1) then monthn=month-1 if (monthn.eq.0) monthn=12 else monthn=month+1 if (monthn.gt.12) monthn=1 endif C Read R12 for current year: R12file='/var/www/izmiran/ionosphere/weather/graf/R12.dat' !Ljuba open(19,file=R12file,status='OLD') 211 read(19,*, end=99) j,(R12y(k),k=1,12) if (j.lt.iyear) then R12y(0)=R12y(12) goto 211 endif R12=R12y(month) if ((monthn.eq.12).and.(month.eq.1)) then R12n=R12y(0) else R12n=R12y(monthn) endif if (R12.gt.150.0) R12=150.0 if (R12n.gt.150.0) R12n=150.0 99 continue close(unit=19) dyri=idy-15.15 R12i=R12+(R12n-R12)*dyri/30.3 ! SSN for current day cov=63.7+(0.728+8.9E-4*R12)*R12 if (cov.gt.193.) cov=193. if (cov.lt.63.0) then cov=63.0 endif covn=63.7+(0.728+8.9E-4*R12n)*R12n if (covn.gt.193.) covn=193. if (covn.lt.63.0) then covn=63.0 endif C call subccir(month,monthn,F2,F2n,FM3,FM3n) !Read CCIR coefficients C call submodip ! read modip for IONEX grids C UT cycle DO 203 iut=0,23 UT=float(iut) C C GLAT cycle C glat=-90. do 301 lat=1,71 glat=glat+2.5 C do k=0,72 av(k)=0. ! enddo C GLONG cycle glon=-185. do lon=0,72 glon=glon+5.0 if (glon.lt.0.) then ylon=glon+360. else ylon=glon endif idmn=im(month) smodip=resmodip(lat,lon) C+ Geomagnetic local time UTS=ut*3600. ! UT /sec. xday=float(ldaynr) gmlt=fmlt(uts,glat,ylon,xday) C call prepccir(glat,ylon,smodip,month,monthn,idy,idmn,UT,R12,R12n +,cov,covn,F2,F2n,FM3,FM3n,hm,fm,foe,foF1,teci) C IMEDF(lat,lon,iut)=nint(fm*10.) ! results of CCIR foF2 map IMEDH(lat,lon,iut)=nint(hm) ! results of CCIR hmF2 map IRIT(lat,lon,iut)=nint(teci) ! results of IRI-Plas CCIR TEC map TECU*10. enddo ! C c write(*,109) (iav(k),k=0,72) 301 continue ! end lat cycle of median map write(*,*) ' UT= ',iut 203 continue ! end ut cycle 108 format(3A2,24(1X,I4)) ! 109 format(73(1X,I4)) ! C GOTO 1234 2 CONTINUE 1234 continue ! c// write(*,*) infile,outfile c// pause ' ' 34 RETURN END C============================================================================== C SUBROUTINE SUBCCIR(IMONTH,IMONTHN,F2,F2n,FM3,FM3n) C Calculation of coefficients C FOF2(76,13) ; SM3000(49,9) C for level of solar activity, MODIP and current month + prec/next month C using CCIR maps DIMENSION F2(13,76,2),FM3(9,49,2),F2n(13,76,2),FM3n(9,49,2) C+ SMI input of CCIR coefficients CHARACTER*46 binfile !Ljuba binfile='/var/www/izmiran/ionosphere/weather/graf/BIN16' !Ljuba !!! C C OPEN (UNIT=16,FILE=binfile,STATUS='OLD', C * ACCESS='DIRECT',FORM='UNFORMATTED',RECL=11432) ! Ljuba C C READ (16,rec=imonth) F2,FM3 C REWIND 16 C READ (16,rec=imonthn) F2n,FM3n C REWIND 16 C CLOSE (UNIT=16) C C C OPEN (UNIT=16,FILE=binfile,STATUS='OLD', * ACCESS='DIRECT',FORM='UNFORMATTED',RECL=11432) ! Ljuba READ (16,rec=imonth) F2,FM3 CLOSE (UNIT=16) ! Ljuba OPEN (UNIT=16,FILE=binfile,STATUS='OLD', * ACCESS='DIRECT',FORM='UNFORMATTED',RECL=11432) READ (16,rec=imonthn) F2n,FM3n CLOSE (UNIT=16) RETURN END C subroutine cciri(xMODIP,mth,UT,R12,alat,along,foF2,M3000,F2x,FM3x) implicit real (a-h,o-z) real M3000 dimension FF0(988),xm0(441) +,F2x(13,76,2) +,FM3x(9,49,2) integer QM(7),QF(9) save data QF/11,11,8,4,1,0,0,0,0/,QM/6,7,5,2,1,0,0/ data montha,monthb,Rga/13,14,-10.0D0/ if (mth.ne.montha) then montha=mth endif if (R12.ne.Rga.or.mth.ne.monthb) then RR2=R12/100.0 RR1=1.0-RR2 do i=1,76 do j=1,13 k=j+13*(i-1) FF0(k)=F2x(j,i,1)*RR1+F2x(j,i,2)*RR2 enddo enddo do i=1,49 do j=1,9 k=j+9*(i-1) xm0(k)=FM3x(j,i,1)*RR1+FM3x(j,i,2)*RR2 enddo enddo Rga=R12 monthb=mth endif foF2= gamma1(xMODIP,alat,along,UT,6,QF,9,76,13,988,FF0) M3000=gamma1(xMODIP,alat,along,UT,4,QM,7,49, 9,441,xm0) return end C-------------------------------------------------- C-------------------------------------------------- real function peakh(foE,foF2,M3000) real MF,M3000 sqM=M3000*M3000 MF=M3000*sqrt((0.0196*sqM+1.)/(1.2967*sqM-1.0)) If(foE.ge.1.0E-30) then ratio=foF2/foE ratio=djoin(ratio,1.75,20.0,ratio-1.75) dM=0.253/(ratio-1.215)-0.012 else dM=-0.012 endif peakh=1490.0*MF/(M3000+dM)-176.0 return end C-------------------------------------------------- REAL FUNCTION GAMMA1(SMODIP,SLAT,SLONG,HOUR,IHARM,NQ, & K1,M,MM,M3,SFE) C CALCULATES GAMMA1=FOF2 OR M3000 USING CCIR NUMERICAL MAP C COEFFICIENTS SFE(M3) FOR MODIFIED DIP LATITUDE (SMODIP/DEG) C GEOGRAPHIC LATITUDE (SLAT/DEG) AND LONGITUDE (SLONG/DEG) C AND UNIVERSIAL TIME (HOUR/DECIMAL HOURS). C NQ(K1) IS AN INTEGER ARRAY GIVING THE HIGHEST DEGREES IN C LATITUDE FOR EACH LONGITUDE HARMONIC. C M=1+NQ1+2(NQ2+1)+2(NQ3+1)+... . C SHEIKH,4.3.77. REAL*8 C(12),S(12),COEF(100),SUM DIMENSION NQ(K1),XSINX(13),SFE(M3) COMMON/CONST/UMR HOU=(15.0*HOUR-180.0)*UMR S(1)=SIN(HOU) C(1)=COS(HOU) DO 250 I=2,IHARM C(I)=C(1)*C(I-1)-S(1)*S(I-1) S(I)=C(1)*S(I-1)+S(1)*C(I-1) 250 CONTINUE DO 300 I=1,M MI=(I-1)*MM COEF(I)=SFE(MI+1) DO 300 J=1,IHARM COEF(I)=COEF(I)+SFE(MI+2*J)*S(J)+SFE(MI+2*J+1)*C(J) 300 CONTINUE SUM=COEF(1) SS=SIN(SMODIP*UMR) S3=SS XSINX(1)=1.0 INDEX=NQ(1) DO 350 J=1,INDEX SUM=SUM+COEF(1+J)*SS XSINX(J+1)=SS SS=SS*S3 350 CONTINUE XSINX(NQ(1)+2)=SS NP=NQ(1)+1 SS=COS(SLAT*UMR) S3=SS DO 400 J=2,K1 S0=SLONG*(J-1.)*UMR S1=COS(S0) S2=SIN(S0) INDEX=NQ(J)+1 DO 450 L=1,INDEX NP=NP+1 SUM=SUM+COEF(NP)*XSINX(L)*SS*S1 NP=NP+1 SUM=SUM+COEF(NP)*XSINX(L)*SS*S2 450 CONTINUE SS=SS*S3 400 CONTINUE GAMMA1=SUM RETURN END C C C C************************************************************ subroutine prepccir(alat,along,xmodip,mth,mth1,idy,idmn,UT,R12 +,R12n,flx,flxn,F2,F2n,FM3,FM3n,hmf2,fof2,foe,foF1,teci) real M3000,M3000a,M3000b C++ integer*2 idy dimension hm(3),aep(3),bb(6) dimension + F2(13,76,2) +,FM3(9,49,2) +,F2n(13,76,2) +,FM3n(9,49,2) save parameter (pi12=0.26179939) parameter (DR=1.745329252E-2,RD=57.29577951) COMMON /bplas/rzs,plat,plong,put,ldaynr,heitop,eltop data UT0,mth0/-100.0,-1/ ut1=mod(UT+24.0,24.0) if (ut1.ne.UT0.or.mth1.ne.mth0) then call sdec(mth1,ut1,sdelta,cdelta) UT0=ut1 mth0=mth1 endif if (idy.lt.idmn/2+1) then C 30.3 days between mid-prec month and mid-current month dyri=15.15+idy else C 30.3 days between mid-current month and mid-next month dyri=idy-15.15 endif xlt=mod(ut1+along/15.0+24.0,24.0) cchi=sin(alat*DR)*sdelta-cos(alat*DR)*cdelta* * cos(pi12*xlt) cchinon=sin(alat*DR)*sdelta-cos(alat*DR)*cdelta* * cos(pi12*12.0) chi=atan2(sqrt(1.0-cchi*cchi),cchi)*RD chinon=atan2(sqrt(1.0-cchi*cchinon),cchinon)*RD call cciri(xMODIP,mth,UT0,R12,alat,along,foF2a,M3000a,F2,FM3) call cciri(xMODIP,mth1,UT0,R12n,alat,along,foF2b,M3000b,F2n,FM3n) fof2=foF2a+(foF2b-foF2a)*dyri/30.3 abslat=abs(alat) foEa=FOEEDI(flx,chi,chinon,abslat) foEb=FOEEDI(flxn,chi,chinon,abslat) foe=foea+(foeb-foea)*dyri/30.3 efoF1a=ef1r1(foF2a,foEa) efoF1b=ef1r1(foF2b,foEb) foF1=efoF1a+(efoF1b-efoF1a)*dyri/30.3 M3000=M3000a+(M3000b-M3000a)*dyri/30.3 call prepmdgr(R12,foF2,foF1,foE,M3000,hm,bb,aep) hmF2=hm(1) C call prepPlas(hm,aep,bb) C call subinteg(20200.0,teci,hm,aep,bb) ! vTEC [80:20000]km return end C----------------------------------------------------------------- subroutine submodip C TLG.............................................March 2014 C Read data from file 'modip2.asc' dimension resmodip(71,0:72) character*51 filenam ! Ljuba заменила 20 на 51 common/BLOKM/resmodip c C filenam='modip2.asc' filenam='/var/www/izmiran/ionosphere/weather/graf/modip2.asc' !Ljuba open(21,file=filenam,form='FORMATTED') do i=0,72 read(21,*) (resmodip(j,i),j=1,71) enddo close(unit=21) return end C----------------------------------------------------------------- subroutine sdec(mth,UT,sdelta,cdelta) parameter (DR=1.745329252E-2) doy=mth*30.5-15.0 t =doy + (18.0-UT)/24.0 amrad=(0.9856*t - 3.289)*DR aLrad = amrad + (1.916*sin(amrad)+0.02*sin(2.0*amrad)+ + 282.634)*DR sdelta=0.397820*sin(aLrad) cdelta=sqrt(1.0-sdelta*sdelta) return end C ---------------------------------------------------- function Fne(x) c- FNe(X)=0.124*X*X FNe=0.124*X*X return end C ---------------------------------------------------- function FEpst(X,Y,Z,W) c- FEpst(X,Y,Z,W)=X*fexp((W-Y)/Z)/(1.0+fexp((W-Y)/Z))**2 FEpst=X*fexp((W-Y)/Z)/(1.0+fexp((W-Y)/Z))**2 return end C ---------------------------------------------------- subroutine prepmdgr(R12,foF2,efoF1,efoE,M3000,hm,bb,aep) real NmE,NmF1,NmF2,M3000 dimension aep(3),hm(3),bb(6) data hmE /120.0/ NmF2=FNe(foF2) NmF1=FNe(efoF1) NmE= FNe(efoE) hmF2=peakh(efoE,foF2,M3000) hmF1=0.5*(hmF2+hmE) hm(1)=hmF2 hm(2)=hmF1 hm(3)=hmE dNdHmx=-3.467+1.714*log(foF2)+2.02*log(M3000) dNdHmx=exp(dNdHmx) B2bot=0.385E2*NmF2/dNdHmx hF2F1E=hmF2-hmF1 B1top=0.3*hF2F1E B1bot=0.5*hF2F1E Betop=B1bot if (Betop.lt.7.0) Betop=7.0 Bebot=5.0 aep(1)=4.0*NmF2 aep(2)=4.0*NmF1 aep(3)=4.0*NmE if (efoF1.le.0.5) then aep(2)=0.0 aep(3)=4.0*(NmE-FEpst(aep(1),hmF2,B2bot,hmE)) else do izml=1,5 aep(2)=4.0*(NmF1-FEpst(aep(1),hmF2,B2bot,hmF1) & -FEpst(aep(3),hmE, Betop,hmF1)) aep(2)=djoin(aep(2),0.8*NmF1,1.0,aep(2) & -0.8*NmF1) aep(3)=4.0*(NmE -FEpst(aep(2),hmF1,B1bot,hmE) & -FEpst(aep(1),hmF2,B2bot,hmE)) enddo if (aep(2).lt.0.0) then aep(2)=0.0 aep(3)=4.0*(NmE-FEpst(aep(1),hmF2,B2bot,hmE)) endif endif aep(3)=djoin(aep(3),0.5E-2,60.0,aep(3)-0.5E-2) b2k = 3.22-0.538E-1*foF2-0.664E-2*hmF2+0.113*hmF2/B2bot + +0.257E-2*R12 b2k=djoin(b2k,1.0,2.0,b2k-1.0) B2top=b2k*B2bot bb(1)=Bebot bb(2)=Betop bb(3)=B1bot bb(4)=B1top bb(5)=B2bot bb(6)=B2top return end C ------------------------------------------------------ REAL FUNCTION FOEEDI(COV,XHI,XHIM,XLATI) C------------------------------------------------------- C CALCULATES FOE/MHZ BY THE EDINBURGH-METHOD. C INPUT: MEAN 10.7CM SOLAR RADIO FLUX (COV), GEOGRAPHIC C LATITUDE (XLATI/DEG), SOLAR ZENITH ANGLE (XHI/DEG AND C XHIM/DEG AT NOON). C REFERENCE: C KOURIS-MUGGELETON, CCIR DOC. 6/3/07, 1973 C TROST, J. GEOPHYS. RES. 84, 2736, 1979 (was used C to improve the nighttime varition) C D.BILITZA--------------------------------- AUGUST 1986. COMMON/CONST/UMR C variation with solar activity (factor A) ............... A=1.0+0.0094*(COV-66.0) C variation with noon solar zenith angle (B) and with latitude (C) SL=COS(XLATI*UMR) IF(XLATI.LT.32.0) THEN SM=-1.93+1.92*SL C=23.0+116.0*SL ELSE SM=0.11-0.49*SL C=92.0+35.0*SL ENDIF if(XHIM.ge.90.) XHIM=89.999 B = COS(XHIM*UMR) ** SM C variation with solar zenith angle (D) .......................... IF(XLATI.GT.12.0) THEN SP=1.2 ELSE SP=1.31 ENDIF C adjusted solar zenith angle during nighttime (XHIC) ............. XHIC=XHI-3.*ALOG(1.+EXP((XHI-89.98)/3.)) D=COS(XHIC*UMR)**SP C determine foE**4 ................................................ R4FOE=A*B*C*D C minimum allowable foE (sqrt[SMIN])............................... SMIN=0.121+0.0015*(COV-60.) SMIN=SMIN*SMIN IF(R4FOE.LT.SMIN) R4FOE=SMIN FOEEDI=R4FOE**0.25 RETURN END c C ---------------------------------------------------- C ---------------------------------------------------- real function ef1r1(foF2,efoE) parameter (DR=1.745329252E-2) parameter (chi0=86.232928) efoF1=djoin(1.4*efoE,0.0,1000.0,efoE-2.0) efoF1=djoin(0.0,efoF1,1000.0,efoE-efoF1) ef1r1=djoin(efoF1,0.85*efoF1,60.0,0.85*foF2-efoF1) if (ef1r1.lt.1.0E-6) ef1r1=0.0 return end C----------------------------------------------------------------- real function djoin(f1,f2,alpha,x) real f1,f2,alpha,x,ee,fexp ee=fexp(alpha*x) djoin=(f1*ee+f2)/(ee+1.0) return end C----------------------------------------------------------------- real function fexp(a) real a if(a.gt.80.0) then fexp=5.5406E34 return endif if(a.lt.-80.0D0) then fexp=1.8049E-35 return endif fexp=exp(a) return end C----------------------------------------------------------------- real function fmlt(ut,xlat,xlong,day) C Calculats geomagnetic local time in hours for given C universal time UT(sec) C XLAT/XLONG=geodetic Latitude/Longitude (deg.) C DAY= day of year C D. Bilitza, COSPAR-1980, paper 7.3.1 C rad=57.29578 flat=xlat/rad flong=xlong/rad delta=0.409207*sin((day-80.0)*3.14159/184.0) beta=(180.-ut/240.0)/rad px=cos(flat)*cos(flong) py=cos(flat)*sin(flong) pz=sin(flat) sx=cos(delta)*cos(beta) sy=cos(delta)*sin(beta) sz=sin(delta) p1=0.35117*px-0.91483*py-0.19937*pz p2=0.93358*px+0.35837*py s1=0.35117*sx-0.91483*sy-0.19937*sz s2=0.93358*sx+0.35837*sy thetp=atan2(p2,p1) thets=atan2(s2,s1) fmlt=(thetp-thets+3.141593)/.2617994 afmlt=anint(fmlt*10.) fmlt=afmlt/10. if(fmlt.ge.24.) fmlt=fmlt-24. if(fmlt.lt.0.) fmlt=fmlt+24. return end C C ---------------------------------------------------------------------