Changeset 3083 for palm/trunk/SOURCE/model_1d_mod.f90
- Timestamp:
- Jun 19, 2018 2:03:12 PM (6 years ago)
- Location:
- palm/trunk
- Files:
-
- 3 edited
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk
- Property svn:mergeinfo changed
/palm/branches/rans merged: 2919,2922,2928-2929,2960,2962,2966,2976,2982,2987-2988,2991,3008-3009,3023,3047,3050,3059,3062,3071-3082
- Property svn:mergeinfo changed
-
palm/trunk/SOURCE
- Property svn:mergeinfo changed
/palm/branches/rans/SOURCE merged: 2919,2922,2928-2929,2960,2962,2966,2976,2982,2987-2988,2991,3008-3009,3023,3047,3050,3059,3062,3071-3082
- Property svn:mergeinfo changed
-
palm/trunk/SOURCE/model_1d_mod.f90
r3049 r3083 25 25 ! ----------------- 26 26 ! $Id$ 27 ! Bugfixes: 28 ! - preset te_diss and te_e to avoid runtime errors 29 ! - implementation of buoyancy term to dissipation 30 ! according to Sogachev et al. (2012) 31 ! - where diss_p < 0 set diss_p = 0.1 diss 32 ! - calculate progn eq(diss) starting from nzb+1 33 ! Changes: 34 ! - add sig_e to TKE equation 35 ! - adjust prognostic equation of diss 36 ! - set model constants according to Koblitz (2013) 37 ! - renamed c_m to c_0 38 ! - rename l_black into l1d_init 39 ! - calculate l_grid within init_1d_model and save it as l1d_init 40 ! - calculate l1d according to DE85 if dissipation is a prognostic value 41 ! - made annotations doxygen-readable 42 ! 43 ! 3049 2018-05-29 13:52:36Z Giersch 27 44 ! Error messages revised 28 45 ! … … 164 181 USE kinds 165 182 166 USE pegrid 183 USE pegrid, & 184 ONLY: myid 167 185 168 186 … … 175 193 LOGICAL :: stop_dt_1d = .FALSE. !< termination flag, used in case of too small timestep (1d-model) 176 194 177 REAL(wp) :: c_1 = 1.44_wp !< model constant178 REAL(wp) :: c_ 2 = 1.92_wp !< model constant179 REAL(wp) :: c_ 3 = 1.44_wp !< model constant180 REAL(wp) :: c_ h = 0.0015_wp !< model constant according to Detering and Etling (1985)181 REAL(wp) :: c_ m = 0.4_wp !< model constant, 0.4 according to Detering and Etling (1985)182 REAL(wp) :: c_mu = 0.09_wp!< model constant195 REAL(wp) :: alpha_buoyancy !< model constant according to Koblitz (2013) 196 REAL(wp) :: c_0 = 0.03_wp**0.25_wp !< model constant according to Koblitz (2013) 197 REAL(wp) :: c_1 = 1.52_wp !< model constant according to Koblitz (2013) 198 REAL(wp) :: c_2 = 1.83_wp !< model constant according to Koblitz (2013) 199 REAL(wp) :: c_3 !< model constant 200 REAL(wp) :: c_mu !< model constant 183 201 REAL(wp) :: damp_level_1d = -1.0_wp !< namelist parameter 184 202 REAL(wp) :: dt_1d = 60.0_wp !< dynamic timestep (1d-model) … … 187 205 REAL(wp) :: dt_run_control_1d = 60.0_wp !< namelist parameter 188 206 REAL(wp) :: end_time_1d = 864000.0_wp !< namelist parameter 207 REAL(wp) :: lambda !< maximum mixing length 189 208 REAL(wp) :: qs1d !< characteristic humidity scale (1d-model) 190 209 REAL(wp) :: simulated_time_1d = 0.0_wp !< updated simulated time (1d-model) 191 REAL(wp) :: sig_diss = 1.3_wp !< model constant 210 REAL(wp) :: sig_diss = 2.95_wp !< model constant according to Koblitz (2013) 211 REAL(wp) :: sig_e = 2.95_wp !< model constant according to Koblitz (2013) 192 212 REAL(wp) :: time_pr_1d = 0.0_wp !< updated simulated time for profile output (1d-model) 193 213 REAL(wp) :: time_run_control_1d = 0.0_wp !< updated simulated time for run-control output (1d-model) … … 198 218 REAL(wp) :: z01d !< roughness length for momentum (1d-model) 199 219 REAL(wp) :: z0h1d !< roughness length for scalars (1d-model) 200 201 220 202 221 REAL(wp), DIMENSION(:), ALLOCATABLE :: diss1d !< tke dissipation rate (1d-model) … … 279 298 280 299 INTEGER(iwp) :: k !< loop index 281 282 REAL(wp) :: lambda !< maximum mixing length283 300 284 301 ! … … 324 341 ! 325 342 !-- Use the same mixing length as in 3D model (LES-mode) 326 ! @todo:rename (delete?) this option327 ! As the mixing length is different between RANS and LES mode, it328 ! must be distinguished here between these modes. For RANS mode,329 ! the mixing length is calculated accoding to Blackadar, which is330 ! the other option at this point.331 ! Maybe delete this option entirely (not appropriate in LES case)332 ! 2018-03-20, gronemeier343 !> @todo rename (delete?) this option 344 !> As the mixing length is different between RANS and LES mode, it 345 !> must be distinguished here between these modes. For RANS mode, 346 !> the mixing length is calculated accoding to Blackadar, which is 347 !> the other option at this point. 348 !> Maybe delete this option entirely (not appropriate in LES case) 349 !> 2018-03-20, gronemeier 333 350 DO k = nzb+1, nzt 334 351 l1d_init(k) = ( dx * dy * dzw(k) )**0.33333333333333_wp … … 359 376 us1d = 0.1_wp ! without initial friction the flow would not change 360 377 ELSE 361 diss1d(nzb+1) = 1.0_wp378 diss1d(nzb+1) = 0.001_wp 362 379 e1d(nzb+1) = 1.0_wp 363 380 km1d(nzb+1) = 1.0_wp … … 372 389 373 390 ! 374 !-- Tendencies must be preset in order to avoid runtime errors within the375 !-- first Runge-Kutta step391 !-- Tendencies must be preset in order to avoid runtime errors 392 te_diss = 0.0_wp 376 393 te_dissm = 0.0_wp 394 te_e = 0.0_wp 377 395 te_em = 0.0_wp 378 396 te_um = 0.0_wp … … 381 399 ! 382 400 !-- Set model constant 383 IF ( dissipation_1d == 'as_in_3d_model' ) c_m = 0.1_wp 401 IF ( dissipation_1d == 'as_in_3d_model' ) c_0 = 0.1_wp 402 c_mu = c_0**4 384 403 385 404 ! … … 423 442 !-- Determine the time step at the start of a 1D-simulation and 424 443 !-- determine and printout quantities used for run control 425 dt_1d = 1.0_wp444 dt_1d = 0.01_wp 426 445 CALL run_control_1d 427 446 … … 483 502 !-- dissipation rate 484 503 IF ( dissipation_1d == 'detering' ) THEN 485 diss1d(k) = c_ m**3 * e1d(k) * SQRT( e1d(k) ) / l1d_diss(k)504 diss1d(k) = c_0**3 * e1d(k) * SQRT( e1d(k) ) / l1d_diss(k) 486 505 ELSEIF ( dissipation_1d == 'as_in_3d_model' ) THEN 487 506 diss1d(k) = ( 0.19_wp + 0.74_wp * l1d_diss(k) / l1d_init(k) & … … 497 516 kmzp * ( e1d(k+1) - e1d(k) ) * ddzu(k+1) & 498 517 - kmzm * ( e1d(k) - e1d(k-1) ) * ddzu(k) & 499 ) * ddzw(k) 518 ) * ddzw(k) / sig_e & 500 519 - diss1d(k) 501 520 … … 503 522 ! 504 523 !-- dissipation rate 505 te_diss(k) = km1d(k) * & 524 IF ( rif1d(k) >= 0.0_wp ) THEN 525 alpha_buoyancy = 1.0_wp - l1d(k) / lambda 526 ELSE 527 alpha_buoyancy = 1.0_wp - ( 1.0_wp + ( c_2 - 1.0_wp ) & 528 / ( c_2 - c_1 ) ) & 529 * l1d(k) / lambda 530 ENDIF 531 c_3 = ( c_1 - c_2 ) * alpha_buoyancy 532 te_diss(k) = ( km1d(k) * & 506 533 ( ( ( u1d(k+1) - u1d(k-1) ) * dd2zu(k) )**2 & 507 534 + ( ( v1d(k+1) - v1d(k-1) ) * dd2zu(k) )**2 & 508 ) * c_1 * c_mu**0.75 / c_h * f / us1d & 509 * SQRT(e1d(k)) & 535 ) * ( c_1 + (c_2 - c_1) * l1d(k) / lambda ) & 510 536 - g / pt_0 * kh1d(k) * flux * c_3 & 511 * diss1d(k) / ( e1d(k) + 1.0E-20_wp ) & 512 + ( kmzp * ( diss1d(k+1) - diss1d(k) ) & 537 - c_2 * diss1d(k) & 538 ) * diss1d(k) / ( e1d(k) + 1.0E-20_wp ) & 539 + ( kmzp * ( diss1d(k+1) - diss1d(k) ) & 513 540 * ddzu(k+1) & 514 541 - kmzm * ( diss1d(k) - diss1d(k-1) ) & 515 542 * ddzu(k) & 516 ) * ddzw(k) / sig_diss & 517 - c_2 * diss1d(k)**2 / ( e1d(k) + 1.0E-20_wp ) 543 ) * ddzw(k) / sig_diss 518 544 519 545 ENDIF … … 546 572 !-- dissipation rate 547 573 IF ( dissipation_1d == 'detering' ) THEN 548 diss1d(k) = c_ m**3 * e1d(k) * SQRT( e1d(k) ) / l1d_diss(k)574 diss1d(k) = c_0**3 * e1d(k) * SQRT( e1d(k) ) / l1d_diss(k) 549 575 ELSEIF ( dissipation_1d == 'as_in_3d_model' ) THEN 550 576 diss1d(k) = ( 0.19_wp + 0.74_wp * l1d_diss(k) / l1d_init(k) ) & … … 565 591 !-- TKE 566 592 IF ( .NOT. dissipation_1d == 'prognostic' ) THEN 593 !> @query why integrate over 2dz 594 !> Why is it allowed to integrate over two delta-z for e 595 !> while for u and v it is not? 596 !> 2018-04-23, gronemeier 567 597 te_e(k) = km1d(k) * ( ( ( u1d(k+1) - u1d(k-1) ) * dd2zu(k) )**2& 568 598 + ( ( v1d(k+1) - v1d(k-1) ) * dd2zu(k) )**2& … … 572 602 kmzp * ( e1d(k+1) - e1d(k) ) * ddzu(k+1) & 573 603 - kmzm * ( e1d(k) - e1d(k-1) ) * ddzu(k) & 574 ) * ddzw(k) 604 ) * ddzw(k) / sig_e & 575 605 - diss1d(k) 576 606 ENDIF … … 595 625 ENDDO 596 626 597 IF ( dissipation_1d == 'prognostic' ) THEN598 DO k = nzb_diff, nzt599 diss1d_p(k) = diss1d(k) + dt_1d * ( tsc(2) * te_diss(k) + &600 tsc(3) * te_dissm(k) )601 ENDDO602 ENDIF603 627 ! 604 628 !-- Eliminate negative TKE values, which can result from the … … 606 630 !-- value is reduced to 10 percent of its old value. 607 631 WHERE ( e1d_p < 0.0_wp ) e1d_p = 0.1_wp * e1d 632 633 IF ( dissipation_1d == 'prognostic' ) THEN 634 DO k = nzb+1, nzt 635 diss1d_p(k) = diss1d(k) + dt_1d * ( tsc(2) * te_diss(k) + & 636 tsc(3) * te_dissm(k) ) 637 ENDDO 638 WHERE ( diss1d_p < 0.0_wp ) diss1d_p = 0.1_wp * diss1d 639 ENDIF 608 640 ENDIF 609 641 … … 643 675 IF ( dissipation_1d == 'prognostic' ) THEN 644 676 DO k = nzb+1, nzt 645 te_dissm(k) = -9.5625_wp * te_diss(k) + 5.3125_wp * te_dissm(k) 677 te_dissm(k) = -9.5625_wp * te_diss(k) & 678 + 5.3125_wp * te_dissm(k) 646 679 ENDDO 647 680 ENDIF … … 650 683 ENDIF 651 684 ENDIF 652 653 685 654 686 ! … … 701 733 702 734 ENDIF ! constant_flux_layer 703 735 !> @todo combine if clauses 736 !> The previous and following if clauses can be combined into a 737 !> single clause 738 !> 2018-04-23, gronemeier 704 739 ! 705 740 !-- Compute the Richardson-flux numbers, … … 789 824 !-- compatibility with the 3D model. 790 825 IF ( ibc_e_b == 2 ) THEN 791 e1d(nzb+1) = ( us1d / c_ m)**2826 e1d(nzb+1) = ( us1d / c_0 )**2 792 827 ENDIF 793 828 IF ( dissipation_1d == 'prognostic' ) THEN 794 e1d(nzb+1) = us1d**2 / SQRT( c_mu )829 e1d(nzb+1) = ( us1d / c_0 )**2 795 830 diss1d(nzb+1) = us1d**3 / ( kappa * zu(nzb+1) ) 796 831 diss1d(nzb) = diss1d(nzb+1) … … 829 864 !-- in the dissipation of TKE via l1d_diss. Otherwise, km1d would be 830 865 !-- too large. 831 IF ( mixing_length_1d == 'blackadar' ) THEN 866 IF ( dissipation_1d /= 'prognostic' ) THEN 867 IF ( mixing_length_1d == 'blackadar' ) THEN 868 DO k = nzb+1, nzt 869 IF ( rif1d(k) >= 0.0_wp ) THEN 870 l1d(k) = l1d_init(k) / ( 1.0_wp + 5.0_wp * rif1d(k) ) 871 l1d_diss(k) = l1d(k) 872 ELSE 873 l1d(k) = l1d_init(k) 874 l1d_diss(k) = l1d_init(k) * & 875 SQRT( 1.0_wp - 16.0_wp * rif1d(k) ) 876 ENDIF 877 ENDDO 878 ELSEIF ( mixing_length_1d == 'as_in_3d_model' ) THEN 879 DO k = nzb+1, nzt 880 dpt_dz = ( pt_init(k+1) - pt_init(k-1) ) * dd2zu(k) 881 IF ( dpt_dz > 0.0_wp ) THEN 882 l_stable = 0.76_wp * SQRT( e1d(k) ) & 883 / SQRT( g / pt_init(k) * dpt_dz ) + 1E-5_wp 884 ELSE 885 l_stable = l1d_init(k) 886 ENDIF 887 l1d(k) = MIN( l1d_init(k), l_stable ) 888 l1d_diss(k) = l1d(k) 889 ENDDO 890 ENDIF 891 ELSE 832 892 DO k = nzb+1, nzt 833 IF ( rif1d(k) >= 0.0_wp ) THEN 834 l1d(k) = l1d_init(k) / ( 1.0_wp + 5.0_wp * rif1d(k) ) 835 l1d_diss(k) = l1d(k) 836 ELSE 837 l1d(k) = l1d_init(k) 838 l1d_diss(k) = l1d_init(k) * & 839 SQRT( 1.0_wp - 16.0_wp * rif1d(k) ) 840 ENDIF 841 ENDDO 842 ELSEIF ( mixing_length_1d == 'as_in_3d_model' ) THEN 843 DO k = nzb+1, nzt 844 dpt_dz = ( pt_init(k+1) - pt_init(k-1) ) * dd2zu(k) 845 IF ( dpt_dz > 0.0_wp ) THEN 846 l_stable = 0.76_wp * SQRT( e1d(k) ) / & 847 SQRT( g / pt_init(k) * dpt_dz ) + 1E-5_wp 848 ELSE 849 l_stable = l1d_init(k) 850 ENDIF 851 l1d(k) = MIN( l1d_init(k), l_stable ) 852 l1d_diss(k) = l1d(k) 893 l1d(k) = c_0**3 * e1d(k) * SQRT( e1d(k) ) & 894 / ( diss1d(k) + 1.0E-30_wp ) 853 895 ENDDO 854 896 ENDIF … … 867 909 ENDIF 868 910 ENDIF 911 869 912 IF ( dissipation_1d == 'prognostic' ) THEN 870 913 DO k = nzb_diff, nzt … … 873 916 ELSE 874 917 DO k = nzb_diff, nzt 875 km1d(k) = c_ m* SQRT( e1d(k) ) * l1d(k)918 km1d(k) = c_0 * SQRT( e1d(k) ) * l1d(k) 876 919 ENDDO 877 920 ENDIF … … 1018 1061 1019 1062 REAL(wp) :: dt_diff !< time step accorind to diffusion criterion 1063 REAL(wp) :: dt_old !< previous time step 1020 1064 REAL(wp) :: fac !< factor of criterion 1021 1065 REAL(wp) :: value !< auxiliary variable 1022 1066 1023 1067 ! 1068 !-- Save previous time step 1069 dt_old = dt_1d 1070 1071 ! 1024 1072 !-- Compute the currently feasible time step according to the diffusion 1025 1073 !-- criterion. At nzb+1 the half grid length is used. 1026 fac = 0.125 !0.35_wp !### changed from 0.351074 fac = 0.125 1027 1075 dt_diff = dt_max_1d 1028 1076 DO k = nzb+2, nzt … … 1034 1082 1035 1083 ! 1084 !-- Limit the new time step to a maximum of 10 times the previous time step 1085 dt_1d = MIN( dt_old * 10.0_wp, dt_1d ) 1086 1087 ! 1036 1088 !-- Set flag when the time step becomes too small 1037 IF ( dt_1d < ( 0.00001_wp * dt_max_1d ) ) THEN1089 IF ( dt_1d < ( 1.0E-15_wp * dt_max_1d ) ) THEN 1038 1090 stop_dt_1d = .TRUE. 1039 1091
Note: See TracChangeset
for help on using the changeset viewer.
