Changeset 2232 for palm/trunk/SOURCE/diffusivities.f90
- Timestamp:
- May 30, 2017 5:47:52 PM (7 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/diffusivities.f90 (modified) (9 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/diffusivities.f90
r2119 r2232 20 20 ! Current revisions: 21 21 ! ----------------- 22 ! 22 ! Adjustments to new topography and surface concept 23 23 ! 24 24 ! Former revisions: … … 89 89 USE control_parameters, & 90 90 ONLY: atmos_ocean_sign, e_min, g, outflow_l, outflow_n, outflow_r, & 91 outflow_s, use_single_reference_value, wall_adjustment 91 outflow_s, use_single_reference_value, wall_adjustment, & 92 wall_adjustment_factor 92 93 93 94 USE indices, & 94 ONLY: nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb_s_inner, nzb, nzt 95 ONLY: nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, nzt, & 96 wall_flags_0 95 97 USE kinds 96 98 … … 100 102 ONLY : rmask, statistic_regions, sums_l_l 101 103 104 USE surface_mod, & 105 ONLY : bc_h 106 102 107 IMPLICIT NONE 103 108 … … 105 110 INTEGER(iwp) :: j !< 106 111 INTEGER(iwp) :: k !< 112 INTEGER(iwp) :: m !< 107 113 INTEGER(iwp) :: omp_get_thread_num !< 108 114 INTEGER(iwp) :: sr !< … … 110 116 111 117 REAL(wp) :: dvar_dz !< 118 REAL(wp) :: flag !< 112 119 REAL(wp) :: l !< 113 120 REAL(wp) :: ll !< … … 129 136 ! 130 137 !-- Compute the turbulent diffusion coefficient for momentum 131 !$OMP PARALLEL PRIVATE (dvar_dz,i,j,k,l,ll,l_stable,sqrt_e,sr,tn )138 !$OMP PARALLEL PRIVATE (dvar_dz,i,j,k,l,ll,l_stable,sqrt_e,sr,tn,flag) 132 139 !$ tn = omp_get_thread_num() 133 140 … … 138 145 DO i = nxlg, nxrg 139 146 DO j = nysg, nyng 140 DO k = 1, nzt 141 IF ( k > nzb_s_inner(j,i) ) THEN 142 e(k,j,i) = MAX( e(k,j,i), e_min ) 143 ENDIF 147 DO k = nzb+1, nzt 148 e(k,j,i) = MAX( e(k,j,i), e_min ) * & 149 MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) ) 144 150 ENDDO 145 151 ENDDO … … 150 156 DO i = nxlg, nxrg 151 157 DO j = nysg, nyng 152 DO k = 1, nzt 153 154 IF ( k > nzb_s_inner(j,i) ) THEN 155 156 sqrt_e = SQRT( e(k,j,i) ) 157 ! 158 !-- Determine the mixing length 159 dvar_dz = atmos_ocean_sign * & ! inverse effect of pt/rho_ocean gradient 160 ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k) 161 IF ( dvar_dz > 0.0_wp ) THEN 162 IF ( use_single_reference_value ) THEN 163 l_stable = 0.76_wp * sqrt_e / & 164 SQRT( g / var_reference * dvar_dz ) + 1E-5_wp 165 ELSE 166 l_stable = 0.76_wp * sqrt_e / & 167 SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5_wp 168 ENDIF 158 DO k = nzb+1, nzt 159 160 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) ) 161 162 sqrt_e = SQRT( e(k,j,i) ) 163 ! 164 !-- Determine the mixing length 165 dvar_dz = atmos_ocean_sign * & ! inverse effect of pt/rho_ocean gradient 166 ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k) 167 IF ( dvar_dz > 0.0_wp ) THEN 168 IF ( use_single_reference_value ) THEN 169 l_stable = 0.76_wp * sqrt_e / & 170 SQRT( g / var_reference * dvar_dz ) + 1E-5_wp 169 171 ELSE 170 l_stable = l_grid(k) 172 l_stable = 0.76_wp * sqrt_e / & 173 SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5_wp 171 174 ENDIF 172 ! 173 !-- Adjustment of the mixing length 174 IF ( wall_adjustment ) THEN 175 l = MIN( l_wall(k,j,i), l_grid(k), l_stable ) 176 ll = MIN( l_wall(k,j,i), l_grid(k) ) 177 ELSE 178 l = MIN( l_grid(k), l_stable ) 179 ll = l_grid(k) 180 ENDIF 181 182 ! 183 !-- Compute diffusion coefficients for momentum and heat 184 km(k,j,i) = 0.1_wp * l * sqrt_e 185 kh(k,j,i) = ( 1.0_wp + 2.0_wp * l / ll ) * km(k,j,i) 186 187 ! 188 !-- Summation for averaged profile (cf. flow_statistics) 189 DO sr = 0, statistic_regions 190 sums_l_l(k,sr,tn) = sums_l_l(k,sr,tn) + l * rmask(j,i,sr) 191 ENDDO 192 175 ELSE 176 l_stable = l_grid(k) 193 177 ENDIF 178 ! 179 !-- Adjustment of the mixing length 180 IF ( wall_adjustment ) THEN 181 l = MIN( wall_adjustment_factor * l_wall(k,j,i), l_grid(k), & 182 l_stable ) 183 ll = MIN( wall_adjustment_factor * l_wall(k,j,i), l_grid(k) ) 184 ELSE 185 l = MIN( l_grid(k), l_stable ) 186 ll = l_grid(k) 187 ENDIF 188 ! 189 !-- Compute diffusion coefficients for momentum and heat 190 km(k,j,i) = 0.1_wp * l * sqrt_e * flag 191 kh(k,j,i) = ( 1.0_wp + 2.0_wp * l / ll ) * km(k,j,i) * flag 192 193 194 ! 195 !-- Summation for averaged profile (cf. flow_statistics) 196 DO sr = 0, statistic_regions 197 sums_l_l(k,sr,tn) = sums_l_l(k,sr,tn) + l * rmask(j,i,sr) & 198 * flag 199 ENDDO 194 200 195 201 ENDDO … … 202 208 203 209 ! 204 !-- Set vertical boundary values (Neumann conditions both at bottom and top). 210 !-- Set vertical boundary values (Neumann conditions both at upward- and 211 !-- downward facing walls. To set wall-boundary values, the surface data type 212 !-- is applied. 205 213 !-- Horizontal boundary conditions at vertical walls are not set because 206 !-- so far vertical walls require usage of a Prandtl-layer where the boundary 207 !-- values of the diffusivities are not needed 214 !-- so far vertical surfaces require usage of a Prandtl-layer where the boundary 215 !-- values of the diffusivities are not needed. 216 !-- Upward-facing 217 !$OMP PARALLEL DO PRIVATE( i, j, k ) 218 DO m = 1, bc_h(0)%ns 219 i = bc_h(0)%i(m) 220 j = bc_h(0)%j(m) 221 k = bc_h(0)%k(m) 222 km(k-1,j,i) = km(k,j,i) 223 kh(k-1,j,i) = kh(k,j,i) 224 ENDDO 225 ! 226 !-- Downward facing surfaces 227 !$OMP PARALLEL DO PRIVATE( i, j, k ) 228 DO m = 1, bc_h(1)%ns 229 i = bc_h(1)%i(m) 230 j = bc_h(1)%j(m) 231 k = bc_h(1)%k(m) 232 km(k+1,j,i) = km(k,j,i) 233 kh(k+1,j,i) = kh(k,j,i) 234 ENDDO 235 ! 236 !-- Model top 208 237 !$OMP PARALLEL DO 209 238 DO i = nxlg, nxrg 210 239 DO j = nysg, nyng 211 km(nzb_s_inner(j,i),j,i) = km(nzb_s_inner(j,i)+1,j,i) 212 km(nzt+1,j,i) = km(nzt,j,i) 213 kh(nzb_s_inner(j,i),j,i) = kh(nzb_s_inner(j,i)+1,j,i) 214 kh(nzt+1,j,i) = kh(nzt,j,i) 240 km(nzt+1,j,i) = km(nzt,j,i) 241 kh(nzt+1,j,i) = kh(nzt,j,i) 215 242 ENDDO 216 243 ENDDO 244 217 245 218 246 !
Note: See TracChangeset
for help on using the changeset viewer.
