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