source: palm/trunk/SOURCE/wall_fluxes.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: 20.8 KB
Line 
1 MODULE wall_fluxes_mod
2!------------------------------------------------------------------------------!
3! Current revisions:
4! -----------------
5!
6!
7! Former revisions:
8! -----------------
9! $Id: wall_fluxes.f90 667 2010-12-23 12:06:00Z suehring $
10!
11! 187 2008-08-06 16:25:09Z letzel
12! Bugfix: Modification of the evaluation of the vertical turbulent momentum
13! fluxes u'w' and v'w (see prandtl_fluxes), this requires the calculation of
14! us_wall (and vel_total, u_i, v_i, ws) also in wall_fluxes_e.
15! Bugfix: change definition of us_wall from 1D to 2D
16! Bugfix: storage of rifs to rifs_wall in wall_fluxes_e removed
17! Change: add 'minus' sign to fluxes produced by subroutine wall_fluxes_e for
18! consistency with subroutine wall_fluxes
19! Change: Modification of the integrated version of the profile function for
20! momentum for unstable stratification
21!
22! Initial version (2007/03/07)
23!
24! Description:
25! ------------
26! Calculates momentum fluxes at vertical walls assuming Monin-Obukhov
27! similarity.
28! Indices: usvs a=1, vsus b=1, wsvs c1=1, wsus c2=1 (other=0).
29! The all-gridpoint version of wall_fluxes_e is not used so far, because
30! it gives slightly different results from the ij-version for some unknown
31! reason.
32!------------------------------------------------------------------------------!
33    PRIVATE
34    PUBLIC wall_fluxes, wall_fluxes_e
35   
36    INTERFACE wall_fluxes
37       MODULE PROCEDURE wall_fluxes
38       MODULE PROCEDURE wall_fluxes_ij
39    END INTERFACE wall_fluxes
40   
41    INTERFACE wall_fluxes_e
42       MODULE PROCEDURE wall_fluxes_e
43       MODULE PROCEDURE wall_fluxes_e_ij
44    END INTERFACE wall_fluxes_e
45 
46 CONTAINS
47
48!------------------------------------------------------------------------------!
49! Call for all grid points
50!------------------------------------------------------------------------------!
51    SUBROUTINE wall_fluxes( wall_flux, a, b, c1, c2, nzb_uvw_inner, &
52                            nzb_uvw_outer, wall )
53
54       USE arrays_3d
55       USE control_parameters
56       USE grid_variables
57       USE indices
58       USE statistics
59
60       IMPLICIT NONE
61
62       INTEGER ::  i, j, k, wall_index
63
64       INTEGER, DIMENSION(nysg:nyng,nxlg:nxrg) ::  nzb_uvw_inner, &
65                                                       nzb_uvw_outer
66       REAL ::  a, b, c1, c2, h1, h2, zp
67       REAL ::  pts, pt_i, rifs, u_i, v_i, us_wall, vel_total, ws, wspts
68
69       REAL, DIMENSION(nysg:nyng,nxlg:nxrg)   ::  wall
70       REAL, DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) ::  wall_flux
71
72
73       zp         = 0.5 * ( (a+c1) * dy + (b+c2) * dx )
74       wall_flux  = 0.0
75       wall_index = NINT( a+ 2*b + 3*c1 + 4*c2 )
76
77       DO  i = nxl, nxr
78          DO  j = nys, nyn
79
80             IF ( wall(j,i) /= 0.0 )  THEN
81!
82!--             All subsequent variables are computed for the respective
83!--             location where the respective flux is defined.
84                DO  k = nzb_uvw_inner(j,i)+1, nzb_uvw_outer(j,i)
85
86!
87!--                (1) Compute rifs, u_i, v_i, ws, pt' and w'pt'
88                   rifs  = rif_wall(k,j,i,wall_index)
89
90                   u_i   = a * u(k,j,i) + c1 * 0.25 * &
91                           ( u(k+1,j,i+1) + u(k+1,j,i) + u(k,j,i+1) + u(k,j,i) )
92
93                   v_i   = b * v(k,j,i) + c2 * 0.25 * &
94                           ( v(k+1,j+1,i) + v(k+1,j,i) + v(k,j+1,i) + v(k,j,i) )
95
96                   ws    = ( c1 + c2 ) * w(k,j,i) + 0.25 * (                   &
97                     a * ( w(k-1,j,i-1) + w(k-1,j,i) + w(k,j,i-1) + w(k,j,i) ) &
98                   + b * ( w(k-1,j-1,i) + w(k-1,j,i) + w(k,j-1,i) + w(k,j,i) ) &
99                                                           )
100                   pt_i  = 0.5 * ( pt(k,j,i) + a *  pt(k,j,i-1) + &
101                                   b * pt(k,j-1,i) + ( c1 + c2 ) * pt(k+1,j,i) )
102
103                   pts   = pt_i - hom(k,1,4,0)
104                   wspts = ws * pts
105
106!
107!--                (2) Compute wall-parallel absolute velocity vel_total
108                   vel_total = SQRT( ws**2 + (a+c1) * u_i**2 + (b+c2) * v_i**2 )
109
110!
111!--                (3) Compute wall friction velocity us_wall
112                   IF ( rifs >= 0.0 )  THEN
113
114!
115!--                   Stable stratification (and neutral)
116                      us_wall = kappa * vel_total / ( LOG( zp / z0(j,i) ) +    &
117                                            5.0 * rifs * ( zp - z0(j,i) ) / zp &
118                                                    )
119                   ELSE
120
121!
122!--                   Unstable stratification
123                      h1 = SQRT( SQRT( 1.0 - 16.0 * rifs ) )
124                      h2 = SQRT( SQRT( 1.0 - 16.0 * rifs * z0(j,i) / zp ) )
125
126                      us_wall = kappa * vel_total / (                          &
127                           LOG( zp / z0(j,i) ) -                               &
128                           LOG( ( 1.0 + h1 )**2 * ( 1.0 + h1**2 ) / (          &
129                                ( 1.0 + h2 )**2 * ( 1.0 + h2**2 )   ) ) +      &
130                                2.0 * ( ATAN( h1 ) - ATAN( h2 ) )              &
131                                                    )
132                   ENDIF
133
134!
135!--                (4) Compute zp/L (corresponds to neutral Richardson flux
136!--                    number rifs)
137                   rifs = -1.0 * zp * kappa * g * wspts / ( pt_i * &
138                                                        ( us_wall**3 + 1E-30 ) )
139
140!
141!--                Limit the value range of the Richardson numbers.
142!--                This is necessary for very small velocities (u,w --> 0),
143!--                because the absolute value of rif can then become very
144!--                large, which in consequence would result in very large
145!--                shear stresses and very small momentum fluxes (both are
146!--                generally unrealistic).
147                   IF ( rifs < rif_min )  rifs = rif_min
148                   IF ( rifs > rif_max )  rifs = rif_max
149
150!
151!--                (5) Compute wall_flux (u'v', v'u', w'v', or w'u')
152                   IF ( rifs >= 0.0 )  THEN
153
154!
155!--                   Stable stratification (and neutral)
156                      wall_flux(k,j,i) = kappa *                               &
157                              ( a*u(k,j,i) + b*v(k,j,i) + (c1+c2)*w(k,j,i) ) / &
158                              (  LOG( zp / z0(j,i) ) +                         &
159                                 5.0 * rifs * ( zp - z0(j,i) ) / zp            &
160                              )
161                   ELSE
162
163!
164!--                   Unstable stratification
165                      h1 = SQRT( SQRT( 1.0 - 16.0 * rifs ) )
166                      h2 = SQRT( SQRT( 1.0 - 16.0 * rifs * z0(j,i) / zp ) )
167
168                      wall_flux(k,j,i) = kappa *                               &
169                           ( a*u(k,j,i) + b*v(k,j,i) + (c1+c2)*w(k,j,i) ) / (  &
170                           LOG( zp / z0(j,i) ) -                               &
171                           LOG( ( 1.0 + h1 )**2 * ( 1.0 + h1**2 ) / (          &
172                                ( 1.0 + h2 )**2 * ( 1.0 + h2**2 )   ) ) +      &
173                                2.0 * ( ATAN( h1 ) - ATAN( h2 ) )              &
174                                                                            )
175                   ENDIF
176                   wall_flux(k,j,i) = -wall_flux(k,j,i) * us_wall
177
178!
179!--                store rifs for next time step
180                   rif_wall(k,j,i,wall_index) = rifs
181
182                ENDDO
183
184             ENDIF
185
186          ENDDO
187       ENDDO
188
189    END SUBROUTINE wall_fluxes
190
191
192
193!------------------------------------------------------------------------------!
194! Call for all grid point i,j
195!------------------------------------------------------------------------------!
196    SUBROUTINE wall_fluxes_ij( i, j, nzb_w, nzt_w, wall_flux, a, b, c1, c2 )
197
198       USE arrays_3d
199       USE control_parameters
200       USE grid_variables
201       USE indices
202       USE statistics
203
204       IMPLICIT NONE
205
206       INTEGER ::  i, j, k, nzb_w, nzt_w, wall_index
207       REAL    ::  a, b, c1, c2, h1, h2, zp
208
209       REAL ::  pts, pt_i, rifs, u_i, v_i, us_wall, vel_total, ws, wspts
210
211       REAL, DIMENSION(nzb:nzt+1) ::  wall_flux
212
213
214       zp         = 0.5 * ( (a+c1) * dy + (b+c2) * dx )
215       wall_flux  = 0.0
216       wall_index = NINT( a+ 2*b + 3*c1 + 4*c2 )
217
218!
219!--    All subsequent variables are computed for the respective location where
220!--    the respective flux is defined.
221       DO  k = nzb_w, nzt_w
222
223!
224!--       (1) Compute rifs, u_i, v_i, ws, pt' and w'pt'
225          rifs  = rif_wall(k,j,i,wall_index)
226
227          u_i   = a * u(k,j,i) + c1 * 0.25 * &
228                  ( u(k+1,j,i+1) + u(k+1,j,i) + u(k,j,i+1) + u(k,j,i) )
229
230          v_i   = b * v(k,j,i) + c2 * 0.25 * &
231                  ( v(k+1,j+1,i) + v(k+1,j,i) + v(k,j+1,i) + v(k,j,i) )
232
233          ws    = ( c1 + c2 ) * w(k,j,i) + 0.25 * (                            &
234                     a * ( w(k-1,j,i-1) + w(k-1,j,i) + w(k,j,i-1) + w(k,j,i) ) &
235                   + b * ( w(k-1,j-1,i) + w(k-1,j,i) + w(k,j-1,i) + w(k,j,i) ) &
236                                                  )
237          pt_i  = 0.5 * ( pt(k,j,i) + a *  pt(k,j,i-1) + b * pt(k,j-1,i)  &
238                          + ( c1 + c2 ) * pt(k+1,j,i) )
239
240          pts   = pt_i - hom(k,1,4,0)
241          wspts = ws * pts
242
243!
244!--       (2) Compute wall-parallel absolute velocity vel_total
245          vel_total = SQRT( ws**2 + ( a+c1 ) * u_i**2 + ( b+c2 ) * v_i**2 )
246
247!
248!--       (3) Compute wall friction velocity us_wall
249          IF ( rifs >= 0.0 )  THEN
250
251!
252!--          Stable stratification (and neutral)
253             us_wall = kappa * vel_total / ( LOG( zp / z0(j,i) ) +             &
254                                            5.0 * rifs * ( zp - z0(j,i) ) / zp &
255                                           )
256          ELSE
257
258!
259!--          Unstable stratification
260             h1 = SQRT( SQRT( 1.0 - 16.0 * rifs ) )
261             h2 = SQRT( SQRT( 1.0 - 16.0 * rifs * z0(j,i) / zp ) )
262
263             us_wall = kappa * vel_total / (                          &
264                  LOG( zp / z0(j,i) ) -                               &
265                  LOG( ( 1.0 + h1 )**2 * ( 1.0 + h1**2 ) / (          &
266                       ( 1.0 + h2 )**2 * ( 1.0 + h2**2 )   ) ) +      &
267                       2.0 * ( ATAN( h1 ) - ATAN( h2 ) )              &
268                                           )
269          ENDIF
270
271!
272!--       (4) Compute zp/L (corresponds to neutral Richardson flux number
273!--           rifs)
274          rifs = -1.0 * zp * kappa * g * wspts / ( pt_i * (us_wall**3 + 1E-30) )
275
276!
277!--       Limit the value range of the Richardson numbers.
278!--       This is necessary for very small velocities (u,w --> 0), because
279!--       the absolute value of rif can then become very large, which in
280!--       consequence would result in very large shear stresses and very
281!--       small momentum fluxes (both are generally unrealistic).
282          IF ( rifs < rif_min )  rifs = rif_min
283          IF ( rifs > rif_max )  rifs = rif_max
284
285!
286!--       (5) Compute wall_flux (u'v', v'u', w'v', or w'u')
287          IF ( rifs >= 0.0 )  THEN
288
289!
290!--          Stable stratification (and neutral)
291             wall_flux(k) = kappa *                                          &
292                            ( a*u(k,j,i) + b*v(k,j,i) + (c1+c2)*w(k,j,i) ) / &
293                            (  LOG( zp / z0(j,i) ) +                         &
294                               5.0 * rifs * ( zp - z0(j,i) ) / zp            &
295                            )
296          ELSE
297
298!
299!--          Unstable stratification
300             h1 = SQRT( SQRT( 1.0 - 16.0 * rifs ) )
301             h2 = SQRT( SQRT( 1.0 - 16.0 * rifs * z0(j,i) / zp ) )
302
303             wall_flux(k) = kappa *                               &
304                  ( a*u(k,j,i) + b*v(k,j,i) + (c1+c2)*w(k,j,i) ) / (  &
305                  LOG( zp / z0(j,i) ) -                               &
306                  LOG( ( 1.0 + h1 )**2 * ( 1.0 + h1**2 ) / (          &
307                       ( 1.0 + h2 )**2 * ( 1.0 + h2**2 )   ) ) +      &
308                       2.0 * ( ATAN( h1 ) - ATAN( h2 ) )              &
309                                                                   )
310          ENDIF
311          wall_flux(k) = -wall_flux(k) * us_wall
312
313!
314!--       store rifs for next time step
315          rif_wall(k,j,i,wall_index) = rifs
316
317       ENDDO
318
319    END SUBROUTINE wall_fluxes_ij
320
321
322
323!------------------------------------------------------------------------------!
324! Call for all grid points
325!------------------------------------------------------------------------------!
326    SUBROUTINE wall_fluxes_e( wall_flux, a, b, c1, c2, wall )
327
328!------------------------------------------------------------------------------!
329! Description:
330! ------------
331! Calculates momentum fluxes at vertical walls for routine production_e
332! assuming Monin-Obukhov similarity.
333! Indices: usvs a=1, vsus b=1, wsvs c1=1, wsus c2=1 (other=0).
334!------------------------------------------------------------------------------!
335
336       USE arrays_3d
337       USE control_parameters
338       USE grid_variables
339       USE indices
340       USE statistics
341
342       IMPLICIT NONE
343
344       INTEGER ::  i, j, k, kk, wall_index
345       REAL    ::  a, b, c1, c2, h1, h2, u_i, v_i, us_wall, vel_total, vel_zp, &
346                   ws, zp
347
348       REAL ::  rifs
349
350       REAL, DIMENSION(nysg:nyng,nxlg:nxrg)   ::  wall
351       REAL, DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) ::  wall_flux
352
353
354       zp         = 0.5 * ( (a+c1) * dy + (b+c2) * dx )
355       wall_flux  = 0.0
356       wall_index = NINT( a+ 2*b + 3*c1 + 4*c2 )
357
358       DO  i = nxl, nxr
359          DO  j = nys, nyn
360
361             IF ( wall(j,i) /= 0.0 )  THEN
362!
363!--             All subsequent variables are computed for scalar locations.
364                DO  k = nzb_diff_s_inner(j,i)-1, nzb_diff_s_outer(j,i)-2
365!
366!--                (1) Compute rifs, u_i, v_i, and ws
367                   IF ( k == nzb_diff_s_inner(j,i)-1 )  THEN
368                      kk = nzb_diff_s_inner(j,i)-1
369                   ELSE
370                      kk = k-1
371                   ENDIF
372                   rifs  = 0.5 * ( rif_wall(k,j,i,wall_index) +                &
373                          a * rif_wall(k,j,i+1,1) +  b * rif_wall(k,j+1,i,2) + &
374                          c1 * rif_wall(kk,j,i,3) + c2 * rif_wall(kk,j,i,4)    &
375                                 )
376
377                   u_i   = 0.5 * ( u(k,j,i) + u(k,j,i+1) )
378                   v_i   = 0.5 * ( v(k,j,i) + v(k,j+1,i) )
379                   ws    = 0.5 * ( w(k,j,i) + w(k-1,j,i) )
380!
381!--                (2) Compute wall-parallel absolute velocity vel_total and
382!--                interpolate appropriate velocity component vel_zp.
383                   vel_total = SQRT( ws**2 + (a+c1) * u_i**2 + (b+c2) * v_i**2 )
384                   vel_zp = 0.5 * ( a * u_i + b * v_i + (c1+c2) * ws )
385!
386!--                (3) Compute wall friction velocity us_wall
387                   IF ( rifs >= 0.0 )  THEN
388
389!
390!--                   Stable stratification (and neutral)
391                      us_wall = kappa * vel_total / ( LOG( zp / z0(j,i) ) +    &
392                                            5.0 * rifs * ( zp - z0(j,i) ) / zp &
393                                                    )
394                   ELSE
395
396!
397!--                   Unstable stratification
398                      h1 = SQRT( SQRT( 1.0 - 16.0 * rifs ) )
399                      h2 = SQRT( SQRT( 1.0 - 16.0 * rifs * z0(j,i) / zp ) )
400
401                      us_wall = kappa * vel_total / (                          &
402                           LOG( zp / z0(j,i) ) -                               &
403                           LOG( ( 1.0 + h1 )**2 * ( 1.0 + h1**2 ) / (          &
404                                ( 1.0 + h2 )**2 * ( 1.0 + h2**2 )   ) ) +      &
405                                2.0 * ( ATAN( h1 ) - ATAN( h2 ) )              &
406                                                    )
407                   ENDIF
408
409!
410!--                Skip step (4) of wall_fluxes, because here rifs is already
411!--                available from (1)
412!
413!--                (5) Compute wall_flux (u'v', v'u', w'v', or w'u')
414
415                   IF ( rifs >= 0.0 )  THEN
416
417!
418!--                   Stable stratification (and neutral)
419                      wall_flux(k,j,i) = kappa *  vel_zp / &
420                          ( LOG( zp/z0(j,i) ) + 5.0*rifs * ( zp-z0(j,i) ) / zp )
421                   ELSE
422
423!
424!--                   Unstable stratification
425                      h1 = SQRT( SQRT( 1.0 - 16.0 * rifs ) )
426                      h2 = SQRT( SQRT( 1.0 - 16.0 * rifs * z0(j,i) / zp ) )
427
428                      wall_flux(k,j,i) = kappa * vel_zp / (                    &
429                           LOG( zp / z0(j,i) ) -                               &
430                           LOG( ( 1.0 + h1 )**2 * ( 1.0 + h1**2 ) / (          &
431                                ( 1.0 + h2 )**2 * ( 1.0 + h2**2 )   ) ) +      &
432                                2.0 * ( ATAN( h1 ) - ATAN( h2 ) )              &
433                                                          )
434                   ENDIF
435                   wall_flux(k,j,i) = - wall_flux(k,j,i) * us_wall
436
437                ENDDO
438
439             ENDIF
440
441          ENDDO
442       ENDDO
443
444    END SUBROUTINE wall_fluxes_e
445
446
447
448!------------------------------------------------------------------------------!
449! Call for grid point i,j
450!------------------------------------------------------------------------------!
451    SUBROUTINE wall_fluxes_e_ij( i, j, nzb_w, nzt_w, wall_flux, a, b, c1, c2 )
452
453       USE arrays_3d
454       USE control_parameters
455       USE grid_variables
456       USE indices
457       USE statistics
458
459       IMPLICIT NONE
460
461       INTEGER ::  i, j, k, kk, nzb_w, nzt_w, wall_index
462       REAL    ::  a, b, c1, c2, h1, h2, u_i, v_i, us_wall, vel_total, vel_zp, &
463                   ws, zp
464
465       REAL ::  rifs
466
467       REAL, DIMENSION(nzb:nzt+1) ::  wall_flux
468
469
470       zp         = 0.5 * ( (a+c1) * dy + (b+c2) * dx )
471       wall_flux  = 0.0
472       wall_index = NINT( a+ 2*b + 3*c1 + 4*c2 )
473
474!
475!--    All subsequent variables are computed for scalar locations.
476       DO  k = nzb_w, nzt_w
477
478!
479!--       (1) Compute rifs, u_i, v_i, and ws
480          IF ( k == nzb_w )  THEN
481             kk = nzb_w
482          ELSE
483             kk = k-1
484          ENDIF
485          rifs  = 0.5 * ( rif_wall(k,j,i,wall_index) +                         &
486                          a * rif_wall(k,j,i+1,1) +  b * rif_wall(k,j+1,i,2) + &
487                          c1 * rif_wall(kk,j,i,3) + c2 * rif_wall(kk,j,i,4)    &
488                        )
489
490          u_i   = 0.5 * ( u(k,j,i) + u(k,j,i+1) )
491          v_i   = 0.5 * ( v(k,j,i) + v(k,j+1,i) )
492          ws    = 0.5 * ( w(k,j,i) + w(k-1,j,i) )
493!
494!--       (2) Compute wall-parallel absolute velocity vel_total and
495!--       interpolate appropriate velocity component vel_zp.
496          vel_total = SQRT( ws**2 + (a+c1) * u_i**2 + (b+c2) * v_i**2 )
497          vel_zp = 0.5 * ( a * u_i + b * v_i + (c1+c2) * ws )
498!
499!--       (3) Compute wall friction velocity us_wall
500          IF ( rifs >= 0.0 )  THEN
501
502!
503!--          Stable stratification (and neutral)
504             us_wall = kappa * vel_total / ( LOG( zp / z0(j,i) ) +             &
505                                            5.0 * rifs * ( zp - z0(j,i) ) / zp &
506                                           )
507          ELSE
508
509!
510!--          Unstable stratification
511             h1 = SQRT( SQRT( 1.0 - 16.0 * rifs ) )
512             h2 = SQRT( SQRT( 1.0 - 16.0 * rifs * z0(j,i) / zp ) )
513
514             us_wall = kappa * vel_total / (                          &
515                  LOG( zp / z0(j,i) ) -                               &
516                  LOG( ( 1.0 + h1 )**2 * ( 1.0 + h1**2 ) / (          &
517                       ( 1.0 + h2 )**2 * ( 1.0 + h2**2 )   ) ) +      &
518                       2.0 * ( ATAN( h1 ) - ATAN( h2 ) )              &
519                                           )
520          ENDIF
521
522!
523!--       Skip step (4) of wall_fluxes, because here rifs is already
524!--       available from (1)
525!
526!--       (5) Compute wall_flux (u'v', v'u', w'v', or w'u')
527!--       First interpolate the velocity (this is different from
528!--       subroutine wall_fluxes because fluxes in subroutine
529!--       wall_fluxes_e are defined at scalar locations).
530          vel_zp = 0.5 * (       a * ( u(k,j,i) + u(k,j,i+1) ) +  &
531                                 b * ( v(k,j,i) + v(k,j+1,i) ) +  &
532                           (c1+c2) * ( w(k,j,i) + w(k-1,j,i) )    &
533                         )
534
535          IF ( rifs >= 0.0 )  THEN
536
537!
538!--          Stable stratification (and neutral)
539             wall_flux(k) = kappa *  vel_zp / &
540                          ( LOG( zp/z0(j,i) ) + 5.0*rifs * ( zp-z0(j,i) ) / zp )
541          ELSE
542
543!
544!--          Unstable stratification
545             h1 = SQRT( SQRT( 1.0 - 16.0 * rifs ) )
546             h2 = SQRT( SQRT( 1.0 - 16.0 * rifs * z0(j,i) / zp ) )
547
548             wall_flux(k) = kappa * vel_zp / (                        &
549                  LOG( zp / z0(j,i) ) -                               &
550                  LOG( ( 1.0 + h1 )**2 * ( 1.0 + h1**2 ) / (          &
551                       ( 1.0 + h2 )**2 * ( 1.0 + h2**2 )   ) ) +      &
552                       2.0 * ( ATAN( h1 ) - ATAN( h2 ) )              &
553                                                 )
554          ENDIF
555          wall_flux(k) = - wall_flux(k) * us_wall
556
557       ENDDO
558
559    END SUBROUTINE wall_fluxes_e_ij
560
561 END MODULE wall_fluxes_mod
Note: See TracBrowser for help on using the repository browser.