Home


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

source: palm/trunk/SOURCE/diffusion_w.f90 @ 357

Last change on this file since 357 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.3 KB
RevLine 
[1]1 MODULE diffusion_w_mod
2
3!------------------------------------------------------------------------------!
4! Actual revisions:
5! -----------------
[77]6!
[1]7!
8! Former revisions:
9! -----------------
[3]10! $Id: diffusion_w.f90 77 2007-03-29 04:26:56Z heinze $
[39]11!
[77]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
15!
[39]16! 20 2007-02-26 00:12:32Z raasch
17! Bugfix: ddzw dimensioned 1:nzt"+1"
18!
[3]19! RCS Log replace by Id keyword, revision history cleaned up
20!
[1]21! Revision 1.12  2006/02/23 10:38:03  raasch
22! nzb_2d replaced by nzb_w_outer, wall functions added for all vertical walls,
23! +z0 in argument list
24! WARNING: loops containing the MAX function are still not properly vectorized!
25!
26! Revision 1.1  1997/09/12 06:24:11  raasch
27! Initial revision
28!
29!
30! Description:
31! ------------
32! Diffusion term of the w-component
33!------------------------------------------------------------------------------!
34
[56]35    USE wall_fluxes_mod
36
[1]37    PRIVATE
38    PUBLIC diffusion_w
39
40    INTERFACE diffusion_w
41       MODULE PROCEDURE diffusion_w
42       MODULE PROCEDURE diffusion_w_ij
43    END INTERFACE diffusion_w
44
45 CONTAINS
46
47
48!------------------------------------------------------------------------------!
49! Call for all grid points
50!------------------------------------------------------------------------------!
51    SUBROUTINE diffusion_w( ddzu, ddzw, km, km_damp_x, km_damp_y, tend, u, v, &
[57]52                            w )
[1]53
54       USE control_parameters
55       USE grid_variables
56       USE indices
57
58       IMPLICIT NONE
59
60       INTEGER ::  i, j, k
61       REAL    ::  kmxm_x, kmxm_z, kmxp_x, kmxp_z, kmym_y, kmym_z, kmyp_y, &
[51]62                   kmyp_z
[20]63       REAL    ::  ddzu(1:nzt+1), ddzw(1:nzt+1), km_damp_x(nxl-1:nxr+1), &
[1]64                   km_damp_y(nys-1:nyn+1)
65       REAL    ::  tend(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1)
66       REAL, DIMENSION(:,:,:), POINTER ::  km, u, v, w
[56]67       REAL, DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) ::  wsus, wsvs
[1]68
69
[56]70!
71!--    First calculate horizontal momentum flux w'u' and/or w'v' at vertical
72!--    walls, if neccessary
73       IF ( topography /= 'flat' )  THEN
[75]74          CALL wall_fluxes( wsus, 0.0, 0.0, 0.0, 1.0, nzb_w_inner, &
[56]75                            nzb_w_outer, wall_w_x )
[75]76          CALL wall_fluxes( wsvs, 0.0, 0.0, 1.0, 0.0, nzb_w_inner, &
[56]77                            nzb_w_outer, wall_w_y )
78       ENDIF
79
[1]80       DO  i = nxl, nxr
81          DO  j = nys, nyn
82             DO  k = nzb_w_outer(j,i)+1, nzt-1
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+1,j,i)+km(k+1,j,i+1) )
87                kmxm_x = 0.25 * &
88                         ( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) )
89                kmxp_z = kmxp_x
90                kmxm_z = kmxm_x
91                kmyp_y = 0.25 * &
92                         ( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) )
93                kmym_y = 0.25 * &
94                         ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
95                kmyp_z = kmyp_y
96                kmym_z = kmym_y
97!
98!--             Increase diffusion at the outflow boundary in case of
99!--             non-cyclic lateral boundaries. Damping is only needed for
100!--             velocity components parallel to the outflow boundary in
101!--             the direction normal to the outflow boundary.
102                IF ( bc_lr /= 'cyclic' )  THEN
103                   kmxp_x = MAX( kmxp_x, km_damp_x(i) )
104                   kmxm_x = MAX( kmxm_x, km_damp_x(i) )
105                ENDIF
106                IF ( bc_ns /= 'cyclic' )  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                      & + ( kmxp_x * ( w(k,j,i+1)   - w(k,j,i)   ) * ddx       &
113                      &   + kmxp_z * ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1) &
114                      &   - kmxm_x * ( w(k,j,i)   - w(k,j,i-1) ) * ddx         &
115                      &   - kmxm_z * ( u(k+1,j,i) - u(k,j,i)   ) * ddzu(k+1)   &
116                      &   ) * ddx                                              &
117                      & + ( kmyp_y * ( w(k,j+1,i)   - w(k,j,i)   ) * ddy       &
118                      &   + kmyp_z * ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1) &
119                      &   - kmym_y * ( w(k,j,i)   - w(k,j-1,i) ) * ddy         &
120                      &   - kmym_z * ( v(k+1,j,i) - v(k,j,i)   ) * ddzu(k+1)   &
121                      &   ) * ddy                                              &
122                      & + 2.0 * (                                              &
123                      &   km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1)  &
124                      & - km(k,j,i)   * ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k)  &
125                      &         ) * ddzu(k+1)
126             ENDDO
127
128!
129!--          Wall functions at all vertical walls, where necessary
130             IF ( wall_w_x(j,i) /= 0.0  .OR.  wall_w_y(j,i) /= 0.0 )  THEN
[51]131
[1]132                DO  k = nzb_w_inner(j,i)+1, nzb_w_outer(j,i)
133!
134!--                Interpolate eddy diffusivities on staggered gridpoints
135                   kmxp_x = 0.25 * &
136                            ( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) )
137                   kmxm_x = 0.25 * &
138                            ( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) )
139                   kmxp_z = kmxp_x
140                   kmxm_z = kmxm_x
141                   kmyp_y = 0.25 * &
142                            ( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) )
143                   kmym_y = 0.25 * &
144                            ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
145                   kmyp_z = kmyp_y
146                   kmym_z = kmym_y
147!
148!--                Increase diffusion at the outflow boundary in case of
149!--                non-cyclic lateral boundaries. Damping is only needed for
150!--                velocity components parallel to the outflow boundary in
151!--                the direction normal to the outflow boundary.
152                   IF ( bc_lr /= 'cyclic' )  THEN
153                      kmxp_x = MAX( kmxp_x, km_damp_x(i) )
154                      kmxm_x = MAX( kmxm_x, km_damp_x(i) )
155                   ENDIF
156                   IF ( bc_ns /= 'cyclic' )  THEN
157                      kmyp_y = MAX( kmyp_y, km_damp_y(j) )
158                      kmym_y = MAX( kmym_y, km_damp_y(j) )
159                   ENDIF
160
161                   tend(k,j,i) = tend(k,j,i)                                   &
162                                 + (   fwxp(j,i) * (                           &
163                            kmxp_x * ( w(k,j,i+1)   - w(k,j,i)   ) * ddx       &
164                          + kmxp_z * ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1) &
165                                                   )                           &
166                                     - fwxm(j,i) * (                           &
167                            kmxm_x * ( w(k,j,i)     - w(k,j,i-1) ) * ddx       &
168                          + kmxm_z * ( u(k+1,j,i)   - u(k,j,i)   ) * ddzu(k+1) &
169                                                   )                           &
[56]170                                     + wall_w_x(j,i) * wsus(k,j,i)             &
[1]171                                   ) * ddx                                     &
172                                 + (   fwyp(j,i) * (                           &
173                            kmyp_y * ( w(k,j+1,i)   - w(k,j,i)   ) * ddy       &
174                          + kmyp_z * ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1) &
175                                                   )                           &
176                                     - fwym(j,i) * (                           &
177                            kmym_y * ( w(k,j,i)     - w(k,j-1,i) ) * ddy       &
178                          + kmym_z * ( v(k+1,j,i)   - v(k,j,i)   ) * ddzu(k+1) &
179                                                   )                           &
[56]180                                     + wall_w_y(j,i) * wsvs(k,j,i)             &
[1]181                                   ) * ddy                                     &
182                                 + 2.0 * (                                     &
183                           km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1) &
184                         - km(k,j,i)   * ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k) &
185                                         ) * ddzu(k+1)
186                ENDDO
187             ENDIF
188
189          ENDDO
190       ENDDO
191
192    END SUBROUTINE diffusion_w
193
194
195!------------------------------------------------------------------------------!
196! Call for grid point i,j
197!------------------------------------------------------------------------------!
198    SUBROUTINE diffusion_w_ij( i, j, ddzu, ddzw, km, km_damp_x, km_damp_y, &
[57]199                               tend, u, v, w )
[1]200
201       USE control_parameters
202       USE grid_variables
203       USE indices
204
205       IMPLICIT NONE
206
207       INTEGER ::  i, j, k
208       REAL    ::  kmxm_x, kmxm_z, kmxp_x, kmxp_z, kmym_y, kmym_z, kmyp_y, &
[51]209                   kmyp_z
[20]210       REAL    ::  ddzu(1:nzt+1), ddzw(1:nzt+1), km_damp_x(nxl-1:nxr+1), &
[1]211                   km_damp_y(nys-1:nyn+1)
212       REAL    ::  tend(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1)
[51]213       REAL, DIMENSION(nzb:nzt+1)      ::  wsus, wsvs
[1]214       REAL, DIMENSION(:,:,:), POINTER ::  km, u, v, w
215
216
217       DO  k = nzb_w_outer(j,i)+1, nzt-1
218!
219!--       Interpolate eddy diffusivities on staggered gridpoints
220          kmxp_x = 0.25 * ( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) )
221          kmxm_x = 0.25 * ( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) )
222          kmxp_z = kmxp_x
223          kmxm_z = kmxm_x
224          kmyp_y = 0.25 * ( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) )
225          kmym_y = 0.25 * ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
226          kmyp_z = kmyp_y
227          kmym_z = kmym_y
228!
229!--       Increase diffusion at the outflow boundary in case of non-cyclic
230!--       lateral boundaries. Damping is only needed for velocity components
231!--       parallel to the outflow boundary in the direction normal to the
232!--       outflow boundary.
233          IF ( bc_lr /= 'cyclic' )  THEN
234             kmxp_x = MAX( kmxp_x, km_damp_x(i) )
235             kmxm_x = MAX( kmxm_x, km_damp_x(i) )
236          ENDIF
237          IF ( bc_ns /= 'cyclic' )  THEN
238             kmyp_y = MAX( kmyp_y, km_damp_y(j) )
239             kmym_y = MAX( kmym_y, km_damp_y(j) )
240          ENDIF
241
242          tend(k,j,i) = tend(k,j,i)                                            &
243                      & + ( kmxp_x * ( w(k,j,i+1)   - w(k,j,i)   ) * ddx       &
244                      &   + kmxp_z * ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1) &
245                      &   - kmxm_x * ( w(k,j,i)   - w(k,j,i-1) ) * ddx         &
246                      &   - kmxm_z * ( u(k+1,j,i) - u(k,j,i)   ) * ddzu(k+1)   &
247                      &   ) * ddx                                              &
248                      & + ( kmyp_y * ( w(k,j+1,i)   - w(k,j,i)   ) * ddy       &
249                      &   + kmyp_z * ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1) &
250                      &   - kmym_y * ( w(k,j,i)   - w(k,j-1,i) ) * ddy         &
251                      &   - kmym_z * ( v(k+1,j,i) - v(k,j,i)   ) * ddzu(k+1)   &
252                      &   ) * ddy                                              &
253                      & + 2.0 * (                                              &
254                      &   km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1)  &
255                      & - km(k,j,i)   * ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k)  &
256                      &         ) * ddzu(k+1)
257       ENDDO
258
259!
260!--    Wall functions at all vertical walls, where necessary
261       IF ( wall_w_x(j,i) /= 0.0  .OR.  wall_w_y(j,i) /= 0.0 )  THEN
[51]262
263!
264!--       Calculate the horizontal momentum fluxes w'u' and/or w'v'
265          IF ( wall_w_x(j,i) /= 0.0 )  THEN
266             CALL wall_fluxes( i, j, nzb_w_inner(j,i)+1, nzb_w_outer(j,i), &
267                               wsus, 0.0, 0.0, 0.0, 1.0 )
268          ELSE
269             wsus = 0.0
270          ENDIF
271
272          IF ( wall_w_y(j,i) /= 0.0 )  THEN
273             CALL wall_fluxes( i, j, nzb_w_inner(j,i)+1, nzb_w_outer(j,i),  &
274                               wsvs, 0.0, 0.0, 1.0, 0.0 )
275          ELSE
276             wsvs = 0.0
277          ENDIF
278
[1]279          DO  k = nzb_w_inner(j,i)+1, nzb_w_outer(j,i)
280!
281!--          Interpolate eddy diffusivities on staggered gridpoints
282             kmxp_x = 0.25 * ( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) )
283             kmxm_x = 0.25 * ( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) )
284             kmxp_z = kmxp_x
285             kmxm_z = kmxm_x
286             kmyp_y = 0.25 * ( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) )
287             kmym_y = 0.25 * ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
288             kmyp_z = kmyp_y
289             kmym_z = kmym_y
290!
291!--          Increase diffusion at the outflow boundary in case of
292!--          non-cyclic lateral boundaries. Damping is only needed for
293!--          velocity components parallel to the outflow boundary in
294!--          the direction normal to the outflow boundary.
295             IF ( bc_lr /= 'cyclic' )  THEN
296                kmxp_x = MAX( kmxp_x, km_damp_x(i) )
297                kmxm_x = MAX( kmxm_x, km_damp_x(i) )
298             ENDIF
299             IF ( bc_ns /= 'cyclic' )  THEN
300                kmyp_y = MAX( kmyp_y, km_damp_y(j) )
301                kmym_y = MAX( kmym_y, km_damp_y(j) )
302             ENDIF
303
304             tend(k,j,i) = tend(k,j,i)                                         &
305                                 + (   fwxp(j,i) * (                           &
306                            kmxp_x * ( w(k,j,i+1)   - w(k,j,i)   ) * ddx       &
307                          + kmxp_z * ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1) &
308                                                   )                           &
309                                     - fwxm(j,i) * (                           &
310                            kmxm_x * ( w(k,j,i)     - w(k,j,i-1) ) * ddx       &
311                          + kmxm_z * ( u(k+1,j,i)   - u(k,j,i)   ) * ddzu(k+1) &
312                                                   )                           &
[51]313                                     + wall_w_x(j,i) * wsus(k)                 &
[1]314                                   ) * ddx                                     &
315                                 + (   fwyp(j,i) * (                           &
316                            kmyp_y * ( w(k,j+1,i)   - w(k,j,i)   ) * ddy       &
317                          + kmyp_z * ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1) &
318                                                   )                           &
319                                     - fwym(j,i) * (                           &
320                            kmym_y * ( w(k,j,i)     - w(k,j-1,i) ) * ddy       &
321                          + kmym_z * ( v(k+1,j,i)   - v(k,j,i)   ) * ddzu(k+1) &
322                                                   )                           &
[51]323                                     + wall_w_y(j,i) * wsvs(k)                 &
[1]324                                   ) * ddy                                     &
325                                 + 2.0 * (                                     &
326                           km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1) &
327                         - km(k,j,i)   * ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k) &
328                                         ) * ddzu(k+1)
329          ENDDO
330       ENDIF
331
332    END SUBROUTINE diffusion_w_ij
333
334 END MODULE diffusion_w_mod
Note: See TracBrowser for help on using the repository browser.
                                                                                                                                                                                                                                                                                                                                                                               
  | Impressum | ©Leibniz Universität Hannover |