Home


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

source: palm/trunk/SOURCE/diffusion_v.f90 @ 91

Last change on this file since 91 was 77, checked in by raasch, 18 years ago

New:
---

particle reflection from vertical walls implemented, particle SGS model adjusted to walls

Wall functions for vertical walls now include diabatic conditions. New subroutines wall_fluxes, wall_fluxes_e. New 4D-array rif_wall.

new d3par-parameter netcdf_64bit_3d to switch on 64bit offset only for 3D files

new d3par-parameter dt_max to define the maximum value for the allowed timestep

new inipar-parameter loop_optimization to control the loop optimization method

new inipar-parameter pt_refrence. If given, this value is used as the reference that in buoyancy terms (otherwise, the instantaneous horizontally averaged temperature is used).

new user interface user_advec_particles

new initializing action "by_user" calls user_init_3d_model and allows the initial setting of all 3d arrays

topography height informations are stored on arrays zu_s_inner and zw_w_inner and output to the 2d/3d NetCDF files

samples added to the user interface which show how to add user-define time series quantities.

calculation/output of precipitation amount, precipitation rate and z0 (by setting "pra*", "prr*", "z0*" with data_output). The time interval on which the precipitation amount is defined is set by new d3par-parameter precipitation_amount_interval

unit 9 opened for debug output (file DEBUG_<pe#>)

Makefile, advec_particles, average_3d_data, buoyancy, calc_precipitation, check_open, check_parameters, data_output_2d, diffusion_e, diffusion_u, diffusion_v, diffusion_w, diffusivities, header, impact_of_latent_heat, init_particles, init_3d_model, modules, netcdf, parin, production_e, read_var_list, read_3d_binary, sum_up_3d_data, user_interface, write_var_list, write_3d_binary

New: wall_fluxes

Changed:

General revision of non-cyclic horizontal boundary conditions: radiation boundary conditions are now used instead of Neumann conditions at the outflow (calculation needs velocity values for t-dt, which are stored on new arrays u_m_l, u_m_r, etc.), calculation of mean outflow is not needed any more, volume flow control is added for the outflow boundary (currently only for the north boundary!!), additional gridpoints along x and y (uxrp, vynp) are not needed any more, routine "boundary_conds" now operates on timelevel t+dt and is not split in two parts (main, uvw_outflow) any more, Neumann boundary conditions at inflow/outflow in case of non-cyclic boundary conditions for all 2d-arrays that are handled by exchange_horiz_2d

The FFT-method for solving the Poisson-equation is now working with Neumann boundary conditions both at the bottom and the top. This requires adjustments of the tridiagonal coefficients and subtracting the horizontally averaged mean from the vertical velocity field.

+age_m in particle_type

Particles-package is now part of the default code ("-p particles" is not needed any more).

Move call of user_actions( 'after_integration' ) below increment of times
and counters. user_actions is now called for each statistic region and has as an argument the number of the respective region (sr)

d3par-parameter data_output_ts removed. Timeseries output for "profil" removed. Timeseries are now switched on by dt_dots. Timeseries data is collected in flow_statistics.

Initial velocities at nzb+1 are regarded for volume flow control in case they have been set zero before (to avoid small timesteps); see new internal parameters u/v_nzb_p1_for_vfc.

q is not allowed to become negative (prognostic_equations).

poisfft_init is only called if fft-solver is switched on (init_pegrid).

d3par-parameter moisture renamed to humidity.

Subversion global revision number is read from mrun and added to the run description header and to the run control (_rc) file.

vtk directives removed from main program.

The uitility routine interpret_config reads PALM environment variables from NAMELIST instead using the system call GETENV.

advec_u_pw, advec_u_up, advec_v_pw, advec_v_up, asselin_filter, check_parameters, coriolis, data_output_dvrp, data_output_ptseries, data_output_ts, data_output_2d, data_output_3d, diffusion_u, diffusion_v, exchange_horiz, exchange_horiz_2d, flow_statistics, header, init_grid, init_particles, init_pegrid, init_rankine, init_pt_anomaly, init_1d_model, init_3d_model, modules, palm, package_parin, parin, poisfft, poismg, prandtl_fluxes, pres, production_e, prognostic_equations, read_var_list, read_3d_binary, sor, swap_timelevel, time_integration, write_var_list, write_3d_binary

Errors:

Bugfix: preset of tendencies te_em, te_um, te_vm in init_1d_model

Bugfix in sample for reading user defined data from restart file (user_init)

Bugfix in setting diffusivities for cases with the outflow damping layer extending over more than one subdomain (init_3d_model)

Check for possible negative humidities in the initial humidity profile.

in Makefile, default suffixes removed from the suffix list to avoid calling of m2c in
# case of .mod files

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