C.......................................................May 2014
C Plas-cor using Geomagnetic Local Time, GMLT
C					                                  Feb. 2014
C
C       Count global CCIR-based median afc(UT) and ahc*UT) from CCIR maps fcF2 and hcF2(from M3000 map) 
C       for given UT 
      subroutine subsmi4(AYEAR,AMN,ADY,IMEDF,IMEDH,IRIT)
C
C
      DIMENSION im(12),av(0:72)
	+,IMEDF(71,0:72,0:23),IMEDH(71,0:72,0:23) ! results of median maps 
	+,resmodip(71,0:72) ! modip data for IONEX map
	+,IRIT(71,0:72,0:23) ! TEC-NeQ-P
      dimension R12y(0:12),F2(13,76,2),FM3(9,49,2)
	+,F2n(13,76,2),FM3n(9,49,2) !current and prev/next month
	integer month
      character*2 ayr,amn,ady
	character*4 AYEAR
C++	integer*2 iyr,imn,idy,iut
	DATA  IM/31,28,31,30,31,30,31,31,30,31,30,31/
      common/BLOKM/resmodip
	COMMON /bplas/r12i,glat,ylon,ut,ldaynr,heitop,eltop
C
c
      PHI=ATAN(1.0)*4.
	  UMR=PHI/180.

      XFOF2=0.
      XHMF2=0.
	XHI=0.0
c
C ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c
	AYR=AYEAR(3:4)
		read(ayear,*) xyear
	iyear=int(xyear)
		read(ayr,*) yr
		read(amn,*) rmn
	imn=int(rmn)
	month=int(rmn)
C
 	iyr=int(yr)
	if(int(yr/4.)*4.eq.iyr) THEN 
	im(2)=29
	                          ELSE
	IM(2)=28
	ENDIF
	ndnr=im(imn)
C
		read(ady,*) rdy
	idy=int(rdy)
	ldaynr=ndoy(iyr,imn,idy)
C 
	mnmidy=im(imn)/2
	if (idy.lt.mnmidy+1) then
	monthn=month-1
	if (monthn.eq.0) monthn=12
	else
	monthn=month+1
	if (monthn.gt.12) monthn=1
	endif
C Read R12 for current year:
	  open(19,file='R12.dat',status='OLD')
  211    read(19,*, end=99) j,(R12y(k),k=1,12)
         if (j.lt.iyear) then
	R12y(0)=R12y(12)
	   goto 211
	endif
         R12=R12y(month)
	if ((monthn.eq.12).and.(month.eq.1)) then
	 R12n=R12y(0)
	else
	 R12n=R12y(monthn)
	endif
         if (R12.gt.150.0) R12=150.0
	         if (R12n.gt.150.0) R12n=150.0
   99 continue
      close(unit=19)

		dyri=idy-15.15
	R12i=R12+(R12n-R12)*dyri/30.3 ! SSN for current day

         cov=63.7+(0.728+8.9E-4*R12)*R12
	if (cov.gt.193.) cov=193.
	      if (cov.lt.63.0) then
	 cov=63.0
      endif

         covn=63.7+(0.728+8.9E-4*R12n)*R12n
	if (covn.gt.193.) covn=193.
	      if (covn.lt.63.0) then
	 covn=63.0
      endif
C
	call subccir(month,monthn,F2,F2n,FM3,FM3n) !Read CCIR coefficients
C
	         call submodip  ! read modip for IONEX grids
C UT cycle
	DO 203 iut=0,23
	UT=float(iut)
C
C GLAT cycle
C
	glat=-90.
	do 301 lat=1,71
	glat=glat+2.5
C
	do k=0,72
	av(k)=0.        !
	enddo

C GLONG cycle
	glon=-185.
	  	do lon=0,72
	glon=glon+5.0
	if (glon.lt.0.) then
	ylon=glon+360.
	else
	ylon=glon
	endif
	idmn=im(month)
	smodip=resmodip(lat,lon)

C+ Geomagnetic local time
	UTS=ut*3600.        ! UT /sec.
	xday=float(ldaynr)
	gmlt=fmlt(uts,glat,ylon,xday)

C
	call prepccir(glat,ylon,smodip,month,monthn,idy,idmn,UT,R12,R12n
	+,cov,covn,F2,F2n,FM3,FM3n,hm,fm,foe,foF1,teci)

C
	IMEDF(lat,lon,iut)=nint(fm*10.) ! results of CCIR foF2 map
	IMEDH(lat,lon,iut)=nint(hm) ! results of CCIR hmF2 map
	IRIT(lat,lon,iut)=nint(teci) ! results of IRI-Plas CCIR TEC map TECU*10.
	    enddo          !  
C
c	   write(*,109) (iav(k),k=0,72)

  301	continue                 ! end lat cycle of median map
	write(*,*) ' UT= ',iut
  203 continue   ! end ut cycle
 
  108 format(3A2,24(1X,I4))	! 
  109 format(73(1X,I4))	! 
      
C
       GOTO 1234
    2 CONTINUE
 1234	continue  !

c//      write(*,*) infile,outfile
c//      pause ' ' 
   34 RETURN
       END
C==============================================================================
C
      SUBROUTINE SUBCCIR(IMONTH,IMONTHN,F2,F2n,FM3,FM3n)
C     Calculation of coefficients
C     FOF2(76,13) ; SM3000(49,9) 
C     for level of solar activity, MODIP and current month + prec/next month
C     using CCIR maps
      DIMENSION F2(13,76,2),FM3(9,49,2),F2n(13,76,2),FM3n(9,49,2)
C+ SMI input of CCIR coefficients
      OPEN (UNIT=16,FILE='BIN16',STATUS='OLD',
     *    ACCESS='DIRECT',FORM='UNFORMATTED',RECL=11432)

      READ (16,rec=imonth) F2,FM3
      REWIND 16
      READ (16,rec=imonthn) F2n,FM3n
      REWIND 16
	      CLOSE (UNIT=16)
      RETURN
      END
C
      subroutine cciri(xMODIP,mth,UT,R12,alat,along,foF2,M3000,F2x,FM3x)
      implicit real (a-h,o-z)
      real M3000
      dimension FF0(988),xm0(441)
	+,F2x(13,76,2)
	+,FM3x(9,49,2)
      integer QM(7),QF(9)
      save
      data QF/11,11,8,4,1,0,0,0,0/,QM/6,7,5,2,1,0,0/
      data montha,monthb,Rga/13,14,-10.0D0/
      if (mth.ne.montha) then
         montha=mth
      endif
      if (R12.ne.Rga.or.mth.ne.monthb) then
         RR2=R12/100.0
         RR1=1.0-RR2
         do i=1,76
         do j=1,13
            k=j+13*(i-1)
            FF0(k)=F2x(j,i,1)*RR1+F2x(j,i,2)*RR2
         enddo
         enddo
         do i=1,49
            do j=1,9
            k=j+9*(i-1)
            xm0(k)=FM3x(j,i,1)*RR1+FM3x(j,i,2)*RR2
		  enddo
         enddo
         Rga=R12
         monthb=mth
      endif
      foF2= gamma1(xMODIP,alat,along,UT,6,QF,9,76,13,988,FF0)
	M3000=gamma1(xMODIP,alat,along,UT,4,QM,7,49, 9,441,xm0)
      return
      end
C--------------------------------------------------

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--------------------------------------------------

      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************************************************************
      subroutine prepccir(alat,along,xmodip,mth,mth1,idy,idmn,UT,R12
	+,R12n,flx,flxn,F2,F2n,FM3,FM3n,hmf2,fof2,foe,foF1,teci)
      real M3000,M3000a,M3000b
C++	integer*2 idy
	dimension hm(3),aep(3),bb(6)
      dimension 
	+ F2(13,76,2)
	+,FM3(9,49,2)
	+,F2n(13,76,2)
	+,FM3n(9,49,2)
      save
      parameter (pi12=0.26179939)
      parameter (DR=1.745329252E-2,RD=57.29577951)
	COMMON /bplas/rzs,plat,plong,put,ldaynr,heitop,eltop 
      data UT0,mth0/-100.0,-1/
      ut1=mod(UT+24.0,24.0)

      
      if (ut1.ne.UT0.or.mth1.ne.mth0) then
         call sdec(mth1,ut1,sdelta,cdelta)
         UT0=ut1
         mth0=mth1
      endif
	if (idy.lt.idmn/2+1) then
C 30.3 days between mid-prec month and mid-current month 
	dyri=15.15+idy
	else
C 30.3 days between mid-current month and mid-next month
	dyri=idy-15.15
	endif
      xlt=mod(ut1+along/15.0+24.0,24.0)
      cchi=sin(alat*DR)*sdelta-cos(alat*DR)*cdelta*
     * cos(pi12*xlt)
	cchinon=sin(alat*DR)*sdelta-cos(alat*DR)*cdelta*
     * cos(pi12*12.0) 
      chi=atan2(sqrt(1.0-cchi*cchi),cchi)*RD
	chinon=atan2(sqrt(1.0-cchi*cchinon),cchinon)*RD
      call cciri(xMODIP,mth,UT0,R12,alat,along,foF2a,M3000a,F2,FM3)
      call cciri(xMODIP,mth1,UT0,R12n,alat,along,foF2b,M3000b,F2n,FM3n)
	fof2=foF2a+(foF2b-foF2a)*dyri/30.3
	abslat=abs(alat)
	foEa=FOEEDI(flx,chi,chinon,abslat)
	foEb=FOEEDI(flxn,chi,chinon,abslat)
	foe=foea+(foeb-foea)*dyri/30.3
	 efoF1a=ef1r1(foF2a,foEa)
	 efoF1b=ef1r1(foF2b,foEb)
	foF1=efoF1a+(efoF1b-efoF1a)*dyri/30.3
      M3000=M3000a+(M3000b-M3000a)*dyri/30.3
	call prepmdgr(R12,foF2,foF1,foE,M3000,hm,bb,aep)
	hmF2=hm(1)
C
	call prepPlas(hm,aep,bb)
C
	call subinteg(20200.0,teci,hm,aep,bb) ! vTEC [80:20000]km
      return
      end
C-----------------------------------------------------------------
      subroutine submodip
C TLG.............................................March 2014
C Read data from file 'modip2.asc' 
      dimension resmodip(71,0:72)
	character*20 filenam
	common/BLOKM/resmodip
c
	filenam='modip2.asc'
      open(21,file=filenam,form='FORMATTED')
      do i=0,72
      read(21,*) (resmodip(j,i),j=1,71)
      enddo
      close(unit=21)
	return
	end

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 ----------------------------------------------------
	function Fne(x)
c-      FNe(X)=0.124*X*X
	FNe=0.124*X*X
	return
	end
C ----------------------------------------------------
	function FEpst(X,Y,Z,W)
c-      FEpst(X,Y,Z,W)=X*fexp((W-Y)/Z)/(1.0+fexp((W-Y)/Z))**2
	FEpst=X*fexp((W-Y)/Z)/(1.0+fexp((W-Y)/Z))**2
	return
	end

C ----------------------------------------------------
      subroutine prepmdgr(R12,foF2,efoF1,efoE,M3000,hm,bb,aep)
      real NmE,NmF1,NmF2,M3000
      dimension aep(3),hm(3),bb(6)
      data hmE /120.0/

      NmF2=FNe(foF2)
      NmF1=FNe(efoF1)
      NmE= FNe(efoE)
      hmF2=peakh(efoE,foF2,M3000)
      hmF1=0.5*(hmF2+hmE)
      hm(1)=hmF2
      hm(2)=hmF1
      hm(3)=hmE

      dNdHmx=-3.467+1.714*log(foF2)+2.02*log(M3000)
      dNdHmx=exp(dNdHmx)
      B2bot=0.385E2*NmF2/dNdHmx
      hF2F1E=hmF2-hmF1
      B1top=0.3*hF2F1E
      B1bot=0.5*hF2F1E
      Betop=B1bot
      if (Betop.lt.7.0) Betop=7.0
      Bebot=5.0
      aep(1)=4.0*NmF2
      aep(2)=4.0*NmF1
      aep(3)=4.0*NmE
      if (efoF1.le.0.5) then
         aep(2)=0.0
         aep(3)=4.0*(NmE-FEpst(aep(1),hmF2,B2bot,hmE))
      else
         do izml=1,5
            aep(2)=4.0*(NmF1-FEpst(aep(1),hmF2,B2bot,hmF1)
     &                      -FEpst(aep(3),hmE, Betop,hmF1))
            aep(2)=djoin(aep(2),0.8*NmF1,1.0,aep(2)
     &                                     -0.8*NmF1)
            aep(3)=4.0*(NmE -FEpst(aep(2),hmF1,B1bot,hmE)
     &                      -FEpst(aep(1),hmF2,B2bot,hmE))
         enddo
         if (aep(2).lt.0.0) then
            aep(2)=0.0
            aep(3)=4.0*(NmE-FEpst(aep(1),hmF2,B2bot,hmE))
         endif
      endif
      aep(3)=djoin(aep(3),0.5E-2,60.0,aep(3)-0.5E-2)
      b2k = 3.22-0.538E-1*foF2-0.664E-2*hmF2+0.113*hmF2/B2bot
     +    +0.257E-2*R12
      b2k=djoin(b2k,1.0,2.0,b2k-1.0)
      B2top=b2k*B2bot
      bb(1)=Bebot
      bb(2)=Betop
      bb(3)=B1bot
      bb(4)=B1top
      bb(5)=B2bot
      bb(6)=B2top
      return
      end
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 ----------------------------------------------------

C ----------------------------------------------------
      real function ef1r1(foF2,efoE)
      parameter (DR=1.745329252E-2)
      parameter (chi0=86.232928)
       efoF1=djoin(1.4*efoE,0.0,1000.0,efoE-2.0)
      efoF1=djoin(0.0,efoF1,1000.0,efoE-efoF1)
      ef1r1=djoin(efoF1,0.85*efoF1,60.0,0.85*foF2-efoF1)
      if (ef1r1.lt.1.0E-6) ef1r1=0.0
      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.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-----------------------------------------------------------------
      real function fmlt(ut,xlat,xlong,day)
C Calculats geomagnetic local time in hours for given
C universal time UT(sec)
C XLAT/XLONG=geodetic Latitude/Longitude (deg.)
C DAY= day of year
C D. Bilitza, COSPAR-1980, paper 7.3.1
C
	rad=57.29578
	flat=xlat/rad
	flong=xlong/rad
	delta=0.409207*sin((day-80.0)*3.14159/184.0)
	beta=(180.-ut/240.0)/rad
	px=cos(flat)*cos(flong)
	py=cos(flat)*sin(flong)
	pz=sin(flat)
	sx=cos(delta)*cos(beta)
	sy=cos(delta)*sin(beta)
	sz=sin(delta)
	p1=0.35117*px-0.91483*py-0.19937*pz
	p2=0.93358*px+0.35837*py
	s1=0.35117*sx-0.91483*sy-0.19937*sz
	s2=0.93358*sx+0.35837*sy
	thetp=atan2(p2,p1)
	thets=atan2(s2,s1)
	fmlt=(thetp-thets+3.141593)/.2617994
	afmlt=anint(fmlt*10.)
	fmlt=afmlt/10.
	if(fmlt.ge.24.) fmlt=fmlt-24.
	if(fmlt.lt.0.) fmlt=fmlt+24.
	return
	end
C
C ---------------------------------------------------------------------

