Line data Source code
1 : #undef OLDLIQSED
2 : module pkg_cld_sediment
3 :
4 : !---------------------------------------------------------------------------------
5 : ! Purpose:
6 : !
7 : ! Contains routines to compute tendencies from sedimentation of cloud liquid and
8 : ! ice particles
9 : !
10 : ! Author: Byron Boville Sept 19, 2002 from code by P. J. Rasch
11 : !
12 : !---------------------------------------------------------------------------------
13 :
14 : use shr_kind_mod, only: r8=>shr_kind_r8
15 : use spmd_utils, only: masterproc
16 : use ppgrid, only: pcols, pver, pverp
17 : use physconst, only: gravit, latvap, latice, rair, rhoh2o
18 : use cldwat, only: icritc
19 : use pkg_cldoptics, only: reitab, reltab
20 : use phys_control, only: use_simple_phys
21 : use cam_abortutils, only: endrun
22 : use cam_logfile, only: iulog
23 :
24 : implicit none
25 : private
26 : save
27 :
28 : public :: cld_sediment_readnl, cld_sediment_vel, cld_sediment_tend
29 :
30 :
31 : real (r8), parameter :: vland = 1.5_r8 ! liquid fall velocity over land (cm/s)
32 : real (r8), parameter :: vocean = 2.8_r8 ! liquid fall velocity over ocean (cm/s)
33 : real (r8), parameter :: mxsedfac = 0.99_r8 ! maximum sedimentation flux factor
34 :
35 : logical, parameter :: stokes = .true. ! use Stokes velocity instead of McFarquhar and Heymsfield
36 :
37 : ! parameter for modified McFarquhar and Heymsfield
38 : real (r8), parameter :: vice_small = 1._r8 ! ice fall velocity for small concentration (cm/s)
39 :
40 : ! parameters for Stokes velocity
41 : real (r8), parameter :: eta = 1.7e-5_r8 ! viscosity of air (kg m / s)
42 : real (r8), parameter :: r40 = 40._r8 ! 40 micron radius
43 : real (r8), parameter :: r400= 400._r8 ! 400 micron radius
44 : real (r8), parameter :: v400= 1.00_r8 ! fall velocity of 400 micron sphere (m/s)
45 : real (r8) :: v40 ! = (2._r8/9._r8) * rhoh2o * gravit/eta * r40**2 * 1.e-12_r8
46 : ! Stokes fall velocity of 40 micron sphere (m/s)
47 : real (r8) :: vslope ! = (v400 - v40)/(r400 -r40) ! linear slope for large particles m/s/micron
48 :
49 : ! namelist variables
50 : real(r8) :: cldsed_ice_stokes_fac = huge(1._r8) ! factor applied to the ice fall velocity computed from
51 : ! stokes terminal velocity
52 :
53 : !===============================================================================
54 : contains
55 : !===============================================================================
56 :
57 1536 : subroutine cld_sediment_readnl(nlfile)
58 :
59 : use namelist_utils, only: find_group_name
60 : use units, only: getunit, freeunit
61 : use mpishorthand
62 :
63 : character(len=*), intent(in) :: nlfile ! filepath for file containing namelist input
64 :
65 : ! Local variables
66 : integer :: unitn, ierr
67 : character(len=*), parameter :: subname = 'cld_sediment_readnl'
68 :
69 : namelist /cldsed_nl/ cldsed_ice_stokes_fac
70 : !-----------------------------------------------------------------------------
71 :
72 1536 : if (masterproc) then
73 2 : unitn = getunit()
74 2 : open( unitn, file=trim(nlfile), status='old' )
75 2 : call find_group_name(unitn, 'cldsed_nl', status=ierr)
76 2 : if (ierr == 0) then
77 2 : read(unitn, cldsed_nl, iostat=ierr)
78 2 : if (ierr /= 0) then
79 0 : call endrun(subname // ':: ERROR reading namelist')
80 : end if
81 : end if
82 2 : close(unitn)
83 2 : call freeunit(unitn)
84 :
85 : end if
86 :
87 : #ifdef SPMD
88 : ! Broadcast namelist variables
89 1536 : call mpibcast(cldsed_ice_stokes_fac, 1, mpir8, 0, mpicom)
90 : #endif
91 :
92 1536 : if (masterproc .and. .not. use_simple_phys) then
93 2 : write(iulog,*) subname//': cldsed_ice_stokes_fac = ', cldsed_ice_stokes_fac
94 : end if
95 :
96 1536 : end subroutine cld_sediment_readnl
97 :
98 : !===============================================================================
99 :
100 1495368 : subroutine cld_sediment_vel (ncol, &
101 : icefrac , landfrac, ocnfrac , pmid , pdel , t , &
102 : cloud , cldliq , cldice , pvliq , pvice , landm, snowh)
103 :
104 : !----------------------------------------------------------------------
105 :
106 : ! Compute gravitational sedimentation velocities for cloud liquid water
107 : ! and ice, based on Lawrence and Crutzen (1998).
108 :
109 : ! LIQUID
110 :
111 : ! The fall velocities assume that droplets have a gamma distribution
112 : ! with effective radii for land and ocean as assessed by Han et al.;
113 : ! see Lawrence and Crutzen (1998) for a derivation.
114 :
115 : ! ICE
116 :
117 : ! The fall velocities are based on data from McFarquhar and Heymsfield
118 : ! or on Stokes terminal velocity for spheres and the effective radius.
119 :
120 : ! NEED TO BE CAREFUL - VELOCITIES SHOULD BE AT THE *LOWER* INTERFACE
121 : ! (THAT IS, FOR K+1), FLUXES ARE ALSO AT THE LOWER INTERFACE (K+1),
122 : ! BUT MIXING RATIOS ARE AT THE MIDPOINTS (K)...
123 :
124 : ! NOTE THAT PVEL IS ON PVERP (INTERFACES), WHEREAS VFALL IS FOR THE CELL
125 : ! AVERAGES (I.E., MIDPOINTS); ASSUME THE FALL VELOCITY APPLICABLE TO THE
126 : ! LOWER INTERFACE (K+1) IS THE SAME AS THAT APPLICABLE FOR THE CELL (V(K))
127 :
128 : !-----------------------------------------------------------------------
129 : ! MATCH-MPIC version 2.0, Author: mgl, March 1998
130 : ! adapted by P. J. Rasch
131 : ! B. A. Boville, September 19, 2002
132 : ! P. J. Rasch May 22, 2003 (added stokes flow calc for liquid
133 : ! drops based on effect radii)
134 : !-----------------------------------------------------------------------
135 :
136 :
137 : ! Arguments
138 : integer, intent(in) :: ncol ! number of colums to process
139 :
140 : real(r8), intent(in) :: icefrac (pcols) ! sea ice fraction (fraction)
141 : real(r8), intent(in) :: landfrac(pcols) ! land fraction (fraction)
142 : real(r8), intent(in) :: ocnfrac (pcols) ! ocean fraction (fraction)
143 : real(r8), intent(in) :: pmid (pcols,pver) ! pressure of midpoint levels (Pa)
144 : real(r8), intent(in) :: pdel (pcols,pver) ! pressure diff across layer (Pa)
145 : real(r8), intent(in) :: cloud (pcols,pver) ! cloud fraction (fraction)
146 : real(r8), intent(in) :: t (pcols,pver) ! temperature (K)
147 : real(r8), intent(in) :: cldliq(pcols,pver) ! cloud water, liquid (kg/kg)
148 : real(r8), intent(in) :: cldice(pcols,pver) ! cloud water, ice (kg/kg)
149 : real(r8), intent(in) :: snowh(pcols) ! Snow depth over land, water equivalent (m)
150 :
151 : real(r8), intent(out) :: pvliq (pcols,pverp) ! vertical velocity of cloud liquid drops (Pa/s)
152 : real(r8), intent(out) :: pvice (pcols,pverp) ! vertical velocity of cloud ice particles (Pa/s)
153 : real(r8), intent(in) :: landm(pcols) ! land fraction ramped over water
154 : ! -> note that pvel is at the interfaces (loss from cell is based on pvel(k+1))
155 :
156 : ! Local variables
157 : real (r8) :: rho(pcols,pver) ! air density in kg/m3
158 : real (r8) :: vfall ! settling velocity of cloud particles (m/s)
159 : real (r8) :: icice ! in cloud ice water content (kg/kg)
160 : real (r8) :: iciwc ! in cloud ice water content in g/m3
161 : real (r8) :: icefac
162 : real (r8) :: logiwc
163 :
164 : real (r8) :: rei(pcols,pver) ! effective radius of ice particles (microns)
165 : real (r8) :: rel(pcols,pver) ! effective radius of liq particles (microns)
166 : real(r8) pvliq2 (pcols,pverp) ! vertical velocity of cloud liquid drops (Pa/s)
167 :
168 : integer i,k
169 :
170 : real (r8) :: lbound, ac, bc, cc
171 :
172 : !-----------------------------------------------------------------------
173 : !------- initialize linear ramp variables for fall velocity ------------
174 : !-----------------------------------------------------------------------
175 :
176 1495368 : v40 = (2._r8/9._r8) * rhoh2o * gravit/eta * r40**2 * 1.e-12_r8
177 1495368 : vslope = (v400 - v40)/(r400 -r40)
178 :
179 : !-----------------------------------------------------------------------
180 : !--------------------- liquid fall velocity ----------------------------
181 : !-----------------------------------------------------------------------
182 :
183 : ! compute air density
184 650693736 : rho(:ncol,:) = pmid(:ncol,:) / (rair * t(:ncol,:))
185 :
186 675662904 : pvliq(:ncol,:) = 0._r8
187 :
188 : ! get effective radius of liquid drop
189 1495368 : call reltab(ncol, t, landfrac, landm, icefrac, rel, snowh)
190 :
191 40374936 : do k = 1,pver
192 650693736 : do i = 1,ncol
193 649198368 : if (cloud(i,k) > 0._r8 .and. cldliq(i,k) > 0._r8) then
194 :
195 : #ifdef OLDLIQSED
196 : ! oldway
197 : ! merge the liquid fall velocities for land and ocean (cm/s)
198 : ! SHOULD ALSO ACCOUNT FOR ICEFRAC
199 : vfall = vland*landfrac(i) + vocean*(1._r8-landfrac(i))
200 : !!$ vfall = vland*landfrac(i) + vocean*ocnfrac(i) + vseaice*icefrac(i)
201 :
202 : ! convert the fall speed to pressure units, but do not apply the traditional
203 : ! negative convention for pvel.
204 : pvliq(i,k+1) = vfall &
205 : * 0.01_r8 & ! cm to meters
206 : * rho(i,k)*gravit ! meters/sec to pascals/sec
207 : #else
208 :
209 : ! newway
210 150493962 : if (rel(i,k) < 40._r8 ) then
211 150493962 : vfall = 2._r8/9._r8 * rhoh2o * gravit * rel(i,k)**2 / eta * 1.e-12_r8 ! micons^2 -> m^2
212 : else
213 0 : vfall = v40 + vslope * (rel(i,k)-r40) ! linear above 40 microns
214 : end if
215 : ! convert the fall speed to pressure units
216 : ! but do not apply the traditional
217 : ! negative convention for pvel.
218 : ! pvliq2(i,k+1) = vfall * rho(i,k)*gravit ! meters/sec to pascals/sec
219 150493962 : pvliq(i,k+1) = vfall * rho(i,k)*gravit ! meters/sec to pascals/sec
220 : #endif
221 : end if
222 : end do
223 : end do
224 :
225 : !-----------------------------------------------------------------------
226 : !--------------------- ice fall velocity -------------------------------
227 : !-----------------------------------------------------------------------
228 :
229 675662904 : pvice(:ncol,:) = 0._r8
230 :
231 : if (stokes) then
232 :
233 : !-----------------------------------------------------------------------
234 : !--------------------- stokes terminal velocity < 40 microns -----------
235 : !-----------------------------------------------------------------------
236 :
237 : ! get effective radius
238 1495368 : call reitab(ncol, t, rei)
239 :
240 40374936 : do k = 1,pver
241 650693736 : do i = 1,ncol
242 649198368 : if (cloud(i,k) > 0._r8 .and. cldice(i,k) > 0._r8) then
243 135492597 : if (rei(i,k) < 40._r8 ) then
244 37825727 : vfall = 2._r8/9._r8 * rhoh2o * gravit * rei(i,k)**2 / eta * 1.e-12_r8 ! micons^2 -> m^2
245 37825727 : vfall = vfall * cldsed_ice_stokes_fac
246 : else
247 97666870 : vfall = v40 + vslope * (rei(i,k)-r40) ! linear above 40 microns
248 : end if
249 :
250 : ! convert the fall speed to pressure units, but do not apply the traditional
251 : ! negative convention for pvel.
252 135492597 : pvice(i,k+1) = vfall * rho(i,k)*gravit ! meters/sec to pascals/sec
253 : end if
254 : end do
255 : end do
256 :
257 : else
258 :
259 : !-----------------------------------------------------------------------
260 : !--------------------- McFarquhar and Heymsfield > icritc --------------
261 : !-----------------------------------------------------------------------
262 :
263 : ! lower bound for iciwc
264 :
265 : cc = 128.64_r8
266 : bc = 53.242_r8
267 : ac = 5.4795_r8
268 : lbound = (-bc + sqrt(bc*bc-4*ac*cc))/(2*ac)
269 : lbound = 10._r8**lbound
270 :
271 : do k = 1,pver
272 : do i = 1,ncol
273 : if (cloud(i,k) > 0._r8 .and. cldice(i,k) > 0._r8) then
274 :
275 : ! compute the in-cloud ice concentration (kg/kg)
276 : icice = cldice(i,k) / cloud(i,k)
277 :
278 : ! compute the ice water content in g/m3
279 : iciwc = icice * rho(i,k) * 1.e3_r8
280 :
281 : ! set the fall velocity (cm/s) to depend on the ice water content in g/m3,
282 : if (iciwc > lbound) then ! need this because of log10
283 : logiwc = log10(iciwc)
284 : ! Median -
285 : vfall = 128.64_r8 + 53.242_r8*logiwc + 5.4795_r8*logiwc**2
286 : ! Average -
287 : !!$ vfall = 122.63 + 44.111*logiwc + 4.2144*logiwc**2
288 : else
289 : vfall = 0._r8
290 : end if
291 :
292 : ! set ice velocity to 1 cm/s if ice mixing ratio < icritc, ramp to value
293 : ! calculated above at 2*icritc
294 : if (icice <= icritc) then
295 : vfall = vice_small
296 : else if(icice < 2*icritc) then
297 : icefac = (icice-icritc)/icritc
298 : vfall = vice_small * (1._r8-icefac) + vfall * icefac
299 : end if
300 :
301 : ! bound the terminal velocity of ice particles at high concentration
302 : vfall = min(100.0_r8, vfall)
303 :
304 : ! convert the fall speed to pressure units, but do not apply the traditional
305 : ! negative convention for pvel.
306 : pvice(i,k+1) = vfall &
307 : * 0.01_r8 & ! cm to meters
308 : * rho(i,k)*gravit ! meters/sec to pascals/sec
309 : end if
310 : end do
311 : end do
312 :
313 : end if
314 :
315 1495368 : end subroutine cld_sediment_vel
316 :
317 :
318 : !===============================================================================
319 1495368 : subroutine cld_sediment_tend (ncol, dtime , &
320 : pint , pmid , pdel , t , &
321 : cloud , cldliq , cldice , pvliq , pvice , &
322 : liqtend, icetend, wvtend , htend , sfliq , sfice )
323 :
324 : !----------------------------------------------------------------------
325 : ! Apply Cloud Particle Gravitational Sedimentation to Condensate
326 : !----------------------------------------------------------------------
327 :
328 :
329 : ! Arguments
330 : integer, intent(in) :: ncol ! number of colums to process
331 :
332 : real(r8), intent(in) :: dtime ! time step
333 : real(r8), intent(in) :: pint (pcols,pverp) ! interfaces pressure (Pa)
334 : real(r8), intent(in) :: pmid (pcols,pver) ! midpoint pressures (Pa)
335 : real(r8), intent(in) :: pdel (pcols,pver) ! pressure diff across layer (Pa)
336 : real(r8), intent(in) :: cloud (pcols,pver) ! cloud fraction (fraction)
337 : real(r8), intent(in) :: t (pcols,pver) ! temperature (K)
338 : real(r8), intent(in) :: cldliq(pcols,pver) ! cloud liquid water (kg/kg)
339 : real(r8), intent(in) :: cldice(pcols,pver) ! cloud ice water (kg/kg)
340 : real(r8), intent(in) :: pvliq (pcols,pverp) ! vertical velocity of liquid drops (Pa/s)
341 : real(r8), intent(in) :: pvice (pcols,pverp) ! vertical velocity of ice particles (Pa/s)
342 : ! -> note that pvel is at the interfaces (loss from cell is based on pvel(k+1))
343 :
344 : real(r8), intent(out) :: liqtend(pcols,pver) ! liquid condensate tend
345 : real(r8), intent(out) :: icetend(pcols,pver) ! ice condensate tend
346 : real(r8), intent(out) :: wvtend (pcols,pver) ! water vapor tend
347 : real(r8), intent(out) :: htend (pcols,pver) ! heating rate
348 : real(r8), intent(out) :: sfliq (pcols) ! surface flux of liquid (rain, kg/m/s)
349 : real(r8), intent(out) :: sfice (pcols) ! surface flux of ice (snow, kg/m/s)
350 :
351 : ! Local variables
352 : real(r8) :: fxliq(pcols,pverp) ! fluxes at the interfaces, liquid (positive = down)
353 : real(r8) :: fxice(pcols,pverp) ! fluxes at the interfaces, ice (positive = down)
354 : real(r8) :: cldab(pcols) ! cloud in layer above
355 : real(r8) :: evapliq ! evaporation of cloud liquid into environment
356 : real(r8) :: evapice ! evaporation of cloud ice into environment
357 : real(r8) :: cldovrl ! cloud overlap factor
358 :
359 : integer :: i,k
360 : !----------------------------------------------------------------------
361 :
362 : ! initialize variables
363 677158272 : fxliq (:ncol,:) = 0._r8 ! flux at interfaces (liquid)
364 675662904 : fxice (:ncol,:) = 0._r8 ! flux at interfaces (ice)
365 650693736 : liqtend(:ncol,:) = 0._r8 ! condensate tend (liquid)
366 650693736 : icetend(:ncol,:) = 0._r8 ! condensate tend (ice)
367 650693736 : wvtend(:ncol,:) = 0._r8 ! environmental moistening
368 650693736 : htend(:ncol,:) = 0._r8 ! evaporative cooling
369 24969168 : sfliq(:ncol) = 0._r8 ! condensate sedimentation flux out bot of column (liquid)
370 24969168 : sfice(:ncol) = 0._r8 ! condensate sedimentation flux out bot of column (ice)
371 :
372 : ! fluxes at interior points
373 1495368 : call getflx(ncol, pint, cldliq, pvliq, dtime, fxliq)
374 1495368 : call getflx(ncol, pint, cldice, pvice, dtime, fxice)
375 :
376 : ! calculate fluxes at boundaries
377 24969168 : do i = 1,ncol
378 23473800 : fxliq(i,1) = 0._r8
379 23473800 : fxice(i,1) = 0._r8
380 : ! surface flux by upstream scheme
381 23473800 : fxliq(i,pverp) = cldliq(i,pver) * pvliq(i,pverp) * dtime
382 24969168 : fxice(i,pverp) = cldice(i,pver) * pvice(i,pverp) * dtime
383 : end do
384 :
385 : ! filter out any negative fluxes from the getflx routine
386 : ! (typical fluxes are of order > 1e-3 when clouds are present)
387 38879568 : do k = 2,pver
388 624229200 : fxliq(:ncol,k) = max(0._r8, fxliq(:ncol,k))
389 625724568 : fxice(:ncol,k) = max(0._r8, fxice(:ncol,k))
390 : end do
391 :
392 : ! Limit the flux out of the bottom of each cell to the water content in each phase.
393 : ! Apply mxsedfac to prevent generating very small negative cloud water/ice
394 : ! NOTE, REMOVED CLOUD FACTOR FROM AVAILABLE WATER. ALL CLOUD WATER IS IN CLOUDS.
395 : ! ***Should we include the flux in the top, to allow for thin surface layers?
396 : ! ***Requires simple treatment of cloud overlap, already included below.
397 40374936 : do k = 1,pver
398 650693736 : do i = 1,ncol
399 610318800 : fxliq(i,k+1) = min( fxliq(i,k+1), mxsedfac * cldliq(i,k) * pdel(i,k) )
400 649198368 : fxice(i,k+1) = min( fxice(i,k+1), mxsedfac * cldice(i,k) * pdel(i,k) )
401 : !!$ fxliq(i,k+1) = min( fxliq(i,k+1), cldliq(i,k) * pdel(i,k) + fxliq(i,k))
402 : !!$ fxice(i,k+1) = min( fxice(i,k+1), cldice(i,k) * pdel(i,k) + fxice(i,k))
403 : !!$ fxliq(i,k+1) = min( fxliq(i,k+1), cloud(i,k) * cldliq(i,k) * pdel(i,k) )
404 : !!$ fxice(i,k+1) = min( fxice(i,k+1), cloud(i,k) * cldice(i,k) * pdel(i,k) )
405 : end do
406 : end do
407 :
408 : ! Now calculate the tendencies assuming that condensate evaporates when
409 : ! it falls into environment, and does not when it falls into cloud.
410 : ! All flux out of the layer comes from the cloudy part.
411 : ! Assume maximum overlap for stratiform clouds
412 : ! if cloud above < cloud, all water falls into cloud below
413 : ! if cloud above > cloud, water split between cloud and environment
414 40374936 : do k = 1,pver
415 649198368 : cldab(:ncol) = 0._r8
416 650693736 : do i = 1,ncol
417 : ! cloud overlap cloud factor
418 610318800 : cldovrl = min( cloud(i,k) / (cldab(i)+.0001_r8), 1._r8 )
419 610318800 : cldab(i) = cloud(i,k)
420 : ! evaporation into environment cause moistening and cooling
421 610318800 : evapliq = fxliq(i,k) * (1._r8-cldovrl) / (dtime * pdel(i,k)) ! into env (kg/kg/s)
422 610318800 : evapice = fxice(i,k) * (1._r8-cldovrl) / (dtime * pdel(i,k)) ! into env (kg/kg/s)
423 610318800 : wvtend(i,k) = evapliq + evapice ! evaporation into environment (kg/kg/s)
424 610318800 : htend (i,k) = -latvap*evapliq -(latvap+latice)*evapice ! evaporation (W/kg)
425 : ! net flux into cloud changes cloud liquid/ice (all flux is out of cloud)
426 610318800 : liqtend(i,k) = (fxliq(i,k)*cldovrl - fxliq(i,k+1)) / (dtime * pdel(i,k))
427 649198368 : icetend(i,k) = (fxice(i,k)*cldovrl - fxice(i,k+1)) / (dtime * pdel(i,k))
428 : end do
429 : end do
430 :
431 : ! convert flux out the bottom to mass units Pa -> kg/m2/s
432 24969168 : sfliq(:ncol) = fxliq(:ncol,pverp) / (dtime*gravit)
433 24969168 : sfice(:ncol) = fxice(:ncol,pverp) / (dtime*gravit)
434 :
435 1495368 : return
436 : end subroutine cld_sediment_tend
437 :
438 : !===============================================================================
439 2990736 : subroutine getflx(ncol, xw, phi, vel, deltat, flux)
440 :
441 : !.....xw1.......xw2.......xw3.......xw4.......xw5.......xw6
442 : !....psiw1.....psiw2.....psiw3.....psiw4.....psiw5.....psiw6
443 : !....velw1.....velw2.....velw3.....velw4.....velw5.....velw6
444 : !.........phi1......phi2.......phi3.....phi4.......phi5.......
445 :
446 :
447 :
448 : integer, intent(in) :: ncol ! number of colums to process
449 :
450 : integer i
451 : integer k
452 :
453 : real (r8), intent(in) :: vel(pcols,pverp)
454 : real (r8) flux(pcols,pverp)
455 : real (r8) xw(pcols,pverp)
456 : real (r8) psi(pcols,pverp)
457 : real (r8), intent(in) :: phi(pcols,pverp-1)
458 : real (r8) fdot(pcols,pverp)
459 : real (r8) xx(pcols)
460 : real (r8) fxdot(pcols)
461 : real (r8) fxdd(pcols)
462 :
463 : real (r8) psistar(pcols)
464 : real (r8) deltat
465 :
466 : real (r8) xxk(pcols,pver)
467 :
468 49938336 : do i = 1,ncol
469 : ! integral of phi
470 46947600 : psi(i,1) = 0._r8
471 : ! fluxes at boundaries
472 46947600 : flux(i,1) = 0._r8
473 49938336 : flux(i,pverp) = 0._r8
474 : end do
475 :
476 : ! integral function
477 80749872 : do k = 2,pverp
478 1301387472 : do i = 1,ncol
479 1298396736 : psi(i,k) = phi(i,k-1)*(xw(i,k)-xw(i,k-1)) + psi(i,k-1)
480 : end do
481 : end do
482 :
483 :
484 : ! calculate the derivatives for the interpolating polynomial
485 2990736 : call cfdotmc_pro (ncol, xw, psi, fdot)
486 :
487 : ! NEW WAY
488 : ! calculate fluxes at interior pts
489 77759136 : do k = 2,pver
490 1251449136 : do i = 1,ncol
491 1248458400 : xxk(i,k) = xw(i,k)-vel(i,k)*deltat
492 : end do
493 : end do
494 77759136 : do k = 2,pver
495 74768400 : call cfint2(ncol, xw, psi, fdot, xxk(1,k), fxdot, fxdd, psistar)
496 1251449136 : do i = 1,ncol
497 1248458400 : flux(i,k) = (psi(i,k)-psistar(i))
498 : end do
499 : end do
500 :
501 :
502 2990736 : return
503 : end subroutine getflx
504 :
505 :
506 :
507 : !##############################################################################
508 :
509 74768400 : subroutine cfint2 (ncol, x, f, fdot, xin, fxdot, fxdd, psistar)
510 :
511 :
512 :
513 : ! input
514 : integer ncol ! number of colums to process
515 :
516 : real (r8) x(pcols, pverp)
517 : real (r8) f(pcols, pverp)
518 : real (r8) fdot(pcols, pverp)
519 : real (r8) xin(pcols)
520 :
521 : ! output
522 : real (r8) fxdot(pcols)
523 : real (r8) fxdd(pcols)
524 : real (r8) psistar(pcols)
525 :
526 : integer i
527 : integer k
528 : integer intz(pcols)
529 : real (r8) dx
530 : real (r8) s
531 : real (r8) c2
532 : real (r8) c3
533 : real (r8) xx
534 : real (r8) xinf
535 : real (r8) psi1, psi2, psi3, psim
536 : real (r8) cfint
537 : real (r8) cfnew
538 : real (r8) xins(pcols)
539 :
540 : ! the minmod function
541 : real (r8) a, b, c
542 : real (r8) minmod
543 : real (r8) medan
544 : logical found_error
545 :
546 : minmod(a,b) = 0.5_r8*(sign(1._r8,a) + sign(1._r8,b))*min(abs(a),abs(b))
547 : medan(a,b,c) = a + minmod(b-a,c-a)
548 :
549 1248458400 : do i = 1,ncol
550 1173690000 : xins(i) = medan(x(i,1), xin(i), x(i,pverp))
551 1248458400 : intz(i) = 0
552 : end do
553 :
554 : ! first find the interval
555 2018746800 : do k = 1,pverp-1
556 32534686800 : do i = 1,ncol
557 32459918400 : if ((xins(i)-x(i,k))*(x(i,k+1)-xins(i)).ge.0) then
558 2056902817 : intz(i) = k
559 : endif
560 : end do
561 : end do
562 :
563 74768400 : found_error=.false.
564 1248458400 : do i = 1,ncol
565 1248458400 : if (intz(i).eq.0._r8) found_error=.true.
566 : end do
567 74768400 : if(found_error) then
568 0 : do i = 1,ncol
569 0 : if (intz(i).eq.0._r8) then
570 0 : write(iulog,*) ' interval was not found for col i ', i
571 0 : call endrun('CFINT2')
572 : endif
573 : end do
574 : endif
575 :
576 : ! now interpolate
577 1248458400 : do i = 1,ncol
578 1173690000 : k = intz(i)
579 1173690000 : dx = (x(i,k+1)-x(i,k))
580 1173690000 : s = (f(i,k+1)-f(i,k))/dx
581 1173690000 : c2 = (3*s-2*fdot(i,k)-fdot(i,k+1))/dx
582 1173690000 : c3 = (fdot(i,k)+fdot(i,k+1)-2*s)/dx**2
583 1173690000 : xx = (xins(i)-x(i,k))
584 1173690000 : fxdot(i) = (3*c3*xx + 2*c2)*xx + fdot(i,k)
585 1173690000 : fxdd(i) = 6*c3*xx + 2*c2
586 1173690000 : cfint = ((c3*xx + c2)*xx + fdot(i,k))*xx + f(i,k)
587 :
588 : ! limit the interpolant
589 1173690000 : psi1 = f(i,k)+(f(i,k+1)-f(i,k))*xx/dx
590 1173690000 : if (k.eq.1) then
591 0 : psi2 = f(i,1)
592 : else
593 1173690000 : psi2 = f(i,k) + (f(i,k)-f(i,k-1))*xx/(x(i,k)-x(i,k-1))
594 : endif
595 1173690000 : if (k+1.eq.pverp) then
596 23602375 : psi3 = f(i,pverp)
597 : else
598 1150087625 : psi3 = f(i,k+1) - (f(i,k+2)-f(i,k+1))*(dx-xx)/(x(i,k+2)-x(i,k+1))
599 : endif
600 1173690000 : psim = medan(psi1, psi2, psi3)
601 1173690000 : cfnew = medan(cfint, psi1, psim)
602 : if (abs(cfnew-cfint)/(abs(cfnew)+abs(cfint)+1.e-36_r8) .gt..03_r8) then
603 : ! CHANGE THIS BACK LATER!!!
604 : ! $ .gt..1) then
605 :
606 :
607 : ! UNCOMMENT THIS LATER!!!
608 : ! write(iulog,*) ' cfint2 limiting important ', cfint, cfnew
609 :
610 :
611 : endif
612 1248458400 : psistar(i) = cfnew
613 : end do
614 :
615 74768400 : return
616 : end subroutine cfint2
617 :
618 :
619 :
620 : !##############################################################################
621 :
622 2990736 : subroutine cfdotmc_pro (ncol, x, f, fdot)
623 :
624 : ! prototype version; eventually replace with final SPITFIRE scheme
625 :
626 : ! calculate the derivative for the interpolating polynomial
627 : ! multi column version
628 :
629 :
630 :
631 : ! input
632 : integer ncol ! number of colums to process
633 :
634 : real (r8) x(pcols, pverp)
635 : real (r8) f(pcols, pverp)
636 : ! output
637 : real (r8) fdot(pcols, pverp) ! derivative at nodes
638 :
639 : ! assumed variable distribution
640 : ! x1.......x2.......x3.......x4.......x5.......x6 1,pverp points
641 : ! f1.......f2.......f3.......f4.......f5.......f6 1,pverp points
642 : ! ...sh1.......sh2......sh3......sh4......sh5.... 1,pver points
643 : ! .........d2.......d3.......d4.......d5......... 2,pver points
644 : ! .........s2.......s3.......s4.......s5......... 2,pver points
645 : ! .............dh2......dh3......dh4............. 2,pver-1 points
646 : ! .............eh2......eh3......eh4............. 2,pver-1 points
647 : ! ..................e3.......e4.................. 3,pver-1 points
648 : ! .................ppl3......ppl4................ 3,pver-1 points
649 : ! .................ppr3......ppr4................ 3,pver-1 points
650 : ! .................t3........t4.................. 3,pver-1 points
651 : ! ................fdot3.....fdot4................ 3,pver-1 points
652 :
653 :
654 : ! work variables
655 :
656 :
657 : integer i
658 : integer k
659 :
660 : real (r8) a ! work var
661 : real (r8) b ! work var
662 : real (r8) c ! work var
663 : real (r8) s(pcols,pverp) ! first divided differences at nodes
664 : real (r8) sh(pcols,pverp) ! first divided differences between nodes
665 : real (r8) d(pcols,pverp) ! second divided differences at nodes
666 : real (r8) dh(pcols,pverp) ! second divided differences between nodes
667 : real (r8) e(pcols,pverp) ! third divided differences at nodes
668 : real (r8) eh(pcols,pverp) ! third divided differences between nodes
669 : real (r8) pp ! p prime
670 : real (r8) ppl(pcols,pverp) ! p prime on left
671 : real (r8) ppr(pcols,pverp) ! p prime on right
672 : real (r8) qpl
673 : real (r8) qpr
674 : real (r8) ttt
675 : real (r8) t
676 : real (r8) tmin
677 : real (r8) tmax
678 : real (r8) delxh(pcols,pverp)
679 :
680 :
681 : ! the minmod function
682 : real (r8) minmod
683 : real (r8) medan
684 : minmod(a,b) = 0.5_r8*(sign(1._r8,a) + sign(1._r8,b))*min(abs(a),abs(b))
685 : medan(a,b,c) = a + minmod(b-a,c-a)
686 :
687 80749872 : do k = 1,pver
688 :
689 :
690 : ! first divided differences between nodes
691 1298396736 : do i = 1, ncol
692 1220637600 : delxh(i,k) = (x(i,k+1)-x(i,k))
693 1298396736 : sh(i,k) = (f(i,k+1)-f(i,k))/delxh(i,k)
694 : end do
695 :
696 : ! first and second divided differences at nodes
697 80749872 : if (k.ge.2) then
698 1248458400 : do i = 1,ncol
699 1173690000 : d(i,k) = (sh(i,k)-sh(i,k-1))/(x(i,k+1)-x(i,k-1))
700 1248458400 : s(i,k) = minmod(sh(i,k),sh(i,k-1))
701 : end do
702 : endif
703 : end do
704 :
705 : ! second and third divided diffs between nodes
706 74768400 : do k = 2,pver-1
707 1201510800 : do i = 1, ncol
708 1126742400 : eh(i,k) = (d(i,k+1)-d(i,k))/(x(i,k+2)-x(i,k-1))
709 1198520064 : dh(i,k) = minmod(d(i,k),d(i,k+1))
710 : end do
711 : end do
712 :
713 : ! treat the boundaries
714 49938336 : do i = 1,ncol
715 46947600 : e(i,2) = eh(i,2)
716 46947600 : e(i,pver) = eh(i,pver-1)
717 : ! outside level
718 : fdot(i,1) = sh(i,1) - d(i,2)*delxh(i,1) &
719 46947600 : - eh(i,2)*delxh(i,1)*(x(i,1)-x(i,3))
720 46947600 : fdot(i,1) = minmod(fdot(i,1),3*sh(i,1))
721 : fdot(i,pverp) = sh(i,pver) + d(i,pver)*delxh(i,pver) &
722 46947600 : + eh(i,pver-1)*delxh(i,pver)*(x(i,pverp)-x(i,pver-1))
723 46947600 : fdot(i,pverp) = minmod(fdot(i,pverp),3*sh(i,pver))
724 : ! one in from boundary
725 46947600 : fdot(i,2) = sh(i,1) + d(i,2)*delxh(i,1) - eh(i,2)*delxh(i,1)*delxh(i,2)
726 46947600 : fdot(i,2) = minmod(fdot(i,2),3*s(i,2))
727 : fdot(i,pver) = sh(i,pver) - d(i,pver)*delxh(i,pver) &
728 46947600 : - eh(i,pver-1)*delxh(i,pver)*delxh(i,pver-1)
729 49938336 : fdot(i,pver) = minmod(fdot(i,pver),3*s(i,pver))
730 : end do
731 :
732 :
733 71777664 : do k = 3,pver-1
734 1151572464 : do i = 1,ncol
735 1148581728 : e(i,k) = minmod(eh(i,k),eh(i,k-1))
736 : end do
737 : end do
738 :
739 :
740 :
741 71777664 : do k = 3,pver-1
742 :
743 1151572464 : do i = 1,ncol
744 :
745 : ! p prime at k-0.5
746 1079794800 : ppl(i,k)=sh(i,k-1) + dh(i,k-1)*delxh(i,k-1)
747 : ! p prime at k+0.5
748 1079794800 : ppr(i,k)=sh(i,k) - dh(i,k) *delxh(i,k)
749 :
750 1079794800 : t = minmod(ppl(i,k),ppr(i,k))
751 :
752 : ! derivate from parabola thru f(i,k-1), f(i,k), and f(i,k+1)
753 1079794800 : pp = sh(i,k-1) + d(i,k)*delxh(i,k-1)
754 :
755 : ! quartic estimate of fdot
756 : fdot(i,k) = pp &
757 : - delxh(i,k-1)*delxh(i,k) &
758 1079794800 : *( eh(i,k-1)*(x(i,k+2)-x(i,k )) &
759 1079794800 : + eh(i,k )*(x(i,k )-x(i,k-2)) &
760 2159589600 : )/(x(i,k+2)-x(i,k-2))
761 :
762 : ! now limit it
763 : qpl = sh(i,k-1) &
764 : + delxh(i,k-1)*minmod(d(i,k-1)+e(i,k-1)*(x(i,k)-x(i,k-2)), &
765 1079794800 : d(i,k) -e(i,k)*delxh(i,k))
766 : qpr = sh(i,k) &
767 : + delxh(i,k )*minmod(d(i,k) +e(i,k)*delxh(i,k-1), &
768 1079794800 : d(i,k+1)+e(i,k+1)*(x(i,k)-x(i,k+2)))
769 :
770 1079794800 : fdot(i,k) = medan(fdot(i,k), qpl, qpr)
771 :
772 1079794800 : ttt = minmod(qpl, qpr)
773 1079794800 : tmin = min(0._r8,3*s(i,k),1.5_r8*t,ttt)
774 1079794800 : tmax = max(0._r8,3*s(i,k),1.5_r8*t,ttt)
775 :
776 1148581728 : fdot(i,k) = fdot(i,k) + minmod(tmin-fdot(i,k), tmax-fdot(i,k))
777 :
778 : end do
779 :
780 : end do
781 :
782 2990736 : return
783 : end subroutine cfdotmc_pro
784 : end module pkg_cld_sediment
|
