Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |