-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathdust_common.F90
253 lines (221 loc) · 11.3 KB
/
dust_common.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
!=============================================================================
! Common dust module
!=============================================================================
module dust_common
use shr_kind_mod, only: r8 => shr_kind_r8, cl => shr_kind_cl
use cam_abortutils, only: endrun
use cam_logfile, only: iulog
implicit none
private
public :: dust_set_params
contains
!=============================================================================
!
! !DESCRIPTION:
!
! Compute source efficiency factor from topography
! Initialize other variables used in subroutine Dust:
! ovr_src_snk_mss(m,n) and tmp1.
! Define particle diameter and density needed by atm model
! as well as by dry dep model
! Source: Paul Ginoux (for source efficiency factor)
! Modifications by C. Zender and later by S. Levis
! Rest of subroutine from C. Zender's dust model
!=============================================================================
subroutine dust_set_params( nbin, dmt_grd, dmt_vwr, stk_crc )
!
! !USES
!
use physconst, only: pi,rair, gravit
use mo_constants, only: dust_density
use infnan, only: nan, assignment(=)
!
! !ARGUMENTS:
!
integer, intent(in) :: nbin
real(r8),intent(in) :: dmt_grd(:)
real(r8),intent(out) :: dmt_vwr(:)
real(r8),intent(out) :: stk_crc(:)
!
! !REVISION HISTORY
! Created by Samual Levis
! Revised for CAM by Natalie Mahowald
!EOP
!------------------------------------------------------------------------
!------------------------------------------------------------------------
!Local Variables
integer, parameter:: dst_src_nbr =3
integer, parameter:: sz_nbr =200
integer :: m,n !indices
real(r8) :: dmt_min(nbin) ![m] Size grid minimum
real(r8) :: dmt_max(nbin) ![m] Size grid maximum
real(r8) :: dmt_ctr(nbin) ![m] Diameter at bin center
real(r8) :: dmt_dlt(nbin) ![m] Width of size bin
real(r8) :: slp_crc(nbin) ![frc] Slip correction factor
real(r8) :: vlm_rsl(nbin) ![m3 m-3] Volume concentration resolved
real(r8) :: vlc_stk(nbin) ![m s-1] Stokes settling velocity
real(r8) :: vlc_grv(nbin) ![m s-1] Settling velocity
real(r8) :: ryn_nbr_grv(nbin) ![frc] Reynolds number at terminal velocity
real(r8) :: cff_drg_grv(nbin) ![frc] Drag coefficient at terminal velocity
real(r8) :: tmp !temporary
real(r8) :: ln_gsd ![frc] ln(gsd)
real(r8) :: gsd_anl ![frc] Geometric standard deviation
real(r8) :: dmt_vma ![m] Mass median diameter analytic She84 p.75 Tabl.1
real(r8) :: dmt_nma ![m] Number median particle diameter
real(r8) :: lgn_dst !Lognormal distribution at sz_ctr
real(r8) :: eps_max ![frc] Relative accuracy for convergence
real(r8) :: eps_crr ![frc] Current relative accuracy
real(r8) :: itr_idx ![idx] Counting index
real(r8) :: dns_mdp ![kg m-3] Midlayer density
real(r8) :: mfp_atm ![m] Mean free path of air
real(r8) :: vsc_dyn_atm ![kg m-1 s-1] Dynamic viscosity of air
real(r8) :: vsc_knm_atm ![kg m-1 s-1] Kinematic viscosity of air
real(r8) :: vlc_grv_old ![m s-1] Previous gravitational settling velocity
real(r8) :: series_ratio !Factor for logarithmic grid
real(r8) :: lngsdsqrttwopi_rcp !Factor in lognormal distribution
real(r8) :: sz_min(sz_nbr) ![m] Size Bin minima
real(r8) :: sz_max(sz_nbr) ![m] Size Bin maxima
real(r8) :: sz_ctr(sz_nbr) ![m] Size Bin centers
real(r8) :: sz_dlt(sz_nbr) ![m] Size Bin widths
stk_crc(:) = nan
dmt_vwr(:) = nan
! Introducing particle diameter. Needed by atm model and by dry dep model.
! Taken from Charlie Zender's subroutines dst_psd_ini, dst_sz_rsl,
! grd_mk (dstpsd.F90) and subroutine lgn_evl (psdlgn.F90)
! Charlie allows logarithmic or linear option for size distribution
! however, he hardwires the distribution to logarithmic in his code
! therefore, I take his logarithmic code only
! furthermore, if dst_nbr == 4, he overrides the automatic grid calculation
! he currently works with dst_nbr = 4, so I only take the relevant code
! if dust_number ever becomes different from 4, must add call grd_mk (dstpsd.F90)
! as done in subroutine dst_psd_ini
! note that here dust_number = dst_nbr
! Override automatic grid with preset grid if available
do n = 1, nbin
dmt_min(n) = dmt_grd(n) ![m] Max diameter in bin
dmt_max(n) = dmt_grd(n+1) ![m] Min diameter in bin
dmt_ctr(n) = 0.5_r8 * (dmt_min(n)+dmt_max(n)) ![m] Diameter at bin ctr
dmt_dlt(n) = dmt_max(n)-dmt_min(n) ![m] Width of size bin
end do
! sets dust_dmt_vwr ....
! Bin physical properties
gsd_anl = 2.0_r8 ! [frc] Geometric std dev PaG77 p. 2080 Table1
ln_gsd = log(gsd_anl)
! Set a fundamental statistic for each bin
dmt_vma = 2.524e-6_r8 ! [m] Mass median diameter analytic She84 p.75 Table1
dmt_vma = 3.5e-6_r8
! Compute analytic size statistics
! Convert mass median diameter to number median diameter (call vma2nma)
dmt_nma = dmt_vma * exp(-3.0_r8*ln_gsd*ln_gsd) ! [m]
! Compute resolved size statistics for each size distribution
! In C. Zender's code call dst_sz_rsl
do n = 1, nbin
series_ratio = (dmt_max(n)/dmt_min(n))**(1.0_r8/sz_nbr)
sz_min(1) = dmt_min(n)
do m = 2, sz_nbr ! Loop starts at 2
sz_min(m) = sz_min(m-1) * series_ratio
end do
! Derived grid values
do m = 1, sz_nbr-1 ! Loop ends at sz_nbr-1
sz_max(m) = sz_min(m+1) ! [m]
end do
sz_max(sz_nbr) = dmt_max(n) ! [m]
! Final derived grid values
do m = 1, sz_nbr
sz_ctr(m) = 0.5_r8 * (sz_min(m)+sz_max(m))
sz_dlt(m) = sz_max(m)-sz_min(m)
end do
lngsdsqrttwopi_rcp = 1.0_r8 / (ln_gsd*sqrt(2.0_r8*pi))
dmt_vwr(n) = 0.0_r8 ! [m] Mass wgted diameter resolved
vlm_rsl(n) = 0.0_r8 ! [m3 m-3] Volume concentration resolved
do m = 1, sz_nbr
! Evaluate lognormal distribution for these sizes (call lgn_evl)
tmp = log(sz_ctr(m)/dmt_nma) / ln_gsd
lgn_dst = lngsdsqrttwopi_rcp * exp(-0.5_r8*tmp*tmp) / sz_ctr(m)
! Integrate moments of size distribution
dmt_vwr(n) = dmt_vwr(n) + sz_ctr(m) * &
pi / 6.0_r8 * (sz_ctr(m)**3.0_r8) * & ![m3] Volume
lgn_dst * sz_dlt(m) ![# m-3] Number concentrn
vlm_rsl(n) = vlm_rsl(n) + &
pi / 6.0_r8 * (sz_ctr(m)**3.0_r8) * & ![m3] Volume
lgn_dst * sz_dlt(m) ![# m-3] Number concentrn
end do
dmt_vwr(n) = dmt_vwr(n) / vlm_rsl(n) ![m] Mass weighted diameter resolved
end do
! sets stk_crc ...
! calculate correction to Stokes' settling velocity (subroutine stk_crc_get)
eps_max = 1.0e-4_r8
dns_mdp = 100000._r8 / (295.0_r8*rair) ![kg m-3] const prs_mdp & tpt_vrt
! Size-independent thermokinetic properties
vsc_dyn_atm = 1.72e-5_r8 * ((295.0_r8/273.0_r8)**1.5_r8) * 393.0_r8 / &
(295.0_r8+120.0_r8) ![kg m-1 s-1] RoY94 p.102 tpt_mdp=295.0
mfp_atm = 2.0_r8 * vsc_dyn_atm / & !SeP97 p. 455 constant prs_mdp, tpt_mdp
(100000._r8*sqrt(8.0_r8/(pi*rair*295.0_r8)))
vsc_knm_atm = vsc_dyn_atm / dns_mdp ![m2 s-1] Kinematic viscosity of air
do m = 1, nbin
slp_crc(m) = 1.0_r8 + 2.0_r8 * mfp_atm * &
(1.257_r8+0.4_r8*exp(-1.1_r8*dmt_vwr(m)/(2.0_r8*mfp_atm))) / &
dmt_vwr(m) ! [frc] Slip correction factor SeP97 p.464
vlc_stk(m) = (1.0_r8/18.0_r8) * dmt_vwr(m) * dmt_vwr(m) * dust_density * &
gravit * slp_crc(m) / vsc_dyn_atm ! [m s-1] SeP97 p.466
end do
! For Reynolds number flows Re < 0.1 Stokes' velocity is valid for
! vlc_grv SeP97 p. 466 (8.42). For larger Re, inertial effects become
! important and empirical drag coefficients must be employed
! Implicit equation for Re, Cd, and Vt is SeP97 p. 467 (8.44)
! Using Stokes' velocity rather than iterative solution with empirical
! drag coefficient causes 60% errors for D = 200 um SeP97 p. 468
! Iterative solution for drag coefficient, Reynolds number, and terminal veloc
do m = 1, nbin
! Initialize accuracy and counter
eps_crr = eps_max + 1.0_r8 ![frc] Current relative accuracy
itr_idx = 0 ![idx] Counting index
! Initial guess for vlc_grv is exact for Re < 0.1
vlc_grv(m) = vlc_stk(m) ![m s-1]
eps_loop: do while(eps_crr > eps_max)
! Save terminal velocity for convergence test
vlc_grv_old = vlc_grv(m) ![m s-1]
ryn_nbr_grv(m) = vlc_grv(m) * dmt_vwr(m) / vsc_knm_atm !SeP97 p.460
! Update drag coefficient based on new Reynolds number
if (ryn_nbr_grv(m) < 0.1_r8) then
cff_drg_grv(m) = 24.0_r8 / ryn_nbr_grv(m) !Stokes' law Sep97 p.463 (8.32)
else if (ryn_nbr_grv(m) < 2.0_r8) then
cff_drg_grv(m) = (24.0_r8/ryn_nbr_grv(m)) * &
(1.0_r8 + 3.0_r8*ryn_nbr_grv(m)/16.0_r8 + &
9.0_r8*ryn_nbr_grv(m)*ryn_nbr_grv(m)* &
log(2.0_r8*ryn_nbr_grv(m))/160.0_r8) !Sep97 p.463 (8.32)
else if (ryn_nbr_grv(m) < 500.0_r8) then
cff_drg_grv(m) = (24.0_r8/ryn_nbr_grv(m)) * &
(1.0_r8 + 0.15_r8*ryn_nbr_grv(m)**0.687_r8) !Sep97 p.463 (8.32)
else if (ryn_nbr_grv(m) < 2.0e5_r8) then
cff_drg_grv(m) = 0.44_r8 !Sep97 p.463 (8.32)
else
write(iulog,'(a,es9.2)') "ryn_nbr_grv(m) = ",ryn_nbr_grv(m)
call endrun ('Dustini error: Reynolds number too large in stk_crc_get()')
endif
! Update terminal velocity based on new Reynolds number and drag coeff
! [m s-1] Terminal veloc SeP97 p.467 (8.44)
vlc_grv(m) = sqrt(4.0_r8 * gravit * dmt_vwr(m) * slp_crc(m) * dust_density / &
(3.0_r8*cff_drg_grv(m)*dns_mdp))
eps_crr = abs((vlc_grv(m)-vlc_grv_old)/vlc_grv(m)) !Relative convergence
if (itr_idx == 12) then
! Numerical pingpong may occur when Re = 0.1, 2.0, or 500.0
! due to discontinuities in derivative of drag coefficient
vlc_grv(m) = 0.5_r8 * (vlc_grv(m)+vlc_grv_old) ! [m s-1]
endif
if (itr_idx > 20) then
write(iulog,*) 'Dustini error: Terminal velocity not converging ',&
' in stk_crc_get(), breaking loop...'
! to next iteration
exit eps_loop
endif
itr_idx = itr_idx + 1
end do eps_loop !end while
end do !end loop over size
! Compute factors to convert Stokes' settling velocities to
! actual settling velocities
do m = 1, nbin
stk_crc(m) = vlc_grv(m) / vlc_stk(m)
end do
end subroutine dust_set_params
end module dust_common