Home

Context Navigation

source: palm/trunk/SOURCE/spline_z.f90 @ 29

Last change on this file since 29 was 19, checked in by raasch, 18 years ago
preliminary version of modified boundary conditions at top
Property svn:keywords set to `Id`
File size: 14.9 KB

Rev	Line
[1]	1	SUBROUTINE spline_z( vad_in_out, ad_v, dz_spline, spline_tri, var_char )
	2
	3	!------------------------------------------------------------------------------!
	4	! Actual revisions:
	5	! -----------------
[19]	6	! Boundary condition for pt at top adjusted
[1]	7	!
	8	! Former revisions:
	9	! -----------------
[3]	10	! $Id: spline_z.f90 19 2007-02-23 04:53:48Z raasch $
	11	! RCS Log replace by Id keyword, revision history cleaned up
	12	!
[1]	13	! Revision 1.9 2005/06/29 08:22:56 steinfeld
	14	! Dependency of ug and vg on height considered in the determination of the
	15	! upper boundary condition for vad
	16	!
	17	! Revision 1.1 1999/02/05 09:17:16 raasch
	18	! Initial revision
	19	!
	20	!
	21	! Description:
	22	! ------------
	23	! Upstream-spline advection along x
	24	!
	25	! Input/output parameters:
	26	! ad_v = advecting wind speed component
	27	! dz_spline = vertical grid spacing (dzu or dzw, depending on quantity to be
	28	! advected)
	29	! spline_tri = grid spacing factors (spl_tri_zu or spl_tri_zw, depending on
	30	! quantity to be advected)
	31	! vad_in_out = quantity to be advected, excluding ghost- or cyclic boundaries
	32	! result is given to the calling routine in this array
	33	! var_char = string which defines the quantity to be advected
	34	!
	35	! Internal arrays:
	36	! r = 2D-working array (right hand side of linear equation, buffer for
	37	! Long filter)
	38	! tf = tendency field (2D), used for long filter
	39	! vad = quantity to be advected (2D), including ghost- or cyclic
	40	! boundarys along the direction of advection
	41	! wrk_long = working array (long coefficients)
	42	! wrk_spline = working array (spline coefficients)
	43	!------------------------------------------------------------------------------!
	44
	45	USE arrays_3d
	46	USE grid_variables
	47	USE indices
	48	USE statistics
	49	USE control_parameters
	50	USE transpose_indices
	51
	52	IMPLICIT NONE
	53
	54	CHARACTER (LEN=*) :: var_char
	55
	56	INTEGER :: component, i, j, k, sr
	57	REAL :: dzwd, dzwu, overshoot_limit, t1, t2, t3, ups_limit
	58	REAL :: dz_spline(1:nzt+1)
	59	REAL :: spline_tri(5,nzb:nzt+1)
	60	REAL :: ad_v(nzb+1:nzta,nys:nyna,nxl:nxra)
	61
	62	REAL, DIMENSION(:,:), ALLOCATABLE :: r, tf, vad, wrk_spline
	63	REAL, DIMENSION(:,:,:), ALLOCATABLE :: wrk_long
	64
	65	#if defined( __parallel )
	66	REAL :: vad_in_out(nzb+1:nzta,nys:nyna,nxl:nxra)
	67	#else
	68	REAL :: vad_in_out(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1)
	69	#endif
	70
	71	!
	72	!-- Set criteria for switching between upstream- and upstream-spline-method
	73	IF ( var_char == 'u' ) THEN
	74	overshoot_limit = overshoot_limit_u
	75	ups_limit = ups_limit_u
	76	component = 1
	77	ELSEIF ( var_char == 'v' ) THEN
	78	overshoot_limit = overshoot_limit_v
	79	ups_limit = ups_limit_v
	80	component = 2
	81	ELSEIF ( var_char == 'w' ) THEN
	82	overshoot_limit = overshoot_limit_w
	83	ups_limit = ups_limit_w
	84	component = 3
	85	ELSEIF ( var_char == 'pt' ) THEN
	86	overshoot_limit = overshoot_limit_pt
	87	ups_limit = ups_limit_pt
	88	component = 4
	89	ELSEIF ( var_char == 'e' ) THEN
	90	overshoot_limit = overshoot_limit_e
	91	ups_limit = ups_limit_e
	92	component = 5
	93	ENDIF
	94
	95	!
	96	!-- Allocate working arrays
	97	ALLOCATE( r(nzb:nzt+1,nys:nyn), vad(nzb:nzt+1,nys:nyn), &
	98	wrk_spline(nzb:nzt+1,nys:nyn) )
	99	IF ( long_filter_factor /= 0.0 ) THEN
	100	ALLOCATE( tf(nzb:nzt+1,nys:nyn), wrk_long(nzb+1:nzt,nys:nyn,1:3) )
	101	ENDIF
	102
	103	!
	104	!-- Initialize calculation of relative upstream fraction
	105	sums_up_fraction_l(component,3,:) = 0.0
	106
	107	!
	108	!-- Loop over all gridpoints along x
	109	DO i = nxl, nxr
	110
	111	!
	112	!-- Store array to be advected on work array
	113	vad(nzb+1:nzt,:) = vad_in_out(nzb+1:nzt,nys:nyn,i)
	114	!
	115	!-- Add boundary conditions along z
	116	IF ( var_char == 'u' .OR. var_char == 'v' ) THEN
	117	!
	118	!-- Bottom boundary
	119	!-- u- and v-component
	120	IF ( ibc_uv_b == 0 ) THEN
	121	vad(nzb,:) = -vad(nzb+1,:)
	122	ELSE
	123	vad(nzb,:) = vad(nzb+1,:)
	124	ENDIF
	125	!
	126	!-- Top boundary
	127	!-- Dirichlet condition
	128	IF ( ibc_uv_t == 0 .AND. var_char == 'u' ) THEN
	129	!
	130	!-- u-component
	131	vad(nzt+1,:) = ug(nzt+1)
	132	ELSEIF ( ibc_uv_t == 0 .AND. var_char == 'v' ) THEN
	133	!
	134	!-- v-component
	135	vad(nzt+1,:) = vg(nzt+1)
	136	ELSE
	137	!
	138	!-- Neumann condition
	139	vad(nzt+1,:) = vad(nzt,:)
	140	ENDIF
	141
	142	ELSEIF ( var_char == 'w' ) THEN
	143	!
	144	!-- Bottom and top boundary for w-component
	145	vad(nzb,:) = 0.0
	146	vad(nzt+1,:) = 0.0
	147
	148	ELSEIF ( var_char == 'pt' ) THEN
	149	!
	150	!-- Bottom boundary for temperature
	151	IF ( ibc_pt_b == 1 ) THEN
	152	vad(nzb,:) = vad(nzb+1,:)
	153	ELSE
	154	vad(nzb,:) = pt(nzb,:,i)
	155	ENDIF
	156	!
	157	!-- Top boundary for temperature
[19]	158	IF ( ibc_pt_t == 0 ) THEN
	159	vad(nzt+1,:) = pt(nzt+1,nys:nyn,i)
	160	ELSEIF ( ibc_pt_t == 1 ) THEN
	161	vad(nzt+1,:) = vad(nzt,:)
	162	ELSEIF ( ibc_pt_t == 2 ) THEN
[1]	163	vad(nzt+1,:) = vad(nzt,:) + bc_pt_t_val * dz_spline(nzt+1)
	164	ENDIF
	165
	166	ELSEIF ( var_char == 'e' ) THEN
	167	!
	168	!-- Boundary conditions for TKE (Neumann in any case)
	169	vad(nzb,:) = vad(nzb+1,:)
	170	vad(nzt,:) = vad(nzt-1,:)
	171	vad(nzt+1,:) = vad(nzt,:)
	172
	173	ENDIF
	174
	175	!
	176	!-- Calculate right hand side
	177	DO j = nys, nyn
	178	r(nzb,j) = 3.0 * ( vad(nzb+1,j)-vad(nzb,j) ) / dz_spline(1)
	179	r(nzt+1,j) = 3.0 * ( vad(nzt+1,j)-vad(nzt,j) ) / dz_spline(nzt+1)
	180	DO k = nzb+1, nzt
	181	r(k,j) = 3.0 * ( &
	182	spline_tri(2,k) * ( vad(k,j)-vad(k-1,j) ) / dz_spline(k) &
	183	+ spline_tri(3,k) * ( vad(k+1,j)-vad(k,j) ) / dz_spline(k+1) &
	184	)
	185	ENDDO
	186	ENDDO
	187
	188	!
	189	!-- Forward substitution
	190	DO j = nys, nyn
	191	wrk_spline(nzb,j) = r(nzb,j)
	192	DO k = nzb+1, nzt+1
	193	wrk_spline(k,j) = r(k,j) - spline_tri(5,k) * r(k-1,j)
	194	ENDDO
	195	ENDDO
	196
	197	!
	198	!-- Backward substitution
	199	DO j = nys, nyn
	200	r(nzt+1,j) = wrk_spline(nzt+1,j) / spline_tri(4,nzt+1)
	201	DO k = nzt, nzb, -1
	202	r(k,j) = ( wrk_spline(k,j) - spline_tri(3,k) * r(k+1,j) ) / &
	203	spline_tri(4,k)
	204	ENDDO
	205	ENDDO
	206
	207	!
	208	!-- Calculate advection along z
	209	DO j = nys, nyn
	210	DO k = nzb+1, nzt
	211
	212	IF ( ad_v(k,j,i) == 0.0 ) THEN
	213
	214	vad_in_out(k,j,i) = vad(k,j)
	215
	216	ELSEIF ( ad_v(k,j,i) > 0.0 ) THEN
	217
	218	IF ( ABS( vad(k,j) - vad(k-1,j) ) <= ups_limit ) THEN
	219	vad_in_out(k,j,i) = vad(k,j) - dt_3d * ad_v(k,j,i) * &
	220	( vad(k,j) - vad(k-1,j) ) * ddzu(k)
	221	!
	222	!-- Calculate upstream fraction in % (s. flow_statistics)
	223	DO sr = 0, statistic_regions
	224	sums_up_fraction_l(component,3,sr) = &
	225	sums_up_fraction_l(component,3,sr) + 1.0
	226	ENDDO
	227	ELSE
	228	t1 = ad_v(k,j,i) * dt_3d / dz_spline(k)
	229	t2 = 3.0 * ( vad(k-1,j) - vad(k,j) ) + &
	230	( 2.0 * r(k,j) + r(k-1,j) ) * dz_spline(k)
	231	t3 = 2.0 * ( vad(k-1,j) - vad(k,j) ) + &
	232	( r(k,j) + r(k-1,j) ) * dz_spline(k)
	233	vad_in_out(k,j,i) = vad(k,j) - r(k,j) * t1* dz_spline(k) + &
	234	t2 * t1*2 - t3 t1**3
	235	IF ( vad(k-1,j) == vad(k,j) ) THEN
	236	vad_in_out(k,j,i) = vad(k,j)
	237	ENDIF
	238	ENDIF
	239
	240	ELSE
	241
	242	IF( ABS( vad(k,j) - vad(k+1,j) ) <= ups_limit ) THEN
	243	vad_in_out(k,j,i) = vad(k,j) - dt_3d * ad_v(k,j,i) * &
	244	( vad(k+1,j) - vad(k,j) ) * ddzu(k+1)
	245	!
	246	!-- Calculate upstream fraction in % (s. flow_statistics)
	247	DO sr = 0, statistic_regions
	248	sums_up_fraction_l(component,3,sr) = &
	249	sums_up_fraction_l(component,3,sr) + 1.0
	250	ENDDO
	251	ELSE
	252	t1 = -ad_v(k,j,i) * dt_3d / dz_spline(k+1)
	253	t2 = 3.0 * ( vad(k,j) - vad(k+1,j) ) + &
	254	( 2.0 * r(k,j) + r(k+1,j) ) * dz_spline(k+1)
	255	t3 = 2.0 * ( vad(k,j) - vad(k+1,j) ) + &
	256	( r(k,j) + r(k+1,j) ) * dz_spline(k+1)
	257	vad_in_out(k,j,i) = vad(k,j) + r(k,j)t1dz_spline(k+1) - &
	258	t2 * t1*2 + t3 t1**3
	259	IF ( vad(k+1,j) == vad(k,j) ) THEN
	260	vad_in_out(k,j,i) = vad(k,j)
	261	ENDIF
	262	ENDIF
	263
	264	ENDIF
	265	ENDDO
	266	ENDDO
	267
	268	!
	269	!-- Limit values in order to prevent overshooting
	270	IF ( cut_spline_overshoot ) THEN
	271
	272	DO j = nys, nyn
	273	DO k = nzb+1, nzt
	274	IF ( ad_v(k,j,i) > 0.0 ) THEN
	275	IF ( vad(k,j) > vad(k-1,j) ) THEN
	276	vad_in_out(k,j,i) = MIN( vad_in_out(k,j,i), &
	277	vad(k,j) + overshoot_limit )
	278	vad_in_out(k,j,i) = MAX( vad_in_out(k,j,i), &
	279	vad(k-1,j) - overshoot_limit )
	280	ELSE
	281	vad_in_out(k,j,i) = MAX( vad_in_out(k,j,i), &
	282	vad(k,j) - overshoot_limit )
	283	vad_in_out(k,j,i) = MIN( vad_in_out(k,j,i), &
	284	vad(k-1,j) + overshoot_limit )
	285	ENDIF
	286	ELSE
	287	IF ( vad(k,j) > vad(k+1,j) ) THEN
	288	vad_in_out(k,j,i) = MIN( vad_in_out(k,j,i), &
	289	vad(k,j) + overshoot_limit )
	290	vad_in_out(k,j,i) = MAX( vad_in_out(k,j,i), &
	291	vad(k+1,j) - overshoot_limit )
	292	ELSE
	293	vad_in_out(k,j,i) = MAX( vad_in_out(k,j,i), &
	294	vad(k,j) - overshoot_limit )
	295	vad_in_out(k,j,i) = MIN( vad_in_out(k,j,i), &
	296	vad(k+1,j) + overshoot_limit )
	297	ENDIF
	298	ENDIF
	299	ENDDO
	300	ENDDO
	301
	302	ENDIF
	303
	304	!
	305	!-- Long-filter (acting on tendency only)
	306	IF ( long_filter_factor /= 0.0 ) THEN
	307
	308	!
	309	!-- Compute tendency
	310	DO j = nys, nyn
	311
	312	!
	313	!-- Depending on the quantity to be advected, the respective vertical
	314	!-- boundary conditions must be applied.
	315	IF ( var_char == 'u' .OR. var_char == 'v' ) THEN
	316
	317	IF ( ibc_uv_b == 0 ) THEN
	318	tf(nzb,j) = - ( vad_in_out(nzb+1,j,i) - vad(nzb+1,j) )
	319	ELSE
	320	tf(nzb,j) = vad_in_out(nzb+1,j,i) - vad(nzb+1,j)
	321	ENDIF
	322
	323	IF ( ibc_uv_t == 0 ) THEN
	324	tf(nzt+1,j) = 0.0
	325	ELSE
	326	tf(nzt+1,j) = vad_in_out(nzt,j,i) - vad(nzt,j)
	327	ENDIF
	328
	329	ELSEIF ( var_char == 'w' ) THEN
	330
	331	tf(nzb,j) = 0.0
	332	tf(nzt+1,j) = 0.0
	333
	334	ELSEIF ( var_char == 'pt' ) THEN
	335
	336	IF ( ibc_pt_b == 1 ) THEN
	337	tf(nzb,j) = vad_in_out(nzb+1,j,i) - vad(nzb+1,j)
	338	ELSE
	339	tf(nzb,j) = 0.0
	340	ENDIF
	341
	342	IF ( ibc_pt_t == 1 ) THEN
	343	vad_in_out(nzt,j,i) = vad_in_out(nzt-1,j,i) + bc_pt_t_val * &
	344	dz_spline(nzt)
	345	tf(nzt+1,j) = vad_in_out(nzt,j,i) + bc_pt_t_val * &
	346	dz_spline(nzt+1) - vad(nzt+1,j)
	347	ELSE
	348	vad_in_out(nzt,j,i) = pt(nzt,j,i)
	349	tf(nzt+1,j) = 0.0
	350	ENDIF
	351
	352	ENDIF
	353
	354	DO k = nzb+1, nzt
	355	tf(k,j) = vad_in_out(k,j,i) - vad(k,j)
	356	ENDDO
	357
	358	ENDDO
	359
	360	!
	361	!-- Apply the filter.
	362	DO j = nys, nyn
	363
	364	dzwd = dz_spline(1) / ( dz_spline(1) + dz_spline(2) )
	365	dzwu = dz_spline(2) / ( dz_spline(1) + dz_spline(2) )
	366
	367	wrk_long(nzb+1,j,1) = 2.0 * ( 1.0 + long_filter_factor )
	368	wrk_long(nzb+1,j,2) = ( 1.0 - long_filter_factor ) * dzwd / &
	369	wrk_long(nzb+1,j,1)
	370	wrk_long(nzb+1,j,3) = ( long_filter_factor * dzwu * tf(nzb,j) + &
	371	2.0 * tf(nzb+1,j) + dzwd * tf(nzb+2,j) &
	372	) / wrk_long(nzb+1,j,1)
	373
	374	DO k = nzb+2, nzt-1
	375
	376	dzwd = dz_spline(k) / ( dz_spline(k) + dz_spline(k+1) )
	377	dzwu = dz_spline(k+1) / ( dz_spline(k) + dz_spline(k+1) )
	378
	379	wrk_long(k,j,1) = 2.0 * ( 1.0 + long_filter_factor ) - &
	380	( 1.0 - long_filter_factor ) * dzwu * &
	381	wrk_long(k-1,j,2)
	382	wrk_long(k,j,2) = ( 1.0 - long_filter_factor ) * dzwd / &
	383	wrk_long(k,j,1)
	384	wrk_long(k,j,3) = ( dzwu * tf(k-1,j) + 2.0 * tf(k,j) + &
	385	dzwd * tf(k+1,j) - &
	386	( 1.0 - long_filter_factor ) * dzwu * &
	387	wrk_long(k-1,j,3) &
	388	) / wrk_long(k,j,1)
	389	ENDDO
	390
	391	dzwd = dz_spline(nzt) / ( dz_spline(nzt) + dz_spline(nzt+1) )
	392	dzwu = dz_spline(nzt+1) / ( dz_spline(nzt) + dz_spline(nzt+1) )
	393
	394	wrk_long(nzt,j,1) = 2.0 * ( 1.0 + long_filter_factor ) - &
	395	( 1.0 - long_filter_factor ) * dzwu * &
	396	wrk_long(nzt-1,j,2)
	397	wrk_long(nzt,j,2) = ( 1.0 - long_filter_factor ) * dzwd / &
	398	wrk_long(nzt,j,1)
	399	wrk_long(nzt,j,3) = ( dzwu * tf(nzt-1,j) + 2.0 * tf(nzt,j) + &
	400	dzwd * long_filter_factor * tf(nzt+1,j) - &
	401	( 1.0 - long_filter_factor ) * dzwu * &
	402	wrk_long(nzt-1,j,3) &
	403	) / wrk_long(nzt,j,1)
	404	r(nzt,j) = wrk_long(nzt,j,3)
	405
	406	ENDDO
	407
	408	DO j = nys, nyn
	409	DO k = nzt-1, nzb+1, -1
	410	r(k,j) = wrk_long(k,j,3) - wrk_long(k,j,2) * r(k+1,j)
	411	ENDDO
	412	ENDDO
	413
	414	DO j = nys, nyn
	415	DO k = nzb+1, nzt
	416	vad_in_out(k,j,i) = vad(k,j) + r(k,j)
	417	ENDDO
	418	ENDDO
	419
	420	ENDIF ! Long filter
	421
	422	ENDDO
	423
	424	DEALLOCATE( r, vad, wrk_spline )
	425	IF ( long_filter_factor /= 0.0 ) DEALLOCATE( tf, wrk_long )
	426
	427	END SUBROUTINE spline_z

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |