iri90.f

C Subroutine IRI90
C
C Special-purpose version of IRI90 for use with the GLOW model.
C This is not the original IRI90.
C
C Adapted 7/93 from 10/91 version of IRIS12 by Stan Solomon.
C Replaces 3/92 adaptation of 12/90 version of IRIS12.
C Replaced height range and interval with array of heights ZKM
C Also supply number of heights NZ; unlimited number of heights allowed.
C Added DIRECT argument to specify location of CCIR and URSI files.
C Uses ASCII versions of CCIR and URSI files.
C Subroutine DFP splices directory and filename together.
C Accepts longitudes -180 to +180; maps internal longitudes 0-360.
C Added SAVE statements for volatile memory machines.
C Changed 'file not found' error to output message on unit 6 and stop.
C Changed JF(12)=.true. to write messages on unit 12 instead of 6.
C Changed call to IONCOM to use floating point variable ZMONTH.
C Note: JF(5)=.false. uses recommended (URSI) coefficients.
C Note: JF(4)=.false. uses recommended Gulyeava B0 coefficients.
C All other JF's usually should be .true.
C
C The following is the original comment from IRIS12:
C
C IRIS12.FOR ---------------------------------------- OCTOBER 1991
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------------- inlcuded on diskette ------------------------------
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
C*****************************************************************
C********* INTERNATIONAL REFERENCE IONOSPHERE (IRI). *************
C*****************************************************************
C****************    OCTOBER 1991     ****************************
C****************     SUBROUTINE      ****************************
C*****************************************************************
C
C
C-----------------------------------------------------------------
C INTERNATIONAL REFERENCE IONOSPHERE 1991
C
C INPUT:  JMAG=0/1      GEODETIC/GEOMAGNETIC LATITUDE AND LONGITUDE
C         ALATI,ALONG   LATITUDE NORTH AND LONGITUDE EAST IN DEGREES
C         RZ12 (-COV)   12-MONTHS-RUNNING MEAN OF SOLAR SUNSPOT NUMBER
C                          (OR EQUIVALENT F10.7 SOLAR RADIO FLUX AS
C                          NEGATIVE NUMBER)
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,       BEGIN, END, AND STEPWIDTH OF HEIGHT RANGE
C          HEIEND,HEISTP   IN KM (MAXIMUM NUMBER OF STEPS IS 50 !!)
C         JF(1:12)      TRUE/FALSE FLAGS FOR SEVERAL OPTIONS
C          JF(1)=.TRUE.[.FALSE.]   ELECTRON DENSITY IS [NOT] CALCULATED
C          JF(2)=T[F]    TEMPERATURES ARE [NOT] CALCULATED
C          JF(3)=T[F]    ION COMPOSITION IS [NOT] CALCULATED
C          JF(4)=T[F]    B0 FROM TABLE [FROM GULYEAVA 1987]
C          JF(5)=T[F]    F2 PEAK FROM CCIR [FROM URSI]
C          JF(6)=T[F]    ION COMP. STANDARD [DANILOV-YAICHNIKOV-1985]
C          JF(7)=T[F]    STAND. IRI TOPSIDE [IRI-79]
C          JF(8)=T[F]    NMF2 PEAK MODEL [INPUT VALUES]
C          JF(9)=T[F]    HMF2 PEAK MODEL [INPUT VALUES]
C          JF(10)=T[F]   TE MODEL [TE-NE MODEL WITH NE INPUT]
C          JF(11)=T[F]   NE STANDARD [LAY-FUNCTIONS VERSION]
C          JF(12)=T[F]   MESSAGE ARE WRITTEN TO UNIT=12 [=6]
C
C  JF(1:11)=.TRUE. GENERATES THE STANDARD IRI-90 PARAMETERS.
C  IF YOU SET JF(8)=.FALSE., THAN YOU HAVE TO PROVIDE THE F2 PEAK
C  NMF2/M-3 OR FOF2/MHZ IN OARR(1). SIMILARLY, IF YOU SET JF(9)=
C  .FALSE., THAN YOU HAVE TO PROVIDE THE F2 PEAK HEIGHT HMF2/KM IN
C  OARR(2). IF YOU SET JF(10)=.FALSE., THAN YOU HAVE TO PROVIDE THE
C  ELECTRON DENSITY IN M-3 AT 300KM AND/OR 400KM AND/OR 600KM IN
C  OARR(3), OARR(4), AND OARR(5). IF YOU WANT TO USE THIS OPTION AT
C  ONLY ONE OF THE THREE ALTITUDES, THAN SET THE DENSITIES AT THE
C  OTHER TWO TO ZERO.
C
C  OUTPUT:  OUTF(1:10,1:50)   IRI PROFILES
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,*)  PERCENTAGE OF O+ IONS IN %
C              OUTF(6,*)  PERCENTAGE OF H+ IONS IN %
C              OUTF(7,*)  PERCENTAGE OF HE+ IONS IN %
C              OUTF(8,*)  PERCENTAGE OF O2+ IONS IN %
C              OUTF(9,*)  PERCENTAGE OF NO+ IONS IN %
C                 AND, IF JF(6)=.FALSE.:
C              OUTF(10,*)  PERCENTAGE OF CLUSTER IONS IN %
C              OUTF(11,*)  PERCENTAGE OF N+ IONS IN %
C
C            OARR(1:30)   ADDITIONAL OUTPUT PARAMETERS
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) = SOLAR ZENITH ANGLE/DEG
C              OARR(24) = SUN DECLINATION/DEG
C              OARR(25) = DIP
C	       OARR(26) = DIP LATITUDE
C              OARR(27) = MODIFIED DIP LATITUDE
C              OARR(28:30) FREE
C-------------------------------------------------------------------
C*** THIS PROGRAM PRODUCES PROFILES OF                         ***
C***      ELECTRON DENSITY                                     ***
C***      NEUTRAL TEMPERATURE (CIRA 86)                        ***
C***      ELECTRON TEMPERATURE                                 ***
C***      ION TEMPERATURE                                      ***
C***      RELATIVE PERCENTAGE DENSITIES OF THE IONS            ***
C***           ATOMIC OXYGEN, HYDROGEN, HELIUM,                ***
C***           MOLECULAR OXYGEN AND NITROGEN OXYD (NO+)        ***
C*****************************************************************
C*** THE ALTITUDE LIMITS ARE:  LOWER (DAY/NIGHT)  UPPER        ***
C***     ELECTRON DENSITY         60/80 KM       1000 KM       ***
C***     TEMPERATURES              120 KM        3000 KM       ***
C***     ION DENSITIES             100 KM        1000 KM       ***
C*****************************************************************
C*     --------------------ADDRESSES------------------------     *
C*     I  PROF. K. RAWER              DR. D. BILITZA       I     *
C*     I  HERRENSTR. 43               GSFC/NSSDC CODE 633  I     *
C*     I  D-7801 MARCH                GREENBELT MD 20771   I     *
C*     I  F.R.G.                      USA                  I     *
C*     -----------------------------------------------------     *
C*****************************************************************
C*****************************************************************
C*****************************************************************
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********************  OPTIONS  **********************************
C*****************************************************************
C* FOR HMF2=0 OR FOF2=0 THE F2 PEAK VALUES ARE CALCULATED WITH   *
C* THE CCIR OR URSI MODELS. THE CCIR COEFFICIENT SET FOR THE     *
C* MONTH "mm" IS EXPECTED IN THE BINARY FILE "CCIRmm.BIN" AND    *
C* THE URSI SET IN "URSImm.BIN". IF YOU USE THE ASCII CODED      *
C* FILES "CCIRmm.ASC", YOU HAVE TO INCORPORATE THE CHANGES       *
C* INDICTED IN PROGRAM SECTION ENTITLED "READ CCIR COEFFICIENT   *
C* SET FOR CHOSEN MONTH."                                        *
C*****************************************************************
C*****************************************************************
C*****************************************************************
C
       SUBROUTINE IRI90(JF,JMAG,ALATI,ALONG,RZ12,MMDD,DHOUR,
     &                  ZKM,NZ,DIRECT,OUTF,OARR)
      use, intrinsic :: iso_fortran_env, only: stdout=>output_unit,
     &  stderr=>error_unit
      dimension zkm(nz), outf(11,nz), oarr(30)
      character(*), intent(in) :: direct
      character(1024) :: path
      character(10) filename
      INTEGER 		EGNR,AGNR,DAYNR,DDO,DO2,SEASON,SEADAY
      REAL 		LATI,LONGI,MO2,MO,MODIP,NMF2,MAGBR
      REAL  		NMF1,NME,NMD,NEI,MM,MLAT,MLONG,NOBO2
      DIMENSION  F(3),RIF(4),E(4),XDELS(4),DNDS(4)
      DIMENSION  FF0(988),XM0(441),F2(13,76,2),FM3(9,49,2)
      DIMENSION  AMP(4),HXL(4),SCL(4),B0B1(5)
      DIMENSION  CTN(3),CTNN(3),XSM(4),MM(5),DTI(4)
      DIMENSION  AHH(7),STTE(6),DTE(5),ATE(7),TEA(6),HOA(3),XNAR(3)
      DIMENSION  PG1O(80),PG2O(32),PG3O(80),PF1O(12),PF2O(4),PF3O(12)
      DIMENSION  HO(4),MO(5),DDO(4),HO2(2),MO2(3),DO2(2),DION(7)
      LOGICAL		EXT,SCHALT,NIGHT,TCON(3)
      LOGICAL		F1REG,FOF2IN,HMF2IN,URSIF2,LAYVER,DY,GULB0
      LOGICAL		NODEN,NOTEM,NOION,TENEOP
      LOGICAL           OLD79,TOPSI,BOTTO,BELOWE,JF(12),URSIFO
      COMMON	/BLOCK1/HMF2,NMF2,HMF1	         /CONST/UMR
     &		/BLOCK2/B0,B1,C1      /BLOCK3/HZ,T,HST,STR
     &  	/BLOCK4/HME,NME,HEF   /BLOCK5/NIGHT,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	 /ARGEXP/ARGMAX
      EXTERNAL 		XE1,XE2,XE3,XE4,XE5,XE6,TEDER
      DATA  HOA  /300.,400.,600./,   XNAR       /3*0.0/,
     &      XDELS   /3*5.,10./,      DNDS   /.016,.01,2*.016/,
     &      DDO	  /9,5,5,25/,        DO2        /5,5/,
     &      B0B1  /.755566,.778596,.797332,.812928,.826146/
      data icalls/0/
C
      SAVE EGNR,AGNR,DAYNR,DO2,SEASON,SEADAY,LATI,LONGI,MO2,
     &     MODIP,MAGBR,NMF1,NEI,MLAT,MLONG,NOBO2,
     &     F,RIF,FF0,XM0,F2,FM3,AMP,HXL,SCL,
     &     CTN,CTNN,ATE,TEA,HOA,XNAR,PG1O,
     &     PG2O,PG3O,PF1O,PF2O,PF3O,HO,MO,DDO,HO2,DION,EXT,
     &     SCHALT,SSIN,TCON,F1REG,FOF2IN,HMF2IN,URSIF2,LAYVER,
     &     DY,GULB0,NODEN,NOTEM,NOION,TENEOP,OLD79,TOPSI,BOTTO,BELOWE,
     &     URSIFO,MONTH,MONTHO,RG,RGO
C
C PROGAM CONSTANTS
C
	icalls=icalls+1
        HHALF = 0.

!       write(6,"('Enter iri90: icalls=',i4,' jf=',12l2,' jmax=',i4,
!    |    ' alati=',f8.2,' along=',f8.2)")
!    |    icalls,jf,jmax,alati,along
!       write(6,"('rz12=',f8.2,' mmdd=',i4,' dhour=',f8.2,' nz=',i4)")
!    |    rz12,mmdd,dhour,nz
!       write(6,"('DIRECT (file path)=',a)") direct
!       write(6,"('zkm=',/,(8f9.2))") zkm

	ARGMAX=88.0
     	UMR=ATAN(1.0)*4./180.
      	ALOG2=ALOG(2.)
	ALG100=ALOG(100.)
 	ISTART=1
        heibeg=zkm(1)
        heiend=zkm(nz)
C
C Code inserted to aleviate block data problem for PC version.
C Thus avoiding DATA statement with parameters from COMMON block.
C
        AHH(1)=120.
        AHH(2)=0.
        AHH(3)=300.
        AHH(4)=400.
        AHH(5)=600.
        AHH(6)=1400.
        AHH(7)=3000.
        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
      MONITO=6
      IUCCIR=10
      KONSOL=stderr
!      IF (JF(12)) KONSOL=12

c
c selection of density and ion composition options ..................
c

      NODEN=(.NOT.JF(1))
      NOTEM=(.NOT.JF(2))
      NOION=(.NOT.JF(3))
      DY=(.NOT.JF(6))
      LAYVER=(.NOT.JF(11))
      OLD79=(.NOT.JF(7))
      GULB0=(.NOT.JF(4))
c
c f peak density ....................................................
c
      FOF2IN=(.NOT.JF(8))
       IF(FOF2IN) THEN
          AFOF2=OARR(1)
          IF(AFOF2.GT.100.) AFOF2=SQRT(AFOF2/1.24E10)
          ENDIF
      URSIF2=(.NOT.JF(5))
c
c f peak altitude ..................................................
c
      HMF2IN=(.NOT.JF(9))
       IF(HMF2IN) AHMF2=OARR(2)
c
C TE-NE MODEL OPTION ..............................................
C
      TENEOP=(.NOT.JF(10))
        IF(TENEOP) THEN
           DO 8154 JXNAR=1,3
              XNAR(JXNAR)=OARR(JXNAR+2)
	      TCON(JXNAR)=.FALSE.
8154	      IF(XNAR(JXNAR).GT.0.) TCON(JXNAR)=.TRUE.
           ENDIF

      if(icalls.gt.1) goto 8201
	write(konsol,*) '*** IRI parameters are being calculated ***'
      if(NODEN) goto 2889
	if(LAYVER) write(konsol,*) 'Ne, E-F: The LAY-Version is ',
     &	  'prelimenary. Erroneous profile features can occur.'
	if(GULB0) write(konsol,*) 'Ne, B0: Bottomside thickness is ',
     &	  'obtained with Gulyaeva-1987 model.'
	if(OLD79) write(konsol,*) 'Ne: Using IRI-79. Correction',
     &	  ' of equatorial topside is not included.'
	if(HMF2IN) write(konsol,*) 'Ne, hmF2: Input values are used.'
	if(FOF2IN) then
	  write(konsol,*) 'Ne, foF2: Input values are used.'
	  goto 2889
	  endif
	if(URSIF2) then
	  write(konsol,*) 'Ne, foF2: URSI model is used.'
	else
	  write(konsol,*) 'Ne, foF2: CCIR model is used.'
	endif
2889  if((.not.NOION).and.(DY))
     &	   write(konsol,*) 'Ion Com.: Using Danilov-Yaichnikov-1985.'
      if((.not.NOTEM).and.(TENEOP))
     &     write(konsol,*) 'Te: Temperature-density correlation is used'
8201	continue
C
C CALCULATION OF MEAN F10.7CM SOLAR RADIO FLUX (COV)................
C CALCULATION OF RESTRICTED SOLAR ACTIVITIES (RG,COVG)..............
C
      IF(RZ12.GE.0.0) THEN
        R=RZ12
        COV=63.75+R*(0.728+R*0.00089)
      ELSE
        COV=-RZ12
        R=33.52*(SQRT(COV+85.12)-12.2)
      ENDIF
      RG=R
      COVG=COV
      IF(R.GT.150.) RG=150.
      IF(COV.GT.193.) COVG=193.
C
C CALCULATION OF GEOG. OR GEOM. COORDINATES IN DEG....................
C CALCULATION OF MAGNETIC INCLINATION (DIP), DECLINATION (DEC)........
C   DIP LATITUDE (MAGBR) AND MODIFIED DIP (MODIP). ALL IN DEGREE......
C
        IF(JMAG.GT.0) THEN
           MLAT=ALATI
           MLONG=ALONG
           if (mlong .lt. 0.) mlong=mlong+360.
        ELSE
           LATI=ALATI
           LONGI=ALONG
           if (longi .lt. 0.) longi=longi+360.
        ENDIF
        CALL GGM(JMAG,LONGI,LATI,MLONG,MLAT)
        ABSLAT=ABS(LATI)
        CALL FIELDG(LATI,LONGI,300.0,XMA,YMA,ZMA,BET,DIP,DEC,MODIP)
        MAGBR=ATAN(0.5*TAN(DIP*UMR))/UMR
	ABSMLT=ABS(MLAT)
	ABSMDP=ABS(MODIP)
	ABSMBR=ABS(MAGBR)
C
C CALCULATION OF SEASON (SUMMER=2, WINTER=4)..........................
C CALCULATION OF DAY OF YEAR AND SUN DECLINATION......................
C
  	if(MMDD.lt.0) then
		DAYNR=-MMDD
		call MODA(1,MONTH,IDAY,DAYNR)
	else
		MONTH=MMDD/100
		IDAY=MMDD-MONTH*100
		call MODA(0,MONTH,IDAY,DAYNR)
	endif
      SEASON=INT((DAYNR+45.0)/92.0)
      IF(SEASON.LT.1) SEASON=4
      NSESON=SEASON
      seaday=daynr
      IF(LATI.GT.0.0) GOTO 5592
   	SEASON=SEASON-2
    	IF(SEASON.LT.1) SEASON=SEASON+4
	seaday=daynr+183
	if(seaday.gt.366) seaday=seaday-366
C
C CALCULATION OF SOLAR ZENITH ANGLE (XHI/DEG).........................
C NOON VALUE (XHINON).................................................
C
5592  IF(DHOUR.GT.24.1) THEN
	UT=DHOUR-25.
	HOUR=UT+LONGI/15.
	IF(HOUR.GT.24.) HOUR=HOUR-24.
      ELSE
	HOUR=DHOUR
        UT=HOUR-LONGI/15.
        IF(UT.LT.0.) UT=UT+24.
      ENDIF

	CALL SOCO(DAYNR,HOUR,LATI,LONGI,SUNDEC,XHI,SAX,SUX)
	CALL SOCO(DAYNR,12.0,LATI,LONGI,SUNDE1,XHINON,SAXNON,SUXNON)

	        NIGHT=.FALSE.
	if(abs(sax).gt.25.0) then
                if(sax.lt.0.0) NIGHT=.TRUE.
		goto 1334
		endif
      	if(SAX.le.SUX) goto 1386
	if((hour.gt.sux).and.(hour.lt.sax)) night=.true.
	goto 1334
1386  	IF((HOUR.GT.SUX).OR.(HOUR.LT.SAX)) NIGHT=.TRUE.
C
C CALCULATION OF ELECTRON DENSITY PARAMETERS................
C
1334  HNEA=65.
      IF(NIGHT) 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!!!!!!! F-REGION PARAMETERS AND E-PEAK !!!!!!!!!!!!!!!!!!!!!!!!!!
      FOE=FOEEDI(COV,XHI,XHINON,ABSLAT)
      NME=1.24E10*FOE*FOE
      HME=105.0
      IF((FOF2IN).AND.(HMF2IN)) GOTO 501
      IF (URSIF2 .NEQV. URSIFO) GOTO 7797
      IF((MONTH.EQ.MONTHO).AND.(RG.EQ.RGO)) GOTO 4292
      IF(MONTH.EQ.MONTHO) GOTO 4291
C
C READ CCIR COEFFICIENT SET FOR CHOSEN MONTH....................
C
7797    WRITE(filename,104) MONTH+10
104     FORMAT('ccir',I2,'.asc')
        call dfp(direct,filename,path)
        OPEN(IUCCIR,FILE=path,STATUS='OLD',ERR=8448)
        READ(IUCCIR,4689) F2,FM3
4689    FORMAT(4E15.8)
        CLOSE(IUCCIR)
C
C READ URSI COEFFICIENT SET FOR CHOSEN MONTH....................
C
	if (URSIF2) then
	  WRITE(filename,1144) MONTH+10
1144      FORMAT('ursi',I2,'.asc')
          call dfp(direct,filename,path)
          OPEN(IUCCIR,FILE=path,STATUS='OLD',ERR=8448)
          READ(IUCCIR,4689) F2
          CLOSE(IUCCIR)
	endif
        URSIFO=URSIF2
        MONTHO=MONTH
	GOTO 4291

8448	write(stderr,*) path // ' IRI90: File not found'
        error stop
C LINEAR INTERPOLATION IN SOLAR ACTIVITY
C
4291    RR2=RG/100.
        RR1=1.-RR2
        DO 20 I=1,76
        DO 20 J=1,13
        K=J+13*(I-1)
20      FF0(K)=F2(J,I,1)*RR1+F2(J,I,2)*RR2
        DO 30 I=1,49
        DO 30 J=1,9
        K=J+9*(I-1)
30      XM0(K)=FM3(J,I,1)*RR1+FM3(J,I,2)*RR2
	RGO=RG

4292  CALL F2OUT(MODIP,LATI,LONGI,FF0,XM0,UT,YFOF2,XM3000)

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

	IF(HMF2IN) THEN
	  HMF2=AHMF2
	ELSE
          HMF2=HMF2ED(MAGBR,RG,FOF2/FOE,XM3000)
	ENDIF

        TOPSI=(HEIEND.GT.HMF2)
	BOTTO=((HEIEND.GE.HME).AND.(HEIBEG.LE.HMF2))
	BELOWE=(HEIBEG.LT.HME)
c
c topside profile parameters .............................
c
	IF(.NOT.TOPSI) GOTO 1501
      COS2=COS(MLAT*UMR)
      COS2=COS2*COS2
      FLU=(COVG-40.0)/30.0
      IF(OLD79) then
        ETA1=-0.0070305*COS2
      else
        EX=EXP(-MLAT/15.)
        EX1=EX+1
        EPIN=4.*EX/(EX1*EX1)
        ETA1=-0.02*EPIN
      endif
      ETA=0.058798+ETA1+FLU*(-0.014065+0.0069724*COS2)+
     &(0.0024287+0.0042810*COS2-0.00015280*FOF2)*FOF2
      ZETA=0.078922-0.0046702*COS2+FLU*(-0.019132+0.0076545*COS2)+
     &(0.0032513+0.0060290*COS2-0.00020872*FOF2)*FOF2
      BETA=-128.03+20.253*COS2+FLU*(-8.0755-0.65896*COS2)+(0.44041
     &+0.71458*COS2-0.042966*FOF2)*FOF2
      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 bottomside profile parameters .............................
C
1501    HMF1=HMF2
	HZ=HMF2
	HEF=HME
	IF(.not.BOTTO) GOTO 2727
	B1=3.0
C!!!!!!! INTERPOLATION FOR B0 OUT OF ARRAY B0F !!!!!!!!!!!!!!!!!!!!!
	if(GULB0) then
	  call ROGUL(SEADAY,XHI,SEAX,GRAT)
	  if(NIGHT) GRAT=0.91-HMF2/4000.
	  B0CNEW=HMF2*(1.-GRAT)
	  B0=B0CNEW/B0B1(1)
	else
          B0 = B0POL(HOUR,SAX,SUX,SEASON,RG,DELA)
	endif
C!!!!!!! F1-REGION PARAMETERS !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      F1REG=.FALSE.
      HMF1=0.
      PNMF1=0.
      C1=0.
      IF(NIGHT.OR.(SEASON.EQ.4)) GOTO 150
        FOF1=FOF1ED(ABSMBR,R,XHI)
        IF(FOF1.LT.1.E-3) GOTO 150
          F1REG=.TRUE.
          C1=.09+.11/DELA
          PNMF1=1.24E10*FOF1*FOF1
150   NMF1=PNMF1
C!!!!!!! PARAMETER FOR E AND VALLEY-REGION !!!!!!!!!!!!!!!!!!!!!
      XDEL=XDELS(SEASON)/DELA
      DNDHBR=DNDS(SEASON)/DELA
      HDEEP=HPOL(HOUR,10.5/DELA,28.,SAX,SUX,1.,1.)
      WIDTH=HPOL(HOUR,17.8/DELA,45.+22./DELA,SAX,SUX,1.,1.)
      DEPTH=HPOL(HOUR,XDEL,81.,SAX,SUX,1.,1.)
      DLNDH=HPOL(HOUR,DNDHBR,.06,SAX,SUX,1.,1.)
      IF(DEPTH.LT.1.0) GOTO 600
	IF(NIGHT) DEPTH=-DEPTH
      	CALL TAL(HDEEP,DEPTH,WIDTH,DLNDH,EXT,E)
      	IF(.NOT.EXT) GOTO 667
      	  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.) * NME
c
c Parameters below E  .............................
c
2727	IF(.not.BELOWE) GOTO 2726
C!!!!!!!D-REGION PARAMETER!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      NMD=XMDED(XHI,R,4.0E8)
      HMD=HPOL(HOUR,81.0,88.0,SAX,SUX,1.,1.)
      F(1)=HPOL(HOUR,0.02+0.03/DELA,0.05,SAX,SUX,1.,1.)
      F(2)=HPOL(HOUR,4.6,4.5,SAX,SUX,1.,1.)
      F(3)=HPOL(HOUR,-11.5,-4.0,SAX,SUX,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)
      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))
      D1=DXDX/(XDX*XKK*X**(XKK-1.0))
C
C SEARCH FOR HMF1 ..................................................
C
2726	IF(.not.BOTTO) GOTO 4933
	if(LAYVER) goto 6153
924	IF(.not.F1REG) GOTO 380
	XE2H=XE2(HEF)
      CALL REGFA1(HEF,HMF2,XE2H,NMF2,0.001,NMF1,XE2,SCHALT,HMF1)
	IF(.not.SCHALT) GOTO 380
	  WRITE(KONSOL,11)
11    FORMAT(1X,'*NE* HMF1 IS NOT EVALUATED BY THE FUNCTION XE2')
	IREGFA=1
c
c change B1 and try again ..........................................
c
9244 	IF(B1.GT.4.5) GOTO (7398,8922) IREGFA
	   	B1=B1+0.5
 		WRITE(KONSOL,902) B1-0.5,B1
902   FORMAT(6X,'CORR.: B1(OLD)=',F4.1,' B1(NEW)=',F4.1)
		IF(GULB0) then
			ib1=int(b1*2.-5.)
			B0=B0CNEW/b0b1(ib1)
			endif
   		GOTO 924
c
c omit F1 feature ....................................................
c
7398  WRITE(KONSOL,9269)
9269  FORMAT(1X,'CORR.: NO F1 REGION, B1=3, C1=0.0')
      	HMF1=0.
      	NMF1=0.
      	C1=0.0
      	B1=3.
      	F1REG=.FALSE.

C
C SEARCH FOR HST [NE3(HST)=NME] ..........................................
C
380	RRRR=0.5
	IF(F1REG) then
		hf1=hmf1
		xf1=nmf1
		GOTO 3972
		ENDIF
	RATHH=0.5
3973	hf1=hef+(hmf2-hef)*RATHH
	xf1=xe3(hf1)
	IF(XF1.LT.NME) THEN
		RATHH=RATHH+.1
		GOTO 3973
		ENDIF
3972	h=hf1
	deh=10.
	XXMIN=XF1
	HHMIN=HF1
3895    h=h-deh
	if(h.lt.HEF) then
	  h=h+2*deh
	  deh=deh/10.
	  if(deh.lt.1.) goto 3885
	  endif
   	XE3H=XE3(h)
	IF(XE3H.LT.XXMIN) then
	  XXMIN=XE3H
	  HHMIN=h
	  endif
	if(XE3H.gt.NME) goto 3895
      CALL REGFA1(h,HF1,XE3H,XF1,0.001,NME,XE3,SCHALT,HST)
	STR=HST
	IF(.not.SCHALT) GOTO 360
3885	WRITE(KONSOL,100)
100   FORMAT(1X,'*NE* HST IS NOT EVALUATED BY THE FUNCTION XE3')
	IREGFA=2
	IF(XXMIN/NME.LT.1.3) GOTO 9244
c
c assume linear interpolation between HZ and HEF ..................
c
8922    HZ=HHMIN+(HF1-HHMIN)*RRRR
        XNEHZ=XE3(HZ)
        if(xnehz-nme.lt.0.001) then
          RRRR=RRRR+.1
          GOTO 8922
          endif
        WRITE(KONSOL,901) HZ,HEF
901   FORMAT(6X,'CORR.: LIN. APP. BETWEEN HZ=',F5.1,
     &          ' AND HEF=',F5.1)
        T=(XNEHZ-NME)/(HZ-HEF)
        HST=-333.
        GOTO 4933
c
c calculate HZ, D and T ............................................
c
360	HZ=(HST+HF1)/2.0
    	D=HZ-HST
    	T=D*D/(HZ-HEF-D)
	GOTO 4933
C
C LAY-functions for middle ionosphere
C
6153	HMF1M=165.+0.6428*XHI
	HHALF = GRAT * HMF2
	HV1R = HME + WIDTH
	HV2R = HME + HDEEP
	HHMF2 = HMF2
	CALL INILAY(NIGHT,NMF2,NMF1,NME,VNER,HHMF2,HMF1M,HME,
     &			HV1R,HV2R,HHALF,HXL,SCL,AMP,IIQU)
	IF(IIQU.EQ.1) WRITE(KONSOL,7733)
7733	FORMAT('*NE* LAY amplitudes found with 2nd choice of HXL(1).')
	IF(IIQU.EQ.2) WRITE(KONSOL,7722)
7722	FORMAT('*NE* LAY amplitudes could not be found.')

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

4933  HTA=120.0
      HTE=3000.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
      SECNI=(24.-LONGI/15.)*3600.
      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) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      ATE(1)=TN120

C !!!!!!!!!! TE-MAXIMUM (JICAMARCA,ARECIBO) !!!!!!!!!!!!!!!!!!!!
      HMAXD=60.*EXP(-(MLAT/22.41)**2)+210.
      HMAXN=150.
      AHH(2)=HPOL(HOUR,HMAXD,HMAXN,SAX,SUX,1.,1.)
      TMAXD=800.*EXP(-(MLAT/33.)**2)+1500.
      TMAXN=TN(HMAXN,TEXNI,TLBDN,SIGNI)+20
      ATE(2)=HPOL(HOUR,TMAXD,TMAXN,SAX,SUX,1.,1.)

C !!!!!!!!!! TE(300,400KM)=TE-AE-C !!!!!!!!!!!!!!!!!!!!!!!!!!!!!
C !!!!!!!!!! TE(1400,3000KM)=TE-ISIS !!!!!!!!!!!!!!!!!!!!!!!!!!!
	DIPLAT=MAGBR
      CALL TEBA(DIPLAT,HOUR,NSESON,TEA)
      ATE(3)=TEA(1)
      ATE(4)=TEA(2)
      ATE(6)=TEA(3)
      ATE(7)=TEA(4)

C !!!!!!!!!! TE(600KM)=TE-AEROS !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      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,SAX,SUX,1.5,1.5)

C !!!!!!!!!! OPTION TO USE TE-NE-RELATION !!!!!!!!!!!!!!!!!!!!!!
C !!!!!!!!!! AT 300, 400 OR 600 KM  !!!!!!!!!!!!!!!!!!!!!!!!!!!!
      IF(TENEOP) THEN
        DO 3395 I=1,3
3395	  IF(TCON(I)) ATE(I+2)=TEDE(HOA(I),XNAR(I),-COV)
	ENDIF
C !!!!!!!!!! TE'S ARE CORRECTED !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
C !!!!!!!!!! ALSO 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
1901  STTE1=STTE2
C !!!!!!!!!! GRADIENTS ARE CALCULATED WITH !!!!!!!!!!!!!!!!!!!!
C !!!!!!!!!! CORRECTED REGION BOUNDARIES !!!!!!!!!!!!!!!!!!!!!!
      DO 1902 I=1,6
1902  STTE(I)=(ATE(I+1)-ATE(I))/(AHH(I+1)-AHH(I))
      ATE1=ATE(1)
887   CONTINUE
C
C------------ CALCULATION OF ION TEMPERATURE PARAMETERS--------
C
C !!!!!!!!!! TI(430KM,DAY)=TI-AEROS !!!!!!!!!!!!!!!!!!!!!!!!!!!
      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,SAX,SUX,1.,1.)
C !!!!!!!!!!  TI < TE   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        TED1=TEA(6)+30.
        IF(TID1.GT.TED1) TID1=TED1

C !!!!!!!!!! TI(430KM,NIGHT)=TI-AEROS !!!!!!!!!!!!!!!!!!!!!!!!!
      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 !!!!!!!!!! TN < TI < TE !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        TEN1=TEA(5)
	TNN1=TN(XSM1,TEXNI,TLBDN,SIGNI)
        IF(TEN1.LT.TNN1) TEN1=TNN1
        IF(TIN1.GT.TEN1) TIN1=TEN1
        IF(TIN1.LT.TNN1) TIN1=TNN1

C !!!!!!!!!! TI(430KM,LT) FROM STEP FUNCTION !!!!!!!!!!!!!!!!!!
	TI1=TIN1
 	IF(TID1.GT.TIN1) TI1=HPOL(HOUR,TID1,TIN1,SAX,SUX,1.,1.)

C !!!!!!!!!! TANGENT ON TN 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 ALTITTUDE 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.
      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(NIGHT) 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(ABS(XHI).LE.90.0) THEN
	ZZZ1=COS(XHI*UMR)
      ELSE
	ZZZ1=0.0
      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
7105  		MO2(L+1)=PF3O(7+I)+PF3O(8+I)*ZZZ1
      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
141   IF(.NOT.F1REG) HMF1=HZ

      DO 7397 KI=1,11
      DO 7397 KK=1,nz
7397    OUTF(KI,KK)=-1.

      do 7118 kk=1,nz
      height=zkm(kk)
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)
	ELSE
	  ELEDE=XE(HEIGHT)
	ENDIF
      	OUTF(1,KK)=ELEDE
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
7108  IF(NOION) GOTO 7118
      IF((HEIGHT.GT.HNIE).OR.(HEIGHT.LT.HNIA)) GOTO 7118
	if(DY) then
      zmonth = month + iday / 30.0
      call IONCOM(HEIGHT,XHI*UMR,LATI*UMR,COV,ZMONTH,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
      OUTF(5,KK)=ROX
      OUTF(6,KK)=RHX
      OUTF(7,KK)=RHEX
      OUTF(8,KK)=RO2X
      OUTF(9,KK)=RNOX
      OUTF(10,KK)=RNX
      OUTF(11,KK)=RCLUST

7118  continue
C
C ADDITIONAL PARAMETER FIELD OARR
C
        IF(NODEN) GOTO 6192
      OARR(1)=NMF2
      OARR(2)=HMF2
      OARR(3)=NMF1
      OARR(4)=HMF1
      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

3330  CONTINUE
      RETURN
      END
C
C
C Subroutine DFP, Stan Solomon, 3/92, splices filename to directory
C
      subroutine dfp(direct,filename,path)
      character*(*) direct,filename,path
      character*50 blanks
      data blanks/'                                                  '/
      path=blanks
      nch=len(direct)
      do 10 i=1,nch
      if (direct(i:i).ne.' ') goto 20
   10 continue
   20 lb=i
      do 30 i=nch,1,-1
      if (direct(i:i).ne.' ') goto 40
   30 continue
   40 le=i
      if (lb.ge.nch .or. le.le.0) then
        path(1:10)=filename(1:10)
      else
        nd=le-lb+1
        path(1:nd)=direct(lb:le)
        path(nd+1:nd+10)=filename(1:10)
      endif
      return
      end
C
C
C
C
C IRIF12.FOR ------------------------------------- OCTOBER 1991
C**************************************************************
c changes from IRIFU9 to IRIF10:
c       SOCO for solar zenith angle
c	ACOS and ASIN argument forced to be within -1 / +1
c	EPSTEIN functions corrected for large arguments
C**************************************************************
c changes from IRIF10 to IRIF11:
c	LAY subroutines introduced
c       TEBA corrected for 1400 km
C**************************************************************
c changes from IRIF11 to IRIF12:
C       Neutral temperature subroutines now in CIRA86.FOR
C       TEDER changed
C       All names with 6 or more characters replaced
C	10/29/91 XEN: 10^ in loop, instead of at the end
C
C**************************************************************
C********** INTERNATIONAL REFERENCE IONOSPHERE ****************
C**************************************************************
C****************  FUNCTIONS,SUBROUTINES  *********************
C**************************************************************
C** NE: 	XE1,DXE1N,XE2,XE3,XE4,XE5,XE6,XE
C** TE/TI: 	TEBA,SPHARM,ELTE,TEDE,TI,TEDER
C** NI:		RPID,RDHHE,RDNO,KOEFP1,KOEFP2,KOEFP3,SUFE
C** PEAKS:	F2OUT,HMF2ED,FOF1ED,FOEEDI,XMDED,GAMMA1
C** MAG. FIELD: GGM,FIELDG
C** FUNCTIONS: 	REGFA1,TAL
C** TIME:	SOCO,HPOL,MODA
C** INTERPOL.:	B0POL
C** EPSTEIN:	RLAY,D1LAY,D2LAY,EPTR,EPST,EPSTEP,EPLA
C** LAY:	XE2TO5,XEN,ROGUL,VALGUL,LNGLSN,LSKNM,INILAY
C** NI-new:	IONCOM, RPDA
C**************************************************************
C
C**************************************************************
C***  -------------------ADDRESSES------------------------  ***
C***  I  PROF. K. RAWER             DR. D. BILITZA       I  ***
C***  I  HERRENSTR. 43              GSFC CODE 933        I  ***
C***  I  7801 MARCH 1               GREENBELT MD 20771   I  ***
C***  I  F.R.G.                     USA                  I  ***
C***  ----------------------------------------------------  ***
C**************************************************************
C**************************************************************
C
C*************************************************************
C*************** ELECTRON DENSITY ****************************
C*************************************************************
C
C
      FUNCTION XE1(H)
c----------------------------------------------------------------
C REPRESENTING ELECTRON DENSITY(M-3) IN THE TOPSIDE IONOSPHERE
C (H=HMF2....1000 KM) BY HARMONIZED BENT-MODEL ADMITTING
C VARIABILITY OFGLOBAL PARAMETER ETA,ZETA,BETA,DELTA WITH
C GEOM. LATITUDE, SMOOTHED SOLAR FLUX AND CRITICAL FREQUENCY
C (SEE MAIN PROGRAM).
C [REF.:K.RAWER,S.RAMAKRISHNAN,1978]
c----------------------------------------------------------------
      COMMON	/BLOCK1/	HMF2,XNMF2,HMF1
     &		/BLO10/		BETA,ETA,DELTA,ZETA
     &		/ARGEXP/	ARGMAX

      DXDH = (1000.-HMF2)/700.
	x0 = 300. - delta
      xmx0 = (H-HMF2)/DXDH
         x = xmx0 + x0
	eptr1 = eptr(x,beta,394.5) - eptr(x0,beta,394.5)
	eptr2 = eptr(x,100.,300.0) - eptr(x0,100.,300.0)

      y = BETA * ETA * eptr1 + ZETA * (100. * eptr2 - xmx0)

 	Y = y * dxdh
	if(abs(Y).gt.argmax) Y = sign(argmax,Y)
      XE1 = XNMF2 * EXP(-Y)
      RETURN
      END
C
C
      FUNCTION DXE1N(H)
C LOGARITHMIC DERIVATIVE OF FUNCTION XE1 (KM-1).
      COMMON	/BLOCK1/	HMF2,XNMF2,HMF1
     &		/BLO10/		BETA,ETA,DELTA,ZETA

	x0 = 300. - delta
      X=(H-HMF2)/(1000.0-HMF2)*700.0 + x0
	epst2 = epst(x,100.0,300.0)
	epst1 = epst(x,beta ,394.5)
      DXE1N = - ETA * epst1 + ZETA * (1. - epst2)
      RETURN
      END
C
C
      REAL FUNCTION XE2(H)
C ELECTRON DENSITY FOR THE BOTTOMSIDE F-REGION (HMF1...HMF2).
      COMMON	/BLOCK1/HMF2,XNMF2,HMF1
     &		/BLOCK2/B0,B1,C1	/ARGEXP/ARGMAX
      X=(HMF2-H)/B0
	z=x**b1
	if(z.gt.argmax) z=argmax
      XE2=XNMF2*EXP(-z)/COSH(X)
      RETURN
      END
C
C
      REAL FUNCTION XE3(H)
C ELECTRON DENSITY FOR THE F1-LAYER (HZ.....HMF1).
      COMMON	/BLOCK1/	HMF2,XNMF2,HMF1
     &		/BLOCK2/	B0,B1,C1
      XE3=XE2(H)+XNMF2*C1*SQRT(ABS(HMF1-H)/B0)
      RETURN
      END
C
C
      REAL FUNCTION XE4(H)
C ELECTRON DENSITY FOR THE INDERMEDIUM REGION (HEF..HZ).
      COMMON	/BLOCK3/	HZ,T,HST,STR
     &		/BLOCK4/	HME,XNME,HEF
      IF(HST.LT.0.) GOTO 100
      XE4=XE3(HZ+T/2.0-SIGN(1.0,T)*SQRT(T*(HZ-H+T/4.0)))
      RETURN
100   XE4=XNME+T*(H-HEF)
      RETURN
      END
C
C
      REAL FUNCTION XE5(H)
C ELECTRON DENSITY FOR THE E AND VALLEY REGION (HME..HEF).
      LOGICAL NIGHT
      COMMON	/BLOCK4/	HME,XNME,HEF
     &  	/BLOCK5/	NIGHT,E(4)
      T3=H-HME
      T1=T3*T3*(E(1)+T3*(E(2)+T3*(E(3)+T3*E(4))))
      IF(NIGHT) GOTO 100
      XE5=XNME*(1+T1)
      RETURN
100   XE5=XNME*EXP(T1)
      RETURN
      END
C
C
      REAL FUNCTION XE6(H)
C ELECTRON DENSITY FOR THE D REGION (HA...HME).
      COMMON	/BLOCK4/	HME,XNME,HEF
     &		/BLOCK6/	HMD,XNMD,HDX
     &		/BLOCK7/	D1,XKK,FP30,FP3U,FP1,FP2
      IF(H.GT.HDX) GOTO 100
      Z=H-HMD
      FP3=FP3U
      IF(Z.GT.0.0) FP3=FP30
      XE6=XNMD*EXP(Z*(FP1+Z*(FP2+Z*FP3)))
      RETURN
100   Z=HME-H
      XE6=XNME*EXP(-D1*Z**XKK)
      RETURN
      END
C
C
      REAL FUNCTION XE(H)
C ELECTRON DENSITY BEETWEEN HA(KM) AND 1000 KM
C SUMMARIZING PROCEDURES  NE1....6;
      COMMON	/BLOCK1/HMF2,XNMF2,HMF1
     &		/BLOCK3/HZ,T,HST,STR
     &  	/BLOCK4/HME,XNME,HEF
      IF(H.LT.HMF2) GOTO 100
      XE=XE1(H)
      RETURN
100   IF(H.LT.HMF1) GOTO 300
      XE=XE2(H)
      RETURN
300   IF(H.LT.HZ) GOTO 400
      XE=XE3(H)
      RETURN
400   IF(H.LT.HEF) GOTO 500
      XE=XE4(H)
      RETURN
500   IF(H.LT.HME) GOTO 600
      XE=XE5(H)
      RETURN
600   XE=XE6(H)
      RETURN
      END
C
C**********************************************************
C***************** ELECTRON TEMPERATURE ********************
C**********************************************************
C
      SUBROUTINE TEBA(DIPL,SLT,NS,TE)
C CALCULATES ELECTRON TEMPERATURES TE(1) TO TE(4) AT ALTITUDES
C 300, 400, 1400 AND 3000 KM FOR DIP-LATITUDE DIPL/DEG AND
C LOCAL SOLAR TIME SLT/H USING THE BRACE-THEIS-MODELS (J. ATMOS.
C TERR. PHYS. 43, 1317, 1981); NS IS SEASON IN NORTHERN
C HEMISOHERE: IS=1 SPRING, IS=2 SUMMER ....
C ALSO CALCULATED ARE THE TEMPERATURES AT 400 KM ALTITUDE FOR
C MIDNIGHT (TE(5)) AND NOON (TE(6)).
      DIMENSION C(4,2,81),A(82),TE(6)
      COMMON/CONST/UMR
      DATA (C(1,1,J),J=1,81)/
     &.3100E1,-.3215E-2,.2440E+0,-.4613E-3,-.1711E-1,.2605E-1,
     &-.9546E-1,.1794E-1,.1270E-1,.2791E-1,.1536E-1,-.6629E-2,
     &-.3616E-2,.1229E-1,.4147E-3,.1447E-2,-.4453E-3,-.1853,
     &-.1245E-1,-.3675E-1,.4965E-2,.5460E-2,.8117E-2,-.1002E-1,
     &.5466E-3,-.3087E-1,-.3435E-2,-.1107E-3,.2199E-2,.4115E-3,
     &.6061E-3,.2916E-3,-.6584E-1,.4729E-2,-.1523E-2,.6689E-3,
     &.1031E-2,.5398E-3,-.1924E-2,-.4565E-1,.7244E-2,-.8543E-4,
     &.1052E-2,-.6696E-3,-.7492E-3,.4405E-1,.3047E-2,.2858E-2,
     &-.1465E-3,.1195E-2,-.1024E-3,.4582E-1,.8749E-3,.3011E-3,
     &.4473E-3,-.2782E-3,.4911E-1,-.1016E-1,.27E-2,-.9304E-3,
     &-.1202E-2,.2210E-1,.2566E-2,-.122E-3,.3987E-3,-.5744E-1,
     &.4408E-2,-.3497E-2,.83E-3,-.3536E-1,-.8813E-2,.2423E-2,
     &-.2994E-1,-.1929E-2,-.5268E-3,-.2228E-1,.3385E-2,
     &.413E-1,.4876E-2,.2692E-1,.1684E-2/
      DATA (C(1,2,J),J=1,81)/.313654E1,.6796E-2,.181413,.8564E-1,
     &-.32856E-1,-.3508E-2,-.1438E-1,-.2454E-1,.2745E-2,.5284E-1,
     &.1136E-1,-.1956E-1,-.5805E-2,.2801E-2,-.1211E-2,.4127E-2,
     &.2909E-2,-.25751,-.37915E-2,-.136E-1,-.13225E-1,.1202E-1,
     &.1256E-1,-.12165E-1,.1326E-1,-.7123E-1,.5793E-3,.1537E-2,
     &.6914E-2,-.4173E-2,.1052E-3,-.5765E-3,-.4041E-1,-.1752E-2,
     &-.542E-2,-.684E-2,.8921E-3,-.2228E-2,.1428E-2,.6635E-2,-.48045E-2,
     &-.1659E-2,-.9341E-3,.223E-3,-.9995E-3,.4285E-1,-.5211E-3,
     &-.3293E-2,.179E-2,.6435E-3,-.1891E-3,.3844E-1,.359E-2,-.8139E-3,
     &-.1996E-2,.2398E-3,.2938E-1,.761E-2,.347655E-2,.1707E-2,.2769E-3,
     &-.157E-1,.983E-3,-.6532E-3,.929E-4,-.2506E-1,.4681E-2,.1461E-2,
     &-.3757E-5,-.9728E-2,.2315E-2,.6377E-3,-.1705E-1,.2767E-2,
     &-.6992E-3,-.115E-1,-.1644E-2,.3355E-2,-.4326E-2,.2035E-1,.2985E-1/
      DATA (C(2,1,J),J=1,81)/.3136E1,.6498E-2,.2289,.1859E-1,-.3328E-1,
     &-.4889E-2,-.3054E-1,-.1773E-1,-.1728E-1,.6555E-1,.1775E-1,
     &-.2488E-1,-.9498E-2,.1493E-1,.281E-2,.2406E-2,.5436E-2,-.2115,
     &.7007E-2,-.5129E-1,-.7327E-2,.2402E-1,.4772E-2,-.7374E-2,
     &-.3835E-3,-.5013E-1,.2866E-2,.2216E-2,.2412E-3,.2094E-2,.122E-2
     &,-.1703E-3,-.1082,-.4992E-2,-.4065E-2,.3615E-2,-.2738E-2,
     &-.7177E-3,.2173E-3,-.4373E-1,-.375E-2,.5507E-2,-.1567E-2,
     &-.1458E-2,-.7397E-3,.7903E-1,.4131E-2,.3714E-2,.1073E-2,
     &-.8991E-3,.2976E-3,.2623E-1,.2344E-2,.5608E-3,.4124E-3,.1509E-3,
     &.5103E-1,.345E-2,.1283E-2,.7238E-3,-.3464E-4,.1663E-1,-.1644E-2,
     &-.71E-3,.5281E-3,-.2729E-1,.3556E-2,-.3391E-2,-.1787E-3,.2154E-2,
     &.6476E-2,-.8282E-3,-.2361E-1,.9557E-3,.3205E-3,-.2301E-1,
     &-.854E-3,-.1126E-1,-.2323E-2,-.8582E-2,.2683E-1/
      DATA (C(2,2,J),J=1,81)/.3144E1,.8571E-2,.2539,.6937E-1,-.1667E-1,
     &.2249E-1,-.4162E-1,.1201E-1,.2435E-1,.5232E-1,.2521E-1,-.199E-1,
     &-.7671E-2,.1264E-1,-.1551E-2,-.1928E-2,.3652E-2,-.2019,.5697E-2,
     &-.3159E-1,-.1451E-1,.2868E-1,.1377E-1,-.4383E-2,.1172E-1,
     &-.5683E-1,.3593E-2,.3571E-2,.3282E-2,.1732E-2,-.4921E-3,-.1165E-2
     &,-.1066,-.1892E-1,.357E-2,-.8631E-3,-.1876E-2,-.8414E-4,.2356E-2,
     &-.4259E-1,-.322E-2,.4641E-2,.6223E-3,-.168E-2,-.1243E-3,.7393E-1,
     &-.3143E-2,-.2362E-2,.1235E-2,-.1551E-2,.2099E-3,.2299E-1,.5301E-2
     &,-.4306E-2,-.1303E-2,.7687E-5,.5305E-1,.6642E-2,-.1686E-2,
     &.1048E-2,.5958E-3,.4341E-1,-.8819E-4,-.333E-3,-.2158E-3,-.4106E-1
     &,.4191E-2,.2045E-2,-.1437E-3,-.1803E-1,-.8072E-3,-.424E-3,
     &-.26E-1,-.2329E-2,.5949E-3,-.1371E-1,-.2188E-2,.1788E-1,
     &.6405E-3,.5977E-2,.1333E-1/
      DATA (C(3,1,J),J=1,81)/.3372E1,.1006E-1,.1436,.2023E-2,-.5166E-1,
     &.9606E-2,-.5596E-1,.4914E-3,-.3124E-2,-.4713E-1,-.7371E-2,
     &-.4823E-2,-.2213E-2,.6569E-2,-.1962E-3,.3309E-3,-.3908E-3,
     &-.2836,.7829E-2,.1175E-1,.9919E-3,.6589E-2,.2045E-2,-.7346E-2
     &,-.89E-3,-.347E-1,-.4977E-2,.147E-2,-.2823E-5,.6465E-3,
     &-.1448E-3,.1401E-2,-.8988E-1,-.3293E-4,-.1848E-2,.4439E-3,
     &-.1263E-2,.317E-3,-.6227E-3,.1721E-1,-.199E-2,-.4627E-3,
     &.2897E-5,-.5454E-3,.3385E-3,.8432E-1,-.1951E-2,.1487E-2,
     &.1042E-2,-.4788E-3,-.1276E-3,.2373E-1,.2409E-2,.5263E-3,
     &.1301E-2,-.4177E-3,.3974E-1,.1418E-3,-.1048E-2,-.2982E-3,
     &-.3396E-4,.131E-1,.1413E-2,-.1373E-3,.2638E-3,-.4171E-1,
     &-.5932E-3,-.7523E-3,-.6883E-3,-.2355E-1,.5695E-3,-.2219E-4,
     &-.2301E-1,-.9962E-4,-.6761E-3,.204E-2,-.5479E-3,.2591E-1,
     &-.2425E-2,.1583E-1,.9577E-2/
      DATA (C(3,2,J),J=1,81)/.3367E1,.1038E-1,.1407,.3622E-1,-.3144E-1,
     &.112E-1,-.5674E-1,.3219E-1,.1288E-2,-.5799E-1,-.4609E-2,
     &.3252E-2,-.2859E-3,.1226E-1,-.4539E-2,.1310E-2,-.5603E-3,
     &-.311,-.1268E-2,.1539E-1,.3146E-2,.7787E-2,-.143E-2,-.482E-2
     &,.2924E-2,-.9981E-1,-.7838E-2,-.1663E-3,.4769E-3,.4148E-2,
     &-.1008E-2,-.979E-3,-.9049E-1,-.2994E-2,-.6748E-2,-.9889E-3,
     &.1488E-2,-.1154E-2,-.8412E-4,-.1302E-1,-.4859E-2,-.7172E-3,
     &-.9401E-3,.9101E-3,-.1735E-3,.7055E-1,.6398E-2,-.3103E-2,
     &-.938E-3,-.4E-3,-.1165E-2,.2713E-1,-.1654E-2,.2781E-2,
     &-.5215E-5,.2258E-3,.5022E-1,.95E-2,.4147E-3,.3499E-3,
     &-.6097E-3,.4118E-1,.6556E-2,.3793E-2,-.1226E-3,-.2517E-1,
     &.1491E-3,.1075E-2,.4531E-3,-.9012E-2,.3343E-2,.3431E-2,
     &-.2519E-1,.3793E-4,.5973E-3,-.1423E-1,-.132E-2,-.6048E-2,
     &-.5005E-2,-.115E-1,.2574E-1/
      DATA (C(4,1,J),J=1,81)/.3574E1,.0,.7537E-1,.0,-.8459E-1,
     &0.,-.294E-1,0.,.4547E-1,-.5321E-1,0.,.4328E-2,0.,.6022E-2,
     &.0,-.9168E-3,.0,-.1768,.0,.294E-1,.0,.5902E-3,.0,-.9047E-2,
     &.0,-.6555E-1,.0,-.1033E-2,.0,.1674E-2,.0,.2802E-3,-.6786E-1
     &,.0,.4193E-2,.0,-.6448E-3,.0,.9277E-3,-.1634E-1,.0,-.2531E-2
     &,.0,.193E-4,.0,.528E-1,.0,.2438E-2,.0,-.5292E-3,.0,.1555E-1
     &,.0,-.3259E-2,.0,-.5998E-3,.3168E-1,.0,.2382E-2,.0,-.4078E-3
     &,.2312E-1,.0,.1481E-3,.0,-.1885E-1,.0,.1144E-2,.0,-.9952E-2
     &,.0,-.551E-3,-.202E-1,.0,-.7283E-4,-.1272E-1,.0,.2224E-2,
     &.0,-.251E-2,.2434E-1/
      DATA (C(4,2,J),J=1,81)/.3574E1,-.5639E-2,.7094E-1,
     &-.3347E-1,-.861E-1,-.2877E-1,-.3154E-1,-.2847E-2,.1235E-1,
     &-.5966E-1,-.3236E-2,.3795E-3,-.8634E-3,.3377E-2,-.1071E-3,
     &-.2151E-2,-.4057E-3,-.1783,.126E-1,.2835E-1,-.242E-2,
     &.3002E-2,-.4684E-2,-.6756E-2,-.7493E-3,-.6147E-1,-.5636E-2
     &,-.1234E-2,-.1613E-2,-.6353E-4,-.2503E-3,-.1729E-3,-.7148E-1
     &,.5326E-2,.4006E-2,.6484E-3,-.1046E-3,-.6034E-3,-.9435E-3,
     &-.2385E-2,.6853E-2,.151E-2,.1319E-2,.9049E-4,-.1999E-3,
     &.3976E-1,.2802E-2,-.103E-2,.5599E-3,-.4791E-3,-.846E-4,
     &.2683E-1,.427E-2,.5911E-3,.2987E-3,-.208E-3,.1396E-1,
     &-.1922E-2,-.1063E-2,.3803E-3,.1343E-3,.1771E-1,-.1038E-2,
     &-.4645E-3,-.2481E-3,-.2251E-1,-.29E-2,-.3977E-3,-.516E-3,
     &-.8079E-2,-.1528E-2,.306E-3,-.1582E-1,-.8536E-3,.1565E-3,
     &-.1252E-1,.2319E-3,.4311E-2,.1024E-2,.1296E-5,.179E-1/
	IF(NS.LT.3) THEN
	   IS=NS
	ELSE IF(NS.GT.3) THEN
	   IS=2
	   DIPL=-DIPL
	ELSE
	   IS=1
	ENDIF
      COLAT=UMR*(90.-DIPL)
      AZ=.2618*SLT
      CALL SPHARM(A,8,8,COLAT,AZ)
	IF(IS.EQ.2) THEN
	   KEND=3
	ELSE
	   KEND=4
	ENDIF
      DO 2 K=1,KEND
      STE=0.
      DO 1 I=1,81
1       STE=STE+A(I)*C(K,IS,I)
2     TE(K)=10.**STE
	IF(IS.EQ.2) THEN
	   DIPL=-DIPL
	   COLAT=UMR*(90.-DIPL)
           CALL SPHARM(A,8,8,COLAT,AZ)
           STE=0.
           DO 11 I=1,81
11            STE=STE+A(I)*C(4,2,I)
           TE(4)=10.**STE
    	ENDIF

C---------- TEMPERATURE AT 400KM AT MIDNIGHT AND NOON
      DO 4 J=1,2
        STE=0.
        AZ=.2618*(J-1)*12.
        CALL SPHARM(A,8,8,COLAT,AZ)
        DO 3 I=1,81
3         STE=STE+A(I)*C(2,IS,I)
4       TE(J+4)=10.**STE
      RETURN
      END
C
      SUBROUTINE SPHARM(C,L,M,COLAT,AZ)
C CALCULATES THE COEFFICIENTS OF THE SPHERICAL HARMONIC
C EXPANSION THAT WAS USED FOR THE BRACE-THEIS-MODELS.
      DIMENSION C(82)
      C(1)=1.
      K=2
      X=COS(COLAT)
      C(K)=X
      K=K+1
      DO 10 I=2,L
      C(K)=((2*I-1)*X*C(K-1)-(I-1)*C(K-2))/I
10    K=K+1
      Y=SIN(COLAT)
      DO 20 MT=1,M
      CAZ=COS(MT*AZ)
      SAZ=SIN(MT*AZ)
      C(K)=Y**MT
      K=K+1
      IF(MT.EQ.L) GOTO 16
      C(K)=C(K-1)*X*(2*MT+1)
      K=K+1
      IF((MT+1).EQ.L) GOTO 16
      DO 15 I=2+MT,L
      C(K)=((2*I-1)*X*C(K-1)-(I+MT-1)*C(K-2))/(I-MT)
15    K=K+1
16    N=L-MT+1
      DO 18 I=1,N
      C(K)=C(K-N)*CAZ
      C(K-N)=C(K-N)*SAZ
18    K=K+1
20    CONTINUE
      RETURN
      END
C
C
      REAL FUNCTION ELTE(H)
c----------------------------------------------------------------
C ELECTRON TEMPERATURE PROFILE BASED ON THE TEMPERATURES AT 120
C HMAX,300,400,600,1400,3000 KM ALTITUDE. INBETWEEN CONSTANT
C GRADIENT IS ASSUMED. ARGMAX IS MAXIMUM ARGUMENT ALLOWED FOR
C EXP-FUNCTION.
c----------------------------------------------------------------
      COMMON /BLOTE/AH(7),ATE1,ST(6),D(5)
C
      SUM=ATE1+ST(1)*(H-AH(1))
      DO 1 I=1,5
	aa = eptr(h    ,d(i),ah(i+1))
	bb = eptr(ah(1),d(i),ah(i+1))
1     SUM=SUM+(ST(I+1)-ST(I))*(AA-BB)*D(I)
      ELTE=SUM
      RETURN
      END
C
C
      FUNCTION TEDE(H,DEN,COV)
C ELECTRON TEMEPERATURE MODEL AFTER BRACE,THEIS .
C FOR NEG. COV THE MEAN COV-INDEX (3 SOLAR ROT.) IS EXPECTED.
C DEN IS THE ELECTRON DENSITY IN M-3.
      Y=1051.+(17.01*H-2746.)*
     &EXP(-5.122E-4*H+(6.094E-12-3.353E-14*H)*DEN)
      ACOV=ABS(COV)
      YC=1.+(.117+2.02E-3*ACOV)/(1.+EXP(-(ACOV-102.5)/5.))
      IF(COV.LT.0.)
     &YC=1.+(.123+1.69E-3*ACOV)/(1.+EXP(-(ACOV-115.)/10.))
      TEDE=Y*YC
      RETURN
      END
C
C
C*************************************************************
C**************** ION TEMPERATURE ****************************
C*************************************************************
C
C
      REAL FUNCTION TI(H)
c----------------------------------------------------------------
C ION TEMPERATURE FOR HEIGHTS NOT GREATER 1000 KM AND NOT LESS HS
C EXPLANATION SEE FUNCTION RPID.
c----------------------------------------------------------------
      REAL 		MM
      COMMON  /BLOCK8/  HS,TNHS,XSM(4),MM(5),G(4),M

      SUM=MM(1)*(H-HS)+TNHS
      DO 100 I=1,M-1
	aa = eptr(h ,g(i),xsm(i))
	bb = eptr(hs,g(i),xsm(i))
100   	SUM=SUM+(MM(I+1)-MM(I))*(AA-BB)*G(I)
      TI=SUM
      RETURN
      END
C
C
      REAL FUNCTION TEDER(H)
C THIS FUNCTION ALONG WITH PROCEDURE REGFA1 ALLOWS TO FIND
C THE  HEIGHT ABOVE WHICH TN BEGINS TO BE DIFFERENT FROM TI
      COMMON	/BLOTN/XSM1,TEX,TLBD,SIG
      TNH = TN(H,TEX,TLBD,SIG)
      DTDX = DTNDH(H,TEX,TLBD,SIG)
      TEDER = DTDX * ( XSM1 - H ) + TNH
      RETURN
      END
C
C
C*************************************************************
C************* ION RELATIVE PRECENTAGE DENSITY *****************
C*************************************************************
C
C
      REAL FUNCTION RPID (H, H0, N0, M, ST, ID, XS)
c------------------------------------------------------------------
C D.BILITZA,1977,THIS ANALYTIC FUNCTION IS USED TO REPRESENT THE
C RELATIVE PRECENTAGE DENSITY OF ATOMAR AND MOLECULAR OXYGEN IONS.
C THE M+1 HEIGHT GRADIENTS ST(M+1) ARE CONNECTED WITH EPSTEIN-
C STEP-FUNCTIONS AT THE STEP HEIGHTS XS(M) WITH TRANSITION
C THICKNESSES ID(M). RPID(H0,H0,N0,....)=N0.
C ARGMAX is the highest allowed argument for EXP in your system.
c------------------------------------------------------------------
      REAL, intent(in) :: H, H0, N0
      integer, intent(in) :: M, ID(M)
      real, intent(in) :: ST(M+1), XS(M)
      COMMON  /ARGEXP/	ARGMAX

      SUM=(H-H0)*ST(1)
      DO 100  I=1,M
	      XI=ID(I)
		aa = eptr(h ,xi,xs(i))
		bb = eptr(h0,xi,xs(i))
100	      SUM=SUM+(ST(I+1)-ST(I))*(AA-BB)*XI
      IF(ABS(SUM).LT.ARGMAX) then
	SM=EXP(SUM)
      else IF(SUM.Gt.0.0) then
	SM=EXP(ARGMAX)
      else
	SM=0.0
      endif
      RPID= n0 * SM
      RETURN
      END
C
c
      SUBROUTINE RDHHE (H,HB,RDOH,RDO2H,RNO,PEHE,RDH,RDHE)
C BILITZA,FEB.82,H+ AND HE+ RELATIVE PERECENTAGE DENSITY BELOW
C 1000 KM. THE O+ AND O2+ REL. PER. DENSITIES SHOULD BE GIVEN
C (RDOH,RDO2H). HB IS THE ALTITUDE OF MAXIMAL O+ DENSITY. PEHE
C IS THE PRECENTAGE OF HE+ IONS COMPARED TO ALL LIGHT IONS.
C RNO IS THE RATIO OF NO+ TO O2+DENSITY AT H=HB.
      RDHE=0.0
      RDH=0.0
      IF(H.LE.HB) GOTO 100
      REST=100.0-RDOH-RDO2H-RNO*RDO2H
      RDH=REST*(1.-PEHE/100.)
      RDHE=REST*PEHE/100.
100   RETURN
      END
C
C
      REAL FUNCTION RDNO(H,HB,RDO2H,RDOH,RNO)
C D.BILITZA, 1978. NO+ RELATIVE PERCENTAGE DENSITY ABOVE 100KM.
C FOR MORE INFORMATION SEE SUBROUTINE RDHHE.
      IF (H.GT.HB) GOTO 200
      RDNO=100.0-RDO2H-RDOH
      RETURN
200   RDNO=RNO*RDO2H
      RETURN
      END
C
C
      SUBROUTINE  KOEFP1(PG1O)
C THIEMANN,1979,COEFFICIENTS PG1O FOR CALCULATING  O+ PROFILES
C BELOW THE F2-MAXIMUM. CHOSEN TO APPROACH DANILOV-
C SEMENOV'S COMPILATION.
      DIMENSION PG1O(80)
      REAL FELD (80)
      DATA FELD/-11.0,-11.0,4.0,-11.0,0.08018,
     &0.13027,0.04216,0.25  ,-0.00686,0.00999,
     &5.113,0.1 ,170.0,180.0,0.1175,0.15,-11.0,
     &1.0 ,2.0,-11.0,0.069,0.161,0.254,0.18,0.0161,
     &0.0216,0.03014,0.1,152.0,167.0,0.04916,
     &0.17,-11.0,2.0,2.0,-11.0,0.072,0.092,0.014,0.21,
     &0.01389,0.03863,0.05762,0.12,165.0,168.0,0.008,
     &0.258,-11.0,1.0,3.0,-11.0,0.091,0.088,
     &0.008,0.34,0.0067,0.0195,0.04,0.1,158.0,172.0,
     &0.01,0.24,-11.0,2.0,3.0, -11.0,0.083,0.102,
     &0.045,0.03,0.00127,0.01,0.05,0.09,167.0,185.0,
     &0.015,0.18/
      K=0
      DO 10 I=1,80
      K=K+1
10    PG1O(K)=FELD(I)
      RETURN
      END
C
C
      SUBROUTINE KOEFP2(PG2O)
C THIEMANN,1979,COEFFICIENTS FOR CALCULATION OF O+ PROFILES
C ABOVE THE F2-MAXIMUM (DUMBS,SPENNER:AEROS-COMPILATION)
      DIMENSION PG2O(32)
      REAL FELD(32)
      DATA FELD/1.0,-11.0,-11.0,1.0,695.0,-.000781,
     &-.00264,2177.0,1.0,-11.0,-11.0,2.0,570.0,
     &-.002,-.0052,1040.0,2.0,-11.0,-11.0,1.0,695.0,
     &-.000786,-.00165,3367.0,2.0,-11.0,-11.0,2.0,
     &575.0,-.00126,-.00524,1380.0/
      K=0
      DO 10 I=1,32
      K=K+1
10    PG2O(K)=FELD(I)
      RETURN
      END
C
C
      SUBROUTINE  KOEFP3(PG3O)
C THIEMANN,1979,COEFFICIENTS FOR CALCULATING O2+ PROFILES.
C CHOSEN AS TO APPROACH DANILOV-SEMENOV'S COMPILATION.
      DIMENSION PG3O(80)
      REAL FELD(80)
      DATA FELD/-11.0,1.0,2.0,-11.0,160.0,31.0,130.0,
     &-10.0,198.0,0.0,0.05922,-0.07983,
     &-0.00397,0.00085,-0.00313,0.0,-11.0,2.0,2.0,-11.0,
     &140.0,30.0,130.0,-10.0,
     &190.0,0.0,0.05107,-0.07964,0.00097,-0.01118,-0.02614,
     &-0.09537,
     &-11.0,1.0,3.0,-11.0,140.0,37.0,125.0,0.0,182.0,
     &0.0,0.0307,-0.04968,-0.00248,
     &-0.02451,-0.00313,0.0,-11.0,2.0,3.0,-11.0,
     &140.0,37.0,125.0,0.0,170.0,0.0,
     &0.02806,-0.04716,0.00066,-0.02763,-0.02247,-0.01919,
     &-11.0,-11.0,4.0,-11.0,140.0,45.0,136.0,-9.0,
     &181.0,-26.0,0.02994,-0.04879,
     &-0.01396,0.00089,-0.09929,0.05589/
      K=0
      DO 10 I=1,80
      K=K+1
10    PG3O(K)=FELD(I)
      RETURN
      END
C
C
      SUBROUTINE SUFE (FIELD,RFE,M,FE)
C SELECTS THE REQUIRED ION DENSITY PARAMETER SET.
C THE INPUT FIELD INCLUDES DIFFERENT SETS OF DIMENSION M EACH
C CARACTERISED BY 4 HEADER NUMBERS. RFE(4) SHOULD CONTAIN THE
C CHOSEN HEADER NUMBERS.FE(M) IS THE CORRESPONDING SET.
      integer, intent(in) :: M
      real, intent(inout) :: FE(M)
      real, intent(in) :: RFE(4),FIELD(4*M)
      real :: EFE(4)
      K=0
100   DO 101 I=1,4
      K=K+1
101   EFE(I)=FIELD(K)
      DO 111 I=1,M
      K=K+1
111   FE(I)=FIELD(K)
      DO 120 I=1,4
      IF((EFE(I).GT.-10.0).AND.(RFE(I).NE.EFE(I))) GOTO 100
120   CONTINUE
      END
C
C
C*************************************************************
C************* PEAK VALUES ELECTRON DENSITY ******************
C*************************************************************
C
C
      SUBROUTINE F2OUT(XMODIP,XLATI,XLONGI,FF0,XM0,UT,
     &                              FOF2,XM3000)
C CALCULATES FOF2/MHZ AND M3000 USING THE CCIR-MAPS.
C INPUT: MODIFIED DIP LATITUDE XMODIP, GEOG. LATITUDE XLATI,
C LONGITUDE XLONGI (ALL IN DEG.), SMOOTHED SUNSPOT NUMBER R,
C MONTH AND UNIVERSAL TIME UT (DEC. HOURS).
C D.BILITZA,JULY 85.
      DIMENSION FF0(988),XM0(441)
      INTEGER QM(7),QF(9)
      DATA QF/11,11,8,4,1,0,0,0,0/,QM/6,7,5,2,1,0,0/
      FOF2=GAMMA1(XMODIP,XLATI,XLONGI,UT,6,QF,9,76,13,988,FF0)
      XM3000=GAMMA1(XMODIP,XLATI,XLONGI,UT,4,QM,7,49,9,441,XM0)
      RETURN
      END
C
C
      REAL FUNCTION HMF2ED(XMAGBR,R,X,XM3)
C CALCULATES THE PEAK HEIGHT HMF2/KM FOR THE MAGNETIC
C LATITUDE XMAGBR/DEG. AND THE SMOOTHED ZUERICH SUNSPOT
C NUMBER R USING CCIR-M3000 XM3 AND THE RATIO X=FOF2/FOE.
C [REF. D.BILITZA ET AL., TELECOMM.J., 46, 549-553, 1979]
C D.BILITZA,1980.
      F1=(2.32E-3)*R+0.222
      F2=1.2-(1.16E-2)*EXP((2.39E-2)*R)
      F3=0.096*(R-25.0)/150.0
      DELM=F1*(1.0-R/150.0*EXP(-XMAGBR*XMAGBR/1600.0))/(X-F2)+F3
      HMF2ED=1490.0/(XM3+DELM)-176.0
      RETURN
      END
C
C
      REAL FUNCTION FOF1ED(YLATI,R,CHI)
c--------------------------------------------------------------
C CALCULATES THE F1 PEAK PLASMA FREQUENCY (FOF1/MHZ)
C FOR 	DIP-LATITUDE (YLATI/DEGREE)
c	SMOOTHED ZURICH SUNSPOT NUMBER (R)
c	SOLAR ZENITH ANGLE (CHI/DEGREE)
C REFERENCE:
c	E.D.DUCHARME ET AL., RADIO SCIENCE 6, 369-378, 1971
C 				       AND 8, 837-839, 1973
c	HOWEVER WITH MAGNETIC DIP LATITUDE INSTEAD OF GEOMAGNETIC
c	DIPOLE LATITUDE, EYFRIG, 1979
C--------------------------------------------- D. BILITZA, 1988.
	COMMON/CONST/UMR
	FOF1 = 0.0
	DLA =  YLATI
		CHI0 = 49.84733 + 0.349504 * DLA
		CHI100 = 38.96113 + 0.509932 * DLA
		CHIM = ( CHI0 + ( CHI100 - CHI0 ) * R / 100. )
		IF(CHI.GT.CHIM) GOTO 1
  	F0 = 4.35 + DLA * ( 0.0058 - 1.2E-4 * DLA )
	F100 = 5.348 + DLA * ( 0.011 - 2.3E-4 * DLA )
	FS = F0 + ( F100 - F0 ) * R / 100.0
	XMUE = 0.093 + DLA * ( 0.0046 - 5.4E-5 * DLA ) + 3.0E-4 * R
	FOF1 = FS * COS( CHI * UMR ) ** XMUE
1	FOF1ED = FOF1
	RETURN
	END
C
C
      REAL FUNCTION FOEEDI(COV,XHI,XHIM,XLATI)
C-------------------------------------------------------
C CALCULATES FOE/MHZ BY THE EDINBURGH-METHOD.
C INPUT: MEAN 10.7CM SOLAR RADIO FLUX (COV), GEOGRAPHIC
C LATITUDE (XLATI/DEG), SOLAR ZENITH ANGLE (XHI/DEG AND
C XHIM/DEG AT NOON).
C REFERENCE:
C 	KOURIS-MUGGELETON, CCIR DOC. 6/3/07, 1973
C 	TROST, J. GEOPHYS. RES. 84, 2736, 1979 (was used
C		to improve the nighttime varition)
C D.BILITZA--------------------------------- AUGUST 1986.
      COMMON/CONST/UMR
C variation with solar activity (factor A) ...............
      A=1.0+0.0094*(COV-66.0)
C variation with noon solar zenith angle (B) and with latitude (C)
      SL=COS(XLATI*UMR)
	IF(XLATI.LT.32.0) THEN
		SM=-1.93+1.92*SL
		C=23.0+116.0*SL
 	ELSE
	 	SM=0.11-0.49*SL
	 	C=92.0+35.0*SL
  	ENDIF
	if(XHIM.ge.90.) XHIM=89.999
	B = COS(XHIM*UMR) ** SM
C variation with solar zenith angle (D) ..........................
 	IF(XLATI.GT.12.0) THEN
		SP=1.2
	ELSE
		SP=1.31
	ENDIF
C adjusted solar zenith angle during nighttime (XHIC) .............
      XHIC=XHI-3.*ALOG(1.+EXP((XHI-89.98)/3.))
      D=COS(XHIC*UMR)**SP
C determine foE**4 ................................................
      R4FOE=A*B*C*D
C minimum allowable foE (sqrt[SMIN])...............................
      SMIN=0.121+0.0015*(COV-60.)
      SMIN=SMIN*SMIN
      IF(R4FOE.LT.SMIN) R4FOE=SMIN
      FOEEDI=R4FOE**0.25
      RETURN
      END
C
C
      REAL FUNCTION XMDED(XHI,R,YW)
C D. BILITZA, 1978, CALCULATES ELECTRON DENSITY OF D MAXIMUM.
C XHI/DEG. IS SOLAR ZENITH ANGLE, R SMOOTHED ZURICH SUNSPOT NUMBER
C AND YW/M-3 THE ASSUMED CONSTANT NIGHT VALUE.
C [REF.: D.BILITZA, WORLD DATA CENTER A REPORT UAG-82,7,
C       BOULDER,1981]
      COMMON/CONST/UMR
      Y=6.05E8+0.088E8*R
      Z=(-0.1/(ALOG(YW/Y)))**0.3704
	if(abs(z).gt.1.) z=sign(1.,z)
      SUXHI=ACOS(Z)
      IF (SUXHI.LT.1.0472) SUXHI=1.0472
      XXHI=XHI*UMR
      IF (XXHI.GT.SUXHI) GOTO 100
      X=COS(XXHI)
      XMDED=Y*EXP(-0.1/X**2.7)
      RETURN
100   XMDED=YW
      RETURN
      END
C
C
      REAL FUNCTION GAMMA1(SMODIP,SLAT,SLONG,HOUR,IHARM,NQ,
     &  			K1,M,MM,M3,SFE)
C CALCULATES GAMMA1=FOF2 OR M3000 USING CCIR NUMERICAL MAP
C COEFFICIENTS SFE(M3) FOR MODIFIED DIP LATITUDE (SMODIP/DEG)
C GEOGRAPHIC LATITUDE (SLAT/DEG) AND LONGITUDE (SLONG/DEG)
C AND UNIVERSIAL TIME (HOUR/DECIMAL HOURS).
C NQ(K1) IS AN INTEGER ARRAY GIVING THE HIGHEST DEGREES IN
C LATITUDE FOR EACH LONGITUDE HARMONIC.
C M=1+NQ1+2(NQ2+1)+2(NQ3+1)+... .
C SHEIKH,4.3.77.
      REAL*8 C(12),S(12),COEF(100),SUM
      DIMENSION NQ(K1),XSINX(13),SFE(M3)
      COMMON/CONST/UMR
      HOU=(15.0*HOUR-180.0)*UMR
      S(1)=SIN(HOU)
      C(1)=COS(HOU)
      DO 250 I=2,IHARM
      C(I)=C(1)*C(I-1)-S(1)*S(I-1)
      S(I)=C(1)*S(I-1)+S(1)*C(I-1)
250   CONTINUE
      DO 300 I=1,M
      MI=(I-1)*MM
      COEF(I)=SFE(MI+1)
      DO 300 J=1,IHARM
      COEF(I)=COEF(I)+SFE(MI+2*J)*S(J)+SFE(MI+2*J+1)*C(J)
300   CONTINUE
      SUM=COEF(1)
      SS=SIN(SMODIP*UMR)
      S3=SS
      XSINX(1)=1.0
      INDEX=NQ(1)
      DO 350 J=1,INDEX
      SUM=SUM+COEF(1+J)*SS
      XSINX(J+1)=SS
      SS=SS*S3
350   CONTINUE
      XSINX(NQ(1)+2)=SS
      NP=NQ(1)+1
      SS=COS(SLAT*UMR)
      S3=SS
      DO 400 J=2,K1
      S0=SLONG*(J-1.)*UMR
      S1=COS(S0)
      S2=SIN(S0)
      INDEX=NQ(J)+1
      DO 450 L=1,INDEX
      NP=NP+1
      SUM=SUM+COEF(NP)*XSINX(L)*SS*S1
      NP=NP+1
      SUM=SUM+COEF(NP)*XSINX(L)*SS*S2
450   CONTINUE
      SS=SS*S3
400   CONTINUE
      GAMMA1=SUM
      RETURN
      END
C
C
C************************************************************
C*************** EARTH MAGNETIC FIELD ***********************
C**************************************************************
C
C
      SUBROUTINE GGM(ART,LONG,LATI,MLONG,MLAT)
C CALCULATES GEOMAGNETIC LONGITUDE (MLONG) AND LATITUDE (MLAT)
C FROM GEOGRAFIC LONGITUDE (LONG) AND LATITUDE (LATI) FOR ART=0
C AND REVERSE FOR ART=1. ALL ANGLES IN DEGREE.
C LATITUDE:-90 TO 90. LONGITUDE:0 TO 360 EAST.
      INTEGER ART
      REAL MLONG,MLAT,LONG,LATI
      COMMON/CONST/FAKTOR
      ZPI=FAKTOR*360.
      CBG=11.4*FAKTOR
      CI=COS(CBG)
      SI=SIN(CBG)
      IF(ART.EQ.0) GOTO 10
      CBM=COS(MLAT*FAKTOR)
      SBM=SIN(MLAT*FAKTOR)
      CLM=COS(MLONG*FAKTOR)
      SLM=SIN(MLONG*FAKTOR)
      SBG=SBM*CI-CBM*CLM*SI
        IF(ABS(SBG).GT.1.) SBG=SIGN(1.,SBG)
      LATI=ASIN(SBG)
      CBG=COS(LATI)
      SLG=(CBM*SLM)/CBG
      CLG=(SBM*SI+CBM*CLM*CI)/CBG
        IF(ABS(CLG).GT.1.) CLG=SIGN(1.,CLG)
      LONG=ACOS(CLG)
      IF(SLG.LT.0.0) LONG=ZPI-LONG
      LATI=LATI/FAKTOR
      LONG=LONG/FAKTOR
      LONG=LONG-69.8
      IF(LONG.LT.0.0) LONG=LONG+360.0
      RETURN
10    YLG=LONG+69.8
      CBG=COS(LATI*FAKTOR)
      SBG=SIN(LATI*FAKTOR)
      CLG=COS(YLG*FAKTOR)
      SLG=SIN(YLG*FAKTOR)
      SBM=SBG*CI+CBG*CLG*SI
        IF(ABS(SBM).GT.1.) SBM=SIGN(1.,SBM)
      MLAT=ASIN(SBM)
      CBM=COS(MLAT)
      SLM=(CBG*SLG)/CBM
      CLM=(-SBG*SI+CBG*CLG*CI)/CBM
        IF(ABS(CLM).GT.1.) CLM=SIGN(1.,CLM)
      MLONG=ACOS(CLM)
      IF(SLM.LT..0) MLONG=ZPI-MLONG
      MLAT=MLAT/FAKTOR
      MLONG=MLONG/FAKTOR
      RETURN
      END
C
C
      SUBROUTINE FIELDG(DLAT,DLONG,ALT,X,Y,Z,F,DIP,DEC,SMODIP)
C THIS IS A SPECIAL VERSION OF THE POGO 68/10 MAGNETIC FIELD
C LEGENDRE MODEL. TRANSFORMATION COEFF. G(144) VALID FOR 1973.
C INPUT: DLAT, DLONG=GEOGRAPHIC COORDINATES/DEG.(-90/90,0/360),
C        ALT=ALTITUDE/KM.
C OUTPUT: F TOTAL FIELD (GAUSS), Z DOWNWARD VERTICAL COMPONENT
C        X,Y COMPONENTS IN THE EQUATORIAL PLANE (X TO ZERO LONGITUDE).
C        DIP INCLINATION ANGLE(DEGREE). SMODIP RAWER'S MODFIED DIP.
C SHEIK,1977.
      DIMENSION H(144),XI(3),G(144),FEL1(72),FEL2(72)
      COMMON/CONST/UMR
      DATA FEL1/0.0, 0.1506723,0.0101742, -0.0286519, 0.0092606,
     & -0.0130846, 0.0089594, -0.0136808,-0.0001508, -0.0093977,
     & 0.0130650, 0.0020520, -0.0121956, -0.0023451, -0.0208555,
     & 0.0068416,-0.0142659, -0.0093322, -0.0021364, -0.0078910,
     & 0.0045586,  0.0128904, -0.0002951, -0.0237245,0.0289493,
     & 0.0074605, -0.0105741, -0.0005116, -0.0105732, -0.0058542,
     &0.0033268, 0.0078164,0.0211234, 0.0099309, 0.0362792,
     &-0.0201070,-0.0046350,-0.0058722,0.0011147,-0.0013949,
     & -0.0108838,  0.0322263, -0.0147390,  0.0031247, 0.0111986,
     & -0.0109394,0.0058112,  0.2739046, -0.0155682, -0.0253272,
     &  0.0163782, 0.0205730,  0.0022081, 0.0112749,-0.0098427,
     & 0.0072705, 0.0195189, -0.0081132, -0.0071889, -0.0579970,
     & -0.0856642, 0.1884260,-0.7391512, 0.1210288, -0.0241888,
     & -0.0052464, -0.0096312, -0.0044834, 0.0201764,  0.0258343,
     &0.0083033,  0.0077187/
      DATA FEL2/0.0586055,0.0102236,-0.0396107,
     & -0.0167860, -0.2019911, -0.5810815,0.0379916,  3.7508268,
     & 1.8133030, -0.0564250, -0.0557352, 0.1335347, -0.0142641,
     & -0.1024618,0.0970994, -0.0751830,-0.1274948, 0.0402073,
     &  0.0386290, 0.1883088,  0.1838960, -0.7848989,0.7591817,
     & -0.9302389,-0.8560960, 0.6633250, -4.6363869, -13.2599277,
     & 0.1002136,  0.0855714,-0.0991981, -0.0765378,-0.0455264,
     &  0.1169326, -0.2604067, 0.1800076, -0.2223685, -0.6347679,
     &0.5334222, -0.3459502,-0.1573697,  0.8589464, 1.7815990,
     &-6.3347645, -3.1513653, -9.9927750,13.3327637, -35.4897308,
     &37.3466339, -0.5257398,  0.0571474, -0.5421217,  0.2404770,
     & -0.1747774,-0.3433644, 0.4829708,0.3935944, 0.4885033,
     &  0.8488121, -0.7640999, -1.8884945, 3.2930784,-7.3497229,
     & 0.1672821,-0.2306652, 10.5782146, 12.6031065, 8.6579742,
     & 215.5209961, -27.1419220,22.3405762,1108.6394043/
      K=0
      DO 10 I=1,72
      K=K+1
      G(K)=FEL1(I)
10    G(72+K)=FEL2(I)
      RLAT=DLAT*UMR
      CT=SIN(RLAT)
      ST=COS(RLAT)
      NMAX=11
      D=SQRT(40680925.0-272336.0*CT*CT)
      RLONG=DLONG*UMR
      CP=COS(RLONG)
      SP=SIN(RLONG)
      ZZZ=(ALT+40408589.0/D)*CT/6371.2
      RHO=(ALT+40680925.0/D)*ST/6371.2
      XXX=RHO*CP
      YYY=RHO*SP
      RQ=1.0/(XXX*XXX+YYY*YYY+ZZZ*ZZZ)
      XI(1)=XXX*RQ
      XI(2)=YYY*RQ
      XI(3)=ZZZ*RQ
      IHMAX=NMAX*NMAX+1
      LAST=IHMAX+NMAX+NMAX
      IMAX=NMAX+NMAX-1
      DO 100 I=IHMAX,LAST
100   H(I)=G(I)
      DO 200 K=1,3,2
      I=IMAX
      IH=IHMAX
300   IL=IH-I
      F1=2./(I-K+2.)
      X1=XI(1)*F1
      Y1=XI(2)*F1
      Z1=XI(3)*(F1+F1)
      I=I-2
      IF((I-1).LT.0) GOTO 400
      IF((I-1).EQ.0) GOTO 500
      DO 600 M=3,I,2
      H(IL+M+1)=G(IL+M+1)+Z1*H(IH+M+1)+X1*(H(IH+M+3)-H(IH+M-1))-
     &Y1*(H(IH+M+2)+H(IH+M-2))
      H(IL+M)=G(IL+M)+Z1*H(IH+M)+X1*(H(IH+M+2)-H(IH+M-2))+
     &Y1*(H(IH+M+3)+H(IH+M-1))
600   CONTINUE
500   H(IL+2)=G(IL+2)+Z1*H(IH+2)+X1*H(IH+4)-Y1*(H(IH+3)+H(IH))
      H(IL+1)=G(IL+1)+Z1*H(IH+1)+Y1*H(IH+4)+X1*(H(IH+3)-H(IH))
400   H(IL)=G(IL)+Z1*H(IH)+2.0*(X1*H(IH+1)+Y1*H(IH+2))
700   IH=IL
      IF(I.GE.K) GOTO 300
200   CONTINUE
      S=0.5*H(1)+2.0*(H(2)*XI(3)+H(3)*XI(1)+H(4)*XI(2))
      XT=(RQ+RQ)*SQRT(RQ)
      X=XT*(H(3)-S*XXX)
      Y=XT*(H(4)-S*YYY)
      Z=XT*(H(2)-S*ZZZ)
      F=SQRT(X*X+Y*Y+Z*Z)
      BRH0=Y*SP+X*CP
      Y=Y*CP-X*SP
      X=Z*ST-BRH0*CT
      Z=-Z*CT-BRH0*ST
	zdivf=z/f
        IF(ABS(zdivf).GT.1.) zdivf=SIGN(1.,zdivf)
      DIP=ASIN(zdivf)
	ydivs=y/sqrt(x*x+y*y)
        IF(ABS(ydivs).GT.1.) ydivs=SIGN(1.,ydivs)
      DEC=ASIN(ydivs)
	dipdiv=DIP/SQRT(DIP*DIP+ST)
        IF(ABS(dipdiv).GT.1.) dipdiv=SIGN(1.,dipdiv)
      SMODIP=ASIN(dipdiv)
      DIP=DIP/UMR
      DEC=DEC/UMR
      SMODIP=SMODIP/UMR
      RETURN
      END
C
C
C************************************************************
C*********** INTERPOLATION AND REST ***************************
C**************************************************************
C
C
      SUBROUTINE REGFA1(X11,X22,FX11,FX22,EPS,FW,F,SCHALT,X)
C REGULA-FALSI-PROCEDURE TO FIND X WITH F(X)-FW=0. X1,X2 ARE THE
C STARTING VALUES. THE COMUTATION ENDS WHEN THE X-INTERVAL
C HAS BECOME LESS THAN EPS . IF SIGN(F(X1)-FW)= SIGN(F(X2)-FW)
C THEN SCHALT=.TRUE.
      LOGICAL L1,LINKS,K,SCHALT
      SCHALT=.FALSE.
      EP=EPS
      X1=X11
      X2=X22
      F1=FX11-FW
      F2=FX22-FW
      K=.FALSE.
      NG=2
      LFD=0
      IF(F1*F2.LE.0.0) GOTO 200
   	X=0.0
    	SCHALT=.TRUE.
      	RETURN
200   X=(X1*F2-X2*F1)/(F2-F1)
      GOTO 400
300   	L1=LINKS
	DX=(X2-X1)/NG
     	IF(.NOT.LINKS) DX=DX*(NG-1)
     	X=X1+DX
400   FX=F(X)-FW
      LFD=LFD+1
      IF(LFD.GT.20) THEN
	EP=EP*10.
	LFD=0
      ENDIF
      LINKS=(F1*FX.GT.0.0)
      K=.NOT.K
      IF(LINKS) THEN
	X1=X
 	F1=FX
      ELSE
	X2=X
	F2=FX
      ENDIF
      IF(ABS(X2-X1).LE.EP) GOTO 800
      IF(K) GOTO 300
      IF((LINKS.AND.(.NOT.L1)).OR.(.NOT.LINKS.AND.L1)) NG=2*NG
      GOTO 200
800   RETURN
      END
C
C
      SUBROUTINE TAL(SHABR,SDELTA,SHBR,SDTDH0,AUS6,SPT)
C CALCULATES THE COEFFICIENTS SPT FOR THE POLYNOMIAL
C Y(X)=1+SPT(1)*X**2+SPT(2)*X**3+SPT(3)*X**4+SPT(4)*X**5
C TO FIT THE VALLEY IN Y, REPRESENTED BY:
C Y(X=0)=1, THE X VALUE OF THE DEEPEST VALLEY POINT (SHABR),
C THE PRECENTAGE DEPTH (SDELTA), THE WIDTH (SHBR) AND THE
C DERIVATIVE DY/DX AT THE UPPER VALLEY BOUNDRY (SDTDH0).
C IF THERE IS AN UNWANTED ADDITIONAL EXTREMUM IN THE VALLEY
C REGION, THEN AUS6=.TRUE., ELSE AUS6=.FALSE..
C FOR -SDELTA THE COEFF. ARE CALCULATED FOR THE FUNCTION
C Y(X)=EXP(SPT(1)*X**2+...+SPT(4)*X**5).
      DIMENSION SPT(4)
      LOGICAL AUS6
      Z1=-SDELTA/(100.0*SHABR*SHABR)
      IF(SDELTA.GT.0.) GOTO 500
      SDELTA=-SDELTA
      Z1=ALOG(1.-SDELTA/100.)/(SHABR*SHABR)
500   Z3=SDTDH0/(2.*SHBR)
      Z4=SHABR-SHBR
      SPT(4)=2.0*(Z1*(SHBR-2.0*SHABR)*SHBR+Z3*Z4*SHABR)/
     &  (SHABR*SHBR*Z4*Z4*Z4)
      SPT(3)=Z1*(2.0*SHBR-3.0*SHABR)/(SHABR*Z4*Z4)-
     &  (2.*SHABR+SHBR)*SPT(4)
      SPT(2)=-2.0*Z1/SHABR-2.0*SHABR*SPT(3)-3.0*SHABR*SHABR*SPT(4)
      SPT(1)=Z1-SHABR*(SPT(2)+SHABR*(SPT(3)+SHABR*SPT(4)))
      AUS6=.FALSE.
      B=4.*SPT(3)/(5.*SPT(4))+SHABR
      C=-2.*SPT(1)/(5*SPT(4)*SHABR)
      Z2=B*B/4.-C
      IF(Z2.LT.0.0) GOTO 300
      Z3=SQRT(Z2)
      Z1=B/2.
      Z2=-Z1+Z3
      IF(Z2.GT.0.0.AND.Z2.LT.SHBR) AUS6=.TRUE.
      IF (ABS(Z3).GT.1.E-15) GOTO 400
      Z2=C/Z2
      IF(Z2.GT.0.0.AND.Z2.LT.SHBR) AUS6=.TRUE.
      RETURN
400   Z2=-Z1-Z3
      IF(Z2.GT.0.0.AND.Z2.LT.SHBR) AUS6=.TRUE.
300   RETURN
      END
C
C
C******************************************************************
C********** ZENITH ANGLE, DAY OF YEAR, TIME ***********************
C******************************************************************
C
C
	subroutine soco (ld,t,flat,Elon,
     &		DECLIN, ZENITH, SUNRSE, SUNSET)
c--------------------------------------------------------------------
c	s/r to calculate the solar declination, zenith angle, and
c	sunrise & sunset times  - based on Newbern Smith's algorithm
c	[leo mcnamara, 1-sep-86, last modified 16-jun-87]
c	{dieter bilitza, 30-oct-89, modified for IRI application}
c
c in:	ld	local day of year
c	t	local hour (decimal)
c	flat	northern latitude in degrees
c	elon	east longitude in degrees
c
c out:	declin      declination of the sun in degrees
c	zenith	    zenith angle of the sun in degrees
c	sunrse	    local time of sunrise in hours
c	sunset	    local time of sunset in hours
c-------------------------------------------------------------------
c
	common/const/	dtr
c amplitudes of Fourier coefficients  --  1955 epoch.................
	data  	p1,p2,p3,p4,p6 /
     & 	0.017203534,0.034407068,0.051610602,0.068814136,0.103221204 /
c
c s/r is formulated in terms of WEST longitude.......................
	wlon = 360. - Elon
c
c time of equinox for 1980...........................................
	td = ld + (t + Wlon/15.) / 24.
	te = td + 0.9369
c
c declination of the sun..............................................
	dcl = 23.256 * sin(p1*(te-82.242)) + 0.381 * sin(p2*(te-44.855))
     &      + 0.167 * sin(p3*(te-23.355)) - 0.013 * sin(p4*(te+11.97))
     &      + 0.011 * sin(p6*(te-10.41)) + 0.339137
	DECLIN = dcl
	dc = dcl * dtr
c
c the equation of time................................................
	tf = te - 0.5
	eqt = -7.38*sin(p1*(tf-4.)) - 9.87*sin(p2*(tf+9.))
     &      + 0.27*sin(p3*(tf-53.)) - 0.2*cos(p4*(tf-17.))
	et = eqt * dtr / 4.
c
	fa = flat * dtr
	phi = 0.26179939 * ( t - 12.) + et
c
	a = sin(fa) * sin(dc)
	b = cos(fa) * cos(dc)
	cosx = a + b * cos(phi)
	if(abs(cosx).gt.1.) cosx=sign(1.,cosx)
	zenith = acos(cosx) / dtr
c
c calculate sunrise and sunset times --  at the ground...........
c see Explanatory Supplement to the Ephemeris (1961) pg 401......
c sunrise at height h metres is at...............................
c	chi(h) = 90.83 + 0.0347 * sqrt(h)........................
c this includes corrections for horizontal refraction and........
c semi-diameter of the solar disk................................
	ch = cos(90.83 * dtr)
	cosphi = (ch -a ) / b
c if abs(secphi) > 1., sun does not rise/set.....................
c allow for sun never setting - high latitude summer.............
	secphi = 999999.
	if(cosphi.ne.0.) secphi = 1./cosphi
	sunset = 99.
	sunrse = 99.
	if(secphi.gt.-1.0.and.secphi.le.0.) return
c allow for sun never rising - high latitude winter..............
	sunset = -99.
	sunrse = -99.
	if(secphi.gt.0.0.and.secphi.lt.1.) return
c
	if(cosphi.gt.1.) cosphi=sign(1.,cosphi)
	phi = acos(cosphi)
	et = et / 0.26179939
	phi = phi / 0.26179939
	sunrse = 12. - phi - et
	sunset = 12. + phi - et
	if(sunrse.lt.0.) sunrse = sunrse + 24.
	if(sunset.ge.24.) sunset = sunset - 24.
c
	return
	end
c
C
      FUNCTION HPOL(HOUR,TW,XNW,SA,SU,DSA,DSU)
C-------------------------------------------------------
C PROCEDURE FOR SMOOTH TIME-INTERPOLATION USING EPSTEIN
C STEP FUNCTION AT SUNRISE (SA) AND SUNSET (SU). THE
C STEP-WIDTH FOR SUNRISE IS DSA AND FOR SUNSET DSU.
C TW,NW ARE THE DAY AND NIGHT VALUE OF THE PARAMETER TO
C BE INTERPOLATED. SA AND SU ARE TIME OF SUNRIES AND
C SUNSET IN DECIMAL HOURS.
C BILITZA----------------------------------------- 1979.
	IF(ABS(SU).GT.25.) THEN
		IF(SU.GT.0.0) THEN
			HPOL=TW
		ELSE
			HPOL=XNW
		ENDIF
		RETURN
	ENDIF
      HPOL=XNW+(TW-XNW)*EPST(HOUR,DSA,SA)+
     &	(XNW-TW)*EPST(HOUR,DSU,SU)
      RETURN
      END
C
C
	SUBROUTINE MODA(IN,MONTH,IDAY,IDOY)
C-------------------------------------------------------------------
C CALCULATES DAY OF YEAR (IDOY) FROM MONTH (MONTH) AND DAY (IDAY)
C IF IN=0, OR MONTH (MONTH) AND DAY (IDAY) FROM DAY OF
C YEAR (IDOY), IF IN=1.
C-------------------------------------------------------------------
	DIMENSION	MO(12)
	DATA		MO/0,31,59,90,120,151,181,212,243,273,304,334/
	IMO=0
	MOBE=0
	IF(IN.GT.0) GOTO 5
		IDOY=MO(MONTH)+IDAY
		RETURN
5	IMO=IMO+1
		MOOLD=MOBE
		IF(IMO.GT.12) GOTO 55
		MOBE=MO(IMO)
		IF(MOBE.LT.IDOY) GOTO 5
55		MONTH=IMO-1
		IDAY=IDOY-MOOLD
	RETURN
	END
c
C
	REAL FUNCTION B0POL ( HOUR, SAX, SUX, ISEASON, R, DELA)
C-----------------------------------------------------------------
C Interpolation procedure for bottomside thickness parameter B0.
C Array B0F(ILT,ISEASON,IR,ILATI) distinguishes between day and
C night (ILT=1,2), four seasons (ISEASON=1 spring), low and high
C solar activity (IR=1,2), and low and middle modified dip
C latitudes (ILATI=1,2). In the DATA statement the first value
C corresponds to B0F(1,1,1,1), the second to B0F(2,1,1,1), the
C third to B0F(1,2,1,1) and so on.
C JUNE 1989 --------------------------------------- Dieter Bilitza
C
      REAL		NITVAL
      DIMENSION 	B0F(2,4,2,2),SIPH(2),SIPL(2)
      DATA 	B0F/114.,64.0,134.,77.0,128.,66.0,75.,73.0,
     &		    113.,115.,150.,116.,138.,123.,94.,132.,
     &		    72.0,84.0,83.0,89.0,75.0,85.0,57.,76.0,
     &		    102.,100.,120.,110.,107.,103.,76.,86.0/

  	DO 7033 ISR=1,2
 		DO 7034 ISL=1,2
			DAYVAL   = B0F(1,ISEASON,ISR,ISL)
			NITVAL = B0F(2,ISEASON,ISR,ISL)

C Interpolation day/night with transitions at SAX (sunrise) and SUX (sunset)
7034  			SIPH(ISL) = HPOL(HOUR,DAYVAL,NITVAL,
     &				SAX,SUX,1.,1.)

C Interpolation low/middle modip with transition at 30 degrees modip
7033		SIPL(ISR) = SIPH(1) + (SIPH(2) - SIPH(1)) / DELA

C Interpolation low/high Rz12: linear from 10 to 100
	B0POL=SIPL(1)+(SIPL(2)-SIPL(1))/90.*(R-10.)
	RETURN
	END
c
C
C *********************************************************************
C ************************ EPSTEIN FUNCTIONS **************************
C *********************************************************************
C REF:	H. G. BOOKER, J. ATMOS. TERR. PHYS. 39, 619-623, 1977
C 	K. RAWER, ADV. SPACE RES. 4, #1, 11-15, 1984
C *********************************************************************
C
C
	REAL FUNCTION  RLAY ( X, XM, SC, HX )
C -------------------------------------------------------- RAWER  LAYER
	Y1  = EPTR ( X , SC, HX )
	Y1M = EPTR ( XM, SC, HX )
	Y2M = EPST ( XM, SC, HX )
	RLAY = Y1 - Y1M - ( X - XM ) * Y2M / SC
	RETURN
	END
C
C
	REAL FUNCTION D1LAY ( X, XM, SC, HX )
C ------------------------------------------------------------ dLAY/dX
	D1LAY = ( EPST(X,SC,HX) - EPST(XM,SC,HX) ) /  SC
	RETURN
	END
C
C
	REAL FUNCTION D2LAY ( X, XM, SC, HX )
C ---------------------------------------------------------- d2LAY/dX2
	D2LAY = EPLA(X,SC,HX) /  (SC * SC)
	RETURN
	END
C
C
	REAL FUNCTION EPTR ( X, SC, HX )
C ------------------------------------------------------------ TRANSITION
	COMMON/ARGEXP/ARGMAX
	D1 = ( X - HX ) / SC
	IF (ABS(D1).LT.ARGMAX) GOTO 1
	IF (D1.GT.0.0) THEN
	  EPTR = D1
	ELSE
	  EPTR = 0.0
	ENDIF
	RETURN
1	EPTR = ALOG ( 1. + EXP( D1 ))
	RETURN
	END
C
C
	REAL FUNCTION EPST ( X, SC, HX )
C -------------------------------------------------------------- STEP
	COMMON/ARGEXP/ARGMAX
	D1 = ( X - HX ) / SC
	IF (ABS(D1).LT.ARGMAX) GOTO 1
	IF (D1.GT.0.0) THEN
	  EPST = 1.
	ELSE
	  EPST = 0.
	ENDIF
	RETURN
1	EPST = 1. / ( 1. + EXP( -D1 ))
	RETURN
	END
C
C
	REAL FUNCTION EPSTEP ( Y2, Y1, SC, HX, X)
C---------------------------------------------- STEP FROM Y1 TO Y2
	EPSTEP = Y1 + ( Y2 - Y1 ) * EPST ( X, SC, HX)
	RETURN
	END
C
C
	REAL FUNCTION EPLA ( X, SC, HX )
C ------------------------------------------------------------ PEAK
	COMMON/ARGEXP/ARGMAX
	D1 = ( X - HX ) / SC
	IF (ABS(D1).LT.ARGMAX) GOTO 1
		EPLA = 0
		RETURN
1	D0 = EXP ( D1 )
	D2 = 1. + D0
	EPLA = D0 / ( D2 * D2 )
	RETURN
	END
c
c
	FUNCTION XE2TO5(H,HMF2,NL,HX,SC,AMP)
C----------------------------------------------------------------------
C NORMALIZED ELECTRON DENSITY (N/NMF2) FOR THE MIDDLE IONOSPHERE FROM
C HME TO HMF2 USING LAY-FUNCTIONS.
C----------------------------------------------------------------------
	DIMENSION	HX(NL),SC(NL),AMP(NL)
	SUM = 1.0
	DO 1 I=1,NL
	   YLAY = AMP(I) * RLAY( H, HMF2, SC(I), HX(I) )
	   zlay=10.**ylay
1	   sum=sum*zlay
	XE2TO5 = sum
	RETURN
	END
C
C
	REAL FUNCTION XEN(H,HMF2,XNMF2,HME,NL,HX,SC,AMP)
C----------------------------------------------------------------------
C ELECTRON DENSITY WITH NEW MIDDLE IONOSPHERE
C----------------------------------------------------------------------
	DIMENSION	HX(NL),SC(NL),AMP(NL)
C
	IF(H.LT.HMF2) GOTO 100
		XEN = XE1(H)
		RETURN
100	IF(H.LT.HME) GOTO 200
		XEN = XNMF2 * XE2TO5(H,HMF2,NL,HX,SC,AMP)
		RETURN
200	XEN = XE6(H)
	RETURN
	END
C
C
	SUBROUTINE VALGUL(XHI,HVB,VWU,VWA,VDP)
C ---------------------------------------------------------------------
C   CALCULATES E-F VALLEY PARAMETERS; T.L. GULYAEVA, ADVANCES IN
C   SPACE RESEARCH 7, #6, 39-48, 1987.
C
C	INPUT:	XHI	SOLAR ZENITH ANGLE [DEGREE]
C
C	OUTPUT:	VDP	VALLEY DEPTH  (NVB/NME)
C		VWU	VALLEY WIDTH  [KM]
C		VWA	VALLEY WIDTH  (SMALLER, CORRECTED BY RAWER)
C		HVB	HEIGHT OF VALLEY BASE [KM]
C -----------------------------------------------------------------------
C
	COMMON	/CONST/UMR
C
	CS = 0.1 + COS(UMR*XHI)
	ABC = ABS(CS)
	VDP = 0.45 * CS / (0.1 + ABC ) + 0.55
	ARL = ( 0.1 + ABC + CS ) / ( 0.1 + ABC - CS)
	ZZZ = ALOG( ARL )
	VWU = 45. - 10. * ZZZ
	VWA = 45. -  5. * ZZZ
	HVB = 1000. / ( 7.024 + 0.224 * CS + 0.966 * ABC )
	RETURN
	END
C
C
	SUBROUTINE ROGUL(IDAY,XHI,SX,GRO)
C ---------------------------------------------------------------------
C   CALCULATES RATIO H0.5/HMF2 FOR HALF-DENSITY POINT (NE(H0.5)=0.5*NMF2)
C   T.L. GULYAEVA, ADVANCES IN SPACE RESEARCH 7, #6, 39-48, 1987.
C
C	INPUT:	IDAY	DAY OF YEAR
C		XHI	SOLAR ZENITH ANGLE [DEGREE]
C
C	OUTPUT:	GRO	RATIO OF HALF DENSITY HEIGHT TO F PEAK HEIGHT
C		SX	SMOOTHLY VARYING SEASON PARAMTER (SX=1 FOR
C			DAY=1; SX=3 FOR DAY=180; SX=2 FOR EQUINOX)
C -----------------------------------------------------------------------
C
	SX = 2. - COS ( IDAY * 0.017214206 )
	XS = ( XHI - 20. * SX) / 15.
	GRO = 0.8 - 0.2 / ( 1. + EXP(XS) )
c same as gro=0.6+0.2/(1+exp(-xs))
	RETURN
	END
C
C
	SUBROUTINE LNGLSN ( N, A, B, AUS)
C --------------------------------------------------------------------
C SOLVES QUADRATIC SYSTEM OF LINEAR EQUATIONS:
C
C	INPUT:	N	NUMBER OF EQUATIONS (= NUMBER OF UNKNOWNS)
C		A(N,N)	MATRIX (LEFT SIDE OF SYSTEM OF EQUATIONS)
C		B(N)	VECTOR (RIGHT SIDE OF SYSTEM)
C
C	OUTPUT:	AUS	=.TRUE.	  NO SOLUTION FOUND
C			=.FALSE.  SOLUTION IS IN  A(N,J) FOR J=1,N
C --------------------------------------------------------------------
C
	DIMENSION	A(5,5), B(5), AZV(10)
	LOGICAL		AUS
C
	NN = N - 1
	AUS = .FALSE.
	DO 1 K=1,N-1
		IMAX = K
		L    = K
		IZG  = 0
		AMAX  = ABS( A(K,K) )
110		L = L + 1
		IF (L.GT.N) GOTO 111
		HSP = ABS( A(L,K) )
		IF (HSP.LT.1.E-8) IZG = IZG + 1
		IF (HSP.LE.AMAX) GOTO 110
111		IF (ABS(AMAX).GE.1.E-10) GOTO 133
			AUS = .TRUE.
			RETURN
133		IF (IMAX.EQ.K) GOTO 112
		DO 2 L=K,N
			AZV(L+1)  = A(IMAX,L)
			A(IMAX,L) = A(K,L)
2			A(K,L)    = AZV(L+1)
		AZV(1)  = B(IMAX)
		B(IMAX) = B(K)
		B(K)    = AZV(1)
112		IF (IZG.EQ.(N-K)) GOTO 1
		AMAX = 1. / A(K,K)
		AZV(1) = B(K) * AMAX
		DO 3 M=K+1,N
3			AZV(M+1) = A(K,M) * AMAX
		DO 4 L=K+1,N
			AMAX = A(L,K)
			IF (ABS(AMAX).LT.1.E-8) GOTO 4
			A(L,K) = 0.0
			B(L) = B(L) - AZV(1) * AMAX
			DO 5 M=K+1,N
5				A(L,M) = A(L,M) - AMAX * AZV(M+1)
4		CONTINUE
1	CONTINUE
	DO 6 K=N,1,-1
		AMAX = 0.0
		IF (K.LT.N) THEN
			DO 7 L=K+1,N
7				AMAX = AMAX + A(K,L) * A(N,L)
			ENDIF
		IF (ABS(A(K,K)).LT.1.E-6) THEN
			A(N,K) = 0.0
		ELSE
			A(N,K) = ( B(K) - AMAX ) / A(K,K)
		ENDIF
6	CONTINUE
	RETURN
	END
C
C
	SUBROUTINE LSKNM ( N, M, M0, M1, HM, SC, HX, W, X, Y, VAR, SING)
C --------------------------------------------------------------------
C   DETERMINES LAY-FUNCTIONS AMPLITUDES FOR A NUMBER OF CONSTRAINTS:
C
C	INPUT:	N	NUMBER OF AMPLITUDES ( LAY-FUNCTIONS)
C		M	NUMBER OF CONSTRAINTS
C		M0	NUMBER OF POINT CONSTRAINTS
C		M1	NUMBER OF FIRST DERIVATIVE CONSTRAINTS
C		HM	F PEAK ALTITUDE  [KM]
C		SC(N)	SCALE PARAMETERS FOR LAY-FUNCTIONS  [KM]
C		HX(N)	HEIGHT PARAMETERS FOR LAY-FUNCTIONS  [KM]
C		W(M)	WEIGHT OF CONSTRAINTS
C		X(M)	ALTITUDES FOR CONSTRAINTS  [KM]
C		Y(M)	LOG(DENSITY/NMF2) FOR CONSTRAINTS
C
C	OUTPUT:	VAR(M)	AMPLITUDES
C		SING	=.TRUE.   NO SOLUTION
C ------------------------------------------------------------------------
C
	LOGICAL		SING
	DIMENSION	VAR(N), HX(N), SC(N), W(M), X(M), Y(M),
     &			BLI(5), ALI(5,5), XLI(5,10)
C
	M01=M0+M1
	SCM=0
	DO 1 J=1,5
		BLI(J) = 0.
		DO 1 I=1,5
1			ALI(J,I) = 0.
	DO 2 I=1,N
		DO 3 K=1,M0
3			XLI(I,K) = RLAY( X(K), HM, SC(I), HX(I) )
		DO 4 K=M0+1,M01
4			XLI(I,K) = D1LAY( X(K), HM, SC(I), HX(I) )
		DO 5 K=M01+1,M
5			XLI(I,K) = D2LAY( X(K), HM, SC(I), HX(I) )
2	CONTINUE
		DO 7 J=1,N
		DO 6 K=1,M
			BLI(J) = BLI(J) + W(K) * Y(K) * XLI(J,K)
			DO 6 I=1,N
6				ALI(J,I) = ALI(J,I) + W(K) * XLI(I,K)
     &					* XLI(J,K)
7	CONTINUE
	CALL LNGLSN( N, ALI, BLI, SING )
	IF (.NOT.SING) THEN
		DO 8 I=1,N
8			VAR(I) = ALI(N,I)
		ENDIF
	RETURN
	END
C
C
	SUBROUTINE INILAY(NIGHT,XNMF2,XNMF1,XNME,VNE,HMF2,HMF1,
     &				HME,HV1,HV2,HHALF,HXL,SCL,AMP,IQUAL)
C-------------------------------------------------------------------
C CALCULATES AMPLITUDES FOR LAY FUNCTIONS
C D. BILITZA, DECEMBER 1988
C
C INPUT:	NIGHT	LOGICAL VARIABLE FOR DAY/NIGHT DISTINCTION
C		XNMF2	F2 PEAK ELECTRON DENSITY [M-3]
C		XNMF1	F1 PEAK ELECTRON DENSITY [M-3]
C		XNME	E  PEAK ELECTRON DENSITY [M-3]
C		VNE	ELECTRON DENSITY AT VALLEY BASE [M-3]
C		HMF2	F2 PEAK ALTITUDE [KM]
C		HMF1	F1 PEAK ALTITUDE [KM]
C		HME	E  PEAK ALTITUDE [KM]
C		HV1	ALTITUDE OF VALLEY TOP [KM]
C		HV2	ALTITUDE OF VALLEY BASE [KM]
C		HHALF	ALTITUDE OF HALF-F2-PEAK-DENSITY [KM]
C
C OUTPUT:	HXL(4)	HEIGHT PARAMETERS FOR LAY FUNCTIONS [KM]
C		SCL(4)	SCALE PARAMETERS FOR LAY FUNCTIONS [KM]
C		AMP(4)	AMPLITUDES FOR LAY FUNCTIONS
C		IQUAL	=0 ok, =1 ok using second choice for HXL(1)
C                       =2 NO SOLUTION
C---------------------------------------------------------------
	DIMENSION	XX(8),YY(8),WW(8),AMP(4),HXL(4),SCL(4)
	LOGICAL		SSIN,NIGHT
c
c constants --------------------------------------------------------
		NUMLAY=4
		NC1 = 2
		ALG102=ALOG10(2.)
c
c constraints: xx == height	yy == log(Ne/NmF2)    ww == weights
c -----------------------------------------------------------------
		ALOGF = ALOG10(XNMF2)
		ALOGEF = ALOG10(XNME) - ALOGF
		XHALF=XNMF2/2.
		XX(1) = HHALF
 		XX(2) = HV1
		XX(3) = HV2
		XX(4) = HME
		XX(5) = HME - ( HV2 - HME )
		YY(1) = -ALG102
		YY(2) = ALOGEF
		YY(3) = ALOG10(VNE) - ALOGF
		YY(4) = ALOGEF
		YY(5) = YY(3)
		YY(7) = 0.0
		WW(2) = 1.
		WW(3) = 2.
		WW(4) = 5.
c
c geometric paramters for LAY -------------------------------------
c difference to earlier version:  HXL(3) = HV2 + SCL(3)
c
		SCL0 = 0.7 * ( 0.216 * ( HMF2 - HHALF ) + 56.8 )
		SCL(1) = 0.8 * SCL0
		SCL(2) = 10.
 		SCL(3) = 9.
	  	SCL(4) = 6.
		HXL(3) = HV2
c
C DAY CONDITION--------------------------------------------------
c earlier tested: 	HXL(2) = HMF1 + SCL(2)
c
	    IF(NIGHT) GOTO 7711
		NUMCON = 8
		HXL(1) = 0.9 * HMF2
		  HXL1T  = HHALF
		HXL(2) = HMF1
		HXL(4) = HME - SCL(4)
		XX(6) = HMF1
 		XX(7) = HV2
		XX(8) = HME
		YY(8) = 0.0
		WW(5) = 1.
	   	WW(7) = 50.
	   	WW(8) = 500.
c without F-region ----------------------------------------------
		IF(XNMF1.GT.0) GOTO 100
			HXL(2)=(HMF2+HHALF)/2.
			YY(6) = 0.
			WW(6) = 0.
			WW(1) = 1.
			GOTO 7722
c with F-region --------------------------------------------
100		YY(6) = ALOG10(XNMF1) - ALOGF
		WW(6) = 3.
		IF((XNMF1-XHALF)*(HMF1-HHALF).LT.0.0) THEN
		  WW(1)=0.5
		ELSE
		  ZET = YY(1) - YY(6)
		  WW(1) = EPST( ZET, 0.1, 0.15)
		ENDIF
		IF(HHALF.GT.HMF1) THEN
		  HFFF=HMF1
		  XFFF=XNMF1
		ELSE
		  HFFF=HHALF
		  XFFF=XHALF
		ENDIF
	        GOTO 7722
c
C NIGHT CONDITION---------------------------------------------------
c different HXL,SCL values were tested including:
c	SCL(1) = HMF2 * 0.15 - 27.1	HXL(2) = 200.
c	HXL(2) = HMF1 + SCL(2)		HXL(3) = 140.
c  	SCL(3) = 5.			HXL(4) = HME + SCL(4)
c	HXL(4) = 105.
c
7711		NUMCON = 7
		HXL(1) = HHALF
		  HXL1T  = 0.4 * HMF2 + 30.
		HXL(2) = ( HMF2 + HV1 ) / 2.
		HXL(4) = HME
		XX(6) = HV2
 		XX(7) = HME
		YY(6) = 0.0
		WW(1) = 1.
		WW(3) = 3.
		WW(5) = 0.5
		WW(6) = 50.
 		WW(7) = 500.
		HFFF=HHALF
		XFFF=XHALF
c
C are valley-top and bottomside point compatible ? -------------
C
7722	IF((HV1-HFFF)*(XNME-XFFF).LT.0.0) WW(2)=0.5
	IF(HV1.LE.HV2+5.0) WW(2)=0.5
c
C DETERMINE AMPLITUDES-----------------------------------------
C
	    NC0=NUMCON-NC1
	    IQUAL=0
2299	    CALL LSKNM(NUMLAY,NUMCON,NC0,NC1,HMF2,SCL,HXL,WW,XX,YY,
     &		AMP,SSIN)
	   	IF(IQUAL.gt.0) GOTO 1937
	    IF((ABS(AMP(1)).GT.10.0).OR.(SSIN)) THEN
		IQUAL=1
		HXL(1)=HXL1T
		GOTO 2299
		ENDIF
1937	    IF(SSIN) IQUAL=2
	    RETURN
	    END
c
c
	subroutine ioncom(h,z,f,fs,t,cn)
c---------------------------------------------------------------
c ion composition model
c A.D. Danilov and A.P. Yaichnikov, A New Model of the Ion
c   Composition at 75 to 1000 km for IRI, Adv. Space Res. 5, #7,
c   75-79, 107-108, 1985
c
c 	h	altitude in km
c	z	solar zenith angle in radians
c	f	latitude in radians
c	fs	10.7cm solar radio flux
c	t	season (month)
c	cn(1)   O+  relative density in percent
c	cn(2)   H+  relative density in percent
c	cn(3)   N+  relative density in percent
c	cn(4)   He+ relative density in percent
c	cn(5)   NO+ relative density in percent
c	cn(6)   O2+ relative density in percent
c	cn(7)   cluster ions  relative density in percent
c---------------------------------------------------------------
c
	dimension	cn(7),cm(7),hm(7),alh(7),all(7),beth(7),
     &			betl(7),p(5,6,7),var(6),po(5,6),ph(5,6),
     &			pn(5,6),phe(5,6),pno(5,6),po2(5,6),pcl(5,6)

	data po/4*0.,98.5,4*0.,320.,4*0.,-2.59E-4,2.79E-4,-3.33E-3,
     &		-3.52E-3,-5.16E-3,-2.47E-2,4*0.,-2.5E-6,1.04E-3,
     &		-1.79E-4,-4.29E-5,1.01E-5,-1.27E-3/
	data ph/-4.97E-7,-1.21E-1,-1.31E-1,0.,98.1,355.,-191.,
     &		-127.,0.,2040.,4*0.,-4.79E-6,-2.E-4,5.67E-4,
     &		2.6E-4,0.,-5.08E-3,10*0./
	data pn/7.6E-1,-5.62,-4.99,0.,5.79,83.,-369.,-324.,0.,593.,
     &		4*0.,-6.3E-5,-6.74E-3,-7.93E-3,-4.65E-3,0.,-3.26E-3,
     &		4*0.,-1.17E-5,4.88E-3,-1.31E-3,-7.03E-4,0.,-2.38E-3/
	data phe/-8.95E-1,6.1,5.39,0.,8.01,4*0.,1200.,4*0.,-1.04E-5,
     &		1.9E-3,9.53E-4,1.06E-3,0.,-3.44E-3,10*0./
	data pno/-22.4,17.7,-13.4,-4.88,62.3,32.7,0.,19.8,2.07,115.,
     &		5*0.,3.94E-3,0.,2.48E-3,2.15E-4,6.67E-3,5*0.,
     &		-8.4E-3,0.,-3.64E-3,2.E-3,-2.59E-2/
	data po2/8.,-12.2,9.9,5.8,53.4,-25.2,0.,-28.5,-6.72,120.,
     &		5*0.,-1.4E-2,0.,-9.3E-3,3.3E-3,2.8E-2,5*0.,4.25E-3,
     &		0.,-6.04E-3,3.85E-3,-3.64E-2/
	data pcl/4*0.,100.,4*0.,75.,10*0.,4*0.,-9.04E-3,-7.28E-3,
     &		2*0.,3.46E-3,-2.11E-2/

	DO 8 I=1,5
	DO 8 J=1,6
		p(i,j,1)=po(i,j)
		p(i,j,2)=ph(i,j)
		p(i,j,3)=pn(i,j)
		p(i,j,4)=phe(i,j)
		p(i,j,5)=pno(i,j)
		p(i,j,6)=po2(i,j)
		p(i,j,7)=pcl(i,j)
8	continue

	s=0.
	do 5 i=1,7
	  do 7 j=1,6
		var(j) = p(1,j,i)*cos(z) + p(2,j,i)*cos(f) +
     &  		 p(3,j,i)*cos(0.013*(300.-fs)) +
     &			 p(4,j,i)*cos(0.52*(t-6.)) + p(5,j,i)
7	  continue
	  cm(i)  = var(1)
	  hm(i)  = var(2)
	  all(i) = var(3)
	  betl(i)= var(4)
	  alh(i) = var(5)
	  beth(i)= var(6)
	  hx=h-hm(i)
 	  if(hx) 1,2,3
1		cn(i) = cm(i) * exp( hx * (hx * all(i) + betl(i)) )
		goto 4
2	  	cn(i) = cm(i)
		goto 4
3		cn(i) = cm(i) * exp( hx * (hx * alh(i) + beth(i)) )
4	  continue
	  if(cn(i).LT.0.005*cm(i)) cn(i)=0.
	  if(cn(i).GT.cm(i)) cn(i)=cm(i)
	  s=s+cn(i)
5	continue
	do 6 i=1,7
6		cn(i)=cn(i)/s*100.
	return
	end
C
C
C
C
      SUBROUTINE CIRA86(IDAY,SEC,GLAT,GLONG,STL,F107A,TINF,TLB,SIGMA)
C*******************************************************************
C   Calculates neutral temperature parameters for IRI using the
C   MSIS-86/CIRA 1986 Neutral Thermosphere Model. The subroutines
C   GTS5, GLOBE5 and GLOBL5 developed by A.E. Hedin (2/26/87) were
C   modified for use in IRI --------- D. Bilitza -------- March 1991
C
C     INPUT:
C        IDAY - DAY OF YEAR
C        SEC - UT(SEC)
C        GLAT - GEODETIC LATITUDE(DEG)
C        GLONG - GEODETIC LONGITUDE(DEG)
C        STL - LOCAL APPARENT SOLAR TIME(HRS)
C        F107A - 3 MONTH AVERAGE OF F10.7 FLUX
C
C     OUTPUT:
C        TINF - EXOSPHERIC TEMPERATURE (K)
C        TLB - TEMPERATURE AT LOWER BOUNDARY (K)
C        SIGMA - SHAPE PARAMETER FOR TEMPERATURE PROFILE
C
C **********************************************************************
      DIMENSION         PLG(9,4)
      DATA DR,DR2/1.72142E-2,0.0344284/,HR/.2618/,SR/7.2722E-5/,
     $     XL/1000./,TLL/1000./,DGTR/1.74533E-2/
      save plg, c, s, c2, c4, s2,
     $     stloc, ctloc, s2tloc, c2tloc, s3tloc, c3tloc
C
C CALCULATE LEGENDRE POLYNOMIALS
C
      IF(XL.EQ.GLAT)   GO TO 15
      C = SIN(GLAT*DGTR)
      S = COS(GLAT*DGTR)
      C2 = C*C
      C4 = C2*C2
      S2 = S*S
      PLG(2,1) = C
      PLG(3,1) = 0.5*(3.*C2 -1.)
      PLG(4,1) = 0.5*(5.*C*C2-3.*C)
      PLG(5,1) = (35.*C4 - 30.*C2 + 3.)/8.
      PLG(6,1) = (63.*C2*C2*C - 70.*C2*C + 15.*C)/8.
      PLG(2,2) = S
      PLG(3,2) = 3.*C*S
      PLG(4,2) = 1.5*(5.*C2-1.)*S
      PLG(5,2) = 2.5*(7.*C2*C-3.*C)*S
      PLG(6,2) = 1.875*(21.*C4 - 14.*C2 +1.)*S
      PLG(7,2) = (11.*C*PLG(6,2)-6.*PLG(5,2))/5.
      PLG(3,3) = 3.*S2
      PLG(4,3) = 15.*S2*C
      PLG(5,3) = 7.5*(7.*C2 -1.)*S2
      PLG(6,3) = 3.*C*PLG(5,3)-2.*PLG(4,3)
      PLG(4,4) = 15.*S2*S
      PLG(5,4) = 105.*S2*S*C
      PLG(6,4)=(9.*C*PLG(5,4)-7.*PLG(4,4))/2.
      PLG(7,4)=(11.*C*PLG(6,4)-8.*PLG(5,4))/3.
      XL=GLAT
   15 CONTINUE
      IF(TLL.EQ.STL)   GO TO 16
      STLOC = SIN(HR*STL)
      CTLOC = COS(HR*STL)
      S2TLOC = SIN(2.*HR*STL)
      C2TLOC = COS(2.*HR*STL)
      S3TLOC = SIN(3.*HR*STL)
      C3TLOC = COS(3.*HR*STL)
      TLL = STL
   16 CONTINUE
C
      DFA=F107A-150.
C
C EXOSPHERIC TEMPERATURE
C
C         F10.7 EFFECT
      T1 =  ( 3.11701E-3 - 0.64111E-5 * DFA ) * DFA
        F1 = 1. + 0.426385E-2 * DFA
        F2 = 1. + 0.511819E-2 * DFA
        F3 = 1. + 0.292246E-2 * DFA
C        TIME INDEPENDENT
      T2 = 0.385528E-1 * PLG(3,1) + 0.303445E-2 * PLG(5,1)
C        SYMMETRICAL ANNUAL AND SEMIANNUAL
	CD14 = COS( DR  * (IDAY+8.45398) )
      	CD18 = COS( DR2 * (IDAY-125.818) )
      	CD32 = COS( DR  * (IDAY-30.0150) )
      	CD39 = COS( DR2 * (IDAY-2.75905) )
      T3 = 0.805486E-2 * CD32 + 0.14237E-1 * CD18
C        ASYMMETRICAL ANNUAL AND SEMIANNUAL
      T5 = F1 * (-0.127371 * PLG(2,1) - 0.302449E-1 * PLG(4,1) ) * CD14
     &       - 0.192645E-1 * PLG(2,1) * CD39
C        DIURNAL
      T71 =  0.123512E-1 * PLG(3,2) * CD14
      T72 = -0.526277E-2 * PLG(3,2) * CD14
      T7 = ( -0.105531 *PLG(2,2) - 0.607134E-2 *PLG(4,2) + T71 ) *CTLOC
     4   + ( -0.115622 *PLG(2,2) + 0.202240E-2 *PLG(4,2) + T72 ) *STLOC
C        SEMIDIURNAL
      T81 = 0.386578E-2 * PLG(4,3) * CD14
      T82 = 0.389146E-2 * PLG(4,3) * CD14
      T8= (-0.516278E-3 *PLG(3,3) - 0.117388E-2 *PLG(5,3) +T81)*C2TLOC
     3   +( 0.990156E-2 *PLG(3,3) - 0.354589E-3 *PLG(5,3) +T82)*S2TLOC
C        TERDIURNAL
      Z1 =  PLG(5,4) * CD14
      Z2 =  PLG(7,4) * CD14
      T14=(0.147284E-2*PLG(4,4)-0.173933E-3*Z1+0.365016E-4*Z2)*S3TLOC
     2   +(0.341345E-3*PLG(4,4)-0.153218E-3*Z1+0.115102E-3*Z2)*C3TLOC
      T7814 = F2 * ( T7 + T8 + T14 )
C        LONGITUDINAL
      T11= F3 * (( 0.562606E-2 * PLG(3,2) + 0.594053E-2 * PLG(5,2) +
     $       0.109358E-2 * PLG(7,2) - 0.301801E-2 * PLG(2,2) -
     $       0.423564E-2 * PLG(4,2) - 0.248289E-2 * PLG(6,2) +
     $      (0.189689E-2 * PLG(2,2) + 0.415654E-2 * PLG(4,2)) * CD14
     $     ) * COS(DGTR*GLONG) +
     $     ( -0.11654E-1 * PLG(3,2) - 0.449173E-2 * PLG(5,2) -
     $       0.353189E-3 * PLG(7,2) + 0.919286E-3 * PLG(2,2) +
     $       0.216372E-2 * PLG(4,2) + 0.863968E-3 * PLG(6,2) +
     $      (0.118068E-1 * PLG(2,2) + 0.331190E-2 * PLG(4,2)) * CD14
     $     ) * SIN(DGTR*GLONG) )
C        UT AND MIXED UT,LONGITUDE
      T12 = ( 1. - 0.565411 * PLG(2,1) ) * COS( SR*(SEC-31137.0) ) *
     $ (-0.13341E-1*PLG(2,1)-0.243409E-1*PLG(4,1)-0.135688E-1*PLG(6,1))
     $ +    ( 0.845583E-3 * PLG(4,3) + 0.538706E-3  * PLG(6,3) ) *
     $      COS( SR * (SEC-247.956) + 2.*DGTR*GLONG )
C  Exospheric temperature TINF/K  [Eq. A7]
      TINF = 1041.3 * ( 1. + T1+T2+T3+T5+T7814+T11+T12 ) * 0.99604
C
C TEMPERATURE DERIVATIVE AT LOWER BOUNDARY
C
C         F10.7 EFFECT
      T1 =  0.252317E-2 * DFA
C        TIME INDEPENDENT
      T2 = -0.467542E-1 * PLG(3,1) + 0.12026 * PLG(5,1)
C        ASYMMETRICAL ANNUAL
      	CD14 = COS( DR  * (IDAY+8.45398) )
      T5 = -0.13324 * PLG(2,1)  * CD14
C        SEMIDIURNAL
      ZZ = PLG(4,3) * CD14
      T81 = -0.973404E-2 * ZZ
      T82 = -0.718482E-3 * ZZ
      T8 =(0.191357E-1 *PLG(3,3) + 0.787683E-2 *PLG(5,3) + T81) *C2TLOC
     3  + (0.125429E-2 *PLG(3,3) - 0.233698E-2 *PLG(5,3) + T82) *S2TLOC
C  dTn/dh at lower boundary  [Eq. A6]
      G0 = 0.166728E2 * ( 1. + T1+T2+T5+T8 ) * 0.951363
C
C NEUTRAL TEMPERATURE AT LOWER BOUNDARY 120KM
C
      	CD9  = COS( DR2 * (IDAY-89.3820) )
      	CD11 = COS( DR  * (IDAY+8.45398) )
      T1 = 0.568478E-3 * DFA
      T4 = 0.107674E-1 * CD9
      T5 =-0.192414E-1 * PLG(2,1) * CD11
      T7 = -0.2002E-1 *PLG(2,2) *CTLOC - 0.195833E-2 *PLG(2,2) *STLOC
      T8 = (-0.938391E-2 * PLG(3,3) - 0.260147E-2 * PLG(5,3)
     $       + 0.511651E-4 * PLG(6,3) * CD11 ) * C2TLOC
     $   + ( 0.131480E-1 * PLG(3,3) - 0.808556E-3 * PLG(5,3)
     $       + 0.255717E-2 * PLG(6,3) * CD11 ) * S2TLOC
C  Tn at lower boundary 120km   [Eq. A8]
      TLB = 386.0 * ( 1. + T1+T4+T5+T7+T8 ) * 0.976619
C  Sigma      [Eq. A5]
      SIGMA = G0 / ( TINF - TLB )
      RETURN
      END
C
C
      FUNCTION TN(H,TINF,TLBD,S)
C--------------------------------------------------------------------
C       Calculate Temperature for MSIS/CIRA-86 model
C--------------------------------------------------------------------
      ZG2 = ( H - 120. ) * 6476.77 / ( 6356.77 + H )
      TN = TINF - TLBD * EXP ( - S * ZG2 )
      RETURN
      END
C
C
      FUNCTION DTNDH(H,TINF,TLBD,S)
C---------------------------------------------------------------------
      ZG1 = 6356.77 + H
      ZG2 = 6476.77 / ZG1
      ZG3 = ( H - 120. ) * ZG2
      DTNDH = - TLBD * EXP ( - S * ZG3 ) * ( S / ZG1 * ( ZG3 - ZG2 ) )
      RETURN
      END