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

Last change on this file since 707 was 707, checked in by raasch, 11 years ago

New:
---

In case of multigrid method, on coarse grid levels, gathered data are
identically processed on all PEs (before, on PE0 only), so that the subsequent
scattering of data is not neccessary any more. (modules, init_pegrid, poismg)

Changed:


Calculation of weighted average of p is now handled in the same way
regardless of the number of ghost layers (advection scheme). (pres)

multigrid and sor method are using p_loc for iterative
advancements of pressure. p_sub removed. (init_3d_model, modules, poismg, pres, sor)

bc_lr and bc_ns replaced by bc_lr_dirrad, bc_lr_raddir, bc_ns_dirrad, bc_ns_raddir
for speed optimization. (calc_spectra, check_parameters, exchange_horiz,
exchange_horiz_2d, header, init_3d_model, init_grid, init_pegrid, modules,
poismg, pres, sor, time_integration, timestep)

grid_level directly used as index for MPI data type arrays. (exchange_horiz,
poismg)

initial assignments of zero to array p for iterative solvers only (init_3d_model)

Errors:


localsum calculation modified for proper OpenMP reduction. (pres)

Bugfix: bottom (nzb) and top (nzt+1) boundary conditions set in routines
resid and restrict. They were missed before, which may have led to
unpredictable results. (poismg)

  • Property svn:keywords set to Id
File size: 55.3 KB
Line 
1#if defined( __ibmy_special )
2@PROCESS NOOPTimize
3#endif
4 SUBROUTINE init_3d_model
5
6!------------------------------------------------------------------------------!
7! Current revisions:
8! -----------------
9! p_sub renamed p_loc and allocated depending on the chosen pressure solver,
10! initial assignments of zero to array p for iterative solvers only,
11! bc_lr/ns replaced by bc_lr/ns_dirrad/raddir
12!
13! Former revisions:
14! -----------------
15! 680 2011-02-04 23:16:06Z gryschka $
16! bugfix: volume_flow_control
17!
18! 673 2011-01-18 16:19:48Z suehring
19! weight_substep (moved from advec_ws) and weight_pres added.
20! Allocate p_sub when using Multigrid or SOR solver.
21! Call of ws_init moved behind the if requests.
22!
23! 667 2010-12-23 12:06:00Z suehring/gryschka
24! nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng in loops and
25! allocation of arrays. Calls of exchange_horiz are modified.
26! Call ws_init to initialize arrays needed for statistical evaluation and
27! optimization when ws-scheme is used.
28! Initial volume flow is now calculated by using the variable hom_sum.
29! Therefore the correction of initial volume flow for non-flat topography
30! removed (removed u_nzb_p1_for_vfc and v_nzb_p1_for_vfc)
31! Changed surface boundary conditions for u and v in case of ibc_uv_b == 0 from
32! mirror bc to dirichlet boundary conditions (u=v=0), so that k=nzb is
33! representative for the height z0
34! Bugfix: type conversion of '1' to 64bit for the MAX function (ngp_3d_inner)
35!
36! 622 2010-12-10 08:08:13Z raasch
37! optional barriers included in order to speed up collective operations
38!
39! 560 2010-09-09 10:06:09Z weinreis
40! bugfix: correction of calculating ngp_3d for 64 bit
41!
42! 485 2010-02-05 10:57:51Z raasch
43! calculation of ngp_3d + ngp_3d_inner changed because they have now 64 bit
44!
45! 407 2009-12-01 15:01:15Z maronga
46! var_ts is replaced by dots_max
47! Enabled passive scalar/humidity wall fluxes for non-flat topography
48!
49! 388 2009-09-23 09:40:33Z raasch
50! Initialization of prho added.
51! bugfix: correction of initial volume flow for non-flat topography
52! bugfix: zero initialization of arrays within buildings for 'cyclic_fill'
53! bugfix: avoid that ngp_2dh_s_inner becomes zero
54! initializing_actions='read_data_for_recycling' renamed to 'cyclic_fill', now
55! independent of turbulent_inflow
56! Output of messages replaced by message handling routine.
57! Set the starting level and the vertical smoothing factor used for
58! the external pressure gradient
59! +conserve_volume_flow_mode: 'default', 'initial_profiles', 'inflow_profile'
60! and 'bulk_velocity'
61! If the inversion height calculated by the prerun is zero,
62! inflow_damping_height must be explicitly specified.
63!
64! 181 2008-07-30 07:07:47Z raasch
65! bugfix: zero assignments to tendency arrays in case of restarts,
66! further extensions and modifications in the initialisation of the plant
67! canopy model,
68! allocation of hom_sum moved to parin, initialization of spectrum_x|y directly
69! after allocating theses arrays,
70! read data for recycling added as new initialization option,
71! dummy allocation for diss
72!
73! 138 2007-11-28 10:03:58Z letzel
74! New counter ngp_2dh_s_inner.
75! Allow new case bc_uv_t = 'dirichlet_0' for channel flow.
76! Corrected calculation of initial volume flow for 'set_1d-model_profiles' and
77! 'set_constant_profiles' in case of buildings in the reference cross-sections.
78!
79! 108 2007-08-24 15:10:38Z letzel
80! Flux initialization in case of coupled runs, +momentum fluxes at top boundary,
81! +arrays for phase speed c_u, c_v, c_w, indices for u|v|w_m_l|r changed
82! +qswst_remote in case of atmosphere model with humidity coupled to ocean
83! Rayleigh damping for ocean, optionally calculate km and kh from initial
84! TKE e_init
85!
86! 97 2007-06-21 08:23:15Z raasch
87! Initialization of salinity, call of init_ocean
88!
89! 87 2007-05-22 15:46:47Z raasch
90! var_hom and var_sum renamed pr_palm
91!
92! 75 2007-03-22 09:54:05Z raasch
93! Arrays for radiation boundary conditions are allocated (u_m_l, u_m_r, etc.),
94! bugfix for cases with the outflow damping layer extending over more than one
95! subdomain, moisture renamed humidity,
96! new initializing action "by_user" calls user_init_3d_model,
97! precipitation_amount/rate, ts_value are allocated, +module netcdf_control,
98! initial velocities at nzb+1 are regarded for volume
99! flow control in case they have been set zero before (to avoid small timesteps)
100! -uvmean_outflow, uxrp, vynp eliminated
101!
102! 19 2007-02-23 04:53:48Z raasch
103! +handling of top fluxes
104!
105! RCS Log replace by Id keyword, revision history cleaned up
106!
107! Revision 1.49  2006/08/22 15:59:07  raasch
108! No optimization of this file on the ibmy (Yonsei Univ.)
109!
110! Revision 1.1  1998/03/09 16:22:22  raasch
111! Initial revision
112!
113!
114! Description:
115! ------------
116! Allocation of arrays and initialization of the 3D model via
117! a) pre-run the 1D model
118! or
119! b) pre-set constant linear profiles
120! or
121! c) read values of a previous run
122!------------------------------------------------------------------------------!
123
124    USE advec_ws
125    USE arrays_3d
126    USE averaging
127    USE cloud_parameters
128    USE constants
129    USE control_parameters
130    USE cpulog
131    USE indices
132    USE interfaces
133    USE model_1d
134    USE netcdf_control
135    USE particle_attributes
136    USE pegrid
137    USE profil_parameter
138    USE random_function_mod
139    USE statistics
140
141    IMPLICIT NONE
142
143    INTEGER ::  i, ind_array(1), j, k, sr
144
145    INTEGER, DIMENSION(:), ALLOCATABLE ::  ngp_2dh_l
146
147    INTEGER, DIMENSION(:,:), ALLOCATABLE ::  ngp_2dh_outer_l,  &
148         ngp_2dh_s_inner_l
149
150    REAL ::  a, b
151
152    REAL, DIMENSION(1:2) ::  volume_flow_area_l, volume_flow_initial_l
153
154    REAL, DIMENSION(:), ALLOCATABLE ::  ngp_3d_inner_l, ngp_3d_inner_tmp
155
156
157!
158!-- Allocate arrays
159    ALLOCATE( ngp_2dh(0:statistic_regions), ngp_2dh_l(0:statistic_regions), &
160              ngp_3d(0:statistic_regions),                                  &
161              ngp_3d_inner(0:statistic_regions),                            &
162              ngp_3d_inner_l(0:statistic_regions),                          &
163              ngp_3d_inner_tmp(0:statistic_regions),                        &
164              sums_divnew_l(0:statistic_regions),                           &
165              sums_divold_l(0:statistic_regions) )
166    ALLOCATE( dp_smooth_factor(nzb:nzt), rdf(nzb+1:nzt) )
167    ALLOCATE( ngp_2dh_outer(nzb:nzt+1,0:statistic_regions),                 &
168              ngp_2dh_outer_l(nzb:nzt+1,0:statistic_regions),               &
169              ngp_2dh_s_inner(nzb:nzt+1,0:statistic_regions),               &
170              ngp_2dh_s_inner_l(nzb:nzt+1,0:statistic_regions),             &
171              rmask(nysg:nyng,nxlg:nxrg,0:statistic_regions),           &
172              sums(nzb:nzt+1,pr_palm+max_pr_user),                          &
173              sums_l(nzb:nzt+1,pr_palm+max_pr_user,0:threads_per_task-1),   &
174              sums_l_l(nzb:nzt+1,0:statistic_regions,0:threads_per_task-1), &
175              sums_up_fraction_l(10,3,0:statistic_regions),                 &
176              sums_wsts_bc_l(nzb:nzt+1,0:statistic_regions),                &
177              ts_value(dots_max,0:statistic_regions) )
178    ALLOCATE( km_damp_x(nxlg:nxrg), km_damp_y(nysg:nyng) )
179
180    ALLOCATE( rif_1(nysg:nyng,nxlg:nxrg), shf_1(nysg:nyng,nxlg:nxrg), &
181              ts(nysg:nyng,nxlg:nxrg), tswst_1(nysg:nyng,nxlg:nxrg),  &
182              us(nysg:nyng,nxlg:nxrg), usws_1(nysg:nyng,nxlg:nxrg),   &
183              uswst_1(nysg:nyng,nxlg:nxrg),                               &
184              vsws_1(nysg:nyng,nxlg:nxrg),                                &
185              vswst_1(nysg:nyng,nxlg:nxrg), z0(nysg:nyng,nxlg:nxrg) )
186
187    IF ( timestep_scheme(1:5) /= 'runge' )  THEN
188!
189!--    Leapfrog scheme needs two timelevels of diffusion quantities
190       ALLOCATE( rif_2(nysg:nyng,nxlg:nxrg),   &
191                 shf_2(nysg:nyng,nxlg:nxrg),   &
192                 tswst_2(nysg:nyng,nxlg:nxrg), &
193                 usws_2(nysg:nyng,nxlg:nxrg),  &
194                 uswst_2(nysg:nyng,nxlg:nxrg), &
195                 vswst_2(nysg:nyng,nxlg:nxrg), &
196                 vsws_2(nysg:nyng,nxlg:nxrg) )
197    ENDIF
198
199    ALLOCATE( d(nzb+1:nzta,nys:nyna,nxl:nxra),         &
200              e_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),  &
201              e_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),  &
202              e_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg),  &
203              kh_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
204              km_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
205              p(nzb:nzt+1,nysg:nyng,nxlg:nxrg),    &
206              pt_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
207              pt_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
208              pt_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
209              tend(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
210              u_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),  &
211              u_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),  &
212              u_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg),  &
213              v_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
214              v_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
215              v_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
216              w_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),  &
217              w_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),  &
218              w_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
219!
220!-- Following array is required for perturbation pressure within the iterative
221!-- pressure solvers. For the multistep schemes (Runge-Kutta), array p holds
222!-- the weighted average of the substeps and cannot be used in the Poisson
223!-- solver.
224    IF ( psolver == 'sor' )  THEN
225       ALLOCATE( p_loc(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
226    ELSEIF ( psolver == 'multigrid' )  THEN
227!
228!--    For performance reasons, multigrid is using one ghost layer only
229       ALLOCATE( p_loc(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
230    ENDIF
231
232    IF ( timestep_scheme(1:5) /= 'runge' )  THEN
233       ALLOCATE( kh_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
234                 km_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
235    ENDIF
236
237    IF ( humidity  .OR.  passive_scalar ) THEN
238!
239!--    2D-humidity/scalar arrays
240       ALLOCATE ( qs(nysg:nyng,nxlg:nxrg),     &
241                  qsws_1(nysg:nyng,nxlg:nxrg), &
242                  qswst_1(nysg:nyng,nxlg:nxrg) )
243
244       IF ( timestep_scheme(1:5) /= 'runge' )  THEN
245          ALLOCATE( qsws_2(nysg:nyng,nxlg:nxrg), &
246                    qswst_2(nysg:nyng,nxlg:nxrg) )
247       ENDIF
248!
249!--    3D-humidity/scalar arrays
250       ALLOCATE( q_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
251                 q_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
252                 q_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
253
254!
255!--    3D-arrays needed for humidity only
256       IF ( humidity )  THEN
257          ALLOCATE( vpt_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
258
259          IF ( timestep_scheme(1:5) /= 'runge' )  THEN
260             ALLOCATE( vpt_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
261          ENDIF
262
263          IF ( cloud_physics ) THEN
264!
265!--          Liquid water content
266             ALLOCATE ( ql_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
267!
268!--          Precipitation amount and rate (only needed if output is switched)
269             ALLOCATE( precipitation_amount(nysg:nyng,nxlg:nxrg), &
270                       precipitation_rate(nysg:nyng,nxlg:nxrg) )
271          ENDIF
272
273          IF ( cloud_droplets )  THEN
274!
275!--          Liquid water content, change in liquid water content,
276!--          real volume of particles (with weighting), volume of particles
277             ALLOCATE ( ql_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
278                        ql_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
279                        ql_v(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
280                        ql_vp(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
281          ENDIF
282
283       ENDIF
284
285    ENDIF
286
287    IF ( ocean )  THEN
288       ALLOCATE( saswsb_1(nysg:nyng,nxlg:nxrg), &
289                 saswst_1(nysg:nyng,nxlg:nxrg) )
290       ALLOCATE( prho_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
291                 rho_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),  &
292                 sa_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),   &
293                 sa_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),   &
294                 sa_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
295       prho => prho_1
296       rho  => rho_1  ! routines calc_mean_profile and diffusion_e require
297                      ! density to be apointer
298       IF ( humidity_remote )  THEN
299          ALLOCATE( qswst_remote(nysg:nyng,nxlg:nxrg))
300          qswst_remote = 0.0
301       ENDIF
302    ENDIF
303
304!
305!-- 3D-array for storing the dissipation, needed for calculating the sgs
306!-- particle velocities
307    IF ( use_sgs_for_particles )  THEN
308       ALLOCATE ( diss(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
309    ELSE
310       ALLOCATE ( diss(2,2,2) )  ! required because diss is used as a
311                                 ! formal parameter
312    ENDIF
313
314    IF ( dt_dosp /= 9999999.9 )  THEN
315       ALLOCATE( spectrum_x( 1:nx/2, 1:10, 1:10 ), &
316                 spectrum_y( 1:ny/2, 1:10, 1:10 ) )
317       spectrum_x = 0.0
318       spectrum_y = 0.0
319    ENDIF
320
321!
322!-- 3D-arrays for the leaf area density and the canopy drag coefficient
323    IF ( plant_canopy ) THEN
324       ALLOCATE ( lad_s(nzb:nzt+1,nysg:nyng,nxlg:nxrg),  &
325                  lad_u(nzb:nzt+1,nysg:nyng,nxlg:nxrg),  &
326                  lad_v(nzb:nzt+1,nysg:nyng,nxlg:nxrg),  &
327                  lad_w(nzb:nzt+1,nysg:nyng,nxlg:nxrg),  &
328                  cdc(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
329
330       IF ( passive_scalar ) THEN
331          ALLOCATE ( sls(nzb:nzt+1,nysg:nyng,nxlg:nxrg),   &
332                     sec(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 
333       ENDIF
334
335       IF ( cthf /= 0.0 ) THEN
336          ALLOCATE ( lai(nzb:nzt+1,nysg:nyng,nxlg:nxrg),   &
337                     canopy_heat_flux(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
338       ENDIF
339
340    ENDIF
341
342!
343!-- 4D-array for storing the Rif-values at vertical walls
344    IF ( topography /= 'flat' )  THEN
345       ALLOCATE( rif_wall(nzb:nzt+1,nysg:nyng,nxlg:nxrg,1:4) )
346       rif_wall = 0.0
347    ENDIF
348
349!
350!-- Arrays to store velocity data from t-dt and the phase speeds which
351!-- are needed for radiation boundary conditions
352    IF ( outflow_l )  THEN
353       ALLOCATE( u_m_l(nzb:nzt+1,nysg:nyng,1:2), &
354                 v_m_l(nzb:nzt+1,nysg:nyng,0:1), &
355                 w_m_l(nzb:nzt+1,nysg:nyng,0:1) )
356    ENDIF
357    IF ( outflow_r )  THEN
358       ALLOCATE( u_m_r(nzb:nzt+1,nysg:nyng,nx-1:nx), &
359                 v_m_r(nzb:nzt+1,nysg:nyng,nx-1:nx), &
360                 w_m_r(nzb:nzt+1,nysg:nyng,nx-1:nx) )
361    ENDIF
362    IF ( outflow_l  .OR.  outflow_r )  THEN
363       ALLOCATE( c_u(nzb:nzt+1,nysg:nyng), c_v(nzb:nzt+1,nysg:nyng), &
364                 c_w(nzb:nzt+1,nysg:nyng) )
365    ENDIF
366    IF ( outflow_s )  THEN
367       ALLOCATE( u_m_s(nzb:nzt+1,0:1,nxlg:nxrg), &
368                 v_m_s(nzb:nzt+1,1:2,nxlg:nxrg), &
369                 w_m_s(nzb:nzt+1,0:1,nxlg:nxrg) )
370    ENDIF
371    IF ( outflow_n )  THEN
372       ALLOCATE( u_m_n(nzb:nzt+1,ny-1:ny,nxlg:nxrg), &
373                 v_m_n(nzb:nzt+1,ny-1:ny,nxlg:nxrg), &
374                 w_m_n(nzb:nzt+1,ny-1:ny,nxlg:nxrg) )
375    ENDIF
376    IF ( outflow_s  .OR.  outflow_n )  THEN
377       ALLOCATE( c_u(nzb:nzt+1,nxlg:nxrg), c_v(nzb:nzt+1,nxlg:nxrg), &
378                 c_w(nzb:nzt+1,nxlg:nxrg) )
379    ENDIF
380
381!
382!-- Initial assignment of the pointers
383    IF ( timestep_scheme(1:5) /= 'runge' )  THEN
384
385       rif_m   => rif_1;    rif   => rif_2
386       shf_m   => shf_1;    shf   => shf_2
387       tswst_m => tswst_1;  tswst => tswst_2
388       usws_m  => usws_1;   usws  => usws_2
389       uswst_m => uswst_1;  uswst => uswst_2
390       vsws_m  => vsws_1;   vsws  => vsws_2
391       vswst_m => vswst_1;  vswst => vswst_2
392       e_m  => e_1;   e  => e_2;   e_p  => e_3;   te_m  => e_3
393       kh_m => kh_1;  kh => kh_2
394       km_m => km_1;  km => km_2
395       pt_m => pt_1;  pt => pt_2;  pt_p => pt_3;  tpt_m => pt_3
396       u_m  => u_1;   u  => u_2;   u_p  => u_3;   tu_m  => u_3
397       v_m  => v_1;   v  => v_2;   v_p  => v_3;   tv_m  => v_3
398       w_m  => w_1;   w  => w_2;   w_p  => w_3;   tw_m  => w_3
399
400       IF ( humidity  .OR.  passive_scalar )  THEN
401          qsws_m  => qsws_1;   qsws  => qsws_2
402          qswst_m => qswst_1;  qswst => qswst_2
403          q_m    => q_1;     q    => q_2;     q_p => q_3;     tq_m => q_3
404          IF ( humidity )        vpt_m  => vpt_1;   vpt  => vpt_2
405          IF ( cloud_physics )   ql   => ql_1
406          IF ( cloud_droplets )  THEN
407             ql   => ql_1
408             ql_c => ql_2
409          ENDIF
410       ENDIF
411
412    ELSE
413
414       rif   => rif_1
415       shf   => shf_1
416       tswst => tswst_1
417       usws  => usws_1
418       uswst => uswst_1
419       vsws  => vsws_1
420       vswst => vswst_1
421       e     => e_1;   e_p  => e_2;   te_m  => e_3;   e_m  => e_3
422       kh    => kh_1
423       km    => km_1
424       pt    => pt_1;  pt_p => pt_2;  tpt_m => pt_3;  pt_m => pt_3
425       u     => u_1;   u_p  => u_2;   tu_m  => u_3;   u_m  => u_3
426       v     => v_1;   v_p  => v_2;   tv_m  => v_3;   v_m  => v_3
427       w     => w_1;   w_p  => w_2;   tw_m  => w_3;   w_m  => w_3
428
429       IF ( humidity  .OR.  passive_scalar )  THEN
430          qsws   => qsws_1
431          qswst  => qswst_1
432          q      => q_1;     q_p  => q_2;     tq_m  => q_3;    q_m => q_3
433          IF ( humidity )        vpt  => vpt_1
434          IF ( cloud_physics )   ql   => ql_1
435          IF ( cloud_droplets )  THEN
436             ql   => ql_1
437             ql_c => ql_2
438          ENDIF
439       ENDIF
440
441       IF ( ocean )  THEN
442          saswsb => saswsb_1
443          saswst => saswst_1
444          sa     => sa_1;    sa_p => sa_2;    tsa_m => sa_3
445       ENDIF
446
447    ENDIF
448   
449!
450!-- Allocate arrays containing the RK coefficient for right evaluation of
451!-- perturbation pressure and turbulent fluxes. At this point it is needed
452!-- for right pressure correction during initialization. Further below
453!-- the real values will be set.
454    ALLOCATE (weight_substep(1:intermediate_timestep_count_max),           &
455              weight_pres(1:intermediate_timestep_count_max) )
456    weight_substep = 1.
457    weight_pres = 1.
458    intermediate_timestep_count = 1 ! needed for simulated_time=0
459       
460!
461!-- Initialize model variables
462    IF ( TRIM( initializing_actions ) /= 'read_restart_data'  .AND.  &
463         TRIM( initializing_actions ) /= 'cyclic_fill' )  THEN
464!
465!--    First model run of a possible job queue.
466!--    Initial profiles of the variables must be computes.
467       IF ( INDEX( initializing_actions, 'set_1d-model_profiles' ) /= 0 )  THEN
468!
469!--       Use solutions of the 1D model as initial profiles,
470!--       start 1D model
471          CALL init_1d_model
472!
473!--       Transfer initial profiles to the arrays of the 3D model
474          DO  i = nxlg, nxrg
475             DO  j = nysg, nyng
476                e(:,j,i)  = e1d
477                kh(:,j,i) = kh1d
478                km(:,j,i) = km1d
479                pt(:,j,i) = pt_init
480                u(:,j,i)  = u1d
481                v(:,j,i)  = v1d
482             ENDDO
483          ENDDO
484
485          IF ( humidity  .OR.  passive_scalar )  THEN
486             DO  i = nxlg, nxrg
487                DO  j = nysg, nyng
488                   q(:,j,i) = q_init
489                ENDDO
490             ENDDO
491          ENDIF
492
493          IF ( .NOT. constant_diffusion )  THEN
494             DO  i = nxlg, nxrg
495                DO  j = nysg, nyng
496                   e(:,j,i)  = e1d
497                ENDDO
498             ENDDO
499!
500!--          Store initial profiles for output purposes etc.
501             hom(:,1,25,:) = SPREAD( l1d, 2, statistic_regions+1 )
502
503             IF ( prandtl_layer )  THEN
504                rif  = rif1d(nzb+1)
505                ts   = 0.0  ! could actually be computed more accurately in the
506                            ! 1D model. Update when opportunity arises.
507                us   = us1d
508                usws = usws1d
509                vsws = vsws1d
510             ELSE
511                ts   = 0.0  ! must be set, because used in
512                rif  = 0.0  ! flowste
513                us   = 0.0
514                usws = 0.0
515                vsws = 0.0
516             ENDIF
517
518          ELSE
519             e    = 0.0  ! must be set, because used in
520             rif  = 0.0  ! flowste
521             ts   = 0.0
522             us   = 0.0
523             usws = 0.0
524             vsws = 0.0
525          ENDIF
526          uswst = top_momentumflux_u
527          vswst = top_momentumflux_v
528
529!
530!--       In every case qs = 0.0 (see also pt)
531!--       This could actually be computed more accurately in the 1D model.
532!--       Update when opportunity arises!
533          IF ( humidity  .OR.  passive_scalar )  qs = 0.0
534
535!
536!--       inside buildings set velocities back to zero
537          IF ( topography /= 'flat' )  THEN
538             DO  i = nxl-1, nxr+1
539                DO  j = nys-1, nyn+1
540                   u(nzb:nzb_u_inner(j,i),j,i) = 0.0
541                   v(nzb:nzb_v_inner(j,i),j,i) = 0.0
542                ENDDO
543             ENDDO
544             
545!
546!--          WARNING: The extra boundary conditions set after running the
547!--          -------  1D model impose an error on the divergence one layer
548!--                   below the topography; need to correct later
549!--          ATTENTION: Provisional correction for Piacsek & Williams
550!--          ---------  advection scheme: keep u and v zero one layer below
551!--                     the topography.
552!
553!--           Following was removed, because mirror boundary condition are
554!--           replaced by dirichlet boundary conditions
555!
556!             IF ( ibc_uv_b == 0 )  THEN
557!!
558!!--             Satisfying the Dirichlet condition with an extra layer below
559!!--             the surface where the u and v component change their sign.
560!                DO  i = nxl-1, nxr+1
561!                   DO  j = nys-1, nyn+1
562!                      IF ( nzb_u_inner(j,i) == 0 ) u(0,j,i) = -u(1,j,i)
563!                      IF ( nzb_v_inner(j,i) == 0 ) v(0,j,i) = -v(1,j,i)
564!                   ENDDO
565!                ENDDO
566!
567!             ELSE
568             IF ( ibc_uv_b == 1 )  THEN
569!
570!--             Neumann condition
571                DO  i = nxl-1, nxr+1
572                   DO  j = nys-1, nyn+1
573                      IF ( nzb_u_inner(j,i) == 0 ) u(0,j,i) = u(1,j,i)
574                      IF ( nzb_v_inner(j,i) == 0 ) v(0,j,i) = v(1,j,i)
575                   ENDDO
576                ENDDO
577
578             ENDIF
579
580          ENDIF
581
582       ELSEIF ( INDEX(initializing_actions, 'set_constant_profiles') /= 0 ) &
583       THEN
584!
585!--       Use constructed initial profiles (velocity constant with height,
586!--       temperature profile with constant gradient)
587          DO  i = nxlg, nxrg
588             DO  j = nysg, nyng
589                pt(:,j,i) = pt_init
590                u(:,j,i)  = u_init
591                v(:,j,i)  = v_init
592             ENDDO
593          ENDDO
594
595!
596!--       Set initial horizontal velocities at the lowest computational grid
597!--       levels to zero in order to avoid too small time steps caused by the
598!--       diffusion limit in the initial phase of a run (at k=1, dz/2 occurs
599!--       in the limiting formula!). The original values are stored to be later
600!--       used for volume flow control.
601          DO  i = nxlg, nxrg
602             DO  j = nysg, nyng
603                u(nzb:nzb_u_inner(j,i)+1,j,i) = 0.0
604                v(nzb:nzb_v_inner(j,i)+1,j,i) = 0.0
605             ENDDO
606          ENDDO
607
608          IF ( humidity  .OR.  passive_scalar )  THEN
609             DO  i = nxlg, nxrg
610                DO  j = nysg, nyng
611                   q(:,j,i) = q_init
612                ENDDO
613             ENDDO
614          ENDIF
615
616          IF ( ocean )  THEN
617             DO  i = nxlg, nxrg
618                DO  j = nysg, nyng
619                   sa(:,j,i) = sa_init
620                ENDDO
621             ENDDO
622          ENDIF
623         
624          IF ( constant_diffusion )  THEN
625             km   = km_constant
626             kh   = km / prandtl_number
627             e    = 0.0
628          ELSEIF ( e_init > 0.0 )  THEN
629             DO  k = nzb+1, nzt
630                km(k,:,:) = 0.1 * l_grid(k) * SQRT( e_init )
631             ENDDO
632             km(nzb,:,:)   = km(nzb+1,:,:)
633             km(nzt+1,:,:) = km(nzt,:,:)
634             kh   = km / prandtl_number
635             e    = e_init
636          ELSE
637             IF ( .NOT. ocean )  THEN
638                kh   = 0.01   ! there must exist an initial diffusion, because
639                km   = 0.01   ! otherwise no TKE would be produced by the
640                              ! production terms, as long as not yet
641                              ! e = (u*/cm)**2 at k=nzb+1
642             ELSE
643                kh   = 0.00001
644                km   = 0.00001
645             ENDIF
646             e    = 0.0
647          ENDIF
648          rif   = 0.0
649          ts    = 0.0
650          us    = 0.0
651          usws  = 0.0
652          uswst = top_momentumflux_u
653          vsws  = 0.0
654          vswst = top_momentumflux_v
655          IF ( humidity  .OR.  passive_scalar ) qs = 0.0
656
657!
658!--       Compute initial temperature field and other constants used in case
659!--       of a sloping surface
660          IF ( sloping_surface )  CALL init_slope
661
662       ELSEIF ( INDEX(initializing_actions, 'by_user') /= 0 ) &
663       THEN
664!
665!--       Initialization will completely be done by the user
666          CALL user_init_3d_model
667
668       ENDIF
669!
670!--    Bottom boundary
671       IF ( ibc_uv_b == 0 .OR. ibc_uv_b == 2  )  THEN
672          u(nzb,:,:) = 0.0
673          v(nzb,:,:) = 0.0
674       ENDIF
675
676!
677!--    Apply channel flow boundary condition
678       IF ( TRIM( bc_uv_t ) == 'dirichlet_0' )  THEN
679
680          u(nzt+1,:,:) = 0.0
681          v(nzt+1,:,:) = 0.0
682
683!--       For the Dirichlet condition to be correctly applied at the top, set
684!--       ug and vg to zero there
685          ug(nzt+1)    = 0.0
686          vg(nzt+1)    = 0.0
687
688       ENDIF
689
690!
691!--    Calculate virtual potential temperature
692       IF ( humidity ) vpt = pt * ( 1.0 + 0.61 * q )
693
694!
695!--    Store initial profiles for output purposes etc.
696       hom(:,1,5,:) = SPREAD( u(:,nys,nxl), 2, statistic_regions+1 )
697       hom(:,1,6,:) = SPREAD( v(:,nys,nxl), 2, statistic_regions+1 )
698       IF ( ibc_uv_b == 0 .OR. ibc_uv_b == 2)  THEN
699          hom(nzb,1,5,:) = 0.0   
700          hom(nzb,1,6,:) = 0.0 
701       ENDIF
702       hom(:,1,7,:)  = SPREAD( pt(:,nys,nxl), 2, statistic_regions+1 )
703       hom(:,1,23,:) = SPREAD( km(:,nys,nxl), 2, statistic_regions+1 )
704       hom(:,1,24,:) = SPREAD( kh(:,nys,nxl), 2, statistic_regions+1 )
705
706       IF ( ocean )  THEN
707!
708!--       Store initial salinity profile
709          hom(:,1,26,:)  = SPREAD( sa(:,nys,nxl), 2, statistic_regions+1 )
710       ENDIF
711
712       IF ( humidity )  THEN
713!
714!--       Store initial profile of total water content, virtual potential
715!--       temperature
716          hom(:,1,26,:) = SPREAD(   q(:,nys,nxl), 2, statistic_regions+1 )
717          hom(:,1,29,:) = SPREAD( vpt(:,nys,nxl), 2, statistic_regions+1 )
718          IF ( cloud_physics  .OR.  cloud_droplets ) THEN
719!
720!--          Store initial profile of specific humidity and potential
721!--          temperature
722             hom(:,1,27,:) = SPREAD(  q(:,nys,nxl), 2, statistic_regions+1 )
723             hom(:,1,28,:) = SPREAD( pt(:,nys,nxl), 2, statistic_regions+1 )
724          ENDIF
725       ENDIF
726
727       IF ( passive_scalar )  THEN
728!
729!--       Store initial scalar profile
730          hom(:,1,26,:) = SPREAD(  q(:,nys,nxl), 2, statistic_regions+1 )
731       ENDIF
732
733!
734!--    Initialize fluxes at bottom surface
735       IF ( use_surface_fluxes )  THEN
736
737          IF ( constant_heatflux )  THEN
738!
739!--          Heat flux is prescribed
740             IF ( random_heatflux )  THEN
741                CALL disturb_heatflux
742             ELSE
743                shf = surface_heatflux
744!
745!--             Over topography surface_heatflux is replaced by wall_heatflux(0)
746                IF ( TRIM( topography ) /= 'flat' )  THEN
747                   DO  i = nxlg, nxrg
748                      DO  j = nysg, nyng
749                         IF ( nzb_s_inner(j,i) /= 0 )  THEN
750                            shf(j,i) = wall_heatflux(0)
751                         ENDIF
752                      ENDDO
753                   ENDDO
754                ENDIF
755             ENDIF
756             IF ( ASSOCIATED( shf_m ) )  shf_m = shf
757          ENDIF
758
759!
760!--       Determine the near-surface water flux
761          IF ( humidity  .OR.  passive_scalar )  THEN
762             IF ( constant_waterflux )  THEN
763                qsws   = surface_waterflux
764!
765!--             Over topography surface_waterflux is replaced by
766!--             wall_humidityflux(0)
767                IF ( TRIM( topography ) /= 'flat' )  THEN
768                   wall_qflux = wall_humidityflux
769                   DO  i = nxlg, nxrg
770                      DO  j = nysg, nyng
771                         IF ( nzb_s_inner(j,i) /= 0 )  THEN
772                            qsws(j,i) = wall_qflux(0)
773                         ENDIF
774                      ENDDO
775                   ENDDO
776                ENDIF
777                IF ( ASSOCIATED( qsws_m ) )  qsws_m = qsws
778             ENDIF
779          ENDIF
780
781       ENDIF
782
783!
784!--    Initialize fluxes at top surface
785!--    Currently, only the heatflux and salinity flux can be prescribed.
786!--    The latent flux is zero in this case!
787       IF ( use_top_fluxes )  THEN
788
789          IF ( constant_top_heatflux )  THEN
790!
791!--          Heat flux is prescribed
792             tswst = top_heatflux
793             IF ( ASSOCIATED( tswst_m ) )  tswst_m = tswst
794
795             IF ( humidity  .OR.  passive_scalar )  THEN
796                qswst = 0.0
797                IF ( ASSOCIATED( qswst_m ) )  qswst_m = qswst
798             ENDIF
799
800             IF ( ocean )  THEN
801                saswsb = bottom_salinityflux
802                saswst = top_salinityflux
803             ENDIF
804          ENDIF
805
806!
807!--       Initialization in case of a coupled model run
808          IF ( coupling_mode == 'ocean_to_atmosphere' )  THEN
809             tswst = 0.0
810             IF ( ASSOCIATED( tswst_m ) )  tswst_m = tswst
811          ENDIF
812
813       ENDIF
814
815!
816!--    Initialize Prandtl layer quantities
817       IF ( prandtl_layer )  THEN
818
819          z0 = roughness_length
820
821          IF ( .NOT. constant_heatflux )  THEN 
822!
823!--          Surface temperature is prescribed. Here the heat flux cannot be
824!--          simply estimated, because therefore rif, u* and theta* would have
825!--          to be computed by iteration. This is why the heat flux is assumed
826!--          to be zero before the first time step. It approaches its correct
827!--          value in the course of the first few time steps.
828             shf   = 0.0
829             IF ( ASSOCIATED( shf_m ) )  shf_m = 0.0
830          ENDIF
831
832          IF ( humidity  .OR.  passive_scalar )  THEN
833             IF ( .NOT. constant_waterflux )  THEN
834                qsws   = 0.0
835                IF ( ASSOCIATED( qsws_m ) )   qsws_m = 0.0
836             ENDIF
837          ENDIF
838
839       ENDIF
840
841
842!
843!--    For the moment, vertical velocity is zero
844       w = 0.0
845
846!
847!--    Initialize array sums (must be defined in first call of pres)
848       sums = 0.0
849
850!
851!--    In case of iterative solvers, p must get an initial value
852       IF ( psolver == 'multigrid'  .OR.  psolver == 'sor' )  p = 0.0
853
854!
855!--    Treating cloud physics, liquid water content and precipitation amount
856!--    are zero at beginning of the simulation
857       IF ( cloud_physics )  THEN
858          ql = 0.0
859          IF ( precipitation )  precipitation_amount = 0.0
860       ENDIF
861       
862!
863!--    Initialize quantities for special advections schemes
864       CALL init_advec
865
866!
867!--    Impose vortex with vertical axis on the initial velocity profile
868       IF ( INDEX( initializing_actions, 'initialize_vortex' ) /= 0 )  THEN
869          CALL init_rankine
870       ENDIF
871
872!
873!--    Impose temperature anomaly (advection test only)
874       IF ( INDEX( initializing_actions, 'initialize_ptanom' ) /= 0 )  THEN
875          CALL init_pt_anomaly
876       ENDIF
877
878!
879!--    If required, change the surface temperature at the start of the 3D run
880       IF ( pt_surface_initial_change /= 0.0 )  THEN
881          pt(nzb,:,:) = pt(nzb,:,:) + pt_surface_initial_change
882       ENDIF
883
884!
885!--    If required, change the surface humidity/scalar at the start of the 3D
886!--    run
887       IF ( ( humidity .OR. passive_scalar ) .AND. &
888            q_surface_initial_change /= 0.0 )  THEN
889          q(nzb,:,:) = q(nzb,:,:) + q_surface_initial_change
890       ENDIF
891
892!
893!--    Initialize the random number generator (from numerical recipes)
894       CALL random_function_ini
895
896!
897!--    Initialize old and new time levels.
898       IF ( timestep_scheme(1:5) /= 'runge' )  THEN
899          e_m = e; pt_m = pt; u_m = u; v_m = v; w_m = w; kh_m = kh; km_m = km
900       ELSE
901          te_m = 0.0; tpt_m = 0.0; tu_m = 0.0; tv_m = 0.0; tw_m = 0.0
902       ENDIF
903       e_p = e; pt_p = pt; u_p = u; v_p = v; w_p = w
904
905       IF ( humidity  .OR.  passive_scalar )  THEN
906          IF ( ASSOCIATED( q_m ) )  q_m = q
907          IF ( timestep_scheme(1:5) == 'runge' )  tq_m = 0.0
908          q_p = q
909          IF ( humidity  .AND.  ASSOCIATED( vpt_m ) )  vpt_m = vpt
910       ENDIF
911
912       IF ( ocean )  THEN
913          tsa_m = 0.0
914          sa_p  = sa
915       ENDIF
916       
917
918    ELSEIF ( TRIM( initializing_actions ) == 'read_restart_data'  .OR.    &
919         TRIM( initializing_actions ) == 'cyclic_fill' )  &
920    THEN
921!
922!--    When reading data for cyclic fill of 3D prerun data, first read
923!--    some of the global variables from restart file
924       IF ( TRIM( initializing_actions ) == 'cyclic_fill' )  THEN
925
926          CALL read_parts_of_var_list
927          CALL close_file( 13 )
928
929!
930!--       Initialization of the turbulence recycling method
931          IF ( turbulent_inflow )  THEN
932!
933!--          Store temporally and horizontally averaged vertical profiles to be
934!--          used as mean inflow profiles
935             ALLOCATE( mean_inflow_profiles(nzb:nzt+1,5) )
936
937             mean_inflow_profiles(:,1) = hom_sum(:,1,0)    ! u
938             mean_inflow_profiles(:,2) = hom_sum(:,2,0)    ! v
939             mean_inflow_profiles(:,4) = hom_sum(:,4,0)    ! pt
940             mean_inflow_profiles(:,5) = hom_sum(:,8,0)    ! e
941
942!
943!--          Use these mean profiles for the inflow (provided that Dirichlet
944!--          conditions are used)
945             IF ( inflow_l )  THEN
946                DO  j = nysg, nyng
947                   DO  k = nzb, nzt+1
948                      u(k,j,nxlg:-1)  = mean_inflow_profiles(k,1)
949                      v(k,j,nxlg:-1)  = mean_inflow_profiles(k,2)
950                      w(k,j,nxlg:-1)  = 0.0
951                      pt(k,j,nxlg:-1) = mean_inflow_profiles(k,4)
952                      e(k,j,nxlg:-1)  = mean_inflow_profiles(k,5)
953                   ENDDO
954                ENDDO
955             ENDIF
956
957!
958!--          Calculate the damping factors to be used at the inflow. For a
959!--          turbulent inflow the turbulent fluctuations have to be limited
960!--          vertically because otherwise the turbulent inflow layer will grow
961!--          in time.
962             IF ( inflow_damping_height == 9999999.9 )  THEN
963!
964!--             Default: use the inversion height calculated by the prerun; if
965!--             this is zero, inflow_damping_height must be explicitly
966!--             specified.
967                IF ( hom_sum(nzb+6,pr_palm,0) /= 0.0 )  THEN
968                   inflow_damping_height = hom_sum(nzb+6,pr_palm,0)
969                ELSE
970                   WRITE( message_string, * ) 'inflow_damping_height must be ',&
971                        'explicitly specified because&the inversion height ', &
972                        'calculated by the prerun is zero.'
973                   CALL message( 'init_3d_model', 'PA0318', 1, 2, 0, 6, 0 )
974                ENDIF
975
976             ENDIF
977
978             IF ( inflow_damping_width == 9999999.9 )  THEN
979!
980!--             Default for the transition range: one tenth of the undamped
981!--             layer
982                inflow_damping_width = 0.1 * inflow_damping_height
983
984             ENDIF
985
986             ALLOCATE( inflow_damping_factor(nzb:nzt+1) )
987
988             DO  k = nzb, nzt+1
989
990                IF ( zu(k) <= inflow_damping_height )  THEN
991                   inflow_damping_factor(k) = 1.0
992                ELSEIF ( zu(k) <= inflow_damping_height +  &
993                                  inflow_damping_width )  THEN
994                   inflow_damping_factor(k) = 1.0 -                            &
995                                           ( zu(k) - inflow_damping_height ) / &
996                                           inflow_damping_width
997                ELSE
998                   inflow_damping_factor(k) = 0.0
999                ENDIF
1000
1001             ENDDO
1002          ENDIF
1003
1004       ENDIF
1005
1006!
1007!--    Read binary data from restart file
1008
1009       CALL read_3d_binary
1010
1011!
1012!--    Inside buildings set velocities and TKE back to zero
1013       IF ( TRIM( initializing_actions ) == 'cyclic_fill' .AND.  &
1014            topography /= 'flat' )  THEN
1015!
1016!--       Inside buildings set velocities and TKE back to zero.
1017!--       Other scalars (pt, q, s, km, kh, p, sa, ...) are ignored at present,
1018!--       maybe revise later.
1019          IF ( timestep_scheme(1:5) == 'runge' )  THEN
1020             DO  i = nxlg, nxrg
1021                DO  j = nysg, nyng
1022                   u  (nzb:nzb_u_inner(j,i),j,i) = 0.0
1023                   v  (nzb:nzb_v_inner(j,i),j,i) = 0.0
1024                   w  (nzb:nzb_w_inner(j,i),j,i) = 0.0
1025                   e  (nzb:nzb_w_inner(j,i),j,i) = 0.0
1026                   u_m(nzb:nzb_u_inner(j,i),j,i) = 0.0
1027                   v_m(nzb:nzb_v_inner(j,i),j,i) = 0.0
1028                   w_m(nzb:nzb_w_inner(j,i),j,i) = 0.0
1029                   e_m(nzb:nzb_w_inner(j,i),j,i) = 0.0
1030                   tu_m(nzb:nzb_u_inner(j,i),j,i) = 0.0
1031                   tv_m(nzb:nzb_v_inner(j,i),j,i) = 0.0
1032                   tw_m(nzb:nzb_w_inner(j,i),j,i) = 0.0
1033                   te_m(nzb:nzb_w_inner(j,i),j,i) = 0.0
1034                   tpt_m(nzb:nzb_w_inner(j,i),j,i) = 0.0
1035                ENDDO
1036             ENDDO
1037          ELSE
1038             DO  i = nxlg, nxrg
1039                DO  j = nysg, nyng
1040                   u  (nzb:nzb_u_inner(j,i),j,i) = 0.0
1041                   v  (nzb:nzb_v_inner(j,i),j,i) = 0.0
1042                   w  (nzb:nzb_w_inner(j,i),j,i) = 0.0
1043                   e  (nzb:nzb_w_inner(j,i),j,i) = 0.0
1044                   u_m(nzb:nzb_u_inner(j,i),j,i) = 0.0
1045                   v_m(nzb:nzb_v_inner(j,i),j,i) = 0.0
1046                   w_m(nzb:nzb_w_inner(j,i),j,i) = 0.0
1047                   e_m(nzb:nzb_w_inner(j,i),j,i) = 0.0
1048                   u_p(nzb:nzb_u_inner(j,i),j,i) = 0.0
1049                   v_p(nzb:nzb_v_inner(j,i),j,i) = 0.0
1050                   w_p(nzb:nzb_w_inner(j,i),j,i) = 0.0
1051                   e_p(nzb:nzb_w_inner(j,i),j,i) = 0.0
1052                ENDDO
1053             ENDDO
1054          ENDIF
1055
1056       ENDIF
1057
1058!
1059!--    Calculate initial temperature field and other constants used in case
1060!--    of a sloping surface
1061       IF ( sloping_surface )  CALL init_slope
1062
1063!
1064!--    Initialize new time levels (only done in order to set boundary values
1065!--    including ghost points)
1066       e_p = e; pt_p = pt; u_p = u; v_p = v; w_p = w
1067       IF ( humidity  .OR.  passive_scalar )  q_p = q
1068       IF ( ocean )  sa_p = sa
1069
1070!
1071!--    Allthough tendency arrays are set in prognostic_equations, they have
1072!--    have to be predefined here because they are used (but multiplied with 0)
1073!--    there before they are set.
1074       IF ( timestep_scheme(1:5) == 'runge' )  THEN
1075          te_m = 0.0; tpt_m = 0.0; tu_m = 0.0; tv_m = 0.0; tw_m = 0.0
1076          IF ( humidity  .OR.  passive_scalar )  tq_m = 0.0
1077          IF ( ocean )  tsa_m = 0.0
1078       ENDIF
1079
1080    ELSE
1081!
1082!--    Actually this part of the programm should not be reached
1083       message_string = 'unknown initializing problem'
1084       CALL message( 'init_3d_model', 'PA0193', 1, 2, 0, 6, 0 )
1085    ENDIF
1086
1087
1088    IF (  TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
1089!
1090!--    Initialize old timelevels needed for radiation boundary conditions
1091       IF ( outflow_l )  THEN
1092          u_m_l(:,:,:) = u(:,:,1:2)
1093          v_m_l(:,:,:) = v(:,:,0:1)
1094          w_m_l(:,:,:) = w(:,:,0:1)
1095       ENDIF
1096       IF ( outflow_r )  THEN
1097          u_m_r(:,:,:) = u(:,:,nx-1:nx)
1098          v_m_r(:,:,:) = v(:,:,nx-1:nx)
1099          w_m_r(:,:,:) = w(:,:,nx-1:nx)
1100       ENDIF
1101       IF ( outflow_s )  THEN
1102          u_m_s(:,:,:) = u(:,0:1,:)
1103          v_m_s(:,:,:) = v(:,1:2,:)
1104          w_m_s(:,:,:) = w(:,0:1,:)
1105       ENDIF
1106       IF ( outflow_n )  THEN
1107          u_m_n(:,:,:) = u(:,ny-1:ny,:)
1108          v_m_n(:,:,:) = v(:,ny-1:ny,:)
1109          w_m_n(:,:,:) = w(:,ny-1:ny,:)
1110       ENDIF
1111       
1112    ENDIF
1113
1114!
1115!-- Calculate the initial volume flow at the right and north boundary
1116    IF ( conserve_volume_flow ) THEN
1117
1118       volume_flow_initial_l = 0.0
1119       volume_flow_area_l    = 0.0
1120 
1121       IF  ( TRIM( initializing_actions ) == 'cyclic_fill' )  THEN
1122
1123          IF ( nxr == nx )  THEN
1124             DO  j = nys, nyn
1125                DO  k = nzb_2d(j,nx) + 1, nzt
1126                   volume_flow_initial_l(1) = volume_flow_initial_l(1) + &
1127                                              hom_sum(k,1,0) * dzw(k)
1128                   volume_flow_area_l(1)    = volume_flow_area_l(1) + dzw(k)
1129                ENDDO
1130             ENDDO
1131          ENDIF
1132         
1133          IF ( nyn == ny )  THEN
1134             DO  i = nxl, nxr
1135                DO  k = nzb_2d(ny,i) + 1, nzt 
1136                   volume_flow_initial_l(2) = volume_flow_initial_l(2) + &
1137                                               hom_sum(k,2,0) * dzw(k)
1138                   volume_flow_area_l(2)    = volume_flow_area_l(2) + dzw(k)
1139                ENDDO
1140             ENDDO
1141          ENDIF
1142
1143       ELSEIF ( TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
1144
1145          IF ( nxr == nx )  THEN
1146             DO  j = nys, nyn
1147                DO  k = nzb_2d(j,nx) + 1, nzt
1148                   volume_flow_initial_l(1) = volume_flow_initial_l(1) + &
1149                                               u(k,j,nx) * dzw(k)
1150                   volume_flow_area_l(1)    = volume_flow_area_l(1) + dzw(k)
1151                ENDDO
1152             ENDDO
1153          ENDIF
1154         
1155          IF ( nyn == ny )  THEN
1156             DO  i = nxl, nxr
1157                DO  k = nzb_2d(ny,i) + 1, nzt 
1158                   volume_flow_initial_l(2) = volume_flow_initial_l(2) + &
1159                                              v(k,ny,i) * dzw(k)
1160                   volume_flow_area_l(2)    = volume_flow_area_l(2) + dzw(k)
1161                ENDDO
1162             ENDDO
1163          ENDIF
1164
1165       ENDIF
1166
1167#if defined( __parallel )
1168       CALL MPI_ALLREDUCE( volume_flow_initial_l(1), volume_flow_initial(1),&
1169                           2, MPI_REAL, MPI_SUM, comm2d, ierr )
1170       CALL MPI_ALLREDUCE( volume_flow_area_l(1), volume_flow_area(1),      &
1171                           2, MPI_REAL, MPI_SUM, comm2d, ierr )
1172
1173#else
1174       volume_flow_initial = volume_flow_initial_l
1175       volume_flow_area    = volume_flow_area_l
1176#endif 
1177
1178!
1179!--       In case of 'bulk_velocity' mode, volume_flow_initial is overridden
1180!--       and calculated from u|v_bulk instead.
1181       IF ( TRIM( conserve_volume_flow_mode ) == 'bulk_velocity' )  THEN
1182          volume_flow_initial(1) = u_bulk * volume_flow_area(1)
1183          volume_flow_initial(2) = v_bulk * volume_flow_area(2)
1184       ENDIF
1185
1186    ENDIF
1187
1188
1189!
1190!-- Impose random perturbation on the horizontal velocity field and then
1191!-- remove the divergences from the velocity field at the initial stage
1192    IF ( create_disturbances .AND. &
1193         TRIM( initializing_actions ) /= 'read_restart_data'  .AND.  &
1194         TRIM( initializing_actions ) /= 'cyclic_fill' )  THEN
1195
1196       CALL disturb_field( nzb_u_inner, tend, u )
1197       CALL disturb_field( nzb_v_inner, tend, v )
1198       n_sor = nsor_ini
1199       CALL pres
1200       n_sor = nsor
1201    ENDIF
1202
1203!
1204!-- Initialization of the leaf area density
1205     IF ( plant_canopy ) THEN
1206 
1207       SELECT CASE ( TRIM( canopy_mode ) )
1208
1209          CASE( 'block' )
1210
1211             DO  i = nxlg, nxrg
1212                DO  j = nysg, nyng
1213                   lad_s(:,j,i) = lad(:)
1214                   cdc(:,j,i)   = drag_coefficient
1215                   IF ( passive_scalar ) THEN
1216                      sls(:,j,i) = leaf_surface_concentration
1217                      sec(:,j,i) = scalar_exchange_coefficient
1218                   ENDIF
1219                ENDDO
1220             ENDDO
1221
1222          CASE DEFAULT
1223
1224!
1225!--          The DEFAULT case is reached either if the parameter
1226!--          canopy mode contains a wrong character string or if the
1227!--          user has coded a special case in the user interface.
1228!--          There, the subroutine user_init_plant_canopy checks
1229!--          which of these two conditions applies.
1230             CALL user_init_plant_canopy
1231 
1232          END SELECT
1233
1234       CALL exchange_horiz( lad_s, nbgp )
1235       CALL exchange_horiz( cdc, nbgp )
1236
1237       IF ( passive_scalar ) THEN
1238          CALL exchange_horiz( sls, nbgp )
1239          CALL exchange_horiz( sec, nbgp )
1240       ENDIF
1241
1242!
1243!--    Sharp boundaries of the plant canopy in horizontal directions
1244!--    In vertical direction the interpolation is retained, as the leaf
1245!--    area density is initialised by prescribing a vertical profile
1246!--    consisting of piecewise linear segments. The upper boundary
1247!--    of the plant canopy is now defined by lad_w(pch_index,:,:) = 0.0.
1248
1249       DO  i = nxl, nxr
1250          DO  j = nys, nyn
1251             DO  k = nzb, nzt+1 
1252                IF ( lad_s(k,j,i) > 0.0 ) THEN
1253                   lad_u(k,j,i)   = lad_s(k,j,i) 
1254                   lad_u(k,j,i+1) = lad_s(k,j,i)
1255                   lad_v(k,j,i)   = lad_s(k,j,i)
1256                   lad_v(k,j+1,i) = lad_s(k,j,i)
1257                ENDIF
1258             ENDDO
1259             DO  k = nzb, nzt
1260                lad_w(k,j,i) = 0.5 * ( lad_s(k+1,j,i) + lad_s(k,j,i) )
1261             ENDDO
1262          ENDDO
1263       ENDDO
1264
1265       lad_w(pch_index,:,:) = 0.0
1266       lad_w(nzt+1,:,:)     = lad_w(nzt,:,:)
1267
1268       CALL exchange_horiz( lad_u, nbgp )
1269       CALL exchange_horiz( lad_v, nbgp )
1270       CALL exchange_horiz( lad_w, nbgp )
1271
1272!
1273!--    Initialisation of the canopy heat source distribution
1274       IF ( cthf /= 0.0 ) THEN
1275!
1276!--       Piecewise evaluation of the leaf area index by
1277!--       integration of the leaf area density
1278          lai(:,:,:) = 0.0
1279          DO  i = nxlg, nxrg
1280             DO  j = nysg, nyng
1281                DO  k = pch_index-1, 0, -1
1282                   lai(k,j,i) = lai(k+1,j,i) +                   &
1283                                ( 0.5 * ( lad_w(k+1,j,i) +       &
1284                                          lad_s(k+1,j,i) ) *     &
1285                                  ( zw(k+1) - zu(k+1) ) )  +     &
1286                                ( 0.5 * ( lad_w(k,j,i)   +       &
1287                                          lad_s(k+1,j,i) ) *     &
1288                                  ( zu(k+1) - zw(k) ) )
1289                ENDDO
1290             ENDDO
1291          ENDDO
1292
1293!
1294!--       Evaluation of the upward kinematic vertical heat flux within the
1295!--       canopy
1296          DO  i = nxlg, nxrg
1297             DO  j = nysg, nyng
1298                DO  k = 0, pch_index
1299                   canopy_heat_flux(k,j,i) = cthf *                    &
1300                                             exp( -0.6 * lai(k,j,i) )
1301                ENDDO
1302             ENDDO
1303          ENDDO
1304
1305!
1306!--       The near surface heat flux is derived from the heat flux
1307!--       distribution within the canopy
1308          shf(:,:) = canopy_heat_flux(0,:,:)
1309
1310          IF ( ASSOCIATED( shf_m ) ) shf_m = shf
1311
1312       ENDIF
1313
1314    ENDIF
1315
1316!
1317!-- If required, initialize dvrp-software
1318    IF ( dt_dvrp /= 9999999.9 )  CALL init_dvrp
1319
1320    IF ( ocean )  THEN
1321!
1322!--    Initialize quantities needed for the ocean model
1323       CALL init_ocean
1324
1325    ELSE
1326!
1327!--    Initialize quantities for handling cloud physics
1328!--    This routine must be called before init_particles, because
1329!--    otherwise, array pt_d_t, needed in data_output_dvrp (called by
1330!--    init_particles) is not defined.
1331       CALL init_cloud_physics
1332    ENDIF
1333
1334!
1335!-- If required, initialize particles
1336    IF ( particle_advection )  CALL init_particles
1337
1338!
1339!-- Initialize the ws-scheme.   
1340    IF ( ws_scheme_sca .OR. ws_scheme_mom )  CALL ws_init       
1341
1342!
1343!-- Setting weighting factors for right evaluation of perturbation pressure
1344!-- and turbulent quantities during the RK substeps.               
1345    IF ( TRIM(timestep_scheme) == 'runge-kutta-3' )  THEN    ! RK3       
1346       weight_substep(1) = 0.166666666666666
1347       weight_substep(2) = 0.3
1348       weight_substep(3) = 0.533333333333333
1349         
1350       weight_pres(1) = 0.333333333333333
1351       weight_pres(2) = 0.416666666666666
1352       weight_pres(3) = 0.25         
1353    ELSEIF ( TRIM(timestep_scheme) == 'runge-kutta-2' )  THEN  ! RK2       
1354       weight_substep(1) = 0.5
1355       weight_substep(2) = 0.5
1356         
1357       weight_pres(1) = 0.5
1358       weight_pres(2) = 0.5         
1359    ELSE                                                  ! Euler and Leapfrog       
1360       weight_substep(1) = 1.0     
1361       weight_pres(1) = 1.0                   
1362    ENDIF
1363
1364!
1365!-- Initialize Rayleigh damping factors
1366    rdf = 0.0
1367    IF ( rayleigh_damping_factor /= 0.0 )  THEN
1368       IF ( .NOT. ocean )  THEN
1369          DO  k = nzb+1, nzt
1370             IF ( zu(k) >= rayleigh_damping_height )  THEN
1371                rdf(k) = rayleigh_damping_factor * &
1372                      ( SIN( pi * 0.5 * ( zu(k) - rayleigh_damping_height )    &
1373                                      / ( zu(nzt) - rayleigh_damping_height ) )&
1374                      )**2
1375             ENDIF
1376          ENDDO
1377       ELSE
1378          DO  k = nzt, nzb+1, -1
1379             IF ( zu(k) <= rayleigh_damping_height )  THEN
1380                rdf(k) = rayleigh_damping_factor * &
1381                      ( SIN( pi * 0.5 * ( rayleigh_damping_height - zu(k) )    &
1382                                      / ( rayleigh_damping_height - zu(nzb+1)))&
1383                      )**2
1384             ENDIF
1385          ENDDO
1386       ENDIF
1387    ENDIF
1388
1389!
1390!-- Initialize the starting level and the vertical smoothing factor used for
1391!-- the external pressure gradient
1392    dp_smooth_factor = 1.0
1393    IF ( dp_external )  THEN
1394!
1395!--    Set the starting level dp_level_ind_b only if it has not been set before
1396!--    (e.g. in init_grid).
1397       IF ( dp_level_ind_b == 0 )  THEN
1398          ind_array = MINLOC( ABS( dp_level_b - zu ) )
1399          dp_level_ind_b = ind_array(1) - 1 + nzb
1400                                        ! MINLOC uses lower array bound 1
1401       ENDIF
1402       IF ( dp_smooth )  THEN
1403          dp_smooth_factor(:dp_level_ind_b) = 0.0
1404          DO  k = dp_level_ind_b+1, nzt
1405             dp_smooth_factor(k) = 0.5 * ( 1.0 + SIN( pi * &
1406                  ( REAL( k - dp_level_ind_b ) /  &
1407                    REAL( nzt - dp_level_ind_b ) - 0.5 ) ) )
1408          ENDDO
1409       ENDIF
1410    ENDIF
1411
1412!
1413!-- Initialize diffusivities used within the outflow damping layer in case of
1414!-- non-cyclic lateral boundaries. A linear increase is assumed over the first
1415!-- half of the width of the damping layer
1416    IF ( bc_lr_dirrad )  THEN
1417
1418       DO  i = nxlg, nxrg
1419          IF ( i >= nx - outflow_damping_width )  THEN
1420             km_damp_x(i) = km_damp_max * MIN( 1.0,                    &
1421                            ( i - ( nx - outflow_damping_width ) ) /   &
1422                            REAL( outflow_damping_width/2 )            &
1423                                             )
1424          ELSE
1425             km_damp_x(i) = 0.0
1426          ENDIF
1427       ENDDO
1428
1429    ELSEIF ( bc_lr_raddir )  THEN
1430
1431       DO  i = nxlg, nxrg
1432          IF ( i <= outflow_damping_width )  THEN
1433             km_damp_x(i) = km_damp_max * MIN( 1.0,                    &
1434                            ( outflow_damping_width - i ) /            &
1435                            REAL( outflow_damping_width/2 )            &
1436                                             )
1437          ELSE
1438             km_damp_x(i) = 0.0
1439          ENDIF
1440       ENDDO
1441
1442    ENDIF
1443
1444    IF ( bc_ns_dirrad )  THEN
1445
1446       DO  j = nysg, nyng
1447          IF ( j >= ny - outflow_damping_width )  THEN
1448             km_damp_y(j) = km_damp_max * MIN( 1.0,                    &
1449                            ( j - ( ny - outflow_damping_width ) ) /   &
1450                            REAL( outflow_damping_width/2 )            &
1451                                             )
1452          ELSE
1453             km_damp_y(j) = 0.0
1454          ENDIF
1455       ENDDO
1456
1457    ELSEIF ( bc_ns_raddir )  THEN
1458
1459       DO  j = nysg, nyng
1460          IF ( j <= outflow_damping_width )  THEN
1461             km_damp_y(j) = km_damp_max * MIN( 1.0,                    &
1462                            ( outflow_damping_width - j ) /            &
1463                            REAL( outflow_damping_width/2 )            &
1464                                             )
1465          ELSE
1466             km_damp_y(j) = 0.0
1467          ENDIF
1468       ENDDO
1469
1470    ENDIF
1471
1472!
1473!-- Initialize local summation arrays for UP flow_statistics. This is necessary
1474!-- because they may not yet have been initialized when they are called from
1475!-- flow_statistics (or - depending on the chosen model run - are never
1476!-- initialized)
1477    sums_divnew_l      = 0.0
1478    sums_divold_l      = 0.0
1479    sums_l_l           = 0.0
1480    sums_up_fraction_l = 0.0
1481    sums_wsts_bc_l     = 0.0
1482
1483!
1484!-- Pre-set masks for regional statistics. Default is the total model domain.
1485    rmask = 1.0
1486
1487!
1488!-- User-defined initializing actions. Check afterwards, if maximum number
1489!-- of allowed timeseries is not exceeded
1490    CALL user_init
1491
1492    IF ( dots_num > dots_max )  THEN
1493       WRITE( message_string, * ) 'number of time series quantities exceeds', &
1494                                  ' its maximum of dots_max = ', dots_max,    &
1495                                  ' &Please increase dots_max in modules.f90.'
1496       CALL message( 'init_3d_model', 'PA0194', 1, 2, 0, 6, 0 )   
1497    ENDIF
1498
1499!
1500!-- Input binary data file is not needed anymore. This line must be placed
1501!-- after call of user_init!
1502    CALL close_file( 13 )
1503
1504!
1505!-- Compute total sum of active mask grid points
1506!-- ngp_2dh: number of grid points of a horizontal cross section through the
1507!--          total domain
1508!-- ngp_3d:  number of grid points of the total domain
1509    ngp_2dh_outer_l   = 0
1510    ngp_2dh_outer     = 0
1511    ngp_2dh_s_inner_l = 0
1512    ngp_2dh_s_inner   = 0
1513    ngp_2dh_l         = 0
1514    ngp_2dh           = 0
1515    ngp_3d_inner_l    = 0.0
1516    ngp_3d_inner      = 0
1517    ngp_3d            = 0
1518    ngp_sums          = ( nz + 2 ) * ( pr_palm + max_pr_user )
1519
1520    DO  sr = 0, statistic_regions
1521       DO  i = nxl, nxr
1522          DO  j = nys, nyn
1523             IF ( rmask(j,i,sr) == 1.0 )  THEN
1524!
1525!--             All xy-grid points
1526                ngp_2dh_l(sr) = ngp_2dh_l(sr) + 1
1527!
1528!--             xy-grid points above topography
1529                DO  k = nzb_s_outer(j,i), nz + 1
1530                   ngp_2dh_outer_l(k,sr) = ngp_2dh_outer_l(k,sr) + 1
1531                ENDDO
1532                DO  k = nzb_s_inner(j,i), nz + 1
1533                   ngp_2dh_s_inner_l(k,sr) = ngp_2dh_s_inner_l(k,sr) + 1
1534                ENDDO
1535!
1536!--             All grid points of the total domain above topography
1537                ngp_3d_inner_l(sr) = ngp_3d_inner_l(sr) + &
1538                                     ( nz - nzb_s_inner(j,i) + 2 )
1539             ENDIF
1540          ENDDO
1541       ENDDO
1542    ENDDO
1543
1544    sr = statistic_regions + 1
1545#if defined( __parallel )
1546    IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
1547    CALL MPI_ALLREDUCE( ngp_2dh_l(0), ngp_2dh(0), sr, MPI_INTEGER, MPI_SUM,   &
1548                        comm2d, ierr )
1549    IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
1550    CALL MPI_ALLREDUCE( ngp_2dh_outer_l(0,0), ngp_2dh_outer(0,0), (nz+2)*sr,  &
1551                        MPI_INTEGER, MPI_SUM, comm2d, ierr )
1552    IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
1553    CALL MPI_ALLREDUCE( ngp_2dh_s_inner_l(0,0), ngp_2dh_s_inner(0,0),         &
1554                        (nz+2)*sr, MPI_INTEGER, MPI_SUM, comm2d, ierr )
1555    IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
1556    CALL MPI_ALLREDUCE( ngp_3d_inner_l(0), ngp_3d_inner_tmp(0), sr, MPI_REAL, &
1557                        MPI_SUM, comm2d, ierr )
1558    ngp_3d_inner = INT( ngp_3d_inner_tmp, KIND = SELECTED_INT_KIND( 18 ) )
1559#else
1560    ngp_2dh         = ngp_2dh_l
1561    ngp_2dh_outer   = ngp_2dh_outer_l
1562    ngp_2dh_s_inner = ngp_2dh_s_inner_l
1563    ngp_3d_inner    = INT( ngp_3d_inner_l, KIND = SELECTED_INT_KIND( 18 ) )
1564#endif
1565
1566    ngp_3d = INT ( ngp_2dh, KIND = SELECTED_INT_KIND( 18 ) ) * &
1567             INT ( (nz + 2 ), KIND = SELECTED_INT_KIND( 18 ) )
1568
1569!
1570!-- Set a lower limit of 1 in order to avoid zero divisions in flow_statistics,
1571!-- buoyancy, etc. A zero value will occur for cases where all grid points of
1572!-- the respective subdomain lie below the surface topography
1573    ngp_2dh_outer   = MAX( 1, ngp_2dh_outer(:,:)   ) 
1574    ngp_3d_inner    = MAX( INT(1, KIND = SELECTED_INT_KIND( 18 )),            &
1575                           ngp_3d_inner(:) )
1576    ngp_2dh_s_inner = MAX( 1, ngp_2dh_s_inner(:,:) ) 
1577
1578    DEALLOCATE( ngp_2dh_l, ngp_2dh_outer_l, ngp_3d_inner_l, ngp_3d_inner_tmp )
1579
1580
1581 END SUBROUTINE init_3d_model
Note: See TracBrowser for help on using the repository browser.