source: palm/trunk/SOURCE/buoyancy.f90 @ 220

Last change on this file since 220 was 139, checked in by raasch, 16 years ago

New:
---

Plant canopy model of Watanabe (2004,BLM 112,307-341) added.
It can be switched on by the inipar parameter plant_canopy.
The inipar parameter canopy_mode can be used to prescribe a
plant canopy type. The default case is a homogeneous plant
canopy. Heterogeneous distributions of the leaf area
density and the canopy drag coefficient can be defined in the
new routine user_init_plant_canopy (user_interface).
The inipar parameters lad_surface, lad_vertical_gradient and
lad_vertical_gradient_level can be used in order to
prescribe the vertical profile of leaf area density. The
inipar parameter drag_coefficient determines the canopy
drag coefficient.
Finally, the inipar parameter pch_index determines the
index of the upper boundary of the plant canopy.

Allow new case bc_uv_t = 'dirichlet_0' for channel flow.

For unknown variables (CASE DEFAULT) call new subroutine user_data_output_dvrp

Pressure boundary conditions for vertical walls added to the multigrid solver.
They are applied using new wall flag arrays (wall_flags_..) which are defined
for each grid level. New argument gls added to routine user_init_grid
(user_interface).

Frequence of sorting particles can be controlled with new particles_par
parameter dt_sort_particles. Sorting is moved from the SGS timestep loop in
advec_particles after the end of this loop.

advec_particles, check_parameters, data_output_dvrp, header, init_3d_model, init_grid, init_particles, init_pegrid, modules, package_parin, parin, plant_canopy_model, read_var_list, read_3d_binary, user_interface, write_var_list, write_3d_binary

Changed:


Redefine initial nzb_local as the actual total size of topography (later the
extent of topography in nzb_local is reduced by 1dx at the E topography walls
and by 1dy at the N topography walls to form the basis for nzb_s_inner);
for consistency redefine 'single_building' case.

Vertical profiles now based on nzb_s_inner; they are divided by
ngp_2dh_s_inner (scalars, procucts of scalars) and ngp_2dh (staggered velocity
components and their products, procucts of scalars and velocity components),
respectively.

Allow two instead of one digit to specify isosurface and slicer variables.

Status of 3D-volume NetCDF data file only depends on switch netcdf_64bit_3d (check_open)

prognostic_equations include the respective wall_*flux in the parameter list of
calls of diffusion_s. Same as before, only the values of wall_heatflux(0:4)
can be assigned. At present, wall_humidityflux, wall_qflux, wall_salinityflux,
and wall_scalarflux are kept zero. diffusion_s uses the respective wall_*flux
instead of wall_heatflux. This update serves two purposes:

  • it avoids errors in calculations with humidity/scalar/salinity and prescribed

non-zero wall_heatflux,

  • it prepares PALM for a possible assignment of wall fluxes of

humidity/scalar/salinity in a future release.

buoyancy, check_open, data_output_dvrp, diffusion_s, diffusivities, flow_statistics, header, init_3d_model, init_dvrp, init_grid, modules, prognostic_equations

Errors:


Bugfix: summation of sums_l_l in diffusivities.

Several bugfixes in the ocean part: Initial density rho is calculated
(init_ocean). Error in initializing u_init and v_init removed
(check_parameters). Calculation of density flux now starts from
nzb+1 (production_e).

Bugfix: pleft/pright changed to pnorth/psouth in sendrecv of particle tail
numbers along y, small bugfixes in the SGS part (advec_particles)

Bugfix: model_string needed a default value (combine_plot_fields)

Bugfix: wavenumber calculation for even nx in routines maketri (poisfft)

Bugfix: assignment of fluxes at walls

Bugfix: absolute value of f must be used when calculating the Blackadar mixing length (init_1d_model)

advec_particles, check_parameters, combine_plot_fields, diffusion_s, diffusivities, init_ocean, init_1d_model, poisfft, production_e

  • Property svn:keywords set to Id
File size: 9.8 KB
Line 
1 MODULE buoyancy_mod
2
3!------------------------------------------------------------------------------!
4! Actual revisions:
5! -----------------
6!
7!
8! Former revisions:
9! -----------------
10! $Id: buoyancy.f90 139 2007-11-29 09:37:41Z raasch $
11!
12! 132 2007-11-20 09:46:11Z letzel
13! Vertical scalar profiles now based on nzb_s_inner and ngp_2dh_s_inner.
14!
15! 106 2007-08-16 14:30:26Z raasch
16! i loop for u-component (sloping surface) is starting from nxlu (needed for
17! non-cyclic boundary conditions)
18!
19! 97 2007-06-21 08:23:15Z raasch
20! Routine reneralized to be used with temperature AND density:
21! argument theta renamed var, new argument var_reference,
22! use_pt_reference renamed use_reference,
23! calc_mean_pt_profile renamed calc_mean_profile
24!
25! 57 2007-03-09 12:05:41Z raasch
26! Reference temperature pt_reference can be used.
27!
28! RCS Log replace by Id keyword, revision history cleaned up
29!
30! Revision 1.19  2006/04/26 12:09:56  raasch
31! OpenMP optimization (one dimension added to sums_l)
32!
33! Revision 1.1  1997/08/29 08:56:48  raasch
34! Initial revision
35!
36!
37! Description:
38! ------------
39! Buoyancy term of the third component of the equation of motion.
40! WARNING: humidity is not regarded when using a sloping surface!
41!------------------------------------------------------------------------------!
42
43    PRIVATE
44    PUBLIC buoyancy, calc_mean_profile
45
46    INTERFACE buoyancy
47       MODULE PROCEDURE buoyancy
48       MODULE PROCEDURE buoyancy_ij
49    END INTERFACE buoyancy
50
51    INTERFACE calc_mean_profile
52       MODULE PROCEDURE calc_mean_profile
53    END INTERFACE calc_mean_profile
54
55 CONTAINS
56
57
58!------------------------------------------------------------------------------!
59! Call for all grid points
60!------------------------------------------------------------------------------!
61    SUBROUTINE buoyancy( var, var_reference, wind_component, pr )
62
63       USE arrays_3d
64       USE control_parameters
65       USE indices
66       USE pegrid
67       USE statistics
68
69       IMPLICIT NONE
70
71       INTEGER ::  i, j, k, pr, wind_component
72       REAL    ::  var_reference
73       REAL, DIMENSION(:,:,:), POINTER  ::  var
74
75
76       IF ( .NOT. sloping_surface )  THEN
77!
78!--       Normal case: horizontal surface
79          IF ( use_reference )  THEN
80             DO  i = nxl, nxr
81                DO  j = nys, nyn
82                   DO  k = nzb_s_inner(j,i)+1, nzt-1
83                      tend(k,j,i) = tend(k,j,i) + atmos_ocean_sign * g * 0.5 * &
84                                                            (                  &
85                          ( var(k,j,i)   - hom(k,1,pr,0)   ) / var_reference + &
86                          ( var(k+1,j,i) - hom(k+1,1,pr,0) ) / var_reference   &
87                                                            )
88                   ENDDO
89                ENDDO
90             ENDDO
91          ELSE
92             DO  i = nxl, nxr
93                DO  j = nys, nyn
94                   DO  k = nzb_s_inner(j,i)+1, nzt-1
95                      tend(k,j,i) = tend(k,j,i) + atmos_ocean_sign * g * 0.5 * &
96                                                            (                  &
97                          ( var(k,j,i)   - hom(k,1,pr,0)   ) / hom(k,1,pr,0) + &
98                          ( var(k+1,j,i) - hom(k+1,1,pr,0) ) / hom(k+1,1,pr,0) &
99                                                            )
100                   ENDDO
101                ENDDO
102             ENDDO
103          ENDIF
104
105       ELSE
106!
107!--       Buoyancy term for a surface with a slope in x-direction. The equations
108!--       for both the u and w velocity-component contain proportionate terms.
109!--       Temperature field at time t=0 serves as environmental temperature.
110!--       Reference temperature (pt_surface) is the one at the lower left corner
111!--       of the total domain.
112          IF ( wind_component == 1 )  THEN
113
114             DO  i = nxlu, nxr
115                DO  j = nys, nyn
116                   DO  k = nzb_s_inner(j,i)+1, nzt-1
117                      tend(k,j,i) = tend(k,j,i) + g * sin_alpha_surface *      &
118                           0.5 * ( ( pt(k,j,i-1)         + pt(k,j,i)         ) &
119                                 - ( pt_slope_ref(k,i-1) + pt_slope_ref(k,i) ) &
120                                 ) / pt_surface
121                   ENDDO
122                ENDDO
123             ENDDO
124
125          ELSEIF ( wind_component == 3 )  THEN
126
127             DO  i = nxl, nxr
128                DO  j = nys, nyn
129                   DO  k = nzb_s_inner(j,i)+1, nzt-1
130                      tend(k,j,i) = tend(k,j,i) + g * cos_alpha_surface *      &
131                           0.5 * ( ( pt(k,j,i)         + pt(k+1,j,i)         ) &
132                                 - ( pt_slope_ref(k,i) + pt_slope_ref(k+1,i) ) &
133                                 ) / pt_surface
134                   ENDDO
135                ENDDO
136            ENDDO
137
138          ELSE
139
140             IF ( myid == 0 )  PRINT*, '+++ buoyancy: no term for component "',&
141                                       wind_component,'"'
142             CALL local_stop
143
144          ENDIF
145
146       ENDIF
147
148    END SUBROUTINE buoyancy
149
150
151!------------------------------------------------------------------------------!
152! Call for grid point i,j
153!------------------------------------------------------------------------------!
154    SUBROUTINE buoyancy_ij( i, j, var, var_reference, wind_component, pr )
155
156       USE arrays_3d
157       USE control_parameters
158       USE indices
159       USE pegrid
160       USE statistics
161
162       IMPLICIT NONE
163
164       INTEGER ::  i, j, k, pr, wind_component
165       REAL    ::  var_reference
166       REAL, DIMENSION(:,:,:), POINTER  ::  var
167
168
169       IF ( .NOT. sloping_surface )  THEN
170!
171!--       Normal case: horizontal surface
172          IF ( use_reference )  THEN
173             DO  k = nzb_s_inner(j,i)+1, nzt-1
174                 tend(k,j,i) = tend(k,j,i) + atmos_ocean_sign * g * 0.5 * (   &
175                         ( var(k,j,i)   - hom(k,1,pr,0)   ) / var_reference + &
176                         ( var(k+1,j,i) - hom(k+1,1,pr,0) ) / var_reference   &
177                                                                          )
178             ENDDO
179          ELSE
180             DO  k = nzb_s_inner(j,i)+1, nzt-1
181                 tend(k,j,i) = tend(k,j,i) + atmos_ocean_sign * g * 0.5 * (    &
182                          ( var(k,j,i)   - hom(k,1,pr,0)   ) / hom(k,1,pr,0) + &
183                          ( var(k+1,j,i) - hom(k+1,1,pr,0) ) / hom(k+1,1,pr,0) &
184                                                                          )
185             ENDDO
186          ENDIF
187
188       ELSE
189!
190!--       Buoyancy term for a surface with a slope in x-direction. The equations
191!--       for both the u and w velocity-component contain proportionate terms.
192!--       Temperature field at time t=0 serves as environmental temperature.
193!--       Reference temperature (pt_surface) is the one at the lower left corner
194!--       of the total domain.
195          IF ( wind_component == 1 )  THEN
196
197             DO  k = nzb_s_inner(j,i)+1, nzt-1
198                tend(k,j,i) = tend(k,j,i) + g * sin_alpha_surface *            &
199                           0.5 * ( ( pt(k,j,i-1)         + pt(k,j,i)         ) &
200                                 - ( pt_slope_ref(k,i-1) + pt_slope_ref(k,i) ) &
201                                 ) / pt_surface
202             ENDDO
203
204          ELSEIF ( wind_component == 3 )  THEN
205
206             DO  k = nzb_s_inner(j,i)+1, nzt-1
207                tend(k,j,i) = tend(k,j,i) + g * cos_alpha_surface *            &
208                           0.5 * ( ( pt(k,j,i)         + pt(k+1,j,i)         ) &
209                                 - ( pt_slope_ref(k,i) + pt_slope_ref(k+1,i) ) &
210                                 ) / pt_surface
211             ENDDO
212
213          ELSE
214
215             IF ( myid == 0 )  PRINT*, '+++ buoyancy: no term for component "',&
216                                       wind_component,'"'
217             CALL local_stop
218
219          ENDIF
220
221       ENDIF
222
223    END SUBROUTINE buoyancy_ij
224
225
226    SUBROUTINE calc_mean_profile( var, pr )
227
228!------------------------------------------------------------------------------!
229! Description:
230! ------------
231! Calculate the horizontally averaged vertical temperature profile (pr=4 in case
232! of potential temperature and 44 in case of virtual potential temperature).
233!------------------------------------------------------------------------------!
234
235       USE control_parameters
236       USE indices
237       USE pegrid
238       USE statistics
239
240       IMPLICIT NONE
241
242       INTEGER ::  i, j, k, omp_get_thread_num, pr, tn
243       REAL, DIMENSION(:,:,:), POINTER  ::  var
244
245!
246!--    Computation of the horizontally averaged profile of variable var, unless
247!--    already done by the relevant call from flow_statistics. The calculation
248!--    is done only for the first respective intermediate timestep in order to
249!--    spare communication time and to produce identical model results with jobs
250!--    which are calling flow_statistics at different time intervals.
251       IF ( .NOT. flow_statistics_called  .AND. &
252            intermediate_timestep_count == 1 )  THEN
253
254!
255!--       Horizontal average of variable var
256          tn           =   0  ! Default thread number in case of one thread
257          !$OMP PARALLEL PRIVATE( i, j, k, tn )
258!$        tn = omp_get_thread_num()
259          sums_l(:,pr,tn) = 0.0
260          !$OMP DO
261          DO  i = nxl, nxr
262             DO  j =  nys, nyn
263                DO  k = nzb_s_inner(j,i), nzt+1
264                   sums_l(k,pr,tn) = sums_l(k,pr,tn) + var(k,j,i)
265                ENDDO
266             ENDDO
267          ENDDO
268          !$OMP END PARALLEL
269
270          DO  i = 1, threads_per_task-1
271             sums_l(:,pr,0) = sums_l(:,pr,0) + sums_l(:,pr,i)
272          ENDDO
273
274#if defined( __parallel )
275
276          CALL MPI_ALLREDUCE( sums_l(nzb,pr,0), sums(nzb,pr), nzt+2-nzb, &
277                              MPI_REAL, MPI_SUM, comm2d, ierr )
278
279#else
280
281          sums(:,pr) = sums_l(:,pr,0)
282
283#endif
284
285          hom(:,1,pr,0) = sums(:,pr) / ngp_2dh_s_inner(:,0)
286
287       ENDIF
288
289    END SUBROUTINE calc_mean_profile
290
291 END MODULE buoyancy_mod
Note: See TracBrowser for help on using the repository browser.