source: palm/trunk/SOURCE/boundary_conds.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:keywords set to Id
File size: 19.4 KB
Line 
1 SUBROUTINE boundary_conds( range )
2
3!------------------------------------------------------------------------------!
4! Current revisions:
5! -----------------
6!
7! nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng
8!
9!
10! Removed mirror boundary conditions for u and v at the bottom in case of
11! ibc_uv_b == 0. Instead, dirichelt boundary conditions (u=v=0) are set
12! in init_3d_model
13
14! Former revisions:
15! -----------------
16! $Id: boundary_conds.f90 667 2010-12-23 12:06:00Z suehring $
17!
18! 107 2007-08-17 13:54:45Z raasch
19! Boundary conditions for temperature adjusted for coupled runs,
20! bugfixes for the radiation boundary conditions at the outflow: radiation
21! conditions are used for every substep, phase speeds are calculated for the
22! first Runge-Kutta substep only and then reused, several index values changed
23!
24! 95 2007-06-02 16:48:38Z raasch
25! Boundary conditions for salinity added
26!
27! 75 2007-03-22 09:54:05Z raasch
28! The "main" part sets conditions for time level t+dt instead of level t,
29! outflow boundary conditions changed from Neumann to radiation condition,
30! uxrp, vynp eliminated, moisture renamed humidity
31!
32! 19 2007-02-23 04:53:48Z raasch
33! Boundary conditions for e(nzt), pt(nzt), and q(nzt) removed because these
34! gridpoints are now calculated by the prognostic equation,
35! Dirichlet and zero gradient condition for pt established at top boundary
36!
37! RCS Log replace by Id keyword, revision history cleaned up
38!
39! Revision 1.15  2006/02/23 09:54:55  raasch
40! Surface boundary conditions in case of topography: nzb replaced by
41! 2d-k-index-arrays (nzb_w_inner, etc.). Conditions for u and v remain
42! unchanged (still using nzb) because a non-flat topography must use a
43! Prandtl-layer, which don't requires explicit setting of the surface values.
44!
45! Revision 1.1  1997/09/12 06:21:34  raasch
46! Initial revision
47!
48!
49! Description:
50! ------------
51! Boundary conditions for the prognostic quantities (range='main').
52! In case of non-cyclic lateral boundaries the conditions for velocities at
53! the outflow are set after the pressure solver has been called (range=
54! 'outflow_uvw').
55! One additional bottom boundary condition is applied for the TKE (=(u*)**2)
56! in prandtl_fluxes. The cyclic lateral boundary conditions are implicitly
57! handled in routine exchange_horiz. Pressure boundary conditions are
58! explicitly set in routines pres, poisfft, poismg and sor.
59!------------------------------------------------------------------------------!
60
61    USE arrays_3d
62    USE control_parameters
63    USE grid_variables
64    USE indices
65    USE pegrid
66
67    IMPLICIT NONE
68
69    CHARACTER (LEN=*) ::  range
70
71    INTEGER ::  i, j, k
72
73    REAL    ::  c_max, denom
74
75
76    IF ( range == 'main')  THEN
77!
78!--    Bottom boundary
79       IF ( ibc_uv_b == 1 )  THEN
80          u_p(nzb,:,:) = u_p(nzb+1,:,:)
81          v_p(nzb,:,:) = v_p(nzb+1,:,:)
82       ENDIF
83       DO  i = nxlg, nxrg
84          DO  j = nysg, nyng
85             w_p(nzb_w_inner(j,i),j,i) = 0.0
86          ENDDO
87       ENDDO
88
89!
90!--    Top boundary
91       IF ( ibc_uv_t == 0 )  THEN
92           u_p(nzt+1,:,:) = ug(nzt+1)
93           v_p(nzt+1,:,:) = vg(nzt+1)
94       ELSE
95           u_p(nzt+1,:,:) = u_p(nzt,:,:)
96           v_p(nzt+1,:,:) = v_p(nzt,:,:)
97       ENDIF
98       w_p(nzt:nzt+1,:,:) = 0.0  ! nzt is not a prognostic level (but cf. pres)
99
100!
101!--    Temperature at bottom boundary.
102!--    In case of coupled runs (ibc_pt_b = 2) the temperature is given by
103!--    the sea surface temperature of the coupled ocean model.
104       IF ( ibc_pt_b == 0 )  THEN
105          DO  i = nxlg, nxrg
106             DO  j = nysg, nyng
107                pt_p(nzb_s_inner(j,i),j,i) = pt(nzb_s_inner(j,i),j,i)
108             ENDDO
109          ENDDO
110       ELSEIF ( ibc_pt_b == 1 )  THEN
111          DO  i = nxlg, nxrg
112             DO  j = nysg, nyng
113                pt_p(nzb_s_inner(j,i),j,i) = pt_p(nzb_s_inner(j,i)+1,j,i)
114             ENDDO
115          ENDDO
116       ENDIF
117
118!
119!--    Temperature at top boundary
120       IF ( ibc_pt_t == 0 )  THEN
121           pt_p(nzt+1,:,:) = pt(nzt+1,:,:)
122       ELSEIF ( ibc_pt_t == 1 )  THEN
123           pt_p(nzt+1,:,:) = pt_p(nzt,:,:)
124       ELSEIF ( ibc_pt_t == 2 )  THEN
125           pt_p(nzt+1,:,:) = pt_p(nzt,:,:)   + bc_pt_t_val * dzu(nzt+1)
126       ENDIF
127
128!
129!--    Boundary conditions for TKE
130!--    Generally Neumann conditions with de/dz=0 are assumed
131       IF ( .NOT. constant_diffusion )  THEN
132          DO  i = nxlg, nxrg
133             DO  j = nysg, nyng
134                e_p(nzb_s_inner(j,i),j,i) = e_p(nzb_s_inner(j,i)+1,j,i)
135             ENDDO
136          ENDDO
137          e_p(nzt+1,:,:) = e_p(nzt,:,:)
138       ENDIF
139
140!
141!--    Boundary conditions for salinity
142       IF ( ocean )  THEN
143!
144!--       Bottom boundary: Neumann condition because salinity flux is always
145!--       given
146          DO  i = nxlg, nxrg
147             DO  j = nysg, nyng
148                sa_p(nzb_s_inner(j,i),j,i) = sa_p(nzb_s_inner(j,i)+1,j,i)
149             ENDDO
150          ENDDO
151
152!
153!--       Top boundary: Dirichlet or Neumann
154          IF ( ibc_sa_t == 0 )  THEN
155              sa_p(nzt+1,:,:) = sa(nzt+1,:,:)
156          ELSEIF ( ibc_sa_t == 1 )  THEN
157              sa_p(nzt+1,:,:) = sa_p(nzt,:,:)
158          ENDIF
159
160       ENDIF
161
162!
163!--    Boundary conditions for total water content or scalar,
164!--    bottom and top boundary (see also temperature)
165       IF ( humidity  .OR.  passive_scalar )  THEN
166!
167!--       Surface conditions for constant_humidity_flux
168          IF ( ibc_q_b == 0 ) THEN
169             DO  i = nxlg, nxrg
170                DO  j = nysg, nyng
171                   q_p(nzb_s_inner(j,i),j,i) = q(nzb_s_inner(j,i),j,i)
172                ENDDO
173             ENDDO
174          ELSE
175             DO  i = nxlg, nxrg
176                DO  j = nysg, nyng
177                   q_p(nzb_s_inner(j,i),j,i) = q_p(nzb_s_inner(j,i)+1,j,i)
178                ENDDO
179             ENDDO
180          ENDIF
181!
182!--       Top boundary
183          q_p(nzt+1,:,:) = q_p(nzt,:,:)   + bc_q_t_val * dzu(nzt+1)
184
185
186       ENDIF
187
188!
189!--    Lateral boundary conditions at the inflow. Quasi Neumann conditions
190!--    are needed for the wall normal velocity in order to ensure zero
191!--    divergence. Dirichlet conditions are used for all other quantities.
192       IF ( inflow_s )  THEN
193          v_p(:,nys,:) = v_p(:,nys-1,:)
194       ELSEIF ( inflow_n ) THEN
195          v_p(:,nyn,:) = v_p(:,nyn+1,:)
196       ELSEIF ( inflow_l ) THEN
197          u_p(:,:,nxl) = u_p(:,:,nxl-1)
198       ELSEIF ( inflow_r ) THEN
199          u_p(:,:,nxr) = u_p(:,:,nxr+1)
200       ENDIF
201
202!
203!--    Lateral boundary conditions for scalar quantities at the outflow
204       IF ( outflow_s )  THEN
205          pt_p(:,nys-1,:)     = pt_p(:,nys,:)
206          IF ( .NOT. constant_diffusion     )  e_p(:,nys-1,:) = e_p(:,nys,:)
207          IF ( humidity .OR. passive_scalar )  q_p(:,nys-1,:) = q_p(:,nys,:)
208       ELSEIF ( outflow_n )  THEN
209          pt_p(:,nyn+1,:)     = pt_p(:,nyn,:)
210          IF ( .NOT. constant_diffusion     )  e_p(:,nyn+1,:) = e_p(:,nyn,:)
211          IF ( humidity .OR. passive_scalar )  q_p(:,nyn+1,:) = q_p(:,nyn,:)
212       ELSEIF ( outflow_l )  THEN
213          pt_p(:,:,nxl-1)     = pt_p(:,:,nxl)
214          IF ( .NOT. constant_diffusion     )  e_p(:,:,nxl-1) = e_p(:,:,nxl)
215          IF ( humidity .OR. passive_scalar )  q_p(:,:,nxl-1) = q_p(:,:,nxl)
216       ELSEIF ( outflow_r )  THEN
217          pt_p(:,:,nxr+1)     = pt_p(:,:,nxr)
218          IF ( .NOT. constant_diffusion     )  e_p(:,:,nxr+1) = e_p(:,:,nxr)
219          IF ( humidity .OR. passive_scalar )  q_p(:,:,nxr+1) = q_p(:,:,nxr)
220       ENDIF
221
222    ENDIF
223
224!
225!-- Radiation boundary condition for the velocities at the respective outflow
226    IF ( outflow_s )  THEN
227
228       c_max = dy / dt_3d
229
230       DO i = nxlg, nxrg
231          DO k = nzb+1, nzt+1
232
233!
234!--          Calculate the phase speeds for u,v, and w. In case of using a
235!--          Runge-Kutta scheme, do this for the first substep only and then
236!--          reuse this values for the further substeps.
237             IF ( intermediate_timestep_count == 1 )  THEN
238
239                denom = u_m_s(k,0,i) - u_m_s(k,1,i)
240
241                IF ( denom /= 0.0 )  THEN
242                   c_u(k,i) = -c_max * ( u(k,0,i) - u_m_s(k,0,i) ) / denom
243                   IF ( c_u(k,i) < 0.0 )  THEN
244                      c_u(k,i) = 0.0
245                   ELSEIF ( c_u(k,i) > c_max )  THEN
246                      c_u(k,i) = c_max
247                   ENDIF
248                ELSE
249                   c_u(k,i) = c_max
250                ENDIF
251
252                denom = v_m_s(k,1,i) - v_m_s(k,2,i)
253
254                IF ( denom /= 0.0 )  THEN
255                   c_v(k,i) = -c_max * ( v(k,1,i) - v_m_s(k,1,i) ) / denom
256                   IF ( c_v(k,i) < 0.0 )  THEN
257                      c_v(k,i) = 0.0
258                   ELSEIF ( c_v(k,i) > c_max )  THEN
259                      c_v(k,i) = c_max
260                   ENDIF
261                ELSE
262                   c_v(k,i) = c_max
263                ENDIF
264
265                denom = w_m_s(k,0,i) - w_m_s(k,1,i)
266
267                IF ( denom /= 0.0 )  THEN
268                   c_w(k,i) = -c_max * ( w(k,0,i) - w_m_s(k,0,i) ) / denom
269                   IF ( c_w(k,i) < 0.0 )  THEN
270                      c_w(k,i) = 0.0
271                   ELSEIF ( c_w(k,i) > c_max )  THEN
272                      c_w(k,i) = c_max
273                   ENDIF
274                ELSE
275                   c_w(k,i) = c_max
276                ENDIF
277
278!
279!--             Save old timelevels for the next timestep
280                u_m_s(k,:,i) = u(k,0:1,i)
281                v_m_s(k,:,i) = v(k,1:2,i)
282                w_m_s(k,:,i) = w(k,0:1,i)
283
284             ENDIF
285
286!
287!--          Calculate the new velocities
288             u_p(k,-1,i) = u(k,-1,i) - dt_3d * tsc(2) * c_u(k,i) * &
289                                       ( u(k,-1,i) - u(k,0,i) ) * ddy
290
291             v_p(k,0,i)  = v(k,0,i)  - dt_3d * tsc(2) * c_v(k,i) * &
292                                       ( v(k,0,i) - v(k,1,i) ) * ddy
293
294             w_p(k,-1,i) = w(k,-1,i) - dt_3d * tsc(2) * c_w(k,i) * &
295                                       ( w(k,-1,i) - w(k,0,i) ) * ddy
296
297          ENDDO
298       ENDDO
299
300!
301!--    Bottom boundary at the outflow
302       IF ( ibc_uv_b == 0 )  THEN
303          u_p(nzb,-1,:) = 0.0 
304          v_p(nzb,0,:)  = 0.0 
305       ELSE                   
306          u_p(nzb,-1,:) =  u_p(nzb+1,-1,:)
307          v_p(nzb,0,:)  =  v_p(nzb+1,0,:)
308       ENDIF
309       w_p(nzb,-1,:) = 0.0
310
311!
312!--    Top boundary at the outflow
313       IF ( ibc_uv_t == 0 )  THEN
314          u_p(nzt+1,-1,:) = ug(nzt+1)
315          v_p(nzt+1,0,:)  = vg(nzt+1)
316       ELSE
317          u_p(nzt+1,-1,:) = u(nzt,-1,:)
318          v_p(nzt+1,0,:)  = v(nzt,0,:)
319       ENDIF
320       w_p(nzt:nzt+1,-1,:) = 0.0
321
322    ENDIF
323
324    IF ( outflow_n )  THEN
325
326       c_max = dy / dt_3d
327
328       DO i = nxlg, nxrg
329          DO k = nzb+1, nzt+1
330
331!
332!--          Calculate the phase speeds for u,v, and w. In case of using a
333!--          Runge-Kutta scheme, do this for the first substep only and then
334!--          reuse this values for the further substeps.
335             IF ( intermediate_timestep_count == 1 )  THEN
336
337                denom = u_m_n(k,ny,i) - u_m_n(k,ny-1,i)
338
339                IF ( denom /= 0.0 )  THEN
340                   c_u(k,i) = -c_max * ( u(k,ny,i) - u_m_n(k,ny,i) ) / denom
341                   IF ( c_u(k,i) < 0.0 )  THEN
342                      c_u(k,i) = 0.0
343                   ELSEIF ( c_u(k,i) > c_max )  THEN
344                      c_u(k,i) = c_max
345                   ENDIF
346                ELSE
347                   c_u(k,i) = c_max
348                ENDIF
349
350                denom = v_m_n(k,ny,i) - v_m_n(k,ny-1,i)
351
352                IF ( denom /= 0.0 )  THEN
353                   c_v(k,i) = -c_max * ( v(k,ny,i) - v_m_n(k,ny,i) ) / denom
354                   IF ( c_v(k,i) < 0.0 )  THEN
355                      c_v(k,i) = 0.0
356                   ELSEIF ( c_v(k,i) > c_max )  THEN
357                      c_v(k,i) = c_max
358                   ENDIF
359                ELSE
360                   c_v(k,i) = c_max
361                ENDIF
362
363                denom = w_m_n(k,ny,i) - w_m_n(k,ny-1,i)
364
365                IF ( denom /= 0.0 )  THEN
366                   c_w(k,i) = -c_max * ( w(k,ny,i) - w_m_n(k,ny,i) ) / denom
367                   IF ( c_w(k,i) < 0.0 )  THEN
368                      c_w(k,i) = 0.0
369                   ELSEIF ( c_w(k,i) > c_max )  THEN
370                      c_w(k,i) = c_max
371                   ENDIF
372                ELSE
373                   c_w(k,i) = c_max
374                ENDIF
375
376!
377!--             Swap timelevels for the next timestep
378                u_m_n(k,:,i) = u(k,ny-1:ny,i)
379                v_m_n(k,:,i) = v(k,ny-1:ny,i)
380                w_m_n(k,:,i) = w(k,ny-1:ny,i)
381
382             ENDIF
383
384!
385!--          Calculate the new velocities
386             u_p(k,ny+1,i) = u(k,ny+1,i) - dt_3d * tsc(2) * c_u(k,i) * &
387                                           ( u(k,ny+1,i) - u(k,ny,i) ) * ddy
388
389             v_p(k,ny+1,i) = v(k,ny+1,i) - dt_3d * tsc(2) * c_v(k,i) * &
390                                           ( v(k,ny+1,i) - v(k,ny,i) ) * ddy
391
392             w_p(k,ny+1,i) = w(k,ny+1,i) - dt_3d * tsc(2) * c_w(k,i) * &
393                                           ( w(k,ny+1,i) - w(k,ny,i) ) * ddy
394
395          ENDDO
396       ENDDO
397
398!
399!--    Bottom boundary at the outflow
400       IF ( ibc_uv_b == 0 )  THEN
401          u_p(nzb,ny+1,:) = 0.0 
402          v_p(nzb,ny+1,:) = 0.0   
403       ELSE                   
404          u_p(nzb,ny+1,:) =  u_p(nzb+1,ny+1,:)
405          v_p(nzb,ny+1,:) =  v_p(nzb+1,ny+1,:)
406       ENDIF
407       w_p(nzb,ny+1,:) = 0.0
408
409!
410!--    Top boundary at the outflow
411       IF ( ibc_uv_t == 0 )  THEN
412          u_p(nzt+1,ny+1,:) = ug(nzt+1)
413          v_p(nzt+1,ny+1,:) = vg(nzt+1)
414       ELSE
415          u_p(nzt+1,ny+1,:) = u_p(nzt,nyn+1,:)
416          v_p(nzt+1,ny+1,:) = v_p(nzt,nyn+1,:)
417       ENDIF
418       w_p(nzt:nzt+1,ny+1,:) = 0.0
419
420    ENDIF
421
422    IF ( outflow_l )  THEN
423
424       c_max = dx / dt_3d
425
426       DO j = nysg, nyng
427          DO k = nzb+1, nzt+1
428
429!
430!--          Calculate the phase speeds for u,v, and w. In case of using a
431!--          Runge-Kutta scheme, do this for the first substep only and then
432!--          reuse this values for the further substeps.
433             IF ( intermediate_timestep_count == 1 )  THEN
434
435                denom = u_m_l(k,j,1) - u_m_l(k,j,2)
436
437                IF ( denom /= 0.0 )  THEN
438                   c_u(k,j) = -c_max * ( u(k,j,1) - u_m_l(k,j,1) ) / denom
439                   IF ( c_u(k,j) < 0.0 )  THEN
440                      c_u(k,j) = 0.0
441                   ELSEIF ( c_u(k,j) > c_max )  THEN
442                      c_u(k,j) = c_max
443                   ENDIF
444                ELSE
445                   c_u(k,j) = c_max
446                ENDIF
447
448                denom = v_m_l(k,j,0) - v_m_l(k,j,1)
449
450                IF ( denom /= 0.0 )  THEN
451                   c_v(k,j) = -c_max * ( v(k,j,0) - v_m_l(k,j,0) ) / denom
452                   IF ( c_v(k,j) < 0.0 )  THEN
453                      c_v(k,j) = 0.0
454                   ELSEIF ( c_v(k,j) > c_max )  THEN
455                      c_v(k,j) = c_max
456                   ENDIF
457                ELSE
458                   c_v(k,j) = c_max
459                ENDIF
460
461                denom = w_m_l(k,j,0) - w_m_l(k,j,1)
462
463                IF ( denom /= 0.0 )  THEN
464                   c_w(k,j) = -c_max * ( w(k,j,0) - w_m_l(k,j,0) ) / denom
465                   IF ( c_w(k,j) < 0.0 )  THEN
466                      c_w(k,j) = 0.0
467                   ELSEIF ( c_w(k,j) > c_max )  THEN
468                      c_w(k,j) = c_max
469                   ENDIF
470                ELSE
471                   c_w(k,j) = c_max
472                ENDIF
473
474!
475!--             Swap timelevels for the next timestep
476                u_m_l(k,j,:) = u(k,j,1:2)
477                v_m_l(k,j,:) = v(k,j,0:1)
478                w_m_l(k,j,:) = w(k,j,0:1)
479
480             ENDIF
481
482!
483!--          Calculate the new velocities
484             u_p(k,j,0)  = u(k,j,0)  - dt_3d * tsc(2) * c_u(k,j) * &
485                                       ( u(k,j,0) - u(k,j,1) ) * ddx
486
487             v_p(k,j,-1) = v(k,j,-1) - dt_3d * tsc(2) * c_v(k,j) * &
488                                       ( v(k,j,-1) - v(k,j,0) ) * ddx
489
490             w_p(k,j,-1) = w(k,j,-1) - dt_3d * tsc(2) * c_w(k,j) * &
491                                       ( w(k,j,-1) - w(k,j,0) ) * ddx
492
493          ENDDO
494       ENDDO
495
496!
497!--    Bottom boundary at the outflow
498       IF ( ibc_uv_b == 0 )  THEN
499          u_p(nzb,:,0)  = 0.0 
500          v_p(nzb,:,-1) = 0.0
501       ELSE                   
502          u_p(nzb,:,0)  =  u_p(nzb+1,:,0)
503          v_p(nzb,:,-1) =  v_p(nzb+1,:,-1)
504       ENDIF
505       w_p(nzb,:,-1) = 0.0
506
507!
508!--    Top boundary at the outflow
509       IF ( ibc_uv_t == 0 )  THEN
510          u_p(nzt+1,:,-1) = ug(nzt+1)
511          v_p(nzt+1,:,-1) = vg(nzt+1)
512       ELSE
513          u_p(nzt+1,:,-1) = u_p(nzt,:,-1)
514          v_p(nzt+1,:,-1) = v_p(nzt,:,-1)
515       ENDIF
516       w_p(nzt:nzt+1,:,-1) = 0.0
517
518    ENDIF
519
520    IF ( outflow_r )  THEN
521
522       c_max = dx / dt_3d
523
524       DO j = nysg, nyng
525          DO k = nzb+1, nzt+1
526
527!
528!--          Calculate the phase speeds for u,v, and w. In case of using a
529!--          Runge-Kutta scheme, do this for the first substep only and then
530!--          reuse this values for the further substeps.
531             IF ( intermediate_timestep_count == 1 )  THEN
532
533                denom = u_m_r(k,j,nx) - u_m_r(k,j,nx-1)
534
535                IF ( denom /= 0.0 )  THEN
536                   c_u(k,j) = -c_max * ( u(k,j,nx) - u_m_r(k,j,nx) ) / denom
537                   IF ( c_u(k,j) < 0.0 )  THEN
538                      c_u(k,j) = 0.0
539                   ELSEIF ( c_u(k,j) > c_max )  THEN
540                      c_u(k,j) = c_max
541                   ENDIF
542                ELSE
543                   c_u(k,j) = c_max
544                ENDIF
545
546                denom = v_m_r(k,j,nx) - v_m_r(k,j,nx-1)
547
548                IF ( denom /= 0.0 )  THEN
549                   c_v(k,j) = -c_max * ( v(k,j,nx) - v_m_r(k,j,nx) ) / denom
550                   IF ( c_v(k,j) < 0.0 )  THEN
551                      c_v(k,j) = 0.0
552                   ELSEIF ( c_v(k,j) > c_max )  THEN
553                      c_v(k,j) = c_max
554                   ENDIF
555                ELSE
556                   c_v(k,j) = c_max
557                ENDIF
558
559                denom = w_m_r(k,j,nx) - w_m_r(k,j,nx-1)
560
561                IF ( denom /= 0.0 )  THEN
562                   c_w(k,j) = -c_max * ( w(k,j,nx) - w_m_r(k,j,nx) ) / denom
563                   IF ( c_w(k,j) < 0.0 )  THEN
564                      c_w(k,j) = 0.0
565                   ELSEIF ( c_w(k,j) > c_max )  THEN
566                      c_w(k,j) = c_max
567                   ENDIF
568                ELSE
569                   c_w(k,j) = c_max
570                ENDIF
571
572!
573!--             Swap timelevels for the next timestep
574                u_m_r(k,j,:) = u(k,j,nx-1:nx)
575                v_m_r(k,j,:) = v(k,j,nx-1:nx)
576                w_m_r(k,j,:) = w(k,j,nx-1:nx)
577
578             ENDIF
579
580!
581!--          Calculate the new velocities
582             u_p(k,j,nx+1) = u(k,j,nx+1) - dt_3d * tsc(2) * c_u(k,j) * &
583                                           ( u(k,j,nx+1) - u(k,j,nx) ) * ddx
584
585             v_p(k,j,nx+1) = v(k,j,nx+1) - dt_3d * tsc(2) * c_v(k,j) * &
586                                           ( v(k,j,nx+1) - v(k,j,nx) ) * ddx
587
588             w_p(k,j,nx+1) = w(k,j,nx+1) - dt_3d * tsc(2) * c_w(k,j) * &
589                                           ( w(k,j,nx+1) - w(k,j,nx) ) * ddx
590
591          ENDDO
592       ENDDO
593
594!
595!--    Bottom boundary at the outflow
596       IF ( ibc_uv_b == 0 )  THEN
597          u_p(nzb,:,nx+1) = 0.0
598          v_p(nzb,:,nx+1) = 0.0 
599       ELSE                   
600          u_p(nzb,:,nx+1) =  u_p(nzb+1,:,nx+1)
601          v_p(nzb,:,nx+1) =  v_p(nzb+1,:,nx+1)
602       ENDIF
603       w_p(nzb,:,nx+1) = 0.0
604
605!
606!--    Top boundary at the outflow
607       IF ( ibc_uv_t == 0 )  THEN
608          u_p(nzt+1,:,nx+1) = ug(nzt+1)
609          v_p(nzt+1,:,nx+1) = vg(nzt+1)
610       ELSE
611          u_p(nzt+1,:,nx+1) = u_p(nzt,:,nx+1)
612          v_p(nzt+1,:,nx+1) = v_p(nzt,:,nx+1)
613       ENDIF
614       w(nzt:nzt+1,:,nx+1) = 0.0
615
616    ENDIF
617
618 
619 END SUBROUTINE boundary_conds
Note: See TracBrowser for help on using the repository browser.