Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |