Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |