c irisub.for, version number can be found at the end of this comment.
c-----------------------------------------------------------------------        
C Includes subroutine IRI_SUB to compute IRI parameters for specified
C location, date, time, and altitude range (Ne: 60/80km day/night to
C 30,000km; Te,Ti:60 - 3,000km; Ni: 75-2,000km) and subroutine
C IRI_WEB to computes IRI parameters for specified location, date, time 
C and variable range; variable can be altitude, latitude, longitude, 
C year, month, day of month, day of year, or hour (UT or LT). 
C IRI_WEB requires IRI_SUB. Both subroutines require linking with the 
c following files: IRIFUN.FOR, IRITEC.FOR, IRIDREG.FOR, 
c IRIFLIP.FOR CIRA.FOR, IGRF.FOR
c
c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
c!!!!!!!!!!!!!!!!!!!!!!!! IMPORTANT !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!        
c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
c
c Programs using subroutine IRI_SUB need to include (see IRITEST.FOR):
c
c	call read_ig_rz
c       call readapf107
c
c Programs using subroutine IRI_WEB need to include (see IRITEST.FOR):
c
c       do i=1,100
c          oar(i,1)=-1.0
c          enddo
c
c Please also make sure to use the default setting for the switches jf
c as noted in this comment section further down:
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 parmeters 
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 i/o units:  
c  KONSOL= 6 IRISUB: Program messages (used when jf(12)=.true. -> konsol)
c  IUCCIR=10 IRISUB: CCIR and URSI coefficients (CCIR%%.ASC, %%=month+10)
c  KONSOL=11 IRISUB: Program messages (used when jf(12)=.false. -> MESSAGES.TXT)
c  KONSOL=6/11 is also used in IRIFUN and IGRF. 
C  COMMON/iounit/konsol,mess 
c    is used to pass the value of KONSOL. If mess=false messages are turned off.
c  UNIT=12 IRIFUN/TCON:  Solar/ionospheric indices IG12, R12 (IG_RZ.DAT) 
c  UNIT=13 IRIFUN/APF..: Magnetic indices and F10.7 (APF107.DAT 
c  UNIT=14 IGRF/GETSHC:  IGRF coeff. (DGRF%%%%.DAT or IGRF%%%%.DAT, %%%%=year)
c  UNIT=15 IRIFUN/read_data_SD: coefficients of Shubin (2015) hmF2 model  
c  UNIT=16 IRIFUN/ibplod: coefficients of Stolle et al. (2023) IBP model  
c
c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
C
C CHANGES FROM  IRIS11.FOR  TO   IRIS12.FOR:
C    - CIRA-1986 INSTEAD OF CIRA-1972 FOR NEUTRAL TEMPERATURE
C    - 10/30/91 VNER FOR NIGHTTIME LAY-VERSION:  ABS(..)
C    - 10/30/91 XNE(..) IN CASE OF LAY-VERSION
C    - 10/30/91 CHANGE SSIN=F/T TO IIQU=0,1,2
C    - 10/30/91 Te > Ti > Tn ENFORCED IN FINAL PROFILE
C    - 10/30/91 SUB ALL NAMES WITH 6 OR MORE CHARACTERS
C    - 10/31/91 CORRECTED HF1 IN HST SEARCH:  NE(HF1)>NME
C    - 11/14/91 C1=0 IF NO F1-REGION
C    - 11/14/91 CORRECTED HHMIN AND HZ FOR LIN. APP.
C    -  1/28/92 RZ12=0 included
C    -  1/29/92 NEQV instead of NE between URSIF2 and URSIFO
C    -  5/ 1/92 CCIR and URSI input as in IRID12
C    -  9/ 2/94 Decimal month (ZMONTH) for IONCOM
C    -  9/ 2/94 Replace B0POL with B0_TAB; better annually
C    -  1/ 4/95 DY for h>hmF2
C    -  2/ 2/95 IG for foF2, topside; RZ for hmF2, B0_TAB, foF1, NmD
C    -  2/ 2/95 winter no longer exclusive for F1 occurrrence
C    -  2/ 2/95 RZ and IG incl as DATA statement; smooth annual var.
C
C CHANGES FROM  IRIS12.FOR  TO   IRIS13.FOR:
C    - 10/26/95 incl year as input and corrected MODA; nrm for zmonth
C    - 10/26/95 use TCON and month-month interpolation in foF2, hmF2
C    - 10/26/95 TCON only if date changes
C    - 11/25/95 take out logicals TOPSI, BOTTO, and BELOWE
C    - 12/ 1/95 UT_LT for (date-)correct UT<->LT conversion
C    - 12/22/95 Change ZETA cov term to cov < 180; use cov inst covsat
C    -  2/23/96 take covmax(R<150) for topside; lyear,.. for lt
C    -  3/26/96 topside: 94.5/BETA inst 94.45/..; cov -> covsat(<=188)
C    -  5/01/96 No longer DY for h>hmF2 (because of discontinuity)
C    - 12/01/96 IRIV13: HOUR for IVAR=1 (height)
C    -  4/25/97 D-region: XKK le 10 with D1 calc accordingly.
C    -  1/12/97 DS model for lower ion compoistion DY model
C    -  5/19/98 seamon=zmonth if lati>0; zmonth= ...(1.0*iday)/..
C    -  5/19/98 DY ion composition model below 300 km now DS model
C    -  5/19/98 DS model includes N+, Cl down to 75 km HNIA changed
C    -  5/28/98 User input for Rz12, foF1/NmF1, hmF1, foE/NmE, hmE
C    -  9/ 2/98 1 instead of 0 in MODA after UT_LT call
C    -  4/30/99 constants moved from DATA statement into program
C    -  4/30/99 changed konsol-unit to 13 (12 is for IG_RZ).
C    -  5/29/99 the limit for IG comp. from Rz12-input is 174 not 274
C    - 11/08/99 jf(18)=t simple UT to LT conversion, otherwise UT_LT
C    - 11/09/99 added COMMON/const1/humr,dumr also for CIRA86
C
C CHANGES FROM  IRIS13.FOR  TO   IRISUB.FOR:
C-Version-MM/DD/YY-Description (person reporting correction)
C 2000.01 05/09/00 B0_98 replaces B0_TAB and B1: 1.9/day to 2.6/night
C 2000.02 06/11/00 including new F1 and indermediate region
C 2000.03 10/15/00 include Scherliess-Fejer drift model
C 2000.04 10/29/00 include special option for D region models
C 2000.05 12/07/00 change name IRIS13 to IRISUB
C 2000.06 12/14/00 jf(30),outf(20,100),oarr(50)
C 2000.07 03/17/01 include Truhlik-Triskova Te model and IGRF
C 2000.08 05/07/01 include Fuller-Rowell-Condrescu storm model 
C 2000.09 07/09/01 LATI instead of LAT1 in F00 call -------- M. Torkar
C 2000.10 07/09/01 sdte instead of dte in ELTEIK call --- P. Wilkinson
C 2000.11 09/18/01 correct computation of foF2 for Rz12 user input
C 2000.12 09/19/01 Call APF only if different date and time -- P. Webb
c 2000.13 10/28/02 replace TAB/6 blanks, enforce 72/line -- D. Simpson
C 2000.14 11/08/02 change unit for message file to 11 (13 is Kp)
C 2000.15 01/27/03 change F1_prob output; Te-IK for fix h and ELTE(h)
C 2000.16 02/04/03 along<0 -> along=along+360; F1 occ for hmf1&foF1
C 2000.17 02/05/03 zyear =12.97 (Dec 31); idayy=#days per year
C 2000.18 02/06/03 jf(27) for IG12 user input; all F1 prob in oar
C 2000.19 07/14/04 covsat<188 instead of covsat=<f(IG)<188
C 2000.19 02/09/05 declare INVDIP as real ------------------ F. Morgan
C 2000.20 11/09/05 replace B0B1 with BCOEF --------------- T. Gulyaeva
C 2005.01 11/09/05 new topside ion composition; F107D from file 
C 2005.02 11/14/05 jf(18)=T: dip,mlat IGRF10 (igrf_dip igrf.for); F:POGO-75;
C 2005.03 11/15/05 sunrise/sunset/night for D,E,F1,F2; UT_LT removed
C 2005.04 05/06/06 FIRI D-region option not tied to peak
C 2005.04 05/06/06 Spread-F included, NeQuick included
C 2005.05 01/15/07 NeQuick uses CCIR-M3000F2 even if user-hmF2  
C 2007.00 05/18/07 Release of IRI-2007
C 2007.01 01/23/08 ryear = .. (daynr-1.0)/idayy ---------- R. Scharroo
C 2007.02 10/31/08 outf(100) -> outf(500), numhei=numstp=500 
C 2007.03 02/12/09 Jf(24)=.false.-> outf(1,60-140km)=FIRI- M. Friedrich
C 2007.04 03/14/09 SOCO(70->80;500->300km) --------------- R. Davidson
C 2007.05 03/26/09 call for APF_ONLY includes F107M
C 2007.09 08/17/09 STROM off if input; fof2in, fof1in,foein corr
C 2007.10 02/03/10 F10.7D = F10.7M = COV if EOF
C 2007.11 04/19/10 Corrections in irifun.for, cira.for 
C 2007.12 11/23/10 FNIGHT computed twice at 8334 --------- C. Vasly 
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.01 11/01/11 delete TEDER from EXTERNAL; GTD7 call 0 to 0.0
C 2012.01 12/12/11 put FMODIP in EXTERNAL; cgn_lon -> cgm_lon
C 2012.01 01/24/12 Change FLAT to LATI in SHAB1D call [D. Altadill]
C 2012.01 08/09/12 add jf(36)=t/f foF2 for hmF2 wout/with storm
C 2012.01 08/09/12 replace foF2_storm with foF2 for topside (NeQ, corr)
C 2012.01 08/09/12 call stormE_ap only if ap available
C 2012.01 08/09/12 If ap not available then auroral boundary for Kp=3
C 2012.02 12/17/12 Add magnetic declination as oarr(84) output
C 2012.03 02/13/13 Move B1 before B0 for Gulyaeva-1987
C 2012.03 02/20/13 Use foot-point for CGM to be closer to AACGM
C 2012.03 02/20/13 DAT(11,*) is UT time of MLT=0
C 2012.04 09/12/13 Replace HOUR with HOURUT in APFMSIS ---- P. Coisson
C 2012.05 01/22/14 TMAXN in GTD7 SEC->SECNI  HOUR->0.0 
C 2012.06 07/17/14 Change estromcor to estormcor -------- A.Shabanloui
C 2012.07 07/24/14 COMMON/iounit/: added 'mess'
C 2012.08 09/18/14 JF(18): FIELDG not UT_LT ............... A.Mazzella
C 2012.08 09/18/14 jf(12)&jf(34): create messages.txt ..... A.Mazzella
C 2012.08 09/18/14 change: icalls.gt.1 to icalls.ge.1  .... A.Mazzella
C 2012.09 09/24/14 added oarr(85)=L and oarr(86)=DIMO
C 2012.10 11/26/14 reading INDAP the first time
C 2012.11 12/22/14 COMMON/CSW/, ISW=0, SW(9)=-1 or =0 (no Ap depend.)
C 2012.12 07/12/15 adapting calls for TCON,APF,APF_ONLY,APFMSIS
C
C 2016.01 08/13/15 moving SWMI(25) from DATA into program
C 2016.01 08/13/15 add PI to COMMON/CONST; delete COMMON/CONST2
C 2016.01 08/23/15 Earth parameters now initialized in IRI_SUB
C 2016.02 09/14/15 JF(41)=t COV=F10.7_12, =f COV=f(IG12)
C 2016.03 09/16/15 observed F10.7 for GTD7 and CHEMION
C 2016.04 09/28/15 Calculate invdip parameter only once for 600 km
C 2016.05 09/30/15 hmF2: AMTB-2013 and SHU-2015 models; JF(39:40)
C 2016.05 09/30/15 revised ELTEIK and CALION calls
C 2016.06 10/14/15 added FELDCOF call for CLCMLT .......... M.Hausman
C 2016.06 10/14/15 COMMON/IGRF1/...,DIMO
C 2016.07 02/01/16 if(hef.le.hme)  no F1 and no valley .... M.Hausman
C 2016.08 06/01/16 User-specified B0 when jf(43)=false
C 2016.09 08/15/16 Corrected input of F10.7D,Y,81, and 365..M.Hausman
C 2016.09 08/15/16 ITOPN=3(Gulyaeva topside) not yet active M.Hausman
C 2016.10 09/08/16 CHEMION call now with n(H) input
C 2016.10 09/08/16 Replace SDMF2 with model_hmF2 (Shubin)
C 2016.11 09/22/16 COMMON: NmF2s,NmEs (STORM foF2, foE for profile)
C 2016.12 10/20/16 IG12_in->R12=f(IG12_in), R12_in->IG12=f(R12_in)
C 2016.12 10/20/16 R12=f_Gulyaeva(IG12_in), IG12=f_Gulyaeva(R12_in)
C 2016.12 10/20/16 F10.7_81_in -> F10.7_365 = F10.7_81_in
C 2016.13 01/26/17 B1 user input; 0.6<B1<6; B1_user only if B0_user
C 2016.14 02/23/17 XM3_CCIR for NeQuick; foF2s option for M3000F2in
C 2016.15 02/27/17 Using XHI1, XHI2, XHI3, XHI4 and related SAX, SUX
C 2016.15 02/27/17 No F1 layer if jf(19) and jf(20) are false
C 2016.16 10/27/17 F10.7Din and not F10.7_81in -> F10.7_81=F10.7Din 
C 2016.16 10/27/17 F10.7_81in and not F10.7Din -> F10.7D=F10.7_81in 
C 2016.17 10/30/17 OARR(87,88)=SAX300,SUX300
C 2016.18 03/22/18 f107in, f107ino, f107_81in, f107_81ino  M. Butala
C 2016.18 03/22/18 f107yo -> f107yobs, f10781o -> f10781obs
C 2016.18 03/22/18 invdip_old for Te elteik() ........... V. Truhlik
C 2016.18 03/22/18 ELTEIK and CALION use PF107OBS ....... V. Truhlik
C 2016.19 03/28/18 OAR(1:100)=-1 corrected
C 2016.20 04/23/18 Versioning now based on year of major releases
C 2016.21 04/25/18 Moved secni to Te calculation .......... C. Vasly
C 2016.21 04/25/18 Deleted arrays ddo and d2o; not used ... C. Vasly
C 2016.22 08/23/18 CNEW! B0-Gulyaeva option revised .... T. Gulyaeva
C 2016.23 08/27/18 Moved FIRI option under ELDE calc ..... P. Sultan
C 2016.23 08/27/18 3-h ap,kp available even if storm models are off
C 2016.23 08/27/18 daily ap avail. even if F10.7din or F10.7_81in
C 2016.24 08/29/18 user input for HNEA and HNEE if jf(45) jf(46)
C 2016.25 06/11/19 comments for OARR and output incl HNEA and HNEE
C
C 2020.01 07/03/19 changed argmax to 87.3 (consistent with IDL)
C 2020.01 07/03/19 itopn=1 now with PF10.7 correction
C 2020.02 07/19/19 itopn=1 cor option, itopn=3 cor2 option
C 2020.03 07/29/19 added 'endif' itopn=3 and declared a01(2,2)        
C 2020.04 08/05/19 itopn=3 requires itopn=1, BLO11 change        
C 2020.05 01/16/20 ion composition topside if h.ge.300km
C 2020.06 06/15/20 Den_NO=0.0,rn=0.0 before CHEMION .... M. Hausman
C 2020.06 06/15/20 IF(..hcor1) after IF(itopn....) ..... M. Hausman 
C 2020.07 08/24/20 changed if(sam_moye ...) ............ A. Rodland
C 2020.07 08/24/20 FELDCOF only when new date .......... A. Rodland             
C 2020.08 09/16/20 added HPOL for topside COR2 option
C 2020.08 09/16/20 not comp'd oarr variables: -1 or -100 for lati
C 2020.08 09/16/20 jf(47) for CGM computation on/off; lati>25
C 2020.09 01/03/21 Jf(19,20)=f no F1 -> C1=0, hmF1=0
C 2020.09 01/03/21 f1_c1 call with ABSMDP
C 2020.09 01/03/21 Intermediate region: hst 100<->hmF1(hmF2+hef)/2
C 2020.09 01/03/21 Intermediate region: T calculated in XE_3
C 2020.10 04/19/22 TEA(6 -> 4) used in TEBA call
C 2020.10 04/20/22 New Ti option Tru-2021 using IONTIF
C 2020.10 04/23/22 funct. ELTE changed to BOOKER1, COMMON/BLOTE del. 
C 2020.10 04/23/22 funct. TI changed to BOOKER1, COMMON/BLOCK8 del. 
C 2020.10 04/29/22 ROCSAT ion drift model of Fejer et al., 2008 
C 2020.10 04/29/22 --- Requires rocdrift.for file  
C 2020.11 08/05/22 F00 cal is now for new FIRI-2018 (IRIDREG.FOR)
C 2020.12 09/28/22 COR2: exp merging from hmF2 to hcor2   
C 2020.12 10/02/22 Changed fill value for vi to 9999 .... I. Girach 
C 2020.12 10/02/22 Fill val: -99 for lat. var, sza, sundec, IG12 
C 2020.13 11/28/22 Addded plasmasphere for IRI2001cor and COR2
C 2020.13 11/28/22 outf(15,*)= ratio plasma to gyro frequency
C 2020.14 02/07/23 Init. D_MSIS(1)=0 for GTD7(cira.for). M. Hausman 
C 2020.14 02/07/23 Ni: RBTT needs IAPO,TEH,TIH jf(2)!... M. Hausman 
C 2020.14 02/07/23 elteik,iontif calls: invdip_old at 600km
C 2020.14 02/07/23 calion call: invdip at height
C 2020.15 02/19/23 plasmapause fixed at L=5
C 2020.16 03/18/23 vfjmodelrocstart: array fjm specified 
C 2020.17 03/23/23 replacing subroutine iri_tec with IRITEC
C 2020.18 04/27/23 w/o plasmapause (fixpt:5000,10000,20000,30000km)
C 2020.19 10/22/23 iri_web: added h_tec_min
C 2020.19 10/22/23 added text if h<HNEA or h>HNEE
C 2020.20 10/25/23 jf(50)=t/f without/with plasmapause
C 2020.21 03/04/24 corrected: do 8429 kind=47,57,1.... J. Lilensten
C 2020.22 04/04/24 including sporadic-E occ. prob. model (ESPROB)
C 2020.22 11/05/24 TCOR2 above peak: 'h' -> 'height'
C 2020.23 11/24/24 Tru-2021: adding XTETI and Tn<Ti<Te requirement
C 2020.23 11/24/24 BIL-1985 until 3000km, TBT-2012 until 2000km
C 2020.23 11/24/24 Ne: HNEE=30000km
C 2020.24 11/29/24 B1.GE.0.6 changed to B1.GE.1.2 ........ B. Reid
C 2020.24 11/29/24 Avoid HF1.lt.hst ...................... B. Reid
C 2020.24 11/29/24 Changed HZ to 2*(HF1+hst)/3
C 2020.24 11/29/24 if h>hmF2 and Ne(h)>NmF2 then Ne(h)=NmF2
C 2020.25 12/04/24 Omit F1 if NmF1<1.2*NmE or NmF1>0.8*NmF2 
C 2020.25 12/04/24 Omit F1 if NmF1 < XE2(HEF) 
C 2020.26 12/22/24 ESPROB_new using F10.7_365 and Kp_sum
C 2020.26 12/22/24 jf(46)=t/f: ESPROB wout/with solar/mag depend.
C 2020.27 01/25/25 ESPROB call: MLONG=360 --> MLONG=0
C 2020.28 03/19/25 D-region: DELA -> DELL
C 2020.28 03/25/25 ******** released on irimodel.org *************
C 2020.29 07/19/25 ESPROB: change to xMLLOO,xKpsum ... K. Johnston
C 2020.30 07/23/25 xKpsum: change i=6,13 to i=5,12 .... V. Truhlik
C 2020.31 08/01/25 Sun declination: oarr(24)=SD300 ... K. Johnston
C 2020.31 08/01/25 ESPROB call: use invdip_old_110 
C 2020.32 08/08/25 spreadf_brazil call: daynr->daynr1. K. Johnston
C 2020.33 08/12/25 CALION: invdip only calculated at 1000km
C 2020.33 08/12/25 IONTIF: use invdip_old_600 
C 2020.33 08/12/25 If foF1in then f1pb=1.0 ........... K. Johnston
C 2020.34 08/30/25 If RZIN or IGIN no TCON 
C 2020.34 08/30/25 Improved Bilitza-1985 Te option
C 2020.35 09/15/25 Declare DAYNR1 as integer
C 2020.35 09/15/25 iri_web: dhour->xhour(=xvar(8)+iut*25.)
C 2020.35 09/15/25 iri_web: including MODA and int for 
C 2020.35 09/25/25 ******** released on irimodel.org *************
C
C 2026.01 10/18/25 xteti=ahh(7) and IF(X.GT.AHH(7)) ... V. Truhlik
C 2026.02 06/11/26 Replacing Te model TBT-2012 with TBPS-2026
C 2026.03 06/23/26 IONTIF: delete SIGTV - not used
C 2026.03 06/23/26 Adding IBP-2023 plasma bubble probability model
C 2026.03 06/23/26 spreadf_brazil: corrected month1 and daynr1
C 2026.03 06/23/26 change start of xteti computation to ahh(8)
C 2026.04 07/08/26 Ti,Te: GTD7(..SECNI->SEC,0->hour..) 
C 2026.04 07/08/26 Ti,Te: added numap, numti and rearranged
C
C*****************************************************************
C********* INTERNATIONAL REFERENCE IONOSPHERE (IRI). *************
C*****************************************************************
C**************** ALL-IN-ONE SUBROUTINE  *************************
C*****************************************************************
C
C
       SUBROUTINE IRI_SUB(JF,JMAG,ALATI,ALONG,IYYYY,MMDD,DHOUR,
     &    HEIBEG,HEIEND,HEISTP,OUTF,OARR)
C-----------------------------------------------------------------
C
C INPUT:  JF(1:50)      true/false switches for several options
C         JMAG          =0 geographic   = 1 geomagnetic coordinates
C         ALATI,ALONG   LATITUDE NORTH AND LONGITUDE EAST IN DEGREES
C         IYYYY         Year as YYYY, e.g. 1985
C         MMDD (-DDD)   DATE (OR DAY OF YEAR AS A NEGATIVE NUMBER)
C         DHOUR         LOCAL TIME (OR UNIVERSAL TIME + 25) IN DECIMAL 
C                          HOURS
C         HEIBEG,       HEIGHT RANGE IN KM; maximal 100 heights, i.e.
C          HEIEND,HEISTP        int((heiend-heibeg)/heistp)+1.le.100
C
C    JF switches to turn off/on (.true./.false.) several options
C
C    i       .true.                  .false.          standard version
C    -----------------------------------------------------------------
C    1    Ne computed            Ne not computed                     t
C    2    Te, Ti computed        Te, Ti not computed                 t
C    3    Ne & Ni computed       Ni not computed                     t
C    4    B0,B1 - Bil-2000       B0,B1 - other models see jf(31) false
C    5    foF2 - CCIR            foF2 - URSI                     false
C    6    Ni - DS-1995 & DY-1985 Ni - RBV-2010 & TBT-2015        false
C    7    Ne - Tops: f10.7<188   f10.7 unlimited                     t            
C    8    foF2 from model        foF2 or NmF2 - user input           t
C    9    hmF2 from model        hmF2 or M3000F2 - user input        t
C   10    Te - Standard          Te - Using Te/Ne correlation        t
C   11    Ne - Standard Profile  Ne - Lay-function formalism         t
C   12    Messages to unit 6     to messages.txt on unit 11          t
C   13    foF1 from model        foF1 or NmF1 - user input           t
C   14    hmF1 from model        hmF1 - user input (only Lay version)t
C   15    foE  from model        foE or NmE - user input             t
C   16    hmE  from model        hmE - user input                    t
C   17    Rz12 from file         Rz12 - user input (oarr(33))        t
C   18    IGRF dip, magbr, modip old FIELDG using POGO68/10 for 1973 t
C   19    F1 probability model           t
C   20     F1            standard F1 plus L condition        t
C   19,20 (t,t) f1-prob Scotto1997, (f,t) Ducharme1973 & not NIGHT t,t
C         (t,f) f1-prob-Scotto with L cond., (f,f) no F1
C   21    ion drift computed     ion drift not computed              t
C   22    ion densities in %     ion densities in m-3                t
C   23    Te_tops (Bil-1985)     Te_topside (TBPS-2026)          false
C   24    D-region: IRI-1990     FT-2001 and DRS-1995                t
C   25    F107D from APF107.DAT  F107D user input (oarr(41))         t
C   26    foF2 storm model       no storm updating                   t
C   27    IG12 from file         IG12 user input (oarr(39))          t
C   28    spread-F probability 	 not computed                        t
C   29,30 (t,t) IRI2001, (f,t) IRIcor, (f,f) NeQuick, (t,f) IRIcor2  t,f
C   31    B0,B1 ABT-2009	 B0 Gulyaeva-1987 h0.5               t   
C   32    F10.7_81 from file     F10.7_81 - user input (oarr(46))    t
C   33    Auroral boundary model on/off  true/false	         false
C   34    Messages on            Messages off                        t
C   35    foE storm model        no foE storm updating           false
C   36    hmF2 w/out foF2_storm  with foF2-storm                     t
C   37    topside w/out foF2-storm  with foF2-storm                  t
C   38    IBP-2023 bubbles model on/off   t/f                        t
C   39    hmF2 (M3000F2)         new models                      false
C   40    hmF2 AMTB-model        Shubin-COSMIC model             false
C   41    Use COV=F10.7_365      COV=f(IG12) (IRI before Oct 2015)   t
C   42    Te-TBPS with PF10.7 	 w/o PF10.7 dependance               t
C   43    B0 from model          B0 - user input (OARR(10))          t
C   44    B1 from model          B1 - user input (OARR(35))          t
C   45    Es occ. prob on        off                                 t
C   46    Es prob without solar and magnetic	Es prob with...      t
C   45,46 (t,t) ESPROB, (t,f) ESPRO with F10.7, (f,f) ESPRO off,   t,f
C         (f,t) ESPROB with F10.7 and kp dependence
C   47    CGM computation on 	 CGM computation off             false
C   48    Ti  Tru-2021           Bil-1981                            t
C   49    Plasmasphere: Ozhogin  Hartley model                       t
C   50    without plasmapause	 with plasmapause                    t
C   ------------------------------------------------------------------
C
C  Depending on the jf() settings additional INPUT parameters may 
c  be required:
C
C       Setting              INPUT parameter
C    -----------------------------------------------------------------
C    jf(8)  =.false.     OARR(1)=user input for foF2/MHz or NmF2/m-3
C    jf(9)  =.false.     OARR(2)=user input for hmF2/km or M(3000)F2
C    jf(10 )=.false.     OARR(15),OARR(16)=user input for Ne(300km),
C       Ne(400km)/m-3. Use OARR()=-1 if one of these values is not 
C       available. If jf(23)=.false. then Ne(300km), Ne(550km)/m-3.
C    jf(13) =.false.     OARR(3)=user input for foF1/MHz or NmF1/m-3 
C    jf(14) =.false.     OARR(4)=user input for hmF1/km
C    jf(15) =.false.     OARR(5)=user input for foE/MHz or NmE/m-3 
C    jf(16) =.false.     OARR(6)=user input for hmE/km
C    jf(17) =.false.     OARR(33)=user input for Rz12
C    jf(25) =.false.     OARR(41)=user input for daily F10.7 index
C    jf(27) =.false.     OARR(39)=user input for IG12
C    jf(32) =.false.     OARR(46)=user input for 81-day avg F10.7
C    jf(43) =.false.     OARR(10)=user input for B0
C    jf(44) =.false.     OARR(35)=user input for B1
C
C
C  OUTPUT:  OUTF(1:20,1:1000)
C               OUTF(1,*)  ELECTRON DENSITY/M-3
C               OUTF(2,*)  NEUTRAL TEMPERATURE/K
C               OUTF(3,*)  ION TEMPERATURE/K
C               OUTF(4,*)  ELECTRON TEMPERATURE/K
C               OUTF(5,*)  O+ ION DENSITY/% or /M-3 if jf(22)=f 
C               OUTF(6,*)  H+ ION DENSITY/% or /M-3 if jf(22)=f
C               OUTF(7,*)  HE+ ION DENSITY/% or /M-3 if jf(22)=f
C               OUTF(8,*)  O2+ ION DENSITY/% or /M-3 if jf(22)=f
C               OUTF(9,*)  NO+ ION DENSITY/% or /M-3 if jf(22)=f
C                 AND, IF JF(6)=.FALSE.:
C               OUTF(10,*)  CLUSTER IONS DEN/% or /M-3 if jf(22)=f
C               OUTF(11,*)  N+ ION DENSITY/% or /M-3 if jf(22)=f
C               OUTF(12,*)  
C               OUTF(13,*)  
C  if(jf(24)    OUTF(14,1:11) standard IRI-Ne for 60,65,..,110km 
C     =.false.)        12:22) Friedrich (FIRI) model at these heights 
C                      23:33) standard Danilov (SW=0, WA=0) 
C                      34:44) for minor Stratospheric Warming (SW=0.5) 
C                      45:55) for major Stratospheric Warming (SW=1) 
C                      56:66) weak Winter Anomaly (WA=0.5) conditions
C                      67:77) strong Winter Anomaly (WA=1) conditions
C               OUTF(15,*) plasma frequency divided by gyro frequency
C               OUTF(16-20,*)  free
c
C            OARR(1:100)   ADDITIONAL OUTPUT PARAMETERS         
C
C      #OARR(1) = NMF2/M-3           #OARR(2) = HMF2/KM
C      #OARR(3) = NMF1/M-3           #OARR(4) = HMF1/KM
C      #OARR(5) = NME/M-3            #OARR(6) = HME/KM
C       OARR(7) = NMD/M-3             OARR(8) = HMD/KM
C       OARR(9) = HHALF/KM           #OARR(10) = B0/KM
C       OARR(11) =VALLEY-BASE/M-3     OARR(12) = VALLEY-TOP/KM
C       OARR(13) = TE-PEAK/K          OARR(14) = TE-PEAK HEIGHT/KM
C      #OARR(15) = TE-MOD(300KM)     #OARR(16) = TE-MOD(400KM)/K
C       OARR(17) = TE-MOD(600KM)      OARR(18) = TE-MOD(1400KM)/K
C       OARR(19) = TE-MOD(3000KM)     OARR(20) = TE(120KM)=TN=TI/K
C       OARR(21) = TI-MOD(430KM)      OARR(22) = X/KM, WHERE TE=TI
C       OARR(23) = SOL ZENITH ANG/DEG OARR(24) = SUN DECLINATION/DEG
C       OARR(25) = DIP/deg            OARR(26) = DIP LATITUDE/deg
C       OARR(27) = MODIFIED DIP LAT.  OARR(28) = Geographic latitude
C       OARR(29) = sunrise/dec. hours OARR(30) = sunset/dec. hours
C       OARR(31) = ISEASON (1=spring) OARR(32) = Geographic longitude
C      #OARR(33) = Rz12               OARR(34) = Covington Index
C      #OARR(35) = B1                 OARR(36) = M(3000)F2
C      $OARR(37) = TEC/m-2           $OARR(38) = TEC_top/TEC*100.
C      #OARR(39) = gind (IG12)        OARR(40) = F1 probability 
C      #OARR(41) = F10.7 daily        OARR(42) = c1 (F1 shape)
C       OARR(43) = daynr              OARR(44) = equatorial vertical 
C       OARR(45) = foF2_storm/foF2_quiet            ion drift in m/s
C      #OARR(46) = F10.7_81           OARR(47) = foE_storm/foE_quiet 
C       OARR(48) = spread-F probability          
C       OARR(49) = Geomag. latitude   OARR(50) = Geomag. longitude  
C       OARR(51) = ap at current time OARR(52) = daily ap
C       OARR(53) = invdip/degree      OARR(54) = MLT-Te
C       OARR(55) = CGM-latitude       OARR(56) = CGM-longitude
C       OARR(57) = CGM-MLT            OARR(58) = CGM lat eq. aurl bodry
C       OARR(59) = CGM-lati(MLT=0)    OARR(60) = CGM-lati for MLT=1
C       OARR(61) = CGM-lati(MLT=2)    OARR(62) = CGM-lati for MLT=3
C       OARR(63) = CGM-lati(MLT=4)    OARR(64) = CGM-lati for MLT=5
C       OARR(65) = CGM-lati(MLT=6)    OARR(66) = CGM-lati for MLT=7
C       OARR(67) = CGM-lati(MLT=8)    OARR(68) = CGM-lati for MLT=9
C       OARR(69) = CGM-lati(MLT=10)   OARR(70) = CGM-lati for MLT=11
C       OARR(71) = CGM-lati(MLT=12)   OARR(72) = CGM-lati for MLT=13
C       OARR(73) = CGM-lati(MLT=14)   OARR(74) = CGM-lati for MLT=15
C       OARR(75) = CGM-lati(MLT=16)   OARR(76) = CGM-lati for MLT=17
C       OARR(77) = CGM-lati(MLT=18)   OARR(78) = CGM-lati for MLT=19
C       OARR(79) = CGM-lati(MLT=20)   OARR(80) = CGM-lati for MLT=21
C       OARR(81) = CGM-lati(MLT=22)   OARR(82) = CGM-lati for MLT=23
C       OARR(83) = Kp at current time OARR(84) = magnetic declination 
C       OARR(85) = L-value            OARR(86) = dipole moment 
C       OARR(87) = SAX300             OARR(88) = SUX300 
C       OARR(89) = HNEA               OARR(90) = HNEE 
C       OARR(91) = Es occ. prob in %  OARR(92) = bubble prob. in %  
C                # INPUT as well as OUTPUT parameter
C                $ special for iri_web (only place-holders)
C		for more details got to end of subroutine
c-----------------------------------------------------------------------        
C*****************************************************************
C*** THE ALTITUDE LIMITS ARE:  LOWER (DAY/NIGHT)  UPPER        ***
C***     ELECTRON DENSITY         60/80 KM       1500 KM       ***
C***     TEMPERATURES               60 KM        2500/3000 KM  ***
C***     ION DENSITIES             100 KM        1500 KM       ***
C*****************************************************************
C*****************************************************************
C*********            INTERNALLY                    **************
C*********       ALL ANGLES ARE IN DEGREE           **************
C*********       ALL DENSITIES ARE IN M-3           **************
C*********       ALL ALTITUDES ARE IN KM            **************
C*********     ALL TEMPERATURES ARE IN KELVIN       **************
C*********     ALL TIMES ARE IN DECIMAL HOURS       **************
C*****************************************************************
C*****************************************************************
C*****************************************************************
      INTEGER    DAYNR,DAYNR1,DDO,DO2,SEASON,SEADAY,MMO(12)
      REAL       LATI,LONGI,MO2,MO,MODIP,NMF2,MAGBR,INVDIP,IAPO,  
     &           NMF1,NME,NMD,MM,MLAT,MLONG,NMF2S,NMES,INVDPC,
     &           INVDIP_OLD,INVDIP_OLD_110,INVDIP_OLD_600,
     &           INVDPC_OLD
      DOUBLE PRECISION XLONGI,DPHOUR,DPF107D,IBPROB
      CHARACTER  FILNAM*12
c-web-for webversion
c      CHARACTER FILNAM*53

      DIMENSION  ARIG(3),RZAR(3),F(3),E(4),XDELS(4),DNDS(4),
     &  FF0(988),XM0(441),F2(13,76,2),FM3(9,49,2),ddens(5,11),
     &  elg(7),FF0N(988),XM0N(441),F2N(13,76,2),FM3N(9,49,2),
     &  INDAP(13),AMP(4),HXL(4),SCL(4),XSM(7),atii(7),MM(6),
     &  DTI(6),AHH(8),STTE(7),DTE(7),ATE(8),TEA(4),XNAR(2),
     &  OUTF(20,1000),DION(7),osfbr(25),D_MSIS(9),T_MSIS(2),
     &  IAPO(7),SWMI(25),ab_mlat(48),DAT(11,4),PLA(4),PLO(4),
     &  a01(2,2),teva(6),tiv(4),FM(59,25,4,11),atei(8),
     &  param(2),OARR(100),xlpl(4)

      DIMENSION palogne(6),dplas(4),pah(6),iap(13),fjm(59,25,4,11)

      LOGICAL  EXT,SCHALT,TECON(2),sam_mon,sam_yea,sam_ut,sam_date,
     &  F1REG,FOF2IN,HMF2IN,URSIF2,LAYVER,RBTT,DREG,rzino,FOF1IN,
     &  HMF1IN,FOEIN,HMEIN,RZIN,sam_doy,F1_OCPRO,F1_L_COND,NODEN,
     &  NOTEM,NOION,TENEOP,JF(50),URSIFO,igin,igino,mess,
     &  dnight,enight,fnight,fstorm_on,estorm_on,B0IN,B1IN,
     &  fof2ino,hmf2ino,f107in,f107ino,f107_81in,f107_81ino,
     &  sam_moye

      COMMON /CONST/UMR,PI  /const1/humr,dumr   /ARGEXP/ARGMAX
     & /IGRF1/ERA,AQUAD,BQUAD,DIMO	/BLOCK2/B0,B1,C1
     & /BLOCK1/HMF2,NMF2S,HMF1,F1REG	/BLOCK3/HZ,T,HST   
     & /BLOCK4/HME,NMES,HEF	/BLOCK5/ENIGHT,E
     & /BLOCK6/HMD,NMD,HDX	/BLOCK7/D1,XKK,FP30,FP3U,FP1,FP2
     & /BLO10/BETA,ETA,DELTA,ZETA /BLO11/B2TOP,itopn,tcor1,tcor2   
      COMMON /findRLAT/FLON,RYEAR  /cotec/hnea,hplpa
     & /iounit/konsol,mess     /CSW/SW(25),ISW,SWC(25)
     & /QTOP/Y05,H05TOP,QF,XNETOP,XM3000,HHALF,TAU 

      COMMON /testtei/numti,xsm,atii,numap,ahh,ate
      COMMON /testne/hpla(4),xnpla(4)

      EXTERNAL  XE1,XE2,XE3_1,XE4_1,XE5,XE6,FMODIP

      DATA icalls/0/, dplas/100,150,10,10/,jfirsta,jfirste/0,0/
c      DATA DTE/5.,5.,5.,5.,10.,20.,20./, DTI/5.,5.,10.,20.,20.,20./
      DATA DTE/10.,10.,10.,10.,20.,30.,30./
      DATA DTI/10.,10.,20.,30.,30.,30./
      DATA MMO/31,28,31,30,31,30,31,31,30,31,30,31/
      DATA hpla/5000.,10000.,20000.,30000./

        save
                
        mess=jf(34)
        
c set switches for NRLMSIS00  
        ISW=0
        do 6492 KI=1,25
6492    SWMI(KI)=1.

        nummax=1000
        DO 7397 KI=1,20
        do 7397 kk=1,nummax
7397    OUTF(KI,kk)=-1.
C
C oarr(1:6,10,15,16,33,35,39,41,46) are used for inputs.
C The fill value is -1 for most oarr output parameters. It is 
C -99 for all latitudinal variables, solar zenith angle, sun
C declination and IG12. It is 9999 for ion drift.
C 
        oarr(7)=-1.
        oarr(8)=-1.
        oarr(9)=-1.
        do 8398 kind=11,14,1
8398    	oarr(kind)=-1.
        do 8378 kind=17,22,1
8378    	oarr(kind)=-1.
        do 8478 kind=23,28,1
8478    	oarr(kind)=-99.
        do 8479 kind=29,32,1
8479    	oarr(kind)=-1.
        oarr(34)=-1.
        oarr(36)=-1.
        oarr(37)=-1.
        oarr(38)=-1.
        oarr(40)=-1.
        oarr(42)=-1.
        oarr(43)=-1.
        oarr(44)=9999.
        oarr(45)=-1.
        do 8429 kind=47,57,1
8429       oarr(kind)=-1.
        oarr(49)=-99.
        oarr(53)=-99.
        oarr(55)=-99.
        do 8428 kind=58,82,1
8428       oarr(kind)=-99.
        oarr(84)=-99.
        do 8427 kind=85,100,1
8427       oarr(kind)=-1.

C
C PROGRAM CONSTANTS AND INITIALIZATION
C
        if(icalls.lt.1) then
        	ARGMAX=87.3
        	pi=ATAN(1.0)*4.
        	UMR=pi/180.
        	humr=pi/12.
        	dumr=pi/182.5
        	ALOG2=ALOG(2.)
        	ALG10=ALOG(10.)
        	ALG100=ALOG(100.)
        	montho=-1
		nmono=-1
		iyearo=-1
		idaynro=-1
		rzino=.true.
		igino=.true.
		f107ino=.true.
		f107_81ino=.true.
		fof2ino=.true.
		hmf2ino=.true.
		ut0=-1
		ursifo=.true.
C Initialize parameters for COMMON/IGRF1/
C   ERA		EARTH RADIUS (WGS-84: 6371.137 KM) 
C   EREQU   MAJOR HALF AXIS FOR EARTH ELLIPSOID (6378.160 KM)
C   ERPOL   MINOR HALF AXIS FOR EARTH ELLIPSOID (6356.775 KM)
C   AQUAD   SQUARE OF MAJOR HALF AXIS FOR EARTH ELLIPSOID
C   BQUAD   SQUARE OF MINOR HALF AXIS FOR EARTH ELLIPSOID
C   EEXC	Eccentricity of Earth's orbit
C   DIMO    Earth's dipole moment in Gauss 
C ERA, EREQU and ERPOL as recommended by the INTERNATIONAL 
C ASTRONOMICAL UNION .
        	ERA=6371.2
        	EREQU=6378.16
        	ERPOL=6356.775
        	AQUAD=EREQU*EREQU
        	BQUAD=ERPOL*ERPOL
        	EEXC=0.01675
        	dimo=0.311653
C Initialize D_MSIS(1) to avoid accidental activation of 
C user input for Tn-exospheric in calls to GTD7 (CIRA.FOR)
                D_MSIS(1) = 0.0
        	endif
 
        numhei=int(abs(heiend-heibeg)/abs(heistp))+1
        if(numhei.gt.nummax) numhei=nummax
C
C Code inserted to aleviate block data problem for PC version.
C Thus avoiding DATA statement with parameters from COMMON block.
C
        XDELS(1)=5.
        XDELS(2)=5.
        XDELS(3)=5.
        XDELS(4)=10.
        DNDS(1)=.016
        DNDS(2)=.01
        DNDS(3)=.016
        DNDS(4)=.016
        XNAR(1)=0.0
        XNAR(2)=0.0
C
C FIRST SPECIFY YOUR COMPUTERS CHANNEL NUMBERS ....................
C AGNR=OUTPUT (OUTPUT IS DISPLAYED OR STORED IN FILE OUTPUT.IRI)...
C IUCCIR=UNIT NUMBER FOR CCIR COEFFICIENTS ........................
c
        IUCCIR=10
c-web- special for web version
c-web- messages should be turned off with mess=jf(34)=.false. 

        KONSOL=6
        if(.not.jf(12).and.mess) then
           konsol=11
           open(11,file='messages.txt')
           endif
c
c selection of density, temperature and ion composition options ......
c

      NODEN=(.not.jf(1))
      NOTEM=(.not.jf(2))
      NOION=(.not.jf(3))
      if(.not.NOION) NODEN=.false.
      RBTT=(.not.jf(6))
      LAYVER=(.not.jf(11))
      F1_OCPRO=(jf(19))
      F1_L_COND=(.not.jf(20))
      DREG=jf(24)
      fstorm_on=jf(26).and.jf(8)
      estorm_on=jf(35).and.jf(15)
c
c rz12, IG12, F10.7D, PF10.7 input option ............................
c
      RZIN=(.not.jf(17))
      IF(RZIN) THEN
          ARZIN=OARR(33)
      else
          oarr(33)=-1.
      ENDIF
      
      IGIN=(.not.jf(27))
      IF(IGIN) THEN
          AIGIN=OARR(39)
      else
          oarr(39)=-99.
      ENDIF

      F107IN=(.not.jf(25))
      IF(F107IN) THEN
          f107din=OARR(41)
      else
          oarr(41)=-1.
      ENDIF

      F107_81IN=(.not.jf(32))
      IF(F107_81IN) THEN
          f10781in=OARR(46)
      else
          oarr(46)=-1.
      ENDIF
c
c Topside density ....................................................
c
        if(jf(29)) then
             if (jf(30)) then    
                 itopn=0		! IRI2001 topside option
             else
                 itopn=3		! IRI-cor2 topside option
             endif
        else 
             if (jf(30)) then
                 itopn=1		! IRI-cor topside option
             else
                 itopn=2		! NeQuick topside option
             endif
        endif
c
c F2 peak density ....................................................
c
      FOF2IN=(.not.jf(8))
       IF(FOF2IN) THEN
          OARR1=OARR(1)
       	  AFOF2=OARR1
       	  ANMF2=OARR1
          IF(OARR1.LT.100.) ANMF2=1.24E10*AFOF2*AFOF2
          IF(OARR1.GE.100.) AFOF2=SQRT(ANMF2/1.24E10)
        else
          oarr(1)=-1.
        ENDIF
      URSIF2=(.not.jf(5))
c
c F2 peak altitude ..................................................
c
      HMF2IN=(.not.jf(9))
       IF(HMF2IN) then
         AHMF2=OARR(2)
       else
         oarr(2)=-1.
       endif
c
c F1 peak density ...................................................
c
       FOF1IN=(.not.jf(13))
       IF(FOF1IN) THEN
          OARR3=OARR(3)
       	  AFOF1=OARR3
       	  ANMF1=OARR3
          IF(OARR3.LT.100.) ANMF1=1.24E10*AFOF1*AFOF1
          IF(OARR3.GE.100.) AFOF1=SQRT(ANMF1/1.24E10)
       else
          oarr(3)=-1.
       ENDIF
c
c F1 peak altitude ..................................................
c
       HMF1IN=(.not.jf(14))
       IF(HMF1IN) then
          AHMF1=OARR(4)
          if(.not.layver.and.mess) write(konsol,1939)
1939      format(' *Ne* User input of hmF1 is only possible',
     &      ' for the LAY-version')
       else
          oarr(4)=-1.
       endif
c
c E peak density ....................................................
c
       FOEIN=(.not.jf(15))
       IF(FOEIN) THEN
          OARR5=OARR(5)
       	  AFOE=OARR5
       	  ANME=OARR5
          IF(OARR5.LT.100.) ANME=1.24E10*AFOE*AFOE
          IF(OARR5.GE.100.) AFOE=SQRT(ANME/1.24E10)
       else
          oarr(5)=-1.
       ENDIF
c
c E peak altitude ..................................................
c
       HMEIN=(.not.jf(16))
       IF(HMEIN) then
          AHME=OARR(6)
       else
          oarr(6)=-1.
       endif
c
c B0 bottomside thickness ...........................................
c
       B0IN=(.not.jf(43))
       IF(B0IN) then
          B0_US=OARR(10)
       else
          oarr(10)=-1.
       endif
c
c B1 bottomside profile shape ......................................
c
       B1IN=(.not.jf(44))
       if(jf(43)) B1IN=.false.
       IF(B1IN) then
          B1_US=OARR(35)
       else
          oarr(35)=-1.
       endif
c
c sporadic-E occurrence probability .................................
c
       if(jf(45)) then
         if (jf(46)) then    
           iesn=0		! ESPROB Es probability model
         else
           iesn=1	        ! ESPROB with F10.7 dependence
         endif
       else 
         if (jf(46)) then
           iesn=2		! ESPROB with F10.7 and ap dep.
         else
           iesn=3		! ESPROB off
         endif
       endif
C
C TE MODEL OPTIONS ..............................................
C
       TENEOP=(.not.jf(10))
       IF(TENEOP) THEN
           DO 8154 JXNAR=1,2
              XNAR(JXNAR)=OARR(JXNAR+14)
              TECON(JXNAR)=.FALSE. 
8154          IF(XNAR(JXNAR).GT.0.) TECON(JXNAR)=.TRUE. 
       else
           oarr(15)=-1.
           oarr(16)=-1.
       ENDIF
c
c lists the selected options before starting the table
c
       if(icalls.ge.1.or.(.not.mess)) goto 8201
          write(konsol,2911) 
          if(NODEN) goto 2889
          if(LAYVER) write(konsol,9012) 
          if(jf(49)) then
             write(konsol,9229) 
          else
             write(konsol,9228) 
          endif
          if(jf(50)) then
             write(konsol,9227) 
          else
             write(konsol,9226) 
          endif
          if(.not.JF(7)) write(konsol,9014) 
          if (itopn.eq.0) write(konsol,9207)
          if (itopn.eq.1) write(konsol,9204)
          if (itopn.eq.2) write(konsol,9205)
          if (itopn.eq.3) write(konsol,9206)
          if(FOF2IN) then
             write(konsol,9015)
          else 
             if(URSIF2) then
                write(konsol,9016) 
             else
                write(konsol,9017) 
             endif
          endif
          if(HMF2IN) then
             write(konsol,9018)
          else 
             if(jf(39)) then
                write(konsol,9020)
             else 
                if(jf(40)) then
                   write(konsol,9717) 
                else
                   write(konsol,9716) 
                endif
             endif
          endif
          if(fof1in) write(konsol,9019) 
          if(HMF1IN.and.LAYVER) write(konsol,9021) 
          if(jf(4)) then
             write(konsol,9214)
          else 
             if(jf(31)) then
                write(konsol,9216)
             else
                write(konsol,9215)
             endif
          endif
          if(foein) write(konsol,9022) 
          if(HMEIN) write(konsol,9023)
          if(jf(35)) then
             write(konsol,9128)
          else
             write(konsol,9129)
          endif
          if(B0IN) write(konsol,9923) 
          if(B1IN) write(konsol,9927) 
          if(F1_OCPRO) then
             write(konsol,9024) 
             if(F1_L_COND) write(konsol,9025)
             endif
          if(.not.f1_ocpro.and.f1_l_cond) write(konsol,2917)
          if(iesn.eq.0) write(konsol,9036)
          if(iesn.eq.1) write(konsol,9037)
          if(iesn.eq.2) write(konsol,9038)
          if(iesn.eq.3) write(konsol,9040)
          if(jf(38)) write(konsol,9939)
          if(DREG) then 
             write(konsol,9026) 
          else
             write(konsol,9027) 
          endif
          if(jf(26)) then
             if(fof2in) then 
                write(konsol,9028) 
                jf(26)=.false.
             else
                write(konsol,9029) 
             endif
          endif

2889    continue
        if(jf(33)) then
            write(konsol,4031)
        else
            write(konsol,4032)
        endif

        if((.not.NOION).and.(RBTT)) write(konsol,9031) 
        if((.not.NOION).and.(.not.RBTT)) write(konsol,9039) 

        if(NOTEM) goto 8201
          if(TENEOP) write(konsol,9032) 
          if(jf(23)) then
            write(konsol,9033)		! Bil-1985 Te option
          else
            if (jf(42)) then
              write(konsol,9034)	! TBPS-2026 with F10.7
            else
              write(konsol,9535)	! TBPS-2026 without
            endif
          endif
        if(jf(48)) then
            write(konsol,9333)
        else
            write(konsol,9335)
        endif

2911    format('*** IRI parameters are being calculated ***')
9012    format('Ne: The LAY-Version is prelimenary.',
     &          ' Erroneous profile features can occur.')
9014    format('Ne: No upper limit for F10.7 in',
     &          ' topside formula.')
9204    format('Ne: IRI-cor for Topside')
9205    format('Ne: NeQuick for Topside')
9206    format('Ne: IRI-cor2 for Topside')
9207    format('Ne: IRI-2001 for Topside')
9214    format('Ne: Bil-2000 for B0 and B1')
9215    format('Ne: Gul-1987 for B0')
9216    format('Ne: ABT-2009 for B0 and B1')
9015    format('Ne: foF2 or NmF2 provided by user')
9016    format('Ne: URSI-1989 model for foF2')
9017    format('Ne: CCIR-1967 model for foF2')
9716    format('Ne: Shubin-2015 model for hmF2')
9717    format('Ne: AMBT-2012 model for hmF2')
9020    format('Ne: hmF2(M3000F2) BSE-1979 ')
9018    format('Ne: hmF2 or M3000F2 provided by user')
9019    format('Ne: foF1 or NmF1 provided by user')
9021    format('Ne: hmF1 provided by user')
9022    format('Ne: foE or NmE provided by user')
9023    format('Ne: hmE provided by user')
9923    format('Ne: B0 provided by user')
9927    format('Ne: B1 provided by user')
9024    format('Ne: foF1 probability function used')
9025    format('Ne: foF1 L condition cases included')
2917    format('Ne: profile without F1 layer')
9036    format('Ne: sporadic-E occ. pro. model on')
9037    format('Ne: Es ESPROB with F10.7 dependence')
9038    format('Ne: Es ESPROB with F10.7 and ap dependence')
9040    format('Ne: sporadic-E occ. prob. model off')
9939	format('Ne: bubble occ. prob. model IBP on ')
9026    format('Ne: IRI1990 is used for D region')
9027    format('Ne: IRI-1990, FT-2001, DRS-1995 available for D region')
9028    format('Ne: Storm model is turned off if foF2 or',
     &          ' NmF2 user input')
9029    format('Ne: foF2 storm model included')
9128    format('Ne: foE storm model included')
9229    format('Ne: Ohzogin-2012 for plasmasphere')
9228    format('Ne: Hartley-2023 for plasmasphere')
9227    format('Ne: profile without plasmapause')
9226    format('Ne: profile with plasmapause')
9129    format('Ne: foE storm model turned off')
9039    format('Ni: DS-1995 below 300km and DY-1985 above')
9031    format('Ni: RBV-2010 below 300km and TBT-2015 above')
9032    format('Te: Temperature-density correlation is used')
9033    format('Te: Bil-1985 model')
9034    format('Te: TBPS-2026 with PF10.7 dependence')
9535    format('Te: TBPS-2026 without PF10.7 dependence')
9333    format('Ti: TBKS-2021 model with PF10.7 dependence')
9335    format('Ti: Bil-1981 model')
4031    format('AB: Auroral boundary model on')
4032    format('AB: Auroral boundary model off')

8201    continue

C
C CALCULATION OF DAY OF YEAR OR MONTH/DAY AND DECIMAL YEAR 
c NRDAYM is the number of days in the current month 
c IDAYY is the number of days in the current year
c
c  leap year rule: years evenly divisible by 4 are leap years, except
c  years also evenly divisible by 100 are not leap years, except years 
c  also evenly divisible by 400 are leap years. The year 2000 is a 100 
c  and 400 year exception and therefore it is a normal leap year. 
c  The next 100 year exception will be in the year 2100!
c

        iyear=iyyyy
        idayy=365
        if(iyear/4*4.eq.iyear) idayy=366    ! leap year

        if(MMDD.lt.0) then
           DAYNR=-MMDD
           call MODA(1,iyear,MONTH,IDAY,DAYNR,nrdaym)
        else
           MONTH=MMDD/100
           IDAY=MMDD-MONTH*100
           call MODA(0,iyear,MONTH,IDAY,DAYNR,nrdaym)
        endif

        ryear = iyear + (daynr-1.0)/idayy
        iyd = iyear *1000 + daynr
        amx = pi*(daynr-3.)/182.6
        radj = 1.-eexc*(cos(amx)+eexc*(cos(2*amx)-1.)/2.)

C
C calculate center height for CGM computation
C

        height_center=(HEIBEG+HEIEND)/2.
        

C
C CALCULATION OF GEODETIC/GEOMAGNETIC COORDINATES (LATI, LONGI AND 
C MLAT, MLONG), MAGNETIC INCLINATION (DIP), DIP LATITUDE (MAGBR) 
C AND MODIFIED DIP (MODIP), ALL IN DEGREES
C

        if(along.lt.0.) along = along + 360. ! -180/180 to 0-360
        
        IF(JMAG.GT.0) THEN
           MLAT=ALATI
           MLONG=ALONG
        ELSE
           LATI=ALATI
           LONGI=ALONG
        ENDIF
        CALL GEODIP(IYEAR,LATI,LONGI,MLAT,MLONG,JMAG)

        if((iyear.ne.iyearo).or.(daynr.ne.idaynro)) CALL FELDCOF(RYEAR)

        if(jf(18)) then
          call igrf_dip(lati,longi,ryear,300.0,dec,dip,magbr,modip)
        else
          CALL FIELDG(LATI,LONGI,300.0,XMA,YMA,ZMA,BET,DIP,DEC,MODIP)
          MAGBR=ATAN(0.5*TAN(DIP*UMR))/UMR
        endif
c
c calculate L-value, dip lati, and B_abs needed for invdip computation
c calculating invdip at 1000 km for Ni-TBT-2015 and invdip_old for
c Te-TBT-2012 at 1000km and for Ti-Tru-2021 at 600km.
c
        invdip=-99.0
        invdip_old=-99.0
        invdip_old_110=-99.0
        invdip_old_600=-99.0
        call igrf_sub(lati,longi,ryear,110.0,fl,icode,dipl,babs)
           if(fl.gt.10.) fl=10.
      	   invdip_old_110=INVDPC_OLD(FL,DIMO,BABS,DIPL)
        if(jf(3).and.(.not.jf(6))) then
           call igrf_sub(lati,longi,ryear,1000.0,fl,icode,dipl,babs)
           if(fl.gt.10.) fl=10.
           invdip=INVDPC(FL,DIMO,BABS,DIPL)
	   endif
        if(jf(2).and.(.not.jf(23))) then
           call igrf_sub(lati,longi,ryear,1000.0,fl,icode,dipl,babs)
           if(fl.gt.10.) fl=10.
      	   invdip_old=INVDPC_OLD(FL,DIMO,BABS,DIPL)
	   endif
        if(jf(2).and.jf(48)) then
           call igrf_sub(lati,longi,ryear,600.0,fl,icode,dipl,babs)
           if(fl.gt.10.) fl=10.
      	   invdip_old_600=INVDPC_OLD(FL,DIMO,BABS,DIPL)
	   endif

        ABSLAT=ABS(LATI)
        ABSMLT=ABS(MLAT)
        ABSMDP=ABS(MODIP)
        ABSMBR=ABS(MAGBR)
c 
C CALCULATION OF UT/LT and XMLT  ...............
c
        IF(DHOUR.le.24.0) then
           HOUR=DHOUR					! dhour =< 24 is LT
           hourut=hour-longi/15.
           if(hourut.lt.0) hourut=hourut+24.
        ELSE
           hourut=DHOUR-25.				 ! dhour>24 is UT+25
           hour=hourut+longi/15.	 		 ! hour becomes LT
           if(hour.gt.24.) hour=hour-24.
        ENDIF
        CALL CLCMLT(IYEAR,DAYNR,HOURUT,LATI,LONGI,XMLT)
c
c calculation of CGM coordinates if abs(latitude) > 25 degrees
c
        cgm_lat=-99.0
        cgm_lon=-99.0
        cgm_mlt=-1.0
c        if(jf(47).and.(abslat.gt.25.0)) then
        if(jf(47)) then
	   DAT(1,1)=lati
	   DAT(2,1)=longi
           call GEOCGM01(1,IYEAR,height_center,DAT,PLA,PLO)
c           cgm_lat=DAT(3,1)
c           cgm_lon=DAT(4,1)
c           cgm_mlt00_ut=DAT(11,1)
           cgm_lat=DAT(3,3)
           cgm_lon=DAT(4,3)
           cgm_mlt00_ut=DAT(11,3)
           cgm_mlt=hourut-cgm_mlt00_ut
           if(cgm_mlt.lt.0.) cgm_mlt=24.+hourut-cgm_mlt00_ut
           endif
c
c SEASON assumes equal length seasons (92 days) with spring 
c (SEASON=1) starting at day-of-year=45; for lati < 0 adjustment 
c for southern hemisphere is made. Some models require the
c seasonal month (ISEAMON) or the seasonal day-of year (SEADAY)
c ZMONTH is decimal month (Jan 1 = 1.0 and Dec 31 = 12.97)
c SDAY is the day number reduced to a 360 day year (TOPH05)
c NRDAYM is the number of days in the current month 
c IDAYY is the number of days in the current year
c 
      
        SEASON=INT((DAYNR+45.0)/92.0)
        IF(SEASON.LT.1) SEASON=4
        NSEASN=SEASON		    ! Northern hemisphere season
        zmonth = month + (iday-1)*1./nrdaym
        seaday=daynr
        iseamon=month
        IF(LATI.lt.0.0) then
          SEASON=SEASON-2			
          IF(SEASON.LT.1) SEASON=SEASON+4
          iseamon=month+6
          if(iseamon.gt.12) iseamon=iseamon-12
          seaday=daynr+idayy/2.
          if(seaday.gt.idayy) seaday=seaday-idayy
          endif
C
C 12-month running mean sunspot number (rssn) and Ionospheric Global 
C index (gind), daily F10.7 cm solar radio flux (f107d) and monthly 
C F10.7 (cov) index   
C
        sam_mon=(month.eq.montho)
        sam_yea=(iyear.eq.iyearo)
        sam_doy=(daynr.eq.idaynro)
        sam_date=(sam_yea.and.sam_doy)
        sam_moye=(sam_yea.and.sam_mon)
        sam_ut=(hourut.eq.ut0)
        
        if(sam_date.and..not.rzin.and..not.rzino.and..not.igin
     &     .and..not.igino.and..not.f107in.and..not.f107ino
     &     .and..not.f107_81in.and..not.f107_81ino) goto 2910     
 
       if(RZIN.or.IGIN) then
           midm=15
           if(month.eq.2) midm=14
           call MODA(0,iyear,month,midm,idd1,nrdyr)
           imm2=month+1
           if(imm2.gt.12) then
              imm2=1
              idd2=380            ! =365+15 mid-January
              if(iyear/4*4.eq.iyear) idd2=381
           else
              call MODA(0,iyear,imm2,midm,IDD2,nrdym)
           endif
           ttt=(daynr-idd1)*1./(idd2-idd1)
           nmonth=imm2
        else
           call tcon(iyear,month,iday,daynr,rzar,arig,ttt,nmonth)
           if(nmonth.lt.0) goto 3330
        endif
	
        if(RZIN) then
           rrr = arzin
           rzar(1) = rrr
           rzar(2) = rrr
           rzar(3) = rrr
           if(IGIN) then
              zi = aigin
           else 
              zi=(-0.0031*ARZIN+1.5332)*ARZIN-11.5634
           endif
           arig(1) = zi
           arig(2) = zi
           arig(3) = zi
           endif
         	
        if(IGIN) then
           zi = aigin
           arig(1) = zi
           arig(2) = zi
           arig(3) = zi
           if(RZIN) then
              rrr = arzin
           else
              if(zi.gt.178.0066) zi=178.0066
              rrr = 247.29 - 17.96*sqrt(178.0066-xigin)
           endif             
           rzar(1) = rrr
           rzar(2) = rrr
           rzar(3) = rrr
           endif

        rssn=rzar(3)
        gind=arig(3)
        COV=63.75+RSSN*(0.728+RSSN*0.00089)
C
C Getting F10.7 index: daily (f107d), previous day (f107y; 
C required by MSIS), 81-day average (f10781), 365-day average
C (f107365), and PF10.7=(F10.7_daily + F10.7_81_day)/2. 
C F10.7 should be adjusted (to top of atmosphere) value not 
C observed (at the ground) value.         
        f107d=cov
        f107y=cov
        f10781=cov
        f107365=cov
        call APF_ONLY(iyear,month,iday,F107_daily,F107PD,
     &      F107_81,F107_365,IAP_daily,isdate)
        if(.not.f107in.or..not.f107_81in) then
          if(F107_daily.gt.-11.1) then
            f107d=f107_daily
            f107y=f107PD
            f10781=f107_81
            f107365=f107_365
            endif
	  endif
        if(f107in) then
          f107d=f107din             ! user input: F10.7 daily 
          f107y=f107din             ! same input: F10.7 previous day
          if(.not.f107_81in) then
            f10781=f107din          ! same input: F10.7 81-day
            f107365=f107din         ! same input: F10.7 yearly average
            endif	        		
          endif
        if(f107_81in) then
          f10781=f10781in 	    ! user input: F10.7 81-day average
          f107365=f10781in	    ! same input: F10.7 yearly average
          if(.not.f107in) then
            f107d=f10781in          ! same input: F10.7 daily	
            f107y=f10781in          ! same input: F10.7 previous day
            endif	        		
          endif		
          pf107=(f107d+f10781)/2.

c Correcting F10.7 adjusted flux from APF107.DAT to flux observed at
c Earth that is expected by NRLMSIS00 (GTD7) and CHEMION AND TBT-2015
C (ION COMPOSITION) AND TBT-2012 (ELECTRON TEMPERATURE)          

        f_adj=radj*radj        
        f107yobs=f107y/f_adj
        f10781obs=f10781/f_adj
        pf107obs=pf107/f_adj
		       
        if(jf(41)) cov=f107365				
        COVSAT=cov
        if(covsat.gt.188.) covsat=188

C
C CALCULATION OF SOLAR ZENITH ANGLE IN DEGREE FOR THE USER-SPECIFIED 
C TIME (HOUR) AND FOR NOON, SUN DECLINATION ANGLE (SUNDEC), AND TIME 
C OF LOCAL SUNRISE/SUNSET (SAX, SUX in DECIMAL HOURS) AT 80 KM (D-
C REGION), 110 KM (E-REGION), 200 KM (F1-REGION), AND 300 KM (F-
C REGION AND TOPSIDE). 
C

2910    continue
        CALL SOCO(daynr,HOUR,LATI,LONGI,80.,SUNDEC,XHI1,SAX80,SUX80)
        CALL SOCO(daynr,HOUR,LATI,LONGI,110.,SUD1,XHI2,SAX110,SUX110)
        CALL SOCO(daynr,HOUR,LATI,LONGI,200.,SUD1,XHI3,SAX200,SUX200)
        CALL SOCO(daynr,HOUR,LATI,LONGI,300.,SD300,XHI4,SAX300,SUX300)
        CALL SOCO(daynr,12.0,LATI,LONGI,110.,SUNDE1,XHINON,SAX1,SUX1)
        CALL SOCO(daynr,12.0,LATI,LONGI,200.,SUNDE2,XHINON2,SAX2,SUX2)
        DNIGHT=.FALSE.
        if(abs(sax80).gt.25.0) then
           if(sax80.lt.0.0) DNIGHT=.TRUE.
           goto 9334
           endif
        if(SAX80.le.SUX80) goto 1386
        if((hour.gt.sux80).and.(hour.lt.sax80)) dnight=.true.
        goto 9334
1386    IF((HOUR.GT.SUX80).OR.(HOUR.LT.SAX80)) DNIGHT=.TRUE.

9334    ENIGHT=.FALSE.
        if(abs(sax110).gt.25.0) then
           if(sax110.lt.0.0) ENIGHT=.TRUE.
           goto 8334
           endif
        if(SAX110.le.SUX110) goto 9386
        if((hour.gt.sux110).and.(hour.lt.sax110)) enight=.true.
        goto 8334
9386    IF((HOUR.GT.SUX110).OR.(HOUR.LT.SAX110)) ENIGHT=.TRUE.

8334    FNIGHT=.FALSE.
        if(abs(sax200).gt.25.0) then
           if(sax200.lt.0.0) FNIGHT=.TRUE.
           goto 1334
           endif
        if(SAX200.le.SUX200) goto 7386
        if((hour.gt.sux200).and.(hour.lt.sax200)) fnight=.true.
        goto 1334
7386    IF((HOUR.GT.SUX200).OR.(HOUR.LT.SAX200)) FNIGHT=.TRUE.

C
C CALCULATION OF ELECTRON DENSITY PARAMETERS................
C lower height boundary (HNEA), upper boundary (HNEE)
C
      
1334  continue

      HNEE = 30000.
      HNEA = 65.
      IF(DNIGHT) HNEA = 80.
      IF(NODEN) GOTO 4933

      DELA = 4.32
      IF(ABSMDP.GE.18.) DELA=1.0+EXP(-(ABSMDP-30.0)/10.0)
      DELL=1+EXP(-(ABSLAT-20.)/10.)
c
c E peak critical frequency (foE), density (NmE), and height (hmE)
c
        IF(FOEIN) THEN
          FOE=AFOE
          NME=ANME
        ELSE
          FOE=FOEEDI(COV,XHI2,XHINON,ABSLAT)
          NME=1.24E10*FOE*FOE
        ENDIF
        IF(HMEIN) THEN
          HME=AHME
        ELSE
          HME=110.0
        ENDIF
c
c sporadic E occurrence probability ESP 
c
        call apf(isdate,hourut,indap)

        if(iesn.lt.3) then
           xKpsum=0
           do i=5,12
              xxkp=ckp(indap(i))
              xKpsum=xKpsum+xxkp
              enddo
           isode=0
           if(iesn.gt.0) isode=1
           imade=0
           if(iesn.gt.1) imade=1
           XMLLOO = MLONG
           if(XMLLOO.eq.360.0) XMLLOO = 0.0
           CALL ESPROB(isode,imade,INVDIP_OLD_110,XMLLOO,XMLT,
     &       DAYNR,F107_365,xKpsum,ESP)          
           endif
c
c ionospheric bubble probability (IBP) model (valid range: 
C 350-510km, LT: 18 to 6, -45 to +45 deg latitude)
c         
         ibprob=-1.0
         if(.not.jf(38).or.abs(LATI).gt.45.0) goto 2356
         if(hour.lt.18.0.and.hour.gt.6.0) goto 2356
           CALL IBPLOD
           xlongi=dble(longi)
           dphour=dble(hour)
           dpf107d=dble(f107d)
           if(longi.gt.180) xlongi=xlongi-360.0
           CALL IBPEMP(MONTH,DAYNR,DPHOUR,XLONGI,DPF107D,IBPROB)
           xibp=sngl(ibprob)
C
C F2 peak critical frequency foF2, density NmF2, and height hmF2
C
C READ CCIR AND URSI COEFFICIENT SET FOR CHOSEN MONTH 
C
2356  IF((FOF2IN).AND.(HMF2IN).and.(itopn.ne.2)) GOTO 501
      IF((FOF2INO).OR.(HMF2INO)) GOTO 7797
      IF(URSIF2.NEQV.URSIFO) GOTO 7797
      IF(sam_moye.AND.(nmonth.NE.nmono)) GOTO 4293
      if(rzin.or.igin.or.rzino.or.igino) then
        IF(sam_moye.AND.(nmonth.EQ.nmono)) GOTO 4291
      else
        IF(sam_moye.AND.(nmonth.EQ.nmono)) GOTO 4292
      endif

7797    URSIFO=URSIF2
        WRITE(FILNAM,104) MONTH+10
104         FORMAT('ccir',I2,'.asc')
c-web-for webversion
c104     FORMAT('/var/www/omniweb/cgi/vitmo/IRI/ccir',I2,'.asc')
        OPEN(IUCCIR,FILE=FILNAM,STATUS='OLD',ERR=8448,
     &          FORM='FORMATTED')
        READ(IUCCIR,4689) F2,FM3
4689    FORMAT(1X,4E15.8)
        CLOSE(IUCCIR)
C
C then URSI if chosen ....................................
C
        if(URSIF2) then
          WRITE(FILNAM,1144) MONTH+10
1144          FORMAT('ursi',I2,'.asc')
c-web-for webversion
c1144    FORMAT('/var/www/omniweb/cgi/vitmo/IRI/ursi',I2,'.asc')
          OPEN(IUCCIR,FILE=FILNAM,STATUS='OLD',ERR=8448,
     &         FORM='FORMATTED')
          READ(IUCCIR,4689) F2
          CLOSE(IUCCIR)
          endif

C
C READ CCIR AND URSI COEFFICIENT SET FOR NMONTH, i.e. previous 
c month if day is less than 15 and following month otherwise 
C

4293    continue

c
c first CCIR ..............................................
c

        WRITE(FILNAM,104) NMONTH+10
        OPEN(IUCCIR,FILE=FILNAM,STATUS='OLD',ERR=8448,
     &          FORM='FORMATTED')
        READ(IUCCIR,4689) F2N,FM3N
        CLOSE(IUCCIR)

C
C then URSI if chosen .....................................
C
        if(URSIF2) then
          WRITE(FILNAM,1144) NMONTH+10
          OPEN(IUCCIR,FILE=FILNAM,STATUS='OLD',ERR=8448,
     &         FORM='FORMATTED')
          READ(IUCCIR,4689) F2N
          CLOSE(IUCCIR)
          endif

        GOTO 4291
        
8448    WRITE(konsol,8449) FILNAM
8449    FORMAT(1X////,
     &    ' The file ',A30,'is not in your directory.')
        GOTO 3330
C
C LINEAR INTERPOLATION IN SOLAR ACTIVITY. IG12 used for foF2
C

4291    continue
        RR2=ARIG(1)/100.
        RR2N=ARIG(2)/100.
        RR1=1.-RR2
        RR1N=1.-RR2N
        DO I=1,76
           DO J=1,13
              K=J+13*(I-1)
              FF0N(K)=F2N(J,I,1)*RR1N+F2N(J,I,2)*RR2N
              FF0(K)=F2(J,I,1)*RR1+F2(J,I,2)*RR2
              enddo
           enddo

        RR2=RZAR(1)/100.
        RR2N=RZAR(2)/100.
        RR1=1.-RR2
        RR1N=1.-RR2N
        DO I=1,49
           DO J=1,9
              K=J+9*(I-1)
              XM0N(K)=FM3N(J,I,1)*RR1N+FM3N(J,I,2)*RR2N
              XM0(K)=FM3(J,I,1)*RR1+FM3(J,I,2)*RR2
              enddo
           enddo

4292    zfof2  =  FOUT(MODIP,LATI,LONGI,HOURUT,FF0)
        fof2n  =  FOUT(MODIP,LATI,LONGI,HOURUT,FF0N)
        zm3000 = XMOUT(MODIP,LATI,LONGI,HOURUT,XM0)
        xm300n = XMOUT(MODIP,LATI,LONGI,HOURUT,XM0N)
        midm=15
        if(month.eq.2) midm=14
        if (iday.lt.midm) then
           yfof2 = fof2n + ttt * (zfof2-fof2n)
           xm3000= xm300n+ ttt * (zm3000-xm300n)
        else
           yfof2 = zfof2 + ttt * (fof2n-zfof2)
           xm3000= zm3000+ ttt * (xm300n-zm3000)
        endif
        XM3_CCIR=XM3000
			
501     IF(FOF2IN) THEN
          FOF2=AFOF2
          NMF2=ANMF2
        ELSE
          FOF2=YFOF2
          NMF2=1.24E10*FOF2*FOF2
        ENDIF
c
c stormtime updating for foF2 (foF2s, NmF2s) and foE (foEs,
c NmEs) and auroral boundary computation.
c
        foF2s=foF2
        foEs=foE
        NmF2s=NmF2
        NmEs=NmE
        stormcorr=-1.
        estormcor=-1.
        index_3h_ap=indap(13)
        if(index_3h_ap.gt.-1) then
 	   xkp=ckp(index_3h_ap)
        else     
           xkp=3.0
	endif   
           
c
c stormtime updating for foF2 (foF2s, NmF2s) 
c
        if(fstorm_on.and.(indap(1).gt.-1)) then    
           kut=int(hourut)
c using the internal corrected geomagnetic latitude
c           icoord=2
c           call STORM(indap,lati,longi,icoord,cgm_lat,kut,
c     &          daynr,stormcorr)
c using the corrected geomagnetic latitude computed from
c latitude and longitude
           icoord=1
           call STORM(indap,lati,longi,icoord,cglat,kut,
     &          daynr,stormcorr)
           fof2s=fof2*stormcorr
           NMF2S=1.24E10*FOF2S*FOF2S
           endif
c
c stormtime updating for foE (foEs, NmEs)
c
        if(estorm_on.and.(index_3h_ap.gt.-1)) then    
           estormcor=STORME_AP(DAYNR,MLAT,index_3h_ap*1.0)
           if(estormcor.gt.-2.0) foes=foe*estormcor
           NMES=1.24E10*FOES*FOES
           endif
c
c calculation of equatorward auroral boundary: corrected magnetic
c latitude (CGM) of the boundary (ab_mlat(1:48)) for MLT=0.0,0.5,
c 1.0 ... 23.5 h and kp=xkp.
c 
        cgmlat=-99.0
        do 2929 i=1,48  
2929       ab_mlat(i)=-99.0
        if(jf(33)) then
           zmlt=xmlt
c           zmlt=cgm_mlt
           if(zmlt.lt.0.0.or.zmlt.gt.24.0) zmlt=-1.0 
           call auroral_boundary(xkp,zmlt,cgmlat,ab_mlat)
           endif
c
c determine latitude RLAT for AMTB-2013 hmF2 model and ABT-2009 B0 model
c            
       IF(((.not.JF(4)).and.JF(31)).or.((.not.JF(39)).and.JF(40))) THEN
	   FLON=LONGI+15.*hourut
           if(FLON.gt.360.) FLON=FLON-360.
           X11=-90
           X22=90
           FX11=fmodip(x11)
           FX22=fmodip(x22)
        CALL REGFA1(X11,X22,FX11,FX22,0.001,MODIP,FMODIP,SCHALT,XRLAT)
           IF(SCHALT) THEN
              XRLAT=LATI
              if(mess) WRITE(KONSOL,656) LATI
656     FORMAT(1X,'*NE* ABT-B0 computed with RLAT=LATI=',F6.2)
              endif
           RLAT=XRLAT
           ENDIF
c          	
c Computation of hmF2: user input or 3 model options		
c
    	IF(HMF2IN) THEN
           IF(AHMF2.LT.50.0) THEN
              XM3000=AHMF2
              ratf=fof2/foe
c if jf(36)=false then foF2_storm in hmF2 formula  
              if(.not.jf(36)) ratf=fof2s/foe
              HMF2=HMF2ED(MAGBR,RSSN,RATF,XM3000)
           ELSE 
              HMF2=AHMF2
c xm3000 calculation from hmF2 user input is no longer needed 
c because NeQuick requires CCIR-M(3000)F2 input always.
	      xm3000=-1.0
           endif
           goto 9917
           ENDIF

        IF(JF(39)) THEN
           ratf=fof2/foe
c if jf(36)=false then foF2_storm in hmF2 formula  
           if(.not.jf(36)) ratf=fof2s/foe
           HMF2=HMF2ED(MAGBR,RSSN,RATF,XM3000)
           goto 9917
           endif

	IF(JF(40)) THEN
c AMTB digisonde model		
 	   CALL SHAMDHMF2(RLAT,FLON,ZMONTH,RSSN,HMF2)
           goto 9917
           endif

c SHUBIN-COSMIC model
	CALL model_hmF2(iday,month,hourut,modip,
     &	   longi,F10781,HMF2)

9917    nmono=nmonth
        MONTHO=MONTH
        iyearo=iyear
        idaynro=daynr
        rzino=rzin
        igino=igin
        f107ino=f107in
        f107_81ino=f107_81in
        ut0=hourut
        foF2ino=foF2in
        hmF2ino=hmF2in
c
c topside profile parameters .............................
c
      if (itopn.eq.2) then     !NeQuick topside parameters 
         fo2=foF2s
         if(jf(37)) fo2=foF2
c NeQuick option requires CCIR-M3000F2, ignores user input.
         dNdHmx=-3.467+1.714*log(fo2)+2.02*log(XM3_CCIR)
         dNdHmx=exp(dNdHmx)*0.01
         B2bot=0.04774*fo2*fo2/dNdHmx
         b2k = 3.22-0.0538*fo2-0.00664*hmF2+0.113*hmF2/B2bot
     &   	+0.00257*rssn
         ee=exp(2.0*(b2k-1.0))
         b2k=(b2k*ee+1.0)/(ee+1.0)
         B2TOP=b2k*B2bot
      else
         COS2=COS(MLAT*UMR)
         COS2=COS2*COS2
         FLU=(COVSAT-40.0)/30.0
c option to use unlimited F10.7M for the topside
         IF(.not.JF(7)) FLU=(COV-40.0)/30.0
         FO1 = FOF2S
         IF(JF(37)) FO1 = FOF2
         EX=EXP(-MLAT/15.)
         EX1=EX+1
         EPIN=4.*EX/(EX1*EX1)
         ETA1=-0.02*EPIN
         ETA = 0.058798 + ETA1 - 
     &     FLU*(0.014065-0.0069724*COS2)+FO1* 
     &     (0.0024287+0.0042810*COS2-0.0001528*FO1)
         ZETA = 0.078922 - 0.0046702 * COS2 -  
     &     FLU*(0.019132-0.0076545*COS2)+FO1*
     &     (0.0032513+0.0060290*COS2-0.00020872*FO1)
         BETA = -128.03 + 20.253 * COS2 -
     &     FLU*(8.0755+0.65896*COS2)+FO1* 
     &     (0.44041+0.71458*COS2-0.042966*FO1)
         Z=EXP(94.5/BETA)
         Z1=Z+1
         Z2=Z/(BETA*Z1*Z1)
         DELTA=(ETA/Z1-ZETA/2.0)/(ETA*Z2+ZETA/400.0)
      endif

c
c Calculation of parameters for 2001cor and COR2 topside
c and extension to plasmasphere  
c
      TCOR1=0.0
      TCOR2=0.0
      if(itopn.eq.1.or.itopn.eq.3) then
        zmp111 = EPLA(modip,10.0,0.0)
        zmp222 = EPLA(modip,19.0,0.0) 
        r2n = -0.84 - 1.60 * zmp111  
        r2d = -0.84 - 0.64 * zmp111 
        x1n = 230. -  700. * zmp222  
        x1d = 550. - 1900. * zmp222
        r2 = HPOL(HOUR,r2d,r2n,SAX300,SUX300,1.,1.)
        x1 = HPOL(HOUR,x1d,x1n,SAX300,SUX300,1.,1.)
        hcor1 = hmF2 + x1
c        x12 = 1500. - x1
c        tc3 = r2 / x12
        hcor2 = (hcor1 + hmF2)/2.0
        shc = (hcor2-hmF2) / alog(2.0)
c
c Assuming a dipole magnetic field for |maglat|<60 degrees. 
c The L-value is related to height h by L = (1 + h/ERA)/(cos(maglat)^2.  
c 
        ppmlat=mlat
        ppb=60.0
        if(abs(ppmlat).gt.ppb) ppmlat=sign(ppb,ppmlat)
        cosmag=cos(ppmlat*umr)
        cosmag2=cosmag*cosmag

c
c Version without Plasmapause (hpla(4):5000,10000,20000,30000km)
c
   
        hplpa=hpla(2)
        if(jf(50)) then
          do i=1,4
            xlpl(i)=(1.0+hpla(i)/ERA)/cosmag2
            enddo
        else

c
c Version with Plasmapause: 
c Plasmapause inner/outer boundary fixed at L=5/5.1 
c

          xlpl(2)=5.0
          xlpl(3)=5.1		
          hpla(2)=ERA*(xlpl(2)*cosmag2-1.0)
          hpla(3)=ERA*(xlpl(3)*cosmag2-1.0)
          hpla(1)=(hpla(2)+hmF2+1500.0)/2.0
          hpla(4)=2*hpla(2)-hpla(1)
          xlpl(1)=(1.0+hpla(1)/ERA)/cosmag2
          xlpl(4)=(1.0+hpla(4)/ERA)/cosmag2
        endif

c		
c plasmasphere models (Hartley-2023, Ohzogin-2012)
c		

        do i=1,4
          if(jf(49)) then
            xnpla(i)=ohzden(xlpl(i),ppmlat)
          else
            xnpla(i)=harden(xlpl(i),xmlt,ppmlat)
          endif
          enddo

c
c Version with Plasmapause: 
c Density decrease by 10 at plasmapause or use
c Carpenter & Anderson (1992) for extended plasma trough  
c

        if(.not.jf(50)) then
          xnpla(3)=xnpla(2)/10.0
c          xnpla(3)=caadenet(xlpl(3),xmlt)
          xnpla(4)=caadenet(xlpl(4),xmlt)
          endif

c		
c Booker parameters for plasmaspheric extension
c		
        PAH(1)=hcor1
        PAH(2)=hmF2 + 1500.0
        PALOGNE(1)=0.0
        PALOGNE(2)=r2*alg10
        do i=1,4
          PAH(i+2)=hpla(i)		
          PALOGNE(i+2)=log(xnpla(i)/nmf2s)
          enddo

        tcor1=0.0
        do i=1,4        
          CALL SOCO(daynr,HOUR,LATI,LONGI,hpla(i),SUC,zxz,sap,sup)
          TCOR2=TCOR2CAL(hpla(i),hour,modip,pf107,sap,sup)
          znep=log(XE_1(hpla(i))/nmf2s)
          palogne(i+2)=palogne(i+2)-znep
          enddo
      endif

c
c Bottomside thickness parameter B0 and shape parameters B1
c                              
        if(jf(4)) then
          B0=B0_98(HOUR,SAX200,SUX200,NSEASN,RSSN,LONGI,MODIP)
	  B1=HPOL(HOUR,1.9,2.6,SAX200,SUX200,1.,1.)
        else if (jf(31)) then
 	  CALL SHAMDB0D (RLAT,FLON,ZMONTH,RSSN,B0)
          CALL SHAB1D (LATI,FLON,ZMONTH,RSSN,B1)
        else
          CALL ROGUL(SEADAY,XHI3,SEAX,GRAT)
cnew!        IF (FNIGHT) GRAT = 0.91D0 - HMF2/4000.D0
	  B1=HPOL(HOUR,1.9,2.6,SAX200,SUX200,1.,1.)
          BCOEF = B1*(B1*(0.0046D0*B1-0.0548D0)+0.2546D0) + 0.3606D0
          B0CNEW = HMF2*(1.D0-GRAT)
          B0 = B0CNEW/BCOEF
        endif
        if(B0IN) B0=B0_US
        if(B1IN) B1=B1_US
        if(B1.gt.6.0) B1=6.0
c        if(B1.lt.0.6) B1=0.6        
        if(B1.lt.1.2) B1=1.2        

c
c F1 layer height hmF1, critical frequency foF1, peak density NmF1
c profile parameter c1, and F1REG (=.false. no F1 region) 
c "No F1 layer" can be enforced by setting JF(19) and JF(20) to .false.
c
        if(.not.jf(19).and..not.jf(20)) then 
          F1REG=.false.
          FOF1=-1.0
          NMF1=-1.0
          hmF1=0.0
          c1=0.0
          if(mess) WRITE(KONSOL,1123)
1123      FORMAT(1X,'*NE* User selected No_F1_region option')
          goto 2918
          endif
        	
c
c F1 layer thickness parameter c1
c
        c1 = f1_c1(absmdp,hour,sax2,sux2)
c
c User provided foF1 or NmF1
c
        IF(FOF1IN) THEN
          FOF1=AFOF1
          NMF1=ANMF1
          F1REG=.TRUE.
          goto 2918
        ELSE
          FOF1=FOF1ED(ABSMBR,RSSN,XHI3)
          NMF1=1.24E10*FOF1*FOF1
        ENDIF
c
c F1 occurrence probability: jf(19), jf(20)
c jf(19),jf(20)=t,t Scotto et al. 1997
c jf(19),jf(20)=f,t Ducharme et al. 1973
c jf(19),jf(20)=t,f Scotto et al. 1997 with L-condition 
c jf(19),jf(20)=f,f no F1 layer 
c
        if(fof1in) then
          f1pb=1.0
        else
          if(f1_ocpro) then
            call f1_prob(xhi3,mlat,rssn,f1pbw,f1pbl)
            f1pb = f1pbw
            if(f1_l_cond) f1pb = f1pbl
          else			
            f1pb = 0.0 
            if((.not.fnight).and.(fof1.gt.0.0)) f1pb=1. 
          endif
        endif
        f1reg=.false.
        if(f1pb.ge.0.5) f1reg=.true.
2918	continue            
c
c E-valley: DEPTH=(NmE-N_deepest)/NmE*100, WIDTH=HEF-HmE, 
c distance of deepest value point above E-peak(HDEEP), 
c derivative at valley top divided by NmE (DLNDH), 
c and height of valley top (HEF)
c
      XDEL=XDELS(SEASON)/DELA
      DNDHBR=DNDS(SEASON)/DELA
      HDEEP=HPOL(HOUR,10.5/DELA,28.,SAX110,SUX110,1.,1.)
      WIDTH=HPOL(HOUR,17.8/DELA,45.+22./DELA,SAX110,SUX110,1.,1.)
      DEPTH=HPOL(HOUR,XDEL,81.,SAX110,SUX110,1.,1.)
      DLNDH=HPOL(HOUR,DNDHBR,.06,SAX110,SUX110,1.,1.)
      IF(DEPTH.LT.1.0) GOTO 600
        IF(ENIGHT) DEPTH=-DEPTH
        CALL TAL(HDEEP,DEPTH,WIDTH,DLNDH,EXT,E)
        IF(.NOT.EXT) GOTO 667
        if(mess) WRITE(KONSOL,650)
650     FORMAT(1X,'*NE* E-REGION VALLEY CAN NOT BE MODELLED')
600   WIDTH=.0
667   HEF=HME+WIDTH
      VNER = (1. - ABS(DEPTH) / 100.) * NMES

c
c Parameters below E  .............................
c

2727  continue
      hmex=hme-9.
      NMD=XMDED(XHI1,RSSN,4.0E8)
      HMD=HPOL(HOUR,81.0,88.0,SAX80,SUX80,1.,1.)
      F(1)=HPOL(HOUR,0.02+0.03/DELL,0.05,SAX80,SUX80,1.,1.)
      F(2)=HPOL(HOUR,4.6,4.5,SAX80,SUX80,1.,1.)
      F(3)=HPOL(HOUR,-11.5,-4.0,SAX80,SUX80,1.,1.)
      FP1=F(1)
      FP2=-FP1*FP1/2.0
      FP30=(-F(2)*FP2-FP1+1.0/F(2))/(F(2)*F(2))
      FP3U=(-F(3)*FP2-FP1-1.0/F(3))/(F(3)*F(3))
      HDX=HMD+F(2)
c
c indermediate region between D and E region; parameters xkk
c and d1 are found such that the function reaches hdx/xdx/dxdh
c
       X=HDX-HMD
       XDX=NMD*EXP(X*(FP1+X*(FP2+X*FP30)))
       DXDX=XDX*(FP1+X*(2.0*FP2+X*3.0*FP30))
       X=HME-HDX
       XKK=-DXDX*X/(XDX*ALOG(XDX/NMES))
c
c if exponent xkk is larger than xkkmax, then xkk will be set to 
c xkkmax and d1 will be determined such that the point hdx/xdx is 
c reached; derivative is no longer continuous.
c
        xkkmax=5.
        if(xkk.gt.xkkmax) then
                xkk=xkkmax
                d1=-alog(xdx/nmes)/(x**xkk)
        else
                D1=DXDX/(XDX*XKK*X**(XKK-1.0))
        endif
c
c compute Danilov et al. (1995) D-region model values
c
      if(.not.dreg) then
          vKp=1.
          f5sw=0.
          f6wa=0.
          call DRegion(xhi1,month,f107d,vKp,f5SW,f6WA,elg)
          do ii=1,11
            ddens(1,ii)=-1.
            if(ii.lt.8) ddens(1,ii)=10**(elg(ii)+6)
            enddo
          f5sw=0.5
          f6wa=0.
          call DRegion(xhi1,month,f107d,vKp,f5SW,f6WA,elg)
          do ii=1,11
            ddens(2,ii)=-1.
            if(ii.lt.8) ddens(2,ii)=10**(elg(ii)+6)
            enddo
          f5sw=1.
          f6wa=0.
          call DRegion(xhi1,month,f107d,vKp,f5SW,f6WA,elg)
          do ii=1,11
            ddens(3,ii)=-1.
            if(ii.lt.8) ddens(3,ii)=10**(elg(ii)+6)
            enddo
          f5sw=0.
          f6wa=0.5
          call DRegion(xhi1,month,f107d,vKp,f5SW,f6WA,elg)
          do ii=1,11
            ddens(4,ii)=-1.
            if(ii.lt.8) ddens(4,ii)=10**(elg(ii)+6)
            enddo          
          f5sw=0.
          f6wa=1.
          call DRegion(xhi1,month,f107d,vKp,f5SW,f6WA,elg)
          do ii=1,11
            ddens(5,ii)=-1.
            if(ii.lt.8) ddens(5,ii)=10**(elg(ii)+6)
            enddo
          endif
C
C SEARCH FOR HMF1 ..................................................
C
       if(LAYVER) goto 6153
       hmf1=0
       IF(.not.F1REG) GOTO 380
c
c Omit F1 feature if NmF1 is smaller than NmE*1.2 or
c greater than NmF2*0.8
c
       bnme=1.2*nmes
       bnmf2=0.8*nmf2
       if(nmf1.lt.bnme.or.nmf1.gt.bnmf2) goto 9427
c
c Omit F1 layer if NmF1 is below Ne given by the bottomside
c function at the E-valley top height. 
c
       XXE1=XE2(HEF)
       if(nmf1.lt.xxe1) goto 9427

       CALL REGFA1(HEF,HMF2,XXE1,NMF2S,0.001,NMF1,XE2,SCHALT,HMF1)
       IF(.not.SCHALT) GOTO 380

c
c omit F1 feature ....................................................
c

9427    if(mess) WRITE(KONSOL,11)
11      FORMAT(1X,'*NE* F1 omitted because NmF1 < 1.2*NmE or NmF1 ', 
     &  '> 0.8*NmF2'/1X,' or NmF1 < Ne(valley-top) or REGFA1 problem.')
        HMF1=0.
        F1REG=.FALSE.
        C1=0.0

C
C SEARCH FOR HST [NE3(HST)=NMEs] ......................................
C

380     continue

        IF(F1REG) then
            hf1=hmf1
            xf1=nmf1
        else
            hf1=2.*(hmf2+hef)/3.
            xf1=xe2(hf1)
        endif

        hf2=100.0
3762    xf2=xe3_1(hf2)
        if(xf2.gt.nmes) then
            hf2=hf2-10.0
            if(hf2.lt.hnea) goto 3885
            goto 3762
            endif
    
3763    if(xf1.lt.nmes) then
            hf1=hf1+10.0
            if(hf1.gt.hmF2-10) goto 3885
            xf1=xe2(hf1)
            goto 3763
            endif

        CALL REGFA1(hf1,HF2,XF1,XF2,0.001,NMES,XE3_1,SCHALT,HST)
        if(schalt) goto 3885
        HZ=(HST+HF1)/2.0
        GOTO 4933

c
c assume linear interpolation between HZ and HEF ..................
c

3885    if(mess) WRITE(KONSOL,100)
100     FORMAT(1X,'*NE* HST IS NOT EVALUATED BY THE FUNCTION XE3')
        HZ=(HEF+HF1)/2.
        xnehz=xe3_1(hz)
        if(mess) WRITE(KONSOL,901) HZ,HEF
901     FORMAT(6X,'CORR.: LIN. APP. BETWEEN HZ=',F5.1,
     &          ' AND HEF=',F5.1)
        T=(XNEHZ-NMES)/(HZ-HEF)
        HST=-333.
        GOTO 4933
C
C LAY-functions for middle ionosphere
C
6153    IF(HMF1IN) THEN
          HMF1M=AHMF1
        ELSE
          HMF1M=165.+0.6428*XHI3
        ENDIF
        HHALF = GRAT * HMF2
        HV1R = HME + WIDTH
        HV2R = HME + HDEEP
        HHMF2 = HMF2
        CALL INILAY(FNIGHT,F1REG,NMF2S,NMF1,NMES,VNER,HHMF2,HMF1M,HME,
     &                  HV1R,HV2R,HHALF,HXL,SCL,AMP,IIQU)
        IF((IIQU.EQ.1).and.mess) WRITE(KONSOL,7733)
7733   FORMAT('*NE* LAY amplitudes found with 2nd choice of HXL(1).')
        IF((IIQU.EQ.2).and.mess) WRITE(KONSOL,7722)
7722   FORMAT('*NE* LAY amplitudes could not be found.')

C
C---------- CALCULATION OF NEUTRAL TEMPERATURE PARAMETER-------
C

4933  HTA=60.0
      HEQUI=120.0
      IF(NOTEM.and.(NOION.or..not.RBTT)) GOTO 240
      SEC=hourut*3600.
      CALL APFMSIS(ISDATE,HOURUT,IAPO)
      if(iapo(2).lt.0.0) then
           SWMI(9)=0.
           IAPO(1)=0.
      else
           SWMI(9)=-1.0
      endif           
      CALL TSELEC(SWMI)
      CALL GTD7(IYD,SEC,HEQUI,LATI,LONGI,HOUR,F10781OBS,
     &        F107YOBS,IAPO,0,D_MSIS,T_MSIS)
      TN120=T_MSIS(2)

C
C--------- CALCULATION OF ELECTRON TEMPERATURE PARAMETER--------
C

881   CONTINUE

c Te(120km) = Tn(120km)

         AHH(1)=120.
         ATE(1)=TN120

c Te-MAXIMUM based on JICAMARCA and ARECIBO data 

        HMAXD=60.*EXP(-(MLAT/22.41)**2)+210.
        HMAXN=150.
        AHH(2)=HPOL(HOUR,HMAXD,HMAXN,SAX200,SUX200,1.,1.)
        TMAXD=800.*EXP(-(MLAT/33.)**2)+1500.
        secni=(24.-longi/15)*3600.
        CALL GTD7(IYD,SECNI,HMAXN,LATI,LONGI,0.0,F10781OBS,
     &        F107YOBS,IAPO,0,D_MSIS,T_MSIS)
        TMAXN=T_MSIS(2)
        ATE(2)=HPOL(HOUR,TMAXD,TMAXN,SAX200,SUX200,1.,1.)

          if(jf(23)) then

c BIL-1985 model: Te(300km), Te(400km) from AE-C, Te(1400km) 
c from ISIS-2, Te(3000km) from ISIS-1 (Brace and Theis, 1981)
              numap=7
              DIPLAT=MAGBR
              CALL TEBA(DIPLAT,HOUR,NSEASN,TEA)
              AHH(3)=300.
              AHH(4)=400.
              AHH(5)=600.
              AHH(6)=1400.
              AHH(7)=3000.
              ATE(3)=TEA(1)
              ATE(4)=TEA(2)
              ATE(6)=TEA(3)
              ATE(7)=TEA(4)

c Te(600km) from AEROS, Spenner and Plugge (1979)
              ETT=EXP(-MLAT/11.35)
              TET=2900.-5600.*ETT/((ETT+1)**2.)
c              TEN=839.+1161./(1.+EXP(-(ABSMLT-45.)/5.))
              TEN=839.0+1161.0/(1.0+EXP(-(ABSMLT-50.)/10.))
              ATE(5)=HPOL(HOUR,TET,TEN,SAX300,SUX300,1.5,1.5)

          else

c TBPS-2026 model:Te at 310, 415, 580, 850, 1400, 2000 km
c with or without solar activity effects included 
              numap=8
              AHH(3)=310.
              AHH(4)=415.
              AHH(5)=580.
              AHH(6)=850.
              AHH(7)=1400.
              AHH(8)=2000.
c isa=0/1 for solar activity correction off/on
              isa=0
              if(jf(42)) isa=1
              call elteo(isa,invdip_old,xmlt,daynr,
     &              pf107obs,teva)
              do ijk=3,numap
                 ate(ijk)=teva(ijk-2)
                 enddo
          endif
          numap1=numap-1
          numap2=numap-2


c Option to use Te = f(Ne) relation at ahh(3), ahh(4)
          IF(TENEOP) THEN
              DO I=1,2 
                 IF(TECON(I)) ATE(I+2)=
     &              TEDE(AHH(I+2),XNAR(I),-COV)
                  enddo
              endif

c The Te values at the fixpoints are slightly different from the actual 
c values, because of the analytical representation (Bilitza, ASR, 1985). 
c They are corrected to virtual values resulting in a better agreement. 
c And Te > Tn is enforced.
      CALL GTD7(IYD,SEC,AHH(2),LATI,LONGI,HOUR,F10781OBS,
     &        F107YOBS,IAPO,0,D_MSIS,T_MSIS)
      TNAHH2=T_MSIS(2)
      IF(ATE(2).LT.TNAHH2) ATE(2)=TNAHH2
      STTE1=(ATE(2)-ATE(1))/(AHH(2)-AHH(1))
      DO I=2,numap1
         CALL GTD7(IYD,SEC,AHH(I+1),LATI,LONGI,HOUR,F10781OBS,
     &        F107YOBS,IAPO,0,D_MSIS,T_MSIS)
         TNAHHI=T_MSIS(2)
         IF(ATE(I+1).LT.TNAHHI) ATE(I+1)=TNAHHI
         STTE2=(ATE(I+1)-ATE(I))/(AHH(I+1)-AHH(I))
c         atei(i)=ate(i)
         ATE(I)=ATE(I)-(STTE2-STTE1)*DTE(I-1)*ALOG2
         STTE1=STTE2
         enddo

c Te gradients STTE are computed for each segment
      DO I=1,numap1
         STTE(I)=(ATE(I+1)-ATE(I))/(AHH(I+1)-AHH(I))
         enddo
      ATE1=ATE(1)
887   CONTINUE
      HTE = AHH(numap)

C
C------------ CALCULATION OF ION TEMPERATURE PARAMETERS--------
C

c Starting height is 200km. Below 200km Ti=Tn
      	HS=200.
      	XSM(1)=HS
      	CALL GTD7(IYD,SEC,HS,LATI,LONGI,HOUR,F10781OBS,F107YOBS,
     &        IAPO,0,D_MSIS,T_MSIS)
      	TNHS=T_MSIS(2)
        atii(1)=tnhs

      if(jf(48)) then
c
c Tru-2021 model: TIV(1:4) ion temperature at 350,430,600,850km
c JSOL=0/1 solar activity dependence excluded/included
        numti=5
        numti2=numti-2
	JSOL=1
	xsm(2)=350.
	xsm(3)=430.
	xsm(4)=600.
	xsm(5)=850.
      	CALL IONTIF(JSOL,INVDIP_OLD_600,XMLT,DAYNR,PF107OBS,TIV)
c Tn < Ti < Te enforced at xsm(2:5)
        do i10=2,5
          HTIX=XSM(i10)
          TEXSM=BOOKER1(HTIX,numap2,ATE1,AHH,STTE,DTE)	
      	  CALL GTD7(IYD,SEC,htix,LATI,LONGI,hour,
     &        F10781OBS,F107YOBS,IAPO,0,D_MSIS,T_MSIS)
      	  TNXSM=T_MSIS(2)
      	  IF(TEXSM.LT.TNXSM) TEXSM=TNXSM
          atiii10=tiv(i10-1)
      	  IF(atiii10.GT.TEXSM) atiii10=TEXSM
          IF(atiii10.LT.TNXSM) atiii10=TNXSM
          atii(i10)=atiii10
          mm(i10-1)=(atii(i10)-atii(i10-1))/(xsm(i10)-xsm(i10-1))
          enddo
c        mm(1)=(TIV(1)-TNHS)/(xsm(2)-xsm(1))
c        mm(2)=(TIV(2)-TIV(1))/(xsm(3)-xsm(2))
c        mm(3)=(TIV(3)-TIV(2))/(xsm(4)-xsm(3))
c        mm(4)=(TIV(4)-TIV(3))/(xsm(5)-xsm(4))

c XTETI is altitude where Te=Ti
        XTETI=AHH(numap)
        XTTS=500.
        X=500.
2397    X=X+XTTS
        IF(X.GT.AHH(numap)) GOTO 240
        TEX=BOOKER1(X,numap2,ATE1,AHH,STTE,DTE)	
        TIX=BOOKER1(X,numti2,TNHS,XSM,MM,DTI)
        IF(TIX.LT.TEX) GOTO 2397
        X=X-XTTS
        XTTS=XTTS/10.
        IF(XTTS.GT.0.1) GOTO 2397
        XTETI=X+XTTS*5.

      else
c
c Bil-1981 model 
        numti=3
        numti2=numti-2

c Ti(430km) during daytime from AEROS data
      	XSM(2)=430.0
      	Z1=EXP(-0.09*MLAT)
      	Z2=Z1+1.
      	TID2=1240.0-1400.0*Z1/(Z2*Z2)

c Ti(430km) during nighttime from AEROS data
      	Z1=ABSMLT
      	Z2=Z1*(0.47+Z1*0.024)*UMR
      	Z3=COS(Z2)
      	TIN2=1200.0-300.0*SIGN(1.0,Z3)*SQRT(ABS(Z3))

c Ti(430km) for specified time using HPOL
      	TI2=TIN2  
      	IF(TID2.GT.TIN2) 
     &     TI2=HPOL(HOUR,TID2,TIN2,SAX300,SUX300,1.,1.)      

c Tn < Ti < Te enforced at XSM(2)=430 km 
        TEXSM2=BOOKER1(xsm(2),numap2,ATE1,AHH,STTE,DTE)	
      	CALL GTD7(IYD,SEC,xsm(2),LATI,LONGI,hour,F10781OBS,
     &        F107YOBS,IAPO,0,D_MSIS,T_MSIS)
      	TNXSM2=T_MSIS(2)
      	IF(TEXSM2.LT.TNXSM2) TEXSM2=TNXSM2
      	IF(TI2.GT.TEXSM2) TI2=TEXSM2
      	IF(TI2.LT.TNXSM2) TI2=TNXSM2
        atii(2)=TI2

c First segment is from 200km to 430km
      	MM(1)=(atii(2)-atii(1))/(XSM(2)-XSM(1))

c Gradient for second segment 
      	MM(2)=HPOL(HOUR,3.0,0.0,SAX300,SUX300,1.,1.)

c XTETI is altitude where Te=Ti
        XTETI=AHH(numap)
        XTTS=500.
        X=500.
2390    X=X+XTTS
        IF(X.GE.AHH(numap)) GOTO 240
        TEX=BOOKER1(X,numap2,ATE1,AHH,STTE,DTE)
        TIX=BOOKER1(X,numti2,TNHS,XSM,MM,DTI)
        IF(TIX.LT.TEX) GOTO 2390
        X=X-XTTS
        XTTS=XTTS/10.
        IF(XTTS.GT.0.1) GOTO 2390
        XTETI=X+XTTS*5.
      	XSM(3)=xteti
        atii(3)=tix

      endif

C
C CALCULATION OF ION DENSITY PARAMETER..................
C

240   IF(NOION) GOTO 141
      HNIA=75.
      if(RBTT) HNIA=80.
      HNIE=2000.
 
C*****************************************************************
C CALCULATION FOR THE REQUIRED HEIGHT RANGE.......................
C In the absence of an F1 layer hmf1=hz because hmf1 is used in XE
C
141     xhmf1=hmf1
        IF(hmf1.le.0.0) HMF1=HZ

        height=heibeg
        kk=1
	xinv=0.0

300   CALL SOCO(daynr,HOUR,LATI,LONGI,height,SUNDEC,XHI,SAX,SUX)

c
c electron density ELEDE in m-3 in outf(1,*)
c
      elede=-1.0  
      IF(NODEN) GOTO 330
      IF(HEIGHT.LT.HNEA) then
        if((jfirsta.lt.1).and.mess) write(konsol,2920) hnea
2920    format(/'At heights below ',F4.1,' km Ne is ',  
     &    'negligible and IRI values are not calculated.')
        jfirsta=2
        GOTO 330
        endif
      IF(HEIGHT.GT.HNEE) then
        if((jfirsta.lt.1).and.mess) write(konsol,2720) int(hnee)
2720    format(/'The upper boundary for Ne profiles is ',
     &    I5,' km.')
        jfirsta=2
        GOTO 330
        endif
      IF((HEIGHT.GT.2000).and.(jfirste.lt.1)) then
	if(mess) write(konsol,2919) int(hnee)
2919    format(/'At heights above 2000 km a first-order ',
     &    'plasmaspheric extension is applied up to ',I5,'km.')
        jfirste=2
        endif

      IF(LAYVER) THEN
          ELEDE=-9.
          IF(IIQU.LT.2) ELEDE=
     &         XEN(HEIGHT,HMF2,NMF2S,HME,4,HXL,SCL,AMP)
          outf(1,kk)=elede
          goto 330
          endif
c
c Calculates height-dependent correction factors for tcor1 and tcor2
c for topside options 2001cor and COR2
c		 
      IF(itopn.eq.1.or.itopn.eq.3) then             
	tcor1 = 0.0
	IF(height.ge.hcor1) tcor1 = 
     &    BOOKER(height,6,pah,palogne,dplas)

        tcor2 = 0.0
        if(itopn.eq.3.and.height.gt.hmF2) then
c        if(itopn.eq.3.and.height.ge.hcor2) then
          CALL SOCO(daynr,HOUR,LATI,LONGI,height,SUC,zxz,sap,sup)
c          tcor2 = 0.0
c          if(height.ge.hcor2) tcor2 = 
c     &      TCOR2CAL(height,hour,modip,pf107,sap,sup)
          tcor2 = TCOR2CAL(height,hour,modip,pf107,sap,sup)
          if(height.lt.hcor2) tcor2=(exp((height-hmF2)/shc)-1)*tcor2 
          endif
        endif
		
      ELEDE=XE_1(HEIGHT)
      if(height.gt.hmF2.and.elede.gt.nmf2s) elede=nmf2s
	  
c
c FIRI D region 
c
      if(.not.dreg.and.height.le.140.) then
            elede=-1.
            call F00(HEIGHT,LATI,DAYNR,XHI,F107D,EDENS,IERROR)
            if(ierror.eq.0.or.ierror.eq.2) elede=edens
            endif
      OUTF(1,kk)=ELEDE
	  
c
c plasma-frequency divided by gyro-frequency
c	
  
      pf_gf = 3.2045e-3 * sqrt(elede/1.e6)/babs	
      OUTF(15,kk)=pf_gf

c
c plasma temperatures in Kelvin
c

330   IF(NOTEM.and.(NOION.or..not.RBTT)) GOTO 7108
      IF((HEIGHT.GT.HTE).OR.(HEIGHT.LT.HTA)) GOTO 7108
      CALL GTD7(IYD,SEC,HEIGHT,LATI,LONGI,HOUR,F10781OBS,
     &        F107YOBS,IAPO,0,D_MSIS,T_MSIS)
      TNH=T_MSIS(2)
      TEH=TNH
      if(HEIGHT.GT.HEQUI) 
     &  TEH=BOOKER1(HEIGHT,numap2,ATE1,AHH,STTE,DTE)
      if(HEIGHT.LE.HS) THEN
         TIH=TNH
      ELSE IF(HEIGHT.LE.XTETI) THEN
         TIH=BOOKER1(HEIGHT,numti2,TNHS,XSM,MM,DTI) 
      ELSE    
         TIH=TEH
      ENDIF     
c
c Tn < Ti < Te enforced
c  
      if(TIH.lt.TNH) TIH=TNH
      if(TEH.lt.TNH) TEH=TNH
      if(TIH.gt.TEH) TIH=TEH

      OUTF(2,kk)=TNH
      OUTF(3,kk)=TIH
      OUTF(4,kk)=TEH

c
c ion composition
c
7108  IF(NOION) GOTO 7118
      IF((HEIGHT.GT.HNIE).OR.(HEIGHT.LT.HNIA)) GOTO 7118
      ROX=-1.
      RHX=-1.
      RNX=-1.
      RHEX=-1.
      RNOX=-1.
      RO2X=-1.
      RCLUST=-1.
      if(RBTT) then
        if(height.ge.300.) then
c Triskova-Truhlik-Smilauer-2003 model
          call CALION(invdip,xmlt,height,daynr,pf107obs,
     &      xic_O,xic_H,xic_He,xic_N)
          rox=xic_O*100.
          rhx=xic_H*100.
          rnx=xic_N*100.
          rhex=xic_He*100.
          rnox=0.
          ro2x=0.
        else
c Richards-Bilitza-Voglozin-2010 IDC model
          CALL GTD7(IYD,SEC,height,lati,longi,HOUR,f10781obs,
     &      f107yobs,IAPO,48,D_MSIS,T_MSIS)
          XN4S=0.5*D_MSIS(8)
          EDENS=ELEDE/1.e6
c         jprint=1      ! turns messages on 
          jprint=0
          Den_NO = 0.0
          rn = 0.0
          CALL CHEMION(jprint,height,F107YOBS,F10781OBS,TEH,TIH,
     &      TNH,D_MSIS(2),D_MSIS(4),D_MSIS(3),D_MSIS(1),
     &      D_MSIS(7),-1.0,XN4S,EDENS,-1.0,xhi,ro,ro2,rno,rn2,
     &      rn,Den_NO,Den_N2D,INEWT)                              
          if(INEWT.gt.0) then
            sumion = edens/100.
            rox=ro/sumion
            rhx=0.
            rhex=0.
            rnx=rn/sumion
            rnox=rno/sumion
            ro2x=ro2/sumion
            endif
        endif
      else
c Danilov-Smirnova-1995 model below 300km and Danilov-Yaichnikov-1985 model above
        call iondani(iday,iseamon,height,xhi,lati,f107365,dion)
        ROX=DION(1)
        RHX=DION(2)
        RNX=DION(3)
        RHEX=DION(4)
        RNOX=DION(5)
        RO2X=DION(6)
        RCLUST=DION(7)
      endif
c
c ion densities are given in percent of total electron density;
c

      if(jf(22)) then 
            xnorm=1
      else
            xnorm=elede/100.
      endif
      OUTF(5,kk)=ROX*xnorm
      OUTF(6,kk)=RHX*xnorm
      OUTF(7,kk)=RHEX*xnorm
      OUTF(8,kk)=RO2X*xnorm
      OUTF(9,kk)=RNOX*xnorm
      OUTF(10,kk)=RCLUST*xnorm
      OUTF(11,kk)=RNX*xnorm

7118  height=height+heistp
      kk=kk+1
      if(kk.le.numhei) goto 300

C
C END OF PARAMETER COMPUTATION LOOP 
C

c
c D region special: densities for 11 heights (60,65,70,..,110km)
c outf(14,1:11)=IRI-07, outf(14,12:22)=FIRI, 
c outf(14,23:33)= Danilov et al.(1995) with SW=0,WA=0 
c outf(14,34:44)= with SW=0.5,WA=0, 
c outf(14,45:55)= with SW=1,WA=0,  
c outf(14,56:66)= with SW=0,WA=0.5, 
c outf(14,67:77)= with SW=0,WA=1,  
c

      if(.not.dreg) then
        do ii=1,11
          Htemp=55+ii*5  
          outf(14,ii)=-1.     
          if(Htemp.ge.65.) outf(14,ii)=XE6(Htemp)     
          outf(14,11+ii)=-1.
          call F00(Htemp,LATI,DAYNR,XHI1,F107D,EDENS,IERROR)
          if(ierror.eq.0.or.ierror.eq.2) outf(14,11+ii)=edens
          outf(14,22+ii)=ddens(1,ii)      
          outf(14,33+ii)=ddens(2,ii)      
          outf(14,44+ii)=ddens(3,ii)      
          outf(14,55+ii)=ddens(4,ii)      
          outf(14,66+ii)=ddens(5,ii)      
          enddo
        endif

c
c equatorial vertical ion drift
c

      drift=-1.
      if(jf(21).and.abs(magbr).lt.25.0) then
        CALL vfjmodelrocstart(FJM)
        CALL vfjmodelrocinit(F107D,DAYNR,JSEA,JF107)
        CALL vfjmodelroc(FJM,HOUR,LONGI,JSEA,JF107,DRIFT)
        endif

c
c spread-F occurrence probability (range: 0-1)
c
      spreadf=-1.
      if(.not.jf(28)) goto 1937
      if(hour.gt.7.25.and.hour.lt.17.75) goto 1937
      if(abs(lati).gt.25.0) goto 1937
        spfhour=hour
        daynr1=daynr
        month1=month
        if(hour.lt.12.0) then
          spfhour=hour+24.0
          daynr1=daynr-1
          if(daynr1.lt.1) then
             month1=month-1
             if(month1.lt.1) then
                month1=12
                daynr1=idayy
              else
                daynr1=mmo(month1)
              endif
            endif  
          endif
        call spreadf_brazil(daynr1,month1,f107d,lati,osfbr)
        ispf=int((spfhour-17.75)/0.5)+1
        if(ispf.gt.0.and.ispf.lt.26) spreadf=osfbr(ispf)
1937   continue

C
C ADDITIONAL PARAMETER FIELD OARR: angles are given in degrees,
C times in decimal hours, altitudes in km, densities in m-3, and
C temperatures in K      
C
      IF(NODEN) GOTO 6192
      OARR(1)=NMF2S	! F2-peak density in m-3
      OARR(2)=HMF2	! F2-peak height in km
      OARR(3)=NMF1	! F1-peak density in m-3
      OARR(4)=XHMF1	! F1-peak height in km
      OARR(5)=NMES	! E-peak density in m-3
      OARR(6)=HME	! E-peak height in km
      OARR(7)=NMD	! density in m-3 of D-region inflection point 
      OARR(8)=HMD	! height in km of D-region inflection point
      OARR(9)=HHALF	! height used by Gulyaeva B0 model
      OARR(10)=B0	! bottomside thickness parameter in km
      OARR(11)=VNER	! density in m-3 at E-valley bottom 
      OARR(12)=HEF	! height in km of E-valley top (Ne(HEF)=NmE)
6192    IF(NOTEM.and.(NOION.or..not.RBTT)) GOTO 6092
      OARR(13)=ATE(2)	! electron temperature Te in K at AHH(2)
      OARR(14)=AHH(2)	! intermediate height between 120km and 300/350km
      OARR(15)=ATE(3)	! Te at 300km/350km for BIL-1995/TBT2012+SA model
      OARR(16)=ATE(4)	! Te at 400km/550km for BIL-1995/TBT2012+SA model
      OARR(17)=ATE(5)	! Te at 600km/850km for BIL-1995/TBT2012+SA model
      OARR(18)=ATE(6)	! Te at 1400km/1400km for BIL-1995/TBT2012+SA model
      OARR(19)=ATE(7)	! Te at 3000km/2000km for BIL-1995/TBT2012+SA model
      OARR(20)=ATE(1)	! Te at 120km = neutral temperature from CIRA
      OARR(21)=TI1	! ion temperature in K at 430km
      OARR(22)=XTETI	! altitude where Te=Ti
6092  OARR(23)=XHI3	! solar zenith angle at 200 km
      OARR(24)=SD300	! sun declination 
      OARR(25)=DIP	! IGRF magnetic inclination (dip)
      OARR(26)=MAGBR	! IGRF dip latitude
      OARR(27)=MODIP	! modified dip latitude
      OARR(28)=LATI	! geographic latitude		
      OARR(29)=SAX200	! time of sunrise at 200 km
      OARR(30)=SUX200	! time of sunset at 200 km
      OARR(31)=SEASON	! =1 spring, 2= summer .. 
c SEASON assumes equal length seasons (92 days) with spring 
c (SEASON=1) starting at day-of-year=45
      OARR(32)=LONGI	! geographic longitude		
      OARR(33)=rssn	! 12-month running mean of sunspot number
      OARR(34)=COV	! 12-month running mean of F10.7
      OARR(35)=B1	! Bottomside shape parameter
      OARR(36)=xm3000	! Propagation factor M(3000)F2
C OARR(37) used for TEC and 38 for TEC-top
      OARR(39)=gind	! 12-month running mean of IG index
      OARR(40)=f1pb	! probability for an F1 layer
      OARR(41)=f107d	! daily solar radio flux at 10.7cm:F10.7
      OARR(42)=c1	! shape parameter for F1 layer
      OARR(43)=daynr	! day of year
      OARR(44)=drift	! vertical ion drift at equator in m/s
      OARR(45)=stormcorr! ratio foF2_storm/foF2_quiet
      OARR(46)=f10781	! 81-day average of F10.7
      OARR(47)=estormcor! ratio foE_storm/foE_quiet
      OARR(48)=spreadf	! probability of spread-F occurrence
      OARR(49)=MLAT	! IGRF magnetic latitude
      OARR(50)=MLONG	! IGRF magnetic longitude
      OARR(51)=index_3h_ap*1.0	! ap index for current UT
      OARR(52)=IAP_daily*1.0	! daily ap index
      OARR(53)=invdip_old	! invariant dip latitude
      OARR(54)=XMLT	! Magnetic Local Time
C Please check subroutine GEOCGM01 in file IGRF.FOR for more
C information on the Corrected Geomagnetic (CGM) coordinates.
C CGM coordinates are only calculated if you select
C AURORAL BOUNDARIES <on>      
      OARR(55)=cgm_lat	! Corrected Geomagnetic (CGM) latitude
      OARR(56)=cgm_lon	! Corrected Geomagnetic (CGM) longitude
      OARR(57)=cgm_mlt	! Magnetic Local Time for CGM coord.
      OARR(58)=cgmlat   ! CGM latitude of equatorward boundary
c include only every second auroral boundary point (MLT=0,1,2..23)
      jjj=58
      do iii=1,47,2	! CGM latitude at MLT=0,1,2 ...23
         jjj=jjj+1 
         oarr(jjj)=ab_mlat(iii)
         enddo
      OARR(83)=xkp	! Kp at the time specified by the user
      OARR(84)=dec	! magnetic declination in degrees
      OARR(85)=fl	! L-value
      OARR(86)=dimo	! Earth's dipole moment
      OARR(87)=SAX300	! sunrise at 300km in decimal hours
      OARR(88)=SUX300	! sunset at 300km in decimal hours		
      OARR(89)=HNEA	! lower boundary in km of Ne valid range	
      OARR(90)=HNEE	! upper boundary in km of Ne valid range	
      OARR(91)=ESP	! sporadic-E occurence probability (0-1)		
      OARR(92)=XIBP	! plasma bubble occ. probability (0-1)		

3330  CONTINUE

       icalls=icalls+1

       RETURN
       END
c
c
        subroutine iri_web(jmag,jf,alati,along,iyyyy,mmdd,iut,dhour,
     &    height,h_tec_min,h_tec_max,ivar,vbeg,vend,vstp,a,b)
c-----------------------------------------------------------------------        
c changes:
c       11/16/99 jf(30) instead of jf(17)
c       10/31/08 outf, a, b (100 -> 500)
c
c-----------------------------------------------------------------------        
c input:   jmag,jf(50),alati,along,iyyyy,mmdd,iut,dhour (see IRI_SUB)          
c          height  height in km
c          h_tec_min lower boundary in km for TEC integral
c          h_tec_max upper boundary in km (=0 TEC not computed) 
c          ivar    parameter that is varying
c                  =1      altitude
c                  =2,3    latitude,longitude
c                  =4,5,6  year,month,day
c                  =7      day of year
c                  =8      hour (UT or LT)
c          vbeg,vend,vstp  variable range (begin,end,step)
c output:  a(20,1000)      contains outf(20) output parameters for all 
c                               maximaly 1000 variable steps 
c          b(100,1000)     contains oar(100) output parameters for all 
c                               maximaly 1000 variable steps 
c
c          numstp  number of steps; maximal 1000
c-----------------------------------------------------------------------        
        dimension   outf(20,1000),oar(100),oarr(100),a(20,1000)
        dimension   xvar(8),b(100,1000)
        logical     jf(50)

        nummax=1000
        numstp=int((vend-vbeg)/vstp)+1
        if(numstp.gt.nummax) numstp=nummax

        do 6249 i=1,100
6249      oar(i)=b(i,1) 

        if(MMDD.lt.0) then
           IDOY=-MMDD
           call MODA(1,iyyyy,MONTH,IDAY,IDOY,nrdaym)
           MMDD=MONTH*100+IDAY
        else
           MONTH=MMDD/100
           IDAY=MMDD-MONTH*100
           call MODA(0,iyyyy,MONTH,IDAY,IDOY,nrdaym)
        endif

        if(ivar.eq.1) then
            do 1249 i=1,100
1249            oarr(i)=oar(i) 
            xhour=dhour+iut*25.
            call IRI_SUB(JF,JMAG,ALATI,ALONG,IYYYY,MMDD,
     &                  XHOUR,VBEG,VEND,VSTP,a,OARR)
            if(h_tec_max.gt.50.) then 
               call IRITEC(ALATI,ALONG,jmag,jf,iyyyy,mmdd,
     &           xhour,h_tec_min,h_tec_max,1.0,oarr,tecbo,tecto)
               oarr(37) = tecbo + tecto
               oarr(38) = tecto / oarr(37) * 100
               endif
            do 1111 i=1,100
1111           b(i,1)=oarr(i)
            return
            endif

        if(height.le.0.0) height=100
        xvar(2)=alati
        xvar(3)=along
        xvar(4)=iyyyy*1.0
        xvar(5)=month*1.0
        xvar(6)=iday*1.0
        xvar(7)=idoy
        xvar(8)=dhour

        xvar(ivar)=vbeg

        alati=xvar(2)
        along=xvar(3)
        iyyyy=int(xvar(4))
        if(ivar.eq.7) then
           mmdd=-int(vbeg)
        else
           mmdd=int(xvar(5)*100+xvar(6))
        endif
        xhour=xvar(8)+iut*25.

        do 1 i=1,numstp	

          do 1349 iii=1,100
1349        oarr(iii)=b(iii,i)

          call IRI_SUB(JF,JMAG,ALATI,ALONG,IYYYY,MMDD,
     &          XHOUR,HEIGHT,HEIGHT,1.,OUTF,OARR)
          if(h_tec_max.gt.50.) then
            call IRITEC(ALATI,ALONG,jmag,jf,iyyyy,mmdd,
     &        xhour,h_tec_min,h_tec_max,1.0,oarr,tecbo,tecto)
	    oarr(37) = tecbo + tecto
            oarr(38) = tecto / oarr(37) * 100
            endif
			
          do 2 ii=1,20
2           a(ii,i)=outf(ii,1)
          do 2222 ii=1,100
2222        b(ii,i)=oarr(ii)

          xvar(ivar)=xvar(ivar)+vstp

          alati=xvar(2)
          along=xvar(3)
          iyyyy=int(xvar(4))
          if(ivar.eq.7) then
            mmdd=-int(xvar(7))
          else
            mmdd=int(xvar(5)*100+xvar(6))
          endif
          xhour=xvar(8)+iut*25.
1       continue

        return
        end

