| blob | 7fc24c8150652e0e990e8cf7529d53732d3076f2 |
1 MODULE qf03\r
2 !\r
3 ! Y. Shao, 29 Jan 2004 \r
4 ! \r
5 ! JY Kang, 01 Dec 2008\r
6 ! Modify the code for WRF_chem\r
7 !\r
8 ! M. Klose, 2010-2015 Modifications\r
9 \r
10 ! A. Ukhov, 11 March 2021, Now "bems" is dust emission (kg/m2) per dt.\r
11 ! Before was instantenious flux (kg m^-2 s^-1).\r
12 \r
13 !-----------------------------------------------------------------------------------\r
14 ! Calculate sediment flux for multi-particle size soils as a weighted average of Q(d) \r
15 ! dust emission F(d) for covered and moisture soil\r
16 !\r
17 ! Options for dust calculation:\r
18 ! 1 Shao (2001)\r
19 ! 2 Shao (2004)\r
20 ! 3 Shao (2011): simplification of 2; added on 26 Sep 2009\r
21 !\r
22 !--------------------------------------------------------------------------------------\r
23 !\r
24 ! input: \r
25 ! n: number of particle size ranges.\r
26 ! dm: median diameter of each particle size. [m]\r
27 ! m_fract: Weight fraction of each particle range. Sum m_fract = 1\r
28 ! ustar: Friciton velocity. [m/s]\r
29 ! cf: fraction area covered by roughness elements/vegetation cover\r
30 ! w: surface soil moisture content [m^3/m^3]\r
31 ! wr: air-dry soil moisture [m^3/m^3] from SOILPARM.TBL\r
32 ! c: Owen's coefficient\r
33 !\r
34 ! output: \r
35 ! ustart: Mean threshold velocity of each particle range.\r
36 ! q: Sand flux from each size range.\r
37 ! ffq: Weighted sand flux from each size range.\r
38 ! qtotal: Total sand flux (weighted average).\r
39 ! f: Weighted dust flux ejected by single sand size range for given d_d. \r
40 ! fff: Weighted dust flux ejected by all sand size range from each dust size range.\r
41 ! ftotal: Total dust flux (weighted average).\r
42 !--------------------------------------------------------------------------------------\r
43 !\r
44 \r
45 CONTAINS\r
46 \r
47 subroutine qf03_driver ( nmx, idst, g, rhop, rho, dt, &\r
48 ustar, w, cf, ust_min, imod, dz_lowest, &\r
49 soilc, tot_soilc, domsoilc, &\r
50 tc, bems, rough_cor_in, smois_cor_in, wr, &\r
51 d0, dd, psd_m, dpsd_m, ppsd_m, psd_f, dpsd_f, ppsd_f, imax, stype)\r
52 \r
53 INTEGER, INTENT(IN ) :: nmx, imod, idst\r
54 REAL , INTENT(IN ) :: g, dt, rho, dz_lowest, ustar\r
55 REAL(8), INTENT(IN ) :: w\r
56 REAL(8), INTENT(INOUT) :: cf\r
57 REAL(8), INTENT( OUT) :: ust_min, rough_cor_in, smois_cor_in\r
58 REAL , INTENT(INOUT), DIMENSION(nmx) :: tc\r
59 REAL , INTENT( OUT), DIMENSION(nmx) :: bems\r
60 REAL , INTENT(IN ), DIMENSION(12) :: wr\r
61 \r
62 !kang [2009/01/07] soil class type\r
63 REAL, INTENT(IN ) :: tot_soilc \r
64 REAL, DIMENSION(16), INTENT(IN) :: soilc \r
65 INTEGER, INTENT(IN ) :: domsoilc\r
66 \r
67 integer :: imax, stype\r
68 real(8), dimension(0:imax), intent(in) :: d0\r
69 real(8), dimension(imax), intent(in) :: dd \r
70 real(8), dimension(imax,stype), intent(in) :: psd_m, dpsd_m, ppsd_m\r
71 real(8), dimension(imax,stype), intent(in) :: psd_f, dpsd_f, ppsd_f\r
72 \r
73 !local variables\r
74 ! particle-size distributions \r
75 real(8), dimension(imax) :: psdm, dpsdm, ppsdm ! minimally dispersed\r
76 real(8), dimension(imax) :: psdf, dpsdf, ppsdf ! fully dispersed\r
77 real(8), dimension(imax) :: psds, dpsds, ppsds ! sediment \r
78 \r
79 real(8) :: smois_correc\r
80 integer :: i, j, n, kk, ij, ffile\r
81 real(8) :: total, qtotal, ftotal\r
82 real(8) :: ftotalb, ftotalp\r
83 !\r
84 real(8), dimension(imax) :: beta_d, beta_s ! beta1 and beta2 used by Shao et al. (1996) dust emission\r
85 real(8), dimension(imax) :: ustart, q, ffq, fff\r
86 real(8), dimension(imax,imax) :: f\r
87 !\r
88 real(8), parameter :: c_lambda=0.35\r
89 real(8) :: h, lambda\r
90 real(8) :: ghl, fc\r
91 real(8) :: phl\r
92 real(8) :: cys, u0, al0, rhos, smass, omega, rys\r
93 real(8) :: ddm, a1, a2, a3 \r
94 real(8) :: zeta, sigma_m ! u*sqrt(rhos/p), bombardment coefficient\r
95 !\r
96 real(8) :: ustart0_out, qwhite_out, f_mb_out, f_hlys_out, pmass_out, vhlys_out\r
97 ! \r
98 integer, parameter :: nbins = 5 \r
99 real(8) :: sigma\r
100 real, dimension(nbins) :: dbin, fbin, cell_fbin\r
101 integer, dimension(nbins) :: ibin\r
102 data dbin/2.,3.6,6.,12.,20./ !size cut diameter (um) consistent with GOCART model\r
103 real(8), intent(in) :: rhop\r
104 real :: cell_ftotal\r
105 integer :: isl, cc\r
106 character :: msg, soil\r
107 !*******************************************************************************************\r
108 !\r
109 ! initialization\r
110 cell_ftotal = 0.\r
111 do n = 1, nbins\r
112 cell_fbin(n) = 0.\r
113 enddo\r
114 ! \r
115 DO cc = 1, 12 ! soil category\r
116 if (soilc(cc).eq.0.) then\r
117 go to 103\r
118 endif\r
119 ! \r
120 psdm(:)=psd_m(:,cc)\r
121 psdf(:)=psd_f(:,cc)\r
122 dpsdm(:)=dpsd_m(:,cc)\r
123 dpsdf(:)=dpsd_f(:,cc)\r
124 ppsdm(:)=ppsd_m(:,cc)\r
125 ppsdf(:)=ppsd_f(:,cc)\r
126 \r
127 ! default settings: \r
128 rhos = 1000.d0 ! bulk density of soil [kg/m3] \r
129 phl = 30000. ! plastic pressure [N/m2]\r
130 cys = 0.00001 ! cys : parameter\r
131 \r
132 sigma = rhop/rho ! particle-air density ratio\r
133 \r
134 !\r
135 !-------------------\r
136 ! frontal area index\r
137 !-------------------\r
138 if (cf .gt. 0.95) then \r
139 write(6,*) 'cover fraction too large, reset'\r
140 cf = 0.95\r
141 endif! as larger values lead to large lambda, which is problematic in Raupach scheme\r
142 \r
143 lambda = - c_lambda*dlog( 1.d0 - cf )\r
144 call r_c(lambda, rough_cor_in)\r
145 \r
146 ! Matching WRF_soil class with Shao's class for moisture correction of Fecan (cc:WRF_soil class, isl:Shao's class)\r
147 ! cc\r
148 ! 1:sand, 2:loamy sand, 3:sandy loam, 4:silt loam, 5:silt, 6:loam, 7:sandy clay loam,\r
149 ! 8:silty clay loam, 9:clay loam, 10:sandy clay, 11:silty clay, 12:clay\r
150 ! isl\r
151 ! 1:sand, 2:loamy sand, 3:sandy loam, 4:loam, 5:silt loam, 6: silt, 7:sandy clay loam, \r
152 ! 8:clay loam, 9:silty clay loam, 10:sandy clay, 11:silty clay, 12:clay\r
153 if (cc.eq.1.or.cc.eq.2.or.cc.eq.3.or.cc.eq.7.or.cc.eq.10.or. & \r
154 & cc.eq.11.or.cc.eq.12) then\r
155 isl = cc\r
156 elseif (cc.eq.4) then\r
157 isl = 5\r
158 elseif (cc.eq.5) then\r
159 isl = 6\r
160 elseif (cc.eq.6) then\r
161 isl = 4\r
162 elseif (cc.eq.8) then\r
163 isl = 9\r
164 elseif (cc.eq.9) then\r
165 isl = 8\r
166 endif\r
167 \r
168 call h_c(w, wr(cc), isl, smois_correc)\r
169 !----------------------------------------------\r
170 ! for each particle size group, estimate ustart\r
171 !----------------------------------------------\r
172 !\r
173 ust_min = 999.0\r
174 do i = 1, imax\r
175 call ustart0(dd(i), sigma, g, rho, ustart0_out)\r
176 ustart(i) = ustart0_out\r
177 ustart(i) = rough_cor_in*smois_correc*ustart(i)\r
178 ust_min = dmin1(ust_min, ustart(i))\r
179 call qwhite(ustart(i), ustar, rho, g, qwhite_out)\r
180 q(i) = qwhite_out\r
181 q(i) = (1.d0-cf)*q(i)\r
182 enddo\r
183 if (cc.eq.domsoilc) then\r
184 smois_cor_in = smois_correc\r
185 endif\r
186 ! \r
187 ! \r
188 IF ( ustar .le. ust_min ) THEN ! no erosion goto 102\r
189 q = 0.d0\r
190 ffq = 0.d0\r
191 qtotal = 0.d0\r
192 fff = 0.d0\r
193 ftotal = 0.d0\r
194 fbin = 0.d0\r
195 goto 102\r
196 ELSE\r
197 ghl = dexp( -(ustar - ust_min)**3.d0 )\r
198 dpsds = ghl*dpsdm + (1.-ghl)*dpsdf\r
199 psds = ghl*psdm + (1.-ghl)*psdf\r
200 ppsds = ghl*ppsdm + (1.-ghl)*ppsdf\r
201 !\r
202 ffq = q*dpsds\r
203 qtotal = sum(ffq)\r
204 \r
205 !--------------\r
206 ! dust emission\r
207 !--------------\r
208 ! \r
209 ! size bin\r
210 do n=1,nbins\r
211 ibin(n)=0\r
212 do i=imax,1,-1\r
213 if(d0(i).ge.dbin(n)) ibin(n)=i\r
214 enddo\r
215 if(ibin(n).eq.0) stop 'wrong dust classes'\r
216 enddo\r
217 ! \r
218 ! \r
219 ! \r
220 !--------------------------------\r
221 ! Shao (2001) dust emission model\r
222 !--------------------------------\r
223 IF (imod .eq. 1) THEN \r
224 do i = idst+1, imax\r
225 ddm = dd(i)*1.d-6\r
226 call pmass(rhop, ddm, pmass_out) ! mass of saltating particles\r
227 smass = pmass_out\r
228 u0 = 10*ustar\r
229 al0 = 13.d0*3.14159d0/180.d0\r
230 call vhlys(phl, 2, smass, al0, u0, ddm, vhlys_out)\r
231 omega = vhlys_out ![m3]\r
232 ! \r
233 do j = 1, idst\r
234 rys = psdm(j)/psdf(j)\r
235 if (rys.gt.1.) then \r
236 rys = 1.\r
237 endif\r
238 a1 = cys*( (1.-ghl) + ghl*rys )\r
239 a2 = ffq(i)*g/ustar**2/smass\r
240 \r
241 if ( dpsdf(j) .lt. dpsdm(j) ) then \r
242 a3 = dpsdf(j)*rhos*omega\r
243 else \r
244 a3 = dpsdf(j)*rhos*omega + (dpsdf(j)-dpsdm(j))*smass\r
245 endif\r
246 f(i,j) = a1*a2*a3 ![kg/m2/s]\r
247 enddo\r
248 enddo\r
249 ! \r
250 ftotal = 0.0\r
251 do j = 1, idst\r
252 fff(j) = 0.\r
253 do i = idst+1, imax\r
254 fff(j) = fff(j) + f(i,j)\r
255 enddo\r
256 fff(j) = (1.d0-cf)*fff(j)\r
257 ftotal = ftotal + fff(j)\r
258 enddo\r
259 ! \r
260 do n=1,nbins\r
261 j0=1\r
262 if(n.gt.1) j0=ibin(n-1)+1\r
263 fbin(n)=0\r
264 do j=j0,ibin(n)\r
265 fbin(n)=fbin(n)+fff(j)\r
266 enddo\r
267 enddo\r
268 ! \r
269 !--------------------------------\r
270 ! Shao (2004) dust emission model\r
271 !--------------------------------\r
272 ELSEIF (imod .eq. 2) THEN\r
273 zeta = ustar*dsqrt( rhos/phl )\r
274 sigma_m = 12.d0*zeta*zeta*(1.d0+14.d0*zeta)\r
275 ! \r
276 do i = idst+1, imax\r
277 do j = 1, idst\r
278 rys = psdm(j)/psdf(j)\r
279 if (rys .gt.1.) then\r
280 rys = 1.\r
281 endif\r
282 a1 = cys*dpsdf(j)*( (1.-ghl) + ghl*rys )\r
283 a2 = (1.+sigma_m)\r
284 a3 = ffq(i)*g/ustar**2\r
285 f(i,j) = a1*a2*a3\r
286 enddo\r
287 enddo\r
288 ! \r
289 ftotal = 0.0\r
290 do j = 1, idst\r
291 fff(j) = 0.\r
292 do i = idst+1, imax\r
293 fff(j) = fff(j) + f(i,j)\r
294 enddo\r
295 fff(j) = (1.d0-cf)*fff(j) \r
296 ftotal = ftotal + fff(j)\r
297 enddo\r
298 ! \r
299 do n=1,nbins\r
300 j0=1\r
301 if(n.gt.1) j0=ibin(n-1)+1\r
302 fbin(n)=0\r
303 do j=j0,ibin(n)\r
304 fbin(n)=fbin(n)+fff(j) \r
305 enddo\r
306 enddo\r
307 ! \r
308 ! \r
309 !--------------------------------------------------------------------------\r
310 ! Shao (2011) minimal version, ghl = 1, Q independent of sand particle size\r
311 ! \r
312 ! See Eq. (34) in \r
313 ! Shao, Y., M. Ishizuka, M. Mikami, J. Leys (2011): Parameterization of size-\r
314 ! resolved dust emission and validation with measurements, JGR, 116, D08203, \r
315 ! doi: 10.1029/2010JD014527\r
316 !--------------------------------------------------------------------------\r
317 ELSEIF (imod .eq. 3) THEN \r
318 zeta = ustar*dsqrt( rhos/phl )\r
319 sigma_m = 12.d0*zeta*zeta*(1.d0+14.d0*zeta)\r
320 ! \r
321 ftotal = 0.0 \r
322 ! \r
323 do j = 1, idst\r
324 a1 = cys*dpsdm(j)\r
325 a2 = (1.+sigma_m)\r
326 a3 = qtotal*g/ustar**2\r
327 fff(j) = a1*a2*a3\r
328 fff(j) = (1.d0-cf)*fff(j)\r
329 ftotal = ftotal + fff(j)\r
330 enddo\r
331 ! \r
332 do n=1,nbins\r
333 j0=1\r
334 if(n.gt.1) j0=ibin(n-1)+1\r
335 fbin(n)=0\r
336 do j=j0,ibin(n)\r
337 fbin(n)=fbin(n)+fff(j) \r
338 enddo\r
339 enddo\r
340 ! \r
341 ENDIF ! dust scheme\r
342 ENDIF ! ust < ust_t\r
343 ! \r
344 ! \r
345 ! \r
346 102 continue\r
347 \r
348 do n = 1, nbins\r
349 cell_fbin(n) = cell_fbin(n) + (soilc(cc)/tot_soilc)*fbin(n)\r
350 enddo\r
351 \r
352 103 continue \r
353 \r
354 ENDDO ! cc, soil category\r
355 \r
356 do n = 1, nbins\r
357 ! fbin : [kg/m2/s], dz_lowest : [m], rho : [kg/m3], dt : [s] -> tc : [kg/kg-dryair]\r
358 tc(n) = tc(n) + cell_fbin(n)/dz_lowest/rho*dt ! (kg/kg)\r
359 bems(n) = cell_fbin(n)*dt ! diagnostic (kg/m2) per dt\r
360 enddo\r
361 \r
362 \r
363 END subroutine qf03_driver\r
364 \r
365 \r
366 !*****************************************************************************\r
367 subroutine ustart0(dum, sigma, g, rho, ustart0_out)\r
368 !\r
369 ! Y. Shao, 13 June 2000\r
370 !\r
371 ! Calculate ustar0(d) using Shao and Lu (2000) for uncovered\r
372 ! dry surface\r
373 !\r
374 ! dum: particle diameter [um]\r
375 ! ustar0: threshold friction velocity [m/s]\r
376 !\r
377 real, intent(in) :: g, rho\r
378 real(8), intent(in) :: dum, sigma\r
379 real(8), intent(out) :: ustart0_out\r
380 real(8) :: dm\r
381 real(8), parameter :: gamma = 1.65d-4 ! a constant\r
382 real(8), parameter :: f = 0.0123\r
383 \r
384 dm = dum*1d-6\r
385 \r
386 ustart0_out = f*(sigma*g*dm + gamma/(rho*dm) )\r
387 ustart0_out = dsqrt( ustart0_out )\r
388 ! end function\r
389 end subroutine ustart0\r
390 !*****************************************************************************\r
391 subroutine qwhite(ust, ustar, rho, g, qwhite_out)\r
392 !\r
393 ! Yaping Shao 17-07-99!\r
394 !\r
395 ! White (1979) Sand Flux Equation\r
396 ! Q = c*rho*u_*^3 over g (1 - u_*t over u_*)(1 + u_*t^2/u_*^2)\r
397 ! qwhite: Streamwise Sand Flux; [kg m-1 s-1]\r
398 ! c : 2.6\r
399 ! ust : threhold friction velocity [m/s]\r
400 ! ustar : friction velocity [m/s]\r
401 !\r
402 real(8) :: c\r
403 real, intent(in) :: ustar, rho, g\r
404 real(8), intent(in) :: ust\r
405 real(8), intent(out) :: qwhite_out\r
406 real(8) :: a, b\r
407 ! default setting: \r
408 c = 2.3d0\r
409 a = rho/g\r
410 \r
411 IF (ustar.lt.ust) THEN \r
412 qwhite_out = 0.d0\r
413 ELSE\r
414 b = ust/ustar \r
415 qwhite_out = c*a*ustar**3.*(1.-b)*(1.+b)**2.\r
416 ENDIF\r
417 \r
418 END subroutine qwhite\r
419 !*****************************************************************************\r
420 subroutine vhlys(p, k, xm, alpha, u, d, vhlys_out)\r
421 !\r
422 ! Volume removal according to Lu and Shao (1999), Equation (8)\r
423 ! alpha: impact angle [^o]\r
424 ! u : impact velocity [m/s]\r
425 ! p : plastic pressure [N/m^2]\r
426 ! xm : particle mass [kg]\r
427 ! d : particle diameter [m]\r
428 !\r
429 !\r
430 REAL(8),intent(in) :: alpha, xm, u, d, p\r
431 REAL(8) :: beta\r
432 REAL(8), PARAMETER :: pi=3.1415927d0\r
433 REAL(8) :: t1, t2, t3\r
434 INTEGER,intent(in) :: k\r
435 real(8),intent(out) :: vhlys_out\r
436 \r
437 beta = dsqrt( p*k*d/xm )\r
438 t1 = u*u/(beta*beta)*( dsin(2.d0*alpha) - 4.d0*dsin(alpha)*dsin(alpha) )\r
439 t2 = u*dsin(alpha)/beta\r
440 t3 = 7.5d0*pi*t2**3.d0/d\r
441 vhlys_out = d*( t1 + t3 )\r
442 \r
443 END subroutine vhlys \r
444 !*****************************************************************************\r
445 ! A routine for correction of ust for soil moisture content\r
446 !\r
447 ! w : volumetric soil moisture\r
448 ! isl: soil texture type, ranging from 1 to 12\r
449 !\r
450 ! Author: Yaping Shao, 5/05/2001\r
451 ! Reference: Fecan et al. (1999), Ann. Geophysicae,17,149-157\r
452 !\r
453 ! Data based on Shao and Jung, 2000, unpublished manuscript\r
454 ! Data invented for sand, loamy sand, sandy loam, loam, clay loam, and clay\r
455 ! Soil classes:\r
456 ! 1 sand, 2 loamy sand, 3 sandy loam, 4 loam, 5 silty loam, 6 silt, 7 sandy clay loam, 8 clay loam, \r
457 ! 9 silty clay loam, 10 sandy clay, 11 silty clay, 12 clay\r
458 !----------------------------------------------------------------------\r
459 subroutine h_c (w, wr, isl, h)\r
460 \r
461 real(8) :: a(12), b(12)!, thr(12) \r
462 real(8), intent(in) :: w\r
463 real(8), intent(out) :: h\r
464 real, intent(in) :: wr\r
465 integer, intent(in) :: isl\r
466 character*100 :: msg ! error message string\r
467 \r
468 ! data thr/0.001, 0.003, 0.037, 0.049, 0.061, 0.072, 0.084, 0.095, 0.110, 0.126, 0.141, 0.156/\r
469 data a /21.19, 33.03, 44.87, 17.79, 20.81, 23.83, 26.84, 29.86, 27.51, 25.17, 22.82, 20.47/\r
470 data b / 0.68, 0.71, 0.85, 0.61, 0.66, 0.71, 0.75, 0.80, 0.75, 0.70, 0.64, 0.59/ \r
471 \r
472 if ( w.lt.0. ) then\r
473 write(msg, *) 'soil moisture correction (h_c): w = ', w, ' < 0'\r
474 call wrf_error_fatal(msg)\r
475 ! stop\r
476 endif\r
477 \r
478 if ( w.le.wr ) then\r
479 h = 1.0\r
480 else\r
481 h = sqrt( 1 + a(isl)*( w-wr )**b(isl) ) \r
482 endif\r
483 \r
484 END subroutine h_c\r
485 \r
486 !*****************************************************************************\r
487 subroutine pmass(rhop, d, pmass_out)\r
488 !\r
489 ! Particle Mass \r
490 ! rhop: particle density [kg m^-3]\r
491 ! d : particle size [m]\r
492 ! \r
493 REAL(8), PARAMETER :: pi=3.1415927d0\r
494 REAL(8),intent(in) :: rhop, d\r
495 real(8),intent(out) :: pmass_out\r
496 \r
497 pmass_out = (pi*rhop*d**3.d0)/6.d0\r
498 \r
499 END subroutine pmass\r
500 \r
501 !*****************************************************************************\r
502 subroutine r_c (x,r)\r
503 !\r
504 ! Y. Shao 17-07-92\r
505 ! CORRECTION FUNCTION FOR UST(D) BASED ON Raupach et al. (1992)\r
506 ! x = frontal area index\r
507 !\r
508 ! R_C = (1 - sig m x)^{1/2} (1 + m beta x)^{1/2}\r
509 ! Note I deife R_C = u_{*tR}/u_{*tS}\r
510 ! While Raupach et al. defined \r
511 ! R_C = u_{*tS}/u_{*tR} and their R function is\r
512 ! R_C = (1 - sig m x)^{-1/2} (1 + m beta x)^{-1/2}\r
513 ! \r
514 ! sig : basal to frontal area; about 1\r
515 ! m : parameter less than 1; about 0.5\r
516 ! beta : a ratio of drag coef.; about 90. ; changed to 200 (recommendation by Y. Shao based on data by Gillies et al.)\r
517 !\r
518 real(8) :: xc\r
519 real(8), intent(in) :: x\r
520 real(8), intent(out) :: r\r
521 real(8), parameter :: sig=1., m=0.5, beta=200.\r
522 !\r
523 xc = 1./(sig*m)\r
524 IF (x.ge.xc) THEN\r
525 r = 999. ! Full covered surface\r
526 ELSE\r
527 r = dsqrt(1.-sig*m*x)*dsqrt(1.+m*beta*x)\r
528 ENDIF\r
529 !\r
530 END subroutine r_c\r
531 \r
532 \r
533 \r
534 END MODULE qf03\r