Line data Source code
1 : module gw_movmtn
2 :
3 : !
4 : ! This module parameterizes gravity waves generated by the obstacle effect produced by
5 : ! boundary layer turbulence for convection.
6 : !
7 :
8 : use gw_utils, only: r8
9 :
10 : implicit none
11 : private
12 : save
13 :
14 : public :: MovMtnSourceDesc
15 : public :: gw_movmtn_src
16 :
17 : type :: MovMtnSourceDesc
18 : ! Whether wind speeds are shifted to be relative to storm cells.
19 : logical :: storm_shift
20 : ! Index for level where wind speed is used as the source speed.
21 : integer :: k
22 : ! Heating depths below this value [m] will be ignored.
23 : real(r8) :: min_hdepth
24 : ! Table bounds, for convenience. (Could be inferred from shape(mfcc).)
25 : integer :: maxh !-bounds of the lookup table heating depths
26 : integer :: maxuh ! bounds of the lookup table wind
27 : ! Heating depths [m].
28 : real(r8), allocatable :: hd(:), uh(:)
29 : ! Table of source spectra.
30 : real(r8), allocatable :: mfcc(:,:,:) !is the lookup table f(depth, wind, phase speed)
31 : end type MovMtnSourceDesc
32 :
33 : contains
34 :
35 : !==========================================================================
36 :
37 2978352 : subroutine gw_movmtn_src(ncol,lchnk, band, desc, u, v, &
38 2978352 : netdt, netdt_shcu, xpwp_shcu, &
39 2978352 : zm, alpha_gw_movmtn, src_level, tend_level, tau, ubm, ubi, xv, yv, &
40 1489176 : c, hdepth)
41 : !-----------------------------------------------------------------------
42 : ! Flexible driver for gravity wave source from obstacle effects produced
43 : ! by boundary layer turbulence or deep convection
44 : !-----------------------------------------------------------------------
45 : use gw_utils, only: get_unit_vector, dot_2d, midpoint_interp
46 : use gw_common, only: GWBand, pver, qbo_hdepth_scaling
47 : use cam_history, only: outfld
48 : use phys_control, only: use_gw_movmtn_pbl
49 : use physconst, only: rair, gravit
50 : !------------------------------Arguments--------------------------------
51 : ! Column dimension.
52 : integer, intent(in) :: ncol , lchnk
53 :
54 : ! Wavelengths triggered by convection.
55 : type(GWBand), intent(in) :: band
56 :
57 : ! Settings for convection type (e.g. deep vs shallow).
58 : type(MovMtnSourceDesc), intent(in) :: desc
59 :
60 : ! Midpoint zonal/meridional winds.
61 : real(r8), intent(in) :: u(ncol,pver), v(ncol,pver)
62 : ! Heating rate due to convection.
63 : real(r8), intent(in) :: netdt(:,:) !from deep scheme
64 : ! Heating rate due to shallow convection and PBL turbulence.
65 : real(r8), intent(in) :: netdt_shcu(:,:)
66 : ! Higher order flux from ShCu/PBL.
67 : real(r8), intent(in) :: xpwp_shcu(ncol,pver+1)
68 : ! Midpoint altitudes.
69 : real(r8), intent(in) :: zm(ncol,pver)
70 : ! tunable parameter controlling proportion of PBL momentum flux emitted as GW
71 : real(r8), intent(in) :: alpha_gw_movmtn
72 :
73 : ! Indices of top gravity wave source level and lowest level where wind
74 : ! tendencies are allowed.
75 : integer, intent(out) :: src_level(ncol)
76 : integer, intent(out) :: tend_level(ncol)
77 :
78 : ! Wave Reynolds stress.
79 : real(r8), intent(out) :: tau(ncol,-band%ngwv:band%ngwv,pver+1) !tau = momentum flux (m2/s2) at interface level ngwv = band of phase speeds
80 : ! Projection of wind at midpoints and interfaces.
81 : real(r8), intent(out) :: ubm(ncol,pver), ubi(ncol,pver+1)
82 : ! Unit vectors of source wind (zonal and meridional components).
83 : real(r8), intent(out) :: xv(ncol), yv(ncol) !determined by vector direction of wind at source
84 : ! Phase speeds.
85 : real(r8), intent(out) :: c(ncol,-band%ngwv:band%ngwv)
86 :
87 : ! Heating depth [m] and maximum heating in each column.
88 : real(r8), intent(out) :: hdepth(ncol) !calculated here in this code
89 :
90 : !---------------------------Local Storage-------------------------------
91 : ! Column and (vertical) level indices.
92 : integer :: i, k
93 :
94 : ! Zonal/meridional wind at steering level, i.e., 'cell speed'.
95 : ! May be later modified by retrograde motion ....
96 2978352 : real(r8) :: usteer(ncol), vsteer(ncol)
97 2978352 : real(r8) :: uwavef(ncol,pver),vwavef(ncol,pver)
98 : ! Steering level (integer converted to real*8)
99 2978352 : real(r8) :: steer_level(ncol)
100 : ! Retrograde motion of Cell
101 2978352 : real(r8) :: Cell_Retro_Speed(ncol)
102 :
103 : ! Maximum heating rate.
104 2978352 : real(r8) :: q0(ncol), qj(ncol)
105 : ! unit vector components at steering level and mag
106 2978352 : real(r8) :: xv_steer(ncol), yv_steer(ncol), umag_steer(ncol)
107 : ! Bottom/top heating range index.
108 2978352 : integer :: boti(ncol), topi(ncol)
109 : ! Index for looking up heating depth dimension in the table.
110 2978352 : integer :: hd_idx(ncol)
111 : ! Mean wind in heating region.
112 2978352 : real(r8) :: uh(ncol)
113 : ! Min/max wavenumber for critical level filtering.
114 : integer :: Umini(ncol), Umaxi(ncol)
115 : ! Source level tau for a column.
116 2978352 : real(r8) :: tau0(-band%ngwv:band%ngwv)
117 : ! Speed of convective cells relative to storm.
118 2978352 : real(r8) :: CS(ncol),CS1(ncol)
119 : ! Wind speeds in wave direction
120 2978352 : real(r8) :: udiff(ncol),vdiff(ncol)
121 : ! "on-crest" source level wind
122 2978352 : real(r8) :: ubmsrc(ncol),ubisrc(ncol)
123 :
124 : ! Index to shift spectra relative to ground.
125 : integer :: shift
126 : ! Other wind quantities
127 2978352 : real(r8) :: ut(ncol),uc(ncol),umm(ncol)
128 : ! Tau from moving mountain lookup table
129 2978352 : real(r8) :: taumm(ncol)
130 : ! Heating rate conversion factor. -> tuning factors
131 : real(r8), parameter :: CF = 20._r8 !(1/ (5%)) -> 5% of grid cell is covered with convection
132 : ! Averaging length.
133 : real(r8), parameter :: AL = 1.0e5_r8
134 : ! Index for moving mountain lookuptable
135 2978352 : integer :: hdmm_idx(ncol), uhmm_idx(ncol)
136 : ! Index for ground based phase speed bin
137 2978352 : real(r8) :: c0(ncol,-band%ngwv:band%ngwv)
138 2978352 : integer :: c_idx(ncol,-band%ngwv:band%ngwv)
139 : ! Flux source from ShCu/PBL
140 1489176 : real(r8) :: xpwp_src(ncol)
141 : ! Manual steering level set
142 : integer :: Steer_k
143 :
144 : !----------------------------------------------------------------------
145 : ! Initialize tau array
146 : !----------------------------------------------------------------------
147 2478854664 : tau = 0.0_r8
148 24865776 : hdepth = 0.0_r8
149 24865776 : q0 = 0.0_r8
150 2978352 : tau0 = 0.0_r8
151 :
152 : !----------------------------------------------------------------------
153 : ! Calculate flux source from ShCu/PBL
154 : !----------------------------------------------------------------------
155 1489176 : xpwp_src = shcu_flux_src( xpwp_shcu, ncol, pver+1, alpha_gw_movmtn )
156 :
157 : !------------------------------------------------------------------------
158 : ! Determine wind and unit vectors approximately at the source (steering level), then
159 : ! project winds.
160 : !------------------------------------------------------------------------
161 :
162 : ! Winds at 'steering level'
163 1489176 : Steer_k = pver-1
164 24865776 : usteer = u(:,Steer_k) !k defined in line21 (at specified altitude)
165 24865776 : vsteer = v(:,Steer_k)
166 24865776 : steer_level = real(Steer_k,r8)
167 :
168 : ! all GW calculations on a plane, which in our case is the wind at source level -> ubi is wind in this plane
169 : ! Get the unit vector components and magnitude at the source level.
170 1489176 : call get_unit_vector(usteer, vsteer, xv_steer, yv_steer, umag_steer)
171 :
172 : !-------------------------------------------------------------------------
173 : ! If we want to account for some retorgrade cell motion,
174 : ! it should be done by vector subtraction from (usteer,vsteer).
175 : ! We assume the retrograde motion is in the same direction as
176 : ! (usteer,vsteer) or the unit vector (xv_steer,yv_steer). Then, the
177 : ! vector retrograde motion is just:
178 : ! = -Cell_Retrograde_Speed * (xv_steer,yv_steer)
179 : ! and we would modify usteer and vsteer
180 : ! usteer = usteer - Cell_Retrograde_Speed * xv_steer
181 : ! vsteer = vsteer - Cell_Retrograde_Speed * yv_steer
182 : !-----------------------------------------------------------------------
183 : ! Cell_Retro_Speed is always =0 for now
184 : !-----------------------------------------------------------------------
185 24865776 : do i=1,ncol
186 24865776 : Cell_Retro_Speed(i) = min( sqrt(usteer(i)**2 + vsteer(i)**2), 0._r8)
187 : end do
188 24865776 : do i=1,ncol
189 23376600 : usteer(i) = usteer(i) - xv_steer(i)*Cell_Retro_Speed(i)
190 24865776 : vsteer(i) = vsteer(i) - yv_steer(i)*Cell_Retro_Speed(i)
191 : end do
192 : !-------------------------------------------------------------------------
193 : ! At this point (usteer,vsteer) is the cell-speed, or equivalently, the 2D
194 : ! ground based wave phase speed for moving mountain GW
195 : !-------------------------------------------------------------------------
196 :
197 :
198 : ! Calculate heating depth.
199 : !
200 : ! Heating depth is defined as the first height range from the bottom in
201 : ! which heating rate is continuously positive.
202 : !-----------------------------------------------------------------------
203 :
204 : ! First find the indices for the top and bottom of the heating range.
205 : !nedt is heating profile from Zhang McFarlane (it's pressure coordinates, therefore k=0 is the top)
206 :
207 24865776 : boti = 0 !bottom
208 24865776 : topi = 0 !top
209 :
210 1489176 : if (use_gw_movmtn_pbl) then
211 24865776 : boti=pver
212 24865776 : topi=Steer_k-10 ! desc%k-5
213 : else
214 0 : do k = pver, 1, -1 !start at surface
215 0 : do i = 1, ncol
216 0 : if (boti(i) == 0) then
217 : ! Detect if we are outside the maximum range (where z = 20 km).
218 0 : if (zm(i,k) >= 20000._r8) then
219 0 : boti(i) = k
220 0 : topi(i) = k
221 : else
222 : ! First spot where heating rate is positive.
223 0 : if (netdt(i,k) > 0.0_r8) boti(i) = k
224 : end if
225 0 : else if (topi(i) == 0) then
226 : ! Detect if we are outside the maximum range (z = 20 km).
227 0 : if (zm(i,k) >= 20000._r8) then
228 0 : topi(i) = k
229 : else
230 : ! First spot where heating rate is no longer positive.
231 0 : if (.not. (netdt(i,k) > 0.0_r8)) topi(i) = k
232 : end if
233 : end if
234 : end do
235 : ! When all done, exit.
236 0 : if (all(topi /= 0)) exit
237 : end do
238 : end if
239 : ! Heating depth in m. (top-bottom altitudes)
240 48242376 : hdepth = [ ( (zm(i,topi(i))-zm(i,boti(i))), i = 1, ncol ) ]
241 1489176 : hd_idx = index_of_nearest(hdepth, desc%hd)
242 :
243 : ! hd_idx=0 signals that a heating depth is too shallow, i.e. that it is
244 : ! either not big enough for the lowest table entry, or it is below the
245 : ! minimum allowed for this convection type.
246 : ! Values above the max in the table still get the highest value, though.
247 :
248 24865776 : where (hdepth < max(desc%min_hdepth, desc%hd(1))) hd_idx = 0
249 :
250 : ! Maximum heating rate.
251 66112488 : do k = minval(topi), maxval(boti)
252 299878488 : where (k >= topi .and. k <= boti)
253 17870112 : q0 = max(q0, netdt(:,k))
254 : end where
255 : end do
256 :
257 : ! Multiply by conversion factor
258 : ! (now 20* larger than what Zhang McFarlane said as they try to describe heating over 100km grid cell)
259 24865776 : q0 = q0 * CF
260 24865776 : qj = gravit/rair*q0 ! unit conversion to m/s3
261 :
262 : !-------------------------------------------------
263 : ! CS1 and CS should be equal in current implemen-
264 : ! tation.
265 : !-------------------------------------------------
266 24865776 : CS1 = sqrt( usteer**2._r8 + vsteer**2._r8 )
267 24865776 : CS = CS1*xv_steer + CS1*yv_steer
268 :
269 : ! -----------------------------------------------------------
270 : ! Calculate winds in reference frame of wave (uwavef,vwavef).
271 : ! This is like "(U-c)" in GW literature, where U and c are in
272 : ! ground-based speeds in a plane perpendicular to wave fronts.
273 : !------------------------------------------------------------
274 24865776 : do i=1,ncol
275 23376600 : udiff(i) = u(i,topi(i)) - usteer(i)
276 23376600 : vdiff(i) = v(i,topi(i)) - vsteer(i)
277 2198889576 : do k=1,pver
278 2174023800 : uwavef(i, k ) = u(i, k ) - usteer(i)
279 2197400400 : vwavef(i, k ) = v(i, k ) - vsteer(i)
280 : end do
281 : end do
282 : !----------------------------------------------------------
283 : ! Wave relative wind at source level. This determines
284 : ! orientation of wave in the XY plane, and therefore the
285 : ! direction in which force from dissipating GW will be
286 : ! applied.
287 : !----------------------------------------------------------
288 24865776 : do i=1,ncol
289 23376600 : udiff(i) = uwavef( i, topi(i) )
290 24865776 : vdiff(i) = vwavef( i, topi(i) )
291 : end do
292 : !-----------------------------------------------------------
293 : ! Unit vector components (xv,yv) in direction of wavevector
294 : ! i.e., in which force will be applied
295 : !-----------------------------------------------------------
296 1489176 : call get_unit_vector(udiff , vdiff , xv, yv, ubisrc )
297 :
298 1489176 : call outfld('UCELL_MOVMTN', usteer, ncol, lchnk)
299 1489176 : call outfld('VCELL_MOVMTN', vsteer, ncol, lchnk)
300 1489176 : call outfld('CS_MOVMTN', CS, ncol, lchnk)
301 1489176 : call outfld('STEER_LEVEL_MOVMTN',steer_level, ncol, lchnk )
302 1489176 : call outfld('XPWP_SRC_MOVMTN', xpwp_src , ncol, lchnk )
303 :
304 : !----------------------------------------------------------
305 : ! Project the local wave relative wind at midpoints onto the
306 : ! direction of the wavevector.
307 : !----------------------------------------------------------
308 139982544 : do k = 1, pver
309 139982544 : ubm(:,k) = dot_2d(uwavef(:,k), vwavef(:,k), xv, yv)
310 : end do
311 : ! Source level on-crest wind
312 24865776 : do i=1,ncol
313 24865776 : ubmsrc(i) = ubm(i,topi(i))
314 : end do
315 :
316 : !---------------------------------------------------------------
317 : ! adjust everything so that source level wave relative on-crest
318 : ! wind is always positive. Also adjust unit vector comps xv,yv
319 : !--------------------------------------------------------------
320 139982544 : do k=1,pver
321 2314006344 : do i=1,ncol
322 2312517168 : ubm(i,k) = sign( 1._r8 , ubmsrc(i) )* ubm(i,k)
323 : end do
324 : end do
325 : !
326 24865776 : do i=1,ncol
327 23376600 : xv(i) = sign( 1._r8 , ubmsrc(i) ) * xv(i)
328 24865776 : yv(i) = sign( 1._r8 , ubmsrc(i) ) * yv(i)
329 : end do
330 :
331 :
332 :
333 : ! Compute the interface wind projection by averaging the midpoint winds. (both same wind profile,
334 : ! just at different points of the grid)
335 :
336 : ! Use the top level wind at the top interface.
337 24865776 : ubi(:,1) = ubm(:,1)
338 :
339 1489176 : ubi(:,2:pver) = midpoint_interp(ubm)
340 :
341 : !-----------------------------------------------------------------------
342 : ! determine wind for lookup table
343 : ! need wind speed at the top of the convecitve cell and at the steering level
344 24865776 : uh = 0._r8
345 24865776 : do i=1,ncol
346 23376600 : ut(i) = ubm(i,topi(i))
347 24865776 : uh(i) = ut(i) - CS(i) ! wind at top in the frame moving with the cell
348 : end do
349 :
350 : ! Set phase speeds; just use reference speeds.
351 24865776 : c(:,0) = 0._r8
352 :
353 : !-----------------------------------------------------------------------
354 : ! Gravity wave sources
355 : !-----------------------------------------------------------------------
356 : ! Start loop over all columns.
357 : !-----------------------------------------------------------------------
358 24865776 : do i=1,ncol
359 :
360 : !---------------------------------------------------------------------
361 : ! Look up spectrum only if the heating depth is large enough, else leave
362 : ! tau = 0.
363 : !---------------------------------------------------------------------
364 24865776 : if (.not. use_gw_movmtn_pbl) then
365 0 : if (hd_idx(i) > 0) then
366 : !------------------------------------------------------------------
367 : ! Look up the spectrum using depth and uh.
368 : !------------------------------------------------------------------
369 : !hdmm_idx = index_of_nearest(hdepth, desc%hd)
370 0 : uhmm_idx = index_of_nearest(uh, desc%uh)
371 0 : taumm(i) = abs(desc%mfcc(uhmm_idx(i),hd_idx(i),0))
372 0 : taumm(i) = taumm(i)*qj(i)*qj(i)/AL/1000._r8
373 : ! assign sign to MF based on the ground based phase speed, ground based phase speed = CS
374 0 : taumm(i) = -1._r8*sign(taumm(i),CS(i))
375 : !find the right phase speed bin
376 0 : c0(i,:) = CS(i)
377 0 : c_idx(i,:) = index_of_nearest(c0(i,:),c(i,:))
378 :
379 : !input tau to top +1 level, interface level just below top of heating, remember it's in pressure
380 : ! everything is upside down (source level of GWs, level where GWs are launched)
381 0 : tau(i,c_idx(i,:),topi(i):topi(i)+1) = taumm(i)
382 :
383 : end if ! heating depth above min and not at the pole
384 : else
385 327272400 : tau(i,0,topi(i):pver+1 ) = xpwp_src(i) ! 0.1_r8/10000._r8
386 : endif
387 :
388 : enddo
389 : !-----------------------------------------------------------------------
390 : ! End loop over all columns.
391 : !-----------------------------------------------------------------------
392 :
393 : ! Output the source level.
394 24865776 : src_level = topi
395 24865776 : tend_level = topi
396 :
397 :
398 1489176 : end subroutine gw_movmtn_src
399 :
400 : ! Short routine to get the indices of a set of values rounded to their
401 : ! nearest points on a grid.
402 1489176 : pure function index_of_nearest(x, grid) result(idx)
403 : real(r8), intent(in) :: x(:)
404 : real(r8), intent(in) :: grid(:)
405 :
406 : integer :: idx(size(x))
407 :
408 2978352 : real(r8) :: interfaces(size(grid)-1)
409 : integer :: i, n
410 :
411 1489176 : n = size(grid)
412 23826816 : interfaces = (grid(:n-1) + grid(2:))/2._r8
413 :
414 24865776 : idx = 1
415 22337640 : do i = 1, n-1
416 349610040 : where (x > interfaces(i)) idx = i + 1
417 : end do
418 :
419 2978352 : end function index_of_nearest
420 :
421 : !!!!!!!!!!!!!!!!!!!!!!!!!!!
422 1489176 : pure function shcu_flux_src (xpwp_shcu , ncol, pverx, alpha_gw_movmtn ) result(xpwp_src)
423 : integer, intent(in) :: ncol,pverx
424 : real(r8), intent(in) :: xpwp_shcu (ncol,pverx)
425 : real(r8), intent(in) :: alpha_gw_movmtn
426 :
427 : real(r8) :: xpwp_src(ncol)
428 :
429 : integer :: k, nlayers
430 :
431 : !-----------------------------------
432 : ! Simple average over layers.
433 : ! Probably can do better
434 : !-----------------------------------
435 : nlayers=5
436 24865776 : xpwp_src(:) =0._r8
437 8935056 : do k = 0, nlayers-1
438 125818056 : xpwp_src(:) = xpwp_src(:) + xpwp_shcu(:,pverx-k)
439 : end do
440 24865776 : xpwp_src(:) = alpha_gw_movmtn * xpwp_src(:)/(1.0_r8*nlayers)
441 :
442 1489176 : end function shcu_flux_src
443 :
444 0 : end module gw_movmtn
|
