SUBROUTINE phytrac (firstcall,lastcall,
     .                   nqmax,nmicro,ptimestep,appkim,dtkim,
     .                   pplev,pplay,delp,ptemp,pmu0,pfract,pdecli,
     .                   lonsol,
     .                   pu,pv,pzlev,pzlay,ftsol,
     .                   tr_seri,qaer,d_tr_mph,d_tr_kim,
     .                   fclat,reservoir)

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 pu-----------input-R-vitesse dans la direction X (de O a E) en m/s (1ere couche)
c pv-----------input-R-vitesse Y (de S a N) en m/s                   (1ere couche)
c pzlev--------input-R-altitude pour chaque inter-couche (en m)
c pzlay--------input-R-altitude pour chaque couche (en m)
c ftsol--------input-R-temperature au sol (en K)
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 fclat--------output-R-flux de chaleur latente d'evaporation du reservoir CH4 (J/m2/s)
c reservoir----outpur-R-un reservoir de surface !!! (m)
c======================================================================
      USE infotrac
      use dimphy
      IMPLICIT none
#include "dimensions.h"
#include "clesphys.h"
#include "YOMCST.h"
#include "microtab.h"
#include "varmuphy.h"
#include "diagmuphy.h"
#include "itemps.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 pu(klon),pv(klon)
      REAL pzlev(klon,klev+1),pzlay(klon,klev)
      REAL ftsol(klon)
      REAL tr_seri(klon,klev,nqmax)
      REAL qaer(klon,klev,nqmax)
      REAL d_tr_mph(klon,klev,nqmax),d_tr_kim(klon,klev,nqmax)
      REAL fclat(klon)
      REAL reservoir(klon)

c======================================================================
c Local variables
      REAL qaer0(klon,klev,nqmax)
      REAL prec(klon,5)

c  ASTUCE POUR EVITER klon... EN ATTENDANT MIEUX
      INTEGER   ngrid,NLAYER
      PARAMETER (ngrid=(jjm-1)*iim+2)  ! = klon
      PARAMETER (NLAYER=llm)           ! = klev
* common relatifs au nuages
      real rmcbar(ngrid,NLAYER),xfbar(ngrid,NLAYER,4) 
      integer ncount(ngrid,NLAYER)
      COMMON/rnuabar/ncount,rmcbar,xfbar

      REAL rcloud(klon,klev,nrad),xfrac(klon,klev,4)

      REAL vcl,nuc,xgsn,xmsn,xesn,xasn


      ReAL gaz1(klon,klev),gaz2(klon,klev),gaz3(klon,klev)

      REAL socccld

c grandeurs en moyennes zonales
      REAL zplev(jjm+1,klev+1),zplay(jjm+1,klev),ztsol(jjm+1)
      REAL zzlev(jjm+1,klev+1),zzlay(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 zdqmufi(jjm+1,klev,nqmax)
      REAL ychim(jjm+1,klev,nqmax-nmicro) 
      REAL zgaz1(jjm+1,klev),zgaz2(jjm+1,klev),zgaz3(jjm+1,klev)
      REAL zgaz10(jjm+1,klev),zgaz20(jjm+1,klev),zgaz30(jjm+1,klev)
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,k,l,iq,ig0
      
      REAL zprec(jjm+1,5),zsolesp(jjm+1,klev,3),
     &     zflxesp_i(jjm+1,klev,3)
      REAL ztau_drop(jjm+1,klev),ztau_aer(jjm+1,klev,nrad)
c
c    indice des esp chimiques utilisees dans la microfi  
      integer icldch4,icldc2h6,icldc2h2
      save icldch4,icldc2h6,icldc2h2
     
      real fte,ftm,Lvch4

      REAL tmp,ex,kmin,kmax,dqsq
      REAL dqch4
c      REAL ch4(jjm+1),ch4b(jjm+1),dch4(jjm+1),ch4c(jjm+1,llm)
c       integer ich4
c       common/ch4ind/ich4

c======================================================================
c======================================================================

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

c  -------- Quelques verifications au demarrage sur les tailles des tableaux.
         IF (microfi.ge.1) then
c        Faire de la microphysique sans traceurs... bon courage !
           if (nmicro.le.0) then
             print*,"aLeRtE cRiTiQuE !!!"
             print*,"Vous faites de la microphysique sans traceurs"
             print*,"microphysique..."
             print*,"Je m'arrete et vous laisse reflechir !"
             stop
           endif
c        Nombre de traceurs incompatibles avec la microphysique. 
           if ((nmicro.ne.ntype*nrad).and.(clouds.eq.1)) then
             print*,"aLeRtE cRiTiQuE !!!"
             print*,"Nb trac imcompatible avec la microphysique."
             print*,nmicro,ntype*nrad
             stop
           endif
           if ((nmicro.ne.nrad).and.(clouds.eq.0)) then
             print*,"aLeRtE cRiTiQuE !!!"
             print*,"Nb trac imcompatible avec la microphysique."
             print*,nmicro,nrad
             stop
           endif
         ENDIF

      endif

c RAZ des sorties : les moyennes se font directement dans IOIPSL :
c
          flxesp_i(:,:,:) = 0.
          tau_drop(:,:)   = 0.
          tau_aer(:,:,:)  = 0.
          solesp(:,:,:)   = 0.
          precip(:,:)     = 0.   ! c'est uniquement une sortie en um/s
c
          prec(:,:)       = 0.   ! c'est la variable temporaire des precipitions de la microfi
                                 ! prec est en m (metre precipitable)

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

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

c   -------------------
c   Gestion de la temperature et de la pression :
c   soit la chimie est active, soit la microphysique se fait en 2D.
      IF (chimi.or.microfi.eq.1) THEN
 
        zplev = 0.0
        zplay = 0.0
        zzlev = 0.0
        zzlay = 0.0
        ztemp = 0.0
        zqaer = 0.0
        ychim = 0.0
        zmu0  = 0.0
        zfract= 0.0
        zgaz1 = 0.0
        zgaz2 = 0.0
        zgaz3 = 0.0
        zprec = 0.0
        zflxesp_i = 0.0
        ztau_drop = 0.0 
        ztau_aer  = 0.0
        zsolesp   = 0.0
   
       do l=1,llm+1
         zplev(1,l) = pplev(1,l)
         zzlev(1,l) = pzlev(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
               zzlev(j,l) = zzlev(j,l) + pzlev(ig0+i,l)/iim
	    enddo
	 enddo
         zplev(jjm+1,l) = pplev(klon,l)
         zzlev(jjm+1,l) = pzlev(klon,l)
       enddo

       do l=1,llm
         ztemp(1,l) = ptemp(1,l)
         zplay(1,l) = pplay(1,l)
         zzlay(1,l) = pzlay(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
               zzlay(j,l) = zzlay(j,l) + pzlay(ig0+i,l)/iim
	    enddo
	 enddo
         ztemp(jjm+1,l) = ptemp(klon,l)
         zplay(jjm+1,l) = pplay(klon,l)
         zzlay(jjm+1,l) = pzlay(klon,l)
	 temp_eq  = ztemp((jjm+1)/2,:)
	 press_eq = zplay((jjm+1)/2,:)/100. ! en mbar
       enddo

      ENDIF   ! chimi or microfi=1

c   -----------------------------
c   Gestion des variables de la microphysique :
c
c   -------------------
      if (microfi.ge.1) then

c   Traceurs microphysiques: passage en extensif: n/kg --> n/m^2 (2D ou 3D passage obligatoire)
      DO iq=1,nmicro
c        print*,tname(iq)
        DO l=1,llm
          DO i = 1, klon
            qaer(i,l,iq) = tr_seri(i,l,iq)*delp(i,l)/RG
          ENDDO
        ENDDO
      ENDDO 
c     copie du tableau de traceur :
      qaer0(:,:,:)=qaer(:,:,:)
c
c   -------------------
c    Extraction des gaz pour les nuages
c
c     recuperation des indices des gaz qui nous interesse       
      if (firstcall) then
        if (clouds.eq.1) then
          icldch4=-1
          icldc2h6=-1
          icldc2h2=-1
          do i=1,nqmax
            if (tname(i).eq."CH4") then
              icldch4=i
c              ich4=i
            elseif (tname(i).eq."C2H6") then
              icldc2h6=i
            elseif (tname(i).eq."C2H2") then
              icldc2h2=i
            endif
          enddo
          if (icldch4 .eq.-1 .or. 
     &        icldc2h6.eq.-1 .or.
     &        icldc2h2.eq.-1 ) then
            print*, "Sacrebleu !!!"
            print*, "Vous voulez faire des nuages sans gaz."
            print*, "Mais vous etes inconscient. Je vais m'arreter la"
            print*, "pour vous laisser reflechir au probleme"
            STOP
          endif
        endif    ! clouds=1
      endif      ! firstcall

c     Saturation et fraction molaire CLOUD 
c     Calcul des saturations pour les esp chimique de la muphy des nuages.
c     On le fait ici pour les sortir dans physiq.F sans avoir a surcharger la routine. 
c     Elles passent ensuite dans un common pour passer dans les I/O.
c
c-------------------------------------------
      IF (clouds.eq.1) THEN
        DO l=1,llm
          DO i = 1, klon
            call ch4sat(ptemp(i,l),pplay(i,l),tmp) !tmp en kg/kg !
            satch4(i,l) = tr_seri(i,l,icldch4)/(tmp*28./16.)

            call c2h6sat(ptemp(i,l),pplay(i,l),tmp)
            satc2h6(i,l) =tr_seri(i,l,icldc2h6)/(tmp*28./30.)

            call c2h2sat(ptemp(i,l),pplay(i,l),tmp)
            satc2h2(i,l) =tr_seri(i,l,icldc2h2)/(tmp*28./26.)
   
          ENDDO
        ENDDO

c   Copie des gaz (en 3D)  <== UNIQUEMENT SI ON FAIT DES NUAGES
        gaz1(:,:) = tr_seri(:,:,icldch4)
        gaz2(:,:) = tr_seri(:,:,icldc2h6)
        gaz3(:,:) = tr_seri(:,:,icldc2h2)

      ENDIF      ! clouds=1
        
      endif      ! microfi.ge.1

c     -------------------
c     Si microfi = 1 on est en 2D :
c     conversion des inputs de muphys
      IF (microfi.eq.1) THEN

         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
c     traceurs 3D --> 2D
c 
      do iq=1,nqmax
       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
c       copie du tableau de traceur
        zqaer0(:,:,:) = zqaer(:,:,:)
c 
c      gaz 3D --> 2D    <=== UNIQUEMENT SI ON FAIT DES NUAGES.
c
        if (clouds.eq.1) then
          do l=1,llm
            zgaz1(1,l) = gaz1(1,l)
            zgaz2(1,l) = gaz2(1,l)
            zgaz3(1,l) = gaz3(1,l)
            do j=2,jjm
              ig0=1+(j-2)*iim
              do i=1,iim
                zgaz1(j,l) = zgaz1(j,l) + gaz1(ig0+i,l)/iim
                zgaz2(j,l) = zgaz2(j,l) + gaz2(ig0+i,l)/iim
                zgaz3(j,l) = zgaz3(j,l) + gaz3(ig0+i,l)/iim
              enddo
            enddo
            zgaz1(jjm+1,l) = gaz1(klon,l)
            zgaz2(jjm+1,l) = gaz2(klon,l)
            zgaz3(jjm+1,l) = gaz3(klon,l)
          enddo
 
          zgaz10=zgaz1
          zgaz20=zgaz2
          zgaz30=zgaz3
        endif ! clouds=1

      endif   ! microfi=1

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   en 2D/3D !!!!!!
c    --------------------------

       IF(firstcall) THEN 
        print*,'MICROPHYSIQUE ',MICROFI
       ENDIF

c       call begintime(tt0)
       IF (MICROFI.eq.0) THEN
c        PAS DE MICROPHYSIQUE :
c        On appelle juste rdf pour creer la grille de rayons.
         IF (firstcall) THEN
          print*,'MICROPHYSIQUE OFF-LINE',MICROFI
           call rdf()
         ENDIF
c        NOTES : 
c        L'appel de rdf ne sert a rien ici mis a part pour le TR. Si cet
c        appel a deja lieu dans le TR inutile de le refaire ici.
c        Je ne sais pas exactement comment marche les modules en F90
c        Mais je recopie les valeurs du common/part/ de rdf pour 
c        les mettre dans un common interne a la microphysique (voir varmuphy.h)
c        DONC J'AI BESOIN D'AVOIR ACCES A L'ANCIEN COMMON !!!
c
       ELSEIF (MICROFI.eq.1) THEN
c      MICROPHYSIQUE 2D :
c      Les input/output comportent le prefixe z pour 2D :)
         zdqmufi = 0.  ! ne sert que pour chimi pour condensation
         call muphys(jjm+1,
     &        zplev,zplay,zzlev,zzlay,
     &        ztemp,zqaer,zgaz1,zgaz2,zgaz3,
     &        nmicro,ptimestep,
     &        zmu0,zfract,
c -------- sorties diagnostiques
     &        zflxesp_i,
     &        ztau_drop,ztau_aer,
     &        zsolesp,zprec)
       ELSE
c      MICROPHYSIQUE 3D :
c      Les input sont des champs 3D directement !
         call muphys(klon,
     &        pplev,pplay,pzlev,pzlay,
     &        ptemp,qaer,gaz1,gaz2,gaz3,
     &        nmicro,ptimestep,
     &        pmu0,pfract,     
c ------ sorties diagnostiques
     &        flxesp_i,
     &        tau_drop,tau_aer,
     &        solesp,prec)
c
c    NOTES :
c    Ici toutes nos sorties sont des champs 3D...(meme les diagnostiques)
c    On a rien a faire mis a part copier les dq dans les d_tr
c
       ENDIF
c       call endtime(tt0,tt1)
c       ttmuphys=ttmuphys+tt1

c-----------------------------------------------------------------------
c     Mise a jour des sorties de muphys
c    -------------
c       En 2D on copie les sorties de muphys de la grille LATxALT
c       sur la grille complete. 
       IF (microfi.eq.1) THEN
c        precipitations
         DO l=1,5
           prec(1,l) = zprec(1,l)
           ig0 = 2
           DO j=2,jjm
             DO i = 1, iim
               prec(ig0,l) = zprec(j,l)
               ig0 = ig0 + 1
             ENDDO
           ENDDO
           prec(ig0,l) = zprec(jjm+1,l)
         ENDDO
c        taux sedimentation
         DO l=1,llm
c          taux sed goutte
           IF (clouds.eq.1) THEN
             tau_drop(1,l) = ztau_drop(1,l)
             ig0 = 2
             DO j=2,jjm
               DO i = 1, iim
                 tau_drop(ig0,l) = ztau_drop(j,l)
                 ig0 = ig0 + 1
               ENDDO
             ENDDO
             tau_drop(ig0,l) = ztau_drop(jjm+1,l)
           ENDIF
c          taux sed aer
           DO iq=1,nrad
             tau_aer(1,l,iq)  = ztau_aer(1,l,iq)
             ig0 = 2
             DO j=2,jjm
               DO i = 1, iim
                 tau_aer(ig0,l,iq)  = ztau_aer(j,l,iq)
                 ig0 = ig0 + 1
               ENDDO
             ENDDO
             tau_aer(ig0,l,iq) = ztau_aer(jjm+1,l,iq)
           ENDDO
         ENDDO
c        flux glace / production glace
         IF (clouds.eq.1) THEN
           DO iq=1,3
             DO l=1,llm
               flxesp_i(1,l,iq) = zflxesp_i(1,l,iq)
               solesp(1,l,iq) = zsolesp(1,l,iq)
               ig0 = 2
               DO j=2,jjm
                 DO i = 1, iim
                   flxesp_i(ig0,l,iq)=zflxesp_i(j,l,iq)
                   solesp(ig0,l,iq) = zsolesp(j,l,iq)    
                   ig0 = ig0 + 1
                 ENDDO
               ENDDO
               flxesp_i(ig0,l,iq)=zflxesp_i(jjm+1,l,iq)
               solesp(ig0,l,iq) = zsolesp(jjm+1,l,iq)
             ENDDO
           ENDDO
         ENDIF
       ENDIF
       
c-----------------------------------------------------------------------
c     Gestion des sources
c    -------------
c
       IF (clouds.eq.1) THEN
         IF (microfi.eq.1) THEN
c          On repasse les gaz en 3D si on a fait de la microphysique en 2D
           DO l=1,llm
             gaz1(1,l) = zgaz1(1,l)
             gaz2(1,l) = zgaz2(1,l)
             gaz3(1,l) = zgaz3(1,l)
             ig0 = 2
             DO j=2,jjm
               DO i = 1, iim
                 gaz1(ig0,l) = zgaz1(j,l)* gaz1(ig0,l) /zgaz10(j,l)
                 gaz2(ig0,l) = zgaz2(j,l)* gaz2(ig0,l) /zgaz20(j,l)
                 gaz3(ig0,l) = zgaz3(j,l)* gaz3(ig0,l) /zgaz30(j,l)
                 ig0 = ig0 + 1
               ENDDO
             ENDDO
             gaz1(ig0,l) = zgaz1(jjm+1,l)
             gaz2(ig0,l) = zgaz2(jjm+1,l)
             gaz3(ig0,l) = zgaz3(jjm+1,l)
           ENDDO
         ENDIF
c        Mise a jour du reservoir de CH4 (ie : seul le CH4 remplit le reservoir)
         DO i=1,klon
            reservoir(i) = reservoir(i)+prec(i,1)
         ENDDO 
c        Calcul des sources :
c        ch4=0.
c        ch4(1) = gaz1(1,1)
c         do j=2,jjm
c           ig0=1+(j-2)*iim
c           do i=1,iim
c             ch4(j)= ch4(j) + gaz1(ig0+i,1)/iim
c           enddo
c         enddo
c         ch4(jjm+1) = gaz1(ig0,1)

         CALL sources(klon,klev,ptimestep,z0,
     &                pu,pv,pplev,pzlay,pzlev,
     &                gaz1,gaz2,gaz3,
     &                ftsol,evapch4,reservoir) 
 
c        ch4b=0.
c        ch4b(1) = gaz1(1,1)
c         do j=2,jjm
c           ig0=1+(j-2)*iim
c           do i=1,iim
c             ch4b(j)= ch4b(j) + gaz1(ig0+i,1)/iim
c           enddo
c         enddo
c         ch4b(jjm+1) = gaz1(ig0,1)
c         do j=1,jjm+1
c           write(499,*) j,ch4(j),ch4b(j)
c         enddo
c         write(499,*) ""
       ENDIF
c-----------------------------------------------------------------------
c     Condensation
c    -------------

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

c   tendance (en /s) passee sur zdqmufi(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
           zdqmufi(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 zdqmufi(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          zdqmufi(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          zdqmufi(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          zdqmufi(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 ----- commentaire de fin (mise a jour des profil de fraction molaire)
          
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
c                                        
c ---> pas de microphysique
       IF (microfi.eq.0) THEN
         DO iq=1,nmicro
           d_tr_mph(:,:,iq)=0.
         ENDDO
       ENDIF
c ---> microphysique 2D
      IF (microfi.eq.1) THEN
c  on repasse le champ de traceurs en 3D (pas les tendances)
         DO iq=1,nmicro
           DO l=1,llm
             qaer(1,l,iq) = zqaer(1,l,iq)
             ig0          = 2
             DO j=2,jjm
               DO i = 1, iim
c    un petit patch : 
c    Si la moyenne zonale au depart est "nulle" :
c    On a quand meme le droit de produire des traceurs dans la cellule.
c    On considere donc que la valeur de sortie 3D correspond a la valeur de sortie 2D.
c    Cela permet aussi entre autre d'eviter les NaN pour les traceurs des nuages !
c    (au dessus de la tropo pas de nuages donc qaer(nrad+1:ntype*nrad) = 0 !!!)
                 IF (zqaer0(j,l,iq).lt.1e-100) THEN
                   qaer(ig0,l,iq) = zqaer(j,l,iq) 
                 ELSE
                   qaer(ig0,l,iq) = zqaer(j,l,iq) *
     &             qaer0(ig0,l,iq)/zqaer0(j,l,iq)
                 ENDIF
                 ig0 = ig0 + 1
               ENDDO
             ENDDO
             qaer(ig0,l,iq) = zqaer(jjm+1,l,iq)
           ENDDO
         ENDDO
c        La tendances correspond a (qaer-qaer0)/ptimestep
         DO iq=1,nmicro
           DO i=1,klon
             DO l=1,llm
               d_tr_mph(i,l,iq) = (qaer(i,l,iq)-qaer0(i,l,iq))/
     &                            ptimestep
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
c ---> microphysique 3D
       ELSEIF(microfi.gt.1) THEN 
         DO iq=1,nmicro
           DO l=1,llm
             DO i = 1, klon
               d_tr_mph(i,l,iq)=(qaer(i,l,iq)-qaer0(i,l,iq))/ptimestep
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 zdqmufi 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) = zdqmufi(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)  = zdqmufi(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) = zdqmufi(jjm+1,l,iq)
         ENDDO
      ENDDO

      endif   ! chimi

c--------------------------------------------------
c  CONDENSATION VIA MICROFI
c----------------------
c La microphysique avec nuages doit se faire obligatoirement en 3D.  (FAUX ACTUELLEMENT)
c Rien n'empeche de faire la chimie en 2D. Cependant pour prendre en compte la 
c condensation due a la microfi (en 3D) on recalcule la tendance finale pour
c les especes concernees (CH4, C2H6 pour le moment). 
       IF (microfi.ge.1.and.clouds.eq.1) THEN
         DO i=1,klon
           DO l=1,klev
c     condensation CH4
             d_tr_mph(i,l,icldch4)=(gaz1(i,l)-tr_seri(i,l,icldch4))
     &                            /ptimestep
c     condensation C2H6
             d_tr_mph(i,l,icldc2h6)=(gaz2(i,l)-tr_seri(i,l,icldc2h6))
     &                             /ptimestep
c     condensation C2H2
             d_tr_mph(i,l,icldc2h2)=(gaz3(i,l)-tr_seri(i,l,icldc2h2))
     &                             /ptimestep
           ENDDO 
         ENDDO
       ENDIF
c        ch4c=0.
c       do l=1,llm
c       ch4c(1,l) = tr_seri(1,l,icldch4)
c       do j=2,jjm
c         ig0=1+(j-2)*iim
c         do i=1,iim
c            ch4c(j,l)= ch4c(j,l)+tr_seri(ig0+i,l,icldch4)/iim
c         enddo
c       enddo
c        ch4c(jjm+1,l) = tr_seri(klon,l,icldch4)
c      enddo
c       do l=1,llm
c         write(500,*) pplay(25,l),ch4c(25,l)
c       enddo
c       write(500,*) ""


c--------------------------------------------------
c  MISE A JOUR CH4 : (pour refixer la fraction 
c                     molaire)
c--------------------------------------------------
c       IF (firstcall) THEN
c         do i=1,klon
c           do j=1,llm
c             call ch4sat(ptemp(i,j),pplay(i,j),tmp) !tmp en kg/kg !
c             tmp=0.95*0.85*tmp*28./16.
c             if (pplay(i,j).lt.20000.) then
c               dqch4 = 1.4e-2          
c             else
c               dqch4 = tmp
c             endif
c             d_tr_mph(i,j,icldch4)=(-tr_seri(i,j,icldch4)+dqch4)/ 
c     &       ptimestep
c           enddo
c         enddo
c         
c       ENDIF

c--------------------------------------------------
c  CONVERSION PRECIPITATION : 
c  en microns/secondes
c--------------------------------------------------
        precip = prec * 1.e6 / ptimestep


c--------------------------------------------------
c CALCUL DU FLUX DE CHALEUR LATENTE D'EVAPORATION 
c DU METHANE 
c--------------------------------------------------
       IF (clouds.eq.1) THEN
         DO i=1,klon
           fte= (1.-ftsol(i)/305.5)
           ftm= (1.-ftsol(i)/190.5)
           if(ftm.le.1.e-3) ftm=1.e-3
           if(fte.le.1.e-3) fte=1.e-3
           Lvch4 =8.314*190.4*
     &     (7.08*ftm**0.354+10.95*1.1e-2*ftm**0.456)
     &     /mch4
           ! evapch4 en m3/m2 {ok}
           ! 425 en kg/m3
           ! Lv en J/kg       {ok}
           ! ptimestep en s   {ok}
           fclat(i)=(evapch4(i)*Lvch4*rhoi_ch4)   ! en J/m2/s
         ENDDO
       ENDIF

c--------------------------------------------------
c  GESTION DES RAYONS DE GOUTTES POUR TR
c--------------------------------------------------
       IF (clouds.eq.1) THEN

c Calcul du rayon des gouttes par bin ...
c----------------------------------------
         DO i=1,klon
           DO j=1,klev 
             DO iq=1,nrad
*      Rayon minimum selon la quantité de noyaux
               IF (qaer(i,j,iq+nrad) .le.   1.e-5) THEN
                  rcloud(i,j,iq) = 1.e-10
               ELSE
                 rcloud(i,j,iq)=
     &           ((qaer(i,j,iq+2*nrad)/qaer(i,j,iq+nrad)+
     &           qaer(i,j,iq+3*nrad)/qaer(i,j,iq+nrad) +
     &           v_e(iq))*0.75/RPI)**(1./3.)
               ENDIF 
             ENDDO
           ENDDO
         ENDDO

c .... et de leur rayon moyen total (tt bins confondu)
c------------------------------------------------------
         DO i=1,klon
           socccld=0.
           DO j=klev,1,-1    !de haut en bas pour le calcul des opacites
             vcl=0.
             nuc=0.
             xgsn=0.
             xmsn=0.
             xesn=0.
             xasn=0.
             DO iq=1,nrad
               vcl=vcl+qaer(i,j,iq+2*nrad)+
     &         qaer(i,j,iq+3*nrad)+ 
     &         qaer(i,j,iq+4*nrad)+ 
     &         v_e(iq)*qaer(i,j,iq+nrad)            ! volume des gouttes
               nuc=nuc+qaer(i,j,iq+nrad)            ! nombre de noyaux
               xgsn=xgsn+qaer(i,j,iq+nrad)*v_e(iq)  ! volume de noyaux
               xmsn=xmsn+qaer(i,j,iq+2*nrad)        ! volume de methane
               xesn=xesn+qaer(i,j,iq+3*nrad)        ! volume d' ethane
               xasn=xasn+qaer(i,j,iq+4*nrad)        ! volume d' acethylene
             ENDDO 
             IF (nuc .le.  1.e-5) THEN
               rmcloud(i,j)=1.0e-10
               xfrac(i,j,:)=0.
             ELSE
               IF(xgsn/vcl.lt.0.  .or. xgsn/vcl.gt.1.001)
     &         print*, 'PB AVEC XFRAC:', i,j,xgsn,vcl 
               rmcloud(i,j)=          ! rayon moyen des gouttes
     &         (vcl/nuc*0.75/RPI)**(1./3.)
               xfrac(i,j,1)=xgsn/vcl         ! fraction volumique noyau/goutte
               xfrac(i,j,2)=xmsn/vcl         ! fraction volumique CH4/goutte
               xfrac(i,j,3)=xesn/vcl         ! fraction volumique C2H6/goutte
               xfrac(i,j,4)=xasn/vcl         ! fraction volumique C2H2/goutte
c              calcul du rayon moyen (moyenne temporelle)
               rmcbar(i,j)=rmcbar(i,j)+rmcloud(i,j)
               xfbar(i,j,:)=xfbar(i,j,:)+xfrac(i,j,:)
               ncount(i,j) = ncount(i,j)+1
             ENDIF
             socccld=socccld+RPI*(rmcloud(i,j)**2.)*nuc
             occcld(i,j)=socccld
           ENDDO
         ENDDO
c
c      OCCCLD
c      Calcul le nombre d'occurence d'un nuage 
c      d'opacité comprise en kmin et kmax
c          k        kmin            kmax
c          1   0.0000000      0.10000000    
c          2  0.10000000      0.17782794    
c          3  0.17782794      0.31622776    
c          4  0.31622776      0.56234139    
c          5  0.56234139       1.0000000    
c          6   1.0000000       1.7782795    
c          7   1.7782795       3.1622777    
c          8   3.1622777       5.6234136    
c          9   5.6234136       10.000000    
c         10   10.000000       17.782795    
c         11   17.782795       31.622778    
c         12   31.622778       100000.00
c
c        mise a zero de occld_m
         occcld_m=0.
         DO i=1,klon
           DO j=1,klev
             DO k=1,12
               ex=10.**(0.25)
               kmin=0.
               kmax=1.e5
               if(k.ne.1)  kmin=0.1*ex**(k-2)
               if(k.ne.12) kmax=0.1*ex**(k-1)
               if(occcld(i,j).ge.kmin .and. occcld(i,j).lt.kmax)
     &         occcld_m(i,j,k)=1.
             ENDDO
           ENDDO
         ENDDO
       ENDIF  ! fin condition clouds => pas besoin de calculer des rayons

      RETURN
      END