Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |