Line data Source code
1 : module gw_convect
2 :
3 : !
4 : ! This module handles gravity waves from convection, and was extracted from
5 : ! gw_drag in May 2013.
6 : !
7 :
8 : use gw_utils, only: r8
9 :
10 : implicit none
11 : private
12 : save
13 :
14 : public :: BeresSourceDesc
15 : public :: gw_beres_src
16 :
17 : type :: BeresSourceDesc
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
26 : integer :: maxuh
27 : ! Heating depths [m].
28 : real(r8), allocatable :: hd(:)
29 : ! Table of source spectra.
30 : real(r8), allocatable :: mfcc(:,:,:)
31 : end type BeresSourceDesc
32 :
33 : contains
34 :
35 : !==========================================================================
36 :
37 1489176 : subroutine gw_beres_src(ncol, band, desc, u, v, &
38 2978352 : netdt, zm, src_level, tend_level, tau, ubm, ubi, xv, yv, &
39 1489176 : c, hdepth, maxq0)
40 : !-----------------------------------------------------------------------
41 : ! Driver for multiple gravity wave drag parameterization.
42 : !
43 : ! The parameterization is assumed to operate only where water vapor
44 : ! concentrations are negligible in determining the density.
45 : !
46 : ! Beres, J.H., M.J. Alexander, and J.R. Holton, 2004: "A method of
47 : ! specifying the gravity wave spectrum above convection based on latent
48 : ! heating properties and background wind". J. Atmos. Sci., Vol 61, No. 3,
49 : ! pp. 324-337.
50 : !
51 : !-----------------------------------------------------------------------
52 : use gw_utils, only: get_unit_vector, dot_2d, midpoint_interp
53 : use gw_common, only: GWBand, pver, qbo_hdepth_scaling
54 :
55 : !------------------------------Arguments--------------------------------
56 : ! Column dimension.
57 : integer, intent(in) :: ncol
58 :
59 : ! Wavelengths triggered by convection.
60 : type(GWBand), intent(in) :: band
61 :
62 : ! Settings for convection type (e.g. deep vs shallow).
63 : type(BeresSourceDesc), intent(in) :: desc
64 :
65 : ! Midpoint zonal/meridional winds.
66 : real(r8), intent(in) :: u(ncol,pver), v(ncol,pver)
67 : ! Heating rate due to convection.
68 : real(r8), intent(in) :: netdt(:,:)
69 : ! Midpoint altitudes.
70 : real(r8), intent(in) :: zm(ncol,pver)
71 :
72 : ! Indices of top gravity wave source level and lowest level where wind
73 : ! tendencies are allowed.
74 : integer, intent(out) :: src_level(ncol)
75 : integer, intent(out) :: tend_level(ncol)
76 :
77 : ! Wave Reynolds stress.
78 : real(r8), intent(out) :: tau(ncol,-band%ngwv:band%ngwv,pver+1)
79 : ! Projection of wind at midpoints and interfaces.
80 : real(r8), intent(out) :: ubm(ncol,pver), ubi(ncol,pver+1)
81 : ! Unit vectors of source wind (zonal and meridional components).
82 : real(r8), intent(out) :: xv(ncol), yv(ncol)
83 : ! Phase speeds.
84 : real(r8), intent(out) :: c(ncol,-band%ngwv:band%ngwv)
85 :
86 : ! Heating depth [m] and maximum heating in each column.
87 : real(r8), intent(out) :: hdepth(ncol), maxq0(ncol)
88 :
89 : !---------------------------Local Storage-------------------------------
90 : ! Column and level indices.
91 : integer :: i, k
92 :
93 : ! Zonal/meridional wind at roughly the level where the convection occurs.
94 2978352 : real(r8) :: uconv(ncol), vconv(ncol)
95 :
96 : ! Maximum heating rate.
97 2978352 : real(r8) :: q0(ncol)
98 :
99 : ! Bottom/top heating range index.
100 2978352 : integer :: boti(ncol), topi(ncol)
101 : ! Index for looking up heating depth dimension in the table.
102 2978352 : integer :: hd_idx(ncol)
103 : ! Mean wind in heating region.
104 2978352 : real(r8) :: uh(ncol)
105 : ! Min/max wavenumber for critical level filtering.
106 2978352 : integer :: Umini(ncol), Umaxi(ncol)
107 : ! Source level tau for a column.
108 2978352 : real(r8) :: tau0(-band%ngwv:band%ngwv)
109 : ! Speed of convective cells relative to storm.
110 2978352 : real(r8) :: CS(ncol)
111 : ! Index to shift spectra relative to ground.
112 : integer :: shift
113 :
114 : ! Heating rate conversion factor.
115 : real(r8), parameter :: CF = 20._r8
116 : ! Averaging length.
117 : real(r8), parameter :: AL = 1.0e5_r8
118 :
119 : !----------------------------------------------------------------------
120 : ! Initialize tau array
121 : !----------------------------------------------------------------------
122 :
123 >15207*10^7 : tau = 0.0_r8
124 24865776 : hdepth = 0.0_r8
125 24865776 : q0 = 0.0_r8
126 98285616 : tau0 = 0.0_r8
127 :
128 : !------------------------------------------------------------------------
129 : ! Determine wind and unit vectors approximately at the source level, then
130 : ! project winds.
131 : !------------------------------------------------------------------------
132 :
133 : ! Source wind speed and direction.
134 24865776 : uconv = u(:,desc%k)
135 24865776 : vconv = v(:,desc%k)
136 :
137 : ! Get the unit vector components and magnitude at the source level.
138 1489176 : call get_unit_vector(uconv, vconv, xv, yv, ubi(:,desc%k+1))
139 :
140 : ! Project the local wind at midpoints onto the source wind.
141 139982544 : do k = 1, pver
142 139982544 : ubm(:,k) = dot_2d(u(:,k), v(:,k), xv, yv)
143 : end do
144 :
145 : ! Compute the interface wind projection by averaging the midpoint winds.
146 : ! Use the top level wind at the top interface.
147 24865776 : ubi(:,1) = ubm(:,1)
148 :
149 1489176 : ubi(:,2:pver) = midpoint_interp(ubm)
150 :
151 : !-----------------------------------------------------------------------
152 : ! Calculate heating depth.
153 : !
154 : ! Heating depth is defined as the first height range from the bottom in
155 : ! which heating rate is continuously positive.
156 : !-----------------------------------------------------------------------
157 :
158 : ! First find the indices for the top and bottom of the heating range.
159 24865776 : boti = 0
160 24865776 : topi = 0
161 84061881 : do k = pver, 1, -1
162 1403715195 : do i = 1, ncol
163 1403715195 : if (boti(i) == 0) then
164 : ! Detect if we are outside the maximum range (where z = 20 km).
165 1115303792 : if (zm(i,k) >= 20000._r8) then
166 18921569 : boti(i) = k
167 18921569 : topi(i) = k
168 : else
169 : ! First spot where heating rate is positive.
170 1096382223 : if (netdt(i,k) > 0.0_r8) boti(i) = k
171 : end if
172 204349522 : else if (topi(i) == 0) then
173 : ! Detect if we are outside the maximum range (z = 20 km).
174 55645860 : if (zm(i,k) >= 20000._r8) then
175 0 : topi(i) = k
176 : else
177 : ! First spot where heating rate is no longer positive.
178 55645860 : if (.not. (netdt(i,k) > 0.0_r8)) topi(i) = k
179 : end if
180 : end if
181 : end do
182 : ! When all done, exit.
183 146927255 : if (all(topi /= 0)) exit
184 : end do
185 :
186 : ! Heating depth in m.
187 48242376 : hdepth = [ ( (zm(i,topi(i))-zm(i,boti(i))), i = 1, ncol ) ]
188 :
189 : ! J. Richter: this is an effective reduction of the GW phase speeds (needed to drive the QBO)
190 24865776 : hdepth = hdepth*qbo_hdepth_scaling
191 :
192 1489176 : hd_idx = index_of_nearest(hdepth, desc%hd)
193 :
194 : ! hd_idx=0 signals that a heating depth is too shallow, i.e. that it is
195 : ! either not big enough for the lowest table entry, or it is below the
196 : ! minimum allowed for this convection type.
197 : ! Values above the max in the table still get the highest value, though.
198 24865776 : where (hdepth < max(desc%min_hdepth, desc%hd(1))) hd_idx = 0
199 :
200 : ! Maximum heating rate.
201 79895164 : do k = minval(topi), maxval(boti)
202 531288868 : where (k >= topi .and. k <= boti)
203 31652788 : q0 = max(q0, netdt(:,k))
204 : end where
205 : end do
206 :
207 : !output max heating rate in K/day
208 24865776 : maxq0 = q0*24._r8*3600._r8
209 :
210 : ! Multipy by conversion factor
211 24865776 : q0 = q0 * CF
212 :
213 1489176 : if (desc%storm_shift) then
214 :
215 : ! Find the cell speed where the storm speed is > 10 m/s.
216 : ! Storm speed is taken to be the source wind speed.
217 24865776 : CS = sign(max(abs(ubm(:,desc%k))-10._r8, 0._r8), ubm(:,desc%k))
218 :
219 : ! Average wind in heating region, relative to storm cells.
220 24865776 : uh = 0._r8
221 79895164 : do k = minval(topi), maxval(boti)
222 531288868 : where (k >= topi .and. k <= boti)
223 31652788 : uh = uh + ubm(:,k)/(boti-topi+1)
224 : end where
225 : end do
226 :
227 24865776 : uh = uh - CS
228 :
229 : else
230 :
231 : ! For shallow convection, wind is relative to ground, and "heating
232 : ! region" wind is just the source level wind.
233 0 : uh = ubm(:,desc%k)
234 :
235 : end if
236 :
237 : ! Limit uh to table range.
238 24865776 : uh = min(uh, real(desc%maxuh, r8))
239 24865776 : uh = max(uh, -real(desc%maxuh, r8))
240 :
241 : ! Speeds for critical level filtering.
242 24865776 : Umini = band%ngwv
243 24865776 : Umaxi = -band%ngwv
244 79895164 : do k = minval(topi), maxval(boti)
245 1559235464 : where (k >= topi .and. k <= boti)
246 31652788 : Umini = min(Umini, nint(ubm(:,k)/band%dc))
247 : Umaxi = max(Umaxi, nint(ubm(:,k)/band%dc))
248 : end where
249 : end do
250 :
251 24865776 : Umini = max(Umini, -band%ngwv)
252 24865776 : Umaxi = min(Umaxi, band%ngwv)
253 :
254 : !-----------------------------------------------------------------------
255 : ! Gravity wave sources
256 : !-----------------------------------------------------------------------
257 : ! Start loop over all columns.
258 : !-----------------------------------------------------------------------
259 24865776 : do i=1,ncol
260 :
261 : !---------------------------------------------------------------------
262 : ! Look up spectrum only if the heating depth is large enough, else set
263 : ! tau0 = 0.
264 : !---------------------------------------------------------------------
265 :
266 24865776 : if (hd_idx(i) > 0) then
267 :
268 : !------------------------------------------------------------------
269 : ! Look up the spectrum using depth and uh.
270 : !------------------------------------------------------------------
271 :
272 67079100 : tau0 = desc%mfcc(hd_idx(i),nint(uh(i)),:)
273 :
274 1016350 : if (desc%storm_shift) then
275 : ! For deep convection, the wind was relative to storm cells, so
276 : ! shift the spectrum so that it is now relative to the ground.
277 1016350 : shift = -nint(CS(i)/band%dc)
278 67079100 : tau0 = eoshift(tau0, shift)
279 : end if
280 :
281 : ! Adjust magnitude.
282 67079100 : tau0 = tau0*q0(i)*q0(i)/AL
283 :
284 : ! Adjust for critical level filtering.
285 6346713 : tau0(Umini(i):Umaxi(i)) = 0.0_r8
286 :
287 67079100 : tau(i,:,topi(i)+1) = tau0
288 :
289 : end if ! heating depth above min and not at the pole
290 :
291 : enddo
292 :
293 : !-----------------------------------------------------------------------
294 : ! End loop over all columns.
295 : !-----------------------------------------------------------------------
296 :
297 : ! Output the source level.
298 24865776 : src_level = topi
299 24865776 : tend_level = topi
300 :
301 : ! Set phase speeds; just use reference speeds.
302 1489176 : c = spread(band%cref, 1, ncol)
303 :
304 1489176 : end subroutine gw_beres_src
305 :
306 : ! Short routine to get the indices of a set of values rounded to their
307 : ! nearest points on a grid.
308 1489176 : function index_of_nearest(x, grid) result(idx)
309 : real(r8), intent(in) :: x(:)
310 : real(r8), intent(in) :: grid(:)
311 :
312 : integer :: idx(size(x))
313 :
314 2978352 : real(r8) :: interfaces(size(grid)-1)
315 : integer :: i, n
316 :
317 1489176 : n = size(grid)
318 31272696 : interfaces = (grid(:n-1) + grid(2:))/2._r8
319 :
320 24865776 : idx = 1
321 29783520 : do i = 1, n-1
322 473938920 : where (x > interfaces(i)) idx = i + 1
323 : end do
324 :
325 1489176 : end function index_of_nearest
326 :
327 0 : end module gw_convect
|
