1 !----------------------------------------------------------------------- 2 SUBROUTINE WAKE_DFSN(LDBHR,XI,SZI,SYI,Z) 3 !----------------------------------------------------------------------- 4 ! 5 ! --- PRIME Version: 1.0 Level: 970812 WAKE_DFSN 6 ! L. Schulman, D. Strimaitis, EARTH TECH 7 ! Prepared for EPRI under contract WO3527-01 8 ! 9 ! --- PURPOSE: Tabulates sigmas and rate of plume growth as function 10 ! of location within the wake from modified Weil (1996) 11 ! analytical expressions 12 ! 13 ! --- MODIFIED: To modify calling arguments in call to subroutine 14 ! WAKE_SIG from turby and turbz to wakiy and wakiz. 15 ! R.W. Brode, PES, Inc. - 07/05/01 16 ! 17 ! --- MODIFIED: For use with the AERMOD model. Use virtual source sigmas 18 ! added in quadrature to ambient sigmas instead of 19 ! virtual source distances. 20 ! R.W. Brode, PES, Inc. - 07/05/01 21 ! 22 ! --- INPUTS: 23 ! ldbhr - logical - Flag for debug write statements 24 ! to upwind bldg wall 25 ! xi - real - distance (m) from upwind bldg wall 26 ! to point where plume intersects wake 27 ! szi - real - sigma-z (m) at xi 28 ! syi - real - sigma-y (m) at xi 29 ! z - real - plume height (m) at xi 30 ! 31 ! Common block /PARAMS/ variables: 32 ! MXNTR, MXNW 33 ! Common block /WAKEDAT/ variables: 34 ! XBADJ, Hb, Wb, xLb, Rb, xLR 35 ! 36 ! --- OUTPUT: 37 ! 38 ! Common block /WAKEDAT/ variables: 39 ! NWAK, XWAK(mxntr), SZWAK(mxntr), SYWAK(mxntr), 40 ! DRWAK(mxntr), XZVWAK, XYVWAK, 41 ! NCAV, XCAV(mxntr), SZCAV(mxntr), SYCAV(mxntr), 42 ! XZVCAV, XYVCAV, LRURL, ISTAB 43 ! 44 ! --- WAKE_DFSN called by: NUMRISE 45 ! --- WAKE_DFSN calls : SIGZPR, SIGYPR, 46 ! WAKE_XA, WAKE_CAV0, WAKE_TURB, WAKE_SIG 47 !---------------------------------------------------------------------- 48 ! 49 INCLUDE 'params.pri' 50 INCLUDE 'numparm.pri' 51 INCLUDE 'wakedat.pri' 52 53 ! --- Define local variable arrays for fine-steps 54 REAL DIST(MXNW) , ASIGZ(MXNW) , ASIGY(MXNW) , DSZ(MXNW) 55 REAL CSIGZ(MXNW) , CSIGY(MXNW) 56 57 LOGICAL LDBHR , LCAV , LWAK , LREVCAV 58 !JRA 30 Sept 2005 - following variables were used without being defined 59 ! this could lead to erroneous results 60 ! spotted by Salford FTN95 version 4.8.0 61 REAL :: YKDUM=0.0 , ZKDUM=0.0 62 63 ! --- Misc. constants 64 DATA ZERO/0.0/ , HALF/0.5/ , ONE/1.0/ 65 DATA RTPIBY2/1.253314/ 66 67 ! --- Compute xa, where turbulent growth in the wake transitions 68 ! --- to ambient growth rate, measured from upwind face of bldg 69 CALL WAKE_XA(ISTAB,LRURL,XLB,RB,XAZ,XAY) ! 5466 70 XAMX = AMAX1(XAZ,XAY) 71 XAMN = AMIN1(XAZ,XAY) 72 73 ! --- Initialize virtual source sigma terms 74 VSIGY = 0.0 75 VSIGZ = 0.0 76 VSIGYC = 0.0 77 VSIGZC = 0.0 78 79 ! --- Initialize CAVITY parameters 80 ! -------------------------------- 81 ! --- Set distance from upwind face of bldg to END of cavity 82 XCAVE = XLB + XLR 83 ! --- Set distance from source to start of cavity 84 DISTC = XLB + XBADJ 85 ! --- Set downwind distance to effective cavity source (when present), 86 ! --- from the upwind face of bldg 87 XBC = AMAX1(XI,XLB) 88 XBC = AMIN1(XBC,XCAVE) 89 ! --- Location of downwind edge of PDF region from effective cavity 90 ! --- source 91 XDC = XBC + XLR 92 ! --- Set initial sigmas for cavity source using sigma-y at xi 93 CALL WAKE_CAV0(SYI,SZCAV0,SYCAV0) 94 ! --- The cavity sigma-y will need to be revised if xi lies upwind of 95 ! --- the downwind face of the bldg. 96 IF ( XI.LT.XLB ) THEN 97 LREVCAV = .TRUE. ! 3710 98 ELSE 99 LREVCAV = .FALSE. ! 1756 100 ENDIF 101 102 ! --- Determine if any plume material in cavity may be modeled 103 ! ------------------------------------------------------------ 104 ! --- Initialize output arrays 105 NCAV = 1 ! 5466 106 XCAV(1) = XBC + XBADJ 107 SZCAV(1) = SZCAV0 108 SYCAV(1) = SYCAV0 109 IF ( XI.GE.XCAVE ) THEN 110 LCAV = .FALSE. ! 780 111 LREVCAV = .FALSE. 112 ELSE 113 LCAV = .TRUE. ! 4686 114 ENDIF 115 116 ! --- Is plume affected by wake turbulence? 117 ! ------------------------------------------- 118 ! --- Initialize output arrays 119 NWAK = 1 ! 5466 120 XWAK(1) = XI + XBADJ 121 SZWAK(1) = SZI 122 SYWAK(1) = SYI 123 DRWAK(1) = ZERO 124 IF ( XI.GE.XAMX ) THEN 125 LWAK = .FALSE. ! 92 126 IF ( LDBHR ) THEN 127 WRITE (IO6,*) '----- WAKE_DFSN: NWAK = ' , NWAK ! 0 128 WRITE (IO6,*) 'Z-dispersion reaches ambient at: ' , & 129 & (XAZ+XBADJ) 130 WRITE (IO6,*) 'Y-dispersion reaches ambient at: ' , & 131 & (XAY+XBADJ) 132 !aer write(io6,*)'z,y virtual distances (m) = ',xzvwak,xyvwak 133 WRITE (IO6,*) 'xadj, yadj, xi (m) = ' , XBADJ , & 134 & YBADJ , XI 135 WRITE (IO6,*) 'Plume NOT altered by wake turbulence!' 136 WRITE (IO6,*) 137 ENDIF 138 ELSE 139 LWAK = .TRUE. ! 5374 140 ENDIF 141 142 ! --- Return now if sigmas in wake do not need to be tabulated 143 IF ( .NOT.LWAK .AND. .NOT.LCAV ) RETURN ! 5466 144 145 ! --- Compute location of downwind edge of PDF region from xi 146 XD = XI + XLR ! 5374 147 148 ! --- Set stepping parameters 149 DX = 2.0 150 ! --- Range of table is from point of entry into wake (xi), to the point 151 ! --- at which ambient growth rate resumes (xa), plus one "ds" so that 152 ! --- both sigmas reach ambient, and virtual distances are computed. 153 ! --- When cavity sigmas are also computed, range may start at the 154 ! --- downwind bldg face, and extend to end of cavity. 155 XLOW = XI 156 XHI = XAMX 157 IF ( LCAV ) THEN 158 XLOW = AMIN1(XI,XBC) ! 4686 159 XHI = AMAX1(XAMX,XCAVE) 160 ENDIF 161 XRANGE = XHI - XLOW + DX ! 5374 162 NP = NINT(XRANGE/DX) + 1 163 NP = MIN(NP,MXNW-1) 164 DX = XRANGE/(FLOAT(NP)-ONE) 165 DXI = ONE/DX 166 NWS = 0 167 NCS = 0 168 169 ! --- Fill first element of marching arrays using values at xlow 170 DIST(1) = XLOW + XBADJ 171 IF ( LWAK ) THEN 172 ASIGZ(1) = SZI ! 5374 173 ASIGY(1) = SYI 174 ! --- Set inital plume growth rate in wake to zero 175 DSZ(1) = ZERO 176 ENDIF 177 IF ( LCAV ) THEN ! 5374 178 CSIGZ(1) = SZCAV0 ! 4686 179 CSIGY(1) = SYCAV0 180 ENDIF 181 182 ! --- Initialize distance (from upwind face of bldg) 183 X = XLOW ! 5374 184 185 ! --- Loop over steps in wake region 186 ! ----------------------------------- 187 DO N = 2 , NP 188 XOLD = X !1921588 189 X = X + DX 190 DIST(N) = DIST(N-1) + DX 191 192 ! --- Check to see if cavity data should be revised based on 193 ! --- data from previous step 194 IF ( LREVCAV .AND. XOLD.GE.XLB ) THEN 195 CALL WAKE_CAV0(ASIGY(N-1),SZCAV0,SYCAV0R) ! 3710 196 IF ( SYCAV0R.GT.SYCAV0 ) THEN 197 SYCAV0 = SYCAV0R ! 0 198 SYCAV(1) = SYCAV0 199 ! --- Replace sigma-y values in stepping arrays 200 DO IR = 1 , N - 1 201 CSIGY(IR) = AMAX1(CSIGY(IR),SYCAV0) ! 0 202 ENDDO 203 ENDIF 204 LREVCAV = .FALSE. ! 3710 205 ENDIF 206 207 ! --- Obtain sigmas for this step 208 209 ! --- First, persist initial values if upwind of starting point 210 IF ( LWAK .AND. (XI.GE.X) ) THEN !1921588 211 ASIGZ(N) = ASIGZ(N-1) ! 0 212 ASIGY(N) = ASIGY(N-1) 213 DSZ(N) = DSZ(N-1) 214 ! --- Set index for skipping entry when filling wake arrays 215 NWS = N 216 ENDIF 217 IF ( LCAV .AND. (XBC.GE.X) ) THEN !1921588 218 CSIGZ(N) = SZCAV0 ! 72000 219 CSIGY(N) = SYCAV0 220 ! --- Set index for skipping entry when filling cav arrays 221 NCS = N 222 ENDIF 223 224 ! --- Now test again and apply full treatment when needed 225 IF ( XOLD.GT.XAMX ) THEN !1921588 226 ! --- Ambient growth region in wake: use virtuals 227 IF ( LWAK .AND. (XI.LT.X) ) THEN ! 5034 228 VSIGZ = MAX(VSIGZ,SZI) ! 5034 229 VSIGY = MAX(VSIGY,SYI) 230 CALL SIGZPR(DIST(N),Z,ASIGZ(N)) 231 CALL SIGYPR(DIST(N),Z,ASIGY(N)) 232 ASIGZ(N) = SQRT(ASIGZ(N)**2+VSIGZ**2) 233 ASIGY(N) = SQRT(ASIGY(N)**2+VSIGY**2) 234 DSZ(N) = (ASIGZ(N)-ASIGZ(N-1))*DXI 235 ENDIF 236 ! --- Cavity source --- 237 IF ( LCAV .AND. (XBC.LT.X) ) THEN ! 5034 238 VSIGZC = MAX(VSIGZC,SZCAV0) ! 4690 239 VSIGYC = MAX(VSIGYC,SYCAV0) 240 CALL SIGZPR(DIST(N),0.0,CSIGZ(N)) 241 CALL SIGYPR(DIST(N),0.0,CSIGY(N)) 242 CSIGZ(N) = SQRT(CSIGZ(N)**2+VSIGZC**2) 243 CSIGY(N) = SQRT(CSIGY(N)**2+VSIGYC**2) 244 ENDIF 245 ELSE 246 IF ( X.LT.XAMN ) THEN !1916554 247 ! --- Wake growth for both sigz and sigy 248 ! --- Set x at mid-point of step 249 XMID = HALF*(X+XOLD) !1755956 250 ! --- Compute turbulence intensities at midpoint 251 CALL WAKE_TURB(ISTAB,LRURL,XMID,XLB,RB,WAKIZ,WAKIY) 252 ! --- Compute sigmas in wake 253 IF ( LWAK .AND. (XI.LE.X) ) & 254 & CALL WAKE_SIG(X,XD,XOLD,WAKIZ,WAKIY,ASIGZ(N-1), & 255 & ASIGY(N-1),HB,WB,RB,ZK,YK,ASIGZ(N),ASIGY(N),DSZ(N)) 256 ! --- Cavity source --- 257 IF ( LCAV .AND. (XBC.LE.X) ) & 258 & CALL WAKE_SIG(X,XDC,XOLD,WAKIZ,WAKIY,CSIGZ(N-1), & 259 & CSIGY(N-1),HB,WB,RB,ZKC,YKC,CSIGZ(N),CSIGY(N), & 260 & DZRATE) 261 ELSE 262 ! --- At least one of the sigmas reaches ambient growth in wake 263 ! --- Process SIGMA-Z 264 IF ( XOLD.GE.XAZ ) THEN ! 160598 265 ! --- Ambient growth region in wake: use virtual x 266 IF ( LWAK .AND. (XI.LE.X) ) THEN ! 358 267 CALL SIGZPR(DIST(N),Z,ASIGZ(N)) ! 358 268 VSIGZ = MAX(VSIGZ,SZI) 269 ASIGZ(N) = SQRT(ASIGZ(N)**2+VSIGZ**2) 270 DSZ(N) = (ASIGZ(N)-ASIGZ(N-1))*DXI 271 ENDIF 272 ! --- Cavity source --- 273 IF ( LCAV .AND. (XBC.LE.X) ) THEN ! 358 274 CALL SIGZPR(DIST(N),0.0,CSIGZ(N)) ! 358 275 CSIGZ(N) = SQRT(CSIGZ(N)**2+VSIGZC**2) 276 ENDIF 277 ELSEIF ( X.GE.XAZ ) THEN 278 ! --- Transition from wake to ambient 279 XNEW = XAZ ! 5372 280 XMID = HALF*(XNEW+XOLD) 281 ! --- Compute turbulence intensities at midpoint 282 CALL WAKE_TURB(ISTAB,LRURL,XMID,XLB,RB,WAKIZ,WAKIY) 283 IF ( LWAK .AND. (XI.LE.XNEW) ) THEN 284 ! --- Compute wake sigma at xaz 285 CALL WAKE_SIG(XNEW,XD,XOLD,WAKIZ,WAKIY,ASIGZ(N-1), & 286 & ASIGY(N-1),HB,WB,RB,ZK,YKDUM,SIGZXA, & 287 & SYDUM,DZRATE) 288 ! --- Get virtual source term as difference in quadrature between 289 ! --- wake and ambient sigmas at transition distance 290 CALL SIGZPR(XAZ+XBADJ,Z,SZ) 291 IF ( SIGZXA.GT.SZ ) THEN 292 VSIGZ = SQRT(SIGZXA**2-SZ**2) ! 5372 293 ELSE 294 VSIGZ = 0.0 ! 0 295 ENDIF 296 ! --- Now compute sigma at dist(n) with virtual source 297 CALL SIGZPR(DIST(N),Z,ASIGZ(N)) ! 5372 298 ASIGZ(N) = SQRT(ASIGZ(N)**2+VSIGZ**2) 299 DSZ(N) = (ASIGZ(N)-ASIGZ(N-1))*DXI 300 ENDIF 301 ! --- Cavity source --- 302 IF ( LCAV .AND. (XBC.LE.XNEW) ) THEN ! 5372 303 ! --- Compute wake sigma at xaz 304 CALL WAKE_SIG(XNEW,XDC,XOLD,WAKIZ,WAKIY,CSIGZ(N-1),& 305 & CSIGY(N-1),HB,WB,RB,ZKC,YKDUM,SIGZXA,& 306 & SYDUM,DZRATE) 307 ! --- Get virtual source term as difference in quadrature between 308 ! --- wake and ambient sigmas at transition distance 309 CALL SIGZPR(XAZ+XBADJ,0.0,SZ) 310 IF ( SIGZXA.GT.SZ ) THEN 311 VSIGZC = SQRT(SIGZXA**2-SZ**2) ! 4684 312 ELSE 313 VSIGZC = 0.0 ! 0 314 ENDIF 315 ! --- Now compute sigma at dist(n) with virtual source 316 CALL SIGZPR(DIST(N),0.0,CSIGZ(N)) ! 4684 317 CSIGZ(N) = SQRT(CSIGZ(N)**2+VSIGZC**2) 318 ENDIF 319 ELSE 320 ! --- Wake growth for sigz 321 ! --- Set x at mid-point of step 322 XMID = HALF*(X+XOLD) ! 154868 323 ! --- Compute turbulence intensities at midpoint 324 CALL WAKE_TURB(ISTAB,LRURL,XMID,XLB,RB,WAKIZ,WAKIY) 325 ! --- Compute sigmaz 326 IF ( LWAK .AND. (XI.LE.X) ) & 327 & CALL WAKE_SIG(X,XD,XOLD,WAKIZ,WAKIY,ASIGZ(N-1), & 328 & ASIGY(N-1),HB,WB,RB,ZK,YKDUM,ASIGZ(N),SYDUM, & 329 & DSZ(N)) 330 ! --- Cavity source --- 331 IF ( LCAV .AND. (XBC.LE.X) ) & 332 & CALL WAKE_SIG(X,XDC,XOLD,WAKIZ,WAKIY,CSIGZ(N-1), & 333 & CSIGY(N-1),HB,WB,RB,ZKC,YKDUM,CSIGZ(N),SYDUM, & 334 & DZRATE) 335 ENDIF 336 ! --- Process SIGMA-Y 337 IF ( XOLD.GE.XAY ) THEN ! 160598 338 ! --- Ambient growth region in wake: use virtual x 339 IF ( LWAK .AND. (XI.LE.X) ) THEN ! 154886 340 CALL SIGYPR(DIST(N),Z,ASIGY(N)) ! 154886 341 VSIGY = MAX(VSIGY,SYI) 342 ASIGY(N) = SQRT(ASIGY(N)**2+VSIGY**2) 343 ENDIF 344 ! --- Cavity source --- 345 IF ( LCAV .AND. (XBC.LE.X) ) THEN ! 154886 346 CALL SIGYPR(DIST(N),0.0,CSIGY(N)) ! 149308 347 VSIGYC = MAX(VSIGYC,SYCAV0) 348 CSIGY(N) = SQRT(CSIGY(N)**2+VSIGYC**2) 349 ENDIF 350 ELSEIF ( X.GE.XAY ) THEN 351 ! --- Transition from wake to ambient 352 XNEW = XAY ! 5372 353 XMID = HALF*(XNEW+XOLD) 354 ! --- Compute turbulence intensities at midpoint 355 CALL WAKE_TURB(ISTAB,LRURL,XMID,XLB,RB,WAKIZ,WAKIY) 356 IF ( LWAK .AND. (XI.LE.XNEW) ) THEN 357 ! --- Compute sigma at xay 358 !RWB call WAKE_SIG(xnew,xd,xold,turbz,turby,asigz(n-1), 359 !RWB turbz and turby appear to be the wrong variables for this 360 !RWB call to WAKE_SIG, try wakiz and wakiy 361 CALL WAKE_SIG(XNEW,XD,XOLD,WAKIZ,WAKIY,ASIGZ(N-1), & 362 & ASIGY(N-1),HB,WB,RB,ZKDUM,YK,SZDUM, & 363 & SIGYXA,DZRATE) 364 ! --- Get virtual source term as difference in quadrature between 365 ! --- wake and ambient sigmas at transition distance 366 CALL SIGYPR(XAY+XBADJ,Z,SY) 367 IF ( SIGYXA.GT.SY ) THEN 368 VSIGY = SQRT(SIGYXA**2-SY**2) ! 5328 369 ELSE 370 VSIGY = 0.0 ! 44 371 ENDIF 372 ! --- Now compute sigma at dist(n) with virtual source 373 CALL SIGYPR(DIST(N),Z,ASIGY(N)) ! 5372 374 ASIGY(N) = SQRT(ASIGY(N)**2+VSIGY**2) 375 ENDIF 376 ! --- Cavity source --- 377 IF ( LCAV .AND. (XBC.LE.XNEW) ) THEN ! 5372 378 CALL WAKE_SIG(XNEW,XDC,XOLD,WAKIZ,WAKIY,CSIGZ(N-1),& 379 & CSIGY(N-1),HB,WB,RB,ZKDUM,YKC,SZDUM, & 380 & SIGYXA,DZRATE) 381 CALL SIGYPR(XAY+XBADJ,0.0,SY) 382 IF ( SIGYXA.GT.SY ) THEN 383 VSIGYC = SQRT(SIGYXA**2-SY**2) ! 4686 384 ELSE 385 VSIGYC = 0.0 ! 0 386 ENDIF 387 CALL SIGYPR(DIST(N),0.0,CSIGY(N)) ! 4686 388 CSIGY(N) = SQRT(CSIGY(N)**2+VSIGYC**2) 389 ENDIF 390 ELSE 391 ! --- Wake growth for sigy 392 ! --- Set x at mid-point of step 393 XMID = HALF*(X+XOLD) ! 340 394 ! --- Compute turbulence intensities at midpoint 395 CALL WAKE_TURB(ISTAB,LRURL,XMID,XLB,RB,WAKIZ,WAKIY) 396 ! --- Compute sigmay 397 IF ( LWAK .AND. (XI.LE.X) ) & 398 & CALL WAKE_SIG(X,XD,XOLD,WAKIZ,WAKIY,ASIGZ(N-1), & 399 & ASIGY(N-1),HB,WB,RB,ZKDUM,YK,SZDUM,ASIGY(N), & 400 & DZRATE) 401 ! --- Cavity source 402 IF ( LCAV .AND. (XBC.LE.X) ) & 403 & CALL WAKE_SIG(X,XDC,XOLD,WAKIZ,WAKIY,CSIGZ(N-1), & 404 & CSIGY(N-1),HB,WB,RB,ZKDUM,YKC,SZDUM,CSIGY(N), & 405 & DZRATE) 406 ENDIF 407 ENDIF 408 ENDIF 409 410 ! --- Next distance 411 ENDDO 412 413 ! --- Construct arrays for /WAKEDAT/ 414 ! ---------------------------------- 415 416 IF ( LWAK ) THEN ! 5374 417 ! --- WAK arrays: 418 NPW = NP - NWS ! 5374 419 420 ! --- Place initial values into first element 421 XWAK(1) = XI + XBADJ 422 SZWAK(1) = SZI 423 SYWAK(1) = SYI 424 DRWAK(1) = ZERO 425 IF ( NPW.GE.MXNTR ) THEN 426 ! --- Sample a subset of the npw points 427 NWAK = MXNTR ! 5324 428 XWAK(NWAK) = DIST(NP) 429 SZWAK(NWAK) = ASIGZ(NP) 430 SYWAK(NWAK) = ASIGY(NP) 431 DRWAK(NWAK) = RTPIBY2*DSZ(NP) 432 IF ( NPW.LE.2*MXNTR ) THEN 433 ! --- Fill elements with nearest values 434 DELN = FLOAT(NPW)/FLOAT(NWAK) ! 150 435 DO IN = 2 , NWAK - 1 436 JN = IN*DELN + NWS ! 7200 437 XWAK(IN) = DIST(JN) 438 SZWAK(IN) = ASIGZ(JN) 439 SYWAK(IN) = ASIGY(JN) 440 DRWAK(IN) = RTPIBY2*DSZ(JN) 441 ENDDO 442 ELSE 443 ! --- Use sliding step-size to sample nearfield more frequently 444 DELN = 2.*FLOAT(NPW-MXNTR)/FLOAT(MXNTR*(MXNTR-1)) ! 5174 445 RN = ONE 446 DO IN = 2 , NWAK - 1 447 RN = RN + ONE + (IN-1)*DELN ! 248352 448 JN = RN + NWS 449 XWAK(IN) = DIST(JN) 450 SZWAK(IN) = ASIGZ(JN) 451 SYWAK(IN) = ASIGY(JN) 452 DRWAK(IN) = RTPIBY2*DSZ(JN) 453 ENDDO 454 ENDIF 455 ELSE 456 ! --- Fill only those elements used 457 NWAK = NPW ! 50 458 DO IN = 2 , NPW 459 INP = IN + NWS ! 1856 460 XWAK(IN) = DIST(INP) 461 SZWAK(IN) = ASIGZ(INP) 462 SYWAK(IN) = ASIGY(INP) 463 DRWAK(IN) = RTPIBY2*DSZ(INP) 464 ENDDO 465 ENDIF 466 ENDIF 467 468 IF ( LCAV ) THEN ! 5374 469 ! --- CAV arrays: 470 NPC = NP - NCS ! 4686 471 472 ! --- Place initial values into first element 473 XCAV(1) = XBC + XBADJ 474 SZCAV(1) = SZCAV0 475 SYCAV(1) = SYCAV0 476 IF ( NPC.GE.MXNTR ) THEN 477 ! --- Sample a subset of the npc points 478 NCAV = MXNTR ! 4606 479 XCAV(NCAV) = DIST(NP) 480 SZCAV(NCAV) = CSIGZ(NP) 481 SYCAV(NCAV) = CSIGY(NP) 482 IF ( NPC.LE.2*MXNTR ) THEN 483 ! --- Fill elements with nearest values 484 DELN = FLOAT(NPC)/FLOAT(NCAV) ! 116 485 DO IN = 2 , NCAV - 1 486 JN = IN*DELN + NCS ! 5568 487 XCAV(IN) = DIST(JN) 488 SZCAV(IN) = CSIGZ(JN) 489 SYCAV(IN) = CSIGY(JN) 490 ENDDO 491 ELSE 492 ! --- Use sliding step-size to sample nearfield more frequently 493 DELN = 2.*FLOAT(NPC-MXNTR)/FLOAT(MXNTR*(MXNTR-1)) ! 4490 494 RN = ONE 495 DO IN = 2 , NCAV - 1 496 RN = RN + ONE + (IN-1)*DELN ! 215520 497 JN = RN + NCS 498 XCAV(IN) = DIST(JN) 499 SZCAV(IN) = CSIGZ(JN) 500 SYCAV(IN) = CSIGY(JN) 501 ENDDO 502 ENDIF 503 ELSE 504 ! --- Fill only those elements used 505 NCAV = NPC ! 80 506 DO IN = 2 , NPC 507 INP = IN + NCS ! 3480 508 XCAV(IN) = DIST(INP) 509 SZCAV(IN) = CSIGZ(INP) 510 SYCAV(IN) = CSIGY(INP) 511 ENDDO 512 ENDIF 513 ENDIF 514 515 IF ( LDBHR ) THEN ! 5374 516 517 WRITE (IO6,*) ! 0 518 WRITE (IO6,*) '----- WAKE_DFSN: NWAK = ' , NWAK 519 WRITE (IO6,*) 'Z-dispersion reaches ambient at: ' , (XAZ+XBADJ) 520 WRITE (IO6,*) 'Y-dispersion reaches ambient at: ' , (XAY+XBADJ) 521 !aer write(io6,*)'z,y virtual dist (m) - WAKE = ',xzvwak,xyvwak 522 !aer write(io6,*)'z,y virtual dist (m) - CAV = ',xzvcav,xyvcav 523 WRITE (IO6,*) 'xadj, yadj, xi (m) = ' , XBADJ , & 524 & YBADJ , XI 525 WRITE (IO6,*) 'xbc,distc,xdc (m) = ' , XBC , DISTC ,& 526 & XDC 527 WRITE (IO6,*) 'lwak, nws, npw = ' , LWAK , NWS , & 528 & NPW 529 WRITE (IO6,*) 'lcav, ncs, npc = ' , LCAV , NCS , & 530 & NPC 531 WRITE (IO6,*) 532 ! 533 ! --- Write the arrays passed back to the calling routine 534 WRITE (IO6,28) 535 28 FORMAT (/4x,'I',9x,'XWAK',6x,'SZWAK',6x,'SYWAK',6x,'DRWAK',/) 536 DO I = 1 , NWAK 537 WRITE (IO6,32) I , XWAK(I) , SZWAK(I) , SYWAK(I) , DRWAK(I) ! 0 538 32 FORMAT (i5,3x,4(f10.4,1x)) 539 ENDDO 540 WRITE (IO6,*) ! 0 541 542 WRITE (IO6,29) 543 29 FORMAT (/4x,'I',9x,'XCAV',6x,'SZCAV',6x,'SYCAV',/) 544 DO I = 1 , NCAV 545 WRITE (IO6,33) I , XCAV(I) , SZCAV(I) , SYCAV(I) ! 0 546 33 FORMAT (i5,3x,3(f10.4,1x)) 547 ENDDO 548 WRITE (IO6,*) ! 0 549 ENDIF 550 551 CONTINUE ! 5374 552 END
HyperKWIC - Version 1.00DD executed at 20:00 on 1 Mar 2018 | Personal or Academic or Evaluation User | Free for Non-Commercial, Non-Government Use