c
c     Global Model of Critical Frequency of the Ionospheric F2 Layer
c
c     Input:
c      dayut - (integer)
c    monthut - (integer)
c         UT - (real)
c       lati - geographic latitude, degrees (real)
c      longi - geographic longitude, degrees (real)
c    SSN_R12 - index of solar activity at a given day (real)
c         RR - integral index of solar activity for 27 days - 1 revolution of the Sun (real)
c 
c    Output:
c     model_GMF2R - Global Model of Critical Frequency foF2 (real)
c
c
      real function model_GMF2Rf3(dayut,monthut,UT,lati,longi,SSN_R12,RR
     &                                                         ,F2param)
	implicit none
c     .. scalar arguments ..
	character(*) F2param
      integer dayut, monthut
      real UT, lati, longi
	real SSN_R12, RR
c     .. function references ..
	real fof2_R_RR, fof2_R_SSN2
c
	select case(F2param)
	  case('R27')	
          model_GMF2Rf3=fof2_R_RR(dayut,monthut,UT,lati,longi,RR)
	  case('SSN2')
	    model_GMF2Rf3=fof2_R_SSN2(dayut,monthut,UT,lati,longi,SSN_R12)
	end select
	
	return
      end function model_GMF2Rf3
c
      real function fof2_R_SSN2(dayut,monthut,UT,lati,longi,index)
	implicit none
c     .. scalar arguments ..
      integer dayut, monthut
      real UT, lati, longi
	real index
c     .. local scalars ..
      integer monthr, montha
	real max_F2_0, max_F2_m, max_F2_p
	real quiet_F2
c     .. function references ..
	real med_max_F2_SSN2f
c
	max_F2_0 = med_max_F2_SSN2f(UT,monthut,index,lati,longi)
      if (dayut.le.15) then
        if (dayut.eq.15) then
          quiet_F2 = max_F2_0
        else
          monthr = monthut-1
          if (monthr.eq.0) monthr = 12
	    max_F2_m = med_max_F2_SSN2f(UT,monthr,index,lati,longi)
          quiet_F2 = (max_F2_0-(dayut-15)*(max_F2_m-max_F2_0)/30.)
        end if
      else
        montha = mod(monthut,12) + 1
	  max_F2_p = med_max_F2_SSN2f(UT,montha,index,lati,longi)
        quiet_F2  = (max_F2_0+(dayut-15)*(max_F2_p-max_F2_0)/30.)
      end if

	fof2_R_SSN2 = quiet_F2

	return
      end function fof2_R_SSN2
c
c      Median of foF2 for a given UT and monthut
c    Input:
c        UT - (real)
c   monthut - (integer)
c       cov - SSN_R12 (real)
c      dlat - geographic latitude, degrees (real)
c     dlong - geographic longitude, degrees (real)
c 
c   Output:
c   med_max_F2_SSN2f - median of foF2 for a given UT and monthut (real)
c
c   Shubin, V.N. Geomagnetism and Aeronomy. V.57, N4. P.414–425. 2017. 
c   DOI: 10.1134/S0016793217040181
c
      real function med_max_F2_SSN2f(UT,monthut,cov,dlat,dlong)
      implicit none
c     .. scalar arguments ..
      integer monthut
      real cov
      real UT
      real dlat, dlong
c     .. local scalars ..
      integer i
      double precision t
c     .. local arrays ..
	double precision xUT(0:23)
c  .. array in common ..
	double precision max_F2_ut(0:23)
	common/max_F2UT/max_F2_ut
c     .. function references ..
      real maxF2_med_SSN2f, fun_maxF2UTf

      max_F2_ut = 0.0
	do i=0,23
         max_F2_ut(i) = maxF2_med_SSN2f(i,monthut,cov,dlat,dlong)
	   xUT(i) = dble(i)
      end do
 
      t = dble(UT)

      med_max_F2_SSN2f = fun_maxF2UTf(t) 
      return
      end function med_max_F2_SSN2f
c
c  Median of foF2 at UT and month (R linear interpolation)
c
c     Input:
c        iUT - (integer)
c    monthut - (integer)
c        cov - SSN_R12 (real)
c       dlat - geographic latitude, degrees (real)
c      dlong - geographic longitude, degrees (real)
c 
c    Output:
c     maxF2_med_SSN2f - foF2 for a given UT and monthut (real) with SSN2
c  
      real function maxF2_med_SSN2f(iUT,monthut,cov,dlat,dlong)
      implicit none
c ..   scalar arguments ..
	integer monthut, iUT
      real cov
	real dlat, dlong
c ..   local scalars ..
	real cov1, cov2, cov_g
      real a, b, maxF2_1, maxF2_2
	real fexp1, fexp2
	double precision teta
c ..   local arrays ..
	double precision coeff_month(0:148,0:47)
      double precision Kf(0:148)
	real ft_fof2_1(12), ft_fof2_2(12)
c
c                                SSN2 2020
c                     Jan   Feb   Mar   Apr   May   Jun
      data ft_fof2_1/ 14.0, 16.1, 17.7, 15.1, 11.2, 12.8,
     *                15.1,  9.9,  9.6, 14.3, 15.0, 12.7/
      data ft_fof2_2/113.8,110.5,114.3,116.4,115.0,114.3,
     *               112.6,108.3,104.6,107.0,106.9,107.6/
c                     Jul   Aug   Sep   Okt   Nov   Dec
c
c     .. local in common ..
	double precision umr
	common/const_r/umr
c     .. function references ..
	real fun_maxF2f
c     .. subroutine references ..
c        read_foF2_SSN2, read_hmF2

	umr=atan(1.0)*4./180

      teta = 90.0-dlat

	call read_foF2(monthut,coeff_month)

      Kf = coeff_month(0:148,iUT)
	maxF2_1 = fun_maxF2f(teta,dlong,Kf)
      Kf = coeff_month(0:148,iUT+24)
	maxF2_2 = fun_maxF2f(teta,dlong,Kf)
c
	cov1 = ft_fof2_1(monthut)
	cov2 = ft_fof2_2(monthut)
c
	select case(1)
	  case(1)
	    a = (maxF2_2 - maxF2_1)/(cov2-cov1)
	    b = maxF2_2- a*cov2
		  if (cov>180.0) then
	        cov_g = 180.0 
	    else
              cov_g = cov 
	    end if
	    maxF2_med_SSN2f = a*cov_g + b
	  case(2)
          fexp1 = exp(-(cov1/(cov2-cov1)))  
          fexp2 = exp(-(cov2/(cov2-cov1)))
	    a = (maxF2_2 - maxF2_1)/(fexp2-fexp1)
	    b = maxF2_2- a*fexp2
	    maxF2_med_SSN2f = a*exp(-(cov/(cov2-cov1))) + b
      end select
c
      return
      end function maxF2_med_SSN2f
c
c   Read median coefficients of foF2
c
c    Input:
c      monthut -(integer)
c   Output:
c .    coeff_month(0:148,0:47) - median coefficients of foF2 (real*8) 
c
      subroutine read_foF2(month,coeff_month)
	implicit none
c     .. scalar arguments ..
	integer month
	double precision coeff_month(0:148,0:47)
c     .. local scalars ..
	character(256) filedata
	integer i, j, kd
c     ..   local arrays ..
	integer coeff_month_read(1:12)
c
	double precision coeff_month_all(0:148,0:47,1:12)
	save coeff_month_all
	data coeff_month_read /12*0/
c
      if (coeff_month_read(month).eq.0) then
         write(filedata,10) month+10
         open(10, File=filedata,Iostat=kd,status='old')
	   if (kd.ne.0) then
	      print 100, filedata
            stop
	   end if
	   do j=0,47
	      read(10,30,err=5) (coeff_month_all(i,j,month),i=0,148)
         end do
	   close(10)
	   coeff_month_read(month) = 1
	endif

      coeff_month = coeff_month_all(0:148,0:47,month)	

	return
   5	print 110, filedata
      stop
c
  10  format('data\foF2\fgssn',i2,'.dat')
  30  format(6(d12.5))
 100  format('Non file',1x,a)
 110  format('Corrupt file',1x,a)
      end subroutine read_foF2
c
c
      real function fof2_R_RR(dayut,monthut,UT,lati,longi,index)
	implicit none
c     .. scalar arguments ..
      integer dayut, monthut
      real UT, lati, longi
	real index
c     .. local scalars ..
      integer monthr, montha
	real max_F2_0, max_F2_m, max_F2_p
	real quiet_F2
c     .. function references ..
	real med_max_F2_RRf
c
	max_F2_0 = med_max_F2_RRf(UT,monthut,index,lati,longi)
      if (dayut.le.15) then
        if (dayut.eq.15) then
          quiet_F2 = max_F2_0
        else
          monthr = monthut-1
          if (monthr.eq.0) monthr = 12
	    max_F2_m = med_max_F2_RRf(UT,monthr,index,lati,longi)
          quiet_F2 = (max_F2_0-(dayut-15)*(max_F2_m-max_F2_0)/30.)
        end if
      else
        montha = mod(monthut,12) + 1
	  max_F2_p = med_max_F2_RRf(UT,montha,index,lati,longi)
        quiet_F2  = (max_F2_0+(dayut-15)*(max_F2_p-max_F2_0)/30.)
      end if

	fof2_R_RR = quiet_F2

	return
      end function fof2_R_RR
c
c      Median of foF2 for a given UT and monthut
c    Input:
c        UT - (real)
c   monthut - (integer)
c       cov - RR (real)
c      dlat - geographic latitude, degrees (real)
c     dlong - geographic longitude, degrees (real)
c 
c   Output:
c   med_max_F2_RRf - median of foF2 for a given UT and monthut (real)
c
      real function med_max_F2_RRf(UT,monthut,cov,dlat,dlong)
      implicit none
c     .. scalar arguments ..
      integer monthut
      real cov
      real UT
      real dlat, dlong
c     .. local scalars ..
      integer i
      double precision t
c     .. local arrays ..
	double precision xUT(0:23)
c  .. array in common ..
	double precision max_F2_ut(0:23)
	common/max_F2UT/max_F2_ut
c     .. function references ..
      real maxF2_med_RRf, fun_maxF2UTf

      max_F2_ut = 0.0
	do i=0,23
         max_F2_ut(i) = maxF2_med_RRf(i,monthut,cov,dlat,dlong)
	   xUT(i) = dble(i)
      end do
 
      t = dble(UT)

      med_max_F2_RRf = fun_maxF2UTf(t) 
      return
      end function med_max_F2_RRf
c
c  Median of foF2 at UT and month (R linear interpolation)
c
c     Input:
c        iUT - (integer)
c    monthut - (integer)
c        cov - RR (real
c       dlat - geographic latitude, degrees (real)
c      dlong - geographic longitude, degrees (real)
c 
c    Output:
c     maxF2_med_RRf - foF2 for a given UT and monthut (real) with SSN2
c  
      real function maxF2_med_RRf(iUT,monthut,cov,dlat,dlong)
      implicit none
c ..   scalar arguments ..
	integer monthut, iUT
      real cov
	real dlat, dlong
c ..   local scalars ..
	real cov1, cov2, cov_g
      real a, b, maxF2_1, maxF2_2
	real fexp1, fexp2
	double precision teta
c ..   local arrays ..
	double precision coeff_month(0:148,0:47)
      double precision Kf(0:148)
	real ft_fof2_1(12), ft_fof2_2(12)
c                                R27 2020
c                     Jan   Feb   Mar   Apr   May   Jun
      data ft_fof2_1/ 10.6, 15.3,  9.0,  5.9, 11.1, 10.9,
     *                12.7,  9.3,  5.1,  8.9, 10.9, 12.2/
      data ft_fof2_2/123.7,125.4,123.4,107.5,115.9,118.1,
     *               112.8,106.9,111.6,116.1,125.1,123.3/
c                     Jul   Aug   Sep   Okt   Nov   Dec
c     .. local in common ..
	double precision umr
	common/const_r/umr
c     .. function references ..
	real fun_maxF2f
c     .. subroutine references ..
c        read_foF2_RR, read_hmF2

	umr=atan(1.0)*4./180

      teta = 90.0-dlat

	call read_foF2(monthut,coeff_month)

      Kf = coeff_month(0:148,iUT)
	maxF2_1 = fun_maxF2f(teta,dlong,Kf)
      Kf = coeff_month(0:148,iUT+24)
	maxF2_2 = fun_maxF2f(teta,dlong,Kf)
c
	cov1 = ft_fof2_1(monthut)
	cov2 = ft_fof2_2(monthut)
c
c
	select case(1)
	  case(1)
	    a = (maxF2_2 - maxF2_1)/(cov2-cov1)
	    b = maxF2_2- a*cov2
		  if (cov>200.0) then
	        cov_g = 200.0 
	    else
              cov_g = cov 
	    end if
	    maxF2_med_RRf = a*cov_g + b
	  case(2)
          fexp1 = exp(-(cov1/(cov2-cov1)))  
          fexp2 = exp(-(cov2/(cov2-cov1)))
	    a = (maxF2_2 - maxF2_1)/(fexp2-fexp1)
	    b = maxF2_2- a*fexp2
	    maxF2_med_RRf = a*exp(-(cov/(cov2-cov1))) + b
      end select
c
      return
      end function maxF2_med_RRf
c      
c
c    Input:
c      year - (integer)
c      nday - day od year (integer)
c
c    Output:
c   SSN_R12 - index of solar activity at a given day (real)
c        Ri - index of solar activity at a given day (real)
c       R27 - integral index of solar activity for 27 days - 1 revolution of the Sun (real)
c
c-	module Mod_SSN
c-	implicit none
c-        integer iny_R12
c-        integer iny_R
        	  
c=	  character(256) file_R12 
c-	  data file_R12/'data/Indices/SSN_R12.dat'/
c-        integer n_R12
c-	  data n_R12/212/

c-	  character(256) file_R 
c-	  data file_R/'data/Indices/SN_R.dat'/
c-        integer n_R
c-	  data n_R/214/

c-	end module Mod_SSN
c	
c
c
c
c
c
c  
c        функция fun_maxF2UTf
c    Разложение по  UT в ряды Фурье медиан foF2 или hmF2
c . Входной параметр:
c.    t - время (real*8)
c .  Выходные значения функции fun_maxF2UTf:
c .   foF2 или hmF2 в произвольный момент UT (real)
c
      real function fun_maxF2UTf(t)
      implicit none
c  .. scalar arguments ..
      double precision t
c  .. local scalars ..
      integer k
      double precision maxF2
c     .. local arrays ..	
	double precision Gk_UT(0:6)
	double precision Kf_UT(0:6)
c  .. local in common ..	
	double precision dtr
	common/radUT/dtr
c     .. subroutine references ..
c        Koeff_UTf, fun_Gk_UTf
c
      dtr=atan(1.0)*4.0/12.0
	call Koeff_UTf(Kf_UT)
      call fun_Gk_UTf(t,Gk_UT)
	maxF2 = 0.d0
	do k=0,6
	   maxF2 = maxF2 + Kf_UT(k)*Gk_UT(k) 
	end do
c
      fun_maxF2UTf = maxF2
      return
      end function fun_maxF2UTf
c
c        fun_maxF2UT
c     Вспомогательная функция Koeff_UTf
c
      subroutine Koeff_UTf(Kf_UT)
      implicit none
c  .. array arguments ..
	double precision Kf_UT(0:6)
c  .. local scalars ..
      integer mm, mk, k, m
      double precision sum_D
c     .. local arrays ..
	double precision Akp_UT(0:6,0:6)
	double precision Dk_UT(0:6)
c     .. subroutine references ..
c        fun_Akp_UTf
c
      mm = 3
	mk = 2*mm
      call fun_Akp_UTf(Akp_UT,Dk_UT)
      Kf_UT = 0.d0
      do k=mk,0,-1
         sum_D = 0.d0
         do m=k+1,mk
            sum_D = sum_D + Akp_UT(m,k)*Kf_UT(m)
         end do
         Kf_UT(k) = sum_D + Dk_UT(k)
      end do   
      return
	end subroutine Koeff_UTf
c
c        fun_maxF2UT
c     Вспомогательная функция fun_Akp_UTf
c
      subroutine fun_Akp_UTf(Akp_UT,Dk_UT)
      implicit none
c  .. array arguments ..
	double precision Akp_UT(0:6,0:6), Dk_UT(0:6)
c  .. local scalars ..
	integer mm, mk, i, k, p
	double precision t
	double precision sum_An, sum_Dn
	double precision sum_Ad, sum_Dd
c     .. local arrays ..
	double precision Gk_UT(0:6), Fk_UT(0:6)
c  .. array in common ..	
	double precision foF2_ut(0:23)
	common/max_F2UT/foF2_ut
c     .. subroutine references ..
c        fun_Gk_UTf, fun_Fk_UTf
c
      mm = 3
      mk = 2*mm
c
      Gk_UT = 0.d0
      Gk_UT(0) = 1.0
	Fk_UT = 0.d0
      Fk_UT(0) = 1.d0
	Akp_UT = 0.d0
	Dk_UT = 0.d0
	do p=0,mk
         sum_Dn=0.d0
         sum_Dd=0.d0
         do k=p+1,mk
            sum_An=0.d0
	      sum_Ad=0.d0
	      do i=0,23
               t = dble(i)
               call fun_Gk_UTf(t,Gk_UT)
	         call fun_Fk_UTf(Gk_UT,Akp_UT,Fk_UT)
		     sum_An = sum_An + Gk_UT(k)*Fk_UT(p)
		     sum_Ad = sum_Ad + Fk_UT(p)*Fk_UT(p)
			 if (p.eq.(k-1)) then
			     sum_Dn = sum_Dn + foF2_ut(i)*Fk_UT(p)
                   sum_Dd = sum_Dd + Fk_UT(p)*Fk_UT(p)
               end if
            end do
            Akp_UT(k,p) = - sum_An/sum_Ad
	   end do
         if (p.lt.mk) then
	      Dk_UT(p) = sum_Dn/sum_Dd
         end if
	end do

      p=mk
      sum_Dn=0.d0
      sum_Dd=0.d0
      do i=0,23
         t = dble(i)
         call fun_Gk_UTf(t,Gk_UT)
         call fun_Fk_UTf(Gk_UT,Akp_UT,Fk_UT)
         sum_Dn = sum_Dn + foF2_ut(i)*Fk_UT(p)
         sum_Dd = sum_Dd + Fk_UT(p)*Fk_UT(p)
      end do
      Dk_UT(p) = sum_Dn/sum_Dd
      return
      end subroutine fun_Akp_UTf
c
c        fun_maxF2UT
c     Вспомогательная функция fun_Fk_UTf
c
      subroutine fun_Fk_UTf(Gk_UT,Akp_UT,Fk_UT)
      implicit none
c  .. array arguments ..
	double precision Gk_UT(0:6)
	double precision Akp_UT(0:6,0:6)
	double precision Fk_UT(0:6)
c  .. local scalars ..
      integer k, p
      double precision sum_G

      Fk_UT = 0.d0
      do k=0,6
         sum_G = 0.d0
         do p=0,k
            if (k.eq.p) cycle
            sum_G = sum_G + Akp_UT(k,p)*Fk_UT(p)
         end do
         Fk_UT(k) = sum_G + Gk_UT(k)
      end do
c
      return
      end subroutine fun_Fk_UTf
c
c        fun_maxF2UT
c     Вспомогательная функция fun_Gk_UTf
c
	subroutine fun_Gk_UTf(t,Gk_UT)
      implicit none
c  .. scalar arguments ..
      double precision t
c  .. array arguments ..
	double precision Gk_UT(0:6)
c  .. local scalars ..
      integer m, mm, k
c  .. local in common ..
	double precision dtr
	common/radUT/dtr
c
      mm = 3
	Gk_UT = 0.d0
      k = 0
      do m=0,mm
         if (m.eq.0) then
            Gk_UT(k) = 1
            k = k + 1
         else 
            Gk_UT(k)   = cos(m*t*dtr)
            Gk_UT(k+1) = sin(m*t*dtr)
		  k = k + 2
         end if
      end do
      return
      end subroutine fun_Gk_UTf
c
c    Пространствеенная интерполяция по сферическим гармоническим функциям
c    медианы foF2 или hmF2 
c . Входные параметры:
c.    teta - коширота (real*8)
c.    long - долгота  (real)
c .  Выходные значения функции fun_maxF2f:
c.      Kf - массив коэффициентов разложения (real*8)
c
      real function fun_maxF2f(teta,long,Kf)
      implicit none
! .. scalar arguments ..
      double precision teta
      real long
      double precision Kf(0:148)
! .. local scalars ..
      integer k
      double precision maxF2
! .. local arrays ..
      double precision Gk(0:148)
c     .. subroutine references ..
c        fun_Gkf

	call fun_Gkf(teta,long,Gk)
	maxF2 = 0.d0
	do k=0,148
	   maxF2 = maxF2 + Kf(k)*Gk(k) 
	end do
      fun_maxF2f = maxF2
      return
      end function fun_maxF2f
c
c       функции Gk
c    Разложение по сферическим гармоническим функциям
c . Входные параметры:
c.    teta - коширота (real*8)
c.    long - долгота  (real)
c . Выходной параметр:
c.      Gk - массив коэффициентов разложения (real*8)
c
      subroutine fun_Gkf(teta,long,Gk)
      implicit none
! .. scalar arguments ..
      double precision teta
      real long
      double precision Gk(0:148)
! .. local scalars ..
      integer mm, nn, m, n, k
! .. local arrays ..
	double precision Pl_mn(0:8,0:12)
!  .. local in common ..
	double precision umr
	common/const_r/umr
c     .. subroutine references ..
c        apoly_Legendre

      Pl_mn = 0.d0
	mm = 8
	nn = 12
      call apoly_Legendre(mm,nn,teta,Pl_mn)
      Gk = 0.d0
      k = 0
      do m=0,mm
         if (m==0) then
	      do n=0,nn
              Gk(k) = Pl_mn(m,n)
			k = k + 1
            end do
         else 
	      do n=m,nn
              Gk(k)   = Pl_mn(m,n)*cos(m*long*umr)
		      Gk(k+1) = Pl_mn(m,n)*sin(m*long*umr)
			k = k + 2
           end do
         end if
      end do
      return
      end subroutine fun_Gkf
c
c        fun_Gk