Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |