Home


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

source: palm/trunk/SOURCE/diffusion_e.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: 13.2 KB
Line 
1 MODULE diffusion_e_mod
2
3!------------------------------------------------------------------------------!
4! Actual revisions:
5! -----------------
6!
7!
8! Former revisions:
9! -----------------
10! $Id: diffusion_e.f90 77 2007-03-29 04:26:56Z raasch $
11!
12! 65 2007-03-13 12:11:43Z raasch
13! Reference temperature pt_reference can be used in buoyancy term
14!
15! 20 2007-02-26 00:12:32Z raasch
16! Bugfix: ddzw dimensioned 1:nzt"+1"
17! Calculation extended for gridpoint nzt
18!
19! RCS Log replace by Id keyword, revision history cleaned up
20!
21! Revision 1.18  2006/08/04 14:29:43  raasch
22! dissipation is stored in extra array diss if needed later on for calculating
23! the sgs particle velocities
24!
25! Revision 1.1  1997/09/19 07:40:24  raasch
26! Initial revision
27!
28!
29! Description:
30! ------------
31! Diffusion- and dissipation terms for the TKE
32!------------------------------------------------------------------------------!
33
34    PRIVATE
35    PUBLIC diffusion_e
36   
37
38    INTERFACE diffusion_e
39       MODULE PROCEDURE diffusion_e
40       MODULE PROCEDURE diffusion_e_ij
41    END INTERFACE diffusion_e
42 
43 CONTAINS
44
45
46!------------------------------------------------------------------------------!
47! Call for all grid points
48!------------------------------------------------------------------------------!
49    SUBROUTINE diffusion_e( ddzu, dd2zu, ddzw, diss, e, km, l_grid, theta, &
50                            rif, tend, zu )
51
52       USE control_parameters
53       USE grid_variables
54       USE indices
55       USE particle_attributes
56
57       IMPLICIT NONE
58
59       INTEGER ::  i, j, k
60       REAL            ::  dpt_dz, l_stable, phi_m
61       REAL            ::  ddzu(1:nzt+1), dd2zu(1:nzt), ddzw(1:nzt+1), &
62                           l_grid(1:nzt), zu(0:nzt+1)
63       REAL, DIMENSION(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) :: diss, tend
64       REAL, DIMENSION(:,:), POINTER   ::  rif
65       REAL, DIMENSION(:,:,:), POINTER ::  e, km, theta
66       REAL, DIMENSION(nzb+1:nzt,nys:nyn) ::  dissipation, l, ll
67 
68
69!
70!--    This if clause must be outside the k-loop because otherwise
71!--    runtime errors occur with -C hopt on NEC
72       IF ( use_pt_reference )  THEN
73
74          DO  i = nxl, nxr
75             DO  j = nys, nyn
76!
77!--             First, calculate phi-function for eventually adjusting the &
78!--             mixing length to the prandtl mixing length
79                IF ( adjust_mixing_length  .AND.  prandtl_layer )  THEN
80                   IF ( rif(j,i) >= 0.0 )  THEN
81                      phi_m = 1.0 + 5.0 * rif(j,i)
82                   ELSE
83                      phi_m = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rif(j,i) ) )
84                   ENDIF
85                ENDIF
86
87                DO  k = nzb_s_inner(j,i)+1, nzt
88!
89!--                Calculate the mixing length (for dissipation)
90                   dpt_dz = ( theta(k+1,j,i) - theta(k-1,j,i) ) * dd2zu(k)
91                   IF ( dpt_dz > 0.0 ) THEN
92                      l_stable = 0.76 * SQRT( e(k,j,i) ) / &
93                                        SQRT( g / pt_reference * dpt_dz ) + 1E-5
94                   ELSE
95                      l_stable = l_grid(k)
96                   ENDIF
97!
98!--                Adjustment of the mixing length
99                   IF ( wall_adjustment )  THEN
100                      l(k,j)  = MIN( wall_adjustment_factor * zu(k), l_grid(k),&
101                                     l_stable )
102                      ll(k,j) = MIN( wall_adjustment_factor * zu(k), l_grid(k) )
103                   ELSE
104                      l(k,j)  = MIN( l_grid(k), l_stable )
105                      ll(k,j) = l_grid(k)
106                   ENDIF
107                   IF ( adjust_mixing_length  .AND.  prandtl_layer )  THEN
108                      l(k,j)  = MIN( l(k,j),  kappa * zu(k) / phi_m )
109                      ll(k,j) = MIN( ll(k,j), kappa * zu(k) / phi_m )
110                   ENDIF
111
112                ENDDO
113             ENDDO
114
115!
116!--          Calculate the tendency terms
117             DO  j = nys, nyn
118                DO  k = nzb_s_inner(j,i)+1, nzt
119
120                    dissipation(k,j) = ( 0.19 + 0.74 * l(k,j) / ll(k,j) ) * &
121                                       e(k,j,i) * SQRT( e(k,j,i) ) / l(k,j)
122
123                    tend(k,j,i) = tend(k,j,i)                                  &
124                                        + (                                    &
125                          ( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)-e(k,j,i) )  &
126                        - ( km(k,j,i)+km(k,j,i-1) ) * ( e(k,j,i)-e(k,j,i-1) )  &
127                                          ) * ddx2                             &
128                                        + (                                    &
129                          ( km(k,j,i)+km(k,j+1,i) ) * ( e(k,j+1,i)-e(k,j,i) )  &
130                        - ( km(k,j,i)+km(k,j-1,i) ) * ( e(k,j,i)-e(k,j-1,i) )  &
131                                          ) * ddy2                             &
132                                        + (                                    &
133               ( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1) &
134             - ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k)   &
135                                          ) * ddzw(k)                          &
136                             - dissipation(k,j)
137
138                ENDDO
139             ENDDO
140
141!
142!--          Store dissipation if needed for calculating the sgs particle
143!--          velocities
144             IF ( use_sgs_for_particles )  THEN
145                DO  j = nys, nyn
146                   DO  k = nzb_s_inner(j,i)+1, nzt
147                      diss(k,j,i) = dissipation(k,j)
148                   ENDDO
149                ENDDO
150             ENDIF
151
152          ENDDO
153
154       ELSE
155
156          DO  i = nxl, nxr
157             DO  j = nys, nyn
158!
159!--             First, calculate phi-function for eventually adjusting the &
160!--             mixing length to the prandtl mixing length
161                IF ( adjust_mixing_length  .AND.  prandtl_layer )  THEN
162                   IF ( rif(j,i) >= 0.0 )  THEN
163                      phi_m = 1.0 + 5.0 * rif(j,i)
164                   ELSE
165                      phi_m = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rif(j,i) ) )
166                   ENDIF
167                ENDIF
168
169                DO  k = nzb_s_inner(j,i)+1, nzt
170!
171!--                Calculate the mixing length (for dissipation)
172                   dpt_dz = ( theta(k+1,j,i) - theta(k-1,j,i) ) * dd2zu(k)
173                   IF ( dpt_dz > 0.0 ) THEN
174                      l_stable = 0.76 * SQRT( e(k,j,i) ) / &
175                                        SQRT( g / theta(k,j,i) * dpt_dz ) + 1E-5
176                   ELSE
177                      l_stable = l_grid(k)
178                   ENDIF
179!
180!--                Adjustment of the mixing length
181                   IF ( wall_adjustment )  THEN
182                      l(k,j)  = MIN( wall_adjustment_factor * zu(k), l_grid(k),&
183                                     l_stable )
184                      ll(k,j) = MIN( wall_adjustment_factor * zu(k), l_grid(k) )
185                   ELSE
186                      l(k,j)  = MIN( l_grid(k), l_stable )
187                      ll(k,j) = l_grid(k)
188                   ENDIF
189                   IF ( adjust_mixing_length  .AND.  prandtl_layer )  THEN
190                      l(k,j)  = MIN( l(k,j),  kappa * zu(k) / phi_m )
191                      ll(k,j) = MIN( ll(k,j), kappa * zu(k) / phi_m )
192                   ENDIF
193
194                ENDDO
195             ENDDO
196
197!
198!--          Calculate the tendency terms
199             DO  j = nys, nyn
200                DO  k = nzb_s_inner(j,i)+1, nzt
201
202                    dissipation(k,j) = ( 0.19 + 0.74 * l(k,j) / ll(k,j) ) * &
203                                       e(k,j,i) * SQRT( e(k,j,i) ) / l(k,j)
204
205                    tend(k,j,i) = tend(k,j,i)                                  &
206                                        + (                                    &
207                          ( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)-e(k,j,i) )  &
208                        - ( km(k,j,i)+km(k,j,i-1) ) * ( e(k,j,i)-e(k,j,i-1) )  &
209                                          ) * ddx2                             &
210                                        + (                                    &
211                          ( km(k,j,i)+km(k,j+1,i) ) * ( e(k,j+1,i)-e(k,j,i) )  &
212                        - ( km(k,j,i)+km(k,j-1,i) ) * ( e(k,j,i)-e(k,j-1,i) )  &
213                                          ) * ddy2                             &
214                                        + (                                    &
215               ( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1) &
216             - ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k)   &
217                                          ) * ddzw(k)                          &
218                             - dissipation(k,j)
219
220                ENDDO
221             ENDDO
222
223!
224!--          Store dissipation if needed for calculating the sgs particle
225!--          velocities
226             IF ( use_sgs_for_particles )  THEN
227                DO  j = nys, nyn
228                   DO  k = nzb_s_inner(j,i)+1, nzt
229                      diss(k,j,i) = dissipation(k,j)
230                   ENDDO
231                ENDDO
232             ENDIF
233
234          ENDDO
235
236       ENDIF
237
238!
239!--    Boundary condition for dissipation
240       IF ( use_sgs_for_particles )  THEN
241          DO  i = nxl, nxr
242             DO  j = nys, nyn
243                diss(nzb_s_inner(j,i),j,i) = diss(nzb_s_inner(j,i)+1,j,i)
244             ENDDO
245          ENDDO
246       ENDIF
247
248    END SUBROUTINE diffusion_e
249
250
251!------------------------------------------------------------------------------!
252! Call for grid point i,j
253!------------------------------------------------------------------------------!
254    SUBROUTINE diffusion_e_ij( i, j, ddzu, dd2zu, ddzw, diss, e, km, l_grid, &
255                               theta, rif, tend, zu )
256
257       USE control_parameters
258       USE grid_variables
259       USE indices
260       USE particle_attributes
261
262       IMPLICIT NONE
263
264       INTEGER         ::  i, j, k
265       REAL            ::  dpt_dz, l_stable, phi_m
266       REAL            ::  ddzu(1:nzt+1), dd2zu(1:nzt), ddzw(1:nzt+1), &
267                           l_grid(1:nzt), zu(0:nzt+1)
268       REAL, DIMENSION(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) ::  diss, tend
269       REAL, DIMENSION(:,:), POINTER   ::  rif
270       REAL, DIMENSION(:,:,:), POINTER ::  e, km, theta
271       REAL, DIMENSION(nzb+1:nzt)    ::  dissipation, l, ll
272
273
274!
275!--    First, calculate phi-function for eventually adjusting the mixing length
276!--    to the prandtl mixing length
277       IF ( adjust_mixing_length  .AND.  prandtl_layer )  THEN
278          IF ( rif(j,i) >= 0.0 )  THEN
279             phi_m = 1.0 + 5.0 * rif(j,i)
280          ELSE
281             phi_m = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rif(j,i) ) )
282          ENDIF
283       ENDIF
284
285!
286!--    Calculate the mixing length (for dissipation)
287       DO  k = nzb_s_inner(j,i)+1, nzt
288          dpt_dz = ( theta(k+1,j,i) - theta(k-1,j,i) ) * dd2zu(k)
289          IF ( dpt_dz > 0.0 ) THEN
290             IF ( use_pt_reference )  THEN
291                l_stable = 0.76 * SQRT( e(k,j,i) ) / &
292                                  SQRT( g / pt_reference * dpt_dz ) + 1E-5
293             ELSE
294                l_stable = 0.76 * SQRT( e(k,j,i) ) / &
295                                  SQRT( g / theta(k,j,i) * dpt_dz ) + 1E-5
296             ENDIF
297          ELSE
298             l_stable = l_grid(k)
299          ENDIF
300!
301!--       Adjustment of the mixing length
302          IF ( wall_adjustment )  THEN
303             l(k)  = MIN( wall_adjustment_factor * zu(k), l_grid(k), l_stable )
304             ll(k) = MIN( wall_adjustment_factor * zu(k), l_grid(k) )
305          ELSE
306             l(k)  = MIN( l_grid(k), l_stable )
307             ll(k) = l_grid(k)
308          ENDIF
309          IF ( adjust_mixing_length  .AND.  prandtl_layer )  THEN
310             l(k)  = MIN( l(k),  kappa * zu(k) / phi_m )
311             ll(k) = MIN( ll(k), kappa * zu(k) / phi_m )
312          ENDIF
313
314!
315!--       Calculate the tendency term
316          dissipation(k) = ( 0.19 + 0.74 * l(k) / ll(k) ) * e(k,j,i) * &
317                           SQRT( e(k,j,i) ) / l(k)
318
319          tend(k,j,i) = tend(k,j,i)                                           &
320                                       + (                                    &
321                         ( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)-e(k,j,i) )  &
322                       - ( km(k,j,i)+km(k,j,i-1) ) * ( e(k,j,i)-e(k,j,i-1) )  &
323                                         ) * ddx2                             &
324                                       + (                                    &
325                         ( km(k,j,i)+km(k,j+1,i) ) * ( e(k,j+1,i)-e(k,j,i) )  &
326                       - ( km(k,j,i)+km(k,j-1,i) ) * ( e(k,j,i)-e(k,j-1,i) )  &
327                                         ) * ddy2                             &
328                                       + (                                    &
329              ( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1) &
330            - ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k)   &
331                                         ) * ddzw(k)                          &
332                                       - dissipation(k)
333
334       ENDDO
335
336!
337!--    Store dissipation if needed for calculating the sgs particle velocities
338       IF ( use_sgs_for_particles )  THEN
339          DO  k = nzb_s_inner(j,i)+1, nzt
340             diss(k,j,i) = dissipation(k)
341          ENDDO
342!
343!--       Boundary condition for dissipation
344          diss(nzb_s_inner(j,i),j,i) = diss(nzb_s_inner(j,i)+1,j,i)
345       ENDIF
346
347    END SUBROUTINE diffusion_e_ij
348
349 END MODULE diffusion_e_mod
Note: See TracBrowser for help on using the repository browser.
                                                                                                                                                                                                                                                                                                                                                                               
  | Impressum | ©Leibniz Universität Hannover |