Changeset 3634 for palm/trunk/SOURCE/turbulence_closure_mod.f90
 Timestamp:
 Dec 18, 2018 12:31:28 PM (3 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

palm/trunk/SOURCE/turbulence_closure_mod.f90
r3550 r3634 25 25 !  26 26 ! $Id$ 27 ! OpenACC port for SPEC 28 ! 29 ! 3550 20181121 16:01:01Z gronemeier 27 30 !  calculate diss production same in vector and cache optimization 28 31 !  move boundary condition for e and diss to boundary_conds … … 1619 1622 ENDDO !k loop 1620 1623 1624 !$ACC ENTER DATA COPYIN(l_black(nzb:nzt+1)) 1625 1621 1626 ENDIF !LES or RANS mode 1622 1627 … … 1624 1629 ! Set lateral boundary conditions for l_wall 1625 1630 CALL exchange_horiz( l_wall, nbgp ) 1631 1632 !$ACC ENTER DATA COPYIN(l_grid(nzb:nzt+1)) & 1633 !$ACC COPYIN(l_wall(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) 1626 1634 1627 1635 CONTAINS … … 1777 1785 ! Default surfaces, upwardfacing 1778 1786 !$OMP PARALLEL DO PRIVATE(i,j,k,m) 1787 !$ACC PARALLEL LOOP PRIVATE(i, j, k, m, km_sfc) & 1788 !$ACC PRESENT(surf_def_h(0), u, v, drho_air_zw, zu) 1779 1789 DO m = 1, surf_def_h(0)%ns 1780 1790 … … 1811 1821 ! Default surfaces, downwardfacing surfaces 1812 1822 !$OMP PARALLEL DO PRIVATE(i,j,k,m) 1823 !$ACC PARALLEL LOOP PRIVATE(i, j, k, m, km_sfc) & 1824 !$ACC PRESENT(surf_def_h(1), u, v, drho_air_zw, zu, km) 1813 1825 DO m = 1, surf_def_h(1)%ns 1814 1826 … … 1842 1854 ! Natural surfaces, upwardfacing 1843 1855 !$OMP PARALLEL DO PRIVATE(i,j,k,m) 1856 !$ACC PARALLEL LOOP PRIVATE(i, j, k, m, km_sfc) & 1857 !$ACC PRESENT(surf_lsm_h, u, v, drho_air_zw, zu) 1844 1858 DO m = 1, surf_lsm_h%ns 1845 1859 … … 1876 1890 ! Urban surfaces, upwardfacing 1877 1891 !$OMP PARALLEL DO PRIVATE(i,j,k,m) 1892 !$ACC PARALLEL LOOP PRIVATE(i, j, k, m, km_sfc) & 1893 !$ACC PRESENT(surf_usm_h, u, v, drho_air_zw, zu) 1878 1894 DO m = 1, surf_usm_h%ns 1879 1895 … … 1970 1986 CALL advec_s_up( e ) 1971 1987 ELSE 1988 !$ACC KERNELS PRESENT(tend) 1972 1989 tend = 0.0_wp 1990 !$ACC END KERNELS 1973 1991 IF ( timestep_scheme(1:5) == 'runge' ) THEN 1974 1992 IF ( ws_scheme_sca ) THEN … … 2006 2024 ! reasons in the course of the integration. In such cases the old TKE 2007 2025 ! value is reduced by 90%. 2026 !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) & 2027 !$ACC PRESENT(e, tend, te_m, wall_flags_0) & 2028 !$ACC PRESENT(tsc(3:3)) & 2029 !$ACC PRESENT(e_p) 2008 2030 DO i = nxl, nxr 2009 2031 DO j = nys, nyn … … 2024 2046 IF ( timestep_scheme(1:5) == 'runge' ) THEN 2025 2047 IF ( intermediate_timestep_count == 1 ) THEN 2048 !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) & 2049 !$ACC PRESENT(tend, te_m) 2026 2050 DO i = nxl, nxr 2027 2051 DO j = nys, nyn … … 2033 2057 ELSEIF ( intermediate_timestep_count < & 2034 2058 intermediate_timestep_count_max ) THEN 2059 !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) & 2060 !$ACC PRESENT(tend, te_m) 2035 2061 DO i = nxl, nxr 2036 2062 DO j = nys, nyn … … 2380 2406 ! points first, gradients at surfacebounded grid points will be 2381 2407 ! overwritten further below. 2408 !$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j, l) & 2409 !$ACC PRIVATE(surf_s, surf_e) & 2410 !$ACC PRIVATE(dudx(:), dudy(:), dudz(:), dvdx(:), dvdy(:), dvdz(:), dwdx(:), dwdy(:), dwdz(:)) & 2411 !$ACC PRESENT(e, u, v, w, diss, dd2zu, ddzw, km, wall_flags_0) & 2412 !$ACC PRESENT(tend) & 2413 !$ACC PRESENT(surf_def_h(0:1), surf_def_v(0:3)) & 2414 !$ACC PRESENT(surf_lsm_h, surf_lsm_v(0:3)) & 2415 !$ACC PRESENT(surf_usm_h, surf_usm_v(0:3)) 2382 2416 DO i = nxl, nxr 2383 2417 DO j = nys, nyn 2418 !$ACC LOOP PRIVATE(k) 2384 2419 DO k = nzb+1, nzt 2385 2420 … … 2423 2458 surf_s = surf_def_v(l)%start_index(j,i) 2424 2459 surf_e = surf_def_v(l)%end_index(j,i) 2460 !$ACC LOOP PRIVATE(m, k, usvs, wsvs, km_neutral, sign_dir) 2425 2461 DO m = surf_s, surf_e 2426 2462 k = surf_def_v(l)%k(m) … … 2441 2477 surf_s = surf_lsm_v(l)%start_index(j,i) 2442 2478 surf_e = surf_lsm_v(l)%end_index(j,i) 2479 !$ACC LOOP PRIVATE(m, k, usvs, wsvs, km_neutral, sign_dir) 2443 2480 DO m = surf_s, surf_e 2444 2481 k = surf_lsm_v(l)%k(m) … … 2459 2496 surf_s = surf_usm_v(l)%start_index(j,i) 2460 2497 surf_e = surf_usm_v(l)%end_index(j,i) 2498 !$ACC LOOP PRIVATE(m, k, usvs, wsvs, km_neutral, sign_dir) 2461 2499 DO m = surf_s, surf_e 2462 2500 k = surf_usm_v(l)%k(m) … … 2479 2517 surf_s = surf_def_v(l)%start_index(j,i) 2480 2518 surf_e = surf_def_v(l)%end_index(j,i) 2519 !$ACC LOOP PRIVATE(m, k, vsus, wsus, km_neutral, sign_dir) 2481 2520 DO m = surf_s, surf_e 2482 2521 k = surf_def_v(l)%k(m) … … 2497 2536 surf_s = surf_lsm_v(l)%start_index(j,i) 2498 2537 surf_e = surf_lsm_v(l)%end_index(j,i) 2538 !$ACC LOOP PRIVATE(m, k, vsus, wsus, km_neutral, sign_dir) 2499 2539 DO m = surf_s, surf_e 2500 2540 k = surf_lsm_v(l)%k(m) … … 2515 2555 surf_s = surf_usm_v(l)%start_index(j,i) 2516 2556 surf_e = surf_usm_v(l)%end_index(j,i) 2557 !$ACC LOOP PRIVATE(m, k, vsus, wsus, km_neutral, sign_dir) 2517 2558 DO m = surf_s, surf_e 2518 2559 k = surf_usm_v(l)%k(m) … … 2534 2575 surf_s = surf_def_h(0)%start_index(j,i) 2535 2576 surf_e = surf_def_h(0)%end_index(j,i) 2577 !$ACC LOOP PRIVATE(m, k) 2536 2578 DO m = surf_s, surf_e 2537 2579 k = surf_def_h(0)%k(m) … … 2550 2592 surf_s = surf_lsm_h%start_index(j,i) 2551 2593 surf_e = surf_lsm_h%end_index(j,i) 2594 !$ACC LOOP PRIVATE(m, k) 2552 2595 DO m = surf_s, surf_e 2553 2596 k = surf_lsm_h%k(m) … … 2561 2604 surf_s = surf_usm_h%start_index(j,i) 2562 2605 surf_e = surf_usm_h%end_index(j,i) 2606 !$ACC LOOP PRIVATE(m, k) 2563 2607 DO m = surf_s, surf_e 2564 2608 k = surf_usm_h%k(m) … … 2573 2617 surf_s = surf_def_h(1)%start_index(j,i) 2574 2618 surf_e = surf_def_h(1)%end_index(j,i) 2619 !$ACC LOOP PRIVATE(m, k) 2575 2620 DO m = surf_s, surf_e 2576 2621 k = surf_def_h(1)%k(m) … … 2585 2630 2586 2631 2632 !$ACC LOOP PRIVATE(k, def, flag) 2587 2633 DO k = nzb+1, nzt 2588 2634 … … 2684 2730 ELSE ! or IF ( .NOT. ocean_mode ) THEN 2685 2731 2732 !$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j) & 2733 !$ACC PRIVATE(surf_s, surf_e) & 2734 !$ACC PRIVATE(tmp_flux(nzb+1:nzt)) & 2735 !$ACC PRESENT(e, diss, kh, pt, dd2zu, drho_air_zw, wall_flags_0) & 2736 !$ACC PRESENT(tend) & 2737 !$ACC PRESENT(surf_def_h(0:2)) & 2738 !$ACC PRESENT(surf_lsm_h) & 2739 !$ACC PRESENT(surf_usm_h) 2686 2740 DO i = nxl, nxr 2687 2741 DO j = nys, nyn 2688 2742 2743 !$ACC LOOP PRIVATE(k) 2689 2744 DO k = nzb+1, nzt 2690 2745 tmp_flux(k) = 1.0_wp * kh(k,j,i) * ( pt(k+1,j,i)  pt(k1,j,i) ) * dd2zu(k) … … 2697 2752 surf_s = surf_def_h(l)%start_index(j,i) 2698 2753 surf_e = surf_def_h(l)%end_index(j,i) 2754 !$ACC LOOP PRIVATE(m, k) 2699 2755 DO m = surf_s, surf_e 2700 2756 k = surf_def_h(l)%k(m) … … 2706 2762 surf_s = surf_lsm_h%start_index(j,i) 2707 2763 surf_e = surf_lsm_h%end_index(j,i) 2764 !$ACC LOOP PRIVATE(m, k) 2708 2765 DO m = surf_s, surf_e 2709 2766 k = surf_lsm_h%k(m) … … 2714 2771 surf_s = surf_usm_h%start_index(j,i) 2715 2772 surf_e = surf_usm_h%end_index(j,i) 2773 !$ACC LOOP PRIVATE(m, k) 2716 2774 DO m = surf_s, surf_e 2717 2775 k = surf_usm_h%k(m) … … 2723 2781 surf_s = surf_def_h(2)%start_index(j,i) 2724 2782 surf_e = surf_def_h(2)%end_index(j,i) 2783 !$ACC LOOP PRIVATE(m, k) 2725 2784 DO m = surf_s, surf_e 2726 2785 k = surf_def_h(2)%k(m) … … 2732 2791 2733 2792 ! Compute tendency for TKEproduction from shear 2793 !$ACC LOOP PRIVATE(k, flag) 2734 2794 DO k = nzb+1, nzt 2735 2795 flag = MERGE( 1.0_wp, 0.0_wp, BTEST(wall_flags_0(k,j,i),0) ) … … 3628 3688 3629 3689 USE arrays_3d, & 3630 ONLY: ddzu, ddzw, drho_air, rho_air_zw 3690 ONLY: ddzu, dd2zu, ddzw, drho_air, rho_air_zw 3691 3692 USE control_parameters, & 3693 ONLY: atmos_ocean_sign, use_single_reference_value, & 3694 wall_adjustment, wall_adjustment_factor 3631 3695 3632 3696 USE grid_variables, & … … 3653 3717 REAL(wp) :: ll !< adjusted l 3654 3718 REAL(wp) :: var_reference !< reference temperature 3719 #ifdef _OPENACC 3720 ! 3721 ! From mixing_length_les: 3722 REAL(wp) :: l_stable !< mixing length according to stratification 3723 ! 3724 ! From mixing_length_rans: 3725 REAL(wp) :: duv2_dz2 !< squared vertical gradient of wind vector 3726 REAL(wp) :: l_diss !< mixing length for dissipation 3727 REAL(wp) :: rif !< Richardson flux number 3728 ! 3729 ! From both: 3730 REAL(wp) :: dvar_dz !< vertical gradient of var 3731 #endif 3655 3732 3656 3733 #if defined( __nopointer ) … … 3659 3736 REAL(wp), DIMENSION(:,:,:), POINTER :: var !< temperature 3660 3737 #endif 3661 REAL(wp) , DIMENSION(nzb+1:nzt,nys:nyn):: dissipation !< TKE dissipation3738 REAL(wp) :: dissipation !< TKE dissipation 3662 3739 3663 3740 3664 3741 ! 3665 3742 ! Calculate the tendency terms 3743 !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) & 3744 !$ACC PRIVATE(flag, l, ll, dissipation) & 3745 !$ACC PRIVATE(l_stable, duv2_dz2, l_diss, rif, dvar_dz) & 3746 !$ACC PRESENT(e, u, v, km, var, wall_flags_0) & 3747 !$ACC PRESENT(ddzu, dd2zu, ddzw, rho_air_zw, drho_air) & 3748 !$ACC PRESENT(l_black, l_grid, l_wall) & 3749 !$ACC PRESENT(diss, tend) 3666 3750 DO i = nxl, nxr 3667 3751 DO j = nys, nyn … … 3675 3759 IF ( les_dynamic .OR. les_mw ) THEN 3676 3760 3761 #ifdef _OPENACC 3762 ! 3763 ! Cannot call subroutine mixing_length_les because after adding all required 3764 ! OpenACC directives the execution crashed reliably when inside the called 3765 ! subroutine. I'm not sure why that is... 3766 dvar_dz = atmos_ocean_sign * ( var(k+1,j,i)  var(k1,j,i) ) * dd2zu(k) 3767 IF ( dvar_dz > 0.0_wp ) THEN 3768 IF ( use_single_reference_value ) THEN 3769 l_stable = 0.76_wp * SQRT( e(k,j,i) ) & 3770 / SQRT( g / var_reference * dvar_dz ) + 1E5_wp 3771 ELSE 3772 l_stable = 0.76_wp * SQRT( e(k,j,i) ) & 3773 / SQRT( g / var(k,j,i) * dvar_dz ) + 1E5_wp 3774 ENDIF 3775 ELSE 3776 l_stable = l_grid(k) 3777 ENDIF 3778 ! 3779 ! Adjustment of the mixing length 3780 IF ( wall_adjustment ) THEN 3781 l = MIN( wall_adjustment_factor * l_wall(k,j,i), l_grid(k), l_stable ) 3782 ll = MIN( wall_adjustment_factor * l_wall(k,j,i), l_grid(k) ) 3783 ELSE 3784 l = MIN( l_grid(k), l_stable ) 3785 ll = l_grid(k) 3786 ENDIF 3787 #else 3677 3788 CALL mixing_length_les( i, j, k, l, ll, var, var_reference ) 3678 3679 dissipation(k,j) = ( 0.19_wp + 0.74_wp * l / ll ) & 3680 * e(k,j,i) * SQRT( e(k,j,i) ) / l 3789 #endif 3790 3791 dissipation = ( 0.19_wp + 0.74_wp * l / ll ) & 3792 * e(k,j,i) * SQRT( e(k,j,i) ) / l 3681 3793 3682 3794 ELSEIF ( rans_tke_l ) THEN 3683 3795 3796 #ifdef _OPENACC 3797 ! 3798 ! Same reason as above... 3799 dvar_dz = atmos_ocean_sign * ( var(k+1,j,i)  var(k1,j,i) ) * dd2zu(k) 3800 3801 duv2_dz2 = ( ( u(k+1,j,i)  u(k1,j,i) ) * dd2zu(k) )**2 & 3802 + ( ( v(k+1,j,i)  v(k1,j,i) ) * dd2zu(k) )**2 & 3803 + 1E30_wp 3804 3805 IF ( use_single_reference_value ) THEN 3806 rif = g / var_reference * dvar_dz / duv2_dz2 3807 ELSE 3808 rif = g / var(k,j,i) * dvar_dz / duv2_dz2 3809 ENDIF 3810 3811 rif = MAX( rif, 5.0_wp ) 3812 rif = MIN( rif, 1.0_wp ) 3813 3814 ! 3815 ! Calculate diabatic mixing length using Dyerprofile functions 3816 IF ( rif >= 0.0_wp ) THEN 3817 l = MIN( l_black(k) / ( 1.0_wp + 5.0_wp * rif ), l_wall(k,j,i) ) 3818 l_diss = l 3819 ELSE 3820 ! 3821 ! In case of unstable stratification, use mixing length of neutral case 3822 ! for l, but consider profile functions for l_diss 3823 l = l_wall(k,j,i) 3824 l_diss = l * SQRT( 1.0_wp  16.0_wp * rif ) 3825 ENDIF 3826 #else 3684 3827 CALL mixing_length_rans( i, j, k, l, ll, var, var_reference ) 3685 3686 dissipation(k,j) = c_0**3 * e(k,j,i) * SQRT( e(k,j,i) ) / ll 3687 3688 diss(k,j,i) = dissipation(k,j) * flag 3828 #endif 3829 3830 dissipation = c_0**3 * e(k,j,i) * SQRT( e(k,j,i) ) / ll 3831 3832 diss(k,j,i) = dissipation * flag 3689 3833 3690 3834 ELSEIF ( rans_tke_e ) THEN 3691 3835 3692 dissipation (k,j)= diss(k,j,i)3836 dissipation = diss(k,j,i) 3693 3837 3694 3838 ENDIF … … 3710 3854 ) * ddzw(k) * drho_air(k) & 3711 3855 ) * flag * dsig_e & 3712  dissipation(k,j) * flag 3856  dissipation * flag 3857 ! 3858 ! Store dissipation if needed for calculating the sgs particle 3859 ! velocities 3860 IF ( .NOT. rans_tke_e .AND. ( use_sgs_for_particles .OR. & 3861 wang_kernel .OR. collision_turbulence ) ) THEN 3862 diss(k,j,i) = dissipation * flag 3863 ENDIF 3713 3864 3714 3865 ENDDO 3715 3866 ENDDO 3716 3717 !3718 ! Store dissipation if needed for calculating the sgs particle3719 ! velocities3720 IF ( .NOT. rans_tke_e .AND. ( use_sgs_for_particles .OR. &3721 wang_kernel .OR. collision_turbulence ) ) THEN3722 DO j = nys, nyn3723 DO k = nzb+1, nzt3724 diss(k,j,i) = dissipation(k,j) * MERGE( 1.0_wp, 0.0_wp, &3725 BTEST( wall_flags_0(k,j,i), 0 ) )3726 ENDDO3727 ENDDO3728 ENDIF3729 3730 3867 ENDDO 3731 3868 … … 4257 4394 ! Downward facing surfaces 4258 4395 !$OMP PARALLEL DO PRIVATE(i,j,k) 4396 !$ACC PARALLEL LOOP PRIVATE(i,j,k) & 4397 !$ACC PRESENT(bc_h(1), kh, km) 4259 4398 DO m = 1, bc_h(1)%ns 4260 4399 i = bc_h(1)%i(m) … … 4267 4406 ! Downward facing surfaces 4268 4407 !$OMP PARALLEL DO PRIVATE(i,j,k) 4408 !$ACC PARALLEL LOOP PRIVATE(i,j,k) & 4409 !$ACC PRESENT(bc_h(0), kh, km) 4269 4410 DO m = 1, bc_h(0)%ns 4270 4411 i = bc_h(0)%i(m) … … 4277 4418 ! Model top 4278 4419 !$OMP PARALLEL DO 4420 !$ACC PARALLEL LOOP COLLAPSE(2) & 4421 !$ACC PRESENT(kh, km) 4279 4422 DO i = nxlg, nxrg 4280 4423 DO j = nysg, nyng … … 4315 4458 SUBROUTINE tcm_diffusivities_default( var, var_reference ) 4316 4459 4460 USE arrays_3d, & 4461 ONLY: dd2zu 4462 4463 USE control_parameters, & 4464 ONLY: atmos_ocean_sign, use_single_reference_value, & 4465 wall_adjustment, wall_adjustment_factor 4466 4317 4467 USE statistics, & 4318 4468 ONLY : rmask, sums_l_l … … 4331 4481 REAL(wp) :: ll !< adjusted mixing length 4332 4482 REAL(wp) :: var_reference !< reference temperature 4483 #ifdef _OPENACC 4484 ! 4485 ! From mixing_length_les: 4486 REAL(wp) :: l_stable !< mixing length according to stratification 4487 REAL(wp) :: dvar_dz !< vertical gradient of var 4488 #endif 4333 4489 4334 4490 #if defined( __nopointer ) … … 4344 4500 ! 4345 4501 ! Initialization for calculation of the mixing length profile 4502 !$ACC KERNELS PRESENT(sums_l_l) 4346 4503 sums_l_l = 0.0_wp 4504 !$ACC END KERNELS 4347 4505 4348 4506 ! … … 4353 4511 IF ( les_dynamic .OR. les_mw ) THEN 4354 4512 !$OMP DO 4513 !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(sr) & 4514 !$ACC PRIVATE(flag, dvar_dz, l_stable, l, ll) & 4515 !$ACC PRESENT(wall_flags_0, var, dd2zu, e, l_wall, l_grid, rmask) & 4516 !$ACC PRESENT(kh, km, sums_l_l) 4355 4517 DO i = nxlg, nxrg 4356 4518 DO j = nysg, nyng … … 4361 4523 ! 4362 4524 ! Determine the mixing length for LES closure 4525 #ifdef _OPENACC 4526 ! 4527 ! Cannot call subroutine mixing_length_les because after adding all required 4528 ! OpenACC directives the execution crashed reliably when inside the called 4529 ! subroutine. I'm not sure why that is... 4530 dvar_dz = atmos_ocean_sign * ( var(k+1,j,i)  var(k1,j,i) ) * dd2zu(k) 4531 IF ( dvar_dz > 0.0_wp ) THEN 4532 IF ( use_single_reference_value ) THEN 4533 l_stable = 0.76_wp * SQRT( e(k,j,i) ) & 4534 / SQRT( g / var_reference * dvar_dz ) + 1E5_wp 4535 ELSE 4536 l_stable = 0.76_wp * SQRT( e(k,j,i) ) & 4537 / SQRT( g / var(k,j,i) * dvar_dz ) + 1E5_wp 4538 ENDIF 4539 ELSE 4540 l_stable = l_grid(k) 4541 ENDIF 4542 ! 4543 ! Adjustment of the mixing length 4544 IF ( wall_adjustment ) THEN 4545 l = MIN( wall_adjustment_factor * l_wall(k,j,i), l_grid(k), l_stable ) 4546 ll = MIN( wall_adjustment_factor * l_wall(k,j,i), l_grid(k) ) 4547 ELSE 4548 l = MIN( l_grid(k), l_stable ) 4549 ll = l_grid(k) 4550 ENDIF 4551 #else 4363 4552 CALL mixing_length_les( i, j, k, l, ll, var, var_reference ) 4553 #endif 4554 4364 4555 ! 4365 4556 ! Compute diffusion coefficients for momentum and heat … … 4429 4620 ENDIF 4430 4621 4622 !$ACC KERNELS PRESENT(sums_l_l) 4431 4623 sums_l_l(nzt+1,:,tn) = sums_l_l(nzt,:,tn) ! quasi boundarycondition for 4432 4624 ! data output 4625 !$ACC END KERNELS 4433 4626 !$OMP END PARALLEL 4434 4627
Note: See TracChangeset
for help on using the changeset viewer.