source: trunk/LMDZ.VENUS/libf/phyvenus/phytrac_chimie.F @ 1443

!
! $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                    pplev,
     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======================================================================
c      USE ioipsl
c      USE infotrac
c      USE control_mod
c      USE dimphy
c      USE comgeomphy
      USE chemparam_mod
      use conc, only: mmean
      IMPLICIT none
      
c#include "dimensions.h"
#include "clesphys.h" 
c#include "temps.h"
c#include "paramet.h"
c#include "comcstfi.h" !me permet de recuperer mugaz et d'autres constantes comme rad,pi etc
#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  pplev(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         if ((n_lon .GT. 1) .AND. ok_chem) then
c !!! DONC 3D !!!       
c           CALL chemparam_ini() 
c         endif
         
c         if ((n_lon .GT. 1) .AND. ok_cloud) then
c !!! DONC 3D !!!
c         CALL cloud_ini(n_lon,n_lev)
c         endif
           
       IF (reinit_trac) THEN
       PRINT*,'REINIT MIXING RATIO TRACEURS'

c       =============================================================
c                                       Passage de Rm à Rv
c       =============================================================
c       Necessaire si on reprend les start.nc qui sont en MMR

         DO iq=1,nqmax
          trac(:,:,iq)=trac(:,:,iq)*mmean(:,:)/M_tr(iq)
         END DO
c       =============================================================
          
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     !!! AVEC NUAGE !!!
      trac(:,1:22,i_ocs)=3.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
        DO iq=2,nqmax
         trac_sum(:,:)= trac_sum(:,:) + trac(:,:,iq) 
        END DO

        trac(:,:,i_co2)= 1-trac_sum(:,:)
        
c=============================================================
      
c       =============================================================
c                                       Passage de Rv à Rm
c       =============================================================
         DO iq=1,nqmax
          trac(:,:,iq)=trac(:,:,iq)*M_tr(iq)/mmean(:,:)
         END DO
c       =============================================================
    
       ENDIF  !FIN REINIT TRAC

        
c-------------
c fin debutphy
c-------------

      ENDIF  ! fin debutphy 

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


c       =============================================================
c                        Appel Microphysique (sans nucleation)
c                  Volume Mixing Ratio
c       =============================================================

         IF (ok_cloud) THEN
               
c      PRINT*,'DEBUT CLOUD'      
c     On remet tout le RM liq dans la partie gaz
c     !!! On reforme un nuage à 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,pplev,
     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 actualisé 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'
      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 Vénus on prend juste l'opposé de la longitude et on a la rotation 
c     de Vénus et donc le midi local qui se déplace 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
       
      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     
    
      IF (ok_chem) THEN
c      PRINT*,'DEBUT CHEMISTRY'
c       =============================================================
c                                       Appel Photochimie
c       =============================================================
c     Pression en hPa => pplev/100.
        
      CALL new_photochemistry_venus(n_lev, n_lon, pdtphys,
     e                         pplev(ilon,:)/100.,
     e                         temp(ilon,:),
     e                         trac(ilon,:,:),
     e                         sza_local, nqmax)
c       =============================================================
c      PRINT*,'FIN CHEMISTRY'
    
        END IF

      END DO
c       =============================================================
c                                       Passage de Rv à Rm
c       =============================================================
        DO iq=1,nqmax
c               trac(:,:,iq)=trac(:,:,iq)*M_tr(iq)/RMD
                trac(:,:,iq)=trac(:,:,iq)*M_tr(iq)/mmean(:,:)

        END DO
c       =============================================================   
C      PRINT*,'FIN PHYTRAC'
      RETURN
      END