Home

Context Navigation

source: palm/trunk/SOURCE/calc_spectra.f90 @ 1336

Last change on this file since 1336 was 1325, checked in by suehring, 11 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 14.8 KB

Rev	Line
[1]	1	SUBROUTINE calc_spectra
	2
[1036]	3	!--------------------------------------------------------------------------------!
	4	! This file is part of PALM.
	5	!
	6	! PALM is free software: you can redistribute it and/or modify it under the terms
	7	! of the GNU General Public License as published by the Free Software Foundation,
	8	! either version 3 of the License, or (at your option) any later version.
	9	!
	10	! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
	11	! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
	12	! A PARTICULAR PURPOSE. See the GNU General Public License for more details.
	13	!
	14	! You should have received a copy of the GNU General Public License along with
	15	! PALM. If not, see <http://www.gnu.org/licenses/>.
	16	!
[1310]	17	! Copyright 1997-2014 Leibniz Universitaet Hannover
[1036]	18	!--------------------------------------------------------------------------------!
	19	!
[247]	20	! Current revisions:
[1]	21	! -----------------
[1321]	22	!
	23	! Former revisions:
	24	! -----------------
	25	! $Id: calc_spectra.f90 1325 2014-03-21 09:21:15Z raasch $
	26	!
[1325]	27	! 1324 2014-03-21 09:13:16Z suehring
	28	! Bugfix: nzb_x, nzb_yd, nyn_x, nyn_x, nzt_x, nzt_yd belong to transpose_indices
	29	!
[1321]	30	! 1320 2014-03-20 08:40:49Z raasch
[1320]	31	! ONLY-attribute added to USE-statements,
	32	! kind-parameters added to all INTEGER and REAL declaration statements,
	33	! kinds are defined in new module kinds,
	34	! revision history before 2012 removed,
	35	! comment fields (!:) to be used for variable explanations added to
	36	! all variable declaration statements
[198]	37	!
[1319]	38	! 1318 2014-03-17 13:35:16Z raasch
	39	! module interfaces removed
	40	!
[1217]	41	! 1216 2013-08-26 09:31:42Z raasch
	42	! resorting of array moved to separate routine resort_for_zx,
	43	! one argument removed from the transpose_..d routines
	44	!
[1121]	45	! 1120 2013-04-05 15:11:35Z raasch
	46	! bugfix: calls of fft_x\|y replaced by fft_x\|y_1d
	47	!
[1037]	48	! 1036 2012-10-22 13:43:42Z raasch
	49	! code put under GPL (PALM 3.9)
	50	!
[1004]	51	! 1003 2012-09-14 14:35:53Z raasch
	52	! adjustment of array tend for cases with unequal subdomain sizes removed
	53	!
[1]	54	! Revision 1.1 2001/01/05 15:08:07 raasch
	55	! Initial revision
	56	!
	57	!
	58	! Description:
	59	! ------------
	60	! Calculate horizontal spectra along x and y.
	61	! ATTENTION: 1d-decomposition along y still needs improvement, because in that
	62	! case the gridpoint number along z still depends on the PE number
	63	! because transpose_xz has to be used (and possibly also
	64	! transpose_zyd needs modification).
	65	!------------------------------------------------------------------------------!
	66
	67	#if defined( __spectra )
[1320]	68	USE arrays_3d, &
	69	ONLY: d, tend
	70
	71	USE control_parameters, &
	72	ONLY: average_count_sp, bc_lr_cyc, bc_ns_cyc, message_string, psolver
	73
	74	USE cpulog, &
	75	ONLY: cpu_log, log_point
	76
	77	USE fft_xy, &
	78	ONLY: fft_init
	79
	80	USE indices, &
[1324]	81	ONLY: nxl, nxr, nyn, nys, nzb, nzt
[1320]	82
	83	USE kinds
	84
[1]	85	USE pegrid
	86
[1320]	87	USE spectrum, &
	88	ONLY: data_output_sp, spectra_direction
	89
	90
[1]	91	IMPLICIT NONE
	92
[1320]	93	INTEGER(iwp) :: m !:
	94	INTEGER(iwp) :: pr !:
[1]	95
	96
	97	CALL cpu_log( log_point(30), 'calc_spectra', 'start' )
	98
	99	!
	100	!-- Initialize ffts
	101	CALL fft_init
	102
	103	!
[225]	104	!-- Reallocate array d in required size
	105	IF ( psolver == 'multigrid' ) THEN
	106	DEALLOCATE( d )
[1003]	107	ALLOCATE( d(nzb+1:nzt,nys:nyn,nxl:nxr) )
[225]	108	ENDIF
	109
[1]	110	m = 1
	111	DO WHILE ( data_output_sp(m) /= ' ' .AND. m <= 10 )
	112	!
	113	!-- Transposition from z --> x ( y --> x in case of a 1d-decomposition
	114	!-- along x)
	115	IF ( INDEX( spectra_direction(m), 'x' ) /= 0 ) THEN
	116
	117	!
	118	!-- Calculation of spectra works for cyclic boundary conditions only
[707]	119	IF ( .NOT. bc_lr_cyc ) THEN
[247]	120
[274]	121	message_string = 'non-cyclic lateral boundaries along x do not'// &
[247]	122	'& allow calculation of spectra along x'
	123	CALL message( 'calc_spectra', 'PA0160', 1, 2, 0, 6, 0 )
[1]	124	ENDIF
	125
	126	CALL preprocess_spectra( m, pr )
	127
	128	#if defined( __parallel )
	129	IF ( pdims(2) /= 1 ) THEN
[1216]	130	CALL resort_for_zx( d, tend )
	131	CALL transpose_zx( tend, d )
[1]	132	ELSE
[1216]	133	CALL transpose_yxd( d, d )
[1]	134	ENDIF
	135	CALL calc_spectra_x( d, pr, m )
	136	#else
[274]	137	message_string = 'sorry, calculation of spectra in non parallel ' // &
	138	'mode& is still not realized'
[247]	139	CALL message( 'calc_spectra', 'PA0161', 1, 2, 0, 6, 0 )
[1]	140	#endif
	141
	142	ENDIF
	143
	144	!
	145	!-- Transposition from z --> y (d is rearranged only in case of a
	146	!-- 1d-decomposition along x)
	147	IF ( INDEX( spectra_direction(m), 'y' ) /= 0 ) THEN
	148
	149	!
	150	!-- Calculation of spectra works for cyclic boundary conditions only
[707]	151	IF ( .NOT. bc_ns_cyc ) THEN
[1]	152	IF ( myid == 0 ) THEN
[274]	153	message_string = 'non-cyclic lateral boundaries along y do' // &
	154	' not & allow calculation of spectra along y'
[247]	155	CALL message( 'calc_spectra', 'PA0162', 1, 2, 0, 6, 0 )
[1]	156	ENDIF
	157	CALL local_stop
	158	ENDIF
	159
	160	CALL preprocess_spectra( m, pr )
	161
	162	#if defined( __parallel )
[1216]	163	CALL transpose_zyd( d, d )
[1]	164	CALL calc_spectra_y( d, pr, m )
	165	#else
[1320]	166	message_string = 'sorry, calculation of spectra in non parallel' // &
[274]	167	'mode& is still not realized'
[247]	168	CALL message( 'calc_spectra', 'PA0161', 1, 2, 0, 6, 0 )
[1]	169	#endif
	170
	171	ENDIF
	172
	173	!
	174	!-- Increase counter for next spectrum
	175	m = m + 1
	176
	177	ENDDO
	178
	179	!
	180	!-- Increase counter for averaging process in routine plot_spectra
	181	average_count_sp = average_count_sp + 1
	182
	183	CALL cpu_log( log_point(30), 'calc_spectra', 'stop' )
	184
	185	#endif
	186	END SUBROUTINE calc_spectra
	187
	188
	189	#if defined( __spectra )
	190	SUBROUTINE preprocess_spectra( m, pr )
	191
[1320]	192	USE arrays_3d, &
	193	ONLY: d, pt, q, u, v, w
	194
	195	USE indices, &
	196	ONLY: nxl, nxr, nyn, nys, nzb, nzt
	197
	198	USE kinds
	199
[1]	200	USE pegrid
	201
[1320]	202	USE spectrum, &
	203	ONLY: data_output_sp
	204
	205	USE statistics, &
	206	ONLY: sums
	207
	208
[1]	209	IMPLICIT NONE
	210
[1320]	211	INTEGER(iwp) :: i !:
	212	INTEGER(iwp) :: j !:
	213	INTEGER(iwp) :: k !:
	214	INTEGER(iwp) :: m !:
	215	INTEGER(iwp) :: pr !:
[1]	216
	217	SELECT CASE ( TRIM( data_output_sp(m) ) )
	218
	219	CASE ( 'u' )
	220	pr = 1
	221	d(nzb+1:nzt,nys:nyn,nxl:nxr) = u(nzb+1:nzt,nys:nyn,nxl:nxr)
	222
	223	CASE ( 'v' )
	224	pr = 2
	225	d(nzb+1:nzt,nys:nyn,nxl:nxr) = v(nzb+1:nzt,nys:nyn,nxl:nxr)
	226
	227	CASE ( 'w' )
	228	pr = 3
	229	d(nzb+1:nzt,nys:nyn,nxl:nxr) = w(nzb+1:nzt,nys:nyn,nxl:nxr)
	230
	231	CASE ( 'pt' )
	232	pr = 4
	233	d(nzb+1:nzt,nys:nyn,nxl:nxr) = pt(nzb+1:nzt,nys:nyn,nxl:nxr)
	234
	235	CASE ( 'q' )
	236	pr = 41
	237	d(nzb+1:nzt,nys:nyn,nxl:nxr) = q(nzb+1:nzt,nys:nyn,nxl:nxr)
	238
	239	CASE DEFAULT
[144]	240	!
	241	!-- The DEFAULT case is reached either if the parameter data_output_sp(m)
	242	!-- contains a wrong character string or if the user has coded a special
	243	!-- case in the user interface. There, the subroutine user_spectra
	244	!-- checks which of these two conditions applies.
	245	CALL user_spectra( 'preprocess', m, pr )
[1]	246
	247	END SELECT
	248
	249	!
	250	!-- Subtract horizontal mean from the array, for which spectra have to be
	251	!-- calculated
	252	DO i = nxl, nxr
	253	DO j = nys, nyn
	254	DO k = nzb+1, nzt
	255	d(k,j,i) = d(k,j,i) - sums(k,pr)
	256	ENDDO
	257	ENDDO
	258	ENDDO
	259
	260	END SUBROUTINE preprocess_spectra
	261
	262
	263	SUBROUTINE calc_spectra_x( ddd, pr, m )
	264
[1320]	265	USE arrays_3d, &
	266	ONLY:
	267
	268	USE control_parameters, &
	269	ONLY: fft_method
	270
	271	USE fft_xy, &
	272	ONLY: fft_x_1d
	273
	274	USE grid_variables, &
	275	ONLY: dx
	276
	277	USE indices, &
[1324]	278	ONLY: nx, ny
[1320]	279
	280	USE kinds
	281
[1]	282	USE pegrid
[1320]	283
	284	USE spectrum, &
	285	ONLY: comp_spectra_level, n_sp_x
	286
	287	USE statistics, &
	288	ONLY: spectrum_x
	289
[1324]	290	USE transpose_indices, &
	291	ONLY: nyn_x, nys_x, nzb_x, nzt_x
[1]	292
[1320]	293
[1]	294	IMPLICIT NONE
	295
[1320]	296	INTEGER(iwp) :: i !:
	297	INTEGER(iwp) :: ishape(1) !:
	298	INTEGER(iwp) :: j !:
	299	INTEGER(iwp) :: k !:
	300	INTEGER(iwp) :: m !:
	301	INTEGER(iwp) :: n !:
	302	INTEGER(iwp) :: pr !:
[1]	303
[1320]	304	REAL(wp) :: fac !:
	305	REAL(wp) :: exponent !:
	306
	307	REAL(wp), DIMENSION(0:nx) :: work !:
	308
	309	REAL(wp), DIMENSION(0:nx/2) :: sums_spectra_l !:
	310
	311	REAL(wp), DIMENSION(0:nx/2,100) :: sums_spectra !:
	312
	313	REAL(wp), DIMENSION(0:nx,nys_x:nyn_x,nzb_x:nzt_x) :: ddd !:
[1]	314
	315	!
	316	!-- Exponent for geometric average
[192]	317	exponent = 1.0 / ( ny + 1.0 )
[1]	318
	319	!
	320	!-- Loop over all levels defined by the user
	321	n = 1
[189]	322	DO WHILE ( comp_spectra_level(n) /= 999999 .AND. n <= 100 )
[1]	323
	324	k = comp_spectra_level(n)
	325
	326	!
	327	!-- Calculate FFT only if the corresponding level is situated on this PE
	328	IF ( k >= nzb_x .AND. k <= nzt_x ) THEN
	329
	330	DO j = nys_x, nyn_x
	331
	332	work = ddd(0:nx,j,k)
[1120]	333	CALL fft_x_1d( work, 'forward' )
[1]	334
	335	ddd(0,j,k) = dx * work(0)**2
	336	DO i = 1, nx/2
	337	ddd(i,j,k) = dx * ( work(i)2 + work(nx+1-i)2 )
	338	ENDDO
	339
	340	ENDDO
	341
	342	!
	343	!-- Local sum and geometric average of these spectra
	344	!-- (WARNING: no global sum should be performed, because floating
	345	!-- point overflow may occur)
	346	DO i = 0, nx/2
	347
	348	sums_spectra_l(i) = 1.0
	349	DO j = nys_x, nyn_x
	350	sums_spectra_l(i) = sums_spectra_l(i) * ddd(i,j,k)**exponent
	351	ENDDO
	352
	353	ENDDO
	354
	355	ELSE
	356
	357	sums_spectra_l = 1.0
	358
	359	ENDIF
	360
	361	!
	362	!-- Global sum of spectra on PE0 (from where they are written on file)
	363	sums_spectra(:,n) = 0.0
	364	#if defined( __parallel )
	365	CALL MPI_BARRIER( comm2d, ierr ) ! Necessary?
[1320]	366	CALL MPI_REDUCE( sums_spectra_l(0), sums_spectra(0,n), nx/2+1, &
[1]	367	MPI_REAL, MPI_PROD, 0, comm2d, ierr )
	368	#else
	369	sums_spectra(:,n) = sums_spectra_l
	370	#endif
	371
	372	n = n + 1
	373
	374	ENDDO
	375	n = n - 1
	376
	377	IF ( myid == 0 ) THEN
	378	!
[146]	379	!-- Sum of spectra for later averaging (see routine data_output_spectra)
[1]	380	!-- Temperton fft results need to be normalized
	381	IF ( fft_method == 'temperton-algorithm' ) THEN
	382	fac = nx + 1.0
	383	ELSE
	384	fac = 1.0
	385	ENDIF
	386	DO i = 1, nx/2
	387	DO k = 1, n
	388	spectrum_x(i,k,m) = spectrum_x(i,k,m) + sums_spectra(i,k) * fac
	389	ENDDO
	390	ENDDO
	391
	392	ENDIF
	393
	394	!
[146]	395	!-- n_sp_x is needed by data_output_spectra_x
[1]	396	n_sp_x = n
	397
	398	END SUBROUTINE calc_spectra_x
	399
	400
	401	SUBROUTINE calc_spectra_y( ddd, pr, m )
	402
[1320]	403	USE arrays_3d, &
	404	ONLY:
	405
	406	USE control_parameters, &
	407	ONLY: fft_method
	408
	409	USE fft_xy, &
	410	ONLY: fft_y_1d
	411
	412	USE grid_variables, &
	413	ONLY: dy
	414
	415	USE indices, &
[1324]	416	ONLY: nx, ny
[1320]	417
	418	USE kinds
	419
[1]	420	USE pegrid
[1320]	421
	422	USE spectrum, &
	423	ONLY: comp_spectra_level, n_sp_y
	424
	425	USE statistics, &
	426	ONLY: spectrum_y
	427
[1324]	428	USE transpose_indices, &
	429	ONLY: nxl_yd, nxr_yd, nzb_yd, nzt_yd
[1]	430
[1320]	431
[1]	432	IMPLICIT NONE
	433
[1320]	434	INTEGER(iwp) :: i !:
	435	INTEGER(iwp) :: j !:
	436	INTEGER(iwp) :: jshape(1) !:
	437	INTEGER(iwp) :: k !:
	438	INTEGER(iwp) :: m !:
	439	INTEGER(iwp) :: n !:
	440	INTEGER(iwp) :: pr !:
[1]	441
[1320]	442	REAL(wp) :: fac !:
	443	REAL(wp) :: exponent !:
	444
	445	REAL(wp), DIMENSION(0:ny) :: work !:
	446
	447	REAL(wp), DIMENSION(0:ny/2) :: sums_spectra_l !:
	448
	449	REAL(wp), DIMENSION(0:ny/2,100) :: sums_spectra !:
	450
	451	REAL(wp), DIMENSION(0:ny,nxl_yd:nxr_yd,nzb_yd:nzt_yd) :: ddd !:
[1]	452
	453
	454	!
	455	!-- Exponent for geometric average
	456	exponent = 1.0 / ( nx + 1.0 )
	457
	458	!
	459	!-- Loop over all levels defined by the user
	460	n = 1
[189]	461	DO WHILE ( comp_spectra_level(n) /= 999999 .AND. n <= 100 )
[1]	462
	463	k = comp_spectra_level(n)
	464
	465	!
	466	!-- Calculate FFT only if the corresponding level is situated on this PE
	467	IF ( k >= nzb_yd .AND. k <= nzt_yd ) THEN
	468
	469	DO i = nxl_yd, nxr_yd
	470
	471	work = ddd(0:ny,i,k)
[1120]	472	CALL fft_y_1d( work, 'forward' )
[1]	473
	474	ddd(0,i,k) = dy * work(0)**2
	475	DO j = 1, ny/2
	476	ddd(j,i,k) = dy * ( work(j)2 + work(ny+1-j)2 )
	477	ENDDO
	478
	479	ENDDO
	480
	481	!
	482	!-- Local sum and geometric average of these spectra
	483	!-- (WARNING: no global sum should be performed, because floating
	484	!-- point overflow may occur)
	485	DO j = 0, ny/2
	486
	487	sums_spectra_l(j) = 1.0
	488	DO i = nxl_yd, nxr_yd
	489	sums_spectra_l(j) = sums_spectra_l(j) * ddd(j,i,k)**exponent
	490	ENDDO
	491
	492	ENDDO
	493
	494	ELSE
	495
	496	sums_spectra_l = 1.0
	497
	498	ENDIF
	499
	500	!
	501	!-- Global sum of spectra on PE0 (from where they are written on file)
	502	sums_spectra(:,n) = 0.0
	503	#if defined( __parallel )
	504	CALL MPI_BARRIER( comm2d, ierr ) ! Necessary?
[1320]	505	CALL MPI_REDUCE( sums_spectra_l(0), sums_spectra(0,n), ny/2+1, &
[1]	506	MPI_REAL, MPI_PROD, 0, comm2d, ierr )
	507	#else
	508	sums_spectra(:,n) = sums_spectra_l
	509	#endif
	510
	511	n = n + 1
	512
	513	ENDDO
	514	n = n - 1
	515
	516
	517	IF ( myid == 0 ) THEN
	518	!
[146]	519	!-- Sum of spectra for later averaging (see routine data_output_spectra)
[1]	520	!-- Temperton fft results need to be normalized
	521	IF ( fft_method == 'temperton-algorithm' ) THEN
	522	fac = ny + 1.0
	523	ELSE
	524	fac = 1.0
	525	ENDIF
	526	DO j = 1, ny/2
	527	DO k = 1, n
	528	spectrum_y(j,k,m) = spectrum_y(j,k,m) + sums_spectra(j,k) * fac
	529	ENDDO
	530	ENDDO
	531
	532	ENDIF
	533
	534	!
[146]	535	!-- n_sp_y is needed by data_output_spectra_y
[1]	536	n_sp_y = n
	537
	538	END SUBROUTINE calc_spectra_y
	539	#endif

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |