source: trunk/LMDZ.MARS/libf/aeronomars/moldiff_MPF.F90

Last change on this file was 3443, checked in by emillour, 3 months ago

Mars PCM:
Add new molecular diffusion routine moldiff_MPF (modified pass flow scheme)
But for intermetiate testing the legacy moldiff_red routine remains
the default. The moldiff_scheme flag (1: legacy, 2: MPF) can be used
to swich from one scheme to the other.
JYC+EM
File size: 42.4 KB
RevLine 
[3443]1MODULE moldiff_MPF_mod
2
3  real*8,parameter :: Pdiff=15. ! pressure (Pa) below which diffusion is computed
4  real*8,parameter :: tdiffmin=5d0
5  real*8,parameter :: dzres=2d0 ! grid resolution (km) for diffusion
6
7CONTAINS
8
9subroutine moldiff_MPF(ngrid,nlayer,nq,pplay,pplev,pt,pdt,pq,pdq,&
10                       ptimestep,zzlay,pdteuv,pdtconduc,pdqdiff,&
11                       PhiEscH,PhiEscH2,PhiEscD)
12
13use tracer_mod, only: noms, mmol
14use geometry_mod, only: cell_area
15use planetwide_mod, only: planetwide_sumval
16USE mod_phys_lmdz_para, only: is_master,bcast
17use moldiffcoeff_red_mod, only: moldiffcoeff_red
18
19
20implicit none
21
22! July 2014 JYC ADD BALISTIC Transport coupling to compute wup for H and H2
23! June 2023 JYC New method based on the modified pass flow (Parshev et al. 1987)
24
25
26!
27! Input/Output
28!
29      integer,intent(in) :: ngrid ! number of atmospheric columns
30      integer,intent(in) :: nlayer ! number of atmospheric layers
31      integer,intent(in) :: nq ! number of advected tracers
32      real ptimestep
33      real pplay(ngrid,nlayer)
34      real zzlay(ngrid,nlayer)
35      real pplev(ngrid,nlayer+1)
36      real pq(ngrid,nlayer,nq)
37      real pdq(ngrid,nlayer,nq)
38      real pt(ngrid,nlayer)
39      real pdt(ngrid,nlayer)
40      real pdteuv(ngrid,nlayer)
41      real pdtconduc(ngrid,nlayer)
42      real pdqdiff(ngrid,nlayer,nq)
43      real*8 PhiEscH,PhiEscH2,PhiEscD
44!
45! Local
46!
47     
48
49!      real hco2(ncompdiff),ho
50
51      integer,dimension(nq) :: indic_diff
52      integer ig,iq,nz,l,k,n,nn,p,ij0
53      integer istep,il,gcn,ntime,nlraf
54      real*8 masse
55      real*8 masseU,kBolt,RgazP,Rmars,Grav,Mmars
56      real*8 rho0,D0,T0,H0,time0,dZ,time,dZraf,tdiff,Zmin,Zmax,K0,Pk
57      real*8 FacEsc,invsgmu,Ueff,alphaTnn
58      real*8 PhiauxH(ngrid),PhiauxH2(ngrid),PhiauxD(ngrid)
59      real*8 hp(nlayer)
60      real*8 pp(nlayer)
61      real*8 pint(nlayer)
62      real*8 tt(nlayer),tnew(nlayer),tint(nlayer)
63      real*8 zz(nlayer)
64      real*8,dimension(:,:),allocatable,save :: qq,qnew,qint,FacMass
65      real*8,dimension(:,:),allocatable,save :: rhoK,rhokinit
66      real*8 rhoT(nlayer)
67      real*8 dmmeandz(nlayer)
68      real*8 massemoy(nlayer)
69      real*8,dimension(:),allocatable :: Praf,Traf,Rraf,Mraf,Nraf,Draf,Kraf,Hraf,Wraf
70      real*8,dimension(:),allocatable :: Zraf,Tdiffraf
71      real*8,dimension(:),allocatable :: Prafold,Mrafold
72      real*8,dimension(:,:),allocatable :: Qraf,Rrafk,Nrafk
73      real*8,dimension(:,:),allocatable :: Rrafkold
74      real*8,dimension(:,:),allocatable :: Drafmol,Hrafmol,Wrafmol,Tdiffrafmol
75      real*8,dimension(:),allocatable :: Atri,Btri,Ctri,Dtri,Xtri,Tad,Dad,Zad,rhoad,Kad
76      real*8,dimension(:),allocatable :: alpha,beta,delta,ksi,eps,zeta,prod,loss
77      real*8,dimension(:),allocatable,save :: wi,Wad,Uthermal,Lambdaexo,Hspecie,alphaT
78      real*8,dimension(:),allocatable,save :: Mtot1,Mtot2,Mraf1,Mraf2
79      integer,parameter :: ListeDiffNb=16
80      character(len=20),dimension(ListeDiffNb) :: ListeDiff
81      real*8,parameter :: pi=3.141592653
82      real*8,parameter :: g=3.72d0
83!ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
84!     tracer numbering in the molecular diffusion
85!ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
86!  We need the index of escaping species
87
88      integer,save :: i_h2 
89      integer,save :: i_h
90      integer,save :: i_d
91      integer,save :: i_hd
92! vertical index limit for the molecular diffusion
93      integer,save :: il0 
94
95!$OMP THREADPRIVATE(qq,qnew,qint,FacMass,rhoK,rhokinit,wi,Wad,Uthermal,Lambdaexo,Hspecie,Mtot1,Mtot2,Mraf1,Mraf2,i_h2,i_h,i_d,il0)
96
97! Tracer indexes in the GCM:
98!      integer,save :: g_co2=0
99!      integer,save :: g_co=0
100!      integer,save :: g_o=0
101!      integer,save :: g_o1d=0
102!      integer,save :: g_o2=0
103!      integer,save :: g_o3=0
104!      integer,save :: g_h=0
105!      integer,save :: g_h2=0
106!      integer,save :: g_oh=0
107!      integer,save :: g_ho2=0
108!      integer,save :: g_h2o2=0
109!      integer,save :: g_n2=0
110!      integer,save :: g_ar=0
111!      integer,save :: g_h2o=0
112!      integer,save :: g_n=0
113!      integer,save :: g_no=0
114!      integer,save :: g_no2=0
115!      integer,save :: g_n2d=0
116!      integer,save :: g_oplus=0
117!      integer,save :: g_hplus=0
118
119      integer,save :: ncompdiff
120      integer,dimension(:),allocatable,save :: gcmind ! array of GCM indexes
121      integer ierr,compteur
122
123      logical,save :: firstcall=.true.
124!      real abfac(ncompdiff)
125      real,dimension(:,:),allocatable,save :: dij
126      real,save :: step
127
128!$OMP THREADPRIVATE(ncompdiff,gcmind,firstcall,dij,step)     
129
130! Initializations at first call
131      if (firstcall) then
132
133! Liste des especes qui sont diffuses si elles sont presentes
134!       ListeDiff=['co2','o','n2','ar','co','h2','h','d2','d','hd','o2','oh','ho2','h2o_vap','h2o2','o1d','n','no','no2']
135        ListeDiff(1)='co2'
136        ListeDiff(2)='o'
137        ListeDiff(3)='n2'
138        ListeDiff(4)='ar'
139        ListeDiff(5)='co'
140        ListeDiff(6)='h2'
141        ListeDiff(7)='h'
142        ListeDiff(8)='d2'
143        ListeDiff(9)='hd'
144        ListeDiff(10)='d'
145        ListeDiff(11)='o2'
146        ListeDiff(12)='h2o_vap'
147        ListeDiff(13)='o3'
148        ListeDiff(14)='n'
149        ListeDiff(15)='he'
150        ListeDiff(16)='hdo_vap'
151        i_h=1000
152        i_h2=1000
153        i_d=1000
154        i_hd=1000
155! On regarde quelles especes diffusables sont presentes
156
157        ncompdiff=0
158        indic_diff(1:nq)=0
159       
160        do nn=1,nq
161        do n=1,ListeDiffNb
162        if (ListeDiff(n) .eq. noms(nn)) then
163        indic_diff(nn)=1
164!       if (noms(nn) .eq. 'h') i_h=n
165!       if (noms(nn) .eq. 'h2') i_h2=n
166        endif
167
168        enddo
169        if (indic_diff(nn) .eq. 1) then
170        print*,'specie ', noms(nn), 'diffused in moldiff_MPF'
171        ncompdiff=ncompdiff+1
172        endif
173        enddo
174        print*,'number of diffused species:',ncompdiff
175        allocate(gcmind(ncompdiff))
176
177! Store gcm indexes in gcmind
178        n=0
179        do nn=1,nq
180        if (indic_diff(nn) .eq. 1) then
181        n=n+1
182        gcmind(n)=nn
183        if (noms(nn) .eq. 'h') i_h=n
184        if (noms(nn) .eq. 'h2') i_h2=n
185        if (noms(nn) .eq. 'd') i_d=n
186        if (noms(nn) .eq. 'hd') i_hd=n
187        endif
188        enddo
189
190!       print*,'gcmind',gcmind,i_h,i_h2
191
192! find vertical index above which diffusion is computed
193        if(is_master) then
194        do l=1,nlayer
195        if (pplay(1,l) .gt. Pdiff) then
196        il0=l
197        endif
198        enddo
199        il0=il0+1
200        endif ! (is_master)
201        CALL bcast(il0)
202        print*,'vertical index for diffusion',il0,pplay(1,il0)
203
204        allocate(dij(ncompdiff,ncompdiff))
205
206        call moldiffcoeff_red(dij,indic_diff,gcmind,ncompdiff)
207        print*,'MOLDIFF  EXO'
208
209! allocatation des tableaux dependants du nombre d especes diffusees
210        allocate(qq(nlayer,ncompdiff))
211        allocate(qnew(nlayer,ncompdiff))
212        allocate(qint(nlayer,ncompdiff))
213        allocate(FacMass(nlayer,ncompdiff))
214        allocate(rhok(nlayer,ncompdiff))
215        allocate(rhokinit(nlayer,ncompdiff))
216
217        allocate(wi(ncompdiff))
218        allocate(wad(ncompdiff))
219        allocate(alphaT(ncompdiff))
220        allocate(uthermal(ncompdiff))
221        allocate(lambdaexo(ncompdiff))
222        allocate(Hspecie(ncompdiff))
223        allocate(Mtot1(ncompdiff))
224        allocate(Mtot2(ncompdiff))
225        allocate(Mraf1(ncompdiff))
226        allocate(Mraf2(ncompdiff))
227
228        firstcall= .false.
229        step=1
230      endif ! of if (firstcall)
231
232!
233!ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
234
235      masseU=1.660538782d-27
236      kbolt=1.3806504d-23
237      RgazP=8.314472   
238!      PI =3.141592653
239!      g=3.72d0
240      Rmars=3390000d0 ! Used to compute escape parameter no need to be more accurate
241        Grav=6.67d-11
242        Mmars=6.4d23
243        ij0=6000 ! For test
244        invsgmu=1d0/g/masseU
245        K0=1.2d11 ! coefficient for eddy diffusion (using n in m-3 and K in m2/s so diff from Krasno2002)
246!        K0=2.4d11
247!        K0=6.0d10
248!        K0=2.4d10
249!        K0=6.0d11
250!        K0=1.2d10
251!        K0=1.2d12
252!        Pk=0.1 ! If P > Pk K = 0 (Pk in Pa)
253
254! Compute the wup(ig) for H and H2 using the balistic code from R Yelle
255
256        PhiEscH=0D0
257        PhiEscH2=0D0
258        PhiEscD=0D0
259
260      do ig=1,ngrid
261        pp=dble(pplay(ig,:))
262
263! Update the temperature
264
265!       CALL TMNEW(pt(ig,:),pdt(ig,:),pdtconduc(ig,:),pdteuv(ig,:)      &
266!     &  ,tt,ptimestep,nlayer,ig)
267        do l=1,nlayer
268          tt(l)=pt(ig,l)*1D0+(pdt(ig,l)*dble(ptimestep)+                &
269                              pdtconduc(ig,l)*dble(ptimestep)+          &
270                              pdteuv(ig,l)*dble(ptimestep))
271          ! to cach Nans...
272          if (tt(l).ne.tt(l)) then
273            print*,'Err TMNEW',ig,l,tt(l),pt(ig,l),                     &
274              pdt(ig,l),pdtconduc(ig,l),pdteuv(ig,l),dble(ptimestep)
275          endif
276        enddo ! of do l=1,nlayer
277
278! Update the mass mixing ratios modified by other processes
279
280!       CALL QMNEW(pq(ig,:,:),pdq(ig,:,:),qq,ptimestep,nlayer,        &
281!     &  ncompdiff,gcmind,ig)
282        do iq=1,ncompdiff
283         do l=1,nlayer
284          qq(l,iq)=pq(ig,l,gcmind(iq))*1D0+(                            &
285                           pdq(ig,l,gcmind(iq))*dble(ptimestep))
286          qq(l,iq)=max(qq(l,iq),1d-30)
287         enddo ! of do l=1,nlayer
288        enddo ! of do iq=1,ncompdiff
289
290! Compute the Pressure scale height
291
292        CALL HSCALE(pp,hp,nlayer)
293
294! Compute the atmospheric mass (in Dalton)
295
296        CALL MMOY(massemoy,mmol,qq,gcmind,nlayer,ncompdiff)
297
298! Compute the vertical gradient of atmospheric mass
299
300        CALL DMMOY(massemoy,hp,dmmeandz,nlayer)
301
302! Compute the altitude of each layer
303
304        CALL ZVERT(pp,tt,massemoy,zz,nlayer,ig)
305
306! Compute the total mass density (kg/m3)
307       
308        CALL RHOTOT(pp,tt,massemoy,qq,RHOT,RHOK,nlayer,ncompdiff)
309        RHOKINIT=RHOK
310
311! Compute total mass of each specie before new grid
312! For conservation tests
313! The conservation is not always fulfilled especially
314! for species very far from diffusion equilibrium "photochemical species"
315! e.g. OH, O(1D)
316
317        Mtot1(1:ncompdiff)=0d0
318
319        do l=il0,nlayer
320        do nn=1,ncompdiff
321        Mtot1(nn)=Mtot1(nn)+1d0/g*qq(l,nn)*                             &
322     &  (dble(pplev(ig,l))-dble(pplev(ig,l+1)))
323        enddo
324        enddo
325
326        Zmin=zz(il0)
327        Zmax=zz(nlayer)
328
329
330! If Zmax > 4000 km there is a problem / stop
331
332        if (Zmax .gt. 4000000.) then
333        Print*,'Zmax too high',ig,zmax,zmin
334        do l=1,nlayer
335        print*,'old',zz(l),pt(ig,l),pdteuv(ig,l),pdq(ig,l,:)
336        print*,'l',l,rhot(l),tt(l),pp(l),massemoy(l),qq(l,:)
337        enddo
338        stop
339        endif
340
341! The number of diffusion layers is variable
342! and fixed by the resolution (dzres) specified in diffusion.h
343! I fix a minimum  number of layers 40
344
345        nlraf=int((Zmax-Zmin)/1000./dzres)+1
346        if (nlraf .le. 40)  nlraf=40
347
348!       if (nlraf .ge. 200) print*,ig,nlraf,Zmin,Zmax
349       
350         ! allocate arrays:
351         allocate(Praf(nlraf),Traf(nlraf),Rraf(nlraf),Mraf(nlraf))
352         allocate(Nraf(nlraf),Draf(nlraf),Hraf(nlraf),Wraf(nlraf))
353         allocate(Zraf(nlraf),Tdiffraf(nlraf),Kraf(nlraf))
354         allocate(Prafold(nlraf),Mrafold(nlraf))
355         allocate(Qraf(nlraf,ncompdiff),Rrafk(nlraf,ncompdiff),Nrafk(nlraf,ncompdiff))
356         allocate(Rrafkold(nlraf,ncompdiff))
357         allocate(Drafmol(nlraf,ncompdiff),Hrafmol(nlraf,ncompdiff))
358         allocate(Wrafmol(nlraf,ncompdiff),Tdiffrafmol(nlraf,ncompdiff))
359         allocate(Atri(nlraf),Btri(nlraf),Ctri(nlraf),Dtri(nlraf),Xtri(nlraf))
360         allocate(Tad(nlraf),Dad(nlraf),Zad(nlraf),rhoad(nlraf),Kad(nlraf))
361         allocate(alpha(nlraf),beta(nlraf),delta(nlraf),eps(nlraf),ksi(nlraf),zeta(nlraf))
362         allocate(prod(nlraf),loss(nlraf))
363
364! before beginning, I use a better vertical resolution above il0,
365! altitude grid reinterpolation
366! The diffusion is solved on an altitude grid with constant step dz
367
368        CALL UPPER_RESOL(pp,tt,zz,massemoy,RHOT,RHOK,                   &
369     &  qq,mmol,gcmind,Praf,Traf,Qraf,Mraf,Zraf,                        &
370     &  Nraf,Nrafk,Rraf,Rrafk,il0,nlraf,ncompdiff,nlayer,ig)
371
372        Prafold=Praf
373        Rrafkold=Rrafk
374        Mrafold=Mraf
375
376! Eddy mixing profile from Krasnopolsky 2002 Kraf is in cm2/s
377     do l=1,nlraf
378      Kraf(l)=K0*sqrt(Traf(l)*Traf(nlraf)*kbolt/Praf(l))
379!      if (Praf(l) .ge. Pk) Kraf(l)=0D0
380!      print*,l,Praf(l),Traf(l),Kraf(l),Traf(nlraf)
381     enddo
382   
383! We reddo interpolation of the gcm grid to estimate mass loss due to interpolation processes.
384
385        CALL GCMGRID_P(Zraf,Praf,Qraf,Traf,Nrafk,Rrafk,qq,qint,tt,tint  &
386     &    ,pp,mmol,gcmind,nlraf,ncompdiff,nlayer,ig)
387
388! We compute the mass correction factor of each specie at each pressure level
389
390        CALL CORRMASS(qq,qint,FacMass,nlayer,ncompdiff)
391
392! Altitude step
393
394        Dzraf=Zraf(2)-Zraf(1)
395
396! Total mass computed on the altitude grid
397
398        Mraf1(1:ncompdiff)=0d0
399        do nn=1,ncompdiff
400        do l=1,nlraf
401        Mraf1(nn)=Mraf1(nn)+Rrafk(l,nn)*Dzraf
402        enddo
403        enddo
404
405! Reupdate values for mass conservation
406
407!       do l=1,nlraf
408!       print*,'test',l,Nraf(l),Praf(l)
409!       do nn=1,ncompdiff
410!       Rrafk(l,nn)=RrafK(l,nn)*Mtot1(nn)/Mraf1(nn)
411!       enddo
412!       Rraf(l)=sum(Rrafk(l,:))
413!       do nn=1,ncompdiff
414!       Qraf(l,nn)=Rrafk(l,nn)/Rraf(l)
415!       Nrafk(l,nn)=Rrafk(l,nn)/dble(mmol(gcmind(nn)))/masseU
416!       enddo
417!       Mraf(l)=1d0/sum(Qraf(l,:)/dble(mmol(gcmind(:))))
418!       Nraf(l)=Rraf(l)/Mraf(l)/masseU
419!       Praf(l)=kbolt*Traf(l)*Nraf(l)
420!       print*,'test',l,Nraf(l),Praf(l)
421!       enddo
422
423!       do l=1,nlayer
424!       print*,'l',l,zz(l),pp(l),tt(l),sum(qq(l,:)),massemoy(l)
425!       enddo
426
427! The diffusion is computed above il0 computed from Pdiff in diffusion.h
428! No change below il0
429       
430        do l=1,nlayer
431         qnew(l,:)=qq(l,:) ! No effet below il0
432       enddo
433
434! all species treated independently
435
436! Upper boundary velocity
437! Jeans escape for H and H2
438! Comparison with and without escape still to be done
439! No escape for other species
440
441
442        do nn=1,ncompdiff
443        Uthermal(nn)=sqrt(2d0*kbolt*Traf(nlraf)/masseU/                 &
444     &  dble(mmol(gcmind(nn))))
445        Lambdaexo(nn)=masseU*dble(mmol(gcmind(nn)))*Grav*Mmars/         &
446     &  (Rmars+Zraf(nlraf))/kbolt/Traf(nlraf)
447        wi(nn)=0D0
448        alphaT(nn)=0D0
449        if (nn .eq. i_h .or. nn .eq. i_h2 .or. nn .eq. i_d) then
450        wi(nn)=Uthermal(nn)/2d0/sqrt(PI)*exp(-Lambdaexo(nn))*           &
451     &  (Lambdaexo(nn)+1d0)
452        endif
453        if (nn .eq. i_h .or. nn .eq. i_h2 .or. nn .eq. i_d .or. nn .eq. i_hd) then
454        alphaT(nn)=-0.25D0
455        endif
456        enddo
457
458!       print*,'wi',wi(i_h),wi(i_h2),wi(i_d),Uthermal,Lambdaexo,mmol(gcmind(:))
459!       print*,'wi',wi
460
461! Compute time step for diffusion
462
463! Loop on species
464
465        T0=Traf(nlraf)
466        rho0=1d0
467
468        do nn=1,ncompdiff
469        masse=dble(mmol(gcmind(nn)))
470! DIffusion coefficient
471        CALL DCOEFF(nn,dij,Praf,Traf,Nraf,Nrafk,Draf,nlraf,ncompdiff)
472        Drafmol(:,nn)=Draf(:)
473! Scale height of the density of the specie
474        CALL HSCALEREAL(nn,Nrafk,Dzraf,Hraf,nlraf,ncompdiff)
475        Hrafmol(:,nn)=Hraf(:)
476!       Hspecie(nn)=kbolt*T0/masse*invsgmu
477! Computation of the diffusion vertical velocity of the specie
478        CALL VELVERT(nn,Traf,Hraf,Draf,Dzraf,masse,Wraf,nlraf)
479        Wrafmol(:,nn)=Wraf(:)
480! Computation of the diffusion time
481        CALL TIMEDIFF(nn,Hraf,Wraf,Tdiffraf,nlraf)
482        Tdiffrafmol(:,nn)=Tdiffraf
483        enddo
484! We use a lower time step
485        Tdiff=minval(Tdiffrafmol)
486        Tdiff=minval(Tdiffrafmol(nlraf,:))*Mraf(nlraf)
487! Some problems when H is dominant
488! The time step is chosen function of atmospheric mass at higher level
489! It is not very nice
490
491!       if (ig .eq. ij0) then
492!       print*,'test',ig,tdiff,tdiffmin,minloc(Tdiffrafmol),minloc(Tdiffrafmol(nlraf,:))
493!       endif
494        if (tdiff .lt. tdiffmin*Mraf(nlraf)) tdiff=tdiffmin*Mraf(nlraf)
495
496        tdiff=ptimestep/5D0
497       
498! Number of time step
499        ntime=anint(dble(ptimestep)/tdiff)
500!       print*,'ptime',ig,ptimestep,tdiff,ntime,tdiffmin,Mraf(nlraf)
501! Adimensionned temperature
502
503        do l=1,nlraf
504        Tad(l)=Traf(l)/T0
505        enddo
506
507        do istep=1,ntime
508        do nn=1,ncompdiff
509        masse=dble(mmol(gcmind(nn)))
510        Draf(1:nlraf)=Drafmol(1:nlraf,nn)
511
512! Parameters to adimension the problem
513
514        H0=kbolt*T0/masse*invsgmu
515        D0=Draf(nlraf)
516        Time0=H0*H0/D0
517        Time=Tdiff/Time0
518
519! Adimensions
520
521        do l=1,nlraf
522        Dad(l)=Draf(l)/D0
523!        print*,dble(mmol(gcmind(:)))
524       
525        Zad(l)=Zraf(l)/H0
526        Kad(l)=Kraf(l)/D0
527!        print*,'l',l,Zraf(l),Draf(l),Kraf(l),Time0/Dad(l),Time0/Kad(l)
528        enddo
529!        STOP
530        Wad(nn)=wi(nn)*Time0/H0
531        DZ=Zad(2)-Zad(1)
532!       FacEsc=exp(-wad(nn)*DZ/Dad(nlraf-1))
533        Ueff=Wad(nn)
534        alphaTnn=alphaT(nn)
535
536        do l=1,nlraf
537        RhoAd(l)=Rrafk(l,nn)/rho0
538        Dad(l)=Dad(l)/dz/dz
539        Kad(l)=Kad(l)/dz/dz
540        enddo
541
542! Compute intermediary coefficients
543
544        CALL DIFFPARAM(Tad,Dad,Kad,DZ,Rhoad,alphaTnn,delta,ksi,eps,zeta,Mraf,masse,prod,loss,nlraf,Time)
545
546! Compute the alpha and beta recurrent sequences
547
548        CALL SEQUENCY(alpha,beta,delta,ksi,eps,zeta,Dad,Kad,rhoAd,Loss,Prod,Ueff,       &
549     &   dz,time,nlraf)
550
551        Xtri(:)=0D0
552
553! COMPUTE THE DENSITY FROM BOTTOM TO TOP
554
555        Xtri(1)=Prod(1)/Loss(1)
556!        if (ig .eq. ij0) print*,nn,masse
557        DO l=2,nlraf
558        Xtri(l)=(-ALPHA(l-1)+eps(l-1)+zeta(l-1))/(Dad(l-1)+Kad(l-1))*Xtri(l-1)+Beta(l-1)/(Dad(l-1)+Kad(l-1))
559!        if (ig .eq. ij0) print*,'l',l,Xtri(l),rhoAd(l),Prod(l)/Loss(l),ALPHA(l-1),BETA(l-1),eps(l-1),zeta(l-1),Dad(l-1),Kad(l-1)
560        ENDDO
561   
562       
563
564!       Xtri=rhoAd
565
566!       if (ig .eq. ij0 .and. (nn .eq. 1 .or. nn .eq. 5 .or. nn .eq. 6 .or. nn .eq. 16)) then
567!       do l=1,nlraf
568!       if (Xtri(l) .lt. 0.) then
569!       print*,'l',l,rhoAd(l)*rho0,Xtri(l)*rho0,nn,Tad(l),Zad(l),Dad(l)
570!       stop
571!       endif
572!       enddo
573!       endif
574
575! Check mass conservation of speci
576
577!       CALL CheckMass(rhoAd,Xtri,nlraf,nn)
578
579! Update mass density of the species
580
581        do l=1,nlraf
582        Rrafk(l,nn)=rho0*Xtri(l)
583        if (Rrafk(l,nn) .ne. Rrafk(l,nn) .or.                           &
584     &  Rrafk(l,nn) .lt. 0 .and. nn .eq. 16) then
585
586! Test if n(CO2) < 0 skip diffusion (should never happen)
587
588        print*,'pb moldiff',istep,ig,l,nn,Rrafk(l,nn),tdiff,&
589     &  rho0*Rhoad(l),Zmin,Zmax,ntime
590        print*,'Atri',Atri
591        print*,'Btri',Btri
592        print*,'Ctri',Ctri
593        print*,'Dtri',Dtri
594        print*,'Xtri',Xtri
595        print*,'alpha',alpha
596        print*,'beta',beta
597        print*,'delta',delta
598        print*,'eps',eps
599        print*,'Dad',Dad
600        print*,'Temp',Traf
601        print*,'alt',Zraf
602        print*,'Mraf',Mraf
603        stop
604!       pdqdiff(1:ngrid,1:nlayer,1:nq)=0.
605!       return
606!       Rrafk(l,nn)=1D-30*Rraf(l)
607        Rrafk(l,nn)=rho0*Rhoad(l)
608
609        endif
610
611        enddo
612
613        enddo ! END Species Loop
614
615! Update total mass density
616
617        do l=1,nlraf
618        Rraf(l)=sum(Rrafk(l,:))
619        enddo
620
621! Compute new mass average at each altitude level and new pressure
622
623        do l=1,nlraf
624        do nn=1,ncompdiff
625        Qraf(l,nn)=Rrafk(l,nn)/Rraf(l)
626        Nrafk(l,nn)=Rrafk(l,nn)/dble(mmol(gcmind(nn)))/masseU
627        enddo
628        Mraf(l)=1d0/sum(Qraf(l,:)/dble(mmol(gcmind(:))))
629        Nraf(l)=Rraf(l)/Mraf(l)/masseU
630        Praf(l)=Nraf(l)*kbolt*Traf(l)
631        enddo
632
633        enddo ! END time Loop
634
635! Compute the total mass of each species to check mass conservation     
636
637        Mraf2(1:ncompdiff)=0d0
638        do nn=1,ncompdiff
639        do l=1,nlraf
640        Mraf2(nn)=Mraf2(nn)+Rrafk(l,nn)*Dzraf
641        enddo
642        enddo
643
644!        print*,'Mraf',Mraf2
645
646! Reinterpolate values on the GCM pressure levels
647
648        CALL GCMGRID_P2(Zraf,Praf,Qraf,Traf,Nrafk,Rrafk,qq,qnew,tt,tnew,&
649     &   pp,mmol,gcmind,nlraf,ncompdiff,nlayer,FacMass,ig)
650
651        CALL RHOTOT(pp,tt,massemoy,qnew,RHOT,RHOK,nlayer,ncompdiff)
652
653! Update total escape flux of H and H2 (if q was really qnew, but not forget we will output
654! the trend only at the end
655
656!       if (i_h .ne. 1000) PhiEscH=PhiEscH+wi(i_h)*Nrafk(nlraf,i_h)*cell_area(ig) ! in s-1
657        if (i_h .ne. 1000) PhiauxH(ig)=wi(i_h)*Nrafk(nlraf,i_h)*cell_area(ig) ! in s-1
658!       if (i_h2 .ne. 1000) PhiEscH2=PhiEscH2+wi(i_h2)*Nrafk(nlraf,i_h2)*cell_area(ig) ! in s-1 (U in m/s, aire in m2, Nrafk in m-3)
659        if (i_h2 .ne. 1000) PhiauxH2(ig)=wi(i_h2)*Nrafk(nlraf,i_h2)*cell_area(ig)
660!        if (i_d .ne. 1000) PhiEscD=PhiEscD+wi(i_d)*Nrafk(nlraf,i_d)*cell_area(ig)
661        if (i_d .ne. 1000) PhiauxD(ig)=wi(i_d)*Nrafk(nlraf,i_d)*cell_area(ig)
662!       print*,'test',ig,wi(i_h),wi(i_h2),Nrafk(nlraf,i_h),Nrafk(nlraf,i_h2),Nrafk(nlraf,i_d),cell_area(ig),PhiEscH,PhiEscH2,i_h,i_h2,i_d,PhiEscD
663!       stop
664
665
666        if (ig .eq. ij0) then
667        do l=il0,nlayer
668        write(*,'(i2,1x,19(e12.4,1x))') l,zz(l),tt(l),RHOK(l,1)/sum(RHOK(l,:)),RHOKINIT(l,1)/sum(RHOKINIT(l,:)),&
669     &  RHOK(l,2)/sum(RHOK(l,:)),RHOKINIT(l,2)/sum(RHOKINIT(l,:)),&
670     &  RHOK(l,6)/sum(RHOK(l,:)),RHOKINIT(l,6)/sum(RHOKINIT(l,:)),&
671     &  RHOK(l,5)/sum(RHOK(l,:)),RHOKINIT(l,5)/sum(RHOKINIT(l,:)),&
672     &  RHOK(l,7)/sum(RHOK(l,:)),RHOKINIT(l,7)/sum(RHOKINIT(l,:))       
673        enddo
674!        STOP
675        endif
676
677! Compute total mass of each specie on the GCM vertical grid
678
679        Mtot2(1:ncompdiff)=0d0
680
681        do l=il0,nlayer
682        do nn=1,ncompdiff
683        Mtot2(nn)=Mtot2(nn)+1d0/g*qnew(l,nn)*                           &
684     &  (dble(pplev(ig,l))-dble(pplev(ig,l+1)))
685        enddo
686        enddo
687
688! Check mass  conservation of each specie on column
689
690!       do nn=1,ncompdiff
691!       CALL CheckMass2(qq,qnew,pplev(ig,:),il0,nlayer,nn,ncompdiff)
692!       enddo
693
694! Compute the diffusion trends du to diffusion
695
696        do l=1,nlayer
697        do nn=1,ncompdiff
698        pdqdiff(ig,l,gcmind(nn))=(qnew(l,nn)-qq(l,nn))/ptimestep
699        enddo
700        enddo
701
702! deallocation des tableaux
703
704        deallocate(Praf,Traf,Rraf,Mraf)
705        deallocate(Nraf,Draf,Hraf,Wraf)
706        deallocate(Zraf,Tdiffraf,Kraf)
707        deallocate(Prafold,Mrafold)
708        deallocate(Qraf,Rrafk,Nrafk)
709        deallocate(Rrafkold)
710        deallocate(Drafmol,Hrafmol)
711        deallocate(Wrafmol,Tdiffrafmol)
712        deallocate(Atri,Btri,Ctri,Dtri,Xtri)
713        deallocate(Tad,Dad,Kad,Zad,rhoad)
714        deallocate(alpha,beta,delta,ksi,eps,zeta)
715        deallocate(prod,loss)
716
717       enddo  ! ig loop         
718!       print*,'Escape flux H, H2,D (s-1)',PhiEscH,PhiEscH2,PhiEscD
719       if (i_h.ne.1000) call planetwide_sumval(PhiauxH,PhiEscH)
720       if (i_h2.ne.1000) call planetwide_sumval(PhiauxH2,PhiEscH2)
721       if (i_d.ne.1000) call planetwide_sumval(PhiauxD,PhiEscD)
722!       print*,'Escape flux H, H2,D (s-1)',PhiEscH,PhiEscH2,PhiEscD
723
724      end subroutine moldiff_MPF
725
726!    ********************************************************************
727!    ********************************************************************
728!    ********************************************************************
729 
730!     JYC subtroutine solving MX = Y where M is defined as a block tridiagonal
731!       matrix (Thomas algorithm), tested on a example
732
733      subroutine tridagbloc(M,F,X,n1,n2)
734      parameter (nmax=100)
735      real*8 M(n1*n2,n1*n2),F(n1*n2),X(n1*n2)
736      real*8 A(n1,n1,n2),B(n1,n1,n2),C(n1,n1,n2),D(n1,n2)
737      real*8 at(n1,n1),bt(n1,n1),ct(n1,n1),dt(n1),gamt(n1,n1),y(n1,n1)
738      real*8 alf(n1,n1),gam(n1,n1,n2),alfinv(n1,n1)
739      real*8 uvec(n1,n2),uvect(n1),vvect(n1),xt(n1)
740      real*8 indx(n1)
741      real*8 rhu
742      integer n1,n2
743      integer i,p,q   
744
745        X(:)=0.
746! Define the bloc matrix A,B,C and the vector D
747      A(1:n1,1:n1,1)=M(1:n1,1:n1)
748      C(1:n1,1:n1,1)=M(1:n1,n1+1:2*n1)
749      D(1:n1,1)=F(1:n1)
750
751      do i=2,n2-1
752      A(1:n1,1:n1,i)=M((i-1)*n1+1:i*n1,(i-1)*n1+1:i*n1)
753      B(1:n1,1:n1,i)=M((i-1)*n1+1:i*n1,(i-2)*n1+1:(i-1)*n1)
754      C(1:n1,1:n1,i)=M((i-1)*n1+1:i*n1,i*n1+1:(i+1)*n1)
755      D(1:n1,i)=F((i-1)*n1+1:i*n1)
756      enddo
757      A(1:n1,1:n1,n2)=M((n2-1)*n1+1:n2*n1,(n2-1)*n1+1:n2*n1)
758      B(1:n1,1:n1,n2)=M((n2-1)*n1+1:n2*n1,(n2-2)*n1+1:(n2-1)*n1)
759      D(1:n1,n2)=F((n2-1)*n1+1:n2*n1)
760
761! Initialization
762        y(:,:)=0.
763        do i=1,n1
764        y(i,i)=1.
765        enddo
766       
767        at(:,:)=A(:,:,1)
768        ct(:,:)=C(:,:,1)
769        dt(:)=D(:,1)
770        call ludcmp(at,n1,n1,indx,rhu,ierr)
771        do p=1,n1
772         call lubksb(at,n1,n1,indx,y(1,p))
773         do q=1,n1
774         alfinv(q,p)=y(q,p)
775         enddo
776        enddo
777      gamt=matmul(alfinv,ct)
778      gam(:,:,1)=gamt(:,:)     
779      uvect=matmul(alfinv,dt)
780      uvec(:,1)=uvect
781
782      do i=2,n2-1
783        y(:,:)=0.
784       do j=1,n1
785         y(j,j)=1.
786       enddo 
787      bt(:,:)=B(:,:,i)
788      at(:,:)=A(:,:,i)-matmul(bt,gamt)
789      ct(:,:)=C(:,:,i)
790      dt(:)=D(:,i)
791      call ludcmp(at,n1,n1,indx,rhu,ierr)
792        do p=1,n1
793         call lubksb(at,n1,n1,indx,y(1,p))
794         do q=1,n1
795         alfinv(q,p)=y(q,p)
796         enddo
797        enddo
798      gamt=matmul(alfinv,ct)
799      gam(:,:,i)=gamt
800      vvect=dt-matmul(bt,uvect)
801      uvect=matmul(alfinv,vvect)
802      uvec(:,i)=uvect
803      enddo
804      bt=B(:,:,n2)
805      dt=D(:,n2)
806      at=A(:,:,n2)-matmul(bt,gamt)
807      vvect=dt-matmul(bt,uvect)
808       y(:,:)=0.
809       do j=1,n1
810         y(j,j)=1.
811       enddo
812       call ludcmp(at,n1,n1,indx,rhu,ierr)
813        do p=1,n1
814         call lubksb(at,n1,n1,indx,y(1,p))
815         do q=1,n1
816         alfinv(q,p)=y(q,p)
817         enddo
818        enddo
819      xt=matmul(alfinv,vvect)
820      X((n2-1)*n1+1 :n1*n2)=xt
821      do i=n2-1,1,-1
822      gamt=gam(:,:,i)
823      xt=X(i*n1+1:n1*n2)
824      uvect=uvec(:,i)
825      vvect=matmul(gamt,xt)
826      X((i-1)*n1+1:i*n1)=uvect-vvect
827      enddo
828
829      end subroutine tridagbloc
830
831      subroutine tridag(a,b,c,r,u,n)
832!      parameter (nmax=4000)
833!      dimension gam(nmax),a(n),b(n),c(n),r(n),u(n)
834      real*8 gam(n),a(n),b(n),c(n),r(n),u(n)
835      if(b(1).eq.0.)then
836        stop 'tridag: error: b(1)=0 !!! '
837      endif
838      bet=b(1)
839      u(1)=r(1)/bet
840      do 11 j=2,n
841        gam(j)=c(j-1)/bet
842        bet=b(j)-a(j)*gam(j)
843        if(bet.eq.0.) then
844          stop 'tridag: error: bet=0 !!! '
845        endif
846        u(j)=(r(j)-a(j)*u(j-1))/bet
84711    continue
848      do 12 j=n-1,1,-1
849        u(j)=u(j)-gam(j+1)*u(j+1)
85012    continue
851     
852      end subroutine tridag
853
854!    ********************************************************************
855!    ********************************************************************
856!    ********************************************************************
857
858      SUBROUTINE LUBKSB(A,N,NP,INDX,B)
859
860      implicit none
861
862      integer i,j,n,np,ii,ll
863      real*8 sum
864      real*8 a(np,np),indx(np),b(np)
865
866!      DIMENSION A(NP,NP),INDX(N),B(N)
867      II=0
868      DO 12 I=1,N
869        LL=INDX(I)
870        SUM=B(LL)
871        B(LL)=B(I)
872        IF (II.NE.0)THEN
873          DO 11 J=II,I-1
874            SUM=SUM-A(I,J)*B(J)
87511        CONTINUE
876        ELSE IF (SUM.NE.0.) THEN
877          II=I
878        ENDIF
879        B(I)=SUM
88012    CONTINUE
881      DO 14 I=N,1,-1
882        SUM=B(I)
883        IF(I.LT.N)THEN
884          DO 13 J=I+1,N
885            SUM=SUM-A(I,J)*B(J)
88613        CONTINUE
887        ENDIF
888        B(I)=SUM/A(I,I)
88914    CONTINUE
890     
891      END SUBROUTINE LUBKSB
892
893!    ********************************************************************
894!    ********************************************************************
895!    ********************************************************************
896
897      SUBROUTINE LUDCMP(A,N,NP,INDX,D,ierr)
898
899      implicit none
900
901      integer n,np,nmax,i,j,k,imax
902      real*8 d,tiny,aamax
903      real*8 a(np,np),indx(np)
904      integer ierr  ! error =0 if OK, =1 if problem
905
906      PARAMETER (NMAX=100,TINY=1.0E-20)
907!      DIMENSION A(NP,NP),INDX(N),VV(NMAX)
908      real*8 sum,vv(nmax),dum
909
910      D=1.
911      DO 12 I=1,N
912        AAMAX=0.
913        DO 11 J=1,N
914          IF (ABS(A(I,J)).GT.AAMAX) AAMAX=ABS(A(I,J))
91511      CONTINUE
916        IF (AAMAX.EQ.0.) then
917            write(*,*) 'In moldiff: Problem in LUDCMP with matrix A'
918            write(*,*) 'Singular matrix ?'
919            write(*,*) 'Matrix A = ', A
920!            stop
921!           TO DEBUG :
922            ierr =1
923            return
924!           stop
925        END IF
926
927        VV(I)=1./AAMAX
92812    CONTINUE
929      DO 19 J=1,N
930        IF (J.GT.1) THEN
931          DO 14 I=1,J-1
932            SUM=A(I,J)
933            IF (I.GT.1)THEN
934              DO 13 K=1,I-1
935                SUM=SUM-A(I,K)*A(K,J)
93613            CONTINUE
937              A(I,J)=SUM
938            ENDIF
93914        CONTINUE
940        ENDIF
941        AAMAX=0.
942        DO 16 I=J,N
943          SUM=A(I,J)
944          IF (J.GT.1)THEN
945            DO 15 K=1,J-1
946              SUM=SUM-A(I,K)*A(K,J)
94715          CONTINUE
948            A(I,J)=SUM
949          ENDIF
950          DUM=VV(I)*ABS(SUM)
951          IF (DUM.GE.AAMAX) THEN
952            IMAX=I
953            AAMAX=DUM
954          ENDIF
95516      CONTINUE
956        IF (J.NE.IMAX)THEN
957          DO 17 K=1,N
958            DUM=A(IMAX,K)
959            A(IMAX,K)=A(J,K)
960            A(J,K)=DUM
96117        CONTINUE
962          D=-D
963          VV(IMAX)=VV(J)
964        ENDIF
965        INDX(J)=IMAX
966        IF(J.NE.N)THEN
967          IF(A(J,J).EQ.0.)A(J,J)=TINY
968          DUM=1./A(J,J)
969          DO 18 I=J+1,N
970            A(I,J)=A(I,J)*DUM
97118        CONTINUE
972        ENDIF
97319    CONTINUE
974      IF(A(N,N).EQ.0.)A(N,N)=TINY
975      ierr =0
976     
977      END SUBROUTINE LUDCMP
978
979        SUBROUTINE TMNEW(T1,DT1,DT2,DT3,T2,dtime,nl,ig)
980        IMPLICIT NONE
981
982        INTEGER,INTENT(IN) :: nl,ig
983        REAL,INTENT(IN),DIMENSION(nl) :: T1,DT1,DT2,DT3
984        REAL*8,INTENT(OUT),DIMENSION(nl) :: T2
985        REAL,INTENT(IN) :: dtime
986        INTEGER :: l
987        DO l=1,nl
988        T2(l)=T1(l)*1D0+(DT1(l)*dble(dtime)+                            &
989     &  DT2(l)*dble(dtime)+                                             &
990     &  DT3(l)*dble(dtime))*1D0
991        if (T2(l) .ne. T2(l)) then
992        print*,'Err TMNEW',ig,l,T2(l),T1(l),dT1(l),DT2(l),              &
993     &  DT3(l),dtime,dble(dtime)
994        endif
995
996        ENDDO
997        END SUBROUTINE TMNEW
998
999        SUBROUTINE QMNEW(Q1,DQ,Q2,dtime,nl,nq,gc,ig)
1000        use tracer_mod, only: nqmx
1001        IMPLICIT NONE
1002
1003        INTEGER,INTENT(IN) :: nl,nq
1004        INTEGER,INTENT(IN) :: ig
1005        INTEGER,INTENT(IN),dimension(nq) :: gc
1006        REAL,INTENT(IN),DIMENSION(nl,nqmx) :: Q1,DQ
1007        REAL*8,INTENT(OUT),DIMENSION(nl,nq) :: Q2
1008        REAL,INTENT(IN) :: dtime
1009        INTEGER :: l,iq
1010        DO l=1,nl
1011        DO iq=1,nq
1012        Q2(l,iq)=Q1(l,gc(iq))*1D0+(DQ(l,gc(iq))*dble(dtime))*1D0
1013        Q2(l,iq)=max(Q2(l,iq),1d-30)
1014        ENDDO
1015        ENDDO
1016        END SUBROUTINE QMNEW
1017
1018        SUBROUTINE HSCALE(p,hp,nl)
1019        IMPLICIT NONE
1020
1021        INTEGER :: nl
1022        INTEGER :: l
1023        REAL*8,dimension(nl) :: P
1024        REAL*8,DIMENSION(nl) :: Hp
1025       
1026        hp(1)=-log(P(2)/P(1))
1027        hp(nl)=-log(P(nl)/P(nl-1))
1028
1029        DO l=2,nl-1
1030        hp(l)=-log(P(l+1)/P(l-1))
1031        ENDDO
1032        END SUBROUTINE HSCALE
1033
1034        SUBROUTINE MMOY(massemoy,mmol,qq,gc,nl,nq)
1035        use tracer_mod, only: nqmx
1036        IMPLICIT NONE
1037
1038        INTEGER :: nl,nq,l
1039        INTEGER,dimension(nq) :: gc
1040        REAL*8,DIMENSION(nl,nq) :: qq
1041        REAL*8,DIMENSION(nl) :: massemoy
1042        REAL,DIMENSION(nqmx) :: MMOL
1043
1044
1045        do l=1,nl
1046        massemoy(l)=1D0/sum(qq(l,:)/dble(mmol(gc(:))))
1047        enddo
1048
1049        END SUBROUTINE MMOY
1050
1051        SUBROUTINE DMMOY(M,H,DM,nl)
1052        IMPLICIT NONE   
1053        INTEGER :: nl,l
1054        REAL*8,DIMENSION(nl) :: M,H,DM
1055
1056        DM(1)=(-3D0*M(1)+4D0*M(2)-M(3))/2D0/H(1)
1057        DM(nl)=(3D0*M(nl)-4D0*M(nl-1)+M(nl-2))/2D0/H(nl)
1058
1059        do l=2,nl-1
1060        DM(l)=(M(l+1)-M(l-1))/H(l)
1061        enddo
1062
1063        END SUBROUTINE DMMOY
1064
1065        SUBROUTINE ZVERT(P,T,M,Z,nl,ig)
1066        IMPLICIT NONE
1067        INTEGER :: nl,l,ig
1068        REAL*8,dimension(nl) :: P,T,M,Z,H
1069        REAL*8 :: z0
1070        REAL*8 :: kbolt,masseU,Konst,g,Hpm
1071        masseU=1.660538782d-27
1072        kbolt=1.3806504d-23
1073        Konst=kbolt/masseU
1074        g=3.72D0
1075
1076        z0=0d0
1077        Z(1)=z0
1078        H(1)=Konst*T(1)/M(1)/g
1079
1080        do l=2,nl
1081        H(l)=Konst*T(l)/M(l)/g
1082        Hpm=H(l-1)
1083        Z(l)=z(l-1)-Hpm*log(P(l)/P(l-1))
1084        if (Z(l) .ne. Z(l)) then
1085        print*,'pb',l,ig
1086        print*,'P',P
1087        print*,'T',T
1088        print*,'M',M
1089        print*,'Z',Z
1090        print*,'Hpm',Hpm
1091        endif
1092        enddo
1093
1094        END SUBROUTINE ZVERT
1095
1096        SUBROUTINE RHOTOT(P,T,M,qq,rhoN,rhoK,nl,nq)
1097        IMPLICIT NONE
1098
1099        REAL*8 :: kbolt,masseU,Konst
1100        INTEGER :: nl,nq,l,iq
1101        REAL*8,DIMENSION(nl) :: P,T,M,RHON
1102        REAL*8,DIMENSION(nl,nq) :: RHOK,qq
1103        masseU=1.660538782d-27
1104        kbolt=1.3806504d-23
1105        Konst=Kbolt/masseU
1106
1107        do l=1,nl
1108        RHON(l)=P(l)*M(l)/T(l)/Konst
1109        do iq=1,nq
1110        RHOK(l,iq)=qq(l,iq)*RHON(l)
1111        enddo
1112        enddo
1113
1114        END SUBROUTINE RHOTOT
1115
1116        SUBROUTINE UPPER_RESOL(P,T,Z,M,R,Rk,                            &
1117     & qq,mmol,gc,Praf,Traf,Qraf,Mraf,Zraf,                             &
1118     & Nraf,Nrafk,Rraf,Rrafk,il,nl,nq,nlx,ig)
1119        use tracer_mod, only: nqmx
1120        IMPLICIT NONE
1121       
1122        INTEGER :: nl,nq,il,l,i,iq,nlx,iz,ig
1123        INTEGER :: gc(nq)
1124        INTEGER,DIMENSION(1) :: indz
1125        REAL*8, DIMENSION(nlx) :: P,T,Z,M,R
1126        REAL*8, DIMENSION(nlx,nq) :: qq,Rk
1127        REAL*8, DIMENSION(nl) :: Praf,Traf,Mraf,Zraf,Nraf,Rraf
1128        REAL*8 :: kbolt,masseU,Konst,g
1129        REAL*8, DIMENSION(nl,nq) :: Qraf,Rrafk,Nrafk
1130        REAL*8 :: facZ,dZ,H
1131        REAL,DIMENSION(nqmx) :: mmol
1132        masseU=1.660538782d-27
1133        kbolt=1.3806504d-23
1134        Konst=Kbolt/masseU
1135        g=3.72d0
1136
1137
1138        Zraf(1)=z(il)
1139        Praf(1)=P(il)
1140        Traf(1)=T(il)
1141        Nraf(1)=Praf(1)/kbolt/Traf(1)
1142        do iq=1,nq
1143        Qraf(1,iq)=qq(il,iq)
1144        enddo
1145        Mraf(1)=1d0/sum(Qraf(1,:)/dble(mmol(gc(:))))
1146        Rraf(1)=Mraf(1)*masseU*Nraf(1)
1147        do iq=1,nq
1148        Rrafk(1,iq)=Rraf(1)*Qraf(1,iq)
1149        Nrafk(1,iq)=Rrafk(1,iq)/masseU/dble(mmol(gc(iq)))
1150        enddo
1151        Zraf(nl)=z(nlx)
1152
1153        do l=2,nl-1
1154        Zraf(l)=Zraf(1)+(Zraf(nl)-Zraf(1))/dble(nl-1)*dble((l-1))
1155        indz=maxloc(z,mask=z <= Zraf(l))
1156        iz=indz(1)
1157        if (iz .lt. 1 .or. iz .gt. nlx) then
1158        print*,'danger',iz,nl,Zraf(l),l,Zraf(1),Zraf(nl),z,P,T,ig
1159        stop
1160        endif
1161        dZ=Zraf(l)-Zraf(l-1)
1162!       dZ=Zraf(l)-z(iz)
1163        facz=(Zraf(l)-z(iz))/(z(iz+1)-z(iz))
1164        Traf(l)=T(iz)+(T(iz+1)-T(iz))*facz
1165        do iq=1,nq
1166!        Qraf(l,iq)=qq(iz,iq)+(qq(iz+1,iq)-qq(iz,iq))*facz
1167        Rrafk(l,iq)=Rk(iz,iq)+(Rk(iz+1,iq)-Rk(iz,iq))*facZ
1168        Rrafk(l,iq)=Rk(iz,iq)*(Rk(iz+1,iq)/Rk(iz,iq))**facZ
1169        enddo
1170!       Mraf(l)=1D0/(sum(qraf(l,:)/dble(mmol(gc(:)))))
1171        Rraf(l)=sum(Rrafk(l,:))
1172        do iq=1,nq
1173        Qraf(l,iq)=Rrafk(l,iq)/Rraf(l)
1174        enddo
1175        Mraf(l)=1D0/(sum(qraf(l,:)/dble(mmol(gc(:)))))
1176!       H=Konst*Traf(l)/Mraf(l)/g
1177!       H=Konst*T(iz)/M(iz)/g
1178!       Praf(l)=P(iz)*exp(-dZ/H)
1179!       Praf(l)=Praf(l-1)*exp(-dZ/H)
1180!       print*,'iz',l,iz,Praf(il-1)*exp(-dZ/H),z(iz),z(iz+1),H
1181        Nraf(l)=Rraf(l)/Mraf(l)/masseU
1182        Praf(l)=Nraf(l)*kbolt*Traf(l)
1183!       Rraf(l)=Nraf(l)*Mraf(l)*masseU
1184        do iq=1,nq
1185!       Rrafk(l,iq)=Rraf(l)*Qraf(l,iq)
1186        Nrafk(l,iq)=Rrafk(l,iq)/dble(mmol(gc(iq)))/masseU
1187        if (Nrafk(l,iq) .lt. 0. .or.                                    &
1188     &  Nrafk(l,iq) .ne. Nrafk(l,iq)) then
1189        print*,'pb interpolation',l,iq,Nrafk(l,iq),Rrafk(l,iq),         &
1190     &  Qraf(l,iq),Rk(iz,iq),Rk(iz+1,iq),facZ,Zraf(l),z(iz)
1191        stop
1192        endif
1193        enddo
1194        enddo
1195        Zraf(nl)=z(nlx)
1196        Traf(nl)=T(nlx)
1197        do iq=1,nq
1198        Rrafk(nl,iq)=Rk(nlx,iq)
1199        Qraf(nl,iq)=Rk(nlx,iq)/R(nlx)
1200        Nrafk(nl,iq)=Rk(nlx,iq)/dble(mmol(gc(iq)))/masseU
1201        enddo
1202        Nraf(nl)=sum(Nrafk(nl,:))
1203        Praf(nl)=Nraf(nl)*kbolt*Traf(nl)
1204        Mraf(nl)=1D0/sum(Qraf(nl,:)/dble(mmol(gc(:))))
1205        END SUBROUTINE UPPER_RESOL
1206
1207        SUBROUTINE CORRMASS(qq,qint,FacMass,nl,nq)
1208        IMPLICIT NONE
1209        INTEGER :: nl,nq,l,nn
1210        REAL*8,DIMENSION(nl,nq) :: qq,qint,FacMass
1211
1212        do nn=1,nq
1213        do l=1,nl
1214        FacMass(l,nn)=qq(l,nn)/qint(l,nn)
1215        enddo
1216        enddo
1217
1218        END     
1219
1220
1221        SUBROUTINE DCOEFF(nn,Dij,P,T,N,Nk,D,nl,nq)
1222        IMPLICIT NONE
1223        REAL*8,DIMENSION(nl) :: N,T,D,P
1224        REAL*8,DIMENSION(nl,nq) :: Nk
1225        INTEGER :: nn,nl,nq,l,iq
1226        REAL,DIMENSION(nq,nq)  :: Dij
1227        REAL*8 :: interm,P0,T0,ptfac,dfac
1228
1229        P0=1D5
1230        T0=273d0
1231       
1232
1233        do l=1,nl
1234        ptfac=(P0/P(l))*(T(l)/T0)**1.75d0
1235        D(l)=0d0
1236        interm=0d0
1237        do iq=1,nq
1238        if (iq .ne. nn) then
1239        dfac=dble(dij(nn,iq))*ptfac
1240        interm=interm+Nk(l,iq)/N(l)/dfac
1241        endif
1242        enddo
1243        !Temporary: eliminate modification to include Wilke's formulation
1244        !back to the old scheme to check effect
1245        !D(l)=1d0/interm
1246        D(l)=(1D0-Nk(l,nn)/N(l))/interm
1247        enddo
1248        END SUBROUTINE DCOEFF
1249 
1250        SUBROUTINE HSCALEREAL(nn,Nk,Dz,H,nl,nq)
1251        IMPLICIT NONE
1252        INTEGER :: nn,nl,nq,l
1253        REAL*8,DIMENSION(nl) :: H
1254        REAL*8,DIMENSION(nl,nq) :: Nk
1255        REAL*8 :: Dz
1256       
1257        H(1)=(-3D0*Nk(1,nn)+4d0*NK(2,nn)-Nk(3,nn))/(2D0*DZ)/            &
1258     &  NK(1,nn)
1259
1260        H(1)=-1D0/H(1)
1261
1262        DO l=2,nl-1
1263        H(l)=(Nk(l+1,nn)-NK(l-1,nn))/(2D0*DZ)/NK(l,nn)
1264        H(l)=-1D0/H(l)
1265        ENDDO
1266
1267        H(nl)=(3D0*Nk(nl,nn)-4D0*Nk(nl-1,nn)+Nk(nl-2,nn))/(2D0*DZ)/     &
1268     &  Nk(nl,nn)
1269        H(nl)=-1D0/H(nl)
1270
1271!       do l=1,nl
1272!       if (abs(H(l)) .lt. 100.) then
1273!       print*,'H',l,H(l),Nk(l,nn),nn
1274!       endif
1275!       enddo
1276
1277        END SUBROUTINE HSCALEREAL
1278       
1279        SUBROUTINE VELVERT(nn,T,H,D,Dz,masse,W,nl)
1280        IMPLICIT NONE
1281        INTEGER :: l,nl,nn
1282        REAL*8,DIMENSION(nl) :: T,H,D,W,DT
1283        REAL*8 :: Dz,Hmol,masse
1284        REAL*8 :: kbolt,masseU,Konst,g
1285        masseU=1.660538782d-27
1286        kbolt=1.3806504d-23
1287        Konst=Kbolt/masseU
1288        g=3.72d0
1289
1290        DT(1)=1D0/T(1)*(-3D0*T(1)+4D0*T(2)-T(3))/(2D0*DZ)
1291        Hmol=Konst*T(1)/masse/g
1292        W(1)=-D(1)*(1D0/H(1)-1D0/Hmol-DT(1))
1293
1294        DO l=2,nl-1
1295        DT(l)=1D0/T(l)*(T(l+1)-T(l-1))/(2D0*DZ)
1296        Hmol=Konst*T(l)/masse/g
1297        W(l)=-D(l)*(1D0/H(l)-1D0/Hmol-DT(l))
1298        ENDDO
1299
1300        DT(nl)=1D0/T(nl)*(3d0*T(nl)-4D0*T(nl-1)+T(nl-2))/(2D0*DZ)
1301        Hmol=Konst*T(nl)/masse/g
1302        W(nl)=-D(nl)*(1D0/H(nl)-1D0/Hmol-DT(nl))
1303
1304!       do l=1,nl
1305!       print*,'W',W(l),D(l),H(l),DT(l)
1306!       enddo
1307
1308        END SUBROUTINE VELVERT
1309
1310        SUBROUTINE TIMEDIFF(nn,H,W,TIME,nl)
1311        IMPLICIT NONE
1312        INTEGER :: nn,nl,l
1313        REAL*8,DIMENSION(nl) :: W,H,TIME
1314
1315        DO l=1,nl
1316        TIME(l)=abs(H(l)/W(l))
1317        if (TIME(l) .lt. 1.D-4) then
1318!       print*,'low tdiff',H(l),W(l),nn,l
1319        endif
1320        ENDDO
1321
1322        END SUBROUTINE TIMEDIFF
1323
1324
1325        SUBROUTINE DIFFPARAM(T,D,K,dz,RHO,alphaTnn,delta,ksi,eps,zeta,Ma,mi,prod,loss,nl,dtime)
1326        IMPLICIT NONE
1327        INTEGER :: nl,l
1328        REAL*8,DIMENSION(nl) :: T,D,K,RHO,Ma
1329        REAL*8 :: DZ,DZinv,dT,dtime,mi,alphaTnn
1330        REAL*8,DIMENSION(nl) :: delta,ksi,eps,zeta,prod,loss
1331
1332! Compute the vectors delta,eps,prod and loss
1333        DO l=1,nl-1
1334        dT=(1D0/T(l)*(T(l+1)-T(l)))/dZ
1335        delta(l)=(1D0/T(l)+(1D0+alphaTnn)*dT)*D(l)*dz
1336        ksi(l)=(1D0/T(l)*Ma(l)/mi+dT)*K(l)*dz
1337        eps(l)=D(l)-delta(l)
1338        zeta(l)=K(l)-ksi(l)
1339        prod(l)=RHO(l)/dtime
1340        loss(l)=1D0/dtime
1341!        print*,l,dT,delta(l),ksi(l),eps(l),zeta(l),D(l),K(l),eps(l)/D(l),zeta(l)/K(l),Ma(l),mi,dZ
1342        ENDDO
1343       
1344! at top assume T = Cste (dT=0)
1345        delta(nl)=1D0/T(nl)*D(nl)*dz
1346        ksi(nl)=1D0/T(nl)*Ma(nl)/mi*K(nl)*dz
1347        eps(nl)=D(nl)-delta(nl)
1348        zeta(nl)=K(nl)-ksi(nl)
1349        prod(nl)=RHO(nl)/dtime
1350        loss(nl)=1D0/dtime       
1351        END SUBROUTINE DIFFPARAM
1352
1353
1354        SUBROUTINE SEQUENCY(alpha,beta,delta,ksi,eps,zeta,Dad,Kad,rhoad,Loss,Prod,        &
1355     &  Ueff,dz,dt,nl)
1356        IMPLICIT NONE
1357        INTEGER :: nl,l
1358        REAL*8, DIMENSION(nl) :: alpha,beta,delta,ksi,eps,zeta,Dad,Kad,RHoad,Loss,Prod
1359        REAL*8 :: dz,dt,del1,del2,del3,Ueff
1360
1361        ALPHA(nl-1)=(eps(nl-1)+zeta(nl-1))*Ueff/dZ/(Dad(nl-1)+Kad(nl-1)+Ueff/dZ)
1362!        ALPHA(nl-1)=Ueff/dZ
1363        BETA(nl-1)=0D0
1364
1365        DO l=nl-2,1,-1
1366!        print*,l,eps(l),zeta(l),Dad(l),Kad(l),ALPHA(l+1),BETA(l+1),Prod(l+1),Loss(l+1)
1367        ALPHA(l)=(eps(l)+zeta(l))*(ALPHA(l+1)+Loss(l+1))/(ALPHA(l+1)+Dad(l)+Kad(l)+Loss(l+1))
1368        BETA(l)=(Dad(l)+Kad(l))*(BETA(l+1)+Prod(l+1))/(ALPHA(l+1)+Dad(l)+Kad(l)+Loss(l+1))       
1369        ENDDO
1370
1371
1372        END SUBROUTINE SEQUENCY
1373
1374        SUBROUTINE Checkmass(X,Y,nl,nn)
1375        IMPLICIT NONE
1376
1377        INTEGER :: nl,nn
1378        REAL*8,DIMENSION(nl) :: X,Y
1379        REAL*8 Xtot,Ytot
1380
1381        Xtot=sum(X)
1382        Ytot=sum(Y)
1383
1384        if (abs((Xtot-Ytot)/Xtot) .gt. 1d-3) then
1385        print*,'no conservation for mass',Xtot,Ytot,nn
1386        endif
1387        END SUBROUTINE Checkmass
1388
1389        SUBROUTINE Checkmass2(qold,qnew,P,il,nl,nn,nq)
1390        IMPLICIT NONE
1391        INTEGER :: nl,nn,l,nq,il
1392        REAL,DIMENSION(nl+1) :: P
1393        REAL*8,DIMENSION(nl,nq) :: qold,qnew
1394        REAL*8 :: DM,Mold,Mnew,g
1395        g=3.72d0
1396        DM=0d0
1397        Mold=0d0
1398        Mnew=0d0
1399        DO l=il,nl
1400        DM=DM+(qnew(l,nn)-qold(l,nn))*(dble(P(l))-dble(P(l+1)))/g
1401        Mold=Mold+qold(l,nn)*(dble(P(l))-dble(P(l+1)))/g
1402        Mnew=Mnew+qnew(l,nn)*(dble(P(l))-dble(P(l+1)))/g
1403!       print*,'l',l,qold(l,nn),qnew(l,nn),Mold,Mnew,DM,P(l),P(l+1)
1404        ENDDO
1405        IF (abs(DM/Mold) .gt. 1d-2) THEN
1406        Print*,'We dont conserve mas',nn,DM,Mold,Mnew,DM/Mold
1407        ENDIF
1408
1409        END SUBROUTINE Checkmass2
1410
1411        SUBROUTINE GCMGRID_P(Z,P,Q,T,Nk,Rk,qq,qnew,tt,tnew,             &
1412     &    pp,M,gc,nl,nq,nlx,ig)
1413        use tracer_mod, only: nqmx
1414        IMPLICIT NONE
1415        INTEGER :: nl,nq,nlx,il,nn,iP,ig,compteur
1416        INTEGER,DIMENSION(1) :: indP
1417        INTEGER,DIMENSION(nq) :: gc
1418        REAL*8,DIMENSION(nl) :: Z,P,T
1419        REAL*8,DIMENSION(nl,nq) :: Q,Nk,Rk
1420        REAL,DIMENSION(nqmx) :: M
1421        REAL*8,DIMENSION(nq) :: nNew
1422        REAL*8,DIMENSION(nlx) :: pp,tt,tnew
1423        REAL*8,DIMENSION(nlx) :: rhonew
1424        REAL*8,DIMENSION(nlx,nq) :: qq,qnew,rhoknew
1425        REAL*8 :: kbolt,masseU,Konst,g,Dz,facP,Hi
1426        REAL*8 :: Znew,Znew2,Pnew,Pnew2
1427        masseU=1.660538782d-27
1428        kbolt=1.3806504d-23
1429        Konst=Kbolt/masseU
1430        g=3.72d0
1431        Dz=Z(2)-Z(1)
1432        Znew=Z(nl)
1433        Znew2=Znew+dz
1434!       print*,'dz',Znew,Znew2,dz
1435        nNew(1:nq)=Nk(nl,1:nq)
1436        Pnew=P(nl)
1437
1438        do il=1,nlx
1439!       print*,'il',il,pp(il),P(1),P(nl)
1440        if (pp(il) .ge. P(1)) then
1441        qnew(il,:)=qq(il,:)
1442        tnew(il)=tt(il)
1443        endif
1444        if (pp(il) .lt. P(1)) then
1445        if (pp(il) .gt. P(nl)) then
1446        indP=maxloc(P,mask=P <  pp(il))
1447        iP=indP(1)-1
1448        if (iP .lt. 1 .or. iP .gt. nl) then
1449        print*,'danger 2',iP,nl,pp(il)
1450        endif
1451        facP=(pp(il)-P(ip))/(P(ip+1)-P(ip))
1452!       print*,'P',P(ip),P(ip+1),facP,indP,iP
1453
1454!       do nn=1,nq
1455!       qnew(il,nn)=Q(iP,nn)+
1456!     &  (Q(ip+1,nn)-Q(ip,nn))*facP
1457!       enddo
1458
1459        do nn=1,nq
1460        rhoknew(il,nn)=Rk(iP,nn)+                                       &
1461     &  (Rk(ip+1,nn)-Rk(ip,nn))*facP
1462        enddo
1463        tnew(il)=T(iP)+(T(iP+1)-T(iP))*facP
1464        rhonew(il)=sum(rhoknew(il,:))
1465        do nn=1,nq
1466        qnew(il,nn)=rhoknew(il,nn)/rhonew(il)
1467        enddo
1468
1469        else ! pp < P(nl) need to extrapolate density of each specie
1470        Pnew2=Pnew
1471
1472        compteur=0
1473        do while (pnew2 .ge. pp(il))   
1474        compteur=compteur+1
1475        do nn=1,nq
1476        Hi=Konst*T(nl)/dble(M(gc(nn)))/g
1477        Nnew(nn)=Nnew(nn)*exp(-dZ/Hi)
1478        enddo
1479        Pnew=Pnew2
1480        Pnew2=kbolt*T(nl)*sum(Nnew(:))
1481        Znew=Znew2
1482        Znew2=Znew2+Dz
1483        if (compteur .ge. 100000) then
1484        print*,'error moldiff_MPF infinite loop'
1485        print*,ig,il,pp(il),tt(nl),pnew2,qnew(il,:),Znew2
1486        stop
1487        endif
1488!       print*,'test',Pnew2,Znew2,Nnew(nq),pp(il)
1489        enddo
1490       
1491        facP=(pp(il)-Pnew)/(Pnew2-Pnew)
1492
1493!       do nn=1,nq
1494!       qnew(il,nn)=dble(M(gc(nn)))*Nnew(nn)
1495!     &  /sum(dble(M(gc(:)))*Nnew(:))
1496!       enddo
1497
1498        do nn=1,nq
1499        rhoknew(il,nn)=dble(M(gc(nn)))*Nnew(nn)
1500        enddo
1501        rhonew(il)=sum(rhoknew(il,:))
1502        do nn=1,nq
1503        qnew(il,nn)=rhoknew(il,nn)/rhonew(il)
1504        enddo
1505        tnew(il)=T(nl)
1506        endif
1507        endif
1508        enddo
1509
1510        END SUBROUTINE GCMGRID_P
1511
1512        SUBROUTINE GCMGRID_P2(Z,P,Q,T,Nk,Rk,qq,qnew,tt,tnew             &
1513     &   ,pp,M,gc,nl,nq,nlx,facM,ig)
1514        use tracer_mod, only: nqmx
1515        IMPLICIT NONE
1516        INTEGER :: nl,nq,nlx,il,nn,iP,ig,compteur
1517        INTEGER,DIMENSION(1) :: indP
1518        INTEGER,DIMENSION(nq) :: gc
1519        REAL*8,DIMENSION(nl) :: Z,P,T
1520        REAL*8,DIMENSION(nl,nq) :: Q,Nk,Rk
1521        REAL,DIMENSION(nqmx) :: M
1522        REAL*8,DIMENSION(nq) :: nNew
1523        REAL*8,DIMENSION(nlx) :: pp,rhonew,tt,tnew
1524        REAL*8,DIMENSION(nlx,nq) :: qq,qnew,facM,rhoknew
1525        REAL*8 :: kbolt,masseU,Konst,g,Dz,facP,Hi
1526        REAL*8 :: Znew,Znew2,Pnew,Pnew2
1527        masseU=1.660538782d-27
1528        kbolt=1.3806504d-23
1529        Konst=Kbolt/masseU
1530        g=3.72d0
1531        Dz=Z(2)-Z(1)
1532        Znew=Z(nl)
1533        Znew2=Znew+dz
1534!       print*,'dz',Znew,Znew2,dz
1535        nNew(1:nq)=Nk(nl,1:nq)
1536        Pnew=P(nl)
1537
1538        do il=1,nlx
1539!       print*,'il',il,pp(il),P(1),P(nl)
1540        if (pp(il) .ge. P(1)) then
1541        qnew(il,:)=qq(il,:)
1542        tnew(il)=tt(il)
1543        endif
1544        if (pp(il) .lt. P(1)) then
1545        if (pp(il) .gt. P(nl)) then
1546        indP=maxloc(P,mask=P <  pp(il))
1547        iP=indP(1)-1
1548        if (iP .lt. 1 .or. iP .gt. nl) then
1549        print*,'danger 3',iP,nl,pp(il)
1550        endif
1551        facP=(pp(il)-P(ip))/(P(ip+1)-P(ip))
1552!       print*,'P',P(ip),P(ip+1),facP,indP,iP
1553
1554!        do nn=1,nq
1555!        qnew(il,nn)=Q(iP,nn)+
1556!     &  (Q(ip+1,nn)-Q(ip,nn))*facP
1557!        enddo
1558
1559        do nn=1,nq
1560        rhoknew(il,nn)=(RK(iP,nn)+                                       &
1561     &  (RK(iP+1,nn)-Rk(iP,nn))*facP)*facM(il,nn)
1562        enddo
1563        tnew(il)=T(iP)+(T(ip+1)-T(iP))*facP
1564        rhonew(il)=sum(rhoknew(il,:))
1565        do nn=1,nq
1566        qnew(il,nn)=rhoknew(il,nn)/rhonew(il)
1567        enddo
1568
1569        else ! pp < P(nl) need to extrapolate density of each specie
1570        Pnew2=Pnew
1571
1572        compteur=0
1573        do while (pnew2 .ge. pp(il))   
1574        compteur=compteur+1
1575        do nn=1,nq
1576        Hi=Konst*T(nl)/dble(M(gc(nn)))/g
1577        Nnew(nn)=Nnew(nn)*exp(-dZ/Hi)
1578        enddo
1579        Pnew=Pnew2
1580        Pnew2=kbolt*T(nl)*sum(Nnew(:))
1581        Znew=Znew2
1582        Znew2=Znew2+Dz
1583        if (compteur .ge. 100000) then
1584        print*,'pb moldiff_MPF infinite loop'
1585        print*,ig,nl,T(nl),pnew2,qnew(il,:),Znew2
1586        stop
1587        endif
1588
1589!       print*,'test',Pnew2,Znew2,Nnew(nq),pp(il)
1590        enddo
1591       
1592        facP=(pp(il)-Pnew)/(Pnew2-Pnew)
1593
1594!        do nn=1,nq
1595!        qnew(il,nn)=dble(M(gc(nn)))*Nnew(nn)
1596!     &  /sum(dble(M(gc(:)))*Nnew(:))
1597!        enddo
1598
1599        do nn=1,nq
1600        rhoknew(il,nn)=dble(M(gc(nn)))*Nnew(nn)*FacM(il,nn)
1601        enddo
1602        rhonew(il)=sum(rhoknew(il,:))
1603        do nn=1,nq
1604        qnew(il,nn)=rhoknew(il,nn)/rhonew(il)
1605        enddo
1606        tnew(il)=T(nl)
1607
1608        endif
1609        endif
1610        enddo
1611
1612        END SUBROUTINE GCMGRID_P2
1613        END MODULE moldiff_MPF_mod       
Note: See TracBrowser for help on using the repository browser.