Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |