Home

Context Navigation

source: palm/trunk/SOURCE/diffusion_s.f90 @ 1116

Last change on this file since 1116 was 1093, checked in by raasch, 12 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 17.7 KB

Rev	Line
[1]	1	MODULE diffusion_s_mod
	2
[1036]	3	!--------------------------------------------------------------------------------!
	4	! This file is part of PALM.
	5	!
	6	! PALM is free software: you can redistribute it and/or modify it under the terms
	7	! of the GNU General Public License as published by the Free Software Foundation,
	8	! either version 3 of the License, or (at your option) any later version.
	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	!
	17	! Copyright 1997-2012 Leibniz University Hannover
	18	!--------------------------------------------------------------------------------!
	19	!
[484]	20	! Current revisions:
[1001]	21	! ------------------
[1093]	22	!
[1]	23	!
	24	! Former revisions:
	25	! -----------------
[3]	26	! $Id: diffusion_s.f90 1093 2013-02-02 12:58:49Z hoffmann $
[39]	27	!
[1093]	28	! 1092 2013-02-02 11:24:22Z raasch
	29	! unused variables removed
	30	!
[1037]	31	! 1036 2012-10-22 13:43:42Z raasch
	32	! code put under GPL (PALM 3.9)
	33	!
[1017]	34	! 1015 2012-09-27 09:23:24Z raasch
	35	! accelerator version (*_acc) added
	36	!
[1011]	37	! 1010 2012-09-20 07:59:54Z raasch
	38	! cpp switch __nopointer added for pointer free version
	39	!
[1002]	40	! 1001 2012-09-13 14:08:46Z raasch
	41	! some arrays comunicated by module instead of parameter list
	42	!
[668]	43	! 667 2010-12-23 12:06:00Z suehring/gryschka
	44	! nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng
	45	!
[198]	46	! 183 2008-08-04 15:39:12Z letzel
	47	! bugfix: calculation of fluxes at vertical surfaces
	48	!
[139]	49	! 129 2007-10-30 12:12:24Z letzel
	50	! replace wall_heatflux by wall_s_flux that is now included in the parameter
	51	! list, bugfix for assignment of fluxes at walls
	52	!
[39]	53	! 20 2007-02-26 00:12:32Z raasch
	54	! Bugfix: ddzw dimensioned 1:nzt"+1"
	55	! Calculation extended for gridpoint nzt, fluxes can be given at top,
	56	! +s_flux_t in parameter list, s_flux renamed s_flux_b
	57	!
[3]	58	! RCS Log replace by Id keyword, revision history cleaned up
	59	!
[1]	60	! Revision 1.8 2006/02/23 10:34:17 raasch
	61	! nzb_2d replaced by nzb_s_outer in horizontal diffusion and by nzb_s_inner
	62	! or nzb_diff_s_inner, respectively, in vertical diffusion, prescribed surface
	63	! fluxes at vertically oriented topography
	64	!
	65	! Revision 1.1 2000/04/13 14:54:02 schroeter
	66	! Initial revision
	67	!
	68	!
	69	! Description:
	70	! ------------
	71	! Diffusion term of scalar quantities (temperature and water content)
	72	!------------------------------------------------------------------------------!
	73
	74	PRIVATE
[1015]	75	PUBLIC diffusion_s, diffusion_s_acc
[1]	76
	77	INTERFACE diffusion_s
	78	MODULE PROCEDURE diffusion_s
	79	MODULE PROCEDURE diffusion_s_ij
	80	END INTERFACE diffusion_s
	81
[1015]	82	INTERFACE diffusion_s_acc
	83	MODULE PROCEDURE diffusion_s_acc
	84	END INTERFACE diffusion_s_acc
	85
[1]	86	CONTAINS
	87
	88
	89	!------------------------------------------------------------------------------!
	90	! Call for all grid points
	91	!------------------------------------------------------------------------------!
[1001]	92	SUBROUTINE diffusion_s( s, s_flux_b, s_flux_t, wall_s_flux )
[1]	93
[1001]	94	USE arrays_3d
[1]	95	USE control_parameters
	96	USE grid_variables
	97	USE indices
	98
	99	IMPLICIT NONE
	100
	101	INTEGER :: i, j, k
[129]	102	REAL :: wall_s_flux(0:4)
[1001]	103	REAL, DIMENSION(nysg:nyng,nxlg:nxrg) :: s_flux_b, s_flux_t
[1010]	104	#if defined( __nopointer )
	105	REAL, DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: s
	106	#else
[1001]	107	REAL, DIMENSION(:,:,:), POINTER :: s
[1010]	108	#endif
[1]	109
	110	DO i = nxl, nxr
	111	DO j = nys,nyn
	112	!
	113	!-- Compute horizontal diffusion
[19]	114	DO k = nzb_s_outer(j,i)+1, nzt
[1]	115
	116	tend(k,j,i) = tend(k,j,i) &
	117	+ 0.5 * ( &
	118	( kh(k,j,i) + kh(k,j,i+1) ) * ( s(k,j,i+1)-s(k,j,i) ) &
	119	- ( kh(k,j,i) + kh(k,j,i-1) ) * ( s(k,j,i)-s(k,j,i-1) ) &
	120	) * ddx2 &
	121	+ 0.5 * ( &
	122	( kh(k,j,i) + kh(k,j+1,i) ) * ( s(k,j+1,i)-s(k,j,i) ) &
	123	- ( kh(k,j,i) + kh(k,j-1,i) ) * ( s(k,j,i)-s(k,j-1,i) ) &
	124	) * ddy2
	125	ENDDO
	126
	127	!
	128	!-- Apply prescribed horizontal wall heatflux where necessary
	129	IF ( ( wall_w_x(j,i) .NE. 0.0 ) .OR. ( wall_w_y(j,i) .NE. 0.0 ) ) &
	130	THEN
	131	DO k = nzb_s_inner(j,i)+1, nzb_s_outer(j,i)
	132
	133	tend(k,j,i) = tend(k,j,i) &
[183]	134	+ ( fwxp(j,i) * 0.5 * &
[1]	135	( kh(k,j,i) + kh(k,j,i+1) ) * ( s(k,j,i+1)-s(k,j,i) ) &
[129]	136	+ ( 1.0 - fwxp(j,i) ) * wall_s_flux(1) &
[183]	137	-fwxm(j,i) * 0.5 * &
[1]	138	( kh(k,j,i) + kh(k,j,i-1) ) * ( s(k,j,i)-s(k,j,i-1) ) &
[129]	139	+ ( 1.0 - fwxm(j,i) ) * wall_s_flux(2) &
[1]	140	) * ddx2 &
[183]	141	+ ( fwyp(j,i) * 0.5 * &
[1]	142	( kh(k,j,i) + kh(k,j+1,i) ) * ( s(k,j+1,i)-s(k,j,i) ) &
[129]	143	+ ( 1.0 - fwyp(j,i) ) * wall_s_flux(3) &
[183]	144	-fwym(j,i) * 0.5 * &
[1]	145	( kh(k,j,i) + kh(k,j-1,i) ) * ( s(k,j,i)-s(k,j-1,i) ) &
[129]	146	+ ( 1.0 - fwym(j,i) ) * wall_s_flux(4) &
[1]	147	) * ddy2
	148	ENDDO
	149	ENDIF
	150
	151	!
	152	!-- Compute vertical diffusion. In case that surface fluxes have been
[19]	153	!-- prescribed or computed at bottom and/or top, index k starts/ends at
	154	!-- nzb+2 or nzt-1, respectively.
	155	DO k = nzb_diff_s_inner(j,i), nzt_diff
[1]	156
	157	tend(k,j,i) = tend(k,j,i) &
	158	+ 0.5 * ( &
	159	( kh(k,j,i) + kh(k+1,j,i) ) * ( s(k+1,j,i)-s(k,j,i) ) * ddzu(k+1) &
	160	- ( kh(k,j,i) + kh(k-1,j,i) ) * ( s(k,j,i)-s(k-1,j,i) ) * ddzu(k) &
	161	) * ddzw(k)
	162	ENDDO
	163
	164	!
[19]	165	!-- Vertical diffusion at the first computational gridpoint along
[1]	166	!-- z-direction
	167	IF ( use_surface_fluxes ) THEN
	168
	169	k = nzb_s_inner(j,i)+1
	170
	171	tend(k,j,i) = tend(k,j,i) &
	172	+ ( 0.5 * ( kh(k,j,i)+kh(k+1,j,i) ) &
	173	* ( s(k+1,j,i)-s(k,j,i) ) &
	174	* ddzu(k+1) &
[19]	175	+ s_flux_b(j,i) &
[1]	176	) * ddzw(k)
	177
	178	ENDIF
	179
[19]	180	!
	181	!-- Vertical diffusion at the last computational gridpoint along
	182	!-- z-direction
	183	IF ( use_top_fluxes ) THEN
	184
	185	k = nzt
	186
	187	tend(k,j,i) = tend(k,j,i) &
	188	+ ( - s_flux_t(j,i) &
[20]	189	- 0.5 * ( kh(k-1,j,i)+kh(k,j,i) ) &
	190	* ( s(k,j,i)-s(k-1,j,i) ) &
	191	* ddzu(k) &
[19]	192	) * ddzw(k)
	193
	194	ENDIF
	195
[1]	196	ENDDO
	197	ENDDO
	198
	199	END SUBROUTINE diffusion_s
	200
	201
	202	!------------------------------------------------------------------------------!
[1015]	203	! Call for all grid points - accelerator version
	204	!------------------------------------------------------------------------------!
	205	SUBROUTINE diffusion_s_acc( s, s_flux_b, s_flux_t, wall_s_flux )
	206
	207	USE arrays_3d
	208	USE control_parameters
	209	USE grid_variables
	210	USE indices
	211
	212	IMPLICIT NONE
	213
	214	INTEGER :: i, j, k
	215	REAL :: wall_s_flux(0:4)
	216	REAL, DIMENSION(nysg:nyng,nxlg:nxrg) :: s_flux_b, s_flux_t
	217	#if defined( __nopointer )
	218	REAL, DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: s
	219	#else
	220	REAL, DIMENSION(:,:,:), POINTER :: s
	221	#endif
	222
	223	!$acc kernels present( ddzu, ddzw, fwxm, fwxp, fwym, fwyp, kh ) &
	224	!$acc present( nzb_diff_s_inner, nzb_s_inner, nzb_s_outer, s ) &
	225	!$acc present( s_flux_b, s_flux_t, tend, wall_s_flux ) &
	226	!$acc present( wall_w_x, wall_w_y )
	227	!$acc loop
	228	DO i = nxl, nxr
	229	DO j = nys,nyn
	230	!
	231	!-- Compute horizontal diffusion
	232	!$acc loop vector( 32 )
	233	DO k = 1, nzt
	234	IF ( k > nzb_s_outer(j,i) ) THEN
	235
	236	tend(k,j,i) = tend(k,j,i) &
	237	+ 0.5 * ( &
	238	( kh(k,j,i) + kh(k,j,i+1) ) * ( s(k,j,i+1)-s(k,j,i) ) &
	239	- ( kh(k,j,i) + kh(k,j,i-1) ) * ( s(k,j,i)-s(k,j,i-1) ) &
	240	) * ddx2 &
	241	+ 0.5 * ( &
	242	( kh(k,j,i) + kh(k,j+1,i) ) * ( s(k,j+1,i)-s(k,j,i) ) &
	243	- ( kh(k,j,i) + kh(k,j-1,i) ) * ( s(k,j,i)-s(k,j-1,i) ) &
	244	) * ddy2
	245	ENDIF
	246	ENDDO
	247
	248	!
	249	!-- Apply prescribed horizontal wall heatflux where necessary
	250	!$acc loop vector(32)
	251	DO k = 1, nzt
	252	IF ( k > nzb_s_inner(j,i) .AND. k <= nzb_s_outer(j,i) .AND. &
	253	( wall_w_x(j,i) /= 0.0 .OR. wall_w_y(j,i) /= 0.0 ) ) &
	254	THEN
	255	tend(k,j,i) = tend(k,j,i) &
	256	+ ( fwxp(j,i) * 0.5 * &
	257	( kh(k,j,i) + kh(k,j,i+1) ) * ( s(k,j,i+1)-s(k,j,i) ) &
	258	+ ( 1.0 - fwxp(j,i) ) * wall_s_flux(1) &
	259	-fwxm(j,i) * 0.5 * &
	260	( kh(k,j,i) + kh(k,j,i-1) ) * ( s(k,j,i)-s(k,j,i-1) ) &
	261	+ ( 1.0 - fwxm(j,i) ) * wall_s_flux(2) &
	262	) * ddx2 &
	263	+ ( fwyp(j,i) * 0.5 * &
	264	( kh(k,j,i) + kh(k,j+1,i) ) * ( s(k,j+1,i)-s(k,j,i) ) &
	265	+ ( 1.0 - fwyp(j,i) ) * wall_s_flux(3) &
	266	-fwym(j,i) * 0.5 * &
	267	( kh(k,j,i) + kh(k,j-1,i) ) * ( s(k,j,i)-s(k,j-1,i) ) &
	268	+ ( 1.0 - fwym(j,i) ) * wall_s_flux(4) &
	269	) * ddy2
	270	ENDIF
	271	ENDDO
	272
	273	!
	274	!-- Compute vertical diffusion. In case that surface fluxes have been
	275	!-- prescribed or computed at bottom and/or top, index k starts/ends at
	276	!-- nzb+2 or nzt-1, respectively.
	277	!$acc loop vector( 32 )
	278	DO k = 1, nzt_diff
	279	IF ( k >= nzb_diff_s_inner(j,i) ) THEN
	280	tend(k,j,i) = tend(k,j,i) &
	281	+ 0.5 * ( &
	282	( kh(k,j,i) + kh(k+1,j,i) ) * ( s(k+1,j,i)-s(k,j,i) ) * ddzu(k+1) &
	283	- ( kh(k,j,i) + kh(k-1,j,i) ) * ( s(k,j,i)-s(k-1,j,i) ) * ddzu(k) &
	284	) * ddzw(k)
	285	ENDIF
	286	ENDDO
	287
	288	!
	289	!-- Vertical diffusion at the first computational gridpoint along
	290	!-- z-direction
	291	!$acc loop vector( 32 )
	292	DO k = 1, nzt
	293	IF ( use_surface_fluxes .AND. k == nzb_s_inner(j,i)+1 ) THEN
	294	tend(k,j,i) = tend(k,j,i) &
	295	+ ( 0.5 * ( kh(k,j,i)+kh(k+1,j,i) ) &
	296	* ( s(k+1,j,i)-s(k,j,i) ) &
	297	* ddzu(k+1) &
	298	+ s_flux_b(j,i) &
	299	) * ddzw(k)
	300	ENDIF
	301
	302	!
	303	!-- Vertical diffusion at the last computational gridpoint along
	304	!-- z-direction
	305	IF ( use_top_fluxes .AND. k == nzt ) THEN
	306	tend(k,j,i) = tend(k,j,i) &
	307	+ ( - s_flux_t(j,i) &
	308	- 0.5 * ( kh(k-1,j,i)+kh(k,j,i) )&
	309	* ( s(k,j,i)-s(k-1,j,i) ) &
	310	* ddzu(k) &
	311	) * ddzw(k)
	312	ENDIF
	313	ENDDO
	314
	315	ENDDO
	316	ENDDO
	317	!$acc end kernels
	318
	319	END SUBROUTINE diffusion_s_acc
	320
	321
	322	!------------------------------------------------------------------------------!
[1]	323	! Call for grid point i,j
	324	!------------------------------------------------------------------------------!
[1001]	325	SUBROUTINE diffusion_s_ij( i, j, s, s_flux_b, s_flux_t, wall_s_flux )
[1]	326
[1001]	327	USE arrays_3d
[1]	328	USE control_parameters
	329	USE grid_variables
	330	USE indices
	331
	332	IMPLICIT NONE
	333
	334	INTEGER :: i, j, k
[129]	335	REAL :: wall_s_flux(0:4)
[1001]	336	REAL, DIMENSION(nysg:nyng,nxlg:nxrg) :: s_flux_b, s_flux_t
[1010]	337	#if defined( __nopointer )
	338	REAL, DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: s
	339	#else
[1001]	340	REAL, DIMENSION(:,:,:), POINTER :: s
[1010]	341	#endif
[1]	342
	343	!
	344	!-- Compute horizontal diffusion
[19]	345	DO k = nzb_s_outer(j,i)+1, nzt
[1]	346
	347	tend(k,j,i) = tend(k,j,i) &
	348	+ 0.5 * ( &
	349	( kh(k,j,i) + kh(k,j,i+1) ) * ( s(k,j,i+1)-s(k,j,i) ) &
	350	- ( kh(k,j,i) + kh(k,j,i-1) ) * ( s(k,j,i)-s(k,j,i-1) ) &
	351	) * ddx2 &
	352	+ 0.5 * ( &
	353	( kh(k,j,i) + kh(k,j+1,i) ) * ( s(k,j+1,i)-s(k,j,i) ) &
	354	- ( kh(k,j,i) + kh(k,j-1,i) ) * ( s(k,j,i)-s(k,j-1,i) ) &
	355	) * ddy2
	356	ENDDO
	357
	358	!
	359	!-- Apply prescribed horizontal wall heatflux where necessary
	360	IF ( ( wall_w_x(j,i) .NE. 0.0 ) .OR. ( wall_w_y(j,i) .NE. 0.0 ) ) &
	361	THEN
	362	DO k = nzb_s_inner(j,i)+1, nzb_s_outer(j,i)
	363
	364	tend(k,j,i) = tend(k,j,i) &
[183]	365	+ ( fwxp(j,i) * 0.5 * &
[1]	366	( kh(k,j,i) + kh(k,j,i+1) ) * ( s(k,j,i+1)-s(k,j,i) ) &
[129]	367	+ ( 1.0 - fwxp(j,i) ) * wall_s_flux(1) &
[183]	368	-fwxm(j,i) * 0.5 * &
[1]	369	( kh(k,j,i) + kh(k,j,i-1) ) * ( s(k,j,i)-s(k,j,i-1) ) &
[129]	370	+ ( 1.0 - fwxm(j,i) ) * wall_s_flux(2) &
[1]	371	) * ddx2 &
[183]	372	+ ( fwyp(j,i) * 0.5 * &
[1]	373	( kh(k,j,i) + kh(k,j+1,i) ) * ( s(k,j+1,i)-s(k,j,i) ) &
[129]	374	+ ( 1.0 - fwyp(j,i) ) * wall_s_flux(3) &
[183]	375	-fwym(j,i) * 0.5 * &
[1]	376	( kh(k,j,i) + kh(k,j-1,i) ) * ( s(k,j,i)-s(k,j-1,i) ) &
[129]	377	+ ( 1.0 - fwym(j,i) ) * wall_s_flux(4) &
[1]	378	) * ddy2
	379	ENDDO
	380	ENDIF
	381
	382	!
	383	!-- Compute vertical diffusion. In case that surface fluxes have been
[19]	384	!-- prescribed or computed at bottom and/or top, index k starts/ends at
	385	!-- nzb+2 or nzt-1, respectively.
	386	DO k = nzb_diff_s_inner(j,i), nzt_diff
[1]	387
	388	tend(k,j,i) = tend(k,j,i) &
	389	+ 0.5 * ( &
	390	( kh(k,j,i) + kh(k+1,j,i) ) * ( s(k+1,j,i)-s(k,j,i) ) * ddzu(k+1) &
	391	- ( kh(k,j,i) + kh(k-1,j,i) ) * ( s(k,j,i)-s(k-1,j,i) ) * ddzu(k) &
	392	) * ddzw(k)
	393	ENDDO
	394
	395	!
[19]	396	!-- Vertical diffusion at the first computational gridpoint along z-direction
[1]	397	IF ( use_surface_fluxes ) THEN
	398
	399	k = nzb_s_inner(j,i)+1
	400
[19]	401	tend(k,j,i) = tend(k,j,i) + ( 0.5 * ( kh(k,j,i)+kh(k+1,j,i) ) &
	402	* ( s(k+1,j,i)-s(k,j,i) ) &
	403	* ddzu(k+1) &
	404	+ s_flux_b(j,i) &
	405	) * ddzw(k)
[1]	406
	407	ENDIF
	408
[19]	409	!
	410	!-- Vertical diffusion at the last computational gridpoint along z-direction
	411	IF ( use_top_fluxes ) THEN
	412
	413	k = nzt
	414
	415	tend(k,j,i) = tend(k,j,i) + ( - s_flux_t(j,i) &
	416	- 0.5 * ( kh(k-1,j,i)+kh(k,j,i) ) &
	417	* ( s(k,j,i)-s(k-1,j,i) ) &
	418	* ddzu(k) &
	419	) * ddzw(k)
	420
	421	ENDIF
	422
[1]	423	END SUBROUTINE diffusion_s_ij
	424
	425	END MODULE diffusion_s_mod

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |