Home

Context Navigation

spectra_mod.f90 @ 4847

Last change on this file since 4847 was 4843, checked in by raasch, 4 years ago
local namelist parameter added to switch off the module although the respective module namelist appears in the namelist file, further copyright updates
Property svn:keywords set to `Id`
File size: 31.7 KB

Rev	Line
[1833]	1	!> @file spectra_mod.f90
[4591]	2	!--------------------------------------------------------------------------------------------------!
[2696]	3	! This file is part of the PALM model system.
[1036]	4	!
[4591]	5	! PALM is free software: you can redistribute it and/or modify it under the terms of the GNU General
	6	! Public License as published by the Free Software Foundation, either version 3 of the License, or
	7	! (at your option) any later version.
[1036]	8	!
[4591]	9	! PALM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the
	10	! implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
	11	! Public License for more details.
[1036]	12	!
[4591]	13	! You should have received a copy of the GNU General Public License along with PALM. If not, see
	14	! <http://www.gnu.org/licenses/>.
[1036]	15	!
[4828]	16	! Copyright 1997-2021 Leibniz Universitaet Hannover
[4591]	17	!--------------------------------------------------------------------------------------------------!
[1036]	18	!
[4591]	19	!
[247]	20	! Current revisions:
[1]	21	! -----------------
[4591]	22	!
	23	!
[1787]	24	! Former revisions:
	25	! -----------------
	26	! $Id: spectra_mod.f90 4843 2021-01-15 15:22:11Z raasch $
[4843]	27	! local namelist parameter added to switch off the module although the respective module namelist
	28	! appears in the namelist file
	29	!
	30	! 4842 2021-01-14 10:42:28Z raasch
[4842]	31	! reading of namelist file and actions in case of namelist errors revised so that statement labels
	32	! and goto statements are not required any more,
	33	! deprecated namelist removed
	34	!
	35	! 4828 2021-01-05 11:21:41Z Giersch
[4629]	36	! support for MPI Fortran77 interface (mpif.h) removed
	37	!
	38	! 4591 2020-07-06 15:56:08Z raasch
[4591]	39	! File re-formatted to follow the PALM coding standard
	40	!
	41	! 4429 2020-02-27 15:24:30Z raasch
	42	! Bugfix: preprocessor directives rearranged for serial mode
	43	!
[4429]	44	! 4360 2020-01-07 11:25:50Z suehring
[4182]	45	! Corrected "Former revisions" section
[4591]	46	!
[4182]	47	! 3805 2019-03-20 15:26:35Z raasch
[4591]	48	! Unused variable removed
	49	!
[3805]	50	! 3655 2019-01-07 16:51:22Z knoop
[3421]	51	! Renamed output variables
[1787]	52	!
[4182]	53	! Revision 1.1 2001/01/05 15:08:07 raasch
	54	! Initial revision
	55	!
[4591]	56	!--------------------------------------------------------------------------------------------------!
[1]	57	! Description:
	58	! ------------
[4591]	59	!> Calculate horizontal spectra along x and y.
	60	!> ATTENTION: 1d-decomposition along y still needs improvement, because in that case the gridpoint
	61	!> number along z still depends on the PE number because transpose_xz has to be used
	62	!> (and possibly also transpose_zyd needs modification).
	63	!--------------------------------------------------------------------------------------------------!
[1833]	64	MODULE spectra_mod
[1]	65
[1786]	66	USE kinds
[1320]	67
[1786]	68	PRIVATE
[1320]	69
[4591]	70	CHARACTER (LEN=10), DIMENSION(10) :: data_output_sp = ' ' !<
	71	CHARACTER (LEN=2), DIMENSION(10) :: spectra_direction = 'x' !<
[1320]	72
[4591]	73	INTEGER(iwp) :: average_count_sp = 0 !<
	74	INTEGER(iwp) :: comp_spectra_level(100) = 999999 !<
	75	INTEGER(iwp) :: dosp_time_count = 0 !<
	76	INTEGER(iwp) :: n_sp_x = 0, n_sp_y = 0 !<
[1320]	77
[1833]	78	LOGICAL :: calculate_spectra = .FALSE. !< internal switch that spectra are calculated
	79	LOGICAL :: spectra_initialized = .FALSE. !< internal switch that spectra related quantities are initialized
[1320]	80
[1833]	81	REAL(wp) :: averaging_interval_sp = 9999999.9_wp !< averaging interval for spectra output
	82	REAL(wp) :: dt_dosp = 9999999.9_wp !< time interval for spectra output
	83	REAL(wp) :: skip_time_dosp = 9999999.9_wp !< no output of spectra data before this interval has passed
[1]	84
[4591]	85	REAL(wp), DIMENSION(:), ALLOCATABLE :: var_d !<
[1833]	86
[4591]	87	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: spectrum_x !<
	88	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: spectrum_y !<
[1833]	89
[1786]	90	SAVE
[1320]	91
[1786]	92	INTERFACE calc_spectra
	93	MODULE PROCEDURE calc_spectra
	94	END INTERFACE calc_spectra
[1320]	95
[1786]	96	INTERFACE preprocess_spectra
	97	MODULE PROCEDURE preprocess_spectra
	98	END INTERFACE preprocess_spectra
[1]	99
[4429]	100	#if defined( __parallel )
[1786]	101	INTERFACE calc_spectra_x
	102	MODULE PROCEDURE calc_spectra_x
	103	END INTERFACE calc_spectra_x
[1]	104
[1786]	105	INTERFACE calc_spectra_y
	106	MODULE PROCEDURE calc_spectra_y
	107	END INTERFACE calc_spectra_y
[4429]	108	#endif
[1]	109
[1833]	110	INTERFACE spectra_check_parameters
	111	MODULE PROCEDURE spectra_check_parameters
	112	END INTERFACE spectra_check_parameters
[1]	113
[1833]	114	INTERFACE spectra_header
	115	MODULE PROCEDURE spectra_header
	116	END INTERFACE spectra_header
[1786]	117
[1833]	118	INTERFACE spectra_init
	119	MODULE PROCEDURE spectra_init
	120	END INTERFACE spectra_init
	121
	122	INTERFACE spectra_parin
	123	MODULE PROCEDURE spectra_parin
	124	END INTERFACE spectra_parin
	125
[4591]	126	PUBLIC average_count_sp, &
	127	averaging_interval_sp, &
	128	calc_spectra, &
	129	calculate_spectra, &
	130	comp_spectra_level, &
	131	data_output_sp, &
	132	dosp_time_count, &
	133	dt_dosp, &
	134	n_sp_x, &
	135	n_sp_y, &
	136	skip_time_dosp, &
	137	spectra_check_parameters, &
	138	spectra_direction, &
	139	spectra_header, &
	140	spectra_init, &
	141	spectra_parin, &
	142	spectrum_x, &
	143	spectrum_y, &
	144	var_d
[1833]	145
	146
[1786]	147	CONTAINS
	148
[4591]	149	!--------------------------------------------------------------------------------------------------!
[1833]	150	! Description:
	151	! ------------
	152	!> Parin for &spectra_par for calculating spectra
[4591]	153	!--------------------------------------------------------------------------------------------------!
	154	SUBROUTINE spectra_parin
[1833]	155
[4591]	156	USE control_parameters, &
[4842]	157	ONLY: dt_data_output
[1833]	158
[4591]	159	IMPLICIT NONE
[1833]	160
[4842]	161	CHARACTER(LEN=100) :: line !< dummy string that contains the current line of the parameter file
[1833]	162
[4842]	163	INTEGER(iwp) :: io_status !< status after reading the namelist file
[1833]	164
[4843]	165	LOGICAL :: switch_off_module = .FALSE. !< local namelist parameter to switch off the module
	166	!< although the respective module namelist appears in
	167	!< the namelist file
[4842]	168
[4591]	169	NAMELIST /spectra_parameters/ &
	170	averaging_interval_sp, &
	171	comp_spectra_level, &
	172	data_output_sp, &
	173	dt_dosp, &
	174	skip_time_dosp, &
[4843]	175	spectra_direction, &
	176	switch_off_module
[1833]	177	!
[4842]	178	!-- Position the namelist-file at the beginning (it was already opened in parin), and try to read
	179	!-- the namelist.
	180	REWIND( 11 )
	181	READ( 11, spectra_parameters, IOSTAT=io_status )
[1833]	182
	183	!
[4842]	184	!-- Action depending on the READ status
	185	IF ( io_status == 0 ) THEN
[1833]	186	!
[4843]	187	!-- spectra_parameters namelist was found and read correctly. Set flag to indicate that spectra
	188	!-- have to be calculated.
	189	IF ( .NOT. switch_off_module ) THEN
	190	calculate_spectra = .TRUE.
	191	!-- Default setting of dt_dosp here (instead of check_parameters), because its current value
	192	!-- is needed in init_pegrid.
	193	IF ( dt_dosp == 9999999.9_wp ) dt_dosp = dt_data_output
	194	ENDIF
[2932]	195
[4842]	196	ELSEIF ( io_status > 0 ) THEN
[2932]	197	!
[4842]	198	!-- spectra_parameters namelist was found but contained errors. Print an error message including
	199	!-- the line that caused the problem.
	200	BACKSPACE( 11 )
	201	READ( 11 , '(A)' ) line
	202	CALL parin_fail_message( 'spectra_parameters', line )
[2932]	203
[4842]	204	ENDIF
[3246]	205
[4591]	206	END SUBROUTINE spectra_parin
	207
	208
	209
	210	!--------------------------------------------------------------------------------------------------!
[1833]	211	! Description:
	212	! ------------
	213	!> Initialization of spectra related variables
[4591]	214	!--------------------------------------------------------------------------------------------------!
	215	SUBROUTINE spectra_init
[1833]	216
[4591]	217	USE indices, &
	218	ONLY: nx, &
	219	ny, &
	220	nzb, &
	221	nzt
[1833]	222
[4591]	223	IMPLICIT NONE
[1833]	224
[4591]	225	IF ( spectra_initialized ) RETURN
[1833]	226
[4591]	227	IF ( dt_dosp /= 9999999.9_wp ) THEN
[1833]	228
[4591]	229	IF ( .NOT. ALLOCATED( spectrum_x ) ) THEN
	230	ALLOCATE( spectrum_x( 1:nx/2, 1:100, 1:10 ) )
	231	spectrum_x = 0.0_wp
	232	ENDIF
[2956]	233
[4591]	234	IF ( .NOT. ALLOCATED( spectrum_y ) ) THEN
	235	ALLOCATE( spectrum_y( 1:ny/2, 1:100, 1:10 ) )
	236	spectrum_y = 0.0_wp
[1833]	237	ENDIF
	238
[4591]	239	ALLOCATE( var_d(nzb:nzt+1) )
	240	var_d = 0.0_wp
	241	ENDIF
[1833]	242
[4591]	243	spectra_initialized = .TRUE.
[1833]	244
[4591]	245	END SUBROUTINE spectra_init
[1833]	246
	247
[4591]	248
	249	!--------------------------------------------------------------------------------------------------!
[1833]	250	! Description:
	251	! ------------
	252	!> Check spectra related quantities
[4591]	253	!--------------------------------------------------------------------------------------------------!
	254	SUBROUTINE spectra_check_parameters
[1833]	255
[4591]	256	USE control_parameters, &
	257	ONLY: averaging_interval, &
	258	message_string, &
	259	skip_time_data_output
[1833]	260
[4591]	261	IMPLICIT NONE
[1833]	262
	263	!
[4591]	264	!-- Check the average interval
	265	IF ( averaging_interval_sp == 9999999.9_wp ) THEN
	266	averaging_interval_sp = averaging_interval
	267	ENDIF
[1833]	268
[4591]	269	IF ( averaging_interval_sp > dt_dosp ) THEN
	270	WRITE( message_string, * ) 'averaging_interval_sp = ', averaging_interval_sp, &
	271	' must be <= dt_dosp = ', dt_dosp
	272	CALL message( 'spectra_check_parameters', 'PA0087', 1, 2, 0, 6, 0 )
	273	ENDIF
[1833]	274
	275	!
[4591]	276	!-- Set the default skip time interval for data output, if necessary
	277	IF ( skip_time_dosp == 9999999.9_wp ) skip_time_dosp = skip_time_data_output
[1833]	278
[4591]	279	END SUBROUTINE spectra_check_parameters
[1833]	280
	281
	282
[4591]	283	!--------------------------------------------------------------------------------------------------!
[1833]	284	! Description:
	285	! ------------
	286	!> Header output for spectra
	287	!>
	288	!> @todo Output of netcdf data format and compression level
[4591]	289	!--------------------------------------------------------------------------------------------------!
	290	SUBROUTINE spectra_header ( io )
[1833]	291
[4591]	292	USE control_parameters, &
	293	ONLY: dt_averaging_input_pr
[1833]	294
[4591]	295	! USE netcdf_interface, &
	296	! ONLY: netcdf_data_format_string, &
	297	! netcdf_deflate
[1833]	298
[4591]	299	IMPLICIT NONE
[1833]	300
[4591]	301	! CHARACTER (LEN=40) :: output_format !< internal string
[1833]	302
[4591]	303	INTEGER(iwp) :: i !< internal counter
	304	INTEGER(iwp), INTENT(IN) :: io !< unit of the output file
[1833]	305
	306	!
[4591]	307	!-- Spectra output
	308	IF ( dt_dosp /= 9999999.9_wp ) THEN
	309	WRITE ( io, 1 )
[1833]	310
[4591]	311	! output_format = netcdf_data_format_string
	312	! IF ( netcdf_deflate == 0 ) THEN
	313	! WRITE ( io, 2 ) output_format
	314	! ELSE
	315	! WRITE ( io, 3 ) TRIM( output_format ), netcdf_deflate
	316	! ENDIF
	317	WRITE ( io, 2 ) 'see profiles or other quantities'
	318	WRITE ( io, 4 ) dt_dosp
	319	IF ( skip_time_dosp /= 0.0_wp ) WRITE ( io, 5 ) skip_time_dosp
	320	WRITE ( io, 6 ) ( data_output_sp(i), i = 1,10 ), &
	321	( spectra_direction(i), i = 1,10 ), &
	322	( comp_spectra_level(i), i = 1,100 ), &
	323	averaging_interval_sp, dt_averaging_input_pr
	324	ENDIF
[1833]	325
[4591]	326	1 FORMAT (' Spectra:')
	327	2 FORMAT (' Output format: ',A/)
	328	! 3 FORMAT (' Output format: ',A, ' compressed with level: ',I1/)
	329	4 FORMAT (' Output every ',F7.1,' s'/)
	330	5 FORMAT (' No output during initial ',F8.2,' s')
	331	6 FORMAT (' Arrays: ', 10(A5,',')/ &
	332	' Directions: ', 10(A5,',')/ &
	333	' height levels k = ', 20(I3,',')/ &
	334	' ', 20(I3,',')/ &
	335	' ', 20(I3,',')/ &
	336	' ', 20(I3,',')/ &
	337	' ', 19(I3,','),I3,'.'/ &
	338	' Time averaged over ', F7.1, ' s,' / &
	339	' Profiles for the time averaging are taken every ', &
	340	F6.1,' s')
[1833]	341
[4591]	342	END SUBROUTINE spectra_header
[1833]	343
	344
[4591]	345	!--------------------------------------------------------------------------------------------------!
	346	! Description:
	347	! ------------
	348	!> @todo Missing subroutine description.
	349	!--------------------------------------------------------------------------------------------------!
[1833]	350
[4591]	351	SUBROUTINE calc_spectra
[1786]	352
[4591]	353	USE arrays_3d, &
	354	ONLY: d
[4429]	355	#if defined( __parallel )
[4591]	356	USE arrays_3d, &
	357	ONLY: tend
[4429]	358	#endif
[1786]	359
[4591]	360	USE control_parameters, &
	361	ONLY: bc_lr_cyc, &
	362	bc_ns_cyc, &
	363	message_string, &
	364	psolver
[1786]	365
[4591]	366	USE cpulog, &
	367	ONLY: cpu_log, &
	368	log_point
[1786]	369
[4591]	370	USE fft_xy, &
	371	ONLY: fft_init
[1786]	372
[4591]	373	USE indices, &
	374	ONLY: nxl, &
	375	nxr, &
	376	nyn, &
	377	nys, &
	378	nzb, &
	379	nzt
[1786]	380
[4591]	381	USE kinds
[1786]	382
[4591]	383	USE pegrid, &
	384	ONLY: myid
[4429]	385	#if defined( __parallel )
[4591]	386	USE pegrid, &
	387	ONLY: pdims
[4429]	388	#endif
[1786]	389
[4591]	390	IMPLICIT NONE
[1786]	391
[4591]	392	INTEGER(iwp) :: m !<
	393	INTEGER(iwp) :: pr !<
[1786]	394
	395
[1]	396	!
[4591]	397	!-- Check if user gave any levels for spectra to be calculated
	398	IF ( comp_spectra_level(1) == 999999 ) RETURN
[1]	399
[4591]	400	CALL cpu_log( log_point(30), 'calc_spectra', 'start' )
[1786]	401
[1]	402	!
[4591]	403	!-- Initialize spectra related quantities
	404	CALL spectra_init
[1833]	405
	406	!
[4591]	407	!-- Initialize ffts
	408	CALL fft_init
[225]	409
[1]	410	!
[4591]	411	!-- Reallocate array d in required size
	412	IF ( psolver(1:9) == 'multigrid' ) THEN
	413	DEALLOCATE( d )
	414	ALLOCATE( d(nzb+1:nzt,nys:nyn,nxl:nxr) )
	415	ENDIF
[1]	416
[4591]	417	m = 1
	418	DO WHILE ( data_output_sp(m) /= ' ' .AND. m <= 10 )
[1]	419	!
[4591]	420	!-- Transposition from z --> x ( y --> x in case of a 1d-decomposition along x)
	421	IF ( INDEX( spectra_direction(m), 'x' ) /= 0 ) THEN
[247]	422
[1786]	423	!
[4591]	424	!-- Calculation of spectra works for cyclic boundary conditions only
	425	IF ( .NOT. bc_lr_cyc ) THEN
[1]	426
[4591]	427	message_string = 'non-cyclic lateral boundaries along x do' // &
	428	' not & allow calculation of spectra along x'
	429	CALL message( 'calc_spectra', 'PA0160', 1, 2, 0, 6, 0 )
	430	ENDIF
[1]	431
[4591]	432	CALL preprocess_spectra( m, pr )
[1786]	433
[1]	434	#if defined( __parallel )
[4591]	435	IF ( pdims(2) /= 1 ) THEN
	436	CALL resort_for_zx( d, tend )
	437	CALL transpose_zx( tend, d )
	438	ELSE
	439	CALL transpose_yxd( d, d )
	440	ENDIF
	441	CALL calc_spectra_x( d, m )
[1]	442	#else
[4591]	443	message_string = 'sorry, calculation of spectra in non parallel mode& is still not ' // &
	444	'realized'
	445	CALL message( 'calc_spectra', 'PA0161', 1, 2, 0, 6, 0 )
[1]	446	#endif
	447
[4591]	448	ENDIF
[1]	449
	450	!
[4591]	451	!-- Transposition from z --> y (d is rearranged only in case of a 1d-decomposition along x)
	452	IF ( INDEX( spectra_direction(m), 'y' ) /= 0 ) THEN
[1]	453
	454	!
[4591]	455	!-- Calculation of spectra works for cyclic boundary conditions only
	456	IF ( .NOT. bc_ns_cyc ) THEN
	457	IF ( myid == 0 ) THEN
	458	message_string = 'non-cyclic lateral boundaries along y' // &
	459	' do not & allow calculation of spectra along y'
	460	CALL message( 'calc_spectra', 'PA0162', 1, 2, 0, 6, 0 )
[1]	461	ENDIF
[4591]	462	CALL local_stop
	463	ENDIF
[1]	464
[4591]	465	CALL preprocess_spectra( m, pr )
[1]	466
	467	#if defined( __parallel )
[4591]	468	CALL transpose_zyd( d, d )
	469	CALL calc_spectra_y( d, m )
[1]	470	#else
[4591]	471	message_string = 'sorry, calculation of spectra in non parallel mode& is still not ' // &
	472	'realized'
	473	CALL message( 'calc_spectra', 'PA0161', 1, 2, 0, 6, 0 )
[1]	474	#endif
	475
[4591]	476	ENDIF
[1]	477
	478	!
[4591]	479	!-- Increase counter for next spectrum
	480	m = m + 1
[1]	481
[4591]	482	ENDDO
	483
[1]	484	!
[4591]	485	!-- Increase counter for averaging process in routine plot_spectra
	486	average_count_sp = average_count_sp + 1
[1]	487
[4591]	488	CALL cpu_log( log_point(30), 'calc_spectra', 'stop' )
[1]	489
[4591]	490	END SUBROUTINE calc_spectra
[1]	491
	492
[4591]	493	!--------------------------------------------------------------------------------------------------!
[1682]	494	! Description:
	495	! ------------
	496	!> @todo Missing subroutine description.
[4591]	497	!--------------------------------------------------------------------------------------------------!
	498	SUBROUTINE preprocess_spectra( m, pr )
[1]	499
[4591]	500	USE arrays_3d, &
	501	ONLY: d, &
	502	pt, &
	503	q, &
	504	s, &
	505	u, &
	506	v, &
	507	w
[1320]	508
[4591]	509	USE indices, &
	510	ONLY: ngp_2dh, &
	511	nxl, &
	512	nxr, &
	513	nyn, &
	514	nys, &
	515	nzb, &
	516	nzt
[1320]	517
[4591]	518	USE kinds
[1320]	519
[1815]	520	#if defined( __parallel )
[4591]	521	USE MPI
[2841]	522
[4591]	523	USE pegrid, &
	524	ONLY: collective_wait, &
	525	comm2d, &
	526	ierr
[4429]	527	#endif
[1]	528
[4591]	529	USE statistics, &
	530	ONLY: hom
[1320]	531
	532
[4591]	533	IMPLICIT NONE
[1320]	534
[4591]	535	INTEGER(iwp) :: i !<
	536	INTEGER(iwp) :: j !<
	537	INTEGER(iwp) :: k !<
	538	INTEGER(iwp) :: m !<
	539	INTEGER(iwp) :: pr !<
[1]	540
[4591]	541	REAL(wp), DIMENSION(nzb:nzt+1) :: var_d_l !<
[1]	542
[4629]	543
[4591]	544	SELECT CASE ( TRIM( data_output_sp(m) ) )
	545
	546	CASE ( 'u' )
	547	pr = 1
	548	d(nzb+1:nzt,nys:nyn,nxl:nxr) = u(nzb+1:nzt,nys:nyn,nxl:nxr)
	549
	550	CASE ( 'v' )
	551	pr = 2
	552	d(nzb+1:nzt,nys:nyn,nxl:nxr) = v(nzb+1:nzt,nys:nyn,nxl:nxr)
	553
	554	CASE ( 'w' )
	555	pr = 3
	556	d(nzb+1:nzt,nys:nyn,nxl:nxr) = w(nzb+1:nzt,nys:nyn,nxl:nxr)
	557
	558	CASE ( 'theta' )
	559	pr = 4
	560	d(nzb+1:nzt,nys:nyn,nxl:nxr) = pt(nzb+1:nzt,nys:nyn,nxl:nxr)
	561
	562	CASE ( 'q' )
	563	pr = 41
	564	d(nzb+1:nzt,nys:nyn,nxl:nxr) = q(nzb+1:nzt,nys:nyn,nxl:nxr)
	565
	566	CASE ( 's' )
	567	pr = 117
	568	d(nzb+1:nzt,nys:nyn,nxl:nxr) = s(nzb+1:nzt,nys:nyn,nxl:nxr)
	569
	570	CASE DEFAULT
[144]	571	!
[4591]	572	!-- The DEFAULT case is reached either if the parameter data_output_sp(m) contains a wrong
	573	!-- character string or if the user has coded a special case in the user interface. There, the
	574	!-- subroutine user_spectra checks which of these two conditions applies.
	575	CALL user_spectra( 'preprocess', m, pr )
[1]	576
[4591]	577	END SELECT
	578
[1]	579	!
[4591]	580	!-- Subtract horizontal mean from the array, for which spectra have to be calculated. Moreover,
	581	!-- calculate variance of the respective quantitiy, later used for normalizing spectra output.
	582	var_d_l(:) = 0.0_wp
	583	DO i = nxl, nxr
	584	DO j = nys, nyn
	585	DO k = nzb+1, nzt
	586	d(k,j,i) = d(k,j,i) - hom(k,1,pr,0)
	587	var_d_l(k) = var_d_l(k) + d(k,j,i) * d(k,j,i)
[1]	588	ENDDO
	589	ENDDO
[4591]	590	ENDDO
[1431]	591	!
[4591]	592	!-- Compute total variance from local variances
	593	var_d(:) = 0.0_wp
	594	#if defined( __parallel )
	595	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
	596	CALL MPI_ALLREDUCE( var_d_l(0), var_d(0), nzt+1-nzb, MPI_REAL, MPI_SUM, comm2d, ierr )
[1431]	597	#else
[4591]	598	var_d(:) = var_d_l(:)
[1431]	599	#endif
[4591]	600	var_d(:) = var_d(:) / ngp_2dh(0)
[1]	601
[4591]	602	END SUBROUTINE preprocess_spectra
[1]	603
	604
[4591]	605	!--------------------------------------------------------------------------------------------------!
[1682]	606	! Description:
	607	! ------------
	608	!> @todo Missing subroutine description.
[4591]	609	!--------------------------------------------------------------------------------------------------!
[4429]	610	#if defined( __parallel )
[4591]	611	SUBROUTINE calc_spectra_x( ddd, m )
[1]	612
[4591]	613	USE fft_xy, &
	614	ONLY: fft_x_1d
[1320]	615
[4591]	616	USE grid_variables, &
	617	ONLY: dx
[1320]	618
[4591]	619	USE indices, &
	620	ONLY: nx, &
	621	ny
[1320]	622
[4591]	623	USE kinds
[1320]	624
[4591]	625	USE MPI
[2841]	626
[4591]	627	USE pegrid, &
	628	ONLY: comm2d, &
	629	ierr, &
	630	myid
[1320]	631
[4591]	632	USE transpose_indices, &
	633	ONLY: nyn_x, &
	634	nys_x, &
	635	nzb_x, &
	636	nzt_x
[1320]	637
	638
[4591]	639	IMPLICIT NONE
[1]	640
[4591]	641	INTEGER(iwp) :: i !<
	642	INTEGER(iwp) :: j !<
	643	INTEGER(iwp) :: k !<
	644	INTEGER(iwp) :: m !<
	645	INTEGER(iwp) :: n !<
[1320]	646
[4591]	647	REAL(wp) :: exponent !<
	648	REAL(wp) :: sum_spec_dum !< wavenumber-integrated spectrum
[1]	649
[4591]	650	REAL(wp), DIMENSION(0:nx) :: work !<
	651
	652	REAL(wp), DIMENSION(0:nx/2) :: sums_spectra_l !<
	653
	654	REAL(wp), DIMENSION(0:nx/2,100) :: sums_spectra !<
	655
	656	REAL(wp), DIMENSION(0:nx,nys_x:nyn_x,nzb_x:nzt_x) :: ddd !<
	657
[1]	658	!
[4591]	659	!-- Exponent for geometric average
	660	exponent = 1.0_wp / ( ny + 1.0_wp )
[1]	661
	662	!
[4591]	663	!-- Loop over all levels defined by the user
	664	n = 1
	665	DO WHILE ( comp_spectra_level(n) /= 999999 .AND. n <= 100 )
[1]	666
[4591]	667	k = comp_spectra_level(n)
[1]	668
	669	!
[4591]	670	!-- Calculate FFT only if the corresponding level is situated on this PE
	671	IF ( k >= nzb_x .AND. k <= nzt_x ) THEN
[1]	672
[4591]	673	DO j = nys_x, nyn_x
[1]	674
[4591]	675	work = ddd(0:nx,j,k)
	676	CALL fft_x_1d( work, 'forward' )
[1786]	677
[4591]	678	ddd(0,j,k) = dx * work(0)**2
	679	DO i = 1, nx/2
	680	ddd(i,j,k) = dx * ( work(i)2 + work(nx+1-i)2 )
[1]	681	ENDDO
	682
[4591]	683	ENDDO
	684
[1]	685	!
[4591]	686	!-- Local sum and geometric average of these spectra (WARNING: no global sum should be
	687	!-- performed, because floating point overflow may occur)
	688	DO i = 0, nx/2
[1]	689
[4591]	690	sums_spectra_l(i) = 1.0_wp
	691	DO j = nys_x, nyn_x
	692	sums_spectra_l(i) = sums_spectra_l(i) * ddd(i,j,k)**exponent
[1]	693	ENDDO
	694
[4591]	695	ENDDO
[1]	696
[4591]	697	ELSE
[1]	698
[4591]	699	sums_spectra_l = 1.0_wp
	700
	701	ENDIF
	702
[1]	703	!
[4591]	704	!-- Global sum of spectra on PE0 (from where they are written on file)
	705	sums_spectra(:,n) = 0.0_wp
	706	#if defined( __parallel )
	707	CALL MPI_BARRIER( comm2d, ierr ) ! Necessary?
	708	CALL MPI_REDUCE( sums_spectra_l(0), sums_spectra(0,n), nx/2+1, MPI_REAL, MPI_PROD, 0, &
	709	comm2d, ierr )
[1]	710	#else
[4591]	711	sums_spectra(:,n) = sums_spectra_l
[1]	712	#endif
[1431]	713	!
[4591]	714	!-- Normalize spectra by variance
	715	sum_spec_dum = SUM( sums_spectra(1:nx/2,n) )
[2215]	716
[4591]	717	IF ( sum_spec_dum /= 0.0_wp ) THEN
	718	sums_spectra(1:nx/2,n) = sums_spectra(1:nx/2,n) * var_d(k) / sum_spec_dum
	719	ENDIF
	720	n = n + 1
[1]	721
[4591]	722	ENDDO
	723	n = n - 1
[1]	724
[4591]	725	IF ( myid == 0 ) THEN
[1]	726	!
[4591]	727	!-- Sum of spectra for later averaging (see routine data_output_spectra)
	728	DO i = 1, nx/2
	729	DO k = 1, n
	730	spectrum_x(i,k,m) = spectrum_x(i,k,m) + sums_spectra(i,k)
[1]	731	ENDDO
[4591]	732	ENDDO
[1]	733
[4591]	734	ENDIF
[1]	735	!
[4591]	736	!-- n_sp_x is needed by data_output_spectra_x
	737	n_sp_x = n
[1]	738
[4591]	739	END SUBROUTINE calc_spectra_x
[4429]	740	#endif
[1]	741
	742
[4591]	743	!--------------------------------------------------------------------------------------------------!
[1682]	744	! Description:
	745	! ------------
	746	!> @todo Missing subroutine description.
[4591]	747	!--------------------------------------------------------------------------------------------------!
[4429]	748	#if defined( __parallel )
[4591]	749	SUBROUTINE calc_spectra_y( ddd, m )
[1]	750
[4591]	751	USE fft_xy, &
	752	ONLY: fft_y_1d
[1320]	753
[4591]	754	USE grid_variables, &
	755	ONLY: dy
[1320]	756
[4591]	757	USE indices, &
	758	ONLY: nx, &
	759	ny
[1320]	760
[4591]	761	USE kinds
[1320]	762
[4591]	763	USE MPI
[2841]	764
[4591]	765	USE pegrid, &
	766	ONLY: comm2d, &
	767	ierr, &
	768	myid
[1320]	769
[4591]	770	USE transpose_indices, &
	771	ONLY: nxl_yd, &
	772	nxr_yd, &
	773	nzb_yd, &
	774	nzt_yd
[1320]	775
	776
[4591]	777	IMPLICIT NONE
[1]	778
[4591]	779	INTEGER(iwp) :: i !<
	780	INTEGER(iwp) :: j !<
	781	INTEGER(iwp) :: k !<
	782	INTEGER(iwp) :: m !<
	783	INTEGER(iwp) :: n !<
[1320]	784
[4591]	785	REAL(wp) :: exponent !<
	786	REAL(wp) :: sum_spec_dum !< wavenumber-integrated spectrum
[1]	787
[4591]	788	REAL(wp), DIMENSION(0:ny) :: work !<
[1]	789
[4591]	790	REAL(wp), DIMENSION(0:ny/2) :: sums_spectra_l !<
	791
	792	REAL(wp), DIMENSION(0:ny/2,100) :: sums_spectra !<
	793
	794	REAL(wp), DIMENSION(0:ny,nxl_yd:nxr_yd,nzb_yd:nzt_yd) :: ddd !<
	795
	796
[1]	797	!
[4591]	798	!-- Exponent for geometric average
	799	exponent = 1.0_wp / ( nx + 1.0_wp )
[1]	800
	801	!
[4591]	802	!-- Loop over all levels defined by the user
	803	n = 1
	804	DO WHILE ( comp_spectra_level(n) /= 999999 .AND. n <= 100 )
[1]	805
[4591]	806	k = comp_spectra_level(n)
[1]	807
	808	!
[4591]	809	!-- Calculate FFT only if the corresponding level is situated on this PE
	810	IF ( k >= nzb_yd .AND. k <= nzt_yd ) THEN
[1]	811
[4591]	812	DO i = nxl_yd, nxr_yd
[1]	813
[4591]	814	work = ddd(0:ny,i,k)
	815	CALL fft_y_1d( work, 'forward' )
[1786]	816
[4591]	817	ddd(0,i,k) = dy * work(0)**2
	818	DO j = 1, ny/2
	819	ddd(j,i,k) = dy * ( work(j)2 + work(ny+1-j)2 )
[1]	820	ENDDO
	821
[4591]	822	ENDDO
	823
[1]	824	!
[4591]	825	!-- Local sum and geometric average of these spectra (WARNING: no global sum should be
	826	!-- performed, because floating point overflow may occur)
	827	DO j = 0, ny/2
[1]	828
[4591]	829	sums_spectra_l(j) = 1.0_wp
	830	DO i = nxl_yd, nxr_yd
	831	sums_spectra_l(j) = sums_spectra_l(j) * ddd(j,i,k)**exponent
[1]	832	ENDDO
	833
[4591]	834	ENDDO
[1]	835
[4591]	836	ELSE
[1]	837
[4591]	838	sums_spectra_l = 1.0_wp
	839
	840	ENDIF
	841
[1]	842	!
[4591]	843	!-- Global sum of spectra on PE0 (from where they are written on file)
	844	sums_spectra(:,n) = 0.0_wp
	845	#if defined( __parallel )
	846	CALL MPI_BARRIER( comm2d, ierr ) ! Necessary?
	847	CALL MPI_REDUCE( sums_spectra_l(0), sums_spectra(0,n), ny/2+1, MPI_REAL, MPI_PROD, 0, &
	848	comm2d, ierr )
[1]	849	#else
[4591]	850	sums_spectra(:,n) = sums_spectra_l
[1]	851	#endif
[1431]	852	!
[4591]	853	!-- Normalize spectra by variance
	854	sum_spec_dum = SUM( sums_spectra(1:ny/2,n) )
	855	IF ( sum_spec_dum /= 0.0_wp ) THEN
	856	sums_spectra(1:ny/2,n) = sums_spectra(1:ny/2,n) * var_d(k) / sum_spec_dum
	857	ENDIF
	858	n = n + 1
[1]	859
[4591]	860	ENDDO
	861	n = n - 1
[1]	862
	863
[4591]	864	IF ( myid == 0 ) THEN
[1]	865	!
[4591]	866	!-- Sum of spectra for later averaging (see routine data_output_spectra)
	867	DO j = 1, ny/2
	868	DO k = 1, n
	869	spectrum_y(j,k,m) = spectrum_y(j,k,m) + sums_spectra(j,k)
[1]	870	ENDDO
[4591]	871	ENDDO
[1]	872
[4591]	873	ENDIF
[1]	874
	875	!
[4591]	876	!-- n_sp_y is needed by data_output_spectra_y
	877	n_sp_y = n
[1]	878
[4591]	879	END SUBROUTINE calc_spectra_y
[4429]	880	#endif
[1786]	881
[1833]	882	END MODULE spectra_mod

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

Context Navigation

source: palm/trunk/SOURCE/spectra_mod.f90 @ 4847

Download in other formats: