Home

Context Navigation

buoyancy.f90 @ 80

Last change on this file since 80 was 77, checked in by raasch, 18 years ago

New:
---

particle reflection from vertical walls implemented, particle SGS model adjusted to walls

Wall functions for vertical walls now include diabatic conditions. New subroutines wall_fluxes, wall_fluxes_e. New 4D-array rif_wall.

new d3par-parameter netcdf_64bit_3d to switch on 64bit offset only for 3D files

new d3par-parameter dt_max to define the maximum value for the allowed timestep

new inipar-parameter loop_optimization to control the loop optimization method

new inipar-parameter pt_refrence. If given, this value is used as the reference that in buoyancy terms (otherwise, the instantaneous horizontally averaged temperature is used).

new user interface user_advec_particles

new initializing action "by_user" calls user_init_3d_model and allows the initial setting of all 3d arrays

topography height informations are stored on arrays zu_s_inner and zw_w_inner and output to the 2d/3d NetCDF files

samples added to the user interface which show how to add user-define time series quantities.

calculation/output of precipitation amount, precipitation rate and z0 (by setting "pra*", "prr*", "z0*" with data_output). The time interval on which the precipitation amount is defined is set by new d3par-parameter precipitation_amount_interval

unit 9 opened for debug output (file DEBUG_<pe#>)

Makefile, advec_particles, average_3d_data, buoyancy, calc_precipitation, check_open, check_parameters, data_output_2d, diffusion_e, diffusion_u, diffusion_v, diffusion_w, diffusivities, header, impact_of_latent_heat, init_particles, init_3d_model, modules, netcdf, parin, production_e, read_var_list, read_3d_binary, sum_up_3d_data, user_interface, write_var_list, write_3d_binary

New: wall_fluxes

Changed:

General revision of non-cyclic horizontal boundary conditions: radiation boundary conditions are now used instead of Neumann conditions at the outflow (calculation needs velocity values for t-dt, which are stored on new arrays u_m_l, u_m_r, etc.), calculation of mean outflow is not needed any more, volume flow control is added for the outflow boundary (currently only for the north boundary!!), additional gridpoints along x and y (uxrp, vynp) are not needed any more, routine "boundary_conds" now operates on timelevel t+dt and is not split in two parts (main, uvw_outflow) any more, Neumann boundary conditions at inflow/outflow in case of non-cyclic boundary conditions for all 2d-arrays that are handled by exchange_horiz_2d

The FFT-method for solving the Poisson-equation is now working with Neumann boundary conditions both at the bottom and the top. This requires adjustments of the tridiagonal coefficients and subtracting the horizontally averaged mean from the vertical velocity field.

+age_m in particle_type

Particles-package is now part of the default code ("-p particles" is not needed any more).

Move call of user_actions( 'after_integration' ) below increment of times
and counters. user_actions is now called for each statistic region and has as an argument the number of the respective region (sr)

d3par-parameter data_output_ts removed. Timeseries output for "profil" removed. Timeseries are now switched on by dt_dots. Timeseries data is collected in flow_statistics.

Initial velocities at nzb+1 are regarded for volume flow control in case they have been set zero before (to avoid small timesteps); see new internal parameters u/v_nzb_p1_for_vfc.

q is not allowed to become negative (prognostic_equations).

poisfft_init is only called if fft-solver is switched on (init_pegrid).

d3par-parameter moisture renamed to humidity.

Subversion global revision number is read from mrun and added to the run description header and to the run control (_rc) file.

vtk directives removed from main program.

The uitility routine interpret_config reads PALM environment variables from NAMELIST instead using the system call GETENV.

advec_u_pw, advec_u_up, advec_v_pw, advec_v_up, asselin_filter, check_parameters, coriolis, data_output_dvrp, data_output_ptseries, data_output_ts, data_output_2d, data_output_3d, diffusion_u, diffusion_v, exchange_horiz, exchange_horiz_2d, flow_statistics, header, init_grid, init_particles, init_pegrid, init_rankine, init_pt_anomaly, init_1d_model, init_3d_model, modules, palm, package_parin, parin, poisfft, poismg, prandtl_fluxes, pres, production_e, prognostic_equations, read_var_list, read_3d_binary, sor, swap_timelevel, time_integration, write_var_list, write_3d_binary

Errors:

Bugfix: preset of tendencies te_em, te_um, te_vm in init_1d_model

Bugfix in sample for reading user defined data from restart file (user_init)

Bugfix in setting diffusivities for cases with the outflow damping layer extending over more than one subdomain (init_3d_model)

Check for possible negative humidities in the initial humidity profile.

in Makefile, default suffixes removed from the suffix list to avoid calling of m2c in
# case of .mod files

Makefile
check_parameters, init_1d_model, init_3d_model, user_interface

Property svn:keywords set to Id

File size: 9.1 KB

Rev	Line
[1]	1	MODULE buoyancy_mod
	2
	3	!------------------------------------------------------------------------------!
	4	! Actual revisions:
	5	! -----------------
[77]	6	!
[1]	7	!
	8	! Former revisions:
	9	! -----------------
[3]	10	! $Id: buoyancy.f90 77 2007-03-29 04:26:56Z raasch $
[77]	11	!
	12	! 57 2007-03-09 12:05:41Z raasch
	13	! Reference temperature pt_reference can be used.
	14	!
[3]	15	! RCS Log replace by Id keyword, revision history cleaned up
	16	!
[1]	17	! Revision 1.19 2006/04/26 12:09:56 raasch
	18	! OpenMP optimization (one dimension added to sums_l)
	19	!
	20	! Revision 1.1 1997/08/29 08:56:48 raasch
	21	! Initial revision
	22	!
	23	!
	24	! Description:
	25	! ------------
	26	! Buoyancy term of the third component of the equation of motion.
	27	! WARNING: humidity is not regarded when using a sloping surface!
	28	!------------------------------------------------------------------------------!
	29
	30	PRIVATE
	31	PUBLIC buoyancy, calc_mean_pt_profile
	32
	33	INTERFACE buoyancy
	34	MODULE PROCEDURE buoyancy
	35	MODULE PROCEDURE buoyancy_ij
	36	END INTERFACE buoyancy
	37
	38	INTERFACE calc_mean_pt_profile
	39	MODULE PROCEDURE calc_mean_pt_profile
	40	END INTERFACE calc_mean_pt_profile
	41
	42	CONTAINS
	43
	44
	45	!------------------------------------------------------------------------------!
	46	! Call for all grid points
	47	!------------------------------------------------------------------------------!
	48	SUBROUTINE buoyancy( theta, wind_component, pr )
	49
	50	USE arrays_3d
	51	USE control_parameters
	52	USE indices
	53	USE pegrid
	54	USE statistics
	55
	56	IMPLICIT NONE
	57
	58	INTEGER :: i, j, k, pr, wind_component
	59	REAL, DIMENSION(:,:,:), POINTER :: theta
	60
	61
	62	IF ( .NOT. sloping_surface ) THEN
	63	!
	64	!-- Normal case: horizontal surface
[57]	65	IF ( use_pt_reference ) THEN
	66	DO i = nxl, nxr
	67	DO j = nys, nyn
	68	DO k = nzb_s_inner(j,i)+1, nzt-1
	69	tend(k,j,i) = tend(k,j,i) + g * 0.5 * ( &
	70	( theta(k,j,i) - hom(k,1,pr,0) ) / pt_reference + &
	71	( theta(k+1,j,i) - hom(k+1,1,pr,0) ) / pt_reference &
	72	)
	73	ENDDO
	74	ENDDO
	75	ENDDO
	76	ELSE
	77	DO i = nxl, nxr
	78	DO j = nys, nyn
	79	DO k = nzb_s_inner(j,i)+1, nzt-1
	80	tend(k,j,i) = tend(k,j,i) + g * 0.5 * ( &
[1]	81	( theta(k,j,i) - hom(k,1,pr,0) ) / hom(k,1,pr,0) + &
	82	( theta(k+1,j,i) - hom(k+1,1,pr,0) ) / hom(k+1,1,pr,0) &
[57]	83	)
	84	ENDDO
[1]	85	ENDDO
	86	ENDDO
[57]	87	ENDIF
[1]	88
	89	ELSE
	90	!
	91	!-- Buoyancy term for a surface with a slope in x-direction. The equations
	92	!-- for both the u and w velocity-component contain proportionate terms.
	93	!-- Temperature field at time t=0 serves as environmental temperature.
	94	!-- Reference temperature (pt_surface) is the one at the lower left corner
	95	!-- of the total domain.
	96	IF ( wind_component == 1 ) THEN
	97
	98	DO i = nxl, nxr
	99	DO j = nys, nyn
	100	DO k = nzb_s_inner(j,i)+1, nzt-1
	101	tend(k,j,i) = tend(k,j,i) + g * sin_alpha_surface * &
	102	0.5 * ( ( pt(k,j,i-1) + pt(k,j,i) ) &
	103	- ( pt_slope_ref(k,i-1) + pt_slope_ref(k,i) ) &
	104	) / pt_surface
	105	ENDDO
	106	ENDDO
	107	ENDDO
	108
	109	ELSEIF ( wind_component == 3 ) THEN
	110
	111	DO i = nxl, nxr
	112	DO j = nys, nyn
	113	DO k = nzb_s_inner(j,i)+1, nzt-1
	114	tend(k,j,i) = tend(k,j,i) + g * cos_alpha_surface * &
	115	0.5 * ( ( pt(k,j,i) + pt(k+1,j,i) ) &
	116	- ( pt_slope_ref(k,i) + pt_slope_ref(k+1,i) ) &
	117	) / pt_surface
	118	ENDDO
	119	ENDDO
	120	ENDDO
	121
	122	ELSE
	123
	124	IF ( myid == 0 ) PRINT*, '+++ buoyancy: no term for component "',&
	125	wind_component,'"'
	126	CALL local_stop
	127
	128	ENDIF
	129
	130	ENDIF
	131
	132	END SUBROUTINE buoyancy
	133
	134
	135	!------------------------------------------------------------------------------!
	136	! Call for grid point i,j
	137	!------------------------------------------------------------------------------!
	138	SUBROUTINE buoyancy_ij( i, j, theta, wind_component, pr )
	139
	140	USE arrays_3d
	141	USE control_parameters
	142	USE indices
	143	USE pegrid
	144	USE statistics
	145
	146	IMPLICIT NONE
	147
	148	INTEGER :: i, j, k, pr, wind_component
	149	REAL, DIMENSION(:,:,:), POINTER :: theta
	150
	151
	152	IF ( .NOT. sloping_surface ) THEN
	153	!
	154	!-- Normal case: horizontal surface
[57]	155	IF ( use_pt_reference ) THEN
	156	DO k = nzb_s_inner(j,i)+1, nzt-1
	157	tend(k,j,i) = tend(k,j,i) + g * 0.5 * ( &
	158	( theta(k,j,i) - hom(k,1,pr,0) ) / pt_reference + &
	159	( theta(k+1,j,i) - hom(k+1,1,pr,0) ) / pt_reference &
	160	)
	161	ENDDO
	162	ELSE
	163	DO k = nzb_s_inner(j,i)+1, nzt-1
	164	tend(k,j,i) = tend(k,j,i) + g * 0.5 * ( &
[1]	165	( theta(k,j,i) - hom(k,1,pr,0) ) / hom(k,1,pr,0) + &
	166	( theta(k+1,j,i) - hom(k+1,1,pr,0) ) / hom(k+1,1,pr,0) &
[57]	167	)
	168	ENDDO
	169	ENDIF
[1]	170
	171	ELSE
	172	!
	173	!-- Buoyancy term for a surface with a slope in x-direction. The equations
	174	!-- for both the u and w velocity-component contain proportionate terms.
	175	!-- Temperature field at time t=0 serves as environmental temperature.
	176	!-- Reference temperature (pt_surface) is the one at the lower left corner
	177	!-- of the total domain.
	178	IF ( wind_component == 1 ) THEN
	179
	180	DO k = nzb_s_inner(j,i)+1, nzt-1
	181	tend(k,j,i) = tend(k,j,i) + g * sin_alpha_surface * &
	182	0.5 * ( ( pt(k,j,i-1) + pt(k,j,i) ) &
	183	- ( pt_slope_ref(k,i-1) + pt_slope_ref(k,i) ) &
	184	) / pt_surface
	185	ENDDO
	186
	187	ELSEIF ( wind_component == 3 ) THEN
	188
	189	DO k = nzb_s_inner(j,i)+1, nzt-1
	190	tend(k,j,i) = tend(k,j,i) + g * cos_alpha_surface * &
	191	0.5 * ( ( pt(k,j,i) + pt(k+1,j,i) ) &
	192	- ( pt_slope_ref(k,i) + pt_slope_ref(k+1,i) ) &
	193	) / pt_surface
	194	ENDDO
	195
	196	ELSE
	197
	198	IF ( myid == 0 ) PRINT*, '+++ buoyancy: no term for component "',&
	199	wind_component,'"'
	200	CALL local_stop
	201
	202	ENDIF
	203
	204	ENDIF
	205
	206	END SUBROUTINE buoyancy_ij
	207
	208
	209	SUBROUTINE calc_mean_pt_profile( theta, pr )
	210
	211	!------------------------------------------------------------------------------!
	212	! Description:
	213	! ------------
	214	! Calculate the horizontally averaged vertical temperature profile (pr=4 in case
	215	! of potential temperature and 44 in case of virtual potential temperature).
	216	!------------------------------------------------------------------------------!
	217
	218	USE control_parameters
	219	USE indices
	220	USE pegrid
	221	USE statistics
	222
	223	IMPLICIT NONE
	224
	225	INTEGER :: i, j, k, omp_get_thread_num, pr, tn
	226	REAL, DIMENSION(:,:,:), POINTER :: theta
	227
	228	!
	229	!-- Computation of the horizontally averaged temperature profile, unless
	230	!-- already done by the relevant call from flow_statistics. The calculation
	231	!-- is done only for the first respective intermediate timestep in order to
	232	!-- spare communication time and to produce identical model results with jobs
	233	!-- which are calling flow_statistics at different time intervals.
	234	!-- Although this calculation is not required for model runs with a slope,
	235	!-- it is nevertheless also computed.
	236	IF ( .NOT. flow_statistics_called .AND. &
	237	intermediate_timestep_count == 1 ) THEN
	238
	239	!
	240	!-- Horizontal average of the potential temperature
	241	tn = 0 ! Default thread number in case of one thread
	242	!$OMP PARALLEL PRIVATE( i, j, k, tn )
	243	!$ tn = omp_get_thread_num()
	244	sums_l(:,pr,tn) = 0.0
	245	!$OMP DO
	246	DO i = nxl, nxr
	247	DO j = nys, nyn
	248	DO k = nzb_s_outer(j,i), nzt+1
	249	sums_l(k,pr,tn) = sums_l(k,pr,tn) + theta(k,j,i)
	250	ENDDO
	251	ENDDO
	252	ENDDO
	253	!$OMP END PARALLEL
	254
	255	DO i = 1, threads_per_task-1
	256	sums_l(:,pr,0) = sums_l(:,pr,0) + sums_l(:,pr,i)
	257	ENDDO
	258
	259	#if defined( __parallel )
	260
	261	CALL MPI_ALLREDUCE( sums_l(nzb,pr,0), sums(nzb,pr), nzt+2-nzb, &
	262	MPI_REAL, MPI_SUM, comm2d, ierr )
	263
	264	#else
	265
	266	sums(:,pr) = sums_l(:,pr,0)
	267
	268	#endif
	269
	270	hom(:,1,pr,0) = sums(:,pr) / ngp_2dh_outer(:,0)
	271
	272	ENDIF
	273
	274	END SUBROUTINE calc_mean_pt_profile
	275
	276	END MODULE buoyancy_mod

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

Context Navigation

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

Download in other formats: