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