Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |