!
! $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                    coefh,
     I                    yu1,
     I                    yv1,
     I                    ftsol,
     I                    xlat,
     I                    xlon,
     I                    zlev, 
     I                    presnivs,
     I                    pphis,
     I                    pphi,
     I                    albsol,
     O                    tr_seri)

      USE ioipsl

      IMPLICIT none
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======================================================================
#include "YOMCST.h"
#include "dimensions.h"
#include "dimphy.h"
#include "clesphys.h" !///utile?
#include "temps.h"
#include "paramet.h"
#include "control.h"
#include "comgeomphy.h"
#include "advtrac.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 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 Rem : nqmax : nombre de vrais traceurs est defini dans dimphy.h
c
c   convection:
c   -----------
c
      REAL pmfu(nlon,nlev)  ! flux de masse dans le panache montant
      REAL pmfd(nlon,nlev)  ! flux de masse dans le panache descendant
      REAL pen_u(nlon,nlev) ! flux entraine dans le panache montant

c
c   Couche limite:
c   --------------
c
      REAL coefh(nlon,nlev) ! coeff melange CL
      REAL yu1(nlon)        ! vents au premier niveau
      REAL yv1(nlon)        ! vents au premier niveau
      REAL xlat(nlon)       ! latitudes pour chaque point 
      REAL xlon(nlon)       ! longitudes pour chaque point 
      REAL zlev(nlon,nlev+1) ! altitude a chaque niveau (interface inferieure de la couche) 
cAA
cAA Arguments necessaires pour les sources et puits de traceur:
cAA ----------------
cAA
      real ftsol(nlon)  ! Temperature du sol (surf)(Kelvin)

cAA ----------------------------
cAA  VARIABLES LOCALES TRACEURS
cAA ----------------------------
cAA

      CHARACTER*2 itn
C maf ioipsl
      CHARACTER*2 str2
      INTEGER nhori, nvert
      REAL zsto, zout, zjulian
      INTEGER nid_tra
      SAVE nid_tra
      INTEGER nid_tra2,nid_tra3
      SAVE nid_tra2,nid_tra3
      INTEGER ndex(1)
      INTEGER ndex2d(iim*(jjm+1)),ndex3d(iim*(jjm+1)*klev)
      REAL zx_tmp_2d(iim,jjm+1), zx_tmp_3d(iim,jjm+1,klev)
      REAL zx_lon(iim,jjm+1), zx_lat(iim,jjm+1)
c
      integer itau_w   ! pas de temps ecriture = nstep + itau_phy
c

C
C Variables liees a l'ecriture de la bande histoire : phytrac.nc
c
c      INTEGER ecrit_tra
c      SAVE ecrit_tra   
      logical ok_sync
      parameter (ok_sync = .true.)
C
C les traceurs
C
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), 
      real zprof(klon,klev,nqmx)
c      real pzero,gamma
c      parameter (pzero=85000.)
c      parameter (gamma=5000.)
      REAL alpha
      real deltatr(klon,klev,nqmx) ! ecart au profil de ref zprof
      real tau(klev,nqmx)          ! temps de relaxation vers le profil (s)
      save zprof,tau
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*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./

c======================================================================
         ps(:)=paprs(:,1)
c---------
c debutphy
c---------
         if (debutphy) then
                 print*,"DEBUT PHYTRAC"
C         
c=============================================================
c   Initialisation des traceurs
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      

c II) Declaration d'un profil vertical de traceur OK
c 
c tr_seri = profil initial
c zprof   = profil de rappel 
c           (identiques a l'initialisation)
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 Constantes de rappel:

       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 
       do k=1,klev
         tr_seri(:,k,it)=0.
c pour l'instant, tau fixe, mais possibilite de le faire varier avec z
        if (pplay(klon/2,k) >= 1.9e6) then
           tr_seri(:,k,it)=20.
        endif
        if ((pplay(klon/2,k)<=1.9e6).and.(pplay(klon/2,k)>=2.73e5)) then
           alpha=(log(pplay(klon/2,k))-log(2.73e5))/
     .     (log(1.9e6)-log(2.73e5))
           tr_seri(:,k,it)=30.*(20./30.)**alpha
        endif
        if ((pplay(klon/2,k)<=2.73e5).and.(pplay(klon/2,k)>=1.1e4)) then
           alpha=(log(pplay(klon/2,k))-log(1.1e4))/
     .     (log(2.73e5)-log(1.1e4))
           tr_seri(:,k,it)=50.*(30./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))
           tr_seri(:,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))
           tr_seri(:,k,it)=1000.*(2./1000.)**alpha
        endif
        if (pplay(klon/2,k) <= 2.4) then
           tr_seri(:,k,it)=1000.
        endif
       enddo 
 
c OCS
       do k=1,klev
         tr_seri(:,k,it+3)=0.
         if (pplay(klon/2,k) >= 9.45e5) then
           tr_seri(:,k,it+3)=16.
         endif
         if ((pplay(klon/2,k)<=9.45e5).and.(pplay(klon/2,k)>=4.724e5))
     *   then
           alpha=(log(pplay(klon/2,k))-log(4.724e5))/
     *     (log(9.45e5)-log(4.724e5))
           tr_seri(:,k,it+3)=0.5*(16/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))
           tr_seri(:,k,it+3)=0.005*(0.5/0.005)**alpha
         endif
         if (pplay(klon/2,k)<=1.1e4) then
           tr_seri(:,k,it+3)=0.
         endif
       end do
      enddo 

       do it=1,nqmax
         do k=1,klev
           do i=1,klon
             zprof(i,k,it) = tr_seri(i,k,it)
           enddo
         enddo
       enddo

c-------------
c fin debutphy
c-------------
      ENDIF  ! fin debutphy 

c======================================================================
c Rappel vers un profil
c======================================================================
         do it=1,nqmax
           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(i,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(i,k,it),tr_seri(i,k,it),pdtphys/tau(k,it)
             stop
           endif
c FIN VERIF
             enddo
           enddo
         enddo

c======================================================================
c   Calcul de l'effet de la couche limite remis directement dans physiq
c======================================================================


      RETURN
      END