Changeset 2337 for palm/trunk/SOURCE/init_1d_model.f90
- Timestamp:
- Aug 7, 2017 8:59:53 AM (7 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/init_1d_model.f90 (modified) (18 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/init_1d_model.f90
r2334 r2337 25 25 ! ----------------- 26 26 ! $Id$ 27 ! revised calculation of mixing length 28 ! removed rounding of time step 29 ! corrected calculation of virtual potential temperature 30 ! 31 ! 2334 2017-08-04 11:57:04Z gronemeier 27 32 ! set c_m = 0.4 according to Detering and Etling (1985) 28 33 ! … … 120 125 121 126 USE model_1d, & 122 ONLY: e1d, e1d_p, kh1d, km1d, l1d, l _black, qs1d, rif1d,&127 ONLY: e1d, e1d_p, kh1d, km1d, l1d, l1d_diss, l_black, qs1d, rif1d, & 123 128 simulated_time_1d, te_e, te_em, te_u, te_um, te_v, te_vm, ts1d, & 124 129 u1d, u1d_p, us1d, usws1d, v1d, v1d_p, vsws1d, z01d, z0h1d … … 141 146 ALLOCATE( e1d(nzb:nzt+1), e1d_p(nzb:nzt+1), & 142 147 kh1d(nzb:nzt+1), km1d(nzb:nzt+1), & 143 l_black(nzb:nzt+1), l1d(nzb:nzt+1), 148 l_black(nzb:nzt+1), l1d(nzb:nzt+1), l1d_diss(nzb:nzt+1), & 144 149 rif1d(nzb:nzt+1), te_e(nzb:nzt+1), & 145 150 te_em(nzb:nzt+1), te_u(nzb:nzt+1), te_um(nzb:nzt+1), & … … 183 188 ENDIF 184 189 ENDIF 185 l1d = l_black 186 u1d = u_init 187 u1d_p = u_init 188 v1d = v_init 189 v1d_p = v_init 190 l1d = l_black 191 l1d_diss = l_black 192 u1d = u_init 193 u1d_p = u_init 194 v1d = v_init 195 v1d_p = v_init 190 196 191 197 ! … … 226 232 227 233 ! 228 !-- Integrate the 1D-model equations using the leap-frogscheme234 !-- Integrate the 1D-model equations using the Runge-Kutta scheme 229 235 CALL time_integration_1d 230 236 … … 261 267 ONLY: current_timestep_number_1d, damp_level_ind_1d, dt_1d, & 262 268 dt_pr_1d, dt_run_control_1d, e1d, e1d_p, end_time_1d, & 263 kh1d, km1d, l1d, l _black, qs1d, rif1d, simulated_time_1d,&269 kh1d, km1d, l1d, l1d_diss, l_black, qs1d, rif1d, simulated_time_1d, & 264 270 stop_dt_1d, te_e, te_em, te_u, te_um, te_v, te_vm, time_pr_1d, & 265 271 ts1d, time_run_control_1d, u1d, u1d_p, us1d, usws1d, v1d, & … … 340 346 pt_0 = pt_init(k) * ( 1.0_wp + 0.61_wp * q_init(k) ) 341 347 flux = ( ( pt_init(k+1) - pt_init(k-1) ) + & 342 0.61_wp * pt_init(k) * & 343 ( q_init(k+1) - q_init(k-1) ) ) * dd2zu(k) 348 0.61_wp * ( pt_init(k+1) * q_init(k+1) - & 349 pt_init(k-1) * q_init(k-1) ) & 350 ) * dd2zu(k) 344 351 ENDIF 345 352 346 353 IF ( dissipation_1d == 'detering' ) THEN 347 354 ! 348 !-- According to Detering, c_e= 0.064349 dissipation = 0.064_wp * e1d(k) * SQRT( e1d(k) ) / l1d(k)355 !-- According to Detering, c_e=c_m^3 356 dissipation = c_m**3 * e1d(k) * SQRT( e1d(k) ) / l1d_diss(k) 350 357 ELSEIF ( dissipation_1d == 'as_in_3d_model' ) THEN 351 dissipation = ( 0.19_wp + 0.74_wp * l1d(k) / l_grid(k) ) & 352 * e1d(k) * SQRT( e1d(k) ) / l1d(k) 358 dissipation = ( 0.19_wp & 359 + 0.74_wp * l1d_diss(k) / l_grid(k) & 360 ) * e1d(k) * SQRT( e1d(k) ) / l1d_diss(k) 353 361 ENDIF 354 362 … … 380 388 pt_0 = pt_init(k) * ( 1.0_wp + 0.61_wp * q_init(k) ) 381 389 flux = ( ( pt_init(k+1) - pt_init(k-1) ) + & 382 0.61_wp * pt_init(k) * ( q_init(k+1) - q_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 392 ) * dd2zu(k) 384 393 ENDIF … … 387 396 ! 388 397 !-- According to Detering, c_e=0.064 389 dissipation = 0.064_wp * e1d(k) * SQRT( e1d(k) ) / l1d (k)398 dissipation = 0.064_wp * e1d(k) * SQRT( e1d(k) ) / l1d_diss(k) 390 399 ELSEIF ( dissipation_1d == 'as_in_3d_model' ) THEN 391 dissipation = ( 0.19_wp + 0.74_wp * l1d (k) / l_grid(k) )&392 * e1d(k) * SQRT( e1d(k) ) / l1d (k)400 dissipation = ( 0.19_wp + 0.74_wp * l1d_diss(k) / l_grid(k) ) & 401 * e1d(k) * SQRT( e1d(k) ) / l1d_diss(k) 393 402 ENDIF 394 403 … … 513 522 b = SQRT( 1.0_wp - 16.0_wp * rif1d(nzb+1) / & 514 523 zu(nzb+1) * z0h1d ) 515 ! 516 !-- In the borderline case the formula for stable stratification 517 !-- must be applied, because otherwise a zero division would 518 !-- occur in the argument of the logarithm. 519 IF ( a == 0.0_wp .OR. b == 0.0_wp ) THEN 520 ts1d = kappa * ( pt_init(nzb+1) - pt_init(nzb) ) / & 521 ( LOG( zu(nzb+1) / z0h1d ) + & 522 5.0_wp * rif1d(nzb+1) * & 523 ( zu(nzb+1) - z0h1d ) / zu(nzb+1) & 524 ) 525 ELSE 526 ts1d = kappa * ( pt_init(nzb+1) - pt_init(nzb) ) / & 527 LOG( (a-1.0_wp) / (a+1.0_wp) * & 528 (b+1.0_wp) / (b-1.0_wp) ) 529 ENDIF 524 525 ts1d = kappa * ( pt_init(nzb+1) - pt_init(nzb) ) / & 526 LOG( (a-1.0_wp) / (a+1.0_wp) * & 527 (b+1.0_wp) / (b-1.0_wp) ) 530 528 ENDIF 531 529 … … 557 555 pt_0 = pt_init(k) * ( 1.0_wp + 0.61_wp * q_init(k) ) 558 556 flux = ( ( pt_init(k+1) - pt_init(k-1) ) & 559 + 0.61_wp * pt_init(k) & 560 * ( q_init(k+1) - q_init(k-1) ) & 557 + 0.61_wp & 558 * ( pt_init(k+1) * q_init(k+1) & 559 - pt_init(k-1) * q_init(k-1) ) & 561 560 ) * dd2zu(k) 562 561 ENDIF … … 600 599 b = 1.0_wp / SQRT( SQRT( 1.0_wp - 16.0_wp * rif1d(nzb+1) / & 601 600 zu(nzb+1) * z01d ) ) 602 ! 603 !-- In the borderline case the formula for stable stratification 604 !-- must be applied, because otherwise a zero division would 605 !-- occur in the argument of the logarithm. 606 IF ( a == 1.0_wp .OR. b == 1.0_wp ) THEN 607 us1d = kappa * uv_total / ( & 608 LOG( zu(nzb+1) / z01d ) + & 609 5.0_wp * rif1d(nzb+1) * ( zu(nzb+1) - z01d ) / & 610 zu(nzb+1) ) 611 ELSE 612 us1d = kappa * uv_total / ( & 613 LOG( (1.0_wp+b) / (1.0_wp-b) * (1.0_wp-a) / & 614 (1.0_wp+a) ) + & 615 2.0_wp * ( ATAN( b ) - ATAN( a ) ) & 616 ) 617 ENDIF 601 us1d = kappa * uv_total / ( & 602 LOG( (1.0_wp+b) / (1.0_wp-b) * (1.0_wp-a) / & 603 (1.0_wp+a) ) + & 604 2.0_wp * ( ATAN( b ) - ATAN( a ) ) & 605 ) 618 606 ENDIF 619 607 … … 649 637 b = SQRT( 1.0_wp - 16.0_wp * rif1d(nzb+1) / & 650 638 zu(nzb+1) * z0h1d ) 651 ! 652 !-- In the borderline case the formula for stable stratification 653 !-- must be applied, because otherwise a zero division would 654 !-- occur in the argument of the logarithm. 655 IF ( a == 1.0_wp .OR. b == 1.0_wp ) THEN 656 qs1d = kappa * ( q_init(nzb+1) - q_init(nzb) ) / & 657 ( LOG( zu(nzb+1) / z0h1d ) + & 658 5.0_wp * rif1d(nzb+1) * & 659 ( zu(nzb+1) - z0h1d ) / zu(nzb+1) & 660 ) 661 ELSE 662 qs1d = kappa * ( q_init(nzb+1) - q_init(nzb) ) / & 663 LOG( (a-1.0_wp) / (a+1.0_wp) * & 664 (b+1.0_wp) / (b-1.0_wp) ) 665 ENDIF 666 ENDIF 639 qs1d = kappa * ( q_init(nzb+1) - q_init(nzb) ) / & 640 LOG( (a-1.0_wp) / (a+1.0_wp) * & 641 (b+1.0_wp) / (b-1.0_wp) ) 642 ENDIF 667 643 ELSE 668 644 qs1d = 0.0_wp 669 ENDIF 645 ENDIF 670 646 671 647 ENDIF ! constant_flux_layer 672 648 673 649 ! 674 !-- Compute the diabatic mixing length 650 !-- Compute the diabatic mixing length. The unstable stratification 651 !-- must not be considered for l1d (km1d) as it is already considered 652 !-- in the dissipation of TKE via l1d_diss. Otherwise, km1d would be 653 !-- too large. 675 654 IF ( mixing_length_1d == 'blackadar' ) THEN 676 655 DO k = nzb+1, nzt 677 656 IF ( rif1d(k) >= 0.0_wp ) THEN 678 657 l1d(k) = l_black(k) / ( 1.0_wp + 5.0_wp * rif1d(k) ) 658 l1d_diss(k) = l1d(k) 679 659 ELSE 680 l1d(k) = l_black(k) * & 681 ( 1.0_wp - 16.0_wp * rif1d(k) )**0.25_wp 660 l1d(k) = l_black(k) 661 l1d_diss(k) = l_black(k) * & 662 SQRT( 1.0_wp - 16.0_wp * rif1d(k) ) 682 663 ENDIF 683 l1d(k) = l_black(k)684 664 ENDDO 685 686 665 ELSEIF ( mixing_length_1d == 'as_in_3d_model' ) THEN 687 666 DO k = nzb+1, nzt … … 694 673 ENDIF 695 674 l1d(k) = MIN( l_grid(k), l_stable ) 675 l1d_diss(k) = l1d(k) 696 676 ENDDO 697 677 ENDIF … … 700 680 !-- Compute the diffusion coefficients for momentum via the 701 681 !-- corresponding Prandtl-layer relationship and according to 702 !-- Prandtl-Kolmogorov, respectively. The unstable stratification is 703 !-- computed via the adiabatic mixing length, for the unstability has 704 !-- already been taken account of via the TKE (cf. also Diss.). 682 !-- Prandtl-Kolmogorov, respectively 705 683 IF ( constant_flux_layer ) THEN 706 684 IF ( rif1d(nzb+1) >= 0.0_wp ) THEN … … 713 691 ENDIF 714 692 DO k = nzb_diff, nzt 715 km1d(k) = c_m * SQRT( e1d(k) ) 716 717 IF ( rif1d(k) >= 0.0_wp ) THEN 718 km1d(k) = km1d(k) * l1d(k) 719 ELSE 720 km1d(k) = km1d(k) * l_black(k) 721 ENDIF 693 km1d(k) = c_m * SQRT( e1d(k) ) * l1d(k) 722 694 ENDDO 723 695 … … 882 854 883 855 USE model_1d, & 884 ONLY: dt_1d, dt_max_1d, km1d, old_dt_1d,stop_dt_1d856 ONLY: dt_1d, dt_max_1d, km1d, stop_dt_1d 885 857 886 858 USE pegrid … … 922 894 ENDIF 923 895 924 !925 !-- A more or less simple new time step value is obtained taking only the926 !-- first two significant digits927 div = 1000.0_wp928 DO WHILE ( dt_1d < div )929 div = div / 10.0_wp930 ENDDO931 dt_1d = NINT( dt_1d * 100.0_wp / div ) * div / 100.0_wp932 933 old_dt_1d = dt_1d934 935 936 896 END SUBROUTINE timestep_1d 937 897
Note: See TracChangeset
for help on using the changeset viewer.
