!
! $Header: /home/cvsroot/LMDZ4/libf/phylmd/phytrac.F,v 1.16 2006/03/24 15:06:23 lmdzadmin Exp $
!
c
c
      SUBROUTINE phytrac (nstep,
     I                    gmtime,
     I                    debutphy,
     I                    lafin,
     I                    nqmax,
     I                    nlon,
     I                    nlev, 
     I                    pdtphys,
     I                    u,
     I                    v,
     I                    t_seri,
     I                    paprs,
     I                    pplay,
     I                    coefh,
     I                    yu1,
     I                    yv1,
     I                    ftsol,
     I                    xlat,
     I                    xlon,
     I                    zlev, 
     I                    presnivs,
     I                    pphis,
     I                    pphi,
     I                    albsol,
     O                    tr_seri)

      USE ioipsl

      IMPLICIT none
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======================================================================
#include "YOMCST.h"
#include "dimensions.h"
#include "dimphy.h"
#include "clesphys.h" !///utile?
#include "temps.h"
#include "paramet.h"
#include "control.h"
#include "comgeomphy.h"
#include "advtrac.h"
c======================================================================

c Arguments:
c
c   EN ENTREE:
c   ==========
c
c   divers:
c   -------
c
      integer nlon  ! nombre de points horizontaux
      integer nlev  ! nombre de couches verticales
      integer nqmax ! nombre de traceurs auxquels on applique la physique
      integer nstep  ! appel physique
      integer nseuil ! numero du premier traceur non CV
c      integer julien !jour julien
c      integer itop_con(nlon)
c      integer ibas_con(nlon)
      real gmtime
      real pdtphys  ! pas d'integration pour la physique (seconde)
      real t_seri(nlon,nlev) ! temperature
      real tr_seri(nlon,nlev,nqmax) ! traceur  
      real u(nlon,nlev)
      real v(nlon,nlev)
      real albsol(nlon)  ! albedo surface
      real paprs(nlon,nlev+1)  ! pression pour chaque inter-couche (en Pa)
      real ps(nlon)  ! pression surface
      real pplay(nlon,nlev)  ! pression pour le mileu de chaque couche (en Pa)
      real pphi(nlon,nlev) ! geopotentiel
      real pphis(nlon)
      REAL presnivs(nlev)
      logical debutphy       ! le flag de l'initialisation de la physique
      logical lafin          ! le flag de la fin de la physique
c      REAL flxmass_w(nlon,nlev)

cAA Rem : nqmax : nombre de vrais traceurs est defini dans dimphy.h
c
c   convection:
c   -----------
c
      REAL pmfu(nlon,nlev)  ! flux de masse dans le panache montant
      REAL pmfd(nlon,nlev)  ! flux de masse dans le panache descendant
      REAL pen_u(nlon,nlev) ! flux entraine dans le panache montant

c
c   Couche limite:
c   --------------
c
      REAL coefh(nlon,nlev) ! coeff melange CL
      REAL yu1(nlon)        ! vents au premier niveau
      REAL yv1(nlon)        ! vents au premier niveau
      REAL xlat(nlon)       ! latitudes pour chaque point 
      REAL xlon(nlon)       ! longitudes pour chaque point 
      REAL zlev(nlon,nlev+1) ! altitude a chaque niveau (interface inferieure de la couche) 
cAA
cAA Arguments necessaires pour les sources et puits de traceur:
cAA ----------------
cAA
      real ftsol(nlon)  ! Temperature du sol (surf)(Kelvin)

cAA ----------------------------
cAA  VARIABLES LOCALES TRACEURS
cAA ----------------------------
cAA
cAA Sources et puits des traceurs:
cAA ------------------------------
cAA
      REAL source(klon)       ! a voir lorsque le flux est prescrit 

      CHARACTER*2 itn
C maf ioipsl
      CHARACTER*2 str2
      INTEGER nhori, nvert
      REAL zsto, zout, zjulian
      INTEGER nid_tra
      SAVE nid_tra
      INTEGER nid_tra2,nid_tra3
      SAVE nid_tra2,nid_tra3
      INTEGER ndex(1)
      INTEGER ndex2d(iim*(jjm+1)),ndex3d(iim*(jjm+1)*klev)
      REAL zx_tmp_2d(iim,jjm+1), zx_tmp_3d(iim,jjm+1,klev)
      REAL zx_lon(iim,jjm+1), zx_lat(iim,jjm+1)
c
      integer itau_w   ! pas de temps ecriture = nstep + itau_phy
c

C
C Variables liees a l'ecriture de la bande histoire : phytrac.nc
c
c      INTEGER ecrit_tra
c      SAVE ecrit_tra   
      logical ok_sync
      parameter (ok_sync = .true.)
C
C nature du traceur
c
      logical aerosol(nqmx)  ! Nature du traceur
c                            ! aerosol(it) = true  => aerosol 
c                            ! aerosol(it) = false => gaz 
c                            ! nat_trac(it) = 1. aerosol
      logical clsol(nqmx)    ! clsol(it) = true => CL sol calculee
      save aerosol,clsol
C
C les traceurs
C
c===================
c it--------indice de traceur
c k,i---------indices long, vert
c===================
c Variables deja declarees dont on a besoin pour traceurs   
c   k,i,it,tr_seri(klon,klev,nqmax),pplay(nlon,nlev), 
      real zprof(klon,klev)
c      real pzero,gamma
c      parameter (pzero=85000.)
c      parameter (gamma=5000.)
      real deltatr(klon,klev,nqmax) ! ecart au profil de ref zprof
      real tau             ! temps de relaxation vers le profil
c tau pourra dependre de it      
      parameter (tau=1.) ! duree de vie du compose en secondes 3.e10
c======================================================================
c
c Declaration des procedures appelees
c
c--modif convection tiedtke
      INTEGER i, k, it
      INTEGER iq, iiq
      REAL delp(klon,klev)
c--end modif
c
c Variables liees a l'ecriture de la bande histoire physique
c
c Variables locales pour effectuer les appels en serie
c----------------------------------------------------
c
      REAL d_tr(klon,klev), d_trs(klon) ! tendances de traceurs 
      REAL d_tr_cl(klon,klev,nqmax) ! tendance de traceurs  couche limite
      REAL d_tr_cv(klon,klev,nqmax) ! tendance de traceurs  conv pour chq traceur
C
      character*80 abort_message
c
c   Controles
c-------------
      logical first,couchelimite,convection
      save first,couchelimite,convection
c Olivia
       data first,couchelimite,convection
     s     /.true.,.false.,.false./

c======================================================================
         ps(:)=paprs(:,1)
c---------
c debutphy
c---------
         if (debutphy) then
                 print*,"DEBUT PHYTRAC"
C         
c=============================================================
c   Initialisation des traceurs
c=============================================================
c
c Initialisation de la nature des traceurs
c
         DO it = 1, nqmax
            aerosol(it) = .FALSE.  ! Tous les traceurs sont des gaz par defaut
            clsol(it) = .FALSE.  ! Par defaut couche limite avec flux prescrit
         ENDDO
c
c===========
c definition de traceurs idealises
c==========
c
c I) Declaration directe du traceur a altitude fixee
c
c a) traceur en carre OK
c
c         do i=1,klon
c         tr_seri(i,:,1)=0.
c        if ((xlat(i)>=0.).and.(xlat(i)<=-30.)) then
c          if ((xlon(i)>=0.).and.(xlon(i)<=40.)) then
c              tr_seri(i,10,1)=1.
c          endif
c        endif
c      end do 
c
c b) traceur a une bande en latitudeOK
c
c        do i=1,klon
c         tr_seri(i,:,1)=0.
c        if ((xlat(i)>=50.).and.(xlat(i)<=60.)) then
c              tr_seri(i,25,1)=1.
c        endif
c      end do  
c
c c) traceur a plusieurs bandes en latitude OK
c
c         do i=1,klon
c        tr_seri(i,:,2)=0.
c        if ((xlat(i)>=50.).and.(xlat(i)<=70.)) then
c             tr_seri(i,10,2)=1.
c        endif
c        if ((xlat(i)>=-10.).and.(xlat(i)<=10.)) then
c              tr_seri(i,10,2)=1.
c        endif
c
c        if ((xlat(i)>=-70.).and.(xlat(i)<=-50.)) then
c              tr_seri(i,10,2)=1. 
c        endif         
c      end do 
c
c d) traceur a une bande en altitude OK
c
c       do k=1,klev+1
c         tr_seri(:,k,1)=0.
c         if ((pplay(klon/2,k)>=1.e5).and.(pplay(klon/2,k)<=1.e6)) then
c             tr_seri(:,k,1)=1.
c         endif
c       end do 
c
c dbis) plusieurs traceurs a une bande en altitude OK
c
       do k=1,klev
         tr_seri(:,k,1)=0.
         if ((pplay(klon/2,k)>=9.e5).and.(pplay(klon/2,k)<=9.e6)) then
             tr_seri(:,k,1)=1.
         endif
       end do 

       do k=1,klev
         tr_seri(:,k,2)=0.
         if ((pplay(klon/2,k)>=1.e4).and.(pplay(klon/2,k)<=1.e5)) then
             tr_seri(:,k,2)=1.
         endif
       end do 

c      do k=1,klev
c        tr_seri(:,k,3)=0.
c        if ((pplay(klon/2,k)>=5.e1).and.(pplay(klon/2,k)<=5.e2)) then
c            tr_seri(:,k,3)=1.
c        endif
c      end do 
c       
c e) plusieurs couches verticales de traceurs, a plusieurs bandes en latitude???
c       
c au sol
c         do i=1,klon
c        tr_seri(i,:,1)=0.
c        if ((xlat(i)>=50.).and.(xlat(i)<=70.)) then
c             tr_seri(i,5,1)=1.
c        endif
c        if ((xlat(i)>=-10.).and.(xlat(i)<=10.)) then
c              tr_seri(i,5,1)=1.
c        endif
c
c        if ((xlat(i)>=-70.).and.(xlat(i)<=-50.)) then
c              tr_seri(i,5,1)=1. 
c        endif         
c      end do 
c 
c         do i=1,klon
c        tr_seri(i,:,2)=0.
c        if ((xlat(i)>=50.).and.(xlat(i)<=70.)) then
c             tr_seri(i,10,2)=1.
c        endif
c        if ((xlat(i)>=-10.).and.(xlat(i)<=10.)) then
c              tr_seri(i,10,2)=1.
c      endif
c
c        if ((xlat(i)>=-70.).and.(xlat(i)<=-50.)) then
c              tr_seri(i,10,2)=1. 
c        endif         
c      end do 
c
c         do i=1,klon
c        tr_seri(i,:,3)=0.
c        if ((xlat(i)>=50.).and.(xlat(i)<=70.)) then
c             tr_seri(i,30,3)=1.
c        endif
c        if ((xlat(i)>=-10.).and.(xlat(i)<=10.)) then
c              tr_seri(i,30,3)=1.
c        endif
c
c        if ((xlat(i)>=-70.).and.(xlat(i)<=-50.)) then
c              tr_seri(i,30,3)=1. 
c        endif         
c      end do 
c        
c        do i=1,klon
c        tr_seri(i,:,4)=0.
c        if ((xlat(i)>=50.).and.(xlat(i)<=70.)) then
c             tr_seri(i,45,4)=1.
c        endif
c        if ((xlat(i)>=-10.).and.(xlat(i)<=10.)) then
c              tr_seri(i,45,4)=1.
c        endif
c
c        if ((xlat(i)>=-70.).and.(xlat(i)<=-50.)) then
c              tr_seri(i,45,4)=1. 
c        endif         
c      end do 
c      
c II) Declaration d'un profil vertical de traceur OK
c 
c      do i=1,klon
c       do k=1,klev 
c         zprof(i,k)=0.5*(1.+tanh((pplay(i,k)-pzero)/gamma))
c--------x=alog10P+b
c        zprof(i,k)=-8.79e-6*(log10(pplay(i,k)))+7.6e-5
c        zprof(i,k)=8.79e-6*(log10(pplay(i,k)))-7.6e-5
c--------alog10x=alog10P+b        
c        zprof(i,k)=-0.31*(log10(pplay(i,k)))+5.2
c        zprof(i,k)=-0.31*(log10(pplay(i,k)))-2.68
c       enddo  
c      enddo
c
c      do i=1,klon
c       do k=1,klev 
c         zprof(i,k)=0.5*(1.+tanh((pplay(i,k)-pzero)/gamma))
c       enddo 
c      enddo 
c-------------
c fin debutphy
c-------------
      ENDIF  ! fin debutphy 

c======================================================================
c Rappel vers un profil
c======================================================================
c         do it=1,nqmax-2
c           do k=1,klev
c            do i=1,klon
c             deltatr(i,k,1) = 0.
c             tr_seri(i,k,1) = 0.
c             deltatr(i,k,1) = (-tr_seri(i,k,1)+zprof(i,k))/tau
c             tr_seri(i,k,1) =  tr_seri(i,k,1) + deltatr(i,k,1)
c            enddo
c            enddo
c         enddo

c======================================================================
c   Calcul de l'effet de la couche limite
c======================================================================

      if (couchelimite) then

      DO k = 1, klev
      DO i = 1, klon
         delp(i,k) = paprs(i,k)-paprs(i,k+1)
      ENDDO
      ENDDO

C maf modif pour tenir compte du cas rnpb + traceur
      DO it=1, nqmax

c     print *,'it',it,clsol(it)
      if (clsol(it)) then  ! couche limite avec quantite dans le sol calculee

      else ! couche limite avec flux prescrit

Cmaf provisoire source / traceur a creer
        DO i=1, klon
          source(i) = 0.0 ! pas de source, pour l'instant
        ENDDO
C

          CALL cltrac(pdtphys, coefh,t_seri,
     s               tr_seri(1,1,it), source,
     e               paprs, pplay, delp,
     s               d_tr_cl(1,1,it))
            DO k = 1, nlev
               DO i = 1, klon
                  tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_cl(i,k,it)
               ENDDO
            ENDDO
Cmaf          WRITE(itn,'(i1)') it
cmaf          CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'cltracn it='//itn)
      endif
      ENDDO
c
      endif ! couche limite

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

      RETURN
      END