Line data Source code
1 : module dust_sediment_mod
2 :
3 : !---------------------------------------------------------------------------------
4 : ! Purpose:
5 : !
6 : ! Contains routines to compute tendencies from sedimentation of dust
7 : !
8 : ! Author: Phil Rasch
9 : !
10 : !---------------------------------------------------------------------------------
11 :
12 : use shr_kind_mod, only: r8=>shr_kind_r8
13 : use ppgrid, only: pcols, pver, pverp
14 : use physconst, only: gravit, rair
15 : use cam_logfile, only: iulog
16 : use cam_abortutils, only: endrun
17 :
18 : private
19 : public :: dust_sediment_vel, dust_sediment_tend
20 :
21 :
22 : real (r8), parameter :: vland = 2.8_r8 ! dust fall velocity over land (cm/s)
23 : real (r8), parameter :: vocean = 1.5_r8 ! dust fall velocity over ocean (cm/s)
24 : real (r8), parameter :: mxsedfac = 0.99_r8 ! maximum sedimentation flux factor
25 :
26 : contains
27 :
28 : !===============================================================================
29 0 : subroutine dust_sediment_vel (ncol, &
30 : icefrac , landfrac, ocnfrac , pmid , pdel , t , &
31 : dustmr , pvdust )
32 :
33 : !----------------------------------------------------------------------
34 :
35 : ! Compute gravitational sedimentation velocities for dust
36 :
37 : implicit none
38 :
39 : ! Arguments
40 : integer, intent(in) :: ncol ! number of colums to process
41 :
42 : real(r8), intent(in) :: icefrac (pcols) ! sea ice fraction (fraction)
43 : real(r8), intent(in) :: landfrac(pcols) ! land fraction (fraction)
44 : real(r8), intent(in) :: ocnfrac (pcols) ! ocean fraction (fraction)
45 : real(r8), intent(in) :: pmid (pcols,pver) ! pressure of midpoint levels (Pa)
46 : real(r8), intent(in) :: pdel (pcols,pver) ! pressure diff across layer (Pa)
47 : real(r8), intent(in) :: t (pcols,pver) ! temperature (K)
48 : real(r8), intent(in) :: dustmr(pcols,pver) ! dust (kg/kg)
49 :
50 : real(r8), intent(out) :: pvdust (pcols,pverp) ! vertical velocity of dust (Pa/s)
51 : ! -> note that pvel is at the interfaces (loss from cell is based on pvel(k+1))
52 :
53 : ! Local variables
54 : real (r8) :: rho(pcols,pver) ! air density in kg/m3
55 : real (r8) :: vfall(pcols) ! settling velocity of dust particles (m/s)
56 :
57 : integer i,k
58 :
59 : real (r8) :: lbound, ac, bc, cc
60 :
61 : !-----------------------------------------------------------------------
62 : !--------------------- dust fall velocity ----------------------------
63 : !-----------------------------------------------------------------------
64 :
65 0 : do k = 1,pver
66 0 : do i = 1,ncol
67 :
68 : ! merge the dust fall velocities for land and ocean (cm/s)
69 : ! SHOULD ALSO ACCOUNT FOR ICEFRAC
70 0 : vfall(i) = vland*landfrac(i) + vocean*(1._r8-landfrac(i))
71 : !! vfall(i) = vland*landfrac(i) + vocean*ocnfrac(i) + vseaice*icefrac(i)
72 :
73 : ! fall velocity (assume positive downward)
74 0 : pvdust(i,k+1) = vfall(i)
75 : end do
76 : end do
77 :
78 0 : return
79 : end subroutine dust_sediment_vel
80 :
81 :
82 : !===============================================================================
83 56588688 : subroutine dust_sediment_tend ( &
84 : ncol, dtime, pint, pmid, pdel, t, &
85 : dustmr ,pvdust, dusttend, sfdust )
86 :
87 : !----------------------------------------------------------------------
88 : ! Apply Particle Gravitational Sedimentation
89 : !----------------------------------------------------------------------
90 :
91 : implicit none
92 :
93 : ! Arguments
94 : integer, intent(in) :: ncol ! number of colums to process
95 :
96 : real(r8), intent(in) :: dtime ! time step
97 : real(r8), intent(in) :: pint (pcols,pverp) ! interfaces pressure (Pa)
98 : real(r8), intent(in) :: pmid (pcols,pver) ! midpoint pressures (Pa)
99 : real(r8), intent(in) :: pdel (pcols,pver) ! pressure diff across layer (Pa)
100 : real(r8), intent(in) :: t (pcols,pver) ! temperature (K)
101 : real(r8), intent(in) :: dustmr(pcols,pver) ! dust (kg/kg)
102 : real(r8), intent(in) :: pvdust (pcols,pverp) ! vertical velocity of dust drops (Pa/s)
103 : ! -> note that pvel is at the interfaces (loss from cell is based on pvel(k+1))
104 :
105 : real(r8), intent(out) :: dusttend(pcols,pver) ! dust tend
106 : real(r8), intent(out) :: sfdust (pcols) ! surface flux of dust (rain, kg/m/s)
107 :
108 : ! Local variables
109 : real(r8) :: fxdust(pcols,pverp) ! fluxes at the interfaces, dust (positive = down)
110 :
111 : integer :: i,k
112 : !----------------------------------------------------------------------
113 :
114 : ! initialize variables
115 88933729248 : fxdust (:ncol,:) = 0._r8 ! flux at interfaces (dust)
116 87932241072 : dusttend(:ncol,:) = 0._r8 ! tend (dust)
117 944899488 : sfdust(:ncol) = 0._r8 ! sedimentation flux out bot of column (dust)
118 :
119 : ! fluxes at interior points
120 56588688 : call getflx(ncol, pint, dustmr, pvdust, dtime, fxdust)
121 :
122 : ! calculate fluxes at boundaries
123 944899488 : do i = 1,ncol
124 888310800 : fxdust(i,1) = 0
125 : ! surface flux by upstream scheme
126 944899488 : fxdust(i,pverp) = dustmr(i,pver) * pvdust(i,pverp) * dtime
127 : end do
128 :
129 : ! filter out any negative fluxes from the getflx routine
130 5262747984 : do k = 2,pver
131 86987341584 : fxdust(:ncol,k) = max(0._r8, fxdust(:ncol,k))
132 : end do
133 :
134 : ! Limit the flux out of the bottom of each cell to the water content in each phase.
135 : ! Apply mxsedfac to prevent generating very small negative cloud water/ice
136 : ! NOTE, REMOVED CLOUD FACTOR FROM AVAILABLE WATER. ALL CLOUD WATER IS IN CLOUDS.
137 : ! ***Should we include the flux in the top, to allow for thin surface layers?
138 : ! ***Requires simple treatment of cloud overlap, already included below.
139 5319336672 : do k = 1,pver
140 87932241072 : do i = 1,ncol
141 87875652384 : fxdust(i,k+1) = min( fxdust(i,k+1), mxsedfac * dustmr(i,k) * pdel(i,k) )
142 : !!$ fxdust(i,k+1) = min( fxdust(i,k+1), dustmr(i,k) * pdel(i,k) + fxdust(i,k))
143 : end do
144 : end do
145 :
146 : ! Now calculate the tendencies
147 5319336672 : do k = 1,pver
148 87932241072 : do i = 1,ncol
149 : ! net flux into cloud changes cloud dust/ice (all flux is out of cloud)
150 87875652384 : dusttend(i,k) = (fxdust(i,k) - fxdust(i,k+1)) / (dtime * pdel(i,k))
151 : end do
152 : end do
153 :
154 : ! convert flux out the bottom to mass units Pa -> kg/m2/s
155 944899488 : sfdust(:ncol) = fxdust(:ncol,pverp) / (dtime*gravit)
156 :
157 56588688 : return
158 : end subroutine dust_sediment_tend
159 :
160 : !===============================================================================
161 56588688 : subroutine getflx(ncol, xw, phi, vel, deltat, flux)
162 :
163 : !.....xw1.......xw2.......xw3.......xw4.......xw5.......xw6
164 : !....psiw1.....psiw2.....psiw3.....psiw4.....psiw5.....psiw6
165 : !....velw1.....velw2.....velw3.....velw4.....velw5.....velw6
166 : !.........phi1......phi2.......phi3.....phi4.......phi5.......
167 :
168 :
169 : implicit none
170 :
171 : integer ncol ! number of colums to process
172 :
173 : integer i
174 : integer k
175 :
176 : real (r8) vel(pcols,pverp)
177 : real (r8) flux(pcols,pverp)
178 : real (r8) xw(pcols,pverp)
179 : real (r8) psi(pcols,pverp)
180 : real (r8) phi(pcols,pverp-1)
181 : real (r8) fdot(pcols,pverp)
182 : real (r8) xx(pcols)
183 : real (r8) fxdot(pcols)
184 : real (r8) fxdd(pcols)
185 :
186 : real (r8) psistar(pcols)
187 : real (r8) deltat
188 :
189 : real (r8) xxk(pcols,pver)
190 :
191 944899488 : do i = 1,ncol
192 : ! integral of phi
193 888310800 : psi(i,1) = 0._r8
194 : ! fluxes at boundaries
195 888310800 : flux(i,1) = 0
196 944899488 : flux(i,pverp) = 0._r8
197 : end do
198 :
199 : ! integral function
200 5319336672 : do k = 2,pverp
201 87932241072 : do i = 1,ncol
202 87875652384 : psi(i,k) = phi(i,k-1)*(xw(i,k)-xw(i,k-1)) + psi(i,k-1)
203 : end do
204 : end do
205 :
206 :
207 : ! calculate the derivatives for the interpolating polynomial
208 56588688 : call cfdotmc_pro (ncol, xw, psi, fdot)
209 :
210 : ! NEW WAY
211 : ! calculate fluxes at interior pts
212 5262747984 : do k = 2,pver
213 86987341584 : do i = 1,ncol
214 86930752896 : xxk(i,k) = xw(i,k)-vel(i,k)*deltat
215 : end do
216 : end do
217 5262747984 : do k = 2,pver
218 5206159296 : call cfint2(ncol, xw, psi, fdot, xxk(1,k), fxdot, fxdd, psistar)
219 86987341584 : do i = 1,ncol
220 86930752896 : flux(i,k) = (psi(i,k)-psistar(i))
221 : end do
222 : end do
223 :
224 :
225 56588688 : return
226 : end subroutine getflx
227 :
228 :
229 :
230 : !##############################################################################
231 :
232 5206159296 : subroutine cfint2 (ncol, x, f, fdot, xin, fxdot, fxdd, psistar)
233 :
234 :
235 : implicit none
236 :
237 : ! input
238 : integer ncol ! number of colums to process
239 :
240 : real (r8) x(pcols, pverp)
241 : real (r8) f(pcols, pverp)
242 : real (r8) fdot(pcols, pverp)
243 : real (r8) xin(pcols)
244 :
245 : ! output
246 : real (r8) fxdot(pcols)
247 : real (r8) fxdd(pcols)
248 : real (r8) psistar(pcols)
249 :
250 : integer i
251 : integer k
252 : integer intz(pcols)
253 : real (r8) dx
254 : real (r8) s
255 : real (r8) c2
256 : real (r8) c3
257 : real (r8) xx
258 : real (r8) xinf
259 : real (r8) psi1, psi2, psi3, psim
260 : real (r8) cfint
261 : real (r8) cfnew
262 : real (r8) xins(pcols)
263 :
264 : ! the minmod function
265 : real (r8) a, b, c
266 : real (r8) minmod
267 : real (r8) medan
268 : minmod(a,b) = 0.5_r8*(sign(1._r8,a) + sign(1._r8,b))*min(abs(a),abs(b))
269 : medan(a,b,c) = a + minmod(b-a,c-a)
270 :
271 86930752896 : do i = 1,ncol
272 81724593600 : xins(i) = medan(x(i,1), xin(i), x(i,pverp))
273 86930752896 : intz(i) = 0
274 : end do
275 :
276 : ! first find the interval
277 >48937*10^7 : do k = 1,pverp-1
278 >80897*10^8 : do i = 1,ncol
279 >80845*10^8 : if ((xins(i)-x(i,k))*(x(i,k+1)-xins(i)).ge.0._r8) then
280 81724593600 : intz(i) = k
281 : endif
282 : end do
283 : end do
284 :
285 86930752896 : do i = 1,ncol
286 86930752896 : if (intz(i).eq.0) then
287 0 : write(iulog,*) ' interval was not found for col i ', i
288 0 : call endrun('DUST_SEDIMENT_MOD:cfint2 -- interval was not found ')
289 : endif
290 : end do
291 :
292 : ! now interpolate
293 86930752896 : do i = 1,ncol
294 81724593600 : k = intz(i)
295 81724593600 : dx = (x(i,k+1)-x(i,k))
296 81724593600 : s = (f(i,k+1)-f(i,k))/dx
297 81724593600 : c2 = (3*s-2*fdot(i,k)-fdot(i,k+1))/dx
298 81724593600 : c3 = (fdot(i,k)+fdot(i,k+1)-2*s)/dx**2
299 81724593600 : xx = (xins(i)-x(i,k))
300 81724593600 : fxdot(i) = (3*c3*xx + 2*c2)*xx + fdot(i,k)
301 81724593600 : fxdd(i) = 6*c3*xx + 2*c2
302 81724593600 : cfint = ((c3*xx + c2)*xx + fdot(i,k))*xx + f(i,k)
303 :
304 : ! limit the interpolant
305 81724593600 : psi1 = f(i,k)+(f(i,k+1)-f(i,k))*xx/dx
306 81724593600 : if (k.eq.1) then
307 1348699599 : psi2 = f(i,1)
308 : else
309 80375894001 : psi2 = f(i,k) + (f(i,k)-f(i,k-1))*xx/(x(i,k)-x(i,k-1))
310 : endif
311 81724593600 : if (k+1.eq.pverp) then
312 0 : psi3 = f(i,pverp)
313 : else
314 81724593600 : psi3 = f(i,k+1) - (f(i,k+2)-f(i,k+1))*(dx-xx)/(x(i,k+2)-x(i,k+1))
315 : endif
316 81724593600 : psim = medan(psi1, psi2, psi3)
317 81724593600 : cfnew = medan(cfint, psi1, psim)
318 : if (abs(cfnew-cfint)/(abs(cfnew)+abs(cfint)+1.e-36_r8) .gt..03_r8) then
319 : ! CHANGE THIS BACK LATER!!!
320 : ! $ .gt..1) then
321 :
322 :
323 : ! UNCOMMENT THIS LATER!!!
324 : ! write(iulog,*) ' cfint2 limiting important ', cfint, cfnew
325 :
326 :
327 : endif
328 86930752896 : psistar(i) = cfnew
329 : end do
330 :
331 5206159296 : return
332 : end subroutine cfint2
333 :
334 :
335 :
336 : !##############################################################################
337 :
338 56588688 : subroutine cfdotmc_pro (ncol, x, f, fdot)
339 :
340 : ! prototype version; eventually replace with final SPITFIRE scheme
341 :
342 : ! calculate the derivative for the interpolating polynomial
343 : ! multi column version
344 :
345 :
346 : implicit none
347 :
348 : ! input
349 : integer ncol ! number of colums to process
350 :
351 : real (r8) x(pcols, pverp)
352 : real (r8) f(pcols, pverp)
353 : ! output
354 : real (r8) fdot(pcols, pverp) ! derivative at nodes
355 :
356 : ! assumed variable distribution
357 : ! x1.......x2.......x3.......x4.......x5.......x6 1,pverp points
358 : ! f1.......f2.......f3.......f4.......f5.......f6 1,pverp points
359 : ! ...sh1.......sh2......sh3......sh4......sh5.... 1,pver points
360 : ! .........d2.......d3.......d4.......d5......... 2,pver points
361 : ! .........s2.......s3.......s4.......s5......... 2,pver points
362 : ! .............dh2......dh3......dh4............. 2,pver-1 points
363 : ! .............eh2......eh3......eh4............. 2,pver-1 points
364 : ! ..................e3.......e4.................. 3,pver-1 points
365 : ! .................ppl3......ppl4................ 3,pver-1 points
366 : ! .................ppr3......ppr4................ 3,pver-1 points
367 : ! .................t3........t4.................. 3,pver-1 points
368 : ! ................fdot3.....fdot4................ 3,pver-1 points
369 :
370 :
371 : ! work variables
372 :
373 :
374 : integer i
375 : integer k
376 :
377 : real (r8) a ! work var
378 : real (r8) b ! work var
379 : real (r8) c ! work var
380 : real (r8) s(pcols,pverp) ! first divided differences at nodes
381 : real (r8) sh(pcols,pverp) ! first divided differences between nodes
382 : real (r8) d(pcols,pverp) ! second divided differences at nodes
383 : real (r8) dh(pcols,pverp) ! second divided differences between nodes
384 : real (r8) e(pcols,pverp) ! third divided differences at nodes
385 : real (r8) eh(pcols,pverp) ! third divided differences between nodes
386 : real (r8) pp ! p prime
387 : real (r8) ppl(pcols,pverp) ! p prime on left
388 : real (r8) ppr(pcols,pverp) ! p prime on right
389 : real (r8) qpl
390 : real (r8) qpr
391 : real (r8) ttt
392 : real (r8) t
393 : real (r8) tmin
394 : real (r8) tmax
395 : real (r8) delxh(pcols,pverp)
396 :
397 :
398 : ! the minmod function
399 : real (r8) minmod
400 : real (r8) medan
401 : minmod(a,b) = 0.5_r8*(sign(1._r8,a) + sign(1._r8,b))*min(abs(a),abs(b))
402 : medan(a,b,c) = a + minmod(b-a,c-a)
403 :
404 5319336672 : do k = 1,pver
405 :
406 :
407 : ! first divided differences between nodes
408 87875652384 : do i = 1, ncol
409 82612904400 : delxh(i,k) = (x(i,k+1)-x(i,k))
410 87875652384 : sh(i,k) = (f(i,k+1)-f(i,k))/delxh(i,k)
411 : end do
412 :
413 : ! first and second divided differences at nodes
414 5319336672 : if (k.ge.2) then
415 86930752896 : do i = 1,ncol
416 81724593600 : d(i,k) = (sh(i,k)-sh(i,k-1))/(x(i,k+1)-x(i,k-1))
417 86930752896 : s(i,k) = minmod(sh(i,k),sh(i,k-1))
418 : end do
419 : endif
420 : end do
421 :
422 : ! second and third divided diffs between nodes
423 5206159296 : do k = 2,pver-1
424 86042442096 : do i = 1, ncol
425 80836282800 : eh(i,k) = (d(i,k+1)-d(i,k))/(x(i,k+2)-x(i,k-1))
426 85985853408 : dh(i,k) = minmod(d(i,k),d(i,k+1))
427 : end do
428 : end do
429 :
430 : ! treat the boundaries
431 944899488 : do i = 1,ncol
432 888310800 : e(i,2) = eh(i,2)
433 888310800 : e(i,pver) = eh(i,pver-1)
434 : ! outside level
435 : fdot(i,1) = sh(i,1) - d(i,2)*delxh(i,1) &
436 888310800 : - eh(i,2)*delxh(i,1)*(x(i,1)-x(i,3))
437 888310800 : fdot(i,1) = minmod(fdot(i,1),3*sh(i,1))
438 : fdot(i,pverp) = sh(i,pver) + d(i,pver)*delxh(i,pver) &
439 888310800 : + eh(i,pver-1)*delxh(i,pver)*(x(i,pverp)-x(i,pver-1))
440 888310800 : fdot(i,pverp) = minmod(fdot(i,pverp),3*sh(i,pver))
441 : ! one in from boundary
442 888310800 : fdot(i,2) = sh(i,1) + d(i,2)*delxh(i,1) - eh(i,2)*delxh(i,1)*delxh(i,2)
443 888310800 : fdot(i,2) = minmod(fdot(i,2),3*s(i,2))
444 : fdot(i,pver) = sh(i,pver) - d(i,pver)*delxh(i,pver) &
445 888310800 : - eh(i,pver-1)*delxh(i,pver)*delxh(i,pver-1)
446 944899488 : fdot(i,pver) = minmod(fdot(i,pver),3*s(i,pver))
447 : end do
448 :
449 :
450 5149570608 : do k = 3,pver-1
451 85097542608 : do i = 1,ncol
452 85040953920 : e(i,k) = minmod(eh(i,k),eh(i,k-1))
453 : end do
454 : end do
455 :
456 :
457 :
458 5149570608 : do k = 3,pver-1
459 :
460 85097542608 : do i = 1,ncol
461 :
462 : ! p prime at k-0.5
463 79947972000 : ppl(i,k)=sh(i,k-1) + dh(i,k-1)*delxh(i,k-1)
464 : ! p prime at k+0.5
465 79947972000 : ppr(i,k)=sh(i,k) - dh(i,k) *delxh(i,k)
466 :
467 79947972000 : t = minmod(ppl(i,k),ppr(i,k))
468 :
469 : ! derivate from parabola thru f(i,k-1), f(i,k), and f(i,k+1)
470 79947972000 : pp = sh(i,k-1) + d(i,k)*delxh(i,k-1)
471 :
472 : ! quartic estimate of fdot
473 : fdot(i,k) = pp &
474 : - delxh(i,k-1)*delxh(i,k) &
475 79947972000 : *( eh(i,k-1)*(x(i,k+2)-x(i,k )) &
476 79947972000 : + eh(i,k )*(x(i,k )-x(i,k-2)) &
477 >15989*10^7 : )/(x(i,k+2)-x(i,k-2))
478 :
479 : ! now limit it
480 : qpl = sh(i,k-1) &
481 : + delxh(i,k-1)*minmod(d(i,k-1)+e(i,k-1)*(x(i,k)-x(i,k-2)), &
482 79947972000 : d(i,k) -e(i,k)*delxh(i,k))
483 : qpr = sh(i,k) &
484 : + delxh(i,k )*minmod(d(i,k) +e(i,k)*delxh(i,k-1), &
485 79947972000 : d(i,k+1)+e(i,k+1)*(x(i,k)-x(i,k+2)))
486 :
487 79947972000 : fdot(i,k) = medan(fdot(i,k), qpl, qpr)
488 :
489 79947972000 : ttt = minmod(qpl, qpr)
490 79947972000 : tmin = min(0._r8,3*s(i,k),1.5_r8*t,ttt)
491 79947972000 : tmax = max(0._r8,3*s(i,k),1.5_r8*t,ttt)
492 :
493 85040953920 : fdot(i,k) = fdot(i,k) + minmod(tmin-fdot(i,k), tmax-fdot(i,k))
494 :
495 : end do
496 :
497 : end do
498 :
499 56588688 : return
500 : end subroutine cfdotmc_pro
501 : end module dust_sediment_mod
|
