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         R27 - integral index of solar activity for 27 days - 1 revolution of the Sun (real)
c 
c    Output:
c     model_TEC - Global Model of Critical Frequency TEC (real)
c
c
      real function model_TEC(dayut,monthut,UT,lati,longi,SSN_R12,R27,
     &                                                       TECparam)
	implicit none
c     .. scalar arguments ..
	character(*) TECparam
      integer dayut, monthut
      real UT, lati, longi
	real SSN_R12, R27
c     .. function references ..
	real TEC_R27, TEC_SSN2
c
	select case(TECparam)
	  case('R27')	
          model_TEC=TEC_R27(dayut,monthut,UT,lati,longi,R27)
	  case('SSN2')
	    model_TEC=TEC_SSN2(dayut,monthut,UT,lati,longi,SSN_R12)
	end select
	
	return
      end function model_TEC
c
      real function TEC_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 TEC_0, TEC_m, TEC_p
	real quiet_TEC
c     .. function references ..
	real med_TEC_SSN2
c
	TEC_0 = med_TEC_SSN2(UT,monthut,index,lati,longi)
      if (dayut.le.15) then
        if (dayut.eq.15) then
          quiet_TEC = TEC_0
        else
          monthr = monthut-1
          if (monthr.eq.0) monthr = 12
	    TEC_m = med_TEC_SSN2(UT,monthr,index,lati,longi)
          quiet_TEC = (TEC_0-(dayut-15)*(TEC_m-TEC_0)/30.)
        end if
      else
        montha = mod(monthut,12) + 1
	  TEC_p = med_TEC_SSN2(UT,montha,index,lati,longi)
        quiet_TEC  = (TEC_0+(dayut-15)*(TEC_p-TEC_0)/30.)
      end if

	TEC_SSN2 = quiet_TEC

	return
      end function TEC_SSN2
c
c      Median of TEC 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_TEC_SSN2 - median of TEC for a given UT and monthut (real)
c
      real function med_TEC_SSN2(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 TEC_ut(0:23)
	common/TECUT/TEC_ut
c     .. function references ..
      real TEC_med_SSN2, fun_TEC_UT

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

      med_TEC_SSN2 = fun_TEC_UT(t) 
      return
      end function med_TEC_SSN2
c
c  Median of TEC 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     TEC_med_SSN2 - TEC for a given UT and monthut (real) with SSN2
c  
      real function TEC_med_SSN2(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, TEC_1, TEC_2
	real cov_Max
	double precision teta
c ..   local arrays ..
	double precision coeff_month(0:225,0:47)
      double precision Kf(0:225)
	real ft_TEC_1(12), ft_TEC_2(12)
c                                SSN2 2020
c                     Jan   Feb   Mar   Apr   May   Jun
      data ft_TEC_1/ 14.1, 13.4, 11.8,  6.4,  6.4,  6.0,
     *                6.2,  7.4,  9.5,  9.9,  9.2,  7.9/
      data ft_TEC_2/144.0,139.1,138.7,142.9,148.1,151.6,
     *              161.5,155.4,150.0,159.3,149.6,147.0/
c                     Jul   Aug   Sep   Okt   Nov   Dec
c
c     .. local in common ..
	double precision umr
	common/const_r/umr
c     .. function references ..
	real fun_TEC
c     .. subroutine references ..
c        read_TEC_SSN2, read_TEC

	umr=atan(1.0)*4./180

      teta = 90.0-dlat

	call read_TEC(monthut,coeff_month)

      Kf = coeff_month(0:225,iUT)
	TEC_1 = fun_TEC(teta,dlong,Kf)
      Kf = coeff_month(0:225,iUT+24)
	TEC_2 = fun_TEC(teta,dlong,Kf)
c
	cov1 = ft_TEC_1(monthut)
	cov2 = ft_TEC_2(monthut)
c
      a = (TEC_2 - TEC_1)/(cov2-cov1)
	b = TEC_2- a*cov2
c?????	
	cov_Max = 200.0
c?????	
	if (cov > cov_Max) then
	    cov_g = cov_Max 
	else
          cov_g = cov 
	end if
	
	TEC_med_SSN2 = a*cov_g + b
c
      return
      end function TEC_med_SSN2
c
c   Read median coefficients of TEC
c
c    Input:
c      monthut -(integer)
c   Output:
c .    coeff_month(0:225,0:47) - median coefficients of TEC (real*8) 
c
      subroutine read_TEC(month,coeff_month)
	implicit none
c     .. scalar arguments ..
	integer month
	double precision coeff_month(0:225,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:225,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,225)
         end do
	   close(10)
	   coeff_month_read(month) = 1
	endif

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

	return
   5	print 110, filedata
      stop
c
  10  format('data\kof_TEC\gtecsn',i2,'.dat')
  30  format(6(d12.5))
 100  format('Non file',1x,a)
 110  format('Corrupt file',1x,a)
      end subroutine read_TEC
c
c
      real function TEC_R27(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 TEC_0, TEC_m, TEC_p
	real quiet_TEC
c     .. function references ..
	real med_TEC_R27
c
	TEC_0 = med_TEC_R27(UT,monthut,index,lati,longi)
      if (dayut.le.15) then
        if (dayut.eq.15) then
          quiet_TEC = TEC_0
        else
          monthr = monthut-1
          if (monthr.eq.0) monthr = 12
	    TEC_m = med_TEC_R27(UT,monthr,index,lati,longi)
          quiet_TEC = (TEC_0-(dayut-15)*(TEC_m-TEC_0)/30.)
        end if
      else
        montha = mod(monthut,12) + 1
	  TEC_p = med_TEC_R27(UT,montha,index,lati,longi)
        quiet_TEC  = (TEC_0+(dayut-15)*(TEC_p-TEC_0)/30.)
      end if

	TEC_R27 = quiet_TEC

	return
      end function TEC_R27
c
c      Median of TEC for a given UT and monthut
c    Input:
c        UT - (real)
c   monthut - (integer)
c       cov - R27 (real)
c      dlat - geographic latitude, degrees (real)
c     dlong - geographic longitude, degrees (real)
c 
c   Output:
c   med_TEC_R27 - median of TEC for a given UT and monthut (real)
c
      real function med_TEC_R27(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 TEC_ut(0:23)
	common/TECUT/TEC_ut
c     .. function references ..
      real TEC_med_R27, fun_TEC_UT

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

      med_TEC_R27 = fun_TEC_UT(t) 
      return
      end function med_TEC_R27
c
c  Median of TEC at UT and month (R linear interpolation)
c
c     Input:
c        iUT - (integer)
c    monthut - (integer)
c        cov - R27 (real
c       dlat - geographic latitude, degrees (real)
c      dlong - geographic longitude, degrees (real)
c 
c    Output:
c     TEC_med_R27 - TEC for a given UT and monthut (real) with SSN2
c  
      real function TEC_med_R27(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, TEC_1, TEC_2
      real cov_Max
	double precision teta
c ..   local arrays ..
	double precision coeff_month(0:225,0:47)
      double precision Kf(0:225)
	real ft_TEC_1(12), ft_TEC_2(12)
c                                R27 2020
c                     Jan   Feb   Mar   Apr   May   Jun
      data ft_TEC_1/ 14.1, 13.4, 11.8,  6.4,  6.4,  6.0,
     *                6.2,  7.4,  9.5,  9.9,  9.2,  7.9/
      data ft_TEC_2/131.2,128.6,130.3,143.4,138.7,158.4,
     *              153.5,156.5,139.6,148.5,154.3,131.6/
c                     Jul   Aug   Sep   Okt   Nov   Dec
c     .. local in common ..
	double precision umr
	common/const_r/umr
c     .. function references ..
	real fun_TEC
c     .. subroutine references ..
c        read_TEC_R27, read_TEC

	umr=atan(1.0)*4./180

      teta = 90.0-dlat

	call read_TEC(monthut,coeff_month)

      Kf = coeff_month(0:225,iUT)
	TEC_1 = fun_TEC(teta,dlong,Kf)
      Kf = coeff_month(0:225,iUT+24)
	TEC_2 = fun_TEC(teta,dlong,Kf)
c
	cov1 = ft_TEC_1(monthut)
	cov2 = ft_TEC_2(monthut)
c
	a = (TEC_2 - TEC_1)/(cov2-cov1)
	b = TEC_2- a*cov2
      cov_Max = 200.0 
	if (cov > cov_Max) then
	    cov_g = cov_Max 
	else
          cov_g = cov 
	end if

	TEC_med_R27 = a*cov_g + b
c
      return
      end function TEC_med_R27
c  
c        function fun_TEC_UT
c    Furie expansion of mediam TEC on  UT 
c . Input:
c.    t - time (real*8)
c .  Output = function fun_TEC_UT:
c .   TEC at arbitrary time UT (real)
c
c
      real function fun_TEC_UT(t)
      implicit none
c  .. scalar arguments ..
      double precision t
c  .. local scalars ..
      integer k
      double precision TEC
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_UTt, fun_Gk_UTt
c
      dtr=atan(1.0)*4.0/12.0
	call Koeff_UTt(Kf_UT)
      call fun_Gk_UTt(t,Gk_UT)
	TEC = 0.d0
	do k=0,6
	   TEC = TEC + Kf_UT(k)*Gk_UT(k) 
	end do
c
      fun_TEC_UT = TEC
      return
      end function fun_TEC_UT
c
c        fun_TEC_UT
c     Auxiliary function Koeff_UTt
c
      subroutine Koeff_UTt(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_UTt
c
      mm = 3
	mk = 2*mm
      call fun_Akp_UTt(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_UTt
c
c        fun_TEC_UT
c     Auxiliary function fun_Akp_UTt
c
      subroutine fun_Akp_UTt(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 TEC_ut(0:23)
	common/TECUT/TEC_ut
c     .. subroutine references ..
c        fun_Gk_UTt, fun_Fk_UTt
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_UTt(t,Gk_UT)
	         call fun_Fk_UTt(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 + TEC_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_UTt(t,Gk_UT)
         call fun_Fk_UTt(Gk_UT,Akp_UT,Fk_UT)
         sum_Dn = sum_Dn + TEC_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_UTt
c
c        fun_TEC_UT
c     Auxiliary function fun_Fk_UTt
c
      subroutine fun_Fk_UTt(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_UTt
c
c        fun_TEC_UT
c     Auxiliary function fun_Gk_UTt
c
	subroutine fun_Gk_UTt(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_UTt
c
c    Spatial interpolation with Spherical Harmonic Functions
c    of TEC median
c . Input:
c.    teta - colatitude (real*8)
c.    long - longitude  (real)
c .  Output function fun_TEC:
c.      Kf - array of expansion coefficients (real*8)
c
      real function fun_TEC(teta,long,Kf)
      implicit none
! .. scalar arguments ..
      double precision teta
      real long
      double precision Kf(0:225)
! .. local scalars ..
      integer k
      double precision TEC
! .. local arrays ..
      double precision Gk(0:225)
c     .. subroutine references ..
c        fun_Gkt

	call fun_Gkt(teta,long,Gk)
	TEC = 0.d0
	do k=0,225
	   TEC = TEC + Kf(k)*Gk(k) 
	end do
      fun_TEC = TEC
      return
      end function fun_TEC
c
c       function Gk
c    Spherical Harmonic Eaxpansion
c . Input:
c.    teta - colatitude (real*8)
c.    long - longitude  (real)
c . Output:
c.      Gk - array of expansion coefficients (real*8)
c
      subroutine fun_Gkt(teta,long,Gk)
      implicit none
! .. scalar arguments ..
      double precision teta
      real long
      double precision Gk(0:225)
! .. local scalars ..
      integer mm, nn, m, n, k
! .. local arrays ..
	double precision Pl_mn(0:10,0:15)
!  .. local in common ..
	double precision umr
	common/const_r/umr
c     .. subroutine references ..
c        apoly_Legendre

      Pl_mn = 0.d0
	mm = 10
	nn = 15
      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_Gkt
c
c        fun_Gk
c
c     Indices
c  
c     Input:
c       dayut - (integer)
c     monthut - (integer)
c      yearut - (integer)
c          UT - (real)
c
c     Output:
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-      subroutine indices(dayut,monthut,yearut,UT,Ri,SSN_R12,RR)
      subroutine indices(dayut,monthut,yearut,Ri,SSN_R12,RR
	+                                   ,INY_R,F2param)		!++TLG		
	implicit none
c     .. scalar arguments ..
      integer dayut, monthut, yearut
		character(4) F2param			  !++TLG
c-	real UT
	real Ri, SSN_R12, RR
c     .. local scalars ..
      integer ndayut,iny_R			   !++TLG
C-	+, iUT
C?      real R27
	      real Rt
c     .. subroutine references ..
c        nnday, indices_R, R_tau 
c
	call nnday(dayut,monthut,yearut,ndayut)
c-      call indices_R(monthut,yearut,SSN_R12)
		select case(F2param)
		  case('SSN2')
C      call indices_R(monthut,yearut,SSN_R12)
      call indices_R(monthut,yearut,SSN_R12,INY_R)  !+TLG
	RR=SSN_R12
	 case('R27 ')	
C	      call R_tau(dayut,monthut,yearut,Ri,Rt)
      call R_tau(dayut,monthut,yearut,Ri,Rt,iny_R) !+TLG
C	      R27 = (2*Rt+Ri)/3
       RR = (2*Rt+Ri)/3
			end select
c-      call indices_R(monthut,yearut,SSN_R12,INY_R)  !+TLG
c-      call R_tau(dayut,monthut,yearut,Ri,R27)
c-      RR = (2*R27+Ri)/3

c-      iUT = int(UT)
c
	return
 	end subroutine indices
c
      subroutine indices_iny_R12(INY_R12,INY_R)  !+TLG
c-	use Mod_SSN
	implicit none
c     .. local scalars ..
      integer kd
	integer w_year,w_month,w_day,i_R
       integer iny_R12					   !+TLG
        integer iny_R						   !+TLG
	real year_decimal
c-	+, w_SSN
	  character(256) file_R12 			   !+TLG
      integer n_R12						   !+TLG
	+,iupd,iupm,iupy						   !+TLG
	+,imst, imend, iyend					   !+TLG
  	  data n_R12/212/					   !+TLG
      data file_R12/'ssn2_12.dat'/		! TLG
        integer n_R						!+TLG
	  data n_R/214/						!+TLG
	  character(256) file_R 			!+TLG
	  data file_R/'SN_R.dat'/           ! TLG
c     .. local arrays ..
c     ..
c-      open(unit=n_R12,file=file_R12,iostat=kd,STATUS='OLD')
      open(unit=n_R12,file=file_R12,STATUS='OLD')
c-	      if (kd.ne.0) then
c-          print 80, 'SSN_R12.dat'
c-          stop
c-      end if
      
c-	read(n_R12,10) w_year,w_month,year_decimal,w_SSN
          read(n_R12,*) iupd,iupm,iupy			   !+TLG
          read(n_R12,*) imst,w_year, imend, iyend	   !+TLG
	iny_R12 = w_year
	close(n_R12)
c
      open(unit=n_R,file=file_R,iostat=kd,STATUS='OLD')
      if (kd.ne.0) then
          print 80, 'SN_R.dat'
          stop
      end if
      
	read(n_R,20) w_year,w_month,w_day,year_decimal,i_R
	iny_R = w_year
	close(n_R)
c
	return
c-  10  format(i4,1x,i2,1x,f8.3,1x,f6.1)
   10  format(1X,I4)  				 !+TLG
  20  format(i4,2(1x,i2),1x,f8.3,1x,i4)
  80  format(4x,'Not file:',a50)
  99  format(2x,'Year = ',i4,2x,'nday = ',i3) 	
 	end subroutine indices_iny_R12
c
c .. Input:
c ..    day                           (integer)
c .    month                          (integer)
c .     year                          (integer)
c .     nday - day number of the year (integer)
c
c .  Output:
c .    Rt - integral index of solar activity for 27 days - 1 revolution of the Sun (real)
c
C      subroutine R_tau(day,month,year,Ri,Rt)
      subroutine R_tau(day,month,year,Ri,Rt,INY_R)
C	      use Mod_SSN
	implicit none
c     .. scalar arguments ..
      integer day,month,year
      real Ri, Rt
c     .. local scalars ..
C	integer kd, i_R
	integer  i_R
	      integer jd
	integer nday
	integer ns, dyear, irec, nd
	integer w_day,w_month,w_year,nn
	real year_decimal, SD
      real Wsum, Wk, wta, wta0, ww
      real Wsumo
      real:: tauR27=0.964  ! 27 day
	integer issn(0:28),k
        integer iny_R
        integer n_R
	  data n_R/214/
	  character(256) file_R 
	  data file_R/'SN_R.dat'/
c      
C      open(unit=n_R,file=file_R,iostat=kd,
C     *    ACCESS='DIRECT',RECL=38,FORM='FORMATTED',STATUS='OLD')
      open(unit=n_R,file=file_R,action='READ')   ! +TLG
C	      if (kd.ne.0) then
C          print 80, 'SN_R.dat'
C          stop
C      end if
c      
	jd=mod(iny_R-1,4)
	wta0=1.0e-03
	NS = int(1.0/(1-tauR27))
      dyear = year - iny_R
	call nnday(day,month,year,nday)
	irec = dyear*365 + (dyear+jd)/4 + nday
      ww = 1.0
	wta = 1.0
	Wsumo = 0.0
	Wsum = 0.0
C      do nd=0,NS
C	   read(n_R,10,rec=irec,err=5) w_year,w_month,w_day,year_decimal,
C     &                                                       i_R,SD,nn
    1	read(n_R,10,err=5) w_year,w_month,w_day,year_decimal,i_R   ! TLG
     +,SD,nn													   ! TLG
C  NEW TLG:
      issn(0)=i_R
	if ((w_year.eq.year).and.(w_month.eq.month).and.(w_day.eq.day)) 
     +	then
	 goto 2
	    else
C Move data forward by 1 day:
	do k=28,1,-1
	issn(k)=issn(k-1)
	enddo
	goto 1
	endif
    2          do nd=0,NS
C             if (nd.eq.0) Ri=i_R
         if (nd.eq.0) Ri=issn(0)
C	            irec = irec - 1
	Wk=issn(nd+1)
C	         Wk = i_R
	   if (tauR27.le.0.) then
            Rt =  Wk
C-            return
       goto 3
         end if
	   Wsumo = Wsumo + wta
	   if (nd==0) then
	      Wsum = Wsum + Wk
	   else
	      Wsum = Wsum + Wk * wta
	   end if
	   wta = wta * tauR27
      enddo
      Rt = Wsum/Wsumo
   3	close(n_R)
c
      return
c
   5 	print 100
	stop
  10  format(i4,2(1x,i2),1x,f8.3,1x,i4,1x,f5.1,1x,i4)
  80  format(4x,'Not file:',a50)
 100  format(4x,'Error in file: SN_R.dat')
      end subroutine R_tau 
C      subroutine indices_R(month,year,SSN_R12)
      subroutine indices_R(month,year,SSN_R12,INY_R12)	 ! +TLG
C		use Mod_SSN
	implicit none
c     .. scalar arguments ..
        integer iny_R12				!+TLG
      integer month, year 
     +,iupd,iupm,iupy,imst,iyst, imend, iyend	 !	+TLG
	real SSN_R12
c     .. local scalars ..
C	integer dyear, ir, kd
	integer dyear, ir , kk       !+TLG
C	integer w_year,w_month,nn
	integer w_year,w_month  !+TLG
C	real year_decimal, Sd
        integer n_R12
	  data n_R12/212/
	real r12mn(12) !+TLG
		  character(256) file_R12 
      data file_R12/'ssn2_12.dat'/		! TLG
c
      dyear=year-iny_R12
      ir = dyear*12 + month
C      open(unit=n_R12,file=file_R12,iostat=kd,
C     *    ACCESS='DIRECT',RECL=38,FORM='FORMATTED',STATUS='OLD')
      open(unit=n_R12,file=file_R12,action='READ')  !+TLG
c Read the update date, the start date and the end date (mm,yyyy)
c 
          read(n_R12,*) iupd,iupm,iupy						  !+TLG
          read(n_R12,*) imst,iyst, imend, iyend				  !+TLG
C      if (kd.ne.0) then
C          print 80, 'SSN_R12.dat'
C          stop
C      end if
c     
   1     	read(n_R12,*,err=40) w_year,(r12mn(kk),kk=1,12)  ! TLG
       w_month=month
      SSN_R12=r12mn(month)  !+TLG

C	read(n_R12,10,rec=ir,err=40) w_year,w_month,year_decimal,SSN_R12,
C     &                                                           SD,nn
C	if (w_year.ne.year.and.w_month.ne.month) then
	if (w_year.ne.year) then
C	   print 99, w_year,w_month
	goto 1
C	   stop
	endif

c	
	close(n_R12)
	return
c
  40  print 100
      stop
C  10  format(i4,1x,i2,1x,f8.3,1x,2(f6.1),2x,i4)
   10 format(1X,I4,1X,12(1X,F5.1))				 !+TLG
  80  format(4x,'Not file:',a50)
  99  format(2x,'Year = ',i4,2x,'nday = ',i3)
 100  format(4x,'Error in file: SSN_R12.dat')  
 	end subroutine indices_R