Home


About us
PALM
PALM-4U
Installation
Release notes
Documentation
Getting Help
Gallery
Developers area

source: palm/trunk/SOURCE/diffusion_u.f90 @ 399

Last change on this file since 399 was 392, checked in by raasch, 15 years ago

New:
---

Adapted for machine lck
(mrun, mbuild, subjob)

bc_lr/bc_ns in most subroutines replaced by LOGICAL variables bc_lr_cyc,
bc_ns_cyc for speed optimization
(check_parameters, diffusion_u, diffusion_v, diffusion_w, modules)

Additional timestep criterion in case of simulations with plant canopy (timestep)

Check for illegal entries in section_xy|xz|yz that exceed nz+1|ny+1|nx+1
(check_parameters)

Clipping of dvrp output implemented. Default colourtable for particles
implemented, particle attributes (color, dvrp_size) can be set with new
parameters particle_color, particle_dvrpsize, color_interval,
dvrpsize_interval (init_dvrp, data_output_dvrp, modules, user_data_output_dvrp).
Slicer attributes (dvrp) are set with new routine set_slicer_attributes_dvrp
and are controlled with existing parameters slicer_range_limits.
(set_slicer_attributes_dvrp)

Ocean atmosphere coupling allows to use independent precursor runs in order
to account for different spin-up times. The time when coupling has to be
started is given by new inipar parameter coupling_start_time. The precursor
ocean run has to be started using new mrun option "-y" in order to add
appendix "_O" to all output files.
(check_for_restart, check_parameters, data_output_2d, data_output_3d,
data_output_profiles, data_output_ptseries, data_output_spectra,
data_output_tseries, header, init_coupling, modules, mrun,
parin, read_var_list, surface_coupler, time_integration, write_var_list)

Polygon reduction for topography and ground plate isosurface. Reduction level
for buildings can be chosen with parameter cluster_size. (init_dvrp)

External pressure gradient (check_parameters, header, init_3d_model, modules,
parin, prognostic_equations, read_var_list, write_var_list)

New topography case 'single_street_canyon' (header, init_grid, modules, parin,
read_var_list, user_check_parameters, user_header, user_init_grid, write_var_list)

Option to predefine a target bulk velocity for conserve_volume_flow
(check_parameters, header, init_3d_model, modules, parin, read_var_list,
write_var_list)

Option for user defined 2D data output in xy cross sections at z=nzb+1
(data_output_2d, user_data_output_2d)

xy cross section output of surface heatfluxes (latent, sensible)
(average_3d_data, check_parameters, data_output_2d, modules, read_3d_binary,
sum_up_3d_data, write_3d_binary)

average_3d_data, check_for_restart, check_parameters, data_output_2d, data_output_3d, data_output_dvrp, data_output_profiles, data_output_ptseries, data_output_spectra, data_output_tseries, init_coupling, init_dvrp, init_grid, init_3d_model, header, mbuild, modules, mrun, package_parin, parin, prognostic_equations, read_3d_binary, read_var_list, subjob, surface_coupler, timestep, time_integration, user_check_parameters, user_data_output_2d, user_data_output_dvrp, user_header, user_init_grid, write_3d_binary, write_var_list

New: set_particle_attributes, set_slicer_attributes_dvrp

Changed:

lcmuk changed to lc to avoid problems with Intel compiler on sgi-ice
(poisfft)

For extended NetCDF files, the updated title attribute includes an update of
time_average_text where appropriate. (netcdf)

In case of restart runs without extension, initial profiles are not written
to NetCDF-file anymore. (data_output_profiles, modules, read_var_list, write_var_list)

Small change in formatting of the message handling routine concering the output in the
job protocoll. (message)

initializing_actions='read_data_for_recycling' renamed to 'cyclic_fill', now
independent of turbulent_inflow (check_parameters, header, init_3d_model)

2 NetCDF error numbers changed. (data_output_3d)

A Link to the website appendix_a.html is printed for further information
about the possible errors. (message)

Temperature gradient criterion for estimating the boundary layer height
replaced by the gradient criterion of Sullivan et al. (1998). (flow_statistics)

NetCDF unit attribute in timeseries output in case of statistic regions added
(netcdf)

Output of NetCDF messages with aid of message handling routine.
(check_open, close_file, data_output_2d, data_output_3d,
data_output_profiles, data_output_ptseries, data_output_spectra,
data_output_tseries, netcdf, output_particles_netcdf)

Output of messages replaced by message handling routine.
(advec_particles, advec_s_bc, buoyancy, calc_spectra, check_for_restart,
check_open, coriolis, cpu_log, data_output_2d, data_output_3d, data_output_dvrp,
data_output_profiles, data_output_spectra, fft_xy, flow_statistics, header,
init_1d_model, init_3d_model, init_dvrp, init_grid, init_particles, init_pegrid,
netcdf, parin, plant_canopy_model, poisfft_hybrid, poismg, read_3d_binary,
read_var_list, surface_coupler, temperton_fft, timestep, user_actions,
user_data_output_dvrp, user_dvrp_coltab, user_init_grid, user_init_plant_canopy,
user_parin, user_read_restart_data, user_spectra )

Maximum number of tails is calculated from maximum number of particles and
skip_particles_for_tail (init_particles)

Value of vertical_particle_advection may differ for each particle group
(advec_particles, header, modules)

First constant in array den also defined as type double. (eqn_state_seawater)

Parameter dvrp_psize moved from particles_par to dvrp_graphics_par. (package_parin)

topography_grid_convention moved from userpar to inipar (check_parameters,
header, parin, read_var_list, user_check_parameters, user_header,
user_init_grid, user_parin, write_var_list)

Default value of grid_matching changed to strict.

Adjustments for runs on lcxt4 (necessary due to an software update on CRAY) and
for coupled runs on ibmy (mrun, subjob)

advec_particles, advec_s_bc, buoyancy, calc_spectra, check_for_restart, check_open, check_parameters, close_file, coriolis, cpu_log, data_output_2d, data_output_3d, data_output_dvrp, data_output_profiles, data_output_ptseries, data_output_spectra, data_output_tseries, eqn_state_seawater, fft_xy, flow_statistics, header, init_1d_model, init_3d_model, init_dvrp, init_grid, init_particles, init_pegrid, message, mrun, netcdf, output_particles_netcdf, package_parin, parin, plant_canopy_model, poisfft, poisfft_hybrid, poismg, read_3d_binary, read_var_list, sort_particles, subjob, user_check_parameters, user_header, user_init_grid, user_parin, surface_coupler, temperton_fft, timestep, user_actions, user_data_output_dvrp, user_dvrp_coltab, user_init_grid, user_init_plant_canopy, user_parin, user_read_restart_data, user_spectra, write_var_list

Errors:

Bugfix: Initial hydrostatic pressure profile in case of ocean runs is now
calculated in 5 iteration steps. (init_ocean)

Bugfix: wrong sign in buoyancy production of ocean part in case of not using
the reference density (only in 3D routine production_e) (production_e)

Bugfix: output of averaged 2d/3d quantities requires that an avaraging
interval has been set, respective error message is included (check_parameters)

Bugfix: Output on unit 14 only if requested by write_binary.
(user_last_actions)

Bugfix to avoid zero division by km_neutral (production_e)

Bugfix for extended NetCDF files: In order to avoid 'data mode' errors if
updated attributes are larger than their original size, NF90_PUT_ATT is called
in 'define mode' enclosed by NF90_REDEF and NF90_ENDDEF calls. This implies a
possible performance loss; an alternative strategy would be to ensure equal
attribute size in a job chain. (netcdf)

Bugfix: correction of initial volume flow for non-flat topography (init_3d_model)
Bugfix: zero initialization of arrays within buildings for 'cyclic_fill' (init_3d_model)

Bugfix: to_be_resorted => s_av for time-averaged scalars (data_output_2d, data_output_3d)

Bugfix: error in formatting the output (message)

Bugfix: avoid that ngp_2dh_s_inner becomes zero (init_3_model)

Typographical error: unit of wpt in dots_unit (modules)

Bugfix: error in check, if particles moved further than one subdomain length.
This check must not be applied for newly released particles. (advec_particles)

Bugfix: several tail counters are initialized, particle_tail_coordinates is
only written to file if its third index is > 0, arrays for tails are allocated
with a minimum size of 10 tails if there is no tail initially (init_particles,
advec_particles)

Bugfix: pressure included for profile output (check_parameters)

Bugfix: Type of count and count_rate changed to default INTEGER on NEC machines
(cpu_log)

Bugfix: output if particle time series only if particle advection is switched
on. (time_integration)

Bugfix: qsws was calculated in case of constant heatflux = .FALSE. (prandtl_fluxes)

Bugfix: averaging along z is not allowed for 2d quantities (e.g. u* and z0) (data_output_2d)

Typographical errors (netcdf)

If the inversion height calculated by the prerun is zero, inflow_damping_height
must be explicitly specified (init_3d_model)

Small bugfix concerning 3d 64bit netcdf output format (header)

Bugfix: dt_fixed removed from the restart file, because otherwise, no change
from a fixed to a variable timestep would be possible in restart runs.
(read_var_list, write_var_list)

Bugfix: initial setting of time_coupling in coupled restart runs (time_integration)

advec_particles, check_parameters, cpu_log, data_output_2d, data_output_3d, header, init_3d_model, init_particles, init_ocean, modules, netcdf, prandtl_fluxes, production_e, read_var_list, time_integration, user_last_actions, write_var_list
  • Property svn:keywords set to Id
File size: 17.7 KB
Line 
1 MODULE diffusion_u_mod
2
3!------------------------------------------------------------------------------!
4! Actual revisions:
5! -----------------
6!
7!
8! Former revisions:
9! -----------------
10! $Id: diffusion_u.f90 392 2009-09-24 10:39:14Z weinreis $
11!
12! 366 2009-08-25 08:06:27Z raasch
13! bc_ns replaced by bc_ns_cyc
14!
15! 106 2007-08-16 14:30:26Z raasch
16! Momentumflux at top (uswst) included as boundary condition,
17! i loop is starting from nxlu (needed for non-cyclic boundary conditions)
18!
19! 75 2007-03-22 09:54:05Z raasch
20! Wall functions now include diabatic conditions, call of routine wall_fluxes,
21! z0 removed from argument list, uxrp eliminated
22!
23! 20 2007-02-26 00:12:32Z raasch
24! Bugfix: ddzw dimensioned 1:nzt"+1"
25!
26! RCS Log replace by Id keyword, revision history cleaned up
27!
28! Revision 1.15  2006/02/23 10:35:35  raasch
29! nzb_2d replaced by nzb_u_outer in horizontal diffusion and by nzb_u_inner
30! or nzb_diff_u, respectively, in vertical diffusion,
31! wall functions added for north and south walls, +z0 in argument list,
32! terms containing w(k-1,..) are removed from the Prandtl-layer equation
33! because they cause errors at the edges of topography
34! WARNING: loops containing the MAX function are still not properly vectorized!
35!
36! Revision 1.1  1997/09/12 06:23:51  raasch
37! Initial revision
38!
39!
40! Description:
41! ------------
42! Diffusion term of the u-component
43! To do: additional damping (needed for non-cyclic bc) causes bad vectorization
44!        and slows down the speed on NEC about 5-10%
45!------------------------------------------------------------------------------!
46
47    USE wall_fluxes_mod
48
49    PRIVATE
50    PUBLIC diffusion_u
51
52    INTERFACE diffusion_u
53       MODULE PROCEDURE diffusion_u
54       MODULE PROCEDURE diffusion_u_ij
55    END INTERFACE diffusion_u
56
57 CONTAINS
58
59
60!------------------------------------------------------------------------------!
61! Call for all grid points
62!------------------------------------------------------------------------------!
63    SUBROUTINE diffusion_u( ddzu, ddzw, km, km_damp_y, tend, u, usws, uswst, &
64                            v, w )
65
66       USE control_parameters
67       USE grid_variables
68       USE indices
69
70       IMPLICIT NONE
71
72       INTEGER ::  i, j, k
73       REAL    ::  kmym_x, kmym_y, kmyp_x, kmyp_y, kmzm, kmzp
74       REAL    ::  ddzu(1:nzt+1), ddzw(1:nzt+1), km_damp_y(nys-1:nyn+1)
75       REAL    ::  tend(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1)
76       REAL, DIMENSION(:,:),   POINTER ::  usws, uswst
77       REAL, DIMENSION(:,:,:), POINTER ::  km, u, v, w
78       REAL, DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) ::  usvs
79
80!
81!--    First calculate horizontal momentum flux u'v' at vertical walls,
82!--    if neccessary
83       IF ( topography /= 'flat' )  THEN
84          CALL wall_fluxes( usvs, 1.0, 0.0, 0.0, 0.0, nzb_u_inner, &
85                            nzb_u_outer, wall_u )
86       ENDIF
87
88       DO  i = nxlu, nxr
89          DO  j = nys,nyn
90!
91!--          Compute horizontal diffusion
92             DO  k = nzb_u_outer(j,i)+1, nzt
93!
94!--             Interpolate eddy diffusivities on staggered gridpoints
95                kmyp_x = 0.25 * &
96                         ( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) )
97                kmym_x = 0.25 * &
98                         ( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) )
99                kmyp_y = kmyp_x
100                kmym_y = kmym_x
101!
102!--             Increase diffusion at the outflow boundary in case of
103!--             non-cyclic lateral boundaries. Damping is only needed for
104!--             velocity components parallel to the outflow boundary in
105!--             the direction normal to the outflow boundary.
106                IF ( .NOT. bc_ns_cyc )  THEN
107                   kmyp_y = MAX( kmyp_y, km_damp_y(j) )
108                   kmym_y = MAX( kmym_y, km_damp_y(j) )
109                ENDIF
110
111                tend(k,j,i) = tend(k,j,i)                                    &
112                      & + 2.0 * (                                            &
113                      &           km(k,j,i)   * ( u(k,j,i+1) - u(k,j,i)   )  &
114                      &         - km(k,j,i-1) * ( u(k,j,i)   - u(k,j,i-1) )  &
115                      &         ) * ddx2                                     &
116                      & + ( kmyp_y * ( u(k,j+1,i) - u(k,j,i)     ) * ddy     &
117                      &   + kmyp_x * ( v(k,j+1,i) - v(k,j+1,i-1) ) * ddx     &
118                      &   - kmym_y * ( u(k,j,i) - u(k,j-1,i) ) * ddy         &
119                      &   - kmym_x * ( v(k,j,i) - v(k,j,i-1) ) * ddx         &
120                      &   ) * ddy
121             ENDDO
122
123!
124!--          Wall functions at the north and south walls, respectively
125             IF ( wall_u(j,i) /= 0.0 )  THEN
126
127                DO  k = nzb_u_inner(j,i)+1, nzb_u_outer(j,i)
128                   kmyp_x = 0.25 * &
129                            ( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) )
130                   kmym_x = 0.25 * &
131                            ( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) )
132                   kmyp_y = kmyp_x
133                   kmym_y = kmym_x
134!
135!--                Increase diffusion at the outflow boundary in case of
136!--                non-cyclic lateral boundaries. Damping is only needed for
137!--                velocity components parallel to the outflow boundary in
138!--                the direction normal to the outflow boundary.
139                   IF ( .NOT. bc_ns_cyc )  THEN
140                      kmyp_y = MAX( kmyp_y, km_damp_y(j) )
141                      kmym_y = MAX( kmym_y, km_damp_y(j) )
142                   ENDIF
143
144                   tend(k,j,i) = tend(k,j,i)                                   &
145                                 + 2.0 * (                                     &
146                                       km(k,j,i)   * ( u(k,j,i+1) - u(k,j,i) ) &
147                                     - km(k,j,i-1) * ( u(k,j,i) - u(k,j,i-1) ) &
148                                         ) * ddx2                              &
149                                 + (   fyp(j,i) * (                            &
150                                  kmyp_y * ( u(k,j+1,i) - u(k,j,i)     ) * ddy &
151                                + kmyp_x * ( v(k,j+1,i) - v(k,j+1,i-1) ) * ddx &
152                                                  )                            &
153                                     - fym(j,i) * (                            &
154                                  kmym_y * ( u(k,j,i) - u(k,j-1,i) ) * ddy     &
155                                + kmym_x * ( v(k,j,i) - v(k,j,i-1) ) * ddx     &
156                                                  )                            &
157                                     + wall_u(j,i) * usvs(k,j,i)               &
158                                   ) * ddy
159                ENDDO
160             ENDIF
161
162!
163!--          Compute vertical diffusion. In case of simulating a Prandtl layer,
164!--          index k starts at nzb_u_inner+2.
165             DO  k = nzb_diff_u(j,i), nzt_diff
166!
167!--             Interpolate eddy diffusivities on staggered gridpoints
168                kmzp = 0.25 * &
169                       ( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) )
170                kmzm = 0.25 * &
171                       ( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
172
173                tend(k,j,i) = tend(k,j,i)                                    &
174                      & + ( kmzp * ( ( u(k+1,j,i) - u(k,j,i)   ) * ddzu(k+1) &
175                      &            + ( w(k,j,i)   - w(k,j,i-1) ) * ddx       &
176                      &            )                                         &
177                      &   - kmzm * ( ( u(k,j,i)   - u(k-1,j,i)   ) * ddzu(k) &
178                      &            + ( w(k-1,j,i) - w(k-1,j,i-1) ) * ddx     &
179                      &            )                                         &
180                      &   ) * ddzw(k)
181             ENDDO
182
183!
184!--          Vertical diffusion at the first grid point above the surface,
185!--          if the momentum flux at the bottom is given by the Prandtl law or
186!--          if it is prescribed by the user.
187!--          Difference quotient of the momentum flux is not formed over half
188!--          of the grid spacing (2.0*ddzw(k)) any more, since the comparison
189!--          with other (LES) modell showed that the values of the momentum
190!--          flux becomes too large in this case.
191!--          The term containing w(k-1,..) (see above equation) is removed here
192!--          because the vertical velocity is assumed to be zero at the surface.
193             IF ( use_surface_fluxes )  THEN
194                k = nzb_u_inner(j,i)+1
195!
196!--             Interpolate eddy diffusivities on staggered gridpoints
197                kmzp = 0.25 * &
198                      ( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) )
199                kmzm = 0.25 * &
200                      ( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
201
202                tend(k,j,i) = tend(k,j,i)                                    &
203                      & + ( kmzp * ( w(k,j,i)   - w(k,j,i-1)   ) * ddx       &
204                      &   ) * ddzw(k)                                        &
205                      & + ( kmzp * ( u(k+1,j,i) - u(k,j,i)     ) * ddzu(k+1) &
206                      &   + usws(j,i)                                        &
207                      &   ) * ddzw(k)
208             ENDIF
209
210!
211!--          Vertical diffusion at the first gridpoint below the top boundary,
212!--          if the momentum flux at the top is prescribed by the user
213             IF ( use_top_fluxes  .AND.  constant_top_momentumflux )  THEN
214                k = nzt
215!
216!--             Interpolate eddy diffusivities on staggered gridpoints
217                kmzp = 0.25 * &
218                       ( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) )
219                kmzm = 0.25 * &
220                       ( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
221
222                tend(k,j,i) = tend(k,j,i)                                    &
223                      & - ( kmzm * ( w(k-1,j,i) - w(k-1,j,i-1) ) * ddx       &
224                      &   ) * ddzw(k)                                        &
225                      & + ( -uswst(j,i)                                      &
226                      &   - kmzm * ( u(k,j,i)   - u(k-1,j,i)   ) * ddzu(k)   &
227                      &   ) * ddzw(k)
228             ENDIF
229
230          ENDDO
231       ENDDO
232
233    END SUBROUTINE diffusion_u
234
235
236!------------------------------------------------------------------------------!
237! Call for grid point i,j
238!------------------------------------------------------------------------------!
239    SUBROUTINE diffusion_u_ij( i, j, ddzu, ddzw, km, km_damp_y, tend, u, usws, &
240                               uswst, v, w )
241
242       USE control_parameters
243       USE grid_variables
244       USE indices
245
246       IMPLICIT NONE
247
248       INTEGER ::  i, j, k
249       REAL    ::  kmym_x, kmym_y, kmyp_x, kmyp_y, kmzm, kmzp
250       REAL    ::  ddzu(1:nzt+1), ddzw(1:nzt+1), km_damp_y(nys-1:nyn+1)
251       REAL    ::  tend(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1)
252       REAL, DIMENSION(nzb:nzt+1)      ::  usvs
253       REAL, DIMENSION(:,:),   POINTER ::  usws, uswst
254       REAL, DIMENSION(:,:,:), POINTER ::  km, u, v, w
255
256!
257!--    Compute horizontal diffusion
258       DO  k = nzb_u_outer(j,i)+1, nzt
259!
260!--       Interpolate eddy diffusivities on staggered gridpoints
261          kmyp_x = 0.25 * ( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) )
262          kmym_x = 0.25 * ( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) )
263          kmyp_y = kmyp_x
264          kmym_y = kmym_x
265
266!
267!--       Increase diffusion at the outflow boundary in case of non-cyclic
268!--       lateral boundaries. Damping is only needed for velocity components
269!--       parallel to the outflow boundary in the direction normal to the
270!--       outflow boundary.
271          IF ( .NOT. bc_ns_cyc )  THEN
272             kmyp_y = MAX( kmyp_y, km_damp_y(j) )
273             kmym_y = MAX( kmym_y, km_damp_y(j) )
274          ENDIF
275
276          tend(k,j,i) = tend(k,j,i)                                          &
277                      & + 2.0 * (                                            &
278                      &           km(k,j,i)   * ( u(k,j,i+1) - u(k,j,i)   )  &
279                      &         - km(k,j,i-1) * ( u(k,j,i)   - u(k,j,i-1) )  &
280                      &         ) * ddx2                                     &
281                      & + ( kmyp_y * ( u(k,j+1,i) - u(k,j,i)     ) * ddy     &
282                      &   + kmyp_x * ( v(k,j+1,i) - v(k,j+1,i-1) ) * ddx     &
283                      &   - kmym_y * ( u(k,j,i) - u(k,j-1,i) ) * ddy         &
284                      &   - kmym_x * ( v(k,j,i) - v(k,j,i-1) ) * ddx         &
285                      &   ) * ddy
286       ENDDO
287
288!
289!--    Wall functions at the north and south walls, respectively
290       IF ( wall_u(j,i) .NE. 0.0 )  THEN
291
292!
293!--       Calculate the horizontal momentum flux u'v'
294          CALL wall_fluxes( i, j, nzb_u_inner(j,i)+1, nzb_u_outer(j,i),  &
295                            usvs, 1.0, 0.0, 0.0, 0.0 )
296
297          DO  k = nzb_u_inner(j,i)+1, nzb_u_outer(j,i)
298             kmyp_x = 0.25 * ( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) )
299             kmym_x = 0.25 * ( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) )
300             kmyp_y = kmyp_x
301             kmym_y = kmym_x
302!
303!--          Increase diffusion at the outflow boundary in case of
304!--          non-cyclic lateral boundaries. Damping is only needed for
305!--          velocity components parallel to the outflow boundary in
306!--          the direction normal to the outflow boundary.
307             IF ( .NOT. bc_ns_cyc )  THEN
308                kmyp_y = MAX( kmyp_y, km_damp_y(j) )
309                kmym_y = MAX( kmym_y, km_damp_y(j) )
310             ENDIF
311
312             tend(k,j,i) = tend(k,j,i)                                         &
313                                 + 2.0 * (                                     &
314                                       km(k,j,i)   * ( u(k,j,i+1) - u(k,j,i) ) &
315                                     - km(k,j,i-1) * ( u(k,j,i) - u(k,j,i-1) ) &
316                                         ) * ddx2                              &
317                                 + (   fyp(j,i) * (                            &
318                                  kmyp_y * ( u(k,j+1,i) - u(k,j,i)     ) * ddy &
319                                + kmyp_x * ( v(k,j+1,i) - v(k,j+1,i-1) ) * ddx &
320                                                  )                            &
321                                     - fym(j,i) * (                            &
322                                  kmym_y * ( u(k,j,i) - u(k,j-1,i) ) * ddy     &
323                                + kmym_x * ( v(k,j,i) - v(k,j,i-1) ) * ddx     &
324                                                  )                            &
325                                     + wall_u(j,i) * usvs(k)                   &
326                                   ) * ddy
327          ENDDO
328       ENDIF
329
330!
331!--    Compute vertical diffusion. In case of simulating a Prandtl layer,
332!--    index k starts at nzb_u_inner+2.
333       DO  k = nzb_diff_u(j,i), nzt_diff
334!
335!--       Interpolate eddy diffusivities on staggered gridpoints
336          kmzp = 0.25 * ( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) )
337          kmzm = 0.25 * ( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
338
339          tend(k,j,i) = tend(k,j,i)                                          &
340                      & + ( kmzp * ( ( u(k+1,j,i) - u(k,j,i)   ) * ddzu(k+1) &
341                      &            + ( w(k,j,i)   - w(k,j,i-1) ) * ddx       &
342                      &            )                                         &
343                      &   - kmzm * ( ( u(k,j,i)   - u(k-1,j,i)   ) * ddzu(k) &
344                      &            + ( w(k-1,j,i) - w(k-1,j,i-1) ) * ddx     &
345                      &            )                                         &
346                      &   ) * ddzw(k)
347       ENDDO
348
349!
350!--    Vertical diffusion at the first grid point above the surface, if the
351!--    momentum flux at the bottom is given by the Prandtl law or if it is
352!--    prescribed by the user.
353!--    Difference quotient of the momentum flux is not formed over half of
354!--    the grid spacing (2.0*ddzw(k)) any more, since the comparison with
355!--    other (LES) modell showed that the values of the momentum flux becomes
356!--    too large in this case.
357!--    The term containing w(k-1,..) (see above equation) is removed here
358!--    because the vertical velocity is assumed to be zero at the surface.
359       IF ( use_surface_fluxes )  THEN
360          k = nzb_u_inner(j,i)+1
361!
362!--       Interpolate eddy diffusivities on staggered gridpoints
363          kmzp = 0.25 * ( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) )
364          kmzm = 0.25 * ( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
365
366          tend(k,j,i) = tend(k,j,i)                                          &
367                      & + ( kmzp * ( w(k,j,i)   - w(k,j,i-1)   ) * ddx       &
368                      &   ) * ddzw(k)                                        &
369                      & + ( kmzp * ( u(k+1,j,i) - u(k,j,i)     ) * ddzu(k+1) &
370                      &   + usws(j,i)                                        &
371                      &   ) * ddzw(k)
372       ENDIF
373
374!
375!--    Vertical diffusion at the first gridpoint below the top boundary,
376!--    if the momentum flux at the top is prescribed by the user
377       IF ( use_top_fluxes  .AND.  constant_top_momentumflux )  THEN
378          k = nzt
379!
380!--       Interpolate eddy diffusivities on staggered gridpoints
381          kmzp = 0.25 * &
382                 ( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) )
383          kmzm = 0.25 * &
384                 ( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
385
386          tend(k,j,i) = tend(k,j,i)                                          &
387                      & - ( kmzm * ( w(k-1,j,i) - w(k-1,j,i-1) ) * ddx       &
388                      &   ) * ddzw(k)                                        &
389                      & + ( -uswst(j,i)                                      &
390                      &   - kmzm * ( u(k,j,i)   - u(k-1,j,i)   ) * ddzu(k)   &
391                      &   ) * ddzw(k)
392       ENDIF
393
394    END SUBROUTINE diffusion_u_ij
395
396 END MODULE diffusion_u_mod
Note: See TracBrowser for help on using the repository browser.
                                                                                                                                                                                                                                                                                                                                                                               
  | Impressum | ©Leibniz Universität Hannover |