Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |