Home

Context Navigation

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

Last change on this file since 50 was 48, checked in by raasch, 18 years ago
preliminary version, several changes to be explained later
Property svn:keywords set to `Id`
File size: 36.6 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |