Changeset 978 for palm/trunk
- Timestamp:
- Aug 9, 2012 8:28:32 AM (12 years ago)
- Location:
- palm/trunk
- Files:
-
- 27 edited
-
. (modified) (1 prop)
-
EXAMPLES/dvr (modified) (1 prop)
-
SOURCE (modified) (1 prop)
-
SOURCE/average_3d_data.f90 (modified) (2 diffs)
-
SOURCE/boundary_conds.f90 (modified) (19 diffs)
-
SOURCE/check_parameters.f90 (modified) (10 diffs)
-
SOURCE/data_output_2d.f90 (modified) (2 diffs)
-
SOURCE/diffusion_u.f90 (modified) (13 diffs)
-
SOURCE/diffusion_v.f90 (modified) (11 diffs)
-
SOURCE/diffusion_w.f90 (modified) (9 diffs)
-
SOURCE/header.f90 (modified) (5 diffs)
-
SOURCE/init_1d_model.f90 (modified) (8 diffs)
-
SOURCE/init_3d_model.f90 (modified) (7 diffs)
-
SOURCE/init_grid.f90 (modified) (9 diffs)
-
SOURCE/init_masks.f90 (modified) (2 diffs)
-
SOURCE/init_pegrid.f90 (modified) (7 diffs)
-
SOURCE/modules.f90 (modified) (10 diffs)
-
SOURCE/parin.f90 (modified) (4 diffs)
-
SOURCE/poismg.f90 (modified) (7 diffs)
-
SOURCE/prandtl_fluxes.f90 (modified) (5 diffs)
-
SOURCE/prognostic_equations.f90 (modified) (14 diffs)
-
SOURCE/read_3d_binary.f90 (modified) (2 diffs)
-
SOURCE/read_var_list.f90 (modified) (5 diffs)
-
SOURCE/sum_up_3d_data.f90 (modified) (3 diffs)
-
SOURCE/timestep.f90 (modified) (2 diffs)
-
SOURCE/write_3d_binary.f90 (modified) (2 diffs)
-
SOURCE/write_var_list.f90 (modified) (5 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk
- Property svn:mergeinfo changed
/palm/branches/fricke (added) merged: 942-944,967-968,971-972,977
- Property svn:mergeinfo changed
-
palm/trunk/EXAMPLES/dvr
- Property svn:mergeinfo changed
/palm/branches/fricke/EXAMPLES/dvr (added) merged: 972,977
- Property svn:mergeinfo changed
-
palm/trunk/SOURCE
- Property svn:mergeinfo changed
/palm/branches/fricke/SOURCE (added) merged: 967-968,971-972,977
- Property svn:mergeinfo changed
-
palm/trunk/SOURCE/average_3d_data.f90
r772 r978 4 4 ! Current revisions: 5 5 ! ----------------- 6 ! +z0h_av 6 7 ! 7 8 ! Former revisions: … … 284 285 ENDDO 285 286 287 CASE ( 'z0h*' ) 288 DO i = nxlg, nxrg 289 DO j = nysg, nyng 290 z0h_av(j,i) = z0h_av(j,i) / REAL( average_count_3d ) 291 ENDDO 292 ENDDO 293 286 294 CASE DEFAULT 287 295 ! -
palm/trunk/SOURCE/boundary_conds.f90
r876 r978 4 4 ! Current revisions: 5 5 ! ----------------- 6 ! 6 ! Neumann boudnary conditions are added at the inflow boundary for the SGS-TKE. 7 ! Outflow boundary conditions for the velocity components can be set to Neumann 8 ! conditions or to radiation conditions with a horizontal averaged phase 9 ! velocity. 7 10 ! 8 11 ! … … 190 193 q_p(nzt+1,:,:) = q_p(nzt,:,:) + bc_q_t_val * dzu(nzt+1) 191 194 192 193 195 ENDIF 194 196 … … 198 200 !-- Since in prognostic_equations (cache optimized version) these levels are 199 201 !-- handled as a prognostic level, boundary values have to be restored here. 202 !-- For the SGS-TKE, Neumann boundary conditions are used at the inflow. 200 203 IF ( inflow_s ) THEN 201 204 v_p(:,nys,:) = v_p(:,nys-1,:) 205 IF ( .NOT. constant_diffusion ) e_p(:,nys-1,:) = e_p(:,nys,:) 206 ELSEIF ( inflow_n ) THEN 207 IF ( .NOT. constant_diffusion ) e_p(:,nyn+1,:) = e_p(:,nyn,:) 202 208 ELSEIF ( inflow_l ) THEN 203 209 u_p(:,:,nxl) = u_p(:,:,nxl-1) 210 IF ( .NOT. constant_diffusion ) e_p(:,:,nxl-1) = e_p(:,:,nxl) 211 ELSEIF ( inflow_r ) THEN 212 IF ( .NOT. constant_diffusion ) e_p(:,:,nxr+1) = e_p(:,:,nxr) 204 213 ENDIF 205 214 … … 227 236 228 237 ! 229 !-- Radiation boundary condition for the velocities at the respective outflow 238 !-- Neumann or Radiation boundary condition for the velocities at the 239 !-- respective outflow 230 240 IF ( outflow_s ) THEN 231 241 232 c_max = dy / dt_3d 233 234 DO i = nxlg, nxrg 242 IF ( bc_ns_dirneu ) THEN 243 u(:,-1,:) = u(:,0,:) 244 v(:,0,:) = v(:,1,:) 245 w(:,-1,:) = w(:,0,:) 246 ELSEIF ( bc_ns_dirrad ) THEN 247 248 c_max = dy / dt_3d 249 250 c_u_m_l = 0.0 251 c_v_m_l = 0.0 252 c_w_m_l = 0.0 253 254 c_u_m = 0.0 255 c_v_m = 0.0 256 c_w_m = 0.0 257 258 ! 259 !-- Calculate the phase speeds for u,v, and w, first local and then 260 !-- average parallel along the outflow boundary. 235 261 DO k = nzb+1, nzt+1 236 237 ! 238 !-- Calculate the phase speeds for u,v, and w. In case of using a 239 !-- Runge-Kutta scheme, do this for the first substep only and then 240 !-- reuse this values for the further substeps. 241 IF ( intermediate_timestep_count == 1 ) THEN 262 DO i = nxl, nxr 242 263 243 264 denom = u_m_s(k,0,i) - u_m_s(k,1,i) 244 265 245 266 IF ( denom /= 0.0 ) THEN 246 c_u(k,i) = -c_max * ( u(k,0,i) - u_m_s(k,0,i) ) / denom 267 c_u(k,i) = -c_max * ( u(k,0,i) - u_m_s(k,0,i) ) & 268 / ( denom * tsc(2) ) 247 269 IF ( c_u(k,i) < 0.0 ) THEN 248 270 c_u(k,i) = 0.0 … … 257 279 258 280 IF ( denom /= 0.0 ) THEN 259 c_v(k,i) = -c_max * ( v(k,1,i) - v_m_s(k,1,i) ) / denom 281 c_v(k,i) = -c_max * ( v(k,1,i) - v_m_s(k,1,i) ) & 282 / ( denom * tsc(2) ) 260 283 IF ( c_v(k,i) < 0.0 ) THEN 261 284 c_v(k,i) = 0.0 … … 270 293 271 294 IF ( denom /= 0.0 ) THEN 272 c_w(k,i) = -c_max * ( w(k,0,i) - w_m_s(k,0,i) ) / denom 295 c_w(k,i) = -c_max * ( w(k,0,i) - w_m_s(k,0,i) ) & 296 / ( denom * tsc(2) ) 273 297 IF ( c_w(k,i) < 0.0 ) THEN 274 298 c_w(k,i) = 0.0 … … 280 304 ENDIF 281 305 282 ! 283 !-- Save old timelevels for the next timestep 284 u_m_s(k,:,i) = u(k,0:1,i) 285 v_m_s(k,:,i) = v(k,1:2,i) 286 w_m_s(k,:,i) = w(k,0:1,i) 287 288 ENDIF 289 290 ! 291 !-- Calculate the new velocities 292 u_p(k,-1,i) = u(k,-1,i) - dt_3d * tsc(2) * c_u(k,i) * & 306 c_u_m_l(k) = c_u_m_l(k) + c_u(k,i) 307 c_v_m_l(k) = c_v_m_l(k) + c_v(k,i) 308 c_w_m_l(k) = c_w_m_l(k) + c_w(k,i) 309 310 ENDDO 311 ENDDO 312 313 #if defined( __parallel ) 314 IF ( collective_wait ) CALL MPI_BARRIER( comm1dx, ierr ) 315 CALL MPI_ALLREDUCE( c_u_m_l(nzb+1), c_u_m(nzb+1), nzt-nzb, MPI_REAL, & 316 MPI_SUM, comm1dx, ierr ) 317 IF ( collective_wait ) CALL MPI_BARRIER( comm1dx, ierr ) 318 CALL MPI_ALLREDUCE( c_v_m_l(nzb+1), c_v_m(nzb+1), nzt-nzb, MPI_REAL, & 319 MPI_SUM, comm1dx, ierr ) 320 IF ( collective_wait ) CALL MPI_BARRIER( comm1dx, ierr ) 321 CALL MPI_ALLREDUCE( c_w_m_l(nzb+1), c_w_m(nzb+1), nzt-nzb, MPI_REAL, & 322 MPI_SUM, comm1dx, ierr ) 323 #else 324 c_u_m = c_u_m_l 325 c_v_m = c_v_m_l 326 c_w_m = c_w_m_l 327 #endif 328 329 c_u_m = c_u_m / (nx+1) 330 c_v_m = c_v_m / (nx+1) 331 c_w_m = c_w_m / (nx+1) 332 333 ! 334 !-- Save old timelevels for the next timestep 335 IF ( intermediate_timestep_count == 1 ) THEN 336 u_m_s(:,:,:) = u(:,0:1,:) 337 v_m_s(:,:,:) = v(:,1:2,:) 338 w_m_s(:,:,:) = w(:,0:1,:) 339 ENDIF 340 341 ! 342 !-- Calculate the new velocities 343 DO k = nzb+1, nzt+1 344 DO i = nxlg, nxrg 345 u_p(k,-1,i) = u(k,-1,i) - dt_3d * tsc(2) * c_u_m(k) * & 293 346 ( u(k,-1,i) - u(k,0,i) ) * ddy 294 347 295 v_p(k,0,i) = v(k,0,i) - dt_3d * tsc(2) * c_v(k,i) *&348 v_p(k,0,i) = v(k,0,i) - dt_3d * tsc(2) * c_v_m(k) * & 296 349 ( v(k,0,i) - v(k,1,i) ) * ddy 297 350 298 w_p(k,-1,i) = w(k,-1,i) - dt_3d * tsc(2) * c_w(k,i) *&351 w_p(k,-1,i) = w(k,-1,i) - dt_3d * tsc(2) * c_w_m(k) * & 299 352 ( w(k,-1,i) - w(k,0,i) ) * ddy 300 301 ENDDO 302 ENDDO 303 304 ! 305 !-- Bottom boundary at the outflow 306 IF ( ibc_uv_b == 0 ) THEN 307 u_p(nzb,-1,:) = 0.0 308 v_p(nzb,0,:) = 0.0 309 ELSE 310 u_p(nzb,-1,:) = u_p(nzb+1,-1,:) 311 v_p(nzb,0,:) = v_p(nzb+1,0,:) 312 ENDIF 313 w_p(nzb,-1,:) = 0.0 314 315 ! 316 !-- Top boundary at the outflow 317 IF ( ibc_uv_t == 0 ) THEN 318 u_p(nzt+1,-1,:) = u_init(nzt+1) 319 v_p(nzt+1,0,:) = v_init(nzt+1) 320 ELSE 321 u_p(nzt+1,-1,:) = u(nzt,-1,:) 322 v_p(nzt+1,0,:) = v(nzt,0,:) 323 ENDIF 324 w_p(nzt:nzt+1,-1,:) = 0.0 353 ENDDO 354 ENDDO 355 356 ! 357 !-- Bottom boundary at the outflow 358 IF ( ibc_uv_b == 0 ) THEN 359 u_p(nzb,-1,:) = 0.0 360 v_p(nzb,0,:) = 0.0 361 ELSE 362 u_p(nzb,-1,:) = u_p(nzb+1,-1,:) 363 v_p(nzb,0,:) = v_p(nzb+1,0,:) 364 ENDIF 365 w_p(nzb,-1,:) = 0.0 366 367 ! 368 !-- Top boundary at the outflow 369 IF ( ibc_uv_t == 0 ) THEN 370 u_p(nzt+1,-1,:) = u_init(nzt+1) 371 v_p(nzt+1,0,:) = v_init(nzt+1) 372 ELSE 373 u_p(nzt+1,-1,:) = u(nzt,-1,:) 374 v_p(nzt+1,0,:) = v(nzt,0,:) 375 ENDIF 376 w_p(nzt:nzt+1,-1,:) = 0.0 377 378 ENDIF 325 379 326 380 ENDIF … … 328 382 IF ( outflow_n ) THEN 329 383 330 c_max = dy / dt_3d 331 332 DO i = nxlg, nxrg 384 IF ( bc_ns_neudir ) THEN 385 u(:,ny+1,:) = u(:,ny,:) 386 v(:,ny+1,:) = v(:,ny,:) 387 w(:,ny+1,:) = w(:,ny,:) 388 ELSEIF ( bc_ns_dirrad ) THEN 389 390 c_max = dy / dt_3d 391 392 c_u_m_l = 0.0 393 c_v_m_l = 0.0 394 c_w_m_l = 0.0 395 396 c_u_m = 0.0 397 c_v_m = 0.0 398 c_w_m = 0.0 399 400 ! 401 !-- Calculate the phase speeds for u,v, and w, first local and then 402 !-- average parallel along the outflow boundary. 333 403 DO k = nzb+1, nzt+1 334 335 ! 336 !-- Calculate the phase speeds for u,v, and w. In case of using a 337 !-- Runge-Kutta scheme, do this for the first substep only and then 338 !-- reuse this values for the further substeps. 339 IF ( intermediate_timestep_count == 1 ) THEN 404 DO i = nxl, nxr 340 405 341 406 denom = u_m_n(k,ny,i) - u_m_n(k,ny-1,i) 342 407 343 408 IF ( denom /= 0.0 ) THEN 344 c_u(k,i) = -c_max * ( u(k,ny,i) - u_m_n(k,ny,i) ) / denom 409 c_u(k,i) = -c_max * ( u(k,ny,i) - u_m_n(k,ny,i) ) & 410 / ( denom * tsc(2) ) 345 411 IF ( c_u(k,i) < 0.0 ) THEN 346 412 c_u(k,i) = 0.0 … … 355 421 356 422 IF ( denom /= 0.0 ) THEN 357 c_v(k,i) = -c_max * ( v(k,ny,i) - v_m_n(k,ny,i) ) / denom 423 c_v(k,i) = -c_max * ( v(k,ny,i) - v_m_n(k,ny,i) ) & 424 / ( denom * tsc(2) ) 358 425 IF ( c_v(k,i) < 0.0 ) THEN 359 426 c_v(k,i) = 0.0 … … 368 435 369 436 IF ( denom /= 0.0 ) THEN 370 c_w(k,i) = -c_max * ( w(k,ny,i) - w_m_n(k,ny,i) ) / denom 437 c_w(k,i) = -c_max * ( w(k,ny,i) - w_m_n(k,ny,i) ) & 438 / ( denom * tsc(2) ) 371 439 IF ( c_w(k,i) < 0.0 ) THEN 372 440 c_w(k,i) = 0.0 … … 378 446 ENDIF 379 447 380 ! 381 !-- Swap timelevels for the next timestep 382 u_m_n(k,:,i) = u(k,ny-1:ny,i) 383 v_m_n(k,:,i) = v(k,ny-1:ny,i) 384 w_m_n(k,:,i) = w(k,ny-1:ny,i) 385 386 ENDIF 387 388 ! 389 !-- Calculate the new velocities 390 u_p(k,ny+1,i) = u(k,ny+1,i) - dt_3d * tsc(2) * c_u(k,i) * & 391 ( u(k,ny+1,i) - u(k,ny,i) ) * ddy 392 393 v_p(k,ny+1,i) = v(k,ny+1,i) - dt_3d * tsc(2) * c_v(k,i) * & 394 ( v(k,ny+1,i) - v(k,ny,i) ) * ddy 395 396 w_p(k,ny+1,i) = w(k,ny+1,i) - dt_3d * tsc(2) * c_w(k,i) * & 397 ( w(k,ny+1,i) - w(k,ny,i) ) * ddy 398 399 ENDDO 400 ENDDO 401 402 ! 403 !-- Bottom boundary at the outflow 404 IF ( ibc_uv_b == 0 ) THEN 405 u_p(nzb,ny+1,:) = 0.0 406 v_p(nzb,ny+1,:) = 0.0 407 ELSE 408 u_p(nzb,ny+1,:) = u_p(nzb+1,ny+1,:) 409 v_p(nzb,ny+1,:) = v_p(nzb+1,ny+1,:) 410 ENDIF 411 w_p(nzb,ny+1,:) = 0.0 412 413 ! 414 !-- Top boundary at the outflow 415 IF ( ibc_uv_t == 0 ) THEN 416 u_p(nzt+1,ny+1,:) = u_init(nzt+1) 417 v_p(nzt+1,ny+1,:) = v_init(nzt+1) 418 ELSE 419 u_p(nzt+1,ny+1,:) = u_p(nzt,nyn+1,:) 420 v_p(nzt+1,ny+1,:) = v_p(nzt,nyn+1,:) 421 ENDIF 422 w_p(nzt:nzt+1,ny+1,:) = 0.0 448 c_u_m_l(k) = c_u_m_l(k) + c_u(k,i) 449 c_v_m_l(k) = c_v_m_l(k) + c_v(k,i) 450 c_w_m_l(k) = c_w_m_l(k) + c_w(k,i) 451 452 ENDDO 453 ENDDO 454 455 #if defined( __parallel ) 456 IF ( collective_wait ) CALL MPI_BARRIER( comm1dx, ierr ) 457 CALL MPI_ALLREDUCE( c_u_m_l(nzb+1), c_u_m(nzb+1), nzt-nzb, MPI_REAL, & 458 MPI_SUM, comm1dx, ierr ) 459 IF ( collective_wait ) CALL MPI_BARRIER( comm1dx, ierr ) 460 CALL MPI_ALLREDUCE( c_v_m_l(nzb+1), c_v_m(nzb+1), nzt-nzb, MPI_REAL, & 461 MPI_SUM, comm1dx, ierr ) 462 IF ( collective_wait ) CALL MPI_BARRIER( comm1dx, ierr ) 463 CALL MPI_ALLREDUCE( c_w_m_l(nzb+1), c_w_m(nzb+1), nzt-nzb, MPI_REAL, & 464 MPI_SUM, comm1dx, ierr ) 465 #else 466 c_u_m = c_u_m_l 467 c_v_m = c_v_m_l 468 c_w_m = c_w_m_l 469 #endif 470 471 c_u_m = c_u_m / (nx+1) 472 c_v_m = c_v_m / (nx+1) 473 c_w_m = c_w_m / (nx+1) 474 475 ! 476 !-- Save old timelevels for the next timestep 477 IF ( intermediate_timestep_count == 1 ) THEN 478 u_m_n(:,:,:) = u(:,ny-1:ny,:) 479 v_m_n(:,:,:) = v(:,ny-1:ny,:) 480 w_m_n(:,:,:) = w(:,ny-1:ny,:) 481 ENDIF 482 483 ! 484 !-- Calculate the new velocities 485 DO k = nzb+1, nzt+1 486 DO i = nxlg, nxrg 487 u_p(k,ny+1,i) = u(k,ny+1,i) - dt_3d * tsc(2) * c_u_m(k) * & 488 ( u(k,ny+1,i) - u(k,ny,i) ) * ddy 489 490 v_p(k,ny+1,i) = v(k,ny+1,i) - dt_3d * tsc(2) * c_v_m(k) * & 491 ( v(k,ny+1,i) - v(k,ny,i) ) * ddy 492 493 w_p(k,ny+1,i) = w(k,ny+1,i) - dt_3d * tsc(2) * c_w_m(k) * & 494 ( w(k,ny+1,i) - w(k,ny,i) ) * ddy 495 ENDDO 496 ENDDO 497 498 ! 499 !-- Bottom boundary at the outflow 500 IF ( ibc_uv_b == 0 ) THEN 501 u_p(nzb,ny+1,:) = 0.0 502 v_p(nzb,ny+1,:) = 0.0 503 ELSE 504 u_p(nzb,ny+1,:) = u_p(nzb+1,ny+1,:) 505 v_p(nzb,ny+1,:) = v_p(nzb+1,ny+1,:) 506 ENDIF 507 w_p(nzb,ny+1,:) = 0.0 508 509 ! 510 !-- Top boundary at the outflow 511 IF ( ibc_uv_t == 0 ) THEN 512 u_p(nzt+1,ny+1,:) = u_init(nzt+1) 513 v_p(nzt+1,ny+1,:) = v_init(nzt+1) 514 ELSE 515 u_p(nzt+1,ny+1,:) = u_p(nzt,nyn+1,:) 516 v_p(nzt+1,ny+1,:) = v_p(nzt,nyn+1,:) 517 ENDIF 518 w_p(nzt:nzt+1,ny+1,:) = 0.0 519 520 ENDIF 423 521 424 522 ENDIF … … 426 524 IF ( outflow_l ) THEN 427 525 428 c_max = dx / dt_3d 429 430 DO j = nysg, nyng 526 IF ( bc_lr_neudir ) THEN 527 u(:,:,-1) = u(:,:,0) 528 v(:,:,0) = v(:,:,1) 529 w(:,:,-1) = w(:,:,0) 530 ELSEIF ( bc_ns_dirrad ) THEN 531 532 c_max = dx / dt_3d 533 534 c_u_m_l = 0.0 535 c_v_m_l = 0.0 536 c_w_m_l = 0.0 537 538 c_u_m = 0.0 539 c_v_m = 0.0 540 c_w_m = 0.0 541 542 ! 543 !-- Calculate the phase speeds for u,v, and w, first local and then 544 !-- average parallel along the outflow boundary. 431 545 DO k = nzb+1, nzt+1 432 433 ! 434 !-- Calculate the phase speeds for u,v, and w. In case of using a 435 !-- Runge-Kutta scheme, do this for the first substep only and then 436 !-- reuse this values for the further substeps. 437 IF ( intermediate_timestep_count == 1 ) THEN 546 DO j = nys, nyn 438 547 439 548 denom = u_m_l(k,j,1) - u_m_l(k,j,2) 440 549 441 550 IF ( denom /= 0.0 ) THEN 442 c_u(k,j) = -c_max * ( u(k,j,1) - u_m_l(k,j,1) ) / denom 551 c_u(k,j) = -c_max * ( u(k,j,1) - u_m_l(k,j,1) ) & 552 / ( denom * tsc(2) ) 443 553 IF ( c_u(k,j) < 0.0 ) THEN 444 554 c_u(k,j) = 0.0 … … 453 563 454 564 IF ( denom /= 0.0 ) THEN 455 c_v(k,j) = -c_max * ( v(k,j,0) - v_m_l(k,j,0) ) / denom 565 c_v(k,j) = -c_max * ( v(k,j,0) - v_m_l(k,j,0) ) & 566 / ( denom * tsc(2) ) 456 567 IF ( c_v(k,j) < 0.0 ) THEN 457 568 c_v(k,j) = 0.0 … … 466 577 467 578 IF ( denom /= 0.0 ) THEN 468 c_w(k,j) = -c_max * ( w(k,j,0) - w_m_l(k,j,0) ) / denom 579 c_w(k,j) = -c_max * ( w(k,j,0) - w_m_l(k,j,0) ) & 580 / ( denom * tsc(2) ) 469 581 IF ( c_w(k,j) < 0.0 ) THEN 470 582 c_w(k,j) = 0.0 … … 476 588 ENDIF 477 589 478 ! 479 !-- Swap timelevels for the next timestep 480 u_m_l(k,j,:) = u(k,j,1:2) 481 v_m_l(k,j,:) = v(k,j,0:1) 482 w_m_l(k,j,:) = w(k,j,0:1) 483 484 ENDIF 485 486 ! 487 !-- Calculate the new velocities 488 u_p(k,j,0) = u(k,j,0) - dt_3d * tsc(2) * c_u(k,j) * & 590 c_u_m_l(k) = c_u_m_l(k) + c_u(k,j) 591 c_v_m_l(k) = c_v_m_l(k) + c_v(k,j) 592 c_w_m_l(k) = c_w_m_l(k) + c_w(k,j) 593 594 ENDDO 595 ENDDO 596 597 #if defined( __parallel ) 598 IF ( collective_wait ) CALL MPI_BARRIER( comm1dy, ierr ) 599 CALL MPI_ALLREDUCE( c_u_m_l(nzb+1), c_u_m(nzb+1), nzt-nzb, MPI_REAL, & 600 MPI_SUM, comm1dy, ierr ) 601 IF ( collective_wait ) CALL MPI_BARRIER( comm1dy, ierr ) 602 CALL MPI_ALLREDUCE( c_v_m_l(nzb+1), c_v_m(nzb+1), nzt-nzb, MPI_REAL, & 603 MPI_SUM, comm1dy, ierr ) 604 IF ( collective_wait ) CALL MPI_BARRIER( comm1dy, ierr ) 605 CALL MPI_ALLREDUCE( c_w_m_l(nzb+1), c_w_m(nzb+1), nzt-nzb, MPI_REAL, & 606 MPI_SUM, comm1dy, ierr ) 607 #else 608 c_u_m = c_u_m_l 609 c_v_m = c_v_m_l 610 c_w_m = c_w_m_l 611 #endif 612 613 c_u_m = c_u_m / (ny+1) 614 c_v_m = c_v_m / (ny+1) 615 c_w_m = c_w_m / (ny+1) 616 617 ! 618 !-- Save old timelevels for the next timestep 619 IF ( intermediate_timestep_count == 1 ) THEN 620 u_m_l(:,:,:) = u(:,:,1:2) 621 v_m_l(:,:,:) = v(:,:,0:1) 622 w_m_l(:,:,:) = w(:,:,0:1) 623 ENDIF 624 625 ! 626 !-- Calculate the new velocities 627 DO k = nzb+1, nzt+1 628 DO i = nxlg, nxrg 629 u_p(k,j,0) = u(k,j,0) - dt_3d * tsc(2) * c_u_m(k) * & 489 630 ( u(k,j,0) - u(k,j,1) ) * ddx 490 631 491 v_p(k,j,-1) = v(k,j,-1) - dt_3d * tsc(2) * c_v(k,j) *&632 v_p(k,j,-1) = v(k,j,-1) - dt_3d * tsc(2) * c_v_m(k) * & 492 633 ( v(k,j,-1) - v(k,j,0) ) * ddx 493 634 494 w_p(k,j,-1) = w(k,j,-1) - dt_3d * tsc(2) * c_w(k,j) *&635 w_p(k,j,-1) = w(k,j,-1) - dt_3d * tsc(2) * c_w_m(k) * & 495 636 ( w(k,j,-1) - w(k,j,0) ) * ddx 496 497 ENDDO 498 ENDDO 499 500 ! 501 !-- Bottom boundary at the outflow 502 IF ( ibc_uv_b == 0 ) THEN 503 u_p(nzb,:,0) = 0.0 504 v_p(nzb,:,-1) = 0.0 505 ELSE 506 u_p(nzb,:,0) = u_p(nzb+1,:,0) 507 v_p(nzb,:,-1) = v_p(nzb+1,:,-1) 508 ENDIF 509 w_p(nzb,:,-1) = 0.0 510 511 ! 512 !-- Top boundary at the outflow 513 IF ( ibc_uv_t == 0 ) THEN 514 u_p(nzt+1,:,-1) = u_init(nzt+1) 515 v_p(nzt+1,:,-1) = v_init(nzt+1) 516 ELSE 517 u_p(nzt+1,:,-1) = u_p(nzt,:,-1) 518 v_p(nzt+1,:,-1) = v_p(nzt,:,-1) 519 ENDIF 520 w_p(nzt:nzt+1,:,-1) = 0.0 637 ENDDO 638 ENDDO 639 640 ! 641 !-- Bottom boundary at the outflow 642 IF ( ibc_uv_b == 0 ) THEN 643 u_p(nzb,:,0) = 0.0 644 v_p(nzb,:,-1) = 0.0 645 ELSE 646 u_p(nzb,:,0) = u_p(nzb+1,:,0) 647 v_p(nzb,:,-1) = v_p(nzb+1,:,-1) 648 ENDIF 649 w_p(nzb,:,-1) = 0.0 650 651 ! 652 !-- Top boundary at the outflow 653 IF ( ibc_uv_t == 0 ) THEN 654 u_p(nzt+1,:,-1) = u_init(nzt+1) 655 v_p(nzt+1,:,-1) = v_init(nzt+1) 656 ELSE 657 u_p(nzt+1,:,-1) = u_p(nzt,:,-1) 658 v_p(nzt+1,:,-1) = v_p(nzt,:,-1) 659 ENDIF 660 w_p(nzt:nzt+1,:,-1) = 0.0 661 662 ENDIF 521 663 522 664 ENDIF … … 524 666 IF ( outflow_r ) THEN 525 667 526 c_max = dx / dt_3d 527 528 DO j = nysg, nyng 668 IF ( bc_lr_dirneu ) THEN 669 u(:,:,nx+1) = u(:,:,nx) 670 v(:,:,nx+1) = v(:,:,nx) 671 w(:,:,nx+1) = w(:,:,nx) 672 ELSEIF ( bc_ns_dirrad ) THEN 673 674 c_max = dx / dt_3d 675 676 c_u_m_l = 0.0 677 c_v_m_l = 0.0 678 c_w_m_l = 0.0 679 680 c_u_m = 0.0 681 c_v_m = 0.0 682 c_w_m = 0.0 683 684 ! 685 !-- Calculate the phase speeds for u,v, and w, first local and then 686 !-- average parallel along the outflow boundary. 529 687 DO k = nzb+1, nzt+1 530 531 ! 532 !-- Calculate the phase speeds for u,v, and w. In case of using a 533 !-- Runge-Kutta scheme, do this for the first substep only and then 534 !-- reuse this values for the further substeps. 535 IF ( intermediate_timestep_count == 1 ) THEN 688 DO j = nys, nyn 536 689 537 690 denom = u_m_r(k,j,nx) - u_m_r(k,j,nx-1) 538 691 539 692 IF ( denom /= 0.0 ) THEN 540 c_u(k,j) = -c_max * ( u(k,j,nx) - u_m_r(k,j,nx) ) / denom 693 c_u(k,j) = -c_max * ( u(k,j,nx) - u_m_r(k,j,nx) ) & 694 / ( denom * tsc(2) ) 541 695 IF ( c_u(k,j) < 0.0 ) THEN 542 696 c_u(k,j) = 0.0 … … 551 705 552 706 IF ( denom /= 0.0 ) THEN 553 c_v(k,j) = -c_max * ( v(k,j,nx) - v_m_r(k,j,nx) ) / denom 707 c_v(k,j) = -c_max * ( v(k,j,nx) - v_m_r(k,j,nx) ) & 708 / ( denom * tsc(2) ) 554 709 IF ( c_v(k,j) < 0.0 ) THEN 555 710 c_v(k,j) = 0.0 … … 564 719 565 720 IF ( denom /= 0.0 ) THEN 566 c_w(k,j) = -c_max * ( w(k,j,nx) - w_m_r(k,j,nx) ) / denom 721 c_w(k,j) = -c_max * ( w(k,j,nx) - w_m_r(k,j,nx) ) & 722 / ( denom * tsc(2) ) 567 723 IF ( c_w(k,j) < 0.0 ) THEN 568 724 c_w(k,j) = 0.0 … … 574 730 ENDIF 575 731 576 ! 577 !-- Swap timelevels for the next timestep 578 u_m_r(k,j,:) = u(k,j,nx-1:nx) 579 v_m_r(k,j,:) = v(k,j,nx-1:nx) 580 w_m_r(k,j,:) = w(k,j,nx-1:nx) 581 582 ENDIF 583 584 ! 585 !-- Calculate the new velocities 586 u_p(k,j,nx+1) = u(k,j,nx+1) - dt_3d * tsc(2) * c_u(k,j) * & 587 ( u(k,j,nx+1) - u(k,j,nx) ) * ddx 588 589 v_p(k,j,nx+1) = v(k,j,nx+1) - dt_3d * tsc(2) * c_v(k,j) * & 590 ( v(k,j,nx+1) - v(k,j,nx) ) * ddx 591 592 w_p(k,j,nx+1) = w(k,j,nx+1) - dt_3d * tsc(2) * c_w(k,j) * & 593 ( w(k,j,nx+1) - w(k,j,nx) ) * ddx 594 595 ENDDO 596 ENDDO 597 598 ! 599 !-- Bottom boundary at the outflow 600 IF ( ibc_uv_b == 0 ) THEN 601 u_p(nzb,:,nx+1) = 0.0 602 v_p(nzb,:,nx+1) = 0.0 603 ELSE 604 u_p(nzb,:,nx+1) = u_p(nzb+1,:,nx+1) 605 v_p(nzb,:,nx+1) = v_p(nzb+1,:,nx+1) 606 ENDIF 607 w_p(nzb,:,nx+1) = 0.0 608 609 ! 610 !-- Top boundary at the outflow 611 IF ( ibc_uv_t == 0 ) THEN 612 u_p(nzt+1,:,nx+1) = u_init(nzt+1) 613 v_p(nzt+1,:,nx+1) = v_init(nzt+1) 614 ELSE 615 u_p(nzt+1,:,nx+1) = u_p(nzt,:,nx+1) 616 v_p(nzt+1,:,nx+1) = v_p(nzt,:,nx+1) 617 ENDIF 618 w(nzt:nzt+1,:,nx+1) = 0.0 732 c_u_m_l(k) = c_u_m_l(k) + c_u(k,j) 733 c_v_m_l(k) = c_v_m_l(k) + c_v(k,j) 734 c_w_m_l(k) = c_w_m_l(k) + c_w(k,j) 735 736 ENDDO 737 ENDDO 738 739 #if defined( __parallel ) 740 IF ( collective_wait ) CALL MPI_BARRIER( comm1dy, ierr ) 741 CALL MPI_ALLREDUCE( c_u_m_l(nzb+1), c_u_m(nzb+1), nzt-nzb, MPI_REAL, & 742 MPI_SUM, comm1dy, ierr ) 743 IF ( collective_wait ) CALL MPI_BARRIER( comm1dy, ierr ) 744 CALL MPI_ALLREDUCE( c_v_m_l(nzb+1), c_v_m(nzb+1), nzt-nzb, MPI_REAL, & 745 MPI_SUM, comm1dy, ierr ) 746 IF ( collective_wait ) CALL MPI_BARRIER( comm1dy, ierr ) 747 CALL MPI_ALLREDUCE( c_w_m_l(nzb+1), c_w_m(nzb+1), nzt-nzb, MPI_REAL, & 748 MPI_SUM, comm1dy, ierr ) 749 #else 750 c_u_m = c_u_m_l 751 c_v_m = c_v_m_l 752 c_w_m = c_w_m_l 753 #endif 754 755 c_u_m = c_u_m / (ny+1) 756 c_v_m = c_v_m / (ny+1) 757 c_w_m = c_w_m / (ny+1) 758 759 ! 760 !-- Save old timelevels for the next timestep 761 IF ( intermediate_timestep_count == 1 ) THEN 762 u_m_r(:,:,:) = u(:,:,nx-1:nx) 763 v_m_r(:,:,:) = v(:,:,nx-1:nx) 764 w_m_r(:,:,:) = w(:,:,nx-1:nx) 765 ENDIF 766 767 ! 768 !-- Calculate the new velocities 769 DO k = nzb+1, nzt+1 770 DO i = nxlg, nxrg 771 u_p(k,j,nx+1) = u(k,j,nx+1) - dt_3d * tsc(2) * c_u_m(k) * & 772 ( u(k,j,nx+1) - u(k,j,nx) ) * ddx 773 774 v_p(k,j,nx+1) = v(k,j,nx+1) - dt_3d * tsc(2) * c_v_m(k) * & 775 ( v(k,j,nx+1) - v(k,j,nx) ) * ddx 776 777 w_p(k,j,nx+1) = w(k,j,nx+1) - dt_3d * tsc(2) * c_w_m(k) * & 778 ( w(k,j,nx+1) - w(k,j,nx) ) * ddx 779 ENDDO 780 ENDDO 781 782 ! 783 !-- Bottom boundary at the outflow 784 IF ( ibc_uv_b == 0 ) THEN 785 u_p(nzb,:,nx+1) = 0.0 786 v_p(nzb,:,nx+1) = 0.0 787 ELSE 788 u_p(nzb,:,nx+1) = u_p(nzb+1,:,nx+1) 789 v_p(nzb,:,nx+1) = v_p(nzb+1,:,nx+1) 790 ENDIF 791 w_p(nzb,:,nx+1) = 0.0 792 793 ! 794 !-- Top boundary at the outflow 795 IF ( ibc_uv_t == 0 ) THEN 796 u_p(nzt+1,:,nx+1) = u_init(nzt+1) 797 v_p(nzt+1,:,nx+1) = v_init(nzt+1) 798 ELSE 799 u_p(nzt+1,:,nx+1) = u_p(nzt,:,nx+1) 800 v_p(nzt+1,:,nx+1) = v_p(nzt,:,nx+1) 801 ENDIF 802 w(nzt:nzt+1,:,nx+1) = 0.0 803 804 ENDIF 619 805 620 806 ENDIF 621 807 622 623 808 END SUBROUTINE boundary_conds -
palm/trunk/SOURCE/check_parameters.f90
r965 r978 4 4 ! Current revisions: 5 5 ! ----------------- 6 ! 6 ! setting of bc_lr/ns_dirneu/neudir 7 ! outflow damping layer removed 8 ! check for z0h* 9 ! check for pt_damping_width 7 10 ! 8 11 ! Former revisions: … … 1379 1382 !-- Lateral boundary conditions 1380 1383 IF ( bc_lr /= 'cyclic' .AND. bc_lr /= 'dirichlet/radiation' .AND. & 1381 bc_lr /= 'radiation/dirichlet' ) THEN 1384 bc_lr /= 'radiation/dirichlet' .AND. bc_lr /= 'dirichlet/neumann' & 1385 .AND. bc_lr /= 'neumann/dirichlet' ) THEN 1382 1386 message_string = 'unknown boundary condition: bc_lr = "' // & 1383 1387 TRIM( bc_lr ) // '"' … … 1385 1389 ENDIF 1386 1390 IF ( bc_ns /= 'cyclic' .AND. bc_ns /= 'dirichlet/radiation' .AND. & 1387 bc_ns /= 'radiation/dirichlet' ) THEN 1391 bc_ns /= 'radiation/dirichlet' .AND. bc_ns /= 'dirichlet/neumann' & 1392 .AND. bc_ns /= 'neumann/dirichlet' ) THEN 1388 1393 message_string = 'unknown boundary condition: bc_ns = "' // & 1389 1394 TRIM( bc_ns ) // '"' … … 1396 1401 IF ( bc_lr == 'dirichlet/radiation' ) bc_lr_dirrad = .TRUE. 1397 1402 IF ( bc_lr == 'radiation/dirichlet' ) bc_lr_raddir = .TRUE. 1403 IF ( bc_lr == 'dirichlet/neumann' ) bc_lr_dirneu = .TRUE. 1404 IF ( bc_lr == 'neumann/dirichlet' ) bc_lr_neudir = .TRUE. 1398 1405 IF ( bc_ns /= 'cyclic' ) bc_ns_cyc = .FALSE. 1399 1406 IF ( bc_ns == 'dirichlet/radiation' ) bc_ns_dirrad = .TRUE. 1400 1407 IF ( bc_ns == 'radiation/dirichlet' ) bc_ns_raddir = .TRUE. 1408 IF ( bc_ns == 'dirichlet/neumann' ) bc_ns_dirneu = .TRUE. 1409 IF ( bc_ns == 'neumann/dirichlet' ) bc_ns_neudir = .TRUE. 1401 1410 1402 1411 ! … … 2659 2668 unit = 'psu' 2660 2669 2661 CASE ( 'u*', 't*', 'lwp*', 'pra*', 'prr*', 'qsws*', 'shf*', 'z0*' )2670 CASE ( 'u*', 't*', 'lwp*', 'pra*', 'prr*', 'qsws*', 'shf*', 'z0*', 'z0h*' ) 2662 2671 IF ( k == 0 .OR. data_output(i)(ilen-2:ilen) /= '_xy' ) THEN 2663 2672 message_string = 'illegal value for data_output: "' // & … … 2700 2709 IF ( TRIM( var ) == 'u*' ) unit = 'm/s' 2701 2710 IF ( TRIM( var ) == 'z0*' ) unit = 'm' 2711 IF ( TRIM( var ) == 'z0h*' ) unit = 'm' 2702 2712 2703 2713 … … 2964 2974 2965 2975 ! 2966 !-- In case of non-cyclic lateral boundaries, set the default maximum value 2967 !-- for the horizontal diffusivity used within the outflow damping layer, 2968 !-- and check/set the width of the damping layer 2976 !-- In case of non-cyclic lateral boundaries and a damping layer for the 2977 !-- potential temperature, check the width of the damping layer 2969 2978 IF ( bc_lr /= 'cyclic' ) THEN 2970 IF ( km_damp_max == -1.0 ) THEN 2971 km_damp_max = 0.5 * dx 2972 ENDIF 2973 IF ( outflow_damping_width == -1.0 ) THEN 2974 outflow_damping_width = MIN( 20, nx/2 ) 2975 ENDIF 2976 IF ( outflow_damping_width <= 0 .OR. outflow_damping_width > nx ) THEN 2977 message_string = 'outflow_damping width out of range' 2979 IF ( pt_damping_width < 0.0 .OR. pt_damping_width > REAL( nx * dx ) ) THEN 2980 message_string = 'pt_damping_width out of range' 2978 2981 CALL message( 'check_parameters', 'PA0124', 1, 2, 0, 6, 0 ) 2979 2982 ENDIF … … 2981 2984 2982 2985 IF ( bc_ns /= 'cyclic' ) THEN 2983 IF ( km_damp_max == -1.0 ) THEN 2984 km_damp_max = 0.5 * dy 2985 ENDIF 2986 IF ( outflow_damping_width == -1.0 ) THEN 2987 outflow_damping_width = MIN( 20, ny/2 ) 2988 ENDIF 2989 IF ( outflow_damping_width <= 0 .OR. outflow_damping_width > ny ) THEN 2990 message_string = 'outflow_damping width out of range' 2986 IF ( pt_damping_width < 0.0 .OR. pt_damping_width > REAL( ny * dy ) ) THEN 2987 message_string = 'pt_damping_width out of range' 2991 2988 CALL message( 'check_parameters', 'PA0124', 1, 2, 0, 6, 0 ) 2992 2989 ENDIF … … 3108 3105 ENDIF 3109 3106 3110 IF ( bc_lr == 'radiation/dirichlet' ) THEN3107 IF ( bc_lr == 'radiation/dirichlet' .OR. bc_lr == 'neumann/dirichlet' ) THEN 3111 3108 dist_nxr = nx - inflow_disturbance_begin 3112 3109 dist_nxl(1) = nx - inflow_disturbance_end 3113 ELSEIF ( bc_lr == 'dirichlet/radiation' ) THEN3110 ELSEIF ( bc_lr == 'dirichlet/radiation' .OR. bc_lr == 'dirichlet/neumann' ) THEN 3114 3111 dist_nxl = inflow_disturbance_begin 3115 3112 dist_nxr(1) = inflow_disturbance_end 3116 3113 ENDIF 3117 IF ( bc_ns == 'dirichlet/radiation' ) THEN3114 IF ( bc_ns == 'dirichlet/radiation' .OR. bc_ns == 'dirichlet/neumann' ) THEN 3118 3115 dist_nyn = ny - inflow_disturbance_begin 3119 3116 dist_nys(1) = ny - inflow_disturbance_end 3120 ELSEIF ( bc_ns == 'radiation/dirichlet' ) THEN3117 ELSEIF ( bc_ns == 'radiation/dirichlet' .OR. bc_ns == 'neumann/dirichlet' ) THEN 3121 3118 dist_nys = inflow_disturbance_begin 3122 3119 dist_nyn(1) = inflow_disturbance_end … … 3126 3123 !-- A turbulent inflow requires Dirichlet conditions at the respective inflow 3127 3124 !-- boundary (so far, a turbulent inflow is realized from the left side only) 3128 IF ( turbulent_inflow .AND. bc_lr /= 'dirichlet/radiation' ) THEN3125 IF ( turbulent_inflow .AND. bc_lr /= 'dirichlet/radiation' .AND. bc_lr /= 'dirichlet/neumann' ) THEN 3129 3126 message_string = 'turbulent_inflow = .T. requires a Dirichlet ' // & 3130 3127 'condition at the inflow boundary' -
palm/trunk/SOURCE/data_output_2d.f90
r791 r978 4 4 ! Current revisions: 5 5 ! ----------------- 6 ! 6 ! +z0h 7 7 ! 8 8 ! Former revisions: … … 594 594 DO j = nysg, nyng 595 595 local_pf(i,j,nzb+1) = z0_av(j,i) 596 ENDDO 597 ENDDO 598 ENDIF 599 resorted = .TRUE. 600 two_d = .TRUE. 601 level_z(nzb+1) = zu(nzb+1) 602 603 CASE ( 'z0h*_xy' ) ! 2d-array 604 IF ( av == 0 ) THEN 605 DO i = nxlg, nxrg 606 DO j = nysg, nyng 607 local_pf(i,j,nzb+1) = z0h(j,i) 608 ENDDO 609 ENDDO 610 ELSE 611 DO i = nxlg, nxrg 612 DO j = nysg, nyng 613 local_pf(i,j,nzb+1) = z0h_av(j,i) 596 614 ENDDO 597 615 ENDDO -
palm/trunk/SOURCE/diffusion_u.f90
r668 r978 4 4 ! Current revisions: 5 5 ! ----------------- 6 ! outflow damping layer removed 7 ! kmym_x/_y and kmyp_x/_y change to kmym and kmyp 6 8 ! 7 9 ! Former revisions: … … 63 65 ! Call for all grid points 64 66 !------------------------------------------------------------------------------! 65 SUBROUTINE diffusion_u( ddzu, ddzw, km, km_damp_y, tend, u, usws, uswst, & 66 v, w ) 67 SUBROUTINE diffusion_u( ddzu, ddzw, km, tend, u, usws, uswst, v, w ) 67 68 68 69 USE control_parameters … … 73 74 74 75 INTEGER :: i, j, k 75 REAL :: kmym _x, kmym_y, kmyp_x, kmyp_y, kmzm, kmzp76 REAL :: ddzu(1:nzt+1), ddzw(1:nzt+1) , km_damp_y(nysg:nyng)76 REAL :: kmym, kmyp, kmzm, kmzp 77 REAL :: ddzu(1:nzt+1), ddzw(1:nzt+1) 77 78 REAL :: tend(nzb:nzt+1,nysg:nyng,nxlg:nxrg) 78 79 REAL, DIMENSION(:,:), POINTER :: usws, uswst … … 95 96 ! 96 97 !-- Interpolate eddy diffusivities on staggered gridpoints 97 kmyp_x = 0.25 * & 98 ( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) ) 99 kmym_x = 0.25 * & 100 ( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) ) 101 kmyp_y = kmyp_x 102 kmym_y = kmym_x 103 ! 104 !-- Increase diffusion at the outflow boundary in case of 105 !-- non-cyclic lateral boundaries. Damping is only needed for 106 !-- velocity components parallel to the outflow boundary in 107 !-- the direction normal to the outflow boundary. 108 IF ( .NOT. bc_ns_cyc ) THEN 109 kmyp_y = MAX( kmyp_y, km_damp_y(j) ) 110 kmym_y = MAX( kmym_y, km_damp_y(j) ) 111 ENDIF 98 kmyp = 0.25 * & 99 ( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) ) 100 kmym = 0.25 * & 101 ( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) ) 112 102 113 103 tend(k,j,i) = tend(k,j,i) & … … 116 106 & - km(k,j,i-1) * ( u(k,j,i) - u(k,j,i-1) ) & 117 107 & ) * ddx2 & 118 & + ( kmyp _y * ( u(k,j+1,i) - u(k,j,i) ) * ddy&119 & + kmyp _x * ( v(k,j+1,i) - v(k,j+1,i-1) ) * ddx&120 & - kmym _y * ( u(k,j,i) - u(k,j-1,i) ) * ddy&121 & - kmym _x * ( v(k,j,i) - v(k,j,i-1) ) * ddx&108 & + ( kmyp * ( u(k,j+1,i) - u(k,j,i) ) * ddy & 109 & + kmyp * ( v(k,j+1,i) - v(k,j+1,i-1) ) * ddx & 110 & - kmym * ( u(k,j,i) - u(k,j-1,i) ) * ddy & 111 & - kmym * ( v(k,j,i) - v(k,j,i-1) ) * ddx & 122 112 & ) * ddy 123 113 ENDDO … … 128 118 129 119 DO k = nzb_u_inner(j,i)+1, nzb_u_outer(j,i) 130 kmyp_x = 0.25 * & 131 ( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) ) 132 kmym_x = 0.25 * & 133 ( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) ) 134 kmyp_y = kmyp_x 135 kmym_y = kmym_x 136 ! 137 !-- Increase diffusion at the outflow boundary in case of 138 !-- non-cyclic lateral boundaries. Damping is only needed for 139 !-- velocity components parallel to the outflow boundary in 140 !-- the direction normal to the outflow boundary. 141 IF ( .NOT. bc_ns_cyc ) THEN 142 kmyp_y = MAX( kmyp_y, km_damp_y(j) ) 143 kmym_y = MAX( kmym_y, km_damp_y(j) ) 144 ENDIF 120 kmyp = 0.25 * & 121 ( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) ) 122 kmym = 0.25 * & 123 ( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) ) 145 124 146 125 tend(k,j,i) = tend(k,j,i) & … … 150 129 ) * ddx2 & 151 130 + ( fyp(j,i) * ( & 152 kmyp _y * ( u(k,j+1,i) - u(k,j,i) ) * ddy&153 + kmyp _x * ( v(k,j+1,i) - v(k,j+1,i-1) ) * ddx&131 kmyp * ( u(k,j+1,i) - u(k,j,i) ) * ddy & 132 + kmyp * ( v(k,j+1,i) - v(k,j+1,i-1) ) * ddx & 154 133 ) & 155 134 - fym(j,i) * ( & 156 kmym _y * ( u(k,j,i) - u(k,j-1,i) ) * ddy&157 + kmym _x * ( v(k,j,i) - v(k,j,i-1) ) * ddx&135 kmym * ( u(k,j,i) - u(k,j-1,i) ) * ddy & 136 + kmym * ( v(k,j,i) - v(k,j,i-1) ) * ddx & 158 137 ) & 159 138 + wall_u(j,i) * usvs(k,j,i) & … … 239 218 ! Call for grid point i,j 240 219 !------------------------------------------------------------------------------! 241 SUBROUTINE diffusion_u_ij( i, j, ddzu, ddzw, km, km_damp_y,tend, u, usws, &220 SUBROUTINE diffusion_u_ij( i, j, ddzu, ddzw, km, tend, u, usws, & 242 221 uswst, v, w ) 243 222 … … 249 228 250 229 INTEGER :: i, j, k 251 REAL :: kmym _x, kmym_y, kmyp_x, kmyp_y, kmzm, kmzp252 REAL :: ddzu(1:nzt+1), ddzw(1:nzt+1) , km_damp_y(nysg:nyng)230 REAL :: kmym, kmyp, kmzm, kmzp 231 REAL :: ddzu(1:nzt+1), ddzw(1:nzt+1) 253 232 REAL :: tend(nzb:nzt+1,nysg:nyng,nxlg:nxrg) 254 233 REAL, DIMENSION(nzb:nzt+1) :: usvs … … 261 240 ! 262 241 !-- Interpolate eddy diffusivities on staggered gridpoints 263 kmyp_x = 0.25 * ( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) ) 264 kmym_x = 0.25 * ( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) ) 265 kmyp_y = kmyp_x 266 kmym_y = kmym_x 267 268 ! 269 !-- Increase diffusion at the outflow boundary in case of non-cyclic 270 !-- lateral boundaries. Damping is only needed for velocity components 271 !-- parallel to the outflow boundary in the direction normal to the 272 !-- outflow boundary. 273 IF ( .NOT. bc_ns_cyc ) THEN 274 kmyp_y = MAX( kmyp_y, km_damp_y(j) ) 275 kmym_y = MAX( kmym_y, km_damp_y(j) ) 276 ENDIF 242 kmyp = 0.25 * ( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) ) 243 kmym = 0.25 * ( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) ) 277 244 278 245 tend(k,j,i) = tend(k,j,i) & … … 281 248 & - km(k,j,i-1) * ( u(k,j,i) - u(k,j,i-1) ) & 282 249 & ) * ddx2 & 283 & + ( kmyp _y * ( u(k,j+1,i) - u(k,j,i) ) * ddy&284 & + kmyp _x * ( v(k,j+1,i) - v(k,j+1,i-1) ) * ddx&285 & - kmym _y * ( u(k,j,i) - u(k,j-1,i) ) * ddy&286 & - kmym _x * ( v(k,j,i) - v(k,j,i-1) ) * ddx&250 & + ( kmyp * ( u(k,j+1,i) - u(k,j,i) ) * ddy & 251 & + kmyp * ( v(k,j+1,i) - v(k,j+1,i-1) ) * ddx & 252 & - kmym * ( u(k,j,i) - u(k,j-1,i) ) * ddy & 253 & - kmym * ( v(k,j,i) - v(k,j,i-1) ) * ddx & 287 254 & ) * ddy 288 255 ENDDO … … 298 265 299 266 DO k = nzb_u_inner(j,i)+1, nzb_u_outer(j,i) 300 kmyp_x = 0.25 * ( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) ) 301 kmym_x = 0.25 * ( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) ) 302 kmyp_y = kmyp_x 303 kmym_y = kmym_x 304 ! 305 !-- Increase diffusion at the outflow boundary in case of 306 !-- non-cyclic lateral boundaries. Damping is only needed for 307 !-- velocity components parallel to the outflow boundary in 308 !-- the direction normal to the outflow boundary. 309 IF ( .NOT. bc_ns_cyc ) THEN 310 kmyp_y = MAX( kmyp_y, km_damp_y(j) ) 311 kmym_y = MAX( kmym_y, km_damp_y(j) ) 312 ENDIF 267 kmyp = 0.25 * ( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) ) 268 kmym = 0.25 * ( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) ) 313 269 314 270 tend(k,j,i) = tend(k,j,i) & … … 318 274 ) * ddx2 & 319 275 + ( fyp(j,i) * ( & 320 kmyp _y * ( u(k,j+1,i) - u(k,j,i) ) * ddy&321 + kmyp _x * ( v(k,j+1,i) - v(k,j+1,i-1) ) * ddx&276 kmyp * ( u(k,j+1,i) - u(k,j,i) ) * ddy & 277 + kmyp * ( v(k,j+1,i) - v(k,j+1,i-1) ) * ddx & 322 278 ) & 323 279 - fym(j,i) * ( & 324 kmym _y * ( u(k,j,i) - u(k,j-1,i) ) * ddy&325 + kmym _x * ( v(k,j,i) - v(k,j,i-1) ) * ddx&280 kmym * ( u(k,j,i) - u(k,j-1,i) ) * ddy & 281 + kmym * ( v(k,j,i) - v(k,j,i-1) ) * ddx & 326 282 ) & 327 283 + wall_u(j,i) * usvs(k) & -
palm/trunk/SOURCE/diffusion_v.f90
r668 r978 4 4 ! Current revisions: 5 5 ! ----------------- 6 ! outflow damping layer removed 7 ! kmxm_x/_y and kmxp_x/_y change to kmxm and kmxp 6 8 ! 7 9 ! Former revisions: … … 61 63 ! Call for all grid points 62 64 !------------------------------------------------------------------------------! 63 SUBROUTINE diffusion_v( ddzu, ddzw, km, km_damp_x, tend, u, v, vsws, & 64 vswst, w ) 65 SUBROUTINE diffusion_v( ddzu, ddzw, km, tend, u, v, vsws, vswst, w ) 65 66 66 67 USE control_parameters … … 71 72 72 73 INTEGER :: i, j, k 73 REAL :: kmxm _x, kmxm_y, kmxp_x, kmxp_y, kmzm, kmzp74 REAL :: ddzu(1:nzt+1), ddzw(1:nzt+1) , km_damp_x(nxlg:nxrg)74 REAL :: kmxm, kmxp, kmzm, kmzp 75 REAL :: ddzu(1:nzt+1), ddzw(1:nzt+1) 75 76 REAL :: tend(nzb:nzt+1,nysg:nyng,nxlg:nxrg) 76 77 REAL, DIMENSION(:,:), POINTER :: vsws, vswst … … 93 94 ! 94 95 !-- Interpolate eddy diffusivities on staggered gridpoints 95 kmxp_x = 0.25 * & 96 ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) ) 97 kmxm_x = 0.25 * & 98 ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) ) 99 kmxp_y = kmxp_x 100 kmxm_y = kmxm_x 101 ! 102 !-- Increase diffusion at the outflow boundary in case of 103 !-- non-cyclic lateral boundaries. Damping is only needed for 104 !-- velocity components parallel to the outflow boundary in 105 !-- the direction normal to the outflow boundary. 106 IF ( .NOT. bc_lr_cyc ) THEN 107 kmxp_x = MAX( kmxp_x, km_damp_x(i) ) 108 kmxm_x = MAX( kmxm_x, km_damp_x(i) ) 109 ENDIF 96 kmxp = 0.25 * & 97 ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) ) 98 kmxm = 0.25 * & 99 ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) ) 110 100 111 101 tend(k,j,i) = tend(k,j,i) & 112 & + ( kmxp _x * ( v(k,j,i+1) - v(k,j,i) ) * ddx&113 & + kmxp _y * ( u(k,j,i+1) - u(k,j-1,i+1) ) * ddy&114 & - kmxm _x * ( v(k,j,i) - v(k,j,i-1) ) * ddx&115 & - kmxm _y * ( u(k,j,i) - u(k,j-1,i) ) * ddy&102 & + ( kmxp * ( v(k,j,i+1) - v(k,j,i) ) * ddx & 103 & + kmxp * ( u(k,j,i+1) - u(k,j-1,i+1) ) * ddy & 104 & - kmxm * ( v(k,j,i) - v(k,j,i-1) ) * ddx & 105 & - kmxm * ( u(k,j,i) - u(k,j-1,i) ) * ddy & 116 106 & ) * ddx & 117 107 & + 2.0 * ( & … … 126 116 127 117 DO k = nzb_v_inner(j,i)+1, nzb_v_outer(j,i) 128 kmxp_x = 0.25 * & 129 ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) ) 130 kmxm_x = 0.25 * & 131 ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) ) 132 kmxp_y = kmxp_x 133 kmxm_y = kmxm_x 134 ! 135 !-- Increase diffusion at the outflow boundary in case of 136 !-- non-cyclic lateral boundaries. Damping is only needed for 137 !-- velocity components parallel to the outflow boundary in 138 !-- the direction normal to the outflow boundary. 139 IF ( .NOT. bc_lr_cyc ) THEN 140 kmxp_x = MAX( kmxp_x, km_damp_x(i) ) 141 kmxm_x = MAX( kmxm_x, km_damp_x(i) ) 142 ENDIF 143 118 kmxp = 0.25 * & 119 ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) ) 120 kmxm = 0.25 * & 121 ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) ) 122 144 123 tend(k,j,i) = tend(k,j,i) & 145 124 + 2.0 * ( & … … 148 127 ) * ddy2 & 149 128 + ( fxp(j,i) * ( & 150 kmxp _x * ( v(k,j,i+1) - v(k,j,i) ) * ddx&151 + kmxp _y * ( u(k,j,i+1) - u(k,j-1,i+1) ) * ddy&129 kmxp * ( v(k,j,i+1) - v(k,j,i) ) * ddx & 130 + kmxp * ( u(k,j,i+1) - u(k,j-1,i+1) ) * ddy & 152 131 ) & 153 132 - fxm(j,i) * ( & 154 kmxm _x * ( v(k,j,i) - v(k,j,i-1) ) * ddx&155 + kmxm _y * ( u(k,j,i) - u(k,j-1,i) ) * ddy&133 kmxm * ( v(k,j,i) - v(k,j,i-1) ) * ddx & 134 + kmxm * ( u(k,j,i) - u(k,j-1,i) ) * ddy & 156 135 ) & 157 136 + wall_v(j,i) * vsus(k,j,i) & … … 237 216 ! Call for grid point i,j 238 217 !------------------------------------------------------------------------------! 239 SUBROUTINE diffusion_v_ij( i, j, ddzu, ddzw, km, km_damp_x,tend, u, v, &218 SUBROUTINE diffusion_v_ij( i, j, ddzu, ddzw, km, tend, u, v, & 240 219 vsws, vswst, w ) 241 220 … … 247 226 248 227 INTEGER :: i, j, k 249 REAL :: kmxm _x, kmxm_y, kmxp_x, kmxp_y, kmzm, kmzp250 REAL :: ddzu(1:nzt+1), ddzw(1:nzt+1) , km_damp_x(nxlg:nxrg)228 REAL :: kmxm, kmxp, kmzm, kmzp 229 REAL :: ddzu(1:nzt+1), ddzw(1:nzt+1) 251 230 REAL :: tend(nzb:nzt+1,nysg:nyng,nxlg:nxrg) 252 231 REAL, DIMENSION(nzb:nzt+1) :: vsus … … 259 238 ! 260 239 !-- Interpolate eddy diffusivities on staggered gridpoints 261 kmxp_x = 0.25 * ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) ) 262 kmxm_x = 0.25 * ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) ) 263 kmxp_y = kmxp_x 264 kmxm_y = kmxm_x 265 ! 266 !-- Increase diffusion at the outflow boundary in case of non-cyclic 267 !-- lateral boundaries. Damping is only needed for velocity components 268 !-- parallel to the outflow boundary in the direction normal to the 269 !-- outflow boundary. 270 IF ( .NOT. bc_lr_cyc ) THEN 271 kmxp_x = MAX( kmxp_x, km_damp_x(i) ) 272 kmxm_x = MAX( kmxm_x, km_damp_x(i) ) 273 ENDIF 240 kmxp = 0.25 * ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) ) 241 kmxm = 0.25 * ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) ) 274 242 275 243 tend(k,j,i) = tend(k,j,i) & 276 & + ( kmxp _x * ( v(k,j,i+1) - v(k,j,i) ) * ddx&277 & + kmxp _y * ( u(k,j,i+1) - u(k,j-1,i+1) ) * ddy&278 & - kmxm _x * ( v(k,j,i) - v(k,j,i-1) ) * ddx&279 & - kmxm _y * ( u(k,j,i) - u(k,j-1,i) ) * ddy&244 & + ( kmxp * ( v(k,j,i+1) - v(k,j,i) ) * ddx & 245 & + kmxp * ( u(k,j,i+1) - u(k,j-1,i+1) ) * ddy & 246 & - kmxm * ( v(k,j,i) - v(k,j,i-1) ) * ddx & 247 & - kmxm * ( u(k,j,i) - u(k,j-1,i) ) * ddy & 280 248 & ) * ddx & 281 249 & + 2.0 * ( & … … 295 263 296 264 DO k = nzb_v_inner(j,i)+1, nzb_v_outer(j,i) 297 kmxp_x = 0.25 * & 298 ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) ) 299 kmxm_x = 0.25 * & 300 ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) ) 301 kmxp_y = kmxp_x 302 kmxm_y = kmxm_x 303 ! 304 !-- Increase diffusion at the outflow boundary in case of 305 !-- non-cyclic lateral boundaries. Damping is only needed for 306 !-- velocity components parallel to the outflow boundary in 307 !-- the direction normal to the outflow boundary. 308 IF ( .NOT. bc_lr_cyc ) THEN 309 kmxp_x = MAX( kmxp_x, km_damp_x(i) ) 310 kmxm_x = MAX( kmxm_x, km_damp_x(i) ) 311 ENDIF 265 kmxp = 0.25 * & 266 ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) ) 267 kmxm = 0.25 * & 268 ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) ) 312 269 313 270 tend(k,j,i) = tend(k,j,i) & … … 317 274 ) * ddy2 & 318 275 + ( fxp(j,i) * ( & 319 kmxp _x * ( v(k,j,i+1) - v(k,j,i) ) * ddx&320 + kmxp _y * ( u(k,j,i+1) - u(k,j-1,i+1) ) * ddy&276 kmxp * ( v(k,j,i+1) - v(k,j,i) ) * ddx & 277 + kmxp * ( u(k,j,i+1) - u(k,j-1,i+1) ) * ddy & 321 278 ) & 322 279 - fxm(j,i) * ( & 323 kmxm _x * ( v(k,j,i) - v(k,j,i-1) ) * ddx&324 + kmxm _y * ( u(k,j,i) - u(k,j-1,i) ) * ddy&280 kmxm * ( v(k,j,i) - v(k,j,i-1) ) * ddx & 281 + kmxm * ( u(k,j,i) - u(k,j-1,i) ) * ddy & 325 282 ) & 326 283 + wall_v(j,i) * vsus(k) & -
palm/trunk/SOURCE/diffusion_w.f90
r668 r978 4 4 ! Current revisions: 5 5 ! ----------------- 6 ! outflow damping layer removed 7 ! kmxm_x/_z and kmxp_x/_z change to kmxm and kmxp 8 ! kmym_y/_z and kmyp_y/_z change to kmym and kmyp 6 9 ! 7 10 ! Former revisions: … … 54 57 ! Call for all grid points 55 58 !------------------------------------------------------------------------------! 56 SUBROUTINE diffusion_w( ddzu, ddzw, km, km_damp_x, km_damp_y, tend, u, v, & 57 w ) 59 SUBROUTINE diffusion_w( ddzu, ddzw, km, tend, u, v, w ) 58 60 59 61 USE control_parameters … … 64 66 65 67 INTEGER :: i, j, k 66 REAL :: kmxm_x, kmxm_z, kmxp_x, kmxp_z, kmym_y, kmym_z, kmyp_y, & 67 kmyp_z 68 REAL :: ddzu(1:nzt+1), ddzw(1:nzt+1), km_damp_x(nxlg:nxrg), & 69 km_damp_y(nysg:nyng) 68 REAL :: kmxm, kmxp, kmym, kmyp 69 REAL :: ddzu(1:nzt+1), ddzw(1:nzt+1) 70 70 REAL :: tend(nzb:nzt+1,nysg:nyng,nxlg:nxrg) 71 71 REAL, DIMENSION(:,:,:), POINTER :: km, u, v, w … … 88 88 ! 89 89 !-- Interpolate eddy diffusivities on staggered gridpoints 90 kmxp_x = 0.25 * & 91 ( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) ) 92 kmxm_x = 0.25 * & 93 ( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) ) 94 kmxp_z = kmxp_x 95 kmxm_z = kmxm_x 96 kmyp_y = 0.25 * & 97 ( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) ) 98 kmym_y = 0.25 * & 99 ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) ) 100 kmyp_z = kmyp_y 101 kmym_z = kmym_y 102 ! 103 !-- Increase diffusion at the outflow boundary in case of 104 !-- non-cyclic lateral boundaries. Damping is only needed for 105 !-- velocity components parallel to the outflow boundary in 106 !-- the direction normal to the outflow boundary. 107 IF ( .NOT. bc_lr_cyc ) THEN 108 kmxp_x = MAX( kmxp_x, km_damp_x(i) ) 109 kmxm_x = MAX( kmxm_x, km_damp_x(i) ) 110 ENDIF 111 IF ( .NOT. bc_ns_cyc ) THEN 112 kmyp_y = MAX( kmyp_y, km_damp_y(j) ) 113 kmym_y = MAX( kmym_y, km_damp_y(j) ) 114 ENDIF 90 kmxp = 0.25 * & 91 ( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) ) 92 kmxm = 0.25 * & 93 ( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) ) 94 kmyp = 0.25 * & 95 ( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) ) 96 kmym = 0.25 * & 97 ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) ) 115 98 116 99 tend(k,j,i) = tend(k,j,i) & 117 & + ( kmxp _x * ( w(k,j,i+1) - w(k,j,i) ) * ddx&118 & + kmxp _z * ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1)&119 & - kmxm _x * ( w(k,j,i) - w(k,j,i-1) ) * ddx&120 & - kmxm _z * ( u(k+1,j,i) - u(k,j,i) ) * ddzu(k+1)&100 & + ( kmxp * ( w(k,j,i+1) - w(k,j,i) ) * ddx & 101 & + kmxp * ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1) & 102 & - kmxm * ( w(k,j,i) - w(k,j,i-1) ) * ddx & 103 & - kmxm * ( u(k+1,j,i) - u(k,j,i) ) * ddzu(k+1) & 121 104 & ) * ddx & 122 & + ( kmyp _y * ( w(k,j+1,i) - w(k,j,i) ) * ddy&123 & + kmyp _z * ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1)&124 & - kmym _y * ( w(k,j,i) - w(k,j-1,i) ) * ddy&125 & - kmym _z * ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1)&105 & + ( kmyp * ( w(k,j+1,i) - w(k,j,i) ) * ddy & 106 & + kmyp * ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1) & 107 & - kmym * ( w(k,j,i) - w(k,j-1,i) ) * ddy & 108 & - kmym * ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1) & 126 109 & ) * ddy & 127 110 & + 2.0 * ( & … … 138 121 ! 139 122 !-- Interpolate eddy diffusivities on staggered gridpoints 140 kmxp_x = 0.25 * & 141 ( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) ) 142 kmxm_x = 0.25 * & 143 ( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) ) 144 kmxp_z = kmxp_x 145 kmxm_z = kmxm_x 146 kmyp_y = 0.25 * & 147 ( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) ) 148 kmym_y = 0.25 * & 149 ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) ) 150 kmyp_z = kmyp_y 151 kmym_z = kmym_y 152 ! 153 !-- Increase diffusion at the outflow boundary in case of 154 !-- non-cyclic lateral boundaries. Damping is only needed for 155 !-- velocity components parallel to the outflow boundary in 156 !-- the direction normal to the outflow boundary. 157 IF ( .NOT. bc_lr_cyc ) THEN 158 kmxp_x = MAX( kmxp_x, km_damp_x(i) ) 159 kmxm_x = MAX( kmxm_x, km_damp_x(i) ) 160 ENDIF 161 IF ( .NOT. bc_ns_cyc ) THEN 162 kmyp_y = MAX( kmyp_y, km_damp_y(j) ) 163 kmym_y = MAX( kmym_y, km_damp_y(j) ) 164 ENDIF 123 kmxp = 0.25 * & 124 ( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) ) 125 kmxm = 0.25 * & 126 ( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) ) 127 kmyp = 0.25 * & 128 ( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) ) 129 kmym = 0.25 * & 130 ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) ) 165 131 166 132 tend(k,j,i) = tend(k,j,i) & 167 133 + ( fwxp(j,i) * ( & 168 kmxp _x * ( w(k,j,i+1) - w(k,j,i) ) * ddx&169 + kmxp _z * ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1)&134 kmxp * ( w(k,j,i+1) - w(k,j,i) ) * ddx & 135 + kmxp * ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1) & 170 136 ) & 171 137 - fwxm(j,i) * ( & 172 kmxm _x * ( w(k,j,i) - w(k,j,i-1) ) * ddx&173 + kmxm _z * ( u(k+1,j,i) - u(k,j,i) ) * ddzu(k+1)&138 kmxm * ( w(k,j,i) - w(k,j,i-1) ) * ddx & 139 + kmxm * ( u(k+1,j,i) - u(k,j,i) ) * ddzu(k+1) & 174 140 ) & 175 141 + wall_w_x(j,i) * wsus(k,j,i) & 176 142 ) * ddx & 177 143 + ( fwyp(j,i) * ( & 178 kmyp _y * ( w(k,j+1,i) - w(k,j,i) ) * ddy&179 + kmyp _z * ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1)&144 kmyp * ( w(k,j+1,i) - w(k,j,i) ) * ddy & 145 + kmyp * ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1) & 180 146 ) & 181 147 - fwym(j,i) * ( & 182 kmym _y * ( w(k,j,i) - w(k,j-1,i) ) * ddy&183 + kmym _z * ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1)&148 kmym * ( w(k,j,i) - w(k,j-1,i) ) * ddy & 149 + kmym * ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1) & 184 150 ) & 185 151 + wall_w_y(j,i) * wsvs(k,j,i) & … … 201 167 ! Call for grid point i,j 202 168 !------------------------------------------------------------------------------! 203 SUBROUTINE diffusion_w_ij( i, j, ddzu, ddzw, km, km_damp_x, km_damp_y, & 204 tend, u, v, w ) 169 SUBROUTINE diffusion_w_ij( i, j, ddzu, ddzw, km, tend, u, v, w ) 205 170 206 171 USE control_parameters … … 211 176 212 177 INTEGER :: i, j, k 213 REAL :: kmxm_x, kmxm_z, kmxp_x, kmxp_z, kmym_y, kmym_z, kmyp_y, & 214 kmyp_z 215 REAL :: ddzu(1:nzt+1), ddzw(1:nzt+1), km_damp_x(nxlg:nxrg), & 216 km_damp_y(nysg:nyng) 178 REAL :: kmxm, kmxp, kmym, kmyp 179 REAL :: ddzu(1:nzt+1), ddzw(1:nzt+1) 217 180 REAL :: tend(nzb:nzt+1,nysg:nyng,nxlg:nxrg) 218 181 REAL, DIMENSION(nzb:nzt+1) :: wsus, wsvs … … 223 186 ! 224 187 !-- Interpolate eddy diffusivities on staggered gridpoints 225 kmxp_x = 0.25 * ( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) ) 226 kmxm_x = 0.25 * ( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) ) 227 kmxp_z = kmxp_x 228 kmxm_z = kmxm_x 229 kmyp_y = 0.25 * ( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) ) 230 kmym_y = 0.25 * ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) ) 231 kmyp_z = kmyp_y 232 kmym_z = kmym_y 233 ! 234 !-- Increase diffusion at the outflow boundary in case of non-cyclic 235 !-- lateral boundaries. Damping is only needed for velocity components 236 !-- parallel to the outflow boundary in the direction normal to the 237 !-- outflow boundary. 238 IF ( .NOT. bc_lr_cyc ) THEN 239 kmxp_x = MAX( kmxp_x, km_damp_x(i) ) 240 kmxm_x = MAX( kmxm_x, km_damp_x(i) ) 241 ENDIF 242 IF ( .NOT. bc_ns_cyc ) THEN 243 kmyp_y = MAX( kmyp_y, km_damp_y(j) ) 244 kmym_y = MAX( kmym_y, km_damp_y(j) ) 245 ENDIF 188 kmxp = 0.25 * ( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) ) 189 kmxm = 0.25 * ( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) ) 190 kmyp = 0.25 * ( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) ) 191 kmym = 0.25 * ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) ) 246 192 247 193 tend(k,j,i) = tend(k,j,i) & 248 & + ( kmxp _x * ( w(k,j,i+1) - w(k,j,i) ) * ddx&249 & + kmxp _z * ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1)&250 & - kmxm _x * ( w(k,j,i) - w(k,j,i-1) ) * ddx&251 & - kmxm _z * ( u(k+1,j,i) - u(k,j,i) ) * ddzu(k+1)&194 & + ( kmxp * ( w(k,j,i+1) - w(k,j,i) ) * ddx & 195 & + kmxp * ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1) & 196 & - kmxm * ( w(k,j,i) - w(k,j,i-1) ) * ddx & 197 & - kmxm * ( u(k+1,j,i) - u(k,j,i) ) * ddzu(k+1) & 252 198 & ) * ddx & 253 & + ( kmyp _y * ( w(k,j+1,i) - w(k,j,i) ) * ddy&254 & + kmyp _z * ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1)&255 & - kmym _y * ( w(k,j,i) - w(k,j-1,i) ) * ddy&256 & - kmym _z * ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1)&199 & + ( kmyp * ( w(k,j+1,i) - w(k,j,i) ) * ddy & 200 & + kmyp * ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1) & 201 & - kmym * ( w(k,j,i) - w(k,j-1,i) ) * ddy & 202 & - kmym * ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1) & 257 203 & ) * ddy & 258 204 & + 2.0 * ( & … … 285 231 ! 286 232 !-- Interpolate eddy diffusivities on staggered gridpoints 287 kmxp_x = 0.25 * ( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) ) 288 kmxm_x = 0.25 * ( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) ) 289 kmxp_z = kmxp_x 290 kmxm_z = kmxm_x 291 kmyp_y = 0.25 * ( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) ) 292 kmym_y = 0.25 * ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) ) 293 kmyp_z = kmyp_y 294 kmym_z = kmym_y 295 ! 296 !-- Increase diffusion at the outflow boundary in case of 297 !-- non-cyclic lateral boundaries. Damping is only needed for 298 !-- velocity components parallel to the outflow boundary in 299 !-- the direction normal to the outflow boundary. 300 IF ( .NOT. bc_lr_cyc ) THEN 301 kmxp_x = MAX( kmxp_x, km_damp_x(i) ) 302 kmxm_x = MAX( kmxm_x, km_damp_x(i) ) 303 ENDIF 304 IF ( .NOT. bc_ns_cyc ) THEN 305 kmyp_y = MAX( kmyp_y, km_damp_y(j) ) 306 kmym_y = MAX( kmym_y, km_damp_y(j) ) 307 ENDIF 233 kmxp = 0.25 * ( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) ) 234 kmxm = 0.25 * ( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) ) 235 kmyp = 0.25 * ( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) ) 236 kmym = 0.25 * ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) ) 308 237 309 238 tend(k,j,i) = tend(k,j,i) & 310 239 + ( fwxp(j,i) * ( & 311 kmxp _x * ( w(k,j,i+1) - w(k,j,i) ) * ddx&312 + kmxp _z * ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1)&240 kmxp * ( w(k,j,i+1) - w(k,j,i) ) * ddx & 241 + kmxp * ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1) & 313 242 ) & 314 243 - fwxm(j,i) * ( & 315 kmxm _x * ( w(k,j,i) - w(k,j,i-1) ) * ddx&316 + kmxm _z * ( u(k+1,j,i) - u(k,j,i) ) * ddzu(k+1)&244 kmxm * ( w(k,j,i) - w(k,j,i-1) ) * ddx & 245 + kmxm * ( u(k+1,j,i) - u(k,j,i) ) * ddzu(k+1) & 317 246 ) & 318 247 + wall_w_x(j,i) * wsus(k) & 319 248 ) * ddx & 320 249 + ( fwyp(j,i) * ( & 321 kmyp _y * ( w(k,j+1,i) - w(k,j,i) ) * ddy&322 + kmyp _z * ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1)&250 kmyp * ( w(k,j+1,i) - w(k,j,i) ) * ddy & 251 + kmyp * ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1) & 323 252 ) & 324 253 - fwym(j,i) * ( & 325 kmym _y * ( w(k,j,i) - w(k,j-1,i) ) * ddy&326 + kmym _z * ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1)&254 kmym * ( w(k,j,i) - w(k,j-1,i) ) * ddy & 255 + kmym * ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1) & 327 256 ) & 328 257 + wall_w_y(j,i) * wsvs(k) & -
palm/trunk/SOURCE/header.f90
r965 r978 4 4 ! Current revisions: 5 5 ! ----------------- 6 ! 6 ! -km_damp_max, outflow_damping_width 7 ! +pt_damping_factor, pt_damping_width 8 ! +z0h 7 9 ! 8 10 ! Former revisions: … … 704 706 705 707 IF ( prandtl_layer ) THEN 706 WRITE ( io, 305 ) (zu(1)-zu(0)), roughness_length, kappa, & 708 WRITE ( io, 305 ) (zu(1)-zu(0)), roughness_length, & 709 z0h_factor*roughness_length, kappa, & 707 710 rif_min, rif_max 708 711 IF ( .NOT. constant_heatflux ) WRITE ( io, 308 ) … … 721 724 WRITE ( io, 317 ) bc_lr, bc_ns 722 725 IF ( .NOT. bc_lr_cyc .OR. .NOT. bc_ns_cyc ) THEN 723 WRITE ( io, 318 ) outflow_damping_width, km_damp_max726 WRITE ( io, 318 ) pt_damping_width, pt_damping_factor 724 727 IF ( turbulent_inflow ) THEN 725 728 WRITE ( io, 319 ) recycling_width, recycling_plane, & … … 1779 1782 305 FORMAT (//' Prandtl-Layer between bottom surface and first ', & 1780 1783 'computational u,v-level:'// & 1781 ' zp = ',F6.2,' m z0 = ',F6.4,' m kappa = ',F4.2/ & 1784 ' zp = ',F6.2,' m z0 = ',F6.4,' m z0h = ',F7.5,& 1785 ' m kappa = ',F4.2/ & 1782 1786 ' Rif value range: ',F6.2,' <= rif <=',F6.2) 1783 1787 306 FORMAT (' Predefined constant heatflux: ',F9.6,' K m/s') … … 1797 1801 ' left/right: ',A/ & 1798 1802 ' north/south: ',A) 1799 318 FORMAT (/' outflow damping layer width: ',I3,' gridpoints with km_', &1800 ' max =',F5.1,' m**2/s')1803 318 FORMAT (/' pt damping layer width = ',F7.2,' m, pt ', & 1804 'damping factor = ',F6.4) 1801 1805 319 FORMAT (' turbulence recycling at inflow switched on'/ & 1802 1806 ' width of recycling domain: ',F7.1,' m grid index: ',I4/ & -
palm/trunk/SOURCE/init_1d_model.f90
r668 r978 4 4 ! Current revisions: 5 5 ! ----------------- 6 ! 6 ! roughness length for scalar quantities z0h1d added 7 7 ! 8 8 ! Former revisions: … … 149 149 vsws1d = 0.0; vsws1d_m = 0.0 150 150 z01d = roughness_length 151 z0h1d = z0h_factor * z01d 151 152 IF ( humidity .OR. passive_scalar ) qs1d = 0.0 152 153 … … 445 446 !-- Stable stratification 446 447 ts1d = kappa * ( pt_init(nzb+1) - pt_init(nzb) ) / & 447 ( LOG( zu(nzb+1) / z0 1d ) + 5.0 * rif1d(nzb+1) * &448 ( zu(nzb+1) - z0 1d ) / zu(nzb+1) &448 ( LOG( zu(nzb+1) / z0h1d ) + 5.0 * rif1d(nzb+1) * & 449 ( zu(nzb+1) - z0h1d ) / zu(nzb+1) & 449 450 ) 450 451 ELSE … … 452 453 !-- Unstable stratification 453 454 a = SQRT( 1.0 - 16.0 * rif1d(nzb+1) ) 454 b = SQRT( 1.0 - 16.0 * rif1d(nzb+1) / zu(nzb+1) * z0 1d )455 b = SQRT( 1.0 - 16.0 * rif1d(nzb+1) / zu(nzb+1) * z0h1d ) 455 456 ! 456 457 !-- In the borderline case the formula for stable stratification … … 459 460 IF ( a == 0.0 .OR. b == 0.0 ) THEN 460 461 ts1d = kappa * ( pt_init(nzb+1) - pt_init(nzb) ) / & 461 ( LOG( zu(nzb+1) / z0 1d ) + 5.0 * rif1d(nzb+1) * &462 ( zu(nzb+1) - z0 1d ) / zu(nzb+1) &462 ( LOG( zu(nzb+1) / z0h1d ) + 5.0 * rif1d(nzb+1) * & 463 ( zu(nzb+1) - z0h1d ) / zu(nzb+1) & 463 464 ) 464 465 ELSE … … 578 579 !-- Stable stratification 579 580 qs1d = kappa * ( q_init(nzb+1) - q_init(nzb) ) / & 580 ( LOG( zu(nzb+1) / z0 1d ) + 5.0 * rif1d(nzb+1) * &581 ( zu(nzb+1) - z0 1d ) / zu(nzb+1) &581 ( LOG( zu(nzb+1) / z0h1d ) + 5.0 * rif1d(nzb+1) * & 582 ( zu(nzb+1) - z0h1d ) / zu(nzb+1) & 582 583 ) 583 584 ELSE … … 585 586 !-- Unstable stratification 586 587 a = SQRT( 1.0 - 16.0 * rif1d(nzb+1) ) 587 b = SQRT( 1.0 - 16.0 * rif1d(nzb+1) / zu(nzb+1) * z0 1d )588 b = SQRT( 1.0 - 16.0 * rif1d(nzb+1) / zu(nzb+1) * z0h1d ) 588 589 ! 589 590 !-- In the borderline case the formula for stable stratification … … 592 593 IF ( a == 1.0 .OR. b == 1.0 ) THEN 593 594 qs1d = kappa * ( q_init(nzb+1) - q_init(nzb) ) / & 594 ( LOG( zu(nzb+1) / z0 1d ) + 5.0 * rif1d(nzb+1) * &595 ( zu(nzb+1) - z0 1d ) / zu(nzb+1) &595 ( LOG( zu(nzb+1) / z0h1d ) + 5.0 * rif1d(nzb+1) * & 596 ( zu(nzb+1) - z0h1d ) / zu(nzb+1) & 596 597 ) 597 598 ELSE -
palm/trunk/SOURCE/init_3d_model.f90
r850 r978 7 7 ! Current revisions: 8 8 ! ------------------ 9 ! 9 ! outflow damping layer removed 10 ! roughness length for scalar quantites z0h added 11 ! damping zone for the potential temperatur in case of non-cyclic lateral 12 ! boundaries added 13 ! initialization of ptdf_x, ptdf_y 14 ! initialization of c_u_m, c_u_m_l, c_v_m, c_v_m_l, c_w_m, c_w_m_l 10 15 ! 11 16 ! Former revisions: … … 164 169 USE control_parameters 165 170 USE cpulog 171 USE grid_variables 166 172 USE indices 167 173 USE interfaces … … 211 217 sums_wsts_bc_l(nzb:nzt+1,0:statistic_regions), & 212 218 ts_value(dots_max,0:statistic_regions) ) 213 ALLOCATE( km_damp_x(nxlg:nxrg), km_damp_y(nysg:nyng) )219 ALLOCATE( ptdf_x(nxlg:nxrg), ptdf_y(nysg:nyng) ) 214 220 215 221 ALLOCATE( rif_1(nysg:nyng,nxlg:nxrg), shf_1(nysg:nyng,nxlg:nxrg), & … … 218 224 uswst_1(nysg:nyng,nxlg:nxrg), & 219 225 vsws_1(nysg:nyng,nxlg:nxrg), & 220 vswst_1(nysg:nyng,nxlg:nxrg), z0(nysg:nyng,nxlg:nxrg) ) 226 vswst_1(nysg:nyng,nxlg:nxrg), z0(nysg:nyng,nxlg:nxrg), & 227 z0h(nysg:nyng,nxlg:nxrg) ) 221 228 222 229 IF ( timestep_scheme(1:5) /= 'runge' ) THEN … … 413 420 c_w(nzb:nzt+1,nxlg:nxrg) ) 414 421 ENDIF 422 IF ( outflow_l .OR. outflow_r .OR. outflow_s .OR. outflow_n ) THEN 423 ALLOCATE( c_u_m_l(nzb:nzt+1), c_v_m_l(nzb:nzt+1), c_w_m_l(nzb:nzt+1) ) 424 ALLOCATE( c_u_m(nzb:nzt+1), c_v_m(nzb:nzt+1), c_w_m(nzb:nzt+1) ) 425 ENDIF 426 415 427 416 428 ! … … 847 859 848 860 z0 = roughness_length 861 z0h = z0h_factor * z0 849 862 850 863 IF ( .NOT. constant_heatflux ) THEN … … 1533 1546 1534 1547 ! 1535 !-- Initialize diffusivities used within the outflow damping layer in case of 1536 !-- non-cyclic lateral boundaries. A linear increase is assumed over the first 1537 !-- half of the width of the damping layer 1538 IF ( bc_lr_dirrad ) THEN 1539 1540 DO i = nxlg, nxrg 1541 IF ( i >= nx - outflow_damping_width ) THEN 1542 km_damp_x(i) = km_damp_max * MIN( 1.0, & 1543 ( i - ( nx - outflow_damping_width ) ) / & 1544 REAL( outflow_damping_width/2 ) & 1545 ) 1546 ELSE 1547 km_damp_x(i) = 0.0 1548 !-- Initialize damping zone for the potential temperature in case of 1549 !-- non-cyclic lateral boundaries. The damping zone has the maximum value 1550 !-- at the inflow boundary and decreases to zero at pt_damping_width. 1551 ptdf_x = 0.0 1552 ptdf_y = 0.0 1553 IF ( bc_lr_dirrad .OR. bc_lr_dirneu ) THEN 1554 DO i = nxl, nxr 1555 IF ( ( i * dx ) < pt_damping_width ) THEN 1556 ptdf_x(i) = pt_damping_factor * ( SIN( pi * 0.5 * & 1557 REAL( pt_damping_width - i * dx ) / ( & 1558 REAL( pt_damping_width ) ) ) )**2 1548 1559 ENDIF 1549 1560 ENDDO 1550 1551 ELSEIF ( bc_lr_raddir ) THEN 1552 1553 DO i = nxlg, nxrg 1554 IF ( i <= outflow_damping_width ) THEN 1555 km_damp_x(i) = km_damp_max * MIN( 1.0, & 1556 ( outflow_damping_width - i ) / & 1557 REAL( outflow_damping_width/2 ) & 1558 ) 1559 ELSE 1560 km_damp_x(i) = 0.0 1561 ELSEIF ( bc_lr_raddir .OR. bc_lr_neudir ) THEN 1562 DO i = nxl, nxr 1563 IF ( ( i * dx ) > ( nx * dx - pt_damping_width ) ) THEN 1564 ptdf_x(i) = pt_damping_factor * ( SIN( pi * 0.5 * & 1565 REAL( ( i - nx ) * dx + pt_damping_width ) / ( & 1566 REAL( pt_damping_width ) ) ) )**2 1567 ENDIF 1568 ENDDO 1569 ELSEIF ( bc_ns_dirrad .OR. bc_ns_dirneu ) THEN 1570 DO j = nys, nyn 1571 IF ( ( j * dy ) > ( ny * dy - pt_damping_width ) ) THEN 1572 ptdf_y(j) = pt_damping_factor * ( SIN( pi * 0.5 * & 1573 REAL( ( j - ny ) * dy + pt_damping_width ) / ( & 1574 REAL( pt_damping_width ) ) ) )**2 1575 ENDIF 1576 ENDDO 1577 ELSEIF ( bc_ns_raddir .OR. bc_ns_neudir ) THEN 1578 DO j = nys, nyn 1579 IF ( ( j * dy ) < pt_damping_width ) THEN 1580 ptdf_y(j) = pt_damping_factor * ( SIN( pi * 0.5 * & 1581 REAL( pt_damping_width - j * dy ) / ( & 1582 REAL( pt_damping_width ) ) ) )**2 1561 1583 ENDIF 1562 1584 ENDDO 1563 1564 ENDIF1565 1566 IF ( bc_ns_dirrad ) THEN1567 1568 DO j = nysg, nyng1569 IF ( j >= ny - outflow_damping_width ) THEN1570 km_damp_y(j) = km_damp_max * MIN( 1.0, &1571 ( j - ( ny - outflow_damping_width ) ) / &1572 REAL( outflow_damping_width/2 ) &1573 )1574 ELSE1575 km_damp_y(j) = 0.01576 ENDIF1577 ENDDO1578 1579 ELSEIF ( bc_ns_raddir ) THEN1580 1581 DO j = nysg, nyng1582 IF ( j <= outflow_damping_width ) THEN1583 km_damp_y(j) = km_damp_max * MIN( 1.0, &1584 ( outflow_damping_width - j ) / &1585 REAL( outflow_damping_width/2 ) &1586 )1587 ELSE1588 km_damp_y(j) = 0.01589 ENDIF1590 ENDDO1591 1592 1585 ENDIF 1593 1586 -
palm/trunk/SOURCE/init_grid.f90
r928 r978 4 4 ! Current revisions: 5 5 ! ----------------- 6 ! 6 ! Bugfix: nzb_max is set to nzt at non-cyclic lateral boundaries 7 ! Bugfix: Set wall_flags_0 for inflow boundary 7 8 ! 8 9 ! Former revisions: … … 608 609 !-- Determine the maximum level of topography. Furthermore it is used for 609 610 !-- steering the degradation of order of the applied advection scheme. 611 !-- In case of non-cyclic lateral boundaries, the order of the advection 612 !-- scheme is reduced at the lateral boundaries up to nzt. 610 613 nzb_max = MAXVAL( nzb_local ) 614 IF ( inflow_l .OR. outflow_l .OR. inflow_r .OR. outflow_r .OR. & 615 inflow_n .OR. outflow_n .OR. inflow_s .OR. outflow_s ) THEN 616 nzb_max = nzt 617 ENDIF 618 611 619 ! 612 620 !-- Consistency checks and index array initialization are only required for … … 1079 1087 !-- scalar - x-direction 1080 1088 !-- WS1 (0), WS3 (1), WS5 (2) 1081 IF ( k <= nzb_s_inner(j,i+1) .OR. ( outflow_l .AND. i == nxl ) & 1082 .OR. ( outflow_r .AND. i == nxr ) ) THEN 1089 IF ( k <= nzb_s_inner(j,i+1) .OR. ( ( inflow_l .OR. outflow_l )& 1090 .AND. i == nxl ) .OR. ( ( inflow_r .OR. outflow_r ) & 1091 .AND. i == nxr ) ) THEN 1083 1092 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 0 ) 1084 1093 ELSEIF ( k <= nzb_s_inner(j,i+2) .OR. k <= nzb_s_inner(j,i-1) & 1085 .OR. ( outflow_r .AND. i == nxr-1 ) & 1086 .OR. ( outflow_l .AND. i == nxlu ) ) THEN 1094 .OR. ( ( inflow_r .OR. outflow_r ) .AND. i == nxr-1 ) & 1095 .OR. ( ( inflow_l .OR. outflow_l ) .AND. i == nxlu ) & 1096 ) THEN 1087 1097 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 1 ) 1088 1098 ELSE … … 1092 1102 !-- scalar - y-direction 1093 1103 !-- WS1 (3), WS3 (4), WS5 (5) 1094 IF ( k <= nzb_s_inner(j+1,i) .OR. ( outflow_s .AND. j == nys ) & 1095 .OR. ( outflow_n .AND. j == nyn ) ) THEN 1104 IF ( k <= nzb_s_inner(j+1,i) .OR. ( ( inflow_s .OR. outflow_s )& 1105 .AND. j == nys ) .OR. ( ( inflow_n .OR. outflow_n ) & 1106 .AND. j == nyn ) ) THEN 1096 1107 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 3 ) 1097 1108 !-- WS3 1098 1109 ELSEIF ( k <= nzb_s_inner(j+2,i) .OR. k <= nzb_s_inner(j-1,i) & 1099 .OR. ( outflow_s .AND. j == nysv ) & 1100 .OR. ( outflow_n .AND. j == nyn-1 ) ) THEN 1110 .OR. ( ( inflow_s .OR. outflow_s ) .AND. j == nysv ) & 1111 .OR. ( ( inflow_n .OR. outflow_n ) .AND. j == nyn-1 ) & 1112 ) THEN 1101 1113 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 4 ) 1102 1114 !-- WS5 … … 1135 1147 !-- WS1 (9), WS3 (10), WS5 (11) 1136 1148 IF ( k <= nzb_u_inner(j,i+1) & 1137 .OR. ( outflow_r .AND. i == nxr ) ) THEN 1149 .OR. ( ( inflow_l .OR. outflow_l ) .AND. i == nxlu ) & 1150 .OR. ( ( inflow_r .OR. outflow_r ) .AND. i == nxr ) & 1151 ) THEN 1138 1152 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 9 ) 1139 1153 ELSEIF ( k <= nzb_u_inner(j,i+2) .OR. k <= nzb_u_inner(j,i-1) & 1140 .OR. ( outflow_r .AND. i == nxr-1 ) & 1141 .OR. ( outflow_l .AND. i == nxlu ) ) THEN 1154 .OR. ( ( inflow_r .OR. outflow_r ) .AND. i == nxr-1 )& 1155 .OR. ( ( inflow_l .OR. outflow_l ) .AND. i == nxlu+1)& 1156 ) THEN 1142 1157 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 10 ) 1143 1158 ELSE 1144 1159 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 11 ) 1145 1160 ENDIF 1161 1146 1162 ! 1147 1163 !-- u component - y-direction 1148 1164 !-- WS1 (12), WS3 (13), WS5 (14) 1149 IF ( k <= nzb_u_inner(j+1,i) .OR. ( outflow_s .AND. j == nys ) & 1150 .OR. ( outflow_n .AND. j == nyn ) ) THEN 1165 IF ( k <= nzb_u_inner(j+1,i) .OR. ( ( inflow_s .OR. outflow_s )& 1166 .AND. j == nys ) .OR. ( ( inflow_n .OR. outflow_n ) & 1167 .AND. j == nyn ) ) THEN 1151 1168 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 12 ) 1152 1169 ELSEIF ( k <= nzb_u_inner(j+2,i) .OR. k <= nzb_u_inner(j-1,i) & 1153 .OR. ( outflow_s .AND. j == nysv ) & 1154 .OR. ( outflow_n .AND. j == nyn-1 ) ) THEN 1170 .OR. ( ( inflow_s .OR. outflow_s ) .AND. j == nysv ) & 1171 .OR. ( ( inflow_n .OR. outflow_n ) .AND. j == nyn-1 ) & 1172 ) THEN 1155 1173 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 13 ) 1156 1174 ELSE … … 1181 1199 !-- v component - x-direction 1182 1200 !-- WS1 (18), WS3 (19), WS5 (20) 1183 IF ( k <= nzb_v_inner(j,i+1) .OR. ( outflow_l .AND. i == nxl ) & 1184 .OR. ( outflow_r .AND. i == nxr ) ) THEN 1201 IF ( k <= nzb_v_inner(j,i+1) .OR. ( ( inflow_l .OR. outflow_l )& 1202 .AND. i == nxl ) .OR. (( inflow_r .OR. outflow_r ) & 1203 .AND. i == nxr ) ) THEN 1185 1204 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 18 ) 1186 1205 !-- WS3 1187 1206 ELSEIF ( k <= nzb_v_inner(j,i+2) .OR. k <= nzb_v_inner(j,i-1) & 1188 .OR. ( outflow_r .AND. i == nxr-1 ) & 1189 .OR. ( outflow_l .AND. i == nxlu ) ) THEN 1207 .OR. ( ( inflow_r .OR. outflow_r ) .AND. i == nxr-1 ) & 1208 .OR. ( ( inflow_l .OR. outflow_l ) .AND. i == nxlu ) & 1209 ) THEN 1190 1210 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 19 ) 1191 1211 ELSE … … 1196 1216 !-- WS1 (21), WS3 (22), WS5 (23) 1197 1217 IF ( k <= nzb_v_inner(j+1,i) & 1198 .OR. ( outflow_n .AND. j == nyn ) ) THEN 1218 .OR. ( ( inflow_s .OR. outflow_s ) .AND. i == nysv ) & 1219 .OR. ( ( inflow_n .OR. outflow_n ) .AND. j == nyn ) & 1220 ) THEN 1199 1221 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 21 ) 1200 1222 ELSEIF ( k <= nzb_v_inner(j+2,i) .OR. k <= nzb_v_inner(j-1,i) & 1201 .OR. ( outflow_s .AND. j == nysv ) & 1202 .OR. ( outflow_n .AND. j == nyn-1 ) ) THEN 1223 .OR. ( ( inflow_s .OR. outflow_s ) .AND. j == nysv+1 )& 1224 .OR. ( ( inflow_n .OR. outflow_n ) .AND. j == nyn-1 )& 1225 ) THEN 1203 1226 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 22 ) 1204 1227 ELSE … … 1228 1251 !-- w component - x-direction 1229 1252 !-- WS1 (27), WS3 (28), WS5 (29) 1230 IF ( k <= nzb_w_inner(j,i+1) .OR. ( outflow_l .AND. i == nxl ) & 1231 .OR. ( outflow_r .AND. i == nxr ) ) THEN 1253 IF ( k <= nzb_w_inner(j,i+1) .OR. ( ( inflow_l .OR. outflow_l )& 1254 .AND. i == nxl ) .OR. ( ( inflow_r .OR. outflow_r ) & 1255 .AND. i == nxr ) ) THEN 1232 1256 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 27 ) 1233 1257 ELSEIF ( k <= nzb_w_inner(j,i+2) .OR. k <= nzb_w_inner(j,i-1) & 1234 .OR. ( outflow_r .AND. i == nxr-1 ) & 1235 .OR. ( outflow_l .AND. i == nxlu ) ) THEN 1258 .OR. ( ( inflow_r .OR. outflow_r ) .AND. i == nxr-1 ) & 1259 .OR. ( ( inflow_l .OR. outflow_l ) .AND. i == nxlu ) & 1260 ) THEN 1236 1261 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 28 ) 1237 1262 ELSE … … 1241 1266 !-- w component - y-direction 1242 1267 !-- WS1 (30), WS3 (31), WS5 (32) 1243 IF ( k <= nzb_w_inner(j+1,i) .OR. ( outflow_s .AND. j == nys ) & 1244 .OR. ( outflow_n .AND. j == nyn ) ) THEN 1268 IF ( k <= nzb_w_inner(j+1,i) .OR. ( ( inflow_s .OR. outflow_s )& 1269 .AND. j == nys ) .OR. ( ( inflow_n .OR. outflow_n ) & 1270 .AND. j == nyn ) ) THEN 1245 1271 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 30 ) 1246 1272 ELSEIF ( k <= nzb_w_inner(j+2,i) .OR. k <= nzb_w_inner(j-1,i) & 1247 .OR. ( outflow_s .AND. j == nysv ) & 1248 .OR. ( outflow_n .AND. j == nyn-1 ) ) THEN 1273 .OR. ( ( inflow_s .OR. outflow_s ) .AND. j == nysv ) & 1274 .OR. ( ( inflow_n .OR. outflow_n ) .AND. j == nyn-1 ) & 1275 ) THEN 1249 1276 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 31 ) 1250 1277 ELSE -
palm/trunk/SOURCE/init_masks.f90
r810 r978 4 4 ! Current revisions: 5 5 ! ----------------- 6 ! 6 ! +z0h* 7 7 ! 8 8 ! Former revisions: … … 257 257 unit = 'psu' 258 258 259 CASE ( 'u*', 't*', 'lwp*', 'pra*', 'prr*', 'z0*' )259 CASE ( 'u*', 't*', 'lwp*', 'pra*', 'prr*', 'z0*', 'z0h*' ) 260 260 WRITE ( message_string, * ) 'illegal value for data_',& 261 261 'output: "', TRIM( var ), '" is only allowed', & -
palm/trunk/SOURCE/init_pegrid.f90
r810 r978 4 4 ! Current revisions: 5 5 ! ----------------- 6 ! dirichlet/neumann and neumann/dirichlet added 7 ! nxlu and nysv are also calculated for inflow boundary 6 8 ! 7 9 ! ATTENTION: nnz_x undefined problem still has to be solved!!!!!!!! … … 1108 1110 !-- horizontal boundary conditions. 1109 1111 IF ( pleft == MPI_PROC_NULL ) THEN 1110 IF ( bc_lr == 'dirichlet/radiation' ) THEN1112 IF ( bc_lr == 'dirichlet/radiation' .OR. bc_lr == 'dirichlet/neumann' ) THEN 1111 1113 inflow_l = .TRUE. 1112 ELSEIF ( bc_lr == 'radiation/dirichlet' ) THEN1114 ELSEIF ( bc_lr == 'radiation/dirichlet' .OR. bc_lr == 'neumann/dirichlet' ) THEN 1113 1115 outflow_l = .TRUE. 1114 1116 ENDIF … … 1116 1118 1117 1119 IF ( pright == MPI_PROC_NULL ) THEN 1118 IF ( bc_lr == 'dirichlet/radiation' ) THEN1120 IF ( bc_lr == 'dirichlet/radiation' .OR. bc_lr == 'dirichlet/neumann' ) THEN 1119 1121 outflow_r = .TRUE. 1120 ELSEIF ( bc_lr == 'radiation/dirichlet' ) THEN1122 ELSEIF ( bc_lr == 'radiation/dirichlet' .OR. bc_lr == 'neumann/dirichlet' ) THEN 1121 1123 inflow_r = .TRUE. 1122 1124 ENDIF … … 1124 1126 1125 1127 IF ( psouth == MPI_PROC_NULL ) THEN 1126 IF ( bc_ns == 'dirichlet/radiation' ) THEN1128 IF ( bc_ns == 'dirichlet/radiation' .OR. bc_ns == 'dirichlet/neumann' ) THEN 1127 1129 outflow_s = .TRUE. 1128 ELSEIF ( bc_ns == 'radiation/dirichlet' ) THEN1130 ELSEIF ( bc_ns == 'radiation/dirichlet' .OR. bc_ns == 'neumann/dirichlet' ) THEN 1129 1131 inflow_s = .TRUE. 1130 1132 ENDIF … … 1132 1134 1133 1135 IF ( pnorth == MPI_PROC_NULL ) THEN 1134 IF ( bc_ns == 'dirichlet/radiation' ) THEN1136 IF ( bc_ns == 'dirichlet/radiation' .OR. bc_ns == 'dirichlet/neumann' ) THEN 1135 1137 inflow_n = .TRUE. 1136 ELSEIF ( bc_ns == 'radiation/dirichlet' ) THEN1138 ELSEIF ( bc_ns == 'radiation/dirichlet' .OR. bc_ns == 'neumann/dirichlet' ) THEN 1137 1139 outflow_n = .TRUE. 1138 1140 ENDIF … … 1164 1166 1165 1167 #elif ! defined ( __parallel ) 1166 IF ( bc_lr == 'dirichlet/radiation' ) THEN1168 IF ( bc_lr == 'dirichlet/radiation' .OR. bc_lr == 'dirichlet/neumann' ) THEN 1167 1169 inflow_l = .TRUE. 1168 1170 outflow_r = .TRUE. 1169 ELSEIF ( bc_lr == 'radiation/dirichlet' ) THEN1171 ELSEIF ( bc_lr == 'radiation/dirichlet' .OR. bc_lr == 'neumann/dirichlet' ) THEN 1170 1172 outflow_l = .TRUE. 1171 1173 inflow_r = .TRUE. 1172 1174 ENDIF 1173 1175 1174 IF ( bc_ns == 'dirichlet/radiation' ) THEN1176 IF ( bc_ns == 'dirichlet/radiation' .OR. bc_ns == 'dirichlet/neumann' ) THEN 1175 1177 inflow_n = .TRUE. 1176 1178 outflow_s = .TRUE. 1177 ELSEIF ( bc_ns == 'radiation/dirichlet' ) THEN1179 ELSEIF ( bc_ns == 'radiation/dirichlet' .OR. bc_ns == 'neumann/dirichlet' ) THEN 1178 1180 outflow_n = .TRUE. 1179 1181 inflow_s = .TRUE. … … 1182 1184 1183 1185 ! 1184 !-- At the outflow, u or v, respectively, have to be calculated for one more1185 !-- grid point.1186 IF ( outflow_l ) THEN1186 !-- At the inflow or outflow, u or v, respectively, have to be calculated for 1187 !-- one more grid point. 1188 IF ( inflow_l .OR. outflow_l ) THEN 1187 1189 nxlu = nxl + 1 1188 1190 ELSE 1189 1191 nxlu = nxl 1190 1192 ENDIF 1191 IF ( outflow_s ) THEN1193 IF ( inflow_s .OR. outflow_s ) THEN 1192 1194 nysv = nys + 1 1193 1195 ELSE -
palm/trunk/SOURCE/modules.f90
r965 r978 4 4 ! Current revisions: 5 5 ! ----------------- 6 ! 6 ! +c_u_m, c_u_m_l, c_v_m, c_v_m_l, c_w_m, c_w_m_l, 7 ! +bc_lr_dirneu, bc_lr_neudir, bc_ns_dirneu, bc_ns_neudir 8 ! -km_damp_x, km_damp_y, km_damp_max, outflow_damping_width 9 ! +z0h, z0h_av, z0h_factor, z0h1d 10 ! +ptdf_x, ptdf_y, pt_damping_width, pt_damping_factor 7 11 ! 8 12 ! Former revisions: … … 324 328 325 329 REAL, DIMENSION(:), ALLOCATABLE :: & 326 ddzu, ddzu_pres, dd2zu, dzu, ddzw, dzw, hyp, inflow_damping_factor,&327 km_damp_x, km_damp_y, lad, l_grid, pt_init, q_init, rdf, rdf_sc,&328 sa_init, ug, u_init, u_nzb_p1_for_vfc, vg, v_init, v_nzb_p1_for_vfc,&329 w_subs, zu, zw330 c_u_m, c_u_m_l, c_v_m, c_v_m_l, c_w_m, c_w_m_l, ddzu, ddzu_pres, & 331 dd2zu, dzu, ddzw, dzw, hyp, inflow_damping_factor, lad, l_grid, & 332 ptdf_x, ptdf_y, pt_init, q_init, rdf, rdf_sc, sa_init, ug, u_init, & 333 u_nzb_p1_for_vfc, vg, v_init, v_nzb_p1_for_vfc, w_subs, zu, zw 330 334 331 335 REAL, DIMENSION(:,:), ALLOCATABLE :: & 332 c_u, c_v, c_w, d zu_mg, dzw_mg, f1_mg, f2_mg, f3_mg,&333 mean_inflow_profiles, pt_slope_ref, qs, qswst_remote, ts, us, z0,&334 flux_s_ pt, diss_s_pt, flux_s_e, diss_s_e, flux_s_q, diss_s_q, &335 flux_s_sa, diss_s_sa, flux_s_u, flux_s_v, flux_s_w, diss_s_u,&336 diss_s_v, diss_s_w, total_2d_o, total_2d_a336 c_u, c_v, c_w, diss_s_e, diss_s_pt, diss_s_q, diss_s_sa, diss_s_u, & 337 diss_s_v, diss_s_w, dzu_mg, dzw_mg, flux_s_e, flux_s_pt, flux_s_q, & 338 flux_s_sa, flux_s_u, flux_s_v, flux_s_w, f1_mg, f2_mg, f3_mg, & 339 mean_inflow_profiles, pt_slope_ref, qs, qswst_remote, total_2d_a, & 340 total_2d_o, ts, us, z0, z0h 337 341 338 342 REAL, DIMENSION(:,:), ALLOCATABLE, TARGET :: & … … 347 351 348 352 REAL, DIMENSION(:,:,:), ALLOCATABLE :: & 349 canopy_heat_flux, cdc, d, de_dx, de_dy, de_dz, diss, lad_s, lad_u,&350 lad_v, lad_w, lai, l_wall, p_loc, sec, sls, tend, u_m_l, u_m_n,&351 u_m_r, u_m_s, v_m_l, v_m_n, v_m_r, v_m_s, w_m_l, w_m_n, w_m_r,&352 w_m_s, flux_l_pt, diss_l_pt, flux_l_e, diss_l_e, flux_l_q, diss_l_q,&353 flux_l_sa, diss_l_sa, flux_l_u, flux_l_v, flux_l_w, diss_l_u,&354 diss_l_v, diss_l_w353 canopy_heat_flux, cdc, d, de_dx, de_dy, de_dz, diss, diss_l_e, & 354 diss_l_pt, diss_l_q, diss_l_sa, diss_l_u, diss_l_v, diss_l_w, & 355 flux_l_e, flux_l_pt, flux_l_q, flux_l_sa, flux_l_u, flux_l_v, & 356 flux_l_w, lad_s, lad_u, lad_v, lad_w, lai, l_wall, p_loc, sec, sls, & 357 tend, u_m_l, u_m_n, u_m_r, u_m_s, v_m_l, v_m_n, v_m_r, v_m_s, w_m_l, & 358 w_m_n, w_m_r, w_m_s 355 359 356 360 REAL, DIMENSION(:,:,:), ALLOCATABLE, TARGET :: & … … 388 392 389 393 REAL, DIMENSION(:,:), ALLOCATABLE :: lwp_av, precipitation_rate_av, & 390 qsws_av, shf_av,ts_av, us_av, z0_av 394 qsws_av, shf_av,ts_av, us_av, z0_av, & 395 z0h_av 391 396 392 397 REAL, DIMENSION(:,:,:), ALLOCATABLE, TARGET :: & … … 549 554 netcdf_data_format = 2, ngsrb = 2, nsor = 20, & 550 555 nsor_ini = 100, n_sor, normalizing_region = 0, & 551 nz_do3d = -9999, outflow_damping_width = -1, &552 pch_index = 0, prt_time_count = 0,recycling_plane, runnr = 0, &556 nz_do3d = -9999, pch_index = 0, prt_time_count = 0, & 557 recycling_plane, runnr = 0, & 553 558 skip_do_avs = 0, subdomain_size, terminate_coupled = 0, & 554 559 terminate_coupled_remote = 0, timestep_count = 0 … … 582 587 583 588 LOGICAL :: adjust_mixing_length = .FALSE., avs_output = .FALSE., & 584 bc_lr_cyc =.TRUE., bc_lr_dirrad = .FALSE., & 589 bc_lr_cyc =.TRUE., bc_lr_dirneu = .FALSE., & 590 bc_lr_dirrad = .FALSE., bc_lr_neudir = .FALSE., & 585 591 bc_lr_raddir = .FALSE., bc_ns_cyc = .TRUE., & 586 bc_ns_dirrad = .FALSE., bc_ns_raddir = .FALSE., & 592 bc_ns_dirneu = .FALSE., bc_ns_dirrad = .FALSE., & 593 bc_ns_neudir = .FALSE., bc_ns_raddir = .FALSE., & 587 594 call_psolver_at_all_substeps = .TRUE., & 588 595 cloud_droplets = .FALSE., cloud_physics = .FALSE., & … … 659 666 f = 0.0, fs = 0.0, g = 9.81, inflow_damping_height = 9999999.9, & 660 667 inflow_damping_width = 9999999.9, kappa = 0.4, km_constant = -1.0,& 661 km_damp_max = -1.0,lad_surface = 0.0, &668 lad_surface = 0.0, & 662 669 leaf_surface_concentration = 0.0, long_filter_factor = 0.0, & 663 670 mask_scale_x = 1.0, mask_scale_y = 1.0, mask_scale_z = 1.0, & … … 670 677 phi = 55.0, prandtl_number = 1.0, & 671 678 precipitation_amount_interval = 9999999.9, prho_reference, & 679 pt_damping_factor = 0.0, pt_damping_width = 0.0, & 672 680 pt_reference = 9999999.9, pt_slope_offset = 0.0, & 673 681 pt_surface = 300.0, pt_surface_initial_change = 0.0, & … … 700 708 ups_limit_v = 0.0, ups_limit_w = 0.0, vg_surface = 0.0, & 701 709 v_bulk = 0.0, v_gtrans = 0.0, wall_adjustment_factor = 1.8, & 702 z_max_do2d = -1.0 710 z_max_do2d = -1.0, z0h_factor = 1.0 703 711 704 712 REAL :: do2d_xy_last_time(0:1) = -1.0, do2d_xz_last_time(0:1) = -1.0, & … … 1057 1065 qs1d, simulated_time_1d = 0.0, time_pr_1d = 0.0, & 1058 1066 time_run_control_1d = 0.0, ts1d, us1d, usws1d, usws1d_m, & 1059 vsws1d, vsws1d_m, z01d 1067 vsws1d, vsws1d_m, z01d, z0h1d 1060 1068 1061 1069 -
palm/trunk/SOURCE/parin.f90
r965 r978 4 4 ! Current revisions: 5 5 ! ----------------- 6 ! 6 ! -km_damp_max, outflow_damping_width 7 ! +pt_damping_factor, pt_damping_width 8 ! +z0h_factor 7 9 ! 8 10 ! Former revisions: … … 168 170 inflow_damping_height, inflow_damping_width, & 169 171 inflow_disturbance_begin, inflow_disturbance_end, & 170 initializing_actions, km_constant, km_damp_max,lad_surface, &172 initializing_actions, km_constant, lad_surface, & 171 173 lad_vertical_gradient, lad_vertical_gradient_level, & 172 174 leaf_surface_concentration, long_filter_factor, & … … 175 177 netcdf_precision, neutral, ngsrb, nsor, & 176 178 nsor_ini, nx, ny, nz, ocean, omega, omega_sor, & 177 o utflow_damping_width, overshoot_limit_e, overshoot_limit_pt, &179 overshoot_limit_e, overshoot_limit_pt, & 178 180 overshoot_limit_u, overshoot_limit_v, overshoot_limit_w, & 179 181 passive_scalar, pch_index, phi, plant_canopy, prandtl_layer, & 180 prandtl_number, precipitation, psolver, pt_reference, pt_surface, & 182 prandtl_number, precipitation, psolver, pt_damping_factor, & 183 pt_damping_width, pt_reference, pt_surface, & 181 184 pt_surface_initial_change, pt_vertical_gradient, & 182 185 pt_vertical_gradient_level, q_surface, q_surface_initial_change, & … … 200 203 uv_heights, u_bulk, u_profile, vg_surface, vg_vertical_gradient, & 201 204 vg_vertical_gradient_level, v_bulk, v_profile, wall_adjustment, & 202 wall_heatflux, wall_humidityflux, wall_scalarflux 205 wall_heatflux, wall_humidityflux, wall_scalarflux, z0h_factor 203 206 204 207 -
palm/trunk/SOURCE/poismg.f90
r881 r978 7 7 ! Current revisions: 8 8 ! ----------------- 9 ! 9 ! bc_lr/ns_dirneu/neudir added 10 ! 10 11 ! 11 12 ! Former revisions: … … 1229 1230 !-- outflow conditions have to be used on all PEs after the switch, 1230 1231 !-- because then they have the total domain. 1231 IF ( bc_lr_dirrad ) THEN1232 IF ( bc_lr_dirrad .OR. bc_lr_dirneu ) THEN 1232 1233 inflow_l = .TRUE. 1233 1234 inflow_r = .FALSE. 1234 1235 outflow_l = .FALSE. 1235 1236 outflow_r = .TRUE. 1236 ELSEIF ( bc_lr_raddir ) THEN1237 ELSEIF ( bc_lr_raddir .OR. bc_lr_neudir ) THEN 1237 1238 inflow_l = .FALSE. 1238 1239 inflow_r = .TRUE. … … 1241 1242 ENDIF 1242 1243 1243 IF ( bc_ns_dirrad ) THEN1244 IF ( bc_ns_dirrad .OR. bc_ns_dirneu ) THEN 1244 1245 inflow_n = .TRUE. 1245 1246 inflow_s = .FALSE. 1246 1247 outflow_n = .FALSE. 1247 1248 outflow_s = .TRUE. 1248 ELSEIF ( bc_ns_raddir ) THEN1249 ELSEIF ( bc_ns_raddir .OR. bc_ns_neudir ) THEN 1249 1250 inflow_n = .FALSE. 1250 1251 inflow_s = .TRUE. … … 1316 1317 1317 1318 IF ( pleft == MPI_PROC_NULL ) THEN 1318 IF ( bc_lr_dirrad ) THEN1319 IF ( bc_lr_dirrad .OR. bc_lr_dirneu ) THEN 1319 1320 inflow_l = .TRUE. 1320 ELSEIF ( bc_lr_raddir ) THEN1321 ELSEIF ( bc_lr_raddir .OR. bc_lr_neudir ) THEN 1321 1322 outflow_l = .TRUE. 1322 1323 ENDIF … … 1324 1325 1325 1326 IF ( pright == MPI_PROC_NULL ) THEN 1326 IF ( bc_lr_dirrad ) THEN1327 IF ( bc_lr_dirrad .OR. bc_lr_dirneu ) THEN 1327 1328 outflow_r = .TRUE. 1328 ELSEIF ( bc_lr_raddir ) THEN1329 ELSEIF ( bc_lr_raddir .OR. bc_lr_neudir ) THEN 1329 1330 inflow_r = .TRUE. 1330 1331 ENDIF … … 1332 1333 1333 1334 IF ( psouth == MPI_PROC_NULL ) THEN 1334 IF ( bc_ns_dirrad ) THEN1335 IF ( bc_ns_dirrad .OR. bc_ns_dirneu ) THEN 1335 1336 outflow_s = .TRUE. 1336 ELSEIF ( bc_ns_raddir ) THEN1337 ELSEIF ( bc_ns_raddir .OR. bc_ns_neudir ) THEN 1337 1338 inflow_s = .TRUE. 1338 1339 ENDIF … … 1340 1341 1341 1342 IF ( pnorth == MPI_PROC_NULL ) THEN 1342 IF ( bc_ns_dirrad ) THEN1343 IF ( bc_ns_dirrad .OR. bc_ns_dirneu ) THEN 1343 1344 inflow_n = .TRUE. 1344 ELSEIF ( bc_ns_raddir ) THEN1345 ELSEIF ( bc_ns_raddir .OR. bc_ns_neudir ) THEN 1345 1346 outflow_n = .TRUE. 1346 1347 ENDIF -
palm/trunk/SOURCE/prandtl_fluxes.f90
r760 r978 4 4 ! Current revisions: 5 5 ! ----------------- 6 ! 6 ! roughness length for scalar quantities z0h added 7 7 ! 8 8 ! Former revisions: … … 96 96 ! 97 97 !-- Stable stratification 98 ts(j,i) = kappa * ( pt(k+1,j,i) - pt(k,j,i) ) / ( &99 LOG( z_p / z0 (j,i) ) + &100 5.0 * rif(j,i) * ( z_p - z0 (j,i) ) / z_p &98 ts(j,i) = kappa * ( pt(k+1,j,i) - pt(k,j,i) ) / ( & 99 LOG( z_p / z0h(j,i) ) + & 100 5.0 * rif(j,i) * ( z_p - z0h(j,i) ) / z_p & 101 101 ) 102 102 ELSE … … 104 104 !-- Unstable stratification 105 105 a = SQRT( 1.0 - 16.0 * rif(j,i) ) 106 b = SQRT( 1.0 - 16.0 * rif(j,i) * z0 (j,i) / z_p )107 108 ts(j,i) = kappa * ( pt(k+1,j,i) - pt(k,j,i) ) / ( &109 LOG( z_p / z0 (j,i) ) - &106 b = SQRT( 1.0 - 16.0 * rif(j,i) * z0h(j,i) / z_p ) 107 108 ts(j,i) = kappa * ( pt(k+1,j,i) - pt(k,j,i) ) / ( & 109 LOG( z_p / z0h(j,i) ) - & 110 110 2.0 * LOG( ( 1.0 + a ) / ( 1.0 + b ) ) ) 111 111 ENDIF … … 308 308 ! 309 309 !-- Stable stratification 310 qs(j,i) = kappa * ( q(k+1,j,i) - q(k,j,i) ) / ( &311 LOG( z_p / z0 (j,i) ) + &312 5.0 * rif(j,i) * ( z_p - z0 (j,i) ) / z_p &310 qs(j,i) = kappa * ( q(k+1,j,i) - q(k,j,i) ) / ( & 311 LOG( z_p / z0h(j,i) ) + & 312 5.0 * rif(j,i) * ( z_p - z0h(j,i) ) / z_p & 313 313 ) 314 314 ELSE … … 316 316 !-- Unstable stratification 317 317 a = SQRT( 1.0 - 16.0 * rif(j,i) ) 318 b = SQRT( 1.0 - 16.0 * rif(j,i) * z0 (j,i) / z_p )318 b = SQRT( 1.0 - 16.0 * rif(j,i) * z0h(j,i) / z_p ) 319 319 320 qs(j,i) = kappa * ( q(k+1,j,i) - q(k,j,i) ) / ( &321 LOG( z_p / z0 (j,i) ) -&320 qs(j,i) = kappa * ( q(k+1,j,i) - q(k,j,i) ) / ( & 321 LOG( z_p / z0h(j,i) ) - & 322 322 2.0 * LOG( (1.0 + a ) / ( 1.0 + b ) ) ) 323 323 ENDIF -
palm/trunk/SOURCE/prognostic_equations.f90
r941 r978 4 4 ! Current revisions: 5 5 ! ----------------- 6 ! 6 ! km_damp_x and km_damp_y removed in calls of diffusion_u and diffusion_v 7 ! add ptdf_x, ptdf_y for damping the potential temperature at the inflow 8 ! boundary in case of non-cyclic lateral boundaries 9 ! Bugfix: first thread index changes for WS-scheme at the inflow 7 10 ! 8 11 ! Former revisions: … … 201 204 ENDIF 202 205 IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) THEN 203 CALL diffusion_u( i, j, ddzu, ddzw, km_m, km_damp_y, tend, u_m,&204 usws _m, uswst_m, v_m, w_m )206 CALL diffusion_u( i, j, ddzu, ddzw, km_m, tend, u_m, usws_m, & 207 uswst_m, v_m, w_m ) 205 208 ELSE 206 CALL diffusion_u( i, j, ddzu, ddzw, km, km_damp_y, tend, u, usws, &207 uswst,v, w )209 CALL diffusion_u( i, j, ddzu, ddzw, km, tend, u, usws, uswst, & 210 v, w ) 208 211 ENDIF 209 212 CALL coriolis( i, j, 1 ) … … 289 292 ENDIF 290 293 IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) THEN 291 CALL diffusion_v( i, j, ddzu, ddzw, km_m, km_damp_x, tend, u_m, &292 v _m, vsws_m, vswst_m, w_m )294 CALL diffusion_v( i, j, ddzu, ddzw, km_m, tend, u_m, v_m, & 295 vsws_m, vswst_m, w_m ) 293 296 ELSE 294 CALL diffusion_v( i, j, ddzu, ddzw, km, km_damp_x, tend, u, v,&295 vsws , vswst, w )297 CALL diffusion_v( i, j, ddzu, ddzw, km, tend, u, v, vsws, & 298 vswst, w ) 296 299 ENDIF 297 300 CALL coriolis( i, j, 2 ) … … 373 376 ENDIF 374 377 IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) THEN 375 CALL diffusion_w( i, j, ddzu, ddzw, km_m, km_damp_x, km_damp_y, & 376 tend, u_m, v_m, w_m ) 378 CALL diffusion_w( i, j, ddzu, ddzw, km_m, tend, u_m, v_m, w_m ) 377 379 ELSE 378 CALL diffusion_w( i, j, ddzu, ddzw, km, km_damp_x, km_damp_y, & 379 tend, u, v, w ) 380 CALL diffusion_w( i, j, ddzu, ddzw, km, tend, u, v, w ) 380 381 ENDIF 381 382 CALL coriolis( i, j, 3 ) … … 524 525 DO k = nzb_s_inner(j,i)+1, nzt 525 526 pt_p(k,j,i) = ( 1 - sat ) * pt_m(k,j,i) + sat * pt(k,j,i) + & 526 dt_3d * ( 527 dt_3d * ( & 527 528 sbt * tend(k,j,i) + tsc(3) * tpt_m(k,j,i) & 528 ) - & 529 tsc(5) * rdf_sc(k) * ( pt(k,j,i) - pt_init(k) ) 529 ) - & 530 tsc(5) * ( pt(k,j,i) - pt_init(k) ) * & 531 ( rdf_sc(k) + ptdf_x(i) + ptdf_y(j) ) 530 532 ENDDO 531 533 … … 983 985 IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN 984 986 IF ( ws_scheme_mom ) THEN 985 IF ( outflow_l.AND. i_omp_start == nxl ) THEN987 IF ( ( inflow_l .OR. outflow_l ) .AND. i_omp_start == nxl ) THEN 986 988 ! CALL local_diss( i, j, u) ! dissipation control 987 989 CALL advec_u_ws( i, j, i_omp_start + 1, tn ) … … 997 999 IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) & 998 1000 THEN 999 CALL diffusion_u( i, j, ddzu, ddzw, km_m, km_damp_y, tend,&1000 u _m, usws_m, uswst_m, v_m, w_m )1001 ELSE 1002 CALL diffusion_u( i, j, ddzu, ddzw, km, km_damp_y, tend, u, &1003 usws , uswst, v, w )1001 CALL diffusion_u( i, j, ddzu, ddzw, km_m, tend, u_m, & 1002 usws_m, uswst_m, v_m, w_m ) 1003 ELSE 1004 CALL diffusion_u( i, j, ddzu, ddzw, km, tend, u, usws, & 1005 uswst, v, w ) 1004 1006 ENDIF 1005 1007 CALL coriolis( i, j, 1 ) … … 1066 1068 IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) & 1067 1069 THEN 1068 CALL diffusion_v( i, j, ddzu, ddzw, km_m, km_damp_x, tend,&1069 u_m, v_m,vsws_m, vswst_m, w_m )1070 ELSE 1071 CALL diffusion_v( i, j, ddzu, ddzw, km, km_damp_x, tend, u, v, &1072 vsws , vswst, w )1070 CALL diffusion_v( i, j, ddzu, ddzw, km_m, tend, u_m, v_m, & 1071 vsws_m, vswst_m, w_m ) 1072 ELSE 1073 CALL diffusion_v( i, j, ddzu, ddzw, km, tend, u, v, vsws, & 1074 vswst, w ) 1073 1075 ENDIF 1074 1076 CALL coriolis( i, j, 2 ) … … 1130 1132 IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) & 1131 1133 THEN 1132 CALL diffusion_w( i, j, ddzu, ddzw, km_m, km_damp_x,&1133 km_damp_y, tend, u_m, v_m,w_m )1134 CALL diffusion_w( i, j, ddzu, ddzw, km_m, tend, u_m, v_m, & 1135 w_m ) 1134 1136 ELSE 1135 CALL diffusion_w( i, j, ddzu, ddzw, km, km_damp_x, km_damp_y, & 1136 tend, u, v, w ) 1137 CALL diffusion_w( i, j, ddzu, ddzw, km, tend, u, v, w ) 1137 1138 ENDIF 1138 1139 CALL coriolis( i, j, 3 ) … … 1240 1241 dt_3d * ( & 1241 1242 tsc(2) * tend(k,j,i) + tsc(3) * tpt_m(k,j,i) & 1242 ) - & 1243 tsc(5) * rdf_sc(k) * ( pt(k,j,i) - pt_init(k) ) 1243 ) - & 1244 tsc(5) * ( pt(k,j,i) - pt_init(k) ) * & 1245 ( rdf_sc(k) + ptdf_x(i) + ptdf_y(j) ) 1244 1246 ENDDO 1245 1247 … … 1540 1542 ENDIF 1541 1543 IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) THEN 1542 CALL diffusion_u( ddzu, ddzw, km_m, km_damp_y, tend, u_m, usws_m, &1543 uswst_m,v_m, w_m )1544 CALL diffusion_u( ddzu, ddzw, km_m, tend, u_m, usws_m, uswst_m, & 1545 v_m, w_m ) 1544 1546 ELSE 1545 CALL diffusion_u( ddzu, ddzw, km, km_damp_y,tend, u, usws, uswst, v, w )1547 CALL diffusion_u( ddzu, ddzw, km, tend, u, usws, uswst, v, w ) 1546 1548 ENDIF 1547 1549 CALL coriolis( 1 ) … … 1634 1636 ENDIF 1635 1637 IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) THEN 1636 CALL diffusion_v( ddzu, ddzw, km_m, km_damp_x, tend, u_m, v_m, vsws_m, &1637 vswst_m,w_m )1638 CALL diffusion_v( ddzu, ddzw, km_m, tend, u_m, v_m, vsws_m, vswst_m, & 1639 w_m ) 1638 1640 ELSE 1639 CALL diffusion_v( ddzu, ddzw, km, km_damp_x,tend, u, v, vsws, vswst, w )1641 CALL diffusion_v( ddzu, ddzw, km, tend, u, v, vsws, vswst, w ) 1640 1642 ENDIF 1641 1643 CALL coriolis( 2 ) … … 1725 1727 ENDIF 1726 1728 IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) THEN 1727 CALL diffusion_w( ddzu, ddzw, km_m, km_damp_x, km_damp_y, tend, u_m, & 1728 v_m, w_m ) 1729 CALL diffusion_w( ddzu, ddzw, km_m, tend, u_m, v_m, w_m ) 1729 1730 ELSE 1730 CALL diffusion_w( ddzu, ddzw, km, km_damp_x, km_damp_y,tend, u, v, w )1731 CALL diffusion_w( ddzu, ddzw, km, tend, u, v, w ) 1731 1732 ENDIF 1732 1733 CALL coriolis( 3 ) … … 1877 1878 DO j = nys, nyn 1878 1879 DO k = nzb_s_inner(j,i)+1, nzt 1879 pt_p(k,j,i) = ( 1 - sat ) * pt_m(k,j,i) + sat * pt(k,j,i) + & 1880 dt_3d * ( & 1881 sbt * tend(k,j,i) + tsc(3) * tpt_m(k,j,i) & 1882 ) - & 1883 tsc(5) * rdf_sc(k) * ( pt(k,j,i) - pt_init(k) ) 1880 pt_p(k,j,i) = ( 1 - sat ) * pt_m(k,j,i) + sat * pt(k,j,i) + & 1881 dt_3d * ( & 1882 sbt * tend(k,j,i) + tsc(3) * tpt_m(k,j,i)& 1883 ) - & 1884 tsc(5) * ( pt(k,j,i) - pt_init(k) ) * & 1885 ( rdf_sc(k) + ptdf_x(i) + ptdf_y(j) ) 1884 1886 ENDDO 1885 1887 ENDDO -
palm/trunk/SOURCE/read_3d_binary.f90
r777 r978 4 4 ! Current revisions: 5 5 ! ----------------- 6 ! +z0h, z0h_av 6 7 ! 7 8 ! Former revisions: … … 963 964 tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp) 964 965 966 CASE ( 'z0h' ) 967 IF ( k == 1 ) READ ( 13 ) tmp_2d 968 z0h(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = & 969 tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp) 970 971 CASE ( 'z0h_av' ) 972 IF ( .NOT. ALLOCATED( z0h_av ) ) THEN 973 ALLOCATE( z0h_av(nysg:nyng,nxlg:nxrg) ) 974 ENDIF 975 IF ( k == 1 ) READ ( 13 ) tmp_2d 976 z0h_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = & 977 tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp) 978 965 979 CASE DEFAULT 966 980 WRITE( message_string, * ) 'unknown field named "', & -
palm/trunk/SOURCE/read_var_list.f90
r941 r978 4 4 ! Current revisions: 5 5 ! ------------------ 6 ! 6 ! -km_damp_max, outflow_damping_width 7 ! +pt_damping_factor, pt_damping_width 8 ! +z0h_factor 7 9 ! 8 10 ! Former revisions: … … 385 387 CASE ( 'km_constant' ) 386 388 READ ( 13 ) km_constant 387 CASE ( 'km_damp_max' )388 READ ( 13 ) km_damp_max389 389 CASE ( 'lad' ) 390 390 READ ( 13 ) lad … … 446 446 CASE ( 'omega_sor' ) 447 447 READ ( 13 ) omega_sor 448 CASE ( 'outflow_damping_width' )449 READ ( 13 ) outflow_damping_width450 448 CASE ( 'output_for_t0' ) 451 449 READ (13) output_for_t0 … … 476 474 CASE ( 'psolver' ) 477 475 READ ( 13 ) psolver 476 CASE ( 'pt_damping_factor' ) 477 READ ( 13 ) pt_damping_factor 478 CASE ( 'pt_damping_width' ) 479 READ ( 13 ) pt_damping_width 478 480 CASE ( 'pt_init' ) 479 481 READ ( 13 ) pt_init … … 684 686 CASE ( 'w_max_ijk' ) 685 687 READ ( 13 ) w_max_ijk 688 CASE ( 'z0h_factor' ) 689 READ ( 13 ) z0h_factor 686 690 687 691 CASE DEFAULT -
palm/trunk/SOURCE/sum_up_3d_data.f90
r791 r978 4 4 ! Current revisions: 5 5 ! ----------------- 6 ! 6 ! +z0h* 7 7 ! 8 8 ! Former revisions: … … 227 227 ENDIF 228 228 z0_av = 0.0 229 230 CASE ( 'z0h*' ) 231 IF ( .NOT. ALLOCATED( z0h_av ) ) THEN 232 ALLOCATE( z0h_av(nysg:nyng,nxlg:nxrg) ) 233 ENDIF 234 z0h_av = 0.0 229 235 230 236 CASE DEFAULT … … 494 500 ENDDO 495 501 502 CASE ( 'z0h*' ) 503 DO i = nxlg, nxrg 504 DO j = nysg, nyng 505 z0h_av(j,i) = z0h_av(j,i) + z0h(j,i) 506 ENDDO 507 ENDDO 508 496 509 CASE DEFAULT 497 510 ! -
palm/trunk/SOURCE/timestep.f90
r867 r978 4 4 ! Current revisions: 5 5 ! ------------------ 6 ! 6 ! restriction of the outflow damping layer in the diffusion criterion removed 7 7 ! 8 8 ! Former revisions: … … 185 185 dt_diff = dt_diff_l 186 186 #endif 187 188 !189 !-- In case of non-cyclic lateral boundaries, the diffusion time step190 !-- may be further restricted by the lateral damping layer (damping only191 !-- along x and y)192 IF ( .NOT. bc_lr_cyc ) THEN193 dt_diff = MIN( dt_diff, 0.125 * dx2 / ( km_damp_max + 1E-20 ) )194 ELSEIF ( .NOT. bc_ns_cyc ) THEN195 dt_diff = MIN( dt_diff, 0.125 * dy2 / ( km_damp_max + 1E-20 ) )196 ENDIF197 187 198 188 ! -
palm/trunk/SOURCE/write_3d_binary.f90
r791 r978 4 4 ! Current revisions: 5 5 ! ----------------- 6 ! 6 ! +z0h, z0h_av 7 7 ! 8 8 ! Former revisions: … … 298 298 WRITE ( 14 ) 'z0_av '; WRITE ( 14 ) z0_av 299 299 ENDIF 300 WRITE ( 14 ) 'z0h '; WRITE ( 14 ) z0h 301 IF ( ALLOCATED( z0h_av ) ) THEN 302 WRITE ( 14 ) 'z0h_av '; WRITE ( 14 ) z0h_av 303 ENDIF 300 304 301 305 ! -
palm/trunk/SOURCE/write_var_list.f90
r941 r978 4 4 ! Current revisions: 5 5 ! ----------------- 6 ! 6 ! -km_damp_max, outflow_damping_width 7 ! +pt_damping_factor, pt_damping_width 8 ! +z0h_factor 7 9 ! 8 10 ! Former revisions: … … 300 302 WRITE ( 14 ) 'km_constant ' 301 303 WRITE ( 14 ) km_constant 302 WRITE ( 14 ) 'km_damp_max '303 WRITE ( 14 ) km_damp_max304 304 WRITE ( 14 ) 'lad ' 305 305 WRITE ( 14 ) lad … … 358 358 WRITE ( 14 ) 'omega_sor ' 359 359 WRITE ( 14 ) omega_sor 360 WRITE ( 14 ) 'outflow_damping_width '361 WRITE ( 14 ) outflow_damping_width362 360 WRITE ( 14 ) 'output_for_t0 ' 363 361 WRITE ( 14 ) output_for_t0 … … 388 386 WRITE ( 14 ) 'psolver ' 389 387 WRITE ( 14 ) psolver 388 WRITE ( 14 ) 'pt_damping_factor ' 389 WRITE ( 14 ) pt_damping_factor 390 WRITE ( 14 ) 'pt_damping_width ' 391 WRITE ( 14 ) pt_damping_width 390 392 WRITE ( 14 ) 'pt_init ' 391 393 WRITE ( 14 ) pt_init … … 596 598 WRITE ( 14 ) 'w_max_ijk ' 597 599 WRITE ( 14 ) w_max_ijk 600 WRITE ( 14 ) 'z0h_factor ' 601 WRITE ( 14 ) z0h_factor 598 602 599 603 !
Note: See TracChangeset
for help on using the changeset viewer.
