Home

Context Navigation

source: palm/trunk/SOURCE/sor.f90 @ 2037

Last change on this file since 2037 was 2037, checked in by knoop, 7 years ago
Anelastic approximation implemented
Property svn:keywords set to `Id`
File size: 8.4 KB

Rev	Line
[1682]	1	!> @file sor.f90
[2000]	2	!------------------------------------------------------------------------------!
[1036]	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.
[1036]	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	!
[1818]	17	! Copyright 1997-2016 Leibniz Universitaet Hannover
[2000]	18	!------------------------------------------------------------------------------!
[1036]	19	!
[484]	20	! Current revisions:
[1]	21	! -----------------
[2037]	22	! Anelastic approximation implemented
[1354]	23	!
[1321]	24	! Former revisions:
	25	! -----------------
	26	! $Id: sor.f90 2037 2016-10-26 11:15:40Z knoop $
	27	!
[2001]	28	! 2000 2016-08-20 18:09:15Z knoop
	29	! Forced header and separation lines into 80 columns
	30	!
[1763]	31	! 1762 2016-02-25 12:31:13Z hellstea
	32	! Introduction of nested domain feature
	33	!
[1683]	34	! 1682 2015-10-07 23:56:08Z knoop
	35	! Code annotations made doxygen readable
	36	!
[1354]	37	! 1353 2014-04-08 15:21:23Z heinze
	38	! REAL constants provided with KIND-attribute
	39	!
[1321]	40	! 1320 2014-03-20 08:40:49Z raasch
[1320]	41	! ONLY-attribute added to USE-statements,
	42	! kind-parameters added to all INTEGER and REAL declaration statements,
	43	! kinds are defined in new module kinds,
	44	! old module precision_kind is removed,
	45	! revision history before 2012 removed,
	46	! comment fields (!:) to be used for variable explanations added to
	47	! all variable declaration statements
[1]	48	!
[1037]	49	! 1036 2012-10-22 13:43:42Z raasch
	50	! code put under GPL (PALM 3.9)
	51	!
[1]	52	! Revision 1.1 1997/08/11 06:25:56 raasch
	53	! Initial revision
	54	!
	55	!
	56	! Description:
	57	! ------------
[1682]	58	!> Solve the Poisson-equation with the SOR-Red/Black-scheme.
[3]	59	!------------------------------------------------------------------------------!
[1682]	60	SUBROUTINE sor( d, ddzu, ddzw, p )
[1]	61
[2037]	62	USE arrays_3d, &
	63	ONLY: rho_air, rho_air_zw
	64
[1320]	65	USE grid_variables, &
	66	ONLY: ddx2, ddy2
[1]	67
[1320]	68	USE indices, &
	69	ONLY: nbgp, nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nz, nzb, nzt
	70
	71	USE kinds
	72
	73	USE control_parameters, &
	74	ONLY: bc_lr_cyc, bc_ns_cyc, ibc_p_b, ibc_p_t, inflow_l, inflow_n, &
[1762]	75	inflow_r, inflow_s, nest_bound_l, nest_bound_n, nest_bound_r, &
	76	nest_bound_s, n_sor, omega_sor, outflow_l, outflow_n, &
[1320]	77	outflow_r, outflow_s
	78
[1]	79	IMPLICIT NONE
	80
[1682]	81	INTEGER(iwp) :: i !<
	82	INTEGER(iwp) :: j !<
	83	INTEGER(iwp) :: k !<
	84	INTEGER(iwp) :: n !<
	85	INTEGER(iwp) :: nxl1 !<
	86	INTEGER(iwp) :: nxl2 !<
	87	INTEGER(iwp) :: nys1 !<
	88	INTEGER(iwp) :: nys2 !<
[1]	89
[1682]	90	REAL(wp) :: ddzu(1:nz+1) !<
	91	REAL(wp) :: ddzw(1:nzt+1) !<
[1320]	92
[1682]	93	REAL(wp) :: d(nzb+1:nzt,nys:nyn,nxl:nxr) !<
	94	REAL(wp) :: p(nzb:nzt+1,nysg:nyng,nxlg:nxrg) !<
[1320]	95
[1682]	96	REAL(wp), DIMENSION(:), ALLOCATABLE :: f1 !<
	97	REAL(wp), DIMENSION(:), ALLOCATABLE :: f2 !<
	98	REAL(wp), DIMENSION(:), ALLOCATABLE :: f3 !<
[1320]	99
[1]	100	ALLOCATE( f1(1:nz), f2(1:nz), f3(1:nz) )
	101
	102	!
	103	!-- Compute pre-factors.
	104	DO k = 1, nz
[2037]	105	f2(k) = ddzu(k+1) * ddzw(k) * rho_air_zw(k)
	106	f3(k) = ddzu(k) * ddzw(k) * rho_air_zw(k-1)
	107	f1(k) = 2.0_wp * ( ddx2 + ddy2 ) * rho_air(k) + f2(k) + f3(k)
[1]	108	ENDDO
	109
	110	!
	111	!-- Limits for RED- and BLACK-part.
	112	IF ( MOD( nxl , 2 ) == 0 ) THEN
	113	nxl1 = nxl
	114	nxl2 = nxl + 1
	115	ELSE
	116	nxl1 = nxl + 1
	117	nxl2 = nxl
	118	ENDIF
	119	IF ( MOD( nys , 2 ) == 0 ) THEN
	120	nys1 = nys
	121	nys2 = nys + 1
	122	ELSE
	123	nys1 = nys + 1
	124	nys2 = nys
	125	ENDIF
	126
	127	DO n = 1, n_sor
	128
	129	!
	130	!-- RED-part
	131	DO i = nxl1, nxr, 2
	132	DO j = nys2, nyn, 2
	133	DO k = nzb+1, nzt
	134	p(k,j,i) = p(k,j,i) + omega_sor / f1(k) * ( &
[2037]	135	rho_air(k) * ddx2 * ( p(k,j,i+1) + p(k,j,i-1) ) + &
	136	rho_air(k) * ddy2 * ( p(k,j+1,i) + p(k,j-1,i) ) + &
	137	f2(k) * p(k+1,j,i) + &
	138	f3(k) * p(k-1,j,i) - &
	139	d(k,j,i) - &
	140	f1(k) * p(k,j,i) )
[1]	141	ENDDO
	142	ENDDO
	143	ENDDO
	144
	145	DO i = nxl2, nxr, 2
	146	DO j = nys1, nyn, 2
	147	DO k = nzb+1, nzt
[2037]	148	p(k,j,i) = p(k,j,i) + omega_sor / f1(k) * ( &
	149	rho_air(k) * ddx2 * ( p(k,j,i+1) + p(k,j,i-1) ) + &
	150	rho_air(k) * ddy2 * ( p(k,j+1,i) + p(k,j-1,i) ) + &
	151	f2(k) * p(k+1,j,i) + &
	152	f3(k) * p(k-1,j,i) - &
	153	d(k,j,i) - &
	154	f1(k) * p(k,j,i) )
[1]	155	ENDDO
	156	ENDDO
	157	ENDDO
	158
	159	!
	160	!-- Exchange of boundary values for p.
[667]	161	CALL exchange_horiz( p, nbgp )
[1]	162
	163	!
	164	!-- Horizontal (Neumann) boundary conditions in case of non-cyclic boundaries
[707]	165	IF ( .NOT. bc_lr_cyc ) THEN
[1762]	166	IF ( inflow_l .OR. outflow_l .OR. nest_bound_l ) p(:,:,nxl-1) = p(:,:,nxl)
	167	IF ( inflow_r .OR. outflow_r .OR. nest_bound_r ) p(:,:,nxr+1) = p(:,:,nxr)
[1]	168	ENDIF
[707]	169	IF ( .NOT. bc_ns_cyc ) THEN
[1762]	170	IF ( inflow_n .OR. outflow_n .OR. nest_bound_n ) p(:,nyn+1,:) = p(:,nyn,:)
	171	IF ( inflow_s .OR. outflow_s .OR. nest_bound_s ) p(:,nys-1,:) = p(:,nys,:)
[1]	172	ENDIF
	173
	174	!
	175	!-- BLACK-part
	176	DO i = nxl1, nxr, 2
	177	DO j = nys1, nyn, 2
	178	DO k = nzb+1, nzt
	179	p(k,j,i) = p(k,j,i) + omega_sor / f1(k) * ( &
[2037]	180	rho_air(k) * ddx2 * ( p(k,j,i+1) + p(k,j,i-1) ) + &
	181	rho_air(k) * ddy2 * ( p(k,j+1,i) + p(k,j-1,i) ) + &
	182	f2(k) * p(k+1,j,i) + &
	183	f3(k) * p(k-1,j,i) - &
	184	d(k,j,i) - &
	185	f1(k) * p(k,j,i) )
[1]	186	ENDDO
	187	ENDDO
	188	ENDDO
	189
	190	DO i = nxl2, nxr, 2
	191	DO j = nys2, nyn, 2
	192	DO k = nzb+1, nzt
	193	p(k,j,i) = p(k,j,i) + omega_sor / f1(k) * ( &
[2037]	194	rho_air(k) * ddx2 * ( p(k,j,i+1) + p(k,j,i-1) ) + &
	195	rho_air(k) * ddy2 * ( p(k,j+1,i) + p(k,j-1,i) ) + &
	196	f2(k) * p(k+1,j,i) + &
	197	f3(k) * p(k-1,j,i) - &
	198	d(k,j,i) - &
	199	f1(k) * p(k,j,i) )
[1]	200	ENDDO
	201	ENDDO
	202	ENDDO
	203
	204	!
	205	!-- Exchange of boundary values for p.
[667]	206	CALL exchange_horiz( p, nbgp )
[1]	207
	208	!
	209	!-- Boundary conditions top/bottom.
	210	!-- Bottom boundary
[667]	211	IF ( ibc_p_b == 1 ) THEN ! Neumann
[1]	212	p(nzb,:,:) = p(nzb+1,:,:)
[667]	213	ELSE ! Dirichlet
[1353]	214	p(nzb,:,:) = 0.0_wp
[1]	215	ENDIF
	216
	217	!
	218	!-- Top boundary
[667]	219	IF ( ibc_p_t == 1 ) THEN ! Neumann
[1]	220	p(nzt+1,:,:) = p(nzt,:,:)
[667]	221	ELSE ! Dirichlet
[1353]	222	p(nzt+1,:,:) = 0.0_wp
[1]	223	ENDIF
	224
	225	!
	226	!-- Horizontal (Neumann) boundary conditions in case of non-cyclic boundaries
[707]	227	IF ( .NOT. bc_lr_cyc ) THEN
[1762]	228	IF ( inflow_l .OR. outflow_l .OR. nest_bound_l ) p(:,:,nxl-1) = p(:,:,nxl)
	229	IF ( inflow_r .OR. outflow_r .OR. nest_bound_r ) p(:,:,nxr+1) = p(:,:,nxr)
[1]	230	ENDIF
[707]	231	IF ( .NOT. bc_ns_cyc ) THEN
[1762]	232	IF ( inflow_n .OR. outflow_n .OR. nest_bound_n ) p(:,nyn+1,:) = p(:,nyn,:)
	233	IF ( inflow_s .OR. outflow_s .OR. nest_bound_s ) p(:,nys-1,:) = p(:,nys,:)
[1]	234	ENDIF
	235
[667]	236
[1]	237	ENDDO
	238
	239	DEALLOCATE( f1, f2, f3 )
	240
	241	END SUBROUTINE sor

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |