Home

Context Navigation

tridia_solver.f90 @ 1827

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

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

Context Navigation

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

Download in other formats: