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 1022976 : new_fieldlist_vdiff%fields = .false.
126 :
127 4608 : end function new_fieldlist_vdiff
128 :
129 : ! =============================================================================== !
130 : ! !
131 : ! =============================================================================== !
132 :
133 72960 : subroutine compute_vdiff( lchnk , &
134 72960 : pcols , pver , ncnst , ncol , tint , &
135 145920 : p , t , rhoi , ztodt , taux , &
136 72960 : tauy , shflx , cflx , &
137 72960 : kvh , kvm , kvq , cgs , cgh , &
138 72960 : zi , ksrftms , dragblj , &
139 72960 : qmincg , fieldlist , fieldlistm , &
140 72960 : u , v , q , dse , &
141 0 : tautmsx , tautmsy , dtk , topflx , errstring , &
142 72960 : tauresx , tauresy , itaures , cpairv , dse_top, &
143 : do_molec_diff , use_temperature_molec_diff, vd_lu_qdecomp, &
144 291840 : ubc_mmr, ubc_flux, kvt, pmid, &
145 218880 : cnst_mw, cnst_fixed_ubc, cnst_fixed_ubflx, nbot_molec, &
146 72960 : 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 145920 : logical :: lqtst(pcols) ! Adjust vertical profiles
312 :
313 : ! LU decomposition information.
314 72960 : type(TriDiagDecomp) :: decomp
315 72960 : type(TriDiagDecomp) :: no_molec_decomp
316 :
317 : ! Square of derivative of pressure with height (on interfaces).
318 145920 : real(r8) :: dpidz_sq(ncol,pver+1)
319 :
320 : ! Pressure coordinates over the molecular diffusion range only.
321 72960 : type(Coords1D) :: p_molec
322 :
323 : ! Boundary layer objects
324 72960 : type(BoundaryType) :: interface_boundary
325 72960 : type(BoundaryType) :: molec_boundary
326 :
327 145920 : real(r8) :: tmp1(pcols) ! Temporary storage
328 145920 : real(r8) :: tmpi1(pcols,pver+1) ! Interface KE dissipation
329 145920 : real(r8) :: tmpi2(pcols,pver+1) ! dt*(g*rho)**2/dp at interfaces
330 145920 : real(r8) :: keg_in(pcols,pver) ! KE on entry to subroutine
331 145920 : real(r8) :: keg_out(pcols,pver) ! KE after U and V dissipation/diffusion
332 145920 : real(r8) :: rrho(pcols) ! 1./bottom level density
333 :
334 145920 : real(r8) :: tautotx(pcols) ! Total surface stress ( zonal )
335 145920 : real(r8) :: tautoty(pcols) ! Total surface stress ( meridional )
336 :
337 145920 : real(r8) :: dinp_u(pcols,pver+1) ! Vertical difference at interfaces, input u
338 145920 : 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 145920 : real(r8) :: qtm(pcols,pver) ! Temporary copy of q
343 :
344 145920 : real(r8) :: ws(pcols) ! Lowest-level wind speed [ m/s ]
345 145920 : real(r8) :: tau(pcols) ! Turbulent surface stress ( not including mountain stress )
346 145920 : 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 145920 : real(r8) :: ksrf(pcols) ! Surface drag coefficient of 'normal' stress +
349 : ! Surface drag coefficient of 'tms' stress. > 0. [ kg/s/m2 ]
350 145920 : 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 145920 : 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 145920 : real(r8) :: usum_mid(pcols) ! Vertical integral of u-momentum after adding explicit residual stress
355 145920 : real(r8) :: vsum_mid(pcols) ! Vertical integral of v-momentum after adding explicit residual stress
356 145920 : real(r8) :: usum_out(pcols) ! Vertical integral of u-momentum after doing implicit diffusion
357 145920 : real(r8) :: vsum_out(pcols) ! Vertical integral of v-momentum after doing implicit diffusion
358 145920 : real(r8) :: tauimpx(pcols) ! Actual net stress added at the current step other than mountain stress
359 145920 : 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 145920 : real(r8) :: taubljx(pcols) ! recomputed explicit/residual beljaars stress
363 145920 : real(r8) :: taubljy(pcols) ! recomputed explicit/residual beljaars stress
364 :
365 : ! Rate at which external (surface) stress damps wind speeds (1/s).
366 145920 : real(r8) :: tau_damp_rate(ncol, pver)
367 :
368 : ! Combined molecular and eddy diffusion.
369 145920 : real(r8) :: kv_total(pcols,pver+1)
370 :
371 : !--------------------------------
372 : ! Variables needed for WACCM-X
373 : !--------------------------------
374 145920 : real(r8) :: ttemp(ncol,pver) ! temporary temperature array
375 72960 : real(r8) :: ttemp0(ncol,pver) ! temporary temperature array
376 :
377 : ! ------------------------------------------------ !
378 : ! Parameters for implicit surface stress treatment !
379 : ! ------------------------------------------------ !
380 :
381 : real(r8), parameter :: wsmin = 1._r8 ! Minimum sfc wind speed for estimating frictional
382 : ! transfer velocity ksrf. [ m/s ]
383 : real(r8), parameter :: ksrfmin = 1.e-4_r8 ! Minimum surface drag coefficient [ kg/s/m^2 ]
384 : real(r8), parameter :: timeres = 7200._r8 ! Relaxation time scale of residual stress ( >= dt ) [ s ]
385 :
386 : ! ----------------------- !
387 : ! Main Computation Begins !
388 : ! ----------------------- !
389 :
390 72960 : errstring = ''
391 72960 : if( ( diffuse(fieldlist,'u') .or. diffuse(fieldlist,'v') ) .and. .not. diffuse(fieldlist,'s') ) then
392 0 : errstring = 'diffusion_solver.compute_vdiff: must diffuse s if diffusing u or v'
393 0 : return
394 : end if
395 :
396 : !--------------------------------------- !
397 : ! Computation of Molecular Diffusivities !
398 : !--------------------------------------- !
399 :
400 : ! Modification : Why 'kvq' is not changed by molecular diffusion ?
401 :
402 72960 : if( do_molec_diff ) then
403 :
404 : if( (.not.present(vd_lu_qdecomp)) .or. (.not.present(kvt)) &
405 72960 : .or. (.not. present(ubc_mmr)) .or. (.not. present(ubc_flux)) ) then
406 0 : errstring = 'compute_vdiff: do_molec_diff true but vd_lu_qdecomp or kvt missing'
407 0 : return
408 : endif
409 :
410 364800 : p_molec = p%section([1, ncol], [1, nbot_molec])
411 72960 : molec_boundary = BoundaryFixedLayer(p%del(:,nbot_molec+1))
412 :
413 : endif
414 :
415 : ! Boundary condition for a fixed concentration directly on a boundary
416 : ! interface (i.e. a boundary layer of size 0).
417 72960 : interface_boundary = BoundaryFixedLayer(spread(0._r8, 1, ncol))
418 :
419 : ! Note that the *derivative* dp/dz is g*rho
420 79847424 : dpidz_sq = gravit*rhoi(:ncol,:)
421 79847424 : dpidz_sq = dpidz_sq * dpidz_sq
422 :
423 1123584 : rrho(:ncol) = rair * t(:ncol,pver) / p%mid(:,pver)
424 :
425 1123584 : tmpi2(:ncol,1) = ztodt * dpidz_sq(:,1) / ( p%mid(:,1) - p%ifc(:,1) )
426 77600256 : tmpi2(:ncol,2:pver) = ztodt * dpidz_sq(:,2:pver) * p%rdst
427 :
428 : ! FIXME: The following four lines are kept in only to preserve answers;
429 : ! they really should be taken out completely.
430 72960 : if (do_molec_diff) &
431 1123584 : tmpi2(:ncol,1) = ztodt * (gravit * rhoi(:ncol,1))**2 / ( pmid(:ncol,1) - p%ifc(:,1) )
432 1123584 : dpidz_sq(:,1) = gravit*(p%ifc(:,1) / (rairv(:ncol,1,lchnk)*t(:ncol,1)))
433 1123584 : dpidz_sq(:,1) = dpidz_sq(:,1)*dpidz_sq(:,1)
434 :
435 1123584 : tmp1(:ncol) = ztodt * gravit * p%rdel(:,pver)
436 :
437 : !---------------------------- !
438 : ! Diffuse Horizontal Momentum !
439 : !---------------------------- !
440 :
441 5180160 : do k = 1, pver
442 78723840 : do i = 1, ncol
443 78650880 : keg_in(i,k) = 0.5_r8 * ( u(i,k)*u(i,k) + v(i,k)*v(i,k) )
444 : end do
445 : end do
446 :
447 72960 : if( diffuse(fieldlist,'u') .or. diffuse(fieldlist,'v') ) then
448 :
449 : ! Compute the vertical upward differences of the input u,v for KE dissipation
450 : ! at each interface.
451 : ! Velocity = 0 at surface, so difference at the bottom interface is -u,v(pver)
452 : ! These 'dinp_u, dinp_v' are computed using the non-diffused input wind.
453 :
454 1123584 : do i = 1, ncol
455 1050624 : dinp_u(i,1) = 0._r8
456 1050624 : dinp_v(i,1) = 0._r8
457 1050624 : dinp_u(i,pver+1) = -u(i,pver)
458 1123584 : dinp_v(i,pver+1) = -v(i,pver)
459 : end do
460 5107200 : do k = 2, pver
461 77600256 : do i = 1, ncol
462 72493056 : dinp_u(i,k) = u(i,k) - u(i,k-1)
463 77527296 : dinp_v(i,k) = v(i,k) - v(i,k-1)
464 : end do
465 : end do
466 :
467 : ! -------------------------------------------------------------- !
468 : ! Do 'Implicit Surface Stress' treatment for numerical stability !
469 : ! in the lowest model layer. !
470 : ! -------------------------------------------------------------- !
471 :
472 72960 : if( do_iss ) then
473 :
474 : ! Compute surface drag coefficient for implicit diffusion
475 : ! including turbulent mountain stress.
476 :
477 1123584 : do i = 1, ncol
478 1050624 : ws(i) = max( sqrt( u(i,pver)**2._r8 + v(i,pver)**2._r8 ), wsmin )
479 1050624 : tau(i) = sqrt( taux(i)**2._r8 + tauy(i)**2._r8 )
480 1123584 : ksrfturb(i) = max( tau(i) / ws(i), ksrfmin )
481 : end do
482 1123584 : ksrf(:ncol) = ksrfturb(:ncol) + ksrftms(:ncol) ! Do all surface stress ( normal + tms ) implicitly
483 :
484 : ! Vertical integration of input momentum.
485 : ! This is total horizontal momentum per unit area [ kg*m/s/m2 ] in each column.
486 : ! Note (u,v) are the raw input to the PBL scheme, not the
487 : ! provisionally-marched ones within the iteration loop of the PBL scheme.
488 :
489 1123584 : do i = 1, ncol
490 1050624 : usum_in(i) = 0._r8
491 1050624 : vsum_in(i) = 0._r8
492 74667264 : do k = 1, pver
493 73543680 : usum_in(i) = usum_in(i) + (1._r8/gravit)*u(i,k)*p%del(i,k)
494 74594304 : vsum_in(i) = vsum_in(i) + (1._r8/gravit)*v(i,k)*p%del(i,k)
495 : end do
496 : end do
497 :
498 : ! Add residual stress of previous time step explicitly into the lowest
499 : ! model layer with a relaxation time scale of 'timeres'.
500 :
501 72960 : if (am_correction) then
502 : ! preserve time-mean torque
503 : ramda = 1._r8
504 : else
505 72960 : ramda = ztodt / timeres
506 : endif
507 :
508 1123584 : u(:ncol,pver) = u(:ncol,pver) + tmp1(:ncol)*tauresx(:ncol)*ramda
509 1123584 : v(:ncol,pver) = v(:ncol,pver) + tmp1(:ncol)*tauresy(:ncol)*ramda
510 :
511 : ! Vertical integration of momentum after adding explicit residual stress
512 : ! into the lowest model layer.
513 :
514 1123584 : do i = 1, ncol
515 1050624 : usum_mid(i) = 0._r8
516 1050624 : vsum_mid(i) = 0._r8
517 74667264 : do k = 1, pver
518 73543680 : usum_mid(i) = usum_mid(i) + (1._r8/gravit)*u(i,k)*p%del(i,k)
519 74594304 : vsum_mid(i) = vsum_mid(i) + (1._r8/gravit)*v(i,k)*p%del(i,k)
520 : end do
521 : end do
522 :
523 : else
524 :
525 : ! In this case, do 'turbulent mountain stress' implicitly,
526 : ! but do 'normal turbulent stress' explicitly.
527 : ! In this case, there is no 'residual stress' as long as 'tms' is
528 : ! treated in a fully implicit way, which is true.
529 :
530 : ! 1. Do 'tms' implicitly
531 :
532 0 : ksrf(:ncol) = ksrftms(:ncol)
533 :
534 : ! 2. Do 'normal stress' explicitly
535 :
536 0 : u(:ncol,pver) = u(:ncol,pver) + tmp1(:ncol)*taux(:ncol)
537 0 : v(:ncol,pver) = v(:ncol,pver) + tmp1(:ncol)*tauy(:ncol)
538 :
539 : end if ! End of 'do iss' ( implicit surface stress )
540 :
541 : ! --------------------------------------------------------------------------------------- !
542 : ! Diffuse horizontal momentum implicitly using tri-diagnonal matrix. !
543 : ! The 'u,v' are input-output: the output 'u,v' are implicitly diffused winds. !
544 : ! For implicit 'normal' stress : ksrf = ksrftms + ksrfturb, !
545 : ! u(pver) : explicitly include 'residual normal' stress !
546 : ! For explicit 'normal' stress : ksrf = ksrftms !
547 : ! u(pver) : explicitly include 'normal' stress !
548 : ! Note that in all the two cases above, 'tms' is fully implicitly treated. !
549 : ! --------------------------------------------------------------------------------------- !
550 :
551 : ! In most layers, no damping at all.
552 78723840 : tau_damp_rate = 0._r8
553 :
554 : ! Physical interpretation:
555 : ! ksrf is stress per unit wind speed.
556 : ! p%del / gravit is approximately the mass in the layer per unit of
557 : ! surface area.
558 : ! Therefore, gravit*ksrf/p%del is the acceleration of wind per unit
559 : ! wind speed, i.e. the rate at which wind is exponentially damped by
560 : ! surface stress.
561 :
562 : ! Beljaars et al SGO scheme incorporated here. It
563 : ! appears as a "3D" tau_damp_rate specification.
564 :
565 1123584 : tau_damp_rate(:,pver) = -gravit*ksrf(:ncol)*p%rdel(:,pver)
566 5180160 : do k=1,pver
567 78723840 : tau_damp_rate(:,k) = tau_damp_rate(:,k) + dragblj(:ncol,k)
568 : end do
569 :
570 : v(:ncol,:) = fin_vol_solve(ztodt, p, v(:ncol,:), ncol, pver, &
571 : coef_q=tau_damp_rate, &
572 237149184 : coef_q_diff=kvm(:ncol,:)*dpidz_sq)
573 :
574 : u(:ncol,:) = fin_vol_solve(ztodt, p, u(:ncol,:), ncol, pver, &
575 : coef_q=tau_damp_rate, &
576 237149184 : coef_q_diff=kvm(:ncol,:)*dpidz_sq)
577 :
578 :
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 1123584 : 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 1050624 : tautmsx(i) = -ksrftms(i)*u(i,pver)
594 1050624 : 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 1050624 : taubljx(i) = 0._r8
600 1050624 : taubljy(i) = 0._r8
601 74594304 : do k = 1, pver
602 73543680 : taubljx(i) = taubljx(i) + (1._r8/gravit)*dragblj(i,k)*u(i,k)*p%del(i,k)
603 74594304 : taubljy(i) = taubljy(i) + (1._r8/gravit)*dragblj(i,k)*v(i,k)*p%del(i,k)
604 : end do
605 :
606 1123584 : if( do_iss ) then
607 :
608 : ! Compute vertical integration of final horizontal momentum
609 :
610 1050624 : usum_out(i) = 0._r8
611 1050624 : vsum_out(i) = 0._r8
612 74594304 : do k = 1, pver
613 73543680 : usum_out(i) = usum_out(i) + (1._r8/gravit)*u(i,k)*p%del(i,k)
614 74594304 : 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 1050624 : tauimpx(i) = ( usum_out(i) - usum_in(i) ) / ztodt
624 1050624 : tauimpy(i) = ( vsum_out(i) - vsum_in(i) ) / ztodt
625 :
626 1050624 : tautotx(i) = tauimpx(i)
627 1050624 : 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 1050624 : if( itaures .eq. 1 ) then
634 1050624 : tauresx(i) = taux(i) + tautmsx(i) + taubljx(i) + tauresx(i)- tauimpx(i)
635 1050624 : 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 1123584 : k = pver + 1
665 1123584 : do i = 1, ncol
666 1050624 : tmpi1(i,1) = 0._r8
667 : tmpi1(i,k) = 0.5_r8 * ztodt * gravit * &
668 1123584 : ( (-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 5107200 : do k = 2, pver
672 77600256 : do i = 1, ncol
673 72493056 : dout_u = u(i,k) - u(i,k-1)
674 72493056 : dout_v = v(i,k) - v(i,k-1)
675 72493056 : tmpi1(i,k) = 0.25_r8 * tmpi2(i,k) * kvm(i,k) * &
676 77527296 : ( 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 72960 : if (do_beljaars) then
681 :
682 : ! 2. Add Kinetic Energy change across dissipation to Static Energy
683 5180160 : do k = 1, pver
684 78723840 : do i = 1, ncol
685 78650880 : 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 5180160 : do k = 1, pver
690 78723840 : do i = 1, ncol
691 73543680 : dtk(i,k) = keg_in(i,k) - keg_out(i,k)
692 78650880 : 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 72960 : if( diffuse(fieldlist,'s') ) then
718 :
719 : ! Add counter-gradient to input static energy profiles
720 5180160 : do k = 1, pver
721 10214400 : dse(:ncol,k) = dse(:ncol,k) + ztodt * p%rdel(:,k) * gravit * &
722 5107200 : ( rhoi(:ncol,k+1) * kvh(:ncol,k+1) * cgh(:ncol,k+1) &
723 94045440 : - rhoi(:ncol,k ) * kvh(:ncol,k ) * cgh(:ncol,k ) )
724 : end do
725 :
726 : ! Add the explicit surface fluxes to the lowest layer
727 1123584 : dse(:ncol,pver) = dse(:ncol,pver) + tmp1(:ncol) * shflx(:ncol)
728 :
729 : ! Diffuse dry static energy
730 :
731 : !---------------------------------------------------
732 : ! Solve for temperature using thermal conductivity
733 : !---------------------------------------------------
734 72960 : if ( use_temperature_molec_diff ) then
735 : !----------------------------------------------------------------------------------------------------
736 : ! In Extended WACCM, kvt is calculated rather kvh. This is because molecular diffusion operates on
737 : ! temperature, while eddy diffusion operates on dse. Also, pass in constituent dependent "constants"
738 : !----------------------------------------------------------------------------------------------------
739 :
740 : ! Boundary layer thickness of "0._r8" signifies that the boundary
741 : ! condition is defined directly on the top interface.
742 :
743 : dse(:ncol,:) = fin_vol_solve(ztodt, p, dse(:ncol,:), ncol, pver, &
744 : coef_q_diff=kvh(:ncol,:)*dpidz_sq, &
745 : upper_bndry=interface_boundary, &
746 0 : l_cond=BoundaryData(dse_top(:ncol)))
747 :
748 : ! Calculate flux at top interface
749 :
750 : ! Modification : Why molecular diffusion does not work for dry static energy in all layers ?
751 :
752 : topflx(:ncol) = - kvh(:ncol,1) * tmpi2(:ncol,1) / (ztodt*gravit) * &
753 0 : ( dse(:ncol,1) - dse_top(:ncol) )
754 :
755 0 : ttemp0 = t(:ncol,:)
756 0 : ttemp = ttemp0
757 :
758 : ! upper boundary is zero flux for extended model
759 : ttemp = fin_vol_solve(ztodt, p, ttemp, ncol, pver, &
760 0 : coef_q_diff=kvt(:ncol,:)*dpidz_sq, &
761 0 : coef_q_weight=cpairv(:ncol,:))
762 :
763 :
764 : !-------------------------------------
765 : ! Update dry static energy
766 : !-------------------------------------
767 0 : do k = 1,pver
768 0 : dse(:ncol,k) = dse(:ncol,k) + &
769 0 : cpairv(:ncol,k)*(ttemp(:,k) - ttemp0(:,k))
770 : enddo
771 :
772 : else
773 :
774 72960 : if (do_molec_diff) then
775 79847424 : kv_total(:ncol,:) = kvh(:ncol,:) + kvt(:ncol,:)/cpair
776 : else
777 0 : kv_total(:ncol,:) = kvh(:ncol,:)
778 : end if
779 :
780 : ! Boundary layer thickness of "0._r8" signifies that the boundary
781 : ! condition is defined directly on the top interface.
782 : dse(:ncol,:) = fin_vol_solve(ztodt, p, dse(:ncol,:), ncol, pver, &
783 : coef_q_diff=kv_total(:ncol,:)*dpidz_sq, &
784 : upper_bndry=interface_boundary, &
785 237149184 : l_cond=BoundaryData(dse_top(:ncol)))
786 :
787 : ! Calculate flux at top interface
788 :
789 : ! Modification : Why molecular diffusion does not work for dry static energy in all layers ?
790 :
791 72960 : if( do_molec_diff ) then
792 : topflx(:ncol) = - kv_total(:ncol,1) * tmpi2(:ncol,1) / (ztodt*gravit) * &
793 1123584 : ( dse(:ncol,1) - dse_top(:ncol) )
794 : else
795 0 : topflx(:ncol) = 0._r8
796 : end if
797 :
798 : endif
799 :
800 : endif
801 :
802 : !---------------------------- !
803 : ! Diffuse Water Vapor Tracers !
804 : !---------------------------- !
805 :
806 : ! Modification : For aerosols, I need to use separate treatment
807 : ! for aerosol mass and aerosol number.
808 :
809 : ! Loop through constituents
810 :
811 : no_molec_decomp = fin_vol_lu_decomp(ztodt, p, &
812 79847424 : coef_q_diff=kvq(:ncol,:)*dpidz_sq)
813 :
814 15978240 : do m = 1, ncnst
815 :
816 15978240 : if( diffuse(fieldlist,'q',m) ) then
817 :
818 : ! Add the nonlocal transport terms to constituents in the PBL.
819 : ! Check for neg q's in each constituent and put the original vertical
820 : ! profile back if a neg value is found. A neg value implies that the
821 : ! quasi-equilibrium conditions assumed for the countergradient term are
822 : ! strongly violated.
823 :
824 6061735680 : qtm(:ncol,:pver) = q(:ncol,:pver,m)
825 :
826 398872320 : do k = 1, pver
827 393254400 : q(:ncol,k,m) = q(:ncol,k,m) + &
828 393254400 : ztodt * p%rdel(:,k) * gravit * ( cflx(:ncol,m) * rrho(:ncol) ) * &
829 393254400 : ( rhoi(:ncol,k+1) * kvh(:ncol,k+1) * cgs(:ncol,k+1) &
830 7241498880 : - rhoi(:ncol,k ) * kvh(:ncol,k ) * cgs(:ncol,k ) )
831 : end do
832 6061735680 : lqtst(:ncol) = all(q(:ncol,1:pver,m) >= qmincg(m), 2)
833 398872320 : do k = 1, pver
834 6061735680 : q(:ncol,k,m) = merge( q(:ncol,k,m), qtm(:ncol,k), lqtst(:ncol) )
835 : end do
836 :
837 : ! Add the explicit surface fluxes to the lowest layer
838 :
839 86515968 : q(:ncol,pver,m) = q(:ncol,pver,m) + tmp1(:ncol) * cflx(:ncol,m)
840 :
841 : ! Diffuse constituents.
842 :
843 : ! This is for solving molecular diffusion of minor species, thus, for WACCM-X, bypass O and O2 (major species)
844 : ! Major species diffusion is calculated separately. -Hanli Liu
845 :
846 5617920 : if( do_molec_diff .and. diffuse(fieldlistm,'q',m)) then
847 :
848 5617920 : decomp = vd_lu_qdecomp( pcols , pver , ncol , cnst_fixed_ubc(m), cnst_mw(m), &
849 : kvq , kq_scal, mw_fac(:,:,m) ,dpidz_sq , p_molec, &
850 : interface_boundary, molec_boundary, &
851 : tint , ztodt , nbot_molec , &
852 5617920 : lchnk , t , m , no_molec_decomp)
853 :
854 : ! This to calculate the upper boundary flux of H. -Hanli Liu
855 5617920 : if ((cnst_fixed_ubflx(m))) then
856 :
857 : ! ubc_flux is a flux of mass density through space, i.e.:
858 : ! ubc_flux = rho_i * dz/dt = q_i * rho * dz/dt
859 : ! For flux of mmr through pressure level, multiply by g:
860 : ! q_i * rho * gravit * dz/dt = q_i * dp/dt
861 :
862 : call decomp%left_div(q(:ncol,:,m), &
863 : l_cond=BoundaryFlux( &
864 0 : -gravit*ubc_flux(:ncol,m), ztodt, &
865 0 : p%del(:,1)))
866 :
867 : else
868 : call decomp%left_div(q(:ncol,:,m), &
869 5617920 : l_cond=BoundaryData(ubc_mmr(:ncol,m)))
870 : end if
871 :
872 5617920 : call decomp%finalize()
873 :
874 : else
875 :
876 0 : if (present(cnst_fixed_ubc)) then
877 : ! explicitly set mmr in top layer for cases where molecular diffusion is not active
878 0 : if (cnst_fixed_ubc(m)) then
879 0 : q(:ncol,1,m) = ubc_mmr(:ncol,m)
880 : endif
881 : end if
882 :
883 0 : call no_molec_decomp%left_div(q(:ncol,:,m))
884 :
885 : end if
886 :
887 : end if
888 : end do
889 :
890 72960 : call no_molec_decomp%finalize()
891 :
892 656640 : end subroutine compute_vdiff
893 :
894 : ! =============================================================================== !
895 : ! !
896 : ! =============================================================================== !
897 :
898 241152 : character(128) function vdiff_select( fieldlist, name, qindex )
899 : ! --------------------------------------------------------------------- !
900 : ! This function sets the field with incoming name as one to be diffused !
901 : ! --------------------------------------------------------------------- !
902 : type(vdiff_selector), intent(inout) :: fieldlist
903 : character(*), intent(in) :: name
904 : integer, intent(in), optional :: qindex
905 :
906 241152 : vdiff_select = ''
907 1536 : select case (name)
908 : case ('u','U')
909 1536 : fieldlist%fields(1) = .true.
910 : case ('v','V')
911 1536 : fieldlist%fields(2) = .true.
912 : case ('s','S')
913 1536 : fieldlist%fields(3) = .true.
914 : case ('q','Q')
915 236544 : if( present(qindex) ) then
916 236544 : fieldlist%fields(3 + qindex) = .true.
917 : else
918 0 : fieldlist%fields(4) = .true.
919 : endif
920 : case default
921 241152 : write(vdiff_select,*) 'Bad argument to vdiff_index: ', name
922 : end select
923 241152 : return
924 :
925 : end function vdiff_select
926 :
927 0 : type(vdiff_selector) function not(a)
928 : ! ------------------------------------------------------------- !
929 : ! This function extends .not. to operate on type vdiff_selector !
930 : ! ------------------------------------------------------------- !
931 : type(vdiff_selector), intent(in) :: a
932 0 : allocate(not%fields(size(a%fields)))
933 0 : not%fields = .not. a%fields
934 0 : end function not
935 :
936 145920 : logical function my_any(a)
937 : ! -------------------------------------------------- !
938 : ! This function extends the intrinsic function 'any' !
939 : ! to operate on type vdiff_selector !
940 : ! -------------------------------------------------- !
941 : type(vdiff_selector), intent(in) :: a
942 16270080 : my_any = any(a%fields)
943 145920 : end function my_any
944 :
945 21960960 : logical function diffuse(fieldlist,name,qindex)
946 : ! ---------------------------------------------------------------------------- !
947 : ! This function reports whether the field with incoming name is to be diffused !
948 : ! ---------------------------------------------------------------------------- !
949 : type(vdiff_selector), intent(in) :: fieldlist
950 : character(*), intent(in) :: name
951 : integer, intent(in), optional :: qindex
952 :
953 145920 : select case (name)
954 : case ('u','U')
955 145920 : diffuse = fieldlist%fields(1)
956 : case ('v','V')
957 145920 : diffuse = fieldlist%fields(2)
958 : case ('s','S')
959 145920 : diffuse = fieldlist%fields(3)
960 : case ('q','Q')
961 21523200 : if( present(qindex) ) then
962 21523200 : diffuse = fieldlist%fields(3 + qindex)
963 : else
964 0 : diffuse = fieldlist%fields(4)
965 : endif
966 : case default
967 21960960 : diffuse = .false.
968 : end select
969 : return
970 : end function diffuse
971 :
972 0 : end module diffusion_solver
|
