C  =====================================================================
C     Downloaded from ftp://nssdcftp.gsfc.nasa.gov/models/ionospheric/iri/
C     iri2007/
C     by B. Rideout on Aug 3, 2007.  Dated 6/5/2007 on ftp site.
C     A description of this code can be found
C     at http://modelweb.gsfc.nasa.gov/ionos/iri.html
C
C     Removed write statements
C     Removed unused iri_web
C     Modifed so all ap and f10.7 data comes from Madrigal, not files
C
C     $Id: irisub.f 7344 2021-03-29 18:27:39Z brideout $
C
c-----------------------------------------------------------------------        
C 
C Includes subroutines IRI_SUB and IRI_WEB to compute IRI parameters 
C for specified location, date, time, and altitude range 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 library files IRIFUN.FOR, IRITEC.FOR, IRIDREG.FOR, 
c CIRA.FOR, IGRF.FOR
c 
c !!!!! Subroutine INITIALIZE has to be called once before calling
c !!!!! iri_sub. This is already included in subroutine IRI_WEB which
c !!!!! calls iri_sub. 
C
c i/o units:  6   messages (during execution) to monitor
c            10   CCIR and URSI coefficients
c            11   alternate unit for message file (if jf(12)=.false.)
c            12   solar/ionospheric indices: ig_rz.dat       
c            13   magnetic indices: ap.dat
c            14   igrf coefficients
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 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 CHANGES FROM  IRIS13.FOR  TO   IRISUB.FOR:
c-----------------------------------------------------------------------        
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) now IGRF/POGO; dip,mlat now IGRF10; 
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
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,IAP3,F107)
C-----------------------------------------------------------------
C
C INPUT:  JF(1:30)      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         IAP3 - Array of 13 Ap3 (values integers) of the 13 3-hour Ap3
C               values leading up to the present time - added by 
C               Bill Rideout to remove dependence on file ap.dat
C         F107 - f10.7 value in 10^22 * W/m2/Hz (float) - added by 
C               Bill Rideout to remove dependence on file ap.dat
C
C    JF switches to turn off/on (.true./.false.) several options
C
C    i       .true.                  .flase.          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 - Table option      B0 - Gulyaeva (1987)                t
C    5    foF2 - CCIR            foF2 - URSI                     false
C    6    Ni - DS-78 & DY-85     Ni - DS-95 & TTS-03             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     no messages                         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                   t
C   18    IGRF dip, magbr, modip old FIELDG using POGO68/10 for 1973 t
C   19    F1 probability model   critical solar zenith angle (old)   t
C   20    standard F1            standard F1 plus L condition        t
C   21    ion drift computed     ion drift not computed          false
C   22    ion densities in %     ion densities in m-3                t
C   23    Te_tops (Aeros,ISIS)   Te_topside (Intercosmos)        false
C   24    D-region: IRI-95       Special: 3 D-region models          t
C   25    F107D from AP.DAT      F107D user input (oarr(41))         t
C   26    foF2 storm model       no storm updating                   t
C   27    IG12 from file         IG12 - user input					 t
C   28    spread-F probability 	 not computed                    false
C   29    IRI01-topside          new options as def. by JF(30)   false
C   30    IRI01-topside corr.    NeQuick topside model   	     false 
C     (29,30) = (t,t) IRIold, (f,t) IRIcor, (f,f) NeQuick, (t,f) TTS   
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) =.flase.     OARR(33)=user input for Rz12
C    jf(21) =.true.      OARR(41)=user input for daily F10.7 index
C    jf(23) =.false.     OARR(41)=user input for daily F10.7 index
C    jf(24) =.false.     OARR(41)=user input for daily F10.7 index
C          optional for jf(21:24); default is F10.7D=COV
C    jf(25) =.false.     OARR(41)=user input for daily F10.7 index
C          if oarr(41).le.0 then 12-month running mean is 
C          taken from internal file]
C    jf(27) =.flase.     OARR(39)=user input for IG12
C    jf(21) =.true.      OARR(41)=user input for daily F10.7 index
C
C
C  OUTPUT:  OUTF(1:20,1:100)
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,*)  D/E-region densities (Friedrich)
C               OUTF(14,*)  1:7 Danilov for 60,65..90km; 8:14 for a
C                  major Stratospheric Warming (SW=1) event; 15:21 
C                  for strong Winter Anomaly (WA=1) conditions
C               OUTF(15-20,*)  free
c
C            OARR(1:50)   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) = DELA
C       OARR(29) = sunrise/dec. hours OARR(30) = sunset/dec. hours
C       OARR(31) = ISEASON (1=spring) OARR(32) = NSEASON (northern)
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 (old)
C      #OARR(41) = F10.7 daily        oarr(42) = c1 (F1 shape)
C       OARR(43) = daynr
C       OARR(44) = equatorial vertical ion drift in m/s
C       OARR(45) = foF2_storm/foF2_quiet
C       OARR(46) = F1 probability without L condition 
C       OARR(47) = F1 probability with L condition incl.
C       OARR(48) = spread-F occurrence probability (Brazilian model)
C                # INPUT as well as OUTPUT parameter
C                $ special for IRIWeb (only place-holders)
c-----------------------------------------------------------------------        
C*****************************************************************
C*** THE ALTITUDE LIMITS ARE:  LOWER (DAY/NIGHT)  UPPER        ***
C***     ELECTRON DENSITY         60/80 KM       1000 KM       ***
C***     TEMPERATURES              120 KM        2500/3000 KM  ***
C***     ION DENSITIES             100 KM        1000 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,DDO,DO2,SEASON,SEADAY
      REAL       LATI,LONGI,MO2,MO,MODIP,NMF2,MAGBR,INVDIP,  
     &               NMF1,NME,NMD,MM,MLAT,MLONG,NOBO2
      CHARACTER  FILNAM*12
c-web-for webversion
c      CHARACTER FILNAM*53

      CHARACTER*128 FULL_FILENAME,MADDIR
      INTEGER IND, IAP3(13)

      DIMENSION  ARIG(3),RZAR(3),F(3),RIF(4),E(4),XDELS(4),DNDS(4),
     &  FF0(988),XM0(441),F2(13,76,2),FM3(9,49,2),ddens(3,7),elg(7),
     &  FF0N(988),XM0N(441),F2N(13,76,2),FM3N(9,49,2),INDAP(13),
     &  AMP(4),HXL(4),SCL(4),XSM(4),MM(5),DTI(4),AHH(7),
     &  STTE(6),DTE(5),ATE(7),TEA(6),XNAR(2),param(2),OARR(50),
     &  OUTF(20,100),PG1O(80),PG2O(32),PG3O(80),PF1O(12),PF2O(4),
     &  PF3O(12),HO(4),MO(5),DDO(4),HO2(2),MO2(3),DO2(2),DION(7),
     &  osfbr(25)

      LOGICAL    EXT,SCHALT,TECON(2),sam_mon,sam_yea,sam_ut,
     &  sam_date,F1REG,FOF2IN,HMF2IN,URSIF2,LAYVER,DY,GULB0,DREG,
     &  rzino,FOF1IN,HMF1IN,FOEIN,HMEIN,RZIN,sam_doy,F1_OCPRO,
     &  F1_L_COND,NODEN,NOTEM,NOION,TENEOP,OLD79,JF(30),URSIFO,
     &  igin,igino,dnight,enight,fnight,TOPO,TOPC

      COMMON /CONST/UMR  /const1/humr,dumr   /ARGEXP/ARGMAX
     &   /BLOCK1/HMF2,NMF2,HMF1,F1REG  /BLOCK2/B0,B1,C1  
     &   /BLOCK3/HZ,T,HST              /BLOCK4/HME,NME,HEF 
     &   /BLOCK5/ENIGHT,E              /BLOCK6/HMD,NMD,HDX
     &   /BLOCK7/D1,XKK,FP30,FP3U,FP1,FP2
     &   /BLOCK8/HS,TNHS,XSM,MM,DTI,MXSM         
     &   /BLOTN/XSM1,TEXOS,TLBDH,SIGMA /BLOTE/AHH,ATE1,STTE,DTE
     &   /BLO10/BETA,ETA,DELTA,ZETA       
     &   /BLO11/B2TOP,TC3,itopn,alg10,hcor1       
     &   /iounit/konsol,monito

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

        save
	
      MADDIR = 'MADROOT'
      IND = index(MADDIR, ' ') - 1        

        DO KI=1,20
        do kk=1,100
        OUTF(KI,kk)=-1.
        end do
        end do

        do 7398 KI=1,13
        INDAP(KI) = IAP3(KI)
7398    continue

        do 8398 kind=7,14,1
        oarr(kind)=-1.
8398    continue

        do 8378 kind=17,32,1
        oarr(kind)=-1.
8378    continue

        do 8478 kind=34,38,1
        oarr(kind)=-1.
8478    continue

        oarr(40)=-1.
        do 8428 kind=42,50,1
        oarr(kind)=-1.
8428    continue

C
C PROGRAM CONSTANTS
C
        icalls=icalls+1
        ARGMAX=88.0
        pi=ATAN(1.0)*4.
        UMR=pi/180.
        humr=pi/12.
        dumr=pi/182.5
        ALOG2=ALOG(2.)
        ALG10=ALOG(10.)
        ALG100=ALOG(100.)
        numhei=int(abs(heiend-heibeg)/abs(heistp))+1
        if(numhei.gt.100) numhei=100
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
        DDO(1)=9
        DDO(2)=5
        DDO(3)=5
        DDO(4)=25
        DO2(1)=5
        DO2(2)=5
        XNAR(1)=0.0
        XNAR(2)=0.0
        DTE(1)=5.
        DTE(2)=5.
        DTE(3)=10.
        DTE(4)=20.
        DTE(5)=20.
        DTI(1)=10.
        DTI(2)=10.
        DTI(3)=20.
        DTI(4)=20.
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 set konsol=1 if you do not want the konsol information
c
        MONITO=6
        IUCCIR=10
c-web- special for web version
c-web-      KONSOL=1
        KONSOL=6
        if(.not.jf(12)) konsol=1
c        if(.not.jf(12)) then
c                konsol=11
c                open(11,file='messages.txt')
c                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.
      DY=(.not.jf(6))
      LAYVER=(.not.jf(11))
      OLD79=(.not.jf(7))
      GULB0=(.not.jf(4))
      F1_OCPRO=(jf(19))
      f1_l_cond=.false.
      if(F1_OCPRO) F1_L_COND=(.not.jf(20))
      DREG=jf(24)
      TOPO=jf(29)
      TOPC=jf(30)
c
c rz12, F10.7D 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)=-1.
      ENDIF

      IF(.not.jf(25)) THEN
          f107d=OARR(41)
      else
          oarr(41)=-1.
      ENDIF
c
c Topside density ....................................................
c
        if(TOPO) then
             if (TOPC) then 
                 itopn=0
             else
                 itopn=3
             endif
        else 
             if (TOPC) then
                 itopn=1
             else
                 itopn=2
             endif
        endif
c
c F2 peak density ....................................................
c
      FOF2IN=(.not.jf(8))
       IF(FOF2IN) THEN
          AFOF2=OARR(1)
          IF(AFOF2.GT.100.) AFOF2=SQRT(AFOF2/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
          AFOF1=OARR(3)
          IF(AFOF1.GT.100.) AFOF1=SQRT(AFOF1/1.24E10)
        else
          oarr(3)=-1.
          ENDIF
c
c F1 peak altitude ..................................................
c
      HMF1IN=(.not.jf(14))
       IF(HMF1IN) then
                AHMF1=OARR(4)
C                if(.not.layver.and.(konsol.gt.1)) write(konsol,1939)
C1939  format(' *Ne* User input of hmF1 is only possible for the LAY-',
C     &          'version')
        else
                oarr(4)=-1.
        endif
c
c E peak density ....................................................
c
      FOEIN=(.not.jf(15))
       IF(FOEIN) THEN
          AFOE=OARR(5)
          IF(AFOE.GT.100.) AFOE=SQRT(AFOE/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 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       continue
           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.gt.1.or.konsol.eq.1) goto 8201
C          write(konsol,2911) 
        if(NODEN) goto 2889
C          if(LAYVER) write(konsol,9012) 
C          if(GULB0) write(konsol,9013) 
C          if(OLD79) write(konsol,9014) 
C          if(TOPO.and.(.not.TOPC)) write(konsol,9206)
C          if(.not.TOPO) then
C                if(TOPC) then
C                      write(konsol,9204) 
C                else
C                      write(konsol,9205) 
C                endif
C                endif
          if(FOF2IN) then
C                write(konsol,9015) 
                goto 2889
                endif
C        if(URSIF2) then
C                write(konsol,9016) 
C        else
C                write(konsol,9017) 
C        endif
C        if(HMF2IN) write(konsol,9018) 
C        if(fof1in) write(konsol,9019) 
C        if(HMF1IN.and.LAYVER) write(konsol,9021) 
C        if(foein) write(konsol,9022) 
C        if(HMEIN) write(konsol,9023) 
C        if(F1_OCPRO) write(konsol,9024) 
C        if(F1_L_COND) write(konsol,9025) 
C        if(DREG) then 
C            write(konsol,9026) 
C        else
C            write(konsol,9027) 
C        endif
        if(jf(26)) then
            if(fof2in) then 
C                  write(konsol,9028) 
                  jf(26)=.false.
C            else
C                  write(konsol,9029) 
            endif
            endif

2889    continue

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

        if(NOTEM) goto 8201
C          if(TENEOP) write(konsol,9032) 
C          if(jf(23)) then 
C            write(konsol,9033) 
C          else
C            write(konsol,9034) 
C          endif
2911    format('*** IRI parameters are being calculated ***')
9012    format('Ne, E-F: The LAY-Version is prelimenary.',
     &          ' Erroneous profile features can occur.')
9013    format('Ne, B0: Bottomside thickness is ',
     &          'obtained with Gulyaeva-1987 model.')
9014    format('Ne: No upper limit for F10.7 in',
     &          ' topside formula.')
9204    format('Ne: Corrected Topside Formula')
9205    format('Ne: NeQuick for Topside')
9206    format('Ne: TTS for Topside')
9015    format('Ne, foF2/NmF2: provided by user.')
9016    format('Ne, foF2: URSI model is used.')
9017    format('Ne, foF2: CCIR model is used.')
9018    format('Ne, hmF2/M3000F2: provided by user.')
9019    format('Ne, foF1/NmF1: provided by user.')
9021    format('Ne, hmF1: provided by user.')
9022    format('Ne, foE/NmE: provided by user.')
9023    format('Ne, hmE: provided by user.')
9024    format('Ne, foF1: probability function used.')
9025    format('Ne, foF1: L condition cases included.')
9026    format('Ne, D: IRI-90')
9027    format('Ne, D: IRI-90, DRS-95,and FIRI-01)')
9028    format('Ne, foF2: Storm model turned off if foF2 or',
     &          ' NmF2 user input')
9029    format('Ne, foF2: storm model included')
9039    format('Ion Com.: DS-78 & DY-85')
9031    format('Ion Com.: DS-95 & TTS-03')
9032    format('Te: Temperature-density correlation is used.')
9033    format('Te: Aeros/ISIS model')
9034    format('Te: Interkosmos model')

8201    continue

C
C CALCULATION OF DAY OF YEAR OR MONTH/DAY AND DECIMAL 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
        if(iyear.lt.100) iyear=iyear+1900
        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/idayy
        
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)

C        CALL GGM(JMAG,XLONGI1,XLATI1,XMLONG1,XMLAT1)

        call igrf_dip(lati,longi,ryear,300.0,dip,magbr,modip)

        if(.not.jf(18)) then
        CALL FIELDG(LATI,LONGI,300.0,XMA,YMA,ZMA,BET,DIP,DEC,MODIP)
        	MAGBR=ATAN(0.5*TAN(DIP*UMR))/UMR
        	endif

        ABSLAT=ABS(LATI)
        ABSMLT=ABS(MLAT)
        ABSMDP=ABS(MODIP)
        ABSMBR=ABS(MAGBR)

c
C CALCULATION OF UT/LT  ...............
c
        IF(DHOUR.le.24.0) then
        	HOUR=DHOUR					! dhour =< 24 is LT
            ut=hour-longi/15.
            if(ut.lt.0) ut=ut+24.
		else
        	UT=DHOUR-25.				 ! dhour>24 is UT+25
        	hour=ut+longi/15.	 		 ! hour becomes LT
        	if(hour.gt.24.) hour=hour-24.
        endif

c
c season assumes equal length seasons (92 days) with spring 
c (season=1) starting at day-of-year=47; for lati<0 adjustment 
c for southern hemisphere is made. some models require the
c seasonal month or the seasonal day-of year
c 
c zmonth is decimal month (Jan 1 = 1.0 and Dec 31 = 12.97)
c		zmonth = month + (iday*1.-1.)/nrdaym
c is not used currently
      
      SEASON=INT((DAYNR+45.0)/92.0)
      IF(SEASON.LT.1) SEASON=4
      NSEASN=SEASON				! Northern hemisphere season
      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
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_ut=(ut.eq.ut0)
        if(sam_date.and..not.rzino.and..not.rzin.
     &                   and..not.igin.and..not.igino) goto 2910
     
        call tcon(iyear,month,iday,daynr,rzar,arig,ttt,nmonth)
        if(nmonth.lt.0) goto 3330		! jump to end of program

        if(RZIN) then
        	rrr = arzin
        	rzar(1) = rrr
        	rzar(2) = rrr
        	rzar(3) = rrr
c       	zi=-12.349154+(1.4683266-2.67690893e-03*rrr)*rrr
c       	if(zi.gt.174.0) zi=174.0
c       	arig(1) = zi
c       	arig(2) = zi
c       	arig(3) = zi
        	endif
        if(IGIN) then
        	zi = aigin
        	arig(1) = zi
        	arig(2) = zi
        	arig(3) = zi
        	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)
        COVSAT=cov
        if(covsat.gt.188.) covsat=188
        
        if(jf(25)) then
        	f107d=cov
C               Modified by Bill Rideout to remove dependence on ap.dat file
C        	call APF_ONLY(iyear,month,iday,F107DX)
                if (F107 > 0.0) then
		    F107DX = F107
		else
		    F107DX = -3.0
		endif
            if(f107dx.gt.-2) f107d=f107dx
			endif

C
C CALCULATION OF SOLAR ZENITH ANGLE (XHI/DEG), SUN DECLINATION ANGLE 
C (SUNDEC),SOLAR ZENITH ANGLE AT NOON (XHINON) AND TIME OF LOCAL 
C SUNRISE/SUNSET (SAX, SUX; dec. hours) AT 70 KM (D-REGION), 110 KM
C (E-REGION), 200 KM (F1-REGION), AND 500 KM (TE, TI).
C

2910    continue
        CALL SOCO(daynr,HOUR,LATI,LONGI,70.,SUNDEC,XHI1,SAX70,SUX70)
        CALL SOCO(daynr,HOUR,LATI,LONGI,110.,SUD1,XHI2,SAX110,SUX110)
        CALL SOCO(daynr,HOUR,LATI,LONGI,200.,SUD1,XHI,SAX200,SUX200)
        CALL SOCO(daynr,HOUR,LATI,LONGI,500.,SUD1,XHI3,SAX500,SUX500)
        CALL SOCO(daynr,12.0,LATI,LONGI,110.,SUNDE1,XHINON,SAX1,SUX1)
        DNIGHT=.FALSE.
        if(abs(sax70).gt.25.0) then
                if(sax70.lt.0.0) DNIGHT=.TRUE.
                goto 9334
                endif
        if(SAX70.le.SUX70) goto 1386
        if((hour.gt.sux70).and.(hour.lt.sax70)) dnight=.true.
        goto 9334
1386    IF((HOUR.GT.SUX70).OR.(HOUR.LT.SAX70)) 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
      HNEA=65.
      IF(DNIGHT) HNEA=80.
      HNEE=2000.
      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 (foF2), density (NmE), and height (hmE)
c
        IF(FOEIN) THEN
          FOE=AFOE
        ELSE
          FOE=FOEEDI(COV,XHI,XHINON,ABSLAT)
        ENDIF
        NME=1.24E10*FOE*FOE

        IF(HMEIN) THEN
          HME=AHME
        ELSE
          HME=110.0
        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(URSIF2.NEQV.URSIFO) GOTO 7797
      if(.not.rzin.and..not.rzino.and..not.igin.and..not.igino) then
          IF(sam_mon.AND.(nmonth.EQ.nmono).and.sam_yea) GOTO 4292
          IF(sam_mon) GOTO 4293
          endif
c
c the program expects the coefficients files in ASCII format; if you
C want to use the binary version of the coefficients, please use the
C the statements that are commented-out below and comment-out the
C ASCII-related statements.
c
7797    URSIFO=URSIF2
        WRITE(FILNAM,104) MONTH+10
104         FORMAT('ccir',I2,'.asc')
c-binary- if binary files than use:
c-binary-104   FORMAT('ccir',I2,'.bin')
c-web- special for web-version:
c104   FORMAT('/usr/local/etc/httpd/cgi-bin/models/IRI/ccir',I2,'.asc')
c-web- special for web-version:
c
      FULL_FILENAME = MADDIR(1:IND) // '/bin/' // FILNAM
        OPEN(IUCCIR,FILE=FULL_FILENAME,STATUS='OLD',ERR=8448,
     &          FORM='FORMATTED')
c-binary- if binary files than use:
c-binary-     &          FORM='UNFORMATTED')
c
        READ(IUCCIR,4689) F2,FM3
4689    FORMAT(1X,4E15.8)
c-binary- if binary files than use:
c-binary-        READ(IUCCIR) F2,FM3

c
        CLOSE(IUCCIR)
C
C then URSI if chosen ....................................
C
        if(URSIF2) then
          WRITE(FILNAM,1144) MONTH+10
1144          FORMAT('ursi',I2,'.asc')
c-web- special for web-version:
c1144  FORMAT('/usr/local/etc/httpd/cgi-bin/models/IRI/ursi',I2,'.asc')
c-binary- if binary files than use:
c-binary-1144          FORMAT('ursi',I2,'.bin')

      FULL_FILENAME = MADDIR(1:IND) // '/bin/' // FILNAM
          OPEN(IUCCIR,FILE=FULL_FILENAME,STATUS='OLD',ERR=8448,
     &         FORM='FORMATTED')
c-binary- if binary files than use:
c-binary-     &         FORM='UNFORMATTED')

          READ(IUCCIR,4689) F2
c-binary- if binary files than use:
c-binary-          READ(IUCCIR) 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
	FULL_FILENAME = MADDIR(1:IND) // '/bin/' // FILNAM
        OPEN(IUCCIR,FILE=FULL_FILENAME,STATUS='OLD',ERR=8448,
     &          FORM='FORMATTED')
c-binary- if binary files than use:
c-binary-     &          FORM='unFORMATTED')

        READ(IUCCIR,4689) F2N,FM3N
c-binary- if binary files than use:
c-binary-        READ(IUCCIR) F2N,FM3N

        CLOSE(IUCCIR)

C
C then URSI if chosen .....................................
C
        if(URSIF2) then
          WRITE(FILNAM,1144) NMONTH+10
	  FULL_FILENAME = MADDIR(1:IND) // '/bin/' // FILNAM
          OPEN(IUCCIR,FILE=FULL_FILENAME,STATUS='OLD',ERR=8448,
     &         FORM='FORMATTED')
c-binary- if binary files than use:
c-binary-     &         FORM='unFORMATTED')

          READ(IUCCIR,4689) F2N
c-binary- if binary files than use:
c-binary-          READ(IUCCIR) F2N

          CLOSE(IUCCIR)
          endif

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

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

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

4292    zfof2  =  FOUT(MODIP,LATI,LONGI,UT,FF0)
        fof2n  =  FOUT(MODIP,LATI,LONGI,UT,FF0N)
        zm3000 = XMOUT(MODIP,LATI,LONGI,UT,XM0)
        xm300n = XMOUT(MODIP,LATI,LONGI,UT,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

501     IF(FOF2IN) THEN
          FOF2=AFOF2
        ELSE
          FOF2=YFOF2
        ENDIF
        NMF2=1.24E10*FOF2*FOF2

        IF(HMF2IN) THEN
          IF(AHMF2.LT.50.0) THEN
          	XM3000=AHMF2
          	HMF2=HMF2ED(MAGBR,RSSN,FOF2/FOE,XM3000)
          ELSE
          	HMF2=AHMF2
c          	XM3000=XM3000HM(MAGBR,RSSN,FOF2/FOE,HMF2)
          ENDIF
        ELSE
          HMF2=HMF2ED(MAGBR,RSSN,FOF2/FOE,XM3000)
        ENDIF

c
c stormtime updating
c

      stormcorr=-1.
      if(jf(26)) then
C            removed by B. Rideout
C           if(.not.sam_date.or..not.sam_ut)  
C     &          call apf(iyear,month,iday,ut,indap)
            if(indap(1).gt.-2) then
                icoord=1
                kut=int(ut)
                call STORM(indap,lati,longi,icoord,cglat,kut,
     &                     daynr,stormcorr)
                fof2=fof2*stormcorr
                NMF2=1.24E10*FOF2*FOF2
                endif
            endif
        nmono=nmonth
        MONTHO=MONTH
        iyearo=iyear
        idaynro=daynr
        rzino=rzin
        igino=igin
        ut0=ut

c
c topside profile parameters .............................
c
        COS2=COS(MLAT*UMR)
        COS2=COS2*COS2
        FLU=(COVSAT-40.0)/30.0
c option to use unlimiited F10.7M for the topside
        IF(OLD79) FLU=(COV-40.0)/30.0
c previously: IF(OLD79) ETA1=-0.0070305*COS2
        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) + 
     &     FOF2* (0.0024287 + 0.0042810 * COS2  - 0.0001528 * FOF2)
     &     

c        fluu=flu
corr    if(fluu.gt.flumax) fluu=flumax

        ZETA = 0.078922 - 0.0046702 * COS2 -  
     &    FLU * (0.019132  - 0.0076545 * COS2) +
     &    FOF2* (0.0032513 + 0.0060290 * COS2  - 0.00020872 * FOF2)

        BETA=-128.03 + 20.253 * COS2 -
     &    FLU * (8.0755  + 0.65896 * COS2) +
     &    FOF2* (0.44041 + 0.71458 * COS2 - 0.042966 * FOF2)

        Z=EXP(94.5/BETA)
corr    Z=EXP(94.45/BETA)
        Z1=Z+1
        Z2=Z/(BETA*Z1*Z1)
        DELTA=(ETA/Z1-ZETA/2.0)/(ETA*Z2+ZETA/400.0)
c
c Correction term for topside
c
      if(itopn.eq.1) then
          zmp1 = exp(modip / 10.)
          zmp11 = 1. + zmp1
          zmp111 = zmp1 / (zmp11 * zmp11)
          zmp2 = exp(modip / 19.)
          zmp22 = 1. + zmp2
          zmp222 = zmp2 / (zmp22 * zmp22) 
          r2n = -0.84 - 1.6 * zmp111  
          r2d = -0.84 - 0.64 * zmp111 
          x1n = 230. - 700. * zmp222  
          x1d = 550. - 1900. * zmp222
          r2 = HPOL(HOUR,r2d,r2n,SAX500,SUX500,1.,1.)
          x1 = HPOL(HOUR,x1d,x1n,SAX500,SUX500,1.,1.)
          hcor1 = hmF2 + x1
          x12 = 1500. - x1
	      tc3 = r2 / x12
          endif
c
c NeQuick topside parameters (use CCIR-M3000F2 even if user-hmF2)
c
      if (itopn.eq.2) then
         dNdHmx=-3.467+1.714*log(foF2)+2.02*log(xM3000)
         dNdHmx=exp(dNdHmx)*0.01
         B2bot=0.04774*FOF2*FOF2/dNdHmx
         b2k = 3.22-0.0538*foF2-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 F layer - bottomside thickness parameter B0 and shape parameters B1
c
                           
	    B1=HPOL(HOUR,1.9,2.6,SAX200,SUX200,1.,1.)
        if(GULB0) then
                call ROGUL(SEADAY,XHI,SEAX,GRAT)
                if(FNIGHT) GRAT=0.91-HMF2/4000.
                BCOEF=B1*(B1*(0.0046*B1-0.0548)+0.2546)+0.3606
                B0CNEW=HMF2*(1.-GRAT)
                B0=B0CNEW/BCOEF
        else
                B0 = B0_98(HOUR,SAX200,SUX200,NSEASN,RSSN,LONGI,MODIP)
        endif
c
c F1 layer height hmF1, critical frequency foF1, peak density NmF1
c
        IF(FOF1IN) THEN
            FOF1=AFOF1
        ELSE
            FOF1=FOF1ED(ABSMBR,RSSN,XHI)
        ENDIF
c F1 layer thickness parameter c1
        c1 = f1_c1(modip,hour,sux200,sax200)
c F1 occurrence probability with Scotto et al. 1997
        call f1_prob(xhi,mlat,rssn,f1pbw,f1pbl)
c F1 occurrence probability Ducharme et al.        
        f1pbo = 0.0 
        if(fof1in.or.((.not.fnight).and.(fof1.gt.0.0))) f1pbo=1. 

        if(f1_ocpro) then
            f1pb = f1pbw
            if(f1_l_cond) f1pb = f1pbl
            f1reg=.false.
            if((fof1in).or.(f1pb.ge.0.5)) f1reg=.true.
        else
            f1pb = f1pbo
            f1reg=.false.
            if(f1pb.gt.0.0) f1reg=.true.
        endif
            
        NMF1=1.24E10*FOF1*FOF1

c
c E-valley: depth, width, height of deepest point (HDEEP),
c 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
C        if(konsol.gt.1) WRITE(KONSOL,650)
C650     FORMAT(1X,'*NE* E-REGION VALLEY CAN NOT BE MODELLED')
600   WIDTH=.0
667   HEF=HME+WIDTH
      hefold=hef
      VNER = (1. - ABS(DEPTH) / 100.) * NME

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

2727  continue
      hmex=hme-9.
      NMD=XMDED(XHI,RSSN,4.0E8)
      HMD=HPOL(HOUR,81.0,88.0,SAX70,SUX70,1.,1.)
      F(1)=HPOL(HOUR,0.02+0.03/DELA,0.05,SAX70,SUX70,1.,1.)
      F(2)=HPOL(HOUR,4.6,4.5,SAX70,SUX70,1.,1.)
      F(3)=HPOL(HOUR,-11.5,-4.0,SAX70,SUX70,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/NME))
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/nme)/(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(xhi,month,f107d,vKp,f5SW,f6WA,elg)
          do ii=1,7
            ddens(1,ii)=10**(elg(ii)+6)
            enddo
          f5sw=1
          f6wa=0
          call DRegion(xhi,month,f107d,vKp,f5SW,f6WA,elg)
          do ii=1,7
            ddens(2,ii)=10**(elg(ii)+6)
            enddo
          f5sw=0
          f6wa=1
          call DRegion(xhi,month,f107d,vKp,f5SW,f6WA,elg)
          do ii=1,7
            ddens(3,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 omit F1 feature if nmf1*0.9 is smaller than nme
       bnmf1=0.9*nmf1
       if(nme.ge.bnmf1) goto 9427

9245   XE2H=XE2(HEF)
       if(xe2h.gt.bnmf1) then
            hef=hef-1
            if(hef.le.hme) goto 9427
            goto 9245
            endif      
        CALL REGFA1(HEF,HMF2,XE2H,NMF2,0.001,NMF1,XE2,SCHALT,HMF1)
        IF(.not.SCHALT) GOTO 3801

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

C9427    if(konsol.gt.1) WRITE(KONSOL,11) 
C11      FORMAT(1X,'*NE* HMF1 IS NOT EVALUATED BY THE FUNCTION XE2'/
C     &        1X,'CORR.: NO F1 REGION, B1=3, C1=0.0')
9427    HMF1=0.
        F1REG=.FALSE.
c        NMF1=0.
c        C1=0.0
c
c Determine E-valley parameters if HEF was changed
c

3801     continue
         if(hef.ne.hefold) then
            width=hef-hme
            IF(ENIGHT) DEPTH=-DEPTH
            CALL TAL(HDEEP,DEPTH,WIDTH,DLNDH,EXT,E)
            IF(.NOT.EXT) GOTO 380
C            if(konsol.gt.1) WRITE(KONSOL,650)
            WIDTH=.0
            hef=hme
            hefold=hef
            goto 9245
            endif

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

380     continue

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

        hf2=hef
        xf2=xe3_1(hf2)
        if(xf2.gt.nme) goto 3885
            
        CALL REGFA1(hf1,HF2,XF1,XF2,0.001,NME,XE3_1,SCHALT,HST)
        if(schalt) goto 3885

        HZ=(HST+HF1)/2.0
        D=HZ-HST
        T=D*D/(HZ-HEF-D)
        GOTO 4933

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

C3885    if(konsol.gt.1) WRITE(KONSOL,100)
C100     FORMAT(1X,'*NE* HST IS NOT EVALUATED BY THE FUNCTION XE3')
3885    HZ=(HEF+HF1)/2.
        xnehz=xe3_1(hz)
C        if(konsol.gt.1) WRITE(KONSOL,901) HZ,HEF
C901     FORMAT(6X,'CORR.: LIN. APP. BETWEEN HZ=',F5.1,
C     &          ' AND HEF=',F5.1)
        T=(XNEHZ-NME)/(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*XHI
        ENDIF
        HHALF = GRAT * HMF2
        HV1R = HME + WIDTH
        HV2R = HME + HDEEP
        HHMF2 = HMF2
        CALL INILAY(FNIGHT,F1REG,NMF2,NMF1,NME,VNER,HHMF2,HMF1M,HME,
     &                  HV1R,HV2R,HHALF,HXL,SCL,AMP,IIQU)
C        IF((IIQU.EQ.1).and.(konsol.gt.1)) WRITE(KONSOL,7733)
C7733   FORMAT('*NE* LAY amplitudes found with 2nd choice of HXL(1).')
C        IF((IIQU.EQ.2).and.(konsol.gt.1)) WRITE(KONSOL,7722)
C7722   FORMAT('*NE* LAY amplitudes could not be found.')

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

4933  HTA=120.0
      IF(NOTEM) GOTO 240
      SEC=UT*3600.
      CALL CIRA86(DAYNR,SEC,LATI,LONGI,HOUR,COV,TEXOS,TN120,SIGMA)
      IF(HOUR.NE.0.0) THEN
         iyz=iyear
         idz=daynr
         if(jf(18)) then
             secni=(24.-longi/15)*3600.
         else
             call ut_lt(1,utni,0.0,longi,iyz,idz)
             SECNI=utni*3600.
         endif
         CALL CIRA86(DAYNR,SECNI,LATI,LONGI,0.,COV,TEXNI,TN1NI,SIGNI)
      ELSE
         TEXNI=TEXOS
         TN1NI=TN120
         SIGNI=SIGMA
      ENDIF
      TLBDH=TEXOS-TN120
      TLBDN=TEXNI-TN1NI

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.
      TMAXN=TN(HMAXN,TEXNI,TLBDN,SIGNI)+20
      ATE(2)=HPOL(HOUR,TMAXD,TMAXN,SAX200,SUX200,1.,1.)

c Te(300km), Te(400km) from AE-C, Te(1400km), Te(3000km) from 
c ISIS, Brace and Theis

              DIPLAT=MAGBR
              CALL TEBA(DIPLAT,HOUR,NSEASN,TEA)

              icd=0              
          if(jf(23)) then

c Te at fixed heights taken from Brace and Theis

              AHH(3)=300.
              AHH(4)=400.
              AHH(5)=600.
              AHH(6)=1400.
              AHH(7)=3000.
              hte=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.)
              TEN=839.+1161./(1.+EXP(-(ABSMLT-45.)/5.))
              ATE(5)=HPOL(HOUR,TET,TEN,SAX500,SUX500,1.5,1.5)
          else

c Te at fixed heights from IK, Truhlik, Triskova, Smilauer

              AHH(3)=300.
              AHH(4)=550.
              AHH(5)=900.
              AHH(6)=1500.
              AHH(7)=2500.
              hte=2500
              dimo=0.311653
c              icd=1    ! compute INVDIP
c              isa=0    ! solar activity correction off
c              ise=0    ! season correction off
              do 3491 ijk=4,7
                 call igrf_sub(lati,longi,ryear,ahh(ijk),
     &                xl,icode,dipl,babs)
                 if(xl.gt.10.) xl=10.
                 call elteik(1,0,0,invdip,xl,dimo,
     &              babs,dipl,hour,ahh(ijk),daynr,f107d,
     &              teh2,sdte)
                 ate(ijk)=teh2
3491          continue
              ATE(3)=TEA(1)
          endif

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

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

c Te corrected and Te > Tn enforced

      TNAHH2=TN(AHH(2),TEXOS,TLBDH,SIGMA)
      IF(ATE(2).LT.TNAHH2) ATE(2)=TNAHH2
      STTE1=(ATE(2)-ATE(1))/(AHH(2)-AHH(1))
      DO 1901 I=2,6
         TNAHHI=TN(AHH(I+1),TEXOS,TLBDH,SIGMA)
         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
      STTE1=STTE2
1901  CONTINUE

c Te gradients STTE are computed for each segment

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

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

c Ti(430km) during daytime from AEROS data

      XSM1=430.0
      XSM(1)=XSM1
      Z1=EXP(-0.09*MLAT)
      Z2=Z1+1.
      TID1 = 1240.0 - 1400.0 * Z1 / ( Z2 * Z2 )
      MM(2)=HPOL(HOUR,3.0,0.0,SAX500,SUX500,1.,1.)

c Ti(430km) duirng 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,SAX500,SUX500,1.,1.)

c Tn < Ti < Te enforced

      TEN1=ELTE(XSM1)
      TNN1=TN(XSM1,TEXNI,TLBDN,SIGNI)
      IF(TEN1.LT.TNN1) TEN1=TNN1
      IF(TI1.GT.TEN1) TI1=TEN1
      IF(TI1.LT.TNN1) TI1=TNN1

c Tangent on Tn profile determines HS

        TI13=TEDER(130.)
        TI50=TEDER(500.)
      CALL REGFA1(130.0,500.0,TI13,TI50,0.01,TI1,TEDER,SCHALT,HS)
        IF(SCHALT) HS=200.
        TNHS=TN(HS,TEXOS,TLBDH,SIGMA)
        MM(1)=DTNDH(HS,TEXOS,TLBDH,SIGMA)
        IF(SCHALT) MM(1)=(TI1-TNHS)/(XSM1-HS)
        MXSM=2

c XTETI is altitude where Te=Ti

2391    XTTS=500.
        X=500.
2390    X=X+XTTS
        IF(X.GE.AHH(7)) GOTO 240
        TEX=ELTE(X)
        TIX=TI(X)
        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 

        MXSM=3
        MM(3)=STTE(6)
        XSM(2)=XTETI
        IF(XTETI.GT.AHH(6)) GOTO 240
        MXSM=4
        MM(3)=STTE(5)
        MM(4)=STTE(6)
        XSM(3)=AHH(6)
        IF(XTETI.GT.AHH(5)) GOTO 240
        MXSM=5
        DTI(1)=5.
        DTI(2)=5.
        MM(3)=STTE(4)
        MM(4)=STTE(5)
        MM(5)=STTE(6)
        XSM(3)=AHH(5)
        XSM(4)=AHH(6)

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

240   IF(NOION) GOTO 141
      HNIA=100.
      if(DY) hnia=75
      HNIE=2000.
      if(dy) goto 141
C
C INPUT OF THE ION DENSITY PARAMETER ARRAYS PF1O,PF2O AND PF3O......
C
      RIF(1)=2.
      IF(ABSLAT.LT.30.0) RIF(1)=1.
      RIF(2)=2.
      IF(COV.LT.100.0) RIF(2)=1.
      RIF(3)=SEASON
      IF(SEASON.EQ.1) RIF(3)=3.
      RIF(4)=1.
      IF(FNIGHT) RIF(4)=2.
      CALL KOEFP1(PG1O)
      CALL KOEFP2(PG2O)
      CALL KOEFP3(PG3O)
      CALL SUFE(PG1O,RIF,12,PF1O)
      CALL SUFE(PG2O,RIF, 4,PF2O)
      CALL SUFE(PG3O,RIF,12,PF3O)

c
c calculate O+ profile parameters
c

      IF(FNIGHT) THEN
        ZZZ1=0.0
      ELSE
        ZZZ1=COS(XHI*UMR)
      ENDIF
      msumo=4
      RDOMAX=100.0
      MO(1)=EPSTEP(PF1O(1),PF1O(2),PF1O(3),PF1O(4),ZZZ1)
      MO(2)=EPSTEP(PF1O(5),PF1O(6),PF1O(7),PF1O(8),ZZZ1)
      MO(3)=0.0
      HO(1)=EPSTEP(PF1O(9),PF1O(10),PF1O(11),PF1O(12),ZZZ1)
      HO(2)=290.0
      IF((RIF(2).EQ.2.).AND.(RIF(3).EQ.2.)) HO(2)=237.0
      HO(4)=PF2O(1)
      ho05=pf2o(4)
      MO(4)=PF2O(2)
      MO(5)=PF2O(3)

c
c adjust gradient MO(4) of O+ profile segment above F peak
c

7100    HO(3)=(ALG100-MO(5)*(HO(4)-ho05))/MO(4)+HO(4)
        IF(HO(3).LE.HO(2)+20.) THEN
                MO(4)=MO(4)-0.001
                GOTO 7100
                endif
        hfixo=(ho(2)+ho(3))/2.

c
c find height H0O of maximum O+ relative density
c

      DELX=5.0
      X=HO(2)
      YMAXX=0.0
7102  X=X+DELX
      Y=RPID(X,HFIXO,RDOMAX,msumo,MO,DDO,HO)
      IF(Y.LE.YMAXX) then
        if(delx.le.0.1) GOTO 7104
        x=x-delx
        delx=delx/5.
      ELSE
        YMAXX=Y
      ENDIF
      GOTO 7102 
7104  H0O=X-DELX/2.
7101  if(y.lt.100.0) goto 7103
      rdomax=rdomax-0.01
      y=rpid(h0o,hfixo,rdomax,msumo,mo,ddo,ho)
      goto 7101
7103  yo2h0o=100.-y
      yoh0o=y

c
c calculate parameters for O2+ profile
c

      hfixo2 = pf3o(1)
      rdo2mx = pf3o(2) 
      DO 7105 L=1,2
          I = L * 2
          HO2(L)=PF3O(1+I)+PF3O(2+I)*ZZZ1
          MO2(L+1)=PF3O(7+I)+PF3O(8+I)*ZZZ1
7105  CONTINUE

        MO2(1)=PF3O(7)+PF3O(8)*ZZZ1
        if(hfixo2.gt.ho2(1)) then
           ymo2z=mo2(2)
        else
           ymo2z=mo2(1)
        endif
        aldo21=alog(rdo2mx)+ymo2z*(ho2(1)-hfixo2)
        hfixo2=(ho2(2)+ho2(1))/2.
        rdo2mx=exp(aldo21+mo2(2)*(hfixo2-ho2(1)))

c
c make sure that rd(O2+) is less or equal 100-rd(O+) at O+ maximum
c

7106  Y=RPID(H0O,hfixo2,rdo2mx,2,MO2,DO2,HO2)
      IF(Y.GT.yo2h0o) then
        MO2(3)=MO2(3)-0.02
        GOTO 7106
        endif

c
C use ratio of NO+ to O2+ density at O+ maximum to calculate
c NO+ density above the O+ maximum (H0O)
c

      IF(y.LT.1.) then
        NOBO2=0.0
      ELSE
        NOBO2= (yo2h0o-y)/y
      ENDIF

C
C CALCULATION FOR THE REQUIRED HEIGHT RANGE.......................
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

300   IF(NODEN) GOTO 330
      IF((HEIGHT.GT.HNEE).OR.(HEIGHT.LT.HNEA)) GOTO 330
      IF(LAYVER) THEN
          ELEDE=-9.
          IF(IIQU.LT.2) ELEDE=
     &            XEN(HEIGHT,HMF2,NMF2,HME,4,HXL,SCL,AMP)
          outf(1,kk)=elede
          goto 330
          endif
      ELEDE=XE_1(HEIGHT)
c
c TTS model option for topside
c
      if((itopn.eq.3).and.(height.ge.400.0)) then
              call igrf_sub(lati,longi,ryear,height,
     &                xl1,icode,dipl1,babs1)
              if(xl1.gt.10.) xl1=10.
c              icd=1    ! compute INVDIP
        	  dimo=0.311653
      	      call CALNE(1,INVDIP,xl1,DIMO,babs1,
     &                  DIPL,hour,height,daynr,F107D,ELEDE)
           	  endif

c
c electron density in m-3 in outf(1,*)
c

      OUTF(1,kk)=ELEDE

c
c plasma temperatures
c

330   IF(NOTEM) GOTO 7108
      IF((HEIGHT.GT.HTE).OR.(HEIGHT.LT.HTA)) GOTO 7108
        TNH=TN(HEIGHT,TEXOS,TLBDH,SIGMA)
        TIH=TNH
        IF(HEIGHT.GE.HS) TIH=TI(HEIGHT)
        TEH=ELTE(HEIGHT)
        IF(TIH.LT.TNH) TIH=TNH
        IF(TEH.LT.TIH) TEH=TIH
        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
      if(DY) then
            call ioncomp(ryear,daynr,iseamon,hour,height,xhi,
     &            lati,longi,cov,dion)
            ROX=DION(1)
            RHX=DION(2)
            RNX=DION(3)
            RHEX=DION(4)
            RNOX=DION(5)
            RO2X=DION(6)
            RCLUST=DION(7)
      else
            ROX=RPID(HEIGHT,HFIXO,RDOMAX,msumo,MO,DDO,HO)
            RO2X=RPID(HEIGHT,HFIXO2,rdo2mx,2,MO2,DO2,HO2)
            CALL RDHHE(HEIGHT,H0O,ROX,RO2X,NOBO2,10.,RHX,RHEX)
            RNOX=RDNO(HEIGHT,H0O,RO2X,ROX,NOBO2)
            RNX=-1.
            RCLUST=-1.
      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

c
c D region special: Friedrich&Torkar model in outf(13,*)
c

7118    if(.not.dreg) then
            outf(13,kk)=-1.
            call F00(HEIGHT,LATI,DAYNR,XHI,F107D,EDENS,IERROR)
            if(ierror.eq.0) outf(13,kk)=edens
            endif

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

c
c D region special: Danilov et al. model in outf(14,*): 
c outf(14,1:7)=Ne(h,SW=0,WA=0) h=60,65,70,75,80,85,90 km; 
C outf(14,8:14)=Ne(h,SW=1,WA=0); outf(14,15:21)=Ne(h,SW=0,WA=1)
c

      if(.not.dreg) then
            do ii=1,7  
                  outf(14,ii)=ddens(1,ii)      
                  outf(14,7+ii)=ddens(2,ii)      
                  outf(14,14+ii)=ddens(3,ii)      
                  enddo
            endif

c
c equatorial vertical ion drift
c

      drift=-1.
      if(jf(21).and.abs(magbr).lt.25.0) then
            param(1)=daynr
            param(2)=f107d
            call vdrift(hour,longi,param,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(daynr,idayy,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
C

        IF(NODEN) GOTO 6192
      OARR(1)=NMF2
      OARR(2)=HMF2
      if(f1reg) OARR(3)=NMF1
      if(f1reg) OARR(4)=XHMF1
      OARR(5)=NME
      OARR(6)=HME
      OARR(7)=NMD
      OARR(8)=HMD
      OARR(9)=HHALF
      OARR(10)=B0
      OARR(11)=VNER
      OARR(12)=HEF
6192    IF(NOTEM) GOTO 6092
      OARR(13)=ATE(2)
      OARR(14)=AHH(2)
      OARR(15)=ATE(3)
      OARR(16)=ATE(4)
      OARR(17)=ATE(5)
      OARR(18)=ATE(6)
      OARR(19)=ATE(7)
      OARR(20)=ATE(1)
      OARR(21)=TI1
      OARR(22)=XTETI
6092  OARR(23)=XHI
      OARR(24)=SUNDEC
      OARR(25)=DIP
      OARR(26)=MAGBR
      OARR(27)=MODIP
      OARR(28)=DELA
      OARR(29)=SAX200
      OARR(30)=SUX200
      OARR(31)=SEASON
      OARR(32)=NSEASN
      OARR(33)=rssn
      OARR(34)=COV
      OARR(35)=B1
      OARR(36)=xm3000
      OARR(39)=gind
      OARR(40)=f1pbo
      OARR(41)=f107d
      OARR(42)=c1
      OARR(43)=daynr
      OARR(44)=drift
      OARR(45)=stormcorr
      OARR(46)=f1pbw
      OARR(47)=f1pbl
      OARR(48)=spreadf
3330  CONTINUE
      RETURN
      END
c
c