C TLG..........................................................Nov 2024 C Make hmF2R12 similar to foF2R12 C c c Global Model of Peak height hmF2 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_hmF2 - Global Model of hmF2 (real) c c real function model_GMF2Rh3(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? real quiet_F2 C-- +dhmF2 c .. function references .. real hmF2_R_RR , hmF2_R_SSN2 C- +,model_GMF2R_hmF2 C-- +, storm_F2 c select case(F2param) case('R27') model_GMF2Rh3=hmF2_R_RR(dayut,monthut,UT,lati,longi,RR) C? quiet_F2 = model_GMF2R_hmF2 case('SSN2') model_GMF2Rh3=hmF2_R_SSN2(dayut,monthut,UT,lati,longi, & SSN_R12) end select return end function model_GMF2Rh3 c real function hmF2_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_SSN2 c max_F2_0 = med_max_F2_SSN2(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_SSN2(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_SSN2(UT,montha,index,lati,longi) quiet_F2 = (max_F2_0+(dayut-15)*(max_F2_p-max_F2_0)/30.) end if hmF2_R_SSN2 = quiet_F2 return end function hmF2_R_SSN2 c c Median of hmF2 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_SSN2 - median of hmF2 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_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 max_F2_ut(0:23) common/max_F2UT/max_F2_ut c .. function references .. real maxF2_med_SSN2, fun_maxF2UT max_F2_ut = 0.0 do i=0,23 max_F2_ut(i) = maxF2_med_SSN2(i,monthut,cov,dlat,dlong) xUT(i) = dble(i) end do t = dble(UT) med_max_F2_SSN2 = fun_maxF2UT(t) return end function med_max_F2_SSN2 c c Median of hmF2 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_SSN2 - hmF2 for a given UT and monthut (real) with SSN2 c real function maxF2_med_SSN2(iUT,monthut,cov,dlat,dlong) implicit none c .. scalar arguments .. integer monthut, iUT real cov real dlat, dlong +,fexp,flin !++ real fexp1, fexp2 !++ c .. local scalars .. real cov1, cov2 c- +, cov_g !++ real a, b, maxF2_1, maxF2_2 double precision teta c .. local arrays .. double precision coeff_month(0:148,0:47) double precision Kf(0:148) real ft_hmF2_1(12), ft_hmF2_2(12) c SSN2 200 km_2018_hmdsat Bins c Jan Feb Mar Apr May Jun data ft_hmF2_1/ 6.6, 5.6, 5.1, 5.1, 14.4, 12.8, * 11.6, 9.9, 9.5, 10.9, 9.2, 7.9/ data ft_hmF2_2/109.3, 98.1,114.3, 84.3, 87.0,114.1, * 112.6,108.3,101.9, 97.3, 94.7, 92.5/ c Jul Aug Sep Okt Nov Dec c c .. local in common .. double precision umr common/const_r/umr c .. function references .. real fun_maxF2 c .. subroutine references .. c read_hmF2_R umr=atan(1.0)*4./180 teta = 90.0-dlat call read_hmF2_R(monthut,coeff_month) Kf = coeff_month(0:148,iUT) maxF2_1 = fun_maxF2(teta,dlong,Kf) Kf = coeff_month(0:148,iUT+24) maxF2_2 = fun_maxF2(teta,dlong,Kf) c cov1 = ft_hmF2_1(monthut) cov2 = ft_hmF2_2(monthut) c a = (maxF2_2 - maxF2_1)/(cov2-cov1) b = maxF2_2- a*cov2 c c- if (cov>210.0) then c- cov_g = 210.0 c- else c- cov_g = cov c- end if c c- maxF2_med_SSN2 = a*cov + b flin= a*cov + b if (cov.lt.cov1) then !++ maxF2_med_SSN2=flin !++ else !++ fexp1 = exp(-(cov1/(cov2-cov1))) !++ fexp2 = exp(-(cov2/(cov2-cov1))) !++ a = (maxF2_2 - maxF2_1)/(fexp2-fexp1) !++ b = maxF2_2- a*fexp2 !++ fexp= a*exp(-(cov/(cov2-cov1))) + b !++ maxF2_med_SSN2=(flin+fexp)/2. !++ endif return end function maxF2_med_SSN2 c c Read median coefficients of hmF2 c c Input: c monthut -(integer) c Output: c . coeff_month(0:148,0:47) - median coefficients of hmF2 (real*8) c subroutine read_hmF2_R(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 10 format('data\hmF2\hmdsat',i2,'.dat') 30 format(6(d12.5)) 100 format('Non file',1x,a) 110 format('Corrupt file',1x,a) end subroutine read_hmF2_R c c real function hmF2_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_RR c max_F2_0 = med_max_F2_RR(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_RR(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_RR(UT,montha,index,lati,longi) quiet_F2 = (max_F2_0+(dayut-15)*(max_F2_p-max_F2_0)/30.) end if hmF2_R_RR = quiet_F2 return end function hmF2_R_RR c c Median of hmF2 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_RR - median of hmF2 for a given UT and monthut (real) c real function med_max_F2_RR(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_RR, fun_maxF2UT max_F2_ut = 0.0 do i=0,23 max_F2_ut(i) = maxF2_med_RR(i,monthut,cov,dlat,dlong) xUT(i) = dble(i) end do t = dble(UT) med_max_F2_RR = fun_maxF2UT(t) return end function med_max_F2_RR c c Median of hmF2 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_RR - hmF2 for a given UT and monthut (real) with SSN2 c real function maxF2_med_RR(iUT,monthut,cov,dlat,dlong) implicit none c .. scalar arguments .. integer monthut, iUT real cov real dlat, dlong c .. local scalars .. real cov1, cov2 c- +, cov_g !++ real a, b, maxF2_1, maxF2_2 double precision teta c .. local arrays .. double precision coeff_month(0:148,0:47) double precision Kf(0:148) real ft_hmF2_1(12), ft_hmF2_2(12) c RR 200 km_2018_hmdsat Bins c Jan Feb Mar Apr May Jun data ft_hmF2_1/ 2.4, 5.8, 4.9, 3.4, 7.3, 15.6, * 13.6, 7.9, 4.9, 5.3, 6.1, 9.1/ Crem data ft_hmF2_2/113.2, 88.3,132.1, 99.9,109.2, 89.2, Crem * 99.0,104.2,112.6, 98.9,119.2,109.5/ data ft_hmf2_2/148.3,156.9,147.4,141.1,140.8,139.1, * 132.4,141.0,157.6,145.1,146.3,151.5/ c Jul Aug Sep Okt Nov Dec c c .. local in common .. double precision umr common/const_r/umr c .. function references .. real fun_maxF2 +,fexp,flin !++ real fexp1, fexp2 !++ c .. subroutine references .. c read_hmF2_R umr=atan(1.0)*4./180 teta = 90.0-dlat call read_hmF2_R(monthut,coeff_month) Kf = coeff_month(0:148,iUT) maxF2_1 = fun_maxF2(teta,dlong,Kf) Kf = coeff_month(0:148,iUT+24) maxF2_2 = fun_maxF2(teta,dlong,Kf) c cov1 = ft_hmF2_1(monthut) cov2 = ft_hmF2_2(monthut) c a = (maxF2_2 - maxF2_1)/(cov2-cov1) b = maxF2_2 - a*cov2 Crem cov_g = 150.0/0.7 Crem if (cov>cov_g) cov = cov_g c- if (cov>210.0) then c- cov_g = 210.0 c- else c- cov_g = cov c- end if Crem maxF2_med_RR = a*cov + b c- maxF2_med_RR = a*cov_g + b flin= a*cov + b if (cov.lt.cov1) then !++ maxF2_med_RR=flin !++ else !++ fexp1 = exp(-(cov1/(cov2-cov1))) !++ fexp2 = exp(-(cov2/(cov2-cov1))) !++ a = (maxF2_2 - maxF2_1)/(fexp2-fexp1) !++ b = maxF2_2- a*fexp2 !++ C- maxF2_med_SSN2h = a*exp(-(cov/(cov2-cov1))) + b !++ fexp= a*exp(-(cov/(cov2-cov1))) + b !++ maxF2_med_RR=(flin+fexp)/2. !++ endif !++ c return end function maxF2_med_RR c c c function fun_maxF2UT c Furie's UT set expansion of median foF2 or hmF2 c . Input parameter: c. t - time (real*8) c . Output function fun_maxF2UT: c . foF2 or hmF2 at any time UT (real) c real function fun_maxF2UT(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_UT, fun_Gk_UT c dtr=atan(1.0)*4.0/12.0 call Koeff_UT(Kf_UT) call fun_Gk_UT(t,Gk_UT) maxF2 = 0.d0 do k=0,6 maxF2 = maxF2 + Kf_UT(k)*Gk_UT(k) end do c fun_maxF2UT = maxF2 return end function fun_maxF2UT c c fun_maxF2UT c Auxiliary function Koeff_UT c subroutine Koeff_UT(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_UT c mm = 3 mk = 2*mm call fun_Akp_UT(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_UT c c fun_maxF2UT c Auxiliar function fun_Akp_UT c subroutine fun_Akp_UT(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_UT, fun_Fk_UT 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_UT(t,Gk_UT) call fun_Fk_UT(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_UT(t,Gk_UT) call fun_Fk_UT(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_UT c c fun_maxF2UT c Auxiliar function fun_Fk_UT c subroutine fun_Fk_UT(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_UT c c fun_maxF2UT c Auxiliar function fun_Gk_UT c subroutine fun_Gk_UT(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_UT c c Spatial interpolation with spherical harmonic functions c medians of foF2 or hmF2 c . Input parameters: c. teta - colatitude (real*8) c. long - longitude (real) c . Output of function fun_maxF2: c. Kf - array of expansion coefficients (real*8) c real function fun_maxF2(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_Gk call fun_Gk(teta,long,Gk) maxF2 = 0.d0 do k=0,148 maxF2 = maxF2 + Kf(k)*Gk(k) end do fun_maxF2 = maxF2 return end function fun_maxF2 c c functions Gk c Spherical harmonic expansion c . Input parameters: c. teta - colatitude (real*8) c. long - longitude (real) c . Output: c. Gk - expansion coefficients array (real*8) c subroutine fun_Gk(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_Gk c c fun_Gk c Auxiliar function c subroutine apoly_Legendre(mm,nn,teta,p) implicit none c .. scalar arguments .. integer mm, nn double precision teta c .. array arguments .. double precision p(0:mm,0:nn) c .. local scalars .. integer j,l,m,n double precision z, s ! .. local in common .. double precision umr common/const_r/umr p = 0.0 z=cos(umr*teta) p(0,0)=1. p(0,1)=z if (mm.ne.0) p(1,1)=sin(umr*teta) do j=2,mm p(j,j)=(2*j-1)*p(j-1,j-1)*p(1,1) end do do m=0,mm do n=1,nn if (m.gt.n) then p(m,n) = 0.0 cycle end if if ((n+1).gt.nn) exit if (n+1.eq.m) cycle if (m.gt.(n-1)) then p(m,n+1)= (2*n+1)*z*p(m,n)/(n+1-m) else p(m,n+1)=((2*n+1)*z*p(m,n)-(n+m)*p(m,n-1))/(n+1-m) end if end do end do do n=1,nn do m=1,mm if (m.gt.n) then p(m,n) = 0.0 exit end if s=1 do l=n-m+1,n+m s=s*l end do p(m,n)=p(m,n)*sqrt(2./s) end do end do return end subroutine apoly_Legendre 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 module Mod_SSN implicit none integer iny_R12 integer iny_R character(256) file_R12 Cold data file_R12/'data/Indices/SSN_R12.dat'/ data file_R12/'data/Indices/ssn2_12.dat'/ ! NEW integer n_R12 data n_R12/212/ character(256) file_R data file_R/'data/Indices/SN_R.dat'/ integer n_R data n_R/214/ end module Mod_SSN c c c c c Calculating the day number of the year c c .. Input: c . day - (integer) c . month - (integer) c . year - (integer) c c . Output: c . nday - day number of the year c subroutine nnday(day,month,year,nday) c .. scalar arguments .. integer day, month, year, nday c .. local scalars .. integer leap c nday=31*(month-1)+day if(month.le.2) return leap=year/4-(year-1)/4 nday=nday-int(0.4*month+2.3)+leap c return end subroutine nnday c c 1 day increase c c .. Input: c . day - (integer) c . month - (integer) c . year - (integer) c c . Output: c . iday - day+1 (integer) c . imonth - (integer) c . iyear - (integer) c subroutine dayadd(day,month,year,iday,imonth,iyear) c .. scalar arguments .. integer day, month, year, iday, imonth, iyear c .. local scalars .. integer daym c .. local arrays .. integer dayinm(12) c .. data dayinm/31,28,31,30,31,30,2*31,30,31,30,31/ c daym = dayinm(month) if (month.eq.2) daym = daym + year/4 - (year-1)/4 iday = day + 1 imonth = month iyear = year if (iday.gt.daym) then iday = 1 imonth = imonth + 1 if (imonth.gt.12) then imonth = 1 iyear = iyear + 1 endif endif c return end subroutine dayadd c c Subtraction 1 day c c .. Input: c . day - äåíü (integer) c . month - ìåñÿö (integer) c . year - ãîä (integer) c c . Output: c . iday - day-1 (integer) c . imonth - (integer) c . iyear - (integer) c subroutine dayres(day,month,year,iday,imonth,iyear) c .. scalar arguments .. integer day, month, year, iday, imonth, iyear c .. local arrays .. integer dayinm(12) c .. data dayinm/31,28,31,30,31,30,2*31,30,31,30,31/ c iday = day - 1 imonth = month iyear = year if (iday.eq.0) then imonth = month - 1 if (imonth.eq.0) then iday = 31 imonth = 12 iyear = year - 1 else iday = dayinm(imonth) if (imonth.eq.2) iday = iday+year/4-(year-1)/4 endif endif c return end subroutine dayres c c Convert UT to LT and vice versa c c .. Input: c . jhour = 0 UT to LT, and 1 LT to UT (integer) c . day - (integer) c . month - (integer) c . year - (integer) c . longi - geographic longitude (real) c UT - (real) or LT - (real) c c . Output: c UT - (real) or LT - (real) c . days - (integer) c . months - (integer) c . years - (integer) c subroutine rconvtime(jhour,day,month,year,longi,LT,UT, * days,months,years) c .. scalar arguments .. integer jhour, day, month, year, days, months, years real longi, LT, UT c .. local scalars .. real delut c .. subroutine references .. c dayres, dayadd, nnday c if (longi.lt.180.0) then delut = longi/15.0 if (jhour.eq.0) then LT = UT + delut if (LT.ge.24.0) then LT = LT - 24.0 call dayadd(day,month,year,days,months,years) else days = day months = month years = year end if else UT = LT - delut if (UT.lt.0.0) then UT = UT + 24.0 call dayres(day,month,year,days,months,years) else days = day months = month years = year end if end if else delut = (360.0-longi)/15.0 if (jhour.eq.0) then LT = UT - delut if (LT.lt.0.0) then LT = LT + 24.0 call dayres(day,month,year,days,months,years) else days = day months = month years = year end if else UT = LT + delut if (UT.ge.24.0) then UT = UT - 24.0 call dayadd(day,month,year,days,months,years) else days = day months = month years = year end if end if end if return end subroutine rconvtime