Home

Context Navigation

source: palm/trunk/SOURCE/tridia_solver_mod.f90 @ 2011

Last change on this file since 2011 was 2001, checked in by knoop, 8 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 22.4 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |