Home

Context Navigation

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

Last change on this file since 1039 was 1037, 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	! ------------------
[1]	22	!
[1017]	23	!
[1]	24	! Former revisions:
	25	! -----------------
[3]	26	! $Id: diffusion_s.f90 1037 2012-10-22 14:10:22Z raasch $
[39]	27	!
[1037]	28	! 1036 2012-10-22 13:43:42Z raasch
	29	! code put under GPL (PALM 3.9)
	30	!
[1017]	31	! 1015 2012-09-27 09:23:24Z raasch
	32	! accelerator version (*_acc) added
	33	!
[1011]	34	! 1010 2012-09-20 07:59:54Z raasch
	35	! cpp switch __nopointer added for pointer free version
	36	!
[1002]	37	! 1001 2012-09-13 14:08:46Z raasch
	38	! some arrays comunicated by module instead of parameter list
	39	!
[668]	40	! 667 2010-12-23 12:06:00Z suehring/gryschka
	41	! nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng
	42	!
[198]	43	! 183 2008-08-04 15:39:12Z letzel
	44	! bugfix: calculation of fluxes at vertical surfaces
	45	!
[139]	46	! 129 2007-10-30 12:12:24Z letzel
	47	! replace wall_heatflux by wall_s_flux that is now included in the parameter
	48	! list, bugfix for assignment of fluxes at walls
	49	!
[39]	50	! 20 2007-02-26 00:12:32Z raasch
	51	! Bugfix: ddzw dimensioned 1:nzt"+1"
	52	! Calculation extended for gridpoint nzt, fluxes can be given at top,
	53	! +s_flux_t in parameter list, s_flux renamed s_flux_b
	54	!
[3]	55	! RCS Log replace by Id keyword, revision history cleaned up
	56	!
[1]	57	! Revision 1.8 2006/02/23 10:34:17 raasch
	58	! nzb_2d replaced by nzb_s_outer in horizontal diffusion and by nzb_s_inner
	59	! or nzb_diff_s_inner, respectively, in vertical diffusion, prescribed surface
	60	! fluxes at vertically oriented topography
	61	!
	62	! Revision 1.1 2000/04/13 14:54:02 schroeter
	63	! Initial revision
	64	!
	65	!
	66	! Description:
	67	! ------------
	68	! Diffusion term of scalar quantities (temperature and water content)
	69	!------------------------------------------------------------------------------!
	70
	71	PRIVATE
[1015]	72	PUBLIC diffusion_s, diffusion_s_acc
[1]	73
	74	INTERFACE diffusion_s
	75	MODULE PROCEDURE diffusion_s
	76	MODULE PROCEDURE diffusion_s_ij
	77	END INTERFACE diffusion_s
	78
[1015]	79	INTERFACE diffusion_s_acc
	80	MODULE PROCEDURE diffusion_s_acc
	81	END INTERFACE diffusion_s_acc
	82
[1]	83	CONTAINS
	84
	85
	86	!------------------------------------------------------------------------------!
	87	! Call for all grid points
	88	!------------------------------------------------------------------------------!
[1001]	89	SUBROUTINE diffusion_s( s, s_flux_b, s_flux_t, wall_s_flux )
[1]	90
[1001]	91	USE arrays_3d
[1]	92	USE control_parameters
	93	USE grid_variables
	94	USE indices
	95
	96	IMPLICIT NONE
	97
	98	INTEGER :: i, j, k
	99	REAL :: vertical_gridspace
[129]	100	REAL :: wall_s_flux(0:4)
[1001]	101	REAL, DIMENSION(nysg:nyng,nxlg:nxrg) :: s_flux_b, s_flux_t
[1010]	102	#if defined( __nopointer )
	103	REAL, DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: s
	104	#else
[1001]	105	REAL, DIMENSION(:,:,:), POINTER :: s
[1010]	106	#endif
[1]	107
	108	DO i = nxl, nxr
	109	DO j = nys,nyn
	110	!
	111	!-- Compute horizontal diffusion
[19]	112	DO k = nzb_s_outer(j,i)+1, nzt
[1]	113
	114	tend(k,j,i) = tend(k,j,i) &
	115	+ 0.5 * ( &
	116	( kh(k,j,i) + kh(k,j,i+1) ) * ( s(k,j,i+1)-s(k,j,i) ) &
	117	- ( kh(k,j,i) + kh(k,j,i-1) ) * ( s(k,j,i)-s(k,j,i-1) ) &
	118	) * ddx2 &
	119	+ 0.5 * ( &
	120	( kh(k,j,i) + kh(k,j+1,i) ) * ( s(k,j+1,i)-s(k,j,i) ) &
	121	- ( kh(k,j,i) + kh(k,j-1,i) ) * ( s(k,j,i)-s(k,j-1,i) ) &
	122	) * ddy2
	123	ENDDO
	124
	125	!
	126	!-- Apply prescribed horizontal wall heatflux where necessary
	127	IF ( ( wall_w_x(j,i) .NE. 0.0 ) .OR. ( wall_w_y(j,i) .NE. 0.0 ) ) &
	128	THEN
	129	DO k = nzb_s_inner(j,i)+1, nzb_s_outer(j,i)
	130
	131	tend(k,j,i) = tend(k,j,i) &
[183]	132	+ ( fwxp(j,i) * 0.5 * &
[1]	133	( kh(k,j,i) + kh(k,j,i+1) ) * ( s(k,j,i+1)-s(k,j,i) ) &
[129]	134	+ ( 1.0 - fwxp(j,i) ) * wall_s_flux(1) &
[183]	135	-fwxm(j,i) * 0.5 * &
[1]	136	( kh(k,j,i) + kh(k,j,i-1) ) * ( s(k,j,i)-s(k,j,i-1) ) &
[129]	137	+ ( 1.0 - fwxm(j,i) ) * wall_s_flux(2) &
[1]	138	) * ddx2 &
[183]	139	+ ( fwyp(j,i) * 0.5 * &
[1]	140	( kh(k,j,i) + kh(k,j+1,i) ) * ( s(k,j+1,i)-s(k,j,i) ) &
[129]	141	+ ( 1.0 - fwyp(j,i) ) * wall_s_flux(3) &
[183]	142	-fwym(j,i) * 0.5 * &
[1]	143	( kh(k,j,i) + kh(k,j-1,i) ) * ( s(k,j,i)-s(k,j-1,i) ) &
[129]	144	+ ( 1.0 - fwym(j,i) ) * wall_s_flux(4) &
[1]	145	) * ddy2
	146	ENDDO
	147	ENDIF
	148
	149	!
	150	!-- Compute vertical diffusion. In case that surface fluxes have been
[19]	151	!-- prescribed or computed at bottom and/or top, index k starts/ends at
	152	!-- nzb+2 or nzt-1, respectively.
	153	DO k = nzb_diff_s_inner(j,i), nzt_diff
[1]	154
	155	tend(k,j,i) = tend(k,j,i) &
	156	+ 0.5 * ( &
	157	( kh(k,j,i) + kh(k+1,j,i) ) * ( s(k+1,j,i)-s(k,j,i) ) * ddzu(k+1) &
	158	- ( kh(k,j,i) + kh(k-1,j,i) ) * ( s(k,j,i)-s(k-1,j,i) ) * ddzu(k) &
	159	) * ddzw(k)
	160	ENDDO
	161
	162	!
[19]	163	!-- Vertical diffusion at the first computational gridpoint along
[1]	164	!-- z-direction
	165	IF ( use_surface_fluxes ) THEN
	166
	167	k = nzb_s_inner(j,i)+1
	168
	169	tend(k,j,i) = tend(k,j,i) &
	170	+ ( 0.5 * ( kh(k,j,i)+kh(k+1,j,i) ) &
	171	* ( s(k+1,j,i)-s(k,j,i) ) &
	172	* ddzu(k+1) &
[19]	173	+ s_flux_b(j,i) &
[1]	174	) * ddzw(k)
	175
	176	ENDIF
	177
[19]	178	!
	179	!-- Vertical diffusion at the last computational gridpoint along
	180	!-- z-direction
	181	IF ( use_top_fluxes ) THEN
	182
	183	k = nzt
	184
	185	tend(k,j,i) = tend(k,j,i) &
	186	+ ( - s_flux_t(j,i) &
[20]	187	- 0.5 * ( kh(k-1,j,i)+kh(k,j,i) ) &
	188	* ( s(k,j,i)-s(k-1,j,i) ) &
	189	* ddzu(k) &
[19]	190	) * ddzw(k)
	191
	192	ENDIF
	193
[1]	194	ENDDO
	195	ENDDO
	196
	197	END SUBROUTINE diffusion_s
	198
	199
	200	!------------------------------------------------------------------------------!
[1015]	201	! Call for all grid points - accelerator version
	202	!------------------------------------------------------------------------------!
	203	SUBROUTINE diffusion_s_acc( s, s_flux_b, s_flux_t, wall_s_flux )
	204
	205	USE arrays_3d
	206	USE control_parameters
	207	USE grid_variables
	208	USE indices
	209
	210	IMPLICIT NONE
	211
	212	INTEGER :: i, j, k
	213	REAL :: vertical_gridspace
	214	REAL :: wall_s_flux(0:4)
	215	REAL, DIMENSION(nysg:nyng,nxlg:nxrg) :: s_flux_b, s_flux_t
	216	#if defined( __nopointer )
	217	REAL, DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: s
	218	#else
	219	REAL, DIMENSION(:,:,:), POINTER :: s
	220	#endif
	221
	222	!$acc kernels present( ddzu, ddzw, fwxm, fwxp, fwym, fwyp, kh ) &
	223	!$acc present( nzb_diff_s_inner, nzb_s_inner, nzb_s_outer, s ) &
	224	!$acc present( s_flux_b, s_flux_t, tend, wall_s_flux ) &
	225	!$acc present( wall_w_x, wall_w_y )
	226	!$acc loop
	227	DO i = nxl, nxr
	228	DO j = nys,nyn
	229	!
	230	!-- Compute horizontal diffusion
	231	!$acc loop vector( 32 )
	232	DO k = 1, nzt
	233	IF ( k > nzb_s_outer(j,i) ) THEN
	234
	235	tend(k,j,i) = tend(k,j,i) &
	236	+ 0.5 * ( &
	237	( kh(k,j,i) + kh(k,j,i+1) ) * ( s(k,j,i+1)-s(k,j,i) ) &
	238	- ( kh(k,j,i) + kh(k,j,i-1) ) * ( s(k,j,i)-s(k,j,i-1) ) &
	239	) * ddx2 &
	240	+ 0.5 * ( &
	241	( kh(k,j,i) + kh(k,j+1,i) ) * ( s(k,j+1,i)-s(k,j,i) ) &
	242	- ( kh(k,j,i) + kh(k,j-1,i) ) * ( s(k,j,i)-s(k,j-1,i) ) &
	243	) * ddy2
	244	ENDIF
	245	ENDDO
	246
	247	!
	248	!-- Apply prescribed horizontal wall heatflux where necessary
	249	!$acc loop vector(32)
	250	DO k = 1, nzt
	251	IF ( k > nzb_s_inner(j,i) .AND. k <= nzb_s_outer(j,i) .AND. &
	252	( wall_w_x(j,i) /= 0.0 .OR. wall_w_y(j,i) /= 0.0 ) ) &
	253	THEN
	254	tend(k,j,i) = tend(k,j,i) &
	255	+ ( fwxp(j,i) * 0.5 * &
	256	( kh(k,j,i) + kh(k,j,i+1) ) * ( s(k,j,i+1)-s(k,j,i) ) &
	257	+ ( 1.0 - fwxp(j,i) ) * wall_s_flux(1) &
	258	-fwxm(j,i) * 0.5 * &
	259	( kh(k,j,i) + kh(k,j,i-1) ) * ( s(k,j,i)-s(k,j,i-1) ) &
	260	+ ( 1.0 - fwxm(j,i) ) * wall_s_flux(2) &
	261	) * ddx2 &
	262	+ ( fwyp(j,i) * 0.5 * &
	263	( kh(k,j,i) + kh(k,j+1,i) ) * ( s(k,j+1,i)-s(k,j,i) ) &
	264	+ ( 1.0 - fwyp(j,i) ) * wall_s_flux(3) &
	265	-fwym(j,i) * 0.5 * &
	266	( kh(k,j,i) + kh(k,j-1,i) ) * ( s(k,j,i)-s(k,j-1,i) ) &
	267	+ ( 1.0 - fwym(j,i) ) * wall_s_flux(4) &
	268	) * ddy2
	269	ENDIF
	270	ENDDO
	271
	272	!
	273	!-- Compute vertical diffusion. In case that surface fluxes have been
	274	!-- prescribed or computed at bottom and/or top, index k starts/ends at
	275	!-- nzb+2 or nzt-1, respectively.
	276	!$acc loop vector( 32 )
	277	DO k = 1, nzt_diff
	278	IF ( k >= nzb_diff_s_inner(j,i) ) THEN
	279	tend(k,j,i) = tend(k,j,i) &
	280	+ 0.5 * ( &
	281	( kh(k,j,i) + kh(k+1,j,i) ) * ( s(k+1,j,i)-s(k,j,i) ) * ddzu(k+1) &
	282	- ( kh(k,j,i) + kh(k-1,j,i) ) * ( s(k,j,i)-s(k-1,j,i) ) * ddzu(k) &
	283	) * ddzw(k)
	284	ENDIF
	285	ENDDO
	286
	287	!
	288	!-- Vertical diffusion at the first computational gridpoint along
	289	!-- z-direction
	290	!$acc loop vector( 32 )
	291	DO k = 1, nzt
	292	IF ( use_surface_fluxes .AND. k == nzb_s_inner(j,i)+1 ) THEN
	293	tend(k,j,i) = tend(k,j,i) &
	294	+ ( 0.5 * ( kh(k,j,i)+kh(k+1,j,i) ) &
	295	* ( s(k+1,j,i)-s(k,j,i) ) &
	296	* ddzu(k+1) &
	297	+ s_flux_b(j,i) &
	298	) * ddzw(k)
	299	ENDIF
	300
	301	!
	302	!-- Vertical diffusion at the last computational gridpoint along
	303	!-- z-direction
	304	IF ( use_top_fluxes .AND. k == nzt ) THEN
	305	tend(k,j,i) = tend(k,j,i) &
	306	+ ( - s_flux_t(j,i) &
	307	- 0.5 * ( kh(k-1,j,i)+kh(k,j,i) )&
	308	* ( s(k,j,i)-s(k-1,j,i) ) &
	309	* ddzu(k) &
	310	) * ddzw(k)
	311	ENDIF
	312	ENDDO
	313
	314	ENDDO
	315	ENDDO
	316	!$acc end kernels
	317
	318	END SUBROUTINE diffusion_s_acc
	319
	320
	321	!------------------------------------------------------------------------------!
[1]	322	! Call for grid point i,j
	323	!------------------------------------------------------------------------------!
[1001]	324	SUBROUTINE diffusion_s_ij( i, j, s, s_flux_b, s_flux_t, wall_s_flux )
[1]	325
[1001]	326	USE arrays_3d
[1]	327	USE control_parameters
	328	USE grid_variables
	329	USE indices
	330
	331	IMPLICIT NONE
	332
	333	INTEGER :: i, j, k
	334	REAL :: vertical_gridspace
[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 |