Home

Context Navigation

source: palm/trunk/SOURCE/tridia_solver.f90 @ 1806

Last change on this file since 1806 was 1805, checked in by maronga, 9 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 23.0 KB

Rev	Line
[1682]	1	!> @file tridia_solver.f90
[1212]	2	!--------------------------------------------------------------------------------!
	3	! This file is part of PALM.
	4	!
	5	! PALM is free software: you can redistribute it and/or modify it under the terms
	6	! of the GNU General Public License as published by the Free Software Foundation,
	7	! either version 3 of the License, or (at your option) any later version.
	8	!
	9	! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
	10	! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
	11	! A PARTICULAR PURPOSE. See the GNU General Public License for more details.
	12	!
	13	! You should have received a copy of the GNU General Public License along with
	14	! PALM. If not, see <http://www.gnu.org/licenses/>.
	15	!
[1310]	16	! Copyright 1997-2014 Leibniz Universitaet Hannover
[1212]	17	!--------------------------------------------------------------------------------!
	18	!
	19	! Current revisions:
	20	! ------------------
[1343]	21	!
[1805]	22	!
[1321]	23	! Former revisions:
	24	! -----------------
	25	! $Id: tridia_solver.f90 1805 2016-04-05 16:32:34Z gronemeier $
	26	!
[1805]	27	! 1804 2016-04-05 16:30:18Z maronga
	28	! Removed code for parameter file check (__check)
	29	!
[1683]	30	! 1682 2015-10-07 23:56:08Z knoop
	31	! Code annotations made doxygen readable
	32	!
[1407]	33	! 1406 2014-05-16 13:47:01Z raasch
	34	! bugfix for pgi 14.4: declare create moved after array declaration
	35	!
[1343]	36	! 1342 2014-03-26 17:04:47Z kanani
	37	! REAL constants defined as wp-kind
	38	!
[1323]	39	! 1322 2014-03-20 16:38:49Z raasch
	40	! REAL functions provided with KIND-attribute
	41	!
[1321]	42	! 1320 2014-03-20 08:40:49Z raasch
[1320]	43	! ONLY-attribute added to USE-statements,
	44	! kind-parameters added to all INTEGER and REAL declaration statements,
	45	! kinds are defined in new module kinds,
	46	! old module precision_kind is removed,
	47	! revision history before 2012 removed,
	48	! comment fields (!:) to be used for variable explanations added to
	49	! all variable declaration statements
[1213]	50	!
[1258]	51	! 1257 2013-11-08 15:18:40Z raasch
	52	! openacc loop and loop vector clauses removed, declare create moved after
	53	! the FORTRAN declaration statement
	54	!
[1222]	55	! 1221 2013-09-10 08:59:13Z raasch
	56	! dummy argument tri in 1d-routines replaced by tri_for_1d because of name
	57	! conflict with arry tri in module arrays_3d
	58	!
[1217]	59	! 1216 2013-08-26 09:31:42Z raasch
	60	! +tridia_substi_overlap for handling overlapping fft / transposition
	61	!
[1213]	62	! 1212 2013-08-15 08:46:27Z raasch
[1212]	63	! Initial revision.
	64	! Routines have been moved to seperate module from former file poisfft to here.
	65	! The tridiagonal matrix coefficients of array tri are calculated only once at
	66	! the beginning, i.e. routine split is called within tridia_init.
	67	!
	68	!
	69	! Description:
	70	! ------------
[1682]	71	!> solves the linear system of equations:
	72	!>
	73	!> -(4 pi^2(i^2/(dx^2nnx^2)+j^2/(dy^2nny^2))+
	74	!> 1/(dzu(k)dzw(k))+1/(dzu(k-1)dzw(k)))*p(i,j,k)+
	75	!> 1/(dzu(k)dzw(k))p(i,j,k+1)+1/(dzu(k-1)dzw(k))p(i,j,k-1)=d(i,j,k)
	76	!>
	77	!> by using the Thomas algorithm
[1212]	78	!------------------------------------------------------------------------------!
[1682]	79	MODULE tridia_solver
	80
[1212]	81
[1320]	82	USE indices, &
	83	ONLY: nx, ny, nz
[1212]	84
[1320]	85	USE kinds
	86
	87	USE transpose_indices, &
	88	ONLY: nxl_z, nyn_z, nxr_z, nys_z
	89
[1212]	90	IMPLICIT NONE
	91
[1682]	92	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: ddzuw !<
[1212]	93
	94	PRIVATE
	95
	96	INTERFACE tridia_substi
	97	MODULE PROCEDURE tridia_substi
	98	END INTERFACE tridia_substi
	99
[1216]	100	INTERFACE tridia_substi_overlap
	101	MODULE PROCEDURE tridia_substi_overlap
	102	END INTERFACE tridia_substi_overlap
[1212]	103
[1216]	104	PUBLIC tridia_substi, tridia_substi_overlap, tridia_init, tridia_1dd
	105
[1212]	106	CONTAINS
	107
	108
[1682]	109	!------------------------------------------------------------------------------!
	110	! Description:
	111	! ------------
	112	!> @todo Missing subroutine description.
	113	!------------------------------------------------------------------------------!
[1212]	114	SUBROUTINE tridia_init
	115
[1320]	116	USE arrays_3d, &
	117	ONLY: ddzu_pres, ddzw
[1212]	118
[1320]	119	USE kinds
	120
[1212]	121	IMPLICIT NONE
	122
[1682]	123	INTEGER(iwp) :: k !<
[1212]	124
	125	ALLOCATE( ddzuw(0:nz-1,3) )
	126
	127	DO k = 0, nz-1
	128	ddzuw(k,1) = ddzu_pres(k+1) * ddzw(k+1)
	129	ddzuw(k,2) = ddzu_pres(k+2) * ddzw(k+1)
[1342]	130	ddzuw(k,3) = -1.0_wp * &
[1212]	131	( ddzu_pres(k+2) * ddzw(k+1) + ddzu_pres(k+1) * ddzw(k+1) )
	132	ENDDO
	133	!
	134	!-- Calculate constant coefficients of the tridiagonal matrix
	135	CALL maketri
	136	CALL split
	137
	138	END SUBROUTINE tridia_init
	139
	140
	141	!------------------------------------------------------------------------------!
[1682]	142	! Description:
	143	! ------------
	144	!> Computes the i- and j-dependent component of the matrix
	145	!> Provide the constant coefficients of the tridiagonal matrix for solution
	146	!> of the Poisson equation in Fourier space.
	147	!> The coefficients are computed following the method of
	148	!> Schmidt et al. (DFVLR-Mitteilung 84-15), which departs from Stephan
	149	!> Siano's original version by discretizing the Poisson equation,
	150	!> before it is Fourier-transformed.
[1212]	151	!------------------------------------------------------------------------------!
[1682]	152	SUBROUTINE maketri
[1212]	153
[1682]	154
[1320]	155	USE arrays_3d, &
	156	ONLY: tric
[1212]	157
[1320]	158	USE constants, &
	159	ONLY: pi
	160
	161	USE control_parameters, &
	162	ONLY: ibc_p_b, ibc_p_t
	163
	164	USE grid_variables, &
	165	ONLY: dx, dy
	166
	167
	168	USE kinds
	169
[1212]	170	IMPLICIT NONE
	171
[1682]	172	INTEGER(iwp) :: i !<
	173	INTEGER(iwp) :: j !<
	174	INTEGER(iwp) :: k !<
	175	INTEGER(iwp) :: nnxh !<
	176	INTEGER(iwp) :: nnyh !<
[1212]	177
[1682]	178	REAL(wp) :: ll(nxl_z:nxr_z,nys_z:nyn_z) !<
[1212]	179	!$acc declare create( ll )
	180
	181
	182	nnxh = ( nx + 1 ) / 2
	183	nnyh = ( ny + 1 ) / 2
	184
	185	!$acc kernels present( tric )
	186	DO j = nys_z, nyn_z
	187	DO i = nxl_z, nxr_z
	188	IF ( j >= 0 .AND. j <= nnyh ) THEN
	189	IF ( i >= 0 .AND. i <= nnxh ) THEN
[1342]	190	ll(i,j) = 2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * i ) / &
	191	REAL( nx+1, KIND=wp ) ) ) / ( dx * dx ) + &
	192	2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * j ) / &
	193	REAL( ny+1, KIND=wp ) ) ) / ( dy * dy )
[1212]	194	ELSE
[1342]	195	ll(i,j) = 2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * ( nx+1-i ) ) / &
	196	REAL( nx+1, KIND=wp ) ) ) / ( dx * dx ) + &
	197	2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * j ) / &
	198	REAL( ny+1, KIND=wp ) ) ) / ( dy * dy )
[1212]	199	ENDIF
	200	ELSE
	201	IF ( i >= 0 .AND. i <= nnxh ) THEN
[1342]	202	ll(i,j) = 2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * i ) / &
	203	REAL( nx+1, KIND=wp ) ) ) / ( dx * dx ) + &
	204	2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * ( ny+1-j ) ) / &
	205	REAL( ny+1, KIND=wp ) ) ) / ( dy * dy )
[1212]	206	ELSE
[1342]	207	ll(i,j) = 2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * ( nx+1-i ) ) / &
	208	REAL( nx+1, KIND=wp ) ) ) / ( dx * dx ) + &
	209	2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * ( ny+1-j ) ) / &
	210	REAL( ny+1, KIND=wp ) ) ) / ( dy * dy )
[1212]	211	ENDIF
	212	ENDIF
	213	ENDDO
	214	ENDDO
	215
	216	DO k = 0, nz-1
	217	DO j = nys_z, nyn_z
	218	DO i = nxl_z, nxr_z
	219	tric(i,j,k) = ddzuw(k,3) - ll(i,j)
	220	ENDDO
	221	ENDDO
	222	ENDDO
	223	!$acc end kernels
	224
	225	IF ( ibc_p_b == 1 ) THEN
	226	!$acc kernels present( tric )
	227	DO j = nys_z, nyn_z
	228	DO i = nxl_z, nxr_z
	229	tric(i,j,0) = tric(i,j,0) + ddzuw(0,1)
	230	ENDDO
	231	ENDDO
	232	!$acc end kernels
	233	ENDIF
	234	IF ( ibc_p_t == 1 ) THEN
	235	!$acc kernels present( tric )
	236	DO j = nys_z, nyn_z
	237	DO i = nxl_z, nxr_z
	238	tric(i,j,nz-1) = tric(i,j,nz-1) + ddzuw(nz-1,2)
	239	ENDDO
	240	ENDDO
	241	!$acc end kernels
	242	ENDIF
	243
	244	END SUBROUTINE maketri
	245
	246
	247	!------------------------------------------------------------------------------!
[1682]	248	! Description:
	249	! ------------
	250	!> Substitution (Forward and Backward) (Thomas algorithm)
[1212]	251	!------------------------------------------------------------------------------!
[1682]	252	SUBROUTINE tridia_substi( ar )
[1212]	253
[1682]	254
[1320]	255	USE arrays_3d, &
	256	ONLY: tri
[1212]	257
[1320]	258	USE control_parameters, &
	259	ONLY: ibc_p_b, ibc_p_t
	260
	261	USE kinds
	262
[1212]	263	IMPLICIT NONE
	264
[1682]	265	INTEGER(iwp) :: i !<
	266	INTEGER(iwp) :: j !<
	267	INTEGER(iwp) :: k !<
[1212]	268
[1682]	269	REAL(wp) :: ar(nxl_z:nxr_z,nys_z:nyn_z,1:nz) !<
[1212]	270
[1682]	271	REAL(wp), DIMENSION(nxl_z:nxr_z,nys_z:nyn_z,0:nz-1) :: ar1 !<
[1212]	272	!$acc declare create( ar1 )
	273
	274	!
	275	!-- Forward substitution
	276	DO k = 0, nz - 1
	277	!$acc kernels present( ar, tri )
	278	DO j = nys_z, nyn_z
	279	DO i = nxl_z, nxr_z
	280
	281	IF ( k == 0 ) THEN
	282	ar1(i,j,k) = ar(i,j,k+1)
	283	ELSE
	284	ar1(i,j,k) = ar(i,j,k+1) - tri(i,j,k,2) * ar1(i,j,k-1)
	285	ENDIF
	286
	287	ENDDO
	288	ENDDO
	289	!$acc end kernels
	290	ENDDO
	291
	292	!
	293	!-- Backward substitution
	294	!-- Note, the 1.0E-20 in the denominator is due to avoid divisions
	295	!-- by zero appearing if the pressure bc is set to neumann at the top of
	296	!-- the model domain.
	297	DO k = nz-1, 0, -1
	298	!$acc kernels present( ar, tri )
	299	DO j = nys_z, nyn_z
	300	DO i = nxl_z, nxr_z
	301
	302	IF ( k == nz-1 ) THEN
[1342]	303	ar(i,j,k+1) = ar1(i,j,k) / ( tri(i,j,k,1) + 1.0E-20_wp )
[1212]	304	ELSE
	305	ar(i,j,k+1) = ( ar1(i,j,k) - ddzuw(k,2) * ar(i,j,k+2) ) &
	306	/ tri(i,j,k,1)
	307	ENDIF
	308	ENDDO
	309	ENDDO
	310	!$acc end kernels
	311	ENDDO
	312
	313	!
	314	!-- Indices i=0, j=0 correspond to horizontally averaged pressure.
	315	!-- The respective values of ar should be zero at all k-levels if
	316	!-- acceleration of horizontally averaged vertical velocity is zero.
	317	IF ( ibc_p_b == 1 .AND. ibc_p_t == 1 ) THEN
	318	IF ( nys_z == 0 .AND. nxl_z == 0 ) THEN
	319	!$acc kernels loop present( ar )
	320	DO k = 1, nz
[1342]	321	ar(nxl_z,nys_z,k) = 0.0_wp
[1212]	322	ENDDO
[1257]	323	!$acc end kernels loop
[1212]	324	ENDIF
	325	ENDIF
	326
	327	END SUBROUTINE tridia_substi
	328
	329
[1216]	330	!------------------------------------------------------------------------------!
[1682]	331	! Description:
	332	! ------------
	333	!> Substitution (Forward and Backward) (Thomas algorithm)
[1216]	334	!------------------------------------------------------------------------------!
[1682]	335	SUBROUTINE tridia_substi_overlap( ar, jj )
[1216]	336
[1682]	337
[1320]	338	USE arrays_3d, &
	339	ONLY: tri
[1216]	340
[1320]	341	USE control_parameters, &
	342	ONLY: ibc_p_b, ibc_p_t
	343
	344	USE kinds
	345
[1216]	346	IMPLICIT NONE
	347
[1682]	348	INTEGER(iwp) :: i !<
	349	INTEGER(iwp) :: j !<
	350	INTEGER(iwp) :: jj !<
	351	INTEGER(iwp) :: k !<
[1216]	352
[1682]	353	REAL(wp) :: ar(nxl_z:nxr_z,nys_z:nyn_z,1:nz) !<
[1216]	354
[1682]	355	REAL(wp), DIMENSION(nxl_z:nxr_z,nys_z:nyn_z,0:nz-1) :: ar1 !<
[1216]	356	!$acc declare create( ar1 )
	357
	358	!
	359	!-- Forward substitution
	360	DO k = 0, nz - 1
	361	!$acc kernels present( ar, tri )
	362	!$acc loop
	363	DO j = nys_z, nyn_z
	364	DO i = nxl_z, nxr_z
	365
	366	IF ( k == 0 ) THEN
	367	ar1(i,j,k) = ar(i,j,k+1)
	368	ELSE
	369	ar1(i,j,k) = ar(i,j,k+1) - tri(i,jj,k,2) * ar1(i,j,k-1)
	370	ENDIF
	371
	372	ENDDO
	373	ENDDO
	374	!$acc end kernels
	375	ENDDO
	376
	377	!
	378	!-- Backward substitution
	379	!-- Note, the 1.0E-20 in the denominator is due to avoid divisions
	380	!-- by zero appearing if the pressure bc is set to neumann at the top of
	381	!-- the model domain.
	382	DO k = nz-1, 0, -1
	383	!$acc kernels present( ar, tri )
	384	!$acc loop
	385	DO j = nys_z, nyn_z
	386	DO i = nxl_z, nxr_z
	387
	388	IF ( k == nz-1 ) THEN
[1342]	389	ar(i,j,k+1) = ar1(i,j,k) / ( tri(i,jj,k,1) + 1.0E-20_wp )
[1216]	390	ELSE
	391	ar(i,j,k+1) = ( ar1(i,j,k) - ddzuw(k,2) * ar(i,j,k+2) ) &
	392	/ tri(i,jj,k,1)
	393	ENDIF
	394	ENDDO
	395	ENDDO
	396	!$acc end kernels
	397	ENDDO
	398
	399	!
	400	!-- Indices i=0, j=0 correspond to horizontally averaged pressure.
	401	!-- The respective values of ar should be zero at all k-levels if
	402	!-- acceleration of horizontally averaged vertical velocity is zero.
	403	IF ( ibc_p_b == 1 .AND. ibc_p_t == 1 ) THEN
	404	IF ( nys_z == 0 .AND. nxl_z == 0 ) THEN
	405	!$acc kernels loop present( ar )
	406	DO k = 1, nz
[1342]	407	ar(nxl_z,nys_z,k) = 0.0_wp
[1216]	408	ENDDO
	409	ENDIF
	410	ENDIF
	411
	412	END SUBROUTINE tridia_substi_overlap
	413
	414
[1212]	415	!------------------------------------------------------------------------------!
[1682]	416	! Description:
	417	! ------------
	418	!> Splitting of the tridiagonal matrix (Thomas algorithm)
[1212]	419	!------------------------------------------------------------------------------!
[1682]	420	SUBROUTINE split
[1212]	421
[1682]	422
[1320]	423	USE arrays_3d, &
	424	ONLY: tri, tric
[1212]	425
[1320]	426	USE kinds
	427
[1212]	428	IMPLICIT NONE
	429
[1682]	430	INTEGER(iwp) :: i !<
	431	INTEGER(iwp) :: j !<
	432	INTEGER(iwp) :: k !<
[1212]	433	!
	434	!-- Splitting
	435	!$acc kernels present( tri, tric )
	436	!$acc loop
	437	DO j = nys_z, nyn_z
	438	!$acc loop vector( 32 )
	439	DO i = nxl_z, nxr_z
	440	tri(i,j,0,1) = tric(i,j,0)
	441	ENDDO
	442	ENDDO
	443	!$acc end kernels
	444
	445	DO k = 1, nz-1
	446	!$acc kernels present( tri, tric )
	447	!$acc loop
	448	DO j = nys_z, nyn_z
	449	!$acc loop vector( 32 )
	450	DO i = nxl_z, nxr_z
	451	tri(i,j,k,2) = ddzuw(k,1) / tri(i,j,k-1,1)
	452	tri(i,j,k,1) = tric(i,j,k) - ddzuw(k-1,2) * tri(i,j,k,2)
	453	ENDDO
	454	ENDDO
	455	!$acc end kernels
	456	ENDDO
	457
	458	END SUBROUTINE split
	459
	460
	461	!------------------------------------------------------------------------------!
[1682]	462	! Description:
	463	! ------------
	464	!> Solves the linear system of equations for a 1d-decomposition along x (see
	465	!> tridia)
	466	!>
	467	!> @attention when using the intel compilers older than 12.0, array tri must
	468	!> be passed as an argument to the contained subroutines. Otherwise
	469	!> addres faults will occur. This feature can be activated with
	470	!> cpp-switch __intel11
	471	!> On NEC, tri should not be passed (except for routine substi_1dd)
	472	!> because this causes very bad performance.
[1212]	473	!------------------------------------------------------------------------------!
[1682]	474
	475	SUBROUTINE tridia_1dd( ddx2, ddy2, nx, ny, j, ar, tri_for_1d )
[1212]	476
[1682]	477
[1320]	478	USE arrays_3d, &
	479	ONLY: ddzu_pres, ddzw
[1212]	480
[1320]	481	USE control_parameters, &
	482	ONLY: ibc_p_b, ibc_p_t
[1212]	483
[1320]	484	USE kinds
	485
[1212]	486	IMPLICIT NONE
	487
[1682]	488	INTEGER(iwp) :: i !<
	489	INTEGER(iwp) :: j !<
	490	INTEGER(iwp) :: k !<
	491	INTEGER(iwp) :: nnyh !<
	492	INTEGER(iwp) :: nx !<
	493	INTEGER(iwp) :: ny !<
	494	INTEGER(iwp) :: omp_get_thread_num !<
	495	INTEGER(iwp) :: tn !<
[1212]	496
[1682]	497	REAL(wp) :: ddx2 !<
	498	REAL(wp) :: ddy2 !<
[1212]	499
[1682]	500	REAL(wp), DIMENSION(0:nx,1:nz) :: ar !<
	501	REAL(wp), DIMENSION(5,0:nx,0:nz-1) :: tri_for_1d !<
[1212]	502
	503
	504	nnyh = ( ny + 1 ) / 2
	505
	506	!
	507	!-- Define constant elements of the tridiagonal matrix.
	508	!-- The compiler on SX6 does loop exchange. If 0:nx is a high power of 2,
	509	!-- the exchanged loops create bank conflicts. The following directive
	510	!-- prohibits loop exchange and the loops perform much better.
	511	! tn = omp_get_thread_num()
	512	! WRITE( 120+tn, * ) '+++ id=',myid,' nx=',nx,' thread=', omp_get_thread_num()
	513	! CALL local_flush( 120+tn )
	514	!CDIR NOLOOPCHG
	515	DO k = 0, nz-1
	516	DO i = 0,nx
[1221]	517	tri_for_1d(2,i,k) = ddzu_pres(k+1) * ddzw(k+1)
	518	tri_for_1d(3,i,k) = ddzu_pres(k+2) * ddzw(k+1)
[1212]	519	ENDDO
	520	ENDDO
	521	! WRITE( 120+tn, * ) '+++ id=',myid,' end of first tridia loop thread=', omp_get_thread_num()
	522	! CALL local_flush( 120+tn )
	523
	524	IF ( j <= nnyh ) THEN
	525	#if defined( __intel11 )
[1221]	526	CALL maketri_1dd( j, tri_for_1d )
[1212]	527	#else
	528	CALL maketri_1dd( j )
	529	#endif
	530	ELSE
	531	#if defined( __intel11 )
[1221]	532	CALL maketri_1dd( ny+1-j, tri_for_1d )
[1212]	533	#else
	534	CALL maketri_1dd( ny+1-j )
	535	#endif
	536	ENDIF
	537	#if defined( __intel11 )
[1221]	538	CALL split_1dd( tri_for_1d )
[1212]	539	#else
	540	CALL split_1dd
	541	#endif
[1221]	542	CALL substi_1dd( ar, tri_for_1d )
[1212]	543
	544	CONTAINS
	545
[1682]	546
	547	!------------------------------------------------------------------------------!
	548	! Description:
	549	! ------------
	550	!> computes the i- and j-dependent component of the matrix
	551	!------------------------------------------------------------------------------!
[1212]	552	#if defined( __intel11 )
[1221]	553	SUBROUTINE maketri_1dd( j, tri_for_1d )
[1212]	554	#else
	555	SUBROUTINE maketri_1dd( j )
	556	#endif
	557
[1320]	558	USE constants, &
	559	ONLY: pi
[1212]	560
[1320]	561	USE kinds
	562
[1212]	563	IMPLICIT NONE
	564
[1682]	565	INTEGER(iwp) :: i !<
	566	INTEGER(iwp) :: j !<
	567	INTEGER(iwp) :: k !<
	568	INTEGER(iwp) :: nnxh !<
[1212]	569
[1682]	570	REAL(wp) :: a !<
	571	REAL(wp) :: c !<
[1212]	572
[1682]	573	REAL(wp), DIMENSION(0:nx) :: l !<
[1320]	574
[1212]	575	#if defined( __intel11 )
[1682]	576	REAL(wp), DIMENSION(5,0:nx,0:nz-1) :: tri_for_1d !<
[1212]	577	#endif
	578
	579
	580	nnxh = ( nx + 1 ) / 2
	581	!
	582	!-- Provide the tridiagonal matrix for solution of the Poisson equation in
	583	!-- Fourier space. The coefficients are computed following the method of
	584	!-- Schmidt et al. (DFVLR-Mitteilung 84-15), which departs from Stephan
	585	!-- Siano's original version by discretizing the Poisson equation,
	586	!-- before it is Fourier-transformed
	587	DO i = 0, nx
	588	IF ( i >= 0 .AND. i <= nnxh ) THEN
[1342]	589	l(i) = 2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * i ) / &
	590	REAL( nx+1, KIND=wp ) ) ) * ddx2 + &
	591	2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * j ) / &
	592	REAL( ny+1, KIND=wp ) ) ) * ddy2
[1212]	593	ELSE
[1342]	594	l(i) = 2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * ( nx+1-i ) ) / &
	595	REAL( nx+1, KIND=wp ) ) ) * ddx2 + &
	596	2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * j ) / &
	597	REAL( ny+1, KIND=wp ) ) ) * ddy2
[1212]	598	ENDIF
	599	ENDDO
	600
	601	DO k = 0, nz-1
	602	DO i = 0, nx
[1342]	603	a = -1.0_wp * ddzu_pres(k+2) * ddzw(k+1)
	604	c = -1.0_wp * ddzu_pres(k+1) * ddzw(k+1)
[1221]	605	tri_for_1d(1,i,k) = a + c - l(i)
[1212]	606	ENDDO
	607	ENDDO
	608	IF ( ibc_p_b == 1 ) THEN
	609	DO i = 0, nx
[1221]	610	tri_for_1d(1,i,0) = tri_for_1d(1,i,0) + tri_for_1d(2,i,0)
[1212]	611	ENDDO
	612	ENDIF
	613	IF ( ibc_p_t == 1 ) THEN
	614	DO i = 0, nx
[1221]	615	tri_for_1d(1,i,nz-1) = tri_for_1d(1,i,nz-1) + tri_for_1d(3,i,nz-1)
[1212]	616	ENDDO
	617	ENDIF
	618
	619	END SUBROUTINE maketri_1dd
	620
	621
[1682]	622	!------------------------------------------------------------------------------!
	623	! Description:
	624	! ------------
	625	!> Splitting of the tridiagonal matrix (Thomas algorithm)
	626	!------------------------------------------------------------------------------!
[1212]	627	#if defined( __intel11 )
[1221]	628	SUBROUTINE split_1dd( tri_for_1d )
[1212]	629	#else
	630	SUBROUTINE split_1dd
	631	#endif
	632
	633
	634	IMPLICIT NONE
	635
[1682]	636	INTEGER(iwp) :: i !<
	637	INTEGER(iwp) :: k !<
[1212]	638
	639	#if defined( __intel11 )
[1682]	640	REAL(wp), DIMENSION(5,0:nx,0:nz-1) :: tri_for_1d !<
[1212]	641	#endif
	642
	643
	644	!
	645	!-- Splitting
	646	DO i = 0, nx
[1221]	647	tri_for_1d(4,i,0) = tri_for_1d(1,i,0)
[1212]	648	ENDDO
	649	DO k = 1, nz-1
	650	DO i = 0, nx
[1221]	651	tri_for_1d(5,i,k) = tri_for_1d(2,i,k) / tri_for_1d(4,i,k-1)
	652	tri_for_1d(4,i,k) = tri_for_1d(1,i,k) - tri_for_1d(3,i,k-1) * tri_for_1d(5,i,k)
[1212]	653	ENDDO
	654	ENDDO
	655
	656	END SUBROUTINE split_1dd
	657
	658
	659	!------------------------------------------------------------------------------!
[1682]	660	! Description:
	661	! ------------
	662	!> Substitution (Forward and Backward) (Thomas algorithm)
[1212]	663	!------------------------------------------------------------------------------!
[1682]	664	SUBROUTINE substi_1dd( ar, tri_for_1d )
[1212]	665
[1682]	666
[1212]	667	IMPLICIT NONE
	668
[1682]	669	INTEGER(iwp) :: i !<
	670	INTEGER(iwp) :: k !<
[1212]	671
[1682]	672	REAL(wp), DIMENSION(0:nx,nz) :: ar !<
	673	REAL(wp), DIMENSION(0:nx,0:nz-1) :: ar1 !<
	674	REAL(wp), DIMENSION(5,0:nx,0:nz-1) :: tri_for_1d !<
[1212]	675
	676	!
	677	!-- Forward substitution
	678	DO i = 0, nx
	679	ar1(i,0) = ar(i,1)
	680	ENDDO
	681	DO k = 1, nz-1
	682	DO i = 0, nx
[1221]	683	ar1(i,k) = ar(i,k+1) - tri_for_1d(5,i,k) * ar1(i,k-1)
[1212]	684	ENDDO
	685	ENDDO
	686
	687	!
	688	!-- Backward substitution
	689	!-- Note, the add of 1.0E-20 in the denominator is due to avoid divisions
	690	!-- by zero appearing if the pressure bc is set to neumann at the top of
	691	!-- the model domain.
	692	DO i = 0, nx
[1342]	693	ar(i,nz) = ar1(i,nz-1) / ( tri_for_1d(4,i,nz-1) + 1.0E-20_wp )
[1212]	694	ENDDO
	695	DO k = nz-2, 0, -1
	696	DO i = 0, nx
[1221]	697	ar(i,k+1) = ( ar1(i,k) - tri_for_1d(3,i,k) * ar(i,k+2) ) &
	698	/ tri_for_1d(4,i,k)
[1212]	699	ENDDO
	700	ENDDO
	701
	702	!
	703	!-- Indices i=0, j=0 correspond to horizontally averaged pressure.
	704	!-- The respective values of ar should be zero at all k-levels if
	705	!-- acceleration of horizontally averaged vertical velocity is zero.
	706	IF ( ibc_p_b == 1 .AND. ibc_p_t == 1 ) THEN
	707	IF ( j == 0 ) THEN
	708	DO k = 1, nz
[1342]	709	ar(0,k) = 0.0_wp
[1212]	710	ENDDO
	711	ENDIF
	712	ENDIF
	713
	714	END SUBROUTINE substi_1dd
	715
	716	END SUBROUTINE tridia_1dd
	717
	718
	719	END MODULE tridia_solver

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |