MODULE advec_ws !--------------------------------------------------------------------------------! ! This file is part of PALM. ! ! PALM is free software: you can redistribute it and/or modify it under the terms ! of the GNU General Public License as published by the Free Software Foundation, ! either version 3 of the License, or (at your option) any later version. ! ! PALM is distributed in the hope that it will be useful, but WITHOUT ANY ! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR ! A PARTICULAR PURPOSE. See the GNU General Public License for more details. ! ! You should have received a copy of the GNU General Public License along with ! PALM. If not, see . ! ! Copyright 1997-2014 Leibniz Universitaet Hannover !--------------------------------------------------------------------------------! ! ! Current revisions: ! ------------------ ! ! ! Former revisions: ! ----------------- ! $Id: advec_ws.f90 1375 2014-04-25 13:07:08Z fricke $ ! ! 1374 2014-04-25 12:55:07Z raasch ! missing variables added to ONLY list ! ! 1361 2014-04-16 15:17:48Z hoffmann ! accelerator and vector version for qr and nr added ! ! 1353 2014-04-08 15:21:23Z heinze ! REAL constants provided with KIND-attribute, ! module kinds added ! some formatting adjustments ! ! 1322 2014-03-20 16:38:49Z raasch ! REAL constants defined as wp-kind ! ! 1320 2014-03-20 08:40:49Z raasch ! ONLY-attribute added to USE-statements, ! kind-parameters added to all INTEGER and REAL declaration statements, ! kinds are defined in new module kinds, ! old module precision_kind is removed, ! revision history before 2012 removed, ! comment fields (!:) to be used for variable explanations added to ! all variable declaration statements ! ! 1257 2013-11-08 15:18:40Z raasch ! accelerator loop directives removed ! ! 1221 2013-09-10 08:59:13Z raasch ! wall_flags_00 introduced, which holds bits 32-... ! ! 1128 2013-04-12 06:19:32Z raasch ! loop index bounds in accelerator version replaced by i_left, i_right, j_south, ! j_north ! ! 1115 2013-03-26 18:16:16Z hoffmann ! calculation of qr and nr is restricted to precipitation ! ! 1053 2012-11-13 17:11:03Z hoffmann ! necessary expansions according to the two new prognostic equations (nr, qr) ! of the two-moment cloud physics scheme: ! +flux_l_*, flux_s_*, diss_l_*, diss_s_*, sums_ws*s_ws_l ! ! 1036 2012-10-22 13:43:42Z raasch ! code put under GPL (PALM 3.9) ! ! 1027 2012-10-15 17:18:39Z suehring ! Bugfix in calculation indices k_mm, k_pp in accelerator version ! ! 1019 2012-09-28 06:46:45Z raasch ! small change in comment lines ! ! 1015 2012-09-27 09:23:24Z raasch ! accelerator versions (*_acc) added ! ! 1010 2012-09-20 07:59:54Z raasch ! cpp switch __nopointer added for pointer free version ! ! 888 2012-04-20 15:03:46Z suehring ! Number of IBITS() calls with identical arguments is reduced. ! ! 862 2012-03-26 14:21:38Z suehring ! ws-scheme also work with topography in combination with vector version. ! ws-scheme also work with outflow boundaries in combination with ! vector version. ! Degradation of the applied order of scheme is now steered by multiplying with ! Integer wall_flags_0. 2nd order scheme, WS3 and WS5 are calculated on each ! grid point and mulitplied with the appropriate flag. ! 2nd order numerical dissipation term changed. Now the appropriate 2nd order ! term derived according to the 4th and 6th order terms is applied. It turns ! out that diss_2nd does not provide sufficient dissipation near walls. ! Therefore, the function diss_2nd is removed. ! Near walls a divergence correction is necessary to overcome numerical ! instabilities due to too less divergence reduction of the velocity field. ! boundary_flags and logicals steering the degradation are removed. ! Empty SUBROUTINE local_diss removed. ! Further formatting adjustments. ! ! 801 2012-01-10 17:30:36Z suehring ! Bugfix concerning OpenMP parallelization. Summation of sums_wsus_ws_l, ! sums_wsvs_ws_l, sums_us2_ws_l, sums_vs2_ws_l, sums_ws2_ws_l, sums_wspts_ws_l, ! sums_wsqs_ws_l, sums_wssas_ws_l is now thread-safe by adding an additional ! dimension. ! ! Initial revision ! ! 411 2009-12-11 12:31:43 Z suehring ! ! Description: ! ------------ ! Advection scheme for scalars and momentum using the flux formulation of ! Wicker and Skamarock 5th order. Additionally the module contains of a ! routine using for initialisation and steering of the statical evaluation. ! The computation of turbulent fluxes takes place inside the advection ! routines. ! Near non-cyclic boundaries the order of the applied advection scheme is ! degraded. ! A divergence correction is applied. It is necessary for topography, since ! the divergence is not sufficiently reduced, resulting in erroneous fluxes and ! partly numerical instabilities. !-----------------------------------------------------------------------------! PRIVATE PUBLIC advec_s_ws, advec_s_ws_acc, advec_u_ws, advec_u_ws_acc, & advec_v_ws, advec_v_ws_acc, advec_w_ws, advec_w_ws_acc, & ws_init, ws_statistics INTERFACE ws_init MODULE PROCEDURE ws_init END INTERFACE ws_init INTERFACE ws_statistics MODULE PROCEDURE ws_statistics END INTERFACE ws_statistics INTERFACE advec_s_ws MODULE PROCEDURE advec_s_ws MODULE PROCEDURE advec_s_ws_ij END INTERFACE advec_s_ws INTERFACE advec_u_ws MODULE PROCEDURE advec_u_ws MODULE PROCEDURE advec_u_ws_ij END INTERFACE advec_u_ws INTERFACE advec_u_ws_acc MODULE PROCEDURE advec_u_ws_acc END INTERFACE advec_u_ws_acc INTERFACE advec_v_ws MODULE PROCEDURE advec_v_ws MODULE PROCEDURE advec_v_ws_ij END INTERFACE advec_v_ws INTERFACE advec_v_ws_acc MODULE PROCEDURE advec_v_ws_acc END INTERFACE advec_v_ws_acc INTERFACE advec_w_ws MODULE PROCEDURE advec_w_ws MODULE PROCEDURE advec_w_ws_ij END INTERFACE advec_w_ws INTERFACE advec_w_ws_acc MODULE PROCEDURE advec_w_ws_acc END INTERFACE advec_w_ws_acc CONTAINS !------------------------------------------------------------------------------! ! Initialization of WS-scheme !------------------------------------------------------------------------------! SUBROUTINE ws_init USE arrays_3d, & ONLY: diss_l_e, diss_l_nr, diss_l_pt, diss_l_q, diss_l_qr, & diss_l_sa, diss_l_u, diss_l_v, diss_l_w, flux_l_e, & flux_l_nr, flux_l_pt, flux_l_q, flux_l_qr, flux_l_sa, & flux_l_u, flux_l_v, flux_l_w, diss_s_e, diss_s_nr, diss_s_pt,& diss_s_q, diss_s_qr, diss_s_sa, diss_s_u, diss_s_v, diss_s_w,& flux_s_e, flux_s_nr, flux_s_pt, flux_s_q, flux_s_qr, & flux_s_sa, flux_s_u, flux_s_v, flux_s_w USE constants, & ONLY: adv_mom_1, adv_mom_3, adv_mom_5, adv_sca_1, adv_sca_3, & adv_sca_5 USE control_parameters, & ONLY: cloud_physics, humidity, icloud_scheme, loop_optimization, & passive_scalar, precipitation, ocean, ws_scheme_mom, & ws_scheme_sca USE indices, & ONLY: nyn, nys, nzb, nzt USE kinds USE pegrid USE statistics, & ONLY: sums_us2_ws_l, sums_vs2_ws_l, sums_ws2_ws_l, sums_wsnrs_ws_l,& sums_wspts_ws_l, sums_wsqrs_ws_l, sums_wsqs_ws_l, & sums_wssas_ws_l, sums_wsus_ws_l, sums_wsvs_ws_l ! !-- Set the appropriate factors for scalar and momentum advection. adv_sca_5 = 1.0_wp / 60.0_wp adv_sca_3 = 1.0_wp / 12.0_wp adv_sca_1 = 1.0_wp / 2.0_wp adv_mom_5 = 1.0_wp / 120.0_wp adv_mom_3 = 1.0_wp / 24.0_wp adv_mom_1 = 1.0_wp / 4.0_wp ! !-- Arrays needed for statical evaluation of fluxes. IF ( ws_scheme_mom ) THEN ALLOCATE( sums_wsus_ws_l(nzb:nzt+1,0:threads_per_task-1), & sums_wsvs_ws_l(nzb:nzt+1,0:threads_per_task-1), & sums_us2_ws_l(nzb:nzt+1,0:threads_per_task-1), & sums_vs2_ws_l(nzb:nzt+1,0:threads_per_task-1), & sums_ws2_ws_l(nzb:nzt+1,0:threads_per_task-1) ) sums_wsus_ws_l = 0.0_wp sums_wsvs_ws_l = 0.0_wp sums_us2_ws_l = 0.0_wp sums_vs2_ws_l = 0.0_wp sums_ws2_ws_l = 0.0_wp ENDIF IF ( ws_scheme_sca ) THEN ALLOCATE( sums_wspts_ws_l(nzb:nzt+1,0:threads_per_task-1) ) sums_wspts_ws_l = 0.0_wp IF ( humidity .OR. passive_scalar ) THEN ALLOCATE( sums_wsqs_ws_l(nzb:nzt+1,0:threads_per_task-1) ) sums_wsqs_ws_l = 0.0_wp ENDIF IF ( cloud_physics .AND. icloud_scheme == 0 .AND. & precipitation ) THEN ALLOCATE( sums_wsqrs_ws_l(nzb:nzt+1,0:threads_per_task-1) ) ALLOCATE( sums_wsnrs_ws_l(nzb:nzt+1,0:threads_per_task-1) ) sums_wsqrs_ws_l = 0.0_wp sums_wsnrs_ws_l = 0.0_wp ENDIF IF ( ocean ) THEN ALLOCATE( sums_wssas_ws_l(nzb:nzt+1,0:threads_per_task-1) ) sums_wssas_ws_l = 0.0_wp ENDIF ENDIF ! !-- Arrays needed for reasons of speed optimization for cache version. !-- For the vector version the buffer arrays are not necessary, !-- because the the fluxes can swapped directly inside the loops of the !-- advection routines. IF ( loop_optimization /= 'vector' ) THEN IF ( ws_scheme_mom ) THEN ALLOCATE( flux_s_u(nzb+1:nzt,0:threads_per_task-1), & flux_s_v(nzb+1:nzt,0:threads_per_task-1), & flux_s_w(nzb+1:nzt,0:threads_per_task-1), & diss_s_u(nzb+1:nzt,0:threads_per_task-1), & diss_s_v(nzb+1:nzt,0:threads_per_task-1), & diss_s_w(nzb+1:nzt,0:threads_per_task-1) ) ALLOCATE( flux_l_u(nzb+1:nzt,nys:nyn,0:threads_per_task-1), & flux_l_v(nzb+1:nzt,nys:nyn,0:threads_per_task-1), & flux_l_w(nzb+1:nzt,nys:nyn,0:threads_per_task-1), & diss_l_u(nzb+1:nzt,nys:nyn,0:threads_per_task-1), & diss_l_v(nzb+1:nzt,nys:nyn,0:threads_per_task-1), & diss_l_w(nzb+1:nzt,nys:nyn,0:threads_per_task-1) ) ENDIF IF ( ws_scheme_sca ) THEN ALLOCATE( flux_s_pt(nzb+1:nzt,0:threads_per_task-1), & flux_s_e(nzb+1:nzt,0:threads_per_task-1), & diss_s_pt(nzb+1:nzt,0:threads_per_task-1), & diss_s_e(nzb+1:nzt,0:threads_per_task-1) ) ALLOCATE( flux_l_pt(nzb+1:nzt,nys:nyn,0:threads_per_task-1), & flux_l_e(nzb+1:nzt,nys:nyn,0:threads_per_task-1), & diss_l_pt(nzb+1:nzt,nys:nyn,0:threads_per_task-1), & diss_l_e(nzb+1:nzt,nys:nyn,0:threads_per_task-1) ) IF ( humidity .OR. passive_scalar ) THEN ALLOCATE( flux_s_q(nzb+1:nzt,0:threads_per_task-1), & diss_s_q(nzb+1:nzt,0:threads_per_task-1) ) ALLOCATE( flux_l_q(nzb+1:nzt,nys:nyn,0:threads_per_task-1), & diss_l_q(nzb+1:nzt,nys:nyn,0:threads_per_task-1) ) ENDIF IF ( cloud_physics .AND. icloud_scheme == 0 .AND. & precipitation ) THEN ALLOCATE( flux_s_qr(nzb+1:nzt,0:threads_per_task-1), & diss_s_qr(nzb+1:nzt,0:threads_per_task-1), & flux_s_nr(nzb+1:nzt,0:threads_per_task-1), & diss_s_nr(nzb+1:nzt,0:threads_per_task-1) ) ALLOCATE( flux_l_qr(nzb+1:nzt,nys:nyn,0:threads_per_task-1), & diss_l_qr(nzb+1:nzt,nys:nyn,0:threads_per_task-1), & flux_l_nr(nzb+1:nzt,nys:nyn,0:threads_per_task-1), & diss_l_nr(nzb+1:nzt,nys:nyn,0:threads_per_task-1) ) ENDIF IF ( ocean ) THEN ALLOCATE( flux_s_sa(nzb+1:nzt,0:threads_per_task-1), & diss_s_sa(nzb+1:nzt,0:threads_per_task-1) ) ALLOCATE( flux_l_sa(nzb+1:nzt,nys:nyn,0:threads_per_task-1), & diss_l_sa(nzb+1:nzt,nys:nyn,0:threads_per_task-1) ) ENDIF ENDIF ENDIF END SUBROUTINE ws_init !------------------------------------------------------------------------------! ! Initialize variables used for storing statistic quantities (fluxes, variances) !------------------------------------------------------------------------------! SUBROUTINE ws_statistics USE control_parameters, & ONLY: cloud_physics, humidity, icloud_scheme, passive_scalar, & precipitation, ocean, ws_scheme_mom, ws_scheme_sca USE kinds USE statistics, & ONLY: sums_us2_ws_l, sums_vs2_ws_l, sums_ws2_ws_l, sums_wsnrs_ws_l,& sums_wspts_ws_l, sums_wsqrs_ws_l, sums_wsqs_ws_l, & sums_wssas_ws_l, sums_wsus_ws_l, sums_wsvs_ws_l IMPLICIT NONE ! !-- The arrays needed for statistical evaluation are set to to 0 at the !-- beginning of prognostic_equations. IF ( ws_scheme_mom ) THEN sums_wsus_ws_l = 0.0_wp sums_wsvs_ws_l = 0.0_wp sums_us2_ws_l = 0.0_wp sums_vs2_ws_l = 0.0_wp sums_ws2_ws_l = 0.0_wp ENDIF IF ( ws_scheme_sca ) THEN sums_wspts_ws_l = 0.0_wp IF ( humidity .OR. passive_scalar ) sums_wsqs_ws_l = 0.0_wp IF ( cloud_physics .AND. icloud_scheme == 0 .AND. & precipitation ) THEN sums_wsqrs_ws_l = 0.0_wp sums_wsnrs_ws_l = 0.0_wp ENDIF IF ( ocean ) sums_wssas_ws_l = 0.0_wp ENDIF END SUBROUTINE ws_statistics !------------------------------------------------------------------------------! ! Scalar advection - Call for grid point i,j !------------------------------------------------------------------------------! SUBROUTINE advec_s_ws_ij( i, j, sk, sk_char, swap_flux_y_local, & swap_diss_y_local, swap_flux_x_local, & swap_diss_x_local, i_omp, tn ) USE arrays_3d, & ONLY: ddzw, tend, u, v, w USE constants, & ONLY: adv_sca_1, adv_sca_3, adv_sca_5 USE control_parameters, & ONLY: intermediate_timestep_count, u_gtrans, v_gtrans USE grid_variables, & ONLY: ddx, ddy USE indices, & ONLY: nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, nzb_max, & nzt, wall_flags_0 USE kinds USE pegrid USE statistics, & ONLY: sums_wsnrs_ws_l, sums_wspts_ws_l, sums_wsqrs_ws_l, & sums_wsqs_ws_l, sums_wssas_ws_l, weight_substep IMPLICIT NONE CHARACTER (LEN = *), INTENT(IN) :: sk_char !: INTEGER(iwp) :: i !: INTEGER(iwp) :: ibit0 !: INTEGER(iwp) :: ibit1 !: INTEGER(iwp) :: ibit2 !: INTEGER(iwp) :: ibit3 !: INTEGER(iwp) :: ibit4 !: INTEGER(iwp) :: ibit5 !: INTEGER(iwp) :: ibit6 !: INTEGER(iwp) :: ibit7 !: INTEGER(iwp) :: ibit8 !: INTEGER(iwp) :: i_omp !: INTEGER(iwp) :: j !: INTEGER(iwp) :: k !: INTEGER(iwp) :: k_mm !: INTEGER(iwp) :: k_pp !: INTEGER(iwp) :: k_ppp !: INTEGER(iwp) :: tn !: REAL(wp) :: diss_d !: REAL(wp) :: div !: REAL(wp) :: flux_d !: REAL(wp) :: u_comp !: REAL(wp) :: v_comp !: #if defined( __nopointer ) REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: sk !: #else REAL(wp), DIMENSION(:,:,:), POINTER :: sk !: #endif REAL(wp), DIMENSION(nzb:nzt+1) :: diss_n !: REAL(wp), DIMENSION(nzb:nzt+1) :: diss_r !: REAL(wp), DIMENSION(nzb:nzt+1) :: diss_t !: REAL(wp), DIMENSION(nzb:nzt+1) :: flux_n !: REAL(wp), DIMENSION(nzb:nzt+1) :: flux_r !: REAL(wp), DIMENSION(nzb:nzt+1) :: flux_t !: REAL(wp), DIMENSION(nzb+1:nzt,0:threads_per_task-1) :: swap_diss_y_local !: REAL(wp), DIMENSION(nzb+1:nzt,0:threads_per_task-1) :: swap_flux_y_local !: REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn,0:threads_per_task-1) :: swap_diss_x_local !: REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn,0:threads_per_task-1) :: swap_flux_x_local !: ! !-- Compute southside fluxes of the respective PE bounds. IF ( j == nys ) THEN ! !-- Up to the top of the highest topography. DO k = nzb+1, nzb_max ibit5 = IBITS(wall_flags_0(k,j,i),5,1) ibit4 = IBITS(wall_flags_0(k,j,i),4,1) ibit3 = IBITS(wall_flags_0(k,j,i),3,1) v_comp = v(k,j,i) - v_gtrans swap_flux_y_local(k,tn) = v_comp * ( & ( 37.0_wp * ibit5 * adv_sca_5 & + 7.0_wp * ibit4 * adv_sca_3 & + ibit3 * adv_sca_1 & ) * & ( sk(k,j,i) + sk(k,j-1,i) ) & - ( 8.0_wp * ibit5 * adv_sca_5 & + ibit4 * adv_sca_3 & ) * & ( sk(k,j+1,i) + sk(k,j-2,i) ) & + ( ibit5 * adv_sca_5 & ) * & ( sk(k,j+2,i) + sk(k,j-3,i) ) & ) swap_diss_y_local(k,tn) = -ABS( v_comp ) * ( & ( 10.0_wp * ibit5 * adv_sca_5 & + 3.0_wp * ibit4 * adv_sca_3 & + ibit3 * adv_sca_1 & ) * & ( sk(k,j,i) - sk(k,j-1,i) ) & - ( 5.0_wp * ibit5 * adv_sca_5 & + ibit4 * adv_sca_3 & ) * & ( sk(k,j+1,i) - sk(k,j-2,i) ) & + ( ibit5 * adv_sca_5 & ) * & ( sk(k,j+2,i) - sk(k,j-3,i) ) & ) ENDDO ! !-- Above to the top of the highest topography. No degradation necessary. DO k = nzb_max+1, nzt v_comp = v(k,j,i) - v_gtrans swap_flux_y_local(k,tn) = v_comp * ( & 37.0_wp * ( sk(k,j,i) + sk(k,j-1,i) ) & - 8.0_wp * ( sk(k,j+1,i) + sk(k,j-2,i) ) & + ( sk(k,j+2,i) + sk(k,j-3,i) ) & ) * adv_sca_5 swap_diss_y_local(k,tn) = -ABS( v_comp ) * ( & 10.0_wp * ( sk(k,j,i) - sk(k,j-1,i) ) & - 5.0_wp * ( sk(k,j+1,i) - sk(k,j-2,i) ) & + sk(k,j+2,i) - sk(k,j-3,i) & ) * adv_sca_5 ENDDO ENDIF ! !-- Compute leftside fluxes of the respective PE bounds. IF ( i == i_omp ) THEN DO k = nzb+1, nzb_max ibit2 = IBITS(wall_flags_0(k,j,i),2,1) ibit1 = IBITS(wall_flags_0(k,j,i),1,1) ibit0 = IBITS(wall_flags_0(k,j,i),0,1) u_comp = u(k,j,i) - u_gtrans swap_flux_x_local(k,j,tn) = u_comp * ( & ( 37.0_wp * ibit2 * adv_sca_5 & + 7.0_wp * ibit1 * adv_sca_3 & + ibit0 * adv_sca_1 & ) * & ( sk(k,j,i) + sk(k,j,i-1) ) & - ( 8.0_wp * ibit2 * adv_sca_5 & + ibit1 * adv_sca_3 & ) * & ( sk(k,j,i+1) + sk(k,j,i-2) ) & + ( ibit2 * adv_sca_5 & ) * & ( sk(k,j,i+2) + sk(k,j,i-3) ) & ) swap_diss_x_local(k,j,tn) = -ABS( u_comp ) * ( & ( 10.0_wp * ibit2 * adv_sca_5 & + 3.0_wp * ibit1 * adv_sca_3 & + ibit0 * adv_sca_1 & ) * & ( sk(k,j,i) - sk(k,j,i-1) ) & - ( 5.0_wp * ibit2 * adv_sca_5 & + ibit1 * adv_sca_3 & ) * & ( sk(k,j,i+1) - sk(k,j,i-2) ) & + ( ibit2 * adv_sca_5 & ) * & ( sk(k,j,i+2) - sk(k,j,i-3) ) & ) ENDDO DO k = nzb_max+1, nzt u_comp = u(k,j,i) - u_gtrans swap_flux_x_local(k,j,tn) = u_comp * ( & 37.0_wp * ( sk(k,j,i) + sk(k,j,i-1) ) & - 8.0_wp * ( sk(k,j,i+1) + sk(k,j,i-2) ) & + ( sk(k,j,i+2) + sk(k,j,i-3) ) & ) * adv_sca_5 swap_diss_x_local(k,j,tn) = -ABS( u_comp ) * ( & 10.0_wp * ( sk(k,j,i) - sk(k,j,i-1) ) & - 5.0_wp * ( sk(k,j,i+1) - sk(k,j,i-2) ) & + ( sk(k,j,i+2) - sk(k,j,i-3) ) & ) * adv_sca_5 ENDDO ENDIF flux_t(0) = 0.0_wp diss_t(0) = 0.0_wp flux_d = 0.0_wp diss_d = 0.0_wp ! !-- Now compute the fluxes and tendency terms for the horizontal and !-- vertical parts up to the top of the highest topography. DO k = nzb+1, nzb_max ! !-- Note: It is faster to conduct all multiplications explicitly, e.g. !-- * adv_sca_5 ... than to determine a factor and multiplicate the !-- flux at the end. ibit2 = IBITS(wall_flags_0(k,j,i),2,1) ibit1 = IBITS(wall_flags_0(k,j,i),1,1) ibit0 = IBITS(wall_flags_0(k,j,i),0,1) u_comp = u(k,j,i+1) - u_gtrans flux_r(k) = u_comp * ( & ( 37.0_wp * ibit2 * adv_sca_5 & + 7.0_wp * ibit1 * adv_sca_3 & + ibit0 * adv_sca_1 & ) * & ( sk(k,j,i+1) + sk(k,j,i) ) & - ( 8.0_wp * ibit2 * adv_sca_5 & + ibit1 * adv_sca_3 & ) * & ( sk(k,j,i+2) + sk(k,j,i-1) ) & + ( ibit2 * adv_sca_5 & ) * & ( sk(k,j,i+3) + sk(k,j,i-2) ) & ) diss_r(k) = -ABS( u_comp ) * ( & ( 10.0_wp * ibit2 * adv_sca_5 & + 3.0_wp * ibit1 * adv_sca_3 & + ibit0 * adv_sca_1 & ) * & ( sk(k,j,i+1) - sk(k,j,i) ) & - ( 5.0_wp * ibit2 * adv_sca_5 & + ibit1 * adv_sca_3 & ) * & ( sk(k,j,i+2) - sk(k,j,i-1) ) & + ( ibit2 * adv_sca_5 & ) * & ( sk(k,j,i+3) - sk(k,j,i-2) ) & ) ibit5 = IBITS(wall_flags_0(k,j,i),5,1) ibit4 = IBITS(wall_flags_0(k,j,i),4,1) ibit3 = IBITS(wall_flags_0(k,j,i),3,1) v_comp = v(k,j+1,i) - v_gtrans flux_n(k) = v_comp * ( & ( 37.0_wp * ibit5 * adv_sca_5 & + 7.0_wp * ibit4 * adv_sca_3 & + ibit3 * adv_sca_1 & ) * & ( sk(k,j+1,i) + sk(k,j,i) ) & - ( 8.0_wp * ibit5 * adv_sca_5 & + ibit4 * adv_sca_3 & ) * & ( sk(k,j+2,i) + sk(k,j-1,i) ) & + ( ibit5 * adv_sca_5 & ) * & ( sk(k,j+3,i) + sk(k,j-2,i) ) & ) diss_n(k) = -ABS( v_comp ) * ( & ( 10.0_wp * ibit5 * adv_sca_5 & + 3.0_wp * ibit4 * adv_sca_3 & + ibit3 * adv_sca_1 & ) * & ( sk(k,j+1,i) - sk(k,j,i) ) & - ( 5.0_wp * ibit5 * adv_sca_5 & + ibit4 * adv_sca_3 & ) * & ( sk(k,j+2,i) - sk(k,j-1,i) ) & + ( ibit5 * adv_sca_5 & ) * & ( sk(k,j+3,i) - sk(k,j-2,i) ) & ) ! !-- k index has to be modified near bottom and top, else array !-- subscripts will be exceeded. ibit8 = IBITS(wall_flags_0(k,j,i),8,1) ibit7 = IBITS(wall_flags_0(k,j,i),7,1) ibit6 = IBITS(wall_flags_0(k,j,i),6,1) k_ppp = k + 3 * ibit8 k_pp = k + 2 * ( 1 - ibit6 ) k_mm = k - 2 * ibit8 flux_t(k) = w(k,j,i) * ( & ( 37.0_wp * ibit8 * adv_sca_5 & + 7.0_wp * ibit7 * adv_sca_3 & + ibit6 * adv_sca_1 & ) * & ( sk(k+1,j,i) + sk(k,j,i) ) & - ( 8.0_wp * ibit8 * adv_sca_5 & + ibit7 * adv_sca_3 & ) * & ( sk(k_pp,j,i) + sk(k-1,j,i) ) & + ( ibit8 * adv_sca_5 & ) * ( sk(k_ppp,j,i)+ sk(k_mm,j,i) ) & ) diss_t(k) = -ABS( w(k,j,i) ) * ( & ( 10.0_wp * ibit8 * adv_sca_5 & + 3.0_wp * ibit7 * adv_sca_3 & + ibit6 * adv_sca_1 & ) * & ( sk(k+1,j,i) - sk(k,j,i) ) & - ( 5.0_wp * ibit8 * adv_sca_5 & + ibit7 * adv_sca_3 & ) * & ( sk(k_pp,j,i) - sk(k-1,j,i) ) & + ( ibit8 * adv_sca_5 & ) * & ( sk(k_ppp,j,i) - sk(k_mm,j,i) ) & ) ! !-- Calculate the divergence of the velocity field. A respective !-- correction is needed to overcome numerical instabilities caused !-- by a not sufficient reduction of divergences near topography. div = ( u(k,j,i+1) - u(k,j,i) ) * ddx & + ( v(k,j+1,i) - v(k,j,i) ) * ddy & + ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) tend(k,j,i) = tend(k,j,i) - ( & ( flux_r(k) + diss_r(k) - swap_flux_x_local(k,j,tn) - & swap_diss_x_local(k,j,tn) ) * ddx & + ( flux_n(k) + diss_n(k) - swap_flux_y_local(k,tn) - & swap_diss_y_local(k,tn) ) * ddy & + ( flux_t(k) + diss_t(k) - flux_d - diss_d & ) * ddzw(k) & ) + sk(k,j,i) * div swap_flux_y_local(k,tn) = flux_n(k) swap_diss_y_local(k,tn) = diss_n(k) swap_flux_x_local(k,j,tn) = flux_r(k) swap_diss_x_local(k,j,tn) = diss_r(k) flux_d = flux_t(k) diss_d = diss_t(k) ENDDO ! !-- Now compute the fluxes and tendency terms for the horizontal and !-- vertical parts above the top of the highest topography. No degradation !-- for the horizontal parts, but for the vertical it is stell needed. DO k = nzb_max+1, nzt u_comp = u(k,j,i+1) - u_gtrans flux_r(k) = u_comp * ( & 37.0_wp * ( sk(k,j,i+1) + sk(k,j,i) ) & - 8.0_wp * ( sk(k,j,i+2) + sk(k,j,i-1) ) & + ( sk(k,j,i+3) + sk(k,j,i-2) ) ) * adv_sca_5 diss_r(k) = -ABS( u_comp ) * ( & 10.0_wp * ( sk(k,j,i+1) - sk(k,j,i) ) & - 5.0_wp * ( sk(k,j,i+2) - sk(k,j,i-1) ) & + ( sk(k,j,i+3) - sk(k,j,i-2) ) ) * adv_sca_5 v_comp = v(k,j+1,i) - v_gtrans flux_n(k) = v_comp * ( & 37.0_wp * ( sk(k,j+1,i) + sk(k,j,i) ) & - 8.0_wp * ( sk(k,j+2,i) + sk(k,j-1,i) ) & + ( sk(k,j+3,i) + sk(k,j-2,i) ) ) * adv_sca_5 diss_n(k) = -ABS( v_comp ) * ( & 10.0_wp * ( sk(k,j+1,i) - sk(k,j,i) ) & - 5.0_wp * ( sk(k,j+2,i) - sk(k,j-1,i) ) & + ( sk(k,j+3,i) - sk(k,j-2,i) ) ) * adv_sca_5 ! !-- k index has to be modified near bottom and top, else array !-- subscripts will be exceeded. ibit8 = IBITS(wall_flags_0(k,j,i),8,1) ibit7 = IBITS(wall_flags_0(k,j,i),7,1) ibit6 = IBITS(wall_flags_0(k,j,i),6,1) k_ppp = k + 3 * ibit8 k_pp = k + 2 * ( 1 - ibit6 ) k_mm = k - 2 * ibit8 flux_t(k) = w(k,j,i) * ( & ( 37.0_wp * ibit8 * adv_sca_5 & + 7.0_wp * ibit7 * adv_sca_3 & + ibit6 * adv_sca_1 & ) * & ( sk(k+1,j,i) + sk(k,j,i) ) & - ( 8.0_wp * ibit8 * adv_sca_5 & + ibit7 * adv_sca_3 & ) * & ( sk(k_pp,j,i) + sk(k-1,j,i) ) & + ( ibit8 * adv_sca_5 & ) * ( sk(k_ppp,j,i)+ sk(k_mm,j,i) ) & ) diss_t(k) = -ABS( w(k,j,i) ) * ( & ( 10.0_wp * ibit8 * adv_sca_5 & + 3.0_wp * ibit7 * adv_sca_3 & + ibit6 * adv_sca_1 & ) * & ( sk(k+1,j,i) - sk(k,j,i) ) & - ( 5.0_wp * ibit8 * adv_sca_5 & + ibit7 * adv_sca_3 & ) * & ( sk(k_pp,j,i) - sk(k-1,j,i) ) & + ( ibit8 * adv_sca_5 & ) * & ( sk(k_ppp,j,i) - sk(k_mm,j,i) ) & ) ! !-- Calculate the divergence of the velocity field. A respective !-- correction is needed to overcome numerical instabilities introduced !-- by a not sufficient reduction of divergences near topography. div = ( u(k,j,i+1) - u(k,j,i) ) * ddx & + ( v(k,j+1,i) - v(k,j,i) ) * ddy & + ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) tend(k,j,i) = tend(k,j,i) - ( & ( flux_r(k) + diss_r(k) - swap_flux_x_local(k,j,tn) - & swap_diss_x_local(k,j,tn) ) * ddx & + ( flux_n(k) + diss_n(k) - swap_flux_y_local(k,tn) - & swap_diss_y_local(k,tn) ) * ddy & + ( flux_t(k) + diss_t(k) - flux_d - diss_d & ) * ddzw(k) & ) + sk(k,j,i) * div swap_flux_y_local(k,tn) = flux_n(k) swap_diss_y_local(k,tn) = diss_n(k) swap_flux_x_local(k,j,tn) = flux_r(k) swap_diss_x_local(k,j,tn) = diss_r(k) flux_d = flux_t(k) diss_d = diss_t(k) ENDDO ! !-- Evaluation of statistics SELECT CASE ( sk_char ) CASE ( 'pt' ) DO k = nzb, nzt sums_wspts_ws_l(k,tn) = sums_wspts_ws_l(k,tn) + & ( flux_t(k) + diss_t(k) ) & * weight_substep(intermediate_timestep_count) ENDDO CASE ( 'sa' ) DO k = nzb, nzt sums_wssas_ws_l(k,tn) = sums_wssas_ws_l(k,tn) + & ( flux_t(k) + diss_t(k) ) & * weight_substep(intermediate_timestep_count) ENDDO CASE ( 'q' ) DO k = nzb, nzt sums_wsqs_ws_l(k,tn) = sums_wsqs_ws_l(k,tn) + & ( flux_t(k) + diss_t(k) ) & * weight_substep(intermediate_timestep_count) ENDDO CASE ( 'qr' ) DO k = nzb, nzt sums_wsqrs_ws_l(k,tn) = sums_wsqrs_ws_l(k,tn) + & ( flux_t(k) + diss_t(k) ) & * weight_substep(intermediate_timestep_count) ENDDO CASE ( 'nr' ) DO k = nzb, nzt sums_wsnrs_ws_l(k,tn) = sums_wsnrs_ws_l(k,tn) + & ( flux_t(k) + diss_t(k) ) & * weight_substep(intermediate_timestep_count) ENDDO END SELECT END SUBROUTINE advec_s_ws_ij !------------------------------------------------------------------------------! ! Advection of u-component - Call for grid point i,j !------------------------------------------------------------------------------! SUBROUTINE advec_u_ws_ij( i, j, i_omp, tn ) USE arrays_3d, & ONLY: ddzw, diss_l_u, diss_s_u, flux_l_u, flux_s_u, tend, u, v, w USE constants, & ONLY: adv_mom_1, adv_mom_3, adv_mom_5 USE control_parameters, & ONLY: intermediate_timestep_count, u_gtrans, v_gtrans USE grid_variables, & ONLY: ddx, ddy USE indices, & ONLY: nxl, nxr, nyn, nys, nzb, nzb_max, nzt, wall_flags_0 USE kinds USE statistics, & ONLY: hom, sums_us2_ws_l, sums_wsus_ws_l, weight_substep IMPLICIT NONE INTEGER(iwp) :: i !: INTEGER(iwp) :: ibit9 !: INTEGER(iwp) :: ibit10 !: INTEGER(iwp) :: ibit11 !: INTEGER(iwp) :: ibit12 !: INTEGER(iwp) :: ibit13 !: INTEGER(iwp) :: ibit14 !: INTEGER(iwp) :: ibit15 !: INTEGER(iwp) :: ibit16 !: INTEGER(iwp) :: ibit17 !: INTEGER(iwp) :: i_omp !: INTEGER(iwp) :: j !: INTEGER(iwp) :: k !: INTEGER(iwp) :: k_mm !: INTEGER(iwp) :: k_pp !: INTEGER(iwp) :: k_ppp !: INTEGER(iwp) :: tn !: REAL(wp) :: diss_d !: REAL(wp) :: div !: REAL(wp) :: flux_d !: REAL(wp) :: gu !: REAL(wp) :: gv !: REAL(wp) :: u_comp_l !: REAL(wp) :: v_comp !: REAL(wp) :: w_comp !: REAL(wp), DIMENSION(nzb:nzt+1) :: diss_n !: REAL(wp), DIMENSION(nzb:nzt+1) :: diss_r !: REAL(wp), DIMENSION(nzb:nzt+1) :: diss_t !: REAL(wp), DIMENSION(nzb:nzt+1) :: flux_n !: REAL(wp), DIMENSION(nzb:nzt+1) :: flux_r !: REAL(wp), DIMENSION(nzb:nzt+1) :: flux_t !: REAL(wp), DIMENSION(nzb:nzt+1) :: u_comp !: gu = 2.0_wp * u_gtrans gv = 2.0_wp * v_gtrans ! !-- Compute southside fluxes for the respective boundary of PE IF ( j == nys ) THEN DO k = nzb+1, nzb_max ibit14 = IBITS(wall_flags_0(k,j,i),14,1) ibit13 = IBITS(wall_flags_0(k,j,i),13,1) ibit12 = IBITS(wall_flags_0(k,j,i),12,1) v_comp = v(k,j,i) + v(k,j,i-1) - gv flux_s_u(k,tn) = v_comp * ( & ( 37.0_wp * ibit14 * adv_mom_5 & + 7.0_wp * ibit13 * adv_mom_3 & + ibit12 * adv_mom_1 & ) * & ( u(k,j,i) + u(k,j-1,i) ) & - ( 8.0_wp * ibit14 * adv_mom_5 & + ibit13 * adv_mom_3 & ) * & ( u(k,j+1,i) + u(k,j-2,i) ) & + ( ibit14 * adv_mom_5 & ) * & ( u(k,j+2,i) + u(k,j-3,i) ) & ) diss_s_u(k,tn) = - ABS ( v_comp ) * ( & ( 10.0_wp * ibit14 * adv_mom_5 & + 3.0_wp * ibit13 * adv_mom_3 & + ibit12 * adv_mom_1 & ) * & ( u(k,j,i) - u(k,j-1,i) ) & - ( 5.0_wp * ibit14 * adv_mom_5 & + ibit13 * adv_mom_3 & ) * & ( u(k,j+1,i) - u(k,j-2,i) ) & + ( ibit14 * adv_mom_5 & ) * & ( u(k,j+2,i) - u(k,j-3,i) ) & ) ENDDO DO k = nzb_max+1, nzt v_comp = v(k,j,i) + v(k,j,i-1) - gv flux_s_u(k,tn) = v_comp * ( & 37.0_wp * ( u(k,j,i) + u(k,j-1,i) ) & - 8.0_wp * ( u(k,j+1,i) + u(k,j-2,i) ) & + ( u(k,j+2,i) + u(k,j-3,i) ) ) * adv_mom_5 diss_s_u(k,tn) = - ABS(v_comp) * ( & 10.0_wp * ( u(k,j,i) - u(k,j-1,i) ) & - 5.0_wp * ( u(k,j+1,i) - u(k,j-2,i) ) & + ( u(k,j+2,i) - u(k,j-3,i) ) ) * adv_mom_5 ENDDO ENDIF ! !-- Compute leftside fluxes for the respective boundary of PE IF ( i == i_omp ) THEN DO k = nzb+1, nzb_max ibit11 = IBITS(wall_flags_0(k,j,i),11,1) ibit10 = IBITS(wall_flags_0(k,j,i),10,1) ibit9 = IBITS(wall_flags_0(k,j,i),9,1) u_comp_l = u(k,j,i) + u(k,j,i-1) - gu flux_l_u(k,j,tn) = u_comp_l * ( & ( 37.0_wp * ibit11 * adv_mom_5 & + 7.0_wp * ibit10 * adv_mom_3 & + ibit9 * adv_mom_1 & ) * & ( u(k,j,i) + u(k,j,i-1) ) & - ( 8.0_wp * ibit11 * adv_mom_5 & + ibit10 * adv_mom_3 & ) * & ( u(k,j,i+1) + u(k,j,i-2) ) & + ( ibit11 * adv_mom_5 & ) * & ( u(k,j,i+2) + u(k,j,i-3) ) & ) diss_l_u(k,j,tn) = - ABS( u_comp_l ) * ( & ( 10.0_wp * ibit11 * adv_mom_5 & + 3.0_wp * ibit10 * adv_mom_3 & + ibit9 * adv_mom_1 & ) * & ( u(k,j,i) - u(k,j,i-1) ) & - ( 5.0_wp * ibit11 * adv_mom_5 & + ibit10 * adv_mom_3 & ) * & ( u(k,j,i+1) - u(k,j,i-2) ) & + ( ibit11 * adv_mom_5 & ) * & ( u(k,j,i+2) - u(k,j,i-3) ) & ) ENDDO DO k = nzb_max+1, nzt u_comp_l = u(k,j,i) + u(k,j,i-1) - gu flux_l_u(k,j,tn) = u_comp_l * ( & 37.0_wp * ( u(k,j,i) + u(k,j,i-1) ) & - 8.0_wp * ( u(k,j,i+1) + u(k,j,i-2) ) & + ( u(k,j,i+2) + u(k,j,i-3) ) ) * adv_mom_5 diss_l_u(k,j,tn) = - ABS(u_comp_l) * ( & 10.0_wp * ( u(k,j,i) - u(k,j,i-1) ) & - 5.0_wp * ( u(k,j,i+1) - u(k,j,i-2) ) & + ( u(k,j,i+2) - u(k,j,i-3) ) ) * adv_mom_5 ENDDO ENDIF flux_t(0) = 0.0_wp diss_t(0) = 0.0_wp flux_d = 0.0_wp diss_d = 0.0_wp ! !-- Now compute the fluxes tendency terms for the horizontal and !-- vertical parts. DO k = nzb+1, nzb_max ibit11 = IBITS(wall_flags_0(k,j,i),11,1) ibit10 = IBITS(wall_flags_0(k,j,i),10,1) ibit9 = IBITS(wall_flags_0(k,j,i),9,1) u_comp(k) = u(k,j,i+1) + u(k,j,i) flux_r(k) = ( u_comp(k) - gu ) * ( & ( 37.0_wp * ibit11 * adv_mom_5 & + 7.0_wp * ibit10 * adv_mom_3 & + ibit9 * adv_mom_1 & ) * & ( u(k,j,i+1) + u(k,j,i) ) & - ( 8.0_wp * ibit11 * adv_mom_5 & + ibit10 * adv_mom_3 & ) * & ( u(k,j,i+2) + u(k,j,i-1) ) & + ( ibit11 * adv_mom_5 & ) * & ( u(k,j,i+3) + u(k,j,i-2) ) & ) diss_r(k) = - ABS( u_comp(k) - gu ) * ( & ( 10.0_wp * ibit11 * adv_mom_5 & + 3.0_wp * ibit10 * adv_mom_3 & + ibit9 * adv_mom_1 & ) * & ( u(k,j,i+1) - u(k,j,i) ) & - ( 5.0_wp * ibit11 * adv_mom_5 & + ibit10 * adv_mom_3 & ) * & ( u(k,j,i+2) - u(k,j,i-1) ) & + ( ibit11 * adv_mom_5 & ) * & ( u(k,j,i+3) - u(k,j,i-2) ) & ) ibit14 = IBITS(wall_flags_0(k,j,i),14,1) ibit13 = IBITS(wall_flags_0(k,j,i),13,1) ibit12 = IBITS(wall_flags_0(k,j,i),12,1) v_comp = v(k,j+1,i) + v(k,j+1,i-1) - gv flux_n(k) = v_comp * ( & ( 37.0_wp * ibit14 * adv_mom_5 & + 7.0_wp * ibit13 * adv_mom_3 & + ibit12 * adv_mom_1 & ) * & ( u(k,j+1,i) + u(k,j,i) ) & - ( 8.0_wp * ibit14 * adv_mom_5 & + ibit13 * adv_mom_3 & ) * & ( u(k,j+2,i) + u(k,j-1,i) ) & + ( ibit14 * adv_mom_5 & ) * & ( u(k,j+3,i) + u(k,j-2,i) ) & ) diss_n(k) = - ABS ( v_comp ) * ( & ( 10.0_wp * ibit14 * adv_mom_5 & + 3.0_wp * ibit13 * adv_mom_3 & + ibit12 * adv_mom_1 & ) * & ( u(k,j+1,i) - u(k,j,i) ) & - ( 5.0_wp * ibit14 * adv_mom_5 & + ibit13 * adv_mom_3 & ) * & ( u(k,j+2,i) - u(k,j-1,i) ) & + ( ibit14 * adv_mom_5 & ) * & ( u(k,j+3,i) - u(k,j-2,i) ) & ) ! !-- k index has to be modified near bottom and top, else array !-- subscripts will be exceeded. ibit17 = IBITS(wall_flags_0(k,j,i),17,1) ibit16 = IBITS(wall_flags_0(k,j,i),16,1) ibit15 = IBITS(wall_flags_0(k,j,i),15,1) k_ppp = k + 3 * ibit17 k_pp = k + 2 * ( 1 - ibit15 ) k_mm = k - 2 * ibit17 w_comp = w(k,j,i) + w(k,j,i-1) flux_t(k) = w_comp * ( & ( 37.0_wp * ibit17 * adv_mom_5 & + 7.0_wp * ibit16 * adv_mom_3 & + ibit15 * adv_mom_1 & ) * & ( u(k+1,j,i) + u(k,j,i) ) & - ( 8.0_wp * ibit17 * adv_mom_5 & + ibit16 * adv_mom_3 & ) * & ( u(k_pp,j,i) + u(k-1,j,i) ) & + ( ibit17 * adv_mom_5 & ) * & ( u(k_ppp,j,i) + u(k_mm,j,i) ) & ) diss_t(k) = - ABS( w_comp ) * ( & ( 10.0_wp * ibit17 * adv_mom_5 & + 3.0_wp * ibit16 * adv_mom_3 & + ibit15 * adv_mom_1 & ) * & ( u(k+1,j,i) - u(k,j,i) ) & - ( 5.0_wp * ibit17 * adv_mom_5 & + ibit16 * adv_mom_3 & ) * & ( u(k_pp,j,i) - u(k-1,j,i) ) & + ( ibit17 * adv_mom_5 & ) * & ( u(k_ppp,j,i) - u(k_mm,j,i) ) & ) ! !-- Calculate the divergence of the velocity field. A respective !-- correction is needed to overcome numerical instabilities introduced !-- by a not sufficient reduction of divergences near topography. div = ( ( u_comp(k) - ( u(k,j,i) + u(k,j,i-1) ) ) * ddx & + ( v_comp + gv - ( v(k,j,i) + v(k,j,i-1 ) ) ) * ddy & + ( w_comp - ( w(k-1,j,i) + w(k-1,j,i-1) ) ) * ddzw(k) & ) * 0.5_wp tend(k,j,i) = tend(k,j,i) - ( & ( flux_r(k) + diss_r(k) & - flux_l_u(k,j,tn) - diss_l_u(k,j,tn) ) * ddx & + ( flux_n(k) + diss_n(k) & - flux_s_u(k,tn) - diss_s_u(k,tn) ) * ddy & + ( flux_t(k) + diss_t(k) & - flux_d - diss_d ) * ddzw(k) & ) + div * u(k,j,i) flux_l_u(k,j,tn) = flux_r(k) diss_l_u(k,j,tn) = diss_r(k) flux_s_u(k,tn) = flux_n(k) diss_s_u(k,tn) = diss_n(k) flux_d = flux_t(k) diss_d = diss_t(k) ! !-- Statistical Evaluation of u'u'. The factor has to be applied for !-- right evaluation when gallilei_trans = .T. . sums_us2_ws_l(k,tn) = sums_us2_ws_l(k,tn) & + ( flux_r(k) * & ( u_comp(k) - 2.0_wp * hom(k,1,1,0) ) & / ( u_comp(k) - gu + 1.0E-20_wp ) & + diss_r(k) * & ABS( u_comp(k) - 2.0_wp * hom(k,1,1,0) ) & / ( ABS( u_comp(k) - gu ) + 1.0E-20_wp ) ) & * weight_substep(intermediate_timestep_count) ! !-- Statistical Evaluation of w'u'. sums_wsus_ws_l(k,tn) = sums_wsus_ws_l(k,tn) & + ( flux_t(k) + diss_t(k) ) & * weight_substep(intermediate_timestep_count) ENDDO DO k = nzb_max+1, nzt u_comp(k) = u(k,j,i+1) + u(k,j,i) flux_r(k) = ( u_comp(k) - gu ) * ( & 37.0_wp * ( u(k,j,i+1) + u(k,j,i) ) & - 8.0_wp * ( u(k,j,i+2) + u(k,j,i-1) ) & + ( u(k,j,i+3) + u(k,j,i-2) ) ) * adv_mom_5 diss_r(k) = - ABS( u_comp(k) - gu ) * ( & 10.0_wp * ( u(k,j,i+1) - u(k,j,i) ) & - 5.0_wp * ( u(k,j,i+2) - u(k,j,i-1) ) & + ( u(k,j,i+3) - u(k,j,i-2) ) ) * adv_mom_5 v_comp = v(k,j+1,i) + v(k,j+1,i-1) - gv flux_n(k) = v_comp * ( & 37.0_wp * ( u(k,j+1,i) + u(k,j,i) ) & - 8.0_wp * ( u(k,j+2,i) + u(k,j-1,i) ) & + ( u(k,j+3,i) + u(k,j-2,i) ) ) * adv_mom_5 diss_n(k) = - ABS( v_comp ) * ( & 10.0_wp * ( u(k,j+1,i) - u(k,j,i) ) & - 5.0_wp * ( u(k,j+2,i) - u(k,j-1,i) ) & + ( u(k,j+3,i) - u(k,j-2,i) ) ) * adv_mom_5 ! !-- k index has to be modified near bottom and top, else array !-- subscripts will be exceeded. ibit17 = IBITS(wall_flags_0(k,j,i),17,1) ibit16 = IBITS(wall_flags_0(k,j,i),16,1) ibit15 = IBITS(wall_flags_0(k,j,i),15,1) k_ppp = k + 3 * ibit17 k_pp = k + 2 * ( 1 - ibit15 ) k_mm = k - 2 * ibit17 w_comp = w(k,j,i) + w(k,j,i-1) flux_t(k) = w_comp * ( & ( 37.0_wp * ibit17 * adv_mom_5 & + 7.0_wp * ibit16 * adv_mom_3 & + ibit15 * adv_mom_1 & ) * & ( u(k+1,j,i) + u(k,j,i) ) & - ( 8.0_wp * ibit17 * adv_mom_5 & + ibit16 * adv_mom_3 & ) * & ( u(k_pp,j,i) + u(k-1,j,i) ) & + ( ibit17 * adv_mom_5 & ) * & ( u(k_ppp,j,i) + u(k_mm,j,i) ) & ) diss_t(k) = - ABS( w_comp ) * ( & ( 10.0_wp * ibit17 * adv_mom_5 & + 3.0_wp * ibit16 * adv_mom_3 & + ibit15 * adv_mom_1 & ) * & ( u(k+1,j,i) - u(k,j,i) ) & - ( 5.0_wp * ibit17 * adv_mom_5 & + ibit16 * adv_mom_3 & ) * & ( u(k_pp,j,i) - u(k-1,j,i) ) & + ( ibit17 * adv_mom_5 & ) * & ( u(k_ppp,j,i) - u(k_mm,j,i) ) & ) ! !-- Calculate the divergence of the velocity field. A respective !-- correction is needed to overcome numerical instabilities introduced !-- by a not sufficient reduction of divergences near topography. div = ( ( u_comp(k) - ( u(k,j,i) + u(k,j,i-1) ) ) * ddx & + ( v_comp + gv - ( v(k,j,i) + v(k,j,i-1 ) ) ) * ddy & + ( w_comp - ( w(k-1,j,i) + w(k-1,j,i-1) ) ) * ddzw(k) & ) * 0.5_wp tend(k,j,i) = tend(k,j,i) - ( & ( flux_r(k) + diss_r(k) & - flux_l_u(k,j,tn) - diss_l_u(k,j,tn) ) * ddx & + ( flux_n(k) + diss_n(k) & - flux_s_u(k,tn) - diss_s_u(k,tn) ) * ddy & + ( flux_t(k) + diss_t(k) & - flux_d - diss_d ) * ddzw(k) & ) + div * u(k,j,i) flux_l_u(k,j,tn) = flux_r(k) diss_l_u(k,j,tn) = diss_r(k) flux_s_u(k,tn) = flux_n(k) diss_s_u(k,tn) = diss_n(k) flux_d = flux_t(k) diss_d = diss_t(k) ! !-- Statistical Evaluation of u'u'. The factor has to be applied for !-- right evaluation when gallilei_trans = .T. . sums_us2_ws_l(k,tn) = sums_us2_ws_l(k,tn) & + ( flux_r(k) * & ( u_comp(k) - 2.0_wp * hom(k,1,1,0) ) & / ( u_comp(k) - gu + 1.0E-20_wp ) & + diss_r(k) * & ABS( u_comp(k) - 2.0_wp * hom(k,1,1,0) ) & / ( ABS( u_comp(k) - gu ) + 1.0E-20_wp ) ) & * weight_substep(intermediate_timestep_count) ! !-- Statistical Evaluation of w'u'. sums_wsus_ws_l(k,tn) = sums_wsus_ws_l(k,tn) & + ( flux_t(k) + diss_t(k) ) & * weight_substep(intermediate_timestep_count) ENDDO sums_us2_ws_l(nzb,tn) = sums_us2_ws_l(nzb+1,tn) END SUBROUTINE advec_u_ws_ij !-----------------------------------------------------------------------------! ! Advection of v-component - Call for grid point i,j !-----------------------------------------------------------------------------! SUBROUTINE advec_v_ws_ij( i, j, i_omp, tn ) USE arrays_3d, & ONLY: ddzw, diss_l_v, diss_s_v, flux_l_v, flux_s_v, tend, u, v, w USE constants, & ONLY: adv_mom_1, adv_mom_3, adv_mom_5 USE control_parameters, & ONLY: intermediate_timestep_count, u_gtrans, v_gtrans USE grid_variables, & ONLY: ddx, ddy USE indices, & ONLY: nxl, nxr, nyn, nys, nysv, nzb, nzb_max, nzt, wall_flags_0 USE kinds USE statistics, & ONLY: hom, sums_vs2_ws_l, sums_wsvs_ws_l, weight_substep IMPLICIT NONE INTEGER(iwp) :: i !: INTEGER(iwp) :: ibit18 !: INTEGER(iwp) :: ibit19 !: INTEGER(iwp) :: ibit20 !: INTEGER(iwp) :: ibit21 !: INTEGER(iwp) :: ibit22 !: INTEGER(iwp) :: ibit23 !: INTEGER(iwp) :: ibit24 !: INTEGER(iwp) :: ibit25 !: INTEGER(iwp) :: ibit26 !: INTEGER(iwp) :: i_omp !: INTEGER(iwp) :: j !: INTEGER(iwp) :: k !: INTEGER(iwp) :: k_mm !: INTEGER(iwp) :: k_pp !: INTEGER(iwp) :: k_ppp !: INTEGER(iwp) :: tn !: REAL(wp) :: diss_d !: REAL(wp) :: div !: REAL(wp) :: flux_d !: REAL(wp) :: gu !: REAL(wp) :: gv !: REAL(wp) :: u_comp !: REAL(wp) :: v_comp_l !: REAL(wp) :: w_comp !: REAL(wp), DIMENSION(nzb:nzt+1) :: diss_n !: REAL(wp), DIMENSION(nzb:nzt+1) :: diss_r !: REAL(wp), DIMENSION(nzb:nzt+1) :: diss_t !: REAL(wp), DIMENSION(nzb:nzt+1) :: flux_n !: REAL(wp), DIMENSION(nzb:nzt+1) :: flux_r !: REAL(wp), DIMENSION(nzb:nzt+1) :: flux_t !: REAL(wp), DIMENSION(nzb:nzt+1) :: v_comp !: gu = 2.0_wp * u_gtrans gv = 2.0_wp * v_gtrans ! !-- Compute leftside fluxes for the respective boundary. IF ( i == i_omp ) THEN DO k = nzb+1, nzb_max ibit20 = IBITS(wall_flags_0(k,j,i),20,1) ibit19 = IBITS(wall_flags_0(k,j,i),19,1) ibit18 = IBITS(wall_flags_0(k,j,i),18,1) u_comp = u(k,j-1,i) + u(k,j,i) - gu flux_l_v(k,j,tn) = u_comp * ( & ( 37.0_wp * ibit20 * adv_mom_5 & + 7.0_wp * ibit19 * adv_mom_3 & + ibit18 * adv_mom_1 & ) * & ( v(k,j,i) + v(k,j,i-1) ) & - ( 8.0_wp * ibit20 * adv_mom_5 & + ibit19 * adv_mom_3 & ) * & ( v(k,j,i+1) + v(k,j,i-2) ) & + ( ibit20 * adv_mom_5 & ) * & ( v(k,j,i+2) + v(k,j,i-3) ) & ) diss_l_v(k,j,tn) = - ABS( u_comp ) * ( & ( 10.0_wp * ibit20 * adv_mom_5 & + 3.0_wp * ibit19 * adv_mom_3 & + ibit18 * adv_mom_1 & ) * & ( v(k,j,i) - v(k,j,i-1) ) & - ( 5.0_wp * ibit20 * adv_mom_5 & + ibit19 * adv_mom_3 & ) * & ( v(k,j,i+1) - v(k,j,i-2) ) & + ( ibit20 * adv_mom_5 & ) * & ( v(k,j,i+2) - v(k,j,i-3) ) & ) ENDDO DO k = nzb_max+1, nzt u_comp = u(k,j-1,i) + u(k,j,i) - gu flux_l_v(k,j,tn) = u_comp * ( & 37.0_wp * ( v(k,j,i) + v(k,j,i-1) ) & - 8.0_wp * ( v(k,j,i+1) + v(k,j,i-2) ) & + ( v(k,j,i+2) + v(k,j,i-3) ) ) * adv_mom_5 diss_l_v(k,j,tn) = - ABS( u_comp ) * ( & 10.0_wp * ( v(k,j,i) - v(k,j,i-1) ) & - 5.0_wp * ( v(k,j,i+1) - v(k,j,i-2) ) & + ( v(k,j,i+2) - v(k,j,i-3) ) ) * adv_mom_5 ENDDO ENDIF ! !-- Compute southside fluxes for the respective boundary. IF ( j == nysv ) THEN DO k = nzb+1, nzb_max ibit23 = IBITS(wall_flags_0(k,j,i),23,1) ibit22 = IBITS(wall_flags_0(k,j,i),22,1) ibit21 = IBITS(wall_flags_0(k,j,i),21,1) v_comp_l = v(k,j,i) + v(k,j-1,i) - gv flux_s_v(k,tn) = v_comp_l * ( & ( 37.0_wp * ibit23 * adv_mom_5 & + 7.0_wp * ibit22 * adv_mom_3 & + ibit21 * adv_mom_1 & ) * & ( v(k,j,i) + v(k,j-1,i) ) & - ( 8.0_wp * ibit23 * adv_mom_5 & + ibit22 * adv_mom_3 & ) * & ( v(k,j+1,i) + v(k,j-2,i) ) & + ( ibit23 * adv_mom_5 & ) * & ( v(k,j+2,i) + v(k,j-3,i) ) & ) diss_s_v(k,tn) = - ABS( v_comp_l ) * ( & ( 10.0_wp * ibit23 * adv_mom_5 & + 3.0_wp * ibit22 * adv_mom_3 & + ibit21 * adv_mom_1 & ) * & ( v(k,j,i) - v(k,j-1,i) ) & - ( 5.0_wp * ibit23 * adv_mom_5 & + ibit22 * adv_mom_3 & ) * & ( v(k,j+1,i) - v(k,j-2,i) ) & + ( ibit23 * adv_mom_5 & ) * & ( v(k,j+2,i) - v(k,j-3,i) ) & ) ENDDO DO k = nzb_max+1, nzt v_comp_l = v(k,j,i) + v(k,j-1,i) - gv flux_s_v(k,tn) = v_comp_l * ( & 37.0_wp * ( v(k,j,i) + v(k,j-1,i) ) & - 8.0_wp * ( v(k,j+1,i) + v(k,j-2,i) ) & + ( v(k,j+2,i) + v(k,j-3,i) ) ) * adv_mom_5 diss_s_v(k,tn) = - ABS( v_comp_l ) * ( & 10.0_wp * ( v(k,j,i) - v(k,j-1,i) ) & - 5.0_wp * ( v(k,j+1,i) - v(k,j-2,i) ) & + ( v(k,j+2,i) - v(k,j-3,i) ) ) * adv_mom_5 ENDDO ENDIF flux_t(0) = 0.0_wp diss_t(0) = 0.0_wp flux_d = 0.0_wp diss_d = 0.0_wp ! !-- Now compute the fluxes and tendency terms for the horizontal and !-- verical parts. DO k = nzb+1, nzb_max ibit20 = IBITS(wall_flags_0(k,j,i),20,1) ibit19 = IBITS(wall_flags_0(k,j,i),19,1) ibit18 = IBITS(wall_flags_0(k,j,i),18,1) u_comp = u(k,j-1,i+1) + u(k,j,i+1) - gu flux_r(k) = u_comp * ( & ( 37.0_wp * ibit20 * adv_mom_5 & + 7.0_wp * ibit19 * adv_mom_3 & + ibit18 * adv_mom_1 & ) * & ( v(k,j,i+1) + v(k,j,i) ) & - ( 8.0_wp * ibit20 * adv_mom_5 & + ibit19 * adv_mom_3 & ) * & ( v(k,j,i+2) + v(k,j,i-1) ) & + ( ibit20 * adv_mom_5 & ) * & ( v(k,j,i+3) + v(k,j,i-2) ) & ) diss_r(k) = - ABS( u_comp ) * ( & ( 10.0_wp * ibit20 * adv_mom_5 & + 3.0_wp * ibit19 * adv_mom_3 & + ibit18 * adv_mom_1 & ) * & ( v(k,j,i+1) - v(k,j,i) ) & - ( 5.0_wp * ibit20 * adv_mom_5 & + ibit19 * adv_mom_3 & ) * & ( v(k,j,i+2) - v(k,j,i-1) ) & + ( ibit20 * adv_mom_5 & ) * & ( v(k,j,i+3) - v(k,j,i-2) ) & ) ibit23 = IBITS(wall_flags_0(k,j,i),23,1) ibit22 = IBITS(wall_flags_0(k,j,i),22,1) ibit21 = IBITS(wall_flags_0(k,j,i),21,1) v_comp(k) = v(k,j+1,i) + v(k,j,i) flux_n(k) = ( v_comp(k) - gv ) * ( & ( 37.0_wp * ibit23 * adv_mom_5 & + 7.0_wp * ibit22 * adv_mom_3 & + ibit21 * adv_mom_1 & ) * & ( v(k,j+1,i) + v(k,j,i) ) & - ( 8.0_wp * ibit23 * adv_mom_5 & + ibit22 * adv_mom_3 & ) * & ( v(k,j+2,i) + v(k,j-1,i) ) & + ( ibit23 * adv_mom_5 & ) * & ( v(k,j+3,i) + v(k,j-2,i) ) & ) diss_n(k) = - ABS( v_comp(k) - gv ) * ( & ( 10.0_wp * ibit23 * adv_mom_5 & + 3.0_wp * ibit22 * adv_mom_3 & + ibit21 * adv_mom_1 & ) * & ( v(k,j+1,i) - v(k,j,i) ) & - ( 5.0_wp * ibit23 * adv_mom_5 & + ibit22 * adv_mom_3 & ) * & ( v(k,j+2,i) - v(k,j-1,i) ) & + ( ibit23 * adv_mom_5 & ) * & ( v(k,j+3,i) - v(k,j-2,i) ) & ) ! !-- k index has to be modified near bottom and top, else array !-- subscripts will be exceeded. ibit26 = IBITS(wall_flags_0(k,j,i),26,1) ibit25 = IBITS(wall_flags_0(k,j,i),25,1) ibit24 = IBITS(wall_flags_0(k,j,i),24,1) k_ppp = k + 3 * ibit26 k_pp = k + 2 * ( 1 - ibit24 ) k_mm = k - 2 * ibit26 w_comp = w(k,j-1,i) + w(k,j,i) flux_t(k) = w_comp * ( & ( 37.0_wp * ibit26 * adv_mom_5 & + 7.0_wp * ibit25 * adv_mom_3 & + ibit24 * adv_mom_1 & ) * & ( v(k+1,j,i) + v(k,j,i) ) & - ( 8.0_wp * ibit26 * adv_mom_5 & + ibit25 * adv_mom_3 & ) * & ( v(k_pp,j,i) + v(k-1,j,i) ) & + ( ibit26 * adv_mom_5 & ) * & ( v(k_ppp,j,i) + v(k_mm,j,i) ) & ) diss_t(k) = - ABS( w_comp ) * ( & ( 10.0_wp * ibit26 * adv_mom_5 & + 3.0_wp * ibit25 * adv_mom_3 & + ibit24 * adv_mom_1 & ) * & ( v(k+1,j,i) - v(k,j,i) ) & - ( 5.0_wp * ibit26 * adv_mom_5 & + ibit25 * adv_mom_3 & ) * & ( v(k_pp,j,i) - v(k-1,j,i) ) & + ( ibit26 * adv_mom_5 & ) * & ( v(k_ppp,j,i) - v(k_mm,j,i) ) & ) ! !-- Calculate the divergence of the velocity field. A respective !-- correction is needed to overcome numerical instabilities introduced !-- by a not sufficient reduction of divergences near topography. div = ( ( u_comp + gu - ( u(k,j-1,i) + u(k,j,i) ) ) * ddx & + ( v_comp(k) - ( v(k,j,i) + v(k,j-1,i) ) ) * ddy & + ( w_comp - ( w(k-1,j-1,i) + w(k-1,j,i) ) ) * ddzw(k) & ) * 0.5_wp tend(k,j,i) = tend(k,j,i) - ( & ( flux_r(k) + diss_r(k) & - flux_l_v(k,j,tn) - diss_l_v(k,j,tn) ) * ddx & + ( flux_n(k) + diss_n(k) & - flux_s_v(k,tn) - diss_s_v(k,tn) ) * ddy & + ( flux_t(k) + diss_t(k) & - flux_d - diss_d ) * ddzw(k) & ) + v(k,j,i) * div flux_l_v(k,j,tn) = flux_r(k) diss_l_v(k,j,tn) = diss_r(k) flux_s_v(k,tn) = flux_n(k) diss_s_v(k,tn) = diss_n(k) flux_d = flux_t(k) diss_d = diss_t(k) ! !-- Statistical Evaluation of v'v'. The factor has to be applied for !-- right evaluation when gallilei_trans = .T. . sums_vs2_ws_l(k,tn) = sums_vs2_ws_l(k,tn) & + ( flux_n(k) & * ( v_comp(k) - 2.0_wp * hom(k,1,2,0) ) & / ( v_comp(k) - gv + 1.0E-20_wp ) & + diss_n(k) & * ABS( v_comp(k) - 2.0_wp * hom(k,1,2,0) ) & / ( ABS( v_comp(k) - gv ) +1.0E-20_wp ) ) & * weight_substep(intermediate_timestep_count) ! !-- Statistical Evaluation of w'v'. sums_wsvs_ws_l(k,tn) = sums_wsvs_ws_l(k,tn) & + ( flux_t(k) + diss_t(k) ) & * weight_substep(intermediate_timestep_count) ENDDO DO k = nzb_max+1, nzt u_comp = u(k,j-1,i+1) + u(k,j,i+1) - gu flux_r(k) = u_comp * ( & 37.0_wp * ( v(k,j,i+1) + v(k,j,i) ) & - 8.0_wp * ( v(k,j,i+2) + v(k,j,i-1) ) & + ( v(k,j,i+3) + v(k,j,i-2) ) ) * adv_mom_5 diss_r(k) = - ABS( u_comp ) * ( & 10.0_wp * ( v(k,j,i+1) - v(k,j,i) ) & - 5.0_wp * ( v(k,j,i+2) - v(k,j,i-1) ) & + ( v(k,j,i+3) - v(k,j,i-2) ) ) * adv_mom_5 v_comp(k) = v(k,j+1,i) + v(k,j,i) flux_n(k) = ( v_comp(k) - gv ) * ( & 37.0_wp * ( v(k,j+1,i) + v(k,j,i) ) & - 8.0_wp * ( v(k,j+2,i) + v(k,j-1,i) ) & + ( v(k,j+3,i) + v(k,j-2,i) ) ) * adv_mom_5 diss_n(k) = - ABS( v_comp(k) - gv ) * ( & 10.0_wp * ( v(k,j+1,i) - v(k,j,i) ) & - 5.0_wp * ( v(k,j+2,i) - v(k,j-1,i) ) & + ( v(k,j+3,i) - v(k,j-2,i) ) ) * adv_mom_5 ! !-- k index has to be modified near bottom and top, else array !-- subscripts will be exceeded. ibit26 = IBITS(wall_flags_0(k,j,i),26,1) ibit25 = IBITS(wall_flags_0(k,j,i),25,1) ibit24 = IBITS(wall_flags_0(k,j,i),24,1) k_ppp = k + 3 * ibit26 k_pp = k + 2 * ( 1 - ibit24 ) k_mm = k - 2 * ibit26 w_comp = w(k,j-1,i) + w(k,j,i) flux_t(k) = w_comp * ( & ( 37.0_wp * ibit26 * adv_mom_5 & + 7.0_wp * ibit25 * adv_mom_3 & + ibit24 * adv_mom_1 & ) * & ( v(k+1,j,i) + v(k,j,i) ) & - ( 8.0_wp * ibit26 * adv_mom_5 & + ibit25 * adv_mom_3 & ) * & ( v(k_pp,j,i) + v(k-1,j,i) ) & + ( ibit26 * adv_mom_5 & ) * & ( v(k_ppp,j,i) + v(k_mm,j,i) ) & ) diss_t(k) = - ABS( w_comp ) * ( & ( 10.0_wp * ibit26 * adv_mom_5 & + 3.0_wp * ibit25 * adv_mom_3 & + ibit24 * adv_mom_1 & ) * & ( v(k+1,j,i) - v(k,j,i) ) & - ( 5.0_wp * ibit26 * adv_mom_5 & + ibit25 * adv_mom_3 & ) * & ( v(k_pp,j,i) - v(k-1,j,i) ) & + ( ibit26 * adv_mom_5 & ) * & ( v(k_ppp,j,i) - v(k_mm,j,i) ) & ) ! !-- Calculate the divergence of the velocity field. A respective !-- correction is needed to overcome numerical instabilities introduced !-- by a not sufficient reduction of divergences near topography. div = ( ( u_comp + gu - ( u(k,j-1,i) + u(k,j,i) ) ) * ddx & + ( v_comp(k) - ( v(k,j,i) + v(k,j-1,i) ) ) * ddy & + ( w_comp - ( w(k-1,j-1,i) + w(k-1,j,i) ) ) * ddzw(k) & ) * 0.5_wp tend(k,j,i) = tend(k,j,i) - ( & ( flux_r(k) + diss_r(k) & - flux_l_v(k,j,tn) - diss_l_v(k,j,tn) ) * ddx & + ( flux_n(k) + diss_n(k) & - flux_s_v(k,tn) - diss_s_v(k,tn) ) * ddy & + ( flux_t(k) + diss_t(k) & - flux_d - diss_d ) * ddzw(k) & ) + v(k,j,i) * div flux_l_v(k,j,tn) = flux_r(k) diss_l_v(k,j,tn) = diss_r(k) flux_s_v(k,tn) = flux_n(k) diss_s_v(k,tn) = diss_n(k) flux_d = flux_t(k) diss_d = diss_t(k) ! !-- Statistical Evaluation of v'v'. The factor has to be applied for !-- right evaluation when gallilei_trans = .T. . sums_vs2_ws_l(k,tn) = sums_vs2_ws_l(k,tn) & + ( flux_n(k) & * ( v_comp(k) - 2.0_wp * hom(k,1,2,0) ) & / ( v_comp(k) - gv + 1.0E-20_wp ) & + diss_n(k) & * ABS( v_comp(k) - 2.0_wp * hom(k,1,2,0) ) & / ( ABS( v_comp(k) - gv ) +1.0E-20_wp ) ) & * weight_substep(intermediate_timestep_count) ! !-- Statistical Evaluation of w'v'. sums_wsvs_ws_l(k,tn) = sums_wsvs_ws_l(k,tn) & + ( flux_t(k) + diss_t(k) ) & * weight_substep(intermediate_timestep_count) ENDDO sums_vs2_ws_l(nzb,tn) = sums_vs2_ws_l(nzb+1,tn) END SUBROUTINE advec_v_ws_ij !------------------------------------------------------------------------------! ! Advection of w-component - Call for grid point i,j !------------------------------------------------------------------------------! SUBROUTINE advec_w_ws_ij( i, j, i_omp, tn ) USE arrays_3d, & ONLY: ddzu, diss_l_w, diss_s_w, flux_l_w, flux_s_w, tend, u, v, w USE constants, & ONLY: adv_mom_1, adv_mom_3, adv_mom_5 USE control_parameters, & ONLY: intermediate_timestep_count, u_gtrans, v_gtrans USE grid_variables, & ONLY: ddx, ddy USE indices, & ONLY: nxl, nxr, nyn, nys, nzb, nzb_max, nzt, wall_flags_0, & wall_flags_00 USE kinds USE statistics, & ONLY: hom, sums_ws2_ws_l, weight_substep IMPLICIT NONE INTEGER(iwp) :: i !: INTEGER(iwp) :: ibit27 !: INTEGER(iwp) :: ibit28 !: INTEGER(iwp) :: ibit29 !: INTEGER(iwp) :: ibit30 !: INTEGER(iwp) :: ibit31 !: INTEGER(iwp) :: ibit32 !: INTEGER(iwp) :: ibit33 !: INTEGER(iwp) :: ibit34 !: INTEGER(iwp) :: ibit35 !: INTEGER(iwp) :: i_omp !: INTEGER(iwp) :: j !: INTEGER(iwp) :: k !: INTEGER(iwp) :: k_mm !: INTEGER(iwp) :: k_pp !: INTEGER(iwp) :: k_ppp !: INTEGER(iwp) :: tn !: REAL(wp) :: diss_d !: REAL(wp) :: div !: REAL(wp) :: flux_d !: REAL(wp) :: gu !: REAL(wp) :: gv !: REAL(wp) :: u_comp !: REAL(wp) :: v_comp !: REAL(wp) :: w_comp !: REAL(wp), DIMENSION(nzb:nzt+1) :: diss_n !: REAL(wp), DIMENSION(nzb:nzt+1) :: diss_r !: REAL(wp), DIMENSION(nzb:nzt+1) :: diss_t !: REAL(wp), DIMENSION(nzb:nzt+1) :: flux_n !: REAL(wp), DIMENSION(nzb:nzt+1) :: flux_r !: REAL(wp), DIMENSION(nzb:nzt+1) :: flux_t !: gu = 2.0_wp * u_gtrans gv = 2.0_wp * v_gtrans ! !-- Compute southside fluxes for the respective boundary. IF ( j == nys ) THEN DO k = nzb+1, nzb_max ibit32 = IBITS(wall_flags_00(k,j,i),0,1) ibit31 = IBITS(wall_flags_0(k,j,i),31,1) ibit30 = IBITS(wall_flags_0(k,j,i),30,1) v_comp = v(k+1,j,i) + v(k,j,i) - gv flux_s_w(k,tn) = v_comp * ( & ( 37.0_wp * ibit32 * adv_mom_5 & + 7.0_wp * ibit31 * adv_mom_3 & + ibit30 * adv_mom_1 & ) * & ( w(k,j,i) + w(k,j-1,i) ) & - ( 8.0_wp * ibit32 * adv_mom_5 & + ibit31 * adv_mom_3 & ) * & ( w(k,j+1,i) + w(k,j-2,i) ) & + ( ibit32 * adv_mom_5 & ) * & ( w(k,j+2,i) + w(k,j-3,i) ) & ) diss_s_w(k,tn) = - ABS( v_comp ) * ( & ( 10.0_wp * ibit32 * adv_mom_5 & + 3.0_wp * ibit31 * adv_mom_3 & + ibit30 * adv_mom_1 & ) * & ( w(k,j,i) - w(k,j-1,i) ) & - ( 5.0_wp * ibit32 * adv_mom_5 & + ibit31 * adv_mom_3 & ) * & ( w(k,j+1,i) - w(k,j-2,i) ) & + ( ibit32 * adv_mom_5 & ) * & ( w(k,j+2,i) - w(k,j-3,i) ) & ) ENDDO DO k = nzb_max+1, nzt v_comp = v(k+1,j,i) + v(k,j,i) - gv flux_s_w(k,tn) = v_comp * ( & 37.0_wp * ( w(k,j,i) + w(k,j-1,i) ) & - 8.0_wp * ( w(k,j+1,i) +w(k,j-2,i) ) & + ( w(k,j+2,i) + w(k,j-3,i) ) ) * adv_mom_5 diss_s_w(k,tn) = - ABS( v_comp ) * ( & 10.0_wp * ( w(k,j,i) - w(k,j-1,i) ) & - 5.0_wp * ( w(k,j+1,i) - w(k,j-2,i) ) & + ( w(k,j+2,i) - w(k,j-3,i) ) ) * adv_mom_5 ENDDO ENDIF ! !-- Compute leftside fluxes for the respective boundary. IF ( i == i_omp ) THEN DO k = nzb+1, nzb_max ibit29 = IBITS(wall_flags_0(k,j,i),29,1) ibit28 = IBITS(wall_flags_0(k,j,i),28,1) ibit27 = IBITS(wall_flags_0(k,j,i),27,1) u_comp = u(k+1,j,i) + u(k,j,i) - gu flux_l_w(k,j,tn) = u_comp * ( & ( 37.0_wp * ibit29 * adv_mom_5 & + 7.0_wp * ibit28 * adv_mom_3 & + ibit27 * adv_mom_1 & ) * & ( w(k,j,i) + w(k,j,i-1) ) & - ( 8.0_wp * ibit29 * adv_mom_5 & + ibit28 * adv_mom_3 & ) * & ( w(k,j,i+1) + w(k,j,i-2) ) & + ( ibit29 * adv_mom_5 & ) * & ( w(k,j,i+2) + w(k,j,i-3) ) & ) diss_l_w(k,j,tn) = - ABS( u_comp ) * ( & ( 10.0_wp * ibit29 * adv_mom_5 & + 3.0_wp * ibit28 * adv_mom_3 & + ibit27 * adv_mom_1 & ) * & ( w(k,j,i) - w(k,j,i-1) ) & - ( 5.0_wp * ibit29 * adv_mom_5 & + ibit28 * adv_mom_3 & ) * & ( w(k,j,i+1) - w(k,j,i-2) ) & + ( ibit29 * adv_mom_5 & ) * & ( w(k,j,i+2) - w(k,j,i-3) ) & ) ENDDO DO k = nzb_max+1, nzt u_comp = u(k+1,j,i) + u(k,j,i) - gu flux_l_w(k,j,tn) = u_comp * ( & 37.0_wp * ( w(k,j,i) + w(k,j,i-1) ) & - 8.0_wp * ( w(k,j,i+1) + w(k,j,i-2) ) & + ( w(k,j,i+2) + w(k,j,i-3) ) ) * adv_mom_5 diss_l_w(k,j,tn) = - ABS( u_comp ) * ( & 10.0_wp * ( w(k,j,i) - w(k,j,i-1) ) & - 5.0_wp * ( w(k,j,i+1) - w(k,j,i-2) ) & + ( w(k,j,i+2) - w(k,j,i-3) ) ) * adv_mom_5 ENDDO ENDIF ! !-- The lower flux has to be calculated explicetely for the tendency at !-- the first w-level. For topography wall this is done implicitely by !-- wall_flags_0. k = nzb + 1 w_comp = w(k,j,i) + w(k-1,j,i) flux_t(0) = w_comp * ( w(k,j,i) + w(k-1,j,i) ) * adv_mom_1 diss_t(0) = -ABS(w_comp) * ( w(k,j,i) - w(k-1,j,i) ) * adv_mom_1 flux_d = flux_t(0) diss_d = diss_t(0) ! !-- Now compute the fluxes and tendency terms for the horizontal !-- and vertical parts. DO k = nzb+1, nzb_max ibit29 = IBITS(wall_flags_0(k,j,i),29,1) ibit28 = IBITS(wall_flags_0(k,j,i),28,1) ibit27 = IBITS(wall_flags_0(k,j,i),27,1) u_comp = u(k+1,j,i+1) + u(k,j,i+1) - gu flux_r(k) = u_comp * ( & ( 37.0_wp * ibit29 * adv_mom_5 & + 7.0_wp * ibit28 * adv_mom_3 & + ibit27 * adv_mom_1 & ) * & ( w(k,j,i+1) + w(k,j,i) ) & - ( 8.0_wp * ibit29 * adv_mom_5 & + ibit28 * adv_mom_3 & ) * & ( w(k,j,i+2) + w(k,j,i-1) ) & + ( ibit29 * adv_mom_5 & ) * & ( w(k,j,i+3) + w(k,j,i-2) ) & ) diss_r(k) = - ABS( u_comp ) * ( & ( 10.0_wp * ibit29 * adv_mom_5 & + 3.0_wp * ibit28 * adv_mom_3 & + ibit27 * adv_mom_1 & ) * & ( w(k,j,i+1) - w(k,j,i) ) & - ( 5.0_wp * ibit29 * adv_mom_5 & + ibit28 * adv_mom_3 & ) * & ( w(k,j,i+2) - w(k,j,i-1) ) & + ( ibit29 * adv_mom_5 & ) * & ( w(k,j,i+3) - w(k,j,i-2) ) & ) ibit32 = IBITS(wall_flags_00(k,j,i),0,1) ibit31 = IBITS(wall_flags_0(k,j,i),31,1) ibit30 = IBITS(wall_flags_0(k,j,i),30,1) v_comp = v(k+1,j+1,i) + v(k,j+1,i) - gv flux_n(k) = v_comp * ( & ( 37.0_wp * ibit32 * adv_mom_5 & + 7.0_wp * ibit31 * adv_mom_3 & + ibit30 * adv_mom_1 & ) * & ( w(k,j+1,i) + w(k,j,i) ) & - ( 8.0_wp * ibit32 * adv_mom_5 & + ibit31 * adv_mom_3 & ) * & ( w(k,j+2,i) + w(k,j-1,i) ) & + ( ibit32 * adv_mom_5 & ) * & ( w(k,j+3,i) + w(k,j-2,i) ) & ) diss_n(k) = - ABS( v_comp ) * ( & ( 10.0_wp * ibit32 * adv_mom_5 & + 3.0_wp * ibit31 * adv_mom_3 & + ibit30 * adv_mom_1 & ) * & ( w(k,j+1,i) - w(k,j,i) ) & - ( 5.0_wp * ibit32 * adv_mom_5 & + ibit31 * adv_mom_3 & ) * & ( w(k,j+2,i) - w(k,j-1,i) ) & + ( ibit32 * adv_mom_5 & ) * & ( w(k,j+3,i) - w(k,j-2,i) ) & ) ! !-- k index has to be modified near bottom and top, else array !-- subscripts will be exceeded. ibit35 = IBITS(wall_flags_00(k,j,i),3,1) ibit34 = IBITS(wall_flags_00(k,j,i),2,1) ibit33 = IBITS(wall_flags_00(k,j,i),1,1) k_ppp = k + 3 * ibit35 k_pp = k + 2 * ( 1 - ibit33 ) k_mm = k - 2 * ibit35 w_comp = w(k+1,j,i) + w(k,j,i) flux_t(k) = w_comp * ( & ( 37.0_wp * ibit35 * adv_mom_5 & + 7.0_wp * ibit34 * adv_mom_3 & + ibit33 * adv_mom_1 & ) * & ( w(k+1,j,i) + w(k,j,i) ) & - ( 8.0_wp * ibit35 * adv_mom_5 & + ibit34 * adv_mom_3 & ) * & ( w(k_pp,j,i) + w(k-1,j,i) ) & + ( ibit35 * adv_mom_5 & ) * & ( w(k_ppp,j,i) + w(k_mm,j,i) ) & ) diss_t(k) = - ABS( w_comp ) * ( & ( 10.0_wp * ibit35 * adv_mom_5 & + 3.0_wp * ibit34 * adv_mom_3 & + ibit33 * adv_mom_1 & ) * & ( w(k+1,j,i) - w(k,j,i) ) & - ( 5.0_wp * ibit35 * adv_mom_5 & + ibit34 * adv_mom_3 & ) * & ( w(k_pp,j,i) - w(k-1,j,i) ) & + ( ibit35 * adv_mom_5 & ) * & ( w(k_ppp,j,i) - w(k_mm,j,i) ) & ) ! !-- Calculate the divergence of the velocity field. A respective !-- correction is needed to overcome numerical instabilities introduced !-- by a not sufficient reduction of divergences near topography. div = ( ( u_comp + gu - ( u(k+1,j,i) + u(k,j,i) ) ) * ddx & + ( v_comp + gv - ( v(k+1,j,i) + v(k,j,i) ) ) * ddy & + ( w_comp - ( w(k,j,i) + w(k-1,j,i) ) ) * ddzu(k+1) & ) * 0.5_wp tend(k,j,i) = tend(k,j,i) - ( & ( flux_r(k) + diss_r(k) & - flux_l_w(k,j,tn) - diss_l_w(k,j,tn) ) * ddx & + ( flux_n(k) + diss_n(k) & - flux_s_w(k,tn) - diss_s_w(k,tn) ) * ddy & + ( flux_t(k) + diss_t(k) & - flux_d - diss_d ) * ddzu(k+1) & ) + div * w(k,j,i) flux_l_w(k,j,tn) = flux_r(k) diss_l_w(k,j,tn) = diss_r(k) flux_s_w(k,tn) = flux_n(k) diss_s_w(k,tn) = diss_n(k) flux_d = flux_t(k) diss_d = diss_t(k) ! !-- Statistical Evaluation of w'w'. sums_ws2_ws_l(k,tn) = sums_ws2_ws_l(k,tn) & + ( flux_t(k) + diss_t(k) ) & * weight_substep(intermediate_timestep_count) ENDDO DO k = nzb_max+1, nzt u_comp = u(k+1,j,i+1) + u(k,j,i+1) - gu flux_r(k) = u_comp * ( & 37.0_wp * ( w(k,j,i+1) + w(k,j,i) ) & - 8.0_wp * ( w(k,j,i+2) + w(k,j,i-1) ) & + ( w(k,j,i+3) + w(k,j,i-2) ) ) * adv_mom_5 diss_r(k) = - ABS( u_comp ) * ( & 10.0_wp * ( w(k,j,i+1) - w(k,j,i) ) & - 5.0_wp * ( w(k,j,i+2) - w(k,j,i-1) ) & + ( w(k,j,i+3) - w(k,j,i-2) ) ) * adv_mom_5 v_comp = v(k+1,j+1,i) + v(k,j+1,i) - gv flux_n(k) = v_comp * ( & 37.0_wp * ( w(k,j+1,i) + w(k,j,i) ) & - 8.0_wp * ( w(k,j+2,i) + w(k,j-1,i) ) & + ( w(k,j+3,i) + w(k,j-2,i) ) ) * adv_mom_5 diss_n(k) = - ABS( v_comp ) * ( & 10.0_wp * ( w(k,j+1,i) - w(k,j,i) ) & - 5.0_wp * ( w(k,j+2,i) - w(k,j-1,i) ) & + ( w(k,j+3,i) - w(k,j-2,i) ) ) * adv_mom_5 ! !-- k index has to be modified near bottom and top, else array !-- subscripts will be exceeded. ibit35 = IBITS(wall_flags_00(k,j,i),3,1) ibit34 = IBITS(wall_flags_00(k,j,i),2,1) ibit33 = IBITS(wall_flags_00(k,j,i),1,1) k_ppp = k + 3 * ibit35 k_pp = k + 2 * ( 1 - ibit33 ) k_mm = k - 2 * ibit35 w_comp = w(k+1,j,i) + w(k,j,i) flux_t(k) = w_comp * ( & ( 37.0_wp * ibit35 * adv_mom_5 & + 7.0_wp * ibit34 * adv_mom_3 & + ibit33 * adv_mom_1 & ) * & ( w(k+1,j,i) + w(k,j,i) ) & - ( 8.0_wp * ibit35 * adv_mom_5 & + ibit34 * adv_mom_3 & ) * & ( w(k_pp,j,i) + w(k-1,j,i) ) & + ( ibit35 * adv_mom_5 & ) * & ( w(k_ppp,j,i) + w(k_mm,j,i) ) & ) diss_t(k) = - ABS( w_comp ) * ( & ( 10.0_wp * ibit35 * adv_mom_5 & + 3.0_wp * ibit34 * adv_mom_3 & + ibit33 * adv_mom_1 & ) * & ( w(k+1,j,i) - w(k,j,i) ) & - ( 5.0_wp * ibit35 * adv_mom_5 & + ibit34 * adv_mom_3 & ) * & ( w(k_pp,j,i) - w(k-1,j,i) ) & + ( ibit35 * adv_mom_5 & ) * & ( w(k_ppp,j,i) - w(k_mm,j,i) ) & ) ! !-- Calculate the divergence of the velocity field. A respective !-- correction is needed to overcome numerical instabilities introduced !-- by a not sufficient reduction of divergences near topography. div = ( ( u_comp + gu - ( u(k+1,j,i) + u(k,j,i) ) ) * ddx & + ( v_comp + gv - ( v(k+1,j,i) + v(k,j,i) ) ) * ddy & + ( w_comp - ( w(k,j,i) + w(k-1,j,i) ) ) * ddzu(k+1) & ) * 0.5_wp tend(k,j,i) = tend(k,j,i) - ( & ( flux_r(k) + diss_r(k) & - flux_l_w(k,j,tn) - diss_l_w(k,j,tn) ) * ddx & + ( flux_n(k) + diss_n(k) & - flux_s_w(k,tn) - diss_s_w(k,tn) ) * ddy & + ( flux_t(k) + diss_t(k) & - flux_d - diss_d ) * ddzu(k+1) & ) + div * w(k,j,i) flux_l_w(k,j,tn) = flux_r(k) diss_l_w(k,j,tn) = diss_r(k) flux_s_w(k,tn) = flux_n(k) diss_s_w(k,tn) = diss_n(k) flux_d = flux_t(k) diss_d = diss_t(k) ! !-- Statistical Evaluation of w'w'. sums_ws2_ws_l(k,tn) = sums_ws2_ws_l(k,tn) & + ( flux_t(k) + diss_t(k) ) & * weight_substep(intermediate_timestep_count) ENDDO END SUBROUTINE advec_w_ws_ij !------------------------------------------------------------------------------! ! Scalar advection - Call for all grid points !------------------------------------------------------------------------------! SUBROUTINE advec_s_ws( sk, sk_char ) USE arrays_3d, & ONLY: ddzw, tend, u, v, w USE constants, & ONLY: adv_sca_1, adv_sca_3, adv_sca_5 USE control_parameters, & ONLY: intermediate_timestep_count, u_gtrans, v_gtrans USE grid_variables, & ONLY: ddx, ddy USE indices, & ONLY: nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, nzb_max, & nzt, wall_flags_0 USE kinds USE statistics, & ONLY: sums_wspts_ws_l, sums_wsqs_ws_l, sums_wssas_ws_l, & sums_wsqrs_ws_l, sums_wsnrs_ws_l, weight_substep IMPLICIT NONE CHARACTER (LEN = *), INTENT(IN) :: sk_char !: INTEGER(iwp) :: i !: INTEGER(iwp) :: ibit0 !: INTEGER(iwp) :: ibit1 !: INTEGER(iwp) :: ibit2 !: INTEGER(iwp) :: ibit3 !: INTEGER(iwp) :: ibit4 !: INTEGER(iwp) :: ibit5 !: INTEGER(iwp) :: ibit6 !: INTEGER(iwp) :: ibit7 !: INTEGER(iwp) :: ibit8 !: INTEGER(iwp) :: j !: INTEGER(iwp) :: k !: INTEGER(iwp) :: k_mm !: INTEGER(iwp) :: k_pp !: INTEGER(iwp) :: k_ppp !: INTEGER(iwp) :: tn = 0 !: #if defined( __nopointer ) REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: sk !: #else REAL(wp), DIMENSION(:,:,:), POINTER :: sk !: #endif REAL(wp) :: diss_d !: REAL(wp) :: div !: REAL(wp) :: flux_d !: REAL(wp) :: u_comp !: REAL(wp) :: v_comp !: REAL(wp), DIMENSION(nzb:nzt) :: diss_n !: REAL(wp), DIMENSION(nzb:nzt) :: diss_r !: REAL(wp), DIMENSION(nzb:nzt) :: diss_t !: REAL(wp), DIMENSION(nzb:nzt) :: flux_n !: REAL(wp), DIMENSION(nzb:nzt) :: flux_r !: REAL(wp), DIMENSION(nzb:nzt) :: flux_t !: REAL(wp), DIMENSION(nzb+1:nzt) :: swap_diss_y_local !: REAL(wp), DIMENSION(nzb+1:nzt) :: swap_flux_y_local !: REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) :: swap_diss_x_local !: REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) :: swap_flux_x_local !: ! !-- Compute the fluxes for the whole left boundary of the processor domain. i = nxl DO j = nys, nyn DO k = nzb+1, nzb_max ibit2 = IBITS(wall_flags_0(k,j,i),2,1) ibit1 = IBITS(wall_flags_0(k,j,i),1,1) ibit0 = IBITS(wall_flags_0(k,j,i),0,1) u_comp = u(k,j,i) - u_gtrans swap_flux_x_local(k,j) = u_comp * ( & ( 37.0_wp * ibit2 * adv_sca_5 & + 7.0_wp * ibit1 * adv_sca_3 & + ibit0 * adv_sca_1 & ) * & ( sk(k,j,i) + sk(k,j,i-1) ) & - ( 8.0_wp * ibit2 * adv_sca_5 & + ibit1 * adv_sca_3 & ) * & ( sk(k,j,i+1) + sk(k,j,i-2) ) & + ( ibit2 * adv_sca_5 & ) * & ( sk(k,j,i+2) + sk(k,j,i-3) ) & ) swap_diss_x_local(k,j) = -ABS( u_comp ) * ( & ( 10.0_wp * ibit2 * adv_sca_5 & + 3.0_wp * ibit1 * adv_sca_3 & + ibit0 * adv_sca_1 & ) * & ( sk(k,j,i) - sk(k,j,i-1) ) & - ( 5.0_wp * ibit2 * adv_sca_5 & + ibit1 * adv_sca_3 & ) * & ( sk(k,j,i+1) - sk(k,j,i-2) ) & + ( ibit2 * adv_sca_5 & ) * & ( sk(k,j,i+2) - sk(k,j,i-3) ) & ) ENDDO DO k = nzb_max+1, nzt u_comp = u(k,j,i) - u_gtrans swap_flux_x_local(k,j) = u_comp * ( & 37.0_wp * ( sk(k,j,i) + sk(k,j,i-1) ) & - 8.0_wp * ( sk(k,j,i+1) + sk(k,j,i-2) ) & + ( sk(k,j,i+2) + sk(k,j,i-3) ) & ) * adv_sca_5 swap_diss_x_local(k,j) = -ABS( u_comp ) * ( & 10.0_wp * ( sk(k,j,i) - sk(k,j,i-1) ) & - 5.0_wp * ( sk(k,j,i+1) - sk(k,j,i-2) ) & + ( sk(k,j,i+2) - sk(k,j,i-3) ) & ) * adv_sca_5 ENDDO ENDDO DO i = nxl, nxr j = nys DO k = nzb+1, nzb_max ibit5 = IBITS(wall_flags_0(k,j,i),5,1) ibit4 = IBITS(wall_flags_0(k,j,i),4,1) ibit3 = IBITS(wall_flags_0(k,j,i),3,1) v_comp = v(k,j,i) - v_gtrans swap_flux_y_local(k) = v_comp * ( & ( 37.0_wp * ibit5 * adv_sca_5 & + 7.0_wp * ibit4 * adv_sca_3 & + ibit3 * adv_sca_1 & ) * & ( sk(k,j,i) + sk(k,j-1,i) ) & - ( 8.0_wp * ibit5 * adv_sca_5 & + ibit4 * adv_sca_3 & ) * & ( sk(k,j+1,i) + sk(k,j-2,i) ) & + ( ibit5 * adv_sca_5 & ) * & ( sk(k,j+2,i) + sk(k,j-3,i) ) & ) swap_diss_y_local(k) = -ABS( v_comp ) * ( & ( 10.0_wp * ibit5 * adv_sca_5 & + 3.0_wp * ibit4 * adv_sca_3 & + ibit3 * adv_sca_1 & ) * & ( sk(k,j,i) - sk(k,j-1,i) ) & - ( 5.0_wp * ibit5 * adv_sca_5 & + ibit4 * adv_sca_3 & ) * & ( sk(k,j+1,i) - sk(k,j-2,i) ) & + ( ibit5 * adv_sca_5 & ) * & ( sk(k,j+2,i) - sk(k,j-3,i) ) & ) ENDDO ! !-- Above to the top of the highest topography. No degradation necessary. DO k = nzb_max+1, nzt v_comp = v(k,j,i) - v_gtrans swap_flux_y_local(k) = v_comp * ( & 37.0_wp * ( sk(k,j,i) + sk(k,j-1,i) ) & - 8.0_wp * ( sk(k,j+1,i) + sk(k,j-2,i) ) & + ( sk(k,j+2,i) + sk(k,j-3,i) ) & ) * adv_sca_5 swap_diss_y_local(k) = -ABS( v_comp ) * ( & 10.0_wp * ( sk(k,j,i) - sk(k,j-1,i) ) & - 5.0_wp * ( sk(k,j+1,i) - sk(k,j-2,i) ) & + sk(k,j+2,i) - sk(k,j-3,i) & ) * adv_sca_5 ENDDO DO j = nys, nyn flux_t(0) = 0.0_wp diss_t(0) = 0.0_wp flux_d = 0.0_wp diss_d = 0.0_wp DO k = nzb+1, nzb_max ibit2 = IBITS(wall_flags_0(k,j,i),2,1) ibit1 = IBITS(wall_flags_0(k,j,i),1,1) ibit0 = IBITS(wall_flags_0(k,j,i),0,1) u_comp = u(k,j,i+1) - u_gtrans flux_r(k) = u_comp * ( & ( 37.0_wp * ibit2 * adv_sca_5 & + 7.0_wp * ibit1 * adv_sca_3 & + ibit0 * adv_sca_1 & ) * & ( sk(k,j,i+1) + sk(k,j,i) ) & - ( 8.0_wp * ibit2 * adv_sca_5 & + ibit1 * adv_sca_3 & ) * & ( sk(k,j,i+2) + sk(k,j,i-1) ) & + ( ibit2 * adv_sca_5 & ) * & ( sk(k,j,i+3) + sk(k,j,i-2) ) & ) diss_r(k) = -ABS( u_comp ) * ( & ( 10.0_wp * ibit2 * adv_sca_5 & + 3.0_wp * ibit1 * adv_sca_3 & + ibit0 * adv_sca_1 & ) * & ( sk(k,j,i+1) - sk(k,j,i) ) & - ( 5.0_wp * ibit2 * adv_sca_5 & + ibit1 * adv_sca_3 & ) * & ( sk(k,j,i+2) - sk(k,j,i-1) ) & + ( ibit2 * adv_sca_5 & ) * & ( sk(k,j,i+3) - sk(k,j,i-2) ) & ) ibit5 = IBITS(wall_flags_0(k,j,i),5,1) ibit4 = IBITS(wall_flags_0(k,j,i),4,1) ibit3 = IBITS(wall_flags_0(k,j,i),3,1) v_comp = v(k,j+1,i) - v_gtrans flux_n(k) = v_comp * ( & ( 37.0_wp * ibit5 * adv_sca_5 & + 7.0_wp * ibit4 * adv_sca_3 & + ibit3 * adv_sca_1 & ) * & ( sk(k,j+1,i) + sk(k,j,i) ) & - ( 8.0_wp * ibit5 * adv_sca_5 & + ibit4 * adv_sca_3 & ) * & ( sk(k,j+2,i) + sk(k,j-1,i) ) & + ( ibit5 * adv_sca_5 & ) * & ( sk(k,j+3,i) + sk(k,j-2,i) ) & ) diss_n(k) = -ABS( v_comp ) * ( & ( 10.0_wp * ibit5 * adv_sca_5 & + 3.0_wp * ibit4 * adv_sca_3 & + ibit3 * adv_sca_1 & ) * & ( sk(k,j+1,i) - sk(k,j,i) ) & - ( 5.0_wp * ibit5 * adv_sca_5 & + ibit4 * adv_sca_3 & ) * & ( sk(k,j+2,i) - sk(k,j-1,i) ) & + ( ibit5 * adv_sca_5 & ) * & ( sk(k,j+3,i) - sk(k,j-2,i) ) & ) ! !-- k index has to be modified near bottom and top, else array !-- subscripts will be exceeded. ibit8 = IBITS(wall_flags_0(k,j,i),8,1) ibit7 = IBITS(wall_flags_0(k,j,i),7,1) ibit6 = IBITS(wall_flags_0(k,j,i),6,1) k_ppp = k + 3 * ibit8 k_pp = k + 2 * ( 1 - ibit6 ) k_mm = k - 2 * ibit8 flux_t(k) = w(k,j,i) * ( & ( 37.0_wp * ibit8 * adv_sca_5 & + 7.0_wp * ibit7 * adv_sca_3 & + ibit6 * adv_sca_1 & ) * & ( sk(k+1,j,i) + sk(k,j,i) ) & - ( 8.0_wp * ibit8 * adv_sca_5 & + ibit7 * adv_sca_3 & ) * & ( sk(k_pp,j,i) + sk(k-1,j,i) ) & + ( ibit8 * adv_sca_5 & ) * ( sk(k_ppp,j,i)+ sk(k_mm,j,i) ) & ) diss_t(k) = -ABS( w(k,j,i) ) * ( & ( 10.0_wp * ibit8 * adv_sca_5 & + 3.0_wp * ibit7 * adv_sca_3 & + ibit6 * adv_sca_1 & ) * & ( sk(k+1,j,i) - sk(k,j,i) ) & - ( 5.0_wp * ibit8 * adv_sca_5 & + ibit7 * adv_sca_3 & ) * & ( sk(k_pp,j,i) - sk(k-1,j,i) ) & + ( ibit8 * adv_sca_5 & ) * & ( sk(k_ppp,j,i) - sk(k_mm,j,i) ) & ) ! !-- Calculate the divergence of the velocity field. A respective !-- correction is needed to overcome numerical instabilities caused !-- by a not sufficient reduction of divergences near topography. div = ( u(k,j,i+1) - u(k,j,i) ) * ddx & + ( v(k,j+1,i) - v(k,j,i) ) * ddy & + ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) tend(k,j,i) = tend(k,j,i) - ( & ( flux_r(k) + diss_r(k) - swap_flux_x_local(k,j) - & swap_diss_x_local(k,j) ) * ddx & + ( flux_n(k) + diss_n(k) - swap_flux_y_local(k) - & swap_diss_y_local(k) ) * ddy & + ( flux_t(k) + diss_t(k) - flux_d - diss_d & ) * ddzw(k) & ) + sk(k,j,i) * div swap_flux_y_local(k) = flux_n(k) swap_diss_y_local(k) = diss_n(k) swap_flux_x_local(k,j) = flux_r(k) swap_diss_x_local(k,j) = diss_r(k) flux_d = flux_t(k) diss_d = diss_t(k) ENDDO DO k = nzb_max+1, nzt u_comp = u(k,j,i+1) - u_gtrans flux_r(k) = u_comp * ( & 37.0_wp * ( sk(k,j,i+1) + sk(k,j,i) ) & - 8.0_wp * ( sk(k,j,i+2) + sk(k,j,i-1) ) & + ( sk(k,j,i+3) + sk(k,j,i-2) ) ) * adv_sca_5 diss_r(k) = -ABS( u_comp ) * ( & 10.0_wp * ( sk(k,j,i+1) - sk(k,j,i) ) & - 5.0_wp * ( sk(k,j,i+2) - sk(k,j,i-1) ) & + ( sk(k,j,i+3) - sk(k,j,i-2) ) ) * adv_sca_5 v_comp = v(k,j+1,i) - v_gtrans flux_n(k) = v_comp * ( & 37.0_wp * ( sk(k,j+1,i) + sk(k,j,i) ) & - 8.0_wp * ( sk(k,j+2,i) + sk(k,j-1,i) ) & + ( sk(k,j+3,i) + sk(k,j-2,i) ) ) * adv_sca_5 diss_n(k) = -ABS( v_comp ) * ( & 10.0_wp * ( sk(k,j+1,i) - sk(k,j,i) ) & - 5.0_wp * ( sk(k,j+2,i) - sk(k,j-1,i) ) & + ( sk(k,j+3,i) - sk(k,j-2,i) ) ) * adv_sca_5 ! !-- k index has to be modified near bottom and top, else array !-- subscripts will be exceeded. ibit8 = IBITS(wall_flags_0(k,j,i),8,1) ibit7 = IBITS(wall_flags_0(k,j,i),7,1) ibit6 = IBITS(wall_flags_0(k,j,i),6,1) k_ppp = k + 3 * ibit8 k_pp = k + 2 * ( 1 - ibit6 ) k_mm = k - 2 * ibit8 flux_t(k) = w(k,j,i) * ( & ( 37.0_wp * ibit8 * adv_sca_5 & + 7.0_wp * ibit7 * adv_sca_3 & + ibit6 * adv_sca_1 & ) * & ( sk(k+1,j,i) + sk(k,j,i) ) & - ( 8.0_wp * ibit8 * adv_sca_5 & + ibit7 * adv_sca_3 & ) * & ( sk(k_pp,j,i) + sk(k-1,j,i) ) & + ( ibit8 * adv_sca_5 & ) * ( sk(k_ppp,j,i)+ sk(k_mm,j,i) ) & ) diss_t(k) = -ABS( w(k,j,i) ) * ( & ( 10.0_wp * ibit8 * adv_sca_5 & + 3.0_wp * ibit7 * adv_sca_3 & + ibit6 * adv_sca_1 & ) * & ( sk(k+1,j,i) - sk(k,j,i) ) & - ( 5.0_wp * ibit8 * adv_sca_5 & + ibit7 * adv_sca_3 & ) * & ( sk(k_pp,j,i) - sk(k-1,j,i) ) & + ( ibit8 * adv_sca_5 & ) * & ( sk(k_ppp,j,i) - sk(k_mm,j,i) ) & ) ! !-- Calculate the divergence of the velocity field. A respective !-- correction is needed to overcome numerical instabilities introduced !-- by a not sufficient reduction of divergences near topography. div = ( u(k,j,i+1) - u(k,j,i) ) * ddx & + ( v(k,j+1,i) - v(k,j,i) ) * ddy & + ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) tend(k,j,i) = tend(k,j,i) - ( & ( flux_r(k) + diss_r(k) - swap_flux_x_local(k,j) - & swap_diss_x_local(k,j) ) * ddx & + ( flux_n(k) + diss_n(k) - swap_flux_y_local(k) - & swap_diss_y_local(k) ) * ddy & + ( flux_t(k) + diss_t(k) - flux_d - diss_d & ) * ddzw(k) & ) + sk(k,j,i) * div swap_flux_y_local(k) = flux_n(k) swap_diss_y_local(k) = diss_n(k) swap_flux_x_local(k,j) = flux_r(k) swap_diss_x_local(k,j) = diss_r(k) flux_d = flux_t(k) diss_d = diss_t(k) ENDDO ! !-- evaluation of statistics SELECT CASE ( sk_char ) CASE ( 'pt' ) DO k = nzb, nzt sums_wspts_ws_l(k,tn) = sums_wspts_ws_l(k,tn) & + ( flux_t(k) + diss_t(k) ) & * weight_substep(intermediate_timestep_count) ENDDO CASE ( 'sa' ) DO k = nzb, nzt sums_wssas_ws_l(k,tn) = sums_wssas_ws_l(k,tn) & + ( flux_t(k) + diss_t(k) ) & * weight_substep(intermediate_timestep_count) ENDDO CASE ( 'q' ) DO k = nzb, nzt sums_wsqs_ws_l(k,tn) = sums_wsqs_ws_l(k,tn) & + ( flux_t(k) + diss_t(k) ) & * weight_substep(intermediate_timestep_count) ENDDO CASE ( 'qr' ) DO k = nzb, nzt sums_wsqrs_ws_l(k,tn) = sums_wsqrs_ws_l(k,tn) & + ( flux_t(k) + diss_t(k) ) & * weight_substep(intermediate_timestep_count) ENDDO CASE ( 'nr' ) DO k = nzb, nzt sums_wsnrs_ws_l(k,tn) = sums_wsnrs_ws_l(k,tn) & + ( flux_t(k) + diss_t(k) ) & * weight_substep(intermediate_timestep_count) ENDDO END SELECT ENDDO ENDDO END SUBROUTINE advec_s_ws !------------------------------------------------------------------------------! ! Scalar advection - Call for all grid points - accelerator version !------------------------------------------------------------------------------! SUBROUTINE advec_s_ws_acc ( sk, sk_char ) USE arrays_3d, & ONLY: ddzw, tend, u, v, w USE constants, & ONLY: adv_sca_1, adv_sca_3, adv_sca_5 USE control_parameters, & ONLY: intermediate_timestep_count, u_gtrans, v_gtrans USE grid_variables, & ONLY: ddx, ddy USE indices, & ONLY: i_left, i_right, j_north, j_south, nxlg, nxrg, nyng, nysg, & nzb, nzb_max, nzt, wall_flags_0 USE kinds ! USE statistics, & ! ONLY: sums_wspts_ws_l, sums_wsqs_ws_l, sums_wssas_ws_l, & ! sums_wsqrs_ws_l, sums_wsnrs_ws_l, weight_substep IMPLICIT NONE CHARACTER (LEN = *), INTENT(IN) :: sk_char !: INTEGER(iwp) :: i !: INTEGER(iwp) :: ibit0 !: INTEGER(iwp) :: ibit1 !: INTEGER(iwp) :: ibit2 !: INTEGER(iwp) :: ibit3 !: INTEGER(iwp) :: ibit4 !: INTEGER(iwp) :: ibit5 !: INTEGER(iwp) :: ibit6 !: INTEGER(iwp) :: ibit7 !: INTEGER(iwp) :: ibit8 !: INTEGER(iwp) :: j !: INTEGER(iwp) :: k !: INTEGER(iwp) :: k_mm !: INTEGER(iwp) :: k_mmm !: INTEGER(iwp) :: k_pp !: INTEGER(iwp) :: k_ppp !: INTEGER(iwp) :: tn = 0 !: REAL(wp) :: diss_d !: REAL(wp) :: diss_l !: REAL(wp) :: diss_n !: REAL(wp) :: diss_r !: REAL(wp) :: diss_s !: REAL(wp) :: diss_t !: REAL(wp) :: div !: REAL(wp) :: flux_d !: REAL(wp) :: flux_l !: REAL(wp) :: flux_n !: REAL(wp) :: flux_r !: REAL(wp) :: flux_s !: REAL(wp) :: flux_t !: REAL(wp) :: u_comp !: REAL(wp) :: v_comp !: REAL(wp), INTENT(IN), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: sk !: ! !-- Computation of fluxes and tendency terms !$acc kernels present( ddzw, sk, tend, u, v, w, wall_flags_0 ) DO i = i_left, i_right DO j = j_south, j_north DO k = nzb+1, nzt ibit2 = IBITS(wall_flags_0(k,j,i),2,1) ibit1 = IBITS(wall_flags_0(k,j,i),1,1) ibit0 = IBITS(wall_flags_0(k,j,i),0,1) u_comp = u(k,j,i) - u_gtrans flux_l = u_comp * ( & ( 37.0_wp * ibit2 * adv_sca_5 & + 7.0_wp * ibit1 * adv_sca_3 & + ibit0 * adv_sca_1 & ) * & ( sk(k,j,i) + sk(k,j,i-1) ) & - ( 8.0_wp * ibit2 * adv_sca_5 & + ibit1 * adv_sca_3 & ) * & ( sk(k,j,i+1) + sk(k,j,i-2) ) & + ( ibit2 * adv_sca_5 & ) * & ( sk(k,j,i+2) + sk(k,j,i-3) ) & ) diss_l = -ABS( u_comp ) * ( & ( 10.0_wp * ibit2 * adv_sca_5 & + 3.0_wp * ibit1 * adv_sca_3 & + ibit0 * adv_sca_1 & ) * & ( sk(k,j,i) - sk(k,j,i-1) ) & - ( 5.0_wp * ibit2 * adv_sca_5 & + ibit1 * adv_sca_3 & ) * & ( sk(k,j,i+1) - sk(k,j,i-2) ) & + ( ibit2 * adv_sca_5 & ) * & ( sk(k,j,i+2) - sk(k,j,i-3) ) & ) u_comp = u(k,j,i+1) - u_gtrans flux_r = u_comp * ( & ( 37.0_wp * ibit2 * adv_sca_5 & + 7.0_wp * ibit1 * adv_sca_3 & + ibit0 * adv_sca_1 & ) * & ( sk(k,j,i+1) + sk(k,j,i) ) & - ( 8.0_wp * ibit2 * adv_sca_5 & + ibit1 * adv_sca_3 & ) * & ( sk(k,j,i+2) + sk(k,j,i-1) ) & + ( ibit2 * adv_sca_5 & ) * & ( sk(k,j,i+3) + sk(k,j,i-2) ) & ) diss_r = -ABS( u_comp ) * ( & ( 10.0_wp * ibit2 * adv_sca_5 & + 3.0_wp * ibit1 * adv_sca_3 & + ibit0 * adv_sca_1 & ) * & ( sk(k,j,i+1) - sk(k,j,i) ) & - ( 5.0_wp * ibit2 * adv_sca_5 & + ibit1 * adv_sca_3 & ) * & ( sk(k,j,i+2) - sk(k,j,i-1) ) & + ( ibit2 * adv_sca_5 & ) * & ( sk(k,j,i+3) - sk(k,j,i-2) ) & ) ibit5 = IBITS(wall_flags_0(k,j,i),5,1) ibit4 = IBITS(wall_flags_0(k,j,i),4,1) ibit3 = IBITS(wall_flags_0(k,j,i),3,1) v_comp = v(k,j,i) - v_gtrans flux_s = v_comp * ( & ( 37.0_wp * ibit5 * adv_sca_5 & + 7.0_wp * ibit4 * adv_sca_3 & + ibit3 * adv_sca_1 & ) * & ( sk(k,j,i) + sk(k,j-1,i) ) & - ( 8.0_wp * ibit5 * adv_sca_5 & + ibit4 * adv_sca_3 & ) * & ( sk(k,j+1,i) + sk(k,j-2,i) ) & + ( ibit5 * adv_sca_5 & ) * & ( sk(k,j+2,i) + sk(k,j-3,i) ) & ) diss_s = -ABS( v_comp ) * ( & ( 10.0_wp * ibit5 * adv_sca_5 & + 3.0_wp * ibit4 * adv_sca_3 & + ibit3 * adv_sca_1 & ) * & ( sk(k,j,i) - sk(k,j-1,i) ) & - ( 5.0_wp * ibit5 * adv_sca_5 & + ibit4 * adv_sca_3 & ) * & ( sk(k,j+1,i) - sk(k,j-2,i) ) & + ( ibit5 * adv_sca_5 & ) * & ( sk(k,j+2,i) - sk(k,j-3,i) ) & ) v_comp = v(k,j+1,i) - v_gtrans flux_n = v_comp * ( & ( 37.0_wp * ibit5 * adv_sca_5 & + 7.0_wp * ibit4 * adv_sca_3 & + ibit3 * adv_sca_1 & ) * & ( sk(k,j+1,i) + sk(k,j,i) ) & - ( 8.0_wp * ibit5 * adv_sca_5 & + ibit4 * adv_sca_3 & ) * & ( sk(k,j+2,i) + sk(k,j-1,i) ) & + ( ibit5 * adv_sca_5 & ) * & ( sk(k,j+3,i) + sk(k,j-2,i) ) & ) diss_n = -ABS( v_comp ) * ( & ( 10.0_wp * ibit5 * adv_sca_5 & + 3.0_wp * ibit4 * adv_sca_3 & + ibit3 * adv_sca_1 & ) * & ( sk(k,j+1,i) - sk(k,j,i) ) & - ( 5.0_wp * ibit5 * adv_sca_5 & + ibit4 * adv_sca_3 & ) * & ( sk(k,j+2,i) - sk(k,j-1,i) ) & + ( ibit5 * adv_sca_5 & ) * & ( sk(k,j+3,i) - sk(k,j-2,i) ) & ) ! !-- indizes k_m, k_mm, ... should be known at these point ibit8 = IBITS(wall_flags_0(k-1,j,i),8,1) ibit7 = IBITS(wall_flags_0(k-1,j,i),7,1) ibit6 = IBITS(wall_flags_0(k-1,j,i),6,1) k_pp = k + 2 * ibit8 k_mm = k - 2 * ( ibit7 + ibit8 ) k_mmm = k - 3 * ibit8 flux_d = w(k-1,j,i) * ( & ( 37.0_wp * ibit8 * adv_sca_5 & + 7.0_wp * ibit7 * adv_sca_3 & + ibit6 * adv_sca_1 & ) * & ( sk(k,j,i) + sk(k-1,j,i) ) & - ( 8.0_wp * ibit8 * adv_sca_5 & + ibit7 * adv_sca_3 & ) * & ( sk(k+1,j,i) + sk(k_mm,j,i) ) & + ( ibit8 * adv_sca_5 & ) * ( sk(k_pp,j,i)+ sk(k_mmm,j,i) ) & ) diss_d = -ABS( w(k-1,j,i) ) * ( & ( 10.0_wp * ibit8 * adv_sca_5 & + 3.0_wp * ibit7 * adv_sca_3 & + ibit6 * adv_sca_1 & ) * & ( sk(k,j,i) - sk(k-1,j,i) ) & - ( 5.0_wp * ibit8 * adv_sca_5 & + ibit7 * adv_sca_3 & ) * & ( sk(k+1,j,i) - sk(k_mm,j,i) ) & + ( ibit8 * adv_sca_5 & ) * & ( sk(k_pp,j,i) - sk(k_mmm,j,i) ) & ) ibit8 = IBITS(wall_flags_0(k,j,i),8,1) ibit7 = IBITS(wall_flags_0(k,j,i),7,1) ibit6 = IBITS(wall_flags_0(k,j,i),6,1) k_ppp = k + 3 * ibit8 k_pp = k + 2 * ( 1 - ibit6 ) k_mm = k - 2 * ibit8 flux_t = w(k,j,i) * ( & ( 37.0_wp * ibit8 * adv_sca_5 & + 7.0_wp * ibit7 * adv_sca_3 & + ibit6 * adv_sca_1 & ) * & ( sk(k+1,j,i) + sk(k,j,i) ) & - ( 8.0_wp * ibit8 * adv_sca_5 & + ibit7 * adv_sca_3 & ) * & ( sk(k_pp,j,i) + sk(k-1,j,i) ) & + ( ibit8 * adv_sca_5 & ) * ( sk(k_ppp,j,i)+ sk(k_mm,j,i) ) & ) diss_t = -ABS( w(k,j,i) ) * ( & ( 10.0_wp * ibit8 * adv_sca_5 & + 3.0_wp * ibit7 * adv_sca_3 & + ibit6 * adv_sca_1 & ) * & ( sk(k+1,j,i) - sk(k,j,i) ) & - ( 5.0_wp * ibit8 * adv_sca_5 & + ibit7 * adv_sca_3 & ) * & ( sk(k_pp,j,i) - sk(k-1,j,i) ) & + ( ibit8 * adv_sca_5 & ) * & ( sk(k_ppp,j,i) - sk(k_mm,j,i) ) & ) ! !-- Calculate the divergence of the velocity field. A respective !-- correction is needed to overcome numerical instabilities caused !-- by a not sufficient reduction of divergences near topography. div = ( u(k,j,i+1) - u(k,j,i) ) * ddx & + ( v(k,j+1,i) - v(k,j,i) ) * ddy & + ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) tend(k,j,i) = - ( & ( flux_r + diss_r - flux_l - diss_l ) * ddx & + ( flux_n + diss_n - flux_s - diss_s ) * ddy & + ( flux_t + diss_t - flux_d - diss_d ) * ddzw(k)& ) + div * sk(k,j,i) !++ !-- Evaluation of statistics ! SELECT CASE ( sk_char ) ! ! CASE ( 'pt' ) ! sums_wspts_ws_l(k,tn) = sums_wspts_ws_l(k,tn) & ! + ( flux_t + diss_t ) & ! * weight_substep(intermediate_timestep_count) ! CASE ( 'sa' ) ! sums_wssas_ws_l(k,tn) = sums_wssas_ws_l(k,tn) & ! + ( flux_t + diss_t ) & ! * weight_substep(intermediate_timestep_count) ! CASE ( 'q' ) ! sums_wsqs_ws_l(k,tn) = sums_wsqs_ws_l(k,tn) & ! + ( flux_t + diss_t ) & ! * weight_substep(intermediate_timestep_count) ! CASE ( 'qr' ) ! sums_wsqrs_ws_l(k,tn) = sums_wsqrs_ws_l(k,tn) & ! + ( flux_t + diss_t ) & ! * weight_substep(intermediate_timestep_count) ! CASE ( 'nr' ) ! sums_wsnrs_ws_l(k,tn) = sums_wsnrs_ws_l(k,tn) & ! + ( flux_t + diss_t ) & ! * weight_substep(intermediate_timestep_count) ! ! END SELECT ENDDO ENDDO ENDDO !$acc end kernels END SUBROUTINE advec_s_ws_acc !------------------------------------------------------------------------------! ! Advection of u - Call for all grid points !------------------------------------------------------------------------------! SUBROUTINE advec_u_ws USE arrays_3d, & ONLY: ddzw, tend, u, v, w USE constants, & ONLY: adv_mom_1, adv_mom_3, adv_mom_5 USE control_parameters, & ONLY: intermediate_timestep_count, u_gtrans, v_gtrans USE grid_variables, & ONLY: ddx, ddy USE indices, & ONLY: nxl, nxlu, nxr, nyn, nys, nzb, nzb_max, nzt, wall_flags_0 USE kinds USE statistics, & ONLY: hom, sums_us2_ws_l, sums_wsus_ws_l, weight_substep IMPLICIT NONE INTEGER(iwp) :: i !: INTEGER(iwp) :: ibit9 !: INTEGER(iwp) :: ibit10 !: INTEGER(iwp) :: ibit11 !: INTEGER(iwp) :: ibit12 !: INTEGER(iwp) :: ibit13 !: INTEGER(iwp) :: ibit14 !: INTEGER(iwp) :: ibit15 !: INTEGER(iwp) :: ibit16 !: INTEGER(iwp) :: ibit17 !: INTEGER(iwp) :: j !: INTEGER(iwp) :: k !: INTEGER(iwp) :: k_mm !: INTEGER(iwp) :: k_pp !: INTEGER(iwp) :: k_ppp !: INTEGER(iwp) :: tn = 0 !: REAL(wp) :: diss_d !: REAL(wp) :: div !: REAL(wp) :: flux_d !: REAL(wp) :: gu !: REAL(wp) :: gv !: REAL(wp) :: v_comp !: REAL(wp) :: w_comp !: REAL(wp), DIMENSION(nzb+1:nzt) :: swap_diss_y_local_u !: REAL(wp), DIMENSION(nzb+1:nzt) :: swap_flux_y_local_u !: REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) :: swap_diss_x_local_u !: REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) :: swap_flux_x_local_u !: REAL(wp), DIMENSION(nzb:nzt) :: diss_n !: REAL(wp), DIMENSION(nzb:nzt) :: diss_r !: REAL(wp), DIMENSION(nzb:nzt) :: diss_t !: REAL(wp), DIMENSION(nzb:nzt) :: flux_n !: REAL(wp), DIMENSION(nzb:nzt) :: flux_r !: REAL(wp), DIMENSION(nzb:nzt) :: flux_t !: REAL(wp), DIMENSION(nzb:nzt) :: u_comp !: gu = 2.0_wp * u_gtrans gv = 2.0_wp * v_gtrans ! !-- Compute the fluxes for the whole left boundary of the processor domain. i = nxlu DO j = nys, nyn DO k = nzb+1, nzb_max ibit11 = IBITS(wall_flags_0(k,j,i),11,1) ibit10 = IBITS(wall_flags_0(k,j,i),10,1) ibit9 = IBITS(wall_flags_0(k,j,i),9,1) u_comp(k) = u(k,j,i) + u(k,j,i-1) - gu swap_flux_x_local_u(k,j) = u_comp(k) * ( & ( 37.0_wp * ibit11 * adv_mom_5 & + 7.0_wp * ibit10 * adv_mom_3 & + ibit9 * adv_mom_1 & ) * & ( u(k,j,i) + u(k,j,i-1) ) & - ( 8.0_wp * ibit11 * adv_mom_5 & + ibit10 * adv_mom_3 & ) * & ( u(k,j,i+1) + u(k,j,i-2) ) & + ( ibit11 * adv_mom_5 & ) * & ( u(k,j,i+2) + u(k,j,i-3) ) & ) swap_diss_x_local_u(k,j) = - ABS( u_comp(k) ) * ( & ( 10.0_wp * ibit11 * adv_mom_5 & + 3.0_wp * ibit10 * adv_mom_3 & + ibit9 * adv_mom_1 & ) * & ( u(k,j,i) - u(k,j,i-1) ) & - ( 5.0_wp * ibit11 * adv_mom_5 & + ibit10 * adv_mom_3 & ) * & ( u(k,j,i+1) - u(k,j,i-2) ) & + ( ibit11 * adv_mom_5 & ) * & ( u(k,j,i+2) - u(k,j,i-3) ) & ) ENDDO DO k = nzb_max+1, nzt u_comp(k) = u(k,j,i) + u(k,j,i-1) - gu swap_flux_x_local_u(k,j) = u_comp(k) * ( & 37.0_wp * ( u(k,j,i) + u(k,j,i-1) ) & - 8.0_wp * ( u(k,j,i+1) + u(k,j,i-2) ) & + ( u(k,j,i+2) + u(k,j,i-3) ) ) * adv_mom_5 swap_diss_x_local_u(k,j) = - ABS(u_comp(k)) * ( & 10.0_wp * ( u(k,j,i) - u(k,j,i-1) ) & - 5.0_wp * ( u(k,j,i+1) - u(k,j,i-2) ) & + ( u(k,j,i+2) - u(k,j,i-3) ) ) * adv_mom_5 ENDDO ENDDO DO i = nxlu, nxr ! !-- The following loop computes the fluxes for the south boundary points j = nys DO k = nzb+1, nzb_max ibit14 = IBITS(wall_flags_0(k,j,i),14,1) ibit13 = IBITS(wall_flags_0(k,j,i),13,1) ibit12 = IBITS(wall_flags_0(k,j,i),12,1) v_comp = v(k,j,i) + v(k,j,i-1) - gv swap_flux_y_local_u(k) = v_comp * ( & ( 37.0_wp * ibit14 * adv_mom_5 & + 7.0_wp * ibit13 * adv_mom_3 & + ibit12 * adv_mom_1 & ) * & ( u(k,j,i) + u(k,j-1,i) ) & - ( 8.0_wp * ibit14 * adv_mom_5 & + ibit13 * adv_mom_3 & ) * & ( u(k,j+1,i) + u(k,j-2,i) ) & + ( ibit14 * adv_mom_5 & ) * & ( u(k,j+2,i) + u(k,j-3,i) ) & ) swap_diss_y_local_u(k) = - ABS ( v_comp ) * ( & ( 10.0_wp * ibit14 * adv_mom_5 & + 3.0_wp * ibit13 * adv_mom_3 & + ibit12 * adv_mom_1 & ) * & ( u(k,j,i) - u(k,j-1,i) ) & - ( 5.0_wp * ibit14 * adv_mom_5 & + ibit13 * adv_mom_3 & ) * & ( u(k,j+1,i) - u(k,j-2,i) ) & + ( ibit14 * adv_mom_5 & ) * & ( u(k,j+2,i) - u(k,j-3,i) ) & ) ENDDO DO k = nzb_max+1, nzt v_comp = v(k,j,i) + v(k,j,i-1) - gv swap_flux_y_local_u(k) = v_comp * ( & 37.0_wp * ( u(k,j,i) + u(k,j-1,i) ) & - 8.0_wp * ( u(k,j+1,i) + u(k,j-2,i) ) & + ( u(k,j+2,i) + u(k,j-3,i) ) ) * adv_mom_5 swap_diss_y_local_u(k) = - ABS(v_comp) * ( & 10.0_wp * ( u(k,j,i) - u(k,j-1,i) ) & - 5.0_wp * ( u(k,j+1,i) - u(k,j-2,i) ) & + ( u(k,j+2,i) - u(k,j-3,i) ) ) * adv_mom_5 ENDDO ! !-- Computation of interior fluxes and tendency terms DO j = nys, nyn flux_t(0) = 0.0_wp diss_t(0) = 0.0_wp flux_d = 0.0_wp diss_d = 0.0_wp DO k = nzb+1, nzb_max ibit11 = IBITS(wall_flags_0(k,j,i),11,1) ibit10 = IBITS(wall_flags_0(k,j,i),10,1) ibit9 = IBITS(wall_flags_0(k,j,i),9,1) u_comp(k) = u(k,j,i+1) + u(k,j,i) flux_r(k) = ( u_comp(k) - gu ) * ( & ( 37.0_wp * ibit11 * adv_mom_5 & + 7.0_wp * ibit10 * adv_mom_3 & + ibit9 * adv_mom_1 & ) * & ( u(k,j,i+1) + u(k,j,i) ) & - ( 8.0_wp * ibit11 * adv_mom_5 & + ibit10 * adv_mom_3 & ) * & ( u(k,j,i+2) + u(k,j,i-1) ) & + ( ibit11 * adv_mom_5 & ) * & ( u(k,j,i+3) + u(k,j,i-2) ) & ) diss_r(k) = - ABS( u_comp(k) - gu ) * ( & ( 10.0_wp * ibit11 * adv_mom_5 & + 3.0_wp * ibit10 * adv_mom_3 & + ibit9 * adv_mom_1 & ) * & ( u(k,j,i+1) - u(k,j,i) ) & - ( 5.0_wp * ibit11 * adv_mom_5 & + ibit10 * adv_mom_3 & ) * & ( u(k,j,i+2) - u(k,j,i-1) ) & + ( ibit11 * adv_mom_5 & ) * & ( u(k,j,i+3) - u(k,j,i-2) ) & ) ibit14 = IBITS(wall_flags_0(k,j,i),14,1) ibit13 = IBITS(wall_flags_0(k,j,i),13,1) ibit12 = IBITS(wall_flags_0(k,j,i),12,1) v_comp = v(k,j+1,i) + v(k,j+1,i-1) - gv flux_n(k) = v_comp * ( & ( 37.0_wp * ibit14 * adv_mom_5 & + 7.0_wp * ibit13 * adv_mom_3 & + ibit12 * adv_mom_1 & ) * & ( u(k,j+1,i) + u(k,j,i) ) & - ( 8.0_wp * ibit14 * adv_mom_5 & + ibit13 * adv_mom_3 & ) * & ( u(k,j+2,i) + u(k,j-1,i) ) & + ( ibit14 * adv_mom_5 & ) * & ( u(k,j+3,i) + u(k,j-2,i) ) & ) diss_n(k) = - ABS ( v_comp ) * ( & ( 10.0_wp * ibit14 * adv_mom_5 & + 3.0_wp * ibit13 * adv_mom_3 & + ibit12 * adv_mom_1 & ) * & ( u(k,j+1,i) - u(k,j,i) ) & - ( 5.0_wp * ibit14 * adv_mom_5 & + ibit13 * adv_mom_3 & ) * & ( u(k,j+2,i) - u(k,j-1,i) ) & + ( ibit14 * adv_mom_5 & ) * & ( u(k,j+3,i) - u(k,j-2,i) ) & ) ! !-- k index has to be modified near bottom and top, else array !-- subscripts will be exceeded. ibit17 = IBITS(wall_flags_0(k,j,i),17,1) ibit16 = IBITS(wall_flags_0(k,j,i),16,1) ibit15 = IBITS(wall_flags_0(k,j,i),15,1) k_ppp = k + 3 * ibit17 k_pp = k + 2 * ( 1 - ibit15 ) k_mm = k - 2 * ibit17 w_comp = w(k,j,i) + w(k,j,i-1) flux_t(k) = w_comp * ( & ( 37.0_wp * ibit17 * adv_mom_5 & + 7.0_wp * ibit16 * adv_mom_3 & + ibit15 * adv_mom_1 & ) * & ( u(k+1,j,i) + u(k,j,i) ) & - ( 8.0_wp * ibit17 * adv_mom_5 & + ibit16 * adv_mom_3 & ) * & ( u(k_pp,j,i) + u(k-1,j,i) ) & + ( ibit17 * adv_mom_5 & ) * & ( u(k_ppp,j,i) + u(k_mm,j,i) ) & ) diss_t(k) = - ABS( w_comp ) * ( & ( 10.0_wp * ibit17 * adv_mom_5 & + 3.0_wp * ibit16 * adv_mom_3 & + ibit15 * adv_mom_1 & ) * & ( u(k+1,j,i) - u(k,j,i) ) & - ( 5.0_wp * ibit17 * adv_mom_5 & + ibit16 * adv_mom_3 & ) * & ( u(k_pp,j,i) - u(k-1,j,i) ) & + ( ibit17 * adv_mom_5 & ) * & ( u(k_ppp,j,i) - u(k_mm,j,i) ) & ) ! !-- Calculate the divergence of the velocity field. A respective !-- correction is needed to overcome numerical instabilities caused !-- by a not sufficient reduction of divergences near topography. div = ( ( u_comp(k) - ( u(k,j,i) + u(k,j,i-1) ) ) * ddx & + ( v_comp + gv - ( v(k,j,i) + v(k,j,i-1 ) ) ) * ddy & + ( w_comp - ( w(k-1,j,i) + w(k-1,j,i-1) ) ) & * ddzw(k) & ) * 0.5_wp tend(k,j,i) = tend(k,j,i) - ( & ( flux_r(k) + diss_r(k) & - swap_flux_x_local_u(k,j) - swap_diss_x_local_u(k,j) ) * ddx & + ( flux_n(k) + diss_n(k) & - swap_flux_y_local_u(k) - swap_diss_y_local_u(k) ) * ddy & + ( flux_t(k) + diss_t(k) & - flux_d - diss_d & ) * ddzw(k) & ) + div * u(k,j,i) swap_flux_x_local_u(k,j) = flux_r(k) swap_diss_x_local_u(k,j) = diss_r(k) swap_flux_y_local_u(k) = flux_n(k) swap_diss_y_local_u(k) = diss_n(k) flux_d = flux_t(k) diss_d = diss_t(k) ! !-- Statistical Evaluation of u'u'. The factor has to be applied !-- for right evaluation when gallilei_trans = .T. . sums_us2_ws_l(k,tn) = sums_us2_ws_l(k,tn) & + ( flux_r(k) * & ( u_comp(k) - 2.0_wp * hom(k,1,1,0) ) & / ( u_comp(k) - gu + 1.0E-20_wp ) & + diss_r(k) * & ABS( u_comp(k) - 2.0_wp * hom(k,1,1,0) ) & / ( ABS( u_comp(k) - gu ) + 1.0E-20_wp ) ) & * weight_substep(intermediate_timestep_count) ! !-- Statistical Evaluation of w'u'. sums_wsus_ws_l(k,tn) = sums_wsus_ws_l(k,tn) & + ( flux_t(k) + diss_t(k) ) & * weight_substep(intermediate_timestep_count) ENDDO DO k = nzb_max+1, nzt u_comp(k) = u(k,j,i+1) + u(k,j,i) flux_r(k) = ( u_comp(k) - gu ) * ( & 37.0_wp * ( u(k,j,i+1) + u(k,j,i) ) & - 8.0_wp * ( u(k,j,i+2) + u(k,j,i-1) ) & + ( u(k,j,i+3) + u(k,j,i-2) ) ) * adv_mom_5 diss_r(k) = - ABS( u_comp(k) - gu ) * ( & 10.0_wp * ( u(k,j,i+1) - u(k,j,i) ) & - 5.0_wp * ( u(k,j,i+2) - u(k,j,i-1) ) & + ( u(k,j,i+3) - u(k,j,i-2) ) ) * adv_mom_5 v_comp = v(k,j+1,i) + v(k,j+1,i-1) - gv flux_n(k) = v_comp * ( & 37.0_wp * ( u(k,j+1,i) + u(k,j,i) ) & - 8.0_wp * ( u(k,j+2,i) + u(k,j-1,i) ) & + ( u(k,j+3,i) + u(k,j-2,i) ) ) * adv_mom_5 diss_n(k) = - ABS( v_comp ) * ( & 10.0_wp * ( u(k,j+1,i) - u(k,j,i) ) & - 5.0_wp * ( u(k,j+2,i) - u(k,j-1,i) ) & + ( u(k,j+3,i) - u(k,j-2,i) ) ) * adv_mom_5 ! !-- k index has to be modified near bottom and top, else array !-- subscripts will be exceeded. ibit17 = IBITS(wall_flags_0(k,j,i),17,1) ibit16 = IBITS(wall_flags_0(k,j,i),16,1) ibit15 = IBITS(wall_flags_0(k,j,i),15,1) k_ppp = k + 3 * ibit17 k_pp = k + 2 * ( 1 - ibit15 ) k_mm = k - 2 * ibit17 w_comp = w(k,j,i) + w(k,j,i-1) flux_t(k) = w_comp * ( & ( 37.0_wp * ibit17 * adv_mom_5 & + 7.0_wp * ibit16 * adv_mom_3 & + ibit15 * adv_mom_1 & ) * & ( u(k+1,j,i) + u(k,j,i) ) & - ( 8.0_wp * ibit17 * adv_mom_5 & + ibit16 * adv_mom_3 & ) * & ( u(k_pp,j,i) + u(k-1,j,i) ) & + ( ibit17 * adv_mom_5 & ) * & ( u(k_ppp,j,i) + u(k_mm,j,i) ) & ) diss_t(k) = - ABS( w_comp ) * ( & ( 10.0_wp * ibit17 * adv_mom_5 & + 3.0_wp * ibit16 * adv_mom_3 & + ibit15 * adv_mom_1 & ) * & ( u(k+1,j,i) - u(k,j,i) ) & - ( 5.0_wp * ibit17 * adv_mom_5 & + ibit16 * adv_mom_3 & ) * & ( u(k_pp,j,i) - u(k-1,j,i) ) & + ( ibit17 * adv_mom_5 & ) * & ( u(k_ppp,j,i) - u(k_mm,j,i) ) & ) ! !-- Calculate the divergence of the velocity field. A respective !-- correction is needed to overcome numerical instabilities caused !-- by a not sufficient reduction of divergences near topography. div = ( ( u_comp(k) - ( u(k,j,i) + u(k,j,i-1) ) ) * ddx & + ( v_comp + gv - ( v(k,j,i) + v(k,j,i-1 ) ) ) * ddy & + ( w_comp - ( w(k-1,j,i) + w(k-1,j,i-1) ) ) & * ddzw(k) & ) * 0.5_wp tend(k,j,i) = tend(k,j,i) - ( & ( flux_r(k) + diss_r(k) & - swap_flux_x_local_u(k,j) - swap_diss_x_local_u(k,j) ) * ddx & + ( flux_n(k) + diss_n(k) & - swap_flux_y_local_u(k) - swap_diss_y_local_u(k) ) * ddy & + ( flux_t(k) + diss_t(k) & - flux_d - diss_d & ) * ddzw(k) & ) + div * u(k,j,i) swap_flux_x_local_u(k,j) = flux_r(k) swap_diss_x_local_u(k,j) = diss_r(k) swap_flux_y_local_u(k) = flux_n(k) swap_diss_y_local_u(k) = diss_n(k) flux_d = flux_t(k) diss_d = diss_t(k) ! !-- Statistical Evaluation of u'u'. The factor has to be applied !-- for right evaluation when gallilei_trans = .T. . sums_us2_ws_l(k,tn) = sums_us2_ws_l(k,tn) & + ( flux_r(k) * & ( u_comp(k) - 2.0_wp * hom(k,1,1,0) ) & / ( u_comp(k) - gu + 1.0E-20_wp ) & + diss_r(k) * & ABS( u_comp(k) - 2.0_wp * hom(k,1,1,0) ) & / ( ABS( u_comp(k) - gu ) + 1.0E-20_wp ) ) & * weight_substep(intermediate_timestep_count) ! !-- Statistical Evaluation of w'u'. sums_wsus_ws_l(k,tn) = sums_wsus_ws_l(k,tn) & + ( flux_t(k) + diss_t(k) ) & * weight_substep(intermediate_timestep_count) ENDDO ENDDO ENDDO sums_us2_ws_l(nzb,tn) = sums_us2_ws_l(nzb+1,tn) END SUBROUTINE advec_u_ws !------------------------------------------------------------------------------! ! Advection of u - Call for all grid points - accelerator version !------------------------------------------------------------------------------! SUBROUTINE advec_u_ws_acc USE arrays_3d, & ONLY: ddzw, tend, u, v, w USE constants, & ONLY: adv_mom_1, adv_mom_3, adv_mom_5 USE control_parameters, & ONLY: intermediate_timestep_count, u_gtrans, v_gtrans USE grid_variables, & ONLY: ddx, ddy USE indices, & ONLY: i_left, i_right, j_north, j_south, nxl, nxr, nyn, nys, nzb, & nzb_max, nzt, wall_flags_0 USE kinds ! USE statistics, & ! ONLY: hom, sums_us2_ws_l, sums_wsus_ws_l, weight_substep IMPLICIT NONE INTEGER(iwp) :: i !: INTEGER(iwp) :: ibit9 !: INTEGER(iwp) :: ibit10 !: INTEGER(iwp) :: ibit11 !: INTEGER(iwp) :: ibit12 !: INTEGER(iwp) :: ibit13 !: INTEGER(iwp) :: ibit14 !: INTEGER(iwp) :: ibit15 !: INTEGER(iwp) :: ibit16 !: INTEGER(iwp) :: ibit17 !: INTEGER(iwp) :: j !: INTEGER(iwp) :: k !: INTEGER(iwp) :: k_mmm !: INTEGER(iwp) :: k_mm !: INTEGER(iwp) :: k_pp !: INTEGER(iwp) :: k_ppp !: INTEGER(iwp) :: tn = 0 !: REAL(wp) :: diss_d !: REAL(wp) :: diss_l !: REAL(wp) :: diss_n !: REAL(wp) :: diss_r !: REAL(wp) :: diss_s !: REAL(wp) :: diss_t !: REAL(wp) :: div !: REAL(wp) :: flux_d !: REAL(wp) :: flux_l !: REAL(wp) :: flux_n !: REAL(wp) :: flux_r !: REAL(wp) :: flux_s !: REAL(wp) :: flux_t !: REAL(wp) :: gu !: REAL(wp) :: gv !: REAL(wp) :: u_comp !: REAL(wp) :: u_comp_l !: REAL(wp) :: v_comp !: REAL(wp) :: v_comp_s !: REAL(wp) :: w_comp !: gu = 2.0_wp * u_gtrans gv = 2.0_wp * v_gtrans ! !-- Computation of fluxes and tendency terms !$acc kernels present( ddzw, tend, u, v, w, wall_flags_0 ) DO i = i_left, i_right DO j = j_south, j_north DO k = nzb+1, nzt ibit11 = IBITS(wall_flags_0(k,j,i),11,1) ibit10 = IBITS(wall_flags_0(k,j,i),10,1) ibit9 = IBITS(wall_flags_0(k,j,i),9,1) u_comp_l = u(k,j,i) + u(k,j,i-1) - gu flux_l = u_comp_l * ( & ( 37.0_wp * ibit11 * adv_mom_5 & + 7.0_wp * ibit10 * adv_mom_3 & + ibit9 * adv_mom_1 & ) * & ( u(k,j,i) + u(k,j,i-1) ) & - ( 8.0_wp * ibit11 * adv_mom_5 & + ibit10 * adv_mom_3 & ) * & ( u(k,j,i+1) + u(k,j,i-2) ) & + ( ibit11 * adv_mom_5 & ) * & ( u(k,j,i+2) + u(k,j,i-3) ) & ) diss_l = - ABS( u_comp_l ) * ( & ( 10.0_wp * ibit11 * adv_mom_5 & + 3.0_wp * ibit10 * adv_mom_3 & + ibit9 * adv_mom_1 & ) * & ( u(k,j,i) - u(k,j,i-1) ) & - ( 5.0_wp * ibit11 * adv_mom_5 & + ibit10 * adv_mom_3 & ) * & ( u(k,j,i+1) - u(k,j,i-2) ) & + ( ibit11 * adv_mom_5 & ) * & ( u(k,j,i+2) - u(k,j,i-3) ) & ) u_comp = u(k,j,i+1) + u(k,j,i) flux_r = ( u_comp - gu ) * ( & ( 37.0_wp * ibit11 * adv_mom_5 & + 7.0_wp * ibit10 * adv_mom_3 & + ibit9 * adv_mom_1 & ) * & ( u(k,j,i+1) + u(k,j,i) ) & - ( 8.0_wp * ibit11 * adv_mom_5 & + ibit10 * adv_mom_3 & ) * & ( u(k,j,i+2) + u(k,j,i-1) ) & + ( ibit11 * adv_mom_5 & ) * & ( u(k,j,i+3) + u(k,j,i-2) ) & ) diss_r = - ABS( u_comp - gu ) * ( & ( 10.0_wp * ibit11 * adv_mom_5 & + 3.0_wp * ibit10 * adv_mom_3 & + ibit9 * adv_mom_1 & ) * & ( u(k,j,i+1) - u(k,j,i) ) & - ( 5.0_wp * ibit11 * adv_mom_5 & + ibit10 * adv_mom_3 & ) * & ( u(k,j,i+2) - u(k,j,i-1) ) & + ( ibit11 * adv_mom_5 & ) * & ( u(k,j,i+3) - u(k,j,i-2) ) & ) ibit14 = IBITS(wall_flags_0(k,j,i),14,1) ibit13 = IBITS(wall_flags_0(k,j,i),13,1) ibit12 = IBITS(wall_flags_0(k,j,i),12,1) v_comp_s = v(k,j,i) + v(k,j,i-1) - gv flux_s = v_comp_s * ( & ( 37.0_wp * ibit14 * adv_mom_5 & + 7.0_wp * ibit13 * adv_mom_3 & + ibit12 * adv_mom_1 & ) * & ( u(k,j,i) + u(k,j-1,i) ) & - ( 8.0_wp * ibit14 * adv_mom_5 & + ibit13 * adv_mom_3 & ) * & ( u(k,j+1,i) + u(k,j-2,i) ) & + ( ibit14 * adv_mom_5 & ) * & ( u(k,j+2,i) + u(k,j-3,i) ) & ) diss_s = - ABS ( v_comp_s ) * ( & ( 10.0_wp * ibit14 * adv_mom_5 & + 3.0_wp * ibit13 * adv_mom_3 & + ibit12 * adv_mom_1 & ) * & ( u(k,j,i) - u(k,j-1,i) ) & - ( 5.0_wp * ibit14 * adv_mom_5 & + ibit13 * adv_mom_3 & ) * & ( u(k,j+1,i) - u(k,j-2,i) ) & + ( ibit14 * adv_mom_5 & ) * & ( u(k,j+2,i) - u(k,j-3,i) ) & ) v_comp = v(k,j+1,i) + v(k,j+1,i-1) - gv flux_n = v_comp * ( & ( 37.0_wp * ibit14 * adv_mom_5 & + 7.0_wp * ibit13 * adv_mom_3 & + ibit12 * adv_mom_1 & ) * & ( u(k,j+1,i) + u(k,j,i) ) & - ( 8.0_wp * ibit14 * adv_mom_5 & + ibit13 * adv_mom_3 & ) * & ( u(k,j+2,i) + u(k,j-1,i) ) & + ( ibit14 * adv_mom_5 & ) * & ( u(k,j+3,i) + u(k,j-2,i) ) & ) diss_n = - ABS ( v_comp ) * ( & ( 10.0_wp * ibit14 * adv_mom_5 & + 3.0_wp * ibit13 * adv_mom_3 & + ibit12 * adv_mom_1 & ) * & ( u(k,j+1,i) - u(k,j,i) ) & - ( 5.0_wp * ibit14 * adv_mom_5 & + ibit13 * adv_mom_3 & ) * & ( u(k,j+2,i) - u(k,j-1,i) ) & + ( ibit14 * adv_mom_5 & ) * & ( u(k,j+3,i) - u(k,j-2,i) ) & ) ibit17 = IBITS(wall_flags_0(k-1,j,i),17,1) ibit16 = IBITS(wall_flags_0(k-1,j,i),16,1) ibit15 = IBITS(wall_flags_0(k-1,j,i),15,1) k_pp = k + 2 * ibit17 k_mm = k - 2 * ( ibit16 + ibit17 ) k_mmm = k - 3 * ibit17 w_comp = w(k-1,j,i) + w(k-1,j,i-1) flux_d = w_comp * ( & ( 37.0_wp * ibit17 * adv_mom_5 & + 7.0_wp * ibit16 * adv_mom_3 & + ibit15 * adv_mom_1 & ) * & ( u(k,j,i) + u(k-1,j,i) ) & - ( 8.0_wp * ibit17 * adv_mom_5 & + ibit16 * adv_mom_3 & ) * & ( u(k+1,j,i) + u(k_mm,j,i) ) & + ( ibit17 * adv_mom_5 & ) * & ( u(k_pp,j,i) + u(k_mmm,j,i) ) & ) diss_d = - ABS( w_comp ) * ( & ( 10.0_wp * ibit17 * adv_mom_5 & + 3.0_wp * ibit16 * adv_mom_3 & + ibit15 * adv_mom_1 & ) * & ( u(k,j,i) - u(k-1,j,i) ) & - ( 5.0_wp * ibit17 * adv_mom_5 & + ibit16 * adv_mom_3 & ) * & ( u(k+1,j,i) - u(k_mm,j,i) ) & + ( ibit17 * adv_mom_5 & ) * & ( u(k_pp,j,i) - u(k_mmm,j,i) ) & ) ! !-- k index has to be modified near bottom and top, else array !-- subscripts will be exceeded. ibit17 = IBITS(wall_flags_0(k,j,i),17,1) ibit16 = IBITS(wall_flags_0(k,j,i),16,1) ibit15 = IBITS(wall_flags_0(k,j,i),15,1) k_ppp = k + 3 * ibit17 k_pp = k + 2 * ( 1 - ibit15 ) k_mm = k - 2 * ibit17 w_comp = w(k,j,i) + w(k,j,i-1) flux_t = w_comp * ( & ( 37.0_wp * ibit17 * adv_mom_5 & + 7.0_wp * ibit16 * adv_mom_3 & + ibit15 * adv_mom_1 & ) * & ( u(k+1,j,i) + u(k,j,i) ) & - ( 8.0_wp * ibit17 * adv_mom_5 & + ibit16 * adv_mom_3 & ) * & ( u(k_pp,j,i) + u(k-1,j,i) ) & + ( ibit17 * adv_mom_5 & ) * & ( u(k_ppp,j,i) + u(k_mm,j,i) ) & ) diss_t = - ABS( w_comp ) * ( & ( 10.0_wp * ibit17 * adv_mom_5 & + 3.0_wp * ibit16 * adv_mom_3 & + ibit15 * adv_mom_1 & ) * & ( u(k+1,j,i) - u(k,j,i) ) & - ( 5.0_wp * ibit17 * adv_mom_5 & + ibit16 * adv_mom_3 & ) * & ( u(k_pp,j,i) - u(k-1,j,i) ) & + ( ibit17 * adv_mom_5 & ) * & ( u(k_ppp,j,i) - u(k_mm,j,i) ) & ) ! !-- Calculate the divergence of the velocity field. A respective !-- correction is needed to overcome numerical instabilities caused !-- by a not sufficient reduction of divergences near topography. div = ( ( u_comp - ( u(k,j,i) + u(k,j,i-1) ) ) * ddx & + ( v_comp + gv - ( v(k,j,i) + v(k,j,i-1 ) ) ) * ddy & + ( w_comp - ( w(k-1,j,i) + w(k-1,j,i-1) ) ) & * ddzw(k) & ) * 0.5_wp tend(k,j,i) = - ( & ( flux_r + diss_r - flux_l - diss_l ) * ddx & + ( flux_n + diss_n - flux_s - diss_s ) * ddy & + ( flux_t + diss_t - flux_d - diss_d ) * ddzw(k) & ) + div * u(k,j,i) !++ !-- Statistical Evaluation of u'u'. The factor has to be applied !-- for right evaluation when gallilei_trans = .T. . ! sums_us2_ws_l(k,tn) = sums_us2_ws_l(k,tn) & ! + ( flux_r * & ! ( u_comp - 2.0_wp * hom(k,1,1,0) ) & ! / ( u_comp - gu + 1.0E-20_wp ) & ! + diss_r * & ! ABS( u_comp - 2.0_wp * hom(k,1,1,0) ) & ! / ( ABS( u_comp - gu ) + 1.0E-20_wp ) ) & ! * weight_substep(intermediate_timestep_count) ! !-- Statistical Evaluation of w'u'. ! sums_wsus_ws_l(k,tn) = sums_wsus_ws_l(k,tn) & ! + ( flux_t + diss_t ) & ! * weight_substep(intermediate_timestep_count) ENDDO ENDDO ENDDO !$acc end kernels !++ ! sums_us2_ws_l(nzb,tn) = sums_us2_ws_l(nzb+1,tn) END SUBROUTINE advec_u_ws_acc !------------------------------------------------------------------------------! ! Advection of v - Call for all grid points !------------------------------------------------------------------------------! SUBROUTINE advec_v_ws USE arrays_3d, & ONLY: ddzw, tend, u, v, w USE constants, & ONLY: adv_mom_1, adv_mom_3, adv_mom_5 USE control_parameters, & ONLY: intermediate_timestep_count, u_gtrans, v_gtrans USE grid_variables, & ONLY: ddx, ddy USE indices, & ONLY: nxl, nxr, nyn, nys, nysv, nzb, nzb_max, nzt, wall_flags_0 USE kinds USE statistics, & ONLY: hom, sums_vs2_ws_l, sums_wsvs_ws_l, weight_substep IMPLICIT NONE INTEGER(iwp) :: i !: INTEGER(iwp) :: ibit18 !: INTEGER(iwp) :: ibit19 !: INTEGER(iwp) :: ibit20 !: INTEGER(iwp) :: ibit21 !: INTEGER(iwp) :: ibit22 !: INTEGER(iwp) :: ibit23 !: INTEGER(iwp) :: ibit24 !: INTEGER(iwp) :: ibit25 !: INTEGER(iwp) :: ibit26 !: INTEGER(iwp) :: j !: INTEGER(iwp) :: k !: INTEGER(iwp) :: k_mm !: INTEGER(iwp) :: k_pp !: INTEGER(iwp) :: k_ppp !: INTEGER(iwp) :: tn = 0 !: REAL(wp) :: diss_d !: REAL(wp) :: div !: REAL(wp) :: flux_d !: REAL(wp) :: gu !: REAL(wp) :: gv !: REAL(wp) :: u_comp !: REAL(wp) :: w_comp !: REAL(wp), DIMENSION(nzb+1:nzt) :: swap_diss_y_local_v !: REAL(wp), DIMENSION(nzb+1:nzt) :: swap_flux_y_local_v !: REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) :: swap_diss_x_local_v !: REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) :: swap_flux_x_local_v !: REAL(wp), DIMENSION(nzb:nzt) :: diss_n !: REAL(wp), DIMENSION(nzb:nzt) :: diss_r !: REAL(wp), DIMENSION(nzb:nzt) :: diss_t !: REAL(wp), DIMENSION(nzb:nzt) :: flux_n !: REAL(wp), DIMENSION(nzb:nzt) :: flux_r !: REAL(wp), DIMENSION(nzb:nzt) :: flux_t !: REAL(wp), DIMENSION(nzb:nzt) :: v_comp !: gu = 2.0_wp * u_gtrans gv = 2.0_wp * v_gtrans ! !-- First compute the whole left boundary of the processor domain i = nxl DO j = nysv, nyn DO k = nzb+1, nzb_max ibit20 = IBITS(wall_flags_0(k,j,i),20,1) ibit19 = IBITS(wall_flags_0(k,j,i),19,1) ibit18 = IBITS(wall_flags_0(k,j,i),18,1) u_comp = u(k,j-1,i) + u(k,j,i) - gu swap_flux_x_local_v(k,j) = u_comp * ( & ( 37.0_wp * ibit20 * adv_mom_5 & + 7.0_wp * ibit19 * adv_mom_3 & + ibit18 * adv_mom_1 & ) * & ( v(k,j,i) + v(k,j,i-1) ) & - ( 8.0_wp * ibit20 * adv_mom_5 & + ibit19 * adv_mom_3 & ) * & ( v(k,j,i+1) + v(k,j,i-2) ) & + ( ibit20 * adv_mom_5 & ) * & ( v(k,j,i+2) + v(k,j,i-3) ) & ) swap_diss_x_local_v(k,j) = - ABS( u_comp ) * ( & ( 10.0_wp * ibit20 * adv_mom_5 & + 3.0_wp * ibit19 * adv_mom_3 & + ibit18 * adv_mom_1 & ) * & ( v(k,j,i) - v(k,j,i-1) ) & - ( 5.0_wp * ibit20 * adv_mom_5 & + ibit19 * adv_mom_3 & ) * & ( v(k,j,i+1) - v(k,j,i-2) ) & + ( ibit20 * adv_mom_5 & ) * & ( v(k,j,i+2) - v(k,j,i-3) ) & ) ENDDO DO k = nzb_max+1, nzt u_comp = u(k,j-1,i) + u(k,j,i) - gu swap_flux_x_local_v(k,j) = u_comp * ( & 37.0_wp * ( v(k,j,i) + v(k,j,i-1) ) & - 8.0_wp * ( v(k,j,i+1) + v(k,j,i-2) ) & + ( v(k,j,i+2) + v(k,j,i-3) ) ) * adv_mom_5 swap_diss_x_local_v(k,j) = - ABS( u_comp ) * ( & 10.0_wp * ( v(k,j,i) - v(k,j,i-1) ) & - 5.0_wp * ( v(k,j,i+1) - v(k,j,i-2) ) & + ( v(k,j,i+2) - v(k,j,i-3) ) ) * adv_mom_5 ENDDO ENDDO DO i = nxl, nxr j = nysv DO k = nzb+1, nzb_max ibit23 = IBITS(wall_flags_0(k,j,i),23,1) ibit22 = IBITS(wall_flags_0(k,j,i),22,1) ibit21 = IBITS(wall_flags_0(k,j,i),21,1) v_comp(k) = v(k,j,i) + v(k,j-1,i) - gv swap_flux_y_local_v(k) = v_comp(k) * ( & ( 37.0_wp * ibit23 * adv_mom_5 & + 7.0_wp * ibit22 * adv_mom_3 & + ibit21 * adv_mom_1 & ) * & ( v(k,j,i) + v(k,j-1,i) ) & - ( 8.0_wp * ibit23 * adv_mom_5 & + ibit22 * adv_mom_3 & ) * & ( v(k,j+1,i) + v(k,j-2,i) ) & + ( ibit23 * adv_mom_5 & ) * & ( v(k,j+2,i) + v(k,j-3,i) ) & ) swap_diss_y_local_v(k) = - ABS( v_comp(k) ) * ( & ( 10.0_wp * ibit23 * adv_mom_5 & + 3.0_wp * ibit22 * adv_mom_3 & + ibit21 * adv_mom_1 & ) * & ( v(k,j,i) - v(k,j-1,i) ) & - ( 5.0_wp * ibit23 * adv_mom_5 & + ibit22 * adv_mom_3 & ) * & ( v(k,j+1,i) - v(k,j-2,i) ) & + ( ibit23 * adv_mom_5 & ) * & ( v(k,j+2,i) - v(k,j-3,i) ) & ) ENDDO DO k = nzb_max+1, nzt v_comp(k) = v(k,j,i) + v(k,j-1,i) - gv swap_flux_y_local_v(k) = v_comp(k) * ( & 37.0_wp * ( v(k,j,i) + v(k,j-1,i) ) & - 8.0_wp * ( v(k,j+1,i) + v(k,j-2,i) ) & + ( v(k,j+2,i) + v(k,j-3,i) ) ) * adv_mom_5 swap_diss_y_local_v(k) = - ABS( v_comp(k) ) * ( & 10.0_wp * ( v(k,j,i) - v(k,j-1,i) ) & - 5.0_wp * ( v(k,j+1,i) - v(k,j-2,i) ) & + ( v(k,j+2,i) - v(k,j-3,i) ) ) * adv_mom_5 ENDDO DO j = nysv, nyn flux_t(0) = 0.0_wp diss_t(0) = 0.0_wp flux_d = 0.0_wp diss_d = 0.0_wp DO k = nzb+1, nzb_max ibit20 = IBITS(wall_flags_0(k,j,i),20,1) ibit19 = IBITS(wall_flags_0(k,j,i),19,1) ibit18 = IBITS(wall_flags_0(k,j,i),18,1) u_comp = u(k,j-1,i+1) + u(k,j,i+1) - gu flux_r(k) = u_comp * ( & ( 37.0_wp * ibit20 * adv_mom_5 & + 7.0_wp * ibit19 * adv_mom_3 & + ibit18 * adv_mom_1 & ) * & ( v(k,j,i+1) + v(k,j,i) ) & - ( 8.0_wp * ibit20 * adv_mom_5 & + ibit19 * adv_mom_3 & ) * & ( v(k,j,i+2) + v(k,j,i-1) ) & + ( ibit20 * adv_mom_5 & ) * & ( v(k,j,i+3) + v(k,j,i-2) ) & ) diss_r(k) = - ABS( u_comp ) * ( & ( 10.0_wp * ibit20 * adv_mom_5 & + 3.0_wp * ibit19 * adv_mom_3 & + ibit18 * adv_mom_1 & ) * & ( v(k,j,i+1) - v(k,j,i) ) & - ( 5.0_wp * ibit20 * adv_mom_5 & + ibit19 * adv_mom_3 & ) * & ( v(k,j,i+2) - v(k,j,i-1) ) & + ( ibit20 * adv_mom_5 & ) * & ( v(k,j,i+3) - v(k,j,i-2) ) & ) ibit23 = IBITS(wall_flags_0(k,j,i),23,1) ibit22 = IBITS(wall_flags_0(k,j,i),22,1) ibit21 = IBITS(wall_flags_0(k,j,i),21,1) v_comp(k) = v(k,j+1,i) + v(k,j,i) flux_n(k) = ( v_comp(k) - gv ) * ( & ( 37.0_wp * ibit23 * adv_mom_5 & + 7.0_wp * ibit22 * adv_mom_3 & + ibit21 * adv_mom_1 & ) * & ( v(k,j+1,i) + v(k,j,i) ) & - ( 8.0_wp * ibit23 * adv_mom_5 & + ibit22 * adv_mom_3 & ) * & ( v(k,j+2,i) + v(k,j-1,i) ) & + ( ibit23 * adv_mom_5 & ) * & ( v(k,j+3,i) + v(k,j-2,i) ) & ) diss_n(k) = - ABS( v_comp(k) - gv ) * ( & ( 10.0_wp * ibit23 * adv_mom_5 & + 3.0_wp * ibit22 * adv_mom_3 & + ibit21 * adv_mom_1 & ) * & ( v(k,j+1,i) - v(k,j,i) ) & - ( 5.0_wp * ibit23 * adv_mom_5 & + ibit22 * adv_mom_3 & ) * & ( v(k,j+2,i) - v(k,j-1,i) ) & + ( ibit23 * adv_mom_5 & ) * & ( v(k,j+3,i) - v(k,j-2,i) ) & ) ! !-- k index has to be modified near bottom and top, else array !-- subscripts will be exceeded. ibit26 = IBITS(wall_flags_0(k,j,i),26,1) ibit25 = IBITS(wall_flags_0(k,j,i),25,1) ibit24 = IBITS(wall_flags_0(k,j,i),24,1) k_ppp = k + 3 * ibit26 k_pp = k + 2 * ( 1 - ibit24 ) k_mm = k - 2 * ibit26 w_comp = w(k,j-1,i) + w(k,j,i) flux_t(k) = w_comp * ( & ( 37.0_wp * ibit26 * adv_mom_5 & + 7.0_wp * ibit25 * adv_mom_3 & + ibit24 * adv_mom_1 & ) * & ( v(k+1,j,i) + v(k,j,i) ) & - ( 8.0_wp * ibit26 * adv_mom_5 & + ibit25 * adv_mom_3 & ) * & ( v(k_pp,j,i) + v(k-1,j,i) ) & + ( ibit26 * adv_mom_5 & ) * & ( v(k_ppp,j,i) + v(k_mm,j,i) ) & ) diss_t(k) = - ABS( w_comp ) * ( & ( 10.0_wp * ibit26 * adv_mom_5 & + 3.0_wp * ibit25 * adv_mom_3 & + ibit24 * adv_mom_1 & ) * & ( v(k+1,j,i) - v(k,j,i) ) & - ( 5.0_wp * ibit26 * adv_mom_5 & + ibit25 * adv_mom_3 & ) * & ( v(k_pp,j,i) - v(k-1,j,i) ) & + ( ibit26 * adv_mom_5 & ) * & ( v(k_ppp,j,i) - v(k_mm,j,i) ) & ) ! !-- Calculate the divergence of the velocity field. A respective !-- correction is needed to overcome numerical instabilities caused !-- by a not sufficient reduction of divergences near topography. div = ( ( u_comp + gu - ( u(k,j-1,i) + u(k,j,i) ) ) * ddx & + ( v_comp(k) - ( v(k,j,i) + v(k,j-1,i) ) ) * ddy & + ( w_comp - ( w(k-1,j-1,i) + w(k-1,j,i) ) & ) * ddzw(k) & ) * 0.5_wp tend(k,j,i) = tend(k,j,i) - ( & ( flux_r(k) + diss_r(k) & - swap_flux_x_local_v(k,j) - swap_diss_x_local_v(k,j) & ) * ddx & + ( flux_n(k) + diss_n(k) & - swap_flux_y_local_v(k) - swap_diss_y_local_v(k) & ) * ddy & + ( flux_t(k) + diss_t(k) & - flux_d - diss_d & ) * ddzw(k) & ) + v(k,j,i) * div swap_flux_x_local_v(k,j) = flux_r(k) swap_diss_x_local_v(k,j) = diss_r(k) swap_flux_y_local_v(k) = flux_n(k) swap_diss_y_local_v(k) = diss_n(k) flux_d = flux_t(k) diss_d = diss_t(k) ! !-- Statistical Evaluation of v'v'. The factor has to be applied !-- for right evaluation when gallilei_trans = .T. . sums_vs2_ws_l(k,tn) = sums_vs2_ws_l(k,tn) & + ( flux_n(k) & * ( v_comp(k) - 2.0_wp * hom(k,1,2,0) ) & / ( v_comp(k) - gv + 1.0E-20_wp ) & + diss_n(k) & * ABS( v_comp(k) - 2.0_wp * hom(k,1,2,0) ) & / ( ABS( v_comp(k) - gv ) +1.0E-20_wp ) ) & * weight_substep(intermediate_timestep_count) ! !-- Statistical Evaluation of w'v'. sums_wsvs_ws_l(k,tn) = sums_wsvs_ws_l(k,tn) & + ( flux_t(k) + diss_t(k) ) & * weight_substep(intermediate_timestep_count) ENDDO DO k = nzb_max+1, nzt u_comp = u(k,j-1,i+1) + u(k,j,i+1) - gu flux_r(k) = u_comp * ( & 37.0_wp * ( v(k,j,i+1) + v(k,j,i) ) & - 8.0_wp * ( v(k,j,i+2) + v(k,j,i-1) ) & + ( v(k,j,i+3) + v(k,j,i-2) ) ) * adv_mom_5 diss_r(k) = - ABS( u_comp ) * ( & 10.0_wp * ( v(k,j,i+1) - v(k,j,i) ) & - 5.0_wp * ( v(k,j,i+2) - v(k,j,i-1) ) & + ( v(k,j,i+3) - v(k,j,i-2) ) ) * adv_mom_5 v_comp(k) = v(k,j+1,i) + v(k,j,i) flux_n(k) = ( v_comp(k) - gv ) * ( & 37.0_wp * ( v(k,j+1,i) + v(k,j,i) ) & - 8.0_wp * ( v(k,j+2,i) + v(k,j-1,i) ) & + ( v(k,j+3,i) + v(k,j-2,i) ) ) * adv_mom_5 diss_n(k) = - ABS( v_comp(k) - gv ) * ( & 10.0_wp * ( v(k,j+1,i) - v(k,j,i) ) & - 5.0_wp * ( v(k,j+2,i) - v(k,j-1,i) ) & + ( v(k,j+3,i) - v(k,j-2,i) ) ) * adv_mom_5 ! !-- k index has to be modified near bottom and top, else array !-- subscripts will be exceeded. ibit26 = IBITS(wall_flags_0(k,j,i),26,1) ibit25 = IBITS(wall_flags_0(k,j,i),25,1) ibit24 = IBITS(wall_flags_0(k,j,i),24,1) k_ppp = k + 3 * ibit26 k_pp = k + 2 * ( 1 - ibit24 ) k_mm = k - 2 * ibit26 w_comp = w(k,j-1,i) + w(k,j,i) flux_t(k) = w_comp * ( & ( 37.0_wp * ibit26 * adv_mom_5 & + 7.0_wp * ibit25 * adv_mom_3 & + ibit24 * adv_mom_1 & ) * & ( v(k+1,j,i) + v(k,j,i) ) & - ( 8.0_wp * ibit26 * adv_mom_5 & + ibit25 * adv_mom_3 & ) * & ( v(k_pp,j,i) + v(k-1,j,i) ) & + ( ibit26 * adv_mom_5 & ) * & ( v(k_ppp,j,i) + v(k_mm,j,i) ) & ) diss_t(k) = - ABS( w_comp ) * ( & ( 10.0_wp * ibit26 * adv_mom_5 & + 3.0_wp * ibit25 * adv_mom_3 & + ibit24 * adv_mom_1 & ) * & ( v(k+1,j,i) - v(k,j,i) ) & - ( 5.0_wp * ibit26 * adv_mom_5 & + ibit25 * adv_mom_3 & ) * & ( v(k_pp,j,i) - v(k-1,j,i) ) & + ( ibit26 * adv_mom_5 & ) * & ( v(k_ppp,j,i) - v(k_mm,j,i) ) & ) ! !-- Calculate the divergence of the velocity field. A respective !-- correction is needed to overcome numerical instabilities caused !-- by a not sufficient reduction of divergences near topography. div = ( ( u_comp + gu - ( u(k,j-1,i) + u(k,j,i) ) ) * ddx & + ( v_comp(k) - ( v(k,j,i) + v(k,j-1,i) ) ) * ddy & + ( w_comp - ( w(k-1,j-1,i) + w(k-1,j,i) ) ) & * ddzw(k) & ) * 0.5_wp tend(k,j,i) = tend(k,j,i) - ( & ( flux_r(k) + diss_r(k) & - swap_flux_x_local_v(k,j) - swap_diss_x_local_v(k,j) & ) * ddx & + ( flux_n(k) + diss_n(k) & - swap_flux_y_local_v(k) - swap_diss_y_local_v(k) & ) * ddy & + ( flux_t(k) + diss_t(k) & - flux_d - diss_d & ) * ddzw(k) & ) + v(k,j,i) * div swap_flux_x_local_v(k,j) = flux_r(k) swap_diss_x_local_v(k,j) = diss_r(k) swap_flux_y_local_v(k) = flux_n(k) swap_diss_y_local_v(k) = diss_n(k) flux_d = flux_t(k) diss_d = diss_t(k) ! !-- Statistical Evaluation of v'v'. The factor has to be applied !-- for right evaluation when gallilei_trans = .T. . sums_vs2_ws_l(k,tn) = sums_vs2_ws_l(k,tn) & + ( flux_n(k) & * ( v_comp(k) - 2.0_wp * hom(k,1,2,0) ) & / ( v_comp(k) - gv + 1.0E-20_wp ) & + diss_n(k) & * ABS( v_comp(k) - 2.0_wp * hom(k,1,2,0) ) & / ( ABS( v_comp(k) - gv ) +1.0E-20_wp ) ) & * weight_substep(intermediate_timestep_count) ! !-- Statistical Evaluation of w'v'. sums_wsvs_ws_l(k,tn) = sums_wsvs_ws_l(k,tn) & + ( flux_t(k) + diss_t(k) ) & * weight_substep(intermediate_timestep_count) ENDDO ENDDO ENDDO sums_vs2_ws_l(nzb,tn) = sums_vs2_ws_l(nzb+1,tn) END SUBROUTINE advec_v_ws !------------------------------------------------------------------------------! ! Advection of v - Call for all grid points - accelerator version !------------------------------------------------------------------------------! SUBROUTINE advec_v_ws_acc USE arrays_3d, & ONLY: ddzw, tend, u, v, w USE constants, & ONLY: adv_mom_1, adv_mom_3, adv_mom_5 USE control_parameters, & ONLY: intermediate_timestep_count, u_gtrans, v_gtrans USE grid_variables, & ONLY: ddx, ddy USE indices, & ONLY: i_left, i_right, j_north, j_south, nxl, nxr, nyn, nys, nzb, & nzb_max, nzt, wall_flags_0 USE kinds ! USE statistics, & ! ONLY: hom, sums_vs2_ws_l, sums_wsvs_ws_l, weight_substep IMPLICIT NONE INTEGER(iwp) :: i !: INTEGER(iwp) :: ibit18 !: INTEGER(iwp) :: ibit19 !: INTEGER(iwp) :: ibit20 !: INTEGER(iwp) :: ibit21 !: INTEGER(iwp) :: ibit22 !: INTEGER(iwp) :: ibit23 !: INTEGER(iwp) :: ibit24 !: INTEGER(iwp) :: ibit25 !: INTEGER(iwp) :: ibit26 !: INTEGER(iwp) :: j !: INTEGER(iwp) :: k !: INTEGER(iwp) :: k_mm !: INTEGER(iwp) :: k_mmm !: INTEGER(iwp) :: k_pp !: INTEGER(iwp) :: k_ppp !: INTEGER(iwp) :: tn = 0 !: REAL(wp) :: diss_d !: REAL(wp) :: diss_l !: REAL(wp) :: diss_n !: REAL(wp) :: diss_r !: REAL(wp) :: diss_s !: REAL(wp) :: diss_t !: REAL(wp) :: div !: REAL(wp) :: flux_d !: REAL(wp) :: flux_l !: REAL(wp) :: flux_n !: REAL(wp) :: flux_r !: REAL(wp) :: flux_s !: REAL(wp) :: flux_t !: REAL(wp) :: gu !: REAL(wp) :: gv !: REAL(wp) :: u_comp !: REAL(wp) :: u_comp_l !: REAL(wp) :: v_comp !: REAL(wp) :: v_comp_s !: REAL(wp) :: w_comp !: gu = 2.0_wp * u_gtrans gv = 2.0_wp * v_gtrans ! !-- Computation of fluxes and tendency terms !$acc kernels present( ddzw, tend, u, v, w, wall_flags_0 ) DO i = i_left, i_right DO j = j_south, j_north DO k = nzb+1, nzt ibit20 = IBITS(wall_flags_0(k,j,i),20,1) ibit19 = IBITS(wall_flags_0(k,j,i),19,1) ibit18 = IBITS(wall_flags_0(k,j,i),18,1) u_comp_l = u(k,j-1,i) + u(k,j,i) - gu flux_l = u_comp_l * ( & ( 37.0_wp * ibit20 * adv_mom_5 & + 7.0_wp * ibit19 * adv_mom_3 & + ibit18 * adv_mom_1 & ) * & ( v(k,j,i) + v(k,j,i-1) ) & - ( 8.0_wp * ibit20 * adv_mom_5 & + ibit19 * adv_mom_3 & ) * & ( v(k,j,i+1) + v(k,j,i-2) ) & + ( ibit20 * adv_mom_5 & ) * & ( v(k,j,i+2) + v(k,j,i-3) ) & ) diss_l = - ABS( u_comp_l ) * ( & ( 10.0_wp * ibit20 * adv_mom_5 & + 3.0_wp * ibit19 * adv_mom_3 & + ibit18 * adv_mom_1 & ) * & ( v(k,j,i) - v(k,j,i-1) ) & - ( 5.0_wp * ibit20 * adv_mom_5 & + ibit19 * adv_mom_3 & ) * & ( v(k,j,i+1) - v(k,j,i-2) ) & + ( ibit20 * adv_mom_5 & ) * & ( v(k,j,i+2) - v(k,j,i-3) ) & ) u_comp = u(k,j-1,i+1) + u(k,j,i+1) - gu flux_r = u_comp * ( & ( 37.0_wp * ibit20 * adv_mom_5 & + 7.0_wp * ibit19 * adv_mom_3 & + ibit18 * adv_mom_1 & ) * & ( v(k,j,i+1) + v(k,j,i) ) & - ( 8.0_wp * ibit20 * adv_mom_5 & + ibit19 * adv_mom_3 & ) * & ( v(k,j,i+2) + v(k,j,i-1) ) & + ( ibit20 * adv_mom_5 & ) * & ( v(k,j,i+3) + v(k,j,i-2) ) & ) diss_r = - ABS( u_comp ) * ( & ( 10.0_wp * ibit20 * adv_mom_5 & + 3.0_wp * ibit19 * adv_mom_3 & + ibit18 * adv_mom_1 & ) * & ( v(k,j,i+1) - v(k,j,i) ) & - ( 5.0_wp * ibit20 * adv_mom_5 & + ibit19 * adv_mom_3 & ) * & ( v(k,j,i+2) - v(k,j,i-1) ) & + ( ibit20 * adv_mom_5 & ) * & ( v(k,j,i+3) - v(k,j,i-2) ) & ) ibit23 = IBITS(wall_flags_0(k,j,i),23,1) ibit22 = IBITS(wall_flags_0(k,j,i),22,1) ibit21 = IBITS(wall_flags_0(k,j,i),21,1) v_comp_s = v(k,j,i) + v(k,j-1,i) - gv flux_s = v_comp_s * ( & ( 37.0_wp * ibit23 * adv_mom_5 & + 7.0_wp * ibit22 * adv_mom_3 & + ibit21 * adv_mom_1 & ) * & ( v(k,j,i) + v(k,j-1,i) ) & - ( 8.0_wp * ibit23 * adv_mom_5 & + ibit22 * adv_mom_3 & ) * & ( v(k,j+1,i) + v(k,j-2,i) ) & + ( ibit23 * adv_mom_5 & ) * & ( v(k,j+2,i) + v(k,j-3,i) ) & ) diss_s = - ABS( v_comp_s ) * ( & ( 10.0_wp * ibit23 * adv_mom_5 & + 3.0_wp * ibit22 * adv_mom_3 & + ibit21 * adv_mom_1 & ) * & ( v(k,j,i) - v(k,j-1,i) ) & - ( 5.0_wp * ibit23 * adv_mom_5 & + ibit22 * adv_mom_3 & ) * & ( v(k,j+1,i) - v(k,j-2,i) ) & + ( ibit23 * adv_mom_5 & ) * & ( v(k,j+2,i) - v(k,j-3,i) ) & ) v_comp = v(k,j+1,i) + v(k,j,i) flux_n = ( v_comp - gv ) * ( & ( 37.0_wp * ibit23 * adv_mom_5 & + 7.0_wp * ibit22 * adv_mom_3 & + ibit21 * adv_mom_1 & ) * & ( v(k,j+1,i) + v(k,j,i) ) & - ( 8.0_wp * ibit23 * adv_mom_5 & + ibit22 * adv_mom_3 & ) * & ( v(k,j+2,i) + v(k,j-1,i) ) & + ( ibit23 * adv_mom_5 & ) * & ( v(k,j+3,i) + v(k,j-2,i) ) & ) diss_n = - ABS( v_comp - gv ) * ( & ( 10.0_wp * ibit23 * adv_mom_5 & + 3.0_wp * ibit22 * adv_mom_3 & + ibit21 * adv_mom_1 & ) * & ( v(k,j+1,i) - v(k,j,i) ) & - ( 5.0_wp * ibit23 * adv_mom_5 & + ibit22 * adv_mom_3 & ) * & ( v(k,j+2,i) - v(k,j-1,i) ) & + ( ibit23 * adv_mom_5 & ) * & ( v(k,j+3,i) - v(k,j-2,i) ) & ) ibit26 = IBITS(wall_flags_0(k-1,j,i),26,1) ibit25 = IBITS(wall_flags_0(k-1,j,i),25,1) ibit24 = IBITS(wall_flags_0(k-1,j,i),24,1) k_pp = k + 2 * ibit26 k_mm = k - 2 * ( ibit25 + ibit26 ) k_mmm = k - 3 * ibit26 w_comp = w(k-1,j-1,i) + w(k-1,j,i) flux_d = w_comp * ( & ( 37.0_wp * ibit26 * adv_mom_5 & + 7.0_wp * ibit25 * adv_mom_3 & + ibit24 * adv_mom_1 & ) * & ( v(k,j,i) + v(k-1,j,i) ) & - ( 8.0_wp * ibit26 * adv_mom_5 & + ibit25 * adv_mom_3 & ) * & ( v(k+1,j,i) + v(k_mm,j,i) ) & + ( ibit26 * adv_mom_5 & ) * & ( v(k_pp,j,i) + v(k_mmm,j,i) ) & ) diss_d = - ABS( w_comp ) * ( & ( 10.0_wp * ibit26 * adv_mom_5 & + 3.0_wp * ibit25 * adv_mom_3 & + ibit24 * adv_mom_1 & ) * & ( v(k,j,i) - v(k-1,j,i) ) & - ( 5.0_wp * ibit26 * adv_mom_5 & + ibit25 * adv_mom_3 & ) * & ( v(k+1,j,i) - v(k_mm,j,i) ) & + ( ibit26 * adv_mom_5 & ) * & ( v(k_pp,j,i) - v(k_mmm,j,i) ) & ) ! !-- k index has to be modified near bottom and top, else array !-- subscripts will be exceeded. ibit26 = IBITS(wall_flags_0(k,j,i),26,1) ibit25 = IBITS(wall_flags_0(k,j,i),25,1) ibit24 = IBITS(wall_flags_0(k,j,i),24,1) k_ppp = k + 3 * ibit26 k_pp = k + 2 * ( 1 - ibit24 ) k_mm = k - 2 * ibit26 w_comp = w(k,j-1,i) + w(k,j,i) flux_t = w_comp * ( & ( 37.0_wp * ibit26 * adv_mom_5 & + 7.0_wp * ibit25 * adv_mom_3 & + ibit24 * adv_mom_1 & ) * & ( v(k+1,j,i) + v(k,j,i) ) & - ( 8.0_wp * ibit26 * adv_mom_5 & + ibit25 * adv_mom_3 & ) * & ( v(k_pp,j,i) + v(k-1,j,i) ) & + ( ibit26 * adv_mom_5 & ) * & ( v(k_ppp,j,i) + v(k_mm,j,i) ) & ) diss_t = - ABS( w_comp ) * ( & ( 10.0_wp * ibit26 * adv_mom_5 & + 3.0_wp * ibit25 * adv_mom_3 & + ibit24 * adv_mom_1 & ) * & ( v(k+1,j,i) - v(k,j,i) ) & - ( 5.0_wp * ibit26 * adv_mom_5 & + ibit25 * adv_mom_3 & ) * & ( v(k_pp,j,i) - v(k-1,j,i) ) & + ( ibit26 * adv_mom_5 & ) * & ( v(k_ppp,j,i) - v(k_mm,j,i) ) & ) ! !-- Calculate the divergence of the velocity field. A respective !-- correction is needed to overcome numerical instabilities caused !-- by a not sufficient reduction of divergences near topography. div = ( ( u_comp + gu - ( u(k,j-1,i) + u(k,j,i) ) ) * ddx & + ( v_comp - ( v(k,j,i) + v(k,j-1,i) ) ) * ddy & + ( w_comp - ( w(k-1,j-1,i) + w(k-1,j,i) ) & ) * ddzw(k) & ) * 0.5_wp tend(k,j,i) = - ( & ( flux_r + diss_r - flux_l - diss_l ) * ddx & + ( flux_n + diss_n - flux_s - diss_s ) * ddy & + ( flux_t + diss_t - flux_d - diss_d ) * ddzw(k) & ) + div * v(k,j,i) !++ !-- Statistical Evaluation of v'v'. The factor has to be applied !-- for right evaluation when gallilei_trans = .T. . ! sums_vs2_ws_l(k,tn) = sums_vs2_ws_l(k,tn) & ! + ( flux_n & ! * ( v_comp - 2.0_wp * hom(k,1,2,0) ) & ! / ( v_comp - gv + 1.0E-20_wp ) & ! + diss_n & ! * ABS( v_comp - 2.0_wp * hom(k,1,2,0) ) & ! / ( ABS( v_comp - gv ) +1.0E-20_wp ) ) & ! * weight_substep(intermediate_timestep_count) ! !-- Statistical Evaluation of w'v'. ! sums_wsvs_ws_l(k,tn) = sums_wsvs_ws_l(k,tn) & ! + ( flux_t + diss_t ) & ! * weight_substep(intermediate_timestep_count) ENDDO ENDDO ENDDO !$acc end kernels !++ ! sums_vs2_ws_l(nzb,tn) = sums_vs2_ws_l(nzb+1,tn) END SUBROUTINE advec_v_ws_acc !------------------------------------------------------------------------------! ! Advection of w - Call for all grid points !------------------------------------------------------------------------------! SUBROUTINE advec_w_ws USE arrays_3d, & ONLY: ddzu, tend, u, v, w USE constants, & ONLY: adv_mom_1, adv_mom_3, adv_mom_5 USE control_parameters, & ONLY: intermediate_timestep_count, u_gtrans, v_gtrans USE grid_variables, & ONLY: ddx, ddy USE indices, & ONLY: nxl, nxr, nyn, nys, nzb, nzb_max, nzt, wall_flags_0, & wall_flags_00 USE kinds USE statistics, & ONLY: hom, sums_ws2_ws_l, weight_substep IMPLICIT NONE INTEGER(iwp) :: i !: INTEGER(iwp) :: ibit27 !: INTEGER(iwp) :: ibit28 !: INTEGER(iwp) :: ibit29 !: INTEGER(iwp) :: ibit30 !: INTEGER(iwp) :: ibit31 !: INTEGER(iwp) :: ibit32 !: INTEGER(iwp) :: ibit33 !: INTEGER(iwp) :: ibit34 !: INTEGER(iwp) :: ibit35 !: INTEGER(iwp) :: j !: INTEGER(iwp) :: k !: INTEGER(iwp) :: k_mm !: INTEGER(iwp) :: k_pp !: INTEGER(iwp) :: k_ppp !: INTEGER(iwp) :: tn = 0 !: REAL(wp) :: diss_d !: REAL(wp) :: div !: REAL(wp) :: flux_d !: REAL(wp) :: gu !: REAL(wp) :: gv !: REAL(wp) :: u_comp !: REAL(wp) :: v_comp !: REAL(wp) :: w_comp !: REAL(wp), DIMENSION(nzb:nzt) :: diss_t !: REAL(wp), DIMENSION(nzb:nzt) :: flux_t !: REAL(wp), DIMENSION(nzb+1:nzt) :: diss_n !: REAL(wp), DIMENSION(nzb+1:nzt) :: diss_r !: REAL(wp), DIMENSION(nzb+1:nzt) :: flux_n !: REAL(wp), DIMENSION(nzb+1:nzt) :: flux_r !: REAL(wp), DIMENSION(nzb+1:nzt) :: swap_diss_y_local_w !: REAL(wp), DIMENSION(nzb+1:nzt) :: swap_flux_y_local_w !: REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) :: swap_diss_x_local_w !: REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) :: swap_flux_x_local_w !: gu = 2.0_wp * u_gtrans gv = 2.0_wp * v_gtrans ! !-- compute the whole left boundary of the processor domain i = nxl DO j = nys, nyn DO k = nzb+1, nzb_max ibit29 = IBITS(wall_flags_0(k,j,i),29,1) ibit28 = IBITS(wall_flags_0(k,j,i),28,1) ibit27 = IBITS(wall_flags_0(k,j,i),27,1) u_comp = u(k+1,j,i) + u(k,j,i) - gu swap_flux_x_local_w(k,j) = u_comp * ( & ( 37.0_wp * ibit29 * adv_mom_5 & + 7.0_wp * ibit28 * adv_mom_3 & + ibit27 * adv_mom_1 & ) * & ( w(k,j,i) + w(k,j,i-1) ) & - ( 8.0_wp * ibit29 * adv_mom_5 & + ibit28 * adv_mom_3 & ) * & ( w(k,j,i+1) + w(k,j,i-2) ) & + ( ibit29 * adv_mom_5 & ) * & ( w(k,j,i+2) + w(k,j,i-3) ) & ) swap_diss_x_local_w(k,j) = - ABS( u_comp ) * ( & ( 10.0_wp * ibit29 * adv_mom_5 & + 3.0_wp * ibit28 * adv_mom_3 & + ibit27 * adv_mom_1 & ) * & ( w(k,j,i) - w(k,j,i-1) ) & - ( 5.0_wp * ibit29 * adv_mom_5 & + ibit28 * adv_mom_3 & ) * & ( w(k,j,i+1) - w(k,j,i-2) ) & + ( ibit29 * adv_mom_5 & ) * & ( w(k,j,i+2) - w(k,j,i-3) ) & ) ENDDO DO k = nzb_max+1, nzt u_comp = u(k+1,j,i) + u(k,j,i) - gu swap_flux_x_local_w(k,j) = u_comp * ( & 37.0_wp * ( w(k,j,i) + w(k,j,i-1) ) & - 8.0_wp * ( w(k,j,i+1) + w(k,j,i-2) ) & + ( w(k,j,i+2) + w(k,j,i-3) ) ) * adv_mom_5 swap_diss_x_local_w(k,j) = - ABS( u_comp ) * ( & 10.0_wp * ( w(k,j,i) - w(k,j,i-1) ) & - 5.0_wp * ( w(k,j,i+1) - w(k,j,i-2) ) & + ( w(k,j,i+2) - w(k,j,i-3) ) ) * adv_mom_5 ENDDO ENDDO DO i = nxl, nxr j = nys DO k = nzb+1, nzb_max ibit32 = IBITS(wall_flags_00(k,j,i),0,1) ibit31 = IBITS(wall_flags_0(k,j,i),31,1) ibit30 = IBITS(wall_flags_0(k,j,i),30,1) v_comp = v(k+1,j,i) + v(k,j,i) - gv swap_flux_y_local_w(k) = v_comp * ( & ( 37.0_wp * ibit32 * adv_mom_5 & + 7.0_wp * ibit31 * adv_mom_3 & + ibit30 * adv_mom_1 & ) * & ( w(k,j,i) + w(k,j-1,i) ) & - ( 8.0_wp * ibit32 * adv_mom_5 & + ibit31 * adv_mom_3 & ) * & ( w(k,j+1,i) + w(k,j-2,i) ) & + ( ibit32 * adv_mom_5 & ) * & ( w(k,j+2,i) + w(k,j-3,i) ) & ) swap_diss_y_local_w(k) = - ABS( v_comp ) * ( & ( 10.0_wp * ibit32 * adv_mom_5 & + 3.0_wp * ibit31 * adv_mom_3 & + ibit30 * adv_mom_1 & ) * & ( w(k,j,i) - w(k,j-1,i) ) & - ( 5.0_wp * ibit32 * adv_mom_5 & + ibit31 * adv_mom_3 & ) * & ( w(k,j+1,i) - w(k,j-2,i) ) & + ( ibit32 * adv_mom_5 & ) * & ( w(k,j+2,i) - w(k,j-3,i) ) & ) ENDDO DO k = nzb_max+1, nzt v_comp = v(k+1,j,i) + v(k,j,i) - gv swap_flux_y_local_w(k) = v_comp * ( & 37.0_wp * ( w(k,j,i) + w(k,j-1,i) ) & - 8.0_wp * ( w(k,j+1,i) +w(k,j-2,i) ) & + ( w(k,j+2,i) + w(k,j-3,i) ) ) * adv_mom_5 swap_diss_y_local_w(k) = - ABS( v_comp ) * ( & 10.0_wp * ( w(k,j,i) - w(k,j-1,i) ) & - 5.0_wp * ( w(k,j+1,i) - w(k,j-2,i) ) & + ( w(k,j+2,i) - w(k,j-3,i) ) ) * adv_mom_5 ENDDO DO j = nys, nyn ! !-- The lower flux has to be calculated explicetely for the tendency !-- at the first w-level. For topography wall this is done implicitely !-- by wall_flags_0. k = nzb + 1 w_comp = w(k,j,i) + w(k-1,j,i) flux_t(0) = w_comp * ( w(k,j,i) + w(k-1,j,i) ) * adv_mom_1 diss_t(0) = -ABS(w_comp) * ( w(k,j,i) - w(k-1,j,i) ) * adv_mom_1 flux_d = flux_t(0) diss_d = diss_t(0) DO k = nzb+1, nzb_max ibit29 = IBITS(wall_flags_0(k,j,i),29,1) ibit28 = IBITS(wall_flags_0(k,j,i),28,1) ibit27 = IBITS(wall_flags_0(k,j,i),27,1) u_comp = u(k+1,j,i+1) + u(k,j,i+1) - gu flux_r(k) = u_comp * ( & ( 37.0_wp * ibit29 * adv_mom_5 & + 7.0_wp * ibit28 * adv_mom_3 & + ibit27 * adv_mom_1 & ) * & ( w(k,j,i+1) + w(k,j,i) ) & - ( 8.0_wp * ibit29 * adv_mom_5 & + ibit28 * adv_mom_3 & ) * & ( w(k,j,i+2) + w(k,j,i-1) ) & + ( ibit29 * adv_mom_5 & ) * & ( w(k,j,i+3) + w(k,j,i-2) ) & ) diss_r(k) = - ABS( u_comp ) * ( & ( 10.0_wp * ibit29 * adv_mom_5 & + 3.0_wp * ibit28 * adv_mom_3 & + ibit27 * adv_mom_1 & ) * & ( w(k,j,i+1) - w(k,j,i) ) & - ( 5.0_wp * ibit29 * adv_mom_5 & + ibit28 * adv_mom_3 & ) * & ( w(k,j,i+2) - w(k,j,i-1) ) & + ( ibit29 * adv_mom_5 & ) * & ( w(k,j,i+3) - w(k,j,i-2) ) & ) ibit32 = IBITS(wall_flags_00(k,j,i),0,1) ibit31 = IBITS(wall_flags_0(k,j,i),31,1) ibit30 = IBITS(wall_flags_0(k,j,i),30,1) v_comp = v(k+1,j+1,i) + v(k,j+1,i) - gv flux_n(k) = v_comp * ( & ( 37.0_wp * ibit32 * adv_mom_5 & + 7.0_wp * ibit31 * adv_mom_3 & + ibit30 * adv_mom_1 & ) * & ( w(k,j+1,i) + w(k,j,i) ) & - ( 8.0_wp * ibit32 * adv_mom_5 & + ibit31 * adv_mom_3 & ) * & ( w(k,j+2,i) + w(k,j-1,i) ) & + ( ibit32 * adv_mom_5 & ) * & ( w(k,j+3,i) + w(k,j-2,i) ) & ) diss_n(k) = - ABS( v_comp ) * ( & ( 10.0_wp * ibit32 * adv_mom_5 & + 3.0_wp * ibit31 * adv_mom_3 & + ibit30 * adv_mom_1 & ) * & ( w(k,j+1,i) - w(k,j,i) ) & - ( 5.0_wp * ibit32 * adv_mom_5 & + ibit31 * adv_mom_3 & ) * & ( w(k,j+2,i) - w(k,j-1,i) ) & + ( ibit32 * adv_mom_5 & ) * & ( w(k,j+3,i) - w(k,j-2,i) ) & ) ! !-- k index has to be modified near bottom and top, else array !-- subscripts will be exceeded. ibit35 = IBITS(wall_flags_00(k,j,i),3,1) ibit34 = IBITS(wall_flags_00(k,j,i),2,1) ibit33 = IBITS(wall_flags_00(k,j,i),1,1) k_ppp = k + 3 * ibit35 k_pp = k + 2 * ( 1 - ibit33 ) k_mm = k - 2 * ibit35 w_comp = w(k+1,j,i) + w(k,j,i) flux_t(k) = w_comp * ( & ( 37.0_wp * ibit35 * adv_mom_5 & + 7.0_wp * ibit34 * adv_mom_3 & + ibit33 * adv_mom_1 & ) * & ( w(k+1,j,i) + w(k,j,i) ) & - ( 8.0_wp * ibit35 * adv_mom_5 & + ibit34 * adv_mom_3 & ) * & ( w(k_pp,j,i) + w(k-1,j,i) ) & + ( ibit35 * adv_mom_5 & ) * & ( w(k_ppp,j,i) + w(k_mm,j,i) ) & ) diss_t(k) = - ABS( w_comp ) * ( & ( 10.0_wp * ibit35 * adv_mom_5 & + 3.0_wp * ibit34 * adv_mom_3 & + ibit33 * adv_mom_1 & ) * & ( w(k+1,j,i) - w(k,j,i) ) & - ( 5.0_wp * ibit35 * adv_mom_5 & + ibit34 * adv_mom_3 & ) * & ( w(k_pp,j,i) - w(k-1,j,i) ) & + ( ibit35 * adv_mom_5 & ) * & ( w(k_ppp,j,i) - w(k_mm,j,i) ) & ) ! !-- Calculate the divergence of the velocity field. A respective !-- correction is needed to overcome numerical instabilities caused !-- by a not sufficient reduction of divergences near topography. div = ( ( u_comp + gu - ( u(k+1,j,i) + u(k,j,i) ) ) * ddx & + ( v_comp + gv - ( v(k+1,j,i) + v(k,j,i) ) ) * ddy & + ( w_comp - ( w(k,j,i) + w(k-1,j,i) ) ) & * ddzu(k+1) & ) * 0.5_wp tend(k,j,i) = tend(k,j,i) - ( & ( flux_r(k) + diss_r(k) & - swap_flux_x_local_w(k,j) - swap_diss_x_local_w(k,j) & ) * ddx & + ( flux_n(k) + diss_n(k) & - swap_flux_y_local_w(k) - swap_diss_y_local_w(k) & ) * ddy & + ( flux_t(k) + diss_t(k) & - flux_d - diss_d & ) * ddzu(k+1) & ) + div * w(k,j,i) swap_flux_x_local_w(k,j) = flux_r(k) swap_diss_x_local_w(k,j) = diss_r(k) swap_flux_y_local_w(k) = flux_n(k) swap_diss_y_local_w(k) = diss_n(k) flux_d = flux_t(k) diss_d = diss_t(k) sums_ws2_ws_l(k,tn) = sums_ws2_ws_l(k,tn) & + ( flux_t(k) + diss_t(k) ) & * weight_substep(intermediate_timestep_count) ENDDO DO k = nzb_max+1, nzt u_comp = u(k+1,j,i+1) + u(k,j,i+1) - gu flux_r(k) = u_comp * ( & 37.0_wp * ( w(k,j,i+1) + w(k,j,i) ) & - 8.0_wp * ( w(k,j,i+2) + w(k,j,i-1) ) & + ( w(k,j,i+3) + w(k,j,i-2) ) ) * adv_mom_5 diss_r(k) = - ABS( u_comp ) * ( & 10.0_wp * ( w(k,j,i+1) - w(k,j,i) ) & - 5.0_wp * ( w(k,j,i+2) - w(k,j,i-1) ) & + ( w(k,j,i+3) - w(k,j,i-2) ) ) * adv_mom_5 v_comp = v(k+1,j+1,i) + v(k,j+1,i) - gv flux_n(k) = v_comp * ( & 37.0_wp * ( w(k,j+1,i) + w(k,j,i) ) & - 8.0_wp * ( w(k,j+2,i) + w(k,j-1,i) ) & + ( w(k,j+3,i) + w(k,j-2,i) ) ) * adv_mom_5 diss_n(k) = - ABS( v_comp ) * ( & 10.0_wp * ( w(k,j+1,i) - w(k,j,i) ) & - 5.0_wp * ( w(k,j+2,i) - w(k,j-1,i) ) & + ( w(k,j+3,i) - w(k,j-2,i) ) ) * adv_mom_5 ! !-- k index has to be modified near bottom and top, else array !-- subscripts will be exceeded. ibit35 = IBITS(wall_flags_00(k,j,i),3,1) ibit34 = IBITS(wall_flags_00(k,j,i),2,1) ibit33 = IBITS(wall_flags_00(k,j,i),1,1) k_ppp = k + 3 * ibit35 k_pp = k + 2 * ( 1 - ibit33 ) k_mm = k - 2 * ibit35 w_comp = w(k+1,j,i) + w(k,j,i) flux_t(k) = w_comp * ( & ( 37.0_wp * ibit35 * adv_mom_5 & + 7.0_wp * ibit34 * adv_mom_3 & + ibit33 * adv_mom_1 & ) * & ( w(k+1,j,i) + w(k,j,i) ) & - ( 8.0_wp * ibit35 * adv_mom_5 & + ibit34 * adv_mom_3 & ) * & ( w(k_pp,j,i) + w(k-1,j,i) ) & + ( ibit35 * adv_mom_5 & ) * & ( w(k_ppp,j,i) + w(k_mm,j,i) ) & ) diss_t(k) = - ABS( w_comp ) * ( & ( 10.0_wp * ibit35 * adv_mom_5 & + 3.0_wp * ibit34 * adv_mom_3 & + ibit33 * adv_mom_1 & ) * & ( w(k+1,j,i) - w(k,j,i) ) & - ( 5.0_wp * ibit35 * adv_mom_5 & + ibit34 * adv_mom_3 & ) * & ( w(k_pp,j,i) - w(k-1,j,i) ) & + ( ibit35 * adv_mom_5 & ) * & ( w(k_ppp,j,i) - w(k_mm,j,i) ) & ) ! !-- Calculate the divergence of the velocity field. A respective !-- correction is needed to overcome numerical instabilities caused !-- by a not sufficient reduction of divergences near topography. div = ( ( u_comp + gu - ( u(k+1,j,i) + u(k,j,i) ) ) * ddx & + ( v_comp + gv - ( v(k+1,j,i) + v(k,j,i) ) ) * ddy & + ( w_comp - ( w(k,j,i) + w(k-1,j,i) ) ) & * ddzu(k+1) & ) * 0.5_wp tend(k,j,i) = tend(k,j,i) - ( & ( flux_r(k) + diss_r(k) & - swap_flux_x_local_w(k,j) - swap_diss_x_local_w(k,j) & ) * ddx & + ( flux_n(k) + diss_n(k) & - swap_flux_y_local_w(k) - swap_diss_y_local_w(k) & ) * ddy & + ( flux_t(k) + diss_t(k) & - flux_d - diss_d & ) * ddzu(k+1) & ) + div * w(k,j,i) swap_flux_x_local_w(k,j) = flux_r(k) swap_diss_x_local_w(k,j) = diss_r(k) swap_flux_y_local_w(k) = flux_n(k) swap_diss_y_local_w(k) = diss_n(k) flux_d = flux_t(k) diss_d = diss_t(k) sums_ws2_ws_l(k,tn) = sums_ws2_ws_l(k,tn) & + ( flux_t(k) + diss_t(k) ) & * weight_substep(intermediate_timestep_count) ENDDO ENDDO ENDDO END SUBROUTINE advec_w_ws !------------------------------------------------------------------------------! ! Advection of w - Call for all grid points - accelerator version !------------------------------------------------------------------------------! SUBROUTINE advec_w_ws_acc USE arrays_3d, & ONLY: ddzu, tend, u, v, w USE constants, & ONLY: adv_mom_1, adv_mom_3, adv_mom_5 USE control_parameters, & ONLY: intermediate_timestep_count, u_gtrans, v_gtrans USE grid_variables, & ONLY: ddx, ddy USE indices, & ONLY: i_left, i_right, j_north, j_south, nxl, nxr, nyn, nys, nzb, & nzb_max, nzt, wall_flags_0, wall_flags_00 USE kinds ! USE statistics, & ! ONLY: hom, sums_ws2_ws_l, weight_substep IMPLICIT NONE INTEGER(iwp) :: i !: INTEGER(iwp) :: ibit27 !: INTEGER(iwp) :: ibit28 !: INTEGER(iwp) :: ibit29 !: INTEGER(iwp) :: ibit30 !: INTEGER(iwp) :: ibit31 !: INTEGER(iwp) :: ibit32 !: INTEGER(iwp) :: ibit33 !: INTEGER(iwp) :: ibit34 !: INTEGER(iwp) :: ibit35 !: INTEGER(iwp) :: j !: INTEGER(iwp) :: k !: INTEGER(iwp) :: k_mmm !: INTEGER(iwp) :: k_mm !: INTEGER(iwp) :: k_pp !: INTEGER(iwp) :: k_ppp !: INTEGER(iwp) :: tn = 0 !: REAL(wp) :: diss_d !: REAL(wp) :: diss_l !: REAL(wp) :: diss_n !: REAL(wp) :: diss_r !: REAL(wp) :: diss_s !: REAL(wp) :: diss_t !: REAL(wp) :: div !: REAL(wp) :: flux_d !: REAL(wp) :: flux_l !: REAL(wp) :: flux_n !: REAL(wp) :: flux_r !: REAL(wp) :: flux_s !: REAL(wp) :: flux_t !: REAL(wp) :: gu !: REAL(wp) :: gv !: REAL(wp) :: u_comp !: REAL(wp) :: u_comp_l !: REAL(wp) :: v_comp !: REAL(wp) :: v_comp_s !: REAL(wp) :: w_comp !: gu = 2.0_wp * u_gtrans gv = 2.0_wp * v_gtrans ! !-- Computation of fluxes and tendency terms !$acc kernels present( ddzu, tend, u, v, w, wall_flags_0, wall_flags_00 ) DO i = i_left, i_right DO j = j_south, j_north DO k = nzb+1, nzt ibit27 = IBITS(wall_flags_0(k,j,i),27,1) ibit28 = IBITS(wall_flags_0(k,j,i),28,1) ibit29 = IBITS(wall_flags_0(k,j,i),29,1) u_comp_l = u(k+1,j,i) + u(k,j,i) - gu flux_l = u_comp_l * ( & ( 37.0_wp * ibit29 * adv_mom_5 & + 7.0_wp * ibit28 * adv_mom_3 & + ibit27 * adv_mom_1 & ) * & ( w(k,j,i) + w(k,j,i-1) ) & - ( 8.0_wp * ibit29 * adv_mom_5 & + ibit28 * adv_mom_3 & ) * & ( w(k,j,i+1) + w(k,j,i-2) ) & + ( ibit29 * adv_mom_5 & ) * & ( w(k,j,i+2) + w(k,j,i-3) ) & ) diss_l = - ABS( u_comp_l ) * ( & ( 10.0_wp * ibit29 * adv_mom_5 & + 3.0_wp * ibit28 * adv_mom_3 & + ibit27 * adv_mom_1 & ) * & ( w(k,j,i) - w(k,j,i-1) ) & - ( 5.0_wp * ibit29 * adv_mom_5 & + ibit28 * adv_mom_3 & ) * & ( w(k,j,i+1) - w(k,j,i-2) ) & + ( ibit29 * adv_mom_5 & ) * & ( w(k,j,i+2) - w(k,j,i-3) ) & ) u_comp = u(k+1,j,i+1) + u(k,j,i+1) - gu flux_r = u_comp * ( & ( 37.0_wp * ibit29 * adv_mom_5 & + 7.0_wp * ibit28 * adv_mom_3 & + ibit27 * adv_mom_1 & ) * & ( w(k,j,i+1) + w(k,j,i) ) & - ( 8.0_wp * ibit29 * adv_mom_5 & + ibit28 * adv_mom_3 & ) * & ( w(k,j,i+2) + w(k,j,i-1) ) & + ( ibit29 * adv_mom_5 & ) * & ( w(k,j,i+3) + w(k,j,i-2) ) & ) diss_r = - ABS( u_comp ) * ( & ( 10.0_wp * ibit29 * adv_mom_5 & + 3.0_wp * ibit28 * adv_mom_3 & + ibit27 * adv_mom_1 & ) * & ( w(k,j,i+1) - w(k,j,i) ) & - ( 5.0_wp * ibit29 * adv_mom_5 & + ibit28 * adv_mom_3 & ) * & ( w(k,j,i+2) - w(k,j,i-1) ) & + ( ibit29 * adv_mom_5 & ) * & ( w(k,j,i+3) - w(k,j,i-2) ) & ) ibit32 = IBITS(wall_flags_00(k,j,i),0,1) ibit31 = IBITS(wall_flags_0(k,j,i),31,1) ibit30 = IBITS(wall_flags_0(k,j,i),30,1) v_comp_s = v(k+1,j,i) + v(k,j,i) - gv flux_s = v_comp_s * ( & ( 37.0_wp * ibit32 * adv_mom_5 & + 7.0_wp * ibit31 * adv_mom_3 & + ibit30 * adv_mom_1 & ) * & ( w(k,j,i) + w(k,j-1,i) ) & - ( 8.0_wp * ibit32 * adv_mom_5 & + ibit31 * adv_mom_3 & ) * & ( w(k,j+1,i) + w(k,j-2,i) ) & + ( ibit32 * adv_mom_5 & ) * & ( w(k,j+2,i) + w(k,j-3,i) ) & ) diss_s = - ABS( v_comp_s ) * ( & ( 10.0_wp * ibit32 * adv_mom_5 & + 3.0_wp * ibit31 * adv_mom_3 & + ibit30 * adv_mom_1 & ) * & ( w(k,j,i) - w(k,j-1,i) ) & - ( 5.0_wp * ibit32 * adv_mom_5 & + ibit31 * adv_mom_3 & ) * & ( w(k,j+1,i) - w(k,j-2,i) ) & + ( ibit32 * adv_mom_5 & ) * & ( w(k,j+2,i) - w(k,j-3,i) ) & ) v_comp = v(k+1,j+1,i) + v(k,j+1,i) - gv flux_n = v_comp * ( & ( 37.0_wp * ibit32 * adv_mom_5 & + 7.0_wp * ibit31 * adv_mom_3 & + ibit30 * adv_mom_1 & ) * & ( w(k,j+1,i) + w(k,j,i) ) & - ( 8.0_wp * ibit32 * adv_mom_5 & + ibit31 * adv_mom_3 & ) * & ( w(k,j+2,i) + w(k,j-1,i) ) & + ( ibit32 * adv_mom_5 & ) * & ( w(k,j+3,i) + w(k,j-2,i) ) & ) diss_n = - ABS( v_comp ) * ( & ( 10.0_wp * ibit32 * adv_mom_5 & + 3.0_wp * ibit31 * adv_mom_3 & + ibit30 * adv_mom_1 & ) * & ( w(k,j+1,i) - w(k,j,i) ) & - ( 5.0_wp * ibit32 * adv_mom_5 & + ibit31 * adv_mom_3 & ) * & ( w(k,j+2,i) - w(k,j-1,i) ) & + ( ibit32 * adv_mom_5 & ) * & ( w(k,j+3,i) - w(k,j-2,i) ) & ) ibit35 = IBITS(wall_flags_00(k-1,j,i),3,1) ibit34 = IBITS(wall_flags_00(k-1,j,i),2,1) ibit33 = IBITS(wall_flags_00(k-1,j,i),1,1) k_pp = k + 2 * ibit35 k_mm = k - 2 * ( ibit34 + ibit35 ) k_mmm = k - 3 * ibit35 w_comp = w(k,j,i) + w(k-1,j,i) flux_d = w_comp * ( & ( 37.0_wp * ibit35 * adv_mom_5 & + 7.0_wp * ibit34 * adv_mom_3 & + ibit33 * adv_mom_1 & ) * & ( w(k,j,i) + w(k-1,j,i) ) & - ( 8.0_wp * ibit35 * adv_mom_5 & + ibit34 * adv_mom_3 & ) * & ( w(k+1,j,i) + w(k_mm,j,i) ) & + ( ibit35 * adv_mom_5 & ) * & ( w(k_pp,j,i) + w(k_mmm,j,i) ) & ) diss_d = - ABS( w_comp ) * ( & ( 10.0_wp * ibit35 * adv_mom_5 & + 3.0_wp * ibit34 * adv_mom_3 & + ibit33 * adv_mom_1 & ) * & ( w(k,j,i) - w(k-1,j,i) ) & - ( 5.0_wp * ibit35 * adv_mom_5 & + ibit34 * adv_mom_3 & ) * & ( w(k+1,j,i) - w(k_mm,j,i) ) & + ( ibit35 * adv_mom_5 & ) * & ( w(k_pp,j,i) - w(k_mmm,j,i) ) & ) ! !-- k index has to be modified near bottom and top, else array !-- subscripts will be exceeded. ibit35 = IBITS(wall_flags_00(k,j,i),3,1) ibit34 = IBITS(wall_flags_00(k,j,i),2,1) ibit33 = IBITS(wall_flags_00(k,j,i),1,1) k_ppp = k + 3 * ibit35 k_pp = k + 2 * ( 1 - ibit33 ) k_mm = k - 2 * ibit35 w_comp = w(k+1,j,i) + w(k,j,i) flux_t = w_comp * ( & ( 37.0_wp * ibit35 * adv_mom_5 & + 7.0_wp * ibit34 * adv_mom_3 & + ibit33 * adv_mom_1 & ) * & ( w(k+1,j,i) + w(k,j,i) ) & - ( 8.0_wp * ibit35 * adv_mom_5 & + ibit34 * adv_mom_3 & ) * & ( w(k_pp,j,i) + w(k-1,j,i) ) & + ( ibit35 * adv_mom_5 & ) * & ( w(k_ppp,j,i) + w(k_mm,j,i) ) & ) diss_t = - ABS( w_comp ) * ( & ( 10.0_wp * ibit35 * adv_mom_5 & + 3.0_wp * ibit34 * adv_mom_3 & + ibit33 * adv_mom_1 & ) * & ( w(k+1,j,i) - w(k,j,i) ) & - ( 5.0_wp * ibit35 * adv_mom_5 & + ibit34 * adv_mom_3 & ) * & ( w(k_pp,j,i) - w(k-1,j,i) ) & + ( ibit35 * adv_mom_5 & ) * & ( w(k_ppp,j,i) - w(k_mm,j,i) ) & ) ! !-- Calculate the divergence of the velocity field. A respective !-- correction is needed to overcome numerical instabilities caused !-- by a not sufficient reduction of divergences near topography. div = ( ( u_comp + gu - ( u(k+1,j,i) + u(k,j,i) ) ) * ddx & + ( v_comp + gv - ( v(k+1,j,i) + v(k,j,i) ) ) * ddy & + ( w_comp - ( w(k,j,i) + w(k-1,j,i) ) ) & * ddzu(k+1) & ) * 0.5_wp tend(k,j,i) = - ( & ( flux_r + diss_r - flux_l - diss_l ) * ddx & + ( flux_n + diss_n - flux_s - diss_s ) * ddy & + ( flux_t + diss_t - flux_d - diss_d ) * ddzu(k+1) & ) + div * w(k,j,i) !++ !-- Statistical Evaluation of w'w'. ! sums_ws2_ws_l(k,tn) = sums_ws2_ws_l(k,tn) & ! + ( flux_t + diss_t ) & ! * weight_substep(intermediate_timestep_count) ENDDO ENDDO ENDDO !$acc end kernels END SUBROUTINE advec_w_ws_acc END MODULE advec_ws