source: trunk/LMDZ.VENUS/libf/phyvenus/Iondiff_red.F @ 3461

Last change on this file since 3461 was 3035, checked in by amartinez, 15 months ago

======= VENUS PCM ===============

COMMIT BY ANTOINE MARTINEZ

Implementation of vertical ambipolar diffusion in physics

=================================

NEW KEYWORD OF PHYSIQ.DEF

=================================

NEW VERSION OF "physiq.def"

  • deftank/physiq-96x96x90-chemistry-IONOSPHERE-IONDIFF.def
  • ok_iondiff: keyword supposed to activate ion ambipolar diffusion
  • nbapp_chem: replaces nbapp_chim, in order to characterize the Number of calls to the chemistry routines (per Venusian day)

================

phyvenus

================

Iondiff_red.F

  • Calculation of the Ion Ambipolar Diffusion for 13 ions on 14:

O2+, O+, C+, N+, CO2+,
NO+, CO+, H+, H2O+, H3O+,

HCO+, N2+, OH+


The ion temperature is fixed as the half of the electron temperature
calculated in ion_h.F for the stability of the program and because the
ion temperature is lower than electron temperature.

plasmaphys_venus_mod.F90

  • parameters of the ambipolar diffusion scheme:

parameter (Pdiffion = 7.0D-4) ! pressure in Pa below which ion diffusion is computed
parameter (naltvert = 168) ! number of level on the altimetric subgrid
parameter (nsubvert = 24000) ! ptimephysiq/nsubvert - minimum sub-timestep allowed
parameter ( nsubmin = 40) ! ptimephysiq/nsubmin - maximum sub-timestep allowed

physic_mod.F

  • nbapp_chem is not fixed anymore here but deftank/in physiq.def
  • Ambipolar diffusion activated if (ok_iondiff) is true

conf_phys.f90

  • add ok_iondiff as parameters to read from physiq.def file set to .false. by default
  • add nbapp_chem as parameters to characterize the Number of calls to the chemistry routines (per Venusian day) instead of be fixed in physic_mod.F

to read from physiq.def file set to 24000 by default

cleshphys.h

  • added ok_iondiff & nbapp_chem in COMMON/clesphys_l/
  • removed nbapp_chim

phytrac_chemistry.F

  • Added security in the calculation of the sza_local in order to avoid the rare case where the |range| is above 1
File size: 64.4 KB
Line 
1!       SUBROUTINE plasmaphys(q,rhoN,P,pk,teta,mmean,mudyn,
2!    & tempk,tetak,tempe,tetae,wcontk,wcovk,Efieldz,wphysk,alt)
3     
4       SUBROUTINE iondiff_red(ngrid,nlayer,nqmx,pplay,pplev,
5!     &                        ptneutre,ptelect,ption,pq,
6     &                        ptneutre,pq,pphis,
7     &                        gmtime,lat,lon,
8     &                        ptimestep,pdteuv,pdtconduc,pdqdiff,
9     &                        pdqiondiff)
10     
11!****************************************************************
12!
13!     Ion ambipolar diffusion routine
14!
15!     Authors: J-Y. Chaufray
16!     -------
17!
18!     J-Y. Chaufray
19!
20!     Vitesse verticale et ses effets sur nk, Tk, uk, vk, etc...
21!     Seule les effets du terme (wk - wn) sont pris en compte ici
22!     Pour l instant seul les effets sur la densite nk sont pris en compte
23!
24!     Version: 14/11/2020
25!
26!     -------
27!
28!     2022/10/22: adapted and adding by Antoine Martinez
29!
30!*****************************************************************
31     
32      USE chemparam_mod
33      USE infotrac_phy, only: tname
34      !USE dimphy   
35     
36      use plasmaphys_venus_mod
37      use iono_h, only: temp_elect, temp_ion
38     
39      IMPLICIT NONE
40     
41! Solve the vertical dynamics equation
42! The equation is solved on a altitude refined grid after interpolation
43! final results are reinterpolated on the atmospheric pressure grids
44! Upper boundary
45! The velocity at upper altitude should be parametrized (escape)
46! Here I assume w2up (nlayer)=0
47! Lower boundaru
48! W2(1) = O (The ion velocity is equal to neutral velocity)
49     
50!#include "dimensions.h"
51!#include "paramet.h"
52!#include "plasmadyn.h"
53!#include "comgeom.h"
54     
55! =========================
56! ==== INPUT / OUTPUT =====
57! =========================
58     
59! ====== INPUT ======
60     
61      INTEGER,intent(in) :: ngrid  ! number of atmospheric columns
62      INTEGER,intent(in) :: nlayer ! number of atmospheric layers
63      INTEGER,intent(in) :: nqmx   ! number of advected tracers
64      REAL,intent(in) :: ptimestep ! physical timestep (s)
65      REAL,intent(in) :: pphis(ngrid) ! geopotential of the surface (m2/s)
66     
67      REAL,intent(in) :: pplay(ngrid,nlayer)        ! mid-layer pressure (Pa)
68      REAL,intent(in) :: pplev(ngrid,nlayer+1)      ! inter-layer pressure (Pa)
69      REAL,intent(in) :: pq(ngrid,nlayer,nqmx)      ! mid-layer mass mixing ratio tracers (kg/kg)
70     
71      REAL,intent(in) :: ptneutre(ngrid,nlayer)     ! mid-layer NEUTRAL  temperature (K)
72!      REAL pdt(ngrid,nlayer)
73!      REAL,intent(in) :: ptelect(ngrid,nlayer)      ! mid-layer ELECTRON temperature (K)
74!      REAL,intent(in) :: ption(ngrid,nlayer)        ! mid-layer ION      temperature (K)
75     
76      REAL,intent(in) :: pdteuv(ngrid,nlayer)       ! Temperature EUV heating tendancy (K/s)
77      REAL,intent(in) :: pdtconduc(ngrid,nlayer)    ! Temperature conduction tendancy  (K/s)
78      REAL,intent(in) :: pdqdiff(ngrid,nlayer,nqmx) ! neutral mmr molecular diffusion tendancy (kg/kg/s)
79!      real pdq(ngrid,nlayer,nq)
80     
81      REAL,intent(in) :: gmtime     ! day fraction   ratio [from 0 to 1]
82      REAL,intent(in) :: lat(ngrid) ! grid latitude  [degree]
83      REAL,intent(in) :: lon(ngrid) ! grid longitude [degree]
84     
85! ====== OUTPUT ======
86     
87      REAL            :: pdqiondiff(ngrid,nlayer,nqmx) ! ion mmr tendancy (kg/kg/s)
88     
89! =========================
90! ======== LOCAL ==========
91! =========================
92     
93      integer, parameter :: t_elec_origin= 2
94      integer, parameter :: t_ion_origin = 1
95      !Electronic temperature. Index 1 -> Theis et al. 1980 - model data ; Index 2-> Theis et al. 1984 - model data
96      !Ion        temperature. Index 1 -> from VIRA Model based on the model of Miller et al. 1980 & 1984.
97     
98      REAL :: lon_sun, sza_local,szad_local,mu_local
99      REAL :: lon_local(ngrid), lat_local(ngrid)
100     
101      REAL :: ptelect(ngrid,nlayer)      ! mid-layer ELECTRON temperature (K)
102      REAL :: ption(ngrid,nlayer)        ! mid-layer ION      temperature (K)
103     
104      REAL :: qnew(ngrid,nlayer,nqmx), qold(ngrid,nlayer,nqmx)
105     
106      REAL :: mmean_new(ngrid,nlayer), vmr_new(ngrid,nlayer,nqmx)
107      REAL :: vmr_ion(ngrid,nlayer)
108     
109      !REAL wcontk(ip1jmp1,llp+1,ndiffions),Efieldz(ip1jmp1,llp+1)     
110      !REAL wcovk(ip1jmp1,llp+1,ndiffions),wphysk(ip1jmp1,llp+1,ndiffions)
111     
112! 1D field I pass to double precision
113      REAL(kind=kind(1.D0)) tdiff1,tdiff2
114      REAL(kind=kind(1.D0)),dimension(nlayer+1) :: P1D,Ez1d
115      REAL(kind=kind(1.D0)),dimension(nlayer) :: Mmean1D,RhoN1D,qe1D
116      REAL(kind=kind(1.D0)),dimension(nlayer) :: Pnlay1D,ZZ1D,TempN1d
117      REAL(kind=kind(1.D0)),dimension(nlayer) :: TempE1D,TempE1dint
118      REAL(kind=kind(1.D0)),dimension(nlayer) :: FacTe,TempE1dnew
119     
120      REAL(kind=kind(1.D0)),dimension(naltvert) :: Zraf,Praf,Tnraf,Mraf
121      REAL(kind=kind(1.D0)),dimension(naltvert) :: Rraf,QEraf,REraf
122      REAL(kind=kind(1.D0)),dimension(naltvert) :: RErafold,TEraf,Ez
123     
124      !REAL(kind=kind(1.D0)) X(lld),Y(lld),a(lld),b(lld),c(lld)
125      !REAL(kind=kind(1.D0)) X2(lld),Y2(lld),a2(lld),b2(lld),c2(lld)
126     
127      INTEGER,dimension(1) :: indZ
128     
129      REAL(kind=kind(1.D0)) :: Rho0,Nu0,T0,H0,D0,dt,dtad,time0
130     
131      REAL(kind=kind(1.D0)) :: tdiff3,tmin,Wmax,ttot,tsim,tdiff
132     
133      REAL(kind=kind(1.D0)) :: pphis1D    ! geopotential of the surface (m2/s)
134     
135      REAL(kind=kind(1.D0)),save :: Wlim
136     
137      REAL(kind=kind(1.D0)),dimension(naltvert) :: RIad,Tnad,TIad,Tead
138      REAL(kind=kind(1.D0)),dimension(naltvert) :: RElad,Dad,Zad
139      REAL(kind=kind(1.D0)) :: Dzraf,Dz,FacEsc,Time
140      REAL(kind=kind(1.D0)),dimension(naltvert) :: alpha,beta,delta
141      REAL(kind=kind(1.D0)),dimension(naltvert) :: eps,gama,dzeta,eta
142      REAL(kind=kind(1.D0)),dimension(naltvert) :: Atri,Btri,Ctri
143      REAL(kind=kind(1.D0)),dimension(naltvert) :: Dtri,Xtri,Ytri
144      INTEGER :: ig,k,l,ln,kn,ij0,l0,ln0,iter,ld,iq,istep,ke,niter,nn
145      INTEGER :: iqelec
146      INTEGER :: nsubreal
147      LOGICAL,SAVE :: firstcall=.true.
148     
149      REAL(kind=kind(1.D0)) masse,invsgmu,Konst
150      REAL(kind=kind(1.D0)) masseU,kBolt,g,RgazP,Rvenus,Grav,Mvenus,PI
151     
152      integer,save :: ntracers
153      integer,dimension(:),allocatable,save :: gcmind ! array of GCM indexes
154      real,dimension(:),allocatable,save :: mol_tr    ! array of mol mass traceurs
155      real,dimension(:),allocatable,save :: Charge_tr ! array of electric charge traceurs
156     
157      integer,dimension(nqmx) :: indic_tracers,indic_iondiff
158      integer,parameter       :: nb_espIon_diff=24
159      character(len=20),dimension(nb_espIon_diff) :: ListeIonDiff
160      character(len=20)                           :: electrans(1)
161     
162      integer,dimension(:),allocatable,save :: itrans    ! array of sub index gcmind of fluid ions
163      integer,dimension(1),save             :: etrans    ! index of gcmind electron   
164     
165      INTEGER,save :: iz_vertplasma ! vertical index above which ion diffusion is computed
166      INTEGER,save :: ndiffions     ! nombre d'ions decrit
167     
168      ! ---------- ALLOCATABLE ARRAY ------------------------------
169     
170      ! (ngrid,nlayer,ntracers)
171      !REAL(kind=kind(1.D0)),dimension(:,:),allocatable :: qnew
172      !REAL(kind=kind(1.D0)),dimension(:,:),allocatable :: qold
173     
174      ! (nlayer,ntracers)
175      REAL(kind=kind(1.D0)),dimension(:,:),allocatable,save :: q1D
176      REAL(kind=kind(1.D0)),dimension(:,:),allocatable,save :: RhoK1D
177     
178      ! (nlayer,ndiffions)           
179      REAL(kind=kind(1.D0)),dimension(:,:),allocatable,save :: qk1d
180      REAL(kind=kind(1.D0)),dimension(:,:),allocatable,save :: qk1dnew
181      REAL(kind=kind(1.D0)),dimension(:,:),allocatable,save :: qk1dint
182      REAL(kind=kind(1.D0)),dimension(:,:),allocatable,save :: TempI1D
183      REAL(kind=kind(1.D0)),dimension(:,:),allocatable,save :: TempI1din
184      REAL(kind=kind(1.D0)),dimension(:,:),allocatable,save :: TempI1dnw
185      REAL(kind=kind(1.D0)),dimension(:,:),allocatable,save :: FacQ
186      REAL(kind=kind(1.D0)),dimension(:,:),allocatable,save :: FacTI
187     
188      ! (nlayer+1,ndiffions)
189      REAL(kind=kind(1.D0)),dimension(:,:),allocatable,save :: Wk1d
190      REAL(kind=kind(1.D0)),dimension(:,:),allocatable,save :: Wk1d2
191      REAL(kind=kind(1.D0)),dimension(:,:),allocatable,save :: Wk1d3
192     
193      ! (ngrid,ndiffions)
194      REAL(kind=kind(1.D0)),dimension(:,:),allocatable,save :: ueff
195     
196      ! (ndiffions)
197      REAL(kind=kind(1.D0)),dimension(:),allocatable,save :: Mtot
198      REAL(kind=kind(1.D0)),dimension(:),allocatable,save :: Mtot2
199     
200      ! ---------- ALLOCATABLE ARRAY ------------------------------
201     
202      ! (naltvert,ntracers)
203      REAL(kind=kind(1.D0)),dimension(:,:),allocatable,save :: Qraf
204      REAL(kind=kind(1.D0)),dimension(:,:),allocatable,save :: Rkraf
205     
206      ! (naltvert,ndiffions)
207      REAL(kind=kind(1.D0)),dimension(:,:),allocatable,save :: QIraf
208      REAL(kind=kind(1.D0)),dimension(:,:),allocatable,save :: RIraf
209      REAL(kind=kind(1.D0)),dimension(:,:),allocatable,save :: RIrafold
210      REAL(kind=kind(1.D0)),dimension(:,:),allocatable,save :: TIraf
211      REAL(kind=kind(1.D0)),dimension(:,:),allocatable,save :: FIraf
212      REAL(kind=kind(1.D0)),dimension(:,:),allocatable,save :: DIraf
213      REAL(kind=kind(1.D0)),dimension(:,:),allocatable,save :: Rhock
214     
215      ! (ndiffions)
216      REAL(kind=kind(1.D0)),dimension(:),allocatable,save :: MtotZ1
217      REAL(kind=kind(1.D0)),dimension(:),allocatable,save :: MtotZ2
218     
219! ==========================================================
220! ==========================================================
221! ========== FIRST CALL AND MAIN PARAMETERS ================
222! ==========================================================
223! ==========================================================
224     
225! ======= Initializations at first call ==========
226     
227      if (firstcall) then
228       
229      ! === On elimine la microphysique ==========
230       
231        ntracers=0
232        indic_tracers(1:nqmx)=0
233       
234        ! type_tr ==> 1: neutral gas, 2: ion gas, 3: liquid, 10: others
235        do iq=1,nqmx
236          if ((type_tr(iq) .le. 2.) .and. (type_tr(iq) .gt. 0.)) then
237            indic_tracers(iq)=1
238            ntracers=ntracers+1
239          endif
240        enddo
241        allocate(gcmind(ntracers))
242        allocate(mol_tr(ntracers))
243        allocate(Charge_tr(ntracers))
244       
245        ! Store gcm indexes in gcmind and molar masses in mol_tr
246        nn=0
247        do iq=1,nqmx
248          if (indic_tracers(iq) .eq. 1) then
249            nn=nn+1
250            gcmind(nn) = iq
251            mol_tr(nn) = M_tr(iq)
252            Charge_tr(nn) = 0.D0
253          endif
254        enddo
255       
256      ! === On selectionne les especes ioniques ==========
257       
258        ndiffions =  0
259        indic_iondiff(1:nqmx) = 0
260       
261        ! define used ions
262        ListeIonDiff(1)="o2plus"
263        ListeIonDiff(2)="oplus"
264        ListeIonDiff(3)="cplus"
265        ListeIonDiff(4)="nplus"
266        ListeIonDiff(5)="co2plus"
267        ListeIonDiff(6)="noplus"
268        ListeIonDiff(7)="coplus"
269        ListeIonDiff(8)="hplus"
270        ListeIonDiff(9)="h2oplus"
271        ListeIonDiff(10)="h3oplus"
272        ListeIonDiff(11)="hcoplus"
273        ListeIonDiff(12)="n2plus"
274        ListeIonDiff(13)="ohplus"
275       
276        ! define electron   
277        electrans(1)="elec"     
278       
279        ! seek the index and number of ion/electron
280        do iq=1,ntracers
281          do nn=1,nb_espIon_diff
282            if (ListeIonDiff(nn) .eq. tname(gcmind(iq))) then
283              indic_iondiff(iq)=iq ! 1
284            endif
285            ! we don't diffuse electron but we readapt electron density at the end
286            if (electrans(1) .eq. tname(gcmind(iq))) then
287              etrans(1)=iq
288            endif
289          enddo
290          if (indic_iondiff(iq) .ge. 1) then
291            print*,'Ion specie ', tname(gcmind(iq)),
292     &               'diffused in Iondiff_red'
293            ndiffions=ndiffions+1
294          endif
295        enddo
296        print*,'number of diffused ion:',ndiffions
297        allocate(itrans(ndiffions))
298       
299        !   Store gcm ion indexes in itrans
300        nn=0
301        do iq=1,ntracers
302          if (indic_iondiff(iq) .ge. 1) then
303            nn=nn+1
304            itrans(nn)=indic_iondiff(iq) ! iq
305            Charge_tr(itrans(nn))=1.!qcharge
306          endif
307        enddo
308        Charge_tr(etrans(1))=-1.!qcharge*(-1.)
309       
310        ! find vertical index above which ion diffusion is computed
311        do l=1,nlayer
312          if (pplay(1,l) .gt. Pdiffion) then
313            iz_vertplasma = l
314          endif
315        enddo
316        iz_vertplasma = iz_vertplasma+1 ! seuil de la diffusion verticale
317        print*,'vertical index for ion diffusion',
318     &         iz_vertplasma,pplay(1,iz_vertplasma)
319        ! iz_plasma     = iz_vertplasma-3  ! seuil de la dynamique horizontale
320       
321      ! allocatation des tableaux dependants du nombre d especes diffusees
322        allocate(ueff(ngrid,ndiffions))
323       
324        ! ----------------------------------- !
325       
326        !allocate-(qnew(ngrid,nlayer,ntracers))
327        !allocate-(qold(ngrid,nlayer,ntracers))
328       
329        allocate(q1D(nlayer,ntracers))
330        allocate(RhoK1D(nlayer,ntracers))
331       
332        allocate(qk1d(nlayer,ndiffions))
333        allocate(qk1dnew(nlayer,ndiffions))
334        allocate(qk1dint(nlayer,ndiffions))
335        allocate(TempI1D(nlayer,ndiffions))
336        allocate(TempI1din(nlayer,ndiffions))
337        allocate(TempI1dnw(nlayer,ndiffions))
338        allocate(FacQ(nlayer,ndiffions))
339        allocate(FacTI(nlayer,ndiffions))
340       
341        allocate(Wk1d(nlayer+1,ndiffions))
342        allocate(Wk1d2(nlayer+1,ndiffions))
343        allocate(Wk1d3(nlayer+1,ndiffions))
344       
345        allocate(Mtot(ndiffions))
346        allocate(Mtot2(ndiffions))
347       
348        ! ----------------------------------- !
349        allocate(Qraf(naltvert,ntracers))
350        allocate(Rkraf(naltvert,ntracers))
351       
352        allocate(QIraf(naltvert,ndiffions))
353        allocate(RIraf(naltvert,ndiffions))
354        allocate(RIrafold(naltvert,ndiffions))
355        allocate(TIraf(naltvert,ndiffions))
356        allocate(FIraf(naltvert,ndiffions))
357        allocate(DIraf(naltvert,ndiffions))
358        allocate(Rhock(naltvert,ndiffions))
359       
360        allocate(MtotZ1(ndiffions))
361        allocate(MtotZ2(ndiffions))
362       
363        ! ----------------------------------- !
364       
365      ! Conditions aux limites
366       
367        Wlim=500. ! Maximal absolute value of vertical velocity (m/s)  (to avoid unstabilities in wdu/dz term)
368        ueff(1:ngrid,1:ndiffions)=0. ! Effusion velocity at top in m/s
369        !ueff(1:ngrid,1:ndiffions)=200.
370       
371        firstcall=.FALSE.
372      endif ! firstcall
373     
374      masseU = 1.660538782d-27
375      kBolt  = 1.3806504d-23
376      Konst  = kBolt/masseU
377      RgazP  = 8.314472
378      PI     = 2.*ASIN(1.) ! 3.141592653
379      g      = 8.87d0
380      Rvenus = 6051800d0 ! Used to compute escape parameter no need to be more accurate
381      Grav   = 6.67d-11
382      Mvenus = 4.86d24
383      ij0    = 6000 ! For test
384      invsgmu= 1d0/g/masseU
385     
386      ke=etrans(1)
387     
388      alpha(1:naltvert)=0D0
389      beta(1:naltvert)=0D0
390      delta(1:naltvert)=0D0
391      eps(1:naltvert)=0D0
392      dzeta(1:naltvert)=0D0
393      eta(1:naltvert)=0D0
394!! Open Vertical file to check it
395! if (firstcall) then
396! open(unit=301,file='Plasma_vertical.dat',status='unknown')
397! firstcall=.FALSE.
398! endif
399!! First we derive only parameters useful for diffusion
400!! Extrapolation of paramaters to all altitudes (we solve the dynamics equation on larger altitude range)
401!!
402!!   ===> ON A PAS BESOIN DE CETTE PARTIE CAR ON EST DEJA DANS LA PHYSIQUE
403!!       CALL DYN2DIFFK(Pk,teta,P,mmean,rhoN,Preff,cpp,kappa,
404!!     &  q,Mtraceur,Tempk,tetak,TempE,TetaE,Preffplasma,itrans,etrans,
405!!     &  llm,nqtot,llp,ndiffions,ip1jmp1,iz_plasma,kappak,kappae,
406!!     &  PNlay,TempN,Rhok,TempI,TetaI,rhotetaI,TempF,TetaF,rhotetaF)
407!!
408     
409! ==========================================================
410! ==========================================================
411! ========== Main Loop on horizontal grids =================
412! ==========================================================
413! ==========================================================
414     
415! ========== Initialization =================
416     
417      dt=dble(ptimestep)      !dtplasma*iappvert
418      nsubreal=nsubvert
419      niter=2
420      qold(1:ngrid,1:nlayer,1:nqmx)=pq(1:ngrid,1:nlayer,1:nqmx)
421      qnew(1:ngrid,1:nlayer,1:nqmx)=pq(1:ngrid,1:nlayer,1:nqmx)
422!      tdiff=dt/nsubreal
423!      print*,'dt',dt
424     
425      lon_sun = (0.5 - gmtime)*2.*PI
426      lon_local(:) = lon(:)*PI/180.
427      lat_local(:) = lat(:)*PI/180.
428     
429! ========== Horizontal loop start ==========
430     
431      DO ig=1,ngrid
432       
433        sza_local  =  cos(lat_local(ig))*cos(lon_local(ig))
434     &               *cos(lon_sun) + cos(lat_local(ig))
435     &               *sin(lon_local(ig))*sin(lon_sun)
436        sza_local  = min(sza_local,1.)
437        sza_local  = max(-1.,sza_local)
438        sza_local  = acos(sza_local)                   ! RAD
439        szad_local = sza_local * 180. / PI             ! Degree
440        mu_local   = cos(sza_local)
441       
442        ttot   = dt
443        tsim   = 0.
444        tdiff1 = dt/nsubvert !
445        tdiff2 = dt/nsubmin  !
446        !if (nsubreal .ne. nsubvert) print*,'ig',ig,tdiff
447       
448        ZZ1D(1:nlayer)     = 0D0
449        P1D(1:nlayer)      = 0D0
450        Pnlay1D(1:nlayer)  = 0D0
451        TempN1D(1:nlayer)  = 0D0
452        Mmean1D(1:nlayer)  = 0D0
453        RHON1D(1:nlayer)   = 0D0
454        RHOK1D(1:nlayer,1:ntracers)=0D0
455        Q1D(1:nlayer,1:ntracers)=0D0
456        TEMPI1D(1:nlayer,1:ndiffions)=0D0
457        !TETAI1D(1:nlayer,1:ndiffions)=0D0
458        TEMPE1D(1:nlayer)=0D0
459        !TETAE1D(1:nlayer)=0D0
460        Zraf(1:naltvert)=0D0
461        Praf(1:naltvert)=0D0
462        Tnraf(1:naltvert)=0D0
463        Mraf(1:naltvert)=0D0
464        Qraf(1:naltvert,1:ntracers)=0D0
465        Rkraf(1:naltvert,1:ntracers)=0D0
466        REraf(1:naltvert)=0D0
467        QIraf(1:naltvert,1:ndiffions)=0D0
468        TIraf(1:naltvert,1:ndiffions)=0D0
469        FIraf(1:naltvert,1:ndiffions)=0D0
470        QEraf(1:naltvert)=0D0
471        TEraf(1:naltvert)=0D0
472        rhock(1:naltvert,1:ndiffions)=0D0
473       
474! print*,'ig',ig,llm
475! First we compute 1d fields in double precision
476! The mixing ratios slightly change here (< 0.1%)
477       
478        pphis1D       = pphis(ig)
479        P1D(nlayer+1) = pplev(ig,nlayer+1)
480        do l=1,nlayer
481          P1D(l)      = pplev(ig,l)     
482          Pnlay1D(l)  = pplay(ig,l)
483          TempN1D(l)  = ptneutre(ig,l)
484          do iq=1,ntracers
485            q1d(l,iq)=qold(ig,l,gcmind(iq))
486          enddo
487        ! compute the mean molecular mass
488          Mmean1d(l)  = (1D0/sum(q1d(l,:)/dble(mol_tr(:))))
489         
490        ! Compute the total mass density (kg/m3)
491          rhok1d(l,:) = 0D0
492          RhoN1D(l)   = Pnlay1D(l)*Mmean1d(l)/TempN1D(l)/Konst
493          do iq=1,ntracers
494            rhok1d(l,iq)=q1d(l,iq)*RhoN1D(l)
495          enddo
496         
497          do k=1,ndiffions
498            kn=itrans(k)
499            qk1d(l,k)=q1D(l,kn)
500            !TempI1D(l,k)=ption(ig,l)
501            !TetaI1D(l,k)=TetaI(ig,l,k)
502            !rhotetaI1D(l,k)=rhotetaI(ig,l,k)
503          enddo
504         
505          !ke=etrans(1)
506          qe1D(l)=q1D(l,ke)
507          !TempE1D(l)=ptelect(ig,l)
508          !TetaE1D(l)=TetaF(ig,l)
509          !rhoTetaE1D(l)=rhoTetaf(ig,l)
510        enddo ! l (vertical levels)
511       
512! Compute colonne mass of each dynamic ion
513       
514        Mtot(1:ndiffions)=0D0
515       
516!        do k=1,ndiffions
517!        kn=itrans(k)
518!        do l=1,lld
519!        ln=l+iz_vertplasma
520!        print*,ln,llm
521!        Mtot(k)=Mtot(k)+qk1d(ln,k)*rhoN1d(ln)
522!        enddo
523!        enddo
524       
525!        print*,'Mtot before',Mtot
526       
527! Compute altitude of the scalar grid pressure level
528       
529        !if (isnan(TempN1d(1))) then
530        !  write (*,*) 'PROBLEME NAN: TITI'
531        !endif
532       
533        CALL ZVERTK(Pnlay1D,TempN1D,mmean1d,ZZ1d,nlayer,pphis1D)
534       
535! Compute lectron and Ion temperature of the scalar grid pressure level
536       
537        do l=1,nlayer
538          ptelect(ig,l) = temp_elect(ZZ1d(l)/1000.,TempN1d(l),
539     &                               szad_local,t_elec_origin) 
540          ption(ig,l)   = temp_ion(ZZ1d(l)/1000.,TempN1d(l),
541     &                             szad_local,t_ion_origin) 
542         
543          TempE1D(l) = ptelect(ig,l)
544          !!! The ion temperature is fixed as the half of the electron temperature
545          !!! calculated in ion_h.F for the stability of the program and because the
546          !!! ion temperature is lower than electron temperature.
547          do k=1,ndiffions
548            !kn=itrans(k)
549            TempI1D(l,k)=max(0.5*ptelect(ig,l),TempN1d(l))
550          enddo 
551        enddo ! l (vertical levels)       
552        !write(*,*) szad_local,ZZ1d(:)/1000.,ptelect(ig,:), ption(ig,:)
553!        do l=1,llm
554!          print*,'Pnlay1D, Plev',Pnlay1D(l),P1D(l),rhok1D(l,1:ntracers)
555!        enddo
556! I use a better vertial resolution above iz_vertplasma and altitude grid interpolation
557       
558        indZ(1)=0
559        CALL UPPER_RESOLK(ZZ1D,Pnlay1D,TempN1d,mmean1d,rhok1d,TempI1d,
560     & TempE1D,mol_tr,nlayer,ndiffions,ntracers,naltvert,iz_vertplasma,
561     & itrans,etrans,Zraf,Tnraf,Praf,Mraf,Qraf,Rkraf,
562     & QIraf,TIraf,FIraf,Qeraf,Teraf,indZ,gcmind)
563       
564!       print*,'z last level, P last level',Zraf(:),Praf(:)
565       
566        DZraf=Zraf(2)-Zraf(1)
567       
568        do l=1,naltvert
569          do k=1,ndiffions
570            kn=itrans(k)
571            RIraf(l,k)=Rkraf(l,kn)
572          enddo
573          !ke=etrans(1)
574          REraf(l)=Rkraf(l,ke)
575        enddo
576       
577! Test
578!       do l=1,naltvert
579!        print*,'l',l,Zraf(l),Praf(l),Tnraf(l),Qiraf(l,1),Tiraf(l,1),
580!     &  Qeraf(l),Teraf(l),Riraf(l,1),Reraf(l)
581!        enddo
582       
583       
584! Before solving the diffusion equation, we reddo interpolation on GCM grid
585! to derive mass correction factor due to interpolation (for mass conservation purpose)
586       
587        CALL GCMGRID_PK(Zraf,Praf,TNraf,Mraf,
588     &  QIraf,RIraf,TIraf,Teraf,
589     &  P1D,Pnlay1D,TempN1d,Mmean1d,
590     &  qk1d,qk1dint,TempI1D,TempI1Din,TempE1d,TempE1dint,
591     &  naltvert,ndiffions,nlayer,iz_vertplasma,itrans,etrans)
592
593        !DO ln=1,nlayer
594          !IF (k.eq.1) THEN
595        !    write(*,*) real(qk1Dint(ln,:)), qold(ig,ln,gcmind(ke))
596          !ENDIF
597        !ENDDO
598       
599        FacQ(1:nlayer,1:ndiffions)=1D0
600        FacTI(1:nlayer,1:ndiffions)=1D0
601        FacTE(1:nlayer)=1D0
602       
603        CALL CORRFACK(qk1d,qk1dint,TempI1d,TempI1din,TempE1d,
604     &  TempE1dint,FacQ,FacTI,FacTE,nlayer,ndiffions)
605       
606        CALL DCOEFFK(TIraf,FIraf,mol_tr,DIraf,
607     &  naltvert,ndiffions,ntracers,itrans)
608       
609        Rho0=1d-20
610       
611        T0=Tnraf(naltvert)!TIraf(naltvert,1)
612       
613        do ln=1,naltvert
614          Tnad(ln)=TNraf(ln)/T0
615        enddo
616       
617!        if(abs(lon(ig)-(48.75)).le.0.1) then
618!          if(abs(lat(ig)-(-80.625)).le.0.1) then
619!            write(*,*) 'ZZ1d= ', ZZ1d(iz_vertplasma:90)
620!            write(*,*) 'Zraf= ', Zraf(1:naltvert)
621!            write(*,*) 'Pnlay1D= ', Pnlay1D(iz_vertplasma:90)
622!            write(*,*) 'Praf= ', Praf(1:naltvert)
623!            DO k=1,ndiffions
624!            kn=itrans(k)
625!            if((tname(gcmind(kn)).eq.'hplus').or.
626!     &         (tname(gcmind(kn)).eq.'1oplus').or.
627!     &         (tname(gcmind(kn)).eq.'1noplus')) then
628!             write(*,*) tname(gcmind(kn))//' DEBUT'
629!             write(*,*) 'avant= ',
630!     &       qk1d(iz_vertplasma:90,k)
631!             write(*,*) 'apres= ',
632!     &       qk1dint(iz_vertplasma:90,k)
633!             write(*,*) 'interm= ',
634!     &       Rkraf(1:naltvert,kn)
635!             write(*,*) 'corrmass= ',
636!     &       real(FacQ(iz_vertplasma:90,k))
637!            endif
638!            ENDDO
639!
640!          endif
641!        endif       
642
643! do l=1,llm
644!        print*,'effect interp',l,qk1d(l,2),qk1dint(l,2)*facQ(l,2),
645!     & TempI1d(l,2),TempI1din(l,2)*facTI(l,2),TempE1D(l),
646!     & TempE1Dint(l)*facTE(l),facQ(l,:),FacTI(l,:),FacTe(l)
647! enddo
648! stop
649       
650! Compute vertical velocity to derive subtime step
651       
652        Tdiff=1.E-9
653        DO k=1,ndiffions
654          kn=itrans(k)
655          H0=kbolt*T0/mol_tr(kn)/g/masseU
656          D0=DIraf(naltvert,k)
657          Time0=H0*H0/D0
658          Time=Tdiff/Time0
659          do ln=1,naltvert
660            RIad(ln)=RIraf(ln,k)/Rho0
661            RElad(ln)=Reraf(ln)/rho0
662            TIad(ln)=TIraf(ln,k)/T0
663            TEad(ln)=TEraf(ln)/T0
664            Zad(ln)=Zraf(ln)/H0
665            Dad(ln)=DIraf(ln,k)/D0
666          enddo
667         
668          DZ=Zad(2)-Zad(1)
669         
670          CALL DIFFPARAMK(Tiad,Tead,RElad,Dad,DZ,
671     &  alpha,beta,gama,delta,eps,dzeta,eta,naltvert,Wmax)
672         
673! if (ig .eq. ij0) then
674! print*,'zeta',dzeta(:)
675! print*,'eta',eta(:)
676! endif
677         
678          FacEsc=exp(-ueff(ig,k)*Time0/H0*DZ/Dad(naltvert-1))
679         
680          CALL MATCOEFFK(alpha,beta,gama,delta,eps,dzeta,eta,Dad,Riad,
681     &  FacEsc,dZ,time,Atri,Btri,Ctri,Dtri,Tiad,Tead,naltvert)
682         
683          rhock(1,k)=0D0
684          rhock(2,k)=0D0
685          do ln=3,naltvert-1
686     
687            if (ln .lt. naltvert-1) then
688              rhock(ln,k)=-Dad(ln)*(-RIraf(ln+2,k)+8.*RIraf(ln+1,k)
689     &  -8.*RIraf(ln-1,k)+RIraf(ln-2,k))/RIraf(ln,k)/12D0/DZ
690     &  -Dad(ln)*(alpha(ln)+beta(ln)+dzeta(ln))
691            endif
692           
693            if (ln .eq. naltvert-1) then
694              rhock(ln,k)=-Dad(ln)*
695     &  (log(RIraf(ln+1,k))-log(RIraf(ln,k)))/Dz
696     &  -Dad(ln)*(alpha(ln)+beta(ln)+dzeta(ln))
697            endif
698           
699            rhock(ln,k)=D0/H0*rhock(ln,k)
700          enddo
701         
702          rhock(naltvert,k)=ueff(ig,k)
703        enddo
704       
705       
706        Wmax=maxval(abs(rhock))
707        if (mu_local .ge. 0.3) tdiff3=h0*Dz/Wmax/300d0 !20d0 !/100d0
708        if (mu_local .le. -0.3) tdiff3=h0*dz/Wmax/300d0 !300d0
709        if (abs(mu_local) .lt. 0.3) tdiff3=h0*dz/Wmax/400d0 !/1000d0
710       
711        tmin=max(tdiff1,tdiff3)
712        tmin=min(tmin,tdiff2)
713         
714        nsubreal=anint(dt/tmin)
715               
716        tdiff=tmin
717
718        !! This is a weird factor to optimize time calculation
719        !! Here, it is to be sure that the first dtime is very low
720        tmin = 2*0.45*(DZ*H0)**2.
721        tmin = tmin/(maxval(abs(DIraf(:,:)))+(DZ*H0)*Wmax)
722        if (tdiff .ge. tmin) then
723          tdiff = tmin
724        endif
725!        if (ig .eq. ij0) then
726!          print*,'nreal',nsubreal,nsubvert,nsubmin,
727!       &  Wmax,tmin,tdiff,tdiff1,tdiff2,tdiff3,h0*dZ,tsim,ttot,h0,dZ
728!          write(301,*) nsubreal
729!        endif
730       
731        !tdiff=8.75D-5 !tmin
732        tdiff=tmin
733       
734! Fix a minimal value for tdiff
735       
736! ========== TEMPORAL Loop =================
737       
738! Loop over dynamic ions to solve the diffusion equation
739       
740! Compute total mass on the Z grid before diffusion
741       
742!        MtotZ1(:)=0D0
743!        do k=1,ndiffions
744!        do l=1,naltvert
745!        MtotZ1(k)=MtotZ1(k)+RIraf(l,k)*Dzraf
746!        enddo
747!        enddo
748!        print*,'MtotZ',MtotZ1
749       
750        !if (szad_local.ge.93.) tsim = ttot
751!    DO istep=1,nsubreal
752        DO WHILE (tsim .lt. ttot)
753! if (ig .eq. ij0) then
754! print*,'tsim',tsim
755!        do l=1,naltvert
756! print*,'density1',l,REraf(l)/Mtraceur(31)/masseU/1d6,
757!     &                   Riraf(l,2)/Mtraceur(22)/masseU/1d6
758! enddo
759! endif
760        RErafold=REraf
761        RIrafold=RIraf
762        DO iter=1,niter ! needed to use approximatively electron density (t+dt)
763          DO ln=1,naltvert
764            RElad(ln)=REraf(ln)/Rho0
765            REraf(ln)=RErafold(ln)   
766         
767          ! if (ig .eq. ij0) then
768          ! print*,'density',RElad(ln)/Mtraceur(31)/masseU/1d6*rho0
769          !     &                  ,Riraf(ln,2)/Mtraceur(22)/masseU/1d6
770          ! endif
771         
772          ENDDO ! ln
773         
774          DO k=1,ndiffions
775            kn=itrans(k)
776            H0=kbolt*T0/mol_tr(kn)/g/masseU
777            D0=DIraf(naltvert,k)
778            Time0=H0*H0/D0
779            Time=Tdiff/Time0
780           
781! Compute the diffusion coefficient using the collion frequency
782           
783!            print*,'now adimension'
784           
785! Adimension all parameters by value at iz_vertplasma+1
786           
787            do ln=1,naltvert
788              RIad(ln)=RIraf(ln,k)/Rho0
789!              RElad(ln)=REraf(ln)/Rho0
790              TIad(ln)=TIraf(ln,k)/T0
791              TEad(ln)=TEraf(ln)/T0
792              Zad(ln)=Zraf(ln)/H0
793              Dad(ln)=DIraf(ln,k)/D0
794            enddo ! ln
795           
796            DZ=Zad(2)-Zad(1)
797           
798!            Tead(:)=0D0
799           
800            CALL DIFFPARAMK(Tiad,Tead,RElad,Dad,DZ,
801     &  alpha,beta,gama,delta,eps,dzeta,eta,naltvert,Wmax)
802           
803!            if (ig .eq. ij0) then
804!        print*,'zeta',dzeta(:)
805!        print*,'eta',eta(:)
806!        endif
807           
808            FacEsc=exp(-ueff(ig,k)*Time0/H0*DZ/Dad(naltvert-1))
809!     & *Tiad(ln-1)/(Tiad(ln-1)+Tead(ln-1)))
810           
811!            eta(:)=0D0
812!            dzeta(:)=0D0
813           
814!            do l=1,naltvert
815!            print*,l,alpha(l),dzeta(l),eta(l),RElad(l),Riad(l)
816!            enddo
817!            stop
818           
819           
820            CALL MATCOEFFK(alpha,beta,gama,delta,eps,dzeta,eta,Dad,
821     &  Riad,FacEsc,dZ,time,Atri,Btri,Ctri,Dtri,Tiad,Tead,naltvert)
822           
823            Xtri(:)=0D0
824           
825            Xtri=Riad
826            Ytri=Xtri
827             
828! Solve system
829            !write(*,*) tsim, iter
830            CALL TridagDP(Atri,Btri,Ctri,Dtri,Xtri,naltvert)
831
832            !STOP
833! Check mass conservation
834           
835! CALL CheckmassK(RIad,Xtri,naltvert)
836           
837! if (ig .eq. 1653) then
838! do l=1,naltvert
839! print*,'lphysnew',l,rho0*Xtri(l),k
840! enddo
841! endif
842           
843            do l=1,naltvert
844              if (iter .eq. niter) RIraf(l,k)=Rho0*Xtri(l)
845              if (RIraf(l,k) .lt. 0. .or. isnan(RIraf(l,k))) then
846              ! print*,'phys q <0',ig,l,k,iter,istep,Xtri(l),Ytri(l),Riraf(l,k),
847              !     &  REraf(l),Tead(l)*T0,Tiad(l)*T0
848                RIraf(l,k)=1D-24
849              ! RIraf(l,k)=RIraf(l-1,k)*exp(-DZ*(alpha(l-1)+beta(l-1)+dzeta(l-1)))
850              ! nsubreal=nsubreal+50
851              ! if (nsubreal .gt. 2000) then
852              ! print*,'oula nsubreal',ig,nsubreal
853              ! endif
854              ! goto 1
855              endif
856            enddo ! l
857           
858! Update electron mass density for next iteration
859           
860            DO l=1,naltvert
861              REraf(l)=REraf(l)+0.5*Charge_tr(kn)*mol_tr(ke)/mol_tr(kn)
862     &                  * rho0*(Xtri(l)-Ytri(l))
863            ENDDO ! l
864           
865!      if (ig .eq. ij0 .and. k .eq. 2) then
866!       do l=1,naltvert
867!         print*,'l1',l,Reraf(l)/mol_tr(ke)/masseU/1d6,
868!     &        Xtri(l)*rho0/mol_tr(kn)/masseU/1d6,
869!     &        Ytri(l)*rho0/mol_tr(kn)/masseU/1d6
870!       enddo ! l
871!      endif
872           
873!       Compute vertical velocity in atmospheric frame
874           
875            IF (iter .eq. niter) THEN
876              rhock(1,k)=0D0
877              rhock(2,k)=0D0
878              do ln=3,naltvert-1
879!       rhock(ln,k)=-Dad(ln)*((RIraf(ln+1,k)-RIraf(ln-1,k))/2D0/Dz
880!     & +RIraf(ln,k)*(alpha(ln)+beta(ln)+dzeta(ln)))
881!       rhock(ln,k)=D0/H0*rhock(ln,k)/RIraf(ln,k)
882             
883             
884! rhock(ln,k)=-Dad(ln)*
885!     &  ((log(RIraf(ln+1,k))-log(RIraf(ln,k)))/Dz
886!     &   +(alpha(ln)+beta(ln))*(Tiad(ln)/(Tead(ln)+Tiad(ln))))
887!       rhock(ln,k)=D0/H0*rhock(ln,k)
888             
889              if (ln .lt. naltvert-1) then
890! rhock(ln,k)=-Dad(ln)*(-log(RIraf(ln+2,k))+8.*log(RIraf(ln+1,k))
891!     &  -8.*log(RIraf(ln-1,k))+log(RIraf(ln-2,k)))
892!     &  /12D0/Dz
893                rhock(ln,k)=-Dad(ln)*(-RIraf(ln+2,k)+8.*RIraf(ln+1,k)
894     &        -8.*RIraf(ln-1,k)+RIraf(ln-2,k))/RIraf(ln,k)/12D0/DZ
895     &        -Dad(ln)*(alpha(ln)+beta(ln)+dzeta(ln))
896              endif
897             
898              if (ln .eq. naltvert-1) then
899! rhock(ln,k)=-Dad(ln)*((Tiad(ln)+Tead(ln))/Tiad(ln)*
900!     &  (log(RIraf(ln+1,k))-log(RIraf(ln,k)))/Dz)
901!     &  -Dad(ln)*(alpha(ln)+beta(ln))
902             
903                rhock(ln,k)=-Dad(ln)*
904     &        (log(RIraf(ln+1,k))-log(RIraf(ln,k)))/Dz
905     &        -Dad(ln)*(alpha(ln)+beta(ln)+dzeta(ln))
906             
907              endif
908             
909              rhock(ln,k)=D0/H0*rhock(ln,k)
910             
911! if (RIraf(ln,k) .eq. 1D-24 .or. RIraf(ln+1,k) .eq. 1D-24)
912!     &  rhock(ln,k)=0.
913! I limit the vertical velocity to 1km/s to avoid instabilities in velocity advection (TO BE IMPROVED)
914             
915! print*,'rhock',istep,ln,k,Zraf(ln),ln,rhock(ln,k),RIraf(ln,k),
916!     & Facesc
917!     &  RIraf(ln,k)*exp(-DZ*(alpha(ln)+beta(ln)+dzeta(ln))),
918!     &  (alpha(ln)+beta(ln)+dzeta(ln)),
919!     & (log(RIraf(ln+1,k))-log(RIraf(ln,k)))/Dz,Atri(ln+1),
920!     &  exp(-DZ*(alpha(ln)+beta(ln)+dzeta(ln)))
921             
922! if (abs(rhock(ln,k)) .gt. 1000. .or. isnan(rhock(ln,k))) then
923! print*,ig,istep
924! print*,'rhock',Zraf(ln),ln,rhock(ln,k),RIraf(ln+1,k)
925!     &  RIraf(ln,k)*exp(-DZ*(alpha(ln)+beta(ln)+dzeta(ln))),
926!     &  (alpha(ln)+beta(ln)+dzeta(ln)),
927!     & (log(RIraf(ln+1,k))-log(RIraf(ln,k)))/Dz,Atri(ln+1),
928!     &  exp(-DZ*(alpha(ln)+beta(ln)+dzeta(ln)))
929! if (abs(rhock(ln,k)) .gt. 10000.) stop
930! rhock(ln,k) =1000.*rhock(ln,k)/abs(rhock(ln,k))
931! endif
932             
933              enddo ! ln
934            rhock(naltvert,k)=ueff(ig,k)
935            ENDIF ! (iter .eq. niter)
936          ENDDO ! dynamic ions
937         
938! Upadte the electron density to assume neutral equilibrium
939         
940          IF (iter .eq. niter) THEN
941            RERaf=Rerafold
942            !ke=etrans(1)
943            DO l=1,naltvert
944              DO k=1,ndiffions
945                kn=itrans(k)
946                REraf(l)=REraf(l)+Charge_tr(kn)*mol_tr(ke)/mol_tr(kn)
947     &            * (RIraf(l,k)-Rirafold(l,k))
948              ENDDO
949            ! if (ig .eq. ij0) then
950            ! DO k=1,ndiffions
951            ! kn=itrans(k)
952            ! print*,'fin cycle',l,Rerafold(l)/Mtraceur(ke)/masseU/1d6,
953            !     &  Reraf(l)/Mtraceur(ke)/masseU/1d6,
954            !     &  Rirafold(l,k)/Mtraceur(kn)/masseU/1d6,
955            !     &  Riraf(l,k)/Mtraceur(kn)/masseU/1d6
956            ! enddo
957            ! endif
958            ENDDO
959          ENDIF ! (iter .eq. niter)
960         
961! Compute vertical component of electric field
962         
963!          DZ=Zraf(2)-Zraf(1)
964          Ez(1)=0D0
965          DO ln=2,naltvert-1
966            Ez(ln)=(log(Teraf(ln+1))-log(Teraf(ln-1)))+
967     &  (log(Reraf(ln+1))-log(Reraf(ln-1)))
968            Ez(ln)=-kbolt/qcharge*Teraf(ln)/2D0/Dz/H0*Ez(ln)
969          ENDDO
970       Ez(naltvert)=3D0*log(Reraf(naltvert))-4D0*log(Reraf(naltvert-1))
971     &  +log(Reraf(naltvert-2))+
972     &  3D0*log(Teraf(naltvert))-4D0*log(Teraf(naltvert-1))
973     &  +log(Teraf(naltvert-2))
974      Ez(naltvert)=-kbolt/qcharge*Teraf(naltvert)/2D0/Dz/H0*Ez(naltvert)
975!      Ez(naltvert)=0D0
976           
977!       if (ig .eq. ij0) then
978!         do ln=1,naltvert
979!          ! print*,'reraf write',ln,Reraf(ln)/Mtraceur(31)/masseU/1d6
980!           write(301,*) ig,istep,ln,Zraf(ln),RIraf(ln,:),rhock(ln,:),
981!     &        REraf(ln),TIraf(ln,:),Teraf(ln),Tnraf(ln)
982!         enddo
983!       endif
984         
985          ENDDO ! iteration
986         
987          tsim=tsim+tdiff
988! compute new time step
989          Wmax=maxval(abs(rhock))
990!          if (mu_local .ge. 0.) tdiff3=Dz/Wmax/300d0
991!          if (mu_local .lt. 0.) tdiff3=dz/Wmax/500d0
992          !if (mu_local .ge. 0.5) tdiff3=h0*Dz/Wmax/20d0!20d0 !20d0 !/100d0
993          !if (mu_local .le. -0.5) tdiff3=2d0*16d0*h0*dz/Wmax/20d0!20d0 !300d0
994          !if (abs(mu_local) .lt. 0.5) tdiff3=2d0*16d0*h0*dz/Wmax/20d0!20d0 !400d0 !/1000d0
995         
996          !! This is a weird factor to optimize time calculation
997          if (mu_local .ge. 0.3) tdiff3=2d0*16d0*h0*Dz/Wmax/20d0!20d0 !20d0 !/100d0
998          if (mu_local .le. -0.3) tdiff3=2d0*16d0*h0*dz/Wmax/20d0!20d0 !300d0
999          if (abs(mu_local) .lt. 0.3) tdiff3=2d0*16d0*h0*dz/Wmax/20d0!20d0 !400d0 !/1000d0
1000         
1001          tmin=max(tdiff1,tdiff3)
1002          tmin=min(tmin,tdiff2)
1003          tdiff=tmin
1004          !! This is a weird factor to optimize time calculation
1005          !!tmin = 2*0.45*(DZ*H0)**2./maxval(abs(DIraf(:,:)))
1006          tmin = 3.*8.*4.*200*2*0.45*(DZ*H0)**2.
1007          tmin = tmin/(maxval(abs(DIraf(:,:)))+(DZ*H0)*Wmax)
1008          if (tdiff .ge. tmin) then
1009            tdiff = tmin
1010          endif
1011          tdiff=max(tdiff1,tdiff)
1012          !if (szad_local .le. 20.) tdiff = tdiff/20.
1013          !!if (tsim .le. 8.75D-2) tdiff = 8.75D-5
1014          if (tsim+tdiff .gt. ttot) tdiff=ttot-tsim
1015         
1016!      if (ig .eq. ij0) then
1017!        print*,'dt',Wmax,tdiff1,tdiff2,tdiff3,tdiff,tsim,ttot,
1018!     &  alpha(naltvert-1),beta(naltvert-1),dzeta(naltvert-1),
1019!     &  Atri(naltvert),Xtri(naltvert),Xtri(naltvert-1)*Atri(naltvert),
1020!     &  h0,dz,h0*dz,mu_local
1021!      endif
1022         
1023        ENDDO ! time step ou while
1024        ! nsubreal=nsubvert
1025        ! end time step here
1026!        if (ig .eq. ij0) close(301)
1027       
1028        MtotZ2(:)=0D0
1029       
1030! Limites pour la vitesse verticale pour eviter problemes numeriques
1031       
1032        DO k=1,ndiffions
1033          DO ln=1,naltvert
1034            IF (abs(rhock(ln,k)).gt.Wlim .or. isnan(rhock(ln,k))) THEN
1035              rhock(ln,k) = Wlim*rhock(ln,k)/abs(rhock(ln,k))
1036            ENDIF
1037            !IF (k.eq.1) THEN
1038            !    write(*,*) real(qk1Dnew(ln,:)), qold(ig,ln,gcmind(ke))
1039            !ENDIF
1040          ENDDO ! ln
1041        ENDDO ! k
1042       
1043! print*,'rhock',l,k,Zraf(l),l,rhock(l,k),RIraf(l,k),Facesc
1044!        MtotZ2(k)=MtotZ2(k)+RIraf(l,k)*Dzraf
1045!        enddo
1046!        enddo
1047       
1048! We reinterpolate results on the Pressure grid and correct with FacQ
1049         
1050        CALL GCMGRID_P2K(Zraf,Praf,TNraf,Mraf,QIraf,RIraf,TIraf,TEraf,
1051     &  rhock,Ez,P1D,Pnlay1D,TempN1d,Mmean1d,qk1d,qk1dnew,
1052     &  TempI1D,TempI1Dnw,TempE1D,TempE1dnew,Wk1D,Wk1D2,Wk1d3,
1053     &  Ez1d,ZZ1d,naltvert,ndiffions,nlayer,
1054     &  iz_vertplasma,itrans,etrans,FacQ,FacTI,FacTe)
1055
1056        !iq = itrans(1)
1057        !write(*,*) tname(gcmind(iq)),ptimestep,qk1dnew(:,1)
1058
1059     
1060! do l=1,llm+1
1061! print*,'Ez,W',P1D(l),Ez1D(l),Wk1D(l,:)
1062! enddo
1063     
1064! Compute the mixing ratio and temperature & Update potential temperature
1065! Here I neglect advection of temperature (later)
1066     
1067!    do ln=1,llm
1068!      l=ln-iz_plasma
1069!      Alt(ig,ln)=ZZ1d(ln)
1070!      do k=1,ndiffions
1071!        kn=itrans(k)
1072!! if (ig .eq. ij0) print*,'qold',ln,q(ig,ln,kn),qk1d(ln,k)
1073!        q(ig,ln,kn)=real(qk1Dnew(ln,k))
1074!! if (ig .eq. ij0) print*,'qnew',q(ig,ln,kn),qk1dnew(ln,k)
1075!        if (l .gt. 0) then
1076!          Tetak(ig,l,k)=real(tempI1Dnw(ln,k))**(1D0-kappak(k))*
1077!     & (masseU*Mtraceur(kn)/kbolt*Preffplasma/RhoN(ig,ln)/q(ig,ln,kn))
1078!     & **(kappak(k))
1079!          wcontk(ig,l,k)=real(Wk1D(ln,k))
1080!          wcovk(ig,l,k)=real(Wk1D2(ln,k))
1081!          wphysk(ig,l,k)=real(Wk1D3(ln,k))
1082!          EfieldZ(ig,l)=real(Ez1D(ln))
1083!        endif
1084!     
1085!        if (q(ig,ln,kn) .lt. 0. .or. isnan(q(ig,ln,kn))) then
1086!          print*,'aie q < 0',q(ig,ln,kn),ln,l,ig,k,kn
1087!          print*,qk1dnew(:,k),qk1d(:,k)
1088!          stop
1089!        endif
1090!       
1091!        if (l .ge. 1) then
1092!          if (isnan(Wcontk(ig,l,k))) then
1093!            Wcontk(ig,l,k)=0.
1094!            Wcovk(ig,l,k)=0.
1095!        ! print*,'aie Wcontk',l,ln,ig,Wcontk(ig,l,k),Wk1D(ln,k)
1096!        ! do l=1,naltvert
1097!        ! print*,'rhock',l,rhock(l,:),RIraf(l,:),Dad(l),alpha(l),beta(l),
1098!        !     &  dzeta(l),Praf(l)
1099!        ! enddo
1100!          endif
1101!         if (isnan(Efieldz(ig,l))) then
1102!           Efieldz(ig,l)=0.
1103!         endif
1104!       endif
1105!     
1106!      enddo
1107!    enddo
1108!    do k=1,ndiffions
1109!      Wcontk(ig,llp+1,k)=real(Wk1D(llm+1,k))
1110!      Wcovk(ig,llp+1,k)=real(Wk1D2(llm+1,k))
1111!      wphysk(ig,llp+1,k)=real(Wk1D3(llm+1,k))
1112!    enddo
1113!   EfieldZ(ig,llp+1)=Ez1D(llm+1)
1114       
1115        DO ln=1,nlayer
1116          l=ln-iz_vertplasma
1117          DO k=1,ndiffions
1118            iq=gcmind(itrans(k))
1119            qnew(ig,ln,iq)=real(qk1Dnew(ln,k))
1120          ENDDO
1121        ENDDO
1122       
1123      ENDDO !  END ig -  Main Loop on horizontal grids - plan horizontal
1124     
1125! ==========================================================
1126! ===== Compute the diffusion trends du to diffusion =======
1127! ================== for the ions ==========================
1128! ==========================================================
1129     
1130      !ke = etrans(1)
1131      DO l=1,nlayer
1132        DO k=1,ndiffions
1133          iq = gcmind(itrans(k))
1134          pdqiondiff(:,l,iq)    = (qnew(:,l,iq)-qold(:,l,iq))/ptimestep
1135!         CE CALCUL EST FAUX, CE N'EST PAS UNE CONSERVATION MMR QU'IL FAUT MAIS UNE CONSERVATION VMR !!!!!!!!!!
1136!          pdqiondiff(:,l,iqelec)= pdqiondiff(:,l,iqelec)
1137!     &                                   + pdqiondiff(:,l,iq)
1138          !pdqiondiff(:,l,iq)    =(qnew(:,l,kn)-qold(:,l,kn))/ptimestep
1139          !pdqiondiff(:,l,iqelec)= pdqiondiff(:,l,ke)+pdqiondiff(:,l,iq)
1140        ENDDO
1141      ENDDO
1142      !iq = gcmind(itrans(1))
1143      !write(*,*) tname(iq),iq,ptimestep,pdqiondiff(1,:,iq),qold(1,:,iq)
1144      !write(*,*) maxval(pdqiondiff(:,:,iq))     
1145! ==========================================================
1146! ===== Compute the diffusion trends du to diffusion =======
1147! ================== for the electron ======================
1148! ==========================================================
1149     
1150      ! ------ update mmean ------ !
1151      mmean_new(:,:) = 0.
1152      do iq = 1,ntracers
1153         mmean_new(:,:) = mmean_new(:,:) +
1154     &             qnew(:,:,gcmind(iq))/M_tr(gcmind(iq))
1155      enddo
1156      mmean_new(:,:) = 1./mmean_new(:,:)
1157     
1158      ! ------ conversion mmr into vmr ------ !
1159      do iq = 1,ntracers
1160        vmr_new(:,:,gcmind(iq)) = qnew(:,:,gcmind(iq)) *
1161     &                            mmean_new(:,:)/M_tr(gcmind(iq))
1162      enddo
1163     
1164      iqelec = gcmind(ke)
1165      ! ------ vmr ion density ------ !
1166      vmr_ion(:,:) = 0.
1167      do iq = 1,ntracers
1168        if ((type_tr(gcmind(iq)) .eq. 2) .and.
1169     &      (  tname(gcmind(iq)) .ne. tname(iqelec))) then
1170          vmr_ion(:,:) = vmr_ion(:,:) + vmr_new(:,:,gcmind(iq))
1171        endif
1172      enddo
1173     
1174      ! ------ force charge neutrality ------ !
1175      ! ------ vmr(ion) = vmr(elec) ------ !
1176      DO ig=1,ngrid
1177        DO l=1,nlayer
1178          if (vmr_new(ig,l,iqelec) .ne. vmr_ion(ig,l)) then
1179            vmr_new(ig,l,iqelec) = vmr_ion(ig,l)
1180!            DO k=1,ndiffions
1181!              vmr_new(ig,l,iqelec) = vmr_new(ig,l,iqelec) +
1182!     &                               vmr_new(ig,l,gcmind(itrans(k)))
1183!            ENDDO
1184          endif
1185        ENDDO
1186      ENDDO
1187     
1188      ! ------ conversion vmr into mmr for electron ------ !
1189      qnew(:,:,iqelec) = vmr_new(:,:,iqelec)*m_tr(iqelec)/mmean_new(:,:)
1190     
1191      ! ------ Compute the diffusion trends du to diffusion ------ !*
1192      ! ------ for the electron ------ !
1193      DO l=1,nlayer
1194        pdqiondiff(:,l,iqelec) =
1195     &            (qnew(:,l,iqelec)-qold(:,l,iqelec))/ptimestep
1196      ENDDO
1197     
1198     
1199!      write(*,*) 'TATA EST PLUS LA'
1200!  do ig=1,ip1jm+1,iip1
1201!    do l=1,llp
1202!      ln=l+iz_plasma
1203!      do k=1,ndiffions
1204!        kn=itrans(k)
1205!        q(ig+iim,ln,kn)=q(ig,ln,kn)
1206!        tetak(ig+iim,l,k)=tetak(ig,l,k)
1207!        wcontk(ig+iim,l,k)=wcontk(ig,l,k)
1208!        if (ig .eq. ij0 .and. l .eq. 2 .and. k .eq.1)
1209!     &  print*,ij0+iim,wcontk(ij0,2,1),wcontk(ij0+iim,2,1)
1210!        wcovk(ig+iim,l,k)=wcovk(ig,l,k)
1211!        wphysk(ig+iim,l,k)=wphysk(ig,l,k)
1212!        Efieldz(ig+iim,l)=Efieldz(ig,l)
1213!        tetae(ig+iim,l)=tetae(ig,l)
1214!      enddo
1215!    enddo
1216!  enddo
1217! print*,'wphys3',ij0,ij0+iim,wcontk(ij0,2,1),wcontk(ij0+iim,2,1)
1218     
1219      RETURN
1220      END
1221     
1222!    ********************************************************************
1223!    ********************************************************************
1224!    ********************************************************************
1225     
1226      SUBROUTINE TridagDP(a,b,c,r,u,n)
1227!      USE mod_phys_lmdz_para, only: mpi_rank
1228!      parameter (nmax=5000)
1229!      dimension gam(nmax),a(n),b(n),c(n),r(n),u(n)
1230      real(kind=kind(1.D0)) gam(n),a(n),b(n),c(n),r(n),u(n)
1231      if(b(1).eq.0.)then
1232        stop 'tridag: error: b(1)=0 !!! '
1233      endif
1234      bet=b(1)
1235      u(1)=r(1)/bet
1236      do 11 j=2,n
1237        gam(j)=c(j-1)/bet
1238        bet=b(j)-a(j)*gam(j)
1239        if(bet.eq.0.) then
1240          stop 'tridag: error: bet=0 !!! '
1241        endif
1242        u(j)=(r(j)-a(j)*u(j-1))/bet       
124311    continue
1244      do 12 j=n-1,1,-1
1245        u(j)=u(j)-gam(j+1)*u(j+1)
124612    continue
1247      return
1248      END
1249     
1250     
1251! SUBROUTINE DYN2DIFFK(Pk,teta,P,mmean,rhoN,Preff,cpp,kappa,
1252!     &  q,Mtraceur,Tempk,tetak,TempE,tetaE,Preffplasma,itrans,etrans,
1253!     &  llm,nqtot,llp,ndiffions,ip1jmp1,iz_plasma,kappak,kappae,
1254!     &  Pnlay,TempN,Rhok,TempI,TetaI,rhotetaI,TempF,TetaF,rhotetaF)
1255!
1256!  USE infotrac, only: tname
1257! 
1258!  IMPLICIT NONE
1259! 
1260!  INTEGER :: llm,llp,nqtot,ndiffions,ip1jmp1,l,ig,k,ln,kn,iz_plasma
1261!  REAL :: Pk(ip1jmp1,llm),teta(ip1jmp1,llm)
1262!  REAL :: mmean(ip1jmp1,llm),rhoN(ip1jmp1,llm)
1263!  REAL :: P(ip1jmp1,llm+1),Pnlay(ip1jmp1,llm)
1264!  REAL :: q(ip1jmp1,llm,nqtot)
1265!  REAL :: Tempk(ip1jmp1,llp,ndiffions),tetak(ip1jmp1,llp,ndiffions)
1266!  REAL :: TempE(ip1jmp1,llp),TetaE(ip1jmp1,llp)
1267!  REAL :: Mtraceur(nqtot)
1268!  INTEGER :: itrans(ndiffions),etrans(1)
1269!  REAL :: Preff,Cpp,Preffplasma,kappa,unskappa
1270!  REAL :: kappak(ndiffions)
1271!  REAL :: kappae
1272! 
1273!  ! Output
1274!  REAL :: TempN(ip1jmp1,llm)
1275!  REAL :: Rhok(ip1jmp1,llm,nqtot)
1276!  REAL :: TempI(ip1jmp1,llm,ndiffions)
1277!  REAL :: TetaI(ip1jmp1,llm,ndiffions)
1278!  REAL :: rhotetaI(ip1jmp1,llm,ndiffions)
1279!  REAL :: TempF(ip1jmp1,llm)
1280!  REAL :: TetaF(ip1jmp1,llm)
1281!  REAL :: rhotetaF(ip1jmp1,llm)
1282!
1283! REAL :: masseU,kbolt
1284! masseU=1.660538782d-27
1285! kbolt=1.3806504d-23
1286!
1287! unskappa=1./kappa
1288!
1289!  Rhok(:,:,:)=0.
1290! DO ig=1,ip1jmp1
1291!   DO ln=1,llm
1292!     l=ln-iz_plasma
1293!     TempN(ig,ln)=Teta(ig,ln)*pk(ig,ln)/Cpp
1294!     Pnlay(ig,ln)=preff*(pk(ig,ln)/cpp)**unskappa
1295!     do k=1,nqtot
1296!        if (tname(k) .ne. "dust_number" .and. tname(k) .ne. "dust_mass"
1297!     &  .and. tname(k).ne."ccn_number".and.tname(k).ne."ccn_mass") then
1298!         if (q(ig,ln,k) .le. 0.) q(ig,ln,k)=1d-30
1299!         Rhok(ig,ln,k)=q(ig,ln,k)*RhoN(ig,ln)
1300!! if (q(ig,ln,k) .le. 0. .and. l .le. 0) Rhok(ig,ln,k)=1d-30*rhoN(ig,ln)
1301!       endif
1302!     enddo
1303!
1304!! Ions parameters
1305!     do k=1,ndiffions
1306!       kn=itrans(k)
1307!       if (ln .le. iz_plasma) TempI(ig,ln,k)=TempN(ig,ln)
1308!       if (ln .gt. iz_plasma) TempI(ig,ln,k)=Tempk(ig,l,k)
1309!     
1310!       if (ln .le. iz_plasma) TetaI(ig,ln,k)=teta(ig,ln)*
1311!     &  (Preffplasma**kappak(k))/(Preff**kappa)*
1312!!     &  (Mtraceur(kn)/mmean(ig,ln)/q(ig,ln,kn))**kappa
1313!     &  (rhoN(ig,ln)/mmean(ig,ln))**kappa*
1314!     & 1D0/(rhok(ig,ln,kn)/Mtraceur(kn))**kappak(k)*
1315!     & (kbolt*TempN(ig,ln)/masseu)**(kappa-kappak(k))
1316!        if (ln .gt. iz_plasma) TetaI(ig,ln,k)=Tetak(ig,l,k)
1317!       
1318!       rhotetaI(ig,ln,k)=rhok(ig,ln,kn)*TetaI(ig,ln,k)
1319!     enddo       
1320!
1321!! Electron parameters
1322!
1323!      kn=etrans(1)
1324!      if (ln .le. iz_plasma) TempF(ig,ln)=TempN(ig,ln)
1325!      if (ln .gt. iz_plasma) TempF(ig,ln)=TempE(ig,l)
1326!     
1327!      if (ln .le. iz_plasma) TetaF(ig,ln)=teta(ig,ln)*
1328!     &  (Preffplasma**kappae)/(Preff**kappa)*
1329!!     &  (Mtraceur(kn)/mmean(ig,l)/q(ig,ln,kn))**kappa
1330!     &  (rhoN(ig,ln)/mmean(ig,ln))**kappa*
1331!     &  1D0/(rhok(ig,ln,kn)/Mtraceur(kn))**kappae*
1332!     &  (kbolt*TempN(ig,ln)/masseU)**(kappa-kappae)
1333!      if (ln .gt. iz_plasma) TetaF(ig,ln)=TetaE(ig,l)
1334!     
1335!      rhotetaF(ig,ln)=rhok(ig,ln,kn)*TetaF(ig,ln)
1336!
1337!
1338!  ENDDO ! l or l=ln-iz_plasma
1339! ENDDO ! ig
1340! END
1341     
1342      SUBROUTINE ZVERTK(P,T,M,Z,nl,gsurf)
1343      IMPLICIT NONE
1344#include "YOMCST.h"
1345      integer :: nl,l
1346      REAL(kind=kind(1.D0)):: P(nl),T(nl),M(nl),Z(nl)
1347      REAL(kind=kind(1.D0)):: masseU,kbolt,Konst,g,z0,Hpm
1348      REAL(kind=kind(1.D0)):: Tmean,Mmean,gsurf
1349      masseU=1.e-3/RNAVO
1350      kbolt=RKBOL
1351      Konst=kbolt/masseU
1352      g=RG
1353     
1354      Z0=gsurf/g!0d0
1355      Z(1)=z0
1356      Hpm=Konst*T(1)/g/M(1)
1357     
1358      do l=2,nl
1359        Tmean=T(l-1)
1360        Mmean=M(l-1)
1361        Hpm=Konst*Tmean/g/Mmean
1362        Z(l)=Z(l-1)-Hpm*log(P(l)/P(l-1))
1363!        print*,'z',l,Z(l)
1364      enddo
1365     
1366      END
1367     
1368      SUBROUTINE UPPER_RESOLK(ZZ1D,P1D,TN1d,M1d,rhok1d,TI1d,Te1d,
1369     & mol_tr,nlayer,ndiffions,ntracers,nlraf,idiffZ,
1370     & itrans,etrans,Zraf,TNraf,Praf,Mraf,Qraf,rhokraf,
1371     & QIraf,TIraf,FIraf,QEraf,Teraf,IndZ,gcmind)
1372      USE infotrac_phy, only: tname
1373      IMPLICIT NONE
1374#include "YOMCST.h"
1375      INTEGER :: nlayer,ndiffions,ntracers,nlraf,idiffZ,iq,k,kn,ke,iz,l
1376      INTEGER :: indZ(1)
1377      INTEGER :: itrans(ndiffions),etrans(1),gcmind(ntracers)
1378      REAL :: mol_tr(ntracers)
1379      REAL(kind=kind(1.D0)),dimension(nlayer) :: P1D,TN1D,M1D,ZZ1D,Te1D
1380      REAL(kind=kind(1.D0)),dimension(nlayer,ntracers) :: rhok1d
1381      REAL(kind=kind(1.D0)),dimension(nlayer,ndiffions) :: TI1d
1382      REAL(kind=kind(1.D0)),dimension(nlraf) :: Zraf,Praf,TNraf
1383      REAL(kind=kind(1.D0)),dimension(nlraf) :: Mraf
1384      REAL(kind=kind(1.D0)),dimension(nlraf,ntracers) :: Qraf,Rhokraf
1385      REAL(kind=kind(1.D0)),dimension(nlraf,ndiffions) :: QIraf,FIraf
1386      REAL(kind=kind(1.D0)),dimension(nlraf,ndiffions) :: TIraf
1387      REAL(kind=kind(1.D0)),dimension(nlraf) :: Qeraf,TEraf
1388      REAL(kind=kind(1.D0)) :: kbolt,masseU,Konst,g,a0,apol
1389      REAL(kind=kind(1.D0)) :: mu,nlay,facZ,dZ
1390      masseU=1.e-3/RNAVO
1391      kbolt=RKBOL
1392      Konst=kbolt/masseU
1393      g=RG
1394       
1395      a0=2.9D-15   ! absolute constant (when densities expressed in /m3)
1396      !apol=2.911D0 ! in 10^-24 cm3 (Withers 2009 for CO2)
1397     
1398      Zraf(1)   = ZZ1D(idiffZ)
1399      Praf(1)   = P1D(idiffZ)
1400      TNraf(1)  = TN1D(idiffZ)
1401      Qraf(:,:) = 0d0
1402      Rhokraf(1:nlraf,1:ntracers)=0d0
1403      do iq=1,ntracers
1404      !        if (tname(iq) .ne. "dust_number" .and. tname(iq) .ne. "dust_mass"
1405      !     & .and. tname(iq).ne."ccn_number".and.tname(iq).ne."ccn_mass") then 
1406        Rhokraf(1,iq)=rhok1d(idiffZ,iq)
1407        Qraf(1,iq)=rhok1d(idiffZ,iq)/sum(rhok1d(idiffZ,1:ntracers))
1408      !        endif
1409      enddo
1410      Mraf(1)=1D0/sum(Qraf(1,1:ntracers)/mol_tr(1:ntracers))
1411      nlay=Praf(1)/kbolt/TNraf(1)
1412      do k=1,ndiffions
1413        kn=itrans(k)
1414        TIraf(1,k)=TI1D(idiffZ,k)
1415        QIraf(1,k)=Qraf(1,kn)
1416!        mu=mol_tr(kn)*Mraf(1)/(mol_tr(kn)+Mraf(1))
1417!        FIraf(1,k)=nlay*a0*Mraf(1)/(Mraf(1)+mol_tr(kn))*sqrt(apol/mu)
1418      enddo
1419      ke=etrans(1)
1420      Teraf(1)=TE1D(idiffZ)
1421      Qeraf(1)=Qraf(1,ke)
1422       
1423      Zraf(nlraf)=ZZ1d(nlayer)
1424      do l=2,nlraf-1
1425        Zraf(l)=Zraf(1)+(Zraf(nlraf)-Zraf(1))/DBLE(nlraf-1)*DBLE(l-1)
1426       
1427        iZ=1
1428        do while (ZZ1D(iz) .le. Zraf(l))
1429          iZ=iZ+1
1430        enddo
1431        Iz=iz-1
1432       
1433      ! indZ=maxloc(ZZ1D,MASK=ZZ1d < Zraf(l))
1434      ! print*,'indZ',indZ
1435      ! iZ=indZ(1)
1436        dZ=Zraf(l)-Zraf(l-1)
1437        facZ=(Zraf(l)-zz1d(iz))/(zz1d(iZ+1)-ZZ1d(iz))
1438        TNraf(l)=TN1D(iz)*(TN1D(iz+1)/TN1d(iz))**facZ
1439        do iq=1,ntracers
1440          if (Rhok1d(iz,iq) .gt. 0.) then
1441             Rhokraf(l,iq)=Rhok1d(iz,iq) *
1442     &                    (Rhok1d(iz+1,iq)/Rhok1d(iz,iq))**facZ
1443          endif
1444        enddo
1445        do iq=1,ntracers
1446          Qraf(l,iq)=Rhokraf(l,iq)/sum(Rhokraf(l,1:ntracers))
1447        enddo
1448        Mraf(l)=1D0/sum(Qraf(l,1:ntracers)/mol_tr(1:ntracers))
1449        nlay=sum(Rhokraf(l,1:ntracers))/Mraf(l)/masseU
1450        Praf(l)=nlay*kbolt*TNraf(l)
1451        do k=1,ndiffions
1452          kn=itrans(k)
1453          TIraf(l,k)=TI1D(iz,k)*(TI1D(iz+1,k)/TI1D(iz,k))**facZ
1454          Qiraf(l,k)=Qraf(l,kn)
1455!          mu=mol_tr(kn)*Mraf(l)/(mol_tr(kn)+Mraf(l))
1456!          FIraf(l,k)=nlay*a0*Mraf(l)/(Mraf(l)+mol_tr(kn))*sqrt(apol/mu)
1457        enddo
1458        Teraf(l)=TE1D(iz)*(TE1D(iz+1)/TE1D(iz))**facZ
1459        Qeraf(l)=Qraf(l,ke)
1460      enddo
1461     
1462      Zraf(nlraf)=ZZ1d(nlayer)
1463      TNraf(nlraf)=TN1D(nlayer)
1464!      Qraf(:,:)=0d0
1465      do iq=1,ntracers
1466        Rhokraf(nlraf,iq)=rhok1d(nlayer,iq)
1467        Qraf(nlraf,iq)=rhok1d(nlayer,iq)/sum(rhok1d(nlayer,1:ntracers))
1468      enddo
1469      Mraf(nlraf)=1D0/sum(Qraf(nlraf,1:ntracers)/mol_tr(1:ntracers))
1470      nlay=sum(rhokraf(nlraf,1:ntracers))/Mraf(nlraf)/masseU
1471      Praf(nlraf)=P1D(nlayer)
1472      do k=1,ndiffions
1473        kn=itrans(k)
1474        TIraf(nlraf,k)=TI1D(nlayer,k)
1475        QIraf(nlraf,k)=Qraf(nlraf,kn)
1476!        mu=mol_tr(kn)*Mraf(nlraf)/(mol_tr(kn)+Mraf(nlraf))
1477!        nlay=Praf(nlraf)/kbolt/TNraf(nlraf)
1478!        FIraf(nlraf,k)=nlay*a0*Mraf(nlraf)/(Mraf(nlraf)+mol_tr(kn))
1479!     &  *sqrt(apol/mu)
1480      enddo
1481      Teraf(nlraf)=TE1D(nlayer)
1482      QEraf(nlraf)=Qraf(nlraf,ke)
1483     
1484      apol=0D0
1485      do iq=1,ntracers
1486       
1487        ! REF: CRC Handbook CHEMISTRY and Physics - 95th EDITION 2014-2015
1488        if (tname(gcmind(iq)).eq.'co2') apol=2.911d0      ! in 10^-24 cm3 (Withers 2009 for CO2)
1489        if (tname(gcmind(iq)).eq.'co')  apol=1.95d0       ! in 10^-24 cm3
1490        if (tname(gcmind(iq)).eq.'o2')  apol=1.5689d0     ! in 10^-24 cm3
1491        if (tname(gcmind(iq)).eq.'n2')  apol=1.7403d0     ! in 10^-24 cm3
1492        if (tname(gcmind(iq)).eq.'h2')  apol=0.8042d0     ! in 10^-24 cm3
1493        if (tname(gcmind(iq)).eq.'he')  apol=0.20550522d0 ! in 10^-24 cm3
1494        if (tname(gcmind(iq)).eq.'h')   apol=0.666793d0   ! in 10^-24 cm3
1495        if (tname(gcmind(iq)).eq.'n')   apol=1.10d0       ! in 10^-24 cm3
1496        if (tname(gcmind(iq)).eq.'o')   apol=0.802d0      ! in 10^-24 cm3
1497       
1498        if(apol.ne.0d0) then
1499         
1500          do l=1,nlraf
1501            nlay=Rhokraf(l,iq)/Mraf(l)/masseU
1502            do k=1,ndiffions
1503              kn=itrans(k)
1504              mu=mol_tr(kn)*mol_tr(iq)/(mol_tr(kn)+mol_tr(iq))
1505              FIraf(l,k) = FIraf(l,k) + nlay * a0
1506     &    * mol_tr(iq)/(mol_tr(iq)+mol_tr(kn)) * sqrt(apol/mu)
1507            enddo
1508          enddo
1509         
1510        endif
1511        apol=0d0
1512       
1513      enddo
1514     
1515      END
1516     
1517      SUBROUTINE GCMGRID_PK(Zraf,Praf,Traf,Mraf,QIraf,RIraf,
1518     &  TIraf,TEraf,P1di,P1D,Tn1d,M1d,qk1d,qk1dnew,TI1D,TI1Dnew,
1519     &  TE1D,TE1Dnew,nlraf,ndiffions,nlayer,izv,itrans,etrans)
1520     
1521      IMPLICIT NONE
1522#include "YOMCST.h"   
1523      INTEGER :: nlraf,ndiffions,nlayer,l,iP,k,izv
1524      INTEGER,dimension(ndiffions) :: itrans
1525      INTEGER,dimension(1) :: etrans
1526      REAL(kind=kind(1.D0)),dimension(nlayer) :: P1d,Tn1d,M1d
1527      REAL(kind=kind(1.D0)),dimension(nlayer+1) :: P1di
1528      REAL(kind=kind(1.D0)),dimension(nlayer,ndiffions) :: Qk1d,Qk1dnew
1529      REAL(kind=kind(1.D0)),dimension(nlayer,ndiffions) :: TI1d,TI1dnew
1530      REAL(kind=kind(1.D0)),dimension(nlayer) :: TE1d,TE1Dnew
1531      REAL(kind=kind(1.D0)),dimension(nlayer,ndiffions) :: Rknew
1532      REAL(kind=kind(1.D0)),dimension(nlraf) :: Zraf,Praf,Traf,Mraf
1533      REAL(kind=kind(1.D0)),dimension(nlraf,ndiffions) :: QIraf,TIraf
1534      REAL(kind=kind(1.D0)),dimension(nlraf) :: TEraf
1535      REAL(kind=kind(1.D0)),dimension(nlraf,ndiffions) :: RIraf
1536      REAL(kind=kind(1.D0)) :: masseU,kbolt,Konst,g
1537      REAL(kind=kind(1.D0)) :: facP,rhoNloc
1538      masseU=1.e-3/RNAVO
1539      kbolt=RKBOL
1540      Konst=kbolt/masseU
1541      g=RG
1542     
1543     
1544      ! below layer iz_vertplasma no effect due to vertical diffusion
1545      do l=1,nlayer
1546        if (P1D(l) .ge. Praf(1)) then
1547          qk1dnew(l,:)=qk1d(l,:)
1548          Ti1dnew(l,:)=Ti1d(l,:)
1549          Te1dnew(l)=TE1d(l)
1550        endif
1551         
1552        if (P1D(l) .lt. Praf(1)) then
1553          iP=1
1554          do while(Praf(iP) .gt. P1D(l))
1555            iP=iP+1
1556          enddo
1557          IP=iP-1
1558          facP=(P1D(l)-Praf(iP))/(Praf(iP+1)-Praf(iP))
1559          rhoNloc=P1D(l)/kbolt/TN1D(l)*M1d(l)*masseU
1560         
1561          do k=1,ndiffions
1562            Rknew(l,k)=RIraf(iP,k)+(RIraf(iP+1,k)-RIraf(iP,k))*facP
1563            TI1dnew(l,k)=TIraf(IP,k)+(TIraf(iP+1,k)-TIraf(iP,k))*facP
1564          enddo
1565          TE1Dnew(l)=TEraf(iP)+(TEraf(iP+1)-TEraf(iP))*facP
1566         
1567          do k=1,ndiffions
1568            Qk1dnew(l,k)=Rknew(l,k)/rhoNloc
1569          enddo
1570         
1571        endif
1572      enddo
1573       
1574      END
1575     
1576      SUBROUTINE GCMGRID_P2K(Zraf,Praf,Traf,Mraf,QIraf,RIraf,
1577     &  Tiraf,Teraf,Wk,Ez,P1Di,P1D,Tn1d,M1d,qk1d,qk1dnew,
1578     &  TI1D,TI1Dnew,TE1D,TE1Dnew,Wk1D,Wk1D2,Wk1D3,Ez1d,ZZ1d,
1579     &  nlraf,ndiffions,nlayer,izv,itrans,etrans,FacQ,FacTI,FacTe)
1580     
1581      IMPLICIT NONE
1582#include "YOMCST.h"   
1583      INTEGER :: nlraf,ndiffions,nlayer,l,iP,k,izv
1584      INTEGER :: itrans(ndiffions)
1585      INTEGER :: etrans(1)
1586      REAL(kind=kind(1.D0)),dimension(nlayer) :: P1d,Tn1d,M1d,ZZ1d
1587      REAL(kind=kind(1.D0)),dimension(nlayer+1) :: P1di
1588      REAL(kind=kind(1.D0)),dimension(nlayer,ndiffions) :: Qk1d,Qk1dnew
1589      REAL(kind=kind(1.D0)),dimension(nlayer,ndiffions) :: TI1d,TI1dnew
1590      REAL(kind=kind(1.D0)),dimension(nlayer) :: TE1D,TE1Dnew
1591      REAL(kind=kind(1.D0)),dimension(nlayer+1,ndiffions) :: Wk1D,WK1d2
1592      REAL(kind=kind(1.D0)),dimension(nlayer+1,ndiffions) :: Wk1D3
1593      REAL(kind=kind(1.D0)),dimension(nlayer+1) :: Ez1D
1594!      REAL(kind=kind(1.D0)),dimension(nlayer,ndiffions) :: Rknew
1595      REAL(kind=kind(1.D0)),dimension(nlayer,ndiffions) :: FacQ,FacTI
1596      REAL(kind=kind(1.D0)),dimension(nlayer) :: FacTE
1597      REAL(kind=kind(1.D0)),dimension(nlraf) :: Zraf,Praf,Traf
1598      REAL(kind=kind(1.D0)),dimension(nlraf) :: Mraf
1599      REAL(kind=kind(1.D0)),dimension(nlraf,ndiffions) :: QIraf
1600      REAL(kind=kind(1.D0)),dimension(nlraf,ndiffions) :: TIraf
1601      REAL(kind=kind(1.D0)),dimension(nlraf) :: TEraf
1602      REAL(kind=kind(1.D0)),dimension(nlraf,ndiffions) :: RIraf
1603      REAL(kind=kind(1.D0)),dimension(nlraf,ndiffions) :: Wk
1604      REAL(kind=kind(1.D0)),dimension(nlraf) :: Ez
1605      REAL(kind=kind(1.D0)) masseU,kbolt,Konst,g
1606      REAL(kind=kind(1.D0)) facP,rhoNloc
1607      masseU=1.e-3/RNAVO
1608      kbolt=RKBOL
1609      Konst=kbolt/masseU
1610      g=RG
1611     
1612     
1613! below layer iz_vertplasma no effect due to vertical diffusion
1614      do l=1,nlayer
1615        if (P1D(l) .ge. Praf(1)) then
1616        qk1dnew(l,1:ndiffions)=qk1d(l,1:ndiffions)
1617        Ti1dnew(l,1:ndiffions)=Ti1d(l,1:ndiffions)
1618        endif
1619       
1620        if (P1di(l) .ge. Praf(1)) then
1621          Ez1D(l)=0D0
1622          Wk1D(l,1:ndiffions)=0D0
1623          Wk1D2(l,1:ndiffions)=0D0
1624          Wk1D3(l,1:ndiffions)=0D0
1625        endif
1626       
1627        if (P1D(l) .lt. Praf(1)) then
1628          iP=1
1629          do while(Praf(iP) .gt. P1D(l))
1630            iP=iP+1
1631          enddo
1632          iP=iP-1
1633         
1634!          indP=maxloc(Praf,mask=Praf < P1D(l))
1635!          IP=indP(1)-1
1636         
1637          facP=(P1D(l)-Praf(iP))/(Praf(iP+1)-Praf(iP))
1638          rhoNloc=P1D(l)/kbolt/TN1D(l)*M1d(l)*masseU
1639          do k=1,ndiffions
1640             Qk1dnew(l,k)=(RIraf(iP,k)+(RIraf(iP+1,k)-RIraf(iP,k))*facP)
1641     &  *facQ(l,k)/rhoNloc
1642             TI1dnew(l,k)=(TIraf(IP,k)+(TIraf(iP+1,k)-TIraf(iP,k))*facP)
1643     &  *facTI(l,k)
1644          enddo
1645          TE1Dnew(l)=(TEraf(iP)+(TEraf(iP+1)-TEraf(iP))*facP)*facTE(l)
1646         
1647!      if (l .eq. nlayer) then
1648!        print*,'iP',iP,Praf(iP),P1D(l),Praf(iP+1),nlraf,facP
1649!          print*,RIraf(iP,1),RIraf(iP+1,1),facQ(l,1),rhoNloc,facTI(l,1)
1650!          print*,Qk1Dnew(l,1),TI1Dnew(l,1)
1651!        endif
1652       
1653        endif
1654       
1655        if (P1Di(l) .lt. Praf(1)) then
1656          iP=1
1657          do while(Praf(iP) .gt. P1Di(l))
1658            iP=iP+1
1659          enddo
1660          iP=iP-1
1661         
1662          facP=(P1Di(l)-Praf(iP))/(Praf(iP+1)-Praf(iP))
1663          do k=1,ndiffions
1664            wk1d(l,k)=(Wk(iP,k)+(Wk(iP+1,k)-Wk(iP,k))*facP)
1665     &  /(ZZ1D(l)-ZZ1D(l-1))
1666            wk1d2(l,k)=(Wk(iP,k)+(Wk(iP+1,k)-Wk(iP,k))*facP)
1667     &  *(ZZ1D(l)-ZZ1D(l-1))
1668            wk1d3(l,k)=(Wk(iP,k)+(Wk(iP+1,k)-Wk(iP,k))*facP)
1669          enddo
1670          Ez1D(l)=(EZ(iP)+(Ez(iP+1)-Ez(iP))*facP)
1671        endif
1672     
1673      enddo ! l
1674      Wk1d(nlayer+1,1:ndiffions)=0D0
1675      Wk1d2(nlayer+1,1:ndiffions)=0D0
1676      Wk1d3(nlayer+1,1:ndiffions)=0D0
1677      Ez1D(nlayer+1)=0D0
1678     
1679     
1680      END
1681     
1682     
1683      SUBROUTINE CORRFACK(qbef,qaft,Tbef,Taft,Tebef,Teaf,
1684     &  Fq,Ft,Fte,nl,nq)
1685     
1686      IMPLICIT NONE
1687   
1688      INTEGER :: nl,nq,l,k
1689      REAL(kind=kind(1.D0)),dimension(nl,nq) :: qbef,qaft,Fq
1690      REAL(kind=kind(1.D0)),dimension(nl,nq) :: Tbef,Taft,Ft
1691      REAL(kind=kind(1.D0)),dimension(nl) :: Tebef,Teaf,Fte
1692     
1693      do l=1,nl
1694       do k=1,nq
1695          Fq(l,k)=qbef(l,k)/qaft(l,k)
1696          Ft(l,k)=Tbef(l,k)/Taft(l,k)
1697       enddo
1698       Fte(l)=Tebef(l)/Teaf(l)
1699      enddo
1700     
1701      END
1702     
1703      SUBROUTINE DCOEFFK(TI,NUI,M,DI,nl,nqdyn,nqtot,ind)
1704     
1705      IMPLICIT NONE
1706     
1707      INTEGER :: nl,l,nqdyn,nqtot,k,kn
1708      INTEGER,dimension(nqdyn) :: ind
1709      REAL(kind=kind(1.D0)),dimension(nl,nqdyn) :: NUI,TI,DI
1710      REAL,dimension(nqtot) :: M
1711      REAL(kind=kind(1.D0)) :: kbolt,masseU
1712      masseU=1.660538782d-27
1713      kbolt=1.3806504d-23
1714     
1715      do k=1,nqdyn
1716        kn=ind(k)
1717        do l=1,nl
1718          DI(l,k)=kbolt/masseU*TI(l,k)/NUI(l,k)/M(kn)
1719        enddo
1720      enddo
1721     
1722      END
1723     
1724      SUBROUTINE DIFFPARAMK(T,Te,RE,D,Dz,
1725     &  alpha,beta,gama,delta,eps,dzeta,eta,nl,Wmax)
1726     
1727      IMPLICIT NONE
1728     
1729      INTEGER :: nl,l
1730      REAL(kind=kind(1.D0)),DIMENSION(nl) :: T,D
1731      REAL(kind=kind(1.D0)),DIMENSION(nl) :: Te,Re
1732      REAL(kind=kind(1.D0)) :: DZ,Wmax
1733      REAL(kind=kind(1.D0)),DIMENSION(nl) :: alpha,beta,gama,delta,eps
1734      REAL(kind=kind(1.D0)),DIMENSION(nl) :: dzeta,eta
1735     
1736      alpha(1)=1D0/T(1)*(-3D0*T(1)+4D0*T(2)-T(3))/(2D0*DZ)+
1737     &  1D0/T(1)*(-3D0*Te(1)+4D0*Te(2)-Te(3))/(2D0*DZ)
1738      delta(1)=(-3D0*D(1)+4D0*D(2)-D(3))/(2D0*Dz)
1739      beta(1)=1D0/T(1)
1740      dzeta(1)=Te(1)/T(1)*
1741     & (-3D0*log(Re(1))+4D0*log(Re(2))-log(Re(3)))/(2D0*DZ)
1742     
1743! Add a limit to dzeta (test to avoid unstabilities)
1744     
1745!      if (abs(dzeta(1)) .ge. 1d0 .and. Wmax .ge. 500d0) then
1746!        dzeta(1)=1d0*dzeta(1)/abs(dzeta(1))
1747!      endif
1748     
1749!      dzeta(1)=Te(1)/T(1)/Re(1)*
1750!     & (-3D0*Re(1)+4D0*Re(2)-Re(3))/(2D0*Dz)
1751     
1752      do l=2,nl-1
1753        alpha(l)=1D0/T(l)*(T(l+1)-T(l-1))/(2D0*DZ)+
1754     &  1D0/T(l)*(Te(l+1)-Te(l-1))/(2D0*DZ)
1755      delta(l)=(D(l+1)-D(l-1))/(2D0*Dz)
1756      beta(l)=1D0/T(l)
1757      dzeta(l)=Te(l)/T(l)*(log(Re(l+1))-log(Re(l-1)))/(2D0*DZ)
1758!      dzeta(l)=Te(l)/T(l)/Re(l)*(Re(l+1)-Re(l-1))/(2D0*DZ)
1759     
1760!      if (abs(dzeta(l)) .ge. 1d0 .and. Wmax .ge. 500d0) then
1761!        dzeta(l)=1d0*dzeta(l)/abs(dzeta(l))
1762!      endif
1763     
1764      enddo
1765     
1766! Upper altitude values
1767     
1768      alpha(nl)=1D0/T(nl)*(3D0*T(nl)-4D0*T(nl-1)+T(nl-2))/(2D0*DZ)+
1769     &  1D0/T(nl)*(3D0*Te(nl)-4D0*Te(nl-1)+Te(nl-2))/(2D0*DZ)
1770        delta(nl)=(3D0*D(nl)-4D0*D(nl-1)+D(nl-2))/(2D0*DZ)
1771        beta(nl)=1D0/T(nl)
1772      dzeta(nl)=Te(nl)/T(nl)*
1773     & (3D0*log(Re(nl))-4D0*log(Re(nl-1))+log(Re(nl-2)))/(2D0*DZ)
1774     
1775!      if (abs(dzeta(nl)) .ge. 1d0 .and. Wmax .ge. 500d0) then
1776!        dzeta(nl)=1d0*dzeta(nl)/abs(dzeta(nl))
1777!      endif
1778     
1779!      dzeta(nl)=Te(nl)/T(nl)/Re(nl)*
1780!     & (3D0*Re(nl)-4D0*Re(nl-1)+Re(nl-2))/(2D0*Dz)
1781     
1782! Compute the gama and eps coefficients
1783! Lower altitude values
1784     
1785      gama(1)=(-3D0*beta(1) +4D0*beta(2) -beta(3) )/(2d0*dz)
1786      eps(1) =(-3D0*alpha(1)+4D0*alpha(2)-alpha(3))/(2d0*dz)
1787      eta(1) =(-3D0*dzeta(1)+4D0*dzeta(2)-dzeta(3))/(2D0*dz)
1788     
1789        do l=2,nl-1
1790          gama(l)=(beta(l+1) -beta(l-1) )/(2d0*Dz)
1791          eps(l) =(alpha(l+1)-alpha(l-1))/(2d0*Dz)
1792          eta(l) =(dzeta(l+1)-dzeta(l-1))/(2D0*Dz)
1793        end do
1794     
1795        gama(nl)=(3D0*beta(nl) -4D0*beta(nl-1) +beta(nl-2) )/(2D0*DZ)
1796        eps(nl) =(3D0*alpha(nl)-4D0*alpha(nl-1)+alpha(nl-2))/(2D0*DZ)
1797        eta(nl) =(3D0*dzeta(nl)-4D0*dzeta(nl-1)+dzeta(nl-2))/(2D0*dz)
1798!        print*,'test',alpha(nl-1),beta(nl-1),dzeta(nl-1)
1799     
1800     
1801        END
1802     
1803     
1804      SUBROUTINE MATCOEFFK(alpha,beta,gama,delta,eps,dzeta,eta,
1805     &  Dad,Rhoad,Facesc,dz,dt,A,B,C,D,Ti,Te,nl)
1806     
1807      IMPLICIT NONE
1808     
1809      INTEGER :: nl,l
1810      REAL*8,DIMENSION(nl) :: alpha,beta,gama,delta,eps,dzeta,eta
1811      REAL*8,DIMENSION(nl) :: Dad,RHoad,Ti,Te
1812      REAL*8 :: dz,dt,del1,del2,del3,FacEsc
1813      REAL*8,DIMENSION(nl) :: A,B,C,D
1814      del1=dt/2d0/dz
1815      del2=dt/dz/dz
1816      del3=dt
1817     
1818! lower boundary coefficients no change
1819      A(1)=0d0
1820      B(1)=1d0
1821      C(1)=0d0
1822      D(1)=rhoAd(1)
1823     
1824      do l=2,nl-1
1825        A(l)=(delta(l)+(alpha(l)+beta(l)+dzeta(l))*Dad(l))*del1
1826     &       -Dad(l)*del2
1827        B(l)=-(delta(l)*(alpha(l)+beta(l)+dzeta(l))
1828     &       +Dad(l)*(gama(l)+eps(l)+eta(l)))*del3
1829        B(l)=B(l)+1d0+2d0*Dad(l)*del2
1830        C(l)=-(delta(l)+(alpha(l)+beta(l)+dzeta(l))*Dad(l))*del1
1831     &       -Dad(l)*del2
1832        D(l)=rhoAD(l)
1833      enddo
1834     
1835! Upper boundary coefficients Diffusion profile
1836      C(nl)=0d0
1837      B(nl)=-1d0
1838      A(nl)=exp(-dZ*(alpha(nl-1)+beta(nl-1)+dzeta(nl-1)))*FacEsc
1839!      A(nl)=exp(-dZ*(alpha(nl)+beta(nl)+dzeta(nl)))*FacEsc
1840!      A(nl)=exp(-dZ*(alpha(nl-1)+beta(nl-1)))*FacEsc
1841!      A(nl)=exp(-dZ*(alpha(nl-1)+beta(nl-1))*
1842!     &  (Ti(nl-1)/(Te(nl-1)+Ti(nl-1))))*Facesc
1843      D(nl)=0D0
1844     
1845      END
1846     
1847      SUBROUTINE CheckmassK(X,Y,nl)
1848     
1849      IMPLICIT NONE
1850     
1851      INTEGER :: nl
1852      REAL*8,DIMENSION(nl) :: X,Y
1853      REAL*8 Xtot,Ytot
1854     
1855      Xtot=sum(X)
1856      Ytot=sum(Y)
1857     
1858      if (abs((Xtot-Ytot)/Xtot) .gt. 1d-4) then
1859        print*,'no conservation for mass',Xtot,Ytot
1860      endif
1861     
1862      END
1863     
1864     
Note: See TracBrowser for help on using the repository browser.