Changeset 1342 for palm/trunk/SOURCE/production_e.f90
- Timestamp:
- Mar 26, 2014 5:04:47 PM (10 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/production_e.f90 (modified) (48 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/production_e.f90
r1321 r1342 20 20 ! Current revisions: 21 21 ! ----------------- 22 ! 22 ! REAL constants defined as wp-kind 23 23 ! 24 24 ! Former revisions: … … 173 173 DO k = nzb_diff_s_outer(j,i), nzt 174 174 175 dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx176 dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &177 u(k,j-1,i) - u(k,j-1,i+1) ) * ddy178 dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &179 u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)180 181 dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &182 v(k,j,i-1) - v(k,j+1,i-1) ) * ddx183 dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy184 dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &185 v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)186 187 dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &188 w(k,j,i-1) - w(k-1,j,i-1) ) * ddx189 dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &190 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy191 dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)192 193 def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) + &175 dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx 176 dudy = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - & 177 u(k,j-1,i) - u(k,j-1,i+1) ) * ddy 178 dudz = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - & 179 u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k) 180 181 dvdx = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - & 182 v(k,j,i-1) - v(k,j+1,i-1) ) * ddx 183 dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy 184 dvdz = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - & 185 v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k) 186 187 dwdx = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - & 188 w(k,j,i-1) - w(k-1,j,i-1) ) * ddx 189 dwdy = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - & 190 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy 191 dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) 192 193 def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 ) + & 194 194 dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + & 195 dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )196 197 IF ( def < 0.0 ) def = 0.0195 dvdz**2 + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz ) 196 197 IF ( def < 0.0_wp ) def = 0.0_wp 198 198 199 199 tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def … … 210 210 DO j = nys, nyn 211 211 212 IF ( ( wall_e_x(j,i) /= 0.0 ) .OR. ( wall_e_y(j,i) /= 0.0) ) &212 IF ( ( wall_e_x(j,i) /= 0.0_wp ) .OR. ( wall_e_y(j,i) /= 0.0_wp ) ) & 213 213 THEN 214 214 215 215 k = nzb_diff_s_inner(j,i) - 1 216 216 dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx 217 dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &218 u_0(j,i) - u_0(j,i+1) ) * dd2zu(k)217 dudz = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - & 218 u_0(j,i) - u_0(j,i+1) ) * dd2zu(k) 219 219 dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy 220 dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &221 v_0(j,i) - v_0(j+1,i) ) * dd2zu(k)220 dvdz = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - & 221 v_0(j,i) - v_0(j+1,i) ) * dd2zu(k) 222 222 dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) 223 223 224 IF ( wall_e_y(j,i) /= 0.0 ) THEN224 IF ( wall_e_y(j,i) /= 0.0_wp ) THEN 225 225 ! 226 226 !-- Inconsistency removed: as the thermal stratification is … … 236 236 CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, & 237 237 wsvs, 0.0_wp, 0.0_wp, 1.0_wp, 0.0_wp ) 238 km_neutral = kappa * ( usvs(k)**2 + wsvs(k)**2 )**0.25 * &239 0.5 * dy240 IF ( km_neutral > 0.0 ) THEN238 km_neutral = kappa * ( usvs(k)**2 + wsvs(k)**2 )**0.25_wp * & 239 0.5_wp * dy 240 IF ( km_neutral > 0.0_wp ) THEN 241 241 dudy = - wall_e_y(j,i) * usvs(k) / km_neutral 242 242 dwdy = - wall_e_y(j,i) * wsvs(k) / km_neutral 243 243 ELSE 244 dudy = 0.0 245 dwdy = 0.0 244 dudy = 0.0_wp 245 dwdy = 0.0_wp 246 246 ENDIF 247 247 ELSE 248 dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &249 u(k,j-1,i) - u(k,j-1,i+1) ) * ddy250 dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &251 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy248 dudy = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - & 249 u(k,j-1,i) - u(k,j-1,i+1) ) * ddy 250 dwdy = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - & 251 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy 252 252 ENDIF 253 253 254 IF ( wall_e_x(j,i) /= 0.0 ) THEN254 IF ( wall_e_x(j,i) /= 0.0_wp ) THEN 255 255 ! 256 256 !-- Inconsistency removed: as the thermal stratification is … … 266 266 CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, & 267 267 wsus, 0.0_wp, 0.0_wp, 0.0_wp, 1.0_wp ) 268 km_neutral = kappa * ( vsus(k)**2 + wsus(k)**2 )**0.25 * &269 0.5 * dx270 IF ( km_neutral > 0.0 ) THEN268 km_neutral = kappa * ( vsus(k)**2 + wsus(k)**2 )**0.25_wp * & 269 0.5_wp * dx 270 IF ( km_neutral > 0.0_wp ) THEN 271 271 dvdx = - wall_e_x(j,i) * vsus(k) / km_neutral 272 272 dwdx = - wall_e_x(j,i) * wsus(k) / km_neutral 273 273 ELSE 274 dvdx = 0.0 275 dwdx = 0.0 274 dvdx = 0.0_wp 275 dwdx = 0.0_wp 276 276 ENDIF 277 277 ELSE 278 dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &279 v(k,j,i-1) - v(k,j+1,i-1) ) * ddx280 dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &281 w(k,j,i-1) - w(k-1,j,i-1) ) * ddx278 dvdx = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - & 279 v(k,j,i-1) - v(k,j+1,i-1) ) * ddx 280 dwdx = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - & 281 w(k,j,i-1) - w(k-1,j,i-1) ) * ddx 282 282 ENDIF 283 283 284 def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) + &284 def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 ) + & 285 285 dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + & 286 dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )287 288 IF ( def < 0.0 ) def = 0.0286 dvdz**2 + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz ) 287 288 IF ( def < 0.0_wp ) def = 0.0_wp 289 289 290 290 tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def … … 301 301 302 302 dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx 303 dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &304 u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)305 dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy306 dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &307 v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)303 dudz = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - & 304 u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k) 305 dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy 306 dvdz = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - & 307 v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k) 308 308 dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) 309 309 310 IF ( wall_e_y(j,i) /= 0.0 ) THEN310 IF ( wall_e_y(j,i) /= 0.0_wp ) THEN 311 311 ! 312 312 !-- Inconsistency removed: as the thermal stratification … … 319 319 !-- validation has been available 320 320 km_neutral = kappa * ( usvs(k)**2 + & 321 wsvs(k)**2 )**0.25 * 0.5* dy322 IF ( km_neutral > 0.0 ) THEN321 wsvs(k)**2 )**0.25_wp * 0.5_wp * dy 322 IF ( km_neutral > 0.0_wp ) THEN 323 323 dudy = - wall_e_y(j,i) * usvs(k) / km_neutral 324 324 dwdy = - wall_e_y(j,i) * wsvs(k) / km_neutral 325 325 ELSE 326 dudy = 0.0 327 dwdy = 0.0 326 dudy = 0.0_wp 327 dwdy = 0.0_wp 328 328 ENDIF 329 329 ELSE 330 dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &331 u(k,j-1,i) - u(k,j-1,i+1) ) * ddy332 dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &333 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy330 dudy = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - & 331 u(k,j-1,i) - u(k,j-1,i+1) ) * ddy 332 dwdy = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - & 333 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy 334 334 ENDIF 335 335 336 IF ( wall_e_x(j,i) /= 0.0 ) THEN336 IF ( wall_e_x(j,i) /= 0.0_wp ) THEN 337 337 ! 338 338 !-- Inconsistency removed: as the thermal stratification … … 345 345 !-- validation has been available 346 346 km_neutral = kappa * ( vsus(k)**2 + & 347 wsus(k)**2 )**0.25 * 0.5* dx348 IF ( km_neutral > 0.0 ) THEN347 wsus(k)**2 )**0.25_wp * 0.5_wp * dx 348 IF ( km_neutral > 0.0_wp ) THEN 349 349 dvdx = - wall_e_x(j,i) * vsus(k) / km_neutral 350 350 dwdx = - wall_e_x(j,i) * wsus(k) / km_neutral 351 351 ELSE 352 dvdx = 0.0 353 dwdx = 0.0 352 dvdx = 0.0_wp 353 dwdx = 0.0_wp 354 354 ENDIF 355 355 ELSE 356 dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &357 v(k,j,i-1) - v(k,j+1,i-1) ) * ddx358 dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &359 w(k,j,i-1) - w(k-1,j,i-1) ) * ddx356 dvdx = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - & 357 v(k,j,i-1) - v(k,j+1,i-1) ) * ddx 358 dwdx = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - & 359 w(k,j,i-1) - w(k-1,j,i-1) ) * ddx 360 360 ENDIF 361 361 362 def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) + &362 def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 ) + & 363 363 dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + & 364 dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )365 366 IF ( def < 0.0 ) def = 0.0364 dvdz**2 + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz ) 365 366 IF ( def < 0.0_wp ) def = 0.0_wp 367 367 368 368 tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def … … 379 379 DO j = nys, nyn 380 380 381 IF ( ( wall_e_x(j,i) /= 0.0 ) .OR. ( wall_e_y(j,i) /= 0.0) ) &381 IF ( ( wall_e_x(j,i) /= 0.0_wp ) .OR. ( wall_e_y(j,i) /= 0.0_wp ) ) & 382 382 THEN 383 383 384 384 k = nzb_diff_s_outer(j,i)-1 385 385 386 dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx387 dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &388 u(k,j-1,i) - u(k,j-1,i+1) ) * ddy389 dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &390 u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)391 392 dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &393 v(k,j,i-1) - v(k,j+1,i-1) ) * ddx394 dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy395 dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &396 v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)397 398 dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &399 w(k,j,i-1) - w(k-1,j,i-1) ) * ddx400 dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &401 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy402 dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)403 404 def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) + &386 dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx 387 dudy = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - & 388 u(k,j-1,i) - u(k,j-1,i+1) ) * ddy 389 dudz = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - & 390 u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k) 391 392 dvdx = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - & 393 v(k,j,i-1) - v(k,j+1,i-1) ) * ddx 394 dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy 395 dvdz = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - & 396 v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k) 397 398 dwdx = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - & 399 w(k,j,i-1) - w(k-1,j,i-1) ) * ddx 400 dwdy = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - & 401 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy 402 dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) 403 404 def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 ) + & 405 405 dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + & 406 dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )407 408 IF ( def < 0.0 ) def = 0.0406 dvdz**2 + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz ) 407 408 IF ( def < 0.0_wp ) def = 0.0_wp 409 409 410 410 tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def … … 420 420 DO j = nys, nyn 421 421 422 IF ( ( wall_e_x(j,i) == 0.0 ) .AND. ( wall_e_y(j,i) == 0.0) ) &422 IF ( ( wall_e_x(j,i) == 0.0_wp ) .AND. ( wall_e_y(j,i) == 0.0_wp ) ) & 423 423 THEN 424 424 425 425 k = nzb_diff_s_inner(j,i)-1 426 426 427 dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx 428 dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - & 427 dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx 428 dudy = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - & 429 u(k,j-1,i) - u(k,j-1,i+1) ) * ddy 430 dudz = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - & 431 u_0(j,i) - u_0(j,i+1) ) * dd2zu(k) 432 433 dvdx = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - & 434 v(k,j,i-1) - v(k,j+1,i-1) ) * ddx 435 dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy 436 dvdz = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - & 437 v_0(j,i) - v_0(j+1,i) ) * dd2zu(k) 438 439 dwdx = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - & 440 w(k,j,i-1) - w(k-1,j,i-1) ) * ddx 441 dwdy = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - & 442 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy 443 dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) 444 445 def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 ) + & 446 dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + & 447 dvdz**2 + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz ) 448 449 IF ( def < 0.0_wp ) def = 0.0_wp 450 451 tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def 452 453 ENDIF 454 455 ENDDO 456 457 ELSEIF ( use_surface_fluxes ) THEN 458 459 DO j = nys, nyn 460 461 k = nzb_diff_s_outer(j,i)-1 462 463 dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx 464 dudy = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - & 429 465 u(k,j-1,i) - u(k,j-1,i+1) ) * ddy 430 dudz = 0.5* ( u(k+1,j,i) + u(k+1,j,i+1) - &431 u_0(j,i) - u_0(j,i+1)) * dd2zu(k)432 433 dvdx = 0.25* ( v(k,j,i+1) + v(k,j+1,i+1) - &466 dudz = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - & 467 u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k) 468 469 dvdx = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - & 434 470 v(k,j,i-1) - v(k,j+1,i-1) ) * ddx 435 dvdy =( v(k,j+1,i) - v(k,j,i) ) * ddy436 dvdz = 0.5* ( v(k+1,j,i) + v(k+1,j+1,i) - &437 v _0(j,i) - v_0(j+1,i)) * dd2zu(k)438 439 dwdx = 0.25* ( w(k,j,i+1) + w(k-1,j,i+1) - &471 dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy 472 dvdz = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - & 473 v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k) 474 475 dwdx = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - & 440 476 w(k,j,i-1) - w(k-1,j,i-1) ) * ddx 441 dwdy = 0.25* ( w(k,j+1,i) + w(k-1,j+1,i) - &477 dwdy = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - & 442 478 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy 443 dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) 444 445 def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) + & 446 dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + & 447 dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz ) 448 449 IF ( def < 0.0 ) def = 0.0 450 451 tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def 452 453 ENDIF 454 455 ENDDO 456 457 ELSEIF ( use_surface_fluxes ) THEN 458 459 DO j = nys, nyn 460 461 k = nzb_diff_s_outer(j,i)-1 462 463 dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx 464 dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - & 465 u(k,j-1,i) - u(k,j-1,i+1) ) * ddy 466 dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - & 467 u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k) 468 469 dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - & 470 v(k,j,i-1) - v(k,j+1,i-1) ) * ddx 471 dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy 472 dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - & 473 v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k) 474 475 dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - & 476 w(k,j,i-1) - w(k-1,j,i-1) ) * ddx 477 dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - & 478 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy 479 dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) 480 481 def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) + & 479 dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) 480 481 def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 ) + & 482 482 dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + & 483 dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )484 485 IF ( def < 0.0 ) def = 0.0483 dvdz**2 + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz ) 484 485 IF ( def < 0.0_wp ) def = 0.0_wp 486 486 487 487 tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def … … 586 586 587 587 IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN 588 k1 = 1.0 + 0.61* q(k,j,i)589 k2 = 0.61 * pt(k,j,i)588 k1 = 1.0_wp + 0.61_wp * q(k,j,i) 589 k2 = 0.61_wp * pt(k,j,i) 590 590 tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * & 591 591 g / vpt(k,j,i) * & … … 594 594 ) * dd2zu(k) 595 595 ELSE IF ( cloud_physics ) THEN 596 IF ( ql(k,j,i) == 0.0 ) THEN597 k1 = 1.0 + 0.61* q(k,j,i)598 k2 = 0.61 * pt(k,j,i)596 IF ( ql(k,j,i) == 0.0_wp ) THEN 597 k1 = 1.0_wp + 0.61_wp * q(k,j,i) 598 k2 = 0.61_wp * pt(k,j,i) 599 599 ELSE 600 600 theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i) 601 601 temp = theta * t_d_pt(k) 602 k1 = ( 1.0 - q(k,j,i) + 1.61 *&603 ( q(k,j,i) - ql(k,j,i) ) * &604 ( 1.0 + 0.622* l_d_r / temp ) ) / &605 ( 1.0 + 0.622* l_d_r * l_d_cp * &602 k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * & 603 ( q(k,j,i) - ql(k,j,i) ) * & 604 ( 1.0_wp + 0.622_wp * l_d_r / temp ) ) / & 605 ( 1.0_wp + 0.622_wp * l_d_r * l_d_cp * & 606 606 ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) ) 607 k2 = theta * ( l_d_cp / temp * k1 - 1.0 )607 k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp ) 608 608 ENDIF 609 609 tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * & … … 613 613 ) * dd2zu(k) 614 614 ELSE IF ( cloud_droplets ) THEN 615 k1 = 1.0 + 0.61* q(k,j,i) - ql(k,j,i)616 k2 = 0.61 * pt(k,j,i)615 k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i) 616 k2 = 0.61_wp * pt(k,j,i) 617 617 tend(k,j,i) = tend(k,j,i) - & 618 618 kh(k,j,i) * g / vpt(k,j,i) * & … … 634 634 635 635 IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN 636 k1 = 1.0 + 0.61* q(k,j,i)637 k2 = 0.61 * pt(k,j,i)636 k1 = 1.0_wp + 0.61_wp * q(k,j,i) 637 k2 = 0.61_wp * pt(k,j,i) 638 638 ELSE IF ( cloud_physics ) THEN 639 IF ( ql(k,j,i) == 0.0 ) THEN640 k1 = 1.0 + 0.61* q(k,j,i)641 k2 = 0.61 * pt(k,j,i)639 IF ( ql(k,j,i) == 0.0_wp ) THEN 640 k1 = 1.0_wp + 0.61_wp * q(k,j,i) 641 k2 = 0.61_wp * pt(k,j,i) 642 642 ELSE 643 643 theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i) 644 644 temp = theta * t_d_pt(k) 645 k1 = ( 1.0 - q(k,j,i) + 1.61 *&646 ( q(k,j,i) - ql(k,j,i) ) * &647 ( 1.0 + 0.622* l_d_r / temp ) ) / &648 ( 1.0 + 0.622* l_d_r * l_d_cp * &645 k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * & 646 ( q(k,j,i) - ql(k,j,i) ) * & 647 ( 1.0_wp + 0.622_wp * l_d_r / temp ) ) / & 648 ( 1.0_wp + 0.622_wp * l_d_r * l_d_cp * & 649 649 ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) ) 650 k2 = theta * ( l_d_cp / temp * k1 - 1.0 )650 k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp ) 651 651 ENDIF 652 652 ELSE IF ( cloud_droplets ) THEN 653 k1 = 1.0 + 0.61* q(k,j,i) - ql(k,j,i)654 k2 = 0.61 * pt(k,j,i)653 k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i) 654 k2 = 0.61_wp * pt(k,j,i) 655 655 ENDIF 656 656 … … 668 668 669 669 IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN 670 k1 = 1.0 + 0.61* q(k,j,i)671 k2 = 0.61 * pt(k,j,i)670 k1 = 1.0_wp + 0.61_wp * q(k,j,i) 671 k2 = 0.61_wp * pt(k,j,i) 672 672 ELSE IF ( cloud_physics ) THEN 673 IF ( ql(k,j,i) == 0.0 ) THEN674 k1 = 1.0 + 0.61* q(k,j,i)675 k2 = 0.61 * pt(k,j,i)673 IF ( ql(k,j,i) == 0.0_wp ) THEN 674 k1 = 1.0_wp + 0.61_wp * q(k,j,i) 675 k2 = 0.61_wp * pt(k,j,i) 676 676 ELSE 677 677 theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i) 678 678 temp = theta * t_d_pt(k) 679 k1 = ( 1.0 - q(k,j,i) + 1.61* &680 ( q(k,j,i) - ql(k,j,i) ) * &681 ( 1.0 + 0.622* l_d_r / temp ) ) / &682 ( 1.0 + 0.622* l_d_r * l_d_cp * &679 k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * & 680 ( q(k,j,i) - ql(k,j,i) ) * & 681 ( 1.0_wp + 0.622_wp * l_d_r / temp ) ) / & 682 ( 1.0_wp + 0.622_wp * l_d_r * l_d_cp * & 683 683 ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) ) 684 k2 = theta * ( l_d_cp / temp * k1 - 1.0 )684 k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp ) 685 685 ENDIF 686 686 ELSE IF ( cloud_droplets ) THEN 687 k1 = 1.0 + 0.61* q(k,j,i) - ql(k,j,i)688 k2 = 0.61 * pt(k,j,i)687 k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i) 688 k2 = 0.61_wp * pt(k,j,i) 689 689 ENDIF 690 690 … … 783 783 IF ( k >= nzb_diff_s_outer(j,i) ) THEN 784 784 785 dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx786 dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &787 u(k,j-1,i) - u(k,j-1,i+1) ) * ddy788 dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &789 u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)790 791 dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &792 v(k,j,i-1) - v(k,j+1,i-1) ) * ddx793 dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy794 dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &795 v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)796 797 dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &798 w(k,j,i-1) - w(k-1,j,i-1) ) * ddx799 dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &800 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy801 dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)802 803 def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) + &785 dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx 786 dudy = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - & 787 u(k,j-1,i) - u(k,j-1,i+1) ) * ddy 788 dudz = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - & 789 u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k) 790 791 dvdx = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - & 792 v(k,j,i-1) - v(k,j+1,i-1) ) * ddx 793 dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy 794 dvdz = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - & 795 v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k) 796 797 dwdx = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - & 798 w(k,j,i-1) - w(k-1,j,i-1) ) * ddx 799 dwdy = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - & 800 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy 801 dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) 802 803 def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 ) + & 804 804 dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + & 805 dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )806 807 IF ( def < 0.0 ) def = 0.0805 dvdz**2 + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz ) 806 807 IF ( def < 0.0_wp ) def = 0.0_wp 808 808 809 809 tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def … … 825 825 DO k = 1, nzt 826 826 827 IF ( ( wall_e_x(j,i) /= 0.0 ).OR.( wall_e_y(j,i) /= 0.0) ) &827 IF ( ( wall_e_x(j,i) /= 0.0_wp ).OR.( wall_e_y(j,i) /= 0.0_wp ) ) & 828 828 THEN 829 829 830 830 IF ( k == nzb_diff_s_inner(j,i) - 1 ) THEN 831 831 dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx 832 dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &833 u_0(j,i) - u_0(j,i+1) ) * dd2zu(k)832 dudz = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - & 833 u_0(j,i) - u_0(j,i+1) ) * dd2zu(k) 834 834 dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy 835 dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &836 v_0(j,i) - v_0(j+1,i) ) * dd2zu(k)835 dvdz = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - & 836 v_0(j,i) - v_0(j+1,i) ) * dd2zu(k) 837 837 dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) 838 838 839 IF ( wall_e_y(j,i) /= 0.0 ) THEN839 IF ( wall_e_y(j,i) /= 0.0_wp ) THEN 840 840 ! 841 841 !-- Inconsistency removed: as the thermal stratification is … … 852 852 ! wsvs, 0.0_wp, 0.0_wp, 1.0_wp, 0.0_wp ) 853 853 km_neutral = kappa * & 854 ( usvs(k,j,i)**2 + wsvs(k,j,i)**2 )**0.25 * &855 0.5 * dy856 IF ( km_neutral > 0.0 ) THEN854 ( usvs(k,j,i)**2 + wsvs(k,j,i)**2 )**0.25_wp * & 855 0.5_wp * dy 856 IF ( km_neutral > 0.0_wp ) THEN 857 857 dudy = - wall_e_y(j,i) * usvs(k,j,i) / km_neutral 858 858 dwdy = - wall_e_y(j,i) * wsvs(k,j,i) / km_neutral 859 859 ELSE 860 dudy = 0.0 861 dwdy = 0.0 860 dudy = 0.0_wp 861 dwdy = 0.0_wp 862 862 ENDIF 863 863 ELSE 864 dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &865 u(k,j-1,i) - u(k,j-1,i+1) ) * ddy866 dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &867 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy864 dudy = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - & 865 u(k,j-1,i) - u(k,j-1,i+1) ) * ddy 866 dwdy = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - & 867 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy 868 868 ENDIF 869 869 870 IF ( wall_e_x(j,i) /= 0.0 ) THEN870 IF ( wall_e_x(j,i) /= 0.0_wp ) THEN 871 871 ! 872 872 !-- Inconsistency removed: as the thermal stratification is … … 883 883 ! wsus, 0.0_wp, 0.0_wp, 0.0_wp, 1.0_wp ) 884 884 km_neutral = kappa * & 885 ( vsus(k,j,i)**2 + wsus(k,j,i)**2 )**0.25 * &886 0.5 * dx887 IF ( km_neutral > 0.0 ) THEN885 ( vsus(k,j,i)**2 + wsus(k,j,i)**2 )**0.25_wp * & 886 0.5_wp * dx 887 IF ( km_neutral > 0.0_wp ) THEN 888 888 dvdx = - wall_e_x(j,i) * vsus(k,j,i) / km_neutral 889 889 dwdx = - wall_e_x(j,i) * wsus(k,j,i) / km_neutral 890 890 ELSE 891 dvdx = 0.0 892 dwdx = 0.0 891 dvdx = 0.0_wp 892 dwdx = 0.0_wp 893 893 ENDIF 894 894 ELSE 895 dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &896 v(k,j,i-1) - v(k,j+1,i-1) ) * ddx897 dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &898 w(k,j,i-1) - w(k-1,j,i-1) ) * ddx895 dvdx = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - & 896 v(k,j,i-1) - v(k,j+1,i-1) ) * ddx 897 dwdx = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - & 898 w(k,j,i-1) - w(k-1,j,i-1) ) * ddx 899 899 ENDIF 900 900 901 def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) + &901 def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 ) + & 902 902 dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + & 903 dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )904 905 IF ( def < 0.0 ) def = 0.0903 dvdz**2 + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz ) 904 905 IF ( def < 0.0_wp ) def = 0.0_wp 906 906 907 907 tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def … … 919 919 920 920 dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx 921 dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &922 u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)923 dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy924 dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &925 v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)921 dudz = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - & 922 u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k) 923 dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy 924 dvdz = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - & 925 v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k) 926 926 dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) 927 927 928 IF ( wall_e_y(j,i) /= 0.0 ) THEN928 IF ( wall_e_y(j,i) /= 0.0_wp ) THEN 929 929 ! 930 930 !-- Inconsistency removed: as the thermal stratification … … 937 937 !-- validation has been available 938 938 km_neutral = kappa * ( usvs(k,j,i)**2 + & 939 wsvs(k,j,i)**2 )**0.25 * 0.5* dy940 IF ( km_neutral > 0.0 ) THEN939 wsvs(k,j,i)**2 )**0.25_wp * 0.5_wp * dy 940 IF ( km_neutral > 0.0_wp ) THEN 941 941 dudy = - wall_e_y(j,i) * usvs(k,j,i) / km_neutral 942 942 dwdy = - wall_e_y(j,i) * wsvs(k,j,i) / km_neutral 943 943 ELSE 944 dudy = 0.0 945 dwdy = 0.0 944 dudy = 0.0_wp 945 dwdy = 0.0_wp 946 946 ENDIF 947 947 ELSE 948 dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &949 u(k,j-1,i) - u(k,j-1,i+1) ) * ddy950 dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &951 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy948 dudy = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - & 949 u(k,j-1,i) - u(k,j-1,i+1) ) * ddy 950 dwdy = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - & 951 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy 952 952 ENDIF 953 953 954 IF ( wall_e_x(j,i) /= 0.0 ) THEN954 IF ( wall_e_x(j,i) /= 0.0_wp ) THEN 955 955 ! 956 956 !-- Inconsistency removed: as the thermal stratification … … 963 963 !-- validation has been available 964 964 km_neutral = kappa * ( vsus(k,j,i)**2 + & 965 wsus(k,j,i)**2 )**0.25 * 0.5* dx966 IF ( km_neutral > 0.0 ) THEN965 wsus(k,j,i)**2 )**0.25_wp * 0.5_wp * dx 966 IF ( km_neutral > 0.0_wp ) THEN 967 967 dvdx = - wall_e_x(j,i) * vsus(k,j,i) / km_neutral 968 968 dwdx = - wall_e_x(j,i) * wsus(k,j,i) / km_neutral 969 969 ELSE 970 dvdx = 0.0 971 dwdx = 0.0 970 dvdx = 0.0_wp 971 dwdx = 0.0_wp 972 972 ENDIF 973 973 ELSE 974 dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &975 v(k,j,i-1) - v(k,j+1,i-1) ) * ddx976 dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &977 w(k,j,i-1) - w(k-1,j,i-1) ) * ddx974 dvdx = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - & 975 v(k,j,i-1) - v(k,j+1,i-1) ) * ddx 976 dwdx = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - & 977 w(k,j,i-1) - w(k-1,j,i-1) ) * ddx 978 978 ENDIF 979 979 980 def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) + &980 def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 ) + & 981 981 dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + & 982 dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )983 984 IF ( def < 0.0 ) def = 0.0982 dvdz**2 + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz ) 983 984 IF ( def < 0.0_wp ) def = 0.0_wp 985 985 986 986 tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def … … 993 993 IF ( k == nzb_diff_s_outer(j,i)-1 ) THEN 994 994 995 dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx996 dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &997 u(k,j-1,i) - u(k,j-1,i+1) ) * ddy998 dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &999 u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)1000 1001 dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &1002 v(k,j,i-1) - v(k,j+1,i-1) ) * ddx1003 dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy1004 dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &1005 v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)1006 1007 dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &1008 w(k,j,i-1) - w(k-1,j,i-1) ) * ddx1009 dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &1010 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy1011 dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)1012 1013 def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) + &995 dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx 996 dudy = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - & 997 u(k,j-1,i) - u(k,j-1,i+1) ) * ddy 998 dudz = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - & 999 u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k) 1000 1001 dvdx = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - & 1002 v(k,j,i-1) - v(k,j+1,i-1) ) * ddx 1003 dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy 1004 dvdz = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - & 1005 v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k) 1006 1007 dwdx = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - & 1008 w(k,j,i-1) - w(k-1,j,i-1) ) * ddx 1009 dwdy = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - & 1010 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy 1011 dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) 1012 1013 def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 ) + & 1014 1014 dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + & 1015 dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )1016 1017 IF ( def < 0.0 ) def = 0.01015 dvdz**2 + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz ) 1016 1017 IF ( def < 0.0_wp ) def = 0.0_wp 1018 1018 1019 1019 tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def … … 1035 1035 DO k = 1, nzt 1036 1036 1037 IF ( ( wall_e_x(j,i) == 0.0 ) .AND. ( wall_e_y(j,i) == 0.0) ) &1037 IF ( ( wall_e_x(j,i) == 0.0_wp ) .AND. ( wall_e_y(j,i) == 0.0_wp ) ) & 1038 1038 THEN 1039 1039 1040 1040 IF ( k == nzb_diff_s_inner(j,i)-1 ) THEN 1041 1041 1042 dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx1043 dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &1044 u(k,j-1,i) - u(k,j-1,i+1) ) * ddy1045 dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &1046 u_0(j,i) - u_0(j,i+1) ) * dd2zu(k)1047 1048 dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &1049 v(k,j,i-1) - v(k,j+1,i-1) ) * ddx1050 dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy1051 dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &1052 v_0(j,i) - v_0(j+1,i) ) * dd2zu(k)1053 1054 dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &1055 w(k,j,i-1) - w(k-1,j,i-1) ) * ddx1056 dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &1057 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy1058 dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)1059 1060 def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) + &1042 dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx 1043 dudy = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - & 1044 u(k,j-1,i) - u(k,j-1,i+1) ) * ddy 1045 dudz = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - & 1046 u_0(j,i) - u_0(j,i+1) ) * dd2zu(k) 1047 1048 dvdx = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - & 1049 v(k,j,i-1) - v(k,j+1,i-1) ) * ddx 1050 dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy 1051 dvdz = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - & 1052 v_0(j,i) - v_0(j+1,i) ) * dd2zu(k) 1053 1054 dwdx = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - & 1055 w(k,j,i-1) - w(k-1,j,i-1) ) * ddx 1056 dwdy = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - & 1057 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy 1058 dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) 1059 1060 def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 ) + & 1061 1061 dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + & 1062 dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )1063 1064 IF ( def < 0.0 ) def = 0.01062 dvdz**2 + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz ) 1063 1064 IF ( def < 0.0_wp ) def = 0.0_wp 1065 1065 1066 1066 tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def … … 1082 1082 IF ( k == nzb_diff_s_outer(j,i)-1 ) THEN 1083 1083 1084 dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx1085 dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &1086 u(k,j-1,i) - u(k,j-1,i+1) ) * ddy1087 dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &1088 u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)1089 1090 dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &1091 v(k,j,i-1) - v(k,j+1,i-1) ) * ddx1092 dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy1093 dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &1094 v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)1095 1096 dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &1097 w(k,j,i-1) - w(k-1,j,i-1) ) * ddx1098 dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &1099 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy1100 dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)1101 1102 def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) + &1084 dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx 1085 dudy = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - & 1086 u(k,j-1,i) - u(k,j-1,i+1) ) * ddy 1087 dudz = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - & 1088 u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k) 1089 1090 dvdx = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - & 1091 v(k,j,i-1) - v(k,j+1,i-1) ) * ddx 1092 dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy 1093 dvdz = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - & 1094 v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k) 1095 1096 dwdx = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - & 1097 w(k,j,i-1) - w(k-1,j,i-1) ) * ddx 1098 dwdy = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - & 1099 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy 1100 dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) 1101 1102 def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 ) + & 1103 1103 dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + & 1104 dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )1105 1106 IF ( def < 0.0 ) def = 0.01104 dvdz**2 + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz ) 1105 1106 IF ( def < 0.0_wp ) def = 0.0_wp 1107 1107 1108 1108 tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def … … 1232 1232 ! 1233 1233 ! IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN 1234 ! k1 = 1.0 + 0.61* q(k,j,i)1235 ! k2 = 0.61 * pt(k,j,i)1234 ! k1 = 1.0_wp + 0.61_wp * q(k,j,i) 1235 ! k2 = 0.61_wp * pt(k,j,i) 1236 1236 ! tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * & 1237 1237 ! g / vpt(k,j,i) * & … … 1240 1240 ! ) * dd2zu(k) 1241 1241 ! ELSE IF ( cloud_physics ) THEN 1242 ! IF ( ql(k,j,i) == 0.0 ) THEN1243 ! k1 = 1.0 + 0.61* q(k,j,i)1244 ! k2 = 0.61 * pt(k,j,i)1242 ! IF ( ql(k,j,i) == 0.0_wp ) THEN 1243 ! k1 = 1.0_wp + 0.61_wp * q(k,j,i) 1244 ! k2 = 0.61_wp * pt(k,j,i) 1245 1245 ! ELSE 1246 1246 ! theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i) 1247 1247 ! temp = theta * t_d_pt(k) 1248 ! k1 = ( 1.0 - q(k,j,i) + 1.61* &1249 ! ( q(k,j,i) - ql(k,j,i) ) * &1250 ! ( 1.0 + 0.622* l_d_r / temp ) ) / &1251 ! ( 1.0 + 0.622* l_d_r * l_d_cp * &1248 ! k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * & 1249 ! ( q(k,j,i) - ql(k,j,i) ) * & 1250 ! ( 1.0_wp + 0.622_wp * l_d_r / temp ) ) / & 1251 ! ( 1.0_wp + 0.622_wp * l_d_r * l_d_cp * & 1252 1252 ! ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) ) 1253 ! k2 = theta * ( l_d_cp / temp * k1 - 1.0 )1253 ! k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp ) 1254 1254 ! ENDIF 1255 1255 ! tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * & … … 1259 1259 ! ) * dd2zu(k) 1260 1260 ! ELSE IF ( cloud_droplets ) THEN 1261 ! k1 = 1.0 + 0.61* q(k,j,i) - ql(k,j,i)1262 ! k2 = 0.61 * pt(k,j,i)1261 ! k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i) 1262 ! k2 = 0.61_wp * pt(k,j,i) 1263 1263 ! tend(k,j,i) = tend(k,j,i) - & 1264 1264 ! kh(k,j,i) * g / vpt(k,j,i) * & … … 1289 1289 ! 1290 1290 ! IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN 1291 ! k1 = 1.0 + 0.61* q(k,j,i)1292 ! k2 = 0.61 * pt(k,j,i)1291 ! k1 = 1.0_wp + 0.61_wp * q(k,j,i) 1292 ! k2 = 0.61_wp * pt(k,j,i) 1293 1293 ! ELSE IF ( cloud_physics ) THEN 1294 ! IF ( ql(k,j,i) == 0.0 ) THEN1295 ! k1 = 1.0 + 0.61* q(k,j,i)1296 ! k2 = 0.61 * pt(k,j,i)1294 ! IF ( ql(k,j,i) == 0.0_wp ) THEN 1295 ! k1 = 1.0_wp + 0.61_wp * q(k,j,i) 1296 ! k2 = 0.61_wp * pt(k,j,i) 1297 1297 ! ELSE 1298 1298 ! theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i) 1299 1299 ! temp = theta * t_d_pt(k) 1300 ! k1 = ( 1.0 - q(k,j,i) + 1.61 *&1301 ! ( q(k,j,i) - ql(k,j,i) ) *&1302 ! ( 1.0 + 0.622 * l_d_r / temp ) ) /&1303 ! ( 1.0 + 0.622 * l_d_r * l_d_cp *&1300 ! k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * & 1301 ! ( q(k,j,i) - ql(k,j,i) ) * & 1302 ! ( 1.0_wp + 0.622_wp * l_d_r / temp ) ) /& 1303 ! ( 1.0_wp + 0.622_wp * l_d_r * l_d_cp * & 1304 1304 ! ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) ) 1305 ! k2 = theta * ( l_d_cp / temp * k1 - 1.0 )1305 ! k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp ) 1306 1306 ! ENDIF 1307 1307 ! ELSE IF ( cloud_droplets ) THEN 1308 ! k1 = 1.0 + 0.61* q(k,j,i) - ql(k,j,i)1309 ! k2 = 0.61 * pt(k,j,i)1308 ! k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i) 1309 ! k2 = 0.61_wp * pt(k,j,i) 1310 1310 ! ENDIF 1311 1311 ! … … 1330 1330 ! 1331 1331 ! IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN 1332 ! k1 = 1.0 + 0.61* q(k,j,i)1333 ! k2 = 0.61 * pt(k,j,i)1332 ! k1 = 1.0_wp + 0.61_wp * q(k,j,i) 1333 ! k2 = 0.61_wp * pt(k,j,i) 1334 1334 ! ELSE IF ( cloud_physics ) THEN 1335 ! IF ( ql(k,j,i) == 0.0 ) THEN1336 ! k1 = 1.0 + 0.61* q(k,j,i)1337 ! k2 = 0.61 * pt(k,j,i)1335 ! IF ( ql(k,j,i) == 0.0_wp ) THEN 1336 ! k1 = 1.0_wp + 0.61_wp * q(k,j,i) 1337 ! k2 = 0.61_wp * pt(k,j,i) 1338 1338 ! ELSE 1339 1339 ! theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i) 1340 1340 ! temp = theta * t_d_pt(k) 1341 ! k1 = ( 1.0 - q(k,j,i) + 1.61 *&1342 ! ( q(k,j,i) - ql(k,j,i) ) *&1343 ! ( 1.0 + 0.622 * l_d_r / temp ) ) /&1344 ! ( 1.0 + 0.622 * l_d_r * l_d_cp *&1341 ! k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * & 1342 ! ( q(k,j,i) - ql(k,j,i) ) * & 1343 ! ( 1.0_wp + 0.622_wp * l_d_r / temp ) ) /& 1344 ! ( 1.0_wp + 0.622_wp * l_d_r * l_d_cp * & 1345 1345 ! ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) ) 1346 ! k2 = theta * ( l_d_cp / temp * k1 - 1.0 )1346 ! k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp ) 1347 1347 ! ENDIF 1348 1348 ! ELSE IF ( cloud_droplets ) THEN 1349 ! k1 = 1.0 + 0.61* q(k,j,i) - ql(k,j,i)1350 ! k2 = 0.61 * pt(k,j,i)1349 ! k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i) 1350 ! k2 = 0.61_wp * pt(k,j,i) 1351 1351 ! ENDIF 1352 1352 ! … … 1426 1426 DO k = nzb_diff_s_outer(j,i), nzt 1427 1427 1428 dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx1429 dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &1430 u(k,j-1,i) - u(k,j-1,i+1) ) * ddy1431 dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &1432 u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)1433 1434 dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &1435 v(k,j,i-1) - v(k,j+1,i-1) ) * ddx1436 dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy1437 dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &1438 v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)1439 1440 dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &1441 w(k,j,i-1) - w(k-1,j,i-1) ) * ddx1442 dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &1443 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy1444 dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)1445 1446 def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) &1447 + dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + dvdz**2 &1448 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )1449 1450 IF ( def < 0.0 ) def = 0.01428 dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx 1429 dudy = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - & 1430 u(k,j-1,i) - u(k,j-1,i+1) ) * ddy 1431 dudz = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - & 1432 u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k) 1433 1434 dvdx = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - & 1435 v(k,j,i-1) - v(k,j+1,i-1) ) * ddx 1436 dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy 1437 dvdz = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - & 1438 v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k) 1439 1440 dwdx = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - & 1441 w(k,j,i-1) - w(k-1,j,i-1) ) * ddx 1442 dwdy = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - & 1443 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy 1444 dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) 1445 1446 def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 ) & 1447 + dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + dvdz**2 & 1448 + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz ) 1449 1450 IF ( def < 0.0_wp ) def = 0.0_wp 1451 1451 1452 1452 tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def … … 1456 1456 IF ( prandtl_layer ) THEN 1457 1457 1458 IF ( ( wall_e_x(j,i) /= 0.0 ) .OR. ( wall_e_y(j,i) /= 0.0) ) THEN1458 IF ( ( wall_e_x(j,i) /= 0.0_wp ) .OR. ( wall_e_y(j,i) /= 0.0_wp ) ) THEN 1459 1459 1460 1460 ! … … 1465 1465 1466 1466 dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx 1467 dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &1468 u_0(j,i) - u_0(j,i+1) ) * dd2zu(k)1469 dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy1470 dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &1471 v_0(j,i) - v_0(j+1,i) ) * dd2zu(k)1467 dudz = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - & 1468 u_0(j,i) - u_0(j,i+1) ) * dd2zu(k) 1469 dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy 1470 dvdz = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - & 1471 v_0(j,i) - v_0(j+1,i) ) * dd2zu(k) 1472 1472 dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) 1473 1473 1474 IF ( wall_e_y(j,i) /= 0.0 ) THEN1474 IF ( wall_e_y(j,i) /= 0.0_wp ) THEN 1475 1475 ! 1476 1476 !-- Inconsistency removed: as the thermal stratification … … 1486 1486 CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, & 1487 1487 wsvs, 0.0_wp, 0.0_wp, 1.0_wp, 0.0_wp ) 1488 km_neutral = kappa * ( usvs(k)**2 + wsvs(k)**2 )**0.25 * &1489 0.5 * dy1490 IF ( km_neutral > 0.0 ) THEN1488 km_neutral = kappa * ( usvs(k)**2 + wsvs(k)**2 )**0.25_wp * & 1489 0.5_wp * dy 1490 IF ( km_neutral > 0.0_wp ) THEN 1491 1491 dudy = - wall_e_y(j,i) * usvs(k) / km_neutral 1492 1492 dwdy = - wall_e_y(j,i) * wsvs(k) / km_neutral 1493 1493 ELSE 1494 dudy = 0.0 1495 dwdy = 0.0 1494 dudy = 0.0_wp 1495 dwdy = 0.0_wp 1496 1496 ENDIF 1497 1497 ELSE 1498 dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &1499 u(k,j-1,i) - u(k,j-1,i+1) ) * ddy1500 dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &1501 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy1498 dudy = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - & 1499 u(k,j-1,i) - u(k,j-1,i+1) ) * ddy 1500 dwdy = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - & 1501 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy 1502 1502 ENDIF 1503 1503 1504 IF ( wall_e_x(j,i) /= 0.0 ) THEN1504 IF ( wall_e_x(j,i) /= 0.0_wp ) THEN 1505 1505 ! 1506 1506 !-- Inconsistency removed: as the thermal stratification … … 1516 1516 CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, & 1517 1517 wsus, 0.0_wp, 0.0_wp, 0.0_wp, 1.0_wp ) 1518 km_neutral = kappa * ( vsus(k)**2 + wsus(k)**2 )**0.25 * &1519 0.5 * dx1520 IF ( km_neutral > 0.0 ) THEN1518 km_neutral = kappa * ( vsus(k)**2 + wsus(k)**2 )**0.25_wp * & 1519 0.5_wp * dx 1520 IF ( km_neutral > 0.0_wp ) THEN 1521 1521 dvdx = - wall_e_x(j,i) * vsus(k) / km_neutral 1522 1522 dwdx = - wall_e_x(j,i) * wsus(k) / km_neutral 1523 1523 ELSE 1524 dvdx = 0.0 1525 dwdx = 0.0 1524 dvdx = 0.0_wp 1525 dwdx = 0.0_wp 1526 1526 ENDIF 1527 1527 ELSE 1528 dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &1529 v(k,j,i-1) - v(k,j+1,i-1) ) * ddx1530 dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &1531 w(k,j,i-1) - w(k-1,j,i-1) ) * ddx1528 dvdx = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - & 1529 v(k,j,i-1) - v(k,j+1,i-1) ) * ddx 1530 dwdx = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - & 1531 w(k,j,i-1) - w(k-1,j,i-1) ) * ddx 1532 1532 ENDIF 1533 1533 1534 def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) + &1534 def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 ) + & 1535 1535 dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + & 1536 dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )1537 1538 IF ( def < 0.0 ) def = 0.01536 dvdz**2 + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz ) 1537 1538 IF ( def < 0.0_wp ) def = 0.0_wp 1539 1539 1540 1540 tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def … … 1549 1549 1550 1550 dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx 1551 dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &1552 u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)1553 dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy1554 dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &1555 v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)1551 dudz = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - & 1552 u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k) 1553 dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy 1554 dvdz = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - & 1555 v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k) 1556 1556 dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) 1557 1557 1558 IF ( wall_e_y(j,i) /= 0.0 ) THEN1558 IF ( wall_e_y(j,i) /= 0.0_wp ) THEN 1559 1559 ! 1560 1560 !-- Inconsistency removed: as the thermal stratification … … 1567 1567 !-- validation has been available 1568 1568 km_neutral = kappa * ( usvs(k)**2 + & 1569 wsvs(k)**2 )**0.25 * 0.5* dy1570 IF ( km_neutral > 0.0 ) THEN1569 wsvs(k)**2 )**0.25_wp * 0.5_wp * dy 1570 IF ( km_neutral > 0.0_wp ) THEN 1571 1571 dudy = - wall_e_y(j,i) * usvs(k) / km_neutral 1572 1572 dwdy = - wall_e_y(j,i) * wsvs(k) / km_neutral 1573 1573 ELSE 1574 dudy = 0.0 1575 dwdy = 0.0 1574 dudy = 0.0_wp 1575 dwdy = 0.0_wp 1576 1576 ENDIF 1577 1577 ELSE 1578 dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &1579 u(k,j-1,i) - u(k,j-1,i+1) ) * ddy1580 dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &1581 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy1578 dudy = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - & 1579 u(k,j-1,i) - u(k,j-1,i+1) ) * ddy 1580 dwdy = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - & 1581 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy 1582 1582 ENDIF 1583 1583 1584 IF ( wall_e_x(j,i) /= 0.0 ) THEN1584 IF ( wall_e_x(j,i) /= 0.0_wp ) THEN 1585 1585 ! 1586 1586 !-- Inconsistency removed: as the thermal stratification … … 1593 1593 !-- validation has been available 1594 1594 km_neutral = kappa * ( vsus(k)**2 + & 1595 wsus(k)**2 )**0.25 * 0.5* dx1596 IF ( km_neutral > 0.0 ) THEN1595 wsus(k)**2 )**0.25_wp * 0.5_wp * dx 1596 IF ( km_neutral > 0.0_wp ) THEN 1597 1597 dvdx = - wall_e_x(j,i) * vsus(k) / km_neutral 1598 1598 dwdx = - wall_e_x(j,i) * wsus(k) / km_neutral 1599 1599 ELSE 1600 dvdx = 0.0 1601 dwdx = 0.0 1600 dvdx = 0.0_wp 1601 dwdx = 0.0_wp 1602 1602 ENDIF 1603 1603 ELSE 1604 dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &1605 v(k,j,i-1) - v(k,j+1,i-1) ) * ddx1606 dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &1607 w(k,j,i-1) - w(k-1,j,i-1) ) * ddx1604 dvdx = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - & 1605 v(k,j,i-1) - v(k,j+1,i-1) ) * ddx 1606 dwdx = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - & 1607 w(k,j,i-1) - w(k-1,j,i-1) ) * ddx 1608 1608 ENDIF 1609 1609 1610 def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) + &1610 def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 ) + & 1611 1611 dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + & 1612 dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )1613 1614 IF ( def < 0.0 ) def = 0.01612 dvdz**2 + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz ) 1613 1614 IF ( def < 0.0_wp ) def = 0.0_wp 1615 1615 1616 1616 tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def … … 1623 1623 k = nzb_diff_s_outer(j,i)-1 1624 1624 1625 dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx1626 dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &1627 u(k,j-1,i) - u(k,j-1,i+1) ) * ddy1628 dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &1629 u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)1630 1631 dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &1632 v(k,j,i-1) - v(k,j+1,i-1) ) * ddx1633 dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy1634 dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &1635 v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)1636 1637 dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &1638 w(k,j,i-1) - w(k-1,j,i-1) ) * ddx1639 dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &1640 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy1641 dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)1625 dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx 1626 dudy = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - & 1627 u(k,j-1,i) - u(k,j-1,i+1) ) * ddy 1628 dudz = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - & 1629 u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k) 1630 1631 dvdx = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - & 1632 v(k,j,i-1) - v(k,j+1,i-1) ) * ddx 1633 dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy 1634 dvdz = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - & 1635 v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k) 1636 1637 dwdx = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - & 1638 w(k,j,i-1) - w(k-1,j,i-1) ) * ddx 1639 dwdy = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - & 1640 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy 1641 dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) 1642 1642 1643 1643 def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) + & … … 1645 1645 dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz ) 1646 1646 1647 IF ( def < 0.0 ) def = 0.01647 IF ( def < 0.0_wp ) def = 0.0_wp 1648 1648 1649 1649 tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def … … 1657 1657 k = nzb_diff_s_inner(j,i)-1 1658 1658 1659 dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx 1660 dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - & 1659 dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx 1660 dudy = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - & 1661 u(k,j-1,i) - u(k,j-1,i+1) ) * ddy 1662 dudz = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - & 1663 u_0(j,i) - u_0(j,i+1) ) * dd2zu(k) 1664 1665 dvdx = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - & 1666 v(k,j,i-1) - v(k,j+1,i-1) ) * ddx 1667 dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy 1668 dvdz = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - & 1669 v_0(j,i) - v_0(j+1,i) ) * dd2zu(k) 1670 1671 dwdx = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - & 1672 w(k,j,i-1) - w(k-1,j,i-1) ) * ddx 1673 dwdy = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - & 1674 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy 1675 dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) 1676 1677 def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 ) & 1678 + dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + dvdz**2 & 1679 + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz ) 1680 1681 IF ( def < 0.0_wp ) def = 0.0_wp 1682 1683 tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def 1684 1685 ENDIF 1686 1687 ELSEIF ( use_surface_fluxes ) THEN 1688 1689 k = nzb_diff_s_outer(j,i)-1 1690 1691 dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx 1692 dudy = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - & 1661 1693 u(k,j-1,i) - u(k,j-1,i+1) ) * ddy 1662 dudz = 0.5* ( u(k+1,j,i) + u(k+1,j,i+1) - &1663 u _0(j,i) - u_0(j,i+1)) * dd2zu(k)1664 1665 dvdx = 0.25* ( v(k,j,i+1) + v(k,j+1,i+1) - &1694 dudz = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - & 1695 u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k) 1696 1697 dvdx = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - & 1666 1698 v(k,j,i-1) - v(k,j+1,i-1) ) * ddx 1667 dvdy =( v(k,j+1,i) - v(k,j,i) ) * ddy1668 dvdz = 0.5* ( v(k+1,j,i) + v(k+1,j+1,i) - &1669 v _0(j,i) - v_0(j+1,i)) * dd2zu(k)1670 1671 dwdx = 0.25* ( w(k,j,i+1) + w(k-1,j,i+1) - &1699 dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy 1700 dvdz = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - & 1701 v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k) 1702 1703 dwdx = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - & 1672 1704 w(k,j,i-1) - w(k-1,j,i-1) ) * ddx 1673 dwdy = 0.25* ( w(k,j+1,i) + w(k-1,j+1,i) - &1705 dwdy = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - & 1674 1706 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy 1675 dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) 1676 1677 def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) & 1678 + dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + dvdz**2 & 1679 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz ) 1680 1681 IF ( def < 0.0 ) def = 0.0 1682 1683 tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def 1684 1685 ENDIF 1686 1687 ELSEIF ( use_surface_fluxes ) THEN 1688 1689 k = nzb_diff_s_outer(j,i)-1 1690 1691 dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx 1692 dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - & 1693 u(k,j-1,i) - u(k,j-1,i+1) ) * ddy 1694 dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - & 1695 u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k) 1696 1697 dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - & 1698 v(k,j,i-1) - v(k,j+1,i-1) ) * ddx 1699 dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy 1700 dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - & 1701 v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k) 1702 1703 dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - & 1704 w(k,j,i-1) - w(k-1,j,i-1) ) * ddx 1705 dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - & 1706 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy 1707 dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) 1708 1709 def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) + & 1707 dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) 1708 1709 def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 ) + & 1710 1710 dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + & 1711 dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )1712 1713 IF ( def < 0.0 ) def = 0.01711 dvdz**2 + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz ) 1712 1713 IF ( def < 0.0_wp ) def = 0.0_wp 1714 1714 1715 1715 tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def … … 1794 1794 1795 1795 IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN 1796 k1 = 1.0 + 0.61* q(k,j,i)1797 k2 = 0.61 * pt(k,j,i)1796 k1 = 1.0_wp + 0.61_wp * q(k,j,i) 1797 k2 = 0.61_wp * pt(k,j,i) 1798 1798 tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / vpt(k,j,i) * & 1799 1799 ( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) + & … … 1801 1801 ) * dd2zu(k) 1802 1802 ELSE IF ( cloud_physics ) THEN 1803 IF ( ql(k,j,i) == 0.0 ) THEN1804 k1 = 1.0 + 0.61* q(k,j,i)1805 k2 = 0.61 * pt(k,j,i)1803 IF ( ql(k,j,i) == 0.0_wp ) THEN 1804 k1 = 1.0_wp + 0.61_wp * q(k,j,i) 1805 k2 = 0.61_wp * pt(k,j,i) 1806 1806 ELSE 1807 1807 theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i) 1808 1808 temp = theta * t_d_pt(k) 1809 k1 = ( 1.0 - q(k,j,i) + 1.61* &1810 ( q(k,j,i) - ql(k,j,i) ) * &1811 ( 1.0 + 0.622* l_d_r / temp ) ) / &1812 ( 1.0 + 0.622* l_d_r * l_d_cp * &1809 k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * & 1810 ( q(k,j,i) - ql(k,j,i) ) * & 1811 ( 1.0_wp + 0.622_wp * l_d_r / temp ) ) / & 1812 ( 1.0_wp + 0.622_wp * l_d_r * l_d_cp * & 1813 1813 ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) ) 1814 k2 = theta * ( l_d_cp / temp * k1 - 1.0 )1814 k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp ) 1815 1815 ENDIF 1816 1816 tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / vpt(k,j,i) * & … … 1819 1819 ) * dd2zu(k) 1820 1820 ELSE IF ( cloud_droplets ) THEN 1821 k1 = 1.0 + 0.61* q(k,j,i) - ql(k,j,i)1822 k2 = 0.61 * pt(k,j,i)1821 k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i) 1822 k2 = 0.61_wp * pt(k,j,i) 1823 1823 tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / vpt(k,j,i) * & 1824 1824 ( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) + & … … 1833 1833 1834 1834 IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN 1835 k1 = 1.0 + 0.61* q(k,j,i)1836 k2 = 0.61 * pt(k,j,i)1835 k1 = 1.0_wp + 0.61_wp * q(k,j,i) 1836 k2 = 0.61_wp * pt(k,j,i) 1837 1837 ELSE IF ( cloud_physics ) THEN 1838 IF ( ql(k,j,i) == 0.0 ) THEN1839 k1 = 1.0 + 0.61* q(k,j,i)1840 k2 = 0.61 * pt(k,j,i)1838 IF ( ql(k,j,i) == 0.0_wp ) THEN 1839 k1 = 1.0_wp + 0.61_wp * q(k,j,i) 1840 k2 = 0.61_wp * pt(k,j,i) 1841 1841 ELSE 1842 1842 theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i) 1843 1843 temp = theta * t_d_pt(k) 1844 k1 = ( 1.0 - q(k,j,i) + 1.61* &1845 ( q(k,j,i) - ql(k,j,i) ) * &1846 ( 1.0 + 0.622* l_d_r / temp ) ) / &1847 ( 1.0 + 0.622* l_d_r * l_d_cp * &1844 k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * & 1845 ( q(k,j,i) - ql(k,j,i) ) * & 1846 ( 1.0_wp + 0.622_wp * l_d_r / temp ) ) / & 1847 ( 1.0_wp + 0.622_wp * l_d_r * l_d_cp * & 1848 1848 ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) ) 1849 k2 = theta * ( l_d_cp / temp * k1 - 1.0 )1849 k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp ) 1850 1850 ENDIF 1851 1851 ELSE IF ( cloud_droplets ) THEN 1852 k1 = 1.0 + 0.61* q(k,j,i) - ql(k,j,i)1853 k2 = 0.61 * pt(k,j,i)1852 k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i) 1853 k2 = 0.61_wp * pt(k,j,i) 1854 1854 ENDIF 1855 1855 … … 1862 1862 1863 1863 IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN 1864 k1 = 1.0 + 0.61* q(k,j,i)1865 k2 = 0.61 * pt(k,j,i)1864 k1 = 1.0_wp + 0.61_wp * q(k,j,i) 1865 k2 = 0.61_wp * pt(k,j,i) 1866 1866 ELSE IF ( cloud_physics ) THEN 1867 IF ( ql(k,j,i) == 0.0 ) THEN1868 k1 = 1.0 + 0.61* q(k,j,i)1869 k2 = 0.61 * pt(k,j,i)1867 IF ( ql(k,j,i) == 0.0_wp ) THEN 1868 k1 = 1.0_wp + 0.61_wp * q(k,j,i) 1869 k2 = 0.61_wp * pt(k,j,i) 1870 1870 ELSE 1871 1871 theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i) 1872 1872 temp = theta * t_d_pt(k) 1873 k1 = ( 1.0 - q(k,j,i) + 1.61* &1874 ( q(k,j,i) - ql(k,j,i) ) * &1875 ( 1.0 + 0.622* l_d_r / temp ) ) / &1876 ( 1.0 + 0.622* l_d_r * l_d_cp * &1873 k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * & 1874 ( q(k,j,i) - ql(k,j,i) ) * & 1875 ( 1.0_wp + 0.622_wp * l_d_r / temp ) ) / & 1876 ( 1.0_wp + 0.622_wp * l_d_r * l_d_cp * & 1877 1877 ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) ) 1878 k2 = theta * ( l_d_cp / temp * k1 - 1.0 )1878 k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp ) 1879 1879 ENDIF 1880 1880 ELSE IF ( cloud_droplets ) THEN 1881 k1 = 1.0 + 0.61* q(k,j,i) - ql(k,j,i)1882 k2 = 0.61 * pt(k,j,i)1881 k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i) 1882 k2 = 0.61_wp * pt(k,j,i) 1883 1883 ENDIF 1884 1884 … … 1917 1917 IF ( first_call ) THEN 1918 1918 ALLOCATE( u_0(nysg:nyng,nxlg:nxrg), v_0(nysg:nyng,nxlg:nxrg) ) 1919 u_0 = 0.0 ! just to avoid access of uninitialized memory1920 v_0 = 0.0 ! within exchange_horiz_2d1919 u_0 = 0.0_wp ! just to avoid access of uninitialized memory 1920 v_0 = 0.0_wp ! within exchange_horiz_2d 1921 1921 first_call = .FALSE. 1922 1922 ENDIF … … 1942 1942 1943 1943 u_0(j,i) = u(ku+1,j,i) + usws(j,i) * ( zu(ku+1) - zu(ku-1) ) / & 1944 ( 0.5 * ( km(ku,j,i) + km(ku,j,i-1) ) +&1945 1.0E-20 )1944 ( 0.5_wp * ( km(ku,j,i) + km(ku,j,i-1) ) + & 1945 1.0E-20_wp ) 1946 1946 ! ( us(j,i) * kappa * zu(1) ) 1947 1947 v_0(j,i) = v(kv+1,j,i) + vsws(j,i) * ( zu(kv+1) - zu(kv-1) ) / & 1948 ( 0.5 * ( km(kv,j,i) + km(kv,j-1,i) ) +&1949 1.0E-20 )1948 ( 0.5_wp * ( km(kv,j,i) + km(kv,j-1,i) ) + & 1949 1.0E-20_wp ) 1950 1950 ! ( us(j,i) * kappa * zu(1) ) 1951 1951
Note: See TracChangeset
for help on using the changeset viewer.
