-
Notifications
You must be signed in to change notification settings - Fork 1
/
NORTRIP_functions.f90
463 lines (360 loc) · 15.1 KB
/
NORTRIP_functions.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
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
!----------------------------------------------------------------------
! Various routines and functions
!----------------------------------------------------------------------
! Calculates running mean temperature
!----------------------------------------------------------------------
subroutine NORTRIP_running_mean_temperature(num_running_hours)
use NORTRIP_definitions
implicit none
!Avergaing time is 1.5 - 3 days
!road_meteo_data(T_sub_index,:,tr,ro) has already been initialised as meteo_data(T_a_index,:,ro) in NORTRIP_initialise_data
real num_running_hours
integer min_running_index
logical :: use_running_mean=.false.
!do ro=n_roads_start,n_roads_end
do tr=1,num_track
if (available_meteo_data(T_sub_input_index)) then
road_meteo_data(T_sub_index,:,tr,ro)=meteo_data(T_sub_input_index,:,ro)
else
do ti=min_time,max_time
!specify the minimum index
min_running_index=max(1,ti-int(num_running_hours))
!write(*,*)min_running_index,tr,ro
if (use_running_mean) then
!Calculate running mean
road_meteo_data(T_sub_index,ti,tr,ro)=sum(meteo_data(T_a_index,min_running_index:ti,ro))/size(meteo_data(T_a_index,min_running_index:ti,ro),1)
else
!Alternative formulation so as to preserve the initial value in min_time index. Assumes a number in min_time already
if (ti.gt.min_time) then
road_meteo_data(T_sub_index,ti,tr,ro)= &
road_meteo_data(T_sub_index,max(1,ti-1),tr,ro)*(1.-dt/num_running_hours) &
+meteo_data(T_a_index,ti,ro)*dt/num_running_hours
!write(*,*) ti,road_meteo_data(T_sub_index,ti,tr,ro),road_meteo_data(T_sub_index,max(1,ti-1),tr,ro),meteo_data(T_a_index,ti,ro)
endif
endif
!write(*,*) ti,road_meteo_data(T_sub_index,ti,tr,ro)
end do
endif
!Treat this differently in a tunnel
if (roadtype_index(ro).eq.tunnel_roadtype.or.roadtype_index(ro).eq.tunnelportal_roadtype) then
!Not certain what to do here as not enough information is available about a tunnels temperature. Set at ambient temperature
!Do nothing for the time being
!road_meteo_data(T_sub_index,:,tr,ro)=meteo_data(T_a_index,:,ro)
endif
!Treat this differently on a bridge. Set sublayer temperature to atmospheric layer
if (roadtype_index(ro).eq.bridge_roadtype) then
road_meteo_data(T_sub_index,:,tr,ro)=meteo_data(T_a_index,:,ro)
endif
end do
if (((ro.eq.1.or.ro.eq.n_roads).and..not.use_single_road_loop_flag).or.((ro_tot.eq.1.or.ro_tot.eq.n_roads_total).and.use_single_road_loop_flag)) then
write(unit_logfile,'(A)') ''
write(unit_logfile,'(A)') 'Calculating subsurface temperature (NORTRIP_running_mean_temperature)'
write(unit_logfile,'(A)') '================================================================'
write(unit_logfile,'(A40,A3,L)') trim(meteo_match_str(T_sub_input_index))//' available',' = ',available_meteo_data(T_sub_input_index)
write(unit_logfile,'(A)') '----------------------------------------------------------------'
write(unit_logfile,'(a32,a14,a14,a14)') 'Temperature parameter','Min value','Max value','Mean value'
write(unit_logfile,'(A)') '----------------------------------------------------------------'
write(unit_logfile,'(a32,f14.2,f14.2,f14.2)') 'T_a',minval(meteo_data(T_a_index,min_time:max_time,ro)),maxval(meteo_data(T_a_index,min_time:max_time,ro)),sum(meteo_data(T_a_index,min_time:max_time,ro)/(max_time-min_time+1))
write(unit_logfile,'(a32,f14.2,f14.2,f14.2)') 'T_sub',minval(road_meteo_data(T_sub_index,min_time:max_time,1,ro)),maxval(road_meteo_data(T_sub_index,min_time:max_time,1,ro)),sum(road_meteo_data(T_sub_index,min_time:max_time,1,ro)/(max_time-min_time+1))
write(unit_logfile,'(A)') '----------------------------------------------------------------'
endif
!(*,*) road_meteo_data(T_sub_index,min_time:max_time,1,ro)
!write(*,*) meteo_data(T_a_index,min_time:max_time,ro)
!write(*,*) num_running_hours
!stop
!enddo
end subroutine NORTRIP_running_mean_temperature
!----------------------------------------------------------------------
!----------------------------------------------------------------------
function e_sat_func(TC,P)
implicit none
!TC: Degrees C
!P: mbar or hPa
!esat: hPa
!qsat: kg/kg
real TC,P
real e_sat_func
real a,b,c
parameter (a=6.1121,b=17.67,c=243.5)
e_sat_func=a*exp(b*TC/(c+TC))
!qsat=0.622*esat./(P-0.378*esat)
!d_esat_dT=esat*b*c./(TC+c).^2
!d_qsat_dT=0.622.*P./(P-0.378*esat).^2.*d_esat_dT
end function e_sat_func
!----------------------------------------------------------------------
!----------------------------------------------------------------------
function q_sat_func(TC,P)
implicit none
!TC: Degrees C
!P: mbar or hPa
!esat: hPa
!qsat: kg/kg
real TC,P
real esat
real q_sat_func
real a,b,c
parameter (a=6.1121,b=17.67,c=243.5)
esat=a*exp(b*TC/(c+TC))
q_sat_func=0.622*esat/(P-0.378*esat)
!d_esat_dT=esat*b*c./(TC+c).^2
!d_qsat_dT=0.622.*P./(P-0.378*esat).^2.*d_esat_dT
end function q_sat_func
!----------------------------------------------------------------------
!----------------------------------------------------------------------
function q_sat_ice_func(TC,P)
implicit none
!TC: Degrees C
!P: mbar or hPa
!esat: hPa
!qsat: kg/kg
real TC,P
real esat
real q_sat_ice_func
real a,b,c
parameter (a=6.1121,b=22.46,c=272.62)
esat=a*exp(b*TC/(c+TC))
q_sat_ice_func=0.622*esat/(P-0.378*esat)
!d_esat_dT=esat*b*c./(TC+c).^2
!d_qsat_dT=0.622.*P./(P-0.378*esat).^2.*d_esat_dT
end function q_sat_ice_func
!----------------------------------------------------------------------
!----------------------------------------------------------------------
function dewpoint_from_RH_func(TC,RH)
implicit none
!TC: Degrees C
!esat: hPa
!RH: !
real dewpoint_from_RH_func
real TC,RH
real esat,eair
real a,b,c
parameter (a=6.1121,b=17.67,c=243.5)
esat=a*exp(b*TC/(c+TC))
!+.01 to avoid a NaN error in the log when RH=0
eair=(RH+.01)/100.*esat
dewpoint_from_RH_func=c*log(eair/a)/(b-log(eair/a))
end function dewpoint_from_RH_func
!----------------------------------------------------------------------
!----------------------------------------------------------------------
function RH_from_dewpoint_func(TC,dewpoint)
implicit none
!TC: Degrees C
!dewpoint: Degrees C
!esat: hPa
!RH: !
real RH_from_dewpoint_func
real TC,dewpoint
real esat,eair
real a,b,c
parameter (a=6.1121,b=17.67,c=243.5)
esat=a*exp(b*TC/(c+TC))
!The original from dewpoint_from_RH_func
!+.01 to avoid a NaN error in the log when RH=0
!eair=(RH+.01)/100.*esat
!dewpoint_from_RH_func=c*log(eair/a)/(b-log(eair/a))
!process of inversion
!dewpoint*(b-log(eair/a))=c*log(eair/a)
!dewpoint*b-dewpoint*log(eair/a)=c*log(eair/a)
!dewpoint*b=dewpoint*log(eair/a)+c*log(eair/a)
!dewpoint*b=(dewpoint+c)*log(eair/a)
!dewpoint*b=(dewpoint+c)*log(eair/a)
!log(eair/a)=dewpoint*b/(dewpoint+c)
eair=a*exp(dewpoint*b/(dewpoint+c))
RH_from_dewpoint_func=eair/esat*100
end function RH_from_dewpoint_func
!----------------------------------------------------------------------
!----------------------------------------------------------------------
function r_aero_func(FF,z_FF,z_T,z0,z0t,V_veh,N_v,num_veh,a_traffic)
implicit none
real r_aero_func
real FF,z_FF,z_T,z0,z0t
real V_veh(num_veh),N_v(num_veh),a_traffic(num_veh)
integer num_veh
real inv_r_wind,inv_r_traffic,inv_r_aero
integer v
real kappa
parameter (kappa=0.4)
inv_r_wind=max(FF,0.2)*kappa**2/(log(z_FF/z0)*log(z_T/z0t))
inv_r_traffic=0
do v=1,num_veh
inv_r_traffic=inv_r_traffic+N_v(v)*V_veh(v)*a_traffic(v)
enddo
inv_r_traffic=max(1e-6,inv_r_traffic/3600./3.6)
inv_r_aero=inv_r_traffic+inv_r_wind
r_aero_func=1./inv_r_aero
end function r_aero_func
!----------------------------------------------------------------------
!----------------------------------------------------------------------
function r_aero_func_with_stability(FF,Tc,Ts,z_FF,z_T,z0,z0t,V_veh,N_v,num_veh,a_traffic)
implicit none
real r_aero_func_with_stability
real FF,Tc,Ts,z_FF,z_T,z0,z0t
real V_veh(num_veh),N_v(num_veh),a_traffic(num_veh)
integer num_veh
real inv_r_wind,inv_r_traffic,inv_r_aero
real phim,phih,phi_m,phi_h,FF_temp,Rib,eps
integer v,iterations,i
real kappa,g,T0K,a,b,p,q,pi
parameter (kappa=0.4,g=9.8,T0K=273.15,a=16.,b=5.,p=-0.25,q=-0.5,pi=3.1415926)
iterations=2
phi_m=0.0
phi_h=0.0
do i=1,iterations
!Calculate wind at temperature height
FF_temp=max(0.2,FF*(log(z_T/z0)-phi_m)/(log(z_FF/z0)-phi_m))
!Set bilk Richardsons number
Rib=g/(Tc+T0K)*z_T*(Tc-Ts)/FF_temp/FF_temp
!Calculate z/L
eps=Rib*(log(z_T/z0)-phi_m)*(log(z_T/z0)-phi_m)/((log(z_T/z0t)-phi_h))
if (eps.ge.0) then
phim=1+b*eps
phih=1+b*eps
phi_m=-b*eps
phi_h=-b*eps
else
phim=exp(p*log((1.-a*eps)))
phih=exp(q*log((1.-a*eps)))
phi_m=2.*log((1.+1./phim)/2.)+log((1.+1./(phim*phim))/2.)-2.*atan(1./phim)+pi/2.
phi_h=2.*log((1.+1./phih)/2.)
endif
enddo
!Calculate bulk exchange coefficient including wind (inverse of aerodynamic resistance)
inv_r_wind=FF_temp*kappa*kappa/((log(z_T/z0)-phi_m)*(log(z_T/z0t)-phi_h));
inv_r_traffic=0
do v=1,num_veh
inv_r_traffic=inv_r_traffic+N_v(v)*V_veh(v)*a_traffic(v)
enddo
inv_r_traffic=max(1e-6,inv_r_traffic/3600./3.6)
inv_r_aero=inv_r_traffic+inv_r_wind
r_aero_func_with_stability=1./inv_r_aero
end function r_aero_func_with_stability
!----------------------------------------------------------------------
!----------------------------------------------------------------------
function f_spray_func(R_0_spray,V_veh,V_ref_spray,V_thresh_spray,a_spray,water_spray_flag)
implicit none
real f_spray_func
real R_0_spray,V_veh,V_ref_spray,V_thresh_spray,a_spray
integer water_spray_flag
f_spray_func=0
if (water_spray_flag.gt.0.and.V_ref_spray.gt.V_thresh_spray) then
f_spray_func=R_0_spray*(max(0.,(V_veh-V_thresh_spray))/(V_ref_spray-V_thresh_spray))**a_spray
endif
end function f_spray_func
!----------------------------------------------------------------------
!----------------------------------------------------------------------
function mass_balance_func(M_0,P,R,dt)
!mass_balance_func: Caclulates temporal mass balance changes
implicit none
real mass_balance_func
real M_0,P,R,dt
if (P.lt.R*1.0E8) then
mass_balance_func=P/R*(1-exp(-R*dt))+M_0*exp(-R*dt)
else
mass_balance_func=M_0+P*dt
endif
end function mass_balance_func
!----------------------------------------------------------------------
!----------------------------------------------------------------------
function W_func(W_0,h_pave_in,h_dc_in,V_veh,a_wear,s_road,s_roadwear_thresh,s,road_index,tyre_index,brake_index)
! W_func: Wear function
! Depends on:
implicit none
!Output variable
real W_func
!Input variables
real h_pave_in,h_dc_in
real W_0,V_veh,a_wear(5),s_road,s_roadwear_thresh
integer s,road_index,tyre_index,brake_index
!Internal variables
real f_V,f_snow
real h_pave,h_dc
h_pave=h_pave_in
h_dc=h_dc_in
!No wear production due to snow on the surface
f_snow=1.
if (s_road.gt.s_roadwear_thresh) then
f_snow=0.
endif
f_V=max(0.,a_wear(1)+a_wear(2)*(max(V_veh,a_wear(5))/a_wear(4))**a_wear(3))
if (s.eq.road_index) then
h_dc=1.
endif
if (s.eq.tyre_index) then
h_dc=1.
h_pave=1.
endif
if (s.eq.brake_index) then
h_pave=1.
f_snow=1.
endif
W_func=W_0*h_pave*h_dc*f_V*f_snow
end function W_func
!----------------------------------------------------------------------
!----------------------------------------------------------------------
function f_abrasion_func(f_sandpaper_0,h_pave,V_veh,s_road,V_ref,s_roadwear_thresh)
!W_func: Abrasion function
implicit none
real f_abrasion_func
real f_sandpaper_0,h_pave,V_veh,s_road,V_ref,s_roadwear_thresh
real f_V,f_snow
if (V_ref.eq.0) then
f_V=1.
else
f_V=(V_veh/V_ref)
endif
!No wear production due to snow on the surface
f_snow=1.
if (s_road.gt.s_roadwear_thresh) then
f_snow=0.
endif
f_abrasion_func=f_sandpaper_0*h_pave*f_V*f_snow
end function f_abrasion_func
!----------------------------------------------------------------------
!----------------------------------------------------------------------
function f_crushing_func(f_crushing_0,V_veh,s_road,V_ref,s_roadwear_thresh)
!W_func: Sandpaper function
implicit none
real f_crushing_func
real f_crushing_0,V_veh,s_road,V_ref,s_roadwear_thresh
real f_V,f_snow
if (V_ref.eq.0) then
f_V=1.
else
f_V=(V_veh/V_ref)
endif
!No wear production due to snow on the surface
f_snow=1.
if (s_road.gt.s_roadwear_thresh) then
f_snow=0.
endif
f_crushing_func=f_crushing_0*f_V*f_snow
end function f_crushing_func
!----------------------------------------------------------------------
!----------------------------------------------------------------------
function f_susroad_func(V_veh,a_sus)
!f_sus_func: Vehicle speed dependence function for suspension
implicit none
real f_susroad_func
integer num_size
real V_veh
real a_sus(5)
real h_V
h_V=max(0.,a_sus(1)+a_sus(2)*(max(V_veh,a_sus(5))/a_sus(4))**a_sus(3))
f_susroad_func=h_V
end function f_susroad_func
!----------------------------------------------------------------------
!----------------------------------------------------------------------
function R_0_wind_func(FF,tau_wind,FF_thresh)
!R_0_wind_func: Wind blown dust wind speed depedency
implicit none
real R_0_wind_func
real FF,tau_wind,FF_thresh
real h_FF
if (FF.gt.FF_thresh) then
h_FF=(FF/FF_thresh-1.)**3
else
h_FF=0.
endif
R_0_wind_func=1./tau_wind*h_FF
end function R_0_wind_func
!----------------------------------------------------------------------