!
! $Header: /home/cvsroot/LMDZ4/libf/phylmd/phytrac.F,v 1.16 2006/03/24 15:06:23 lmdzadmin Exp $
!
c
c
      SUBROUTINE phytrac_emiss (timesimu,
     I                    debutphy,
     I                    lafin,
     I                    nqmax,
     I                    nlon,
     I                    nlev, 
     I                    pdtphys,
     I                    paprs,
     I                    xlat,xlon,
     O                    tr_seri)

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 SL: Janvier 2014
c Version developed for surface emission
c Maybe could be used just to compute the 'source' variable from physiq
c
c======================================================================
      USE infotrac_phy, ONLY: nqtot
      use dimphy
      USE geometry_mod, only: cell_area
      USE chemparam_mod,only:M_tr
      USE mod_grid_phy_lmdz, ONLY: nbp_lon, nbp_lat
      IMPLICIT none
#include "YOMCST.h"
#include "clesphys.h"
c======================================================================

c Arguments:

c   EN ENTREE:
c   ==========

      real timesimu   ! duree depuis debut simu (s)
      logical debutphy       ! le flag de l'initialisation de la physique
      logical lafin          ! le flag de la fin de la physique
      integer nqmax ! nombre de traceurs auxquels on applique la physique
      integer nlon  ! nombre de points horizontaux
      integer nlev  ! nombre de couches verticales
      real pdtphys  ! pas d'integration pour la physique (seconde)
      real paprs(nlon,nlev+1)  ! pression pour chaque inter-couche (en Pa)
      REAL xlat(nlon)       ! latitudes pour chaque point 
      REAL xlon(nlon)       ! longitudes pour chaque point 

c   EN ENTREE/SORTIE:
c   =================

      real tr_seri(nlon,nlev,nqmax) ! traceur  

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

c pour emission volcan
      real :: deltatr(klon,klev,nqtot)

      integer,parameter :: nblat=5,nblon=4,nbz=3
      integer,parameter :: Nemiss=1     ! duree emission (Ed)
      integer,save :: Nheight(nbz)      ! layer emission
      real,save :: so2_quantity         ! quantity so2 (kg)
      real,save :: lat_volcan(nblat),lon_volcan(nblon)
      real,save :: area_emiss(nblat,nblon)
      integer,save :: ig_volcan(nblat,nblon)

      INTEGER i, k, it
      integer ilat,ilon,iz
      real    deltalat,deltalon
c======================================================================

c EMISSION TRACEURS 

c---------
c debutphy
c---------      
      if (debutphy) then

        print*,"DEBUT PHYTRAC"
        print*,"PHYTRAC: EMISSION"

        ALLOCATE(M_tr(nqtot))
        M_tr(:)=64.                 ! SO2
	
C=========================================================================
c Caracteristiques des traceurs emis:
C=========================================================================

c nombre total de traceur
         if (nbz*nblat*nblon .gt. nqtot) then
            print*, nbz*nblat*nblon, nqtot
            write(*,*) "Attention, pas assez de traceurs"
            write(*,*) "le dernier sera bien le dernier"
         endif

c quantite en kg
         so2_quantity = 20.*10.**9.

c height (in layer index)
         Nheight(1) =  6  ! ~ 1 km
         Nheight(2) = 16  ! ~ 25 km
         Nheight(3) = 24  ! ~ 50 km

c localisation volcan
         lat_volcan(1) =  70.
         lat_volcan(2) =  35.
         lat_volcan(3) =   0.
         lat_volcan(4) = -35.
         lat_volcan(5) = -70.
         lon_volcan(1) = -125.
         lon_volcan(2) =  -35.
         lon_volcan(3) =   55.
         lon_volcan(4) =  145.
	 
	 ig_volcan(ilat,ilon)= 0
         if ((nbp_lon*nbp_lat)==1) then ! running a 1D simulation
           deltalat=180.
           deltalon=360.
         else
           deltalat = 180./(nbp_lat-1)
           deltalon = 360./nbp_lon
         endif

         do i=1,nlon
          do ilat=1,nblat
           do ilon=1,nblon
            if ((xlat(i).ge.lat_volcan(ilat))
     &     .and.((xlat(i)-deltalat).lt.lat_volcan(ilat))
     &     .and.(xlon(i).le.lon_volcan(ilon))
     &     .and.((xlon(i)+deltalon).gt.lon_volcan(ilon)) ) then
             ig_volcan(ilat,ilon)= i
             area_emiss(ilat,ilon) = cell_area(i)
	     print*,"Lat,lon=",ilat,ilon," OK"
            end if
           end do
          end do
         end do

c Reinit des traceurs si necessaire
         if (reinit_trac) then
	   tr_seri(:,:,:)=0.
	 endif
	 
C=========================================================================
C=========================================================================
      ENDIF  ! fin debutphy 
c-------------
c fin debutphy
c-------------

c======================================================================
c Emission d'un traceur pendant un certain temps
c necessite raz_date=1 dans run.def
c et reinit_trac=y
c======================================================================
       deltatr(:,:,:) = 0.

c source appliquee pendant Nemiss Ed
       if (timesimu .lt. 86400*Nemiss) then

c emet les traceurs qui sont presents sur la grille
        do ilat  = 1,nblat
        do ilon  = 1,nblon
         if (ig_volcan(ilat,ilon).ne.0) then
	 
          do iz = 1,nbz
           it=min( (iz-1)*nblat*nblon+(ilat-1)*nblon+ilon , nqtot )	    
           i=ig_volcan(ilat,ilon)
	   
c injection dans une seule cellule:
c source en kg/kg/s
c            deltatr(i,Nheight(iz),it) = so2_quantity/(86400.*Nemiss) ! kg/s
c     $ *RG/( area_emiss(ilat,ilon)
c     $      *(paprs(i,Nheight(iz))-paprs(i,Nheight(iz)+1)) )    ! /kg (masse cellule)
     
c            tr_seri(i,Nheight(iz),it) = tr_seri(i,Nheight(iz),it)
c     $      + deltatr(i,Nheight(iz),it)*pdtphys

c injection dans toute la colonne (a faire):
            do k=1,Nheight(iz)
             deltatr(i,k,it) = so2_quantity/(86400.*Nemiss) ! kg/s
     $  *RG/( area_emiss(ilat,ilon)
     $       *(paprs(i,1)-paprs(i,Nheight(iz)+1)) )    ! /kg (masse colonne)
     
             tr_seri(i,k,it) = tr_seri(i,k,it)+deltatr(i,k,it)*pdtphys
	    end do
	    
          end do
	  
         endif  ! ig_volcan!=0
        end do
        end do

       end if  ! duree emission
       
c======================================================================
c======================================================================

      RETURN
      END