Home

Context Navigation

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

Last change on this file since 4648 was 4591, checked in by raasch, 4 years ago
files re-formatted to follow the PALM coding standard
Property svn:keywords set to `Id`
File size: 10.9 KB

Rev	Line
[1682]	1	!> @file sor.f90
[4591]	2	!--------------------------------------------------------------------------------------------------!
[2696]	3	! This file is part of the PALM model system.
[1036]	4	!
[4591]	5	! PALM is free software: you can redistribute it and/or modify it under the terms of the GNU General
	6	! Public License as published by the Free Software Foundation, either version 3 of the License, or
	7	! (at your option) any later version.
[1036]	8	!
[4591]	9	! PALM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the
	10	! implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
	11	! Public License for more details.
[1036]	12	!
[4591]	13	! You should have received a copy of the GNU General Public License along with PALM. If not, see
	14	! <http://www.gnu.org/licenses/>.
[1036]	15	!
[4360]	16	! Copyright 1997-2020 Leibniz Universitaet Hannover
[4591]	17	!--------------------------------------------------------------------------------------------------!
[1036]	18	!
[4591]	19	!
[484]	20	! Current revisions:
[1]	21	! -----------------
[4591]	22	!
	23	!
[1321]	24	! Former revisions:
	25	! -----------------
	26	! $Id: sor.f90 4591 2020-07-06 15:56:08Z raasch $
[4591]	27	! File re-formatted to follow the PALM coding standard
	28	!
	29	!
	30	! 4457 2020-03-11 14:20:43Z raasch
	31	! Use statement for exchange horiz added
	32	!
[4457]	33	! 4360 2020-01-07 11:25:50Z suehring
[4182]	34	! Corrected "Former revisions" section
[4591]	35	!
[4182]	36	! 3655 2019-01-07 16:51:22Z knoop
[3183]	37	! Rename variables in mesoscale-offline nesting mode
[1321]	38	!
[4182]	39	! Revision 1.1 1997/08/11 06:25:56 raasch
	40	! Initial revision
	41	!
[4591]	42	!--------------------------------------------------------------------------------------------------!
[1]	43	! Description:
	44	! ------------
[1682]	45	!> Solve the Poisson-equation with the SOR-Red/Black-scheme.
[4591]	46	!--------------------------------------------------------------------------------------------------!
[1682]	47	SUBROUTINE sor( d, ddzu, ddzw, p )
[1]	48
[4591]	49	USE arrays_3d, &
	50	ONLY: rho_air, &
	51	rho_air_zw
[2037]	52
[4591]	53	USE control_parameters, &
	54	ONLY: bc_dirichlet_l, &
	55	bc_dirichlet_n, &
	56	bc_dirichlet_r, &
	57	bc_dirichlet_s, &
	58	bc_lr_cyc, &
	59	bc_ns_cyc, &
	60	bc_radiation_l, &
	61	bc_radiation_n, &
	62	bc_radiation_r, &
	63	bc_radiation_s, &
	64	ibc_p_b, &
	65	ibc_p_t, &
	66	n_sor, &
	67	omega_sor
[4457]	68
[4591]	69	USE exchange_horiz_mod, &
[4457]	70	ONLY: exchange_horiz
	71
[4591]	72	USE grid_variables, &
	73	ONLY: ddx2, &
	74	ddy2
[1]	75
[4591]	76	USE indices, &
	77	ONLY: nbgp, &
	78	nxl, &
	79	nxlg, &
	80	nxr, &
	81	nxrg, &
	82	nyn, &
	83	nyng, &
	84	nys, &
	85	nysg, &
	86	nz, &
	87	nzb, &
	88	nzt
[1320]	89
	90	USE kinds
	91
[1]	92	IMPLICIT NONE
	93
[4591]	94	INTEGER(iwp) :: i !<
	95	INTEGER(iwp) :: j !<
	96	INTEGER(iwp) :: k !<
	97	INTEGER(iwp) :: n !<
	98	INTEGER(iwp) :: nxl1 !<
	99	INTEGER(iwp) :: nxl2 !<
	100	INTEGER(iwp) :: nys1 !<
	101	INTEGER(iwp) :: nys2 !<
[1]	102
[4591]	103	REAL(wp) :: ddzu(1:nz+1) !<
	104	REAL(wp) :: ddzw(1:nzt+1) !<
[1320]	105
[4591]	106	REAL(wp) :: d(nzb+1:nzt,nys:nyn,nxl:nxr) !<
	107	REAL(wp) :: p(nzb:nzt+1,nysg:nyng,nxlg:nxrg) !<
[1320]	108
[4591]	109	REAL(wp), DIMENSION(:), ALLOCATABLE :: f1 !<
	110	REAL(wp), DIMENSION(:), ALLOCATABLE :: f2 !<
	111	REAL(wp), DIMENSION(:), ALLOCATABLE :: f3 !<
[1320]	112
[1]	113	ALLOCATE( f1(1:nz), f2(1:nz), f3(1:nz) )
	114
	115	!
	116	!-- Compute pre-factors.
	117	DO k = 1, nz
[2037]	118	f2(k) = ddzu(k+1) * ddzw(k) * rho_air_zw(k)
	119	f3(k) = ddzu(k) * ddzw(k) * rho_air_zw(k-1)
	120	f1(k) = 2.0_wp * ( ddx2 + ddy2 ) * rho_air(k) + f2(k) + f3(k)
[1]	121	ENDDO
	122
	123	!
	124	!-- Limits for RED- and BLACK-part.
	125	IF ( MOD( nxl , 2 ) == 0 ) THEN
	126	nxl1 = nxl
	127	nxl2 = nxl + 1
	128	ELSE
	129	nxl1 = nxl + 1
	130	nxl2 = nxl
	131	ENDIF
	132	IF ( MOD( nys , 2 ) == 0 ) THEN
	133	nys1 = nys
	134	nys2 = nys + 1
	135	ELSE
	136	nys1 = nys + 1
	137	nys2 = nys
	138	ENDIF
	139
	140	DO n = 1, n_sor
	141
	142	!
	143	!-- RED-part
	144	DO i = nxl1, nxr, 2
	145	DO j = nys2, nyn, 2
	146	DO k = nzb+1, nzt
[4591]	147	p(k,j,i) = p(k,j,i) + omega_sor / f1(k) * ( &
	148	rho_air(k) * ddx2 * ( p(k,j,i+1) + p(k,j,i-1) ) + &
	149	rho_air(k) * ddy2 * ( p(k,j+1,i) + p(k,j-1,i) ) + &
	150	f2(k) * p(k+1,j,i) + &
	151	f3(k) * p(k-1,j,i) - &
	152	d(k,j,i) - &
	153	f1(k) * p(k,j,i) )
[1]	154	ENDDO
	155	ENDDO
	156	ENDDO
	157
	158	DO i = nxl2, nxr, 2
	159	DO j = nys1, nyn, 2
	160	DO k = nzb+1, nzt
[4591]	161	p(k,j,i) = p(k,j,i) + omega_sor / f1(k) * ( &
	162	rho_air(k) * ddx2 * ( p(k,j,i+1) + p(k,j,i-1) ) + &
	163	rho_air(k) * ddy2 * ( p(k,j+1,i) + p(k,j-1,i) ) + &
	164	f2(k) * p(k+1,j,i) + &
	165	f3(k) * p(k-1,j,i) - &
	166	d(k,j,i) - &
	167	f1(k) * p(k,j,i) )
[1]	168	ENDDO
	169	ENDDO
	170	ENDDO
	171
	172	!
	173	!-- Exchange of boundary values for p.
[667]	174	CALL exchange_horiz( p, nbgp )
[1]	175
[4591]	176
[1]	177	!
	178	!-- Horizontal (Neumann) boundary conditions in case of non-cyclic boundaries
[707]	179	IF ( .NOT. bc_lr_cyc ) THEN
[3182]	180	IF ( bc_dirichlet_l .OR. bc_radiation_l ) p(:,:,nxl-1) = p(:,:,nxl)
	181	IF ( bc_dirichlet_r .OR. bc_radiation_r ) p(:,:,nxr+1) = p(:,:,nxr)
[1]	182	ENDIF
[707]	183	IF ( .NOT. bc_ns_cyc ) THEN
[3182]	184	IF ( bc_dirichlet_n .OR. bc_radiation_n ) p(:,nyn+1,:) = p(:,nyn,:)
	185	IF ( bc_dirichlet_s .OR. bc_radiation_s ) p(:,nys-1,:) = p(:,nys,:)
[1]	186	ENDIF
	187
	188	!
	189	!-- BLACK-part
	190	DO i = nxl1, nxr, 2
	191	DO j = nys1, nyn, 2
	192	DO k = nzb+1, nzt
[4591]	193	p(k,j,i) = p(k,j,i) + omega_sor / f1(k) * ( &
	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	DO i = nxl2, nxr, 2
	205	DO j = nys2, nyn, 2
	206	DO k = nzb+1, nzt
[4591]	207	p(k,j,i) = p(k,j,i) + omega_sor / f1(k) * ( &
	208	rho_air(k) * ddx2 * ( p(k,j,i+1) + p(k,j,i-1) ) + &
	209	rho_air(k) * ddy2 * ( p(k,j+1,i) + p(k,j-1,i) ) + &
	210	f2(k) * p(k+1,j,i) + &
	211	f3(k) * p(k-1,j,i) - &
	212	d(k,j,i) - &
	213	f1(k) * p(k,j,i) )
[1]	214	ENDDO
	215	ENDDO
	216	ENDDO
	217
	218	!
	219	!-- Exchange of boundary values for p.
[667]	220	CALL exchange_horiz( p, nbgp )
[1]	221
	222	!
	223	!-- Boundary conditions top/bottom.
	224	!-- Bottom boundary
[4591]	225	IF ( ibc_p_b == 1 ) THEN ! Neumann
[1]	226	p(nzb,:,:) = p(nzb+1,:,:)
[4591]	227	ELSE ! Dirichlet
[1353]	228	p(nzb,:,:) = 0.0_wp
[1]	229	ENDIF
	230
	231	!
	232	!-- Top boundary
[4591]	233	IF ( ibc_p_t == 1 ) THEN ! Neumann
[1]	234	p(nzt+1,:,:) = p(nzt,:,:)
[4591]	235	ELSE ! Dirichlet
[1353]	236	p(nzt+1,:,:) = 0.0_wp
[1]	237	ENDIF
	238
	239	!
	240	!-- Horizontal (Neumann) boundary conditions in case of non-cyclic boundaries
[707]	241	IF ( .NOT. bc_lr_cyc ) THEN
[3182]	242	IF ( bc_dirichlet_l .OR. bc_radiation_l ) p(:,:,nxl-1) = p(:,:,nxl)
	243	IF ( bc_dirichlet_r .OR. bc_radiation_r ) p(:,:,nxr+1) = p(:,:,nxr)
[1]	244	ENDIF
[707]	245	IF ( .NOT. bc_ns_cyc ) THEN
[3182]	246	IF ( bc_dirichlet_n .OR. bc_radiation_n ) p(:,nyn+1,:) = p(:,nyn,:)
	247	IF ( bc_dirichlet_s .OR. bc_radiation_s ) p(:,nys-1,:) = p(:,nys,:)
[1]	248	ENDIF
	249
[667]	250
[1]	251	ENDDO
	252
	253	DEALLOCATE( f1, f2, f3 )
	254
[4457]	255	END SUBROUTINE sor

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |