source: palm/tags/release-3.2a/SOURCE/prandtl_fluxes.f90 @ 3984

Last change on this file since 3984 was 77, checked in by raasch, 17 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.8 KB
Line 
1 SUBROUTINE prandtl_fluxes
2
3!------------------------------------------------------------------------------!
4! Actual revisions:
5! -----------------
6!
7!
8! Former revisions:
9! -----------------
10! $Id: prandtl_fluxes.f90 77 2007-03-29 04:26:56Z hellstea $
11!
12! 75 2007-03-22 09:54:05Z raasch
13! moisture renamed humidity
14!
15! RCS Log replace by Id keyword, revision history cleaned up
16!
17! Revision 1.19  2006/04/26 12:24:35  raasch
18! +OpenMP directives and optimization (array assignments replaced by DO loops)
19!
20! Revision 1.1  1998/01/23 10:06:06  raasch
21! Initial revision
22!
23!
24! Description:
25! ------------
26! Diagnostic computation of vertical fluxes in the Prandtl layer from the
27! values of the variables at grid point k=1
28!------------------------------------------------------------------------------!
29
30    USE arrays_3d
31    USE control_parameters
32    USE grid_variables
33    USE indices
34
35    IMPLICIT NONE
36
37    INTEGER ::  i, j, k
38    REAL    ::  a, b, rifm, uv_total, z_p
39
40!
41!-- Compute theta*
42    IF ( constant_heatflux )  THEN
43!
44!--    For a given heat flux in the Prandtl layer:
45!--    for u* use the value from the previous time step
46       !$OMP PARALLEL DO
47       DO  i = nxl-1, nxr+1
48          DO  j = nys-1, nyn+1
49             ts(j,i) = -shf(j,i) / ( us(j,i) + 1E-30 )
50!
51!--          ts must be limited, because otherwise overflow may occur in case of
52!--          us=0 when computing rif further below
53             IF ( ts(j,i) < -1.05E5 )  ts = -1.0E5
54             IF ( ts(j,i) >   1.0E5 )  ts =  1.0E5
55          ENDDO
56       ENDDO
57
58    ELSE
59!
60!--    For a given surface temperature:
61!--    (the Richardson number is still the one from the previous time step)
62       !$OMP PARALLEL DO PRIVATE( a, b, k, z_p )
63       DO  i = nxl-1, nxr+1
64          DO  j = nys-1, nyn+1
65
66             k   = nzb_s_inner(j,i)
67             z_p = zu(k+1) - zw(k)
68
69             IF ( rif(j,i) >= 0.0 )  THEN
70!
71!--             Stable stratification
72                ts(j,i) = kappa * ( pt(k+1,j,i) - pt(k,j,i) ) / (          &
73                                  LOG( z_p / z0(j,i) ) +                   &
74                                  5.0 * rif(j,i) * ( z_p - z0(j,i) ) / z_p &
75                                                                )
76             ELSE
77!
78!--             Unstable stratification
79                a = SQRT( 1.0 - 16.0 * rif(j,i) )
80                b = SQRT( 1.0 - 16.0 * rif(j,i) / z_p * z0(j,i) )
81!
82!--             If a borderline case occurs, the formula for stable
83!--             stratification must be used anyway, or else a zero division
84!--             would occur in the argument of the logarithm
85                IF ( a == 1.0  .OR.  b == 1.0 )  THEN
86                   ts(j,i) = kappa * ( pt(k+1,j,i) - pt(k,j,i) ) / (          &
87                                     LOG( z_p / z0(j,i) ) +                   &
88                                     5.0 * rif(j,i) * ( z_p - z0(j,i) ) / z_p &
89                                                                   )
90                ELSE
91                   ts(j,i) = kappa * ( pt(k+1,j,i) - pt(k,j,i) ) / (          &
92                                 LOG( (1.0+b) / (1.0-b) * (1.0-a) / (1.0+a) ) &
93                                                                   )
94                ENDIF
95             ENDIF
96
97          ENDDO
98       ENDDO
99    ENDIF
100
101!
102!-- Compute z_p/L (corresponds to the Richardson-flux number)
103    IF ( .NOT. humidity )  THEN
104       !$OMP PARALLEL DO PRIVATE( k, z_p )
105       DO  i = nxl-1, nxr+1
106          DO  j = nys-1, nyn+1
107             k   = nzb_s_inner(j,i)
108             z_p = zu(k+1) - zw(k)
109             rif(j,i) = z_p * kappa * g * ts(j,i) / &
110                        ( pt(k+1,j,i) * ( us(j,i)**2 + 1E-30 ) )
111!
112!--          Limit the value range of the Richardson numbers.
113!--          This is necessary for very small velocities (u,v --> 0), because
114!--          the absolute value of rif can then become very large, which in
115!--          consequence would result in very large shear stresses and very
116!--          small momentum fluxes (both are generally unrealistic).
117             IF ( rif(j,i) < rif_min )  rif(j,i) = rif_min
118             IF ( rif(j,i) > rif_max )  rif(j,i) = rif_max
119          ENDDO
120       ENDDO
121    ELSE
122       !$OMP PARALLEL DO PRIVATE( k, z_p )
123       DO  i = nxl-1, nxr+1
124          DO  j = nys-1, nyn+1
125             k   = nzb_s_inner(j,i)
126             z_p = zu(k+1) - zw(k)
127             rif(j,i) = z_p * kappa * g *                            &
128                        ( ts(j,i) + 0.61 * pt(k+1,j,i) * qs(j,i) ) / &
129                        ( vpt(k+1,j,i) * ( us(j,i)**2 + 1E-30 ) )
130!
131!--          Limit the value range of the Richardson numbers.
132!--          This is necessary for very small velocities (u,v --> 0), because
133!--          the absolute value of rif can then become very large, which in
134!--          consequence would result in very large shear stresses and very
135!--          small momentum fluxes (both are generally unrealistic).
136             IF ( rif(j,i) < rif_min )  rif(j,i) = rif_min
137             IF ( rif(j,i) > rif_max )  rif(j,i) = rif_max
138          ENDDO
139       ENDDO       
140    ENDIF
141
142!
143!-- Compute u* at the scalars' grid points
144    !$OMP PARALLEL DO PRIVATE( a, b, k, uv_total, z_p )
145    DO  i = nxl, nxr
146       DO  j = nys, nyn
147
148          k   = nzb_s_inner(j,i)
149          z_p = zu(k+1) - zw(k)
150
151!
152!--       Compute the absolute value of the horizontal velocity
153          uv_total = SQRT( ( 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) ) )**2 + &
154                           ( 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) ) )**2   &
155                         )
156
157          IF ( rif(j,i) >= 0.0 )  THEN
158!
159!--          Stable stratification
160             us(j,i) = kappa * uv_total / (                                &
161                                  LOG( z_p / z0(j,i) ) +                   &
162                                  5.0 * rif(j,i) * ( z_p - z0(j,i) ) / z_p &
163                                          )
164          ELSE
165!
166!--          Unstable stratification
167             a = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rif(j,i) ) )
168             b = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rif(j,i) / z_p * z0(j,i) ) )
169!
170!--          If a borderline case occurs, the formula for stable stratification
171!--          must be used anyway, or else a zero division would occur in the
172!--          argument of the logarithm.
173             IF ( a == 1.0  .OR.  b == 1.0 )  THEN
174                us(j,i) = kappa * uv_total / (                                &
175                                     LOG( z_p / z0(j,i) ) +                   &
176                                     5.0 * rif(j,i) * ( z_p - z0(j,i) ) / z_p &
177                                             )
178             ELSE
179                us(j,i) = kappa * uv_total / (                               &
180                              LOG( (1.0+b) / (1.0-b) * (1.0-a) / (1.0+a) ) + &
181                              2.0 * ( ATAN( b ) - ATAN( a ) )                &
182                                             )
183             ENDIF
184          ENDIF
185       ENDDO
186    ENDDO
187
188!
189!-- Compute u'w' for the total model domain.
190!-- First compute the corresponding component of u* and square it.
191    !$OMP PARALLEL DO PRIVATE( a, b, k, rifm, z_p )
192    DO  i = nxl, nxr
193       DO  j = nys, nyn
194
195          k   = nzb_u_inner(j,i)
196          z_p = zu(k+1) - zw(k)
197
198!
199!--       Compute Richardson-flux number for this point
200          rifm = 0.5 * ( rif(j,i-1) + rif(j,i) )
201          IF ( rifm >= 0.0 )  THEN
202!
203!--          Stable stratification
204             usws(j,i) = kappa * u(k+1,j,i) / (                           &
205                                     LOG( z_p / z0(j,i) ) +               &
206                                     5.0 * rifm * ( z_p - z0(j,i) ) / z_p &
207                                              )
208          ELSE
209!
210!--          Unstable stratification
211             a = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifm ) )
212             b = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifm / z_p * z0(j,i) ) )
213!
214!--          If a borderline case occurs, the formula for stable stratification
215!--          must be used anyway, or else a zero division would occur in the
216!--          argument of the logarithm.
217             IF ( a == 1.0  .OR.  B == 1.0 )  THEN
218                usws(j,i) = kappa * u(k+1,j,i) / (                           &
219                                        LOG( z_p / z0(j,i) ) +               &
220                                        5.0 * rifm * ( z_p - z0(j,i) ) / z_p &
221                                                 )
222             ELSE
223                usws(j,i) = kappa * u(k+1,j,i) / (                           &
224                              LOG( (1.0+b) / (1.0-b) * (1.0-a) / (1.0+a) ) + &
225                              2.0 * ( ATAN( b ) - ATAN( a ) )                &
226                                                 )
227             ENDIF
228          ENDIF
229          usws(j,i) = -usws(j,i) * ABS( usws(j,i) )
230       ENDDO
231    ENDDO
232
233!
234!-- Compute v'w' for the total model domain.
235!-- First compute the corresponding component of u* and square it.
236    !$OMP PARALLEL DO PRIVATE( a, b, k, rifm, z_p )
237    DO  i = nxl, nxr
238       DO  j = nys, nyn
239
240          k   = nzb_v_inner(j,i)
241          z_p = zu(k+1) - zw(k)
242
243!
244!--       Compute Richardson-flux number for this point
245          rifm = 0.5 * ( rif(j-1,i) + rif(j,i) )
246          IF ( rifm >= 0.0 )  THEN
247!
248!--          Stable stratification
249             vsws(j,i) = kappa * v(k+1,j,i) / (                           &
250                                     LOG( z_p / z0(j,i) ) +               &
251                                     5.0 * rifm * ( z_p - z0(j,i) ) / z_p &
252                                              )
253          ELSE
254!
255!--          Unstable stratification
256             a = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifm ) )
257             b = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifm / z_p * z0(j,i) ) )
258!
259!--          If a borderline case occurs, the formula for stable stratification
260!--          must be used anyway, or else a zero division would occur in the
261!--          argument of the logarithm.
262             IF ( a == 1.0  .OR.  b == 1.0 )  THEN
263                vsws(j,i) = kappa * v(k+1,j,i) / (                           &
264                                        LOG( z_p / z0(j,i) ) +               &
265                                        5.0 * rifm * ( z_p - z0(j,i) ) / z_p &
266                                                 )
267             ELSE
268                vsws(j,i) = kappa * v(k+1,j,i) / (                           &
269                              LOG( (1.0+b) / (1.0-b) * (1.0-a) / (1.0+a) ) + &
270                              2.0 * ( ATAN( b ) - ATAN( a ) )                &
271                                                 )
272             ENDIF
273          ENDIF
274          vsws(j,i) = -vsws(j,i) * ABS( vsws(j,i) )
275       ENDDO
276    ENDDO
277
278!
279!-- If required compute q*
280    IF ( humidity  .OR.  passive_scalar )  THEN
281       IF ( constant_waterflux )  THEN
282!
283!--       For a given water flux in the Prandtl layer:
284          !$OMP PARALLEL DO
285          DO  i = nxl-1, nxr+1
286             DO  j = nys-1, nyn+1
287                qs(j,i) = -qsws(j,i) / ( us(j,i) + 1E-30 )
288             ENDDO
289          ENDDO
290         
291       ELSE         
292          !$OMP PARALLEL DO PRIVATE( a, b, k, z_p )
293          DO  i = nxl-1, nxr+1
294             DO  j = nys-1, nyn+1
295
296                k   = nzb_s_inner(j,i)
297                z_p = zu(k+1) - zw(k)
298
299                IF ( rif(j,i) >= 0.0 )  THEN
300!
301!--                Stable stratification
302                   qs(j,i) = kappa * ( q(k+1,j,i) - q(k,j,i) ) / (         &
303                                  LOG( z_p / z0(j,i) ) +                   &
304                                  5.0 * rif(j,i) * ( z_p - z0(j,i) ) / z_p &
305                                                                 )
306                ELSE
307!
308!--                Unstable stratification
309                   a = SQRT( 1.0 - 16.0 * rif(j,i) )
310                   b = SQRT( 1.0 - 16.0 * rif(j,i) / z_p * z0(j,i) )
311!
312!--                If a borderline case occurs, the formula for stable
313!--                stratification must be used anyway, or else a zero division
314!--                would occur in the argument of the logarithm.
315                   IF ( a == 1.0  .OR.  b == 1.0 )  THEN
316                      qs(j,i) = kappa * ( q(k+1,j,i) - q(k,j,i) ) / (         &
317                                     LOG( z_p / z0(j,i) ) +                   &
318                                     5.0 * rif(j,i) * ( z_p - z0(j,i) ) / z_p &
319                                                                    )
320                   ELSE
321                      qs(j,i) = kappa * ( q(k+1,j,i) - q(k,j,i) ) / (          &
322                                  LOG( (1.0+b) / (1.0-b) * (1.0-a) / (1.0+a) ) &
323                                                                    )
324                   ENDIF
325                ENDIF
326
327             ENDDO
328          ENDDO
329       ENDIF
330    ENDIF
331
332!
333!-- Exchange the boundaries for u* and the momentum fluxes (fluxes only for
334!-- completeness's sake).
335    CALL exchange_horiz_2d( us )
336    CALL exchange_horiz_2d( usws )
337    CALL exchange_horiz_2d( vsws )
338    IF ( humidity  .OR.  passive_scalar )  CALL exchange_horiz_2d( qsws )
339
340!
341!-- Compute the vertical kinematic heat flux
342    IF ( .NOT. constant_heatflux )  THEN
343       !$OMP PARALLEL DO
344       DO  i = nxl-1, nxr+1
345          DO  j = nys-1, nyn+1
346             shf(j,i) = -ts(j,i) * us(j,i)
347          ENDDO
348       ENDDO
349    ENDIF
350
351!
352!-- Compute the vertical water/scalar flux
353    IF ( .NOT. constant_heatflux .AND. ( humidity .OR. passive_scalar ) ) THEN
354       !$OMP PARALLEL DO
355       DO  i = nxl-1, nxr+1
356          DO  j = nys-1, nyn+1
357             qsws(j,i) = -qs(j,i) * us(j,i)
358          ENDDO
359       ENDDO
360    ENDIF
361
362!
363!-- Bottom boundary condition for the TKE
364    IF ( ibc_e_b == 2 )  THEN
365       !$OMP PARALLEL DO
366       DO  i = nxl-1, nxr+1
367          DO  j = nys-1, nyn+1
368             e(nzb_s_inner(j,i)+1,j,i) = ( us(j,i) / 0.1 )**2
369!
370!--          As a test: cm = 0.4
371!            e(nzb_s_inner(j,i)+1,j,i) = ( us(j,i) / 0.4 )**2
372             e(nzb_s_inner(j,i),j,i)   = e(nzb_s_inner(j,i)+1,j,i)
373          ENDDO
374       ENDDO
375    ENDIF
376
377
378 END SUBROUTINE prandtl_fluxes
Note: See TracBrowser for help on using the repository browser.