Home

Context Navigation

source: palm/trunk/SOURCE/spline_y.f90 @ 354

Last change on this file since 354 was 4, checked in by raasch, 18 years ago
Id keyword set as property for all *.f90 files
Property svn:keywords set to `Id`
File size: 19.0 KB

Rev	Line
[1]	1	SUBROUTINE spline_y( vad_in_out, ad_v, var_char )
	2
	3	!------------------------------------------------------------------------------!
	4	! Actual revisions:
	5	! -----------------
	6	!
	7	!
	8	! Former revisions:
	9	! -----------------
[3]	10	! $Id: spline_y.f90 4 2007-02-13 11:33:16Z maronga $
	11	! RCS Log replace by Id keyword, revision history cleaned up
	12	!
[1]	13	! Revision 1.9 2004/04/30 12:54:37 raasch
	14	! Names of transpose indices changed, enlarged transposition arrays introduced
	15	!
	16	! Revision 1.1 1999/02/05 09:16:31 raasch
	17	! Initial revision
	18	!
	19	!
	20	! Description:
	21	! ------------
	22	! Upstream-spline advection along x
	23	!
	24	! Input/output parameters:
	25	! ad_v = advecting wind speed component
	26	! vad_in_out = quantity to be advected, excluding ghost- or cyclic boundaries
	27	! result is given to the calling routine in this array
	28	! var_char = string which defines the quantity to be advected
	29	!
	30	! Internal arrays:
	31	! r = 2D-working array (right hand side of linear equation, buffer for
	32	! Long filter)
	33	! tf = tendency field (2D), used for long filter
	34	! vad = quantity to be advected (2D), including ghost- or cyclic
	35	! boundarys along the direction of advection
	36	! wrk_long = working array (long coefficients)
	37	! wrk_spline = working array (spline coefficients)
	38	!------------------------------------------------------------------------------!
	39
	40	USE advection
	41	USE grid_variables
	42	USE indices
	43	USE statistics
	44	USE control_parameters
	45	USE transpose_indices
	46
	47	IMPLICIT NONE
	48
	49	CHARACTER (LEN=*) :: var_char
	50
	51	INTEGER :: component, i, j, k, sr
	52	REAL :: overshoot_limit, sm_faktor, t1, t2, t3, ups_limit
	53	REAL, DIMENSION(:,:), ALLOCATABLE :: r, tf, vad, wrk_spline
	54	REAL, DIMENSION(:,:,:), ALLOCATABLE :: wrk_long
	55
	56	#if defined( __parallel )
	57	REAL :: ad_v(0:nya,nxl_y:nxr_ya,nzb_y:nzt_ya), &
	58	vad_in_out(0:nya,nxl_y:nxr_ya,nzb_y:nzt_ya)
	59	#else
	60	REAL :: ad_v(nzb+1:nzt,nys:nyn,nxl:nxr), &
	61	vad_in_out(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1)
	62	#endif
	63
	64	!
	65	!-- Set criteria for switching between upstream- and upstream-spline-method
	66	IF ( var_char == 'u' ) THEN
	67	overshoot_limit = overshoot_limit_u
	68	ups_limit = ups_limit_u
	69	component = 1
	70	ELSEIF ( var_char == 'v' ) THEN
	71	overshoot_limit = overshoot_limit_v
	72	ups_limit = ups_limit_v
	73	component = 2
	74	ELSEIF ( var_char == 'w' ) THEN
	75	overshoot_limit = overshoot_limit_w
	76	ups_limit = ups_limit_w
	77	component = 3
	78	ELSEIF ( var_char == 'pt' ) THEN
	79	overshoot_limit = overshoot_limit_pt
	80	ups_limit = ups_limit_pt
	81	component = 4
	82	ELSEIF ( var_char == 'e' ) THEN
	83	overshoot_limit = overshoot_limit_e
	84	ups_limit = ups_limit_e
	85	component = 5
	86	ENDIF
	87
	88	!
	89	!-- Initialize calculation of relative upstream fraction
	90	sums_up_fraction_l(component,2,:) = 0.0
	91
	92	#if defined( __parallel )
	93
	94	!
	95	!-- Allocate working arrays
	96	ALLOCATE( r(-1:ny+1,nxl_y:nxr_y), &
	97	vad(-1:ny+1,nxl_y:nxr_y), &
	98	wrk_spline(0:ny,nxl_y:nxr_y) )
	99	IF ( long_filter_factor /= 0.0 ) THEN
	100	ALLOCATE( tf(0:ny,nxl_y:nxr_y), &
	101	wrk_long(0:ny,nxl_y:nxr_y,1:3) )
	102	ENDIF
	103
	104	!
	105	!-- Loop over all gridpoints along z
	106	DO k = nzb_y, nzt_y
	107
	108	!
	109	!-- Store array to be advected on work array and add cyclic boundary along y
	110	vad(0:ny,nxl_y:nxr_y) = vad_in_out(0:ny,nxl_y:nxr_y,k)
	111	vad(-1,:) = vad(ny,:)
	112	vad(ny+1,:) = vad(0,:)
	113
	114	!
	115	!-- Calculate right hand side
	116	DO i = nxl_y, nxr_y
	117	DO j = 0, ny
	118	r(j,i) = 3.0 * ( &
	119	spl_tri_y(2,j) * ( vad(j,i) - vad(j-1,i) ) * ddy + &
	120	spl_tri_y(3,j) * ( vad(j+1,i) - vad(j,i) ) * ddy &
	121	)
	122	ENDDO
	123	ENDDO
	124
	125	!
	126	!-- Forward substitution
	127	DO i = nxl_y, nxr_y
	128	wrk_spline(0,i) = r(0,i)
	129	DO j = 1, ny
	130	wrk_spline(j,i) = r(j,i) - spl_tri_y(5,j) * wrk_spline(j-1,i)
	131	ENDDO
	132	ENDDO
	133
	134	!
	135	!-- Backward substitution (sherman-Morrison-formula)
	136	DO i = nxl_y, nxr_y
	137	r(ny,i) = wrk_spline(ny,i) / spl_tri_y(4,ny)
	138	DO j = ny-1, 0, -1
	139	r(j,i) = ( wrk_spline(j,i) - spl_tri_y(3,j) * r(j+1,i) ) / &
	140	spl_tri_y(4,j)
	141	ENDDO
	142	sm_faktor = ( r(0,i) + 0.5 * r(ny,i) / spl_gamma_y ) / &
	143	( 1.0 + spl_z_y(0) + 0.5 * spl_z_y(ny) / spl_gamma_y )
	144	DO j = 0, ny
	145	r(j,i) = r(j,i) - sm_faktor * spl_z_y(j)
	146	ENDDO
	147	ENDDO
	148
	149	!
	150	!-- Add cyclic boundary conditions to right hand side
	151	r(-1,:) = r(ny,:)
	152	r(ny+1,:) = r(0,:)
	153
	154	!
	155	!-- Calculate advection along y
	156	DO i = nxl_y, nxr_y
	157	DO j = 0, ny
	158
	159	IF ( ad_v(j,i,k) == 0.0 ) THEN
	160
	161	vad_in_out(j,i,k) = vad(j,i)
	162
	163	ELSEIF ( ad_v(j,i,k) > 0.0 ) THEN
	164
	165	IF ( ABS( vad(j,i) - vad(j-1,i) ) <= ups_limit ) THEN
	166	vad_in_out(j,i,k) = vad(j,i) - dt_3d * ad_v(j,i,k) * &
	167	( vad(j,i) - vad(j-1,i) ) * ddy
	168	!
	169	!-- Calculate upstream fraction in % (s. flow_statistics)
	170	DO sr = 0, statistic_regions
	171	sums_up_fraction_l(component,2,sr) = &
	172	sums_up_fraction_l(component,2,sr) + 1.0
	173	ENDDO
	174	ELSE
	175	t1 = ad_v(j,i,k) * dt_3d * ddy
	176	t2 = 3.0 * ( vad(j-1,i) - vad(j,i) ) + &
	177	( 2.0 * r(j,i) + r(j-1,i) ) * dy
	178	t3 = 2.0 * ( vad(j-1,i) - vad(j,i) ) + &
	179	( r(j,i) + r(j-1,i) ) * dy
	180	vad_in_out(j,i,k) = vad(j,i) - r(j,i) * t1 * dy + &
	181	t2 * t1*2 - t3 t1**3
	182	IF ( vad(j-1,i) == vad(j,i) ) THEN
	183	vad_in_out(j,i,k) = vad(j,i)
	184	ENDIF
	185	ENDIF
	186
	187	ELSE
	188
	189	IF ( ABS( vad(j,i) - vad(j+1,i) ) <= ups_limit ) THEN
	190	vad_in_out(j,i,k) = vad(j,i) - dt_3d * ad_v(j,i,k) * &
	191	( vad(j+1,i) - vad(j,i) ) * ddy
	192	!
	193	!-- Calculate upstream fraction in % (s. flow_statistics)
	194	DO sr = 0, statistic_regions
	195	sums_up_fraction_l(component,2,sr) = &
	196	sums_up_fraction_l(component,2,sr) + 1.0
	197	ENDDO
	198	ELSE
	199	t1 = -ad_v(j,i,k) * dt_3d * ddy
	200	t2 = 3.0 * ( vad(j,i) - vad(j+1,i) ) + &
	201	( 2.0 * r(j,i) + r(j+1,i) ) * dy
	202	t3 = 2.0 * ( vad(j,i) - vad(j+1,i) ) + &
	203	( r(j,i) + r(j+1,i) ) * dy
	204	vad_in_out(j,i,k) = vad(j,i) + r(j,i) * t1 * dy - &
	205	t2 * t1*2 + t3 t1**3
	206	IF ( vad(j+1,i) == vad(j,i) ) THEN
	207	vad_in_out(j,i,k) = vad(j,i)
	208	ENDIF
	209	ENDIF
	210
	211	ENDIF
	212	ENDDO
	213	ENDDO
	214
	215	!
	216	!-- Limit values in order to prevent overshooting
	217	IF ( cut_spline_overshoot ) THEN
	218
	219	DO i = nxl_y, nxr_y
	220	DO j = 0, ny
	221	IF ( ad_v(j,i,k) > 0.0 ) THEN
	222	IF ( vad(j,i) > vad(j-1,i) ) THEN
	223	vad_in_out(j,i,k) = MIN( vad_in_out(j,i,k), &
	224	vad(j,i) + overshoot_limit )
	225	vad_in_out(j,i,k) = MAX( vad_in_out(j,i,k), &
	226	vad(j-1,i) - overshoot_limit )
	227	ELSE
	228	vad_in_out(j,i,k) = MAX( vad_in_out(j,i,k), &
	229	vad(j,i) - overshoot_limit )
	230	vad_in_out(j,i,k) = MIN( vad_in_out(j,i,k), &
	231	vad(j-1,i) + overshoot_limit )
	232	ENDIF
	233	ELSE
	234	IF ( vad(j,i) > vad(j+1,i) ) THEN
	235	vad_in_out(j,i,k) = MIN( vad_in_out(j,i,k), &
	236	vad(j,i) + overshoot_limit )
	237	vad_in_out(j,i,k) = MAX( vad_in_out(j,i,k), &
	238	vad(j+1,i) - overshoot_limit )
	239	ELSE
	240	vad_in_out(j,i,k) = MAX( vad_in_out(j,i,k), &
	241	vad(j,i) - overshoot_limit )
	242	vad_in_out(j,i,k) = MIN( vad_in_out(j,i,k), &
	243	vad(j+1,i) + overshoot_limit )
	244	ENDIF
	245	ENDIF
	246	ENDDO
	247	ENDDO
	248
	249	ENDIF
	250
	251	!
	252	!-- Long-filter (acting on tendency only)
	253	IF ( long_filter_factor /= 0.0 ) THEN
	254
	255	!
	256	!-- Compute tendency. Filter only acts on this quantity.
	257	DO i = nxl_y, nxr_y
	258	DO j = 0, ny
	259	tf(j,i) = vad_in_out(j,i,k) - vad(j,i)
	260	ENDDO
	261	ENDDO
	262
	263	!
	264	!-- Apply the filter.
	265	DO i = nxl_y, nxr_y
	266	wrk_long(0,i,1) = 2.0 * ( 1.0 + long_filter_factor )
	267	wrk_long(0,i,2) = ( 1.0 - long_filter_factor ) / wrk_long(0,i,1)
	268	wrk_long(0,i,3) = ( long_filter_factor * tf(ny,i) + &
	269	2.0 * tf(0,i) + tf(1,i) &
	270	) / wrk_long(0,i,1)
	271	DO j = 1, ny-1
	272	wrk_long(j,i,1) = 2.0 * ( 1.0 + long_filter_factor ) - &
	273	( 1.0 - long_filter_factor ) * wrk_long(j-1,i,2)
	274	wrk_long(j,i,2) = ( 1.0 - long_filter_factor ) / wrk_long(j,i,1)
	275	wrk_long(j,i,3) = ( tf(j-1,i) + 2.0 * tf(j,i) + &
	276	tf(j+1,i) - ( 1.0 - long_filter_factor ) * &
	277	wrk_long(j-1,i,3) ) / wrk_long(j,i,1)
	278	ENDDO
	279	wrk_long(ny,i,1) = 2.0 * ( 1.0 + long_filter_factor ) - &
	280	( 1.0 - long_filter_factor ) * wrk_long(ny-1,i,2)
	281	wrk_long(ny,i,2) = ( 1.0 - long_filter_factor ) / wrk_long(ny,i,1)
	282	wrk_long(ny,i,3) = ( tf(ny-1,i) + 2.0 * tf(ny,i) + &
	283	long_filter_factor * tf(0,i) - &
	284	( 1.0 - long_filter_factor ) * &
	285	wrk_long(ny-1,i,3) &
	286	) / wrk_long(ny,i,1)
	287	r(ny,i) = wrk_long(ny,i,3)
	288	ENDDO
	289
	290	DO j = ny-1, 0, -1
	291	DO i = nxl_y, nxr_y
	292	r(j,i) = wrk_long(j,i,3) - wrk_long(j,i,2) * r(j+1,i)
	293	ENDDO
	294	ENDDO
	295
	296	DO i = nxl_y, nxr_y
	297	DO j = 0, ny
	298	vad_in_out(j,i,k) = vad(j,i) + r(j,i)
	299	ENDDO
	300	ENDDO
	301
	302	ENDIF ! Long filter
	303
	304	ENDDO
	305
	306	#else
	307
	308	!
	309	!-- Allocate working arrays
	310	ALLOCATE( r(nzb+1:nzt,nys-1:nyn+1), vad(nzb:nzt+1,nys-1:nyn+1), &
	311	wrk_spline(nzb+1:nzt,nys-1:nyn+1) )
	312	IF ( long_filter_factor /= 0.0 ) THEN
	313	ALLOCATE( tf(nzb+1:nzt,nys-1:nyn+1), wrk_long(nzb+1:nzt,0:ny,1:3) )
	314	ENDIF
	315
	316	!
	317	!-- Loop over all gridpoints along x
	318	DO i = nxl, nxr
	319
	320	!
	321	!-- Store array to be advected on work array and add cyclic boundary along x
	322	vad(:,:) = vad_in_out(:,:,i)
	323	vad(:,-1) = vad(:,ny)
	324	vad(:,ny+1) = vad(:,0)
	325
	326	!
	327	!-- Calculate right hand side
	328	DO j = 0, ny
	329	DO k = nzb+1, nzt
	330	r(k,j) = 3.0 * ( &
	331	spl_tri_y(2,j) * ( vad(k,j) - vad(k,j-1) ) * ddy + &
	332	spl_tri_y(3,j) * ( vad(k,j+1) - vad(k,j) ) * ddy &
	333	)
	334	ENDDO
	335	ENDDO
	336
	337	!
	338	!-- Forward substitution
	339	DO k = nzb+1, nzt
	340	wrk_spline(k,0) = r(k,0)
	341	ENDDO
	342
	343	DO j = 1, ny
	344	DO k = nzb+1, nzt
	345	wrk_spline(k,j) = r(k,j) - spl_tri_y(5,j) * wrk_spline(k,j-1)
	346	ENDDO
	347	ENDDO
	348
	349	!
	350	!-- Backward substitution (Sherman-Morrison-formula)
	351	DO k = nzb+1, nzt
	352	r(k,ny) = wrk_spline(k,ny) / spl_tri_y(4,ny)
	353	ENDDO
	354
	355	DO k = nzb+1, nzt
	356	DO j = ny-1, 0, -1
	357	r(k,j) = ( wrk_spline(k,j) - spl_tri_y(3,j) * r(k,j+1) ) / &
	358	spl_tri_y(4,j)
	359	ENDDO
	360	sm_faktor = ( r(k,0) + 0.5 * r(k,ny) / spl_gamma_y ) / &
	361	( 1.0 + spl_z_y(0) + 0.5 * spl_z_y(ny) / spl_gamma_y )
	362	DO j = 0, ny
	363	r(k,j) = r(k,j) - sm_faktor * spl_z_y(j)
	364	ENDDO
	365	ENDDO
	366
	367	!
	368	!-- Add cyclic boundary to the right hand side
	369	r(:,-1) = r(:,ny)
	370	r(:,ny+1) = r(:,0)
	371
	372	!
	373	!-- Calculate advection along y
	374	DO j = 0, ny
	375	DO k = nzb+1, nzt
	376
	377	IF ( ad_v(k,j,i) == 0.0 ) THEN
	378
	379	vad_in_out(k,j,i) = vad(k,j)
	380
	381	ELSEIF ( ad_v(k,j,i) > 0.0 ) THEN
	382
	383	IF ( ABS( vad(k,j) - vad(k,j-1) ) <= ups_limit ) THEN
	384	vad_in_out(k,j,i) = vad(k,j) - dt_3d * ad_v(k,j,i) * &
	385	( vad(k,j) - vad(k,j-1) ) * ddy
	386	!
	387	!-- Calculate upstream fraction in % (s. flow_statistics)
	388	DO sr = 0, statistic_regions
	389	sums_up_fraction_l(component,2,sr) = &
	390	sums_up_fraction_l(component,2,sr) + 1.0
	391	ENDDO
	392	ELSE
	393	t1 = ad_v(k,j,i) * dt_3d * ddy
	394	t2 = 3.0 * ( vad(k,j-1) - vad(k,j) ) + &
	395	( 2.0 * r(k,j) + r(k,j-1) ) * dy
	396	t3 = 2.0 * ( vad(k,j-1) - vad(k,j) ) + &
	397	( r(k,j) + r(k,j-1) ) * dy
	398	vad_in_out(k,j,i) = vad(k,j) - r(k,j) * t1 * dy + &
	399	t2 * t1*2 - t3 t1**3
	400	IF ( vad(k,j-1) == vad(k,j) ) THEN
	401	vad_in_out(k,j,i) = vad(k,j)
	402	ENDIF
	403	ENDIF
	404
	405	ELSE
	406
	407	IF ( ABS( vad(k,j) - vad(k,j+1) ) <= ups_limit ) THEN
	408	vad_in_out(k,j,i) = vad(k,j) - dt_3d * ad_v(k,j,i) * &
	409	( vad(k,j+1) - vad(k,j) ) * ddy
	410	!
	411	!-- Calculate upstream fraction in % (s. flow_statistics)
	412	DO sr = 0, statistic_regions
	413	sums_up_fraction_l(component,2,sr) = &
	414	sums_up_fraction_l(component,2,sr) + 1.0
	415	ENDDO
	416	ELSE
	417	t1 = -ad_v(k,j,i) * dt_3d * ddy
	418	t2 = 3.0 * ( vad(k,j) - vad(k,j+1) ) + &
	419	( 2.0 * r(k,j) + r(k,j+1) ) * dy
	420	t3 = 2.0 * ( vad(k,j) - vad(k,j+1) ) + &
	421	( r(k,j) + r(k,j+1) ) * dy
	422	vad_in_out(k,j,i) = vad(k,j) + r(k,j) * t1 * dy - &
	423	t2 * t1*2 + t3 t1**3
	424	IF ( vad(k,j+1) == vad(k,j) ) THEN
	425	vad_in_out(k,j,i) = vad(k,j)
	426	ENDIF
	427	ENDIF
	428
	429	ENDIF
	430	ENDDO
	431	ENDDO
	432
	433	!
	434	!-- Limit values in order to prevent overshooting
	435	IF ( cut_spline_overshoot ) THEN
	436
	437	DO j = 0, ny
	438	DO k = nzb+1, nzt
	439	IF ( ad_v(k,j,i) > 0.0 ) THEN
	440	IF ( vad(k,j) > vad(k,j-1) ) THEN
	441	vad_in_out(k,j,i) = MIN( vad_in_out(k,j,i), &
	442	vad(k,j) + overshoot_limit )
	443	vad_in_out(k,j,i) = MAX( vad_in_out(k,j,i), &
	444	vad(k,j-1) - overshoot_limit )
	445	ELSE
	446	vad_in_out(k,j,i) = MAX( vad_in_out(k,j,i), &
	447	vad(k,j) - overshoot_limit )
	448	vad_in_out(k,j,i) = MIN( vad_in_out(k,j,i), &
	449	vad(k,j-1) + overshoot_limit )
	450	ENDIF
	451	ELSE
	452	IF ( vad(k,j) > vad(k,j+1) ) THEN
	453	vad_in_out(k,j,i) = MIN( vad_in_out(k,j,i), &
	454	vad(k,j) + overshoot_limit )
	455	vad_in_out(k,j,i) = MAX( vad_in_out(k,j,i), &
	456	vad(k,j+1) - overshoot_limit )
	457	ELSE
	458	vad_in_out(k,j,i) = MAX( vad_in_out(k,j,i), &
	459	vad(k,j) - overshoot_limit )
	460	vad_in_out(k,j,i) = MIN( vad_in_out(k,j,i), &
	461	vad(k,j+1) + overshoot_limit )
	462	ENDIF
	463	ENDIF
	464	ENDDO
	465	ENDDO
	466
	467	ENDIF
	468
	469	!
	470	!-- Long filter (acting on tendency only)
	471	IF ( long_filter_factor /= 0.0 ) THEN
	472
	473	!
	474	!-- Compute tendency
	475	DO j = nys, nyn
	476	DO k = nzb+1, nzt
	477	tf(k,j) = vad_in_out(k,j,i) - vad(k,j)
	478	ENDDO
	479	ENDDO
	480
	481	!
	482	!-- Apply the filter
	483	wrk_long(:,0,1) = 2.0 * ( 1.0 + long_filter_factor )
	484	wrk_long(:,0,2) = ( 1.0 - long_filter_factor ) / wrk_long(:,0,1)
	485	wrk_long(:,0,3) = ( long_filter_factor * tf(:,ny) + &
	486	2.0 * tf(:,0) + tf(:,1) ) / wrk_long(:,0,1)
	487
	488	DO j = 1, ny-1
	489	DO k = nzb+1, nzt
	490	wrk_long(k,j,1) = 2.0 * ( 1.0 + long_filter_factor ) - &
	491	( 1.0 - long_filter_factor ) * &
	492	wrk_long(k,j-1,2)
	493	wrk_long(k,j,2) = ( 1.0 - long_filter_factor ) / wrk_long(k,j,1)
	494	wrk_long(k,j,3) = ( tf(k,j-1) + 2.0 * tf(k,j) + &
	495	tf(k,j+1) - ( 1.0 - long_filter_factor ) * &
	496	wrk_long(k,j-1,3) ) / wrk_long(k,j,1)
	497	ENDDO
	498	wrk_long(:,ny,1) = 2.0 * ( 1.0 + long_filter_factor ) - &
	499	( 1.0 - long_filter_factor ) * &
	500	wrk_long(:,ny-1,2)
	501	wrk_long(:,ny,2) = ( 1.0 - long_filter_factor ) / wrk_long(:,ny,1)
	502	wrk_long(:,ny,3) = ( tf(:,ny-1) + 2.0 * tf(:,ny) + &
	503	long_filter_factor * tf(:,0) - &
	504	( 1.0 - long_filter_factor ) * &
	505	wrk_long(:,ny-1,3) ) / wrk_long(:,ny,1)
	506	r(:,ny) = wrk_long(:,ny,3)
	507	ENDDO
	508	DO j = ny-1, 0, -1
	509	DO k = nzb+1, nzt
	510	r(k,j) = wrk_long(k,j,3) - wrk_long(k,j,2) * r(k,j+1)
	511	ENDDO
	512	ENDDO
	513	DO j = 0, ny
	514	DO k = nzb+1, nzt
	515	vad_in_out(k,j,i) = vad(k,j) + r(k,j)
	516	ENDDO
	517	ENDDO
	518
	519	ENDIF ! Long filter
	520
	521	ENDDO
	522	#endif
	523
	524	IF ( long_filter_factor /= 0.0 ) DEALLOCATE( tf, wrk_long )
	525	DEALLOCATE( r, vad, wrk_spline )
	526
	527	END SUBROUTINE spline_y

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |