Index: /trunk/LMDZ.VENUS/deftank/physiq-96x96x90-chemistry-IONOSPHERE-IONDIFF.def
===================================================================
--- /trunk/LMDZ.VENUS/deftank/physiq-96x96x90-chemistry-IONOSPHERE-IONDIFF.def	(revision 3035)
+++ /trunk/LMDZ.VENUS/deftank/physiq-96x96x90-chemistry-IONOSPHERE-IONDIFF.def	(revision 3035)
@@ -0,0 +1,157 @@
+## $Header$
+#
+##  Include diurnal cycle or not
+cycle_diurne=y
+##  Include soil model or not
+soil_model=y
+## Use orodr or not for orography
+ok_orodr=n
+## Use orolf or not for orography
+ok_orolf=n
+## Use non-orographic Gravity Waves of not
+ok_gw_nonoro=y
+## Number of calls to the radiative transfer (per day)
+nbapp_rad=24000
+## Number of calls to the chemistry routines (per day)
+nbapp_chem=24000
+##  Flag for convection : 1 pour LMD, 2 pour Tiedtke, 3 KE(new version JYG), 30 KE(version IPCC AR4), 4 KE vect
+iflag_con=0
+#
+# orbital parameters
+##~~~~~~~~~~~~~~~~~~
+#
+# VENUS ##
+# R_ecc = 0.006787
+R_ecc = 0.
+R_peri = 0.
+# R_incl = 177.4
+R_incl = 0.
+# solar: effective flux, given at 1 UA
+solaire = 328.
+#
+# parameters for the tracers
+##~~~~~~~~~~~~~~~~~~~~~~~~~
+#
+# 0: nothing (passive tracers)
+# 1: pseudo-chemistry relaxation (phytrac_relax)
+# 2: surface emission (phytrac_emiss)
+# 3: full chemistry (phytrac_chem)
+tr_scheme = 3
+#
+# Reinitialization of tracer abundances 
+reinit_trac=n
+#
+# parameters for chemistry and microphysics
+##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+#
+# use chemistry ?
+ok_chem=y
+# use clouds ? (needed for chemistry)
+ok_cloud=y
+# if yes, which scheme:
+#  1 => simple scheme (Aurelien Stolzenbach)
+#  2 => full microphysical scheme (Sabrina Guilbon)
+cl_scheme=1
+# use sedimentation (goes with clouds)
+ok_sedim=y
+nb_mode=3
+#
+# parameters for the boundary layer
+##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+#
+iflag_pbl=8
+z0 = 0.01
+lmixmin = 35.
+ksta = 1.e-7
+ok_kzmin=n
+#
+# Surface thermal inertia
+inertie=2000
+#
+# dry convection parameters
+##~~~~~~~~~~~~~~~~~~~~~~~~~
+#
+iflag_ajs = 1
+
+# Solar radiation module
+##~~~~~~~~~~~~~~~~~~~~~~
+## 1 = Rainer Haus Tables
+## 2 = Generic solar module
+#
+solarchoice=1
+
+# Thermosphere and nlte parameters
+##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+#
+## Compute non-LTE or not
+callnlte=y
+## choice of species profile to use
+nltemodel=2
+## Compute CO2 IR absorption or not
+callnirco2=y
+##CO2 IR Absorption model
+nircorr=0
+##include thermosphere or not
+callthermos=y
+
+## Thermospheric options 
+##~~~~~~~~~~~~~~~~~~~~~~
+## Tuning of photochemistry for oxygen production ?
+tuneupperatm = n
+
+## include photolysis j online calculated or not
+ok_jonline = y
+
+## include ion chemistry or not
+ok_ionchem = y
+
+## include ion diffusion or not | WARNING: it is 2.5-3. time slower with
+#  BE CAREFULL IF (ok_iondiff = y)
+#  ==> nbapp_chem need to be equal to day_step in gcm.def
+ok_iondiff = n
+
+#Method to include solar variability?
+#0-> Old method   1-> Variability with E10.7 as observed
+solvarmod=0
+
+## (Solar min=70 ave=140 max=300)
+fixed_euv_value = 200.
+
+# value for the UV heating efficiency
+##(experimental values between 0.19 and 0.23, lower values may
+## be used to compensate for low 15 um cooling)
+euveff = 0.195       
+#
+#
+#
+#
+# Parameters for IOIPSL output files
+##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+## OLD. Now we use XIOS => see context_lmdz_physics.xml to taylor the output files
+#
+### OK_journe= y for daily output file histday.nc, =n no histday.nc output
+### Meaningless for Venus
+OK_journe=n
+### OK_mensuel= y for monthly output file histmth.nc, =n no histmth.nc
+### For Venus, only these averaged outputs
+OK_mensuel=n
+## rate (in days) at which the Venus histmth file is to be written               
+# sets the output rate in histmth and/or histins 
+ecritphy=0.1
+### OK_instan=y, make some "instantaneous" outputs (same rate as histmth)
+OK_instan=n
+# 
+# Output levels for the various output files
+#
+# output level for  "day" lev_histday
+# - lev_hist*=1 => baseline 2D fields
+# - lev_hist*=2 => baseline 3D fields (default)
+# - lev_hist*=3 => radiative transfert
+# - lev_hist*=4 => 3D tendencies
+# - lev_hist*=5 => tracers and others
+lev_histday=2
+#output level for "mth" lev_histmth 
+lev_histmth=2
+#output level for "ins" lev_histins 
+lev_histins=2
+
Index: /trunk/LMDZ.VENUS/libf/phyvenus/Iondiff_red.F
===================================================================
--- /trunk/LMDZ.VENUS/libf/phyvenus/Iondiff_red.F	(revision 3035)
+++ /trunk/LMDZ.VENUS/libf/phyvenus/Iondiff_red.F	(revision 3035)
@@ -0,0 +1,1864 @@
+!       SUBROUTINE plasmaphys(q,rhoN,P,pk,teta,mmean,mudyn,
+!    & tempk,tetak,tempe,tetae,wcontk,wcovk,Efieldz,wphysk,alt)
+      
+       SUBROUTINE iondiff_red(ngrid,nlayer,nqmx,pplay,pplev,
+!     &                        ptneutre,ptelect,ption,pq,
+     &                        ptneutre,pq,pphis,
+     &                        gmtime,lat,lon,
+     &                        ptimestep,pdteuv,pdtconduc,pdqdiff,
+     &                        pdqiondiff)
+      
+!****************************************************************
+!
+!     Ion ambipolar diffusion routine 
+!
+!     Authors: J-Y. Chaufray
+!     -------
+!
+!     J-Y. Chaufray 
+! 
+!     Vitesse verticale et ses effets sur nk, Tk, uk, vk, etc...
+!     Seule les effets du terme (wk - wn) sont pris en compte ici 
+!     Pour l instant seul les effets sur la densite nk sont pris en compte
+!
+!     Version: 14/11/2020
+!
+!     -------
+!
+!     2022/10/22: adapted and adding by Antoine Martinez
+!
+!*****************************************************************
+      
+      USE chemparam_mod
+      USE infotrac_phy, only: tname
+      !USE dimphy    
+      
+      use plasmaphys_venus_mod
+      use iono_h, only: temp_elect, temp_ion
+      
+      IMPLICIT NONE
+      
+! Solve the vertical dynamics equation
+! The equation is solved on a altitude refined grid after interpolation
+! final results are reinterpolated on the atmospheric pressure grids
+! Upper boundary
+! The velocity at upper altitude should be parametrized (escape)
+! Here I assume w2up (nlayer)=0
+! Lower boundaru 
+! W2(1) = O (The ion velocity is equal to neutral velocity)
+      
+!#include "dimensions.h"
+!#include "paramet.h"
+!#include "plasmadyn.h"
+!#include "comgeom.h"
+      
+! =========================
+! ==== INPUT / OUTPUT =====
+! =========================
+      
+! ====== INPUT ======
+      
+      INTEGER,intent(in) :: ngrid  ! number of atmospheric columns
+      INTEGER,intent(in) :: nlayer ! number of atmospheric layers
+      INTEGER,intent(in) :: nqmx   ! number of advected tracers
+      REAL,intent(in) :: ptimestep ! physical timestep (s)
+      REAL,intent(in) :: pphis(ngrid) ! geopotential of the surface (m2/s)
+      
+      REAL,intent(in) :: pplay(ngrid,nlayer)        ! mid-layer pressure (Pa)
+      REAL,intent(in) :: pplev(ngrid,nlayer+1)      ! inter-layer pressure (Pa)
+      REAL,intent(in) :: pq(ngrid,nlayer,nqmx)      ! mid-layer mass mixing ratio tracers (kg/kg)
+      
+      REAL,intent(in) :: ptneutre(ngrid,nlayer)     ! mid-layer NEUTRAL  temperature (K)
+!      REAL pdt(ngrid,nlayer)
+!      REAL,intent(in) :: ptelect(ngrid,nlayer)      ! mid-layer ELECTRON temperature (K)
+!      REAL,intent(in) :: ption(ngrid,nlayer)        ! mid-layer ION      temperature (K)
+      
+      REAL,intent(in) :: pdteuv(ngrid,nlayer)       ! Temperature EUV heating tendancy (K/s)
+      REAL,intent(in) :: pdtconduc(ngrid,nlayer)    ! Temperature conduction tendancy  (K/s)
+      REAL,intent(in) :: pdqdiff(ngrid,nlayer,nqmx) ! neutral mmr molecular diffusion tendancy (kg/kg/s)
+!      real pdq(ngrid,nlayer,nq)
+      
+      REAL,intent(in) :: gmtime     ! day fraction   ratio [from 0 to 1]
+      REAL,intent(in) :: lat(ngrid) ! grid latitude  [degree]
+      REAL,intent(in) :: lon(ngrid) ! grid longitude [degree]
+      
+! ====== OUTPUT ======
+      
+      REAL            :: pdqiondiff(ngrid,nlayer,nqmx) ! ion mmr tendancy (kg/kg/s)
+      
+! =========================
+! ======== LOCAL ==========
+! =========================
+      
+      integer, parameter :: t_elec_origin= 2
+      integer, parameter :: t_ion_origin = 1
+      !Electronic temperature. Index 1 -> Theis et al. 1980 - model data ; Index 2-> Theis et al. 1984 - model data 
+      !Ion        temperature. Index 1 -> from VIRA Model based on the model of Miller et al. 1980 & 1984. 
+      
+      REAL :: lon_sun, sza_local,szad_local,mu_local
+      REAL :: lon_local(ngrid), lat_local(ngrid) 
+      
+      REAL :: ptelect(ngrid,nlayer)      ! mid-layer ELECTRON temperature (K)
+      REAL :: ption(ngrid,nlayer)        ! mid-layer ION      temperature (K)
+      
+      REAL :: qnew(ngrid,nlayer,nqmx), qold(ngrid,nlayer,nqmx)
+      
+      REAL :: mmean_new(ngrid,nlayer), vmr_new(ngrid,nlayer,nqmx)
+      REAL :: vmr_ion(ngrid,nlayer)
+      
+      !REAL wcontk(ip1jmp1,llp+1,ndiffions),Efieldz(ip1jmp1,llp+1)      
+      !REAL wcovk(ip1jmp1,llp+1,ndiffions),wphysk(ip1jmp1,llp+1,ndiffions)
+      
+! 1D field I pass to double precision
+      REAL(kind=kind(1.D0)) tdiff1,tdiff2
+      REAL(kind=kind(1.D0)),dimension(nlayer+1) :: P1D,Ez1d
+      REAL(kind=kind(1.D0)),dimension(nlayer) :: Mmean1D,RhoN1D,qe1D
+      REAL(kind=kind(1.D0)),dimension(nlayer) :: Pnlay1D,ZZ1D,TempN1d
+      REAL(kind=kind(1.D0)),dimension(nlayer) :: TempE1D,TempE1dint
+      REAL(kind=kind(1.D0)),dimension(nlayer) :: FacTe,TempE1dnew
+      
+      REAL(kind=kind(1.D0)),dimension(naltvert) :: Zraf,Praf,Tnraf,Mraf
+      REAL(kind=kind(1.D0)),dimension(naltvert) :: Rraf,QEraf,REraf
+      REAL(kind=kind(1.D0)),dimension(naltvert) :: RErafold,TEraf,Ez
+      
+      !REAL(kind=kind(1.D0)) X(lld),Y(lld),a(lld),b(lld),c(lld)
+      !REAL(kind=kind(1.D0)) X2(lld),Y2(lld),a2(lld),b2(lld),c2(lld)
+      
+      INTEGER,dimension(1) :: indZ
+      
+      REAL(kind=kind(1.D0)) :: Rho0,Nu0,T0,H0,D0,dt,dtad,time0
+      
+      REAL(kind=kind(1.D0)) :: tdiff3,tmin,Wmax,ttot,tsim,tdiff
+      
+      REAL(kind=kind(1.D0)) :: pphis1D    ! geopotential of the surface (m2/s)
+      
+      REAL(kind=kind(1.D0)),save :: Wlim
+      
+      REAL(kind=kind(1.D0)),dimension(naltvert) :: RIad,Tnad,TIad,Tead
+      REAL(kind=kind(1.D0)),dimension(naltvert) :: RElad,Dad,Zad
+      REAL(kind=kind(1.D0)) :: Dzraf,Dz,FacEsc,Time
+      REAL(kind=kind(1.D0)),dimension(naltvert) :: alpha,beta,delta
+      REAL(kind=kind(1.D0)),dimension(naltvert) :: eps,gama,dzeta,eta
+      REAL(kind=kind(1.D0)),dimension(naltvert) :: Atri,Btri,Ctri
+      REAL(kind=kind(1.D0)),dimension(naltvert) :: Dtri,Xtri,Ytri
+      INTEGER :: ig,k,l,ln,kn,ij0,l0,ln0,iter,ld,iq,istep,ke,niter,nn
+      INTEGER :: iqelec
+      INTEGER :: nsubreal
+      LOGICAL,SAVE :: firstcall=.true.
+      
+      REAL(kind=kind(1.D0)) masse,invsgmu,Konst
+      REAL(kind=kind(1.D0)) masseU,kBolt,g,RgazP,Rvenus,Grav,Mvenus,PI
+      
+      integer,save :: ntracers
+      integer,dimension(:),allocatable,save :: gcmind ! array of GCM indexes
+      real,dimension(:),allocatable,save :: mol_tr    ! array of mol mass traceurs
+      real,dimension(:),allocatable,save :: Charge_tr ! array of electric charge traceurs
+      
+      integer,dimension(nqmx) :: indic_tracers,indic_iondiff
+      integer,parameter       :: nb_espIon_diff=24
+      character(len=20),dimension(nb_espIon_diff) :: ListeIonDiff
+      character(len=20)                           :: electrans(1)
+      
+      integer,dimension(:),allocatable,save :: itrans    ! array of sub index gcmind of fluid ions
+      integer,dimension(1),save             :: etrans    ! index of gcmind electron   
+      
+      INTEGER,save :: iz_vertplasma ! vertical index above which ion diffusion is computed
+      INTEGER,save :: ndiffions     ! nombre d'ions decrit
+      
+      ! ---------- ALLOCATABLE ARRAY ------------------------------
+      
+      ! (ngrid,nlayer,ntracers)
+      !REAL(kind=kind(1.D0)),dimension(:,:),allocatable :: qnew
+      !REAL(kind=kind(1.D0)),dimension(:,:),allocatable :: qold
+      
+      ! (nlayer,ntracers)
+      REAL(kind=kind(1.D0)),dimension(:,:),allocatable,save :: q1D
+      REAL(kind=kind(1.D0)),dimension(:,:),allocatable,save :: RhoK1D
+      
+      ! (nlayer,ndiffions)            
+      REAL(kind=kind(1.D0)),dimension(:,:),allocatable,save :: qk1d
+      REAL(kind=kind(1.D0)),dimension(:,:),allocatable,save :: qk1dnew
+      REAL(kind=kind(1.D0)),dimension(:,:),allocatable,save :: qk1dint
+      REAL(kind=kind(1.D0)),dimension(:,:),allocatable,save :: TempI1D
+      REAL(kind=kind(1.D0)),dimension(:,:),allocatable,save :: TempI1din
+      REAL(kind=kind(1.D0)),dimension(:,:),allocatable,save :: TempI1dnw
+      REAL(kind=kind(1.D0)),dimension(:,:),allocatable,save :: FacQ
+      REAL(kind=kind(1.D0)),dimension(:,:),allocatable,save :: FacTI
+      
+      ! (nlayer+1,ndiffions) 
+      REAL(kind=kind(1.D0)),dimension(:,:),allocatable,save :: Wk1d
+      REAL(kind=kind(1.D0)),dimension(:,:),allocatable,save :: Wk1d2
+      REAL(kind=kind(1.D0)),dimension(:,:),allocatable,save :: Wk1d3
+      
+      ! (ngrid,ndiffions)
+      REAL(kind=kind(1.D0)),dimension(:,:),allocatable,save :: ueff
+      
+      ! (ndiffions) 
+      REAL(kind=kind(1.D0)),dimension(:),allocatable,save :: Mtot
+      REAL(kind=kind(1.D0)),dimension(:),allocatable,save :: Mtot2
+      
+      ! ---------- ALLOCATABLE ARRAY ------------------------------
+      
+      ! (naltvert,ntracers)
+      REAL(kind=kind(1.D0)),dimension(:,:),allocatable,save :: Qraf
+      REAL(kind=kind(1.D0)),dimension(:,:),allocatable,save :: Rkraf
+      
+      ! (naltvert,ndiffions) 
+      REAL(kind=kind(1.D0)),dimension(:,:),allocatable,save :: QIraf
+      REAL(kind=kind(1.D0)),dimension(:,:),allocatable,save :: RIraf
+      REAL(kind=kind(1.D0)),dimension(:,:),allocatable,save :: RIrafold
+      REAL(kind=kind(1.D0)),dimension(:,:),allocatable,save :: TIraf
+      REAL(kind=kind(1.D0)),dimension(:,:),allocatable,save :: FIraf
+      REAL(kind=kind(1.D0)),dimension(:,:),allocatable,save :: DIraf
+      REAL(kind=kind(1.D0)),dimension(:,:),allocatable,save :: Rhock
+      
+      ! (ndiffions) 
+      REAL(kind=kind(1.D0)),dimension(:),allocatable,save :: MtotZ1
+      REAL(kind=kind(1.D0)),dimension(:),allocatable,save :: MtotZ2
+      
+! ==========================================================
+! ==========================================================
+! ========== FIRST CALL AND MAIN PARAMETERS ================
+! ==========================================================
+! ==========================================================
+      
+! ======= Initializations at first call ==========
+      
+      if (firstcall) then
+        
+      ! === On elimine la microphysique ==========
+        
+        ntracers=0
+        indic_tracers(1:nqmx)=0
+        
+        ! type_tr ==> 1: neutral gas, 2: ion gas, 3: liquid, 10: others
+        do iq=1,nqmx
+          if ((type_tr(iq) .le. 2.) .and. (type_tr(iq) .gt. 0.)) then
+            indic_tracers(iq)=1
+            ntracers=ntracers+1
+          endif
+        enddo
+        allocate(gcmind(ntracers))
+        allocate(mol_tr(ntracers))
+        allocate(Charge_tr(ntracers))
+        
+        ! Store gcm indexes in gcmind and molar masses in mol_tr
+        nn=0
+        do iq=1,nqmx 
+          if (indic_tracers(iq) .eq. 1) then
+            nn=nn+1
+            gcmind(nn) = iq
+            mol_tr(nn) = M_tr(iq)
+            Charge_tr(nn) = 0.D0
+          endif
+        enddo
+        
+      ! === On selectionne les especes ioniques ==========
+        
+        ndiffions =  0
+        indic_iondiff(1:nqmx) = 0
+        
+        ! define used ions
+        ListeIonDiff(1)="o2plus"
+        ListeIonDiff(2)="oplus"
+        ListeIonDiff(3)="cplus"
+        ListeIonDiff(4)="nplus"
+        ListeIonDiff(5)="co2plus"
+        ListeIonDiff(6)="noplus"
+        ListeIonDiff(7)="coplus"
+        ListeIonDiff(8)="hplus"
+        ListeIonDiff(9)="h2oplus"
+        ListeIonDiff(10)="h3oplus"
+        ListeIonDiff(11)="hcoplus"
+        ListeIonDiff(12)="n2plus"
+        ListeIonDiff(13)="ohplus"
+        
+        ! define electron   
+        electrans(1)="elec"     
+        
+        ! seek the index and number of ion/electron
+        do iq=1,ntracers
+          do nn=1,nb_espIon_diff
+            if (ListeIonDiff(nn) .eq. tname(gcmind(iq))) then
+              indic_iondiff(iq)=iq ! 1
+            endif
+            ! we don't diffuse electron but we readapt electron density at the end
+            if (electrans(1) .eq. tname(gcmind(iq))) then
+              etrans(1)=iq
+            endif
+          enddo
+          if (indic_iondiff(iq) .ge. 1) then
+            print*,'Ion specie ', tname(gcmind(iq)),
+     &               'diffused in Iondiff_red'
+            ndiffions=ndiffions+1
+          endif
+        enddo
+        print*,'number of diffused ion:',ndiffions
+        allocate(itrans(ndiffions))
+        
+        !   Store gcm ion indexes in itrans
+        nn=0
+        do iq=1,ntracers
+          if (indic_iondiff(iq) .ge. 1) then
+            nn=nn+1
+            itrans(nn)=indic_iondiff(iq) ! iq
+            Charge_tr(itrans(nn))=1.!qcharge
+          endif
+        enddo
+        Charge_tr(etrans(1))=-1.!qcharge*(-1.)
+        
+        ! find vertical index above which ion diffusion is computed
+        do l=1,nlayer
+          if (pplay(1,l) .gt. Pdiffion) then
+            iz_vertplasma = l
+          endif
+        enddo
+        iz_vertplasma = iz_vertplasma+1 ! seuil de la diffusion verticale
+        print*,'vertical index for ion diffusion',
+     &         iz_vertplasma,pplay(1,iz_vertplasma)
+        ! iz_plasma     = iz_vertplasma-3  ! seuil de la dynamique horizontale
+        
+      ! allocatation des tableaux dependants du nombre d especes diffusees
+        allocate(ueff(ngrid,ndiffions))
+        
+        ! ----------------------------------- ! 
+        
+        !allocate-(qnew(ngrid,nlayer,ntracers))
+        !allocate-(qold(ngrid,nlayer,ntracers))
+        
+        allocate(q1D(nlayer,ntracers))
+        allocate(RhoK1D(nlayer,ntracers))
+        
+        allocate(qk1d(nlayer,ndiffions))
+        allocate(qk1dnew(nlayer,ndiffions))
+        allocate(qk1dint(nlayer,ndiffions))
+        allocate(TempI1D(nlayer,ndiffions))
+        allocate(TempI1din(nlayer,ndiffions))
+        allocate(TempI1dnw(nlayer,ndiffions))
+        allocate(FacQ(nlayer,ndiffions))
+        allocate(FacTI(nlayer,ndiffions))
+        
+        allocate(Wk1d(nlayer+1,ndiffions))
+        allocate(Wk1d2(nlayer+1,ndiffions))
+        allocate(Wk1d3(nlayer+1,ndiffions))
+        
+        allocate(Mtot(ndiffions))
+        allocate(Mtot2(ndiffions))
+        
+        ! ----------------------------------- ! 
+        allocate(Qraf(naltvert,ntracers))
+        allocate(Rkraf(naltvert,ntracers))
+        
+        allocate(QIraf(naltvert,ndiffions))
+        allocate(RIraf(naltvert,ndiffions))
+        allocate(RIrafold(naltvert,ndiffions))
+        allocate(TIraf(naltvert,ndiffions))
+        allocate(FIraf(naltvert,ndiffions))
+        allocate(DIraf(naltvert,ndiffions))
+        allocate(Rhock(naltvert,ndiffions))
+        
+        allocate(MtotZ1(ndiffions))
+        allocate(MtotZ2(ndiffions))
+        
+        ! ----------------------------------- ! 
+        
+      ! Conditions aux limites
+        
+        Wlim=500. ! Maximal absolute value of vertical velocity (m/s)  (to avoid unstabilities in wdu/dz term)
+        ueff(1:ngrid,1:ndiffions)=0. ! Effusion velocity at top in m/s
+        !ueff(1:ngrid,1:ndiffions)=200.
+        
+        firstcall=.FALSE.
+      endif ! firstcall
+      
+      masseU = 1.660538782d-27
+      kBolt  = 1.3806504d-23
+      Konst  = kBolt/masseU
+      RgazP  = 8.314472 
+      PI     = 2.*ASIN(1.) ! 3.141592653
+      g      = 8.87d0
+      Rvenus = 6051800d0 ! Used to compute escape parameter no need to be more accurate
+      Grav   = 6.67d-11
+      Mvenus = 4.86d24
+      ij0    = 6000 ! For test
+      invsgmu= 1d0/g/masseU
+      
+      ke=etrans(1)
+      
+      alpha(1:naltvert)=0D0
+      beta(1:naltvert)=0D0
+      delta(1:naltvert)=0D0
+      eps(1:naltvert)=0D0
+      dzeta(1:naltvert)=0D0
+      eta(1:naltvert)=0D0
+!! Open Vertical file to check it
+! if (firstcall) then
+! open(unit=301,file='Plasma_vertical.dat',status='unknown')
+! firstcall=.FALSE.
+! endif
+!! First we derive only parameters useful for diffusion
+!! Extrapolation of paramaters to all altitudes (we solve the dynamics equation on larger altitude range)
+!!
+!!   ===> ON A PAS BESOIN DE CETTE PARTIE CAR ON EST DEJA DANS LA PHYSIQUE
+!!       CALL DYN2DIFFK(Pk,teta,P,mmean,rhoN,Preff,cpp,kappa,
+!!     &  q,Mtraceur,Tempk,tetak,TempE,TetaE,Preffplasma,itrans,etrans,
+!!     &  llm,nqtot,llp,ndiffions,ip1jmp1,iz_plasma,kappak,kappae,
+!!     &  PNlay,TempN,Rhok,TempI,TetaI,rhotetaI,TempF,TetaF,rhotetaF)
+!!
+      
+! ==========================================================
+! ==========================================================
+! ========== Main Loop on horizontal grids =================
+! ==========================================================
+! ==========================================================
+      
+! ========== Initialization =================
+      
+      dt=dble(ptimestep)      !dtplasma*iappvert
+      nsubreal=nsubvert
+      niter=2
+      qold(1:ngrid,1:nlayer,1:nqmx)=pq(1:ngrid,1:nlayer,1:nqmx)
+      qnew(1:ngrid,1:nlayer,1:nqmx)=pq(1:ngrid,1:nlayer,1:nqmx)
+!      tdiff=dt/nsubreal
+!      print*,'dt',dt
+      
+      lon_sun = (0.5 - gmtime)*2.*PI
+      lon_local(:) = lon(:)*PI/180.
+      lat_local(:) = lat(:)*PI/180.
+      
+! ========== Horizontal loop start ==========
+      
+      DO ig=1,ngrid
+        
+        sza_local  =  cos(lat_local(ig))*cos(lon_local(ig))
+     &               *cos(lon_sun) + cos(lat_local(ig))
+     &               *sin(lon_local(ig))*sin(lon_sun)
+        sza_local  = min(sza_local,1.)
+        sza_local  = max(-1.,sza_local)
+        sza_local  = acos(sza_local)                   ! RAD
+        szad_local = sza_local * 180. / PI             ! Degree
+        mu_local   = cos(sza_local) 
+        
+        ttot   = dt
+        tsim   = 0.
+        tdiff1 = dt/nsubvert ! 
+        tdiff2 = dt/nsubmin  !
+        !if (nsubreal .ne. nsubvert) print*,'ig',ig,tdiff
+        
+        ZZ1D(1:nlayer)     = 0D0
+        P1D(1:nlayer)      = 0D0
+        Pnlay1D(1:nlayer)  = 0D0
+        TempN1D(1:nlayer)  = 0D0
+        Mmean1D(1:nlayer)  = 0D0
+        RHON1D(1:nlayer)   = 0D0
+        RHOK1D(1:nlayer,1:ntracers)=0D0
+        Q1D(1:nlayer,1:ntracers)=0D0
+        TEMPI1D(1:nlayer,1:ndiffions)=0D0
+        !TETAI1D(1:nlayer,1:ndiffions)=0D0
+        TEMPE1D(1:nlayer)=0D0
+        !TETAE1D(1:nlayer)=0D0
+        Zraf(1:naltvert)=0D0
+        Praf(1:naltvert)=0D0
+        Tnraf(1:naltvert)=0D0
+        Mraf(1:naltvert)=0D0
+        Qraf(1:naltvert,1:ntracers)=0D0
+        Rkraf(1:naltvert,1:ntracers)=0D0
+        REraf(1:naltvert)=0D0
+        QIraf(1:naltvert,1:ndiffions)=0D0
+        TIraf(1:naltvert,1:ndiffions)=0D0
+        FIraf(1:naltvert,1:ndiffions)=0D0
+        QEraf(1:naltvert)=0D0
+        TEraf(1:naltvert)=0D0
+        rhock(1:naltvert,1:ndiffions)=0D0
+        
+! print*,'ig',ig,llm
+! First we compute 1d fields in double precision 
+! The mixing ratios slightly change here (< 0.1%)
+        
+        pphis1D       = pphis(ig)
+        P1D(nlayer+1) = pplev(ig,nlayer+1)
+        do l=1,nlayer 
+          P1D(l)      = pplev(ig,l)     
+          Pnlay1D(l)  = pplay(ig,l) 
+          TempN1D(l)  = ptneutre(ig,l) 
+          do iq=1,ntracers
+            q1d(l,iq)=qold(ig,l,gcmind(iq))
+          enddo
+        ! compute the mean molecular mass
+          Mmean1d(l)  = (1D0/sum(q1d(l,:)/dble(mol_tr(:))))
+          
+        ! Compute the total mass density (kg/m3)
+          rhok1d(l,:) = 0D0
+          RhoN1D(l)   = Pnlay1D(l)*Mmean1d(l)/TempN1D(l)/Konst
+          do iq=1,ntracers
+            rhok1d(l,iq)=q1d(l,iq)*RhoN1D(l)
+          enddo
+          
+          do k=1,ndiffions
+            kn=itrans(k)
+            qk1d(l,k)=q1D(l,kn)
+            !TempI1D(l,k)=ption(ig,l)
+            !TetaI1D(l,k)=TetaI(ig,l,k)
+            !rhotetaI1D(l,k)=rhotetaI(ig,l,k)
+          enddo
+          
+          !ke=etrans(1)
+          qe1D(l)=q1D(l,ke)
+          !TempE1D(l)=ptelect(ig,l)
+          !TetaE1D(l)=TetaF(ig,l)
+          !rhoTetaE1D(l)=rhoTetaf(ig,l)
+        enddo ! l (vertical levels)
+        
+! Compute colonne mass of each dynamic ion
+        
+        Mtot(1:ndiffions)=0D0
+        
+!        do k=1,ndiffions
+!        kn=itrans(k)
+!        do l=1,lld
+!        ln=l+iz_vertplasma
+!        print*,ln,llm
+!        Mtot(k)=Mtot(k)+qk1d(ln,k)*rhoN1d(ln)
+!        enddo
+!        enddo
+        
+!        print*,'Mtot before',Mtot
+        
+! Compute altitude of the scalar grid pressure level
+        
+        !if (isnan(TempN1d(1))) then
+        !  write (*,*) 'PROBLEME NAN: TITI'
+        !endif
+        
+        CALL ZVERTK(Pnlay1D,TempN1D,mmean1d,ZZ1d,nlayer,pphis1D)
+        
+! Compute lectron and Ion temperature of the scalar grid pressure level 
+        
+        do l=1,nlayer
+          ptelect(ig,l) = temp_elect(ZZ1d(l)/1000.,TempN1d(l),
+     &                               szad_local,t_elec_origin)  
+          ption(ig,l)   = temp_ion(ZZ1d(l)/1000.,TempN1d(l),
+     &                             szad_local,t_ion_origin)  
+          
+          TempE1D(l) = ptelect(ig,l)
+          !!! 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.
+          do k=1,ndiffions
+            !kn=itrans(k)
+            TempI1D(l,k)=max(0.5*ptelect(ig,l),TempN1d(l))
+          enddo  
+        enddo ! l (vertical levels)        
+        !write(*,*) szad_local,ZZ1d(:)/1000.,ptelect(ig,:), ption(ig,:)
+!        do l=1,llm
+!          print*,'Pnlay1D, Plev',Pnlay1D(l),P1D(l),rhok1D(l,1:ntracers)
+!        enddo
+! I use a better vertial resolution above iz_vertplasma and altitude grid interpolation
+        
+        indZ(1)=0
+        CALL UPPER_RESOLK(ZZ1D,Pnlay1D,TempN1d,mmean1d,rhok1d,TempI1d,
+     & TempE1D,mol_tr,nlayer,ndiffions,ntracers,naltvert,iz_vertplasma,
+     & itrans,etrans,Zraf,Tnraf,Praf,Mraf,Qraf,Rkraf,
+     & QIraf,TIraf,FIraf,Qeraf,Teraf,indZ,gcmind)
+        
+!       print*,'z last level, P last level',Zraf(:),Praf(:)
+        
+        DZraf=Zraf(2)-Zraf(1)
+        
+        do l=1,naltvert
+          do k=1,ndiffions
+            kn=itrans(k)
+            RIraf(l,k)=Rkraf(l,kn)
+          enddo
+          !ke=etrans(1)
+          REraf(l)=Rkraf(l,ke)
+        enddo
+        
+! Test
+!       do l=1,naltvert
+!        print*,'l',l,Zraf(l),Praf(l),Tnraf(l),Qiraf(l,1),Tiraf(l,1),
+!     &  Qeraf(l),Teraf(l),Riraf(l,1),Reraf(l)
+!        enddo
+        
+        
+! Before solving the diffusion equation, we reddo interpolation on GCM grid 
+! to derive mass correction factor due to interpolation (for mass conservation purpose)
+        
+        CALL GCMGRID_PK(Zraf,Praf,TNraf,Mraf,
+     &  QIraf,RIraf,TIraf,Teraf,
+     &  P1D,Pnlay1D,TempN1d,Mmean1d,
+     &  qk1d,qk1dint,TempI1D,TempI1Din,TempE1d,TempE1dint,
+     &  naltvert,ndiffions,nlayer,iz_vertplasma,itrans,etrans)
+
+        !DO ln=1,nlayer
+          !IF (k.eq.1) THEN
+        !    write(*,*) real(qk1Dint(ln,:)), qold(ig,ln,gcmind(ke))
+          !ENDIF
+        !ENDDO
+        
+        FacQ(1:nlayer,1:ndiffions)=1D0
+        FacTI(1:nlayer,1:ndiffions)=1D0
+        FacTE(1:nlayer)=1D0
+        
+        CALL CORRFACK(qk1d,qk1dint,TempI1d,TempI1din,TempE1d,
+     &  TempE1dint,FacQ,FacTI,FacTE,nlayer,ndiffions)
+        
+        CALL DCOEFFK(TIraf,FIraf,mol_tr,DIraf,
+     &  naltvert,ndiffions,ntracers,itrans)
+        
+        Rho0=1d-20
+        
+        T0=Tnraf(naltvert)!TIraf(naltvert,1)
+        
+        do ln=1,naltvert
+          Tnad(ln)=TNraf(ln)/T0
+        enddo
+        
+!        if(abs(lon(ig)-(48.75)).le.0.1) then
+!          if(abs(lat(ig)-(-80.625)).le.0.1) then
+!            write(*,*) 'ZZ1d= ', ZZ1d(iz_vertplasma:90)
+!            write(*,*) 'Zraf= ', Zraf(1:naltvert)
+!            write(*,*) 'Pnlay1D= ', Pnlay1D(iz_vertplasma:90)
+!            write(*,*) 'Praf= ', Praf(1:naltvert)
+!            DO k=1,ndiffions
+!            kn=itrans(k)
+!            if((tname(gcmind(kn)).eq.'hplus').or.
+!     &         (tname(gcmind(kn)).eq.'1oplus').or.
+!     &         (tname(gcmind(kn)).eq.'1noplus')) then
+!             write(*,*) tname(gcmind(kn))//' DEBUT'
+!             write(*,*) 'avant= ',
+!     &       qk1d(iz_vertplasma:90,k)
+!             write(*,*) 'apres= ',
+!     &       qk1dint(iz_vertplasma:90,k)
+!             write(*,*) 'interm= ', 
+!     &       Rkraf(1:naltvert,kn)
+!             write(*,*) 'corrmass= ',
+!     &       real(FacQ(iz_vertplasma:90,k))
+!            endif
+!            ENDDO
+!
+!          endif
+!        endif        
+
+! do l=1,llm
+!        print*,'effect interp',l,qk1d(l,2),qk1dint(l,2)*facQ(l,2),
+!     & TempI1d(l,2),TempI1din(l,2)*facTI(l,2),TempE1D(l),
+!     & TempE1Dint(l)*facTE(l),facQ(l,:),FacTI(l,:),FacTe(l)
+! enddo
+! stop
+        
+! Compute vertical velocity to derive subtime step
+        
+        Tdiff=1.E-9
+        DO k=1,ndiffions
+          kn=itrans(k)
+          H0=kbolt*T0/mol_tr(kn)/g/masseU
+          D0=DIraf(naltvert,k)
+          Time0=H0*H0/D0
+          Time=Tdiff/Time0
+          do ln=1,naltvert
+            RIad(ln)=RIraf(ln,k)/Rho0
+            RElad(ln)=Reraf(ln)/rho0
+            TIad(ln)=TIraf(ln,k)/T0
+            TEad(ln)=TEraf(ln)/T0
+            Zad(ln)=Zraf(ln)/H0
+            Dad(ln)=DIraf(ln,k)/D0
+          enddo
+          
+          DZ=Zad(2)-Zad(1)
+          
+          CALL DIFFPARAMK(Tiad,Tead,RElad,Dad,DZ,
+     &  alpha,beta,gama,delta,eps,dzeta,eta,naltvert,Wmax)
+          
+! if (ig .eq. ij0) then
+! print*,'zeta',dzeta(:)
+! print*,'eta',eta(:)
+! endif
+          
+          FacEsc=exp(-ueff(ig,k)*Time0/H0*DZ/Dad(naltvert-1))
+          
+          CALL MATCOEFFK(alpha,beta,gama,delta,eps,dzeta,eta,Dad,Riad,
+     &  FacEsc,dZ,time,Atri,Btri,Ctri,Dtri,Tiad,Tead,naltvert)
+          
+          rhock(1,k)=0D0
+          rhock(2,k)=0D0
+          do ln=3,naltvert-1
+      
+            if (ln .lt. naltvert-1) then
+              rhock(ln,k)=-Dad(ln)*(-RIraf(ln+2,k)+8.*RIraf(ln+1,k)
+     &  -8.*RIraf(ln-1,k)+RIraf(ln-2,k))/RIraf(ln,k)/12D0/DZ
+     &  -Dad(ln)*(alpha(ln)+beta(ln)+dzeta(ln))
+            endif
+            
+            if (ln .eq. naltvert-1) then
+              rhock(ln,k)=-Dad(ln)*
+     &  (log(RIraf(ln+1,k))-log(RIraf(ln,k)))/Dz
+     &  -Dad(ln)*(alpha(ln)+beta(ln)+dzeta(ln))
+            endif
+            
+            rhock(ln,k)=D0/H0*rhock(ln,k)
+          enddo
+          
+          rhock(naltvert,k)=ueff(ig,k)
+        enddo
+        
+        
+        Wmax=maxval(abs(rhock))
+        if (mu_local .ge. 0.3) tdiff3=h0*Dz/Wmax/300d0 !20d0 !/100d0
+        if (mu_local .le. -0.3) tdiff3=h0*dz/Wmax/300d0 !300d0
+        if (abs(mu_local) .lt. 0.3) tdiff3=h0*dz/Wmax/400d0 !/1000d0
+        
+        tmin=max(tdiff1,tdiff3)
+        tmin=min(tmin,tdiff2)
+         
+        nsubreal=anint(dt/tmin)
+                
+        tdiff=tmin
+
+        !! This is a weird factor to optimize time calculation
+        !! Here, it is to be sure that the first dtime is very low
+        tmin = 2*0.45*(DZ*H0)**2.
+        tmin = tmin/(maxval(abs(DIraf(:,:)))+(DZ*H0)*Wmax)
+        if (tdiff .ge. tmin) then
+          tdiff = tmin
+        endif
+!        if (ig .eq. ij0) then
+!          print*,'nreal',nsubreal,nsubvert,nsubmin,
+!       &  Wmax,tmin,tdiff,tdiff1,tdiff2,tdiff3,h0*dZ,tsim,ttot,h0,dZ
+!          write(301,*) nsubreal
+!        endif
+        
+        !tdiff=8.75D-5 !tmin
+        tdiff=tmin
+        
+! Fix a minimal value for tdiff
+        
+! ========== TEMPORAL Loop =================
+        
+! Loop over dynamic ions to solve the diffusion equation
+        
+! Compute total mass on the Z grid before diffusion
+        
+!        MtotZ1(:)=0D0
+!        do k=1,ndiffions
+!        do l=1,naltvert
+!        MtotZ1(k)=MtotZ1(k)+RIraf(l,k)*Dzraf
+!        enddo
+!        enddo
+!        print*,'MtotZ',MtotZ1
+        
+        !if (szad_local.ge.93.) tsim = ttot
+!    DO istep=1,nsubreal
+        DO WHILE (tsim .lt. ttot)
+! if (ig .eq. ij0) then 
+! print*,'tsim',tsim
+!        do l=1,naltvert
+! print*,'density1',l,REraf(l)/Mtraceur(31)/masseU/1d6,
+!     &                   Riraf(l,2)/Mtraceur(22)/masseU/1d6
+! enddo
+! endif
+        RErafold=REraf
+        RIrafold=RIraf
+        DO iter=1,niter ! needed to use approximatively electron density (t+dt)
+          DO ln=1,naltvert
+            RElad(ln)=REraf(ln)/Rho0
+            REraf(ln)=RErafold(ln)   
+          
+          ! if (ig .eq. ij0) then
+          ! print*,'density',RElad(ln)/Mtraceur(31)/masseU/1d6*rho0
+          !     &                  ,Riraf(ln,2)/Mtraceur(22)/masseU/1d6
+          ! endif
+          
+          ENDDO ! ln
+          
+          DO k=1,ndiffions
+            kn=itrans(k)
+            H0=kbolt*T0/mol_tr(kn)/g/masseU
+            D0=DIraf(naltvert,k)
+            Time0=H0*H0/D0
+            Time=Tdiff/Time0
+            
+! Compute the diffusion coefficient using the collion frequency
+            
+!            print*,'now adimension'
+            
+! Adimension all parameters by value at iz_vertplasma+1
+            
+            do ln=1,naltvert
+              RIad(ln)=RIraf(ln,k)/Rho0
+!              RElad(ln)=REraf(ln)/Rho0
+              TIad(ln)=TIraf(ln,k)/T0
+              TEad(ln)=TEraf(ln)/T0
+              Zad(ln)=Zraf(ln)/H0
+              Dad(ln)=DIraf(ln,k)/D0
+            enddo ! ln
+            
+            DZ=Zad(2)-Zad(1)
+            
+!            Tead(:)=0D0
+            
+            CALL DIFFPARAMK(Tiad,Tead,RElad,Dad,DZ,
+     &  alpha,beta,gama,delta,eps,dzeta,eta,naltvert,Wmax)
+            
+!            if (ig .eq. ij0) then
+!        print*,'zeta',dzeta(:)
+!        print*,'eta',eta(:)
+!        endif
+            
+            FacEsc=exp(-ueff(ig,k)*Time0/H0*DZ/Dad(naltvert-1))
+!     & *Tiad(ln-1)/(Tiad(ln-1)+Tead(ln-1)))
+            
+!            eta(:)=0D0
+!            dzeta(:)=0D0
+            
+!            do l=1,naltvert
+!            print*,l,alpha(l),dzeta(l),eta(l),RElad(l),Riad(l)
+!            enddo
+!            stop
+            
+            
+            CALL MATCOEFFK(alpha,beta,gama,delta,eps,dzeta,eta,Dad,
+     &  Riad,FacEsc,dZ,time,Atri,Btri,Ctri,Dtri,Tiad,Tead,naltvert)
+            
+            Xtri(:)=0D0
+            
+            Xtri=Riad
+            Ytri=Xtri
+             
+! Solve system
+            !write(*,*) tsim, iter
+            CALL TridagDP(Atri,Btri,Ctri,Dtri,Xtri,naltvert)
+
+            !STOP
+! Check mass conservation
+            
+! CALL CheckmassK(RIad,Xtri,naltvert)
+            
+! if (ig .eq. 1653) then
+! do l=1,naltvert
+! print*,'lphysnew',l,rho0*Xtri(l),k
+! enddo
+! endif
+            
+            do l=1,naltvert
+              if (iter .eq. niter) RIraf(l,k)=Rho0*Xtri(l)
+              if (RIraf(l,k) .lt. 0. .or. isnan(RIraf(l,k))) then 
+              ! print*,'phys q <0',ig,l,k,iter,istep,Xtri(l),Ytri(l),Riraf(l,k),
+              !     &  REraf(l),Tead(l)*T0,Tiad(l)*T0
+                RIraf(l,k)=1D-24
+              ! RIraf(l,k)=RIraf(l-1,k)*exp(-DZ*(alpha(l-1)+beta(l-1)+dzeta(l-1)))
+              ! nsubreal=nsubreal+50
+              ! if (nsubreal .gt. 2000) then
+              ! print*,'oula nsubreal',ig,nsubreal
+              ! endif
+              ! goto 1
+              endif
+            enddo ! l
+            
+! Update electron mass density for next iteration
+            
+            DO l=1,naltvert
+              REraf(l)=REraf(l)+0.5*Charge_tr(kn)*mol_tr(ke)/mol_tr(kn)
+     &                  * rho0*(Xtri(l)-Ytri(l))
+            ENDDO ! l
+            
+!      if (ig .eq. ij0 .and. k .eq. 2) then 
+!       do l=1,naltvert
+!         print*,'l1',l,Reraf(l)/mol_tr(ke)/masseU/1d6,
+!     &        Xtri(l)*rho0/mol_tr(kn)/masseU/1d6,
+!     &        Ytri(l)*rho0/mol_tr(kn)/masseU/1d6
+!       enddo ! l
+!      endif 
+            
+!       Compute vertical velocity in atmospheric frame
+            
+            IF (iter .eq. niter) THEN
+              rhock(1,k)=0D0
+              rhock(2,k)=0D0
+              do ln=3,naltvert-1
+!       rhock(ln,k)=-Dad(ln)*((RIraf(ln+1,k)-RIraf(ln-1,k))/2D0/Dz
+!     & +RIraf(ln,k)*(alpha(ln)+beta(ln)+dzeta(ln)))
+!       rhock(ln,k)=D0/H0*rhock(ln,k)/RIraf(ln,k)
+              
+              
+! rhock(ln,k)=-Dad(ln)*
+!     &  ((log(RIraf(ln+1,k))-log(RIraf(ln,k)))/Dz
+!     &   +(alpha(ln)+beta(ln))*(Tiad(ln)/(Tead(ln)+Tiad(ln))))
+!       rhock(ln,k)=D0/H0*rhock(ln,k)
+              
+              if (ln .lt. naltvert-1) then
+! rhock(ln,k)=-Dad(ln)*(-log(RIraf(ln+2,k))+8.*log(RIraf(ln+1,k))
+!     &  -8.*log(RIraf(ln-1,k))+log(RIraf(ln-2,k)))
+!     &  /12D0/Dz
+                rhock(ln,k)=-Dad(ln)*(-RIraf(ln+2,k)+8.*RIraf(ln+1,k)
+     &        -8.*RIraf(ln-1,k)+RIraf(ln-2,k))/RIraf(ln,k)/12D0/DZ
+     &        -Dad(ln)*(alpha(ln)+beta(ln)+dzeta(ln))
+              endif
+              
+              if (ln .eq. naltvert-1) then
+! rhock(ln,k)=-Dad(ln)*((Tiad(ln)+Tead(ln))/Tiad(ln)*
+!     &  (log(RIraf(ln+1,k))-log(RIraf(ln,k)))/Dz)
+!     &  -Dad(ln)*(alpha(ln)+beta(ln))
+              
+                rhock(ln,k)=-Dad(ln)*
+     &        (log(RIraf(ln+1,k))-log(RIraf(ln,k)))/Dz
+     &        -Dad(ln)*(alpha(ln)+beta(ln)+dzeta(ln))
+              
+              endif
+              
+              rhock(ln,k)=D0/H0*rhock(ln,k)
+              
+! if (RIraf(ln,k) .eq. 1D-24 .or. RIraf(ln+1,k) .eq. 1D-24)
+!     &  rhock(ln,k)=0.
+! I limit the vertical velocity to 1km/s to avoid instabilities in velocity advection (TO BE IMPROVED)
+              
+! print*,'rhock',istep,ln,k,Zraf(ln),ln,rhock(ln,k),RIraf(ln,k),
+!     & Facesc
+!     &  RIraf(ln,k)*exp(-DZ*(alpha(ln)+beta(ln)+dzeta(ln))),
+!     &  (alpha(ln)+beta(ln)+dzeta(ln)),
+!     & (log(RIraf(ln+1,k))-log(RIraf(ln,k)))/Dz,Atri(ln+1),
+!     &  exp(-DZ*(alpha(ln)+beta(ln)+dzeta(ln)))
+              
+! if (abs(rhock(ln,k)) .gt. 1000. .or. isnan(rhock(ln,k))) then
+! print*,ig,istep
+! print*,'rhock',Zraf(ln),ln,rhock(ln,k),RIraf(ln+1,k)
+!     &  RIraf(ln,k)*exp(-DZ*(alpha(ln)+beta(ln)+dzeta(ln))),
+!     &  (alpha(ln)+beta(ln)+dzeta(ln)),
+!     & (log(RIraf(ln+1,k))-log(RIraf(ln,k)))/Dz,Atri(ln+1),
+!     &  exp(-DZ*(alpha(ln)+beta(ln)+dzeta(ln)))
+! if (abs(rhock(ln,k)) .gt. 10000.) stop
+! rhock(ln,k) =1000.*rhock(ln,k)/abs(rhock(ln,k))
+! endif
+              
+              enddo ! ln
+            rhock(naltvert,k)=ueff(ig,k)
+            ENDIF ! (iter .eq. niter)
+          ENDDO ! dynamic ions
+          
+! Upadte the electron density to assume neutral equilibrium
+          
+          IF (iter .eq. niter) THEN
+            RERaf=Rerafold
+            !ke=etrans(1)
+            DO l=1,naltvert
+              DO k=1,ndiffions
+                kn=itrans(k)
+                REraf(l)=REraf(l)+Charge_tr(kn)*mol_tr(ke)/mol_tr(kn)
+     &            * (RIraf(l,k)-Rirafold(l,k))
+              ENDDO
+            ! if (ig .eq. ij0) then
+            ! DO k=1,ndiffions
+            ! kn=itrans(k)
+            ! print*,'fin cycle',l,Rerafold(l)/Mtraceur(ke)/masseU/1d6,
+            !     &  Reraf(l)/Mtraceur(ke)/masseU/1d6,
+            !     &  Rirafold(l,k)/Mtraceur(kn)/masseU/1d6,
+            !     &  Riraf(l,k)/Mtraceur(kn)/masseU/1d6
+            ! enddo
+            ! endif
+            ENDDO
+          ENDIF ! (iter .eq. niter)
+          
+! Compute vertical component of electric field
+          
+!          DZ=Zraf(2)-Zraf(1)
+          Ez(1)=0D0
+          DO ln=2,naltvert-1
+            Ez(ln)=(log(Teraf(ln+1))-log(Teraf(ln-1)))+
+     &  (log(Reraf(ln+1))-log(Reraf(ln-1)))
+            Ez(ln)=-kbolt/qcharge*Teraf(ln)/2D0/Dz/H0*Ez(ln)
+          ENDDO
+       Ez(naltvert)=3D0*log(Reraf(naltvert))-4D0*log(Reraf(naltvert-1))
+     &  +log(Reraf(naltvert-2))+
+     &  3D0*log(Teraf(naltvert))-4D0*log(Teraf(naltvert-1))
+     &  +log(Teraf(naltvert-2))
+      Ez(naltvert)=-kbolt/qcharge*Teraf(naltvert)/2D0/Dz/H0*Ez(naltvert)
+!      Ez(naltvert)=0D0
+           
+!       if (ig .eq. ij0) then
+!         do ln=1,naltvert
+!          ! print*,'reraf write',ln,Reraf(ln)/Mtraceur(31)/masseU/1d6
+!           write(301,*) ig,istep,ln,Zraf(ln),RIraf(ln,:),rhock(ln,:),
+!     &        REraf(ln),TIraf(ln,:),Teraf(ln),Tnraf(ln)
+!         enddo
+!       endif
+          
+          ENDDO ! iteration
+          
+          tsim=tsim+tdiff
+! compute new time step
+          Wmax=maxval(abs(rhock))
+!          if (mu_local .ge. 0.) tdiff3=Dz/Wmax/300d0
+!          if (mu_local .lt. 0.) tdiff3=dz/Wmax/500d0
+          !if (mu_local .ge. 0.5) tdiff3=h0*Dz/Wmax/20d0!20d0 !20d0 !/100d0
+          !if (mu_local .le. -0.5) tdiff3=2d0*16d0*h0*dz/Wmax/20d0!20d0 !300d0
+          !if (abs(mu_local) .lt. 0.5) tdiff3=2d0*16d0*h0*dz/Wmax/20d0!20d0 !400d0 !/1000d0
+          
+          !! This is a weird factor to optimize time calculation
+          if (mu_local .ge. 0.3) tdiff3=2d0*16d0*h0*Dz/Wmax/20d0!20d0 !20d0 !/100d0
+          if (mu_local .le. -0.3) tdiff3=2d0*16d0*h0*dz/Wmax/20d0!20d0 !300d0
+          if (abs(mu_local) .lt. 0.3) tdiff3=2d0*16d0*h0*dz/Wmax/20d0!20d0 !400d0 !/1000d0
+          
+          tmin=max(tdiff1,tdiff3)
+          tmin=min(tmin,tdiff2)
+          tdiff=tmin
+          !! This is a weird factor to optimize time calculation
+          !!tmin = 2*0.45*(DZ*H0)**2./maxval(abs(DIraf(:,:)))
+          tmin = 3.*8.*4.*200*2*0.45*(DZ*H0)**2.
+          tmin = tmin/(maxval(abs(DIraf(:,:)))+(DZ*H0)*Wmax)
+          if (tdiff .ge. tmin) then
+            tdiff = tmin 
+          endif
+          tdiff=max(tdiff1,tdiff)
+          !if (szad_local .le. 20.) tdiff = tdiff/20.
+          !!if (tsim .le. 8.75D-2) tdiff = 8.75D-5
+          if (tsim+tdiff .gt. ttot) tdiff=ttot-tsim
+          
+!      if (ig .eq. ij0) then
+!        print*,'dt',Wmax,tdiff1,tdiff2,tdiff3,tdiff,tsim,ttot,
+!     &  alpha(naltvert-1),beta(naltvert-1),dzeta(naltvert-1),
+!     &  Atri(naltvert),Xtri(naltvert),Xtri(naltvert-1)*Atri(naltvert),
+!     &  h0,dz,h0*dz,mu_local
+!      endif
+          
+        ENDDO ! time step ou while
+        ! nsubreal=nsubvert
+        ! end time step here
+!        if (ig .eq. ij0) close(301)
+        
+        MtotZ2(:)=0D0
+        
+! Limites pour la vitesse verticale pour eviter problemes numeriques
+        
+        DO k=1,ndiffions
+          DO ln=1,naltvert
+            IF (abs(rhock(ln,k)).gt.Wlim .or. isnan(rhock(ln,k))) THEN
+              rhock(ln,k) = Wlim*rhock(ln,k)/abs(rhock(ln,k))
+            ENDIF
+            !IF (k.eq.1) THEN
+            !    write(*,*) real(qk1Dnew(ln,:)), qold(ig,ln,gcmind(ke))
+            !ENDIF
+          ENDDO ! ln
+        ENDDO ! k
+        
+! print*,'rhock',l,k,Zraf(l),l,rhock(l,k),RIraf(l,k),Facesc
+!        MtotZ2(k)=MtotZ2(k)+RIraf(l,k)*Dzraf
+!        enddo
+!        enddo
+        
+! We reinterpolate results on the Pressure grid and correct with FacQ
+         
+        CALL GCMGRID_P2K(Zraf,Praf,TNraf,Mraf,QIraf,RIraf,TIraf,TEraf,
+     &  rhock,Ez,P1D,Pnlay1D,TempN1d,Mmean1d,qk1d,qk1dnew,
+     &  TempI1D,TempI1Dnw,TempE1D,TempE1dnew,Wk1D,Wk1D2,Wk1d3,
+     &  Ez1d,ZZ1d,naltvert,ndiffions,nlayer,
+     &  iz_vertplasma,itrans,etrans,FacQ,FacTI,FacTe)
+
+        !iq = itrans(1)
+        !write(*,*) tname(gcmind(iq)),ptimestep,qk1dnew(:,1)
+
+      
+! do l=1,llm+1
+! print*,'Ez,W',P1D(l),Ez1D(l),Wk1D(l,:)
+! enddo
+      
+! Compute the mixing ratio and temperature & Update potential temperature
+! Here I neglect advection of temperature (later)
+      
+!    do ln=1,llm
+!      l=ln-iz_plasma
+!      Alt(ig,ln)=ZZ1d(ln)
+!      do k=1,ndiffions
+!        kn=itrans(k)
+!! if (ig .eq. ij0) print*,'qold',ln,q(ig,ln,kn),qk1d(ln,k)
+!        q(ig,ln,kn)=real(qk1Dnew(ln,k))
+!! if (ig .eq. ij0) print*,'qnew',q(ig,ln,kn),qk1dnew(ln,k)
+!        if (l .gt. 0) then
+!          Tetak(ig,l,k)=real(tempI1Dnw(ln,k))**(1D0-kappak(k))*
+!     & (masseU*Mtraceur(kn)/kbolt*Preffplasma/RhoN(ig,ln)/q(ig,ln,kn))
+!     & **(kappak(k))
+!          wcontk(ig,l,k)=real(Wk1D(ln,k))
+!          wcovk(ig,l,k)=real(Wk1D2(ln,k))
+!          wphysk(ig,l,k)=real(Wk1D3(ln,k))
+!          EfieldZ(ig,l)=real(Ez1D(ln))
+!        endif
+!      
+!        if (q(ig,ln,kn) .lt. 0. .or. isnan(q(ig,ln,kn))) then
+!          print*,'aie q < 0',q(ig,ln,kn),ln,l,ig,k,kn
+!          print*,qk1dnew(:,k),qk1d(:,k)
+!          stop
+!        endif
+!        
+!        if (l .ge. 1) then
+!          if (isnan(Wcontk(ig,l,k))) then
+!            Wcontk(ig,l,k)=0.
+!            Wcovk(ig,l,k)=0.
+!        ! print*,'aie Wcontk',l,ln,ig,Wcontk(ig,l,k),Wk1D(ln,k)
+!        ! do l=1,naltvert
+!        ! print*,'rhock',l,rhock(l,:),RIraf(l,:),Dad(l),alpha(l),beta(l),
+!        !     &  dzeta(l),Praf(l)
+!        ! enddo
+!          endif
+!         if (isnan(Efieldz(ig,l))) then
+!           Efieldz(ig,l)=0.
+!         endif
+!       endif
+!      
+!      enddo
+!    enddo
+!    do k=1,ndiffions
+!      Wcontk(ig,llp+1,k)=real(Wk1D(llm+1,k))
+!      Wcovk(ig,llp+1,k)=real(Wk1D2(llm+1,k))
+!      wphysk(ig,llp+1,k)=real(Wk1D3(llm+1,k))
+!    enddo
+!   EfieldZ(ig,llp+1)=Ez1D(llm+1)
+        
+        DO ln=1,nlayer
+          l=ln-iz_vertplasma
+          DO k=1,ndiffions
+            iq=gcmind(itrans(k))
+            qnew(ig,ln,iq)=real(qk1Dnew(ln,k))
+          ENDDO
+        ENDDO
+        
+      ENDDO !  END ig -  Main Loop on horizontal grids - plan horizontal
+      
+! ==========================================================
+! ===== Compute the diffusion trends du to diffusion =======
+! ================== for the ions ==========================
+! ==========================================================
+      
+      !ke = etrans(1)
+      DO l=1,nlayer
+        DO k=1,ndiffions
+          iq = gcmind(itrans(k))
+          pdqiondiff(:,l,iq)    = (qnew(:,l,iq)-qold(:,l,iq))/ptimestep
+!         CE CALCUL EST FAUX, CE N'EST PAS UNE CONSERVATION MMR QU'IL FAUT MAIS UNE CONSERVATION VMR !!!!!!!!!!
+!          pdqiondiff(:,l,iqelec)= pdqiondiff(:,l,iqelec)
+!     &                                   + pdqiondiff(:,l,iq)
+          !pdqiondiff(:,l,iq)    =(qnew(:,l,kn)-qold(:,l,kn))/ptimestep
+          !pdqiondiff(:,l,iqelec)= pdqiondiff(:,l,ke)+pdqiondiff(:,l,iq)
+        ENDDO
+      ENDDO
+      !iq = gcmind(itrans(1))
+      !write(*,*) tname(iq),iq,ptimestep,pdqiondiff(1,:,iq),qold(1,:,iq)
+      !write(*,*) maxval(pdqiondiff(:,:,iq))      
+! ==========================================================
+! ===== Compute the diffusion trends du to diffusion =======
+! ================== for the electron ======================
+! ==========================================================
+      
+      ! ------ update mmean ------ !
+      mmean_new(:,:) = 0.
+      do iq = 1,ntracers
+         mmean_new(:,:) = mmean_new(:,:) + 
+     &             qnew(:,:,gcmind(iq))/M_tr(gcmind(iq))
+      enddo
+      mmean_new(:,:) = 1./mmean_new(:,:)
+      
+      ! ------ conversion mmr into vmr ------ !
+      do iq = 1,ntracers
+        vmr_new(:,:,gcmind(iq)) = qnew(:,:,gcmind(iq)) * 
+     &                            mmean_new(:,:)/M_tr(gcmind(iq))
+      enddo
+      
+      iqelec = gcmind(ke)
+      ! ------ vmr ion density ------ !
+      vmr_ion(:,:) = 0.
+      do iq = 1,ntracers
+        if ((type_tr(gcmind(iq)) .eq. 2) .and. 
+     &      (  tname(gcmind(iq)) .ne. tname(iqelec))) then
+          vmr_ion(:,:) = vmr_ion(:,:) + vmr_new(:,:,gcmind(iq))
+        endif
+      enddo
+      
+      ! ------ force charge neutrality ------ !
+      ! ------ vmr(ion) = vmr(elec) ------ !
+      DO ig=1,ngrid
+        DO l=1,nlayer
+          if (vmr_new(ig,l,iqelec) .ne. vmr_ion(ig,l)) then
+            vmr_new(ig,l,iqelec) = vmr_ion(ig,l)
+!            DO k=1,ndiffions
+!              vmr_new(ig,l,iqelec) = vmr_new(ig,l,iqelec) + 
+!     &                               vmr_new(ig,l,gcmind(itrans(k)))
+!            ENDDO
+          endif
+        ENDDO
+      ENDDO
+      
+      ! ------ conversion vmr into mmr for electron ------ !
+      qnew(:,:,iqelec) = vmr_new(:,:,iqelec)*m_tr(iqelec)/mmean_new(:,:)
+      
+      ! ------ Compute the diffusion trends du to diffusion ------ !*
+      ! ------ for the electron ------ !
+      DO l=1,nlayer
+        pdqiondiff(:,l,iqelec) = 
+     &            (qnew(:,l,iqelec)-qold(:,l,iqelec))/ptimestep
+      ENDDO
+      
+      
+!      write(*,*) 'TATA EST PLUS LA'
+!  do ig=1,ip1jm+1,iip1
+!    do l=1,llp
+!      ln=l+iz_plasma
+!      do k=1,ndiffions
+!        kn=itrans(k)
+!        q(ig+iim,ln,kn)=q(ig,ln,kn)
+!        tetak(ig+iim,l,k)=tetak(ig,l,k)
+!        wcontk(ig+iim,l,k)=wcontk(ig,l,k)
+!        if (ig .eq. ij0 .and. l .eq. 2 .and. k .eq.1) 
+!     &  print*,ij0+iim,wcontk(ij0,2,1),wcontk(ij0+iim,2,1)
+!        wcovk(ig+iim,l,k)=wcovk(ig,l,k)
+!        wphysk(ig+iim,l,k)=wphysk(ig,l,k)
+!        Efieldz(ig+iim,l)=Efieldz(ig,l)
+!        tetae(ig+iim,l)=tetae(ig,l)
+!      enddo
+!    enddo
+!  enddo
+! print*,'wphys3',ij0,ij0+iim,wcontk(ij0,2,1),wcontk(ij0+iim,2,1)
+      
+      RETURN
+      END
+      
+!    ********************************************************************
+!    ********************************************************************
+!    ********************************************************************
+      
+      SUBROUTINE TridagDP(a,b,c,r,u,n)
+!      USE mod_phys_lmdz_para, only: mpi_rank
+!      parameter (nmax=5000)
+!      dimension gam(nmax),a(n),b(n),c(n),r(n),u(n)
+      real(kind=kind(1.D0)) gam(n),a(n),b(n),c(n),r(n),u(n)
+      if(b(1).eq.0.)then
+        stop 'tridag: error: b(1)=0 !!! '
+      endif
+      bet=b(1)
+      u(1)=r(1)/bet
+      do 11 j=2,n
+        gam(j)=c(j-1)/bet
+        bet=b(j)-a(j)*gam(j)
+        if(bet.eq.0.) then
+          stop 'tridag: error: bet=0 !!! '
+        endif
+        u(j)=(r(j)-a(j)*u(j-1))/bet        
+11    continue
+      do 12 j=n-1,1,-1
+        u(j)=u(j)-gam(j+1)*u(j+1)
+12    continue
+      return
+      END
+      
+      
+! SUBROUTINE DYN2DIFFK(Pk,teta,P,mmean,rhoN,Preff,cpp,kappa,
+!     &  q,Mtraceur,Tempk,tetak,TempE,tetaE,Preffplasma,itrans,etrans,
+!     &  llm,nqtot,llp,ndiffions,ip1jmp1,iz_plasma,kappak,kappae,
+!     &  Pnlay,TempN,Rhok,TempI,TetaI,rhotetaI,TempF,TetaF,rhotetaF)
+!
+!  USE infotrac, only: tname
+!  
+!  IMPLICIT NONE
+!  
+!  INTEGER :: llm,llp,nqtot,ndiffions,ip1jmp1,l,ig,k,ln,kn,iz_plasma
+!  REAL :: Pk(ip1jmp1,llm),teta(ip1jmp1,llm)
+!  REAL :: mmean(ip1jmp1,llm),rhoN(ip1jmp1,llm)
+!  REAL :: P(ip1jmp1,llm+1),Pnlay(ip1jmp1,llm)
+!  REAL :: q(ip1jmp1,llm,nqtot)
+!  REAL :: Tempk(ip1jmp1,llp,ndiffions),tetak(ip1jmp1,llp,ndiffions)
+!  REAL :: TempE(ip1jmp1,llp),TetaE(ip1jmp1,llp)
+!  REAL :: Mtraceur(nqtot)
+!  INTEGER :: itrans(ndiffions),etrans(1)
+!  REAL :: Preff,Cpp,Preffplasma,kappa,unskappa
+!  REAL :: kappak(ndiffions)
+!  REAL :: kappae
+!  
+!  ! Output
+!  REAL :: TempN(ip1jmp1,llm)
+!  REAL :: Rhok(ip1jmp1,llm,nqtot)
+!  REAL :: TempI(ip1jmp1,llm,ndiffions)
+!  REAL :: TetaI(ip1jmp1,llm,ndiffions)
+!  REAL :: rhotetaI(ip1jmp1,llm,ndiffions)
+!  REAL :: TempF(ip1jmp1,llm)
+!  REAL :: TetaF(ip1jmp1,llm)
+!  REAL :: rhotetaF(ip1jmp1,llm)
+!
+! REAL :: masseU,kbolt
+! masseU=1.660538782d-27
+! kbolt=1.3806504d-23
+!
+! unskappa=1./kappa
+!
+!  Rhok(:,:,:)=0.
+! DO ig=1,ip1jmp1
+!   DO ln=1,llm
+!     l=ln-iz_plasma
+!     TempN(ig,ln)=Teta(ig,ln)*pk(ig,ln)/Cpp
+!     Pnlay(ig,ln)=preff*(pk(ig,ln)/cpp)**unskappa
+!     do k=1,nqtot
+!        if (tname(k) .ne. "dust_number" .and. tname(k) .ne. "dust_mass"
+!     &  .and. tname(k).ne."ccn_number".and.tname(k).ne."ccn_mass") then
+!         if (q(ig,ln,k) .le. 0.) q(ig,ln,k)=1d-30
+!         Rhok(ig,ln,k)=q(ig,ln,k)*RhoN(ig,ln)
+!! if (q(ig,ln,k) .le. 0. .and. l .le. 0) Rhok(ig,ln,k)=1d-30*rhoN(ig,ln)
+!       endif
+!     enddo
+!
+!! Ions parameters 
+!     do k=1,ndiffions
+!       kn=itrans(k)
+!       if (ln .le. iz_plasma) TempI(ig,ln,k)=TempN(ig,ln)
+!       if (ln .gt. iz_plasma) TempI(ig,ln,k)=Tempk(ig,l,k)
+!      
+!       if (ln .le. iz_plasma) TetaI(ig,ln,k)=teta(ig,ln)*
+!     &  (Preffplasma**kappak(k))/(Preff**kappa)*
+!!     &  (Mtraceur(kn)/mmean(ig,ln)/q(ig,ln,kn))**kappa
+!     &  (rhoN(ig,ln)/mmean(ig,ln))**kappa*
+!     & 1D0/(rhok(ig,ln,kn)/Mtraceur(kn))**kappak(k)*
+!     & (kbolt*TempN(ig,ln)/masseu)**(kappa-kappak(k))
+!        if (ln .gt. iz_plasma) TetaI(ig,ln,k)=Tetak(ig,l,k)
+!        
+!       rhotetaI(ig,ln,k)=rhok(ig,ln,kn)*TetaI(ig,ln,k)
+!     enddo        
+!
+!! Electron parameters
+!
+!      kn=etrans(1)
+!      if (ln .le. iz_plasma) TempF(ig,ln)=TempN(ig,ln)
+!      if (ln .gt. iz_plasma) TempF(ig,ln)=TempE(ig,l)
+!      
+!      if (ln .le. iz_plasma) TetaF(ig,ln)=teta(ig,ln)*
+!     &  (Preffplasma**kappae)/(Preff**kappa)*
+!!     &  (Mtraceur(kn)/mmean(ig,l)/q(ig,ln,kn))**kappa
+!     &  (rhoN(ig,ln)/mmean(ig,ln))**kappa*
+!     &  1D0/(rhok(ig,ln,kn)/Mtraceur(kn))**kappae*
+!     &  (kbolt*TempN(ig,ln)/masseU)**(kappa-kappae)
+!      if (ln .gt. iz_plasma) TetaF(ig,ln)=TetaE(ig,l)
+!      
+!      rhotetaF(ig,ln)=rhok(ig,ln,kn)*TetaF(ig,ln)
+!
+!
+!  ENDDO ! l or l=ln-iz_plasma
+! ENDDO ! ig
+! END
+      
+      SUBROUTINE ZVERTK(P,T,M,Z,nl,gsurf)
+      IMPLICIT NONE
+#include "YOMCST.h"
+      integer :: nl,l
+      REAL(kind=kind(1.D0)):: P(nl),T(nl),M(nl),Z(nl)
+      REAL(kind=kind(1.D0)):: masseU,kbolt,Konst,g,z0,Hpm
+      REAL(kind=kind(1.D0)):: Tmean,Mmean,gsurf
+      masseU=1.e-3/RNAVO
+      kbolt=RKBOL
+      Konst=kbolt/masseU
+      g=RG
+      
+      Z0=gsurf/g!0d0
+      Z(1)=z0
+      Hpm=Konst*T(1)/g/M(1)
+      
+      do l=2,nl
+        Tmean=T(l-1)
+        Mmean=M(l-1)
+        Hpm=Konst*Tmean/g/Mmean
+        Z(l)=Z(l-1)-Hpm*log(P(l)/P(l-1))
+!        print*,'z',l,Z(l)
+      enddo
+      
+      END
+      
+      SUBROUTINE UPPER_RESOLK(ZZ1D,P1D,TN1d,M1d,rhok1d,TI1d,Te1d,
+     & mol_tr,nlayer,ndiffions,ntracers,nlraf,idiffZ,
+     & itrans,etrans,Zraf,TNraf,Praf,Mraf,Qraf,rhokraf,
+     & QIraf,TIraf,FIraf,QEraf,Teraf,IndZ,gcmind)
+      USE infotrac_phy, only: tname
+      IMPLICIT NONE
+#include "YOMCST.h"
+      INTEGER :: nlayer,ndiffions,ntracers,nlraf,idiffZ,iq,k,kn,ke,iz,l
+      INTEGER :: indZ(1)
+      INTEGER :: itrans(ndiffions),etrans(1),gcmind(ntracers)
+      REAL :: mol_tr(ntracers)
+      REAL(kind=kind(1.D0)),dimension(nlayer) :: P1D,TN1D,M1D,ZZ1D,Te1D
+      REAL(kind=kind(1.D0)),dimension(nlayer,ntracers) :: rhok1d
+      REAL(kind=kind(1.D0)),dimension(nlayer,ndiffions) :: TI1d
+      REAL(kind=kind(1.D0)),dimension(nlraf) :: Zraf,Praf,TNraf
+      REAL(kind=kind(1.D0)),dimension(nlraf) :: Mraf
+      REAL(kind=kind(1.D0)),dimension(nlraf,ntracers) :: Qraf,Rhokraf
+      REAL(kind=kind(1.D0)),dimension(nlraf,ndiffions) :: QIraf,FIraf
+      REAL(kind=kind(1.D0)),dimension(nlraf,ndiffions) :: TIraf
+      REAL(kind=kind(1.D0)),dimension(nlraf) :: Qeraf,TEraf
+      REAL(kind=kind(1.D0)) :: kbolt,masseU,Konst,g,a0,apol
+      REAL(kind=kind(1.D0)) :: mu,nlay,facZ,dZ
+      masseU=1.e-3/RNAVO
+      kbolt=RKBOL
+      Konst=kbolt/masseU
+      g=RG
+       
+      a0=2.9D-15   ! absolute constant (when densities expressed in /m3)
+      !apol=2.911D0 ! in 10^-24 cm3 (Withers 2009 for CO2)
+      
+      Zraf(1)   = ZZ1D(idiffZ)
+      Praf(1)   = P1D(idiffZ)
+      TNraf(1)  = TN1D(idiffZ)
+      Qraf(:,:) = 0d0
+      Rhokraf(1:nlraf,1:ntracers)=0d0
+      do iq=1,ntracers
+      !        if (tname(iq) .ne. "dust_number" .and. tname(iq) .ne. "dust_mass"
+      !     & .and. tname(iq).ne."ccn_number".and.tname(iq).ne."ccn_mass") then  
+        Rhokraf(1,iq)=rhok1d(idiffZ,iq)
+        Qraf(1,iq)=rhok1d(idiffZ,iq)/sum(rhok1d(idiffZ,1:ntracers))
+      !        endif
+      enddo
+      Mraf(1)=1D0/sum(Qraf(1,1:ntracers)/mol_tr(1:ntracers))
+      nlay=Praf(1)/kbolt/TNraf(1)
+      do k=1,ndiffions
+        kn=itrans(k)
+        TIraf(1,k)=TI1D(idiffZ,k)
+        QIraf(1,k)=Qraf(1,kn)
+!        mu=mol_tr(kn)*Mraf(1)/(mol_tr(kn)+Mraf(1))
+!        FIraf(1,k)=nlay*a0*Mraf(1)/(Mraf(1)+mol_tr(kn))*sqrt(apol/mu)
+      enddo
+      ke=etrans(1)
+      Teraf(1)=TE1D(idiffZ) 
+      Qeraf(1)=Qraf(1,ke)
+       
+      Zraf(nlraf)=ZZ1d(nlayer)
+      do l=2,nlraf-1
+        Zraf(l)=Zraf(1)+(Zraf(nlraf)-Zraf(1))/DBLE(nlraf-1)*DBLE(l-1)
+       
+        iZ=1
+        do while (ZZ1D(iz) .le. Zraf(l))
+          iZ=iZ+1
+        enddo
+        Iz=iz-1 
+       
+      ! indZ=maxloc(ZZ1D,MASK=ZZ1d < Zraf(l))
+      ! print*,'indZ',indZ
+      ! iZ=indZ(1)
+        dZ=Zraf(l)-Zraf(l-1)
+        facZ=(Zraf(l)-zz1d(iz))/(zz1d(iZ+1)-ZZ1d(iz))
+        TNraf(l)=TN1D(iz)*(TN1D(iz+1)/TN1d(iz))**facZ
+        do iq=1,ntracers
+          if (Rhok1d(iz,iq) .gt. 0.) then
+             Rhokraf(l,iq)=Rhok1d(iz,iq) * 
+     &                    (Rhok1d(iz+1,iq)/Rhok1d(iz,iq))**facZ
+          endif
+        enddo
+        do iq=1,ntracers
+          Qraf(l,iq)=Rhokraf(l,iq)/sum(Rhokraf(l,1:ntracers))
+        enddo
+        Mraf(l)=1D0/sum(Qraf(l,1:ntracers)/mol_tr(1:ntracers))
+        nlay=sum(Rhokraf(l,1:ntracers))/Mraf(l)/masseU
+        Praf(l)=nlay*kbolt*TNraf(l)
+        do k=1,ndiffions
+          kn=itrans(k)
+          TIraf(l,k)=TI1D(iz,k)*(TI1D(iz+1,k)/TI1D(iz,k))**facZ
+          Qiraf(l,k)=Qraf(l,kn)
+!          mu=mol_tr(kn)*Mraf(l)/(mol_tr(kn)+Mraf(l))
+!          FIraf(l,k)=nlay*a0*Mraf(l)/(Mraf(l)+mol_tr(kn))*sqrt(apol/mu)
+        enddo
+        Teraf(l)=TE1D(iz)*(TE1D(iz+1)/TE1D(iz))**facZ
+        Qeraf(l)=Qraf(l,ke)
+      enddo
+      
+      Zraf(nlraf)=ZZ1d(nlayer)
+      TNraf(nlraf)=TN1D(nlayer)
+!      Qraf(:,:)=0d0
+      do iq=1,ntracers
+        Rhokraf(nlraf,iq)=rhok1d(nlayer,iq)
+        Qraf(nlraf,iq)=rhok1d(nlayer,iq)/sum(rhok1d(nlayer,1:ntracers))
+      enddo
+      Mraf(nlraf)=1D0/sum(Qraf(nlraf,1:ntracers)/mol_tr(1:ntracers))
+      nlay=sum(rhokraf(nlraf,1:ntracers))/Mraf(nlraf)/masseU
+      Praf(nlraf)=P1D(nlayer)
+      do k=1,ndiffions
+        kn=itrans(k)
+        TIraf(nlraf,k)=TI1D(nlayer,k)
+        QIraf(nlraf,k)=Qraf(nlraf,kn)
+!        mu=mol_tr(kn)*Mraf(nlraf)/(mol_tr(kn)+Mraf(nlraf))
+!        nlay=Praf(nlraf)/kbolt/TNraf(nlraf)
+!        FIraf(nlraf,k)=nlay*a0*Mraf(nlraf)/(Mraf(nlraf)+mol_tr(kn))
+!     &  *sqrt(apol/mu)
+      enddo
+      Teraf(nlraf)=TE1D(nlayer)
+      QEraf(nlraf)=Qraf(nlraf,ke)
+      
+      apol=0D0
+      do iq=1,ntracers
+        
+        ! REF: CRC Handbook CHEMISTRY and Physics - 95th EDITION 2014-2015
+        if (tname(gcmind(iq)).eq.'co2') apol=2.911d0      ! in 10^-24 cm3 (Withers 2009 for CO2)
+        if (tname(gcmind(iq)).eq.'co')  apol=1.95d0       ! in 10^-24 cm3 
+        if (tname(gcmind(iq)).eq.'o2')  apol=1.5689d0     ! in 10^-24 cm3 
+        if (tname(gcmind(iq)).eq.'n2')  apol=1.7403d0     ! in 10^-24 cm3 
+        if (tname(gcmind(iq)).eq.'h2')  apol=0.8042d0     ! in 10^-24 cm3 
+        if (tname(gcmind(iq)).eq.'he')  apol=0.20550522d0 ! in 10^-24 cm3 
+        if (tname(gcmind(iq)).eq.'h')   apol=0.666793d0   ! in 10^-24 cm3
+        if (tname(gcmind(iq)).eq.'n')   apol=1.10d0       ! in 10^-24 cm3 
+        if (tname(gcmind(iq)).eq.'o')   apol=0.802d0      ! in 10^-24 cm3 
+        
+        if(apol.ne.0d0) then
+          
+          do l=1,nlraf
+            nlay=Rhokraf(l,iq)/Mraf(l)/masseU
+            do k=1,ndiffions
+              kn=itrans(k)
+              mu=mol_tr(kn)*mol_tr(iq)/(mol_tr(kn)+mol_tr(iq))
+              FIraf(l,k) = FIraf(l,k) + nlay * a0 
+     &    * mol_tr(iq)/(mol_tr(iq)+mol_tr(kn)) * sqrt(apol/mu) 
+            enddo
+          enddo
+          
+        endif
+        apol=0d0
+        
+      enddo
+      
+      END
+      
+      SUBROUTINE GCMGRID_PK(Zraf,Praf,Traf,Mraf,QIraf,RIraf,
+     &  TIraf,TEraf,P1di,P1D,Tn1d,M1d,qk1d,qk1dnew,TI1D,TI1Dnew,
+     &  TE1D,TE1Dnew,nlraf,ndiffions,nlayer,izv,itrans,etrans)
+      
+      IMPLICIT NONE
+#include "YOMCST.h"    
+      INTEGER :: nlraf,ndiffions,nlayer,l,iP,k,izv
+      INTEGER,dimension(ndiffions) :: itrans
+      INTEGER,dimension(1) :: etrans
+      REAL(kind=kind(1.D0)),dimension(nlayer) :: P1d,Tn1d,M1d
+      REAL(kind=kind(1.D0)),dimension(nlayer+1) :: P1di
+      REAL(kind=kind(1.D0)),dimension(nlayer,ndiffions) :: Qk1d,Qk1dnew
+      REAL(kind=kind(1.D0)),dimension(nlayer,ndiffions) :: TI1d,TI1dnew 
+      REAL(kind=kind(1.D0)),dimension(nlayer) :: TE1d,TE1Dnew
+      REAL(kind=kind(1.D0)),dimension(nlayer,ndiffions) :: Rknew
+      REAL(kind=kind(1.D0)),dimension(nlraf) :: Zraf,Praf,Traf,Mraf
+      REAL(kind=kind(1.D0)),dimension(nlraf,ndiffions) :: QIraf,TIraf
+      REAL(kind=kind(1.D0)),dimension(nlraf) :: TEraf
+      REAL(kind=kind(1.D0)),dimension(nlraf,ndiffions) :: RIraf
+      REAL(kind=kind(1.D0)) :: masseU,kbolt,Konst,g
+      REAL(kind=kind(1.D0)) :: facP,rhoNloc
+      masseU=1.e-3/RNAVO
+      kbolt=RKBOL
+      Konst=kbolt/masseU
+      g=RG
+      
+      
+      ! below layer iz_vertplasma no effect due to vertical diffusion
+      do l=1,nlayer
+        if (P1D(l) .ge. Praf(1)) then
+          qk1dnew(l,:)=qk1d(l,:)
+          Ti1dnew(l,:)=Ti1d(l,:)
+          Te1dnew(l)=TE1d(l)
+        endif
+          
+        if (P1D(l) .lt. Praf(1)) then
+          iP=1
+          do while(Praf(iP) .gt. P1D(l))
+            iP=iP+1
+          enddo
+          IP=iP-1
+          facP=(P1D(l)-Praf(iP))/(Praf(iP+1)-Praf(iP))
+          rhoNloc=P1D(l)/kbolt/TN1D(l)*M1d(l)*masseU
+          
+          do k=1,ndiffions
+            Rknew(l,k)=RIraf(iP,k)+(RIraf(iP+1,k)-RIraf(iP,k))*facP
+            TI1dnew(l,k)=TIraf(IP,k)+(TIraf(iP+1,k)-TIraf(iP,k))*facP
+          enddo
+          TE1Dnew(l)=TEraf(iP)+(TEraf(iP+1)-TEraf(iP))*facP
+          
+          do k=1,ndiffions
+            Qk1dnew(l,k)=Rknew(l,k)/rhoNloc
+          enddo
+          
+        endif
+      enddo
+       
+      END
+      
+      SUBROUTINE GCMGRID_P2K(Zraf,Praf,Traf,Mraf,QIraf,RIraf,
+     &  Tiraf,Teraf,Wk,Ez,P1Di,P1D,Tn1d,M1d,qk1d,qk1dnew,
+     &  TI1D,TI1Dnew,TE1D,TE1Dnew,Wk1D,Wk1D2,Wk1D3,Ez1d,ZZ1d,
+     &  nlraf,ndiffions,nlayer,izv,itrans,etrans,FacQ,FacTI,FacTe)
+      
+      IMPLICIT NONE
+#include "YOMCST.h"    
+      INTEGER :: nlraf,ndiffions,nlayer,l,iP,k,izv
+      INTEGER :: itrans(ndiffions)
+      INTEGER :: etrans(1)
+      REAL(kind=kind(1.D0)),dimension(nlayer) :: P1d,Tn1d,M1d,ZZ1d
+      REAL(kind=kind(1.D0)),dimension(nlayer+1) :: P1di
+      REAL(kind=kind(1.D0)),dimension(nlayer,ndiffions) :: Qk1d,Qk1dnew
+      REAL(kind=kind(1.D0)),dimension(nlayer,ndiffions) :: TI1d,TI1dnew
+      REAL(kind=kind(1.D0)),dimension(nlayer) :: TE1D,TE1Dnew
+      REAL(kind=kind(1.D0)),dimension(nlayer+1,ndiffions) :: Wk1D,WK1d2
+      REAL(kind=kind(1.D0)),dimension(nlayer+1,ndiffions) :: Wk1D3
+      REAL(kind=kind(1.D0)),dimension(nlayer+1) :: Ez1D
+!      REAL(kind=kind(1.D0)),dimension(nlayer,ndiffions) :: Rknew
+      REAL(kind=kind(1.D0)),dimension(nlayer,ndiffions) :: FacQ,FacTI
+      REAL(kind=kind(1.D0)),dimension(nlayer) :: FacTE
+      REAL(kind=kind(1.D0)),dimension(nlraf) :: Zraf,Praf,Traf
+      REAL(kind=kind(1.D0)),dimension(nlraf) :: Mraf
+      REAL(kind=kind(1.D0)),dimension(nlraf,ndiffions) :: QIraf
+      REAL(kind=kind(1.D0)),dimension(nlraf,ndiffions) :: TIraf
+      REAL(kind=kind(1.D0)),dimension(nlraf) :: TEraf
+      REAL(kind=kind(1.D0)),dimension(nlraf,ndiffions) :: RIraf
+      REAL(kind=kind(1.D0)),dimension(nlraf,ndiffions) :: Wk
+      REAL(kind=kind(1.D0)),dimension(nlraf) :: Ez
+      REAL(kind=kind(1.D0)) masseU,kbolt,Konst,g
+      REAL(kind=kind(1.D0)) facP,rhoNloc
+      masseU=1.e-3/RNAVO
+      kbolt=RKBOL
+      Konst=kbolt/masseU
+      g=RG
+      
+      
+! below layer iz_vertplasma no effect due to vertical diffusion
+      do l=1,nlayer
+        if (P1D(l) .ge. Praf(1)) then
+        qk1dnew(l,1:ndiffions)=qk1d(l,1:ndiffions)
+        Ti1dnew(l,1:ndiffions)=Ti1d(l,1:ndiffions)
+        endif
+        
+        if (P1di(l) .ge. Praf(1)) then
+          Ez1D(l)=0D0
+          Wk1D(l,1:ndiffions)=0D0
+          Wk1D2(l,1:ndiffions)=0D0
+          Wk1D3(l,1:ndiffions)=0D0
+        endif
+        
+        if (P1D(l) .lt. Praf(1)) then
+          iP=1
+          do while(Praf(iP) .gt. P1D(l))
+            iP=iP+1
+          enddo
+          iP=iP-1
+          
+!          indP=maxloc(Praf,mask=Praf < P1D(l))
+!          IP=indP(1)-1
+          
+          facP=(P1D(l)-Praf(iP))/(Praf(iP+1)-Praf(iP))
+          rhoNloc=P1D(l)/kbolt/TN1D(l)*M1d(l)*masseU
+          do k=1,ndiffions
+             Qk1dnew(l,k)=(RIraf(iP,k)+(RIraf(iP+1,k)-RIraf(iP,k))*facP)
+     &  *facQ(l,k)/rhoNloc
+             TI1dnew(l,k)=(TIraf(IP,k)+(TIraf(iP+1,k)-TIraf(iP,k))*facP)
+     &  *facTI(l,k)
+          enddo
+          TE1Dnew(l)=(TEraf(iP)+(TEraf(iP+1)-TEraf(iP))*facP)*facTE(l)
+          
+!      if (l .eq. nlayer) then
+!        print*,'iP',iP,Praf(iP),P1D(l),Praf(iP+1),nlraf,facP
+!          print*,RIraf(iP,1),RIraf(iP+1,1),facQ(l,1),rhoNloc,facTI(l,1)
+!          print*,Qk1Dnew(l,1),TI1Dnew(l,1)
+!        endif
+        
+        endif
+        
+        if (P1Di(l) .lt. Praf(1)) then
+          iP=1
+          do while(Praf(iP) .gt. P1Di(l))
+            iP=iP+1
+          enddo
+          iP=iP-1
+          
+          facP=(P1Di(l)-Praf(iP))/(Praf(iP+1)-Praf(iP))
+          do k=1,ndiffions
+            wk1d(l,k)=(Wk(iP,k)+(Wk(iP+1,k)-Wk(iP,k))*facP)
+     &  /(ZZ1D(l)-ZZ1D(l-1))
+            wk1d2(l,k)=(Wk(iP,k)+(Wk(iP+1,k)-Wk(iP,k))*facP)
+     &  *(ZZ1D(l)-ZZ1D(l-1))
+            wk1d3(l,k)=(Wk(iP,k)+(Wk(iP+1,k)-Wk(iP,k))*facP)
+          enddo
+          Ez1D(l)=(EZ(iP)+(Ez(iP+1)-Ez(iP))*facP)
+        endif
+      
+      enddo ! l
+      Wk1d(nlayer+1,1:ndiffions)=0D0
+      Wk1d2(nlayer+1,1:ndiffions)=0D0
+      Wk1d3(nlayer+1,1:ndiffions)=0D0
+      Ez1D(nlayer+1)=0D0 
+      
+      
+      END
+      
+      
+      SUBROUTINE CORRFACK(qbef,qaft,Tbef,Taft,Tebef,Teaf,
+     &  Fq,Ft,Fte,nl,nq)
+     
+      IMPLICIT NONE
+    
+      INTEGER :: nl,nq,l,k
+      REAL(kind=kind(1.D0)),dimension(nl,nq) :: qbef,qaft,Fq
+      REAL(kind=kind(1.D0)),dimension(nl,nq) :: Tbef,Taft,Ft
+      REAL(kind=kind(1.D0)),dimension(nl) :: Tebef,Teaf,Fte
+      
+      do l=1,nl
+       do k=1,nq
+          Fq(l,k)=qbef(l,k)/qaft(l,k)
+          Ft(l,k)=Tbef(l,k)/Taft(l,k)
+       enddo
+       Fte(l)=Tebef(l)/Teaf(l)
+      enddo
+      
+      END
+      
+      SUBROUTINE DCOEFFK(TI,NUI,M,DI,nl,nqdyn,nqtot,ind)
+      
+      IMPLICIT NONE
+      
+      INTEGER :: nl,l,nqdyn,nqtot,k,kn
+      INTEGER,dimension(nqdyn) :: ind
+      REAL(kind=kind(1.D0)),dimension(nl,nqdyn) :: NUI,TI,DI
+      REAL,dimension(nqtot) :: M
+      REAL(kind=kind(1.D0)) :: kbolt,masseU
+      masseU=1.660538782d-27
+      kbolt=1.3806504d-23
+      
+      do k=1,nqdyn
+        kn=ind(k)
+        do l=1,nl
+          DI(l,k)=kbolt/masseU*TI(l,k)/NUI(l,k)/M(kn)
+        enddo
+      enddo
+      
+      END
+      
+      SUBROUTINE DIFFPARAMK(T,Te,RE,D,Dz,
+     &  alpha,beta,gama,delta,eps,dzeta,eta,nl,Wmax)
+     
+      IMPLICIT NONE
+      
+      INTEGER :: nl,l
+      REAL(kind=kind(1.D0)),DIMENSION(nl) :: T,D
+      REAL(kind=kind(1.D0)),DIMENSION(nl) :: Te,Re
+      REAL(kind=kind(1.D0)) :: DZ,Wmax
+      REAL(kind=kind(1.D0)),DIMENSION(nl) :: alpha,beta,gama,delta,eps
+      REAL(kind=kind(1.D0)),DIMENSION(nl) :: dzeta,eta
+      
+      alpha(1)=1D0/T(1)*(-3D0*T(1)+4D0*T(2)-T(3))/(2D0*DZ)+
+     &  1D0/T(1)*(-3D0*Te(1)+4D0*Te(2)-Te(3))/(2D0*DZ)
+      delta(1)=(-3D0*D(1)+4D0*D(2)-D(3))/(2D0*Dz)
+      beta(1)=1D0/T(1)
+      dzeta(1)=Te(1)/T(1)*
+     & (-3D0*log(Re(1))+4D0*log(Re(2))-log(Re(3)))/(2D0*DZ)
+      
+! Add a limit to dzeta (test to avoid unstabilities)
+      
+!      if (abs(dzeta(1)) .ge. 1d0 .and. Wmax .ge. 500d0) then
+!        dzeta(1)=1d0*dzeta(1)/abs(dzeta(1))
+!      endif
+      
+!      dzeta(1)=Te(1)/T(1)/Re(1)*
+!     & (-3D0*Re(1)+4D0*Re(2)-Re(3))/(2D0*Dz)
+      
+      do l=2,nl-1 
+        alpha(l)=1D0/T(l)*(T(l+1)-T(l-1))/(2D0*DZ)+
+     &  1D0/T(l)*(Te(l+1)-Te(l-1))/(2D0*DZ)
+      delta(l)=(D(l+1)-D(l-1))/(2D0*Dz)
+      beta(l)=1D0/T(l)
+      dzeta(l)=Te(l)/T(l)*(log(Re(l+1))-log(Re(l-1)))/(2D0*DZ)
+!      dzeta(l)=Te(l)/T(l)/Re(l)*(Re(l+1)-Re(l-1))/(2D0*DZ)
+      
+!      if (abs(dzeta(l)) .ge. 1d0 .and. Wmax .ge. 500d0) then
+!        dzeta(l)=1d0*dzeta(l)/abs(dzeta(l))
+!      endif
+      
+      enddo
+      
+! Upper altitude values
+      
+      alpha(nl)=1D0/T(nl)*(3D0*T(nl)-4D0*T(nl-1)+T(nl-2))/(2D0*DZ)+
+     &  1D0/T(nl)*(3D0*Te(nl)-4D0*Te(nl-1)+Te(nl-2))/(2D0*DZ)
+        delta(nl)=(3D0*D(nl)-4D0*D(nl-1)+D(nl-2))/(2D0*DZ)
+        beta(nl)=1D0/T(nl)
+      dzeta(nl)=Te(nl)/T(nl)*
+     & (3D0*log(Re(nl))-4D0*log(Re(nl-1))+log(Re(nl-2)))/(2D0*DZ)
+      
+!      if (abs(dzeta(nl)) .ge. 1d0 .and. Wmax .ge. 500d0) then
+!        dzeta(nl)=1d0*dzeta(nl)/abs(dzeta(nl))
+!      endif
+      
+!      dzeta(nl)=Te(nl)/T(nl)/Re(nl)*
+!     & (3D0*Re(nl)-4D0*Re(nl-1)+Re(nl-2))/(2D0*Dz)
+      
+! Compute the gama and eps coefficients
+! Lower altitude values
+      
+      gama(1)=(-3D0*beta(1) +4D0*beta(2) -beta(3) )/(2d0*dz)
+      eps(1) =(-3D0*alpha(1)+4D0*alpha(2)-alpha(3))/(2d0*dz)
+      eta(1) =(-3D0*dzeta(1)+4D0*dzeta(2)-dzeta(3))/(2D0*dz)
+      
+        do l=2,nl-1
+          gama(l)=(beta(l+1) -beta(l-1) )/(2d0*Dz)
+          eps(l) =(alpha(l+1)-alpha(l-1))/(2d0*Dz)
+          eta(l) =(dzeta(l+1)-dzeta(l-1))/(2D0*Dz)
+        end do
+      
+        gama(nl)=(3D0*beta(nl) -4D0*beta(nl-1) +beta(nl-2) )/(2D0*DZ)
+        eps(nl) =(3D0*alpha(nl)-4D0*alpha(nl-1)+alpha(nl-2))/(2D0*DZ)
+        eta(nl) =(3D0*dzeta(nl)-4D0*dzeta(nl-1)+dzeta(nl-2))/(2D0*dz)
+!        print*,'test',alpha(nl-1),beta(nl-1),dzeta(nl-1)
+      
+      
+        END
+      
+      
+      SUBROUTINE MATCOEFFK(alpha,beta,gama,delta,eps,dzeta,eta,
+     &  Dad,Rhoad,Facesc,dz,dt,A,B,C,D,Ti,Te,nl)
+      
+      IMPLICIT NONE
+      
+      INTEGER :: nl,l
+      REAL*8,DIMENSION(nl) :: alpha,beta,gama,delta,eps,dzeta,eta
+      REAL*8,DIMENSION(nl) :: Dad,RHoad,Ti,Te
+      REAL*8 :: dz,dt,del1,del2,del3,FacEsc
+      REAL*8,DIMENSION(nl) :: A,B,C,D
+      del1=dt/2d0/dz
+      del2=dt/dz/dz
+      del3=dt
+      
+! lower boundary coefficients no change
+      A(1)=0d0
+      B(1)=1d0
+      C(1)=0d0
+      D(1)=rhoAd(1)
+      
+      do l=2,nl-1
+        A(l)=(delta(l)+(alpha(l)+beta(l)+dzeta(l))*Dad(l))*del1
+     &       -Dad(l)*del2
+        B(l)=-(delta(l)*(alpha(l)+beta(l)+dzeta(l))
+     &       +Dad(l)*(gama(l)+eps(l)+eta(l)))*del3
+        B(l)=B(l)+1d0+2d0*Dad(l)*del2
+        C(l)=-(delta(l)+(alpha(l)+beta(l)+dzeta(l))*Dad(l))*del1
+     &       -Dad(l)*del2
+        D(l)=rhoAD(l)
+      enddo
+      
+! Upper boundary coefficients Diffusion profile
+      C(nl)=0d0
+      B(nl)=-1d0
+      A(nl)=exp(-dZ*(alpha(nl-1)+beta(nl-1)+dzeta(nl-1)))*FacEsc
+!      A(nl)=exp(-dZ*(alpha(nl)+beta(nl)+dzeta(nl)))*FacEsc
+!      A(nl)=exp(-dZ*(alpha(nl-1)+beta(nl-1)))*FacEsc 
+!      A(nl)=exp(-dZ*(alpha(nl-1)+beta(nl-1))*
+!     &  (Ti(nl-1)/(Te(nl-1)+Ti(nl-1))))*Facesc
+      D(nl)=0D0
+      
+      END
+      
+      SUBROUTINE CheckmassK(X,Y,nl)
+      
+      IMPLICIT NONE
+      
+      INTEGER :: nl
+      REAL*8,DIMENSION(nl) :: X,Y
+      REAL*8 Xtot,Ytot
+      
+      Xtot=sum(X)
+      Ytot=sum(Y)
+      
+      if (abs((Xtot-Ytot)/Xtot) .gt. 1d-4) then
+        print*,'no conservation for mass',Xtot,Ytot
+      endif
+      
+      END
+      
+      
Index: /trunk/LMDZ.VENUS/libf/phyvenus/clesphys.h
===================================================================
--- /trunk/LMDZ.VENUS/libf/phyvenus/clesphys.h	(revision 3034)
+++ /trunk/LMDZ.VENUS/libf/phyvenus/clesphys.h	(revision 3035)
@@ -11,8 +11,9 @@
        LOGICAL ok_orodr,ok_orolf,ok_gw_nonoro
        LOGICAL ok_kzmin,tuneupperatm, ok_ionchem, ok_jonline
+       LOGICAL ok_iondiff
        LOGICAL callnlte,callnirco2,callthermos
        LOGICAL ok_cloud, ok_chem, reinit_trac, ok_sedim
        LOGICAL ok_clmain, physideal, startphy_file
-       INTEGER nbapp_rad, nbapp_chim, iflag_con, iflag_ajs
+       INTEGER nbapp_rad, nbapp_chem, iflag_con, iflag_ajs
        INTEGER lev_histins, lev_histday, lev_histmth
        INTEGER tr_scheme, cl_scheme
@@ -29,7 +30,8 @@
      &     tuneupperatm,callnlte,callnirco2,callthermos,                &
      &     ok_cloud, ok_chem, reinit_trac, ok_sedim,                    &
-     &     ok_clmain, physideal, startphy_file, ok_ionchem, ok_jonline
+     &     ok_clmain, physideal, startphy_file, ok_ionchem, ok_jonline, &
+     &     ok_iondiff
 
-       COMMON/clesphys_i/ nbapp_rad, nbapp_chim,                        &
+       COMMON/clesphys_i/ nbapp_rad, nbapp_chem,                        &
      &     iflag_con, iflag_ajs,                                        &
      &     lev_histins, lev_histday, lev_histmth, tr_scheme,            &
Index: /trunk/LMDZ.VENUS/libf/phyvenus/conf_phys.F90
===================================================================
--- /trunk/LMDZ.VENUS/libf/phyvenus/conf_phys.F90	(revision 3034)
+++ /trunk/LMDZ.VENUS/libf/phyvenus/conf_phys.F90	(revision 3035)
@@ -92,11 +92,11 @@
        call getin('nbapp_rad',nbapp_rad)
        print*,"nbapp_rad",nbapp_rad
-!Config  Key  = nbapp_chim
+!Config  Key  = nbapp_chem
 !Config  Desc = Frequence d'appel a la chimie
-!Config  Def  = 1
+!Config  Def  = 24000
 !Config  Help = Nombre  d'appels des routines de chimie
 !Config         par jour.
-       nbapp_chim = 1
-       call getin('nbapp_chim',nbapp_chim)
+       nbapp_chem = 24000
+       call getin('nbapp_chem',nbapp_chem)
 
 !Config  Key  = iflag_con
@@ -481,5 +481,5 @@
 !Config Desc =
 !Config Def  = .false.
-!Config Help =
+!Config Help = Activate thermosphere tuning for oxygen compensation
 !
   tuneupperatm = .false.
@@ -499,5 +499,5 @@
 !Config Desc =
 !Config Def  = .false.
-!Config Help =
+!Config Help = Activate the ionosphere chemistry
 !
   ok_ionchem = .false.
@@ -514,4 +514,24 @@
 
 !
+!Config Key  = ok_iondiff
+!Config Desc =
+!Config Def  = .false.
+!Config Help = Activate the ambipolar ion diffusion in physic
+!
+  ok_iondiff = .false.
+  call getin('ok_iondiff',ok_iondiff)
+
+  if ((ok_ionchem.eqv..false.).and.(.true..eqv.ok_iondiff)) then
+    write(*,*) "Error incoherent flags :"
+    write(*,*) "ok_ionchem=",ok_ionchem," / ok_iondiff=",ok_iondiff
+    write(*,*) "If you include ions diffusion,"
+    write(*,*) "   you need also ok_ionchem==.true."
+    write(*,*) "If you do not include ions diffusion,"
+    write(*,*) "   you need ok_iondiff==.false."
+    write(*,*) "Check physiq.def"
+    stop
+  endif
+
+!
 !
 !Config Key  = 
@@ -535,5 +555,5 @@
   write(lunout,*)' ok_gw_nonoro = ', ok_gw_nonoro
   write(lunout,*)' nbapp_rad = ', nbapp_rad
-  write(lunout,*)' nbapp_chim = ', nbapp_chim
+  write(lunout,*)' nbapp_chem = ', nbapp_chem
   write(lunout,*)' iflag_con = ', iflag_con
   write(lunout,*)' Sortie journaliere = ', ok_journe
@@ -576,5 +596,6 @@
   write(lunout,*)' ok_jonline = ',ok_jonline
   write(lunout,*)' ok_ionchem = ',ok_ionchem
-  
+  write(lunout,*)' ok_iondiff = ',ok_iondiff
+
   end subroutine conf_phys
 
Index: /trunk/LMDZ.VENUS/libf/phyvenus/nonoro_gwd_ran_mod.F90
===================================================================
--- /trunk/LMDZ.VENUS/libf/phyvenus/nonoro_gwd_ran_mod.F90	(revision 3034)
+++ /trunk/LMDZ.VENUS/libf/phyvenus/nonoro_gwd_ran_mod.F90	(revision 3035)
@@ -264,4 +264,6 @@
              zp(jw, ii) = (sign(1., 0.5 - ran_num_1) &
                         + 1.) * RPI / 2.
+           !   zp(jw, ii) = atan2(4.*vh(ii, launch),uh(ii, launch)) &
+           !            + (sign(1., 0.5 - ran_num_1) + 1.) * RPI / 2.
            ! --------- TRY 00 -----------------------
            !IF((abs(latitude_deg(ii)) .le. 55.)) THEN
@@ -301,8 +303,12 @@
                      sqrt(3.) / sqrt(na * 1.)
            end DO
+           !cpha = cpha - uh(ii, launch)
            IF (cpha < 0.) THEN
               cpha = - 1. * cpha
               zp (jw, ii) = zp(jw, ii) + RPI
            ENDIF
+           !IF (cpha < 1.) THEN
+           !  cpha = 1.
+           !ENDIF
 ! otherwise, with the computation below, we get a uniform distribution between cmin and cmax.
 !           ran_num_3 = mod(tt(ii, jw)**2, 1.)
Index: /trunk/LMDZ.VENUS/libf/phyvenus/photochemistry_venus.F90
===================================================================
--- /trunk/LMDZ.VENUS/libf/phyvenus/photochemistry_venus.F90	(revision 3034)
+++ /trunk/LMDZ.VENUS/libf/phyvenus/photochemistry_venus.F90	(revision 3035)
@@ -2507,5 +2507,5 @@
 !-- TuneE
 ! VCD 2.0 tuning
-    v_phot(65:90,4) = v_phot(65:90,4)*10. ! CO2 + hv ==> O(1D) + CO
+    v_phot(65:nz,4) = v_phot(65:nz,4)*10. ! CO2 + hv ==> O(1D) + CO
 !--
 !   v_phot(:,4) = v_phot(:,4)*10.
@@ -4590,4 +4590,7 @@
           v_phot(65:nz,4) = v_phot(65:nz,4)*8. ! CO2 + hv ==> O(1D) + CO
       !--
+      ! TEST TUNEH
+      !     v_phot(62:74,4) = v_phot(62:74,4)*6.5 ! CO2 + hv ==> O(1D) + CO
+      !--
       !   v_phot(:,4) = v_phot(:,4)*10.
       !do ij=3,4
Index: /trunk/LMDZ.VENUS/libf/phyvenus/physiq_mod.F
===================================================================
--- /trunk/LMDZ.VENUS/libf/phyvenus/physiq_mod.F	(revision 3034)
+++ /trunk/LMDZ.VENUS/libf/phyvenus/physiq_mod.F	(revision 3035)
@@ -298,5 +298,5 @@
 
       integer :: chempas
-      real :: nbapp_chem, zctime
+      real :: zctime
 
 ! sedimentation
@@ -340,4 +340,7 @@
 c Tendencies due to molecular diffusion
       real d_q_moldif(klon,klev,nqmax)
+
+c Tendencies due to ambipolar ion diffusion
+      real d_q_iondif(klon,klev,nqmax)
 
 c
@@ -671,4 +674,14 @@
          ENDIF
 
+c Verification synchronize AMBIPOLAR DIFFUSION & CHEMISTRY
+
+         IF ((ok_iondiff) .and. (NINT(RDAY/dtime).ne.nbapp_chem)) THEN
+          WRITE(lunout,*)'nbapp_chem .NE. day_step'
+          WRITE(lunout,*)'nbapp_chem ', nbapp_chem
+          WRITE(lunout,*)'day_step ', NINT(RDAY/dtime)
+          WRITE(lunout,*)'nbapp_chem must be equal to day_step' 
+          STOP
+         ENDIF
+        
 c Initialisation des sorties 
 c===========================
@@ -1050,5 +1063,5 @@
 !=========
 
-         nbapp_chem = 24000
+         
          chempas = nint(rday/pdtphys/nbapp_chem)
          zctime = dtime*real(chempas)             ! chemical timestep
@@ -1659,4 +1672,37 @@
 
 c====================================================================
+
+c==================================
+!  --  ION AMBIPOLAR DIFFUSION ---
+c==================================
+
+         d_q_iondif(:,:,:)=0
+
+         IF (callthermos .and. ok_chem .and. ok_ionchem) THEN
+           IF (ok_iondiff) THEN
+
+            call iondiff_red(klon,klev,nqmax,pplay,paprs,
+     &                        t_seri,tr_seri,pphis,
+     &                        gmtime,latitude_deg,longitude_deg,
+     &                        pdtphys,d_t_euv,d_t_conduc,d_q_moldif,
+     &                        d_q_iondif)
+
+            !write (*,*) 'TITI EST PASSE PAR LA'
+! --- update tendencies tracers ---
+
+            DO iq = 1, nqmax
+              IF (type_tr(iq) .eq. 2) THEN
+                DO k=1,klev
+                  DO ig=1,klon
+                    tr_seri(ig,k,iq)= max(tr_seri(ig,k,iq)+
+     &                           d_q_iondif(ig,k,iq)*dtime,1.e-30)
+                  ENDDO
+                ENDDO
+              ENDIF
+            ENDDO
+           ENDIF ! ok_iondiff
+         ENDIF  ! callthermos & ok_chem & ok_ionchem
+
+c====================================================================
       ! tests: output tendencies
 !      call writefield_phy('physiq_dtrad',dtrad,klev)
Index: /trunk/LMDZ.VENUS/libf/phyvenus/phytrac_chimie.F
===================================================================
--- /trunk/LMDZ.VENUS/libf/phyvenus/phytrac_chimie.F	(revision 3034)
+++ /trunk/LMDZ.VENUS/libf/phyvenus/phytrac_chimie.F	(revision 3035)
@@ -391,8 +391,16 @@
                         
 !           solar zenith angle
-
-            sza_local = acos(cos(lat_local(ilon))*cos(lon_local(ilon))
+!            sza_local = acos(cos(lat_local(ilon))*cos(lon_local(ilon))
+!     $                 *cos(lon_sun) + cos(lat_local(ilon))
+!     $                 *sin(lon_local(ilon))*sin(lon_sun))*180./rpi
+
+            sza_local = cos(lat_local(ilon))*cos(lon_local(ilon))
      $                 *cos(lon_sun) + cos(lat_local(ilon))
-     $                 *sin(lon_local(ilon))*sin(lon_sun))*180./rpi 
+     $                 *sin(lon_local(ilon))*sin(lon_sun)
+
+            ! Security - Handle rare cases where |sza_local| > 1
+            sza_local = min(sza_local,1.)
+            sza_local = max(-1.,sza_local)
+            sza_local = acos(sza_local)*180./rpi
 
 !           electron temperature 
Index: /trunk/LMDZ.VENUS/libf/phyvenus/plasmaphys_venus_mod.F90
===================================================================
--- /trunk/LMDZ.VENUS/libf/phyvenus/plasmaphys_venus_mod.F90	(revision 3035)
+++ /trunk/LMDZ.VENUS/libf/phyvenus/plasmaphys_venus_mod.F90	(revision 3035)
@@ -0,0 +1,88 @@
+MODULE plasmaphys_venus_mod
+
+!       COMMON/com_plasma/
+ 
+IMPLICIT NONE
+
+!       LOGICAL :: plasma 
+!       LOGICAL :: isoT      
+       !INTEGER :: ndiffions    ! nombre d'ions decrit
+       INTEGER :: nsubstep     ! nombre de sous pas de temps pour le calcul des tendances dynamiques
+       INTEGER :: naltvert     ! nombre de points sur la grille altitude pour dynamique verticale
+       INTEGER :: nsubvert     ! nombre de sous pas de temps maximal pour dynamique verticale
+       INTEGER :: nsubmin      ! nombre de sous pas de temps minimal pour dynamique verticale
+!       INTEGER :: iappvert     ! periode appel de dynamique verticale
+       !INTEGER :: iz_vertplasma! vertical index above which ion diffusion is computed
+!       INTEGER :: irun         ! numero de simulation
+!       INTEGER nspecies     ! number of molecular species (different from nqmx if active dust included)
+!       REAL :: qcharge,muvide,preffplasma,masseU,melec                     &
+!     & ,dtplasma,Freqklim,kbolt
+       REAL*8 :: Pdiffion 
+       REAL*8 :: qcharge,muvide
+
+!       parameter (ndiffions = 8)
+!       parameter (Pdiffion  = 7.D-4) ! pressure in Pa below which ion diffusion is computed
+!        ORIGINE -------------------------------
+       !parameter (naltvert  = 28)
+       !parameter (nsubvert  = 12000)   ! ptimestep/nsubvert
+       !parameter ( nsubmin  = 80)    ! ptimestep/nsubmin
+!        STABILITY MORE5  ----------------------------------
+       !parameter (Pdiffion  = 7.0D-4) ! pressure in Pa below which ion diffusion is computed
+       !parameter (naltvert  = 28*3)
+       !parameter (nsubvert  = 12000)   ! ptimestep/nsubvert
+       !parameter ( nsubmin  = 40)    ! ptimestep/nsubmin
+!        STABILITY MORE6  ----------------------------------
+       !parameter (Pdiffion  = 7.0D-4) ! pressure in Pa below which ion diffusion is computed
+       !parameter (naltvert  = 28*3*2)
+       !parameter (nsubvert  = 12000*3)   ! ptimestep/nsubvert
+       !parameter ( nsubmin  = 40)    ! ptimestep/nsubmin
+!       CURRENT BEST TIME CALCULATION  ----------------------------------
+       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
+
+
+       !parameter (iappvert  = 500)
+       parameter (qcharge=1.602176487D-19)
+       !parameter (kbolt=1.3806504d-23)
+       parameter (muvide=1.2566371D-6)
+       !parameter (preffplasma=1.e-10)
+       !parameter (masseU=1.660538782e-27)
+       !parameter (melec=9.1093826e-31) 
+!       parameter (isoT=.TRUE.)
+
+!! OLDER CODE BELOW ======================================================
+       
+! Ions parameter 
+        
+       !INTEGER itrans(ndynions)                  ! global index of fluid ions
+       !INTEGER etrans(1)                         ! index of electron    
+       !INTEGER,DIMENSION(:),ALLOCATABLE :: iplasma ! = 0 ou 1 si traceur ou fluide 
+       !CHARACTER(LEN=20) ionstrans(ndynions)
+       !CHARACTER(LEN=20) electrans(1)  
+       !REAL,DIMENSION(:),ALLOCATABLE ::  Mtraceur
+       !REAL, DIMENSION(:),ALLOCATABLE :: Charge
+       !REAL Cpkabs(ndynions),cpk(ndynions),kappak(ndynions)
+       !REAL cpeabs,cpe,kappae
+
+!	parameter (ionstrans=["o2plus","oplus"])
+!	parameter (electrans=["elec"])
+
+!
+!! Dissipation parameters
+!
+!!       REAL tetaudivk,tetaurotk,tetahk,cdivuk,crotk,cdivhk	
+!        REAL,DIMENSION(:),ALLOCATABLE :: tetaudivk,tetaurotk,tetahk
+!
+!!       COMMON/com_plasma/                                               &
+!!     & ionstrans,electrans,irun,                                        &
+!!     & plasma,itrans,etrans,iplasma,Mtraceur,dtplasma,Charge,           &
+!!     & cpkabs,cpk,kappak,cpeabs,cpe,kappae
+!
+!!       COMMON/diss_plasma/tetaudivk(llp),tetaurotk(llp),tetahk(llp)
+!       	
+
+END MODULE
+
+!-----------------------------------------------------------
