-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathNRNG.f90
169 lines (154 loc) · 6.86 KB
/
NRNG.f90
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
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
program nrng
implicit none
integer, parameter :: n0 = 24, n1 = 50
real(8) :: tide_height(n0), glacier_size(n1), perturbation, total_sum, c2, y2, t2, carry_over2
real(8) :: mean, sum_sq, std_dev
integer :: i, io_status0, io_status1, count
real(8) :: min_height, max_height, normalized_height, rand_num0, total_tide, c0, y0, t0, carry_over0, c3, y3, t3
real(8) :: min_size, max_size, normalized_size, rand_num1, total_glacier, c1, y1, t1, carry_over1, carry_over3
character(len=100) :: line
character(len=20) :: datetime, year
character(len=20) :: tide_height_str
character(len=20) :: glacier_size_str
! Set a small offset to avoid normalization issues
perturbation = 1.0e-6
total_tide = 0.0
total_glacier = 0.0
total_sum = 0.0
carry_over0 = 0.0
carry_over1 = 0.0
carry_over2 = 0.0
carry_over3 = 0.0
c0 = 0.0
c1 = 0.0
c2 = 0.0
c3 = 0.0
count = 0
! Open the tide data file
open(unit=20, file='mock_tide_data.csv', status='old', action='read')
! Skip the header line
read(20, '(A)')
! Read the tide heights
do i = 1, n0
read(20, '(A)', iostat=io_status0) line
if (io_status0 == 0) then
read(line, '(A20, 1X, A)', iostat=io_status0) datetime, tide_height_str
if (io_status0 == 0) then
read(tide_height_str, *) tide_height(i)
else
print *, "Error parsing line ", i
exit
end if
else
print *, "Error reading line ", i
exit
end if
end do
close(20)
! Find the min and max tide heights
min_height = minval(tide_height)
max_height = maxval(tide_height)
! Open the glacier data file
open(unit=20, file='simulated_glacier_data.csv', status='old', action='read')
! Skip the header line
read(20, '(A)')
! Read the glacier sizes
do i = 1, n1
read(20, '(A)', iostat=io_status1) line
if (io_status1 == 0) then
read(line, '(A20, 1X, A)', iostat=io_status1) year, glacier_size_str
if (io_status1 == 0) then
read(glacier_size_str, *) glacier_size(i)
else
print *, "Error parsing line ", i
exit
end if
else
print *, "Error reading line ", i
exit
end if
end do
close(20)
! Find the min and max glacier sizes
min_size = minval(glacier_size)
max_size = maxval(glacier_size)
! Generate random numbers based on normalized tide heights
do i = 1, 10
call random_number(normalized_height)
normalized_height = ((tide_height(mod(i, n0) + 1) - min_height) / (max_height - min_height)) + perturbation * rand()
rand_num0 = normalized_height * rand()
!total_tide = total_tide + rand_num0
carry_over0 = carry_over0 + rand_num0
y0 = rand_num0 - c0 ! So far, so good: c0 is 0
t0 = total_tide + y0 ! Alas, sum is big, y0 small, so low-order digits of y are lost.
c0 = (t0 - total_tide) - y0 ! (t0-total_tide) recovers the high part of y; subtracting y recovers -(low part of y)
total_tide = t0 ! Algebraically, c0 should always be zero. Beware overly-aggressive optimizing compilers!
! Next time around, the lost low part will be added to y in a fresh attempt.
! Print the current total_tide and the total sum with more decimal places
total_tide = total_tide + carry_over0
carry_over0 = 0.0
! Debugging print statements
!print *, 'Normalized simulated tide height:', normalized_height
!print *, 'TRNG: ', rand_num0
! Print statements with higher precision
!print '(A, F18.15)', 'TRNG: ', rand_num0
! Generate random numbers based on normalized glacier sizes
call random_number(normalized_size)
normalized_size = ((glacier_size(mod(i, n1) + 1) - min_size) / (max_size - min_size)) + perturbation * rand()
rand_num1 = normalized_size * rand()
!total_glacier = total_glacier + rand_num1
carry_over1 = carry_over1 + rand_num1
y1 = rand_num1 - c1 ! So far, so good: c1 is 0
t1 = total_glacier + y1 ! Alas, sum is big, y1 small, so low-order digits of y are lost.
c1 = (t1 - total_glacier) - y1 ! (t1-total_glacier) recovers the high part of y; subtracting y recovers -(low part of y)
total_glacier = t1 ! Algebraically, c1 should always be zero. Beware overly-aggressive optimizing compilers!
! Next time around, the lost low part will be added to y in a fresh attempt.
! Print the current total_tide and the total sum with more decimal places
total_glacier = total_glacier + carry_over1
carry_over1 = 0.0
! Debugging print statements
!print *, 'Simulate Glacier Decay:', normalized_size
!print *, 'GRNG: ', rand_num1
!
print '(A, F24.18)', 'TRNG: ', rand_num0
print '(A, F24.18)', 'GRNG: ', rand_num1
!total_sum = total_tide + total_glacier
!total_sum = total_sum + total_tide
carry_over3 = carry_over3 + total_tide
y3 = total_tide - c3 ! So far, so good: c3 is 0
t3 = total_sum + y3 ! Alas, sum is big, y3 small, so low-order digits of y are lost.
c3 = (t3 - total_sum) - y3 ! (t3-total_sum) recovers the high part of y; subtracting y recovers -(low part of y)
total_sum = t3 ! Algebraically, c should always be zero. Beware overly-aggressive optimizing compilers!
! Next time around, the lost low part will be added to y in a fresh attempt.
! Print the current total_tide and the total sum with more decimal places
total_sum = total_sum + carry_over3
carry_over3 = 0.0
carry_over2 = carry_over2 + total_glacier
y2 = total_glacier - c2 ! So far, so good: c2 is 0
t2 = total_sum + y2 ! Alas, sum is big, y2 small, so low-order digits of y are lost.
c2 = (t2 - total_sum) - y2 ! (t2-total_sum) recovers the high part of y; subtracting y recovers -(low part of y)
total_sum = t2 ! Algebraically, c should always be zero. Beware overly-aggressive optimizing compilers!
! Next time around, the lost low part will be added to y in a fresh attempt.
! Print the current total_tide and the total sum with more decimal places
total_sum = total_sum + carry_over2
carry_over2 = 0.0
count = count + 1
end do
!print '(A, F24.18)', 'Total TRNG: ', total_tide
!print '(A, F24.18)', 'Total GRNG: ', total_glacier
print '(A, F30.20)', 'Sum: ', total_sum
mean = total_sum / count
print '(A, F30.20)', 'Mean:', mean
! Calculate standard deviation (sigma)
sum_sq = 0.0
do i = 1, count
sum_sq = sum_sq + (total_sum - mean)**2
end do
std_dev = sqrt(sum_sq / (count - 1))
print '(A, F30.20)', 'Std :', std_dev
contains
function rand() result(r)
real(8) :: r
call random_number(r)
end function rand
end program nrng