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 58824 : subroutine gw_beres_src(ncol, band, desc, u, v, &
38 117648 : netdt, zm, src_level, tend_level, tau, ubm, ubi, xv, yv, &
39 58824 : 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 117648 : real(r8) :: uconv(ncol), vconv(ncol)
95 :
96 : ! Maximum heating rate.
97 117648 : real(r8) :: q0(ncol)
98 :
99 : ! Bottom/top heating range index.
100 117648 : integer :: boti(ncol), topi(ncol)
101 : ! Index for looking up heating depth dimension in the table.
102 117648 : integer :: hd_idx(ncol)
103 : ! Mean wind in heating region.
104 117648 : real(r8) :: uh(ncol)
105 : ! Min/max wavenumber for critical level filtering.
106 117648 : integer :: Umini(ncol), Umaxi(ncol)
107 : ! Source level tau for a column.
108 117648 : real(r8) :: tau0(-band%ngwv:band%ngwv)
109 : ! Speed of convective cells relative to storm.
110 117648 : 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 6006976920 : tau = 0.0_r8
124 982224 : hdepth = 0.0_r8
125 982224 : q0 = 0.0_r8
126 3882384 : 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 982224 : uconv = u(:,desc%k)
135 982224 : vconv = v(:,desc%k)
136 :
137 : ! Get the unit vector components and magnitude at the source level.
138 58824 : 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 5529456 : do k = 1, pver
142 5529456 : 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 982224 : ubi(:,1) = ubm(:,1)
148 :
149 58824 : 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 982224 : boti = 0
160 982224 : topi = 0
161 3302182 : do k = pver, 1, -1
162 55146871 : do i = 1, ncol
163 55146871 : if (boti(i) == 0) then
164 : ! Detect if we are outside the top of range (where z = 20 km).
165 44832201 : if (zm(i,k) >= 20000._r8) then
166 763752 : boti(i) = k
167 763752 : topi(i) = k
168 : else
169 : ! First spot where heating rate is positive.
170 44068449 : if (netdt(i,k) > 0.0_r8) boti(i) = k
171 : end if
172 : end if
173 : end do
174 : ! When all done, exit
175 6181257 : if (all(boti /= 0)) exit
176 : end do
177 :
178 1727181 : do k = 1, pver
179 28899489 : do i = 1, ncol
180 28899489 : if (topi(i) == 0) then
181 : ! First spot where heating rate is positive.
182 11731624 : if ((netdt(i,k) > 0.0_r8) .AND. (zm(i,k) <= 20000._r8)) topi(i) = k-1
183 : end if
184 : end do
185 : ! When all done, exit
186 7800910 : if (all(topi /= 0)) exit
187 : end do
188 :
189 : ! Heating depth in m.
190 1905624 : hdepth = [ ( (zm(i,topi(i))-zm(i,boti(i))), i = 1, ncol ) ]
191 :
192 : ! J. Richter: this is an effective reduction of the GW phase speeds (needed to drive the QBO)
193 982224 : hdepth = hdepth*qbo_hdepth_scaling
194 :
195 58824 : hd_idx = index_of_nearest(hdepth, desc%hd)
196 :
197 : ! hd_idx=0 signals that a heating depth is too shallow, i.e. that it is
198 : ! either not big enough for the lowest table entry, or it is below the
199 : ! minimum allowed for this convection type.
200 : ! Values above the max in the table still get the highest value, though.
201 982224 : where (hdepth < max(desc%min_hdepth, desc%hd(1))) hd_idx = 0
202 :
203 : ! Maximum heating rate.
204 3060879 : do k = minval(topi), maxval(boti)
205 19394952 : where (k >= topi .and. k <= boti)
206 1155255 : q0 = max(q0, netdt(:,k))
207 : end where
208 : end do
209 :
210 : !output max heating rate in K/day
211 982224 : maxq0 = q0*24._r8*3600._r8
212 :
213 : ! Multipy by conversion factor
214 982224 : q0 = q0 * CF
215 :
216 58824 : if (desc%storm_shift) then
217 :
218 : ! Find the cell speed where the storm speed is > 10 m/s.
219 : ! Storm speed is taken to be the source wind speed.
220 982224 : CS = sign(max(abs(ubm(:,desc%k))-10._r8, 0._r8), ubm(:,desc%k))
221 :
222 : ! Average wind in heating region, relative to storm cells.
223 982224 : uh = 0._r8
224 3060879 : do k = minval(topi), maxval(boti)
225 19394952 : where (k >= topi .and. k <= boti)
226 1155255 : uh = uh + ubm(:,k)/(boti-topi+1)
227 : end where
228 : end do
229 :
230 982224 : uh = uh - CS
231 :
232 : else
233 :
234 : ! For shallow convection, wind is relative to ground, and "heating
235 : ! region" wind is just the source level wind.
236 0 : uh = ubm(:,desc%k)
237 :
238 : end if
239 :
240 : ! Limit uh to table range.
241 982224 : uh = min(uh, real(desc%maxuh, r8))
242 982224 : uh = max(uh, -real(desc%maxuh, r8))
243 :
244 : ! Speeds for critical level filtering.
245 982224 : Umini = band%ngwv
246 982224 : Umaxi = -band%ngwv
247 3060879 : do k = minval(topi), maxval(boti)
248 56911953 : where (k >= topi .and. k <= boti)
249 1155255 : Umini = min(Umini, nint(ubm(:,k)/band%dc))
250 : Umaxi = max(Umaxi, nint(ubm(:,k)/band%dc))
251 : end where
252 : end do
253 :
254 982224 : Umini = max(Umini, -band%ngwv)
255 982224 : Umaxi = min(Umaxi, band%ngwv)
256 :
257 : !-----------------------------------------------------------------------
258 : ! Gravity wave sources
259 : !-----------------------------------------------------------------------
260 : ! Start loop over all columns.
261 : !-----------------------------------------------------------------------
262 982224 : do i=1,ncol
263 :
264 : !---------------------------------------------------------------------
265 : ! Look up spectrum only if the heating depth is large enough, else set
266 : ! tau0 = 0.
267 : !---------------------------------------------------------------------
268 :
269 982224 : if (hd_idx(i) > 0) then
270 :
271 : !------------------------------------------------------------------
272 : ! Look up the spectrum using depth and uh.
273 : !------------------------------------------------------------------
274 :
275 3423024 : tau0 = desc%mfcc(hd_idx(i),nint(uh(i)),:)
276 :
277 51864 : if (desc%storm_shift) then
278 : ! For deep convection, the wind was relative to storm cells, so
279 : ! shift the spectrum so that it is now relative to the ground.
280 51864 : shift = -nint(CS(i)/band%dc)
281 3423024 : tau0 = eoshift(tau0, shift)
282 : end if
283 :
284 : ! Adjust magnitude.
285 3423024 : tau0 = tau0*q0(i)*q0(i)/AL
286 :
287 : ! Adjust for critical level filtering.
288 330826 : tau0(Umini(i):Umaxi(i)) = 0.0_r8
289 :
290 3423024 : tau(i,:,topi(i)+1) = tau0
291 :
292 : end if ! heating depth above min and not at the pole
293 :
294 : enddo
295 :
296 : !-----------------------------------------------------------------------
297 : ! End loop over all columns.
298 : !-----------------------------------------------------------------------
299 :
300 : ! Output the source level.
301 982224 : src_level = topi
302 982224 : tend_level = topi
303 :
304 : ! Set phase speeds; just use reference speeds.
305 58824 : c = spread(band%cref, 1, ncol)
306 :
307 58824 : end subroutine gw_beres_src
308 :
309 : ! Short routine to get the indices of a set of values rounded to their
310 : ! nearest points on a grid.
311 58824 : function index_of_nearest(x, grid) result(idx)
312 : real(r8), intent(in) :: x(:)
313 : real(r8), intent(in) :: grid(:)
314 :
315 : integer :: idx(size(x))
316 :
317 117648 : real(r8) :: interfaces(size(grid)-1)
318 : integer :: i, n
319 :
320 58824 : n = size(grid)
321 1235304 : interfaces = (grid(:n-1) + grid(2:))/2._r8
322 :
323 982224 : idx = 1
324 1176480 : do i = 1, n-1
325 18721080 : where (x > interfaces(i)) idx = i + 1
326 : end do
327 :
328 58824 : end function index_of_nearest
329 :
330 0 : end module gw_convect
|
