MODULE advec_ws !-----------------------------------------------------------------------------! ! Current revisions: ! ------------------ ! ! Former revisions: ! ----------------- ! $Id: advec_ws.f90 889 2012-04-20 15:06:33Z maronga $ ! ! 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. ! ! 743 2011-08-18 16:10:16Z suehring ! Evaluation of turbulent fluxes with WS-scheme only for the whole model ! domain. Therefor dimension of arrays needed for statistical evaluation ! decreased. ! ! 736 2011-08-17 14:13:26Z suehring ! Bugfix concerning OpenMP parallelization. i_omp introduced, because first ! index where fluxes on left side have to be calculated explicitly is ! different on each thread. Furthermore the swapping of fluxes is now ! thread-safe by adding an additional dimension. ! ! 713 2011-03-30 14:21:21Z suehring ! File reformatted. ! Bugfix in vertical advection of w concerning the optimized version for ! vector architecture. ! Constants adv_mom_3, adv_mom_5, adv_sca_5, adv_sca_3 reformulated as ! broken numbers. ! ! 709 2011-03-30 09:31:40Z raasch ! formatting adjustments ! ! 705 2011-03-25 11:21:43 Z suehring $ ! Bugfix in declaration of logicals concerning outflow boundaries. ! ! 411 2009-12-11 12:31:43 Z suehring ! Allocation of weight_substep moved to init_3d_model. ! Declaration of ws_scheme_sca and ws_scheme_mom moved to check_parameters. ! Setting bc for the horizontal velocity variances added (moved from ! flow_statistics). ! ! 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_u_ws, advec_v_ws, advec_w_ws, & 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_v_ws MODULE PROCEDURE advec_v_ws MODULE PROCEDURE advec_v_ws_ij END INTERFACE advec_v_ws INTERFACE advec_w_ws MODULE PROCEDURE advec_w_ws MODULE PROCEDURE advec_w_ws_ij END INTERFACE advec_w_ws CONTAINS !------------------------------------------------------------------------------! ! Initialization of WS-scheme !------------------------------------------------------------------------------! SUBROUTINE ws_init USE arrays_3d USE constants USE control_parameters USE indices USE pegrid USE statistics ! !-- Set the appropriate factors for scalar and momentum advection. adv_sca_5 = 1./60. adv_sca_3 = 1./12. adv_sca_1 = 1./2. adv_mom_5 = 1./120. adv_mom_3 = 1./24. adv_mom_1 = 1./4. ! !-- 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 sums_wsvs_ws_l = 0.0 sums_us2_ws_l = 0.0 sums_vs2_ws_l = 0.0 sums_ws2_ws_l = 0.0 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 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 ENDIF IF ( ocean ) THEN ALLOCATE( sums_wssas_ws_l(nzb:nzt+1,0:threads_per_task-1) ) sums_wssas_ws_l = 0.0 ENDIF ENDIF ! !-- Arrays needed for reasons of speed optimization for cache and noopt !-- 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 ( 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 qauntities (fluxes, variances) !------------------------------------------------------------------------------! SUBROUTINE ws_statistics USE control_parameters USE statistics 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 sums_wsvs_ws_l = 0.0 sums_us2_ws_l = 0.0 sums_vs2_ws_l = 0.0 sums_ws2_ws_l = 0.0 ENDIF IF ( ws_scheme_sca ) THEN sums_wspts_ws_l = 0.0 IF ( humidity .OR. passive_scalar ) sums_wsqs_ws_l = 0.0 IF ( ocean ) sums_wssas_ws_l = 0.0 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 USE constants USE control_parameters USE grid_variables USE indices USE pegrid USE statistics IMPLICIT NONE INTEGER :: i, ibit0, ibit1, ibit2, ibit3, ibit4, ibit5, ibit6, & ibit7, ibit8, i_omp, j, k, k_mm, k_pp, k_ppp, tn REAL :: diss_d, div, flux_d, u_comp, v_comp REAL, DIMENSION(:,:,:), POINTER :: sk REAL, DIMENSION(nzb:nzt+1) :: diss_n, diss_r, diss_t, flux_n, & flux_r, flux_t REAL, DIMENSION(nzb+1:nzt,0:threads_per_task-1) :: swap_diss_y_local, & swap_flux_y_local REAL, DIMENSION(nzb+1:nzt,nys:nyn,0:threads_per_task-1) :: & swap_diss_x_local, & swap_flux_x_local CHARACTER (LEN = *), INTENT(IN) :: sk_char ! !-- 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 * ibit5 * adv_sca_5 & + 7.0 * ibit4 * adv_sca_3 & + ibit3 * adv_sca_1 & ) * & ( sk(k,j,i) + sk(k,j-1,i) ) & - ( 8.0 * 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 * ibit5 * adv_sca_5 & + 3.0 * ibit4 * adv_sca_3 & + ibit3 * adv_sca_1 & ) * & ( sk(k,j,i) - sk(k,j-1,i) ) & - ( 5.0 * 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 * ( sk(k,j,i) + sk(k,j-1,i) ) & - 8.0 * ( 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 * ( sk(k,j,i) - sk(k,j-1,i) ) & - 5.0 * ( 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 * ibit2 * adv_sca_5 & + 7.0 * ibit1 * adv_sca_3 & + ibit0 * adv_sca_1 & ) * & ( sk(k,j,i) + sk(k,j,i-1) ) & - ( 8.0 * 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 * ibit2 * adv_sca_5 & + 3.0 * ibit1 * adv_sca_3 & + ibit0 * adv_sca_1 & ) * & ( sk(k,j,i) - sk(k,j,i-1) ) & - ( 5.0 * 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 * ( sk(k,j,i) + sk(k,j,i-1) ) & - 8.0 * ( 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 * ( sk(k,j,i) - sk(k,j,i-1) ) & - 5.0 * ( 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 diss_t(0) = 0.0 flux_d = 0.0 diss_d = 0.0 ! !-- 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 * ibit2 * adv_sca_5 & + 7.0 * ibit1 * adv_sca_3 & + ibit0 * adv_sca_1 & ) * & ( sk(k,j,i+1) + sk(k,j,i) ) & - ( 8.0 * 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 * ibit2 * adv_sca_5 & + 3.0 * ibit1 * adv_sca_3 & + ibit0 * adv_sca_1 & ) * & ( sk(k,j,i+1) - sk(k,j,i) ) & - ( 5.0 * 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 * ibit5 * adv_sca_5 & + 7.0 * ibit4 * adv_sca_3 & + ibit3 * adv_sca_1 & ) * & ( sk(k,j+1,i) + sk(k,j,i) ) & - ( 8.0 * 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 * ibit5 * adv_sca_5 & + 3.0 * ibit4 * adv_sca_3 & + ibit3 * adv_sca_1 & ) * & ( sk(k,j+1,i) - sk(k,j,i) ) & - ( 5.0 * 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 * ibit8 * adv_sca_5 & + 7.0 * ibit7 * adv_sca_3 & + ibit6 * adv_sca_1 & ) * & ( sk(k+1,j,i) + sk(k,j,i) ) & - ( 8.0 * 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 * ibit8 * adv_sca_5 & + 3.0 * ibit7 * adv_sca_3 & + ibit6 * adv_sca_1 & ) * & ( sk(k+1,j,i) - sk(k,j,i) ) & - ( 5.0 * 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 * ( sk(k,j,i+1) + sk(k,j,i) ) & - 8.0 * ( 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 * ( sk(k,j,i+1) - sk(k,j,i) ) & - 5.0 * ( 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 * ( sk(k,j+1,i) + sk(k,j,i) ) & - 8.0 * ( 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 * ( sk(k,j+1,i) - sk(k,j,i) ) & - 5.0 * ( 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 * ibit8 * adv_sca_5 & + 7.0 * ibit7 * adv_sca_3 & + ibit6 * adv_sca_1 & ) * & ( sk(k+1,j,i) + sk(k,j,i) ) & - ( 8.0 * 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 * ibit8 * adv_sca_5 & + 3.0 * ibit7 * adv_sca_3 & + ibit6 * adv_sca_1 & ) * & ( sk(k+1,j,i) - sk(k,j,i) ) & - ( 5.0 * 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 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 USE constants USE control_parameters USE grid_variables USE indices USE statistics IMPLICIT NONE INTEGER :: i, ibit9, ibit10, ibit11, ibit12, ibit13, ibit14, ibit15, & ibit16, ibit17, i_omp, j, k, k_mm, k_pp, k_ppp, tn REAL :: diss_d, div, flux_d, gu, gv, u_comp_l, v_comp, w_comp REAL, DIMENSION(nzb:nzt+1) :: diss_n, diss_r, diss_t, flux_n, flux_r, & flux_t, u_comp gu = 2.0 * u_gtrans gv = 2.0 * 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 * ibit14 * adv_mom_5 & + 7.0 * ibit13 * adv_mom_3 & + ibit12 * adv_mom_1 & ) * & ( u(k,j,i) + u(k,j-1,i) ) & - ( 8.0 * 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 * ibit14 * adv_mom_5 & + 3.0 * ibit13 * adv_mom_3 & + ibit12 * adv_mom_1 & ) * & ( u(k,j,i) - u(k,j-1,i) ) & - ( 5.0 * 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 * ( u(k,j,i) + u(k,j-1,i) ) & - 8.0 * ( 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 * ( u(k,j,i) - u(k,j-1,i) ) & - 5.0 * ( 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 * ibit11 * adv_mom_5 & + 7.0 * ibit10 * adv_mom_3 & + ibit9 * adv_mom_1 & ) * & ( u(k,j,i) + u(k,j,i-1) ) & - ( 8.0 * 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 * ibit11 * adv_mom_5 & + 3.0 * ibit10 * adv_mom_3 & + ibit9 * adv_mom_1 & ) * & ( u(k,j,i) - u(k,j,i-1) ) & - ( 5.0 * 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 * ( u(k,j,i) + u(k,j,i-1) ) & - 8.0 * ( 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 * ( u(k,j,i) - u(k,j,i-1) ) & - 5.0 * ( 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 diss_t(0) = 0.0 flux_d = 0.0 diss_d = 0.0 ! !-- 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 * ibit11 * adv_mom_5 & + 7.0 * ibit10 * adv_mom_3 & + ibit9 * adv_mom_1 & ) * & ( u(k,j,i+1) + u(k,j,i) ) & - ( 8.0 * 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 * ibit11 * adv_mom_5 & + 3.0 * ibit10 * adv_mom_3 & + ibit9 * adv_mom_1 & ) * & ( u(k,j,i+1) - u(k,j,i) ) & - ( 5.0 * 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 * ibit14 * adv_mom_5 & + 7.0 * ibit13 * adv_mom_3 & + ibit12 * adv_mom_1 & ) * & ( u(k,j+1,i) + u(k,j,i) ) & - ( 8.0 * 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 * ibit14 * adv_mom_5 & + 3.0 * ibit13 * adv_mom_3 & + ibit12 * adv_mom_1 & ) * & ( u(k,j+1,i) - u(k,j,i) ) & - ( 5.0 * 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 * ibit17 * adv_mom_5 & + 7.0 * ibit16 * adv_mom_3 & + ibit15 * adv_mom_1 & ) * & ( u(k+1,j,i) + u(k,j,i) ) & - ( 8.0 * 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 * ibit17 * adv_mom_5 & + 3.0 * ibit16 * adv_mom_3 & + ibit15 * adv_mom_1 & ) * & ( u(k+1,j,i) - u(k,j,i) ) & - ( 5.0 * 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 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 * hom(k,1,1,0) ) & / ( u_comp(k) - gu + 1.0E-20 ) & + diss_r(k) * & ABS( u_comp(k) - 2.0 * hom(k,1,1,0) ) & / ( ABS( u_comp(k) - gu ) + 1.0E-20 ) ) & * 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 * ( u(k,j,i+1) + u(k,j,i) ) & - 8.0 * ( 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 * ( u(k,j,i+1) - u(k,j,i) ) & - 5.0 * ( 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 * ( u(k,j+1,i) + u(k,j,i) ) & - 8.0 * ( 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 * ( u(k,j+1,i) - u(k,j,i) ) & - 5.0 * ( 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 * ibit17 * adv_mom_5 & + 7.0 * ibit16 * adv_mom_3 & + ibit15 * adv_mom_1 & ) * & ( u(k+1,j,i) + u(k,j,i) ) & - ( 8.0 * 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 * ibit17 * adv_mom_5 & + 3.0 * ibit16 * adv_mom_3 & + ibit15 * adv_mom_1 & ) * & ( u(k+1,j,i) - u(k,j,i) ) & - ( 5.0 * 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 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 * hom(k,1,1,0) ) & / ( u_comp(k) - gu + 1.0E-20 ) & + diss_r(k) * & ABS( u_comp(k) - 2.0 * hom(k,1,1,0) ) & / ( ABS( u_comp(k) - gu ) + 1.0E-20 ) ) & * 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 USE constants USE control_parameters USE grid_variables USE indices USE statistics IMPLICIT NONE INTEGER :: i, ibit18, ibit19, ibit20, ibit21, ibit22, ibit23, ibit24, & ibit25, ibit26, i_omp, j, k, k_mm, k_pp, k_ppp, tn REAL :: diss_d, div, flux_d, gu, gv, u_comp, v_comp_l, w_comp REAL, DIMENSION(nzb:nzt+1) :: diss_n, diss_r, diss_t, flux_n, flux_r, & flux_t, v_comp gu = 2.0 * u_gtrans gv = 2.0 * 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 * ibit20 * adv_mom_5 & + 7.0 * ibit19 * adv_mom_3 & + ibit18 * adv_mom_1 & ) * & ( v(k,j,i) + v(k,j,i-1) ) & - ( 8.0 * 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 * ibit20 * adv_mom_5 & + 3.0 * ibit19 * adv_mom_3 & + ibit18 * adv_mom_1 & ) * & ( v(k,j,i) - v(k,j,i-1) ) & - ( 5.0 * 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 * ( v(k,j,i) + v(k,j,i-1) ) & - 8.0 * ( 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 * ( v(k,j,i) - v(k,j,i-1) ) & - 5.0 * ( 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 * ibit23 * adv_mom_5 & + 7.0 * ibit22 * adv_mom_3 & + ibit21 * adv_mom_1 & ) * & ( v(k,j,i) + v(k,j-1,i) ) & - ( 8.0 * 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 * ibit23 * adv_mom_5 & + 3.0 * ibit22 * adv_mom_3 & + ibit21 * adv_mom_1 & ) * & ( v(k,j,i) - v(k,j-1,i) ) & - ( 5.0 * 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 * ( v(k,j,i) + v(k,j-1,i) ) & - 8.0 * ( 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 * ( v(k,j,i) - v(k,j-1,i) ) & - 5.0 * ( 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 diss_t(0) = 0.0 flux_d = 0.0 diss_d = 0.0 ! !-- 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 * ibit20 * adv_mom_5 & + 7.0 * ibit19 * adv_mom_3 & + ibit18 * adv_mom_1 & ) * & ( v(k,j,i+1) + v(k,j,i) ) & - ( 8.0 * 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 * ibit20 * adv_mom_5 & + 3.0 * ibit19 * adv_mom_3 & + ibit18 * adv_mom_1 & ) * & ( v(k,j,i+1) - v(k,j,i) ) & - ( 5.0 * 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 * ibit23 * adv_mom_5 & + 7.0 * ibit22 * adv_mom_3 & + ibit21 * adv_mom_1 & ) * & ( v(k,j+1,i) + v(k,j,i) ) & - ( 8.0 * 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 * ibit23 * adv_mom_5 & + 3.0 * ibit22 * adv_mom_3 & + ibit21 * adv_mom_1 & ) * & ( v(k,j+1,i) - v(k,j,i) ) & - ( 5.0 * 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 * ibit26 * adv_mom_5 & + 7.0 * ibit25 * adv_mom_3 & + ibit24 * adv_mom_1 & ) * & ( v(k+1,j,i) + v(k,j,i) ) & - ( 8.0 * 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 * ibit26 * adv_mom_5 & + 3.0 * ibit25 * adv_mom_3 & + ibit24 * adv_mom_1 & ) * & ( v(k+1,j,i) - v(k,j,i) ) & - ( 5.0 * 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 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 * hom(k,1,2,0) ) & / ( v_comp(k) - gv + 1.0E-20 ) & + diss_n(k) & * ABS( v_comp(k) - 2.0 * hom(k,1,2,0) ) & / ( ABS( v_comp(k) - gv ) +1.0E-20 ) ) & * 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 * ( v(k,j,i+1) + v(k,j,i) ) & - 8.0 * ( 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 * ( v(k,j,i+1) - v(k,j,i) ) & - 5.0 * ( 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 * ( v(k,j+1,i) + v(k,j,i) ) & - 8.0 * ( 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 * ( v(k,j+1,i) - v(k,j,i) ) & - 5.0 * ( 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 * ibit26 * adv_mom_5 & + 7.0 * ibit25 * adv_mom_3 & + ibit24 * adv_mom_1 & ) * & ( v(k+1,j,i) + v(k,j,i) ) & - ( 8.0 * 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 * ibit26 * adv_mom_5 & + 3.0 * ibit25 * adv_mom_3 & + ibit24 * adv_mom_1 & ) * & ( v(k+1,j,i) - v(k,j,i) ) & - ( 5.0 * 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 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 * hom(k,1,2,0) ) & / ( v_comp(k) - gv + 1.0E-20 ) & + diss_n(k) & * ABS( v_comp(k) - 2.0 * hom(k,1,2,0) ) & / ( ABS( v_comp(k) - gv ) +1.0E-20 ) ) & * 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 USE constants USE control_parameters USE grid_variables USE indices USE statistics IMPLICIT NONE INTEGER :: i, ibit27, ibit28, ibit29, ibit30, ibit31, ibit32, ibit33, & ibit34, ibit35, i_omp, j, k, k_mm, k_pp, k_ppp, tn REAL :: diss_d, div, flux_d, gu, gv, u_comp, v_comp, w_comp REAL, DIMENSION(nzb:nzt+1) :: diss_n, diss_r, diss_t, flux_n, flux_r, & flux_t gu = 2.0 * u_gtrans gv = 2.0 * v_gtrans ! !-- Compute southside fluxes for the respective boundary. IF ( j == nys ) THEN DO k = nzb+1, nzb_max ibit32 = IBITS(wall_flags_0(k,j,i),32,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 * ibit32 * adv_mom_5 & + 7.0 * ibit31 * adv_mom_3 & + ibit30 * adv_mom_1 & ) * & ( w(k,j,i) + w(k,j-1,i) ) & - ( 8.0 * 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 * ibit32 * adv_mom_5 & + 3.0 * ibit31 * adv_mom_3 & + ibit30 * adv_mom_1 & ) * & ( w(k,j,i) - w(k,j-1,i) ) & - ( 5.0 * 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 * ( w(k,j,i) + w(k,j-1,i) ) & - 8.0 * ( 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 * ( w(k,j,i) - w(k,j-1,i) ) & - 5.0 * ( 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 * ibit29 * adv_mom_5 & + 7.0 * ibit28 * adv_mom_3 & + ibit27 * adv_mom_1 & ) * & ( w(k,j,i) + w(k,j,i-1) ) & - ( 8.0 * 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 * ibit29 * adv_mom_5 & + 3.0 * ibit28 * adv_mom_3 & + ibit27 * adv_mom_1 & ) * & ( w(k,j,i) - w(k,j,i-1) ) & - ( 5.0 * 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 * ( w(k,j,i) + w(k,j,i-1) ) & - 8.0 * ( 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 * ( w(k,j,i) - w(k,j,i-1) ) & - 5.0 * ( 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 * ibit29 * adv_mom_5 & + 7.0 * ibit28 * adv_mom_3 & + ibit27 * adv_mom_1 & ) * & ( w(k,j,i+1) + w(k,j,i) ) & - ( 8.0 * 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 * ibit29 * adv_mom_5 & + 3.0 * ibit28 * adv_mom_3 & + ibit27 * adv_mom_1 & ) * & ( w(k,j,i+1) - w(k,j,i) ) & - ( 5.0 * 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_0(k,j,i),32,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 * ibit32 * adv_mom_5 & + 7.0 * ibit31 * adv_mom_3 & + ibit30 * adv_mom_1 & ) * & ( w(k,j+1,i) + w(k,j,i) ) & - ( 8.0 * 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 * ibit32 * adv_mom_5 & + 3.0 * ibit31 * adv_mom_3 & + ibit30 * adv_mom_1 & ) * & ( w(k,j+1,i) - w(k,j,i) ) & - ( 5.0 * 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_0(k,j,i),35,1) ibit34 = IBITS(wall_flags_0(k,j,i),34,1) ibit33 = IBITS(wall_flags_0(k,j,i),33,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 * ibit35 * adv_mom_5 & + 7.0 * ibit34 * adv_mom_3 & + ibit33 * adv_mom_1 & ) * & ( w(k+1,j,i) + w(k,j,i) ) & - ( 8.0 * 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 * ibit35 * adv_mom_5 & + 3.0 * ibit34 * adv_mom_3 & + ibit33 * adv_mom_1 & ) * & ( w(k+1,j,i) - w(k,j,i) ) & - ( 5.0 * 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 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 * ( w(k,j,i+1) + w(k,j,i) ) & - 8.0 * ( 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 * ( w(k,j,i+1) - w(k,j,i) ) & - 5.0 * ( 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 * ( w(k,j+1,i) + w(k,j,i) ) & - 8.0 * ( 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 * ( w(k,j+1,i) - w(k,j,i) ) & - 5.0 * ( 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_0(k,j,i),35,1) ibit34 = IBITS(wall_flags_0(k,j,i),34,1) ibit33 = IBITS(wall_flags_0(k,j,i),33,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 * ibit35 * adv_mom_5 & + 7.0 * ibit34 * adv_mom_3 & + ibit33 * adv_mom_1 & ) * & ( w(k+1,j,i) + w(k,j,i) ) & - ( 8.0 * 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 * ibit35 * adv_mom_5 & + 3.0 * ibit34 * adv_mom_3 & + ibit33 * adv_mom_1 & ) * & ( w(k+1,j,i) - w(k,j,i) ) & - ( 5.0 * 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 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 USE constants USE control_parameters USE grid_variables USE indices USE statistics IMPLICIT NONE INTEGER :: i, ibit0, ibit1, ibit2, ibit3, ibit4, ibit5, ibit6, & ibit7, ibit8, j, k, k_mm, k_pp, k_ppp, tn = 0 REAL, DIMENSION(:,:,:), POINTER :: sk REAL :: diss_d, div, flux_d, u_comp, v_comp REAL, DIMENSION(nzb:nzt) :: diss_n, diss_r, diss_t, flux_n, flux_r, & flux_t REAL, DIMENSION(nzb+1:nzt) :: swap_diss_y_local, swap_flux_y_local REAL, DIMENSION(nzb+1:nzt,nys:nyn) :: swap_diss_x_local, & swap_flux_x_local CHARACTER (LEN = *), INTENT(IN) :: sk_char ! !-- 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 * ibit2 * adv_sca_5 & + 7.0 * ibit1 * adv_sca_3 & + ibit0 * adv_sca_1 & ) * & ( sk(k,j,i) + sk(k,j,i-1) ) & - ( 8.0 * 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 * ibit2 * adv_sca_5 & + 3.0 * ibit1 * adv_sca_3 & + ibit0 * adv_sca_1 & ) * & ( sk(k,j,i) - sk(k,j,i-1) ) & - ( 5.0 * 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 * ( sk(k,j,i) + sk(k,j,i-1) ) & - 8.0 * ( 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 * ( sk(k,j,i) - sk(k,j,i-1) ) & - 5.0 * ( 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 * ibit5 * adv_sca_5 & + 7.0 * ibit4 * adv_sca_3 & + ibit3 * adv_sca_1 & ) * & ( sk(k,j,i) + sk(k,j-1,i) ) & - ( 8.0 * 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 * ibit5 * adv_sca_5 & + 3.0 * ibit4 * adv_sca_3 & + ibit3 * adv_sca_1 & ) * & ( sk(k,j,i) - sk(k,j-1,i) ) & - ( 5.0 * 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 * ( sk(k,j,i) + sk(k,j-1,i) ) & - 8.0 * ( 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 * ( sk(k,j,i) - sk(k,j-1,i) ) & - 5.0 * ( 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 diss_t(0) = 0.0 flux_d = 0.0 diss_d = 0.0 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 * ibit2 * adv_sca_5 & + 7.0 * ibit1 * adv_sca_3 & + ibit0 * adv_sca_1 & ) * & ( sk(k,j,i+1) + sk(k,j,i) ) & - ( 8.0 * 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 * ibit2 * adv_sca_5 & + 3.0 * ibit1 * adv_sca_3 & + ibit0 * adv_sca_1 & ) * & ( sk(k,j,i+1) - sk(k,j,i) ) & - ( 5.0 * 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 * ibit5 * adv_sca_5 & + 7.0 * ibit4 * adv_sca_3 & + ibit3 * adv_sca_1 & ) * & ( sk(k,j+1,i) + sk(k,j,i) ) & - ( 8.0 * 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 * ibit5 * adv_sca_5 & + 3.0 * ibit4 * adv_sca_3 & + ibit3 * adv_sca_1 & ) * & ( sk(k,j+1,i) - sk(k,j,i) ) & - ( 5.0 * 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 * ibit8 * adv_sca_5 & + 7.0 * ibit7 * adv_sca_3 & + ibit6 * adv_sca_1 & ) * & ( sk(k+1,j,i) + sk(k,j,i) ) & - ( 8.0 * 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 * ibit8 * adv_sca_5 & + 3.0 * ibit7 * adv_sca_3 & + ibit6 * adv_sca_1 & ) * & ( sk(k+1,j,i) - sk(k,j,i) ) & - ( 5.0 * 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 * ( sk(k,j,i+1) + sk(k,j,i) ) & - 8.0 * ( 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 * ( sk(k,j,i+1) - sk(k,j,i) ) & - 5.0 * ( 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 * ( sk(k,j+1,i) + sk(k,j,i) ) & - 8.0 * ( 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 * ( sk(k,j+1,i) - sk(k,j,i) ) & - 5.0 * ( 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 * ibit8 * adv_sca_5 & + 7.0 * ibit7 * adv_sca_3 & + ibit6 * adv_sca_1 & ) * & ( sk(k+1,j,i) + sk(k,j,i) ) & - ( 8.0 * 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 * ibit8 * adv_sca_5 & + 3.0 * ibit7 * adv_sca_3 & + ibit6 * adv_sca_1 & ) * & ( sk(k+1,j,i) - sk(k,j,i) ) & - ( 5.0 * 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 END SELECT ENDDO ENDDO END SUBROUTINE advec_s_ws !------------------------------------------------------------------------------! ! Advection of u - Call for all grid points !------------------------------------------------------------------------------! SUBROUTINE advec_u_ws USE arrays_3d USE constants USE control_parameters USE grid_variables USE indices USE statistics IMPLICIT NONE INTEGER :: i, ibit9, ibit10, ibit11, ibit12, ibit13, ibit14, ibit15, & ibit16, ibit17, j, k, k_mm, k_pp, k_ppp, tn = 0 REAL :: diss_d, div, flux_d, gu, gv, v_comp, w_comp REAL, DIMENSION(nzb+1:nzt) :: swap_diss_y_local_u, swap_flux_y_local_u REAL, DIMENSION(nzb+1:nzt,nys:nyn) :: swap_diss_x_local_u, & swap_flux_x_local_u REAL, DIMENSION(nzb:nzt) :: diss_n, diss_r, diss_t, flux_n, flux_r, & flux_t, u_comp gu = 2.0 * u_gtrans gv = 2.0 * 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 * ibit11 * adv_mom_5 & + 7.0 * ibit10 * adv_mom_3 & + ibit9 * adv_mom_1 & ) * & ( u(k,j,i) + u(k,j,i-1) ) & - ( 8.0 * 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 * ibit11 * adv_mom_5 & + 3.0 * ibit10 * adv_mom_3 & + ibit9 * adv_mom_1 & ) * & ( u(k,j,i) - u(k,j,i-1) ) & - ( 5.0 * 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 * ( u(k,j,i) + u(k,j,i-1) ) & - 8.0 * ( 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 * ( u(k,j,i) - u(k,j,i-1) ) & - 5.0 * ( 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 * ibit14 * adv_mom_5 & + 7.0 * ibit13 * adv_mom_3 & + ibit12 * adv_mom_1 & ) * & ( u(k,j,i) + u(k,j-1,i) ) & - ( 8.0 * 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 * ibit14 * adv_mom_5 & + 3.0 * ibit13 * adv_mom_3 & + ibit12 * adv_mom_1 & ) * & ( u(k,j,i) - u(k,j-1,i) ) & - ( 5.0 * 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 * ( u(k,j,i) + u(k,j-1,i) ) & - 8.0 * ( 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 * ( u(k,j,i) - u(k,j-1,i) ) & - 5.0 * ( 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 diss_t(0) = 0.0 flux_d = 0.0 diss_d = 0.0 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 * ibit11 * adv_mom_5 & + 7.0 * ibit10 * adv_mom_3 & + ibit9 * adv_mom_1 & ) * & ( u(k,j,i+1) + u(k,j,i) ) & - ( 8.0 * 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 * ibit11 * adv_mom_5 & + 3.0 * ibit10 * adv_mom_3 & + ibit9 * adv_mom_1 & ) * & ( u(k,j,i+1) - u(k,j,i) ) & - ( 5.0 * 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 * ibit14 * adv_mom_5 & + 7.0 * ibit13 * adv_mom_3 & + ibit12 * adv_mom_1 & ) * & ( u(k,j+1,i) + u(k,j,i) ) & - ( 8.0 * 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 * ibit14 * adv_mom_5 & + 3.0 * ibit13 * adv_mom_3 & + ibit12 * adv_mom_1 & ) * & ( u(k,j+1,i) - u(k,j,i) ) & - ( 5.0 * 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 * ibit17 * adv_mom_5 & + 7.0 * ibit16 * adv_mom_3 & + ibit15 * adv_mom_1 & ) * & ( u(k+1,j,i) + u(k,j,i) ) & - ( 8.0 * 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 * ibit17 * adv_mom_5 & + 3.0 * ibit16 * adv_mom_3 & + ibit15 * adv_mom_1 & ) * & ( u(k+1,j,i) - u(k,j,i) ) & - ( 5.0 * 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 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 * hom(k,1,1,0) ) & / ( u_comp(k) - gu + 1.0E-20 ) & + diss_r(k) * & ABS( u_comp(k) - 2.0 * hom(k,1,1,0) ) & / ( ABS( u_comp(k) - gu ) + 1.0E-20 ) ) & * 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 * ( u(k,j,i+1) + u(k,j,i) ) & - 8.0 * ( 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 * ( u(k,j,i+1) - u(k,j,i) ) & - 5.0 * ( 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 * ( u(k,j+1,i) + u(k,j,i) ) & - 8.0 * ( 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 * ( u(k,j+1,i) - u(k,j,i) ) & - 5.0 * ( 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 * ibit17 * adv_mom_5 & + 7.0 * ibit16 * adv_mom_3 & + ibit15 * adv_mom_1 & ) * & ( u(k+1,j,i) + u(k,j,i) ) & - ( 8.0 * 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 * ibit17 * adv_mom_5 & + 3.0 * ibit16 * adv_mom_3 & + ibit15 * adv_mom_1 & ) * & ( u(k+1,j,i) - u(k,j,i) ) & - ( 5.0 * 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 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 * hom(k,1,1,0) ) & / ( u_comp(k) - gu + 1.0E-20 ) & + diss_r(k) * & ABS( u_comp(k) - 2.0 * hom(k,1,1,0) ) & / ( ABS( u_comp(k) - gu ) + 1.0E-20 ) ) & * 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 v - Call for all grid points !------------------------------------------------------------------------------! SUBROUTINE advec_v_ws USE arrays_3d USE constants USE control_parameters USE grid_variables USE indices USE statistics IMPLICIT NONE INTEGER :: i, ibit18, ibit19, ibit20, ibit21, ibit22, ibit23, ibit24, & ibit25, ibit26, j, k, k_mm, k_pp, k_ppp, tn = 0 REAL :: diss_d, div, flux_d, gu, gv, u_comp, w_comp REAL, DIMENSION(nzb+1:nzt) :: swap_diss_y_local_v, swap_flux_y_local_v REAL, DIMENSION(nzb+1:nzt,nys:nyn) :: swap_diss_x_local_v, & swap_flux_x_local_v REAL, DIMENSION(nzb:nzt) :: diss_n, diss_r, diss_t, flux_n, flux_r, & flux_t, v_comp gu = 2.0 * u_gtrans gv = 2.0 * 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 * ibit20 * adv_mom_5 & + 7.0 * ibit19 * adv_mom_3 & + ibit18 * adv_mom_1 & ) * & ( v(k,j,i) + v(k,j,i-1) ) & - ( 8.0 * 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 * ibit20 * adv_mom_5 & + 3.0 * ibit19 * adv_mom_3 & + ibit18 * adv_mom_1 & ) * & ( v(k,j,i) - v(k,j,i-1) ) & - ( 5.0 * 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 * ( v(k,j,i) + v(k,j,i-1) ) & - 8.0 * ( 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 * ( v(k,j,i) - v(k,j,i-1) ) & - 5.0 * ( 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 * ibit23 * adv_mom_5 & + 7.0 * ibit22 * adv_mom_3 & + ibit21 * adv_mom_1 & ) * & ( v(k,j,i) + v(k,j-1,i) ) & - ( 8.0 * 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 * ibit23 * adv_mom_5 & + 3.0 * ibit22 * adv_mom_3 & + ibit21 * adv_mom_1 & ) * & ( v(k,j,i) - v(k,j-1,i) ) & - ( 5.0 * 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 * ( v(k,j,i) + v(k,j-1,i) ) & - 8.0 * ( 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 * ( v(k,j,i) - v(k,j-1,i) ) & - 5.0 * ( 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 diss_t(0) = 0.0 flux_d = 0.0 diss_d = 0.0 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 * ibit20 * adv_mom_5 & + 7.0 * ibit19 * adv_mom_3 & + ibit18 * adv_mom_1 & ) * & ( v(k,j,i+1) + v(k,j,i) ) & - ( 8.0 * 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 * ibit20 * adv_mom_5 & + 3.0 * ibit19 * adv_mom_3 & + ibit18 * adv_mom_1 & ) * & ( v(k,j,i+1) - v(k,j,i) ) & - ( 5.0 * 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 * ibit23 * adv_mom_5 & + 7.0 * ibit22 * adv_mom_3 & + ibit21 * adv_mom_1 & ) * & ( v(k,j+1,i) + v(k,j,i) ) & - ( 8.0 * 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 * ibit23 * adv_mom_5 & + 3.0 * ibit22 * adv_mom_3 & + ibit21 * adv_mom_1 & ) * & ( v(k,j+1,i) - v(k,j,i) ) & - ( 5.0 * 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 * ibit26 * adv_mom_5 & + 7.0 * ibit25 * adv_mom_3 & + ibit24 * adv_mom_1 & ) * & ( v(k+1,j,i) + v(k,j,i) ) & - ( 8.0 * 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 * ibit26 * adv_mom_5 & + 3.0 * ibit25 * adv_mom_3 & + ibit24 * adv_mom_1 & ) * & ( v(k+1,j,i) - v(k,j,i) ) & - ( 5.0 * 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 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 * hom(k,1,2,0) ) & / ( v_comp(k) - gv + 1.0E-20 ) & + diss_n(k) & * ABS( v_comp(k) - 2.0 * hom(k,1,2,0) ) & / ( ABS( v_comp(k) - gv ) +1.0E-20 ) ) & * 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 * ( v(k,j,i+1) + v(k,j,i) ) & - 8.0 * ( 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 * ( v(k,j,i+1) - v(k,j,i) ) & - 5.0 * ( 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 * ( v(k,j+1,i) + v(k,j,i) ) & - 8.0 * ( 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 * ( v(k,j+1,i) - v(k,j,i) ) & - 5.0 * ( 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 * ibit26 * adv_mom_5 & + 7.0 * ibit25 * adv_mom_3 & + ibit24 * adv_mom_1 & ) * & ( v(k+1,j,i) + v(k,j,i) ) & - ( 8.0 * 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 * ibit26 * adv_mom_5 & + 3.0 * ibit25 * adv_mom_3 & + ibit24 * adv_mom_1 & ) * & ( v(k+1,j,i) - v(k,j,i) ) & - ( 5.0 * 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 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 * hom(k,1,2,0) ) & / ( v_comp(k) - gv + 1.0E-20 ) & + diss_n(k) & * ABS( v_comp(k) - 2.0 * hom(k,1,2,0) ) & / ( ABS( v_comp(k) - gv ) +1.0E-20 ) ) & * 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 w - Call for all grid points !------------------------------------------------------------------------------! SUBROUTINE advec_w_ws USE arrays_3d USE constants USE control_parameters USE grid_variables USE indices USE statistics IMPLICIT NONE INTEGER :: i, ibit27, ibit28, ibit29, ibit30, ibit31, ibit32, ibit33, & ibit34, ibit35, j, k, k_mm, k_pp, k_ppp, tn = 0 REAL :: diss_d, div, flux_d, gu, gv, u_comp, v_comp, w_comp REAL, DIMENSION(nzb:nzt) :: diss_t, flux_t REAL, DIMENSION(nzb+1:nzt) :: diss_n, diss_r, flux_n, flux_r, & swap_diss_y_local_w, & swap_flux_y_local_w REAL, DIMENSION(nzb+1:nzt,nys:nyn) :: swap_diss_x_local_w, & swap_flux_x_local_w gu = 2.0 * u_gtrans gv = 2.0 * 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 * ibit29 * adv_mom_5 & + 7.0 * ibit28 * adv_mom_3 & + ibit27 * adv_mom_1 & ) * & ( w(k,j,i) + w(k,j,i-1) ) & - ( 8.0 * 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 * ibit29 * adv_mom_5 & + 3.0 * ibit28 * adv_mom_3 & + ibit27 * adv_mom_1 & ) * & ( w(k,j,i) - w(k,j,i-1) ) & - ( 5.0 * 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 * ( w(k,j,i) + w(k,j,i-1) ) & - 8.0 * ( 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 * ( w(k,j,i) - w(k,j,i-1) ) & - 5.0 * ( 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_0(k,j,i),32,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 * ibit32 * adv_mom_5 & + 7.0 * ibit31 * adv_mom_3 & + ibit30 * adv_mom_1 & ) * & ( w(k,j,i) + w(k,j-1,i) ) & - ( 8.0 * 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 * ibit32 * adv_mom_5 & + 3.0 * ibit31 * adv_mom_3 & + ibit30 * adv_mom_1 & ) * & ( w(k,j,i) - w(k,j-1,i) ) & - ( 5.0 * 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 * ( w(k,j,i) + w(k,j-1,i) ) & - 8.0 * ( 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 * ( w(k,j,i) - w(k,j-1,i) ) & - 5.0 * ( 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 * ibit29 * adv_mom_5 & + 7.0 * ibit28 * adv_mom_3 & + ibit27 * adv_mom_1 & ) * & ( w(k,j,i+1) + w(k,j,i) ) & - ( 8.0 * 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 * ibit29 * adv_mom_5 & + 3.0 * ibit28 * adv_mom_3 & + ibit27 * adv_mom_1 & ) * & ( w(k,j,i+1) - w(k,j,i) ) & - ( 5.0 * 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_0(k,j,i),32,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 * ibit32 * adv_mom_5 & + 7.0 * ibit31 * adv_mom_3 & + ibit30 * adv_mom_1 & ) * & ( w(k,j+1,i) + w(k,j,i) ) & - ( 8.0 * 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 * ibit32 * adv_mom_5 & + 3.0 * ibit31 * adv_mom_3 & + ibit30 * adv_mom_1 & ) * & ( w(k,j+1,i) - w(k,j,i) ) & - ( 5.0 * 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_0(k,j,i),35,1) ibit34 = IBITS(wall_flags_0(k,j,i),34,1) ibit33 = IBITS(wall_flags_0(k,j,i),33,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 * ibit35 * adv_mom_5 & + 7.0 * ibit34 * adv_mom_3 & + ibit33 * adv_mom_1 & ) * & ( w(k+1,j,i) + w(k,j,i) ) & - ( 8.0 * 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 * ibit35 * adv_mom_5 & + 3.0 * ibit34 * adv_mom_3 & + ibit33 * adv_mom_1 & ) * & ( w(k+1,j,i) - w(k,j,i) ) & - ( 5.0 * 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 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 * ( w(k,j,i+1) + w(k,j,i) ) & - 8.0 * ( 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 * ( w(k,j,i+1) - w(k,j,i) ) & - 5.0 * ( 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 * ( w(k,j+1,i) + w(k,j,i) ) & - 8.0 * ( 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 * ( w(k,j+1,i) - w(k,j,i) ) & - 5.0 * ( 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_0(k,j,i),35,1) ibit34 = IBITS(wall_flags_0(k,j,i),34,1) ibit33 = IBITS(wall_flags_0(k,j,i),33,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 * ibit35 * adv_mom_5 & + 7.0 * ibit34 * adv_mom_3 & + ibit33 * adv_mom_1 & ) * & ( w(k+1,j,i) + w(k,j,i) ) & - ( 8.0 * 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 * ibit35 * adv_mom_5 & + 3.0 * ibit34 * adv_mom_3 & + ibit33 * adv_mom_1 & ) * & ( w(k+1,j,i) - w(k,j,i) ) & - ( 5.0 * 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 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 END MODULE advec_ws