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. 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
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.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.34 09/15/25 Declare DAYNR1 as integer
C 2020.34 09/15/25 iri_web: dhour->xhour(=xvar(8)+iut*25.)
C 2020.34 09/15/25 iri_web: including MODA and int for 
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 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   only if foF1>0 and not NIGHT        t
C   20    standard F1            standard F1 plus L condition        t
C (19,20) = (t,t) f1-prob, (t,f) f1-prob-L, (f,t) old F1, (f,f) no F1
C   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 (TBT-2012)           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    IRI01-topside          new options as def. by JF(30)       t
C   30    IRI01-topside corr.    NeQuick topside model   	 false 
C (29,30) = (t,t) IRIold, (f,t) IRIcor, (f,f) NeQuick, (t,f) IRIcor2
C   31    B0,B1 ABT-2009	 B0 Gulyaeva-1987 h0.5               t   
C (4,31) = (t,t) Bil-00, (f,t) ABT-09, (f,f) Gul-87, (t,f) not used
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    turn WRITEs off in IRIFLIP   turn WRITEs on                t
C   39    hmF2 (M3000F2)         new models                      false
C   40    hmF2 AMTB-model        Shubin-COSMIC model             false
C (39,40) = (t,t) hmF2-old, (f,t) AMTB, (f,f) Shubin, (t,f) not used
C   41    Use COV=F10.7_365      COV=f(IG12) (IRI before Oct 2015)   t
C   42    Te with PF10.7 dep.	 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   47    CGM computation on 	 CGM computation off             false
C   48    Ti  Tru-2021           Bil-1981                            t
C   49    Plasmasphere: Ozhogin  Gallagher 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 %   
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
      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
      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),MM(6),DTI(6),AHH(7),
     &  STTE(6),DTE(5),ATE(7),TEA(4),XNAR(2),param(2),OARR(100),
     &  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(5),sdteva(5),tiv(4),sigtv(4),FM(59,25,4,11)

      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,OLD79,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
     &   /iounit/konsol,mess     /CSW/SW(25),ISW,SWC(25)
     &   /QTOP/Y05,H05TOP,QF,XNETOP,XM3000,HHALF,TAU 
     &   /cotec/hnea,hpp
      EXTERNAL          XE1,XE2,XE3_1,XE4_1,XE5,XE6,FMODIP

      DATA icalls/0/, dplas/100,150,10,10/,jfirsta,jfirste/0,0/
      DATA DTE/5.,5.,10.,20.,20./, DTI/5.,5.,10.,20.,20.,20./

        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,91,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 NEW-GUL------------------------------
c         Y05=.6931473
c         QF=1.
c         h05top=0.
C NEW-GUL------------------------------

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))
      OLD79=(.not.jf(7))
      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 TE-NE MODEL OPTION ..............................................
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(OLD79) 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) 
             goto 2889
             endif
          if(URSIF2) then
             write(konsol,9016) 
          else
             write(konsol,9017) 
          endif
          if(jf(40)) then
             write(konsol,9717) 
          else
             write(konsol,9716) 
          endif
          if(HMF2IN) write(konsol,9018) 
          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(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) 
          else
            if(jf(42)) then
              write(konsol,9034)
            else
              write(konsol,9534)
            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')
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')
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: OTSR2012 for plasmasphere')
9228    format('Ne: GCC-2000 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 is used')
9034    format('Te: TBT-2012 with PF10.7 depend. is used')
9534    format('Te: TBT-2012 without PF10.7 depend. is used')
9333    format('Ti: TBKS2021 model is used')
9335    format('Ti: Bil-1981 model is used')
4031    format('AB: Auroral boundary model on')
4032    format('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
C NEW-GUL------------------------------
c        sday=daynr/idayy*360.			 
C NEW-GUL------------------------------
        seaday=daynr
        iseamon=month
        IF(LATI.GE.0.0) GOTO 5592
          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
C NEW-GUL------------------------------
c            sday=sday+180.						
c            if (sday.gt.360.) sday=sday-360.	
C NEW-GUL------------------------------

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

5592    continue
        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        rlimit=gind
c        COVSAT=63.75+rlimit*(0.728+rlimit*0.00089)

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 (if PF107Y=1 then PF107OBS)
c
        call apf(isdate,hourut,indap)

        if(jf(45)) then
           xKpsum=0
           do i=5,12
              xxkp=ckp(indap(i))
              xKpsum=xKpsum+xxkp
              enddo
           isoma=1
           if(jf(46)) isoma=0
           xMLLOO = MLONG
           if(xMLLOO.eq.360.0) xMLLOO = 0.0
           CALL ESPROB(isoma,isoma,INVDIP_OLD_110,xMLLOO,XMLT,
     &       DAYNR,F107_365,xKpsum,ESP)
           endif
C
C F2 peak critical frequency foF2, density NmF2, and height hmF2
C
C READ CCIR AND URSI COEFFICIENT SET FOR CHOSEN MONTH 
C
      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 20 I=1,76
        DO 20 J=1,13
              K=J+13*(I-1)
              FF0N(K)=F2N(J,I,1)*RR1N+F2N(J,I,2)*RR2N
20            FF0(K)=F2(J,I,1)*RR1+F2(J,I,2)*RR2

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

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    
           icoord=1
           kut=int(hourut)
           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
        COS2=COS(MLAT*UMR)
        COS2=COS2*COS2
        FLU=(COVSAT-40.0)/30.0
c
c option to use unlimited F10.7M for the topside
c previously: IF(OLD79) ETA1=-0.0070305*COS2
c
        IF(OLD79) 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)
c
c NEW-GUL--------------------------------
c
c Correction term for topside (Gulyaeva)
c
c      if(itopn.eq.3) then
c          hei05=0.
c          CALL TOPH05(COV,MLAT,HOUR,HMF2,HEI05,SDAY)
c          h05top=hei05
c          xnetop=XE_1(H05TOP)
c          endif
c NEW-GUL--------------------------------

c
c NeQuick topside parameters 
C Use CCIR-M3000F2 even if user_hmF2 or user_M(3000)F2
c
      if (itopn.eq.2) then
         fo2=foF2s
         if(jf(37)) fo2=foF2
         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
      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 (fixpt:5000,10000,20000,30000km)
c

        if(jf(50)) then
          hmid=5000.0
          hpp =10000.0
          hppo=20000.0
          hpt =30000.0
          xlmid=(1.0+hmid/ERA)/cosmag2
          xlpp =(1.0+hpp/ERA)/cosmag2
          xlppo=(1.0+hppo/ERA)/cosmag2
          xlpt =(1.0+hpt/ERA)/cosmag2
        else

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

          xlpp =5.0
          xlppo=5.1		
          hpp =ERA*(xlpp*cosmag2-1.0)
          hppo=ERA*(xlppo*cosmag2-1.0)
          hmid=(hpp+hmF2+1500.0)/2.0
          hpt=2*hpp-hmid
          xlmid=(1.0+hmid/ERA)/cosmag2
          xlpt =(1.0+hpt/ERA)/cosmag2
        endif

c		
c plasmasphere models (Gallagher-2000, Ohzogin-2012)
c		

        if(jf(49)) then
          xnepp =ohzden(xlpp,ppmlat)
          xnemid=ohzden(xlmid,ppmlat)
          xneppo=ohzden(xlppo,ppmlat)
          xnept =ohzden(xlpt,ppmlat)
        else
          xnepp =gallden(xlpp,daynr,rssn)
          xnemid=gallden(xlmid,daynr,rssn)
          xneppo=gallden(xlppo,daynr,rssn)
          xnept =gallden(xlpt,daynr,rssn)
        endif

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
          xneppo=xnepp/10.0
c          xneppo=caadenet(xlppo,xmlt)
          xnept=caadenet(xlpt,xmlt)
          endif

c		
c Booker parameters for plasmaspheric extension
c		
        PAH(1)=hcor1
        PAH(2)=hmF2 + 1500.0
        PAH(3)=hmid		
        PAH(4)=hpp
        PAH(5)=hppo
        PAH(6)=hpt      
        PALOGNE(1)=0.0
        PALOGNE(2)=r2*alg10
        PALOGNE(3)=log(xnemid/nmf2s)
        PALOGNE(4)=log(xnepp/nmf2s)
        PALOGNE(5)=log(xneppo/nmf2s)
        PALOGNE(6)=log(xnept/nmf2s)

        tcor1=0.0        
        CALL SOCO(daynr,HOUR,LATI,LONGI,hmid,SUC,zxz,sap,sup)
        TCOR2 = TCOR2CAL(hmid,hour,modip,pf107,sap,sup)
        znemid = log(XE_1(hmid)/nmf2s)
        palogne(3)=palogne(3)-znemid
        CALL SOCO(daynr,HOUR,LATI,LONGI,hpp,SUC,zxz,sap,sup)
        TCOR2 = TCOR2CAL(hpp,hour,modip,pf107,sap,sup)
        znepp = log(XE_1(hpp)/nmf2s)
        palogne(4)=palogne(4)-znepp
        CALL SOCO(daynr,HOUR,LATI,LONGI,hppo,SUC,zxz,sap,sup)
        TCOR2 = TCOR2CAL(hppo,hour,modip,pf107,sap,sup)
        zneppo = log(XE_1(hppo)/nmf2s)
        palogne(5)=palogne(5)-zneppo
        CALL SOCO(daynr,HOUR,LATI,LONGI,hpt,SUC,zxz,sap,sup)
        TCOR2 = TCOR2CAL(hpt,hour,modip,pf107,sap,sup)
        znept = log(XE_1(hpt)/nmf2s)
        palogne(6)=palogne(6)-znept
      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:
c Te(300km), Te(400km) from AE-C, Te(1400km) from ISIS-2, and 
C Te(3000km) from ISIS-1 (Brace and Theis, 1981)
              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 TBT-2012 model:
c Te at fixed heights 350, 550, 850, 1400, and 2000 km with and
c without solar activity effects included (Truhlik et al., 2012)
              AHH(3)=350.
              AHH(4)=550.
              AHH(5)=850.
              AHH(6)=1400.
              AHH(7)=2000.
c isa for solar activity correction: isa=0 sol activity corr off
              isa=0
              if(jf(42)) isa=1
              call elteik(isa,invdip_old,xmlt,daynr,
     &              pf107obs,teva,sdteva)
              do ijk=3,7
                 ate(ijk)=teva(ijk-2)
                 enddo
          endif

c Option to use Te = f(Ne) relation at ahh(3), ahh(4)

          IF(TENEOP) THEN
              DO 3395 I=1,2
3395              IF(TECON(I)) ATE(I+2)=TEDE(AHH(I+2),XNAR(I),-COV)
              endif

c Te corrected and Te > Tn 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 1901 I=2,6
         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))
         ATE(I)=ATE(I)-(STTE2-STTE1)*DTE(I-1)*ALOG2
1901     STTE1=STTE2

c Te gradients STTE are computed for each segment

      DO 1902 I=1,6
1902     STTE(I)=(ATE(I+1)-ATE(I))/(AHH(I+1)-AHH(I))
      ATE1=ATE(1)
887   CONTINUE

      HTE = AHH(7)

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)

      if(jf(48)) then

c Tru-2021 model: TIV(1:4) ion temperature at 350,430,600,850km
c JSOL=1 solar activity dependence included
c JSOL=0 solar activity dependence excluded
 
	JSOL=1
	xsm(1)=200
	xsm(2)=350
	xsm(3)=430
	xsm(4)=600
	xsm(5)=850
      	CALL IONTIF(JSOL,INVDIP_OLD_600,XMLT,DAYNR,PF107OBS,TIV,SIGTV)      
c Tn < Ti < Te enforced at xsm(2:5)
        do 2391 i10=2,5
          HTIX=XSM(i10)
          TEXSM=BOOKER1(HTIX,5,ATE1,AHH,STTE,DTE)	
      	  CALL GTD7(IYD,SECNI,XSM(i10),LATI,LONGI,0.0,
     &        F10781OBS,F107YOBS,IAPO,0,D_MSIS,T_MSIS)
      	  TNXSM=T_MSIS(2)
      	  IF(TEXSM.LT.TNXSM) TEXSM=TNXSM
      	  IF(TIV(i10-1).GT.TEXSM) TIV(i10-1)=TEXSM
2391      IF(TIV(i10-1).LT.TNXSM) TIV(i10-1)=TNXSM

        mm(1)=(TIV(1)-TNHS)/(xsm(2)-xsm(1))
        mm(2)=(TIV(2)-TIV(1))/(xsm(3)-xsm(2))
        mm(3)=(TIV(3)-TIV(2))/(xsm(4)-xsm(3))
        mm(4)=(TIV(4)-TIV(3))/(xsm(5)-xsm(4))
        MXSM=3

c XTETI is altitude where Te=Ti
        XTTS=500.
        X=500.
2397    X=X+XTTS
        IF(X.GE.AHH(7)) GOTO 240
        TEX=BOOKER1(X,5,ATE1,AHH,STTE,DTE)	
        TIX=BOOKER1(X,MXSM,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.
        if(xteti.gt.xsm(5)+100) then
           xsm(6)=xteti
           TEX5=BOOKER1(XTETI,5,ATE1,AHH,STTE,DTE)
           mm(5)=(TEX5-TIV(4))/(xsm(6)-xsm(5))
           MXSM=4
        else if(xteti.gt.xsm(4)) then
           xsm(5)=xteti           
           TEX5=BOOKER1(XTETI,5,ATE1,AHH,STTE,DTE)
           mm(4)=(TEX5-TIV(3))/(xsm(5)-xsm(4))
        else 
           xsm(4)=xteti
           TEX5=BOOKER1(XTETI,5,ATE1,AHH,STTE,DTE)
           mm(3)=(TEX5-TIV(2))/(xsm(4)-xsm(3))
           MXSM=2
 	   endif        
      else
c
c Bil-1981 model: 
c Ti(430km) during daytime from AEROS data
c
      	XSM1=430.0
      	XSM(2)=XSM1
      	Z1=EXP(-0.09*MLAT)
      	Z2=Z1+1.
      	TID1 = 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)
      	TIN1=1200.0-300.0*SIGN(1.0,Z3)*SQRT(ABS(Z3))

c Ti(430km) for specified time using HPOL
      	TI1=TIN1  
      	IF(TID1.GT.TIN1) TI1=HPOL(HOUR,TID1,TIN1,SAX300,SUX300,1.,1.)
      
c Tn < Ti < Te enforced at 430 km 
        TEN1=BOOKER1(XSM1,5,ATE1,AHH,STTE,DTE)	
      	CALL GTD7(IYD,SECNI,XSM1,LATI,LONGI,0.0,F10781OBS,
     &        F107YOBS,IAPO,0,D_MSIS,T_MSIS)
      	TNN1=T_MSIS(2)
      	IF(TEN1.LT.TNN1) TEN1=TNN1
      	IF(TI1.GT.TEN1) TI1=TEN1
      	IF(TI1.LT.TNN1) TI1=TNN1

c First segment is from 200km to 430km
      	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)
      	MM(1)=(TI1-TNHS)/(XSM1-HS)
      	MXSM=1

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

c XTETI is altitude where Te=Ti
        XTTS=500.
        X=500.
2390    X=X+XTTS
        IF(X.GE.AHH(7)) GOTO 240
        TEX=BOOKER1(X,5,ATE1,AHH,STTE,DTE)	
        TIX=BOOKER1(X,MXSM,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.

c Ti=Te above XTETI 
c
c        MXSM=2
c        MM(3)=STTE(6)
c        XSM(3)=XTETI
c        IF(XTETI.GT.AHH(6)) GOTO 240
c        MXSM=3
c        MM(3)=STTE(5)
c        MM(4)=STTE(6)
c        XSM(4)=AHH(6)
c        IF(XTETI.GT.AHH(5)) GOTO 240
c        MXSM=4
c        MM(3)=STTE(4)
c        MM(4)=STTE(5)
c        MM(5)=STTE(6)
c        XSM(4)=AHH(5)
c        XSM(5)=AHH(6)
      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
C In the absence of an F1 layer hmf1=hz since 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 no longer calculating invdip for each height
c       call igrf_sub(lati,longi,ryear,height,fl,icode,dipl,babs)
c       if(fl.gt.10.) fl=10.
c      	invdip=INVDPC(FL,DIMO,BABS,DIPL)
c	yinv=invdip-xinv
c      	invdip_old=INVDPC_OLD(FL,DIMO,BABS,DIPL)
c       if(height.eq.heibeg) xinv=invdip


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,5,ATE1,AHH,STTE,DTE)
      if(HEIGHT.LE.HS) THEN
         TIH=TNH
      ELSE IF(HEIGHT.LT.XTETI) THEN
         TIH=BOOKER1(HEIGHT,MXSM,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
          jprint=1
          if(jf(38)) 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
c
c old drift model of Scherliess & Fejer (1999) 
c      param(1)=daynr
c      param(2)=f107d
c      call vdrift(hour,longi,param,drift)
c
        CALL vfjmodelrocstart(FJM)
        CALL vfjmodelrocinit(F107D,DAYNR,JSEA,JF107)
        CALL vfjmodelroc(FJM,HOUR,LONGI,JSEA,JF107,DRIFT)
        endif

c
c spread-F occurrence probability
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
        if(hour.lt.12.0) then
          spfhour=hour+24.0
          daynr1=daynr-1
          if(daynr1.lt.1) daynr1=idayy
          endif
        call spreadf_brazil(daynr1,month,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 in %		

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

