Changeset 1320 for palm/trunk/SOURCE/diffusion_v.f90
- Timestamp:
- Mar 20, 2014 8:40:49 AM (10 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/diffusion_v.f90
r1310 r1320 20 20 ! Current revisions: 21 21 ! ----------------- 22 ! 22 ! ONLY-attribute added to USE-statements, 23 ! kind-parameters added to all INTEGER and REAL declaration statements, 24 ! kinds are defined in new module kinds, 25 ! old module precision_kind is removed, 26 ! revision history before 2012 removed, 27 ! comment fields (!:) to be used for variable explanations added to 28 ! all variable declaration statements 23 29 ! 24 30 ! Former revisions: … … 46 52 ! outflow damping layer removed 47 53 ! kmxm_x/_y and kmxp_x/_y change to kmxm and kmxp 48 !49 ! 667 2010-12-23 12:06:00Z suehring/gryschka50 ! nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng51 !52 ! 366 2009-08-25 08:06:27Z raasch53 ! bc_lr replaced by bc_lr_cyc54 !55 ! 106 2007-08-16 14:30:26Z raasch56 ! Momentumflux at top (vswst) included as boundary condition,57 ! j loop is starting from nysv (needed for non-cyclic boundary conditions)58 !59 ! 75 2007-03-22 09:54:05Z raasch60 ! Wall functions now include diabatic conditions, call of routine wall_fluxes,61 ! z0 removed from argument list, vynp eliminated62 !63 ! 20 2007-02-26 00:12:32Z raasch64 ! Bugfix: ddzw dimensioned 1:nzt"+1"65 !66 ! RCS Log replace by Id keyword, revision history cleaned up67 !68 ! Revision 1.15 2006/02/23 10:36:00 raasch69 ! nzb_2d replaced by nzb_v_outer in horizontal diffusion and by nzb_v_inner70 ! or nzb_diff_v, respectively, in vertical diffusion,71 ! wall functions added for north and south walls, +z0 in argument list,72 ! terms containing w(k-1,..) are removed from the Prandtl-layer equation73 ! because they cause errors at the edges of topography74 ! WARNING: loops containing the MAX function are still not properly vectorized!75 54 ! 76 55 ! Revision 1.1 1997/09/12 06:24:01 raasch … … 105 84 SUBROUTINE diffusion_v 106 85 107 USE arrays_3d 108 USE control_parameters 109 USE grid_variables 110 USE indices 86 USE arrays_3d, & 87 ONLY: ddzu, ddzw, km, tend, u, v, vsws, vswst, w 88 89 USE control_parameters, & 90 ONLY: constant_top_momentumflux, topography, use_surface_fluxes, & 91 use_top_fluxes 92 93 USE grid_variables, & 94 ONLY: ddx, ddy, ddy2, fxm, fxp, wall_v 95 96 USE indices, & 97 ONLY: nxl, nxr, nyn, nys, nysv, nzb, nzb_diff_v, nzb_v_inner, & 98 nzb_v_outer, nzt, nzt_diff 99 100 USE kinds 111 101 112 102 IMPLICIT NONE 113 103 114 INTEGER :: i, j, k 115 REAL :: kmxm, kmxp, kmzm, kmzp 116 117 REAL, DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) :: vsus 104 INTEGER(iwp) :: i !: 105 INTEGER(iwp) :: j !: 106 INTEGER(iwp) :: k !: 107 REAL(wp) :: kmxm !: 108 REAL(wp) :: kmxp !: 109 REAL(wp) :: kmzm !: 110 REAL(wp) :: kmzp !: 111 112 REAL(wp), DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) :: vsus !: 118 113 119 114 ! … … 121 116 !-- if neccessary 122 117 IF ( topography /= 'flat' ) THEN 123 CALL wall_fluxes( vsus, 0.0 , 1.0, 0.0, 0.0, nzb_v_inner, &118 CALL wall_fluxes( vsus, 0.0_wp, 1.0_wp, 0.0_wp, 0.0_wp, nzb_v_inner, & 124 119 nzb_v_outer, wall_v ) 125 120 ENDIF … … 137 132 ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) ) 138 133 139 tend(k,j,i) = tend(k,j,i) &140 & + ( kmxp * ( v(k,j,i+1) - v(k,j,i) ) * ddx &141 & + kmxp * ( u(k,j,i+1) - u(k,j-1,i+1) ) * ddy &142 & - kmxm * ( v(k,j,i) - v(k,j,i-1) ) * ddx &143 & - kmxm * ( u(k,j,i) - u(k,j-1,i) ) * ddy &144 & ) * ddx &145 & + 2.0 * ( &146 & km(k,j,i) * ( v(k,j+1,i) - v(k,j,i) ) &147 & - km(k,j-1,i) * ( v(k,j,i) - v(k,j-1,i) ) &134 tend(k,j,i) = tend(k,j,i) & 135 & + ( kmxp * ( v(k,j,i+1) - v(k,j,i) ) * ddx & 136 & + kmxp * ( u(k,j,i+1) - u(k,j-1,i+1) ) * ddy & 137 & - kmxm * ( v(k,j,i) - v(k,j,i-1) ) * ddx & 138 & - kmxm * ( u(k,j,i) - u(k,j-1,i) ) * ddy & 139 & ) * ddx & 140 & + 2.0 * ( & 141 & km(k,j,i) * ( v(k,j+1,i) - v(k,j,i) ) & 142 & - km(k,j-1,i) * ( v(k,j,i) - v(k,j-1,i) ) & 148 143 & ) * ddy2 149 144 ENDDO … … 154 149 155 150 DO k = nzb_v_inner(j,i)+1, nzb_v_outer(j,i) 156 kmxp = 0.25 * &151 kmxp = 0.25 * & 157 152 ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) ) 158 kmxm = 0.25 * &153 kmxm = 0.25 * & 159 154 ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) ) 160 155 … … 188 183 ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) ) 189 184 190 tend(k,j,i) = tend(k,j,i) &191 & + ( kmzp * ( ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1) &192 & + ( w(k,j,i) - w(k,j-1,i) ) * ddy &193 & ) &194 & - kmzm * ( ( v(k,j,i) - v(k-1,j,i) ) * ddzu(k) &195 & + ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy &196 & ) &185 tend(k,j,i) = tend(k,j,i) & 186 & + ( kmzp * ( ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1) & 187 & + ( w(k,j,i) - w(k,j-1,i) ) * ddy & 188 & ) & 189 & - kmzm * ( ( v(k,j,i) - v(k-1,j,i) ) * ddzu(k) & 190 & + ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy & 191 & ) & 197 192 & ) * ddzw(k) 198 193 ENDDO … … 204 199 !-- Difference quotient of the momentum flux is not formed over 205 200 !-- half of the grid spacing (2.0*ddzw(k)) any more, since the 206 !-- comparison with other (LES) model lshowed that the values of201 !-- comparison with other (LES) models showed that the values of 207 202 !-- the momentum flux becomes too large in this case. 208 203 !-- The term containing w(k-1,..) (see above equation) is removed here … … 212 207 ! 213 208 !-- Interpolate eddy diffusivities on staggered gridpoints 214 kmzp = 0.25 * &209 kmzp = 0.25 * & 215 210 ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) ) 216 kmzm = 0.25 * &211 kmzm = 0.25 * & 217 212 ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) ) 218 213 219 tend(k,j,i) = tend(k,j,i) &220 & + ( kmzp * ( w(k,j,i) - w(k,j-1,i) ) * ddy &221 & ) * ddzw(k) &222 & + ( kmzp * ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1) &223 & + vsws(j,i) &214 tend(k,j,i) = tend(k,j,i) & 215 & + ( kmzp * ( w(k,j,i) - w(k,j-1,i) ) * ddy & 216 & ) * ddzw(k) & 217 & + ( kmzp * ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1) & 218 & + vsws(j,i) & 224 219 & ) * ddzw(k) 225 220 ENDIF … … 232 227 ! 233 228 !-- Interpolate eddy diffusivities on staggered gridpoints 234 kmzp = 0.25 * &229 kmzp = 0.25 * & 235 230 ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) ) 236 kmzm = 0.25 * &231 kmzm = 0.25 * & 237 232 ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) ) 238 233 239 tend(k,j,i) = tend(k,j,i) &240 & - ( kmzm * ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy &241 & ) * ddzw(k) &242 & + ( -vswst(j,i) &243 & - kmzm * ( v(k,j,i) - v(k-1,j,i) ) * ddzu(k) &234 tend(k,j,i) = tend(k,j,i) & 235 & - ( kmzm * ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy & 236 & ) * ddzw(k) & 237 & + ( -vswst(j,i) & 238 & - kmzm * ( v(k,j,i) - v(k-1,j,i) ) * ddzu(k) & 244 239 & ) * ddzw(k) 245 240 ENDIF … … 256 251 SUBROUTINE diffusion_v_acc 257 252 258 USE arrays_3d 259 USE control_parameters 260 USE grid_variables 261 USE indices 253 USE arrays_3d, & 254 ONLY: ddzu, ddzw, km, tend, u, v, vsws, vswst, w 255 256 USE control_parameters, & 257 ONLY: constant_top_momentumflux, topography, use_surface_fluxes, & 258 use_top_fluxes 259 260 USE grid_variables, & 261 ONLY: ddx, ddy, ddy2, fxm, fxp, wall_v 262 263 USE indices, & 264 ONLY: i_left, i_right, j_north, j_south, nxl, nxr, nyn, nys, nzb, & 265 nzb_diff_v, nzb_v_inner, nzb_v_outer, nzt, nzt_diff 266 267 USE kinds 262 268 263 269 IMPLICIT NONE 264 270 265 INTEGER :: i, j, k 266 REAL :: kmxm, kmxp, kmzm, kmzp 267 268 REAL, DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) :: vsus 271 INTEGER(iwp) :: i !: 272 INTEGER(iwp) :: j !: 273 INTEGER(iwp) :: k !: 274 REAL(wp) :: kmxm !: 275 REAL(wp) :: kmxp !: 276 REAL(wp) :: kmzm !: 277 REAL(wp) :: kmzp !: 278 279 REAL(wp), DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) :: vsus !: 269 280 !$acc declare create ( vsus ) 270 281 … … 273 284 !-- if neccessary 274 285 IF ( topography /= 'flat' ) THEN 275 CALL wall_fluxes_acc( vsus, 0.0 , 1.0, 0.0, 0.0, nzb_v_inner,&276 nzb_v_ outer, wall_v )277 ENDIF 278 279 !$acc kernels present ( u, v, w, km, tend, vsws, vswst ) &280 !$acc present ( ddzu, ddzw, fxm, fxp, wall_v ) &286 CALL wall_fluxes_acc( vsus, 0.0_wp, 1.0_wp, 0.0_wp, 0.0_wp, & 287 nzb_v_inner, nzb_v_outer, wall_v ) 288 ENDIF 289 290 !$acc kernels present ( u, v, w, km, tend, vsws, vswst ) & 291 !$acc present ( ddzu, ddzw, fxm, fxp, wall_v ) & 281 292 !$acc present ( nzb_v_inner, nzb_v_outer, nzb_diff_v ) 282 293 DO i = i_left, i_right … … 288 299 ! 289 300 !-- Interpolate eddy diffusivities on staggered gridpoints 290 kmxp = 0.25 * &301 kmxp = 0.25 * & 291 302 ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) ) 292 kmxm = 0.25 * &303 kmxm = 0.25 * & 293 304 ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) ) 294 305 … … 309 320 !-- Wall functions at the left and right walls, respectively 310 321 DO k = 1, nzt 311 IF( k > nzb_v_inner(j,i) .AND. k <= nzb_v_outer(j,i) .AND. &322 IF( k > nzb_v_inner(j,i) .AND. k <= nzb_v_outer(j,i) .AND. & 312 323 wall_v(j,i) /= 0.0 ) THEN 313 324 314 kmxp = 0.25 * &325 kmxp = 0.25 * & 315 326 ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) ) 316 kmxm = 0.25 * &327 kmxm = 0.25 * & 317 328 ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) ) 318 329 … … 342 353 ! 343 354 !-- Interpolate eddy diffusivities on staggered gridpoints 344 kmzp = 0.25 * &355 kmzp = 0.25 * & 345 356 ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) ) 346 kmzm = 0.25 * &357 kmzm = 0.25 * & 347 358 ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) ) 348 359 … … 367 378 !-- Difference quotient of the momentum flux is not formed over 368 379 !-- half of the grid spacing (2.0*ddzw(k)) any more, since the 369 !-- comparison with other (LES) model lshowed that the values of380 !-- comparison with other (LES) models showed that the values of 370 381 !-- the momentum flux becomes too large in this case. 371 382 !-- The term containing w(k-1,..) (see above equation) is removed here … … 379 390 ! 380 391 !-- Interpolate eddy diffusivities on staggered gridpoints 381 kmzp = 0.25 * &392 kmzp = 0.25 * & 382 393 ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) ) 383 kmzm = 0.25 * &394 kmzm = 0.25 * & 384 395 ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) ) 385 396 386 tend(k,j,i) = tend(k,j,i) &387 & + ( kmzp * ( w(k,j,i) - w(k,j-1,i) ) * ddy &388 & ) * ddzw(k) &389 & + ( kmzp * ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1) &390 & + vsws(j,i) &397 tend(k,j,i) = tend(k,j,i) & 398 & + ( kmzp * ( w(k,j,i) - w(k,j-1,i) ) * ddy & 399 & ) * ddzw(k) & 400 & + ( kmzp * ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1) & 401 & + vsws(j,i) & 391 402 & ) * ddzw(k) 392 403 ENDDO … … 407 418 ! 408 419 !-- Interpolate eddy diffusivities on staggered gridpoints 409 kmzp = 0.25 * &420 kmzp = 0.25 * & 410 421 ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) ) 411 kmzm = 0.25 * &422 kmzm = 0.25 * & 412 423 ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) ) 413 424 414 tend(k,j,i) = tend(k,j,i) &415 & - ( kmzm * ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy &416 & ) * ddzw(k) &417 & + ( -vswst(j,i) &418 & - kmzm * ( v(k,j,i) - v(k-1,j,i) ) * ddzu(k) &425 tend(k,j,i) = tend(k,j,i) & 426 & - ( kmzm * ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy & 427 & ) * ddzw(k) & 428 & + ( -vswst(j,i) & 429 & - kmzm * ( v(k,j,i) - v(k-1,j,i) ) * ddzu(k) & 419 430 & ) * ddzw(k) 420 431 ENDDO … … 432 443 SUBROUTINE diffusion_v_ij( i, j ) 433 444 434 USE arrays_3d 435 USE control_parameters 436 USE grid_variables 437 USE indices 445 USE arrays_3d, & 446 ONLY: ddzu, ddzw, km, tend, u, v, vsws, vswst, w 447 448 USE control_parameters, & 449 ONLY: constant_top_momentumflux, use_surface_fluxes, use_top_fluxes 450 451 USE grid_variables, & 452 ONLY: ddx, ddy, ddy2, fxm, fxp, wall_v 453 454 USE indices, & 455 ONLY: nzb, nzb_diff_v, nzb_v_inner, nzb_v_outer, nzt, nzt_diff 456 457 USE kinds 438 458 439 459 IMPLICIT NONE 440 460 441 INTEGER :: i, j, k 442 REAL :: kmxm, kmxp, kmzm, kmzp 443 444 REAL, DIMENSION(nzb:nzt+1) :: vsus 461 INTEGER(iwp) :: i !: 462 INTEGER(iwp) :: j !: 463 INTEGER(iwp) :: k !: 464 REAL(wp) :: kmxm !: 465 REAL(wp) :: kmxp !: 466 REAL(wp) :: kmzm !: 467 REAL(wp) :: kmzp !: 468 469 REAL(wp), DIMENSION(nzb:nzt+1) :: vsus !: 445 470 446 471 ! … … 452 477 kmxm = 0.25 * ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) ) 453 478 454 tend(k,j,i) = tend(k,j,i) &455 & + ( kmxp * ( v(k,j,i+1) - v(k,j,i) ) * ddx &456 & + kmxp * ( u(k,j,i+1) - u(k,j-1,i+1) ) * ddy &457 & - kmxm * ( v(k,j,i) - v(k,j,i-1) ) * ddx &458 & - kmxm * ( u(k,j,i) - u(k,j-1,i) ) * ddy &459 & ) * ddx &460 & + 2.0 * ( &461 & km(k,j,i) * ( v(k,j+1,i) - v(k,j,i) ) &462 & - km(k,j-1,i) * ( v(k,j,i) - v(k,j-1,i) ) &479 tend(k,j,i) = tend(k,j,i) & 480 & + ( kmxp * ( v(k,j,i+1) - v(k,j,i) ) * ddx & 481 & + kmxp * ( u(k,j,i+1) - u(k,j-1,i+1) ) * ddy & 482 & - kmxm * ( v(k,j,i) - v(k,j,i-1) ) * ddx & 483 & - kmxm * ( u(k,j,i) - u(k,j-1,i) ) * ddy & 484 & ) * ddx & 485 & + 2.0 * ( & 486 & km(k,j,i) * ( v(k,j+1,i) - v(k,j,i) ) & 487 & - km(k,j-1,i) * ( v(k,j,i) - v(k,j-1,i) ) & 463 488 & ) * ddy2 464 489 ENDDO … … 470 495 ! 471 496 !-- Calculate the horizontal momentum flux v'u' 472 CALL wall_fluxes( i, j, nzb_v_inner(j,i)+1, nzb_v_outer(j,i), &473 vsus, 0.0 , 1.0, 0.0, 0.0)497 CALL wall_fluxes( i, j, nzb_v_inner(j,i)+1, nzb_v_outer(j,i), & 498 vsus, 0.0_wp, 1.0_wp, 0.0_wp, 0.0_wp ) 474 499 475 500 DO k = nzb_v_inner(j,i)+1, nzb_v_outer(j,i) 476 kmxp = 0.25 * &501 kmxp = 0.25 * & 477 502 ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) ) 478 kmxm = 0.25 * &503 kmxm = 0.25 * & 479 504 ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) ) 480 505 … … 506 531 kmzm = 0.25 * ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) ) 507 532 508 tend(k,j,i) = tend(k,j,i) &509 & + ( kmzp * ( ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1) &510 & + ( w(k,j,i) - w(k,j-1,i) ) * ddy &511 & ) &512 & - kmzm * ( ( v(k,j,i) - v(k-1,j,i) ) * ddzu(k) &513 & + ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy &514 & ) &533 tend(k,j,i) = tend(k,j,i) & 534 & + ( kmzp * ( ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1) & 535 & + ( w(k,j,i) - w(k,j-1,i) ) * ddy & 536 & ) & 537 & - kmzm * ( ( v(k,j,i) - v(k-1,j,i) ) * ddzu(k) & 538 & + ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy & 539 & ) & 515 540 & ) * ddzw(k) 516 541 ENDDO … … 522 547 !-- Difference quotient of the momentum flux is not formed over half of 523 548 !-- the grid spacing (2.0*ddzw(k)) any more, since the comparison with 524 !-- other (LES) model lshowed that the values of the momentum flux becomes549 !-- other (LES) models showed that the values of the momentum flux becomes 525 550 !-- too large in this case. 526 551 !-- The term containing w(k-1,..) (see above equation) is removed here … … 533 558 kmzm = 0.25 * ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) ) 534 559 535 tend(k,j,i) = tend(k,j,i) &536 & + ( kmzp * ( w(k,j,i) - w(k,j-1,i) ) * ddy &537 & ) * ddzw(k) &538 & + ( kmzp * ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1) &539 & + vsws(j,i) &560 tend(k,j,i) = tend(k,j,i) & 561 & + ( kmzp * ( w(k,j,i) - w(k,j-1,i) ) * ddy & 562 & ) * ddzw(k) & 563 & + ( kmzp * ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1) & 564 & + vsws(j,i) & 540 565 & ) * ddzw(k) 541 566 ENDIF … … 553 578 ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) ) 554 579 555 tend(k,j,i) = tend(k,j,i) &556 & - ( kmzm * ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy &557 & ) * ddzw(k) &558 & + ( -vswst(j,i) &559 & - kmzm * ( v(k,j,i) - v(k-1,j,i) ) * ddzu(k) &580 tend(k,j,i) = tend(k,j,i) & 581 & - ( kmzm * ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy & 582 & ) * ddzw(k) & 583 & + ( -vswst(j,i) & 584 & - kmzm * ( v(k,j,i) - v(k-1,j,i) ) * ddzu(k) & 560 585 & ) * ddzw(k) 561 586 ENDIF
Note: See TracChangeset
for help on using the changeset viewer.