Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |