Changeset 2232 for palm/trunk/SOURCE/production_e.f90
- Timestamp:
- May 30, 2017 5:47:52 PM (7 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/production_e.f90 (modified) (19 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/production_e.f90
r2127 r2232 20 20 ! Current revisions: 21 21 ! ----------------- 22 ! 22 ! Adjustments to new surface concept 23 23 ! 24 24 ! Former revisions: … … 106 106 !------------------------------------------------------------------------------! 107 107 MODULE production_e_mod 108 109 110 USE wall_fluxes_mod, &111 ONLY: wall_fluxes_e112 108 113 109 USE kinds … … 115 111 PRIVATE 116 112 PUBLIC production_e, production_e_init 117 118 LOGICAL, SAVE :: first_call = .TRUE. !<119 120 REAL(wp), DIMENSION(:,:), ALLOCATABLE, SAVE :: u_0 !<121 REAL(wp), DIMENSION(:,:), ALLOCATABLE, SAVE :: v_0 !<122 113 123 114 INTERFACE production_e … … 141 132 142 133 USE arrays_3d, & 143 ONLY: ddzw, dd2zu, kh, km, prho, pt, q, ql, qsws, qswst, shf, & 144 tend, tswst, u, v, vpt, w 134 ONLY: ddzw, dd2zu, kh, km, prho, pt, q, ql, tend, u, v, vpt, w 145 135 146 136 USE cloud_parameters, & … … 154 144 155 145 USE grid_variables, & 156 ONLY: ddx, dx, ddy, dy , wall_e_x, wall_e_y146 ONLY: ddx, dx, ddy, dy 157 147 158 148 USE indices, & 159 ONLY: nxl, nxr, nys, nyn, nzb, nzb_diff_s_inner, & 160 nzb_diff_s_outer, nzb_s_inner, nzt, nzt_diff 149 ONLY: nxl, nxr, nys, nyn, nzb, nzt, wall_flags_0 150 151 USE surface_mod, & 152 ONLY : surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, & 153 surf_usm_v 161 154 162 155 IMPLICIT NONE 163 156 164 INTEGER(iwp) :: i !< 165 INTEGER(iwp) :: j !< 166 INTEGER(iwp) :: k !< 167 168 REAL(wp) :: def !< 169 REAL(wp) :: dudx !< 170 REAL(wp) :: dudy !< 171 REAL(wp) :: dudz !< 172 REAL(wp) :: dvdx !< 173 REAL(wp) :: dvdy !< 174 REAL(wp) :: dvdz !< 175 REAL(wp) :: dwdx !< 176 REAL(wp) :: dwdy !< 177 REAL(wp) :: dwdz !< 157 INTEGER(iwp) :: i !< running index x-direction 158 INTEGER(iwp) :: j !< running index y-direction 159 INTEGER(iwp) :: k !< running index z-direction 160 INTEGER(iwp) :: l !< running index for different surface type orientation 161 INTEGER(iwp) :: m !< running index surface elements 162 INTEGER(iwp) :: surf_e !< end index of surface elements at given i-j position 163 INTEGER(iwp) :: surf_s !< start index of surface elements at given i-j position 164 165 REAL(wp) :: def !< 166 REAL(wp) :: flag !< flag to mask topography 178 167 REAL(wp) :: k1 !< 179 168 REAL(wp) :: k2 !< 180 REAL(wp) :: km_neutral !< 169 REAL(wp) :: km_neutral !< diffusion coefficient assuming neutral conditions - used to compute shear production at surfaces 181 170 REAL(wp) :: theta !< 182 171 REAL(wp) :: temp !< 183 184 ! REAL(wp), DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) :: usvs, vsus, wsus, wsvs 185 REAL(wp), DIMENSION(nzb:nzt+1) :: usvs !< 186 REAL(wp), DIMENSION(nzb:nzt+1) :: vsus !< 187 REAL(wp), DIMENSION(nzb:nzt+1) :: wsus !< 188 REAL(wp), DIMENSION(nzb:nzt+1) :: wsvs !< 189 190 ! 191 !-- First calculate horizontal momentum flux u'v', w'v', v'u', w'u' at 192 !-- vertical walls, if neccessary 193 !-- So far, results are slightly different from the ij-Version. 194 !-- Therefore, ij-Version is called further below within the ij-loops. 195 ! IF ( topography /= 'flat' ) THEN 196 ! CALL wall_fluxes_e( usvs, 1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, wall_e_y ) 197 ! CALL wall_fluxes_e( wsvs, 0.0_wp, 0.0_wp, 1.0_wp, 0.0_wp, wall_e_y ) 198 ! CALL wall_fluxes_e( vsus, 0.0_wp, 1.0_wp, 0.0_wp, 0.0_wp, wall_e_x ) 199 ! CALL wall_fluxes_e( wsus, 0.0_wp, 0.0_wp, 0.0_wp, 1.0_wp, wall_e_x ) 200 ! ENDIF 201 172 REAL(wp) :: sign_dir !< sign of wall-tke flux, depending on wall orientation 173 REAL(wp) :: usvs !< momentum flux u"v" 174 REAL(wp) :: vsus !< momentum flux v"u" 175 REAL(wp) :: wsus !< momentum flux w"u" 176 REAL(wp) :: wsvs !< momentum flux w"v" 177 178 REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) :: dudx !< Gradient of u-component in x-direction 179 REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) :: dudy !< Gradient of u-component in y-direction 180 REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) :: dudz !< Gradient of u-component in z-direction 181 REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) :: dvdx !< Gradient of v-component in x-direction 182 REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) :: dvdy !< Gradient of v-component in y-direction 183 REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) :: dvdz !< Gradient of v-component in z-direction 184 REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) :: dwdx !< Gradient of w-component in x-direction 185 REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) :: dwdy !< Gradient of w-component in y-direction 186 REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) :: dwdz !< Gradient of w-component in z-direction 202 187 203 188 DO i = nxl, nxr 204 189 205 ! 206 !-- Calculate TKE production by shear 207 DO j = nys, nyn 208 DO k = nzb_diff_s_outer(j,i), nzt 209 210 dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx 211 dudy = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - & 212 u(k,j-1,i) - u(k,j-1,i+1) ) * ddy 213 dudz = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - & 214 u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k) 215 216 dvdx = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - & 217 v(k,j,i-1) - v(k,j+1,i-1) ) * ddx 218 dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy 219 dvdz = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - & 220 v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k) 221 222 dwdx = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - & 223 w(k,j,i-1) - w(k-1,j,i-1) ) * ddx 224 dwdy = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - & 225 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy 226 dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) 227 228 def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 ) + & 229 dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + & 230 dvdz**2 + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz ) 231 232 IF ( def < 0.0_wp ) def = 0.0_wp 233 234 tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def 235 190 IF ( constant_flux_layer ) THEN 191 192 ! 193 !-- Calculate TKE production by shear. Calculate gradients at all grid 194 !-- points first, gradients at surface-bounded grid points will be 195 !-- overwritten further below. 196 DO j = nys, nyn 197 DO k = nzb+1, nzt 198 199 dudx(k,j) = ( u(k,j,i+1) - u(k,j,i) ) * ddx 200 dudy(k,j) = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - & 201 u(k,j-1,i) - u(k,j-1,i+1) ) * ddy 202 dudz(k,j) = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - & 203 u(k-1,j,i) - u(k-1,j,i+1) ) * & 204 dd2zu(k) 205 206 dvdx(k,j) = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - & 207 v(k,j,i-1) - v(k,j+1,i-1) ) * ddx 208 dvdy(k,j) = ( v(k,j+1,i) - v(k,j,i) ) * ddy 209 dvdz(k,j) = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - & 210 v(k-1,j,i) - v(k-1,j+1,i) ) * & 211 dd2zu(k) 212 213 dwdx(k,j) = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - & 214 w(k,j,i-1) - w(k-1,j,i-1) ) * ddx 215 dwdy(k,j) = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - & 216 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy 217 dwdz(k,j) = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) 218 219 ENDDO 236 220 ENDDO 237 ENDDO238 239 IF ( constant_flux_layer ) THEN240 221 241 222 ! … … 244 225 !-- 'bottom and wall: use u_0,v_0 and wall functions' 245 226 DO j = nys, nyn 246 247 IF ( ( wall_e_x(j,i) /= 0.0_wp ) .OR. ( wall_e_y(j,i) /= 0.0_wp ) ) & 248 THEN 249 250 k = nzb_diff_s_inner(j,i) - 1 251 dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx 252 dudz = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - & 253 u_0(j,i) - u_0(j,i+1) ) * dd2zu(k) 254 dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy 255 dvdz = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - & 256 v_0(j,i) - v_0(j+1,i) ) * dd2zu(k) 257 dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) 258 259 IF ( wall_e_y(j,i) /= 0.0_wp ) THEN 260 ! 261 !-- Inconsistency removed: as the thermal stratification is 262 !-- not taken into account for the evaluation of the wall 263 !-- fluxes at vertical walls, the eddy viscosity km must not 264 !-- be used for the evaluation of the velocity gradients dudy 265 !-- and dwdy 266 !-- Note: The validity of the new method has not yet been 267 !-- shown, as so far no suitable data for a validation 268 !-- has been available 269 CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, & 270 usvs, 1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp ) 271 CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, & 272 wsvs, 0.0_wp, 0.0_wp, 1.0_wp, 0.0_wp ) 273 km_neutral = kappa * ( usvs(k)**2 + wsvs(k)**2 )**0.25_wp * & 274 0.5_wp * dy 275 IF ( km_neutral > 0.0_wp ) THEN 276 dudy = - wall_e_y(j,i) * usvs(k) / km_neutral 277 dwdy = - wall_e_y(j,i) * wsvs(k) / km_neutral 278 ELSE 279 dudy = 0.0_wp 280 dwdy = 0.0_wp 281 ENDIF 282 ELSE 283 dudy = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - & 284 u(k,j-1,i) - u(k,j-1,i+1) ) * ddy 285 dwdy = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - & 286 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy 287 ENDIF 288 289 IF ( wall_e_x(j,i) /= 0.0_wp ) THEN 290 ! 291 !-- Inconsistency removed: as the thermal stratification is 292 !-- not taken into account for the evaluation of the wall 293 !-- fluxes at vertical walls, the eddy viscosity km must not 294 !-- be used for the evaluation of the velocity gradients dvdx 295 !-- and dwdx 296 !-- Note: The validity of the new method has not yet been 297 !-- shown, as so far no suitable data for a validation 298 !-- has been available 299 CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, & 300 vsus, 0.0_wp, 1.0_wp, 0.0_wp, 0.0_wp ) 301 CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, & 302 wsus, 0.0_wp, 0.0_wp, 0.0_wp, 1.0_wp ) 303 km_neutral = kappa * ( vsus(k)**2 + wsus(k)**2 )**0.25_wp * & 304 0.5_wp * dx 305 IF ( km_neutral > 0.0_wp ) THEN 306 dvdx = - wall_e_x(j,i) * vsus(k) / km_neutral 307 dwdx = - wall_e_x(j,i) * wsus(k) / km_neutral 308 ELSE 309 dvdx = 0.0_wp 310 dwdx = 0.0_wp 311 ENDIF 312 ELSE 313 dvdx = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - & 314 v(k,j,i-1) - v(k,j+1,i-1) ) * ddx 315 dwdx = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - & 316 w(k,j,i-1) - w(k-1,j,i-1) ) * ddx 317 ENDIF 318 319 def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 ) + & 320 dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + & 321 dvdz**2 + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz ) 227 ! 228 !-- Compute gradients at north- and south-facing surfaces. 229 !-- First, for default surfaces, then for urban surfaces. 230 !-- Note, so far no natural vertical surfaces implemented 231 DO l = 0, 1 232 surf_s = surf_def_v(l)%start_index(j,i) 233 surf_e = surf_def_v(l)%end_index(j,i) 234 DO m = surf_s, surf_e 235 k = surf_def_v(l)%k(m) 236 usvs = surf_def_v(l)%mom_flux_tke(0,m) 237 wsvs = surf_def_v(l)%mom_flux_tke(1,m) 238 239 km_neutral = kappa * ( usvs**2 + wsvs**2 )**0.25_wp & 240 * 0.5_wp * dy 241 ! 242 !-- -1.0 for right-facing wall, 1.0 for left-facing wall 243 sign_dir = MERGE( 1.0_wp, -1.0_wp, & 244 BTEST( wall_flags_0(k,j-1,i), 0 ) ) 245 dudy(k,j) = sign_dir * usvs / ( km_neutral + 1E-10_wp ) 246 dwdy(k,j) = sign_dir * wsvs / ( km_neutral + 1E-10_wp ) 247 ENDDO 248 ! 249 !-- Natural surfaces 250 surf_s = surf_lsm_v(l)%start_index(j,i) 251 surf_e = surf_lsm_v(l)%end_index(j,i) 252 DO m = surf_s, surf_e 253 k = surf_lsm_v(l)%k(m) 254 usvs = surf_lsm_v(l)%mom_flux_tke(0,m) 255 wsvs = surf_lsm_v(l)%mom_flux_tke(1,m) 256 257 km_neutral = kappa * ( usvs**2 + wsvs**2 )**0.25_wp & 258 * 0.5_wp * dy 259 ! 260 !-- -1.0 for right-facing wall, 1.0 for left-facing wall 261 sign_dir = MERGE( 1.0_wp, -1.0_wp, & 262 BTEST( wall_flags_0(k,j-1,i), 0 ) ) 263 dudy(k,j) = sign_dir * usvs / ( km_neutral + 1E-10_wp ) 264 dwdy(k,j) = sign_dir * wsvs / ( km_neutral + 1E-10_wp ) 265 ENDDO 266 ! 267 !-- Urban surfaces 268 surf_s = surf_usm_v(l)%start_index(j,i) 269 surf_e = surf_usm_v(l)%end_index(j,i) 270 DO m = surf_s, surf_e 271 k = surf_usm_v(l)%k(m) 272 usvs = surf_usm_v(l)%mom_flux_tke(0,m) 273 wsvs = surf_usm_v(l)%mom_flux_tke(1,m) 274 275 km_neutral = kappa * ( usvs**2 + wsvs**2 )**0.25_wp & 276 * 0.5_wp * dy 277 ! 278 !-- -1.0 for right-facing wall, 1.0 for left-facing wall 279 sign_dir = MERGE( 1.0_wp, -1.0_wp, & 280 BTEST( wall_flags_0(k,j-1,i), 0 ) ) 281 dudy(k,j) = sign_dir * usvs / ( km_neutral + 1E-10_wp ) 282 dwdy(k,j) = sign_dir * wsvs / ( km_neutral + 1E-10_wp ) 283 ENDDO 284 ENDDO 285 ! 286 !-- Compute gradients at east- and west-facing walls 287 DO l = 2, 3 288 surf_s = surf_def_v(l)%start_index(j,i) 289 surf_e = surf_def_v(l)%end_index(j,i) 290 DO m = surf_s, surf_e 291 k = surf_def_v(l)%k(m) 292 vsus = surf_def_v(l)%mom_flux_tke(0,m) 293 wsus = surf_def_v(l)%mom_flux_tke(1,m) 294 295 km_neutral = kappa * ( vsus**2 + wsus**2 )**0.25_wp & 296 * 0.5_wp * dx 297 ! 298 !-- -1.0 for right-facing wall, 1.0 for left-facing wall 299 sign_dir = MERGE( 1.0_wp, -1.0_wp, & 300 BTEST( wall_flags_0(k,j,i-1), 0 ) ) 301 dvdx(k,j) = sign_dir * vsus / ( km_neutral + 1E-10_wp ) 302 dwdx(k,j) = sign_dir * wsus / ( km_neutral + 1E-10_wp ) 303 ENDDO 304 ! 305 !-- Natural surfaces 306 surf_s = surf_lsm_v(l)%start_index(j,i) 307 surf_e = surf_lsm_v(l)%end_index(j,i) 308 DO m = surf_s, surf_e 309 k = surf_lsm_v(l)%k(m) 310 vsus = surf_lsm_v(l)%mom_flux_tke(0,m) 311 wsus = surf_lsm_v(l)%mom_flux_tke(1,m) 312 313 km_neutral = kappa * ( vsus**2 + wsus**2 )**0.25_wp & 314 * 0.5_wp * dx 315 ! 316 !-- -1.0 for right-facing wall, 1.0 for left-facing wall 317 sign_dir = MERGE( 1.0_wp, -1.0_wp, & 318 BTEST( wall_flags_0(k,j,i-1), 0 ) ) 319 dvdx(k,j) = sign_dir * vsus / ( km_neutral + 1E-10_wp ) 320 dwdx(k,j) = sign_dir * wsus / ( km_neutral + 1E-10_wp ) 321 ENDDO 322 ! 323 !-- Urban surfaces 324 surf_s = surf_usm_v(l)%start_index(j,i) 325 surf_e = surf_usm_v(l)%end_index(j,i) 326 DO m = surf_s, surf_e 327 k = surf_usm_v(l)%k(m) 328 vsus = surf_usm_v(l)%mom_flux_tke(0,m) 329 wsus = surf_usm_v(l)%mom_flux_tke(1,m) 330 331 km_neutral = kappa * ( vsus**2 + wsus**2 )**0.25_wp & 332 * 0.5_wp * dx 333 ! 334 !-- -1.0 for right-facing wall, 1.0 for left-facing wall 335 sign_dir = MERGE( 1.0_wp, -1.0_wp, & 336 BTEST( wall_flags_0(k,j,i-1), 0 ) ) 337 dvdx(k,j) = sign_dir * vsus / ( km_neutral + 1E-10_wp ) 338 dwdx(k,j) = sign_dir * wsus / ( km_neutral + 1E-10_wp ) 339 ENDDO 340 ENDDO 341 ! 342 !-- Compute gradients at upward-facing surfaces 343 surf_s = surf_def_h(0)%start_index(j,i) 344 surf_e = surf_def_h(0)%end_index(j,i) 345 DO m = surf_s, surf_e 346 k = surf_def_h(0)%k(m) 347 ! 348 !-- Please note, actually, an interpolation of u_0 and v_0 349 !-- onto the grid center would be required. However, this 350 !-- would require several data transfers between 2D-grid and 351 !-- wall type. The effect of this missing interpolation is 352 !-- negligible. (See also production_e_init). 353 dudz(k,j) = ( u(k+1,j,i) - surf_def_h(0)%u_0(m) ) * dd2zu(k) 354 dvdz(k,j) = ( v(k+1,j,i) - surf_def_h(0)%v_0(m) ) * dd2zu(k) 355 356 ENDDO 357 ! 358 !-- Natural surfaces 359 surf_s = surf_lsm_h%start_index(j,i) 360 surf_e = surf_lsm_h%end_index(j,i) 361 DO m = surf_s, surf_e 362 k = surf_lsm_h%k(m) 363 ! 364 !-- Please note, actually, an interpolation of u_0 and v_0 365 !-- onto the grid center would be required. However, this 366 !-- would require several data transfers between 2D-grid and 367 !-- wall type. The effect of this missing interpolation is 368 !-- negligible. (See also production_e_init). 369 dudz(k,j) = ( u(k+1,j,i) - surf_lsm_h%u_0(m) ) * dd2zu(k) 370 dvdz(k,j) = ( v(k+1,j,i) - surf_lsm_h%v_0(m) ) * dd2zu(k) 371 372 ENDDO 373 ! 374 !-- Urban surfaces 375 surf_s = surf_usm_h%start_index(j,i) 376 surf_e = surf_usm_h%end_index(j,i) 377 DO m = surf_s, surf_e 378 k = surf_usm_h%k(m) 379 ! 380 !-- Please note, actually, an interpolation of u_0 and v_0 381 !-- onto the grid center would be required. However, this 382 !-- would require several data transfers between 2D-grid and 383 !-- wall type. The effect of this missing interpolation is 384 !-- negligible. (See also production_e_init). 385 dudz(k,j) = ( u(k+1,j,i) - surf_usm_h%u_0(m) ) * dd2zu(k) 386 dvdz(k,j) = ( v(k+1,j,i) - surf_usm_h%v_0(m) ) * dd2zu(k) 387 388 ENDDO 389 ! 390 !-- Compute gradients at downward-facing walls, only for 391 !-- non-natural default surfaces 392 surf_s = surf_def_h(1)%start_index(j,i) 393 surf_e = surf_def_h(1)%end_index(j,i) 394 DO m = surf_s, surf_e 395 k = surf_def_h(1)%k(m) 396 ! 397 !-- Please note, actually, an interpolation of u_0 and v_0 398 !-- onto the grid center would be required. However, this 399 !-- would require several data transfers between 2D-grid and 400 !-- wall type. The effect of this missing interpolation is 401 !-- negligible. (See also production_e_init). 402 dudz(k,j) = ( surf_def_h(1)%u_0(m) - u(k-1,j,i) ) * dd2zu(k) 403 dvdz(k,j) = ( surf_def_h(1)%v_0(m) - v(k-1,j,i) ) * dd2zu(k) 404 405 ENDDO 406 407 ENDDO 408 409 DO j = nys, nyn 410 DO k = nzb+1, nzt 411 412 def = 2.0_wp * ( dudx(k,j)**2 + dvdy(k,j)**2 + dwdz(k,j)**2 ) + & 413 dudy(k,j)**2 + dvdx(k,j)**2 + dwdx(k,j)**2 + & 414 dwdy(k,j)**2 + dudz(k,j)**2 + dvdz(k,j)**2 + & 415 2.0_wp * ( dvdx(k,j)*dudy(k,j) + dwdx(k,j)*dudz(k,j) + & 416 dwdy(k,j)*dvdz(k,j) ) 322 417 323 418 IF ( def < 0.0_wp ) def = 0.0_wp 324 419 325 tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def 326 327 328 ! 329 !-- (3) - will be executed only, if there is at least one level 330 !-- between (2) and (4), i.e. the topography must have a 331 !-- minimum height of 2 dz. Wall fluxes for this case have 332 !-- already been calculated for (2). 333 !-- 'wall only: use wall functions' 334 335 DO k = nzb_diff_s_inner(j,i), nzb_diff_s_outer(j,i)-2 336 337 dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx 338 dudz = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - & 339 u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k) 340 dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy 341 dvdz = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - & 342 v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k) 343 dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) 344 345 IF ( wall_e_y(j,i) /= 0.0_wp ) THEN 346 ! 347 !-- Inconsistency removed: as the thermal stratification 348 !-- is not taken into account for the evaluation of the 349 !-- wall fluxes at vertical walls, the eddy viscosity km 350 !-- must not be used for the evaluation of the velocity 351 !-- gradients dudy and dwdy 352 !-- Note: The validity of the new method has not yet 353 !-- been shown, as so far no suitable data for a 354 !-- validation has been available 355 km_neutral = kappa * ( usvs(k)**2 + & 356 wsvs(k)**2 )**0.25_wp * 0.5_wp * dy 357 IF ( km_neutral > 0.0_wp ) THEN 358 dudy = - wall_e_y(j,i) * usvs(k) / km_neutral 359 dwdy = - wall_e_y(j,i) * wsvs(k) / km_neutral 360 ELSE 361 dudy = 0.0_wp 362 dwdy = 0.0_wp 363 ENDIF 364 ELSE 365 dudy = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - & 420 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) ) 421 422 tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def * flag 423 424 ENDDO 425 ENDDO 426 427 ELSEIF ( use_surface_fluxes ) THEN 428 429 DO j = nys, nyn 430 ! 431 !-- Calculate TKE production by shear. Here, no additional 432 !-- wall-bounded code is considered. 433 !-- Why? 434 DO k = nzb+1, nzt 435 436 dudx(k,j) = ( u(k,j,i+1) - u(k,j,i) ) * ddx 437 dudy(k,j) = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - & 366 438 u(k,j-1,i) - u(k,j-1,i+1) ) * ddy 367 dwdy = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - & 439 dudz(k,j) = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - & 440 u(k-1,j,i) - u(k-1,j,i+1) ) * & 441 dd2zu(k) 442 443 dvdx(k,j) = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - & 444 v(k,j,i-1) - v(k,j+1,i-1) ) * ddx 445 dvdy(k,j) = ( v(k,j+1,i) - v(k,j,i) ) * ddy 446 dvdz(k,j) = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - & 447 v(k-1,j,i) - v(k-1,j+1,i) ) * & 448 dd2zu(k) 449 450 dwdx(k,j) = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - & 451 w(k,j,i-1) - w(k-1,j,i-1) ) * ddx 452 dwdy(k,j) = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - & 368 453 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy 369 ENDIF 370 371 IF ( wall_e_x(j,i) /= 0.0_wp ) THEN 372 ! 373 !-- Inconsistency removed: as the thermal stratification 374 !-- is not taken into account for the evaluation of the 375 !-- wall fluxes at vertical walls, the eddy viscosity km 376 !-- must not be used for the evaluation of the velocity 377 !-- gradients dvdx and dwdx 378 !-- Note: The validity of the new method has not yet 379 !-- been shown, as so far no suitable data for a 380 !-- validation has been available 381 km_neutral = kappa * ( vsus(k)**2 + & 382 wsus(k)**2 )**0.25_wp * 0.5_wp * dx 383 IF ( km_neutral > 0.0_wp ) THEN 384 dvdx = - wall_e_x(j,i) * vsus(k) / km_neutral 385 dwdx = - wall_e_x(j,i) * wsus(k) / km_neutral 386 ELSE 387 dvdx = 0.0_wp 388 dwdx = 0.0_wp 389 ENDIF 390 ELSE 391 dvdx = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - & 392 v(k,j,i-1) - v(k,j+1,i-1) ) * ddx 393 dwdx = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - & 394 w(k,j,i-1) - w(k-1,j,i-1) ) * ddx 395 ENDIF 396 397 def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 ) + & 398 dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + & 399 dvdz**2 + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz ) 400 401 IF ( def < 0.0_wp ) def = 0.0_wp 402 403 tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def 404 405 ENDDO 406 407 ENDIF 408 409 ENDDO 410 411 ! 412 !-- (4) - will allways be executed. 413 !-- 'special case: free atmosphere' (as for case (0)) 414 DO j = nys, nyn 415 416 IF ( ( wall_e_x(j,i) /= 0.0_wp ) .OR. ( wall_e_y(j,i) /= 0.0_wp ) ) & 417 THEN 418 419 k = nzb_diff_s_outer(j,i)-1 420 421 dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx 422 dudy = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - & 423 u(k,j-1,i) - u(k,j-1,i+1) ) * ddy 424 dudz = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - & 425 u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k) 426 427 dvdx = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - & 428 v(k,j,i-1) - v(k,j+1,i-1) ) * ddx 429 dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy 430 dvdz = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - & 431 v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k) 432 433 dwdx = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - & 434 w(k,j,i-1) - w(k-1,j,i-1) ) * ddx 435 dwdy = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - & 436 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy 437 dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) 438 439 def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 ) + & 440 dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + & 441 dvdz**2 + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz ) 454 dwdz(k,j) = ( w(k,j,i) - w(k-1,j,i) ) * & 455 ddzw(k) 456 457 def = 2.0_wp * ( & 458 dudx(k,j)**2 + dvdy(k,j)**2 + dwdz(k,j)**2 & 459 ) + & 460 dudy(k,j)**2 + dvdx(k,j)**2 + dwdx(k,j)**2 + & 461 dwdy(k,j)**2 + dudz(k,j)**2 + dvdz(k,j)**2 + & 462 2.0_wp * ( & 463 dvdx(k,j)*dudy(k,j) + dwdx(k,j)*dudz(k,j) + & 464 dwdy(k,j)*dvdz(k,j) & 465 ) 442 466 443 467 IF ( def < 0.0_wp ) def = 0.0_wp 444 468 445 tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def 446 447 ENDIF 448 449 ENDDO 450 451 ! 452 !-- Position without adjacent wall 453 !-- (1) - will allways be executed. 454 !-- 'bottom only: use u_0,v_0' 455 DO j = nys, nyn 456 457 IF ( ( wall_e_x(j,i) == 0.0_wp ) .AND. ( wall_e_y(j,i) == 0.0_wp ) ) & 458 THEN 459 460 k = nzb_diff_s_inner(j,i)-1 461 462 dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx 463 dudy = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - & 464 u(k,j-1,i) - u(k,j-1,i+1) ) * ddy 465 dudz = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - & 466 u_0(j,i) - u_0(j,i+1) ) * dd2zu(k) 467 468 dvdx = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - & 469 v(k,j,i-1) - v(k,j+1,i-1) ) * ddx 470 dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy 471 dvdz = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - & 472 v_0(j,i) - v_0(j+1,i) ) * dd2zu(k) 473 474 dwdx = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - & 475 w(k,j,i-1) - w(k-1,j,i-1) ) * ddx 476 dwdy = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - & 477 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy 478 dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) 479 480 def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 ) + & 481 dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + & 482 dvdz**2 + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz ) 483 484 IF ( def < 0.0_wp ) def = 0.0_wp 485 486 tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def 487 488 ENDIF 489 490 ENDDO 491 492 ELSEIF ( use_surface_fluxes ) THEN 493 494 DO j = nys, nyn 495 496 k = nzb_diff_s_outer(j,i)-1 497 498 dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx 499 dudy = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - & 500 u(k,j-1,i) - u(k,j-1,i+1) ) * ddy 501 dudz = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - & 502 u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k) 503 504 dvdx = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - & 505 v(k,j,i-1) - v(k,j+1,i-1) ) * ddx 506 dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy 507 dvdz = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - & 508 v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k) 509 510 dwdx = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - & 511 w(k,j,i-1) - w(k-1,j,i-1) ) * ddx 512 dwdy = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - & 513 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy 514 dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) 515 516 def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 ) + & 517 dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + & 518 dvdz**2 + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz ) 519 520 IF ( def < 0.0_wp ) def = 0.0_wp 521 522 tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def 523 469 flag = MERGE( 1.0_wp, 0.0_wp, & 470 BTEST( wall_flags_0(k,j,i), 29 ) ) 471 tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def * flag 472 473 ENDDO 524 474 ENDDO 525 475 … … 539 489 !-- in the bottom and top surface layer 540 490 DO j = nys, nyn 541 DO k = nzb _s_inner(j,i)+1, nzt491 DO k = nzb+1, nzt 542 492 tend(k,j,i) = tend(k,j,i) + & 543 493 kh(k,j,i) * g / rho_reference * & 544 494 ( prho(k+1,j,i) - prho(k-1,j,i) ) * & 545 dd2zu(k) 495 dd2zu(k) * & 496 MERGE( 1.0_wp, 0.0_wp, & 497 BTEST( wall_flags_0(k,j,i), 0 ) & 498 ) 546 499 ENDDO 547 500 ENDDO … … 550 503 551 504 DO j = nys, nyn 552 DO k = nzb_diff_s_inner(j,i), nzt_diff 553 tend(k,j,i) = tend(k,j,i) - & 554 kh(k,j,i) * g / pt_reference * & 555 ( pt(k+1,j,i) - pt(k-1,j,i) ) * & 556 dd2zu(k) 505 DO k = nzb+1, nzt 506 ! 507 !-- Flag 9 is used to mask top fluxes, flag 30 to mask 508 !-- surface fluxes 509 tend(k,j,i) = tend(k,j,i) - & 510 kh(k,j,i) * g / pt_reference * & 511 ( pt(k+1,j,i) - pt(k-1,j,i) ) * & 512 dd2zu(k) * & 513 MERGE( 1.0_wp, 0.0_wp, & 514 BTEST( wall_flags_0(k,j,i), 30 ) & 515 ) * & 516 MERGE( 1.0_wp, 0.0_wp, & 517 BTEST( wall_flags_0(k,j,i), 9 ) & 518 ) 557 519 ENDDO 558 520 559 521 IF ( use_surface_fluxes ) THEN 560 k = nzb_diff_s_inner(j,i)-1 561 tend(k,j,i) = tend(k,j,i) + g / pt_reference * & 562 shf(j,i) 522 ! 523 !-- Default surfaces, up- and downward-facing 524 DO l = 0, 1 525 surf_s = surf_def_h(l)%start_index(j,i) 526 surf_e = surf_def_h(l)%end_index(j,i) 527 DO m = surf_s, surf_e 528 k = surf_def_h(l)%k(m) 529 tend(k,j,i) = tend(k,j,i) + g / pt_reference & 530 * surf_def_h(l)%shf(m) 531 ENDDO 532 ENDDO 533 ! 534 !-- Natural surfaces 535 surf_s = surf_lsm_h%start_index(j,i) 536 surf_e = surf_lsm_h%end_index(j,i) 537 DO m = surf_s, surf_e 538 k = surf_lsm_h%k(m) 539 tend(k,j,i) = tend(k,j,i) + g / pt_reference & 540 * surf_lsm_h%shf(m) 541 ENDDO 542 ! 543 !-- Urban surfaces 544 surf_s = surf_usm_h%start_index(j,i) 545 surf_e = surf_usm_h%end_index(j,i) 546 DO m = surf_s, surf_e 547 k = surf_usm_h%k(m) 548 tend(k,j,i) = tend(k,j,i) + g / pt_reference & 549 * surf_usm_h%shf(m) 550 ENDDO 563 551 ENDIF 564 552 565 553 IF ( use_top_fluxes ) THEN 566 k = nzt 567 tend(k,j,i) = tend(k,j,i) + g / pt_reference * & 568 tswst(j,i) 554 surf_s = surf_def_h(2)%start_index(j,i) 555 surf_e = surf_def_h(2)%end_index(j,i) 556 DO m = surf_s, surf_e 557 k = surf_def_h(2)%k(m) 558 tend(k,j,i) = tend(k,j,i) + g / pt_reference * & 559 surf_def_h(2)%shf(m) 560 ENDDO 569 561 ENDIF 570 562 ENDDO … … 579 571 !-- in the bottom and top surface layer 580 572 DO j = nys, nyn 581 DO k = nzb _s_inner(j,i)+1, nzt573 DO k = nzb+1, nzt 582 574 tend(k,j,i) = tend(k,j,i) + & 583 575 kh(k,j,i) * g / prho(k,j,i) * & 584 576 ( prho(k+1,j,i) - prho(k-1,j,i) ) * & 585 dd2zu(k) 577 dd2zu(k) * & 578 MERGE( 1.0_wp, 0.0_wp, & 579 BTEST( wall_flags_0(k,j,i), 0 ) & 580 ) 586 581 ENDDO 587 582 ENDDO … … 590 585 591 586 DO j = nys, nyn 592 DO k = nzb_diff_s_inner(j,i), nzt_diff 593 tend(k,j,i) = tend(k,j,i) - & 594 kh(k,j,i) * g / pt(k,j,i) * & 595 ( pt(k+1,j,i) - pt(k-1,j,i) ) * & 596 dd2zu(k) 587 DO k = nzb+1, nzt 588 ! 589 !-- Flag 9 is used to mask top fluxes, flag 30 to mask 590 !-- surface fluxes 591 tend(k,j,i) = tend(k,j,i) - & 592 kh(k,j,i) * g / pt(k,j,i) * & 593 ( pt(k+1,j,i) - pt(k-1,j,i) ) * & 594 dd2zu(k) * & 595 MERGE( 1.0_wp, 0.0_wp, & 596 BTEST( wall_flags_0(k,j,i), 30 ) & 597 ) * & 598 MERGE( 1.0_wp, 0.0_wp, & 599 BTEST( wall_flags_0(k,j,i), 9 ) & 600 ) 597 601 ENDDO 598 602 599 603 IF ( use_surface_fluxes ) THEN 600 k = nzb_diff_s_inner(j,i)-1 601 tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * & 602 shf(j,i) 604 ! 605 !-- Default surfaces, up- and downwrd-facing 606 DO l = 0, 1 607 surf_s = surf_def_h(l)%start_index(j,i) 608 surf_e = surf_def_h(l)%end_index(j,i) 609 DO m = surf_s, surf_e 610 k = surf_def_h(l)%k(m) 611 tend(k,j,i) = tend(k,j,i) + g / pt_reference & 612 * surf_def_h(l)%shf(m) 613 ENDDO 614 ENDDO 615 ! 616 !-- Natural surfaces 617 surf_s = surf_lsm_h%start_index(j,i) 618 surf_e = surf_lsm_h%end_index(j,i) 619 DO m = surf_s, surf_e 620 k = surf_lsm_h%k(m) 621 tend(k,j,i) = tend(k,j,i) + g / pt_reference & 622 * surf_lsm_h%shf(m) 623 ENDDO 624 ! 625 !-- Urban surfaces 626 surf_s = surf_usm_h%start_index(j,i) 627 surf_e = surf_usm_h%end_index(j,i) 628 DO m = surf_s, surf_e 629 k = surf_usm_h%k(m) 630 tend(k,j,i) = tend(k,j,i) + g / pt_reference & 631 * surf_usm_h%shf(m) 632 ENDDO 603 633 ENDIF 604 634 605 635 IF ( use_top_fluxes ) THEN 606 k = nzt 607 tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * & 608 tswst(j,i) 636 surf_s = surf_def_h(2)%start_index(j,i) 637 surf_e = surf_def_h(2)%end_index(j,i) 638 DO m = surf_s, surf_e 639 k = surf_def_h(2)%k(m) 640 tend(k,j,i) = tend(k,j,i) + g / pt_reference * & 641 surf_def_h(2)%shf(m) 642 ENDDO 609 643 ENDIF 610 644 ENDDO … … 618 652 DO j = nys, nyn 619 653 620 DO k = nzb_diff_s_inner(j,i), nzt_diff 621 654 DO k = nzb+1, nzt 655 ! 656 !-- Flag 9 is used to mask top fluxes, flag 30 to mask 657 !-- surface fluxes 622 658 IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN 623 659 k1 = 1.0_wp + 0.61_wp * q(k,j,i) … … 627 663 ( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) + & 628 664 k2 * ( q(k+1,j,i) - q(k-1,j,i) ) & 629 ) * dd2zu(k) 630 ELSE IF ( cloud_physics ) THEN 631 IF ( ql(k,j,i) == 0.0_wp ) THEN 632 k1 = 1.0_wp + 0.61_wp * q(k,j,i) 633 k2 = 0.61_wp * pt(k,j,i) 634 ELSE 635 theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i) 636 temp = theta * t_d_pt(k) 637 k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * & 638 ( q(k,j,i) - ql(k,j,i) ) * & 639 ( 1.0_wp + 0.622_wp * l_d_r / temp ) ) / & 640 ( 1.0_wp + 0.622_wp * l_d_r * l_d_cp * & 641 ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) ) 642 k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp ) 643 ENDIF 644 tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * & 645 g / vpt(k,j,i) * & 646 ( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) + & 647 k2 * ( q(k+1,j,i) - q(k-1,j,i) ) & 648 ) * dd2zu(k) 649 ELSE IF ( cloud_droplets ) THEN 650 k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i) 651 k2 = 0.61_wp * pt(k,j,i) 652 tend(k,j,i) = tend(k,j,i) - & 653 kh(k,j,i) * g / vpt(k,j,i) * & 654 ( k1 * ( pt(k+1,j,i)- pt(k-1,j,i) ) + & 655 k2 * ( q(k+1,j,i) - q(k-1,j,i) ) - & 656 pt(k,j,i) * ( ql(k+1,j,i) - & 657 ql(k-1,j,i) ) ) * dd2zu(k) 658 ENDIF 659 660 ENDDO 661 662 ENDDO 663 664 IF ( use_surface_fluxes ) THEN 665 666 DO j = nys, nyn 667 668 k = nzb_diff_s_inner(j,i)-1 669 670 IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN 671 k1 = 1.0_wp + 0.61_wp * q(k,j,i) 672 k2 = 0.61_wp * pt(k,j,i) 665 ) * dd2zu(k) * & 666 MERGE( 1.0_wp, 0.0_wp, & 667 BTEST( wall_flags_0(k,j,i), 30 ) & 668 ) * & 669 MERGE( 1.0_wp, 0.0_wp, & 670 BTEST( wall_flags_0(k,j,i), 9 ) & 671 ) 673 672 ELSE IF ( cloud_physics ) THEN 674 673 IF ( ql(k,j,i) == 0.0_wp ) THEN … … 685 684 k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp ) 686 685 ENDIF 686 tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * & 687 g / vpt(k,j,i) * & 688 ( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) + & 689 k2 * ( q(k+1,j,i) - q(k-1,j,i) ) & 690 ) * dd2zu(k) * & 691 MERGE( 1.0_wp, 0.0_wp, & 692 BTEST( wall_flags_0(k,j,i), 30 ) & 693 ) * & 694 MERGE( 1.0_wp, 0.0_wp, & 695 BTEST( wall_flags_0(k,j,i), 9 ) & 696 ) 687 697 ELSE IF ( cloud_droplets ) THEN 688 698 k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i) 689 699 k2 = 0.61_wp * pt(k,j,i) 700 tend(k,j,i) = tend(k,j,i) - & 701 kh(k,j,i) * g / vpt(k,j,i) * & 702 ( k1 * ( pt(k+1,j,i)- pt(k-1,j,i) ) + & 703 k2 * ( q(k+1,j,i) - q(k-1,j,i) ) - & 704 pt(k,j,i) * ( ql(k+1,j,i) - & 705 ql(k-1,j,i) ) ) * dd2zu(k) * & 706 MERGE( 1.0_wp, 0.0_wp, & 707 BTEST( wall_flags_0(k,j,i), 30 ) & 708 ) * & 709 MERGE( 1.0_wp, 0.0_wp, & 710 BTEST( wall_flags_0(k,j,i), 9 ) & 711 ) 690 712 ENDIF 691 713 692 tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) * &693 ( k1* shf(j,i) + k2 * qsws(j,i) )694 714 ENDDO 695 715 716 ENDDO 717 718 IF ( use_surface_fluxes ) THEN 719 720 DO j = nys, nyn 721 ! 722 !-- Treat horizontal default surfaces, up- and downward-facing 723 DO l = 0, 1 724 surf_s = surf_def_h(l)%start_index(j,i) 725 surf_e = surf_def_h(l)%end_index(j,i) 726 DO m = surf_s, surf_e 727 k = surf_def_h(l)%k(m) 728 729 IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN 730 k1 = 1.0_wp + 0.61_wp * q(k,j,i) 731 k2 = 0.61_wp * pt(k,j,i) 732 ELSE IF ( cloud_physics ) THEN 733 IF ( ql(k,j,i) == 0.0_wp ) THEN 734 k1 = 1.0_wp + 0.61_wp * q(k,j,i) 735 k2 = 0.61_wp * pt(k,j,i) 736 ELSE 737 theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i) 738 temp = theta * t_d_pt(k) 739 k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * & 740 ( q(k,j,i) - ql(k,j,i) ) * & 741 ( 1.0_wp + 0.622_wp * l_d_r / temp ) ) / & 742 ( 1.0_wp + 0.622_wp * l_d_r * l_d_cp * & 743 ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) ) 744 k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp ) 745 ENDIF 746 ELSE IF ( cloud_droplets ) THEN 747 k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i) 748 k2 = 0.61_wp * pt(k,j,i) 749 ENDIF 750 751 tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) * & 752 ( k1 * surf_def_h(l)%shf(m) + & 753 k2 * surf_def_h(l)%qsws(m) ) 754 ENDDO 755 ENDDO 756 ! 757 !-- Treat horizontal natural surfaces 758 surf_s = surf_lsm_h%start_index(j,i) 759 surf_e = surf_lsm_h%end_index(j,i) 760 DO m = surf_s, surf_e 761 k = surf_lsm_h%k(m) 762 763 IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN 764 k1 = 1.0_wp + 0.61_wp * q(k,j,i) 765 k2 = 0.61_wp * pt(k,j,i) 766 ELSE IF ( cloud_physics ) THEN 767 IF ( ql(k,j,i) == 0.0_wp ) THEN 768 k1 = 1.0_wp + 0.61_wp * q(k,j,i) 769 k2 = 0.61_wp * pt(k,j,i) 770 ELSE 771 theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i) 772 temp = theta * t_d_pt(k) 773 k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * & 774 ( q(k,j,i) - ql(k,j,i) ) * & 775 ( 1.0_wp + 0.622_wp * l_d_r / temp ) ) / & 776 ( 1.0_wp + 0.622_wp * l_d_r * l_d_cp * & 777 ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) ) 778 k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp ) 779 ENDIF 780 ELSE IF ( cloud_droplets ) THEN 781 k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i) 782 k2 = 0.61_wp * pt(k,j,i) 783 ENDIF 784 785 tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) * & 786 ( k1 * surf_lsm_h%shf(m) + & 787 k2 * surf_lsm_h%qsws(m) ) 788 ENDDO 789 ! 790 !-- Treat horizontal urban surfaces 791 surf_s = surf_usm_h%start_index(j,i) 792 surf_e = surf_usm_h%end_index(j,i) 793 DO m = surf_s, surf_e 794 k = surf_lsm_h%k(m) 795 796 IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN 797 k1 = 1.0_wp + 0.61_wp * q(k,j,i) 798 k2 = 0.61_wp * pt(k,j,i) 799 ELSE IF ( cloud_physics ) THEN 800 IF ( ql(k,j,i) == 0.0_wp ) THEN 801 k1 = 1.0_wp + 0.61_wp * q(k,j,i) 802 k2 = 0.61_wp * pt(k,j,i) 803 ELSE 804 theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i) 805 temp = theta * t_d_pt(k) 806 k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * & 807 ( q(k,j,i) - ql(k,j,i) ) * & 808 ( 1.0_wp + 0.622_wp * l_d_r / temp ) ) / & 809 ( 1.0_wp + 0.622_wp * l_d_r * l_d_cp * & 810 ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) ) 811 k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp ) 812 ENDIF 813 ELSE IF ( cloud_droplets ) THEN 814 k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i) 815 k2 = 0.61_wp * pt(k,j,i) 816 ENDIF 817 818 tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) * & 819 ( k1 * surf_usm_h%shf(m) + & 820 k2 * surf_usm_h%qsws(m) ) 821 ENDDO 822 823 ENDDO 824 696 825 ENDIF 697 826 … … 700 829 DO j = nys, nyn 701 830 702 k = nzt 831 surf_s = surf_def_h(2)%start_index(j,i) 832 surf_e = surf_def_h(2)%end_index(j,i) 833 DO m = surf_s, surf_e 834 k = surf_def_h(2)%k(m) 835 836 IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN 837 k1 = 1.0_wp + 0.61_wp * q(k,j,i) 838 k2 = 0.61_wp * pt(k,j,i) 839 ELSE IF ( cloud_physics ) THEN 840 IF ( ql(k,j,i) == 0.0_wp ) THEN 841 k1 = 1.0_wp + 0.61_wp * q(k,j,i) 842 k2 = 0.61_wp * pt(k,j,i) 843 ELSE 844 theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i) 845 temp = theta * t_d_pt(k) 846 k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * & 847 ( q(k,j,i) - ql(k,j,i) ) * & 848 ( 1.0_wp + 0.622_wp * l_d_r / temp ) ) / & 849 ( 1.0_wp + 0.622_wp * l_d_r * l_d_cp * & 850 ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) ) 851 k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp ) 852 ENDIF 853 ELSE IF ( cloud_droplets ) THEN 854 k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i) 855 k2 = 0.61_wp * pt(k,j,i) 856 ENDIF 857 858 tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) * & 859 ( k1 * surf_def_h(2)%shf(m) + & 860 k2 * surf_def_h(2)%qsws(m) ) 861 862 ENDDO 863 864 ENDDO 865 866 ENDIF 867 868 ENDIF 869 870 ENDIF 871 872 ENDDO 873 874 END SUBROUTINE production_e 875 876 877 !------------------------------------------------------------------------------! 878 ! Description: 879 ! ------------ 880 !> Call for grid point i,j 881 !------------------------------------------------------------------------------! 882 SUBROUTINE production_e_ij( i, j ) 883 884 USE arrays_3d, & 885 ONLY: ddzw, dd2zu, kh, km, prho, pt, q, ql, tend, u, v, vpt, w 886 887 USE cloud_parameters, & 888 ONLY: l_d_cp, l_d_r, pt_d_t, t_d_pt 889 890 USE control_parameters, & 891 ONLY: cloud_droplets, cloud_physics, constant_flux_layer, g, & 892 humidity, kappa, neutral, ocean, pt_reference, & 893 rho_reference, use_single_reference_value, & 894 use_surface_fluxes, use_top_fluxes 895 896 USE grid_variables, & 897 ONLY: ddx, dx, ddy, dy 898 899 USE indices, & 900 ONLY: nxl, nxr, nys, nyn, nzb, nzb, nzt, wall_flags_0 901 902 USE surface_mod, & 903 ONLY : surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, & 904 surf_usm_v 905 906 IMPLICIT NONE 907 908 INTEGER(iwp) :: i !< running index x-direction 909 INTEGER(iwp) :: j !< running index y-direction 910 INTEGER(iwp) :: k !< running index z-direction 911 INTEGER(iwp) :: l !< running index for different surface type orientation 912 INTEGER(iwp) :: m !< running index surface elements 913 INTEGER(iwp) :: surf_e !< end index of surface elements at given i-j position 914 INTEGER(iwp) :: surf_s !< start index of surface elements at given i-j position 915 916 REAL(wp) :: def !< 917 REAL(wp) :: flag !< flag to mask topography 918 REAL(wp) :: k1 !< 919 REAL(wp) :: k2 !< 920 REAL(wp) :: km_neutral !< diffusion coefficient assuming neutral conditions - used to compute shear production at surfaces 921 REAL(wp) :: theta !< 922 REAL(wp) :: temp !< 923 REAL(wp) :: sign_dir !< sign of wall-tke flux, depending on wall orientation 924 REAL(wp) :: usvs !< momentum flux u"v" 925 REAL(wp) :: vsus !< momentum flux v"u" 926 REAL(wp) :: wsus !< momentum flux w"u" 927 REAL(wp) :: wsvs !< momentum flux w"v" 928 929 930 REAL(wp), DIMENSION(nzb+1:nzt) :: dudx !< Gradient of u-component in x-direction 931 REAL(wp), DIMENSION(nzb+1:nzt) :: dudy !< Gradient of u-component in y-direction 932 REAL(wp), DIMENSION(nzb+1:nzt) :: dudz !< Gradient of u-component in z-direction 933 REAL(wp), DIMENSION(nzb+1:nzt) :: dvdx !< Gradient of v-component in x-direction 934 REAL(wp), DIMENSION(nzb+1:nzt) :: dvdy !< Gradient of v-component in y-direction 935 REAL(wp), DIMENSION(nzb+1:nzt) :: dvdz !< Gradient of v-component in z-direction 936 REAL(wp), DIMENSION(nzb+1:nzt) :: dwdx !< Gradient of w-component in x-direction 937 REAL(wp), DIMENSION(nzb+1:nzt) :: dwdy !< Gradient of w-component in y-direction 938 REAL(wp), DIMENSION(nzb+1:nzt) :: dwdz !< Gradient of w-component in z-direction 939 940 941 IF ( constant_flux_layer ) THEN 942 ! 943 !-- Calculate TKE production by shear. Calculate gradients at all grid 944 !-- points first, gradients at surface-bounded grid points will be 945 !-- overwritten further below. 946 DO k = nzb+1, nzt 947 948 dudx(k) = ( u(k,j,i+1) - u(k,j,i) ) * ddx 949 dudy(k) = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - & 950 u(k,j-1,i) - u(k,j-1,i+1) ) * ddy 951 dudz(k) = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - & 952 u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k) 953 954 dvdx(k) = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - & 955 v(k,j,i-1) - v(k,j+1,i-1) ) * ddx 956 dvdy(k) = ( v(k,j+1,i) - v(k,j,i) ) * ddy 957 dvdz(k) = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - & 958 v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k) 959 960 dwdx(k) = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - & 961 w(k,j,i-1) - w(k-1,j,i-1) ) * ddx 962 dwdy(k) = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - & 963 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy 964 dwdz(k) = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) 965 966 ENDDO 967 ! 968 !-- Compute gradients at north- and south-facing surfaces. 969 !-- Note, no vertical natural surfaces so far. 970 DO l = 0, 1 971 ! 972 !-- Default surfaces 973 surf_s = surf_def_v(l)%start_index(j,i) 974 surf_e = surf_def_v(l)%end_index(j,i) 975 DO m = surf_s, surf_e 976 k = surf_def_v(l)%k(m) 977 usvs = surf_def_v(l)%mom_flux_tke(0,m) 978 wsvs = surf_def_v(l)%mom_flux_tke(1,m) 979 980 km_neutral = kappa * ( usvs**2 + wsvs**2 )**0.25_wp & 981 * 0.5_wp * dy 982 ! 983 !-- -1.0 for right-facing wall, 1.0 for left-facing wall 984 sign_dir = MERGE( 1.0_wp, -1.0_wp, & 985 BTEST( wall_flags_0(k,j-1,i), 0 ) ) 986 dudy(k) = sign_dir * usvs / ( km_neutral + 1E-10_wp ) 987 dwdy(k) = sign_dir * wsvs / ( km_neutral + 1E-10_wp ) 988 ENDDO 989 ! 990 !-- Natural surfaces 991 surf_s = surf_lsm_v(l)%start_index(j,i) 992 surf_e = surf_lsm_v(l)%end_index(j,i) 993 DO m = surf_s, surf_e 994 k = surf_lsm_v(l)%k(m) 995 usvs = surf_lsm_v(l)%mom_flux_tke(0,m) 996 wsvs = surf_lsm_v(l)%mom_flux_tke(1,m) 997 998 km_neutral = kappa * ( usvs**2 + wsvs**2 )**0.25_wp & 999 * 0.5_wp * dy 1000 ! 1001 !-- -1.0 for right-facing wall, 1.0 for left-facing wall 1002 sign_dir = MERGE( 1.0_wp, -1.0_wp, & 1003 BTEST( wall_flags_0(k,j-1,i), 0 ) ) 1004 dudy(k) = sign_dir * usvs / ( km_neutral + 1E-10_wp ) 1005 dwdy(k) = sign_dir * wsvs / ( km_neutral + 1E-10_wp ) 1006 ENDDO 1007 ! 1008 !-- Urban surfaces 1009 surf_s = surf_usm_v(l)%start_index(j,i) 1010 surf_e = surf_usm_v(l)%end_index(j,i) 1011 DO m = surf_s, surf_e 1012 k = surf_usm_v(l)%k(m) 1013 usvs = surf_usm_v(l)%mom_flux_tke(0,m) 1014 wsvs = surf_usm_v(l)%mom_flux_tke(1,m) 1015 1016 km_neutral = kappa * ( usvs**2 + wsvs**2 )**0.25_wp & 1017 * 0.5_wp * dy 1018 ! 1019 !-- -1.0 for right-facing wall, 1.0 for left-facing wall 1020 sign_dir = MERGE( 1.0_wp, -1.0_wp, & 1021 BTEST( wall_flags_0(k,j-1,i), 0 ) ) 1022 dudy(k) = sign_dir * usvs / ( km_neutral + 1E-10_wp ) 1023 dwdy(k) = sign_dir * wsvs / ( km_neutral + 1E-10_wp ) 1024 ENDDO 1025 ENDDO 1026 ! 1027 !-- Compute gradients at east- and west-facing walls 1028 DO l = 2, 3 1029 ! 1030 !-- Default surfaces 1031 surf_s = surf_def_v(l)%start_index(j,i) 1032 surf_e = surf_def_v(l)%end_index(j,i) 1033 DO m = surf_s, surf_e 1034 k = surf_def_v(l)%k(m) 1035 vsus = surf_def_v(l)%mom_flux_tke(0,m) 1036 wsus = surf_def_v(l)%mom_flux_tke(1,m) 1037 1038 km_neutral = kappa * ( vsus**2 + wsus**2 )**0.25_wp & 1039 * 0.5_wp * dx 1040 ! 1041 !-- -1.0 for right-facing wall, 1.0 for left-facing wall 1042 sign_dir = MERGE( 1.0_wp, -1.0_wp, & 1043 BTEST( wall_flags_0(k,j,i-1), 0 ) ) 1044 dvdx(k) = sign_dir * vsus / ( km_neutral + 1E-10_wp ) 1045 dwdx(k) = sign_dir * wsus / ( km_neutral + 1E-10_wp ) 1046 ENDDO 1047 ! 1048 !-- Natural surfaces 1049 surf_s = surf_lsm_v(l)%start_index(j,i) 1050 surf_e = surf_lsm_v(l)%end_index(j,i) 1051 DO m = surf_s, surf_e 1052 k = surf_lsm_v(l)%k(m) 1053 vsus = surf_lsm_v(l)%mom_flux_tke(0,m) 1054 wsus = surf_lsm_v(l)%mom_flux_tke(1,m) 1055 1056 km_neutral = kappa * ( vsus**2 + wsus**2 )**0.25_wp & 1057 * 0.5_wp * dx 1058 ! 1059 !-- -1.0 for right-facing wall, 1.0 for left-facing wall 1060 sign_dir = MERGE( 1.0_wp, -1.0_wp, & 1061 BTEST( wall_flags_0(k,j,i-1), 0 ) ) 1062 dvdx(k) = sign_dir * vsus / ( km_neutral + 1E-10_wp ) 1063 dwdx(k) = sign_dir * wsus / ( km_neutral + 1E-10_wp ) 1064 ENDDO 1065 ! 1066 !-- Urban surfaces 1067 surf_s = surf_usm_v(l)%start_index(j,i) 1068 surf_e = surf_usm_v(l)%end_index(j,i) 1069 DO m = surf_s, surf_e 1070 k = surf_usm_v(l)%k(m) 1071 vsus = surf_usm_v(l)%mom_flux_tke(0,m) 1072 wsus = surf_usm_v(l)%mom_flux_tke(1,m) 1073 1074 km_neutral = kappa * ( vsus**2 + wsus**2 )**0.25_wp & 1075 * 0.5_wp * dx 1076 ! 1077 !-- -1.0 for right-facing wall, 1.0 for left-facing wall 1078 sign_dir = MERGE( 1.0_wp, -1.0_wp, & 1079 BTEST( wall_flags_0(k,j,i-1), 0 ) ) 1080 dvdx(k) = sign_dir * vsus / ( km_neutral + 1E-10_wp ) 1081 dwdx(k) = sign_dir * wsus / ( km_neutral + 1E-10_wp ) 1082 ENDDO 1083 ENDDO 1084 ! 1085 !-- Compute gradients at upward-facing walls, first for 1086 !-- non-natural default surfaces 1087 surf_s = surf_def_h(0)%start_index(j,i) 1088 surf_e = surf_def_h(0)%end_index(j,i) 1089 DO m = surf_s, surf_e 1090 k = surf_def_h(0)%k(m) 1091 ! 1092 !-- Please note, actually, an interpolation of u_0 and v_0 1093 !-- onto the grid center would be required. However, this 1094 !-- would require several data transfers between 2D-grid and 1095 !-- wall type. The effect of this missing interpolation is 1096 !-- negligible. (See also production_e_init). 1097 dudz(k) = ( u(k+1,j,i) - surf_def_h(0)%u_0(m) ) * dd2zu(k) 1098 dvdz(k) = ( v(k+1,j,i) - surf_def_h(0)%v_0(m) ) * dd2zu(k) 1099 1100 ENDDO 1101 ! 1102 !-- Natural surfaces 1103 surf_s = surf_lsm_h%start_index(j,i) 1104 surf_e = surf_lsm_h%end_index(j,i) 1105 DO m = surf_s, surf_e 1106 k = surf_lsm_h%k(m) 1107 ! 1108 !-- Please note, actually, an interpolation of u_0 and v_0 1109 !-- onto the grid center would be required. However, this 1110 !-- would require several data transfers between 2D-grid and 1111 !-- wall type. The effect of this missing interpolation is 1112 !-- negligible. (See also production_e_init). 1113 dudz(k) = ( u(k+1,j,i) - surf_lsm_h%u_0(m) ) * dd2zu(k) 1114 dvdz(k) = ( v(k+1,j,i) - surf_lsm_h%v_0(m) ) * dd2zu(k) 1115 ENDDO 1116 ! 1117 !-- Urban surfaces 1118 surf_s = surf_usm_h%start_index(j,i) 1119 surf_e = surf_usm_h%end_index(j,i) 1120 DO m = surf_s, surf_e 1121 k = surf_usm_h%k(m) 1122 ! 1123 !-- Please note, actually, an interpolation of u_0 and v_0 1124 !-- onto the grid center would be required. However, this 1125 !-- would require several data transfers between 2D-grid and 1126 !-- wall type. The effect of this missing interpolation is 1127 !-- negligible. (See also production_e_init). 1128 dudz(k) = ( u(k+1,j,i) - surf_usm_h%u_0(m) ) * dd2zu(k) 1129 dvdz(k) = ( v(k+1,j,i) - surf_usm_h%v_0(m) ) * dd2zu(k) 1130 ENDDO 1131 ! 1132 !-- Compute gradients at downward-facing walls, only for 1133 !-- non-natural default surfaces 1134 surf_s = surf_def_h(1)%start_index(j,i) 1135 surf_e = surf_def_h(1)%end_index(j,i) 1136 DO m = surf_s, surf_e 1137 k = surf_def_h(1)%k(m) 1138 ! 1139 !-- Please note, actually, an interpolation of u_0 and v_0 1140 !-- onto the grid center would be required. However, this 1141 !-- would require several data transfers between 2D-grid and 1142 !-- wall type. The effect of this missing interpolation is 1143 !-- negligible. (See also production_e_init). 1144 dudz(k) = ( surf_def_h(1)%u_0(m) - u(k-1,j,i) ) * dd2zu(k) 1145 dvdz(k) = ( surf_def_h(1)%v_0(m) - v(k-1,j,i) ) * dd2zu(k) 1146 1147 ENDDO 1148 1149 DO k = nzb+1, nzt 1150 1151 def = 2.0_wp * ( dudx(k)**2 + dvdy(k)**2 + dwdz(k)**2 ) + & 1152 dudy(k)**2 + dvdx(k)**2 + dwdx(k)**2 + & 1153 dwdy(k)**2 + dudz(k)**2 + dvdz(k)**2 + & 1154 2.0_wp * ( dvdx(k)*dudy(k) + dwdx(k)*dudz(k) + dwdy(k)*dvdz(k) ) 1155 1156 IF ( def < 0.0_wp ) def = 0.0_wp 1157 1158 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) ) 1159 1160 tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def * flag 1161 1162 ENDDO 1163 1164 ELSEIF ( use_surface_fluxes ) THEN 1165 ! 1166 !-- Calculate TKE production by shear. Here, no additional 1167 !-- wall-bounded code is considered. 1168 !-- Why? 1169 DO k = nzb+1, nzt 1170 1171 dudx(k) = ( u(k,j,i+1) - u(k,j,i) ) * ddx 1172 dudy(k) = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - & 1173 u(k,j-1,i) - u(k,j-1,i+1) ) * ddy 1174 dudz(k) = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - & 1175 u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k) 1176 1177 dvdx(k) = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - & 1178 v(k,j,i-1) - v(k,j+1,i-1) ) * ddx 1179 dvdy(k) = ( v(k,j+1,i) - v(k,j,i) ) * ddy 1180 dvdz(k) = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - & 1181 v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k) 1182 1183 dwdx(k) = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - & 1184 w(k,j,i-1) - w(k-1,j,i-1) ) * ddx 1185 dwdy(k) = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - & 1186 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy 1187 dwdz(k) = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) 1188 1189 def = 2.0_wp * ( dudx(k)**2 + dvdy(k)**2 + dwdz(k)**2 ) + & 1190 dudy(k)**2 + dvdx(k)**2 + dwdx(k)**2 + & 1191 dwdy(k)**2 + dudz(k)**2 + dvdz(k)**2 + & 1192 2.0_wp * ( dvdx(k)*dudy(k) + dwdx(k)*dudz(k) + dwdy(k)*dvdz(k) ) 1193 1194 IF ( def < 0.0_wp ) def = 0.0_wp 1195 1196 flag = MERGE( 1.0_wp, 0.0_wp, & 1197 BTEST( wall_flags_0(k,j,i), 29 ) ) 1198 tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def * flag 1199 1200 ENDDO 1201 1202 ENDIF 1203 1204 ! 1205 !-- If required, calculate TKE production by buoyancy 1206 IF ( .NOT. neutral ) THEN 1207 1208 IF ( .NOT. humidity ) THEN 1209 1210 IF ( use_single_reference_value ) THEN 1211 1212 IF ( ocean ) THEN 1213 ! 1214 !-- So far in the ocean no special treatment of density flux in 1215 !-- the bottom and top surface layer 1216 DO k = nzb+1, nzt 1217 tend(k,j,i) = tend(k,j,i) + & 1218 kh(k,j,i) * g / rho_reference * & 1219 ( prho(k+1,j,i) - prho(k-1,j,i) ) * & 1220 dd2zu(k) * & 1221 MERGE( 1.0_wp, 0.0_wp, & 1222 BTEST( wall_flags_0(k,j,i), 0 ) & 1223 ) 1224 ENDDO 1225 1226 ELSE 1227 1228 DO k = nzb+1, nzt 1229 ! 1230 !-- Flag 9 is used to mask top fluxes, flag 30 to mask 1231 !-- surface fluxes 1232 tend(k,j,i) = tend(k,j,i) - & 1233 kh(k,j,i) * g / pt_reference * & 1234 ( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k) * & 1235 MERGE( 1.0_wp, 0.0_wp, & 1236 BTEST( wall_flags_0(k,j,i), 30 ) & 1237 ) * & 1238 MERGE( 1.0_wp, 0.0_wp, & 1239 BTEST( wall_flags_0(k,j,i), 9 ) & 1240 ) 1241 1242 ENDDO 1243 1244 IF ( use_surface_fluxes ) THEN 1245 ! 1246 !-- Default surfaces, up- and downward-facing 1247 DO l = 0, 1 1248 surf_s = surf_def_h(l)%start_index(j,i) 1249 surf_e = surf_def_h(l)%end_index(j,i) 1250 DO m = surf_s, surf_e 1251 k = surf_def_h(l)%k(m) 1252 tend(k,j,i) = tend(k,j,i) + g / pt_reference * & 1253 surf_def_h(l)%shf(m) 1254 ENDDO 1255 ENDDO 1256 ! 1257 !-- Natural surfaces 1258 surf_s = surf_lsm_h%start_index(j,i) 1259 surf_e = surf_lsm_h%end_index(j,i) 1260 DO m = surf_s, surf_e 1261 k = surf_lsm_h%k(m) 1262 tend(k,j,i) = tend(k,j,i) + g / pt_reference * & 1263 surf_lsm_h%shf(m) 1264 ENDDO 1265 ! 1266 !-- Urban surfaces 1267 surf_s = surf_usm_h%start_index(j,i) 1268 surf_e = surf_usm_h%end_index(j,i) 1269 DO m = surf_s, surf_e 1270 k = surf_usm_h%k(m) 1271 tend(k,j,i) = tend(k,j,i) + g / pt_reference * & 1272 surf_usm_h%shf(m) 1273 ENDDO 1274 ENDIF 1275 1276 IF ( use_top_fluxes ) THEN 1277 surf_s = surf_def_h(2)%start_index(j,i) 1278 surf_e = surf_def_h(2)%end_index(j,i) 1279 DO m = surf_s, surf_e 1280 k = surf_def_h(2)%k(m) 1281 tend(k,j,i) = tend(k,j,i) + g / pt_reference * & 1282 surf_def_h(2)%shf(m) 1283 ENDDO 1284 ENDIF 1285 1286 ENDIF 1287 1288 ELSE 1289 1290 IF ( ocean ) THEN 1291 ! 1292 !-- So far in the ocean no special treatment of density flux in 1293 !-- the bottom and top surface layer 1294 DO k = nzb+1, nzt 1295 tend(k,j,i) = tend(k,j,i) + & 1296 kh(k,j,i) * g / prho(k,j,i) * & 1297 ( prho(k+1,j,i) - prho(k-1,j,i) ) * & 1298 dd2zu(k) * & 1299 MERGE( 1.0_wp, 0.0_wp, & 1300 BTEST( wall_flags_0(k,j,i), 0 ) & 1301 ) 1302 ENDDO 1303 1304 ELSE 1305 1306 DO k = nzb+1, nzt 1307 ! 1308 !-- Flag 9 is used to mask top fluxes, flag 30 to mask 1309 !-- surface fluxes 1310 tend(k,j,i) = tend(k,j,i) - & 1311 kh(k,j,i) * g / pt(k,j,i) * & 1312 ( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k) * & 1313 MERGE( 1.0_wp, 0.0_wp, & 1314 BTEST( wall_flags_0(k,j,i), 30 ) & 1315 ) * & 1316 MERGE( 1.0_wp, 0.0_wp, & 1317 BTEST( wall_flags_0(k,j,i), 9 ) & 1318 ) 1319 ENDDO 1320 1321 IF ( use_surface_fluxes ) THEN 1322 ! 1323 !-- Default surfaces, up- and downward-facing 1324 DO l = 0, 1 1325 surf_s = surf_def_h(l)%start_index(j,i) 1326 surf_e = surf_def_h(l)%end_index(j,i) 1327 DO m = surf_s, surf_e 1328 k = surf_def_h(l)%k(m) 1329 tend(k,j,i) = tend(k,j,i) + g / pt_reference & 1330 * surf_def_h(l)%shf(m) 1331 ENDDO 1332 ENDDO 1333 ! 1334 !-- Natural surfaces 1335 surf_s = surf_lsm_h%start_index(j,i) 1336 surf_e = surf_lsm_h%end_index(j,i) 1337 DO m = surf_s, surf_e 1338 k = surf_lsm_h%k(m) 1339 tend(k,j,i) = tend(k,j,i) + g / pt_reference & 1340 * surf_lsm_h%shf(m) 1341 ENDDO 1342 ! 1343 !-- Urban surfaces 1344 surf_s = surf_usm_h%start_index(j,i) 1345 surf_e = surf_usm_h%end_index(j,i) 1346 DO m = surf_s, surf_e 1347 k = surf_usm_h%k(m) 1348 tend(k,j,i) = tend(k,j,i) + g / pt_reference & 1349 * surf_usm_h%shf(m) 1350 ENDDO 1351 ENDIF 1352 1353 IF ( use_top_fluxes ) THEN 1354 surf_s = surf_def_h(2)%start_index(j,i) 1355 surf_e = surf_def_h(2)%end_index(j,i) 1356 DO m = surf_s, surf_e 1357 k = surf_def_h(2)%k(m) 1358 tend(k,j,i) = tend(k,j,i) + g / pt_reference * & 1359 surf_def_h(2)%shf(m) 1360 ENDDO 1361 ENDIF 1362 1363 ENDIF 1364 1365 ENDIF 1366 1367 ELSE 1368 1369 DO k = nzb+1, nzt 1370 ! 1371 !-- Flag 9 is used to mask top fluxes, flag 30 to mask surface fluxes 1372 IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN 1373 k1 = 1.0_wp + 0.61_wp * q(k,j,i) 1374 k2 = 0.61_wp * pt(k,j,i) 1375 tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / vpt(k,j,i) * & 1376 ( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) + & 1377 k2 * ( q(k+1,j,i) - q(k-1,j,i) ) & 1378 ) * dd2zu(k) * & 1379 MERGE( 1.0_wp, 0.0_wp, & 1380 BTEST( wall_flags_0(k,j,i), 30 ) & 1381 ) * & 1382 MERGE( 1.0_wp, 0.0_wp, & 1383 BTEST( wall_flags_0(k,j,i), 9 ) & 1384 ) 1385 1386 ELSE IF ( cloud_physics ) THEN 1387 IF ( ql(k,j,i) == 0.0_wp ) THEN 1388 k1 = 1.0_wp + 0.61_wp * q(k,j,i) 1389 k2 = 0.61_wp * pt(k,j,i) 1390 ELSE 1391 theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i) 1392 temp = theta * t_d_pt(k) 1393 k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * & 1394 ( q(k,j,i) - ql(k,j,i) ) * & 1395 ( 1.0_wp + 0.622_wp * l_d_r / temp ) ) / & 1396 ( 1.0_wp + 0.622_wp * l_d_r * l_d_cp * & 1397 ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) ) 1398 k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp ) 1399 ENDIF 1400 tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / vpt(k,j,i) * & 1401 ( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) + & 1402 k2 * ( q(k+1,j,i) - q(k-1,j,i) ) & 1403 ) * dd2zu(k) * & 1404 MERGE( 1.0_wp, 0.0_wp, & 1405 BTEST( wall_flags_0(k,j,i), 30 ) & 1406 ) * & 1407 MERGE( 1.0_wp, 0.0_wp, & 1408 BTEST( wall_flags_0(k,j,i), 9 ) & 1409 ) 1410 ELSE IF ( cloud_droplets ) THEN 1411 k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i) 1412 k2 = 0.61_wp * pt(k,j,i) 1413 tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / vpt(k,j,i) * & 1414 ( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) + & 1415 k2 * ( q(k+1,j,i) - q(k-1,j,i) ) - & 1416 pt(k,j,i) * ( ql(k+1,j,i) - & 1417 ql(k-1,j,i) ) ) * dd2zu(k)& 1418 * MERGE( 1.0_wp, 0.0_wp, & 1419 BTEST( wall_flags_0(k,j,i), 30 ) & 1420 ) & 1421 * MERGE( 1.0_wp, 0.0_wp, & 1422 BTEST( wall_flags_0(k,j,i), 9 ) & 1423 ) 1424 ENDIF 1425 ENDDO 1426 1427 IF ( use_surface_fluxes ) THEN 1428 ! 1429 !-- Treat horizontal default surfaces, up- and downward-facing 1430 DO l = 0, 1 1431 surf_s = surf_def_h(l)%start_index(j,i) 1432 surf_e = surf_def_h(l)%end_index(j,i) 1433 DO m = surf_s, surf_e 1434 k = surf_def_h(l)%k(m) 703 1435 704 1436 IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN … … 710 1442 k2 = 0.61_wp * pt(k,j,i) 711 1443 ELSE 712 713 714 k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp *&715 ( q(k,j,i) - ql(k,j,i) ) *&716 ( 1.0_wp + 0.622_wp * l_d_r / temp ) ) /&717 ( 1.0_wp + 0.622_wp * l_d_r * l_d_cp * &718 719 1444 theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i) 1445 temp = theta * t_d_pt(k) 1446 k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * & 1447 ( q(k,j,i) - ql(k,j,i) ) * & 1448 ( 1.0_wp + 0.622_wp * l_d_r / temp ) ) / & 1449 ( 1.0_wp + 0.622_wp * l_d_r * l_d_cp * & 1450 ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) ) 1451 k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp ) 720 1452 ENDIF 721 1453 ELSE IF ( cloud_droplets ) THEN … … 724 1456 ENDIF 725 1457 726 tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) * & 727 ( k1* tswst(j,i) + k2 * qswst(j,i) ) 1458 tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) * & 1459 ( k1 * surf_def_h(l)%shf(m) + & 1460 k2 * surf_def_h(l)%qsws(m) ) 728 1461 ENDDO 729 730 ENDIF 1462 ENDDO 1463 ! 1464 !-- Treat horizontal natural surfaces 1465 surf_s = surf_lsm_h%start_index(j,i) 1466 surf_e = surf_lsm_h%end_index(j,i) 1467 DO m = surf_s, surf_e 1468 k = surf_lsm_h%k(m) 1469 1470 IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN 1471 k1 = 1.0_wp + 0.61_wp * q(k,j,i) 1472 k2 = 0.61_wp * pt(k,j,i) 1473 ELSE IF ( cloud_physics ) THEN 1474 IF ( ql(k,j,i) == 0.0_wp ) THEN 1475 k1 = 1.0_wp + 0.61_wp * q(k,j,i) 1476 k2 = 0.61_wp * pt(k,j,i) 1477 ELSE 1478 theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i) 1479 temp = theta * t_d_pt(k) 1480 k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * & 1481 ( q(k,j,i) - ql(k,j,i) ) * & 1482 ( 1.0_wp + 0.622_wp * l_d_r / temp ) ) / & 1483 ( 1.0_wp + 0.622_wp * l_d_r * l_d_cp * & 1484 ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) ) 1485 k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp ) 1486 ENDIF 1487 ELSE IF ( cloud_droplets ) THEN 1488 k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i) 1489 k2 = 0.61_wp * pt(k,j,i) 1490 ENDIF 1491 1492 tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) * & 1493 ( k1 * surf_lsm_h%shf(m) + & 1494 k2 * surf_lsm_h%qsws(m) ) 1495 ENDDO 1496 ! 1497 !-- Treat horizontal urban surfaces 1498 surf_s = surf_usm_h%start_index(j,i) 1499 surf_e = surf_usm_h%end_index(j,i) 1500 DO m = surf_s, surf_e 1501 k = surf_usm_h%k(m) 1502 1503 IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN 1504 k1 = 1.0_wp + 0.61_wp * q(k,j,i) 1505 k2 = 0.61_wp * pt(k,j,i) 1506 ELSE IF ( cloud_physics ) THEN 1507 IF ( ql(k,j,i) == 0.0_wp ) THEN 1508 k1 = 1.0_wp + 0.61_wp * q(k,j,i) 1509 k2 = 0.61_wp * pt(k,j,i) 1510 ELSE 1511 theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i) 1512 temp = theta * t_d_pt(k) 1513 k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * & 1514 ( q(k,j,i) - ql(k,j,i) ) * & 1515 ( 1.0_wp + 0.622_wp * l_d_r / temp ) ) / & 1516 ( 1.0_wp + 0.622_wp * l_d_r * l_d_cp * & 1517 ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) ) 1518 k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp ) 1519 ENDIF 1520 ELSE IF ( cloud_droplets ) THEN 1521 k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i) 1522 k2 = 0.61_wp * pt(k,j,i) 1523 ENDIF 1524 1525 tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) * & 1526 ( k1 * surf_usm_h%shf(m) + & 1527 k2 * surf_usm_h%qsws(m) ) 1528 ENDDO 731 1529 732 1530 ENDIF 733 1531 734 ENDIF 735 736 ENDDO 737 738 END SUBROUTINE production_e 739 740 741 !------------------------------------------------------------------------------! 742 ! Description: 743 ! ------------ 744 !> Call for grid point i,j 745 !------------------------------------------------------------------------------! 746 SUBROUTINE production_e_ij( i, j ) 747 748 USE arrays_3d, & 749 ONLY: ddzw, dd2zu, kh, km, prho, pt, q, ql, qsws, qswst, shf, & 750 tend, tswst, u, v, vpt, w 751 752 USE cloud_parameters, & 753 ONLY: l_d_cp, l_d_r, pt_d_t, t_d_pt 754 755 USE control_parameters, & 756 ONLY: cloud_droplets, cloud_physics, constant_flux_layer, g, & 757 humidity, kappa, neutral, ocean, pt_reference, & 758 rho_reference, use_single_reference_value, & 759 use_surface_fluxes, use_top_fluxes 760 761 USE grid_variables, & 762 ONLY: ddx, dx, ddy, dy, wall_e_x, wall_e_y 763 764 USE indices, & 765 ONLY: nxl, nxr, nys, nyn, nzb, nzb_diff_s_inner, & 766 nzb_diff_s_outer, nzb_s_inner, nzt, nzt_diff 767 768 IMPLICIT NONE 769 770 INTEGER(iwp) :: i !< 771 INTEGER(iwp) :: j !< 772 INTEGER(iwp) :: k !< 773 774 REAL(wp) :: def !< 775 REAL(wp) :: dudx !< 776 REAL(wp) :: dudy !< 777 REAL(wp) :: dudz !< 778 REAL(wp) :: dvdx !< 779 REAL(wp) :: dvdy !< 780 REAL(wp) :: dvdz !< 781 REAL(wp) :: dwdx !< 782 REAL(wp) :: dwdy !< 783 REAL(wp) :: dwdz !< 784 REAL(wp) :: k1 !< 785 REAL(wp) :: k2 !< 786 REAL(wp) :: km_neutral !< 787 REAL(wp) :: theta !< 788 REAL(wp) :: temp !< 789 790 REAL(wp), DIMENSION(nzb:nzt+1) :: usvs !< 791 REAL(wp), DIMENSION(nzb:nzt+1) :: vsus !< 792 REAL(wp), DIMENSION(nzb:nzt+1) :: wsus !< 793 REAL(wp), DIMENSION(nzb:nzt+1) :: wsvs !< 794 795 ! 796 !-- Calculate TKE production by shear 797 DO k = nzb_diff_s_outer(j,i), nzt 798 799 dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx 800 dudy = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - & 801 u(k,j-1,i) - u(k,j-1,i+1) ) * ddy 802 dudz = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - & 803 u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k) 804 805 dvdx = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - & 806 v(k,j,i-1) - v(k,j+1,i-1) ) * ddx 807 dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy 808 dvdz = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - & 809 v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k) 810 811 dwdx = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - & 812 w(k,j,i-1) - w(k-1,j,i-1) ) * ddx 813 dwdy = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - & 814 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy 815 dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) 816 817 def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 ) & 818 + dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + dvdz**2 & 819 + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz ) 820 821 IF ( def < 0.0_wp ) def = 0.0_wp 822 823 tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def 824 825 ENDDO 826 827 IF ( constant_flux_layer ) THEN 828 829 IF ( ( wall_e_x(j,i) /= 0.0_wp ) .OR. ( wall_e_y(j,i) /= 0.0_wp ) ) THEN 830 831 ! 832 !-- Position beneath wall 833 !-- (2) - Will allways be executed. 834 !-- 'bottom and wall: use u_0,v_0 and wall functions' 835 k = nzb_diff_s_inner(j,i)-1 836 837 dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx 838 dudz = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - & 839 u_0(j,i) - u_0(j,i+1) ) * dd2zu(k) 840 dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy 841 dvdz = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - & 842 v_0(j,i) - v_0(j+1,i) ) * dd2zu(k) 843 dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) 844 845 IF ( wall_e_y(j,i) /= 0.0_wp ) THEN 846 ! 847 !-- Inconsistency removed: as the thermal stratification 848 !-- is not taken into account for the evaluation of the 849 !-- wall fluxes at vertical walls, the eddy viscosity km 850 !-- must not be used for the evaluation of the velocity 851 !-- gradients dudy and dwdy 852 !-- Note: The validity of the new method has not yet 853 !-- been shown, as so far no suitable data for a 854 !-- validation has been available 855 CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, & 856 usvs, 1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp ) 857 CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, & 858 wsvs, 0.0_wp, 0.0_wp, 1.0_wp, 0.0_wp ) 859 km_neutral = kappa * ( usvs(k)**2 + wsvs(k)**2 )**0.25_wp * & 860 0.5_wp * dy 861 IF ( km_neutral > 0.0_wp ) THEN 862 dudy = - wall_e_y(j,i) * usvs(k) / km_neutral 863 dwdy = - wall_e_y(j,i) * wsvs(k) / km_neutral 864 ELSE 865 dudy = 0.0_wp 866 dwdy = 0.0_wp 867 ENDIF 868 ELSE 869 dudy = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - & 870 u(k,j-1,i) - u(k,j-1,i+1) ) * ddy 871 dwdy = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - & 872 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy 873 ENDIF 874 875 IF ( wall_e_x(j,i) /= 0.0_wp ) THEN 876 ! 877 !-- Inconsistency removed: as the thermal stratification 878 !-- is not taken into account for the evaluation of the 879 !-- wall fluxes at vertical walls, the eddy viscosity km 880 !-- must not be used for the evaluation of the velocity 881 !-- gradients dvdx and dwdx 882 !-- Note: The validity of the new method has not yet 883 !-- been shown, as so far no suitable data for a 884 !-- validation has been available 885 CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, & 886 vsus, 0.0_wp, 1.0_wp, 0.0_wp, 0.0_wp ) 887 CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, & 888 wsus, 0.0_wp, 0.0_wp, 0.0_wp, 1.0_wp ) 889 km_neutral = kappa * ( vsus(k)**2 + wsus(k)**2 )**0.25_wp * & 890 0.5_wp * dx 891 IF ( km_neutral > 0.0_wp ) THEN 892 dvdx = - wall_e_x(j,i) * vsus(k) / km_neutral 893 dwdx = - wall_e_x(j,i) * wsus(k) / km_neutral 894 ELSE 895 dvdx = 0.0_wp 896 dwdx = 0.0_wp 897 ENDIF 898 ELSE 899 dvdx = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - & 900 v(k,j,i-1) - v(k,j+1,i-1) ) * ddx 901 dwdx = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - & 902 w(k,j,i-1) - w(k-1,j,i-1) ) * ddx 903 ENDIF 904 905 def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 ) + & 906 dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + & 907 dvdz**2 + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz ) 908 909 IF ( def < 0.0_wp ) def = 0.0_wp 910 911 tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def 912 913 ! 914 !-- (3) - will be executed only, if there is at least one level 915 !-- between (2) and (4), i.e. the topography must have a 916 !-- minimum height of 2 dz. Wall fluxes for this case have 917 !-- already been calculated for (2). 918 !-- 'wall only: use wall functions' 919 DO k = nzb_diff_s_inner(j,i), nzb_diff_s_outer(j,i)-2 920 921 dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx 922 dudz = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - & 923 u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k) 924 dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy 925 dvdz = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - & 926 v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k) 927 dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) 928 929 IF ( wall_e_y(j,i) /= 0.0_wp ) THEN 930 ! 931 !-- Inconsistency removed: as the thermal stratification 932 !-- is not taken into account for the evaluation of the 933 !-- wall fluxes at vertical walls, the eddy viscosity km 934 !-- must not be used for the evaluation of the velocity 935 !-- gradients dudy and dwdy 936 !-- Note: The validity of the new method has not yet 937 !-- been shown, as so far no suitable data for a 938 !-- validation has been available 939 km_neutral = kappa * ( usvs(k)**2 + & 940 wsvs(k)**2 )**0.25_wp * 0.5_wp * dy 941 IF ( km_neutral > 0.0_wp ) THEN 942 dudy = - wall_e_y(j,i) * usvs(k) / km_neutral 943 dwdy = - wall_e_y(j,i) * wsvs(k) / km_neutral 944 ELSE 945 dudy = 0.0_wp 946 dwdy = 0.0_wp 947 ENDIF 948 ELSE 949 dudy = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - & 950 u(k,j-1,i) - u(k,j-1,i+1) ) * ddy 951 dwdy = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - & 952 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy 953 ENDIF 954 955 IF ( wall_e_x(j,i) /= 0.0_wp ) THEN 956 ! 957 !-- Inconsistency removed: as the thermal stratification 958 !-- is not taken into account for the evaluation of the 959 !-- wall fluxes at vertical walls, the eddy viscosity km 960 !-- must not be used for the evaluation of the velocity 961 !-- gradients dvdx and dwdx 962 !-- Note: The validity of the new method has not yet 963 !-- been shown, as so far no suitable data for a 964 !-- validation has been available 965 km_neutral = kappa * ( vsus(k)**2 + & 966 wsus(k)**2 )**0.25_wp * 0.5_wp * dx 967 IF ( km_neutral > 0.0_wp ) THEN 968 dvdx = - wall_e_x(j,i) * vsus(k) / km_neutral 969 dwdx = - wall_e_x(j,i) * wsus(k) / km_neutral 970 ELSE 971 dvdx = 0.0_wp 972 dwdx = 0.0_wp 973 ENDIF 974 ELSE 975 dvdx = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - & 976 v(k,j,i-1) - v(k,j+1,i-1) ) * ddx 977 dwdx = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - & 978 w(k,j,i-1) - w(k-1,j,i-1) ) * ddx 979 ENDIF 980 981 def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 ) + & 982 dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + & 983 dvdz**2 + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz ) 984 985 IF ( def < 0.0_wp ) def = 0.0_wp 986 987 tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def 988 989 ENDDO 990 991 ! 992 !-- (4) - will allways be executed. 993 !-- 'special case: free atmosphere' (as for case (0)) 994 k = nzb_diff_s_outer(j,i)-1 995 996 dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx 997 dudy = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - & 998 u(k,j-1,i) - u(k,j-1,i+1) ) * ddy 999 dudz = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - & 1000 u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k) 1001 1002 dvdx = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - & 1003 v(k,j,i-1) - v(k,j+1,i-1) ) * ddx 1004 dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy 1005 dvdz = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - & 1006 v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k) 1007 1008 dwdx = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - & 1009 w(k,j,i-1) - w(k-1,j,i-1) ) * ddx 1010 dwdy = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - & 1011 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy 1012 dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) 1013 1014 def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 ) + & 1015 dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + & 1016 dvdz**2 + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz ) 1017 1018 IF ( def < 0.0_wp ) def = 0.0_wp 1019 1020 tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def 1021 1022 ELSE 1023 1024 ! 1025 !-- Position without adjacent wall 1026 !-- (1) - will allways be executed. 1027 !-- 'bottom only: use u_0,v_0' 1028 k = nzb_diff_s_inner(j,i)-1 1029 1030 dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx 1031 dudy = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - & 1032 u(k,j-1,i) - u(k,j-1,i+1) ) * ddy 1033 dudz = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - & 1034 u_0(j,i) - u_0(j,i+1) ) * dd2zu(k) 1035 1036 dvdx = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - & 1037 v(k,j,i-1) - v(k,j+1,i-1) ) * ddx 1038 dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy 1039 dvdz = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - & 1040 v_0(j,i) - v_0(j+1,i) ) * dd2zu(k) 1041 1042 dwdx = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - & 1043 w(k,j,i-1) - w(k-1,j,i-1) ) * ddx 1044 dwdy = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - & 1045 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy 1046 dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) 1047 1048 def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 ) & 1049 + dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + dvdz**2 & 1050 + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz ) 1051 1052 IF ( def < 0.0_wp ) def = 0.0_wp 1053 1054 tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def 1055 1056 ENDIF 1057 1058 ELSEIF ( use_surface_fluxes ) THEN 1059 1060 k = nzb_diff_s_outer(j,i)-1 1061 1062 dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx 1063 dudy = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - & 1064 u(k,j-1,i) - u(k,j-1,i+1) ) * ddy 1065 dudz = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - & 1066 u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k) 1067 1068 dvdx = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - & 1069 v(k,j,i-1) - v(k,j+1,i-1) ) * ddx 1070 dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy 1071 dvdz = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - & 1072 v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k) 1073 1074 dwdx = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - & 1075 w(k,j,i-1) - w(k-1,j,i-1) ) * ddx 1076 dwdy = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - & 1077 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy 1078 dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) 1079 1080 def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 ) + & 1081 dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + & 1082 dvdz**2 + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz ) 1083 1084 IF ( def < 0.0_wp ) def = 0.0_wp 1085 1086 tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def 1087 1088 ENDIF 1089 1090 ! 1091 !-- If required, calculate TKE production by buoyancy 1092 IF ( .NOT. neutral ) THEN 1093 1094 IF ( .NOT. humidity ) THEN 1095 1096 IF ( use_single_reference_value ) THEN 1097 1098 IF ( ocean ) THEN 1099 ! 1100 !-- So far in the ocean no special treatment of density flux in 1101 !-- the bottom and top surface layer 1102 DO k = nzb_s_inner(j,i)+1, nzt 1103 tend(k,j,i) = tend(k,j,i) + & 1104 kh(k,j,i) * g / rho_reference * & 1105 ( prho(k+1,j,i) - prho(k-1,j,i) ) * dd2zu(k) 1106 ENDDO 1107 1108 ELSE 1109 1110 DO k = nzb_diff_s_inner(j,i), nzt_diff 1111 tend(k,j,i) = tend(k,j,i) - & 1112 kh(k,j,i) * g / pt_reference * & 1113 ( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k) 1114 ENDDO 1115 1116 IF ( use_surface_fluxes ) THEN 1117 k = nzb_diff_s_inner(j,i)-1 1118 tend(k,j,i) = tend(k,j,i) + g / pt_reference * shf(j,i) 1119 ENDIF 1120 1121 IF ( use_top_fluxes ) THEN 1122 k = nzt 1123 tend(k,j,i) = tend(k,j,i) + g / pt_reference * tswst(j,i) 1124 ENDIF 1125 1126 ENDIF 1127 1128 ELSE 1129 1130 IF ( ocean ) THEN 1131 ! 1132 !-- So far in the ocean no special treatment of density flux in 1133 !-- the bottom and top surface layer 1134 DO k = nzb_s_inner(j,i)+1, nzt 1135 tend(k,j,i) = tend(k,j,i) + & 1136 kh(k,j,i) * g / prho(k,j,i) * & 1137 ( prho(k+1,j,i) - prho(k-1,j,i) ) * dd2zu(k) 1138 ENDDO 1139 1140 ELSE 1141 1142 DO k = nzb_diff_s_inner(j,i), nzt_diff 1143 tend(k,j,i) = tend(k,j,i) - & 1144 kh(k,j,i) * g / pt(k,j,i) * & 1145 ( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k) 1146 ENDDO 1147 1148 IF ( use_surface_fluxes ) THEN 1149 k = nzb_diff_s_inner(j,i)-1 1150 tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * shf(j,i) 1151 ENDIF 1152 1153 IF ( use_top_fluxes ) THEN 1154 k = nzt 1155 tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * tswst(j,i) 1156 ENDIF 1157 1158 ENDIF 1159 1160 ENDIF 1161 1162 ELSE 1163 1164 DO k = nzb_diff_s_inner(j,i), nzt_diff 1165 1166 IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN 1167 k1 = 1.0_wp + 0.61_wp * q(k,j,i) 1168 k2 = 0.61_wp * pt(k,j,i) 1169 tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / vpt(k,j,i) * & 1170 ( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) + & 1171 k2 * ( q(k+1,j,i) - q(k-1,j,i) ) & 1172 ) * dd2zu(k) 1173 ELSE IF ( cloud_physics ) THEN 1174 IF ( ql(k,j,i) == 0.0_wp ) THEN 1532 IF ( use_top_fluxes ) THEN 1533 surf_s = surf_def_h(2)%start_index(j,i) 1534 surf_e = surf_def_h(2)%end_index(j,i) 1535 DO m = surf_s, surf_e 1536 k = surf_def_h(2)%k(m) 1537 1538 1539 1540 IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN 1175 1541 k1 = 1.0_wp + 0.61_wp * q(k,j,i) 1176 1542 k2 = 0.61_wp * pt(k,j,i) 1177 ELSE 1178 theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i) 1179 temp = theta * t_d_pt(k) 1180 k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * & 1181 ( q(k,j,i) - ql(k,j,i) ) * & 1182 ( 1.0_wp + 0.622_wp * l_d_r / temp ) ) / & 1183 ( 1.0_wp + 0.622_wp * l_d_r * l_d_cp * & 1543 ELSE IF ( cloud_physics ) THEN 1544 IF ( ql(k,j,i) == 0.0_wp ) THEN 1545 k1 = 1.0_wp + 0.61_wp * q(k,j,i) 1546 k2 = 0.61_wp * pt(k,j,i) 1547 ELSE 1548 theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i) 1549 temp = theta * t_d_pt(k) 1550 k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * & 1551 ( q(k,j,i) - ql(k,j,i) ) * & 1552 ( 1.0_wp + 0.622_wp * l_d_r / temp ) ) / & 1553 ( 1.0_wp + 0.622_wp * l_d_r * l_d_cp * & 1184 1554 ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) ) 1185 k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp ) 1555 k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp ) 1556 ENDIF 1557 ELSE IF ( cloud_droplets ) THEN 1558 k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i) 1559 k2 = 0.61_wp * pt(k,j,i) 1186 1560 ENDIF 1187 tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / vpt(k,j,i) * & 1188 ( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) + & 1189 k2 * ( q(k+1,j,i) - q(k-1,j,i) ) & 1190 ) * dd2zu(k) 1191 ELSE IF ( cloud_droplets ) THEN 1192 k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i) 1193 k2 = 0.61_wp * pt(k,j,i) 1194 tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / vpt(k,j,i) * & 1195 ( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) + & 1196 k2 * ( q(k+1,j,i) - q(k-1,j,i) ) - & 1197 pt(k,j,i) * ( ql(k+1,j,i) - & 1198 ql(k-1,j,i) ) ) * dd2zu(k) 1199 ENDIF 1200 ENDDO 1201 1202 IF ( use_surface_fluxes ) THEN 1203 k = nzb_diff_s_inner(j,i)-1 1204 1205 IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN 1206 k1 = 1.0_wp + 0.61_wp * q(k,j,i) 1207 k2 = 0.61_wp * pt(k,j,i) 1208 ELSE IF ( cloud_physics ) THEN 1209 IF ( ql(k,j,i) == 0.0_wp ) THEN 1210 k1 = 1.0_wp + 0.61_wp * q(k,j,i) 1211 k2 = 0.61_wp * pt(k,j,i) 1212 ELSE 1213 theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i) 1214 temp = theta * t_d_pt(k) 1215 k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * & 1216 ( q(k,j,i) - ql(k,j,i) ) * & 1217 ( 1.0_wp + 0.622_wp * l_d_r / temp ) ) / & 1218 ( 1.0_wp + 0.622_wp * l_d_r * l_d_cp * & 1219 ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) ) 1220 k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp ) 1221 ENDIF 1222 ELSE IF ( cloud_droplets ) THEN 1223 k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i) 1224 k2 = 0.61_wp * pt(k,j,i) 1225 ENDIF 1226 1227 tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) * & 1228 ( k1* shf(j,i) + k2 * qsws(j,i) ) 1229 ENDIF 1230 1231 IF ( use_top_fluxes ) THEN 1232 k = nzt 1233 1234 IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN 1235 k1 = 1.0_wp + 0.61_wp * q(k,j,i) 1236 k2 = 0.61_wp * pt(k,j,i) 1237 ELSE IF ( cloud_physics ) THEN 1238 IF ( ql(k,j,i) == 0.0_wp ) THEN 1239 k1 = 1.0_wp + 0.61_wp * q(k,j,i) 1240 k2 = 0.61_wp * pt(k,j,i) 1241 ELSE 1242 theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i) 1243 temp = theta * t_d_pt(k) 1244 k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * & 1245 ( q(k,j,i) - ql(k,j,i) ) * & 1246 ( 1.0_wp + 0.622_wp * l_d_r / temp ) ) / & 1247 ( 1.0_wp + 0.622_wp * l_d_r * l_d_cp * & 1248 ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) ) 1249 k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp ) 1250 ENDIF 1251 ELSE IF ( cloud_droplets ) THEN 1252 k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i) 1253 k2 = 0.61_wp * pt(k,j,i) 1254 ENDIF 1255 1256 tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) * & 1257 ( k1* tswst(j,i) + k2 * qswst(j,i) ) 1561 1562 tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) * & 1563 ( k1* surf_def_h(2)%shf(m) + & 1564 k2 * surf_def_h(2)%qsws(m) ) 1565 ENDDO 1566 1258 1567 ENDIF 1259 1568 … … 1273 1582 1274 1583 USE arrays_3d, & 1275 ONLY: kh, km, u, us, usws, v, vsws, zu1584 ONLY: kh, km, u, v, zu 1276 1585 1277 1586 USE control_parameters, & … … 1282 1591 nzb_v_inner 1283 1592 1593 USE surface_mod, & 1594 ONLY : surf_def_h, surf_def_v, surf_lsm_h, surf_usm_h 1595 1284 1596 IMPLICIT NONE 1285 1597 1286 INTEGER(iwp) :: i !< 1287 INTEGER(iwp) :: j !< 1288 INTEGER(iwp) :: ku !< 1289 INTEGER(iwp) :: kv !< 1598 INTEGER(iwp) :: i !< grid index x-direction 1599 INTEGER(iwp) :: j !< grid index y-direction 1600 INTEGER(iwp) :: k !< grid index z-direction 1601 INTEGER(iwp) :: l !< running index surface type (up- or downward-facing) 1602 INTEGER(iwp) :: m !< running index surface elements 1290 1603 1291 1604 IF ( constant_flux_layer ) THEN 1292 1293 IF ( first_call ) THEN1294 ALLOCATE( u_0(nysg:nyng,nxlg:nxrg), v_0(nysg:nyng,nxlg:nxrg) )1295 u_0 = 0.0_wp ! just to avoid access of uninitialized memory1296 v_0 = 0.0_wp ! within exchange_horiz_2d1297 first_call = .FALSE.1298 ENDIF1299 1300 1605 ! 1301 1606 !-- Calculate a virtual velocity at the surface in a way that the … … 1306 1611 !-- (otherwise the timestep may be significantly reduced by large 1307 1612 !-- surface winds). 1308 !-- Upper bounds are nxr+1 and nyn+1 because otherwise these values are1309 1613 !-- not available in case of non-cyclic boundary conditions. 1310 1614 !-- WARNING: the exact analytical solution would require the determination 1311 1615 !-- of the eddy diffusivity by km = u* * kappa * zp / phi_m. 1312 !$OMP PARALLEL DO PRIVATE( ku, kv ) 1313 DO i = nxl, nxr+1 1314 DO j = nys, nyn+1 1315 1316 ku = nzb_u_inner(j,i)+1 1317 kv = nzb_v_inner(j,i)+1 1318 1319 u_0(j,i) = u(ku+1,j,i) + usws(j,i) * ( zu(ku+1) - zu(ku-1) ) / & 1320 ( 0.5_wp * ( km(ku,j,i) + km(ku,j,i-1) ) + & 1321 1.0E-20_wp ) 1322 ! ( us(j,i) * kappa * zu(1) ) 1323 v_0(j,i) = v(kv+1,j,i) + vsws(j,i) * ( zu(kv+1) - zu(kv-1) ) / & 1324 ( 0.5_wp * ( km(kv,j,i) + km(kv,j-1,i) ) + & 1325 1.0E-20_wp ) 1326 ! ( us(j,i) * kappa * zu(1) ) 1327 1328 IF ( ABS( u(ku+1,j,i) - u_0(j,i) ) > & 1329 ABS( u(ku+1,j,i) - u(ku-1,j,i) ) ) u_0(j,i) = u(ku-1,j,i) 1330 IF ( ABS( v(kv+1,j,i) - v_0(j,i) ) > & 1331 ABS( v(kv+1,j,i) - v(kv-1,j,i) ) ) v_0(j,i) = v(kv-1,j,i) 1332 1333 ENDDO 1616 !-- Default surfaces, upward-facing 1617 !$OMP PARALLEL DO PRIVATE(i,j,k,m) 1618 DO m = 1, surf_def_h(0)%ns 1619 1620 i = surf_def_h(0)%i(m) 1621 j = surf_def_h(0)%j(m) 1622 k = surf_def_h(0)%k(m) 1623 ! 1624 !-- Note, calculatione of u_0 and v_0 is not fully accurate, as u/v 1625 !-- and km are not on the same grid. Actually, a further 1626 !-- interpolation of km onto the u/v-grid is necessary. However, the 1627 !-- effect of this error is negligible. 1628 surf_def_h(0)%u_0(m) = u(k+1,j,i) + surf_def_h(0)%usws(m) * & 1629 ( zu(k+1) - zu(k-1) ) / & 1630 ( km(k,j,i) + 1.0E-20_wp ) 1631 surf_def_h(0)%v_0(m) = v(k+1,j,i) + surf_def_h(0)%vsws(m) * & 1632 ( zu(k+1) - zu(k-1) ) / & 1633 ( km(k,j,i) + 1.0E-20_wp ) 1634 1635 IF ( ABS( u(k+1,j,i) - surf_def_h(0)%u_0(m) ) > & 1636 ABS( u(k+1,j,i) - u(k-1,j,i) ) & 1637 ) surf_def_h(0)%u_0(m) = u(k-1,j,i) 1638 1639 IF ( ABS( v(k+1,j,i) - surf_def_h(0)%v_0(m) ) > & 1640 ABS( v(k+1,j,i) - v(k-1,j,i) ) & 1641 ) surf_def_h(0)%v_0(m) = v(k-1,j,i) 1642 1334 1643 ENDDO 1335 1336 CALL exchange_horiz_2d( u_0 ) 1337 CALL exchange_horiz_2d( v_0 ) 1644 ! 1645 !-- Default surfaces, downward-facing 1646 !$OMP PARALLEL DO PRIVATE(i,j,k,m) 1647 DO m = 1, surf_def_h(1)%ns 1648 1649 i = surf_def_h(1)%i(m) 1650 j = surf_def_h(1)%j(m) 1651 k = surf_def_h(1)%k(m) 1652 ! 1653 !-- Note, calculatione of u_0 and v_0 is not fully accurate, as u/v 1654 !-- and km are not on the same grid. Actually, a further 1655 !-- interpolation of km onto the u/v-grid is necessary. However, the 1656 !-- effect of this error is negligible. 1657 !-- In case of downward-facing surfaces, gradient is calculated 1658 !-- between u_0 and u(k-1). 1659 surf_def_h(1)%u_0(m) = u(k-1,j,i) - surf_def_h(1)%usws(m) * & 1660 ( zu(k+1) - zu(k-1) ) / & 1661 ( km(k,j,i) + 1.0E-20_wp ) 1662 surf_def_h(1)%v_0(m) = v(k-1,j,i) - surf_def_h(1)%vsws(m) * & 1663 ( zu(k+1) - zu(k-1) ) / & 1664 ( km(k,j,i) + 1.0E-20_wp ) 1665 1666 IF ( ABS( surf_def_h(1)%u_0(m) - u(k-1,j,i) ) > & 1667 ABS( u(k+1,j,i) - u(k-1,j,i) ) & 1668 ) surf_def_h(1)%u_0(m) = u(k+1,j,i) 1669 1670 IF ( ABS( surf_def_h(1)%v_0(m) - v(k-1,j,i) ) > & 1671 ABS( v(k+1,j,i) - v(k-1,j,i) ) & 1672 ) surf_def_h(1)%v_0(m) = v(k+1,j,i) 1673 1674 ENDDO 1675 ! 1676 !-- Natural surfaces, upward-facing 1677 !$OMP PARALLEL DO PRIVATE(i,j,k,m) 1678 DO m = 1, surf_lsm_h%ns 1679 1680 i = surf_lsm_h%i(m) 1681 j = surf_lsm_h%j(m) 1682 k = surf_lsm_h%k(m) 1683 ! 1684 !-- Note, calculatione of u_0 and v_0 is not fully accurate, as u/v 1685 !-- and km are not on the same grid. Actually, a further 1686 !-- interpolation of km onto the u/v-grid is necessary. However, the 1687 !-- effect of this error is negligible. 1688 surf_lsm_h%u_0(m) = u(k+1,j,i) + surf_lsm_h%usws(m) * & 1689 ( zu(k+1) - zu(k-1) ) / & 1690 ( km(k,j,i) + 1.0E-20_wp ) 1691 surf_lsm_h%v_0(m) = v(k+1,j,i) + surf_lsm_h%vsws(m) * & 1692 ( zu(k+1) - zu(k-1) ) / & 1693 ( km(k,j,i) + 1.0E-20_wp ) 1694 1695 IF ( ABS( u(k+1,j,i) - surf_lsm_h%u_0(m) ) > & 1696 ABS( u(k+1,j,i) - u(k-1,j,i) ) & 1697 ) surf_lsm_h%u_0(m) = u(k-1,j,i) 1698 1699 IF ( ABS( v(k+1,j,i) - surf_lsm_h%v_0(m) ) > & 1700 ABS( v(k+1,j,i) - v(k-1,j,i) ) & 1701 ) surf_lsm_h%v_0(m) = v(k-1,j,i) 1702 1703 ENDDO 1704 ! 1705 !-- Urban surfaces, upward-facing 1706 !$OMP PARALLEL DO PRIVATE(i,j,k,m) 1707 DO m = 1, surf_usm_h%ns 1708 1709 i = surf_usm_h%i(m) 1710 j = surf_usm_h%j(m) 1711 k = surf_usm_h%k(m) 1712 ! 1713 !-- Note, calculatione of u_0 and v_0 is not fully accurate, as u/v 1714 !-- and km are not on the same grid. Actually, a further 1715 !-- interpolation of km onto the u/v-grid is necessary. However, the 1716 !-- effect of this error is negligible. 1717 surf_usm_h%u_0(m) = u(k+1,j,i) + surf_usm_h%usws(m) * & 1718 ( zu(k+1) - zu(k-1) ) / & 1719 ( km(k,j,i) + 1.0E-20_wp ) 1720 surf_usm_h%v_0(m) = v(k+1,j,i) + surf_usm_h%vsws(m) * & 1721 ( zu(k+1) - zu(k-1) ) / & 1722 ( km(k,j,i) + 1.0E-20_wp ) 1723 1724 IF ( ABS( u(k+1,j,i) - surf_usm_h%u_0(m) ) > & 1725 ABS( u(k+1,j,i) - u(k-1,j,i) ) & 1726 ) surf_usm_h%u_0(m) = u(k-1,j,i) 1727 1728 IF ( ABS( v(k+1,j,i) - surf_usm_h%v_0(m) ) > & 1729 ABS( v(k+1,j,i) - v(k-1,j,i) ) & 1730 ) surf_usm_h%v_0(m) = v(k-1,j,i) 1731 1732 ENDDO 1338 1733 1339 1734 ENDIF
Note: See TracChangeset
for help on using the changeset viewer.
