Line data Source code
1 : module gw_rdg
2 :
3 : !
4 : ! This module handles gravity waves from orographic sources, and was
5 : ! extracted from gw_drag in May 2013.
6 : !
7 : use shr_const_mod, only: pii => shr_const_pi
8 : use gw_utils, only: r8
9 : use gw_common, only: pver
10 : use coords_1d, only: Coords1D
11 : use spmd_utils,only: masterproc
12 : use cam_abortutils, only: endrun
13 :
14 :
15 : implicit none
16 : private
17 : save
18 :
19 : ! Public interface
20 : public :: gw_rdg_readnl
21 : public :: gw_rdg_src
22 : public :: gw_rdg_resid_src
23 : public :: gw_rdg_belowpeak
24 : public :: gw_rdg_break_trap
25 : public :: gw_rdg_do_vdiff
26 :
27 : ! Tunable Parameters
28 : !--------------------
29 : logical :: do_divstream
30 :
31 : !===========================================
32 : ! Parameters for DS2017 (do_divstream=.T.)
33 : !===========================================
34 : ! Amplification factor - 1.0 for
35 : ! high-drag/windstorm regime
36 : real(r8), protected :: C_BetaMax_DS
37 :
38 : ! Max Ratio Fr2:Fr1 - 1.0
39 : real(r8), protected :: C_GammaMax
40 :
41 : ! Normalized limits for Fr2(Frx) function
42 : real(r8), protected :: Frx0
43 : real(r8), protected :: Frx1
44 :
45 :
46 : !===========================================
47 : ! Parameters for SM2000
48 : !===========================================
49 : ! Amplification factor - 1.0 for
50 : ! high-drag/windstorm regime
51 : real(r8), protected :: C_BetaMax_SM
52 :
53 :
54 :
55 : ! NOTE: Critical inverse Froude number Fr_c is
56 : ! 1./(SQRT(2.)~0.707 in SM2000
57 : ! (should be <= 1)
58 : real(r8), protected :: Fr_c
59 :
60 : logical, protected :: gw_rdg_do_vdiff=.true.
61 :
62 : logical :: do_smooth_regimes
63 : logical :: do_adjust_tauoro
64 : logical :: do_backward_compat
65 :
66 :
67 : ! Limiters (min/max values)
68 : ! min surface displacement height for orographic waves (m)
69 : real(r8), protected :: orohmin
70 : ! min wind speed for orographic waves
71 : real(r8), protected :: orovmin
72 : ! min stratification allowing wave behavior
73 : real(r8), protected :: orostratmin
74 : ! min stratification allowing wave behavior
75 : real(r8), protected :: orom2min
76 :
77 : !==========================================================================
78 : contains
79 : !==========================================================================
80 :
81 2304 : subroutine gw_rdg_readnl(nlfile)
82 : use namelist_utils, only: find_group_name
83 : use units, only: getunit, freeunit
84 : use spmd_utils, only: mpicom, mstrid=>masterprocid, mpi_real8, mpi_logical
85 :
86 : ! File containing namelist input.
87 : character(len=*), intent(in) :: nlfile
88 :
89 : ! Local variables
90 : integer :: unitn, ierr
91 : character(len=*), parameter :: sub = 'gw_rdg_readnl'
92 :
93 : logical :: gw_rdg_do_divstream, gw_rdg_do_smooth_regimes, gw_rdg_do_adjust_tauoro, &
94 : gw_rdg_do_backward_compat
95 :
96 :
97 : real(r8) :: gw_rdg_C_BetaMax_DS, gw_rdg_C_GammaMax, &
98 : gw_rdg_Frx0, gw_rdg_Frx1, gw_rdg_C_BetaMax_SM, gw_rdg_Fr_c, &
99 : gw_rdg_orohmin, gw_rdg_orovmin, gw_rdg_orostratmin, gw_rdg_orom2min
100 :
101 : namelist /gw_rdg_nl/ gw_rdg_do_divstream, gw_rdg_C_BetaMax_DS, gw_rdg_C_GammaMax, &
102 : gw_rdg_Frx0, gw_rdg_Frx1, gw_rdg_C_BetaMax_SM, gw_rdg_Fr_c, &
103 : gw_rdg_do_smooth_regimes, gw_rdg_do_adjust_tauoro, &
104 : gw_rdg_do_backward_compat, gw_rdg_orohmin, gw_rdg_orovmin, &
105 : gw_rdg_orostratmin, gw_rdg_orom2min, gw_rdg_do_vdiff
106 :
107 : !----------------------------------------------------------------------
108 :
109 2304 : if (masterproc) then
110 3 : unitn = getunit()
111 3 : open( unitn, file=trim(nlfile), status='old' )
112 3 : call find_group_name(unitn, 'gw_rdg_nl', status=ierr)
113 3 : if (ierr == 0) then
114 3 : read(unitn, gw_rdg_nl, iostat=ierr)
115 3 : if (ierr /= 0) then
116 0 : call endrun(sub // ':: ERROR reading namelist')
117 : end if
118 : end if
119 3 : close(unitn)
120 3 : call freeunit(unitn)
121 :
122 : ! Set the local variables
123 3 : do_divstream = gw_rdg_do_divstream
124 3 : C_BetaMax_DS = gw_rdg_C_BetaMax_DS
125 3 : C_GammaMax = gw_rdg_C_GammaMax
126 3 : Frx0 = gw_rdg_Frx0
127 3 : Frx1 = gw_rdg_Frx1
128 3 : C_BetaMax_SM = gw_rdg_C_BetaMax_SM
129 3 : Fr_c = gw_rdg_Fr_c
130 3 : do_smooth_regimes = gw_rdg_do_smooth_regimes
131 3 : do_adjust_tauoro = gw_rdg_do_adjust_tauoro
132 3 : do_backward_compat = gw_rdg_do_backward_compat
133 3 : orohmin = gw_rdg_orohmin
134 3 : orovmin = gw_rdg_orovmin
135 3 : orostratmin = gw_rdg_orostratmin
136 3 : orom2min = gw_rdg_orom2min
137 : end if
138 :
139 : ! Broadcast the local variables
140 :
141 2304 : call mpi_bcast(do_divstream, 1, mpi_logical, mstrid, mpicom, ierr)
142 2304 : if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: do_divstream")
143 2304 : call mpi_bcast(do_smooth_regimes, 1, mpi_logical, mstrid, mpicom, ierr)
144 2304 : if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: do_smooth_regimes")
145 2304 : call mpi_bcast(do_adjust_tauoro, 1, mpi_logical, mstrid, mpicom, ierr)
146 2304 : if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: do_adjust_tauoro")
147 2304 : call mpi_bcast(do_backward_compat, 1, mpi_logical, mstrid, mpicom, ierr)
148 2304 : if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: do_backward_compat")
149 :
150 2304 : call mpi_bcast(C_BetaMax_DS, 1, mpi_real8, mstrid, mpicom, ierr)
151 2304 : if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: C_BetaMax_DS")
152 2304 : call mpi_bcast(C_GammaMax, 1, mpi_real8, mstrid, mpicom, ierr)
153 2304 : if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: C_GammaMax")
154 2304 : call mpi_bcast(Frx0, 1, mpi_real8, mstrid, mpicom, ierr)
155 2304 : if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: Frx0")
156 2304 : call mpi_bcast(Frx1, 1, mpi_real8, mstrid, mpicom, ierr)
157 2304 : if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: Frx1")
158 2304 : call mpi_bcast(C_BetaMax_SM, 1, mpi_real8, mstrid, mpicom, ierr)
159 2304 : if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: C_BetaMax_SM")
160 2304 : call mpi_bcast(Fr_c, 1, mpi_real8, mstrid, mpicom, ierr)
161 2304 : if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: Fr_c")
162 2304 : call mpi_bcast(orohmin, 1, mpi_real8, mstrid, mpicom, ierr)
163 2304 : if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: orohmin")
164 2304 : call mpi_bcast(orovmin, 1, mpi_real8, mstrid, mpicom, ierr)
165 2304 : if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: orovmin")
166 2304 : call mpi_bcast(orostratmin, 1, mpi_real8, mstrid, mpicom, ierr)
167 2304 : if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: orostratmin")
168 2304 : call mpi_bcast(orom2min, 1, mpi_real8, mstrid, mpicom, ierr)
169 2304 : if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: orom2min")
170 :
171 2304 : call mpi_bcast(gw_rdg_do_vdiff, 1, mpi_logical, mstrid, mpicom, ierr)
172 2304 : if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: gw_rdg_do_vdiff")
173 :
174 2304 : if (Fr_c > 1.0_r8) call endrun(sub//": FATAL: Fr_c must be <= 1")
175 :
176 2304 : end subroutine gw_rdg_readnl
177 :
178 :
179 : !==========================================================================
180 0 : subroutine gw_rdg_resid_src(ncol, band, p, &
181 0 : u, v, t, mxdis, kwvrdg, zi, nm, &
182 0 : src_level, tend_level, tau, ubm, ubi, xv, yv, &
183 0 : ubmsrc, usrc, vsrc, nsrc, rsrc, m2src, c, tauoro )
184 : use gw_common, only: rair, GWBand
185 : use gw_utils, only: dot_2d, midpoint_interp, get_unit_vector
186 : !-----------------------------------------------------------------------
187 : ! Orographic source for multiple gravity wave drag parameterization.
188 : !
189 : ! The stress is returned for a single wave with c=0, over orography.
190 : ! For points where the orographic variance is small (including ocean),
191 : ! the returned stress is zero.
192 : !------------------------------Arguments--------------------------------
193 : ! Column dimension.
194 : integer, intent(in) :: ncol
195 :
196 : ! Band to emit orographic waves in.
197 : ! Regardless, we will only ever emit into l = 0.
198 : type(GWBand), intent(in) :: band
199 : ! Pressure coordinates.
200 : type(Coords1D), intent(in) :: p
201 :
202 :
203 : ! Midpoint zonal/meridional winds. ( m s-1)
204 : real(r8), intent(in) :: u(ncol,pver), v(ncol,pver)
205 : ! Midpoint temperatures. (K)
206 : real(r8), intent(in) :: t(ncol,pver)
207 : ! Height estimate for ridge (m) [anisotropic orography].
208 : real(r8), intent(in) :: mxdis(ncol)
209 : ! horiz wavenumber [anisotropic orography].
210 : real(r8), intent(in) :: kwvrdg(ncol)
211 : ! Interface altitudes above ground (m).
212 : real(r8), intent(in) :: zi(ncol,pver+1)
213 : ! Midpoint Brunt-Vaisalla frequencies (s-1).
214 : real(r8), intent(in) :: nm(ncol,pver)
215 :
216 : ! Indices of top gravity wave source level and lowest level where wind
217 : ! tendencies are allowed.
218 : integer, intent(out) :: src_level(ncol)
219 : integer, intent(out) :: tend_level(ncol)
220 :
221 : ! Averages over source region.
222 : real(r8), intent(out) :: nsrc(ncol) ! B-V frequency.
223 : real(r8), intent(out) :: rsrc(ncol) ! Density.
224 : real(r8), intent(out) :: usrc(ncol) ! Zonal wind.
225 : real(r8), intent(out) :: vsrc(ncol) ! Meridional wind.
226 : real(r8), intent(out) :: ubmsrc(ncol) ! On-obstacle wind.
227 : ! normalized wavenumber
228 : real(r8), intent(out) :: m2src(ncol)
229 :
230 :
231 : ! Wave Reynolds stress.
232 : real(r8), intent(out) :: tau(ncol,-band%ngwv:band%ngwv,pver+1)
233 : ! Projection of wind at midpoints and interfaces.
234 : real(r8), intent(out) :: ubm(ncol,pver), ubi(ncol,pver+1)
235 : ! Unit vectors of source wind (zonal and meridional components).
236 : real(r8), intent(out) :: xv(ncol), yv(ncol)
237 : ! Phase speeds.
238 : real(r8), intent(out) :: c(ncol,-band%ngwv:band%ngwv)
239 : ! source level mom. flux
240 : real(r8), intent(out) :: tauoro(ncol)
241 :
242 : !---------------------------Local Storage-------------------------------
243 : ! Column and level indices.
244 : integer :: i, k
245 :
246 : ! Surface streamline displacement height (2*sgh).
247 0 : real(r8) :: hdsp(ncol)
248 :
249 : ! Difference in interface pressure across source region.
250 0 : real(r8) :: dpsrc(ncol)
251 : ! Thickness of downslope wind region.
252 : real(r8) :: ddw(ncol)
253 : ! Thickness of linear wave region.
254 : real(r8) :: dwv(ncol)
255 : ! Wind speed in source region.
256 0 : real(r8) :: wmsrc(ncol)
257 :
258 : real(r8) :: ragl(ncol)
259 : real(r8) :: Fcrit_res,sghmax
260 :
261 : !--------------------------------------------------------------------------
262 : ! Check that ngwav is equal to zero, otherwise end the job
263 : !--------------------------------------------------------------------------
264 0 : if (band%ngwv /= 0) call endrun(' gw_rdg_src :: ERROR - band%ngwv must be zero and it is not')
265 :
266 : !--------------------------------------------------------------------------
267 : ! Average the basic state variables for the wave source over the depth of
268 : ! the orographic standard deviation. Here we assume that the appropiate
269 : ! values of wind, stability, etc. for determining the wave source are
270 : ! averages over the depth of the atmosphere penterated by the typical
271 : ! mountain.
272 : ! Reduces to the bottom midpoint values when mxdis=0, such as over ocean.
273 : !--------------------------------------------------------------------------
274 :
275 0 : Fcrit_res = 1.0_r8
276 0 : hdsp = mxdis ! no longer multipied by 2
277 0 : where(hdsp < 10._r8)
278 : hdsp = 0._r8
279 : end where
280 :
281 0 : src_level = pver+1
282 :
283 0 : tau(:,0,:) = 0.0_r8
284 :
285 : ! Find depth of "source layer" for mountain waves
286 : ! i.e., between ground and mountain top
287 0 : do k = pver, 1, -1
288 0 : do i = 1, ncol
289 : ! Need to have h >= z(k+1) here or code will bomb when h=0.
290 0 : if ( (hdsp(i) >= zi(i,k+1)) .and. (hdsp(i) < zi(i,k)) ) then
291 0 : src_level(i) = k
292 : end if
293 : end do
294 : end do
295 :
296 0 : rsrc = 0._r8
297 0 : usrc = 0._r8
298 0 : vsrc = 0._r8
299 0 : nsrc = 0._r8
300 0 : do i = 1, ncol
301 0 : do k = pver, src_level(i), -1
302 0 : rsrc(i) = rsrc(i) + p%mid(i,k) / (rair*t(i,k))* p%del(i,k)
303 0 : usrc(i) = usrc(i) + u(i,k) * p%del(i,k)
304 0 : vsrc(i) = vsrc(i) + v(i,k) * p%del(i,k)
305 0 : nsrc(i) = nsrc(i) + nm(i,k)* p%del(i,k)
306 : end do
307 : end do
308 :
309 :
310 0 : do i = 1, ncol
311 0 : dpsrc(i) = p%ifc(i,pver+1) - p%ifc(i,src_level(i))
312 : end do
313 :
314 0 : rsrc = rsrc / dpsrc
315 0 : usrc = usrc / dpsrc
316 0 : vsrc = vsrc / dpsrc
317 0 : nsrc = nsrc / dpsrc
318 :
319 : ! Get the unit vector components and magnitude at the surface.
320 0 : call get_unit_vector(usrc, vsrc, xv, yv, wmsrc )
321 :
322 0 : ubmsrc = wmsrc
323 :
324 : ! Project the local wind at midpoints onto the source wind.
325 0 : do k = 1, pver
326 0 : ubm(:,k) = dot_2d(u(:,k), v(:,k), xv, yv)
327 : end do
328 :
329 : ! Compute the interface wind projection by averaging the midpoint winds.
330 : ! Use the top level wind at the top interface.
331 0 : ubi(:,1) = ubm(:,1)
332 :
333 0 : ubi(:,2:pver) = midpoint_interp(ubm)
334 :
335 : ! The minimum stratification allowing GW behavior
336 : ! should really depend on horizontal scale since
337 : !
338 : ! m^2 ~ (N/U)^2 - k^2
339 : !
340 :
341 0 : m2src = ( (nsrc/(ubmsrc+0.01_r8))**2 - kwvrdg**2 ) /((nsrc/(ubmsrc+0.01_r8))**2)
342 :
343 : ! Compute the interface wind projection by averaging the midpoint winds.
344 : ! Use the top level wind at the top interface.
345 0 : ubi(:,1) = ubm(:,1)
346 0 : ubi(:,2:pver) = midpoint_interp(ubm)
347 0 : ubi(:,pver+1) = ubm(:,pver)
348 :
349 :
350 :
351 : ! Determine the orographic c=0 source term following McFarlane (1987).
352 : ! (DOI: https://doi.org/10.1175/1520-0469(1987)044<1775:TEOOEG>2.0.CO;2)
353 : ! Set the source top interface index to pver, if the orographic term is
354 : ! zero.
355 0 : do i = 1, ncol
356 0 : if ( ( src_level(i) > 0 ) .and. ( m2src(i) > orom2min ) ) then
357 0 : sghmax = Fcrit_res * (ubmsrc(i) / nsrc(i))**2
358 : tauoro(i) = 0.5_r8 * kwvrdg(i) * min(hdsp(i)**2, sghmax) * &
359 0 : rsrc(i) * nsrc(i) * ubmsrc(i)
360 : else
361 0 : tauoro(i) = 0._r8
362 : end if
363 : end do
364 :
365 0 : do i = 1, ncol
366 0 : do k=src_level(i),pver+1
367 0 : tau(i,0,k) = tauoro(i)
368 : end do
369 : end do
370 :
371 :
372 : ! Allow wind tendencies all the way to the model bottom.
373 0 : tend_level = pver
374 :
375 : ! No spectrum; phase speed is just 0.
376 0 : c = 0._r8
377 :
378 0 : end subroutine gw_rdg_resid_src
379 :
380 :
381 : !==========================================================================
382 :
383 897840 : subroutine gw_rdg_src(ncol, band, p, &
384 897840 : u, v, t, mxdis, angxy, anixy, kwvrdg, iso, zi, nm, &
385 897840 : src_level, tend_level, bwv_level ,tlb_level , tau, ubm, ubi, xv, yv, &
386 897840 : ubmsrc, usrc, vsrc, nsrc, rsrc, m2src, tlb, bwv, Fr1, Fr2, Frx, c)
387 : use gw_common, only: rair, GWBand
388 : use gw_utils, only: dot_2d, midpoint_interp
389 : !-----------------------------------------------------------------------
390 : ! Orographic source for multiple gravity wave drag parameterization.
391 : !
392 : ! The stress is returned for a single wave with c=0, over orography.
393 : ! For points where the orographic variance is small (including ocean),
394 : ! the returned stress is zero.
395 : !------------------------------Arguments--------------------------------
396 : ! Column dimension.
397 : integer, intent(in) :: ncol
398 :
399 : ! Band to emit orographic waves in.
400 : ! Regardless, we will only ever emit into l = 0.
401 : type(GWBand), intent(in) :: band
402 : ! Pressure coordinates.
403 : type(Coords1D), intent(in) :: p
404 :
405 :
406 : ! Midpoint zonal/meridional winds. ( m s-1)
407 : real(r8), intent(in) :: u(ncol,pver), v(ncol,pver)
408 : ! Midpoint temperatures. (K)
409 : real(r8), intent(in) :: t(ncol,pver)
410 : ! Height estimate for ridge (m) [anisotropic orography].
411 : real(r8), intent(in) :: mxdis(ncol)
412 : ! Angle of ridge axis w/resp to north (degrees) [anisotropic orography].
413 : real(r8), intent(in) :: angxy(ncol)
414 : ! Anisotropy parameter [anisotropic orography].
415 : real(r8), intent(in) :: anixy(ncol)
416 : ! horiz wavenumber [anisotropic orography].
417 : real(r8), intent(in) :: kwvrdg(ncol)
418 : ! Isotropic source flag [anisotropic orography].
419 : integer, intent(in) :: iso(ncol)
420 : ! Interface altitudes above ground (m).
421 : real(r8), intent(in) :: zi(ncol,pver+1)
422 : ! Midpoint Brunt-Vaisalla frequencies (s-1).
423 : real(r8), intent(in) :: nm(ncol,pver)
424 :
425 : ! Indices of top gravity wave source level and lowest level where wind
426 : ! tendencies are allowed.
427 : integer, intent(out) :: src_level(ncol)
428 : integer, intent(out) :: tend_level(ncol)
429 : integer, intent(out) :: bwv_level(ncol),tlb_level(ncol)
430 :
431 : ! Averages over source region.
432 : real(r8), intent(out) :: nsrc(ncol) ! B-V frequency.
433 : real(r8), intent(out) :: rsrc(ncol) ! Density.
434 : real(r8), intent(out) :: usrc(ncol) ! Zonal wind.
435 : real(r8), intent(out) :: vsrc(ncol) ! Meridional wind.
436 : real(r8), intent(out) :: ubmsrc(ncol) ! On-ridge wind.
437 : ! Top of low-level flow layer.
438 : real(r8), intent(out) :: tlb(ncol)
439 : ! Bottom of linear wave region.
440 : real(r8), intent(out) :: bwv(ncol)
441 : ! normalized wavenumber
442 : real(r8), intent(out) :: m2src(ncol)
443 :
444 :
445 : ! Wave Reynolds stress.
446 : real(r8), intent(out) :: tau(ncol,-band%ngwv:band%ngwv,pver+1)
447 : ! Projection of wind at midpoints and interfaces.
448 : real(r8), intent(out) :: ubm(ncol,pver), ubi(ncol,pver+1)
449 : ! Unit vectors of source wind (zonal and meridional components).
450 : real(r8), intent(out) :: xv(ncol), yv(ncol)
451 : ! Phase speeds.
452 : real(r8), intent(out) :: c(ncol,-band%ngwv:band%ngwv)
453 : ! Froude numbers for flow/drag regimes
454 : real(r8), intent(out) :: Fr1(ncol), Fr2(ncol), Frx(ncol)
455 :
456 : !---------------------------Local Storage-------------------------------
457 : ! Column and level indices.
458 : integer :: i, k
459 :
460 : ! Surface streamline displacement height (2*sgh).
461 1795680 : real(r8) :: hdsp(ncol)
462 :
463 : ! Difference in interface pressure across source region.
464 1795680 : real(r8) :: dpsrc(ncol)
465 : ! Thickness of downslope wind region.
466 1795680 : real(r8) :: ddw(ncol)
467 : ! Thickness of linear wave region.
468 1795680 : real(r8) :: dwv(ncol)
469 : ! Wind speed in source region.
470 1795680 : real(r8) :: wmsrc(ncol)
471 :
472 1795680 : real(r8) :: ragl(ncol)
473 :
474 : !--------------------------------------------------------------------------
475 : ! Check that ngwav is equal to zero, otherwise end the job
476 : !--------------------------------------------------------------------------
477 897840 : if (band%ngwv /= 0) call endrun(' gw_rdg_src :: ERROR - band%ngwv must be zero and it is not')
478 :
479 : !--------------------------------------------------------------------------
480 : ! Average the basic state variables for the wave source over the depth of
481 : ! the orographic standard deviation. Here we assume that the appropiate
482 : ! values of wind, stability, etc. for determining the wave source are
483 : ! averages over the depth of the atmosphere penterated by the typical
484 : ! mountain.
485 : ! Reduces to the bottom midpoint values when mxdis=0, such as over ocean.
486 : !--------------------------------------------------------------------------
487 :
488 14991840 : hdsp = mxdis ! no longer multipied by 2
489 14991840 : src_level = pver+1
490 14991840 : bwv_level = -1
491 14991840 : tlb_level = -1
492 :
493 1410130800 : tau(:,0,:) = 0.0_r8
494 :
495 : ! Find depth of "source layer" for mountain waves
496 : ! i.e., between ground and mountain top
497 84396960 : do k = pver, 1, -1
498 1395138960 : do i = 1, ncol
499 : ! Need to have h >= z(k+1) here or code will bomb when h=0.
500 1394241120 : if ( (hdsp(i) >= zi(i,k+1)) .and. (hdsp(i) < zi(i,k)) ) then
501 14094000 : src_level(i) = k
502 : end if
503 : end do
504 : end do
505 :
506 14991840 : rsrc = 0._r8
507 14991840 : usrc = 0._r8
508 14991840 : vsrc = 0._r8
509 14991840 : nsrc = 0._r8
510 14991840 : do i = 1, ncol
511 33126218 : do k = pver, src_level(i), -1
512 18134378 : rsrc(i) = rsrc(i) + p%mid(i,k) / (rair*t(i,k))* p%del(i,k)
513 18134378 : usrc(i) = usrc(i) + u(i,k) * p%del(i,k)
514 18134378 : vsrc(i) = vsrc(i) + v(i,k) * p%del(i,k)
515 32228378 : nsrc(i) = nsrc(i) + nm(i,k)* p%del(i,k)
516 : end do
517 : end do
518 :
519 :
520 14991840 : do i = 1, ncol
521 14991840 : dpsrc(i) = p%ifc(i,pver+1) - p%ifc(i,src_level(i))
522 : end do
523 :
524 14991840 : rsrc = rsrc / dpsrc
525 14991840 : usrc = usrc / dpsrc
526 14991840 : vsrc = vsrc / dpsrc
527 14991840 : nsrc = nsrc / dpsrc
528 :
529 14991840 : wmsrc = sqrt( usrc**2 + vsrc**2 )
530 :
531 :
532 : ! Get the unit vector components
533 : ! Want agl=0 with U>0 to give xv=1
534 :
535 14991840 : ragl = angxy * pii/180._r8
536 :
537 : ! protect from wierd "bad" angles
538 : ! that may occur if hdsp is zero
539 14991840 : where( hdsp <= orohmin )
540 : ragl = 0._r8
541 : end where
542 :
543 14991840 : yv =-sin( ragl )
544 14991840 : xv = cos( ragl )
545 :
546 :
547 : ! Kluge in possible "isotropic" obstacle.
548 43179840 : where( ( iso == 1 ) .and. (wmsrc > orovmin) )
549 : xv = usrc/wmsrc
550 : yv = vsrc/wmsrc
551 : end where
552 :
553 :
554 : ! Project the local wind at midpoints into the on-ridge direction
555 84396960 : do k = 1, pver
556 84396960 : ubm(:,k) = dot_2d(u(:,k), v(:,k), xv, yv)
557 : end do
558 897840 : ubmsrc = dot_2d(usrc , vsrc , xv, yv)
559 :
560 : ! Ensure on-ridge wind is positive at source level
561 84396960 : do k = 1, pver
562 1395138960 : ubm(:,k) = sign( ubmsrc*0._r8+1._r8 , ubmsrc ) * ubm(:,k)
563 : end do
564 :
565 : ! Sean says just use 1._r8 as
566 : ! first argument
567 14991840 : xv = sign( ubmsrc*0._r8+1._r8 , ubmsrc ) * xv
568 14991840 : yv = sign( ubmsrc*0._r8+1._r8 , ubmsrc ) * yv
569 :
570 : ! Now make ubmsrc positive and protect
571 : ! against zero
572 14991840 : ubmsrc = abs(ubmsrc)
573 14991840 : ubmsrc = max( 0.01_r8 , ubmsrc )
574 :
575 :
576 : ! The minimum stratification allowing GW behavior
577 : ! should really depend on horizontal scale since
578 : !
579 : ! m^2 ~ (N/U)^2 - k^2
580 : !
581 : ! Should also think about parameterizing
582 : ! trapped lee-waves.
583 :
584 :
585 : ! This needs to be made constistent with later
586 : ! treatment of nonhydrostatic effects.
587 14991840 : m2src = ( (nsrc/(ubmsrc+0.01_r8))**2 - kwvrdg**2 ) /((nsrc/(ubmsrc+0.01_r8))**2)
588 :
589 :
590 : !-------------------------------------------------------------
591 : ! Calculate provisional limits (in Z [m]) for 3 regimes. This
592 : ! will modified later if wave breaking or trapping are
593 : ! diagnosed
594 : !
595 : ! ^
596 : ! | *** linear propagation ***
597 : ! (H) -------- mountain top ------------- | *** or wave breaking ****
598 : ! | *** regimes *************
599 : ! (BWV)------ bottom of linear waves ---- |
600 : ! : |
601 : ! ******* |
602 : ! : |
603 : ! (TLB)--- top of flow diversion layer--- '
604 : ! :
605 : ! **** flow diversion *****
606 : ! :
607 : !============================================
608 :
609 : !============================================
610 : ! For Dividing streamline para (DS2017)
611 : !--------------------------------------------
612 : ! High-drag downslope wind regime exists
613 : ! between bottom of linear waves and top of
614 : ! flow diversion. Linear waves can only
615 : ! attain vertical displacment of f1*U/N. So,
616 : ! bottom of linear waves is given by
617 : !
618 : ! BWV = H - Fr1*U/N
619 : !
620 : ! Downslope wind layer begins at BWV and
621 : ! extends below it until some maximum high
622 : ! drag obstacle height Fr2*U/N is attained
623 : ! (where Fr2 >= f1). Below downslope wind
624 : ! there is flow diversion, so top of
625 : ! diversion layer (TLB) is equivalent to
626 : ! bottom of downslope wind layer and is;
627 : !
628 : ! TLB = H - Fr2*U/N
629 : !
630 : !-----------------------------------------
631 :
632 : ! Critical inverse Froude number
633 : !-----------------------------------------------
634 14991840 : Fr1(:) = Fr_c * 1.00_r8
635 14991840 : Frx(:) = hdsp(:)*nsrc(:)/abs( ubmsrc(:) ) / Fr_c
636 :
637 897840 : if ( do_divstream ) then
638 : !------------------------------------------------
639 : ! Calculate Fr2(Frx) for DS2017
640 : !------------------------------------------------
641 85461840 : where(Frx <= Frx0)
642 : Fr2(:) = Fr1(:) + Fr1(:)* C_GammaMax * anixy(:)
643 : elsewhere((Frx > Frx0).and.(Frx <= Frx1) )
644 : Fr2(:) = Fr1(:) + Fr1(:)* C_GammaMax * anixy(:) &
645 : * (Frx1 - Frx(:))/(Frx1-Frx0)
646 : elsewhere(Frx > Frx1)
647 : Fr2(:)=Fr1(:)
648 : endwhere
649 : else
650 : !------------------------------------------
651 : ! Regime distinctions entirely carried by
652 : ! amplification of taudsw (next subr)
653 : !------------------------------------------
654 0 : Fr2(:)=Fr1(:)
655 : end if
656 :
657 :
658 :
659 14991840 : where( m2src > orom2min )
660 : ddw = Fr2 * ( abs(ubmsrc) )/nsrc
661 : elsewhere
662 : ddw = 0._r8
663 : endwhere
664 :
665 :
666 : ! If TLB is less than zero then obstacle is not
667 : ! high enough to produce an low-level diversion layer
668 14991840 : tlb = mxdis - ddw
669 14991840 : where( tlb < 0._r8)
670 : tlb = 0._r8
671 : endwhere
672 43994160 : do k = pver, pver/2, -1
673 720506160 : do i = 1, ncol
674 719608320 : if ( (tlb(i) > zi(i,k+1)) .and. (tlb(i) <= zi(i,k)) ) then
675 363338 : tlb_level(i) = k
676 : end if
677 : end do
678 : end do
679 :
680 :
681 : ! Find *BOTTOM* of linear wave layer (BWV)
682 : !where ( nsrc > orostratmin )
683 14991840 : where( m2src > orom2min )
684 : dwv = Fr1 * ( abs(ubmsrc) )/nsrc
685 : elsewhere
686 : dwv = -9.999e9_r8 ! if weak strat - no waves
687 : endwhere
688 :
689 14991840 : bwv = mxdis - dwv
690 14991840 : where(( bwv < 0._r8) .or. (dwv < 0._r8) )
691 : bwv = 0._r8
692 : endwhere
693 84396960 : do k = pver,1, -1
694 1395138960 : do i = 1, ncol
695 1394241120 : if ( (bwv(i) > zi(i,k+1)) .and. (bwv(i) <= zi(i,k)) ) then
696 581501 : bwv_level(i) = k+1
697 : end if
698 : end do
699 : end do
700 :
701 :
702 :
703 : ! Compute the interface wind projection by averaging the midpoint winds.
704 : ! Use the top level wind at the top interface.
705 14991840 : ubi(:,1) = ubm(:,1)
706 897840 : ubi(:,2:pver) = midpoint_interp(ubm)
707 14991840 : ubi(:,pver+1) = ubm(:,pver)
708 :
709 : ! Allow wind tendencies all the way to the model bottom.
710 14991840 : tend_level = pver
711 :
712 : ! No spectrum; phase speed is just 0.
713 15889680 : c = 0._r8
714 :
715 43179840 : where( m2src < orom2min )
716 : tlb = mxdis
717 : tlb_level = src_level
718 : endwhere
719 :
720 :
721 897840 : end subroutine gw_rdg_src
722 :
723 :
724 : !==========================================================================
725 :
726 897840 : subroutine gw_rdg_belowpeak(ncol, band, rdg_cd_llb, &
727 897840 : t, mxdis, anixy, kwvrdg, zi, nm, ni, rhoi, &
728 897840 : src_level , tau, &
729 897840 : ubmsrc, nsrc, rsrc, m2src,tlb,bwv,Fr1,Fr2,Frx, &
730 897840 : tauoro,taudsw, hdspwv,hdspdw )
731 :
732 : use gw_common, only: GWBand
733 : !-----------------------------------------------------------------------
734 : ! Orographic source for multiple gravity wave drag parameterization.
735 : !
736 : ! The stress is returned for a single wave with c=0, over orography.
737 : ! For points where the orographic variance is small (including ocean),
738 : ! the returned stress is zero.
739 : !------------------------------Arguments--------------------------------
740 : ! Column dimension.
741 : integer, intent(in) :: ncol
742 : ! Band to emit orographic waves in.
743 : ! Regardless, we will only ever emit into l = 0.
744 : type(GWBand), intent(in) :: band
745 : ! Drag coefficient for low-level flow
746 : real(r8), intent(in) :: rdg_cd_llb
747 :
748 :
749 : ! Midpoint temperatures. (K)
750 : real(r8), intent(in) :: t(ncol,pver)
751 : ! Height estimate for ridge (m) [anisotropic orography].
752 : real(r8), intent(in) :: mxdis(ncol)
753 : ! Anisotropy parameter [0-1] [anisotropic orography].
754 : real(r8), intent(in) :: anixy(ncol)
755 : ! Inverse cross-ridge lengthscale (m-1) [anisotropic orography].
756 : real(r8), intent(inout) :: kwvrdg(ncol)
757 : ! Interface altitudes above ground (m).
758 : real(r8), intent(in) :: zi(ncol,pver+1)
759 : ! Midpoint Brunt-Vaisalla frequencies (s-1).
760 : real(r8), intent(in) :: nm(ncol,pver)
761 : ! Interface Brunt-Vaisalla frequencies (s-1).
762 : real(r8), intent(in) :: ni(ncol,pver+1)
763 : ! Interface density (kg m-3).
764 : real(r8), intent(in) :: rhoi(ncol,pver+1)
765 :
766 : ! Indices of top gravity wave source level
767 : integer, intent(inout) :: src_level(ncol)
768 :
769 : ! Wave Reynolds stress.
770 : real(r8), intent(inout) :: tau(ncol,-band%ngwv:band%ngwv,pver+1)
771 : ! Top of low-level flow layer.
772 : real(r8), intent(in) :: tlb(ncol)
773 : ! Bottom of linear wave region.
774 : real(r8), intent(in) :: bwv(ncol)
775 : ! surface stress from linear waves.
776 : real(r8), intent(out) :: tauoro(ncol)
777 : ! surface stress for downslope wind regime.
778 : real(r8), intent(out) :: taudsw(ncol)
779 :
780 : ! Surface streamline displacement height for linear waves.
781 : real(r8), intent(out) :: hdspwv(ncol)
782 : ! Surface streamline displacement height for downslope wind regime.
783 : real(r8), intent(out) :: hdspdw(ncol)
784 :
785 :
786 :
787 : ! Froude numbers for flow/drag regimes
788 : real(r8), intent(in) :: Fr1(ncol), Fr2(ncol),Frx(ncol)
789 :
790 : ! Averages over source region.
791 : real(r8), intent(in) :: m2src(ncol) ! normalized non-hydro wavenumber
792 : real(r8), intent(in) :: nsrc(ncol) ! B-V frequency.
793 : real(r8), intent(in) :: rsrc(ncol) ! Density.
794 : real(r8), intent(in) :: ubmsrc(ncol) ! On-ridge wind.
795 :
796 :
797 : !logical, intent(in), optional :: forcetlb
798 :
799 : !---------------------------Local Storage-------------------------------
800 : ! Column and level indices.
801 : integer :: i, k
802 :
803 1795680 : real(r8) :: Coeff_LB(ncol),tausat,ubsrcx(ncol),dswamp
804 1795680 : real(r8) :: taulin(ncol),BetaMax
805 :
806 : ! ubsrcx introduced to account for situations with high shear, strong strat.
807 14991840 : do i = 1, ncol
808 14991840 : ubsrcx(i) = max( ubmsrc(i) , 0._r8 )
809 : end do
810 :
811 14991840 : do i = 1, ncol
812 14991840 : if ( m2src(i) > orom2min ) then
813 1020198 : hdspwv(i) = min( mxdis(i) , Fr1(i) * ubsrcx(i) / nsrc(i) )
814 : else
815 13073802 : hdspwv(i) = 0._r8
816 : end if
817 : end do
818 :
819 897840 : if (do_divstream) then
820 14991840 : do i = 1, ncol
821 14991840 : if ( m2src(i) > orom2min ) then
822 1020198 : hdspdw(i) = min( mxdis(i) , Fr2(i) * ubsrcx(i) / nsrc(i) )
823 : else
824 13073802 : hdspdw(i) = 0._r8
825 : end if
826 : end do
827 : else
828 0 : do i = 1, ncol
829 : ! Needed only to mark where a DSW occurs
830 0 : if ( m2src(i) > orom2min ) then
831 0 : hdspdw(i) = mxdis(i)
832 : else
833 0 : hdspdw(i) = 0._r8
834 : end if
835 : end do
836 : end if
837 :
838 : ! Calculate form drag coefficient ("CD")
839 : !--------------------------------------
840 14991840 : Coeff_LB = rdg_cd_llb*anixy
841 :
842 : ! Determine the orographic c=0 source term following McFarlane (1987).
843 : ! Set the source top interface index to pver, if the orographic term is
844 : ! zero.
845 : !
846 : ! This formula is basically from
847 : !
848 : ! tau(src) = rho * u' * w'
849 : ! where
850 : ! u' ~ N*h' and w' ~ U*h'/b (b="breite")
851 : !
852 : ! and 1/b has been replaced with k (kwvrdg)
853 : !
854 14991840 : do i = 1, ncol
855 14991840 : if ( ( src_level(i) > 0 ) .and. ( m2src(i) > orom2min ) ) then
856 : tauoro(i) = kwvrdg(i) * ( hdspwv(i)**2 ) * rsrc(i) * nsrc(i) &
857 1020198 : * ubsrcx(i)
858 : taudsw(i) = kwvrdg(i) * ( hdspdw(i)**2 ) * rsrc(i) * nsrc(i) &
859 1020198 : * ubsrcx(i)
860 : else
861 13073802 : tauoro(i) = 0._r8
862 13073802 : taudsw(i) = 0._r8
863 : end if
864 : end do
865 :
866 897840 : if (do_divstream) then
867 14991840 : do i = 1, ncol
868 14991840 : taulin(i) = 0._r8
869 : end do
870 : !---------------------------------------
871 : ! Need linear drag when divstream is not used is used
872 : !---------------------------------------
873 : else
874 0 : do i = 1, ncol
875 0 : if ( ( src_level(i) > 0 ) .and. ( m2src(i) > orom2min ) ) then
876 : taulin(i) = kwvrdg(i) * ( mxdis(i)**2 ) * rsrc(i) * nsrc(i) &
877 0 : * ubsrcx(i)
878 : else
879 0 : taulin(i) = 0._r8
880 : end if
881 : end do
882 : end if
883 :
884 897840 : if ( do_divstream ) then
885 : ! Amplify DSW between Frx=1. and Frx=Frx1
886 14991840 : do i = 1,ncol
887 14094000 : dswamp=0._r8
888 14094000 : BetaMax = C_BetaMax_DS * anixy(i)
889 14094000 : if ( (Frx(i)>1._r8).and.(Frx(i)<=Frx1)) then
890 288383 : dswamp = (Frx(i)-1._r8)*(Frx1-Frx(i))/(0.25_r8*(Frx1-1._r8)**2)
891 : end if
892 14991840 : taudsw(i) = (1._r8 + BetaMax*dswamp)*taudsw(i)
893 : end do
894 : else
895 : !-------------------
896 : ! Scinocca&McFarlane
897 : !--------------------
898 0 : do i = 1, ncol
899 0 : BetaMax = C_BetaMax_SM * anixy(i)
900 0 : if ( (Frx(i) >=1._r8) .and. (Frx(i) < 1.5_r8) ) then
901 0 : dswamp = 2._r8 * BetaMax * (Frx(i) -1._r8)
902 0 : else if ( ( Frx(i) >= 1.5_r8 ) .and. (Frx(i) < 3._r8 ) ) then
903 : dswamp = ( 1._r8 + BetaMax - (0.666_r8**2) ) * ( 0.666_r8*(3._r8 - Frx(i) ))**2 &
904 0 : + ( 1._r8 / Frx(i) )**2 -1._r8
905 : else
906 : dswamp = 0._r8
907 : end if
908 0 : if ( (Frx(i) >=1._r8) .and. (Frx(i) < 3._r8) ) then
909 0 : taudsw(i) = (1._r8 + dswamp )*taulin(i) - tauoro(i)
910 : else
911 0 : taudsw(i) = 0._r8
912 : endif
913 : ! This code defines "taudsw" as SUM of freely-propagating
914 : ! DSW enhancement. Different than in SM2000
915 0 : taudsw(i) = taudsw(i) + tauoro(i)
916 : end do
917 : !----------------------------------------------------
918 : end if
919 :
920 :
921 14991840 : do i = 1, ncol
922 14991840 : if ( m2src(i) > orom2min ) then
923 385634844 : where ( ( zi(i,:) < mxdis(i) ) .and. ( zi(i,:) >= bwv(i) ) )
924 1020198 : tau(i,0,:) = tauoro(i)
925 : else where ( ( zi(i,:) < bwv(i) ) .and. ( zi(i,:) >= tlb(i) ) )
926 : tau(i,0,:) = tauoro(i) +( taudsw(i)-tauoro(i) )* &
927 : ( bwv(i) - zi(i,:) ) / &
928 : ( bwv(i) - tlb(i) )
929 : endwhere
930 : ! low-level form drag on obstacle. Quantity kwvrdg (~1/b) appears for consistency
931 : ! with tauoro and taudsw forms. Should be weighted by L*b/A_g before applied to flow.
932 96918810 : where ( ( zi(i,:) < tlb(i) ) .and. ( zi(i,:) >= 0._r8 ) )
933 : tau(i,0,:) = taudsw(i) + &
934 : Coeff_LB(i) * kwvrdg(i) * rsrc(i) * 0.5_r8 * (ubsrcx(i)**2) * ( tlb(i) - zi(i,:) )
935 : endwhere
936 :
937 1020198 : if (do_smooth_regimes) then
938 : ! This blocks accounts for case where both mxdis and tlb fall
939 : ! between adjacent edges
940 0 : do k=1,pver
941 0 : if ( (zi(i,k) >= tlb(i)).and.(zi(i,k+1) < tlb(i)).and. &
942 0 : (zi(i,k) >= mxdis(i)).and.(zi(i,k+1) < mxdis(i)) ) then
943 0 : src_level(i) = src_level(i)-1
944 0 : tau(i,0,k) = tauoro(i)
945 : end if
946 : end do
947 : end if
948 :
949 : else !----------------------------------------------
950 : ! This block allows low-level dynamics to occur
951 : ! even if m2 is less than orom2min
952 1242011190 : where ( ( zi(i,:) < tlb(i) ) .and. ( zi(i,:) >= 0._r8 ) )
953 13073802 : tau(i,0,:) = taudsw(i) + &
954 : Coeff_LB(i) * kwvrdg(i) * rsrc(i) * 0.5_r8 * &
955 : (ubsrcx(i)**2) * ( tlb(i) - zi(i,:) )
956 : endwhere
957 : endif
958 : end do
959 :
960 : ! This may be redundant with newest version of gw_drag_prof.
961 : ! That code reaches down to level k=src_level+1. (jtb 1/5/16)
962 14991840 : do i = 1, ncol
963 14094000 : k=src_level(i)
964 14094000 : if ( ni(i,k) > orostratmin ) then
965 : tausat = (Fr_c**2) * kwvrdg(i) * rhoi(i,k) * ubsrcx(i)**3 / &
966 14094000 : (1._r8*ni(i,k))
967 : else
968 : tausat = 0._r8
969 : endif
970 14991840 : tau(i,0,src_level(i)) = min( tauoro(i), tausat )
971 : end do
972 :
973 :
974 :
975 : ! Final clean-up. Do nothing if obstacle less than orohmin
976 14991840 : do i = 1, ncol
977 14991840 : if ( mxdis(i) < orohmin ) then
978 1244056065 : tau(i,0,:) = 0._r8
979 13095327 : tauoro(i) = 0._r8
980 13095327 : taudsw(i) = 0._r8
981 : endif
982 : end do
983 :
984 : ! Disable vertical propagation if Scorer param is
985 : ! too small.
986 14991840 : do i = 1, ncol
987 14991840 : if ( m2src(i) <= orom2min ) then
988 13073802 : src_level(i)=1
989 : endif
990 : end do
991 :
992 :
993 :
994 897840 : end subroutine gw_rdg_belowpeak
995 :
996 : !==========================================================================
997 897840 : subroutine gw_rdg_break_trap(ncol, band, &
998 897840 : zi, nm, ni, ubm, ubi, rhoi, kwvrdg, bwv, tlb, wbr, &
999 897840 : src_level, tlb_level, &
1000 897840 : hdspwv, hdspdw, mxdis, &
1001 897840 : tauoro, taudsw, tau, &
1002 : ldo_trapped_waves, wdth_kwv_scale_in )
1003 : use gw_common, only: GWBand
1004 : !-----------------------------------------------------------------------
1005 : ! Parameterization of high-drag regimes and trapped lee-waves for CAM
1006 : !
1007 : !------------------------------Arguments--------------------------------
1008 : ! Column dimension.
1009 : integer, intent(in) :: ncol
1010 : ! Band to emit orographic waves in.
1011 : ! Regardless, we will only ever emit into l = 0.
1012 : type(GWBand), intent(in) :: band
1013 :
1014 :
1015 : ! Height estimate for ridge (m) [anisotropic orography].
1016 : !real(r8), intent(in) :: mxdis(ncol)
1017 : ! Horz wavenumber for ridge (1/m) [anisotropic orography].
1018 : real(r8), intent(in) :: kwvrdg(ncol)
1019 : ! Interface altitudes above ground (m).
1020 : real(r8), intent(in) :: zi(ncol,pver+1)
1021 : ! Midpoint Brunt-Vaisalla frequencies (s-1).
1022 : real(r8), intent(in) :: nm(ncol,pver)
1023 : ! Interface Brunt-Vaisalla frequencies (s-1).
1024 : real(r8), intent(in) :: ni(ncol,pver+1)
1025 :
1026 : ! Indices of gravity wave sources.
1027 : integer, intent(inout) :: src_level(ncol), tlb_level(ncol)
1028 :
1029 : ! Wave Reynolds stress.
1030 : real(r8), intent(inout) :: tau(ncol,-band%ngwv:band%ngwv,pver+1)
1031 : ! Wave Reynolds stresses at source.
1032 : real(r8), intent(in) :: taudsw(ncol)
1033 : real(r8), intent(inout) :: tauoro(ncol)
1034 : ! Projection of wind at midpoints and interfaces.
1035 : real(r8), intent(in) :: ubm(ncol,pver)
1036 : real(r8), intent(in) :: ubi(ncol,pver+1)
1037 : ! Interface density (kg m-3).
1038 : real(r8), intent(in) :: rhoi(ncol,pver+1)
1039 :
1040 : ! Top of low-level flow layer.
1041 : real(r8), intent(in) :: tlb(ncol)
1042 : ! Bottom of linear wave region.
1043 : real(r8), intent(in) :: bwv(ncol)
1044 :
1045 : ! Surface streamline displacement height for linear waves.
1046 : real(r8), intent(in) :: hdspwv(ncol)
1047 : ! Surface streamline displacement height for downslope wind regime.
1048 : real(r8), intent(in) :: hdspdw(ncol)
1049 : ! Ridge height.
1050 : real(r8), intent(in) :: mxdis(ncol)
1051 :
1052 :
1053 : ! Wave breaking level
1054 : real(r8), intent(out) :: wbr(ncol)
1055 :
1056 : logical, intent(in), optional :: ldo_trapped_waves
1057 : real(r8), intent(in), optional :: wdth_kwv_scale_in
1058 :
1059 : !---------------------------Local Storage-------------------------------
1060 : ! Column and level indices.
1061 : integer :: i, k, kp1, non_hydro
1062 1795680 : real(r8):: m2(ncol,pver),delz(ncol),tausat(ncol),trn(ncol)
1063 1795680 : real(r8):: wbrx(ncol)
1064 1795680 : real(r8):: phswkb(ncol,pver+1)
1065 : logical :: lldo_trapped_waves
1066 : real(r8):: wdth_kwv_scale
1067 : ! Indices of important levels.
1068 1795680 : integer :: trn_level(ncol)
1069 :
1070 897840 : if (present(ldo_trapped_waves)) then
1071 897840 : lldo_trapped_waves = ldo_trapped_waves
1072 897840 : if(lldo_trapped_waves) then
1073 : non_hydro = 1
1074 : else
1075 897840 : non_hydro = 0
1076 : endif
1077 : else
1078 : lldo_trapped_waves = .false.
1079 : non_hydro = 0
1080 : endif
1081 :
1082 897840 : if (present(wdth_kwv_scale_in)) then
1083 0 : wdth_kwv_scale = wdth_kwv_scale_in
1084 : else
1085 : wdth_kwv_scale = 1._r8
1086 : endif
1087 :
1088 : ! Calculate vertical wavenumber**2
1089 : !---------------------------------
1090 1395138960 : m2 = (nm / (abs(ubm)+.01_r8))**2
1091 84396960 : do k=pver,1,-1
1092 1394241120 : m2(:,k) = m2(:,k) - non_hydro*(wdth_kwv_scale*kwvrdg)**2
1093 : ! sweeping up, zero out m2 above first occurence
1094 : ! of m2(:,k)<=0
1095 83499120 : kp1=min( k+1, pver )
1096 1395138960 : where( (m2(:,k) <= 0.0_r8 ).or.(m2(:,kp1) <= 0.0_r8 ) )
1097 : m2(:,k) = 0._r8
1098 : endwhere
1099 : end do
1100 :
1101 : ! Take square root of m**2 and
1102 : ! do vertical integral to find
1103 : ! WKB phase.
1104 : !-----------------------------
1105 1395138960 : m2 = SQRT( m2 )
1106 1410130800 : phswkb(:,:)=0
1107 84396960 : do k=pver,1,-1
1108 4183621200 : where( zi(:,k) > tlb(:) )
1109 83499120 : delz(:) = min( zi(:,k)-zi(:,k+1) , zi(:,k)-tlb(:) )
1110 83499120 : phswkb(:,k) = phswkb(:,k+1) + m2(:,k)*delz(:)
1111 : endwhere
1112 : end do
1113 :
1114 : ! Identify top edge of layer in which phswkb reaches 3*pi/2
1115 : ! - approximately the "breaking level"
1116 : !----------------------------------------------------------
1117 14991840 : wbr(:)=0._r8
1118 14991840 : wbrx(:)=0._r8
1119 897840 : if (do_smooth_regimes) then
1120 0 : do k=pver,1,-1
1121 0 : where( (phswkb(:,k+1)<1.5_r8*pii).and.(phswkb(:,k)>=1.5_r8*pii) &
1122 0 : .and.(hdspdw(:)>hdspwv(:)) )
1123 : wbr(:) = zi(:,k)
1124 : ! Extrapolation to make regime
1125 : ! transitions smoother
1126 : wbrx(:) = zi(:,k) - ( phswkb(:,k) - 1.5_r8*pii ) &
1127 0 : / ( m2(:,k) + 1.e-6_r8 )
1128 : src_level(:) = k-1
1129 : endwhere
1130 : end do
1131 : else
1132 84396960 : do k=pver,1,-1
1133 166998240 : where( (phswkb(:,k+1)<1.5_r8*pii).and.(phswkb(:,k)>=1.5_r8*pii) &
1134 4183621200 : .and.(hdspdw(:)>hdspwv(:)) )
1135 : wbr(:) = zi(:,k)
1136 : src_level(:) = k
1137 : endwhere
1138 : end do
1139 : end if
1140 :
1141 : ! Adjust tauoro at new source levels if needed.
1142 : ! This is problematic if Fr_c<1.0. Not sure why.
1143 : !----------------------------------------------------------
1144 897840 : if (do_adjust_tauoro) then
1145 14991840 : do i = 1,ncol
1146 14991840 : if (wbr(i) > 0._r8 ) then
1147 288383 : tausat(i) = (Fr_c**2) * kwvrdg(i) * rhoi( i, src_level(i) ) &
1148 : * abs(ubi(i , src_level(i) ))**3 &
1149 288383 : / ni( i , src_level(i) )
1150 288383 : tauoro(i) = min( tauoro(i), tausat(i) )
1151 : end if
1152 : end do
1153 : end if
1154 :
1155 897840 : if (do_smooth_regimes) then
1156 0 : do i = 1, ncol
1157 0 : do k=1,pver+1
1158 0 : if ( ( zi(i,k) <= wbr(i) ) .and. ( zi(i,k) > tlb(i) ) ) then
1159 0 : tau(i,0,k) = tauoro(i) + (taudsw(i)-tauoro(i)) * &
1160 : ( wbrx(i) - zi(i,k) ) / &
1161 0 : ( wbrx(i) - tlb(i) )
1162 0 : tau(i,0,k) = max( tau(i,0,k), tauoro(i) )
1163 : endif
1164 : end do
1165 : end do
1166 : else
1167 : ! Following is for backwards B4B compatibility with earlier versions
1168 : ! ("N1" and "N5" -- Note: "N5" used do_backward_compat=.true.)
1169 897840 : if (.not.do_backward_compat) then
1170 14991840 : do i = 1, ncol
1171 1339827840 : do k=1,pver+1
1172 1338930000 : if ( ( zi(i,k) < wbr(i) ) .and. ( zi(i,k) >= tlb(i) ) ) then
1173 2907092 : tau(i,0,k) = tauoro(i) + (taudsw(i)-tauoro(i)) * &
1174 : ( wbr(i) - zi(i,k) ) / &
1175 2907092 : ( wbr(i) - tlb(i) )
1176 : endif
1177 : end do
1178 : end do
1179 : else
1180 0 : do i = 1, ncol
1181 0 : do k=1,pver+1
1182 0 : if ( ( zi(i,k) <= wbr(i) ) .and. ( zi(i,k) > tlb(i) ) ) then
1183 0 : tau(i,0,k) = tauoro(i) + (taudsw(i)-tauoro(i)) * &
1184 : ( wbr(i) - zi(i,k) ) / &
1185 0 : ( wbr(i) - tlb(i) )
1186 : endif
1187 : end do
1188 : end do
1189 : end if
1190 : end if
1191 :
1192 897840 : if (lldo_trapped_waves) then
1193 :
1194 : ! Identify top edge of layer in which Scorer param drops below 0
1195 : ! - approximately the "turning level"
1196 : !----------------------------------------------------------
1197 0 : trn(:)=1.e8_r8
1198 0 : trn_level(:) = 0 ! pver+1
1199 0 : where( m2(:,pver)<= 0._r8 )
1200 0 : trn(:) = zi(:,pver)
1201 : trn_level(:) = pver
1202 : endwhere
1203 0 : do k=pver-1,1,-1
1204 0 : where( (m2(:,k+1)> 0._r8).and.(m2(:,k)<= 0._r8) )
1205 0 : trn(:) = zi(:,k)
1206 : trn_level(:) = k
1207 : endwhere
1208 : end do
1209 :
1210 0 : do i = 1,ncol
1211 : ! Case: Turning below mountain top
1212 0 : if ( (trn(i) < mxdis(i)).and.(trn_level(i)>=1) ) then
1213 0 : tau(i,0,:) = tau(i,0,:) - max( tauoro(i),taudsw(i) )
1214 0 : tau(i,0,:) = max( tau(i,0,:) , 0._r8 )
1215 0 : tau(i,0,1:tlb_level(i))=0._r8
1216 0 : src_level(i) = 1 ! disable any more tau calculation
1217 : end if
1218 : ! Case: Turning but no breaking
1219 0 : if ( (wbr(i) == 0._r8 ).and.(trn(i)>mxdis(i)).and.(trn_level(i)>=1) ) then
1220 0 : where ( ( zi(i,:) <= trn(i) ) .and. ( zi(i,:) >= bwv(i) ) )
1221 0 : tau(i,0,:) = tauoro(i) * &
1222 : ( trn(i) - zi(i,:) ) / &
1223 : ( trn(i) - bwv(i) )
1224 : end where
1225 0 : src_level(i) = 1 ! disable any more tau calculation
1226 : end if
1227 : ! Case: Turning AND breaking. Turning ABOVE breaking
1228 0 : if ( (wbr(i) > 0._r8 ).and.(trn(i) >= wbr(i)).and.(trn_level(i)>=1) ) then
1229 0 : where ( ( zi(i,:) <= trn(i) ) .and. ( zi(i,:) >= wbr(i) ) )
1230 0 : tau(i,0,:) = tauoro(i) * &
1231 : ( trn(i) - zi(i,:) ) / &
1232 : ( trn(i) - wbr(i) )
1233 : endwhere
1234 0 : src_level(i) = 1 ! disable any more tau calculation
1235 : end if
1236 : ! Case: Turning AND breaking. Turning BELOW breaking
1237 0 : if ( (wbr(i) > 0._r8 ).and.(trn(i) < wbr(i)).and.(trn_level(i)>=1) ) then
1238 0 : tauoro(i) = 0._r8
1239 0 : where ( ( zi(i,:) < wbr(i) ) .and. ( zi(i,:) >= tlb(i) ) )
1240 0 : tau(i,0,:) = tauoro(i) + (taudsw(i)-tauoro(i)) * &
1241 : ( wbr(i) - zi(i,:) ) / &
1242 : ( wbr(i) - tlb(i) )
1243 : endwhere
1244 0 : src_level(i) = 1 ! disable any more tau calculation
1245 : end if
1246 : end do
1247 : end if
1248 :
1249 897840 : end subroutine gw_rdg_break_trap
1250 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1251 :
1252 :
1253 :
1254 : end module gw_rdg
|
