[1992] | 1 | |
---|
[1279] | 2 | ! $Id: hines_gwd.F90 3102 2017-12-03 20:27:42Z lebasn $ |
---|
[1001] | 3 | |
---|
[1992] | 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) |
---|
[1001] | 6 | |
---|
[1992] | 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 | ! ######################################################################## |
---|
[1001] | 13 | |
---|
| 14 | |
---|
[1992] | 15 | USE dimphy |
---|
| 16 | IMPLICIT NONE |
---|
[1001] | 17 | |
---|
[1992] | 18 | include "YOEGWD.h" |
---|
| 19 | include "YOMCST.h" |
---|
[1001] | 20 | |
---|
[1992] | 21 | INTEGER nazmth |
---|
| 22 | PARAMETER (nazmth=8) |
---|
[1001] | 23 | |
---|
[1992] | 24 | ! INPUT ARGUMENTS. |
---|
| 25 | ! ----- ---------- |
---|
[1001] | 26 | |
---|
[1992] | 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) |
---|
[1001] | 33 | |
---|
| 34 | |
---|
[1992] | 35 | ! REFERENCE. |
---|
| 36 | ! ---------- |
---|
| 37 | ! SEE MODEL DOCUMENTATION |
---|
[1001] | 38 | |
---|
[1992] | 39 | ! AUTHOR. |
---|
| 40 | ! ------- |
---|
[1001] | 41 | |
---|
[1992] | 42 | ! N. MCFARLANE DKRZ-HAMBURG MAY 1995 |
---|
| 43 | ! STAND ALONE E. MANZINI MPI-HAMBURG FEBRUARY 1997 |
---|
[1001] | 44 | |
---|
[1992] | 45 | ! BASED ON A COMBINATION OF THE OROGRAPHIC SCHEME BY N.MCFARLANE 1987 |
---|
| 46 | ! AND THE HINES SCHEME AS CODED BY C. MCLANDRESS 1995. |
---|
[1001] | 47 | |
---|
| 48 | |
---|
| 49 | |
---|
[1992] | 50 | ! ym INTEGER KLEVM1 |
---|
[1001] | 51 | |
---|
[1992] | 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) |
---|
[1001] | 59 | |
---|
[1992] | 60 | REAL th(klon, klev), utendgw(klon, klev), vtendgw(klon, klev), & |
---|
| 61 | pressg(klon), uhs(klon, klev), vhs(klon, klev), zpr(klon) |
---|
[1001] | 62 | |
---|
[1992] | 63 | ! * VERTICAL POSITIONING ARRAYS. |
---|
[1001] | 64 | |
---|
[1992] | 65 | REAL sgj(klon, klev), shj(klon, klev), shxkj(klon, klev), dsgj(klon, klev) |
---|
[1001] | 66 | |
---|
[1992] | 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 |
---|
[1001] | 70 | |
---|
| 71 | |
---|
[1992] | 72 | ! * WORK ARRAYS. |
---|
[1001] | 73 | |
---|
[1992] | 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) |
---|
[1001] | 86 | |
---|
[1992] | 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. |
---|
[1001] | 93 | |
---|
| 94 | |
---|
[1992] | 95 | REAL rmscon |
---|
| 96 | INTEGER nmessg, iprint, ilrms |
---|
| 97 | INTEGER ifl |
---|
[1001] | 98 | |
---|
[1992] | 99 | INTEGER naz, icutoff, nsmax, iheatcal |
---|
| 100 | REAL slope, f1, f2, f3, f5, f6, kstar(klon), alt_cutoff, smco |
---|
[1001] | 101 | |
---|
[1992] | 102 | ! PROVIDED AS INPUT |
---|
[1001] | 103 | |
---|
[1992] | 104 | INTEGER nlon, nlev |
---|
[1001] | 105 | |
---|
[1992] | 106 | REAL dtime |
---|
| 107 | REAL paphm1x(nlon, nlev+1), papm1x(nlon, nlev) |
---|
| 108 | REAL ux(nlon, nlev), vx(nlon, nlev), tx(nlon, nlev) |
---|
[1001] | 109 | |
---|
[1992] | 110 | ! VARIABLES FOR OUTPUT |
---|
[1001] | 111 | |
---|
| 112 | |
---|
[1992] | 113 | REAL d_t_hin(nlon, nlev), d_u_hin(nlon, nlev), d_v_hin(nlon, nlev) |
---|
| 114 | REAL zustrhi(nlon), zvstrhi(nlon) |
---|
[1001] | 115 | |
---|
| 116 | |
---|
[1992] | 117 | ! * LOGICAL SWITCHES TO CONTROL PRECIP ENHANCEMENT AND |
---|
| 118 | ! * HINES' DOPPLER SPREADING EXTROWAVE GW DRAG. |
---|
| 119 | ! * LOZPR IS TRUE FOR ZPR ENHANCEMENT |
---|
[1001] | 120 | |
---|
[1992] | 121 | LOGICAL lozpr, lorms(klon) |
---|
[1001] | 122 | |
---|
[1992] | 123 | ! LOCAL PARAMETERS TO MAKE THINGS WORK (TEMPORARY VARIABLE) |
---|
[1001] | 124 | |
---|
[1992] | 125 | REAL rhoh2o, zpcons, rgocp, zlat, dttdsf, ratio, hscal |
---|
| 126 | INTEGER i, j, l, jl, jk, le, lref, lrefp, levbot |
---|
[1001] | 127 | |
---|
[1992] | 128 | ! DATA PARAMETERS NEEDED, EXPLAINED LATER |
---|
[1001] | 129 | |
---|
[1992] | 130 | REAL v0, vmin, dmpscal, taufac, hmin, apibt, cpart, fcrit |
---|
| 131 | REAL pcrit, pcons |
---|
| 132 | INTEGER iplev, ierror |
---|
[1001] | 133 | |
---|
| 134 | |
---|
| 135 | |
---|
[1992] | 136 | ! PRINT *,' IT IS STARTED HINES GOING ON...' |
---|
[1001] | 137 | |
---|
| 138 | |
---|
| 139 | |
---|
| 140 | |
---|
[1992] | 141 | ! * COMPUTATIONAL CONSTANTS. |
---|
| 142 | ! ------------- ---------- |
---|
[1001] | 143 | |
---|
| 144 | |
---|
[1992] | 145 | d_t_hin(:, :) = 0. |
---|
[1001] | 146 | |
---|
[1992] | 147 | rhoh2o = 1000. |
---|
| 148 | zpcons = (1000.*86400.)/rhoh2o |
---|
| 149 | ! ym KLEVM1=KLEV-1 |
---|
[1001] | 150 | |
---|
| 151 | |
---|
[1992] | 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 |
---|
[1001] | 163 | |
---|
[1992] | 164 | ! Define constants and arrays needed for the ccc/mam gwd scheme |
---|
| 165 | ! *Constants: |
---|
[1001] | 166 | |
---|
[1992] | 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) |
---|
[1001] | 307 | END IF |
---|
[1992] | 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) |
---|
[1001] | 311 | END IF |
---|
[1992] | 312 | END DO |
---|
| 313 | END DO |
---|
[1001] | 314 | |
---|
[1992] | 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. |
---|
[1001] | 319 | |
---|
[1992] | 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 |
---|
[1001] | 331 | |
---|
[1992] | 332 | ! * ALTITUDE AND DENSITY AT BOTTOM. |
---|
[1001] | 333 | |
---|
[1992] | 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 |
---|
[1001] | 383 | END IF |
---|
[1992] | 384 | END DO |
---|
| 385 | END IF |
---|
[1001] | 386 | |
---|
[1992] | 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 |
---|
[1001] | 393 | |
---|
[1992] | 394 | ! * CALCULATE GWD (NOTE THAT DIFFUSION COEFFICIENT AND |
---|
| 395 | ! * HEATING RATE ONLY CALCULATED IF IHEATCAL = 1). |
---|
[1001] | 396 | |
---|
[1992] | 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) |
---|
[2197] | 714 | IMPLICIT NONE |
---|
[1992] | 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 |
---|
[2197] | 765 | REAL slope, kstar(nlons), f1, f2, f3, f2mfac |
---|
[1992] | 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) |
---|
[1001] | 839 | END IF |
---|
[1992] | 840 | END DO |
---|
| 841 | END DO |
---|
[1001] | 842 | |
---|
[1992] | 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 |
---|
[1001] | 872 | END IF |
---|
| 873 | END IF |
---|
[1992] | 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 | |
---|
[1001] | 901 | END IF |
---|
[1992] | 902 | END DO |
---|
| 903 | END DO |
---|
[1001] | 904 | |
---|
[1992] | 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) |
---|
[2197] | 937 | IMPLICIT NONE |
---|
[1992] | 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) |
---|
[2197] | 1033 | IMPLICIT NONE |
---|
[1992] | 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 | |
---|
[2197] | 1088 | INTEGER i, l, lev1p, lev2m, lev2p |
---|
[1992] | 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 | |
---|
[1001] | 1196 | END IF |
---|
[1992] | 1197 | END DO |
---|
| 1198 | END DO |
---|
[1001] | 1199 | |
---|
[1992] | 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 |
---|
[1001] | 1222 | END IF |
---|
[1992] | 1223 | END DO |
---|
| 1224 | END IF |
---|
[1001] | 1225 | |
---|
[1992] | 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) |
---|
[2197] | 1233 | IMPLICIT NONE |
---|
[1992] | 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) |
---|
[2197] | 1354 | IMPLICIT NONE |
---|
[1992] | 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) |
---|
[2197] | 1449 | IMPLICIT NONE |
---|
[1992] | 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. |
---|
[1001] | 1523 | ELSE |
---|
[1992] | 1524 | i_alpha(i, n) = (qm**2/2.+qm**3/3.+qm**4/4.+qm**5/5.)/ & |
---|
| 1525 | q_alpha**2 |
---|
[1001] | 1526 | END IF |
---|
[1992] | 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. |
---|
[1001] | 1556 | ELSE |
---|
[1992] | 1557 | i_alpha(i, n) = (qm**3/3.+qm**4/4.+qm**5/5.+qm**6/6.)/ & |
---|
| 1558 | q_alpha**3 |
---|
[1001] | 1559 | END IF |
---|
[1992] | 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 |
---|
[1001] | 1589 | ELSE |
---|
[1992] | 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 |
---|
[1001] | 1592 | END IF |
---|
[1992] | 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 |
---|
[1001] | 1604 | |
---|
[1992] | 1605 | END IF |
---|
| 1606 | END DO |
---|
| 1607 | END DO |
---|
[1001] | 1608 | |
---|
[1992] | 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) |
---|
[2197] | 1642 | IMPLICIT NONE |
---|
[1992] | 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) |
---|
[2197] | 1771 | IMPLICIT NONE |
---|
[1992] | 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) |
---|
[2197] | 1863 | IMPLICIT NONE |
---|
[1992] | 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) |
---|
[1001] | 1920 | END IF |
---|
[1992] | 1921 | END DO |
---|
| 1922 | END DO |
---|
[1001] | 1923 | |
---|
[1992] | 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) |
---|
[1001] | 1930 | END IF |
---|
[1992] | 1931 | END DO |
---|
| 1932 | END DO |
---|
[1001] | 1933 | |
---|
[1992] | 1934 | RETURN |
---|
| 1935 | ! ----------------------------------------------------------------------- |
---|
| 1936 | END SUBROUTINE hines_exp |
---|
[1001] | 1937 | |
---|
[1992] | 1938 | SUBROUTINE vert_smooth(data, work, coeff, nsmooth, il1, il2, lev1, lev2, & |
---|
| 1939 | nlons, nlevs) |
---|
[2197] | 1940 | IMPLICIT NONE |
---|
[1992] | 1941 | ! Smooth a longitude by altitude array in the vertical over a |
---|
| 1942 | ! specified number of levels using a three point smoother. |
---|
[1001] | 1943 | |
---|
[1992] | 1944 | ! NOTE: input array DATA is modified on output! |
---|
[1001] | 1945 | |
---|
[1992] | 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 | |
---|