Changeset 4677 for palm/trunk/SOURCE/model_1d_mod.f90
- Timestamp:
- Sep 14, 2020 7:55:28 AM (4 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/model_1d_mod.f90 (modified) (47 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/model_1d_mod.f90
r4586 r4677 1 1 !> @file model_1d_mod.f90 2 !------------------------------------------------------------------------------ !2 !--------------------------------------------------------------------------------------------------! 3 3 ! This file is part of the PALM model system. 4 4 ! 5 ! PALM is free software: you can redistribute it and/or modify it under the 6 ! terms of the GNU General Public License as published by the Free Software 7 ! Foundation, either version 3 of the License, or (at your option) any later 8 ! version. 9 ! 10 ! PALM is distributed in the hope that it will be useful, but WITHOUT ANY 11 ! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR 12 ! A PARTICULAR PURPOSE. See the GNU General Public License for more details. 13 ! 14 ! You should have received a copy of the GNU General Public License along with 15 ! PALM. If not, see <http://www.gnu.org/licenses/>. 5 ! PALM is free software: you can redistribute it and/or modify it under the terms of the GNU General 6 ! Public License as published by the Free Software Foundation, either version 3 of the License, or 7 ! (at your option) any later version. 8 ! 9 ! PALM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the 10 ! implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General 11 ! Public License for more details. 12 ! 13 ! You should have received a copy of the GNU General Public License along with PALM. If not, see 14 ! <http://www.gnu.org/licenses/>. 16 15 ! 17 16 ! Copyright 1997-2020 Leibniz Universitaet Hannover 18 !------------------------------------------------------------------------------ !17 !--------------------------------------------------------------------------------------------------! 19 18 ! 20 19 ! Current revisions: … … 25 24 ! ----------------- 26 25 ! $Id$ 26 ! file re-formatted to follow the PALM coding standard 27 ! 28 ! 4586 2020-07-01 16:16:43Z gronemeier 27 29 ! renamed Richardson flux number into gradient Richardson number 28 30 ! … … 53 55 !> the 1D model uses different turbulence closure approaches at least if 54 56 !> the 3D model is set to LES-mode. 55 !------------------------------------------------------------------------------ !57 !--------------------------------------------------------------------------------------------------! 56 58 MODULE model_1d_mod 57 59 58 USE arrays_3d, & 59 ONLY: dd2zu, ddzu, ddzw, dzu, dzw, pt_init, q_init, ug, u_init, & 60 vg, v_init, zu 61 62 USE basic_constants_and_equations_mod, & 60 USE arrays_3d, & 61 ONLY: dd2zu, ddzu, ddzw, dzu, dzw, pt_init, q_init, ug, u_init, vg, v_init, zu 62 63 USE basic_constants_and_equations_mod, & 63 64 ONLY: g, kappa, pi 64 65 65 USE control_parameters, & 66 ONLY: constant_diffusion, constant_flux_layer, dissipation_1d, f, & 67 humidity, ibc_e_b, intermediate_timestep_count, & 68 intermediate_timestep_count_max, km_constant, & 69 message_string, mixing_length_1d, prandtl_number, & 70 roughness_length, run_description_header, simulated_time_chr, & 71 timestep_scheme, tsc, z0h_factor 72 73 USE indices, & 66 USE control_parameters, & 67 ONLY: constant_diffusion, constant_flux_layer, dissipation_1d, f, humidity, ibc_e_b, & 68 intermediate_timestep_count, intermediate_timestep_count_max, km_constant, & 69 message_string, mixing_length_1d, prandtl_number, roughness_length, & 70 run_description_header, simulated_time_chr, timestep_scheme, tsc, z0h_factor 71 72 USE indices, & 74 73 ONLY: nzb, nzb_diff, nzt 75 74 76 75 USE kinds 77 76 78 USE pegrid, &77 USE pegrid, & 79 78 ONLY: myid 80 79 … … 176 175 ! 177 176 !-- Public variables 178 PUBLIC damp_level_1d, damp_level_ind_1d, diss1d, dt_pr_1d, & 179 dt_run_control_1d, e1d, end_time_1d, kh1d, km1d, l1d, ri1d, u1d, & 180 us1d, usws1d, v1d, vsws1d 177 PUBLIC damp_level_1d, damp_level_ind_1d, diss1d, dt_pr_1d, dt_run_control_1d, e1d, & 178 end_time_1d, kh1d, km1d, l1d, ri1d, u1d, us1d, usws1d, v1d, vsws1d 181 179 182 180 … … 185 183 SUBROUTINE init_1d_model 186 184 187 USE grid_variables, &185 USE grid_variables, & 188 186 ONLY: dx, dy 189 187 … … 196 194 ! 197 195 !-- Allocate required 1D-arrays 198 ALLOCATE( diss1d(nzb:nzt+1), diss1d_p(nzb:nzt+1), &199 e1d(nzb:nzt+1), e1d_p(nzb:nzt+1), kh1d(nzb:nzt+1), &200 km1d(nzb:nzt+1), l1d(nzb:nzt+1), l1d_init(nzb:nzt+1), &201 l1d_diss(nzb:nzt+1), ri1d(nzb:nzt+1), te_diss(nzb:nzt+1), &202 te_dissm(nzb:nzt+1), te_e(nzb:nzt+1), &203 te_em(nzb:nzt+1), te_u(nzb:nzt+1), te_um(nzb:nzt+1), &204 te_v(nzb:nzt+1), te_vm(nzb:nzt+1), u1d(nzb:nzt+1), &196 ALLOCATE( diss1d(nzb:nzt+1), diss1d_p(nzb:nzt+1), & 197 e1d(nzb:nzt+1), e1d_p(nzb:nzt+1), kh1d(nzb:nzt+1), & 198 km1d(nzb:nzt+1), l1d(nzb:nzt+1), l1d_init(nzb:nzt+1), & 199 l1d_diss(nzb:nzt+1), ri1d(nzb:nzt+1), te_diss(nzb:nzt+1), & 200 te_dissm(nzb:nzt+1), te_e(nzb:nzt+1), & 201 te_em(nzb:nzt+1), te_u(nzb:nzt+1), te_um(nzb:nzt+1), & 202 te_v(nzb:nzt+1), te_vm(nzb:nzt+1), u1d(nzb:nzt+1), & 205 203 u1d_p(nzb:nzt+1), v1d(nzb:nzt+1), v1d_p(nzb:nzt+1) ) 206 204 … … 223 221 !-- Blackadar mixing length 224 222 IF ( f /= 0.0_wp ) THEN 225 lambda = 2.7E-4_wp * SQRT( ug(nzt+1)**2 + vg(nzt+1)**2 ) / & 226 ABS( f ) + 1E-10_wp 223 lambda = 2.7E-4_wp * SQRT( ug(nzt+1)**2 + vg(nzt+1)**2 ) / ABS( f ) + 1E-10_wp 227 224 ELSE 228 225 lambda = 30.0_wp … … 258 255 259 256 ! 260 !-- Set initial horizontal velocities at the lowest grid levels to a very small 261 !-- value in order to avoid too small time steps caused by the diffusion limit262 !-- in the initial phase of a run (at k=1,dz/2 occurs in the limiting formula!)257 !-- Set initial horizontal velocities at the lowest grid levels to a very small value in order to 258 !-- avoid too small time steps caused by the diffusion limit in the initial phase of a run (at k=1, 259 !-- dz/2 occurs in the limiting formula!) 263 260 u1d(0:1) = 0.1_wp 264 261 u1d_p(0:1) = 0.1_wp … … 310 307 311 308 312 !------------------------------------------------------------------------------ !309 !--------------------------------------------------------------------------------------------------! 313 310 ! Description: 314 311 ! ------------ 315 312 !> Runge-Kutta time differencing scheme for the 1D-model. 316 !------------------------------------------------------------------------------ !313 !--------------------------------------------------------------------------------------------------! 317 314 318 315 SUBROUTINE time_integration_1d … … 335 332 336 333 ! 337 !-- Determine the time step at the start of a 1D-simulation and 338 !-- determine and printout quantitiesused for run control334 !-- Determine the time step at the start of a 1D-simulation and determine and printout quantities 335 !-- used for run control 339 336 dt_1d = 0.01_wp 340 337 CALL run_control_1d … … 345 342 346 343 ! 347 !-- Depending on the timestep scheme, carry out one or more intermediate 348 !-- timesteps 344 !-- Depending on the timestep scheme, carry out one or more intermediate timesteps 349 345 350 346 intermediate_timestep_count = 0 351 DO WHILE ( intermediate_timestep_count < & 352 intermediate_timestep_count_max ) 347 DO WHILE ( intermediate_timestep_count < intermediate_timestep_count_max ) 353 348 354 349 intermediate_timestep_count = intermediate_timestep_count + 1 … … 357 352 358 353 ! 359 !-- Compute all tendency terms. If a constant-flux layer is simulated, 360 !-- k starts at nzb+2. 354 !-- Compute all tendency terms. If a constant-flux layer is simulated, k starts at nzb+2. 361 355 DO k = nzb_diff, nzt 362 356 … … 365 359 ! 366 360 !-- u-component 367 te_u(k) = f * ( v1d(k) - vg(k) ) + ( &368 kmzp * ( u1d(k+1) - u1d(k) ) * ddzu(k+1) &369 - kmzm * ( u1d(k) - u1d(k-1) ) * ddzu(k) &361 te_u(k) = f * ( v1d(k) - vg(k) ) + ( & 362 kmzp * ( u1d(k+1) - u1d(k) ) * ddzu(k+1) & 363 - kmzm * ( u1d(k) - u1d(k-1) ) * ddzu(k) & 370 364 ) * ddzw(k) 371 365 ! 372 366 !-- v-component 373 te_v(k) = -f * ( u1d(k) - ug(k) ) + ( &374 kmzp * ( v1d(k+1) - v1d(k) ) * ddzu(k+1) &375 - kmzm * ( v1d(k) - v1d(k-1) ) * ddzu(k) &367 te_v(k) = -f * ( u1d(k) - ug(k) ) + ( & 368 kmzp * ( v1d(k+1) - v1d(k) ) * ddzu(k+1) & 369 - kmzm * ( v1d(k) - v1d(k-1) ) * ddzu(k) & 376 370 ) * ddzw(k) 377 371 ENDDO … … 387 381 ELSE 388 382 pt_0 = pt_init(k) * ( 1.0_wp + 0.61_wp * q_init(k) ) 389 flux = ( ( pt_init(k+1) - pt_init(k-1) ) + &390 0.61_wp * ( pt_init(k+1) * q_init(k+1) - &391 pt_init(k-1) * q_init(k-1) ) &383 flux = ( ( pt_init(k+1) - pt_init(k-1) ) + & 384 0.61_wp * ( pt_init(k+1) * q_init(k+1) - & 385 pt_init(k-1) * q_init(k-1) ) & 392 386 ) * dd2zu(k) 393 387 ENDIF 394 388 395 389 ! 396 !-- Calculate dissipation rate if no prognostic equation is used for 397 !-- dissipation rate 390 !-- Calculate dissipation rate if no prognostic equation is used for dissipation rate. 398 391 IF ( dissipation_1d == 'detering' ) THEN 399 392 diss1d(k) = c_0**3 * e1d(k) * SQRT( e1d(k) ) / l1d_diss(k) 400 393 ELSEIF ( dissipation_1d == 'as_in_3d_model' ) THEN 401 diss1d(k) = ( 0.19_wp + 0.74_wp * l1d_diss(k) / l1d_init(k) &402 )* e1d(k) * SQRT( e1d(k) ) / l1d_diss(k)394 diss1d(k) = ( 0.19_wp + 0.74_wp * l1d_diss(k) / l1d_init(k) ) & 395 * e1d(k) * SQRT( e1d(k) ) / l1d_diss(k) 403 396 ENDIF 404 397 ! 405 398 !-- TKE 406 te_e(k) = km1d(k) * ( ( ( u1d(k+1) - u1d(k-1) ) * dd2zu(k) )**2 &407 + ( ( v1d(k+1) - v1d(k-1) ) * dd2zu(k) )**2 &408 ) &409 - g / pt_0 * kh1d(k) * flux &410 + ( &411 kmzp * ( e1d(k+1) - e1d(k) ) * ddzu(k+1) &412 - kmzm * ( e1d(k) - e1d(k-1) ) * ddzu(k) &413 ) * ddzw(k) / sig_e &399 te_e(k) = km1d(k) * ( ( ( u1d(k+1) - u1d(k-1) ) * dd2zu(k) )**2 & 400 + ( ( v1d(k+1) - v1d(k-1) ) * dd2zu(k) )**2 & 401 ) & 402 - g / pt_0 * kh1d(k) * flux & 403 + ( & 404 kmzp * ( e1d(k+1) - e1d(k) ) * ddzu(k+1) & 405 - kmzm * ( e1d(k) - e1d(k-1) ) * ddzu(k) & 406 ) * ddzw(k) / sig_e & 414 407 - diss1d(k) 415 408 … … 420 413 alpha_buoyancy = 1.0_wp - l1d(k) / lambda 421 414 ELSE 422 alpha_buoyancy = 1.0_wp - ( 1.0_wp + ( c_2 - 1.0_wp ) &423 / ( c_2 - c_1 ) ) &415 alpha_buoyancy = 1.0_wp - ( 1.0_wp + ( c_2 - 1.0_wp ) & 416 / ( c_2 - c_1 ) ) & 424 417 * l1d(k) / lambda 425 418 ENDIF 426 419 c_3 = ( c_1 - c_2 ) * alpha_buoyancy + 1.0_wp 427 te_diss(k) = ( km1d(k) * &428 ( ( ( u1d(k+1) - u1d(k-1) ) * dd2zu(k) )**2 &429 + ( ( v1d(k+1) - v1d(k-1) ) * dd2zu(k) )**2 &430 ) * ( c_1 + (c_2 - c_1) * l1d(k) / lambda ) &431 - g / pt_0 * kh1d(k) * flux * c_3 &432 - c_2 * diss1d(k) &433 ) * diss1d(k) / ( e1d(k) + 1.0E-20_wp ) &434 + ( kmzp * ( diss1d(k+1) - diss1d(k) ) &435 * ddzu(k+1) &436 - kmzm * ( diss1d(k) - diss1d(k-1) ) &437 * ddzu(k) &420 te_diss(k) = ( km1d(k) * & 421 ( ( ( u1d(k+1) - u1d(k-1) ) * dd2zu(k) )**2 & 422 + ( ( v1d(k+1) - v1d(k-1) ) * dd2zu(k) )**2 & 423 ) * ( c_1 + (c_2 - c_1) * l1d(k) / lambda ) & 424 - g / pt_0 * kh1d(k) * flux * c_3 & 425 - c_2 * diss1d(k) & 426 ) * diss1d(k) / ( e1d(k) + 1.0E-20_wp ) & 427 + ( kmzp * ( diss1d(k+1) - diss1d(k) ) & 428 * ddzu(k+1) & 429 - kmzm * ( diss1d(k) - diss1d(k-1) ) & 430 * ddzu(k) & 438 431 ) * ddzw(k) / sig_diss 439 432 … … 445 438 ! 446 439 !-- Tendency terms at the top of the constant-flux layer. 447 !-- Finite differences of the momentum fluxes are computed using half the 448 !-- normal grid length(2.0*ddzw(k)) for the sake of enhanced accuracy440 !-- Finite differences of the momentum fluxes are computed using half the normal grid length 441 !-- (2.0*ddzw(k)) for the sake of enhanced accuracy 449 442 IF ( constant_flux_layer ) THEN 450 443 … … 457 450 ELSE 458 451 pt_0 = pt_init(k) * ( 1.0_wp + 0.61_wp * q_init(k) ) 459 flux = ( ( pt_init(k+1) - pt_init(k-1) ) + &460 0.61_wp * ( pt_init(k+1) * q_init(k+1) - &461 pt_init(k-1) * q_init(k-1) ) &452 flux = ( ( pt_init(k+1) - pt_init(k-1) ) + & 453 0.61_wp * ( pt_init(k+1) * q_init(k+1) - & 454 pt_init(k-1) * q_init(k-1) ) & 462 455 ) * dd2zu(k) 463 456 ENDIF 464 457 465 458 ! 466 !-- Calculate dissipation rate if no prognostic equation is used for 467 !-- dissipation rate 459 !-- Calculate dissipation rate if no prognostic equation is used for dissipation rate. 468 460 IF ( dissipation_1d == 'detering' ) THEN 469 461 diss1d(k) = c_0**3 * e1d(k) * SQRT( e1d(k) ) / l1d_diss(k) 470 462 ELSEIF ( dissipation_1d == 'as_in_3d_model' ) THEN 471 diss1d(k) = ( 0.19_wp + 0.74_wp * l1d_diss(k) / l1d_init(k) ) &463 diss1d(k) = ( 0.19_wp + 0.74_wp * l1d_diss(k) / l1d_init(k) ) & 472 464 * e1d(k) * SQRT( e1d(k) ) / l1d_diss(k) 473 465 ENDIF … … 475 467 ! 476 468 !-- u-component 477 te_u(k) = f * ( v1d(k) - vg(k) ) + ( &478 kmzp * ( u1d(k+1) - u1d(k) ) * ddzu(k+1) + usws1d &469 te_u(k) = f * ( v1d(k) - vg(k) ) + ( & 470 kmzp * ( u1d(k+1) - u1d(k) ) * ddzu(k+1) + usws1d & 479 471 ) * 2.0_wp * ddzw(k) 480 472 ! 481 473 !-- v-component 482 te_v(k) = -f * ( u1d(k) - ug(k) ) + ( &483 kmzp * ( v1d(k+1) - v1d(k) ) * ddzu(k+1) + vsws1d &474 te_v(k) = -f * ( u1d(k) - ug(k) ) + ( & 475 kmzp * ( v1d(k+1) - v1d(k) ) * ddzu(k+1) + vsws1d & 484 476 ) * 2.0_wp * ddzw(k) 485 477 ! … … 490 482 !> while for u and v it is not? 491 483 !> 2018-04-23, gronemeier 492 te_e(k) = km1d(k) * ( ( ( u1d(k+1) - u1d(k-1) ) * dd2zu(k) )**2 &493 + ( ( v1d(k+1) - v1d(k-1) ) * dd2zu(k) )**2 &494 ) &495 - g / pt_0 * kh1d(k) * flux &496 + ( &497 kmzp * ( e1d(k+1) - e1d(k) ) * ddzu(k+1) &498 - kmzm * ( e1d(k) - e1d(k-1) ) * ddzu(k) &499 ) * ddzw(k) / sig_e &484 te_e(k) = km1d(k) * ( ( ( u1d(k+1) - u1d(k-1) ) * dd2zu(k) )**2 & 485 + ( ( v1d(k+1) - v1d(k-1) ) * dd2zu(k) )**2 & 486 ) & 487 - g / pt_0 * kh1d(k) * flux & 488 + ( & 489 kmzp * ( e1d(k+1) - e1d(k) ) * ddzu(k+1) & 490 - kmzm * ( e1d(k) - e1d(k-1) ) * ddzu(k) & 491 ) * ddzw(k) / sig_e & 500 492 - diss1d(k) 501 493 ENDIF … … 507 499 DO k = nzb+1, nzt 508 500 509 u1d_p(k) = u1d(k) + dt_1d * ( tsc(2) * te_u(k) + & 510 tsc(3) * te_um(k) ) 511 v1d_p(k) = v1d(k) + dt_1d * ( tsc(2) * te_v(k) + & 512 tsc(3) * te_vm(k) ) 501 u1d_p(k) = u1d(k) + dt_1d * ( tsc(2) * te_u(k) + tsc(3) * te_um(k) ) 502 v1d_p(k) = v1d(k) + dt_1d * ( tsc(2) * te_v(k) + tsc(3) * te_vm(k) ) 513 503 514 504 ENDDO … … 516 506 517 507 DO k = nzb+1, nzt 518 e1d_p(k) = e1d(k) + dt_1d * ( tsc(2) * te_e(k) + & 519 tsc(3) * te_em(k) ) 508 e1d_p(k) = e1d(k) + dt_1d * ( tsc(2) * te_e(k) + tsc(3) * te_em(k) ) 520 509 ENDDO 521 510 522 511 ! 523 !-- Eliminate negative TKE values, which can result from the 524 !-- integration due to numerical inaccuracies. In such cases the TKE 525 !-- value is reduced to 10 percent of its old value. 512 !-- Eliminate negative TKE values, which can result from the integration due to numerical 513 !-- inaccuracies. In such cases the TKE value is reduced to 10 percent of its old value. 526 514 WHERE ( e1d_p < 0.0_wp ) e1d_p = 0.1_wp * e1d 527 515 528 516 IF ( dissipation_1d == 'prognostic' ) THEN 529 517 DO k = nzb+1, nzt 530 diss1d_p(k) = diss1d(k) + dt_1d * ( tsc(2) * te_diss(k) + & 531 tsc(3) * te_dissm(k) ) 518 diss1d_p(k) = diss1d(k) + dt_1d * ( tsc(2) * te_diss(k) + tsc(3) * te_dissm(k) ) 532 519 ENDDO 533 520 WHERE ( diss1d_p < 0.0_wp ) diss1d_p = 0.1_wp * diss1d … … 556 543 ENDIF 557 544 558 ELSEIF ( intermediate_timestep_count < & 559 intermediate_timestep_count_max ) THEN 545 ELSEIF ( intermediate_timestep_count < intermediate_timestep_count_max ) THEN 560 546 561 547 DO k = nzb+1, nzt … … 570 556 IF ( dissipation_1d == 'prognostic' ) THEN 571 557 DO k = nzb+1, nzt 572 te_dissm(k) = -9.5625_wp * te_diss(k) & 573 + 5.3125_wp * te_dissm(k) 558 te_dissm(k) = -9.5625_wp * te_diss(k) + 5.3125_wp * te_dissm(k) 574 559 ENDDO 575 560 ENDIF … … 581 566 ! 582 567 !-- Boundary conditions for the prognostic variables. 583 !-- At the top boundary (nzt+1) u, v, e, and diss keep their initial 584 !-- v alues (ug(nzt+1), vg(nzt+1), 0, 0).585 !-- At the bottom boundary, Dirichlet condition is used for u and v (0) 586 !-- and Neumann conditionfor e and diss (e(nzb)=e(nzb+1)).568 !-- At the top boundary (nzt+1) u, v, e, and diss keep their initial values (ug(nzt+1), 569 !-- vg(nzt+1), 0, 0). 570 !-- At the bottom boundary, Dirichlet condition is used for u and v (0) and Neumann condition 571 !-- for e and diss (e(nzb)=e(nzb+1)). 587 572 u1d_p(nzb) = 0.0_wp 588 573 v1d_p(nzb) = 0.0_wp … … 611 596 ! 612 597 !-- Stable stratification 613 ts1d = kappa * ( pt_init(nzb+1) - pt_init(nzb) ) / &614 ( LOG( zu(nzb+1) / z0h1d ) + 5.0_wp * ri1d(nzb+1) * &615 ( zu(nzb+1) - z0h1d ) / zu(nzb+1) &598 ts1d = kappa * ( pt_init(nzb+1) - pt_init(nzb) ) / & 599 ( LOG( zu(nzb+1) / z0h1d ) + 5.0_wp * ri1d(nzb+1) * & 600 ( zu(nzb+1) - z0h1d ) / zu(nzb+1) & 616 601 ) 617 602 ELSE … … 619 604 !-- Unstable stratification 620 605 a = SQRT( 1.0_wp - 16.0_wp * ri1d(nzb+1) ) 621 b = SQRT( 1.0_wp - 16.0_wp * ri1d(nzb+1) / & 622 zu(nzb+1) * z0h1d ) 623 624 ts1d = kappa * ( pt_init(nzb+1) - pt_init(nzb) ) / & 625 LOG( (a-1.0_wp) / (a+1.0_wp) * & 626 (b+1.0_wp) / (b-1.0_wp) ) 606 b = SQRT( 1.0_wp - 16.0_wp * ri1d(nzb+1) / zu(nzb+1) * z0h1d ) 607 608 ts1d = kappa * ( pt_init(nzb+1) - pt_init(nzb) ) / & 609 LOG( (a-1.0_wp) / (a+1.0_wp) * (b+1.0_wp) / (b-1.0_wp) ) 627 610 ENDIF 628 611 … … 634 617 ! 635 618 !-- Compute the gradient Richardson numbers, 636 !-- first at the top of the constant-flux layer using u* of the 637 !-- previous time step(+1E-30, if u* = 0), then in the remaining area.619 !-- first at the top of the constant-flux layer using u* of the previous time step 620 !-- (+1E-30, if u* = 0), then in the remaining area. 638 621 !-- There, the Ri numbers of the previous time step are used. 639 622 … … 646 629 flux = ts1d + 0.61_wp * pt_init(k) * qs1d 647 630 ENDIF 648 ri1d(nzb+1) = zu(nzb+1) * kappa * g * flux / & 649 ( pt_0 * ( us1d**2 + 1E-30_wp ) ) 631 ri1d(nzb+1) = zu(nzb+1) * kappa * g * flux / ( pt_0 * ( us1d**2 + 1E-30_wp ) ) 650 632 ENDIF 651 633 … … 656 638 ELSE 657 639 pt_0 = pt_init(k) * ( 1.0_wp + 0.61_wp * q_init(k) ) 658 flux = ( ( pt_init(k+1) - pt_init(k-1) ) &659 + 0.61_wp &660 * ( pt_init(k+1) * q_init(k+1) &661 - pt_init(k-1) * q_init(k-1) ) &640 flux = ( ( pt_init(k+1) - pt_init(k-1) ) & 641 + 0.61_wp & 642 * ( pt_init(k+1) * q_init(k+1) & 643 - pt_init(k-1) * q_init(k-1) ) & 662 644 ) * dd2zu(k) 663 645 ENDIF 664 646 IF ( ri1d(k) >= 0.0_wp ) THEN 665 ri1d(k) = g / pt_0 * flux / &666 ( ( ( u1d(k+1) - u1d(k-1) ) * dd2zu(k) )**2 &667 + ( ( v1d(k+1) - v1d(k-1) ) * dd2zu(k) )**2 &668 + 1E-30_wp &647 ri1d(k) = g / pt_0 * flux / & 648 ( ( ( u1d(k+1) - u1d(k-1) ) * dd2zu(k) )**2 & 649 + ( ( v1d(k+1) - v1d(k-1) ) * dd2zu(k) )**2 & 650 + 1E-30_wp & 669 651 ) 670 652 ELSE 671 ri1d(k) = g / pt_0 * flux / &672 ( ( ( u1d(k+1) - u1d(k-1) ) * dd2zu(k) )**2 &673 + ( ( v1d(k+1) - v1d(k-1) ) * dd2zu(k) )**2 &674 + 1E-30_wp &653 ri1d(k) = g / pt_0 * flux / & 654 ( ( ( u1d(k+1) - u1d(k-1) ) * dd2zu(k) )**2 & 655 + ( ( v1d(k+1) - v1d(k-1) ) * dd2zu(k) )**2 & 656 + 1E-30_wp & 675 657 ) * ( 1.0_wp - 16.0_wp * ri1d(k) )**0.25_wp 676 658 ENDIF 677 659 ENDDO 678 660 ! 679 !-- Richardson numbers must remain restricted to a realistic value 680 !-- range. It is exceeded excessively for very small velocities 681 !-- (u,v --> 0). 661 !-- Richardson numbers must remain restricted to a realistic value range. It is exceeded 662 !-- excessively for very small velocities (u,v --> 0). 682 663 WHERE ( ri1d < -5.0_wp ) ri1d = -5.0_wp 683 664 WHERE ( ri1d > 1.0_wp ) ri1d = 1.0_wp … … 691 672 ! 692 673 !-- Stable stratification 693 us1d = kappa * uv_total / ( 694 LOG( zu(nzb+1) / z01d ) + 5.0_wp * ri1d(nzb+1) *&695 ( zu(nzb+1) - z01d ) / zu(nzb+1)&674 us1d = kappa * uv_total / ( LOG( zu(nzb+1) / z01d ) & 675 + 5.0_wp * ri1d(nzb+1) * ( zu(nzb+1) - z01d ) & 676 / zu(nzb+1) & 696 677 ) 697 678 ELSE … … 699 680 !-- Unstable stratification 700 681 a = 1.0_wp / SQRT( SQRT( 1.0_wp - 16.0_wp * ri1d(nzb+1) ) ) 701 b = 1.0_wp / SQRT( SQRT( 1.0_wp - 16.0_wp * ri1d(nzb+1) / & 702 zu(nzb+1) * z01d ) ) 703 us1d = kappa * uv_total / ( & 704 LOG( (1.0_wp+b) / (1.0_wp-b) * (1.0_wp-a) / & 705 (1.0_wp+a) ) + & 706 2.0_wp * ( ATAN( b ) - ATAN( a ) ) & 682 b = 1.0_wp / SQRT( SQRT( 1.0_wp - 16.0_wp * ri1d(nzb+1) / zu(nzb+1) * z01d ) ) 683 us1d = kappa * uv_total / ( LOG( (1.0_wp+b) / (1.0_wp-b) * (1.0_wp-a) / & 684 (1.0_wp+a) ) + & 685 2.0_wp * ( ATAN( b ) - ATAN( a ) ) & 707 686 ) 708 687 ENDIF … … 714 693 715 694 ! 716 !-- Boundary condition for the turbulent kinetic energy and 717 !-- dissipation rate at the topof the constant-flux layer.718 !-- Additional Neumann condition de/dz = 0 at nzb is set to ensure 719 !-- compatibility withthe 3D model.695 !-- Boundary condition for the turbulent kinetic energy and dissipation rate at the top 696 !-- of the constant-flux layer. 697 !-- Additional Neumann condition de/dz = 0 at nzb is set to ensure compatibility with 698 !-- the 3D model. 720 699 IF ( ibc_e_b == 2 ) THEN 721 700 e1d(nzb+1) = ( us1d / c_0 )**2 … … 734 713 ! 735 714 !-- Stable stratification 736 qs1d = kappa * ( q_init(nzb+1) - q_init(nzb) ) / &737 ( LOG( zu(nzb+1) / z0h1d ) + 5.0_wp * ri1d(nzb+1) *&738 ( zu(nzb+1) - z0h1d ) / zu(nzb+1)&739 )715 qs1d = kappa * ( q_init(nzb+1) - q_init(nzb) ) / & 716 ( LOG( zu(nzb+1) / z0h1d ) + 5.0_wp * ri1d(nzb+1) * & 717 ( zu(nzb+1) - z0h1d ) / zu(nzb+1) & 718 ) 740 719 ELSE 741 720 ! 742 721 !-- Unstable stratification 743 722 a = SQRT( 1.0_wp - 16.0_wp * ri1d(nzb+1) ) 744 b = SQRT( 1.0_wp - 16.0_wp * ri1d(nzb+1) / & 745 zu(nzb+1) * z0h1d ) 746 qs1d = kappa * ( q_init(nzb+1) - q_init(nzb) ) / & 747 LOG( (a-1.0_wp) / (a+1.0_wp) * & 748 (b+1.0_wp) / (b-1.0_wp) ) 723 b = SQRT( 1.0_wp - 16.0_wp * ri1d(nzb+1) / zu(nzb+1) * z0h1d ) 724 qs1d = kappa * ( q_init(nzb+1) - q_init(nzb) ) / & 725 LOG( (a-1.0_wp) / (a+1.0_wp) * (b+1.0_wp) / (b-1.0_wp) ) 749 726 ENDIF 750 727 ELSE … … 755 732 756 733 ! 757 !-- Compute the diabatic mixing length. The unstable stratification 758 !-- must not be considered for l1d (km1d) as it is already considered 759 !-- in the dissipation of TKE via l1d_diss. Otherwise, km1d would be 760 !-- too large. 734 !-- Compute the diabatic mixing length. The unstable stratification must not be considered 735 !-- for l1d (km1d) as it is already considered in the dissipation of TKE via l1d_diss. 736 !-- Otherwise, km1d would be too large. 761 737 IF ( dissipation_1d /= 'prognostic' ) THEN 762 738 IF ( mixing_length_1d == 'blackadar' ) THEN … … 767 743 ELSE 768 744 l1d(k) = l1d_init(k) 769 l1d_diss(k) = l1d_init(k) * & 770 SQRT( 1.0_wp - 16.0_wp * ri1d(k) ) 745 l1d_diss(k) = l1d_init(k) * SQRT( 1.0_wp - 16.0_wp * ri1d(k) ) 771 746 ENDIF 772 747 ENDDO … … 775 750 dpt_dz = ( pt_init(k+1) - pt_init(k-1) ) * dd2zu(k) 776 751 IF ( dpt_dz > 0.0_wp ) THEN 777 l_stable = 0.76_wp * SQRT( e1d(k) ) &778 / SQRT( g / pt_init(k) * dpt_dz ) + 1E-5_wp752 l_stable = 0.76_wp * SQRT( e1d(k) ) & 753 / SQRT( g / pt_init(k) * dpt_dz ) + 1E-5_wp 779 754 ELSE 780 755 l_stable = l1d_init(k) … … 786 761 ELSE 787 762 DO k = nzb+1, nzt 788 l1d(k) = c_0**3 * e1d(k) * SQRT( e1d(k) ) & 789 / ( diss1d(k) + 1.0E-30_wp ) 763 l1d(k) = c_0**3 * e1d(k) * SQRT( e1d(k) ) / ( diss1d(k) + 1.0E-30_wp ) 790 764 ENDDO 791 765 ENDIF 792 766 793 767 ! 794 !-- Compute the diffusion coefficients for momentum via the 795 !-- corresponding Prandtl-layer relationship and according to 796 !-- Prandtl-Kolmogorov, respectively 768 !-- Compute the diffusion coefficients for momentum via the corresponding Prandtl-layer 769 !-- relationship and according to Prandtl-Kolmogorov, respectively 797 770 IF ( constant_flux_layer ) THEN 798 771 IF ( ri1d(nzb+1) >= 0.0_wp ) THEN 799 km1d(nzb+1) = us1d * kappa * zu(nzb+1) / &772 km1d(nzb+1) = us1d * kappa * zu(nzb+1) / & 800 773 ( 1.0_wp + 5.0_wp * ri1d(nzb+1) ) 801 774 ELSE 802 km1d(nzb+1) = us1d * kappa * zu(nzb+1) * &775 km1d(nzb+1) = us1d * kappa * zu(nzb+1) * & 803 776 ( 1.0_wp - 16.0_wp * ri1d(nzb+1) )**0.25_wp 804 777 ENDIF … … 823 796 824 797 ! 825 !-- Compute the diffusion coefficient for heat via the relationship 826 !-- kh = phim / phih * km 798 !-- Compute the diffusion coefficient for heat via the relationship kh = phim / phih * km 827 799 DO k = nzb+1, nzt 828 800 IF ( ri1d(k) >= 0.0_wp ) THEN … … 865 837 ENDDO ! time loop 866 838 ! 867 !-- Set intermediate_timestep_count back to zero. This is required e.g. for 868 !-- initial calls ofcalc_mean_profile.839 !-- Set intermediate_timestep_count back to zero. This is required e.g. for initial calls of 840 !-- calc_mean_profile. 869 841 intermediate_timestep_count = 0 870 842 … … 872 844 873 845 874 !------------------------------------------------------------------------------ !846 !--------------------------------------------------------------------------------------------------! 875 847 ! Description: 876 848 ! ------------ 877 849 !> Compute and print out quantities for run control of the 1D model. 878 !------------------------------------------------------------------------------ !850 !--------------------------------------------------------------------------------------------------! 879 851 880 852 SUBROUTINE run_control_1d … … 922 894 alpha = alpha / ( 2.0_wp * pi ) * 360.0_wp 923 895 924 WRITE ( 15, 101 ) current_timestep_number_1d, simulated_time_chr, &925 dt_1d, umax, vmax, us1d,alpha, energy896 WRITE ( 15, 101 ) current_timestep_number_1d, simulated_time_chr, dt_1d, umax, vmax, us1d, & 897 alpha, energy 926 898 ! 927 899 !-- Write buffer contents to disc immediately … … 932 904 ! 933 905 !-- formats 934 100 FORMAT (///'1D run control output:'/ &935 &'------------------------------'//&936 &'ITER. HH:MM:SS DT UMAX VMAX U* ALPHA ENERG.'/&937 &'-------------------------------------------------------------')906 100 FORMAT (///'1D run control output:'/ & 907 '------------------------------'// & 908 'ITER. HH:MM:SS DT UMAX VMAX U* ALPHA ENERG.'/ & 909 '-------------------------------------------------------------') 938 910 101 FORMAT (I7,1X,A9,1X,F6.2,2X,F6.2,1X,F6.2,1X,F6.3,2X,F5.1,2X,F7.2) 939 911 … … 943 915 944 916 945 !------------------------------------------------------------------------------ !917 !--------------------------------------------------------------------------------------------------! 946 918 ! Description: 947 919 ! ------------ 948 920 !> Compute the time step w.r.t. the diffusion criterion 949 !------------------------------------------------------------------------------ !921 !--------------------------------------------------------------------------------------------------! 950 922 951 923 SUBROUTINE timestep_1d … … 965 937 966 938 ! 967 !-- Compute the currently feasible time step according to the diffusion 968 !-- criterion. At nzb+1 the halfgrid length is used.939 !-- Compute the currently feasible time step according to the diffusion criterion. At nzb+1 the half 940 !-- grid length is used. 969 941 fac = 0.125 970 942 dt_diff = dt_max_1d … … 985 957 stop_dt_1d = .TRUE. 986 958 987 WRITE( message_string, * ) 'timestep has exceeded the lower limit&', &988 ' dt_1d = ',dt_1d,'s simulation stopped!'959 WRITE( message_string, * ) 'timestep has exceeded the lower limit&', 'dt_1d = ',dt_1d, & 960 ' s simulation stopped!' 989 961 CALL message( 'timestep_1d', 'PA0192', 1, 2, 0, 6, 0 ) 990 962 … … 995 967 996 968 997 !------------------------------------------------------------------------------ !969 !--------------------------------------------------------------------------------------------------! 998 970 ! Description: 999 971 ! ------------ 1000 972 !> List output of profiles from the 1D-model 1001 !------------------------------------------------------------------------------ !973 !--------------------------------------------------------------------------------------------------! 1002 974 1003 975 SUBROUTINE print_1d_model … … 1029 1001 WRITE ( 17, 101 ) 1030 1002 DO k = nzt+1, nzb, -1 1031 WRITE ( 17, 103) k, zu(k), u1d(k), v1d(k), pt_init(k), e1d(k), &1032 ri1d(k), km1d(k),kh1d(k), l1d(k), diss1d(k)1003 WRITE ( 17, 103) k, zu(k), u1d(k), v1d(k), pt_init(k), e1d(k), ri1d(k), km1d(k), & 1004 kh1d(k), l1d(k), diss1d(k) 1033 1005 ENDDO 1034 1006 WRITE ( 17, 101 )
Note: See TracChangeset
for help on using the changeset viewer.
