Changeset 4586 for palm/trunk/SOURCE/model_1d_mod.f90
- Timestamp:
- Jul 1, 2020 4:16:43 PM (4 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/model_1d_mod.f90 (modified) (21 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/model_1d_mod.f90
r4449 r4586 20 20 ! Current revisions: 21 21 ! ----------------- 22 ! 23 ! 22 ! 23 ! 24 24 ! Former revisions: 25 25 ! ----------------- 26 26 ! $Id$ 27 ! renamed Richardson flux number into gradient Richardson number 28 ! 29 ! 4449 2020-03-09 14:43:16Z suehring 27 30 ! Set intermediate_timestep_count back to zero after 1D-model integration. 28 31 ! This is required e.g. for initial calls of calc_mean_profile. … … 120 123 REAL(wp), DIMENSION(:), ALLOCATABLE :: l1d_init !< initial mixing length (1d-model) 121 124 REAL(wp), DIMENSION(:), ALLOCATABLE :: l1d_diss !< mixing length for dissipation (1d-model) 122 REAL(wp), DIMENSION(:), ALLOCATABLE :: ri f1d !< Richardson fluxnumber (1d-model)125 REAL(wp), DIMENSION(:), ALLOCATABLE :: ri1d !< gradient Richardson number (1d-model) 123 126 REAL(wp), DIMENSION(:), ALLOCATABLE :: te_diss !< tendency of diss (1d-model) 124 127 REAL(wp), DIMENSION(:), ALLOCATABLE :: te_dissm !< weighted tendency of diss for previous sub-timestep (1d-model) … … 174 177 !-- Public variables 175 178 PUBLIC damp_level_1d, damp_level_ind_1d, diss1d, dt_pr_1d, & 176 dt_run_control_1d, e1d, end_time_1d, kh1d, km1d, l1d, ri f1d, u1d,&179 dt_run_control_1d, e1d, end_time_1d, kh1d, km1d, l1d, ri1d, u1d, & 177 180 us1d, usws1d, v1d, vsws1d 178 181 … … 196 199 e1d(nzb:nzt+1), e1d_p(nzb:nzt+1), kh1d(nzb:nzt+1), & 197 200 km1d(nzb:nzt+1), l1d(nzb:nzt+1), l1d_init(nzb:nzt+1), & 198 l1d_diss(nzb:nzt+1), ri f1d(nzb:nzt+1), te_diss(nzb:nzt+1),&201 l1d_diss(nzb:nzt+1), ri1d(nzb:nzt+1), te_diss(nzb:nzt+1), & 199 202 te_dissm(nzb:nzt+1), te_e(nzb:nzt+1), & 200 203 te_em(nzb:nzt+1), te_u(nzb:nzt+1), te_um(nzb:nzt+1), & … … 211 214 e1d = 0.0_wp; e1d_p = 0.0_wp 212 215 kh1d = 0.0_wp; km1d = 0.0_wp 213 ri f1d = 0.0_wp216 ri1d = 0.0_wp 214 217 ! 215 218 !-- Compute the mixing length … … 414 417 ! 415 418 !-- dissipation rate 416 IF ( ri f1d(k) >= 0.0_wp ) THEN419 IF ( ri1d(k) >= 0.0_wp ) THEN 417 420 alpha_buoyancy = 1.0_wp - l1d(k) / lambda 418 421 ELSE … … 604 607 IF ( constant_flux_layer ) THEN 605 608 ! 606 !-- Compute theta* using Ri fnumbers of the previous time step607 IF ( ri f1d(nzb+1) >= 0.0_wp ) THEN609 !-- Compute theta* using Ri numbers of the previous time step 610 IF ( ri1d(nzb+1) >= 0.0_wp ) THEN 608 611 ! 609 612 !-- Stable stratification 610 613 ts1d = kappa * ( pt_init(nzb+1) - pt_init(nzb) ) / & 611 ( LOG( zu(nzb+1) / z0h1d ) + 5.0_wp * ri f1d(nzb+1) *&614 ( LOG( zu(nzb+1) / z0h1d ) + 5.0_wp * ri1d(nzb+1) * & 612 615 ( zu(nzb+1) - z0h1d ) / zu(nzb+1) & 613 616 ) … … 615 618 ! 616 619 !-- Unstable stratification 617 a = SQRT( 1.0_wp - 16.0_wp * ri f1d(nzb+1) )618 b = SQRT( 1.0_wp - 16.0_wp * ri f1d(nzb+1) /&620 a = SQRT( 1.0_wp - 16.0_wp * ri1d(nzb+1) ) 621 b = SQRT( 1.0_wp - 16.0_wp * ri1d(nzb+1) / & 619 622 zu(nzb+1) * z0h1d ) 620 623 … … 630 633 !> 2018-04-23, gronemeier 631 634 ! 632 !-- Compute the Richardson-fluxnumbers,635 !-- Compute the gradient Richardson numbers, 633 636 !-- first at the top of the constant-flux layer using u* of the 634 637 !-- previous time step (+1E-30, if u* = 0), then in the remaining area. 635 !-- There the rif-numbers of the previous time step are used.638 !-- There, the Ri numbers of the previous time step are used. 636 639 637 640 IF ( constant_flux_layer ) THEN … … 643 646 flux = ts1d + 0.61_wp * pt_init(k) * qs1d 644 647 ENDIF 645 ri f1d(nzb+1) = zu(nzb+1) * kappa * g * flux /&648 ri1d(nzb+1) = zu(nzb+1) * kappa * g * flux / & 646 649 ( pt_0 * ( us1d**2 + 1E-30_wp ) ) 647 650 ENDIF … … 659 662 ) * dd2zu(k) 660 663 ENDIF 661 IF ( ri f1d(k) >= 0.0_wp ) THEN662 ri f1d(k) = g / pt_0 * flux /&664 IF ( ri1d(k) >= 0.0_wp ) THEN 665 ri1d(k) = g / pt_0 * flux / & 663 666 ( ( ( u1d(k+1) - u1d(k-1) ) * dd2zu(k) )**2 & 664 667 + ( ( v1d(k+1) - v1d(k-1) ) * dd2zu(k) )**2 & … … 666 669 ) 667 670 ELSE 668 ri f1d(k) = g / pt_0 * flux /&671 ri1d(k) = g / pt_0 * flux / & 669 672 ( ( ( u1d(k+1) - u1d(k-1) ) * dd2zu(k) )**2 & 670 673 + ( ( v1d(k+1) - v1d(k-1) ) * dd2zu(k) )**2 & 671 674 + 1E-30_wp & 672 ) * ( 1.0_wp - 16.0_wp * ri f1d(k) )**0.25_wp675 ) * ( 1.0_wp - 16.0_wp * ri1d(k) )**0.25_wp 673 676 ENDIF 674 677 ENDDO 675 678 ! 676 !-- Richardson -numbers must remain restricted to a realistic value679 !-- Richardson numbers must remain restricted to a realistic value 677 680 !-- range. It is exceeded excessively for very small velocities 678 681 !-- (u,v --> 0). 679 WHERE ( ri f1d < -5.0_wp ) rif1d = -5.0_wp680 WHERE ( ri f1d > 1.0_wp ) rif1d = 1.0_wp682 WHERE ( ri1d < -5.0_wp ) ri1d = -5.0_wp 683 WHERE ( ri1d > 1.0_wp ) ri1d = 1.0_wp 681 684 682 685 ! … … 685 688 uv_total = SQRT( u1d(nzb+1)**2 + v1d(nzb+1)**2 ) 686 689 687 IF ( ri f1d(nzb+1) >= 0.0_wp ) THEN690 IF ( ri1d(nzb+1) >= 0.0_wp ) THEN 688 691 ! 689 692 !-- Stable stratification 690 693 us1d = kappa * uv_total / ( & 691 LOG( zu(nzb+1) / z01d ) + 5.0_wp * ri f1d(nzb+1) *&694 LOG( zu(nzb+1) / z01d ) + 5.0_wp * ri1d(nzb+1) * & 692 695 ( zu(nzb+1) - z01d ) / zu(nzb+1) & 693 696 ) … … 695 698 ! 696 699 !-- Unstable stratification 697 a = 1.0_wp / SQRT( SQRT( 1.0_wp - 16.0_wp * ri f1d(nzb+1) ) )698 b = 1.0_wp / SQRT( SQRT( 1.0_wp - 16.0_wp * ri f1d(nzb+1) /&700 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) / & 699 702 zu(nzb+1) * z01d ) ) 700 703 us1d = kappa * uv_total / ( & … … 728 731 ! 729 732 !-- Compute q* 730 IF ( ri f1d(nzb+1) >= 0.0_wp ) THEN733 IF ( ri1d(nzb+1) >= 0.0_wp ) THEN 731 734 ! 732 735 !-- Stable stratification 733 736 qs1d = kappa * ( q_init(nzb+1) - q_init(nzb) ) / & 734 ( LOG( zu(nzb+1) / z0h1d ) + 5.0_wp * ri f1d(nzb+1) *&737 ( LOG( zu(nzb+1) / z0h1d ) + 5.0_wp * ri1d(nzb+1) * & 735 738 ( zu(nzb+1) - z0h1d ) / zu(nzb+1) & 736 739 ) … … 738 741 ! 739 742 !-- Unstable stratification 740 a = SQRT( 1.0_wp - 16.0_wp * ri f1d(nzb+1) )741 b = SQRT( 1.0_wp - 16.0_wp * ri f1d(nzb+1) /&743 a = SQRT( 1.0_wp - 16.0_wp * ri1d(nzb+1) ) 744 b = SQRT( 1.0_wp - 16.0_wp * ri1d(nzb+1) / & 742 745 zu(nzb+1) * z0h1d ) 743 746 qs1d = kappa * ( q_init(nzb+1) - q_init(nzb) ) / & … … 759 762 IF ( mixing_length_1d == 'blackadar' ) THEN 760 763 DO k = nzb+1, nzt 761 IF ( ri f1d(k) >= 0.0_wp ) THEN762 l1d(k) = l1d_init(k) / ( 1.0_wp + 5.0_wp * ri f1d(k) )764 IF ( ri1d(k) >= 0.0_wp ) THEN 765 l1d(k) = l1d_init(k) / ( 1.0_wp + 5.0_wp * ri1d(k) ) 763 766 l1d_diss(k) = l1d(k) 764 767 ELSE 765 768 l1d(k) = l1d_init(k) 766 769 l1d_diss(k) = l1d_init(k) * & 767 SQRT( 1.0_wp - 16.0_wp * ri f1d(k) )770 SQRT( 1.0_wp - 16.0_wp * ri1d(k) ) 768 771 ENDIF 769 772 ENDDO … … 793 796 !-- Prandtl-Kolmogorov, respectively 794 797 IF ( constant_flux_layer ) THEN 795 IF ( ri f1d(nzb+1) >= 0.0_wp ) THEN798 IF ( ri1d(nzb+1) >= 0.0_wp ) THEN 796 799 km1d(nzb+1) = us1d * kappa * zu(nzb+1) / & 797 ( 1.0_wp + 5.0_wp * ri f1d(nzb+1) )800 ( 1.0_wp + 5.0_wp * ri1d(nzb+1) ) 798 801 ELSE 799 802 km1d(nzb+1) = us1d * kappa * zu(nzb+1) * & 800 ( 1.0_wp - 16.0_wp * ri f1d(nzb+1) )**0.25_wp803 ( 1.0_wp - 16.0_wp * ri1d(nzb+1) )**0.25_wp 801 804 ENDIF 802 805 ENDIF … … 823 826 !-- kh = phim / phih * km 824 827 DO k = nzb+1, nzt 825 IF ( ri f1d(k) >= 0.0_wp ) THEN828 IF ( ri1d(k) >= 0.0_wp ) THEN 826 829 kh1d(k) = km1d(k) 827 830 ELSE 828 kh1d(k) = km1d(k) * ( 1.0_wp - 16.0_wp * ri f1d(k) )**0.25_wp831 kh1d(k) = km1d(k) * ( 1.0_wp - 16.0_wp * ri1d(k) )**0.25_wp 829 832 ENDIF 830 833 ENDDO … … 1026 1029 WRITE ( 17, 101 ) 1027 1030 DO k = nzt+1, nzb, -1 1028 WRITE ( 17, 103) k, zu(k), u1d(k), v1d(k), pt_init(k), e1d(k), &1029 ri f1d(k), km1d(k), kh1d(k), l1d(k), diss1d(k)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) 1030 1033 ENDDO 1031 1034 WRITE ( 17, 101 ) … … 1045 1048 101 FORMAT ('#',111('-')) 1046 1049 102 FORMAT ('# k zu u v pt e ', & 1047 ' rifKm Kh l diss ')1050 'Ri Km Kh l diss ') 1048 1051 103 FORMAT (1X,I4,1X,F7.1,9(1X,E10.3)) 1049 1052
Note: See TracChangeset
for help on using the changeset viewer.
