Changeset 1015 for palm/trunk/SOURCE/diffusion_e.f90
- Timestamp:
- Sep 27, 2012 9:23:24 AM (12 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/diffusion_e.f90 (modified) (12 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/diffusion_e.f90
r1011 r1015 4 4 ! Current revisions: 5 5 ! ----------------- 6 ! 6 ! accelerator version (*_acc) added, 7 ! adjustment of mixing length to the Prandtl mixing length at first grid point 8 ! above ground removed 7 9 ! 8 10 ! Former revisions: … … 56 58 57 59 PRIVATE 58 PUBLIC diffusion_e 60 PUBLIC diffusion_e, diffusion_e_acc 59 61 60 62 … … 64 66 END INTERFACE diffusion_e 65 67 68 INTERFACE diffusion_e_acc 69 MODULE PROCEDURE diffusion_e_acc 70 END INTERFACE diffusion_e_acc 71 66 72 CONTAINS 67 73 … … 81 87 82 88 INTEGER :: i, j, k 83 REAL :: dvar_dz, l_stable, phi_m,var_reference89 REAL :: dvar_dz, l_stable, var_reference 84 90 85 91 #if defined( __nopointer ) … … 98 104 DO i = nxl, nxr 99 105 DO j = nys, nyn 100 !101 !-- First, calculate phi-function for eventually adjusting the &102 !-- mixing length to the prandtl mixing length103 IF ( adjust_mixing_length .AND. prandtl_layer ) THEN104 IF ( rif(j,i) >= 0.0 ) THEN105 phi_m = 1.0 + 5.0 * rif(j,i)106 ELSE107 phi_m = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rif(j,i) ) )108 ENDIF109 ENDIF110 111 106 DO k = nzb_s_inner(j,i)+1, nzt 112 107 ! … … 133 128 ll(k,j) = l_grid(k) 134 129 ENDIF 135 IF ( adjust_mixing_length .AND. prandtl_layer ) THEN136 l(k,j) = MIN( l(k,j), kappa * &137 ( zu(k) - zw(nzb_s_inner(j,i)) ) &138 / phi_m )139 ll(k,j) = MIN( ll(k,j), kappa * &140 ( zu(k) - zw(nzb_s_inner(j,i)) ) &141 / phi_m )142 ENDIF143 130 144 131 ENDDO … … 188 175 DO i = nxl, nxr 189 176 DO j = nys, nyn 190 !191 !-- First, calculate phi-function for eventually adjusting the &192 !-- mixing length to the prandtl mixing length193 IF ( adjust_mixing_length .AND. prandtl_layer ) THEN194 IF ( rif(j,i) >= 0.0 ) THEN195 phi_m = 1.0 + 5.0 * rif(j,i)196 ELSE197 phi_m = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rif(j,i) ) )198 ENDIF199 ENDIF200 201 177 DO k = nzb_s_inner(j,i)+1, nzt 202 178 ! … … 223 199 ll(k,j) = l_grid(k) 224 200 ENDIF 225 IF ( adjust_mixing_length .AND. prandtl_layer ) THEN226 l(k,j) = MIN( l(k,j), kappa * &227 ( zu(k) - zw(nzb_s_inner(j,i)) ) &228 / phi_m )229 ll(k,j) = MIN( ll(k,j), kappa * &230 ( zu(k) - zw(nzb_s_inner(j,i)) ) &231 / phi_m )232 ENDIF233 201 234 202 ENDDO … … 290 258 291 259 !------------------------------------------------------------------------------! 292 ! Call for grid point i,j293 !------------------------------------------------------------------------------! 294 SUBROUTINE diffusion_e_ ij( i, j,var, var_reference )260 ! Call for all grid points - accelerator version 261 !------------------------------------------------------------------------------! 262 SUBROUTINE diffusion_e_acc( var, var_reference ) 295 263 296 264 USE arrays_3d … … 303 271 304 272 INTEGER :: i, j, k 305 REAL :: dvar_dz, l_stable, phi_m, var_reference 273 REAL :: dissipation, dvar_dz, l, ll, l_stable, var_reference 274 275 #if defined( __nopointer ) 276 REAL, DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: var 277 #else 278 REAL, DIMENSION(:,:,:), POINTER :: var 279 #endif 280 281 282 ! 283 !-- This if clause must be outside the k-loop because otherwise 284 !-- runtime errors occur with -C hopt on NEC 285 IF ( use_reference ) THEN 286 STOP '+++ use_reference in diffusion_e not implemented' 287 ! DO i = nxl, nxr 288 ! DO j = nys, nyn 289 ! DO k = nzb_s_inner(j,i)+1, nzt 290 ! 291 !-- Calculate the mixing length (for dissipation) 292 ! dvar_dz = atmos_ocean_sign * & 293 ! ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k) 294 ! IF ( dvar_dz > 0.0 ) THEN 295 ! l_stable = 0.76 * SQRT( e(k,j,i) ) / & 296 ! SQRT( g / var_reference * dvar_dz ) + 1E-5 297 ! ELSE 298 ! l_stable = l_grid(k) 299 ! ENDIF 300 ! 301 !-- Adjustment of the mixing length 302 ! IF ( wall_adjustment ) THEN 303 ! l(k,j) = MIN( wall_adjustment_factor * & 304 ! ( zu(k) - zw(nzb_s_inner(j,i)) ), & 305 ! l_grid(k), l_stable ) 306 ! ll(k,j) = MIN( wall_adjustment_factor * & 307 ! ( zu(k) - zw(nzb_s_inner(j,i)) ), & 308 ! l_grid(k) ) 309 ! ELSE 310 ! l(k,j) = MIN( l_grid(k), l_stable ) 311 ! ll(k,j) = l_grid(k) 312 ! ENDIF 313 ! 314 ! ENDDO 315 ! ENDDO 316 ! 317 ! 318 !-- Calculate the tendency terms 319 ! DO j = nys, nyn 320 ! DO k = nzb_s_inner(j,i)+1, nzt 321 ! 322 ! dissipation(k,j) = ( 0.19 + 0.74 * l(k,j) / ll(k,j) ) * & 323 ! e(k,j,i) * SQRT( e(k,j,i) ) / l(k,j) 324 ! 325 ! tend(k,j,i) = tend(k,j,i) & 326 ! + ( & 327 ! ( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)-e(k,j,i) ) & 328 ! - ( km(k,j,i)+km(k,j,i-1) ) * ( e(k,j,i)-e(k,j,i-1) ) & 329 ! ) * ddx2 & 330 ! + ( & 331 ! ( km(k,j,i)+km(k,j+1,i) ) * ( e(k,j+1,i)-e(k,j,i) ) & 332 ! - ( km(k,j,i)+km(k,j-1,i) ) * ( e(k,j,i)-e(k,j-1,i) ) & 333 ! ) * ddy2 & 334 ! + ( & 335 ! ( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1) & 336 ! - ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k) & 337 ! ) * ddzw(k) & 338 ! - dissipation(k,j) 339 ! 340 ! ENDDO 341 ! ENDDO 342 ! 343 ! 344 !-- Store dissipation if needed for calculating the sgs particle 345 !-- velocities 346 ! IF ( use_sgs_for_particles .OR. wang_kernel ) THEN 347 ! DO j = nys, nyn 348 ! DO k = nzb_s_inner(j,i)+1, nzt 349 ! diss(k,j,i) = dissipation(k,j) 350 ! ENDDO 351 ! ENDDO 352 ! ENDIF 353 ! 354 ! ENDDO 355 ! 356 ELSE 357 358 !$acc kernels present( ddzu, ddzw, dd2zu, diss, e, km, l_grid ) & 359 !$acc present( nzb_s_inner, rif, tend, var, zu, zw ) 360 !$acc loop 361 DO i = nxl, nxr 362 DO j = nys, nyn 363 !$acc loop vector( 32 ) 364 DO k = 1, nzt 365 366 IF ( k > nzb_s_inner(j,i) ) THEN 367 ! 368 !-- Calculate the mixing length (for dissipation) 369 dvar_dz = atmos_ocean_sign * & 370 ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k) 371 IF ( dvar_dz > 0.0 ) THEN 372 l_stable = 0.76 * SQRT( e(k,j,i) ) / & 373 SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5 374 ELSE 375 l_stable = l_grid(k) 376 ENDIF 377 ! 378 !-- Adjustment of the mixing length 379 IF ( wall_adjustment ) THEN 380 l = MIN( wall_adjustment_factor * & 381 ( zu(k) - zw(nzb_s_inner(j,i)) ), & 382 l_grid(k), l_stable ) 383 ll = MIN( wall_adjustment_factor * & 384 ( zu(k) - zw(nzb_s_inner(j,i)) ), & 385 l_grid(k) ) 386 ELSE 387 l = MIN( l_grid(k), l_stable ) 388 ll = l_grid(k) 389 ENDIF 390 ! 391 !-- Calculate the tendency terms 392 dissipation = ( 0.19 + 0.74 * l / ll ) * & 393 e(k,j,i) * SQRT( e(k,j,i) ) / l 394 395 tend(k,j,i) = tend(k,j,i) & 396 + ( & 397 ( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)-e(k,j,i) ) & 398 - ( km(k,j,i)+km(k,j,i-1) ) * ( e(k,j,i)-e(k,j,i-1) ) & 399 ) * ddx2 & 400 + ( & 401 ( km(k,j,i)+km(k,j+1,i) ) * ( e(k,j+1,i)-e(k,j,i) ) & 402 - ( km(k,j,i)+km(k,j-1,i) ) * ( e(k,j,i)-e(k,j-1,i) ) & 403 ) * ddy2 & 404 + ( & 405 ( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1) & 406 - ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k) & 407 ) * ddzw(k) & 408 - dissipation 409 410 ! 411 !-- Store dissipation if needed for calculating the sgs 412 !-- particle velocities 413 IF ( use_sgs_for_particles .OR. wang_kernel ) THEN 414 diss(k,j,i) = dissipation 415 ENDIF 416 417 ENDIF 418 419 ENDDO 420 ENDDO 421 ENDDO 422 !$acc end kernels 423 424 ENDIF 425 426 ! 427 !-- Boundary condition for dissipation 428 IF ( use_sgs_for_particles .OR. wang_kernel ) THEN 429 !$acc kernels present( diss, nzb_s_inner ) 430 !$acc loop 431 DO i = nxl, nxr 432 !$acc loop vector( 32 ) 433 DO j = nys, nyn 434 diss(nzb_s_inner(j,i),j,i) = diss(nzb_s_inner(j,i)+1,j,i) 435 ENDDO 436 ENDDO 437 !$acc end kernels 438 ENDIF 439 440 END SUBROUTINE diffusion_e_acc 441 442 443 !------------------------------------------------------------------------------! 444 ! Call for grid point i,j 445 !------------------------------------------------------------------------------! 446 SUBROUTINE diffusion_e_ij( i, j, var, var_reference ) 447 448 USE arrays_3d 449 USE control_parameters 450 USE grid_variables 451 USE indices 452 USE particle_attributes 453 454 IMPLICIT NONE 455 456 INTEGER :: i, j, k 457 REAL :: dvar_dz, l_stable, var_reference 306 458 307 459 #if defined( __nopointer ) … … 312 464 REAL, DIMENSION(nzb+1:nzt) :: dissipation, l, ll 313 465 314 315 !316 !-- First, calculate phi-function for eventually adjusting the mixing length317 !-- to the prandtl mixing length318 IF ( adjust_mixing_length .AND. prandtl_layer ) THEN319 IF ( rif(j,i) >= 0.0 ) THEN320 phi_m = 1.0 + 5.0 * rif(j,i)321 ELSE322 phi_m = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rif(j,i) ) )323 ENDIF324 ENDIF325 466 326 467 ! … … 352 493 ll(k) = l_grid(k) 353 494 ENDIF 354 IF ( adjust_mixing_length .AND. prandtl_layer ) THEN355 l(k) = MIN( l(k), kappa * &356 ( zu(k) - zw(nzb_s_inner(j,i)) ) / phi_m )357 ll(k) = MIN( ll(k), kappa * &358 ( zu(k) - zw(nzb_s_inner(j,i)) ) / phi_m )359 ENDIF360 361 495 ! 362 496 !-- Calculate the tendency term
Note: See TracChangeset
for help on using the changeset viewer.
