Changeset 3634 for palm/trunk/SOURCE/advec_ws.f90
- Timestamp:
- Dec 18, 2018 12:31:28 PM (6 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/advec_ws.f90 (modified) (80 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/advec_ws.f90
r3589 r3634 25 25 ! ----------------- 26 26 ! $Id$ 27 ! OpenACC port for SPEC 28 ! 29 ! 3589 2018-11-30 15:09:51Z suehring 27 30 ! Move the control parameter "salsa" from salsa_mod to control_parameters 28 31 ! (M. Kurppa) … … 1143 1146 !-- beginning of prognostic_equations. 1144 1147 IF ( ws_scheme_mom ) THEN 1148 !$ACC KERNELS PRESENT(sums_wsus_ws_l, sums_wsvs_ws_l) & 1149 !$ACC PRESENT(sums_us2_ws_l, sums_vs2_ws_l, sums_ws2_ws_l) 1145 1150 sums_wsus_ws_l = 0.0_wp 1146 1151 sums_wsvs_ws_l = 0.0_wp … … 1148 1153 sums_vs2_ws_l = 0.0_wp 1149 1154 sums_ws2_ws_l = 0.0_wp 1155 !$ACC END KERNELS 1150 1156 ENDIF 1151 1157 1152 1158 IF ( ws_scheme_sca ) THEN 1159 !$ACC KERNELS PRESENT(sums_wspts_ws_l) 1153 1160 sums_wspts_ws_l = 0.0_wp 1161 !$ACC END KERNELS 1154 1162 IF ( humidity ) sums_wsqs_ws_l = 0.0_wp 1155 1163 IF ( passive_scalar ) sums_wsss_ws_l = 0.0_wp … … 3300 3308 3301 3309 CHARACTER (LEN = *), INTENT(IN) :: sk_char !< string identifier, used for assign fluxes to the correct dimension in the analysis array 3310 INTEGER(iwp) :: sk_num !< integer identifier, used for assign fluxes to the correct dimension in the analysis array 3302 3311 3303 3312 INTEGER(iwp) :: i !< grid index along x-direction … … 3320 3329 REAL(wp) :: ibit1 !< flag indicating 3rd-order scheme along x-direction 3321 3330 REAL(wp) :: ibit2 !< flag indicating 5th-order scheme along x-direction 3331 #ifdef _OPENACC 3332 REAL(wp) :: ibit0_l !< flag indicating 1st-order scheme along x-direction 3333 REAL(wp) :: ibit1_l !< flag indicating 3rd-order scheme along x-direction 3334 REAL(wp) :: ibit2_l !< flag indicating 5th-order scheme along x-direction 3335 #endif 3322 3336 REAL(wp) :: ibit3 !< flag indicating 1st-order scheme along y-direction 3323 3337 REAL(wp) :: ibit4 !< flag indicating 3rd-order scheme along y-direction 3324 3338 REAL(wp) :: ibit5 !< flag indicating 5th-order scheme along y-direction 3339 #ifdef _OPENACC 3340 REAL(wp) :: ibit3_s !< flag indicating 1st-order scheme along y-direction 3341 REAL(wp) :: ibit4_s !< flag indicating 3rd-order scheme along y-direction 3342 REAL(wp) :: ibit5_s !< flag indicating 5th-order scheme along y-direction 3343 #endif 3325 3344 REAL(wp) :: ibit6 !< flag indicating 1st-order scheme along z-direction 3326 3345 REAL(wp) :: ibit7 !< flag indicating 3rd-order scheme along z-direction … … 3330 3349 REAL(wp) :: flux_d !< 6th-order flux at grid box bottom 3331 3350 REAL(wp) :: u_comp !< advection velocity along x-direction 3351 #ifdef _OPENACC 3352 REAL(wp) :: u_comp_l !< advection velocity along x-direction 3353 #endif 3332 3354 REAL(wp) :: v_comp !< advection velocity along y-direction 3355 #ifdef _OPENACC 3356 REAL(wp) :: v_comp_s !< advection velocity along y-direction 3357 #endif 3333 3358 3334 REAL(wp) , DIMENSION(nzb:nzt):: diss_n !< discretized artificial dissipation at northward-side of the grid box3335 REAL(wp) , DIMENSION(nzb:nzt):: diss_r !< discretized artificial dissipation at rightward-side of the grid box3336 REAL(wp) , DIMENSION(nzb:nzt):: diss_t !< discretized artificial dissipation at rightward-side of the grid box3337 REAL(wp) , DIMENSION(nzb:nzt):: flux_n !< discretized 6th-order flux at northward-side of the grid box3338 REAL(wp) , DIMENSION(nzb:nzt):: flux_r !< discretized 6th-order flux at rightward-side of the grid box3339 REAL(wp) , DIMENSION(nzb:nzt):: flux_t !< discretized 6th-order flux at rightward-side of the grid box3359 REAL(wp) :: diss_n !< discretized artificial dissipation at northward-side of the grid box 3360 REAL(wp) :: diss_r !< discretized artificial dissipation at rightward-side of the grid box 3361 REAL(wp) :: diss_t !< discretized artificial dissipation at rightward-side of the grid box 3362 REAL(wp) :: flux_n !< discretized 6th-order flux at northward-side of the grid box 3363 REAL(wp) :: flux_r !< discretized 6th-order flux at rightward-side of the grid box 3364 REAL(wp) :: flux_t !< discretized 6th-order flux at rightward-side of the grid box 3340 3365 3366 REAL(wp) :: diss_s !< discretized artificial dissipation term at southward-side of the grid box 3367 REAL(wp) :: flux_s !< discretized 6th-order flux at northward-side of the grid box 3368 #ifndef _OPENACC 3341 3369 REAL(wp), DIMENSION(nzb+1:nzt) :: swap_diss_y_local !< discretized artificial dissipation term at southward-side of the grid box 3342 3370 REAL(wp), DIMENSION(nzb+1:nzt) :: swap_flux_y_local !< discretized 6th-order flux at northward-side of the grid box 3371 #endif 3343 3372 3373 REAL(wp) :: diss_l !< discretized artificial dissipation term at leftward-side of the grid box 3374 REAL(wp) :: flux_l !< discretized 6th-order flux at leftward-side of the grid box 3375 #ifndef _OPENACC 3344 3376 REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) :: swap_diss_x_local !< discretized artificial dissipation term at leftward-side of the grid box 3345 3377 REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) :: swap_flux_x_local !< discretized 6th-order flux at leftward-side of the grid box 3378 #endif 3346 3379 3347 3380 SELECT CASE ( sk_char ) 3381 3382 CASE ( 'pt' ) 3383 sk_num = 1 3384 CASE ( 'sa' ) 3385 sk_num = 2 3386 CASE ( 'q' ) 3387 sk_num = 3 3388 CASE ( 'qc' ) 3389 sk_num = 4 3390 CASE ( 'qr' ) 3391 sk_num = 5 3392 CASE ( 'nc' ) 3393 sk_num = 6 3394 CASE ( 'nr' ) 3395 sk_num = 7 3396 CASE ( 's' ) 3397 sk_num = 8 3398 CASE ( 'aerosol_mass', 'aerosol_number', 'salsa_gas' ) 3399 sk_num = 9 3400 3401 END SELECT 3402 3403 #ifndef _OPENACC 3348 3404 ! 3349 3405 !-- Compute the fluxes for the whole left boundary of the processor domain. … … 3408 3464 3409 3465 ENDDO 3410 3466 #endif 3467 3468 !$ACC PARALLEL LOOP COLLAPSE(2) FIRSTPRIVATE(tn, sk_num) & 3469 !$ACC PRIVATE(i, j, k, k_mm, k_pp, k_ppp) & 3470 !$ACC PRIVATE(ibit0, ibit1, ibit2, ibit3, ibit4, ibit5) & 3471 !$ACC PRIVATE(ibit0_l, ibit1_l, ibit2_l) & 3472 !$ACC PRIVATE(ibit3_s, ibit4_s, ibit5_s) & 3473 !$ACC PRIVATE(ibit6, ibit7, ibit8) & 3474 !$ACC PRIVATE(flux_r, diss_r, flux_l, diss_l) & 3475 !$ACC PRIVATE(flux_n, diss_n, flux_s, diss_s) & 3476 !$ACC PRIVATE(flux_t, diss_t, flux_d, diss_d) & 3477 !$ACC PRIVATE(div, u_comp, u_comp_l, v_comp, v_comp_s) & 3478 !$ACC PRESENT(advc_flags_1) & 3479 !$ACC PRESENT(sk, u, v, w, u_stokes_zu, v_stokes_zu) & 3480 !$ACC PRESENT(drho_air, rho_air_zw, ddzw) & 3481 !$ACC PRESENT(tend) & 3482 !$ACC PRESENT(hom(nzb+1:nzb_max,1,1:3,0)) & 3483 !$ACC PRESENT(weight_substep(intermediate_timestep_count)) & 3484 !$ACC PRESENT(sums_wspts_ws_l, sums_wssas_ws_l) & 3485 !$ACC PRESENT(sums_wsqs_ws_l, sums_wsqcs_ws_l) & 3486 !$ACC PRESENT(sums_wsqrs_ws_l, sums_wsncs_ws_l) & 3487 !$ACC PRESENT(sums_wsnrs_ws_l, sums_wsss_ws_l) & 3488 !$ACC PRESENT(sums_salsa_ws_l) 3411 3489 DO i = nxl, nxr 3412 3490 3491 #ifndef _OPENACC 3413 3492 j = nys 3414 3493 DO k = nzb+1, nzb_max … … 3467 3546 3468 3547 ENDDO 3548 #endif 3469 3549 3470 3550 DO j = nys, nyn 3471 3551 3472 flux_t(0) = 0.0_wp3473 diss_t(0) = 0.0_wp3474 3552 flux_d = 0.0_wp 3475 3553 diss_d = 0.0_wp … … 3481 3559 ibit0 = REAL( IBITS(advc_flags_1(k,j,i),0,1), KIND = wp ) 3482 3560 3483 u_comp 3484 flux_r (k) = u_comp * (&3561 u_comp = u(k,j,i+1) - u_gtrans + u_stokes_zu(k) 3562 flux_r = u_comp * ( & 3485 3563 ( 37.0_wp * ibit2 * adv_sca_5 & 3486 3564 + 7.0_wp * ibit1 * adv_sca_3 & … … 3497 3575 ) 3498 3576 3499 diss_r (k) = -ABS( u_comp ) * (&3577 diss_r = -ABS( u_comp ) * ( & 3500 3578 ( 10.0_wp * ibit2 * adv_sca_5 & 3501 3579 + 3.0_wp * ibit1 * adv_sca_3 & … … 3512 3590 ) 3513 3591 3592 #ifdef _OPENACC 3593 ! 3594 !-- Recompute the left fluxes. 3595 ibit2_l = REAL( IBITS(advc_flags_1(k,j,i-1),2,1), KIND = wp ) 3596 ibit1_l = REAL( IBITS(advc_flags_1(k,j,i-1),1,1), KIND = wp ) 3597 ibit0_l = REAL( IBITS(advc_flags_1(k,j,i-1),0,1), KIND = wp ) 3598 3599 u_comp_l = u(k,j,i) - u_gtrans + u_stokes_zu(k) 3600 flux_l = u_comp_l * ( & 3601 ( 37.0_wp * ibit2_l * adv_sca_5 & 3602 + 7.0_wp * ibit1_l * adv_sca_3 & 3603 + ibit0_l * adv_sca_1 & 3604 ) * & 3605 ( sk(k,j,i) + sk(k,j,i-1) ) & 3606 - ( 8.0_wp * ibit2_l * adv_sca_5 & 3607 + ibit1_l * adv_sca_3 & 3608 ) * & 3609 ( sk(k,j,i+1) + sk(k,j,i-2) ) & 3610 + ( ibit2_l * adv_sca_5 & 3611 ) * & 3612 ( sk(k,j,i+2) + sk(k,j,i-3) ) & 3613 ) 3614 3615 diss_l = -ABS( u_comp_l ) * ( & 3616 ( 10.0_wp * ibit2_l * adv_sca_5 & 3617 + 3.0_wp * ibit1_l * adv_sca_3 & 3618 + ibit0_l * adv_sca_1 & 3619 ) * & 3620 ( sk(k,j,i) - sk(k,j,i-1) ) & 3621 - ( 5.0_wp * ibit2_l * adv_sca_5 & 3622 + ibit1_l * adv_sca_3 & 3623 ) * & 3624 ( sk(k,j,i+1) - sk(k,j,i-2) ) & 3625 + ( ibit2_l * adv_sca_5 & 3626 ) * & 3627 ( sk(k,j,i+2) - sk(k,j,i-3) ) & 3628 ) 3629 #else 3630 flux_l = swap_flux_x_local(k,j) 3631 diss_l = swap_diss_x_local(k,j) 3632 #endif 3633 3514 3634 ibit5 = REAL( IBITS(advc_flags_1(k,j,i),5,1), KIND = wp ) 3515 3635 ibit4 = REAL( IBITS(advc_flags_1(k,j,i),4,1), KIND = wp ) 3516 3636 ibit3 = REAL( IBITS(advc_flags_1(k,j,i),3,1), KIND = wp ) 3517 3637 3518 v_comp 3519 flux_n (k) = v_comp * (&3638 v_comp = v(k,j+1,i) - v_gtrans + v_stokes_zu(k) 3639 flux_n = v_comp * ( & 3520 3640 ( 37.0_wp * ibit5 * adv_sca_5 & 3521 3641 + 7.0_wp * ibit4 * adv_sca_3 & … … 3532 3652 ) 3533 3653 3534 diss_n (k) = -ABS( v_comp ) * (&3654 diss_n = -ABS( v_comp ) * ( & 3535 3655 ( 10.0_wp * ibit5 * adv_sca_5 & 3536 3656 + 3.0_wp * ibit4 * adv_sca_3 & … … 3546 3666 ( sk(k,j+3,i) - sk(k,j-2,i) ) & 3547 3667 ) 3668 3669 #ifdef _OPENACC 3670 ! 3671 !-- Recompute the south fluxes. 3672 ibit5_s = REAL( IBITS(advc_flags_1(k,j-1,i),5,1), KIND = wp ) 3673 ibit4_s = REAL( IBITS(advc_flags_1(k,j-1,i),4,1), KIND = wp ) 3674 ibit3_s = REAL( IBITS(advc_flags_1(k,j-1,i),3,1), KIND = wp ) 3675 3676 v_comp_s = v(k,j,i) - v_gtrans + v_stokes_zu(k) 3677 flux_s = v_comp_s * ( & 3678 ( 37.0_wp * ibit5_s * adv_sca_5 & 3679 + 7.0_wp * ibit4_s * adv_sca_3 & 3680 + ibit3_s * adv_sca_1 & 3681 ) * & 3682 ( sk(k,j,i) + sk(k,j-1,i) ) & 3683 - ( 8.0_wp * ibit5_s * adv_sca_5 & 3684 + ibit4_s * adv_sca_3 & 3685 ) * & 3686 ( sk(k,j+1,i) + sk(k,j-2,i) ) & 3687 + ( ibit5_s * adv_sca_5 & 3688 ) * & 3689 ( sk(k,j+2,i) + sk(k,j-3,i) ) & 3690 ) 3691 3692 diss_s = -ABS( v_comp_s ) * ( & 3693 ( 10.0_wp * ibit5_s * adv_sca_5 & 3694 + 3.0_wp * ibit4_s * adv_sca_3 & 3695 + ibit3_s * adv_sca_1 & 3696 ) * & 3697 ( sk(k,j,i) - sk(k,j-1,i) ) & 3698 - ( 5.0_wp * ibit5_s * adv_sca_5 & 3699 + ibit4_s * adv_sca_3 & 3700 ) * & 3701 ( sk(k,j+1,i) - sk(k,j-2,i) ) & 3702 + ( ibit5_s * adv_sca_5 & 3703 ) * & 3704 ( sk(k,j+2,i) - sk(k,j-3,i) ) & 3705 ) 3706 #else 3707 flux_s = swap_flux_y_local(k) 3708 diss_s = swap_diss_y_local(k) 3709 #endif 3710 3548 3711 ! 3549 3712 !-- k index has to be modified near bottom and top, else array … … 3558 3721 3559 3722 3560 flux_t (k) = w(k,j,i) * rho_air_zw(k) * (&3723 flux_t = w(k,j,i) * rho_air_zw(k) * ( & 3561 3724 ( 37.0_wp * ibit8 * adv_sca_5 & 3562 3725 + 7.0_wp * ibit7 * adv_sca_3 & … … 3572 3735 ) 3573 3736 3574 diss_t (k) = -ABS( w(k,j,i) ) * rho_air_zw(k) * (&3737 diss_t = -ABS( w(k,j,i) ) * rho_air_zw(k) * ( & 3575 3738 ( 10.0_wp * ibit8 * adv_sca_5 & 3576 3739 + 3.0_wp * ibit7 * adv_sca_3 & … … 3616 3779 3617 3780 tend(k,j,i) = tend(k,j,i) - ( & 3618 ( flux_r(k) + diss_r(k) - swap_flux_x_local(k,j) - & 3619 swap_diss_x_local(k,j) ) * ddx & 3620 + ( flux_n(k) + diss_n(k) - swap_flux_y_local(k) - & 3621 swap_diss_y_local(k) ) * ddy & 3622 + ( ( flux_t(k) + diss_t(k) ) - & 3623 ( flux_d + diss_d ) & 3624 ) * drho_air(k) * ddzw(k) & 3781 ( ( flux_r + diss_r ) & 3782 - ( flux_l + diss_l ) ) * ddx & 3783 + ( ( flux_n + diss_n ) & 3784 - ( flux_s + diss_s ) ) * ddy & 3785 + ( ( flux_t + diss_t ) & 3786 - ( flux_d + diss_d ) ) * drho_air(k) * ddzw(k) & 3625 3787 ) + sk(k,j,i) * div 3626 3788 3627 swap_flux_y_local(k) = flux_n(k) 3628 swap_diss_y_local(k) = diss_n(k) 3629 swap_flux_x_local(k,j) = flux_r(k) 3630 swap_diss_x_local(k,j) = diss_r(k) 3631 flux_d = flux_t(k) 3632 diss_d = diss_t(k) 3789 #ifndef _OPENACC 3790 swap_flux_y_local(k) = flux_n 3791 swap_diss_y_local(k) = diss_n 3792 swap_flux_x_local(k,j) = flux_r 3793 swap_diss_x_local(k,j) = diss_r 3794 #endif 3795 flux_d = flux_t 3796 diss_d = diss_t 3797 3798 ! 3799 !-- Evaluation of statistics. 3800 SELECT CASE ( sk_num ) 3801 3802 CASE ( 1 ) 3803 !$ACC ATOMIC 3804 sums_wspts_ws_l(k,tn) = sums_wspts_ws_l(k,tn) & 3805 + ( flux_t & 3806 / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) ) & 3807 * ( w(k,j,i) - hom(k,1,3,0) ) & 3808 + diss_t & 3809 / ( ABS(w(k,j,i)) + 1.0E-20_wp ) & 3810 * ABS(w(k,j,i) - hom(k,1,3,0) ) & 3811 ) * weight_substep(intermediate_timestep_count) 3812 CASE ( 2 ) 3813 !$ACC ATOMIC 3814 sums_wssas_ws_l(k,tn) = sums_wssas_ws_l(k,tn) & 3815 + ( flux_t & 3816 / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) ) & 3817 * ( w(k,j,i) - hom(k,1,3,0) ) & 3818 + diss_t & 3819 / ( ABS(w(k,j,i)) + 1.0E-20_wp ) & 3820 * ABS(w(k,j,i) - hom(k,1,3,0) ) & 3821 ) * weight_substep(intermediate_timestep_count) 3822 CASE ( 3 ) 3823 !$ACC ATOMIC 3824 sums_wsqs_ws_l(k,tn) = sums_wsqs_ws_l(k,tn) & 3825 + ( flux_t & 3826 / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) ) & 3827 * ( w(k,j,i) - hom(k,1,3,0) ) & 3828 + diss_t & 3829 / ( ABS(w(k,j,i)) + 1.0E-20_wp ) & 3830 * ABS(w(k,j,i) - hom(k,1,3,0) ) & 3831 ) * weight_substep(intermediate_timestep_count) 3832 CASE ( 4 ) 3833 !$ACC ATOMIC 3834 sums_wsqcs_ws_l(k,tn) = sums_wsqcs_ws_l(k,tn) & 3835 + ( flux_t & 3836 / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) ) & 3837 * ( w(k,j,i) - hom(k,1,3,0) ) & 3838 + diss_t & 3839 / ( ABS(w(k,j,i)) + 1.0E-20_wp ) & 3840 * ABS(w(k,j,i) - hom(k,1,3,0) ) & 3841 ) * weight_substep(intermediate_timestep_count) 3842 CASE ( 5 ) 3843 !$ACC ATOMIC 3844 sums_wsqrs_ws_l(k,tn) = sums_wsqrs_ws_l(k,tn) & 3845 + ( flux_t & 3846 / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) ) & 3847 * ( w(k,j,i) - hom(k,1,3,0) ) & 3848 + diss_t & 3849 / ( ABS(w(k,j,i)) + 1.0E-20_wp ) & 3850 * ABS(w(k,j,i) - hom(k,1,3,0) ) & 3851 ) * weight_substep(intermediate_timestep_count) 3852 CASE ( 6 ) 3853 !$ACC ATOMIC 3854 sums_wsncs_ws_l(k,tn) = sums_wsncs_ws_l(k,tn) & 3855 + ( flux_t & 3856 / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) ) & 3857 * ( w(k,j,i) - hom(k,1,3,0) ) & 3858 + diss_t & 3859 / ( ABS(w(k,j,i)) + 1.0E-20_wp ) & 3860 * ABS(w(k,j,i) - hom(k,1,3,0) ) & 3861 ) * weight_substep(intermediate_timestep_count) 3862 CASE ( 7 ) 3863 !$ACC ATOMIC 3864 sums_wsnrs_ws_l(k,tn) = sums_wsnrs_ws_l(k,tn) & 3865 + ( flux_t & 3866 / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) ) & 3867 * ( w(k,j,i) - hom(k,1,3,0) ) & 3868 + diss_t & 3869 / ( ABS(w(k,j,i)) + 1.0E-20_wp ) & 3870 * ABS(w(k,j,i) - hom(k,1,3,0) ) & 3871 ) * weight_substep(intermediate_timestep_count) 3872 CASE ( 8 ) 3873 !$ACC ATOMIC 3874 sums_wsss_ws_l(k,tn) = sums_wsss_ws_l(k,tn) & 3875 + ( flux_t & 3876 / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) ) & 3877 * ( w(k,j,i) - hom(k,1,3,0) ) & 3878 + diss_t & 3879 / ( ABS(w(k,j,i)) + 1.0E-20_wp ) & 3880 * ABS(w(k,j,i) - hom(k,1,3,0) ) & 3881 ) * weight_substep(intermediate_timestep_count) 3882 CASE ( 9 ) 3883 !$ACC ATOMIC 3884 sums_salsa_ws_l(k,tn) = sums_salsa_ws_l(k,tn) & 3885 + ( flux_t & 3886 / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) ) & 3887 * ( w(k,j,i) - hom(k,1,3,0) ) & 3888 + diss_t & 3889 / ( ABS(w(k,j,i)) + 1.0E-20_wp ) & 3890 * ABS(w(k,j,i) - hom(k,1,3,0) ) & 3891 ) * weight_substep(intermediate_timestep_count) 3892 3893 END SELECT 3633 3894 3634 3895 ENDDO … … 3636 3897 DO k = nzb_max+1, nzt 3637 3898 3638 u_comp 3639 flux_r (k) = u_comp * (&3899 u_comp = u(k,j,i+1) - u_gtrans + u_stokes_zu(k) 3900 flux_r = u_comp * ( & 3640 3901 37.0_wp * ( sk(k,j,i+1) + sk(k,j,i) ) & 3641 3902 - 8.0_wp * ( sk(k,j,i+2) + sk(k,j,i-1) ) & 3642 3903 + ( sk(k,j,i+3) + sk(k,j,i-2) ) ) * adv_sca_5 3643 diss_r (k) = -ABS( u_comp ) * (&3904 diss_r = -ABS( u_comp ) * ( & 3644 3905 10.0_wp * ( sk(k,j,i+1) - sk(k,j,i) ) & 3645 3906 - 5.0_wp * ( sk(k,j,i+2) - sk(k,j,i-1) ) & 3646 3907 + ( sk(k,j,i+3) - sk(k,j,i-2) ) ) * adv_sca_5 3647 3908 3648 v_comp = v(k,j+1,i) - v_gtrans + v_stokes_zu(k) 3649 flux_n(k) = v_comp * ( & 3909 #ifdef _OPENACC 3910 ! 3911 !-- Recompute the left fluxes. 3912 u_comp_l = u(k,j,i) - u_gtrans + u_stokes_zu(k) 3913 flux_l = u_comp_l * ( & 3914 37.0_wp * ( sk(k,j,i) + sk(k,j,i-1) ) & 3915 - 8.0_wp * ( sk(k,j,i+1) + sk(k,j,i-2) ) & 3916 + ( sk(k,j,i+2) + sk(k,j,i-3) ) & 3917 ) * adv_sca_5 3918 3919 diss_l = -ABS( u_comp_l ) * ( & 3920 10.0_wp * ( sk(k,j,i) - sk(k,j,i-1) ) & 3921 - 5.0_wp * ( sk(k,j,i+1) - sk(k,j,i-2) ) & 3922 + ( sk(k,j,i+2) - sk(k,j,i-3) ) & 3923 ) * adv_sca_5 3924 #else 3925 flux_l = swap_flux_x_local(k,j) 3926 diss_l = swap_diss_x_local(k,j) 3927 #endif 3928 3929 v_comp = v(k,j+1,i) - v_gtrans + v_stokes_zu(k) 3930 flux_n = v_comp * ( & 3650 3931 37.0_wp * ( sk(k,j+1,i) + sk(k,j,i) ) & 3651 3932 - 8.0_wp * ( sk(k,j+2,i) + sk(k,j-1,i) ) & 3652 3933 + ( sk(k,j+3,i) + sk(k,j-2,i) ) ) * adv_sca_5 3653 diss_n (k) = -ABS( v_comp ) * (&3934 diss_n = -ABS( v_comp ) * ( & 3654 3935 10.0_wp * ( sk(k,j+1,i) - sk(k,j,i) ) & 3655 3936 - 5.0_wp * ( sk(k,j+2,i) - sk(k,j-1,i) ) & 3656 3937 + ( sk(k,j+3,i) - sk(k,j-2,i) ) ) * adv_sca_5 3938 3939 #ifdef _OPENACC 3940 ! 3941 !-- Recompute the south fluxes. 3942 v_comp_s = v(k,j,i) - v_gtrans + v_stokes_zu(k) 3943 flux_s = v_comp_s * ( & 3944 37.0_wp * ( sk(k,j,i) + sk(k,j-1,i) ) & 3945 - 8.0_wp * ( sk(k,j+1,i) + sk(k,j-2,i) ) & 3946 + ( sk(k,j+2,i) + sk(k,j-3,i) ) & 3947 ) * adv_sca_5 3948 diss_s = -ABS( v_comp_s ) * ( & 3949 10.0_wp * ( sk(k,j,i) - sk(k,j-1,i) ) & 3950 - 5.0_wp * ( sk(k,j+1,i) - sk(k,j-2,i) ) & 3951 + sk(k,j+2,i) - sk(k,j-3,i) & 3952 ) * adv_sca_5 3953 #else 3954 flux_s = swap_flux_y_local(k) 3955 diss_s = swap_diss_y_local(k) 3956 #endif 3957 3657 3958 ! 3658 3959 !-- k index has to be modified near bottom and top, else array … … 3667 3968 3668 3969 3669 flux_t (k)= w(k,j,i) * rho_air_zw(k) * ( &3970 flux_t = w(k,j,i) * rho_air_zw(k) * ( & 3670 3971 ( 37.0_wp * ibit8 * adv_sca_5 & 3671 3972 + 7.0_wp * ibit7 * adv_sca_3 & … … 3681 3982 ) 3682 3983 3683 diss_t (k)= -ABS( w(k,j,i) ) * rho_air_zw(k) * ( &3984 diss_t = -ABS( w(k,j,i) ) * rho_air_zw(k) * ( & 3684 3985 ( 10.0_wp * ibit8 * adv_sca_5 & 3685 3986 + 3.0_wp * ibit7 * adv_sca_3 & … … 3706 4007 3707 4008 tend(k,j,i) = tend(k,j,i) - ( & 3708 ( flux_r(k) + diss_r(k) - swap_flux_x_local(k,j) - & 3709 swap_diss_x_local(k,j) ) * ddx & 3710 + ( flux_n(k) + diss_n(k) - swap_flux_y_local(k) - & 3711 swap_diss_y_local(k) ) * ddy & 3712 + ( ( flux_t(k) + diss_t(k) ) - & 3713 ( flux_d + diss_d ) & 3714 ) * drho_air(k) * ddzw(k) & 4009 ( ( flux_r + diss_r ) & 4010 - ( flux_l + diss_l ) ) * ddx & 4011 + ( ( flux_n + diss_n ) & 4012 - ( flux_s + diss_s ) ) * ddy & 4013 + ( ( flux_t + diss_t ) & 4014 - ( flux_d + diss_d ) ) * drho_air(k) * ddzw(k) & 3715 4015 ) + sk(k,j,i) * div 3716 4016 3717 swap_flux_y_local(k) = flux_n(k) 3718 swap_diss_y_local(k) = diss_n(k) 3719 swap_flux_x_local(k,j) = flux_r(k) 3720 swap_diss_x_local(k,j) = diss_r(k) 3721 flux_d = flux_t(k) 3722 diss_d = diss_t(k) 3723 3724 ENDDO 3725 ! 3726 !-- Evaluation of statistics. 3727 SELECT CASE ( sk_char ) 3728 3729 CASE ( 'pt' ) 3730 DO k = nzb, nzt 4017 #ifndef _OPENACC 4018 swap_flux_y_local(k) = flux_n 4019 swap_diss_y_local(k) = diss_n 4020 swap_flux_x_local(k,j) = flux_r 4021 swap_diss_x_local(k,j) = diss_r 4022 #endif 4023 flux_d = flux_t 4024 diss_d = diss_t 4025 4026 ! 4027 !-- Evaluation of statistics. 4028 SELECT CASE ( sk_num ) 4029 4030 CASE ( 1 ) 4031 !$ACC ATOMIC 3731 4032 sums_wspts_ws_l(k,tn) = sums_wspts_ws_l(k,tn) & 3732 + ( flux_t (k)&4033 + ( flux_t & 3733 4034 / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) ) & 3734 4035 * ( w(k,j,i) - hom(k,1,3,0) ) & 3735 + diss_t (k)&4036 + diss_t & 3736 4037 / ( ABS(w(k,j,i)) + 1.0E-20_wp ) & 3737 4038 * ABS(w(k,j,i) - hom(k,1,3,0) ) & 3738 4039 ) * weight_substep(intermediate_timestep_count) 3739 ENDDO 3740 CASE ( 'sa' ) 3741 DO k = nzb, nzt 4040 CASE ( 2 ) 4041 !$ACC ATOMIC 3742 4042 sums_wssas_ws_l(k,tn) = sums_wssas_ws_l(k,tn) & 3743 + ( flux_t (k)&4043 + ( flux_t & 3744 4044 / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) ) & 3745 4045 * ( w(k,j,i) - hom(k,1,3,0) ) & 3746 + diss_t (k)&4046 + diss_t & 3747 4047 / ( ABS(w(k,j,i)) + 1.0E-20_wp ) & 3748 4048 * ABS(w(k,j,i) - hom(k,1,3,0) ) & 3749 4049 ) * weight_substep(intermediate_timestep_count) 3750 ENDDO 3751 CASE ( 'q' ) 3752 DO k = nzb, nzt 4050 CASE ( 3 ) 4051 !$ACC ATOMIC 3753 4052 sums_wsqs_ws_l(k,tn) = sums_wsqs_ws_l(k,tn) & 3754 + ( flux_t (k)&4053 + ( flux_t & 3755 4054 / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) ) & 3756 4055 * ( w(k,j,i) - hom(k,1,3,0) ) & 3757 + diss_t (k)&4056 + diss_t & 3758 4057 / ( ABS(w(k,j,i)) + 1.0E-20_wp ) & 3759 4058 * ABS(w(k,j,i) - hom(k,1,3,0) ) & 3760 4059 ) * weight_substep(intermediate_timestep_count) 3761 ENDDO 3762 CASE ( 'qc' ) 3763 DO k = nzb, nzt 4060 CASE ( 4 ) 4061 !$ACC ATOMIC 3764 4062 sums_wsqcs_ws_l(k,tn) = sums_wsqcs_ws_l(k,tn) & 3765 + ( flux_t (k)&4063 + ( flux_t & 3766 4064 / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) ) & 3767 4065 * ( w(k,j,i) - hom(k,1,3,0) ) & 3768 + diss_t (k)&4066 + diss_t & 3769 4067 / ( ABS(w(k,j,i)) + 1.0E-20_wp ) & 3770 4068 * ABS(w(k,j,i) - hom(k,1,3,0) ) & 3771 4069 ) * weight_substep(intermediate_timestep_count) 3772 ENDDO 3773 CASE ( 'qr' ) 3774 DO k = nzb, nzt 4070 CASE ( 5 ) 4071 !$ACC ATOMIC 3775 4072 sums_wsqrs_ws_l(k,tn) = sums_wsqrs_ws_l(k,tn) & 3776 + ( flux_t (k)&4073 + ( flux_t & 3777 4074 / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) ) & 3778 4075 * ( w(k,j,i) - hom(k,1,3,0) ) & 3779 + diss_t (k)&4076 + diss_t & 3780 4077 / ( ABS(w(k,j,i)) + 1.0E-20_wp ) & 3781 4078 * ABS(w(k,j,i) - hom(k,1,3,0) ) & 3782 4079 ) * weight_substep(intermediate_timestep_count) 3783 ENDDO 3784 CASE ( 'nc' ) 3785 DO k = nzb, nzt 4080 CASE ( 6 ) 4081 !$ACC ATOMIC 3786 4082 sums_wsncs_ws_l(k,tn) = sums_wsncs_ws_l(k,tn) & 3787 + ( flux_t (k)&4083 + ( flux_t & 3788 4084 / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) ) & 3789 4085 * ( w(k,j,i) - hom(k,1,3,0) ) & 3790 + diss_t (k)&4086 + diss_t & 3791 4087 / ( ABS(w(k,j,i)) + 1.0E-20_wp ) & 3792 4088 * ABS(w(k,j,i) - hom(k,1,3,0) ) & 3793 4089 ) * weight_substep(intermediate_timestep_count) 3794 ENDDO 3795 CASE ( 'nr' ) 3796 DO k = nzb, nzt 4090 CASE ( 7 ) 4091 !$ACC ATOMIC 3797 4092 sums_wsnrs_ws_l(k,tn) = sums_wsnrs_ws_l(k,tn) & 3798 + ( flux_t (k)&4093 + ( flux_t & 3799 4094 / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) ) & 3800 4095 * ( w(k,j,i) - hom(k,1,3,0) ) & 3801 + diss_t (k)&4096 + diss_t & 3802 4097 / ( ABS(w(k,j,i)) + 1.0E-20_wp ) & 3803 4098 * ABS(w(k,j,i) - hom(k,1,3,0) ) & 3804 4099 ) * weight_substep(intermediate_timestep_count) 3805 ENDDO 3806 CASE ( 's' ) 3807 DO k = nzb, nzt 4100 CASE ( 8 ) 4101 !$ACC ATOMIC 3808 4102 sums_wsss_ws_l(k,tn) = sums_wsss_ws_l(k,tn) & 3809 + ( flux_t (k)&4103 + ( flux_t & 3810 4104 / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) ) & 3811 4105 * ( w(k,j,i) - hom(k,1,3,0) ) & 3812 + diss_t (k)&4106 + diss_t & 3813 4107 / ( ABS(w(k,j,i)) + 1.0E-20_wp ) & 3814 4108 * ABS(w(k,j,i) - hom(k,1,3,0) ) & 3815 4109 ) * weight_substep(intermediate_timestep_count) 3816 ENDDO 3817 3818 CASE ( 'aerosol_mass', 'aerosol_number', 'salsa_gas' ) 3819 DO k = nzb, nzt 4110 CASE ( 9 ) 4111 !$ACC ATOMIC 3820 4112 sums_salsa_ws_l(k,tn) = sums_salsa_ws_l(k,tn) & 3821 + ( flux_t (k)&4113 + ( flux_t & 3822 4114 / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) ) & 3823 4115 * ( w(k,j,i) - hom(k,1,3,0) ) & 3824 + diss_t (k)&4116 + diss_t & 3825 4117 / ( ABS(w(k,j,i)) + 1.0E-20_wp ) & 3826 4118 * ABS(w(k,j,i) - hom(k,1,3,0) ) & 3827 4119 ) * weight_substep(intermediate_timestep_count) 3828 ENDDO 3829 3830 3831 END SELECT4120 4121 END SELECT 4122 4123 ENDDO 3832 4124 3833 4125 ENDDO … … 3874 4166 REAL(wp) :: ibit10 !< flag indicating 3rd-order scheme along x-direction 3875 4167 REAL(wp) :: ibit11 !< flag indicating 5th-order scheme along x-direction 4168 #ifdef _OPENACC 4169 REAL(wp) :: ibit9_l !< flag indicating 1st-order scheme along x-direction 4170 REAL(wp) :: ibit10_l !< flag indicating 3rd-order scheme along x-direction 4171 REAL(wp) :: ibit11_l !< flag indicating 5th-order scheme along x-direction 4172 #endif 3876 4173 REAL(wp) :: ibit12 !< flag indicating 1st-order scheme along y-direction 3877 4174 REAL(wp) :: ibit13 !< flag indicating 3rd-order scheme along y-direction 3878 4175 REAL(wp) :: ibit14 !< flag indicating 5th-order scheme along y-direction 4176 #ifdef _OPENACC 4177 REAL(wp) :: ibit12_s !< flag indicating 1st-order scheme along y-direction 4178 REAL(wp) :: ibit13_s !< flag indicating 3rd-order scheme along y-direction 4179 REAL(wp) :: ibit14_s !< flag indicating 5th-order scheme along y-direction 4180 #endif 3879 4181 REAL(wp) :: ibit15 !< flag indicating 1st-order scheme along z-direction 3880 4182 REAL(wp) :: ibit16 !< flag indicating 3rd-order scheme along z-direction … … 3886 4188 REAL(wp) :: gv !< Galilei-transformation velocity along y 3887 4189 REAL(wp) :: v_comp !< advection velocity along y 4190 #ifdef _OPENACC 4191 REAL(wp) :: v_comp_s !< advection velocity along y 4192 #endif 3888 4193 REAL(wp) :: w_comp !< advection velocity along z 3889 4194 4195 REAL(wp) :: diss_s !< discretized artificial dissipation at southward-side of the grid box 4196 REAL(wp) :: flux_s !< discretized 6th-order flux at southward-side of the grid box 4197 #ifndef _OPENACC 3890 4198 REAL(wp), DIMENSION(nzb+1:nzt) :: swap_diss_y_local_u !< discretized artificial dissipation at southward-side of the grid box 3891 4199 REAL(wp), DIMENSION(nzb+1:nzt) :: swap_flux_y_local_u !< discretized 6th-order flux at southward-side of the grid box 4200 #endif 3892 4201 4202 REAL(wp) :: diss_l !< discretized artificial dissipation at leftward-side of the grid box 4203 REAL(wp) :: flux_l !< discretized 6th-order flux at leftward-side of the grid box 4204 #ifndef _OPENACC 3893 4205 REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) :: swap_diss_x_local_u !< discretized artificial dissipation at leftward-side of the grid box 3894 4206 REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) :: swap_flux_x_local_u !< discretized 6th-order flux at leftward-side of the grid box 4207 #endif 3895 4208 3896 REAL(wp), DIMENSION(nzb:nzt) :: diss_n !< discretized artificial dissipation at northward-side of the grid box 3897 REAL(wp), DIMENSION(nzb:nzt) :: diss_r !< discretized artificial dissipation at leftward-side of the grid box 3898 REAL(wp), DIMENSION(nzb:nzt) :: diss_t !< discretized artificial dissipation at top of the grid box 3899 REAL(wp), DIMENSION(nzb:nzt) :: flux_n !< discretized 6th-order flux at northward-side of the grid box 3900 REAL(wp), DIMENSION(nzb:nzt) :: flux_r !< discretized 6th-order flux at rightward-side of the grid box 3901 REAL(wp), DIMENSION(nzb:nzt) :: flux_t !< discretized 6th-order flux at top of the grid box 3902 REAL(wp), DIMENSION(nzb:nzt) :: u_comp !< advection velocity along x 4209 REAL(wp) :: diss_n !< discretized artificial dissipation at northward-side of the grid box 4210 REAL(wp) :: diss_r !< discretized artificial dissipation at leftward-side of the grid box 4211 REAL(wp) :: diss_t !< discretized artificial dissipation at top of the grid box 4212 REAL(wp) :: flux_n !< discretized 6th-order flux at northward-side of the grid box 4213 REAL(wp) :: flux_r !< discretized 6th-order flux at rightward-side of the grid box 4214 REAL(wp) :: flux_t !< discretized 6th-order flux at top of the grid box 4215 REAL(wp) :: u_comp !< advection velocity along x 4216 #ifdef _OPENACC 4217 REAL(wp) :: u_comp_l !< 4218 #endif 3903 4219 3904 4220 gu = 2.0_wp * u_gtrans 3905 4221 gv = 2.0_wp * v_gtrans 3906 4222 4223 #ifndef _OPENACC 3907 4224 ! 3908 4225 !-- Compute the fluxes for the whole left boundary of the processor domain. … … 3915 4232 ibit9 = REAL( IBITS(advc_flags_1(k,j,i-1),9,1), KIND = wp ) 3916 4233 3917 u_comp (k)= u(k,j,i) + u(k,j,i-1) - gu3918 swap_flux_x_local_u(k,j) = u_comp (k) * (&4234 u_comp = u(k,j,i) + u(k,j,i-1) - gu 4235 swap_flux_x_local_u(k,j) = u_comp * ( & 3919 4236 ( 37.0_wp * ibit11 * adv_mom_5 & 3920 4237 + 7.0_wp * ibit10 * adv_mom_3 & … … 3931 4248 ) 3932 4249 3933 swap_diss_x_local_u(k,j) = - ABS( u_comp (k) ) * (&4250 swap_diss_x_local_u(k,j) = - ABS( u_comp ) * ( & 3934 4251 ( 10.0_wp * ibit11 * adv_mom_5 & 3935 4252 + 3.0_wp * ibit10 * adv_mom_3 & … … 3950 4267 DO k = nzb_max+1, nzt 3951 4268 3952 u_comp (k)= u(k,j,i) + u(k,j,i-1) - gu3953 swap_flux_x_local_u(k,j) = u_comp (k) * (&4269 u_comp = u(k,j,i) + u(k,j,i-1) - gu 4270 swap_flux_x_local_u(k,j) = u_comp * ( & 3954 4271 37.0_wp * ( u(k,j,i) + u(k,j,i-1) ) & 3955 4272 - 8.0_wp * ( u(k,j,i+1) + u(k,j,i-2) ) & 3956 4273 + ( u(k,j,i+2) + u(k,j,i-3) ) ) * adv_mom_5 3957 swap_diss_x_local_u(k,j) = - ABS(u_comp (k)) * (&4274 swap_diss_x_local_u(k,j) = - ABS(u_comp) * ( & 3958 4275 10.0_wp * ( u(k,j,i) - u(k,j,i-1) ) & 3959 4276 - 5.0_wp * ( u(k,j,i+1) - u(k,j,i-2) ) & … … 3962 4279 ENDDO 3963 4280 ENDDO 3964 4281 #endif 4282 4283 !$ACC PARALLEL LOOP COLLAPSE(2) FIRSTPRIVATE(tn, gu, gv) & 4284 !$ACC PRIVATE(i, j, k, k_mm, k_pp, k_ppp) & 4285 !$ACC PRIVATE(ibit9, ibit10, ibit11, ibit12, ibit13, ibit14) & 4286 !$ACC PRIVATE(ibit9_l, ibit10_l, ibit11_l) & 4287 !$ACC PRIVATE(ibit12_s, ibit13_s, ibit14_s) & 4288 !$ACC PRIVATE(ibit15, ibit16, ibit17) & 4289 !$ACC PRIVATE(flux_r, diss_r, flux_l, diss_l) & 4290 !$ACC PRIVATE(flux_n, diss_n, flux_s, diss_s) & 4291 !$ACC PRIVATE(flux_t, diss_t, flux_d, diss_d) & 4292 !$ACC PRIVATE(div, u_comp, u_comp_l, v_comp, v_comp_s, w_comp) & 4293 !$ACC PRESENT(advc_flags_1) & 4294 !$ACC PRESENT(u, v, w) & 4295 !$ACC PRESENT(drho_air, rho_air_zw, ddzw) & 4296 !$ACC PRESENT(tend) & 4297 !$ACC PRESENT(hom(nzb+1:nzb_max,1,1:3,0)) & 4298 !$ACC PRESENT(weight_substep(intermediate_timestep_count)) & 4299 !$ACC PRESENT(sums_us2_ws_l, sums_wsus_ws_l) 3965 4300 DO i = nxlu, nxr 4301 #ifndef _OPENACC 3966 4302 ! 3967 4303 !-- The following loop computes the fluxes for the south boundary points … … 4019 4355 4020 4356 ENDDO 4357 #endif 4358 4021 4359 ! 4022 4360 !-- Computation of interior fluxes and tendency terms 4023 4361 DO j = nys, nyn 4024 4362 4025 flux_t(0) = 0.0_wp4026 diss_t(0) = 0.0_wp4027 4363 flux_d = 0.0_wp 4028 4364 diss_d = 0.0_wp … … 4034 4370 ibit9 = REAL( IBITS(advc_flags_1(k,j,i),9,1), KIND = wp ) 4035 4371 4036 u_comp (k)= u(k,j,i+1) + u(k,j,i)4037 flux_r (k) = ( u_comp(k) - gu ) * (&4372 u_comp = u(k,j,i+1) + u(k,j,i) 4373 flux_r = ( u_comp - gu ) * ( & 4038 4374 ( 37.0_wp * ibit11 * adv_mom_5 & 4039 4375 + 7.0_wp * ibit10 * adv_mom_3 & … … 4050 4386 ) 4051 4387 4052 diss_r (k) = - ABS( u_comp(k) - gu ) * (&4388 diss_r = - ABS( u_comp - gu ) * ( & 4053 4389 ( 10.0_wp * ibit11 * adv_mom_5 & 4054 4390 + 3.0_wp * ibit10 * adv_mom_3 & … … 4065 4401 ) 4066 4402 4403 #ifdef _OPENACC 4404 ! 4405 !-- Recompute the left fluxes. 4406 ibit11_l = REAL( IBITS(advc_flags_1(k,j,i-1),11,1), KIND = wp ) 4407 ibit10_l = REAL( IBITS(advc_flags_1(k,j,i-1),10,1), KIND = wp ) 4408 ibit9_l = REAL( IBITS(advc_flags_1(k,j,i-1),9,1), KIND = wp ) 4409 4410 u_comp_l = u(k,j,i) + u(k,j,i-1) - gu 4411 flux_l = u_comp_l * ( & 4412 ( 37.0_wp * ibit11_l * adv_mom_5 & 4413 + 7.0_wp * ibit10_l * adv_mom_3 & 4414 + ibit9_l * adv_mom_1 & 4415 ) * & 4416 ( u(k,j,i) + u(k,j,i-1) ) & 4417 - ( 8.0_wp * ibit11_l * adv_mom_5 & 4418 + ibit10_l * adv_mom_3 & 4419 ) * & 4420 ( u(k,j,i+1) + u(k,j,i-2) ) & 4421 + ( ibit11_l * adv_mom_5 & 4422 ) * & 4423 ( u(k,j,i+2) + u(k,j,i-3) ) & 4424 ) 4425 4426 diss_l = - ABS( u_comp_l ) * ( & 4427 ( 10.0_wp * ibit11_l * adv_mom_5 & 4428 + 3.0_wp * ibit10_l * adv_mom_3 & 4429 + ibit9_l * adv_mom_1 & 4430 ) * & 4431 ( u(k,j,i) - u(k,j,i-1) ) & 4432 - ( 5.0_wp * ibit11_l * adv_mom_5 & 4433 + ibit10_l * adv_mom_3 & 4434 ) * & 4435 ( u(k,j,i+1) - u(k,j,i-2) ) & 4436 + ( ibit11_l * adv_mom_5 & 4437 ) * & 4438 ( u(k,j,i+2) - u(k,j,i-3) ) & 4439 ) 4440 #else 4441 flux_l = swap_flux_x_local_u(k,j) 4442 diss_l = swap_diss_x_local_u(k,j) 4443 #endif 4444 4067 4445 ibit14 = REAL( IBITS(advc_flags_1(k,j,i),14,1), KIND = wp ) 4068 4446 ibit13 = REAL( IBITS(advc_flags_1(k,j,i),13,1), KIND = wp ) 4069 4447 ibit12 = REAL( IBITS(advc_flags_1(k,j,i),12,1), KIND = wp ) 4070 4448 4071 v_comp 4072 flux_n (k) = v_comp * (&4449 v_comp = v(k,j+1,i) + v(k,j+1,i-1) - gv 4450 flux_n = v_comp * ( & 4073 4451 ( 37.0_wp * ibit14 * adv_mom_5 & 4074 4452 + 7.0_wp * ibit13 * adv_mom_3 & … … 4085 4463 ) 4086 4464 4087 diss_n (k) = - ABS ( v_comp ) * (&4465 diss_n = - ABS ( v_comp ) * ( & 4088 4466 ( 10.0_wp * ibit14 * adv_mom_5 & 4089 4467 + 3.0_wp * ibit13 * adv_mom_3 & … … 4099 4477 ( u(k,j+3,i) - u(k,j-2,i) ) & 4100 4478 ) 4479 4480 #ifdef _OPENACC 4481 ! 4482 !-- Recompute the south fluxes. 4483 ibit14_s = REAL( IBITS(advc_flags_1(k,j-1,i),14,1), KIND = wp ) 4484 ibit13_s = REAL( IBITS(advc_flags_1(k,j-1,i),13,1), KIND = wp ) 4485 ibit12_s = REAL( IBITS(advc_flags_1(k,j-1,i),12,1), KIND = wp ) 4486 4487 v_comp_s = v(k,j,i) + v(k,j,i-1) - gv 4488 flux_s = v_comp_s * ( & 4489 ( 37.0_wp * ibit14_s * adv_mom_5 & 4490 + 7.0_wp * ibit13_s * adv_mom_3 & 4491 + ibit12_s * adv_mom_1 & 4492 ) * & 4493 ( u(k,j,i) + u(k,j-1,i) ) & 4494 - ( 8.0_wp * ibit14_s * adv_mom_5 & 4495 + ibit13_s * adv_mom_3 & 4496 ) * & 4497 ( u(k,j+1,i) + u(k,j-2,i) ) & 4498 + ( ibit14_s * adv_mom_5 & 4499 ) * & 4500 ( u(k,j+2,i) + u(k,j-3,i) ) & 4501 ) 4502 4503 diss_s = - ABS ( v_comp_s ) * ( & 4504 ( 10.0_wp * ibit14_s * adv_mom_5 & 4505 + 3.0_wp * ibit13_s * adv_mom_3 & 4506 + ibit12_s * adv_mom_1 & 4507 ) * & 4508 ( u(k,j,i) - u(k,j-1,i) ) & 4509 - ( 5.0_wp * ibit14_s * adv_mom_5 & 4510 + ibit13_s * adv_mom_3 & 4511 ) * & 4512 ( u(k,j+1,i) - u(k,j-2,i) ) & 4513 + ( ibit14_s * adv_mom_5 & 4514 ) * & 4515 ( u(k,j+2,i) - u(k,j-3,i) ) & 4516 ) 4517 #else 4518 flux_s = swap_flux_y_local_u(k) 4519 diss_s = swap_diss_y_local_u(k) 4520 #endif 4521 4101 4522 ! 4102 4523 !-- k index has to be modified near bottom and top, else array … … 4110 4531 k_mm = k - 2 * ibit17 4111 4532 4112 w_comp 4113 flux_t (k) = w_comp * rho_air_zw(k) * (&4533 w_comp = w(k,j,i) + w(k,j,i-1) 4534 flux_t = w_comp * rho_air_zw(k) * ( & 4114 4535 ( 37.0_wp * ibit17 * adv_mom_5 & 4115 4536 + 7.0_wp * ibit16 * adv_mom_3 & … … 4126 4547 ) 4127 4548 4128 diss_t (k) = - ABS( w_comp ) * rho_air_zw(k) * (&4549 diss_t = - ABS( w_comp ) * rho_air_zw(k) * ( & 4129 4550 ( 10.0_wp * ibit17 * adv_mom_5 & 4130 4551 + 3.0_wp * ibit16 * adv_mom_3 & … … 4144 4565 !-- correction is needed to overcome numerical instabilities caused 4145 4566 !-- by a not sufficient reduction of divergences near topography. 4146 div = ( ( u_comp (k) * ( ibit9 + ibit10 + ibit11 )&4567 div = ( ( u_comp * ( ibit9 + ibit10 + ibit11 ) & 4147 4568 - ( u(k,j,i) + u(k,j,i-1) ) & 4148 4569 * ( & … … 4173 4594 4174 4595 tend(k,j,i) = tend(k,j,i) - ( & 4175 ( flux_r(k) + diss_r(k) & 4176 - swap_flux_x_local_u(k,j) - swap_diss_x_local_u(k,j) ) * ddx & 4177 + ( flux_n(k) + diss_n(k) & 4178 - swap_flux_y_local_u(k) - swap_diss_y_local_u(k) ) * ddy & 4179 + ( ( flux_t(k) + diss_t(k) ) & 4180 - ( flux_d + diss_d ) & 4181 ) * drho_air(k) * ddzw(k) & 4596 ( ( flux_r + diss_r ) & 4597 - ( flux_l + diss_l ) ) * ddx & 4598 + ( ( flux_n + diss_n ) & 4599 - ( flux_s + diss_s ) ) * ddy & 4600 + ( ( flux_t + diss_t ) & 4601 - ( flux_d + diss_d ) ) * drho_air(k) * ddzw(k) & 4182 4602 ) + div * u(k,j,i) 4183 4603 4184 swap_flux_x_local_u(k,j) = flux_r(k) 4185 swap_diss_x_local_u(k,j) = diss_r(k) 4186 swap_flux_y_local_u(k) = flux_n(k) 4187 swap_diss_y_local_u(k) = diss_n(k) 4188 flux_d = flux_t(k) 4189 diss_d = diss_t(k) 4604 #ifndef _OPENACC 4605 swap_flux_x_local_u(k,j) = flux_r 4606 swap_diss_x_local_u(k,j) = diss_r 4607 swap_flux_y_local_u(k) = flux_n 4608 swap_diss_y_local_u(k) = diss_n 4609 #endif 4610 flux_d = flux_t 4611 diss_d = diss_t 4190 4612 ! 4191 4613 !-- Statistical Evaluation of u'u'. The factor has to be applied 4192 4614 !-- for right evaluation when gallilei_trans = .T. . 4615 !$ACC ATOMIC 4193 4616 sums_us2_ws_l(k,tn) = sums_us2_ws_l(k,tn) & 4194 + ( flux_r (k)&4195 * ( u_comp (k) - 2.0_wp * hom(k,1,1,0) )&4196 / ( u_comp (k) - gu + SIGN( 1.0E-20_wp, u_comp(k) - gu ) )&4197 + diss_r (k)&4198 * ABS( u_comp (k) - 2.0_wp * hom(k,1,1,0) )&4199 / ( ABS( u_comp (k) - gu ) + 1.0E-20_wp )&4617 + ( flux_r & 4618 * ( u_comp - 2.0_wp * hom(k,1,1,0) ) & 4619 / ( u_comp - gu + SIGN( 1.0E-20_wp, u_comp - gu ) ) & 4620 + diss_r & 4621 * ABS( u_comp - 2.0_wp * hom(k,1,1,0) ) & 4622 / ( ABS( u_comp - gu ) + 1.0E-20_wp ) & 4200 4623 ) * weight_substep(intermediate_timestep_count) 4201 4624 ! 4202 4625 !-- Statistical Evaluation of w'u'. 4626 !$ACC ATOMIC 4203 4627 sums_wsus_ws_l(k,tn) = sums_wsus_ws_l(k,tn) & 4204 + ( flux_t (k)&4628 + ( flux_t & 4205 4629 * ( w_comp - 2.0_wp * hom(k,1,3,0) ) & 4206 4630 / ( w_comp + SIGN( 1.0E-20_wp, w_comp ) ) & 4207 + diss_t (k)&4631 + diss_t & 4208 4632 * ABS( w_comp - 2.0_wp * hom(k,1,3,0) ) & 4209 4633 / ( ABS( w_comp ) + 1.0E-20_wp ) & … … 4214 4638 DO k = nzb_max+1, nzt 4215 4639 4216 u_comp (k)= u(k,j,i+1) + u(k,j,i)4217 flux_r (k) = ( u_comp(k) - gu ) * (&4640 u_comp = u(k,j,i+1) + u(k,j,i) 4641 flux_r = ( u_comp - gu ) * ( & 4218 4642 37.0_wp * ( u(k,j,i+1) + u(k,j,i) ) & 4219 4643 - 8.0_wp * ( u(k,j,i+2) + u(k,j,i-1) ) & 4220 4644 + ( u(k,j,i+3) + u(k,j,i-2) ) ) * adv_mom_5 4221 diss_r (k) = - ABS( u_comp(k) - gu ) * (&4645 diss_r = - ABS( u_comp - gu ) * ( & 4222 4646 10.0_wp * ( u(k,j,i+1) - u(k,j,i) ) & 4223 4647 - 5.0_wp * ( u(k,j,i+2) - u(k,j,i-1) ) & 4224 4648 + ( u(k,j,i+3) - u(k,j,i-2) ) ) * adv_mom_5 4225 4649 4226 v_comp = v(k,j+1,i) + v(k,j+1,i-1) - gv 4227 flux_n(k) = v_comp * ( & 4650 #ifdef _OPENACC 4651 ! 4652 !-- Recompute the left fluxes. 4653 u_comp_l = u(k,j,i) + u(k,j,i-1) - gu 4654 flux_l = u_comp_l * ( & 4655 37.0_wp * ( u(k,j,i) + u(k,j,i-1) ) & 4656 - 8.0_wp * ( u(k,j,i+1) + u(k,j,i-2) ) & 4657 + ( u(k,j,i+2) + u(k,j,i-3) ) ) * adv_mom_5 4658 diss_l = - ABS(u_comp_l) * ( & 4659 10.0_wp * ( u(k,j,i) - u(k,j,i-1) ) & 4660 - 5.0_wp * ( u(k,j,i+1) - u(k,j,i-2) ) & 4661 + ( u(k,j,i+2) - u(k,j,i-3) ) ) * adv_mom_5 4662 #else 4663 flux_l = swap_flux_x_local_u(k,j) 4664 diss_l = swap_diss_x_local_u(k,j) 4665 #endif 4666 4667 v_comp = v(k,j+1,i) + v(k,j+1,i-1) - gv 4668 flux_n = v_comp * ( & 4228 4669 37.0_wp * ( u(k,j+1,i) + u(k,j,i) ) & 4229 4670 - 8.0_wp * ( u(k,j+2,i) + u(k,j-1,i) ) & 4230 4671 + ( u(k,j+3,i) + u(k,j-2,i) ) ) * adv_mom_5 4231 diss_n (k) = - ABS( v_comp ) * (&4672 diss_n = - ABS( v_comp ) * ( & 4232 4673 10.0_wp * ( u(k,j+1,i) - u(k,j,i) ) & 4233 4674 - 5.0_wp * ( u(k,j+2,i) - u(k,j-1,i) ) & 4234 4675 + ( u(k,j+3,i) - u(k,j-2,i) ) ) * adv_mom_5 4676 4677 #ifdef _OPENACC 4678 ! 4679 !-- Recompute the south fluxes. 4680 v_comp_s = v(k,j,i) + v(k,j,i-1) - gv 4681 flux_s = v_comp_s * ( & 4682 37.0_wp * ( u(k,j,i) + u(k,j-1,i) ) & 4683 - 8.0_wp * ( u(k,j+1,i) + u(k,j-2,i) ) & 4684 + ( u(k,j+2,i) + u(k,j-3,i) ) ) * adv_mom_5 4685 diss_s = - ABS( v_comp_s ) * ( & 4686 10.0_wp * ( u(k,j,i) - u(k,j-1,i) ) & 4687 - 5.0_wp * ( u(k,j+1,i) - u(k,j-2,i) ) & 4688 + ( u(k,j+2,i) - u(k,j-3,i) ) ) * adv_mom_5 4689 #else 4690 flux_s = swap_flux_y_local_u(k) 4691 diss_s = swap_diss_y_local_u(k) 4692 #endif 4693 4235 4694 ! 4236 4695 !-- k index has to be modified near bottom and top, else array … … 4244 4703 k_mm = k - 2 * ibit17 4245 4704 4246 w_comp 4247 flux_t (k) = w_comp * rho_air_zw(k) * (&4705 w_comp = w(k,j,i) + w(k,j,i-1) 4706 flux_t = w_comp * rho_air_zw(k) * ( & 4248 4707 ( 37.0_wp * ibit17 * adv_mom_5 & 4249 4708 + 7.0_wp * ibit16 * adv_mom_3 & … … 4260 4719 ) 4261 4720 4262 diss_t (k) = - ABS( w_comp ) * rho_air_zw(k) * (&4721 diss_t = - ABS( w_comp ) * rho_air_zw(k) * ( & 4263 4722 ( 10.0_wp * ibit17 * adv_mom_5 & 4264 4723 + 3.0_wp * ibit16 * adv_mom_3 & … … 4278 4737 !-- correction is needed to overcome numerical instabilities caused 4279 4738 !-- by a not sufficient reduction of divergences near topography. 4280 div = ( ( u_comp (k)- ( u(k,j,i) + u(k,j,i-1) ) ) * ddx &4739 div = ( ( u_comp - ( u(k,j,i) + u(k,j,i-1) ) ) * ddx & 4281 4740 + ( v_comp + gv - ( v(k,j,i) + v(k,j,i-1 ) ) ) * ddy & 4282 4741 + ( w_comp * rho_air_zw(k) - & … … 4286 4745 4287 4746 tend(k,j,i) = tend(k,j,i) - ( & 4288 ( flux_r(k) + diss_r(k) & 4289 - swap_flux_x_local_u(k,j) - swap_diss_x_local_u(k,j) ) * ddx & 4290 + ( flux_n(k) + diss_n(k) & 4291 - swap_flux_y_local_u(k) - swap_diss_y_local_u(k) ) * ddy & 4292 + ( ( flux_t(k) + diss_t(k) ) & 4293 - ( flux_d + diss_d ) & 4294 ) * drho_air(k) * ddzw(k) & 4747 ( ( flux_r + diss_r ) & 4748 - ( flux_l + diss_l ) ) * ddx & 4749 + ( ( flux_n + diss_n ) & 4750 - ( flux_s + diss_s ) ) * ddy & 4751 + ( ( flux_t + diss_t ) & 4752 - ( flux_d + diss_d ) ) * drho_air(k) * ddzw(k) & 4295 4753 ) + div * u(k,j,i) 4296 4754 4297 swap_flux_x_local_u(k,j) = flux_r(k) 4298 swap_diss_x_local_u(k,j) = diss_r(k) 4299 swap_flux_y_local_u(k) = flux_n(k) 4300 swap_diss_y_local_u(k) = diss_n(k) 4301 flux_d = flux_t(k) 4302 diss_d = diss_t(k) 4755 #ifndef _OPENACC 4756 swap_flux_x_local_u(k,j) = flux_r 4757 swap_diss_x_local_u(k,j) = diss_r 4758 swap_flux_y_local_u(k) = flux_n 4759 swap_diss_y_local_u(k) = diss_n 4760 #endif 4761 flux_d = flux_t 4762 diss_d = diss_t 4303 4763 ! 4304 4764 !-- Statistical Evaluation of u'u'. The factor has to be applied 4305 4765 !-- for right evaluation when gallilei_trans = .T. . 4766 !$ACC ATOMIC 4306 4767 sums_us2_ws_l(k,tn) = sums_us2_ws_l(k,tn) & 4307 + ( flux_r (k)&4308 * ( u_comp (k) - 2.0_wp * hom(k,1,1,0) )&4309 / ( u_comp (k) - gu + SIGN( 1.0E-20_wp, u_comp(k) - gu ) )&4310 + diss_r (k)&4311 * ABS( u_comp (k) - 2.0_wp * hom(k,1,1,0) )&4312 / ( ABS( u_comp (k) - gu ) + 1.0E-20_wp )&4768 + ( flux_r & 4769 * ( u_comp - 2.0_wp * hom(k,1,1,0) ) & 4770 / ( u_comp - gu + SIGN( 1.0E-20_wp, u_comp - gu ) ) & 4771 + diss_r & 4772 * ABS( u_comp - 2.0_wp * hom(k,1,1,0) ) & 4773 / ( ABS( u_comp - gu ) + 1.0E-20_wp ) & 4313 4774 ) * weight_substep(intermediate_timestep_count) 4314 4775 ! 4315 4776 !-- Statistical Evaluation of w'u'. 4777 !$ACC ATOMIC 4316 4778 sums_wsus_ws_l(k,tn) = sums_wsus_ws_l(k,tn) & 4317 + ( flux_t (k)&4779 + ( flux_t & 4318 4780 * ( w_comp - 2.0_wp * hom(k,1,3,0) ) & 4319 4781 / ( w_comp + SIGN( 1.0E-20_wp, w_comp ) ) & 4320 + diss_t (k)&4782 + diss_t & 4321 4783 * ABS( w_comp - 2.0_wp * hom(k,1,3,0) ) & 4322 4784 / ( ABS( w_comp ) + 1.0E-20_wp ) & … … 4369 4831 REAL(wp) :: ibit19 !< flag indicating 3rd-order scheme along x-direction 4370 4832 REAL(wp) :: ibit20 !< flag indicating 5th-order scheme along x-direction 4833 #ifdef _OPENACC 4834 REAL(wp) :: ibit18_l !< flag indicating 1st-order scheme along x-direction 4835 REAL(wp) :: ibit19_l !< flag indicating 3rd-order scheme along x-direction 4836 REAL(wp) :: ibit20_l !< flag indicating 5th-order scheme along x-direction 4837 #endif 4371 4838 REAL(wp) :: ibit21 !< flag indicating 1st-order scheme along y-direction 4372 4839 REAL(wp) :: ibit22 !< flag indicating 3rd-order scheme along y-direction 4373 4840 REAL(wp) :: ibit23 !< flag indicating 5th-order scheme along y-direction 4841 #ifdef _OPENACC 4842 REAL(wp) :: ibit21_s !< flag indicating 1st-order scheme along y-direction 4843 REAL(wp) :: ibit22_s !< flag indicating 3rd-order scheme along y-direction 4844 REAL(wp) :: ibit23_s !< flag indicating 5th-order scheme along y-direction 4845 #endif 4374 4846 REAL(wp) :: ibit24 !< flag indicating 1st-order scheme along z-direction 4375 4847 REAL(wp) :: ibit25 !< flag indicating 3rd-order scheme along z-direction … … 4381 4853 REAL(wp) :: gv !< Galilei-transformation velocity along y 4382 4854 REAL(wp) :: u_comp !< advection velocity along x 4855 #ifdef _OPENACC 4856 REAL(wp) :: u_comp_l !< advection velocity along x 4857 #endif 4383 4858 REAL(wp) :: w_comp !< advection velocity along z 4384 4859 4860 REAL(wp) :: diss_s !< discretized artificial dissipation at southward-side of the grid box 4861 REAL(wp) :: flux_s !< discretized 6th-order flux at southward-side of the grid box 4862 #ifndef _OPENACC 4385 4863 REAL(wp), DIMENSION(nzb+1:nzt) :: swap_diss_y_local_v !< discretized artificial dissipation at southward-side of the grid box 4386 4864 REAL(wp), DIMENSION(nzb+1:nzt) :: swap_flux_y_local_v !< discretized 6th-order flux at southward-side of the grid box 4865 #endif 4387 4866 4867 REAL(wp) :: diss_l !< discretized artificial dissipation at leftward-side of the grid box 4868 REAL(wp) :: flux_l !< discretized 6th-order flux at leftward-side of the grid box 4869 #ifndef _OPENACC 4388 4870 REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) :: swap_diss_x_local_v !< discretized artificial dissipation at leftward-side of the grid box 4389 4871 REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) :: swap_flux_x_local_v !< discretized 6th-order flux at leftward-side of the grid box 4872 #endif 4390 4873 4391 REAL(wp), DIMENSION(nzb:nzt) :: diss_n !< discretized artificial dissipation at northward-side of the grid box 4392 REAL(wp), DIMENSION(nzb:nzt) :: diss_r !< discretized artificial dissipation at rightward-side of the grid box 4393 REAL(wp), DIMENSION(nzb:nzt) :: diss_t !< discretized artificial dissipation at top of the grid box 4394 REAL(wp), DIMENSION(nzb:nzt) :: flux_n !< discretized 6th-order flux at northward-side of the grid box 4395 REAL(wp), DIMENSION(nzb:nzt) :: flux_r !< discretized 6th-order flux at rightward-side of the grid box 4396 REAL(wp), DIMENSION(nzb:nzt) :: flux_t !< discretized 6th-order flux at top of the grid box 4397 REAL(wp), DIMENSION(nzb:nzt) :: v_comp !< advection velocity along y 4874 REAL(wp) :: diss_n !< discretized artificial dissipation at northward-side of the grid box 4875 REAL(wp) :: diss_r !< discretized artificial dissipation at rightward-side of the grid box 4876 REAL(wp) :: diss_t !< discretized artificial dissipation at top of the grid box 4877 REAL(wp) :: flux_n !< discretized 6th-order flux at northward-side of the grid box 4878 REAL(wp) :: flux_r !< discretized 6th-order flux at rightward-side of the grid box 4879 REAL(wp) :: flux_t !< discretized 6th-order flux at top of the grid box 4880 REAL(wp) :: v_comp !< advection velocity along y 4881 #ifdef _OPENACC 4882 REAL(wp) :: v_comp_s !< 4883 #endif 4398 4884 4399 4885 gu = 2.0_wp * u_gtrans 4400 4886 gv = 2.0_wp * v_gtrans 4887 4888 #ifndef _OPENACC 4401 4889 ! 4402 4890 !-- First compute the whole left boundary of the processor domain … … 4457 4945 4458 4946 ENDDO 4459 4947 #endif 4948 4949 !$ACC PARALLEL LOOP COLLAPSE(2) FIRSTPRIVATE(tn, gu, gv) & 4950 !$ACC PRIVATE(i, j, k, k_mm, k_pp, k_ppp) & 4951 !$ACC PRIVATE(ibit18, ibit19, ibit20, ibit21, ibit22, ibit23) & 4952 !$ACC PRIVATE(ibit18_l, ibit19_l, ibit20_l) & 4953 !$ACC PRIVATE(ibit21_s, ibit22_s, ibit23_s) & 4954 !$ACC PRIVATE(ibit24, ibit25, ibit26) & 4955 !$ACC PRIVATE(flux_r, diss_r, flux_l, diss_l) & 4956 !$ACC PRIVATE(flux_n, diss_n, flux_s, diss_s) & 4957 !$ACC PRIVATE(flux_t, diss_t, flux_d, diss_d) & 4958 !$ACC PRIVATE(div, u_comp, u_comp_l, v_comp, v_comp_s, w_comp) & 4959 !$ACC PRESENT(advc_flags_1) & 4960 !$ACC PRESENT(u, v, w) & 4961 !$ACC PRESENT(drho_air, rho_air_zw, ddzw) & 4962 !$ACC PRESENT(tend) & 4963 !$ACC PRESENT(hom(nzb+1:nzb_max,1,2:3,0)) & 4964 !$ACC PRESENT(weight_substep(intermediate_timestep_count)) & 4965 !$ACC PRESENT(sums_vs2_ws_l, sums_wsvs_ws_l) 4460 4966 DO i = nxl, nxr 4461 4967 4968 #ifndef _OPENACC 4462 4969 j = nysv 4463 4970 DO k = nzb+1, nzb_max … … 4467 4974 ibit21 = REAL( IBITS(advc_flags_1(k,j-1,i),21,1), KIND = wp ) 4468 4975 4469 v_comp (k)= v(k,j,i) + v(k,j-1,i) - gv4470 swap_flux_y_local_v(k) = v_comp (k) * (&4976 v_comp = v(k,j,i) + v(k,j-1,i) - gv 4977 swap_flux_y_local_v(k) = v_comp * ( & 4471 4978 ( 37.0_wp * ibit23 * adv_mom_5 & 4472 4979 + 7.0_wp * ibit22 * adv_mom_3 & … … 4483 4990 ) 4484 4991 4485 swap_diss_y_local_v(k) = - ABS( v_comp (k) ) * (&4992 swap_diss_y_local_v(k) = - ABS( v_comp ) * ( & 4486 4993 ( 10.0_wp * ibit23 * adv_mom_5 & 4487 4994 + 3.0_wp * ibit22 * adv_mom_3 & … … 4502 5009 DO k = nzb_max+1, nzt 4503 5010 4504 v_comp (k)= v(k,j,i) + v(k,j-1,i) - gv4505 swap_flux_y_local_v(k) = v_comp (k) * (&5011 v_comp = v(k,j,i) + v(k,j-1,i) - gv 5012 swap_flux_y_local_v(k) = v_comp * ( & 4506 5013 37.0_wp * ( v(k,j,i) + v(k,j-1,i) ) & 4507 5014 - 8.0_wp * ( v(k,j+1,i) + v(k,j-2,i) ) & 4508 5015 + ( v(k,j+2,i) + v(k,j-3,i) ) ) * adv_mom_5 4509 swap_diss_y_local_v(k) = - ABS( v_comp (k) ) * (&5016 swap_diss_y_local_v(k) = - ABS( v_comp ) * ( & 4510 5017 10.0_wp * ( v(k,j,i) - v(k,j-1,i) ) & 4511 5018 - 5.0_wp * ( v(k,j+1,i) - v(k,j-2,i) ) & … … 4513 5020 4514 5021 ENDDO 5022 #endif 4515 5023 4516 5024 DO j = nysv, nyn 4517 5025 4518 flux_t(0) = 0.0_wp4519 diss_t(0) = 0.0_wp4520 5026 flux_d = 0.0_wp 4521 5027 diss_d = 0.0_wp … … 4527 5033 ibit18 = REAL( IBITS(advc_flags_1(k,j,i),18,1), KIND = wp ) 4528 5034 4529 u_comp 4530 flux_r (k) = u_comp * (&5035 u_comp = u(k,j-1,i+1) + u(k,j,i+1) - gu 5036 flux_r = u_comp * ( & 4531 5037 ( 37.0_wp * ibit20 * adv_mom_5 & 4532 5038 + 7.0_wp * ibit19 * adv_mom_3 & … … 4543 5049 ) 4544 5050 4545 diss_r (k) = - ABS( u_comp ) * (&5051 diss_r = - ABS( u_comp ) * ( & 4546 5052 ( 10.0_wp * ibit20 * adv_mom_5 & 4547 5053 + 3.0_wp * ibit19 * adv_mom_3 & … … 4558 5064 ) 4559 5065 5066 #ifdef _OPENACC 5067 ! 5068 !-- Recompute the left fluxes. 5069 ibit20_l = REAL( IBITS(advc_flags_1(k,j,i-1),20,1), KIND = wp ) 5070 ibit19_l = REAL( IBITS(advc_flags_1(k,j,i-1),19,1), KIND = wp ) 5071 ibit18_l = REAL( IBITS(advc_flags_1(k,j,i-1),18,1), KIND = wp ) 5072 5073 u_comp_l = u(k,j-1,i) + u(k,j,i) - gu 5074 flux_l = u_comp_l * ( & 5075 ( 37.0_wp * ibit20_l * adv_mom_5 & 5076 + 7.0_wp * ibit19_l * adv_mom_3 & 5077 + ibit18_l * adv_mom_1 & 5078 ) * & 5079 ( v(k,j,i) + v(k,j,i-1) ) & 5080 - ( 8.0_wp * ibit20_l * adv_mom_5 & 5081 + ibit19_l * adv_mom_3 & 5082 ) * & 5083 ( v(k,j,i+1) + v(k,j,i-2) ) & 5084 + ( ibit20_l * adv_mom_5 & 5085 ) * & 5086 ( v(k,j,i+2) + v(k,j,i-3) ) & 5087 ) 5088 5089 diss_l = - ABS( u_comp_l ) * ( & 5090 ( 10.0_wp * ibit20_l * adv_mom_5 & 5091 + 3.0_wp * ibit19_l * adv_mom_3 & 5092 + ibit18_l * adv_mom_1 & 5093 ) * & 5094 ( v(k,j,i) - v(k,j,i-1) ) & 5095 - ( 5.0_wp * ibit20_l * adv_mom_5 & 5096 + ibit19_l * adv_mom_3 & 5097 ) * & 5098 ( v(k,j,i+1) - v(k,j,i-2) ) & 5099 + ( ibit20_l * adv_mom_5 & 5100 ) * & 5101 ( v(k,j,i+2) - v(k,j,i-3) ) & 5102 ) 5103 #else 5104 flux_l = swap_flux_x_local_v(k,j) 5105 diss_l = swap_diss_x_local_v(k,j) 5106 #endif 5107 4560 5108 ibit23 = REAL( IBITS(advc_flags_1(k,j,i),23,1), KIND = wp ) 4561 5109 ibit22 = REAL( IBITS(advc_flags_1(k,j,i),22,1), KIND = wp ) 4562 5110 ibit21 = REAL( IBITS(advc_flags_1(k,j,i),21,1), KIND = wp ) 4563 5111 4564 v_comp (k)= v(k,j+1,i) + v(k,j,i)4565 flux_n (k) = ( v_comp(k) - gv ) * (&5112 v_comp = v(k,j+1,i) + v(k,j,i) 5113 flux_n = ( v_comp - gv ) * ( & 4566 5114 ( 37.0_wp * ibit23 * adv_mom_5 & 4567 5115 + 7.0_wp * ibit22 * adv_mom_3 & … … 4578 5126 ) 4579 5127 4580 diss_n (k) = - ABS( v_comp(k) - gv ) * (&5128 diss_n = - ABS( v_comp - gv ) * ( & 4581 5129 ( 10.0_wp * ibit23 * adv_mom_5 & 4582 5130 + 3.0_wp * ibit22 * adv_mom_3 & … … 4592 5140 ( v(k,j+3,i) - v(k,j-2,i) ) & 4593 5141 ) 5142 5143 #ifdef _OPENACC 5144 ! 5145 !-- Recompute the south fluxes. 5146 ibit23_s = REAL( IBITS(advc_flags_1(k,j-1,i),23,1), KIND = wp ) 5147 ibit22_s = REAL( IBITS(advc_flags_1(k,j-1,i),22,1), KIND = wp ) 5148 ibit21_s = REAL( IBITS(advc_flags_1(k,j-1,i),21,1), KIND = wp ) 5149 5150 v_comp_s = v(k,j,i) + v(k,j-1,i) - gv 5151 flux_s = v_comp_s * ( & 5152 ( 37.0_wp * ibit23_s * adv_mom_5 & 5153 + 7.0_wp * ibit22_s * adv_mom_3 & 5154 + ibit21_s * adv_mom_1 & 5155 ) * & 5156 ( v(k,j,i) + v(k,j-1,i) ) & 5157 - ( 8.0_wp * ibit23_s * adv_mom_5 & 5158 + ibit22_s * adv_mom_3 & 5159 ) * & 5160 ( v(k,j+1,i) + v(k,j-2,i) ) & 5161 + ( ibit23_s * adv_mom_5 & 5162 ) * & 5163 ( v(k,j+2,i) + v(k,j-3,i) ) & 5164 ) 5165 5166 diss_s = - ABS( v_comp_s ) * ( & 5167 ( 10.0_wp * ibit23_s * adv_mom_5 & 5168 + 3.0_wp * ibit22_s * adv_mom_3 & 5169 + ibit21_s * adv_mom_1 & 5170 ) * & 5171 ( v(k,j,i) - v(k,j-1,i) ) & 5172 - ( 5.0_wp * ibit23_s * adv_mom_5 & 5173 + ibit22_s * adv_mom_3 & 5174 ) * & 5175 ( v(k,j+1,i) - v(k,j-2,i) ) & 5176 + ( ibit23_s * adv_mom_5 & 5177 ) * & 5178 ( v(k,j+2,i) - v(k,j-3,i) ) & 5179 ) 5180 #else 5181 flux_s = swap_flux_y_local_v(k) 5182 diss_s = swap_diss_y_local_v(k) 5183 #endif 5184 4594 5185 ! 4595 5186 !-- k index has to be modified near bottom and top, else array … … 4603 5194 k_mm = k - 2 * ibit26 4604 5195 4605 w_comp 4606 flux_t (k) = w_comp * rho_air_zw(k) * (&5196 w_comp = w(k,j-1,i) + w(k,j,i) 5197 flux_t = w_comp * rho_air_zw(k) * ( & 4607 5198 ( 37.0_wp * ibit26 * adv_mom_5 & 4608 5199 + 7.0_wp * ibit25 * adv_mom_3 & … … 4619 5210 ) 4620 5211 4621 diss_t (k) = - ABS( w_comp ) * rho_air_zw(k) * (&5212 diss_t = - ABS( w_comp ) * rho_air_zw(k) * ( & 4622 5213 ( 10.0_wp * ibit26 * adv_mom_5 & 4623 5214 + 3.0_wp * ibit25 * adv_mom_3 & … … 4646 5237 ) & 4647 5238 ) * ddx & 4648 + ( v_comp (k)&5239 + ( v_comp & 4649 5240 * ( ibit21 + ibit22 + ibit23 ) & 4650 5241 - ( v(k,j,i) + v(k,j-1,i) ) & … … 4668 5259 4669 5260 tend(k,j,i) = tend(k,j,i) - ( & 4670 ( flux_r(k) + diss_r(k) & 4671 - swap_flux_x_local_v(k,j) - swap_diss_x_local_v(k,j) & 4672 ) * ddx & 4673 + ( flux_n(k) + diss_n(k) & 4674 - swap_flux_y_local_v(k) - swap_diss_y_local_v(k) & 4675 ) * ddy & 4676 + ( ( flux_t(k) + diss_t(k) ) & 4677 - ( flux_d + diss_d ) & 4678 ) * drho_air(k) * ddzw(k) & 5261 ( ( flux_r + diss_r ) & 5262 - ( flux_l + diss_l ) ) * ddx & 5263 + ( ( flux_n + diss_n ) & 5264 - ( flux_s + diss_s ) ) * ddy & 5265 + ( ( flux_t + diss_t ) & 5266 - ( flux_d + diss_d ) ) * drho_air(k) * ddzw(k) & 4679 5267 ) + v(k,j,i) * div 4680 5268 4681 swap_flux_x_local_v(k,j) = flux_r(k) 4682 swap_diss_x_local_v(k,j) = diss_r(k) 4683 swap_flux_y_local_v(k) = flux_n(k) 4684 swap_diss_y_local_v(k) = diss_n(k) 4685 flux_d = flux_t(k) 4686 diss_d = diss_t(k) 5269 #ifndef _OPENACC 5270 swap_flux_x_local_v(k,j) = flux_r 5271 swap_diss_x_local_v(k,j) = diss_r 5272 swap_flux_y_local_v(k) = flux_n 5273 swap_diss_y_local_v(k) = diss_n 5274 #endif 5275 flux_d = flux_t 5276 diss_d = diss_t 4687 5277 4688 5278 ! 4689 5279 !-- Statistical Evaluation of v'v'. The factor has to be applied 4690 5280 !-- for right evaluation when gallilei_trans = .T. . 5281 !$ACC ATOMIC 4691 5282 sums_vs2_ws_l(k,tn) = sums_vs2_ws_l(k,tn) & 4692 + ( flux_n (k)&4693 * ( v_comp (k) - 2.0_wp * hom(k,1,2,0) )&4694 / ( v_comp (k) - gv + SIGN( 1.0E-20_wp, v_comp(k) - gv ) )&4695 + diss_n (k)&4696 * ABS( v_comp (k) - 2.0_wp * hom(k,1,2,0) )&4697 / ( ABS( v_comp (k) - gv ) + 1.0E-20_wp )&5283 + ( flux_n & 5284 * ( v_comp - 2.0_wp * hom(k,1,2,0) ) & 5285 / ( v_comp - gv + SIGN( 1.0E-20_wp, v_comp - gv ) ) & 5286 + diss_n & 5287 * ABS( v_comp - 2.0_wp * hom(k,1,2,0) ) & 5288 / ( ABS( v_comp - gv ) + 1.0E-20_wp ) & 4698 5289 ) * weight_substep(intermediate_timestep_count) 4699 5290 ! 4700 5291 !-- Statistical Evaluation of w'u'. 5292 !$ACC ATOMIC 4701 5293 sums_wsvs_ws_l(k,tn) = sums_wsvs_ws_l(k,tn) & 4702 + ( flux_t (k)&5294 + ( flux_t & 4703 5295 * ( w_comp - 2.0_wp * hom(k,1,3,0) ) & 4704 5296 / ( w_comp + SIGN( 1.0E-20_wp, w_comp ) ) & 4705 + diss_t (k)&5297 + diss_t & 4706 5298 * ABS( w_comp - 2.0_wp * hom(k,1,3,0) ) & 4707 5299 / ( ABS( w_comp ) + 1.0E-20_wp ) & … … 4712 5304 DO k = nzb_max+1, nzt 4713 5305 4714 u_comp 4715 flux_r (k) = u_comp * (&5306 u_comp = u(k,j-1,i+1) + u(k,j,i+1) - gu 5307 flux_r = u_comp * ( & 4716 5308 37.0_wp * ( v(k,j,i+1) + v(k,j,i) ) & 4717 5309 - 8.0_wp * ( v(k,j,i+2) + v(k,j,i-1) ) & 4718 5310 + ( v(k,j,i+3) + v(k,j,i-2) ) ) * adv_mom_5 4719 5311 4720 diss_r (k) = - ABS( u_comp ) * (&5312 diss_r = - ABS( u_comp ) * ( & 4721 5313 10.0_wp * ( v(k,j,i+1) - v(k,j,i) ) & 4722 5314 - 5.0_wp * ( v(k,j,i+2) - v(k,j,i-1) ) & 4723 5315 + ( v(k,j,i+3) - v(k,j,i-2) ) ) * adv_mom_5 4724 5316 4725 4726 v_comp(k) = v(k,j+1,i) + v(k,j,i) 4727 flux_n(k) = ( v_comp(k) - gv ) * ( & 5317 #ifdef _OPENACC 5318 ! 5319 !-- Recompute the left fluxes. 5320 u_comp_l = u(k,j-1,i) + u(k,j,i) - gu 5321 flux_l = u_comp_l * ( & 5322 37.0_wp * ( v(k,j,i) + v(k,j,i-1) ) & 5323 - 8.0_wp * ( v(k,j,i+1) + v(k,j,i-2) ) & 5324 + ( v(k,j,i+2) + v(k,j,i-3) ) ) * adv_mom_5 5325 diss_l = - ABS( u_comp_l ) * ( & 5326 10.0_wp * ( v(k,j,i) - v(k,j,i-1) ) & 5327 - 5.0_wp * ( v(k,j,i+1) - v(k,j,i-2) ) & 5328 + ( v(k,j,i+2) - v(k,j,i-3) ) ) * adv_mom_5 5329 #else 5330 flux_l = swap_flux_x_local_v(k,j) 5331 diss_l = swap_diss_x_local_v(k,j) 5332 #endif 5333 5334 v_comp = v(k,j+1,i) + v(k,j,i) 5335 flux_n = ( v_comp - gv ) * ( & 4728 5336 37.0_wp * ( v(k,j+1,i) + v(k,j,i) ) & 4729 5337 - 8.0_wp * ( v(k,j+2,i) + v(k,j-1,i) ) & 4730 5338 + ( v(k,j+3,i) + v(k,j-2,i) ) ) * adv_mom_5 4731 5339 4732 diss_n (k) = - ABS( v_comp(k) - gv ) * (&5340 diss_n = - ABS( v_comp - gv ) * ( & 4733 5341 10.0_wp * ( v(k,j+1,i) - v(k,j,i) ) & 4734 5342 - 5.0_wp * ( v(k,j+2,i) - v(k,j-1,i) ) & 4735 5343 + ( v(k,j+3,i) - v(k,j-2,i) ) ) * adv_mom_5 5344 5345 #ifdef _OPENACC 5346 ! 5347 !-- Recompute the south fluxes. 5348 v_comp_s = v(k,j,i) + v(k,j-1,i) - gv 5349 flux_s = v_comp_s * ( & 5350 37.0_wp * ( v(k,j,i) + v(k,j-1,i) ) & 5351 - 8.0_wp * ( v(k,j+1,i) + v(k,j-2,i) ) & 5352 + ( v(k,j+2,i) + v(k,j-3,i) ) ) * adv_mom_5 5353 diss_s = - ABS( v_comp_s ) * ( & 5354 10.0_wp * ( v(k,j,i) - v(k,j-1,i) ) & 5355 - 5.0_wp * ( v(k,j+1,i) - v(k,j-2,i) ) & 5356 + ( v(k,j+2,i) - v(k,j-3,i) ) ) * adv_mom_5 5357 #else 5358 flux_s = swap_flux_y_local_v(k) 5359 diss_s = swap_diss_y_local_v(k) 5360 #endif 5361 4736 5362 ! 4737 5363 !-- k index has to be modified near bottom and top, else array … … 4745 5371 k_mm = k - 2 * ibit26 4746 5372 4747 w_comp 4748 flux_t (k) = w_comp * rho_air_zw(k) * (&5373 w_comp = w(k,j-1,i) + w(k,j,i) 5374 flux_t = w_comp * rho_air_zw(k) * ( & 4749 5375 ( 37.0_wp * ibit26 * adv_mom_5 & 4750 5376 + 7.0_wp * ibit25 * adv_mom_3 & … … 4761 5387 ) 4762 5388 4763 diss_t (k) = - ABS( w_comp ) * rho_air_zw(k) * (&5389 diss_t = - ABS( w_comp ) * rho_air_zw(k) * ( & 4764 5390 ( 10.0_wp * ibit26 * adv_mom_5 & 4765 5391 + 3.0_wp * ibit25 * adv_mom_3 & … … 4780 5406 !-- by a not sufficient reduction of divergences near topography. 4781 5407 div = ( ( u_comp + gu - ( u(k,j-1,i) + u(k,j,i) ) ) * ddx & 4782 + ( v_comp (k)- ( v(k,j,i) + v(k,j-1,i) ) ) * ddy &5408 + ( v_comp - ( v(k,j,i) + v(k,j-1,i) ) ) * ddy & 4783 5409 + ( w_comp * rho_air_zw(k) - & 4784 5410 ( w(k-1,j-1,i) + w(k-1,j,i) ) * rho_air_zw(k-1) & … … 4787 5413 4788 5414 tend(k,j,i) = tend(k,j,i) - ( & 4789 ( flux_r(k) + diss_r(k) & 4790 - swap_flux_x_local_v(k,j) - swap_diss_x_local_v(k,j) & 4791 ) * ddx & 4792 + ( flux_n(k) + diss_n(k) & 4793 - swap_flux_y_local_v(k) - swap_diss_y_local_v(k) & 4794 ) * ddy & 4795 + ( ( flux_t(k) + diss_t(k) ) & 4796 - ( flux_d + diss_d ) & 4797 ) * drho_air(k) * ddzw(k) & 5415 ( ( flux_r + diss_r ) & 5416 - ( flux_l + diss_l ) ) * ddx & 5417 + ( ( flux_n + diss_n ) & 5418 - ( flux_s + diss_s ) ) * ddy & 5419 + ( ( flux_t + diss_t ) & 5420 - ( flux_d + diss_d ) ) * drho_air(k) * ddzw(k) & 4798 5421 ) + v(k,j,i) * div 4799 5422 4800 swap_flux_x_local_v(k,j) = flux_r(k) 4801 swap_diss_x_local_v(k,j) = diss_r(k) 4802 swap_flux_y_local_v(k) = flux_n(k) 4803 swap_diss_y_local_v(k) = diss_n(k) 4804 flux_d = flux_t(k) 4805 diss_d = diss_t(k) 5423 #ifndef _OPENACC 5424 swap_flux_x_local_v(k,j) = flux_r 5425 swap_diss_x_local_v(k,j) = diss_r 5426 swap_flux_y_local_v(k) = flux_n 5427 swap_diss_y_local_v(k) = diss_n 5428 #endif 5429 flux_d = flux_t 5430 diss_d = diss_t 4806 5431 4807 5432 ! 4808 5433 !-- Statistical Evaluation of v'v'. The factor has to be applied 4809 5434 !-- for right evaluation when gallilei_trans = .T. . 5435 !$ACC ATOMIC 4810 5436 sums_vs2_ws_l(k,tn) = sums_vs2_ws_l(k,tn) & 4811 + ( flux_n (k)&4812 * ( v_comp (k) - 2.0_wp * hom(k,1,2,0) )&4813 / ( v_comp (k) - gv + SIGN( 1.0E-20_wp, v_comp(k) - gv ) )&4814 + diss_n (k)&4815 * ABS( v_comp (k) - 2.0_wp * hom(k,1,2,0) )&4816 / ( ABS( v_comp (k) - gv ) + 1.0E-20_wp )&5437 + ( flux_n & 5438 * ( v_comp - 2.0_wp * hom(k,1,2,0) ) & 5439 / ( v_comp - gv + SIGN( 1.0E-20_wp, v_comp - gv ) ) & 5440 + diss_n & 5441 * ABS( v_comp - 2.0_wp * hom(k,1,2,0) ) & 5442 / ( ABS( v_comp - gv ) + 1.0E-20_wp ) & 4817 5443 ) * weight_substep(intermediate_timestep_count) 4818 5444 ! 4819 5445 !-- Statistical Evaluation of w'u'. 5446 !$ACC ATOMIC 4820 5447 sums_wsvs_ws_l(k,tn) = sums_wsvs_ws_l(k,tn) & 4821 + ( flux_t (k)&5448 + ( flux_t & 4822 5449 * ( w_comp - 2.0_wp * hom(k,1,3,0) ) & 4823 5450 / ( w_comp + SIGN( 1.0E-20_wp, w_comp ) ) & 4824 + diss_t (k)&5451 + diss_t & 4825 5452 * ABS( w_comp - 2.0_wp * hom(k,1,3,0) ) & 4826 5453 / ( ABS( w_comp ) + 1.0E-20_wp ) & … … 4830 5457 ENDDO 4831 5458 ENDDO 5459 !$ACC UPDATE HOST(sums_vs2_ws_l(nzb+1,tn)) 4832 5460 sums_vs2_ws_l(nzb,tn) = sums_vs2_ws_l(nzb+1,tn) 5461 !$ACC UPDATE DEVICE(sums_vs2_ws_l(nzb,tn)) 4833 5462 4834 5463 … … 4874 5503 REAL(wp) :: ibit28 !< flag indicating 3rd-order scheme along x-direction 4875 5504 REAL(wp) :: ibit29 !< flag indicating 5th-order scheme along x-direction 5505 #ifdef _OPENACC 5506 REAL(wp) :: ibit27_l !< flag indicating 1st-order scheme along x-direction 5507 REAL(wp) :: ibit28_l !< flag indicating 3rd-order scheme along x-direction 5508 REAL(wp) :: ibit29_l !< flag indicating 5th-order scheme along x-direction 5509 #endif 4876 5510 REAL(wp) :: ibit30 !< flag indicating 1st-order scheme along y-direction 4877 5511 REAL(wp) :: ibit31 !< flag indicating 3rd-order scheme along y-direction 4878 5512 REAL(wp) :: ibit32 !< flag indicating 5th-order scheme along y-direction 5513 #ifdef _OPENACC 5514 REAL(wp) :: ibit30_s !< flag indicating 1st-order scheme along y-direction 5515 REAL(wp) :: ibit31_s !< flag indicating 3rd-order scheme along y-direction 5516 REAL(wp) :: ibit32_s !< flag indicating 5th-order scheme along y-direction 5517 #endif 4879 5518 REAL(wp) :: ibit33 !< flag indicating 1st-order scheme along z-direction 4880 5519 REAL(wp) :: ibit34 !< flag indicating 3rd-order scheme along z-direction … … 4886 5525 REAL(wp) :: gv !< Galilei-transformation velocity along y 4887 5526 REAL(wp) :: u_comp !< advection velocity along x 5527 #ifdef _OPENACC 5528 REAL(wp) :: u_comp_l !< advection velocity along x 5529 #endif 4888 5530 REAL(wp) :: v_comp !< advection velocity along y 5531 #ifdef _OPENACC 5532 REAL(wp) :: v_comp_s !< advection velocity along y 5533 #endif 4889 5534 REAL(wp) :: w_comp !< advection velocity along z 4890 5535 4891 REAL(wp) , DIMENSION(nzb:nzt):: diss_t !< discretized artificial dissipation at top of the grid box4892 REAL(wp) , DIMENSION(nzb:nzt):: flux_t !< discretized 6th-order flux at top of the grid box5536 REAL(wp) :: diss_t !< discretized artificial dissipation at top of the grid box 5537 REAL(wp) :: flux_t !< discretized 6th-order flux at top of the grid box 4893 5538 4894 REAL(wp), DIMENSION(nzb+1:nzt) :: diss_n !< discretized artificial dissipation at northward-side of the grid box 4895 REAL(wp), DIMENSION(nzb+1:nzt) :: diss_r !< discretized artificial dissipation at rightward-side of the grid box 4896 REAL(wp), DIMENSION(nzb+1:nzt) :: flux_n !< discretized 6th-order flux at northward-side of the grid box 4897 REAL(wp), DIMENSION(nzb+1:nzt) :: flux_r !< discretized 6th-order flux at rightward-side of the grid box 5539 REAL(wp) :: diss_n !< discretized artificial dissipation at northward-side of the grid box 5540 REAL(wp) :: diss_r !< discretized artificial dissipation at rightward-side of the grid box 5541 REAL(wp) :: flux_n !< discretized 6th-order flux at northward-side of the grid box 5542 REAL(wp) :: flux_r !< discretized 6th-order flux at rightward-side of the grid box 5543 5544 REAL(wp) :: diss_s !< discretized artificial dissipation at southward-side of the grid box 5545 REAL(wp) :: flux_s !< discretized 6th-order flux at southward-side of the grid box 5546 #ifndef _OPENACC 4898 5547 REAL(wp), DIMENSION(nzb+1:nzt) :: swap_diss_y_local_w !< discretized artificial dissipation at southward-side of the grid box 4899 5548 REAL(wp), DIMENSION(nzb+1:nzt) :: swap_flux_y_local_w !< discretized 6th-order flux at southward-side of the grid box 5549 #endif 4900 5550 5551 REAL(wp) :: diss_l !< discretized artificial dissipation at leftward-side of the grid box 5552 REAL(wp) :: flux_l !< discretized 6th-order flux at leftward-side of the grid box 5553 #ifndef _OPENACC 4901 5554 REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) :: swap_diss_x_local_w !< discretized artificial dissipation at leftward-side of the grid box 4902 5555 REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) :: swap_flux_x_local_w !< discretized 6th-order flux at leftward-side of the grid box 5556 #endif 4903 5557 4904 5558 gu = 2.0_wp * u_gtrans 4905 5559 gv = 2.0_wp * v_gtrans 5560 5561 #ifndef _OPENACC 4906 5562 ! 4907 5563 !-- compute the whole left boundary of the processor domain … … 4962 5618 4963 5619 ENDDO 4964 5620 #endif 5621 5622 !$ACC PARALLEL LOOP COLLAPSE(2) FIRSTPRIVATE(tn, gu, gv) & 5623 !$ACC PRIVATE(i, j, k, k_mm, k_pp, k_ppp) & 5624 !$ACC PRIVATE(ibit27, ibit28, ibit29, ibit30, ibit31, ibit32) & 5625 !$ACC PRIVATE(ibit27_l, ibit28_l, ibit29_l) & 5626 !$ACC PRIVATE(ibit30_s, ibit31_s, ibit32_s) & 5627 !$ACC PRIVATE(ibit33, ibit34, ibit35) & 5628 !$ACC PRIVATE(flux_r, diss_r, flux_l, diss_l) & 5629 !$ACC PRIVATE(flux_n, diss_n, flux_s, diss_s) & 5630 !$ACC PRIVATE(flux_t, diss_t, flux_d, diss_d) & 5631 !$ACC PRIVATE(div, u_comp, u_comp_l, v_comp, v_comp_s, w_comp) & 5632 !$ACC PRESENT(advc_flags_1, advc_flags_2) & 5633 !$ACC PRESENT(u, v, w) & 5634 !$ACC PRESENT(rho_air, drho_air_zw, ddzu) & 5635 !$ACC PRESENT(tend) & 5636 !$ACC PRESENT(hom(nzb+1:nzb_max,1,3,0)) & 5637 !$ACC PRESENT(weight_substep(intermediate_timestep_count)) & 5638 !$ACC PRESENT(sums_ws2_ws_l(nzb+1:nzb_max,0)) 4965 5639 DO i = nxl, nxr 4966 5640 5641 #ifndef _OPENACC 4967 5642 j = nys 4968 5643 DO k = nzb+1, nzb_max … … 5018 5693 5019 5694 ENDDO 5695 #endif 5020 5696 5021 5697 DO j = nys, nyn … … 5025 5701 !-- at the first w-level. For topography wall this is done implicitely 5026 5702 !-- by advc_flags_1. 5027 k = nzb + 1 5028 w_comp = w(k,j,i) + w(k-1,j,i) 5029 flux_t(0) = w_comp * ( w(k,j,i) + w(k-1,j,i) ) * adv_mom_1 5030 diss_t(0) = -ABS(w_comp) * ( w(k,j,i) - w(k-1,j,i) ) * adv_mom_1 5031 flux_d = flux_t(0) 5032 diss_d = diss_t(0) 5703 k = nzb + 1 5704 w_comp = w(k,j,i) + w(k-1,j,i) 5705 flux_d = w_comp * ( w(k,j,i) + w(k-1,j,i) ) * adv_mom_1 5706 diss_d = -ABS(w_comp) * ( w(k,j,i) - w(k-1,j,i) ) * adv_mom_1 5033 5707 5034 5708 DO k = nzb+1, nzb_max … … 5038 5712 ibit27 = REAL( IBITS(advc_flags_1(k,j,i),27,1), KIND = wp ) 5039 5713 5040 u_comp 5041 flux_r (k) = u_comp * (&5714 u_comp = u(k+1,j,i+1) + u(k,j,i+1) - gu 5715 flux_r = u_comp * ( & 5042 5716 ( 37.0_wp * ibit29 * adv_mom_5 & 5043 5717 + 7.0_wp * ibit28 * adv_mom_3 & … … 5054 5728 ) 5055 5729 5056 diss_r (k) = - ABS( u_comp ) * (&5730 diss_r = - ABS( u_comp ) * ( & 5057 5731 ( 10.0_wp * ibit29 * adv_mom_5 & 5058 5732 + 3.0_wp * ibit28 * adv_mom_3 & … … 5069 5743 ) 5070 5744 5745 #ifdef _OPENACC 5746 ! 5747 !-- Recompute the left fluxes. 5748 ibit29_l = REAL( IBITS(advc_flags_1(k,j,i-1),29,1), KIND = wp ) 5749 ibit28_l = REAL( IBITS(advc_flags_1(k,j,i-1),28,1), KIND = wp ) 5750 ibit27_l = REAL( IBITS(advc_flags_1(k,j,i-1),27,1), KIND = wp ) 5751 5752 u_comp_l = u(k+1,j,i) + u(k,j,i) - gu 5753 flux_l = u_comp_l * ( & 5754 ( 37.0_wp * ibit29_l * adv_mom_5 & 5755 + 7.0_wp * ibit28_l * adv_mom_3 & 5756 + ibit27_l * adv_mom_1 & 5757 ) * & 5758 ( w(k,j,i) + w(k,j,i-1) ) & 5759 - ( 8.0_wp * ibit29_l * adv_mom_5 & 5760 + ibit28_l * adv_mom_3 & 5761 ) * & 5762 ( w(k,j,i+1) + w(k,j,i-2) ) & 5763 + ( ibit29_l * adv_mom_5 & 5764 ) * & 5765 ( w(k,j,i+2) + w(k,j,i-3) ) & 5766 ) 5767 5768 diss_l = - ABS( u_comp_l ) * ( & 5769 ( 10.0_wp * ibit29_l * adv_mom_5 & 5770 + 3.0_wp * ibit28_l * adv_mom_3 & 5771 + ibit27_l * adv_mom_1 & 5772 ) * & 5773 ( w(k,j,i) - w(k,j,i-1) ) & 5774 - ( 5.0_wp * ibit29_l * adv_mom_5 & 5775 + ibit28_l * adv_mom_3 & 5776 ) * & 5777 ( w(k,j,i+1) - w(k,j,i-2) ) & 5778 + ( ibit29_l * adv_mom_5 & 5779 ) * & 5780 ( w(k,j,i+2) - w(k,j,i-3) ) & 5781 ) 5782 #else 5783 flux_l = swap_flux_x_local_w(k,j) 5784 diss_l = swap_diss_x_local_w(k,j) 5785 #endif 5786 5787 5071 5788 ibit32 = REAL( IBITS(advc_flags_2(k,j,i),0,1), KIND = wp ) 5072 5789 ibit31 = REAL( IBITS(advc_flags_1(k,j,i),31,1), KIND = wp ) 5073 5790 ibit30 = REAL( IBITS(advc_flags_1(k,j,i),30,1), KIND = wp ) 5074 5791 5075 v_comp 5076 flux_n (k) = v_comp * (&5792 v_comp = v(k+1,j+1,i) + v(k,j+1,i) - gv 5793 flux_n = v_comp * ( & 5077 5794 ( 37.0_wp * ibit32 * adv_mom_5 & 5078 5795 + 7.0_wp * ibit31 * adv_mom_3 & … … 5089 5806 ) 5090 5807 5091 diss_n (k) = - ABS( v_comp ) * (&5808 diss_n = - ABS( v_comp ) * ( & 5092 5809 ( 10.0_wp * ibit32 * adv_mom_5 & 5093 5810 + 3.0_wp * ibit31 * adv_mom_3 & … … 5103 5820 ( w(k,j+3,i) - w(k,j-2,i) ) & 5104 5821 ) 5822 5823 #ifdef _OPENACC 5824 ! 5825 !-- Recompute the south fluxes. 5826 ibit32_s = REAL( IBITS(advc_flags_2(k,j-1,i),0,1), KIND = wp ) 5827 ibit31_s = REAL( IBITS(advc_flags_1(k,j-1,i),31,1), KIND = wp ) 5828 ibit30_s = REAL( IBITS(advc_flags_1(k,j-1,i),30,1), KIND = wp ) 5829 5830 v_comp_s = v(k+1,j,i) + v(k,j,i) - gv 5831 flux_s = v_comp_s * ( & 5832 ( 37.0_wp * ibit32_s * adv_mom_5 & 5833 + 7.0_wp * ibit31_s * adv_mom_3 & 5834 + ibit30_s * adv_mom_1 & 5835 ) * & 5836 ( w(k,j,i) + w(k,j-1,i) ) & 5837 - ( 8.0_wp * ibit32_s * adv_mom_5 & 5838 + ibit31_s * adv_mom_3 & 5839 ) * & 5840 ( w(k,j+1,i) + w(k,j-2,i) ) & 5841 + ( ibit32_s * adv_mom_5 & 5842 ) * & 5843 ( w(k,j+2,i) + w(k,j-3,i) ) & 5844 ) 5845 5846 diss_s = - ABS( v_comp_s ) * ( & 5847 ( 10.0_wp * ibit32_s * adv_mom_5 & 5848 + 3.0_wp * ibit31_s * adv_mom_3 & 5849 + ibit30_s * adv_mom_1 & 5850 ) * & 5851 ( w(k,j,i) - w(k,j-1,i) ) & 5852 - ( 5.0_wp * ibit32_s * adv_mom_5 & 5853 + ibit31_s * adv_mom_3 & 5854 ) * & 5855 ( w(k,j+1,i) - w(k,j-2,i) ) & 5856 + ( ibit32_s * adv_mom_5 & 5857 ) * & 5858 ( w(k,j+2,i) - w(k,j-3,i) ) & 5859 ) 5860 #else 5861 flux_s = swap_flux_y_local_w(k) 5862 diss_s = swap_diss_y_local_w(k) 5863 #endif 5864 5105 5865 ! 5106 5866 !-- k index has to be modified near bottom and top, else array … … 5114 5874 k_mm = k - 2 * ibit35 5115 5875 5116 w_comp 5117 flux_t (k) = w_comp * rho_air(k+1) * (&5876 w_comp = w(k+1,j,i) + w(k,j,i) 5877 flux_t = w_comp * rho_air(k+1) * ( & 5118 5878 ( 37.0_wp * ibit35 * adv_mom_5 & 5119 5879 + 7.0_wp * ibit34 * adv_mom_3 & … … 5130 5890 ) 5131 5891 5132 diss_t (k) = - ABS( w_comp ) * rho_air(k+1) * (&5892 diss_t = - ABS( w_comp ) * rho_air(k+1) * ( & 5133 5893 ( 10.0_wp * ibit35 * adv_mom_5 & 5134 5894 + 3.0_wp * ibit34 * adv_mom_3 & … … 5177 5937 5178 5938 tend(k,j,i) = tend(k,j,i) - ( & 5179 ( flux_r(k) + diss_r(k) & 5180 - swap_flux_x_local_w(k,j) - swap_diss_x_local_w(k,j) & 5181 ) * ddx & 5182 + ( flux_n(k) + diss_n(k) & 5183 - swap_flux_y_local_w(k) - swap_diss_y_local_w(k) & 5184 ) * ddy & 5185 + ( ( flux_t(k) + diss_t(k) ) & 5186 - ( flux_d + diss_d ) & 5187 ) * drho_air_zw(k) * ddzu(k+1) & 5939 ( ( flux_r + diss_r ) & 5940 - ( flux_l + diss_l ) ) * ddx & 5941 + ( ( flux_n + diss_n ) & 5942 - ( flux_s + diss_s ) ) * ddy & 5943 + ( ( flux_t + diss_t ) & 5944 - ( flux_d + diss_d ) ) * drho_air_zw(k) * ddzu(k+1) & 5188 5945 ) + div * w(k,j,i) 5189 5946 5190 swap_flux_x_local_w(k,j) = flux_r(k) 5191 swap_diss_x_local_w(k,j) = diss_r(k) 5192 swap_flux_y_local_w(k) = flux_n(k) 5193 swap_diss_y_local_w(k) = diss_n(k) 5194 flux_d = flux_t(k) 5195 diss_d = diss_t(k) 5196 5947 #ifndef _OPENACC 5948 swap_flux_x_local_w(k,j) = flux_r 5949 swap_diss_x_local_w(k,j) = diss_r 5950 swap_flux_y_local_w(k) = flux_n 5951 swap_diss_y_local_w(k) = diss_n 5952 #endif 5953 flux_d = flux_t 5954 diss_d = diss_t 5955 5956 !$ACC ATOMIC 5197 5957 sums_ws2_ws_l(k,tn) = sums_ws2_ws_l(k,tn) & 5198 + ( flux_t (k)&5958 + ( flux_t & 5199 5959 * ( w_comp - 2.0_wp * hom(k,1,3,0) ) & 5200 5960 / ( w_comp + SIGN( 1.0E-20_wp, w_comp ) ) & 5201 + diss_t (k)&5961 + diss_t & 5202 5962 * ABS( w_comp - 2.0_wp * hom(k,1,3,0) ) & 5203 5963 / ( ABS( w_comp ) + 1.0E-20_wp ) & … … 5208 5968 DO k = nzb_max+1, nzt 5209 5969 5210 u_comp 5211 flux_r (k) = u_comp * (&5970 u_comp = u(k+1,j,i+1) + u(k,j,i+1) - gu 5971 flux_r = u_comp * ( & 5212 5972 37.0_wp * ( w(k,j,i+1) + w(k,j,i) ) & 5213 5973 - 8.0_wp * ( w(k,j,i+2) + w(k,j,i-1) ) & 5214 5974 + ( w(k,j,i+3) + w(k,j,i-2) ) ) * adv_mom_5 5215 5975 5216 diss_r (k) = - ABS( u_comp ) * (&5976 diss_r = - ABS( u_comp ) * ( & 5217 5977 10.0_wp * ( w(k,j,i+1) - w(k,j,i) ) & 5218 5978 - 5.0_wp * ( w(k,j,i+2) - w(k,j,i-1) ) & 5219 5979 + ( w(k,j,i+3) - w(k,j,i-2) ) ) * adv_mom_5 5220 5980 5221 v_comp = v(k+1,j+1,i) + v(k,j+1,i) - gv 5222 flux_n(k) = v_comp * ( & 5981 #ifdef _OPENACC 5982 ! 5983 !-- Recompute the left fluxes. 5984 u_comp_l = u(k+1,j,i) + u(k,j,i) - gu 5985 flux_l = u_comp_l * ( & 5986 37.0_wp * ( w(k,j,i) + w(k,j,i-1) ) & 5987 - 8.0_wp * ( w(k,j,i+1) + w(k,j,i-2) ) & 5988 + ( w(k,j,i+2) + w(k,j,i-3) ) ) * adv_mom_5 5989 diss_l = - ABS( u_comp_l ) * ( & 5990 10.0_wp * ( w(k,j,i) - w(k,j,i-1) ) & 5991 - 5.0_wp * ( w(k,j,i+1) - w(k,j,i-2) ) & 5992 + ( w(k,j,i+2) - w(k,j,i-3) ) ) * adv_mom_5 5993 #else 5994 flux_l = swap_flux_x_local_w(k,j) 5995 diss_l = swap_diss_x_local_w(k,j) 5996 #endif 5997 5998 v_comp = v(k+1,j+1,i) + v(k,j+1,i) - gv 5999 flux_n = v_comp * ( & 5223 6000 37.0_wp * ( w(k,j+1,i) + w(k,j,i) ) & 5224 6001 - 8.0_wp * ( w(k,j+2,i) + w(k,j-1,i) ) & 5225 6002 + ( w(k,j+3,i) + w(k,j-2,i) ) ) * adv_mom_5 5226 6003 5227 diss_n (k) = - ABS( v_comp ) * (&6004 diss_n = - ABS( v_comp ) * ( & 5228 6005 10.0_wp * ( w(k,j+1,i) - w(k,j,i) ) & 5229 6006 - 5.0_wp * ( w(k,j+2,i) - w(k,j-1,i) ) & 5230 6007 + ( w(k,j+3,i) - w(k,j-2,i) ) ) * adv_mom_5 6008 6009 #ifdef _OPENACC 6010 ! 6011 !-- Recompute the south fluxes. 6012 v_comp_s = v(k+1,j,i) + v(k,j,i) - gv 6013 flux_s = v_comp_s * ( & 6014 37.0_wp * ( w(k,j,i) + w(k,j-1,i) ) & 6015 - 8.0_wp * ( w(k,j+1,i) +w(k,j-2,i) ) & 6016 + ( w(k,j+2,i) + w(k,j-3,i) ) ) * adv_mom_5 6017 diss_s = - ABS( v_comp_s ) * ( & 6018 10.0_wp * ( w(k,j,i) - w(k,j-1,i) ) & 6019 - 5.0_wp * ( w(k,j+1,i) - w(k,j-2,i) ) & 6020 + ( w(k,j+2,i) - w(k,j-3,i) ) ) * adv_mom_5 6021 #else 6022 flux_s = swap_flux_y_local_w(k) 6023 diss_s = swap_diss_y_local_w(k) 6024 #endif 6025 5231 6026 ! 5232 6027 !-- k index has to be modified near bottom and top, else array … … 5240 6035 k_mm = k - 2 * ibit35 5241 6036 5242 w_comp 5243 flux_t (k) = w_comp * rho_air(k+1) * (&6037 w_comp = w(k+1,j,i) + w(k,j,i) 6038 flux_t = w_comp * rho_air(k+1) * ( & 5244 6039 ( 37.0_wp * ibit35 * adv_mom_5 & 5245 6040 + 7.0_wp * ibit34 * adv_mom_3 & … … 5256 6051 ) 5257 6052 5258 diss_t (k) = - ABS( w_comp ) * rho_air(k+1) * (&6053 diss_t = - ABS( w_comp ) * rho_air(k+1) * ( & 5259 6054 ( 10.0_wp * ibit35 * adv_mom_5 & 5260 6055 + 3.0_wp * ibit34 * adv_mom_3 & … … 5282 6077 5283 6078 tend(k,j,i) = tend(k,j,i) - ( & 5284 ( flux_r(k) + diss_r(k) & 5285 - swap_flux_x_local_w(k,j) - swap_diss_x_local_w(k,j) & 5286 ) * ddx & 5287 + ( flux_n(k) + diss_n(k) & 5288 - swap_flux_y_local_w(k) - swap_diss_y_local_w(k) & 5289 ) * ddy & 5290 + ( ( flux_t(k) + diss_t(k) ) & 5291 - ( flux_d + diss_d ) & 5292 ) * drho_air_zw(k) * ddzu(k+1) & 6079 ( ( flux_r + diss_r ) & 6080 - ( flux_l + diss_l ) ) * ddx & 6081 + ( ( flux_n + diss_n ) & 6082 - ( flux_s + diss_s ) ) * ddy & 6083 + ( ( flux_t + diss_t ) & 6084 - ( flux_d + diss_d ) ) * drho_air_zw(k) * ddzu(k+1) & 5293 6085 ) + div * w(k,j,i) 5294 6086 5295 swap_flux_x_local_w(k,j) = flux_r(k) 5296 swap_diss_x_local_w(k,j) = diss_r(k) 5297 swap_flux_y_local_w(k) = flux_n(k) 5298 swap_diss_y_local_w(k) = diss_n(k) 5299 flux_d = flux_t(k) 5300 diss_d = diss_t(k) 5301 6087 #ifndef _OPENACC 6088 swap_flux_x_local_w(k,j) = flux_r 6089 swap_diss_x_local_w(k,j) = diss_r 6090 swap_flux_y_local_w(k) = flux_n 6091 swap_diss_y_local_w(k) = diss_n 6092 #endif 6093 flux_d = flux_t 6094 diss_d = diss_t 6095 6096 !$ACC ATOMIC 5302 6097 sums_ws2_ws_l(k,tn) = sums_ws2_ws_l(k,tn) & 5303 + ( flux_t (k)&6098 + ( flux_t & 5304 6099 * ( w_comp - 2.0_wp * hom(k,1,3,0) ) & 5305 6100 / ( w_comp + SIGN( 1.0E-20_wp, w_comp ) ) & 5306 + diss_t (k)&6101 + diss_t & 5307 6102 * ABS( w_comp - 2.0_wp * hom(k,1,3,0) ) & 5308 6103 / ( ABS( w_comp ) + 1.0E-20_wp ) &
Note: See TracChangeset
for help on using the changeset viewer.
