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