Home

Context Navigation

flow_statistics.f90 @ 86

Last change on this file since 86 was 83, checked in by raasch, 18 years ago

New:
---

Changed:

PALM can be generally installed on any kind of Linux-, IBM-AIX-, or NEC-SX-system by adding appropriate settings to the configuration file.

Scripts are also running under the public domain ksh.

All system relevant compile and link options as well as the host identifier (local_host) are specified in the configuration file.

Filetransfer by ftp removed (options -f removed from mrun and mbuild).

Call of (system-)FLUSH routine moved to new routine local_flush.

return_addres and return_username are read from ENVPAR-NAMELIST-file instead of using local_getenv.

Preprocessor strings for different linux clusters changed to "lc", some preprocessor directives renamed (new: intel_openmp_bug), preprocessor directives for old systems removed

advec_particles, check_open, cpu_log, cpu_statistics, data_output_dvrp, flow_statistics, header, init_dvrp, init_particles, init_1d_model, init_dvrp, init_pegrid, local_getenv, local_system, local_tremain, local_tremain_ini, modules, palm, parin, run_control

new:
local_flush

mbuild, mrun

Errors:

Property svn:keywords set to Id

File size: 36.9 KB

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

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

Context Navigation

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

Download in other formats: