Context Navigation

moldiff.F @ 2325

Last change on this file since 2325 was 1266, checked in by aslmd, 11 years ago

LMDZ.MARS
IMPORTANT CHANGE

Remove all reference/use of nlayermx and dimphys.h
Made use of automatic arrays whenever arrays are needed with dimension nlayer
Remove lots of obsolete reference to dimensions.h
Converted iono.h and param_v4.h into corresponding modules

(with embedded subroutine to allocate arrays)
(no arrays allocated if thermosphere not used)

Deleted param.h and put contents into module param_v4_h
Adapted testphys1d, newstart, etc...
Made DATA arrays in param_read to be initialized by subroutine

fill_data_thermos in module param_v4_h

Optimized computations in paramfoto_compact (twice less dlog10 calculations)
Checked consistency before/after modification in debug mode
Checked performance is not impacted (same as before)

File size: 17.8 KB

Line
1	subroutine moldiff(ngrid,nlayer,nq,
2	& pplay,pplev,pt,pdt,pq,pdq,ptimestep,
3	& zzlay,pdteuv,pdtconduc,pdqdiff)
4
5	use tracer_mod, only: igcm_co2, igcm_co, igcm_o, igcm_o1d,
6	& igcm_o2, igcm_o3, igcm_h, igcm_h2, igcm_oh,
7	& igcm_ho2, igcm_h2o2, igcm_n2, igcm_ar,
8	& igcm_h2o_vap, mmol
9	use conc_mod, only: rnew, mmean
10	USE comcstfi_h
11	implicit none
12
13	c
14	c Input/Output
15	c
16	integer,intent(in) :: ngrid ! number of atmospheric columns
17	integer,intent(in) :: nlayer ! number of atmospheric layers
18	integer,intent(in) :: nq ! number of advected tracers
19	real ptimestep
20	real pplay(ngrid,nlayer)
21	real zzlay(ngrid,nlayer)
22	real pplev(ngrid,nlayer+1)
23	real pq(ngrid,nlayer,nq)
24	real pdq(ngrid,nlayer,nq)
25	real pt(ngrid,nlayer)
26	real pdt(ngrid,nlayer)
27	real pdteuv(ngrid,nlayer)
28	real pdtconduc(ngrid,nlayer)
29	real pdqdiff(ngrid,nlayer,nq)
30	c
31	c Local
32	c
33
34	integer,parameter :: ncompmoldiff = 14
35	real hco2(ncompmoldiff),ho
36
37	integer ig,nz,l,n,nn,iq
38	real del1,del2, tmean ,dalfinvdz, d
39	real hh,dcoef,dcoef1,ptfac, ntot, dens, dens2, dens3
40	real hp(nlayer)
41	real tt(nlayer)
42	real qq(nlayer,ncompmoldiff)
43	real dmmeandz(nlayer)
44	real qnew(nlayer,ncompmoldiff)
45	real zlocal(nlayer)
46	real alf(ncompmoldiff-1,ncompmoldiff-1)
47	real alfinv(nlayer,ncompmoldiff-1,ncompmoldiff-1)
48	real indx(ncompmoldiff-1)
49	real b(nlayer,ncompmoldiff-1)
50	real y(ncompmoldiff-1,ncompmoldiff-1)
51	real aa(nlayer,ncompmoldiff-1,ncompmoldiff-1)
52	real bb(nlayer,ncompmoldiff-1,ncompmoldiff-1)
53	real cc(nlayer,ncompmoldiff-1,ncompmoldiff-1)
54	real atri(nlayer-2)
55	real btri(nlayer-2)
56	real ctri(nlayer-2)
57	real rtri(nlayer-2)
58	real qtri(nlayer-2)
59	real alfdiag(ncompmoldiff-1)
60	real wi(ncompmoldiff), flux(ncompmoldiff), pote
61
62	cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
63	c tracer numbering in the molecular diffusion
64	cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
65	c Atomic oxygen must always be the LAST species of the list as
66	c it is the dominant species at high altitudes.
67	integer,parameter :: i_co = 1
68	integer,parameter :: i_n2 = 2
69	integer,parameter :: i_o2 = 3
70	integer,parameter :: i_co2 = 4
71	integer,parameter :: i_h2 = 5
72	integer,parameter :: i_h = 6
73	integer,parameter :: i_oh = 7
74	integer,parameter :: i_ho2 = 8
75	integer,parameter :: i_h2o = 9
76	integer,parameter :: i_h2o2 = 10
77	integer,parameter :: i_o1d = 11
78	integer,parameter :: i_o3 = 12
79	integer,parameter :: i_ar = 13
80	integer,parameter :: i_o = 14
81
82	! Tracer indexes in the GCM:
83	integer,save :: g_co2=0
84	integer,save :: g_co=0
85	integer,save :: g_o=0
86	integer,save :: g_o1d=0
87	integer,save :: g_o2=0
88	integer,save :: g_o3=0
89	integer,save :: g_h=0
90	integer,save :: g_h2=0
91	integer,save :: g_oh=0
92	integer,save :: g_ho2=0
93	integer,save :: g_h2o2=0
94	integer,save :: g_n2=0
95	integer,save :: g_ar=0
96	integer,save :: g_h2o=0
97
98	integer,save :: gcmind(ncompmoldiff) ! array of GCM indexes
99	integer ierr
100
101	logical,save :: firstcall=.true.
102	real abfac(ncompmoldiff)
103	real,save :: dij(ncompmoldiff,ncompmoldiff)
104
105	! Initializations at first call
106	if (firstcall) then
107	call moldiffcoeff(dij)
108	print*,'MOLDIFF EXO'
109
110	! identify the indexes of the tracers we'll need
111	g_co2=igcm_co2
112	if (g_co2.eq.0) then
113	write(,) "moldiff: Error; no CO2 tracer !!!"
114	stop
115	endif
116	g_co=igcm_co
117	if (g_co.eq.0) then
118	write(,) "moldiff: Error; no CO tracer !!!"
119	stop
120	endif
121	g_o=igcm_o
122	if (g_o.eq.0) then
123	write(,) "moldiff: Error; no O tracer !!!"
124	stop
125	endif
126	g_o1d=igcm_o1d
127	if (g_o1d.eq.0) then
128	write(,) "moldiff: Error; no O1D tracer !!!"
129	stop
130	endif
131	g_o2=igcm_o2
132	if (g_o2.eq.0) then
133	write(,) "moldiff: Error; no O2 tracer !!!"
134	stop
135	endif
136	g_o3=igcm_o3
137	if (g_o3.eq.0) then
138	write(,) "moldiff: Error; no O3 tracer !!!"
139	stop
140	endif
141	g_h=igcm_h
142	if (g_h.eq.0) then
143	write(,) "moldiff: Error; no H tracer !!!"
144	stop
145	endif
146	g_h2=igcm_h2
147	if (g_h2.eq.0) then
148	write(,) "moldiff: Error; no H2 tracer !!!"
149	stop
150	endif
151	g_oh=igcm_oh
152	if (g_oh.eq.0) then
153	write(,) "moldiff: Error; no OH tracer !!!"
154	stop
155	endif
156	g_ho2=igcm_ho2
157	if (g_ho2.eq.0) then
158	write(,) "moldiff: Error; no HO2 tracer !!!"
159	stop
160	endif
161	g_h2o2=igcm_h2o2
162	if (g_h2o2.eq.0) then
163	write(,) "moldiff: Error; no H2O2 tracer !!!"
164	stop
165	endif
166	g_n2=igcm_n2
167	if (g_n2.eq.0) then
168	write(,) "moldiff: Error; no N2 tracer !!!"
169	stop
170	endif
171	g_ar=igcm_ar
172	if (g_ar.eq.0) then
173	write(,) "moldiff: Error; no AR tracer !!!"
174	stop
175	endif
176	g_h2o=igcm_h2o_vap
177	if (g_h2o.eq.0) then
178	write(,) "moldiff: Error; no water vapor tracer !!!"
179	stop
180	endif
181
182	cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
183	c fill array to relate local indexes to gcm indexes
184	cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
185
186	gcmind(i_co) = g_co
187	gcmind(i_n2) = g_n2
188	gcmind(i_o2) = g_o2
189	gcmind(i_co2) = g_co2
190	gcmind(i_h2) = g_h2
191	gcmind(i_h) = g_h
192	gcmind(i_oh) = g_oh
193	gcmind(i_ho2) = g_ho2
194	gcmind(i_h2o) = g_h2o
195	gcmind(i_h2o2) = g_h2o2
196	gcmind(i_o1d) = g_o1d
197	gcmind(i_o3) = g_o3
198	gcmind(i_o) = g_o
199	gcmind(i_ar) = g_ar
200
201	firstcall= .false.
202	endif ! of if (firstcall)
203
204
205
206	c
207	cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
208
209	nz=nlayer
210
211	do ig=1,ngrid
212
213	do l=2,nz-1
214	tt(l)=pt(ig,l)+pdt(ig,l)*ptimestep
215	& +pdteuv(ig,l)*ptimestep
216	& +pdtconduc(ig,l)*ptimestep
217
218	do nn=1,ncompmoldiff
219	qq(l,nn)=pq(ig,l,gcmind(nn))+pdq(ig,l,gcmind(nn))*ptimestep
220	qq(l,nn)=max(qq(l,nn),1.e-30)
221	enddo
222	hp(l)=-log(pplay(ig,l+1)/pplay(ig,l-1))
223	dmmeandz(l)=(mmean(ig,l+1)-mmean(ig,l-1))/hp(l)
224	enddo
225
226	tt(1)=pt(ig,1) +pdt(ig,1)*ptimestep
227	& +pdteuv(ig,1)*ptimestep
228	& +pdtconduc(ig,1)*ptimestep
229	tt(nz)=pt(ig,nz)+pdt(ig,nz)*ptimestep
230	& +pdteuv(ig,nz)*ptimestep
231	& +pdtconduc(ig,nz)*ptimestep
232
233	do nn=1,ncompmoldiff
234	qq(1,nn)=pq(ig,1,gcmind(nn))+pdq(ig,1,gcmind(nn))*ptimestep
235	qq(nz,nn)=pq(ig,nz,gcmind(nn))+pdq(ig,nz,gcmind(nn))*ptimestep
236	qq(1,nn)=max(qq(1,nn),1.e-30)
237	qq(nz,nn)=max(qq(nz,nn),1.e-30)
238	enddo
239	hp(1)=-log(pplay(ig,2)/pplay(ig,1))
240	dmmeandz(1)=(-3.mmean(ig,1)+4.mmean(ig,2)
241	& -mmean(ig,3))/(2.*hp(1))
242	hp(nz)=-log(pplay(ig,nz)/pplay(ig,nz-1))
243	dmmeandz(nz)=(3.mmean(ig,nz)-4.mmean(ig,nz-1)
244	& +mmean(ig,nz-2))/(2.*hp(nz))
245	c
246	c Setting-up matrix of alfa coefficients from Dickinson and Ridley 1972
247	c
248	do l=1,nz
249	if(abs(dmmeandz(l)) .lt. 1.e-5) dmmeandz(l)=0.0
250	hh=rnew(ig,l)*tt(l)/g
251	ptfac=(1.e5/pplay(ig,l))(tt(l)/273)*1.75
252	ntot=pplay(ig,l)/(1.381e-23*tt(l)) ! in #/m3
253
254	do nn=1,ncompmoldiff-1 ! rows
255	alfdiag(nn)=0.
256	dcoef1=dij(nn,i_o)*ptfac
257	do n=1,ncompmoldiff-1 ! columns
258	y(nn,n)=0.
259	dcoef=dij(nn,n)*ptfac
260	alf(nn,n)=qq(l,nn)/ntot/1.66e-27
261	& (1./(mmol(gcmind(n))dcoef)-1./(mmol(g_o)*dcoef1))
262	alfdiag(nn)=alfdiag(nn)
263	& +(1./(mmol(gcmind(n))dcoef)-1./(mmol(g_o)dcoef1))
264	& *qq(l,n)
265	enddo
266	dcoef=dij(nn,nn)*ptfac
267	alfdiag(nn)=alfdiag(nn)
268	& -(1./(mmol(gcmind(nn))dcoef)-1./(mmol(g_o)dcoef1))
269	& *qq(l,nn)
270	alf(nn,nn)=-(alfdiag(nn)
271	& +1./(mmol(g_o)*dcoef1))/ntot/1.66e-27
272	y(nn,nn)=1.
273	b(l,nn)=-(dmmeandz(l)/mmean(ig,l)
274	& +mmol(gcmind(nn))/mmean(ig,l)-1.)
275	enddo
276
277	c
278	c Inverting the alfa matrix
279	c
280	call ludcmp_sp(alf,ncompmoldiff-1,ncompmoldiff-1,indx,d,ierr)
281
282	c TEMPORAIRE *****************************
283	if (ierr.ne.0) then
284	write(,)'In moldiff: Problem in LUDCMP_SP with matrix alf'
285	write(,) 'Singular matrix ?'
286	c write(,) 'Matrix alf = ', alf
287	write(,) 'ig, l=',ig, l
288	write(,) 'No molecular diffusion this time !'
289	pdqdiff(1:ngrid,1:nlayer,1:nq)=0
290	return
291	c stop
292	end if
293	c *******************************************
294	do n=1,ncompmoldiff-1
295	call lubksb_sp(alf,ncompmoldiff-1,ncompmoldiff-1,indx,y(1,n))
296	do nn=1,ncompmoldiff-1
297	alfinv(l,nn,n)=y(nn,n)/hh
298	enddo
299	enddo
300	enddo
301
302	c
303	c Calculating coefficients of the system
304	c
305
306	c zlocal(1)=-log(pplay(ig,1)/pplev(ig,1))* Rnew(ig,1)*tt(1)/g
307	zlocal(1)=zzlay(ig,1)
308
309	do l=2,nz-1
310	del1=hp(l)pplay(ig,l)/(gptimestep)
311	del2=(hp(l)/2)*2pplay(ig,l)/(g*ptimestep)
312	do nn=1,ncompmoldiff-1
313	do n=1,ncompmoldiff-1
314	dalfinvdz=(alfinv(l+1,nn,n)-alfinv(l-1,nn,n))/hp(l)
315	aa(l,nn,n)=-dalfinvdz/del1+alfinv(l,nn,n)/del2
316	& +alfinv(l-1,nn,n)*b(l-1,n)/del1
317	bb(l,nn,n)=-2.*alfinv(l,nn,n)/del2
318	cc(l,nn,n)=dalfinvdz/del1+alfinv(l,nn,n)/del2
319	& -alfinv(l+1,nn,n)*b(l+1,n)/del1
320	enddo
321	enddo
322
323	c tmean=tt(l)
324	c if(tt(l).ne.tt(l-1))
325	c & tmean=(tt(l)-tt(l-1))/log(tt(l)/tt(l-1))
326	c zlocal(l)= zlocal(l-1)
327	c & -log(pplay(ig,l)/pplay(ig,l-1))rnew(ig,l)tmean/g
328	zlocal(l)=zzlay(ig,l)
329	enddo
330
331	c zlocal(nz)= zlocal(nz-1)
332	c & -log(pplay(ig,nz)/pplay(ig,nz-1))rnew(ig,nz)tmean/g
333	zlocal(nz)=zzlay(ig,nz)
334
335	ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
336	c Escape velocity from Jeans equation for the flux of H and H2
337	c (Hunten 1973, eq. 5)
338
339	do n=1,ncompmoldiff
340	wi(n)=1.
341	flux(n)=0.
342	abfac(n)=1.
343	enddo
344
345	dens=pplay(ig,nz)/(rnew(ig,nz)*tt(nz))
346	c
347	c For H:
348	c
349	pote=(3398000.+zlocal(nz))/
350	& (1.381e-23tt(nz)/(1.6605e-27mmol(g_h)*g))
351	wi(i_h)=sqrt(2.1.381e-23tt(nz)/(1.6605e-27*mmol(g_h)))
352	& /(2.sqrt(3.1415))(1.+pote)*exp(-pote)
353	flux(i_h)=qq(nz,i_h)dens/(1.6605e-27mmol(g_h))*wi(i_h)
354	flux(i_h)=flux(i_h)*1.6606e-27
355	abfac(i_h)=0.
356	c
357	c For H2:
358	c
359	pote=(3398000.+zlocal(nz))/
360	& (1.381e-23tt(nz)/(1.6605e-27mmol(g_h2)*g))
361	wi(i_h2)=sqrt(2.1.381e-23tt(nz)/(1.6605e-27*mmol(g_h2)))
362	& /(2.sqrt(3.1415))(1.+pote)*exp(-pote)
363	flux(i_h2)=qq(nz,i_h2)dens/(1.6605e-27mmol(g_h2))*wi(i_h2)
364	flux(i_h2)=flux(i_h2)*1.6606e-27
365	abfac(i_h2)=0.
366
367	c ********* TEMPORAIRE : no escape for h and h2
368	c do n=1,ncomptot
369	c wi(n)=1.
370	c flux(n)=0.
371	c abfac(n)=1.
372	c enddo
373	c ********************************************
374
375
376	ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
377
378	c
379	c Setting coefficients for tridiagonal matrix and solving the system
380	c
381
382	do nn=1,ncompmoldiff-1
383	do l=2,nz-1
384	atri(l-1)=aa(l,nn,nn)
385	btri(l-1)=bb(l,nn,nn)+1.
386	ctri(l-1)=cc(l,nn,nn)
387	rtri(l-1)=qq(l,nn)
388	do n=1,ncompmoldiff-1
389	rtri(l-1)=rtri(l-1)-(aa(l,nn,n)*qq(l-1,n)
390	& +bb(l,nn,n)*qq(l,n)
391	& +cc(l,nn,n)*qq(l+1,n))
392	enddo
393	rtri(l-1)=rtri(l-1)+(aa(l,nn,nn)*qq(l-1,nn)
394	& +bb(l,nn,nn)*qq(l,nn)
395	& +cc(l,nn,nn)*qq(l+1,nn))
396	enddo
397
398	c
399	c Boundary conditions:
400	c Escape flux for H and H2 at top
401	c Diffusive equilibrium for the other species at top
402	c Perfect mixing for all at bottom
403	c
404
405	rtri(nz-2)=rtri(nz-2)
406	& -ctri(nz-2)flux(nn)mmol(gcmind(nn))/(dens*wi(nn))
407
408	atri(nz-2)=atri(nz-2)
409	& -abfac(nn)ctri(nz-2)/(3.-2.hp(nz)*b(nz,nn))
410	btri(nz-2)=btri(nz-2)
411	& +abfac(nn)4.ctri(nz-2)/(3.-2.hp(nz)b(nz,nn))
412
413	c rtri(1)=rtri(1)-atri(1)*qq(1,nn)
414	btri(1)=btri(1)+atri(1)
415
416	call tridag_sp(atri,btri,ctri,rtri,qtri,nz-2)
417
418	do l=2,nz-1
419	c qnew(l,nn)=qtri(l-1)
420	qnew(l,nn)=max(qtri(l-1),1.e-30)
421	enddo
422
423	qnew(nz,nn)=flux(nn)mmol(gcmind(nn))/(denswi(nn))
424	& +abfac(nn)(4.qnew(nz-1,nn)-qnew(nz-2,nn))
425	& /(3.-2.hp(nz)b(nz,nn))
426	c qnew(1,nn)=qq(1,nn)
427	qnew(1,nn)=qnew(2,nn)
428
429	qnew(nz,nn)=max(qnew(nz,nn),1.e-30)
430	qnew(1,nn)=max(qnew(1,nn),1.e-30)
431
432	enddo ! loop on species
433
434	DO l=1,nz
435	if(zlocal(l).gt.65000.)then
436	pdqdiff(ig,l,g_o)=0.
437	do n=1,ncompmoldiff-1
438	pdqdiff(ig,l,gcmind(n))=(qnew(l,n)-qq(l,n))
439	pdqdiff(ig,l,g_o)=pdqdiff(ig,l,g_o)-(qnew(l,n)-qq(l,n))
440	pdqdiff(ig,l,gcmind(n))=pdqdiff(ig,l,gcmind(n))/ptimestep
441	enddo
442	pdqdiff(ig,l,g_o)=pdqdiff(ig,l,g_o)/ptimestep
443	endif
444	ENDDO
445
446	c do l=2,nz
447	c do n=1,ncomptot-1
448	c hco2(n)=qnew(l,n)pplay(ig,l)/(rnew(ig,l)tt(l)) /
449	c & (qnew(l-1,n)pplay(ig,l-1)/(rnew(ig,l-1)tt(l-1)))
450	c hco2(n)=-(zlocal(l)-zlocal(l-1))/log(hco2(n))/1000.
451	c enddo
452	c write(225,*),l,pt(1,l),(hco2(n),n=1,6)
453	c write(226,*),l,pt(1,l),(hco2(n),n=7,12)
454	c enddo
455
456	enddo ! ig loop
457
458	return
459	end
460
461	c ********************************************************************
462	c ********************************************************************
463	c ********************************************************************
464
465	subroutine tridag_sp(a,b,c,r,u,n)
466	c parameter (nmax=100)
467	c dimension gam(nmax),a(n),b(n),c(n),r(n),u(n)
468	real gam(n),a(n),b(n),c(n),r(n),u(n)
469	if(b(1).eq.0.)then
470	stop 'tridag_sp: error: b(1)=0 !!! '
471	endif
472	bet=b(1)
473	u(1)=r(1)/bet
474	do 11 j=2,n
475	gam(j)=c(j-1)/bet
476	bet=b(j)-a(j)*gam(j)
477	if(bet.eq.0.) then
478	stop 'tridag_sp: error: bet=0 !!! '
479	endif
480	u(j)=(r(j)-a(j)*u(j-1))/bet
481	11 continue
482	do 12 j=n-1,1,-1
483	u(j)=u(j)-gam(j+1)*u(j+1)
484	12 continue
485	return
486	end
487
488	c ********************************************************************
489	c ********************************************************************
490	c ********************************************************************
491
492	SUBROUTINE LUBKSB_SP(A,N,NP,INDX,B)
493
494	implicit none
495
496	integer i,j,n,np,ii,ll
497	real sum
498	real a(np,np),indx(np),b(np)
499
500	c DIMENSION A(NP,NP),INDX(N),B(N)
501	II=0
502	DO 12 I=1,N
503	LL=INDX(I)
504	SUM=B(LL)
505	B(LL)=B(I)
506	IF (II.NE.0)THEN
507	DO 11 J=II,I-1
508	SUM=SUM-A(I,J)*B(J)
509	11 CONTINUE
510	ELSE IF (SUM.NE.0.) THEN
511	II=I
512	ENDIF
513	B(I)=SUM
514	12 CONTINUE
515	DO 14 I=N,1,-1
516	SUM=B(I)
517	IF(I.LT.N)THEN
518	DO 13 J=I+1,N
519	SUM=SUM-A(I,J)*B(J)
520	13 CONTINUE
521	ENDIF
522	B(I)=SUM/A(I,I)
523	14 CONTINUE
524	RETURN
525	END
526
527	c ********************************************************************
528	c ********************************************************************
529	c ********************************************************************
530
531	SUBROUTINE LUDCMP_SP(A,N,NP,INDX,D,ierr)
532
533	implicit none
534
535	integer n,np,nmax,i,j,k,imax
536	real d,tiny,aamax
537	real a(np,np),indx(np)
538	integer ierr ! error =0 if OK, =1 if problem
539
540	PARAMETER (NMAX=100,TINY=1.0E-20)
541	c DIMENSION A(NP,NP),INDX(N),VV(NMAX)
542	real sum,vv(nmax),dum
543
544	D=1.
545	DO 12 I=1,N
546	AAMAX=0.
547	DO 11 J=1,N
548	IF (ABS(A(I,J)).GT.AAMAX) AAMAX=ABS(A(I,J))
549	11 CONTINUE
550	IF (AAMAX.EQ.0.) then
551	write(,) 'In moldiff: Problem in LUDCMP_SP with matrix A'
552	write(,) 'Singular matrix ?'
553	c write(,) 'Matrix A = ', A
554	c TO DEBUG :
555	ierr =1
556	return
557	c stop
558	END IF
559
560	VV(I)=1./AAMAX
561	12 CONTINUE
562	DO 19 J=1,N
563	IF (J.GT.1) THEN
564	DO 14 I=1,J-1
565	SUM=A(I,J)
566	IF (I.GT.1)THEN
567	DO 13 K=1,I-1
568	SUM=SUM-A(I,K)*A(K,J)
569	13 CONTINUE
570	A(I,J)=SUM
571	ENDIF
572	14 CONTINUE
573	ENDIF
574	AAMAX=0.
575	DO 16 I=J,N
576	SUM=A(I,J)
577	IF (J.GT.1)THEN
578	DO 15 K=1,J-1
579	SUM=SUM-A(I,K)*A(K,J)
580	15 CONTINUE
581	A(I,J)=SUM
582	ENDIF
583	DUM=VV(I)*ABS(SUM)
584	IF (DUM.GE.AAMAX) THEN
585	IMAX=I
586	AAMAX=DUM
587	ENDIF
588	16 CONTINUE
589	IF (J.NE.IMAX)THEN
590	DO 17 K=1,N
591	DUM=A(IMAX,K)
592	A(IMAX,K)=A(J,K)
593	A(J,K)=DUM
594	17 CONTINUE
595	D=-D
596	VV(IMAX)=VV(J)
597	ENDIF
598	INDX(J)=IMAX
599	IF(J.NE.N)THEN
600	IF(A(J,J).EQ.0.)A(J,J)=TINY
601	DUM=1./A(J,J)
602	DO 18 I=J+1,N
603	A(I,J)=A(I,J)*DUM
604	18 CONTINUE
605	ENDIF
606	19 CONTINUE
607	IF(A(N,N).EQ.0.)A(N,N)=TINY
608	ierr =0
609	RETURN
610	END
611

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

Download in other formats:

Original Format