Home

Context Navigation

source: palm/trunk/SOURCE/flow_statistics.f90 @ 88

Last change on this file since 88 was 87, checked in by raasch, 18 years ago
Preliminary update for user defined profiles
Property svn:keywords set to `Id`
File size: 37.9 KB

Rev	Line
[1]	1	SUBROUTINE flow_statistics
	2
	3	!------------------------------------------------------------------------------!
	4	! Actual revisions:
	5	! -----------------
[87]	6	! Two more arguments added to user_statistics, which is now also called for
	7	! user-defined profiles,
	8	! var_hom and var_sum renamed pr_palm
[1]	9	!
	10	! Former revisions:
	11	! -----------------
[3]	12	! $Id: flow_statistics.f90 87 2007-05-22 15:46:47Z raasch $
[39]	13	!
[83]	14	! 82 2007-04-16 15:40:52Z raasch
	15	! Cpp-directive lcmuk changed to intel_openmp_bug
	16	!
[77]	17	! 75 2007-03-22 09:54:05Z raasch
	18	! Collection of time series quantities moved from routine flow_statistics to
	19	! here, routine user_statistics is called for each statistic region,
	20	! moisture renamed humidity
	21	!
[39]	22	! 19 2007-02-23 04:53:48Z raasch
[77]	23	! fluxes at top modified (tswst, qswst)
[39]	24	!
[3]	25	! RCS Log replace by Id keyword, revision history cleaned up
	26	!
[1]	27	! Revision 1.41 2006/08/04 14:37:50 raasch
	28	! Error removed in non-parallel part (sums_l)
	29	!
	30	! Revision 1.1 1997/08/11 06:15:17 raasch
	31	! Initial revision
	32	!
	33	!
	34	! Description:
	35	! ------------
	36	! Compute average profiles and further average flow quantities for the different
	37	! user-defined (sub-)regions. The region indexed 0 is the total model domain.
	38	!
	39	! NOTE: For simplicity, nzb_s_outer and nzb_diff_s_outer are being used as a
	40	! ---- lower vertical index for k-loops for all variables so that regardless
	41	! of the variable and its respective staggered grid always the same number of
	42	! grid points is used for 2D averages. The disadvantage: depending on the
	43	! variable, up to one grid layer adjacent to the (vertical walls of the)
	44	! topography is missed out by this simplification.
	45	!------------------------------------------------------------------------------!
	46
	47	USE arrays_3d
	48	USE cloud_parameters
	49	USE cpulog
	50	USE grid_variables
	51	USE indices
	52	USE interfaces
	53	USE pegrid
	54	USE statistics
	55	USE control_parameters
	56
	57	IMPLICIT NONE
	58
	59	INTEGER :: i, j, k, omp_get_thread_num, sr, tn
	60	LOGICAL :: first
	61	REAL :: height, pts, sums_l_eper, sums_l_etot, ust, ust2, u2, vst, &
	62	vst2, v2, w2, z_i(2)
	63	REAL :: sums_ll(nzb:nzt+1,2)
	64
	65
	66	CALL cpu_log( log_point(10), 'flow_statistics', 'start' )
	67
	68	!
	69	!-- To be on the safe side, check whether flow_statistics has already been
	70	!-- called once after the current time step
	71	IF ( flow_statistics_called ) THEN
	72	IF ( myid == 0 ) PRINT*, '+++ WARNING: flow_statistics is called two', &
	73	' times within one timestep'
	74	CALL local_stop
	75	ENDIF
	76
	77	!
	78	!-- Compute statistics for each (sub-)region
	79	DO sr = 0, statistic_regions
	80
	81	!
	82	!-- Initialize (local) summation array
	83	sums_l = 0.0
	84
	85	!
	86	!-- Store sums that have been computed in other subroutines in summation
	87	!-- array
	88	sums_l(:,11,:) = sums_l_l(:,sr,:) ! mixing length from diffusivities
	89	!-- WARNING: next line still has to be adjusted for OpenMP
	90	sums_l(:,21,0) = sums_wsts_bc_l(:,sr) ! heat flux from advec_s_bc
[87]	91	sums_l(nzb+9,pr_palm,0) = sums_divold_l(sr) ! old divergence from pres
	92	sums_l(nzb+10,pr_palm,0) = sums_divnew_l(sr) ! new divergence from pres
[1]	93	!-- WARNING: next four lines still may have to be adjusted for OpenMP
[87]	94	sums_l(nzb:nzb+2,pr_palm-1,0) = sums_up_fraction_l(1,1:3,sr)! upstream
	95	sums_l(nzb+3:nzb+5,pr_palm-1,0) = sums_up_fraction_l(2,1:3,sr)! parts
	96	sums_l(nzb+6:nzb+8,pr_palm-1,0) = sums_up_fraction_l(3,1:3,sr)! from
	97	sums_l(nzb+9:nzb+11,pr_palm-1,0) = sums_up_fraction_l(4,1:3,sr)! spline
[1]	98
	99	!
	100	!-- Horizontally averaged profiles of horizontal velocities and temperature.
	101	!-- They must have been computed before, because they are already required
	102	!-- for other horizontal averages.
	103	tn = 0
	104	!$OMP PARALLEL PRIVATE( i, j, k, tn )
[82]	105	#if defined( __intel_openmp_bug )
[1]	106	tn = omp_get_thread_num()
	107	#else
	108	!$ tn = omp_get_thread_num()
	109	#endif
	110
	111	!$OMP DO
	112	DO i = nxl, nxr
	113	DO j = nys, nyn
	114	DO k = nzb_s_outer(j,i), nzt+1
	115	sums_l(k,1,tn) = sums_l(k,1,tn) + u(k,j,i) * rmask(j,i,sr)
	116	sums_l(k,2,tn) = sums_l(k,2,tn) + v(k,j,i) * rmask(j,i,sr)
	117	sums_l(k,4,tn) = sums_l(k,4,tn) + pt(k,j,i) * rmask(j,i,sr)
	118	ENDDO
	119	ENDDO
	120	ENDDO
	121
	122	!
	123	!-- Horizontally averaged profiles of virtual potential temperature,
	124	!-- total water content, specific humidity and liquid water potential
	125	!-- temperature
[75]	126	IF ( humidity ) THEN
[1]	127	!$OMP DO
	128	DO i = nxl, nxr
	129	DO j = nys, nyn
	130	DO k = nzb_s_outer(j,i), nzt+1
	131	sums_l(k,44,tn) = sums_l(k,44,tn) + &
	132	vpt(k,j,i) * rmask(j,i,sr)
	133	sums_l(k,41,tn) = sums_l(k,41,tn) + &
	134	q(k,j,i) * rmask(j,i,sr)
	135	ENDDO
	136	ENDDO
	137	ENDDO
	138	IF ( cloud_physics ) THEN
	139	!$OMP DO
	140	DO i = nxl, nxr
	141	DO j = nys, nyn
	142	DO k = nzb_s_outer(j,i), nzt+1
	143	sums_l(k,42,tn) = sums_l(k,42,tn) + &
	144	( q(k,j,i) - ql(k,j,i) ) * rmask(j,i,sr)
	145	sums_l(k,43,tn) = sums_l(k,43,tn) + ( &
	146	pt(k,j,i) + l_d_cppt_d_t(k) ql(k,j,i) &
	147	) * rmask(j,i,sr)
	148	ENDDO
	149	ENDDO
	150	ENDDO
	151	ENDIF
	152	ENDIF
	153
	154	!
	155	!-- Horizontally averaged profiles of passive scalar
	156	IF ( passive_scalar ) THEN
	157	!$OMP DO
	158	DO i = nxl, nxr
	159	DO j = nys, nyn
	160	DO k = nzb_s_outer(j,i), nzt+1
	161	sums_l(k,41,tn) = sums_l(k,41,tn) + q(k,j,i) * rmask(j,i,sr)
	162	ENDDO
	163	ENDDO
	164	ENDDO
	165	ENDIF
	166	!$OMP END PARALLEL
	167
	168	!
	169	!-- Summation of thread sums
	170	IF ( threads_per_task > 1 ) THEN
	171	DO i = 1, threads_per_task-1
	172	sums_l(:,1,0) = sums_l(:,1,0) + sums_l(:,1,i)
	173	sums_l(:,2,0) = sums_l(:,2,0) + sums_l(:,2,i)
	174	sums_l(:,4,0) = sums_l(:,4,0) + sums_l(:,4,i)
[75]	175	IF ( humidity ) THEN
[1]	176	sums_l(:,41,0) = sums_l(:,41,0) + sums_l(:,41,i)
	177	sums_l(:,44,0) = sums_l(:,44,0) + sums_l(:,44,i)
	178	IF ( cloud_physics ) THEN
	179	sums_l(:,42,0) = sums_l(:,42,0) + sums_l(:,42,i)
	180	sums_l(:,43,0) = sums_l(:,43,0) + sums_l(:,43,i)
	181	ENDIF
	182	ENDIF
	183	IF ( passive_scalar ) THEN
	184	sums_l(:,41,0) = sums_l(:,41,0) + sums_l(:,41,i)
	185	ENDIF
	186	ENDDO
	187	ENDIF
	188
	189	#if defined( __parallel )
	190	!
	191	!-- Compute total sum from local sums
	192	CALL MPI_ALLREDUCE( sums_l(nzb,1,0), sums(nzb,1), nzt+2-nzb, MPI_REAL, &
	193	MPI_SUM, comm2d, ierr )
	194	CALL MPI_ALLREDUCE( sums_l(nzb,2,0), sums(nzb,2), nzt+2-nzb, MPI_REAL, &
	195	MPI_SUM, comm2d, ierr )
	196	CALL MPI_ALLREDUCE( sums_l(nzb,4,0), sums(nzb,4), nzt+2-nzb, MPI_REAL, &
	197	MPI_SUM, comm2d, ierr )
[75]	198	IF ( humidity ) THEN
[1]	199	CALL MPI_ALLREDUCE( sums_l(nzb,44,0), sums(nzb,44), nzt+2-nzb, &
	200	MPI_REAL, MPI_SUM, comm2d, ierr )
	201	CALL MPI_ALLREDUCE( sums_l(nzb,41,0), sums(nzb,41), nzt+2-nzb, &
	202	MPI_REAL, MPI_SUM, comm2d, ierr )
	203	IF ( cloud_physics ) THEN
	204	CALL MPI_ALLREDUCE( sums_l(nzb,42,0), sums(nzb,42), nzt+2-nzb, &
	205	MPI_REAL, MPI_SUM, comm2d, ierr )
	206	CALL MPI_ALLREDUCE( sums_l(nzb,43,0), sums(nzb,43), nzt+2-nzb, &
	207	MPI_REAL, MPI_SUM, comm2d, ierr )
	208	ENDIF
	209	ENDIF
	210
	211	IF ( passive_scalar ) THEN
	212	CALL MPI_ALLREDUCE( sums_l(nzb,41,0), sums(nzb,41), nzt+2-nzb, &
	213	MPI_REAL, MPI_SUM, comm2d, ierr )
	214	ENDIF
	215	#else
	216	sums(:,1) = sums_l(:,1,0)
	217	sums(:,2) = sums_l(:,2,0)
	218	sums(:,4) = sums_l(:,4,0)
[75]	219	IF ( humidity ) THEN
[1]	220	sums(:,44) = sums_l(:,44,0)
	221	sums(:,41) = sums_l(:,41,0)
	222	IF ( cloud_physics ) THEN
	223	sums(:,42) = sums_l(:,42,0)
	224	sums(:,43) = sums_l(:,43,0)
	225	ENDIF
	226	ENDIF
	227	IF ( passive_scalar ) sums(:,41) = sums_l(:,41,0)
	228	#endif
	229
	230	!
	231	!-- Final values are obtained by division by the total number of grid points
	232	!-- used for summation. After that store profiles.
	233	sums(:,1) = sums(:,1) / ngp_2dh_outer(:,sr)
	234	sums(:,2) = sums(:,2) / ngp_2dh_outer(:,sr)
	235	sums(:,4) = sums(:,4) / ngp_2dh_outer(:,sr)
	236	hom(:,1,1,sr) = sums(:,1) ! u
	237	hom(:,1,2,sr) = sums(:,2) ! v
	238	hom(:,1,4,sr) = sums(:,4) ! pt
	239
	240	!
	241	!-- Humidity and cloud parameters
[75]	242	IF ( humidity ) THEN
[1]	243	sums(:,44) = sums(:,44) / ngp_2dh_outer(:,sr)
	244	sums(:,41) = sums(:,41) / ngp_2dh_outer(:,sr)
	245	hom(:,1,44,sr) = sums(:,44) ! vpt
	246	hom(:,1,41,sr) = sums(:,41) ! qv (q)
	247	IF ( cloud_physics ) THEN
	248	sums(:,42) = sums(:,42) / ngp_2dh_outer(:,sr)
	249	sums(:,43) = sums(:,43) / ngp_2dh_outer(:,sr)
	250	hom(:,1,42,sr) = sums(:,42) ! qv
	251	hom(:,1,43,sr) = sums(:,43) ! pt
	252	ENDIF
	253	ENDIF
	254
	255	!
	256	!-- Passive scalar
	257	IF ( passive_scalar ) hom(:,1,41,sr) = sums(:,41) / ngp_2dh_outer(:,sr)
	258
	259	!
	260	!-- Horizontally averaged profiles of the remaining prognostic variables,
	261	!-- variances, the total and the perturbation energy (single values in last
	262	!-- column of sums_l) and some diagnostic quantities.
	263	!-- NOTE: for simplicity, nzb_s_outer is used below, although strictly
	264	!-- ---- speaking the following k-loop would have to be split up and
	265	!-- rearranged according to the staggered grid.
	266	tn = 0
[82]	267	#if defined( __intel_openmp_bug )
[1]	268	!$OMP PARALLEL PRIVATE( i, j, k, pts, sums_ll, sums_l_eper, sums_l_etot, &
	269	!$OMP tn, ust, ust2, u2, vst, vst2, v2, w2 )
	270	tn = omp_get_thread_num()
	271	#else
	272	!$OMP PARALLEL PRIVATE( i, j, k, pts, sums_ll, sums_l_eper, sums_l_etot, tn, ust, ust2, u2, vst, vst2, v2, w2 )
	273	!$ tn = omp_get_thread_num()
	274	#endif
	275	!$OMP DO
	276	DO i = nxl, nxr
	277	DO j = nys, nyn
	278	sums_l_etot = 0.0
	279	sums_l_eper = 0.0
	280	DO k = nzb_s_outer(j,i), nzt+1
	281	u2 = u(k,j,i)**2
	282	v2 = v(k,j,i)**2
	283	w2 = w(k,j,i)**2
	284	ust2 = ( u(k,j,i) - hom(k,1,1,sr) )**2
	285	vst2 = ( v(k,j,i) - hom(k,1,2,sr) )**2
	286	!
	287	!-- Prognostic and diagnostic variables
	288	sums_l(k,3,tn) = sums_l(k,3,tn) + w(k,j,i) * rmask(j,i,sr)
	289	sums_l(k,8,tn) = sums_l(k,8,tn) + e(k,j,i) * rmask(j,i,sr)
	290	sums_l(k,9,tn) = sums_l(k,9,tn) + km(k,j,i) * rmask(j,i,sr)
	291	sums_l(k,10,tn) = sums_l(k,10,tn) + kh(k,j,i) * rmask(j,i,sr)
	292	sums_l(k,40,tn) = sums_l(k,40,tn) + p(k,j,i)
	293
	294	!
	295	!-- Variances
	296	sums_l(k,30,tn) = sums_l(k,30,tn) + ust2 * rmask(j,i,sr)
	297	sums_l(k,31,tn) = sums_l(k,31,tn) + vst2 * rmask(j,i,sr)
	298	sums_l(k,32,tn) = sums_l(k,32,tn) + w2 * rmask(j,i,sr)
	299	sums_l(k,33,tn) = sums_l(k,33,tn) + &
	300	( pt(k,j,i)-hom(k,1,4,sr) )*2 rmask(j,i,sr)
	301	!
	302	!-- Higher moments
	303	!-- (Computation of the skewness of w further below)
	304	sums_l(k,38,tn) = sums_l(k,38,tn) + w(k,j,i) * w2 * &
	305	rmask(j,i,sr)
	306	!
	307	!-- Perturbation energy
	308	sums_l(k,34,tn) = sums_l(k,34,tn) + 0.5 * ( ust2 + vst2 + w2 ) &
	309	* rmask(j,i,sr)
	310	sums_l_etot = sums_l_etot + &
	311	0.5 * ( u2 + v2 + w2 ) * rmask(j,i,sr)
	312	sums_l_eper = sums_l_eper + &
	313	0.5 * ( ust2+vst2+w2 ) * rmask(j,i,sr)
	314	ENDDO
	315	!
	316	!-- Total and perturbation energy for the total domain (being
	317	!-- collected in the last column of sums_l). Summation of these
	318	!-- quantities is seperated from the previous loop in order to
	319	!-- allow vectorization of that loop.
[87]	320	sums_l(nzb+4,pr_palm,tn) = sums_l(nzb+4,pr_palm,tn) + sums_l_etot
	321	sums_l(nzb+5,pr_palm,tn) = sums_l(nzb+5,pr_palm,tn) + sums_l_eper
[1]	322	!
	323	!-- 2D-arrays (being collected in the last column of sums_l)
[87]	324	sums_l(nzb,pr_palm,tn) = sums_l(nzb,pr_palm,tn) + &
[1]	325	us(j,i) * rmask(j,i,sr)
[87]	326	sums_l(nzb+1,pr_palm,tn) = sums_l(nzb+1,pr_palm,tn) + &
[1]	327	usws(j,i) * rmask(j,i,sr)
[87]	328	sums_l(nzb+2,pr_palm,tn) = sums_l(nzb+2,pr_palm,tn) + &
[1]	329	vsws(j,i) * rmask(j,i,sr)
[87]	330	sums_l(nzb+3,pr_palm,tn) = sums_l(nzb+3,pr_palm,tn) + &
[1]	331	ts(j,i) * rmask(j,i,sr)
	332	ENDDO
	333	ENDDO
	334
	335	!
	336	!-- Horizontally averaged profiles of the vertical fluxes
	337	!$OMP DO
	338	DO i = nxl, nxr
	339	DO j = nys, nyn
	340	!
	341	!-- Subgridscale fluxes (without Prandtl layer from k=nzb,
	342	!-- oterwise from k=nzb+1)
	343	!-- NOTE: for simplicity, nzb_diff_s_outer is used below, although
	344	!-- ---- strictly speaking the following k-loop would have to be
	345	!-- split up according to the staggered grid.
	346	DO k = nzb_diff_s_outer(j,i)-1, nzt
	347	!
	348	!-- Momentum flux w"u"
	349	sums_l(k,12,tn) = sums_l(k,12,tn) - 0.25 * ( &
	350	km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) &
	351	) * ( &
	352	( u(k+1,j,i) - u(k,j,i) ) * ddzu(k+1) &
	353	+ ( w(k,j,i) - w(k,j,i-1) ) * ddx &
	354	) * rmask(j,i,sr)
	355	!
	356	!-- Momentum flux w"v"
	357	sums_l(k,14,tn) = sums_l(k,14,tn) - 0.25 * ( &
	358	km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) &
	359	) * ( &
	360	( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1) &
	361	+ ( w(k,j,i) - w(k,j-1,i) ) * ddy &
	362	) * rmask(j,i,sr)
[19]	363	ENDDO
	364
	365	DO k = nzb_diff_s_outer(j,i)-1, nzt_diff
[1]	366	!
	367	!-- Heat flux w"pt"
	368	sums_l(k,16,tn) = sums_l(k,16,tn) &
	369	- 0.5 * ( kh(k,j,i) + kh(k+1,j,i) ) &
	370	* ( pt(k+1,j,i) - pt(k,j,i) ) &
	371	* ddzu(k+1) * rmask(j,i,sr)
	372
	373
	374	!
	375	!-- Buoyancy flux, water flux (humidity flux) w"q"
[75]	376	IF ( humidity ) THEN
[1]	377	sums_l(k,45,tn) = sums_l(k,45,tn) &
	378	- 0.5 * ( kh(k,j,i) + kh(k+1,j,i) ) &
	379	* ( vpt(k+1,j,i) - vpt(k,j,i) ) &
	380	* ddzu(k+1) * rmask(j,i,sr)
	381	sums_l(k,48,tn) = sums_l(k,48,tn) &
	382	- 0.5 * ( kh(k,j,i) + kh(k+1,j,i) ) &
	383	* ( q(k+1,j,i) - q(k,j,i) ) &
	384	* ddzu(k+1) * rmask(j,i,sr)
	385	IF ( cloud_physics ) THEN
	386	sums_l(k,51,tn) = sums_l(k,51,tn) &
	387	- 0.5 * ( kh(k,j,i) + kh(k+1,j,i) ) &
	388	* ( ( q(k+1,j,i) - ql(k+1,j,i) )&
	389	- ( q(k,j,i) - ql(k,j,i) ) ) &
	390	* ddzu(k+1) * rmask(j,i,sr)
	391	ENDIF
	392	ENDIF
	393
	394	!
	395	!-- Passive scalar flux
	396	IF ( passive_scalar ) THEN
	397	sums_l(k,48,tn) = sums_l(k,48,tn) &
	398	- 0.5 * ( kh(k,j,i) + kh(k+1,j,i) ) &
	399	* ( q(k+1,j,i) - q(k,j,i) ) &
	400	* ddzu(k+1) * rmask(j,i,sr)
	401	ENDIF
	402
	403	ENDDO
	404
	405	!
	406	!-- Subgridscale fluxes in the Prandtl layer
	407	IF ( use_surface_fluxes ) THEN
	408	sums_l(nzb,12,tn) = sums_l(nzb,12,tn) + &
	409	usws(j,i) * rmask(j,i,sr) ! w"u"
	410	sums_l(nzb,14,tn) = sums_l(nzb,14,tn) + &
	411	vsws(j,i) * rmask(j,i,sr) ! w"v"
	412	sums_l(nzb,16,tn) = sums_l(nzb,16,tn) + &
	413	shf(j,i) * rmask(j,i,sr) ! w"pt"
	414	sums_l(nzb,58,tn) = sums_l(nzb,58,tn) + &
	415	0.0 * rmask(j,i,sr) ! u"pt"
	416	sums_l(nzb,61,tn) = sums_l(nzb,61,tn) + &
	417	0.0 * rmask(j,i,sr) ! v"pt"
[75]	418	IF ( humidity ) THEN
[1]	419	sums_l(nzb,48,tn) = sums_l(nzb,48,tn) + &
	420	qsws(j,i) * rmask(j,i,sr) ! w"q" (w"qv")
	421	IF ( cloud_physics ) THEN
	422	sums_l(nzb,45,tn) = sums_l(nzb,45,tn) + ( &
	423	( 1.0 + 0.61 * q(nzb,j,i) ) * &
	424	shf(j,i) + 0.61 * pt(nzb,j,i) * &
	425	qsws(j,i) &
	426	)
	427	!
	428	!-- Formula does not work if ql(nzb) /= 0.0
	429	sums_l(nzb,51,tn) = sums_l(nzb,51,tn) + & ! w"q" (w"qv")
	430	qsws(j,i) * rmask(j,i,sr)
	431	ENDIF
	432	ENDIF
	433	IF ( passive_scalar ) THEN
	434	sums_l(nzb,48,tn) = sums_l(nzb,48,tn) + &
	435	qsws(j,i) * rmask(j,i,sr) ! w"q" (w"qv")
	436	ENDIF
	437	ENDIF
	438
	439	!
[19]	440	!-- Subgridscale fluxes at the top surface
	441	IF ( use_top_fluxes ) THEN
	442	sums_l(nzt,16,tn) = sums_l(nzt,16,tn) + &
	443	tswst(j,i) * rmask(j,i,sr) ! w"pt"
	444	sums_l(nzt,58,tn) = sums_l(nzt,58,tn) + &
	445	0.0 * rmask(j,i,sr) ! u"pt"
	446	sums_l(nzt,61,tn) = sums_l(nzt,61,tn) + &
	447	0.0 * rmask(j,i,sr) ! v"pt"
[75]	448	IF ( humidity ) THEN
[19]	449	sums_l(nzt,48,tn) = sums_l(nzt,48,tn) + &
	450	qswst(j,i) * rmask(j,i,sr) ! w"q" (w"qv")
	451	IF ( cloud_physics ) THEN
	452	sums_l(nzt,45,tn) = sums_l(nzt,45,tn) + ( &
	453	( 1.0 + 0.61 * q(nzt,j,i) ) * &
	454	tswst(j,i) + 0.61 * pt(nzt,j,i) * &
	455	qsws(j,i) &
	456	)
	457	!
	458	!-- Formula does not work if ql(nzb) /= 0.0
	459	sums_l(nzt,51,tn) = sums_l(nzt,51,tn) + & ! w"q" (w"qv")
	460	qswst(j,i) * rmask(j,i,sr)
	461	ENDIF
	462	ENDIF
	463	IF ( passive_scalar ) THEN
	464	sums_l(nzt,48,tn) = sums_l(nzt,48,tn) + &
	465	qswst(j,i) * rmask(j,i,sr) ! w"q" (w"qv")
	466	ENDIF
	467	ENDIF
	468
	469	!
[1]	470	!-- Resolved fluxes (can be computed for all horizontal points)
	471	!-- NOTE: for simplicity, nzb_s_outer is used below, although strictly
	472	!-- ---- speaking the following k-loop would have to be split up and
	473	!-- rearranged according to the staggered grid.
	474	DO k = nzb_s_outer(j,i), nzt
	475	ust = 0.5 * ( u(k,j,i) - hom(k,1,1,sr) + &
	476	u(k+1,j,i) - hom(k+1,1,1,sr) )
	477	vst = 0.5 * ( v(k,j,i) - hom(k,1,2,sr) + &
	478	v(k+1,j,i) - hom(k+1,1,2,sr) )
	479	pts = 0.5 * ( pt(k,j,i) - hom(k,1,4,sr) + &
	480	pt(k+1,j,i) - hom(k+1,1,4,sr) )
	481	!
	482	!-- Momentum flux wu
	483	sums_l(k,13,tn) = sums_l(k,13,tn) + 0.5 * &
	484	( w(k,j,i-1) + w(k,j,i) ) &
	485	* ust * rmask(j,i,sr)
	486	!
	487	!-- Momentum flux wv
	488	sums_l(k,15,tn) = sums_l(k,15,tn) + 0.5 * &
	489	( w(k,j-1,i) + w(k,j,i) ) &
	490	* vst * rmask(j,i,sr)
	491	!
	492	!-- Heat flux wpt
	493	!-- The following formula (comment line, not executed) does not
	494	!-- work if applied to subregions
	495	! sums_l(k,17,tn) = sums_l(k,17,tn) + 0.5 * &
	496	! ( pt(k,j,i)+pt(k+1,j,i) ) &
	497	! * w(k,j,i) * rmask(j,i,sr)
	498	sums_l(k,17,tn) = sums_l(k,17,tn) + pts * w(k,j,i) * &
	499	rmask(j,i,sr)
	500	!
	501	!-- Higher moments
	502	sums_l(k,35,tn) = sums_l(k,35,tn) + pts * w(k,j,i)*2 &
	503	rmask(j,i,sr)
	504	sums_l(k,36,tn) = sums_l(k,36,tn) + pts*2 w(k,j,i) * &
	505	rmask(j,i,sr)
	506
	507	!
	508	!-- Buoyancy flux, water flux, humidity flux and liquid water
	509	!-- content
[75]	510	IF ( humidity ) THEN
[1]	511	pts = 0.5 * ( vpt(k,j,i) - hom(k,1,44,sr) + &
	512	vpt(k+1,j,i) - hom(k+1,1,44,sr) )
	513	sums_l(k,46,tn) = sums_l(k,46,tn) + pts * w(k,j,i) * &
	514	rmask(j,i,sr)
	515	pts = 0.5 * ( q(k,j,i) - hom(k,1,41,sr) + &
	516	q(k+1,j,i) - hom(k+1,1,41,sr) )
	517	sums_l(k,49,tn) = sums_l(k,49,tn) + pts * w(k,j,i) * &
	518	rmask(j,i,sr)
	519	IF ( cloud_physics .OR. cloud_droplets ) THEN
	520	pts = 0.5 * &
	521	( ( q(k,j,i) - ql(k,j,i) ) - hom(k,1,42,sr) &
	522	+ ( q(k+1,j,i) - ql(k+1,j,i) ) - hom(k+1,1,42,sr) )
	523	sums_l(k,52,tn) = sums_l(k,52,tn) + pts * w(k,j,i) * &
	524	rmask(j,i,sr)
	525	sums_l(k,54,tn) = sums_l(k,54,tn) + ql(k,j,i) * &
	526	rmask(j,i,sr)
	527	ENDIF
	528	ENDIF
	529
	530	!
	531	!-- Passive scalar flux
	532	IF ( passive_scalar ) THEN
	533	pts = 0.5 * ( q(k,j,i) - hom(k,1,41,sr) + &
	534	q(k+1,j,i) - hom(k+1,1,41,sr) )
	535	sums_l(k,49,tn) = sums_l(k,49,tn) + pts * w(k,j,i) * &
	536	rmask(j,i,sr)
	537	ENDIF
	538
	539	!
	540	!-- Energy flux we
	541	sums_l(k,37,tn) = sums_l(k,37,tn) + w(k,j,i) * 0.5 * &
	542	( ust2 + vst2 + w(k,j,i)**2 )&
	543	* rmask(j,i,sr)
	544
	545	ENDDO
	546	ENDDO
	547	ENDDO
	548
	549	!
	550	!-- Divergence of vertical flux of resolved scale energy and pressure
	551	!-- fluctuations. First calculate the products, then the divergence.
	552	!-- Calculation is time consuming. Do it only, if profiles shall be plotted.
	553	IF ( hom(nzb+1,2,55,0) /= 0.0 ) THEN
	554
	555	sums_ll = 0.0 ! local array
	556
	557	!$OMP DO
	558	DO i = nxl, nxr
	559	DO j = nys, nyn
	560	DO k = nzb_s_outer(j,i)+1, nzt
	561
	562	sums_ll(k,1) = sums_ll(k,1) + 0.5 * w(k,j,i) * ( &
	563	( 0.25 * ( u(k,j,i)+u(k+1,j,i)+u(k,j,i+1)+u(k+1,j,i+1) &
	564	- 2.0 * ( hom(k,1,1,sr) + hom(k+1,1,1,sr) ) &
	565	) )**2 &
	566	+ ( 0.25 * ( v(k,j,i)+v(k+1,j,i)+v(k,j+1,i)+v(k+1,j+1,i) &
	567	- 2.0 * ( hom(k,1,2,sr) + hom(k+1,1,2,sr) ) &
	568	) )**2 &
	569	+ w(k,j,i)**2 )
	570
	571	sums_ll(k,2) = sums_ll(k,2) + 0.5 * w(k,j,i) &
	572	* ( p(k,j,i) + p(k+1,j,i) )
	573
	574	ENDDO
	575	ENDDO
	576	ENDDO
	577	sums_ll(0,1) = 0.0 ! because w is zero at the bottom
	578	sums_ll(nzt+1,1) = 0.0
	579	sums_ll(0,2) = 0.0
	580	sums_ll(nzt+1,2) = 0.0
	581
	582	DO k = nzb_s_outer(j,i)+1, nzt
	583	sums_l(k,55,tn) = ( sums_ll(k,1) - sums_ll(k-1,1) ) * ddzw(k)
	584	sums_l(k,56,tn) = ( sums_ll(k,2) - sums_ll(k-1,2) ) * ddzw(k)
	585	ENDDO
	586	sums_l(nzb,55,tn) = sums_l(nzb+1,55,tn)
	587	sums_l(nzb,56,tn) = sums_l(nzb+1,56,tn)
	588
	589	ENDIF
	590
	591	!
	592	!-- Divergence of vertical flux of SGS TKE
	593	IF ( hom(nzb+1,2,57,0) /= 0.0 ) THEN
	594
	595	!$OMP DO
	596	DO i = nxl, nxr
	597	DO j = nys, nyn
	598	DO k = nzb_s_outer(j,i)+1, nzt
	599
	600	sums_l(k,57,tn) = sums_l(k,57,tn) + ( &
	601	(km(k,j,i)+km(k+1,j,i)) * (e(k+1,j,i)-e(k,j,i)) * ddzu(k+1) &
	602	- (km(k-1,j,i)+km(k,j,i)) * (e(k,j,i)-e(k-1,j,i)) * ddzu(k) &
	603	) * ddzw(k)
	604
	605	ENDDO
	606	ENDDO
	607	ENDDO
	608	sums_l(nzb,57,tn) = sums_l(nzb+1,57,tn)
	609
	610	ENDIF
	611
	612	!
	613	!-- Horizontal heat fluxes (subgrid, resolved, total).
	614	!-- Do it only, if profiles shall be plotted.
	615	IF ( hom(nzb+1,2,58,0) /= 0.0 ) THEN
	616
	617	!$OMP DO
	618	DO i = nxl, nxr
	619	DO j = nys, nyn
	620	DO k = nzb_s_outer(j,i)+1, nzt
	621	!
	622	!-- Subgrid horizontal heat fluxes u"pt", v"pt"
	623	sums_l(k,58,tn) = sums_l(k,58,tn) - 0.5 * &
	624	( kh(k,j,i) + kh(k,j,i-1) ) &
	625	* ( pt(k,j,i-1) - pt(k,j,i) ) &
	626	* ddx * rmask(j,i,sr)
	627	sums_l(k,61,tn) = sums_l(k,61,tn) - 0.5 * &
	628	( kh(k,j,i) + kh(k,j-1,i) ) &
	629	* ( pt(k,j-1,i) - pt(k,j,i) ) &
	630	* ddy * rmask(j,i,sr)
	631	!
	632	!-- Resolved horizontal heat fluxes upt, vpt
	633	sums_l(k,59,tn) = sums_l(k,59,tn) + &
	634	( u(k,j,i) - hom(k,1,1,sr) ) &
	635	* 0.5 * ( pt(k,j,i-1) - hom(k,1,4,sr) + &
	636	pt(k,j,i) - hom(k,1,4,sr) )
	637	pts = 0.5 * ( pt(k,j-1,i) - hom(k,1,4,sr) + &
	638	pt(k,j,i) - hom(k,1,4,sr) )
	639	sums_l(k,62,tn) = sums_l(k,62,tn) + &
	640	( v(k,j,i) - hom(k,1,2,sr) ) &
	641	* 0.5 * ( pt(k,j-1,i) - hom(k,1,4,sr) + &
	642	pt(k,j,i) - hom(k,1,4,sr) )
	643	ENDDO
	644	ENDDO
	645	ENDDO
	646	!
	647	!-- Fluxes at the surface must be zero (e.g. due to the Prandtl-layer)
	648	sums(nzb,58) = 0.0
	649	sums(nzb,59) = 0.0
	650	sums(nzb,60) = 0.0
	651	sums(nzb,61) = 0.0
	652	sums(nzb,62) = 0.0
	653	sums(nzb,63) = 0.0
	654
	655	ENDIF
[87]	656
	657	!
	658	!-- Calculate the user-defined profiles
	659	CALL user_statistics( 'profiles', sr, tn )
[1]	660	!$OMP END PARALLEL
	661
	662	!
	663	!-- Summation of thread sums
	664	IF ( threads_per_task > 1 ) THEN
	665	DO i = 1, threads_per_task-1
	666	sums_l(:,3,0) = sums_l(:,3,0) + sums_l(:,3,i)
	667	sums_l(:,4:40,0) = sums_l(:,4:40,0) + sums_l(:,4:40,i)
[87]	668	sums_l(:,45:pr_palm,0) = sums_l(:,45:pr_palm,0) + &
	669	sums_l(:,45:pr_palm,i)
	670	IF ( max_pr_user > 0 ) THEN
	671	sums_l(:,pr_palm+1:pr_palm+max_pr_user,0) = &
	672	sums_l(:,pr_palm+1:pr_palm+max_pr_user,0) + &
	673	sums_l(:,pr_palm+1:pr_palm+max_pr_user,i)
	674	ENDIF
[1]	675	ENDDO
	676	ENDIF
	677
	678	#if defined( __parallel )
	679	!
	680	!-- Compute total sum from local sums
	681	CALL MPI_ALLREDUCE( sums_l(nzb,1,0), sums(nzb,1), ngp_sums, MPI_REAL, &
	682	MPI_SUM, comm2d, ierr )
	683	#else
	684	sums = sums_l(:,:,0)
	685	#endif
	686
	687	!
	688	!-- Final values are obtained by division by the total number of grid points
	689	!-- used for summation. After that store profiles.
	690	!-- Profiles:
	691	DO k = nzb, nzt+1
[87]	692	sums(k,:pr_palm-2) = sums(k,:pr_palm-2) / ngp_2dh_outer(k,sr)
[1]	693	ENDDO
	694	!-- Upstream-parts
[87]	695	sums(nzb:nzb+11,pr_palm-1) = sums(nzb:nzb+11,pr_palm-1) / ngp_3d(sr)
[1]	696	!-- u* and so on
[87]	697	!-- As sums(nzb:nzb+3,pr_palm) are full 2D arrays (us, usws, vsws, ts) whose
[1]	698	!-- size is always ( nx + 1 ) * ( ny + 1 ), defined at the first grid layer
	699	!-- above the topography, they are being divided by ngp_2dh(sr)
[87]	700	sums(nzb:nzb+3,pr_palm) = sums(nzb:nzb+3,pr_palm) / &
[1]	701	ngp_2dh(sr)
	702	!-- eges, e*
[87]	703	sums(nzb+4:nzb+5,pr_palm) = sums(nzb+4:nzb+5,pr_palm) / &
[1]	704	ngp_3d_inner(sr)
	705	!-- Old and new divergence
[87]	706	sums(nzb+9:nzb+10,pr_palm) = sums(nzb+9:nzb+10,pr_palm) / &
[1]	707	ngp_3d_inner(sr)
	708
[87]	709	!-- User-defined profiles
	710	IF ( max_pr_user > 0 ) THEN
	711	DO k = nzb, nzt+1
	712	sums(k,pr_palm+1:pr_palm+max_pr_user) = &
	713	sums(k,pr_palm+1:pr_palm+max_pr_user) / &
	714	ngp_2dh_outer(k,sr)
	715	ENDDO
	716	ENDIF
	717
[1]	718	!
	719	!-- Collect horizontal average in hom.
	720	!-- Compute deduced averages (e.g. total heat flux)
	721	hom(:,1,3,sr) = sums(:,3) ! w
	722	hom(:,1,8,sr) = sums(:,8) ! e profiles 5-7 are initial profiles
	723	hom(:,1,9,sr) = sums(:,9) ! km
	724	hom(:,1,10,sr) = sums(:,10) ! kh
	725	hom(:,1,11,sr) = sums(:,11) ! l
	726	hom(:,1,12,sr) = sums(:,12) ! w"u"
	727	hom(:,1,13,sr) = sums(:,13) ! wu
	728	hom(:,1,14,sr) = sums(:,14) ! w"v"
	729	hom(:,1,15,sr) = sums(:,15) ! wv
	730	hom(:,1,16,sr) = sums(:,16) ! w"pt"
	731	hom(:,1,17,sr) = sums(:,17) ! wpt
	732	hom(:,1,18,sr) = sums(:,16) + sums(:,17) ! wpt
	733	hom(:,1,19,sr) = sums(:,12) + sums(:,13) ! wu
	734	hom(:,1,20,sr) = sums(:,14) + sums(:,15) ! wv
	735	hom(:,1,21,sr) = sums(:,21) ! wptBC
	736	hom(:,1,22,sr) = sums(:,16) + sums(:,21) ! wptBC
	737	! profiles 23-29 left empty for initial
	738	! profiles
	739	hom(:,1,30,sr) = sums(:,30) ! u*2
	740	hom(:,1,31,sr) = sums(:,31) ! v*2
	741	hom(:,1,32,sr) = sums(:,32) ! w*2
	742	hom(:,1,33,sr) = sums(:,33) ! pt*2
	743	hom(:,1,34,sr) = sums(:,34) ! e*
	744	hom(:,1,35,sr) = sums(:,35) ! w2pt
	745	hom(:,1,36,sr) = sums(:,36) ! wpt2
	746	hom(:,1,37,sr) = sums(:,37) ! we
	747	hom(:,1,38,sr) = sums(:,38) ! w*3
	748	hom(:,1,39,sr) = sums(:,38) / ( sums(:,32) + 1E-20 )**1.5 ! Sw
	749	hom(:,1,40,sr) = sums(:,40) ! p
	750	hom(:,1,45,sr) = sums(:,45) ! w"q"
	751	hom(:,1,46,sr) = sums(:,46) ! wvpt
	752	hom(:,1,47,sr) = sums(:,45) + sums(:,46) ! wvpt
	753	hom(:,1,48,sr) = sums(:,48) ! w"q" (w"qv")
	754	hom(:,1,49,sr) = sums(:,49) ! wq (wqv)
	755	hom(:,1,50,sr) = sums(:,48) + sums(:,49) ! wq (wqv)
	756	hom(:,1,51,sr) = sums(:,51) ! w"qv"
	757	hom(:,1,52,sr) = sums(:,52) ! wqv
	758	hom(:,1,53,sr) = sums(:,52) + sums(:,51) ! wq (wqv)
	759	hom(:,1,54,sr) = sums(:,54) ! ql
	760	hom(:,1,55,sr) = sums(:,55) ! wuu*/dz
	761	hom(:,1,56,sr) = sums(:,56) ! wp/dz
	762	hom(:,1,57,sr) = sums(:,57) ! w"e/dz
	763	hom(:,1,58,sr) = sums(:,58) ! u"pt"
	764	hom(:,1,59,sr) = sums(:,59) ! upt
	765	hom(:,1,60,sr) = sums(:,58) + sums(:,59) ! upt_t
	766	hom(:,1,61,sr) = sums(:,61) ! v"pt"
	767	hom(:,1,62,sr) = sums(:,62) ! vpt
	768	hom(:,1,63,sr) = sums(:,61) + sums(:,62) ! vpt_t
	769
[87]	770	hom(:,1,pr_palm-1,sr) = sums(:,pr_palm-1)
[1]	771	! upstream-parts u_x, u_y, u_z, v_x,
	772	! v_y, usw. (in last but one profile)
[87]	773	hom(:,1,pr_palm,sr) = sums(:,pr_palm)
[1]	774	! u, w'u', w'v', t (in last profile)
	775
[87]	776	IF ( max_pr_user > 0 ) THEN ! user-defined profiles
	777	hom(:,1,pr_palm+1:pr_palm+max_pr_user,sr) = &
	778	sums(:,pr_palm+1:pr_palm+max_pr_user)
	779	ENDIF
	780
[1]	781	!
	782	!-- Determine the boundary layer height using two different schemes.
	783	!-- First scheme: Starting from the Earth's surface, look for the first
	784	!-- relative minimum of the total heat flux. The corresponding height is
	785	!-- accepted as the boundary layer height, if it is less than 1.5 times the
	786	!-- height where the heat flux becomes negative for the first time.
	787	!-- NOTE: This criterion is still capable of improving!
	788	z_i(1) = 0.0
	789	first = .TRUE.
	790	DO k = nzb, nzt-1
	791	IF ( first .AND. hom(k,1,18,sr) < 0.0 ) THEN
	792	first = .FALSE.
	793	height = zw(k)
	794	ENDIF
	795	IF ( hom(k,1,18,sr) < 0.0 .AND. &
	796	hom(k+1,1,18,sr) > hom(k,1,18,sr) ) THEN
	797	IF ( zw(k) < 1.5 * height ) THEN
	798	z_i(1) = zw(k)
	799	ELSE
	800	z_i(1) = height
	801	ENDIF
	802	EXIT
	803	ENDIF
	804	ENDDO
	805
	806	!
	807	!-- Second scheme: Starting from the top model boundary, look for the first
	808	!-- characteristic kink in the temperature profile, where the originally
	809	!-- stable stratification notably weakens.
	810	z_i(2) = 0.0
	811	DO k = nzt-1, nzb+1, -1
	812	IF ( ( hom(k+1,1,4,sr) - hom(k,1,4,sr) ) > &
	813	2.0 * ( hom(k,1,4,sr) - hom(k-1,1,4,sr) ) ) THEN
	814	z_i(2) = zu(k)
	815	EXIT
	816	ENDIF
	817	ENDDO
	818
[87]	819	hom(nzb+6,1,pr_palm,sr) = z_i(1)
	820	hom(nzb+7,1,pr_palm,sr) = z_i(2)
[1]	821
	822	!
	823	!-- Computation of both the characteristic vertical velocity and
	824	!-- the characteristic convective boundary layer temperature.
	825	!-- The horizontal average at nzb+1 is input for the average temperature.
	826	IF ( hom(nzb,1,18,sr) > 0.0 .AND. z_i(1) /= 0.0 ) THEN
[87]	827	hom(nzb+8,1,pr_palm,sr) = ( g / hom(nzb+1,1,4,sr) * &
[1]	828	hom(nzb,1,18,sr) * z_i(1) )**0.333333333
	829	!-- so far this only works if Prandtl layer is used
[87]	830	hom(nzb+11,1,pr_palm,sr) = hom(nzb,1,16,sr) / hom(nzb+8,1,pr_palm,sr)
[1]	831	ELSE
[87]	832	hom(nzb+8,1,pr_palm,sr) = 0.0
	833	hom(nzb+11,1,pr_palm,sr) = 0.0
[1]	834	ENDIF
	835
[48]	836	!
	837	!-- Collect the time series quantities
[87]	838	ts_value(1,sr) = hom(nzb+4,1,pr_palm,sr) ! E
	839	ts_value(2,sr) = hom(nzb+5,1,pr_palm,sr) ! E*
[48]	840	ts_value(3,sr) = dt_3d
[87]	841	ts_value(4,sr) = hom(nzb,1,pr_palm,sr) ! u*
	842	ts_value(5,sr) = hom(nzb+3,1,pr_palm,sr) ! th*
[48]	843	ts_value(6,sr) = u_max
	844	ts_value(7,sr) = v_max
	845	ts_value(8,sr) = w_max
[87]	846	ts_value(9,sr) = hom(nzb+10,1,pr_palm,sr) ! new divergence
	847	ts_value(10,sr) = hom(nzb+9,1,pr_palm,sr) ! old Divergence
	848	ts_value(11,sr) = hom(nzb+6,1,pr_palm,sr) ! z_i(1)
	849	ts_value(12,sr) = hom(nzb+7,1,pr_palm,sr) ! z_i(2)
	850	ts_value(13,sr) = hom(nzb+8,1,pr_palm,sr) ! w*
[48]	851	ts_value(14,sr) = hom(nzb,1,16,sr) ! w'pt' at k=0
	852	ts_value(15,sr) = hom(nzb+1,1,16,sr) ! w'pt' at k=1
	853	ts_value(16,sr) = hom(nzb+1,1,18,sr) ! wpt at k=1
	854	ts_value(17,sr) = hom(nzb,1,4,sr) ! pt(0)
	855	ts_value(18,sr) = hom(nzb+1,1,4,sr) ! pt(zp)
[87]	856	ts_value(19,sr) = hom(nzb+9,1,pr_palm-1,sr) ! splptx
	857	ts_value(20,sr) = hom(nzb+10,1,pr_palm-1,sr) ! splpty
	858	ts_value(21,sr) = hom(nzb+11,1,pr_palm-1,sr) ! splptz
[48]	859	IF ( ts_value(5,sr) /= 0.0 ) THEN
	860	ts_value(22,sr) = ts_value(4,sr)**2 / &
	861	( kappa * g * ts_value(5,sr) / ts_value(18,sr) ) ! L
	862	ELSE
	863	ts_value(22,sr) = 10000.0
	864	ENDIF
[1]	865
	866	!
[48]	867	!-- Calculate additional statistics provided by the user interface
[87]	868	CALL user_statistics( 'time_series', sr, 0 )
[1]	869
[48]	870	ENDDO ! loop of the subregions
	871
[1]	872	!
	873	!-- If required, sum up horizontal averages for subsequent time averaging
	874	IF ( do_sum ) THEN
	875	IF ( average_count_pr == 0 ) hom_sum = 0.0
	876	hom_sum = hom_sum + hom(:,1,:,:)
	877	average_count_pr = average_count_pr + 1
	878	do_sum = .FALSE.
	879	ENDIF
	880
	881	!
	882	!-- Set flag for other UPs (e.g. output routines, but also buoyancy).
	883	!-- This flag is reset after each time step in time_integration.
	884	flow_statistics_called = .TRUE.
	885
	886	CALL cpu_log( log_point(10), 'flow_statistics', 'stop' )
	887
	888
	889	END SUBROUTINE flow_statistics
	890
	891
	892

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |