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