-
Notifications
You must be signed in to change notification settings - Fork 7
/
Copy pathhgpt.f90
335 lines (299 loc) · 13.1 KB
/
hgpt.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
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
! hgpt.f90
!
!Mateus, P.; Catalão, J.; Mendes, V.B.; Nico, G. An ERA5-Based Hourly Global Pressure and Temperature (HGPT) Model.
!Remote Sens. 2020, 12, 1098; https://doi.org/10.3390/rs12071098
!
! This routine determines the surface pressure (P), surface air temperature (T),
! weighed mean temperature (Tm), and zenith hydrostatic delay (ZHD) from binary coefficient files
! As available from:
! https://github.com/pjmateus/hgpt_model (release v1.0)
! press_grid.bin; temp_grid.bin; and tm_grid.bin
!
! It is admitted that the binary files with the coefficients are in the same directory as this script
! If you don't want it you can change it in the open statement
!
! The epoch can be an real*8 array of size 1, and in this case is the Modified Julian Date (MJD)
! or can be an integer array of size 6, with the Gregorian Calendar in the following format (year, month, day, hour, min, sec)
!
! All parameters are bilinear interpolated to the input ellipsoidal longitude and latitude
!
! Reference for HGPT:
! An ERA5-based hourly global temperature and pressure (HGTP) model (submitted to Remote Sensing, MDPI)
!
! INPUT:
! dt : if size(dt)=1 => modified julian date
! if size(dt)=6 => year, month, day, hour, min, sec
! ndt : size of dt
! x0 : ellipsoidal longitude (degrees)
! y0 : ellipsoidal latitude (degrees)
! z0 : height (m)
! z0_type : ‘orth’ for orthometric height or ‘elli’ for ellipsoidal height
!
! OUTPUT:
! P : surface pressure valid at (x0, y0, z0), in hPa
! T : surface air temperature valid at (x0, y0, z0), in Kelvins
! Tm : weighed mean temperature valid at (x0, y0, z0), in Kelvins
! ZHD : zenith hydrostatic delay, valid at (x0, y0, z0), in meters
!
! Set the following variables in the script from which you call hgpt.f90:
!--------------------------------------------------------------------------!
! real :: x0, y0, z0, P, T, Tm, ZHD !
! real*8, dimension(6) :: dt ! Gregorian Calendar !
! !
! Define variables values and call hgtp: !
! call hgpt(dt, 6, x0, y0, z0, 'orth', P, T, Tm, ZHD) !
! or !
! real :: x0, y0, z0, P, T, Tm, ZHD !
! real*8, dimension(1) :: dt ! Modified Julian Date (MJD) !
! !
! Define variables values and call hgtp in the same way: !
! call hgpt(dt, 1, x0, y0, z0, 'orth', P, T, Tm, ZHD) !
!--------------------------------------------------------------------------!
!Example: !
!program call_hgpt !
! implicit None !
! real :: x0, y0, z0, P, T, Tm, ZHD !
! real*8, dimension(1) :: dt !
! y0 = 38.5519 !
! x0 = -9.0147 !
! z0 = 25 !
! dt(1) = 58119.5 !
! call hgpt(dt, size(dt), x0, y0, z0, 'orth', P, T, Tm, ZHD) !
! write(*,*) P, T, Tm, ZHD !
!end program call_hgpt !
!--------------------------------------------------------------------------!
! written by Pedro Mateus (2020/01/15)
! Instituto Dom Luiz (IDL), Faculdade de Ciências, Universidade de Lisboa, 1749-016 Lisboa, Portugal
! pjmateus@fc.ul.pt
!
!
subroutine hgpt(dt, ndt, x0, y0, z0, z0_type, P, T, Tm, ZHD)
implicit None
real, intent(in) :: x0, y0, z0
integer, intent(in) :: ndt
real*8, intent(in) :: dt(ndt)
character(len=4), intent(in) :: z0_type
real, intent(out) :: T, P, Tm, ZHD
real*8 :: mjd
integer :: ios, i, j, year, month, day, hour, minute, sec
integer, parameter :: row = 721, col = 1440
real, parameter :: pi = 4*atan(1.0), p1 = 365.25, p2 = 182.6250, p3 = 91.3125, deg2rad = pi/180.0
real, dimension(:) :: lon(col), lat(row)
real, dimension(:,:) :: y_intercept(row, col), slope(row, col), a1(row, col), a2(row, col), a3(row, col), &
orography(row, col), undu(row, col)
integer*2, dimension(:,:) :: f1(row, col), f2(row, col), f3(row, col)
real :: a, b, amp1, pha1, amp2, pha2, amp3, pha3, geo_height, N, H_orth, H_ellip
! Input datetime format
if ( ndt == 6) then
year = idint(dt(1)); month = idint(dt(2)); day = idint(dt(3))
hour = idint(dt(4)); minute = idint(dt(5)); sec = idint(dt(6))
call greg2modjul(year, month, day, hour, minute, sec, mjd)
else if ( ndt == 1 ) then
mjd = dt(1)
call modjul2greg(mjd, year, month, day, hour, minute, sec)
else
stop 'Use 1) Modified Julian Date (MJD) or 2) Gregorian date (y,m,d,HH,MM,SS).'
end if
! Geographic latitude and longitude (equal to ERA5 grid)
lon = (/(( -180+0.25*i), i=1, col)/)
lat = (/((-90.25+0.25*i), i=1, row)/)
! Locate i, j indices
call locate_ij(row, lat, y0, i)
call locate_ij(col, lon, x0, j)
! Open and read the surface air temperature coefficients file
open(10, file='temp_grid.bin', status='old', action='read', &
form='unformatted', access='stream', iostat=ios)
if ( ios /= 0 ) stop 'Error opening file temp_grid.bin'
! hour varies from 0 to 23 UTC
call fseek(10, (row*col*26)*hour, 0)
read(10) y_intercept
read(10) slope
read(10) a1
read(10) f1
read(10) a2
read(10) f2
read(10) a3
read(10) f3
close(10)
! Bilinear interpolation
call interpolate(lon, lat, y_intercept, x0, y0, i, j, a)
call interpolate(lon, lat, slope, x0, y0, i, j, b)
call interpolate(lon, lat, a1, x0, y0, i, j, amp1)
call interpolate(lon, lat, float(f1)/10000.0, x0, y0, i, j, pha1)
call interpolate(lon, lat, a2, x0, y0, i, j, amp2)
call interpolate(lon, lat, float(f2)/10000.0, x0, y0, i, j, pha2)
call interpolate(lon, lat, a3, x0, y0, i, j, amp3)
call interpolate(lon, lat, float(f3)/10000.0, x0, y0, i, j, pha3)
! Surface air temperature model
T = a + b*(mjd - 51178) + amp1*cos(2*pi*(mjd - 51178)/p1+pha1) + &
amp2*cos(2*pi*(mjd - 51178)/p2+pha2) + &
amp3*cos(2*pi*(mjd - 51178)/p3+pha3)
! Open and read the surface pressure coefficients file
open(11, file='press_grid.bin', status='old', action='read', &
form='unformatted', access='stream', iostat=ios)
if ( ios /= 0 ) stop 'Error opening file press_grid.bin'
! hour varies from 0 to 23 UTC
call fseek(11, (row*col*20)*hour, 0)
read(11) y_intercept
read(11) slope
read(11) a1
read(11) f1
read(11) a2
read(11) f2
close(11)
! Bilinear interpolation
call interpolate(lon, lat, y_intercept, x0, y0, i, j, a)
call interpolate(lon, lat, slope, x0, y0, i, j, b)
call interpolate(lon, lat, a1, x0, y0, i, j, amp1)
call interpolate(lon, lat, float(f1)/10000.0, x0, y0, i, j, pha1)
call interpolate(lon, lat, a2, x0, y0, i, j, amp2)
call interpolate(lon, lat, float(f2)/10000.0, x0, y0, i, j, pha2)
! Surface pressure model
P = a + b*(mjd - 51178) + amp1*cos(2*pi*(mjd - 51178)/p1+pha1) + &
amp2*cos(2*pi*(mjd - 51178)/p2+pha2)
! Open and read the weight mean temperature and undulation coefficients file
open(12, file='tm_grid.bin', status='old', action='read', &
form='unformatted', access='stream', iostat=ios)
if ( ios /= 0 ) stop 'Error opening file tm_grid.bin'
read(12) y_intercept
read(12) slope
read(12) orography
read(12) undu
close(12)
! Bilinear interpolation
call interpolate(lon, lat, y_intercept, x0, y0, i, j, a)
call interpolate(lon, lat, slope, x0, y0, i, j, b)
call interpolate(lon, lat, orography, x0, y0, i, j, geo_height)
call interpolate(lon, lat, undu, x0, y0, i, j, N)
! Zenith hydrostatic delay (ZHD), Saastamoinen model
if ( z0_type == 'orth' ) then
H_orth = z0
else if ( z0_type == 'elli' ) then
H_orth = z0 - N
else
stop 'Use 1) <<orth>> for Orthometric height or 2) <<elli>> for Ellipsoidal height (in m).'
end if
! Correction to P and T (see Guochanf Xu, GPS Theory, Algorithms and Applications, 2nd Edition, page 56)
P = (P*100.0 * (1.0 - 0.0065/T * (H_orth - geo_height))**5.2559)/100.0
T = T - 0.0065*(H_orth - geo_height)
! Weight mean temperature (Tm) linear model
Tm = a + b * T
! ZHD using the Saastamoinen Model (see Saastamoinen, 1973)
ZHD = (0.0022768 * P)/(1 - 0.0026*cos(2*deg2rad*y0)-0.00000028*H_orth)
return
end subroutine hgpt
subroutine locate_ij(length, array, value, idx)
! Given an array and a value, returns the index of the element that
! is closest to, but less than, the given value.
! Uses a binary search algorithm.
implicit none
integer, intent(in) :: length
real, dimension(length), intent(in) :: array
real, intent(in) :: value
integer, intent(out) :: idx
integer :: left, middle, right
left = 1
right = length
do
if (left > right) then
exit
end if
middle = nint((left+right) / 2.0)
if ( abs(array(middle) - value) <= 1e-9) then
idx = middle
return
else if (array(middle) > value) then
right = middle - 1
else
left = middle + 1
end if
end do
idx = right
return
end subroutine locate_ij
subroutine interpolate(lon, lat, f, x0, y0, i, j, value)
! This function uses bilinear interpolation to estimate the value
! of a function f at point (x0,y0) using the point grid (i, j), (i+1,j), (i,j+1)
! and (i+1,j+1), see locate_ij subroutine
! In the limits (borders) the nearest-neighbor interpolation is used
! f is assumed to be sampled on a regular grid, with the grid x values specified
! by lon array and the grid y values specified by lat array
implicit none
integer, parameter :: col = 1440, row = 721
real, dimension(:), intent(in) :: lon(col), lat(row)
real, dimension(row, col), intent(in) :: f
integer, intent(in) :: i, j
real, intent(in) :: x0, y0
integer :: ix, jy
real, intent(out) :: value
real :: x1, x2, y1, y2, xy1, xy2
ix = i
jy = j
if ( ix == 0 .or. ix >= row .or. jy == 0 .or. jy >= col ) then
! Nearest-neighbor interpolation
if ( ix == 0 ) ix = ix + 1
if ( ix >= row ) ix = row
if ( jy == 0 ) jy = jy + 1
if ( jy >= col ) jy = col
value = f(ix, jy)
else
! Bilinear interpolation
x1 = lon( jy )
x2 = lon( jy + 1 )
y1 = lat( ix )
y2 = lat( ix + 1 )
xy1 = (x2 - x0)/(x2 - x1)*f(ix, jy) + (x0 - x1)/(x2 - x1)*f(ix, jy+1)
xy2 = (x2 - x0)/(x2 - x1)*f(ix+1, jy) + (x0 - x1)/(x2 - x1)*f(ix+1, jy+1)
value = (y2 - y0)/(y2 - y1)*xy1 + (y0 - y1)/(y2 - y1)*xy2
end if
return
end subroutine interpolate
subroutine greg2modjul(year, month, day, hour, minute, sec, mjd)
! greg2modjul(..., mjd) returns the
! Modified Julian Date (MJD) corresponding to the
! Gregorian calendar date (year, month, day, hour, minute, and second).
implicit none
integer, intent(in) :: year, month, day, hour, minute, sec
real*8, intent(out) :: mjd
real*8 :: yyyy, mm, jd
! January & February
if ( month <= 2 ) then
yyyy = dfloat(year) - 1_8
mm = dfloat(month) + 12_8
else
yyyy = dfloat(year)
mm = dfloat(month)
end if
jd = floor( 365.25_8*(yyyy + 4716.0_8)) + floor( 30.6001_8*( mm + 1.0_8)) + 2.0_8 - &
floor( yyyy/100.0_8 ) + floor( floor( yyyy/100.0_8 )/4.0_8 ) + dfloat(day) - 1524.5_8 + &
(dfloat(hour) + dfloat(minute)/60_8 + dfloat(sec)/3600_8)/24_8
mjd = jd - 2400000.5_8
return
end subroutine greg2modjul
subroutine modjul2greg(mjd, year, month, day, hour, minute, sec)
! modjul2greg(mjd, ...) returns the
! Gregorian calendar date (year, month, day, hour, minute, and second)
! corresponding to the Modified Julian Date (MJD).
implicit none
real*8, intent(in) :: mjd
real*8 :: jd, ijd, fjd, a, b, c, d, e, m
integer, intent(out):: year, month, day, hour, minute, sec
jd = mjd + 2400000.5_8
ijd = floor(jd + 0.5_8)
fjd = jd - ijd + 0.5_8
sec = 86400_8 * fjd
hour = ifix( sec / 3600.0 )
sec = sec - 3600_8 * hour
minute = ifix( sec / 60.0 )
sec = sec - 60_8 * minute
a = ijd + 32044_8
b = floor((4_8 * a + 3_8) / 146097_8)
c = a - floor((b * 146097_8) / 4_8)
d = floor((4_8 * c + 3_8) / 1461_8)
e = c - floor((1461_8 * d) / 4_8)
m = floor((5_8 * e + 2_8) / 153_8)
day = e - floor((153_8 * m + 2_8) / 5_8) + 1_8
month = m + 3_8 - 12_8 * floor(m / 10_8)
year = b * 100_8 + d - 4800_8 + floor(m / 10_8)
return
end subroutine modjul2greg