c iritest.for, version number can be found at the end of this comment. c----------------------------------------------------------------------- c c test program for the iri_web subroutine c c*********************************************************************** c!********************* IMPORTANT PLEASE READ ************************** c*********************************************************************** c c This test program shows how to initialize and call the IRI subroutine c IRI_SUB. Any program using IRI_SUB needs to include the following c statements (this iritest.for program shows you how and where to include c these statements): c c call read_ig_rz c call readapf107 c c The indices values are stored in the arrays ionoindx, indrz, aap, af107 c and are passed to other subroutines via common/igrz/ and common/apfa/ c c do i=1,100 c oar(i,1)=-1.0 c enddo c c If you want to use the standard version of IRI you should make sure c that the JF(50) logical values are set to the recommended default c values: c jf(4,5,6,23,30,33,35,39,40,47)=.false. all others =.true. c You can turn off (jf( )=.false.) the computation of certain parameters c if you do not need these parameters: c jf(1): Ne; jf(2): Te Ti Tn; jf(3): Ni; jf(21): ion drift; c jf(28): spread-F probability c For some parameters the default is already .false.: c jf(33): auroral boundaries; jf(35): foE storm model; c jf(47): CGM coordinates c If you use IRI with JF values other than the default values please c make sure to mention this in any publication that results from your c research. c c Required files: irisub.for, irifun.for, iritec.for, iridreg.for, c iriflip.for, igrf.for, cira.for, rocdrift.for c c Required i/o units: c IRISUB.FOR, IRIFUN.FOR, IGRF.FOR: c KONSOL= 6: Program messages to konsol (if jf(12)=.true.) c KONSOL=11: Program messages to MESSAGES.TXT (if jf(12)=.false.) c (COMMON/iounit/konsol,mess is used to pass the values from c IRISUB to IRIFUN and IGRF. If mess=false messages are turned off) c IRISUB: c IUCCIR=10: CCIR and URSI coefficients (CCIR%%.ASC, %%=month+10) c IRIFUN: c read_data_SD: UNIT=10 coefficients of Shubin (2015) hmF2 model c read_ig_rz: UNIT=12 Solar/ionospheric indices IG12, R12 (IG_RZ.DAT) c readapf107: UNIT=13 Magnetic indices and F10.7 (APF107.DAT) c IGRF: c GETSHC: UNIT=14 IGRF coeff. (DGRF%%%%.DAT, IGRF%%%%.DAT, %%%%=year) c c*********************************************************************** c*********************************************************************** c*********************************************************************** c c-version-mm/dd/yy ----------corrections-------------------------- c 2000.01 05/07/01 initial version c 2000.02 07/11/01 line 210: do i=1,100 instead of i=2,100 (K. Tokar) c 2000.03 28/12/01 output oar(39) for IG12 (R. Conkright, NGDC, NOAA) c 2000.04 28/10/02 replace TAB/6 blanks, enforce 72/line (D. Simpson) c 2000.05 02/06/03 Ne(Te) only 300,400; foF1 and hmF1 output corr. c 2000.06 01/19/05 (.not.jf(20)) instead of (..jf(2)) (G. Schiralli) c 2005.01 05/06/06 included spread-F (jf(28)) and topside (jf(29)) options C 2007.00 05/18/07 Release of IRI-2007 c 2007.02 10/31/08 outf(100) -> outf(500), numhei=numstp=500 c 2007.03 02/12/09 added new D-region option (h=-3) c 2007.11 04/19/10 correct TEC for normal output [Shunrong Zhang] c C 2012.00 10/05/11 IRI-2012: bottomside B0 B1 model (SHAMDB0D, SHAB1D), C 2012.00 10/05/11 bottomside Ni model (iriflip.for), auroral foE C 2012.00 10/05/11 storm model (storme_ap), Te with PF10.7 (elteik), C 2012.00 10/05/11 oval kp model (auroral_boundary), IGRF-11(igrf.for), C 2012.00 10/05/11 NRLMSIS00 (cira.for), CGM coordinates, F10.7 daily C 2012.00 10/05/11 81-day 365-day indices (apf107.dat), ap->kp (ckp), C 2012.00 10/05/11 array size change jf(50) outf(20,1000), oarr(100). C 2012.00 03/21/12 PIKTAB=4 output for D-region C 2012.01 09/16/12 Corrected UT output (hour-25) C 2012.03 09/18/14 input JF(18): FIELDG not UT_LT (A. Mazzella) C 2012.03 09/18/14 rzino, igino logical not real; and more (A. Mazzella) C 2012.03 09/18/14 jf defaults and special for piktab=3 (A. Mazzella) C 2012.03 09/23/14 jf(22): output option Ni in [m-3]/1.e9 C 2012.03 09/23/14 jf(36:38) include in input choices C 2012.03 09/23/14 added output options for parameters oar(59:86) C 2012.05 04/27/15 delete double line CHARACTER*9 pname(6) C 2015.01 06/30/15 scid=1.0E-9 (A. Charisi) C 2015.01 06/30/15 jino-outf(9)&jio2-outf(8) write.jio2,jino(A. Charisi) C 2015.01 07/12/15 output of ion density/composition now with 3 digits C 2015.01 07/12/15 call read_ig_rz readapf107 C 2015.02 08/13/15 delete COMMON/CONST2 C 2015.02 08/13/15 ursifo=jf(5) just before IRI_SUB call C 2016.01 06/01/16 User specified B0 when jf(43)=false C 2016.02 08/22/16 Allow user input options for F10.7D and F10.7_81 C 2016.03 04/25/18 cleaned up user input options: B0,B1,F10.7D,F10.7_81 C 2016.04 09/05/18 user input options: HNEA and HNEE C 2020.01 09/08/20 deleted extra jf(43) input (C. Vasyl) C 2020.02 09/08/20 re-arranged jf inputs - if user input no options C 2020.03 09/16/20 messages if auroral boundary or CGM are off C 2020.04 10/03/21 hmF2: Shubin model default, jf(40)=.false C 2020.05 10/22/23 user input for htec_min for iri_web C 2020.05 10/22/23 delete user input option for HNEA and HNEE C 2020.06 12/10/24 Ne-topside: COR2 default C 2020.06 12/10/24 Included output for Sporadic-E Occ. Prob. oar(91,*) C 2020.06 09/25/25 ******** released on irimodel.org ************* C C 2026.01 06/15/26 Change to Te-TBPS-2026 model C 2026.02 06/20/26 Added online/table text for Te,Ti,Es-prob,IBP C 2026.02 06/20/26 Added ESPROB, IBP to 6-parameter output option C INTEGER pad1(6),jdprof(77),piktab DIMENSION outf(20,1000),oar(100,1000),jfi(6) LOGICAL jf(50),rzino,igino CHARACTER*2 timev(2) CHARACTER*3 uni(92),sopt,seopt,ibpopt CHARACTER*4 IMZ(8),MAP,xtex,coorv(2),plpa CHARACTER*5 ITEXT(8),tsopt CHARACTER*6 pna(92) CHARACTER*7 popt CHARACTER*8 bopt,topt,tiopt,pplas CHARACTER*9 pname(7) CHARACTER*10 dopt,hopt CHARACTER*11 iopt,rzopt,igopt,f8opt,fdopt CHARACTER*16 f1opt CHARACTER*21 esopt DATA IMZ /' km ','GEOD','GEOD','yyyy',' mm ',' dd ','YEAR', & 'L.T.'/, ITEXT/' H ',' LATI', & ' LONG',' YEAR','MONTH',' DAY ','DAYOF',' HOUR'/ data pna/'NmF2','hmF2','NmF1','hmF1','NmE','hmE','NmD','hmD', & 'h05','B0','NVmin','hVtop','Tpeak','hTpek','T300','T400','T600', & 'T1400','T3000','T120','Ti450','hTeTi','sza','sundec','dip', & 'diplat','modip','Lati','Srise','Sset','season','Longi', & 'Rz12','cov','B1','M3000','TEC','TECtop','IG12','F1prb','F107d', & 'C1','daynr','vdrft','foF2r','F10781','foEr','sprd_F','MLAT', & 'MLON','Ap_t','Ap_d','invdip','MLTinv','CGMlat','CGMlon', & 'CGMmlt','CGM_AB','CGMm0','CGMm1','CGMm2','CGMm3','CGMm4', & 'CGMm5','CGMm6','CGMm7','CGMm8','CGMm9','CGMm10','CGMm11', & 'CGMm12','CGMm13','CGMm14','CGMm15','CGMm16','CGMm17','CGMm18', & 'CGMm19','CGMm20','CGMm21','CGMm22','CGMm23','kp_t','dec','L', & 'DIMO','SR300','SS300','HNEA','HNEE','spor_E',' IBP '/ data uni/'m-3','km','m-3','km','m-3','km','m-3','km','km','km', & 'm-3','km','K','km',7*'K','km',6*'deg',2*'h',' ','deg',4*' ', & 'm-2','%',5*' ','m/s',3*' ','0-1',2*'deg',2*' ','deg','h', & 2*'deg','h',25*'deg',' ','deg',' ','Gau','h','h','km','km', & 2*'0-1'/,timev/'LT','UT'/,coorv/'geog','geom'/ data jfi/8,9,13,14,15,16/ c COMMON/const2/icalls,montho,nmono,iyearo,idaynro,ursifo,rzino, c & igino,ut0 c c icalls=0 c montho=-1 c nmono=-1 c iyearo=-1 c idaynro=-1 c rzino=.true. c igino=.true. c ut0=-1 call read_ig_rz call readapf107 nummax=1000 do 6249 i=1,100 6249 oar(i,1)=-1.0 c user input of IRI input parameters c 1 print *,'jmag(=0/1,geog/geom),lati/deg,long/deg' read(5,*) jm,xlat,xlon print *,'year(yyyy),mmdd(or -ddd),iut(=0/1,LT/UT),hour' read(5,*) iy,imd,iut,hour print *,'height/km' read(5,*) hx print *,'variable? (1/2/../8 for height/lat/long/year/month/', & 'day/day of year/hour)' read(5,*) ivar print *,'begin, end, and stepsize for the selected variable' read(5,*) vbeg,vend,vstp print *,'output-option (if variable=height then choose 0, 3,', & '4, or 5)' print *,'(enter 0 for standard table of IRI parameters)' print *,'(enter 1 for list of peak heights and densities)' print *,'(enter 2 for plasma frequencies, B0, M3000, ', & 'valley width, and depth,)' c outf(15,*)= ratio plasma to gyro frequency print *,'(enter 3 for 6 parameters of your choice)' print *,'(enter 4 for D-region models at 60,65,..,110 km)' print *,'(enter 5 special test output)' read(5,*) piktab print *,'lower and upper height [km] for TEC integration', & ' (valid range: 65-30000 km) (0,0 for no TEC)' read(5,*) htec_min,htec_max print *,'Options (50): t(rue) or f(alse)' print *,'Enter 0 to use standard or 1 to enter your own' read(5,*) jchoice do i=1,50 jf(i)=.true. enddo if(piktab.eq.4) jf(24)=.false. if(jchoice.eq.0) then c defaults for jf(1:50) c jf(1)=.false. ! f=no electron densities (t) c jf(2)=.false. ! f=no temperatures (t) c jf(3)=.false. ! f=no ion composition (t) jf(4)=.false. ! t=B0table f=other models (f) jf(5)=.false. ! t=CCIR f=URSI foF2 model (f) jf(6)=.false. ! t=DS95+DY85 f=RBV10+TBT15 (f) c jf(7)=.false. ! t=tops f10.7<188 f=unlimited (t) c jf(8)=.false. ! t=foF2 f=user input (t) c jf(9)=.false. ! t=hmF2 from model f=user input (t) c jf(10)=.false. ! t=Te standard f=using Ne user input (t) c jf(11)=.false. ! t=Ne standard f=Lay-function (t) c jf(12)=.false. ! t=messages to unit 6 f=messages.txt on 11 (t) c jf(13)=.false. ! t=foF1 f=foF1 or NmF1 user input (t) c jf(14)=.false. ! t=hmF1 f=user input (only Lay version) (t) c jf(15)=.false. ! t=foE f= foE or NmE user input (t) c jf(16)=.false. ! t=hmE f=hmE user input (t) c jf(17)=.false. ! t=Rz12 from file f=user input (only Lay version) (t) c jf(18)=.false. ! t=IGRF dip,magbr,modip f=POGO68/10 for 1973(t) c jf(19)=.false. !F1 prob model only if foF1>0 and not NIGHT (t) c jf(20)=.false. !standard F1 standard F1 plus L condition (t) c (19,20) = (t,t) f1-prob, (t,f) f1-prob-L, (f,t) old F1, (f,f) no F1 c jf(21)=.false. ! t=ion drift computed f=not comp.(t) c jf(22)=.false. ! ion densities in m-3 (t) jf(23)=.false. ! t=Te-AEROS/ISIS f=Te-TBPS-2026 (f) c jf(24)=.false. ! t=D-reg-IRI-1990 f=FT-2001 (t) c jf(25)=.false. ! t=F107D from APF107.DAT f=user (t) c jf(26)=.false. ! t=STORM model on f= off (t) c jf(28)=.false. ! t=spread-F computed f=not comp. (f) c jf(29)=.false. ! t=old f=New Topside options (t) jf(30)=.false. ! t=corr f=NeQuick topside (f) c (29,30) = (t,t) IRIold, (f,t) IRIcor, (f,f) NeQuick, (t,f) COR2 c jf(31)=.false. ! t=B0ABT f=Gulyaeva (t) c jf(32)=.false. ! t=F107_81 from APF107.DAT f=user (t) jf(33)=.false. ! t=auroral boundary f=off (f) c jf(34)=.false. ! t=messages on f= off (t) jf(35)=.false. ! t=auroral E-storm model on f=off (f) c jf(36)=.false. ! t=hmF2 w/out foF2_storm f=with (t) c jf(37)=.false. ! t=topside w/out foF2_storm f=with (t) c jf(38)=.false. ! t=bubble prob. IBP on f=off (t) jf(39)=.false. ! t=M3000F2 model f=new hmF2 models (f) jf(40)=.false. ! t=AMTB-model, f=Shubin-COSMIC model (f) c jf(41)=.false. ! t:COV=F10.7_386 f:COV=f(IG12) (t) c jf(42)=.false. ! t/f=TBPS-2026 w/o PF10.7 dependance (t) c jf(43)=.false. ! t= B0 model f= B0 user input (t) c jf(44)=.false. ! t= B1 model f= B1 user input (t) C jf(45)=.false. ! t=Es occ. prob on f=off (t) C jf(46)=.false. ! t=Es prob w/o f=Es prob with F10.7 dep.(t) C (45,46) = (t,t) ESPROB, (t,f) ESPRO with F10.7, (f,f) ESPRO off, C (f,t) ESPROB with F10.7 and kp dependence jf(47)=.false. ! t=CGM on f=CGM off (f) c jf(48)=.false. ! t=Ozhogin plasmasphere f=Gallagher (t) c jf(49)=.false. ! t=Ti-Tru2021 f=Ti-Bil1981 (t) c jf(50)=.false. ! t=without plasmapause f=with plasmapause (t) else print *,'Compute Ne, T, Ni? (enter: t,t,t if you want all)' read(5,*) jf(1),jf(2),jf(3) if(jf(1)) then c print *,'Ne lower boundary: t=65/80km day/night, ', c & 'f=user input {t}' c read(5,*) jf(45) c print *,'Ne upper boundary: t=2000km day/night, ', c & 'f=user input {t}' c read(5,*) jf(46) print *,'Sporadic E Occ. Prob.: jf(45),jf(46)' print *,'(t,t) ESPROB on, (t,f) with F10.7 dep., ', & '(f,t)=with F10.7&ap dep., f,f=ESPROB off ' read(5,*) jf(45),jf(46) print *,'Plasma Bubble Prob. model IBP: t=on, f=off', & ' {standard:t}' read(5,*) jf(38) print *,'LAY version: t=standard ver., f=LAY version.', & ' {standard:t}' read(5,*) jf(11) print *,'Ne Topside: jf(29),jf(30) {t,f}' print *,'(t,t) IRI-2001, (t,f) COR2, ', & '(f,t)=IRIcor, f,f=NeQuick ' read(5,*) jf(29),jf(30) print *,'Ne Topside: t=F10.7<188, f=unlimited {t}' read(5,*) jf(7) print *,'F2 peak density or foF2: t=model, ', & 'f=user input {t}' read(5,*) jf(8) if(jf(8)) then print *,'foF2 model: t=CCIR-67, f=URSI-89 {f}' read(5,*) jf(5) print *,'foF2: t=with storm model, f=without {t}' read(5,*) jf(26) print *,'Ne Topside: t=w/o foF2 storm model, f=with {t}' read(5,*) jf(37) endif print *,'F2 peak height or M3000F2: t=model, ', & 'f=user input {t}' read(5,*) jf(9) if(jf(9)) then print *,'hmF2: t=f(M3000F2), f=new models {f}' read(5,*) jf(39) print *,'hmF2: t=AMTB-model, f=Shubin-COSMIC model {f}' read(5,*) jf(40) if(jf(39)) then print *,'hmF2: t=w/o foF2 storm model, f=with {t}' read(5,*) jf(36) endif endif print *,'B0 bottomside thickness: t=model, ', & 'f=user input {t}' read(5,*) jf(43) print *,'B1 bottomside shape: t=model, ', & 'f=user input {t}' read(5,*) jf(44) if(jf(43)) then print *,'Bottomside thickness B0: t=Bil-2000, ', & 'f=other options {f}.' read(5,*) jf(4) print *,'Bottomside thickness B0: t=ABT-2009, ', & 'f= Gul-1987 {t}.' read(5,*) jf(31) endif print *,'F1 peak density or foF1: t=model, ', & 'f=user input {t}' read(5,*) jf(13) if(.not.jf(11)) then print *,'F1 peak height: t=model, f=user input {t}' read(5,*) jf(14) endif print *,'F1: jf(19),jf(20) {t,t}' print *,'(t,t) stdard f1 prob, (t,f) f1-prob with L-cond' print *,'(f,t) only if foF1>0 and not NIGHT, (f,f) no F1' read(5,*) jf(19),jf(20) print *,'E peak density or foE: t=model, f=user input {t}' read(5,*) jf(15) print *,'E peak height: t=model, f=user input {t}' read(5,*) jf(16) print *,'E peak auroral storm model: t=on, f=off {f}' read(5,*) jf(35) print *,'D: t=IRI-1990, f= FT-2001 {t}' read(5,*) jf(24) endif if(jf(2)) then print *,'Te(Ne) model: t=not used, f=correlation is', & ' used. {t}' read(5,*) jf(10) print *,'Te: t=Bil-1985, f=TBPS-2026 {f}' read(5,*) jf(23) print *,'Te: t=TBPS-2026 with PF107 dep., f=w/o {t}' read(5,*) jf(42) print *,'Ti: t=Tru-2021, f=Bil-1981 {t}' read(5,*) jf(48) endif if(jf(3)) then print *,'Ion comp. model: t=DS95/DY85, f=RBV10/TBT15 {f}' read(5,*) jf(6) print *,'Ni: t=ion composition in %, f=ion densities', & 'in cm-3 {t}' read(5,*) jf(22) endif print *,'Equat. Vert. Ion Drift: t=computed, ', & 'f=not computed {t}' read(5,*) jf(21) print *,'Spread-F probability: t=computed, ', & 'f=not computed {t}' read(5,*) jf(28) print *,'Auroral boundary model: t=on, f=off {f}' read(5,*) jf(33) print *,'COV: t: COV=F10.7_365, f: COV=func(IG12). {t}' read(5,*) jf(41) print *,'Sunspot index: t=from file, f=user input. {t}' read(5,*) jf(17) print *,'Ionospheric index: t=from file, f=user input. {t}' read(5,*) jf(27) print *,'F10.7D Index: t=from file, f=user input {t}' read(5,*) jf(25) print *,'F10.7_81 Index: t=from file, f=user input {t}' read(5,*) jf(32) print *,'dip, magbr, modip: t=IGRF, f=old FIELDG using ', & 'POGO68/10 for 1973 {t}' read(5,*) jf(18) print *,'CGM computation: t=on, f=off {f}' read(5,*) jf(47) print *,'Messages on (t) off (f) {t}' read(5,*) jf(34) print *,'Message output unit: t=(UNIT=6), f=(UNIT=11). {t}' read(5,*) jf(12) endif c if(piktab.gt.3) jf(24)=.false. c option to enter six additional parameters c if(PIKTAB.eq.3) then print *,'6 Parameters of your choice (number:1-91)' print *,(pna(j),j=1,10) print *,(pna(j),j=11,20) print *,(pna(j),j=21,30) print *,(pna(j),j=31,40) print *,(pna(j),j=41,50) print *,(pna(j),j=51,60) print *,(pna(j),j=61,70) print *,(pna(j),j=71,80) print *,(pna(j),j=81,90) print *,(pna(91)) print *,'e.g. 48,44,45,91,92,58 for:' print *,' spread-F probability [48]' print *,' equatorial vertical ion drift [44]' print *,' foF2_storm/foF2_quiet [45]' print *,' sporadic-E probability [91]' print *,' plasma bubble probability [92]' print *,' CGM-Lat of equatorw. auroral boundary [58]' print *,'or 23,27,58,85,51,53 for:' print *,' solar zenith angle [23]' print *,' modified dip latitude [27]' print *,' L-value [85]' print *,' Ap for current UT [51]' print *,' magnetic local time (MLT) [54]' print *,' invdip_old [53]' c change defaults for computation of specific parameters jf(21)=.true. ! spread-F prob. computed jf(28)=.true. ! vertical ion drift computed jf(33)=.true. ! auroral boundary computed jf(45)=.true. ! sporadic-E prob. computed read(5,*) (pad1(jjj),jjj=1,6) do jjj=1,6 if(pad1(jjj).eq.44) JF(21)=.true. if(pad1(jjj).eq.45) JF(26)=.true. if(pad1(jjj).eq.47) JF(35)=.true. if(pad1(jjj).eq.48) JF(28)=.true. if(pad1(jjj).eq.55) JF(47)=.true. if(pad1(jjj).eq.58) JF(33)=.true. if(pad1(jjj).eq.91) JF(45)=.true. if(pad1(jjj).eq.92) JF(38)=.true. enddo endif c option to enter measured values for NmF2, hmF2, NmF1, hmF1, NmE, hmE, c B0, N(300), N(400), N(600) if available; c print *,' ' print *,' ' print *,' ' numstp=int((vend-vbeg)/vstp)+1 if(ivar.eq.1) numstp=1 if(jf(1)) then if(.not.jf(8).or..not.jf(9).or..not.jf(13).or..not.jf(14).or. & .not.jf(15).or..not.jf(16).or..not.jf(43).or..not.jf(44)) then var=vbeg i=1 c --------- user input: foF2 or NmF2 2234 if(.not.jf(8)) then jf(26)=.false. ! storm model off, if user input print *,'foF2/Mhz or NmF2/m-3 for ',itext(ivar), & '=',var read(5,*) oar(1,i) pname(1)='foF2/MHz' if(oar(1,i).gt.30.) pname(1)='NmF2/m-3' endif c --------- user input: hmf2 or M(3000)F2 if(.not.jf(9)) then print *,'hmF2/km or M3000F2 for ',itext(ivar), & '=',var read(5,*) oar(2,i) pname(2)='M(3000)F2' if(oar(2,i).gt.50.) pname(2)='hmF2/km' endif c --------- user input: foF1 or NmF1 if(.not.jf(13)) then print *,'foF1/MHz or NmF1/m-3 for ',itext(ivar), & '=',var read(5,*) oar(3,i) pname(3)='foF1/MHz' if(oar(3,i).gt.30.) pname(3)='NmF1/m-3' endif c --------- user input: hmF1 if(.not.jf(14)) then print *,'hmF1/km for ',itext(ivar),'=',var read(5,*) oar(4,i) pname(4)='hmF1/km' endif c --------- user input: foE or NmE if(.not.jf(15)) then print *,'foE/MHz or NmE/m-3 for ',itext(ivar), & '=',var read(5,*) oar(5,i) pname(5)='foE/MHz' if(oar(5,i).gt.30.) pname(5)='NmE/m-3' endif c --------- user input: hmE if(.not.jf(16)) then print *,'hmE/km for ',itext(ivar),'=',var read(5,*) oar(6,i) pname(6)='hmE/km' endif c --------- user input: B0 if(.not.jf(43)) then print *,'B0/km for ',itext(ivar),'=',var read(5,*) oar(10,i) pname(7)='B0/km ' endif c --------- user input: B1 if(.not.jf(44)) then print *,'B1 for ',itext(ivar),'=',var read(5,*) oar(35,i) pname(7)='B1 ' endif i=i+1 var=var+vstp if(ivar.gt.1.and.var.le.vend) goto 2234 endif endif c option to enter Ne for Te-Ne relationship c if(jf(2).and..not.jf(10)) then var=vbeg do 1235 i=1,numstp print *,'Ne(300km),Ne(400km)/m-3', & ' for ',itext(ivar),'=',var,' [-1 if not]' read(5,*) oar(15,i),oar(16,i) 1235 var=var+vstp endif c option to enter F107D and/or F107_81 c if(.not.jf(25)) then print *,'User input for F10.7D' read(5,*) oar(41,1) do i=2,100 oar(41,i)=oar(41,1) enddo endif if(.not.jf(32)) then print *,'User input for F10.7_81avg' read(5,*) oar(46,1) do i=2,100 oar(46,i)=oar(46,1) enddo endif c option to enter Rz12 and/or IG12 c if(.not.jf(17)) then print *,'User input for Rz12' read(5,*) oar(33,1) do i=2,100 oar(33,i)=oar(33,1) enddo endif if(.not.jf(27)) then print *,'User input for IG12' read(5,*) oar(39,1) do i=2,100 oar(39,i)=oar(39,1) enddo endif c end of user input c num1=(vend-vbeg)/vstp+1 numstp=iabs(num1) if(numstp.GT.nummax) numstp=nummax pplas='Har-2023' if(jf(49)) pplas='Ozh-2012' plpa='with' if(jf(50)) plpa='w/o ' if(jf(29)) then if(jf(30)) then popt='IRI2001' else popt='IRIcor2' endif else if(jf(30)) then popt='IRIcorr' else popt='NeQuick' endif endif map='URSI' if(jf(5)) map='CCIR' if(jf(39)) then hopt='CCIR-M3000' else if(jf(40)) then hopt='AMTB-2013' else hopt='Shubin2015' endif endif if(jf(4)) then bopt='BIl-2000' else if(jf(31)) then bopt='ABT-2009' else bopt='Gul-1987' endif endif if(jf(45)) then if(jf(46)) then esopt='without F10.7,ap dep.' else esopt='with F10.7 dependence' endif else if(jf(46)) then esopt='with F10.7 and ap dep.' else esopt='model turned off ' endif endif ibpopt='off' if(jf(38)) ibpopt='on ' iopt='RBV10+TBT15' if(jf(6)) iopt='DS95 + DY85' dopt='FT01+DRS95' if(jf(24)) dopt='IRI-1990' sopt='off' if(jf(26)) sopt='on ' seopt='off' if(jf(35)) seopt='on ' topt='TBPS2026' if(jf(23)) topt='Bil-1985' tsopt=' with' if(jf(23)) tsopt=' w/o' tiopt='Bil-1981' if(jf(48)) tiopt='TBKS2021' if(jf(19)) then f1opt='Scotto-97 no L' if(.not.jf(20)) f1opt='Scotto-97 with L' else f1opt='IRI-95' if(.not.jf(20)) f1opt='no F1 region' endif rzopt=' user input' if(jf(17)) rzopt=' ' igopt=' user input' if(jf(27)) igopt=' ' fdopt=' user input' if(jf(25)) fdopt=' ' f8opt=' user input' if(jf(32)) f8opt=' ' hxx=hx jmag=jm mmdd=imd c calling IRI subroutine c phour=hour call iri_web(jmag,jf,xlat,xlon,iy,mmdd,iut,hour, & hxx,htec_min,htec_max,ivar,vbeg,vend,vstp,outf,oar) c preparation of results page c write(7,3991) iy,mmdd,phour,timev(iut+1), & coorv(jmag+1),xlat,xlon,hxx if(jf(1)) then write(7,3314) popt write(7,4314) pplas,plpa if(jf(8)) then write(7,301) map write(7,3291) sopt endif if(jf(9)) write(7,303) hopt write(7,309) bopt write(7,3295) f1opt write(7,3299) seopt write(7,3081) dopt numi=numstp if(ivar.eq.1) numi=1 do j=1,6 ij=jfi(j) if(.not.jf(ij)) then write(7,302) pname(j) write(7,402) (oar(j,i),i=1,numi) endif enddo if(.not.jf(43)) then write(7,302) pname(7) write(7,402) (oar(10,i),i=1,numi) endif endif write(7,3277) esopt write(7,3278) ibpopt if(jf(2)) then write(7,3292) topt,tsopt write(7,3297) tiopt endif if(jf(3)) write(7,329) iopt if(.not.jf(47)) write(7,7234) if(.not.jf(33)) write(7,7235) if(ivar.eq.1) then if(oar(3,1).lt.1.) oar(4,1)=0. yp2=0 if(oar(3,1).gt.0.0) yp2=oar(3,1)/1.e6 write(7,213) oar(1,1)/1.E6,yp2,oar(5,1)/1.E6 write(7,214) oar(2,1),oar(4,1),oar(6,1) else write(7,307) endif write(7,211) oar(23,1),oar(25,1),oar(27,1) write(7,223) oar(33,1),rzopt write(7,2231) oar(39,1),igopt write(7,2238) oar(41,1),fdopt write(7,2237) oar(46,1),f8opt if(htec_max.gt.50.0) write(7,3914) htec_min,htec_max 3991 format(///'yyyy/mmdd(or -ddd)/hh.h):',I4,'/',I4,'/',F4.1, & A2,2X,A4,' Lat/Long/Alt=',F5.1,'/',F6.1,'/',F6.1/) 3914 format(/'TEC [1.E16 m-2] is obtained by numerical integ', & 'ration in 1km steps'/' from ',F7.1,' to ',F7.1,' km. ', & 't is the percentage of TEC above the F peak.') 3916 format(/'M3000F2: Propagation factor related to hmF2'/ & 'B0: bottomside thickness parameter.') 301 format(A4,' maps are used for the F2 peak density (NmF2)') 302 format(A9,' provided by user:') 402 format(7(1PE10.3)) 303 format(A10,' model is used for F2 peak height (hmF2)') 304 format(A25,'=',F5.1,') provided by user') 307 format(1x/'Solar and magnetic parameter for the 1st profile', & ' point:') 309 format(A8,' option is used for the bottomside thickness ', & 'parameter B0') 3081 format(A10,' option is used for D-region') 329 format(A11,' option is used for ion composition') 3291 format('foF2 STORM model is turned ',A3) 3299 format('foE auroral storm model is turned ',A3) 3277 format('sporadic-E occ. prob. ',A21) 3278 format('bubble probability model (IBP) ',A3) 3292 format(A8,A5,' solar dependence is used for the electron', & ' temperature') 3297 format(A8,' option is used for the ion temperature') 3293 format(A8,' option is used for the D-region Ne') 3295 format(A16,' option is used for the F1 occurrence probability') 211 format('Solar Zenith Angle/degree',28X,F6.1/ & 'Dip (Magnetic Inclination)/degree',21X, & F6.2/'Modip (Modified Dip)/degree',27X,F6.2) 223 format('Solar Sunspot Number (12-months running mean) Rz12', & 4X,F5.1,A11) 2231 format('Ionospheric-Effective Solar Index IG12',16X, & F5.1,A11) 2238 format('Solar radio flux F10.7 (daily)',24X, & F5.1,A11) 2237 format('Solar radio flux F10.7 (81-day average)',15X, & F5.1,A11) 213 format(/'Peak Densities/cm-3: NmF2=',F9.1,' NmF1=',F9.1, & ' NmE=',F9.1) 214 format('Peak Heights/km: hmF2=',F9.2,' hmF1=',F9.2, & ' hmE=',F9.2/) 3314 format(A7,' is used for topside Ne profile') 7234 format('CGM coordinate computation is turned off') 7235 format('Auroral boundary computation is turned off') 4314 format(A8,' plasmasphere model ',A4,' plasmapause') c c table head ....................................................... c agnr=7 !output unit number xtex=imz(ivar) if(jmag.gt.0.and.(ivar.eq.2.or.ivar.eq.3)) xtex='GEOM' if(iut.gt.0.and.ivar.eq.8) xtex='U.T.' IF(PIKTAB.EQ.4) WRITE(7,8199) IF(PIKTAB.EQ.3) WRITE(7,8191) ITEXT(IVAR), & (pna(pad1(j)),j=1,6),xtex,(uni(pad1(j)),j=1,6) IF(PIKTAB.EQ.2) WRITE(7,8194) ITEXT(IVAR),xtex IF(PIKTAB.EQ.1) WRITE(7,8192) ITEXT(IVAR),xtex IF(PIKTAB.EQ.0) THEN if(jf(22)) then WRITE(7,8193) ITEXT(IVAR),xtex else WRITE(7,9193) ITEXT(IVAR),xtex endif ENDIF 8191 FORMAT(/'-'/2X,A5,6A10/3X,A4,6A10) 8192 FORMAT(/'-'/2X,A5,6X,'PEAK ALTITUDES IN KM',8X,'PEAK DEN', & 'SITIES IN cm-3 TEC top/%'/3X,A4,' hmF2 hmF1 hmE ', & 'hmD NmF2 NmF1 NmE NmD 1E16m-2') 8194 FORMAT(/'-'/2X,A5,3X,'M3000 B0',3X,'B1 E-VALLEY',7X, & 'PLASMA ','FREQUENCIES / MHz'/3X,A4,11X,'km W/km ', & ' Depth',5X,'foF2 foF1 foE foD') 8193 FORMAT(/'-'/1X,A5,' ELECTRON DENSITY TEMPERATURES ', & 5X,'ION PERCENTAGES[%]*10',4X,'1E16m-2'/2X,A4,' Ne/cm-3 Ne/', & 'NmF2 Tn/K Ti/K Te/K O+ N+ H+ He+ O2+ NO+ Clust TEC t/%') 9193 FORMAT(/'-'/1X,A5,' ELECTRON DENSITY TEMPERATURES ', & 4X,'ION DENSITIES[cm-3]/100',4x,'1E16m-2'/2X,A4,' Ne/cm-3 Ne/', & 'NmF2 Tn/K Ti/K Te/K O+ N+ H+ He+ O2+ NO+ Clust TEC t/%') 8199 FORMAT(/'-'/1X,'h',8X,' D-REGION ELECTRON DENSITY IN CM-3'/ & 1X,'km',18X,'DRS-95: Stratos Warming/Winter Anomaly'/5X, & 'IRI-07',4x,'FIRI SW/WA=0/0 0.5/0 1/0 0/0.5 0/1') c c output: D-region PIKTAB=4 c c D-REGION ELECTRON DENSITY IN CM-3: c IRI-07 FIRI Danilov:SW/WA=0/0 0.5/0 1/0 0/0.5 0/1 c DRS-95: Stratos Warming/Winter Anomaly c if(piktab.eq.4) then do 2591 lix=1,77 jdprof(lix)=-1 dichte=outf(14,lix) 2591 if(dichte.gt.0.) jdprof(lix)=int(dichte/1.e6+0.5) do 2592 lix=1,11 ihtemp=55+lix*5 WRITE(7,3810) ihtemp,jdprof(lix),jdprof(lix+11), & jdprof(lix+22),jdprof(lix+33),jdprof(lix+44), & jdprof(lix+55),jdprof(lix+66) 2592 continue 3810 FORMAT(I3,7I8) goto 2357 endif xcor=vbeg do 1234 li=1,numstp c c output: peak densities and altitudes PIKTAB=1 c IF(PIKTAB.eq.1) THEN if(oar(3,li).lt.1.) oar(4,li)=0. iyp1=int(oar(1,li)/1.e6+.5) iyp2=0 if(oar(3,li).gt.0.0) iyp2=int(oar(3,li)/1.e6+.5) iyp3=int(oar(5,li)/1.e6+.5) iyp4=int(oar(7,li)/1.e6+.5) tec=oar(37,li) if(tec.gt.0.0) then tec=tec/1.e16 itopp=int(oar(38,li)+.5) else tec=-1.0 itopp=-1 endif WRITE(7,3910) XCOR,oar(2,li),oar(4,li),oar(6,li),oar(8,li), & iyp1,iyp2,iyp3,iyp4,tec,itopp 3910 FORMAT(F7.1,2X,4F6.1,1X,I9,3I7,1X,F6.2,I4) GOTO 1234 ENDIF c c output: plasma frequencies and profile parameters PIKTAB=2 c IF(PIKTAB.eq.2) THEN if(oar(3,li).lt.1.) oar(4,li)=0. fyp1=SQRT(oar(1,li)/1.24E10) fyp2=0 if(oar(3,li).gt.0.0) fyp2=SQRT(oar(3,li)/1.24E10) fyp3=SQRT(oar(5,li)/1.24E10) fyp4=SQRT(oar(7,li)/1.24E10) tec=oar(37,li) if(tec.gt.0.0) then tec=tec/1.e16 itopp=int(oar(38,li)+.5) else tec=-1.0 itopp=-1 endif wvalley=oar(12,li)-oar(6,li) dvalley=0.0 if(oar(5,li).gt.0.0) dvalley=oar(11,li)/oar(5,li) WRITE(7,3950) XCOR,oar(36,li),oar(10,li),oar(35,li), & wvalley,dvalley,fyp1,fyp2,fyp3,fyp4 3950 FORMAT(F7.1,2X,F6.4,F6.1,F4.1,F6.1,F8.4,1X,4F7.3) GOTO 1234 ENDIF c c output: 6 parameters of your choice PIKTAB=3 c IF(PIKTAB.eq.3) THEN c if(pad1.eq.45.and.oar(pad1,li).le.0.0) oar(pad1,li)=-1. c if(pad2.eq.45.and.oar(pad2,li).le.0.0) oar(pad2,li)=-1. c if(pad3.eq.45.and.oar(pad3,li).le.0.0) oar(pad3,li)=-1. WRITE(7,3919) XCOR,oar(pad1(1),li),oar(pad1(2),li), & oar(pad1(3),li),oar(pad1(4),li),oar(pad1(5),li), & oar(pad1(6),li) 3919 FORMAT(F7.1,6(1X,1PE9.2)) GOTO 1234 ENDIF c c output: special for test purposes PIKTAB=5 c IF(PIKTAB.eq.5) THEN c ----------- B0, B1 ---------------- c WRITE(8,4919) XCOR,jf(4),jf(31),jf(43),oar(10,li),oar(35,li) c4919 FORMAT(F7.1,2X,F6.2,2X,F5.3) c print*,iy,imd,hour,xlat,xlon,jf(4),jf(31),jf(43),oar(10,li) c ----------- Tn, Ti, Te ---------------- c print*,XCOR,outf(2,li),outf(3,li),outf(4,li) WRITE(8,4919) XCOR,outf(2,li),outf(3,li),outf(4,li) 4919 FORMAT(F7.1,2X,F6.1,2X,F6.1,2X,F6.1) c ----------- SZA, Ne ---------------- c WRITE(8,4919) oar(23,li),outf(1,li) c4919 FORMAT(F7.1,2X,E12.5) c ----------- Ne, TEC ---------------- c WRITE(8,4919) XCOR,outf(1,li),oar(37,li) c4919 FORMAT(F7.1,2X,E12.5,2X,E12.5) c ----------- Ni ---------------- c type*,XCOR,outf(1,li),outf(5,li),outf(6,li), c & outf(7,li),outf(8,li),outf(9,li),outf(10,li),outf(11,li) c WRITE(8,4919) XCOR,outf(1,li),outf(5,li),outf(6,li), c & outf(7,li),outf(8,li),outf(9,li),outf(10,li),outf(11,li) c4919 FORMAT(F7.1,2X,E12.5,7F10.4) c ----------- ion drift ---------------- c type*,xcor,oar(44,li),oar(41,li),oar(43,li) c WRITE(8,4919) XCOR,oar(44,li) c4919 FORMAT(F7.1,1X,F6.2) c ----------- hmF2 ---------------- c type*,XCOR,oar(26,li),oar(2,li),oar(36,li), c & sqrt(oar(1,li)/oar(5,li)),oar(33,li),oar(39,li) c WRITE(8,4919) XCOR,outf(1,li),outf(5,li),outf(6,li), c & outf(7,li),outf(8,li),outf(9,li),outf(10,li),outf(11,li) c4919 FORMAT(F7.1,1X,E12.5,7F10.4) c ----------- auroral boundary ---------------- c print *,XCOR,oar(55,li),oar(56,li),oar(54,li), c & oar(57,li),oar(54,li)-oar(57,li),oar(58,li) cc WRITE(8,4919) XCOR,outf(1,li),outf(5,li),outf(6,li), cc & outf(7,li),outf(8,li),outf(9,li),outf(10,li),outf(11,li) cc4919 FORMAT(F7.1,1X,E12.5,7F10.4) c ----------- Sporadic E occurrence probability --------- c WRITE(8,4919) XCOR,oar(91,li) c4919 FORMAT(F7.1,2X,F6.1) c ----------- Spread F occurrence probability --------- c WRITE(8,4919) XCOR,oar(48,li) c4919 FORMAT(F7.1,2X,F6.1) c ----------- CGM Coordinates --------- c WRITE(8,4919) XCOR,(oar(jk,li),jk=54,62), c & oar(54,li)-oar(57,li) c4919 FORMAT(F7.1,2X,9F6.1) c ----------- NmF2, hmF2, vi, spread_F, MLAT, Ap_daily --------- c WRITE(8,4919) XCOR,oar(1,li),oar(2,li),oar(44,li), c & oar(48,li),oar(49,li), oar(52,li) c4919 FORMAT(F7.1,2X,1PE8.2,2X,1PE8.2,2X,1PE8.2,2X,1PE8.2,2X, c & 1PE8.2) GOTO 1234 ENDIF c c output: standard c if(ivar.eq.1) then oar(1,li)=oar(1,1) oar(37,li)=oar(37,1) oar(38,li)=oar(38,1) endif jne=int(outf(1,li)/1.e6+.5) xner=outf(1,li)/oar(1,li) jtn=int(outf(2,li)+.5) jti=int(outf(3,li)+.5) jte=int(outf(4,li)+.5) scid=1.0E-8 if(jf(22)) scid=10. jio=INT(OUTF(5,li)*scid+.5) jih=INT(OUTF(6,li)*scid+.5) jihe=INT(OUTF(7,li)*scid+.5) jio2=INT(OUTF(8,li)*scid+.5) jino=INT(OUTF(9,li)*scid+.5) jicl=INT(OUTF(10,li)*scid+.5) jin=INT(OUTF(11,li)*scid+.5) c print *,'O+,N+,O2+,NO+=',outf(5,li), c & outf(11,li),outf(9,li),outf(8,li) if(outf(1,li).lt.0) jne=-1 if(outf(1,li).lt.0) xner=-1. if(outf(2,li).lt.0) jtn=-1 if(outf(3,li).lt.0) jti=-1 if(outf(4,li).lt.0) jte=-1 if(outf(5,li).lt.0) jio=-1 if(outf(6,li).lt.0) jih=-1 if(outf(7,li).lt.0) jihe=-1 if(outf(8,li).lt.0) jio2=-1 if(outf(9,li).lt.0) jino=-1 if(outf(10,li).lt.0) jicl=-1 if(outf(11,li).lt.0) jin=-1 tec=oar(37,li) if(tec.gt.0.0) then tec=tec/1.e16 itopp=int(oar(38,li)+.5) else tec=-1.0 itopp=-1 endif c print *, XCOR,jne,xner,jtn,jti,jte,jio,jin, c & jih,jihe,jino,jio2,jicl,tec,itopp WRITE(7,7117) XCOR,jne,xner,jtn,jti,jte,jio,jin, & jih,jihe,jio2,jino,jicl,tec,itopp 7117 FORMAT(F7.1,I7,1x,F6.3,3I6,7I4,f6.1,i4) 1234 xcor=xcor+vstp 2357 print *,'Enter 0 to exit or 1 to generate another profile?' read(5,*) icontinue if (icontinue.gt.0) goto 1 c print *,oar(51,1),oar(52,1),oar(83,1) stop end