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

Last change on this file since 3567 was 3466, checked in by emillour, 2 months ago

Mars PCM:
More tidying in aeronomars:

  • remove unused "inv.F" and remove "dtridgl.F" which is not used here and is a duplicate of the "dtridgl" routine in phymars/swr_toon.F
  • turn aeronomars routines to modules, for those which aren't in modules yet.

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