Changeset 106 for palm/trunk/SOURCE/prognostic_equations.f90
- Timestamp:
- Aug 16, 2007 2:30:26 PM (17 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/prognostic_equations.f90 (modified) (13 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/prognostic_equations.f90
r102 r106 4 4 ! Actual revisions: 5 5 ! ----------------- 6 ! +uswst, vswst as arguments in calls of diffusion_u|v 6 ! +uswst, vswst as arguments in calls of diffusion_u|v, 7 ! loops for u and v are starting from index nxlu, nysv, respectively (needed 8 ! for non-cyclic boundary conditions) 7 9 ! 8 10 ! Former revisions: … … 133 135 ! 134 136 !-- u-tendency terms with no communication 135 DO i = nxl , nxr137 DO i = nxlu, nxr 136 138 DO j = nys, nyn 137 139 ! … … 202 204 !-- v-tendency terms with no communication 203 205 DO i = nxl, nxr 204 DO j = nys , nyn206 DO j = nysv, nyn 205 207 ! 206 208 !-- Tendency terms … … 806 808 ! 807 809 !-- Tendency terms for u-velocity component 808 IF ( j < nyn+1) THEN810 IF ( .NOT. outflow_l .OR. i > nxl ) THEN 809 811 810 812 tend(:,j,i) = 0.0 … … 857 859 ! 858 860 !-- Tendency terms for v-velocity component 859 IF ( i < nxr+1) THEN861 IF ( .NOT. outflow_s .OR. j > nys ) THEN 860 862 861 863 tend(:,j,i) = 0.0 … … 906 908 ! 907 909 !-- Tendency terms for w-velocity component 908 IF ( i < nxr+1 .AND. j < nyn+1 ) THEN 909 910 tend(:,j,i) = 0.0 911 IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN 912 CALL advec_w_pw( i, j ) 913 ELSE 914 CALL advec_w_up( i, j ) 915 ENDIF 916 IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) & 917 THEN 918 CALL diffusion_w( i, j, ddzu, ddzw, km_m, km_damp_x, & 919 km_damp_y, tend, u_m, v_m, w_m ) 920 ELSE 921 CALL diffusion_w( i, j, ddzu, ddzw, km, km_damp_x, km_damp_y, & 922 tend, u, v, w ) 923 ENDIF 924 CALL coriolis( i, j, 3 ) 925 IF ( ocean ) THEN 926 CALL buoyancy( i, j, rho, prho_reference, 3, 64 ) 927 ELSE 928 IF ( .NOT. humidity ) THEN 929 CALL buoyancy( i, j, pt, pt_reference, 3, 4 ) 930 ELSE 931 CALL buoyancy( i, j, vpt, pt_reference, 3, 44 ) 932 ENDIF 933 ENDIF 934 CALL user_actions( i, j, 'w-tendency' ) 935 936 ! 937 !-- Prognostic equation for w-velocity component 938 DO k = nzb_w_inner(j,i)+1, nzt-1 939 w_p(k,j,i) = ( 1.0-tsc(1) ) * w_m(k,j,i) + tsc(1) * w(k,j,i) + & 940 dt_3d * ( & 941 tsc(2) * tend(k,j,i) + tsc(3) * tw_m(k,j,i) & 942 - tsc(4) * ( p(k+1,j,i) - p(k,j,i) ) * ddzu(k+1) & 943 ) - & 944 tsc(5) * rdf(k) * w(k,j,i) 945 ENDDO 946 947 ! 948 !-- Calculate tendencies for the next Runge-Kutta step 949 IF ( timestep_scheme(1:5) == 'runge' ) THEN 950 IF ( intermediate_timestep_count == 1 ) THEN 951 DO k = nzb_w_inner(j,i)+1, nzt-1 952 tw_m(k,j,i) = tend(k,j,i) 953 ENDDO 954 ELSEIF ( intermediate_timestep_count < & 955 intermediate_timestep_count_max ) THEN 956 DO k = nzb_w_inner(j,i)+1, nzt-1 957 tw_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tw_m(k,j,i) 958 ENDDO 959 ENDIF 960 ENDIF 961 962 ! 963 !-- Tendency terms for potential temperature 964 tend(:,j,i) = 0.0 965 IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN 966 CALL advec_s_pw( i, j, pt ) 967 ELSE 968 CALL advec_s_up( i, j, pt ) 969 ENDIF 970 IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) & 971 THEN 972 CALL diffusion_s( i, j, ddzu, ddzw, kh_m, pt_m, shf_m, & 973 tswst_m, tend ) 974 ELSE 975 CALL diffusion_s( i, j, ddzu, ddzw, kh, pt, shf, tswst, tend ) 976 ENDIF 977 978 ! 979 !-- If required compute heating/cooling due to long wave radiation 980 !-- processes 981 IF ( radiation ) THEN 982 CALL calc_radiation( i, j ) 983 ENDIF 984 985 ! 986 !-- If required compute impact of latent heat due to precipitation 987 IF ( precipitation ) THEN 988 CALL impact_of_latent_heat( i, j ) 989 ENDIF 990 CALL user_actions( i, j, 'pt-tendency' ) 991 992 ! 993 !-- Prognostic equation for potential temperature 994 DO k = nzb_s_inner(j,i)+1, nzt 995 pt_p(k,j,i) = ( 1.0-tsc(1) ) * pt_m(k,j,i) + tsc(1)*pt(k,j,i) +& 996 dt_3d * ( & 997 tsc(2) * tend(k,j,i) + tsc(3) * tpt_m(k,j,i) & 998 ) - & 999 tsc(5) * rdf(k) * ( pt(k,j,i) - pt_init(k) ) 1000 ENDDO 1001 1002 ! 1003 !-- Calculate tendencies for the next Runge-Kutta step 1004 IF ( timestep_scheme(1:5) == 'runge' ) THEN 1005 IF ( intermediate_timestep_count == 1 ) THEN 1006 DO k = nzb_s_inner(j,i)+1, nzt 1007 tpt_m(k,j,i) = tend(k,j,i) 1008 ENDDO 1009 ELSEIF ( intermediate_timestep_count < & 1010 intermediate_timestep_count_max ) THEN 1011 DO k = nzb_s_inner(j,i)+1, nzt 1012 tpt_m(k,j,i) = -9.5625 * tend(k,j,i) + & 1013 5.3125 * tpt_m(k,j,i) 1014 ENDDO 1015 ENDIF 1016 ENDIF 1017 1018 ! 1019 !-- If required, compute prognostic equation for salinity 1020 IF ( ocean ) THEN 1021 1022 ! 1023 !-- Tendency-terms for salinity 910 1024 tend(:,j,i) = 0.0 911 IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN 912 CALL advec_w_pw( i, j ) 913 ELSE 914 CALL advec_w_up( i, j ) 915 ENDIF 916 IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) & 1025 IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) & 917 1026 THEN 918 CALL diffusion_w( i, j, ddzu, ddzw, km_m, km_damp_x, & 919 km_damp_y, tend, u_m, v_m, w_m ) 920 ELSE 921 CALL diffusion_w( i, j, ddzu, ddzw, km, km_damp_x, km_damp_y, & 922 tend, u, v, w ) 923 ENDIF 924 CALL coriolis( i, j, 3 ) 925 IF ( ocean ) THEN 926 CALL buoyancy( i, j, rho, prho_reference, 3, 64 ) 927 ELSE 928 IF ( .NOT. humidity ) THEN 929 CALL buoyancy( i, j, pt, pt_reference, 3, 4 ) 930 ELSE 931 CALL buoyancy( i, j, vpt, pt_reference, 3, 44 ) 932 ENDIF 933 ENDIF 934 CALL user_actions( i, j, 'w-tendency' ) 935 936 ! 937 !-- Prognostic equation for w-velocity component 938 DO k = nzb_w_inner(j,i)+1, nzt-1 939 w_p(k,j,i) = ( 1.0-tsc(1) ) * w_m(k,j,i) + tsc(1) * w(k,j,i) + & 940 dt_3d * ( & 941 tsc(2) * tend(k,j,i) + tsc(3) * tw_m(k,j,i) & 942 - tsc(4) * ( p(k+1,j,i) - p(k,j,i) ) * ddzu(k+1) & 943 ) - & 944 tsc(5) * rdf(k) * w(k,j,i) 945 ENDDO 946 947 ! 948 !-- Calculate tendencies for the next Runge-Kutta step 949 IF ( timestep_scheme(1:5) == 'runge' ) THEN 950 IF ( intermediate_timestep_count == 1 ) THEN 951 DO k = nzb_w_inner(j,i)+1, nzt-1 952 tw_m(k,j,i) = tend(k,j,i) 953 ENDDO 954 ELSEIF ( intermediate_timestep_count < & 955 intermediate_timestep_count_max ) THEN 956 DO k = nzb_w_inner(j,i)+1, nzt-1 957 tw_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tw_m(k,j,i) 958 ENDDO 959 ENDIF 960 ENDIF 961 962 ! 963 !-- Tendency terms for potential temperature 964 tend(:,j,i) = 0.0 965 IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN 966 CALL advec_s_pw( i, j, pt ) 967 ELSE 968 CALL advec_s_up( i, j, pt ) 969 ENDIF 970 IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) & 971 THEN 972 CALL diffusion_s( i, j, ddzu, ddzw, kh_m, pt_m, shf_m, & 973 tswst_m, tend ) 974 ELSE 975 CALL diffusion_s( i, j, ddzu, ddzw, kh, pt, shf, tswst, tend ) 976 ENDIF 977 978 ! 979 !-- If required compute heating/cooling due to long wave radiation 980 !-- processes 981 IF ( radiation ) THEN 982 CALL calc_radiation( i, j ) 983 ENDIF 984 985 ! 986 !-- If required compute impact of latent heat due to precipitation 987 IF ( precipitation ) THEN 988 CALL impact_of_latent_heat( i, j ) 989 ENDIF 990 CALL user_actions( i, j, 'pt-tendency' ) 991 992 ! 993 !-- Prognostic equation for potential temperature 1027 CALL advec_s_pw( i, j, sa ) 1028 ELSE 1029 CALL advec_s_up( i, j, sa ) 1030 ENDIF 1031 CALL diffusion_s( i, j, ddzu, ddzw, kh, sa, saswsb, saswst, & 1032 tend ) 1033 1034 CALL user_actions( i, j, 'sa-tendency' ) 1035 1036 ! 1037 !-- Prognostic equation for salinity 994 1038 DO k = nzb_s_inner(j,i)+1, nzt 995 pt_p(k,j,i) = ( 1.0-tsc(1) ) * pt_m(k,j,i) + tsc(1)*pt(k,j,i) +& 996 dt_3d * ( & 997 tsc(2) * tend(k,j,i) + tsc(3) * tpt_m(k,j,i) & 998 ) - & 999 tsc(5) * rdf(k) * ( pt(k,j,i) - pt_init(k) ) 1039 sa_p(k,j,i) = tsc(1) * sa(k,j,i) + & 1040 dt_3d * ( & 1041 tsc(2) * tend(k,j,i) + tsc(3) * tsa_m(k,j,i) & 1042 ) - & 1043 tsc(5) * rdf(k) * ( sa(k,j,i) - sa_init(k) ) 1044 IF ( sa_p(k,j,i) < 0.0 ) sa_p(k,j,i) = 0.1 * sa(k,j,i) 1000 1045 ENDDO 1001 1046 … … 1005 1050 IF ( intermediate_timestep_count == 1 ) THEN 1006 1051 DO k = nzb_s_inner(j,i)+1, nzt 1007 t pt_m(k,j,i) = tend(k,j,i)1052 tsa_m(k,j,i) = tend(k,j,i) 1008 1053 ENDDO 1009 1054 ELSEIF ( intermediate_timestep_count < & 1010 1055 intermediate_timestep_count_max ) THEN 1011 1056 DO k = nzb_s_inner(j,i)+1, nzt 1012 tpt_m(k,j,i) = -9.5625 * tend(k,j,i) + & 1013 5.3125 * tpt_m(k,j,i) 1014 ENDDO 1015 ENDIF 1016 ENDIF 1017 1018 ! 1019 !-- If required, compute prognostic equation for salinity 1020 IF ( ocean ) THEN 1021 1022 ! 1023 !-- Tendency-terms for salinity 1024 tend(:,j,i) = 0.0 1025 IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) & 1026 THEN 1027 CALL advec_s_pw( i, j, sa ) 1057 tsa_m(k,j,i) = -9.5625 * tend(k,j,i) + & 1058 5.3125 * tsa_m(k,j,i) 1059 ENDDO 1060 ENDIF 1061 ENDIF 1062 1063 ! 1064 !-- Calculate density by the equation of state for seawater 1065 CALL eqn_state_seawater( i, j ) 1066 1067 ENDIF 1068 1069 ! 1070 !-- If required, compute prognostic equation for total water content / 1071 !-- scalar 1072 IF ( humidity .OR. passive_scalar ) THEN 1073 1074 ! 1075 !-- Tendency-terms for total water content / scalar 1076 tend(:,j,i) = 0.0 1077 IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) & 1078 THEN 1079 CALL advec_s_pw( i, j, q ) 1080 ELSE 1081 CALL advec_s_up( i, j, q ) 1082 ENDIF 1083 IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' )& 1084 THEN 1085 CALL diffusion_s( i, j, ddzu, ddzw, kh_m, q_m, qsws_m, & 1086 qswst_m, tend ) 1087 ELSE 1088 CALL diffusion_s( i, j, ddzu, ddzw, kh, q, qsws, qswst, & 1089 tend ) 1090 ENDIF 1091 1092 ! 1093 !-- If required compute decrease of total water content due to 1094 !-- precipitation 1095 IF ( precipitation ) THEN 1096 CALL calc_precipitation( i, j ) 1097 ENDIF 1098 CALL user_actions( i, j, 'q-tendency' ) 1099 1100 ! 1101 !-- Prognostic equation for total water content / scalar 1102 DO k = nzb_s_inner(j,i)+1, nzt 1103 q_p(k,j,i) = ( 1.0-tsc(1) ) * q_m(k,j,i) + tsc(1)*q(k,j,i) +& 1104 dt_3d * ( & 1105 tsc(2) * tend(k,j,i) + tsc(3) * tq_m(k,j,i) & 1106 ) - & 1107 tsc(5) * rdf(k) * ( q(k,j,i) - q_init(k) ) 1108 IF ( q_p(k,j,i) < 0.0 ) q_p(k,j,i) = 0.1 * q(k,j,i) 1109 ENDDO 1110 1111 ! 1112 !-- Calculate tendencies for the next Runge-Kutta step 1113 IF ( timestep_scheme(1:5) == 'runge' ) THEN 1114 IF ( intermediate_timestep_count == 1 ) THEN 1115 DO k = nzb_s_inner(j,i)+1, nzt 1116 tq_m(k,j,i) = tend(k,j,i) 1117 ENDDO 1118 ELSEIF ( intermediate_timestep_count < & 1119 intermediate_timestep_count_max ) THEN 1120 DO k = nzb_s_inner(j,i)+1, nzt 1121 tq_m(k,j,i) = -9.5625 * tend(k,j,i) + & 1122 5.3125 * tq_m(k,j,i) 1123 ENDDO 1124 ENDIF 1125 ENDIF 1126 1127 ENDIF 1128 1129 ! 1130 !-- If required, compute prognostic equation for turbulent kinetic 1131 !-- energy (TKE) 1132 IF ( .NOT. constant_diffusion ) THEN 1133 1134 ! 1135 !-- Tendency-terms for TKE 1136 tend(:,j,i) = 0.0 1137 IF ( ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) & 1138 .AND. .NOT. use_upstream_for_tke ) THEN 1139 CALL advec_s_pw( i, j, e ) 1140 ELSE 1141 CALL advec_s_up( i, j, e ) 1142 ENDIF 1143 IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' )& 1144 THEN 1145 IF ( .NOT. humidity ) THEN 1146 CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e_m, & 1147 km_m, l_grid, pt_m, pt_reference, & 1148 rif_m, tend, zu, zw ) 1028 1149 ELSE 1029 CALL advec_s_up( i, j, sa ) 1030 ENDIF 1031 CALL diffusion_s( i, j, ddzu, ddzw, kh, sa, saswsb, saswst, & 1032 tend ) 1033 1034 CALL user_actions( i, j, 'sa-tendency' ) 1035 1036 ! 1037 !-- Prognostic equation for salinity 1038 DO k = nzb_s_inner(j,i)+1, nzt 1039 sa_p(k,j,i) = tsc(1) * sa(k,j,i) + & 1040 dt_3d * ( & 1041 tsc(2) * tend(k,j,i) + tsc(3) * tsa_m(k,j,i) & 1042 ) - & 1043 tsc(5) * rdf(k) * ( sa(k,j,i) - sa_init(k) ) 1044 IF ( sa_p(k,j,i) < 0.0 ) sa_p(k,j,i) = 0.1 * sa(k,j,i) 1045 ENDDO 1046 1047 ! 1048 !-- Calculate tendencies for the next Runge-Kutta step 1049 IF ( timestep_scheme(1:5) == 'runge' ) THEN 1050 IF ( intermediate_timestep_count == 1 ) THEN 1051 DO k = nzb_s_inner(j,i)+1, nzt 1052 tsa_m(k,j,i) = tend(k,j,i) 1053 ENDDO 1054 ELSEIF ( intermediate_timestep_count < & 1055 intermediate_timestep_count_max ) THEN 1056 DO k = nzb_s_inner(j,i)+1, nzt 1057 tsa_m(k,j,i) = -9.5625 * tend(k,j,i) + & 1058 5.3125 * tsa_m(k,j,i) 1059 ENDDO 1060 ENDIF 1061 ENDIF 1062 1063 ! 1064 !-- Calculate density by the equation of state for seawater 1065 CALL eqn_state_seawater( i, j ) 1066 1067 ENDIF 1068 1069 ! 1070 !-- If required, compute prognostic equation for total water content / 1071 !-- scalar 1072 IF ( humidity .OR. passive_scalar ) THEN 1073 1074 ! 1075 !-- Tendency-terms for total water content / scalar 1076 tend(:,j,i) = 0.0 1077 IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) & 1078 THEN 1079 CALL advec_s_pw( i, j, q ) 1080 ELSE 1081 CALL advec_s_up( i, j, q ) 1082 ENDIF 1083 IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' )& 1084 THEN 1085 CALL diffusion_s( i, j, ddzu, ddzw, kh_m, q_m, qsws_m, & 1086 qswst_m, tend ) 1087 ELSE 1088 CALL diffusion_s( i, j, ddzu, ddzw, kh, q, qsws, qswst, & 1089 tend ) 1090 ENDIF 1091 1092 ! 1093 !-- If required compute decrease of total water content due to 1094 !-- precipitation 1095 IF ( precipitation ) THEN 1096 CALL calc_precipitation( i, j ) 1097 ENDIF 1098 CALL user_actions( i, j, 'q-tendency' ) 1099 1100 ! 1101 !-- Prognostic equation for total water content / scalar 1102 DO k = nzb_s_inner(j,i)+1, nzt 1103 q_p(k,j,i) = ( 1.0-tsc(1) ) * q_m(k,j,i) + tsc(1)*q(k,j,i) +& 1104 dt_3d * ( & 1105 tsc(2) * tend(k,j,i) + tsc(3) * tq_m(k,j,i) & 1106 ) - & 1107 tsc(5) * rdf(k) * ( q(k,j,i) - q_init(k) ) 1108 IF ( q_p(k,j,i) < 0.0 ) q_p(k,j,i) = 0.1 * q(k,j,i) 1109 ENDDO 1110 1111 ! 1112 !-- Calculate tendencies for the next Runge-Kutta step 1113 IF ( timestep_scheme(1:5) == 'runge' ) THEN 1114 IF ( intermediate_timestep_count == 1 ) THEN 1115 DO k = nzb_s_inner(j,i)+1, nzt 1116 tq_m(k,j,i) = tend(k,j,i) 1117 ENDDO 1118 ELSEIF ( intermediate_timestep_count < & 1119 intermediate_timestep_count_max ) THEN 1120 DO k = nzb_s_inner(j,i)+1, nzt 1121 tq_m(k,j,i) = -9.5625 * tend(k,j,i) + & 1122 5.3125 * tq_m(k,j,i) 1123 ENDDO 1124 ENDIF 1125 ENDIF 1126 1127 ENDIF 1128 1129 ! 1130 !-- If required, compute prognostic equation for turbulent kinetic 1131 !-- energy (TKE) 1132 IF ( .NOT. constant_diffusion ) THEN 1133 1134 ! 1135 !-- Tendency-terms for TKE 1136 tend(:,j,i) = 0.0 1137 IF ( ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) & 1138 .AND. .NOT. use_upstream_for_tke ) THEN 1139 CALL advec_s_pw( i, j, e ) 1140 ELSE 1141 CALL advec_s_up( i, j, e ) 1142 ENDIF 1143 IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' )& 1144 THEN 1145 IF ( .NOT. humidity ) THEN 1146 CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e_m, & 1147 km_m, l_grid, pt_m, pt_reference, & 1148 rif_m, tend, zu, zw ) 1150 CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e_m, & 1151 km_m, l_grid, vpt_m, pt_reference, & 1152 rif_m, tend, zu, zw ) 1153 ENDIF 1154 ELSE 1155 IF ( .NOT. humidity ) THEN 1156 IF ( ocean ) THEN 1157 CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, & 1158 km, l_grid, rho, prho_reference, & 1159 rif, tend, zu, zw ) 1149 1160 ELSE 1150 CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e _m,&1151 km _m, l_grid, vpt_m, pt_reference,&1152 rif_m,tend, zu, zw )1161 CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, & 1162 km, l_grid, pt, pt_reference, rif, & 1163 tend, zu, zw ) 1153 1164 ENDIF 1154 1165 ELSE 1155 IF ( .NOT. humidity ) THEN 1156 IF ( ocean ) THEN 1157 CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, & 1158 km, l_grid, rho, prho_reference, & 1159 rif, tend, zu, zw ) 1160 ELSE 1161 CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, & 1162 km, l_grid, pt, pt_reference, rif, & 1163 tend, zu, zw ) 1164 ENDIF 1165 ELSE 1166 CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, km, & 1167 l_grid, vpt, pt_reference, rif, tend, & 1168 zu, zw ) 1169 ENDIF 1170 ENDIF 1171 CALL production_e( i, j ) 1172 CALL user_actions( i, j, 'e-tendency' ) 1173 1174 ! 1175 !-- Prognostic equation for TKE. 1176 !-- Eliminate negative TKE values, which can occur due to numerical 1177 !-- reasons in the course of the integration. In such cases the old 1178 !-- TKE value is reduced by 90%. 1179 DO k = nzb_s_inner(j,i)+1, nzt 1180 e_p(k,j,i) = ( 1.0-tsc(1) ) * e_m(k,j,i) + tsc(1)*e(k,j,i) +& 1181 dt_3d * ( & 1182 tsc(2) * tend(k,j,i) + tsc(3) * te_m(k,j,i) & 1183 ) 1184 IF ( e_p(k,j,i) < 0.0 ) e_p(k,j,i) = 0.1 * e(k,j,i) 1185 ENDDO 1186 1187 ! 1188 !-- Calculate tendencies for the next Runge-Kutta step 1189 IF ( timestep_scheme(1:5) == 'runge' ) THEN 1190 IF ( intermediate_timestep_count == 1 ) THEN 1191 DO k = nzb_s_inner(j,i)+1, nzt 1192 te_m(k,j,i) = tend(k,j,i) 1193 ENDDO 1194 ELSEIF ( intermediate_timestep_count < & 1195 intermediate_timestep_count_max ) THEN 1196 DO k = nzb_s_inner(j,i)+1, nzt 1197 te_m(k,j,i) = -9.5625 * tend(k,j,i) + & 1198 5.3125 * te_m(k,j,i) 1199 ENDDO 1200 ENDIF 1201 ENDIF 1202 1203 ENDIF ! TKE equation 1204 1205 ENDIF ! Gridpoints excluding the non-cyclic wall 1166 CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, km, & 1167 l_grid, vpt, pt_reference, rif, tend, & 1168 zu, zw ) 1169 ENDIF 1170 ENDIF 1171 CALL production_e( i, j ) 1172 CALL user_actions( i, j, 'e-tendency' ) 1173 1174 ! 1175 !-- Prognostic equation for TKE. 1176 !-- Eliminate negative TKE values, which can occur due to numerical 1177 !-- reasons in the course of the integration. In such cases the old 1178 !-- TKE value is reduced by 90%. 1179 DO k = nzb_s_inner(j,i)+1, nzt 1180 e_p(k,j,i) = ( 1.0-tsc(1) ) * e_m(k,j,i) + tsc(1)*e(k,j,i) +& 1181 dt_3d * ( & 1182 tsc(2) * tend(k,j,i) + tsc(3) * te_m(k,j,i) & 1183 ) 1184 IF ( e_p(k,j,i) < 0.0 ) e_p(k,j,i) = 0.1 * e(k,j,i) 1185 ENDDO 1186 1187 ! 1188 !-- Calculate tendencies for the next Runge-Kutta step 1189 IF ( timestep_scheme(1:5) == 'runge' ) THEN 1190 IF ( intermediate_timestep_count == 1 ) THEN 1191 DO k = nzb_s_inner(j,i)+1, nzt 1192 te_m(k,j,i) = tend(k,j,i) 1193 ENDDO 1194 ELSEIF ( intermediate_timestep_count < & 1195 intermediate_timestep_count_max ) THEN 1196 DO k = nzb_s_inner(j,i)+1, nzt 1197 te_m(k,j,i) = -9.5625 * tend(k,j,i) + & 1198 5.3125 * te_m(k,j,i) 1199 ENDDO 1200 ENDIF 1201 ENDIF 1202 1203 ENDIF ! TKE equation 1206 1204 1207 1205 ENDDO … … 1268 1266 ! 1269 1267 !-- Prognostic equation for u-velocity component 1270 DO i = nxl , nxr1268 DO i = nxlu, nxr 1271 1269 DO j = nys, nyn 1272 1270 DO k = nzb_u_inner(j,i)+1, nzt … … 1285 1283 IF ( timestep_scheme(1:5) == 'runge' ) THEN 1286 1284 IF ( intermediate_timestep_count == 1 ) THEN 1287 DO i = nxl , nxr1285 DO i = nxlu, nxr 1288 1286 DO j = nys, nyn 1289 1287 DO k = nzb_u_inner(j,i)+1, nzt … … 1294 1292 ELSEIF ( intermediate_timestep_count < & 1295 1293 intermediate_timestep_count_max ) THEN 1296 DO i = nxl , nxr1294 DO i = nxlu, nxr 1297 1295 DO j = nys, nyn 1298 1296 DO k = nzb_u_inner(j,i)+1, nzt … … 1340 1338 !-- Prognostic equation for v-velocity component 1341 1339 DO i = nxl, nxr 1342 DO j = nys , nyn1340 DO j = nysv, nyn 1343 1341 DO k = nzb_v_inner(j,i)+1, nzt 1344 1342 v_p(k,j,i) = ( 1.0-tsc(1) ) * v_m(k,j,i) + tsc(1) * v(k,j,i) + & … … 1357 1355 IF ( intermediate_timestep_count == 1 ) THEN 1358 1356 DO i = nxl, nxr 1359 DO j = nys , nyn1357 DO j = nysv, nyn 1360 1358 DO k = nzb_v_inner(j,i)+1, nzt 1361 1359 tv_m(k,j,i) = tend(k,j,i) … … 1366 1364 intermediate_timestep_count_max ) THEN 1367 1365 DO i = nxl, nxr 1368 DO j = nys , nyn1366 DO j = nysv, nyn 1369 1367 DO k = nzb_v_inner(j,i)+1, nzt 1370 1368 tv_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tv_m(k,j,i)
Note: See TracChangeset
for help on using the changeset viewer.
