MODULE advec_ws !-----------------------------------------------------------------------------! ! Current revisions: ! ----------------- ! ! Former revisions: ! ----------------- ! $Id: advec_s_ws.f90 411 2009-12-11 12:31:43 Z suehring $ ! ! 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. ! In case of vector architectures Dirichlet and Radiation boundary conditions ! are outstanding and not available. ! A further routine local_diss_ij is available (next weeks) and is used if a ! control of dissipative fluxes is desired. ! In case of vertical grid stretching the order of advection scheme is ! degraded. This is also realized for the ocean where the stretching is ! downwards. ! ! OUTSTANDING: - Dirichlet and Radiation boundary conditions for ! vector architectures ! - dissipation control for cache architectures ( next weeks ) ! - Topography ( next weeks ) !-----------------------------------------------------------------------------! PRIVATE PUBLIC advec_s_ws, advec_u_ws, advec_v_ws, advec_w_ws, & local_diss, 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 INTERFACE local_diss MODULE PROCEDURE local_diss MODULE PROCEDURE local_diss_ij END INTERFACE local_diss CONTAINS ! !-- Initialisation SUBROUTINE ws_init USE arrays_3d USE constants USE control_parameters USE indices USE statistics ! !-- Allocate arrays needed for dissipation control. IF ( dissipation_control ) THEN ! ALLOCATE(var_x(nzb+1:nzt,nys:nyn,nxl:nxr), & ! var_y(nzb+1:nzt,nys:nyn,nxl:nxr), & ! var_z(nzb+1:nzt,nys:nyn,nxl:nxr), & ! gamma_x(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & ! gamma_y(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & ! gamma_z(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) ENDIF !-- Set the appropriate factors for scalar and momentum advection. adv_sca_5 = 0.016666666666666 adv_sca_3 = 0.083333333333333 adv_mom_5 = 0.0083333333333333 adv_mom_3 = 0.041666666666666 ! !-- Arrays needed for statical evaluation of fluxes. IF ( ws_scheme_mom ) THEN ALLOCATE( sums_wsus_ws_l(nzb:nzt+1,0:statistic_regions), & sums_wsvs_ws_l(nzb:nzt+1,0:statistic_regions), & sums_us2_ws_l(nzb:nzt+1,0:statistic_regions), & sums_vs2_ws_l(nzb:nzt+1,0:statistic_regions), & sums_ws2_ws_l(nzb:nzt+1,0:statistic_regions)) 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:statistic_regions)) sums_wspts_ws_l = 0.0 IF ( humidity .OR. passive_scalar ) THEN ALLOCATE( sums_wsqs_ws_l(nzb:nzt+1,0:statistic_regions)) sums_wsqs_ws_l = 0.0 ENDIF IF ( ocean ) THEN ALLOCATE( sums_wssas_ws_l(nzb:nzt+1,0:statistic_regions)) 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), flux_s_v(nzb+1:nzt), & flux_s_w(nzb+1:nzt), diss_s_u(nzb+1:nzt), & diss_s_v(nzb+1:nzt), diss_s_w(nzb+1:nzt)) ALLOCATE(flux_l_u(nzb+1:nzt,nys:nyn), & flux_l_v(nzb+1:nzt,nys:nyn), flux_l_w(nzb+1:nzt,nys:nyn), & diss_l_u(nzb+1:nzt,nys:nyn), diss_l_v(nzb+1:nzt,nys:nyn), & diss_l_w(nzb+1:nzt,nys:nyn)) ENDIF IF ( ws_scheme_sca ) THEN ALLOCATE(flux_s_pt(nzb+1:nzt), flux_s_e(nzb+1:nzt), & diss_s_pt(nzb+1:nzt), diss_s_e(nzb+1:nzt)) ALLOCATE(flux_l_pt(nzb+1:nzt,nys:nyn), & flux_l_e(nzb+1:nzt,nys:nyn), diss_l_pt(nzb+1:nzt,nys:nyn), & diss_l_e(nzb+1:nzt,nys:nyn)) IF ( humidity .OR. passive_scalar ) THEN ALLOCATE(flux_s_q(nzb+1:nzt), diss_s_q(nzb+1:nzt)) ALLOCATE(flux_l_q(nzb+1:nzt,nys:nyn), & diss_l_q(nzb+1:nzt,nys:nyn)) ENDIF IF ( ocean ) THEN ALLOCATE(flux_s_sa(nzb+1:nzt), diss_s_sa(nzb+1:nzt)) ALLOCATE(flux_l_sa(nzb+1:nzt,nys:nyn), & diss_l_sa(nzb+1:nzt,nys:nyn)) ENDIF ENDIF ENDIF ! !-- Determine the flags where the order of the scheme for horizontal !-- advection should be degraded. ALLOCATE( boundary_flags(nys:nyn,nxl:nxr) ) DO i = nxl, nxr DO j = nys, nyn boundary_flags(j,i) = 0 IF ( outflow_l ) THEN IF ( i == nxlu ) boundary_flags(j,i) = 5 IF ( i == nxl ) boundary_flags(j,i) = 6 ELSEIF ( outflow_r ) THEN IF ( i == nxr-1 ) boundary_flags(j,i) = 1 IF ( i == nxr ) boundary_flags(j,i) = 2 ELSEIF ( outflow_n ) THEN IF ( j == nyn-1 ) boundary_flags(j,i) = 3 IF ( j == nyn ) boundary_flags(j,i) = 4 ELSEIF ( outflow_s ) THEN IF ( j == nysv ) boundary_flags(j,i) = 7 IF ( j == nys ) boundary_flags(j,i) = 8 ENDIF ENDDO ENDDO END SUBROUTINE ws_init SUBROUTINE ws_statistics USE control_parameters USE statistics ! !-- The arrays needed for statistical evaluation are set to to 0 at the !-- begin 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 ) USE arrays_3d USE constants USE control_parameters USE grid_variables USE indices USE statistics IMPLICIT NONE INTEGER :: i, j, k LOGICAL :: degraded_l = .FALSE., degraded_s = .FALSE. REAL :: flux_d, diss_d, u_comp, v_comp REAL, DIMENSION(:,:,:), POINTER :: sk REAL, DIMENSION(nzb:nzt+1) :: flux_t, diss_t, flux_r, diss_r, & flux_n, diss_n REAL, DIMENSION(nzb+1:nzt) :: swap_flux_y_local, & swap_diss_y_local REAL, DIMENSION(nzb+1:nzt,nys:nyn) :: swap_flux_x_local, & swap_diss_x_local CHARACTER (LEN = *), INTENT(IN) :: sk_char IF ( boundary_flags(j,i) /= 0 ) THEN ! !-- Degrade the order for Dirichlet bc. at the outflow boundary. SELECT CASE ( boundary_flags(j,i) ) CASE ( 1 ) DO k = nzb_s_inner(j,i) + 1, nzt u_comp = u(k,j,i+1) - u_gtrans flux_r(k) = u_comp * ( & 7. * ( sk(k,j,i+1) + sk(k,j,i) ) & - ( sk(k,j,i+2) + sk(k,j,i-1) ) ) * adv_sca_3 diss_r(k) = - abs(u_comp) * ( & 3. * ( sk(k,j,i+1)-sk(k,j,i) ) & - ( sk(k,j,i+2) - sk(k,j,i-1) ) ) * adv_sca_3 ENDDO CASE ( 2 ) DO k = nzb_s_inner(j,i) + 1, nzt u_comp = u(k,j,i+1) - u_gtrans flux_r(k) = u_comp * ( sk(k,j,i+1) + sk(k,j,i) ) * 0.5 diss_r(k) = diss_2nd( sk(k,j,i+1), sk(k,j,i+1), sk(k,j,i), & sk(k,j,i-1), sk(k,j,i-2), u_comp, & 0.5, ddx ) ENDDO CASE ( 3 ) DO k = nzb_s_inner(j,i) + 1, nzt v_comp = v(k,j+1,i) - v_gtrans flux_n(k) = v_comp * ( & 7. * ( sk(k,j+1,i) + sk(k,j,i) ) & - ( sk(k,j+2,i) + sk(k,j-1,i) ) ) * adv_sca_3 diss_n(k) = - abs(v_comp) * ( & 3. * ( sk(k,j+1,i)-sk(k,j,i) ) & - (sk(k,j+2,i) - sk(k,j-1,i) ) ) * adv_sca_3 ENDDO CASE ( 4 ) DO k = nzb_s_inner(j,i) + 1, nzt v_comp = v(k,j+1,i) - v_gtrans flux_n(k) = v_comp* ( sk(k,j+1,i) + sk(k,j,i) ) * 0.5 diss_n(k) = diss_2nd( sk(k,j+1,i), sk(k,j+1,i), sk(k,j,i), & sk(k,j-1,i), sk(k,j-2,i), v_comp, & 0.5, ddy ) ENDDO CASE ( 5 ) DO k = nzb_w_inner(j,i) + 1, nzt u_comp = u(k,j,i+1) - u_gtrans flux_r(k) = u_comp * ( & 7. * ( sk(k,j,i+1) + sk(k,j,i) ) & - ( sk(k,j,i+2) + sk(k,j,i-1) ) ) * adv_sca_3 diss_r(k) = - abs(u_comp) * ( & 3. * ( sk(k,j,i+1)-sk(k,j,i) ) & - ( sk(k,j,i+2) - sk(k,j,i-1) ) ) * adv_sca_3 ENDDO CASE ( 6 ) DO k = nzb_s_inner(j,i) + 1, nzt u_comp = u(k,j,i+1) - u_gtrans flux_r(k) = u_comp * ( & 7. * ( sk(k,j,i+1) + sk(k,j,i) ) & - ( sk(k,j,i+2) + sk(k,j,i-1) ) ) * adv_sca_3 diss_r(k) = - abs(u_comp) * ( & 3. * ( sk(k,j,i+1)-sk(k,j,i) ) & - ( sk(k,j,i+2) - sk(k,j,i-1) ) ) * adv_sca_3 ! !-- Compute leftside fluxes for the left boundary of PE domain. u_comp = u(k,j,i) - u_gtrans swap_flux_x_local(k,j) = u_comp * ( & sk(k,j,i) + sk(k,j,i-1) ) * 0.5 swap_diss_x_local(k,j) = diss_2nd( sk(k,j,i+2),sk(k,j,i+1), & sk(k,j,i), sk(k,j,i-1), & sk(k,j,i-1), u_comp, & 0.5, ddx ) ENDDO degraded_l = .TRUE. CASE ( 7 ) DO k = nzb_s_inner(j,i) + 1, nzt v_comp = v(k,j+1,i)-v_gtrans flux_n(k) = v_comp * ( & 7. * ( sk(k,j+1,i) + sk(k,j,i) ) & - ( sk(k,j+2,i) + sk(k,j-1,i) ) ) * adv_sca_3 diss_n(k) = - abs(v_comp) * ( & 3. * ( sk(k,j+1,i) - sk(k,j,i) ) & - ( sk(k,j+2,i) - sk(k,j-1,i) ) ) * adv_sca_3 ENDDO CASE ( 8 ) DO k = nzb_s_inner(j,i) + 1, nzt v_comp = v(k,j+1,i) - v_gtrans flux_n(k) = v_comp * ( & 7. * ( sk(k,j+1,i) + sk(k,j,i) ) & - ( sk(k,j+2,i) + sk(k,j-1,i) ) ) * adv_sca_3 diss_n(k) = - abs(v_comp) * ( & 3. * ( sk(k,j+1,i) - sk(k,j,i) ) & - ( sk(k,j+2,i) - sk(k,j-1,i) ) ) * adv_sca_3 ! !-- Compute southside fluxes for the south boundary of PE domain v_comp = v(k,j,i) - v_gtrans swap_flux_y_local(k) = v_comp * ( & sk(k,j,i)+ sk(k,j-1,i) ) * 0.5 swap_diss_y_local(k) = diss_2nd( sk(k,j+2,i), sk(k,j+1,i), & sk(k,j,i), sk(k,j-1,i), & sk(k,j-1,i), v_comp, & 0.5, ddy ) ENDDO degraded_s = .TRUE. CASE DEFAULT END SELECT ! !-- Compute the crosswise 5th order fluxes at the outflow IF ( boundary_flags(j,i) == 1 .OR. boundary_flags(j,i) == 2 & .OR. boundary_flags(j,i) == 5 .OR. boundary_flags(j,i) == 6 ) THEN DO k = nzb_s_inner(j,i) + 1, nzt v_comp = v(k,j+1,i) - v_gtrans flux_n(k) = v_comp * ( & 37. * ( sk(k,j+1,i) + sk(k,j,i) ) & - 8. * ( 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. * ( sk(k,j+1,i) - sk(k,j,i) ) & - 5. * ( sk(k,j+2,i) - sk(k,j-1,i) ) & + ( sk(k,j+3,i) - sk(k,j-2,i) ) ) * adv_sca_5 ENDDO ELSE DO k = nzb_s_inner(j,i) + 1, nzt u_comp = u(k,j,i+1) - u_gtrans flux_r(k) = u_comp * ( & 37. * ( sk(k,j,i+1) + sk(k,j,i) ) & - 8. * ( 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. * ( sk(k,j,i+1) - sk(k,j,i) ) & - 5. * ( sk(k,j,i+2) - sk(k,j,i-1) ) & + ( sk(k,j,i+3) - sk(k,j,i-2) ) ) * adv_sca_5 ENDDO ENDIF ! !-- Compute the fifth order fluxes for the interior of PE domain. ELSE DO k = nzb_u_inner(j,i) + 1, nzt u_comp = u(k,j,i+1) - u_gtrans flux_r(k) = u_comp * ( & 37. * ( sk(k,j,i+1) + sk(k,j,i) ) & - 8. * ( 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. * ( sk(k,j,i+1) - sk(k,j,i) ) & - 5. * ( 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. * ( sk(k,j+1,i) + sk(k,j,i) ) & - 8. * ( 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. * ( sk(k,j+1,i) - sk(k,j,i) ) & - 5. * ( sk(k,j+2,i) - sk(k,j-1,i) ) & + ( sk(k,j+3,i) - sk(k,j-2,i) ) ) * adv_sca_5 ENDDO ENDIF ! !-- Compute left- and southside fluxes of the respective PE bounds. IF ( j == nys .AND. ( .NOT. degraded_s ) ) THEN DO k = nzb_s_inner(j,i) + 1, nzt v_comp = v(k,j,i) - v_gtrans swap_flux_y_local(k) = v_comp * ( & 37. * ( sk(k,j,i) + sk(k,j-1,i) ) & - 8. * ( 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. * ( sk(k,j,i) - sk(k,j-1,i) ) & - 5. * ( 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 IF ( i == nxl .AND. ( .NOT. degraded_l ) ) THEN DO k = nzb_s_inner(j,i) + 1, nzt u_comp = u(k,j,i) - u_gtrans swap_flux_x_local(k,j) = u_comp * ( & 37. * ( sk(k,j,i) + sk(k,j,i-1) ) & - 8. * ( 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. * ( sk(k,j,i) - sk(k,j,i-1) ) & - 5. * ( 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 ! !-- Now compute the tendency terms for the horizontal parts. DO k = nzb_s_inner(j,i) + 1, nzt 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 ) 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) ENDDO ! !-- Vertical advection, degradation of order near surface and top. !-- The fluxes flux_d and diss_d at the surface are 0. Due to reasons of !-- statistical evaluation the top flux at the surface should be 0 flux_t(nzb_s_inner(j,i)) = 0. diss_t(nzb_s_inner(j,i)) = 0. ! !-- 2nd order scheme k = nzb_s_inner(j,i) + 1 flux_d = flux_t(k-1) diss_d = diss_t(k-1) flux_t(k) = w(k,j,i) * ( sk(k+1,j,i) + sk(k,j,i) ) * 0.5 ! !-- sk(k,j,i) is referenced three times to avoid a access below surface. diss_t(k) = diss_2nd( sk(k+2,j,i), sk(k+1,j,i), sk(k,j,i), sk(k,j,i), & sk(k,j,i), w(k,j,i), 0.5, ddzw(k) ) tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) & - ( flux_d + diss_d ) ) * ddzw(k) ! !-- WS3 as an intermediate step k = nzb_s_inner(j,i) + 2 flux_d = flux_t(k-1) diss_d = diss_t(k-1) flux_t(k) = w(k,j,i) * ( & 7. * ( sk(k+1,j,i) + sk(k,j,i) ) & - ( sk(k+2,j,i) + sk(k-1,j,i) ) ) * adv_sca_3 diss_t(k) = - abs(w(k,j,i)) * ( & 3. * ( sk(k+1,j,i) - sk(k,j,i) ) & - ( sk(k+2,j,i) - sk(k-1,j,i) ) ) * adv_sca_3 tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) & - ( flux_d + diss_d ) ) * ddzw(k) ! !-- WS5 DO k = nzb_s_inner(j,i) + 3, nzt - 2 flux_d = flux_t(k-1) diss_d = diss_t(k-1) flux_t(k) = w(k,j,i) * ( & 37. * ( sk(k+1,j,i) + sk(k,j,i) ) & - 8. * ( sk(k+2,j,i) + sk(k-1,j,i) ) & + ( sk(k+3,j,i) + sk(k-2,j,i) ) ) *adv_sca_5 diss_t(k) = - abs(w(k,j,i)) * ( & 10. * ( sk(k+1,j,i) - sk(k,j,i) ) & - 5. * ( sk(k+2,j,i) - sk(k-1,j,i) ) & + ( sk(k+3,j,i) - sk(k-2,j,i) ) ) * adv_sca_5 tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) & - ( flux_d + diss_d ) ) * ddzw(k) ENDDO ! !-- WS3 as an intermediate step k = nzt - 1 flux_d = flux_t(k-1) diss_d = diss_t(k-1) flux_t(k) = w(k,j,i) * ( & 7. * ( sk(k+1,j,i) + sk(k,j,i) ) & - ( sk(k+2,j,i) + sk(k-1,j,i) ) ) * adv_sca_3 diss_t(k) = - abs(w(k,j,i)) * ( & 3. * ( sk(k+1,j,i) - sk(k,j,i) ) & - ( sk(k+2,j,i) - sk(k-1,j,i) ) ) * adv_sca_3 tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) & - ( flux_d + diss_d ) ) * ddzw(k) ! !-- 2nd k = nzt flux_d = flux_t(k-1) diss_d = diss_t(k-1) flux_t(k) = w(k,j,i) * ( sk(k+1,j,i) + sk(k,j,i) ) *0.5 ! !-- sk(k+1) is referenced two times to avoid a segmentation fault at top. diss_t(k) = diss_2nd( sk(k+1,j,i), sk(k+1,j,i), sk(k,j,i), sk(k-1,j,i),& sk(k-2,j,i), w(k,j,i), 0.5, ddzw(k) ) tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) & - ( flux_d + diss_d ) ) * ddzw(k) ! !-- evaluation of statistics SELECT CASE ( sk_char ) CASE ( 'pt' ) DO k = nzb_s_inner(j,i), nzt sums_wspts_ws_l(k,:) = sums_wspts_ws_l(k,:) & + (flux_t(k)+diss_t(k)) & * weight_substep(intermediate_timestep_count) * rmask(j,i,:) ENDDO CASE ( 'sa' ) DO k = nzb_s_inner(j,i), nzt sums_wssas_ws_l(k,:) = sums_wssas_ws_l(k,:) & +(flux_t(k)+diss_t(k)) & * weight_substep(intermediate_timestep_count) * rmask(j,i,:) ENDDO CASE ( 'q' ) DO k = nzb_s_inner(j,i), nzt sums_wsqs_ws_l(k,:) = sums_wsqs_ws_l(k,:) & + (flux_t(k)+diss_t(k)) & * weight_substep(intermediate_timestep_count) * rmask(j,i,:) 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 ) USE arrays_3d USE constants USE control_parameters USE grid_variables USE indices USE statistics IMPLICIT NONE INTEGER :: i, j, k LOGICAL :: degraded_l = .FALSE., degraded_s = .FALSE. REAL :: gu, gv, flux_d, diss_d, u_comp_l, v_comp, w_comp REAL, DIMENSION(nzb:nzt+1) :: flux_t, diss_t, flux_r, diss_r, & flux_n, diss_n, u_comp gu = 2.0 * u_gtrans gv = 2.0 * v_gtrans IF ( boundary_flags(j,i) /= 0 ) THEN ! !-- Degrade the order for Dirichlet bc. at the outflow boundary SELECT CASE ( boundary_flags(j,i) ) CASE ( 1 ) DO k = nzb_u_inner(j,i) + 1, nzt u_comp(k) = u(k,j,i+1) + u(k,j,i) flux_r(k) = ( u_comp(k) - gu ) * ( & 7. * (u(k,j,i+1) + u(k,j,i) ) & - ( u(k,j,i+2) + u(k,j,i-1) ) ) * adv_mom_3 diss_r(k) = - abs(u_comp(k) - gu) * ( & 3. * (u(k,j,i+1) - u(k,j,i) ) & - ( u(k,j,i+2) - u(k,j,i-1) ) ) * adv_mom_3 ENDDO CASE ( 2 ) DO k = nzb_u_inner(j,i) + 1, nzt u_comp(k) = u(k,j,i+1) + u(k,j,i) flux_r(k) = (u_comp(k) - gu) * ( u(k,j,i+1) + u(k,j,i) )*0.25 diss_r(k) = diss_2nd( u(k,j,i+1) ,u(k,j,i+1), u(k,j,i), & u(k,j,i-1), u(k,j,i-2), u_comp(k), & 0.25, ddx ) ENDDO CASE ( 3 ) DO k = nzb_u_inner(j,i) + 1, nzt v_comp = v(k,j+1,i) + v(k,j+1,i-1) - gv flux_n(k) = v_comp * ( & 7. * ( u(k,j+1,i) + u(k,j,i) ) & - ( u(k,j+2,i) + u(k,j-1,i) ) ) * adv_mom_3 diss_n(k) = - abs(v_comp) * ( & 3. * ( u(k,j+1,i) - u(k,j,i) ) & - ( u(k,j+2,i) - u(k,j-1,i) ) ) * adv_mom_3 ENDDO CASE ( 4 ) DO k = nzb_u_inner(j,i) + 1, nzt v_comp = v(k,j+1,i) + v(k,j+1,i-1) - gv flux_n(k) = v_comp * ( u(k,j+1,i) + u(k,j,i) ) * 0.25 diss_n(k) = diss_2nd( u(k,j+1,i), u(k,j+1,i), u(k,j,i), & u(k,j-1,i), u(k,j-2,i), v_comp, & 0.25, ddy ) ENDDO CASE ( 5 ) DO k = nzb_u_inner(j,i) + 1, nzt ! !-- Compute leftside fluxes for the left boundary of PE domain u_comp(k) = u(k,j,i) + u(k,j,i-1) - gu flux_l_u(k,j) = u_comp(k) * ( u(k,j,i) + u(k,j,i-1) ) * 0.25 diss_l_u(k,j) = diss_2nd(u(k,j,i+2), u(k,j,i+1), u(k,j,i), & u(k,j,i-1), u(k,j,i-1), u_comp(k), & 0.25, ddx ) u_comp(k) = u(k,j,i+1) + u(k,j,i) flux_r(k) = ( u_comp(k) - gu ) * ( & 7. * (u(k,j,i+1) + u(k,j,i) ) & - ( u(k,j,i+2) + u(k,j,i-1) ) ) * adv_mom_3 diss_r(k) = - abs( u_comp(k) -gu ) * ( & 3. * ( u(k,j,i+1) - u(k,j,i) ) & - ( u(k,j,i+2) - u(k,j,i-1) ) ) * adv_mom_3 ENDDO degraded_l = .TRUE. CASE ( 7 ) DO k = nzb_u_inner(j,i) + 1, nzt v_comp = v(k,j+1,i) + v(k,j+1,i-1) - gv flux_n(k) = v_comp * ( & 7. * ( u(k,j+1,i) + u(k,j,i) ) & - ( u(k,j+2,i) + u(k,j-1,i) ) ) * adv_mom_3 diss_n(k) = - abs(v_comp) * ( & 3. * ( u(k,j+1,i) - u(k,j,i) ) & - ( u(k,j+2,i) - u(k,j-1,i) ) ) * adv_mom_3 ENDDO CASE ( 8 ) DO k = nzb_u_inner(j,i) + 1, nzt ! !-- Compute southside fluxes for the south boundary of PE domain v_comp = v(k,j,i) + v(k,j,i-1) - gv flux_s_u(k) = v_comp * ( u(k,j,i) + u(k,j-1,i) ) * 0.25 diss_s_u(k) = diss_2nd( u(k,j+2,i), u(k,j+1,i), u(k,j,i), & u(k,j-1,i), u(k,j-1,i), v_comp, & 0.25, ddy ) v_comp = v(k,j+1,i) + v(k,j+1,i-1) - gv flux_n(k) = v_comp * ( & 7. * ( u(k,j+1,i) + u(k,j,i) ) & - ( u(k,j+2,i) + u(k,j-1,i) ) ) * adv_mom_3 diss_n(k) = - abs(v_comp) * ( & 3. * ( u(k,j+1,i) - u(k,j,i) ) & - ( u(k,j+2,i) - u(k,j-1,i) ) ) * adv_mom_3 ENDDO degraded_s = .TRUE. CASE DEFAULT END SELECT ! !-- Compute the crosswise 5th order fluxes at the outflow IF ( boundary_flags(j,i) == 1 .OR. boundary_flags(j,i) == 2 & .OR. boundary_flags(j,i) == 5 ) THEN DO k = nzb_u_inner(j,i)+1, nzt v_comp = v(k,j+1,i) + v(k,j+1,i-1) - gv flux_n(k) = v_comp * ( & 37. * ( u(k,j+1,i) + u(k,j,i) ) & - 8. * ( 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. * ( u(k,j+1,i) - u(k,j,i) ) & - 5. * ( u(k,j+2,i) - u(k,j-1,i) ) & + ( u(k,j+3,i) - u(k,j-2,i) ) ) * adv_mom_5 ENDDO ELSE DO k = nzb_u_inner(j,i) + 1, nzt u_comp(k) = u(k,j,i+1) + u(k,j,i) flux_r(k) = ( u_comp(k) - gu ) * ( & 37. * ( u(k,j,i+1) + u(k,j,i) ) & - 8. * ( 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. * ( u(k,j,i+1) - u(k,j,i) ) & - 5. * ( u(k,j,i+2) - u(k,j,i-1) ) & + ( u(k,j,i+3) - u(k,j,i-2) ) ) *adv_mom_5 ENDDO ENDIF ELSE ! !-- Compute the fifth order fluxes for the interior of PE domain. DO k = nzb_u_inner(j,i) + 1, nzt u_comp(k) = u(k,j,i+1) + u(k,j,i) flux_r(k) = ( u_comp(k) - gu ) * ( & 37. * ( u(k,j,i+1) + u(k,j,i) ) & - 8. * ( 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. * ( u(k,j,i+1) - u(k,j,i) ) & - 5. * ( 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. * ( u(k,j+1,i) + u(k,j,i) ) & - 8. * ( 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. * ( u(k,j+1,i) - u(k,j,i) ) & - 5. * ( u(k,j+2,i) - u(k,j-1,i) ) & + ( u(k,j+3,i) - u(k,j-2,i) ) ) * adv_mom_5 ENDDO ENDIF ! !-- Compute left- and southside fluxes for the respective boundary of PE IF ( j == nys .AND. ( .NOT. degraded_s ) ) THEN DO k = nzb_u_inner(j,i) + 1, nzt v_comp = v(k,j,i) + v(k,j,i-1) - gv flux_s_u(k) = v_comp * ( & 37. * ( u(k,j,i) + u(k,j-1,i) ) & - 8. * ( 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) = - abs(v_comp) * ( & 10. * ( u(k,j,i) - u(k,j-1,i) ) & - 5. * ( 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 IF ( i == nxlu .AND. ( .NOT. degraded_l ) ) THEN DO k = nzb_u_inner(j,i)+1, nzt u_comp_l = u(k,j,i)+u(k,j,i-1)-gu flux_l_u(k,j) = u_comp_l * ( & 37. * ( u(k,j,i) + u(k,j,i-1) ) & - 8. * ( 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) = - abs(u_comp_l) * ( & 10. * ( u(k,j,i) - u(k,j,i-1) ) & - 5. * ( 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 ! !-- Now compute the tendency terms for the horizontal parts. DO k = nzb_u_inner(j,i) + 1, nzt tend(k,j,i) = tend(k,j,i) - ( & ( flux_r(k) + diss_r(k) & - ( flux_l_u(k,j) + diss_l_u(k,j) ) ) * ddx & + ( flux_n(k) + diss_n(k) & - ( flux_s_u(k) + diss_s_u(k) ) ) * ddy ) flux_l_u(k,j) = flux_r(k) diss_l_u(k,j) = diss_r(k) flux_s_u(k) = flux_n(k) diss_s_u(k) = diss_n(k) ! !-- Statistical Evaluation of u'u'. The factor has to be applied for !-- right evaluation when gallilei_trans = .T. . sums_us2_ws_l(k,:) = sums_us2_ws_l(k,:) & + ( flux_r(k) * & ( u_comp(k) - 2. * hom(k,1,1,:) ) / ( u_comp(k) - gu + 1.0E-20 ) & + diss_r(k) & * abs(u_comp(k) - 2. * hom(k,1,1,:) ) & / (abs(u_comp(k) - gu) + 1.0E-20) ) & * weight_substep(intermediate_timestep_count) * rmask(j,i,:) ENDDO sums_us2_ws_l(nzb_u_inner(j,i),:) = & sums_us2_ws_l(nzb_u_inner(j,i)+1,:) ! !-- Vertical advection, degradation of order near surface and top. !-- The fluxes flux_d and diss_d at the surface are 0. Due to reasons of !-- statistical evaluation the top flux at the surface should be 0 flux_t(nzb_u_inner(j,i)) = 0. !statistical reasons diss_t(nzb_u_inner(j,i)) = 0. ! !-- 2nd order scheme k = nzb_u_inner(j,i) + 1 flux_d = flux_t(k-1) diss_d = diss_t(k-1) w_comp = w(k,j,i) + w(k,j,i-1) flux_t(k) = w_comp * ( u(k+1,j,i) + u(k,j,i) ) *0.25 diss_t(k) = diss_2nd( u(k+2,j,i), u(k+1,j,i), u(k,j,i), 0., 0., & w_comp, 0.25, ddzw(k) ) tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) & - ( flux_d + diss_d ) ) * ddzw(k) ! !-- WS3 as an intermediate step k = nzb_u_inner(j,i) + 2 flux_d = flux_t(k-1) diss_d = diss_t(k-1) w_comp = w(k,j,i) + w(k,j,i-1) flux_t(k) = w_comp*( & 7. * ( u(k+1,j,i) + u(k,j,i) ) & - ( u(k+2,j,i) + u(k-1,j,i) ) ) * adv_mom_3 diss_t(k) = -abs(w_comp)*( & 3. * ( u(k+1,j,i) - u(k,j,i) ) & - ( u(k+2,j,i) - u(k-1,j,i) ) ) * adv_mom_3 tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) & - ( flux_d + diss_d ) ) * ddzw(k) DO k = nzb_u_inner(j,i) + 3, nzt - 2 flux_d = flux_t(k-1) diss_d = diss_t(k-1) w_comp = w(k,j,i) + w(k,j,i-1) flux_t(k) = w_comp*( & 37.* ( u(k+1,j,i) + u(k,j,i) ) & - 8. * ( u(k+2,j,i) + u(k-1,j,i) ) & + ( u(k+3,j,i) + u(k-2,j,i) ) ) * adv_mom_5 diss_t(k) = - abs(w_comp) * ( & 10. * ( u(k+1,j,i) - u(k,j,i) ) & - 5. * ( u(k+2,j,i) - u(k-1,j,i) ) & + ( u(k+3,j,i) - u(k-2,j,i) ) ) * adv_mom_5 tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) & - ( flux_d + diss_d ) ) * ddzw(k) ENDDO ! !-- WS3 as an intermediate step k = nzt - 1 flux_d = flux_t(k-1) diss_d = diss_t(k-1) w_comp = w(k,j,i) + w(k,j,i-1) flux_t(k) = w_comp * ( & 7. * ( u(k+1,j,i) + u(k,j,i) ) & - ( u(k+2,j,i) + u(k-1,j,i) ) ) * adv_mom_3 diss_t(k) = - abs(w_comp) * ( & 3. * ( u(k+1,j,i) - u(k,j,i) ) & - ( u(k+2,j,i) - u(k-1,j,i) ) ) * adv_mom_3 tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) & - ( flux_d + diss_d ) ) * ddzw(k) ! !-- 2nd order scheme k = nzt flux_d = flux_t(k-1) diss_d = diss_t(k-1) w_comp = w(k,j,i)+w(k,j,i-1) flux_t(k) = w_comp * ( u(k+1,j,i) + u(k,j,i) ) * 0.25 diss_t(k) = diss_2nd( u(k+1,j,i), u(k+1,j,i), u(k,j,i), u(k-1,j,i), & u(k-2,j,i), w_comp, 0.25, ddzw(k) ) tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) & - ( flux_d + diss_d ) ) * ddzw(k) ! !-- sum up the vertical momentum fluxes DO k = nzb_u_inner(j,i), nzt sums_wsus_ws_l(k,:) = sums_wsus_ws_l(k,:) & + ( flux_t(k) + diss_t(k) ) & * weight_substep(intermediate_timestep_count) * rmask(j,i,:) ENDDO END SUBROUTINE advec_u_ws_ij !------------------------------------------------------------------------------! ! Advection of v-component - Call for grid point i,j !------------------------------------------------------------------------------! SUBROUTINE advec_v_ws_ij( i, j ) USE arrays_3d USE constants USE control_parameters USE grid_variables USE indices USE statistics IMPLICIT NONE INTEGER :: i, j, k LOGICAL :: degraded_l = .FALSE., degraded_s = .FALSE. REAL :: gu, gv, flux_d, diss_d, u_comp, v_comp_l, w_comp REAL, DIMENSION(nzb:nzt+1) :: flux_t, diss_t, flux_n, & diss_n, flux_r, diss_r, v_comp gu = 2.0 * u_gtrans gv = 2.0 * v_gtrans IF ( boundary_flags(j,i) /= 0 ) THEN ! !-- Degrade the order for Dirichlet bc. at the outflow boundary SELECT CASE ( boundary_flags(j,i) ) CASE ( 1 ) DO k = nzb_v_inner(j,i) + 1, nzt u_comp = u(k,j-1,i+1) + u(k,j,i+1) - gu flux_r(k) = u_comp * ( & 7. * (v(k,j,i+1) + v(k,j,i) ) & - ( v(k,j,i+2) + v(k,j,i-1) ) ) * adv_mom_3 diss_r(k) = - abs(u_comp) * ( & 3. * ( v(k,j,i+1) - v(k,j,i) ) & - ( v(k,j,i+2) - v(k,j,i-1) ) ) * adv_mom_3 ENDDO CASE ( 2 ) DO k = nzb_v_inner(j,i) + 1, nzt u_comp = u(k,j-1,i+1) + u(k,j,i+1) - gu flux_r(k) = u_comp * ( v(k,j,i+1) + v(k,j,i) ) * 0.25 diss_r(k) = diss_2nd( v(k,j,i+1), v(k,j,i+1), v(k,j,i), & v(k,j,i-1), v(k,j,i-2), u_comp, & 0.25, ddx ) ENDDO CASE ( 3 ) DO k = nzb_v_inner(j,i) + 1, nzt v_comp(k) = v(k,j+1,i) + v(k,j,i) flux_n(k) = ( v_comp(k)- gv ) * ( & 7. * ( v(k,j+1,i) + v(k,j,i) ) & - ( v(k,j+2,i) + v(k,j-1,i) ) ) * adv_mom_3 diss_n(k) = - abs(v_comp(k) - gv) * ( & 3. * ( v(k,j+1,i) - v(k,j,i) ) & - ( v(k,j+2,i) - v(k,j-1,i) ) ) * adv_mom_3 ENDDO CASE ( 4 ) DO k = nzb_v_inner(j,i) + 1, nzt v_comp(k) = v(k,j+1,i) + v(k,j,i) flux_n(k) = ( v_comp(k) - gv ) * & ( v(k,j+1,i) + v(k,j,i) ) * 0.25 diss_n(k) = diss_2nd( v(k,j+1,i), v(k,j+1,i), v(k,j,i), & v(k,j-1,i), v(k,j-2,i), v_comp(k), & 0.25, ddy ) ENDDO CASE ( 5 ) DO k = nzb_w_inner(j,i) + 1, nzt u_comp = u(k,j-1,i) + u(k,j,i) - gu flux_r(k) = u_comp * ( & 7. * ( v(k,j,i+1) + v(k,j,i) ) & - ( v(k,j,i+2) + v(k,j,i-1) ) ) * adv_mom_3 diss_r(k) = - abs(u_comp) * ( & 3. * ( w(k,j,i+1) - w(k,j,i) ) & - ( v(k,j,i+2) - v(k,j,i-1) ) ) * adv_mom_3 ENDDO CASE ( 6 ) DO k = nzb_v_inner(j,i) + 1, nzt u_comp = u(k,j-1,i) + u(k,j,i) - gu flux_l_v(k,j) = u_comp * ( v(k,j,i) + v(k,j,i-1) ) * 0.25 diss_l_v(k,j) = diss_2nd( v(k,j,i+2), v(k,j,i+1), v(k,j,i),& v(k,j,i-1), v(k,j,i-1), u_comp, & 0.25, ddx ) u_comp = u(k,j-1,i+1) + u(k,j,i+1) - gu flux_r(k) = u_comp * ( & 7. * ( v(k,j,i+1) + v(k,j,i) ) & - ( v(k,j,i+2) + v(k,j,i-1) ) ) * adv_mom_3 diss_r(k) = - abs(u_comp) * ( & 3. * ( v(k,j,i+1) - v(k,j,i) ) & - ( v(k,j,i+2) - v(k,j,i-1) ) ) * adv_mom_3 ENDDO degraded_l = .TRUE. CASE ( 7 ) DO k = nzb_v_inner(j,i) + 1, nzt v_comp(k) = v(k,j,i) + v(k,j-1,i) - gv flux_s_v(k) = v_comp(k) * ( v(k,j,i) + v(k,j-1,i) ) * 0.25 diss_s_v(k) = diss_2nd( v(k,j+2,i), v(k,j+1,i), v(k,j,i), & v(k,j-1,i), v(k,j-1,i), v_comp(k), & 0.25, ddy ) v_comp(k) = v(k,j+1,i) + v(k,j,i) flux_n(k) = ( v_comp(k) - gv ) * ( & 7. * ( v(k,j+1,i) + v(k,j,i) ) & - ( v(k,j+2,i) + v(k,j-1,i) ) ) * adv_mom_3 diss_n(k) = - abs(v_comp(k) - gv) * ( & 3. * ( v(k,j+1,i) - v(k,j,i) ) & - ( v(k,j+2,i) - v(k,j-1,i) ) ) * adv_mom_3 ENDDO degraded_s = .TRUE. CASE DEFAULT END SELECT ! !-- Compute the crosswise 5th order fluxes at the outflow IF ( boundary_flags(j,i) == 1 .OR. boundary_flags(j,i) == 2 & .OR. boundary_flags(j,i) == 5 .OR. boundary_flags(j,i) == 6 ) THEN DO k = nzb_v_inner(j,i) + 1, nzt v_comp(k) = v(k,j+1,i) + v(k,j,i) flux_n(k) = ( v_comp(k) - gv ) * ( & 37. * ( v(k,j+1,i) + v(k,j,i) ) & - 8. * ( 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. * ( v(k,j+1,i) - v(k,j,i) ) & - 5. * ( v(k,j+2,i) - v(k,j-1,i) ) & + ( v(k,j+3,i) - v(k,j-2,i) ) ) *adv_mom_5 ENDDO ELSE DO k = nzb_v_inner(j,i) + 1, nzt u_comp = u(k,j-1,i+1) + u(k,j,i+1) - gu flux_r(k) = u_comp * ( & 37. * ( v(k,j,i+1) + v(k,j,i) ) & - 8. * ( 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. * ( v(k,j,i+1) - v(k,j,i) ) & -5. * ( v(k,j,i+2) - v(k,j,i-1) ) & + ( v(k,j,i+3) - v(k,j,i-2) ) ) * adv_mom_5 ENDDO ENDIF ELSE ! !-- Compute the fifth order fluxes for the interior of PE domain. DO k = nzb_v_inner(j,i) + 1, nzt u_comp = u(k,j-1,i+1) + u(k,j,i+1) - gu flux_r(k) = u_comp * ( & 37. * ( v(k,j,i+1) + v(k,j,i) ) & - 8. * ( 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. * ( v(k,j,i+1) - v(k,j,i) ) & - 5. * ( 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. * ( v(k,j+1,i) + v(k,j,i) ) & - 8. * ( 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. * ( v(k,j+1,i) - v(k,j,i) ) & - 5. * ( v(k,j+2,i) - v(k,j-1,i) ) & + ( v(k,j+3,i) - v(k,j-2,i) ) ) * adv_mom_5 ENDDO ENDIF ! !-- Compute left- and southside fluxes for the respective boundary IF ( i == nxl .AND. ( .NOT. degraded_l ) ) THEN DO k = nzb_v_inner(j,i) + 1, nzt u_comp = u(k,j-1,i) + u(k,j,i) - gu flux_l_v(k,j) = u_comp * ( & 37. * ( v(k,j,i) + v(k,j,i-1) ) & - 8. * ( 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) = - abs(u_comp) * ( & 10. * ( v(k,j,i) - v(k,j,i-1) ) & - 5. * ( 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 IF ( j == nysv .AND. ( .NOT. degraded_s ) ) THEN DO k = nzb_v_inner(j,i) + 1, nzt v_comp_l = v(k,j,i) + v(k,j-1,i) - gv flux_s_v(k) = v_comp_l * ( & 37. * ( v(k,j,i) + v(k,j-1,i) ) & - 8. * ( 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) = - abs(v_comp_l) * ( & 10. * ( v(k,j,i) - v(k,j-1,i) ) & - 5. * ( 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 ! !-- Now compute the tendency terms for the horizontal parts. DO k = nzb_v_inner(j,i) + 1, nzt tend(k,j,i) = tend(k,j,i) - ( & ( flux_r(k) + diss_r(k) & - ( flux_l_v(k,j) + diss_l_v(k,j) ) ) * ddx & + ( flux_n(k) + diss_n(k) & - ( flux_s_v(k) + diss_s_v(k) ) ) * ddy ) flux_l_v(k,j) = flux_r(k) diss_l_v(k,j) = diss_r(k) flux_s_v(k) = flux_n(k) diss_s_v(k) = diss_n(k) ! !-- Statistical Evaluation of v'v'. The factor has to be applied for !-- right evaluation when gallilei_trans = .T. . sums_vs2_ws_l(k,:) = sums_vs2_ws_l(k,:) & + ( flux_n(k) & * ( v_comp(k) - 2. * hom(k,1,2,:) ) & / ( v_comp(k) - gv + 1.0E-20 ) & + diss_n(k) & * abs( v_comp(k) - 2. * hom(k,1,2,:) ) & / ( abs( v_comp(k) - gv ) +1.0E-20 ) ) & * weight_substep(intermediate_timestep_count) * rmask(j,i,:) ENDDO sums_vs2_ws_l(nzb_v_inner(j,i),:) = & sums_vs2_ws_l(nzb_v_inner(j,i)+1,:) ! !-- Vertical advection, degradation of order near surface and top. !-- The fluxes flux_d and diss_d at the surface are 0. Due to reasons of !-- statistical evaluation the top flux at the surface should be 0 flux_t(nzb_v_inner(j,i)) = 0. !statistical reasons diss_t(nzb_v_inner(j,i)) = 0. ! !-- 2nd order scheme k = nzb_v_inner(j,i) + 1 flux_d = flux_t(k-1) diss_d = diss_t(k-1) w_comp = w(k,j-1,i) + w(k,j,i) flux_t(k) = w_comp * ( v(k+1,j,i) + v(k,j,i) ) * 0.25 diss_t(k) = diss_2nd( v(k+2,j,i), v(k+1,j,i), v(k,j,i), 0., 0., w_comp,& 0.25, ddzw(k) ) tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) & - ( flux_d + diss_d ) ) * ddzw(k) ! !-- WS3 as an intermediate step k = nzb_v_inner(j,i) + 2 flux_d = flux_t(k-1) diss_d = diss_t(k-1) w_comp = w(k,j-1,i) + w(k,j,i) flux_t(k) = w_comp * ( & 7. * ( v(k+1,j,i) + v(k,j,i) ) & - ( v(k+2,j,i) + v(k-1,j,i) ) ) * adv_mom_3 diss_t(k) = - abs(w_comp) * ( & 3. * ( v(k+1,j,i) - v(k,j,i) ) & - ( v(k+2,j,i) - v(k-1,j,i) ) ) * adv_mom_3 tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) & - ( flux_d + diss_d ) ) * ddzw(k) !-- WS5 DO k = nzb_v_inner(j,i) + 3, nzt - 2 flux_d = flux_t(k-1) diss_d = diss_t(k-1) w_comp = w(k,j-1,i) + w(k,j,i) flux_t(k) = w_comp * ( & 37. * ( v(k+1,j,i) + v(k,j,i) ) & - 8. * ( v(k+2,j,i) + v(k-1,j,i) ) & + ( v(k+3,j,i) + v(k-2,j,i) ) ) * adv_mom_5 diss_t(k) = - abs(w_comp) * ( & 10. * ( v(k+1,j,i) - v(k,j,i) ) & - 5. * ( v(k+2,j,i) - v(k-1,j,i) ) & + ( v(k+3,j,i) - v(k-2,j,i) ) ) * adv_mom_5 tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) & - ( flux_d + diss_d ) ) * ddzw(k) ENDDO ! !-- WS3 as an intermediate step k = nzt - 1 flux_d = flux_t(k-1) diss_d = diss_t(k-1) w_comp = w(k,j-1,i) + w(k,j,i) flux_t(k) = w_comp * ( & 7. * ( v(k+1,j,i) + v(k,j,i) ) & - ( v(k+2,j,i) + v(k-1,j,i) ) ) * adv_mom_3 diss_t(k) = - abs(w_comp) * ( & 3. * ( v(k+1,j,i) - v(k,j,i) ) & - ( v(k+2,j,i) - v(k-1,j,i) ) ) * adv_mom_3 tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) & - ( flux_d + diss_d ) ) * ddzw(k) ! !-- 2nd order scheme k = nzt flux_d = flux_t(k-1) diss_d = diss_t(k-1) w_comp = w(k,j-1,i)+w(k,j,i) flux_t(k) = w_comp * ( v(k+1,j,i) + v(k,j,i) ) * 0.25 diss_t(k) = diss_2nd( v(k+1,j,i), v(k+1,j,i), v(k,j,i), v(k-1,j,i), & v(k-2,j,i), w_comp, 0.25, ddzw(k) ) tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) & - ( flux_d + diss_d ) ) * ddzw(k) DO k = nzb_v_inner(j,i), nzt sums_wsvs_ws_l(k,:) = sums_wsvs_ws_l(k,:) & + ( flux_t(k) + diss_t(k) ) & * weight_substep(intermediate_timestep_count) * rmask(j,i,:) ENDDO END SUBROUTINE advec_v_ws_ij !------------------------------------------------------------------------------! ! Advection of w-component - Call for grid point i,j !------------------------------------------------------------------------------! SUBROUTINE advec_w_ws_ij( i, j ) USE arrays_3d USE constants USE control_parameters USE grid_variables USE indices USE statistics IMPLICIT NONE INTEGER :: i, j, k LOGICAL :: degraded_l = .FALSE., degraded_s = .FALSE. REAL :: gu, gv, flux_d, diss_d, u_comp, v_comp, w_comp REAL, DIMENSION(nzb:nzt+1) :: flux_t, diss_t, flux_r, diss_r, flux_n, & diss_n gu = 2.0 * u_gtrans gv = 2.0 * v_gtrans IF ( boundary_flags(j,i) /= 0 ) THEN ! !-- Degrade the order for Dirichlet bc. at the outflow boundary SELECT CASE ( boundary_flags(j,i) ) CASE ( 1 ) DO k = nzb_w_inner(j,i) + 1, nzt u_comp = u(k+1,j,i+1) + u(k,j,i+1) - gu flux_r(k) = u_comp * ( & 7. * ( w(k,j,i+1) + w(k,j,i) ) & - ( w(k,j,i+2) + w(k,j,i-1) ) ) * adv_mom_3 diss_r(k) = -abs(u_comp) * ( & 3. * ( w(k,j,i+1) - w(k,j,i) ) & - ( w(k,j,i+2) - w(k,j,i-1) ) ) * adv_mom_3 ENDDO CASE ( 2 ) DO k = nzb_w_inner(j,i) + 1, nzt u_comp = u(k+1,j,i+1) + u(k,j,i+1) - gu flux_r(k) = u_comp * ( w(k,j,i+1) + w(k,j,i) ) * 0.25 diss_r(k) = diss_2nd( w(k,j,i+1), w(k,j,i+1), w(k,j,i), & w(k,j,i-1), w(k,j,i-2), u_comp, & 0.25, ddx ) ENDDO CASE ( 3 ) DO k = nzb_w_inner(j,i) + 1, nzt v_comp = v(k+1,j+1,i) + v(k,j+1,i) - gv flux_n(k) = v_comp * ( & 7. * ( w(k,j+1,i) + w(k,j,i) ) & - ( w(k,j+2,i) + w(k,j-1,i) ) ) * adv_mom_3 diss_n(k) = -abs(v_comp) * ( & 3. * ( w(k,j+1,i) - w(k,j,i) ) & - ( w(k,j+2,i) - w(k,j-1,i) ) ) * adv_mom_3 ENDDO CASE ( 4 ) DO k = nzb_w_inner(j,i) + 1, nzt v_comp = v(k+1,j+1,i) + v(k,j+1,i) - gv flux_n(k) = v_comp * ( w(k,j+1,i) + w(k,j,i) ) * 0.25 diss_n(k) = diss_2nd( w(k,j+1,i), w(k,j+1,i), w(k,j,i), & w(k,j-1,i), w(k,j-2,i), v_comp, & 0.25, ddy ) ENDDO CASE ( 5 ) DO k = nzb_w_inner(j,i) + 1, nzt u_comp = u(k+1,j,i+1) + u(k,j,i+1) - gu flux_r(k) = u_comp * ( & 7. * ( w(k,j,i+1) + w(k,j,i) ) & - ( w(k,j,i+2) + w(k,j,i-1) ) ) * adv_mom_3 diss_r(k) = - abs(u_comp) * ( & 3. * ( w(k,j,i+1) - w(k,j,i) ) & - ( w(k,j,i+2) - w(k,j,i-1) ) ) * adv_mom_3 ENDDO CASE ( 6 ) DO k = nzb_w_inner(j,i) + 1, nzt ! !-- Compute leftside fluxes for the left boundary of PE domain u_comp = u(k+1,j,i+1) + u(k,j,i+1) - gu flux_r(k) = u_comp *( & 7. * ( w(k,j,i+1) + w(k,j,i) ) & - ( w(k,j,i+2) + w(k,j,i-1) ) ) * adv_mom_3 diss_r(k) = - abs(u_comp) * ( & 3. * ( w(k,j,i+1) - w(k,j,i) ) & - ( w(k,j,i+2) - w(k,j,i-1) ) ) * adv_mom_3 u_comp = u(k+1,j,i) + u(k,j,i) - gu flux_l_w(k,j) = u_comp * ( w(k,j,i) + w(k,j,i-1) ) * 0.25 diss_l_w(k,j) = diss_2nd( w(k,j,i+2), w(k,j,i+1), w(k,j,i), & w(k,j,i-1), w(k,j,i-1), u_comp, & 0.25,ddx) ENDDO degraded_l = .TRUE. CASE ( 7 ) DO k = nzb_w_inner(j,i) + 1, nzt v_comp = v(k+1,j+1,i) + v(k,j+1,i) - gv flux_n(k) = v_comp *( & 7. * ( w(k,j+1,i) + w(k,j,i) ) & - ( w(k,j+2,i) + w(k,j-1,i) ) ) * adv_mom_3 diss_n(k) = - abs(v_comp) * ( & 3. * ( w(k,j+1,i) - w(k,j,i) ) & - ( w(k,j+2,i) - w(k,j-1,i) ) ) * adv_mom_3 ENDDO CASE ( 8 ) DO k = nzb_w_inner(j,i) + 1, nzt v_comp = v(k+1,j+1,i) + v(k,j+1,i) - gv flux_n(k) = v_comp * ( & 7. * ( w(k,j+1,i) + w(k,j,i) ) & - ( w(k,j+2,i) + w(k,j-1,i) ) ) * adv_mom_3 diss_n(k) = - abs(v_comp) * ( & 3. * ( w(k,j+1,i) - w(k,j,i) ) & - ( w(k,j+2,i) - w(k,j-1,i) ) ) * adv_mom_3 ! !-- Compute southside fluxes for the south boundary of PE domain v_comp = v(k+1,j,i) + v(k,j,i) - gv flux_s_w(k) = v_comp * ( w(k,j,i) + w(k,j-1,i) ) * 0.25 diss_s_w(k) = diss_2nd( w(k,j+2,i), w(k,j+1,i), w(k,j,i), & w(k,j-1,i), w(k,j-1,i), v_comp, & 0.25, ddy ) ENDDO degraded_s = .TRUE. CASE DEFAULT END SELECT ! !-- Compute the crosswise 5th order fluxes at the outflow IF ( boundary_flags(j,i) == 1 .OR. boundary_flags(j,i) == 2 & .OR. boundary_flags(j,i) == 5 .OR. boundary_flags(j,i) == 6 ) THEN DO k = nzb_w_inner(j,i) + 1, nzt v_comp = v(k+1,j+1,i) + v(k,j+1,i) - gv flux_n(k) = v_comp * ( & 37. * ( w(k,j+1,i) + w(k,j,i) ) & - 8. * ( 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. * ( w(k,j+1,i) - w(k,j,i) ) & - 5. * ( w(k,j+2,i) - w(k,j-1,i) ) & + ( w(k,j+3,i) - w(k,j-2,i) ) ) * adv_mom_5 ENDDO ELSE DO k = nzb_w_inner(j,i) + 1, nzt u_comp = u(k+1,j,i+1) + u(k,j,i+1) - gu flux_r(k) = u_comp * ( & 37. * ( w(k,j,i+1) + w(k,j,i) ) & - 8. * ( 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. * ( w(k,j,i+1) - w(k,j,i) ) & - 5. * ( w(k,j,i+2) - w(k,j,i-1) ) & + ( w(k,j,i+3) - w(k,j,i-2) ) ) * adv_mom_5 ENDDO ENDIF ELSE ! !-- Compute the fifth order fluxes for the interior of PE domain. DO k = nzb_w_inner(j,i) + 1, nzt u_comp = u(k+1,j,i+1) + u(k,j,i+1) - gu flux_r(k) = u_comp * ( & 37. * ( w(k,j,i+1) + w(k,j,i) ) & - 8. * ( 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. * ( w(k,j,i+1) - w(k,j,i) ) & - 5. * ( 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. * ( w(k,j+1,i) + w(k,j,i) ) & - 8. * ( 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. * ( w(k,j+1,i) - w(k,j,i) ) & - 5. * ( w(k,j+2,i) - w(k,j-1,i) ) & + ( w(k,j+3,i) - w(k,j-2,i) ) ) * adv_mom_5 ENDDO ENDIF ! !-- Compute left- and southside fluxes for the respective boundary IF ( j == nys .AND. ( .NOT. degraded_s ) ) THEN DO k = nzb_w_inner(j,i) + 1, nzt v_comp = v(k+1,j,i) + v(k,j,i) - gv flux_s_w(k) = v_comp * ( & 37. * ( w(k,j,i) + w(k,j-1,i) ) & - 8. * ( 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) = - abs(v_comp) * ( & 10. * ( w(k,j,i) - w(k,j-1,i) ) & - 5. * ( 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 IF ( i == nxl .AND. ( .NOT. degraded_l ) ) THEN DO k = nzb_w_inner(j,i) + 1, nzt u_comp = u(k+1,j,i) + u(k,j,i) - gu flux_l_w(k,j) = u_comp * ( & 37. * ( w(k,j,i) + w(k,j,i-1) ) & - 8. * ( 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) = - abs(u_comp) * ( & 10. * ( w(k,j,i) - w(k,j,i-1) ) & - 5. * ( 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 ! !-- Now compute the tendency terms for the horizontal parts. DO k = nzb_w_inner(j,i) + 1, nzt tend(k,j,i) = tend(k,j,i) - ( & ( flux_r(k) + diss_r(k) & - ( flux_l_w(k,j) + diss_l_w(k,j) ) ) * ddx & + ( flux_n(k) + diss_n(k) & - ( flux_s_w(k) + diss_s_w(k) ) ) * ddy ) flux_l_w(k,j) = flux_r(k) diss_l_w(k,j) = diss_r(k) flux_s_w(k) = flux_n(k) diss_s_w(k) = diss_n(k) ENDDO ! !-- Vertical advection, degradation of order near surface and top. !-- The fluxes flux_d and diss_d at the surface are 0. Due to reasons of !-- statistical evaluation the top flux at the surface should be 0 flux_t(nzb_w_inner(j,i)) = 0. !statistical reasons diss_t(nzb_w_inner(j,i)) = 0. ! !-- 2nd order scheme k = nzb_w_inner(j,i) + 1 flux_d = flux_t(k-1) diss_d = diss_t(k-1) w_comp = w(k+1,j,i) + w(k,j,i) flux_t(k) = w_comp * ( w(k+1,j,i) + w(k,j,i) ) * 0.25 diss_t(k) = diss_2nd( w(k+2,j,i), w(k+1,j,i), w(k,j,i), 0., 0., & w_comp, 0.25, ddzu(k+1) ) tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) & - ( flux_d + diss_d ) ) * ddzu(k+1) ! !-- WS3 as an intermediate step k = nzb_w_inner(j,i) + 2 flux_d = flux_t(k-1) diss_d = diss_t(k-1) w_comp = w(k+1,j,i) + w(k,j,i) flux_t(k) = w_comp * ( & 7. * ( w(k+1,j,i) + w(k,j,i) ) & - ( w(k+2,j,i) + w(k-1,j,i) ) ) * adv_mom_3 diss_t(k) = - abs(w_comp) * ( & 3. * ( w(k+1,j,i) - w(k,j,i) ) & - ( w(k+2,j,i) - w(k-1,j,i) ) ) * adv_mom_3 tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) & - ( flux_d + diss_d ) ) * ddzu(k+1) ! !-- WS5 DO k = nzb_w_inner(j,i) + 3, nzt -2 flux_d = flux_t(k-1) diss_d = diss_t(k-1) w_comp = w(k+1,j,i) + w(k,j,i) flux_t(k) = w_comp * ( & 37. * ( w(k+1,j,i) + w(k,j,i) ) & - 8. * ( w(k+2,j,i) + w(k-1,j,i) ) & + ( w(k+3,j,i) + w(k-2,j,i) ) ) * adv_mom_5 diss_t(k) = - abs(w_comp) * ( & 10. * ( w(k+1,j,i) - w(k,j,i) ) & - 5. * ( w(k+2,j,i) - w(k-1,j,i) ) & + ( w(k+3,j,i) - w(k-2,j,i) ) ) * adv_mom_5 tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) & - ( flux_d + diss_d ) ) * ddzu(k+1) ENDDO !-- WS3 as an intermediate step k = nzt - 1 flux_d = flux_t(k-1) diss_d = diss_t(k-1) w_comp = w(k+1,j,i) + w(k,j,i) flux_t(k) = w_comp * ( & 7. * ( w(k+1,j,i) + w(k,j,i) ) & - ( w(k+2,j,i) + w(k-1,j,i) ) ) *adv_mom_3 diss_t(k) = - abs(w_comp) * ( & 3. * ( w(k+1,j,i) - w(k,j,i) ) & - ( w(k+2,j,i) - w(k-1,j,i) ) ) * adv_mom_3 tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) & - ( flux_d + diss_d ) ) * ddzu(k+1) ! !-- 2nd order scheme k = nzt flux_d = flux_t(k-1) diss_d = diss_t(k-1) w_comp = w(k+1,j,i) + w(k,j,i) flux_t(k) = w_comp * ( w(k+1,j,i) + w(k,j,i) ) * 0.25 diss_t(k) = diss_2nd( w(k+1,j,i), w(k+1,j,i), w(k,j,i), w(k-1,j,i), & w(k-2,j,i), w_comp, 0.25, ddzu(k+1) ) tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) & - ( flux_d + diss_d ) ) * ddzu(k+1) DO k = nzb_w_inner(j,i), nzt sums_ws2_ws_l(k,:) = sums_ws2_ws_l(k,:) & + ( flux_t(k) + diss_t(k) ) & * weight_substep(intermediate_timestep_count) * rmask(j,i,:) 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, j, k REAL, DIMENSION(:,:,:), POINTER :: sk REAL :: flux_d, diss_d, u_comp, v_comp REAL, DIMENSION(nzb:nzt+1) :: flux_r, diss_r, flux_n, diss_n REAL, DIMENSION(nzb+1:nzt) :: swap_flux_y_local, swap_diss_y_local, & flux_t, diss_t REAL, DIMENSION(nzb+1:nzt,nys:nyn) :: swap_flux_x_local, & swap_diss_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 IF ( boundary_flags(j,i) == 6 ) THEN DO k = nzb_s_inner(j,i) + 1, nzt u_comp = u(k,j,i) - u_gtrans swap_flux_x_local(k,j) = u_comp * ( & sk(k,j,i) + sk(k,j,i-1)) * 0.5 swap_diss_x_local(k,j) = diss_2nd( sk(k,j,i+2), sk(k,j,i+1), & sk(k,j,i), sk(k,j,i-1), & sk(k,j,i-1), u_comp, & 0.5, ddx ) ENDDO ELSE DO k = nzb_s_inner(j,i) + 1, nzt u_comp = u(k,j,i) - u_gtrans swap_flux_x_local(k,j) = u_comp*( & 37. * (sk(k,j,i)+sk(k,j,i-1) ) & -8. * ( 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. * (sk(k,j,i) - sk(k,j,i-1) ) & -5. * ( 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 ENDDO ! !-- The following loop computes the horizontal fluxes for the interior of the !-- processor domain plus south boundary points. Furthermore tendency terms !-- are computed. DO i = nxl, nxr j = nys IF ( boundary_flags(j,i) == 8 ) THEN DO k = nzb_s_inner(j,i) + 1, nzt v_comp = v(k,j,i) - v_gtrans swap_flux_y_local(k) = v_comp * & ( sk(k,j,i) + sk(k,j-1,i) ) * 0.5 swap_diss_y_local(k) = diss_2nd( sk(k,j+2,i), sk(k,j+1,i), & sk(k,j,i), sk(k,j-1,i), & sk(k,j-1,i), v_comp, 0.5, ddy) ENDDO ELSE DO k = nzb_s_inner(j,i) + 1, nzt v_comp = v(k,j,i) - v_gtrans swap_flux_y_local(k) = v_comp * ( & 37. * ( sk(k,j,i) + sk(k,j-1,i) ) & - 8. * ( 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. * ( sk(k,j,i) - sk(k,j-1,i) ) & - 5. * ( 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 DO j = nys, nyn IF ( boundary_flags(j,i) /= 0 ) THEN ! !-- Degrade the order for Dirichlet bc. at the outflow boundary SELECT CASE ( boundary_flags(j,i) ) CASE ( 1 ) DO k = nzb_s_inner(j,i) + 1, nzt u_comp = u(k,j,i+1) - u_gtrans flux_r(k) = u_comp * ( & 7. * ( sk(k,j,i+1) + sk(k,j,i) ) & - ( sk(k,j,i+2) + sk(k,j,i-1) ) ) * adv_sca_3 diss_r(k) = - abs(u_comp) * ( & 3. * ( sk(k,j,i+1) - sk(k,j,i) ) & - ( sk(k,j,i+2) - sk(k,j,i-1) ) ) * adv_sca_3 ENDDO CASE ( 2 ) DO k = nzb_s_inner(j,i) + 1, nzt u_comp = u(k,j,i+1) - u_gtrans flux_r(k) = u_comp * ( sk(k,j,i+1) + sk(k,j,i) ) * 0.5 diss_r(k) = diss_2nd( sk(k,j,i+1), sk(k,j,i+1), & sk(k,j,i), sk(k,j,i-1), & sk(k,j,i-2), u_comp, 0.5, ddx ) ENDDO CASE ( 3 ) DO k = nzb_s_inner(j,i) + 1, nzt v_comp = v(k,j+1,i) - v_gtrans flux_n(k) = v_comp * ( & 7. * ( sk(k,j+1,i) + sk(k,j,i) ) & - ( sk(k,j+2,i) + sk(k,j-1,i) ) ) * adv_sca_3 diss_n(k) = - abs(v_comp) * ( & 3. * ( sk(k,j+1,i) - sk(k,j,i) ) & - ( sk(k,j+2,i) - sk(k,j-1,i) ) ) * adv_sca_3 ENDDO CASE ( 4 ) DO k = nzb_s_inner(j,i) + 1, nzt v_comp = v(k,j+1,i) - v_gtrans flux_n(k) = v_comp * ( sk(k,j+1,i) + sk(k,j,i) ) * 0.5 diss_n(k) = diss_2nd( sk(k,j+1,i), sk(k,j+1,i), & sk(k,j,i), sk(k,j-1,i), & sk(k,j-2,i), v_comp, 0.5, ddy ) ENDDO CASE ( 5 ) DO k = nzb_w_inner(j,i) + 1, nzt u_comp = u(k,j,i+1) - u_gtrans flux_r(k) = u_comp * ( & 7. * ( sk(k,j,i+1) + sk(k,j,i) ) & - ( sk(k,j,i+2) + sk(k,j,i-1) ) ) * adv_sca_3 diss_r(k) = - abs(u_comp) * ( & 3. * ( sk(k,j,i+1) - sk(k,j,i) ) & - ( sk(k,j,i+2) - sk(k,j,i-1) ) ) * adv_sca_3 ENDDO CASE ( 6 ) DO k = nzb_s_inner(j,i) + 1, nzt u_comp = u(k,j,i+1) - u_gtrans flux_r(k) = u_comp * ( & 7. * ( sk(k,j,i+1) + sk(k,j,i) ) & - ( sk(k,j,i+2) + sk(k,j,i-1) ) ) * adv_sca_3 diss_r(k) = - abs(u_comp) * ( & 3. * ( sk(k,j,i+1) - sk(k,j,i) ) & - ( sk(k,j,i+2) - sk(k,j,i-1) ) ) * adv_sca_3 ENDDO CASE ( 7 ) DO k = nzb_s_inner(j,i) + 1, nzt v_comp = v(k,j+1,i) - v_gtrans flux_n(k) = v_comp * ( & 7. * ( sk(k,j+1,i) + sk(k,j,i) ) & - ( sk(k,j+2,i) + sk(k,j-1,i) ) ) * adv_sca_3 diss_n(k) = - abs(v_comp) * ( & 3. * ( sk(k,j+1,i) - sk(k,j,i) ) & - ( sk(k,j+2,i) - sk(k,j-1,i) ) ) * adv_sca_3 ENDDO CASE ( 8 ) DO k = nzb_s_inner(j,i) + 1, nzt v_comp = v(k,j+1,i) - v_gtrans flux_n(k) = v_comp * ( & 7. * ( sk(k,j+1,i) + sk(k,j,i) ) & - ( sk(k,j+2,i) + sk(k,j-1,i) ) ) * adv_sca_3 diss_n(k) = - abs(v_comp) * ( & 3. * ( sk(k,j+1,i) - sk(k,j,i) ) & - ( sk(k,j+2,i) - sk(k,j-1,i) ) ) * adv_sca_3 ENDDO CASE DEFAULT END SELECT IF (boundary_flags(j,i) == 1 .OR. boundary_flags(j,i) == 2 & .OR. boundary_flags(j,i) == 5 .OR. boundary_flags(j,i) == 6) THEN DO k = nzb_s_inner(j,i) + 1, nzt v_comp = v(k,j+1,i) - v_gtrans flux_n(k) = v_comp * ( & 37. * ( sk(k,j+1,i) + sk(k,j,i) ) & - 8. * (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. * ( sk(k,j+1,i) - sk(k,j,i) ) & - 5. * ( sk(k,j+2,i) - sk(k,j-1,i) ) & + ( sk(k,j+3,i) - sk(k,j-2,i) ) ) & * adv_sca_5 ENDDO ELSE DO k = nzb_s_inner(j,i) + 1, nzt u_comp = u(k,j,i+1) - u_gtrans flux_r(k) = u_comp * ( & 37. * ( sk(k,j,i+1) + sk(k,j,i) ) & - 8. * ( 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. * ( sk(k,j,i+1) - sk(k,j,i) ) & - 5. * ( sk(k,j,i+2) - sk(k,j,i-1) ) & + ( sk(k,j,i+3) - sk(k,j,i-2) ) ) & * adv_sca_5 ENDDO ENDIF ELSE DO k = nzb_s_inner(j,i) + 1, nzt u_comp = u(k,j,i+1) - u_gtrans flux_r(k) = u_comp * ( & 37. * ( sk(k,j,i+1) + sk(k,j,i) ) & - 8. * ( 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. * ( sk(k,j,i+1) - sk(k,j,i) ) & - 5. * ( 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. * ( sk(k,j+1,i) + sk(k,j,i) ) & - 8. * ( 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. * ( sk(k,j+1,i) - sk(k,j,i) ) & - 5. * ( sk(k,j+2,i) - sk(k,j-1,i) ) & + ( sk(k,j+3,i) - sk(k,j-2,i) ) ) & * adv_sca_5 ENDDO ENDIF DO k = nzb_s_inner(j,i) + 1, nzt 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) swap_flux_x_local(k,j) = flux_r(k) swap_diss_x_local(k,j) = diss_r(k) swap_flux_y_local(k) = flux_n(k) swap_diss_y_local(k) = diss_n(k) ENDDO ENDDO ENDDO ! !-- Vertical advection, degradation of order near surface and top. !-- The fluxes flux_d and diss_d at the surface are 0. Due to reasons of !-- statistical evaluation the top flux at the surface should be 0 DO i = nxl, nxr DO j = nys, nyn ! !-- 2nd order scheme k=nzb_s_inner(j,i) + 1 ! !-- The fluxes flux_d and diss_d at the surface are 0. Due to static !-- reasons the top flux at the surface should be 0. flux_t(nzb_s_inner(j,i)) = 0. diss_t(nzb_s_inner(j,i)) = 0. flux_d = flux_t(k-1) diss_d = diss_t(k-1) flux_t(k) = w(k,j,i) * (sk(k+1,j,i) + sk(k,j,i) ) *0.5 diss_t(k) = diss_2nd( sk(k+2,j,i), sk(k+1,j,i), sk(k,j,i), & sk(k,j,i), sk(k,j,i), w(k,j,i), & 0.5, ddzw(k) ) tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) & - ( flux_d + diss_d ) ) * ddzw(k) ! !-- WS3 as an intermediate step k = nzb_s_inner(j,i) + 2 flux_d = flux_t(k-1) diss_d = diss_t(k-1) flux_t(k) = w(k,j,i) * ( & 7. * ( sk(k+1,j,i) + sk(k,j,i) ) & - ( sk(k+2,j,i) + sk(k-1,j,i) ) ) * adv_sca_3 diss_t(k) = - abs(w(k,j,i)) * ( & 3. * ( sk(k+1,j,i) - sk(k,j,i) ) & - ( sk(k+2,j,i) - sk(k-1,j,i) ) ) * adv_sca_3 tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) & - ( flux_d + diss_d ) ) * ddzw(k) ! !-- WS5 DO k = nzb_s_inner(j,i) + 3, nzt - 2 flux_d = flux_t(k-1) diss_d = diss_t(k-1) flux_t(k) = w(k,j,i) * ( & 37. * ( sk(k+1,j,i) + sk(k,j,i) ) & - 8. * ( sk(k+2,j,i) + sk(k-1,j,i) ) & + ( sk(k+3,j,i) + sk(k-2,j,i) ) ) * adv_sca_5 diss_t(k) = - abs(w(k,j,i)) * ( & 10. * ( sk(k+1,j,i) -sk(k,j,i) ) & - 5. * ( sk(k+2,j,i) - sk(k-1,j,i) ) & + ( sk(k+3,j,i) - sk(k-2,j,i) ) ) * adv_sca_5 tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) & - ( flux_d + diss_d ) ) * ddzw(k) ENDDO ! !-- WS3 as an intermediate step k = nzt - 1 flux_d = flux_t(k-1) diss_d = diss_t(k-1) flux_t(k) = w(k,j,i) * ( & 7. * ( sk(k+1,j,i) + sk(k,j,i) ) & - ( sk(k+2,j,i) + sk(k-1,j,i) ) ) * adv_sca_3 diss_t(k) = - abs(w(k,j,i)) * ( & 3. * ( sk(k+1,j,i) - sk(k,j,i) ) & - ( sk(k+2,j,i) - sk(k-1,j,i) ) ) * adv_sca_3 tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) & - ( flux_d + diss_d ) ) * ddzw(k) ! !-- 2nd order scheme k = nzt flux_d = flux_t(k-1) diss_d = diss_t(k-1) flux_t(k) = w(k,j,i) * ( sk(k+1,j,i) + sk(k,j,i) ) * 0.5 diss_t(k) = diss_2nd( sk(k+1,j,i), sk(k+1,j,i), sk(k,j,i), & sk(k-1,j,i), sk(k-2,j,i), w(k,j,i), & 0.5, ddzw(k) ) tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) & - ( flux_d + diss_d ) ) * ddzw(k) ! !-- evaluation of statistics SELECT CASE ( sk_char ) CASE ( 'pt' ) DO k = nzb_s_inner(j,i), nzt sums_wspts_ws_l(k,:) = sums_wspts_ws_l(k,:) & + ( flux_t(k) + diss_t(k) ) & * weight_substep(intermediate_timestep_count) & * rmask(j,i,:) ENDDO CASE ( 'sa' ) DO k = nzb_s_inner(j,i), nzt sums_wssas_ws_l(k,:) = sums_wssas_ws_l(k,:) & + ( flux_t(k) + diss_t(k) ) & * weight_substep(intermediate_timestep_count) & * rmask(j,i,:) ENDDO CASE ( 'q' ) DO k = nzb_s_inner(j,i), nzt sums_wsqs_ws_l(k,:) = sums_wsqs_ws_l(k,:) & + ( flux_t(k) + diss_t(k) ) & * weight_substep(intermediate_timestep_count) & * rmask(j,i,:) 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, j, k REAL :: gu, gv, flux_d, diss_d, v_comp, w_comp REAL, DIMENSION(nzb+1:nzt) :: swap_flux_y_local_u, swap_diss_y_local_u REAL, DIMENSION(nzb+1:nzt,nys:nyn) :: swap_flux_x_local_u, & swap_diss_x_local_u REAL, DIMENSION(nzb:nzt+1) :: flux_t, diss_t, flux_r, diss_r, flux_n, & diss_n, 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 IF( boundary_flags(j,i) == 5 ) THEN DO k = nzb_u_inner(j,i) + 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) * & ( u(k,j,i) + u(k,j,i-1) ) * 0.25 swap_diss_x_local_u(k,j) = diss_2nd( u(k,j,i+2), u(k,j,i+1), & u(k,j,i), u(k,j,i-1), & u(k,j,i-1), u_comp(k), & 0.25, ddx ) ENDDO ELSE DO k = nzb_u_inner(j,i) + 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. * ( u(k,j,i) + u(k,j,i-1) ) & - 8. * ( 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. * ( u(k,j,i) - u(k,j,i-1) ) & - 5. * ( 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 ENDDO DO i = nxlu, nxr ! !-- The following loop computes the fluxes for the south boundary points j = nys IF ( boundary_flags(j,i) == 8 ) THEN ! !-- Compute southside fluxes for the south boundary of PE domain DO k = nzb_u_inner(j,i) + 1, nzt v_comp = v(k,j,i) + v(k,j,i-1) - gv swap_flux_y_local_u(k) = v_comp * & ( u(k,j,i) + u(k,j-1,i) ) * 0.25 swap_diss_y_local_u(k) = diss_2nd( u(k,j+2,i), u(k,j+1,i), & u(k,j,i), u(k,j-1,i), & u(k,j-1,i), v_comp, & 0.25, ddy ) ENDDO ELSE DO k = nzb_u_inner(j,i) + 1, nzt v_comp = v(k,j,i) + v(k,j,i-1) - gv swap_flux_y_local_u(k) = v_comp * ( & 37. * ( u(k,j,i) + u(k,j-1,i) ) & - 8. * ( 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. * ( u(k,j,i) - u(k,j-1,i) ) & - 5. * ( 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 ! !-- Computation of interior fluxes and tendency terms DO j = nys, nyn IF ( boundary_flags(j,i) /= 0 ) THEN ! !-- Degrade the order for Dirichlet bc. at the outflow boundary SELECT CASE ( boundary_flags(j,i) ) CASE ( 1 ) DO k = nzb_u_inner(j,i) + 1, nzt u_comp(k) = u(k,j,i+1) + u(k,j,i) flux_r(k) = ( u_comp(k) - gu ) * ( & 7. * ( u(k,j,i+1) + u(k,j,i) ) & - ( u(k,j,i+2) + u(k,j,i-1) ) ) * adv_mom_3 diss_r(k) = - abs(u_comp(k) - gu) * ( & 3. * ( u(k,j,i+1) - u(k,j,i) ) & - ( u(k,j,i+2) - u(k,j,i-1) ) ) * adv_mom_3 ENDDO CASE ( 2 ) DO k = nzb_u_inner(j,i) + 1, nzt u_comp(k) = u(k,j,i+1) + u(k,j,i) flux_r(k) = ( u_comp(k) - gu ) * & ( u(k,j,i+1) + u(k,j,i) ) * 0.25 diss_r(k) = diss_2nd( u(k,j,i+1), u(k,j,i+1), & u(k,j,i), u(k,j,i-1), & u(k,j,i-2), u_comp(k) ,0.25 ,ddx) ENDDO CASE ( 3 ) DO k = nzb_u_inner(j,i) + 1, nzt v_comp = v(k,j+1,i) + v(k,j+1,i-1) - gv flux_n(k) = v_comp * ( & 7. * ( u(k,j+1,i) + u(k,j,i) ) & - ( u(k,j+2,i) + u(k,j-1,i) ) ) * adv_mom_3 diss_n(k) = - abs(v_comp) * ( & 3. * ( u(k,j+1,i) - u(k,j,i) ) & - ( u(k,j+2,i) - u(k,j-1,i) ) ) * adv_mom_3 ENDDO CASE ( 4 ) DO k = nzb_u_inner(j,i) + 1, nzt v_comp = v(k,j+1,i) + v(k,j+1,i-1) - gv flux_n(k) = v_comp * ( u(k,j+1,i) + u(k,j,i) ) * 0.25 diss_n(k) = diss_2nd( u(k,j+1,i), u(k,j+1,i), & u(k,j,i), u(k,j-1,i), & u(k,j-2,i), v_comp, 0.25, ddy ) ENDDO CASE ( 5 ) DO k = nzb_u_inner(j,i) + 1, nzt u_comp(k) = u(k,j,i+1) + u(k,j,i) flux_r(k) = ( u_comp(k) - gu ) * ( & 7. * ( u(k,j,i+1) + u(k,j,i) ) & - ( u(k,j,i+2) + u(k,j,i-1) ) ) * adv_mom_3 diss_r(k) = - abs(u_comp(k) - gu) * ( & 3. * ( u(k,j,i+1) - u(k,j,i) ) & - ( u(k,j,i+2) - u(k,j,i-1) ) ) * adv_mom_3 ENDDO CASE ( 7 ) DO k = nzb_u_inner(j,i) + 1, nzt v_comp = v(k,j+1,i) + v(k,j+1,i-1) - gv flux_n(k) = v_comp * ( & 7. * ( u(k,j+1,i) + u(k,j,i) ) & - ( u(k,j+2,i) + u(k,j-1,i) ) ) * adv_mom_3 diss_n(k) = - abs(v_comp) * ( & 3. * ( u(k,j+1,i) - u(k,j,i) ) & - ( u(k,j+2,i) - u(k,j-1,i) ) ) * adv_mom_3 ENDDO CASE ( 8 ) DO k = nzb_u_inner(j,i) + 1, nzt v_comp = v(k,j+1,i) + v(k,j+1,i-1) - gv flux_n(k) = v_comp * ( & 7. * ( u(k,j+1,i) + u(k,j,i) ) & - ( u(k,j+2,i) + u(k,j-1,i) ) ) * adv_mom_3 diss_n(k) = - abs(v_comp) * ( & 3. * ( u(k,j+1,i) - u(k,j,i) ) & - ( u(k,j+2,i) - u(k,j-1,i) ) ) * adv_mom_3 ENDDO CASE DEFAULT END SELECT ! !-- Compute the crosswise 5th order fluxes at the outflow IF (boundary_flags(j,i) == 1 .OR. boundary_flags(j,i) == 2 & .OR. boundary_flags(j,i) == 5) THEN DO k = nzb_u_inner(j,i) + 1, nzt v_comp = v(k,j+1,i) + v(k,j+1,i-1) - gv flux_n(k) = v_comp * ( & 37. * ( u(k,j+1,i) + u(k,j,i) ) & - 8. * ( 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. * ( u(k,j+1,i) - u(k,j,i) ) & - 5. * ( u(k,j+2,i) - u(k,j-1,i) ) & + ( u(k,j+3,i) - u(k,j-2,i) ) ) * adv_mom_5 ENDDO ELSE DO k = nzb_u_inner(j,i) + 1, nzt u_comp(k) = u(k,j,i+1) + u(k,j,i) flux_r(k) = ( u_comp(k) - gu ) * ( & 37. * ( u(k,j,i+1) + u(k,j,i) ) & - 8. * ( 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. * ( u(k,j,i+1) - u(k,j,i) ) & - 5. * ( u(k,j,i+2) - u(k,j,i-1) ) & + ( u(k,j,i+3) - u(k,j,i-2) ) ) * adv_mom_5 ENDDO ENDIF ELSE DO k = nzb_u_inner(j,i) + 1, nzt u_comp(k) = u(k,j,i+1) + u(k,j,i) flux_r(k) = ( u_comp(k) - gu ) * ( & 37. * ( u(k,j,i+1) + u(k,j,i) ) & - 8. * ( 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. * ( u(k,j,i+1) - u(k,j,i) ) & - 5. * ( 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. * ( u(k,j+1,i) + u(k,j,i) ) & - 8. * ( 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. * ( u(k,j+1,i) - u(k,j,i) ) & - 5. * ( u(k,j+2,i) - u(k,j-1,i) ) & + ( u(k,j+3,i) - u(k,j-2,i) ) ) * adv_mom_5 ENDDO ENDIF DO k = nzb_u_inner(j,i) + 1, nzt 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 ) 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) sums_us2_ws_l(k,:) = sums_us2_ws_l(k,:) & + ( flux_r(k) & * ( u_comp(k) - 2. * hom(k,1,1,:) ) & / ( u_comp(k) - gu + 1.0E-20 ) & + diss_r(k) & * abs(u_comp(k) - 2. * hom(k,1,1,:) ) & / (abs(u_comp(k) - gu) + 1.0E-20) ) & * weight_substep(intermediate_timestep_count) * rmask(j,i,:) ENDDO sums_us2_ws_l(nzb_u_inner(j,i),:) = & sums_us2_ws_l(nzb_u_inner(j,i)+1,:) ENDDO ENDDO ! !-- Vertical advection, degradation of order near surface and top. !-- The fluxes flux_d and diss_d at the surface are 0. Due to reasons of !-- statistical evaluation the top flux at the surface should be 0 DO i = nxlu, nxr DO j = nys, nyn k = nzb_u_inner(j,i) + 1 ! !-- The fluxes flux_d and diss_d at the surface are 0. Due to static !-- reasons the top flux at the surface should be 0. flux_t(nzb_u_inner(j,i)) = 0. diss_t(nzb_u_inner(j,i)) = 0. flux_d = flux_t(k-1) diss_d = diss_t(k-1) ! !-- 2nd order scheme w_comp = w(k,j,i) + w(k,j,i-1) flux_t(k) = w_comp * ( u(k+1,j,i) + u(k,j,i) ) * 0.25 diss_t(k) = diss_2nd( u(k+2,j,i), u(k+1,j,i), u(k,j,i), & 0., 0., w_comp, 0.25, ddzw(k) ) tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) & - ( flux_d + diss_d ) ) * ddzw(k) ! !-- WS3 as an intermediate step k = nzb_u_inner(j,i) + 2 flux_d = flux_t(k-1) diss_d = diss_t(k-1) w_comp = w(k,j,i) + w(k,j,i-1) flux_t(k) = w_comp * ( & 7. * (u(k+1,j,i) + u(k,j,i) ) & - ( u(k+2,j,i) + u(k-1,j,i) ) ) * adv_mom_3 diss_t(k) = - abs(w_comp) * ( & 3. * ( u(k+1,j,i) - u(k,j,i) ) & - ( u(k+2,j,i) - u(k-1,j,i) ) ) * adv_mom_3 tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) & - ( flux_d + diss_d ) ) * ddzw(k) ! !WS5 DO k = nzb_u_inner(j,i) + 3, nzt - 3 flux_d = flux_t(k-1) diss_d = diss_t(k-1) w_comp = w(k,j,i) + w(k,j,i-1) flux_t(k) = w_comp * ( & 37. * ( u(k+1,j,i) + u(k,j,i) ) & - 8. * ( u(k+2,j,i) + u(k-1,j,i) ) & + ( u(k+3,j,i) + u(k-2,j,i) ) ) * adv_mom_5 diss_t(k) = - abs(w_comp) * ( & 10. * ( u(k+1,j,i) - u(k,j,i) ) & - 5. * ( u(k+2,j,i) - u(k-1,j,i) ) & + ( u(k+3,j,i) - u(k-2,j,i) ) ) * adv_mom_5 tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) & - ( flux_d + diss_d ) ) * ddzw(k) ENDDO ! !-- WS3 as an intermediate step DO k = nzt - 2, nzt - 1 flux_d = flux_t(k-1) diss_d = diss_t(k-1) w_comp = w(k,j,i) + w(k,j,i-1) flux_t(k) = w_comp * ( & 7. * ( u(k+1,j,i) + u(k,j,i) ) & - ( u(k+2,j,i) + u(k-1,j,i) ) ) * adv_mom_3 diss_t(k) = - abs(w_comp) * ( & 3. * ( u(k+1,j,i) - u(k,j,i) ) & - ( u(k+2,j,i) - u(k-1,j,i) ) ) * adv_mom_3 tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) & - ( flux_d + diss_d ) ) * ddzw(k) ENDDO ! !-- 2nd order scheme k = nzt flux_d = flux_t(k-1) diss_d = diss_t(k-1) w_comp = w(k,j,i) + w(k,j,i-1) flux_t(k) = w_comp * ( u(k+1,j,i) + u(k,j,i) ) * 0.25 diss_t(k) = diss_2nd( u(nzt+1,j,i), u(nzt+1,j,i), u(k,j,i), & u(k-1,j,i), u(k-2,j,i), w_comp, & 0.25, ddzw(k)) tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) & - ( flux_d + diss_d ) ) * ddzw(k) ! !-- at last vertical momentum flux is accumulated DO k = nzb_u_inner(j,i), nzt sums_wsus_ws_l(k,:) = sums_wsus_ws_l(k,:) & + ( flux_t(k) + diss_t(k) ) & * weight_substep(intermediate_timestep_count) & * rmask(j,i,:) ENDDO ENDDO ENDDO 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, j, k REAL :: gu, gv, flux_l, flux_s, flux_d, diss_l, diss_s, diss_d, & u_comp, w_comp REAL, DIMENSION(nzb+1:nzt) :: swap_flux_y_local_v, swap_diss_y_local_v REAL, DIMENSION(nzb+1:nzt,nys:nyn) :: swap_flux_x_local_v, & swap_diss_x_local_v REAL, DIMENSION(nzb:nzt+1) :: flux_t, diss_t, flux_n, diss_n, flux_r, & diss_r, 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 IF ( boundary_flags(j,i) == 6 ) THEN DO k = nzb_v_inner(j,i) + 1, nzt u_comp = u(k,j-1,i) + u(k,j,i) - gu swap_flux_x_local_v(k,j) = u_comp * & ( v(k,j,i) + v(k,j,i-1)) * 0.25 swap_diss_x_local_v(k,j) = diss_2nd( v(k,j,i+2), v(k,j,i+1), & v(k,j,i), v(k,j,i-1), & v(k,j,i-1), u_comp, & 0.25, ddx ) ENDDO ELSE DO k = nzb_v_inner(j,i)+1, nzt u_comp = u(k,j-1,i) + u(k,j,i) - gu swap_flux_x_local_v(k,j) = u_comp * ( & 37. * ( v(k,j,i) + v(k,j,i-1) ) & - 8. * ( 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. * ( v(k,j,i) - v(k,j,i-1) ) & -5. * ( 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 ENDDO DO i = nxl, nxr j = nysv IF ( boundary_flags(j,i) == 7 ) THEN DO k = nzb_v_inner(j,i) + 1, nzt v_comp(k) = v(k,j,i) + v(k,j-1,i) - gv swap_flux_y_local_v(k) = v_comp(k) * & ( v(k,j,i) + v(k,j-1,i) ) * 0.25 swap_diss_y_local_v(k) = diss_2nd( v(k,j+2,i), v(k,j+1,i), & v(k,j,i), v(k,j-1,i), & v(k,j-1,i), v_comp(k), & 0.25, ddy ) ENDDO ELSE DO k = nzb_v_inner(j,i) + 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. * ( v(k,j,i) + v(k,j-1,i) ) & - 8. * ( 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. * ( v(k,j,i) - v(k,j-1,i) ) & - 5. * ( 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 DO j = nysv, nyn IF ( boundary_flags(j,i) /= 0 ) THEN ! !-- Degrade the order for Dirichlet bc. at the outflow boundary SELECT CASE ( boundary_flags(j,i) ) CASE ( 1 ) DO k = nzb_v_inner(j,i) + 1, nzt u_comp = u(k,j-1,i+1) + u(k,j,i+1) - gu flux_r(k) = u_comp * ( & 7. * (v(k,j,i+1) + v(k,j,i) ) & - ( v(k,j,i+2) + v(k,j,i-1) ) ) * adv_mom_3 diss_r(k) = - abs(u_comp) * ( & 3. * ( v(k,j,i+1) - v(k,j,i) ) & - ( v(k,j,i+2) - v(k,j,i-1) ) ) * adv_mom_3 ENDDO CASE ( 2 ) DO k = nzb_v_inner(j,i) + 1, nzt u_comp = u(k,j-1,i+1) + u(k,j,i+1) - gu flux_r(k) = u_comp * ( v(k,j,i+1) + v(k,j,i) ) * 0.25 diss_r(k) = diss_2nd( v(k,j,i+1), v(k,j,i+1), & v(k,j,i), v(k,j,i-1), & v(k,j,i-2), u_comp, 0.25, ddx ) ENDDO CASE ( 3 ) DO k = nzb_v_inner(j,i) + 1, nzt v_comp(k) = v(k,j+1,i) + v(k,j,i) flux_n(k) = ( v_comp(k)- gv ) * ( & 7. * ( v(k,j+1,i) + v(k,j,i) ) & - ( v(k,j+2,i) + v(k,j-1,i) ) ) * adv_mom_3 diss_n(k) = - abs(v_comp(k) - gv) * ( & 3. * ( v(k,j+1,i) - v(k,j,i) ) & - ( v(k,j+2,i) - v(k,j-1,i) ) ) * adv_mom_3 ENDDO CASE ( 4 ) DO k = nzb_v_inner(j,i) + 1, nzt v_comp(k) = v(k,j+1,i) + v(k,j,i) flux_n(k) = ( v_comp(k) - gv ) * & ( v(k,j+1,i) + v(k,j,i) ) * 0.25 diss_n(k) = diss_2nd( v(k,j+1,i), v(k,j+1,i), & v(k,j,i), v(k,j-1,i), & v(k,j-2,i), v_comp(k), 0.25, ddy) ENDDO CASE ( 5 ) DO k = nzb_w_inner(j,i) + 1, nzt u_comp = u(k,j-1,i) + u(k,j,i) - gu flux_r(k) = u_comp * ( & 7. * ( v(k,j,i+1) + v(k,j,i) ) & - ( v(k,j,i+2) + v(k,j,i-1) ) ) * adv_mom_3 diss_r(k) = - abs(u_comp) * ( & 3. * ( w(k,j,i+1) - w(k,j,i) ) & - ( v(k,j,i+2) - v(k,j,i-1) ) ) * adv_mom_3 ENDDO CASE ( 6 ) DO k = nzb_v_inner(j,i) + 1, nzt u_comp = u(k,j-1,i+1) + u(k,j,i+1) - gu flux_r(k) = u_comp * ( & 7. * ( v(k,j,i+1) + v(k,j,i) ) & - ( v(k,j,i+2) + v(k,j,i-1) ) ) * adv_mom_3 diss_r(k) = - abs(u_comp) * ( & 3. * ( v(k,j,i+1) - v(k,j,i) ) & - ( v(k,j,i+2) - v(k,j,i-1) ) ) * adv_mom_3 ENDDO CASE ( 7 ) DO k = nzb_v_inner(j,i) + 1, nzt v_comp(k) = v(k,j+1,i) + v(k,j,i) flux_n(k) = ( v_comp(k) - gv ) * ( & 7. * ( v(k,j+1,i) + v(k,j,i) ) & - ( v(k,j+2,i) + v(k,j-1,i) ) ) * adv_mom_3 diss_n(k) = - abs(v_comp(k) - gv) * ( & 3. * ( v(k,j+1,i) - v(k,j,i) ) & - ( v(k,j+2,i) - v(k,j-1,i) ) ) * adv_mom_3 ENDDO CASE DEFAULT END SELECT ! !-- Compute the crosswise 5th order fluxes at the outflow IF ( boundary_flags(j,i) == 1 .OR. boundary_flags(j,i) == 2 & .OR. boundary_flags(j,i) == 5 .OR. boundary_flags(j,i) == 6 ) THEN DO k = nzb_v_inner(j,i) + 1, nzt v_comp(k) = v(k,j+1,i) + v(k,j,i) flux_n(k) = ( v_comp(k) - gv ) * ( & 37. * ( v(k,j+1,i) + v(k,j,i) ) & - 8. * ( 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. * ( v(k,j+1,i) - v(k,j,i) ) & - 5. * ( v(k,j+2,i) - v(k,j-1,i) ) & + ( v(k,j+3,i) - v(k,j-2,i) ) ) *adv_mom_5 ENDDO ELSE DO k = nzb_v_inner(j,i) + 1, nzt u_comp = u(k,j-1,i+1) + u(k,j,i+1) - gu flux_r(k) = u_comp * ( & 37. * ( v(k,j,i+1) + v(k,j,i) ) & - 8. * ( 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. * ( v(k,j,i+1) - v(k,j,i) ) & -5. * ( v(k,j,i+2) - v(k,j,i-1) ) & + ( v(k,j,i+3) - v(k,j,i-2) ) ) * adv_mom_5 ENDDO ENDIF ELSE DO k = nzb_v_inner(j,i) + 1, nzt u_comp = u(k,j-1,i+1) + u(k,j,i+1) - gu flux_r(k) = u_comp * ( & 37. * ( v(k,j,i+1) + v(k,j,i) ) & - 8. * ( 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. * ( v(k,j,i+1) - v(k,j,i) ) & -5. * ( 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. * ( v(k,j+1,i) + v(k,j,i) ) & - 8. * ( 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. * ( v(k,j+1,i) - v(k,j,i) ) & - 5. * ( v(k,j+2,i) - v(k,j-1,i) ) & + ( v(k,j+3,i) - v(k,j-2,i) ) ) *adv_mom_5 ENDDO ENDIF DO k = nzb_v_inner(j,i) + 1, nzt 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 ) 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) sums_vs2_ws_l(k,:) = sums_vs2_ws_l(k,:) & + (flux_n(k) *( v_comp(k) - 2. * hom(k,1,2,:) ) & / ( v_comp(k) - gv + 1.0E-20) & + diss_n(k) * abs(v_comp(k) - 2. * hom(k,1,2,:) ) & /(abs(v_comp(k) - gv) + 1.0E-20 ) ) & * weight_substep(intermediate_timestep_count) * rmask(j,i,:) ENDDO sums_vs2_ws_l(nzb_v_inner(j,i),:) = & sums_vs2_ws_l(nzb_v_inner(j,i)+1,:) ENDDO ENDDO ! !-- Vertical advection, degradation of order near surface and top. !-- The fluxes flux_d and diss_d at the surface are 0. Due to reasons of !-- statistical evaluation the top flux at the surface should be 0 DO i = nxl, nxr DO j = nysv, nyn ! !-- The fluxes flux_d and diss_d at the surface are 0. Due to static !-- reasons the top flux at the surface should be 0. flux_t(nzb_v_inner(j,i)) = 0. diss_t(nzb_v_inner(j,i)) = 0. k = nzb_v_inner(j,i) + 1 flux_d = flux_t(k-1) diss_d = diss_t(k-1) ! !-- 2nd order scheme w_comp = w(k,j-1,i) + w(k,j,i) flux_t(k) = w_comp * ( v(k+1,j,i) + v(k,j,i) ) * 0.25 diss_t(k) = diss_2nd( v(k+2,j,i), v(k+1,j,i), v(k,j,i), & 0., 0., w_comp, 0.25, ddzw(k) ) tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) & - ( flux_d + diss_d ) ) * ddzw(k) ! !-- WS3 as an intermediate step k = nzb_v_inner(j,i) + 2 flux_d = flux_t(k-1) diss_d = diss_t(k-1) w_comp = w(k,j-1,i) + w(k,j,i) flux_t(k) = w_comp * ( & 7. * ( v(k+1,j,i) + v(k,j,i) ) & - ( v(k+2,j,i) + v(k-1,j,i) ) ) * adv_mom_3 diss_t(k) = - abs(w_comp) * ( & 3. * ( v(k+1,j,i) - v(k,j,i) ) & - ( v(k+2,j,i) - v(k-1,j,i) ) ) * adv_mom_3 tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) & - ( flux_d + diss_d ) ) * ddzw(k) ! !-- WS5 DO k = nzb_v_inner(j,i) + 3, nzt - 3 flux_d = flux_t(k-1) diss_d = diss_t(k-1) w_comp = w(k,j-1,i) + w(k,j,i) flux_t(k) = w_comp * ( & 37. * ( v(k+1,j,i) + v(k,j,i) ) & - 8. * ( v(k+2,j,i) + v(k-1,j,i) ) & + ( v(k+3,j,i) + v(k-2,j,i) ) ) * adv_mom_5 diss_t(k) = - abs(w_comp) * ( & 10. * ( v(k+1,j,i) - v(k,j,i) ) & -5. * ( v(k+2,j,i) - v(k-1,j,i) ) & + ( v(k+3,j,i) - v(k-2,j,i) ) ) * adv_mom_5 tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) & - ( flux_d + diss_d ) ) * ddzw(k) ENDDO ! !-- WS3 as an intermediate step DO k = nzt - 2, nzt - 1 flux_d = flux_t(k-1) diss_d = diss_t(k-1) w_comp = w(k,j-1,i) + w(k,j,i) flux_t(k) = w_comp * ( & 7. * ( v(k+1,j,i) + v(k,j,i) ) & - ( v(k+2,j,i) + v(k-1,j,i) ) ) * adv_mom_3 diss_t(k) = - abs(w_comp) * ( & 3. * ( v(k+1,j,i) - v(k,j,i) ) & - ( v(k+2,j,i) - v(k-1,j,i) ) ) * adv_mom_3 tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) & - ( flux_d + diss_d ) ) * ddzw(k) ENDDO ! !-- 2nd order scheme k = nzt flux_d = flux_t(k-1) diss_d = diss_t(k-1) w_comp = w(k,j-1,i) + w(k,j,i) flux_t(k) = w_comp * ( v(k+1,j,i) + v(k,j,i) ) * 0.25 diss_t(k) = diss_2nd( v(nzt+1,j,i), v(nzt+1,j,i), v(k,j,i), & v(k-1,j,i), v(k-2,j,i), w_comp, 0.25,ddzw(k)) tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) & - ( flux_d + diss_d ) ) * ddzw(k) ! !- at last vertical momentum flux is accumulated DO k = nzb_v_inner(j,i), nzt sums_wsvs_ws_l(k,:) = sums_wsvs_ws_l(k,:) & + ( flux_t(k) + diss_t(k) ) & * weight_substep(intermediate_timestep_count) & * rmask(j,i,:) ENDDO sums_vs2_ws_l(nzb_v_inner(j,i),:) = & sums_vs2_ws_l(nzb_v_inner(j,i)+1,:) ENDDO ENDDO 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, j, k REAL :: gu, gv, flux_d, diss_d, u_comp, v_comp, w_comp REAL :: flux_t(nzb:nzt+1), diss_t(nzb:nzt+1) REAL, DIMENSION(nzb:nzt+1) :: flux_r, diss_r, flux_n, diss_n REAL, DIMENSION(nzb+1:nzt) :: swap_flux_y_local_w, swap_diss_y_local_w REAL, DIMENSION(nzb+1:nzt,nys:nyn) :: swap_flux_x_local_w, & swap_diss_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 IF ( boundary_flags(j,i) == 6 ) THEN DO k = nzb_v_inner(j,i) + 1, nzt u_comp = u(k+1,j,i) + u(k,j,i) - gu swap_flux_x_local_w(k,j) = u_comp * & (w(k,j,i)+w(k,j,i-1)) * 0.25 swap_flux_x_local_w(k,j) = diss_2nd( w(k,j,i+2), w(k,j,i+1), & w(k,j,i), w(k,j,i-1), & w(k,j,i-1), u_comp, & 0.25, ddx ) ENDDO ELSE DO k = nzb_v_inner(j,i) + 1, nzt u_comp = u(k+1,j,i) + u(k,j,i) - gu swap_flux_x_local_w(k,j) = u_comp * ( & 37. * ( w(k,j,i) + w(k,j,i-1) ) & - 8. * ( 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. * ( w(k,j,i) - w(k,j,i-1) ) & - 5. * ( 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 ENDDO DO i = nxl, nxr j = nys IF ( boundary_flags(j,i) == 8 ) THEN DO k = nzb_v_inner(j,i) + 1, nzt v_comp = v(k+1,j,i) + v(k,j,i) - gv swap_flux_y_local_w(k) = v_comp * & ( w(k,j,i) + w(k,j-1,i) ) * 0.25 swap_diss_y_local_w(k) = diss_2nd( w(k,j+2,i), w(k,j+1,i), & w(k,j,i), w(k,j-1,i), & w(k,j-1,i), v_comp, 0.25,ddy) ENDDO ELSE DO k = nzb_v_inner(j,i) + 1, nzt v_comp = v(k+1,j,i) + v(k,j,i) - gv swap_flux_y_local_w(k) = v_comp * ( & 37. * ( w(k,j,i) + w(k,j-1,i) ) & - 8. * ( 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. * ( w(k,j,i) - w(k,j-1,i) ) & - 5. * ( 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 DO j = nys, nyn IF ( boundary_flags(j,i) /= 0 ) THEN ! !-- Degrade the order for Dirichlet bc. at the outflow boundary SELECT CASE ( boundary_flags(j,i) ) CASE ( 1 ) DO k = nzb_w_inner(j,i) + 1, nzt u_comp = u(k+1,j,i+1) + u(k,j,i+1) - gu flux_r(k) = u_comp * ( & 7. * ( w(k,j,i+1) + w(k,j,i) ) & - ( w(k,j,i+2) + w(k,j,i-1) ) ) * adv_mom_3 diss_r(k) = -abs(u_comp) * ( & 3. * ( w(k,j,i+1) - w(k,j,i) ) & - ( w(k,j,i+2) - w(k,j,i-1) ) ) * adv_mom_3 ENDDO CASE ( 2 ) DO k = nzb_w_inner(j,i) + 1, nzt u_comp = u(k+1,j,i+1) + u(k,j,i+1) - gu flux_r(k) = u_comp * ( w(k,j,i+1) + w(k,j,i) ) * 0.25 diss_r(k) = diss_2nd( w(k,j,i+1), w(k,j,i+1), & w(k,j,i), w(k,j,i-1), & w(k,j,i-2), u_comp, 0.25, ddx ) ENDDO CASE ( 3 ) DO k = nzb_w_inner(j,i) + 1, nzt v_comp = v(k+1,j+1,i) + v(k,j+1,i) - gv flux_n(k) = v_comp * ( & 7. * ( w(k,j+1,i) + w(k,j,i) ) & - ( w(k,j+2,i) + w(k,j-1,i) ) ) * adv_mom_3 diss_n(k) = -abs(v_comp) * ( & 3. * ( w(k,j+1,i) - w(k,j,i) ) & - ( w(k,j+2,i) - w(k,j-1,i) ) ) * adv_mom_3 ENDDO CASE ( 4 ) DO k = nzb_w_inner(j,i) + 1, nzt v_comp = v(k+1,j+1,i) + v(k,j+1,i) - gv flux_n(k) = v_comp * ( w(k,j+1,i) + w(k,j,i) ) * 0.25 diss_n(k) = diss_2nd( w(k,j+1,i), w(k,j+1,i), & w(k,j,i), w(k,j-1,i), & w(k,j-2,i), v_comp, 0.25, ddy ) ENDDO CASE ( 5 ) DO k = nzb_w_inner(j,i) + 1, nzt u_comp = u(k+1,j,i+1) + u(k,j,i+1) - gu flux_r(k) = u_comp * ( & 7. * ( w(k,j,i+1) + w(k,j,i) ) & - ( w(k,j,i+2) + w(k,j,i-1) ) ) * adv_mom_3 diss_r(k) = - abs(u_comp) * ( & 3. * ( w(k,j,i+1) - w(k,j,i) ) & - ( w(k,j,i+2) - w(k,j,i-1) ) ) * adv_mom_3 ENDDO CASE ( 6 ) DO k = nzb_w_inner(j,i) + 1, nzt u_comp = u(k+1,j,i+1) + u(k,j,i+1) - gu flux_r(k) = u_comp *( & 7. * ( w(k,j,i+1) + w(k,j,i) ) & - ( w(k,j,i+2) + w(k,j,i-1) ) ) * adv_mom_3 diss_r(k) = - abs(u_comp) * ( & 3. * ( w(k,j,i+1) - w(k,j,i) ) & - ( w(k,j,i+2) - w(k,j,i-1) ) ) * adv_mom_3 ENDDO CASE ( 7 ) DO k = nzb_w_inner(j,i) + 1, nzt v_comp = v(k+1,j+1,i) + v(k,j+1,i) - gv flux_n(k) = v_comp *( & 7. * ( w(k,j+1,i) + w(k,j,i) ) & - ( w(k,j+2,i) + w(k,j-1,i) ) ) * adv_mom_3 diss_n(k) = - abs(v_comp) * ( & 3. * ( w(k,j+1,i) - w(k,j,i) ) & - ( w(k,j+2,i) - w(k,j-1,i) ) ) * adv_mom_3 ENDDO CASE ( 8 ) DO k = nzb_w_inner(j,i) + 1, nzt v_comp = v(k+1,j+1,i) + v(k,j+1,i) - gv flux_n(k) = v_comp * ( & 7. * ( w(k,j+1,i) + w(k,j,i) ) & - ( w(k,j+2,i) + w(k,j-1,i) ) ) * adv_mom_3 diss_n(k) = - abs(v_comp) * ( & 3. * ( w(k,j+1,i) - w(k,j,i) ) & - ( w(k,j+2,i) - w(k,j-1,i) ) ) * adv_mom_3 ENDDO CASE DEFAULT END SELECT ! !-- Compute the crosswise 5th order fluxes at the outflow IF ( boundary_flags(j,i) == 1 .OR. boundary_flags(j,i) == 2 & .OR. boundary_flags(j,i) == 5 .OR. boundary_flags(j,i) == 6 ) THEN DO k = nzb_w_inner(j,i) + 1, nzt v_comp = v(k+1,j+1,i) + v(k,j+1,i) - gv flux_n(k) = v_comp * ( & 37. * ( w(k,j+1,i) + w(k,j,i) ) & - 8. * ( 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. * ( w(k,j+1,i) - w(k,j,i) ) & - 5. * ( w(k,j+2,i) - w(k,j-1,i) ) & + ( w(k,j+3,i) - w(k,j-2,i) ) ) * adv_mom_5 ENDDO ELSE DO k = nzb_w_inner(j,i) + 1, nzt u_comp = u(k+1,j,i+1) + u(k,j,i+1) - gu flux_r(k) = u_comp * ( & 37. * ( w(k,j,i+1) + w(k,j,i) ) & - 8. * ( 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. * ( w(k,j,i+1) - w(k,j,i) ) & - 5. * ( w(k,j,i+2) - w(k,j,i-1) ) & + ( w(k,j,i+3) - w(k,j,i-2) ) ) * adv_mom_5 ENDDO ENDIF ELSE ! !-- Compute the fifth order fluxes for the interior of PE domain. DO k = nzb_w_inner(j,i) + 1, nzt u_comp = u(k+1,j,i+1) + u(k,j,i+1) - gu flux_r(k) = u_comp * ( & 37. * ( w(k,j,i+1) + w(k,j,i) ) & - 8. * ( 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. * ( w(k,j,i+1) - w(k,j,i) ) & - 5. * ( 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. * ( w(k,j+1,i) + w(k,j,i) ) & - 8. * ( 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. * ( w(k,j+1,i) - w(k,j,i) ) & - 5. * ( w(k,j+2,i) - w(k,j-1,i) ) & + ( w(k,j+3,i) - w(k,j-2,i) ) ) * adv_mom_5 ENDDO ENDIF DO k = nzb_v_inner(j,i)+1, nzt 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 ) 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) ENDDO ENDDO ENDDO ! !-- Vertical advection, degradation of order near surface and top. !-- The fluxes flux_d and diss_d at the surface are 0. Due to reasons of !-- statistical evaluation the top flux at the surface should be 0 DO i = nxl, nxr DO j = nys, nyn k = nzb_w_inner(j,i) + 1 flux_d = w(k-1,j,i) * ( w(k,j,i) + w(k-1,j,i)) flux_t(k-1) = flux_d diss_d = 0.0 diss_t(k-1) = diss_d ! !-- 2nd order scheme w_comp = w(k+1,j,i) + w(k,j,i) flux_t(k) = w_comp * ( w(k+1,j,i) + w(k,j,i) ) * 0.25 diss_t(k) = diss_2nd( w(k+2,j,i), w(k+1,j,i), w(k,j,i), 0., 0., & w_comp, 0.25, ddzu(k+1) ) tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) & - ( flux_d + diss_d ) ) * ddzu(k+1) ! !-- WS3 as an intermediate step k = nzb_w_inner(j,i) + 2 flux_d = flux_t(k-1) diss_d = diss_t(k-1) w_comp = w(k+1,j,i)+w(k,j,i) flux_t(k) = w_comp * ( & 7. * ( w(k+1,j,i) + w(k,j,i) ) & - ( w(k+2,j,i) + w(k-1,j,i) ) ) * adv_mom_3 diss_t(k) = - abs(w_comp) * ( & 3. * ( w(k+1,j,i) - w(k,j,i) ) & - ( w(k+2,j,i) - w(k-1,j,i) ) ) * adv_mom_3 tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) & - ( flux_d + diss_d ) ) * ddzu(k+1) ! !-- WS5 DO k = nzb_v_inner(j,i) + 3, nzt flux_d = flux_t(k-1) diss_d = diss_t(k-1) w_comp = w(k+1,j,i) + w(k,j,i) flux_t(k) = w_comp * ( & 37. * ( w(k+1,j,i) + w(k,j,i) ) & - 8. * ( w(k+2,j,i) + w(k-1,j,i) ) & + ( w(k+3,j,i) + w(k-2,j,i) ) ) * adv_mom_5 diss_t(k) = - abs(w_comp) * ( & 10. * ( w(k+1,j,i) - w(k,j,i) ) & - 5. * ( w(k+2,j,i) - w(k-1,j,i) ) & + ( w(k+3,j,i) - w(k-2,j,i) ) ) * adv_mom_5 tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) & - ( flux_d + diss_d ) ) * ddzu(k+1) ENDDO ! !-- WS3 as an intermediate step DO k = nzt - 2, nzt - 1 flux_d = flux_t(k-1) diss_d = diss_t(k-1) w_comp = w(k+1,j,i)+w(k,j,i) flux_t(k) = w_comp * ( & 7. * ( w(k+1,j,i) + w(k,j,i) ) & - ( w(k+2,j,i) + w(k-1,j,i) ) ) * adv_mom_3 diss_t(k) = - abs(w_comp) * ( & 3. * ( w(k+1,j,i) - w(k,j,i) ) & - ( w(k+2,j,i) - w(k-1,j,i) ) ) * adv_mom_3 tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) & - ( flux_d + diss_d ) ) * ddzu(k+1) ENDDO ! !-- 2nd order scheme k = nzt flux_d = flux_t(k-1) diss_d = diss_t(k-1) w_comp = w(k+1,j,i) + w(k,j,i) flux_t(k) = w_comp * ( w(k+1,j,i) + w(k,j,i) ) * 0.25 diss_t(k) = diss_2nd( w(nzt+1,j,i), w(nzt+1,j,i), w(k,j,i), & w(k-1,j,i), w(k-2,j,i), w_comp, & 0.25, ddzu(k+1) ) tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) & - ( flux_d + diss_d ) ) * ddzu(k+1) ! !-- at last vertical momentum flux is accumulated DO k = nzb_w_inner(j,i), nzt sums_ws2_ws_l(k,:) = sums_ws2_ws_l(k,:) & + ( flux_t(k) + diss_t(k) ) & * weight_substep(intermediate_timestep_count) & * rmask(j,i,:) ENDDO ENDDO ENDDO END SUBROUTINE advec_w_ws SUBROUTINE local_diss_ij( i, j, ar, var_char ) END SUBROUTINE local_diss_ij SUBROUTINE local_diss( ar ) END SUBROUTINE local_diss ! !-- Computation of 2nd order dissipation. This numerical dissipation is !-- necessary to keep a stable numerical solution in region where the !-- order of the schemes degraded. REAL FUNCTION diss_2nd( v2, v1, v0, vm1, vm2, vel_comp, factor, grid ) & RESULT(dissip) IMPLICIT NONE REAL, INTENT(IN) :: v2, v1, v0, vm1, vm2, vel_comp, factor, & grid REAL :: value_min_m, value_max_m, value_min, value_max, & value_min_p, value_max_p, diss_m, diss_0, diss_p value_min_m = MIN(v0, vm1, vm2) value_max_m = MAX(v0, vm1, vm2) value_min = MIN(v1, v0, vm2) value_max = MAX(v1, v0, vm2) value_min_p = MIN(v2, v1, v0) value_max_p = MAX(v2, v1, v0) diss_m = MAX(0.,vm1 - value_max_m) + MIN(0.,vm1 - value_min_m) diss_0 = MAX(0.,v0 - value_max) + MIN(0.,v0 - value_min) diss_p = MAX(0.,v1 - value_max_p) + MIN(0.,v1 - value_min_p) dissip = abs(vel_comp) * (diss_p - 2.*diss_0 + diss_m) & * grid**2 * factor END FUNCTION diss_2nd END MODULE advec_ws