Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |