Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |