C........................................................Nov. 2021 C Use R12 instead GEC12 as solar proxy C........................................................July 2018 C Add map GIM-Hc center-of-mass production C........................................................Dec. 2015 C Produce GIM-TEC, GIM-foF2, GIM-hmF2, GIM-TECeff. in magn. coords. C IONEX mlats=-87.5:2.5:87.5, mlongs=0:5:360. C........................................................Apr. 2014 C replace hmF2 calculation with C hmF2=peakh(foE,foF2,XM3000) C........................................................Oct. 2013 C Exclude Hakura Tables C Use GEC array instead of IG12 array C C........................................................Nov.2012 C NEW!!!!!!!!! C Produce maps foF2, hmF2, TEC, TECeff in IONEX style format C Add calculation of TEC(h) profile for [0 20200]km in step of 200 km C multiplied by Cv(h) coefficients to produce GEC C Add output of TECeff map with TECeff equal to sum of TEC(h)*Cv(h) C C input i=1,2,...71 C longi:-180.0, -175.0,... 0...180.0 ==> j=0,...,72 C output maps at C lats=[-87.5, -85.0,...,0,...,87.5 for i= 71, 70,...1 C longi:-180.0, -175.0,... 0...180.0 ==> j=0,...,72 C C C C .......................................................Feb. 2012 C C NO!!!!!!!!! Include GEC calculation for UT=0,1,...,24 hrs UT C NO!!!!!!!!! outfile='gecyr.br' C C C T.L.Gulyaeva...........................................June 2011 C CInclude input of IONEX_1h files C Produce fcF2 & hcF2 numerical maps with IRI-Plas code C for glats [-87.5:2.5:87.5], glons [-180:5:180] C in 48 daily files: , for ut=0,1,...,23hUT C PROGRAM MAPIRIPC9jm C---------------------------------------------------------------- C Select NEW+++:Input 71 (instead of 35) geographic latitudes: res1=87.5,85,,.0..-87.5N C at: NEW+++: 0,1,...,72 geographic longitudes [-180,-175,...,180] = IONEX step C Produce maps fc & hc for 0, 1,..., 23hr UT.syr C 2 output arrays for each hourUT =48 daily files C ,,, C send results to d:/web/dfc/yr/mn/fcmndyut.syr; C d:/web/dhc/yr/mn/hcmndyut.syr; C program for EXTRACTING FROM "JPLGdoy0.iYR" ionex_TEC data C +++++++++ C Add+Input: selected glats& glons for mlati=+60,+55,...,-60 deg. at MLT=0,1,...,23h MLT. C++ Add+preprocessing results: TEC(glats=87.5:-2.5-87.5,glons=-180:15:165) Array TECPRE(72,73) C C*********TEST dimension outf(0:11,0:120),oarr(50),indap(13),akp(13) +,iprekp(0:7),icurkp(0:7) REAL tech(200,2) Cold +,rkp(0:366,0:7) logical jf(30),F1REG integer daynr,numhei,ifcf2(71,0:72),ihcf2(71,0:72),itec(71,0:72) +,ieff(71,0:72),ihcm(71,0:72) CHARACTER*2 YY,AMN,DY1,DY2,CMN CHARACTER*4 TTT DIMENSION GLATS(71,0:72),GLONS(71,0:72),resgec(72) ! GEC results for each latitude INTEGER*2 iyr,imn,idy,nndy,iut,ndyend +,iut0 C? character*4 xxpr,yypr !* C? COMMON /BXY/xxpr,yypr !* common /block1/hmf2,xnmf2,hmf1,F1REG &/PLAS/HPL,XNEPL,XKP,ICALLS,RUT,RLT,OARR,TECI,TCB,TCT,TCPL,TEC COMMON /CONST/UMR/ARGEXP/ARGMAX & /CONST1/HUMR,DUMR,FAP,INDAP,AKP,COV,RZS,HMLON &/dem/W,glong,hour,daynr,cn1000,alon,blon,ENRE,HSC C# new common block: & /QTOP/Y05,H05TOP,QF,XNETOP,XM3000,HHALF,TAU,HTOP,HCM,DHTOP,TCM !Hcm -center mass,Tcm=TEC(Hcm) common /temod/month ! added in IRIS2006.for to call Temodel C? common /tecprof/tech COMMON /BL1/AYEAR,YY,AMN,DY1,DY2,DD1,DD2,ID1,ID2,ADY COMMON /BL2/ DD,TTT,ZZ Cold COMMON /FIKP/rkp,prekp COMMON /FIKP/icurkp,iprekp,iyear_cur,imn_cur,idy_cur !*TEST=new common/BLMAG/glats,glons,itec common /tecprof/tech C C**********TEST CHARACTER*128 OUTFILEF,OUTFILEH,OUTFILET,OUTFILEC !Hc center-of-mass +,OUTFILEE ! TECeff EXTERNAL XE1,XE2,XE3,XE3_1,XE4,XE5,XE6,XXE6,TEDER,ABLON CHARACTER(10) DD CHARACTER(10) TT CHARACTER(5) ZZ CHARACTER*4 AYEAR CHARACTER*60 inglats,inglons character*2 ayr,ady,aut,ady0 character*3 aday dimension im(12) C? +,teceff(71,0:72) ! TECeff DATA IM/31,28,31,30,31,30,31,31,30,31,30,31/ C++ C++ ICALLS=0 CALL DATE_AND_TIME(DATE=DD,TIME=TT,ZONE=ZZ) TT(5:)=' ' WRITE(*,*) 'PC Date: Year,Month,Day = ',DD,'Time = ',TT C WRITE(*,'(\A\)') ' ENTER YEAR =' READ(*,'(A)') AYEAR ctemp AYEAR='2012' !ctemp ayr=ayear(3:4) WRITE(*,'(\A\)') ' ENTER MN =' READ(*,'(A)') AMN ctemp amn='01' WRITE(*,*) 'Enter ndy1, ndy2 : ' read(*,987) nndy,ndyend 987 format(I2) call blet2(nndy,DY1) call blet2(ndyend,DY2) ady0=DY1 ctemp nndy=2 ! temp for test ctemp ndyend=5 !temp for test C c// outfile='ghYRMN.txt' irunkp=0 C 1001 DO I=1,30 JF(I)=.TRUE. ENDDO C C Read geogr. coords. for magn. grids inglats='glatgrid1.txt' inglons='glongrid1.txt' 171 format(73(1X,F5.1)) OPEN(10,file=inglats) OPEN(11,file=inglons) C do i=1,71 read(10,171) (glats(i,k),k=0,72) do k=0,72 if (glats(i,k).gt.87.5) glats(i,k)=87.5 if (glats(i,k).lt.-87.5) glats(i,k)=-87.5 enddo read(11,171) (glons(i,k),k=0,72) enddo close(unit=10) close(unit=11) C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C Cd UMR=ATAN(1.0)*4./180. IRUN=0 ! NORMAL COMMAND iut0=0 ! NORMAL COMMAND CSTART TEMP: Next two command for starting UT (arbitrary) C/ iut0=23 !TEMP C/ irun=23 !TEMP ady0=ady CEND TEMP 118 format(A4) cc pause ' ' C C ========================================= C? do k=0,72 C? do i=1,71 C? teceff(i,k)=0. ! arrange TECeff for mlong 0, 5..,360 deg.E C? enddo C? enddo c read(ayear,*) xyear IYYYY=int(xyear) yy=ayr read(yy,*) ryr iyr=int(ryr) c read(amn,*) xmn imn=int(xmn) CC call blet2(nndy,dy1) ady0=dy1 write(*,199) dy1 199 format(' day=',A2) read(dy1,*) xdy idy1=int(xdy) if (irunkp.eq.0) then CMN=AMN Cold CALL SUBKP CALL SUBKP2018 AMN=CMN irunkp=1 endif idy2=idy1 c nrday=int(daynr) c outfiles ================================= outfilef='d:\web\dfs\yr\mn\fomndyut.jyr' ! ionex type map in magn. coord. outfilef(12:13)=ayr outfilef(15:16)=amn outfilef(20:21)=amn C outfilef(22:23)=ady1 outfilef(28:29)=ayr outfileh=outfilef outfileh(9:9)='h' outfileh(18:18)='h' outfilet=outfilef outfilet(9:9)='t' outfilet(18:18)='t' C>> outfilet(18:18)='d' ! Residual error: devTEC=TECobs-TECmod outfilee=outfilet outfilee(9:9)='e' outfilee(18:18)='e' ! eomndyut.jyr ! TECeff for GEC calculation C outfilec=outfileh outfilec(10:10)='c' outfilec(18:19)='co' ! \comndyut.jyr = center-of-mass in magn. coord C++ ist=0 JOUT=0 JSTR=0 JUC=0 C> JMAG=0 JMAG=1 ! in magn. coords RZS=-1. ! temp HBEG=80. jf(25)=.true. oarr(33)=rzs JF(3)=.FALSE. JF(4)=.FALSE. IF (JUC.EQ.2) JF(5)=.FALSE. C IF (JOUT.EQ.0) JF(18)=.FALSE. IF (JSTR.EQ.0) JF(26)=.FALSE. iiyr=int(xyear) if(int(xyear/4.)*4.eq.iiyr) THEN im(2)=29 ELSE IM(2)=28 ENDIF c xdy=1. ! start day-of-month CC CC GOTO 801 !TEMPPPPPPPPPPPPPPPPPPP C Cycle day-to-day +++++++++++++++++++++++++++++ 700 DO 800 idy=nndy,ndyend call blet2(idy,ADY) outfilef(22:23)=ady outfileh(22:23)=ady outfilet(22:23)=ady outfilee(22:23)=ady outfilec(22:23)=ady lda=ndoy(iyr,imn,idy) daynr=lda 178 call blet3(lda,ADAY) C 113 CONTINUE 27 format(5X,F3.1) C C Prepare results of GEC for UT=0,...,23 + daily mean 114 CONTINUE C C Cycle on UT: C DO 775 iut=iut0,23 cc write(*,89) ryear,rmn,rdy,iut rut=float(iut) call blet2(iut,AUT) C Transform GIM-TEC to magn. coords for given day and UT : call subexmag4(AYR,AMN,ADY,AUT) outfilef(24:25)=aut ! true outfileh(24:25)=aut ! true outfilet(24:25)=aut ! true outfilee(24:25)=aut ! true outfilec(24:25)=aut ! true C C>Move next commands to subkp2018: C>prepare extracting UKP: C> do i=0,7 C> curkp(i)=rkp(lda,i) C> prekp(i)=rkp(lda-1,i) C> enddo C extract UKP=maxkp for 6 prec.hrs from rkp array: ukp=0. iut1=iut/3-1 iut2=iut1-1 if (iut1.ge.0) then ukp1=icurkp(iut1)/10. else ukp1=iprekp(iut1+8)/10. endif if (iut2.ge.0) then ukp2=icurkp(iut2)/10. else ukp2=iprekp(iut2+8)/10. endif UKP=MAX(ukp1,ukp2) XKP=UKP C +++++++++++ jt=0 C C ilat=0 C start cycle reading TEC on glatitudes: gmlati=-90.0 do k=1,72 resgec(k)=0. enddo C DO 777 lat=1,71 ! Start cycle on mlati ++++++++++++++++++++++++++++++++ ilat=ilat+1 ! ilat=1 for -87.5 gmlati=gmlati+2.5 ALATI=gmlati ! for iriplas 91 j1=0 C ilatgec=ilat CC C Select values for MLong=0,5,...360!!!!!!!!!!!!: C ii=-1 jj=71 C> 9 continue C==================================================== do 177 il=0,360,5 ! for mlong-cycle ii=ii+1 ! ii=0 for mlong=0,...,ii=72 for mlong=360 ires=ii c if ((alati.eq.-10.0).or.(alati.eq.-7.5)) then !TEMP+++++++++++++ C/ if (alati.eq.-17.5) then !TEMP+++++++++++++ C/ ifcf2(lat,ires)=ifcf2(lat-1,ires) !TEMP+++++++++++++ C/ ihcf2(lat,ires)=ihcf2(lat-1,ires) !TEMP+++++++++++++ C/ itec(lat,ires)=itec(lat-1,ires) !TEMP+++++++++++++ C/ ieff(lat,ires)=ieff(lat-1,ires) !TEMP+++++++++++++ C/ goto 177 !TEMP+++++++++++++ C/ endif !TEMP+++++++++++++ gmlong=float(il) CC++ ALONG=gmlong if (along.lt.0.) along=along+360. MMDD=imn*100+idy HOURS=IUT+25. HMF2I=0. FOF2=0. JF(8)=.TRUE. JF(9)=.TRUE. TECI=itec(lat,ii)/10. IF (TECI.LT.2.) TECI=0. IF (TECI.GT.0.) JF(30)=.FALSE. tec = -111. tect= -111. tecb= -111. TCPL= -111. hmf2=0. xm3000=0. CC Call iriplas for fcF2 & hcF2 calculations c c initialize IRI parameter in COMMON blocks c abeg=hbeg C aend=hend AEND=1336. astp=aend-hbeg numhei=101 ! 101 ranges = 200 km CREM jind=5 ! TECgn index C_PRE jind=0 ! gec_rz.dat file input jind=1 !SSN1 C jf(30)=.FALSE. ! TECI>0. C call IRIS2018c(JF,JMAG,ALATI,ALONG,IYYYY,MMDD,HOURS, & abeg,aend,numhei,OUTF,JIND) C Add: if (JF(30)) GOTO 123 ! No TECI input C C Add: 2nd call IRIS2015c: replace TEC input with 'input' of updated foF2, hmF2: C jf(30)=.TRUE. ! no TECI option TECI=0. JF(8)=.FALSE. ! foF2>0. JF(9)=.FALSE. ! hmF2>0. JF(30)=.TRUE. !TECI=0. tec = -111. tect= -111. tecb= -111. TCPL= -111. icalls=icalls-1 call IRIS2018c(JF,JMAG,ALATI,ALONG,IYYYY,MMDD,HOURS, !* C_pre & abeg,aend,numhei,OUTF,JINP) !* & abeg,aend,numhei,OUTF,JIND) !* C End Add commands C 123 IF(JF(8)) XNMF2=OARR(1) IF(JF(9)) HMF2=OARR(2) IF(JF(17)) RZS=OARR(33) AL=OARR(32) ! L value IF (AL.GT.99.99) AL=99.99 XHI=OARR(23) YMODI=OARR(27) CCC GLATI=OARR(37) GLONG=OARR(38) rlt=iut+glong/15. if (rlt.lt.0.) rlt=rlt+24. XMLAT=OARR(42) XMLON=OARR(40) GMLT=OARR(44) HTOP=OARR(46) XNMHC=0.5*XNMF2 ! NEW: Ne ~ 0.5*NMF2 for Hc limits Test XNMF2=XNMF2/1.0E+06 FOF2=SQRT(XNMF2/1.24E+04) TOT=TCB+TCT+TCPL TAU=(TOT/XNMF2)*1.0E+07 C++++++++++++++++++++++++++++++++++++++++++++++++++ C Center-of-mass Hc calculation C Obtain: hbothc(at XNMHC below hmF2) and htophc(at XNMHC above hmF2) C? h1=120. C? tec1=0. C TCM=TCM*1.0E+10 C? TEC_01=0. C? DO 33 i=50,0,-1 C? h0=OUTF(0,I+1) C? XNE0=OUTF(1,I+1) C? h1=OUTF(0,I) C? XNE1=OUTF(1,I) C? TEC_01=TEC_01+(XNE0+XNE1)/2.*(h1-h0)*1.0E+03 C? i1=50-i C? OUTF(11,I)=TEC_01 C? 33 CONTINUE ifcf2(lat,ires)=nint(foF2*10.) ihcf2(lat,ires)=nint(hmF2) ieff(lat,ires)=nint(OARR(13)*10.) ihcm(lat,ires)=nint(Hcm) C>> TEMP C>> output TEC residual error = TECobs-TECmod C>> itec(lat,ires)=nint((TECI-TOT)*100.) 177 continue C test write(*,181) iut,ALATI cc write(*,182) (ifcf2(lat,k),k=0,24) cc write(*,182) (ihcf2(lat,k),k=0,24) cc write(*,182) (itec(lat,k),k=0,24) !TEMP 777 continue ! End cycle on MLATs CREM rgec(iut)=rgec(iut)+GECUT/1.0E+32 181 format(2X,I2,2X,F5.1) C C ================================================== C 1611 continue ! cycle C --------------- C Prepare output of results for nnlong=0,1,...,24 OPEN(12,FILE=OUTFILEF) OPEN(13,FILE=OUTFILEH) OPEN(14,FILE=OUTFILET) OPEN(15,FILE=OUTFILEE) OPEN(16,FILE=OUTFILEC) C OUTPUT OF RESULTS DO nnlat=1,71 write(12,182) (ifcf2(nnlat,k),k=0,72) write(13,182) (ihcf2(nnlat,k),k=0,72) crem write(*,182) (itec(nnlat,k),k=0,72) write(14,182) (itec(nnlat,k),k=0,72) !OUT of TEC hourly maps write(15,182) (ieff(nnlat,k),k=0,72) write(16,182) (ihcm(nnlat,k),k=0,72) ENDDO 182 format(73(1X,I4)) close(unit=12) close(unit=13) close(unit=14) close(unit=15) close(unit=16) 775 continue ! End calculation for given UT C C> C Continue to the next UT hr : 32 irun=irun+1 C> ady0=rdy C> if (irun.lt.iutend) goto 80 ! continue to the next UT hr =>>>> c if (iut.lt.22) goto 80 ! continue to the next UT hr =>>>> C pause ' ' C++++++++++++++++++++++++++++ C C C C 101 irun=0 c pause ' ' 800 continue ! ! continue for the next day of month 801 do idy=nndy,ndyend call blet2(idy,ADY) CALL subgec8o(AYEAR,AMN,ADY,gecdy) CALL submag14m1(AYEAR,AMN,idy) ! WU, WL, WE m.n,s enddo goto 33 ! end of calculations C c write (*,888) alat,alon,aday c 888 format(1X,' alat=',A3,' alon=',A3,' day=',A3) 2 CONTINUE if (irun.eq.0) then write(*,*) 'INPUT FILE IS NOT IN YOUR DIRECTORY ' pause ' ' goto 35 endif 33 CONTINUE 35 write(*,*) outfilef,outfileh,outfilet,outfilee write(*,*) irun C temp C> nndy=nndy+1 crem if (nndy.le.im(mn)) then C> if (nndy.le.ndyend) then C> goto 1001 C> endif pause ' ' STOP END C***********************************************************