Home

Context Navigation

init_3d_model.f90 @ 91

Last change on this file since 91 was 90, checked in by raasch, 17 years ago

New:
---
Calculation and output of user-defined profiles. New &userpar parameters data_output_pr_user and max_pr_user.

check_parameters, flow_statistics, modules, parin, read_var_list, user_interface, write_var_list

Changed:

Division through dt_3d replaced by multiplication of the inverse. For performance optimisation, this is done in the loop calculating the divergence instead of using a seperate loop. (pres.f90) var_hom and var_sum renamed pr_palm.

data_output_profiles, flow_statistics, init_3d_model, modules, parin, pres, read_var_list, run_control, time_integration

Errors:

Bugfix: work_fft*_vec removed from some PRIVATE-declarations (poisfft).

Bugfix: field_chr renamed field_char (user_interface).

Bugfix: output of use_upstream_for_tke (header).

header, poisfft, user_interface

Property svn:keywords set to Id

File size: 33.5 KB

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

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

Context Navigation

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

Download in other formats: