Home

Context Navigation

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

Last change on this file since 883 was 484, checked in by raasch, 15 years ago
typo in file headers removed
Property svn:keywords set to `Id`
File size: 15.0 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |