Changeset 3083 for palm/trunk/SOURCE/turbulence_closure_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/turbulence_closure_mod.f90
r3045 r3083 25 25 ! ----------------- 26 26 ! $Id$ 27 ! - set limits of diss at the end of prognostic equations 28 ! - call production_e to calculate production term of diss 29 ! - limit change of diss to -90% to +100% 30 ! - remove factor 0.5 from diffusion_diss_ij 31 ! - rename c_m into c_0, and c_h into c_4 32 ! - add rans_const_c and rans_const_sigma as namelist parameters 33 ! - add calculation of mixing length for profile output in case of rans_tke_e 34 ! - changed format of annotations to comply with doxygen standards 35 ! - calculate and save dissipation rate during rans_tke_l mode 36 ! - set bc at vertical walls for e, diss, km, kh 37 ! - bugfix: set l_wall = 0.0 within buildings 38 ! - set l_wall at bottom and top boundary (rans-mode) 39 ! - bugfix in production term for dissipation rate 40 ! - bugfix in diffusion of dissipation rate 41 ! - disable check for 1D model if rans_tke_e is used 42 ! - bugfixes for initialization (rans-mode): 43 ! - correction of dissipation-rate formula 44 ! - calculate km based on l_wall 45 ! - initialize diss if 1D model is not used 46 ! 47 ! 3045 2018-05-28 07:55:41Z Giersch 27 48 ! Error message revised 28 49 ! … … 36 57 ! Further todo's 37 58 ! 38 ! 2936 2018-03-27 14:49:27Z suehring59 ! 2936 2018-03-27 14:49:27Z gronemeier 39 60 ! - defined l_grid only within this module 40 61 ! - Moved l_wall definition from modules.f90 … … 81 102 !> add OpenMP directives whereever possible 82 103 !> remove debug output variables (dummy1, dummy2, dummy3) 83 !> @todo Move initialization of wall-mixing length from init_grid84 104 !> @todo Check for random disturbances 85 105 !> @note <Enter notes on the module> … … 146 166 147 167 148 REAL(wp) :: c_1 = 1.44_wp !< model constant for RANS mode 149 REAL(wp) :: c_2 = 1.92_wp !< model constant for RANS mode 150 REAL(wp) :: c_3 = 1.44_wp !< model constant for RANS mode 151 REAL(wp) :: c_h = 0.0015_wp !< model constant for RANS mode 152 REAL(wp) :: c_m !< constant used for diffusion coefficient and dissipation (dependent on mode RANS/LES) 153 REAL(wp) :: c_mu = 0.09_wp !< model constant for RANS mode 154 REAL(wp) :: l_max !< maximum length scale for Blackadar mixing length 155 REAL(wp) :: sig_e = 1.0_wp !< factor to calculate Ke from Km 156 REAL(wp) :: sig_diss = 1.3_wp !< factor to calculate K_diss from Km 157 INTEGER(iwp) :: surf_e !< end index of surface elements at given i-j position 158 INTEGER(iwp) :: surf_s !< start index of surface elements at given i-j position 168 REAL(wp) :: c_0 !< constant used for diffusion coefficient and dissipation (dependent on mode RANS/LES) 169 REAL(wp) :: c_1 !< model constant for RANS mode 170 REAL(wp) :: c_2 !< model constant for RANS mode 171 REAL(wp) :: c_3 !< model constant for RANS mode 172 REAL(wp) :: c_4 !< model constant for RANS mode 173 REAL(wp) :: l_max !< maximum length scale for Blackadar mixing length 174 REAL(wp) :: dsig_e = 1.0_wp !< factor to calculate Ke from Km (1/sigma_e) 175 REAL(wp) :: dsig_diss = 1.0_wp !< factor to calculate K_diss from Km (1/sigma_diss) 176 INTEGER(iwp) :: surf_e !< end index of surface elements at given i-j position 177 INTEGER(iwp) :: surf_s !< start index of surface elements at given i-j position 178 179 REAL(wp), DIMENSION(0:4) :: rans_const_c = & !< model constants for RANS mode (namelist param) 180 (/ 0.55_wp, 1.44_wp, 1.92_wp, 0.0_wp, 0.0_wp /) !> default values fit for standard-tke-e closure 181 182 REAL(wp), DIMENSION(2) :: rans_const_sigma = & !< model constants for RANS mode, sigma values (sigma_e, sigma_diss) (namelist param) 183 (/ 1.0_wp, 0.77_wp /) 159 184 160 185 REAL(wp), DIMENSION(:), ALLOCATABLE :: l_black !< mixing length according to Blackadar … … 163 188 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: l_wall !< near-wall mixing length 164 189 165 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: dummy1 !< debug output variable 166 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: dummy2 !< debug output variable 167 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: dummy3 !< debug output variable 168 169 170 PUBLIC c_m, dummy1, dummy2, dummy3 190 !> @todo remove debug variables 191 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: diss_prod1, diss_adve1, diss_diff1, & 192 diss_prod2, diss_adve2, diss_diff2, & 193 diss_prod3, diss_adve3, diss_diff3, & 194 dummy1, dummy2, dummy3 195 196 197 PUBLIC c_0, rans_const_c, rans_const_sigma 171 198 172 199 ! … … 294 321 ! ------------ 295 322 !> Check parameters routine for turbulence closure module. 323 !> @todo remove rans_mode from initialization namelist and rework checks 324 !> The way it is implemented at the moment, the user has to set two variables 325 !> so that the RANS mode is working. It would be better if only one parameter 326 !> has to be set. 327 !> 2018-06-18, gronemeier 296 328 !------------------------------------------------------------------------------! 297 329 SUBROUTINE tcm_check_parameters … … 307 339 IF ( rans_mode ) THEN 308 340 309 c_m = 0.4_wp !according to Detering and Etling (1985) 341 ! 342 !-- Assign values to constants for RANS mode 343 dsig_e = 1.0_wp / rans_const_sigma(1) 344 dsig_diss = 1.0_wp / rans_const_sigma(2) 345 346 c_0 = rans_const_c(0) 347 c_1 = rans_const_c(1) 348 c_2 = rans_const_c(2) 349 !c_3 = rans_const_c(3) !> @todo clarify how to switch between different models 350 c_4 = rans_const_c(4) 310 351 311 352 SELECT CASE ( TRIM( turbulence_closure ) ) … … 316 357 CASE ( 'TKE-e' ) 317 358 rans_tke_e = .TRUE. 318 319 IF ( INDEX( initializing_actions, 'set_1d-model_profiles' ) == 0 &320 .AND. .NOT. nest_domain ) THEN321 message_string = 'Initializing without 1D model while ' // &322 'using TKE-e closure&' // &323 'is not possible at the moment!'324 CALL message( 'tcm_check_parameters', 'TG0005', 1, 2, 0, 6, 0 )325 ENDIF326 359 327 360 CASE DEFAULT 328 361 message_string = 'Unknown turbulence closure: ' // & 329 362 TRIM( turbulence_closure ) 330 CALL message( 'tcm_check_parameters', ' TG0001', 1, 2, 0, 6, 0 )363 CALL message( 'tcm_check_parameters', 'PA0500', 1, 2, 0, 6, 0 ) 331 364 332 365 END SELECT 366 367 IF ( turbulent_inflow .OR. turbulent_outflow ) THEN 368 message_string = 'turbulent inflow/outflow is not yet '// & 369 'implemented for RANS mode' 370 CALL message( 'tcm_check_parameters', 'PA0501', 1, 2, 0, 6, 0 ) 371 ENDIF 333 372 334 373 message_string = 'RANS mode is still in development! ' // & 335 374 '&Not all features of PALM are yet compatible '// & 336 375 'with RANS mode. &Use at own risk!' 337 CALL message( 'tcm_check_parameters', ' TG0003', 0, 1, 0, 6, 0 )376 CALL message( 'tcm_check_parameters', 'PA0502', 0, 1, 0, 6, 0 ) 338 377 339 378 ELSE 340 379 341 c_m = 0.1_wp !according to Lilly (1967) and Deardorff (1980) 380 c_0 = 0.1_wp !according to Lilly (1967) and Deardorff (1980) 381 382 dsig_e = 1.0_wp !assure to use K_m to calculate TKE instead 383 !of K_e which is used in RANS mode 342 384 343 385 SELECT CASE ( TRIM( turbulence_closure ) ) … … 347 389 348 390 CASE DEFAULT 391 !> @todo rework this part so that only one call of this error exists 349 392 message_string = 'Unknown turbulence closure: ' // & 350 393 TRIM( turbulence_closure ) 351 CALL message( 'tcm_check_parameters', ' TG0001', 1, 2, 0, 6, 0 )394 CALL message( 'tcm_check_parameters', 'PA0500', 1, 2, 0, 6, 0 ) 352 395 353 396 END SELECT 354 355 ENDIF356 357 IF ( rans_tke_e ) THEN358 359 IF ( turbulent_inflow .OR. turbulent_outflow ) THEN360 message_string = 'turbulent inflow/outflow is not yet '// &361 'implemented for TKE-e closure'362 CALL message( 'tcm_check_parameters', 'TG0002', 1, 2, 0, 6, 0 )363 ENDIF364 397 365 398 ENDIF … … 379 412 IMPLICIT NONE 380 413 381 CHARACTER (LEN=*) :: unit !< 382 CHARACTER (LEN=*) :: var !< 383 384 INTEGER(iwp) :: i !< 385 INTEGER(iwp) :: ilen !< 386 INTEGER(iwp) :: k !< 414 CHARACTER (LEN=*) :: unit !< unit of output variable 415 CHARACTER (LEN=*) :: var !< name of output variable 416 417 INTEGER(iwp) :: i !< index of var in data_output 418 INTEGER(iwp) :: ilen !< length of var string 419 INTEGER(iwp) :: k !< flag if var contains one of '_xy', '_xz' or '_yz' 387 420 388 421 SELECT CASE ( TRIM( var ) ) 389 422 390 423 CASE ( 'diss' ) 391 IF ( .NOT. rans_tke_e ) THEN392 message_string = 'output of "' // TRIM( var ) // '" requi' // &393 'res TKE-e closure for RANS mode.'394 CALL message( 'tcm_check_data_output', 'TG0101', 1, 2, 0, 6, 0 )395 ENDIF396 424 unit = 'm2/s3' 397 425 398 CASE ( 'dummy2', 'dummy3', 'dummy1' ) 399 unit = 'mixing length' 426 CASE ( 'diss1', 'diss2', & !> @todo remove later 427 'diss_prod1', 'diss_adve1', 'diss_diff1', & 428 'diss_prod2', 'diss_adve2', 'diss_diff2', & 429 'diss_prod3', 'diss_adve3', 'diss_diff3', 'dummy3' ) 430 unit = 'debug output' 400 431 401 432 CASE ( 'kh', 'km' ) … … 420 451 IMPLICIT NONE 421 452 422 CHARACTER (LEN=*), INTENT(OUT) :: grid_x !< 423 CHARACTER (LEN=*), INTENT(OUT) :: grid_y !< 424 CHARACTER (LEN=*), INTENT(OUT) :: grid_z !< 425 CHARACTER (LEN=*), INTENT(IN) :: var !< 426 427 LOGICAL, INTENT(OUT) :: found !< 428 453 CHARACTER (LEN=*), INTENT(OUT) :: grid_x !< x grid of output variable 454 CHARACTER (LEN=*), INTENT(OUT) :: grid_y !< y grid of output variable 455 CHARACTER (LEN=*), INTENT(OUT) :: grid_z !< z grid of output variable 456 CHARACTER (LEN=*), INTENT(IN) :: var !< name of output variable 457 458 LOGICAL, INTENT(OUT) :: found !< flag if output variable is found 459 429 460 found = .TRUE. 430 461 … … 438 469 grid_z = 'zu' 439 470 440 CASE ( 'dummy2', 'dummy3', 'dummy1' ) !### remove later 471 CASE ( 'diss1', 'diss2', & !> @todo remove later 472 'diss_prod1', 'diss_adve1', 'diss_diff1', & 473 'diss_prod2', 'diss_adve2', 'diss_diff2', & 474 'diss_prod3', 'diss_adve3', 'diss_diff3', 'dummy3' ) 441 475 grid_x = 'x' 442 476 grid_y = 'y' … … 480 514 IMPLICIT NONE 481 515 482 CHARACTER (LEN=*) :: mode !< 483 CHARACTER (LEN=*) :: variable !< 484 485 INTEGER(iwp) :: i !< 486 INTEGER(iwp) :: j !< 487 INTEGER(iwp) :: k !< 516 CHARACTER (LEN=*) :: mode !< flag defining mode 'allocate', 'sum' or 'average' 517 CHARACTER (LEN=*) :: variable !< name of variable 518 519 INTEGER(iwp) :: i !< loop index 520 INTEGER(iwp) :: j !< loop index 521 INTEGER(iwp) :: k !< loop index 488 522 489 523 IF ( mode == 'allocate' ) THEN … … 616 650 IMPLICIT NONE 617 651 618 CHARACTER (LEN=*) :: grid !< 619 CHARACTER (LEN=*) :: mode !< 620 CHARACTER (LEN=*) :: variable !< 621 622 INTEGER(iwp) :: av !< 623 INTEGER(iwp) :: i !< 624 INTEGER(iwp) :: j !< 625 INTEGER(iwp) :: k !< 626 INTEGER(iwp) :: nzb_do !< 627 INTEGER(iwp) :: nzt_do !< 628 629 LOGICAL :: found !< 652 CHARACTER (LEN=*) :: grid !< name of vertical grid 653 CHARACTER (LEN=*) :: mode !< either 'xy', 'xz' or 'yz' 654 CHARACTER (LEN=*) :: variable !< name of variable 655 656 INTEGER(iwp) :: av !< flag for (non-)average output 657 INTEGER(iwp) :: i !< loop index 658 INTEGER(iwp) :: j !< loop index 659 INTEGER(iwp) :: k !< loop index 660 INTEGER(iwp) :: nzb_do !< vertical output index (bottom) 661 INTEGER(iwp) :: nzt_do !< vertical output index (top) 662 663 LOGICAL :: found !< flag if output variable is found 630 664 LOGICAL :: two_d !< flag parameter that indicates 2D variables (horizontal cross sections) 631 665 … … 737 771 IMPLICIT NONE 738 772 739 CHARACTER (LEN=*) :: variable !< 740 741 INTEGER(iwp) :: av !< 742 INTEGER(iwp) :: i !< 743 INTEGER(iwp) :: j !< 744 INTEGER(iwp) :: k !< 773 CHARACTER (LEN=*) :: variable !< name of variable 774 775 INTEGER(iwp) :: av !< flag for (non-)average output 776 INTEGER(iwp) :: i !< loop index 777 INTEGER(iwp) :: j !< loop index 778 INTEGER(iwp) :: k !< loop index 745 779 INTEGER(iwp) :: nzb_do !< lower limit of the data output (usually 0) 746 780 INTEGER(iwp) :: nzt_do !< vertical upper limit of the data output (usually nz_do3d) 747 781 748 LOGICAL :: found !< 782 LOGICAL :: found !< flag if output variable is found 749 783 750 784 REAL(wp) :: fill_value = -999.0_wp !< value for the _FillValue attribute … … 829 863 ENDIF 830 864 831 CASE ( 'dummy1' ) !### remove later 865 CASE ( 'dummy3' ) !> @todo remove later 866 IF ( av == 0 ) THEN 867 DO i = nxl, nxr 868 DO j = nys, nyn 869 DO k = nzb_do, nzt_do 870 local_pf(i,j,k) = dummy3(k,j,i) 871 ENDDO 872 ENDDO 873 ENDDO 874 ENDIF 875 876 CASE ( 'diss1' ) !> @todo remove later 832 877 IF ( av == 0 ) THEN 833 878 DO i = nxl, nxr … … 840 885 ENDIF 841 886 842 CASE ( 'd ummy2' ) !###remove later887 CASE ( 'diss2' ) !> @todo remove later 843 888 IF ( av == 0 ) THEN 844 889 DO i = nxl, nxr … … 851 896 ENDIF 852 897 853 CASE ( 'd ummy3' ) !###remove later898 CASE ( 'diss_prod1' ) !> @todo remove later 854 899 IF ( av == 0 ) THEN 855 900 DO i = nxl, nxr 856 901 DO j = nys, nyn 857 902 DO k = nzb_do, nzt_do 858 local_pf(i,j,k) = d ummy3(k,j,i)903 local_pf(i,j,k) = diss_prod1(k,j,i) 859 904 ENDDO 860 905 ENDDO … … 862 907 ENDIF 863 908 909 CASE ( 'diss_adve1' ) !> @todo remove later 910 IF ( av == 0 ) THEN 911 DO i = nxl, nxr 912 DO j = nys, nyn 913 DO k = nzb_do, nzt_do 914 local_pf(i,j,k) = diss_adve1(k,j,i) 915 ENDDO 916 ENDDO 917 ENDDO 918 ENDIF 919 920 CASE ( 'diss_diff1' ) !> @todo remove later 921 IF ( av == 0 ) THEN 922 DO i = nxl, nxr 923 DO j = nys, nyn 924 DO k = nzb_do, nzt_do 925 local_pf(i,j,k) = diss_diff1(k,j,i) 926 ENDDO 927 ENDDO 928 ENDDO 929 ENDIF 930 931 CASE ( 'diss_prod2' ) !> @todo remove later 932 IF ( av == 0 ) THEN 933 DO i = nxl, nxr 934 DO j = nys, nyn 935 DO k = nzb_do, nzt_do 936 local_pf(i,j,k) = diss_prod2(k,j,i) 937 ENDDO 938 ENDDO 939 ENDDO 940 ENDIF 941 942 CASE ( 'diss_adve2' ) !> @todo remove later 943 IF ( av == 0 ) THEN 944 DO i = nxl, nxr 945 DO j = nys, nyn 946 DO k = nzb_do, nzt_do 947 local_pf(i,j,k) = diss_adve2(k,j,i) 948 ENDDO 949 ENDDO 950 ENDDO 951 ENDIF 952 953 CASE ( 'diss_diff2' ) !> @todo remove later 954 IF ( av == 0 ) THEN 955 DO i = nxl, nxr 956 DO j = nys, nyn 957 DO k = nzb_do, nzt_do 958 local_pf(i,j,k) = diss_diff2(k,j,i) 959 ENDDO 960 ENDDO 961 ENDDO 962 ENDIF 963 964 CASE ( 'diss_prod3' ) !> @todo remove later 965 IF ( av == 0 ) THEN 966 DO i = nxl, nxr 967 DO j = nys, nyn 968 DO k = nzb_do, nzt_do 969 local_pf(i,j,k) = diss_prod3(k,j,i) 970 ENDDO 971 ENDDO 972 ENDDO 973 ENDIF 974 975 CASE ( 'diss_adve3' ) !> @todo remove later 976 IF ( av == 0 ) THEN 977 DO i = nxl, nxr 978 DO j = nys, nyn 979 DO k = nzb_do, nzt_do 980 local_pf(i,j,k) = diss_adve3(k,j,i) 981 ENDDO 982 ENDDO 983 ENDDO 984 ENDIF 985 986 CASE ( 'diss_diff3' ) !> @todo remove later 987 IF ( av == 0 ) THEN 988 DO i = nxl, nxr 989 DO j = nys, nyn 990 DO k = nzb_do, nzt_do 991 local_pf(i,j,k) = diss_diff3(k,j,i) 992 ENDDO 993 ENDDO 994 ENDDO 995 ENDIF 996 864 997 CASE DEFAULT 865 998 found = .FALSE. … … 891 1024 ALLOCATE( km(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 892 1025 893 ALLOCATE( dummy1(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) ! ###remove later1026 ALLOCATE( dummy1(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) !> @todo remove later 894 1027 ALLOCATE( dummy2(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 895 1028 ALLOCATE( dummy3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 896 897 IF ( rans_mode ) ALLOCATE( l_black(nzb:nzt+1) ) 1029 ALLOCATE( diss_adve1(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 1030 ALLOCATE( diss_adve2(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 1031 ALLOCATE( diss_adve3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 1032 ALLOCATE( diss_prod1(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 1033 ALLOCATE( diss_prod2(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 1034 ALLOCATE( diss_prod3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 1035 ALLOCATE( diss_diff1(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 1036 ALLOCATE( diss_diff2(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 1037 ALLOCATE( diss_diff3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 1038 dummy1 = 0.0_wp 1039 dummy2 = 0.0_wp 1040 dummy3 = 0.0_wp 1041 diss_adve1 = 0.0_wp 1042 diss_adve2 = 0.0_wp 1043 diss_adve3 = 0.0_wp 1044 diss_prod1 = 0.0_wp 1045 diss_prod2 = 0.0_wp 1046 diss_prod3 = 0.0_wp 1047 diss_diff1 = 0.0_wp 1048 diss_diff2 = 0.0_wp 1049 diss_diff3 = 0.0_wp 898 1050 899 1051 #if defined( __nopointer ) … … 912 1064 !-- they do not necessarily need to be transferred, which is attributed to 913 1065 !-- the design of the model coupler which allocates memory for each variable. 914 IF ( rans_ tke_e .OR. use_sgs_for_particles .OR. wang_kernel .OR.&1066 IF ( rans_mode .OR. use_sgs_for_particles .OR. wang_kernel .OR. & 915 1067 collision_turbulence .OR. nested_run ) THEN 916 1068 #if defined( __nopointer ) … … 934 1086 e => e_1; e_p => e_2; te_m => e_3 935 1087 936 IF ( rans_ tke_e .OR. use_sgs_for_particles .OR. &1088 IF ( rans_mode .OR. use_sgs_for_particles .OR. & 937 1089 wang_kernel .OR. collision_turbulence .OR. nested_run ) THEN 938 1090 diss => diss_1 … … 967 1119 INTEGER(iwp) :: j !< loop index 968 1120 INTEGER(iwp) :: k !< loop index 969 INTEGER(iwp) :: nz_s_shift !< 970 INTEGER(iwp) :: nz_s_shift_l !< 1121 INTEGER(iwp) :: nz_s_shift !< lower shift index for scalars 1122 INTEGER(iwp) :: nz_s_shift_l !< local lower shift index in case of turbulent inflow 971 1123 972 1124 ! 973 1125 !-- Initialize mixing length 974 1126 CALL tcm_init_mixing_length 1127 dummy3 = l_wall !> @todo remove later 975 1128 976 1129 ! … … 995 1148 996 1149 IF ( rans_tke_e ) THEN 997 IF ( dissipation_1d == 'prognostic' ) THEN ! ###Why must this be checked?998 DO i = nxlg, nxrg ! ###Should 'diss' not always999 DO j = nysg, nyng ! ###be prognostic in case rans_tke_e?1150 IF ( dissipation_1d == 'prognostic' ) THEN !> @query Why must this be checked? 1151 DO i = nxlg, nxrg !> Should 'diss' not always 1152 DO j = nysg, nyng !> be prognostic in case rans_tke_e? 1000 1153 diss(:,j,i) = diss1d 1001 1154 ENDDO … … 1005 1158 DO j = nysg, nyng 1006 1159 DO k = nzb+1, nzt 1007 diss(k,j,i) = e(k,j,i) * SQRT( e(k,j,i) ) / l1d(k)1160 diss(k,j,i) = c_0**4 * e(k,j,i)**2 / km1d(k) 1008 1161 ENDDO 1009 1162 ENDDO … … 1016 1169 1017 1170 IF ( constant_diffusion ) THEN 1018 km 1019 kh 1020 e 1171 km = km_constant 1172 kh = km / prandtl_number 1173 e = 0.0_wp 1021 1174 ELSEIF ( e_init > 0.0_wp ) THEN 1022 DO k = nzb+1, nzt 1023 km(k,:,:) = c_m * l_grid(k) * SQRT( e_init ) 1175 DO i = nxlg, nxrg 1176 DO j = nysg, nyng 1177 DO k = nzb+1, nzt 1178 km(k,j,i) = c_0 * l_wall(k,j,i) * SQRT( e_init ) 1179 ENDDO 1180 ENDDO 1024 1181 ENDDO 1025 1182 km(nzb,:,:) = km(nzb+1,:,:) 1026 1183 km(nzt+1,:,:) = km(nzt,:,:) 1027 kh 1028 e 1184 kh = km / prandtl_number 1185 e = e_init 1029 1186 ELSE 1030 1187 IF ( .NOT. ocean ) THEN … … 1040 1197 ENDIF 1041 1198 1199 IF ( rans_tke_e ) THEN 1200 DO i = nxlg, nxrg 1201 DO j = nysg, nyng 1202 DO k = nzb+1, nzt 1203 diss(k,j,i) = c_0**4 * e(k,j,i)**2 / km(k,j,i) 1204 ENDDO 1205 ENDDO 1206 ENDDO 1207 diss(nzb,:,:) = diss(nzb+1,:,:) 1208 diss(nzt+1,:,:) = diss(nzt,:,:) 1209 ENDIF 1210 1042 1211 ENDIF 1043 1212 ! … … 1073 1242 ENDDO 1074 1243 ENDDO 1244 IF ( rans_tke_e ) THEN 1245 DO i = nxlg, nxrg 1246 DO j = nysg, nyng 1247 nz_s_shift = get_topography_top_index_ji( j, i, 's' ) 1248 1249 diss(nz_s_shift:nzt+1,j,i) = diss(0:nzt+1-nz_s_shift,j,i) 1250 ENDDO 1251 ENDDO 1252 ENDIF 1075 1253 ENDIF 1076 1254 … … 1084 1262 !-- boundary and adjust mean inflow profiles 1085 1263 IF ( complex_terrain ) THEN 1086 IF ( nxlg <= 0 .AND. nxrg >= 0 .AND. nysg <= 0 .AND. nyng >= 0 ) THEN 1264 IF ( nxlg <= 0 .AND. nxrg >= 0 .AND. & 1265 nysg <= 0 .AND. nyng >= 0 ) THEN 1087 1266 nz_s_shift_l = get_topography_top_index_ji( 0, 0, 's' ) 1088 1267 ELSE … … 1095 1274 nz_s_shift = nz_s_shift_l 1096 1275 #endif 1097 mean_inflow_profiles(nz_s_shift:nzt+1,5) = hom_sum(0:nzt+1-nz_s_shift,8,0) ! e 1276 mean_inflow_profiles(nz_s_shift:nzt+1,5) = & 1277 hom_sum(0:nzt+1-nz_s_shift,8,0) ! e 1098 1278 ENDIF 1099 1279 ! … … 1114 1294 ! 1115 1295 !-- Inside buildings set TKE back to zero. 1116 !-- Other scalars (km, kh, diss,...) are ignored at present,1296 !-- Other scalars (km, kh,...) are ignored at present, 1117 1297 !-- maybe revise later. 1118 1298 DO i = nxlg, nxrg … … 1121 1301 e(k,j,i) = MERGE( e(k,j,i), 0.0_wp, & 1122 1302 BTEST( wall_flags_0(k,j,i), 0 ) ) 1123 te_m(k,j,i) = MERGE( te_m(k,j,i), 0.0_wp, &1124 BTEST( wall_flags_0(k,j,i), 0 ) )1125 1303 ENDDO 1126 1304 ENDDO 1127 1305 ENDDO 1128 1306 1307 IF ( rans_tke_e ) THEN 1308 DO i = nxlg, nxrg 1309 DO j = nysg, nyng 1310 DO k = nzb, nzt 1311 diss(k,j,i) = MERGE( diss(k,j,i), 0.0_wp, & 1312 BTEST( wall_flags_0(k,j,i), 0 ) ) 1313 ENDDO 1314 ENDDO 1315 ENDDO 1316 ENDIF 1129 1317 ENDIF 1130 1318 ! … … 1134 1322 ! 1135 1323 !-- Allthough tendency arrays are set in prognostic_equations, they have 1136 !-- have to be predefined here because they are used (but multiplied with 0)1137 !-- there before they are set.1324 !-- to be predefined here because there they are used (but multiplied with 0) 1325 !-- before they are set. 1138 1326 te_m = 0.0_wp 1327 1328 IF ( rans_tke_e ) THEN 1329 diss_p = diss 1330 tdiss_m = 0.0_wp 1331 ENDIF 1139 1332 1140 1333 ENDIF … … 1215 1408 1216 1409 DO k = 1, nzt 1217 IF ( l_grid(k) > 1.5_wp * dx * wall_adjustment_factor .OR. &1410 IF ( l_grid(k) > 1.5_wp * dx * wall_adjustment_factor .OR. & 1218 1411 l_grid(k) > 1.5_wp * dy * wall_adjustment_factor ) THEN 1219 WRITE( message_string, * ) 'grid anisotropy exceeds ', &1220 'threshold given by only local', &1221 ' horizontal reduction of near_wall ', &1222 'mixing length l_wall', &1223 ' starting from height level k = ', k, &1412 WRITE( message_string, * ) 'grid anisotropy exceeds ', & 1413 'threshold given by only local', & 1414 ' &horizontal reduction of near_wall ', & 1415 'mixing length l_wall', & 1416 ' &starting from height level k = ', k, & 1224 1417 '.' 1225 1418 CALL message( 'init_grid', 'PA0202', 0, 1, 0, 6, 0 ) … … 1320 1513 ! 1321 1514 !-- Initialize the mixing length in case of a RANS simulation 1515 ALLOCATE( l_black(nzb:nzt+1) ) 1322 1516 1323 1517 ! 1324 1518 !-- Calculate mixing length according to Blackadar (1962) 1325 1519 IF ( f /= 0.0_wp ) THEN 1326 l_max = 2.7E-4 * SQRT( ug(nzt+1)**2 + vg(nzt+1)**2 ) /&1327 ABS( f ) + 1 E-10_wp1520 l_max = 2.7E-4_wp * SQRT( ug(nzt+1)**2 + vg(nzt+1)**2 ) / & 1521 ABS( f ) + 1.0E-10_wp 1328 1522 ELSE 1329 1523 l_max = 30.0_wp … … 1338 1532 ! 1339 1533 !-- Gather topography information of whole domain 1340 ! ## TODO:reduce amount of data sent by MPI call1341 ! ##By now, a whole global 3D-array is sent and received with1342 ! ##MPI_ALLREDUCE although most of the array is 0. This can be1343 ! ##drastically reduced if only the local subarray is sent and stored1344 ! ##in a global array. For that, an MPI data type or subarray must be1345 ! ##defined.1346 ! ##2018-03-19, gronemeier1534 !> @todo reduce amount of data sent by MPI call 1535 !> By now, a whole global 3D-array is sent and received with 1536 !> MPI_ALLREDUCE although most of the array is 0. This can be 1537 !> drastically reduced if only the local subarray is sent and stored 1538 !> in a global array. For that, an MPI data type or subarray must be 1539 !> defined. 1540 !> 2018-03-19, gronemeier 1347 1541 ALLOCATE( wall_flags_0_global(nzb:nzt+1,0:ny,0:nx) ) 1348 1542 … … 1373 1567 ENDDO 1374 1568 ENDDO 1569 1570 l_wall(nzb,:,:) = l_black(nzb) 1571 l_wall(nzt+1,:,:) = l_black(nzt+1) 1375 1572 ! 1376 1573 !-- Limit mixing length to either nearest wall or Blackadar mixing length. … … 1415 1612 IF ( rad_k_b /= 0 .OR. rad_k_t /= 0 ) THEN 1416 1613 1417 ! ## NOTE:shape of vicinity is larger in z direction1418 ! ##Shape of vicinity is two grid points larger than actual search1419 ! ##radius in vertical direction. The first and last grid point is1420 ! ##always set to 1 to asure correct detection of topography. See1421 ! ##function "shortest_distance" for details.1422 ! ##2018-03-16, gronemeier1614 !> @note shape of vicinity is larger in z direction 1615 !> Shape of vicinity is two grid points larger than actual search 1616 !> radius in vertical direction. The first and last grid point is 1617 !> always set to 1 to asure correct detection of topography. See 1618 !> function "shortest_distance" for details. 1619 !> 2018-03-16, gronemeier 1423 1620 ALLOCATE( vicinity(-rad_k-1:rad_k+1,-rad_j:rad_j,-rad_i:rad_i) ) 1424 1621 ALLOCATE( vic_yz(0:rad_k+1,0:rad_j) ) … … 1601 1798 ELSE !Check if (i,j,k) belongs to atmosphere 1602 1799 1603 l_wall(k,j,i) = -999.01800 l_wall(k,j,i) = l_black(k) 1604 1801 1605 1802 ENDIF … … 1630 1827 !> (pos_i/jj/kk), where (jj/kk) is the position of the maximum of 'array' 1631 1828 !> closest to the origin (0/0) of 'array'. 1829 !> @todo this part of PALM does not reproduce the same results for optimized 1830 !> and debug options for the compiler. This should be fixed 1632 1831 !------------------------------------------------------------------------------! 1633 1832 REAL FUNCTION shortest_distance( array, orientation, pos_i ) … … 1659 1858 IF ( orientation ) THEN !if array is oriented upwards 1660 1859 DO jj = 0, rad_j 1661 shortest_distance = MIN( shortest_distance, & 1662 SQRT( MAX(pos_i*dx-0.5*dx,0.0)**2 & 1663 + MAX(jj*dy-0.5*dy,0.0)**2 & 1664 + MAX(zw(loc_k(jj)+k-1)-zu(k),0.0)**2 & 1665 ) & 1666 ) 1860 shortest_distance = & 1861 MIN( shortest_distance, & 1862 SQRT( MAX(REAL(pos_i, KIND=wp)*dx-0.5_wp*dx, 0.0_wp)**2 & 1863 + MAX(REAL(jj, KIND=wp)*dy-0.5_wp*dy, 0.0_wp)**2 & 1864 + MAX(zw(loc_k(jj)+k-1)-zu(k), 0.0_wp)**2 & 1865 ) & 1866 ) 1667 1867 ENDDO 1668 1868 ELSE !if array is oriented downwards 1669 ! ## NOTE:MAX within zw required to circumvent error at domain border1670 ! ##At the domain border, if non-cyclic boundary is present, the1671 ! ##index for zw could be -1, which will be errorneous (zw(-1) does1672 ! ##not exist). The MAX function limits the index to be at least 0.1869 !> @note MAX within zw required to circumvent error at domain border 1870 !> At the domain border, if non-cyclic boundary is present, the 1871 !> index for zw could be -1, which will be errorneous (zw(-1) does 1872 !> not exist). The MAX function limits the index to be at least 0. 1673 1873 DO jj = 0, rad_j 1674 shortest_distance = MIN( shortest_distance, & 1675 SQRT( MAX(pos_i*dx-0.5*dx,0.0)**2 & 1676 + MAX(jj*dy-0.5*dy,0.0)**2 & 1677 + MAX(zu(k)-zw(MAX(k-loc_k(jj), & 1678 0_iwp)), & 1679 0.0)**2 & 1680 ) & 1681 ) 1874 shortest_distance = & 1875 MIN( shortest_distance, & 1876 SQRT( MAX(REAL(pos_i, KIND=wp)*dx-0.5_wp*dx, 0.0_wp)**2 & 1877 + MAX(REAL(jj, KIND=wp)*dy-0.5_wp*dy, 0.0_wp)**2 & 1878 + MAX(zu(k)-zw(MAX(k-loc_k(jj),0_iwp)), 0.0_wp)**2 & 1879 ) & 1880 ) 1682 1881 ENDDO 1683 1882 ENDIF 1684 1883 1685 1884 END FUNCTION 1686 1885 … … 2004 2203 ! 2005 2204 !-- Use special boundary condition in case of TKE-e closure 2205 !> @todo do the same for usm and lsm surfaces 2206 !> 2018-06-05, gronemeier 2006 2207 IF ( rans_tke_e ) THEN 2007 2208 DO i = nxl, nxr … … 2011 2212 DO m = surf_s, surf_e 2012 2213 k = surf_def_h(0)%k(m) 2013 e_p(k,j,i) = surf_def_h(0)%us(m)**2 / c_ m**22214 e_p(k,j,i) = surf_def_h(0)%us(m)**2 / c_0**2 2014 2215 ENDDO 2015 2216 ENDDO … … 2092 2293 DO k = nzb+1, nzt 2093 2294 ! tend(k,j,i) = tend(k,j,i) + c_1 * diss(k,j,i) / ( e(k,j,i) + 1.0E-20_wp ) * produc(k) 2094 tend(k,j,i) = tend(k,j,i) + c_1 * c_ mu * f / c_h & !###needs revision2295 tend(k,j,i) = tend(k,j,i) + c_1 * c_0**4 * f / c_4 & !> @todo needs revision 2095 2296 / surf_def_h(0)%us(surf_def_h(0)%start_index(j,i)) & 2096 2297 * SQRT(e(k,j,i)) * produc(k,j,i) … … 2103 2304 ! 2104 2305 !-- Additional sink term for flows through plant canopies 2105 ! IF ( plant_canopy ) CALL pcm_tendency( ? ) ! ###what to do with this?2106 2107 ! CALL user_actions( 'diss-tendency' ) ! ###not yet implemented2306 ! IF ( plant_canopy ) CALL pcm_tendency( ? ) !> @query what to do with this? 2307 2308 ! CALL user_actions( 'diss-tendency' ) !> @todo not yet implemented 2108 2309 2109 2310 ! … … 2181 2382 USE arrays_3d, & 2182 2383 ONLY: ddzu, diss_l_diss, diss_l_e, diss_s_diss, diss_s_e, & 2183 flux_l_diss, flux_l_e, flux_s_diss, flux_s_e 2384 flux_l_diss, flux_l_e, flux_s_diss, flux_s_e,& 2385 u_p,v_p,w_p 2184 2386 2185 2387 USE control_parameters, & 2186 2388 ONLY: f, tsc 2389 2390 USE grid_variables, & 2391 ONLY: dx, dy 2187 2392 2188 2393 USE surface_mod, & … … 2190 2395 surf_usm_v 2191 2396 2397 use indices, only: nx, ny 2398 2192 2399 IMPLICIT NONE 2193 2400 2194 2401 INTEGER(iwp) :: i !< loop index x direction 2195 INTEGER(iwp) :: i_omp !< 2402 INTEGER(iwp) :: i_omp !< first loop index of i-loop in prognostic_equations 2196 2403 INTEGER(iwp) :: j !< loop index y direction 2197 2404 INTEGER(iwp) :: k !< loop index z direction 2405 INTEGER(iwp) :: l !< loop index 2198 2406 INTEGER(iwp) :: m !< loop index 2199 2407 INTEGER(iwp) :: surf_e !< end index of surface elements at given i-j position 2200 2408 INTEGER(iwp) :: surf_s !< start index of surface elements at given i-j position 2201 INTEGER(iwp) :: tn !< 2202 2203 REAL(wp), DIMENSION(nzb:nzt+1) :: advec !< advection term of TKE tendency 2204 REAL(wp), DIMENSION(nzb:nzt+1) :: produc !< production term of TKE tendency 2409 INTEGER(iwp) :: tn !< task number of openmp task 2410 2411 INTEGER(iwp) :: pis = 32 !< debug variable, print from i=pis !> @todo remove later 2412 INTEGER(iwp) :: pie = 32 !< debug variable, print until i=pie !> @todo remove later 2413 INTEGER(iwp) :: pjs = 26 !< debug variable, print from j=pjs !> @todo remove later 2414 INTEGER(iwp) :: pje = 26 !< debug variable, print until j=pje !> @todo remove later 2415 INTEGER(iwp) :: pkb = 1 !< debug variable, print from k=pkb !> @todo remove later 2416 INTEGER(iwp) :: pkt = 7 !< debug variable, print until k=pkt !> @todo remove later 2417 2418 REAL(wp), DIMENSION(nzb:nzt+1) :: dum_adv !< debug variable !> @todo remove later 2419 REAL(wp), DIMENSION(nzb:nzt+1) :: dum_pro !< debug variable !> @todo remove later 2420 REAL(wp), DIMENSION(nzb:nzt+1) :: dum_dif !< debug variable !> @todo remove later 2421 2422 5555 FORMAT(A,7(1X,E12.5)) !> @todo remove later 2205 2423 2206 2424 ! … … 2224 2442 ENDIF 2225 2443 2226 advec(:) = tend(:,j,i)2227 2228 CALL production_e( i, j )2229 2230 produc(:) = tend(:,j,i) - advec(:)2444 dum_adv = tend(:,j,i) !> @todo remove later 2445 2446 CALL production_e( i, j, .FALSE. ) 2447 2448 dum_pro = tend(:,j,i) - dum_adv !> @todo remove later 2231 2449 2232 2450 IF ( .NOT. humidity ) THEN … … 2239 2457 CALL diffusion_e( i, j, vpt, pt_reference ) 2240 2458 ENDIF 2459 2460 dum_dif = tend(:,j,i) - dum_adv - dum_pro !> @todo remove later 2241 2461 2242 2462 ! … … 2268 2488 DO m = surf_s, surf_e 2269 2489 k = surf_def_h(0)%k(m) 2270 e_p(k,j,i) = surf_def_h(0)%us(m)**2 / c_m**2 2490 e_p(k,j,i) = ( surf_def_h(0)%us(m) / c_0 )**2 2491 ENDDO 2492 2493 DO l = 0, 3 2494 surf_s = surf_def_v(l)%start_index(j,i) 2495 surf_e = surf_def_v(l)%end_index(j,i) 2496 DO m = surf_s, surf_e 2497 k = surf_def_v(l)%k(m) 2498 e_p(k,j,i) = ( surf_def_v(l)%us(m) / c_0 )**2 2499 ENDDO 2271 2500 ENDDO 2272 2501 ENDIF … … 2287 2516 ENDIF 2288 2517 ENDIF 2518 2519 ! if ( i >= pis .and. i <= pie .and. j >= pjs .and. j <= pje ) then !> @todo remove later 2520 ! WRITE(9, *) '------' 2521 ! WRITE(9, '(A,F8.3,1X,F8.3,1X,I2)') 't, dt, int_ts:', simulated_time, dt_3d, intermediate_timestep_count 2522 ! WRITE(9, *) 'i:', i 2523 ! WRITE(9, *) 'j:', j 2524 ! WRITE(9, *) 'k:', pkb, ' - ', pkt 2525 ! WRITE(9, *) '---' 2526 ! WRITE(9, *) 'e:' 2527 ! WRITE(9, 5555) 'adv :', dum_adv(pkb:pkt) 2528 ! WRITE(9, 5555) 'pro :', dum_pro(pkb:pkt) 2529 ! WRITE(9, 5555) 'dif :', dum_dif(pkb:pkt) 2530 ! WRITE(9, 5555) 'tend:', tend(pkb:pkt,j,i) 2531 ! WRITE(9, 5555) 'e_p :', e_p(pkb:pkt,j,i) 2532 ! WRITE(9, 5555) 'e :', e(pkb:pkt,j,i) 2533 ! FLUSH(9) 2534 ! endif 2289 2535 2290 2536 ENDIF ! TKE equation … … 2309 2555 ENDIF 2310 2556 2557 IF ( intermediate_timestep_count == 1 ) diss_adve1(:,j,i) = tend(:,j,i) !> @todo remove later 2558 IF ( intermediate_timestep_count == 2 ) diss_adve2(:,j,i) = tend(:,j,i) 2559 IF ( intermediate_timestep_count == 3 ) diss_adve3(:,j,i) = tend(:,j,i) 2560 2311 2561 ! 2312 2562 !-- Production of TKE dissipation rate 2313 DO k = nzb+1, nzt 2314 ! tend(k,j,i) = tend(k,j,i) + c_1 * diss(k,j,i) / ( e(k,j,i) + 1.0E-20_wp ) * produc(k) 2315 tend(k,j,i) = tend(k,j,i) + c_1 * c_mu * f / c_h & !### needs revision 2316 / surf_def_h(0)%us(surf_def_h(0)%start_index(j,i)) & 2317 * SQRT(e(k,j,i)) * produc(k) 2318 ENDDO 2319 2563 CALL production_e( i, j, .TRUE. ) 2564 2565 IF ( intermediate_timestep_count == 1 ) diss_prod1(:,j,i) = tend(:,j,i) - diss_adve1(:,j,i) !> @todo remove later 2566 IF ( intermediate_timestep_count == 2 ) diss_prod2(:,j,i) = tend(:,j,i) - diss_adve2(:,j,i) 2567 IF ( intermediate_timestep_count == 3 ) diss_prod3(:,j,i) = tend(:,j,i) - diss_adve3(:,j,i) 2568 2569 dum_pro = tend(:,j,i) - dum_adv !> @todo remove later 2570 2571 ! 2572 !-- Diffusion term of TKE dissipation rate 2320 2573 CALL diffusion_diss( i, j ) 2321 2574 2575 IF ( intermediate_timestep_count == 1 ) diss_diff1(:,j,i) = tend(:,j,i) - diss_adve1(:,j,i) - diss_prod1(:,j,i) !> @todo remove later 2576 IF ( intermediate_timestep_count == 2 ) diss_diff2(:,j,i) = tend(:,j,i) - diss_adve2(:,j,i) - diss_prod2(:,j,i) 2577 IF ( intermediate_timestep_count == 3 ) diss_diff3(:,j,i) = tend(:,j,i) - diss_adve3(:,j,i) - diss_prod3(:,j,i) 2578 IF ( intermediate_timestep_count == 3 ) dummy3(:,j,i) = km(:,j,i) 2579 2580 dum_dif = tend(:,j,i) - dum_adv - dum_pro !> @todo remove later 2581 2322 2582 ! 2323 2583 !-- Additional sink term for flows through plant canopies 2324 ! IF ( plant_canopy ) CALL pcm_tendency( i, j, ? ) ! ###not yet implemented2325 2326 ! CALL user_actions( i, j, 'diss-tendency' ) ! ###not yet implemented2584 ! IF ( plant_canopy ) CALL pcm_tendency( i, j, ? ) !> @todo not yet implemented 2585 2586 ! CALL user_actions( i, j, 'diss-tendency' ) !> @todo not yet implemented 2327 2587 2328 2588 ! … … 2338 2598 BTEST( wall_flags_0(k,j,i), 0 )& 2339 2599 ) 2340 IF ( diss_p(k,j,i) <= 0.0_wp ) diss_p(k,j,i) = 0.1_wp * diss(k,j,i)2341 2600 ENDDO 2342 2601 … … 2350 2609 ENDDO 2351 2610 2611 DO l = 0, 1 2612 surf_s = surf_def_v(l)%start_index(j,i) 2613 surf_e = surf_def_v(l)%end_index(j,i) 2614 DO m = surf_s, surf_e 2615 k = surf_def_v(l)%k(m) 2616 diss_p(k,j,i) = surf_def_v(l)%us(m)**3 / ( kappa * 0.5_wp * dy ) 2617 ENDDO 2618 ENDDO 2619 2620 DO l = 2, 3 2621 surf_s = surf_def_v(l)%start_index(j,i) 2622 surf_e = surf_def_v(l)%end_index(j,i) 2623 DO m = surf_s, surf_e 2624 k = surf_def_v(l)%k(m) 2625 diss_p(k,j,i) = surf_def_v(l)%us(m)**3 / ( kappa * 0.5_wp * dx ) 2626 ENDDO 2627 ENDDO 2352 2628 ! 2353 2629 !-- Calculate tendencies for the next Runge-Kutta step … … 2365 2641 ENDIF 2366 2642 ENDIF 2367 2368 ! IF ( intermediate_timestep_count == 1 ) dummy1(:,j,i) = e_p(:,j,i) 2369 ! IF ( intermediate_timestep_count == 2 ) dummy2(:,j,i) = e_p(:,j,i) 2370 ! IF ( intermediate_timestep_count == 3 ) dummy3(:,j,i) = e_p(:,j,i) 2643 ! 2644 !-- Limit change of diss to be between -90% and +100%. Also, set an absolute 2645 !-- minimum value 2646 DO k = nzb, nzt+1 2647 diss_p(k,j,i) = MIN( MAX( diss_p(k,j,i), & 2648 0.1_wp * diss(k,j,i), & 2649 0.0001_wp ), & 2650 2.0_wp * diss(k,j,i) ) 2651 ENDDO 2652 2653 IF ( intermediate_timestep_count == 1 ) dummy1(:,j,i) = diss_p(:,j,i) !> @todo remove later 2654 IF ( intermediate_timestep_count == 2 ) dummy2(:,j,i) = diss_p(:,j,i) 2655 2656 ! if ( i >= pis .and. i <= pie .and. j >= pjs .and. j <= pje ) then !> @todo remove later 2657 ! WRITE(9, *) '---' 2658 ! WRITE(9, *) 'diss:' 2659 ! WRITE(9, 5555) 'adv :', dum_adv(pkb:pkt) 2660 ! WRITE(9, 5555) 'pro :', dum_pro(pkb:pkt) 2661 ! WRITE(9, 5555) 'dif :', dum_dif(pkb:pkt) 2662 ! WRITE(9, 5555) 'tend :', tend(pkb:pkt,j,i) 2663 ! WRITE(9, 5555) 'diss_p:', diss_p(pkb:pkt,j,i) 2664 ! WRITE(9, 5555) 'diss :', diss(pkb:pkt,j,i) 2665 ! WRITE(9, *) '---' 2666 ! WRITE(9, 5555) 'km :', km(pkb:pkt,j,i) 2667 ! flush(9) 2668 ! endif 2371 2669 2372 2670 ENDIF ! dissipation equation … … 2382 2680 !> @warning The case with constant_flux_layer = F and use_surface_fluxes = T is 2383 2681 !> not considered well! 2682 !> @todo Adjust production term in case of rans_tke_e simulation 2384 2683 !------------------------------------------------------------------------------! 2385 2684 SUBROUTINE production_e … … 3099 3398 !> @warning The case with constant_flux_layer = F and use_surface_fluxes = T is 3100 3399 !> not considered well! 3400 !> @todo non-neutral case is not yet considered for RANS mode 3101 3401 !------------------------------------------------------------------------------! 3102 SUBROUTINE production_e_ij( i, j )3402 SUBROUTINE production_e_ij( i, j, diss_production ) 3103 3403 3104 3404 USE arrays_3d, & … … 3121 3421 3122 3422 IMPLICIT NONE 3423 3424 LOGICAL :: diss_production 3123 3425 3124 3426 INTEGER(iwp) :: i !< running index x-direction … … 3153 3455 REAL(wp), DIMENSION(nzb+1:nzt) :: dwdy !< Gradient of w-component in y-direction 3154 3456 REAL(wp), DIMENSION(nzb+1:nzt) :: dwdz !< Gradient of w-component in z-direction 3155 3457 REAL(wp), DIMENSION(nzb+1:nzt) :: tend_temp !< temporal tendency 3156 3458 3157 3459 IF ( constant_flux_layer ) THEN … … 3348 3650 ENDDO 3349 3651 3652 ! IF ( .NOT. rans_tke_e ) THEN 3653 3350 3654 DO k = nzb+1, nzt 3351 3655 3352 def = 2.0_wp * ( dudx(k)**2 + dvdy(k)**2 + dwdz(k)**2 ) + & 3353 dudy(k)**2 + dvdx(k)**2 + dwdx(k)**2 + & 3354 dwdy(k)**2 + dudz(k)**2 + dvdz(k)**2 + & 3355 2.0_wp * ( dvdx(k)*dudy(k) + dwdx(k)*dudz(k) + dwdy(k)*dvdz(k) ) 3356 3357 ! 3358 !-- Production term according to Kato and Launder (1993) 3359 ! def = SQRT( ( dudx(k) + dudy(k) + dudz(k) + & 3360 ! dvdx(k) + dvdy(k) + dvdz(k) + & 3361 ! dwdx(k) + dwdy(k) + dwdz(k) & 3362 ! )**4 - & 3363 ! ( dudx(k)**2 + dvdy(k)**2 + dwdz(k)**2 + & 3364 ! 2.0_wp * ( dudy(k) * dvdx(k) + & 3365 ! dudz(k) * dwdx(k) + & 3366 ! dvdz(k) * dwdy(k) ) & 3367 ! )**2 & 3368 ! ) 3656 def = 2.0_wp * ( dudx(k)**2 + dvdy(k)**2 + dwdz(k)**2 ) + & 3657 dudy(k)**2 + dvdx(k)**2 + dwdx(k)**2 + & 3658 dwdy(k)**2 + dudz(k)**2 + dvdz(k)**2 + & 3659 2.0_wp * ( dvdx(k)*dudy(k) + & 3660 dwdx(k)*dudz(k) + & 3661 dwdy(k)*dvdz(k) ) 3369 3662 3370 3663 IF ( def < 0.0_wp ) def = 0.0_wp … … 3372 3665 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) ) 3373 3666 3374 tend (k,j,i) = tend(k,j,i) +km(k,j,i) * def * flag3667 tend_temp(k) = km(k,j,i) * def * flag 3375 3668 3376 3669 ENDDO 3377 3670 3671 ! ELSE 3672 ! 3673 ! DO k = nzb+1, nzt 3674 ! ! 3675 ! !-- Production term according to Kato and Launder (1993) 3676 ! def = SQRT( ( dudy(k)**2 + dvdz(k)**2 + dwdx(k)**2 + & 3677 ! dudz(k)**2 + dvdx(k)**2 + dwdy(k)**2 + & 3678 ! 2.0_wp * ( dudy(k) * dvdx(k) + & 3679 ! dvdz(k) * dwdy(k) + & 3680 ! dwdx(k) * dudz(k) ) ) & 3681 ! * ( dudy(k)**2 + dvdz(k)**2 + dwdx(k)**2 + & 3682 ! dudz(k)**2 + dvdx(k)**2 + dwdy(k)**2 - & 3683 ! 2.0_wp * ( dudy(k) * dvdx(k) + & 3684 ! dvdz(k) * dwdy(k) + & 3685 ! dwdx(k) * dudz(k) ) ) & 3686 ! ) 3687 ! 3688 ! IF ( def < 0.0_wp ) def = 0.0_wp 3689 ! 3690 ! flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) ) 3691 ! 3692 ! tend_temp(k) = km(k,j,i) * def * flag 3693 ! 3694 ! ENDDO 3695 ! 3696 ! ENDIF 3697 3698 ELSE ! not constant_flux_layer 3699 3700 ! IF ( .NOT. rans_tke_e ) THEN 3701 ! 3702 !-- Calculate TKE production by shear. Here, no additional 3703 !-- wall-bounded code is considered. 3704 !-- Why? 3705 DO k = nzb+1, nzt 3706 3707 dudx(k) = ( u(k,j,i+1) - u(k,j,i) ) * ddx 3708 dudy(k) = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - & 3709 u(k,j-1,i) - u(k,j-1,i+1) ) * ddy 3710 dudz(k) = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - & 3711 u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k) 3712 3713 dvdx(k) = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - & 3714 v(k,j,i-1) - v(k,j+1,i-1) ) * ddx 3715 dvdy(k) = ( v(k,j+1,i) - v(k,j,i) ) * ddy 3716 dvdz(k) = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - & 3717 v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k) 3718 3719 dwdx(k) = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - & 3720 w(k,j,i-1) - w(k-1,j,i-1) ) * ddx 3721 dwdy(k) = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - & 3722 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy 3723 dwdz(k) = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) 3724 3725 def = 2.0_wp * ( dudx(k)**2 + dvdy(k)**2 + dwdz(k)**2 ) + & 3726 dudy(k)**2 + dvdx(k)**2 + dwdx(k)**2 + & 3727 dwdy(k)**2 + dudz(k)**2 + dvdz(k)**2 + & 3728 2.0_wp * ( dvdx(k)*dudy(k) + & 3729 dwdx(k)*dudz(k) + & 3730 dwdy(k)*dvdz(k) ) 3731 3732 IF ( def < 0.0_wp ) def = 0.0_wp 3733 3734 flag = MERGE( 1.0_wp, 0.0_wp, & 3735 BTEST( wall_flags_0(k,j,i), 29 ) ) 3736 tend_temp(k) = km(k,j,i) * def * flag 3737 3738 ENDDO 3739 3740 ! ELSE 3741 ! 3742 ! DO k = nzb+1, nzt 3743 ! 3744 ! dudx(k) = ( u(k,j,i+1) - u(k,j,i) ) * ddx 3745 ! dudy(k) = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - & 3746 ! u(k,j-1,i) - u(k,j-1,i+1) ) * ddy 3747 ! dudz(k) = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - & 3748 ! u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k) 3749 ! 3750 ! dvdx(k) = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - & 3751 ! v(k,j,i-1) - v(k,j+1,i-1) ) * ddx 3752 ! dvdy(k) = ( v(k,j+1,i) - v(k,j,i) ) * ddy 3753 ! dvdz(k) = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - & 3754 ! v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k) 3755 ! 3756 ! dwdx(k) = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - & 3757 ! w(k,j,i-1) - w(k-1,j,i-1) ) * ddx 3758 ! dwdy(k) = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - & 3759 ! w(k,j-1,i) - w(k-1,j-1,i) ) * ddy 3760 ! dwdz(k) = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) 3761 ! ! 3762 ! !-- Production term according to Kato and Launder (1993) 3763 ! def = SQRT( ( dudy(k)**2 + dvdz(k)**2 + dwdx(k)**2 + & 3764 ! dudz(k)**2 + dvdx(k)**2 + dwdy(k)**2 + & 3765 ! 2.0_wp * ( dudy(k) * dvdx(k) + & 3766 ! dvdz(k) * dwdy(k) + & 3767 ! dwdx(k) * dudz(k) ) ) & 3768 ! * ( dudy(k)**2 + dvdz(k)**2 + dwdx(k)**2 + & 3769 ! dudz(k)**2 + dvdx(k)**2 + dwdy(k)**2 - & 3770 ! 2.0_wp * ( dudy(k) * dvdx(k) + & 3771 ! dvdz(k) * dwdy(k) + & 3772 ! dwdx(k) * dudz(k) ) ) & 3773 ! ) 3774 ! 3775 ! IF ( def < 0.0_wp ) def = 0.0_wp 3776 ! 3777 ! flag = MERGE( 1.0_wp, 0.0_wp, & 3778 ! BTEST( wall_flags_0(k,j,i), 29 ) ) 3779 ! tend_temp(k) = km(k,j,i) * def * flag 3780 ! 3781 ! ENDDO 3782 ! 3783 ! ENDIF 3784 3785 ENDIF 3786 3787 IF ( .NOT. diss_production ) THEN 3788 ! 3789 !-- Production term in case of TKE production 3790 DO k = nzb+1, nzt 3791 tend(k,j,i) = tend(k,j,i) + tend_temp(k) 3792 ENDDO 3378 3793 ELSE 3379 3794 ! 3380 !-- Calculate TKE production by shear. Here, no additional 3381 !-- wall-bounded code is considered. 3382 !-- Why? 3795 !-- Production term in case of dissipation-rate production (rans_tke_e) 3383 3796 DO k = nzb+1, nzt 3384 3797 3385 dudx(k) = ( u(k,j,i+1) - u(k,j,i) ) * ddx 3386 dudy(k) = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - & 3387 u(k,j-1,i) - u(k,j-1,i+1) ) * ddy 3388 dudz(k) = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - & 3389 u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k) 3390 3391 dvdx(k) = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - & 3392 v(k,j,i-1) - v(k,j+1,i-1) ) * ddx 3393 dvdy(k) = ( v(k,j+1,i) - v(k,j,i) ) * ddy 3394 dvdz(k) = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - & 3395 v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k) 3396 3397 dwdx(k) = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - & 3398 w(k,j,i-1) - w(k-1,j,i-1) ) * ddx 3399 dwdy(k) = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - & 3400 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy 3401 dwdz(k) = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) 3402 3403 def = 2.0_wp * ( dudx(k)**2 + dvdy(k)**2 + dwdz(k)**2 ) + & 3404 dudy(k)**2 + dvdx(k)**2 + dwdx(k)**2 + & 3405 dwdy(k)**2 + dudz(k)**2 + dvdz(k)**2 + & 3406 2.0_wp * ( dvdx(k)*dudy(k) + dwdx(k)*dudz(k) + dwdy(k)*dvdz(k) ) 3407 3408 ! 3409 !-- Production term according to Kato and Launder (1993) 3410 ! def = SQRT( ( dudx(k) + dudy(k) + dudz(k) + & 3411 ! dvdx(k) + dvdy(k) + dvdz(k) + & 3412 ! dwdx(k) + dwdy(k) + dwdz(k) & 3413 ! )**4 - & 3414 ! ( dudx(k)**2 + dvdy(k)**2 + dwdz(k)**2 + & 3415 ! 2.0_wp * ( dudy(k) * dvdx(k) + & 3416 ! dudz(k) * dwdx(k) + & 3417 ! dvdz(k) * dwdy(k) ) & 3418 ! )**2 & 3419 ! ) 3420 3421 IF ( def < 0.0_wp ) def = 0.0_wp 3422 3423 flag = MERGE( 1.0_wp, 0.0_wp, & 3424 BTEST( wall_flags_0(k,j,i), 29 ) ) 3425 tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def * flag 3426 3798 ! Standard TKE-e closure 3799 tend(k,j,i) = tend(k,j,i) + tend_temp(k) * diss(k,j,i) & 3800 /( e(k,j,i) + 1.0E-20_wp ) & 3801 * c_1 3802 ! ! Production according to Koblitz (2013) 3803 ! tend(k,j,i) = tend(k,j,i) + tend_temp(k) * diss(k,j,i) & 3804 ! /( e(k,j,i) + 1.0E-20_wp ) & 3805 ! * ( c_1 + ( c_2 - c_1 ) & 3806 ! * l_wall(k,j,i) / l_max ) 3807 ! ! Production according to Detering and Etling (1985) 3808 ! !> @todo us is not correct if there are vertical walls 3809 ! tend(k,j,i) = tend(k,j,i) + tend_temp(k) * SQRT(e(k,j,i)) & 3810 ! * c_1 * c_0**3 / c_4 * f & 3811 ! / surf_def_h(0)%us(surf_def_h(0)%start_index(j,i)) 3427 3812 ENDDO 3428 3429 3813 ENDIF 3430 3814 … … 3822 4206 REAL(wp) :: l !< mixing length 3823 4207 REAL(wp) :: ll !< adjusted l 3824 REAL(wp) :: var_reference !< 4208 REAL(wp) :: var_reference !< reference temperature 3825 4209 3826 4210 #if defined( __nopointer ) 3827 REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: var !< 4211 REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: var !< temperature 3828 4212 #else 3829 REAL(wp), DIMENSION(:,:,:), POINTER :: var !< 4213 REAL(wp), DIMENSION(:,:,:), POINTER :: var !< temperature 3830 4214 #endif 3831 4215 REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) :: dissipation !< TKE dissipation … … 3837 4221 DO j = nys, nyn 3838 4222 DO k = nzb+1, nzt 4223 ! 4224 !-- Predetermine flag to mask topography 4225 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) ) 3839 4226 3840 4227 ! … … 3851 4238 CALL mixing_length_rans( i, j, k, l, ll, var, var_reference ) 3852 4239 3853 dissipation(k,j) = c_m**3 * e(k,j,i) * SQRT( e(k,j,i) ) / ll 4240 dissipation(k,j) = c_0**3 * e(k,j,i) * SQRT( e(k,j,i) ) / ll 4241 4242 diss(k,j,i) = dissipation(k,j) * flag 3854 4243 3855 4244 ELSEIF ( rans_tke_e ) THEN … … 3859 4248 ENDIF 3860 4249 3861 ! 3862 !-- Predetermine flag to mask topography 3863 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) ) 3864 3865 tend(k,j,i) = tend(k,j,i) & 3866 + ( & 4250 tend(k,j,i) = tend(k,j,i) + ( & 4251 ( & 3867 4252 ( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)-e(k,j,i) ) & 3868 4253 - ( km(k,j,i)+km(k,j,i-1) ) * ( e(k,j,i)-e(k,j,i-1) ) & 3869 ) * ddx2 * flag&3870 + (&4254 ) * ddx2 * flag & 4255 + ( & 3871 4256 ( km(k,j,i)+km(k,j+1,i) ) * ( e(k,j+1,i)-e(k,j,i) ) & 3872 4257 - ( km(k,j,i)+km(k,j-1,i) ) * ( e(k,j,i)-e(k,j-1,i) ) & 3873 ) * ddy2 * flag&3874 + (&4258 ) * ddy2 * flag & 4259 + ( & 3875 4260 ( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1) & 3876 4261 * rho_air_zw(k) & 3877 4262 - ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k) & 3878 4263 * rho_air_zw(k-1) & 3879 ) * ddzw(k) * drho_air(k) * flag & 4264 ) * ddzw(k) * drho_air(k) & 4265 ) * flag * dsig_e & 3880 4266 - dissipation(k,j) * flag 3881 4267 … … 3959 4345 REAL(wp) :: l !< mixing length 3960 4346 REAL(wp) :: ll !< adjusted l 3961 REAL(wp) :: var_reference !< 4347 REAL(wp) :: var_reference !< reference temperature 3962 4348 3963 4349 #if defined( __nopointer ) 3964 REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: var !< 4350 REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: var !< temperature 3965 4351 #else 3966 REAL(wp), DIMENSION(:,:,:), POINTER :: var !< 4352 REAL(wp), DIMENSION(:,:,:), POINTER :: var !< temperature 3967 4353 #endif 3968 4354 REAL(wp), DIMENSION(nzb+1:nzt) :: dissipation !< dissipation of TKE … … 3991 4377 CALL mixing_length_rans( i, j, k, l, ll, var, var_reference ) 3992 4378 3993 dissipation(k) = c_m**3 * e(k,j,i) * SQRT( e(k,j,i) ) / ll 4379 dissipation(k) = c_0**3 * e(k,j,i) * SQRT( e(k,j,i) ) / ll 4380 4381 diss(k,j,i) = dissipation(k) * flag 3994 4382 3995 4383 ELSEIF ( rans_tke_e ) THEN … … 4001 4389 ! 4002 4390 !-- Calculate the tendency term 4003 tend(k,j,i) = tend(k,j,i) 4004 +( &4391 tend(k,j,i) = tend(k,j,i) + ( & 4392 ( & 4005 4393 ( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)-e(k,j,i) ) & 4006 4394 - ( km(k,j,i)+km(k,j,i-1) ) * ( e(k,j,i)-e(k,j,i-1) ) & 4007 ) * ddx2 * flag / sig_e&4395 ) * ddx2 & 4008 4396 + ( & 4009 4397 ( km(k,j,i)+km(k,j+1,i) ) * ( e(k,j+1,i)-e(k,j,i) ) & 4010 4398 - ( km(k,j,i)+km(k,j-1,i) ) * ( e(k,j,i)-e(k,j-1,i) ) & 4011 ) * ddy2 * flag / sig_e&4399 ) * ddy2 & 4012 4400 + ( & 4013 4401 ( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1) & … … 4015 4403 - ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k) & 4016 4404 * rho_air_zw(k-1) & 4017 ) * ddzw(k) * drho_air(k) * flag / sig_e & 4018 - dissipation(k) * flag 4405 ) * ddzw(k) * drho_air(k) & 4406 ) * flag * dsig_e & 4407 - dissipation(k) * flag 4019 4408 4020 4409 ENDDO … … 4082 4471 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) ) 4083 4472 4084 tend(k,j,i) = tend(k,j,i) 4085 +( &4473 tend(k,j,i) = tend(k,j,i) + & 4474 ( ( & 4086 4475 ( km(k,j,i)+km(k,j,i+1) ) * ( diss(k,j,i+1)-diss(k,j,i) ) & 4087 4476 - ( km(k,j,i)+km(k,j,i-1) ) * ( diss(k,j,i)-diss(k,j,i-1) ) & 4088 ) * ddx2 * flag&4477 ) * ddx2 & 4089 4478 + ( & 4090 4479 ( km(k,j,i)+km(k,j+1,i) ) * ( diss(k,j+1,i)-diss(k,j,i) ) & 4091 4480 - ( km(k,j,i)+km(k,j-1,i) ) * ( diss(k,j,i)-diss(k,j-1,i) ) & 4092 ) * ddy2 * flag&4481 ) * ddy2 & 4093 4482 + ( & 4094 4483 ( km(k,j,i)+km(k+1,j,i) ) * ( diss(k+1,j,i)-diss(k,j,i) ) * ddzu(k+1) & … … 4096 4485 - ( km(k,j,i)+km(k-1,j,i) ) * ( diss(k,j,i)-diss(k-1,j,i) ) * ddzu(k) & 4097 4486 * rho_air_zw(k-1) & 4098 ) * ddzw(k) * drho_air(k) * flag & 4099 - c_2 * diss(k,j,i)**2 & 4100 / ( e(k,j,i) + 1.0E-20_wp ) * flag 4487 ) * ddzw(k) * drho_air(k) & 4488 ) * flag * dsig_diss & 4489 - c_2 * diss(k,j,i)**2 & 4490 / ( e(k,j,i) + 1.0E-20_wp ) * flag 4101 4491 4102 4492 ENDDO … … 4123 4513 IMPLICIT NONE 4124 4514 4125 INTEGER(iwp) :: i !< running index x direction 4126 INTEGER(iwp) :: j !< running index y direction 4127 INTEGER(iwp) :: k !< running index z direction 4128 4129 REAL(wp) :: flag !< flag to mask topography 4130 4131 REAL(wp), DIMENSION(nzb+1:nzt) :: tend_temp 4515 INTEGER(iwp) :: i !< running index x direction 4516 INTEGER(iwp) :: j !< running index y direction 4517 INTEGER(iwp) :: k !< running index z direction 4518 4519 REAL(wp) :: flag !< flag to mask topography 4132 4520 4133 4521 ! … … 4141 4529 ! 4142 4530 !-- Calculate the tendency term 4143 tend_temp(k) = ( & 4531 tend(k,j,i) = tend(k,j,i) + & 4532 ( ( & 4144 4533 ( km(k,j,i)+km(k,j,i+1) ) * ( diss(k,j,i+1)-diss(k,j,i) ) & 4145 4534 - ( km(k,j,i)+km(k,j,i-1) ) * ( diss(k,j,i)-diss(k,j,i-1) ) & 4146 ) * ddx2 * flag / sig_diss&4535 ) * ddx2 & 4147 4536 + ( & 4148 4537 ( km(k,j,i)+km(k,j+1,i) ) * ( diss(k,j+1,i)-diss(k,j,i) ) & 4149 4538 - ( km(k,j,i)+km(k,j-1,i) ) * ( diss(k,j,i)-diss(k,j-1,i) ) & 4150 ) * ddy2 * flag / sig_diss&4539 ) * ddy2 & 4151 4540 + ( & 4152 4541 ( km(k,j,i)+km(k+1,j,i) ) * ( diss(k+1,j,i)-diss(k,j,i) ) * ddzu(k+1) & … … 4154 4543 - ( km(k,j,i)+km(k-1,j,i) ) * ( diss(k,j,i)-diss(k-1,j,i) ) * ddzu(k) & 4155 4544 * rho_air_zw(k-1) & 4156 ) * ddzw(k) * drho_air(k) * flag / sig_diss & 4157 - c_2 * diss(k,j,i)**2 & 4158 / ( e(k,j,i) + 1.0E-20_wp ) * flag 4159 4160 tend(k,j,i) = tend(k,j,i) + tend_temp(k) 4545 ) * ddzw(k) * drho_air(k) & 4546 ) * flag * dsig_diss & 4547 - c_2 * diss(k,j,i)**2 / ( e(k,j,i) + 1.0E-20_wp ) * flag 4161 4548 4162 4549 ENDDO … … 4290 4677 !> Computation of the turbulent diffusion coefficients for momentum and heat 4291 4678 !> according to Prandtl-Kolmogorov. 4679 !> @todo consider non-default surfaces 4292 4680 !------------------------------------------------------------------------------! 4293 4681 SUBROUTINE tcm_diffusivities( var, var_reference ) … … 4297 4685 ONLY: e_min, outflow_l, outflow_n, outflow_r, outflow_s 4298 4686 4687 USE grid_variables, & 4688 ONLY: dx, dy 4689 4299 4690 USE statistics, & 4300 4691 ONLY : rmask, sums_l_l 4301 4692 4302 4693 USE surface_mod, & 4303 ONLY : bc_h, surf_def_h 4694 ONLY : bc_h, surf_def_h, surf_def_v 4304 4695 4305 4696 IMPLICIT NONE 4306 4697 4307 INTEGER(iwp) :: i !< 4308 INTEGER(iwp) :: j !< 4309 INTEGER(iwp) :: k !< 4310 INTEGER(iwp) :: m !< 4311 INTEGER(iwp) :: n !< 4312 INTEGER(iwp) :: omp_get_thread_num !< 4313 INTEGER(iwp) :: sr !< 4314 INTEGER(iwp) :: tn !< 4315 4316 REAL(wp) :: flag !< 4317 REAL(wp) :: l !< 4318 REAL(wp) :: ll !< 4319 REAL(wp) :: var_reference !< 4698 INTEGER(iwp) :: i !< loop index 4699 INTEGER(iwp) :: j !< loop index 4700 INTEGER(iwp) :: k !< loop index 4701 INTEGER(iwp) :: m !< loop index 4702 INTEGER(iwp) :: n !< loop index 4703 INTEGER(iwp) :: omp_get_thread_num !< opemmp function to get thread number 4704 INTEGER(iwp) :: sr !< statistic region 4705 INTEGER(iwp) :: tn !< thread number 4706 4707 REAL(wp) :: flag !< topography flag 4708 REAL(wp) :: l !< mixing length 4709 REAL(wp) :: ll !< adjusted mixing length 4710 REAL(wp) :: var_reference !< reference temperature 4320 4711 4321 4712 #if defined( __nopointer ) 4322 REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: var !< 4713 REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: var !< temperature 4323 4714 #else 4324 REAL(wp), DIMENSION(:,:,:), POINTER :: var !< 4715 REAL(wp), DIMENSION(:,:,:), POINTER :: var !< temperature 4325 4716 #endif 4326 4717 … … 4365 4756 ! 4366 4757 !-- Compute diffusion coefficients for momentum and heat 4367 km(k,j,i) = c_ m* l * SQRT( e(k,j,i) ) * flag4758 km(k,j,i) = c_0 * l * SQRT( e(k,j,i) ) * flag 4368 4759 kh(k,j,i) = ( 1.0_wp + 2.0_wp * l / ll ) * km(k,j,i) * flag 4369 4760 ! … … 4391 4782 ! 4392 4783 !-- Compute diffusion coefficients for momentum and heat 4393 km(k,j,i) = c_ m* l * SQRT( e(k,j,i) ) * flag4784 km(k,j,i) = c_0 * l * SQRT( e(k,j,i) ) * flag 4394 4785 kh(k,j,i) = km(k,j,i) / prandtl_number * flag 4395 4786 ! … … 4414 4805 ! 4415 4806 !-- Compute diffusion coefficients for momentum and heat 4416 km(k,j,i) = c_ mu * e(k,j,i)**2 / ( diss(k,j,i) + 1.E-10) * flag4807 km(k,j,i) = c_0**4 * e(k,j,i)**2 / ( diss(k,j,i) + 1.0E-30_wp ) * flag 4417 4808 kh(k,j,i) = km(k,j,i) / prandtl_number * flag 4809 ! 4810 !-- Summation for averaged profile of mixing length (cf. flow_statistics) 4811 DO sr = 0, statistic_regions 4812 sums_l_l(k,sr,tn) = sums_l_l(k,sr,tn) + & 4813 c_0**3 * e(k,j,i) * SQRT(e(k,j,i)) / & 4814 ( diss(k,j,i) + 1.0E-30_wp ) * rmask(j,i,sr) * flag 4815 ENDDO 4418 4816 4419 4817 ENDDO … … 4434 4832 !-- so far vertical surfaces require usage of a Prandtl-layer where the boundary 4435 4833 !-- values of the diffusivities are not needed. 4436 4437 4834 IF ( .NOT. rans_tke_e ) THEN 4438 4835 ! … … 4468 4865 ENDDO 4469 4866 ENDDO 4867 ! 4868 !-- North- and southward facing surfaces 4869 DO n = 0, 1 4870 DO m = 1, surf_def_v(n)%ns 4871 i = surf_def_v(n)%i(m) 4872 j = surf_def_v(n)%j(m) 4873 k = surf_def_v(n)%k(m) 4874 km(k,j,i) = kappa * surf_def_v(n)%us(m) * 0.5_wp * dy 4875 kh(k,j,i) = 1.35_wp * km(k,j,i) 4876 ENDDO 4877 ENDDO 4878 ! 4879 !-- West- and eastward facing surfaces 4880 DO n = 2, 3 4881 DO m = 1, surf_def_v(n)%ns 4882 i = surf_def_v(n)%i(m) 4883 j = surf_def_v(n)%j(m) 4884 k = surf_def_v(n)%k(m) 4885 km(k,j,i) = kappa * surf_def_v(n)%us(m) * 0.5_wp * dx 4886 kh(k,j,i) = 1.35_wp * km(k,j,i) 4887 ENDDO 4888 ENDDO 4470 4889 4471 4890 CALL exchange_horiz( km, nbgp ) … … 4473 4892 4474 4893 ENDIF 4894 4475 4895 ! 4476 4896 !-- Model top … … 4518 4938 INTEGER(iwp) :: j !< loop index y direction 4519 4939 INTEGER(iwp) :: k !< loop index z direction 4520 INTEGER, INTENT(IN) :: mod_count !< 4940 INTEGER, INTENT(IN) :: mod_count !< flag defining where pointers point to 4521 4941 4522 4942 #if defined( __nopointer )
Note: See TracChangeset
for help on using the changeset viewer.
