subroutine subhmf2c(alati,along,idy,xf2,xm3,xhm)
C................................................................Dec. 2023
C
C................................................................Oct.2013
C Use subinter1.for
C ...............................................................Jan. 2012
C
C Produce hmF2 from xf2, xm3 for one day:
C xf2(0:23),xm3(0:23),xhf(0:23) ! daily sets for UT=0,...,23
C
c "TAB-hmF2 FROM foF2,foE,M3F2 FOR MONTH (31,24)
c
c THE MEAN COV-INDEX (3 SOLAR ROT.) IS EXPECTED
c DEN IS THE ELECTRON DENSITY IN M-3 H=HEIGHT, KM
c  LATI,LONG,DIP,NMF2,MAGBR,MLAT,MLONG,FOF2,HMF2,M3F2,FOE,RZ,COV
      DIMENSION HF2(31,0:23),IM(0:12)
	+,GLATS(59),GLONS(59),JFS(59),JHS(59),JTS(59)
c-	+,rzar(3),arig(3)
	+,xf2(0:23),xm3(0:23),xhm(0:23) ! daily sets for UT=0,...,23
	 INTEGER*2 ID1,ID2
	+,iyri,idyi
	CHARACTER*2	sta,ayr,amn,dy1,dy2,DD1,DD2,YY
	+,STS(59),AMNS(59),DY22,ADY,AUT
	CHARACTER*4	YEAR,TT
	CHARACTER*13 STATS(59),STAT
	CHARACTER*1 hhh
		CHARACTER(10) DD
      CHARACTER(5) ZZ
	integer*2 iyr,imn,idy,jut
	parameter (pi12=0.26179939)
	parameter (DR=1.745329252E-2,RD=57.29577951)
	parameter (xhi0=86.232928)
	COMMON   /CONST/UMR
      COMMON /BL1/STA,YEAR,YY,AMN,DY1,DY2,DD1,DD2,ID1,ID2,JF,JH,JT,DY22
	COMMON /BL2/ STAT,GLAT,GLON,DD,TT,ZZ,R12
	COMMON /BL3/ STATS,STS,GLATS,GLONS,JFS,JHS,II,AMNS,JTS
      DATA IM/31,31,28,31,30,31,30,31,31,30,31,30,31/

	      UMR=.0174533
      PHI=3.141593
	call blet2(idy,ADY)
C
	do i=1,31
	do k=0,23
	hf2(i,k)=0.
	enddo
	enddo
c Test
	read(year,*) ryear
	iyear=int(ryear)
      ayr=year(3:4)
	read(ayr,*) yr
	iyr=int(yr)
	read(amn,*) rmn
	mn=int(rmn)

	IF(mn.EQ.0) GOTO 2

c 1110	write (*,*) 'Enter DY1,DY2:'
c	read (*,1109) DY1,DY2
c 1109 format(A2)     
c==      read(dy1,*) rdy1
c==	read(dy2,*) rdy2
c==	idy1=int(rdy1)
c==	idy2=int(rdy2)

c==	     outfile1='c:\web\idce\styrh.br'
c==	outfile1(13:14)=sta
C ===========================================
C++ 
C Add extract Rzs from ig_rz.dat file
	imn=mn
	iyri=iyr
	idyi=15
	ldai=ndoy(iyri,imn,idyi)
CREM	call tcon(iyear,mn,15,ldai,rzar,arig,ttt,nmonth)
c-	rz=rzar(3)        
c- 	COV=63.75+RZ*(0.728+RZ*0.00089) ! COV index F10.7
	 cov=63.7+(0.728+8.9E-4*R12)*R12
		if (cov.gt.193.) cov=193.
	      if (cov.lt.63.0) cov=63.0

	IF(IYR.LT.90) THEN
         IYYYY=2000+IYR
      ELSE
         IYYYY=1900+IYR
      ENDIF
      IF (YR/4.eq.INT(YR/4)) THEN
 	 IM(2)=29 
	      IDNR=366
	 ELSE 
	 IM(2)=28
	      IDNR=365
	endif

	ABSLAT=ABS(ALATI)
C
crem	numd=idy2-idy1+1       ! number of days for hmF2 calculations
		numd=1


C ++++++++++++++++++ Input of foF2 cloned data and M3000 data:++++++++++++++++++++++
	ihflag=0
	goto 119
  118	continue  ! no data h-input
	goto 117
C
  119	continue 

  117	continue
	jd=idy ! one given day
CC CALCULATION OF DAY OF YEAR	(LDA)
	      mosum=0
      if(mn.gt.1) then
         do 1234 i=1,mn-1
1234   mosum=mosum+im(i)
         endif
      lda=mosum+jd          ! daynr

C !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 	do 1111 jut=0,23               ! start cycle on UT 
	fof2=xf2(jut)
	xm3000=xm3(jut)
	if (xhm(jut).gt.0.) goto 1111
	if ((xm3000.lt.2.).or.(fof2.eq.0.)) then
	hhh='h'
	call blet2(jut,aut)
CREM	call subinterj(hhh,ayr,amn,ady,aut,resm)
	call subinter1(hhh,ayr,amn,ady,aut,resm)
	xhm(jut)=resm
	goto 1111
	endif
	ut1=float(iut)
      hour=mod(ut1+along/15.0+24.0,24.0)	 ! Local time
	   call sdec(mtn,ut1,sdelta,cdelta)
      cchi=sin(alati*DR)*sdelta-cos(alati*DR)*cdelta*
     * cos(pi12*hour)
      xhi=atan2(sqrt(1.0-cchi*cchi),cchi)*RD
	cchinon=sin(alati*DR)*sdelta-cos(alati*DR)*cdelta*
     * cos(pi12*12.0)
c-      xhinon=djoin(90.0-0.24*fexp(20.0-0.20*xhi),xhi,12.0,
c-     & xhi-xhi0)
	xhinon=atan2(sqrt(1.0-cchinon*cchinon),cchinon)*RD

	FOE=FOEEDI(COV,XHI,XHINON,ABSLAT)
c
		fof2=xf2(jut)
		xm3000=xm3(jut)
C++
	HMF2=peakh(FOE,FOF2,XM3000)
 1112	if (hmf2.lt.170.) then
 	hmf2=hmf2+10.
	goto 1112
           endif
	if (hmf2.gt.600.) then
	hmf2=0.
	endif
c 350	hmf(jj,jut)=hmf2
 350	xhm(jut)=hmf2
	if (xhm(jut).gt.900.) then
Crem	call subinterj(hhh,ayr,amn,ady,aut,resm)
	call subinter1(hhh,ayr,amn,ady,aut,resm)
	xhm(jut)=resm
	endif
c==	ihmf(jj,jut)=nint(hmf2)
 1111	continue          ! end UT-cycle


    2	continue
c==      write(*,*) sta
C      pause ' '
      
      return
      END
C--------------------------------------------------
      real function peakh(foE,foF2,M3000)
      real MF,M3000
      sqM=M3000*M3000
      MF=M3000*sqrt((0.0196*sqM+1.)/(1.2967*sqM-1.0))
      If(foE.ge.1.0E-30) then
         ratio=foF2/foE
         ratio=djoin(ratio,1.75,20.0,ratio-1.75)
         dM=0.253/(ratio-1.215)-0.012
      else
         dM=-0.012
      endif
      peakh=1490.0*MF/(M3000+dM)-176.0
      return
      end
C--------------------------------------------------
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.0)
      return
      end
C-----------------------------------------------------------------
       real function fexp(a)
      real a
      if(a.gt.80.0) then
         fexp=5.5406E34
         return
      endif
      if(a.lt.-80.0D0) then
         fexp=1.8049E-35
         return
      endif
      fexp=exp(a)
      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 ----------------------------------------------------
      subroutine sdec(mth,UT,sdelta,cdelta)
      parameter (DR=1.745329252E-2)

      doy=mth*30.5-15.0
      t =doy + (18.0-UT)/24.0
      amrad=(0.9856*t - 3.289)*DR
      aLrad = amrad + (1.916*sin(amrad)+0.02*sin(2.0*amrad)+
     + 282.634)*DR
      sdelta=0.397820*sin(aLrad)
      cdelta=sqrt(1.0-sdelta*sdelta)
      return
      end
C ----------------------------------------------------