Changeset 3129 for palm/trunk/SOURCE/turbulence_closure_mod.f90
- Timestamp:
- Jul 16, 2018 7:45:13 AM (6 years ago)
- Location:
- palm/trunk
- Files:
-
- 3 edited
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk
- Property svn:mergeinfo changed
/palm/branches/rans merged: 3087,3100-3102,3128
- Property svn:mergeinfo changed
-
palm/trunk/SOURCE
- Property svn:mergeinfo changed
/palm/branches/rans/SOURCE merged: 3087,3100-3102,3128
- Property svn:mergeinfo changed
-
palm/trunk/SOURCE/turbulence_closure_mod.f90
r3121 r3129 25 25 ! ----------------- 26 26 ! $Id$ 27 ! - move limitation of diss to boundary_conds 28 ! - move boundary conditions for e and diss to boundary_conds 29 ! - consider non-default surfaces in tcm_diffusivities 30 ! - use z_mo within surface layer instead of calculating it 31 ! - resort output after case select -> reduced code duplication 32 ! - when using rans_tke_e and 1d-model, do not use e1d, km1d and diss1d 33 ! 34 ! 3121 2018-07-11 18:46:49Z gronemeier 27 35 ! - created the function phi_m 28 36 ! - implemented km = u* * kappa * zp / phi_m in production_e_init for all … … 117 125 !> @todo Check for random disturbances 118 126 !> @note <Enter notes on the module> 119 !> @bug TKE-e closure still crashes due to too small dt120 127 !------------------------------------------------------------------------------! 121 128 MODULE turbulence_closure_mod … … 202 209 203 210 !> @todo remove debug variables 204 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: diss_prod1, diss_adve1, diss_diff1, & 205 diss_prod2, diss_adve2, diss_diff2, & 206 diss_prod3, diss_adve3, diss_diff3, & 207 dummy1, dummy2, dummy3 211 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: & 212 diss_prod1, diss_adve1, diss_diff1, & 213 diss_prod2, diss_adve2, diss_diff2, & 214 diss_prod3, diss_adve3, diss_diff3, & 215 dummy1, dummy2, dummy3 208 216 209 217 … … 689 697 CHARACTER (LEN=*) :: variable !< name of variable 690 698 691 INTEGER(iwp) :: av !< flag for (non-)average output 692 INTEGER(iwp) :: i !< loop index 693 INTEGER(iwp) :: j !< loop index 694 INTEGER(iwp) :: k !< loop index 695 INTEGER(iwp) :: nzb_do !< vertical output index (bottom) 696 INTEGER(iwp) :: nzt_do !< vertical output index (top) 699 INTEGER(iwp) :: av !< flag for (non-)average output 700 INTEGER(iwp) :: flag_nr !< number of masking flag 701 INTEGER(iwp) :: i !< loop index 702 INTEGER(iwp) :: j !< loop index 703 INTEGER(iwp) :: k !< loop index 704 INTEGER(iwp) :: nzb_do !< vertical output index (bottom) 705 INTEGER(iwp) :: nzt_do !< vertical output index (top) 697 706 698 707 LOGICAL :: found !< flag if output variable is found 699 708 LOGICAL :: two_d !< flag parameter that indicates 2D variables (horizontal cross sections) 709 LOGICAL :: resorted !< flag if output is already resorted 700 710 701 711 REAL(wp) :: fill_value = -999.0_wp !< value for the _FillValue attribute … … 704 714 !< array to which output data is resorted to 705 715 716 REAL(wp), DIMENSION(:,:,:), POINTER :: to_be_resorted !< points to selected output variable 717 706 718 found = .TRUE. 719 resorted = .FALSE. 720 ! 721 !-- Set masking flag for topography for not resorted arrays 722 flag_nr = 0 707 723 708 724 SELECT CASE ( TRIM( variable ) ) 709 710 725 711 726 CASE ( 'diss_xy', 'diss_xz', 'diss_yz' ) 712 727 IF ( av == 0 ) THEN 713 DO i = nxl, nxr 714 DO j = nys, nyn 715 DO k = nzb_do, nzt_do 716 local_pf(i,j,k) = diss(k,j,i) 717 ENDDO 718 ENDDO 719 ENDDO 728 to_be_resorted => diss 720 729 ELSE 721 730 IF ( .NOT. ALLOCATED( diss_av ) ) THEN … … 723 732 diss_av = REAL( fill_value, KIND = wp ) 724 733 ENDIF 725 DO i = nxl, nxr 726 DO j = nys, nyn 727 DO k = nzb_do, nzt_do 728 local_pf(i,j,k) = diss_av(k,j,i) 729 ENDDO 730 ENDDO 731 ENDDO 734 to_be_resorted => diss_av 732 735 ENDIF 733 734 736 IF ( mode == 'xy' ) grid = 'zu' 735 737 736 738 CASE ( 'kh_xy', 'kh_xz', 'kh_yz' ) 737 739 IF ( av == 0 ) THEN 738 DO i = nxl, nxr 739 DO j = nys, nyn 740 DO k = nzb_do, nzt_do 741 local_pf(i,j,k) = kh(k,j,i) 742 ENDDO 743 ENDDO 740 to_be_resorted => kh 741 ELSE 742 IF ( .NOT. ALLOCATED( kh_av ) ) THEN 743 ALLOCATE( kh_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 744 kh_av = REAL( fill_value, KIND = wp ) 745 ENDIF 746 to_be_resorted => kh_av 747 ENDIF 748 IF ( mode == 'xy' ) grid = 'zu' 749 750 CASE ( 'km_xy', 'km_xz', 'km_yz' ) 751 IF ( av == 0 ) THEN 752 to_be_resorted => km 753 ELSE 754 IF ( .NOT. ALLOCATED( km_av ) ) THEN 755 ALLOCATE( km_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 756 km_av = REAL( fill_value, KIND = wp ) 757 ENDIF 758 to_be_resorted => km_av 759 ENDIF 760 IF ( mode == 'xy' ) grid = 'zu' 761 762 CASE DEFAULT 763 found = .FALSE. 764 grid = 'none' 765 766 END SELECT 767 768 IF ( found .AND. .NOT. resorted ) THEN 769 DO i = nxl, nxr 770 DO j = nys, nyn 771 DO k = nzb_do, nzt_do 772 local_pf(i,j,k) = MERGE( to_be_resorted(k,j,i), & 773 REAL( fill_value, KIND = wp ), & 774 BTEST( wall_flags_0(k,j,i), flag_nr ) ) 744 775 ENDDO 776 ENDDO 777 ENDDO 778 ENDIF 779 780 END SUBROUTINE tcm_data_output_2d 781 782 783 !------------------------------------------------------------------------------! 784 ! Description: 785 ! ------------ 786 !> Define 3D output variables. 787 !------------------------------------------------------------------------------! 788 SUBROUTINE tcm_data_output_3d( av, variable, found, local_pf, nzb_do, nzt_do ) 789 790 791 USE averaging, & 792 ONLY: diss_av, kh_av, km_av 793 794 IMPLICIT NONE 795 796 CHARACTER (LEN=*) :: variable !< name of variable 797 798 INTEGER(iwp) :: av !< flag for (non-)average output 799 INTEGER(iwp) :: flag_nr !< number of masking flag 800 INTEGER(iwp) :: i !< loop index 801 INTEGER(iwp) :: j !< loop index 802 INTEGER(iwp) :: k !< loop index 803 INTEGER(iwp) :: nzb_do !< lower limit of the data output (usually 0) 804 INTEGER(iwp) :: nzt_do !< vertical upper limit of the data output (usually nz_do3d) 805 806 LOGICAL :: found !< flag if output variable is found 807 LOGICAL :: resorted !< flag if output is already resorted 808 809 REAL(wp) :: fill_value = -999.0_wp !< value for the _FillValue attribute 810 811 REAL(sp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) :: local_pf !< local 812 !< array to which output data is resorted to 813 814 REAL(wp), DIMENSION(:,:,:), POINTER :: to_be_resorted !< points to selected output variable 815 816 found = .TRUE. 817 resorted = .FALSE. 818 ! 819 !-- Set masking flag for topography for not resorted arrays 820 flag_nr = 0 821 822 SELECT CASE ( TRIM( variable ) ) 823 824 CASE ( 'diss' ) 825 IF ( av == 0 ) THEN 826 to_be_resorted => diss 745 827 ELSE 746 828 IF ( .NOT. ALLOCATED( diss_av ) ) THEN … … 748 830 diss_av = REAL( fill_value, KIND = wp ) 749 831 ENDIF 750 DO i = nxl, nxr 751 DO j = nys, nyn 752 DO k = nzb_do, nzt_do 753 local_pf(i,j,k) = kh_av(k,j,i) 754 ENDDO 755 ENDDO 756 ENDDO 757 ENDIF 758 759 IF ( mode == 'xy' ) grid = 'zu' 760 761 CASE ( 'km_xy', 'km_xz', 'km_yz' ) 762 IF ( av == 0 ) THEN 763 DO i = nxl, nxr 764 DO j = nys, nyn 765 DO k = nzb_do, nzt_do 766 local_pf(i,j,k) = km(k,j,i) 767 ENDDO 768 ENDDO 769 ENDDO 770 ELSE 771 IF ( .NOT. ALLOCATED( diss_av ) ) THEN 772 ALLOCATE( diss_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 773 diss_av = REAL( fill_value, KIND = wp ) 774 ENDIF 775 DO i = nxl, nxr 776 DO j = nys, nyn 777 DO k = nzb_do, nzt_do 778 local_pf(i,j,k) = km_av(k,j,i) 779 ENDDO 780 ENDDO 781 ENDDO 782 ENDIF 783 784 IF ( mode == 'xy' ) grid = 'zu' 785 786 CASE DEFAULT 787 found = .FALSE. 788 grid = 'none' 789 790 END SELECT 791 792 END SUBROUTINE tcm_data_output_2d 793 794 795 !------------------------------------------------------------------------------! 796 ! Description: 797 ! ------------ 798 !> Define 3D output variables. 799 !------------------------------------------------------------------------------! 800 SUBROUTINE tcm_data_output_3d( av, variable, found, local_pf, nzb_do, nzt_do ) 801 802 803 USE averaging, & 804 ONLY: diss_av, kh_av, km_av 805 806 IMPLICIT NONE 807 808 CHARACTER (LEN=*) :: variable !< name of variable 809 810 INTEGER(iwp) :: av !< flag for (non-)average output 811 INTEGER(iwp) :: i !< loop index 812 INTEGER(iwp) :: j !< loop index 813 INTEGER(iwp) :: k !< loop index 814 INTEGER(iwp) :: nzb_do !< lower limit of the data output (usually 0) 815 INTEGER(iwp) :: nzt_do !< vertical upper limit of the data output (usually nz_do3d) 816 817 LOGICAL :: found !< flag if output variable is found 818 819 REAL(wp) :: fill_value = -999.0_wp !< value for the _FillValue attribute 820 821 REAL(sp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) :: local_pf !< local 822 !< array to which output data is resorted to 823 824 825 found = .TRUE. 826 827 828 SELECT CASE ( TRIM( variable ) ) 829 830 831 CASE ( 'diss' ) 832 IF ( av == 0 ) THEN 833 DO i = nxl, nxr 834 DO j = nys, nyn 835 DO k = nzb_do, nzt_do 836 local_pf(i,j,k) = diss(k,j,i) 837 ENDDO 838 ENDDO 839 ENDDO 840 ELSE 841 IF ( .NOT. ALLOCATED( diss_av ) ) THEN 842 ALLOCATE( diss_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 843 diss_av = REAL( fill_value, KIND = wp ) 844 ENDIF 845 DO i = nxl, nxr 846 DO j = nys, nyn 847 DO k = nzb_do, nzt_do 848 local_pf(i,j,k) = diss_av(k,j,i) 849 ENDDO 850 ENDDO 851 ENDDO 832 to_be_resorted => diss_av 852 833 ENDIF 853 834 854 835 CASE ( 'kh' ) 855 836 IF ( av == 0 ) THEN 856 DO i = nxl, nxr 857 DO j = nys, nyn 858 DO k = nzb_do, nzt_do 859 local_pf(i,j,k) = kh(k,j,i) 860 ENDDO 861 ENDDO 862 ENDDO 837 to_be_resorted => kh 863 838 ELSE 864 839 IF ( .NOT. ALLOCATED( kh_av ) ) THEN … … 866 841 kh_av = REAL( fill_value, KIND = wp ) 867 842 ENDIF 868 DO i = nxl, nxr 869 DO j = nys, nyn 870 DO k = nzb_do, nzt_do 871 local_pf(i,j,k) = kh_av(k,j,i) 872 ENDDO 873 ENDDO 874 ENDDO 843 to_be_resorted => kh_av 875 844 ENDIF 876 845 877 846 CASE ( 'km' ) 878 847 IF ( av == 0 ) THEN 879 DO i = nxl, nxr 880 DO j = nys, nyn 881 DO k = nzb_do, nzt_do 882 local_pf(i,j,k) = km(k,j,i) 883 ENDDO 884 ENDDO 885 ENDDO 848 to_be_resorted => km 886 849 ELSE 887 850 IF ( .NOT. ALLOCATED( km_av ) ) THEN … … 889 852 km_av = REAL( fill_value, KIND = wp ) 890 853 ENDIF 891 DO i = nxl, nxr 892 DO j = nys, nyn 893 DO k = nzb_do, nzt_do 894 local_pf(i,j,k) = km_av(k,j,i) 895 ENDDO 896 ENDDO 897 ENDDO 854 to_be_resorted => km_av 898 855 ENDIF 899 856 900 857 CASE ( 'dummy3' ) !> @todo remove later 901 858 IF ( av == 0 ) THEN 902 DO i = nxl, nxr 903 DO j = nys, nyn 904 DO k = nzb_do, nzt_do 905 local_pf(i,j,k) = dummy3(k,j,i) 906 ENDDO 907 ENDDO 908 ENDDO 859 to_be_resorted => dummy3 909 860 ENDIF 910 861 911 862 CASE ( 'diss1' ) !> @todo remove later 912 863 IF ( av == 0 ) THEN 913 DO i = nxl, nxr 914 DO j = nys, nyn 915 DO k = nzb_do, nzt_do 916 local_pf(i,j,k) = dummy1(k,j,i) 917 ENDDO 918 ENDDO 919 ENDDO 864 to_be_resorted => dummy1 920 865 ENDIF 921 866 922 867 CASE ( 'diss2' ) !> @todo remove later 923 868 IF ( av == 0 ) THEN 924 DO i = nxl, nxr 925 DO j = nys, nyn 926 DO k = nzb_do, nzt_do 927 local_pf(i,j,k) = dummy2(k,j,i) 928 ENDDO 929 ENDDO 930 ENDDO 869 to_be_resorted => dummy2 931 870 ENDIF 932 871 933 872 CASE ( 'diss_prod1' ) !> @todo remove later 934 873 IF ( av == 0 ) THEN 935 DO i = nxl, nxr 936 DO j = nys, nyn 937 DO k = nzb_do, nzt_do 938 local_pf(i,j,k) = diss_prod1(k,j,i) 939 ENDDO 940 ENDDO 941 ENDDO 874 to_be_resorted => diss_prod1 942 875 ENDIF 943 876 944 877 CASE ( 'diss_adve1' ) !> @todo remove later 945 878 IF ( av == 0 ) THEN 946 DO i = nxl, nxr 947 DO j = nys, nyn 948 DO k = nzb_do, nzt_do 949 local_pf(i,j,k) = diss_adve1(k,j,i) 950 ENDDO 951 ENDDO 952 ENDDO 879 to_be_resorted => diss_adve1 953 880 ENDIF 954 881 955 882 CASE ( 'diss_diff1' ) !> @todo remove later 956 883 IF ( av == 0 ) THEN 957 DO i = nxl, nxr 958 DO j = nys, nyn 959 DO k = nzb_do, nzt_do 960 local_pf(i,j,k) = diss_diff1(k,j,i) 961 ENDDO 962 ENDDO 963 ENDDO 884 to_be_resorted => diss_diff1 964 885 ENDIF 965 886 966 887 CASE ( 'diss_prod2' ) !> @todo remove later 967 888 IF ( av == 0 ) THEN 968 DO i = nxl, nxr 969 DO j = nys, nyn 970 DO k = nzb_do, nzt_do 971 local_pf(i,j,k) = diss_prod2(k,j,i) 972 ENDDO 973 ENDDO 974 ENDDO 889 to_be_resorted => diss_prod2 975 890 ENDIF 976 891 977 892 CASE ( 'diss_adve2' ) !> @todo remove later 978 893 IF ( av == 0 ) THEN 979 DO i = nxl, nxr 980 DO j = nys, nyn 981 DO k = nzb_do, nzt_do 982 local_pf(i,j,k) = diss_adve2(k,j,i) 983 ENDDO 984 ENDDO 985 ENDDO 894 to_be_resorted => diss_adve2 986 895 ENDIF 987 896 988 897 CASE ( 'diss_diff2' ) !> @todo remove later 989 898 IF ( av == 0 ) THEN 990 DO i = nxl, nxr 991 DO j = nys, nyn 992 DO k = nzb_do, nzt_do 993 local_pf(i,j,k) = diss_diff2(k,j,i) 994 ENDDO 995 ENDDO 996 ENDDO 899 to_be_resorted => diss_diff2 997 900 ENDIF 998 901 999 902 CASE ( 'diss_prod3' ) !> @todo remove later 1000 903 IF ( av == 0 ) THEN 1001 DO i = nxl, nxr 1002 DO j = nys, nyn 1003 DO k = nzb_do, nzt_do 1004 local_pf(i,j,k) = diss_prod3(k,j,i) 1005 ENDDO 1006 ENDDO 1007 ENDDO 904 to_be_resorted => diss_prod3 1008 905 ENDIF 1009 906 1010 907 CASE ( 'diss_adve3' ) !> @todo remove later 1011 908 IF ( av == 0 ) THEN 1012 DO i = nxl, nxr 1013 DO j = nys, nyn 1014 DO k = nzb_do, nzt_do 1015 local_pf(i,j,k) = diss_adve3(k,j,i) 1016 ENDDO 1017 ENDDO 1018 ENDDO 909 to_be_resorted => diss_adve3 1019 910 ENDIF 1020 911 1021 912 CASE ( 'diss_diff3' ) !> @todo remove later 1022 913 IF ( av == 0 ) THEN 1023 DO i = nxl, nxr 1024 DO j = nys, nyn 1025 DO k = nzb_do, nzt_do 1026 local_pf(i,j,k) = diss_diff3(k,j,i) 1027 ENDDO 1028 ENDDO 1029 ENDDO 914 to_be_resorted => diss_diff3 1030 915 ENDIF 1031 916 … … 1034 919 1035 920 END SELECT 921 922 923 IF ( found .AND. .NOT. resorted ) THEN 924 DO i = nxl, nxr 925 DO j = nys, nyn 926 DO k = nzb_do, nzt_do 927 local_pf(i,j,k) = MERGE( & 928 to_be_resorted(k,j,i), & 929 REAL( fill_value, KIND = wp ), & 930 BTEST( wall_flags_0(k,j,i), flag_nr ) ) 931 ENDDO 932 ENDDO 933 ENDDO 934 resorted = .TRUE. 935 ENDIF 1036 936 1037 937 END SUBROUTINE tcm_data_output_3d … … 1168 1068 1169 1069 IF ( INDEX( initializing_actions, 'set_1d-model_profiles' ) /= 0 ) THEN 1170 ! 1171 !-- Transfer initial profiles to the arrays of the 3D model 1172 DO i = nxlg, nxrg 1173 DO j = nysg, nyng 1174 e(:,j,i) = e1d 1175 kh(:,j,i) = kh1d 1176 km(:,j,i) = km1d 1070 1071 IF ( .NOT. rans_tke_e ) THEN 1072 ! 1073 !-- Transfer initial profiles to the arrays of the 3D model 1074 DO i = nxlg, nxrg 1075 DO j = nysg, nyng 1076 e(:,j,i) = e1d 1077 kh(:,j,i) = kh1d 1078 km(:,j,i) = km1d 1079 ENDDO 1177 1080 ENDDO 1178 ENDDO 1179 1180 IF ( constant_diffusion ) THEN1181 e = 0.0_wp1182 ENDIF 1183 1184 IF ( rans_tke_e ) THEN 1185 IF ( dissipation_1d == 'prognostic' ) THEN !> @query Why must this be checked? 1186 DO i = nxlg, nxrg !> Should 'diss' not always1187 DO j = nysg, nyng !> be prognostic in case rans_tke_e?1188 diss(:,j,i) = diss1d1189 ENDDO1190 ENDDO1191 ELSE 1081 1082 IF ( constant_diffusion ) THEN 1083 e = 0.0_wp 1084 ENDIF 1085 1086 ELSE 1087 ! 1088 !-- In case of TKE-e closure in RANS mode, do not use e, diss, and km 1089 !-- profiles from 1D model. Instead, initialize with constant profiles 1090 IF ( constant_diffusion ) THEN 1091 km = km_constant 1092 kh = km / prandtl_number 1093 e = 0.0_wp 1094 ELSEIF ( e_init > 0.0_wp ) THEN 1192 1095 DO i = nxlg, nxrg 1193 1096 DO j = nysg, nyng 1194 1097 DO k = nzb+1, nzt 1195 diss(k,j,i) = c_0**4 * e(k,j,i)**2 / km1d(k)1098 km(k,j,i) = c_0 * l_wall(k,j,i) * SQRT( e_init ) 1196 1099 ENDDO 1197 1100 ENDDO 1198 1101 ENDDO 1102 km(nzb,:,:) = km(nzb+1,:,:) 1103 km(nzt+1,:,:) = km(nzt,:,:) 1104 kh = km / prandtl_number 1105 e = e_init 1106 ELSE 1107 IF ( .NOT. ocean ) THEN 1108 kh = 0.01_wp ! there must exist an initial diffusion, because 1109 km = 0.01_wp ! otherwise no TKE would be produced by the 1110 ! production terms, as long as not yet 1111 ! e = (u*/cm)**2 at k=nzb+1 1112 ELSE 1113 kh = 0.00001_wp 1114 km = 0.00001_wp 1115 ENDIF 1116 e = 0.0_wp 1199 1117 ENDIF 1118 1119 DO i = nxlg, nxrg 1120 DO j = nysg, nyng 1121 DO k = nzb+1, nzt 1122 diss(k,j,i) = c_0**4 * e(k,j,i)**2 / km(k,j,i) 1123 ENDDO 1124 ENDDO 1125 ENDDO 1126 diss(nzb,:,:) = diss(nzb+1,:,:) 1127 diss(nzt+1,:,:) = diss(nzt,:,:) 1128 1200 1129 ENDIF 1201 1130 … … 1381 1310 1382 1311 USE control_parameters, & 1383 ONLY: bc_lr_cyc, bc_ns_cyc, f, kappa, message_string, & 1384 wall_adjustment_factor 1312 ONLY: bc_lr_cyc, bc_ns_cyc, f, message_string, wall_adjustment_factor 1385 1313 1386 1314 USE grid_variables, & … … 1973 1901 1974 1902 USE surface_mod, & 1975 ONLY : surf_def_h, surf_ def_v, surf_lsm_h, surf_usm_h1903 ONLY : surf_def_h, surf_lsm_h, surf_usm_h 1976 1904 1977 1905 IMPLICIT NONE … … 2556 2484 2557 2485 ! 2558 !-- Use special boundary condition in case of TKE-e closure2559 IF ( rans_tke_e ) THEN2560 surf_s = surf_def_h(0)%start_index(j,i)2561 surf_e = surf_def_h(0)%end_index(j,i)2562 DO m = surf_s, surf_e2563 k = surf_def_h(0)%k(m)2564 e_p(k,j,i) = ( surf_def_h(0)%us(m) / c_0 )**22565 ENDDO2566 2567 DO l = 0, 32568 surf_s = surf_def_v(l)%start_index(j,i)2569 surf_e = surf_def_v(l)%end_index(j,i)2570 DO m = surf_s, surf_e2571 k = surf_def_v(l)%k(m)2572 e_p(k,j,i) = ( surf_def_v(l)%us(m) / c_0 )**22573 ENDDO2574 ENDDO2575 ENDIF2576 2577 !2578 2486 !-- Calculate tendencies for the next Runge-Kutta step 2579 2487 IF ( timestep_scheme(1:5) == 'runge' ) THEN … … 2675 2583 2676 2584 ! 2677 !-- Use special boundary condition in case of TKE-e closure2678 surf_s = surf_def_h(0)%start_index(j,i)2679 surf_e = surf_def_h(0)%end_index(j,i)2680 DO m = surf_s, surf_e2681 k = surf_def_h(0)%k(m)2682 diss_p(k,j,i) = surf_def_h(0)%us(m)**3 / kappa * ddzu(k)2683 ENDDO2684 2685 DO l = 0, 12686 surf_s = surf_def_v(l)%start_index(j,i)2687 surf_e = surf_def_v(l)%end_index(j,i)2688 DO m = surf_s, surf_e2689 k = surf_def_v(l)%k(m)2690 diss_p(k,j,i) = surf_def_v(l)%us(m)**3 / ( kappa * 0.5_wp * dy )2691 ENDDO2692 ENDDO2693 2694 DO l = 2, 32695 surf_s = surf_def_v(l)%start_index(j,i)2696 surf_e = surf_def_v(l)%end_index(j,i)2697 DO m = surf_s, surf_e2698 k = surf_def_v(l)%k(m)2699 diss_p(k,j,i) = surf_def_v(l)%us(m)**3 / ( kappa * 0.5_wp * dx )2700 ENDDO2701 ENDDO2702 !2703 2585 !-- Calculate tendencies for the next Runge-Kutta step 2704 2586 IF ( timestep_scheme(1:5) == 'runge' ) THEN … … 2715 2597 ENDIF 2716 2598 ENDIF 2717 !2718 !-- Limit change of diss to be between -90% and +100%. Also, set an absolute2719 !-- minimum value2720 DO k = nzb, nzt+12721 diss_p(k,j,i) = MIN( MAX( diss_p(k,j,i), &2722 0.1_wp * diss(k,j,i), &2723 0.0001_wp ), &2724 2.0_wp * diss(k,j,i) )2725 ENDDO2726 2599 2727 2600 IF ( intermediate_timestep_count == 1 ) dummy1(:,j,i) = diss_p(:,j,i) !> @todo remove later … … 4777 4650 !> Computation of the turbulent diffusion coefficients for momentum and heat 4778 4651 !> according to Prandtl-Kolmogorov. 4779 !> @todo consider non-default surfaces4780 4652 !------------------------------------------------------------------------------! 4781 4653 SUBROUTINE tcm_diffusivities_default( var, var_reference ) … … 4792 4664 4793 4665 USE surface_mod, & 4794 ONLY : bc_h, surf_def_h, surf_def_v 4666 ONLY : bc_h, surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, & 4667 surf_usm_h, surf_usm_v 4795 4668 4796 4669 IMPLICIT NONE … … 4956 4829 ! 4957 4830 !-- Up- and downward facing surfaces 4831 !-- Default surfaces 4958 4832 DO n = 0, 1 4959 4833 DO m = 1, surf_def_h(n)%ns … … 4961 4835 j = surf_def_h(n)%j(m) 4962 4836 k = surf_def_h(n)%k(m) 4963 km(k,j,i) = kappa * surf_def_h(n)%us(m) * dzu(k)4837 km(k,j,i) = kappa * surf_def_h(n)%us(m) * surf_def_h(n)%z_mo(m) 4964 4838 kh(k,j,i) = 1.35_wp * km(k,j,i) 4965 4839 ENDDO 4966 4840 ENDDO 4967 4841 ! 4968 !-- North- and southward facing surfaces 4969 DO n = 0, 1 4842 !-- Upward facing surfaces 4843 !-- Natural surfaces 4844 DO m = 1, surf_lsm_h%ns 4845 i = surf_lsm_h%i(m) 4846 j = surf_lsm_h%j(m) 4847 k = surf_lsm_h%k(m) 4848 km(k,j,i) = kappa * surf_lsm_h%us(m) * surf_lsm_h%z_mo(m) 4849 kh(k,j,i) = 1.35_wp * km(k,j,i) 4850 ENDDO 4851 ! 4852 !-- Urban surfaces 4853 DO m = 1, surf_usm_h%ns 4854 i = surf_usm_h%i(m) 4855 j = surf_usm_h%j(m) 4856 k = surf_usm_h%k(m) 4857 km(k,j,i) = kappa * surf_usm_h%us(m) * surf_usm_h%z_mo(m) 4858 kh(k,j,i) = 1.35_wp * km(k,j,i) 4859 ENDDO 4860 4861 ! 4862 !-- North-, south-, west and eastward facing surfaces 4863 DO n = 0, 3 4864 ! 4865 !-- Default surfaces 4970 4866 DO m = 1, surf_def_v(n)%ns 4971 4867 i = surf_def_v(n)%i(m) 4972 4868 j = surf_def_v(n)%j(m) 4973 4869 k = surf_def_v(n)%k(m) 4974 km(k,j,i) = kappa * surf_def_v(n)%us(m) * 0.5_wp * dy4870 km(k,j,i) = kappa * surf_def_v(n)%us(m) * surf_def_v(n)%z_mo(m) 4975 4871 kh(k,j,i) = 1.35_wp * km(k,j,i) 4976 4872 ENDDO 4977 ENDDO 4978 ! 4979 !-- West- and eastward facing surfaces 4980 DO n = 2, 3 4981 DO m = 1, surf_def_v(n)%ns 4982 i = surf_def_v(n)%i(m) 4983 j = surf_def_v(n)%j(m) 4984 k = surf_def_v(n)%k(m) 4985 km(k,j,i) = kappa * surf_def_v(n)%us(m) * 0.5_wp * dx 4873 ! 4874 !-- Natural surfaces 4875 DO m = 1, surf_lsm_v(n)%ns 4876 i = surf_lsm_v(n)%i(m) 4877 j = surf_lsm_v(n)%j(m) 4878 k = surf_lsm_v(n)%k(m) 4879 km(k,j,i) = kappa * surf_lsm_v(n)%us(m) * surf_lsm_v(n)%z_mo(m) 4880 kh(k,j,i) = 1.35_wp * km(k,j,i) 4881 ENDDO 4882 ! 4883 !-- Urban surfaces 4884 DO m = 1, surf_usm_v(n)%ns 4885 i = surf_usm_v(n)%i(m) 4886 j = surf_usm_v(n)%j(m) 4887 k = surf_usm_v(n)%k(m) 4888 km(k,j,i) = kappa * surf_usm_v(n)%us(m) * surf_usm_v(n)%z_mo(m) 4986 4889 kh(k,j,i) = 1.35_wp * km(k,j,i) 4987 4890 ENDDO
Note: See TracChangeset
for help on using the changeset viewer.
