Home

Context Navigation

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

Last change on this file since 4231 was 4182, checked in by scharf, 5 years ago
corrected "Former revisions" section minor formatting in "Former revisions" section added "Author" section
Property svn:keywords set to `Id`
File size: 20.9 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |