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

Rev	Line
[1]	1	#if defined( __ibmy_special )
	2	@PROCESS NOOPTimize
	3	#endif
	4	SUBROUTINE init_3d_model
	5
	6	!------------------------------------------------------------------------------!
	7	! Actual revisions:
	8	! -----------------
[98]	9	!
[77]	10	!
	11	! Former revisions:
	12	! -----------------
	13	! $Id: init_3d_model.f90 98 2007-06-21 09:36:33Z raasch $
	14	!
[98]	15	! 97 2007-06-21 08:23:15Z raasch
	16	! Initialization of salinity, call of init_ocean
	17	!
[90]	18	! 87 2007-05-22 15:46:47Z raasch
	19	! var_hom and var_sum renamed pr_palm
	20	!
[77]	21	! 75 2007-03-22 09:54:05Z raasch
[73]	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
[75]	24	! subdomain, moisture renamed humidity,
	25	! new initializing action "by_user" calls user_init_3d_model,
[72]	26	! precipitation_amount/rate, ts_value are allocated, +module netcdf_control,
[51]	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)
[75]	29	! -uvmean_outflow, uxrp, vynp eliminated
[1]	30	!
[39]	31	! 19 2007-02-23 04:53:48Z raasch
	32	! +handling of top fluxes
	33	!
[3]	34	! RCS Log replace by Id keyword, revision history cleaned up
	35	!
[1]	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
[72]	55	USE cloud_parameters
[1]	56	USE constants
	57	USE control_parameters
	58	USE cpulog
	59	USE indices
	60	USE interfaces
	61	USE model_1d
[51]	62	USE netcdf_control
[1]	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) )
[75]	88	ALLOCATE( rdf(nzb+1:nzt) )
[87]	89	ALLOCATE( hom_sum(nzb:nzt+1,pr_palm+max_pr_user,0:statistic_regions), &
[1]	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), &
[87]	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), &
[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), &
[48]	97	sums_wsts_bc_l(nzb:nzt+1,0:statistic_regions), &
	98	ts_value(var_ts,0:statistic_regions) )
[1]	99	ALLOCATE( km_damp_x(nxl-1:nxr+1), km_damp_y(nys-1:nyn+1) )
	100
[19]	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) )
[1]	105
	106	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
	107	!
	108	!-- Leapfrog scheme needs two timelevels of diffusion quantities
[19]	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), &
[1]	113	vsws_2(nys-1:nyn+1,nxl-1:nxr+1) )
	114	ENDIF
	115
[75]	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), &
[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
[75]	142	IF ( humidity .OR. passive_scalar ) THEN
[1]	143	!
[75]	144	!-- 2D-humidity/scalar arrays
[1]	145	ALLOCATE ( qs(nys-1:nyn+1,nxl-1:nxr+1), &
[19]	146	qsws_1(nys-1:nyn+1,nxl-1:nxr+1), &
	147	qswst_1(nys-1:nyn+1,nxl-1:nxr+1) )
[1]	148
	149	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
[19]	150	ALLOCATE( qsws_2(nys-1:nyn+1,nxl-1:nxr+1), &
	151	qswst_2(nys-1:nyn+1,nxl-1:nxr+1) )
[1]	152	ENDIF
	153	!
[75]	154	!-- 3D-humidity/scalar arrays
[1]	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	!
[75]	160	!-- 3D-arrays needed for humidity only
	161	IF ( humidity ) THEN
[1]	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) )
[72]	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) )
[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
[94]	192	IF ( ocean ) THEN
[95]	193	ALLOCATE( saswsb_1(nys-1:nyn+1,nxl-1:nxr+1), &
	194	saswst_1(nys-1:nyn+1,nxl-1:nxr+1) )
[96]	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), &
[94]	198	sa_3(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
[96]	199	rho => rho_1 ! routine calc_mean_profile requires density to be a
	200	! pointer
[94]	201	ENDIF
	202
[1]	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	!
[51]	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	!
[73]	231	!-- Arrays to store velocity data from t-dt needed for radiation boundary
	232	!-- conditions
	233	IF ( outflow_l ) THEN
[75]	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) )
[73]	237	ENDIF
	238	IF ( outflow_r ) THEN
[75]	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) )
[73]	242	ENDIF
	243	IF ( outflow_s ) THEN
[75]	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) )
[73]	247	ENDIF
	248	IF ( outflow_n ) THEN
[75]	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) )
[73]	252	ENDIF
	253
	254	!
[1]	255	!-- Initial assignment of the pointers
	256	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
	257
[19]	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
[1]	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
[75]	271	IF ( humidity .OR. passive_scalar ) THEN
[19]	272	qsws_m => qsws_1; qsws => qsws_2
	273	qswst_m => qswst_1; qswst => qswst_2
[1]	274	q_m => q_1; q => q_2; q_p => q_3; tq_m => q_3
[75]	275	IF ( humidity ) vpt_m => vpt_1; vpt => vpt_2
[1]	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
[19]	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
[1]	297
[75]	298	IF ( humidity .OR. passive_scalar ) THEN
[1]	299	qsws => qsws_1
[19]	300	qswst => qswst_1
[94]	301	q => q_1; q_p => q_2; tq_m => q_3; q_m => q_3
[75]	302	IF ( humidity ) vpt => vpt_1
[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
[94]	310	IF ( ocean ) THEN
[95]	311	saswsb => saswsb_1
[94]	312	saswst => saswst_1
	313	sa => sa_1; sa_p => sa_2; tsa_m => sa_3
	314	ENDIF
	315
[1]	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
[75]	342	IF ( humidity .OR. passive_scalar ) THEN
[1]	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!
[75]	388	IF ( humidity .OR. passive_scalar ) qs = 0.0
[1]	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
[75]	443
[1]	444	!
[51]	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
[1]	447	!-- limit in the initial phase of a run (at k=1, dz/2 occurs in the
[51]	448	!-- limiting formula!). The original values are stored to be later used for
	449	!-- volume flow control.
[1]	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
[51]	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
[1]	470
[75]	471	IF ( humidity .OR. passive_scalar ) THEN
[1]	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
[94]	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
[1]	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
[75]	502	IF ( humidity .OR. passive_scalar ) qs = 0.0
[1]	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
[46]	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
[1]	515	ENDIF
	516
	517	!
	518	!-- Calculate virtual potential temperature
[75]	519	IF ( humidity ) vpt = pt * ( 1.0 + 0.61 * q )
[1]	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
[97]	534	IF ( ocean ) THEN
	535	!
	536	!-- Store initial salinity profile
	537	hom(:,1,26,:) = SPREAD( sa(:,nys,nxl), 2, statistic_regions+1 )
	538	ENDIF
[1]	539
[75]	540	IF ( humidity ) THEN
[1]	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	!
[19]	562	!-- Initialize fluxes at bottom surface
[1]	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
[75]	589	IF ( humidity .OR. passive_scalar ) THEN
[1]	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	!
[19]	599	!-- Initialize fluxes at top surface
[94]	600	!-- Currently, only the heatflux and salinity flux can be prescribed.
	601	!-- The latent flux is zero in this case!
[19]	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
[75]	610	IF ( humidity .OR. passive_scalar ) THEN
[19]	611	qswst = 0.0
	612	IF ( ASSOCIATED( qswst_m ) ) qswst_m = qswst
	613	ENDIF
[94]	614
	615	IF ( ocean ) THEN
[95]	616	saswsb = bottom_salinityflux
[94]	617	saswst = top_salinityflux
	618	ENDIF
[19]	619	ENDIF
	620
	621	ENDIF
	622
	623	!
[1]	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
[75]	640	IF ( humidity .OR. passive_scalar ) THEN
[1]	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
[51]	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)
[1]	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
[51]	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)
[1]	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	!
[72]	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
[1]	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
[75]	739	IF ( ( humidity .OR. passive_scalar ) .AND. &
[1]	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
[75]	752	CALL disturb_field( nzb_u_inner, tend, u )
	753	CALL disturb_field( nzb_v_inner, tend, v )
[1]	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
[75]	776	IF ( humidity .OR. passive_scalar ) THEN
[1]	777	IF ( ASSOCIATED( q_m ) ) q_m = q
	778	IF ( timestep_scheme(1:5) == 'runge' ) tq_m = 0.0
	779	q_p = q
[75]	780	IF ( humidity .AND. ASSOCIATED( vpt_m ) ) vpt_m = vpt
[1]	781	ENDIF
	782
[94]	783	IF ( ocean ) THEN
	784	tsa_m = 0.0
	785	sa_p = sa
	786	ENDIF
	787
[73]	788	!
	789	!-- Initialize old timelevels needed for radiation boundary conditions
	790	IF ( outflow_l ) THEN
[75]	791	u_m_l(:,:,:) = u(:,:,-1:1)
	792	v_m_l(:,:,:) = v(:,:,-1:1)
	793	w_m_l(:,:,:) = w(:,:,-1:1)
[73]	794	ENDIF
	795	IF ( outflow_r ) THEN
[75]	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)
[73]	799	ENDIF
	800	IF ( outflow_s ) THEN
[75]	801	u_m_s(:,:,:) = u(:,-1:1,:)
	802	v_m_s(:,:,:) = v(:,-1:1,:)
	803	w_m_s(:,:,:) = w(:,-1:1,:)
[73]	804	ENDIF
	805	IF ( outflow_n ) THEN
[75]	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,:)
[73]	809	ENDIF
	810
[1]	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
[75]	826	IF ( humidity .OR. passive_scalar ) q_p = q
[94]	827	IF ( ocean ) sa_p = sa
[1]	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
[96]	841	IF ( ocean ) THEN
[1]	842	!
[96]	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
[1]	853
	854	!
	855	!-- If required, initialize particles
[63]	856	IF ( particle_advection ) CALL init_particles
[1]	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
[73]	880	IF ( bc_lr == 'dirichlet/radiation' ) THEN
[1]	881
	882	DO i = nxl-1, nxr+1
[73]	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
[1]	892
[73]	893	ELSEIF ( bc_lr == 'radiation/dirichlet' ) THEN
[1]	894
[73]	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
[1]	905
[73]	906	ENDIF
[1]	907
[73]	908	IF ( bc_ns == 'dirichlet/radiation' ) THEN
[1]	909
[73]	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
[1]	918	ENDIF
	919	ENDDO
	920
[73]	921	ELSEIF ( bc_ns == 'radiation/dirichlet' ) THEN
[1]	922
	923	DO j = nys-1, nyn+1
[73]	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
[1]	931	ENDIF
[73]	932	ENDDO
[1]	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	!
[51]	952	!-- User-defined initializing actions. Check afterwards, if maximum number
	953	!-- of allowed timeseries is not exceeded
[1]	954	CALL user_init
	955
[51]	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
[1]	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
[87]	985	ngp_sums = ( nz + 2 ) * ( pr_palm + max_pr_user )
[1]	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 |