Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |