1 | |
---|
2 | ! $Id: hines_gwd.F90 2408 2015-12-14 10:43:09Z fcheruy $ |
---|
3 | |
---|
4 | SUBROUTINE hines_gwd(nlon, nlev, dtime, paphm1x, papm1x, rlat, tx, ux, vx, & |
---|
5 | zustrhi, zvstrhi, d_t_hin, d_u_hin, d_v_hin) |
---|
6 | |
---|
7 | ! ######################################################################## |
---|
8 | ! Parametrization of the momentum flux deposition due to a broad band |
---|
9 | ! spectrum of gravity waves, following Hines (1997a,b), as coded by |
---|
10 | ! McLANDRESS (1995). Modified by McFARLANE and MANZINI (1995-1997) |
---|
11 | ! MAECHAM model stand alone version |
---|
12 | ! ######################################################################## |
---|
13 | |
---|
14 | |
---|
15 | USE dimphy |
---|
16 | IMPLICIT NONE |
---|
17 | |
---|
18 | include "YOEGWD.h" |
---|
19 | include "YOMCST.h" |
---|
20 | |
---|
21 | INTEGER nazmth |
---|
22 | PARAMETER (nazmth=8) |
---|
23 | |
---|
24 | ! INPUT ARGUMENTS. |
---|
25 | ! ----- ---------- |
---|
26 | |
---|
27 | ! - 2D |
---|
28 | ! PAPHM1 : HALF LEVEL PRESSURE (T-DT) |
---|
29 | ! PAPM1 : FULL LEVEL PRESSURE (T-DT) |
---|
30 | ! PTM1 : TEMPERATURE (T-DT) |
---|
31 | ! PUM1 : ZONAL WIND (T-DT) |
---|
32 | ! PVM1 : MERIDIONAL WIND (T-DT) |
---|
33 | |
---|
34 | |
---|
35 | ! REFERENCE. |
---|
36 | ! ---------- |
---|
37 | ! SEE MODEL DOCUMENTATION |
---|
38 | |
---|
39 | ! AUTHOR. |
---|
40 | ! ------- |
---|
41 | |
---|
42 | ! N. MCFARLANE DKRZ-HAMBURG MAY 1995 |
---|
43 | ! STAND ALONE E. MANZINI MPI-HAMBURG FEBRUARY 1997 |
---|
44 | |
---|
45 | ! BASED ON A COMBINATION OF THE OROGRAPHIC SCHEME BY N.MCFARLANE 1987 |
---|
46 | ! AND THE HINES SCHEME AS CODED BY C. MCLANDRESS 1995. |
---|
47 | |
---|
48 | |
---|
49 | |
---|
50 | ! ym INTEGER KLEVM1 |
---|
51 | |
---|
52 | REAL paphm1(klon, klev+1), papm1(klon, klev) |
---|
53 | REAL ptm1(klon, klev), pum1(klon, klev), pvm1(klon, klev) |
---|
54 | REAL prflux(klon) |
---|
55 | ! 1 |
---|
56 | ! 1 |
---|
57 | ! 1 |
---|
58 | REAL rlat(klon), coslat(klon) |
---|
59 | |
---|
60 | REAL th(klon, klev), utendgw(klon, klev), vtendgw(klon, klev), & |
---|
61 | pressg(klon), uhs(klon, klev), vhs(klon, klev), zpr(klon) |
---|
62 | |
---|
63 | ! * VERTICAL POSITIONING ARRAYS. |
---|
64 | |
---|
65 | REAL sgj(klon, klev), shj(klon, klev), shxkj(klon, klev), dsgj(klon, klev) |
---|
66 | |
---|
67 | ! * LOGICAL SWITCHES TO CONTROL ROOF DRAG, ENVELOP GW DRAG AND |
---|
68 | ! * HINES' DOPPLER SPREADING EXTROWAVE GW DRAG. |
---|
69 | ! * LOZPR IS TRUE FOR ZPR ENHANCEMENT |
---|
70 | |
---|
71 | |
---|
72 | ! * WORK ARRAYS. |
---|
73 | |
---|
74 | REAL m_alpha(klon, klev, nazmth), v_alpha(klon, klev, nazmth), & |
---|
75 | sigma_alpha(klon, klev, nazmth), sigsqh_alpha(klon, klev, nazmth), & |
---|
76 | drag_u(klon, klev), drag_v(klon, klev), flux_u(klon, klev), & |
---|
77 | flux_v(klon, klev), heat(klon, klev), diffco(klon, klev), & |
---|
78 | bvfreq(klon, klev), density(klon, klev), sigma_t(klon, klev), & |
---|
79 | visc_mol(klon, klev), alt(klon, klev), sigsqmcw(klon, klev, nazmth), & |
---|
80 | sigmatm(klon, klev), ak_alpha(klon, nazmth), k_alpha(klon, nazmth), & |
---|
81 | mmin_alpha(klon, nazmth), i_alpha(klon, nazmth), rmswind(klon), & |
---|
82 | bvfbot(klon), densbot(klon) |
---|
83 | REAL smoothr1(klon, klev), smoothr2(klon, klev) |
---|
84 | REAL sigalpmc(klon, klev, nazmth) |
---|
85 | REAL f2mod(klon, klev) |
---|
86 | |
---|
87 | ! * THES ARE THE INPUT PARAMETERS FOR HINES ROUTINE AND |
---|
88 | ! * ARE SPECIFIED IN ROUTINE HINES_SETUP. SINCE THIS IS CALLED |
---|
89 | ! * ONLY AT FIRST CALL TO THIS ROUTINE THESE VARIABLES MUST BE SAVED |
---|
90 | ! * FOR USE AT SUBSEQUENT CALLS. THIS CAN BE AVOIDED BY CALLING |
---|
91 | ! * HINES_SETUP IN MAIN PROGRAM AND PASSING THE PARAMETERS AS |
---|
92 | ! * SUBROUTINE ARGUEMENTS. |
---|
93 | |
---|
94 | |
---|
95 | REAL rmscon |
---|
96 | INTEGER nmessg, iprint, ilrms |
---|
97 | INTEGER ifl |
---|
98 | |
---|
99 | INTEGER naz, icutoff, nsmax, iheatcal |
---|
100 | REAL slope, f1, f2, f3, f5, f6, kstar(klon), alt_cutoff, smco |
---|
101 | |
---|
102 | ! PROVIDED AS INPUT |
---|
103 | |
---|
104 | INTEGER nlon, nlev |
---|
105 | |
---|
106 | REAL dtime |
---|
107 | REAL paphm1x(nlon, nlev+1), papm1x(nlon, nlev) |
---|
108 | REAL ux(nlon, nlev), vx(nlon, nlev), tx(nlon, nlev) |
---|
109 | |
---|
110 | ! VARIABLES FOR OUTPUT |
---|
111 | |
---|
112 | |
---|
113 | REAL d_t_hin(nlon, nlev), d_u_hin(nlon, nlev), d_v_hin(nlon, nlev) |
---|
114 | REAL zustrhi(nlon), zvstrhi(nlon) |
---|
115 | |
---|
116 | |
---|
117 | ! * LOGICAL SWITCHES TO CONTROL PRECIP ENHANCEMENT AND |
---|
118 | ! * HINES' DOPPLER SPREADING EXTROWAVE GW DRAG. |
---|
119 | ! * LOZPR IS TRUE FOR ZPR ENHANCEMENT |
---|
120 | |
---|
121 | LOGICAL lozpr, lorms(klon) |
---|
122 | |
---|
123 | ! LOCAL PARAMETERS TO MAKE THINGS WORK (TEMPORARY VARIABLE) |
---|
124 | |
---|
125 | REAL rhoh2o, zpcons, rgocp, zlat, dttdsf, ratio, hscal |
---|
126 | INTEGER i, j, l, jl, jk, le, lref, lrefp, levbot |
---|
127 | |
---|
128 | ! DATA PARAMETERS NEEDED, EXPLAINED LATER |
---|
129 | |
---|
130 | REAL v0, vmin, dmpscal, taufac, hmin, apibt, cpart, fcrit |
---|
131 | REAL pcrit, pcons |
---|
132 | INTEGER iplev, ierror |
---|
133 | |
---|
134 | |
---|
135 | |
---|
136 | ! PRINT *,' IT IS STARTED HINES GOING ON...' |
---|
137 | |
---|
138 | |
---|
139 | |
---|
140 | |
---|
141 | ! * COMPUTATIONAL CONSTANTS. |
---|
142 | ! ------------- ---------- |
---|
143 | |
---|
144 | |
---|
145 | d_t_hin(:, :) = 0. |
---|
146 | |
---|
147 | rhoh2o = 1000. |
---|
148 | zpcons = (1000.*86400.)/rhoh2o |
---|
149 | ! ym KLEVM1=KLEV-1 |
---|
150 | |
---|
151 | |
---|
152 | DO jl = kidia, kfdia |
---|
153 | paphm1(jl, 1) = paphm1x(jl, klev+1) |
---|
154 | DO jk = 1, klev |
---|
155 | le = klev + 1 - jk |
---|
156 | paphm1(jl, jk+1) = paphm1x(jl, le) |
---|
157 | papm1(jl, jk) = papm1x(jl, le) |
---|
158 | ptm1(jl, jk) = tx(jl, le) |
---|
159 | pum1(jl, jk) = ux(jl, le) |
---|
160 | pvm1(jl, jk) = vx(jl, le) |
---|
161 | END DO |
---|
162 | END DO |
---|
163 | |
---|
164 | ! Define constants and arrays needed for the ccc/mam gwd scheme |
---|
165 | ! *Constants: |
---|
166 | |
---|
167 | rgocp = rd/rcpd |
---|
168 | lrefp = klev - 1 |
---|
169 | lref = klev - 2 |
---|
170 | ! 1 |
---|
171 | ! 1 *Arrays |
---|
172 | ! 1 |
---|
173 | DO jk = 1, klev |
---|
174 | DO jl = kidia, kfdia |
---|
175 | shj(jl, jk) = papm1(jl, jk)/paphm1(jl, klev+1) |
---|
176 | sgj(jl, jk) = papm1(jl, jk)/paphm1(jl, klev+1) |
---|
177 | dsgj(jl, jk) = (paphm1(jl,jk+1)-paphm1(jl,jk))/paphm1(jl, klev+1) |
---|
178 | shxkj(jl, jk) = (papm1(jl,jk)/paphm1(jl,klev+1))**rgocp |
---|
179 | th(jl, jk) = ptm1(jl, jk) |
---|
180 | END DO |
---|
181 | END DO |
---|
182 | |
---|
183 | ! C |
---|
184 | DO jl = kidia, kfdia |
---|
185 | pressg(jl) = paphm1(jl, klev+1) |
---|
186 | END DO |
---|
187 | |
---|
188 | |
---|
189 | DO jl = kidia, kfdia |
---|
190 | prflux(jl) = 0.0 |
---|
191 | zpr(jl) = zpcons*prflux(jl) |
---|
192 | zlat = (rlat(jl)/180.)*rpi |
---|
193 | coslat(jl) = cos(zlat) |
---|
194 | END DO |
---|
195 | |
---|
196 | ! /######################################################################### |
---|
197 | ! / |
---|
198 | ! / |
---|
199 | |
---|
200 | ! * AUG. 14/95 - C. MCLANDRESS. |
---|
201 | ! * SEP. 95 N. MCFARLANE. |
---|
202 | |
---|
203 | ! * THIS ROUTINE CALCULATES THE HORIZONTAL WIND TENDENCIES |
---|
204 | ! * DUE TO MCFARLANE'S OROGRAPHIC GW DRAG SCHEME, HINES' |
---|
205 | ! * DOPPLER SPREAD SCHEME FOR "EXTROWAVES" AND ADDS ON |
---|
206 | ! * ROOF DRAG. IT IS BASED ON THE ROUTINE GWDFLX8. |
---|
207 | |
---|
208 | ! * LREFP IS THE INDEX OF THE MODEL LEVEL BELOW THE REFERENCE LEVEL |
---|
209 | ! * I/O ARRAYS PASSED FROM MAIN. |
---|
210 | ! * (PRESSG = SURFACE PRESSURE) |
---|
211 | |
---|
212 | |
---|
213 | |
---|
214 | |
---|
215 | ! * CONSTANTS VALUES DEFINED IN DATA STATEMENT ARE : |
---|
216 | ! * VMIN = MIMINUM WIND IN THE DIRECTION OF REFERENCE LEVEL |
---|
217 | ! * WIND BEFORE WE CONSIDER BREAKING TO HAVE OCCURED. |
---|
218 | ! * DMPSCAL = DAMPING TIME FOR GW DRAG IN SECONDS. |
---|
219 | ! * TAUFAC = 1/(LENGTH SCALE). |
---|
220 | ! * HMIN = MIMINUM ENVELOPE HEIGHT REQUIRED TO PRODUCE GW DRAG. |
---|
221 | ! * V0 = VALUE OF WIND THAT APPROXIMATES ZERO. |
---|
222 | |
---|
223 | |
---|
224 | DATA vmin/5.0/, v0/1.E-10/, taufac/5.E-6/, hmin/40000./, dmpscal/6.5E+6/, & |
---|
225 | apibt/1.5708/, cpart/0.7/, fcrit/1./ |
---|
226 | |
---|
227 | ! * HINES EXTROWAVE GWD CONSTANTS DEFINED IN DATA STATEMENT ARE: |
---|
228 | ! * RMSCON = ROOT MEAN SQUARE GRAVITY WAVE WIND AT LOWEST LEVEL (M/S). |
---|
229 | ! * NMESSG = UNIT NUMBER FOR PRINTED MESSAGES. |
---|
230 | ! * IPRINT = 1 TO DO PRINT OUT SOME HINES ARRAYS. |
---|
231 | ! * IFL = FIRST CALL FLAG TO HINES_SETUP ("SAVE" IT) |
---|
232 | ! * PCRIT = CRITICAL VALUE OF ZPR (MM/D) |
---|
233 | ! * IPLEV = LEVEL OF APPLICATION OF PRCIT |
---|
234 | ! * PCONS = FACTOR OF ZPR ENHANCEMENT |
---|
235 | |
---|
236 | |
---|
237 | DATA pcrit/5./, pcons/4.75/ |
---|
238 | |
---|
239 | iplev = lrefp - 1 |
---|
240 | |
---|
241 | DATA rmscon/1.00/iprint/2/, nmessg/6/ |
---|
242 | DATA ifl/0/ |
---|
243 | |
---|
244 | lozpr = .FALSE. |
---|
245 | |
---|
246 | ! ----------------------------------------------------------------------- |
---|
247 | |
---|
248 | |
---|
249 | |
---|
250 | ! * SET ERROR FLAG |
---|
251 | |
---|
252 | ierror = 0 |
---|
253 | |
---|
254 | ! * SPECIFY VARIOUS PARAMETERS FOR HINES ROUTINE AT VERY FIRST CALL. |
---|
255 | ! * (NOTE THAT ARRAY K_ALPHA IS SPECIFIED SO MAKE SURE THAT |
---|
256 | ! * IT IS NOT OVERWRITTEN LATER ON). |
---|
257 | |
---|
258 | CALL hines_setup(naz, slope, f1, f2, f3, f5, f6, kstar, icutoff, & |
---|
259 | alt_cutoff, smco, nsmax, iheatcal, k_alpha, ierror, nmessg, klon, nazmth, & |
---|
260 | coslat) |
---|
261 | IF (ierror/=0) GO TO 999 |
---|
262 | |
---|
263 | ! * START GWD CALCULATIONS. |
---|
264 | |
---|
265 | lref = lrefp - 1 |
---|
266 | |
---|
267 | |
---|
268 | DO j = 1, nazmth |
---|
269 | DO l = 1, klev |
---|
270 | DO i = kidia, klon |
---|
271 | sigsqmcw(i, l, j) = 0. |
---|
272 | END DO |
---|
273 | END DO |
---|
274 | END DO |
---|
275 | |
---|
276 | |
---|
277 | |
---|
278 | ! * INITIALIZE NECESSARY ARRAYS. |
---|
279 | |
---|
280 | DO l = 1, klev |
---|
281 | DO i = kidia, kfdia |
---|
282 | utendgw(i, l) = 0. |
---|
283 | vtendgw(i, l) = 0. |
---|
284 | |
---|
285 | uhs(i, l) = 0. |
---|
286 | vhs(i, l) = 0. |
---|
287 | |
---|
288 | END DO |
---|
289 | END DO |
---|
290 | |
---|
291 | ! * IF USING HINES SCHEME THEN CALCULATE B V FREQUENCY AT ALL POINTS |
---|
292 | ! * AND SMOOTH BVFREQ. |
---|
293 | |
---|
294 | DO l = 2, klev |
---|
295 | DO i = kidia, kfdia |
---|
296 | dttdsf = (th(i,l)/shxkj(i,l)-th(i,l-1)/shxkj(i,l-1))/ & |
---|
297 | (shj(i,l)-shj(i,l-1)) |
---|
298 | dttdsf = min(dttdsf, -5./sgj(i,l)) |
---|
299 | bvfreq(i, l) = sqrt(-dttdsf*sgj(i,l)*(sgj(i,l)**rgocp)/rd)*rg/ptm1(i, l & |
---|
300 | ) |
---|
301 | END DO |
---|
302 | END DO |
---|
303 | DO l = 1, klev |
---|
304 | DO i = kidia, kfdia |
---|
305 | IF (l==1) THEN |
---|
306 | bvfreq(i, l) = bvfreq(i, l+1) |
---|
307 | END IF |
---|
308 | IF (l>1) THEN |
---|
309 | ratio = 5.*log(sgj(i,l)/sgj(i,l-1)) |
---|
310 | bvfreq(i, l) = (bvfreq(i,l-1)+ratio*bvfreq(i,l))/(1.+ratio) |
---|
311 | END IF |
---|
312 | END DO |
---|
313 | END DO |
---|
314 | |
---|
315 | ! * CALCULATE GW DRAG DUE TO HINES' EXTROWAVES |
---|
316 | ! * SET MOLECULAR VISCOSITY TO A VERY SMALL VALUE. |
---|
317 | ! * IF THE MODEL TOP IS GREATER THAN 100 KM THEN THE ACTUAL |
---|
318 | ! * VISCOSITY COEFFICIENT COULD BE SPECIFIED HERE. |
---|
319 | |
---|
320 | DO l = 1, klev |
---|
321 | DO i = kidia, kfdia |
---|
322 | visc_mol(i, l) = 1.5E-5 |
---|
323 | drag_u(i, l) = 0. |
---|
324 | drag_v(i, l) = 0. |
---|
325 | flux_u(i, l) = 0. |
---|
326 | flux_v(i, l) = 0. |
---|
327 | heat(i, l) = 0. |
---|
328 | diffco(i, l) = 0. |
---|
329 | END DO |
---|
330 | END DO |
---|
331 | |
---|
332 | ! * ALTITUDE AND DENSITY AT BOTTOM. |
---|
333 | |
---|
334 | DO i = kidia, kfdia |
---|
335 | hscal = rd*ptm1(i, klev)/rg |
---|
336 | density(i, klev) = sgj(i, klev)*pressg(i)/(rg*hscal) |
---|
337 | alt(i, klev) = 0. |
---|
338 | END DO |
---|
339 | |
---|
340 | ! * ALTITUDE AND DENSITY AT REMAINING LEVELS. |
---|
341 | |
---|
342 | DO l = klev - 1, 1, -1 |
---|
343 | DO i = kidia, kfdia |
---|
344 | hscal = rd*ptm1(i, l)/rg |
---|
345 | alt(i, l) = alt(i, l+1) + hscal*dsgj(i, l)/sgj(i, l) |
---|
346 | density(i, l) = sgj(i, l)*pressg(i)/(rg*hscal) |
---|
347 | END DO |
---|
348 | END DO |
---|
349 | |
---|
350 | |
---|
351 | ! * INITIALIZE SWITCHES FOR HINES GWD CALCULATION |
---|
352 | |
---|
353 | ilrms = 0 |
---|
354 | |
---|
355 | DO i = kidia, kfdia |
---|
356 | lorms(i) = .FALSE. |
---|
357 | END DO |
---|
358 | |
---|
359 | |
---|
360 | ! * DEFILE BOTTOM LAUNCH LEVEL |
---|
361 | |
---|
362 | levbot = iplev |
---|
363 | |
---|
364 | ! * BACKGROUND WIND MINUS VALUE AT BOTTOM LAUNCH LEVEL. |
---|
365 | |
---|
366 | DO l = 1, levbot |
---|
367 | DO i = kidia, kfdia |
---|
368 | uhs(i, l) = pum1(i, l) - pum1(i, levbot) |
---|
369 | vhs(i, l) = pvm1(i, l) - pvm1(i, levbot) |
---|
370 | END DO |
---|
371 | END DO |
---|
372 | |
---|
373 | ! * SPECIFY ROOT MEAN SQUARE WIND AT BOTTOM LAUNCH LEVEL. |
---|
374 | |
---|
375 | DO i = kidia, kfdia |
---|
376 | rmswind(i) = rmscon |
---|
377 | END DO |
---|
378 | |
---|
379 | IF (lozpr) THEN |
---|
380 | DO i = kidia, kfdia |
---|
381 | IF (zpr(i)>pcrit) THEN |
---|
382 | rmswind(i) = rmscon + ((zpr(i)-pcrit)/zpr(i))*pcons |
---|
383 | END IF |
---|
384 | END DO |
---|
385 | END IF |
---|
386 | |
---|
387 | DO i = kidia, kfdia |
---|
388 | IF (rmswind(i)>0.0) THEN |
---|
389 | ilrms = ilrms + 1 |
---|
390 | lorms(i) = .TRUE. |
---|
391 | END IF |
---|
392 | END DO |
---|
393 | |
---|
394 | ! * CALCULATE GWD (NOTE THAT DIFFUSION COEFFICIENT AND |
---|
395 | ! * HEATING RATE ONLY CALCULATED IF IHEATCAL = 1). |
---|
396 | |
---|
397 | IF (ilrms>0) THEN |
---|
398 | |
---|
399 | CALL hines_extro0(drag_u, drag_v, heat, diffco, flux_u, flux_v, uhs, vhs, & |
---|
400 | bvfreq, density, visc_mol, alt, rmswind, k_alpha, m_alpha, v_alpha, & |
---|
401 | sigma_alpha, sigsqh_alpha, ak_alpha, mmin_alpha, i_alpha, sigma_t, & |
---|
402 | densbot, bvfbot, 1, iheatcal, icutoff, iprint, nsmax, smco, alt_cutoff, & |
---|
403 | kstar, slope, f1, f2, f3, f5, f6, naz, sigsqmcw, sigmatm, kidia, klon, & |
---|
404 | 1, levbot, klon, klev, nazmth, lorms, smoothr1, smoothr2, sigalpmc, & |
---|
405 | f2mod) |
---|
406 | |
---|
407 | ! * ADD ON HINES' GWD TENDENCIES TO OROGRAPHIC TENDENCIES AND |
---|
408 | ! * APPLY HINES' GW DRAG ON (UROW,VROW) WORK ARRAYS. |
---|
409 | |
---|
410 | DO l = 1, klev |
---|
411 | DO i = kidia, kfdia |
---|
412 | utendgw(i, l) = utendgw(i, l) + drag_u(i, l) |
---|
413 | vtendgw(i, l) = vtendgw(i, l) + drag_v(i, l) |
---|
414 | END DO |
---|
415 | END DO |
---|
416 | |
---|
417 | |
---|
418 | ! * END OF HINES CALCULATIONS. |
---|
419 | |
---|
420 | END IF |
---|
421 | |
---|
422 | ! ----------------------------------------------------------------------- |
---|
423 | |
---|
424 | DO jl = kidia, kfdia |
---|
425 | zustrhi(jl) = flux_u(jl, 1) |
---|
426 | zvstrhi(jl) = flux_v(jl, 1) |
---|
427 | DO jk = 1, klev |
---|
428 | le = klev - jk + 1 |
---|
429 | d_u_hin(jl, jk) = utendgw(jl, le)*dtime |
---|
430 | d_v_hin(jl, jk) = vtendgw(jl, le)*dtime |
---|
431 | END DO |
---|
432 | END DO |
---|
433 | |
---|
434 | ! PRINT *,'UTENDGW:',UTENDGW |
---|
435 | |
---|
436 | ! PRINT *,' HINES HAS BEEN COMPLETED (LONG ISNT IT...)' |
---|
437 | |
---|
438 | RETURN |
---|
439 | 999 CONTINUE |
---|
440 | |
---|
441 | ! * IF ERROR DETECTED THEN ABORT. |
---|
442 | |
---|
443 | WRITE (nmessg, 6000) |
---|
444 | WRITE (nmessg, 6010) ierror |
---|
445 | 6000 FORMAT (/' EXECUTION ABORTED IN GWDOREXV') |
---|
446 | 6010 FORMAT (' ERROR FLAG =', I4) |
---|
447 | |
---|
448 | |
---|
449 | RETURN |
---|
450 | END SUBROUTINE hines_gwd |
---|
451 | ! / |
---|
452 | ! / |
---|
453 | |
---|
454 | |
---|
455 | SUBROUTINE hines_extro0(drag_u, drag_v, heat, diffco, flux_u, flux_v, vel_u, & |
---|
456 | vel_v, bvfreq, density, visc_mol, alt, rmswind, k_alpha, m_alpha, & |
---|
457 | v_alpha, sigma_alpha, sigsqh_alpha, ak_alpha, mmin_alpha, i_alpha, & |
---|
458 | sigma_t, densb, bvfb, iorder, iheatcal, icutoff, iprint, nsmax, smco, & |
---|
459 | alt_cutoff, kstar, slope, f1, f2, f3, f5, f6, naz, sigsqmcw, sigmatm, & |
---|
460 | il1, il2, lev1, lev2, nlons, nlevs, nazmth, lorms, smoothr1, smoothr2, & |
---|
461 | sigalpmc, f2mod) |
---|
462 | |
---|
463 | IMPLICIT NONE |
---|
464 | |
---|
465 | ! Main routine for Hines' "extrowave" gravity wave parameterization based |
---|
466 | ! on Hines' Doppler spread theory. This routine calculates zonal |
---|
467 | ! and meridional components of gravity wave drag, heating rates |
---|
468 | ! and diffusion coefficient on a longitude by altitude grid. |
---|
469 | ! No "mythical" lower boundary region calculation is made so it |
---|
470 | ! is assumed that lowest level winds are weak (i.e, approximately zero). |
---|
471 | |
---|
472 | ! Aug. 13/95 - C. McLandress |
---|
473 | ! SEPT. /95 - N.McFarlane |
---|
474 | |
---|
475 | ! Modifications: |
---|
476 | |
---|
477 | ! Output arguements: |
---|
478 | |
---|
479 | ! * DRAG_U = zonal component of gravity wave drag (m/s^2). |
---|
480 | ! * DRAG_V = meridional component of gravity wave drag (m/s^2). |
---|
481 | ! * HEAT = gravity wave heating (K/sec). |
---|
482 | ! * DIFFCO = diffusion coefficient (m^2/sec) |
---|
483 | ! * FLUX_U = zonal component of vertical momentum flux (Pascals) |
---|
484 | ! * FLUX_V = meridional component of vertical momentum flux (Pascals) |
---|
485 | |
---|
486 | ! Input arguements: |
---|
487 | |
---|
488 | ! * VEL_U = background zonal wind component (m/s). |
---|
489 | ! * VEL_V = background meridional wind component (m/s). |
---|
490 | ! * BVFREQ = background Brunt Vassala frequency (radians/sec). |
---|
491 | ! * DENSITY = background density (kg/m^3) |
---|
492 | ! * VISC_MOL = molecular viscosity (m^2/s) |
---|
493 | ! * ALT = altitude of momentum, density, buoyancy levels (m) |
---|
494 | ! * (NOTE: levels ordered so that ALT(I,1) > ALT(I,2), etc.) |
---|
495 | ! * RMSWIND = root mean square gravity wave wind at lowest level (m/s). |
---|
496 | ! * K_ALPHA = horizontal wavenumber of each azimuth (1/m). |
---|
497 | ! * IORDER = 1 means vertical levels are indexed from top down |
---|
498 | ! * (i.e., highest level indexed 1 and lowest level NLEVS); |
---|
499 | ! * .NE. 1 highest level is index NLEVS. |
---|
500 | ! * IHEATCAL = 1 to calculate heating rates and diffusion coefficient. |
---|
501 | ! * IPRINT = 1 to print out various arrays. |
---|
502 | ! * ICUTOFF = 1 to exponentially damp GWD, heating and diffusion |
---|
503 | ! * arrays above ALT_CUTOFF; otherwise arrays not modified. |
---|
504 | ! * ALT_CUTOFF = altitude in meters above which exponential decay applied. |
---|
505 | ! * SMCO = smoothing factor used to smooth cutoff vertical |
---|
506 | ! * wavenumbers and total rms winds in vertical direction |
---|
507 | ! * before calculating drag or heating |
---|
508 | ! * (SMCO >= 1 ==> 1:SMCO:1 stencil used). |
---|
509 | ! * NSMAX = number of times smoother applied ( >= 1), |
---|
510 | ! * = 0 means no smoothing performed. |
---|
511 | ! * KSTAR = typical gravity wave horizontal wavenumber (1/m). |
---|
512 | ! * SLOPE = slope of incident vertical wavenumber spectrum |
---|
513 | ! * (SLOPE must equal 1., 1.5 or 2.). |
---|
514 | ! * F1 to F6 = Hines's fudge factors (F4 not needed since used for |
---|
515 | ! * vertical flux of vertical momentum). |
---|
516 | ! * NAZ = actual number of horizontal azimuths used. |
---|
517 | ! * IL1 = first longitudinal index to use (IL1 >= 1). |
---|
518 | ! * IL2 = last longitudinal index to use (IL1 <= IL2 <= NLONS). |
---|
519 | ! * LEV1 = index of first level for drag calculation. |
---|
520 | ! * LEV2 = index of last level for drag calculation |
---|
521 | ! * (i.e., LEV1 < LEV2 <= NLEVS). |
---|
522 | ! * NLONS = number of longitudes. |
---|
523 | ! * NLEVS = number of vertical levels. |
---|
524 | ! * NAZMTH = azimuthal array dimension (NAZMTH >= NAZ). |
---|
525 | |
---|
526 | ! Work arrays. |
---|
527 | |
---|
528 | ! * M_ALPHA = cutoff vertical wavenumber (1/m). |
---|
529 | ! * V_ALPHA = wind component at each azimuth (m/s) and if IHEATCAL=1 |
---|
530 | ! * holds vertical derivative of cutoff wavenumber. |
---|
531 | ! * SIGMA_ALPHA = total rms wind in each azimuth (m/s). |
---|
532 | ! * SIGSQH_ALPHA = portion of wind variance from waves having wave |
---|
533 | ! * normals in the alpha azimuth (m/s). |
---|
534 | ! * SIGMA_T = total rms horizontal wind (m/s). |
---|
535 | ! * AK_ALPHA = spectral amplitude factor at each azimuth |
---|
536 | ! * (i.e.,{AjKj}) in m^4/s^2. |
---|
537 | ! * I_ALPHA = Hines' integral. |
---|
538 | ! * MMIN_ALPHA = minimum value of cutoff wavenumber. |
---|
539 | ! * DENSB = background density at bottom level. |
---|
540 | ! * BVFB = buoyancy frequency at bottom level and |
---|
541 | ! * work array for ICUTOFF = 1. |
---|
542 | |
---|
543 | ! * LORMS = .TRUE. for drag computation |
---|
544 | |
---|
545 | INTEGER naz, nlons, nlevs, nazmth, il1, il2, lev1, lev2 |
---|
546 | INTEGER icutoff, nsmax, iorder, iheatcal, iprint |
---|
547 | REAL kstar(nlons), f1, f2, f3, f5, f6, slope |
---|
548 | REAL alt_cutoff, smco |
---|
549 | REAL drag_u(nlons, nlevs), drag_v(nlons, nlevs) |
---|
550 | REAL heat(nlons, nlevs), diffco(nlons, nlevs) |
---|
551 | REAL flux_u(nlons, nlevs), flux_v(nlons, nlevs) |
---|
552 | REAL vel_u(nlons, nlevs), vel_v(nlons, nlevs) |
---|
553 | REAL bvfreq(nlons, nlevs), density(nlons, nlevs) |
---|
554 | REAL visc_mol(nlons, nlevs), alt(nlons, nlevs) |
---|
555 | REAL rmswind(nlons), bvfb(nlons), densb(nlons) |
---|
556 | REAL sigma_t(nlons, nlevs), sigsqmcw(nlons, nlevs, nazmth) |
---|
557 | REAL sigma_alpha(nlons, nlevs, nazmth), sigmatm(nlons, nlevs) |
---|
558 | REAL sigsqh_alpha(nlons, nlevs, nazmth) |
---|
559 | REAL m_alpha(nlons, nlevs, nazmth), v_alpha(nlons, nlevs, nazmth) |
---|
560 | REAL ak_alpha(nlons, nazmth), k_alpha(nlons, nazmth) |
---|
561 | REAL mmin_alpha(nlons, nazmth), i_alpha(nlons, nazmth) |
---|
562 | REAL smoothr1(nlons, nlevs), smoothr2(nlons, nlevs) |
---|
563 | REAL sigalpmc(nlons, nlevs, nazmth) |
---|
564 | REAL f2mod(nlons, nlevs) |
---|
565 | |
---|
566 | LOGICAL lorms(nlons) |
---|
567 | |
---|
568 | ! Internal variables. |
---|
569 | |
---|
570 | INTEGER levbot, levtop, i, n, l, lev1p, lev2m |
---|
571 | INTEGER ilprt1, ilprt2 |
---|
572 | ! ----------------------------------------------------------------------- |
---|
573 | |
---|
574 | ! PRINT *,' IN HINES_EXTRO0' |
---|
575 | lev1p = lev1 + 1 |
---|
576 | lev2m = lev2 - 1 |
---|
577 | |
---|
578 | ! Index of lowest altitude level (bottom of drag calculation). |
---|
579 | |
---|
580 | levbot = lev2 |
---|
581 | levtop = lev1 |
---|
582 | IF (iorder/=1) THEN |
---|
583 | WRITE (6, 1) |
---|
584 | 1 FORMAT (2X, ' error: IORDER NOT ONE! ') |
---|
585 | END IF |
---|
586 | |
---|
587 | ! Buoyancy and density at bottom level. |
---|
588 | |
---|
589 | DO i = il1, il2 |
---|
590 | bvfb(i) = bvfreq(i, levbot) |
---|
591 | densb(i) = density(i, levbot) |
---|
592 | END DO |
---|
593 | |
---|
594 | ! initialize some variables |
---|
595 | |
---|
596 | DO n = 1, naz |
---|
597 | DO l = lev1, lev2 |
---|
598 | DO i = il1, il2 |
---|
599 | m_alpha(i, l, n) = 0.0 |
---|
600 | END DO |
---|
601 | END DO |
---|
602 | END DO |
---|
603 | DO l = lev1, lev2 |
---|
604 | DO i = il1, il2 |
---|
605 | sigma_t(i, l) = 0.0 |
---|
606 | END DO |
---|
607 | END DO |
---|
608 | DO n = 1, naz |
---|
609 | DO i = il1, il2 |
---|
610 | i_alpha(i, n) = 0.0 |
---|
611 | END DO |
---|
612 | END DO |
---|
613 | |
---|
614 | ! Compute azimuthal wind components from zonal and meridional winds. |
---|
615 | |
---|
616 | CALL hines_wind(v_alpha, vel_u, vel_v, naz, il1, il2, lev1, lev2, nlons, & |
---|
617 | nlevs, nazmth) |
---|
618 | |
---|
619 | ! Calculate cutoff vertical wavenumber and velocity variances. |
---|
620 | |
---|
621 | CALL hines_wavnum(m_alpha, sigma_alpha, sigsqh_alpha, sigma_t, ak_alpha, & |
---|
622 | v_alpha, visc_mol, density, densb, bvfreq, bvfb, rmswind, i_alpha, & |
---|
623 | mmin_alpha, kstar, slope, f1, f2, f3, naz, levbot, levtop, il1, il2, & |
---|
624 | nlons, nlevs, nazmth, sigsqmcw, sigmatm, lorms, sigalpmc, f2mod) |
---|
625 | ! Smooth cutoff wavenumbers and total rms velocity in the vertical |
---|
626 | ! direction NSMAX times, using FLUX_U as temporary work array. |
---|
627 | |
---|
628 | IF (nsmax>0) THEN |
---|
629 | DO n = 1, naz |
---|
630 | DO l = lev1, lev2 |
---|
631 | DO i = il1, il2 |
---|
632 | smoothr1(i, l) = m_alpha(i, l, n) |
---|
633 | END DO |
---|
634 | END DO |
---|
635 | CALL vert_smooth(smoothr1, smoothr2, smco, nsmax, il1, il2, lev1, lev2, & |
---|
636 | nlons, nlevs) |
---|
637 | DO l = lev1, lev2 |
---|
638 | DO i = il1, il2 |
---|
639 | m_alpha(i, l, n) = smoothr1(i, l) |
---|
640 | END DO |
---|
641 | END DO |
---|
642 | END DO |
---|
643 | CALL vert_smooth(sigma_t, smoothr2, smco, nsmax, il1, il2, lev1, lev2, & |
---|
644 | nlons, nlevs) |
---|
645 | END IF |
---|
646 | |
---|
647 | ! Calculate zonal and meridional components of the |
---|
648 | ! momentum flux and drag. |
---|
649 | |
---|
650 | CALL hines_flux(flux_u, flux_v, drag_u, drag_v, alt, density, densb, & |
---|
651 | m_alpha, ak_alpha, k_alpha, slope, naz, il1, il2, lev1, lev2, nlons, & |
---|
652 | nlevs, nazmth, lorms) |
---|
653 | |
---|
654 | ! Cutoff drag above ALT_CUTOFF, using BVFB as temporary work array. |
---|
655 | |
---|
656 | IF (icutoff==1) THEN |
---|
657 | CALL hines_exp(drag_u, bvfb, alt, alt_cutoff, iorder, il1, il2, lev1, & |
---|
658 | lev2, nlons, nlevs) |
---|
659 | CALL hines_exp(drag_v, bvfb, alt, alt_cutoff, iorder, il1, il2, lev1, & |
---|
660 | lev2, nlons, nlevs) |
---|
661 | END IF |
---|
662 | |
---|
663 | ! Print out various arrays for diagnostic purposes. |
---|
664 | |
---|
665 | IF (iprint==1) THEN |
---|
666 | ilprt1 = 15 |
---|
667 | ilprt2 = 16 |
---|
668 | CALL hines_print(flux_u, flux_v, drag_u, drag_v, alt, sigma_t, & |
---|
669 | sigma_alpha, v_alpha, m_alpha, 1, 1, 6, ilprt1, ilprt2, lev1, lev2, & |
---|
670 | naz, nlons, nlevs, nazmth) |
---|
671 | END IF |
---|
672 | |
---|
673 | ! If not calculating heating rate and diffusion coefficient then finished. |
---|
674 | |
---|
675 | IF (iheatcal/=1) RETURN |
---|
676 | |
---|
677 | ! Calculate vertical derivative of cutoff wavenumber (store |
---|
678 | ! in array V_ALPHA) using centered differences at interior gridpoints |
---|
679 | ! and one-sided differences at first and last levels. |
---|
680 | |
---|
681 | DO n = 1, naz |
---|
682 | DO l = lev1p, lev2m |
---|
683 | DO i = il1, il2 |
---|
684 | v_alpha(i, l, n) = (m_alpha(i,l+1,n)-m_alpha(i,l-1,n))/ & |
---|
685 | (alt(i,l+1)-alt(i,l-1)) |
---|
686 | END DO |
---|
687 | END DO |
---|
688 | DO i = il1, il2 |
---|
689 | v_alpha(i, lev1, n) = (m_alpha(i,lev1p,n)-m_alpha(i,lev1,n))/ & |
---|
690 | (alt(i,lev1p)-alt(i,lev1)) |
---|
691 | END DO |
---|
692 | DO i = il1, il2 |
---|
693 | v_alpha(i, lev2, n) = (m_alpha(i,lev2,n)-m_alpha(i,lev2m,n))/ & |
---|
694 | (alt(i,lev2)-alt(i,lev2m)) |
---|
695 | END DO |
---|
696 | END DO |
---|
697 | |
---|
698 | ! Heating rate and diffusion coefficient. |
---|
699 | |
---|
700 | CALL hines_heat(heat, diffco, m_alpha, v_alpha, ak_alpha, k_alpha, bvfreq, & |
---|
701 | density, densb, sigma_t, visc_mol, kstar, slope, f2, f3, f5, f6, naz, & |
---|
702 | il1, il2, lev1, lev2, nlons, nlevs, nazmth) |
---|
703 | |
---|
704 | ! Finished. |
---|
705 | |
---|
706 | RETURN |
---|
707 | ! ----------------------------------------------------------------------- |
---|
708 | END SUBROUTINE hines_extro0 |
---|
709 | |
---|
710 | SUBROUTINE hines_wavnum(m_alpha, sigma_alpha, sigsqh_alpha, sigma_t, & |
---|
711 | ak_alpha, v_alpha, visc_mol, density, densb, bvfreq, bvfb, rms_wind, & |
---|
712 | i_alpha, mmin_alpha, kstar, slope, f1, f2, f3, naz, levbot, levtop, il1, & |
---|
713 | il2, nlons, nlevs, nazmth, sigsqmcw, sigmatm, lorms, sigalpmc, f2mod) |
---|
714 | IMPLICIT NONE |
---|
715 | ! This routine calculates the cutoff vertical wavenumber and velocity |
---|
716 | ! variances on a longitude by altitude grid for the Hines' Doppler |
---|
717 | ! spread gravity wave drag parameterization scheme. |
---|
718 | ! NOTE: (1) only values of four or eight can be used for # azimuths (NAZ). |
---|
719 | ! (2) only values of 1.0, 1.5 or 2.0 can be used for slope (SLOPE). |
---|
720 | |
---|
721 | ! Aug. 10/95 - C. McLandress |
---|
722 | |
---|
723 | ! Output arguements: |
---|
724 | |
---|
725 | ! * M_ALPHA = cutoff wavenumber at each azimuth (1/m). |
---|
726 | ! * SIGMA_ALPHA = total rms wind in each azimuth (m/s). |
---|
727 | ! * SIGSQH_ALPHA = portion of wind variance from waves having wave |
---|
728 | ! * normals in the alpha azimuth (m/s). |
---|
729 | ! * SIGMA_T = total rms horizontal wind (m/s). |
---|
730 | ! * AK_ALPHA = spectral amplitude factor at each azimuth |
---|
731 | ! * (i.e.,{AjKj}) in m^4/s^2. |
---|
732 | |
---|
733 | ! Input arguements: |
---|
734 | |
---|
735 | ! * V_ALPHA = wind component at each azimuth (m/s). |
---|
736 | ! * VISC_MOL = molecular viscosity (m^2/s) |
---|
737 | ! * DENSITY = background density (kg/m^3). |
---|
738 | ! * DENSB = background density at model bottom (kg/m^3). |
---|
739 | ! * BVFREQ = background Brunt Vassala frequency (radians/sec). |
---|
740 | ! * BVFB = background Brunt Vassala frequency at model bottom. |
---|
741 | ! * RMS_WIND = root mean square gravity wave wind at lowest level (m/s). |
---|
742 | ! * KSTAR = typical gravity wave horizontal wavenumber (1/m). |
---|
743 | ! * SLOPE = slope of incident vertical wavenumber spectrum |
---|
744 | ! * (SLOPE = 1., 1.5 or 2.). |
---|
745 | ! * F1,F2,F3 = Hines's fudge factors. |
---|
746 | ! * NAZ = actual number of horizontal azimuths used (4 or 8). |
---|
747 | ! * LEVBOT = index of lowest vertical level. |
---|
748 | ! * LEVTOP = index of highest vertical level |
---|
749 | ! * (NOTE: if LEVTOP < LEVBOT then level index |
---|
750 | ! * increases from top down). |
---|
751 | ! * IL1 = first longitudinal index to use (IL1 >= 1). |
---|
752 | ! * IL2 = last longitudinal index to use (IL1 <= IL2 <= NLONS). |
---|
753 | ! * NLONS = number of longitudes. |
---|
754 | ! * NLEVS = number of vertical levels. |
---|
755 | ! * NAZMTH = azimuthal array dimension (NAZMTH >= NAZ). |
---|
756 | |
---|
757 | ! * LORMS = .TRUE. for drag computation |
---|
758 | |
---|
759 | ! Input work arrays: |
---|
760 | |
---|
761 | ! * I_ALPHA = Hines' integral at a single level. |
---|
762 | ! * MMIN_ALPHA = minimum value of cutoff wavenumber. |
---|
763 | |
---|
764 | INTEGER naz, levbot, levtop, il1, il2, nlons, nlevs, nazmth |
---|
765 | REAL slope, kstar(nlons), f1, f2, f3, f2mfac |
---|
766 | REAL m_alpha(nlons, nlevs, nazmth) |
---|
767 | REAL sigma_alpha(nlons, nlevs, nazmth) |
---|
768 | REAL sigalpmc(nlons, nlevs, nazmth) |
---|
769 | REAL sigsqh_alpha(nlons, nlevs, nazmth) |
---|
770 | REAL sigsqmcw(nlons, nlevs, nazmth) |
---|
771 | REAL sigma_t(nlons, nlevs) |
---|
772 | REAL sigmatm(nlons, nlevs) |
---|
773 | REAL ak_alpha(nlons, nazmth) |
---|
774 | REAL v_alpha(nlons, nlevs, nazmth) |
---|
775 | REAL visc_mol(nlons, nlevs) |
---|
776 | REAL f2mod(nlons, nlevs) |
---|
777 | REAL density(nlons, nlevs), densb(nlons) |
---|
778 | REAL bvfreq(nlons, nlevs), bvfb(nlons), rms_wind(nlons) |
---|
779 | REAL i_alpha(nlons, nazmth), mmin_alpha(nlons, nazmth) |
---|
780 | |
---|
781 | LOGICAL lorms(nlons) |
---|
782 | |
---|
783 | ! Internal variables. |
---|
784 | |
---|
785 | INTEGER i, l, n, lstart, lend, lincr, lbelow |
---|
786 | REAL m_sub_m_turb, m_sub_m_mol, m_trial |
---|
787 | REAL visc, visc_min, azfac, sp1 |
---|
788 | |
---|
789 | ! c REAL N_OVER_M(1000), SIGFAC(1000) |
---|
790 | |
---|
791 | REAL n_over_m(nlons), sigfac(nlons) |
---|
792 | DATA visc_min/1.E-10/ |
---|
793 | ! ----------------------------------------------------------------------- |
---|
794 | |
---|
795 | |
---|
796 | ! PRINT *,'IN HINES_WAVNUM' |
---|
797 | sp1 = slope + 1. |
---|
798 | |
---|
799 | ! Indices of levels to process. |
---|
800 | |
---|
801 | IF (levbot>levtop) THEN |
---|
802 | lstart = levbot - 1 |
---|
803 | lend = levtop |
---|
804 | lincr = -1 |
---|
805 | ELSE |
---|
806 | WRITE (6, 1) |
---|
807 | 1 FORMAT (2X, ' error: IORDER NOT ONE! ') |
---|
808 | END IF |
---|
809 | |
---|
810 | ! Use horizontal isotropy to calculate azimuthal variances at bottom level. |
---|
811 | |
---|
812 | azfac = 1./real(naz) |
---|
813 | DO n = 1, naz |
---|
814 | DO i = il1, il2 |
---|
815 | sigsqh_alpha(i, levbot, n) = azfac*rms_wind(i)**2 |
---|
816 | END DO |
---|
817 | END DO |
---|
818 | |
---|
819 | ! Velocity variances at bottom level. |
---|
820 | |
---|
821 | CALL hines_sigma(sigma_t, sigma_alpha, sigsqh_alpha, naz, levbot, il1, il2, & |
---|
822 | nlons, nlevs, nazmth) |
---|
823 | |
---|
824 | CALL hines_sigma(sigmatm, sigalpmc, sigsqmcw, naz, levbot, il1, il2, nlons, & |
---|
825 | nlevs, nazmth) |
---|
826 | |
---|
827 | ! Calculate cutoff wavenumber and spectral amplitude factor |
---|
828 | ! at bottom level where it is assumed that background winds vanish |
---|
829 | ! and also initialize minimum value of cutoff wavnumber. |
---|
830 | |
---|
831 | DO n = 1, naz |
---|
832 | DO i = il1, il2 |
---|
833 | IF (lorms(i)) THEN |
---|
834 | m_alpha(i, levbot, n) = bvfb(i)/(f1*sigma_alpha(i,levbot,n)+f2* & |
---|
835 | sigma_t(i,levbot)) |
---|
836 | ak_alpha(i, n) = sigsqh_alpha(i, levbot, n)/ & |
---|
837 | (m_alpha(i,levbot,n)**sp1/sp1) |
---|
838 | mmin_alpha(i, n) = m_alpha(i, levbot, n) |
---|
839 | END IF |
---|
840 | END DO |
---|
841 | END DO |
---|
842 | |
---|
843 | ! Calculate quantities from the bottom upwards, |
---|
844 | ! starting one level above bottom. |
---|
845 | |
---|
846 | DO l = lstart, lend, lincr |
---|
847 | |
---|
848 | ! Level beneath present level. |
---|
849 | |
---|
850 | lbelow = l - lincr |
---|
851 | |
---|
852 | ! Calculate N/m_M where m_M is maximum permissible value of the vertical |
---|
853 | ! wavenumber (i.e., m > m_M are obliterated) and N is buoyancy frequency. |
---|
854 | ! m_M is taken as the smaller of the instability-induced |
---|
855 | ! wavenumber (M_SUB_M_TURB) and that imposed by molecular viscosity |
---|
856 | ! (M_SUB_M_MOL). Since variance at this level is not yet known |
---|
857 | ! use value at level below. |
---|
858 | |
---|
859 | DO i = il1, il2 |
---|
860 | IF (lorms(i)) THEN |
---|
861 | |
---|
862 | f2mfac = sigmatm(i, lbelow)**2 |
---|
863 | f2mod(i, lbelow) = 1. + 2.*f2mfac/(f2mfac+sigma_t(i,lbelow)**2) |
---|
864 | |
---|
865 | visc = amax1(visc_mol(i,l), visc_min) |
---|
866 | m_sub_m_turb = bvfreq(i, l)/(f2*f2mod(i,lbelow)*sigma_t(i,lbelow)) |
---|
867 | m_sub_m_mol = (bvfreq(i,l)*kstar(i)/visc)**0.33333333/f3 |
---|
868 | IF (m_sub_m_turb<m_sub_m_mol) THEN |
---|
869 | n_over_m(i) = f2*f2mod(i, lbelow)*sigma_t(i, lbelow) |
---|
870 | ELSE |
---|
871 | n_over_m(i) = bvfreq(i, l)/m_sub_m_mol |
---|
872 | END IF |
---|
873 | END IF |
---|
874 | END DO |
---|
875 | |
---|
876 | ! Calculate cutoff wavenumber at this level. |
---|
877 | |
---|
878 | DO n = 1, naz |
---|
879 | DO i = il1, il2 |
---|
880 | IF (lorms(i)) THEN |
---|
881 | |
---|
882 | ! Calculate trial value (since variance at this level is not yet |
---|
883 | ! known |
---|
884 | ! use value at level below). If trial value is negative or if it |
---|
885 | ! exceeds |
---|
886 | ! minimum value (not permitted) then set it to minimum value. |
---|
887 | |
---|
888 | m_trial = bvfb(i)/(f1*(sigma_alpha(i,lbelow,n)+sigalpmc(i,lbelow, & |
---|
889 | n))+n_over_m(i)+v_alpha(i,l,n)) |
---|
890 | IF (m_trial<=0. .OR. m_trial>mmin_alpha(i,n)) THEN |
---|
891 | m_trial = mmin_alpha(i, n) |
---|
892 | END IF |
---|
893 | m_alpha(i, l, n) = m_trial |
---|
894 | |
---|
895 | ! Reset minimum value of cutoff wavenumber if necessary. |
---|
896 | |
---|
897 | IF (m_alpha(i,l,n)<mmin_alpha(i,n)) THEN |
---|
898 | mmin_alpha(i, n) = m_alpha(i, l, n) |
---|
899 | END IF |
---|
900 | |
---|
901 | END IF |
---|
902 | END DO |
---|
903 | END DO |
---|
904 | |
---|
905 | ! Calculate the Hines integral at this level. |
---|
906 | |
---|
907 | CALL hines_intgrl(i_alpha, v_alpha, m_alpha, bvfb, slope, naz, l, il1, & |
---|
908 | il2, nlons, nlevs, nazmth, lorms) |
---|
909 | |
---|
910 | |
---|
911 | ! Calculate the velocity variances at this level. |
---|
912 | |
---|
913 | DO i = il1, il2 |
---|
914 | sigfac(i) = densb(i)/density(i, l)*bvfreq(i, l)/bvfb(i) |
---|
915 | END DO |
---|
916 | DO n = 1, naz |
---|
917 | DO i = il1, il2 |
---|
918 | sigsqh_alpha(i, l, n) = sigfac(i)*ak_alpha(i, n)*i_alpha(i, n) |
---|
919 | END DO |
---|
920 | END DO |
---|
921 | CALL hines_sigma(sigma_t, sigma_alpha, sigsqh_alpha, naz, l, il1, il2, & |
---|
922 | nlons, nlevs, nazmth) |
---|
923 | |
---|
924 | CALL hines_sigma(sigmatm, sigalpmc, sigsqmcw, naz, l, il1, il2, nlons, & |
---|
925 | nlevs, nazmth) |
---|
926 | |
---|
927 | ! End of level loop. |
---|
928 | |
---|
929 | END DO |
---|
930 | |
---|
931 | RETURN |
---|
932 | ! ----------------------------------------------------------------------- |
---|
933 | END SUBROUTINE hines_wavnum |
---|
934 | |
---|
935 | SUBROUTINE hines_wind(v_alpha, vel_u, vel_v, naz, il1, il2, lev1, lev2, & |
---|
936 | nlons, nlevs, nazmth) |
---|
937 | IMPLICIT NONE |
---|
938 | ! This routine calculates the azimuthal horizontal background wind |
---|
939 | ! components |
---|
940 | ! on a longitude by altitude grid for the case of 4 or 8 azimuths for |
---|
941 | ! the Hines' Doppler spread GWD parameterization scheme. |
---|
942 | |
---|
943 | ! Aug. 7/95 - C. McLandress |
---|
944 | |
---|
945 | ! Output arguement: |
---|
946 | |
---|
947 | ! * V_ALPHA = background wind component at each azimuth (m/s). |
---|
948 | ! * (note: first azimuth is in eastward direction |
---|
949 | ! * and rotate in counterclockwise direction.) |
---|
950 | |
---|
951 | ! Input arguements: |
---|
952 | |
---|
953 | ! * VEL_U = background zonal wind component (m/s). |
---|
954 | ! * VEL_V = background meridional wind component (m/s). |
---|
955 | ! * NAZ = actual number of horizontal azimuths used (must be 4 or 8). |
---|
956 | ! * IL1 = first longitudinal index to use (IL1 >= 1). |
---|
957 | ! * IL2 = last longitudinal index to use (IL1 <= IL2 <= NLONS). |
---|
958 | ! * LEV1 = first altitude level to use (LEV1 >=1). |
---|
959 | ! * LEV2 = last altitude level to use (LEV1 < LEV2 <= NLEVS). |
---|
960 | ! * NLONS = number of longitudes. |
---|
961 | ! * NLEVS = number of vertical levels. |
---|
962 | ! * NAZMTH = azimuthal array dimension (NAZMTH >= NAZ). |
---|
963 | |
---|
964 | ! Constants in DATA statements. |
---|
965 | |
---|
966 | ! * COS45 = cosine of 45 degrees. |
---|
967 | ! * UMIN = minimum allowable value for zonal or meridional |
---|
968 | ! * wind component (m/s). |
---|
969 | |
---|
970 | ! Subroutine arguements. |
---|
971 | |
---|
972 | INTEGER naz, il1, il2, lev1, lev2 |
---|
973 | INTEGER nlons, nlevs, nazmth |
---|
974 | REAL v_alpha(nlons, nlevs, nazmth) |
---|
975 | REAL vel_u(nlons, nlevs), vel_v(nlons, nlevs) |
---|
976 | |
---|
977 | ! Internal variables. |
---|
978 | |
---|
979 | INTEGER i, l |
---|
980 | REAL u, v, cos45, umin |
---|
981 | |
---|
982 | DATA cos45/0.7071068/ |
---|
983 | DATA umin/0.001/ |
---|
984 | ! ----------------------------------------------------------------------- |
---|
985 | |
---|
986 | ! Case with 4 azimuths. |
---|
987 | |
---|
988 | |
---|
989 | ! PRINT *,'IN HINES_WIND' |
---|
990 | IF (naz==4) THEN |
---|
991 | DO l = lev1, lev2 |
---|
992 | DO i = il1, il2 |
---|
993 | u = vel_u(i, l) |
---|
994 | v = vel_v(i, l) |
---|
995 | IF (abs(u)<umin) u = umin |
---|
996 | IF (abs(v)<umin) v = umin |
---|
997 | v_alpha(i, l, 1) = u |
---|
998 | v_alpha(i, l, 2) = v |
---|
999 | v_alpha(i, l, 3) = -u |
---|
1000 | v_alpha(i, l, 4) = -v |
---|
1001 | END DO |
---|
1002 | END DO |
---|
1003 | END IF |
---|
1004 | |
---|
1005 | ! Case with 8 azimuths. |
---|
1006 | |
---|
1007 | IF (naz==8) THEN |
---|
1008 | DO l = lev1, lev2 |
---|
1009 | DO i = il1, il2 |
---|
1010 | u = vel_u(i, l) |
---|
1011 | v = vel_v(i, l) |
---|
1012 | IF (abs(u)<umin) u = umin |
---|
1013 | IF (abs(v)<umin) v = umin |
---|
1014 | v_alpha(i, l, 1) = u |
---|
1015 | v_alpha(i, l, 2) = cos45*(v+u) |
---|
1016 | v_alpha(i, l, 3) = v |
---|
1017 | v_alpha(i, l, 4) = cos45*(v-u) |
---|
1018 | v_alpha(i, l, 5) = -u |
---|
1019 | v_alpha(i, l, 6) = -v_alpha(i, l, 2) |
---|
1020 | v_alpha(i, l, 7) = -v |
---|
1021 | v_alpha(i, l, 8) = -v_alpha(i, l, 4) |
---|
1022 | END DO |
---|
1023 | END DO |
---|
1024 | END IF |
---|
1025 | |
---|
1026 | RETURN |
---|
1027 | ! ----------------------------------------------------------------------- |
---|
1028 | END SUBROUTINE hines_wind |
---|
1029 | |
---|
1030 | SUBROUTINE hines_flux(flux_u, flux_v, drag_u, drag_v, alt, density, densb, & |
---|
1031 | m_alpha, ak_alpha, k_alpha, slope, naz, il1, il2, lev1, lev2, nlons, & |
---|
1032 | nlevs, nazmth, lorms) |
---|
1033 | IMPLICIT NONE |
---|
1034 | ! Calculate zonal and meridional components of the vertical flux |
---|
1035 | ! of horizontal momentum and corresponding wave drag (force per unit mass) |
---|
1036 | ! on a longitude by altitude grid for the Hines' Doppler spread |
---|
1037 | ! GWD parameterization scheme. |
---|
1038 | ! NOTE: only 4 or 8 azimuths can be used. |
---|
1039 | |
---|
1040 | ! Aug. 6/95 - C. McLandress |
---|
1041 | |
---|
1042 | ! Output arguements: |
---|
1043 | |
---|
1044 | ! * FLUX_U = zonal component of vertical momentum flux (Pascals) |
---|
1045 | ! * FLUX_V = meridional component of vertical momentum flux (Pascals) |
---|
1046 | ! * DRAG_U = zonal component of drag (m/s^2). |
---|
1047 | ! * DRAG_V = meridional component of drag (m/s^2). |
---|
1048 | |
---|
1049 | ! Input arguements: |
---|
1050 | |
---|
1051 | ! * ALT = altitudes (m). |
---|
1052 | ! * DENSITY = background density (kg/m^3). |
---|
1053 | ! * DENSB = background density at bottom level (kg/m^3). |
---|
1054 | ! * M_ALPHA = cutoff vertical wavenumber (1/m). |
---|
1055 | ! * AK_ALPHA = spectral amplitude factor (i.e., {AjKj} in m^4/s^2). |
---|
1056 | ! * K_ALPHA = horizontal wavenumber (1/m). |
---|
1057 | ! * SLOPE = slope of incident vertical wavenumber spectrum. |
---|
1058 | ! * NAZ = actual number of horizontal azimuths used (must be 4 or 8). |
---|
1059 | ! * IL1 = first longitudinal index to use (IL1 >= 1). |
---|
1060 | ! * IL2 = last longitudinal index to use (IL1 <= IL2 <= NLONS). |
---|
1061 | ! * LEV1 = first altitude level to use (LEV1 >=1). |
---|
1062 | ! * LEV2 = last altitude level to use (LEV1 < LEV2 <= NLEVS). |
---|
1063 | ! * NLONS = number of longitudes. |
---|
1064 | ! * NLEVS = number of vertical levels. |
---|
1065 | ! * NAZMTH = azimuthal array dimension (NAZMTH >= NAZ). |
---|
1066 | |
---|
1067 | ! * LORMS = .TRUE. for drag computation |
---|
1068 | |
---|
1069 | ! Constant in DATA statement. |
---|
1070 | |
---|
1071 | ! * COS45 = cosine of 45 degrees. |
---|
1072 | |
---|
1073 | ! Subroutine arguements. |
---|
1074 | |
---|
1075 | INTEGER naz, il1, il2, lev1, lev2 |
---|
1076 | INTEGER nlons, nlevs, nazmth |
---|
1077 | REAL slope |
---|
1078 | REAL flux_u(nlons, nlevs), flux_v(nlons, nlevs) |
---|
1079 | REAL drag_u(nlons, nlevs), drag_v(nlons, nlevs) |
---|
1080 | REAL alt(nlons, nlevs), density(nlons, nlevs), densb(nlons) |
---|
1081 | REAL m_alpha(nlons, nlevs, nazmth) |
---|
1082 | REAL ak_alpha(nlons, nazmth), k_alpha(nlons, nazmth) |
---|
1083 | |
---|
1084 | LOGICAL lorms(nlons) |
---|
1085 | |
---|
1086 | ! Internal variables. |
---|
1087 | |
---|
1088 | INTEGER i, l, lev1p, lev2m, lev2p |
---|
1089 | REAL cos45, prod2, prod4, prod6, prod8, dendz, dendz2 |
---|
1090 | DATA cos45/0.7071068/ |
---|
1091 | ! ----------------------------------------------------------------------- |
---|
1092 | |
---|
1093 | lev1p = lev1 + 1 |
---|
1094 | lev2m = lev2 - 1 |
---|
1095 | lev2p = lev2 + 1 |
---|
1096 | |
---|
1097 | ! Sum over azimuths for case where SLOPE = 1. |
---|
1098 | |
---|
1099 | IF (slope==1.) THEN |
---|
1100 | |
---|
1101 | ! Case with 4 azimuths. |
---|
1102 | |
---|
1103 | IF (naz==4) THEN |
---|
1104 | DO l = lev1, lev2 |
---|
1105 | DO i = il1, il2 |
---|
1106 | flux_u(i, l) = ak_alpha(i, 1)*k_alpha(i, 1)*m_alpha(i, l, 1) - & |
---|
1107 | ak_alpha(i, 3)*k_alpha(i, 3)*m_alpha(i, l, 3) |
---|
1108 | flux_v(i, l) = ak_alpha(i, 2)*k_alpha(i, 2)*m_alpha(i, l, 2) - & |
---|
1109 | ak_alpha(i, 4)*k_alpha(i, 4)*m_alpha(i, l, 4) |
---|
1110 | END DO |
---|
1111 | END DO |
---|
1112 | END IF |
---|
1113 | |
---|
1114 | ! Case with 8 azimuths. |
---|
1115 | |
---|
1116 | IF (naz==8) THEN |
---|
1117 | DO l = lev1, lev2 |
---|
1118 | DO i = il1, il2 |
---|
1119 | prod2 = ak_alpha(i, 2)*k_alpha(i, 2)*m_alpha(i, l, 2) |
---|
1120 | prod4 = ak_alpha(i, 4)*k_alpha(i, 4)*m_alpha(i, l, 4) |
---|
1121 | prod6 = ak_alpha(i, 6)*k_alpha(i, 6)*m_alpha(i, l, 6) |
---|
1122 | prod8 = ak_alpha(i, 8)*k_alpha(i, 8)*m_alpha(i, l, 8) |
---|
1123 | flux_u(i, l) = ak_alpha(i, 1)*k_alpha(i, 1)*m_alpha(i, l, 1) - & |
---|
1124 | ak_alpha(i, 5)*k_alpha(i, 5)*m_alpha(i, l, 5) + & |
---|
1125 | cos45*(prod2-prod4-prod6+prod8) |
---|
1126 | flux_v(i, l) = ak_alpha(i, 3)*k_alpha(i, 3)*m_alpha(i, l, 3) - & |
---|
1127 | ak_alpha(i, 7)*k_alpha(i, 7)*m_alpha(i, l, 7) + & |
---|
1128 | cos45*(prod2+prod4-prod6-prod8) |
---|
1129 | END DO |
---|
1130 | END DO |
---|
1131 | END IF |
---|
1132 | |
---|
1133 | END IF |
---|
1134 | |
---|
1135 | ! Sum over azimuths for case where SLOPE not equal to 1. |
---|
1136 | |
---|
1137 | IF (slope/=1.) THEN |
---|
1138 | |
---|
1139 | ! Case with 4 azimuths. |
---|
1140 | |
---|
1141 | IF (naz==4) THEN |
---|
1142 | DO l = lev1, lev2 |
---|
1143 | DO i = il1, il2 |
---|
1144 | flux_u(i, l) = ak_alpha(i, 1)*k_alpha(i, 1)* & |
---|
1145 | m_alpha(i, l, 1)**slope - ak_alpha(i, 3)*k_alpha(i, 3)*m_alpha(i, & |
---|
1146 | l, 3)**slope |
---|
1147 | flux_v(i, l) = ak_alpha(i, 2)*k_alpha(i, 2)* & |
---|
1148 | m_alpha(i, l, 2)**slope - ak_alpha(i, 4)*k_alpha(i, 4)*m_alpha(i, & |
---|
1149 | l, 4)**slope |
---|
1150 | END DO |
---|
1151 | END DO |
---|
1152 | END IF |
---|
1153 | |
---|
1154 | ! Case with 8 azimuths. |
---|
1155 | |
---|
1156 | IF (naz==8) THEN |
---|
1157 | DO l = lev1, lev2 |
---|
1158 | DO i = il1, il2 |
---|
1159 | prod2 = ak_alpha(i, 2)*k_alpha(i, 2)*m_alpha(i, l, 2)**slope |
---|
1160 | prod4 = ak_alpha(i, 4)*k_alpha(i, 4)*m_alpha(i, l, 4)**slope |
---|
1161 | prod6 = ak_alpha(i, 6)*k_alpha(i, 6)*m_alpha(i, l, 6)**slope |
---|
1162 | prod8 = ak_alpha(i, 8)*k_alpha(i, 8)*m_alpha(i, l, 8)**slope |
---|
1163 | flux_u(i, l) = ak_alpha(i, 1)*k_alpha(i, 1)* & |
---|
1164 | m_alpha(i, l, 1)**slope - ak_alpha(i, 5)*k_alpha(i, 5)*m_alpha(i, & |
---|
1165 | l, 5)**slope + cos45*(prod2-prod4-prod6+prod8) |
---|
1166 | flux_v(i, l) = ak_alpha(i, 3)*k_alpha(i, 3)* & |
---|
1167 | m_alpha(i, l, 3)**slope - ak_alpha(i, 7)*k_alpha(i, 7)*m_alpha(i, & |
---|
1168 | l, 7)**slope + cos45*(prod2+prod4-prod6-prod8) |
---|
1169 | END DO |
---|
1170 | END DO |
---|
1171 | END IF |
---|
1172 | |
---|
1173 | END IF |
---|
1174 | |
---|
1175 | ! Calculate flux from sum. |
---|
1176 | |
---|
1177 | DO l = lev1, lev2 |
---|
1178 | DO i = il1, il2 |
---|
1179 | flux_u(i, l) = flux_u(i, l)*densb(i)/slope |
---|
1180 | flux_v(i, l) = flux_v(i, l)*densb(i)/slope |
---|
1181 | END DO |
---|
1182 | END DO |
---|
1183 | |
---|
1184 | ! Calculate drag at intermediate levels using centered differences |
---|
1185 | |
---|
1186 | DO l = lev1p, lev2m |
---|
1187 | DO i = il1, il2 |
---|
1188 | IF (lorms(i)) THEN |
---|
1189 | ! cc DENDZ2 = DENSITY(I,L) * ( ALT(I,L+1) - ALT(I,L-1) ) |
---|
1190 | dendz2 = density(i, l)*(alt(i,l-1)-alt(i,l)) |
---|
1191 | ! cc DRAG_U(I,L) = - ( FLUX_U(I,L+1) - FLUX_U(I,L-1) ) / DENDZ2 |
---|
1192 | drag_u(i, l) = -(flux_u(i,l-1)-flux_u(i,l))/dendz2 |
---|
1193 | ! cc DRAG_V(I,L) = - ( FLUX_V(I,L+1) - FLUX_V(I,L-1) ) / DENDZ2 |
---|
1194 | drag_v(i, l) = -(flux_v(i,l-1)-flux_v(i,l))/dendz2 |
---|
1195 | |
---|
1196 | END IF |
---|
1197 | END DO |
---|
1198 | END DO |
---|
1199 | |
---|
1200 | ! Drag at first and last levels using one-side differences. |
---|
1201 | |
---|
1202 | DO i = il1, il2 |
---|
1203 | IF (lorms(i)) THEN |
---|
1204 | dendz = density(i, lev1)*(alt(i,lev1)-alt(i,lev1p)) |
---|
1205 | drag_u(i, lev1) = flux_u(i, lev1)/dendz |
---|
1206 | drag_v(i, lev1) = flux_v(i, lev1)/dendz |
---|
1207 | END IF |
---|
1208 | END DO |
---|
1209 | DO i = il1, il2 |
---|
1210 | IF (lorms(i)) THEN |
---|
1211 | dendz = density(i, lev2)*(alt(i,lev2m)-alt(i,lev2)) |
---|
1212 | drag_u(i, lev2) = -(flux_u(i,lev2m)-flux_u(i,lev2))/dendz |
---|
1213 | drag_v(i, lev2) = -(flux_v(i,lev2m)-flux_v(i,lev2))/dendz |
---|
1214 | END IF |
---|
1215 | END DO |
---|
1216 | IF (nlevs>lev2) THEN |
---|
1217 | DO i = il1, il2 |
---|
1218 | IF (lorms(i)) THEN |
---|
1219 | dendz = density(i, lev2p)*(alt(i,lev2)-alt(i,lev2p)) |
---|
1220 | drag_u(i, lev2p) = -flux_u(i, lev2)/dendz |
---|
1221 | drag_v(i, lev2p) = -flux_v(i, lev2)/dendz |
---|
1222 | END IF |
---|
1223 | END DO |
---|
1224 | END IF |
---|
1225 | |
---|
1226 | RETURN |
---|
1227 | ! ----------------------------------------------------------------------- |
---|
1228 | END SUBROUTINE hines_flux |
---|
1229 | |
---|
1230 | SUBROUTINE hines_heat(heat, diffco, m_alpha, dmdz_alpha, ak_alpha, k_alpha, & |
---|
1231 | bvfreq, density, densb, sigma_t, visc_mol, kstar, slope, f2, f3, f5, f6, & |
---|
1232 | naz, il1, il2, lev1, lev2, nlons, nlevs, nazmth) |
---|
1233 | IMPLICIT NONE |
---|
1234 | ! This routine calculates the gravity wave induced heating and |
---|
1235 | ! diffusion coefficient on a longitude by altitude grid for |
---|
1236 | ! the Hines' Doppler spread gravity wave drag parameterization scheme. |
---|
1237 | |
---|
1238 | ! Aug. 6/95 - C. McLandress |
---|
1239 | |
---|
1240 | ! Output arguements: |
---|
1241 | |
---|
1242 | ! * HEAT = gravity wave heating (K/sec). |
---|
1243 | ! * DIFFCO = diffusion coefficient (m^2/sec) |
---|
1244 | |
---|
1245 | ! Input arguements: |
---|
1246 | |
---|
1247 | ! * M_ALPHA = cutoff vertical wavenumber (1/m). |
---|
1248 | ! * DMDZ_ALPHA = vertical derivative of cutoff wavenumber. |
---|
1249 | ! * AK_ALPHA = spectral amplitude factor of each azimuth |
---|
1250 | ! (i.e., {AjKj} in m^4/s^2). |
---|
1251 | ! * K_ALPHA = horizontal wavenumber of each azimuth (1/m). |
---|
1252 | ! * BVFREQ = background Brunt Vassala frequency (rad/sec). |
---|
1253 | ! * DENSITY = background density (kg/m^3). |
---|
1254 | ! * DENSB = background density at bottom level (kg/m^3). |
---|
1255 | ! * SIGMA_T = total rms horizontal wind (m/s). |
---|
1256 | ! * VISC_MOL = molecular viscosity (m^2/s). |
---|
1257 | ! * KSTAR = typical gravity wave horizontal wavenumber (1/m). |
---|
1258 | ! * SLOPE = slope of incident vertical wavenumber spectrum. |
---|
1259 | ! * F2,F3,F5,F6 = Hines's fudge factors. |
---|
1260 | ! * NAZ = actual number of horizontal azimuths used. |
---|
1261 | ! * IL1 = first longitudinal index to use (IL1 >= 1). |
---|
1262 | ! * IL2 = last longitudinal index to use (IL1 <= IL2 <= NLONS). |
---|
1263 | ! * LEV1 = first altitude level to use (LEV1 >=1). |
---|
1264 | ! * LEV2 = last altitude level to use (LEV1 < LEV2 <= NLEVS). |
---|
1265 | ! * NLONS = number of longitudes. |
---|
1266 | ! * NLEVS = number of vertical levels. |
---|
1267 | ! * NAZMTH = azimuthal array dimension (NAZMTH >= NAZ). |
---|
1268 | |
---|
1269 | INTEGER naz, il1, il2, lev1, lev2, nlons, nlevs, nazmth |
---|
1270 | REAL kstar(nlons), slope, f2, f3, f5, f6 |
---|
1271 | REAL heat(nlons, nlevs), diffco(nlons, nlevs) |
---|
1272 | REAL m_alpha(nlons, nlevs, nazmth), dmdz_alpha(nlons, nlevs, nazmth) |
---|
1273 | REAL ak_alpha(nlons, nazmth), k_alpha(nlons, nazmth) |
---|
1274 | REAL bvfreq(nlons, nlevs), density(nlons, nlevs), densb(nlons) |
---|
1275 | REAL sigma_t(nlons, nlevs), visc_mol(nlons, nlevs) |
---|
1276 | |
---|
1277 | ! Internal variables. |
---|
1278 | |
---|
1279 | INTEGER i, l, n |
---|
1280 | REAL m_sub_m_turb, m_sub_m_mol, m_sub_m, heatng |
---|
1281 | REAL visc, visc_min, cpgas, sm1 |
---|
1282 | |
---|
1283 | ! specific heat at constant pressure |
---|
1284 | |
---|
1285 | DATA cpgas/1004./ |
---|
1286 | |
---|
1287 | ! minimum permissible viscosity |
---|
1288 | |
---|
1289 | DATA visc_min/1.E-10/ |
---|
1290 | ! ----------------------------------------------------------------------- |
---|
1291 | |
---|
1292 | ! Initialize heating array. |
---|
1293 | |
---|
1294 | DO l = 1, nlevs |
---|
1295 | DO i = 1, nlons |
---|
1296 | heat(i, l) = 0. |
---|
1297 | END DO |
---|
1298 | END DO |
---|
1299 | |
---|
1300 | ! Perform sum over azimuths for case where SLOPE = 1. |
---|
1301 | |
---|
1302 | IF (slope==1.) THEN |
---|
1303 | DO n = 1, naz |
---|
1304 | DO l = lev1, lev2 |
---|
1305 | DO i = il1, il2 |
---|
1306 | heat(i, l) = heat(i, l) + ak_alpha(i, n)*k_alpha(i, n)*dmdz_alpha(i & |
---|
1307 | , l, n) |
---|
1308 | END DO |
---|
1309 | END DO |
---|
1310 | END DO |
---|
1311 | END IF |
---|
1312 | |
---|
1313 | ! Perform sum over azimuths for case where SLOPE not 1. |
---|
1314 | |
---|
1315 | IF (slope/=1.) THEN |
---|
1316 | sm1 = slope - 1. |
---|
1317 | DO n = 1, naz |
---|
1318 | DO l = lev1, lev2 |
---|
1319 | DO i = il1, il2 |
---|
1320 | heat(i, l) = heat(i, l) + ak_alpha(i, n)*k_alpha(i, n)*m_alpha(i, l & |
---|
1321 | , n)**sm1*dmdz_alpha(i, l, n) |
---|
1322 | END DO |
---|
1323 | END DO |
---|
1324 | END DO |
---|
1325 | END IF |
---|
1326 | |
---|
1327 | ! Heating and diffusion. |
---|
1328 | |
---|
1329 | DO l = lev1, lev2 |
---|
1330 | DO i = il1, il2 |
---|
1331 | |
---|
1332 | ! Maximum permissible value of cutoff wavenumber is the smaller |
---|
1333 | ! of the instability-induced wavenumber (M_SUB_M_TURB) and |
---|
1334 | ! that imposed by molecular viscosity (M_SUB_M_MOL). |
---|
1335 | |
---|
1336 | visc = amax1(visc_mol(i,l), visc_min) |
---|
1337 | m_sub_m_turb = bvfreq(i, l)/(f2*sigma_t(i,l)) |
---|
1338 | m_sub_m_mol = (bvfreq(i,l)*kstar(i)/visc)**0.33333333/f3 |
---|
1339 | m_sub_m = amin1(m_sub_m_turb, m_sub_m_mol) |
---|
1340 | |
---|
1341 | heatng = -heat(i, l)*f5*bvfreq(i, l)/m_sub_m*densb(i)/density(i, l) |
---|
1342 | diffco(i, l) = f6*heatng**0.33333333/m_sub_m**1.33333333 |
---|
1343 | heat(i, l) = heatng/cpgas |
---|
1344 | |
---|
1345 | END DO |
---|
1346 | END DO |
---|
1347 | |
---|
1348 | RETURN |
---|
1349 | ! ----------------------------------------------------------------------- |
---|
1350 | END SUBROUTINE hines_heat |
---|
1351 | |
---|
1352 | SUBROUTINE hines_sigma(sigma_t, sigma_alpha, sigsqh_alpha, naz, lev, il1, & |
---|
1353 | il2, nlons, nlevs, nazmth) |
---|
1354 | IMPLICIT NONE |
---|
1355 | ! This routine calculates the total rms and azimuthal rms horizontal |
---|
1356 | ! velocities at a given level on a longitude by altitude grid for |
---|
1357 | ! the Hines' Doppler spread GWD parameterization scheme. |
---|
1358 | ! NOTE: only four or eight azimuths can be used. |
---|
1359 | |
---|
1360 | ! Aug. 7/95 - C. McLandress |
---|
1361 | |
---|
1362 | ! Output arguements: |
---|
1363 | |
---|
1364 | ! * SIGMA_T = total rms horizontal wind (m/s). |
---|
1365 | ! * SIGMA_ALPHA = total rms wind in each azimuth (m/s). |
---|
1366 | |
---|
1367 | ! Input arguements: |
---|
1368 | |
---|
1369 | ! * SIGSQH_ALPHA = portion of wind variance from waves having wave |
---|
1370 | ! * normals in the alpha azimuth (m/s). |
---|
1371 | ! * NAZ = actual number of horizontal azimuths used (must be 4 or 8). |
---|
1372 | ! * LEV = altitude level to process. |
---|
1373 | ! * IL1 = first longitudinal index to use (IL1 >= 1). |
---|
1374 | ! * IL2 = last longitudinal index to use (IL1 <= IL2 <= NLONS). |
---|
1375 | ! * NLONS = number of longitudes. |
---|
1376 | ! * NLEVS = number of vertical levels. |
---|
1377 | ! * NAZMTH = azimuthal array dimension (NAZMTH >= NAZ). |
---|
1378 | |
---|
1379 | ! Subroutine arguements. |
---|
1380 | |
---|
1381 | INTEGER lev, naz, il1, il2 |
---|
1382 | INTEGER nlons, nlevs, nazmth |
---|
1383 | REAL sigma_t(nlons, nlevs) |
---|
1384 | REAL sigma_alpha(nlons, nlevs, nazmth) |
---|
1385 | REAL sigsqh_alpha(nlons, nlevs, nazmth) |
---|
1386 | |
---|
1387 | ! Internal variables. |
---|
1388 | |
---|
1389 | INTEGER i, n |
---|
1390 | REAL sum_even, sum_odd |
---|
1391 | ! ----------------------------------------------------------------------- |
---|
1392 | |
---|
1393 | ! Calculate azimuthal rms velocity for the 4 azimuth case. |
---|
1394 | |
---|
1395 | IF (naz==4) THEN |
---|
1396 | DO i = il1, il2 |
---|
1397 | sigma_alpha(i, lev, 1) = sqrt(sigsqh_alpha(i,lev,1)+sigsqh_alpha(i,lev, & |
---|
1398 | 3)) |
---|
1399 | sigma_alpha(i, lev, 2) = sqrt(sigsqh_alpha(i,lev,2)+sigsqh_alpha(i,lev, & |
---|
1400 | 4)) |
---|
1401 | sigma_alpha(i, lev, 3) = sigma_alpha(i, lev, 1) |
---|
1402 | sigma_alpha(i, lev, 4) = sigma_alpha(i, lev, 2) |
---|
1403 | END DO |
---|
1404 | END IF |
---|
1405 | |
---|
1406 | ! Calculate azimuthal rms velocity for the 8 azimuth case. |
---|
1407 | |
---|
1408 | IF (naz==8) THEN |
---|
1409 | DO i = il1, il2 |
---|
1410 | sum_odd = (sigsqh_alpha(i,lev,1)+sigsqh_alpha(i,lev,3)+ & |
---|
1411 | sigsqh_alpha(i,lev,5)+sigsqh_alpha(i,lev,7))/2. |
---|
1412 | sum_even = (sigsqh_alpha(i,lev,2)+sigsqh_alpha(i,lev,4)+ & |
---|
1413 | sigsqh_alpha(i,lev,6)+sigsqh_alpha(i,lev,8))/2. |
---|
1414 | sigma_alpha(i, lev, 1) = sqrt(sigsqh_alpha(i,lev,1)+sigsqh_alpha(i,lev, & |
---|
1415 | 5)+sum_even) |
---|
1416 | sigma_alpha(i, lev, 2) = sqrt(sigsqh_alpha(i,lev,2)+sigsqh_alpha(i,lev, & |
---|
1417 | 6)+sum_odd) |
---|
1418 | sigma_alpha(i, lev, 3) = sqrt(sigsqh_alpha(i,lev,3)+sigsqh_alpha(i,lev, & |
---|
1419 | 7)+sum_even) |
---|
1420 | sigma_alpha(i, lev, 4) = sqrt(sigsqh_alpha(i,lev,4)+sigsqh_alpha(i,lev, & |
---|
1421 | 8)+sum_odd) |
---|
1422 | sigma_alpha(i, lev, 5) = sigma_alpha(i, lev, 1) |
---|
1423 | sigma_alpha(i, lev, 6) = sigma_alpha(i, lev, 2) |
---|
1424 | sigma_alpha(i, lev, 7) = sigma_alpha(i, lev, 3) |
---|
1425 | sigma_alpha(i, lev, 8) = sigma_alpha(i, lev, 4) |
---|
1426 | END DO |
---|
1427 | END IF |
---|
1428 | |
---|
1429 | ! Calculate total rms velocity. |
---|
1430 | |
---|
1431 | DO i = il1, il2 |
---|
1432 | sigma_t(i, lev) = 0. |
---|
1433 | END DO |
---|
1434 | DO n = 1, naz |
---|
1435 | DO i = il1, il2 |
---|
1436 | sigma_t(i, lev) = sigma_t(i, lev) + sigsqh_alpha(i, lev, n) |
---|
1437 | END DO |
---|
1438 | END DO |
---|
1439 | DO i = il1, il2 |
---|
1440 | sigma_t(i, lev) = sqrt(sigma_t(i,lev)) |
---|
1441 | END DO |
---|
1442 | |
---|
1443 | RETURN |
---|
1444 | ! ----------------------------------------------------------------------- |
---|
1445 | END SUBROUTINE hines_sigma |
---|
1446 | |
---|
1447 | SUBROUTINE hines_intgrl(i_alpha, v_alpha, m_alpha, bvfb, slope, naz, lev, & |
---|
1448 | il1, il2, nlons, nlevs, nazmth, lorms) |
---|
1449 | IMPLICIT NONE |
---|
1450 | ! This routine calculates the vertical wavenumber integral |
---|
1451 | ! for a single vertical level at each azimuth on a longitude grid |
---|
1452 | ! for the Hines' Doppler spread GWD parameterization scheme. |
---|
1453 | ! NOTE: (1) only spectral slopes of 1, 1.5 or 2 are permitted. |
---|
1454 | ! (2) the integral is written in terms of the product QM |
---|
1455 | ! which by construction is always less than 1. Series |
---|
1456 | ! solutions are used for small |QM| and analytical solutions |
---|
1457 | ! for remaining values. |
---|
1458 | |
---|
1459 | ! Aug. 8/95 - C. McLandress |
---|
1460 | |
---|
1461 | ! Output arguement: |
---|
1462 | |
---|
1463 | ! * I_ALPHA = Hines' integral. |
---|
1464 | |
---|
1465 | ! Input arguements: |
---|
1466 | |
---|
1467 | ! * V_ALPHA = azimuthal wind component (m/s). |
---|
1468 | ! * M_ALPHA = azimuthal cutoff vertical wavenumber (1/m). |
---|
1469 | ! * BVFB = background Brunt Vassala frequency at model bottom. |
---|
1470 | ! * SLOPE = slope of initial vertical wavenumber spectrum |
---|
1471 | ! * (must use SLOPE = 1., 1.5 or 2.) |
---|
1472 | ! * NAZ = actual number of horizontal azimuths used. |
---|
1473 | ! * LEV = altitude level to process. |
---|
1474 | ! * IL1 = first longitudinal index to use (IL1 >= 1). |
---|
1475 | ! * IL2 = last longitudinal index to use (IL1 <= IL2 <= NLONS). |
---|
1476 | ! * NLONS = number of longitudes. |
---|
1477 | ! * NLEVS = number of vertical levels. |
---|
1478 | ! * NAZMTH = azimuthal array dimension (NAZMTH >= NAZ). |
---|
1479 | |
---|
1480 | ! * LORMS = .TRUE. for drag computation |
---|
1481 | |
---|
1482 | ! Constants in DATA statements: |
---|
1483 | |
---|
1484 | ! * QMIN = minimum value of Q_ALPHA (avoids indeterminant form of integral) |
---|
1485 | ! * QM_MIN = minimum value of Q_ALPHA * M_ALPHA (used to avoid numerical |
---|
1486 | ! * problems). |
---|
1487 | |
---|
1488 | INTEGER lev, naz, il1, il2, nlons, nlevs, nazmth |
---|
1489 | REAL i_alpha(nlons, nazmth) |
---|
1490 | REAL v_alpha(nlons, nlevs, nazmth) |
---|
1491 | REAL m_alpha(nlons, nlevs, nazmth) |
---|
1492 | REAL bvfb(nlons), slope |
---|
1493 | |
---|
1494 | LOGICAL lorms(nlons) |
---|
1495 | |
---|
1496 | ! Internal variables. |
---|
1497 | |
---|
1498 | INTEGER i, n |
---|
1499 | REAL q_alpha, qm, sqrtqm, q_min, qm_min |
---|
1500 | |
---|
1501 | DATA q_min/1.0/, qm_min/0.01/ |
---|
1502 | ! ----------------------------------------------------------------------- |
---|
1503 | |
---|
1504 | ! For integer value SLOPE = 1. |
---|
1505 | |
---|
1506 | IF (slope==1.) THEN |
---|
1507 | |
---|
1508 | DO n = 1, naz |
---|
1509 | DO i = il1, il2 |
---|
1510 | IF (lorms(i)) THEN |
---|
1511 | |
---|
1512 | q_alpha = v_alpha(i, lev, n)/bvfb(i) |
---|
1513 | qm = q_alpha*m_alpha(i, lev, n) |
---|
1514 | |
---|
1515 | ! If |QM| is small then use first 4 terms series of Taylor series |
---|
1516 | ! expansion of integral in order to avoid indeterminate form of |
---|
1517 | ! integral, |
---|
1518 | ! otherwise use analytical form of integral. |
---|
1519 | |
---|
1520 | IF (abs(q_alpha)<q_min .OR. abs(qm)<qm_min) THEN |
---|
1521 | IF (q_alpha==0.) THEN |
---|
1522 | i_alpha(i, n) = m_alpha(i, lev, n)**2/2. |
---|
1523 | ELSE |
---|
1524 | i_alpha(i, n) = (qm**2/2.+qm**3/3.+qm**4/4.+qm**5/5.)/ & |
---|
1525 | q_alpha**2 |
---|
1526 | END IF |
---|
1527 | ELSE |
---|
1528 | i_alpha(i, n) = -(alog(1.-qm)+qm)/q_alpha**2 |
---|
1529 | END IF |
---|
1530 | |
---|
1531 | END IF |
---|
1532 | END DO |
---|
1533 | END DO |
---|
1534 | |
---|
1535 | END IF |
---|
1536 | |
---|
1537 | ! For integer value SLOPE = 2. |
---|
1538 | |
---|
1539 | IF (slope==2.) THEN |
---|
1540 | |
---|
1541 | DO n = 1, naz |
---|
1542 | DO i = il1, il2 |
---|
1543 | IF (lorms(i)) THEN |
---|
1544 | |
---|
1545 | q_alpha = v_alpha(i, lev, n)/bvfb(i) |
---|
1546 | qm = q_alpha*m_alpha(i, lev, n) |
---|
1547 | |
---|
1548 | ! If |QM| is small then use first 4 terms series of Taylor series |
---|
1549 | ! expansion of integral in order to avoid indeterminate form of |
---|
1550 | ! integral, |
---|
1551 | ! otherwise use analytical form of integral. |
---|
1552 | |
---|
1553 | IF (abs(q_alpha)<q_min .OR. abs(qm)<qm_min) THEN |
---|
1554 | IF (q_alpha==0.) THEN |
---|
1555 | i_alpha(i, n) = m_alpha(i, lev, n)**3/3. |
---|
1556 | ELSE |
---|
1557 | i_alpha(i, n) = (qm**3/3.+qm**4/4.+qm**5/5.+qm**6/6.)/ & |
---|
1558 | q_alpha**3 |
---|
1559 | END IF |
---|
1560 | ELSE |
---|
1561 | i_alpha(i, n) = -(alog(1.-qm)+qm+qm**2/2.)/q_alpha**3 |
---|
1562 | END IF |
---|
1563 | |
---|
1564 | END IF |
---|
1565 | END DO |
---|
1566 | END DO |
---|
1567 | |
---|
1568 | END IF |
---|
1569 | |
---|
1570 | ! For real value SLOPE = 1.5 |
---|
1571 | |
---|
1572 | IF (slope==1.5) THEN |
---|
1573 | |
---|
1574 | DO n = 1, naz |
---|
1575 | DO i = il1, il2 |
---|
1576 | IF (lorms(i)) THEN |
---|
1577 | |
---|
1578 | q_alpha = v_alpha(i, lev, n)/bvfb(i) |
---|
1579 | qm = q_alpha*m_alpha(i, lev, n) |
---|
1580 | |
---|
1581 | ! If |QM| is small then use first 4 terms series of Taylor series |
---|
1582 | ! expansion of integral in order to avoid indeterminate form of |
---|
1583 | ! integral, |
---|
1584 | ! otherwise use analytical form of integral. |
---|
1585 | |
---|
1586 | IF (abs(q_alpha)<q_min .OR. abs(qm)<qm_min) THEN |
---|
1587 | IF (q_alpha==0.) THEN |
---|
1588 | i_alpha(i, n) = m_alpha(i, lev, n)**2.5/2.5 |
---|
1589 | ELSE |
---|
1590 | i_alpha(i, n) = (qm/2.5+qm**2/3.5+qm**3/4.5+qm**4/5.5)* & |
---|
1591 | m_alpha(i, lev, n)**1.5/q_alpha |
---|
1592 | END IF |
---|
1593 | ELSE |
---|
1594 | qm = abs(qm) |
---|
1595 | sqrtqm = sqrt(qm) |
---|
1596 | IF (q_alpha>=0.) THEN |
---|
1597 | i_alpha(i, n) = (alog((1.+sqrtqm)/(1.-sqrtqm))-2.*sqrtqm*(1.+qm & |
---|
1598 | /3.))/q_alpha**2.5 |
---|
1599 | ELSE |
---|
1600 | i_alpha(i, n) = 2.*(atan(sqrtqm)+sqrtqm*(qm/3.-1.))/ & |
---|
1601 | abs(q_alpha)**2.5 |
---|
1602 | END IF |
---|
1603 | END IF |
---|
1604 | |
---|
1605 | END IF |
---|
1606 | END DO |
---|
1607 | END DO |
---|
1608 | |
---|
1609 | END IF |
---|
1610 | |
---|
1611 | ! If integral is negative (which in principal should not happen) then |
---|
1612 | ! print a message and some info since execution will abort when calculating |
---|
1613 | ! the variances. |
---|
1614 | |
---|
1615 | ! DO 80 N = 1,NAZ |
---|
1616 | ! DO 70 I = IL1,IL2 |
---|
1617 | ! IF (I_ALPHA(I,N).LT.0.) THEN |
---|
1618 | ! WRITE (6,*) |
---|
1619 | ! WRITE (6,*) '******************************' |
---|
1620 | ! WRITE (6,*) 'Hines integral I_ALPHA < 0 ' |
---|
1621 | ! WRITE (6,*) ' longitude I=',I |
---|
1622 | ! WRITE (6,*) ' azimuth N=',N |
---|
1623 | ! WRITE (6,*) ' level LEV=',LEV |
---|
1624 | ! WRITE (6,*) ' I_ALPHA =',I_ALPHA(I,N) |
---|
1625 | ! WRITE (6,*) ' V_ALPHA =',V_ALPHA(I,LEV,N) |
---|
1626 | ! WRITE (6,*) ' M_ALPHA =',M_ALPHA(I,LEV,N) |
---|
1627 | ! WRITE (6,*) ' Q_ALPHA =',V_ALPHA(I,LEV,N) / BVFB(I) |
---|
1628 | ! WRITE (6,*) ' QM =',V_ALPHA(I,LEV,N) / BVFB(I) |
---|
1629 | ! ^ * M_ALPHA(I,LEV,N) |
---|
1630 | ! WRITE (6,*) '******************************' |
---|
1631 | ! END IF |
---|
1632 | ! 70 CONTINUE |
---|
1633 | ! 80 CONTINUE |
---|
1634 | |
---|
1635 | RETURN |
---|
1636 | ! ----------------------------------------------------------------------- |
---|
1637 | END SUBROUTINE hines_intgrl |
---|
1638 | |
---|
1639 | SUBROUTINE hines_setup(naz, slope, f1, f2, f3, f5, f6, kstar, icutoff, & |
---|
1640 | alt_cutoff, smco, nsmax, iheatcal, k_alpha, ierror, nmessg, nlons, & |
---|
1641 | nazmth, coslat) |
---|
1642 | IMPLICIT NONE |
---|
1643 | ! This routine specifies various parameters needed for the |
---|
1644 | ! the Hines' Doppler spread gravity wave drag parameterization scheme. |
---|
1645 | |
---|
1646 | ! Aug. 8/95 - C. McLandress |
---|
1647 | |
---|
1648 | ! Output arguements: |
---|
1649 | |
---|
1650 | ! * NAZ = actual number of horizontal azimuths used |
---|
1651 | ! * (code set up presently for only NAZ = 4 or 8). |
---|
1652 | ! * SLOPE = slope of incident vertical wavenumber spectrum |
---|
1653 | ! * (code set up presently for SLOPE 1., 1.5 or 2.). |
---|
1654 | ! * F1 = "fudge factor" used in calculation of trial value of |
---|
1655 | ! * azimuthal cutoff wavenumber M_ALPHA (1.2 <= F1 <= 1.9). |
---|
1656 | ! * F2 = "fudge factor" used in calculation of maximum |
---|
1657 | ! * permissible instabiliy-induced cutoff wavenumber |
---|
1658 | ! * M_SUB_M_TURB (0.1 <= F2 <= 1.4). |
---|
1659 | ! * F3 = "fudge factor" used in calculation of maximum |
---|
1660 | ! * permissible molecular viscosity-induced cutoff wavenumber |
---|
1661 | ! * M_SUB_M_MOL (0.1 <= F2 <= 1.4). |
---|
1662 | ! * F5 = "fudge factor" used in calculation of heating rate |
---|
1663 | ! * (1 <= F5 <= 3). |
---|
1664 | ! * F6 = "fudge factor" used in calculation of turbulent |
---|
1665 | ! * diffusivity coefficient. |
---|
1666 | ! * KSTAR = typical gravity wave horizontal wavenumber (1/m) |
---|
1667 | ! * used in calculation of M_SUB_M_TURB. |
---|
1668 | ! * ICUTOFF = 1 to exponentially damp off GWD, heating and diffusion |
---|
1669 | ! * arrays above ALT_CUTOFF; otherwise arrays not modified. |
---|
1670 | ! * ALT_CUTOFF = altitude in meters above which exponential decay applied. |
---|
1671 | ! * SMCO = smoother used to smooth cutoff vertical wavenumbers |
---|
1672 | ! * and total rms winds before calculating drag or heating. |
---|
1673 | ! * (==> a 1:SMCO:1 stencil used; SMCO >= 1.). |
---|
1674 | ! * NSMAX = number of times smoother applied ( >= 1), |
---|
1675 | ! * = 0 means no smoothing performed. |
---|
1676 | ! * IHEATCAL = 1 to calculate heating rates and diffusion coefficient. |
---|
1677 | ! * = 0 means only drag and flux calculated. |
---|
1678 | ! * K_ALPHA = horizontal wavenumber of each azimuth (1/m) which |
---|
1679 | ! * is set here to KSTAR. |
---|
1680 | ! * IERROR = error flag. |
---|
1681 | ! * = 0 no errors. |
---|
1682 | ! * = 10 ==> NAZ > NAZMTH |
---|
1683 | ! * = 20 ==> invalid number of azimuths (NAZ must be 4 or 8). |
---|
1684 | ! * = 30 ==> invalid slope (SLOPE must be 1., 1.5 or 2.). |
---|
1685 | ! * = 40 ==> invalid smoother (SMCO must be >= 1.) |
---|
1686 | |
---|
1687 | ! Input arguements: |
---|
1688 | |
---|
1689 | ! * NMESSG = output unit number where messages to be printed. |
---|
1690 | ! * NLONS = number of longitudes. |
---|
1691 | ! * NAZMTH = azimuthal array dimension (NAZMTH >= NAZ). |
---|
1692 | |
---|
1693 | INTEGER naz, nlons, nazmth, iheatcal, icutoff |
---|
1694 | INTEGER nmessg, nsmax, ierror |
---|
1695 | REAL kstar(nlons), slope, f1, f2, f3, f5, f6, alt_cutoff, smco |
---|
1696 | REAL k_alpha(nlons, nazmth), coslat(nlons) |
---|
1697 | REAL ksmin, ksmax |
---|
1698 | |
---|
1699 | ! Internal variables. |
---|
1700 | |
---|
1701 | INTEGER i, n |
---|
1702 | ! ----------------------------------------------------------------------- |
---|
1703 | |
---|
1704 | ! Specify constants. |
---|
1705 | |
---|
1706 | naz = 8 |
---|
1707 | slope = 1. |
---|
1708 | f1 = 1.5 |
---|
1709 | f2 = 0.3 |
---|
1710 | f3 = 1.0 |
---|
1711 | f5 = 3.0 |
---|
1712 | f6 = 1.0 |
---|
1713 | ksmin = 1.E-5 |
---|
1714 | ksmax = 1.E-4 |
---|
1715 | DO i = 1, nlons |
---|
1716 | kstar(i) = ksmin/(coslat(i)+(ksmin/ksmax)) |
---|
1717 | END DO |
---|
1718 | icutoff = 1 |
---|
1719 | alt_cutoff = 105.E3 |
---|
1720 | smco = 2.0 |
---|
1721 | ! SMCO = 1.0 |
---|
1722 | nsmax = 5 |
---|
1723 | ! NSMAX = 2 |
---|
1724 | iheatcal = 0 |
---|
1725 | |
---|
1726 | ! Print information to output file. |
---|
1727 | |
---|
1728 | ! WRITE (NMESSG,6000) |
---|
1729 | ! 6000 FORMAT (/' Subroutine HINES_SETUP:') |
---|
1730 | ! WRITE (NMESSG,*) ' SLOPE = ', SLOPE |
---|
1731 | ! WRITE (NMESSG,*) ' NAZ = ', NAZ |
---|
1732 | ! WRITE (NMESSG,*) ' F1,F2,F3 = ', F1, F2, F3 |
---|
1733 | ! WRITE (NMESSG,*) ' F5,F6 = ', F5, F6 |
---|
1734 | ! WRITE (NMESSG,*) ' KSTAR = ', KSTAR |
---|
1735 | ! > ,' COSLAT = ', COSLAT |
---|
1736 | ! IF (ICUTOFF .EQ. 1) THEN |
---|
1737 | ! WRITE (NMESSG,*) ' Drag exponentially damped above ', |
---|
1738 | ! & ALT_CUTOFF/1.E3 |
---|
1739 | ! END IF |
---|
1740 | ! IF (NSMAX.LT.1 ) THEN |
---|
1741 | ! WRITE (NMESSG,*) ' No smoothing of cutoff wavenumbers, etc' |
---|
1742 | ! ELSE |
---|
1743 | ! WRITE (NMESSG,*) ' Cutoff wavenumbers and sig_t smoothed:' |
---|
1744 | ! WRITE (NMESSG,*) ' SMCO =', SMCO |
---|
1745 | ! WRITE (NMESSG,*) ' NSMAX =', NSMAX |
---|
1746 | ! END IF |
---|
1747 | |
---|
1748 | ! Check that things are setup correctly and log error if not |
---|
1749 | |
---|
1750 | ierror = 0 |
---|
1751 | IF (naz>nazmth) ierror = 10 |
---|
1752 | IF (naz/=4 .AND. naz/=8) ierror = 20 |
---|
1753 | IF (slope/=1. .AND. slope/=1.5 .AND. slope/=2.) ierror = 30 |
---|
1754 | IF (smco<1.) ierror = 40 |
---|
1755 | |
---|
1756 | ! Use single value for azimuthal-dependent horizontal wavenumber. |
---|
1757 | |
---|
1758 | DO n = 1, naz |
---|
1759 | DO i = 1, nlons |
---|
1760 | k_alpha(i, n) = kstar(i) |
---|
1761 | END DO |
---|
1762 | END DO |
---|
1763 | |
---|
1764 | RETURN |
---|
1765 | ! ----------------------------------------------------------------------- |
---|
1766 | END SUBROUTINE hines_setup |
---|
1767 | |
---|
1768 | SUBROUTINE hines_print(flux_u, flux_v, drag_u, drag_v, alt, sigma_t, & |
---|
1769 | sigma_alpha, v_alpha, m_alpha, iu_print, iv_print, nmessg, ilprt1, & |
---|
1770 | ilprt2, levprt1, levprt2, naz, nlons, nlevs, nazmth) |
---|
1771 | IMPLICIT NONE |
---|
1772 | ! Print out altitude profiles of various quantities from |
---|
1773 | ! Hines' Doppler spread gravity wave drag parameterization scheme. |
---|
1774 | ! (NOTE: only for NAZ = 4 or 8). |
---|
1775 | |
---|
1776 | ! Aug. 8/95 - C. McLandress |
---|
1777 | |
---|
1778 | ! Input arguements: |
---|
1779 | |
---|
1780 | ! * IU_PRINT = 1 to print out values in east-west direction. |
---|
1781 | ! * IV_PRINT = 1 to print out values in north-south direction. |
---|
1782 | ! * NMESSG = unit number for printed output. |
---|
1783 | ! * ILPRT1 = first longitudinal index to print. |
---|
1784 | ! * ILPRT2 = last longitudinal index to print. |
---|
1785 | ! * LEVPRT1 = first altitude level to print. |
---|
1786 | ! * LEVPRT2 = last altitude level to print. |
---|
1787 | |
---|
1788 | INTEGER naz, ilprt1, ilprt2, levprt1, levprt2 |
---|
1789 | INTEGER nlons, nlevs, nazmth |
---|
1790 | INTEGER iu_print, iv_print, nmessg |
---|
1791 | REAL flux_u(nlons, nlevs), flux_v(nlons, nlevs) |
---|
1792 | REAL drag_u(nlons, nlevs), drag_v(nlons, nlevs) |
---|
1793 | REAL alt(nlons, nlevs), sigma_t(nlons, nlevs) |
---|
1794 | REAL sigma_alpha(nlons, nlevs, nazmth) |
---|
1795 | REAL v_alpha(nlons, nlevs, nazmth), m_alpha(nlons, nlevs, nazmth) |
---|
1796 | |
---|
1797 | ! Internal variables. |
---|
1798 | |
---|
1799 | INTEGER n_east, n_west, n_north, n_south |
---|
1800 | INTEGER i, l |
---|
1801 | ! ----------------------------------------------------------------------- |
---|
1802 | |
---|
1803 | ! Azimuthal indices of cardinal directions. |
---|
1804 | |
---|
1805 | n_east = 1 |
---|
1806 | IF (naz==4) THEN |
---|
1807 | n_west = 3 |
---|
1808 | n_north = 2 |
---|
1809 | n_south = 4 |
---|
1810 | ELSE IF (naz==8) THEN |
---|
1811 | n_west = 5 |
---|
1812 | n_north = 3 |
---|
1813 | n_south = 7 |
---|
1814 | END IF |
---|
1815 | |
---|
1816 | ! Print out values for range of longitudes. |
---|
1817 | |
---|
1818 | DO i = ilprt1, ilprt2 |
---|
1819 | |
---|
1820 | ! Print east-west wind, sigmas, cutoff wavenumbers, flux and drag. |
---|
1821 | |
---|
1822 | IF (iu_print==1) THEN |
---|
1823 | WRITE (nmessg, *) |
---|
1824 | WRITE (nmessg, 6001) i |
---|
1825 | WRITE (nmessg, 6005) |
---|
1826 | 6001 FORMAT ('Hines GW (east-west) at longitude I =', I3) |
---|
1827 | 6005 FORMAT (15X, ' U ', 2X, 'sig_E', 2X, 'sig_T', 3X, 'm_E', 4X, 'm_W', 4X, & |
---|
1828 | 'fluxU', 5X, 'gwdU') |
---|
1829 | DO l = levprt1, levprt2 |
---|
1830 | WRITE (nmessg, 6701) alt(i, l)/1.E3, v_alpha(i, l, n_east), & |
---|
1831 | sigma_alpha(i, l, n_east), sigma_t(i, l), & |
---|
1832 | m_alpha(i, l, n_east)*1.E3, m_alpha(i, l, n_west)*1.E3, & |
---|
1833 | flux_u(i, l)*1.E5, drag_u(i, l)*24.*3600. |
---|
1834 | END DO |
---|
1835 | 6701 FORMAT (' z=', F7.2, 1X, 3F7.1, 2F7.3, F9.4, F9.3) |
---|
1836 | END IF |
---|
1837 | |
---|
1838 | ! Print north-south winds, sigmas, cutoff wavenumbers, flux and drag. |
---|
1839 | |
---|
1840 | IF (iv_print==1) THEN |
---|
1841 | WRITE (nmessg, *) |
---|
1842 | WRITE (nmessg, 6002) i |
---|
1843 | 6002 FORMAT ('Hines GW (north-south) at longitude I =', I3) |
---|
1844 | WRITE (nmessg, 6006) |
---|
1845 | 6006 FORMAT (15X, ' V ', 2X, 'sig_N', 2X, 'sig_T', 3X, 'm_N', 4X, 'm_S', 4X, & |
---|
1846 | 'fluxV', 5X, 'gwdV') |
---|
1847 | DO l = levprt1, levprt2 |
---|
1848 | WRITE (nmessg, 6701) alt(i, l)/1.E3, v_alpha(i, l, n_north), & |
---|
1849 | sigma_alpha(i, l, n_north), sigma_t(i, l), & |
---|
1850 | m_alpha(i, l, n_north)*1.E3, m_alpha(i, l, n_south)*1.E3, & |
---|
1851 | flux_v(i, l)*1.E5, drag_v(i, l)*24.*3600. |
---|
1852 | END DO |
---|
1853 | END IF |
---|
1854 | |
---|
1855 | END DO |
---|
1856 | |
---|
1857 | RETURN |
---|
1858 | ! ----------------------------------------------------------------------- |
---|
1859 | END SUBROUTINE hines_print |
---|
1860 | |
---|
1861 | SUBROUTINE hines_exp(data, data_zmax, alt, alt_exp, iorder, il1, il2, lev1, & |
---|
1862 | lev2, nlons, nlevs) |
---|
1863 | IMPLICIT NONE |
---|
1864 | ! This routine exponentially damps a longitude by altitude array |
---|
1865 | ! of data above a specified altitude. |
---|
1866 | |
---|
1867 | ! Aug. 13/95 - C. McLandress |
---|
1868 | |
---|
1869 | ! Output arguements: |
---|
1870 | |
---|
1871 | ! * DATA = modified data array. |
---|
1872 | |
---|
1873 | ! Input arguements: |
---|
1874 | |
---|
1875 | ! * DATA = original data array. |
---|
1876 | ! * ALT = altitudes. |
---|
1877 | ! * ALT_EXP = altitude above which exponential decay applied. |
---|
1878 | ! * IORDER = 1 means vertical levels are indexed from top down |
---|
1879 | ! * (i.e., highest level indexed 1 and lowest level NLEVS); |
---|
1880 | ! * .NE. 1 highest level is index NLEVS. |
---|
1881 | ! * IL1 = first longitudinal index to use (IL1 >= 1). |
---|
1882 | ! * IL2 = last longitudinal index to use (IL1 <= IL2 <= NLONS). |
---|
1883 | ! * LEV1 = first altitude level to use (LEV1 >=1). |
---|
1884 | ! * LEV2 = last altitude level to use (LEV1 < LEV2 <= NLEVS). |
---|
1885 | ! * NLONS = number of longitudes. |
---|
1886 | ! * NLEVS = number of vertical |
---|
1887 | |
---|
1888 | ! Input work arrays: |
---|
1889 | |
---|
1890 | ! * DATA_ZMAX = data values just above altitude ALT_EXP. |
---|
1891 | |
---|
1892 | INTEGER iorder, il1, il2, lev1, lev2, nlons, nlevs |
---|
1893 | REAL alt_exp |
---|
1894 | REAL data(nlons, nlevs), data_zmax(nlons), alt(nlons, nlevs) |
---|
1895 | |
---|
1896 | ! Internal variables. |
---|
1897 | |
---|
1898 | INTEGER levbot, levtop, lincr, i, l |
---|
1899 | REAL hscale |
---|
1900 | DATA hscale/5.E3/ |
---|
1901 | ! ----------------------------------------------------------------------- |
---|
1902 | |
---|
1903 | ! Index of lowest altitude level (bottom of drag calculation). |
---|
1904 | |
---|
1905 | levbot = lev2 |
---|
1906 | levtop = lev1 |
---|
1907 | lincr = 1 |
---|
1908 | IF (iorder/=1) THEN |
---|
1909 | levbot = lev1 |
---|
1910 | levtop = lev2 |
---|
1911 | lincr = -1 |
---|
1912 | END IF |
---|
1913 | |
---|
1914 | ! Data values at first level above ALT_EXP. |
---|
1915 | |
---|
1916 | DO i = il1, il2 |
---|
1917 | DO l = levtop, levbot, lincr |
---|
1918 | IF (alt(i,l)>=alt_exp) THEN |
---|
1919 | data_zmax(i) = data(i, l) |
---|
1920 | END IF |
---|
1921 | END DO |
---|
1922 | END DO |
---|
1923 | |
---|
1924 | ! Exponentially damp field above ALT_EXP to model top at L=1. |
---|
1925 | |
---|
1926 | DO l = 1, lev2 |
---|
1927 | DO i = il1, il2 |
---|
1928 | IF (alt(i,l)>=alt_exp) THEN |
---|
1929 | data(i, l) = data_zmax(i)*exp((alt_exp-alt(i,l))/hscale) |
---|
1930 | END IF |
---|
1931 | END DO |
---|
1932 | END DO |
---|
1933 | |
---|
1934 | RETURN |
---|
1935 | ! ----------------------------------------------------------------------- |
---|
1936 | END SUBROUTINE hines_exp |
---|
1937 | |
---|
1938 | SUBROUTINE vert_smooth(data, work, coeff, nsmooth, il1, il2, lev1, lev2, & |
---|
1939 | nlons, nlevs) |
---|
1940 | IMPLICIT NONE |
---|
1941 | ! Smooth a longitude by altitude array in the vertical over a |
---|
1942 | ! specified number of levels using a three point smoother. |
---|
1943 | |
---|
1944 | ! NOTE: input array DATA is modified on output! |
---|
1945 | |
---|
1946 | ! Aug. 3/95 - C. McLandress |
---|
1947 | |
---|
1948 | ! Output arguement: |
---|
1949 | |
---|
1950 | ! * DATA = smoothed array (on output). |
---|
1951 | |
---|
1952 | ! Input arguements: |
---|
1953 | |
---|
1954 | ! * DATA = unsmoothed array of data (on input). |
---|
1955 | ! * WORK = work array of same dimension as DATA. |
---|
1956 | ! * COEFF = smoothing coefficient for a 1:COEFF:1 stencil. |
---|
1957 | ! * (e.g., COEFF = 2 will result in a smoother which |
---|
1958 | ! * weights the level L gridpoint by two and the two |
---|
1959 | ! * adjecent levels (L+1 and L-1) by one). |
---|
1960 | ! * NSMOOTH = number of times to smooth in vertical. |
---|
1961 | ! * (e.g., NSMOOTH=1 means smoothed only once, |
---|
1962 | ! * NSMOOTH=2 means smoothing repeated twice, etc.) |
---|
1963 | ! * IL1 = first longitudinal index to use (IL1 >= 1). |
---|
1964 | ! * IL2 = last longitudinal index to use (IL1 <= IL2 <= NLONS). |
---|
1965 | ! * LEV1 = first altitude level to use (LEV1 >=1). |
---|
1966 | ! * LEV2 = last altitude level to use (LEV1 < LEV2 <= NLEVS). |
---|
1967 | ! * NLONS = number of longitudes. |
---|
1968 | ! * NLEVS = number of vertical levels. |
---|
1969 | |
---|
1970 | ! Subroutine arguements. |
---|
1971 | |
---|
1972 | INTEGER nsmooth, il1, il2, lev1, lev2, nlons, nlevs |
---|
1973 | REAL coeff |
---|
1974 | REAL data(nlons, nlevs), work(nlons, nlevs) |
---|
1975 | |
---|
1976 | ! Internal variables. |
---|
1977 | |
---|
1978 | INTEGER i, l, ns, lev1p, lev2m |
---|
1979 | REAL sum_wts |
---|
1980 | ! ----------------------------------------------------------------------- |
---|
1981 | |
---|
1982 | ! Calculate sum of weights. |
---|
1983 | |
---|
1984 | sum_wts = coeff + 2. |
---|
1985 | |
---|
1986 | lev1p = lev1 + 1 |
---|
1987 | lev2m = lev2 - 1 |
---|
1988 | |
---|
1989 | ! Smooth NSMOOTH times |
---|
1990 | |
---|
1991 | DO ns = 1, nsmooth |
---|
1992 | |
---|
1993 | ! Copy data into work array. |
---|
1994 | |
---|
1995 | DO l = lev1, lev2 |
---|
1996 | DO i = il1, il2 |
---|
1997 | work(i, l) = data(i, l) |
---|
1998 | END DO |
---|
1999 | END DO |
---|
2000 | |
---|
2001 | ! Smooth array WORK in vertical direction and put into DATA. |
---|
2002 | |
---|
2003 | DO l = lev1p, lev2m |
---|
2004 | DO i = il1, il2 |
---|
2005 | data(i, l) = (work(i,l+1)+coeff*work(i,l)+work(i,l-1))/sum_wts |
---|
2006 | END DO |
---|
2007 | END DO |
---|
2008 | |
---|
2009 | END DO |
---|
2010 | |
---|
2011 | RETURN |
---|
2012 | ! ----------------------------------------------------------------------- |
---|
2013 | END SUBROUTINE vert_smooth |
---|
2014 | |
---|
2015 | |
---|
2016 | |
---|
2017 | |
---|
2018 | |
---|
2019 | |
---|