Changeset 1557 for palm/trunk/SOURCE
- Timestamp:
- Mar 5, 2015 4:43:04 PM (10 years ago)
- Location:
- palm/trunk/SOURCE
- Files:
-
- 7 edited
-
advec_ws.f90 (modified) (15 diffs)
-
check_parameters.f90 (modified) (7 diffs)
-
flow_statistics.f90 (modified) (13 diffs)
-
header.f90 (modified) (3 diffs)
-
init_grid.f90 (modified) (2 diffs)
-
init_pegrid.f90 (modified) (2 diffs)
-
modules.f90 (modified) (2 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/advec_ws.f90
r1375 r1557 20 20 ! Current revisions: 21 21 ! ------------------ 22 ! 22 ! Enable monotone advection for scalars using monotonic limiter 23 23 ! 24 24 ! Former revisions: … … 194 194 USE control_parameters, & 195 195 ONLY: cloud_physics, humidity, icloud_scheme, loop_optimization, & 196 passive_scalar, precipitation, ocean, ws_scheme_mom,&197 ws_scheme_ sca196 monotonic_adjustment, passive_scalar, precipitation, ocean, & 197 ws_scheme_mom, ws_scheme_sca 198 198 199 199 USE indices, & … … 386 386 387 387 USE control_parameters, & 388 ONLY: intermediate_timestep_count, u_gtrans, v_gtrans 388 ONLY: intermediate_timestep_count, monotonic_adjustment, u_gtrans, & 389 v_gtrans 389 390 390 391 USE grid_variables, & … … 421 422 INTEGER(iwp) :: k !: 422 423 INTEGER(iwp) :: k_mm !: 424 INTEGER(iwp) :: k_mmm !: 423 425 INTEGER(iwp) :: k_pp !: 424 426 INTEGER(iwp) :: k_ppp !: … … 428 430 REAL(wp) :: div !: 429 431 REAL(wp) :: flux_d !: 432 REAL(wp) :: fd_1 !: 433 REAL(wp) :: fl_1 !: 434 REAL(wp) :: fn_1 !: 435 REAL(wp) :: fr_1 !: 436 REAL(wp) :: fs_1 !: 437 REAL(wp) :: ft_1 !: 438 REAL(wp) :: phi_d !: 439 REAL(wp) :: phi_l !: 440 REAL(wp) :: phi_n !: 441 REAL(wp) :: phi_r !: 442 REAL(wp) :: phi_s !: 443 REAL(wp) :: phi_t !: 444 REAL(wp) :: rd !: 445 REAL(wp) :: rl !: 446 REAL(wp) :: rn !: 447 REAL(wp) :: rr !: 448 REAL(wp) :: rs !: 449 REAL(wp) :: rt !: 430 450 REAL(wp) :: u_comp !: 431 451 REAL(wp) :: v_comp !: … … 697 717 ) 698 718 ! 719 !-- Apply monotonic adjustment. 720 IF ( monotonic_adjustment ) THEN 721 ! 722 !-- At first, calculate first order fluxes. 723 u_comp = u(k,j,i) - u_gtrans 724 fl_1 = ( u_comp * ( sk(k,j,i) + sk(k,j,i-1) ) & 725 -ABS( u_comp ) * ( sk(k,j,i) - sk(k,j,i-1) ) & 726 ) * adv_sca_1 727 728 u_comp = u(k,j,i+1) - u_gtrans 729 fr_1 = ( u_comp * ( sk(k,j,i+1) + sk(k,j,i) ) & 730 -ABS( u_comp ) * ( sk(k,j,i+1) - sk(k,j,i) ) & 731 ) * adv_sca_1 732 733 v_comp = v(k,j,i) - v_gtrans 734 fs_1 = ( v_comp * ( sk(k,j,i) + sk(k,j-1,i) ) & 735 -ABS( v_comp ) * ( sk(k,j,i) - sk(k,j-1,i) ) & 736 ) * adv_sca_1 737 738 v_comp = v(k,j+1,i) - v_gtrans 739 fn_1 = ( v_comp * ( sk(k,j+1,i) + sk(k,j,i) ) & 740 -ABS( v_comp ) * ( sk(k,j+1,i) - sk(k,j,i) ) & 741 ) * adv_sca_1 742 743 fd_1 = ( w(k-1,j,i) * ( sk(k,j,i) + sk(k-1,j,i) ) & 744 -ABS( w(k-1,j,i) ) * ( sk(k,j,i) - sk(k-1,j,i) ) & 745 ) * adv_sca_1 746 747 ft_1 = ( w(k,j,i) * ( sk(k+1,j,i) + sk(k,j,i) ) & 748 -ABS( w(k,j,i) ) * ( sk(k+1,j,i) - sk(k,j,i) ) & 749 ) * adv_sca_1 750 ! 751 !-- Calculate ratio of upwind gradients. Note, Min/Max is just to 752 !-- avoid if statements. 753 rl = ( MAX( 0.0, u(k,j,i) - u_gtrans ) * & 754 ABS( ( sk(k,j,i-1) - sk(k,j,i-2) ) / & 755 ( sk(k,j,i) - sk(k,j,i-1) + 1E-20_wp ) & 756 ) + & 757 MIN( 0.0, u(k,j,i) - u_gtrans ) * & 758 ABS( ( sk(k,j,i) - sk(k,j,i+1) ) / & 759 ( sk(k,j,i-1) - sk(k,j,i) + 1E-20_wp ) & 760 ) & 761 ) / ABS( u(k,j,i) - u_gtrans + 1E-20_wp ) 762 763 rr = ( MAX( 0.0, u(k,j,i+1) - u_gtrans ) * & 764 ABS( ( sk(k,j,i) - sk(k,j,i-1) ) / & 765 ( sk(k,j,i+1) - sk(k,j,i) + 1E-20_wp ) & 766 ) + & 767 MIN( 0.0, u(k,j,i+1) - u_gtrans ) * & 768 ABS( ( sk(k,j,i+1) - sk(k,j,i+2) ) / & 769 ( sk(k,j,i) - sk(k,j,i+1) + 1E-20_wp ) & 770 ) & 771 ) / ABS( u(k,j,i+1) - u_gtrans + 1E-20_wp ) 772 773 rs = ( MAX( 0.0, v(k,j,i) - v_gtrans ) * & 774 ABS( ( sk(k,j-1,i) - sk(k,j-2,i) ) / & 775 ( sk(k,j,i) - sk(k,j-1,i) + 1E-20_wp ) & 776 ) + & 777 MIN( 0.0, v(k,j,i) - v_gtrans ) * & 778 ABS( ( sk(k,j,i) - sk(k,j+1,i) ) / & 779 ( sk(k,j-1,i) - sk(k,j,i) + 1E-20_wp ) & 780 ) & 781 ) / ABS( v(k,j,i) - v_gtrans + 1E-20_wp ) 782 783 rn = ( MAX( 0.0, v(k,j+1,i) - v_gtrans ) * & 784 ABS( ( sk(k,j,i) - sk(k,j-1,i) ) / & 785 ( sk(k,j+1,i) - sk(k,j,i) + 1E-20_wp ) & 786 ) + & 787 MIN( 0.0, v(k,j+1,i) - v_gtrans ) * & 788 ABS( ( sk(k,j+1,i) - sk(k,j+2,i) ) / & 789 ( sk(k,j,i) - sk(k,j+1,i) + 1E-20_wp ) & 790 ) & 791 ) / ABS( v(k,j+1,i) - v_gtrans + 1E-20_wp ) 792 ! 793 !-- Reuse k_mm and compute k_mmm for the vertical gradient ratios. 794 !-- Note, for vertical advection below the third grid point above 795 !-- surface ( or below the model top) rd and rt are set to 0, i.e. 796 !-- use of first order scheme is enforced. 797 k_mmm = k - 3 * ibit8 798 799 rd = ( MAX( 0.0, w(k-1,j,i) ) * & 800 ABS( ( sk(k_mm,j,i) - sk(k_mmm,j,i) ) / & 801 ( sk(k-1,j,i) - sk(k_mm,j,i) + 1E-20_wp ) & 802 ) + & 803 MIN( 0.0, w(k-1,j,i) ) * & 804 ABS( ( sk(k-1,j,i) - sk(k,j,i) ) / & 805 ( sk(k_mm,j,i) - sk(k-1,j,i) + 1E-20_wp ) & 806 ) & 807 ) * ibit8 / ABS( w(k-1,j,i) + 1E-20_wp ) 808 809 rt = ( MAX( 0.0, w(k,j,i) ) * & 810 ABS( ( sk(k,j,i) - sk(k-1,j,i) ) / & 811 ( sk(k+1,j,i) - sk(k,j,i) + 1E-20_wp ) & 812 ) + & 813 MIN( 0.0, w(k,j,i) ) * & 814 ABS( ( sk(k+1,j,i) - sk(k_pp,j,i) ) / & 815 ( sk(k,j,i) - sk(k+1,j,i) + 1E-20_wp ) & 816 ) & 817 ) * ibit8 / ABS( w(k,j,i) + 1E-20_wp ) 818 ! 819 !-- Calculate empirical limiter function (van Albada2 limiter). 820 phi_l = MIN( 1.0_wp, ( 2.0_wp * rl ) / & 821 ( rl**2 + 1.0_wp ) ) 822 phi_r = MIN( 1.0_wp, ( 2.0_wp * rr ) / & 823 ( rr**2 + 1.0_wp ) ) 824 phi_s = MIN( 1.0_wp, ( 2.0_wp * rs ) / & 825 ( rs**2 + 1.0_wp ) ) 826 phi_n = MIN( 1.0_wp, ( 2.0_wp * rn ) / & 827 ( rn**2 + 1.0_wp ) ) 828 phi_d = MIN( 1.0_wp, ( 2.0_wp * rd ) / & 829 ( rd**2 + 1.0_wp ) ) 830 phi_t = MIN( 1.0_wp, ( 2.0_wp * rt ) / & 831 ( rt**2 + 1.0_wp ) ) 832 ! 833 !-- Calculate the resulting monotone flux. 834 swap_flux_x_local(k,j,tn) = fl_1 - phi_l * & 835 ( fl_1 - swap_flux_x_local(k,j,tn) ) 836 flux_r(k) = fr_1 - phi_r * & 837 ( fr_1 - flux_r(k) ) 838 swap_flux_y_local(k,tn) = fs_1 - phi_s * & 839 ( fs_1 - swap_flux_y_local(k,tn) ) 840 flux_n(k) = fn_1 - phi_n * & 841 ( fn_1 - flux_n(k) ) 842 flux_d = fd_1 !- phi_d * & 843 ! ( fd_1 - flux_d ) 844 flux_t(k) = ft_1 !- phi_t * & 845 ! ( ft_1 - flux_t(k) ) 846 ! 847 !-- Moreover, modify dissipation flux according to the limiter. 848 swap_diss_x_local(k,j,tn) = swap_diss_x_local(k,j,tn) * phi_l 849 diss_r(k) = diss_r(k) * phi_r 850 swap_diss_y_local(k,tn) = swap_diss_y_local(k,tn) * phi_s 851 diss_n(k) = diss_n(k) * phi_n 852 diss_d = 0.0 !diss_d * phi_d 853 diss_t(k) = 0.0 !diss_t(k) * phi_t 854 855 ENDIF 856 ! 699 857 !-- Calculate the divergence of the velocity field. A respective 700 858 !-- correction is needed to overcome numerical instabilities caused … … 756 914 k_pp = k + 2 * ( 1 - ibit6 ) 757 915 k_mm = k - 2 * ibit8 758 759 916 760 917 flux_t(k) = w(k,j,i) * ( & … … 786 943 ( sk(k_ppp,j,i) - sk(k_mm,j,i) ) & 787 944 ) 945 946 947 ! 948 !-- Apply monotonic adjustment. 949 IF ( monotonic_adjustment ) THEN 950 ! 951 !-- At first, calculate first order fluxes. 952 u_comp = u(k,j,i) - u_gtrans 953 fl_1 = ( u_comp * ( sk(k,j,i) + sk(k,j,i-1) ) & 954 -ABS( u_comp ) * ( sk(k,j,i) - sk(k,j,i-1) ) & 955 ) * adv_sca_1 956 957 u_comp = u(k,j,i+1) - u_gtrans 958 fr_1 = ( u_comp * ( sk(k,j,i+1) + sk(k,j,i) ) & 959 -ABS( u_comp ) * ( sk(k,j,i+1) - sk(k,j,i) ) & 960 ) * adv_sca_1 961 962 v_comp = v(k,j,i) - v_gtrans 963 fs_1 = ( v_comp * ( sk(k,j,i) + sk(k,j-1,i) ) & 964 -ABS( v_comp ) * ( sk(k,j,i) - sk(k,j-1,i) ) & 965 ) * adv_sca_1 966 967 v_comp = v(k,j+1,i) - v_gtrans 968 fn_1 = ( v_comp * ( sk(k,j+1,i) + sk(k,j,i) ) & 969 -ABS( v_comp ) * ( sk(k,j+1,i) - sk(k,j,i) ) & 970 ) * adv_sca_1 971 972 fd_1 = ( w(k-1,j,i) * ( sk(k,j,i) + sk(k-1,j,i) ) & 973 -ABS( w(k-1,j,i) ) * ( sk(k,j,i) - sk(k-1,j,i) ) & 974 ) * adv_sca_1 975 976 ft_1 = ( w(k,j,i) * ( sk(k+1,j,i) + sk(k,j,i) ) & 977 -ABS( w(k,j,i) ) * ( sk(k+1,j,i) - sk(k,j,i) ) & 978 ) * adv_sca_1 979 ! 980 !-- Calculate ratio of upwind gradients. Note, Min/Max is just to 981 !-- avoid if statements. 982 rl = ( MAX( 0.0, u(k,j,i) - u_gtrans ) * & 983 ABS( ( sk(k,j,i-1) - sk(k,j,i-2) ) / & 984 ( sk(k,j,i) - sk(k,j,i-1) + 1E-20_wp ) & 985 ) + & 986 MIN( 0.0, u(k,j,i) - u_gtrans ) * & 987 ABS( ( sk(k,j,i) - sk(k,j,i+1) ) / & 988 ( sk(k,j,i-1) - sk(k,j,i) + 1E-20_wp ) & 989 ) & 990 ) / ABS( u(k,j,i) - u_gtrans + 1E-20_wp ) 991 992 rr = ( MAX( 0.0, u(k,j,i+1) - u_gtrans ) * & 993 ABS( ( sk(k,j,i) - sk(k,j,i-1) ) / & 994 ( sk(k,j,i+1) - sk(k,j,i) + 1E-20_wp ) & 995 ) + & 996 MIN( 0.0, u(k,j,i+1) - u_gtrans ) * & 997 ABS( ( sk(k,j,i+1) - sk(k,j,i+2) ) / & 998 ( sk(k,j,i) - sk(k,j,i+1) + 1E-20_wp ) & 999 ) & 1000 ) / ABS( u(k,j,i+1) - u_gtrans + 1E-20_wp ) 1001 1002 rs = ( MAX( 0.0, v(k,j,i) - v_gtrans ) * & 1003 ABS( ( sk(k,j-1,i) - sk(k,j-2,i) ) / & 1004 ( sk(k,j,i) - sk(k,j-1,i) + 1E-20_wp ) & 1005 ) + & 1006 MIN( 0.0, v(k,j,i) - v_gtrans ) * & 1007 ABS( ( sk(k,j,i) - sk(k,j+1,i) ) / & 1008 ( sk(k,j-1,i) - sk(k,j,i) + 1E-20_wp ) & 1009 ) & 1010 ) / ABS( v(k,j,i) - v_gtrans + 1E-20_wp ) 1011 1012 rn = ( MAX( 0.0, v(k,j+1,i) - v_gtrans ) * & 1013 ABS( ( sk(k,j,i) - sk(k,j-1,i) ) / & 1014 ( sk(k,j+1,i) - sk(k,j,i) + 1E-20_wp ) & 1015 ) + & 1016 MIN( 0.0, v(k,j+1,i) - v_gtrans ) * & 1017 ABS( ( sk(k,j+1,i) - sk(k,j+2,i) ) / & 1018 ( sk(k,j,i) - sk(k,j+1,i) + 1E-20_wp ) & 1019 ) & 1020 ) / ABS( v(k,j+1,i) - v_gtrans + 1E-20_wp ) 1021 ! 1022 !-- Reuse k_mm and compute k_mmm for the vertical gradient ratios. 1023 !-- Note, for vertical advection below the third grid point above 1024 !-- surface ( or below the model top) rd and rt are set to 0, i.e. 1025 !-- use of first order scheme is enforced. 1026 k_mmm = k - 3 * ibit8 1027 1028 rd = ( MAX( 0.0, w(k-1,j,i) ) * & 1029 ABS( ( sk(k_mm,j,i) - sk(k_mmm,j,i) ) / & 1030 ( sk(k-1,j,i) - sk(k_mm,j,i) + 1E-20_wp ) & 1031 ) + & 1032 MIN( 0.0, w(k-1,j,i) ) * & 1033 ABS( ( sk(k-1,j,i) - sk(k,j,i) ) / & 1034 ( sk(k_mm,j,i) - sk(k-1,j,i) + 1E-20_wp ) & 1035 ) & 1036 ) * ibit8 / ABS( w(k-1,j,i) + 1E-20_wp ) 1037 1038 rt = ( MAX( 0.0, w(k,j,i) ) * & 1039 ABS( ( sk(k,j,i) - sk(k-1,j,i) ) / & 1040 ( sk(k+1,j,i) - sk(k,j,i) + 1E-20_wp ) & 1041 ) + & 1042 MIN( 0.0, w(k,j,i) ) * & 1043 ABS( ( sk(k+1,j,i) - sk(k_pp,j,i) ) / & 1044 ( sk(k,j,i) - sk(k+1,j,i) + 1E-20_wp ) & 1045 ) & 1046 ) * ibit8 / ABS( w(k,j,i) + 1E-20_wp ) 1047 ! 1048 !-- Calculate empirical limiter function (van Albada2 limiter). 1049 phi_l = MIN( 1.0_wp, ( 2.0_wp * rl ) / & 1050 ( rl**2 + 1.0_wp ) ) 1051 phi_r = MIN( 1.0_wp, ( 2.0_wp * rr ) / & 1052 ( rr**2 + 1.0_wp ) ) 1053 phi_s = MIN( 1.0_wp, ( 2.0_wp * rs ) / & 1054 ( rs**2 + 1.0_wp ) ) 1055 phi_n = MIN( 1.0_wp, ( 2.0_wp * rn ) / & 1056 ( rn**2 + 1.0_wp ) ) 1057 phi_d = MIN( 1.0_wp, ( 2.0_wp * rd ) / & 1058 ( rd**2 + 1.0_wp ) ) 1059 phi_t = MIN( 1.0_wp, ( 2.0_wp * rt ) / & 1060 ( rt**2 + 1.0_wp ) ) 1061 ! 1062 !-- Calculate the resulting monotone flux. 1063 swap_flux_x_local(k,j,tn) = fl_1 - phi_l * & 1064 ( fl_1 - swap_flux_x_local(k,j,tn) ) 1065 flux_r(k) = fr_1 - phi_r * & 1066 ( fr_1 - flux_r(k) ) 1067 swap_flux_y_local(k,tn) = fs_1 - phi_s * & 1068 ( fs_1 - swap_flux_y_local(k,tn) ) 1069 flux_n(k) = fn_1 - phi_n * & 1070 ( fn_1 - flux_n(k) ) 1071 flux_d = fd_1 - phi_d * & 1072 ( fd_1 - flux_d ) 1073 flux_t(k) = ft_1 - phi_t * & 1074 ( ft_1 - flux_t(k) ) 1075 ! 1076 !-- Moreover, modify dissipation flux according to the limiter. 1077 swap_diss_x_local(k,j,tn) = swap_diss_x_local(k,j,tn) * phi_l 1078 diss_r(k) = diss_r(k) * phi_r 1079 swap_diss_y_local(k,tn) = swap_diss_y_local(k,tn) * phi_s 1080 diss_n(k) = diss_n(k) * phi_n 1081 diss_d = diss_d * phi_d 1082 diss_t(k) = diss_t(k) * phi_t 1083 1084 ENDIF 788 1085 ! 789 1086 !-- Calculate the divergence of the velocity field. A respective … … 803 1100 ) + sk(k,j,i) * div 804 1101 1102 805 1103 swap_flux_y_local(k,tn) = flux_n(k) 806 1104 swap_diss_y_local(k,tn) = diss_n(k) … … 813 1111 814 1112 ! 815 !-- Evaluation of statistics 816 SELECT CASE ( sk_char ) 817 818 CASE ( 'pt' ) 819 820 DO k = nzb, nzt 821 sums_wspts_ws_l(k,tn) = sums_wspts_ws_l(k,tn) + & 822 ( flux_t(k) + diss_t(k) ) & 823 * weight_substep(intermediate_timestep_count) 1113 !-- Evaluation of statistics. Please note, in case of using monotone limiter 1114 !-- vertical fluxes will be not calculated by the advective fluxes. 1115 IF ( .NOT. monotonic_adjustment ) THEN 1116 1117 SELECT CASE ( sk_char ) 1118 1119 CASE ( 'pt' ) 1120 1121 DO k = nzb, nzt 1122 sums_wspts_ws_l(k,tn) = sums_wspts_ws_l(k,tn) + & 1123 ( flux_t(k) + diss_t(k) ) & 1124 * weight_substep(intermediate_timestep_count) 824 1125 ENDDO 825 1126 826 CASE ( 'sa' )827 828 DO k = nzb, nzt829 sums_wssas_ws_l(k,tn) = sums_wssas_ws_l(k,tn) +&830 ( flux_t(k) + diss_t(k) )&831 * weight_substep(intermediate_timestep_count)832 ENDDO1127 CASE ( 'sa' ) 1128 1129 DO k = nzb, nzt 1130 sums_wssas_ws_l(k,tn) = sums_wssas_ws_l(k,tn) + & 1131 ( flux_t(k) + diss_t(k) ) & 1132 * weight_substep(intermediate_timestep_count) 1133 ENDDO 833 1134 834 CASE ( 'q' ) 835 836 DO k = nzb, nzt 837 sums_wsqs_ws_l(k,tn) = sums_wsqs_ws_l(k,tn) + & 838 ( flux_t(k) + diss_t(k) ) & 839 * weight_substep(intermediate_timestep_count) 840 ENDDO 841 842 CASE ( 'qr' ) 843 844 DO k = nzb, nzt 845 sums_wsqrs_ws_l(k,tn) = sums_wsqrs_ws_l(k,tn) + & 846 ( flux_t(k) + diss_t(k) ) & 847 * weight_substep(intermediate_timestep_count) 848 ENDDO 849 850 CASE ( 'nr' ) 851 852 DO k = nzb, nzt 853 sums_wsnrs_ws_l(k,tn) = sums_wsnrs_ws_l(k,tn) + & 854 ( flux_t(k) + diss_t(k) ) & 855 * weight_substep(intermediate_timestep_count) 856 ENDDO 857 858 END SELECT 1135 CASE ( 'q' ) 1136 1137 DO k = nzb, nzt 1138 sums_wsqs_ws_l(k,tn) = sums_wsqs_ws_l(k,tn) + & 1139 ( flux_t(k) + diss_t(k) ) & 1140 * weight_substep(intermediate_timestep_count) 1141 ENDDO 1142 1143 CASE ( 'qr' ) 1144 1145 DO k = nzb, nzt 1146 sums_wsqrs_ws_l(k,tn) = sums_wsqrs_ws_l(k,tn) + & 1147 ( flux_t(k) + diss_t(k) ) & 1148 * weight_substep(intermediate_timestep_count) 1149 ENDDO 1150 1151 CASE ( 'nr' ) 1152 1153 DO k = nzb, nzt 1154 sums_wsnrs_ws_l(k,tn) = sums_wsnrs_ws_l(k,tn) + & 1155 ( flux_t(k) + diss_t(k) ) & 1156 * weight_substep(intermediate_timestep_count) 1157 ENDDO 1158 1159 END SELECT 1160 ENDIF 859 1161 860 1162 END SUBROUTINE advec_s_ws_ij … … 2228 2530 2229 2531 USE control_parameters, & 2230 ONLY: intermediate_timestep_count, u_gtrans, v_gtrans 2532 ONLY: intermediate_timestep_count, monotonic_adjustment, u_gtrans,& 2533 v_gtrans 2231 2534 2232 2535 USE grid_variables, & … … 2247 2550 CHARACTER (LEN = *), INTENT(IN) :: sk_char !: 2248 2551 2249 INTEGER(iwp) :: i !:2250 INTEGER(iwp) :: ibit0 !:2251 INTEGER(iwp) :: ibit1 !:2252 INTEGER(iwp) :: ibit2 !:2253 INTEGER(iwp) :: ibit3 !:2254 INTEGER(iwp) :: ibit4 !:2255 INTEGER(iwp) :: ibit5 !:2256 INTEGER(iwp) :: ibit6 !:2257 INTEGER(iwp) :: ibit7 !:2258 INTEGER(iwp) :: ibit8 !:2259 INTEGER(iwp) :: j !:2260 INTEGER(iwp) :: k !:2261 INTEGER(iwp) :: k_mm !:2262 INTEGER(iwp) :: k_pp !:2263 INTEGER(iwp) :: k_ppp !:2264 INTEGER(iwp) :: tn = 0 !:2265 2266 #if defined( __nopointer )2267 REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: sk !:2268 #else2269 REAL(wp), DIMENSION(:,:,:), POINTER :: sk !:2270 #endif2271 2272 REAL(wp) :: diss_d !:2273 REAL(wp) :: div !:2274 REAL(wp) :: flux_d !:2275 REAL(wp) :: u_comp !:2276 REAL(wp) :: v_comp !:2277 2278 REAL(wp), DIMENSION(nzb:nzt) :: diss_n !:2279 REAL(wp), DIMENSION(nzb:nzt) :: diss_r !:2280 REAL(wp), DIMENSION(nzb:nzt) :: diss_t !:2281 REAL(wp), DIMENSION(nzb:nzt) :: flux_n !:2282 REAL(wp), DIMENSION(nzb:nzt) :: flux_r !:2283 REAL(wp), DIMENSION(nzb:nzt) :: flux_t !:2284 2285 REAL(wp), DIMENSION(nzb+1:nzt) :: swap_diss_y_local !:2286 REAL(wp), DIMENSION(nzb+1:nzt) :: swap_flux_y_local !:2287 2288 REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) :: swap_diss_x_local !:2289 REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) :: swap_flux_x_local !:2290 2291 2292 !2293 !-- Compute the fluxes for the whole left boundary of the processor domain.2294 i = nxl2295 DO j = nys, nyn2296 2297 DO k = nzb+1, nzb_max2298 2299 ibit2 = IBITS(wall_flags_0(k,j,i),2,1)2300 ibit1 = IBITS(wall_flags_0(k,j,i),1,1)2301 ibit0 = IBITS(wall_flags_0(k,j,i),0,1)2302 2303 u_comp = u(k,j,i) - u_gtrans2304 swap_flux_x_local(k,j) = u_comp * ( &2305 ( 37.0_wp * ibit2 * adv_sca_5 &2306 + 7.0_wp * ibit1 * adv_sca_3 &2307 + ibit0 * adv_sca_1 &2308 ) * &2309 ( sk(k,j,i) + sk(k,j,i-1) ) &2310 - ( 8.0_wp * ibit2 * adv_sca_5 &2311 + ibit1 * adv_sca_3 &2312 ) * &2313 ( sk(k,j,i+1) + sk(k,j,i-2) ) &2314 + ( ibit2 * adv_sca_5 &2315 ) * &2316 ( sk(k,j,i+2) + sk(k,j,i-3) ) &2317 )2318 2319 swap_diss_x_local(k,j) = -ABS( u_comp ) * ( &2320 ( 10.0_wp * ibit2 * adv_sca_5 &2321 + 3.0_wp * ibit1 * adv_sca_3 &2322 + ibit0 * adv_sca_1 &2323 ) * &2324 ( sk(k,j,i) - sk(k,j,i-1) ) &2325 - ( 5.0_wp * ibit2 * adv_sca_5 &2326 + ibit1 * adv_sca_3 &2327 ) * &2328 ( sk(k,j,i+1) - sk(k,j,i-2) ) &2329 + ( ibit2 * adv_sca_5 &2330 ) * &2331 ( sk(k,j,i+2) - sk(k,j,i-3) ) &2332 )2333 2334 ENDDO2335 2336 DO k = nzb_max+1, nzt2337 2338 u_comp = u(k,j,i) - u_gtrans2339 swap_flux_x_local(k,j) = u_comp * ( &2340 37.0_wp * ( sk(k,j,i) + sk(k,j,i-1) ) &2341 - 8.0_wp * ( sk(k,j,i+1) + sk(k,j,i-2) ) &2342 + ( sk(k,j,i+2) + sk(k,j,i-3) ) &2343 ) * adv_sca_52344 2345 swap_diss_x_local(k,j) = -ABS( u_comp ) * ( &2346 10.0_wp * ( sk(k,j,i) - sk(k,j,i-1) ) &2347 - 5.0_wp * ( sk(k,j,i+1) - sk(k,j,i-2) ) &2348 + ( sk(k,j,i+2) - sk(k,j,i-3) ) &2349 ) * adv_sca_52350 2351 ENDDO2352 2353 ENDDO2354 2355 DO i = nxl, nxr2356 2357 j = nys2358 DO k = nzb+1, nzb_max2359 2360 ibit5 = IBITS(wall_flags_0(k,j,i),5,1)2361 ibit4 = IBITS(wall_flags_0(k,j,i),4,1)2362 ibit3 = IBITS(wall_flags_0(k,j,i),3,1)2363 2364 v_comp = v(k,j,i) - v_gtrans2365 swap_flux_y_local(k) = v_comp * ( &2366 ( 37.0_wp * ibit5 * adv_sca_5 &2367 + 7.0_wp * ibit4 * adv_sca_3 &2368 + ibit3 * adv_sca_1 &2369 ) * &2370 ( sk(k,j,i) + sk(k,j-1,i) ) &2371 - ( 8.0_wp * ibit5 * adv_sca_5 &2372 + ibit4 * adv_sca_3 &2373 ) * &2374 ( sk(k,j+1,i) + sk(k,j-2,i) ) &2375 + ( ibit5 * adv_sca_5 &2376 ) * &2377 ( sk(k,j+2,i) + sk(k,j-3,i) ) &2378 )2379 2380 swap_diss_y_local(k) = -ABS( v_comp ) * ( &2381 ( 10.0_wp * ibit5 * adv_sca_5 &2382 + 3.0_wp * ibit4 * adv_sca_3 &2383 + ibit3 * adv_sca_1 &2384 ) * &2385 ( sk(k,j,i) - sk(k,j-1,i) ) &2386 - ( 5.0_wp * ibit5 * adv_sca_5 &2387 + ibit4 * adv_sca_3 &2388 ) * &2389 ( sk(k,j+1,i) - sk(k,j-2,i) ) &2390 + ( ibit5 * adv_sca_5 &2391 ) * &2392 ( sk(k,j+2,i) - sk(k,j-3,i) ) &2393 )2394 2395 ENDDO2396 !2397 !-- Above to the top of the highest topography. No degradation necessary.2398 DO k = nzb_max+1, nzt2399 2400 v_comp = v(k,j,i) - v_gtrans2401 swap_flux_y_local(k) = v_comp * ( &2402 37.0_wp * ( sk(k,j,i) + sk(k,j-1,i) ) &2403 - 8.0_wp * ( sk(k,j+1,i) + sk(k,j-2,i) ) &2404 + ( sk(k,j+2,i) + sk(k,j-3,i) ) &2405 ) * adv_sca_52406 swap_diss_y_local(k) = -ABS( v_comp ) * ( &2407 10.0_wp * ( sk(k,j,i) - sk(k,j-1,i) ) &2408 - 5.0_wp * ( sk(k,j+1,i) - sk(k,j-2,i) ) &2409 + sk(k,j+2,i) - sk(k,j-3,i) &2410 ) * adv_sca_52411 2412 ENDDO2413 2414 DO j = nys, nyn2415 2416 flux_t(0) = 0.0_wp2417 diss_t(0) = 0.0_wp2418 flux_d = 0.0_wp2419 diss_d = 0.0_wp2420 2421 DO k = nzb+1, nzb_max2422 2423 ibit2 = IBITS(wall_flags_0(k,j,i),2,1)2424 ibit1 = IBITS(wall_flags_0(k,j,i),1,1)2425 ibit0 = IBITS(wall_flags_0(k,j,i),0,1)2426 2427 u_comp = u(k,j,i+1) - u_gtrans2428 flux_r(k) = u_comp * ( &2429 ( 37.0_wp * ibit2 * adv_sca_5 &2430 + 7.0_wp * ibit1 * adv_sca_3 &2431 + ibit0 * adv_sca_1 &2432 ) * &2433 ( sk(k,j,i+1) + sk(k,j,i) ) &2434 - ( 8.0_wp * ibit2 * adv_sca_5 &2435 + ibit1 * adv_sca_3 &2436 ) * &2437 ( sk(k,j,i+2) + sk(k,j,i-1) ) &2438 + ( ibit2 * adv_sca_5 &2439 ) * &2440 ( sk(k,j,i+3) + sk(k,j,i-2) ) &2441 )2442 2443 diss_r(k) = -ABS( u_comp ) * ( &2444 ( 10.0_wp * ibit2 * adv_sca_5 &2445 + 3.0_wp * ibit1 * adv_sca_3 &2446 + ibit0 * adv_sca_1 &2447 ) * &2448 ( sk(k,j,i+1) - sk(k,j,i) ) &2449 - ( 5.0_wp * ibit2 * adv_sca_5 &2450 + ibit1 * adv_sca_3 &2451 ) * &2452 ( sk(k,j,i+2) - sk(k,j,i-1) ) &2453 + ( ibit2 * adv_sca_5 &2454 ) * &2455 ( sk(k,j,i+3) - sk(k,j,i-2) ) &2456 )2457 2458 ibit5 = IBITS(wall_flags_0(k,j,i),5,1)2459 ibit4 = IBITS(wall_flags_0(k,j,i),4,1)2460 ibit3 = IBITS(wall_flags_0(k,j,i),3,1)2461 2462 v_comp = v(k,j+1,i) - v_gtrans2463 flux_n(k) = v_comp * ( &2464 ( 37.0_wp * ibit5 * adv_sca_5 &2465 + 7.0_wp * ibit4 * adv_sca_3 &2466 + ibit3 * adv_sca_1 &2467 ) * &2468 ( sk(k,j+1,i) + sk(k,j,i) ) &2469 - ( 8.0_wp * ibit5 * adv_sca_5 &2470 + ibit4 * adv_sca_3 &2471 ) * &2472 ( sk(k,j+2,i) + sk(k,j-1,i) ) &2473 + ( ibit5 * adv_sca_5 &2474 ) * &2475 ( sk(k,j+3,i) + sk(k,j-2,i) ) &2476 )2477 2478 diss_n(k) = -ABS( v_comp ) * ( &2479 ( 10.0_wp * ibit5 * adv_sca_5 &2480 + 3.0_wp * ibit4 * adv_sca_3 &2481 + ibit3 * adv_sca_1 &2482 ) * &2483 ( sk(k,j+1,i) - sk(k,j,i) ) &2484 - ( 5.0_wp * ibit5 * adv_sca_5 &2485 + ibit4 * adv_sca_3 &2486 ) * &2487 ( sk(k,j+2,i) - sk(k,j-1,i) ) &2488 + ( ibit5 * adv_sca_5 &2489 ) * &2490 ( sk(k,j+3,i) - sk(k,j-2,i) ) &2491 )2492 !2493 !-- k index has to be modified near bottom and top, else array2494 !-- subscripts will be exceeded.2495 ibit8 = IBITS(wall_flags_0(k,j,i),8,1)2496 ibit7 = IBITS(wall_flags_0(k,j,i),7,1)2497 ibit6 = IBITS(wall_flags_0(k,j,i),6,1)2498 2499 k_ppp = k + 3 * ibit82500 k_pp = k + 2 * ( 1 - ibit6 )2501 k_mm = k - 2 * ibit82502 2503 2504 flux_t(k) = w(k,j,i) * ( &2505 ( 37.0_wp * ibit8 * adv_sca_5 &2506 + 7.0_wp * ibit7 * adv_sca_3 &2507 + ibit6 * adv_sca_1 &2508 ) * &2509 ( sk(k+1,j,i) + sk(k,j,i) ) &2510 - ( 8.0_wp * ibit8 * adv_sca_5 &2511 + ibit7 * adv_sca_3 &2512 ) * &2513 ( sk(k_pp,j,i) + sk(k-1,j,i) ) &2514 + ( ibit8 * adv_sca_5 &2515 ) * ( sk(k_ppp,j,i)+ sk(k_mm,j,i) ) &2516 )2517 2518 diss_t(k) = -ABS( w(k,j,i) ) * ( &2519 ( 10.0_wp * ibit8 * adv_sca_5 &2520 + 3.0_wp * ibit7 * adv_sca_3 &2521 + ibit6 * adv_sca_1 &2522 ) * &2523 ( sk(k+1,j,i) - sk(k,j,i) ) &2524 - ( 5.0_wp * ibit8 * adv_sca_5 &2525 + ibit7 * adv_sca_3 &2526 ) * &2527 ( sk(k_pp,j,i) - sk(k-1,j,i) ) &2528 + ( ibit8 * adv_sca_5 &2529 ) * &2530 ( sk(k_ppp,j,i) - sk(k_mm,j,i) ) &2531 )2532 !2533 !-- Calculate the divergence of the velocity field. A respective2534 !-- correction is needed to overcome numerical instabilities caused2535 !-- by a not sufficient reduction of divergences near topography.2536 div = ( u(k,j,i+1) - u(k,j,i) ) * ddx &2537 + ( v(k,j+1,i) - v(k,j,i) ) * ddy &2538 + ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)2539 2540 tend(k,j,i) = tend(k,j,i) - ( &2541 ( flux_r(k) + diss_r(k) - swap_flux_x_local(k,j) - &2542 swap_diss_x_local(k,j) ) * ddx &2543 + ( flux_n(k) + diss_n(k) - swap_flux_y_local(k) - &2544 swap_diss_y_local(k) ) * ddy &2545 + ( flux_t(k) + diss_t(k) - flux_d - diss_d &2546 ) * ddzw(k) &2547 ) + sk(k,j,i) * div2548 2549 swap_flux_y_local(k) = flux_n(k)2550 swap_diss_y_local(k) = diss_n(k)2551 swap_flux_x_local(k,j) = flux_r(k)2552 swap_diss_x_local(k,j) = diss_r(k)2553 flux_d = flux_t(k)2554 diss_d = diss_t(k)2555 2556 ENDDO2557 2558 DO k = nzb_max+1, nzt2559 2560 u_comp = u(k,j,i+1) - u_gtrans2561 flux_r(k) = u_comp * ( &2562 37.0_wp * ( sk(k,j,i+1) + sk(k,j,i) ) &2563 - 8.0_wp * ( sk(k,j,i+2) + sk(k,j,i-1) ) &2564 + ( sk(k,j,i+3) + sk(k,j,i-2) ) ) * adv_sca_52565 diss_r(k) = -ABS( u_comp ) * ( &2566 10.0_wp * ( sk(k,j,i+1) - sk(k,j,i) ) &2567 - 5.0_wp * ( sk(k,j,i+2) - sk(k,j,i-1) ) &2568 + ( sk(k,j,i+3) - sk(k,j,i-2) ) ) * adv_sca_52569 2570 v_comp = v(k,j+1,i) - v_gtrans2571 flux_n(k) = v_comp * ( &2572 37.0_wp * ( sk(k,j+1,i) + sk(k,j,i) ) &2573 - 8.0_wp * ( sk(k,j+2,i) + sk(k,j-1,i) ) &2574 + ( sk(k,j+3,i) + sk(k,j-2,i) ) ) * adv_sca_52575 diss_n(k) = -ABS( v_comp ) * ( &2576 10.0_wp * ( sk(k,j+1,i) - sk(k,j,i) ) &2577 - 5.0_wp * ( sk(k,j+2,i) - sk(k,j-1,i) ) &2578 + ( sk(k,j+3,i) - sk(k,j-2,i) ) ) * adv_sca_52579 !2580 !-- k index has to be modified near bottom and top, else array2581 !-- subscripts will be exceeded.2582 ibit8 = IBITS(wall_flags_0(k,j,i),8,1)2583 ibit7 = IBITS(wall_flags_0(k,j,i),7,1)2584 ibit6 = IBITS(wall_flags_0(k,j,i),6,1)2585 2586 k_ppp = k + 3 * ibit82587 k_pp = k + 2 * ( 1 - ibit6 )2588 k_mm = k - 2 * ibit82589 2590 2591 flux_t(k) = w(k,j,i) * ( &2592 ( 37.0_wp * ibit8 * adv_sca_5 &2593 + 7.0_wp * ibit7 * adv_sca_3 &2594 + ibit6 * adv_sca_1 &2595 ) * &2596 ( sk(k+1,j,i) + sk(k,j,i) ) &2597 - ( 8.0_wp * ibit8 * adv_sca_5 &2598 + ibit7 * adv_sca_3 &2599 ) * &2600 ( sk(k_pp,j,i) + sk(k-1,j,i) ) &2601 + ( ibit8 * adv_sca_5 &2602 ) * ( sk(k_ppp,j,i)+ sk(k_mm,j,i) ) &2603 )2604 2605 diss_t(k) = -ABS( w(k,j,i) ) * ( &2606 ( 10.0_wp * ibit8 * adv_sca_5 &2607 + 3.0_wp * ibit7 * adv_sca_3 &2608 + ibit6 * adv_sca_1 &2609 ) * &2610 ( sk(k+1,j,i) - sk(k,j,i) ) &2611 - ( 5.0_wp * ibit8 * adv_sca_5 &2612 + ibit7 * adv_sca_3 &2613 ) * &2614 ( sk(k_pp,j,i) - sk(k-1,j,i) ) &2615 + ( ibit8 * adv_sca_5 &2616 ) * &2617 ( sk(k_ppp,j,i) - sk(k_mm,j,i) ) &2618 )2619 !2620 !-- Calculate the divergence of the velocity field. A respective2621 !-- correction is needed to overcome numerical instabilities introduced2622 !-- by a not sufficient reduction of divergences near topography.2623 div = ( u(k,j,i+1) - u(k,j,i) ) * ddx &2624 + ( v(k,j+1,i) - v(k,j,i) ) * ddy &2625 + ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)2626 2627 tend(k,j,i) = tend(k,j,i) - ( &2628 ( flux_r(k) + diss_r(k) - swap_flux_x_local(k,j) - &2629 swap_diss_x_local(k,j) ) * ddx &2630 + ( flux_n(k) + diss_n(k) - swap_flux_y_local(k) - &2631 swap_diss_y_local(k) ) * ddy &2632 + ( flux_t(k) + diss_t(k) - flux_d - diss_d &2633 ) * ddzw(k) &2634 ) + sk(k,j,i) * div2635 2636 swap_flux_y_local(k) = flux_n(k)2637 swap_diss_y_local(k) = diss_n(k)2638 swap_flux_x_local(k,j) = flux_r(k)2639 swap_diss_x_local(k,j) = diss_r(k)2640 flux_d = flux_t(k)2641 diss_d = diss_t(k)2642 2643 ENDDO2644 !2645 !-- evaluation of statistics2646 SELECT CASE ( sk_char )2647 2648 CASE ( 'pt' )2649 DO k = nzb, nzt2650 sums_wspts_ws_l(k,tn) = sums_wspts_ws_l(k,tn) &2651 + ( flux_t(k) + diss_t(k) ) &2652 * weight_substep(intermediate_timestep_count)2653 ENDDO2654 CASE ( 'sa' )2655 DO k = nzb, nzt2656 sums_wssas_ws_l(k,tn) = sums_wssas_ws_l(k,tn) &2657 + ( flux_t(k) + diss_t(k) ) &2658 * weight_substep(intermediate_timestep_count)2659 ENDDO2660 CASE ( 'q' )2661 DO k = nzb, nzt2662 sums_wsqs_ws_l(k,tn) = sums_wsqs_ws_l(k,tn) &2663 + ( flux_t(k) + diss_t(k) ) &2664 * weight_substep(intermediate_timestep_count)2665 ENDDO2666 CASE ( 'qr' )2667 DO k = nzb, nzt2668 sums_wsqrs_ws_l(k,tn) = sums_wsqrs_ws_l(k,tn) &2669 + ( flux_t(k) + diss_t(k) ) &2670 * weight_substep(intermediate_timestep_count)2671 ENDDO2672 CASE ( 'nr' )2673 DO k = nzb, nzt2674 sums_wsnrs_ws_l(k,tn) = sums_wsnrs_ws_l(k,tn) &2675 + ( flux_t(k) + diss_t(k) ) &2676 * weight_substep(intermediate_timestep_count)2677 ENDDO2678 2679 END SELECT2680 2681 ENDDO2682 ENDDO2683 2684 END SUBROUTINE advec_s_ws2685 2686 2687 !------------------------------------------------------------------------------!2688 ! Scalar advection - Call for all grid points - accelerator version2689 !------------------------------------------------------------------------------!2690 SUBROUTINE advec_s_ws_acc ( sk, sk_char )2691 2692 USE arrays_3d, &2693 ONLY: ddzw, tend, u, v, w2694 2695 USE constants, &2696 ONLY: adv_sca_1, adv_sca_3, adv_sca_52697 2698 USE control_parameters, &2699 ONLY: intermediate_timestep_count, u_gtrans, v_gtrans2700 2701 USE grid_variables, &2702 ONLY: ddx, ddy2703 2704 USE indices, &2705 ONLY: i_left, i_right, j_north, j_south, nxlg, nxrg, nyng, nysg, &2706 nzb, nzb_max, nzt, wall_flags_02707 2708 USE kinds2709 2710 ! USE statistics, &2711 ! ONLY: sums_wspts_ws_l, sums_wsqs_ws_l, sums_wssas_ws_l, &2712 ! sums_wsqrs_ws_l, sums_wsnrs_ws_l, weight_substep2713 2714 IMPLICIT NONE2715 2716 CHARACTER (LEN = *), INTENT(IN) :: sk_char !:2717 2718 2552 INTEGER(iwp) :: i !: 2719 2553 INTEGER(iwp) :: ibit0 !: … … 2733 2567 INTEGER(iwp) :: k_ppp !: 2734 2568 INTEGER(iwp) :: tn = 0 !: 2569 2570 #if defined( __nopointer ) 2571 REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: sk !: 2572 #else 2573 REAL(wp), DIMENSION(:,:,:), POINTER :: sk !: 2574 #endif 2575 2576 REAL(wp) :: diss_d !: 2577 REAL(wp) :: div !: 2578 REAL(wp) :: flux_d !: 2579 REAL(wp) :: fd_1 !: 2580 REAL(wp) :: fl_1 !: 2581 REAL(wp) :: fn_1 !: 2582 REAL(wp) :: fr_1 !: 2583 REAL(wp) :: fs_1 !: 2584 REAL(wp) :: ft_1 !: 2585 REAL(wp) :: phi_d !: 2586 REAL(wp) :: phi_l !: 2587 REAL(wp) :: phi_n !: 2588 REAL(wp) :: phi_r !: 2589 REAL(wp) :: phi_s !: 2590 REAL(wp) :: phi_t !: 2591 REAL(wp) :: rd !: 2592 REAL(wp) :: rl !: 2593 REAL(wp) :: rn !: 2594 REAL(wp) :: rr !: 2595 REAL(wp) :: rs !: 2596 REAL(wp) :: rt !: 2597 REAL(wp) :: u_comp !: 2598 REAL(wp) :: v_comp !: 2599 2600 REAL(wp), DIMENSION(nzb:nzt) :: diss_n !: 2601 REAL(wp), DIMENSION(nzb:nzt) :: diss_r !: 2602 REAL(wp), DIMENSION(nzb:nzt) :: diss_t !: 2603 REAL(wp), DIMENSION(nzb:nzt) :: flux_n !: 2604 REAL(wp), DIMENSION(nzb:nzt) :: flux_r !: 2605 REAL(wp), DIMENSION(nzb:nzt) :: flux_t !: 2606 2607 REAL(wp), DIMENSION(nzb+1:nzt) :: swap_diss_y_local !: 2608 REAL(wp), DIMENSION(nzb+1:nzt) :: swap_flux_y_local !: 2609 2610 REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) :: swap_diss_x_local !: 2611 REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) :: swap_flux_x_local !: 2612 2613 2614 ! 2615 !-- Compute the fluxes for the whole left boundary of the processor domain. 2616 i = nxl 2617 DO j = nys, nyn 2618 2619 DO k = nzb+1, nzb_max 2620 2621 ibit2 = IBITS(wall_flags_0(k,j,i),2,1) 2622 ibit1 = IBITS(wall_flags_0(k,j,i),1,1) 2623 ibit0 = IBITS(wall_flags_0(k,j,i),0,1) 2624 2625 u_comp = u(k,j,i) - u_gtrans 2626 swap_flux_x_local(k,j) = u_comp * ( & 2627 ( 37.0_wp * ibit2 * adv_sca_5 & 2628 + 7.0_wp * ibit1 * adv_sca_3 & 2629 + ibit0 * adv_sca_1 & 2630 ) * & 2631 ( sk(k,j,i) + sk(k,j,i-1) ) & 2632 - ( 8.0_wp * ibit2 * adv_sca_5 & 2633 + ibit1 * adv_sca_3 & 2634 ) * & 2635 ( sk(k,j,i+1) + sk(k,j,i-2) ) & 2636 + ( ibit2 * adv_sca_5 & 2637 ) * & 2638 ( sk(k,j,i+2) + sk(k,j,i-3) ) & 2639 ) 2640 2641 swap_diss_x_local(k,j) = -ABS( u_comp ) * ( & 2642 ( 10.0_wp * ibit2 * adv_sca_5 & 2643 + 3.0_wp * ibit1 * adv_sca_3 & 2644 + ibit0 * adv_sca_1 & 2645 ) * & 2646 ( sk(k,j,i) - sk(k,j,i-1) ) & 2647 - ( 5.0_wp * ibit2 * adv_sca_5 & 2648 + ibit1 * adv_sca_3 & 2649 ) * & 2650 ( sk(k,j,i+1) - sk(k,j,i-2) ) & 2651 + ( ibit2 * adv_sca_5 & 2652 ) * & 2653 ( sk(k,j,i+2) - sk(k,j,i-3) ) & 2654 ) 2655 2656 ENDDO 2657 2658 DO k = nzb_max+1, nzt 2659 2660 u_comp = u(k,j,i) - u_gtrans 2661 swap_flux_x_local(k,j) = u_comp * ( & 2662 37.0_wp * ( sk(k,j,i) + sk(k,j,i-1) ) & 2663 - 8.0_wp * ( sk(k,j,i+1) + sk(k,j,i-2) ) & 2664 + ( sk(k,j,i+2) + sk(k,j,i-3) ) & 2665 ) * adv_sca_5 2666 2667 swap_diss_x_local(k,j) = -ABS( u_comp ) * ( & 2668 10.0_wp * ( sk(k,j,i) - sk(k,j,i-1) ) & 2669 - 5.0_wp * ( sk(k,j,i+1) - sk(k,j,i-2) ) & 2670 + ( sk(k,j,i+2) - sk(k,j,i-3) ) & 2671 ) * adv_sca_5 2672 2673 ENDDO 2674 2675 ENDDO 2676 2677 DO i = nxl, nxr 2678 2679 j = nys 2680 DO k = nzb+1, nzb_max 2681 2682 ibit5 = IBITS(wall_flags_0(k,j,i),5,1) 2683 ibit4 = IBITS(wall_flags_0(k,j,i),4,1) 2684 ibit3 = IBITS(wall_flags_0(k,j,i),3,1) 2685 2686 v_comp = v(k,j,i) - v_gtrans 2687 swap_flux_y_local(k) = v_comp * ( & 2688 ( 37.0_wp * ibit5 * adv_sca_5 & 2689 + 7.0_wp * ibit4 * adv_sca_3 & 2690 + ibit3 * adv_sca_1 & 2691 ) * & 2692 ( sk(k,j,i) + sk(k,j-1,i) ) & 2693 - ( 8.0_wp * ibit5 * adv_sca_5 & 2694 + ibit4 * adv_sca_3 & 2695 ) * & 2696 ( sk(k,j+1,i) + sk(k,j-2,i) ) & 2697 + ( ibit5 * adv_sca_5 & 2698 ) * & 2699 ( sk(k,j+2,i) + sk(k,j-3,i) ) & 2700 ) 2701 2702 swap_diss_y_local(k) = -ABS( v_comp ) * ( & 2703 ( 10.0_wp * ibit5 * adv_sca_5 & 2704 + 3.0_wp * ibit4 * adv_sca_3 & 2705 + ibit3 * adv_sca_1 & 2706 ) * & 2707 ( sk(k,j,i) - sk(k,j-1,i) ) & 2708 - ( 5.0_wp * ibit5 * adv_sca_5 & 2709 + ibit4 * adv_sca_3 & 2710 ) * & 2711 ( sk(k,j+1,i) - sk(k,j-2,i) ) & 2712 + ( ibit5 * adv_sca_5 & 2713 ) * & 2714 ( sk(k,j+2,i) - sk(k,j-3,i) ) & 2715 ) 2716 2717 ENDDO 2718 ! 2719 !-- Above to the top of the highest topography. No degradation necessary. 2720 DO k = nzb_max+1, nzt 2721 2722 v_comp = v(k,j,i) - v_gtrans 2723 swap_flux_y_local(k) = v_comp * ( & 2724 37.0_wp * ( sk(k,j,i) + sk(k,j-1,i) ) & 2725 - 8.0_wp * ( sk(k,j+1,i) + sk(k,j-2,i) ) & 2726 + ( sk(k,j+2,i) + sk(k,j-3,i) ) & 2727 ) * adv_sca_5 2728 swap_diss_y_local(k) = -ABS( v_comp ) * ( & 2729 10.0_wp * ( sk(k,j,i) - sk(k,j-1,i) ) & 2730 - 5.0_wp * ( sk(k,j+1,i) - sk(k,j-2,i) ) & 2731 + sk(k,j+2,i) - sk(k,j-3,i) & 2732 ) * adv_sca_5 2733 2734 ENDDO 2735 2736 DO j = nys, nyn 2737 2738 flux_t(0) = 0.0_wp 2739 diss_t(0) = 0.0_wp 2740 flux_d = 0.0_wp 2741 diss_d = 0.0_wp 2742 2743 DO k = nzb+1, nzb_max 2744 2745 ibit2 = IBITS(wall_flags_0(k,j,i),2,1) 2746 ibit1 = IBITS(wall_flags_0(k,j,i),1,1) 2747 ibit0 = IBITS(wall_flags_0(k,j,i),0,1) 2748 2749 u_comp = u(k,j,i+1) - u_gtrans 2750 flux_r(k) = u_comp * ( & 2751 ( 37.0_wp * ibit2 * adv_sca_5 & 2752 + 7.0_wp * ibit1 * adv_sca_3 & 2753 + ibit0 * adv_sca_1 & 2754 ) * & 2755 ( sk(k,j,i+1) + sk(k,j,i) ) & 2756 - ( 8.0_wp * ibit2 * adv_sca_5 & 2757 + ibit1 * adv_sca_3 & 2758 ) * & 2759 ( sk(k,j,i+2) + sk(k,j,i-1) ) & 2760 + ( ibit2 * adv_sca_5 & 2761 ) * & 2762 ( sk(k,j,i+3) + sk(k,j,i-2) ) & 2763 ) 2764 2765 diss_r(k) = -ABS( u_comp ) * ( & 2766 ( 10.0_wp * ibit2 * adv_sca_5 & 2767 + 3.0_wp * ibit1 * adv_sca_3 & 2768 + ibit0 * adv_sca_1 & 2769 ) * & 2770 ( sk(k,j,i+1) - sk(k,j,i) ) & 2771 - ( 5.0_wp * ibit2 * adv_sca_5 & 2772 + ibit1 * adv_sca_3 & 2773 ) * & 2774 ( sk(k,j,i+2) - sk(k,j,i-1) ) & 2775 + ( ibit2 * adv_sca_5 & 2776 ) * & 2777 ( sk(k,j,i+3) - sk(k,j,i-2) ) & 2778 ) 2779 2780 ibit5 = IBITS(wall_flags_0(k,j,i),5,1) 2781 ibit4 = IBITS(wall_flags_0(k,j,i),4,1) 2782 ibit3 = IBITS(wall_flags_0(k,j,i),3,1) 2783 2784 v_comp = v(k,j+1,i) - v_gtrans 2785 flux_n(k) = v_comp * ( & 2786 ( 37.0_wp * ibit5 * adv_sca_5 & 2787 + 7.0_wp * ibit4 * adv_sca_3 & 2788 + ibit3 * adv_sca_1 & 2789 ) * & 2790 ( sk(k,j+1,i) + sk(k,j,i) ) & 2791 - ( 8.0_wp * ibit5 * adv_sca_5 & 2792 + ibit4 * adv_sca_3 & 2793 ) * & 2794 ( sk(k,j+2,i) + sk(k,j-1,i) ) & 2795 + ( ibit5 * adv_sca_5 & 2796 ) * & 2797 ( sk(k,j+3,i) + sk(k,j-2,i) ) & 2798 ) 2799 2800 diss_n(k) = -ABS( v_comp ) * ( & 2801 ( 10.0_wp * ibit5 * adv_sca_5 & 2802 + 3.0_wp * ibit4 * adv_sca_3 & 2803 + ibit3 * adv_sca_1 & 2804 ) * & 2805 ( sk(k,j+1,i) - sk(k,j,i) ) & 2806 - ( 5.0_wp * ibit5 * adv_sca_5 & 2807 + ibit4 * adv_sca_3 & 2808 ) * & 2809 ( sk(k,j+2,i) - sk(k,j-1,i) ) & 2810 + ( ibit5 * adv_sca_5 & 2811 ) * & 2812 ( sk(k,j+3,i) - sk(k,j-2,i) ) & 2813 ) 2814 ! 2815 !-- k index has to be modified near bottom and top, else array 2816 !-- subscripts will be exceeded. 2817 ibit8 = IBITS(wall_flags_0(k,j,i),8,1) 2818 ibit7 = IBITS(wall_flags_0(k,j,i),7,1) 2819 ibit6 = IBITS(wall_flags_0(k,j,i),6,1) 2820 2821 k_ppp = k + 3 * ibit8 2822 k_pp = k + 2 * ( 1 - ibit6 ) 2823 k_mm = k - 2 * ibit8 2824 2825 2826 flux_t(k) = w(k,j,i) * ( & 2827 ( 37.0_wp * ibit8 * adv_sca_5 & 2828 + 7.0_wp * ibit7 * adv_sca_3 & 2829 + ibit6 * adv_sca_1 & 2830 ) * & 2831 ( sk(k+1,j,i) + sk(k,j,i) ) & 2832 - ( 8.0_wp * ibit8 * adv_sca_5 & 2833 + ibit7 * adv_sca_3 & 2834 ) * & 2835 ( sk(k_pp,j,i) + sk(k-1,j,i) ) & 2836 + ( ibit8 * adv_sca_5 & 2837 ) * ( sk(k_ppp,j,i)+ sk(k_mm,j,i) ) & 2838 ) 2839 2840 diss_t(k) = -ABS( w(k,j,i) ) * ( & 2841 ( 10.0_wp * ibit8 * adv_sca_5 & 2842 + 3.0_wp * ibit7 * adv_sca_3 & 2843 + ibit6 * adv_sca_1 & 2844 ) * & 2845 ( sk(k+1,j,i) - sk(k,j,i) ) & 2846 - ( 5.0_wp * ibit8 * adv_sca_5 & 2847 + ibit7 * adv_sca_3 & 2848 ) * & 2849 ( sk(k_pp,j,i) - sk(k-1,j,i) ) & 2850 + ( ibit8 * adv_sca_5 & 2851 ) * & 2852 ( sk(k_ppp,j,i) - sk(k_mm,j,i) ) & 2853 ) 2854 ! 2855 !-- Apply monotonic adjustment. 2856 IF ( monotonic_adjustment ) THEN 2857 ! 2858 !-- At first, calculate first order fluxes. 2859 u_comp = u(k,j,i) - u_gtrans 2860 fl_1 = ( u_comp * ( sk(k,j,i) + sk(k,j,i-1) ) & 2861 -ABS( u_comp ) * ( sk(k,j,i) - sk(k,j,i-1) ) & 2862 ) * adv_sca_1 2863 2864 u_comp = u(k,j,i+1) - u_gtrans 2865 fr_1 = ( u_comp * ( sk(k,j,i+1) + sk(k,j,i) ) & 2866 -ABS( u_comp ) * ( sk(k,j,i+1) - sk(k,j,i) ) & 2867 ) * adv_sca_1 2868 2869 v_comp = v(k,j,i) - v_gtrans 2870 fs_1 = ( v_comp * ( sk(k,j,i) + sk(k,j-1,i) ) & 2871 -ABS( v_comp ) * ( sk(k,j,i) - sk(k,j-1,i) ) & 2872 ) * adv_sca_1 2873 2874 v_comp = v(k,j+1,i) - v_gtrans 2875 fn_1 = ( v_comp * ( sk(k,j+1,i) + sk(k,j,i) ) & 2876 -ABS( v_comp ) * ( sk(k,j+1,i) - sk(k,j,i) ) & 2877 ) * adv_sca_1 2878 2879 fd_1 = ( w(k-1,j,i) * ( sk(k,j,i) + sk(k-1,j,i) ) & 2880 -ABS( w(k-1,j,i) ) * ( sk(k,j,i) - sk(k-1,j,i) ) & 2881 ) * adv_sca_1 2882 2883 ft_1 = ( w(k,j,i) * ( sk(k+1,j,i) + sk(k,j,i) ) & 2884 -ABS( w(k,j,i) ) * ( sk(k+1,j,i) - sk(k,j,i) ) & 2885 ) * adv_sca_1 2886 ! 2887 !-- Calculate ratio of upwind gradients. Note, Min/Max is just 2888 !-- to avoid if statements. 2889 rl = ( MAX( 0.0, u(k,j,i) - u_gtrans ) * & 2890 ABS( ( sk(k,j,i-1) - sk(k,j,i-2) ) /& 2891 ( sk(k,j,i) - sk(k,j,i-1) + 1E-20_wp ) & 2892 ) + & 2893 MIN( 0.0, u(k,j,i) - u_gtrans ) * & 2894 ABS( ( sk(k,j,i) - sk(k,j,i+1) ) /& 2895 ( sk(k,j,i-1) - sk(k,j,i) + 1E-20_wp ) & 2896 ) & 2897 ) / ABS( u(k,j,i) - u_gtrans + 1E-20_wp ) 2898 2899 rr = ( MAX( 0.0, u(k,j,i+1) - u_gtrans ) * & 2900 ABS( ( sk(k,j,i) - sk(k,j,i-1) ) /& 2901 ( sk(k,j,i+1) - sk(k,j,i) + 1E-20_wp ) & 2902 ) + & 2903 MIN( 0.0, u(k,j,i+1) - u_gtrans ) * & 2904 ABS( ( sk(k,j,i+1) - sk(k,j,i+2) ) /& 2905 ( sk(k,j,i) - sk(k,j,i+1) + 1E-20_wp ) & 2906 ) & 2907 ) / ABS( u(k,j,i+1) - u_gtrans + 1E-20_wp ) 2908 2909 rs = ( MAX( 0.0, v(k,j,i) - v_gtrans ) * & 2910 ABS( ( sk(k,j-1,i) - sk(k,j-2,i) ) /& 2911 ( sk(k,j,i) - sk(k,j-1,i) + 1E-20_wp ) & 2912 ) + & 2913 MIN( 0.0, v(k,j,i) - v_gtrans ) * & 2914 ABS( ( sk(k,j,i) - sk(k,j+1,i) ) /& 2915 ( sk(k,j-1,i) - sk(k,j,i) + 1E-20_wp ) & 2916 ) & 2917 ) / ABS( v(k,j,i) - v_gtrans + 1E-20_wp ) 2918 2919 rn = ( MAX( 0.0, v(k,j+1,i) - v_gtrans ) * & 2920 ABS( ( sk(k,j,i) - sk(k,j-1,i) ) /& 2921 ( sk(k,j+1,i) - sk(k,j,i) + 1E-20_wp ) & 2922 ) + & 2923 MIN( 0.0, v(k,j+1,i) - v_gtrans ) * & 2924 ABS( ( sk(k,j+1,i) - sk(k,j+2,i) ) /& 2925 ( sk(k,j,i) - sk(k,j+1,i) + 1E-20_wp ) & 2926 ) & 2927 ) / ABS( v(k,j+1,i) - v_gtrans + 1E-20_wp ) 2928 ! 2929 !-- Reuse k_mm and compute k_mmm for the vertical gradient ratios. 2930 !-- Note, for vertical advection below the third grid point above 2931 !-- surface ( or below the model top) rd and rt are set to 0, i.e. 2932 !-- use of first order scheme is enforced. 2933 k_mmm = k - 3 * ibit8 2934 2935 rd = ( MAX( 0.0, w(k-1,j,i) ) * & 2936 ABS( ( sk(k_mm,j,i) - sk(k_mmm,j,i) ) / & 2937 ( sk(k-1,j,i) - sk(k_mm,j,i) + 1E-20_wp ) & 2938 ) + & 2939 MIN( 0.0, w(k-1,j,i) ) * & 2940 ABS( ( sk(k-1,j,i) - sk(k,j,i) ) / & 2941 ( sk(k_mm,j,i) - sk(k-1,j,i) + 1E-20_wp ) & 2942 ) & 2943 ) * ibit8 / ABS( w(k-1,j,i) + 1E-20_wp ) 2944 2945 rt = ( MAX( 0.0, w(k,j,i) ) * & 2946 ABS( ( sk(k,j,i) - sk(k-1,j,i) ) / & 2947 ( sk(k+1,j,i) - sk(k,j,i) + 1E-20_wp ) & 2948 ) + & 2949 MIN( 0.0, w(k,j,i) ) * & 2950 ABS( ( sk(k+1,j,i) - sk(k_pp,j,i) ) / & 2951 ( sk(k,j,i) - sk(k+1,j,i) + 1E-20_wp ) & 2952 ) & 2953 ) * ibit8 / ABS( w(k,j,i) + 1E-20_wp ) 2954 ! 2955 !-- Calculate empirical limiter function (van Albada2 limiter). 2956 phi_l = MIN( 1.0_wp, ( 2.0_wp * rl ) / & 2957 ( rl**2 + 1.0_wp ) ) 2958 phi_r = MIN( 1.0_wp, ( 2.0_wp * rr ) / & 2959 ( rr**2 + 1.0_wp ) ) 2960 phi_s = MIN( 1.0_wp, ( 2.0_wp * rs ) / & 2961 ( rs**2 + 1.0_wp ) ) 2962 phi_n = MIN( 1.0_wp, ( 2.0_wp * rn ) / & 2963 ( rn**2 + 1.0_wp ) ) 2964 phi_d = MIN( 1.0_wp, ( 2.0_wp * rd ) / & 2965 ( rd**2 + 1.0_wp ) ) 2966 phi_t = MIN( 1.0_wp, ( 2.0_wp * rt ) / & 2967 ( rt**2 + 1.0_wp ) ) 2968 ! 2969 !-- Calculate the resulting monotone flux. 2970 swap_flux_x_local(k,j) = fl_1 - phi_l * & 2971 ( fl_1 - swap_flux_x_local(k,j) ) 2972 flux_r(k) = fr_1 - phi_r * & 2973 ( fr_1 - flux_r(k) ) 2974 swap_flux_y_local(k) = fs_1 - phi_s * & 2975 ( fs_1 - swap_flux_y_local(k) ) 2976 flux_n(k) = fn_1 - phi_n * & 2977 ( fn_1 - flux_n(k) ) 2978 flux_d = fd_1 - phi_d * & 2979 ( fd_1 - flux_d ) 2980 flux_t(k) = ft_1 - phi_t * & 2981 ( ft_1 - flux_t(k) ) 2982 ! 2983 !-- Moreover, modify dissipation flux according to the limiter. 2984 swap_diss_x_local(k,j) = swap_diss_x_local(k,j) * phi_l 2985 diss_r(k) = diss_r(k) * phi_r 2986 swap_diss_y_local(k) = swap_diss_y_local(k) * phi_s 2987 diss_n(k) = diss_n(k) * phi_n 2988 diss_d = diss_d * phi_d 2989 diss_t(k) = diss_t(k) * phi_t 2990 2991 ENDIF 2992 ! 2993 !-- Calculate the divergence of the velocity field. A respective 2994 !-- correction is needed to overcome numerical instabilities caused 2995 !-- by a not sufficient reduction of divergences near topography. 2996 div = ( u(k,j,i+1) - u(k,j,i) ) * ddx & 2997 + ( v(k,j+1,i) - v(k,j,i) ) * ddy & 2998 + ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) 2999 3000 tend(k,j,i) = tend(k,j,i) - ( & 3001 ( flux_r(k) + diss_r(k) - swap_flux_x_local(k,j) - & 3002 swap_diss_x_local(k,j) ) * ddx & 3003 + ( flux_n(k) + diss_n(k) - swap_flux_y_local(k) - & 3004 swap_diss_y_local(k) ) * ddy & 3005 + ( flux_t(k) + diss_t(k) - flux_d - diss_d & 3006 ) * ddzw(k) & 3007 ) + sk(k,j,i) * div 3008 3009 swap_flux_y_local(k) = flux_n(k) 3010 swap_diss_y_local(k) = diss_n(k) 3011 swap_flux_x_local(k,j) = flux_r(k) 3012 swap_diss_x_local(k,j) = diss_r(k) 3013 flux_d = flux_t(k) 3014 diss_d = diss_t(k) 3015 3016 ENDDO 3017 3018 DO k = nzb_max+1, nzt 3019 3020 u_comp = u(k,j,i+1) - u_gtrans 3021 flux_r(k) = u_comp * ( & 3022 37.0_wp * ( sk(k,j,i+1) + sk(k,j,i) ) & 3023 - 8.0_wp * ( sk(k,j,i+2) + sk(k,j,i-1) ) & 3024 + ( sk(k,j,i+3) + sk(k,j,i-2) ) ) * adv_sca_5 3025 diss_r(k) = -ABS( u_comp ) * ( & 3026 10.0_wp * ( sk(k,j,i+1) - sk(k,j,i) ) & 3027 - 5.0_wp * ( sk(k,j,i+2) - sk(k,j,i-1) ) & 3028 + ( sk(k,j,i+3) - sk(k,j,i-2) ) ) * adv_sca_5 3029 3030 v_comp = v(k,j+1,i) - v_gtrans 3031 flux_n(k) = v_comp * ( & 3032 37.0_wp * ( sk(k,j+1,i) + sk(k,j,i) ) & 3033 - 8.0_wp * ( sk(k,j+2,i) + sk(k,j-1,i) ) & 3034 + ( sk(k,j+3,i) + sk(k,j-2,i) ) ) * adv_sca_5 3035 diss_n(k) = -ABS( v_comp ) * ( & 3036 10.0_wp * ( sk(k,j+1,i) - sk(k,j,i) ) & 3037 - 5.0_wp * ( sk(k,j+2,i) - sk(k,j-1,i) ) & 3038 + ( sk(k,j+3,i) - sk(k,j-2,i) ) ) * adv_sca_5 3039 ! 3040 !-- k index has to be modified near bottom and top, else array 3041 !-- subscripts will be exceeded. 3042 ibit8 = IBITS(wall_flags_0(k,j,i),8,1) 3043 ibit7 = IBITS(wall_flags_0(k,j,i),7,1) 3044 ibit6 = IBITS(wall_flags_0(k,j,i),6,1) 3045 3046 k_ppp = k + 3 * ibit8 3047 k_pp = k + 2 * ( 1 - ibit6 ) 3048 k_mm = k - 2 * ibit8 3049 3050 3051 flux_t(k) = w(k,j,i) * ( & 3052 ( 37.0_wp * ibit8 * adv_sca_5 & 3053 + 7.0_wp * ibit7 * adv_sca_3 & 3054 + ibit6 * adv_sca_1 & 3055 ) * & 3056 ( sk(k+1,j,i) + sk(k,j,i) ) & 3057 - ( 8.0_wp * ibit8 * adv_sca_5 & 3058 + ibit7 * adv_sca_3 & 3059 ) * & 3060 ( sk(k_pp,j,i) + sk(k-1,j,i) ) & 3061 + ( ibit8 * adv_sca_5 & 3062 ) * ( sk(k_ppp,j,i)+ sk(k_mm,j,i) ) & 3063 ) 3064 3065 diss_t(k) = -ABS( w(k,j,i) ) * ( & 3066 ( 10.0_wp * ibit8 * adv_sca_5 & 3067 + 3.0_wp * ibit7 * adv_sca_3 & 3068 + ibit6 * adv_sca_1 & 3069 ) * & 3070 ( sk(k+1,j,i) - sk(k,j,i) ) & 3071 - ( 5.0_wp * ibit8 * adv_sca_5 & 3072 + ibit7 * adv_sca_3 & 3073 ) * & 3074 ( sk(k_pp,j,i) - sk(k-1,j,i) ) & 3075 + ( ibit8 * adv_sca_5 & 3076 ) * & 3077 ( sk(k_ppp,j,i) - sk(k_mm,j,i) ) & 3078 ) 3079 ! 3080 !-- Apply monotonic adjustment. 3081 IF ( monotonic_adjustment ) THEN 3082 ! 3083 !-- At first, calculate first order fluxes. 3084 u_comp = u(k,j,i) - u_gtrans 3085 fl_1 = ( u_comp * ( sk(k,j,i) + sk(k,j,i-1) ) & 3086 -ABS( u_comp ) * ( sk(k,j,i) - sk(k,j,i-1) ) & 3087 ) * adv_sca_1 3088 3089 u_comp = u(k,j,i+1) - u_gtrans 3090 fr_1 = ( u_comp * ( sk(k,j,i+1) + sk(k,j,i) ) & 3091 -ABS( u_comp ) * ( sk(k,j,i+1) - sk(k,j,i) ) & 3092 ) * adv_sca_1 3093 3094 v_comp = v(k,j,i) - v_gtrans 3095 fs_1 = ( v_comp * ( sk(k,j,i) + sk(k,j-1,i) ) & 3096 -ABS( v_comp ) * ( sk(k,j,i) - sk(k,j-1,i) ) & 3097 ) * adv_sca_1 3098 3099 v_comp = v(k,j+1,i) - v_gtrans 3100 fn_1 = ( v_comp * ( sk(k,j+1,i) + sk(k,j,i) ) & 3101 -ABS( v_comp ) * ( sk(k,j+1,i) - sk(k,j,i) ) & 3102 ) * adv_sca_1 3103 3104 fd_1 = ( w(k-1,j,i) * ( sk(k,j,i) + sk(k-1,j,i) ) & 3105 -ABS( w(k-1,j,i) ) * ( sk(k,j,i) - sk(k-1,j,i) ) & 3106 ) * adv_sca_1 3107 3108 ft_1 = ( w(k,j,i) * ( sk(k+1,j,i) + sk(k,j,i) ) & 3109 -ABS( w(k,j,i) ) * ( sk(k+1,j,i) - sk(k,j,i) ) & 3110 ) * adv_sca_1 3111 ! 3112 !-- Calculate ratio of upwind gradients. Note, Min/Max is just 3113 !-- to avoid if statements. 3114 rl = ( MAX( 0.0, u(k,j,i) - u_gtrans ) * & 3115 ABS( ( sk(k,j,i-1) - sk(k,j,i-2) ) /& 3116 ( sk(k,j,i) - sk(k,j,i-1) + 1E-20_wp ) & 3117 ) + & 3118 MIN( 0.0, u(k,j,i) - u_gtrans ) * & 3119 ABS( ( sk(k,j,i) - sk(k,j,i+1) ) /& 3120 ( sk(k,j,i-1) - sk(k,j,i) + 1E-20_wp ) & 3121 ) & 3122 ) / ABS( u(k,j,i) - u_gtrans + 1E-20_wp ) 3123 3124 rr = ( MAX( 0.0, u(k,j,i+1) - u_gtrans ) * & 3125 ABS( ( sk(k,j,i) - sk(k,j,i-1) ) /& 3126 ( sk(k,j,i+1) - sk(k,j,i) + 1E-20_wp ) & 3127 ) + & 3128 MIN( 0.0, u(k,j,i+1) - u_gtrans ) * & 3129 ABS( ( sk(k,j,i+1) - sk(k,j,i+2) ) /& 3130 ( sk(k,j,i) - sk(k,j,i+1) + 1E-20_wp ) & 3131 ) & 3132 ) / ABS( u(k,j,i+1) - u_gtrans + 1E-20_wp ) 3133 3134 rs = ( MAX( 0.0, v(k,j,i) - v_gtrans ) * & 3135 ABS( ( sk(k,j-1,i) - sk(k,j-2,i) ) /& 3136 ( sk(k,j,i) - sk(k,j-1,i) + 1E-20_wp ) & 3137 ) + & 3138 MIN( 0.0, v(k,j,i) - v_gtrans ) * & 3139 ABS( ( sk(k,j,i) - sk(k,j+1,i) ) /& 3140 ( sk(k,j-1,i) - sk(k,j,i) + 1E-20_wp ) & 3141 ) & 3142 ) / ABS( v(k,j,i) - v_gtrans + 1E-20_wp ) 3143 3144 rn = ( MAX( 0.0, v(k,j+1,i) - v_gtrans ) * & 3145 ABS( ( sk(k,j,i) - sk(k,j-1,i) ) /& 3146 ( sk(k,j+1,i) - sk(k,j,i) + 1E-20_wp ) & 3147 ) + & 3148 MIN( 0.0, v(k,j+1,i) - v_gtrans ) * & 3149 ABS( ( sk(k,j+1,i) - sk(k,j+2,i) ) /& 3150 ( sk(k,j,i) - sk(k,j+1,i) + 1E-20_wp ) & 3151 ) & 3152 ) / ABS( v(k,j+1,i) - v_gtrans + 1E-20_wp ) 3153 ! 3154 !-- Reuse k_mm and compute k_mmm for the vertical gradient ratios. 3155 !-- Note, for vertical advection below the third grid point above 3156 !-- surface ( or below the model top) rd and rt are set to 0, i.e. 3157 !-- use of first order scheme is enforced. 3158 k_mmm = k - 3 * ibit8 3159 3160 rd = ( MAX( 0.0, w(k-1,j,i) ) * & 3161 ABS( ( sk(k_mm,j,i) - sk(k_mmm,j,i) ) / & 3162 ( sk(k-1,j,i) - sk(k_mm,j,i) + 1E-20_wp ) & 3163 ) + & 3164 MIN( 0.0, w(k-1,j,i) ) * & 3165 ABS( ( sk(k-1,j,i) - sk(k,j,i) ) / & 3166 ( sk(k_mm,j,i) - sk(k-1,j,i) + 1E-20_wp ) & 3167 ) & 3168 ) * ibit8 / ABS( w(k-1,j,i) + 1E-20_wp ) 3169 3170 rt = ( MAX( 0.0, w(k,j,i) ) * & 3171 ABS( ( sk(k,j,i) - sk(k-1,j,i) ) / & 3172 ( sk(k+1,j,i) - sk(k,j,i) + 1E-20_wp ) & 3173 ) + & 3174 MIN( 0.0, w(k,j,i) ) * & 3175 ABS( ( sk(k+1,j,i) - sk(k_pp,j,i) ) / & 3176 ( sk(k,j,i) - sk(k+1,j,i) + 1E-20_wp ) & 3177 ) & 3178 ) * ibit8 / ABS( w(k,j,i) + 1E-20_wp ) 3179 ! 3180 !-- Calculate empirical limiter function (van Albada2 limiter). 3181 phi_l = MIN( 1.0_wp, ( 2.0_wp * rl ) / & 3182 ( rl**2 + 1.0_wp ) ) 3183 phi_r = MIN( 1.0_wp, ( 2.0_wp * rr ) / & 3184 ( rr**2 + 1.0_wp ) ) 3185 phi_s = MIN( 1.0_wp, ( 2.0_wp * rs ) / & 3186 ( rs**2 + 1.0_wp ) ) 3187 phi_n = MIN( 1.0_wp, ( 2.0_wp * rn ) / & 3188 ( rn**2 + 1.0_wp ) ) 3189 phi_d = MIN( 1.0_wp, ( 2.0_wp * rd ) / & 3190 ( rd**2 + 1.0_wp ) ) 3191 phi_t = MIN( 1.0_wp, ( 2.0_wp * rt ) / & 3192 ( rt**2 + 1.0_wp ) ) 3193 ! 3194 !-- Calculate the resulting monotone flux. 3195 swap_flux_x_local(k,j) = fl_1 - phi_l * & 3196 ( fl_1 - swap_flux_x_local(k,j) ) 3197 flux_r(k) = fr_1 - phi_r * & 3198 ( fr_1 - flux_r(k) ) 3199 swap_flux_y_local(k) = fs_1 - phi_s * & 3200 ( fs_1 - swap_flux_y_local(k) ) 3201 flux_n(k) = fn_1 - phi_n * & 3202 ( fn_1 - flux_n(k) ) 3203 flux_d = fd_1 - phi_d * & 3204 ( fd_1 - flux_d ) 3205 flux_t(k) = ft_1 - phi_t * & 3206 ( ft_1 - flux_t(k) ) 3207 ! 3208 !-- Moreover, modify dissipation flux according to the limiter. 3209 swap_diss_x_local(k,j) = swap_diss_x_local(k,j) * phi_l 3210 diss_r(k) = diss_r(k) * phi_r 3211 swap_diss_y_local(k) = swap_diss_y_local(k) * phi_s 3212 diss_n(k) = diss_n(k) * phi_n 3213 diss_d = diss_d * phi_d 3214 diss_t(k) = diss_t(k) * phi_t 3215 3216 ENDIF 3217 ! 3218 !-- Calculate the divergence of the velocity field. A respective 3219 !-- correction is needed to overcome numerical instabilities introduced 3220 !-- by a not sufficient reduction of divergences near topography. 3221 div = ( u(k,j,i+1) - u(k,j,i) ) * ddx & 3222 + ( v(k,j+1,i) - v(k,j,i) ) * ddy & 3223 + ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) 3224 3225 tend(k,j,i) = tend(k,j,i) - ( & 3226 ( flux_r(k) + diss_r(k) - swap_flux_x_local(k,j) - & 3227 swap_diss_x_local(k,j) ) * ddx & 3228 + ( flux_n(k) + diss_n(k) - swap_flux_y_local(k) - & 3229 swap_diss_y_local(k) ) * ddy & 3230 + ( flux_t(k) + diss_t(k) - flux_d - diss_d & 3231 ) * ddzw(k) & 3232 ) + sk(k,j,i) * div 3233 3234 swap_flux_y_local(k) = flux_n(k) 3235 swap_diss_y_local(k) = diss_n(k) 3236 swap_flux_x_local(k,j) = flux_r(k) 3237 swap_diss_x_local(k,j) = diss_r(k) 3238 flux_d = flux_t(k) 3239 diss_d = diss_t(k) 3240 3241 ENDDO 3242 ! 3243 !-- Evaluation of statistics. Please note, in case of using monotone 3244 !-- limiter vertical fluxes will be not calculated by the advective 3245 !-- fluxes. 3246 IF ( .NOT. monotonic_adjustment ) THEN 3247 3248 SELECT CASE ( sk_char ) 3249 3250 CASE ( 'pt' ) 3251 DO k = nzb, nzt 3252 sums_wspts_ws_l(k,tn) = sums_wspts_ws_l(k,tn) & 3253 + ( flux_t(k) + diss_t(k) ) & 3254 * weight_substep(intermediate_timestep_count) 3255 ENDDO 3256 CASE ( 'sa' ) 3257 DO k = nzb, nzt 3258 sums_wssas_ws_l(k,tn) = sums_wssas_ws_l(k,tn) & 3259 + ( flux_t(k) + diss_t(k) ) & 3260 * weight_substep(intermediate_timestep_count) 3261 ENDDO 3262 CASE ( 'q' ) 3263 DO k = nzb, nzt 3264 sums_wsqs_ws_l(k,tn) = sums_wsqs_ws_l(k,tn) & 3265 + ( flux_t(k) + diss_t(k) ) & 3266 * weight_substep(intermediate_timestep_count) 3267 ENDDO 3268 CASE ( 'qr' ) 3269 DO k = nzb, nzt 3270 sums_wsqrs_ws_l(k,tn) = sums_wsqrs_ws_l(k,tn) & 3271 + ( flux_t(k) + diss_t(k) ) & 3272 * weight_substep(intermediate_timestep_count) 3273 ENDDO 3274 CASE ( 'nr' ) 3275 DO k = nzb, nzt 3276 sums_wsnrs_ws_l(k,tn) = sums_wsnrs_ws_l(k,tn) & 3277 + ( flux_t(k) + diss_t(k) ) & 3278 * weight_substep(intermediate_timestep_count) 3279 ENDDO 3280 3281 END SELECT 3282 3283 ENDIF 3284 3285 ENDDO 3286 ENDDO 3287 3288 END SUBROUTINE advec_s_ws 3289 3290 3291 !------------------------------------------------------------------------------! 3292 ! Scalar advection - Call for all grid points - accelerator version 3293 !------------------------------------------------------------------------------! 3294 SUBROUTINE advec_s_ws_acc ( sk, sk_char ) 3295 3296 USE arrays_3d, & 3297 ONLY: ddzw, tend, u, v, w 3298 3299 USE constants, & 3300 ONLY: adv_sca_1, adv_sca_3, adv_sca_5 3301 3302 USE control_parameters, & 3303 ONLY: intermediate_timestep_count, monotonic_adjustment, u_gtrans,& 3304 v_gtrans 3305 3306 USE grid_variables, & 3307 ONLY: ddx, ddy 3308 3309 USE indices, & 3310 ONLY: i_left, i_right, j_north, j_south, nxlg, nxrg, nyng, nysg, & 3311 nzb, nzb_max, nzt, wall_flags_0 3312 3313 USE kinds 3314 3315 ! USE statistics, & 3316 ! ONLY: sums_wspts_ws_l, sums_wsqs_ws_l, sums_wssas_ws_l, & 3317 ! sums_wsqrs_ws_l, sums_wsnrs_ws_l, weight_substep 3318 3319 IMPLICIT NONE 3320 3321 CHARACTER (LEN = *), INTENT(IN) :: sk_char !: 3322 3323 INTEGER(iwp) :: i !: 3324 INTEGER(iwp) :: ibit0 !: 3325 INTEGER(iwp) :: ibit1 !: 3326 INTEGER(iwp) :: ibit2 !: 3327 INTEGER(iwp) :: ibit3 !: 3328 INTEGER(iwp) :: ibit4 !: 3329 INTEGER(iwp) :: ibit5 !: 3330 INTEGER(iwp) :: ibit6 !: 3331 INTEGER(iwp) :: ibit7 !: 3332 INTEGER(iwp) :: ibit8 !: 3333 INTEGER(iwp) :: j !: 3334 INTEGER(iwp) :: k !: 3335 INTEGER(iwp) :: k_mm !: 3336 INTEGER(iwp) :: k_mmm !: 3337 INTEGER(iwp) :: k_pp !: 3338 INTEGER(iwp) :: k_ppp !: 3339 INTEGER(iwp) :: tn = 0 !: 2735 3340 2736 3341 REAL(wp) :: diss_d !: … … 2747 3352 REAL(wp) :: flux_s !: 2748 3353 REAL(wp) :: flux_t !: 3354 REAL(wp) :: fd_1 !: 3355 REAL(wp) :: fl_1 !: 3356 REAL(wp) :: fn_1 !: 3357 REAL(wp) :: fr_1 !: 3358 REAL(wp) :: fs_1 !: 3359 REAL(wp) :: ft_1 !: 3360 REAL(wp) :: phi_d !: 3361 REAL(wp) :: phi_l !: 3362 REAL(wp) :: phi_n !: 3363 REAL(wp) :: phi_r !: 3364 REAL(wp) :: phi_s !: 3365 REAL(wp) :: phi_t !: 3366 REAL(wp) :: rd !: 3367 REAL(wp) :: rl !: 3368 REAL(wp) :: rn !: 3369 REAL(wp) :: rr !: 3370 REAL(wp) :: rs !: 3371 REAL(wp) :: rt !: 2749 3372 REAL(wp) :: u_comp !: 2750 3373 REAL(wp) :: v_comp !: … … 2967 3590 ( sk(k_ppp,j,i) - sk(k_mm,j,i) ) & 2968 3591 ) 3592 ! 3593 !-- Apply monotonic adjustment. 3594 IF ( monotonic_adjustment ) THEN 3595 ! 3596 !-- At first, calculate first order fluxes. 3597 u_comp = u(k,j,i) - u_gtrans 3598 fl_1 = ( u_comp * ( sk(k,j,i) + sk(k,j,i-1) ) & 3599 -ABS( u_comp ) * ( sk(k,j,i) - sk(k,j,i-1) ) & 3600 ) * adv_sca_1 3601 3602 u_comp = u(k,j,i+1) - u_gtrans 3603 fr_1 = ( u_comp * ( sk(k,j,i+1) + sk(k,j,i) ) & 3604 -ABS( u_comp ) * ( sk(k,j,i+1) - sk(k,j,i) ) & 3605 ) * adv_sca_1 3606 3607 v_comp = v(k,j,i) - v_gtrans 3608 fs_1 = ( v_comp * ( sk(k,j,i) + sk(k,j-1,i) ) & 3609 -ABS( v_comp ) * ( sk(k,j,i) - sk(k,j-1,i) ) & 3610 ) * adv_sca_1 3611 3612 v_comp = v(k,j+1,i) - v_gtrans 3613 fn_1 = ( v_comp * ( sk(k,j+1,i) + sk(k,j,i) ) & 3614 -ABS( v_comp ) * ( sk(k,j+1,i) - sk(k,j,i) ) & 3615 ) * adv_sca_1 3616 3617 fd_1 = ( w(k-1,j,i) * ( sk(k,j,i) + sk(k-1,j,i) ) & 3618 -ABS( w(k-1,j,i) ) * ( sk(k,j,i) - sk(k-1,j,i) ) & 3619 ) * adv_sca_1 3620 3621 ft_1 = ( w(k,j,i) * ( sk(k+1,j,i) + sk(k,j,i) ) & 3622 -ABS( w(k,j,i) ) * ( sk(k+1,j,i) - sk(k,j,i) ) & 3623 ) * adv_sca_1 3624 ! 3625 !-- Calculate ratio of upwind gradients. Note, Min/Max is just 3626 !-- to avoid if statements. 3627 rl = ( MAX( 0.0, u(k,j,i) - u_gtrans ) * & 3628 ABS( ( sk(k,j,i-1) - sk(k,j,i-2) ) /& 3629 ( sk(k,j,i) - sk(k,j,i-1) + 1E-20_wp ) & 3630 ) + & 3631 MIN( 0.0, u(k,j,i) - u_gtrans ) * & 3632 ABS( ( sk(k,j,i) - sk(k,j,i+1) ) /& 3633 ( sk(k,j,i-1) - sk(k,j,i) + 1E-20_wp ) & 3634 ) & 3635 ) / ABS( u(k,j,i) - u_gtrans + 1E-20_wp ) 3636 3637 rr = ( MAX( 0.0, u(k,j,i+1) - u_gtrans ) * & 3638 ABS( ( sk(k,j,i) - sk(k,j,i-1) ) /& 3639 ( sk(k,j,i+1) - sk(k,j,i) + 1E-20_wp ) & 3640 ) + & 3641 MIN( 0.0, u(k,j,i+1) - u_gtrans ) * & 3642 ABS( ( sk(k,j,i+1) - sk(k,j,i+2) ) /& 3643 ( sk(k,j,i) - sk(k,j,i+1) + 1E-20_wp ) & 3644 ) & 3645 ) / ABS( u(k,j,i+1) - u_gtrans + 1E-20_wp ) 3646 3647 rs = ( MAX( 0.0, v(k,j,i) - v_gtrans ) * & 3648 ABS( ( sk(k,j-1,i) - sk(k,j-2,i) ) /& 3649 ( sk(k,j,i) - sk(k,j-1,i) + 1E-20_wp ) & 3650 ) + & 3651 MIN( 0.0, v(k,j,i) - v_gtrans ) * & 3652 ABS( ( sk(k,j,i) - sk(k,j+1,i) ) /& 3653 ( sk(k,j-1,i) - sk(k,j,i) + 1E-20_wp ) & 3654 ) & 3655 ) / ABS( v(k,j,i) - v_gtrans + 1E-20_wp ) 3656 3657 rn = ( MAX( 0.0, v(k,j+1,i) - v_gtrans ) * & 3658 ABS( ( sk(k,j,i) - sk(k,j-1,i) ) /& 3659 ( sk(k,j+1,i) - sk(k,j,i) + 1E-20_wp ) & 3660 ) + & 3661 MIN( 0.0, v(k,j+1,i) - v_gtrans ) * & 3662 ABS( ( sk(k,j+1,i) - sk(k,j+2,i) ) /& 3663 ( sk(k,j,i) - sk(k,j+1,i) + 1E-20_wp ) & 3664 ) & 3665 ) / ABS( v(k,j+1,i) - v_gtrans + 1E-20_wp ) 3666 ! 3667 !-- Reuse k_mm and compute k_mmm for the vertical gradient ratios. 3668 !-- Note, for vertical advection below the third grid point above 3669 !-- surface ( or below the model top) rd and rt are set to 0, i.e. 3670 !-- use of first order scheme is enforced. 3671 k_mmm = k - 3 * ibit8 3672 3673 rd = ( MAX( 0.0, w(k-1,j,i) ) * & 3674 ABS( ( sk(k_mm,j,i) - sk(k_mmm,j,i) ) / & 3675 ( sk(k-1,j,i) - sk(k_mm,j,i) + 1E-20_wp ) & 3676 ) + & 3677 MIN( 0.0, w(k-1,j,i) ) * & 3678 ABS( ( sk(k-1,j,i) - sk(k,j,i) ) / & 3679 ( sk(k_mm,j,i) - sk(k-1,j,i) + 1E-20_wp ) & 3680 ) & 3681 ) * ibit8 / ABS( w(k-1,j,i) + 1E-20_wp ) 3682 3683 rt = ( MAX( 0.0, w(k,j,i) ) * & 3684 ABS( ( sk(k,j,i) - sk(k-1,j,i) ) / & 3685 ( sk(k+1,j,i) - sk(k,j,i) + 1E-20_wp ) & 3686 ) + & 3687 MIN( 0.0, w(k,j,i) ) * & 3688 ABS( ( sk(k+1,j,i) - sk(k_pp,j,i) ) / & 3689 ( sk(k,j,i) - sk(k+1,j,i) + 1E-20_wp ) & 3690 ) & 3691 ) * ibit8 / ABS( w(k,j,i) + 1E-20_wp ) 3692 ! 3693 !-- Calculate empirical limiter function (van Albada2 limiter). 3694 phi_l = MIN( 1.0_wp, ( 2.0_wp * rl ) / & 3695 ( rl**2 + 1.0_wp ) ) 3696 phi_r = MIN( 1.0_wp, ( 2.0_wp * rr ) / & 3697 ( rr**2 + 1.0_wp ) ) 3698 phi_s = MIN( 1.0_wp, ( 2.0_wp * rs ) / & 3699 ( rs**2 + 1.0_wp ) ) 3700 phi_n = MIN( 1.0_wp, ( 2.0_wp * rn ) / & 3701 ( rn**2 + 1.0_wp ) ) 3702 phi_d = MIN( 1.0_wp, ( 2.0_wp * rd ) / & 3703 ( rd**2 + 1.0_wp ) ) 3704 phi_t = MIN( 1.0_wp, ( 2.0_wp * rt ) / & 3705 ( rt**2 + 1.0_wp ) ) 3706 ! 3707 !-- Calculate the resulting monotone flux. 3708 flux_l = fl_1 - phi_l * ( fl_1 - flux_l ) 3709 flux_r = fr_1 - phi_r * ( fr_1 - flux_r ) 3710 flux_s = fs_1 - phi_s * ( fs_1 - flux_s ) 3711 flux_n = fn_1 - phi_n * ( fn_1 - flux_n ) 3712 flux_d = fd_1 - phi_d * ( fd_1 - flux_d ) 3713 flux_t = ft_1 - phi_t * ( ft_1 - flux_t ) 3714 ! 3715 !-- Moreover, modify dissipation flux according to the limiter. 3716 diss_l = diss_l * phi_l 3717 diss_r = diss_r * phi_r 3718 diss_s = diss_s * phi_s 3719 diss_n = diss_n * phi_n 3720 diss_d = diss_d * phi_d 3721 diss_t = diss_t * phi_t 3722 3723 ENDIF 2969 3724 ! 2970 3725 !-- Calculate the divergence of the velocity field. A respective -
palm/trunk/SOURCE/check_parameters.f90
r1556 r1557 20 20 ! Current revisions: 21 21 ! ----------------- 22 ! 22 ! Added checks for monotonic limiter 23 23 ! 24 24 ! Former revisions: … … 605 605 IF ( topography /= 'flat' ) THEN 606 606 action = ' ' 607 IF ( scalar_advec /= 'pw-scheme' .AND. scalar_advec /= 'ws-scheme') THEN 607 IF ( scalar_advec /= 'pw-scheme' .AND. scalar_advec /= 'ws-scheme' & 608 .AND. scalar_advec /= 'ws-scheme-mono' ) THEN 608 609 WRITE( action, '(A,A)' ) 'scalar_advec = ', scalar_advec 609 610 ENDIF … … 751 752 CALL message( 'check_parameters', 'PA0022', 1, 2, 0, 6, 0 ) 752 753 ENDIF 753 IF ( ( momentum_advec == 'ws-scheme' .OR. scalar_advec == 'ws-scheme' ) & 754 IF ( ( momentum_advec == 'ws-scheme' .OR. scalar_advec == 'ws-scheme' & 755 .OR. scalar_advec == 'ws-scheme-mono' ) & 754 756 .AND. ( timestep_scheme == 'euler' .OR. & 755 757 timestep_scheme == 'runge-kutta-2' ) ) & … … 761 763 ENDIF 762 764 IF ( scalar_advec /= 'pw-scheme' .AND. scalar_advec /= 'ws-scheme' .AND. & 763 scalar_advec /= ' bc-scheme' )&765 scalar_advec /= 'ws-scheme-mono' .AND. scalar_advec /= 'bc-scheme' ) & 764 766 THEN 765 767 message_string = 'unknown advection scheme: scalar_advec = "' // & … … 776 778 777 779 IF ( use_sgs_for_particles .AND. .NOT. use_upstream_for_tke .AND. & 778 scalar_advec /= 'ws-scheme' ) THEN 780 ( scalar_advec /= 'ws-scheme' .OR. & 781 scalar_advec /= 'ws-scheme-mono' ) & 782 ) THEN 779 783 use_upstream_for_tke = .TRUE. 780 784 message_string = 'use_upstream_for_tke set .TRUE. because ' // & … … 792 796 ! 793 797 !-- Set LOGICAL switches to enhance performance 794 IF ( momentum_advec == 'ws-scheme' ) ws_scheme_mom = .TRUE. 795 IF ( scalar_advec == 'ws-scheme' ) ws_scheme_sca = .TRUE. 798 IF ( momentum_advec == 'ws-scheme' ) ws_scheme_mom = .TRUE. 799 IF ( scalar_advec == 'ws-scheme' .OR. & 800 scalar_advec == 'ws-scheme-mono' ) ws_scheme_sca = .TRUE. 801 IF ( scalar_advec == 'ws-scheme-mono' ) monotonic_adjustment = .TRUE. 802 796 803 797 804 ! … … 1707 1714 CALL message( 'check_parameters', 'PA0051', 1, 2, 0, 6, 0 ) 1708 1715 ENDIF 1709 IF ( momentum_advec /= 'pw-scheme' .AND. & 1710 momentum_advec /= 'ws-scheme') THEN 1716 IF ( scalar_advec /= 'pw-scheme' .AND. & 1717 ( scalar_advec /= 'ws-scheme' .OR. & 1718 scalar_advec /= 'ws-scheme-mono' ) & 1719 ) THEN 1720 1711 1721 message_string = 'non-cyclic lateral boundaries do not allow ' // & 1712 1722 'momentum_advec = "' // TRIM( momentum_advec ) // '"' 1713 1723 CALL message( 'check_parameters', 'PA0052', 1, 2, 0, 6, 0 ) 1714 1724 ENDIF 1715 IF ( scalar_advec /= 'pw-scheme' .AND. & 1716 scalar_advec /= 'ws-scheme' ) THEN 1725 IF ( scalar_advec /= 'pw-scheme' .AND. & 1726 ( scalar_advec /= 'ws-scheme' .OR. & 1727 scalar_advec /= 'ws-scheme-mono' ) & 1728 ) THEN 1717 1729 message_string = 'non-cyclic lateral boundaries do not allow ' // & 1718 1730 'scalar_advec = "' // TRIM( scalar_advec ) // '"' -
palm/trunk/SOURCE/flow_statistics.f90
r1556 r1557 21 21 ! Current revisions: 22 22 ! ----------------- 23 ! 23 ! Adjustments for monotonic limiter 24 24 ! 25 25 ! Former revisions: … … 143 143 ONLY: average_count_pr, cloud_droplets, cloud_physics, do_sum, & 144 144 dt_3d, g, humidity, icloud_scheme, kappa, large_scale_forcing, & 145 large_scale_subsidence, max_pr_user, message_string, ocean, & 145 large_scale_subsidence, max_pr_user, message_string, & 146 monotonic_adjustment, ocean, & 146 147 passive_scalar, precipitation, simulated_time, & 147 148 use_subsidence_tendencies, use_surface_fluxes, use_top_fluxes, & … … 281 282 ENDIF 282 283 283 IF ( ws_scheme_sca .AND. sr == 0 ) THEN 284 IF ( ws_scheme_sca .AND. .NOT. monotonic_adjustment & 285 .AND. sr == 0 ) THEN 284 286 285 287 DO i = 0, threads_per_task-1 … … 802 804 !-- but so far there is no other suitable place to calculate) 803 805 IF ( ocean ) THEN 804 IF( .NOT. ws_scheme_sca .OR. sr /= 0 ) THEN 806 IF( .NOT. ws_scheme_sca .OR. monotonic_adjustment .OR. & 807 sr /= 0 ) THEN 805 808 pts = 0.5_wp * ( sa(k,j,i) - hom(k,1,23,sr) + & 806 809 sa(k+1,j,i) - hom(k+1,1,23,sr) ) … … 853 856 ENDIF 854 857 ELSE 855 IF( .NOT. ws_scheme_sca .OR. sr /= 0 ) THEN 858 IF( .NOT. ws_scheme_sca .OR. monotonic_adjustment .OR. & 859 sr /= 0 ) THEN 856 860 pts = 0.5_wp * ( vpt(k,j,i) - hom(k,1,44,sr) + & 857 861 vpt(k+1,j,i) - hom(k+1,1,44,sr) ) 858 862 sums_l(k,46,tn) = sums_l(k,46,tn) + pts * w(k,j,i) * & 859 863 rmask(j,i,sr) 860 ELSE IF ( ws_scheme_sca .AND. sr == 0 ) THEN 864 ELSE IF ( ws_scheme_sca .AND. .NOT. monotonic_adjustment & 865 .AND. sr == 0 ) THEN 861 866 sums_l(k,46,tn) = ( 1.0_wp + 0.61_wp * & 862 867 hom(k,1,41,sr) ) * & … … 870 875 !-- Passive scalar flux 871 876 IF ( passive_scalar .AND. ( .NOT. ws_scheme_sca & 872 .OR. sr /= 0 ) ) THEN877 .OR. monotonic_adjustment .OR. sr /= 0 ) ) THEN 873 878 pts = 0.5_wp * ( q(k,j,i) - hom(k,1,41,sr) + & 874 879 q(k+1,j,i) - hom(k+1,1,41,sr) ) … … 914 919 915 920 ENDIF 916 IF ( .NOT. ws_scheme_sca .OR. sr /= 0 ) THEN921 IF ( .NOT. ws_scheme_sca .OR. monotonic_adjustment .OR. sr /= 0 ) THEN 917 922 !$OMP DO 918 923 DO i = nxl, nxr … … 1525 1530 ONLY : average_count_pr, cloud_droplets, cloud_physics, do_sum, & 1526 1531 dt_3d, g, humidity, icloud_scheme, kappa, large_scale_forcing, & 1527 large_scale_subsidence, max_pr_user, message_string, ocean, & 1532 large_scale_subsidence, max_pr_user, message_string, & 1533 monotonic_adjustment, ocean, & 1528 1534 passive_scalar, precipitation, simulated_time, & 1529 1535 use_subsidence_tendencies, use_surface_fluxes, use_top_fluxes, & … … 1663 1669 ENDIF 1664 1670 1665 IF ( ws_scheme_sca .AND. sr == 0 ) THEN1671 IF ( ws_scheme_sca .AND. .NOT. monotonic_adjustment .AND. sr == 0 ) THEN 1666 1672 1667 1673 DO i = 0, threads_per_task-1 … … 2486 2492 IF ( ocean ) THEN 2487 2493 2488 IF( .NOT. ws_scheme_sca .OR. sr /= 0 ) THEN2494 IF( .NOT. ws_scheme_sca .OR. monotonic_adjustment .OR. sr /= 0 ) THEN 2489 2495 2490 2496 !$acc parallel loop gang present( hom, rflags_invers, rmask, sa, sums_l, w ) create( s1 ) … … 2628 2634 ELSE 2629 2635 2630 IF( .NOT. ws_scheme_sca .OR. sr /= 0 ) THEN2636 IF( .NOT. ws_scheme_sca .OR. monotonic_adjustment .OR. sr /= 0 ) THEN 2631 2637 2632 2638 !$acc parallel loop gang present( hom, rflags_invers, rmask, sums_l, vpt, w ) create( s1 ) … … 2660 2666 ! 2661 2667 !-- Passive scalar flux 2662 IF ( passive_scalar .AND. ( .NOT. ws_scheme_sca .OR. sr /= 0 ) ) THEN 2668 IF ( passive_scalar .AND. ( .NOT. ws_scheme_sca .OR. monotonic_adjustment & 2669 .OR. sr /= 0 ) ) THEN 2663 2670 2664 2671 !$acc parallel loop gang present( hom, q, rflags_invers, rmask, sums_l, w ) create( s1 ) … … 2714 2721 ENDIF 2715 2722 2716 IF ( .NOT. ws_scheme_sca .OR. sr /= 0 ) THEN2723 IF ( .NOT. ws_scheme_sca .OR. monotonic_adjustment .OR. sr /= 0 ) THEN 2717 2724 2718 2725 !$OMP DO -
palm/trunk/SOURCE/header.f90
r1552 r1557 20 20 ! Current revisions: 21 21 ! ----------------- 22 ! 22 ! output for monotonic limiter 23 23 ! 24 24 ! Former revisions: … … 452 452 ELSEIF ( scalar_advec == 'ws-scheme' ) THEN 453 453 WRITE ( io, 504 ) 454 ELSEIF ( scalar_advec == 'ws-scheme-mono' ) THEN 455 WRITE ( io, 513 ) 454 456 ELSE 455 457 WRITE ( io, 118 ) … … 2330 2332 ' Time: ',A8,6X,'Run-No.: ',I2.2/ & 2331 2333 ' Run on host: ',A10,6X,'En-No.: ',I2.2) 2334 513 FORMAT (' --> Scalar advection via Wicker-Skamarock-Scheme 5th order ' // & 2335 '+ monotonic adjustment') 2336 2332 2337 2333 2338 END SUBROUTINE header -
palm/trunk/SOURCE/init_grid.f90
r1419 r1557 20 20 ! Current revisions: 21 21 ! ----------------- 22 ! 22 ! Adjustment for monotoinic limiter 23 23 ! 24 24 ! Former revisions: … … 1169 1169 wall_flags_00 = 0 1170 1170 1171 IF ( scalar_advec == 'ws-scheme' ) THEN 1171 IF ( scalar_advec == 'ws-scheme' .OR. & 1172 scalar_advec == 'ws-scheme-mono' ) THEN 1172 1173 ! 1173 1174 !-- Set flags to steer the degradation of the advection scheme in advec_ws -
palm/trunk/SOURCE/init_pegrid.f90
r1469 r1557 20 20 ! Current revisions: 21 21 ! ------------------ 22 ! 22 ! Adjustment for monotonic limiter 23 23 ! 24 24 ! Former revisions: … … 591 591 ! 592 592 !-- Determine the number of ghost point layers 593 IF ( scalar_advec == 'ws-scheme' .OR. momentum_advec == 'ws-scheme' ) THEN 593 IF ( scalar_advec == 'ws-scheme' .OR. & 594 scalar_advec == 'ws-scheme-mono' .OR. & 595 momentum_advec == 'ws-scheme' ) THEN 594 596 nbgp = 3 595 597 ELSE -
palm/trunk/SOURCE/modules.f90
r1552 r1557 20 20 ! Current revisions: 21 21 ! ------------------ 22 ! 22 ! +monotonic_adjustment 23 23 ! 24 24 ! Former revisions: … … 659 659 lunudge = .FALSE., lvnudge = .FALSE., lwnudge = .FALSE., & 660 660 masking_method = .FALSE., mg_switch_to_pe0 = .FALSE., & 661 monotonic_adjustment = .FALSE., & 661 662 neutral = .FALSE., nudging = .FALSE., & 662 663 ocean = .FALSE., on_device = .FALSE., &
Note: See TracChangeset
for help on using the changeset viewer.
