Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |