Home

Context Navigation

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

Last change on this file since 1350 was 1343, checked in by kanani, 11 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 22.3 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |