Home


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

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

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