Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |