Home

Context Navigation

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

Last change on this file since 72 was 72, checked in by raasch, 17 years ago
preliminary changes for precipitation output
Property svn:keywords set to `Id`
File size: 32.2 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |