source: palm/trunk/SOURCE/init_3d_model.f90 @ 77

Last change on this file since 77 was 77, checked in by raasch, 14 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: 33.4 KB
Line 
1#if defined( __ibmy_special )
2@PROCESS NOOPTimize
3#endif
4 SUBROUTINE init_3d_model
5
6!------------------------------------------------------------------------------!
7! Actual revisions:
8! -----------------
9!
10!
11! Former revisions:
12! -----------------
13! $Id: init_3d_model.f90 77 2007-03-29 04:26:56Z raasch $
14!
15! 75 2007-03-22 09:54:05Z raasch
16! Arrays for radiation boundary conditions are allocated (u_m_l, u_m_r, etc.),
17! bugfix for cases with the outflow damping layer extending over more than one
18! subdomain, moisture renamed humidity,
19! new initializing action "by_user" calls user_init_3d_model,
20! precipitation_amount/rate, ts_value are allocated, +module netcdf_control,
21! initial velocities at nzb+1 are regarded for volume
22! flow control in case they have been set zero before (to avoid small timesteps)
23! -uvmean_outflow, uxrp, vynp eliminated
24!
25! 19 2007-02-23 04:53:48Z raasch
26! +handling of top fluxes
27!
28! RCS Log replace by Id keyword, revision history cleaned up
29!
30! Revision 1.49  2006/08/22 15:59:07  raasch
31! No optimization of this file on the ibmy (Yonsei Univ.)
32!
33! Revision 1.1  1998/03/09 16:22:22  raasch
34! Initial revision
35!
36!
37! Description:
38! ------------
39! Allocation of arrays and initialization of the 3D model via
40! a) pre-run the 1D model
41! or
42! b) pre-set constant linear profiles
43! or
44! c) read values of a previous run
45!------------------------------------------------------------------------------!
46
47    USE arrays_3d
48    USE averaging
49    USE cloud_parameters
50    USE constants
51    USE control_parameters
52    USE cpulog
53    USE indices
54    USE interfaces
55    USE model_1d
56    USE netcdf_control
57    USE particle_attributes
58    USE pegrid
59    USE profil_parameter
60    USE random_function_mod
61    USE statistics
62
63    IMPLICIT NONE
64
65    INTEGER ::  i, j, k, sr
66
67    INTEGER, DIMENSION(:), ALLOCATABLE ::  ngp_2dh_l, ngp_3d_inner_l
68
69    INTEGER, DIMENSION(:,:), ALLOCATABLE ::  ngp_2dh_outer_l
70
71    REAL, DIMENSION(1:2) ::  volume_flow_area_l, volume_flow_initial_l
72
73
74!
75!-- Allocate arrays
76    ALLOCATE( ngp_2dh(0:statistic_regions), ngp_2dh_l(0:statistic_regions), &
77              ngp_3d(0:statistic_regions),                                  &
78              ngp_3d_inner(0:statistic_regions),                            &
79              ngp_3d_inner_l(0:statistic_regions),                          &
80              sums_divnew_l(0:statistic_regions),                           &
81              sums_divold_l(0:statistic_regions) )
82    ALLOCATE( rdf(nzb+1:nzt) )
83    ALLOCATE( hom_sum(nzb:nzt+1,var_hom,0:statistic_regions),               &
84              ngp_2dh_outer(nzb:nzt+1,0:statistic_regions),                 &
85              ngp_2dh_outer_l(nzb:nzt+1,0:statistic_regions),               &
86              rmask(nys-1:nyn+1,nxl-1:nxr+1,0:statistic_regions),           &
87              sums(nzb:nzt+1,var_sum),                                      &
88              sums_l(nzb:nzt+1,var_sum,0:threads_per_task-1),               &
89              sums_l_l(nzb:nzt+1,0:statistic_regions,0:threads_per_task-1), &
90              sums_up_fraction_l(10,3,0:statistic_regions),                 &
91              sums_wsts_bc_l(nzb:nzt+1,0:statistic_regions),                &
92              ts_value(var_ts,0:statistic_regions) )
93    ALLOCATE( km_damp_x(nxl-1:nxr+1), km_damp_y(nys-1:nyn+1) )
94
95    ALLOCATE( rif_1(nys-1:nyn+1,nxl-1:nxr+1), shf_1(nys-1:nyn+1,nxl-1:nxr+1), &
96              ts(nys-1:nyn+1,nxl-1:nxr+1), tswst_1(nys-1:nyn+1,nxl-1:nxr+1),  &
97              us(nys-1:nyn+1,nxl-1:nxr+1), usws_1(nys-1:nyn+1,nxl-1:nxr+1),   &
98              vsws_1(nys-1:nyn+1,nxl-1:nxr+1), z0(nys-1:nyn+1,nxl-1:nxr+1) )
99
100    IF ( timestep_scheme(1:5) /= 'runge' )  THEN
101!
102!--    Leapfrog scheme needs two timelevels of diffusion quantities
103       ALLOCATE( rif_2(nys-1:nyn+1,nxl-1:nxr+1),   &
104                 shf_2(nys-1:nyn+1,nxl-1:nxr+1),   &
105                 tswst_2(nys-1:nyn+1,nxl-1:nxr+1), &
106                 usws_2(nys-1:nyn+1,nxl-1:nxr+1),  &
107                 vsws_2(nys-1:nyn+1,nxl-1:nxr+1) )
108    ENDIF
109
110    ALLOCATE( d(nzb+1:nzta,nys:nyna,nxl:nxra),         &
111              e_1(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1),  &
112              e_2(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1),  &
113              e_3(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1),  &
114              kh_1(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
115              km_1(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
116              p(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1),    &
117              pt_1(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
118              pt_2(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
119              pt_3(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
120              tend(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
121              u_1(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1),  &
122              u_2(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1),  &
123              u_3(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1),  &
124              v_1(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
125              v_2(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
126              v_3(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
127              w_1(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1),  &
128              w_2(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1),  &
129              w_3(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
130
131    IF ( timestep_scheme(1:5) /= 'runge' )  THEN
132       ALLOCATE( kh_2(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
133                 km_2(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
134    ENDIF
135
136    IF ( humidity  .OR.  passive_scalar ) THEN
137!
138!--    2D-humidity/scalar arrays
139       ALLOCATE ( qs(nys-1:nyn+1,nxl-1:nxr+1),     &
140                  qsws_1(nys-1:nyn+1,nxl-1:nxr+1), &
141                  qswst_1(nys-1:nyn+1,nxl-1:nxr+1) )
142
143       IF ( timestep_scheme(1:5) /= 'runge' )  THEN
144          ALLOCATE( qsws_2(nys-1:nyn+1,nxl-1:nxr+1), &
145                    qswst_2(nys-1:nyn+1,nxl-1:nxr+1) )
146       ENDIF
147!
148!--    3D-humidity/scalar arrays
149       ALLOCATE( q_1(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
150                 q_2(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
151                 q_3(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
152
153!
154!--    3D-arrays needed for humidity only
155       IF ( humidity )  THEN
156          ALLOCATE( vpt_1(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
157
158          IF ( timestep_scheme(1:5) /= 'runge' )  THEN
159             ALLOCATE( vpt_2(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
160          ENDIF
161
162          IF ( cloud_physics ) THEN
163!
164!--          Liquid water content
165             ALLOCATE ( ql_1(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
166!
167!--          Precipitation amount and rate (only needed if output is switched)
168             ALLOCATE( precipitation_amount(nys-1:nyn+1,nxl-1:nxr+1), &
169                       precipitation_rate(nys-1:nyn+1,nxl-1:nxr+1) )
170          ENDIF
171
172          IF ( cloud_droplets )  THEN
173!
174!--          Liquid water content, change in liquid water content,
175!--          real volume of particles (with weighting), volume of particles
176             ALLOCATE ( ql_1(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
177                        ql_2(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
178                        ql_v(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
179                        ql_vp(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
180          ENDIF
181
182       ENDIF
183
184    ENDIF
185
186!
187!-- 3D-array for storing the dissipation, needed for calculating the sgs
188!-- particle velocities
189    IF ( use_sgs_for_particles )  THEN
190       ALLOCATE ( diss(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
191    ENDIF
192
193    IF ( dt_dosp /= 9999999.9 )  THEN
194       ALLOCATE( spectrum_x( 1:nx/2, 1:10, 1:10 ), &
195                 spectrum_y( 1:ny/2, 1:10, 1:10 ) )
196    ENDIF
197
198!
199!-- 4D-array for storing the Rif-values at vertical walls
200    IF ( topography /= 'flat' )  THEN
201       ALLOCATE( rif_wall(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1,1:4) )
202       rif_wall = 0.0
203    ENDIF
204
205!
206!-- Velocities at nzb+1 needed for volume flow control
207    IF ( conserve_volume_flow )  THEN
208       ALLOCATE( u_nzb_p1_for_vfc(nys:nyn), v_nzb_p1_for_vfc(nxl:nxr) )
209       u_nzb_p1_for_vfc = 0.0
210       v_nzb_p1_for_vfc = 0.0
211    ENDIF
212
213!
214!-- Arrays to store velocity data from t-dt needed for radiation boundary
215!-- conditions
216    IF ( outflow_l )  THEN
217       ALLOCATE( u_m_l(nzb:nzt+1,nys-1:nyn+1,-1:1), &
218                 v_m_l(nzb:nzt+1,nys-1:nyn+1,-1:1), &
219                 w_m_l(nzb:nzt+1,nys-1:nyn+1,-1:1) )
220    ENDIF
221    IF ( outflow_r )  THEN
222       ALLOCATE( u_m_r(nzb:nzt+1,nys-1:nyn+1,nx-1:nx+1), &
223                 v_m_r(nzb:nzt+1,nys-1:nyn+1,nx-1:nx+1), &
224                 w_m_r(nzb:nzt+1,nys-1:nyn+1,nx-1:nx+1) )
225    ENDIF
226    IF ( outflow_s )  THEN
227       ALLOCATE( u_m_s(nzb:nzt+1,-1:1,nxl-1:nxr+1), &
228                 v_m_s(nzb:nzt+1,-1:1,nxl-1:nxr+1), &
229                 w_m_s(nzb:nzt+1,-1:1,nxl-1:nxr+1) )
230    ENDIF
231    IF ( outflow_n )  THEN
232       ALLOCATE( u_m_n(nzb:nzt+1,ny-1:ny+1,nxl-1:nxr+1), &
233                 v_m_n(nzb:nzt+1,ny-1:ny+1,nxl-1:nxr+1), &
234                 w_m_n(nzb:nzt+1,ny-1:ny+1,nxl-1:nxr+1) )
235    ENDIF
236
237!
238!-- Initial assignment of the pointers
239    IF ( timestep_scheme(1:5) /= 'runge' )  THEN
240
241       rif_m   => rif_1;    rif   => rif_2
242       shf_m   => shf_1;    shf   => shf_2
243       tswst_m => tswst_1;  tswst => tswst_2
244       usws_m  => usws_1;   usws  => usws_2
245       vsws_m  => vsws_1;   vsws  => vsws_2
246       e_m  => e_1;   e  => e_2;   e_p  => e_3;   te_m  => e_3
247       kh_m => kh_1;  kh => kh_2
248       km_m => km_1;  km => km_2
249       pt_m => pt_1;  pt => pt_2;  pt_p => pt_3;  tpt_m => pt_3
250       u_m  => u_1;   u  => u_2;   u_p  => u_3;   tu_m  => u_3
251       v_m  => v_1;   v  => v_2;   v_p  => v_3;   tv_m  => v_3
252       w_m  => w_1;   w  => w_2;   w_p  => w_3;   tw_m  => w_3
253
254       IF ( humidity  .OR.  passive_scalar )  THEN
255          qsws_m  => qsws_1;   qsws  => qsws_2
256          qswst_m => qswst_1;  qswst => qswst_2
257          q_m    => q_1;     q    => q_2;     q_p => q_3;     tq_m => q_3
258          IF ( humidity )        vpt_m  => vpt_1;   vpt  => vpt_2
259          IF ( cloud_physics )   ql   => ql_1
260          IF ( cloud_droplets )  THEN
261             ql   => ql_1
262             ql_c => ql_2
263          ENDIF
264       ENDIF
265
266    ELSE
267
268       rif   => rif_1
269       shf   => shf_1
270       tswst => tswst_1
271       usws  => usws_1
272       vsws  => vsws_1
273       e     => e_1;   e_p  => e_2;   te_m  => e_3;   e_m  => e_3
274       kh    => kh_1
275       km    => km_1
276       pt    => pt_1;  pt_p => pt_2;  tpt_m => pt_3;  pt_m => pt_3
277       u     => u_1;   u_p  => u_2;   tu_m  => u_3;   u_m  => u_3
278       v     => v_1;   v_p  => v_2;   tv_m  => v_3;   v_m  => v_3
279       w     => w_1;   w_p  => w_2;   tw_m  => w_3;   w_m  => w_3
280
281       IF ( humidity  .OR.  passive_scalar )  THEN
282          qsws   => qsws_1
283          qswst  => qswst_1
284          q      => q_1;     q_p  => q_2;     tq_m => q_3;    q_m => q_3
285          IF ( humidity )        vpt  => vpt_1
286          IF ( cloud_physics )   ql   => ql_1
287          IF ( cloud_droplets )  THEN
288             ql   => ql_1
289             ql_c => ql_2
290          ENDIF
291       ENDIF
292
293    ENDIF
294
295!
296!-- Initialize model variables
297    IF ( TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
298!
299!--    First model run of a possible job queue.
300!--    Initial profiles of the variables must be computes.
301       IF ( INDEX( initializing_actions, 'set_1d-model_profiles' ) /= 0 )  THEN
302!
303!--       Use solutions of the 1D model as initial profiles,
304!--       start 1D model
305          CALL init_1d_model
306!
307!--       Transfer initial profiles to the arrays of the 3D model
308          DO  i = nxl-1, nxr+1
309             DO  j = nys-1, nyn+1
310                e(:,j,i)  = e1d
311                kh(:,j,i) = kh1d
312                km(:,j,i) = km1d
313                pt(:,j,i) = pt_init
314                u(:,j,i)  = u1d
315                v(:,j,i)  = v1d
316             ENDDO
317          ENDDO
318
319          IF ( humidity  .OR.  passive_scalar )  THEN
320             DO  i = nxl-1, nxr+1
321                DO  j = nys-1, nyn+1
322                   q(:,j,i) = q_init
323                ENDDO
324             ENDDO
325          ENDIF
326
327          IF ( .NOT. constant_diffusion )  THEN
328             DO  i = nxl-1, nxr+1
329                DO  j = nys-1, nyn+1
330                   e(:,j,i)  = e1d
331                ENDDO
332             ENDDO
333!
334!--          Store initial profiles for output purposes etc.
335             hom(:,1,25,:) = SPREAD( l1d, 2, statistic_regions+1 )
336
337             IF ( prandtl_layer )  THEN
338                rif  = rif1d(nzb+1)
339                ts   = 0.0  ! could actually be computed more accurately in the
340                            ! 1D model. Update when opportunity arises.
341                us   = us1d
342                usws = usws1d
343                vsws = vsws1d
344             ELSE
345                ts   = 0.0  ! must be set, because used in
346                rif  = 0.0  ! flowste
347                us   = 0.0
348                usws = 0.0
349                vsws = 0.0
350             ENDIF
351
352          ELSE
353             e    = 0.0  ! must be set, because used in
354             rif  = 0.0  ! flowste
355             ts   = 0.0
356             us   = 0.0
357             usws = 0.0
358             vsws = 0.0
359          ENDIF
360
361!
362!--       In every case qs = 0.0 (see also pt)
363!--       This could actually be computed more accurately in the 1D model.
364!--       Update when opportunity arises!
365          IF ( humidity  .OR.  passive_scalar )  qs = 0.0
366
367!
368!--       inside buildings set velocities back to zero
369          IF ( topography /= 'flat' )  THEN
370             DO  i = nxl-1, nxr+1
371                DO  j = nys-1, nyn+1
372                   u(nzb:nzb_u_inner(j,i),j,i) = 0.0
373                   v(nzb:nzb_v_inner(j,i),j,i) = 0.0
374                ENDDO
375             ENDDO
376!
377!--          WARNING: The extra boundary conditions set after running the
378!--          -------  1D model impose an error on the divergence one layer
379!--                   below the topography; need to correct later
380!--          ATTENTION: Provisional correction for Piacsek & Williams
381!--          ---------  advection scheme: keep u and v zero one layer below
382!--                     the topography.
383             IF ( ibc_uv_b == 0 )  THEN
384!
385!--             Satisfying the Dirichlet condition with an extra layer below
386!--             the surface where the u and v component change their sign.
387                DO  i = nxl-1, nxr+1
388                   DO  j = nys-1, nyn+1
389                      IF ( nzb_u_inner(j,i) == 0 ) u(0,j,i) = -u(1,j,i)
390                      IF ( nzb_v_inner(j,i) == 0 ) v(0,j,i) = -v(1,j,i)
391                   ENDDO
392                ENDDO
393
394             ELSE
395!
396!--             Neumann condition
397                DO  i = nxl-1, nxr+1
398                   DO  j = nys-1, nyn+1
399                      IF ( nzb_u_inner(j,i) == 0 ) u(0,j,i) = u(1,j,i)
400                      IF ( nzb_v_inner(j,i) == 0 ) v(0,j,i) = v(1,j,i)
401                   ENDDO
402                ENDDO
403
404             ENDIF
405
406          ENDIF
407
408       ELSEIF ( INDEX(initializing_actions, 'set_constant_profiles') /= 0 ) &
409       THEN
410!
411!--       Use constructed initial profiles (velocity constant with height,
412!--       temperature profile with constant gradient)
413          DO  i = nxl-1, nxr+1
414             DO  j = nys-1, nyn+1
415                pt(:,j,i) = pt_init
416                u(:,j,i)  = u_init
417                v(:,j,i)  = v_init
418             ENDDO
419          ENDDO
420
421!
422!--       Set initial horizontal velocities at the lowest computational grid levels
423!--       to zero in order to avoid too small time steps caused by the diffusion
424!--       limit in the initial phase of a run (at k=1, dz/2 occurs in the
425!--       limiting formula!). The original values are stored to be later used for
426!--       volume flow control.
427          DO  i = nxl-1, nxr+1
428             DO  j = nys-1, nyn+1
429                u(nzb:nzb_u_inner(j,i)+1,j,i) = 0.0
430                v(nzb:nzb_v_inner(j,i)+1,j,i) = 0.0
431             ENDDO
432          ENDDO
433          IF ( conserve_volume_flow )  THEN
434             IF ( nxr == nx )  THEN
435                DO  j = nys, nyn
436                   k = nzb_u_inner(j,nx) + 1
437                   u_nzb_p1_for_vfc(j) = u_init(k) * dzu(k)
438                ENDDO
439             ENDIF
440             IF ( nyn == ny )  THEN
441                DO  i = nxl, nxr
442                   k = nzb_v_inner(ny,i) + 1
443                   v_nzb_p1_for_vfc(i) = v_init(k) * dzu(k)
444                ENDDO
445             ENDIF
446          ENDIF
447
448          IF ( humidity  .OR.  passive_scalar )  THEN
449             DO  i = nxl-1, nxr+1
450                DO  j = nys-1, nyn+1
451                   q(:,j,i) = q_init
452                ENDDO
453             ENDDO
454          ENDIF
455
456         
457          IF ( constant_diffusion )  THEN
458             km   = km_constant
459             kh   = km / prandtl_number
460          ELSE
461             kh   = 0.01   ! there must exist an initial diffusion, because
462             km   = 0.01   ! otherwise no TKE would be produced by the
463                           ! production terms, as long as not yet
464                           ! e = (u*/cm)**2 at k=nzb+1
465          ENDIF
466          e    = 0.0
467          rif  = 0.0
468          ts   = 0.0
469          us   = 0.0
470          usws = 0.0
471          vsws = 0.0
472          IF ( humidity  .OR.  passive_scalar ) qs = 0.0
473
474!
475!--       Compute initial temperature field and other constants used in case
476!--       of a sloping surface
477          IF ( sloping_surface )  CALL init_slope
478
479       ELSEIF ( INDEX(initializing_actions, 'by_user') /= 0 ) &
480       THEN
481!
482!--       Initialization will completely be done by the user
483          CALL user_init_3d_model
484
485       ENDIF
486
487!
488!--    Calculate virtual potential temperature
489       IF ( humidity ) vpt = pt * ( 1.0 + 0.61 * q )
490
491!
492!--    Store initial profiles for output purposes etc.
493       hom(:,1,5,:) = SPREAD( u(:,nys,nxl), 2, statistic_regions+1 )
494       hom(:,1,6,:) = SPREAD( v(:,nys,nxl), 2, statistic_regions+1 )
495       IF ( ibc_uv_b == 0 )  THEN
496          hom(nzb,1,5,:) = -hom(nzb+1,1,5,:)  ! due to satisfying the Dirichlet
497          hom(nzb,1,6,:) = -hom(nzb+1,1,6,:)  ! condition with an extra layer
498              ! below the surface where the u and v component change their sign
499       ENDIF
500       hom(:,1,7,:)  = SPREAD( pt(:,nys,nxl), 2, statistic_regions+1 )
501       hom(:,1,23,:) = SPREAD( km(:,nys,nxl), 2, statistic_regions+1 )
502       hom(:,1,24,:) = SPREAD( kh(:,nys,nxl), 2, statistic_regions+1 )
503
504
505       IF ( humidity )  THEN
506!
507!--       Store initial profile of total water content, virtual potential
508!--       temperature
509          hom(:,1,26,:) = SPREAD(   q(:,nys,nxl), 2, statistic_regions+1 )
510          hom(:,1,29,:) = SPREAD( vpt(:,nys,nxl), 2, statistic_regions+1 )
511          IF ( cloud_physics  .OR.  cloud_droplets ) THEN
512!
513!--          Store initial profile of specific humidity and potential
514!--          temperature
515             hom(:,1,27,:) = SPREAD(  q(:,nys,nxl), 2, statistic_regions+1 )
516             hom(:,1,28,:) = SPREAD( pt(:,nys,nxl), 2, statistic_regions+1 )
517          ENDIF
518       ENDIF
519
520       IF ( passive_scalar )  THEN
521!
522!--       Store initial scalar profile
523          hom(:,1,26,:) = SPREAD(  q(:,nys,nxl), 2, statistic_regions+1 )
524       ENDIF
525
526!
527!--    Initialize fluxes at bottom surface
528       IF ( use_surface_fluxes )  THEN
529
530          IF ( constant_heatflux )  THEN
531!
532!--          Heat flux is prescribed
533             IF ( random_heatflux )  THEN
534                CALL disturb_heatflux
535             ELSE
536                shf = surface_heatflux
537!
538!--             Over topography surface_heatflux is replaced by wall_heatflux(0)
539                IF ( TRIM( topography ) /= 'flat' )  THEN
540                   DO  i = nxl-1, nxr+1
541                      DO  j = nys-1, nyn+1
542                         IF ( nzb_s_inner(j,i) /= 0 )  THEN
543                            shf(j,i) = wall_heatflux(0)
544                         ENDIF
545                      ENDDO
546                   ENDDO
547                ENDIF
548             ENDIF
549             IF ( ASSOCIATED( shf_m ) )  shf_m = shf
550          ENDIF
551
552!
553!--       Determine the near-surface water flux
554          IF ( humidity  .OR.  passive_scalar )  THEN
555             IF ( constant_waterflux )  THEN
556                qsws   = surface_waterflux
557                IF ( ASSOCIATED( qsws_m ) )  qsws_m = qsws
558             ENDIF
559          ENDIF
560
561       ENDIF
562
563!
564!--    Initialize fluxes at top surface
565!--    Currently, only the heatflux can be prescribed. The latent flux is
566!--    zero in this case!
567       IF ( use_top_fluxes )  THEN
568
569          IF ( constant_top_heatflux )  THEN
570!
571!--          Heat flux is prescribed
572             tswst = top_heatflux
573             IF ( ASSOCIATED( tswst_m ) )  tswst_m = tswst
574
575             IF ( humidity  .OR.  passive_scalar )  THEN
576                qswst = 0.0
577                IF ( ASSOCIATED( qswst_m ) )  qswst_m = qswst
578             ENDIF
579         ENDIF
580
581       ENDIF
582
583!
584!--    Initialize Prandtl layer quantities
585       IF ( prandtl_layer )  THEN
586
587          z0 = roughness_length
588
589          IF ( .NOT. constant_heatflux )  THEN 
590!
591!--          Surface temperature is prescribed. Here the heat flux cannot be
592!--          simply estimated, because therefore rif, u* and theta* would have
593!--          to be computed by iteration. This is why the heat flux is assumed
594!--          to be zero before the first time step. It approaches its correct
595!--          value in the course of the first few time steps.
596             shf   = 0.0
597             IF ( ASSOCIATED( shf_m ) )  shf_m = 0.0
598          ENDIF
599
600          IF ( humidity  .OR.  passive_scalar )  THEN
601             IF ( .NOT. constant_waterflux )  THEN
602                qsws   = 0.0
603                IF ( ASSOCIATED( qsws_m ) )   qsws_m = 0.0
604             ENDIF
605          ENDIF
606
607       ENDIF
608
609!
610!--    Calculate the initial volume flow at the right and north boundary
611       IF ( conserve_volume_flow )  THEN
612
613          volume_flow_initial_l = 0.0
614          volume_flow_area_l    = 0.0
615 
616          IF ( nxr == nx )  THEN
617             DO  j = nys, nyn
618                DO  k = nzb_2d(j,nx) + 1, nzt
619                   volume_flow_initial_l(1) = volume_flow_initial_l(1) + &
620                                              u(k,j,nx) * dzu(k)
621                   volume_flow_area_l(1)    = volume_flow_area_l(1) + dzu(k)
622                ENDDO
623!
624!--             Correction if velocity at nzb+1 has been set zero further above
625                volume_flow_initial_l(1) = volume_flow_initial_l(1) + &
626                                           u_nzb_p1_for_vfc(j)
627             ENDDO
628          ENDIF
629
630          IF ( nyn == ny )  THEN
631             DO  i = nxl, nxr
632                DO  k = nzb_2d(ny,i) + 1, nzt 
633                   volume_flow_initial_l(2) = volume_flow_initial_l(2) + &
634                                              v(k,ny,i) * dzu(k)
635                   volume_flow_area_l(2)    = volume_flow_area_l(2) + dzu(k)
636                ENDDO
637!
638!--             Correction if velocity at nzb+1 has been set zero further above
639                volume_flow_initial_l(2) = volume_flow_initial_l(2) + &
640                                           v_nzb_p1_for_vfc(i)
641             ENDDO
642          ENDIF
643
644#if defined( __parallel )
645          CALL MPI_ALLREDUCE( volume_flow_initial_l(1), volume_flow_initial(1),&
646                              2, MPI_REAL, MPI_SUM, comm2d, ierr )
647          CALL MPI_ALLREDUCE( volume_flow_area_l(1), volume_flow_area(1),      &
648                              2, MPI_REAL, MPI_SUM, comm2d, ierr )
649#else
650          volume_flow_initial = volume_flow_initial_l
651          volume_flow_area    = volume_flow_area_l
652#endif 
653       ENDIF
654
655!
656!--    For the moment, perturbation pressure and vertical velocity are zero
657       p = 0.0; w = 0.0
658
659!
660!--    Initialize array sums (must be defined in first call of pres)
661       sums = 0.0
662
663!
664!--    Treating cloud physics, liquid water content and precipitation amount
665!--    are zero at beginning of the simulation
666       IF ( cloud_physics )  THEN
667          ql = 0.0
668          IF ( precipitation )  precipitation_amount = 0.0
669       ENDIF
670
671!
672!--    Initialize spectra
673       IF ( dt_dosp /= 9999999.9 )  THEN
674          spectrum_x = 0.0
675          spectrum_y = 0.0
676       ENDIF
677
678!
679!--    Impose vortex with vertical axis on the initial velocity profile
680       IF ( INDEX( initializing_actions, 'initialize_vortex' ) /= 0 )  THEN
681          CALL init_rankine
682       ENDIF
683
684!
685!--    Impose temperature anomaly (advection test only)
686       IF ( INDEX( initializing_actions, 'initialize_ptanom' ) /= 0 )  THEN
687          CALL init_pt_anomaly
688       ENDIF
689
690!
691!--    If required, change the surface temperature at the start of the 3D run
692       IF ( pt_surface_initial_change /= 0.0 )  THEN
693          pt(nzb,:,:) = pt(nzb,:,:) + pt_surface_initial_change
694       ENDIF
695
696!
697!--    If required, change the surface humidity/scalar at the start of the 3D
698!--    run
699       IF ( ( humidity .OR. passive_scalar ) .AND. &
700            q_surface_initial_change /= 0.0 )  THEN
701          q(nzb,:,:) = q(nzb,:,:) + q_surface_initial_change
702       ENDIF
703
704!
705!--    Initialize the random number generator (from numerical recipes)
706       CALL random_function_ini
707
708!
709!--    Impose random perturbation on the horizontal velocity field and then
710!--    remove the divergences from the velocity field
711       IF ( create_disturbances )  THEN
712          CALL disturb_field( nzb_u_inner, tend, u )
713          CALL disturb_field( nzb_v_inner, tend, v )
714          n_sor = nsor_ini
715          CALL pres
716          n_sor = nsor
717       ENDIF
718
719!
720!--    Once again set the perturbation pressure explicitly to zero in order to
721!--    assure that it does not generate any divergences in the first time step.
722!--    At t=0 the velocity field is free of divergence (as constructed above).
723!--    Divergences being created during a time step are not yet known and thus
724!--    cannot be corrected during the time step yet.
725       p = 0.0
726
727!
728!--    Initialize old and new time levels.
729       IF ( timestep_scheme(1:5) /= 'runge' )  THEN
730          e_m = e; pt_m = pt; u_m = u; v_m = v; w_m = w; kh_m = kh; km_m = km
731       ELSE
732          te_m = 0.0; tpt_m = 0.0; tu_m = 0.0; tv_m = 0.0; tw_m = 0.0
733       ENDIF
734       e_p = e; pt_p = pt; u_p = u; v_p = v; w_p = w
735
736       IF ( humidity  .OR.  passive_scalar )  THEN
737          IF ( ASSOCIATED( q_m ) )  q_m = q
738          IF ( timestep_scheme(1:5) == 'runge' )  tq_m = 0.0
739          q_p = q
740          IF ( humidity  .AND.  ASSOCIATED( vpt_m ) )  vpt_m = vpt
741       ENDIF
742
743!
744!--    Initialize old timelevels needed for radiation boundary conditions
745       IF ( outflow_l )  THEN
746          u_m_l(:,:,:) = u(:,:,-1:1)
747          v_m_l(:,:,:) = v(:,:,-1:1)
748          w_m_l(:,:,:) = w(:,:,-1:1)
749       ENDIF
750       IF ( outflow_r )  THEN
751          u_m_r(:,:,:) = u(:,:,nx-1:nx+1)
752          v_m_r(:,:,:) = v(:,:,nx-1:nx+1)
753          w_m_r(:,:,:) = w(:,:,nx-1:nx+1)
754       ENDIF
755       IF ( outflow_s )  THEN
756          u_m_s(:,:,:) = u(:,-1:1,:)
757          v_m_s(:,:,:) = v(:,-1:1,:)
758          w_m_s(:,:,:) = w(:,-1:1,:)
759       ENDIF
760       IF ( outflow_n )  THEN
761          u_m_n(:,:,:) = u(:,ny-1:ny+1,:)
762          v_m_n(:,:,:) = v(:,ny-1:ny+1,:)
763          w_m_n(:,:,:) = w(:,ny-1:ny+1,:)
764       ENDIF
765
766    ELSEIF ( TRIM( initializing_actions ) == 'read_restart_data' ) &
767    THEN
768!
769!--    Read binary data from restart file
770       CALL read_3d_binary
771
772!
773!--    Calculate initial temperature field and other constants used in case
774!--    of a sloping surface
775       IF ( sloping_surface )  CALL init_slope
776
777!
778!--    Initialize new time levels (only done in order to set boundary values
779!--    including ghost points)
780       e_p = e; pt_p = pt; u_p = u; v_p = v; w_p = w
781       IF ( humidity  .OR.  passive_scalar )  q_p = q
782
783    ELSE
784!
785!--    Actually this part of the programm should not be reached
786       IF ( myid == 0 )  PRINT*,'+++ init_3d_model: unknown initializing ', &
787                                                    'problem'
788       CALL local_stop
789    ENDIF
790
791!
792!-- If required, initialize dvrp-software
793    IF ( dt_dvrp /= 9999999.9 )  CALL init_dvrp
794
795!
796!-- If required, initialize quantities for handling cloud physics
797!-- This routine must be called before init_particles, because
798!-- otherwise, array pt_d_t, needed in data_output_dvrp (called by
799!-- init_particles) is not defined.
800    CALL init_cloud_physics
801
802!
803!-- If required, initialize particles
804    IF ( particle_advection )  CALL init_particles
805
806!
807!-- Initialize quantities for special advections schemes
808    CALL init_advec
809
810!
811!-- Initialize Rayleigh damping factors
812    rdf = 0.0
813    IF ( rayleigh_damping_factor /= 0.0 )  THEN
814       DO  k = nzb+1, nzt
815          IF ( zu(k) >= rayleigh_damping_height )  THEN
816             rdf(k) = rayleigh_damping_factor * &
817                      ( SIN( pi * 0.5 * ( zu(k) - rayleigh_damping_height )    &
818                                      / ( zu(nzt) - rayleigh_damping_height ) )&
819                      )**2
820          ENDIF
821       ENDDO
822    ENDIF
823
824!
825!-- Initialize diffusivities used within the outflow damping layer in case of
826!-- non-cyclic lateral boundaries. A linear increase is assumed over the first
827!-- half of the width of the damping layer
828    IF ( bc_lr == 'dirichlet/radiation' )  THEN
829
830       DO  i = nxl-1, nxr+1
831          IF ( i >= nx - outflow_damping_width )  THEN
832             km_damp_x(i) = km_damp_max * MIN( 1.0,                    &
833                            ( i - ( nx - outflow_damping_width ) ) /   &
834                            REAL( outflow_damping_width/2 )            &
835                                             )
836          ELSE
837             km_damp_x(i) = 0.0
838          ENDIF
839       ENDDO
840
841    ELSEIF ( bc_lr == 'radiation/dirichlet' )  THEN
842
843       DO  i = nxl-1, nxr+1
844          IF ( i <= outflow_damping_width )  THEN
845             km_damp_x(i) = km_damp_max * MIN( 1.0,                    &
846                            ( outflow_damping_width - i ) /            &
847                            REAL( outflow_damping_width/2 )            &
848                                             )
849          ELSE
850             km_damp_x(i) = 0.0
851          ENDIF
852       ENDDO
853
854    ENDIF
855
856    IF ( bc_ns == 'dirichlet/radiation' )  THEN
857
858       DO  j = nys-1, nyn+1
859          IF ( j >= ny - outflow_damping_width )  THEN
860             km_damp_y(j) = km_damp_max * MIN( 1.0,                    &
861                            ( j - ( ny - outflow_damping_width ) ) /   &
862                            REAL( outflow_damping_width/2 )            &
863                                             )
864          ELSE
865             km_damp_y(j) = 0.0
866          ENDIF
867       ENDDO
868
869    ELSEIF ( bc_ns == 'radiation/dirichlet' )  THEN
870
871       DO  j = nys-1, nyn+1
872          IF ( j <= outflow_damping_width )  THEN
873             km_damp_y(j) = km_damp_max * MIN( 1.0,                    &
874                            ( outflow_damping_width - j ) /            &
875                            REAL( outflow_damping_width/2 )            &
876                                             )
877          ELSE
878             km_damp_y(j) = 0.0
879          ENDIF
880       ENDDO
881
882    ENDIF
883
884!
885!-- Initialize local summation arrays for UP flow_statistics. This is necessary
886!-- because they may not yet have been initialized when they are called from
887!-- flow_statistics (or - depending on the chosen model run - are never
888!-- initialized)
889    sums_divnew_l      = 0.0
890    sums_divold_l      = 0.0
891    sums_l_l           = 0.0
892    sums_up_fraction_l = 0.0
893    sums_wsts_bc_l     = 0.0
894
895!
896!-- Pre-set masks for regional statistics. Default is the total model domain.
897    rmask = 1.0
898
899!
900!-- User-defined initializing actions. Check afterwards, if maximum number
901!-- of allowed timeseries is not exceeded
902    CALL user_init
903
904    IF ( dots_num > dots_max )  THEN
905       IF ( myid == 0 )  THEN
906          PRINT*, '+++ user_init: number of time series quantities exceeds', &
907                  ' its maximum of dots_max = ', dots_max
908          PRINT*, '    Please increase dots_max in modules.f90.'
909       ENDIF
910       CALL local_stop
911    ENDIF
912
913!
914!-- Input binary data file is not needed anymore. This line must be placed
915!-- after call of user_init!
916    CALL close_file( 13 )
917
918!
919!-- Compute total sum of active mask grid points
920!-- ngp_2dh: number of grid points of a horizontal cross section through the
921!--          total domain
922!-- ngp_3d:  number of grid points of the total domain
923!-- Note: The lower vertical index nzb_s_outer imposes a small error on the 2D
924!-- ----  averages of staggered variables such as u and v due to the topography
925!--       arrangement on the staggered grid. Maybe revise later.
926    ngp_2dh_outer_l = 0
927    ngp_2dh_outer   = 0
928    ngp_2dh_l       = 0
929    ngp_2dh         = 0
930    ngp_3d_inner_l  = 0
931    ngp_3d_inner    = 0
932    ngp_3d          = 0
933    ngp_sums        = ( nz + 2 ) * var_sum
934
935    DO  sr = 0, statistic_regions
936       DO  i = nxl, nxr
937          DO  j = nys, nyn
938             IF ( rmask(j,i,sr) == 1.0 )  THEN
939!
940!--             All xy-grid points
941                ngp_2dh_l(sr) = ngp_2dh_l(sr) + 1
942!
943!--             xy-grid points above topography
944                DO  k = nzb_s_outer(j,i), nz + 1
945                   ngp_2dh_outer_l(k,sr) = ngp_2dh_outer_l(k,sr) + 1
946                ENDDO
947!
948!--             All grid points of the total domain above topography
949                ngp_3d_inner_l(sr) = ngp_3d_inner_l(sr) + &
950                                     ( nz - nzb_s_inner(j,i) + 2 )
951             ENDIF
952          ENDDO
953       ENDDO
954    ENDDO
955
956    sr = statistic_regions + 1
957#if defined( __parallel )
958    CALL MPI_ALLREDUCE( ngp_2dh_l(0), ngp_2dh(0), sr, MPI_INTEGER, MPI_SUM,  &
959                        comm2d, ierr )
960    CALL MPI_ALLREDUCE( ngp_2dh_outer_l(0,0), ngp_2dh_outer(0,0), (nz+2)*sr, &
961                        MPI_INTEGER, MPI_SUM, comm2d, ierr )
962    CALL MPI_ALLREDUCE( ngp_3d_inner_l(0), ngp_3d_inner(0), sr, MPI_INTEGER, &
963                        MPI_SUM, comm2d, ierr )
964#else
965    ngp_2dh       = ngp_2dh_l
966    ngp_2dh_outer = ngp_2dh_outer_l
967    ngp_3d_inner  = ngp_3d_inner_l
968#endif
969
970    ngp_3d = ngp_2dh * ( nz + 2 )
971
972!
973!-- Set a lower limit of 1 in order to avoid zero divisions in flow_statistics,
974!-- buoyancy, etc. A zero value will occur for cases where all grid points of
975!-- the respective subdomain lie below the surface topography
976    ngp_2dh_outer = MAX( 1, ngp_2dh_outer(:,:) ) 
977    ngp_3d_inner  = MAX( 1, ngp_3d_inner(:)    )
978
979    DEALLOCATE( ngp_2dh_l, ngp_2dh_outer_l, ngp_3d_inner_l )
980
981
982 END SUBROUTINE init_3d_model
Note: See TracBrowser for help on using the repository browser.