C+++++++++++++++++++++++++++++++++++++++++++++++++++++
      subroutine subxhi(xlati,xlongi,nrday,ut,XHI,XHINON)
C
      PHI=3.1415927
       PI=ATAN(1.0)*4.
        UMR=PI/180.
C-      UT=HOUR-XLONGI/15.
	HOUR=UT+XLONGI/15.
	IF (hour.ge.24.) hour=hour-24.
C
      SUNDEC=-0.40915*COS(2.0*PHI/365.25*(NRDAY+8))
      Z1=SIN(SUNDEC)*SIN(XLATI*UMR)
      Z2=COS(SUNDEC)*COS(XLATI*UMR)
      IF (ABS(Z2).GT.0.0) GOTO 120
      SAX=24.0
      IF (Z1.GE.0.0) SAX=0.0
      GOTO 140
  120 IF (ABS(Z1/Z2).LE.1.0) GOTO 510
      SAX=24.0
      IF ((Z1+Z2).GE.0.0) SAX=0.0
      GOTO 140
  510 SAX=12.0-ACOS(-Z1/Z2)/(UMR*15.0)
  140 SUX=24.0-SAX
      XLSTA=15.0*(HOUR-12.0)
      COSXHI=Z1+Z2*COS(XLSTA*UMR)
      XHI=ACOS(COSXHI)/UMR
      XLSTA12=0.0
      COSXHI12=Z1+Z2*COS(XLSTA12*UMR) ! remove?
      COSXHI12=Z1+Z2
      XHINON=ACOS(COSXHI12)/UMR
	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.
C+      COMMON/CONST/UMR,PI
C variation with solar activity (factor A) ...............
       PI=ATAN(1.0)*4.
        UMR=PI/180.
C+

      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---------------------------------------------------------------------
      real function peakh(foE,foF2,M3000)
C     alat:  gg. latitude  (degrees N)
C     along: gg. longitude (degrees E)
C     mth:   month (1 .. 12)
C     flx:   10.7 cm solar radio flux (flux units)
      implicit real (a-h,o-z)
      real MF,M3000
      sqM=M3000*M3000
      MF=M3000*sqrt((0.0196E0*sqM+1.E0)/(1.2967E0*sqM-1.0E0))
      If(foE.ge.1.0E-30) then
         ratio=foF2/foE
         ratio=djoin(ratio,1.75E0,20.0E0,ratio-1.75E0)
         dM=0.253E0/(ratio-1.215E0)-0.012E0
      else
         dM=-0.012E0
      endif
      peakh=1490.0E0*MF/(M3000+dM)-176.0E0
      return
      end
C-------------------------------------------
      real function djoin(f1,f2,alpha,x)
      real f1,f2,alpha,x,ee,fexp
      ee=fexp(alpha*x)
      djoin=(f1*ee+f2)/(ee+1.0E0)
      return
      end
C------------------------------------------------
      real function fexp(a)
      real a
      if(a.gt.80.0E0) then
         fexp=5.5406E34
         return
      endif
      if(a.lt.-80.0E0) then
         fexp=1.8049E-35
         return
      endif
      fexp=exp(a)
      return
      end
C-------------------------------------------------