Home

Context Navigation

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

Last change on this file since 1935 was 1851, checked in by maronga, 9 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 22.3 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |