Home

Context Navigation

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

Last change on this file since 100 was 98, checked in by raasch, 17 years ago
updating comments and rc-file
Property svn:keywords set to `Id`
File size: 34.9 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |