Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |