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