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
