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

Last change on this file since 667 was 667, checked in by suehring, 13 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
RevLine 
[1]1 SUBROUTINE boundary_conds( range )
2
3!------------------------------------------------------------------------------!
[484]4! Current revisions:
[1]5! -----------------
[667]6!
7! nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng
8!
[106]9!
[667]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
[1]14! Former revisions:
15! -----------------
[3]16! $Id: boundary_conds.f90 667 2010-12-23 12:06:00Z suehring $
[39]17!
[110]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!
[98]24! 95 2007-06-02 16:48:38Z raasch
25! Boundary conditions for salinity added
26!
[77]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!
[39]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!
[3]37! RCS Log replace by Id keyword, revision history cleaned up
38!
[1]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
[106]73    REAL    ::  c_max, denom
[1]74
[73]75
[1]76    IF ( range == 'main')  THEN
77!
[667]78!--    Bottom boundary
79       IF ( ibc_uv_b == 1 )  THEN
[73]80          u_p(nzb,:,:) = u_p(nzb+1,:,:)
81          v_p(nzb,:,:) = v_p(nzb+1,:,:)
[1]82       ENDIF
[667]83       DO  i = nxlg, nxrg
84          DO  j = nysg, nyng
[73]85             w_p(nzb_w_inner(j,i),j,i) = 0.0
[1]86          ENDDO
87       ENDDO
88
89!
90!--    Top boundary
91       IF ( ibc_uv_t == 0 )  THEN
[667]92           u_p(nzt+1,:,:) = ug(nzt+1)
93           v_p(nzt+1,:,:) = vg(nzt+1)
[1]94       ELSE
[667]95           u_p(nzt+1,:,:) = u_p(nzt,:,:)
96           v_p(nzt+1,:,:) = v_p(nzt,:,:)
[1]97       ENDIF
[73]98       w_p(nzt:nzt+1,:,:) = 0.0  ! nzt is not a prognostic level (but cf. pres)
[1]99
100!
[102]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.
[1]104       IF ( ibc_pt_b == 0 )  THEN
[667]105          DO  i = nxlg, nxrg
106             DO  j = nysg, nyng
[73]107                pt_p(nzb_s_inner(j,i),j,i) = pt(nzb_s_inner(j,i),j,i)
[1]108             ENDDO
[73]109          ENDDO
[102]110       ELSEIF ( ibc_pt_b == 1 )  THEN
[667]111          DO  i = nxlg, nxrg
112             DO  j = nysg, nyng
[73]113                pt_p(nzb_s_inner(j,i),j,i) = pt_p(nzb_s_inner(j,i)+1,j,i)
[1]114             ENDDO
115          ENDDO
116       ENDIF
117
118!
119!--    Temperature at top boundary
[19]120       IF ( ibc_pt_t == 0 )  THEN
[667]121           pt_p(nzt+1,:,:) = pt(nzt+1,:,:)
[19]122       ELSEIF ( ibc_pt_t == 1 )  THEN
[667]123           pt_p(nzt+1,:,:) = pt_p(nzt,:,:)
[19]124       ELSEIF ( ibc_pt_t == 2 )  THEN
[667]125           pt_p(nzt+1,:,:) = pt_p(nzt,:,:)   + bc_pt_t_val * dzu(nzt+1)
[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
[667]132          DO  i = nxlg, nxrg
133             DO  j = nysg, nyng
[73]134                e_p(nzb_s_inner(j,i),j,i) = e_p(nzb_s_inner(j,i)+1,j,i)
[1]135             ENDDO
136          ENDDO
[73]137          e_p(nzt+1,:,:) = e_p(nzt,:,:)
[1]138       ENDIF
139
140!
[95]141!--    Boundary conditions for salinity
142       IF ( ocean )  THEN
143!
144!--       Bottom boundary: Neumann condition because salinity flux is always
145!--       given
[667]146          DO  i = nxlg, nxrg
147             DO  j = nysg, nyng
[95]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
[667]155              sa_p(nzt+1,:,:) = sa(nzt+1,:,:)
[95]156          ELSEIF ( ibc_sa_t == 1 )  THEN
[667]157              sa_p(nzt+1,:,:) = sa_p(nzt,:,:)
[95]158          ENDIF
159
160       ENDIF
161
162!
[1]163!--    Boundary conditions for total water content or scalar,
[95]164!--    bottom and top boundary (see also temperature)
[75]165       IF ( humidity  .OR.  passive_scalar )  THEN
[1]166!
[75]167!--       Surface conditions for constant_humidity_flux
[1]168          IF ( ibc_q_b == 0 ) THEN
[667]169             DO  i = nxlg, nxrg
170                DO  j = nysg, nyng
[73]171                   q_p(nzb_s_inner(j,i),j,i) = q(nzb_s_inner(j,i),j,i)
[1]172                ENDDO
[73]173             ENDDO
[1]174          ELSE
[667]175             DO  i = nxlg, nxrg
176                DO  j = nysg, nyng
[73]177                   q_p(nzb_s_inner(j,i),j,i) = q_p(nzb_s_inner(j,i)+1,j,i)
[1]178                ENDDO
179             ENDDO
180          ENDIF
181!
182!--       Top boundary
[73]183          q_p(nzt+1,:,:) = q_p(nzt,:,:)   + bc_q_t_val * dzu(nzt+1)
[667]184
185
[1]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
[73]193          v_p(:,nys,:) = v_p(:,nys-1,:)
[1]194       ELSEIF ( inflow_n ) THEN
[75]195          v_p(:,nyn,:) = v_p(:,nyn+1,:)
[1]196       ELSEIF ( inflow_l ) THEN
[73]197          u_p(:,:,nxl) = u_p(:,:,nxl-1)
[1]198       ELSEIF ( inflow_r ) THEN
[75]199          u_p(:,:,nxr) = u_p(:,:,nxr+1)
[1]200       ENDIF
201
202!
203!--    Lateral boundary conditions for scalar quantities at the outflow
204       IF ( outflow_s )  THEN
[73]205          pt_p(:,nys-1,:)     = pt_p(:,nys,:)
206          IF ( .NOT. constant_diffusion     )  e_p(:,nys-1,:) = e_p(:,nys,:)
[75]207          IF ( humidity .OR. passive_scalar )  q_p(:,nys-1,:) = q_p(:,nys,:)
[1]208       ELSEIF ( outflow_n )  THEN
[73]209          pt_p(:,nyn+1,:)     = pt_p(:,nyn,:)
210          IF ( .NOT. constant_diffusion     )  e_p(:,nyn+1,:) = e_p(:,nyn,:)
[75]211          IF ( humidity .OR. passive_scalar )  q_p(:,nyn+1,:) = q_p(:,nyn,:)
[1]212       ELSEIF ( outflow_l )  THEN
[73]213          pt_p(:,:,nxl-1)     = pt_p(:,:,nxl)
214          IF ( .NOT. constant_diffusion     )  e_p(:,:,nxl-1) = e_p(:,:,nxl)
[75]215          IF ( humidity .OR. passive_scalar )  q_p(:,:,nxl-1) = q_p(:,:,nxl)
[1]216       ELSEIF ( outflow_r )  THEN
[73]217          pt_p(:,:,nxr+1)     = pt_p(:,:,nxr)
218          IF ( .NOT. constant_diffusion     )  e_p(:,:,nxr+1) = e_p(:,:,nxr)
[75]219          IF ( humidity .OR. passive_scalar )  q_p(:,:,nxr+1) = q_p(:,:,nxr)
[1]220       ENDIF
221
222    ENDIF
223
224!
[75]225!-- Radiation boundary condition for the velocities at the respective outflow
[106]226    IF ( outflow_s )  THEN
[75]227
228       c_max = dy / dt_3d
229
[667]230       DO i = nxlg, nxrg
[75]231          DO k = nzb+1, nzt+1
232
233!
[106]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
[75]238
[106]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
[75]250                ENDIF
251
[106]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
[75]263                ENDIF
264
[106]265                denom = w_m_s(k,0,i) - w_m_s(k,1,i)
[75]266
[106]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
[75]276                ENDIF
[106]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
[75]284             ENDIF
285
286!
287!--          Calculate the new velocities
[106]288             u_p(k,-1,i) = u(k,-1,i) - dt_3d * tsc(2) * c_u(k,i) * &
[75]289                                       ( u(k,-1,i) - u(k,0,i) ) * ddy
290
[107]291             v_p(k,0,i)  = v(k,0,i)  - dt_3d * tsc(2) * c_v(k,i) * &
[106]292                                       ( v(k,0,i) - v(k,1,i) ) * ddy
[75]293
[106]294             w_p(k,-1,i) = w(k,-1,i) - dt_3d * tsc(2) * c_w(k,i) * &
[75]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
[667]303          u_p(nzb,-1,:) = 0.0 
304          v_p(nzb,0,:)  = 0.0 
[75]305       ELSE                   
306          u_p(nzb,-1,:) =  u_p(nzb+1,-1,:)
[106]307          v_p(nzb,0,:)  =  v_p(nzb+1,0,:)
[73]308       ENDIF
[106]309       w_p(nzb,-1,:) = 0.0
[73]310
[75]311!
312!--    Top boundary at the outflow
313       IF ( ibc_uv_t == 0 )  THEN
314          u_p(nzt+1,-1,:) = ug(nzt+1)
[106]315          v_p(nzt+1,0,:)  = vg(nzt+1)
[75]316       ELSE
317          u_p(nzt+1,-1,:) = u(nzt,-1,:)
[106]318          v_p(nzt+1,0,:)  = v(nzt,0,:)
[75]319       ENDIF
320       w_p(nzt:nzt+1,-1,:) = 0.0
[73]321
[75]322    ENDIF
[73]323
[106]324    IF ( outflow_n )  THEN
[73]325
[75]326       c_max = dy / dt_3d
327
[667]328       DO i = nxlg, nxrg
[75]329          DO k = nzb+1, nzt+1
330
[1]331!
[106]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
[73]336
[106]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
[73]348                ENDIF
349
[106]350                denom = v_m_n(k,ny,i) - v_m_n(k,ny-1,i)
[73]351
[106]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
[73]361                ENDIF
362
[106]363                denom = w_m_n(k,ny,i) - w_m_n(k,ny-1,i)
[73]364
[106]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
[73]374                ENDIF
[106]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
[75]382             ENDIF
[73]383
384!
[75]385!--          Calculate the new velocities
[106]386             u_p(k,ny+1,i) = u(k,ny+1,i) - dt_3d * tsc(2) * c_u(k,i) * &
[75]387                                           ( u(k,ny+1,i) - u(k,ny,i) ) * ddy
[73]388
[106]389             v_p(k,ny+1,i) = v(k,ny+1,i) - dt_3d * tsc(2) * c_v(k,i) * &
[75]390                                           ( v(k,ny+1,i) - v(k,ny,i) ) * ddy
[73]391
[106]392             w_p(k,ny+1,i) = w(k,ny+1,i) - dt_3d * tsc(2) * c_w(k,i) * &
[75]393                                           ( w(k,ny+1,i) - w(k,ny,i) ) * ddy
[73]394
[1]395          ENDDO
[75]396       ENDDO
[1]397
398!
[75]399!--    Bottom boundary at the outflow
400       IF ( ibc_uv_b == 0 )  THEN
[667]401          u_p(nzb,ny+1,:) = 0.0 
402          v_p(nzb,ny+1,:) = 0.0   
[75]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
[73]408
409!
[75]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,:)
[1]417       ENDIF
[75]418       w_p(nzt:nzt+1,ny+1,:) = 0.0
[1]419
[75]420    ENDIF
421
[106]422    IF ( outflow_l )  THEN
[75]423
424       c_max = dx / dt_3d
425
[667]426       DO j = nysg, nyng
[75]427          DO k = nzb+1, nzt+1
428
[1]429!
[106]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
[75]434
[106]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
[107]439                   IF ( c_u(k,j) < 0.0 )  THEN
[106]440                      c_u(k,j) = 0.0
[107]441                   ELSEIF ( c_u(k,j) > c_max )  THEN
442                      c_u(k,j) = c_max
[106]443                   ENDIF
444                ELSE
[107]445                   c_u(k,j) = c_max
[75]446                ENDIF
447
[106]448                denom = v_m_l(k,j,0) - v_m_l(k,j,1)
[75]449
[106]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
[75]459                ENDIF
460
[106]461                denom = w_m_l(k,j,0) - w_m_l(k,j,1)
[75]462
[106]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
[75]472                ENDIF
[106]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
[75]480             ENDIF
481
[73]482!
[75]483!--          Calculate the new velocities
[106]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
[75]486
[106]487             v_p(k,j,-1) = v(k,j,-1) - dt_3d * tsc(2) * c_v(k,j) * &
[75]488                                       ( v(k,j,-1) - v(k,j,0) ) * ddx
489
[106]490             w_p(k,j,-1) = w(k,j,-1) - dt_3d * tsc(2) * c_w(k,j) * &
[75]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
[667]499          u_p(nzb,:,0)  = 0.0 
500          v_p(nzb,:,-1) = 0.0
[75]501       ELSE                   
[667]502          u_p(nzb,:,0)  =  u_p(nzb+1,:,0)
[75]503          v_p(nzb,:,-1) =  v_p(nzb+1,:,-1)
[1]504       ENDIF
[75]505       w_p(nzb,:,-1) = 0.0
[1]506
[75]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
[73]517
[75]518    ENDIF
[73]519
[106]520    IF ( outflow_r )  THEN
[73]521
[75]522       c_max = dx / dt_3d
523
[667]524       DO j = nysg, nyng
[75]525          DO k = nzb+1, nzt+1
526
[1]527!
[106]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
[73]532
[106]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
[73]544                ENDIF
545
[106]546                denom = v_m_r(k,j,nx) - v_m_r(k,j,nx-1)
[73]547
[106]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
[73]557                ENDIF
558
[106]559                denom = w_m_r(k,j,nx) - w_m_r(k,j,nx-1)
[73]560
[106]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
[73]570                ENDIF
[106]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
[75]578             ENDIF
[73]579
580!
[75]581!--          Calculate the new velocities
[106]582             u_p(k,j,nx+1) = u(k,j,nx+1) - dt_3d * tsc(2) * c_u(k,j) * &
[75]583                                           ( u(k,j,nx+1) - u(k,j,nx) ) * ddx
[73]584
[106]585             v_p(k,j,nx+1) = v(k,j,nx+1) - dt_3d * tsc(2) * c_v(k,j) * &
[75]586                                           ( v(k,j,nx+1) - v(k,j,nx) ) * ddx
[73]587
[106]588             w_p(k,j,nx+1) = w(k,j,nx+1) - dt_3d * tsc(2) * c_w(k,j) * &
[75]589                                           ( w(k,j,nx+1) - w(k,j,nx) ) * ddx
[73]590
591          ENDDO
[75]592       ENDDO
[73]593
594!
[75]595!--    Bottom boundary at the outflow
596       IF ( ibc_uv_b == 0 )  THEN
[667]597          u_p(nzb,:,nx+1) = 0.0
598          v_p(nzb,:,nx+1) = 0.0 
[75]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
[73]604
605!
[75]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)
[1]613       ENDIF
[75]614       w(nzt:nzt+1,:,nx+1) = 0.0
[1]615
616    ENDIF
617
618 
619 END SUBROUTINE boundary_conds
Note: See TracBrowser for help on using the repository browser.