Line data Source code
1 :
2 : module diffusion_solver
3 :
4 : !------------------------------------------------------------------------------------ !
5 : ! Module to solve vertical diffusion equations using a tri-diagonal solver. !
6 : ! The module will also apply countergradient fluxes, and apply molecular !
7 : ! diffusion for constituents. !
8 : ! !
9 : ! Public interfaces : !
10 : ! init_vdiff initializes time independent coefficients !
11 : ! compute_vdiff solves diffusion equations !
12 : ! vdiff_selector type for storing fields selected to be diffused !
13 : ! vdiff_select selects fields to be diffused !
14 : ! operator(.not.) extends .not. to operate on type vdiff_selector !
15 : ! any provides functionality of intrinsic any for type vdiff_selector !
16 : ! !
17 : !------------------------------------ Code History ---------------------------------- !
18 : ! Initial subroutines : B. Boville and others, 1991-2004 !
19 : ! Modularization : J. McCaa, September 2004 !
20 : ! Most Recent Code : Sungsu Park, Aug. 2006, Dec. 2008, Jan. 2010. !
21 : !------------------------------------------------------------------------------------ !
22 :
23 : implicit none
24 : private
25 : save
26 :
27 : integer, parameter :: r8 = selected_real_kind(12) ! 8 byte real
28 :
29 : ! ----------------- !
30 : ! Public interfaces !
31 : ! ----------------- !
32 :
33 : public init_vdiff ! Initialization
34 : public new_fieldlist_vdiff ! Returns an empty fieldlist
35 : public compute_vdiff ! Full routine
36 : public vdiff_selector ! Type for storing fields selected to be diffused
37 : public vdiff_select ! Selects fields to be diffused
38 : public operator(.not.) ! Extends .not. to operate on type vdiff_selector
39 : public any ! Provides functionality of intrinsic any for type vdiff_selector
40 :
41 : ! Below stores logical array of fields to be diffused
42 :
43 : type vdiff_selector
44 : private
45 : logical, allocatable, dimension(:) :: fields
46 : end type vdiff_selector
47 :
48 : ! Below extends .not. to operate on type vdiff_selector
49 :
50 : interface operator(.not.)
51 : module procedure not
52 : end interface
53 :
54 : ! Below provides functionality of intrinsic any for type vdiff_selector
55 :
56 : interface any
57 : module procedure my_any
58 : end interface
59 :
60 : ! ------------ !
61 : ! Private data !
62 : ! ------------ !
63 :
64 : ! Unit number for log output
65 : integer :: iulog = -1
66 :
67 : real(r8), private :: cpair ! Specific heat of dry air
68 : real(r8), private :: gravit ! Acceleration due to gravity
69 : real(r8), private :: rair ! Gas constant for dry air
70 :
71 : logical, private :: do_iss ! Use implicit turbulent surface stress computation
72 :
73 : ! Parameters used for Turbulent Mountain Stress
74 :
75 : real(r8), parameter :: z0fac = 0.025_r8 ! Factor determining z_0 from orographic standard deviation
76 : real(r8), parameter :: z0max = 100._r8 ! Max value of z_0 for orography
77 : real(r8), parameter :: horomin = 10._r8 ! Min value of subgrid orographic height for mountain stress
78 : real(r8), parameter :: dv2min = 0.01_r8 ! Minimum shear squared
79 :
80 : logical :: am_correction ! logical switch for AM correction
81 :
82 : contains
83 :
84 : ! =============================================================================== !
85 : ! !
86 : ! =============================================================================== !
87 :
88 1536 : subroutine init_vdiff( kind, iulog_in, rair_in, cpair_in, gravit_in, do_iss_in, &
89 1536 : am_correction_in, errstring )
90 :
91 : integer, intent(in) :: kind ! Kind used for reals
92 : integer, intent(in) :: iulog_in ! Unit number for log output.
93 : real(r8), intent(in) :: rair_in ! Input gas constant for dry air
94 : real(r8), intent(in) :: cpair_in ! Input heat capacity for dry air
95 : real(r8), intent(in) :: gravit_in ! Input gravitational acceleration
96 : logical, intent(in) :: do_iss_in ! Input ISS flag
97 : logical, intent(in) :: am_correction_in! for angular momentum conservation
98 : character(128), intent(out) :: errstring ! Output status
99 :
100 1536 : errstring = ''
101 1536 : iulog = iulog_in
102 1536 : if( kind .ne. r8 ) then
103 0 : write(iulog,*) 'KIND of reals passed to init_vdiff -- exiting.'
104 0 : errstring = 'init_vdiff'
105 0 : return
106 : endif
107 :
108 1536 : rair = rair_in
109 1536 : cpair = cpair_in
110 1536 : gravit = gravit_in
111 1536 : do_iss = do_iss_in
112 1536 : am_correction = am_correction_in
113 :
114 : end subroutine init_vdiff
115 :
116 : ! =============================================================================== !
117 : ! !
118 : ! =============================================================================== !
119 :
120 4608 : type(vdiff_selector) pure function new_fieldlist_vdiff(ncnst)
121 :
122 : integer, intent(in) :: ncnst ! Number of constituents
123 :
124 13824 : allocate( new_fieldlist_vdiff%fields( 3 + ncnst ) )
125 207360 : new_fieldlist_vdiff%fields = .false.
126 :
127 4608 : end function new_fieldlist_vdiff
128 :
129 : ! =============================================================================== !
130 : ! !
131 : ! =============================================================================== !
132 :
133 1489176 : subroutine compute_vdiff( lchnk , &
134 1489176 : pcols , pver , ncnst , ncol , tint , &
135 2978352 : p , t , rhoi , ztodt , taux , &
136 1489176 : tauy , shflx , cflx , &
137 1489176 : kvh , kvm , kvq , cgs , cgh , &
138 1489176 : zi , ksrftms , dragblj , &
139 1489176 : qmincg , fieldlist , fieldlistm , &
140 1489176 : u , v , q , dse , &
141 0 : tautmsx , tautmsy , dtk , topflx , errstring , &
142 1489176 : tauresx , tauresy , itaures , cpairv , dse_top, &
143 : do_molec_diff , use_temperature_molec_diff, vd_lu_qdecomp, &
144 5956704 : ubc_mmr, ubc_flux, kvt, pmid, &
145 4467528 : cnst_mw, cnst_fixed_ubc, cnst_fixed_ubflx, nbot_molec, &
146 1489176 : kq_scal, mw_fac)
147 :
148 : !-------------------------------------------------------------------------- !
149 : ! Driver routine to compute vertical diffusion of momentum, moisture, trace !
150 : ! constituents and dry static energy. The new temperature is computed from !
151 : ! the diffused dry static energy. !
152 : ! Turbulent diffusivities and boundary layer nonlocal transport terms are !
153 : ! obtained from the turbulence module. !
154 : !-------------------------------------------------------------------------- !
155 :
156 : ! Used for CAM debugging.
157 : ! use phys_debug_util, only : phys_debug_col
158 : ! use time_manager, only : is_first_step, get_nstep
159 :
160 : use coords_1d, only: Coords1D
161 : use linear_1d_operators, only : BoundaryType, BoundaryFixedLayer, &
162 : BoundaryData, BoundaryFlux, TriDiagDecomp
163 : use vdiff_lu_solver, only : fin_vol_lu_decomp
164 : use vertical_diffusion_solver, only : fin_vol_solve
165 : use beljaars_drag_cam, only : do_beljaars
166 : ! FIXME: This should not be needed
167 : use air_composition, only: rairv
168 :
169 : use phys_control, only : phys_getopts
170 :
171 : ! Modification : Ideally, we should diffuse 'liquid-ice static energy' (sl), not the dry static energy.
172 : ! Also, vertical diffusion of cloud droplet number concentration and aerosol number
173 : ! concentration should be done very carefully in the future version.
174 :
175 : ! --------------- !
176 : ! Input Arguments !
177 : ! --------------- !
178 :
179 : integer, intent(in) :: lchnk
180 : integer, intent(in) :: pcols
181 : integer, intent(in) :: pver
182 : integer, intent(in) :: ncnst
183 : integer, intent(in) :: ncol ! Number of atmospheric columns
184 : integer, intent(in) :: itaures ! Indicator determining whether 'tauresx,tauresy'
185 : ! is updated (1) or non-updated (0) in this subroutine.
186 :
187 : type(Coords1D), intent(in) :: p ! Pressure coordinates [ Pa ]
188 : real(r8), intent(in) :: tint(pcols,pver+1) ! Temperature [ K ]
189 : real(r8), intent(in) :: t(pcols,pver) ! Temperature [ K ]
190 : real(r8), intent(in) :: rhoi(pcols,pver+1) ! Density of air at interfaces [ kg/m3 ]
191 : real(r8), intent(in) :: ztodt ! 2 delta-t [ s ]
192 : real(r8), intent(in) :: taux(pcols) ! Surface zonal stress.
193 : ! Input u-momentum per unit time per unit area into the atmosphere [ N/m2 ]
194 : real(r8), intent(in) :: tauy(pcols) ! Surface meridional stress.
195 : ! Input v-momentum per unit time per unit area into the atmosphere [ N/m2 ]
196 : real(r8), intent(in) :: shflx(pcols) ! Surface sensible heat flux [ W/m2 ]
197 : real(r8), intent(in) :: cflx(pcols,ncnst) ! Surface constituent flux [ kg/m2/s ]
198 : real(r8), intent(in) :: zi(pcols,pver+1) ! Interface heights [ m ]
199 : real(r8), intent(in) :: ksrftms(pcols) ! Surface drag coefficient for turbulent mountain stress. > 0. [ kg/s/m2 ]
200 : real(r8), intent(in) :: dragblj(pcols,pver) ! Drag profile from Beljaars SGO form drag > 0. [ 1/s ]
201 : real(r8), intent(in) :: qmincg(ncnst) ! Minimum constituent mixing ratios from cg fluxes
202 : real(r8), intent(in) :: cpairv(pcols,pver) ! Specific heat at constant pressure
203 : real(r8), intent(in) :: kvh(pcols,pver+1) ! Eddy diffusivity for heat [ m2/s ]
204 :
205 : logical, intent(in) :: do_molec_diff ! Flag indicating multiple constituent diffusivities
206 : logical, intent(in) :: use_temperature_molec_diff! Flag indicating that molecular diffusion should apply to temperature, not
207 : ! dry static energy.
208 :
209 : type(vdiff_selector), intent(in) :: fieldlist ! Array of flags selecting which fields to diffuse
210 : type(vdiff_selector), intent(in) :: fieldlistm ! Array of flags selecting which fields for molecular diffusion
211 :
212 : ! Dry static energy top boundary condition.
213 : real(r8), intent(in) :: dse_top(pcols)
214 :
215 : real(r8), intent(in) :: kvm(pcols,pver+1) ! Eddy viscosity ( Eddy diffusivity for momentum ) [ m2/s ]
216 : real(r8), intent(in) :: kvq(pcols,pver+1) ! Eddy diffusivity for constituents
217 : real(r8), intent(in) :: cgs(pcols,pver+1) ! Counter-gradient star [ cg/flux ]
218 : real(r8), intent(in) :: cgh(pcols,pver+1) ! Counter-gradient term for heat
219 :
220 : ! ---------------------- !
221 : ! Input-Output Arguments !
222 : ! ---------------------- !
223 :
224 : real(r8), intent(inout) :: u(pcols,pver) ! U wind. This input is the 'raw' input wind to
225 : ! PBL scheme without iterative provisional update. [ m/s ]
226 : real(r8), intent(inout) :: v(pcols,pver) ! V wind. This input is the 'raw' input wind to PBL scheme
227 : ! without iterative provisional update. [ m/s ]
228 : real(r8), intent(inout) :: q(pcols,pver,ncnst) ! Moisture and trace constituents [ kg/kg, #/kg ? ]
229 : real(r8), intent(inout) :: dse(pcols,pver) ! Dry static energy [ J/kg ]
230 :
231 : real(r8), intent(inout) :: tauresx(pcols) ! Input : Reserved surface stress at previous time step
232 : real(r8), intent(inout) :: tauresy(pcols) ! Output : Reserved surface stress at current time step
233 :
234 : ! ---------------- !
235 : ! Output Arguments !
236 : ! ---------------- !
237 :
238 : real(r8), intent(out) :: dtk(pcols,pver) ! T tendency from KE dissipation
239 : real(r8), intent(out) :: tautmsx(pcols) ! Implicit zonal turbulent mountain surface stress
240 : ! [ N/m2 = kg m/s /s/m2 ]
241 : real(r8), intent(out) :: tautmsy(pcols) ! Implicit meridional turbulent mountain surface stress
242 : ! [ N/m2 = kg m/s /s/m2 ]
243 : real(r8), intent(out) :: topflx(pcols) ! Molecular heat flux at the top interface
244 : character(128), intent(out) :: errstring ! Output status
245 :
246 : ! ------------------ !
247 : ! Optional Arguments !
248 : ! ------------------ !
249 :
250 : ! The molecular diffusion module will likely change significantly in
251 : ! the future, and this module may directly depend on it after that.
252 : ! Until then, we have these highly specific interfaces hard-coded.
253 :
254 : optional :: vd_lu_qdecomp ! Constituent-dependent molecular diffusivity routine
255 :
256 : interface
257 : function vd_lu_qdecomp( &
258 : pcols , pver , ncol , fixed_ubc , mw , &
259 : kv , kq_scal, mw_facm , dpidz_sq , coords , &
260 : interface_boundary, molec_boundary, &
261 : tint , ztodt , nbot_molec , &
262 : lchnk , t , m , no_molec_decomp) result(decomp)
263 : import
264 : integer, intent(in) :: pcols
265 : integer, intent(in) :: pver
266 : integer, intent(in) :: ncol
267 : integer, intent(in) :: nbot_molec
268 : logical, intent(in) :: fixed_ubc
269 : real(r8), intent(in) :: kv(pcols,pver+1)
270 : real(r8), intent(in) :: kq_scal(pcols,pver+1)
271 : real(r8), intent(in) :: mw
272 : real(r8), intent(in) :: mw_facm(pcols,pver+1)
273 : real(r8), intent(in) :: dpidz_sq(ncol,pver+1)
274 : type(Coords1D), intent(in) :: coords
275 : type(BoundaryType), intent(in) :: interface_boundary
276 : type(BoundaryType), intent(in) :: molec_boundary
277 : real(r8), intent(in) :: tint(pcols,pver+1)
278 : real(r8), intent(in) :: ztodt
279 : integer, intent(in) :: lchnk
280 : real(r8), intent(in) :: t(pcols,pver)
281 : integer, intent(in) :: m
282 : type(TriDiagDecomp), intent(in) :: no_molec_decomp
283 : type(TriDiagDecomp) :: decomp
284 : end function vd_lu_qdecomp
285 : end interface
286 :
287 : real(r8), intent(in), optional :: ubc_mmr(pcols,ncnst) ! Upper boundary mixing ratios [ kg/kg ]
288 : real(r8), intent(in), optional :: ubc_flux(pcols,ncnst) ! Upper boundary flux [ kg/s/m^2 ]
289 :
290 : real(r8), intent(in), optional :: kvt(pcols,pver+1) ! Kinematic molecular conductivity
291 :
292 : ! FIXME: This input should not be needed (and should not be passed in in vertical_diffusion).
293 : real(r8), intent(in), optional :: pmid(pcols,pver)
294 :
295 : real(r8), intent(in), optional :: cnst_mw(ncnst) ! Molecular weight [ kg/kmole ]
296 : logical, intent(in), optional :: cnst_fixed_ubc(ncnst) ! Whether upper boundary condition is fixed
297 : logical, intent(in), optional :: cnst_fixed_ubflx(ncnst) ! Whether upper boundary flux is a fixed non-zero value
298 :
299 : integer, intent(in), optional :: nbot_molec ! Bottom level where molecular diffusivity is applied
300 :
301 : ! kq_fac*sqrt(T)*m_d/rho for molecular diffusivity
302 : real(r8), intent(in), optional :: kq_scal(pcols,pver+1)
303 : ! Local sqrt(1/M_q + 1/M_d) for each constituent
304 : real(r8), intent(in), optional :: mw_fac(pcols,pver+1,ncnst)
305 :
306 : ! --------------- !
307 : ! Local Variables !
308 : ! --------------- !
309 :
310 : integer :: i, k, m ! Longitude, level, constituent indices
311 2978352 : logical :: lqtst(pcols) ! Adjust vertical profiles
312 :
313 : ! LU decomposition information.
314 1489176 : type(TriDiagDecomp) :: decomp
315 1489176 : type(TriDiagDecomp) :: no_molec_decomp
316 :
317 : ! Square of derivative of pressure with height (on interfaces).
318 2978352 : real(r8) :: dpidz_sq(ncol,pver+1)
319 :
320 : ! Pressure coordinates over the molecular diffusion range only.
321 1489176 : type(Coords1D) :: p_molec
322 :
323 : ! Boundary layer objects
324 1489176 : type(BoundaryType) :: interface_boundary
325 1489176 : type(BoundaryType) :: molec_boundary
326 :
327 2978352 : real(r8) :: tmp1(pcols) ! Temporary storage
328 2978352 : real(r8) :: tmpi1(pcols,pver+1) ! Interface KE dissipation
329 2978352 : real(r8) :: tmpi2(pcols,pver+1) ! dt*(g*rho)**2/dp at interfaces
330 2978352 : real(r8) :: keg_in(pcols,pver) ! KE on entry to subroutine
331 2978352 : real(r8) :: keg_out(pcols,pver) ! KE after U and V dissipation/diffusion
332 2978352 : real(r8) :: rrho(pcols) ! 1./bottom level density
333 :
334 2978352 : real(r8) :: tautotx(pcols) ! Total surface stress ( zonal )
335 2978352 : real(r8) :: tautoty(pcols) ! Total surface stress ( meridional )
336 :
337 2978352 : real(r8) :: dinp_u(pcols,pver+1) ! Vertical difference at interfaces, input u
338 2978352 : real(r8) :: dinp_v(pcols,pver+1) ! Vertical difference at interfaces, input v
339 : real(r8) :: dout_u ! Vertical difference at interfaces, output u
340 : real(r8) :: dout_v ! Vertical difference at interfaces, output v
341 :
342 2978352 : real(r8) :: qtm(pcols,pver) ! Temporary copy of q
343 :
344 2978352 : real(r8) :: ws(pcols) ! Lowest-level wind speed [ m/s ]
345 2978352 : real(r8) :: tau(pcols) ! Turbulent surface stress ( not including mountain stress )
346 2978352 : real(r8) :: ksrfturb(pcols) ! Surface drag coefficient of 'normal' stress. > 0.
347 : ! Virtual mass input per unit time per unit area [ kg/s/m2 ]
348 2978352 : real(r8) :: ksrf(pcols) ! Surface drag coefficient of 'normal' stress +
349 : ! Surface drag coefficient of 'tms' stress. > 0. [ kg/s/m2 ]
350 2978352 : real(r8) :: usum_in(pcols) ! Vertical integral of input u-momentum. Total zonal
351 : ! momentum per unit area in column [ sum of u*dp/g = kg m/s m-2 ]
352 2978352 : real(r8) :: vsum_in(pcols) ! Vertical integral of input v-momentum. Total meridional
353 : ! momentum per unit area in column [ sum of v*dp/g = kg m/s m-2 ]
354 2978352 : real(r8) :: usum_mid(pcols) ! Vertical integral of u-momentum after adding explicit residual stress
355 2978352 : real(r8) :: vsum_mid(pcols) ! Vertical integral of v-momentum after adding explicit residual stress
356 2978352 : real(r8) :: usum_out(pcols) ! Vertical integral of u-momentum after doing implicit diffusion
357 2978352 : real(r8) :: vsum_out(pcols) ! Vertical integral of v-momentum after doing implicit diffusion
358 2978352 : real(r8) :: tauimpx(pcols) ! Actual net stress added at the current step other than mountain stress
359 2978352 : real(r8) :: tauimpy(pcols) ! Actual net stress added at the current step other than mountain stress
360 : real(r8) :: ramda ! dt/timeres [ no unit ]
361 :
362 2978352 : real(r8) :: taubljx(pcols) ! recomputed explicit/residual beljaars stress
363 2978352 : real(r8) :: taubljy(pcols) ! recomputed explicit/residual beljaars stress
364 :
365 : ! Rate at which external (surface) stress damps wind speeds (1/s).
366 2978352 : real(r8) :: tau_damp_rate(ncol, pver)
367 :
368 : ! Combined molecular and eddy diffusion.
369 2978352 : real(r8) :: kv_total(pcols,pver+1)
370 :
371 : logical :: use_spcam
372 :
373 : !--------------------------------
374 : ! Variables needed for WACCM-X
375 : !--------------------------------
376 2978352 : real(r8) :: ttemp(ncol,pver) ! temporary temperature array
377 1489176 : real(r8) :: ttemp0(ncol,pver) ! temporary temperature array
378 :
379 : ! ------------------------------------------------ !
380 : ! Parameters for implicit surface stress treatment !
381 : ! ------------------------------------------------ !
382 :
383 : real(r8), parameter :: wsmin = 1._r8 ! Minimum sfc wind speed for estimating frictional
384 : ! transfer velocity ksrf. [ m/s ]
385 : real(r8), parameter :: ksrfmin = 1.e-4_r8 ! Minimum surface drag coefficient [ kg/s/m^2 ]
386 : real(r8), parameter :: timeres = 7200._r8 ! Relaxation time scale of residual stress ( >= dt ) [ s ]
387 :
388 : ! ----------------------- !
389 : ! Main Computation Begins !
390 : ! ----------------------- !
391 :
392 1489176 : call phys_getopts(use_spcam_out = use_spcam)
393 :
394 1489176 : errstring = ''
395 1489176 : if( ( diffuse(fieldlist,'u') .or. diffuse(fieldlist,'v') ) .and. .not. diffuse(fieldlist,'s') ) then
396 0 : errstring = 'diffusion_solver.compute_vdiff: must diffuse s if diffusing u or v'
397 0 : return
398 : end if
399 :
400 : !--------------------------------------- !
401 : ! Computation of Molecular Diffusivities !
402 : !--------------------------------------- !
403 :
404 : ! Modification : Why 'kvq' is not changed by molecular diffusion ?
405 :
406 1489176 : if( do_molec_diff ) then
407 :
408 : if( (.not.present(vd_lu_qdecomp)) .or. (.not.present(kvt)) &
409 0 : .or. (.not. present(ubc_mmr)) .or. (.not. present(ubc_flux)) ) then
410 0 : errstring = 'compute_vdiff: do_molec_diff true but vd_lu_qdecomp or kvt missing'
411 0 : return
412 : endif
413 :
414 0 : p_molec = p%section([1, ncol], [1, nbot_molec])
415 0 : molec_boundary = BoundaryFixedLayer(p%del(:,nbot_molec+1))
416 :
417 : endif
418 :
419 : ! Boundary condition for a fixed concentration directly on a boundary
420 : ! interface (i.e. a boundary layer of size 0).
421 1489176 : interface_boundary = BoundaryFixedLayer(spread(0._r8, 1, ncol))
422 :
423 : ! Note that the *derivative* dp/dz is g*rho
424 2338872120 : dpidz_sq = gravit*rhoi(:ncol,:)
425 2338872120 : dpidz_sq = dpidz_sq * dpidz_sq
426 :
427 24865776 : rrho(:ncol) = rair * t(:ncol,pver) / p%mid(:,pver)
428 :
429 24865776 : tmpi2(:ncol,1) = ztodt * dpidz_sq(:,1) / ( p%mid(:,1) - p%ifc(:,1) )
430 2289140568 : tmpi2(:ncol,2:pver) = ztodt * dpidz_sq(:,2:pver) * p%rdst
431 :
432 : ! FIXME: The following four lines are kept in only to preserve answers;
433 : ! they really should be taken out completely.
434 1489176 : if (do_molec_diff) &
435 0 : tmpi2(:ncol,1) = ztodt * (gravit * rhoi(:ncol,1))**2 / ( pmid(:ncol,1) - p%ifc(:,1) )
436 24865776 : dpidz_sq(:,1) = gravit*(p%ifc(:,1) / (rairv(:ncol,1,lchnk)*t(:ncol,1)))
437 24865776 : dpidz_sq(:,1) = dpidz_sq(:,1)*dpidz_sq(:,1)
438 :
439 24865776 : tmp1(:ncol) = ztodt * gravit * p%rdel(:,pver)
440 :
441 : !---------------------------- !
442 : ! Diffuse Horizontal Momentum !
443 : !---------------------------- !
444 :
445 139982544 : do k = 1, pver
446 2314006344 : do i = 1, ncol
447 2312517168 : keg_in(i,k) = 0.5_r8 * ( u(i,k)*u(i,k) + v(i,k)*v(i,k) )
448 : end do
449 : end do
450 :
451 1489176 : if( diffuse(fieldlist,'u') .or. diffuse(fieldlist,'v') ) then
452 :
453 : ! Compute the vertical upward differences of the input u,v for KE dissipation
454 : ! at each interface.
455 : ! Velocity = 0 at surface, so difference at the bottom interface is -u,v(pver)
456 : ! These 'dinp_u, dinp_v' are computed using the non-diffused input wind.
457 :
458 24865776 : do i = 1, ncol
459 23376600 : dinp_u(i,1) = 0._r8
460 23376600 : dinp_v(i,1) = 0._r8
461 23376600 : dinp_u(i,pver+1) = -u(i,pver)
462 24865776 : dinp_v(i,pver+1) = -v(i,pver)
463 : end do
464 138493368 : do k = 2, pver
465 2289140568 : do i = 1, ncol
466 2150647200 : dinp_u(i,k) = u(i,k) - u(i,k-1)
467 2287651392 : dinp_v(i,k) = v(i,k) - v(i,k-1)
468 : end do
469 : end do
470 :
471 : ! -------------------------------------------------------------- !
472 : ! Do 'Implicit Surface Stress' treatment for numerical stability !
473 : ! in the lowest model layer. !
474 : ! -------------------------------------------------------------- !
475 :
476 1489176 : if( do_iss ) then
477 :
478 : ! Compute surface drag coefficient for implicit diffusion
479 : ! including turbulent mountain stress.
480 :
481 24865776 : do i = 1, ncol
482 23376600 : ws(i) = max( sqrt( u(i,pver)**2._r8 + v(i,pver)**2._r8 ), wsmin )
483 23376600 : tau(i) = sqrt( taux(i)**2._r8 + tauy(i)**2._r8 )
484 24865776 : ksrfturb(i) = max( tau(i) / ws(i), ksrfmin )
485 : end do
486 24865776 : ksrf(:ncol) = ksrfturb(:ncol) + ksrftms(:ncol) ! Do all surface stress ( normal + tms ) implicitly
487 :
488 : ! Vertical integration of input momentum.
489 : ! This is total horizontal momentum per unit area [ kg*m/s/m2 ] in each column.
490 : ! Note (u,v) are the raw input to the PBL scheme, not the
491 : ! provisionally-marched ones within the iteration loop of the PBL scheme.
492 :
493 24865776 : do i = 1, ncol
494 23376600 : usum_in(i) = 0._r8
495 23376600 : vsum_in(i) = 0._r8
496 2198889576 : do k = 1, pver
497 2174023800 : usum_in(i) = usum_in(i) + (1._r8/gravit)*u(i,k)*p%del(i,k)
498 2197400400 : vsum_in(i) = vsum_in(i) + (1._r8/gravit)*v(i,k)*p%del(i,k)
499 : end do
500 : end do
501 :
502 : ! Add residual stress of previous time step explicitly into the lowest
503 : ! model layer with a relaxation time scale of 'timeres'.
504 :
505 1489176 : if (am_correction) then
506 : ! preserve time-mean torque
507 : ramda = 1._r8
508 : else
509 1489176 : ramda = ztodt / timeres
510 : endif
511 :
512 24865776 : u(:ncol,pver) = u(:ncol,pver) + tmp1(:ncol)*tauresx(:ncol)*ramda
513 24865776 : v(:ncol,pver) = v(:ncol,pver) + tmp1(:ncol)*tauresy(:ncol)*ramda
514 :
515 : ! Vertical integration of momentum after adding explicit residual stress
516 : ! into the lowest model layer.
517 :
518 24865776 : do i = 1, ncol
519 23376600 : usum_mid(i) = 0._r8
520 23376600 : vsum_mid(i) = 0._r8
521 2198889576 : do k = 1, pver
522 2174023800 : usum_mid(i) = usum_mid(i) + (1._r8/gravit)*u(i,k)*p%del(i,k)
523 2197400400 : vsum_mid(i) = vsum_mid(i) + (1._r8/gravit)*v(i,k)*p%del(i,k)
524 : end do
525 : end do
526 :
527 : else
528 :
529 : ! In this case, do 'turbulent mountain stress' implicitly,
530 : ! but do 'normal turbulent stress' explicitly.
531 : ! In this case, there is no 'residual stress' as long as 'tms' is
532 : ! treated in a fully implicit way, which is true.
533 :
534 : ! 1. Do 'tms' implicitly
535 :
536 0 : ksrf(:ncol) = ksrftms(:ncol)
537 :
538 : ! 2. Do 'normal stress' explicitly
539 :
540 0 : u(:ncol,pver) = u(:ncol,pver) + tmp1(:ncol)*taux(:ncol)
541 0 : v(:ncol,pver) = v(:ncol,pver) + tmp1(:ncol)*tauy(:ncol)
542 :
543 : end if ! End of 'do iss' ( implicit surface stress )
544 :
545 : ! --------------------------------------------------------------------------------------- !
546 : ! Diffuse horizontal momentum implicitly using tri-diagnonal matrix. !
547 : ! The 'u,v' are input-output: the output 'u,v' are implicitly diffused winds. !
548 : ! For implicit 'normal' stress : ksrf = ksrftms + ksrfturb, !
549 : ! u(pver) : explicitly include 'residual normal' stress !
550 : ! For explicit 'normal' stress : ksrf = ksrftms !
551 : ! u(pver) : explicitly include 'normal' stress !
552 : ! Note that in all the two cases above, 'tms' is fully implicitly treated. !
553 : ! --------------------------------------------------------------------------------------- !
554 :
555 : ! In most layers, no damping at all.
556 2314006344 : tau_damp_rate = 0._r8
557 :
558 : ! Physical interpretation:
559 : ! ksrf is stress per unit wind speed.
560 : ! p%del / gravit is approximately the mass in the layer per unit of
561 : ! surface area.
562 : ! Therefore, gravit*ksrf/p%del is the acceleration of wind per unit
563 : ! wind speed, i.e. the rate at which wind is exponentially damped by
564 : ! surface stress.
565 :
566 : ! Beljaars et al SGO scheme incorporated here. It
567 : ! appears as a "3D" tau_damp_rate specification.
568 :
569 24865776 : tau_damp_rate(:,pver) = -gravit*ksrf(:ncol)*p%rdel(:,pver)
570 139982544 : do k=1,pver
571 2314006344 : tau_damp_rate(:,k) = tau_damp_rate(:,k) + dragblj(:ncol,k)
572 : end do
573 :
574 : v(:ncol,:) = fin_vol_solve(ztodt, p, v(:ncol,:), ncol, pver, &
575 : coef_q=tau_damp_rate, &
576 5193553248 : coef_q_diff=kvm(:ncol,:)*dpidz_sq)
577 :
578 : u(:ncol,:) = fin_vol_solve(ztodt, p, u(:ncol,:), ncol, pver, &
579 : coef_q=tau_damp_rate, &
580 5193553248 : coef_q_diff=kvm(:ncol,:)*dpidz_sq)
581 :
582 :
583 :
584 : ! ---------------------------------------------------------------------- !
585 : ! Calculate 'total' ( tautotx ) and 'tms' ( tautmsx ) stresses that !
586 : ! have been actually added into the atmosphere at the current time step. !
587 : ! Also, update residual stress, if required. !
588 : ! ---------------------------------------------------------------------- !
589 :
590 24865776 : do i = 1, ncol
591 :
592 : ! Compute the implicit 'tms' using the updated winds.
593 : ! Below 'tautmsx(i),tautmsy(i)' are pure implicit mountain stresses
594 : ! that has been actually added into the atmosphere both for explicit
595 : ! and implicit approach.
596 :
597 23376600 : tautmsx(i) = -ksrftms(i)*u(i,pver)
598 23376600 : tautmsy(i) = -ksrftms(i)*v(i,pver)
599 :
600 : ! We want to add vertically-integrated Beljaars drag to residual stress.
601 : ! So this has to be calculated locally.
602 : ! We may want to rethink the residual drag calculation performed here on. (jtb)
603 23376600 : taubljx(i) = 0._r8
604 23376600 : taubljy(i) = 0._r8
605 2197400400 : do k = 1, pver
606 2174023800 : taubljx(i) = taubljx(i) + (1._r8/gravit)*dragblj(i,k)*u(i,k)*p%del(i,k)
607 2197400400 : taubljy(i) = taubljy(i) + (1._r8/gravit)*dragblj(i,k)*v(i,k)*p%del(i,k)
608 : end do
609 :
610 24865776 : if( do_iss ) then
611 :
612 : ! Compute vertical integration of final horizontal momentum
613 :
614 23376600 : usum_out(i) = 0._r8
615 23376600 : vsum_out(i) = 0._r8
616 2197400400 : do k = 1, pver
617 2174023800 : usum_out(i) = usum_out(i) + (1._r8/gravit)*u(i,k)*p%del(i,k)
618 2197400400 : vsum_out(i) = vsum_out(i) + (1._r8/gravit)*v(i,k)*p%del(i,k)
619 : end do
620 :
621 : ! Compute net stress added into the atmosphere at the current time step.
622 : ! Note that the difference between 'usum_in' and 'usum_out' are induced
623 : ! by 'explicit residual stress + implicit total stress' for implicit case, while
624 : ! by 'explicit normal stress + implicit tms stress' for explicit case.
625 : ! Here, 'tautotx(i)' is net stress added into the air at the current time step.
626 :
627 23376600 : tauimpx(i) = ( usum_out(i) - usum_in(i) ) / ztodt
628 23376600 : tauimpy(i) = ( vsum_out(i) - vsum_in(i) ) / ztodt
629 :
630 23376600 : tautotx(i) = tauimpx(i)
631 23376600 : tautoty(i) = tauimpy(i)
632 :
633 : ! Compute residual stress and update if required.
634 : ! Note that the total stress we should have added at the current step is
635 : ! the sum of 'taux(i) - ksrftms(i)*u(i,pver) + tauresx(i)'.
636 :
637 23376600 : if( itaures .eq. 1 ) then
638 23376600 : tauresx(i) = taux(i) + tautmsx(i) + taubljx(i) + tauresx(i)- tauimpx(i)
639 23376600 : tauresy(i) = tauy(i) + tautmsy(i) + taubljy(i) + tauresy(i)- tauimpy(i)
640 : endif
641 :
642 : else
643 :
644 0 : tautotx(i) = tautmsx(i) + taux(i)
645 0 : tautoty(i) = tautmsy(i) + tauy(i)
646 0 : tauresx(i) = 0._r8
647 0 : tauresy(i) = 0._r8
648 :
649 : end if ! End of 'do_iss' if
650 :
651 : end do ! End of 'do i = 1, ncol' loop
652 :
653 : ! ------------------------------------ !
654 : ! Calculate kinetic energy dissipation !
655 : ! ------------------------------------ !
656 :
657 : ! Modification : In future, this should be set exactly same as
658 : ! the ones in the convection schemes
659 :
660 : ! 1. Compute dissipation term at interfaces
661 : ! Note that 'u,v' are already diffused wind, and 'tautotx,tautoty' are
662 : ! implicit stress that has been actually added. On the other hand,
663 : ! 'dinp_u, dinp_v' were computed using non-diffused input wind.
664 :
665 : ! Modification : I should check whether non-consistency between 'u' and 'dinp_u'
666 : ! is correctly intended approach. I think so.
667 :
668 24865776 : k = pver + 1
669 24865776 : do i = 1, ncol
670 23376600 : tmpi1(i,1) = 0._r8
671 : tmpi1(i,k) = 0.5_r8 * ztodt * gravit * &
672 24865776 : ( (-u(i,k-1) + dinp_u(i,k))*tautotx(i) + (-v(i,k-1) + dinp_v(i,k))*tautoty(i) )
673 : end do
674 :
675 138493368 : do k = 2, pver
676 2289140568 : do i = 1, ncol
677 2150647200 : dout_u = u(i,k) - u(i,k-1)
678 2150647200 : dout_v = v(i,k) - v(i,k-1)
679 2150647200 : tmpi1(i,k) = 0.25_r8 * tmpi2(i,k) * kvm(i,k) * &
680 2287651392 : ( dout_u**2 + dout_v**2 + dout_u*dinp_u(i,k) + dout_v*dinp_v(i,k) )
681 : end do
682 : end do
683 :
684 1489176 : if (do_beljaars) then
685 :
686 : ! 2. Add Kinetic Energy change across dissipation to Static Energy
687 139982544 : do k = 1, pver
688 2314006344 : do i = 1, ncol
689 2312517168 : keg_out(i,k) = 0.5_r8 * ( u(i,k)*u(i,k) + v(i,k)*v(i,k) )
690 : end do
691 : end do
692 :
693 139982544 : do k = 1, pver
694 2314006344 : do i = 1, ncol
695 2174023800 : dtk(i,k) = keg_in(i,k) - keg_out(i,k)
696 2312517168 : dse(i,k) = dse(i,k) + dtk(i,k) ! + dkeblj(i,k)
697 : end do
698 : end do
699 :
700 : else
701 :
702 : ! 2. Compute dissipation term at midpoints, add to dry static energy
703 0 : do k = 1, pver
704 0 : do i = 1, ncol
705 0 : dtk(i,k) = ( tmpi1(i,k+1) + tmpi1(i,k) ) * p%rdel(i,k)
706 0 : dse(i,k) = dse(i,k) + dtk(i,k)
707 : end do
708 : end do
709 :
710 : end if
711 :
712 : end if ! End of diffuse horizontal momentum, diffuse(fieldlist,'u') routine
713 :
714 : !-------------------------- !
715 : ! Diffuse Dry Static Energy !
716 : !-------------------------- !
717 :
718 : ! Modification : In future, we should diffuse the fully conservative
719 : ! moist static energy,not the dry static energy.
720 :
721 1489176 : if( diffuse(fieldlist,'s') ) then
722 1489176 : if (.not. use_spcam) then
723 :
724 : ! Add counter-gradient to input static energy profiles
725 :
726 139982544 : do k = 1, pver
727 276986736 : dse(:ncol,k) = dse(:ncol,k) + ztodt * p%rdel(:,k) * gravit * &
728 138493368 : ( rhoi(:ncol,k+1) * kvh(:ncol,k+1) * cgh(:ncol,k+1) &
729 2729486448 : - rhoi(:ncol,k ) * kvh(:ncol,k ) * cgh(:ncol,k ) )
730 : end do
731 : endif
732 : ! Add the explicit surface fluxes to the lowest layer
733 24865776 : dse(:ncol,pver) = dse(:ncol,pver) + tmp1(:ncol) * shflx(:ncol)
734 :
735 : ! Diffuse dry static energy
736 :
737 : !---------------------------------------------------
738 : ! Solve for temperature using thermal conductivity
739 : !---------------------------------------------------
740 1489176 : if ( use_temperature_molec_diff ) then
741 : !----------------------------------------------------------------------------------------------------
742 : ! In Extended WACCM, kvt is calculated rather kvh. This is because molecular diffusion operates on
743 : ! temperature, while eddy diffusion operates on dse. Also, pass in constituent dependent "constants"
744 : !----------------------------------------------------------------------------------------------------
745 :
746 : ! Boundary layer thickness of "0._r8" signifies that the boundary
747 : ! condition is defined directly on the top interface.
748 :
749 0 : if (.not. use_spcam) then
750 : dse(:ncol,:) = fin_vol_solve(ztodt, p, dse(:ncol,:), ncol, pver, &
751 : coef_q_diff=kvh(:ncol,:)*dpidz_sq, &
752 : upper_bndry=interface_boundary, &
753 0 : l_cond=BoundaryData(dse_top(:ncol)))
754 : endif
755 :
756 : ! Calculate flux at top interface
757 :
758 : ! Modification : Why molecular diffusion does not work for dry static energy in all layers ?
759 :
760 : topflx(:ncol) = - kvh(:ncol,1) * tmpi2(:ncol,1) / (ztodt*gravit) * &
761 0 : ( dse(:ncol,1) - dse_top(:ncol) )
762 :
763 0 : ttemp0 = t(:ncol,:)
764 0 : ttemp = ttemp0
765 :
766 : ! upper boundary is zero flux for extended model
767 0 : if (.not. use_spcam) then
768 : ttemp = fin_vol_solve(ztodt, p, ttemp, ncol, pver, &
769 0 : coef_q_diff=kvt(:ncol,:)*dpidz_sq, &
770 0 : coef_q_weight=cpairv(:ncol,:))
771 : end if
772 :
773 :
774 : !-------------------------------------
775 : ! Update dry static energy
776 : !-------------------------------------
777 0 : do k = 1,pver
778 0 : dse(:ncol,k) = dse(:ncol,k) + &
779 0 : cpairv(:ncol,k)*(ttemp(:,k) - ttemp0(:,k))
780 : enddo
781 :
782 : else
783 :
784 1489176 : if (do_molec_diff) then
785 0 : kv_total(:ncol,:) = kvh(:ncol,:) + kvt(:ncol,:)/cpair
786 : else
787 2338872120 : kv_total(:ncol,:) = kvh(:ncol,:)
788 : end if
789 :
790 : ! Boundary layer thickness of "0._r8" signifies that the boundary
791 : ! condition is defined directly on the top interface.
792 1489176 : if (.not. use_spcam) then
793 : dse(:ncol,:) = fin_vol_solve(ztodt, p, dse(:ncol,:), ncol, pver, &
794 : coef_q_diff=kv_total(:ncol,:)*dpidz_sq, &
795 : upper_bndry=interface_boundary, &
796 5193553248 : l_cond=BoundaryData(dse_top(:ncol)))
797 : end if
798 :
799 : ! Calculate flux at top interface
800 :
801 : ! Modification : Why molecular diffusion does not work for dry static energy in all layers ?
802 :
803 1489176 : if( do_molec_diff ) then
804 : topflx(:ncol) = - kv_total(:ncol,1) * tmpi2(:ncol,1) / (ztodt*gravit) * &
805 0 : ( dse(:ncol,1) - dse_top(:ncol) )
806 : else
807 24865776 : topflx(:ncol) = 0._r8
808 : end if
809 :
810 : endif
811 :
812 : endif
813 :
814 : !---------------------------- !
815 : ! Diffuse Water Vapor Tracers !
816 : !---------------------------- !
817 :
818 : ! Modification : For aerosols, I need to use separate treatment
819 : ! for aerosol mass and aerosol number.
820 :
821 : ! Loop through constituents
822 :
823 : no_molec_decomp = fin_vol_lu_decomp(ztodt, p, &
824 2338872120 : coef_q_diff=kvq(:ncol,:)*dpidz_sq)
825 :
826 62545392 : do m = 1, ncnst
827 :
828 62545392 : if( diffuse(fieldlist,'q',m) ) then
829 31272696 : if (.not. use_spcam) then
830 :
831 : ! Add the nonlocal transport terms to constituents in the PBL.
832 : ! Check for neg q's in each constituent and put the original vertical
833 : ! profile back if a neg value is found. A neg value implies that the
834 : ! quasi-equilibrium conditions assumed for the countergradient term are
835 : ! strongly violated.
836 :
837 48594133224 : qtm(:ncol,:pver) = q(:ncol,:pver,m)
838 :
839 2939633424 : do k = 1, pver
840 2908360728 : q(:ncol,k,m) = q(:ncol,k,m) + &
841 2908360728 : ztodt * p%rdel(:,k) * gravit * ( cflx(:ncol,m) * rrho(:ncol) ) * &
842 2908360728 : ( rhoi(:ncol,k+1) * kvh(:ncol,k+1) * cgs(:ncol,k+1) &
843 57319215408 : - rhoi(:ncol,k ) * kvh(:ncol,k ) * cgs(:ncol,k ) )
844 : end do
845 48594133224 : lqtst(:ncol) = all(q(:ncol,1:pver,m) >= qmincg(m), 2)
846 2939633424 : do k = 1, pver
847 48594133224 : q(:ncol,k,m) = merge( q(:ncol,k,m), qtm(:ncol,k), lqtst(:ncol) )
848 : end do
849 : endif
850 :
851 : ! Add the explicit surface fluxes to the lowest layer
852 :
853 522181296 : q(:ncol,pver,m) = q(:ncol,pver,m) + tmp1(:ncol) * cflx(:ncol,m)
854 :
855 : ! Diffuse constituents.
856 :
857 : ! This is for solving molecular diffusion of minor species, thus, for WACCM-X, bypass O and O2 (major species)
858 : ! Major species diffusion is calculated separately. -Hanli Liu
859 :
860 31272696 : if( do_molec_diff .and. diffuse(fieldlistm,'q',m)) then
861 :
862 0 : decomp = vd_lu_qdecomp( pcols , pver , ncol , cnst_fixed_ubc(m), cnst_mw(m), &
863 : kvq , kq_scal, mw_fac(:,:,m) ,dpidz_sq , p_molec, &
864 : interface_boundary, molec_boundary, &
865 : tint , ztodt , nbot_molec , &
866 0 : lchnk , t , m , no_molec_decomp)
867 :
868 : ! This to calculate the upper boundary flux of H. -Hanli Liu
869 0 : if ((cnst_fixed_ubflx(m))) then
870 :
871 : ! ubc_flux is a flux of mass density through space, i.e.:
872 : ! ubc_flux = rho_i * dz/dt = q_i * rho * dz/dt
873 : ! For flux of mmr through pressure level, multiply by g:
874 : ! q_i * rho * gravit * dz/dt = q_i * dp/dt
875 :
876 : call decomp%left_div(q(:ncol,:,m), &
877 : l_cond=BoundaryFlux( &
878 0 : -gravit*ubc_flux(:ncol,m), ztodt, &
879 0 : p%del(:,1)))
880 :
881 : else
882 : call decomp%left_div(q(:ncol,:,m), &
883 0 : l_cond=BoundaryData(ubc_mmr(:ncol,m)))
884 : end if
885 :
886 0 : call decomp%finalize()
887 :
888 : else
889 :
890 31272696 : if (present(cnst_fixed_ubc)) then
891 : ! explicitly set mmr in top layer for cases where molecular diffusion is not active
892 31272696 : if (cnst_fixed_ubc(m)) then
893 0 : q(:ncol,1,m) = ubc_mmr(:ncol,m)
894 : endif
895 : end if
896 :
897 31272696 : if (.not. use_spcam) then
898 31272696 : call no_molec_decomp%left_div(q(:ncol,:,m))
899 : end if
900 :
901 : end if
902 :
903 : end if
904 : end do
905 :
906 1489176 : call no_molec_decomp%finalize()
907 :
908 13402584 : end subroutine compute_vdiff
909 :
910 : ! =============================================================================== !
911 : ! !
912 : ! =============================================================================== !
913 :
914 69120 : character(128) function vdiff_select( fieldlist, name, qindex )
915 : ! --------------------------------------------------------------------- !
916 : ! This function sets the field with incoming name as one to be diffused !
917 : ! --------------------------------------------------------------------- !
918 : type(vdiff_selector), intent(inout) :: fieldlist
919 : character(*), intent(in) :: name
920 : integer, intent(in), optional :: qindex
921 :
922 69120 : vdiff_select = ''
923 1536 : select case (name)
924 : case ('u','U')
925 1536 : fieldlist%fields(1) = .true.
926 : case ('v','V')
927 1536 : fieldlist%fields(2) = .true.
928 : case ('s','S')
929 1536 : fieldlist%fields(3) = .true.
930 : case ('q','Q')
931 64512 : if( present(qindex) ) then
932 64512 : fieldlist%fields(3 + qindex) = .true.
933 : else
934 0 : fieldlist%fields(4) = .true.
935 : endif
936 : case default
937 69120 : write(vdiff_select,*) 'Bad argument to vdiff_index: ', name
938 : end select
939 69120 : return
940 :
941 : end function vdiff_select
942 :
943 0 : type(vdiff_selector) function not(a)
944 : ! ------------------------------------------------------------- !
945 : ! This function extends .not. to operate on type vdiff_selector !
946 : ! ------------------------------------------------------------- !
947 : type(vdiff_selector), intent(in) :: a
948 0 : allocate(not%fields(size(a%fields)))
949 0 : not%fields = .not. a%fields
950 0 : end function not
951 :
952 2978352 : logical function my_any(a)
953 : ! -------------------------------------------------- !
954 : ! This function extends the intrinsic function 'any' !
955 : ! to operate on type vdiff_selector !
956 : ! -------------------------------------------------- !
957 : type(vdiff_selector), intent(in) :: a
958 68502096 : my_any = any(a%fields)
959 2978352 : end function my_any
960 :
961 101263968 : logical function diffuse(fieldlist,name,qindex)
962 : ! ---------------------------------------------------------------------------- !
963 : ! This function reports whether the field with incoming name is to be diffused !
964 : ! ---------------------------------------------------------------------------- !
965 : type(vdiff_selector), intent(in) :: fieldlist
966 : character(*), intent(in) :: name
967 : integer, intent(in), optional :: qindex
968 :
969 2978352 : select case (name)
970 : case ('u','U')
971 2978352 : diffuse = fieldlist%fields(1)
972 : case ('v','V')
973 2978352 : diffuse = fieldlist%fields(2)
974 : case ('s','S')
975 2978352 : diffuse = fieldlist%fields(3)
976 : case ('q','Q')
977 92328912 : if( present(qindex) ) then
978 92328912 : diffuse = fieldlist%fields(3 + qindex)
979 : else
980 0 : diffuse = fieldlist%fields(4)
981 : endif
982 : case default
983 101263968 : diffuse = .false.
984 : end select
985 : return
986 : end function diffuse
987 :
988 0 : end module diffusion_solver
|
