Context Navigation

moldiff.F @ 1119

Last change on this file since 1119 was 1047, checked in by emillour, 11 years ago

Mars GCM:

IMPORTANT CHANGE: Removed all reference/use of ngridmx (dimphys.h) in routines (necessary prerequisite to using parallel dynamics); in most cases this just means adding 'ngrid' as routine argument, and making local saved variables allocatable (and allocated at first call). In the process, had to convert many *.h files to equivalent modules: yomaer.h => yomaer_h.F90 , surfdat.h => surfdat_h.F90 , comsaison.h => comsaison_h.F90 , yomlw.h => yomlw_h.F90 , comdiurn.h => comdiurn_h.F90 , dimradmars.h => dimradmars_mod.F90 , comgeomfi.h => comgeomfi_h.F90, comsoil.h => comsoil_h.F90 , slope.h => slope_mod.F90
Also updated EOF routines, everything is now in eofdump_mod.F90
Removed unused routine lectfux.F (in dyn3d)

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

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

Download in other formats:

Original Format