source: palm/trunk/SOURCE/advec_ws.f90 @ 667

Last change on this file since 667 was 667, checked in by suehring, 11 years ago

summary:


Gryschka:

  • Coupling with different resolution and different numbers of PEs in ocean and atmosphere is available
  • Exchange of u and v from ocean surface to atmosphere surface
  • Mirror boundary condition for u and v at the bottom are replaced by dirichlet boundary conditions
  • Inflow turbulence is now defined by flucuations around spanwise mean
  • Bugfixes for cyclic_fill and constant_volume_flow

Suehring:

  • New advection added ( Wicker and Skamarock 5th order ), therefore:
    • New module advec_ws.f90
    • Modified exchange of ghost boundaries.
    • Modified evaluation of turbulent fluxes
    • New index bounds nxlg, nxrg, nysg, nyng

advec_ws.f90


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. Furthermore simulations within topography are
not possible so far. A further routine local_diss_ij is available and is used
if a control of dissipative fluxes is desired.

check_parameters.f90


Exchange of parameters between ocean and atmosphere via PE0
Check for illegal combination of ws-scheme and timestep scheme.
Check for topography and ws-scheme.
Check for not cyclic boundary conditions in combination with ws-scheme and
loop_optimization = 'vector'.
Check for call_psolver_at_all_substeps and ws-scheme for momentum_advec.

Different processor/grid topology in atmosphere and ocean is now allowed!
Bugfixes in checking for conserve_volume_flow_mode.

exchange_horiz.f90


Dynamic exchange of ghost points with nbgp_local to ensure that no useless
ghost points exchanged in case of multigrid. type_yz(0) and type_xz(0) used for
normal grid, the remaining types used for the several grid levels.
Exchange is done via MPI-Vectors with a dynamic value of ghost points which
depend on the advection scheme. Exchange of left and right PEs is 10% faster
with MPI-Vectors than without.

flow_statistics.f90


When advection is computed with ws-scheme, turbulent fluxes are already
computed in the respective advection routines and buffered in arrays
sums_xxxx_ws_l(). This is due to a consistent treatment of statistics
with the numerics and to avoid unphysical kinks near the surface. So some if-
requests has to be done to dicern between fluxes from ws-scheme other advection
schemes. Furthermore the computation of z_i is only done if the heat flux
exceeds a minimum value. This affects only simulations of a neutral boundary
layer and is due to reasons of computations in the advection scheme.

inflow_turbulence.f90


Using nbgp recycling planes for a better resolution of the turbulent flow near
the inflow.

init_grid.f90


Definition of new array bounds nxlg, nxrg, nysg, nyng on each PE.
Furthermore the allocation of arrays and steering of loops is done with these
parameters. Call of exchange_horiz are modified.
In case of dirichlet bounday condition at the bottom zu(0)=0.0
dzu_mg has to be set explicitly for a equally spaced grid near bottom.
ddzu_pres added to use a equally spaced grid near bottom.

init_pegrid.f90


Moved determination of target_id's from init_coupling
Determination of parameters needed for coupling (coupling_topology, ngp_a, ngp_o)
with different grid/processor-topology in ocean and atmosphere

Adaption of ngp_xy, ngp_y to a dynamic number of ghost points.
The maximum_grid_level changed from 1 to 0. 0 is the normal grid, 1 to
maximum_grid_level the grids for multigrid, in which 0 and 1 are normal grids.
This distinction is due to reasons of data exchange and performance for the
normal grid and grids in poismg.
The definition of MPI-Vectors adapted to a dynamic numer of ghost points.
New MPI-Vectors for data exchange between left and right boundaries added.
This is due to reasons of performance (10% faster).

ATTENTION: nnz_x undefined problem still has to be solved!!!!!!!!
TEST OUTPUT (TO BE REMOVED) logging mpi2 ierr values

parin.f90


Steering parameter dissipation_control added in inipar.

Makefile


Module advec_ws added.

Modules


Removed u_nzb_p1_for_vfc and v_nzb_p1_for_vfc

For coupling with different resolution in ocean and atmophere:
+nx_a, +nx_o, ny_a, +ny_o, ngp_a, ngp_o, +total_2d_o, +total_2d_a,
+coupling_topology

Buffer arrays for the left sided advective fluxes added in arrays_3d.
+flux_s_u, +flux_s_v, +flux_s_w, +diss_s_u, +diss_s_v, +diss_s_w,
+flux_s_pt, +diss_s_pt, +flux_s_e, +diss_s_e, +flux_s_q, +diss_s_q,
+flux_s_sa, +diss_s_sa
3d arrays for dissipation control added. (only necessary for vector arch.)
+var_x, +var_y, +var_z, +gamma_x, +gamma_y, +gamma_z
Default of momentum_advec and scalar_advec changed to 'ws-scheme' .
+exchange_mg added in control_parameters to steer the data exchange.
Parameters +nbgp, +nxlg, +nxrg, +nysg, +nyng added in indices.
flag array +boundary_flags added in indices to steer the degradation of order
of the advective fluxes when non-cyclic boundaries are used.
MPI-datatypes +type_y, +type_y_int and +type_yz for data_exchange added in
pegrid.
+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_wssas_ws_l, +sums_wsqs_ws_l
and +weight_substep added in statistics to steer the statistical evaluation
of turbulent fluxes in the advection routines.
LOGICALS +ws_scheme_sca and +ws_scheme_mom added to get a better performance
in prognostic_equations.
LOGICAL +dissipation_control control added to steer numerical dissipation
in ws-scheme.

Changed length of string run_description_header

pres.f90


New allocation of tend when ws-scheme and multigrid is used. This is due to
reasons of perforance of the data_exchange. The same is done with p after
poismg is called.
nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng when no
multigrid is used. Calls of exchange_horiz are modified.

bugfix: After pressure correction no volume flow correction in case of
non-cyclic boundary conditions
(has to be done only before pressure correction)

Call of SOR routine is referenced with ddzu_pres.

prognostic_equations.f90


Calls of the advection routines with WS5 added.
Calls of ws_statistics added to set the statistical arrays to zero after each
time step.

advec_particles.f90


Declaration of de_dx, de_dy, de_dz adapted to additional ghost points.
Furthermore the calls of exchange_horiz were modified.

asselin_filter.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng

average_3d_data.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng

boundary_conds.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng
Removed mirror boundary conditions for u and v at the bottom in case of
ibc_uv_b == 0. Instead, dirichelt boundary conditions (u=v=0) are set
in init_3d_model

calc_liquid_water_content.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng

calc_spectra.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng for
allocation of tend.

check_open.f90


Output of total array size was adapted to nbgp.

data_output_2d.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng in loops and
allocation of arrays local_2d and total_2d.
Calls of exchange_horiz are modified.

data_output_2d.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng in loops and
allocation of arrays. Calls of exchange_horiz are modified.
Skip-value skip_do_avs changed to a dynamic adaption of ghost points.

data_output_mask.f90


Calls of exchange_horiz are modified.

diffusion_e.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng

diffusion_s.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng

diffusion_u.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng

diffusion_v.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng

diffusion_w.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng

diffusivities.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng

diffusivities.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng.
Calls of exchange_horiz are modified.

exchange_horiz_2d.f90


Dynamic exchange of ghost points with nbgp, which depends on the advection
scheme. Exchange between left and right PEs is now done with MPI-vectors.

global_min_max.f90


Adapting of the index arrays, because MINLOC assumes lowerbound
at 1 and not at nbgp.

init_3d_model.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng in loops and
allocation of arrays. Calls of exchange_horiz are modified.
Call ws_init to initialize arrays needed for statistical evaluation and
optimization when ws-scheme is used.
Initial volume flow is now calculated by using the variable hom_sum.
Therefore the correction of initial volume flow for non-flat topography
removed (removed u_nzb_p1_for_vfc and v_nzb_p1_for_vfc)
Changed surface boundary conditions for u and v in case of ibc_uv_b == 0 from
mirror bc to dirichlet boundary conditions (u=v=0), so that k=nzb is
representative for the height z0

Bugfix: type conversion of '1' to 64bit for the MAX function (ngp_3d_inner)

init_coupling.f90


determination of target_id's moved to init_pegrid

init_pt_anomaly.f90


Call of exchange_horiz are modified.

init_rankine.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng.
Calls of exchange_horiz are modified.

init_slope.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng.

header.f90


Output of advection scheme.

poismg.f90


Calls of exchange_horiz are modified.

prandtl_fluxes.f90


Changed surface boundary conditions for u and v from mirror bc to dirichelt bc,
therefore u(uzb,:,:) and v(nzb,:,:) is now representative for the height z0
nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng

production_e.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng

read_3d_binary.f90


+/- 1 replaced with +/- nbgp when swapping and allocating variables.

sor.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng.
Call of exchange_horiz are modified.
bug removed in declaration of ddzw(), nz replaced by nzt+1

subsidence.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng.

sum_up_3d_data.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng.

surface_coupler.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng in
MPI_SEND() and MPI_RECV.
additional case for nonequivalent processor and grid topopolgy in ocean and
atmosphere added (coupling_topology = 1)

Added exchange of u and v from Ocean to Atmosphere

time_integration.f90


Calls of exchange_horiz are modified.
Adaption to slooping surface.

timestep.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng.

user_3d_data_averaging.f90, user_data_output_2d.f90, user_data_output_3d.f90,
user_actions.f90, user_init.f90, user_init_plant_canopy.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng.

user_read_restart_data.f90


Allocation with nbgp.

wall_fluxes.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng.

write_compressed.f90


Array bounds and nx, ny adapted with nbgp.

sor.f90


bug removed in declaration of ddzw(), nz replaced by nzt+1

  • Property svn:executable set to *
File size: 151.2 KB
Line 
1 MODULE advec_ws
2
3
4!-----------------------------------------------------------------------------!
5! Current revisions:
6! -----------------
7!
8!
9! Former revisions:
10! -----------------
11! $Id: advec_s_ws.f90 411 2009-12-11 12:31:43 Z suehring $
12!
13!
14! Initial revision
15!
16! Description:
17! ------------
18! Advection scheme for scalars and momentum using the flux formulation of
19! Wicker and Skamarock 5th order. Additionally the module contains of a
20! routine using for initialisation and steering of the statical evaluation.
21! The computation of turbulent fluxes takes place inside the advection
22! routines.
23! In case of vector architectures Dirichlet and Radiation boundary conditions
24! are outstanding and not available.
25! A further routine local_diss_ij is available (next weeks) and is used if a
26! control of dissipative fluxes is desired.
27! In case of vertical grid stretching the order of advection scheme is
28! degraded. This is also realized for the ocean where the stretching is
29! downwards.
30!
31! OUTSTANDING: - Dirichlet and Radiation boundary conditions for
32!                vector architectures
33!              - dissipation control for cache architectures ( next weeks )
34!              - Topography ( next weeks )
35!-----------------------------------------------------------------------------!
36
37
38    PRIVATE
39    PUBLIC   advec_s_ws, advec_u_ws, advec_v_ws, advec_w_ws, &
40             local_diss, ws_init, ws_statistics
41
42    INTERFACE ws_init
43       MODULE PROCEDURE ws_init
44    END INTERFACE ws_init
45    INTERFACE ws_statistics
46       MODULE PROCEDURE ws_statistics
47    END INTERFACE ws_statistics
48    INTERFACE advec_s_ws
49       MODULE PROCEDURE advec_s_ws
50       MODULE PROCEDURE advec_s_ws_ij
51    END INTERFACE advec_s_ws
52    INTERFACE advec_u_ws
53       MODULE PROCEDURE advec_u_ws
54       MODULE PROCEDURE advec_u_ws_ij
55    END INTERFACE advec_u_ws
56    INTERFACE advec_v_ws
57       MODULE PROCEDURE advec_v_ws
58       MODULE PROCEDURE advec_v_ws_ij
59    END INTERFACE advec_v_ws
60    INTERFACE advec_w_ws
61       MODULE PROCEDURE advec_w_ws
62       MODULE PROCEDURE advec_w_ws_ij
63    END INTERFACE advec_w_ws
64    INTERFACE local_diss
65       MODULE PROCEDURE local_diss
66       MODULE PROCEDURE local_diss_ij
67    END INTERFACE local_diss
68
69 CONTAINS
70
71!
72!-- Initialisation
73
74    SUBROUTINE ws_init
75       USE arrays_3d
76       USE constants
77       USE control_parameters
78       USE indices
79       USE statistics
80
81!       
82!--    Set the LOGICALS to enhance the performance.
83       IF ( momentum_advec == 'ws-scheme' )  ws_scheme_mom = .TRUE.
84       IF ( scalar_advec == 'ws-scheme' )  ws_scheme_sca = .TRUE.
85!
86!--    Allocate arrays needed for dissipation control.
87       IF ( dissipation_control )  THEN
88   !          ALLOCATE(var_x(nzb+1:nzt,nys:nyn,nxl:nxr),  &
89    !                  var_y(nzb+1:nzt,nys:nyn,nxl:nxr),  &
90     !                 var_z(nzb+1:nzt,nys:nyn,nxl:nxr),  &
91      !                gamma_x(nzb:nzt+1,nysg:nyng,nxlg:nxrg),  &
92       !               gamma_y(nzb:nzt+1,nysg:nyng,nxlg:nxrg),  &
93        !              gamma_z(nzb:nzt+1,nysg:nyng,nxlg:nxrg))
94       ENDIF
95             
96!--    Set the appropriate factors for scalar and momentum advection.
97       adv_sca_5 = 0.016666666666666
98       adv_sca_3 = 0.083333333333333 
99       adv_mom_5 = 0.0083333333333333
100       adv_mom_3 = 0.041666666666666
101!
102!--    Allocate arrays needed for statistical evaluation of fluxes when
103!--    ws-scheme is used.
104!--    The following array contains the weighting factors of the RK3 average
105!--    of tendecies.
106       ALLOCATE (weight_substep(1:intermediate_timestep_count_max) )
107       IF ( intermediate_timestep_count_max == 3)  THEN   !RK3
108          weight_substep(1) = 0.166666666666666
109          weight_substep(2) = 0.3
110          weight_substep(3) = 0.533333333333333
111       ENDIF
112!         
113!--    Arrays needed for statical evaluation of fluxes.
114       IF ( ws_scheme_mom )  THEN
115          ALLOCATE( sums_wsus_ws_l(nzb:nzt+1,0:statistic_regions),         &
116              sums_wsvs_ws_l(nzb:nzt+1,0:statistic_regions),               &
117              sums_us2_ws_l(nzb:nzt+1,0:statistic_regions),                &
118              sums_vs2_ws_l(nzb:nzt+1,0:statistic_regions),                &
119              sums_ws2_ws_l(nzb:nzt+1,0:statistic_regions))
120
121          sums_wsus_ws_l = 0.0
122          sums_wsvs_ws_l = 0.0
123          sums_us2_ws_l  = 0.0
124          sums_vs2_ws_l  = 0.0
125          sums_ws2_ws_l  = 0.0
126       ENDIF
127
128       IF ( ws_scheme_sca )  THEN
129          ALLOCATE( sums_wspts_ws_l(nzb:nzt+1,0:statistic_regions))
130          sums_wspts_ws_l = 0.0
131          IF ( humidity  .OR.  passive_scalar )  THEN
132             ALLOCATE( sums_wsqs_ws_l(nzb:nzt+1,0:statistic_regions))
133             sums_wsqs_ws_l = 0.0
134          ENDIF
135          IF ( ocean )  THEN
136             ALLOCATE( sums_wssas_ws_l(nzb:nzt+1,0:statistic_regions))
137             sums_wssas_ws_l = 0.0
138          ENDIF
139       ENDIF
140!
141!--    Arrays needed for reasons of speed optimization for cache and noopt
142!--    version. For the vector version the buffer arrays are not necessary,
143!--    because the the fluxes can swapped directly inside the loops of the
144!--    advection routines.
145       IF ( loop_optimization /= 'vector' )  THEN
146          IF ( ws_scheme_mom )  THEN
147             ALLOCATE(flux_s_u(nzb+1:nzt), flux_s_v(nzb+1:nzt),            &
148                flux_s_w(nzb+1:nzt), diss_s_u(nzb+1:nzt),                  &
149                diss_s_v(nzb+1:nzt), diss_s_w(nzb+1:nzt))
150             ALLOCATE(flux_l_u(nzb+1:nzt,nys:nyn),                         &
151                flux_l_v(nzb+1:nzt,nys:nyn), flux_l_w(nzb+1:nzt,nys:nyn),  & 
152                diss_l_u(nzb+1:nzt,nys:nyn), diss_l_v(nzb+1:nzt,nys:nyn),  & 
153                diss_l_w(nzb+1:nzt,nys:nyn)) 
154          ENDIF
155          IF ( ws_scheme_sca )  THEN
156             ALLOCATE(flux_s_pt(nzb+1:nzt), flux_s_e(nzb+1:nzt),           &
157                diss_s_pt(nzb+1:nzt), diss_s_e(nzb+1:nzt))
158             ALLOCATE(flux_l_pt(nzb+1:nzt,nys:nyn),                        &
159                flux_l_e(nzb+1:nzt,nys:nyn), diss_l_pt(nzb+1:nzt,nys:nyn), & 
160                diss_l_e(nzb+1:nzt,nys:nyn))
161             IF ( humidity  .OR.  passive_scalar )  THEN
162                ALLOCATE(flux_s_q(nzb+1:nzt), diss_s_q(nzb+1:nzt))
163                ALLOCATE(flux_l_q(nzb+1:nzt,nys:nyn),                      &
164                   diss_l_q(nzb+1:nzt,nys:nyn))
165             ENDIF
166             IF ( ocean )  THEN
167                ALLOCATE(flux_s_sa(nzb+1:nzt), diss_s_sa(nzb+1:nzt))
168                ALLOCATE(flux_l_sa(nzb+1:nzt,nys:nyn),                     &
169                   diss_l_sa(nzb+1:nzt,nys:nyn))
170             ENDIF
171          ENDIF
172       ENDIF
173!
174!--    Determine the flags where the order of the scheme for horizontal
175!--    advection should be degraded.
176       ALLOCATE( boundary_flags(nys:nyn,nxl:nxr) )
177       DO  i = nxl, nxr
178          DO  j = nys, nyn
179             boundary_flags(j,i) = 0
180             IF ( outflow_l )  THEN
181                IF ( i == nxlu )  boundary_flags(j,i) = 5
182                IF ( i == nxl )  boundary_flags(j,i) = 6
183             ELSEIF ( outflow_r )  THEN
184                IF ( i == nxr-1 )  boundary_flags(j,i) = 1
185                IF ( i == nxr )  boundary_flags(j,i) = 2
186             ELSEIF ( outflow_n )  THEN
187                IF ( j == nyn-1 )  boundary_flags(j,i) = 3
188                IF ( j == nyn )  boundary_flags(j,i) = 4
189             ELSEIF ( outflow_s )  THEN
190                IF ( j == nysv )  boundary_flags(j,i) = 7
191                IF ( j == nys )  boundary_flags(j,i) = 8
192             ENDIF
193          ENDDO
194       ENDDO
195       
196    END SUBROUTINE ws_init
197   
198    SUBROUTINE ws_statistics
199   
200       USE control_parameters
201       USE statistics
202
203!       
204!--      The arrays needed for statistical evaluation are set to to 0 at the
205!--      begin of prognostic_equations.
206          IF ( ws_scheme_mom )  THEN
207             sums_wsus_ws_l = 0.0
208             sums_wsvs_ws_l = 0.0
209             sums_us2_ws_l = 0.0
210             sums_vs2_ws_l = 0.0
211             sums_ws2_ws_l = 0.0
212          ENDIF
213          IF ( ws_scheme_sca )  THEN
214             sums_wspts_ws_l = 0.0
215             IF ( humidity  .OR.  passive_scalar )  sums_wsqs_ws_l = 0.0
216             IF ( ocean )  sums_wssas_ws_l = 0.0
217          ENDIF
218
219    END SUBROUTINE ws_statistics
220
221
222
223!------------------------------------------------------------------------------!
224! Scalar advection - Call for grid point i,j
225!------------------------------------------------------------------------------!
226    SUBROUTINE advec_s_ws_ij( i, j, sk, sk_char,swap_flux_y_local,  &
227               swap_diss_y_local, swap_flux_x_local, swap_diss_x_local  )
228
229       USE arrays_3d
230       USE constants
231       USE control_parameters
232       USE grid_variables
233       USE indices
234       USE statistics
235
236       IMPLICIT NONE
237
238       INTEGER ::  i, j, k
239       LOGICAL :: degraded_l = .FALSE., degraded_s = .FALSE.
240       REAL    ::  flux_d, diss_d, u_comp, v_comp
241       REAL, DIMENSION(:,:,:), POINTER    :: sk
242       REAL, DIMENSION(nzb:nzt+1)         :: flux_t, diss_t, flux_r, diss_r,  &
243                                             flux_n, diss_n
244       REAL, DIMENSION(nzb+1:nzt)         :: swap_flux_y_local,               &       
245                                             swap_diss_y_local
246       REAL, DIMENSION(nzb+1:nzt,nys:nyn) :: swap_flux_x_local,               &
247                                             swap_diss_x_local
248       CHARACTER (LEN = *), INTENT(IN)    :: sk_char
249
250       IF ( boundary_flags(j,i) /= 0 )  THEN
251!       
252!--      Degrade the order for Dirichlet bc. at the outflow boundary.
253         SELECT CASE ( boundary_flags(j,i) )
254         
255            CASE ( 1 )
256               DO  k = nzb_s_inner(j,i) + 1, nzt
257                  u_comp    = u(k,j,i+1) - u_gtrans
258                  flux_r(k) = u_comp * (                                      &
259                              7. * ( sk(k,j,i+1) + sk(k,j,i)   )              &
260                              -    ( sk(k,j,i+2) + sk(k,j,i-1) ) ) * adv_sca_3
261                  diss_r(k) = - abs(u_comp) * (                               &
262                              3. * ( sk(k,j,i+1)-sk(k,j,i)     )              & 
263                              -    ( sk(k,j,i+2) - sk(k,j,i-1) ) ) * adv_sca_3
264               ENDDO
265               
266            CASE ( 2 )
267               DO  k = nzb_s_inner(j,i) + 1, nzt
268                  u_comp    = u(k,j,i+1) - u_gtrans
269                  flux_r(k) = u_comp * ( sk(k,j,i+1) + sk(k,j,i) ) * 0.5
270                  diss_r(k) = diss_2nd( sk(k,j,i+1), sk(k,j,i+1), sk(k,j,i),  &
271                                        sk(k,j,i-1), sk(k,j,i-2), u_comp,     &
272                                        0.5, ddx )
273              ENDDO
274             
275            CASE ( 3 )
276               DO  k = nzb_s_inner(j,i) + 1, nzt
277                  v_comp = v(k,j+1,i) - v_gtrans
278                  flux_n(k) = v_comp * (                                      &
279                              7. * ( sk(k,j+1,i) + sk(k,j,i)   )              &
280                              -    ( sk(k,j+2,i) + sk(k,j-1,i) ) ) * adv_sca_3
281                  diss_n(k) = - abs(v_comp) * (                               &
282                              3. * ( sk(k,j+1,i)-sk(k,j,i)    )               & 
283                              -    (sk(k,j+2,i) - sk(k,j-1,i) ) ) * adv_sca_3
284              ENDDO
285            CASE ( 4 )
286               DO  k = nzb_s_inner(j,i) + 1, nzt
287                  v_comp    = v(k,j+1,i) - v_gtrans
288                  flux_n(k) = v_comp* ( sk(k,j+1,i) + sk(k,j,i) ) * 0.5
289                  diss_n(k) =  diss_2nd( sk(k,j+1,i), sk(k,j+1,i), sk(k,j,i), &
290                                         sk(k,j-1,i), sk(k,j-2,i), v_comp,    &
291                                         0.5, ddy )     
292               ENDDO
293               
294            CASE ( 5 )
295               DO  k = nzb_w_inner(j,i) + 1, nzt
296                  u_comp    = u(k,j,i+1) - u_gtrans
297                  flux_r(k) = u_comp  * (                                     &
298                              7. * ( sk(k,j,i+1) + sk(k,j,i)   )              &
299                              -    ( sk(k,j,i+2) + sk(k,j,i-1) ) ) * adv_sca_3
300                  diss_r(k) = - abs(u_comp) * (                               &
301                              3. * ( sk(k,j,i+1)-sk(k,j,i) )                  & 
302                              -    ( sk(k,j,i+2) - sk(k,j,i-1) ) ) * adv_sca_3
303               ENDDO 
304               
305            CASE ( 6 )
306               DO  k = nzb_s_inner(j,i) + 1, nzt
307                  u_comp    = u(k,j,i+1) - u_gtrans
308                  flux_r(k) = u_comp * (                                      &
309                              7. * ( sk(k,j,i+1) + sk(k,j,i)   )              &
310                              -    ( sk(k,j,i+2) + sk(k,j,i-1) ) ) * adv_sca_3
311                  diss_r(k) = - abs(u_comp) * (                               &
312                              3. * ( sk(k,j,i+1)-sk(k,j,i)     )              & 
313                              -    ( sk(k,j,i+2) - sk(k,j,i-1) ) ) * adv_sca_3
314!                             
315!--               Compute leftside fluxes for the left boundary of PE domain.
316                  u_comp                 = u(k,j,i) - u_gtrans
317                  swap_flux_x_local(k,j) = u_comp * (                         &
318                                           sk(k,j,i) + sk(k,j,i-1) ) * 0.5
319                  swap_diss_x_local(k,j) = diss_2nd( sk(k,j,i+2),sk(k,j,i+1), &
320                                                     sk(k,j,i), sk(k,j,i-1),  &
321                                                     sk(k,j,i-1), u_comp,     &
322                                                     0.5, ddx ) 
323               ENDDO
324               degraded_l = .TRUE.
325               
326            CASE ( 7 )
327               DO  k = nzb_s_inner(j,i) + 1, nzt
328                  v_comp    = v(k,j+1,i)-v_gtrans
329                  flux_n(k) = v_comp * (                                      &
330                              7. * ( sk(k,j+1,i) + sk(k,j,i)   )              &
331                              -    ( sk(k,j+2,i) + sk(k,j-1,i) ) ) * adv_sca_3
332                  diss_n(k) = - abs(v_comp) * (                               &
333                              3. * ( sk(k,j+1,i) - sk(k,j,i)   )              & 
334                              -    ( sk(k,j+2,i) - sk(k,j-1,i) ) ) * adv_sca_3
335               ENDDO
336               
337            CASE ( 8 )
338               DO  k = nzb_s_inner(j,i) + 1, nzt
339                  v_comp    = v(k,j+1,i) - v_gtrans
340                  flux_n(k) = v_comp * (                                      &
341                              7. * ( sk(k,j+1,i) + sk(k,j,i)   )              &
342                              -    ( sk(k,j+2,i) + sk(k,j-1,i) ) ) * adv_sca_3
343                  diss_n(k) = - abs(v_comp) * (                               &
344                              3. * ( sk(k,j+1,i) - sk(k,j,i) )                & 
345                              -    ( sk(k,j+2,i) - sk(k,j-1,i) ) ) * adv_sca_3
346!                             
347!--               Compute southside fluxes for the south boundary of PE domain           
348                  v_comp               = v(k,j,i) - v_gtrans
349                  swap_flux_y_local(k) = v_comp * (                           &
350                                         sk(k,j,i)+ sk(k,j-1,i) ) * 0.5 
351                  swap_diss_y_local(k) = diss_2nd( sk(k,j+2,i), sk(k,j+1,i),  &
352                                                   sk(k,j,i), sk(k,j-1,i),    &
353                                                   sk(k,j-1,i), v_comp,       &
354                                                   0.5, ddy )
355               ENDDO
356               degraded_s = .TRUE.
357               
358            CASE DEFAULT
359           
360         END SELECT
361!         
362!--      Compute the crosswise 5th order fluxes at the outflow
363         IF ( boundary_flags(j,i) == 1 .OR. boundary_flags(j,i) == 2          &
364         .OR. boundary_flags(j,i) == 5 .OR. boundary_flags(j,i) == 6 )  THEN
365            DO  k = nzb_s_inner(j,i) + 1, nzt
366               v_comp    = v(k,j+1,i) - v_gtrans
367               flux_n(k) = v_comp * (                                         &
368                           37.  * ( sk(k,j+1,i) + sk(k,j,i)   )               &
369                           - 8. * ( sk(k,j+2,i) + sk(k,j-1,i) )               &
370                           +      ( sk(k,j+3,i) + sk(k,j-2,i) ) ) * adv_sca_5
371               diss_n(k) = - abs(v_comp) * (                                  &
372                           10. *  ( sk(k,j+1,i) - sk(k,j,i)   )               &
373                           - 5. * ( sk(k,j+2,i) - sk(k,j-1,i) )               &
374                           +      ( sk(k,j+3,i) - sk(k,j-2,i) ) ) * adv_sca_5
375            ENDDO
376           
377         ELSE
378         
379            DO  k = nzb_s_inner(j,i) + 1, nzt
380               u_comp    = u(k,j,i+1) - u_gtrans 
381               flux_r(k) = u_comp * (                                         &
382                           37.  * ( sk(k,j,i+1) + sk(k,j,i)   )               &
383                           - 8. * ( sk(k,j,i+2) + sk(k,j,i-1) )               &
384                           +      ( sk(k,j,i+3) + sk(k,j,i-2) ) ) * adv_sca_5
385               diss_r(k) = - abs(u_comp) * (                                  &
386                           10.  * ( sk(k,j,i+1) - sk(k,j,i)   )               &
387                           - 5. * ( sk(k,j,i+2) - sk(k,j,i-1) )               &
388                           +      ( sk(k,j,i+3) - sk(k,j,i-2) ) ) * adv_sca_5
389            ENDDO
390           
391         ENDIF
392!
393!--    Compute the fifth order fluxes for the interior of PE domain.
394       ELSE
395               
396          DO  k = nzb_u_inner(j,i) + 1, nzt
397             u_comp    = u(k,j,i+1) - u_gtrans
398             flux_r(k) = u_comp * (                                           &
399                         37.  * ( sk(k,j,i+1) + sk(k,j,i)   )                 &
400                         - 8. * ( sk(k,j,i+2) + sk(k,j,i-1) )                 &
401                         +      ( sk(k,j,i+3) + sk(k,j,i-2) ) ) * adv_sca_5
402             diss_r(k) = - abs(u_comp) * (                                    &
403                         10.  * ( sk(k,j,i+1) - sk(k,j,i)   )                 &
404                         - 5. * ( sk(k,j,i+2) - sk(k,j,i-1) )                 &
405                         +      ( sk(k,j,i+3) - sk(k,j,i-2) ) ) * adv_sca_5
406
407             v_comp    = v(k,j+1,i) - v_gtrans
408             flux_n(k) = v_comp * (                                           &
409                         37.  * ( sk(k,j+1,i) + sk(k,j,i)   )                 &
410                         - 8. * ( sk(k,j+2,i) + sk(k,j-1,i) )                 &
411                         +      ( sk(k,j+3,i) + sk(k,j-2,i) ) ) * adv_sca_5
412             diss_n(k) = - abs(v_comp) * (                                    &
413                         10.  * ( sk(k,j+1,i) - sk(k,j,i)   )                 &                     
414                         - 5. * ( sk(k,j+2,i) - sk(k,j-1,i) )                 &
415                         +      ( sk(k,j+3,i) - sk(k,j-2,i) ) ) * adv_sca_5               
416          ENDDO
417         
418       ENDIF
419!       
420!--    Compute left- and southside fluxes of the respective PE bounds.       
421       IF ( j == nys .AND. ( .NOT. degraded_s ) )  THEN
422       
423           DO  k = nzb_s_inner(j,i) + 1, nzt
424              v_comp               = v(k,j,i) - v_gtrans
425              swap_flux_y_local(k) = v_comp * (                               &
426                                     37. *    ( sk(k,j,i) + sk(k,j-1,i)   )   &
427                                     - 8. *   ( sk(k,j+1,i) + sk(k,j-2,i) )   &
428                                     +        ( sk(k,j+2,i) + sk(k,j-3,i) )   &
429                                     ) * adv_sca_5
430              swap_diss_y_local(k) = - abs(v_comp) * (                        &
431                                     10. *  ( sk(k,j,i) - sk(k,j-1,i)   )     &
432                                     - 5. * ( sk(k,j+1,i) - sk(k,j-2,i) )     &
433                                     +        sk(k,j+2,i) - sk(k,j-3,i)       &
434                                     ) * adv_sca_5
435           ENDDO           
436         
437        ENDIF
438       
439        IF ( i == nxl .AND. ( .NOT. degraded_l ) )  THEN
440       
441           DO  k = nzb_s_inner(j,i) + 1, nzt
442              u_comp                 = u(k,j,i) - u_gtrans
443              swap_flux_x_local(k,j) = u_comp * (                             &
444                                       37.  * ( sk(k,j,i) + sk(k,j,i-1)   )   &
445                                       - 8. * ( sk(k,j,i+1) + sk(k,j,i-2) )   &
446                                       +      ( sk(k,j,i+2) + sk(k,j,i-3) )   &
447                                       ) * adv_sca_5
448              swap_diss_x_local(k,j) = - abs(u_comp) * (                      &
449                                       10.  * ( sk(k,j,i) - sk(k,j,i-1)   )   &
450                                       - 5. * ( sk(k,j,i+1) - sk(k,j,i-2) )   &
451                                       +      ( sk(k,j,i+2) - sk(k,j,i-3) )   &
452                                       ) * adv_sca_5
453           ENDDO
454           
455        ENDIF
456!       
457!--    Now compute the tendency terms for the horizontal parts.
458       DO  k = nzb_s_inner(j,i) + 1, nzt
459         
460          tend(k,j,i) = tend(k,j,i) - (                                       &
461               ( flux_r(k) + diss_r(k)                                        &
462             - ( swap_flux_x_local(k,j) + swap_diss_x_local(k,j) ) ) * ddx    &
463             + (flux_n(k)+diss_n(k)                                           &
464             - (swap_flux_y_local(k)    + swap_diss_y_local(k) ) )   * ddy )
465
466          swap_flux_y_local(k)   = flux_n(k)
467          swap_diss_y_local(k)   = diss_n(k)
468          swap_flux_x_local(k,j) = flux_r(k)
469          swap_diss_x_local(k,j) = diss_r(k)
470         
471       ENDDO
472!
473!--    Vertical advection, degradation of order near surface and top.
474!--    The fluxes flux_d and diss_d at the surface are 0. Due to reasons of
475!--    statistical evaluation the top flux at the surface should be 0
476
477       flux_t(nzb_s_inner(j,i)) = 0.
478       diss_t(nzb_s_inner(j,i)) = 0.
479!       
480!--    2nd order scheme       
481       k         = nzb_s_inner(j,i) + 1
482       flux_d    = flux_t(k-1)
483       diss_d    = diss_t(k-1)
484       flux_t(k) = w(k,j,i) * ( sk(k+1,j,i) + sk(k,j,i) ) * 0.5
485!       
486!--    sk(k,j,i) is referenced three times to avoid a access below surface.
487       diss_t(k) = diss_2nd( sk(k+2,j,i), sk(k+1,j,i), sk(k,j,i), sk(k,j,i),  &
488                             sk(k,j,i), w(k,j,i), 0.5, ddzw(k) )
489                   
490       tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)                    &
491                                 - ( flux_d + diss_d )     ) * ddzw(k)
492!
493!--    WS3 as an intermediate step       
494       k         = nzb_s_inner(j,i) + 2
495       flux_d    = flux_t(k-1)
496       diss_d    = diss_t(k-1)
497       flux_t(k) = w(k,j,i) * (                                               &
498                   7. * ( sk(k+1,j,i) + sk(k,j,i)   )                         &
499                   -    ( sk(k+2,j,i) + sk(k-1,j,i) ) ) * adv_sca_3
500       diss_t(k) = - abs(w(k,j,i)) * (                                        &
501                   3. * ( sk(k+1,j,i) - sk(k,j,i)   )                         & 
502                   -    ( sk(k+2,j,i) - sk(k-1,j,i) ) ) * adv_sca_3
503
504       tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)                    &
505                                 - ( flux_d + diss_d ) ) * ddzw(k)
506!
507!--    WS5
508       DO  k = nzb_s_inner(j,i) + 3, nzt - 2
509       
510          flux_d    = flux_t(k-1)
511          diss_d    = diss_t(k-1)
512         
513          flux_t(k) = w(k,j,i) * (                                            & 
514                      37.  * ( sk(k+1,j,i) + sk(k,j,i)   )                    &
515                      - 8. * ( sk(k+2,j,i) + sk(k-1,j,i) )                    &
516                      +      ( sk(k+3,j,i) + sk(k-2,j,i) ) ) *adv_sca_5
517          diss_t(k) = - abs(w(k,j,i)) * (                                     & 
518                      10.  * ( sk(k+1,j,i) - sk(k,j,i)   )                    &     
519                      - 5. * ( sk(k+2,j,i) - sk(k-1,j,i) )                    &
520                      +      ( sk(k+3,j,i) - sk(k-2,j,i) ) ) * adv_sca_5
521
522          tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)                 &
523                                    - ( flux_d + diss_d ) ) * ddzw(k)
524       ENDDO
525!
526!--    WS3 as an intermediate step
527       k         = nzt - 1
528       flux_d    = flux_t(k-1)
529       diss_d    = diss_t(k-1)
530       flux_t(k) = w(k,j,i) * (                                               &
531                   7. * ( sk(k+1,j,i) + sk(k,j,i) )                           &
532                   -    ( sk(k+2,j,i) + sk(k-1,j,i) ) ) * adv_sca_3
533       diss_t(k) = - abs(w(k,j,i)) * (                                        &
534                   3. * ( sk(k+1,j,i) - sk(k,j,i) )                           & 
535                   -    ( sk(k+2,j,i) - sk(k-1,j,i) ) ) * adv_sca_3
536
537       tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)                    &
538                                 - ( flux_d + diss_d ) ) * ddzw(k)
539!                                 
540!--    2nd
541       k         = nzt
542       flux_d    = flux_t(k-1)
543       diss_d    = diss_t(k-1)
544       flux_t(k) = w(k,j,i) * ( sk(k+1,j,i) + sk(k,j,i) ) *0.5
545!       
546!--    sk(k+1) is referenced two times to avoid a segmentation fault at top.
547       diss_t(k) = diss_2nd( sk(k+1,j,i), sk(k+1,j,i), sk(k,j,i), sk(k-1,j,i),&
548                             sk(k-2,j,i), w(k,j,i), 0.5, ddzw(k) )
549
550       tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)                    &
551                                 - ( flux_d + diss_d ) ) * ddzw(k)
552!       
553!--    evaluation of statistics
554       SELECT CASE ( sk_char )
555
556          CASE ( 'pt' )
557            DO  k = nzb_s_inner(j,i), nzt
558               sums_wspts_ws_l(k,:) = sums_wspts_ws_l(k,:)                     &
559                 + (flux_t(k)+diss_t(k))                                       &
560                 * weight_substep(intermediate_timestep_count) * rmask(j,i,:)
561             ENDDO
562             
563          CASE ( 'sa' )
564            DO  k = nzb_s_inner(j,i), nzt
565               sums_wssas_ws_l(k,:) = sums_wssas_ws_l(k,:)                     &
566                 +(flux_t(k)+diss_t(k))                                        &
567                 * weight_substep(intermediate_timestep_count) * rmask(j,i,:)
568             ENDDO
569             
570          CASE ( 'q' )
571             DO  k = nzb_s_inner(j,i), nzt
572                sums_wsqs_ws_l(k,:) = sums_wsqs_ws_l(k,:)                      &
573                  + (flux_t(k)+diss_t(k))                                      &
574                  * weight_substep(intermediate_timestep_count) * rmask(j,i,:)
575             ENDDO
576
577         END SELECT
578
579    END SUBROUTINE advec_s_ws_ij
580
581
582
583
584!------------------------------------------------------------------------------!
585! Advection of u-component - Call for grid point i,j
586!------------------------------------------------------------------------------!
587    SUBROUTINE advec_u_ws_ij( i, j )
588
589       USE arrays_3d
590       USE constants
591       USE control_parameters
592       USE grid_variables
593       USE indices
594       USE statistics
595
596       IMPLICIT NONE
597
598       INTEGER ::  i, j, k
599       LOGICAL ::  degraded_l = .FALSE., degraded_s = .FALSE.
600       REAL    ::  gu, gv, flux_d, diss_d, u_comp_l, v_comp, w_comp
601       REAL, DIMENSION(nzb:nzt+1)  ::  flux_t, diss_t, flux_r, diss_r,        &
602                                       flux_n, diss_n, u_comp
603
604       gu = 2.0 * u_gtrans
605       gv = 2.0 * v_gtrans
606         
607       IF ( boundary_flags(j,i) /= 0 )  THEN
608!       
609!--      Degrade the order for Dirichlet bc. at the outflow boundary
610         SELECT CASE ( boundary_flags(j,i) )
611         
612            CASE ( 1 )
613               DO  k = nzb_u_inner(j,i) + 1, nzt
614                  u_comp(k) = u(k,j,i+1) + u(k,j,i)
615                  flux_r(k) = ( u_comp(k) - gu ) * (                          &
616                              7. * (u(k,j,i+1) + u(k,j,i)    )                &
617                              -    ( u(k,j,i+2) + u(k,j,i-1) ) ) * adv_mom_3
618                  diss_r(k) = - abs(u_comp(k) - gu) * (                       &
619                              3. * (u(k,j,i+1) - u(k,j,i)    )                & 
620                              -    ( u(k,j,i+2) - u(k,j,i-1) ) ) * adv_mom_3
621              ENDDO
622             
623            CASE ( 2 )
624               DO  k = nzb_u_inner(j,i) + 1, nzt
625                  u_comp(k) = u(k,j,i+1) + u(k,j,i)
626                  flux_r(k) = (u_comp(k) - gu) * ( u(k,j,i+1) + u(k,j,i) )*0.25 
627                  diss_r(k) = diss_2nd( u(k,j,i+1) ,u(k,j,i+1), u(k,j,i),     &
628                                        u(k,j,i-1), u(k,j,i-2), u_comp(k),    &
629                                        0.25, ddx )
630               ENDDO
631               
632            CASE ( 3 )
633               DO  k = nzb_u_inner(j,i) + 1, nzt
634                  v_comp    = v(k,j+1,i) + v(k,j+1,i-1) - gv
635                  flux_n(k) = v_comp * (                                      &
636                              7. * ( u(k,j+1,i) + u(k,j,i) )                  &
637                              -    ( u(k,j+2,i) + u(k,j-1,i) ) ) * adv_mom_3
638                  diss_n(k) = - abs(v_comp) * (                               &
639                              3. * ( u(k,j+1,i) - u(k,j,i)   )                & 
640                              -    ( u(k,j+2,i) - u(k,j-1,i) ) ) * adv_mom_3
641               ENDDO
642               
643            CASE ( 4 )
644               DO  k = nzb_u_inner(j,i) + 1, nzt
645                  v_comp    = v(k,j+1,i) + v(k,j+1,i-1) - gv
646                  flux_n(k) = v_comp * ( u(k,j+1,i) + u(k,j,i) ) * 0.25
647                  diss_n(k) = diss_2nd( u(k,j+1,i), u(k,j+1,i), u(k,j,i),     &
648                                        u(k,j-1,i), u(k,j-2,i), v_comp,       &
649                                        0.25, ddy )
650               ENDDO
651               
652            CASE ( 5 )
653               DO  k = nzb_u_inner(j,i) + 1, nzt
654!               
655!--               Compute leftside fluxes for the left boundary of PE domain
656                  u_comp(k)     = u(k,j,i) + u(k,j,i-1) - gu
657                  flux_l_u(k,j) = u_comp(k) * ( u(k,j,i) + u(k,j,i-1) ) * 0.25
658                  diss_l_u(k,j) = diss_2nd(u(k,j,i+2), u(k,j,i+1), u(k,j,i),  &
659                                           u(k,j,i-1), u(k,j,i-1), u_comp(k), &
660                                           0.25, ddx )
661                 
662                  u_comp(k) = u(k,j,i+1) + u(k,j,i)
663                  flux_r(k) = ( u_comp(k) - gu ) * (                          &
664                              7. * (u(k,j,i+1) + u(k,j,i)    )                &
665                              -    ( u(k,j,i+2) + u(k,j,i-1) ) ) * adv_mom_3
666                  diss_r(k) = - abs( u_comp(k) -gu ) * (                      &
667                              3. * ( u(k,j,i+1) - u(k,j,i)     )              & 
668                              -    ( u(k,j,i+2) - u(k,j,i-1) ) ) * adv_mom_3
669               ENDDO
670               degraded_l = .TRUE.
671               
672            CASE ( 7 )
673               DO  k = nzb_u_inner(j,i) + 1, nzt                           
674                  v_comp    = v(k,j+1,i) + v(k,j+1,i-1) - gv
675                  flux_n(k) = v_comp * (                                      &
676                              7. * ( u(k,j+1,i) + u(k,j,i)   )                &
677                              -    ( u(k,j+2,i) + u(k,j-1,i) ) ) * adv_mom_3
678                  diss_n(k) = - abs(v_comp) * (                               &
679                              3. * ( u(k,j+1,i) - u(k,j,i)   )                & 
680                              -    ( u(k,j+2,i) - u(k,j-1,i) ) ) * adv_mom_3
681               ENDDO
682               
683            CASE ( 8 )
684               DO  k = nzb_u_inner(j,i) + 1, nzt
685!               
686!--              Compute southside fluxes for the south boundary of PE domain           
687                  v_comp      = v(k,j,i) + v(k,j,i-1) - gv
688                  flux_s_u(k) = v_comp * ( u(k,j,i) + u(k,j-1,i) ) * 0.25 
689                  diss_s_u(k) = diss_2nd( u(k,j+2,i), u(k,j+1,i), u(k,j,i),   &
690                                          u(k,j-1,i), u(k,j-1,i), v_comp,     &
691                                          0.25, ddy )
692                                 
693                  v_comp    = v(k,j+1,i) + v(k,j+1,i-1) - gv
694                  flux_n(k) = v_comp * (                                      &
695                              7. * ( u(k,j+1,i) + u(k,j,i)   )                &
696                              -    ( u(k,j+2,i) + u(k,j-1,i) ) ) * adv_mom_3
697                  diss_n(k) = - abs(v_comp) * (                               &
698                              3. * ( u(k,j+1,i) - u(k,j,i)   )                & 
699                              -    ( u(k,j+2,i) - u(k,j-1,i) ) ) * adv_mom_3
700               ENDDO 
701               degraded_s = .TRUE.
702               
703            CASE DEFAULT
704           
705         END SELECT
706!         
707!--      Compute the crosswise 5th order fluxes at the outflow
708         IF ( boundary_flags(j,i) == 1 .OR. boundary_flags(j,i) == 2          &
709         .OR. boundary_flags(j,i) == 5 )  THEN
710         
711             DO  k = nzb_u_inner(j,i)+1, nzt
712               v_comp    = v(k,j+1,i) + v(k,j+1,i-1) - gv
713               flux_n(k) = v_comp * (                                         &
714                           37.  * ( u(k,j+1,i) + u(k,j,i)   )                 &
715                           - 8. * ( u(k,j+2,i) + u(k,j-1,i) )                 &
716                           +      ( u(k,j+3,i) + u(k,j-2,i) ) ) *adv_mom_5
717               diss_n(k) = - abs(v_comp) * (                                  &
718                           10.  * ( u(k,j+1,i) - u(k,j,i)   )                 &
719                           - 5. * ( u(k,j+2,i) - u(k,j-1,i) )                 &
720                           +      ( u(k,j+3,i) - u(k,j-2,i) ) ) * adv_mom_5
721             ENDDO
722             
723         ELSE
724         
725            DO  k = nzb_u_inner(j,i) + 1, nzt
726               u_comp(k) = u(k,j,i+1) + u(k,j,i)
727               flux_r(k) = ( u_comp(k) - gu ) * (                             &
728                           37.  * ( u(k,j,i+1) + u(k,j,i)   )                 &
729                           - 8. * ( u(k,j,i+2) + u(k,j,i-1) )                 &
730                           +      ( u(k,j,i+3) + u(k,j,i-2) ) ) * adv_mom_5
731               diss_r(k) = - abs(u_comp(k) - gu) * (                          &
732                           10.  * ( u(k,j,i+1) - u(k,j,i) )                   &
733                           - 5. * ( u(k,j,i+2) - u(k,j,i-1) )                 &
734                           +      ( u(k,j,i+3) - u(k,j,i-2) ) ) *adv_mom_5
735            ENDDO
736           
737         ENDIF
738                 
739       ELSE
740!       
741!--       Compute the fifth order fluxes for the interior of PE domain.                 
742          DO  k = nzb_u_inner(j,i) + 1, nzt
743             u_comp(k) = u(k,j,i+1) + u(k,j,i)
744             flux_r(k) = ( u_comp(k) - gu ) * (                               &
745                         37.  * ( u(k,j,i+1) + u(k,j,i)   )                   &
746                         - 8. * ( u(k,j,i+2) + u(k,j,i-1) )                   &
747                         +      ( u(k,j,i+3) + u(k,j,i-2) ) ) * adv_mom_5
748             diss_r(k) = - abs(u_comp(k) - gu) * (                            &
749                         10.  * ( u(k,j,i+1) - u(k,j,i)   )                   &
750                         - 5. * ( u(k,j,i+2) - u(k,j,i-1) )                   &
751                         +      ( u(k,j,i+3) - u(k,j,i-2) ) ) * adv_mom_5
752
753             v_comp    = v(k,j+1,i) + v(k,j+1,i-1) - gv
754             flux_n(k) = v_comp * (                                           &
755                         37.  * ( u(k,j+1,i) + u(k,j,i)   )                   &
756                         - 8. * ( u(k,j+2,i) + u(k,j-1,i) )                   &
757                         +      ( u(k,j+3,i) + u(k,j-2,i) ) ) * adv_mom_5
758             diss_n(k) = - abs(v_comp) * (                                    &
759                         10.  * ( u(k,j+1,i) - u(k,j,i)   )                   &
760                         - 5. * ( u(k,j+2,i) - u(k,j-1,i) )                   &
761                         +      ( u(k,j+3,i) - u(k,j-2,i) ) ) * adv_mom_5
762          ENDDO
763         
764       ENDIF
765!       
766!--    Compute left- and southside fluxes for the respective boundary of PE     
767       IF ( j == nys .AND. ( .NOT. degraded_s ) )  THEN
768       
769          DO  k = nzb_u_inner(j,i) + 1, nzt
770             v_comp      = v(k,j,i) + v(k,j,i-1) - gv
771             flux_s_u(k) = v_comp * (                                         &
772                           37.  * ( u(k,j,i) + u(k,j-1,i)   )                 &
773                           - 8. * ( u(k,j+1,i) + u(k,j-2,i) )                 &
774                           +      ( u(k,j+2,i) + u(k,j-3,i) ) ) * adv_mom_5
775             diss_s_u(k) = - abs(v_comp) * (                                  &
776                           10.  * ( u(k,j,i) - u(k,j-1,i)    )                &
777                           - 5. * ( u(k,j+1,i) - u(k,j-2,i)  )                &
778                           +      (  u(k,j+2,i) - u(k,j-3,i) ) ) * adv_mom_5
779          ENDDO
780         
781       ENDIF
782       
783       IF ( i == nxlu .AND. ( .NOT. degraded_l ) )  THEN
784       
785          DO  k = nzb_u_inner(j,i)+1, nzt
786             u_comp_l      = u(k,j,i)+u(k,j,i-1)-gu
787             flux_l_u(k,j) = u_comp_l * (                                     &
788                             37.  * ( u(k,j,i) + u(k,j,i-1)   )               &
789                             - 8. * ( u(k,j,i+1) + u(k,j,i-2) )               &
790                             +      ( u(k,j,i+2) + u(k,j,i-3) ) ) * adv_mom_5
791             diss_l_u(k,j) = - abs(u_comp_l) * (                              &
792                             10.  * ( u(k,j,i) - u(k,j,i-1)   )               &
793                             - 5. * ( u(k,j,i+1) - u(k,j,i-2) )               &
794                             +      ( u(k,j,i+2) - u(k,j,i-3) ) ) * adv_mom_5
795          ENDDO
796         
797       ENDIF
798!       
799!--    Now compute the tendency terms for the horizontal parts.
800       DO  k = nzb_u_inner(j,i) + 1, nzt                   
801          tend(k,j,i) = tend(k,j,i) - (                                       &
802                            ( flux_r(k) + diss_r(k)                           &
803                          - ( flux_l_u(k,j) + diss_l_u(k,j) ) ) * ddx         &
804                          + ( flux_n(k) + diss_n(k)                           &
805                          - ( flux_s_u(k) + diss_s_u(k) )     ) * ddy  )
806
807           flux_l_u(k,j) = flux_r(k)
808           diss_l_u(k,j) = diss_r(k)
809           flux_s_u(k)   = flux_n(k)
810           diss_s_u(k)   = diss_n(k)
811!
812!--        Statistical Evaluation of u'u'. The factor has to be applied for
813!--        right evaluation when gallilei_trans = .T. .
814           sums_us2_ws_l(k,:) = sums_us2_ws_l(k,:)                            &
815             + ( flux_r(k) *                                                  &
816             ( u_comp(k) - 2. * hom(k,1,1,:) ) / ( u_comp(k) - gu + 1.0E-20 ) &
817             + diss_r(k)                                                      &
818             * abs(u_comp(k) - 2. * hom(k,1,1,:) )                            & 
819             / (abs(u_comp(k) - gu) + 1.0E-20)         )                      &
820             * weight_substep(intermediate_timestep_count) * rmask(j,i,:)
821       ENDDO   
822!
823!--    Vertical advection, degradation of order near surface and top.
824!--    The fluxes flux_d and diss_d at the surface are 0. Due to reasons of
825!--    statistical evaluation the top flux at the surface should be 0
826       flux_t(nzb_u_inner(j,i)) = 0. !statistical reasons
827       diss_t(nzb_u_inner(j,i)) = 0.
828!
829!--    2nd order scheme
830       k         = nzb_u_inner(j,i) + 1
831       flux_d    = flux_t(k-1)
832       diss_d    = diss_t(k-1)
833       w_comp    = w(k,j,i) + w(k,j,i-1)
834       flux_t(k) = w_comp * ( u(k+1,j,i) + u(k,j,i) ) *0.25
835       diss_t(k) = diss_2nd( u(k+2,j,i), u(k+1,j,i), u(k,j,i), 0., 0.,        &
836                             w_comp, 0.25, ddzw(k) )
837             
838       tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)                    &
839                                 - ( flux_d + diss_d )     ) * ddzw(k)
840!
841!--    WS3 as an intermediate step
842       k         = nzb_u_inner(j,i) + 2
843       flux_d    = flux_t(k-1)
844       diss_d    = diss_t(k-1)
845       w_comp    = w(k,j,i) + w(k,j,i-1)
846       flux_t(k) = w_comp*(                                                   &
847                   7. * ( u(k+1,j,i) + u(k,j,i)   )                           &
848                   -    ( u(k+2,j,i) + u(k-1,j,i) ) ) * adv_mom_3
849       diss_t(k) = -abs(w_comp)*(                                             &
850                   3. * ( u(k+1,j,i) - u(k,j,i)   )                           & 
851                   -    ( u(k+2,j,i) - u(k-1,j,i) ) ) * adv_mom_3
852
853       tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)                    &
854                                 - ( flux_d + diss_d )     ) * ddzw(k)
855
856       DO  k = nzb_u_inner(j,i) + 3, nzt - 2
857          flux_d    = flux_t(k-1)
858          diss_d    = diss_t(k-1)
859          w_comp    = w(k,j,i) + w(k,j,i-1)
860          flux_t(k) = w_comp*(                                                &
861                      37.*   ( u(k+1,j,i) + u(k,j,i)   )                      &
862                      - 8. * ( u(k+2,j,i) + u(k-1,j,i) )                      &
863                      +      ( u(k+3,j,i) + u(k-2,j,i) ) ) * adv_mom_5
864          diss_t(k) = - abs(w_comp) * (                                       &
865                      10. *  ( u(k+1,j,i) - u(k,j,i)   )                      & 
866                      - 5. * ( u(k+2,j,i) - u(k-1,j,i) )                      &
867                      +      ( u(k+3,j,i) - u(k-2,j,i) ) ) * adv_mom_5
868
869          tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)                 &
870                                 - ( flux_d + diss_d )     ) * ddzw(k)
871       ENDDO
872!
873!--    WS3 as an intermediate step
874       k         = nzt - 1
875       flux_d    = flux_t(k-1)
876       diss_d    = diss_t(k-1)
877       w_comp    = w(k,j,i) + w(k,j,i-1)
878       flux_t(k) = w_comp * (                                                 &
879                   7. * ( u(k+1,j,i) + u(k,j,i)   )                           &
880                   -    ( u(k+2,j,i) + u(k-1,j,i) ) ) * adv_mom_3
881       diss_t(k) = - abs(w_comp) * (                                          &
882                   3. * ( u(k+1,j,i) - u(k,j,i)   )                           & 
883                   -    ( u(k+2,j,i) - u(k-1,j,i) ) ) * adv_mom_3
884                   
885       tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)                    &
886                                 - ( flux_d + diss_d )     ) * ddzw(k)
887       
888!
889!--    2nd order scheme
890       k         = nzt
891       flux_d    = flux_t(k-1)
892       diss_d    = diss_t(k-1)
893       w_comp    = w(k,j,i)+w(k,j,i-1)
894       flux_t(k) = w_comp * ( u(k+1,j,i) + u(k,j,i) ) * 0.25
895       diss_t(k) = diss_2nd( u(k+1,j,i), u(k+1,j,i), u(k,j,i), u(k-1,j,i),    &
896                             u(k-2,j,i), w_comp, 0.25, ddzw(k) )
897
898       tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)                    &
899                                 - ( flux_d + diss_d )     ) * ddzw(k)
900
901!
902!--    sum up the vertical momentum fluxes
903       DO  k = nzb_u_inner(j,i), nzt
904          sums_wsus_ws_l(k,:) = sums_wsus_ws_l(k,:)                           &
905              + ( flux_t(k) + diss_t(k) )                                     &
906              * weight_substep(intermediate_timestep_count) * rmask(j,i,:)
907       ENDDO
908
909
910    END SUBROUTINE advec_u_ws_ij
911
912
913
914
915!------------------------------------------------------------------------------!
916! Advection of v-component - Call for grid point i,j
917!------------------------------------------------------------------------------!
918   SUBROUTINE advec_v_ws_ij( i, j )
919
920       USE arrays_3d
921       USE constants
922       USE control_parameters
923       USE grid_variables
924       USE indices
925       USE statistics
926
927       IMPLICIT NONE
928
929       INTEGER  ::  i, j, k
930       LOGICAL  ::  degraded_l = .FALSE., degraded_s = .FALSE.
931       REAL     ::  gu, gv, flux_d, diss_d, u_comp, v_comp_l, w_comp
932       REAL, DIMENSION(nzb:nzt+1)  ::  flux_t, diss_t, flux_n,                &
933                                       diss_n, flux_r, diss_r, v_comp
934     
935       gu = 2.0 * u_gtrans
936       gv = 2.0 * v_gtrans
937       
938       IF ( boundary_flags(j,i) /= 0 )  THEN
939!       
940!--       Degrade the order for Dirichlet bc. at the outflow boundary
941          SELECT CASE ( boundary_flags(j,i) )
942         
943             CASE ( 1 )
944                DO  k = nzb_v_inner(j,i) + 1, nzt
945                   u_comp    = u(k,j-1,i+1) + u(k,j,i+1) - gu
946                   flux_r(k) = u_comp * (                                     &
947                               7. * (v(k,j,i+1) + v(k,j,i)    )               &
948                               -    ( v(k,j,i+2) + v(k,j,i-1) ) ) * adv_mom_3
949                   diss_r(k) = - abs(u_comp) * (                              &
950                               3. * ( v(k,j,i+1) - v(k,j,i)   )               & 
951                               -    ( v(k,j,i+2) - v(k,j,i-1) ) ) * adv_mom_3
952                ENDDO
953               
954             CASE ( 2 )
955                DO  k = nzb_v_inner(j,i) + 1, nzt
956                   u_comp    = u(k,j-1,i+1) + u(k,j,i+1) - gu
957                   flux_r(k) = u_comp * ( v(k,j,i+1) + v(k,j,i) ) * 0.25 
958                   diss_r(k) = diss_2nd( v(k,j,i+1), v(k,j,i+1), v(k,j,i),    &
959                                         v(k,j,i-1), v(k,j,i-2), u_comp,      &
960                                         0.25, ddx )
961                ENDDO
962               
963             CASE ( 3 )
964                DO  k = nzb_v_inner(j,i) + 1, nzt
965                   v_comp(k) = v(k,j+1,i) + v(k,j,i)
966                   flux_n(k) = ( v_comp(k)- gv ) * (                          &
967                               7. * ( v(k,j+1,i) + v(k,j,i)   )               &
968                               -    ( v(k,j+2,i) + v(k,j-1,i) ) ) * adv_mom_3
969                   diss_n(k) = - abs(v_comp(k) - gv) * (                      &
970                               3. * ( v(k,j+1,i) - v(k,j,i)   )               & 
971                               -    ( v(k,j+2,i) - v(k,j-1,i) ) ) * adv_mom_3
972                ENDDO
973               
974             CASE ( 4 )
975                DO  k = nzb_v_inner(j,i) + 1, nzt
976                   v_comp(k) = v(k,j+1,i) + v(k,j,i)
977                   flux_n(k) = ( v_comp(k) - gv ) *                           &
978                               ( v(k,j+1,i) + v(k,j,i) ) * 0.25 
979                   diss_n(k) = diss_2nd( v(k,j+1,i), v(k,j+1,i), v(k,j,i),    &
980                                         v(k,j-1,i), v(k,j-2,i), v_comp(k),   &
981                                         0.25, ddy )
982                ENDDO
983               
984             CASE ( 5 )
985                DO  k = nzb_w_inner(j,i) + 1, nzt
986                   u_comp    = u(k,j-1,i) + u(k,j,i) - gu
987                   flux_r(k) = u_comp * (                                     &
988                               7. * ( v(k,j,i+1) + v(k,j,i)   )               &
989                               -    ( v(k,j,i+2) + v(k,j,i-1) ) ) * adv_mom_3
990                   diss_r(k) = - abs(u_comp) * (                              &
991                               3. * ( w(k,j,i+1) - w(k,j,i)   )               & 
992                               -    ( v(k,j,i+2) - v(k,j,i-1) ) ) * adv_mom_3
993                ENDDO
994               
995             CASE ( 6 )
996                DO  k = nzb_v_inner(j,i) + 1, nzt
997                   u_comp        = u(k,j-1,i) + u(k,j,i) - gu
998                   flux_l_v(k,j) = u_comp * ( v(k,j,i) + v(k,j,i-1) ) * 0.25
999                   diss_l_v(k,j) = diss_2nd( v(k,j,i+2), v(k,j,i+1), v(k,j,i),&
1000                                             v(k,j,i-1), v(k,j,i-1), u_comp,  &
1001                                             0.25, ddx )
1002                                 
1003                   u_comp = u(k,j-1,i+1) + u(k,j,i+1) - gu
1004                   flux_r(k) = u_comp * (                                     &
1005                               7. * ( v(k,j,i+1) + v(k,j,i)   )               &
1006                               -    ( v(k,j,i+2) + v(k,j,i-1) ) ) * adv_mom_3
1007                   diss_r(k) = - abs(u_comp) * (                              &
1008                               3. * ( v(k,j,i+1) - v(k,j,i)   )               & 
1009                               -    ( v(k,j,i+2) - v(k,j,i-1) ) ) * adv_mom_3
1010                ENDDO
1011                degraded_l = .TRUE.
1012               
1013             CASE ( 7 )
1014                DO  k = nzb_v_inner(j,i) + 1, nzt
1015                   v_comp(k)   = v(k,j,i) + v(k,j-1,i) - gv
1016                   flux_s_v(k) = v_comp(k) * ( v(k,j,i) + v(k,j-1,i) ) * 0.25
1017                   diss_s_v(k) = diss_2nd( v(k,j+2,i), v(k,j+1,i), v(k,j,i),  &
1018                                           v(k,j-1,i), v(k,j-1,i), v_comp(k), &
1019                                           0.25, ddy )
1020                                 
1021                   v_comp(k) = v(k,j+1,i) + v(k,j,i)
1022                   flux_n(k) = ( v_comp(k) - gv ) * (                         &
1023                               7. * ( v(k,j+1,i) + v(k,j,i)   )               &
1024                               -    ( v(k,j+2,i) + v(k,j-1,i) ) ) * adv_mom_3
1025                   diss_n(k) = - abs(v_comp(k) - gv) * (                      &
1026                               3. * ( v(k,j+1,i) - v(k,j,i)   )               & 
1027                               -    ( v(k,j+2,i) - v(k,j-1,i) ) ) * adv_mom_3
1028                ENDDO
1029                degraded_s = .TRUE.
1030               
1031              CASE DEFAULT
1032             
1033          END SELECT
1034!         
1035!--       Compute the crosswise 5th order fluxes at the outflow
1036          IF ( boundary_flags(j,i) == 1 .OR. boundary_flags(j,i) == 2         &
1037          .OR. boundary_flags(j,i) == 5 .OR. boundary_flags(j,i) == 6 )  THEN
1038         
1039             DO  k = nzb_v_inner(j,i) + 1, nzt
1040                v_comp(k) = v(k,j+1,i) + v(k,j,i)
1041                flux_n(k) = ( v_comp(k) - gv ) * (                            &
1042                            37.  * ( v(k,j+1,i) + v(k,j,i)   )                &
1043                            - 8. * ( v(k,j+2,i) + v(k,j-1,i) )                &
1044                            +      ( v(k,j+3,i) + v(k,j-2,i) ) ) * adv_mom_5
1045                diss_n(k) = - abs(v_comp(k) - gv ) * (                        &
1046                             10.  * ( v(k,j+1,i) - v(k,j,i)   )               &
1047                             - 5. * ( v(k,j+2,i) - v(k,j-1,i) )               &
1048                             +      ( v(k,j+3,i) - v(k,j-2,i) ) ) *adv_mom_5
1049              ENDDO
1050             
1051          ELSE
1052         
1053             DO  k = nzb_v_inner(j,i) + 1, nzt
1054                u_comp    = u(k,j-1,i+1) + u(k,j,i+1) - gu
1055                flux_r(k) = u_comp * (                                        &
1056                            37.  * ( v(k,j,i+1) + v(k,j,i)   )                &
1057                            - 8. * ( v(k,j,i+2) + v(k,j,i-1) )                &
1058                            +      ( v(k,j,i+3) + v(k,j,i-2) ) ) * adv_mom_5
1059                diss_r(k) = - abs(u_comp) * (                                 &
1060                            10. * ( v(k,j,i+1) - v(k,j,i)   )                 &
1061                            -5. * ( v(k,j,i+2) - v(k,j,i-1) )                 &
1062                            +     ( v(k,j,i+3) - v(k,j,i-2) ) ) * adv_mom_5
1063             ENDDO
1064             
1065          ENDIF
1066         
1067       ELSE
1068!       
1069!--       Compute the fifth order fluxes for the interior of PE domain.                 
1070          DO  k = nzb_v_inner(j,i) + 1, nzt
1071             u_comp    = u(k,j-1,i+1) + u(k,j,i+1) - gu
1072             flux_r(k) = u_comp * (                                           &
1073                         37.  * ( v(k,j,i+1) + v(k,j,i)   )                   &
1074                         - 8. * ( v(k,j,i+2) + v(k,j,i-1) )                   &
1075                         +      ( v(k,j,i+3) + v(k,j,i-2) ) ) * adv_mom_5
1076             diss_r(k) = - abs(u_comp) * (                                    &
1077                         10.  * ( v(k,j,i+1) - v(k,j,i) )                     &
1078                         - 5. * ( v(k,j,i+2) - v(k,j,i-1) )                   &
1079                         +      ( v(k,j,i+3) - v(k,j,i-2) ) ) * adv_mom_5
1080
1081             v_comp(k) = v(k,j+1,i) + v(k,j,i)
1082             flux_n(k) = ( v_comp(k) - gv ) * (                               &
1083                         37.  * ( v(k,j+1,i) + v(k,j,i)   )                   &
1084                         - 8. * ( v(k,j+2,i) + v(k,j-1,i) )                   &
1085                         +      ( v(k,j+3,i) + v(k,j-2,i) ) ) * adv_mom_5
1086             diss_n(k) = - abs(v_comp(k) - gv) * (                            &
1087                         10.  * ( v(k,j+1,i) - v(k,j,i)   )                   &
1088                         - 5. * ( v(k,j+2,i) - v(k,j-1,i) )                   &
1089                         +      ( v(k,j+3,i) - v(k,j-2,i) ) ) * adv_mom_5
1090          ENDDO
1091         
1092       ENDIF
1093!       
1094!--    Compute left- and southside fluxes for the respective boundary       
1095       IF ( i == nxl .AND. ( .NOT. degraded_l ) )  THEN
1096          DO  k = nzb_v_inner(j,i) + 1, nzt
1097             u_comp        = u(k,j-1,i) + u(k,j,i) - gu
1098             flux_l_v(k,j) = u_comp * (                                       &
1099                             37.  * ( v(k,j,i) + v(k,j,i-1)   )               &
1100                             - 8. * ( v(k,j,i+1) + v(k,j,i-2) )               &
1101                             +      ( v(k,j,i+2) + v(k,j,i-3) ) ) * adv_mom_5
1102             diss_l_v(k,j) = - abs(u_comp) * (                                &
1103                             10.  * ( v(k,j,i) - v(k,j,i-1)   )               &
1104                             - 5. * ( v(k,j,i+1) - v(k,j,i-2) )               &
1105                             +      ( v(k,j,i+2) - v(k,j,i-3) ) ) * adv_mom_5
1106          ENDDO
1107         
1108       ENDIF
1109       
1110       IF ( j == nysv .AND. ( .NOT. degraded_s ) )  THEN
1111       
1112          DO  k = nzb_v_inner(j,i) + 1, nzt
1113             v_comp_l    = v(k,j,i) + v(k,j-1,i) - gv
1114             flux_s_v(k) = v_comp_l * (                                       &
1115                           37.  * ( v(k,j,i) + v(k,j-1,i)   )                 &
1116                           - 8. * ( v(k,j+1,i) + v(k,j-2,i) )                 &
1117                           +      ( v(k,j+2,i) + v(k,j-3,i) ) ) * adv_mom_5
1118             diss_s_v(k) = - abs(v_comp_l) * (                                &
1119                           10.  * ( v(k,j,i) - v(k,j-1,i)   )                 &
1120                           - 5. * ( v(k,j+1,i) - v(k,j-2,i) )                 &
1121                           +      ( v(k,j+2,i) - v(k,j-3,i) ) ) * adv_mom_5
1122          ENDDO
1123         
1124       ENDIF
1125!       
1126!--    Now compute the tendency terms for the horizontal parts.         
1127       DO  k = nzb_v_inner(j,i) + 1, nzt                 
1128          tend(k,j,i) = tend(k,j,i) - (                                       &
1129                         ( flux_r(k) + diss_r(k)                              &
1130                       - ( flux_l_v(k,j) + diss_l_v(k,j) ) ) * ddx            &
1131                       + ( flux_n(k) + diss_n(k)                              &
1132                       - ( flux_s_v(k) + diss_s_v(k)     ) ) * ddy     )
1133       
1134           flux_l_v(k,j) = flux_r(k)
1135           diss_l_v(k,j) = diss_r(k)
1136           flux_s_v(k)   = flux_n(k)
1137           diss_s_v(k)   = diss_n(k)
1138
1139!
1140!--        Statistical Evaluation of v'v'. The factor has to be applied for
1141!--        right evaluation when gallilei_trans = .T. .
1142
1143           sums_vs2_ws_l(k,:) = sums_vs2_ws_l(k,:)                            &
1144             + ( flux_n(k)                                                    &
1145             *  ( v_comp(k) - 2. * hom(k,1,2,:) )                             &
1146             / ( v_comp(k) - gv + 1.0E-20 )                                   &
1147             + diss_n(k)                                                      &
1148             * abs( v_comp(k) - 2. * hom(k,1,2,:) )                           &
1149             / ( abs( v_comp(k) - gv ) +1.0E-20 ) )                           &
1150             * weight_substep(intermediate_timestep_count) * rmask(j,i,:)
1151
1152       ENDDO
1153!
1154!--    Vertical advection, degradation of order near surface and top.
1155!--    The fluxes flux_d and diss_d at the surface are 0. Due to reasons of
1156!--    statistical evaluation the top flux at the surface should be 0
1157       flux_t(nzb_v_inner(j,i)) = 0.  !statistical reasons
1158       diss_t(nzb_v_inner(j,i)) = 0.
1159!
1160!--    2nd order scheme       
1161       k         = nzb_v_inner(j,i) + 1
1162       flux_d    = flux_t(k-1)
1163       diss_d    = diss_t(k-1) 
1164       w_comp    = w(k,j-1,i) + w(k,j,i)
1165       flux_t(k) = w_comp * ( v(k+1,j,i) + v(k,j,i) ) * 0.25
1166       diss_t(k) = diss_2nd( v(k+2,j,i), v(k+1,j,i), v(k,j,i), 0., 0., w_comp,&
1167                             0.25, ddzw(k) ) 
1168
1169       tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)                    &
1170                                 - ( flux_d + diss_d )     ) * ddzw(k)
1171       
1172!
1173!--    WS3 as an intermediate step
1174       k         = nzb_v_inner(j,i) + 2
1175       flux_d    = flux_t(k-1)
1176       diss_d    = diss_t(k-1)
1177       w_comp    = w(k,j-1,i) + w(k,j,i)
1178       flux_t(k) = w_comp * (                                                 &
1179                   7. * ( v(k+1,j,i) + v(k,j,i)   )                           &
1180                   -    ( v(k+2,j,i) + v(k-1,j,i) ) ) * adv_mom_3
1181       diss_t(k) = - abs(w_comp) * (                                          &
1182                   3. * ( v(k+1,j,i) - v(k,j,i)   )                           & 
1183                   -    ( v(k+2,j,i) - v(k-1,j,i) ) ) * adv_mom_3
1184
1185       tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)                    &
1186                                 - ( flux_d + diss_d )     ) * ddzw(k)
1187!--    WS5
1188       DO  k = nzb_v_inner(j,i) + 3, nzt - 2
1189          flux_d    = flux_t(k-1)
1190          diss_d    = diss_t(k-1)
1191          w_comp    = w(k,j-1,i) + w(k,j,i)
1192          flux_t(k) = w_comp * (                                              &
1193                      37.  * ( v(k+1,j,i) + v(k,j,i)   )                      &           
1194                      - 8. * ( v(k+2,j,i) + v(k-1,j,i) )                      &
1195                      +      ( v(k+3,j,i) + v(k-2,j,i) ) ) * adv_mom_5
1196          diss_t(k) = - abs(w_comp) * (                                       &
1197                      10.  * ( v(k+1,j,i) - v(k,j,i)   )                      & 
1198                      - 5. * ( v(k+2,j,i) - v(k-1,j,i) )                      &
1199                      +      ( v(k+3,j,i) - v(k-2,j,i) ) ) * adv_mom_5
1200
1201          tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)                 &
1202                                 - ( flux_d + diss_d )     ) * ddzw(k)
1203       ENDDO
1204!
1205!--    WS3 as an intermediate step
1206       k         = nzt - 1
1207       flux_d    = flux_t(k-1)
1208       diss_d    = diss_t(k-1)
1209       w_comp    = w(k,j-1,i) + w(k,j,i)
1210       flux_t(k) = w_comp * (                                                 &
1211                   7. * ( v(k+1,j,i) + v(k,j,i)   )                           &   
1212                   -    ( v(k+2,j,i) + v(k-1,j,i) ) ) * adv_mom_3
1213       diss_t(k) = - abs(w_comp) * (                                          &
1214                   3. * ( v(k+1,j,i) - v(k,j,i) )                             & 
1215                   -    ( v(k+2,j,i) - v(k-1,j,i) ) ) * adv_mom_3
1216       tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)                    &
1217                                 - ( flux_d + diss_d )     ) * ddzw(k)
1218!
1219!--    2nd order scheme
1220       k         = nzt
1221       flux_d    = flux_t(k-1)
1222       diss_d    = diss_t(k-1)
1223       w_comp    = w(k,j-1,i)+w(k,j,i)
1224       flux_t(k) = w_comp * ( v(k+1,j,i) + v(k,j,i) ) * 0.25
1225       diss_t(k) = diss_2nd( v(k+1,j,i), v(k+1,j,i), v(k,j,i), v(k-1,j,i),    &
1226                             v(k-2,j,i), w_comp, 0.25, ddzw(k) )
1227 
1228       tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)                    &
1229                                 - ( flux_d + diss_d )     ) * ddzw(k)
1230             
1231       DO  k = nzb_v_inner(j,i), nzt
1232          sums_wsvs_ws_l(k,:) = sums_wsvs_ws_l(k,:)                           &
1233            + ( flux_t(k) + diss_t(k) )                                       &
1234            * weight_substep(intermediate_timestep_count) * rmask(j,i,:)
1235       ENDDO
1236
1237    END SUBROUTINE advec_v_ws_ij
1238
1239
1240
1241!------------------------------------------------------------------------------!
1242! Advection of w-component - Call for grid point i,j
1243!------------------------------------------------------------------------------!
1244    SUBROUTINE advec_w_ws_ij( i, j )
1245
1246       USE arrays_3d
1247       USE constants
1248       USE control_parameters
1249       USE grid_variables
1250       USE indices
1251       USE statistics
1252
1253       IMPLICIT NONE
1254
1255       INTEGER ::  i, j, k
1256       LOGICAL ::  degraded_l = .FALSE., degraded_s = .FALSE.
1257       REAL    ::  gu, gv, flux_d, diss_d, u_comp, v_comp, w_comp
1258       REAL, DIMENSION(nzb:nzt+1)  ::  flux_t, diss_t, flux_r, diss_r, flux_n, &
1259                                       diss_n
1260
1261       gu = 2.0 * u_gtrans
1262       gv = 2.0 * v_gtrans
1263       
1264       IF ( boundary_flags(j,i) /= 0 )  THEN
1265!       
1266!--      Degrade the order for Dirichlet bc. at the outflow boundary
1267         SELECT CASE ( boundary_flags(j,i) )
1268         
1269            CASE ( 1 )
1270               DO  k = nzb_w_inner(j,i) + 1, nzt
1271                  u_comp    = u(k+1,j,i+1) + u(k,j,i+1) - gu
1272                  flux_r(k) = u_comp * (                                      &
1273                              7. * ( w(k,j,i+1) + w(k,j,i)   )                &
1274                              -    ( w(k,j,i+2) + w(k,j,i-1) ) ) * adv_mom_3
1275                  diss_r(k) = -abs(u_comp) * (                                &
1276                              3. * ( w(k,j,i+1) - w(k,j,i)   )                & 
1277                              -    ( w(k,j,i+2) - w(k,j,i-1) ) ) * adv_mom_3 
1278               ENDDO
1279               
1280            CASE ( 2 )
1281               DO  k = nzb_w_inner(j,i) + 1, nzt
1282                  u_comp    = u(k+1,j,i+1) + u(k,j,i+1) - gu
1283                  flux_r(k) = u_comp * ( w(k,j,i+1) + w(k,j,i) ) * 0.25 
1284                  diss_r(k) = diss_2nd( w(k,j,i+1), w(k,j,i+1), w(k,j,i),     &
1285                                        w(k,j,i-1), w(k,j,i-2), u_comp,       &
1286                                        0.25, ddx )
1287               ENDDO
1288               
1289            CASE ( 3 )
1290               DO  k = nzb_w_inner(j,i) + 1, nzt
1291                  v_comp    = v(k+1,j+1,i) + v(k,j+1,i) - gv
1292                  flux_n(k) = v_comp * (                                      &
1293                              7. * ( w(k,j+1,i) + w(k,j,i)   )                &
1294                              -    ( w(k,j+2,i) + w(k,j-1,i) ) ) * adv_mom_3
1295                  diss_n(k) = -abs(v_comp) * (                                &
1296                              3. * ( w(k,j+1,i) - w(k,j,i)   )                & 
1297                              -    ( w(k,j+2,i) - w(k,j-1,i) ) ) * adv_mom_3
1298               ENDDO
1299               
1300            CASE ( 4 )
1301               DO  k = nzb_w_inner(j,i) + 1, nzt
1302                  v_comp    = v(k+1,j+1,i) + v(k,j+1,i) - gv
1303                  flux_n(k) = v_comp * ( w(k,j+1,i) + w(k,j,i) ) * 0.25 
1304                  diss_n(k) = diss_2nd( w(k,j+1,i), w(k,j+1,i), w(k,j,i),     &
1305                                        w(k,j-1,i), w(k,j-2,i), v_comp,       &
1306                                        0.25, ddy )
1307               ENDDO
1308               
1309            CASE ( 5 )
1310               DO  k = nzb_w_inner(j,i) + 1, nzt
1311                  u_comp    = u(k+1,j,i+1) + u(k,j,i+1) - gu
1312                  flux_r(k) = u_comp * (                                      &
1313                              7. * ( w(k,j,i+1) + w(k,j,i)   )                &
1314                              -    ( w(k,j,i+2) + w(k,j,i-1) ) ) * adv_mom_3
1315                  diss_r(k) = - abs(u_comp) * (                               &
1316                              3. * ( w(k,j,i+1) - w(k,j,i) )                  & 
1317                              -    ( w(k,j,i+2) - w(k,j,i-1) ) ) * adv_mom_3
1318               ENDDO
1319               
1320            CASE ( 6 )
1321               DO  k = nzb_w_inner(j,i) + 1, nzt
1322!               
1323!--               Compute leftside fluxes for the left boundary of PE domain
1324                  u_comp    = u(k+1,j,i+1) + u(k,j,i+1) - gu
1325                  flux_r(k) = u_comp *(                                       &
1326                              7. * ( w(k,j,i+1) + w(k,j,i)   )                &
1327                              -    ( w(k,j,i+2) + w(k,j,i-1) ) )  * adv_mom_3
1328                  diss_r(k) = - abs(u_comp) * (                               &
1329                              3. * ( w(k,j,i+1) - w(k,j,i) )                  & 
1330                              -    ( w(k,j,i+2) - w(k,j,i-1) ) ) * adv_mom_3
1331                 
1332                  u_comp        = u(k+1,j,i) + u(k,j,i) - gu
1333                  flux_l_w(k,j) = u_comp * ( w(k,j,i) + w(k,j,i-1) ) * 0.25
1334                  diss_l_w(k,j) = diss_2nd( w(k,j,i+2), w(k,j,i+1), w(k,j,i), &
1335                                            w(k,j,i-1), w(k,j,i-1), u_comp,   &
1336                                            0.25,ddx)
1337               ENDDO
1338               degraded_l = .TRUE.
1339               
1340            CASE ( 7 )
1341               DO  k = nzb_w_inner(j,i) + 1, nzt
1342                  v_comp    = v(k+1,j+1,i) + v(k,j+1,i) - gv
1343                  flux_n(k) = v_comp *(                                       &
1344                              7. * ( w(k,j+1,i) + w(k,j,i)   )                &
1345                              -    ( w(k,j+2,i) + w(k,j-1,i) ) ) * adv_mom_3
1346                  diss_n(k) = - abs(v_comp) * (                               &
1347                              3. * ( w(k,j+1,i) - w(k,j,i)   )                & 
1348                              -    ( w(k,j+2,i) - w(k,j-1,i) ) ) * adv_mom_3
1349                ENDDO
1350               
1351            CASE ( 8 )
1352               DO  k = nzb_w_inner(j,i) + 1, nzt
1353                  v_comp    = v(k+1,j+1,i) + v(k,j+1,i) - gv
1354                  flux_n(k) = v_comp * (                                      &
1355                              7. * ( w(k,j+1,i) + w(k,j,i)   )                &
1356                             -     ( w(k,j+2,i) + w(k,j-1,i) ) ) * adv_mom_3
1357                  diss_n(k) = - abs(v_comp) * (                               &
1358                              3. * ( w(k,j+1,i) - w(k,j,i) )                  & 
1359                              -    ( w(k,j+2,i) - w(k,j-1,i) ) ) * adv_mom_3
1360!                   
1361!--              Compute southside fluxes for the south boundary of PE domain           
1362                  v_comp      = v(k+1,j,i) + v(k,j,i) - gv
1363                  flux_s_w(k) = v_comp * ( w(k,j,i) + w(k,j-1,i) ) * 0.25
1364                  diss_s_w(k) = diss_2nd( w(k,j+2,i), w(k,j+1,i), w(k,j,i),   &
1365                                          w(k,j-1,i), w(k,j-1,i), v_comp,     &
1366                                          0.25, ddy )                 
1367               ENDDO 
1368               degraded_s = .TRUE.
1369               
1370            CASE DEFAULT
1371           
1372         END SELECT
1373!         
1374!--      Compute the crosswise 5th order fluxes at the outflow
1375         IF ( boundary_flags(j,i) == 1 .OR. boundary_flags(j,i) == 2          &
1376         .OR. boundary_flags(j,i) == 5 .OR. boundary_flags(j,i) == 6 )  THEN
1377         
1378            DO  k = nzb_w_inner(j,i) + 1, nzt
1379               v_comp    = v(k+1,j+1,i) + v(k,j+1,i) - gv
1380               flux_n(k) = v_comp * (                                         &
1381                           37.  * ( w(k,j+1,i) + w(k,j,i)   )                 &
1382                           - 8. * ( w(k,j+2,i) + w(k,j-1,i) )                 &
1383                           +      ( w(k,j+3,i) + w(k,j-2,i) ) ) * adv_mom_5
1384               diss_n(k) = - abs(v_comp) * (                                  &
1385                           10.  * ( w(k,j+1,i) - w(k,j,i)   )                 &
1386                           - 5. * ( w(k,j+2,i) - w(k,j-1,i) )                 &
1387                           +      ( w(k,j+3,i) - w(k,j-2,i) ) ) * adv_mom_5
1388            ENDDO
1389           
1390         ELSE
1391         
1392            DO  k = nzb_w_inner(j,i) + 1, nzt         
1393               u_comp    = u(k+1,j,i+1) + u(k,j,i+1) - gu
1394               flux_r(k) = u_comp * (                                         &
1395                           37.  * ( w(k,j,i+1) + w(k,j,i) )                   &
1396                           - 8. * ( w(k,j,i+2) + w(k,j,i-1) )                 &
1397                           +      ( w(k,j,i+3) + w(k,j,i-2) ) ) * adv_mom_5
1398               diss_r(k) = - abs(u_comp) * (                                  &
1399                           10.  * ( w(k,j,i+1) - w(k,j,i)   )                 &
1400                           - 5. * ( w(k,j,i+2) - w(k,j,i-1) )                 &
1401                           +      ( w(k,j,i+3) - w(k,j,i-2) ) ) * adv_mom_5
1402            ENDDO
1403           
1404         ENDIF
1405         
1406       ELSE
1407!       
1408!--      Compute the fifth order fluxes for the interior of PE domain.               
1409         DO  k = nzb_w_inner(j,i) + 1, nzt
1410            u_comp    = u(k+1,j,i+1) + u(k,j,i+1) - gu
1411            flux_r(k) = u_comp * (                                            &
1412                        37.  * ( w(k,j,i+1) + w(k,j,i)   )                    &
1413                        - 8. * ( w(k,j,i+2) + w(k,j,i-1) )                    &
1414                        +      ( w(k,j,i+3) + w(k,j,i-2) ) ) * adv_mom_5
1415            diss_r(k) = - abs(u_comp) * (                                     &
1416                        10.  * ( w(k,j,i+1) - w(k,j,i)   )                    &
1417                        - 5. * ( w(k,j,i+2) - w(k,j,i-1) )                    &
1418                        +      ( w(k,j,i+3) - w(k,j,i-2) ) ) * adv_mom_5
1419                 
1420            v_comp    = v(k+1,j+1,i) + v(k,j+1,i) - gv
1421            flux_n(k) = v_comp * (                                            &
1422                        37.  * ( w(k,j+1,i) + w(k,j,i)   )                    &
1423                        - 8. * ( w(k,j+2,i) + w(k,j-1,i) )                    &
1424                        +      ( w(k,j+3,i) + w(k,j-2,i) ) ) * adv_mom_5
1425            diss_n(k) = - abs(v_comp) * (                                     &
1426                        10.  * ( w(k,j+1,i) - w(k,j,i)   )                    &
1427                        - 5. * ( w(k,j+2,i) - w(k,j-1,i) )                    &
1428                        +      ( w(k,j+3,i) - w(k,j-2,i) ) ) * adv_mom_5
1429         ENDDO
1430         
1431       ENDIF
1432!       
1433!--    Compute left- and southside fluxes for the respective boundary     
1434       IF ( j == nys .AND. ( .NOT. degraded_s ) )  THEN
1435       
1436          DO  k = nzb_w_inner(j,i) + 1, nzt
1437             v_comp      = v(k+1,j,i) + v(k,j,i) - gv
1438             flux_s_w(k) = v_comp * (                                         &
1439                           37.  * ( w(k,j,i) + w(k,j-1,i)   )                 &
1440                           - 8. * ( w(k,j+1,i) +w(k,j-2,i)  )                 &
1441                           +      ( w(k,j+2,i) + w(k,j-3,i) ) ) * adv_mom_5
1442             diss_s_w(k) = - abs(v_comp) * (                                  &
1443                           10.  * ( w(k,j,i) - w(k,j-1,i)   )                 &
1444                           - 5. * ( w(k,j+1,i) - w(k,j-2,i) )                 &
1445                           +      ( w(k,j+2,i) - w(k,j-3,i) ) ) * adv_mom_5
1446          ENDDO
1447         
1448       ENDIF
1449       
1450       IF ( i == nxl .AND. ( .NOT. degraded_l ) ) THEN
1451       
1452          DO  k = nzb_w_inner(j,i) + 1, nzt
1453            u_comp        = u(k+1,j,i) + u(k,j,i) - gu
1454            flux_l_w(k,j) = u_comp * (                                        &
1455                            37.  * ( w(k,j,i) + w(k,j,i-1)   )                &
1456                            - 8. * ( w(k,j,i+1) + w(k,j,i-2) )                &
1457                            +      ( w(k,j,i+2) + w(k,j,i-3) ) ) * adv_mom_5
1458            diss_l_w(k,j) = - abs(u_comp) * (                                 &
1459                            10.  * ( w(k,j,i) - w(k,j,i-1)   )                &
1460                            - 5. * ( w(k,j,i+1) - w(k,j,i-2) )                &
1461                            +      ( w(k,j,i+2) - w(k,j,i-3) ) ) * adv_mom_5
1462          ENDDO
1463         
1464       ENDIF
1465!       
1466!--    Now compute the tendency terms for the horizontal parts.         
1467       DO  k = nzb_w_inner(j,i) + 1, nzt
1468          tend(k,j,i) = tend(k,j,i) - (                                       &
1469                      ( flux_r(k) + diss_r(k)                                 &
1470                    - ( flux_l_w(k,j) + diss_l_w(k,j) ) ) * ddx               &
1471                    + ( flux_n(k) + diss_n(k)                                 &
1472                    - ( flux_s_w(k) + diss_s_w(k)     ) ) * ddy  )
1473
1474          flux_l_w(k,j) = flux_r(k)
1475          diss_l_w(k,j) = diss_r(k)
1476          flux_s_w(k) = flux_n(k)
1477          diss_s_w(k) = diss_n(k)
1478       ENDDO
1479
1480!
1481!--    Vertical advection, degradation of order near surface and top.
1482!--    The fluxes flux_d and diss_d at the surface are 0. Due to reasons of
1483!--    statistical evaluation the top flux at the surface should be 0
1484       flux_t(nzb_w_inner(j,i)) = 0.  !statistical reasons
1485       diss_t(nzb_w_inner(j,i)) = 0.
1486!
1487!--    2nd order scheme       
1488       k         = nzb_w_inner(j,i) + 1
1489       flux_d    = flux_t(k-1)
1490       diss_d    = diss_t(k-1)
1491       w_comp    = w(k+1,j,i) + w(k,j,i)
1492       flux_t(k) = w_comp * ( w(k+1,j,i) + w(k,j,i) ) * 0.25
1493       diss_t(k) = diss_2nd( w(k+2,j,i), w(k+1,j,i), w(k,j,i), 0., 0.,        &
1494                             w_comp, 0.25, ddzu(k+1) )
1495                     
1496       tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)                    &
1497                                 - ( flux_d + diss_d )     ) * ddzu(k+1)               
1498!
1499!--    WS3 as an intermediate step
1500       k         = nzb_w_inner(j,i) + 2
1501       flux_d    = flux_t(k-1)
1502       diss_d    = diss_t(k-1)
1503       w_comp    = w(k+1,j,i) + w(k,j,i)
1504       flux_t(k) = w_comp * (                                                 &
1505                   7. * ( w(k+1,j,i) + w(k,j,i) )                             &
1506                   -    ( w(k+2,j,i) + w(k-1,j,i) ) ) * adv_mom_3
1507       diss_t(k) = - abs(w_comp) * (                                          &
1508                   3. * ( w(k+1,j,i) - w(k,j,i) )                             & 
1509                   -    ( w(k+2,j,i) - w(k-1,j,i) ) ) * adv_mom_3
1510
1511       tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)                    &
1512                                 - ( flux_d + diss_d )     ) * ddzu(k+1)
1513!                                 
1514!--    WS5
1515       DO  k = nzb_w_inner(j,i) + 3, nzt -2 
1516          flux_d    = flux_t(k-1)
1517          diss_d    = diss_t(k-1)
1518          w_comp    = w(k+1,j,i) + w(k,j,i)
1519          flux_t(k) = w_comp * (                                              &
1520                      37.  * ( w(k+1,j,i) + w(k,j,i)   )                      &
1521                      - 8. * ( w(k+2,j,i) + w(k-1,j,i) )                      &
1522                      +      ( w(k+3,j,i) + w(k-2,j,i) ) ) * adv_mom_5
1523          diss_t(k) = - abs(w_comp) * (                                       &
1524                      10.  * ( w(k+1,j,i) - w(k,j,i)   )                      & 
1525                      - 5. * ( w(k+2,j,i) - w(k-1,j,i) )                      & 
1526                      +      ( w(k+3,j,i) - w(k-2,j,i) ) ) * adv_mom_5
1527
1528          tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)                 &
1529                                 - ( flux_d + diss_d )     ) * ddzu(k+1)
1530       ENDDO
1531!--    WS3 as an intermediate step
1532       k         = nzt - 1
1533       flux_d    = flux_t(k-1)
1534       diss_d    = diss_t(k-1)
1535       w_comp    = w(k+1,j,i) + w(k,j,i)
1536       flux_t(k) = w_comp * (                                                 &
1537                   7. * ( w(k+1,j,i) + w(k,j,i)   )                           &
1538                   -    ( w(k+2,j,i) + w(k-1,j,i) ) ) *adv_mom_3
1539       diss_t(k) = - abs(w_comp) * (                                          &
1540                   3. * ( w(k+1,j,i) - w(k,j,i)     )                         & 
1541                   -    ( w(k+2,j,i) - w(k-1,j,i) ) ) * adv_mom_3
1542                   
1543       tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)                    &
1544                                 - ( flux_d + diss_d )     ) * ddzu(k+1)
1545!
1546!--    2nd order scheme
1547       k         = nzt
1548       flux_d    = flux_t(k-1)
1549       diss_d    = diss_t(k-1)
1550       w_comp    = w(k+1,j,i) + w(k,j,i)
1551       flux_t(k) = w_comp * ( w(k+1,j,i) + w(k,j,i) ) * 0.25
1552       diss_t(k) = diss_2nd( w(k+1,j,i), w(k+1,j,i), w(k,j,i), w(k-1,j,i),    &
1553                             w(k-2,j,i), w_comp, 0.25, ddzu(k+1) )
1554
1555       tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)                    &
1556                                 - ( flux_d + diss_d )     ) * ddzu(k+1)
1557       
1558       DO  k = nzb_w_inner(j,i), nzt
1559           sums_ws2_ws_l(k,:)  = sums_ws2_ws_l(k,:)                           &
1560               + ( flux_t(k) + diss_t(k) )                                    &
1561               * weight_substep(intermediate_timestep_count) * rmask(j,i,:)
1562       ENDDO
1563
1564    END SUBROUTINE advec_w_ws_ij
1565   
1566
1567!
1568! Scalar advection - Call for all grid points
1569!------------------------------------------------------------------------------!
1570    SUBROUTINE advec_s_ws( sk, sk_char )
1571
1572       USE arrays_3d
1573       USE constants
1574       USE control_parameters
1575       USE grid_variables
1576       USE indices
1577       USE statistics
1578
1579       IMPLICIT NONE
1580
1581       INTEGER ::  i, j, k
1582
1583       REAL, DIMENSION(:,:,:), POINTER ::  sk
1584       REAL    :: flux_d, diss_d, u_comp, v_comp
1585       REAL, DIMENSION(nzb:nzt+1)  ::  flux_r, diss_r, flux_n, diss_n
1586       REAL, DIMENSION(nzb+1:nzt) :: swap_flux_y_local, swap_diss_y_local,    &
1587                                     flux_t, diss_t
1588       REAL, DIMENSION(nzb+1:nzt,nys:nyn) :: swap_flux_x_local,               &
1589                                             swap_diss_x_local
1590       CHARACTER (LEN = *), INTENT(IN) :: sk_char
1591
1592!
1593!--    Compute the fluxes for the whole left boundary of the processor domain.
1594       i = nxl
1595       DO  j = nys, nyn
1596          IF ( boundary_flags(j,i) == 6 )  THEN
1597         
1598             DO  k = nzb_s_inner(j,i) + 1, nzt
1599                u_comp                 = u(k,j,i) - u_gtrans
1600                swap_flux_x_local(k,j) = u_comp * (                           &
1601                                         sk(k,j,i) + sk(k,j,i-1)) * 0.5
1602                swap_diss_x_local(k,j) = diss_2nd( sk(k,j,i+2), sk(k,j,i+1),  &
1603                                                   sk(k,j,i), sk(k,j,i-1),    &
1604                                                   sk(k,j,i-1), u_comp,       &
1605                                                   0.5, ddx ) 
1606             ENDDO
1607             
1608          ELSE
1609         
1610             DO  k = nzb_s_inner(j,i) + 1, nzt
1611                u_comp                 = u(k,j,i) - u_gtrans
1612                swap_flux_x_local(k,j) = u_comp*(                             &
1613                                         37.  * (sk(k,j,i)+sk(k,j,i-1)    )   &
1614                                         -8. * ( sk(k,j,i+1) + sk(k,j,i-2) )  &
1615                                         +     ( sk(k,j,i+2) + sk(k,j,i-3) ) )&
1616                                         * adv_sca_5
1617                swap_diss_x_local(k,j) = - abs(u_comp) * (                    &
1618                                         10. * (sk(k,j,i) - sk(k,j,i-1)    )  &
1619                                         -5. * ( sk(k,j,i+1) - sk(k,j,i-2) )  &
1620                                         +     ( sk(k,j,i+2) - sk(k,j,i-3) ) )&
1621                                         * adv_sca_5
1622             ENDDO
1623          ENDIF
1624       ENDDO
1625!
1626!--    The following loop computes the horizontal fluxes for the interior of the
1627!--    processor domain plus south boundary points. Furthermore tendency terms
1628!--    are computed.
1629       DO  i = nxl, nxr
1630          j = nys
1631          IF ( boundary_flags(j,i) == 8 )  THEN
1632         
1633             DO  k = nzb_s_inner(j,i) + 1, nzt
1634                v_comp               = v(k,j,i) - v_gtrans
1635                swap_flux_y_local(k) = v_comp *                               &
1636                                       ( sk(k,j,i) + sk(k,j-1,i) ) * 0.5 
1637                swap_diss_y_local(k) = diss_2nd( sk(k,j+2,i), sk(k,j+1,i),    &
1638                                                 sk(k,j,i), sk(k,j-1,i),      &
1639                                                 sk(k,j-1,i), v_comp, 0.5, ddy)
1640             ENDDO
1641             
1642          ELSE
1643         
1644             DO  k = nzb_s_inner(j,i) + 1, nzt
1645                v_comp               = v(k,j,i) - v_gtrans
1646                swap_flux_y_local(k) = v_comp * (                             &
1647                                       37.  * ( sk(k,j,i) + sk(k,j-1,i)   )   &
1648                                       - 8. * ( sk(k,j+1,i) + sk(k,j-2,i) )   &
1649                                       +      ( sk(k,j+2,i) + sk(k,j-3,i) ) ) & 
1650                                       * adv_sca_5
1651                swap_diss_y_local(k)= - abs(v_comp) * (                       &
1652                                      10.  * ( sk(k,j,i) - sk(k,j-1,i)   )    &
1653                                      - 5. * ( sk(k,j+1,i) - sk(k,j-2,i) )    &
1654                                      +      ( sk(k,j+2,i)-sk(k,j-3,i) ) )    &
1655                                      * adv_sca_5
1656             ENDDO
1657             
1658          ENDIF
1659
1660          DO  j = nys, nyn
1661            IF ( boundary_flags(j,i) /= 0 )  THEN 
1662!
1663!--            Degrade the order for Dirichlet bc. at the outflow boundary
1664               SELECT CASE ( boundary_flags(j,i) )
1665               
1666                  CASE ( 1 )
1667                     DO  k = nzb_s_inner(j,i) + 1, nzt
1668                        u_comp    = u(k,j,i+1) - u_gtrans
1669                        flux_r(k) = u_comp * (                                &
1670                                  7. * ( sk(k,j,i+1) + sk(k,j,i) )            &
1671                                  - ( sk(k,j,i+2) + sk(k,j,i-1) ) ) * adv_sca_3
1672                        diss_r(k) = - abs(u_comp) * (                         &
1673                                  3. * ( sk(k,j,i+1) - sk(k,j,i) )            & 
1674                                  - ( sk(k,j,i+2) - sk(k,j,i-1) ) ) * adv_sca_3
1675                     ENDDO
1676                     
1677                  CASE ( 2 )
1678                     DO  k = nzb_s_inner(j,i) + 1, nzt
1679                        u_comp    = u(k,j,i+1) - u_gtrans
1680                        flux_r(k) = u_comp * ( sk(k,j,i+1) + sk(k,j,i) ) * 0.5
1681                        diss_r(k) = diss_2nd( sk(k,j,i+1), sk(k,j,i+1),       &
1682                                              sk(k,j,i), sk(k,j,i-1),         &
1683                                              sk(k,j,i-2), u_comp, 0.5, ddx )
1684                     ENDDO
1685                     
1686                  CASE ( 3 )
1687                     DO  k = nzb_s_inner(j,i) + 1, nzt
1688                        v_comp    = v(k,j+1,i) - v_gtrans
1689                        flux_n(k) = v_comp * (                                &
1690                                  7. * ( sk(k,j+1,i) + sk(k,j,i) )            &
1691                                  - ( sk(k,j+2,i) + sk(k,j-1,i) ) ) * adv_sca_3
1692                        diss_n(k) = - abs(v_comp) * (                         &
1693                                  3. * ( sk(k,j+1,i) - sk(k,j,i) )            & 
1694                                  - ( sk(k,j+2,i) - sk(k,j-1,i) ) ) * adv_sca_3
1695                     ENDDO
1696                     
1697                  CASE ( 4 )
1698                     DO  k = nzb_s_inner(j,i) + 1, nzt
1699                        v_comp    = v(k,j+1,i) - v_gtrans
1700                        flux_n(k) = v_comp * ( sk(k,j+1,i) + sk(k,j,i) ) * 0.5
1701                        diss_n(k) = diss_2nd( sk(k,j+1,i), sk(k,j+1,i),       &
1702                                              sk(k,j,i), sk(k,j-1,i),         &
1703                                              sk(k,j-2,i), v_comp, 0.5, ddy )     
1704                     ENDDO
1705                     
1706                  CASE ( 5 )
1707                     DO  k = nzb_w_inner(j,i) + 1, nzt
1708                        u_comp    = u(k,j,i+1) - u_gtrans
1709                        flux_r(k) = u_comp * (                                &
1710                                  7. * ( sk(k,j,i+1) + sk(k,j,i) )            &
1711                                  - ( sk(k,j,i+2) + sk(k,j,i-1) ) ) * adv_sca_3
1712                        diss_r(k) = - abs(u_comp) * (                         &
1713                                  3. * ( sk(k,j,i+1) - sk(k,j,i) )            & 
1714                                  - ( sk(k,j,i+2) - sk(k,j,i-1) ) ) * adv_sca_3
1715                     ENDDO 
1716                     
1717                  CASE ( 6 )
1718                     DO  k = nzb_s_inner(j,i) + 1, nzt
1719                        u_comp    = u(k,j,i+1) - u_gtrans
1720                        flux_r(k) = u_comp * (                                &
1721                                  7. * ( sk(k,j,i+1) + sk(k,j,i) )            &
1722                                  - ( sk(k,j,i+2) + sk(k,j,i-1) ) ) * adv_sca_3
1723                        diss_r(k) = - abs(u_comp) * (                         &
1724                                  3. * ( sk(k,j,i+1) - sk(k,j,i) )            & 
1725                                  - ( sk(k,j,i+2) - sk(k,j,i-1) ) ) * adv_sca_3
1726                     ENDDO
1727                     
1728                  CASE ( 7 )
1729                     DO  k = nzb_s_inner(j,i) + 1, nzt
1730                        v_comp    = v(k,j+1,i) - v_gtrans
1731                        flux_n(k) = v_comp * (                                &
1732                                  7. * ( sk(k,j+1,i) + sk(k,j,i) )            &
1733                                  - ( sk(k,j+2,i) + sk(k,j-1,i) ) ) * adv_sca_3
1734                        diss_n(k) = - abs(v_comp) * (                         &
1735                                  3. * ( sk(k,j+1,i) - sk(k,j,i) )            & 
1736                                  - ( sk(k,j+2,i) - sk(k,j-1,i) ) ) * adv_sca_3
1737                     ENDDO
1738                     
1739                  CASE ( 8 )
1740                     DO  k = nzb_s_inner(j,i) + 1, nzt
1741                        v_comp    = v(k,j+1,i) - v_gtrans
1742                        flux_n(k) = v_comp * (                                &
1743                                  7. * ( sk(k,j+1,i) + sk(k,j,i) )            &
1744                                  - ( sk(k,j+2,i) + sk(k,j-1,i) ) ) * adv_sca_3
1745                        diss_n(k) = -  abs(v_comp) * (                        &
1746                                  3. * ( sk(k,j+1,i) - sk(k,j,i) )            & 
1747                                  - ( sk(k,j+2,i) - sk(k,j-1,i) ) ) * adv_sca_3
1748                     ENDDO 
1749                     
1750                  CASE DEFAULT
1751                 
1752               END SELECT
1753               
1754               IF (boundary_flags(j,i) == 1 .OR. boundary_flags(j,i) == 2     &
1755             .OR. boundary_flags(j,i) == 5 .OR. boundary_flags(j,i) == 6)  THEN
1756             
1757                  DO  k = nzb_s_inner(j,i) + 1, nzt
1758                     v_comp    = v(k,j+1,i) - v_gtrans
1759                     flux_n(k) = v_comp * (                                   &
1760                                 37.  * ( sk(k,j+1,i) + sk(k,j,i)  )          &
1761                                 - 8. * (sk(k,j+2,i) + sk(k,j-1,i) )          &
1762                                 +      ( sk(k,j+3,i) + sk(k,j-2,i) ) )       &
1763                                 * adv_sca_5
1764                     diss_n(k) = - abs(v_comp) * (                            &
1765                                 10.  * ( sk(k,j+1,i) - sk(k,j,i)   )         &
1766                                 - 5. * ( sk(k,j+2,i) - sk(k,j-1,i) )         &
1767                                 +      ( sk(k,j+3,i) - sk(k,j-2,i) ) )       &
1768                                 * adv_sca_5
1769                  ENDDO
1770                 
1771               ELSE
1772               
1773                  DO  k = nzb_s_inner(j,i) + 1, nzt
1774                     u_comp    = u(k,j,i+1) - u_gtrans 
1775                     flux_r(k) = u_comp * (                                   & 
1776                                 37.  * ( sk(k,j,i+1) + sk(k,j,i)   )         &
1777                                 - 8. * ( sk(k,j,i+2) + sk(k,j,i-1) )         &
1778                                 +      ( sk(k,j,i+3) + sk(k,j,i-2) ) )       &
1779                                 * adv_sca_5
1780                     diss_r(k) = - abs(u_comp) * (                            &
1781                                 10.  * ( sk(k,j,i+1) - sk(k,j,i)   )         &
1782                                 - 5. * ( sk(k,j,i+2) - sk(k,j,i-1) )         &
1783                                 +      ( sk(k,j,i+3) - sk(k,j,i-2) ) )       &
1784                                 * adv_sca_5
1785                  ENDDO
1786                 
1787               ENDIF     
1788           
1789            ELSE
1790           
1791               DO  k = nzb_s_inner(j,i) + 1, nzt
1792                  u_comp    = u(k,j,i+1) - u_gtrans
1793                  flux_r(k) = u_comp * (                                      &
1794                              37.  * ( sk(k,j,i+1) + sk(k,j,i)   )            &
1795                              - 8. * ( sk(k,j,i+2) + sk(k,j,i-1) )            &
1796                              +      ( sk(k,j,i+3) + sk(k,j,i-2) ) )          &
1797                              * adv_sca_5
1798                  diss_r(k) = - abs(u_comp) * (                               &
1799                              10.  * ( sk(k,j,i+1) - sk(k,j,i)   )            &
1800                              - 5. * ( sk(k,j,i+2) - sk(k,j,i-1) )            &
1801                              +      ( sk(k,j,i+3) - sk(k,j,i-2) ) )          &
1802                              * adv_sca_5
1803                 
1804                  v_comp    = v(k,j+1,i) - v_gtrans
1805                  flux_n(k) = v_comp * (                                      &
1806                              37.  * ( sk(k,j+1,i) + sk(k,j,i)   )            &
1807                              - 8. * ( sk(k,j+2,i) + sk(k,j-1,i) )            &
1808                              +      ( sk(k,j+3,i) + sk(k,j-2,i) ) )          &
1809                              * adv_sca_5
1810                  diss_n(k) = - abs(v_comp) * (                               &
1811                              10.  * ( sk(k,j+1,i) - sk(k,j,i)   )            &
1812                              - 5. * ( sk(k,j+2,i) - sk(k,j-1,i) )            &
1813                              +      ( sk(k,j+3,i) - sk(k,j-2,i) ) )          &
1814                              * adv_sca_5                 
1815               ENDDO
1816               
1817            ENDIF
1818           
1819            DO  k = nzb_s_inner(j,i) + 1, nzt
1820               tend(k,j,i) = tend(k,j,i) - (                                  &
1821                  ( flux_r(k) + diss_r(k)                                     &
1822                - ( swap_flux_x_local(k,j) + swap_diss_x_local(k,j) ) ) * ddx &
1823                + ( flux_n(k) + diss_n(k)                                     &
1824                - ( swap_flux_y_local(k) + swap_diss_y_local(k) )     ) * ddy)
1825                   
1826                swap_flux_x_local(k,j) = flux_r(k)
1827                swap_diss_x_local(k,j) = diss_r(k)
1828                swap_flux_y_local(k)   = flux_n(k)
1829                swap_diss_y_local(k)   = diss_n(k)
1830            ENDDO
1831         ENDDO
1832      ENDDO
1833
1834!
1835!--   Vertical advection, degradation of order near surface and top.
1836!--   The fluxes flux_d and diss_d at the surface are 0. Due to reasons of
1837!--   statistical evaluation the top flux at the surface should be 0
1838      DO  i = nxl, nxr
1839         DO  j = nys, nyn
1840!
1841!--         2nd order scheme
1842            k=nzb_s_inner(j,i) + 1
1843!           
1844!--         The fluxes flux_d and diss_d at the surface are 0. Due to static
1845!--         reasons the top flux at the surface should be 0.
1846            flux_t(nzb_s_inner(j,i)) = 0.
1847            diss_t(nzb_s_inner(j,i)) = 0.
1848            flux_d = flux_t(k-1)
1849            diss_d = diss_t(k-1)
1850            flux_t(k) = w(k,j,i) * (sk(k+1,j,i) + sk(k,j,i) ) *0.5
1851            diss_t(k) = diss_2nd( sk(k+2,j,i), sk(k+1,j,i), sk(k,j,i),        &
1852                                  sk(k,j,i), sk(k,j,i), w(k,j,i),             &
1853                                  0.5, ddzw(k) )
1854            tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)               &
1855                                      - ( flux_d + diss_d )     ) * ddzw(k)   
1856!
1857!--         WS3 as an intermediate step
1858            k         = nzb_s_inner(j,i) + 2
1859            flux_d    = flux_t(k-1)
1860            diss_d    = diss_t(k-1)
1861            flux_t(k) = w(k,j,i) * (                                          &
1862                        7. * ( sk(k+1,j,i) + sk(k,j,i) )                      &
1863                       - ( sk(k+2,j,i) + sk(k-1,j,i) ) ) * adv_sca_3
1864            diss_t(k) = - abs(w(k,j,i)) * (                                   &
1865                        3. * ( sk(k+1,j,i) - sk(k,j,i) )                      & 
1866                        - ( sk(k+2,j,i) - sk(k-1,j,i) ) ) * adv_sca_3
1867            tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)               &
1868                                      - ( flux_d + diss_d )     ) * ddzw(k)
1869!
1870!--         WS5
1871            DO  k = nzb_s_inner(j,i) + 3, nzt - 2
1872               flux_d    = flux_t(k-1)
1873               diss_d    = diss_t(k-1)
1874               flux_t(k) = w(k,j,i) * (                                       &
1875                           37.  * ( sk(k+1,j,i) + sk(k,j,i)   )               &
1876                           - 8. * ( sk(k+2,j,i) + sk(k-1,j,i) )               &
1877                           +      ( sk(k+3,j,i) + sk(k-2,j,i) ) ) * adv_sca_5
1878               diss_t(k) = - abs(w(k,j,i)) * (                                &
1879                           10.  * ( sk(k+1,j,i) -sk(k,j,i) )                  &
1880                           - 5. * ( sk(k+2,j,i) - sk(k-1,j,i) )               &
1881                           +      ( sk(k+3,j,i) - sk(k-2,j,i) ) ) * adv_sca_5
1882
1883               tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)            &
1884                                      - ( flux_d + diss_d )     ) * ddzw(k)
1885            ENDDO
1886!
1887!--         WS3 as an intermediate step
1888            k         = nzt - 1
1889            flux_d    = flux_t(k-1)
1890            diss_d    = diss_t(k-1)
1891            flux_t(k) = w(k,j,i) * (                                          &
1892                        7. * ( sk(k+1,j,i) + sk(k,j,i)  )                     &
1893                        - ( sk(k+2,j,i) + sk(k-1,j,i) ) ) * adv_sca_3
1894            diss_t(k) = - abs(w(k,j,i)) * (                                   &
1895                        3. * ( sk(k+1,j,i) - sk(k,j,i)   )                    & 
1896                        -    ( sk(k+2,j,i) - sk(k-1,j,i) ) ) * adv_sca_3
1897                       
1898            tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)               &
1899                                      - ( flux_d + diss_d )     ) * ddzw(k)
1900!
1901!--         2nd order scheme
1902            k         = nzt
1903            flux_d    = flux_t(k-1)
1904            diss_d    = diss_t(k-1)
1905            flux_t(k) = w(k,j,i) * ( sk(k+1,j,i) + sk(k,j,i) ) * 0.5
1906            diss_t(k) = diss_2nd( sk(k+1,j,i), sk(k+1,j,i), sk(k,j,i),        &
1907                                  sk(k-1,j,i), sk(k-2,j,i), w(k,j,i),         &
1908                                  0.5, ddzw(k) )
1909
1910            tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)               &
1911                                      - ( flux_d + diss_d )     ) * ddzw(k)
1912!
1913!--         evaluation of statistics
1914            SELECT CASE ( sk_char )
1915
1916               CASE ( 'pt' )
1917                 DO  k = nzb_s_inner(j,i), nzt
1918                   sums_wspts_ws_l(k,:) = sums_wspts_ws_l(k,:)                 &
1919                      + ( flux_t(k) + diss_t(k) )                              &
1920                      * weight_substep(intermediate_timestep_count)            &
1921                      * rmask(j,i,:)
1922                 ENDDO
1923               CASE ( 'sa' )
1924                 DO  k = nzb_s_inner(j,i), nzt
1925                   sums_wssas_ws_l(k,:) = sums_wssas_ws_l(k,:)                 &
1926                      + ( flux_t(k) + diss_t(k) )                              &
1927                      * weight_substep(intermediate_timestep_count)            &
1928                      * rmask(j,i,:)
1929                 ENDDO
1930               CASE ( 'q' )
1931                 DO  k = nzb_s_inner(j,i), nzt
1932                   sums_wsqs_ws_l(k,:) = sums_wsqs_ws_l(k,:)                   &
1933                      + ( flux_t(k) + diss_t(k) )                              &
1934                      * weight_substep(intermediate_timestep_count)            &
1935                      * rmask(j,i,:)
1936                 ENDDO
1937
1938            END SELECT
1939         ENDDO
1940      ENDDO
1941
1942
1943    END SUBROUTINE advec_s_ws
1944
1945
1946!
1947! Advection of u - Call for all grid points
1948!------------------------------------------------------------------------------!
1949    SUBROUTINE advec_u_ws
1950
1951       USE arrays_3d
1952       USE constants
1953       USE control_parameters
1954       USE grid_variables
1955       USE indices
1956       USE statistics
1957
1958       IMPLICIT NONE
1959
1960       INTEGER ::  i, j, k
1961       REAL    ::  gu, gv, flux_d, diss_d, v_comp, w_comp
1962       REAL, DIMENSION(nzb+1:nzt) :: swap_flux_y_local_u, swap_diss_y_local_u
1963       REAL, DIMENSION(nzb+1:nzt,nys:nyn) :: swap_flux_x_local_u,             &
1964                                             swap_diss_x_local_u
1965       REAL, DIMENSION(nzb:nzt+1) :: flux_t, diss_t, flux_r, diss_r, flux_n,  &
1966                                     diss_n, u_comp
1967 
1968       gu = 2.0 * u_gtrans
1969       gv = 2.0 * v_gtrans
1970
1971!
1972!--    Compute the fluxes for the whole left boundary of the processor domain.
1973       i = nxlu
1974       DO  j = nys, nyn
1975          IF( boundary_flags(j,i) == 5 )  THEN
1976         
1977             DO  k = nzb_u_inner(j,i) + 1, nzt
1978                u_comp(k)                = u(k,j,i) + u(k,j,i-1) - gu
1979                swap_flux_x_local_u(k,j) = u_comp(k) *                        &
1980                                           ( u(k,j,i) + u(k,j,i-1) ) * 0.25
1981                swap_diss_x_local_u(k,j) = diss_2nd( u(k,j,i+2), u(k,j,i+1),  &
1982                                                     u(k,j,i), u(k,j,i-1),    &
1983                                                     u(k,j,i-1), u_comp(k),   &
1984                                                     0.25, ddx )
1985             ENDDO
1986             
1987          ELSE
1988         
1989             DO  k = nzb_u_inner(j,i) + 1, nzt
1990                u_comp(k)                = u(k,j,i) + u(k,j,i-1) - gu
1991                swap_flux_x_local_u(k,j) = u_comp(k) * (                      &
1992                                         37.  * ( u(k,j,i) + u(k,j,i-1)   )   &
1993                                         - 8. * ( u(k,j,i+1) + u(k,j,i-2) )   &
1994                                         +      (u(k,j,i+2)+u(k,j,i-3) )  )   &
1995                                         * adv_mom_5
1996                swap_diss_x_local_u(k,j) = - abs(u_comp(k)) * (               &
1997                                         10.  * ( u(k,j,i) - u(k,j,i-1)   )   &
1998                                         - 5. * ( u(k,j,i+1) - u(k,j,i-2) )   &
1999                                         +      ( u(k,j,i+2) - u(k,j,i-3) ) ) &
2000                                         * adv_mom_5
2001             ENDDO
2002             
2003          ENDIF
2004         
2005       ENDDO
2006
2007       DO i = nxlu, nxr
2008!       
2009!--      The following loop computes the fluxes for the south boundary points
2010         j = nys
2011         IF ( boundary_flags(j,i) == 8 )  THEN
2012!         
2013!--         Compute southside fluxes for the south boundary of PE domain
2014            DO  k = nzb_u_inner(j,i) + 1, nzt
2015               v_comp                 = v(k,j,i) + v(k,j,i-1) - gv
2016               swap_flux_y_local_u(k) = v_comp *                              &
2017                                        ( u(k,j,i) + u(k,j-1,i) ) * 0.25 
2018               swap_diss_y_local_u(k) = diss_2nd( u(k,j+2,i), u(k,j+1,i),     &
2019                                                  u(k,j,i), u(k,j-1,i),       &
2020                                                  u(k,j-1,i), v_comp,         &
2021                                                  0.25, ddy )
2022            ENDDO
2023           
2024         ELSE
2025         
2026            DO  k = nzb_u_inner(j,i) + 1, nzt
2027               v_comp                 = v(k,j,i) + v(k,j,i-1) - gv
2028               swap_flux_y_local_u(k) = v_comp * (                            &
2029                                      37.  * ( u(k,j,i) + u(k,j-1,i)   )      &
2030                                      - 8. * ( u(k,j+1,i) + u(k,j-2,i) )      &
2031                                      +      ( u(k,j+2,i) + u(k,j-3,i) ) )    &
2032                                      * adv_mom_5
2033               swap_diss_y_local_u(k) = - abs(v_comp) * (                     &
2034                                      10.  * ( u(k,j,i) - u(k,j-1,i)   )      &
2035                                      - 5. * ( u(k,j+1,i) - u(k,j-2,i) )      &
2036                                      +      ( u(k,j+2,i) - u(k,j-3,i) ) )    &
2037                                      * adv_mom_5
2038            ENDDO
2039           
2040         ENDIF
2041!         
2042!--      Computation of interior fluxes and tendency terms
2043         DO  j = nys, nyn
2044            IF ( boundary_flags(j,i) /= 0 )  THEN
2045!           
2046!--            Degrade the order for Dirichlet bc. at the outflow boundary
2047               SELECT CASE ( boundary_flags(j,i) )
2048               
2049                  CASE ( 1 )
2050                     DO  k = nzb_u_inner(j,i) + 1, nzt
2051                        u_comp(k) = u(k,j,i+1) + u(k,j,i)
2052                        flux_r(k) = ( u_comp(k) - gu ) * (                    &
2053                                    7. * ( u(k,j,i+1) + u(k,j,i) )            &
2054                                    - ( u(k,j,i+2) + u(k,j,i-1) ) ) * adv_mom_3
2055                        diss_r(k) = - abs(u_comp(k) - gu) * (                 &
2056                                    3. * ( u(k,j,i+1) - u(k,j,i) )            & 
2057                                    - ( u(k,j,i+2) - u(k,j,i-1) ) ) * adv_mom_3
2058                     ENDDO
2059                     
2060                  CASE ( 2 )
2061                     DO  k = nzb_u_inner(j,i) + 1, nzt
2062                        u_comp(k) = u(k,j,i+1) + u(k,j,i)
2063                        flux_r(k) = ( u_comp(k) - gu ) *                      &
2064                                    ( u(k,j,i+1) + u(k,j,i) ) * 0.25 
2065                        diss_r(k) = diss_2nd( u(k,j,i+1), u(k,j,i+1),         &
2066                                              u(k,j,i), u(k,j,i-1),           &
2067                                              u(k,j,i-2), u_comp(k) ,0.25 ,ddx)
2068                     ENDDO
2069                     
2070                  CASE ( 3 )
2071                     DO  k = nzb_u_inner(j,i) + 1, nzt
2072                        v_comp   = v(k,j+1,i) + v(k,j+1,i-1) - gv
2073                        flux_n(k) = v_comp * (                                &
2074                                    7. * ( u(k,j+1,i) + u(k,j,i) )            &
2075                                    - ( u(k,j+2,i) + u(k,j-1,i) ) ) * adv_mom_3
2076                        diss_n(k) = - abs(v_comp) * (                         &
2077                                    3. * ( u(k,j+1,i) - u(k,j,i) )            & 
2078                                    - ( u(k,j+2,i) - u(k,j-1,i) ) ) * adv_mom_3
2079                     ENDDO
2080                     
2081                  CASE ( 4 )
2082                     DO  k = nzb_u_inner(j,i) + 1, nzt
2083                        v_comp    = v(k,j+1,i) + v(k,j+1,i-1) - gv
2084                        flux_n(k) = v_comp * ( u(k,j+1,i) + u(k,j,i) ) * 0.25
2085                        diss_n(k) = diss_2nd( u(k,j+1,i), u(k,j+1,i),         &
2086                                              u(k,j,i), u(k,j-1,i),           &
2087                                              u(k,j-2,i), v_comp, 0.25, ddy )
2088                     ENDDO
2089                     
2090                  CASE ( 5 )
2091                     DO  k = nzb_u_inner(j,i) + 1, nzt       
2092                        u_comp(k) = u(k,j,i+1) + u(k,j,i)
2093                        flux_r(k) = ( u_comp(k) - gu ) * (                    &
2094                                    7. * ( u(k,j,i+1) + u(k,j,i) )            &
2095                                    - ( u(k,j,i+2) + u(k,j,i-1) ) ) * adv_mom_3
2096                        diss_r(k) = - abs(u_comp(k) - gu) * (                 &
2097                                    3. * ( u(k,j,i+1) - u(k,j,i) )            & 
2098                                    - ( u(k,j,i+2) - u(k,j,i-1) ) ) * adv_mom_3
2099                     ENDDO
2100                     
2101                  CASE ( 7 )
2102                     DO  k = nzb_u_inner(j,i) + 1, nzt                           
2103                        v_comp    = v(k,j+1,i) + v(k,j+1,i-1) - gv
2104                        flux_n(k) = v_comp * (                                &
2105                                    7. * ( u(k,j+1,i) + u(k,j,i) )            &
2106                                    - ( u(k,j+2,i) + u(k,j-1,i) ) ) * adv_mom_3
2107                        diss_n(k) = - abs(v_comp) * (                         &
2108                                    3. * ( u(k,j+1,i) - u(k,j,i) )            & 
2109                                    - ( u(k,j+2,i) - u(k,j-1,i) ) ) * adv_mom_3
2110                     ENDDO
2111                     
2112                  CASE ( 8 )
2113                     DO  k = nzb_u_inner(j,i) + 1, nzt
2114                        v_comp    = v(k,j+1,i) + v(k,j+1,i-1) - gv
2115                        flux_n(k) = v_comp * (                                &
2116                                    7. * ( u(k,j+1,i) + u(k,j,i) )            &
2117                                    - ( u(k,j+2,i) + u(k,j-1,i) ) ) * adv_mom_3
2118                        diss_n(k) = - abs(v_comp) * (                         &
2119                                    3. * ( u(k,j+1,i) - u(k,j,i) )            & 
2120                                    - ( u(k,j+2,i) - u(k,j-1,i) ) ) * adv_mom_3
2121                     ENDDO 
2122                     
2123                  CASE DEFAULT
2124                 
2125               END SELECT
2126!               
2127!--            Compute the crosswise 5th order fluxes at the outflow
2128               IF (boundary_flags(j,i) == 1 .OR. boundary_flags(j,i) == 2     &
2129              .OR. boundary_flags(j,i) == 5)  THEN
2130             
2131                  DO  k = nzb_u_inner(j,i) + 1, nzt
2132                     v_comp = v(k,j+1,i) + v(k,j+1,i-1) - gv
2133                     flux_n(k) = v_comp * (                                   &
2134                               37.  * ( u(k,j+1,i) + u(k,j,i)   )             &
2135                               - 8. * ( u(k,j+2,i) +u(k,j-1,i)  )             &
2136                               +      ( u(k,j+3,i) + u(k,j-2,i) ) ) * adv_mom_5
2137                     diss_n(k) = - abs(v_comp) * (                            &
2138                               10.  * ( u(k,j+1,i) - u(k,j,i) )               &
2139                               - 5. * ( u(k,j+2,i) - u(k,j-1,i) )             &
2140                               +      ( u(k,j+3,i) - u(k,j-2,i) ) ) * adv_mom_5
2141                  ENDDO
2142                 
2143               ELSE
2144               
2145                  DO  k = nzb_u_inner(j,i) + 1, nzt
2146                     u_comp(k) = u(k,j,i+1) + u(k,j,i)
2147                     flux_r(k) = ( u_comp(k) - gu ) * (                       &
2148                               37.  * ( u(k,j,i+1) + u(k,j,i)   )             &
2149                               - 8. * ( u(k,j,i+2) + u(k,j,i-1) )             &
2150                               +      ( u(k,j,i+3) + u(k,j,i-2) ) ) * adv_mom_5
2151                     diss_r(k) = - abs(u_comp(k) - gu) * (                    &
2152                               10.  * ( u(k,j,i+1) - u(k,j,i)   )             &
2153                               - 5. * ( u(k,j,i+2) - u(k,j,i-1) )             &
2154                               +      ( u(k,j,i+3) - u(k,j,i-2) ) ) * adv_mom_5
2155                  ENDDO
2156                 
2157               ENDIF
2158               
2159            ELSE
2160           
2161               DO  k = nzb_u_inner(j,i) + 1, nzt
2162                  u_comp(k) = u(k,j,i+1) + u(k,j,i)
2163                  flux_r(k) = ( u_comp(k) - gu ) * (                          &
2164                            37.  * ( u(k,j,i+1) + u(k,j,i)   )                &
2165                            - 8. * ( u(k,j,i+2) + u(k,j,i-1) )                &
2166                            +      ( u(k,j,i+3) + u(k,j,i-2) ) ) * adv_mom_5
2167                  diss_r(k) = - abs(u_comp(k) - gu) * (                       &
2168                            10.  * ( u(k,j,i+1) - u(k,j,i)   )                &
2169                            - 5. * ( u(k,j,i+2) - u(k,j,i-1) )                &
2170                            +      ( u(k,j,i+3) - u(k,j,i-2) ) ) * adv_mom_5
2171
2172                  v_comp = v(k,j+1,i) + v(k,j+1,i-1) - gv
2173                  flux_n(k) = v_comp * (                                      &
2174                            37.  * ( u(k,j+1,i) + u(k,j,i)   )                &
2175                            - 8. * ( u(k,j+2,i) + u(k,j-1,i) )                &
2176                            +      ( u(k,j+3,i) + u(k,j-2,i) ) ) * adv_mom_5
2177                  diss_n(k) = - abs(v_comp) * (                               &
2178                            10.  * ( u(k,j+1,i) - u(k,j,i)   )                &
2179                            - 5. * ( u(k,j+2,i) - u(k,j-1,i) )                &
2180                            +      ( u(k,j+3,i) - u(k,j-2,i) ) ) * adv_mom_5
2181                               
2182               ENDDO
2183               
2184            ENDIF
2185           
2186            DO  k = nzb_u_inner(j,i) + 1, nzt
2187           
2188               tend(k,j,i) = tend(k,j,i) - (                                  &
2189            ( flux_r(k) + diss_r(k)               &
2190            - ( swap_flux_x_local_u(k,j) + swap_diss_x_local_u(k,j) ) ) * ddx &
2191            + ( flux_n(k) + diss_n(k)               &
2192            - ( swap_flux_y_local_u(k) + swap_diss_y_local_u(k)     ) ) * ddy )
2193                 
2194               swap_flux_x_local_u(k,j) = flux_r(k)   
2195               swap_diss_x_local_u(k,j) = diss_r(k)
2196               swap_flux_y_local_u(k)   = flux_n(k)     
2197               swap_diss_y_local_u(k)   = diss_n(k)
2198                     
2199               sums_us2_ws_l(k,:)  = sums_us2_ws_l(k,:)                       &
2200                 + ( flux_r(k)                                                &
2201                 * ( u_comp(k) - 2. * hom(k,1,1,:) )                          &
2202                 / ( u_comp(k) - gu + 1.0E-20 )                               & 
2203                 + diss_r(k)                                                  &
2204                 * abs(u_comp(k) - 2. * hom(k,1,1,:) )                        &
2205                 / (abs(u_comp(k) - gu) + 1.0E-20) )                          &
2206                 * weight_substep(intermediate_timestep_count) * rmask(j,i,:)
2207            ENDDO
2208         ENDDO
2209       ENDDO
2210
2211!
2212!--   Vertical advection, degradation of order near surface and top.
2213!--   The fluxes flux_d and diss_d at the surface are 0. Due to reasons of
2214!--   statistical evaluation the top flux at the surface should be 0
2215       DO i = nxlu, nxr
2216          DO  j = nys, nyn
2217             k = nzb_u_inner(j,i) + 1
2218!             
2219!--         The fluxes flux_d and diss_d at the surface are 0. Due to static
2220!--         reasons the top flux at the surface should be 0.
2221            flux_t(nzb_u_inner(j,i)) = 0.
2222            diss_t(nzb_u_inner(j,i)) = 0.
2223            flux_d = flux_t(k-1)
2224            diss_d = diss_t(k-1)             
2225!
2226!--         2nd order scheme
2227             w_comp    = w(k,j,i) + w(k,j,i-1)
2228             flux_t(k) = w_comp * ( u(k+1,j,i) + u(k,j,i) ) * 0.25
2229             diss_t(k) = diss_2nd( u(k+2,j,i), u(k+1,j,i), u(k,j,i),          &
2230                                   0., 0., w_comp, 0.25, ddzw(k) )
2231                                   
2232             tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)               &
2233                                      - ( flux_d + diss_d )     ) * ddzw(k)
2234!
2235!--         WS3 as an intermediate step
2236            k         = nzb_u_inner(j,i) + 2
2237            flux_d    = flux_t(k-1)
2238            diss_d    = diss_t(k-1)
2239            w_comp    = w(k,j,i) + w(k,j,i-1)
2240            flux_t(k) = w_comp * (                                            &
2241                        7. * (u(k+1,j,i) + u(k,j,i)    )                      &
2242                        -    ( u(k+2,j,i) + u(k-1,j,i) ) ) * adv_mom_3
2243            diss_t(k) = - abs(w_comp) * (                                     &
2244                        3. * ( u(k+1,j,i) - u(k,j,i) )                        & 
2245                        -    ( u(k+2,j,i) - u(k-1,j,i) ) ) * adv_mom_3
2246
2247            tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)               &
2248                                      - ( flux_d + diss_d )     ) * ddzw(k)
2249!
2250!WS5
2251             DO  k = nzb_u_inner(j,i) + 3, nzt - 3
2252
2253                flux_d    = flux_t(k-1)
2254                diss_d    = diss_t(k-1)
2255                w_comp    = w(k,j,i) + w(k,j,i-1)
2256                flux_t(k) = w_comp * (                                        &
2257                            37.  * ( u(k+1,j,i) + u(k,j,i)   )                &
2258                            - 8. * ( u(k+2,j,i) + u(k-1,j,i) )                &
2259                            +      ( u(k+3,j,i) + u(k-2,j,i) ) ) * adv_mom_5
2260                diss_t(k) = - abs(w_comp) * (                                 &
2261                            10.  * ( u(k+1,j,i) - u(k,j,i)   )                & 
2262                            - 5. * ( u(k+2,j,i) - u(k-1,j,i) )                & 
2263                            +      ( u(k+3,j,i) - u(k-2,j,i) ) ) * adv_mom_5
2264
2265                tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)           &
2266                                      - ( flux_d + diss_d )     ) * ddzw(k)
2267
2268             ENDDO
2269!
2270!--         WS3 as an intermediate step
2271             DO k = nzt - 2, nzt - 1
2272                flux_d = flux_t(k-1)
2273                diss_d = diss_t(k-1)
2274                w_comp = w(k,j,i) + w(k,j,i-1)
2275                flux_t(k) = w_comp * (                                        &
2276                            7. * ( u(k+1,j,i) + u(k,j,i) )                    &
2277                            -    ( u(k+2,j,i) + u(k-1,j,i) ) ) * adv_mom_3
2278                diss_t(k) = - abs(w_comp) * (                                 &
2279                            3. * ( u(k+1,j,i) - u(k,j,i) )                    & 
2280                            -    ( u(k+2,j,i) - u(k-1,j,i) ) ) * adv_mom_3
2281                           
2282                tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)           &
2283                                      - ( flux_d + diss_d )     ) * ddzw(k)
2284             ENDDO
2285!
2286!--         2nd order scheme
2287             k         = nzt
2288             flux_d    = flux_t(k-1)
2289             diss_d    = diss_t(k-1)
2290             w_comp    = w(k,j,i) + w(k,j,i-1)
2291             flux_t(k) = w_comp * ( u(k+1,j,i) + u(k,j,i) ) * 0.25
2292             diss_t(k) = diss_2nd( u(nzt+1,j,i), u(nzt+1,j,i), u(k,j,i),      &
2293                                   u(k-1,j,i), u(k-2,j,i), w_comp,            &
2294                                   0.25, ddzw(k)) 
2295
2296             tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)              &
2297                                      - ( flux_d + diss_d )     ) * ddzw(k)
2298!
2299!-- at last vertical momentum flux is accumulated
2300            DO  k = nzb_u_inner(j,i), nzt
2301               sums_wsus_ws_l(k,:) = sums_wsus_ws_l(k,:)                       &
2302                   + ( flux_t(k) + diss_t(k) )                                 &
2303                   * weight_substep(intermediate_timestep_count)               &
2304                   * rmask(j,i,:)
2305            ENDDO
2306          ENDDO
2307       ENDDO
2308
2309
2310    END SUBROUTINE advec_u_ws
2311   
2312   
2313!
2314! Advection of v - Call for all grid points
2315!------------------------------------------------------------------------------!
2316    SUBROUTINE advec_v_ws
2317
2318       USE arrays_3d
2319       USE constants
2320       USE control_parameters
2321       USE grid_variables
2322       USE indices
2323       USE statistics
2324
2325       IMPLICIT NONE
2326
2327
2328       INTEGER ::  i, j, k
2329       REAL    ::  gu, gv, flux_l, flux_s, flux_d, diss_l, diss_s, diss_d,    &
2330                   u_comp, w_comp
2331       REAL, DIMENSION(nzb+1:nzt) :: swap_flux_y_local_v, swap_diss_y_local_v
2332       REAL, DIMENSION(nzb+1:nzt,nys:nyn) :: swap_flux_x_local_v,             &
2333                                             swap_diss_x_local_v
2334       REAL, DIMENSION(nzb:nzt+1) :: flux_t, diss_t, flux_n, diss_n, flux_r,  &
2335                                     diss_r, v_comp
2336
2337       gu = 2.0 * u_gtrans
2338       gv = 2.0 * v_gtrans
2339!
2340!-- First compute the whole left boundary of the processor domain
2341       i = nxl
2342       DO  j = nysv, nyn
2343       
2344          IF ( boundary_flags(j,i) == 6 )  THEN
2345             DO  k = nzb_v_inner(j,i) + 1, nzt
2346                u_comp                   = u(k,j-1,i) + u(k,j,i) - gu
2347                swap_flux_x_local_v(k,j) = u_comp *                           &
2348                                           ( v(k,j,i) + v(k,j,i-1)) * 0.25
2349                swap_diss_x_local_v(k,j) = diss_2nd( v(k,j,i+2), v(k,j,i+1),  &
2350                                                     v(k,j,i), v(k,j,i-1),    &
2351                                                     v(k,j,i-1), u_comp,      &
2352                                                     0.25, ddx )
2353             ENDDO
2354             
2355          ELSE
2356         
2357             DO  k = nzb_v_inner(j,i)+1, nzt
2358                u_comp                   = u(k,j-1,i) + u(k,j,i) - gu
2359                swap_flux_x_local_v(k,j) = u_comp * (                         &
2360                                         37.  * ( v(k,j,i) + v(k,j,i-1)   )   &
2361                                         - 8. * ( v(k,j,i+1) + v(k,j,i-2) )   &
2362                                         +      ( v(k,j,i+2) + v(k,j,i-3) ) ) &
2363                                         * adv_mom_5
2364                swap_diss_x_local_v(k,j) = - abs(u_comp) * (                  &
2365                                         10. * ( v(k,j,i) - v(k,j,i-1)   )    &
2366                                         -5. * ( v(k,j,i+1) - v(k,j,i-2) )    &
2367                                         +     ( v(k,j,i+2) - v(k,j,i-3) ) )  &
2368                                         * adv_mom_5
2369             ENDDO
2370             
2371          ENDIF
2372         
2373       ENDDO
2374
2375       DO i = nxl, nxr
2376         j = nysv
2377         IF ( boundary_flags(j,i) == 7 )  THEN
2378         
2379            DO  k = nzb_v_inner(j,i) + 1, nzt
2380               v_comp(k)              = v(k,j,i) + v(k,j-1,i) - gv
2381               swap_flux_y_local_v(k) = v_comp(k) *                           &
2382                                        ( v(k,j,i) + v(k,j-1,i) ) * 0.25
2383               swap_diss_y_local_v(k) = diss_2nd( v(k,j+2,i), v(k,j+1,i),     &
2384                                                  v(k,j,i), v(k,j-1,i),       &
2385                                                  v(k,j-1,i), v_comp(k),      &
2386                                                  0.25, ddy )                               
2387            ENDDO
2388           
2389         ELSE
2390         
2391            DO  k = nzb_v_inner(j,i) + 1, nzt
2392               v_comp(k)              = v(k,j,i) + v(k,j-1,i) - gv
2393               swap_flux_y_local_v(k) = v_comp(k) * (                         &
2394                                      37.  * ( v(k,j,i) + v(k,j-1,i)   )      &
2395                                      - 8. * ( v(k,j+1,i) + v(k,j-2,i) )      &
2396                                      +      ( v(k,j+2,i) + v(k,j-3,i) ) )    &
2397                                      * adv_mom_5
2398               swap_diss_y_local_v(k) = - abs(v_comp(k)) * (                  &
2399                                      10.  * ( v(k,j,i) - v(k,j-1,i)   )      &
2400                                      - 5. * ( v(k,j+1,i) - v(k,j-2,i) )      &
2401                                      +      ( v(k,j+2,i) - v(k,j-3,i) ) )    &
2402                                      * adv_mom_5
2403            ENDDO
2404           
2405         ENDIF
2406         
2407         DO  j = nysv, nyn
2408            IF ( boundary_flags(j,i) /= 0 )  THEN
2409!       
2410!--       Degrade the order for Dirichlet bc. at the outflow boundary
2411               SELECT CASE ( boundary_flags(j,i) )
2412         
2413                  CASE ( 1 )
2414                     DO  k = nzb_v_inner(j,i) + 1, nzt
2415                        u_comp    = u(k,j-1,i+1) + u(k,j,i+1) - gu
2416                        flux_r(k) = u_comp * (                                &
2417                                 7. * (v(k,j,i+1) + v(k,j,i)    )             &
2418                                 -    ( v(k,j,i+2) + v(k,j,i-1) ) ) * adv_mom_3
2419                        diss_r(k) = - abs(u_comp) * (                         &
2420                               3. * ( v(k,j,i+1) - v(k,j,i)   )               & 
2421                               -    ( v(k,j,i+2) - v(k,j,i-1) ) ) * adv_mom_3
2422                     ENDDO
2423               
2424                  CASE ( 2 )
2425                     DO  k = nzb_v_inner(j,i) + 1, nzt
2426                        u_comp    = u(k,j-1,i+1) + u(k,j,i+1) - gu
2427                        flux_r(k) = u_comp * ( v(k,j,i+1) + v(k,j,i) ) * 0.25 
2428                        diss_r(k) = diss_2nd( v(k,j,i+1), v(k,j,i+1),         &
2429                                              v(k,j,i), v(k,j,i-1),           &
2430                                              v(k,j,i-2), u_comp, 0.25, ddx )
2431                     ENDDO
2432               
2433                  CASE ( 3 )
2434                     DO  k = nzb_v_inner(j,i) + 1, nzt
2435                        v_comp(k) = v(k,j+1,i) + v(k,j,i)
2436                        flux_n(k) = ( v_comp(k)- gv ) * (                     &
2437                                 7. * ( v(k,j+1,i) + v(k,j,i)   )             &
2438                                 -    ( v(k,j+2,i) + v(k,j-1,i) ) ) * adv_mom_3
2439                        diss_n(k) = - abs(v_comp(k) - gv) * (                 &
2440                                 3. * ( v(k,j+1,i) - v(k,j,i)   )             & 
2441                                 -    ( v(k,j+2,i) - v(k,j-1,i) ) ) * adv_mom_3
2442                     ENDDO
2443               
2444                  CASE ( 4 )
2445                     DO  k = nzb_v_inner(j,i) + 1, nzt
2446                        v_comp(k) = v(k,j+1,i) + v(k,j,i)
2447                        flux_n(k) = ( v_comp(k) - gv ) *                      &
2448                                   ( v(k,j+1,i) + v(k,j,i) ) * 0.25 
2449                        diss_n(k) = diss_2nd( v(k,j+1,i), v(k,j+1,i),         &
2450                                              v(k,j,i), v(k,j-1,i),           &
2451                                              v(k,j-2,i), v_comp(k), 0.25, ddy)
2452                     ENDDO
2453               
2454                  CASE ( 5 )
2455                     DO  k = nzb_w_inner(j,i) + 1, nzt
2456                        u_comp    = u(k,j-1,i) + u(k,j,i) - gu
2457                        flux_r(k) = u_comp * (                                &
2458                               7. * ( v(k,j,i+1) + v(k,j,i)   )               &
2459                               -    ( v(k,j,i+2) + v(k,j,i-1) ) ) * adv_mom_3
2460                        diss_r(k) = - abs(u_comp) * (                         &
2461                               3. * ( w(k,j,i+1) - w(k,j,i)   )               & 
2462                               -    ( v(k,j,i+2) - v(k,j,i-1) ) ) * adv_mom_3
2463                     ENDDO
2464                 
2465                  CASE ( 6 )
2466                     DO  k = nzb_v_inner(j,i) + 1, nzt
2467                                 
2468                        u_comp = u(k,j-1,i+1) + u(k,j,i+1) - gu
2469                        flux_r(k) = u_comp * (                                &
2470                                 7. * ( v(k,j,i+1) + v(k,j,i)   )             &
2471                                 -    ( v(k,j,i+2) + v(k,j,i-1) ) ) * adv_mom_3
2472                        diss_r(k) = - abs(u_comp) * (                         &
2473                                 3. * ( v(k,j,i+1) - v(k,j,i)   )             & 
2474                                 -    ( v(k,j,i+2) - v(k,j,i-1) ) ) * adv_mom_3
2475                     ENDDO
2476               
2477                  CASE ( 7 )
2478                     DO  k = nzb_v_inner(j,i) + 1, nzt                                 
2479                        v_comp(k) = v(k,j+1,i) + v(k,j,i)
2480                        flux_n(k) = ( v_comp(k) - gv ) * (                    &
2481                                 7. * ( v(k,j+1,i) + v(k,j,i)   )             &
2482                                 -    ( v(k,j+2,i) + v(k,j-1,i) ) ) * adv_mom_3
2483                        diss_n(k) = - abs(v_comp(k) - gv) * (                 &
2484                                 3. * ( v(k,j+1,i) - v(k,j,i)   )             & 
2485                                 -    ( v(k,j+2,i) - v(k,j-1,i) ) ) * adv_mom_3
2486                     ENDDO
2487               
2488                  CASE DEFAULT
2489             
2490               END SELECT
2491!         
2492!--         Compute the crosswise 5th order fluxes at the outflow
2493               IF ( boundary_flags(j,i) == 1 .OR. boundary_flags(j,i) == 2    &
2494            .OR. boundary_flags(j,i) == 5 .OR. boundary_flags(j,i) == 6 )  THEN
2495         
2496                  DO  k = nzb_v_inner(j,i) + 1, nzt
2497                     v_comp(k) = v(k,j+1,i) + v(k,j,i)
2498                     flux_n(k) = ( v_comp(k) - gv ) * (                       &
2499                               37.  * ( v(k,j+1,i) + v(k,j,i)   )             &
2500                               - 8. * ( v(k,j+2,i) + v(k,j-1,i) )             &
2501                               +      ( v(k,j+3,i) + v(k,j-2,i) ) ) * adv_mom_5
2502                     diss_n(k) = - abs(v_comp(k) - gv ) * (                   &
2503                               10.  * ( v(k,j+1,i) - v(k,j,i)   )             &
2504                               - 5. * ( v(k,j+2,i) - v(k,j-1,i) )             &
2505                               +      ( v(k,j+3,i) - v(k,j-2,i) ) ) *adv_mom_5
2506                   ENDDO
2507             
2508               ELSE
2509           
2510                  DO  k = nzb_v_inner(j,i) + 1, nzt
2511                     u_comp    = u(k,j-1,i+1) + u(k,j,i+1) - gu
2512                     flux_r(k) = u_comp * (                                   &
2513                               37.  * ( v(k,j,i+1) + v(k,j,i)   )             &
2514                              - 8. * ( v(k,j,i+2) + v(k,j,i-1) )              &
2515                            +      ( v(k,j,i+3) + v(k,j,i-2) ) ) * adv_mom_5
2516                     diss_r(k) = - abs(u_comp) * (                            &
2517                               10. * ( v(k,j,i+1) - v(k,j,i)   )              &
2518                               -5. * ( v(k,j,i+2) - v(k,j,i-1) )              &
2519                               +     ( v(k,j,i+3) - v(k,j,i-2) ) ) * adv_mom_5
2520                  ENDDO
2521               
2522               ENDIF
2523         
2524         
2525            ELSE
2526         
2527               DO  k = nzb_v_inner(j,i) + 1, nzt
2528                  u_comp    = u(k,j-1,i+1) + u(k,j,i+1) - gu
2529                  flux_r(k) = u_comp * (                                      &
2530                            37.  * ( v(k,j,i+1) + v(k,j,i)   )                &
2531                            - 8. * ( v(k,j,i+2) + v(k,j,i-1) )                &
2532                            +      ( v(k,j,i+3) + v(k,j,i-2) ) ) * adv_mom_5
2533                  diss_r(k) = - abs(u_comp) * (                               &
2534                            10. * ( v(k,j,i+1) - v(k,j,i)   )                 &
2535                            -5. * ( v(k,j,i+2) - v(k,j,i-1) )                 &
2536                            +     ( v(k,j,i+3) - v(k,j,i-2) ) ) * adv_mom_5
2537
2538                  v_comp(k) = v(k,j+1,i) + v(k,j,i)
2539                  flux_n(k) = ( v_comp(k) - gv ) * (                          &
2540                            37.  * ( v(k,j+1,i) + v(k,j,i)   )                &
2541                            - 8. * ( v(k,j+2,i) + v(k,j-1,i) )                &
2542                            +      ( v(k,j+3,i) + v(k,j-2,i) ) ) * adv_mom_5
2543                  diss_n(k) = - abs(v_comp(k) - gv ) * (                      &
2544                            10.  * ( v(k,j+1,i) - v(k,j,i)   )                &
2545                            - 5. * ( v(k,j+2,i) - v(k,j-1,i) )                &
2546                            +      ( v(k,j+3,i) - v(k,j-2,i) ) ) *adv_mom_5
2547
2548               ENDDO
2549            ENDIF
2550           
2551            DO  k = nzb_v_inner(j,i) + 1, nzt
2552               tend(k,j,i) = tend(k,j,i) - (                                  &
2553              ( flux_r(k) + diss_r(k)                                         &
2554            - ( swap_flux_x_local_v(k,j) + swap_diss_x_local_v(k,j) ) ) * ddx &
2555            + ( flux_n(k) + diss_n(k)                                         &
2556            - ( swap_flux_y_local_v(k) + swap_diss_y_local_v(k) )     ) * ddy )
2557                 
2558               swap_flux_x_local_v(k,j) = flux_r(k)   
2559               swap_diss_x_local_v(k,j) = diss_r(k)
2560               swap_flux_y_local_v(k)   = flux_n(k)     
2561               swap_diss_y_local_v(k)   = diss_n(k)   
2562
2563               sums_vs2_ws_l(k,:) = sums_vs2_ws_l(k,:)                        &
2564                  + (flux_n(k) *( v_comp(k) - 2. * hom(k,1,2,:) )             &
2565                  / ( v_comp(k) - gv + 1.0E-20)                               &
2566                  + diss_n(k) * abs(v_comp(k) - 2. * hom(k,1,2,:) )           &
2567                  /(abs(v_comp(k) - gv) + 1.0E-20 ) )                         &
2568                  * weight_substep(intermediate_timestep_count) * rmask(j,i,:)
2569            ENDDO
2570         ENDDO
2571       ENDDO
2572!
2573!--   Vertical advection, degradation of order near surface and top.
2574!--   The fluxes flux_d and diss_d at the surface are 0. Due to reasons of
2575!--   statistical evaluation the top flux at the surface should be 0
2576       DO i = nxl, nxr
2577          DO  j = nysv, nyn
2578!         
2579!--         The fluxes flux_d and diss_d at the surface are 0. Due to static
2580!--         reasons the top flux at the surface should be 0.
2581             flux_t(nzb_v_inner(j,i)) = 0.
2582             diss_t(nzb_v_inner(j,i)) = 0.
2583             k      = nzb_v_inner(j,i) + 1             
2584             flux_d = flux_t(k-1)
2585             diss_d = diss_t(k-1)
2586!
2587!--          2nd order scheme
2588             w_comp    = w(k,j-1,i) + w(k,j,i)
2589             flux_t(k) = w_comp * ( v(k+1,j,i) + v(k,j,i) ) * 0.25
2590             diss_t(k) = diss_2nd( v(k+2,j,i), v(k+1,j,i), v(k,j,i),          &
2591                                   0., 0., w_comp, 0.25, ddzw(k) )
2592
2593             tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)              &
2594                                      - ( flux_d + diss_d )     ) * ddzw(k)
2595!
2596!--         WS3 as an intermediate step
2597             k          = nzb_v_inner(j,i) + 2
2598             flux_d     = flux_t(k-1)
2599             diss_d     = diss_t(k-1)
2600             w_comp     = w(k,j-1,i) + w(k,j,i)
2601             flux_t(k) = w_comp * (                                           &
2602                         7. * ( v(k+1,j,i) + v(k,j,i) )                       &
2603                         -    ( v(k+2,j,i) + v(k-1,j,i) ) ) * adv_mom_3
2604             diss_t(k) = - abs(w_comp) * (                                    &
2605                         3. * ( v(k+1,j,i) - v(k,j,i) )                       & 
2606                         -    ( v(k+2,j,i) - v(k-1,j,i) ) ) * adv_mom_3
2607                         
2608             tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)              &
2609                                      - ( flux_d + diss_d )     ) * ddzw(k)
2610
2611!
2612!-- WS5
2613             DO  k = nzb_v_inner(j,i) + 3, nzt - 3
2614                flux_d = flux_t(k-1)
2615                diss_d = diss_t(k-1)
2616                w_comp = w(k,j-1,i) + w(k,j,i)
2617                flux_t(k) = w_comp * (                                        &
2618                            37.  * ( v(k+1,j,i) + v(k,j,i)   )                &
2619                            - 8. * ( v(k+2,j,i) + v(k-1,j,i) )                &
2620                            +      ( v(k+3,j,i) + v(k-2,j,i) ) ) * adv_mom_5
2621                diss_t(k) = - abs(w_comp) * (                                 &
2622                            10. * ( v(k+1,j,i) - v(k,j,i)   )                 & 
2623                            -5. * ( v(k+2,j,i) - v(k-1,j,i) )                 &
2624                            +     ( v(k+3,j,i) - v(k-2,j,i) ) ) * adv_mom_5
2625                           
2626                tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)              &
2627                                      - ( flux_d + diss_d )     ) * ddzw(k)
2628            ENDDO
2629!
2630!--         WS3 as an intermediate step
2631            DO k = nzt - 2, nzt - 1 
2632               flux_d     = flux_t(k-1)
2633               diss_d     = diss_t(k-1)
2634               w_comp     = w(k,j-1,i) + w(k,j,i)
2635               flux_t(k) = w_comp * (                                         &
2636                           7. * ( v(k+1,j,i) + v(k,j,i) )                     &
2637                           -    ( v(k+2,j,i) + v(k-1,j,i) ) ) * adv_mom_3
2638               diss_t(k) = - abs(w_comp) * (                                  &
2639                           3. * ( v(k+1,j,i) - v(k,j,i) )                     & 
2640                           -    ( v(k+2,j,i) - v(k-1,j,i) ) ) * adv_mom_3
2641                         
2642             tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)              &
2643                                      - ( flux_d + diss_d )     ) * ddzw(k)
2644            ENDDO
2645!
2646!--          2nd order scheme
2647            k         = nzt
2648            flux_d    = flux_t(k-1)
2649            diss_d    = diss_t(k-1)
2650            w_comp    = w(k,j-1,i) + w(k,j,i)
2651            flux_t(k) = w_comp * ( v(k+1,j,i) + v(k,j,i) ) * 0.25
2652            diss_t(k) = diss_2nd( v(nzt+1,j,i), v(nzt+1,j,i), v(k,j,i),       &
2653                                  v(k-1,j,i), v(k-2,j,i), w_comp, 0.25,ddzw(k))
2654
2655            tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)               &
2656                                      - ( flux_d + diss_d )     ) * ddzw(k)
2657!
2658!- at last vertical momentum flux is accumulated
2659            DO  k = nzb_v_inner(j,i), nzt
2660               sums_wsvs_ws_l(k,:) = sums_wsvs_ws_l(k,:)                       &
2661                   + ( flux_t(k) + diss_t(k) )                                 &
2662                   * weight_substep(intermediate_timestep_count)               &
2663                   * rmask(j,i,:)
2664            ENDDO
2665            sums_vs2_ws_l(nzb_v_inner(j,i),:) =                                &
2666                   sums_vs2_ws_l(nzb_v_inner(j,i)+1,:)
2667          ENDDO
2668       ENDDO
2669
2670    END SUBROUTINE advec_v_ws
2671   
2672   
2673!
2674! Advection of w - Call for all grid points
2675!------------------------------------------------------------------------------!
2676    SUBROUTINE advec_w_ws
2677
2678       USE arrays_3d
2679       USE constants
2680       USE control_parameters
2681       USE grid_variables
2682       USE indices
2683       USE statistics
2684
2685       IMPLICIT NONE
2686
2687       INTEGER ::  i, j, k
2688       REAL    ::  gu, gv, flux_d, diss_d, u_comp, v_comp, w_comp
2689       REAL    ::  flux_t(nzb:nzt+1), diss_t(nzb:nzt+1)
2690       REAL, DIMENSION(nzb:nzt+1)  ::  flux_r, diss_r, flux_n, diss_n
2691       REAL, DIMENSION(nzb+1:nzt) :: swap_flux_y_local_w, swap_diss_y_local_w
2692       REAL, DIMENSION(nzb+1:nzt,nys:nyn) :: swap_flux_x_local_w,             &
2693                                             swap_diss_x_local_w
2694 
2695       gu = 2.0 * u_gtrans
2696       gv = 2.0 * v_gtrans
2697!
2698!--   compute the whole left boundary of the processor domain
2699       i = nxl
2700       DO  j = nys, nyn
2701       
2702          IF ( boundary_flags(j,i) == 6 )  THEN
2703         
2704             DO  k = nzb_v_inner(j,i) + 1, nzt
2705                u_comp                   = u(k+1,j,i) + u(k,j,i) - gu
2706                swap_flux_x_local_w(k,j) = u_comp *                           &
2707                                           (w(k,j,i)+w(k,j,i-1)) * 0.25
2708                swap_flux_x_local_w(k,j) = diss_2nd( w(k,j,i+2), w(k,j,i+1),  &
2709                                                     w(k,j,i), w(k,j,i-1),    &
2710                                                     w(k,j,i-1), u_comp,      &
2711                                                     0.25, ddx )
2712             ENDDO
2713             
2714          ELSE
2715         
2716             DO  k = nzb_v_inner(j,i) + 1, nzt
2717                u_comp                   = u(k+1,j,i) + u(k,j,i) - gu
2718                swap_flux_x_local_w(k,j) = u_comp * (                         &
2719                                         37.  * ( w(k,j,i) + w(k,j,i-1)   )   &
2720                                         - 8. * ( w(k,j,i+1) + w(k,j,i-2) )   &
2721                                         +      ( w(k,j,i+2) + w(k,j,i-3) ) ) &
2722                                         * adv_mom_5
2723                swap_diss_x_local_w(k,j) = - abs(u_comp) * (                  &
2724                                         10.  * ( w(k,j,i) - w(k,j,i-1)   )   &
2725                                         - 5. * ( w(k,j,i+1) - w(k,j,i-2) )   &
2726                                         +      ( w(k,j,i+2) - w(k,j,i-3) ) ) &
2727                                         * adv_mom_5
2728             ENDDO
2729             
2730          ENDIF
2731         
2732       ENDDO
2733
2734       DO i = nxl, nxr
2735         j = nys
2736         IF ( boundary_flags(j,i) == 8 )  THEN
2737         
2738            DO  k = nzb_v_inner(j,i) + 1, nzt
2739               v_comp                 = v(k+1,j,i) + v(k,j,i) - gv
2740               swap_flux_y_local_w(k) = v_comp *                              &
2741                                        ( w(k,j,i) + w(k,j-1,i) ) * 0.25
2742               swap_diss_y_local_w(k) = diss_2nd( w(k,j+2,i), w(k,j+1,i),     &
2743                                                  w(k,j,i), w(k,j-1,i),       &
2744                                                  w(k,j-1,i), v_comp, 0.25,ddy)
2745            ENDDO
2746           
2747         ELSE
2748         
2749            DO  k = nzb_v_inner(j,i) + 1, nzt           
2750               v_comp                 = v(k+1,j,i) + v(k,j,i) - gv
2751               swap_flux_y_local_w(k) = v_comp * (                            &
2752                                      37. *  ( w(k,j,i) + w(k,j-1,i)  )       &
2753                                      - 8. * ( w(k,j+1,i) +w(k,j-2,i) )       &
2754                                      +      ( w(k,j+2,i) +w(k,j-3,i) ) )     &
2755                                      * adv_mom_5
2756               swap_diss_y_local_w(k) = - abs(v_comp) * (                     &
2757                                      10.  * ( w(k,j,i) - w(k,j-1,i)   )      &
2758                                      - 5. * ( w(k,j+1,i) - w(k,j-2,i) )      &
2759                                      +      ( w(k,j+2,i) - w(k,j-3,i) ) )    &
2760                                      * adv_mom_5
2761            ENDDO
2762           
2763         ENDIF
2764         
2765         DO  j = nys, nyn
2766
2767            IF ( boundary_flags(j,i) /= 0 )  THEN
2768!       
2769!--            Degrade the order for Dirichlet bc. at the outflow boundary
2770               SELECT CASE ( boundary_flags(j,i) )
2771         
2772                  CASE ( 1 )
2773                     DO  k = nzb_w_inner(j,i) + 1, nzt
2774                        u_comp    = u(k+1,j,i+1) + u(k,j,i+1) - gu
2775                        flux_r(k) = u_comp * (                                &
2776                              7. * ( w(k,j,i+1) + w(k,j,i)   )                &
2777                              -    ( w(k,j,i+2) + w(k,j,i-1) ) ) * adv_mom_3
2778                        diss_r(k) = -abs(u_comp) * (                          &
2779                              3. * ( w(k,j,i+1) - w(k,j,i)   )                & 
2780                              -    ( w(k,j,i+2) - w(k,j,i-1) ) ) * adv_mom_3 
2781                     ENDDO
2782               
2783                  CASE ( 2 )
2784                     DO  k = nzb_w_inner(j,i) + 1, nzt
2785                        u_comp    = u(k+1,j,i+1) + u(k,j,i+1) - gu
2786                        flux_r(k) = u_comp * ( w(k,j,i+1) + w(k,j,i) ) * 0.25 
2787                        diss_r(k) = diss_2nd( w(k,j,i+1), w(k,j,i+1),         &
2788                                              w(k,j,i), w(k,j,i-1),           &
2789                                              w(k,j,i-2), u_comp, 0.25, ddx )
2790                     ENDDO
2791               
2792                  CASE ( 3 )
2793                     DO  k = nzb_w_inner(j,i) + 1, nzt
2794                        v_comp    = v(k+1,j+1,i) + v(k,j+1,i) - gv
2795                        flux_n(k) = v_comp * (                                &
2796                              7. * ( w(k,j+1,i) + w(k,j,i)   )                &
2797                              -    ( w(k,j+2,i) + w(k,j-1,i) ) ) * adv_mom_3
2798                        diss_n(k) = -abs(v_comp) * (                          &
2799                              3. * ( w(k,j+1,i) - w(k,j,i)   )                & 
2800                              -    ( w(k,j+2,i) - w(k,j-1,i) ) ) * adv_mom_3
2801                     ENDDO
2802               
2803                  CASE ( 4 )
2804                     DO  k = nzb_w_inner(j,i) + 1, nzt
2805                        v_comp    = v(k+1,j+1,i) + v(k,j+1,i) - gv
2806                        flux_n(k) = v_comp * ( w(k,j+1,i) + w(k,j,i) ) * 0.25 
2807                        diss_n(k) = diss_2nd( w(k,j+1,i), w(k,j+1,i),         &
2808                                              w(k,j,i), w(k,j-1,i),           &
2809                                              w(k,j-2,i), v_comp, 0.25, ddy )
2810                     ENDDO
2811               
2812                  CASE ( 5 )
2813                     DO  k = nzb_w_inner(j,i) + 1, nzt
2814                        u_comp    = u(k+1,j,i+1) + u(k,j,i+1) - gu
2815                        flux_r(k) = u_comp * (                                &
2816                              7. * ( w(k,j,i+1) + w(k,j,i)   )                &
2817                              -    ( w(k,j,i+2) + w(k,j,i-1) ) ) * adv_mom_3
2818                        diss_r(k) = - abs(u_comp) * (                         &
2819                              3. * ( w(k,j,i+1) - w(k,j,i) )                  & 
2820                              -    ( w(k,j,i+2) - w(k,j,i-1) ) ) * adv_mom_3
2821                     ENDDO
2822               
2823                  CASE ( 6 )
2824                     DO  k = nzb_w_inner(j,i) + 1, nzt
2825                        u_comp    = u(k+1,j,i+1) + u(k,j,i+1) - gu
2826                        flux_r(k) = u_comp *(                                 &
2827                              7. * ( w(k,j,i+1) + w(k,j,i)   )                &
2828                              -    ( w(k,j,i+2) + w(k,j,i-1) ) )  * adv_mom_3
2829                        diss_r(k) = - abs(u_comp) * (                         &
2830                              3. * ( w(k,j,i+1) - w(k,j,i) )                  & 
2831                              -    ( w(k,j,i+2) - w(k,j,i-1) ) ) * adv_mom_3
2832                     ENDDO
2833               
2834                  CASE ( 7 )
2835                     DO  k = nzb_w_inner(j,i) + 1, nzt
2836                        v_comp    = v(k+1,j+1,i) + v(k,j+1,i) - gv
2837                        flux_n(k) = v_comp *(                                 &
2838                              7. * ( w(k,j+1,i) + w(k,j,i)   )                &
2839                              -    ( w(k,j+2,i) + w(k,j-1,i) ) ) * adv_mom_3
2840                        diss_n(k) = - abs(v_comp) * (                         &
2841                              3. * ( w(k,j+1,i) - w(k,j,i)   )                & 
2842                              -    ( w(k,j+2,i) - w(k,j-1,i) ) ) * adv_mom_3
2843                      ENDDO
2844               
2845                  CASE ( 8 )
2846                     DO  k = nzb_w_inner(j,i) + 1, nzt
2847                        v_comp    = v(k+1,j+1,i) + v(k,j+1,i) - gv
2848                        flux_n(k) = v_comp * (                                &
2849                              7. * ( w(k,j+1,i) + w(k,j,i)   )                &
2850                             -     ( w(k,j+2,i) + w(k,j-1,i) ) ) * adv_mom_3
2851                        diss_n(k) = - abs(v_comp) * (                         &
2852                              3. * ( w(k,j+1,i) - w(k,j,i) )                  & 
2853                              -    ( w(k,j+2,i) - w(k,j-1,i) ) ) * adv_mom_3
2854               
2855                     ENDDO 
2856               
2857                  CASE DEFAULT
2858           
2859               END SELECT
2860!         
2861!--      Compute the crosswise 5th order fluxes at the outflow
2862               IF ( boundary_flags(j,i) == 1 .OR. boundary_flags(j,i) == 2    &
2863            .OR. boundary_flags(j,i) == 5 .OR. boundary_flags(j,i) == 6 )  THEN
2864         
2865                  DO  k = nzb_w_inner(j,i) + 1, nzt
2866                     v_comp    = v(k+1,j+1,i) + v(k,j+1,i) - gv
2867                     flux_n(k) = v_comp * (                                   &
2868                           37.  * ( w(k,j+1,i) + w(k,j,i)   )                 &
2869                           - 8. * ( w(k,j+2,i) + w(k,j-1,i) )                 &
2870                           +      ( w(k,j+3,i) + w(k,j-2,i) ) ) * adv_mom_5
2871                     diss_n(k) = - abs(v_comp) * (                            &
2872                           10.  * ( w(k,j+1,i) - w(k,j,i)   )                 &
2873                           - 5. * ( w(k,j+2,i) - w(k,j-1,i) )                 &
2874                           +      ( w(k,j+3,i) - w(k,j-2,i) ) ) * adv_mom_5
2875                  ENDDO
2876           
2877               ELSE
2878         
2879                  DO  k = nzb_w_inner(j,i) + 1, nzt         
2880                     u_comp    = u(k+1,j,i+1) + u(k,j,i+1) - gu
2881                     flux_r(k) = u_comp * (                                   &
2882                           37.  * ( w(k,j,i+1) + w(k,j,i) )                   &
2883                           - 8. * ( w(k,j,i+2) + w(k,j,i-1) )                 &
2884                           +      ( w(k,j,i+3) + w(k,j,i-2) ) ) * adv_mom_5
2885                     diss_r(k) = - abs(u_comp) * (                            &
2886                           10.  * ( w(k,j,i+1) - w(k,j,i)   )                 &
2887                           - 5. * ( w(k,j,i+2) - w(k,j,i-1) )                 &
2888                            +      ( w(k,j,i+3) - w(k,j,i-2) ) ) * adv_mom_5
2889                  ENDDO
2890           
2891               ENDIF
2892         
2893             ELSE
2894!       
2895!--            Compute the fifth order fluxes for the interior of PE domain.               
2896               DO  k = nzb_w_inner(j,i) + 1, nzt
2897                  u_comp    = u(k+1,j,i+1) + u(k,j,i+1) - gu
2898                  flux_r(k) = u_comp * (                                      &
2899                        37.  * ( w(k,j,i+1) + w(k,j,i)   )                    &
2900                        - 8. * ( w(k,j,i+2) + w(k,j,i-1) )                    &
2901                        +      ( w(k,j,i+3) + w(k,j,i-2) ) ) * adv_mom_5
2902                  diss_r(k) = - abs(u_comp) * (                               &
2903                        10.  * ( w(k,j,i+1) - w(k,j,i)   )                    &
2904                        - 5. * ( w(k,j,i+2) - w(k,j,i-1) )                    &
2905                        +      ( w(k,j,i+3) - w(k,j,i-2) ) ) * adv_mom_5
2906                 
2907                  v_comp    = v(k+1,j+1,i) + v(k,j+1,i) - gv
2908                  flux_n(k) = v_comp * (                                      &
2909                        37.  * ( w(k,j+1,i) + w(k,j,i)   )                    &
2910                        - 8. * ( w(k,j+2,i) + w(k,j-1,i) )                    &
2911                        +      ( w(k,j+3,i) + w(k,j-2,i) ) ) * adv_mom_5
2912                  diss_n(k) = - abs(v_comp) * (                               &
2913                        10.  * ( w(k,j+1,i) - w(k,j,i)   )                    &
2914                        - 5. * ( w(k,j+2,i) - w(k,j-1,i) )                    &
2915                        +      ( w(k,j+3,i) - w(k,j-2,i) ) ) * adv_mom_5
2916               ENDDO
2917         
2918            ENDIF
2919           
2920            DO  k = nzb_v_inner(j,i)+1, nzt
2921               tend(k,j,i) = tend(k,j,i) - (                                  &
2922              ( flux_r(k) + diss_r(k)                                         &
2923            - ( swap_flux_x_local_w(k,j) + swap_diss_x_local_w(k,j) ) ) * ddx &
2924            + ( flux_n(k) + diss_n(k)                                         &
2925            - ( swap_flux_y_local_w(k) + swap_diss_y_local_w(k)     ) ) * ddy )
2926               
2927               swap_flux_x_local_w(k,j) = flux_r(k)
2928               swap_diss_x_local_w(k,j) = diss_r(k)
2929               swap_flux_y_local_w(k)   = flux_n(k) 
2930               swap_diss_y_local_w(k)   = diss_n(k)
2931            ENDDO
2932         ENDDO
2933       ENDDO
2934
2935!
2936!--    Vertical advection, degradation of order near surface and top.
2937!--    The fluxes flux_d and diss_d at the surface are 0. Due to reasons of
2938!--    statistical evaluation the top flux at the surface should be 0
2939       DO i = nxl, nxr
2940          DO  j = nys, nyn
2941            k      = nzb_w_inner(j,i) + 1
2942            flux_d = w(k-1,j,i) * ( w(k,j,i) + w(k-1,j,i))
2943            flux_t(k-1) = flux_d
2944            diss_d = 0.0
2945            diss_t(k-1) = diss_d
2946!
2947!--         2nd order scheme           
2948            w_comp = w(k+1,j,i) + w(k,j,i)
2949            flux_t(k) = w_comp * ( w(k+1,j,i) + w(k,j,i) ) * 0.25
2950            diss_t(k) = diss_2nd( w(k+2,j,i), w(k+1,j,i), w(k,j,i), 0., 0.,   &
2951                                  w_comp, 0.25, ddzu(k+1) )
2952
2953            tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)               &
2954                                 - ( flux_d + diss_d )     ) * ddzu(k+1)
2955!                                 
2956!--        WS3 as an intermediate step
2957            k         = nzb_w_inner(j,i) + 2
2958            flux_d    = flux_t(k-1)
2959            diss_d    = diss_t(k-1)
2960            w_comp    = w(k+1,j,i)+w(k,j,i)
2961            flux_t(k) = w_comp * (                                            &
2962                        7. * ( w(k+1,j,i) + w(k,j,i) )                        &
2963                        - ( w(k+2,j,i) + w(k-1,j,i) ) ) * adv_mom_3
2964            diss_t(k) = - abs(w_comp) * (                                     &
2965                        3. * ( w(k+1,j,i) - w(k,j,i) )                        & 
2966                        - ( w(k+2,j,i) - w(k-1,j,i) ) ) * adv_mom_3
2967                       
2968            tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)               &
2969                                 - ( flux_d + diss_d )     ) * ddzu(k+1)
2970!
2971!--         WS5
2972            DO  k = nzb_v_inner(j,i) + 3, nzt
2973               flux_d = flux_t(k-1)
2974               diss_d = diss_t(k-1)
2975
2976               w_comp = w(k+1,j,i) + w(k,j,i)
2977               flux_t(k) = w_comp * (                                         &
2978                           37.  * ( w(k+1,j,i) + w(k,j,i)   )                 &
2979                           - 8. * ( w(k+2,j,i) + w(k-1,j,i) )                 &
2980                           +      ( w(k+3,j,i) + w(k-2,j,i) ) ) * adv_mom_5
2981               diss_t(k) = - abs(w_comp) * (                                  &
2982                           10.  * ( w(k+1,j,i) - w(k,j,i)   )                 & 
2983                           - 5. * ( w(k+2,j,i) - w(k-1,j,i) )                 &
2984                           +      ( w(k+3,j,i) - w(k-2,j,i) ) ) * adv_mom_5
2985
2986               tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)            &
2987                                 - ( flux_d + diss_d )     ) * ddzu(k+1)
2988            ENDDO
2989!                                 
2990!--        WS3 as an intermediate step
2991            DO k = nzt - 2, nzt - 1
2992               flux_d    = flux_t(k-1)
2993               diss_d    = diss_t(k-1)
2994               w_comp    = w(k+1,j,i)+w(k,j,i)
2995               flux_t(k) = w_comp * (                                         &
2996                           7. * ( w(k+1,j,i) + w(k,j,i) )                     &
2997                           - ( w(k+2,j,i) + w(k-1,j,i) ) ) * adv_mom_3
2998               diss_t(k) = - abs(w_comp) * (                                  &
2999                           3. * ( w(k+1,j,i) - w(k,j,i) )                     & 
3000                           - ( w(k+2,j,i) - w(k-1,j,i) ) ) * adv_mom_3
3001                       
3002               tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)            &
3003                                 - ( flux_d + diss_d )     ) * ddzu(k+1)
3004            ENDDO
3005!
3006!--         2nd order scheme 
3007            k = nzt
3008            flux_d = flux_t(k-1)
3009            diss_d = diss_t(k-1)
3010            w_comp = w(k+1,j,i) + w(k,j,i)
3011            flux_t(k) = w_comp * ( w(k+1,j,i) + w(k,j,i) ) * 0.25
3012            diss_t(k) = diss_2nd( w(nzt+1,j,i), w(nzt+1,j,i), w(k,j,i),       &
3013                                  w(k-1,j,i), w(k-2,j,i), w_comp,             &
3014                                  0.25, ddzu(k+1) ) 
3015
3016            tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k)               &
3017                                 - ( flux_d + diss_d )     ) * ddzu(k+1)
3018!
3019!--         at last vertical momentum flux is accumulated
3020            DO  k = nzb_w_inner(j,i), nzt
3021               sums_ws2_ws_l(k,:)  = sums_ws2_ws_l(k,:)                        &
3022                   + ( flux_t(k) + diss_t(k) )                                 &
3023                   * weight_substep(intermediate_timestep_count)               &
3024                   * rmask(j,i,:)
3025            ENDDO
3026
3027          ENDDO
3028       ENDDO
3029
3030    END SUBROUTINE advec_w_ws
3031   
3032 
3033 
3034    SUBROUTINE local_diss_ij( i, j, ar, var_char )
3035   
3036    END SUBROUTINE local_diss_ij
3037
3038    SUBROUTINE local_diss( ar )
3039
3040
3041    END SUBROUTINE local_diss
3042   
3043!   
3044!-- Computation of 2nd order dissipation. This numerical dissipation is
3045!-- necessary to keep a stable numerical solution in region where the
3046!-- order of the schemes degraded.
3047 
3048     REAL FUNCTION diss_2nd( v2, v1, v0, vm1, vm2, vel_comp, factor, grid )   &
3049                            RESULT(dissip)
3050       
3051       IMPLICIT NONE
3052       
3053       REAL, INTENT(IN)  :: v2, v1, v0, vm1, vm2, vel_comp, factor,            &
3054                            grid
3055       REAL :: value_min_m, value_max_m, value_min, value_max,                 &
3056               value_min_p, value_max_p, diss_m, diss_0, diss_p
3057             
3058       value_min_m = MIN(v0, vm1, vm2)
3059       value_max_m = MAX(v0, vm1, vm2)
3060       value_min = MIN(v1, v0, vm2)
3061       value_max = MAX(v1, v0, vm2)
3062       value_min_p = MIN(v2, v1, v0)
3063       value_max_p = MAX(v2, v1, v0)
3064       
3065       diss_m = MAX(0.,vm1 - value_max_m) + MIN(0.,vm1 - value_min_m)
3066       diss_0 = MAX(0.,v0 - value_max) + MIN(0.,v0 - value_min)
3067       diss_p = MAX(0.,v1 - value_max_p) + MIN(0.,v1 - value_min_p)
3068       
3069       dissip = abs(vel_comp) * (diss_p - 2.*diss_0 + diss_m)         &
3070                         * grid**2 * factor
3071             
3072    END FUNCTION diss_2nd
3073
3074 END MODULE advec_ws
Note: See TracBrowser for help on using the repository browser.