Line data Source code
1 : module gw_common
2 :
3 : !
4 : ! This module contains code common to different gravity wave
5 : ! parameterizations.
6 : !
7 : use gw_utils, only: r8
8 : use coords_1d, only: Coords1D
9 :
10 :
11 : implicit none
12 : private
13 : save
14 :
15 : ! Public interface.
16 :
17 : public :: GWBand
18 :
19 : public :: gw_common_init
20 : public :: gw_prof
21 : public :: gw_drag_prof
22 : public :: qbo_hdepth_scaling
23 : public :: calc_taucd, momentum_flux, momentum_fixer
24 : public :: energy_change, energy_fixer
25 : public :: coriolis_speed, adjust_inertial
26 :
27 : public :: pver
28 : public :: west, east, north, south
29 : public :: pi
30 : public :: gravit
31 : public :: rair
32 :
33 : ! Number of levels in the atmosphere.
34 : integer, protected :: pver = 0
35 :
36 : ! Whether or not to enforce an upper boundary condition of tau = 0.
37 : logical :: tau_0_ubc = .false.
38 :
39 : ! Index the cardinal directions.
40 : integer, parameter :: west = 1
41 : integer, parameter :: east = 2
42 : integer, parameter :: south = 3
43 : integer, parameter :: north = 4
44 :
45 : ! Scaling factor for generating QBO
46 : real(r8), protected :: qbo_hdepth_scaling
47 :
48 : ! 3.14159...
49 : real(r8), parameter :: pi = acos(-1._r8)
50 :
51 : ! Acceleration due to gravity.
52 : real(r8), protected :: gravit = huge(1._r8)
53 :
54 : ! Gas constant for dry air.
55 : real(r8), protected :: rair = huge(1._r8)
56 :
57 : !
58 : ! Private variables
59 : !
60 :
61 : ! Interface levels for gravity wave sources.
62 : integer :: ktop = huge(1)
63 :
64 : ! Background diffusivity.
65 : real(r8), parameter :: dback = 0.05_r8
66 :
67 : ! rair/gravit
68 : real(r8) :: rog = huge(1._r8)
69 :
70 : ! Newtonian cooling coefficients.
71 : real(r8), allocatable :: alpha(:)
72 :
73 : ! Inverse Prandtl number.
74 : real(r8) :: prndl
75 :
76 : !
77 : ! Limits to keep values reasonable.
78 : !
79 :
80 : ! Minimum non-zero stress.
81 : real(r8), parameter :: taumin = 1.e-10_r8
82 : ! Maximum wind tendency from stress divergence (before efficiency applied).
83 : ! 400 m/s/day
84 : real(r8), parameter :: tndmax = 400._r8 / 86400._r8
85 : ! Maximum allowed change in u-c (before efficiency applied).
86 : real(r8), parameter :: umcfac = 0.5_r8
87 : ! Minimum value of (u-c)**2.
88 : real(r8), parameter :: ubmc2mn = 0.01_r8
89 :
90 : ! Type describing a band of wavelengths into which gravity waves can be
91 : ! emitted.
92 : ! Currently this has to have uniform spacing (i.e. adjacent elements of
93 : ! cref are exactly dc apart).
94 : type :: GWBand
95 : ! Dimension of the spectrum.
96 : integer :: ngwv
97 : ! Delta between nearest phase speeds [m/s].
98 : real(r8) :: dc
99 : ! Reference speeds [m/s].
100 : real(r8), allocatable :: cref(:)
101 : ! Critical Froude number, squared
102 : real(r8) :: fcrit2
103 : ! Horizontal wave number [1/m].
104 : real(r8) :: kwv
105 : ! Effective horizontal wave number [1/m] (fcrit2*kwv).
106 : real(r8) :: effkwv
107 : end type GWBand
108 :
109 : interface GWBand
110 : module procedure new_GWBand
111 : end interface
112 :
113 : contains
114 :
115 : !==========================================================================
116 :
117 : ! Constructor for a GWBand that calculates derived components.
118 6144 : function new_GWBand(ngwv, dc, fcrit2, wavelength) result(band)
119 : ! Used directly to set the type's components.
120 : integer, intent(in) :: ngwv
121 : real(r8), intent(in) :: dc
122 : real(r8), intent(in) :: fcrit2
123 : ! Wavelength in meters.
124 : real(r8), intent(in) :: wavelength
125 :
126 : ! Output.
127 : type(GWBand) :: band
128 :
129 : ! Wavenumber index.
130 : integer :: l
131 :
132 : ! Simple assignments.
133 6144 : band%ngwv = ngwv
134 6144 : band%dc = dc
135 6144 : band%fcrit2 = fcrit2
136 :
137 : ! Uniform phase speed reference grid.
138 18432 : allocate(band%cref(-ngwv:ngwv))
139 30720 : band%cref = [( dc * l, l = -ngwv, ngwv )]
140 :
141 : ! Wavenumber and effective wavenumber come from the wavelength.
142 6144 : band%kwv = 2._r8*pi / wavelength
143 6144 : band%effkwv = band%fcrit2 * band%kwv
144 :
145 6144 : end function new_GWBand
146 :
147 : !==========================================================================
148 :
149 1536 : subroutine gw_common_init(pver_in, &
150 3072 : tau_0_ubc_in, ktop_in, gravit_in, rair_in, alpha_in, &
151 : prndl_in, qbo_hdepth_scaling_in, errstring)
152 :
153 : integer, intent(in) :: pver_in
154 : logical, intent(in) :: tau_0_ubc_in
155 : integer, intent(in) :: ktop_in
156 : real(r8), intent(in) :: gravit_in
157 : real(r8), intent(in) :: rair_in
158 : real(r8), intent(in) :: alpha_in(:)
159 : real(r8), intent(in) :: prndl_in
160 : real(r8), intent(in) :: qbo_hdepth_scaling_in
161 : ! Report any errors from this routine.
162 : character(len=*), intent(out) :: errstring
163 :
164 : integer :: ierr
165 :
166 1536 : errstring = ""
167 :
168 1536 : pver = pver_in
169 1536 : tau_0_ubc = tau_0_ubc_in
170 1536 : ktop = ktop_in
171 1536 : gravit = gravit_in
172 1536 : rair = rair_in
173 4608 : allocate(alpha(pver+1), stat=ierr, errmsg=errstring)
174 0 : if (ierr /= 0) return
175 44544 : alpha = alpha_in
176 1536 : prndl = prndl_in
177 1536 : qbo_hdepth_scaling = qbo_hdepth_scaling_in
178 :
179 1536 : rog = rair/gravit
180 :
181 1536 : end subroutine gw_common_init
182 :
183 : !==========================================================================
184 :
185 1489176 : subroutine gw_prof (ncol, p, cpair, t, rhoi, nm, ni)
186 : !-----------------------------------------------------------------------
187 : ! Compute profiles of background state quantities for the multiple
188 : ! gravity wave drag parameterization.
189 : !
190 : ! The parameterization is assumed to operate only where water vapor
191 : ! concentrations are negligible in determining the density.
192 : !-----------------------------------------------------------------------
193 : use gw_utils, only: midpoint_interp
194 : !------------------------------Arguments--------------------------------
195 : ! Column dimension.
196 : integer, intent(in) :: ncol
197 : ! Pressure coordinates.
198 : type(Coords1D), intent(in) :: p
199 :
200 : ! Specific heat of dry air, constant pressure.
201 : real(r8), intent(in) :: cpair
202 : ! Midpoint temperatures.
203 : real(r8), intent(in) :: t(ncol,pver)
204 :
205 : ! Interface density.
206 : real(r8), intent(out) :: rhoi(ncol,pver+1)
207 : ! Midpoint and interface Brunt-Vaisalla frequencies.
208 : real(r8), intent(out) :: nm(ncol,pver), ni(ncol,pver+1)
209 :
210 : !---------------------------Local Storage-------------------------------
211 : ! Column and level indices.
212 : integer :: i,k
213 :
214 : ! dt/dp
215 : real(r8) :: dtdp
216 : ! Brunt-Vaisalla frequency squared.
217 : real(r8) :: n2
218 :
219 : ! Interface temperature.
220 2978352 : real(r8) :: ti(ncol,pver+1)
221 :
222 : ! Minimum value of Brunt-Vaisalla frequency squared.
223 : real(r8), parameter :: n2min = 5.e-5_r8
224 :
225 : !------------------------------------------------------------------------
226 : ! Determine the interface densities and Brunt-Vaisala frequencies.
227 : !------------------------------------------------------------------------
228 :
229 : ! The top interface values are calculated assuming an isothermal
230 : ! atmosphere above the top level.
231 1489176 : k = 1
232 24865776 : do i = 1, ncol
233 23376600 : ti(i,k) = t(i,k)
234 23376600 : rhoi(i,k) = p%ifc(i,k) / (rair*ti(i,k))
235 24865776 : ni(i,k) = sqrt(gravit*gravit / (cpair*ti(i,k)))
236 : end do
237 :
238 : ! Interior points use centered differences.
239 1489176 : ti(:,2:pver) = midpoint_interp(t)
240 38718576 : do k = 2, pver
241 623133576 : do i = 1, ncol
242 584415000 : rhoi(i,k) = p%ifc(i,k) / (rair*ti(i,k))
243 584415000 : dtdp = (t(i,k)-t(i,k-1)) * p%rdst(i,k-1)
244 584415000 : n2 = gravit*gravit/ti(i,k) * (1._r8/cpair - rhoi(i,k)*dtdp)
245 621644400 : ni(i,k) = sqrt(max(n2min, n2))
246 : end do
247 : end do
248 :
249 : ! Bottom interface uses bottom level temperature, density; next interface
250 : ! B-V frequency.
251 1489176 : k = pver+1
252 24865776 : do i = 1, ncol
253 23376600 : ti(i,k) = t(i,k-1)
254 23376600 : rhoi(i,k) = p%ifc(i,k) / (rair*ti(i,k))
255 24865776 : ni(i,k) = ni(i,k-1)
256 : end do
257 :
258 : !------------------------------------------------------------------------
259 : ! Determine the midpoint Brunt-Vaisala frequencies.
260 : !------------------------------------------------------------------------
261 1489176 : nm = midpoint_interp(ni)
262 :
263 1489176 : end subroutine gw_prof
264 :
265 : !==========================================================================
266 :
267 1489176 : subroutine gw_drag_prof(ncol, band, p, src_level, tend_level, dt, &
268 1489176 : t, vramp, &
269 5956704 : piln, rhoi, nm, ni, ubm, ubi, xv, yv, &
270 7445880 : effgw, c, kvtt, q, dse, tau, utgw, vtgw, &
271 2978352 : ttgw, qtgw, egwdffi, gwut, dttdf, dttke, ro_adjust, &
272 0 : kwvrdg, satfac_in, lapply_effgw_in, lapply_vdiff, tau_diag )
273 :
274 : !-----------------------------------------------------------------------
275 : ! Solve for the drag profile from the multiple gravity wave drag
276 : ! parameterization.
277 : ! 1. scan up from the wave source to determine the stress profile
278 : ! 2. scan down the stress profile to determine the tendencies
279 : ! => apply bounds to the tendency
280 : ! a. from wkb solution
281 : ! b. from computational stability constraints
282 : ! => adjust stress on interface below to reflect actual bounded
283 : ! tendency
284 : !-----------------------------------------------------------------------
285 :
286 : use gw_diffusion, only: gw_ediff, gw_diff_tend
287 : use linear_1d_operators, only: TriDiagDecomp
288 :
289 : !------------------------------Arguments--------------------------------
290 : ! Column dimension.
291 : integer, intent(in) :: ncol
292 : ! Wavelengths.
293 : type(GWBand), intent(in) :: band
294 : ! Pressure coordinates.
295 : type(Coords1D), intent(in) :: p
296 : ! Level from which gravity waves are propagated upward.
297 : integer, intent(in) :: src_level(ncol)
298 : ! Lowest level where wind tendencies are calculated.
299 : integer, intent(in) :: tend_level(ncol)
300 : ! Using tend_level > src_level allows the orographic waves to prescribe
301 : ! wave propagation up to a certain level, but then allow wind tendencies
302 : ! and adjustments to tau below that level.
303 :
304 : ! Time step.
305 : real(r8), intent(in) :: dt
306 :
307 : ! Midpoint and interface temperatures.
308 : real(r8), intent(in) :: t(ncol,pver)
309 : ! Log of interface pressures.
310 : real(r8), intent(in) :: piln(ncol,pver+1)
311 : ! Interface densities.
312 : real(r8), intent(in) :: rhoi(ncol,pver+1)
313 : ! Midpoint and interface Brunt-Vaisalla frequencies.
314 : real(r8), intent(in) :: nm(ncol,pver), ni(ncol,pver+1)
315 : ! Projection of wind at midpoints and interfaces.
316 : real(r8), intent(in) :: ubm(ncol,pver), ubi(ncol,pver+1)
317 : ! Unit vectors of source wind (zonal and meridional components).
318 : real(r8), intent(in) :: xv(ncol), yv(ncol)
319 : ! Tendency efficiency.
320 : real(r8), intent(in) :: effgw(ncol)
321 : ! Wave phase speeds for each column.
322 : real(r8), intent(in) :: c(ncol,-band%ngwv:band%ngwv)
323 : ! Molecular thermal diffusivity.
324 : real(r8), intent(in) :: kvtt(ncol,pver+1)
325 : ! Constituent array.
326 : real(r8), intent(in) :: q(:,:,:)
327 : ! Dry static energy.
328 : real(r8), intent(in) :: dse(ncol,pver)
329 : ! Coefficient to ramp down diffusion coeff.
330 : real(r8), pointer, intent(in) :: vramp(:)
331 :
332 : ! Wave Reynolds stress.
333 : real(r8), intent(inout) :: tau(ncol,-band%ngwv:band%ngwv,pver+1)
334 : ! Zonal/meridional wind tendencies.
335 : real(r8), intent(out) :: utgw(ncol,pver), vtgw(ncol,pver)
336 : ! Gravity wave heating tendency.
337 : real(r8), intent(out) :: ttgw(ncol,pver)
338 : ! Gravity wave constituent tendency.
339 : real(r8), intent(out) :: qtgw(:,:,:)
340 :
341 : ! Effective gravity wave diffusivity at interfaces.
342 : real(r8), intent(out) :: egwdffi(ncol,pver+1)
343 :
344 : ! Gravity wave wind tendency for each wave.
345 : real(r8), intent(out) :: gwut(ncol,pver,-band%ngwv:band%ngwv)
346 :
347 : ! Temperature tendencies from diffusion and kinetic energy.
348 : real(r8), intent(out) :: dttdf(ncol,pver)
349 : real(r8), intent(out) :: dttke(ncol,pver)
350 :
351 : ! Adjustment parameter for IGWs.
352 : real(r8), intent(in), optional :: &
353 : ro_adjust(ncol,-band%ngwv:band%ngwv,pver+1)
354 :
355 : ! Diagnosed horizontal wavenumber for ridges.
356 : real(r8), intent(in), optional :: &
357 : kwvrdg(ncol)
358 :
359 : ! Factor for saturation calculation. Here backwards
360 : ! compatibility. I believe it should be 1.0 (jtb).
361 : ! Looks like it has been 2.0 for a while in CAM.
362 : real(r8), intent(in), optional :: &
363 : satfac_in
364 :
365 : logical, intent(in), optional :: lapply_effgw_in, lapply_vdiff
366 : ! Provisional Wave Reynolds stress.
367 : real(r8), intent(out), optional :: tau_diag(ncol,pver+1)
368 :
369 : !---------------------------Local storage-------------------------------
370 :
371 : ! Level, wavenumber, constituent and column loop indices.
372 : integer :: k, l, m, i
373 :
374 : ! Lowest tendency and source levels.
375 : integer :: kbot_tend, kbot_src
376 :
377 : ! "Total" and saturation diffusivity.
378 2978352 : real(r8) :: d(ncol)
379 : ! Imaginary part of vertical wavenumber.
380 2978352 : real(r8) :: mi(ncol)
381 : ! Stress after damping.
382 2978352 : real(r8) :: taudmp(ncol)
383 : ! Saturation stress.
384 2978352 : real(r8) :: tausat(ncol)
385 : ! (ub-c) and (ub-c)**2
386 2978352 : real(r8) :: ubmc(ncol), ubmc2(ncol)
387 : ! Temporary ubar tendencies (overall, and at wave l).
388 2978352 : real(r8) :: ubt(ncol,pver), ubtl(ncol)
389 2978352 : real(r8) :: wrk(ncol)
390 : ! Ratio used for ubt tndmax limiting.
391 2978352 : real(r8) :: ubt_lim_ratio(ncol)
392 :
393 : ! saturation factor. Defaults to 2.0
394 : ! unless overidden by satfac_in
395 : real(r8) :: satfac
396 :
397 : logical :: lapply_effgw,do_vertical_diffusion
398 :
399 : ! LU decomposition.
400 1489176 : type(TriDiagDecomp) :: decomp
401 :
402 : !------------------------------------------------------------------------
403 :
404 1489176 : if (present(satfac_in)) then
405 0 : satfac = satfac_in
406 : else
407 : satfac = 2._r8
408 : endif
409 :
410 : ! Default behavior is to apply vertical diffusion.
411 : ! The user has the option to turn off vert diffusion
412 1489176 : do_vertical_diffusion = .true.
413 1489176 : if (present(lapply_vdiff)) then
414 0 : do_vertical_diffusion = lapply_vdiff
415 : endif
416 :
417 : ! Default behavior is to apply effgw and
418 : ! tendency limiters as designed by Sean
419 : ! Santos (lapply_effgw=.TRUE.). However,
420 : ! WACCM non-oro GW need to be retuned before
421 : ! this can done to them. --jtb 03/02/16
422 1489176 : if (present(lapply_effgw_in)) then
423 1489176 : lapply_effgw = lapply_effgw_in
424 : else
425 : lapply_effgw = .TRUE.
426 : endif
427 :
428 :
429 : ! Lowest levels that loops need to iterate over.
430 24865776 : kbot_tend = maxval(tend_level)
431 24865776 : kbot_src = maxval(src_level)
432 :
433 : ! Initialize gravity wave drag tendencies to zero.
434 :
435 647999352 : utgw = 0._r8
436 647999352 : vtgw = 0._r8
437 :
438 649488528 : gwut = 0._r8
439 :
440 647999352 : dttke = 0._r8
441 647999352 : ttgw = 0._r8
442 :
443 647999352 : dttdf = 0._r8
444 1945487232 : qtgw = 0._r8
445 :
446 : ! Workaround floating point exception issues on Intel by initializing
447 : ! everything that's first set in a where block.
448 24865776 : mi = 0._r8
449 24865776 : taudmp = 0._r8
450 24865776 : tausat = 0._r8
451 24865776 : ubmc = 0._r8
452 24865776 : ubmc2 = 0._r8
453 24865776 : wrk = 0._r8
454 :
455 : !------------------------------------------------------------------------
456 : ! Compute the stress profiles and diffusivities
457 : !------------------------------------------------------------------------
458 :
459 : ! Loop from bottom to top to get stress profiles.
460 : ! do k = kbot_src-1, ktop, -1 !++jtb I think this is right
461 40207752 : do k = kbot_src, ktop, -1 !++ but this is in model now
462 :
463 : ! Determine the diffusivity for each column.
464 :
465 646510176 : d = dback + kvtt(:,k)
466 :
467 78926328 : do l = -band%ngwv, band%ngwv
468 :
469 : ! Determine the absolute value of the saturation stress.
470 : ! Define critical levels where the sign of (u-c) changes between
471 : ! interfaces.
472 646510176 : ubmc = ubi(:,k) - c(:,l)
473 :
474 646510176 : tausat = 0.0_r8
475 :
476 38718576 : if (present(kwvrdg)) then
477 0 : where (src_level >= k)
478 : ! Test to see if u-c has the same sign here as the level below.
479 0 : where (ubmc > 0.0_r8 .eqv. ubi(:,k+1) > c(:,l))
480 0 : tausat = abs( kwvrdg * rhoi(:,k) * ubmc**3 / &
481 : (satfac*ni(:,k)))
482 : end where
483 : end where
484 : else
485 1900811952 : where (src_level >= k)
486 : ! Test to see if u-c has the same sign here as the level below.
487 38718576 : where (ubmc > 0.0_r8 .eqv. ubi(:,k+1) > c(:,l))
488 : tausat = abs(band%effkwv * rhoi(:,k) * ubmc**3 / &
489 : (satfac*ni(:,k)))
490 : end where
491 : end where
492 : end if
493 :
494 38718576 : if (present(ro_adjust)) then
495 0 : where (src_level >= k)
496 0 : tausat = tausat * sqrt(ro_adjust(:,l,k))
497 : end where
498 : end if
499 :
500 77437152 : if (present(kwvrdg)) then
501 0 : where (src_level >= k)
502 : ! Compute stress for each wave. The stress at this level is the
503 : ! min of the saturation stress and the stress at the level below
504 : ! reduced by damping. The sign of the stress must be the same as
505 : ! at the level below.
506 :
507 : ubmc2 = max(ubmc**2, ubmc2mn)
508 0 : mi = ni(:,k) / (2._r8 * kwvrdg * ubmc2) * & ! Is this 2._r8 related to satfac?
509 0 : (alpha(k) + ni(:,k)**2/ubmc2 * d)
510 0 : wrk = -2._r8*mi*rog*t(:,k)*(piln(:,k+1) - piln(:,k))
511 :
512 : taudmp = tau(:,l,k+1)
513 :
514 : ! For some reason, PGI 14.1 loses bit-for-bit reproducibility if
515 : ! we limit tau, so instead limit the arrays used to set it.
516 : where (tausat <= taumin) tausat = 0._r8
517 : where (taudmp <= taumin) taudmp = 0._r8
518 :
519 : tau(:,l,k) = min(taudmp, tausat)
520 : end where
521 :
522 : else
523 :
524 6194071728 : where (src_level >= k)
525 :
526 : ! Compute stress for each wave. The stress at this level is the
527 : ! min of the saturation stress and the stress at the level below
528 : ! reduced by damping. The sign of the stress must be the same as
529 : ! at the level below.
530 :
531 : ubmc2 = max(ubmc**2, ubmc2mn)
532 : mi = ni(:,k) / (2._r8 * band%kwv * ubmc2) * &
533 38718576 : (alpha(k) + ni(:,k)**2/ubmc2 * d)
534 38718576 : wrk = -2._r8*mi*rog*t(:,k)*(piln(:,k+1) - piln(:,k))
535 :
536 : taudmp = tau(:,l,k+1) * exp(wrk)
537 :
538 : ! For some reason, PGI 14.1 loses bit-for-bit reproducibility if
539 : ! we limit tau, so instead limit the arrays used to set it.
540 : where (tausat <= taumin) tausat = 0._r8
541 : where (taudmp <= taumin) taudmp = 0._r8
542 :
543 : tau(:,l,k) = min(taudmp, tausat)
544 : end where
545 : endif
546 :
547 : end do
548 : end do
549 :
550 : ! Force tau at the top of the model to zero, if requested.
551 1489176 : if (tau_0_ubc) tau(:,:,ktop) = 0._r8
552 :
553 : ! Write out pre-adjustment tau profile for diagnostc purposes.
554 : ! Current implementation only makes sense for orographic waves.
555 : ! Fix later.
556 1489176 : if (PRESENT(tau_diag)) then
557 0 : tau_diag(:,:) = tau(:,0,:)
558 : end if
559 :
560 : ! Apply efficiency to completed stress profile.
561 1489176 : if (lapply_effgw) then
562 41696928 : do k = ktop, kbot_tend+1
563 81904680 : do l = -band%ngwv, band%ngwv
564 711583704 : where (k-1 <= tend_level)
565 40207752 : tau(:,l,k) = tau(:,l,k) * effgw
566 : end where
567 : end do
568 : end do
569 : end if
570 :
571 : !------------------------------------------------------------------------
572 : ! Compute the tendencies from the stress divergence.
573 : !------------------------------------------------------------------------
574 :
575 : ! Loop over levels from top to bottom
576 40207752 : do k = ktop, kbot_tend
577 :
578 : ! Accumulate the mean wind tendency over wavenumber.
579 646510176 : ubt(:,k) = 0.0_r8
580 :
581 77437152 : do l = -band%ngwv, band%ngwv ! loop over wave
582 :
583 : ! Determine the wind tendency, including excess stress carried down
584 : ! from above.
585 646510176 : ubtl = gravit * (tau(:,l,k+1)-tau(:,l,k)) * p%rdel(:,k)
586 :
587 : ! Apply first tendency limit to maintain numerical stability.
588 : ! Enforce du/dt < |c-u|/dt so u-c cannot change sign
589 : ! (u^n+1 = u^n + du/dt * dt)
590 : ! The limiter is somewhat stricter, so that we don't come anywhere
591 : ! near reversing c-u.
592 646510176 : ubtl = min(ubtl, umcfac * abs(c(:,l)-ubm(:,k)) / dt)
593 :
594 38718576 : if (.not. lapply_effgw) ubtl = min(ubtl, tndmax)
595 :
596 2508603552 : where (k <= tend_level)
597 :
598 : ! Save tendency for each wave (for later computation of kzz).
599 : ! sign function returns magnitude of ubtl with sign of c-ubm
600 : ! Renders ubt/ubm check for mountain waves unecessary
601 : gwut(:,k,l) = sign(ubtl, c(:,l)-ubm(:,k))
602 : ubt(:,k) = ubt(:,k) + gwut(:,k,l)
603 :
604 : end where
605 :
606 : end do
607 :
608 38718576 : if (lapply_effgw) then
609 : ! Apply second tendency limit to maintain numerical stability.
610 : ! Enforce du/dt < tndmax so that ridicuously large tendencies are not
611 : ! permitted.
612 : ! This can only happen above tend_level, so don't bother checking the
613 : ! level explicitly.
614 3077676576 : where (abs(ubt(:,k)) > tndmax)
615 : ubt_lim_ratio = tndmax/abs(ubt(:,k))
616 : ubt(:,k) = ubt_lim_ratio * ubt(:,k)
617 : elsewhere
618 : ubt_lim_ratio = 1._r8
619 : end where
620 : else
621 0 : ubt_lim_ratio = 1._r8
622 : end if
623 :
624 77437152 : do l = -band%ngwv, band%ngwv
625 646510176 : gwut(:,k,l) = ubt_lim_ratio*gwut(:,k,l)
626 : ! Redetermine the effective stress on the interface below from the
627 : ! wind tendency. If the wind tendency was limited above, then the
628 : ! new stress will be smaller than the old stress, causing stress
629 : ! divergence in the next layer down. This smoothes large stress
630 : ! divergences downward while conserving total stress.
631 :
632 : ! Protection on SMALL gwut to prevent floating point
633 : ! issues.
634 : !--------------------------------------------------
635 646510176 : where( abs(gwut(:,k,l)) < 1.e-15_r8 )
636 : gwut(:,k,l) = 0._r8
637 : endwhere
638 :
639 685228752 : where (k <= tend_level)
640 77437152 : tau(:,l,k+1) = tau(:,l,k) + &
641 77437152 : abs(gwut(:,k,l)) * p%del(:,k) / gravit
642 : end where
643 : end do
644 :
645 : ! Project the mean wind tendency onto the components.
646 1862093376 : where (k <= tend_level)
647 : utgw(:,k) = ubt(:,k) * xv
648 : vtgw(:,k) = ubt(:,k) * yv
649 : end where
650 :
651 40207752 : if (associated(vramp)) then
652 0 : utgw(:,k) = utgw(:,k) * vramp(k)
653 0 : vtgw(:,k) = vtgw(:,k) * vramp(k)
654 : endif
655 :
656 : ! End of level loop.
657 : end do
658 :
659 :
660 : ! Block to undo Sean Santos mods to effgw and limiters.
661 : ! Here because non-oro GW in WACCM need extensive re-tuning
662 : ! before Sean's mods can be adopted. --jtb 03/02/16
663 : !==========================================
664 1489176 : if (.not.(lapply_effgw)) then
665 0 : do k = ktop, kbot_tend+1
666 0 : do l = -band%ngwv, band%ngwv
667 0 : where (k-1 <= tend_level)
668 0 : tau(:,l,k) = tau(:,l,k) * effgw
669 : end where
670 : end do
671 : end do
672 0 : do k = ktop, kbot_tend
673 0 : do l = -band%ngwv, band%ngwv
674 0 : gwut(:,k,l) = gwut(:,k,l) * effgw
675 : end do
676 0 : utgw(:,k) = utgw(:,k) * effgw
677 0 : vtgw(:,k) = vtgw(:,k) * effgw
678 : end do
679 : end if
680 : !===========================================
681 :
682 1489176 : if (do_vertical_diffusion) then
683 :
684 : ! Calculate effective diffusivity and LU decomposition for the
685 : ! vertical diffusion solver.
686 : call gw_ediff (ncol, pver, band%ngwv, kbot_tend, ktop, tend_level, &
687 : gwut, ubm, nm, rhoi, dt, prndl, gravit, p, c, vramp, &
688 1489176 : egwdffi, decomp, ro_adjust=ro_adjust)
689 :
690 : ! Calculate tendency on each constituent.
691 5956704 : do m = 1, size(q,3)
692 :
693 0 : call gw_diff_tend(ncol, pver, kbot_tend, ktop, q(:,:,m), &
694 5956704 : dt, decomp, qtgw(:,:,m))
695 :
696 : enddo
697 :
698 : ! Calculate tendency from diffusing dry static energy (dttdf).
699 1489176 : call gw_diff_tend(ncol, pver, kbot_tend, ktop, dse, dt, decomp, dttdf)
700 :
701 : endif
702 :
703 : ! Evaluate second temperature tendency term: Conversion of kinetic
704 : ! energy into thermal.
705 2978352 : do l = -band%ngwv, band%ngwv
706 41696928 : do k = ktop, kbot_tend
707 647999352 : dttke(:,k) = dttke(:,k) - (ubm(:,k) - c(:,l)) * gwut(:,k,l)
708 : end do
709 : end do
710 :
711 647999352 : ttgw = dttke + dttdf
712 :
713 1489176 : if (associated(vramp)) then
714 0 : do k = ktop, kbot_tend
715 0 : ttgw(:,k) = ttgw(:,k) * vramp(k)
716 : enddo
717 : endif
718 :
719 : ! Deallocate decomp.
720 1489176 : call decomp%finalize()
721 :
722 2978352 : end subroutine gw_drag_prof
723 :
724 : !==========================================================================
725 :
726 : ! Calculate Reynolds stress for waves propagating in each cardinal
727 : ! direction.
728 :
729 0 : function calc_taucd(ncol, ngwv, tend_level, tau, c, xv, yv, ubi) &
730 0 : result(taucd)
731 :
732 : ! Column and gravity wave wavenumber dimensions.
733 : integer, intent(in) :: ncol, ngwv
734 : ! Lowest level where wind tendencies are calculated.
735 : integer, intent(in) :: tend_level(:)
736 : ! Wave Reynolds stress.
737 : real(r8), intent(in) :: tau(:,-ngwv:,:)
738 : ! Wave phase speeds for each column.
739 : real(r8), intent(in) :: c(:,-ngwv:)
740 : ! Unit vectors of source wind (zonal and meridional components).
741 : real(r8), intent(in) :: xv(:), yv(:)
742 : ! Projection of wind at interfaces.
743 : real(r8), intent(in) :: ubi(:,:)
744 :
745 : real(r8) :: taucd(ncol,pver+1,4)
746 :
747 : ! Indices.
748 : integer :: i, k, l
749 :
750 : ! ubi at tend_level.
751 0 : real(r8) :: ubi_tend(ncol)
752 :
753 : ! Signed wave Reynolds stress.
754 0 : real(r8) :: tausg(ncol)
755 :
756 : ! Reynolds stress for waves propagating behind and forward of the wind.
757 0 : real(r8) :: taub(ncol)
758 0 : real(r8) :: tauf(ncol)
759 :
760 0 : taucd = 0._r8
761 0 : tausg = 0._r8
762 :
763 0 : ubi_tend = (/ (ubi(i,tend_level(i)+1), i = 1, ncol) /)
764 :
765 0 : do k = ktop, maxval(tend_level)+1
766 :
767 0 : taub = 0._r8
768 0 : tauf = 0._r8
769 :
770 0 : do l = -ngwv, ngwv
771 0 : where (k-1 <= tend_level)
772 :
773 0 : tausg = sign(tau(:,l,k), c(:,l)-ubi(:,k))
774 :
775 0 : where ( c(:,l) < ubi_tend )
776 : taub = taub + tausg
777 : elsewhere
778 : tauf = tauf + tausg
779 : end where
780 :
781 : end where
782 : end do
783 :
784 0 : where (k-1 <= tend_level)
785 : where (xv > 0._r8)
786 0 : taucd(:,k,east) = tauf * xv
787 : taucd(:,k,west) = taub * xv
788 : elsewhere
789 : taucd(:,k,east) = taub * xv
790 : taucd(:,k,west) = tauf * xv
791 : end where
792 :
793 : where ( yv > 0._r8)
794 0 : taucd(:,k,north) = tauf * yv
795 : taucd(:,k,south) = taub * yv
796 : elsewhere
797 : taucd(:,k,north) = taub * yv
798 : taucd(:,k,south) = tauf * yv
799 : end where
800 : end where
801 :
802 : end do
803 :
804 0 : end function calc_taucd
805 :
806 : !==========================================================================
807 :
808 : ! Calculate the amount of momentum conveyed from below the gravity wave
809 : ! region, to the region where gravity waves are calculated.
810 0 : subroutine momentum_flux(tend_level, taucd, um_flux, vm_flux)
811 :
812 : ! Bottom stress level.
813 : integer, intent(in) :: tend_level(:)
814 : ! Projected stresses.
815 : real(r8), intent(in) :: taucd(:,:,:)
816 : ! Components of momentum change sourced from the bottom.
817 : real(r8), intent(out) :: um_flux(:), vm_flux(:)
818 :
819 : integer :: i
820 :
821 : ! Tendency for U & V below source level.
822 0 : do i = 1, size(tend_level)
823 0 : um_flux(i) = taucd(i,tend_level(i)+1, east) + &
824 0 : taucd(i,tend_level(i)+1, west)
825 0 : vm_flux(i) = taucd(i,tend_level(i)+1,north) + &
826 0 : taucd(i,tend_level(i)+1,south)
827 : end do
828 :
829 0 : end subroutine momentum_flux
830 :
831 : !==========================================================================
832 :
833 : ! Subtracts a change in momentum in the gravity wave levels from wind
834 : ! tendencies in lower levels, ensuring momentum conservation.
835 0 : subroutine momentum_fixer(tend_level, p, um_flux, vm_flux, utgw, vtgw)
836 :
837 : ! Bottom stress level.
838 : integer, intent(in) :: tend_level(:)
839 : ! Pressure coordinates.
840 : type(Coords1D), intent(in) :: p
841 : ! Components of momentum change sourced from the bottom.
842 : real(r8), intent(in) :: um_flux(:), vm_flux(:)
843 : ! Wind tendencies.
844 : real(r8), intent(inout) :: utgw(:,:), vtgw(:,:)
845 :
846 : ! Indices.
847 : integer :: i, k
848 : ! Reciprocal of total mass.
849 0 : real(r8) :: rdm(size(tend_level))
850 : ! Average changes in velocity from momentum change being spread over
851 : ! total mass.
852 0 : real(r8) :: du(size(tend_level)), dv(size(tend_level))
853 :
854 : ! Total mass from ground to source level: rho*dz = dp/gravit
855 0 : do i = 1, size(tend_level)
856 0 : rdm(i) = gravit/(p%ifc(i,pver+1)-p%ifc(i,tend_level(i)+1))
857 : end do
858 :
859 : ! Average velocity changes.
860 0 : du = -um_flux*rdm
861 0 : dv = -vm_flux*rdm
862 :
863 0 : do k = minval(tend_level)+1, pver
864 0 : where (k > tend_level)
865 0 : utgw(:,k) = utgw(:,k) + du
866 0 : vtgw(:,k) = vtgw(:,k) + dv
867 : end where
868 : end do
869 :
870 0 : end subroutine momentum_fixer
871 :
872 : !==========================================================================
873 :
874 : ! Calculate the change in total energy from tendencies up to this point.
875 1489176 : subroutine energy_change(dt, p, u, v, dudt, dvdt, dsdt, de)
876 :
877 : ! Time step.
878 : real(r8), intent(in) :: dt
879 : ! Pressure coordinates.
880 : type(Coords1D), intent(in) :: p
881 : ! Winds at start of time step.
882 : real(r8), intent(in) :: u(:,:), v(:,:)
883 : ! Wind tendencies.
884 : real(r8), intent(in) :: dudt(:,:), dvdt(:,:)
885 : ! Heating tendency.
886 : real(r8), intent(in) :: dsdt(:,:)
887 : ! Change in energy.
888 : real(r8), intent(out) :: de(:)
889 :
890 : ! Level index.
891 : integer :: k
892 :
893 : ! Net gain/loss of total energy in the column.
894 24865776 : de = 0.0_r8
895 40207752 : do k = 1, pver
896 38718576 : de = de + p%del(:,k)/gravit * (dsdt(:,k) + &
897 0 : dudt(:,k)*(u(:,k)+dudt(:,k)*0.5_r8*dt) + &
898 686717928 : dvdt(:,k)*(v(:,k)+dvdt(:,k)*0.5_r8*dt) )
899 : end do
900 :
901 1489176 : end subroutine energy_change
902 :
903 : !==========================================================================
904 :
905 : ! Subtract change in energy from the heating tendency in the levels below
906 : ! the gravity wave region.
907 0 : subroutine energy_fixer(tend_level, p, de, ttgw)
908 :
909 : ! Bottom stress level.
910 : integer, intent(in) :: tend_level(:)
911 : ! Pressure coordinates.
912 : type(Coords1D), intent(in) :: p
913 : ! Change in energy.
914 : real(r8), intent(in) :: de(:)
915 : ! Heating tendency.
916 : real(r8), intent(inout) :: ttgw(:,:)
917 :
918 : ! Column/level indices.
919 : integer :: i, k
920 : ! Energy change to apply divided by all the mass it is spread across.
921 0 : real(r8) :: de_dm(size(tend_level))
922 :
923 0 : do i = 1, size(tend_level)
924 0 : de_dm(i) = -de(i)*gravit/(p%ifc(i,pver+1)-p%ifc(i,tend_level(i)+1))
925 : end do
926 :
927 : ! Subtract net gain/loss of total energy below tend_level.
928 0 : do k = minval(tend_level)+1, pver
929 0 : where (k > tend_level)
930 0 : ttgw(:,k) = ttgw(:,k) + de_dm
931 : end where
932 : end do
933 :
934 0 : end subroutine energy_fixer
935 :
936 : !==========================================================================
937 :
938 : ! Calculates absolute value of the local Coriolis frequency divided by the
939 : ! spatial frequency kwv, which gives a characteristic speed in m/s.
940 0 : function coriolis_speed(band, lat)
941 : ! Inertial gravity wave lengths.
942 : type(GWBand), intent(in) :: band
943 : ! Latitude in radians.
944 : real(r8), intent(in) :: lat(:)
945 :
946 : real(r8) :: coriolis_speed(size(lat))
947 :
948 : ! 24*3600 = 86400 seconds in a day.
949 : real(r8), parameter :: omega_earth = 2._r8*pi/86400._r8
950 :
951 0 : coriolis_speed = abs(sin(lat) * 2._r8 * omega_earth / band%kwv)
952 :
953 : end function coriolis_speed
954 :
955 : !==========================================================================
956 :
957 0 : subroutine adjust_inertial(band, tend_level, &
958 0 : u_coriolis, c, ubi, tau, ro_adjust)
959 : ! Inertial gravity wave lengths.
960 : type(GWBand), intent(in) :: band
961 : ! Levels above which tau is calculated.
962 : integer, intent(in) :: tend_level(:)
963 : ! Absolute value of the Coriolis frequency for each column,
964 : ! divided by kwv [m/s].
965 : real(r8), intent(in) :: u_coriolis(:)
966 : ! Wave propagation speed.
967 : real(r8), intent(in) :: c(:,-band%ngwv:)
968 : ! Wind speed in the direction of wave propagation.
969 : real(r8), intent(in) :: ubi(:,:)
970 :
971 : ! Tau will be adjusted by blocking wave propagation through cells where
972 : ! the Coriolis effect prevents it.
973 : real(r8), intent(inout) :: tau(:,-band%ngwv:,:)
974 : ! Dimensionless Coriolis term used to reduce gravity wave strength.
975 : ! Equal to max(0, 1 - (1/ro)^2), where ro is the Rossby number of the
976 : ! wind with respect to inertial waves.
977 : real(r8), intent(out) :: ro_adjust(:,-band%ngwv:,:)
978 :
979 : ! Column/level/wavenumber indices.
980 : integer :: i, k, l
981 :
982 : ! For each column and wavenumber, are we clear of levels that block
983 : ! upward propagation?
984 0 : logical :: unblocked_mask(size(tend_level),-band%ngwv:band%ngwv)
985 :
986 0 : unblocked_mask = .true.
987 0 : ro_adjust = 0._r8
988 :
989 : ! Iterate from the bottom up, through every interface level where tau is
990 : ! set.
991 0 : do k = maxval(tend_level)+1, ktop, -1
992 0 : do l = -band%ngwv, band%ngwv
993 0 : do i = 1, size(tend_level)
994 : ! Only operate on valid levels for this column.
995 0 : if (k <= tend_level(i) + 1) then
996 : ! Block waves if Coriolis is too strong.
997 : ! By setting the mask in this way, we avoid division by zero.
998 0 : unblocked_mask(i,l) = unblocked_mask(i,l) .and. &
999 0 : (abs(ubi(i,k) - c(i,l)) > u_coriolis(i))
1000 0 : if (unblocked_mask(i,l)) then
1001 0 : ro_adjust(i,l,k) = &
1002 0 : 1._r8 - (u_coriolis(i)/(ubi(i,k)-c(i,l)))**2
1003 : else
1004 0 : tau(i,l,k) = 0._r8
1005 : end if
1006 : end if
1007 : end do
1008 : end do
1009 : end do
1010 :
1011 0 : end subroutine adjust_inertial
1012 :
1013 0 : end module gw_common
|
