-
Notifications
You must be signed in to change notification settings - Fork 37
Expand file tree
/
Copy pathtest_fang.f90
More file actions
107 lines (89 loc) · 2.35 KB
/
Copy pathtest_fang.f90
File metadata and controls
107 lines (89 loc) · 2.35 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
Program test_fang
!! Need program statement for FORD
!! Reproduces data of:
!! * Figure 3 in https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2008JA013384
!! * Figure 2 in Fang 2010
use, intrinsic:: iso_fortran_env, only: real32, real64
use phys_consts, only: wp
use ionrate, only: ionization_fang2008, ionization_fang2010
implicit none
integer :: i, argc
real(wp) :: Q0_erg, f107, f107a, Ap, glat, glon, UTsec, altrange(3), a
real(wp), allocatable :: alt_km(:), E0_keV(:), Qtot08(:,:), Qtot10(:,:)
integer :: doy
real(wp) :: massden_gcm3, meanmass_g
argc = command_argument_count()
if (argc < 1) then
Q0_erg = 1
else
call parse_arg(1, Q0_erg)
endif
if (argc < 2) then
E0_keV = [0.1, 1., 10., 100., 1000.]
else
allocate(E0_keV(1))
call parse_arg(2, E0_keV(1))
endif
if (argc<5) then
altrange = [20, 400, 2]
else
do i = 3,5
call parse_arg(i, altrange(i-2))
enddo
endif
if (argc<12) then
f107=50
f107a=50
Ap=5
glat=60
glon=0
doy=1
UTsec=0
else
call parse_arg(6, f107)
call parse_arg(7, f107a)
call parse_arg(8, Ap)
call parse_arg(9, glat)
call parse_arg(10, glon)
call parse_arg(11, doy)
call parse_arg(12, UTsec)
endif
i = int((altrange(2)-altrange(1)) / altrange(3)) + 1
allocate(alt_km(i))
alt_km(1) = altrange(1)
do i = 2,size(alt_km)
alt_km(i) = alt_km(i-1) + altrange(3)
enddo
allocate(Qtot08(size(alt_km), size(E0_keV)), Qtot10(size(alt_km), size(E0_keV)))
do i = 1, size(alt_km)
Qtot08(i, :) = ionization_fang2008(Q0_erg, E0_keV, alt_km(i), f107, f107a, Ap, glat, glon, doy, UTsec)
Qtot10(i, :) = ionization_fang2010(Q0_erg, E0_keV, alt_km(i), f107, f107a, Ap, glat, glon, doy, UTsec)
enddo
print '(A,25F10.1)', 'alt[km]/E0[keV]', E0_keV
do i = 1,size(alt_km)
print '(F7.1,5F15.3)', alt_km(i), Qtot08(i,:)
enddo
print '(/,A,25F10.1)', 'alt[km]/Emono[keV]', E0_keV
do i = 1,size(alt_km)
print '(F7.1,5F15.3)', alt_km(i), Qtot10(i,:)
enddo
contains
subroutine parse_arg(iarg, val)
integer, intent(in) :: iarg
class(*), intent(out) :: val
character(80) :: argv
call get_command_argument(iarg, argv)
select type(val)
type is (integer)
read(argv,*) val
type is (real(real64))
read(argv,*) val
type is (real(real32))
read(argv,*) val
type is (character(*))
read(argv,*) val
class default
error stop 'unknown type to parse'
end select
end subroutine parse_arg
end program