Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |