Changeset 106
- Timestamp:
- Aug 16, 2007 2:30:26 PM (17 years ago)
- Location:
- palm/trunk/SOURCE
- Files:
-
- 22 edited
-
CURRENT_MODIFICATIONS (modified) (2 diffs)
-
advec_particles.f90 (modified) (4 diffs)
-
advec_u_pw.f90 (modified) (2 diffs)
-
advec_u_up.f90 (modified) (2 diffs)
-
advec_v_pw.f90 (modified) (2 diffs)
-
advec_v_up.f90 (modified) (2 diffs)
-
boundary_conds.f90 (modified) (12 diffs)
-
buoyancy.f90 (modified) (2 diffs)
-
check_parameters.f90 (modified) (3 diffs)
-
coriolis.f90 (modified) (3 diffs)
-
data_output_dvrp.f90 (modified) (3 diffs)
-
diffusion_u.f90 (modified) (2 diffs)
-
diffusion_v.f90 (modified) (2 diffs)
-
flow_statistics.f90 (modified) (7 diffs)
-
init_3d_model.f90 (modified) (3 diffs)
-
init_particles.f90 (modified) (4 diffs)
-
init_pegrid.f90 (modified) (3 diffs)
-
modules.f90 (modified) (3 diffs)
-
pres.f90 (modified) (3 diffs)
-
production_e.f90 (modified) (3 diffs)
-
prognostic_equations.f90 (modified) (13 diffs)
-
time_integration.f90 (modified) (3 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/CURRENT_MODIFICATIONS
r105 r106 13 13 inipar parameter top_momentumflux_u|v. 14 14 15 boundary_conds, check_open, check_parameters, diffusion_u, diffusion_v, flow_statistics, header, init_pegrid, init_3d_model, modules, palm, parin, prognostic_equations, read_var_list, read_3d_binary, swap_timelevel, time_integration, write_var_list, write_3d_binary 15 Non-cyclic boundary conditions can be used along all horizontal directions. 16 17 Quantities w*p* and w"e can be output as vertical profiles. 18 19 boundary_conds, check_open, check_parameters, diffusion_u, diffusion_v, flow_statistics, header, init_pegrid, init_3d_model, modules, palm, parin, pres, prognostic_equations, read_var_list, read_3d_binary, swap_timelevel, time_integration, write_var_list, write_3d_binary 16 20 17 21 New: … … 21 25 Changed: 22 26 ------- 27 Remaining variables iran changed to iran_part (advec_particles, init_particles). 28 29 In case that the presure solver is not called for every Runge-Kutta substep 30 (call_psolver_at_all_substeps = .F.), it is called after the first substep 31 instead of the last. In that case, random perturbations are also added to the 32 velocity field after the first substep. 33 34 advec_particles, init_particles, time_integration 23 35 24 36 25 37 Errors: 26 38 ------ 27 +dots_num_palm in module user, +module netcdf_control in user_init (both in user_interface.f90) 39 Bugs from code parts for non-cyclic boundary conditions are removed: loops for 40 u and v are starting from index nxlu, nysv, respectively. The radiation boundary 41 condition is used for every Runge-Kutta substep. Velocity phase speeds for 42 the radiation boundary conditions are calculated for the first Runge-Kutta 43 substep only and reused for the further substeps. New arrays c_u, c_v, and c_w 44 are defined for this purpose. Several index errors are removed from the 45 radiation boundary condition code parts. Upper bounds for calculating 46 u_0 and v_0 (in production_e) are nxr+1 and nyn+1 because otherwise these 47 values are not available in case of non-cyclic boundary conditions. 48 49 +dots_num_palm in module user, +module netcdf_control in user_init (both in user_interface) 50 51 Bugfix: wrong sign removed from the buoyancy production term in the case use_reference = .T. (production_e) 28 52 29 53 Bugfix: Error message concerning output of particle concentration (pc) modified (check_parameters). 30 54 31 check_parameters, user_interface55 advec_u_pw, advec_u_up, advec_v_pw, advec_v_up, boundary_conds, buoyancy, check_parameters, coriolis, diffusion_u, diffusion_v, init_pegrid, init_3d_model, modules, production_e, prognostic_equations, user_interface -
palm/trunk/SOURCE/advec_particles.f90
r98 r106 4 4 ! Actual revisions: 5 5 ! ----------------- 6 ! 6 ! remaining variables iran changed to iran_part 7 7 ! TEST: PRINT statements on unit 9 (commented out) 8 8 ! … … 1955 1955 THEN 1956 1956 particles(n)%x = particles(n)%x + & 1957 ( random_function( iran ) - 0.5 ) *&1957 ( random_function( iran_part ) - 0.5 ) *& 1958 1958 pdx(particles(n)%group) 1959 1959 IF ( particles(n)%x <= ( nxl - 0.5 ) * dx ) THEN … … 1966 1966 THEN 1967 1967 particles(n)%y = particles(n)%y + & 1968 ( random_function( iran ) - 0.5 ) *&1968 ( random_function( iran_part ) - 0.5 ) *& 1969 1969 pdy(particles(n)%group) 1970 1970 IF ( particles(n)%y <= ( nys - 0.5 ) * dy ) THEN … … 1977 1977 THEN 1978 1978 particles(n)%z = particles(n)%z + & 1979 ( random_function( iran ) - 0.5 ) *&1979 ( random_function( iran_part ) - 0.5 ) *& 1980 1980 pdz(particles(n)%group) 1981 1981 ENDIF -
palm/trunk/SOURCE/advec_u_pw.f90
r77 r106 4 4 ! Actual revisions: 5 5 ! ----------------- 6 ! 6 ! i loop is starting from nxlu (needed for non-cyclic boundary conditions) 7 7 ! 8 8 ! Former revisions: … … 58 58 gu = 2.0 * u_gtrans 59 59 gv = 2.0 * v_gtrans 60 DO i = nxl , nxr60 DO i = nxlu, nxr 61 61 DO j = nys, nyn 62 62 DO k = nzb_u_inner(j,i)+1, nzt -
palm/trunk/SOURCE/advec_u_up.f90
r77 r106 4 4 ! Actual revisions: 5 5 ! ----------------- 6 ! 6 ! i loop is starting from nxlu (needed for non-cyclic boundary conditions) 7 7 ! 8 8 ! Former revisions: … … 57 57 58 58 59 DO i = nxl , nxr59 DO i = nxlu, nxr 60 60 DO j = nys, nyn 61 61 DO k = nzb_u_inner(j,i)+1, nzt -
palm/trunk/SOURCE/advec_v_pw.f90
r77 r106 4 4 ! Actual revisions: 5 5 ! ----------------- 6 ! 6 ! j loop is starting from nysv (needed for non-cyclic boundary conditions) 7 7 ! 8 8 ! Former revisions: … … 60 60 gv = 2.0 * v_gtrans 61 61 DO i = nxl, nxr 62 DO j = nys , nyn62 DO j = nysv, nyn 63 63 DO k = nzb_v_inner(j,i)+1, nzt 64 64 tend(k,j,i) = tend(k,j,i) - 0.25 * ( & -
palm/trunk/SOURCE/advec_v_up.f90
r77 r106 4 4 ! Actual revisions: 5 5 ! ----------------- 6 ! 6 ! j loop is starting from nysv (needed for non-cyclic boundary conditions) 7 7 ! 8 8 ! Former revisions: … … 57 57 58 58 DO i = nxl, nxr 59 DO j = nys , nyn59 DO j = nysv, nyn 60 60 DO k = nzb_v_inner(j,i)+1, nzt 61 61 ! -
palm/trunk/SOURCE/boundary_conds.f90
r102 r106 4 4 ! Actual revisions: 5 5 ! ----------------- 6 ! Boundary conditions for temperature adjusted for coupled runs 6 ! Boundary conditions for temperature adjusted for coupled runs, 7 ! bugfixes for the radiation boundary conditions at the outflow: radiation 8 ! conditions are used for every substep, phase speeds are calculated for the 9 ! first Runge-Kutta substep only and then reused, several index values changed 10 ! 7 11 ! 8 12 ! Former revisions: … … 59 63 INTEGER :: i, j, k 60 64 61 REAL :: c_max, c_u, c_v, c_w,denom65 REAL :: c_max, denom 62 66 63 67 … … 216 220 ! 217 221 !-- Radiation boundary condition for the velocities at the respective outflow 218 IF ( outflow_s .AND. & 219 intermediate_timestep_count == intermediate_timestep_count_max ) & 220 THEN 222 IF ( outflow_s ) THEN 221 223 222 224 c_max = dy / dt_3d … … 226 228 227 229 ! 228 !-- First calculate the phase speeds for u,v, and w 229 denom = u_m_s(k,0,i) - u_m_s(k,1,i) 230 231 IF ( denom /= 0.0 ) THEN 232 c_u = -c_max * ( u(k,0,i) - u_m_s(k,0,i) ) / denom 233 IF ( c_u < 0.0 ) THEN 234 c_u = 0.0 235 ELSEIF ( c_u > c_max ) THEN 236 c_u = c_max 237 ENDIF 238 ELSE 239 c_u = c_max 230 !-- Calculate the phase speeds for u,v, and w. In case of using a 231 !-- Runge-Kutta scheme, do this for the first substep only and then 232 !-- reuse this values for the further substeps. 233 IF ( intermediate_timestep_count == 1 ) THEN 234 235 denom = u_m_s(k,0,i) - u_m_s(k,1,i) 236 237 IF ( denom /= 0.0 ) THEN 238 c_u(k,i) = -c_max * ( u(k,0,i) - u_m_s(k,0,i) ) / denom 239 IF ( c_u(k,i) < 0.0 ) THEN 240 c_u(k,i) = 0.0 241 ELSEIF ( c_u(k,i) > c_max ) THEN 242 c_u(k,i) = c_max 243 ENDIF 244 ELSE 245 c_u(k,i) = c_max 246 ENDIF 247 248 denom = v_m_s(k,1,i) - v_m_s(k,2,i) 249 250 IF ( denom /= 0.0 ) THEN 251 c_v(k,i) = -c_max * ( v(k,1,i) - v_m_s(k,1,i) ) / denom 252 IF ( c_v(k,i) < 0.0 ) THEN 253 c_v(k,i) = 0.0 254 ELSEIF ( c_v(k,i) > c_max ) THEN 255 c_v(k,i) = c_max 256 ENDIF 257 ELSE 258 c_v(k,i) = c_max 259 ENDIF 260 261 denom = w_m_s(k,0,i) - w_m_s(k,1,i) 262 263 IF ( denom /= 0.0 ) THEN 264 c_w(k,i) = -c_max * ( w(k,0,i) - w_m_s(k,0,i) ) / denom 265 IF ( c_w(k,i) < 0.0 ) THEN 266 c_w(k,i) = 0.0 267 ELSEIF ( c_w(k,i) > c_max ) THEN 268 c_w(k,i) = c_max 269 ENDIF 270 ELSE 271 c_w(k,i) = c_max 272 ENDIF 273 274 ! 275 !-- Save old timelevels for the next timestep 276 u_m_s(k,:,i) = u(k,0:1,i) 277 v_m_s(k,:,i) = v(k,1:2,i) 278 w_m_s(k,:,i) = w(k,0:1,i) 279 240 280 ENDIF 241 denom = v_m_s(k,0,i) - v_m_s(k,1,i)242 243 IF ( denom /= 0.0 ) THEN244 c_v = -c_max * ( v(k,0,i) - v_m_s(k,0,i) ) / denom245 IF ( c_v < 0.0 ) THEN246 c_v = 0.0247 ELSEIF ( c_v > c_max ) THEN248 c_v = c_max249 ENDIF250 ELSE251 c_v = c_max252 ENDIF253 254 denom = w_m_s(k,0,i) - w_m_s(k,1,i)255 256 IF ( denom /= 0.0 ) THEN257 c_w = -c_max * ( w(k,0,i) - w_m_s(k,0,i) ) / denom258 IF ( c_w < 0.0 ) THEN259 c_w = 0.0260 ELSEIF ( c_w > c_max ) THEN261 c_w = c_max262 ENDIF263 ELSE264 c_w = c_max265 ENDIF266 281 267 282 ! 268 283 !-- Calculate the new velocities 269 u_p(k,-1,i) = u(k,-1,i) + dt_3d * c_u* &284 u_p(k,-1,i) = u(k,-1,i) - dt_3d * tsc(2) * c_u(k,i) * & 270 285 ( u(k,-1,i) - u(k,0,i) ) * ddy 271 286 272 v_p(k,-1,i) = v(k,-1,i) + dt_3d * c_v* &273 ( v(k, -1,i) - v_m_s(k,0,i) ) * ddy274 275 w_p(k,-1,i) = w(k,-1,i) + dt_3d * c_w* &287 v_p(k,-1,i) = v(k,-1,i) - dt_3d * tsc(2) * c_v(k,i) * & 288 ( v(k,0,i) - v(k,1,i) ) * ddy 289 290 w_p(k,-1,i) = w(k,-1,i) - dt_3d * tsc(2) * c_w(k,i) * & 276 291 ( w(k,-1,i) - w(k,0,i) ) * ddy 277 278 !279 !-- Save old timelevels for the next timestep280 u_m_s(k,:,i) = u(k,-1:1,i)281 v_m_s(k,:,i) = v(k,-1:1,i)282 w_m_s(k,:,i) = w(k,-1:1,i)283 292 284 293 ENDDO … … 289 298 IF ( ibc_uv_b == 0 ) THEN 290 299 u_p(nzb,-1,:) = -u_p(nzb+1,-1,:) 291 v_p(nzb, -1,:) = -v_p(nzb+1,-1,:)300 v_p(nzb,0,:) = -v_p(nzb+1,0,:) 292 301 ELSE 293 302 u_p(nzb,-1,:) = u_p(nzb+1,-1,:) 294 v_p(nzb, -1,:) = v_p(nzb+1,-1,:)295 ENDIF 296 w_p(nzb, ny+1,:) = 0.0303 v_p(nzb,0,:) = v_p(nzb+1,0,:) 304 ENDIF 305 w_p(nzb,-1,:) = 0.0 297 306 298 307 ! … … 300 309 IF ( ibc_uv_t == 0 ) THEN 301 310 u_p(nzt+1,-1,:) = ug(nzt+1) 302 v_p(nzt+1, -1,:)= vg(nzt+1)311 v_p(nzt+1,0,:) = vg(nzt+1) 303 312 ELSE 304 313 u_p(nzt+1,-1,:) = u(nzt,-1,:) 305 v_p(nzt+1, -1,:) = v(nzt,-1,:)314 v_p(nzt+1,0,:) = v(nzt,0,:) 306 315 ENDIF 307 316 w_p(nzt:nzt+1,-1,:) = 0.0 … … 309 318 ENDIF 310 319 311 IF ( outflow_n .AND. & 312 intermediate_timestep_count == intermediate_timestep_count_max ) & 313 THEN 320 IF ( outflow_n ) THEN 314 321 315 322 c_max = dy / dt_3d … … 319 326 320 327 ! 321 !-- First calculate the phase speeds for u,v, and w 322 denom = u_m_n(k,ny,i) - u_m_n(k,ny-1,i) 323 324 IF ( denom /= 0.0 ) THEN 325 c_u = -c_max * ( u(k,ny,i) - u_m_n(k,ny,i) ) / denom 326 IF ( c_u < 0.0 ) THEN 327 c_u = 0.0 328 ELSEIF ( c_u > c_max ) THEN 329 c_u = c_max 330 ENDIF 331 ELSE 332 c_u = c_max 328 !-- Calculate the phase speeds for u,v, and w. In case of using a 329 !-- Runge-Kutta scheme, do this for the first substep only and then 330 !-- reuse this values for the further substeps. 331 IF ( intermediate_timestep_count == 1 ) THEN 332 333 denom = u_m_n(k,ny,i) - u_m_n(k,ny-1,i) 334 335 IF ( denom /= 0.0 ) THEN 336 c_u(k,i) = -c_max * ( u(k,ny,i) - u_m_n(k,ny,i) ) / denom 337 IF ( c_u(k,i) < 0.0 ) THEN 338 c_u(k,i) = 0.0 339 ELSEIF ( c_u(k,i) > c_max ) THEN 340 c_u(k,i) = c_max 341 ENDIF 342 ELSE 343 c_u(k,i) = c_max 344 ENDIF 345 346 denom = v_m_n(k,ny,i) - v_m_n(k,ny-1,i) 347 348 IF ( denom /= 0.0 ) THEN 349 c_v(k,i) = -c_max * ( v(k,ny,i) - v_m_n(k,ny,i) ) / denom 350 IF ( c_v(k,i) < 0.0 ) THEN 351 c_v(k,i) = 0.0 352 ELSEIF ( c_v(k,i) > c_max ) THEN 353 c_v(k,i) = c_max 354 ENDIF 355 ELSE 356 c_v(k,i) = c_max 357 ENDIF 358 359 denom = w_m_n(k,ny,i) - w_m_n(k,ny-1,i) 360 361 IF ( denom /= 0.0 ) THEN 362 c_w(k,i) = -c_max * ( w(k,ny,i) - w_m_n(k,ny,i) ) / denom 363 IF ( c_w(k,i) < 0.0 ) THEN 364 c_w(k,i) = 0.0 365 ELSEIF ( c_w(k,i) > c_max ) THEN 366 c_w(k,i) = c_max 367 ENDIF 368 ELSE 369 c_w(k,i) = c_max 370 ENDIF 371 372 ! 373 !-- Swap timelevels for the next timestep 374 u_m_n(k,:,i) = u(k,ny-1:ny,i) 375 v_m_n(k,:,i) = v(k,ny-1:ny,i) 376 w_m_n(k,:,i) = w(k,ny-1:ny,i) 377 333 378 ENDIF 334 379 335 denom = v_m_n(k,ny,i) - v_m_n(k,ny-1,i)336 337 IF ( denom /= 0.0 ) THEN338 c_v = -c_max * ( v(k,ny,i) - v_m_n(k,ny,i) ) / denom339 IF ( c_v < 0.0 ) THEN340 c_v = 0.0341 ELSEIF ( c_v > c_max ) THEN342 c_v = c_max343 ENDIF344 ELSE345 c_v = c_max346 ENDIF347 348 denom = w_m_n(k,ny,i) - w_m_n(k,ny-1,i)349 350 IF ( denom /= 0.0 ) THEN351 c_w = -c_max * ( w(k,ny,i) - w_m_n(k,ny,i) ) / denom352 IF ( c_w < 0.0 ) THEN353 c_w = 0.0354 ELSEIF ( c_w > c_max ) THEN355 c_w = c_max356 ENDIF357 ELSE358 c_w = c_max359 ENDIF360 361 380 ! 362 381 !-- Calculate the new velocities 363 u_p(k,ny+1,i) = u(k,ny+1,i) - dt_3d * c_u* &382 u_p(k,ny+1,i) = u(k,ny+1,i) - dt_3d * tsc(2) * c_u(k,i) * & 364 383 ( u(k,ny+1,i) - u(k,ny,i) ) * ddy 365 384 366 v_p(k,ny+1,i) = v(k,ny+1,i) - dt_3d * c_v* &385 v_p(k,ny+1,i) = v(k,ny+1,i) - dt_3d * tsc(2) * c_v(k,i) * & 367 386 ( v(k,ny+1,i) - v(k,ny,i) ) * ddy 368 387 369 w_p(k,ny+1,i) = w(k,ny+1,i) - dt_3d * c_w* &388 w_p(k,ny+1,i) = w(k,ny+1,i) - dt_3d * tsc(2) * c_w(k,i) * & 370 389 ( w(k,ny+1,i) - w(k,ny,i) ) * ddy 371 372 !373 !-- Swap timelevels for the next timestep374 u_m_n(k,:,i) = u(k,ny-1:ny+1,i)375 v_m_n(k,:,i) = v(k,ny-1:ny+1,i)376 w_m_n(k,:,i) = w(k,ny-1:ny+1,i)377 390 378 391 ENDDO … … 403 416 ENDIF 404 417 405 IF ( outflow_l .AND. & 406 intermediate_timestep_count == intermediate_timestep_count_max ) & 407 THEN 418 IF ( outflow_l ) THEN 408 419 409 420 c_max = dx / dt_3d … … 413 424 414 425 ! 415 !-- First calculate the phase speeds for u,v, and w 416 denom = u_m_l(k,j,0) - u_m_l(k,j,1) 417 418 IF ( denom /= 0.0 ) THEN 419 c_u = -c_max * ( u(k,j,0) - u_m_r(k,j,0) ) / denom 420 IF ( c_u > 0.0 ) THEN 421 c_u = 0.0 422 ELSEIF ( c_u < -c_max ) THEN 423 c_u = -c_max 424 ENDIF 425 ELSE 426 c_u = -c_max 426 !-- Calculate the phase speeds for u,v, and w. In case of using a 427 !-- Runge-Kutta scheme, do this for the first substep only and then 428 !-- reuse this values for the further substeps. 429 IF ( intermediate_timestep_count == 1 ) THEN 430 431 denom = u_m_l(k,j,1) - u_m_l(k,j,2) 432 433 IF ( denom /= 0.0 ) THEN 434 c_u(k,j) = -c_max * ( u(k,j,1) - u_m_l(k,j,1) ) / denom 435 IF ( c_u(k,j) > 0.0 ) THEN 436 c_u(k,j) = 0.0 437 ELSEIF ( c_u(k,j) < -c_max ) THEN 438 c_u(k,j) = -c_max 439 ENDIF 440 ELSE 441 c_u(k,j) = -c_max 442 ENDIF 443 444 denom = v_m_l(k,j,0) - v_m_l(k,j,1) 445 446 IF ( denom /= 0.0 ) THEN 447 c_v(k,j) = -c_max * ( v(k,j,0) - v_m_l(k,j,0) ) / denom 448 IF ( c_v(k,j) < 0.0 ) THEN 449 c_v(k,j) = 0.0 450 ELSEIF ( c_v(k,j) > c_max ) THEN 451 c_v(k,j) = c_max 452 ENDIF 453 ELSE 454 c_v(k,j) = c_max 455 ENDIF 456 457 denom = w_m_l(k,j,0) - w_m_l(k,j,1) 458 459 IF ( denom /= 0.0 ) THEN 460 c_w(k,j) = -c_max * ( w(k,j,0) - w_m_l(k,j,0) ) / denom 461 IF ( c_w(k,j) < 0.0 ) THEN 462 c_w(k,j) = 0.0 463 ELSEIF ( c_w(k,j) > c_max ) THEN 464 c_w(k,j) = c_max 465 ENDIF 466 ELSE 467 c_w(k,j) = c_max 468 ENDIF 469 470 ! 471 !-- Swap timelevels for the next timestep 472 u_m_l(k,j,:) = u(k,j,1:2) 473 v_m_l(k,j,:) = v(k,j,0:1) 474 w_m_l(k,j,:) = w(k,j,0:1) 475 427 476 ENDIF 428 477 429 denom = v_m_l(k,j,0) - v_m_l(k,j,1)430 431 IF ( denom /= 0.0 ) THEN432 c_v = -c_max * ( v(k,j,0) - v_m_l(k,j,0) ) / denom433 IF ( c_v < 0.0 ) THEN434 c_v = 0.0435 ELSEIF ( c_v > c_max ) THEN436 c_v = c_max437 ENDIF438 ELSE439 c_v = c_max440 ENDIF441 442 denom = w_m_l(k,j,0) - w_m_l(k,j,1)443 444 IF ( denom /= 0.0 ) THEN445 c_w = -c_max * ( w(k,j,0) - w_m_l(k,j,0) ) / denom446 IF ( c_w < 0.0 ) THEN447 c_w = 0.0448 ELSEIF ( c_w > c_max ) THEN449 c_w = c_max450 ENDIF451 ELSE452 c_w = c_max453 ENDIF454 455 478 ! 456 479 !-- Calculate the new velocities 457 u_p(k,j, -1) = u(k,j,-1) + dt_3d * c_u* &458 ( u(k,j, -1) - u(k,j,0) ) * ddx459 460 v_p(k,j,-1) = v(k,j,-1) + dt_3d * c_v* &480 u_p(k,j,0) = u(k,j,0) - dt_3d * tsc(2) * c_u(k,j) * & 481 ( u(k,j,0) - u(k,j,1) ) * ddx 482 483 v_p(k,j,-1) = v(k,j,-1) - dt_3d * tsc(2) * c_v(k,j) * & 461 484 ( v(k,j,-1) - v(k,j,0) ) * ddx 462 485 463 w_p(k,j,-1) = w(k,j,-1) + dt_3d * c_w* &486 w_p(k,j,-1) = w(k,j,-1) - dt_3d * tsc(2) * c_w(k,j) * & 464 487 ( w(k,j,-1) - w(k,j,0) ) * ddx 465 466 !467 !-- Swap timelevels for the next timestep468 u_m_l(k,j,:) = u(k,j,-1:1)469 v_m_l(k,j,:) = v(k,j,-1:1)470 w_m_l(k,j,:) = w(k,j,-1:1)471 488 472 489 ENDDO … … 497 514 ENDIF 498 515 499 IF ( outflow_r .AND. & 500 intermediate_timestep_count == intermediate_timestep_count_max ) & 501 THEN 516 IF ( outflow_r ) THEN 502 517 503 518 c_max = dx / dt_3d … … 507 522 508 523 ! 509 !-- First calculate the phase speeds for u,v, and w 510 denom = u_m_r(k,j,nx) - u_m_r(k,j,nx-1) 511 512 IF ( denom /= 0.0 ) THEN 513 c_u = -c_max * ( u(k,j,nx) - u_m_r(k,j,nx) ) / denom 514 IF ( c_u < 0.0 ) THEN 515 c_u = 0.0 516 ELSEIF ( c_u > c_max ) THEN 517 c_u = c_max 518 ENDIF 519 ELSE 520 c_u = c_max 524 !-- Calculate the phase speeds for u,v, and w. In case of using a 525 !-- Runge-Kutta scheme, do this for the first substep only and then 526 !-- reuse this values for the further substeps. 527 IF ( intermediate_timestep_count == 1 ) THEN 528 529 denom = u_m_r(k,j,nx) - u_m_r(k,j,nx-1) 530 531 IF ( denom /= 0.0 ) THEN 532 c_u(k,j) = -c_max * ( u(k,j,nx) - u_m_r(k,j,nx) ) / denom 533 IF ( c_u(k,j) < 0.0 ) THEN 534 c_u(k,j) = 0.0 535 ELSEIF ( c_u(k,j) > c_max ) THEN 536 c_u(k,j) = c_max 537 ENDIF 538 ELSE 539 c_u(k,j) = c_max 540 ENDIF 541 542 denom = v_m_r(k,j,nx) - v_m_r(k,j,nx-1) 543 544 IF ( denom /= 0.0 ) THEN 545 c_v(k,j) = -c_max * ( v(k,j,nx) - v_m_r(k,j,nx) ) / denom 546 IF ( c_v(k,j) < 0.0 ) THEN 547 c_v(k,j) = 0.0 548 ELSEIF ( c_v(k,j) > c_max ) THEN 549 c_v(k,j) = c_max 550 ENDIF 551 ELSE 552 c_v(k,j) = c_max 553 ENDIF 554 555 denom = w_m_r(k,j,nx) - w_m_r(k,j,nx-1) 556 557 IF ( denom /= 0.0 ) THEN 558 c_w(k,j) = -c_max * ( w(k,j,nx) - w_m_r(k,j,nx) ) / denom 559 IF ( c_w(k,j) < 0.0 ) THEN 560 c_w(k,j) = 0.0 561 ELSEIF ( c_w(k,j) > c_max ) THEN 562 c_w(k,j) = c_max 563 ENDIF 564 ELSE 565 c_w(k,j) = c_max 566 ENDIF 567 568 ! 569 !-- Swap timelevels for the next timestep 570 u_m_r(k,j,:) = u(k,j,nx-1:nx) 571 v_m_r(k,j,:) = v(k,j,nx-1:nx) 572 w_m_r(k,j,:) = w(k,j,nx-1:nx) 573 521 574 ENDIF 522 575 523 denom = v_m_r(k,j,nx) - v_m_r(k,j,nx-1)524 525 IF ( denom /= 0.0 ) THEN526 c_v = -c_max * ( v(k,j,nx) - v_m_r(k,j,nx) ) / denom527 IF ( c_v < 0.0 ) THEN528 c_v = 0.0529 ELSEIF ( c_v > c_max ) THEN530 c_v = c_max531 ENDIF532 ELSE533 c_v = c_max534 ENDIF535 536 denom = w_m_r(k,j,nx) - w_m_r(k,j,nx-1)537 538 IF ( denom /= 0.0 ) THEN539 c_w = -c_max * ( w(k,j,nx) - w_m_r(k,j,nx) ) / denom540 IF ( c_w < 0.0 ) THEN541 c_w = 0.0542 ELSEIF ( c_w > c_max ) THEN543 c_w = c_max544 ENDIF545 ELSE546 c_w = c_max547 ENDIF548 549 576 ! 550 577 !-- Calculate the new velocities 551 u_p(k,j,nx+1) = u(k,j,nx+1) - dt_3d * c_u* &578 u_p(k,j,nx+1) = u(k,j,nx+1) - dt_3d * tsc(2) * c_u(k,j) * & 552 579 ( u(k,j,nx+1) - u(k,j,nx) ) * ddx 553 580 554 v_p(k,j,nx+1) = v(k,j,nx+1) - dt_3d * c_v* &581 v_p(k,j,nx+1) = v(k,j,nx+1) - dt_3d * tsc(2) * c_v(k,j) * & 555 582 ( v(k,j,nx+1) - v(k,j,nx) ) * ddx 556 583 557 w_p(k,j,nx+1) = w(k,j,nx+1) - dt_3d * c_w* &584 w_p(k,j,nx+1) = w(k,j,nx+1) - dt_3d * tsc(2) * c_w(k,j) * & 558 585 ( w(k,j,nx+1) - w(k,j,nx) ) * ddx 559 560 !561 !-- Swap timelevels for the next timestep562 u_m_r(k,j,:) = u(k,j,nx-1:nx+1)563 v_m_r(k,j,:) = v(k,j,nx-1:nx+1)564 w_m_r(k,j,:) = w(k,j,nx-1:nx+1)565 586 566 587 ENDDO -
palm/trunk/SOURCE/buoyancy.f90
r98 r106 4 4 ! Actual revisions: 5 5 ! ----------------- 6 ! 6 ! i loop for u-component (sloping surface) is starting from nxlu (needed for 7 ! non-cyclic boundary conditions) 7 8 ! 8 9 ! Former revisions: … … 105 106 IF ( wind_component == 1 ) THEN 106 107 107 DO i = nxl , nxr108 DO i = nxlu, nxr 108 109 DO j = nys, nyn 109 110 DO k = nzb_s_inner(j,i)+1, nzt-1 -
palm/trunk/SOURCE/check_parameters.f90
r104 r106 5 5 ! ----------------- 6 6 ! Check coupling_mode and set default (obligatory) values (like boundary 7 ! conditions for temperature and fluxes) in case of coupled runs 7 ! conditions for temperature and fluxes) in case of coupled runs. 8 ! +profiles for w*p* and w"e 8 9 ! Bugfix: Error message concerning output of particle concentration (pc) 9 10 ! modified … … 1968 1969 dopr_index(i) = 56 1969 1970 dopr_unit(i) = 'm2/s3' 1970 hom(:,2,56,:) = SPREAD( z u, 2, statistic_regions+1 )1971 hom(:,2,56,:) = SPREAD( zw, 2, statistic_regions+1 ) 1971 1972 1972 1973 CASE ( 'w"e/dz' ) … … 2051 2052 hom(:,2,67,:) = SPREAD( zw, 2, statistic_regions+1 ) 2052 2053 ENDIF 2054 2055 CASE ( 'w*p*' ) 2056 dopr_index(i) = 68 2057 dopr_unit(i) = 'm3/s3' 2058 hom(:,2,68,:) = SPREAD( zu, 2, statistic_regions+1 ) 2059 2060 CASE ( 'w"e' ) 2061 dopr_index(i) = 69 2062 dopr_unit(i) = 'm3/s3' 2063 hom(:,2,69,:) = SPREAD( zu, 2, statistic_regions+1 ) 2053 2064 2054 2065 -
palm/trunk/SOURCE/coriolis.f90
r77 r106 4 4 ! Actual revisions: 5 5 ! ----------------- 6 ! 6 ! loops for u and v are starting from index nxlu, nysv, respectively (needed 7 ! for non-cyclic boundary conditions) 7 8 ! 8 9 ! Former revisions: … … 60 61 !-- u-component 61 62 CASE ( 1 ) 62 DO i = nxl , nxr63 DO i = nxlu, nxr 63 64 DO j = nys, nyn 64 65 DO k = nzb_u_inner(j,i)+1, nzt … … 78 79 CASE ( 2 ) 79 80 DO i = nxl, nxr 80 DO j = nys , nyn81 DO j = nysv, nyn 81 82 DO k = nzb_v_inner(j,i)+1, nzt 82 83 tend(k,j,i) = tend(k,j,i) - f * ( 0.25 * & -
palm/trunk/SOURCE/data_output_dvrp.f90
r86 r106 151 151 ! 152 152 !-- Definition of characteristics of particle material 153 CALL DVRP_MATERIAL_RGB( m-1, 1, 0.1, 0.7, 0.1, 0.0 ) 153 ! CALL DVRP_MATERIAL_RGB( m-1, 1, 0.1, 0.7, 0.1, 0.0 ) 154 CALL DVRP_MATERIAL_RGB( m-1, 1, 0.0, 0.0, 0.0, 0.0 ) 154 155 155 156 ! … … 318 319 ENDDO 319 320 ENDDO 320 DO k = nzb, nz_do3d 321 DO j = nys+1, nyn 322 DO i = nxl, nxr+1 323 local_pf(i,j,k) = 0.25 * local_pf(i,j-1,k) + & 324 0.50 * local_pf(i,j,k) + & 325 0.25 * local_pf(i,j+1,k) 326 ENDDO 327 ENDDO 328 ENDDO 321 ! Averaging for Langmuir circulation 322 ! DO k = nzb, nz_do3d 323 ! DO j = nys+1, nyn 324 ! DO i = nxl, nxr+1 325 ! local_pf(i,j,k) = 0.25 * local_pf(i,j-1,k) + & 326 ! 0.50 * local_pf(i,j,k) + & 327 ! 0.25 * local_pf(i,j+1,k) 328 ! ENDDO 329 ! ENDDO 330 ! ENDDO 329 331 330 332 CASE ( 'p', 'p_xy', 'p_xz', 'p_yz' ) … … 475 477 CALL DVRP_DATA( m-1, local_pf, 1, nx_dvrp, ny_dvrp, nz_dvrp, & 476 478 cyclic_dvrp, cyclic_dvrp, cyclic_dvrp ) 477 CALL DVRP_SLICER( m-1, section_mode, slicer_position ) 479 ! CALL DVRP_SLICER( m-1, section_mode, slicer_position ) 480 CALL DVRP_SLICER( m-1, 2, 1.0 ) 481 WRITE (9,*) 'nx_dvrp=', nx_dvrp 482 WRITE (9,*) 'ny_dvrp=', ny_dvrp 483 WRITE (9,*) 'nz_dvrp=', nz_dvrp 484 WRITE (9,*) 'section_mode=', section_mode 485 WRITE (9,*) 'slicer_position=', slicer_position 486 CALL local_flush( 9 ) 487 478 488 CALL DVRP_VISUALIZE( m-1, 2, dvrp_filecount ) 479 489 -
palm/trunk/SOURCE/diffusion_u.f90
r103 r106 4 4 ! Actual revisions: 5 5 ! ----------------- 6 ! Momentumflux at top (uswst) included as boundary condition 7 ! 6 ! Momentumflux at top (uswst) included as boundary condition, 7 ! i loop is starting from nxlu (needed for non-cyclic boundary conditions) 8 ! 8 9 ! Former revisions: 9 10 ! ----------------- … … 79 80 ENDIF 80 81 81 DO i = nxl , nxr82 DO i = nxlu, nxr 82 83 DO j = nys,nyn 83 84 ! -
palm/trunk/SOURCE/diffusion_v.f90
r103 r106 4 4 ! Actual revisions: 5 5 ! ----------------- 6 ! Momentumflux at top (vswst) included as boundary condition 6 ! Momentumflux at top (vswst) included as boundary condition, 7 ! j loop is starting from nysv (needed for non-cyclic boundary conditions) 7 8 ! 8 9 ! Former revisions: … … 78 79 79 80 DO i = nxl, nxr 80 DO j = nys , nyn81 DO j = nysv, nyn 81 82 ! 82 83 !-- Compute horizontal diffusion -
palm/trunk/SOURCE/flow_statistics.f90
r102 r106 4 4 ! Actual revisions: 5 5 ! ----------------- 6 ! Prescribed momentum fluxes at the top surface are used 6 ! Prescribed momentum fluxes at the top surface are used, 7 ! profiles for w*p* and w"e are calculated 7 8 ! 8 9 ! Former revisions: … … 620 621 ! 621 622 !-- Divergence of vertical flux of resolved scale energy and pressure 622 !-- fluctuations. First calculate the products, then the divergence. 623 !-- fluctuations as well as flux of pressure fluctuation itself (68). 624 !-- First calculate the products, then the divergence. 623 625 !-- Calculation is time consuming. Do it only, if profiles shall be plotted. 624 IF ( hom(nzb+1,2,55,0) /= 0.0 ) THEN626 IF ( hom(nzb+1,2,55,0) /= 0.0 .OR. hom(nzb+1,2,68,0) /= 0.0 ) THEN 625 627 626 628 sums_ll = 0.0 ! local array … … 654 656 sums_l(k,55,tn) = ( sums_ll(k,1) - sums_ll(k-1,1) ) * ddzw(k) 655 657 sums_l(k,56,tn) = ( sums_ll(k,2) - sums_ll(k-1,2) ) * ddzw(k) 658 sums_l(k,68,tn) = sums_ll(k,2) 656 659 ENDDO 657 660 sums_l(nzb,55,tn) = sums_l(nzb+1,55,tn) 658 661 sums_l(nzb,56,tn) = sums_l(nzb+1,56,tn) 659 660 ENDIF 661 662 ! 663 !-- Divergence of vertical flux of SGS TKE 664 IF ( hom(nzb+1,2,57,0) /= 0.0 ) THEN 662 sums_l(nzb,68,tn) = 0.0 ! because w* = 0 at nzb 663 664 ENDIF 665 666 ! 667 !-- Divergence of vertical flux of SGS TKE and the flux itself (69) 668 IF ( hom(nzb+1,2,57,0) /= 0.0 .OR. hom(nzb+1,2,69,0) /= 0.0 ) THEN 665 669 666 670 !$OMP DO … … 669 673 DO k = nzb_s_outer(j,i)+1, nzt 670 674 671 sums_l(k,57,tn) = sums_l(k,57,tn) + (&675 sums_l(k,57,tn) = sums_l(k,57,tn) - 0.5 * ( & 672 676 (km(k,j,i)+km(k+1,j,i)) * (e(k+1,j,i)-e(k,j,i)) * ddzu(k+1) & 673 677 - (km(k-1,j,i)+km(k,j,i)) * (e(k,j,i)-e(k-1,j,i)) * ddzu(k) & 674 ) * ddzw(k) 678 ) * ddzw(k) 679 680 sums_l(k,69,tn) = sums_l(k,69,tn) - 0.5 * ( & 681 (km(k,j,i)+km(k+1,j,i)) * (e(k+1,j,i)-e(k,j,i)) * ddzu(k+1) & 682 ) 675 683 676 684 ENDDO … … 678 686 ENDDO 679 687 sums_l(nzb,57,tn) = sums_l(nzb+1,57,tn) 688 sums_l(nzb,69,tn) = sums_l(nzb+1,69,tn) 680 689 681 690 ENDIF … … 832 841 hom(:,1,55,sr) = sums(:,55) ! w*u*u*/dz 833 842 hom(:,1,56,sr) = sums(:,56) ! w*p*/dz 834 hom(:,1,57,sr) = sums(:,57) ! w"e/dz843 hom(:,1,57,sr) = sums(:,57) ! ( w"e + w"p"/rho )/dz 835 844 hom(:,1,58,sr) = sums(:,58) ! u"pt" 836 845 hom(:,1,59,sr) = sums(:,59) ! u*pt* … … 843 852 hom(:,1,66,sr) = sums(:,66) ! w*sa* 844 853 hom(:,1,67,sr) = sums(:,65) + sums(:,66) ! wsa 854 hom(:,1,68,sr) = sums(:,68) ! w*p* 855 hom(:,1,69,sr) = sums(:,69) ! w"e + w"p"/rho 845 856 846 857 hom(:,1,pr_palm-1,sr) = sums(:,pr_palm-1) -
palm/trunk/SOURCE/init_3d_model.f90
r102 r106 7 7 ! Actual revisions: 8 8 ! ----------------- 9 ! Flux initialization in case of coupled runs, +momentum fluxes at top boundary 9 ! Flux initialization in case of coupled runs, +momentum fluxes at top boundary, 10 ! +arrays for phase speed c_u, c_v, c_w, indices for u|v|w_m_l|r changed 10 11 ! 11 12 ! Former revisions: … … 233 234 234 235 ! 235 !-- Arrays to store velocity data from t-dt needed for radiation boundary236 !-- conditions236 !-- Arrays to store velocity data from t-dt and the phase speeds which 237 !-- are needed for radiation boundary conditions 237 238 IF ( outflow_l ) THEN 238 ALLOCATE( u_m_l(nzb:nzt+1,nys-1:nyn+1, -1:1), &239 v_m_l(nzb:nzt+1,nys-1:nyn+1, -1:1), &240 w_m_l(nzb:nzt+1,nys-1:nyn+1, -1:1) )239 ALLOCATE( u_m_l(nzb:nzt+1,nys-1:nyn+1,1:2), & 240 v_m_l(nzb:nzt+1,nys-1:nyn+1,0:1), & 241 w_m_l(nzb:nzt+1,nys-1:nyn+1,0:1) ) 241 242 ENDIF 242 243 IF ( outflow_r ) THEN 243 ALLOCATE( u_m_r(nzb:nzt+1,nys-1:nyn+1,nx-1:nx+1), & 244 v_m_r(nzb:nzt+1,nys-1:nyn+1,nx-1:nx+1), & 245 w_m_r(nzb:nzt+1,nys-1:nyn+1,nx-1:nx+1) ) 244 ALLOCATE( u_m_r(nzb:nzt+1,nys-1:nyn+1,nx-1:nx), & 245 v_m_r(nzb:nzt+1,nys-1:nyn+1,nx-1:nx), & 246 w_m_r(nzb:nzt+1,nys-1:nyn+1,nx-1:nx) ) 247 ENDIF 248 IF ( outflow_l .OR. outflow_r ) THEN 249 ALLOCATE( c_u(nzb:nzt+1,nys-1:nyn+1), c_v(nzb:nzt+1,nys-1:nyn+1), & 250 c_w(nzb:nzt+1,nys-1:nyn+1) ) 246 251 ENDIF 247 252 IF ( outflow_s ) THEN 248 ALLOCATE( u_m_s(nzb:nzt+1, -1:1,nxl-1:nxr+1), &249 v_m_s(nzb:nzt+1, -1:1,nxl-1:nxr+1), &250 w_m_s(nzb:nzt+1, -1:1,nxl-1:nxr+1) )253 ALLOCATE( u_m_s(nzb:nzt+1,0:1,nxl-1:nxr+1), & 254 v_m_s(nzb:nzt+1,1:2,nxl-1:nxr+1), & 255 w_m_s(nzb:nzt+1,0:1,nxl-1:nxr+1) ) 251 256 ENDIF 252 257 IF ( outflow_n ) THEN 253 ALLOCATE( u_m_n(nzb:nzt+1,ny-1:ny+1,nxl-1:nxr+1), & 254 v_m_n(nzb:nzt+1,ny-1:ny+1,nxl-1:nxr+1), & 255 w_m_n(nzb:nzt+1,ny-1:ny+1,nxl-1:nxr+1) ) 258 ALLOCATE( u_m_n(nzb:nzt+1,ny-1:ny,nxl-1:nxr+1), & 259 v_m_n(nzb:nzt+1,ny-1:ny,nxl-1:nxr+1), & 260 w_m_n(nzb:nzt+1,ny-1:ny,nxl-1:nxr+1) ) 261 ENDIF 262 IF ( outflow_s .OR. outflow_n ) THEN 263 ALLOCATE( c_u(nzb:nzt+1,nxl-1:nxr+1), c_v(nzb:nzt+1,nxl-1:nxr+1), & 264 c_w(nzb:nzt+1,nxl-1:nxr+1) ) 256 265 ENDIF 257 266 … … 808 817 !-- Initialize old timelevels needed for radiation boundary conditions 809 818 IF ( outflow_l ) THEN 810 u_m_l(:,:,:) = u(:,:, -1:1)811 v_m_l(:,:,:) = v(:,:, -1:1)812 w_m_l(:,:,:) = w(:,:, -1:1)819 u_m_l(:,:,:) = u(:,:,1:2) 820 v_m_l(:,:,:) = v(:,:,0:1) 821 w_m_l(:,:,:) = w(:,:,0:1) 813 822 ENDIF 814 823 IF ( outflow_r ) THEN 815 u_m_r(:,:,:) = u(:,:,nx-1:nx +1)816 v_m_r(:,:,:) = v(:,:,nx-1:nx +1)817 w_m_r(:,:,:) = w(:,:,nx-1:nx +1)824 u_m_r(:,:,:) = u(:,:,nx-1:nx) 825 v_m_r(:,:,:) = v(:,:,nx-1:nx) 826 w_m_r(:,:,:) = w(:,:,nx-1:nx) 818 827 ENDIF 819 828 IF ( outflow_s ) THEN 820 u_m_s(:,:,:) = u(:, -1:1,:)821 v_m_s(:,:,:) = v(:, -1:1,:)822 w_m_s(:,:,:) = w(:, -1:1,:)829 u_m_s(:,:,:) = u(:,0:1,:) 830 v_m_s(:,:,:) = v(:,1:2,:) 831 w_m_s(:,:,:) = w(:,0:1,:) 823 832 ENDIF 824 833 IF ( outflow_n ) THEN 825 u_m_n(:,:,:) = u(:,ny-1:ny +1,:)826 v_m_n(:,:,:) = v(:,ny-1:ny +1,:)827 w_m_n(:,:,:) = w(:,ny-1:ny +1,:)834 u_m_n(:,:,:) = u(:,ny-1:ny,:) 835 v_m_n(:,:,:) = v(:,ny-1:ny,:) 836 w_m_n(:,:,:) = w(:,ny-1:ny,:) 828 837 ENDIF 829 838 -
palm/trunk/SOURCE/init_particles.f90
r83 r106 4 4 ! Actual revisions: 5 5 ! ----------------- 6 ! 6 ! variable iran replaced by iran_part 7 7 ! 8 8 ! Former revisions: … … 354 354 IF ( random_start_position ) THEN 355 355 356 iran = iran + myid ! Random positions should be different on357 ! different PEs358 359 356 DO n = 1, number_of_initial_particles 360 357 IF ( psl(particles(n)%group) /= psr(particles(n)%group) ) THEN 361 358 particles(n)%x = particles(n)%x + & 362 ( random_function( iran ) - 0.5 ) * &359 ( random_function( iran_part ) - 0.5 ) * & 363 360 pdx(particles(n)%group) 364 361 IF ( particles(n)%x <= ( nxl - 0.5 ) * dx ) THEN … … 370 367 IF ( pss(particles(n)%group) /= psn(particles(n)%group) ) THEN 371 368 particles(n)%y = particles(n)%y + & 372 ( random_function( iran ) - 0.5 ) * &369 ( random_function( iran_part ) - 0.5 ) * & 373 370 pdy(particles(n)%group) 374 371 IF ( particles(n)%y <= ( nys - 0.5 ) * dy ) THEN … … 380 377 IF ( psb(particles(n)%group) /= pst(particles(n)%group) ) THEN 381 378 particles(n)%z = particles(n)%z + & 382 ( random_function( iran ) - 0.5 ) * &379 ( random_function( iran_part ) - 0.5 ) * & 383 380 pdz(particles(n)%group) 384 381 ENDIF -
palm/trunk/SOURCE/init_pegrid.f90
r104 r106 5 5 ! ----------------- 6 6 ! Intercommunicator (comm_inter) and derived data type (type_xy) for 7 ! coupled model runs created 7 ! coupled model runs created, 8 ! indices nxlu and nysv are calculated (needed for non-cyclic boundary 9 ! conditions) 8 10 ! 9 11 ! Former revisions: … … 818 820 ! 819 821 !-- Setting of flags for inflow/outflow conditions in case of non-cyclic 820 !-- horizontal boundary conditions. Set variables for extending array u (v) 821 !-- by one gridpoint on the left/rightmost (northest/southest) processor 822 !-- horizontal boundary conditions. 822 823 IF ( pleft == MPI_PROC_NULL ) THEN 823 824 IF ( bc_lr == 'dirichlet/radiation' ) THEN … … 869 870 ENDIF 870 871 #endif 872 ! 873 !-- At the outflow, u or v, respectively, have to be calculated for one grid 874 !-- point less. 875 IF ( outflow_l ) THEN 876 nxlu = nxl + 1 877 ELSE 878 nxlu = nxl 879 ENDIF 880 IF ( outflow_s ) THEN 881 nysv = nys + 1 882 ELSE 883 nysv = nys 884 ENDIF 871 885 872 886 IF ( psolver == 'poisfft_hybrid' ) THEN -
palm/trunk/SOURCE/modules.f90
r103 r106 6 6 ! ----------------- 7 7 ! +comm_inter, constant_top_momentumflux, coupling_char, coupling_mode, 8 ! dt_coupling, ngp_xy, port_name, time_coupling, top_momentumflux_u|v,9 ! t ype_xy, uswst*, vswst*8 ! c_u, c_v, c_w, dt_coupling, ngp_xy, nxlu, nysv, port_name, time_coupling, 9 ! top_momentumflux_u|v, type_xy, uswst*, vswst* 10 10 ! 11 11 ! Former revisions: … … 105 105 106 106 REAL, DIMENSION(:,:), ALLOCATABLE :: & 107 dzu_mg, dzw_mg, f1_mg, f2_mg, f3_mg, pt_slope_ref, qs, ts, us, z0 107 c_u, c_v, c_w, dzu_mg, dzw_mg, f1_mg, f2_mg, f3_mg, pt_slope_ref, & 108 qs, ts, us, z0 108 109 109 110 REAL, DIMENSION(:,:), ALLOCATABLE, TARGET :: & … … 565 566 !------------------------------------------------------------------------------! 566 567 567 INTEGER :: ngp_sums, nnx, nx = 0, nxa, nxl, nx r, nxra, nny, ny = 0, nya,&568 ny n, nyna, nys, nnz, nz = 0, nza, nzb, nzb_diff, nzt, nzta, &569 nzt _diff568 INTEGER :: ngp_sums, nnx, nx = 0, nxa, nxl, nxlu, nxr, nxra, nny, ny = 0, & 569 nya, nyn, nyna, nys, nysv, nnz, nz = 0, nza, nzb, nzb_diff, & 570 nzt, nzta, nzt_diff 570 571 571 572 INTEGER, DIMENSION(:), ALLOCATABLE :: & -
palm/trunk/SOURCE/pres.f90
r90 r106 4 4 ! Actual revisions: 5 5 ! ----------------- 6 ! 6 ! Volume flow conservation added for the remaining three outflow boundaries 7 7 ! 8 8 ! Former revisions: … … 74 74 !-- Conserve the volume flow at the outflow in case of non-cyclic lateral 75 75 !-- boundary conditions 76 IF ( conserve_volume_flow .AND. bc_ns == 'radiation/dirichlet') THEN 76 !-- WARNING: so far, this conservation does not work at the left/south 77 !-- boundary if the topography at the inflow differs from that at the 78 !-- outflow! For this case, volume_flow_area needs adjustment! 79 ! 80 !-- Left/right 81 IF ( conserve_volume_flow .AND. ( outflow_l .OR. outflow_r ) ) THEN 82 83 volume_flow(1) = 0.0 84 volume_flow_l(1) = 0.0 85 86 IF ( outflow_l ) THEN 87 i = 0 88 ELSEIF ( outflow_r ) THEN 89 i = nx+1 90 ENDIF 91 92 DO j = nys, nyn 93 ! 94 !-- Sum up the volume flow through the south/north boundary 95 DO k = nzb_2d(j,i) + 1, nzt 96 volume_flow_l(1) = volume_flow_l(1) + u(k,j,i) * dzu(k) 97 ENDDO 98 ENDDO 99 100 #if defined( __parallel ) 101 CALL MPI_ALLREDUCE( volume_flow_l(1), volume_flow(1), 1, MPI_REAL, & 102 MPI_SUM, comm1dy, ierr ) 103 #else 104 volume_flow = volume_flow_l 105 #endif 106 volume_flow_offset(1) = ( volume_flow_initial(1) - volume_flow(1) ) & 107 / volume_flow_area(1) 108 109 DO j = nys, nyn 110 DO k = nzb_v_inner(j,i) + 1, nzt 111 u(k,j,i) = u(k,j,i) + volume_flow_offset(1) 112 ENDDO 113 ENDDO 114 115 CALL exchange_horiz( u ) 116 117 ENDIF 118 119 ! 120 !-- South/north 121 IF ( conserve_volume_flow .AND. ( outflow_n .OR. outflow_s ) ) THEN 77 122 78 123 volume_flow(2) = 0.0 79 124 volume_flow_l(2) = 0.0 80 125 81 IF ( nyn == ny ) THEN 126 IF ( outflow_s ) THEN 127 j = 0 128 ELSEIF ( outflow_n ) THEN 82 129 j = ny+1 83 DO i = nxl, nxr 84 ! 85 !-- Sum up the volume flow through the north boundary 86 DO k = nzb_2d(j,i) + 1, nzt 87 volume_flow_l(2) = volume_flow_l(2) + v(k,j,i) * dzu(k) 88 ENDDO 89 ENDDO 90 ENDIF 130 ENDIF 131 132 DO i = nxl, nxr 133 ! 134 !-- Sum up the volume flow through the south/north boundary 135 DO k = nzb_2d(j,i) + 1, nzt 136 volume_flow_l(2) = volume_flow_l(2) + v(k,j,i) * dzu(k) 137 ENDDO 138 ENDDO 139 91 140 #if defined( __parallel ) 92 141 CALL MPI_ALLREDUCE( volume_flow_l(2), volume_flow(2), 1, MPI_REAL, & … … 96 145 #endif 97 146 volume_flow_offset(2) = ( volume_flow_initial(2) - volume_flow(2) ) & 98 / volume_flow_area(2) 99 100 IF ( outflow_n ) THEN 101 j = nyn+1 102 DO i = nxl, nxr 103 DO k = nzb_v_inner(j,i) + 1, nzt 104 v(k,j,i) = v(k,j,i) + volume_flow_offset(2) 105 ENDDO 106 ENDDO 107 ENDIF 147 / volume_flow_area(2) 148 149 DO i = nxl, nxr 150 DO k = nzb_v_inner(j,i) + 1, nzt 151 v(k,j,i) = v(k,j,i) + volume_flow_offset(2) 152 ENDDO 153 ENDDO 108 154 109 155 CALL exchange_horiz( v ) -
palm/trunk/SOURCE/production_e.f90
r98 r106 4 4 ! Actual revisions: 5 5 ! ----------------- 6 ! 6 ! Bugfix: wrong sign removed from the buoyancy production term in the case 7 ! use_reference = .T., 8 ! u_0 and v_0 are calculated for nxr+1, nyn+1 also (otherwise these values are 9 ! not available in case of non-cyclic boundary conditions) 7 10 ! 8 11 ! Former revisions: … … 386 389 DO j = nys, nyn 387 390 DO k = nzb_diff_s_inner(j,i), nzt_diff 388 tend(k,j,i) = tend(k,j,i) +&391 tend(k,j,i) = tend(k,j,i) - & 389 392 kh(k,j,i) * g / pt_reference * & 390 393 ( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k) … … 976 979 !-- (otherwise the timestep may be significantly reduced by large 977 980 !-- surface winds). 981 !-- Upper bounds are nxr+1 and nyn+1 because otherwise these values are 982 !-- not available in case of non-cyclic boundary conditions. 978 983 !-- WARNING: the exact analytical solution would require the determination 979 984 !-- of the eddy diffusivity by km = u* * kappa * zp / phi_m. 980 985 !$OMP PARALLEL DO PRIVATE( ku, kv ) 981 DO i = nxl, nxr 982 DO j = nys, nyn 986 DO i = nxl, nxr+1 987 DO j = nys, nyn+1 983 988 984 989 ku = nzb_u_inner(j,i)+1 -
palm/trunk/SOURCE/prognostic_equations.f90
r102 r106 4 4 ! Actual revisions: 5 5 ! ----------------- 6 ! +uswst, vswst as arguments in calls of diffusion_u|v 6 ! +uswst, vswst as arguments in calls of diffusion_u|v, 7 ! loops for u and v are starting from index nxlu, nysv, respectively (needed 8 ! for non-cyclic boundary conditions) 7 9 ! 8 10 ! Former revisions: … … 133 135 ! 134 136 !-- u-tendency terms with no communication 135 DO i = nxl , nxr137 DO i = nxlu, nxr 136 138 DO j = nys, nyn 137 139 ! … … 202 204 !-- v-tendency terms with no communication 203 205 DO i = nxl, nxr 204 DO j = nys , nyn206 DO j = nysv, nyn 205 207 ! 206 208 !-- Tendency terms … … 806 808 ! 807 809 !-- Tendency terms for u-velocity component 808 IF ( j < nyn+1) THEN810 IF ( .NOT. outflow_l .OR. i > nxl ) THEN 809 811 810 812 tend(:,j,i) = 0.0 … … 857 859 ! 858 860 !-- Tendency terms for v-velocity component 859 IF ( i < nxr+1) THEN861 IF ( .NOT. outflow_s .OR. j > nys ) THEN 860 862 861 863 tend(:,j,i) = 0.0 … … 906 908 ! 907 909 !-- Tendency terms for w-velocity component 908 IF ( i < nxr+1 .AND. j < nyn+1 ) THEN 909 910 tend(:,j,i) = 0.0 911 IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN 912 CALL advec_w_pw( i, j ) 913 ELSE 914 CALL advec_w_up( i, j ) 915 ENDIF 916 IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) & 917 THEN 918 CALL diffusion_w( i, j, ddzu, ddzw, km_m, km_damp_x, & 919 km_damp_y, tend, u_m, v_m, w_m ) 920 ELSE 921 CALL diffusion_w( i, j, ddzu, ddzw, km, km_damp_x, km_damp_y, & 922 tend, u, v, w ) 923 ENDIF 924 CALL coriolis( i, j, 3 ) 925 IF ( ocean ) THEN 926 CALL buoyancy( i, j, rho, prho_reference, 3, 64 ) 927 ELSE 928 IF ( .NOT. humidity ) THEN 929 CALL buoyancy( i, j, pt, pt_reference, 3, 4 ) 930 ELSE 931 CALL buoyancy( i, j, vpt, pt_reference, 3, 44 ) 932 ENDIF 933 ENDIF 934 CALL user_actions( i, j, 'w-tendency' ) 935 936 ! 937 !-- Prognostic equation for w-velocity component 938 DO k = nzb_w_inner(j,i)+1, nzt-1 939 w_p(k,j,i) = ( 1.0-tsc(1) ) * w_m(k,j,i) + tsc(1) * w(k,j,i) + & 940 dt_3d * ( & 941 tsc(2) * tend(k,j,i) + tsc(3) * tw_m(k,j,i) & 942 - tsc(4) * ( p(k+1,j,i) - p(k,j,i) ) * ddzu(k+1) & 943 ) - & 944 tsc(5) * rdf(k) * w(k,j,i) 945 ENDDO 946 947 ! 948 !-- Calculate tendencies for the next Runge-Kutta step 949 IF ( timestep_scheme(1:5) == 'runge' ) THEN 950 IF ( intermediate_timestep_count == 1 ) THEN 951 DO k = nzb_w_inner(j,i)+1, nzt-1 952 tw_m(k,j,i) = tend(k,j,i) 953 ENDDO 954 ELSEIF ( intermediate_timestep_count < & 955 intermediate_timestep_count_max ) THEN 956 DO k = nzb_w_inner(j,i)+1, nzt-1 957 tw_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tw_m(k,j,i) 958 ENDDO 959 ENDIF 960 ENDIF 961 962 ! 963 !-- Tendency terms for potential temperature 964 tend(:,j,i) = 0.0 965 IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN 966 CALL advec_s_pw( i, j, pt ) 967 ELSE 968 CALL advec_s_up( i, j, pt ) 969 ENDIF 970 IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) & 971 THEN 972 CALL diffusion_s( i, j, ddzu, ddzw, kh_m, pt_m, shf_m, & 973 tswst_m, tend ) 974 ELSE 975 CALL diffusion_s( i, j, ddzu, ddzw, kh, pt, shf, tswst, tend ) 976 ENDIF 977 978 ! 979 !-- If required compute heating/cooling due to long wave radiation 980 !-- processes 981 IF ( radiation ) THEN 982 CALL calc_radiation( i, j ) 983 ENDIF 984 985 ! 986 !-- If required compute impact of latent heat due to precipitation 987 IF ( precipitation ) THEN 988 CALL impact_of_latent_heat( i, j ) 989 ENDIF 990 CALL user_actions( i, j, 'pt-tendency' ) 991 992 ! 993 !-- Prognostic equation for potential temperature 994 DO k = nzb_s_inner(j,i)+1, nzt 995 pt_p(k,j,i) = ( 1.0-tsc(1) ) * pt_m(k,j,i) + tsc(1)*pt(k,j,i) +& 996 dt_3d * ( & 997 tsc(2) * tend(k,j,i) + tsc(3) * tpt_m(k,j,i) & 998 ) - & 999 tsc(5) * rdf(k) * ( pt(k,j,i) - pt_init(k) ) 1000 ENDDO 1001 1002 ! 1003 !-- Calculate tendencies for the next Runge-Kutta step 1004 IF ( timestep_scheme(1:5) == 'runge' ) THEN 1005 IF ( intermediate_timestep_count == 1 ) THEN 1006 DO k = nzb_s_inner(j,i)+1, nzt 1007 tpt_m(k,j,i) = tend(k,j,i) 1008 ENDDO 1009 ELSEIF ( intermediate_timestep_count < & 1010 intermediate_timestep_count_max ) THEN 1011 DO k = nzb_s_inner(j,i)+1, nzt 1012 tpt_m(k,j,i) = -9.5625 * tend(k,j,i) + & 1013 5.3125 * tpt_m(k,j,i) 1014 ENDDO 1015 ENDIF 1016 ENDIF 1017 1018 ! 1019 !-- If required, compute prognostic equation for salinity 1020 IF ( ocean ) THEN 1021 1022 ! 1023 !-- Tendency-terms for salinity 910 1024 tend(:,j,i) = 0.0 911 IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN 912 CALL advec_w_pw( i, j ) 913 ELSE 914 CALL advec_w_up( i, j ) 915 ENDIF 916 IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) & 1025 IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) & 917 1026 THEN 918 CALL diffusion_w( i, j, ddzu, ddzw, km_m, km_damp_x, & 919 km_damp_y, tend, u_m, v_m, w_m ) 920 ELSE 921 CALL diffusion_w( i, j, ddzu, ddzw, km, km_damp_x, km_damp_y, & 922 tend, u, v, w ) 923 ENDIF 924 CALL coriolis( i, j, 3 ) 925 IF ( ocean ) THEN 926 CALL buoyancy( i, j, rho, prho_reference, 3, 64 ) 927 ELSE 928 IF ( .NOT. humidity ) THEN 929 CALL buoyancy( i, j, pt, pt_reference, 3, 4 ) 930 ELSE 931 CALL buoyancy( i, j, vpt, pt_reference, 3, 44 ) 932 ENDIF 933 ENDIF 934 CALL user_actions( i, j, 'w-tendency' ) 935 936 ! 937 !-- Prognostic equation for w-velocity component 938 DO k = nzb_w_inner(j,i)+1, nzt-1 939 w_p(k,j,i) = ( 1.0-tsc(1) ) * w_m(k,j,i) + tsc(1) * w(k,j,i) + & 940 dt_3d * ( & 941 tsc(2) * tend(k,j,i) + tsc(3) * tw_m(k,j,i) & 942 - tsc(4) * ( p(k+1,j,i) - p(k,j,i) ) * ddzu(k+1) & 943 ) - & 944 tsc(5) * rdf(k) * w(k,j,i) 945 ENDDO 946 947 ! 948 !-- Calculate tendencies for the next Runge-Kutta step 949 IF ( timestep_scheme(1:5) == 'runge' ) THEN 950 IF ( intermediate_timestep_count == 1 ) THEN 951 DO k = nzb_w_inner(j,i)+1, nzt-1 952 tw_m(k,j,i) = tend(k,j,i) 953 ENDDO 954 ELSEIF ( intermediate_timestep_count < & 955 intermediate_timestep_count_max ) THEN 956 DO k = nzb_w_inner(j,i)+1, nzt-1 957 tw_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tw_m(k,j,i) 958 ENDDO 959 ENDIF 960 ENDIF 961 962 ! 963 !-- Tendency terms for potential temperature 964 tend(:,j,i) = 0.0 965 IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN 966 CALL advec_s_pw( i, j, pt ) 967 ELSE 968 CALL advec_s_up( i, j, pt ) 969 ENDIF 970 IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) & 971 THEN 972 CALL diffusion_s( i, j, ddzu, ddzw, kh_m, pt_m, shf_m, & 973 tswst_m, tend ) 974 ELSE 975 CALL diffusion_s( i, j, ddzu, ddzw, kh, pt, shf, tswst, tend ) 976 ENDIF 977 978 ! 979 !-- If required compute heating/cooling due to long wave radiation 980 !-- processes 981 IF ( radiation ) THEN 982 CALL calc_radiation( i, j ) 983 ENDIF 984 985 ! 986 !-- If required compute impact of latent heat due to precipitation 987 IF ( precipitation ) THEN 988 CALL impact_of_latent_heat( i, j ) 989 ENDIF 990 CALL user_actions( i, j, 'pt-tendency' ) 991 992 ! 993 !-- Prognostic equation for potential temperature 1027 CALL advec_s_pw( i, j, sa ) 1028 ELSE 1029 CALL advec_s_up( i, j, sa ) 1030 ENDIF 1031 CALL diffusion_s( i, j, ddzu, ddzw, kh, sa, saswsb, saswst, & 1032 tend ) 1033 1034 CALL user_actions( i, j, 'sa-tendency' ) 1035 1036 ! 1037 !-- Prognostic equation for salinity 994 1038 DO k = nzb_s_inner(j,i)+1, nzt 995 pt_p(k,j,i) = ( 1.0-tsc(1) ) * pt_m(k,j,i) + tsc(1)*pt(k,j,i) +& 996 dt_3d * ( & 997 tsc(2) * tend(k,j,i) + tsc(3) * tpt_m(k,j,i) & 998 ) - & 999 tsc(5) * rdf(k) * ( pt(k,j,i) - pt_init(k) ) 1039 sa_p(k,j,i) = tsc(1) * sa(k,j,i) + & 1040 dt_3d * ( & 1041 tsc(2) * tend(k,j,i) + tsc(3) * tsa_m(k,j,i) & 1042 ) - & 1043 tsc(5) * rdf(k) * ( sa(k,j,i) - sa_init(k) ) 1044 IF ( sa_p(k,j,i) < 0.0 ) sa_p(k,j,i) = 0.1 * sa(k,j,i) 1000 1045 ENDDO 1001 1046 … … 1005 1050 IF ( intermediate_timestep_count == 1 ) THEN 1006 1051 DO k = nzb_s_inner(j,i)+1, nzt 1007 t pt_m(k,j,i) = tend(k,j,i)1052 tsa_m(k,j,i) = tend(k,j,i) 1008 1053 ENDDO 1009 1054 ELSEIF ( intermediate_timestep_count < & 1010 1055 intermediate_timestep_count_max ) THEN 1011 1056 DO k = nzb_s_inner(j,i)+1, nzt 1012 tpt_m(k,j,i) = -9.5625 * tend(k,j,i) + & 1013 5.3125 * tpt_m(k,j,i) 1014 ENDDO 1015 ENDIF 1016 ENDIF 1017 1018 ! 1019 !-- If required, compute prognostic equation for salinity 1020 IF ( ocean ) THEN 1021 1022 ! 1023 !-- Tendency-terms for salinity 1024 tend(:,j,i) = 0.0 1025 IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) & 1026 THEN 1027 CALL advec_s_pw( i, j, sa ) 1057 tsa_m(k,j,i) = -9.5625 * tend(k,j,i) + & 1058 5.3125 * tsa_m(k,j,i) 1059 ENDDO 1060 ENDIF 1061 ENDIF 1062 1063 ! 1064 !-- Calculate density by the equation of state for seawater 1065 CALL eqn_state_seawater( i, j ) 1066 1067 ENDIF 1068 1069 ! 1070 !-- If required, compute prognostic equation for total water content / 1071 !-- scalar 1072 IF ( humidity .OR. passive_scalar ) THEN 1073 1074 ! 1075 !-- Tendency-terms for total water content / scalar 1076 tend(:,j,i) = 0.0 1077 IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) & 1078 THEN 1079 CALL advec_s_pw( i, j, q ) 1080 ELSE 1081 CALL advec_s_up( i, j, q ) 1082 ENDIF 1083 IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' )& 1084 THEN 1085 CALL diffusion_s( i, j, ddzu, ddzw, kh_m, q_m, qsws_m, & 1086 qswst_m, tend ) 1087 ELSE 1088 CALL diffusion_s( i, j, ddzu, ddzw, kh, q, qsws, qswst, & 1089 tend ) 1090 ENDIF 1091 1092 ! 1093 !-- If required compute decrease of total water content due to 1094 !-- precipitation 1095 IF ( precipitation ) THEN 1096 CALL calc_precipitation( i, j ) 1097 ENDIF 1098 CALL user_actions( i, j, 'q-tendency' ) 1099 1100 ! 1101 !-- Prognostic equation for total water content / scalar 1102 DO k = nzb_s_inner(j,i)+1, nzt 1103 q_p(k,j,i) = ( 1.0-tsc(1) ) * q_m(k,j,i) + tsc(1)*q(k,j,i) +& 1104 dt_3d * ( & 1105 tsc(2) * tend(k,j,i) + tsc(3) * tq_m(k,j,i) & 1106 ) - & 1107 tsc(5) * rdf(k) * ( q(k,j,i) - q_init(k) ) 1108 IF ( q_p(k,j,i) < 0.0 ) q_p(k,j,i) = 0.1 * q(k,j,i) 1109 ENDDO 1110 1111 ! 1112 !-- Calculate tendencies for the next Runge-Kutta step 1113 IF ( timestep_scheme(1:5) == 'runge' ) THEN 1114 IF ( intermediate_timestep_count == 1 ) THEN 1115 DO k = nzb_s_inner(j,i)+1, nzt 1116 tq_m(k,j,i) = tend(k,j,i) 1117 ENDDO 1118 ELSEIF ( intermediate_timestep_count < & 1119 intermediate_timestep_count_max ) THEN 1120 DO k = nzb_s_inner(j,i)+1, nzt 1121 tq_m(k,j,i) = -9.5625 * tend(k,j,i) + & 1122 5.3125 * tq_m(k,j,i) 1123 ENDDO 1124 ENDIF 1125 ENDIF 1126 1127 ENDIF 1128 1129 ! 1130 !-- If required, compute prognostic equation for turbulent kinetic 1131 !-- energy (TKE) 1132 IF ( .NOT. constant_diffusion ) THEN 1133 1134 ! 1135 !-- Tendency-terms for TKE 1136 tend(:,j,i) = 0.0 1137 IF ( ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) & 1138 .AND. .NOT. use_upstream_for_tke ) THEN 1139 CALL advec_s_pw( i, j, e ) 1140 ELSE 1141 CALL advec_s_up( i, j, e ) 1142 ENDIF 1143 IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' )& 1144 THEN 1145 IF ( .NOT. humidity ) THEN 1146 CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e_m, & 1147 km_m, l_grid, pt_m, pt_reference, & 1148 rif_m, tend, zu, zw ) 1028 1149 ELSE 1029 CALL advec_s_up( i, j, sa ) 1030 ENDIF 1031 CALL diffusion_s( i, j, ddzu, ddzw, kh, sa, saswsb, saswst, & 1032 tend ) 1033 1034 CALL user_actions( i, j, 'sa-tendency' ) 1035 1036 ! 1037 !-- Prognostic equation for salinity 1038 DO k = nzb_s_inner(j,i)+1, nzt 1039 sa_p(k,j,i) = tsc(1) * sa(k,j,i) + & 1040 dt_3d * ( & 1041 tsc(2) * tend(k,j,i) + tsc(3) * tsa_m(k,j,i) & 1042 ) - & 1043 tsc(5) * rdf(k) * ( sa(k,j,i) - sa_init(k) ) 1044 IF ( sa_p(k,j,i) < 0.0 ) sa_p(k,j,i) = 0.1 * sa(k,j,i) 1045 ENDDO 1046 1047 ! 1048 !-- Calculate tendencies for the next Runge-Kutta step 1049 IF ( timestep_scheme(1:5) == 'runge' ) THEN 1050 IF ( intermediate_timestep_count == 1 ) THEN 1051 DO k = nzb_s_inner(j,i)+1, nzt 1052 tsa_m(k,j,i) = tend(k,j,i) 1053 ENDDO 1054 ELSEIF ( intermediate_timestep_count < & 1055 intermediate_timestep_count_max ) THEN 1056 DO k = nzb_s_inner(j,i)+1, nzt 1057 tsa_m(k,j,i) = -9.5625 * tend(k,j,i) + & 1058 5.3125 * tsa_m(k,j,i) 1059 ENDDO 1060 ENDIF 1061 ENDIF 1062 1063 ! 1064 !-- Calculate density by the equation of state for seawater 1065 CALL eqn_state_seawater( i, j ) 1066 1067 ENDIF 1068 1069 ! 1070 !-- If required, compute prognostic equation for total water content / 1071 !-- scalar 1072 IF ( humidity .OR. passive_scalar ) THEN 1073 1074 ! 1075 !-- Tendency-terms for total water content / scalar 1076 tend(:,j,i) = 0.0 1077 IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) & 1078 THEN 1079 CALL advec_s_pw( i, j, q ) 1080 ELSE 1081 CALL advec_s_up( i, j, q ) 1082 ENDIF 1083 IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' )& 1084 THEN 1085 CALL diffusion_s( i, j, ddzu, ddzw, kh_m, q_m, qsws_m, & 1086 qswst_m, tend ) 1087 ELSE 1088 CALL diffusion_s( i, j, ddzu, ddzw, kh, q, qsws, qswst, & 1089 tend ) 1090 ENDIF 1091 1092 ! 1093 !-- If required compute decrease of total water content due to 1094 !-- precipitation 1095 IF ( precipitation ) THEN 1096 CALL calc_precipitation( i, j ) 1097 ENDIF 1098 CALL user_actions( i, j, 'q-tendency' ) 1099 1100 ! 1101 !-- Prognostic equation for total water content / scalar 1102 DO k = nzb_s_inner(j,i)+1, nzt 1103 q_p(k,j,i) = ( 1.0-tsc(1) ) * q_m(k,j,i) + tsc(1)*q(k,j,i) +& 1104 dt_3d * ( & 1105 tsc(2) * tend(k,j,i) + tsc(3) * tq_m(k,j,i) & 1106 ) - & 1107 tsc(5) * rdf(k) * ( q(k,j,i) - q_init(k) ) 1108 IF ( q_p(k,j,i) < 0.0 ) q_p(k,j,i) = 0.1 * q(k,j,i) 1109 ENDDO 1110 1111 ! 1112 !-- Calculate tendencies for the next Runge-Kutta step 1113 IF ( timestep_scheme(1:5) == 'runge' ) THEN 1114 IF ( intermediate_timestep_count == 1 ) THEN 1115 DO k = nzb_s_inner(j,i)+1, nzt 1116 tq_m(k,j,i) = tend(k,j,i) 1117 ENDDO 1118 ELSEIF ( intermediate_timestep_count < & 1119 intermediate_timestep_count_max ) THEN 1120 DO k = nzb_s_inner(j,i)+1, nzt 1121 tq_m(k,j,i) = -9.5625 * tend(k,j,i) + & 1122 5.3125 * tq_m(k,j,i) 1123 ENDDO 1124 ENDIF 1125 ENDIF 1126 1127 ENDIF 1128 1129 ! 1130 !-- If required, compute prognostic equation for turbulent kinetic 1131 !-- energy (TKE) 1132 IF ( .NOT. constant_diffusion ) THEN 1133 1134 ! 1135 !-- Tendency-terms for TKE 1136 tend(:,j,i) = 0.0 1137 IF ( ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) & 1138 .AND. .NOT. use_upstream_for_tke ) THEN 1139 CALL advec_s_pw( i, j, e ) 1140 ELSE 1141 CALL advec_s_up( i, j, e ) 1142 ENDIF 1143 IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' )& 1144 THEN 1145 IF ( .NOT. humidity ) THEN 1146 CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e_m, & 1147 km_m, l_grid, pt_m, pt_reference, & 1148 rif_m, tend, zu, zw ) 1150 CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e_m, & 1151 km_m, l_grid, vpt_m, pt_reference, & 1152 rif_m, tend, zu, zw ) 1153 ENDIF 1154 ELSE 1155 IF ( .NOT. humidity ) THEN 1156 IF ( ocean ) THEN 1157 CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, & 1158 km, l_grid, rho, prho_reference, & 1159 rif, tend, zu, zw ) 1149 1160 ELSE 1150 CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e _m,&1151 km _m, l_grid, vpt_m, pt_reference,&1152 rif_m,tend, zu, zw )1161 CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, & 1162 km, l_grid, pt, pt_reference, rif, & 1163 tend, zu, zw ) 1153 1164 ENDIF 1154 1165 ELSE 1155 IF ( .NOT. humidity ) THEN 1156 IF ( ocean ) THEN 1157 CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, & 1158 km, l_grid, rho, prho_reference, & 1159 rif, tend, zu, zw ) 1160 ELSE 1161 CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, & 1162 km, l_grid, pt, pt_reference, rif, & 1163 tend, zu, zw ) 1164 ENDIF 1165 ELSE 1166 CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, km, & 1167 l_grid, vpt, pt_reference, rif, tend, & 1168 zu, zw ) 1169 ENDIF 1170 ENDIF 1171 CALL production_e( i, j ) 1172 CALL user_actions( i, j, 'e-tendency' ) 1173 1174 ! 1175 !-- Prognostic equation for TKE. 1176 !-- Eliminate negative TKE values, which can occur due to numerical 1177 !-- reasons in the course of the integration. In such cases the old 1178 !-- TKE value is reduced by 90%. 1179 DO k = nzb_s_inner(j,i)+1, nzt 1180 e_p(k,j,i) = ( 1.0-tsc(1) ) * e_m(k,j,i) + tsc(1)*e(k,j,i) +& 1181 dt_3d * ( & 1182 tsc(2) * tend(k,j,i) + tsc(3) * te_m(k,j,i) & 1183 ) 1184 IF ( e_p(k,j,i) < 0.0 ) e_p(k,j,i) = 0.1 * e(k,j,i) 1185 ENDDO 1186 1187 ! 1188 !-- Calculate tendencies for the next Runge-Kutta step 1189 IF ( timestep_scheme(1:5) == 'runge' ) THEN 1190 IF ( intermediate_timestep_count == 1 ) THEN 1191 DO k = nzb_s_inner(j,i)+1, nzt 1192 te_m(k,j,i) = tend(k,j,i) 1193 ENDDO 1194 ELSEIF ( intermediate_timestep_count < & 1195 intermediate_timestep_count_max ) THEN 1196 DO k = nzb_s_inner(j,i)+1, nzt 1197 te_m(k,j,i) = -9.5625 * tend(k,j,i) + & 1198 5.3125 * te_m(k,j,i) 1199 ENDDO 1200 ENDIF 1201 ENDIF 1202 1203 ENDIF ! TKE equation 1204 1205 ENDIF ! Gridpoints excluding the non-cyclic wall 1166 CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, km, & 1167 l_grid, vpt, pt_reference, rif, tend, & 1168 zu, zw ) 1169 ENDIF 1170 ENDIF 1171 CALL production_e( i, j ) 1172 CALL user_actions( i, j, 'e-tendency' ) 1173 1174 ! 1175 !-- Prognostic equation for TKE. 1176 !-- Eliminate negative TKE values, which can occur due to numerical 1177 !-- reasons in the course of the integration. In such cases the old 1178 !-- TKE value is reduced by 90%. 1179 DO k = nzb_s_inner(j,i)+1, nzt 1180 e_p(k,j,i) = ( 1.0-tsc(1) ) * e_m(k,j,i) + tsc(1)*e(k,j,i) +& 1181 dt_3d * ( & 1182 tsc(2) * tend(k,j,i) + tsc(3) * te_m(k,j,i) & 1183 ) 1184 IF ( e_p(k,j,i) < 0.0 ) e_p(k,j,i) = 0.1 * e(k,j,i) 1185 ENDDO 1186 1187 ! 1188 !-- Calculate tendencies for the next Runge-Kutta step 1189 IF ( timestep_scheme(1:5) == 'runge' ) THEN 1190 IF ( intermediate_timestep_count == 1 ) THEN 1191 DO k = nzb_s_inner(j,i)+1, nzt 1192 te_m(k,j,i) = tend(k,j,i) 1193 ENDDO 1194 ELSEIF ( intermediate_timestep_count < & 1195 intermediate_timestep_count_max ) THEN 1196 DO k = nzb_s_inner(j,i)+1, nzt 1197 te_m(k,j,i) = -9.5625 * tend(k,j,i) + & 1198 5.3125 * te_m(k,j,i) 1199 ENDDO 1200 ENDIF 1201 ENDIF 1202 1203 ENDIF ! TKE equation 1206 1204 1207 1205 ENDDO … … 1268 1266 ! 1269 1267 !-- Prognostic equation for u-velocity component 1270 DO i = nxl , nxr1268 DO i = nxlu, nxr 1271 1269 DO j = nys, nyn 1272 1270 DO k = nzb_u_inner(j,i)+1, nzt … … 1285 1283 IF ( timestep_scheme(1:5) == 'runge' ) THEN 1286 1284 IF ( intermediate_timestep_count == 1 ) THEN 1287 DO i = nxl , nxr1285 DO i = nxlu, nxr 1288 1286 DO j = nys, nyn 1289 1287 DO k = nzb_u_inner(j,i)+1, nzt … … 1294 1292 ELSEIF ( intermediate_timestep_count < & 1295 1293 intermediate_timestep_count_max ) THEN 1296 DO i = nxl , nxr1294 DO i = nxlu, nxr 1297 1295 DO j = nys, nyn 1298 1296 DO k = nzb_u_inner(j,i)+1, nzt … … 1340 1338 !-- Prognostic equation for v-velocity component 1341 1339 DO i = nxl, nxr 1342 DO j = nys , nyn1340 DO j = nysv, nyn 1343 1341 DO k = nzb_v_inner(j,i)+1, nzt 1344 1342 v_p(k,j,i) = ( 1.0-tsc(1) ) * v_m(k,j,i) + tsc(1) * v(k,j,i) + & … … 1357 1355 IF ( intermediate_timestep_count == 1 ) THEN 1358 1356 DO i = nxl, nxr 1359 DO j = nys , nyn1357 DO j = nysv, nyn 1360 1358 DO k = nzb_v_inner(j,i)+1, nzt 1361 1359 tv_m(k,j,i) = tend(k,j,i) … … 1366 1364 intermediate_timestep_count_max ) THEN 1367 1365 DO i = nxl, nxr 1368 DO j = nys , nyn1366 DO j = nysv, nyn 1369 1367 DO k = nzb_v_inner(j,i)+1, nzt 1370 1368 tv_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tv_m(k,j,i) -
palm/trunk/SOURCE/time_integration.f90
r102 r106 4 4 ! Actual revisions: 5 5 ! ----------------- 6 ! Call of new routine surface_coupler 6 ! Call of new routine surface_coupler, 7 ! presure solver is called after the first Runge-Kutta substep instead of the 8 ! last in case that call_psolver_at_all_substeps = .F.; for this case, the 9 ! random perturbation has to be added to the velocity fields also after the 10 ! first substep 7 11 ! 8 12 ! Former revisions: … … 194 198 ! 195 199 !-- Impose a random perturbation on the horizontal velocity field 196 IF ( create_disturbances .AND. & 200 IF ( create_disturbances .AND. & 201 ( call_psolver_at_all_substeps .AND. & 197 202 intermediate_timestep_count == intermediate_timestep_count_max )& 203 .OR. ( .NOT. call_psolver_at_all_substeps .AND. & 204 intermediate_timestep_count == 1 ) ) & 198 205 THEN 199 206 time_disturb = time_disturb + dt_3d … … 218 225 !-- Reduce the velocity divergence via the equation for perturbation 219 226 !-- pressure. 220 IF ( intermediate_timestep_count == intermediate_timestep_count_max&221 .OR.call_psolver_at_all_substeps ) THEN227 IF ( intermediate_timestep_count == 1 .OR. & 228 call_psolver_at_all_substeps ) THEN 222 229 CALL pres 223 230 ENDIF 224 225 !226 !-- In case of a non-cyclic lateral wall, set the boundary conditions for227 !-- the velocities at the outflow228 ! IF ( bc_lr /= 'cyclic' .OR. bc_ns /= 'cyclic' ) THEN229 ! CALL boundary_conds( 'outflow_uvw' )230 ! ENDIF231 231 232 232 !
Note: See TracChangeset
for help on using the changeset viewer.
