Home

Context Navigation

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

Last change on this file since 4476 was 4457, checked in by raasch, 5 years ago
ghost point exchange modularized, bugfix for wrong 2d-exchange
Property svn:keywords set to `Id`
File size: 7.9 KB

Rev	Line
[1682]	1	!> @file sor.f90
[2000]	2	!------------------------------------------------------------------------------!
[2696]	3	! This file is part of the PALM model system.
[1036]	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	!
[4360]	17	! Copyright 1997-2020 Leibniz Universitaet Hannover
[2000]	18	!------------------------------------------------------------------------------!
[1036]	19	!
[484]	20	! Current revisions:
[1]	21	! -----------------
[1354]	22	!
[3183]	23	!
[1321]	24	! Former revisions:
	25	! -----------------
	26	! $Id: sor.f90 4457 2020-03-11 14:20:43Z maronga $
[4457]	27	! use statement for exchange horiz added
	28	!
	29	! 4360 2020-01-07 11:25:50Z suehring
[4182]	30	! Corrected "Former revisions" section
	31	!
	32	! 3655 2019-01-07 16:51:22Z knoop
[3183]	33	! Rename variables in mesoscale-offline nesting mode
[1321]	34	!
[4182]	35	! Revision 1.1 1997/08/11 06:25:56 raasch
	36	! Initial revision
	37	!
	38	!
[1]	39	! Description:
	40	! ------------
[1682]	41	!> Solve the Poisson-equation with the SOR-Red/Black-scheme.
[3]	42	!------------------------------------------------------------------------------!
[1682]	43	SUBROUTINE sor( d, ddzu, ddzw, p )
[1]	44
[2037]	45	USE arrays_3d, &
	46	ONLY: rho_air, rho_air_zw
	47
[4457]	48	USE control_parameters, &
	49	ONLY: bc_dirichlet_l, bc_dirichlet_n, bc_dirichlet_r, &
	50	bc_dirichlet_s, bc_lr_cyc, bc_ns_cyc, bc_radiation_l, &
	51	bc_radiation_n, bc_radiation_r, bc_radiation_s, ibc_p_b, &
	52	ibc_p_t, n_sor, omega_sor
	53
	54	USE exchange_horiz_mod, &
	55	ONLY: exchange_horiz
	56
[1320]	57	USE grid_variables, &
	58	ONLY: ddx2, ddy2
[1]	59
[1320]	60	USE indices, &
	61	ONLY: nbgp, nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nz, nzb, nzt
	62
	63	USE kinds
	64
[1]	65	IMPLICIT NONE
	66
[1682]	67	INTEGER(iwp) :: i !<
	68	INTEGER(iwp) :: j !<
	69	INTEGER(iwp) :: k !<
	70	INTEGER(iwp) :: n !<
	71	INTEGER(iwp) :: nxl1 !<
	72	INTEGER(iwp) :: nxl2 !<
	73	INTEGER(iwp) :: nys1 !<
	74	INTEGER(iwp) :: nys2 !<
[1]	75
[1682]	76	REAL(wp) :: ddzu(1:nz+1) !<
	77	REAL(wp) :: ddzw(1:nzt+1) !<
[1320]	78
[1682]	79	REAL(wp) :: d(nzb+1:nzt,nys:nyn,nxl:nxr) !<
	80	REAL(wp) :: p(nzb:nzt+1,nysg:nyng,nxlg:nxrg) !<
[1320]	81
[1682]	82	REAL(wp), DIMENSION(:), ALLOCATABLE :: f1 !<
	83	REAL(wp), DIMENSION(:), ALLOCATABLE :: f2 !<
	84	REAL(wp), DIMENSION(:), ALLOCATABLE :: f3 !<
[1320]	85
[1]	86	ALLOCATE( f1(1:nz), f2(1:nz), f3(1:nz) )
	87
	88	!
	89	!-- Compute pre-factors.
	90	DO k = 1, nz
[2037]	91	f2(k) = ddzu(k+1) * ddzw(k) * rho_air_zw(k)
	92	f3(k) = ddzu(k) * ddzw(k) * rho_air_zw(k-1)
	93	f1(k) = 2.0_wp * ( ddx2 + ddy2 ) * rho_air(k) + f2(k) + f3(k)
[1]	94	ENDDO
	95
	96	!
	97	!-- Limits for RED- and BLACK-part.
	98	IF ( MOD( nxl , 2 ) == 0 ) THEN
	99	nxl1 = nxl
	100	nxl2 = nxl + 1
	101	ELSE
	102	nxl1 = nxl + 1
	103	nxl2 = nxl
	104	ENDIF
	105	IF ( MOD( nys , 2 ) == 0 ) THEN
	106	nys1 = nys
	107	nys2 = nys + 1
	108	ELSE
	109	nys1 = nys + 1
	110	nys2 = nys
	111	ENDIF
	112
	113	DO n = 1, n_sor
	114
	115	!
	116	!-- RED-part
	117	DO i = nxl1, nxr, 2
	118	DO j = nys2, nyn, 2
	119	DO k = nzb+1, nzt
	120	p(k,j,i) = p(k,j,i) + omega_sor / f1(k) * ( &
[2037]	121	rho_air(k) * ddx2 * ( p(k,j,i+1) + p(k,j,i-1) ) + &
	122	rho_air(k) * ddy2 * ( p(k,j+1,i) + p(k,j-1,i) ) + &
	123	f2(k) * p(k+1,j,i) + &
	124	f3(k) * p(k-1,j,i) - &
	125	d(k,j,i) - &
	126	f1(k) * p(k,j,i) )
[1]	127	ENDDO
	128	ENDDO
	129	ENDDO
	130
	131	DO i = nxl2, nxr, 2
	132	DO j = nys1, nyn, 2
	133	DO k = nzb+1, nzt
[2037]	134	p(k,j,i) = p(k,j,i) + omega_sor / f1(k) * ( &
	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	!
	146	!-- Exchange of boundary values for p.
[667]	147	CALL exchange_horiz( p, nbgp )
[1]	148
	149	!
	150	!-- Horizontal (Neumann) boundary conditions in case of non-cyclic boundaries
[707]	151	IF ( .NOT. bc_lr_cyc ) THEN
[3182]	152	IF ( bc_dirichlet_l .OR. bc_radiation_l ) p(:,:,nxl-1) = p(:,:,nxl)
	153	IF ( bc_dirichlet_r .OR. bc_radiation_r ) p(:,:,nxr+1) = p(:,:,nxr)
[1]	154	ENDIF
[707]	155	IF ( .NOT. bc_ns_cyc ) THEN
[3182]	156	IF ( bc_dirichlet_n .OR. bc_radiation_n ) p(:,nyn+1,:) = p(:,nyn,:)
	157	IF ( bc_dirichlet_s .OR. bc_radiation_s ) p(:,nys-1,:) = p(:,nys,:)
[1]	158	ENDIF
	159
	160	!
	161	!-- BLACK-part
	162	DO i = nxl1, nxr, 2
	163	DO j = nys1, nyn, 2
	164	DO k = nzb+1, nzt
	165	p(k,j,i) = p(k,j,i) + omega_sor / f1(k) * ( &
[2037]	166	rho_air(k) * ddx2 * ( p(k,j,i+1) + p(k,j,i-1) ) + &
	167	rho_air(k) * ddy2 * ( p(k,j+1,i) + p(k,j-1,i) ) + &
	168	f2(k) * p(k+1,j,i) + &
	169	f3(k) * p(k-1,j,i) - &
	170	d(k,j,i) - &
	171	f1(k) * p(k,j,i) )
[1]	172	ENDDO
	173	ENDDO
	174	ENDDO
	175
	176	DO i = nxl2, nxr, 2
	177	DO j = nys2, 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	!
	191	!-- Exchange of boundary values for p.
[667]	192	CALL exchange_horiz( p, nbgp )
[1]	193
	194	!
	195	!-- Boundary conditions top/bottom.
	196	!-- Bottom boundary
[667]	197	IF ( ibc_p_b == 1 ) THEN ! Neumann
[1]	198	p(nzb,:,:) = p(nzb+1,:,:)
[667]	199	ELSE ! Dirichlet
[1353]	200	p(nzb,:,:) = 0.0_wp
[1]	201	ENDIF
	202
	203	!
	204	!-- Top boundary
[667]	205	IF ( ibc_p_t == 1 ) THEN ! Neumann
[1]	206	p(nzt+1,:,:) = p(nzt,:,:)
[667]	207	ELSE ! Dirichlet
[1353]	208	p(nzt+1,:,:) = 0.0_wp
[1]	209	ENDIF
	210
	211	!
	212	!-- Horizontal (Neumann) boundary conditions in case of non-cyclic boundaries
[707]	213	IF ( .NOT. bc_lr_cyc ) THEN
[3182]	214	IF ( bc_dirichlet_l .OR. bc_radiation_l ) p(:,:,nxl-1) = p(:,:,nxl)
	215	IF ( bc_dirichlet_r .OR. bc_radiation_r ) p(:,:,nxr+1) = p(:,:,nxr)
[1]	216	ENDIF
[707]	217	IF ( .NOT. bc_ns_cyc ) THEN
[3182]	218	IF ( bc_dirichlet_n .OR. bc_radiation_n ) p(:,nyn+1,:) = p(:,nyn,:)
	219	IF ( bc_dirichlet_s .OR. bc_radiation_s ) p(:,nys-1,:) = p(:,nys,:)
[1]	220	ENDIF
	221
[667]	222
[1]	223	ENDDO
	224
	225	DEALLOCATE( f1, f2, f3 )
	226
[4457]	227	END SUBROUTINE sor

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |