source: trunk/LMDZ.VENUS/libf/phyvenus/phytrac_emiss.F @ 3556

!
! $Header: /home/cvsroot/LMDZ4/libf/phylmd/phytrac.F,v 1.16 2006/03/24 15:06:23 lmdzadmin Exp $
!
c
c
      SUBROUTINE phytrac_emiss (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: tname, nqtot
#ifdef CPP_XIOS      
      use xios_output_mod, only:  send_xios_field
#endif
      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   ==========

      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 type effusion
      real :: deltatr(klon,klev,nqtot)




      integer,parameter :: nblat=3,nblon=3,nbaire=3,nbflux=2,maxcell=16 
      integer,parameter :: Nheight=3 ! layer emission (150m)
      integer,save :: Ncell(nbaire) 
      real,save :: flux_surface_co2(nbflux) ! flux de CO2 emis (kg/m2/s) 
      real :: flux(nlon,nqtot) 
      real,save :: lat_zone(nblat),lon_zone(nblon) 
      integer,save :: ig_zone(nblat,nblon,nbaire,nbflux,maxcell) 
      integer,save :: numcell(nblat,nblon,nbaire,nbflux)
        

      INTEGER i, k, it
      integer ilat,ilon,iaire,iflux,ipos
      real    deltalat,deltalon
      


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

c EMISSIONS TRACEURS 

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

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

        ALLOCATE(M_tr(nqtot))
        M_tr(:)=44.                ! CO2

C=========================================================================
c Caracteristiques des traceurs emis:
C=========================================================================

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



c flux de CO2 (kg/s/m2)
         flux_surface_co2(1) = 5.*10.**-9.
         flux_surface_co2(2) = 5.*10.**-15.

c nombre de cellule pour le cote du carre d'aire 
         Ncell(1)= 2
         Ncell(2)= 3
         Ncell(3)= 4
      


c localisation zone emission
         lat_zone(1) = 08.
         lat_zone(2) = -50.
         lat_zone(3) = 35.
         lon_zone(1) = -172.
         lon_zone(2) = -20.
         lon_zone(3) = 70.

 
         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

         numcell(:,:,:,:)=0
         ig_zone(:,:,:,:,:)=0
         do i=1,nlon
          do ilat=1,nblat
           do ilon=1,nblon
            do iaire=1,nbaire

             if ((xlat(i).ge.lat_zone(ilat))
     &      .and.((xlat(i)-Ncell(iaire)*deltalat)
     &      .lt.lat_zone(ilat))
     &      .and.(xlon(i).le.lon_zone(ilon))
     &      .and.((xlon(i)+Ncell(iaire)*deltalon)
     &      .gt.lon_zone(ilon))) then

              do iflux=1,nbflux
         numcell(ilat,ilon,iaire,iflux)=numcell(ilat,ilon,iaire,iflux)+1
        ig_zone(ilat,ilon,iaire,iflux,numcell(ilat,ilon,iaire,iflux))= i
             print*,"Lat,lon,naire,nflux,nlon=",ilat,ilon,iaire,iflux
     &        ,i," OK"
              end do

             end if

            end do
           end do
          end do
         end do

c Reinit des traceurs si necessaire
         if (reinit_trac) then
          tr_seri(:,:,:)=0.


           do i=1,klon
            do k=1,klev
              tr_seri(i,k,:)=1.-28/43.44*max(min(0.035,
     &           0.035*(1.-log(paprs(i,k)/6.e6)/log(9.e6/6.e6))),0.)
            end do
           end do
         endif
 
C=========================================================================
C=========================================================================
      ENDIF  ! fin debutphy 
c-------------
c fin debutphy
c-------------

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

c emet les traceurs qui sont presents sur la grille
      do ilat  = 1,nblat
       do ilon  = 1,nblon
        do iaire = 1,nbaire
         do iflux = 1,nbflux
     
              it=min( (ilat-1)*nblon*nbflux*nbaire+(iaire-1)*nbflux
     &         +(ilon-1)*nbaire*nbflux+iflux , nqtot )   
 

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 ipos=1,maxcell
 
           if (ig_zone(ilat,ilon,iaire,iflux,ipos).ne.0) then

             i=ig_zone(ilat,ilon,iaire,iflux,ipos)         
             flux(i,it)=flux_surface_co2(iflux)  
              
            do k=1,Nheight
             deltatr(i,k,it) = flux_surface_co2(iflux) ! kg/s/m2
     $         *RG/(paprs(i,1)-paprs(i,Nheight+1))    ! /kg (masse colonne)
     
               tr_seri(i,k,it) = tr_seri(i,k,it)+deltatr(i,k,it)*pdtphys
            end do

           end if 
          end do
 
        end do
        end do
        end do
        end do

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

   
      

     
#ifdef CPP_XIOS      
       do it=1,nqtot
       CALL  send_xios_field("flux_"//tname(it),
     &                     flux(:,it)) 
      end do 
#endif



      RETURN
      END