SUBROUTINE phytrac (firstcall,lastcall,
     .                   nqmax,nmicro,ptimestep,appkim,dtkim,
     .                   pplev,pplay,delp,ptemp,pmu0,pfract,pdecli,
     .                   lonsol,tr_seri,qaer,d_tr_mph,d_tr_kim)

c======================================================================
c    S. Lebonnois, mai 2008
c
c  Arguments:
c
c firstcall----input-L-variable logique indiquant le premier passage
c lastcall-----input-L-variable logique indiquant le dernier passage
c nqmax--------input-I-nombre de traceurs (total)
c nmicro-------input-I-nombre de traceurs microphysiques !! doivent etre toujours en premiers!!
c ptimestep----input-R-pas d'integration pour la physique (seconde)
c appkim-------input-I-appel a la chimie
c dtkim--------input-R-pas de temps chimique (seconde)
c pplev--------input-R-pression pour chaque inter-couche (en Pa)
c pplay--------input-R-pression pour chaque couche (en Pa)
c delp---------input-R-epaisseur d'une couche (en Pa)
c ptemp--------input-R-temperature (K)
c pmu0---------input-R-cos angle zenithal
c pfract-------input-R-fractional day
c pdecli-------input-R-declinaison en radian
c lonsol-------input-R-longitude solaire en radian
c tr_seri------input-R-mass mixing ratio traceurs (kg/kg) 
c d_tr_mph----output-R-tendance microphysique de "qx" (kg/kg/s)
c d_tr_kim----output-R-tendance chimique de "qx" (kg/kg/s)
c======================================================================
      USE infotrac
      use dimphy
      IMPLICIT none
#include "dimensions.h"
#include "clesphys.h"
#include "YOMCST.h"

c======================================================================
c Variables argument:
c
      LOGICAL firstcall,lastcall
      INTEGER nqmax,nmicro,nlat,appkim
      REAL ptimestep,dtkim
      REAL pplev(klon,klev+1),pplay(klon,klev+1),delp(klon,klev)
      REAL ptemp(klon,klev)
      REAL pmu0(klon), pfract(klon), pdecli, lonsol
      REAL tr_seri(klon,klev,nqmax)
      REAL qaer(klon,klev,nqmax)
      REAL d_tr_mph(klon,klev,nqmax),d_tr_kim(klon,klev,nqmax)

c======================================================================
c Local variables

c grandeurs en moyennes zonales
      REAL zplev(jjm+1,klev+1),zplay(jjm+1,klev)
      REAL ztemp(jjm+1,klev),zmu0(jjm+1),zfract(jjm+1)
      real temp_eq(klev),press_eq(klev)
      REAL zqaer(jjm+1,klev,nqmax)    ! et non nmicro... Permet nmicro=0.
      REAL zqaer0(jjm+1,klev,nqmax) 
      REAL pdqmfi(jjm+1,klev,nqmax)
      REAL ychim(jjm+1,klev,nqmax-nmicro) 
c La saturation n est calculee qu une seule fois: sauvegarde qysat
c La chimie n est pas calculee tous les pas, il faut donc
c                      sauvegarder les sorties de la chimie
      REAL,save,allocatable :: qysat(:,:),pdyfi(:,:,:) 
      
      character*10 nomqy(nqmax-nmicro+1)
      integer      i,j,l,iq,ig0
      
c======================================================================
c======================================================================

      if (firstcall) then
       allocate(qysat(klev,nqmax-nmicro),pdyfi(jjm+1,klev,nqmax-nmicro))
      endif

c-----------------------------------------------------------------------
c   convertion moyennes zonales et changement d unites pour microphy
c   ---------------------------------

c     print*,'CONVERSION 2D ET CHANGEMENT UNITES (PHYTRAC)'

      zplev = 0.0
      zplay = 0.0
      ztemp = 0.0
      zqaer = 0.0
      ychim = 0.0
      zmu0  = 0.0
      zfract= 0.0
      
       do l=1,llm+1
         zplev(1,l) = pplev(1,l)
	 do j=2,jjm
	    ig0=1+(j-2)*iim
	    do i=1,iim
	       zplev(j,l) = zplev(j,l) + pplev(ig0+i,l)/iim
	    enddo
	 enddo
         zplev(jjm+1,l) = pplev(klon,l)
       enddo

       do l=1,llm
         ztemp(1,l) = ptemp(1,l)
         zplay(1,l) = pplay(1,l)
	 do j=2,jjm
	    ig0=1+(j-2)*iim
	    do i=1,iim
	       ztemp(j,l) = ztemp(j,l) + ptemp(ig0+i,l)/iim
	       zplay(j,l) = zplay(j,l) + pplay(ig0+i,l)/iim
	    enddo
	 enddo
         ztemp(jjm+1,l) = ptemp(klon,l)
         zplay(jjm+1,l) = pplay(klon,l)
	 temp_eq  = ztemp((jjm+1)/2,:)
	 press_eq = zplay((jjm+1)/2,:)/100. ! en mbar
       enddo

         zmu0(1)   = pmu0(1)
         zfract(1) = pfract(1)
	 do j=2,jjm
	    ig0=1+(j-2)*iim
	    do i=1,iim
	       zmu0(j)   = zmu0(j)   + pmu0(ig0+i)/iim
	       zfract(j) = zfract(j) + pfract(ig0+i)/iim
	    enddo
	 enddo
         zmu0(jjm+1)   = pmu0(klon)
         zfract(jjm+1) = pfract(klon)

c TRACEURS MICROPHYSIQUES

      if (microfi.eq.1) then

      do iq=1,nmicro
c        print*,tname(iq)

c       Traceurs microphysiques: passage en extensif: n/kg --> n/m^2
         DO l=1,llm
          DO i = 1, klon
            qaer(i,l,iq) = tr_seri(i,l,iq)*delp(i,l)/RG
          ENDDO
         ENDDO

       do l=1,llm
         zqaer(1,l,iq) = qaer(1,l,iq)
	 do j=2,jjm
	    ig0=1+(j-2)*iim
	    do i=1,iim
	       zqaer(j,l,iq) = zqaer(j,l,iq) + qaer(ig0+i,l,iq)/iim
	    enddo
	 enddo
         zqaer(jjm+1,l,iq) = qaer(klon,l,iq)
       enddo
      enddo

      endif   ! microfi

c AUTRES TRACEURS
      
      if (nqmax.gt.nmicro) then
      do iq=nmicro+1,nqmax
       do l=1,llm
         ychim(1,l,iq-nmicro) = tr_seri(1,l,iq)
	 do j=2,jjm
	    ig0=1+(j-2)*iim
	    do i=1,iim
	       ychim(j,l,iq-nmicro) = ychim(j,l,iq-nmicro) 
     .                              + tr_seri(ig0+i,l,iq)/iim
	    enddo
	 enddo
         ychim(jjm+1,l,iq-nmicro) = tr_seri(klon,l,iq)
       enddo
         nomqy(iq-nmicro) = tname(iq)

c       print*,iq-nmicro,nomqy(iq-nmicro)

      enddo
      nomqy(nqmax-nmicro+1) = "HV"
      endif

c-----------------------------------------------------------------------
c   initialisation des qysat au premier appel:
c   ---------------------------------

c!! ATTENTION, qysat pris uniquement a l'equateur
c!!  justifie puisque dans cette region, les var de t et p sont faibles...

      if (firstcall .and. chimi .and.(nqmax.gt.nmicro)) then
           call inicondens(nqmax-nmicro,press_eq,temp_eq,nomqy,qysat)
      endif

c-----------------------------------------------------------------------
c     Appel de la microphysique
c    --------------------------

       pdqmfi = 0.
       zqaer0 = zqaer
       
       IF(firstcall) THEN 
        print*,'MICROPHYSIQUE ',MICROFI
       ENDIF

       IF(MICROFI.eq.1) THEN

         IF(firstcall) THEN 
          print*,'OH! On passe dans la microphysique'
         ENDIF

         CALL pg2(zplev,ztemp,zqaer,pdqmfi    ! tendances 2D, /s
     &   ,nmicro,ptimestep,zmu0,zfract,lastcall)

       ELSE

        IF(firstcall) THEN 
          print*,'MICROPHYSIQUE OFF-LINE',MICROFI
          if (nmicro.gt.0) then
            CALL pg2(zplev,ztemp,zqaer,pdqmfi
     &    ,nmicro,ptimestep,zmu0,zfract,lastcall)
          endif
        ENDIF

       ENDIF

       zqaer  = zqaer+pdqmfi*ptimestep
       pdqmfi = (zqaer-zqaer0)/ptimestep
       
c-----------------------------------------------------------------------
c     Condensation
c    -------------

      if ((chimi).and.(nqmax.gt.nmicro)) then

c   tendance (en /s) passee sur pdqmfi(nmicro+1 a nqmax)
c        print*,'Condensation'

        do iq=1,nqmax-nmicro
           do l=1,llm
              do j=1,jjm+1
                 if (ychim(j,l,iq).gt.qysat(l,iq)) then
           pdqmfi(j,l,nmicro+iq)= (-ychim(j,l,iq)+qysat(l,iq)) !delta y
     .                            / ptimestep                  ! / dt
                 endif
              enddo
           enddo
        enddo

      ENDIF

c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
c eventuellement, modif initiale de la compo
c
c   tendance (en /s) passee sur pdqmfi(nmicro+1 a nqmax)
c
c     if (firstcall .and. chimi .and.(nqmax.gt.nmicro)) then
c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
c!!!remise de CH4 a 1.5%!!!!!!!!!!!!!!!!!!!!!!
c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
c       do iq=1,nqmax-nmicro
c         if (nomqy(iq).eq."CH4") then
c          do l=1,llm
c             do j=1,jjm+1
c                if (ychim(j,l,iq).le.0.015) then
c          pdqmfi(j,l,nmicro+iq)= (-ychim(j,l,iq)+0.015) !delta y
c    .                            / ptimestep                  ! / dt
c                endif
c             enddo
c          enddo
c         endif
c       enddo
c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
c         
c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
c!!!remise de C2H2 a 1.e-5 max !!!!!!!!!!!!!
c!!!remise de C2H6 a 3.e-5 max !!!!!!!!!!!!!
c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
c       do iq=1,nqmax-nmicro
c         if (nomqy(iq).eq."C2H2") then
c          do l=1,llm
c             do j=1,jjm+1
c                if (ychim(j,l,iq).gt.1.e-5) then
c          pdqmfi(j,l,nmicro+iq)= (-ychim(j,l,iq)+1.e-5) !delta y
c    .                            / ptimestep                  ! / dt
c                endif
c             enddo
c          enddo
c         endif
c         if (nomqy(iq).eq."C2H6") then
c          do l=1,llm
c             do j=1,jjm+1
c                if (ychim(j,l,iq).gt.3.e-5) then
c          pdqmfi(j,l,nmicro+iq)= (-ychim(j,l,iq)+3.e-5) !delta y
c    .                            / ptimestep                  ! / dt
c                endif
c             enddo
c          enddo
c         endif
c       enddo
c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
c     endif
c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
          
c-----------------------------------------------------------------------
c     Appel de la chimie
c    --------------------------

      if((appkim.eq.1).and.(chimi)) then
        print*,'On passe dans la CHIMIE'

c       do iq=1,nqmax-nmicro
c         if (nomqy(iq).eq."C2H2") then
c           print*,"C2H2top=",ychim(:,klev,iq)
c         endif
c       enddo

c Appel Chimie
c ------------
       CALL calchim(nqmax-nmicro,ychim,nomqy,pdecli,lonsol,dtkim,
     .              ztemp,zplay,zplev,
     .              pdyfi)   
c ychim ne doit pas etre modifie, pdyfi en /s
	 
      endif
      
c-----------------------------------------------------------------------
c   retour des tendances vers 3D
c   ---------------------------------

c TRACEURS MICROPHYSIQUES

      if (microfi.eq.1) then

      DO iq=1,nmicro
         DO l=1,llm
            d_tr_mph(1,l,iq) = pdqmfi(1,l,iq)
            ig0          = 2
            DO j=2,jjm
               DO i = 1, iim
                  d_tr_mph(ig0,l,iq)  = pdqmfi(j,l,iq)
     &                 *qaer(ig0,l,iq)/zqaer0(j,l,iq)
                  ig0             = ig0 + 1
               ENDDO
            ENDDO
            d_tr_mph(ig0,l,iq) = pdqmfi(jjm+1,l,iq)
         ENDDO
      ENDDO

      do iq=1,nmicro
         DO l=1,llm
          DO i = 1, klon
c  Traceurs microphysiques: passage en intensif: n/m^2 --> n/kg
            d_tr_mph(i,l,iq) = d_tr_mph(i,l,iq)*RG/delp(i,l)
          ENDDO
         ENDDO
      enddo

      endif   ! microfi

c AUTRES TRACEURS

      if ((chimi).and.(nqmax.gt.nmicro)) then
c on passe de pdyfi (tendance chimique en /s calculee quand chimie appelee)
c          a  d_tr_kim (tendance chimique 3D en /s, passee a physiq)
c  et de pdqmfi a d_tr_mph (tendance condensation 3D en /s passee a physiq)

      DO iq=nmicro+1,nqmax
         DO l=1,llm
            d_tr_kim(1,l,iq) = pdyfi(1,l,iq-nmicro)
            d_tr_mph(1,l,iq) = pdqmfi(1,l,iq)
            ig0          = 2
            DO j=2,jjm
               DO i = 1, iim
                  d_tr_kim(ig0,l,iq)  = pdyfi(j,l,iq-nmicro)
     &                 *tr_seri(ig0,l,iq)/ychim(j,l,iq-nmicro)
                  d_tr_mph(ig0,l,iq)  = pdqmfi(j,l,iq)
     &                 *tr_seri(ig0,l,iq)/ychim(j,l,iq-nmicro)
                  ig0             = ig0 + 1
               ENDDO
            ENDDO
            d_tr_kim(ig0,l,iq) = pdyfi(jjm+1,l,iq-nmicro)
            d_tr_mph(ig0,l,iq) = pdqmfi(jjm+1,l,iq)
         ENDDO
      ENDDO

      endif   ! chimi

      RETURN
      END