Home

Context Navigation

source: palm/trunk/SOURCE/buoyancy.f90 @ 1034

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

Rev	Line
[1]	1	MODULE buoyancy_mod
	2
	3	!------------------------------------------------------------------------------!
[247]	4	! Currrent revisions:
[1]	5	! -----------------
[98]	6	!
[1017]	7	!
[98]	8	! Former revisions:
	9	! -----------------
	10	! $Id: buoyancy.f90 1017 2012-09-27 11:28:50Z raasch $
	11	!
[1017]	12	! 1015 2012-09-27 09:23:24Z raasch
	13	! accelerator version (*_acc) added
	14	!
[1011]	15	! 1010 2012-09-20 07:59:54Z raasch
	16	! cpp switch __nopointer added for pointer free version
	17	!
[623]	18	! 622 2010-12-10 08:08:13Z raasch
	19	! optional barriers included in order to speed up collective operations
	20	!
[516]	21	! 515 2010-03-18 02:30:38Z raasch
	22	! PGI-compiler creates SIGFPE in routine buoyancy, if opt>1 is used! Therefore,
	23	! opt=1 is forced by compiler-directive.
	24	!
[392]	25	! 247 2009-02-27 14:01:30Z heinze
	26	! Output of messages replaced by message handling routine
	27	!
[139]	28	! 132 2007-11-20 09:46:11Z letzel
	29	! Vertical scalar profiles now based on nzb_s_inner and ngp_2dh_s_inner.
	30	!
[110]	31	! 106 2007-08-16 14:30:26Z raasch
	32	! i loop for u-component (sloping surface) is starting from nxlu (needed for
	33	! non-cyclic boundary conditions)
	34	!
[98]	35	! 97 2007-06-21 08:23:15Z raasch
[97]	36	! Routine reneralized to be used with temperature AND density:
	37	! argument theta renamed var, new argument var_reference,
	38	! use_pt_reference renamed use_reference,
[96]	39	! calc_mean_pt_profile renamed calc_mean_profile
[1]	40	!
[77]	41	! 57 2007-03-09 12:05:41Z raasch
	42	! Reference temperature pt_reference can be used.
	43	!
[3]	44	! RCS Log replace by Id keyword, revision history cleaned up
	45	!
[1]	46	! Revision 1.19 2006/04/26 12:09:56 raasch
	47	! OpenMP optimization (one dimension added to sums_l)
	48	!
	49	! Revision 1.1 1997/08/29 08:56:48 raasch
	50	! Initial revision
	51	!
	52	!
	53	! Description:
	54	! ------------
	55	! Buoyancy term of the third component of the equation of motion.
	56	! WARNING: humidity is not regarded when using a sloping surface!
	57	!------------------------------------------------------------------------------!
	58
	59	PRIVATE
[1015]	60	PUBLIC buoyancy, buoyancy_acc, calc_mean_profile
[1]	61
	62	INTERFACE buoyancy
	63	MODULE PROCEDURE buoyancy
	64	MODULE PROCEDURE buoyancy_ij
	65	END INTERFACE buoyancy
	66
[1015]	67	INTERFACE buoyancy_acc
	68	MODULE PROCEDURE buoyancy_acc
	69	END INTERFACE buoyancy_acc
	70
[96]	71	INTERFACE calc_mean_profile
	72	MODULE PROCEDURE calc_mean_profile
	73	END INTERFACE calc_mean_profile
[1]	74
	75	CONTAINS
	76
	77
	78	!------------------------------------------------------------------------------!
	79	! Call for all grid points
	80	!------------------------------------------------------------------------------!
[97]	81	SUBROUTINE buoyancy( var, var_reference, wind_component, pr )
[1]	82
	83	USE arrays_3d
	84	USE control_parameters
	85	USE indices
	86	USE pegrid
	87	USE statistics
	88
	89	IMPLICIT NONE
	90
	91	INTEGER :: i, j, k, pr, wind_component
[97]	92	REAL :: var_reference
[1010]	93	#if defined( __nopointer )
	94	REAL, DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: var
	95	#else
	96	REAL, DIMENSION(:,:,:), POINTER :: var
	97	#endif
[1]	98
	99
	100	IF ( .NOT. sloping_surface ) THEN
	101	!
	102	!-- Normal case: horizontal surface
[97]	103	IF ( use_reference ) THEN
[57]	104	DO i = nxl, nxr
	105	DO j = nys, nyn
	106	DO k = nzb_s_inner(j,i)+1, nzt-1
[97]	107	tend(k,j,i) = tend(k,j,i) + atmos_ocean_sign * g * 0.5 * &
	108	( &
	109	( var(k,j,i) - hom(k,1,pr,0) ) / var_reference + &
	110	( var(k+1,j,i) - hom(k+1,1,pr,0) ) / var_reference &
[57]	111	)
	112	ENDDO
	113	ENDDO
	114	ENDDO
	115	ELSE
	116	DO i = nxl, nxr
	117	DO j = nys, nyn
	118	DO k = nzb_s_inner(j,i)+1, nzt-1
[97]	119	tend(k,j,i) = tend(k,j,i) + atmos_ocean_sign * g * 0.5 * &
	120	( &
	121	( var(k,j,i) - hom(k,1,pr,0) ) / hom(k,1,pr,0) + &
	122	( var(k+1,j,i) - hom(k+1,1,pr,0) ) / hom(k+1,1,pr,0) &
[57]	123	)
	124	ENDDO
[1]	125	ENDDO
	126	ENDDO
[57]	127	ENDIF
[1]	128
	129	ELSE
	130	!
	131	!-- Buoyancy term for a surface with a slope in x-direction. The equations
	132	!-- for both the u and w velocity-component contain proportionate terms.
	133	!-- Temperature field at time t=0 serves as environmental temperature.
	134	!-- Reference temperature (pt_surface) is the one at the lower left corner
	135	!-- of the total domain.
	136	IF ( wind_component == 1 ) THEN
	137
[106]	138	DO i = nxlu, nxr
[1]	139	DO j = nys, nyn
	140	DO k = nzb_s_inner(j,i)+1, nzt-1
	141	tend(k,j,i) = tend(k,j,i) + g * sin_alpha_surface * &
	142	0.5 * ( ( pt(k,j,i-1) + pt(k,j,i) ) &
	143	- ( pt_slope_ref(k,i-1) + pt_slope_ref(k,i) ) &
	144	) / pt_surface
	145	ENDDO
	146	ENDDO
	147	ENDDO
	148
	149	ELSEIF ( wind_component == 3 ) THEN
	150
	151	DO i = nxl, nxr
	152	DO j = nys, nyn
	153	DO k = nzb_s_inner(j,i)+1, nzt-1
	154	tend(k,j,i) = tend(k,j,i) + g * cos_alpha_surface * &
	155	0.5 * ( ( pt(k,j,i) + pt(k+1,j,i) ) &
	156	- ( pt_slope_ref(k,i) + pt_slope_ref(k+1,i) ) &
	157	) / pt_surface
	158	ENDDO
	159	ENDDO
	160	ENDDO
	161
	162	ELSE
[247]	163
	164	WRITE( message_string, * ) 'no term for component "',&
[1]	165	wind_component,'"'
[247]	166	CALL message( 'buoyancy', 'PA0159', 1, 2, 0, 6, 0 )
[1]	167
	168	ENDIF
	169
	170	ENDIF
	171
	172	END SUBROUTINE buoyancy
	173
	174
	175	!------------------------------------------------------------------------------!
[1015]	176	! Call for all grid points - accelerator version
	177	!------------------------------------------------------------------------------!
	178	SUBROUTINE buoyancy_acc( var, var_reference, wind_component, pr )
	179
	180	USE arrays_3d
	181	USE control_parameters
	182	USE indices
	183	USE pegrid
	184	USE statistics
	185
	186	IMPLICIT NONE
	187
	188	INTEGER :: i, j, k, pr, wind_component
	189	REAL :: var_reference
	190	#if defined( __nopointer )
	191	REAL, DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: var
	192	#else
	193	REAL, DIMENSION(:,:,:), POINTER :: var
	194	#endif
	195
	196
	197	IF ( .NOT. sloping_surface ) THEN
	198	!
	199	!-- Normal case: horizontal surface
	200	IF ( use_reference ) THEN
	201	DO i = nxl, nxr
	202	DO j = nys, nyn
	203	DO k = nzb_s_inner(j,i)+1, nzt-1
	204	tend(k,j,i) = tend(k,j,i) + atmos_ocean_sign * g * 0.5 * &
	205	( &
	206	( var(k,j,i) - hom(k,1,pr,0) ) / var_reference + &
	207	( var(k+1,j,i) - hom(k+1,1,pr,0) ) / var_reference &
	208	)
	209	ENDDO
	210	ENDDO
	211	ENDDO
	212	ELSE
	213	!$acc kernels present( nzb_s_inner, hom, tend, var )
	214	!$acc loop
	215	DO i = nxl, nxr
	216	DO j = nys, nyn
	217	!$acc loop vector(32)
	218	DO k = 1, nzt-1
	219	IF ( k > nzb_s_inner(j,i) ) THEN
	220	tend(k,j,i) = tend(k,j,i) + atmos_ocean_sign * g * 0.5 * &
	221	( &
	222	( var(k,j,i) - hom(k,1,pr,0) ) / hom(k,1,pr,0) + &
	223	( var(k+1,j,i) - hom(k+1,1,pr,0) ) / hom(k+1,1,pr,0) &
	224	)
	225	ENDIF
	226	ENDDO
	227	ENDDO
	228	ENDDO
	229	!$acc end kernels
	230	ENDIF
	231
	232	ELSE
	233	!
	234	!-- Buoyancy term for a surface with a slope in x-direction. The equations
	235	!-- for both the u and w velocity-component contain proportionate terms.
	236	!-- Temperature field at time t=0 serves as environmental temperature.
	237	!-- Reference temperature (pt_surface) is the one at the lower left corner
	238	!-- of the total domain.
	239	IF ( wind_component == 1 ) THEN
	240
	241	DO i = nxlu, nxr
	242	DO j = nys, nyn
	243	DO k = nzb_s_inner(j,i)+1, nzt-1
	244	tend(k,j,i) = tend(k,j,i) + g * sin_alpha_surface * &
	245	0.5 * ( ( pt(k,j,i-1) + pt(k,j,i) ) &
	246	- ( pt_slope_ref(k,i-1) + pt_slope_ref(k,i) ) &
	247	) / pt_surface
	248	ENDDO
	249	ENDDO
	250	ENDDO
	251
	252	ELSEIF ( wind_component == 3 ) THEN
	253
	254	DO i = nxl, nxr
	255	DO j = nys, nyn
	256	DO k = nzb_s_inner(j,i)+1, nzt-1
	257	tend(k,j,i) = tend(k,j,i) + g * cos_alpha_surface * &
	258	0.5 * ( ( pt(k,j,i) + pt(k+1,j,i) ) &
	259	- ( pt_slope_ref(k,i) + pt_slope_ref(k+1,i) ) &
	260	) / pt_surface
	261	ENDDO
	262	ENDDO
	263	ENDDO
	264
	265	ELSE
	266
	267	WRITE( message_string, * ) 'no term for component "',&
	268	wind_component,'"'
	269	CALL message( 'buoyancy', 'PA0159', 1, 2, 0, 6, 0 )
	270
	271	ENDIF
	272
	273	ENDIF
	274
	275	END SUBROUTINE buoyancy_acc
	276
	277
	278	!------------------------------------------------------------------------------!
[1]	279	! Call for grid point i,j
[515]	280	! ATTENTION: PGI-compiler creates SIGFPE if opt>1 is used! Therefore, opt=1 is
	281	! forced by compiler-directive.
[1]	282	!------------------------------------------------------------------------------!
[515]	283	!pgi$r opt=1
[97]	284	SUBROUTINE buoyancy_ij( i, j, var, var_reference, wind_component, pr )
[1]	285
	286	USE arrays_3d
	287	USE control_parameters
	288	USE indices
	289	USE pegrid
	290	USE statistics
	291
	292	IMPLICIT NONE
	293
	294	INTEGER :: i, j, k, pr, wind_component
[97]	295	REAL :: var_reference
[1010]	296	#if defined( __nopointer )
	297	REAL, DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: var
	298	#else
	299	REAL, DIMENSION(:,:,:), POINTER :: var
	300	#endif
[1]	301
	302
	303	IF ( .NOT. sloping_surface ) THEN
	304	!
	305	!-- Normal case: horizontal surface
[97]	306	IF ( use_reference ) THEN
[57]	307	DO k = nzb_s_inner(j,i)+1, nzt-1
[97]	308	tend(k,j,i) = tend(k,j,i) + atmos_ocean_sign * g * 0.5 * ( &
	309	( var(k,j,i) - hom(k,1,pr,0) ) / var_reference + &
	310	( var(k+1,j,i) - hom(k+1,1,pr,0) ) / var_reference &
	311	)
[57]	312	ENDDO
	313	ELSE
	314	DO k = nzb_s_inner(j,i)+1, nzt-1
[97]	315	tend(k,j,i) = tend(k,j,i) + atmos_ocean_sign * g * 0.5 * ( &
	316	( var(k,j,i) - hom(k,1,pr,0) ) / hom(k,1,pr,0) + &
	317	( var(k+1,j,i) - hom(k+1,1,pr,0) ) / hom(k+1,1,pr,0) &
	318	)
[57]	319	ENDDO
	320	ENDIF
[1]	321
	322	ELSE
	323	!
	324	!-- Buoyancy term for a surface with a slope in x-direction. The equations
	325	!-- for both the u and w velocity-component contain proportionate terms.
	326	!-- Temperature field at time t=0 serves as environmental temperature.
	327	!-- Reference temperature (pt_surface) is the one at the lower left corner
	328	!-- of the total domain.
	329	IF ( wind_component == 1 ) THEN
	330
	331	DO k = nzb_s_inner(j,i)+1, nzt-1
	332	tend(k,j,i) = tend(k,j,i) + g * sin_alpha_surface * &
	333	0.5 * ( ( pt(k,j,i-1) + pt(k,j,i) ) &
	334	- ( pt_slope_ref(k,i-1) + pt_slope_ref(k,i) ) &
	335	) / pt_surface
	336	ENDDO
	337
	338	ELSEIF ( wind_component == 3 ) THEN
	339
	340	DO k = nzb_s_inner(j,i)+1, nzt-1
	341	tend(k,j,i) = tend(k,j,i) + g * cos_alpha_surface * &
	342	0.5 * ( ( pt(k,j,i) + pt(k+1,j,i) ) &
	343	- ( pt_slope_ref(k,i) + pt_slope_ref(k+1,i) ) &
	344	) / pt_surface
	345	ENDDO
	346
	347	ELSE
	348
[247]	349	WRITE( message_string, * ) 'no term for component "',&
[1]	350	wind_component,'"'
[247]	351	CALL message( 'buoyancy', 'PA0159', 1, 2, 0, 6, 0 )
[1]	352
	353	ENDIF
	354
	355	ENDIF
	356
	357	END SUBROUTINE buoyancy_ij
	358
	359
[96]	360	SUBROUTINE calc_mean_profile( var, pr )
[1]	361
	362	!------------------------------------------------------------------------------!
	363	! Description:
	364	! ------------
	365	! Calculate the horizontally averaged vertical temperature profile (pr=4 in case
	366	! of potential temperature and 44 in case of virtual potential temperature).
	367	!------------------------------------------------------------------------------!
	368
	369	USE control_parameters
	370	USE indices
	371	USE pegrid
	372	USE statistics
	373
	374	IMPLICIT NONE
	375
	376	INTEGER :: i, j, k, omp_get_thread_num, pr, tn
[1010]	377	#if defined( __nopointer )
	378	REAL, DIMENSION(:,:,:) :: var
	379	#else
	380	REAL, DIMENSION(:,:,:), POINTER :: var
	381	#endif
[1]	382
	383	!
[96]	384	!-- Computation of the horizontally averaged profile of variable var, unless
[1]	385	!-- already done by the relevant call from flow_statistics. The calculation
	386	!-- is done only for the first respective intermediate timestep in order to
	387	!-- spare communication time and to produce identical model results with jobs
	388	!-- which are calling flow_statistics at different time intervals.
	389	IF ( .NOT. flow_statistics_called .AND. &
	390	intermediate_timestep_count == 1 ) THEN
	391
	392	!
[96]	393	!-- Horizontal average of variable var
[1]	394	tn = 0 ! Default thread number in case of one thread
	395	!$OMP PARALLEL PRIVATE( i, j, k, tn )
	396	!$ tn = omp_get_thread_num()
	397	sums_l(:,pr,tn) = 0.0
	398	!$OMP DO
	399	DO i = nxl, nxr
	400	DO j = nys, nyn
[132]	401	DO k = nzb_s_inner(j,i), nzt+1
[96]	402	sums_l(k,pr,tn) = sums_l(k,pr,tn) + var(k,j,i)
[1]	403	ENDDO
	404	ENDDO
	405	ENDDO
	406	!$OMP END PARALLEL
	407
	408	DO i = 1, threads_per_task-1
	409	sums_l(:,pr,0) = sums_l(:,pr,0) + sums_l(:,pr,i)
	410	ENDDO
	411
	412	#if defined( __parallel )
	413
[622]	414	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
[1]	415	CALL MPI_ALLREDUCE( sums_l(nzb,pr,0), sums(nzb,pr), nzt+2-nzb, &
	416	MPI_REAL, MPI_SUM, comm2d, ierr )
	417
	418	#else
	419
	420	sums(:,pr) = sums_l(:,pr,0)
	421
	422	#endif
	423
[132]	424	hom(:,1,pr,0) = sums(:,pr) / ngp_2dh_s_inner(:,0)
[1]	425
	426	ENDIF
	427
[96]	428	END SUBROUTINE calc_mean_profile
[1]	429
	430	END MODULE buoyancy_mod

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |