Home

Context Navigation

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

Last change on this file since 3182 was 3182, checked in by suehring, 6 years ago
New Inifor features: grid stretching, improved command-interface, support start dates in different formats in both YYYYMMDD and YYYYMMDDHH, Ability to manually control input file prefixes (--radiation-prefix, --soil-preifx, --flow-prefix, --soilmoisture-prefix) for compatiblity with DWD forcast naming scheme; GNU-style short and long option; Prepared output of large-scale forcing profiles (no computation yet); Added preprocessor flag netcdf4 to switch output format between netCDF 3 and 4; Updated netCDF variable names and attributes to comply with PIDS v1.9; Inifor bugfixes: Improved compatibility with older Intel Intel compilers by avoiding implicit array allocation; Added origin_lon/_lat values and correct reference time in dynamic driver global attributes; corresponding PALM changes: adjustments to revised Inifor; variables names in dynamic driver adjusted; enable geostrophic forcing also in offline nested mode; variable names in LES-LES and COSMO offline nesting changed; lateral boundary flags for nesting, in- and outflow conditions renamed
Property svn:keywords set to `Id`
File size: 8.8 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	!
[2718]	17	! Copyright 1997-2018 Leibniz Universitaet Hannover
[2000]	18	!------------------------------------------------------------------------------!
[1036]	19	!
[484]	20	! Current revisions:
[1]	21	! -----------------
[3182]	22	! Rename variables in mesoscale-offline nesting mode
[1354]	23	!
[1321]	24	! Former revisions:
	25	! -----------------
	26	! $Id: sor.f90 3182 2018-07-27 13:36:03Z suehring $
[2716]	27	! Corrected "Former revisions" section
	28	!
[3182]	29	! 2718 2018-01-02 08:49:38Z maronga
	30	! Corrected "Former revisions" section
	31	!
[2716]	32	! 2696 2017-12-14 17:12:51Z kanani
	33	! - Change in file header (GPL part)
[2696]	34	! - Large-scale forcing implemented (MS)
	35	!
	36	! 2101 2017-01-05 16:42:31Z suehring
[1321]	37	!
[2038]	38	! 2037 2016-10-26 11:15:40Z knoop
	39	! Anelastic approximation implemented
	40	!
[2001]	41	! 2000 2016-08-20 18:09:15Z knoop
	42	! Forced header and separation lines into 80 columns
	43	!
[1763]	44	! 1762 2016-02-25 12:31:13Z hellstea
	45	! Introduction of nested domain feature
	46	!
[1683]	47	! 1682 2015-10-07 23:56:08Z knoop
	48	! Code annotations made doxygen readable
	49	!
[1354]	50	! 1353 2014-04-08 15:21:23Z heinze
	51	! REAL constants provided with KIND-attribute
	52	!
[1321]	53	! 1320 2014-03-20 08:40:49Z raasch
[1320]	54	! ONLY-attribute added to USE-statements,
	55	! kind-parameters added to all INTEGER and REAL declaration statements,
	56	! kinds are defined in new module kinds,
	57	! old module precision_kind is removed,
	58	! revision history before 2012 removed,
	59	! comment fields (!:) to be used for variable explanations added to
	60	! all variable declaration statements
[1]	61	!
[1037]	62	! 1036 2012-10-22 13:43:42Z raasch
	63	! code put under GPL (PALM 3.9)
	64	!
[1]	65	! Revision 1.1 1997/08/11 06:25:56 raasch
	66	! Initial revision
	67	!
	68	!
	69	! Description:
	70	! ------------
[1682]	71	!> Solve the Poisson-equation with the SOR-Red/Black-scheme.
[3]	72	!------------------------------------------------------------------------------!
[1682]	73	SUBROUTINE sor( d, ddzu, ddzw, p )
[1]	74
[2037]	75	USE arrays_3d, &
	76	ONLY: rho_air, rho_air_zw
	77
[1320]	78	USE grid_variables, &
	79	ONLY: ddx2, ddy2
[1]	80
[1320]	81	USE indices, &
	82	ONLY: nbgp, nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nz, nzb, nzt
	83
	84	USE kinds
	85
	86	USE control_parameters, &
[3182]	87	ONLY: bc_dirichlet_l, bc_dirichlet_n, bc_dirichlet_r, &
	88	bc_dirichlet_s, bc_lr_cyc, bc_ns_cyc, bc_radiation_l, &
	89	bc_radiation_n, bc_radiation_r, bc_radiation_s, ibc_p_b, &
	90	ibc_p_t, n_sor, omega_sor
[1320]	91
[1]	92	IMPLICIT NONE
	93
[1682]	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
[1682]	103	REAL(wp) :: ddzu(1:nz+1) !<
	104	REAL(wp) :: ddzw(1:nzt+1) !<
[1320]	105
[1682]	106	REAL(wp) :: d(nzb+1:nzt,nys:nyn,nxl:nxr) !<
	107	REAL(wp) :: p(nzb:nzt+1,nysg:nyng,nxlg:nxrg) !<
[1320]	108
[1682]	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
	147	p(k,j,i) = p(k,j,i) + omega_sor / f1(k) * ( &
[2037]	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
[2037]	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
	176	!
	177	!-- Horizontal (Neumann) boundary conditions in case of non-cyclic boundaries
[707]	178	IF ( .NOT. bc_lr_cyc ) THEN
[3182]	179	IF ( bc_dirichlet_l .OR. bc_radiation_l ) p(:,:,nxl-1) = p(:,:,nxl)
	180	IF ( bc_dirichlet_r .OR. bc_radiation_r ) p(:,:,nxr+1) = p(:,:,nxr)
[1]	181	ENDIF
[707]	182	IF ( .NOT. bc_ns_cyc ) THEN
[3182]	183	IF ( bc_dirichlet_n .OR. bc_radiation_n ) p(:,nyn+1,:) = p(:,nyn,:)
	184	IF ( bc_dirichlet_s .OR. bc_radiation_s ) p(:,nys-1,:) = p(:,nys,:)
[1]	185	ENDIF
	186
	187	!
	188	!-- BLACK-part
	189	DO i = nxl1, nxr, 2
	190	DO j = nys1, nyn, 2
	191	DO k = nzb+1, nzt
	192	p(k,j,i) = p(k,j,i) + omega_sor / f1(k) * ( &
[2037]	193	rho_air(k) * ddx2 * ( p(k,j,i+1) + p(k,j,i-1) ) + &
	194	rho_air(k) * ddy2 * ( p(k,j+1,i) + p(k,j-1,i) ) + &
	195	f2(k) * p(k+1,j,i) + &
	196	f3(k) * p(k-1,j,i) - &
	197	d(k,j,i) - &
	198	f1(k) * p(k,j,i) )
[1]	199	ENDDO
	200	ENDDO
	201	ENDDO
	202
	203	DO i = nxl2, nxr, 2
	204	DO j = nys2, nyn, 2
	205	DO k = nzb+1, nzt
	206	p(k,j,i) = p(k,j,i) + omega_sor / f1(k) * ( &
[2037]	207	rho_air(k) * ddx2 * ( p(k,j,i+1) + p(k,j,i-1) ) + &
	208	rho_air(k) * ddy2 * ( p(k,j+1,i) + p(k,j-1,i) ) + &
	209	f2(k) * p(k+1,j,i) + &
	210	f3(k) * p(k-1,j,i) - &
	211	d(k,j,i) - &
	212	f1(k) * p(k,j,i) )
[1]	213	ENDDO
	214	ENDDO
	215	ENDDO
	216
	217	!
	218	!-- Exchange of boundary values for p.
[667]	219	CALL exchange_horiz( p, nbgp )
[1]	220
	221	!
	222	!-- Boundary conditions top/bottom.
	223	!-- Bottom boundary
[667]	224	IF ( ibc_p_b == 1 ) THEN ! Neumann
[1]	225	p(nzb,:,:) = p(nzb+1,:,:)
[667]	226	ELSE ! Dirichlet
[1353]	227	p(nzb,:,:) = 0.0_wp
[1]	228	ENDIF
	229
	230	!
	231	!-- Top boundary
[667]	232	IF ( ibc_p_t == 1 ) THEN ! Neumann
[1]	233	p(nzt+1,:,:) = p(nzt,:,:)
[667]	234	ELSE ! Dirichlet
[1353]	235	p(nzt+1,:,:) = 0.0_wp
[1]	236	ENDIF
	237
	238	!
	239	!-- Horizontal (Neumann) boundary conditions in case of non-cyclic boundaries
[707]	240	IF ( .NOT. bc_lr_cyc ) THEN
[3182]	241	IF ( bc_dirichlet_l .OR. bc_radiation_l ) p(:,:,nxl-1) = p(:,:,nxl)
	242	IF ( bc_dirichlet_r .OR. bc_radiation_r ) p(:,:,nxr+1) = p(:,:,nxr)
[1]	243	ENDIF
[707]	244	IF ( .NOT. bc_ns_cyc ) THEN
[3182]	245	IF ( bc_dirichlet_n .OR. bc_radiation_n ) p(:,nyn+1,:) = p(:,nyn,:)
	246	IF ( bc_dirichlet_s .OR. bc_radiation_s ) p(:,nys-1,:) = p(:,nys,:)
[1]	247	ENDIF
	248
[667]	249
[1]	250	ENDDO
	251
	252	DEALLOCATE( f1, f2, f3 )
	253
	254	END SUBROUTINE sor

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |