!
! $Header: /home/cvsroot/LMDZ4/libf/phylmd/phytrac.F,v 1.16 2006/03/24 15:06:23 lmdzadmin Exp $
!
c
c
      SUBROUTINE  phytrac_chimie (
     I                    debutphy,
     I                    gmtime,
     I                    nqmax,
     I                    n_lon,
     I                    lat,
     I                    lon,
     I                    n_lev, 
     I                    pdtphys,
     I                    temp,
     I                    pplay,
     O                    trac)
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 chemparam_mod
      use conc, only: mmean
      IMPLICIT none
      
#include "clesphys.h" 
#include "YOMCST.h"
c======================================================================
c Arguments:
c
c
c   EN ENTREE:
c   ==========
c
c   divers:
c   -------
c

      REAL  sza_local 
      REAL  gmtime
c      INTEGER, SAVE :: cpt_cloudIO  !un compteur pour fichier sortie cloud_parameter en 1D
      INTEGER  iq
      INTEGER  i
      INTEGER  ilon, ilev
      integer  n_lon  ! nombre de points horizontaux
      INTEGER  n_lev  ! nombre de couches verticales
      INTEGER  nqmax ! nombre de traceurs auxquels on applique la physique

      real  pdtphys  ! pas d'integration pour la physique (seconde)
      real  lat(n_lon), lat_local(n_lon)
      real  lon(n_lon), lon_local(n_lon)
      real  temp(n_lon,n_lev) ! temp
      real  trac(n_lon,n_lev,nqmax) ! traceur
      real  trac_sav(n_lon,n_lev,nqmax)
      real  trac_sum(n_lon,n_lev)
      real  pplay(n_lon,n_lev)  ! pression pour le mileu de chaque couche (en Pa)
      real  lon_sun

      logical debutphy       ! le flag de l'initialisation de la physique

C     Auxilary variables:
 
      REAL, DIMENSION(n_lon,n_lev) :: mrtwv,mrtsa,
     +                                mrwv,mrsa
     
C    ps_sa: satur pressure pure SA 
C    satps_sa: satur pres over mixture in dyne/cm2=Pa/10
C----------------------------------------------------------------------------
        
C Variables liees a l'ecriture de la bande histoire : phytrac.nc
 
      logical ok_sync
      parameter (ok_sync = .true.)

c      modname = 'phytrac'

c      PRINT*,'DEBUT subroutine PHYTRAC'

c----------------------------------
c debutphy: Initiation des traceurs
c----------------------------------

      if (debutphy) then
      
         if (n_lon .EQ. 1) then           
         PRINT*,'n_lon 1D: ',n_lon
         end if

c=============================================================
c == CASE INIT PHOTOCHEM ==
c=============================================================
                             
       IF (reinit_trac.and.ok_chem) THEN
       PRINT*,'REINIT MIXING RATIO TRACEURS'

c					Passage de Rm a Rv
c	=============================================================
c       Necessaire si on reprend les start.nc qui sont en MMR
c SEULEMENT LES GAZ
         DO iq=1,nqmax-nmicro
          trac(:,:,iq)=trac(:,:,iq)*mmean(:,:)/M_tr(iq)
         END DO
c	=============================================================
          
c		Initialisation des profils traceurs en Rv
c=============================================================
c initialisation sert a mettre les valeurs voulues par utilisateur pour
c chaque traceur
c exemple: trac(ilon,ilev,q)=xx

c     trac_sav sert a sauver les valeurs initiales du start.nc      
      trac_sav=trac

c     On initialise les traceurs a zero obligatoire pour la chimie
      trac(:,:,:)=1.0E-30

c     !!! Verif traceurs !!!
      if (   (i_ocs.ne.0).and.(i_co.ne.0).and.(i_hcl.ne.0) 
     &  .and.(i_so2.ne.0).and.(i_h2o.ne.0).and.(i_n2.ne.0)
     &  .and.(i_co2.ne.0) ) then
       trac(:,1:22,i_ocs)=3.E-6
       trac(:,1:22,i_co)=25.E-6
       trac(:,:,i_hcl)=0.4E-6
       trac(:,1:22,i_so2)=10.E-6
       trac(:,1:22,i_h2o)=30.0E-6

c     remettre tous les traceurs du start => trac(:,:,:)=trac_sav(:,:,:)

c     N2 n est pas encore une espece chimique du modele chimique
c     traceur passif pour la chimie-transport
       trac(:,:,i_n2)=0.35d-1
   
!!!! GG:   Initialization CO2 = 1 - qtot 
!!    It assures that vmr_tot = 1
c     On a donc le CO2 qui est le restant d atmosphere Venus  
        trac_sum(:,:)=0.0
c SEULEMENT LES GAZ
        DO iq=2,nqmax-nmicro
         trac_sum(:,:)= trac_sum(:,:) + trac(:,:,iq) 
        END DO

        trac(:,:,i_co2)= 1-trac_sum(:,:)
      else
        write(*,*) "Réinitialisation des traceurs avec chimie "
        write(*,*) "IL manque un traceur pour les options choisies"
        stop
      endif ! verif traceurs
        
c=============================================================
      
c					Passage de Rv a Rm
c	=============================================================
c SEULEMENT LES GAZ
         DO iq=1,nqmax-nmicro
          trac(:,:,iq)=trac(:,:,iq)*M_tr(iq)/mmean(:,:)
         END DO
c	=============================================================
    
       ENDIF  !FIN REINIT TRAC

c=============================================================
c == CASE INIT GAZ when muphy without chemistry ==
c=============================================================
                             
       if (.not.ok_chem.and.ok_cloud.and.(cl_scheme.eq.2)) then

c					Passage de Rm a Rv
c	=============================================================
c       Necessaire si on reprend les start.nc qui sont en MMR
c SEULEMENT LES GAZ
         DO iq=1,nqmax-nmicro
          trac(:,:,iq)=trac(:,:,iq)*mmean(:,:)/M_tr(iq)
         END DO
c	=============================================================
          
        call vapors4muphy_ini(n_lon,n_lev,trac)

c					Passage de Rv a Rm
c	=============================================================
c SEULEMENT LES GAZ
         DO iq=1,nqmax-nmicro
          trac(:,:,iq)=trac(:,:,iq)*M_tr(iq)/mmean(:,:)
         END DO
c	=============================================================
    
      endif
	
c-------------
c fin debutphy
c-------------

      ENDIF  ! fin debutphy 

c	=============================================================
c					Passage de Rm a Rv
c                                        DES GAZ
c	=============================================================
      DO iq=1,nqmax-nmicro
        trac(:,:,iq)=MAX(trac(:,:,iq)*mmean(:,:)/M_tr(iq),1.E-30)
      END DO
c	=============================================================


c	=============================================================
c	    Appel Microphysique (schema simple, these Aurelien Stolzenbach)
c	=============================================================

      IF (ok_cloud .AND. cl_scheme.eq.1) THEN

c	  Passage de Rm a Rv pour les liq
         trac(:,:,i_h2so4liq)=MAX(trac(:,:,i_h2so4liq)
     &                          *mmean(:,:)/M_tr(i_h2so4liq),1.E-30)
         trac(:,:,i_h2oliq)=MAX(trac(:,:,i_h2oliq)
     &                          *mmean(:,:)/M_tr(i_h2oliq),1.E-30)
              
c      PRINT*,'DEBUT CLOUD'      
c     On remet tout le RM liq dans la partie gaz
c     !!! On reforme un nuage a chaque fois !!!

      DO ilev=1, n_lev
      DO ilon=1, n_lon         
      mrtwv(ilon,ilev)=trac(ilon,ilev,i_h2o) +
     &  trac(ilon,ilev,i_h2oliq)
      mrtsa(ilon,ilev)=trac(ilon,ilev,i_h2so4) +
     &  trac(ilon,ilev,i_h2so4liq)
      mrwv(ilon,ilev)=mrtwv(ilon,ilev)
      mrsa(ilon,ilev)=mrtsa(ilon,ilev)
      ENDDO
      ENDDO
	            
      CALL new_cloud_venus(n_lev, n_lon,
     e temp,pplay,
     e mrtwv,mrtsa,
     e mrwv,mrsa)

c	=========================================                
c	Actualisation des mixing ratio liq et gaz
c	=========================================
c	Si la routine new_cloud_venus n'a pas actualise mrwv et mrsa
c	on a alors bien mr=mrt pour sa et wv, donc les parties liq sont=0 hors du nuage
c	ou si on ne condense pas

c      PRINT*,'DEBUT ACTUALISATION OUTPUT CLOUD'
c    si tout se passe bien, mrtwv et mrtsa ne changent pas
      DO ilev=1, n_lev
      DO ilon=1, n_lon       
      trac(ilon,ilev,i_h2o) = mrwv(ilon,ilev)
      trac(ilon,ilev,i_h2oliq) = mrtwv(ilon,ilev) -
     &  trac(ilon,ilev,i_h2o)
     
      trac(ilon,ilev,i_h2so4) = mrsa(ilon,ilev)
      trac(ilon,ilev,i_h2so4liq) = mrtsa(ilon,ilev) -
     &  trac(ilon,ilev,i_h2so4)
      ENDDO
      ENDDO

c	=============================================================
c      PRINT*,'FIN CLOUD, scheme 1'
      ENDIF  

c	=============================================================
c	    Appel Microphysique (schema complet, these Sabrina Guilbon)
c	=============================================================

      IF (ok_cloud .AND. cl_scheme.eq.2) THEN

c	Boucle sur grille (n_lon) et niveaux (n_lev)
        DO ilon=1, n_lon       
         DO ilev=1, n_lev

           CALL MAD_MUPHY(pdtphys,                               ! Timestep
     &  temp(ilon,ilev),pplay(ilon,ilev),                        ! Temperature and pressure
     &  trac(ilon,ilev,i_h2o),trac(ilon,ilev,i_h2so4),           ! Mixing ratio of SA and W
     &  trac(ilon,ilev,i_m0_aer),trac(ilon,ilev,i_m3_aer),       ! Moments of aerosols
     &  trac(ilon,ilev,i_m0_mode1drop),trac(ilon,ilev,i_m0_mode1ccn),  ! Moments of mode 1
     &  trac(ilon,ilev,i_m3_mode1sa),trac(ilon,ilev,i_m3_mode1w),      ! Moments of mode 1
     &  trac(ilon,ilev,i_m3_mode1ccn),                                 ! Moments of mode 1
     &  trac(ilon,ilev,i_m0_mode2drop),trac(ilon,ilev,i_m0_mode2ccn),  ! Moments of mode 2
     &  trac(ilon,ilev,i_m3_mode2sa),trac(ilon,ilev,i_m3_mode2w),      ! Moments of mode 2
     &  trac(ilon,ilev,i_m3_mode2ccn))                                 ! Moments of mode 2

         ENDDO
        ENDDO


c	=============================================================
c      PRINT*,'FIN CLOUD, scheme 2'
      ENDIF  

            
c=============================================================
c		CHIMIE: Boucle sur les lon, lat (n_lon)
c=============================================================

c     AS:
c     Ici, la longitude au midi local se deplace vers l'Ouest
c     c'est le sens terrestre
c     pour Venus on prend juste l'oppose de la longitude et on a la rotation 
c     de Venus et donc le midi local qui se deplace vers l'Est
     
      lon_sun = (0.5 - gmtime) * 2.0 * RPI
      lon_local = lon * RPI/180.0E+0
      lat_local = lat * RPI/180.0E+0
       
      IF (ok_chem) THEN

c      PRINT*,'DEBUT CHEMISTRY'
      DO ilon=1, n_lon

c     calcul sza_local pour obtenir des sza_local > 90, utile pour la chimie
      sza_local = acos(cos(lat_local(ilon))*cos(lon_local(ilon))*
     & cos(lon_sun) + cos(lat_local(ilon))*sin(lon_local(ilon))
     & *sin(lon_sun))* 180.0E+0/RPI 
      
c      PRINT*,'sza_local :', sza_local     
    
c	=============================================================
c					Appel Photochimie
c	=============================================================
c     Pression en hPa => pplay/100.
	
      CALL new_photochemistry_venus(n_lev, n_lon, pdtphys,
     e                         pplay(ilon,:)/100.,
     e                         temp(ilon,:),
     e                         trac(ilon,:,:),
     e                         mmean(ilon,:),
     e                         sza_local, nqmax)
c	=============================================================

      END DO
c      PRINT*,'FIN CHEMISTRY'
    
      END IF
c=============================================================

c	=============================================================
c					Passage de Rv a Rm
c	=============================================================
c                                        DES GAZ
      DO iq=1,nqmax-nmicro
		trac(:,:,iq)=trac(:,:,iq)*M_tr(iq)/mmean(:,:)

      END DO

      IF (ok_cloud .AND. cl_scheme.eq.1) THEN
c	  Passage de Rm a Rv pour les liq
         trac(:,:,i_h2so4liq)=trac(:,:,i_h2so4liq)*M_tr(i_h2so4liq)
     &                          /mmean(:,:)
         trac(:,:,i_h2oliq)=trac(:,:,i_h2oliq)*M_tr(i_h2oliq)
     &                          /mmean(:,:)
      ENDIF
      
c	=============================================================	
C      PRINT*,'FIN PHYTRAC'
      RETURN
      END