source: palm/tags/release-3.7a/SOURCE/init_3d_model.f90 @ 3986

Last change on this file since 3986 was 482, checked in by raasch, 14 years ago

New:
---
compare_palm_logs is additionally compiled with mbuild -u (Makefile in trunk/UTIL)

make options (mopts) to be set by configuration file implemented (mrun, mbuild)

humidity=.T. is now usable for runs with topography. wall_humidityflux and
wall_scalarflux are the corresponding new parin arrays.
(check_parameters, init_3d_model, parin)

Large scale vertical motion (subsidence/ascent) can be applied to the
prognostic equation for the potential temperature. (check_parameters, header,
Makefile, modules, parin, prognostic_equations, read_var_list, subsidence,
write_var_list)

A simple method for installing and running palm (with limited features)
has been added. (Makefile, palm_simple_install, palm_simple_run)

Changed:


2d-decomposition is default for Cray-XT machines (init_pegrid)

var_ts is replaced by dots_max (modules,init_3d_model)

Every cloud droplet has now an own weighting factor and can be deleted due to
collisions. Condensation and collision of cloud droplets are adjusted
accordingly. (advec_particles)

Collision efficiency for large cloud droplets has changed according to table of
Rogers and Yau. (collision_efficiency)

Errors:


Bugfix for generating serial jobs (subjob)

Bugfix: index problem concerning gradient_level indices removed (header)

Dimension of array stat in cascade change to prevent type problems with
mpi2 libraries (poisfft_hybrid)

Loop was split to make runs reproducible when using ifort compiler.
(disturb_field)

Bugfix: exchange of ghost points for prho included (time_integration)

Bugfix: calculation of time-averaged surface heatfluxes (sum_up_3d_data)

Bugfix: calculation of precipitation_rate (calc_precipitation)

Bugfix: initial data assignments to some dvrp arrays changed due to error
messages from gfortran compiler (modules)

Bugfix: calculation of cloud droplet velocity (advec_particles)

Bugfix: transfer of particles at south/left edge (advec_particles)

Bugfix: calculation of collision_efficiency (collision_efficiency)

Bugfix: initialisation of var_mod (subsidence)

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