Index: trunk/LMDZ.VENUS/libf/phyvenus/clesphys.h
===================================================================
--- trunk/LMDZ.VENUS/libf/phyvenus/clesphys.h	(revision 1121)
+++ trunk/LMDZ.VENUS/libf/phyvenus/clesphys.h	(revision 1160)
@@ -10,17 +10,21 @@
        LOGICAL cycle_diurne,soil_model
        LOGICAL ok_orodr,ok_orolf,ok_gw_nonoro
+       LOGICAL ok_kzmin
        INTEGER nbapp_rad, nbapp_chim, iflag_con, iflag_ajs
+       INTEGER lev_histhf, lev_histday, lev_histmth
+       INTEGER tr_scheme
        REAL    ecriphy
        REAL    solaire
        REAL    z0, lmixmin
        REAL    ksta, inertie
-       LOGICAL ok_kzmin
-       INTEGER lev_histhf, lev_histday, lev_histmth
 
-       COMMON/clesphys/cycle_diurne, soil_model,                        &
-     &     ok_orodr, ok_orolf, ok_gw_nonoro, nbapp_rad, nbapp_chim      &
-     &     , ecriphy                                                    &
-     &     , iflag_con, iflag_ajs, solaire, z0, lmixmin, ksta           &
-     &     , ok_kzmin, lev_histhf, lev_histday, lev_histmth             &
-     &     , inertie
+       COMMON/clesphys_l/cycle_diurne, soil_model,                      &
+     &     ok_orodr, ok_orolf, ok_gw_nonoro, ok_kzmin
 
+       COMMON/clesphys_i/nbapp_rad, nbapp_chim,                         &
+     &     iflag_con, iflag_ajs,                                        &
+     &     lev_histhf, lev_histday, lev_histmth, tr_scheme
+
+       COMMON/clesphys_r/ecriphy, solaire, z0, lmixmin,                 &
+     &     ksta, inertie
+
Index: trunk/LMDZ.VENUS/libf/phyvenus/conf_phys.F90
===================================================================
--- trunk/LMDZ.VENUS/libf/phyvenus/conf_phys.F90	(revision 1121)
+++ trunk/LMDZ.VENUS/libf/phyvenus/conf_phys.F90	(revision 1160)
@@ -190,4 +190,25 @@
     call getin('solaire', solaire)
 !
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+! PARAMETER FOR THE TRACERS
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!
+!Config Key  = tr_scheme
+!Config Desc =
+!Config Def  = 0
+!Config Help =
+!
+! 0   = Nothing is done (passive tracers)
+! 1   = pseudo-chemistry with relaxation toward fixed profile
+!       See Marcq&Lebonnois 2013
+! 2   = surface emission
+!       For the moment, inspired from Mars version
+!       However, the variable 'source' could be used in physiq
+!       so the call to phytrac_emiss could be to initialise it.
+! 3   = Full chemistry
+!       To be added by Aurelien
+  tr_scheme = 0
+  call getin('tr_scheme',tr_scheme)
+
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 ! PARAMETER FOR THE PLANETARY BOUNDARY LAYER AND SOIL
@@ -323,4 +344,5 @@
   write(numout,*)' Equinoxe = ',R_peri
   write(numout,*)' Inclinaison =',R_incl
+  write(numout,*)' tr_scheme = ', tr_scheme
   write(numout,*)' iflag_pbl = ', iflag_pbl
   write(numout,*)' z0 = ',z0 
Index: trunk/LMDZ.VENUS/libf/phyvenus/physiq.F
===================================================================
--- trunk/LMDZ.VENUS/libf/phyvenus/physiq.F	(revision 1121)
+++ trunk/LMDZ.VENUS/libf/phyvenus/physiq.F	(revision 1160)
@@ -669,14 +669,27 @@
 
       if (iflag_trac.eq.1) then
-      call phytrac (
-     I                   itap, gmtime,
-     I                   debut,lafin,
-     I                   nqmax,
-     I                   nlon,nlev,dtime,
-     I                   u,v,t,paprs,pplay,
-     I                   rlatd,
-     I                   rlond,presnivs,pphis,pphi,
-     I                   falbe,
+
+       if (tr_scheme.eq.1) then
+! Case 1: pseudo-chemistry with relaxation toward fixed profile
+         call phytrac_relax (debut,lafin,nqmax,
+     I                   nlon,nlev,dtime,pplay,
      O                   tr_seri)
+
+       elseif (tr_scheme.eq.2) then
+! Case 2: surface emission
+! For the moment, inspired from Mars version
+! However, the variable 'source' could be used in physiq
+! so the call to phytrac_emiss could be to initialise it.
+         call phytrac_emiss ( (rjourvrai+gmtime)*RDAY,
+     I                   debut,lafin,nqmax,
+     I                   nlon,nlev,dtime,paprs,
+     I                   rlatd,rlond,
+     O                   tr_seri)
+       elseif (tr_scheme.eq.3) then
+! Case 3: Full chemistry
+!        call phytrac_chem ( ?? )
+         print*,"Chemistry not yet implemented..." 
+         print*,"See Aurelien Stolzenbach" 
+       endif
       endif
 
Index: trunk/LMDZ.VENUS/libf/phyvenus/phytrac.F
===================================================================
--- trunk/LMDZ.VENUS/libf/phyvenus/phytrac.F	(revision 1121)
+++ 	(revision )
@@ -1,886 +1,0 @@
-!
-! $Header: /home/cvsroot/LMDZ4/libf/phylmd/phytrac.F,v 1.16 2006/03/24 15:06:23 lmdzadmin Exp $
-!
-c
-c
-      SUBROUTINE phytrac (nstep,
-     I                    gmtime,
-     I                    debutphy,
-     I                    lafin,
-     I                    nqmax,
-     I                    nlon,
-     I                    nlev, 
-     I                    pdtphys,
-     I                    u,
-     I                    v,
-     I                    t_seri,
-     I                    paprs,
-     I                    pplay,
-     I                    xlat,
-     I                    xlon,
-     I                    presnivs,
-     I                    pphis,
-     I                    pphi,
-     I                    albsol,
-     O                    tr_seri)
-
-c======================================================================
-c Auteur(s) FH
-c Objet: Moniteur general des tendances traceurs
-c
-cAA Remarques en vrac:
-cAA--------------------
-cAA 1/ le call phytrac se fait avec nqmax 
-c======================================================================
-      USE ioipsl
-      USE infotrac
-      USE control_mod
-      use dimphy
-      USE comgeomphy
-      IMPLICIT none
-#include "YOMCST.h"
-#include "dimensions.h"
-#include "clesphys.h" !///utile?
-#include "temps.h"
-#include "paramet.h"
-c======================================================================
-
-c Arguments:
-c
-c   EN ENTREE:
-c   ==========
-c
-c   divers:
-c   -------
-c
-      integer nlon  ! nombre de points horizontaux
-      integer nlev  ! nombre de couches verticales
-      integer nqmax ! nombre de traceurs auxquels on applique la physique
-      integer nstep  ! appel physique
-      integer nseuil ! numero du premier traceur non CV
-c      integer julien !jour julien
-c      integer itop_con(nlon)
-c      integer ibas_con(nlon)
-      real gmtime
-      real pdtphys  ! pas d'integration pour la physique (seconde)
-      real t_seri(nlon,nlev) ! temperature
-      real tr_seri(nlon,nlev,nqmax) ! traceur  
-      real u(nlon,nlev)
-      real v(nlon,nlev)
-      real albsol(nlon)  ! albedo surface
-      real paprs(nlon,nlev+1)  ! pression pour chaque inter-couche (en Pa)
-      real ps(nlon)  ! pression surface
-      real pplay(nlon,nlev)  ! pression pour le mileu de chaque couche (en Pa)
-      real pphi(nlon,nlev) ! geopotentiel
-      real pphis(nlon)
-      REAL xlat(nlon)       ! latitudes pour chaque point 
-      REAL xlon(nlon)       ! longitudes pour chaque point 
-      REAL presnivs(nlev)
-      logical debutphy       ! le flag de l'initialisation de la physique
-      logical lafin          ! le flag de la fin de la physique
-c      REAL flxmass_w(nlon,nlev)
-
-cAA ----------------------------
-cAA  VARIABLES LOCALES TRACEURS
-cAA ----------------------------
-cAA
-
-C les traceurs
-C
-      logical flagCO_OCS
-c===================
-c it--------indice de traceur
-c k,i---------indices long, vert
-c===================
-c Variables deja declarees dont on a besoin pour traceurs   
-c   k,i,it,tr_seri(klon,klev,nqmax),pplay(nlon,nlev), 
-      integer nqCO_OCS
-c      real pzero,gamma
-c      parameter (pzero=85000.)
-c      parameter (gamma=5000.)
-      REAL alpha
-      real deltatr(klon,klev,nqtot) ! ecart au profil de ref zprof
-      real,save,allocatable :: zprof(:,:)
-      real,save,allocatable :: tau(:,:) ! temps de relaxation vers le profil (s)
-c======================================================================
-c
-c Declaration des procedures appelees
-c
-c--modif convection tiedtke
-      INTEGER i, k, it
-      INTEGER iq, iiq
-      REAL delp(klon,klev)
-c--end modif
-c
-c Variables liees a l'ecriture de la bande histoire physique
-c
-c Variables locales pour effectuer les appels en serie
-c----------------------------------------------------
-c
-      REAL d_tr(klon,klev), d_trs(klon) ! tendances de traceurs 
-      REAL d_tr_cl(klon,klev,nqmax) ! tendance de traceurs  couche limite
-      REAL d_tr_cv(klon,klev,nqmax) ! tendance de traceurs  conv pour chq traceur
-C
-      character*20 modname
-      character*80 abort_message
-c
-c   Controles
-c-------------
-      logical first,couchelimite,convection
-      save first,couchelimite,convection
-c Olivia
-       data first,couchelimite,convection
-     s     /.true.,.false.,.false./
-
-      modname = 'phytrac'
-c======================================================================
-
-       if(first) then
-         allocate(zprof(klev,nqtot),tau(klev,nqtot))
-         first = .false.
-       endif
-
-         ps(:)=paprs(:,1)
-c TRACEURS TYPE CO ET OCS
-         flagCO_OCS = .true.
-      if (flagCO_OCS) then
-         nqCO_OCS   = 6
-      else
-         nqCO_OCS   = 0
-      endif  ! flagCO_OCS
-
-c---------
-c debutphy
-c---------
-         if (debutphy) then
-                 print*,"DEBUT PHYTRAC"
-C         
-c=============================================================
-c=============================================================
-c=============================================================
-c   Initialisation des traceurs
-c=============================================================
-c=============================================================
-c=============================================================
-c
-c=============================================================
-c=============================================================
-
-C=========================================================================
-C=========================================================================
-      if (flagCO_OCS) then
-c II) Declaration d'un profil vertical de traceur OK
-c 
-c zprof   = profil de rappel 
-c
-c 1 -> CO ; 2 -> OCS 
-c def des profils en log(a) = a * log(P) + b par morceaux, cf. pollack et al
-c tr_seri en ppm
-c (initialisation seulement si ceux-ci sont nuls)
-
-c ICI, ON UTILISE 3 CONSTANTES DE TEMPS DIFFERENTES POUR CHAQUE,
-c DONC TRACEURS 1 A 3 POUR CO ET 4 A 6 POUR OCS
-C=========================================================================
-
-
-c Constantes de rappel:
-
-       print*,"INIT TAU"
-       do k=1,klev
-         tau(k,1)=1.e6
-         tau(k,2)=1.e7
-         tau(k,3)=1.e8
-         tau(k,4)=1.e6
-         tau(k,5)=1.e7
-         tau(k,6)=1.e8
-       enddo
-
-c CO
-
-      do it=1,3 
-       print*,"INIT ZPROF ",tname(it)
-       do k=1,klev
-         zprof(k,it)=0.
-c pour l'instant, tau fixe, mais possibilite de le faire varier avec z
-        if (pplay(klon/2,k) >= 4.8e6) then
-           zprof(k,it)=14.
-        endif
-        if ((pplay(klon/2,k)<=4.8e6).and.(pplay(klon/2,k)>=1.9e6)) then
-           alpha=(log(pplay(klon/2,k))-log(1.9e6))/
-     .     (log(4.8e6)-log(1.9e6))
-           zprof(k,it)=20.*(14./20.)**alpha
-        endif
-        if ((pplay(klon/2,k)<=1.9e6).and.(pplay(klon/2,k)>=1.5e5)) then
-           alpha=(log(pplay(klon/2,k))-log(1.5e5))/
-     .     (log(1.9e6)-log(1.5e5))
-           zprof(k,it)=39.*(20./39.)**alpha
-        endif
-        if ((pplay(klon/2,k)<=1.5e5).and.(pplay(klon/2,k)>=1.1e4)) then
-           alpha=(log(pplay(klon/2,k))-log(1.1e4))/
-     .     (log(2.73e5)-log(1.1e4))
-           zprof(k,it)=50.*(39./50.)**alpha
-        endif
-        if ((pplay(klon/2,k)<=1.1e4).and.(pplay(klon/2,k)>=1.3e3)) then
-           alpha=(log(pplay(klon/2,k))-log(1.3e3))/
-     .     (log(1.1e4)-log(1.3e3))
-           zprof(k,it)=2.*(50./2.)**alpha
-        endif
-        if ((pplay(klon/2,k)<=1.3e3).and.(pplay(klon/2,k)>=2.4)) then
-           alpha=(log(pplay(klon/2,k))-log(2.4))/
-     .     (log(1.3e3)-log(2.4))
-           zprof(k,it)=1000.*(2./1000.)**alpha
-        endif
-        if (pplay(klon/2,k) <= 2.4) then
-           zprof(k,it)=1000.
-        endif
-       enddo 
-       print*,zprof(:,it)
- 
-c OCS
-       print*,"INIT ZPROF ",tname(it+3)
-       do k=1,klev
-         zprof(k,it+3)=0.
-         if (pplay(klon/2,k) >= 4.8e6) then
-           zprof(k,it+3)=30.
-         endif
-         if ((pplay(klon/2,k)<=4.8e6).and.(pplay(klon/2,k)>=9.4e5))
-     *   then
-           alpha=(log(pplay(klon/2,k))-log(9.4e5))/
-     *     (log(4.8e6)-log(9.4e5))
-           zprof(k,it+3)=20.*(30/20.)**alpha
-         endif
-         if ((pplay(klon/2,k)<=9.4e5).and.(pplay(klon/2,k)>=4.724e5))
-     *   then
-           alpha=(log(pplay(klon/2,k))-log(4.724e5))/
-     *     (log(9.4e5)-log(4.724e5))
-           zprof(k,it+3)=0.5*(20/0.5)**alpha
-         endif
-         if ((pplay(klon/2,k)<=4.724e5).and.(pplay(klon/2,k)>=1.1e4))
-     *   then
-           alpha=(log(pplay(klon/2,k))-log(1.1e4))/
-     *     (log(4.724e5)-log(1.1e4))
-           zprof(k,it+3)=0.005*(0.5/0.005)**alpha
-         endif
-         if (pplay(klon/2,k)<=1.1e4) then
-           zprof(k,it+3)=0.
-         endif
-       end do
-       print*,zprof(:,it+3)
-      enddo 
-
-c Initialisation du traceur s'il est nul:
-       do it=1,nqCO_OCS
-        if ((tr_seri(klon/2,1,it).eq.0.).and.
-     .      (tr_seri(klon/2,klev/2,it).eq.0.).and.
-     .      (tr_seri(klon/2,klev,it).eq.0.)) then
-         print*,"INITIALISATION DE ",tname(it)
-         do k=1,klev
-           do i=1,klon
-             tr_seri(i,k,it) = zprof(k,it)
-           enddo
-         enddo
-        endif
-       enddo
-
-C=========================================================================
-      endif  ! flagCO_OCS
-C=========================================================================
-C=========================================================================
-
-c-------------
-c fin debutphy
-c-------------
-      ENDIF  ! fin debutphy 
-
-c======================================================================
-      if (flagCO_OCS) then
-c Rappel vers un profil
-c======================================================================
-         do it=1,nqCO_OCS
-           do k=1,klev
-             do i=1,klon
-c VERIF
-           if (tr_seri(i,k,it).lt.0) then
-             print*,"Traceur negatif AVANT rappel:",i,k,it
-             stop
-           endif
-c FIN VERIF
-
-           deltatr(i,k,it) = (-tr_seri(i,k,it)+zprof(k,it))/tau(k,it)
-           tr_seri(i,k,it) =  tr_seri(i,k,it) + deltatr(i,k,it)*pdtphys
-
-c VERIF
-           if (tr_seri(i,k,it).lt.0) then
-             print*,"APRES rappel:",i,k,it,
-     .  deltatr(i,k,it),zprof(k,it),tr_seri(i,k,it),pdtphys/tau(k,it)
-             stop
-           endif
-c FIN VERIF
-             enddo
-           enddo
-         enddo
-
-c======================================================================
-      endif  ! flagCO_OCS
-c======================================================================
-
-c======================================================================
-c   Calcul de l'effet de la couche limite remis directement dans physiq
-c======================================================================
-
-
-      RETURN
-      END
-      
-      
-c=========================================================================
-c=========================================================================
-c=========================================================================
-c ARCHIVES ===============================================================
-c=========================================================================
-c=========================================================================
-c=========================================================================
-
-c===========
-c definition de traceurs idealises
-c==========
-c
-c I) Declaration directe du traceur a altitude fixee
-c
-c a) traceur en carre OK
-c
-c         do i=1,klon
-c         tr_seri(i,:,1)=0.
-c        if ((xlat(i)>=0.).and.(xlat(i)<=-30.)) then
-c          if ((xlon(i)>=0.).and.(xlon(i)<=40.)) then
-c              tr_seri(i,10,1)=1.
-c          endif
-c        endif
-c      end do 
-c
-c a bis) 2 traceurs en carre lat/alt, uniforme en longitude OK 
-c
-C entre 45-55 km
-c
-c         do i=1,klon
-c         do k=1,klev+1
-cc         tr_seri(i,k,1)=0.
-c           if ((xlat(i)>=60.).and.(xlat(i)<=80.)) then
-c           if ((xlon(i)>=-180.).and.(xlon(i)<=180.)) then
-c           if ((pplay(klon/2,k)>=5.e4).and.(pplay(klon/2,k)<=4.e5)) then
-c               tr_seri(i,k,1)=1.
-c           endif
-c           endif
-c           endif
-c         else
-c            tr_seri(i,k,1)=0.
-c         end do 
-c         end do
-cc
-c         do i=1,klon
-c         do k=1,klev+1
-cc         tr_seri(i,k,2)=0.
-c           if ((xlat(i)>=-60.).and.(xlat(i)<=-80.)) then
-c           if ((xlon(i)>=-180.).and.(xlon(i)<=180.)) then
-c           if ((pplay(klon/2,k)>=5.e4).and.(pplay(klon/2,k)<=4.e5)) then
-c               tr_seri(i,k,2)=1.
-c           endif
-c           endif
-c           endif
-c         else
-c            tr_seri(i,k,2)=0.
-c         end do 
-c         end do
-cc
-c         do i=1,klon
-c         do k=1,klev+1
-cc         tr_seri(i,k,3)=0.
-c           if ((xlat(i)>=40.).and.(xlat(i)<=60.)) then
-c           if ((xlon(i)>=-180.).and.(xlon(i)<=180.)) then
-c           if ((pplay(klon/2,k)>=5.e4).and.(pplay(klon/2,k)<=4.e5)) then
-c               tr_seri(i,k,3)=1.
-c           endif
-c           endif
-c           endif
-c         else
-c            tr_seri(i,k,3)=0.
-c         end do 
-c         end do
-cc
-c         do i=1,klon
-c         do k=1,klev+1
-cc         tr_seri(i,k,4)=0.
-c           if ((xlat(i)>=-40.).and.(xlat(i)<=-60.)) then
-c           if ((xlon(i)>=-180.).and.(xlon(i)<=180.)) then
-c           if ((pplay(klon/2,k)>=5.e4).and.(pplay(klon/2,k)<=4.e5)) then
-c               tr_seri(i,k,4)=1.
-c           endif
-c           endif
-c           endif
-c         else
-c            tr_seri(i,k,4)=0.
-c         end do 
-c         end do
-cc
-c         do i=1,klon
-c         do k=1,klev+1
-cc         tr_seri(i,k,5)=0.
-c           if ((xlat(i)>=-20.).and.(xlat(i)<=20.)) then
-c           if ((xlon(i)>=-180.).and.(xlon(i)<=180.)) then
-c           if ((pplay(klon/2,k)>=5.e4).and.(pplay(klon/2,k)<=4.e5)) then
-c              tr_seri(i,k,5)=1.
-c           endif
-c           endif
-c           endif
-c         else
-c            tr_seri(i,k,5)=0.
-c         end do 
-c         end do
-c
-c entre 35-45 km
-c
-c         do i=1,klon
-c         do k=1,klev+1
-cc         tr_seri(i,k,6)=0.
-c           if ((xlat(i)>=60.).and.(xlat(i)<=80.)) then
-c           if ((xlon(i)>=-180.).and.(xlon(i)<=180.)) then
-c           if ((pplay(klon/2,k)>=4.e5).and.(pplay(klon/2,k)<=8.e6)) then
-c               tr_seri(i,k,6)=1.
-c           endif
-c           endif
-c           endif
-c         else
-c            tr_seri(i,k,6)=0.
-c         end do 
-c         end do
-c
-c         do i=1,klon
-c         do k=1,klev+1
-cc         tr_seri(i,k,7)=0.
-c           if ((xlat(i)>=-60.).and.(xlat(i)<=-80.)) then
-c           if ((xlon(i)>=-180.).and.(xlon(i)<=180.)) then
-c           if ((pplay(klon/2,k)>=4.e5).and.(pplay(klon/2,k)<=8.e6)) then
-c               tr_seri(i,k,7)=1.
-c           endif
-c           endif
-c           endif
-c         else
-c            tr_seri(i,k,7)=0.
-c         end do 
-c         end do
-c
-C entre 50-60 km
-c
-c         do i=1,klon
-c         do k=1,klev+1
-cc         tr_seri(i,k,8)=0.
-c           if ((xlat(i)>=60.).and.(xlat(i)<=80.)) then
-c           if ((xlon(i)>=-180.).and.(xlon(i)<=180.)) then
-c           if ((pplay(klon/2,k)>=1.e4).and.(pplay(klon/2,k)<=1.e5)) then
-c               tr_seri(i,k,8)=1.
-c           endif
-c           endif
-c          endif
-c         else
-c            tr_seri(i,k,8)=0.
-c         end do 
-c         end do
-c
-c         do i=1,klon
-c         do k=1,klev+1
-cc         tr_seri(i,k,9)=0.
-c           if ((xlat(i)>=-80.).and.(xlat(i)<=-60.)) then
-c           if ((xlon(i)>=-180.).and.(xlon(i)<=180.)) then
-c           if ((pplay(klon/2,k)>=1.e4).and.(pplay(klon/2,k)<=1.e5)) then
-c               tr_seri(i,k,9)=1.
-c           endif
-c           endif
-c           endif
-c         else
-c            tr_seri(i,k,9)=0.
-c         end do 
-c         end do
-c
-c         do i=1,klon
-c         do k=1,klev+1
-cc         tr_seri(i,k,10)=0.
-c           if ((xlat(i)>=40.).and.(xlat(i)<=60.)) then
-c           if ((xlon(i)>=-180.).and.(xlon(i)<=180.)) then
-c           if ((pplay(klon/2,k)>=1.e4).and.(pplay(klon/2,k)<=1.e5)) then
-c               tr_seri(i,k,10)=1.
-c           endif
-c           endif
-c           endif
-c         else
-c            tr_seri(i,k,10)=0.
-c         end do 
-c         end do
-c
-c         do i=1,klon
-c         do k=1,klev+1
-cc         tr_seri(i,k,11)=0.
-c           if ((xlat(i)>=-60.).and.(xlat(i)<=-40.)) then
-c           if ((xlon(i)>=-180.).and.(xlon(i)<=180.)) then
-c           if ((pplay(klon/2,k)>=1.e4).and.(pplay(klon/2,k)<=1.e5)) then
-c               tr_seri(i,k,11)=1.
-c           endif
-c           endif
-c           endif
-c         else
-c            tr_seri(i,k,11)=0.
-c         end do 
-c         end do
-c
-c         do i=1,klon
-c         do k=1,klev+1
-cc         tr_seri(i,k,12)=0.
-c           if ((xlat(i)>=-20.).and.(xlat(i)<=20.)) then
-c           if ((xlon(i)>=-180.).and.(xlon(i)<=180.)) then
-c           if ((pplay(klon/2,k)>=1.e4).and.(pplay(klon/2,k)<=1.e5)) then
-c               tr_seri(i,k,12)=1.
-c           endif
-c           endif
-c           endif
-c         else
-c            tr_seri(i,k,12)=0.
-c         end do 
-c         end do
-c
-c entre 20-30 km
-c
-c         do i=1,klon
-c         do k=1,klev+1
-cc         tr_seri(i,k,13)=0.
-c           if ((xlat(i)>=60.).and.(xlat(i)<=80.)) then
-c           if ((xlon(i)>=-180.).and.(xlon(i)<=180.)) then
-c           if ((pplay(klon/2,k)>=1.e6).and.(pplay(klon/2,k)<=2.e6)) then
-c               tr_seri(i,k,13)=1.
-c           endif
-c           endif
-c           endif
-c         else
-c            tr_seri(i,k,13)=0.
-c         end do 
-c         end do
-c
-c         do i=1,klon
-c         do k=1,klev+1
-cc         tr_seri(i,k,14)=0.
-c           if ((xlat(i)>=-80.).and.(xlat(i)<=-60.)) then
-c           if ((xlon(i)>=-180.).and.(xlon(i)<=180.)) then
-c           if ((pplay(klon/2,k)>=1.e6).and.(pplay(klon/2,k)<=2.e6)) then
-c               tr_seri(i,k,14)=1.
-c           endif
-c           endif
-c           endif
-c         else
-c            tr_seri(i,k,14)=0.
-c         end do 
-c         end do
-c
-c         do i=1,klon
-c         do k=1,klev+1
-cc         tr_seri(i,k,15)=0.
-c           if ((xlat(i)>=-20.).and.(xlat(i)<=20.)) then
-c           if ((xlon(i)>=-180.).and.(xlon(i)<=180.)) then
-c           if ((pplay(klon/2,k)>=1.e6).and.(pplay(klon/2,k)<=2.e6)) then
-c               tr_seri(i,k,15)=1.
-c           endif
-c           endif
-c           endif
-c         else
-c            tr_seri(i,k,15)=0.
-c         end do 
-c         end do
-c
-c entre 55-65 km
-c
-c         do i=1,klon
-c         do k=1,klev+1
-cc         tr_seri(i,k,16)=0.
-c           if ((xlat(i)>=60.).and.(xlat(i)<=80.)) then
-c           if ((xlon(i)>=-180.).and.(xlon(i)<=180.)) then
-c           if ((pplay(klon/2,k)>=1.e4).and.(pplay(klon/2,k)<=5.e4)) then
-c               tr_seri(i,k,16)=1.
-c           endif
-c           endif
-c           endif
-c           endif
-c         else
-c            tr_seri(i,k,16)=0.
-c         end do 
-c         end do
-c
-c         do i=1,klon
-c         do k=1,klev+1
-cc         tr_seri(i,k,17)=0.
-c           if ((xlat(i)>=-80.).and.(xlat(i)<=-60.)) then
-c           if ((xlon(i)>=-180.).and.(xlon(i)<=180.)) then
-c           if ((pplay(klon/2,k)>=1.e4).and.(pplay(klon/2,k)<=5.e4)) then
-c               tr_seri(i,k,17)=1.
-c          endif
-c          endif
-c           endif
-c           endif
-c         else
-c            tr_seri(i,k,17)=0.
-c         end do 
-c         end do
-c
-c         do i=1,klon
-c         do k=1,klev+1
-cc         tr_seri(i,k,18)=0.
-c           if ((xlat(i)>=-20.).and.(xlat(i)<=20.)) then
-c           if ((xlon(i)>=-180.).and.(xlon(i)<=180.)) then
-c           if ((pplay(klon/2,k)>=1.e4).and.(pplay(klon/2,k)<=5.e4)) then
-c               tr_seri(i,k,18)=1.
-c           endif
-c           endif
-c           endif
-c           endif
-c         else
-c            tr_seri(i,k,18)=0.
-c         end do 
-c         end do
-c
-c b) traceur a une bande en latitudeOK
-c
-c a 65km
-c
-c        do i=1,klon
-c         tr_seri(i,:,1)=0.
-c        if ((xlat(i)>=60.).and.(xlat(i)<=80.)) then
-c              tr_seri(i,20,1)=1.
-c        endif
-c      end do  
-c
-c        do i=1,klon
-c         tr_seri(i,:,2)=0.
-c        if ((xlat(i)>=40.).and.(xlat(i)<=60.)) then
-c              tr_seri(i,20,2)=1.
-c        endif
-c      end do  
-c
-c        do i=1,klon
-c         tr_seri(i,:,3)=0.
-c        if ((xlat(i)>=20.).and.(xlat(i)<=40.)) then
-c              tr_seri(i,20,3)=1.
-c        endif
-c      end do  
-c
-c        do i=1,klon
-c         tr_seri(i,:,4)=0.
-c        if ((xlat(i)>=0.).and.(xlat(i)<=20.)) then
-c              tr_seri(i,20,4)=1.
-c        endif
-c      end do  
-c
-c        do i=1,klon
-c         tr_seri(i,:,5)=0.
-c        if ((xlat(i)>=-20.).and.(xlat(i)<=0.)) then
-c              tr_seri(i,20,5)=1.
-c        endif
-c      end do  
-c
-c        do i=1,klon
-c         tr_seri(i,:,6)=0.
-c        if ((xlat(i)>=-40.).and.(xlat(i)<=-20.)) then
-c              tr_seri(i,20,6)=1.
-c        endif
-c      end do  
-c
-c        do i=1,klon
-c         tr_seri(i,:,7)=0.
-c        if ((xlat(i)>=-60.).and.(xlat(i)<=-40.)) then
-c              tr_seri(i,20,7)=1.
-c        endif
-c      end do  
-c
-c        do i=1,klon
-c         tr_seri(i,:,8)=0.
-c        if ((xlat(i)>=-80.).and.(xlat(i)<=-60.)) then
-c              tr_seri(i,20,8)=1.
-c        endif
-c      end do  
-c
-c a 50km
-c
-c        do i=1,klon
-c        tr_seri(i,:,1)=0.
-c        if ((xlat(i)>=40.).and.(xlat(i)<=60.)) then
-c              tr_seri(i,27,1)=1.
-c        endif
-c      end do  
-c
-c        do i=1,klon
-c         tr_seri(i,:,2)=0.
-c        if ((xlat(i)>=60.).and.(xlat(i)<=80.)) then
-c              tr_seri(i,27,2)=1.
-c        endif
-c      end do  
-c
-c        do i=1,klon
-c         tr_seri(i,:,3)=0.
-c        if ((xlat(i)>=20.).and.(xlat(i)<=40.)) then
-c              tr_seri(i,27,3)=1.
-c        endif
-c      end do  
-c
-c        do i=1,klon
-c         tr_seri(i,:4)=0.
-c        if ((xlat(i)>=0.).and.(xlat(i)<=20.)) then
-c              tr_seri(i,27,4)=1.
-c       endif
-c      end do  
-c
-c        do i=1,klon
-c         tr_seri(i,:,5)=0.
-c        if ((xlat(i)>=-20.).and.(xlat(i)<=0.)) then
-c              tr_seri(i,27,5)=1.
-c        endif
-c      end do  
-c
-c        do i=1,klon
-c         tr_seri(i,:,6)=0.
-c        if ((xlat(i)>=-40.).and.(xlat(i)<=-20.)) then
-c              tr_seri(i,27,6)=1.
-c        endif
-c      end do  
-c
-c        do i=1,klon
-c         tr_seri(i,:,7)=0.
-c        if ((xlat(i)>=-60.).and.(xlat(i)<=-40.)) then
-c              tr_seri(i,27,7)=1.
-c        endif
-c      end do  
-c
-c        do i=1,klon
-c         tr_seri(i,:,8)=0.
-c        if ((xlat(i)>=-80.).and.(xlat(i)<=-60.)) then
-c              tr_seri(i,27,8)=1.
-c        endif
-c      end do  
-c
-c c) traceur a plusieurs bandes en latitude OK
-c
-c         do i=1,klon
-c        tr_seri(i,:,2)=0.
-c        if ((xlat(i)>=50.).and.(xlat(i)<=70.)) then
-c             tr_seri(i,10,2)=1.
-c        endif
-c        if ((xlat(i)>=-10.).and.(xlat(i)<=10.)) then
-c              tr_seri(i,10,2)=1.
-c        endif
-c
-c        if ((xlat(i)>=-70.).and.(xlat(i)<=-50.)) then
-c              tr_seri(i,10,2)=1. 
-c        endif         
-c      end do 
-c
-c d) traceur a une bande en altitude OK
-c
-c       do k=1,klev+1
-c         tr_seri(:,k,1)=0.
-c         if ((pplay(klon/2,k)>=1.e5).and.(pplay(klon/2,k)<=1.e6)) then
-c             tr_seri(:,k,1)=1.
-c         endif
-c       end do 
-c
-c dbis) plusieurs traceurs a une bande en altitude OK
-c
-c bande tres basse tropo
-c      do k=1,klev
-c        tr_seri(:,k,1)=0.
-c        if ((pplay(klon/2,k)>=5.e5).and.(pplay(klon/2,k)<=5.e6)) then
-c            tr_seri(:,k,1)=1.
-c        endif
-c      end do 
-c bande dans les nuages et un peu en-dessous    
-c      do k=1,klev
-c         tr_seri(:,k,2)=0.
-c         if ((pplay(klon/2,k)>=5.e4).and.(pplay(klon/2,k)<=5.e5)) then
-c             tr_seri(:,k,2)=1.
-c         endif
-c       end do 
-cune grosse epaisseur: inclue toute la circulation meridienne
-c      do k=1,klev
-c        tr_seri(:,k,1)=0.
-c        if ((pplay(klon/2,k)>=1.e4).and.(pplay(klon/2,k)<=1.e6)) then
-c            tr_seri(:,k,1)=1.
-c        endif
-c      end do
-cune grosse epaisseur: inclue la mesosphere
-c      do k=1,klev
-c         tr_seri(:,k,2)=0.
-c         if ((pplay(klon/2,k)>=2.e2).and.(pplay(klon/2,k)<=1.e4)) then
-c             tr_seri(:,k,2)=1.
-c         endif
-c       end do 
-c
-c      do k=1,klev
-c        tr_seri(:,k,3)=0.
-c        if ((pplay(klon/2,k)>=5.e1).and.(pplay(klon/2,k)<=5.e2)) then
-c            tr_seri(:,k,3)=1.
-c        endif
-c      end do 
-c
-c e) plusieurs couches verticales de traceurs, a plusieurs bandes en latitude???
-c       
-c au sol
-c         do i=1,klon
-c        tr_seri(i,:,1)=0.
-c        if ((xlat(i)>=50.).and.(xlat(i)<=70.)) then
-c             tr_seri(i,5,1)=1.
-c        endif
-c        if ((xlat(i)>=-10.).and.(xlat(i)<=10.)) then
-c              tr_seri(i,5,1)=1.
-c        endif
-c
-c        if ((xlat(i)>=-70.).and.(xlat(i)<=-50.)) then
-c              tr_seri(i,5,1)=1. 
-c        endif         
-c      end do 
-c 
-c         do i=1,klon
-c        tr_seri(i,:,2)=0.
-c        if ((xlat(i)>=50.).and.(xlat(i)<=70.)) then
-c             tr_seri(i,10,2)=1.
-c        endif
-c        if ((xlat(i)>=-10.).and.(xlat(i)<=10.)) then
-c              tr_seri(i,10,2)=1.
-c      endif
-c
-c        if ((xlat(i)>=-70.).and.(xlat(i)<=-50.)) then
-c              tr_seri(i,10,2)=1. 
-c        endif         
-c      end do 
-c
-c         do i=1,klon
-c        tr_seri(i,:,3)=0.
-c        if ((xlat(i)>=50.).and.(xlat(i)<=70.)) then
-c             tr_seri(i,30,3)=1.
-c        endif
-c        if ((xlat(i)>=-10.).and.(xlat(i)<=10.)) then
-c              tr_seri(i,30,3)=1.
-c        endif
-c
-c        if ((xlat(i)>=-70.).and.(xlat(i)<=-50.)) then
-c              tr_seri(i,30,3)=1. 
-c        endif         
-c      end do 
-c        
-c        do i=1,klon
-c        tr_seri(i,:,4)=0.
-c        if ((xlat(i)>=50.).and.(xlat(i)<=70.)) then
-c             tr_seri(i,45,4)=1.
-c        endif
-c        if ((xlat(i)>=-10.).and.(xlat(i)<=10.)) then
-c              tr_seri(i,45,4)=1.
-c        endif
-c
-c        if ((xlat(i)>=-70.).and.(xlat(i)<=-50.)) then
-c              tr_seri(i,45,4)=1. 
-c        endif         
-c      end do 
-c      
-
Index: trunk/LMDZ.VENUS/libf/phyvenus/phytrac_emiss.F
===================================================================
--- trunk/LMDZ.VENUS/libf/phyvenus/phytrac_emiss.F	(revision 1160)
+++ trunk/LMDZ.VENUS/libf/phyvenus/phytrac_emiss.F	(revision 1160)
@@ -0,0 +1,200 @@
+!
+! $Header: /home/cvsroot/LMDZ4/libf/phylmd/phytrac.F,v 1.16 2006/03/24 15:06:23 lmdzadmin Exp $
+!
+c
+c
+      SUBROUTINE phytrac_emiss (timesimu,
+     I                    debutphy,
+     I                    lafin,
+     I                    nqmax,
+     I                    nlon,
+     I                    nlev, 
+     I                    pdtphys,
+     I                    paprs,
+     I                    xlat,xlon,
+     O                    tr_seri)
+
+c======================================================================
+c Auteur(s) FH
+c Objet: Moniteur general des tendances traceurs
+c
+cAA Remarques en vrac:
+cAA--------------------
+cAA 1/ le call phytrac se fait avec nqmax 
+c
+c SL: Janvier 2014
+c Version developed for surface emission
+c Maybe could be used just to compute the 'source' variable from physiq
+c
+c======================================================================
+      USE ioipsl
+      USE infotrac
+      USE control_mod
+      use dimphy
+      USE comgeomphy
+      IMPLICIT none
+#include "YOMCST.h"
+#include "dimensions.h"
+#include "clesphys.h"
+#include "temps.h"
+#include "paramet.h"
+c======================================================================
+
+c Arguments:
+
+c   EN ENTREE:
+c   ==========
+
+      real timesimu   ! duree depuis debut simu (s)
+      logical debutphy       ! le flag de l'initialisation de la physique
+      logical lafin          ! le flag de la fin de la physique
+      integer nqmax ! nombre de traceurs auxquels on applique la physique
+      integer nlon  ! nombre de points horizontaux
+      integer nlev  ! nombre de couches verticales
+      real pdtphys  ! pas d'integration pour la physique (seconde)
+      real paprs(nlon,nlev+1)  ! pression pour chaque inter-couche (en Pa)
+      REAL xlat(nlon)       ! latitudes pour chaque point 
+      REAL xlon(nlon)       ! longitudes pour chaque point 
+
+c   EN ENTREE/SORTIE:
+c   =================
+
+      real tr_seri(nlon,nlev,nqmax) ! traceur  
+
+cAA ----------------------------
+cAA  VARIABLES LOCALES TRACEURS
+cAA ----------------------------
+
+C les traceurs
+
+c===================
+c it--------indice de traceur
+c k,i---------indices long, vert
+c===================
+c pour emission volcan
+      real :: deltatr(klon,klev,nqtot)
+
+      integer,parameter :: nbsrc=2,nblat=5,nblon=4
+      integer,parameter :: Nemiss=1   ! duree emission (Ed)
+      real,save :: source_volcan(nbsrc)
+      real,save :: lat_volcan(nblat),lon_volcan(nblon)
+      real,save :: area_emiss(nblat,nblon)
+      integer,save :: ig_volcan(nblat,nblon)
+
+c======================================================================
+
+      INTEGER i, k, it
+      integer ilat,ilon,iemiss
+      real    deltalat,deltalon
+
+c Variables liees a l'ecriture de la bande histoire physique
+
+c Variables locales pour effectuer les appels en serie
+c----------------------------------------------------
+
+      REAL d_tr(nlon,nlev) ! tendances de traceurs 
+
+      character*20 modname
+      character*80 abort_message
+
+c======================================================================
+
+      modname = 'phytrac_emiss'
+c EMISSION TRACEURS 
+
+c---------
+c debutphy
+c---------
+      if (debutphy) then
+         print*,"DEBUT PHYTRAC"
+         print*,"PHYTRAC: EMISSION"
+
+c=============================================================
+c=============================================================
+c=============================================================
+c   Initialisation des traceurs
+c=============================================================
+c=============================================================
+c=============================================================
+
+C=========================================================================
+C=========================================================================
+c Caracteristiques des traceurs emis:
+
+c nombre total de traceur
+         if (nbsrc*nblat*nblon .gt. nqtot) then
+            write(*,*) "Attention, pas assez de traceurs"
+            write(*,*) "le dernier sera bien le dernier"
+         endif
+
+c source en kg/s
+         source_volcan(1) = 1.
+         source_volcan(2) = 1000.
+c localisation volcan
+         lat_volcan(1) =  70.
+         lat_volcan(2) =  35.
+         lat_volcan(3) =   0.
+         lat_volcan(4) = -35.
+         lat_volcan(5) = -70.
+         lon_volcan(1) = -120.
+         lon_volcan(2) =  -30.
+         lon_volcan(3) =   60.
+         lon_volcan(4) =  150.
+
+         deltalat = 180./jjm
+         deltalon = 360./jjm
+         do i=1,nlon
+          do ilat=1,nblat
+           do ilon=1,nblon
+            if ((xlat(i).ge.lat_volcan(ilat))
+     &     .and.((xlat(i)-deltalat).lt.lat_volcan(ilat))
+     &     .and.(xlon(i).le.lon_volcan(ilon))
+     &     .and.((xlon(i)+deltalon).gt.lon_volcan(ilon)) ) then
+             ig_volcan(ilat,ilon)= i
+             area_emiss(ilat,ilon) = airephy(i)
+            endif
+           enddo
+          enddo
+         enddo
+
+C=========================================================================
+C=========================================================================
+
+c-------------
+c fin debutphy
+c-------------
+      ENDIF  ! fin debutphy 
+
+c======================================================================
+c Emission d'un traceur pendant un certain temps
+c======================================================================
+         do i = 1,nlon
+          do iemiss = 1,nbsrc
+           do ilat  = 1,nblat
+            do ilon  = 1,nblon
+             it=min(iemiss*ilat*ilon,nqtot)
+
+             if (i .eq. ig_volcan(ilat,ilon)) then
+
+c source appliquee pendant Nemiss Ed
+               if (timesimu .lt. 86400.*Nemiss) then
+
+c source en kg/kg/s
+           deltatr(i,1,it) = source_volcan(iemiss)*RG
+     $     /(area_emiss(ilat,ilon)*(paprs(i,1) - paprs(i,2)))
+           tr_seri(i,1,it) = tr_seri(i,1,it) + deltatr(i,1,it)*pdtphys
+
+               else
+           deltatr(i,1,it) = 0.
+               end if  ! duree emission
+
+             end if ! i localisation
+            end do
+           end do
+          end do
+         end do
+c======================================================================
+c======================================================================
+
+      RETURN
+      END
Index: trunk/LMDZ.VENUS/libf/phyvenus/phytrac_relax.F
===================================================================
--- trunk/LMDZ.VENUS/libf/phyvenus/phytrac_relax.F	(revision 1160)
+++ trunk/LMDZ.VENUS/libf/phyvenus/phytrac_relax.F	(revision 1160)
@@ -0,0 +1,272 @@
+!
+! $Header: /home/cvsroot/LMDZ4/libf/phylmd/phytrac.F,v 1.16 2006/03/24 15:06:23 lmdzadmin Exp $
+!
+c
+c
+      SUBROUTINE phytrac_relax (debutphy,lafin,
+     I                    nqmax,
+     I                    nlon,
+     I                    nlev, 
+     I                    pdtphys,
+     I                    pplay,
+     O                    tr_seri)
+
+c======================================================================
+c Auteur(s) FH
+c Objet: Moniteur general des tendances traceurs
+c
+cAA Remarques en vrac:
+cAA--------------------
+cAA 1/ le call phytrac se fait avec nqmax 
+c
+c SL: Janvier 2014
+c Version developed by E. Marcq for pseudo-chemistry relaxation
+c See Marcq&Lebonnois 2013.
+c
+c======================================================================
+      USE ioipsl
+      USE infotrac
+      USE control_mod
+      use dimphy
+      USE comgeomphy
+      IMPLICIT none
+#include "YOMCST.h"
+#include "dimensions.h"
+#include "clesphys.h"
+#include "temps.h"
+#include "paramet.h"
+c======================================================================
+
+c Arguments:
+
+c   EN ENTREE:
+c   ==========
+
+      logical debutphy       ! le flag de l'initialisation de la physique
+      logical lafin          ! le flag de la fin de la physique
+      integer nqmax ! nombre de traceurs auxquels on applique la physique
+      integer nlon  ! nombre de points horizontaux
+      integer nlev  ! nombre de couches verticales
+      real pdtphys  ! pas d'integration pour la physique (seconde)
+      real pplay(nlon,nlev)  ! pression pour le mileu de chaque couche (en Pa)
+
+c   EN ENTREE/SORTIE:
+c   =================
+
+      real tr_seri(nlon,nlev,nqmax) ! traceur  
+
+cAA ----------------------------
+cAA  VARIABLES LOCALES TRACEURS
+cAA ----------------------------
+
+C les traceurs
+
+c===================
+c it--------indice de traceur
+c k,i---------indices long, vert
+c===================
+c Variables deja declarees dont on a besoin pour traceurs   
+c   k,i,it,tr_seri(nlon,nlev,nqmax),pplay(nlon,nlev), 
+      integer nqCO_OCS
+c      real pzero,gamma
+c      parameter (pzero=85000.)
+c      parameter (gamma=5000.)
+      REAL alpha
+      real deltatr(nlon,nlev,nqtot) ! ecart au profil de ref zprof
+      real,save,allocatable :: zprof(:,:)
+      real,save,allocatable :: tau(:,:) ! temps de relaxation vers le profil (s)
+c======================================================================
+
+      INTEGER i, k, it
+
+c Variables liees a l'ecriture de la bande histoire physique
+
+c Variables locales pour effectuer les appels en serie
+c----------------------------------------------------
+
+      REAL d_tr(nlon,nlev) ! tendances de traceurs 
+
+      character*20 modname
+      character*80 abort_message
+
+c======================================================================
+
+      modname = 'phytrac_relax'
+c TRACEURS TYPE CO ET OCS
+      nqCO_OCS   = 6
+
+c---------
+c debutphy
+c---------
+      if (debutphy) then
+         print*,"DEBUT PHYTRAC"
+         print*,"PHYTRAC: RELAXATION"
+         allocate(zprof(nlev,nqtot),tau(nlev,nqtot))
+
+c=============================================================
+c=============================================================
+c=============================================================
+c   Initialisation des traceurs
+c=============================================================
+c=============================================================
+c=============================================================
+
+C=========================================================================
+C=========================================================================
+
+c II) Declaration d'un profil vertical de traceur OK
+c 
+c zprof   = profil de rappel 
+c
+c 1 -> CO ; 2 -> OCS 
+c def des profils en log(a) = a * log(P) + b par morceaux, cf. pollack et al
+c tr_seri en ppm
+c (initialisation seulement si ceux-ci sont nuls)
+
+c ICI, ON UTILISE 3 CONSTANTES DE TEMPS DIFFERENTES POUR CHAQUE,
+c DONC TRACEURS 1 A 3 POUR CO ET 4 A 6 POUR OCS
+C=========================================================================
+
+
+c Constantes de rappel:
+
+       print*,"INIT TAU"
+       do k=1,nlev
+         tau(k,1)=1.e6
+         tau(k,2)=1.e7
+         tau(k,3)=1.e8
+         tau(k,4)=1.e6
+         tau(k,5)=1.e7
+         tau(k,6)=1.e8
+       enddo
+
+c CO
+
+      do it=1,3 
+       print*,"INIT ZPROF ",tname(it)
+       do k=1,nlev
+         zprof(k,it)=0.
+c pour l'instant, tau fixe, mais possibilite de le faire varier avec z
+        if (pplay(nlon/2,k) >= 4.8e6) then
+           zprof(k,it)=14.
+        endif
+        if ((pplay(nlon/2,k)<=4.8e6).and.(pplay(nlon/2,k)>=1.9e6)) then
+           alpha=(log(pplay(nlon/2,k))-log(1.9e6))/
+     .     (log(4.8e6)-log(1.9e6))
+           zprof(k,it)=20.*(14./20.)**alpha
+        endif
+        if ((pplay(nlon/2,k)<=1.9e6).and.(pplay(nlon/2,k)>=1.5e5)) then
+           alpha=(log(pplay(nlon/2,k))-log(1.5e5))/
+     .     (log(1.9e6)-log(1.5e5))
+           zprof(k,it)=39.*(20./39.)**alpha
+        endif
+        if ((pplay(nlon/2,k)<=1.5e5).and.(pplay(nlon/2,k)>=1.1e4)) then
+           alpha=(log(pplay(nlon/2,k))-log(1.1e4))/
+     .     (log(2.73e5)-log(1.1e4))
+           zprof(k,it)=50.*(39./50.)**alpha
+        endif
+        if ((pplay(nlon/2,k)<=1.1e4).and.(pplay(nlon/2,k)>=1.3e3)) then
+           alpha=(log(pplay(nlon/2,k))-log(1.3e3))/
+     .     (log(1.1e4)-log(1.3e3))
+           zprof(k,it)=2.*(50./2.)**alpha
+        endif
+        if ((pplay(nlon/2,k)<=1.3e3).and.(pplay(nlon/2,k)>=2.4)) then
+           alpha=(log(pplay(nlon/2,k))-log(2.4))/
+     .     (log(1.3e3)-log(2.4))
+           zprof(k,it)=1000.*(2./1000.)**alpha
+        endif
+        if (pplay(nlon/2,k) <= 2.4) then
+           zprof(k,it)=1000.
+        endif
+       enddo 
+       print*,zprof(:,it)
+ 
+c OCS
+       print*,"INIT ZPROF ",tname(it+3)
+       do k=1,nlev
+         zprof(k,it+3)=0.
+         if (pplay(nlon/2,k) >= 4.8e6) then
+           zprof(k,it+3)=30.
+         endif
+         if ((pplay(nlon/2,k)<=4.8e6).and.(pplay(nlon/2,k)>=9.4e5))
+     *   then
+           alpha=(log(pplay(nlon/2,k))-log(9.4e5))/
+     *     (log(4.8e6)-log(9.4e5))
+           zprof(k,it+3)=20.*(30/20.)**alpha
+         endif
+         if ((pplay(nlon/2,k)<=9.4e5).and.(pplay(nlon/2,k)>=4.724e5))
+     *   then
+           alpha=(log(pplay(nlon/2,k))-log(4.724e5))/
+     *     (log(9.4e5)-log(4.724e5))
+           zprof(k,it+3)=0.5*(20/0.5)**alpha
+         endif
+         if ((pplay(nlon/2,k)<=4.724e5).and.(pplay(nlon/2,k)>=1.1e4))
+     *   then
+           alpha=(log(pplay(nlon/2,k))-log(1.1e4))/
+     *     (log(4.724e5)-log(1.1e4))
+           zprof(k,it+3)=0.005*(0.5/0.005)**alpha
+         endif
+         if (pplay(nlon/2,k)<=1.1e4) then
+           zprof(k,it+3)=0.
+         endif
+       end do
+       print*,zprof(:,it+3)
+      enddo 
+
+c Initialisation du traceur s'il est nul:
+       do it=1,nqCO_OCS
+        if ((tr_seri(nlon/2,1,it).eq.0.).and.
+     .      (tr_seri(nlon/2,nlev/2,it).eq.0.).and.
+     .      (tr_seri(nlon/2,nlev,it).eq.0.)) then
+         print*,"INITIALISATION DE ",tname(it)
+         do k=1,nlev
+           do i=1,nlon
+             tr_seri(i,k,it) = zprof(k,it)
+           enddo
+         enddo
+        endif
+       enddo
+
+C=========================================================================
+C=========================================================================
+
+c-------------
+c fin debutphy
+c-------------
+      ENDIF  ! fin debutphy 
+
+c======================================================================
+c Rappel vers un profil
+c======================================================================
+         do it=1,nqCO_OCS
+           do k=1,nlev
+             do i=1,nlon
+c VERIF
+           if (tr_seri(i,k,it).lt.0) then
+             print*,"Traceur negatif AVANT rappel:",i,k,it
+             stop
+           endif
+c FIN VERIF
+
+           deltatr(i,k,it) = (-tr_seri(i,k,it)+zprof(k,it))/tau(k,it)
+           tr_seri(i,k,it) =  tr_seri(i,k,it) + deltatr(i,k,it)*pdtphys
+
+c VERIF
+           if (tr_seri(i,k,it).lt.0) then
+             print*,"APRES rappel:",i,k,it,
+     .  deltatr(i,k,it),zprof(k,it),tr_seri(i,k,it),pdtphys/tau(k,it)
+             stop
+           endif
+c FIN VERIF
+             enddo
+           enddo
+         enddo
+
+c======================================================================
+c======================================================================
+
+
+      RETURN
+      END
+      
+      
