source: trunk/LMDZ.MARS/libf/aeronomars/moldiff.F @ 2599

Last change on this file since 2599 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
13c
14c Input/Output
15c
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)
30c
31c Local
32c
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
62cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
63c     tracer numbering in the molecular diffusion
64cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
65c  Atomic oxygen must always be the LAST species of the list as
66c 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
182cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
183c   fill array to relate local indexes to gcm indexes
184cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
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
206c
207cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
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))
245c
246c Setting-up matrix of alfa coefficients from Dickinson and Ridley 1972
247c
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
277c
278c Inverting the alfa matrix
279c
280        call ludcmp_sp(alf,ncompmoldiff-1,ncompmoldiff-1,indx,d,ierr)
281
282c       TEMPORAIRE *****************************
283        if (ierr.ne.0) then
284            write(*,*)'In moldiff: Problem in LUDCMP_SP with matrix alf'
285            write(*,*) 'Singular matrix ?'
286c           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
291c           stop
292        end if
293c       *******************************************
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
302c
303c Calculating coefficients of the system
304c
305
306c      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)/(g*ptimestep)
311        del2=(hp(l)/2)**2*pplay(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
323c        tmean=tt(l)
324c        if(tt(l).ne.tt(l-1))
325c     &        tmean=(tt(l)-tt(l-1))/log(tt(l)/tt(l-1))
326c        zlocal(l)= zlocal(l-1)
327c     &         -log(pplay(ig,l)/pplay(ig,l-1))*rnew(ig,l)*tmean/g
328      zlocal(l)=zzlay(ig,l)
329      enddo
330
331c      zlocal(nz)= zlocal(nz-1)
332c     &         -log(pplay(ig,nz)/pplay(ig,nz-1))*rnew(ig,nz)*tmean/g
333      zlocal(nz)=zzlay(ig,nz)
334       
335ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
336c Escape velocity from Jeans equation for the flux of H and H2
337c (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))
346c
347c For H:
348c
349       pote=(3398000.+zlocal(nz))/
350     &         (1.381e-23*tt(nz)/(1.6605e-27*mmol(g_h)*g))
351       wi(i_h)=sqrt(2.*1.381e-23*tt(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-27*mmol(g_h))*wi(i_h)
354       flux(i_h)=flux(i_h)*1.6606e-27
355       abfac(i_h)=0.
356c
357c For H2:
358c
359       pote=(3398000.+zlocal(nz))/
360     &         (1.381e-23*tt(nz)/(1.6605e-27*mmol(g_h2)*g))
361       wi(i_h2)=sqrt(2.*1.381e-23*tt(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-27*mmol(g_h2))*wi(i_h2)
364       flux(i_h2)=flux(i_h2)*1.6606e-27
365       abfac(i_h2)=0.
366
367c ********* TEMPORAIRE : no escape for h and h2
368c     do n=1,ncomptot
369c       wi(n)=1.
370c       flux(n)=0.
371c       abfac(n)=1.
372c     enddo
373c ********************************************
374
375
376ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
377
378c
379c Setting coefficients for tridiagonal matrix and solving the system
380c
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
398c
399c Boundary conditions:
400c       Escape flux for H and H2 at top
401c       Diffusive equilibrium for the other species at top
402c       Perfect mixing for all at bottom
403c
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
413c        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
419c          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))/(dens*wi(nn))
424     &               +abfac(nn)*(4.*qnew(nz-1,nn)-qnew(nz-2,nn))
425     &                /(3.-2.*hp(nz)*b(nz,nn))
426c        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
446c      do l=2,nz
447c        do n=1,ncomptot-1
448c          hco2(n)=qnew(l,n)*pplay(ig,l)/(rnew(ig,l)*tt(l)) /
449c     &      (qnew(l-1,n)*pplay(ig,l-1)/(rnew(ig,l-1)*tt(l-1)))
450c          hco2(n)=-(zlocal(l)-zlocal(l-1))/log(hco2(n))/1000.
451c        enddo
452c        write(225,*),l,pt(1,l),(hco2(n),n=1,6)
453c        write(226,*),l,pt(1,l),(hco2(n),n=7,12)
454c      enddo
455
456      enddo             ! ig loop
457
458      return
459      end
460
461c    ********************************************************************
462c    ********************************************************************
463c    ********************************************************************
464 
465      subroutine tridag_sp(a,b,c,r,u,n)
466c      parameter (nmax=100)
467c      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
48111    continue
482      do 12 j=n-1,1,-1
483        u(j)=u(j)-gam(j+1)*u(j+1)
48412    continue
485      return
486      end
487
488c    ********************************************************************
489c    ********************************************************************
490c    ********************************************************************
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
500c      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)
50911        CONTINUE
510        ELSE IF (SUM.NE.0.) THEN
511          II=I
512        ENDIF
513        B(I)=SUM
51412    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)
52013        CONTINUE
521        ENDIF
522        B(I)=SUM/A(I,I)
52314    CONTINUE
524      RETURN
525      END
526
527c    ********************************************************************
528c    ********************************************************************
529c    ********************************************************************
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)
541c      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))
54911      CONTINUE
550        IF (AAMAX.EQ.0.) then
551            write(*,*) 'In moldiff: Problem in LUDCMP_SP with matrix A'
552            write(*,*) 'Singular matrix ?'
553c           write(*,*) 'Matrix A = ', A
554c           TO DEBUG :
555            ierr =1
556            return
557c           stop
558        END IF
559
560        VV(I)=1./AAMAX
56112    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)
56913            CONTINUE
570              A(I,J)=SUM
571            ENDIF
57214        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)
58015          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
58816      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
59417        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
60418        CONTINUE
605        ENDIF
60619    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.