Changeset 2232 for palm/trunk/SOURCE/diffusion_e.f90
- Timestamp:
- May 30, 2017 5:47:52 PM (7 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/diffusion_e.f90 (modified) (22 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/diffusion_e.f90
r2119 r2232 20 20 ! Current revisions: 21 21 ! ----------------- 22 ! 22 ! Adjustments to new topography and surface concept 23 23 ! 24 24 ! Former revisions: … … 129 129 130 130 USE arrays_3d, & 131 ONLY: dd2zu, ddzu, ddzw, diss, e, km, l_grid, tend, zu, zw,&131 ONLY: dd2zu, ddzu, ddzw, diss, e, km, l_grid, l_wall, tend, zu, zw,& 132 132 drho_air, rho_air_zw 133 133 … … 140 140 141 141 USE indices, & 142 ONLY: nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, nzb_ s_inner,&143 nzt 142 ONLY: nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, nzb_max, & 143 nzt, wall_flags_0 144 144 145 145 USE kinds … … 151 151 ONLY: use_sgs_for_particles, wang_kernel 152 152 153 USE surface_mod, & 154 ONLY : bc_h 155 153 156 IMPLICIT NONE 154 157 155 INTEGER(iwp) :: i !< 156 INTEGER(iwp) :: j !< 157 INTEGER(iwp) :: k !< 158 INTEGER(iwp) :: i !< running index x direction 159 INTEGER(iwp) :: j !< running index y direction 160 INTEGER(iwp) :: k !< running index z direction 161 INTEGER(iwp) :: m !< running index surface elements 162 158 163 REAL(wp) :: dvar_dz !< 164 REAL(wp) :: flag !< flag to mask topography 159 165 REAL(wp) :: l_stable !< 160 166 REAL(wp) :: var_reference !< … … 177 183 DO i = nxl, nxr 178 184 DO j = nys, nyn 179 DO k = nzb _s_inner(j,i)+1, nzt185 DO k = nzb+1, nzt 180 186 ! 181 187 !-- Calculate the mixing length (for dissipation) … … 191 197 !-- Adjustment of the mixing length 192 198 IF ( wall_adjustment ) THEN 193 l(k,j) = MIN( wall_adjustment_factor * & 194 ( zu(k) - zw(nzb_s_inner(j,i)) ), & 199 l(k,j) = MIN( wall_adjustment_factor * l_wall(k,j,i), & 195 200 l_grid(k), l_stable ) 196 ll(k,j) = MIN( wall_adjustment_factor * & 197 ( zu(k) - zw(nzb_s_inner(j,i)) ), & 201 ll(k,j) = MIN( wall_adjustment_factor * l_wall(k,j,i), & 198 202 l_grid(k) ) 199 203 ELSE … … 208 212 !-- Calculate the tendency terms 209 213 DO j = nys, nyn 210 DO k = nzb_s_inner(j,i)+1, nzt 211 212 dissipation(k,j) = ( 0.19_wp + 0.74_wp * l(k,j) / ll(k,j) ) * & 214 DO k = nzb+1, nzt 215 ! 216 !-- Predetermine flag to mask topography 217 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) ) 218 219 dissipation(k,j) = ( 0.19_wp + 0.74_wp * l(k,j) / ll(k,j) ) * & 213 220 e(k,j,i) * SQRT( e(k,j,i) ) / l(k,j) 214 221 215 tend(k,j,i) = tend(k,j,i)&222 tend(k,j,i) = tend(k,j,i) & 216 223 + ( & 217 224 ( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)-e(k,j,i) ) & 218 225 - ( km(k,j,i)+km(k,j,i-1) ) * ( e(k,j,i)-e(k,j,i-1) ) & 219 ) * ddx2 226 ) * ddx2 * flag & 220 227 + ( & 221 228 ( km(k,j,i)+km(k,j+1,i) ) * ( e(k,j+1,i)-e(k,j,i) ) & 222 229 - ( km(k,j,i)+km(k,j-1,i) ) * ( e(k,j,i)-e(k,j-1,i) ) & 223 ) * ddy2 230 ) * ddy2 * flag & 224 231 + ( & 225 232 ( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1) & … … 227 234 - ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k) & 228 235 * rho_air_zw(k-1) & 229 ) * ddzw(k) * drho_air(k) 230 - dissipation(k,j) 236 ) * ddzw(k) * drho_air(k) * flag & 237 - dissipation(k,j) * flag 231 238 232 239 ENDDO … … 239 246 collision_turbulence ) THEN 240 247 DO j = nys, nyn 241 DO k = nzb_s_inner(j,i)+1, nzt 242 diss(k,j,i) = dissipation(k,j) 248 DO k = nzb+1, nzt 249 diss(k,j,i) = dissipation(k,j) * & 250 MERGE( 1.0_wp, 0.0_wp, & 251 BTEST( wall_flags_0(k,j,i), 0 ) ) 243 252 ENDDO 244 253 ENDDO … … 251 260 DO i = nxl, nxr 252 261 DO j = nys, nyn 253 DO k = nzb _s_inner(j,i)+1, nzt262 DO k = nzb+1, nzt 254 263 ! 255 264 !-- Calculate the mixing length (for dissipation) … … 265 274 !-- Adjustment of the mixing length 266 275 IF ( wall_adjustment ) THEN 267 l(k,j) = MIN( wall_adjustment_factor * & 268 ( zu(k) - zw(nzb_s_inner(j,i)) ), & 276 l(k,j) = MIN( wall_adjustment_factor * l_wall(k,j,i), & 269 277 l_grid(k), l_stable ) 270 ll(k,j) = MIN( wall_adjustment_factor * & 271 ( zu(k) - zw(nzb_s_inner(j,i)) ), & 278 ll(k,j) = MIN( wall_adjustment_factor * l_wall(k,j,i), & 272 279 l_grid(k) ) 273 280 ELSE … … 282 289 !-- Calculate the tendency terms 283 290 DO j = nys, nyn 284 DO k = nzb_s_inner(j,i)+1, nzt 285 286 dissipation(k,j) = ( 0.19_wp + 0.74_wp * l(k,j) / ll(k,j) ) * & 291 DO k = nzb+1, nzt 292 ! 293 !-- Predetermine flag to mask topography 294 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) ) 295 296 dissipation(k,j) = ( 0.19_wp + 0.74_wp * l(k,j) / ll(k,j) ) *& 287 297 e(k,j,i) * SQRT( e(k,j,i) ) / l(k,j) 288 298 289 tend(k,j,i) = tend(k,j,i)&299 tend(k,j,i) = tend(k,j,i) & 290 300 + ( & 291 301 ( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)-e(k,j,i) ) & 292 302 - ( km(k,j,i)+km(k,j,i-1) ) * ( e(k,j,i)-e(k,j,i-1) ) & 293 ) * ddx2 303 ) * ddx2 * flag & 294 304 + ( & 295 305 ( km(k,j,i)+km(k,j+1,i) ) * ( e(k,j+1,i)-e(k,j,i) ) & 296 306 - ( km(k,j,i)+km(k,j-1,i) ) * ( e(k,j,i)-e(k,j-1,i) ) & 297 ) * ddy2 307 ) * ddy2 * flag & 298 308 + ( & 299 309 ( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1) & … … 301 311 - ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k) & 302 312 * rho_air_zw(k-1) & 303 ) * ddzw(k) * drho_air(k) 304 - dissipation(k,j) 313 ) * ddzw(k) * drho_air(k) * flag & 314 - dissipation(k,j) * flag 305 315 306 316 ENDDO … … 313 323 collision_turbulence ) THEN 314 324 DO j = nys, nyn 315 DO k = nzb_s_inner(j,i)+1, nzt 316 diss(k,j,i) = dissipation(k,j) 325 DO k = nzb+1, nzt 326 diss(k,j,i) = dissipation(k,j) * & 327 MERGE( 1.0_wp, 0.0_wp, & 328 BTEST( wall_flags_0(k,j,i), 0 ) ) 317 329 ENDDO 318 330 ENDDO … … 324 336 325 337 ! 326 !-- Boundary condition for dissipation338 !-- Neumann boundary condition for dissipation diss(nzb,:,:) = diss(nzb+1,:,:) 327 339 IF ( use_sgs_for_particles .OR. wang_kernel .OR. & 328 340 collision_turbulence ) THEN 329 DO i = nxl, nxr 330 DO j = nys, nyn 331 diss(nzb_s_inner(j,i),j,i) = diss(nzb_s_inner(j,i)+1,j,i) 332 ENDDO 333 ENDDO 341 ! 342 !-- Upward facing surfaces 343 DO m = 1, bc_h(0)%ns 344 i = bc_h(0)%i(m) 345 j = bc_h(0)%j(m) 346 k = bc_h(0)%k(m) 347 diss(k-1,j,i) = diss(k,j,i) 348 ENDDO 349 ! 350 !-- Downward facing surfaces 351 DO m = 1, bc_h(1)%ns 352 i = bc_h(1)%i(m) 353 j = bc_h(1)%j(m) 354 k = bc_h(1)%k(m) 355 diss(k+1,j,i) = diss(k,j,i) 356 ENDDO 357 334 358 ENDIF 335 359 … … 345 369 346 370 USE arrays_3d, & 347 ONLY: dd2zu, ddzu, ddzw, diss, e, km, l_grid, tend, zu, zw,&371 ONLY: dd2zu, ddzu, ddzw, diss, e, km, l_grid, l_wall, tend, zu, zw,& 348 372 drho_air, rho_air_zw 349 373 … … 356 380 357 381 USE indices, & 358 ONLY: nxlg, nxrg, nyng, nysg, nzb, nzb_ s_inner, nzt382 ONLY: nxlg, nxrg, nyng, nysg, nzb, nzb_max, nzt, wall_flags_0 359 383 360 384 USE kinds … … 366 390 ONLY: use_sgs_for_particles, wang_kernel 367 391 392 USE surface_mod, & 393 ONLY : bc_h 394 368 395 IMPLICIT NONE 369 396 370 INTEGER(iwp) :: i !< 371 INTEGER(iwp) :: j !< 372 INTEGER(iwp) :: k !< 397 INTEGER(iwp) :: i !< running index x direction 398 INTEGER(iwp) :: j !< running index y direction 399 INTEGER(iwp) :: k !< running index z direction 400 INTEGER(iwp) :: m !< running index surface elements 401 INTEGER(iwp) :: surf_e !< End index of surface elements at (j,i)-gridpoint 402 INTEGER(iwp) :: surf_s !< Start index of surface elements at (j,i)-gridpoint 403 373 404 REAL(wp) :: dvar_dz !< 405 REAL(wp) :: flag !< flag to mask topography 374 406 REAL(wp) :: l_stable !< 375 407 REAL(wp) :: var_reference !< … … 384 416 REAL(wp), DIMENSION(nzb+1:nzt) :: ll !< 385 417 386 387 418 ! 388 419 !-- Calculate the mixing length (for dissipation) 389 DO k = nzb_s_inner(j,i)+1, nzt 420 DO k = nzb+1, nzt 421 ! 422 !-- Predetermine flag to mask topography 423 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) ) 424 390 425 dvar_dz = atmos_ocean_sign * & 391 426 ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k) … … 404 439 !-- Adjustment of the mixing length 405 440 IF ( wall_adjustment ) THEN 406 l(k) = MIN( wall_adjustment_factor * & 407 ( zu(k) - zw(nzb_s_inner(j,i)) ), l_grid(k), & 408 l_stable ) 409 ll(k) = MIN( wall_adjustment_factor * & 410 ( zu(k) - zw(nzb_s_inner(j,i)) ), l_grid(k) ) 441 l(k) = MIN( wall_adjustment_factor * l_wall(k,j,i), & 442 l_grid(k), l_stable ) 443 ll(k) = MIN( wall_adjustment_factor * l_wall(k,j,i), l_grid(k) ) 411 444 ELSE 412 445 l(k) = MIN( l_grid(k), l_stable ) … … 415 448 ! 416 449 !-- Calculate the tendency term 417 dissipation(k) = ( 0.19_wp + 0.74_wp * l(k) / ll(k) ) * e(k,j,i) * &450 dissipation(k) = ( 0.19_wp + 0.74_wp * l(k) / ll(k) ) * e(k,j,i) * & 418 451 SQRT( e(k,j,i) ) / l(k) 419 452 420 tend(k,j,i) = tend(k,j,i) &421 + ( &422 ( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)-e(k,j,i) ) &423 - ( km(k,j,i)+km(k,j,i-1) ) * ( e(k,j,i)-e(k,j,i-1) ) &424 ) * ddx2 425 + ( &426 ( km(k,j,i)+km(k,j+1,i) ) * ( e(k,j+1,i)-e(k,j,i) ) &427 - ( km(k,j,i)+km(k,j-1,i) ) * ( e(k,j,i)-e(k,j-1,i) ) &428 ) * ddy2 429 + ( &430 ( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1) &431 * rho_air_zw(k) &432 - ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k) &433 * rho_air_zw(k-1) &434 ) * ddzw(k) * drho_air(k) 435 - dissipation(k) 453 tend(k,j,i) = tend(k,j,i) & 454 + ( & 455 ( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)-e(k,j,i) ) & 456 - ( km(k,j,i)+km(k,j,i-1) ) * ( e(k,j,i)-e(k,j,i-1) ) & 457 ) * ddx2 * flag & 458 + ( & 459 ( km(k,j,i)+km(k,j+1,i) ) * ( e(k,j+1,i)-e(k,j,i) ) & 460 - ( km(k,j,i)+km(k,j-1,i) ) * ( e(k,j,i)-e(k,j-1,i) ) & 461 ) * ddy2 * flag & 462 + ( & 463 ( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1) & 464 * rho_air_zw(k) & 465 - ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k) & 466 * rho_air_zw(k-1) & 467 ) * ddzw(k) * drho_air(k) * flag & 468 - dissipation(k) * flag 436 469 437 470 ENDDO … … 441 474 IF ( use_sgs_for_particles .OR. wang_kernel .OR. & 442 475 collision_turbulence ) THEN 443 DO k = nzb_s_inner(j,i)+1, nzt 444 diss(k,j,i) = dissipation(k) 445 ENDDO 446 ! 447 !-- Boundary condition for dissipation 448 diss(nzb_s_inner(j,i),j,i) = diss(nzb_s_inner(j,i)+1,j,i) 476 DO k = nzb+1, nzt 477 diss(k,j,i) = dissipation(k) * & 478 MERGE( 1.0_wp, 0.0_wp, & 479 BTEST( wall_flags_0(k,j,i), 0 ) ) 480 ENDDO 481 ! 482 !-- Neumann boundary condition for dissipation diss(nzb,:,:) = diss(nzb+1,:,:) 483 !-- For each surface type determine start and end index (in case of elevated 484 !-- toopography several up/downward facing surfaces may exist. 485 surf_s = bc_h(0)%start_index(j,i) 486 surf_e = bc_h(0)%end_index(j,i) 487 DO m = surf_s, surf_e 488 k = bc_h(0)%k(m) 489 diss(k-1,j,i) = diss(k,j,i) 490 ENDDO 491 ! 492 !-- Downward facing surfaces 493 surf_s = bc_h(1)%start_index(j,i) 494 surf_e = bc_h(1)%end_index(j,i) 495 DO m = surf_s, surf_e 496 k = bc_h(1)%k(m) 497 diss(k+1,j,i) = diss(k,j,i) 498 ENDDO 449 499 ENDIF 450 500
Note: See TracChangeset
for help on using the changeset viewer.
