Index: /LMDZ4/branches/LMDZ4-dev/libf/dyn3d/infotrac.F90
===================================================================
--- /LMDZ4/branches/LMDZ4-dev/libf/dyn3d/infotrac.F90	(revision 1190)
+++ /LMDZ4/branches/LMDZ4-dev/libf/dyn3d/infotrac.F90	(revision 1191)
@@ -21,8 +21,11 @@
   INTEGER, ALLOCATABLE, DIMENSION(:), SAVE    :: niadv ! equivalent dyn / physique
 
-! Variables for INCA
+! conv_flg(it)=0 : convection desactivated for tracer number it 
   INTEGER, ALLOCATABLE, DIMENSION(:), SAVE  :: conv_flg
+! pbl_flg(it)=0  : boundary layer diffusion desactivaded for tracer number it 
   INTEGER, ALLOCATABLE, DIMENSION(:), SAVE  :: pbl_flg
 
+  CHARACTER(len=4),SAVE :: type_trac
+ 
 CONTAINS
 
@@ -80,4 +83,11 @@
     descrq(20)='SLP'
     descrq(30)='PRA'
+    
+
+    IF (config_inca=='none') THEN
+       type_trac='lmdz'
+    ELSE
+       type_trac='inca'
+    END IF
 
 !-----------------------------------------------------------------------
@@ -87,5 +97,5 @@
 !
 !-----------------------------------------------------------------------
-    IF (config_inca == 'none') THEN
+    IF (type_trac == 'lmdz') THEN
        OPEN(90,file='traceur.def',form='formatted',status='old', iostat=ierr)
        IF(ierr.EQ.0) THEN
@@ -109,13 +119,11 @@
     END IF
 !
-! Allocate variables depending on nqtrue
+! Allocate variables depending on nqtrue and nbtr
 !
     ALLOCATE(tnom_0(nqtrue), hadv(nqtrue), vadv(nqtrue))
-
-    IF (config_inca /= 'none') THEN
-       ! Varaibles only needed in case of INCA
-       ALLOCATE(conv_flg(nbtr), pbl_flg(nbtr), tracnam(nbtr))
-    END IF
-       
+    ALLOCATE(conv_flg(nbtr), pbl_flg(nbtr), tracnam(nbtr))
+    conv_flg(:) = 1 ! convection activated for all tracers
+    pbl_flg(:)  = 1 ! boundary layer activated for all tracers
+
 !-----------------------------------------------------------------------
 ! 2)     Choix  des schemas d'advection pour l'eau et les traceurs
@@ -144,5 +152,5 @@
 !    Get choice of advection schema from file tracer.def or from INCA
 !---------------------------------------------------------------------
-    IF (config_inca == 'none') THEN
+    IF (type_trac == 'lmdz') THEN
        IF(ierr.EQ.0) THEN
           ! Continue to read tracer.def
@@ -172,5 +180,5 @@
        END DO
 
-    ELSE  ! config_inca='aero' ou 'chem'
+    ELSE  ! type_trac=inca : config_inca='aero' ou 'chem'
 ! le module de chimie fournit les noms des traceurs
 ! et les schemas d'advection associes.
@@ -191,5 +199,5 @@
        END DO
 
-    END IF ! config_inca
+    END IF ! type_trac
 
 !-----------------------------------------------------------------------
@@ -319,5 +327,5 @@
 !
     DEALLOCATE(tnom_0, hadv, vadv)
-    IF (config_inca /= 'none') DEALLOCATE(tracnam)
+    DEALLOCATE(tracnam)
 
 999 FORMAT (i2,1x,i2,1x,a8)
Index: /LMDZ4/branches/LMDZ4-dev/libf/dyn3dpar/infotrac.F90
===================================================================
--- /LMDZ4/branches/LMDZ4-dev/libf/dyn3dpar/infotrac.F90	(revision 1190)
+++ /LMDZ4/branches/LMDZ4-dev/libf/dyn3dpar/infotrac.F90	(revision 1191)
@@ -21,8 +21,11 @@
   INTEGER, ALLOCATABLE, DIMENSION(:), SAVE    :: niadv ! equivalent dyn / physique
 
-! Variables for INCA
+! conv_flg(it)=0 : convection desactivated for tracer number it 
   INTEGER, ALLOCATABLE, DIMENSION(:), SAVE  :: conv_flg
+! pbl_flg(it)=0  : boundary layer diffusion desactivaded for tracer number it 
   INTEGER, ALLOCATABLE, DIMENSION(:), SAVE  :: pbl_flg
 
+  CHARACTER(len=4),SAVE :: type_trac
+ 
 CONTAINS
 
@@ -80,4 +83,11 @@
     descrq(20)='SLP'
     descrq(30)='PRA'
+    
+
+    IF (config_inca=='none') THEN
+       type_trac='lmdz'
+    ELSE
+       type_trac='inca'
+    END IF
 
 !-----------------------------------------------------------------------
@@ -87,5 +97,5 @@
 !
 !-----------------------------------------------------------------------
-    IF (config_inca == 'none') THEN
+    IF (type_trac == 'lmdz') THEN
        OPEN(90,file='traceur.def',form='formatted',status='old', iostat=ierr)
        IF(ierr.EQ.0) THEN
@@ -109,13 +119,11 @@
     END IF
 !
-! Allocate variables depending on nqtrue
+! Allocate variables depending on nqtrue and nbtr
 !
     ALLOCATE(tnom_0(nqtrue), hadv(nqtrue), vadv(nqtrue))
-
-    IF (config_inca /= 'none') THEN
-       ! Varaibles only needed in case of INCA
-       ALLOCATE(conv_flg(nbtr), pbl_flg(nbtr), tracnam(nbtr))
-    END IF
-       
+    ALLOCATE(conv_flg(nbtr), pbl_flg(nbtr), tracnam(nbtr))
+    conv_flg(:) = 1 ! convection activated for all tracers
+    pbl_flg(:)  = 1 ! boundary layer activated for all tracers
+
 !-----------------------------------------------------------------------
 ! 2)     Choix  des schemas d'advection pour l'eau et les traceurs
@@ -144,5 +152,5 @@
 !    Get choice of advection schema from file tracer.def or from INCA
 !---------------------------------------------------------------------
-    IF (config_inca == 'none') THEN
+    IF (type_trac == 'lmdz') THEN
        IF(ierr.EQ.0) THEN
           ! Continue to read tracer.def
@@ -172,5 +180,5 @@
        END DO
 
-    ELSE  ! config_inca='aero' ou 'chem'
+    ELSE  ! type_trac=inca : config_inca='aero' ou 'chem'
 ! le module de chimie fournit les noms des traceurs
 ! et les schemas d'advection associes.
@@ -191,5 +199,5 @@
        END DO
 
-    END IF ! config_inca
+    END IF ! type_trac
 
 !-----------------------------------------------------------------------
@@ -319,5 +327,5 @@
 !
     DEALLOCATE(tnom_0, hadv, vadv)
-    IF (config_inca /= 'none') DEALLOCATE(tracnam)
+    DEALLOCATE(tracnam)
 
 999 FORMAT (i2,1x,i2,1x,a8)
Index: /LMDZ4/branches/LMDZ4-dev/libf/phylmd/YOECUMF.h
===================================================================
--- /LMDZ4/branches/LMDZ4-dev/libf/phylmd/YOECUMF.h	(revision 1190)
+++ /LMDZ4/branches/LMDZ4-dev/libf/phylmd/YOECUMF.h	(revision 1191)
@@ -1,43 +1,48 @@
 !
-! $Header$
+! $Id$
 !
-C     ----------------------------------------------------------------
-C*    *COMMON* *YOECUMF* - PARAMETERS FOR CUMULUS MASSFLUX SCHEME
-C     ----------------------------------------------------------------
-C
-      COMMON /YOECUMF/
-     L                 LMFPEN,LMFSCV,LMFMID,LMFDD,LMFDUDV,
-     R                 ENTRPEN,ENTRSCV,ENTRMID,ENTRDD,CMFCTOP,
-     R                 CMFCMAX,CMFCMIN,CMFDEPS,RHCDD,CPRCON
-C
+!  ATTENTION!!!!: ce fichier include est compatible format fixe/format libre
+!                 veillez n'utiliser que des ! pour les commentaires
+!                 et bien positionner les & des lignes de continuation
+!                 (les placer en colonne 6 et en colonne 73)
+!
+!     ----------------------------------------------------------------
+!*    *COMMON* *YOECUMF* - PARAMETERS FOR CUMULUS MASSFLUX SCHEME
+!     ----------------------------------------------------------------
+!
+      COMMON /YOECUMF/                                                  &
+     &                 LMFPEN,LMFSCV,LMFMID,LMFDD,LMFDUDV,              &
+     &                 ENTRPEN,ENTRSCV,ENTRMID,ENTRDD,CMFCTOP,          &
+     &                 CMFCMAX,CMFCMIN,CMFDEPS,RHCDD,CPRCON
+
       LOGICAL          LMFPEN,LMFSCV,LMFMID,LMFDD,LMFDUDV
       REAL ENTRPEN, ENTRSCV, ENTRMID, ENTRDD
       REAL CMFCTOP, CMFCMAX, CMFCMIN, CMFDEPS, RHCDD, CPRCON
-c$OMP THREADPRIVATE(/YOECUMF/)
-C
-*if (DOC,declared) <> 'UNKNOWN'
-C*    *COMMON* *YOECUMF* - PARAMETERS FOR CUMULUS MASSFLUX SCHEME
-C
-C     M.TIEDTKE       E. C. M. W. F.      18/1/89
-C
-C     NAME      TYPE      PURPOSE
-C     ----      ----      -------
-C
-C     LMFPEN    LOGICAL  TRUE IF PENETRATIVE CONVECTION IS SWITCHED ON
-C     LMFSCV    LOGICAL  TRUE IF SHALLOW     CONVECTION IS SWITCHED ON
-C     LMFMID    LOGICAL  TRUE IF MIDLEVEL    CONVECTION IS SWITCHED ON
-C     LMFDD     LOGICAL  TRUE IF CUMULUS DOWNDRAFT      IS SWITCHED ON
-C     LMFDUDV   LOGICAL  TRUE IF CUMULUS FRICTION       IS SWITCHED ON
-C     ENTRPEN   REAL     ENTRAINMENT RATE FOR PENETRATIVE CONVECTION
-C     ENTRSCV   REAL     ENTRAINMENT RATE FOR SHALLOW CONVECTION
-C     ENTRMID   REAL     ENTRAINMENT RATE FOR MIDLEVEL CONVECTION
-C     ENTRDD    REAL     ENTRAINMENT RATE FOR CUMULUS DOWNDRAFTS
-C     CMFCTOP   REAL     RELAT. CLOUD MASSFLUX AT LEVEL ABOVE NONBUOYANC
-C     CMFCMAX   REAL     MAXIMUM MASSFLUX VALUE ALLOWED FOR
-C     CMFCMIN   REAL     MINIMUM MASSFLUX VALUE (FOR SAFETY)
-C     CMFDEPS   REAL     FRACTIONAL MASSFLUX FOR DOWNDRAFTS AT LFS
-C     RHCDD     REAL     RELATIVE SATURATION IN DOWNDRAFTS
-C     CPRCON    REAL     COEFFICIENTS FOR DETERMINING CONVERSION
-C                        FROM CLOUD WATER TO RAIN
-*ifend
-C     ----------------------------------------------------------------
+!$OMP THREADPRIVATE(/YOECUMF/)
+!
+!*if (DOC,declared) <> 'UNKNOWN'
+!*    *COMMON* *YOECUMF* - PARAMETERS FOR CUMULUS MASSFLUX SCHEME
+!
+!     M.TIEDTKE       E. C. M. W. F.      18/1/89
+!
+!     NAME      TYPE      PURPOSE
+!     ----      ----      -------
+!
+!     LMFPEN    LOGICAL  TRUE IF PENETRATIVE CONVECTION IS SWITCHED ON
+!     LMFSCV    LOGICAL  TRUE IF SHALLOW     CONVECTION IS SWITCHED ON
+!     LMFMID    LOGICAL  TRUE IF MIDLEVEL    CONVECTION IS SWITCHED ON
+!     LMFDD     LOGICAL  TRUE IF CUMULUS DOWNDRAFT      IS SWITCHED ON
+!     LMFDUDV   LOGICAL  TRUE IF CUMULUS FRICTION       IS SWITCHED ON
+!     ENTRPEN   REAL     ENTRAINMENT RATE FOR PENETRATIVE CONVECTION
+!     ENTRSCV   REAL     ENTRAINMENT RATE FOR SHALLOW CONVECTION
+!     ENTRMID   REAL     ENTRAINMENT RATE FOR MIDLEVEL CONVECTION
+!     ENTRDD    REAL     ENTRAINMENT RATE FOR CUMULUS DOWNDRAFTS
+!     CMFCTOP   REAL     RELAT. CLOUD MASSFLUX AT LEVEL ABOVE NONBUOYANC
+!     CMFCMAX   REAL     MAXIMUM MASSFLUX VALUE ALLOWED FOR
+!     CMFCMIN   REAL     MINIMUM MASSFLUX VALUE (FOR SAFETY)
+!     CMFDEPS   REAL     FRACTIONAL MASSFLUX FOR DOWNDRAFTS AT LFS
+!     RHCDD     REAL     RELATIVE SATURATION IN DOWNDRAFTS
+!     CPRCON    REAL     COEFFICIENTS FOR DETERMINING CONVERSION
+!                        FROM CLOUD WATER TO RAIN
+!*ifend
+!     ----------------------------------------------------------------
Index: DZ4/branches/LMDZ4-dev/libf/phylmd/cltrac.F
===================================================================
--- /LMDZ4/branches/LMDZ4-dev/libf/phylmd/cltrac.F	(revision 1190)
+++ 	(revision )
@@ -1,122 +1,0 @@
-!
-! $Header$
-!
-      SUBROUTINE cltrac(dtime,coef,t,tr,flux,paprs,pplay,delp,
-     s                  d_tr)
-      USE dimphy
-      IMPLICIT none
-c======================================================================
-c Auteur(s): O. Boucher (LOA/LMD) date: 19961127
-c            inspire de clvent
-c Objet: diffusion verticale de traceurs avec flux fixe a la surface
-c        ou/et flux du type c-drag
-c======================================================================
-c Arguments:
-c dtime----input-R- intervalle du temps (en second)
-c coef-----input-R- le coefficient d'echange (m**2/s) l>1
-c t--------input-R- temperature (K)
-c tr-------input-R- la q. de traceurs
-c flux-----input-R- le flux de traceurs a la surface
-c paprs----input-R- pression a inter-couche (Pa)
-c pplay----input-R- pression au milieu de couche (Pa)
-c delp-----input-R- epaisseur de couche (Pa)
-c cdrag----input-R- cdrag pour le flux de surface (non active)
-c tr0------input-R- traceurs a la surface ou dans l'ocean (non active)
-c d_tr-----output-R- le changement de tr
-c flux_tr--output-R- flux de tr
-c======================================================================
-cym#include "dimensions.h"
-cym#include "dimphy.h"
-      REAL dtime
-      REAL coef(klon,klev)
-      REAL t(klon,klev), tr(klon,klev)
-      REAL paprs(klon,klev+1), pplay(klon,klev), delp(klon,klev)
-      REAL d_tr(klon,klev)
-      REAL flux(klon), cdrag(klon), tr0(klon)
-c      REAL flux_tr(klon,klev)
-c======================================================================
-#include "YOMCST.h"
-c======================================================================
-      INTEGER i, k
-      REAL zx_ctr(klon,2:klev)
-      REAL zx_dtr(klon,2:klev)
-      REAL zx_buf(klon)
-      REAL zx_coef(klon,klev)
-      REAL local_tr(klon,klev)
-      REAL zx_alf1(klon), zx_alf2(klon), zx_flux(klon)
-c======================================================================
-      DO k = 1, klev
-      DO i = 1, klon
-         local_tr(i,k) = tr(i,k)
-      ENDDO
-      ENDDO
-c
-
-c======================================================================
-      DO i = 1, klon
-         zx_alf1(i) = (paprs(i,1)-pplay(i,2))/(pplay(i,1)-pplay(i,2))
-         zx_alf2(i) = 1.0 - zx_alf1(i)
-         zx_flux(i) =  -flux(i)*dtime*RG
-c--pour le moment le flux est prescrit
-c--cdrag et zx_coef(1) vaut 0
-         cdrag(i) = 0.0 
-         tr0(i) = 0.0
-         zx_coef(i,1) = cdrag(i)*dtime*RG 
-      ENDDO
-c======================================================================
-      DO k = 2, klev
-      DO i = 1, klon
-         zx_coef(i,k) = coef(i,k)*RG/(pplay(i,k-1)-pplay(i,k))
-     .                  *(paprs(i,k)*2/(t(i,k)+t(i,k-1))/RD)**2
-         zx_coef(i,k) = zx_coef(i,k)*dtime*RG
-      ENDDO
-      ENDDO
-c======================================================================
-      DO i = 1, klon
-         zx_buf(i) = delp(i,1) + zx_coef(i,1)*zx_alf1(i) + zx_coef(i,2)
-         zx_ctr(i,2) = (local_tr(i,1)*delp(i,1)+
-     .                  zx_coef(i,1)*tr0(i)-zx_flux(i))/zx_buf(i)
-         zx_dtr(i,2) = (zx_coef(i,2)-zx_alf2(i)*zx_coef(i,1)) / 
-     .                  zx_buf(i)
-      ENDDO
-c
-      DO k = 3, klev
-      DO i = 1, klon
-         zx_buf(i) = delp(i,k-1) + zx_coef(i,k)
-     .                  + zx_coef(i,k-1)*(1.-zx_dtr(i,k-1))
-         zx_ctr(i,k) = (local_tr(i,k-1)*delp(i,k-1)
-     .                  +zx_coef(i,k-1)*zx_ctr(i,k-1) )/zx_buf(i)
-         zx_dtr(i,k) = zx_coef(i,k)/zx_buf(i)
-      ENDDO
-      ENDDO
-      DO i = 1, klon
-         local_tr(i,klev) = ( local_tr(i,klev)*delp(i,klev)
-     .                        +zx_coef(i,klev)*zx_ctr(i,klev) )
-     .                   / ( delp(i,klev) + zx_coef(i,klev)
-     .                       -zx_coef(i,klev)*zx_dtr(i,klev) )
-      ENDDO
-      DO k = klev-1, 1, -1
-      DO i = 1, klon
-         local_tr(i,k) = zx_ctr(i,k+1) + zx_dtr(i,k+1)*local_tr(i,k+1)
-      ENDDO
-      ENDDO
-c======================================================================
-c== flux_tr est le flux de traceur (positif vers bas)
-c      DO i = 1, klon
-c         flux_tr(i,1) = zx_coef(i,1)/(RG*dtime)
-c      ENDDO
-c      DO k = 2, klev
-c      DO i = 1, klon
-c         flux_tr(i,k) = zx_coef(i,k)/(RG*dtime)
-c     .               * (local_tr(i,k)-local_tr(i,k-1))
-c      ENDDO
-c      ENDDO
-c======================================================================
-      DO k = 1, klev
-      DO i = 1, klon
-         d_tr(i,k) = local_tr(i,k) - tr(i,k)
-      ENDDO
-      ENDDO
-c
-      RETURN
-      END
Index: /LMDZ4/branches/LMDZ4-dev/libf/phylmd/cltrac.F90
===================================================================
--- /LMDZ4/branches/LMDZ4-dev/libf/phylmd/cltrac.F90	(revision 1191)
+++ /LMDZ4/branches/LMDZ4-dev/libf/phylmd/cltrac.F90	(revision 1191)
@@ -0,0 +1,140 @@
+!
+! $Id $
+!
+SUBROUTINE cltrac(dtime,coef,t,tr,flux,paprs,pplay,delp,d_tr)
+  USE dimphy
+  IMPLICIT NONE
+!======================================================================
+! Auteur(s): O. Boucher (LOA/LMD) date: 19961127
+!            inspire de clvent
+! Objet: diffusion verticale de traceurs avec flux fixe a la surface
+!        ou/et flux du type c-drag
+!
+! Arguments:
+!-----------
+! dtime....input-R- intervalle du temps (en secondes)
+! coef.....input-R- le coefficient d'echange (m**2/s) l>1
+! t........input-R- temperature (K)
+! tr.......input-R- la q. de traceurs
+! flux.....input-R- le flux de traceurs a la surface
+! paprs....input-R- pression a inter-couche (Pa)
+! pplay....input-R- pression au milieu de couche (Pa)
+! delp.....input-R- epaisseur de couche (Pa)
+! cdrag....input-R- cdrag pour le flux de surface (non active)
+! tr0......input-R- traceurs a la surface ou dans l'ocean (non active)
+! d_tr.....output-R- le changement de tr
+! flux_tr..output-R- flux de tr
+!======================================================================
+  include "YOMCST.h"
+!
+! Entree
+! 
+  REAL,INTENT(IN)                        :: dtime
+  REAL,DIMENSION(klon,klev),INTENT(IN)   :: coef
+  REAL,DIMENSION(klon,klev),INTENT(IN)   :: t, tr
+  REAL,DIMENSION(klon),INTENT(IN)        :: flux !(at/s/m2)
+  REAL,DIMENSION(klon,klev+1),INTENT(IN) :: paprs 
+  REAL,DIMENSION(klon,klev),INTENT(IN)   :: pplay, delp
+!
+! Sorties
+!
+  REAL ,DIMENSION(klon,klev),INTENT(OUT) :: d_tr
+!  REAL ,DIMENSION(klon,klev),INTENT(OUT) :: flux_tr
+!
+! Local
+! 
+  INTEGER                   :: i, k
+  REAL,DIMENSION(klon)      :: cdrag, tr0
+  REAL,DIMENSION(klon,klev) :: zx_ctr
+  REAL,DIMENSION(klon,klev) :: zx_dtr
+  REAL,DIMENSION(klon)      :: zx_buf
+  REAL,DIMENSION(klon,klev) :: zx_coef
+  REAL,DIMENSION(klon,klev) :: local_tr
+  REAL,DIMENSION(klon)      :: zx_alf1,zx_alf2,zx_flux
+
+!======================================================================
+
+  DO k = 1, klev
+     DO i = 1, klon
+        local_tr(i,k) = tr(i,k)
+     ENDDO
+  ENDDO
+
+!======================================================================
+
+  DO i = 1, klon
+     zx_alf1(i) = (paprs(i,1)-pplay(i,2))/(pplay(i,1)-pplay(i,2))
+     zx_alf2(i) = 1.0 - zx_alf1(i)
+     zx_flux(i) =  -flux(i)*dtime*RG
+! Pour le moment le flux est prescrit cdrag et zx_coef(1) vaut 0
+     cdrag(i) = 0.0 
+     tr0(i) = 0.0
+     zx_coef(i,1) = cdrag(i)*dtime*RG 
+     zx_ctr(i,1)=0.
+     zx_dtr(i,1)=0.
+  ENDDO
+
+!======================================================================
+
+  DO k = 2, klev
+     DO i = 1, klon
+        zx_coef(i,k) = coef(i,k)*RG/(pplay(i,k-1)-pplay(i,k))   &
+             *(paprs(i,k)*2/(t(i,k)+t(i,k-1))/RD)**2
+        zx_coef(i,k) = zx_coef(i,k)*dtime*RG  
+     ENDDO
+  ENDDO
+
+!======================================================================
+
+  DO i = 1, klon
+     zx_buf(i) = delp(i,1) + zx_coef(i,1)*zx_alf1(i) + zx_coef(i,2)
+     !
+     zx_ctr(i,2) = (local_tr(i,1)*delp(i,1)+                  &
+          zx_coef(i,1)*tr0(i)-zx_flux(i))/zx_buf(i)
+     !
+     zx_dtr(i,2) = (zx_coef(i,2)-zx_alf2(i)*zx_coef(i,1)) /   & 
+          zx_buf(i)
+  ENDDO
+
+  DO k = 3, klev
+     DO i = 1, klon
+        zx_buf(i) = delp(i,k-1) + zx_coef(i,k)      &
+             + zx_coef(i,k-1)*(1.-zx_dtr(i,k-1))
+        zx_ctr(i,k) = (local_tr(i,k-1)*delp(i,k-1)  & 
+             +zx_coef(i,k-1)*zx_ctr(i,k-1) )/zx_buf(i)
+        zx_dtr(i,k) = zx_coef(i,k)/zx_buf(i)
+     ENDDO
+  ENDDO
+
+  DO i = 1, klon
+     local_tr(i,klev) = ( local_tr(i,klev)*delp(i,klev) &
+          +zx_coef(i,klev)*zx_ctr(i,klev) )             &
+          / ( delp(i,klev) + zx_coef(i,klev)            &
+          -zx_coef(i,klev)*zx_dtr(i,klev) )
+  ENDDO
+
+  DO k = klev-1, 1, -1
+     DO i = 1, klon
+        local_tr(i,k) = zx_ctr(i,k+1) + zx_dtr(i,k+1)*local_tr(i,k+1)
+     ENDDO
+  ENDDO
+
+!======================================================================
+!== flux_tr est le flux de traceur (positif vers bas)
+!      DO i = 1, klon
+!         flux_tr(i,1) = zx_coef(i,1)/(RG*dtime)
+!      ENDDO
+!      DO k = 2, klev
+!      DO i = 1, klon
+!         flux_tr(i,k) = zx_coef(i,k)/(RG*dtime)
+!     .               * (local_tr(i,k)-local_tr(i,k-1))
+!      ENDDO
+!      ENDDO
+!======================================================================
+  DO k = 1, klev
+     DO i = 1, klon
+        d_tr(i,k) = local_tr(i,k) - tr(i,k)
+     ENDDO
+  ENDDO
+  
+END SUBROUTINE cltrac
Index: DZ4/branches/LMDZ4-dev/libf/phylmd/cltracrn.F
===================================================================
--- /LMDZ4/branches/LMDZ4-dev/libf/phylmd/cltracrn.F	(revision 1190)
+++ 	(revision )
@@ -1,294 +1,0 @@
-!
-      SUBROUTINE cltracrn( itr, dtime,u1lay, v1lay,
-     e              cdrag,coef,t,ftsol,pctsrf,
-     e              tr,trs,paprs,pplay,delp,
-     e              masktr,fshtr,hsoltr,tautr,vdeptr,
-     e              lat,
-     s              d_tr,d_trs )
-
-      USE dimphy
-      IMPLICIT none
-c======================================================================
-c Auteur(s): Alex/LMD) date:  fev 99
-c            inspire de clqh + clvent
-c Objet: diffusion verticale de traceurs avec quantite de traceur ds 
-c        le sol ( reservoir de sol de radon ) 
-c        
-c note : pour l'instant le traceur dans le sol et le flux sont
-c        calcules mais ils ne servent que de diagnostiques
-c        seule la tendance sur le traceur est sortie (d_tr)
-c======================================================================
-c Arguments:
-c itr---  -input-R- le type de traceur 1- Rn 2 - Pb
-c dtime----input-R- intervalle du temps (en second)
-c u1lay----input-R- vent u de la premiere couche (m/s)
-c v1lay----input-R- vent v de la premiere couche (m/s)
-c cdrag----input-R- cdrag
-c coef-----input-R- le coefficient d'echange (m**2/s) l>1, valable uniquement pour k entre 2 et klev
-c t--------input-R- temperature (K)
-c paprs----input-R- pression a inter-couche (Pa)
-c pplay----input-R- pression au milieu de couche (Pa)
-c delp-----input-R- epaisseur de couche (Pa)
-c ftsol----input-R- temperature du sol (en Kelvin)
-c tr-------input-R- traceurs
-c trs------input-R- traceurs dans le sol
-c masktr---input-R- Masque reservoir de sol traceur (1 = reservoir)
-c fshtr----input-R- Flux surfacique de production dans le sol
-c tautr----input-R- Constante de decroissance du traceur
-c vdeptr---input-R- Vitesse de depot sec dans la couche brownienne
-c hsoltr---input-R- Epaisseur equivalente du reservoir de sol
-c lat-----input-R- latitude en degree
-c d_tr-----output-R- le changement de "tr"
-c d_trs----output-R- le changement de "trs"
-c======================================================================
-cym#include "dimensions.h"
-cym#include "dimphy.h"
-#include "YOMCST.h"
-#include "indicesol.h"
-c======================================================================
-      REAL dtime
-      REAL u1lay(klon), v1lay(klon)
-      REAL cdrag(klon)
-      REAL coef(klon,klev)
-      REAL t(klon,klev), ftsol(klon,nbsrf), pctsrf(klon,nbsrf) 
-      REAL tr(klon,klev), trs(klon)
-      REAL paprs(klon,klev+1), pplay(klon,klev), delp(klon,klev)
-      REAL masktr(klon) 
-      REAL fshtr(klon) 
-      REAL hsoltr
-      REAL tautr
-      REAL vdeptr
-      REAL lat(klon)   
-      REAL d_tr(klon,klev)
-c======================================================================
-      REAL flux_tr(klon,klev)  ! (diagnostic) flux de traceur
-      REAL d_trs(klon)         ! (diagnostic) traceur ds le sol
-c======================================================================
-      INTEGER i, k, itr, n, l
-      REAL rotrhi(klon)
-      REAL zx_coef(klon,klev)
-      REAL zx_buf(klon)
-      REAL zx_ctr(klon,klev)
-      REAL zx_dtr(klon,klev)
-      REAL zx_trs(klon)
-      REAL zx_a, zx_b
-
-      REAL local_tr(klon,klev)
-      REAL local_trs(klon)
-      REAL zts(klon)
-      REAL zx_alpha1(klon), zx_alpha2(klon)
-c======================================================================
-cAA Pour l'instant les 4 types de surface ne sont pas pris en compte
-cAA On fabrique avec zts un champ de temperature de sol  
-cAA que le pondere par la fraction de nature de sol.
-c 
-c      print*,'PASSAGE DANS CLTRACRN'
-
-      DO i = 1,klon
-            zts(i) = 0. 
-      ENDDO
-c
-      DO n=1,nbsrf
-         DO i = 1,klon
-            zts(i) = zts(i) + ftsol(i,n)*pctsrf(i,n)
-         ENDDO
-      ENDDO
-c
-      DO i = 1,klon
-          rotrhi(i) = RD * zts(i) / hsoltr 
-      END DO
-c
-      DO k = 1, klev
-      DO i = 1, klon
-         local_tr(i,k) = tr(i,k)
-      ENDDO
-      ENDDO
-c
-      DO i = 1, klon
-         local_trs(i) = trs(i)
-      ENDDO
-c======================================================================
-cAA   Attention si dans clmain zx_alf1(i) = 1.0 
-cAA   Il doit y avoir coherence (dc la meme chose ici)
-
-      DO i = 1, klon
-cAA         zx_alpha1(i) = (paprs(i,1)-pplay(i,2))/(pplay(i,1)-pplay(i,2))
-         zx_alpha1(i) = 1.0
-         zx_alpha2(i) = 1.0 - zx_alpha1(i)
-      ENDDO
-c======================================================================
-      DO i = 1, klon
-         zx_coef(i,1) = cdrag(i)
-     .                 * (1.0+SQRT(u1lay(i)**2+v1lay(i)**2))
-     .                 * pplay(i,1)/(RD*t(i,1))
-         zx_coef(i,1) = zx_coef(i,1) * dtime*RG
-      ENDDO
-c
-      DO k = 2, klev
-      DO i = 1, klon
-         zx_coef(i,k) = coef(i,k)*RG/(pplay(i,k-1)-pplay(i,k))
-     .                  *(paprs(i,k)*2/(t(i,k)+t(i,k-1))/RD)**2
-         zx_coef(i,k) = zx_coef(i,k) * dtime*RG
-      ENDDO
-      ENDDO
-c======================================================================
-      DO i = 1, klon
-         zx_buf(i)      = delp(i,klev) + zx_coef(i,klev)
-         zx_ctr(i,klev) = local_tr(i,klev)*delp(i,klev)/zx_buf(i)
-         zx_dtr(i,klev) = zx_coef(i,klev) / zx_buf(i)
-      ENDDO
-c
-      DO l = klev-1, 2 , -1
-      DO i = 1, klon
-         zx_buf(i) = delp(i,l)+zx_coef(i,l)
-     .              +zx_coef(i,l+1)*(1.-zx_dtr(i,l+1))
-         zx_ctr(i,l) = ( local_tr(i,l)*delp(i,l)
-     .                   + zx_coef(i,l+1)*zx_ctr(i,l+1) )/zx_buf(i)
-         zx_dtr(i,l) = zx_coef(i,l) / zx_buf(i)
-      ENDDO
-      ENDDO
-c
-      DO i = 1, klon
-         zx_buf(i) = delp(i,1) + zx_coef(i,2)*(1.-zx_dtr(i,2))
-     .               + masktr(i) * zx_coef(i,1)
-     .               *( zx_alpha1(i)+zx_alpha2(i)*zx_dtr(i,2) )
-         zx_ctr(i,1) = ( local_tr(i,1)*delp(i,1)
-     .                  + zx_ctr(i,2)
-     .                  *(zx_coef(i,2)
-     .                    - masktr(i) * zx_coef(i,1)
-     .                    *zx_alpha2(i) ) ) / zx_buf(i)
-         zx_dtr(i,1) = masktr(i) * zx_coef(i,1) / zx_buf(i)
-      ENDDO
-c======================================================================
-c Calculer d'abord local_trs nouvelle quantite dans le reservoir
-c de sol
-c
-c-------------------------
-c Au dessus des continents
-c-------------------------
-c Le pb peut se deposer partout : vdeptr = 10-3 m/s
-c Le Rn est traiter commme une couche Brownienne puisque vdeptr = 0.
-c
-      DO i = 1, klon
-c
-        IF ( NINT(masktr(i)) .EQ. 1 ) THEN
-         zx_trs(i) = local_trs(i)
-         zx_a = zx_trs(i)
-     .         +fshtr(i)*dtime*rotrhi(i)
-     .         +rotrhi(i)*masktr(i)*zx_coef(i,1)/RG
-     .         *(zx_ctr(i,1)*(zx_alpha1(i)+zx_alpha2(i)*zx_dtr(i,2))
-     .         +zx_alpha2(i)*zx_ctr(i,2))
-         zx_b = 1. + rotrhi(i)*masktr(i)*zx_coef(i,1)/RG
-     .              * (1.-zx_dtr(i,1)
-     .              *(zx_alpha1(i)+zx_alpha2(i)*zx_dtr(i,2)))
-     .              + dtime / tautr
-cAA: Pour l'instant, pour aller vite, le depot sec est traite
-C comme une decroissance
-     .              + dtime * vdeptr / hsoltr
-         zx_trs(i) = zx_a / zx_b
-         local_trs(i) = zx_trs(i)
-        ENDIF
-c
-c Si on est entre 60N et 70N on divise par 2 l'emanation
-c--------------------------------------------------------
-c
-        IF
-     .  ( (itr.eq.1.AND.NINT(masktr(i)).EQ.1.AND.lat(i).GE.60.
-     .     .AND.lat(i).LE.70.)
-     .                      .OR.
-     .     (itr.eq.2.AND.NINT(masktr(i)).EQ.1.AND.lat(i).GE.60.
-     .     .AND.lat(i).LE.70.) )
-     .  THEN
-         zx_trs(i) = local_trs(i)
-         zx_a = zx_trs(i)
-     .         +(fshtr(i)/2.)*dtime*rotrhi(i)
-     .         +rotrhi(i)*masktr(i)*zx_coef(i,1)/RG
-     .         *(zx_ctr(i,1)*(zx_alpha1(i)+zx_alpha2(i)*zx_dtr(i,2))
-     .         +zx_alpha2(i)*zx_ctr(i,2))
-         zx_b = 1. + rotrhi(i)*masktr(i)*zx_coef(i,1)/RG
-     .              * (1.-zx_dtr(i,1)
-     .              *(zx_alpha1(i)+zx_alpha2(i)*zx_dtr(i,2)))
-     .              + dtime / tautr
-     .              + dtime * vdeptr / hsoltr
-         zx_trs(i) = zx_a / zx_b
-         local_trs(i) = zx_trs(i)
-        ENDIF
-c
-c----------------------------------------------
-c Au dessus des oceans et aux hautes latitudes
-c----------------------------------------------
-c
-c au dessous de -60S  pas d'emission de radon au dessus 
-c des oceans et des continents
-c---------------------------------------------------------------
-
-       IF ( (itr.EQ.1.AND.NINT(masktr(i)).EQ.0)
-     .                 .OR.
-     . (itr.EQ.1.AND.NINT(masktr(i)).EQ.1.AND.lat(i).LT.-60.))
-     . THEN
-         zx_trs(i) = 0.
-         local_trs(i) = 0.
-        END IF
-
-c au dessus de 70 N pas d'emission de radon au dessus 
-c des oceans et des continents
-c--------------------------------------------------------------
-       IF ( (itr.EQ.1.AND.NINT(masktr(i)).EQ.0)
-     .                 .OR.
-     . (itr.EQ.1.AND.NINT(masktr(i)).EQ.1.AND.lat(i).GT.70.))
-     . THEN
-         zx_trs(i) = 0.
-         local_trs(i) = 0.
-        END IF
-
-c Au dessus des oceans la source est nulle
-c-----------------------------------------
-c
-        IF (itr.eq.1.AND.NINT(masktr(i)).EQ.0) THEN
-         zx_trs(i) = 0.
-         local_trs(i) = 0.
-        END IF
-c
-      ENDDO    ! sur le i=1,klon
-c
-c======================================================================
-c==== une fois on a zx_trs, on peut faire l'iteration ========
-c
-      DO i = 1, klon
-         local_tr(i,1) = zx_ctr(i,1)+zx_dtr(i,1)*zx_trs(i)
-      ENDDO
-      DO l = 2, klev
-      DO i = 1, klon
-        local_tr(i,l)
-     .    = zx_ctr(i,l) + zx_dtr(i,l)*local_tr(i,l-1)
-      ENDDO
-      ENDDO
-c======================================================================
-c== Calcul du flux de traceur (flux_tr): UA/(m**2 s)
-c
-      DO i = 1, klon
-        flux_tr(i,1) = masktr(i)*zx_coef(i,1)/RG
-     .      * (zx_alpha1(i)*local_tr(i,1)+zx_alpha2(i)*local_tr(i,2)
-     .                   -zx_trs(i)) / dtime
-      ENDDO
-      DO l = 2, klev
-      DO i = 1, klon
-        flux_tr(i,l) = zx_coef(i,l)/RG
-     .       * (local_tr(i,l)-local_tr(i,l-1)) / dtime
-      ENDDO
-      ENDDO
-c======================================================================
-c== Calcul des tendances du traceur ds le sol et dans l'atmosphere
-c
-      DO l = 1, klev
-      DO i = 1, klon
-         d_tr(i,l) = local_tr(i,l) - tr(i,l)
-      ENDDO
-      ENDDO
-      DO i = 1, klon
-         d_trs(i) = local_trs(i) - trs(i)
-      ENDDO
-c======================================================================
-c
-      RETURN
-      END
Index: /LMDZ4/branches/LMDZ4-dev/libf/phylmd/cltracrn.F90
===================================================================
--- /LMDZ4/branches/LMDZ4-dev/libf/phylmd/cltracrn.F90	(revision 1191)
+++ /LMDZ4/branches/LMDZ4-dev/libf/phylmd/cltracrn.F90	(revision 1191)
@@ -0,0 +1,286 @@
+!$Id $
+
+SUBROUTINE cltracrn( itr, dtime,u1lay, v1lay, &
+     cdrag,coef,t,ftsol,pctsrf,               &
+     tr,trs,paprs,pplay,delp,                 &
+     masktr,fshtr,hsoltr,tautr,vdeptr,        &
+     lat,d_tr,d_trs )
+  
+  USE dimphy
+  IMPLICIT NONE
+!======================================================================
+! Auteur(s): Alex/LMD) date:  fev 99
+!            inspire de clqh + clvent
+! Objet: diffusion verticale de traceurs avec quantite de traceur ds 
+!        le sol ( reservoir de sol de radon ) 
+!        
+! note : pour l'instant le traceur dans le sol et le flux sont
+!        calcules mais ils ne servent que de diagnostiques
+!        seule la tendance sur le traceur est sortie (d_tr)
+!---------------------------------------------------------------------
+! Arguments:
+! itr......input-R-  le type de traceur 1- Rn 2 - Pb
+! dtime....input-R-  intervalle du temps (en secondes) ~ pdtphys
+! u1lay....input-R-  vent u de la premiere couche (m/s)
+! v1lay....input-R-  vent v de la premiere couche (m/s)
+! cdrag....input-R-  cdrag
+! coef.....input-R-  le coefficient d'echange (m**2/s) l>1, valable uniquement pour k entre 2 et klev
+! t........input-R-  temperature (K)
+! paprs....input-R-  pression a inter-couche (Pa)
+! pplay....input-R-  pression au milieu de couche (Pa)
+! delp.....input-R-  epaisseur de couche (Pa)
+! ftsol....input-R-  temperature du sol (en Kelvin)
+! tr.......input-R-  traceurs
+! trs......input-R-  traceurs dans le sol
+! masktr...input-R-  Masque reservoir de sol traceur (1 = reservoir)
+! fshtr....input-R-  Flux surfacique de production dans le sol
+! tautr....input-R-  Constante de decroissance du traceur
+! vdeptr...input-R-  Vitesse de depot sec dans la couche brownienne
+! hsoltr...input-R-  Epaisseur equivalente du reservoir de sol
+! lat......input-R-  latitude en degree
+! d_tr.....output-R- le changement de "tr"
+! d_trs....output-R- le changement de "trs"
+!======================================================================
+  include "YOMCST.h"
+  include "indicesol.h"
+!
+!Entrees
+  INTEGER,INTENT(IN)                     :: itr
+  REAL,INTENT(IN)                        :: dtime
+  REAL,DIMENSION(klon),INTENT(IN)        :: u1lay, v1lay
+  REAL,DIMENSION(klon),INTENT(IN)        :: cdrag
+  REAL,DIMENSION(klon,klev),INTENT(IN)   :: coef, t
+  REAL,DIMENSION(klon,nbsrf),INTENT(IN)  :: ftsol, pctsrf 
+  REAL,DIMENSION(klon,klev),INTENT(IN)   :: tr 
+  REAL,DIMENSION(klon),INTENT(IN)        :: trs
+  REAL,DIMENSION(klon,klev+1),INTENT(IN) :: paprs 
+  REAL,DIMENSION(klon,klev),INTENT(IN)   :: pplay, delp
+  REAL,DIMENSION(klon),INTENT(IN)        :: masktr 
+  REAL,DIMENSION(klon),INTENT(IN)        :: fshtr 
+  REAL,INTENT(IN)                        :: hsoltr
+  REAL,INTENT(IN)                        :: tautr
+  REAL,INTENT(IN)                        :: vdeptr
+  REAL,DIMENSION(klon),INTENT(IN)        :: lat   
+
+!Sorties
+  REAL,DIMENSION(klon,klev),INTENT(OUT) :: d_tr
+  REAL,DIMENSION(klon),INTENT(OUT) :: d_trs  ! (diagnostic) traceur ds le sol
+
+!Locales
+  REAL,DIMENSION(klon,klev) :: flux_tr  ! (diagnostic) flux de traceur
+  INTEGER                   :: i, k, n, l
+  REAL,DIMENSION(klon)      :: rotrhi
+  REAL,DIMENSION(klon,klev) :: zx_coef
+  REAL,DIMENSION(klon)      :: zx_buf
+  REAL,DIMENSION(klon,klev) :: zx_ctr
+  REAL,DIMENSION(klon,klev) :: zx_dtr
+  REAL,DIMENSION(klon)      :: zx_trs
+  REAL                      :: zx_a, zx_b
+  
+  REAL,DIMENSION(klon,klev) :: local_tr
+  REAL,DIMENSION(klon)      :: local_trs
+  REAL,DIMENSION(klon)      :: zts      ! champ de temperature du sol
+  REAL,DIMENSION(klon)      :: zx_alpha1, zx_alpha2
+!======================================================================
+!AA Pour l'instant les 4 types de surface ne sont pas pris en compte
+!AA On fabrique avec zts un champ de temperature de sol  
+!AA que le pondere par la fraction de nature de sol.
+ 
+  DO i = 1,klon
+     zts(i) = 0. 
+  ENDDO
+
+  DO n=1,nbsrf
+     DO i = 1,klon
+        zts(i) = zts(i) + ftsol(i,n)*pctsrf(i,n)
+     ENDDO
+  ENDDO
+
+  DO i = 1,klon
+     rotrhi(i) = RD * zts(i) / hsoltr 
+  ENDDO
+
+  DO k = 1, klev
+     DO i = 1, klon
+        local_tr(i,k) = tr(i,k)
+     ENDDO
+  ENDDO
+
+  DO i = 1, klon
+     local_trs(i) = trs(i)
+  ENDDO
+!======================================================================
+!AA   Attention si dans clmain zx_alf1(i) = 1.0 
+!AA   Il doit y avoir coherence (dc la meme chose ici)
+
+  DO i = 1, klon
+!AA         zx_alpha1(i) = (paprs(i,1)-pplay(i,2))/(pplay(i,1)-pplay(i,2))
+     zx_alpha1(i) = 1.0
+     zx_alpha2(i) = 1.0 - zx_alpha1(i)
+  ENDDO
+!======================================================================
+  DO i = 1, klon
+     zx_coef(i,1) = cdrag(i)*(1.0+SQRT(u1lay(i)**2+v1lay(i)**2)) &
+          *pplay(i,1)/(RD*t(i,1))
+     zx_coef(i,1) = zx_coef(i,1) * dtime*RG
+  ENDDO
+
+  DO k = 2, klev
+     DO i = 1, klon
+        zx_coef(i,k) = coef(i,k)*RG/(pplay(i,k-1)-pplay(i,k)) &
+             *(paprs(i,k)*2/(t(i,k)+t(i,k-1))/RD)**2
+        zx_coef(i,k) = zx_coef(i,k) * dtime*RG
+     ENDDO
+  ENDDO
+!======================================================================
+  DO i = 1, klon
+     zx_buf(i)      = delp(i,klev) + zx_coef(i,klev)
+     zx_ctr(i,klev) = local_tr(i,klev)*delp(i,klev)/zx_buf(i)
+     zx_dtr(i,klev) = zx_coef(i,klev) / zx_buf(i)
+  ENDDO
+
+  DO l = klev-1, 2 , -1
+     DO i = 1, klon
+        zx_buf(i) = delp(i,l)+zx_coef(i,l)      &
+             +zx_coef(i,l+1)*(1.-zx_dtr(i,l+1))
+ 
+        zx_ctr(i,l) = ( local_tr(i,l)*delp(i,l) &
+             + zx_coef(i,l+1)*zx_ctr(i,l+1) )/zx_buf(i)
+        zx_dtr(i,l) = zx_coef(i,l) / zx_buf(i)
+     ENDDO
+  ENDDO
+
+  DO i = 1, klon
+     zx_buf(i) = delp(i,1) + zx_coef(i,2)*(1.-zx_dtr(i,2))  &
+          + masktr(i) * zx_coef(i,1)                        &
+          *( zx_alpha1(i)+zx_alpha2(i)*zx_dtr(i,2) )
+
+     zx_ctr(i,1) = ( local_tr(i,1)*delp(i,1)                &
+          + zx_ctr(i,2)                                     &
+          *(zx_coef(i,2)                                    &
+          - masktr(i) * zx_coef(i,1)                        &
+          *zx_alpha2(i) ) ) / zx_buf(i)
+     zx_dtr(i,1) = masktr(i) * zx_coef(i,1) / zx_buf(i)
+  ENDDO
+!======================================================================
+! Calculer d'abord local_trs nouvelle quantite dans le reservoir
+! de sol
+!=====================================================================
+
+  DO i = 1, klon
+!-------------------------
+! Au dessus des continents
+!--
+! Le pb peut se deposer partout : vdeptr = 10-3 m/s
+! Le Rn est traiter commme une couche Brownienne puisque vdeptr = 0.
+!-------------------------------------------------------------------
+     IF ( NINT(masktr(i)) .EQ. 1 ) THEN
+        zx_trs(i) = local_trs(i)
+        zx_a = zx_trs(i)                                           &
+             +fshtr(i)*dtime*rotrhi(i)                             &
+             +rotrhi(i)*masktr(i)*zx_coef(i,1)/RG                  &
+             *(zx_ctr(i,1)*(zx_alpha1(i)+zx_alpha2(i)*zx_dtr(i,2)) &
+             +zx_alpha2(i)*zx_ctr(i,2))
+! Pour l'instant, pour aller vite, le depot sec est traite comme une decroissance
+        zx_b = 1. + rotrhi(i)*masktr(i)*zx_coef(i,1)/RG            &
+             * (1.-zx_dtr(i,1)                                     &
+             *(zx_alpha1(i)+zx_alpha2(i)*zx_dtr(i,2)))             &
+             + dtime / tautr                                       & 
+             + dtime * vdeptr / hsoltr 
+        zx_trs(i) = zx_a / zx_b
+        local_trs(i) = zx_trs(i)
+     ENDIF
+!--------------------------------------------------------
+! Si on est entre 60N et 70N on divise par 2 l'emanation
+!--------------------------------------------------------
+
+     IF ( (itr.eq.1.AND.NINT(masktr(i)).EQ.1.AND.lat(i).GE.60..AND.lat(i).LE.70.).OR.      &
+          (itr.eq.2.AND.NINT(masktr(i)).EQ.1.AND.lat(i).GE.60..AND.lat(i).LE.70.) ) THEN
+        zx_trs(i) = local_trs(i)
+        zx_a = zx_trs(i)                                           &
+             +(fshtr(i)/2.)*dtime*rotrhi(i)                        & 
+             +rotrhi(i)*masktr(i)*zx_coef(i,1)/RG                  &
+             *(zx_ctr(i,1)*(zx_alpha1(i)+zx_alpha2(i)*zx_dtr(i,2)) &
+             +zx_alpha2(i)*zx_ctr(i,2))
+        !
+        zx_b = 1. + rotrhi(i)*masktr(i)*zx_coef(i,1)/RG  &
+             * (1.-zx_dtr(i,1)                           &
+             *(zx_alpha1(i)+zx_alpha2(i)*zx_dtr(i,2)))   &
+             + dtime / tautr                             &
+             + dtime * vdeptr / hsoltr
+        ! 
+        zx_trs(i) = zx_a / zx_b
+        local_trs(i) = zx_trs(i)
+     ENDIF
+
+!----------------------------------------------
+! Au dessus des oceans et aux hautes latitudes
+!--
+! au dessous de -60S  pas d'emission de radon au dessus 
+! des oceans et des continents
+!---------------------------------------------------------------
+
+     IF ( (itr.EQ.1.AND.NINT(masktr(i)).EQ.0).OR.       &
+          (itr.EQ.1.AND.NINT(masktr(i)).EQ.1.AND.lat(i).LT.-60.)) THEN
+        zx_trs(i) = 0.
+        local_trs(i) = 0.
+     END IF
+!--
+! au dessus de 70 N pas d'emission de radon au dessus 
+! des oceans et des continents
+!--------------------------------------------------------------
+     IF ( (itr.EQ.1.AND.NINT(masktr(i)).EQ.0).OR.    &
+          (itr.EQ.1.AND.NINT(masktr(i)).EQ.1.AND.lat(i).GT.70.)) THEN
+        zx_trs(i) = 0.
+        local_trs(i) = 0.
+     END IF
+!---------------------------------------------
+! Au dessus des oceans la source est nulle
+!--------------------------------------------
+
+     IF (itr.eq.1.AND.NINT(masktr(i)).EQ.0) THEN
+        zx_trs(i) = 0.
+        local_trs(i) = 0.
+     END IF
+
+  ENDDO    ! sur le i=1,klon
+!
+!======================================================================
+! Une fois on a zx_trs, on peut faire l'iteration 
+!====================================================================== 
+
+  DO i = 1, klon
+     local_tr(i,1) = zx_ctr(i,1)+zx_dtr(i,1)*zx_trs(i)
+  ENDDO
+  DO l = 2, klev
+     DO i = 1, klon
+        local_tr(i,l) = zx_ctr(i,l) + zx_dtr(i,l)*local_tr(i,l-1)
+     ENDDO
+  ENDDO
+!======================================================================
+! Calcul du flux de traceur (flux_tr): UA/(m**2 s)
+!======================================================================
+  DO i = 1, klon
+     flux_tr(i,1) = masktr(i)*zx_coef(i,1)/RG                      &
+          * (zx_alpha1(i)*local_tr(i,1)+zx_alpha2(i)*local_tr(i,2) &
+          -zx_trs(i)) / dtime
+  ENDDO
+  DO l = 2, klev
+     DO i = 1, klon
+        flux_tr(i,l) = zx_coef(i,l)/RG                    &
+             * (local_tr(i,l)-local_tr(i,l-1)) / dtime
+     ENDDO
+  ENDDO
+!======================================================================
+! Calcul des tendances du traceur ds le sol et dans l'atmosphere
+!======================================================================
+  DO l = 1, klev
+     DO i = 1, klon
+        d_tr(i,l) = local_tr(i,l) - tr(i,l)
+     ENDDO
+  ENDDO
+  DO i = 1, klon
+     d_trs(i) = local_trs(i) - trs(i)
+  ENDDO
+
+END SUBROUTINE cltracrn
Index: DZ4/branches/LMDZ4-dev/libf/phylmd/cvltr.F
===================================================================
--- /LMDZ4/branches/LMDZ4-dev/libf/phylmd/cvltr.F	(revision 1190)
+++ 	(revision )
@@ -1,145 +1,0 @@
-c
-c $Header$
-c
-      SUBROUTINE cvltr(pdtime,da, phi, mp,paprs,pplay,x,upd,dnd,dx)
-      USE dimphy
-      IMPLICIT NONE 
-c=====================================================================
-c Objet : convection des traceurs / KE
-c Auteurs: M-A Filiberti and J-Y Grandpeix
-c=====================================================================
-c
-cym#include "dimensions.h"
-cym#include "dimphy.h"
-#include "YOMCST.h"
-#include "YOECUMF.h" 
-c
-      REAL pdtime
-      REAL paprs(klon,klev+1) ! pression aux 1/2 couches (bas en haut)
-      REAL pplay(klon,klev)  ! pression pour le milieu de chaque couche
-      REAL x(klon,klev)        ! q de traceur (bas en haut) 
-      REAL dx(klon,klev)     ! tendance de traceur  (bas en haut)
-      real da(klon,klev),phi(klon,klev,klev),mp(klon,klev)
-      REAL upd(klon,klev)      ! saturated updraft mass flux
-      REAL dnd(klon,klev)      ! saturated downdraft mass flux
-c
-c--variables locales      
-      real zed(klon,klev),zmd(klon,klev,klev)
-      real za(klon,klev,klev)
-      real zmfd(klon,klev),zmfa(klon,klev)
-      real zmfp(klon,klev),zmfu(klon,klev)
-      integer i,k,j 
-c test conservation
-c      real conserv
-c =========================================
-c calcul des tendances liees au downdraft
-c =========================================
-      zed(:,:)=0.
-      zmfd(:,:)=0.
-      zmfa(:,:)=0.
-      zmfu(:,:)=0.
-      zmfp(:,:)=0.
-      zmd(:,:,:)=0.
-      za(:,:,:)=0.
-c entrainement
-      do k=1,klev-1
-        do i=1,klon
-          zed(i,k)=max(0.,mp(i,k)-mp(i,k+1))
-        end do
-      end do
-c
-c calcul de la matrice d echange
-c matrice de distribution de la masse entrainee en k
-c
-      do k=1,klev
-        do i=1,klon
-          zmd(i,k,k)=zed(i,k)
-        end do
-      end do
-      do k=2,klev
-        do j=k-1,1,-1
-          do i=1,klon
-          if(mp(i,j+1).ne.0) then
-          zmd(i,j,k)=zmd(i,j+1,k)*min(1.,mp(i,j)/mp(i,j+1))
-          endif
-          end do
-        end do
-      end do
-      do k=1,klev
-        do j=1,klev-1
-          do i=1,klon
-          za(i,j,k)=max(0.,zmd(i,j+1,k)-zmd(i,j,k))
-          end do
-        end do
-      end do
-c
-c rajout du terme lie a l ascendance induite
-c
-        do j=2,klev
-         do i=1,klon
-          za(i,j,j-1)=za(i,j,j-1)+mp(i,j)
-         end do
-        end do
-C
-c tendances
-c            
-      do k=1,klev
-        do j=1,klev
-          do i=1,klon
-          zmfd(i,j)=zmfd(i,j)+za(i,j,k)*(x(i,k)-x(i,j))
-          end do
-        end do
-      end do
-c
-c =========================================
-c calcul des tendances liees aux flux satures
-c =========================================
-      do j=1,klev
-        do i=1,klon
-          zmfa(i,j)=da(i,j)*(x(i,1)-x(i,j))
-        end do
-      end do
-      do k=1,klev
-        do j=1,klev
-          do i=1,klon
-          zmfp(i,j)=zmfp(i,j)+phi(i,j,k)*(x(i,k)-x(i,j))
-          end do
-        end do
-      end do
-      do j=1,klev-1
-        do i=1,klon
-          zmfu(i,j)=max(0.,upd(i,j+1)+dnd(i,j+1))*(x(i,j+1)-x(i,j))
-        end do
-      end do
-      do j=2,klev
-        do i=1,klon
-          zmfu(i,j)=zmfu(i,j)
-     .             +min(0.,upd(i,j)+dnd(i,j))*(x(i,j)-x(i,j-1))
-        end do
-      end do
-
-c =========================================
-c--calcul final des tendances
-c =========================================
-      do k=1, klev
-        do i=1, klon
-          dx(i,k)=(zmfd(i,k)+zmfu(i,k)
-     .      +zmfa(i,k)+zmfp(i,k))*pdtime
-     .      *RG/(paprs(i,k)-paprs(i,k+1))
-c          print*,'dx',k,dx(i,k)
-        enddo
-      enddo
-c
-c test de conservation du traceur
-c      conserv=0.
-c      do k=1, klev
-c        do i=1, klon
-c         conserv=conserv+dx(i,k)*
-c     .     (paprs(i,k)-paprs(i,k+1))/RG
-C
-c        enddo
-c      enddo
-c      print *,'conserv',conserv
-     
-      return
-      end
Index: /LMDZ4/branches/LMDZ4-dev/libf/phylmd/cvltr.F90
===================================================================
--- /LMDZ4/branches/LMDZ4-dev/libf/phylmd/cvltr.F90	(revision 1191)
+++ /LMDZ4/branches/LMDZ4-dev/libf/phylmd/cvltr.F90	(revision 1191)
@@ -0,0 +1,145 @@
+!
+! $Id $
+!
+SUBROUTINE cvltr(pdtime,da, phi, mp,paprs,pplay,x,upd,dnd,dx)
+  USE dimphy
+  IMPLICIT NONE 
+!=====================================================================
+! Objet : convection des traceurs / KE
+! Auteurs: M-A Filiberti and J-Y Grandpeix
+!=====================================================================
+
+  include "YOMCST.h"
+  include "YOECUMF.h" 
+
+! Entree
+  REAL,INTENT(IN)                           :: pdtime
+  REAL,DIMENSION(klon,klev),INTENT(IN)      :: da
+  REAL,DIMENSION(klon,klev,klev),INTENT(IN) :: phi
+  REAL,DIMENSION(klon,klev),INTENT(IN)      :: mp
+  REAL,DIMENSION(klon,klev+1),INTENT(IN)    :: paprs ! pression aux 1/2 couches (bas en haut)
+  REAL,DIMENSION(klon,klev),INTENT(IN)      :: pplay ! pression pour le milieu de chaque couche
+  REAL,DIMENSION(klon,klev),INTENT(IN)      :: x     ! q de traceur (bas en haut) 
+  REAL,DIMENSION(klon,klev),INTENT(IN)      :: upd   ! saturated updraft mass flux
+  REAL,DIMENSION(klon,klev),INTENT(IN)      :: dnd   ! saturated downdraft mass flux
+
+! Sortie
+  REAL,DIMENSION(klon,klev),INTENT(OUT) :: dx ! tendance de traceur  (bas en haut)
+
+! Variables locales     
+  REAL,DIMENSION(klon,klev)       :: zed
+  REAL,DIMENSION(klon,klev,klev)  :: zmd
+  REAL,DIMENSION(klon,klev,klev)  :: za
+  REAL,DIMENSION(klon,klev)       :: zmfd,zmfa
+  REAL,DIMENSION(klon,klev)       :: zmfp,zmfu
+  INTEGER                         :: i,k,j 
+
+! =========================================
+! calcul des tendances liees au downdraft
+! =========================================
+  zed(:,:)=0.
+  zmfd(:,:)=0.
+  zmfa(:,:)=0.
+  zmfu(:,:)=0.
+  zmfp(:,:)=0.
+  zmd(:,:,:)=0.
+  za(:,:,:)=0.
+! entrainement
+  DO k=1,klev-1
+     DO i=1,klon
+        zed(i,k)=max(0.,mp(i,k)-mp(i,k+1))
+     END DO
+  END DO
+
+! calcul de la matrice d echange
+! matrice de distribution de la masse entrainee en k
+
+  DO k=1,klev
+     DO i=1,klon
+        zmd(i,k,k)=zed(i,k)
+     END DO
+  END DO
+  DO k=2,klev
+     DO j=k-1,1,-1
+        DO i=1,klon
+           if(mp(i,j+1).ne.0) then
+              zmd(i,j,k)=zmd(i,j+1,k)*min(1.,mp(i,j)/mp(i,j+1))
+           ENDif
+        END DO
+     END DO
+  END DO
+  DO k=1,klev
+     DO j=1,klev-1
+        DO i=1,klon
+           za(i,j,k)=max(0.,zmd(i,j+1,k)-zmd(i,j,k))
+        END DO
+     END DO
+  END DO
+!
+! rajout du terme lie a l ascendance induite
+!
+  DO j=2,klev
+     DO i=1,klon
+        za(i,j,j-1)=za(i,j,j-1)+mp(i,j)
+     END DO
+  END DO
+!
+! tendances
+!            
+  DO k=1,klev
+     DO j=1,klev
+        DO i=1,klon
+           zmfd(i,j)=zmfd(i,j)+za(i,j,k)*(x(i,k)-x(i,j))
+        END DO
+     END DO
+  END DO
+!
+! =========================================
+! calcul des tendances liees aux flux satures
+! =========================================
+  DO j=1,klev
+     DO i=1,klon
+        zmfa(i,j)=da(i,j)*(x(i,1)-x(i,j))
+     END DO
+  END DO
+  DO k=1,klev
+     DO j=1,klev
+        DO i=1,klon
+           zmfp(i,j)=zmfp(i,j)+phi(i,j,k)*(x(i,k)-x(i,j))
+        END DO
+     END DO
+  END DO
+  DO j=1,klev-1
+     DO i=1,klon
+        zmfu(i,j)=max(0.,upd(i,j+1)+dnd(i,j+1))*(x(i,j+1)-x(i,j))
+     END DO
+  END DO
+  DO j=2,klev
+     DO i=1,klon
+        zmfu(i,j)=zmfu(i,j)+min(0.,upd(i,j)+dnd(i,j))*(x(i,j)-x(i,j-1))
+     END DO
+  END DO
+
+! =========================================
+! calcul final des tendances
+! =========================================
+  DO k=1, klev
+     DO i=1, klon
+        dx(i,k)=(zmfd(i,k)+zmfu(i,k)       &
+             +zmfa(i,k)+zmfp(i,k))*pdtime  &
+             *RG/(paprs(i,k)-paprs(i,k+1)) 
+        !          print*,'dx',k,dx(i,k)
+     ENDDO
+  ENDDO
+
+! test de conservation du traceur
+!      conserv=0.
+!      DO k=1, klev
+!        DO i=1, klon
+!         conserv=conserv+dx(i,k)*   &
+!        (paprs(i,k)-paprs(i,k+1))/RG
+!        ENDDO
+!      ENDDO
+!      print *,'conserv',conserv
+     
+END SUBROUTINE cvltr
Index: /LMDZ4/branches/LMDZ4-dev/libf/phylmd/ini_histrac.h
===================================================================
--- /LMDZ4/branches/LMDZ4-dev/libf/phylmd/ini_histrac.h	(revision 1190)
+++ /LMDZ4/branches/LMDZ4-dev/libf/phylmd/ini_histrac.h	(revision 1191)
@@ -1,123 +1,128 @@
 !
-! $Header$
+! $Id $
 !
-      IF (ecrit_tra>0. .AND. config_inca == 'none') THEN
-c$OMP MASTER 
-         CALL ymds2ju(annee_ref, 1, day_ref, 0.0, zjulian)
-c
-         CALL histbeg_phy("histrac", itau_phy, zjulian, pdtphys,
-     .                 nhori, nid_tra)
-         CALL histvert(nid_tra, "presnivs", "Vertical levels", "mb",
-     .                 klev, presnivs, nvert)
+#ifdef CPP_IOIPSL
 
+  IF (ecrit_tra>0. .AND. config_inca == 'none') THEN
+!$OMP MASTER 
+     CALL ymds2ju(annee_ref, 1, day_ref, 0.0, zjulian)
+     CALL histbeg_phy("histrac", itau_phy, zjulian, pdtphys,nhori, nid_tra)
+     CALL histvert(nid_tra, "presnivs", "Vertical levels", "mb",klev, presnivs, nvert)
 
+     zsto = pdtphys
+     zout = ecrit_tra
+     CALL histdef(nid_tra, "phis", "Surface geop. height", "-",   &
+          iim,jj_nb,nhori, 1,1,1, -99, 32,"once",  zsto,zout)
+     CALL histdef(nid_tra, "aire", "Grid area", "-",              &
+          iim,jj_nb,nhori, 1,1,1, -99, 32,"once",  zsto,zout)
 
-         zsto = pdtphys
-         zout = ecrit_tra
-c
-         CALL histdef(nid_tra, "phis", "Surface geop. height", "-",
-     .                iim,jj_nb,nhori, 1,1,1, -99, 32,
-     .                "once",  zsto,zout)
-c
-         CALL histdef(nid_tra, "aire", "Grid area", "-",
-     .                iim,jj_nb,nhori, 1,1,1, -99, 32,
-     .                "once",  zsto,zout)
-         DO it=1,nbtr
-C champ 2D
-         iq=it+2
-         iiq=niadv(iq)
-         CALL histdef(nid_tra, tname(iiq), ttext(iiq), "U/kga",
-     .                iim,jj_nb,nhori, klev,1,klev,nvert, 32,
-     .                "ave(X)", zsto,zout)
-         if (lessivage) THEN
-         CALL histdef(nid_tra, "fl"//tname(iiq),"Flux "//ttext(iiq),
-     .              "U/m2/s",iim,jj_nb,nhori, klev,1,klev,nvert, 32,
-     .              "ave(X)", zsto,zout)
-         endif
+!TRACEURS
+!----------------
+     DO it = 1,nbtr
+        iiq = niadv(it+2)
 
-c---Ajout Olivia
-         CALL histdef(nid_tra, "d_tr_th_"//tname(iiq),
-     .	              "tendance thermique"// ttext(iiq), "?",
-     .                iim,jj_nb,nhori, klev,1,klev,nvert, 32,
-     .                "ave(X)", zsto,zout)
-c
-         if(iflag_con.GE.2) then
-         CALL histdef(nid_tra, "d_tr_cv_"//tname(iiq),
-     .	              "tendance convection"// ttext(iiq), "?",
-     .                iim,jj_nb,nhori, klev,1,klev,nvert, 32,
-     .                "ave(X)", zsto,zout)
-	     endif !(iflag_con.GE.2) then
-         CALL histdef(nid_tra, "d_tr_cl_"//tname(iiq),
-     .	              "tendance couche limite"// ttext(iiq), "?",
-     .                iim,jj_nb,nhori, klev,1,klev,nvert, 32,
-     .                "ave(X)", zsto,zout)
-c---fin Olivia	 
+! CONCENTRATIONS
+        CALL histdef(nid_tra, tname(iiq), ttext(iiq), "U/kga",    &
+             iim,jj_nb,nhori, klev,1,klev,nvert, 32,"ave(X)", zsto,zout)
 
-         ENDDO
+! TD LESSIVAGE
+        IF (lessivage .AND. aerosol(it)) THEN
+           CALL histdef(nid_tra, "fl"//tname(iiq),"Flux "//ttext(iiq), &
+                "U/m2/s",iim,jj_nb,nhori, klev,1,klev,nvert, 32,       &
+                "ave(X)", zsto,zout)
+        END IF
 
-         CALL histdef(nid_tra, "pyu1", "Vent niv 1", "-",
-     .                iim,jj_nb,nhori, 1,1,1, -99, 32,
-     .                "inst(X)",  zout,zout)
+! TD THERMIQUES
+        IF (iflag_thermals.gt.0) THEN
+           CALL histdef(nid_tra, "d_tr_th_"//tname(iiq),      &
+                "tendance thermique"// ttext(iiq), "?",       &
+                iim,jj_nb,nhori, klev,1,klev,nvert, 32,       &
+                "ave(X)", zsto,zout)
+        ENDIF
 
-         CALL histdef(nid_tra, "pyv1", "Vent niv 1", "-",
-     .                iim,jj_nb,nhori, 1,1,1, -99, 32,
-     .                "inst(X)",  zout,zout)
-         CALL histdef(nid_tra, "psrf1", "nature sol", "-",
-     .                iim,jj_nb,nhori, 1,1,1, -99, 32,
-     .                "inst(X)",  zout,zout)
-         CALL histdef(nid_tra, "psrf2", "nature sol", "-",
-     .                iim,jj_nb,nhori, 1,1,1, -99, 32,
-     .                "inst(X)",  zout,zout)
-         CALL histdef(nid_tra, "psrf3", "nature sol", "-",
-     .                iim,jj_nb,nhori, 1,1,1, -99, 32,
-     .                "inst(X)",  zout,zout)
-         CALL histdef(nid_tra, "psrf4", "nature sol", "-",
-     .                iim,jj_nb,nhori, 1,1,1, -99, 32,
-     .                "inst(X)",  zout,zout)
-         CALL histdef(nid_tra, "ftsol1", "temper sol", "-",
-     .                iim,jj_nb,nhori, 1,1,1, -99, 32,
-     .                "inst(X)",  zout,zout)
-         CALL histdef(nid_tra, "ftsol2", "temper sol", "-",
-     .                iim,jj_nb,nhori, 1,1,1, -99, 32,
-     .                "inst(X)",  zout,zout)
-         CALL histdef(nid_tra, "ftsol3", "temper sol", "-",
-     .                iim,jj_nb,nhori, 1,1,1, -99, 32,
-     .                "inst",  zout,zout)
-         CALL histdef(nid_tra, "ftsol4", "temper sol", "-",
-     .                iim,jj_nb,nhori, 1,1,1, -99, 32,
-     .                "inst(X)",  zout,zout)
-         CALL histdef(nid_tra, "pplay", "flux u mont","-",
-     .                iim,jj_nb,nhori, klev,1,klev,nvert, 32,
-     .                "inst(X)", zout,zout)
-         CALL histdef(nid_tra, "t", "flux u mont","-",
-     .                iim,jj_nb,nhori, klev,1,klev,nvert, 32,
-     .                "inst(X)", zout,zout)
-         CALL histdef(nid_tra, "mfu", "flux u mont","-",
-     .                iim,jj_nb,nhori, klev,1,klev,nvert, 32,
-     .                "ave(X)", zsto,zout)
-         CALL histdef(nid_tra, "mfd", "flux u decen","-",
-     .                iim,jj_nb,nhori, klev,1,klev,nvert, 32,
-     .                "ave(X)", zsto,zout)
-         CALL histdef(nid_tra, "en_u", "flux u mont","-",
-     .                iim,jj_nb,nhori, klev,1,klev,nvert, 32,
-     .                "ave(X)", zsto,zout)
-         CALL histdef(nid_tra, "en_d", "flux u mont","-",
-     .                iim,jj_nb,nhori, klev,1,klev,nvert, 32,
-     .                "ave(X)", zsto,zout)
-         CALL histdef(nid_tra, "de_d", "flux u mont","-",
-     .                iim,jj_nb,nhori, klev,1,klev,nvert, 32,
-     .                "ave(X)", zsto,zout)
-         CALL histdef(nid_tra, "de_u", "flux u decen","-",
-     .                iim,jj_nb,nhori, klev,1,klev,nvert, 32,
-     .                "ave(X)", zsto,zout)
-         CALL histdef(nid_tra, "coefh", "turbulent coef","-",
-     .                iim,jj_nb,nhori, klev,1,klev,nvert, 32,
-     .                "ave(X)", zsto,zout)
+! TD CONVECTION
+        IF (iflag_con.GE.2) THEN
+           CALL histdef(nid_tra, "d_tr_cv_"//tname(iiq),   &
+                "tendance convection"// ttext(iiq), "?",   &
+                iim,jj_nb,nhori, klev,1,klev,nvert, 32,    &
+                "ave(X)", zsto,zout)
+        ENDIF
 
-c
-         CALL histend(nid_tra)
-         ndex2d = 0
-         ndex3d = 0
-         ndex = 0
-c$OMP END MASTER
-      END IF ! ecrit_tra>0. .AND. config_inca == 'none'
+! TD COUCHE-LIMITE
+        CALL histdef(nid_tra, "d_tr_cl_"//tname(iiq),      &
+             "tendance couche limite"// ttext(iiq), "?",   &
+             iim,jj_nb,nhori, klev,1,klev,nvert, 32,       &
+             "ave(X)", zsto,zout)
+     ENDDO
+!---------------   
+!
+! VENT (niveau 1)
+     CALL histdef(nid_tra, "pyu1", "Vent niv 1", "-",      &
+          iim,jj_nb,nhori, 1,1,1, -99, 32,                 &
+          "inst(X)",  zout,zout)     
+     CALL histdef(nid_tra, "pyv1", "Vent niv 1", "-",      &
+          iim,jj_nb,nhori, 1,1,1, -99, 32,                 &
+          "inst(X)",  zout,zout)
+
+! TEMPERATURE DU SOL
+     CALL histdef(nid_tra, "ftsol1", "temper sol", "-",    &
+          iim,jj_nb,nhori, 1,1,1, -99, 32,                 &
+          "inst(X)",  zout,zout)
+     CALL histdef(nid_tra, "ftsol2", "temper sol", "-",    &
+          iim,jj_nb,nhori, 1,1,1, -99, 32,                 &
+          "inst(X)",  zout,zout)
+     CALL histdef(nid_tra, "ftsol3", "temper sol", "-",    &
+          iim,jj_nb,nhori, 1,1,1, -99, 32,                 &
+          "inst",  zout,zout)
+     CALL histdef(nid_tra, "ftsol4", "temper sol", "-",    &
+          iim,jj_nb,nhori, 1,1,1, -99, 32,                 &
+          "inst(X)",  zout,zout)
+
+! NATURE DU SOL
+     CALL histdef(nid_tra, "psrf1", "nature sol", "-",     &
+          iim,jj_nb,nhori, 1,1,1, -99, 32,                 &
+          "inst(X)",  zout,zout)
+     CALL histdef(nid_tra, "psrf2", "nature sol", "-",     &
+          iim,jj_nb,nhori, 1,1,1, -99, 32,                 &
+          "inst(X)",  zout,zout)
+     CALL histdef(nid_tra, "psrf3", "nature sol", "-",     &
+          iim,jj_nb,nhori, 1,1,1, -99, 32,                 &
+          "inst(X)",  zout,zout)
+     CALL histdef(nid_tra, "psrf4", "nature sol", "-",     &
+          iim,jj_nb,nhori, 1,1,1, -99, 32,                 & 
+          "inst(X)",  zout,zout)
+! DIVERS
+     CALL histdef(nid_tra, "pplay", "flux u mont","-",     &
+          iim,jj_nb,nhori, klev,1,klev,nvert, 32,          &
+          "inst(X)", zout,zout)
+     CALL histdef(nid_tra, "t", "flux u mont","-",         &
+          iim,jj_nb,nhori, klev,1,klev,nvert, 32,          &
+          "inst(X)", zout,zout)
+     CALL histdef(nid_tra, "mfu", "flux u mont","-",       &
+          iim,jj_nb,nhori, klev,1,klev,nvert, 32,          &
+          "ave(X)", zsto,zout)
+     CALL histdef(nid_tra, "mfd", "flux u decen","-",      &
+          iim,jj_nb,nhori, klev,1,klev,nvert, 32,          &
+          "ave(X)", zsto,zout)
+     CALL histdef(nid_tra, "en_u", "flux u mont","-",      &
+          iim,jj_nb,nhori, klev,1,klev,nvert, 32,          &
+          "ave(X)", zsto,zout)
+     CALL histdef(nid_tra, "en_d", "flux u mont","-",      &
+          iim,jj_nb,nhori, klev,1,klev,nvert, 32,          &
+          "ave(X)", zsto,zout)
+     CALL histdef(nid_tra, "de_d", "flux u mont","-",      &
+          iim,jj_nb,nhori, klev,1,klev,nvert, 32,          &
+          "ave(X)", zsto,zout)
+     CALL histdef(nid_tra, "de_u", "flux u decen","-",     &
+          iim,jj_nb,nhori, klev,1,klev,nvert, 32,          &
+          "ave(X)", zsto,zout)
+     CALL histdef(nid_tra, "coefh", "turbulent coef","-",  &
+          iim,jj_nb,nhori, klev,1,klev,nvert, 32,          &
+          "ave(X)", zsto,zout)   
+     
+     CALL histend(nid_tra)
+!$OMP END MASTER
+  END IF ! ecrit_tra>0. .AND. config_inca == 'none'
+
+#endif
+  
Index: /LMDZ4/branches/LMDZ4-dev/libf/phylmd/init_be.F90
===================================================================
--- /LMDZ4/branches/LMDZ4-dev/libf/phylmd/init_be.F90	(revision 1191)
+++ /LMDZ4/branches/LMDZ4-dev/libf/phylmd/init_be.F90	(revision 1191)
@@ -0,0 +1,510 @@
+!$Id $
+
+SUBROUTINE init_be(pctsrf,masktr,tautr,vdeptr,scavtr,srcbe)
+
+  USE dimphy
+  USE comgeomphy
+  USE infotrac, ONLY : nbtr
+    
+  IMPLICIT NONE 
+!=====================================================================
+! Objet : prescription d'une source de Beryllium 7 
+!         pour 19 niveaux verticaux
+!        (d'apres le diagramme de Lal and Peters, 1967)
+!
+!
+! written by : O. Coindreau (CEA/LDG) 05/2005
+! last modified by : A. Jamelot (LMD/CEA)  04/03/2009 
+!=====================================================================
+
+  INCLUDE "YOMCST.h"
+  INCLUDE "YOECUMF.h" 
+  INCLUDE "indicesol.h"
+
+!
+! Input Arguments
+!
+  REAL,DIMENSION(klon,nbsrf),INTENT(IN) :: pctsrf !Pourcentage de sol (f(nature du sol))
+!
+! Output Arguments
+!
+  REAL,DIMENSION(klon),INTENT(OUT)      :: masktr ! Masque de l'echange avec la surface (possible => 1 )
+  REAL,INTENT(OUT)                      :: tautr  ! Constante de decroissance radioactive
+  REAL,INTENT(OUT)                      :: vdeptr ! Vitesse de depot sec dans la couche Brownienne
+  REAL,INTENT(OUT)                      :: scavtr ! Coefficient de lessivage
+  REAL,DIMENSION(klon,klev),INTENT(OUT) :: srcbe  ! source volumique de 7Be      
+!
+! Local Variables
+!
+  REAL,DIMENSION(klon) :: rlatgeo   ! latitudes geomagnetiques de la grille
+  REAL                 :: glt       ! latitude du pole geomagnetique
+  REAL                 :: glg       ! longitude du pole geomagnetique
+  REAL                 :: latgeo,qcos
+  INTEGER              :: k,i
+
+  WRITE(*,*)'PASSAGE init_be ...'
+
+! Source actuellement definie pour klev = 19 et klev >= 39
+  IF (klev /= 19 .AND. klev<39) CALL abort_gcm("init_be","Source du be7 necessite klev=19 ou klev>=39",1)
+!
+! Definition des constantes
+! -------------------------
+  tautr = 6645000.
+  vdeptr = 1.E-3 
+  scavtr = 0.5 
+
+  WRITE(*,*) '-------------- SOURCE DE BERYLLIUM ------------------- '
+  WRITE(*,*)'Decroissance (s): ', tautr
+  WRITE(*,*)'Vitesse de depot sec: ',vdeptr
+  WRITE(*,*)'Facteur de lessivage: ',scavtr
+
+  DO i = 1,klon
+     masktr(i) = 0.
+     IF ( NINT(pctsrf(i,1)) .EQ. 1 ) masktr(i) = 1.
+  END DO
+
+! Premiers niveaux: source nulle
+! ------------------------------
+  DO k = 1,6
+     DO i = 1,klon
+        srcbe(i,k) = 0.
+     END DO
+  END DO
+!
+! Pour les autres niveaux:
+! 1-passer des coordonnees geographiques a la latitude geomagnetique
+! 2-prescrire la source de Be (en 10exp5 at/g/s) dans ce repere
+! 3-mettre la source de Be ds la bonne unite (en at/kgA/s)
+!
+  glt=78.5*rpi/180.
+  glg=69.0*rpi/180.
+
+  DO i = 1,klon
+     qcos=sin(glt)*sin(rlatd(i))
+     qcos=qcos+cos(glt)*cos(rlatd(i))*cos(rlond(i)+glg)
+     IF ( qcos .LT. -1.) qcos = -1.
+     IF ( qcos .GT. 1.)  qcos = 1.
+     rlatgeo(i)=rpi/2.-acos(qcos)
+  ENDDO
+
+!===========================
+!  Cas 19 niveaux verticaux
+!===========================
+  IF (klev.eq.19) then
+     DO k = 1,klev
+        DO i = 1,klon
+           latgeo=(180./rpi)*abs(rlatgeo(i))
+           IF ( k .EQ. 1 ) THEN
+              IF (latgeo.GE.50.0) srcbe(i,k)=0.1
+              IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=0.09
+              IF (latgeo.GE.30.0 .and. latgeo.LT.40.0) srcbe(i,k)=0.07
+           END IF
+           IF ( k .EQ. 2 ) THEN
+              IF (latgeo.GE.50.0) srcbe(i,k)=0.12
+              IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=0.1
+              IF (latgeo.GE.30.0 .and. latgeo.LT.40.0) srcbe(i,k)=0.09
+              IF (latgeo.GE.20.0 .and. latgeo.LT.30.0) srcbe(i,k)=0.07
+           END IF
+           IF ( k .EQ. 3 ) THEN
+              IF (latgeo.GE.50.0) srcbe(i,k)=0.14
+              IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=0.12
+              IF (latgeo.GE.30.0 .and. latgeo.LT.40.0) srcbe(i,k)=0.1
+              IF (latgeo.GE.20.0 .and. latgeo.LT.30.0) srcbe(i,k)=0.09
+              IF (latgeo.GE.10.0 .and. latgeo.LT.20.0) srcbe(i,k)=0.07
+           END IF
+           IF ( k .EQ. 4 ) THEN
+              IF (latgeo.GE.50.0) srcbe(i,k)=0.175
+              IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=0.16
+              IF (latgeo.GE.30.0 .and. latgeo.LT.40.0) srcbe(i,k)=0.14
+              IF (latgeo.GE.20.0 .and. latgeo.LT.30.0) srcbe(i,k)=0.12
+              IF (latgeo.GE.10.0 .and. latgeo.LT.20.0) srcbe(i,k)=0.1
+              IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=0.09
+           END IF
+           IF ( k .EQ. 5 ) THEN
+              IF (latgeo.GE.50.0) srcbe(i,k)=0.28
+              IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=0.26
+              IF (latgeo.GE.30.0 .and. latgeo.LT.40.0) srcbe(i,k)=0.23
+              IF (latgeo.GE.20.0 .and. latgeo.LT.30.0) srcbe(i,k)=0.175
+              IF (latgeo.GE.10.0 .and. latgeo.LT.20.0) srcbe(i,k)=0.14
+              IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=0.12
+           END IF
+           IF ( k .EQ. 6 ) THEN
+              IF (latgeo.GE.50.0) srcbe(i,k)=0.56
+              IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=0.49
+              IF (latgeo.GE.30.0 .and. latgeo.LT.40.0) srcbe(i,k)=0.42
+              IF (latgeo.GE.20.0 .and. latgeo.LT.30.0) srcbe(i,k)=0.28
+              IF (latgeo.GE.10.0 .and. latgeo.LT.20.0) srcbe(i,k)=0.26
+              IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=0.245
+           END IF
+           IF ( k .EQ. 7 ) THEN
+              IF (latgeo.GE.50.0) srcbe(i,k)=1.05
+              IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=0.875
+              IF (latgeo.GE.30.0 .and. latgeo.LT.40.0) srcbe(i,k)=0.7
+              IF (latgeo.GE.20.0 .and. latgeo.LT.30.0) srcbe(i,k)=0.52
+              IF (latgeo.GE.10.0 .and. latgeo.LT.20.0) srcbe(i,k)=0.44
+              IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=0.385
+           END IF
+           IF ( k .EQ. 8 ) THEN
+              IF (latgeo.GE.50.0) srcbe(i,k)=2.
+              IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=1.8
+              IF (latgeo.GE.30.0 .and. latgeo.LT.40.0) srcbe(i,k)=1.5
+              IF (latgeo.GE.20.0 .and. latgeo.LT.30.0) srcbe(i,k)=1.
+              IF (latgeo.GE.10.0 .and. latgeo.LT.20.0) srcbe(i,k)=0.8
+              IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=0.75
+           END IF
+           IF ( k .EQ. 9 ) THEN
+              IF (latgeo.GE.50.0) srcbe(i,k)=4.
+              IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=3.5
+              IF (latgeo.GE.30.0 .and. latgeo.LT.40.0) srcbe(i,k)=3.
+              IF (latgeo.GE.20.0 .and. latgeo.LT.30.0) srcbe(i,k)=2.5
+              IF (latgeo.GE.10.0 .and. latgeo.LT.20.0) srcbe(i,k)=1.8
+              IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=1.4
+           END IF
+           IF ( k .EQ. 10 ) THEN
+              IF (latgeo.GE.50.0) srcbe(i,k)=8.5
+              IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=8.
+              IF (latgeo.GE.30.0 .and. latgeo.LT.40.0) srcbe(i,k)=7.
+              IF (latgeo.GE.20.0 .and. latgeo.LT.30.0) srcbe(i,k)=4.5
+              IF (latgeo.GE.10.0 .and. latgeo.LT.20.0) srcbe(i,k)=3.5
+              IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=3.
+           END IF
+           IF ( k .EQ. 11 ) THEN
+              IF (latgeo.GE.50.0) srcbe(i,k)=17.
+              IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=15.
+              IF (latgeo.GE.30.0 .and. latgeo.LT.40.0) srcbe(i,k)=11.
+              IF (latgeo.GE.20.0 .and. latgeo.LT.30.0) srcbe(i,k)=8.
+              IF (latgeo.GE.10.0 .and. latgeo.LT.20.0) srcbe(i,k)=5.
+              IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=4.
+           END IF
+           IF ( k .EQ. 12 ) THEN
+              IF (latgeo.GE.50.0) srcbe(i,k)=25.
+              IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=22.
+              IF (latgeo.GE.30.0 .and. latgeo.LT.40.0) srcbe(i,k)=17.
+              IF (latgeo.GE.20.0 .and. latgeo.LT.30.0) srcbe(i,k)=11.
+              IF (latgeo.GE.10.0 .and. latgeo.LT.20.0) srcbe(i,k)=7.5
+              IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=7.
+           END IF
+           IF ( k .EQ. 13 ) THEN
+              IF (latgeo.GE.60.0) srcbe(i,k)=33.
+              IF (latgeo.GE.50.0 .and. latgeo.LT.60.0) srcbe(i,k)=32.
+              IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=30.
+              IF (latgeo.GE.30.0 .and. latgeo.LT.40.0) srcbe(i,k)=22.
+              IF (latgeo.GE.20.0 .and. latgeo.LT.30.0) srcbe(i,k)=15.
+              IF (latgeo.GE.10.0 .and. latgeo.LT.20.0) srcbe(i,k)=11.
+              IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=8.
+           END IF
+           IF ( k .EQ. 14 ) THEN
+              IF (latgeo.GE.60.0) srcbe(i,k)=48.
+              IF (latgeo.GE.50.0 .and. latgeo.LT.60.0) srcbe(i,k)=45.
+              IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=36.
+              IF (latgeo.GE.30.0 .and. latgeo.LT.40.0) srcbe(i,k)=26.
+              IF (latgeo.GE.20.0 .and. latgeo.LT.30.0) srcbe(i,k)=17.5
+              IF (latgeo.GE.10.0 .and. latgeo.LT.20.0) srcbe(i,k)=12.5
+              IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=10.
+           END IF
+           IF ( k .EQ. 15 ) THEN
+              IF (latgeo.GE.70.0) srcbe(i,k)=58.
+              IF (latgeo.GE.60.0 .and. latgeo.LT.70.0) srcbe(i,k)=57.
+              IF (latgeo.GE.50.0 .and. latgeo.LT.60.0) srcbe(i,k)=50.
+              IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=38.
+              IF (latgeo.GE.30.0 .and. latgeo.LT.40.0) srcbe(i,k)=25.
+              IF (latgeo.GE.20.0 .and. latgeo.LT.30.0) srcbe(i,k)=15.
+              IF (latgeo.GE.10.0 .and. latgeo.LT.20.0) srcbe(i,k)=12.5
+              IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=10.
+           END IF
+           IF ( k .EQ. 16 ) THEN
+              IF (latgeo.GE.70.0) srcbe(i,k)=70.
+              IF (latgeo.GE.60.0 .and. latgeo.LT.70.0) srcbe(i,k)=65.
+              IF (latgeo.GE.50.0 .and. latgeo.LT.60.0) srcbe(i,k)=50.
+              IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=32.
+              IF (latgeo.GE.30.0 .and. latgeo.LT.40.0) srcbe(i,k)=20.
+              IF (latgeo.GE.20.0 .and. latgeo.LT.30.0) srcbe(i,k)=13.
+              IF (latgeo.GE.10.0 .and. latgeo.LT.20.0) srcbe(i,k)=9.
+              IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=7.5
+           END IF
+           IF ( k .GE. 17 ) THEN
+              IF (latgeo.GE.70.0) srcbe(i,k)=80.
+              IF (latgeo.GE.60.0 .and. latgeo.LT.70.0) srcbe(i,k)=70.
+              IF (latgeo.GE.50.0 .and. latgeo.LT.60.0) srcbe(i,k)=45.
+              IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=27.
+              IF (latgeo.GE.30.0 .and. latgeo.LT.40.0) srcbe(i,k)=17.5
+              IF (latgeo.GE.20.0 .and. latgeo.LT.30.0) srcbe(i,k)=12.
+              IF (latgeo.GE.10.0 .and. latgeo.LT.20.0) srcbe(i,k)=8.
+              IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=7.
+           END IF
+        END DO
+     END DO
+  END IF ! fin de 19 niveaux verticaux
+
+!================================
+!  Cas 39 niveaux verticaux
+!================================
+  IF (klev .ge. 39) then
+     DO k = 1,klev
+        DO i = 1,klon
+           latgeo=(180./rpi)*abs(rlatgeo(i))
+           IF ( k .LE. 4 ) THEN
+              IF (latgeo.GE.50.0) srcbe(i,k)=0.07
+           END IF
+           IF ( k .EQ. 5 ) THEN
+              IF (latgeo.GE.50.0) srcbe(i,k)=0.1
+              IF (latgeo.GE.20.0 .and. latgeo.LT.50.0) srcbe(i,k)=0.09
+              IF (latgeo.GE.20.0 .and. latgeo.LT.30.0) srcbe(i,k)=0.07
+           END IF
+           IF ( k .EQ. 6 ) THEN
+              IF (latgeo.GE.50.0) srcbe(i,k)=0.14
+              IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=0.12
+              IF (latgeo.GE.30.0 .and. latgeo.LT.40.0) srcbe(i,k)=0.1
+              IF (latgeo.GE.20.0 .and. latgeo.LT.30.0) srcbe(i,k)=0.09
+              IF (latgeo.GE.10.0 .and. latgeo.LT.20.0) srcbe(i,k)=0.07
+           END IF
+           IF ( k .EQ. 7 ) THEN
+              IF (latgeo.GE.50.0) srcbe(i,k)=0.16
+              IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=0.16
+              IF (latgeo.GE.30.0 .and. latgeo.LT.40.0) srcbe(i,k)=0.14
+              IF (latgeo.GE.20.0 .and. latgeo.LT.30.0) srcbe(i,k)=0.12
+              IF (latgeo.GE.10.0 .and. latgeo.LT.20.0) srcbe(i,k)=0.1
+              IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=0.09
+           END IF
+           IF ( k .EQ. 8 ) THEN
+              IF (latgeo.GE.50.0) srcbe(i,k)=0.175
+              IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=0.16
+              IF (latgeo.GE.30.0 .and. latgeo.LT.40.0) srcbe(i,k)=0.14
+              IF (latgeo.GE.20.0 .and. latgeo.LT.30.0) srcbe(i,k)=0.12
+              IF (latgeo.GE.10.0 .and. latgeo.LT.20.0) srcbe(i,k)=0.1
+              IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=0.1
+           END IF
+           IF ( k .EQ. 9 ) THEN
+              IF (latgeo.GE.50.0) srcbe(i,k)=0.245
+              IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=0.21
+              IF (latgeo.GE.30.0 .and. latgeo.LT.40.0) srcbe(i,k)=0.175
+              IF (latgeo.GE.20.0 .and. latgeo.LT.30.0) srcbe(i,k)=0.14
+              IF (latgeo.GE.10.0 .and. latgeo.LT.20.0) srcbe(i,k)=0.12
+              IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=0.12
+           END IF
+           IF ( k .EQ. 10 ) THEN
+              IF (latgeo.GE.50.0) srcbe(i,k)=0.31
+              IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=0.28
+              IF (latgeo.GE.30.0 .and. latgeo.LT.40.0) srcbe(i,k)=0.245
+              IF (latgeo.GE.20.0 .and. latgeo.LT.30.0) srcbe(i,k)=0.21
+              IF (latgeo.GE.10.0 .and. latgeo.LT.20.0) srcbe(i,k)=0.16
+              IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=0.14
+           END IF
+           IF ( k .EQ. 11 ) THEN
+              IF (latgeo.GE.50.0) srcbe(i,k)=0.35
+              IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=0.3
+              IF (latgeo.GE.30.0 .and. latgeo.LT.40.0) srcbe(i,k)=0.3
+              IF (latgeo.GE.20.0 .and. latgeo.LT.30.0) srcbe(i,k)=0.2
+              IF (latgeo.GE.10.0 .and. latgeo.LT.20.0) srcbe(i,k)=0.18
+              IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=0.16
+           END IF
+           IF ( k .EQ. 12 ) THEN
+              IF (latgeo.GE.40.0) srcbe(i,k)=0.5
+              IF (latgeo.GE.30.0 .and. latgeo.LT.40.0) srcbe(i,k)=0.4
+              IF (latgeo.GE.20.0 .and. latgeo.LT.30.0) srcbe(i,k)=0.35
+              IF (latgeo.GE.10.0 .and. latgeo.LT.20.0) srcbe(i,k)=0.3
+              IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=0.25
+           END IF
+           IF ( k .EQ. 13 ) THEN
+              IF (latgeo.GE.50.0) srcbe(i,k)=0.8
+              IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=0.7
+              IF (latgeo.GE.30.0 .and. latgeo.LT.40.0) srcbe(i,k)=0.6
+              IF (latgeo.GE.20.0 .and. latgeo.LT.30.0) srcbe(i,k)=0.5
+              IF (latgeo.GE.10.0 .and. latgeo.LT.20.0) srcbe(i,k)=0.4
+              IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=0.35
+           END IF
+           IF ( k .EQ. 14 ) THEN
+              IF (latgeo.GE.50.0) srcbe(i,k)=1.2
+              IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=1.
+              IF (latgeo.GE.30.0 .and. latgeo.LT.40.0) srcbe(i,k)=0.75
+              IF (latgeo.GE.20.0 .and. latgeo.LT.30.0) srcbe(i,k)=0.6
+              IF (latgeo.GE.10.0 .and. latgeo.LT.20.0) srcbe(i,k)=0.5
+              IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=0.4
+           END IF
+           IF ( k .EQ. 15 ) THEN
+              IF (latgeo.GE.60.0) srcbe(i,k)=1.75
+              IF (latgeo.GE.50.0 .and. latgeo.LT.60.0) srcbe(i,k)=1.8 
+              IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=1.6
+              IF (latgeo.GE.30.0 .and. latgeo.LT.40.0) srcbe(i,k)=1.4
+              IF (latgeo.GE.20.0 .and. latgeo.LT.30.0) srcbe(i,k)=0.9
+              IF (latgeo.GE.10.0 .and. latgeo.LT.20.0) srcbe(i,k)=0.75
+              IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=0.65
+           END IF
+           IF ( k .EQ. 16 ) THEN
+              IF (latgeo.GE.50.0) srcbe(i,k)=3.
+              IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=2.5
+              IF (latgeo.GE.30.0 .and. latgeo.LT.40.0) srcbe(i,k)=1.8
+              IF (latgeo.GE.20.0 .and. latgeo.LT.30.0) srcbe(i,k)=1.5
+              IF (latgeo.GE.10.0 .and. latgeo.LT.20.0) srcbe(i,k)=1.2
+              IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=0.9
+           END IF
+           IF ( k .EQ. 17 ) THEN
+              IF (latgeo.GE.50.0) srcbe(i,k)=4.
+              IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=3.
+              IF (latgeo.GE.30.0 .and. latgeo.LT.40.0) srcbe(i,k)=2.5
+              IF (latgeo.GE.20.0 .and. latgeo.LT.30.0) srcbe(i,k)=2.
+              IF (latgeo.GE.10.0 .and. latgeo.LT.20.0) srcbe(i,k)=1.6
+              IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=1.4
+           END IF
+           IF ( k .EQ. 18 ) THEN
+              IF (latgeo.GE.50.0) srcbe(i,k)=7.
+              IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=6.
+              IF (latgeo.GE.30.0 .and. latgeo.LT.40.0) srcbe(i,k)=4.5
+              IF (latgeo.GE.20.0 .and. latgeo.LT.30.0) srcbe(i,k)=3.5
+              IF (latgeo.GE.10.0 .and. latgeo.LT.20.0) srcbe(i,k)=3.
+              IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=2.
+           END IF
+           IF ( k .EQ. 19 ) THEN
+              IF (latgeo.GE.50.0) srcbe(i,k)=8.5
+              IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=8.
+              IF (latgeo.GE.30.0 .and. latgeo.LT.40.0) srcbe(i,k)=7.
+              IF (latgeo.GE.20.0 .and. latgeo.LT.30.0) srcbe(i,k)=4.
+              IF (latgeo.GE.10.0 .and. latgeo.LT.20.0) srcbe(i,k)=3.5
+              IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=3.
+           END IF
+           IF ( k .EQ. 20 ) THEN
+              IF (latgeo.GE.50.0) srcbe(i,k)=12.5
+              IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=12.
+              IF (latgeo.GE.30.0 .and. latgeo.LT.40.0) srcbe(i,k)=8.
+              IF (latgeo.GE.20.0 .and. latgeo.LT.30.0) srcbe(i,k)=6.
+              IF (latgeo.GE.10.0 .and. latgeo.LT.20.0) srcbe(i,k)=4.
+              IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=3.5
+           END IF
+           IF ( k .EQ. 21 ) THEN
+              IF (latgeo.GE.50.0) srcbe(i,k)=16.
+              IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=13.
+              IF (latgeo.GE.30.0 .and. latgeo.LT.40.0) srcbe(i,k)=10.
+              IF (latgeo.GE.20.0 .and. latgeo.LT.30.0) srcbe(i,k)=7.5
+              IF (latgeo.GE.10.0 .and. latgeo.LT.20.0) srcbe(i,k)=4.5
+              IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=4.
+           END IF
+           IF ( k .EQ. 22 ) THEN
+              IF (latgeo.GE.50.0) srcbe(i,k)=20.
+              IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=17.5
+              IF (latgeo.GE.30.0 .and. latgeo.LT.40.0) srcbe(i,k)=12.5
+              IF (latgeo.GE.20.0 .and. latgeo.LT.30.0) srcbe(i,k)=9.
+              IF (latgeo.GE.10.0 .and. latgeo.LT.20.0) srcbe(i,k)=6.
+              IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=4.5
+           END IF
+           IF ( k .EQ. 23 ) THEN
+              IF (latgeo.GE.50.0) srcbe(i,k)=25.
+              IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=22.
+              IF (latgeo.GE.30.0 .and. latgeo.LT.40.0) srcbe(i,k)=15.
+              IF (latgeo.GE.20.0 .and. latgeo.LT.30.0) srcbe(i,k)=10.
+              IF (latgeo.GE.10.0 .and. latgeo.LT.20.0) srcbe(i,k)=7.5
+              IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=6.
+           END IF
+           IF ( k .EQ. 24 ) THEN
+              IF (latgeo.GE.50.0) srcbe(i,k)=28.
+              IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=26.
+              IF (latgeo.GE.30.0 .and. latgeo.LT.40.0) srcbe(i,k)=18.
+              IF (latgeo.GE.20.0 .and. latgeo.LT.30.0) srcbe(i,k)=12.
+              IF (latgeo.GE.10.0 .and. latgeo.LT.20.0) srcbe(i,k)=8.5
+              IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=7.
+           END IF
+           IF ( k .EQ. 25 ) THEN
+              IF (latgeo.GE.50.0) srcbe(i,k)=33.
+              IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=28.
+              IF (latgeo.GE.30.0 .and. latgeo.LT.40.0) srcbe(i,k)=20.
+              IF (latgeo.GE.20.0 .and. latgeo.LT.30.0) srcbe(i,k)=14.
+              IF (latgeo.GE.10.0 .and. latgeo.LT.20.0) srcbe(i,k)=10.
+              IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=8.5
+           END IF
+           IF ( k .EQ. 26 ) THEN
+              IF (latgeo.GE.60.0) srcbe(i,k)=38.
+              IF (latgeo.GE.50.0 .and. latgeo.LT.60.0) srcbe(i,k)=36.
+              IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=32.
+              IF (latgeo.GE.30.0 .and. latgeo.LT.40.0) srcbe(i,k)=24.
+              IF (latgeo.GE.20.0 .and. latgeo.LT.30.0) srcbe(i,k)=15.
+              IF (latgeo.GE.10.0 .and. latgeo.LT.20.0) srcbe(i,k)=11.5
+              IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=10.
+              IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=6.
+           END IF
+           IF ( k .EQ. 27 ) THEN
+              IF (latgeo.GE.60.0) srcbe(i,k)=46.
+              IF (latgeo.GE.50.0 .and. latgeo.LT.60.0) srcbe(i,k)=44.
+              IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=35.
+              IF (latgeo.GE.30.0 .and. latgeo.LT.40.0) srcbe(i,k)=25.
+              IF (latgeo.GE.20.0 .and. latgeo.LT.30.0) srcbe(i,k)=16.
+              IF (latgeo.GE.10.0 .and. latgeo.LT.20.0) srcbe(i,k)=12.5
+              IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=10.
+           END IF
+           IF ( k .EQ. 28 ) THEN
+              IF (latgeo.GE.60.0) srcbe(i,k)=53.
+              IF (latgeo.GE.50.0 .and. latgeo.LT.60.0) srcbe(i,k)=48.
+              IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=37.
+              IF (latgeo.GE.30.0 .and. latgeo.LT.40.0) srcbe(i,k)=25.
+              IF (latgeo.GE.20.0 .and. latgeo.LT.30.0) srcbe(i,k)=16.
+              IF (latgeo.GE.10.0 .and. latgeo.LT.20.0) srcbe(i,k)=12.5
+              IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=10.
+           END IF
+           IF ( k .EQ. 29 ) THEN
+              IF (latgeo.GE.70.0) srcbe(i,k)=58.
+              IF (latgeo.GE.60.0 .and. latgeo.LT.70.0) srcbe(i,k)=56.
+              IF (latgeo.GE.50.0 .and. latgeo.LT.60.0) srcbe(i,k)=50.
+              IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=36.
+              IF (latgeo.GE.30.0 .and. latgeo.LT.40.0) srcbe(i,k)=24.
+              IF (latgeo.GE.20.0 .and. latgeo.LT.30.0) srcbe(i,k)=15.
+              IF (latgeo.GE.10.0 .and. latgeo.LT.20.0) srcbe(i,k)=11.5
+              IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=10.
+           END IF
+           IF ( k .EQ. 30 ) THEN
+              IF (latgeo.GE.70.0) srcbe(i,k)=65.
+              IF (latgeo.GE.60.0 .and. latgeo.LT.70.0) srcbe(i,k)=60.
+              IF (latgeo.GE.50.0 .and. latgeo.LT.60.0) srcbe(i,k)=50.
+              IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=35.
+              IF (latgeo.GE.30.0 .and. latgeo.LT.40.0) srcbe(i,k)=22.
+              IF (latgeo.GE.20.0 .and. latgeo.LT.30.0) srcbe(i,k)=14.
+              IF (latgeo.GE.10.0 .and. latgeo.LT.20.0) srcbe(i,k)=10.
+              IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=9.
+           END IF
+           IF ( k .EQ. 31 ) THEN
+              IF (latgeo.GE.70.0) srcbe(i,k)=70.
+              IF (latgeo.GE.60.0 .and. latgeo.LT.70.0) srcbe(i,k)=62.
+              IF (latgeo.GE.50.0 .and. latgeo.LT.60.0) srcbe(i,k)=48.
+              IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=32.
+              IF (latgeo.GE.30.0 .and. latgeo.LT.40.0) srcbe(i,k)=21.
+              IF (latgeo.GE.20.0 .and. latgeo.LT.30.0) srcbe(i,k)=13.
+              IF (latgeo.GE.10.0 .and. latgeo.LT.20.0) srcbe(i,k)=9.
+              IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=7.6
+           END IF
+           IF ( k .EQ. 32 ) THEN
+              IF (latgeo.GE.70.0) srcbe(i,k)=80.
+              IF (latgeo.GE.60.0 .and. latgeo.LT.70.0) srcbe(i,k)=60.
+              IF (latgeo.GE.50.0 .and. latgeo.LT.60.0) srcbe(i,k)=46.
+              IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=30.
+              IF (latgeo.GE.30.0 .and. latgeo.LT.40.0) srcbe(i,k)=17.5
+              IF (latgeo.GE.20.0 .and. latgeo.LT.30.0) srcbe(i,k)=11.
+              IF (latgeo.GE.10.0 .and. latgeo.LT.20.0) srcbe(i,k)=8.
+              IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=7.4
+           END IF
+           IF ( k .GE. 33 ) THEN
+              IF (latgeo.GE.70.0) srcbe(i,k)=80.
+              IF (latgeo.GE.60.0 .and. latgeo.LT.70.0) srcbe(i,k)=70.
+              IF (latgeo.GE.50.0 .and. latgeo.LT.60.0) srcbe(i,k)=45.
+              IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=27.
+              IF (latgeo.GE.30.0 .and. latgeo.LT.40.0) srcbe(i,k)=15.
+              IF (latgeo.GE.20.0 .and. latgeo.LT.30.0) srcbe(i,k)=10.
+              IF (latgeo.GE.10.0 .and. latgeo.LT.20.0) srcbe(i,k)=7.6
+              IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=7.
+           END IF
+        END DO
+     END DO
+  END IF ! fin de 39 niveaux verticaux
+
+
+!====================================
+! Conversion de la source en U/s/kgA
+!====================================
+  DO k = 1,klev
+     DO i = 1,klon
+       ! La source est  at/min/m3 -> at/s/kgA
+       ! avec une masse volumique de l'air = 1.295 kg/m3
+       ! 1/(60*1.295) = 0.01287
+       srcbe(i,k)=srcbe(i,k)*0.01287
+       ! La source est  at/min/m3 -> at/s/m3
+       ! srcbe(i,k)=srcbe(i,k)*0.0166667
+    END DO
+ END DO
+
+END SUBROUTINE init_be
Index: DZ4/branches/LMDZ4-dev/libf/phylmd/initrrnpb.F
===================================================================
--- /LMDZ4/branches/LMDZ4-dev/libf/phylmd/initrrnpb.F	(revision 1190)
+++ 	(revision )
@@ -1,125 +1,0 @@
-!
-! $Header$
-!
-      SUBROUTINE  initrrnpb(ftsol,pctsrf,masktr,fshtr,hsoltr,tautr
-     .                   ,vdeptr,scavtr)
-      USE dimphy
-      USE infotrac, ONLY : nbtr
-      IMPLICIT none
-c======================================================================
-c Auteur(s): AA + CG (LGGE/CNRS) Date 24-06-94
-c Objet: initialisation des constantes des traceurs
-CAA Revison pour le controle avec la temperature du sol
-cAA
-CAA   it = 1 radon ss controle de ts
-cAA   it = 2 plomb ss controle de ts  
-c======================================================================
-c Arguments:
-c nbtr------input-I- nombre de vrais traceurs (sans l'eau)
-c ftsol-------input-R- Temperature du sol (Kelvin)
-c pctsrf-----input-R-  Nature de sol (pourcentage de sol)
-c masktr---output-R- Masque reservoir de sol traceur (1 = reservoir)
-c fshtr----output-R- Flux surfacique de production dans le sol
-c hsoltr---output-R- Epaisseur du reservoir de sol
-c tautr----output-R- Constante de decroissance du traceur
-c vdeptr---output-R- Vitesse de depot sec dans la couche Brownienne
-c scavtr---output-R- Coefficient de lessivage
-c======================================================================
-cym#include "dimensions.h"
-cym#include "dimphy.h"
-#include "indicesol.h"
-c======================================================================
-C
-      INTEGER i, it
-      REAL pctsrf(klon,nbsrf) !Pourcentage de sol (f(nature du sol))
-      REAL ftsol(klon,nbsrf)  ! Temperature du sol pour le controle Rn
-c                             ! le cas echeant
-      REAL masktr(klon,nbtr)  ! Masque de l'echange avec la surface
-c                                 (possible => 1 )
-      REAL fshtr(klon,nbtr)  ! Flux surfacique dans le reservoir de sol
-      REAL hsoltr(nbtr)      ! Epaisseur equivalente du reservoir de sol
-      REAL tautr(nbtr)       ! Constante de decroissance radioactive
-      REAL vdeptr(nbtr)      ! Vitesse de depot sec dans la couche Brownienne
-      REAL scavtr(nbtr)      ! Coefficient de lessivage
-      REAL s
-C
-      WRITE(*,'(''PASSAGE initrrnpb ...'',$)')
-      print*,'nbtr= ',nbtr
-      print*,'nbsrf= ',nbsrf
-      print*,'klon= ',klon
-C
-C Puis les initialisation specifiques a chaque traceur (pour le moment, Rn222)
-C
-C
-C Radon it = 1
-c
-      IF ( nbtr .LE. 0 ) STOP 'initrrnpb pas glop pas glop' 
-      it = 1
-      s = 1.E4  !  Source: atome par m2
-      hsoltr(it) = 0.1      ! Hauteur equivalente du reservoir : 
-c                              1 m * porosite 0.1
-      tautr(it) = 4.765E5  ! Decroissance du radon, secondes
-cAA
-c      tautr(it) = 4.765E55  ! Decroissance du radon,infinie
-cAA
-      vdeptr(it) = 0. ! Pas de depot sec pour le radon
-      scavtr(it) = 0. ! Pas de lessivage pour le radon
-
-      print*, '-------------- SOURCE DU RADON ------------------------ '
-      print*,'it = ',it
-      print*,'Source : ', s
-      print*,'Hauteur equivalente du reservoir de sol: ',hsoltr(it) 
-      print*,'Decroissance (s): ', tautr(it)
-      print*,'Vitesse de depot sec: ',vdeptr(it) 
-      print*,'Facteur de lessivage: ',scavtr(it)
-
-      DO i = 1,klon
-        masktr(i,it) = 0.
-        IF ( NINT(pctsrf(i,1)) .EQ. 1 ) masktr(i,it) = 1.
-        fshtr(i,it) = s * masktr(i,it)
-
-cAA
-cAA quelques tests
-cAA POur l'instant le pctsrf(i,3) = 1.0 
-cAA lorsqu'il ya de la terre mias ne prend aucune autre valeur
-cAA il n'est donc pas necessaire de multiplier fshtr par pctsrf
-cAA 
-c       print*, '------------------------------------------ '
-c        print*, 'masktr(',i,it,')= ',masktr(i,it)
-c        print*, 'fshtr(',i,it,')= ',fshtr(i,it)
-c        print*, 'pctsrf(',i,',1)= ',pctsrf(i,1)
-c        print*, 'pctsrf(',i,',2)= ',pctsrf(i,2)
-c        print*, 'pctsrf(',i,',3)= ',pctsrf(i,3)
-c        print*, 'pctsrf(',i,',4)= ',pctsrf(i,4)
-c        print*, 's = ',s
-c        print*, '------------------------------------------ '
-
-      END DO
-C
-C 210Pb it = 2
-C
-      IF ( nbtr .LE. 1 ) STOP 'initrrnpb pas glop pas glop' 
-      it = 2
-      s = 0. !  Pas de source !!!
-      hsoltr(it) = 10.     ! Hauteur equivalente du reservoir 
-c                              a partir duquel le
-c                              depot Brownien a lieu
-      tautr(it) = 1.028E9 ! Decroissance du Pb210, secondes
-      vdeptr(it) = 1.E-3 ! 1 mm/s pour le 210Pb
-      scavtr(it) =  .5   ! Lessivage du Pb210
-      DO i = 1,klon
-        masktr(i,it) = 1. ! Le depot sec peut avoir lieu partout
-        fshtr(i,it) = s * masktr(i,it)
-      END DO
-      print*, '-------------- SOURCE DU PLOMB ------------------------ '
-      print*,'it = ',it
-      print*,'Source : ', s
-      print*,'Hauteur equivalente du reservoir : ',hsoltr(it) 
-      print*,'Decroissance (s): ', tautr(it)
-      print*,'Vitesse de depot sec: ',vdeptr(it) 
-      print*,'Facteur de lessivage: ',scavtr(it)
-c
-      WRITE(*,*) 'initialisation rnpb ok'
-c
-      RETURN
-      END
Index: /LMDZ4/branches/LMDZ4-dev/libf/phylmd/initrrnpb.F90
===================================================================
--- /LMDZ4/branches/LMDZ4-dev/libf/phylmd/initrrnpb.F90	(revision 1191)
+++ /LMDZ4/branches/LMDZ4-dev/libf/phylmd/initrrnpb.F90	(revision 1191)
@@ -0,0 +1,92 @@
+!
+! $Id $
+!
+SUBROUTINE  initrrnpb(ftsol,pctsrf,masktr,fshtr,hsoltr,tautr,vdeptr,scavtr)
+  USE dimphy
+  USE infotrac, ONLY : nbtr
+  IMPLICIT NONE
+!======================================================================
+! Auteur(s): AA + CG (LGGE/CNRS) Date 24-06-94
+! Objet: initialisation des constantes des traceurs
+!AA Revison pour le controle avec la temperature du sol
+!AA
+!AA   it = 1 radon ss controle de ts
+!AA   it = 2 plomb ss controle de ts  
+!======================================================================
+! Arguments:
+! nbtr.............. nombre de vrais traceurs (sans l'eau)
+! ftsol....input-R-  Temperature du sol (Kelvin)
+! pctsrf...input-R-  Nature de sol (pourcentage de sol)
+! masktr...output-R- Masque reservoir de sol traceur (1 = reservoir)
+! fshtr....output-R- Flux surfacique de production dans le reservoir de sol
+! hsoltr...output-R- Epaisseur equivalente du reservoir de sol
+! tautr....output-R- Constante de decroissance radioactive du traceur
+! vdeptr...output-R- Vitesse de depot sec dans la couche Brownienne
+! scavtr...output-R- Coefficient de lessivage
+!======================================================================
+  INCLUDE "indicesol.h"
+!======================================================================
+
+  REAL,DIMENSION(klon,nbsrf),INTENT(IN) :: pctsrf
+  REAL,DIMENSION(klon,nbsrf),INTENT(IN) :: ftsol
+  REAL,DIMENSION(klon,nbtr),INTENT(OUT) :: masktr
+  REAL,DIMENSION(klon,nbtr),INTENT(OUT) :: fshtr
+  REAL,DIMENSION(nbtr),INTENT(OUT)      :: hsoltr
+  REAL,DIMENSION(nbtr),INTENT(OUT)      :: tautr
+  REAL,DIMENSION(nbtr),INTENT(OUT)      :: vdeptr
+  REAL,DIMENSION(nbtr),INTENT(OUT)      :: scavtr
+  INTEGER                               :: i, it
+  REAL                                  :: s
+
+  WRITE(*,*)'PASSAGE initrrnpb ...'
+!
+! Radon it = 1
+!----------------
+  IF ( nbtr .LE. 0 ) STOP '**PHYTRAC:initrrnpb:** nbtr < 0; verifier RN dans traceur.def' 
+  it = 1
+  s = 1.E4             ! Source: atome par m2
+  hsoltr(it) = 0.1     ! Hauteur equivalente du reservoir : 
+                       ! 1 m * porosite 0.1
+  tautr(it) = 4.765E5  ! Decroissance du radon, secondes
+  vdeptr(it) = 0.      ! Pas de depot sec pour le radon
+  scavtr(it) = 0.      ! Pas de lessivage pour le radon
+  
+  WRITE(*,*)'-------------- SOURCE DU RADON ------------------------ '
+  WRITE(*,*)'it = ',it
+  WRITE(*,*)'Source : ', s
+  WRITE(*,*)'Hauteur equivalente du reservoir de sol: ',hsoltr(it) 
+  WRITE(*,*)'Decroissance (s): ', tautr(it)
+  WRITE(*,*)'Vitesse de depot sec: ',vdeptr(it) 
+  WRITE(*,*)'Facteur de lessivage: ',scavtr(it)
+
+  DO i = 1,klon
+     masktr(i,it) = 0.
+     IF ( NINT(pctsrf(i,1)) .EQ. 1 ) masktr(i,it) = 1.
+     fshtr(i,it) = s * masktr(i,it)
+  END DO
+!
+! 210Pb it = 2
+!----------------
+  IF ( nbtr .LE. 1 ) STOP '**PHYTRAC**:initrrnpb:** nbtr <= 1; verifier PB dans traceur.def' 
+  it = 2
+  s = 0.                ! Pas de source 
+  hsoltr(it) = 10.      ! Hauteur equivalente du reservoir 
+                        ! a partir duquel le depot Brownien a lieu
+  tautr(it) = 1.028E9   ! Decroissance du Pb210, secondes
+  vdeptr(it) = 1.E-3    ! 1 mm/s pour le 210Pb
+  scavtr(it) =  .5      ! Lessivage du Pb210
+  DO i = 1,klon
+     masktr(i,it) = 1.  ! Le depot sec peut avoir lieu partout
+     fshtr(i,it) = s * masktr(i,it)
+  END DO
+  WRITE(*,*)'-------------- SOURCE DU PLOMB ------------------------ '
+  WRITE(*,*)'it = ',it
+  WRITE(*,*)'Source : ', s
+  WRITE(*,*)'Hauteur equivalente du reservoir : ',hsoltr(it) 
+  WRITE(*,*)'Decroissance (s): ', tautr(it)
+  WRITE(*,*)'Vitesse de depot sec: ',vdeptr(it) 
+  WRITE(*,*)'Facteur de lessivage: ',scavtr(it)
+
+  WRITE(*,*) 'Initialisation RN et PB ok'
+
+END SUBROUTINE initrrnpb
Index: DZ4/branches/LMDZ4-dev/libf/phylmd/minmaxqfi.F
===================================================================
--- /LMDZ4/branches/LMDZ4-dev/libf/phylmd/minmaxqfi.F	(revision 1190)
+++ 	(revision )
@@ -1,34 +1,0 @@
-!
-! $Header$
-!
-      SUBROUTINE minmaxqfi(zq,qmin,qmax,comment)
-      USE dimphy
-      IMPLICIT none
-
-cym#include "dimensions.h"
-cym#include "dimphy.h"
-
-      CHARACTER*(*) comment
-      real qmin,qmax
-      real zq(klon,klev)
-
-      INTEGER jadrs(klon), jbad, k, i
-
-      DO k = 1, klev
-         jbad = 0
-         DO i = 1, klon
-         IF (zq(i,k).GT.qmax .OR. zq(i,k).LT.qmin) THEN
-            jbad = jbad + 1
-            jadrs(jbad) = i
-         ENDIF
-         ENDDO
-         IF (jbad.GT.0) THEN
-         PRINT*, comment
-         DO i = 1, jbad
-            PRINT*, "i,k,q=", jadrs(i),k,zq(jadrs(i),k)
-         ENDDO
-         ENDIF
-      ENDDO
-
-      return
-      end
Index: /LMDZ4/branches/LMDZ4-dev/libf/phylmd/minmaxqfi.F90
===================================================================
--- /LMDZ4/branches/LMDZ4-dev/libf/phylmd/minmaxqfi.F90	(revision 1191)
+++ /LMDZ4/branches/LMDZ4-dev/libf/phylmd/minmaxqfi.F90	(revision 1191)
@@ -0,0 +1,33 @@
+!
+! $Id$
+!
+SUBROUTINE minmaxqfi(zq,qmin,qmax,comment)
+  USE dimphy
+  IMPLICIT NONE
+
+! Entrees
+  REAL,DIMENSION(klon,klev), INTENT(IN)   :: zq
+  REAL,INTENT(IN)                         :: qmin,qmax
+  CHARACTER(LEN=*),INTENT(IN)             :: comment
+
+! Local  
+  INTEGER,DIMENSION(klon)     :: jadrs 
+  INTEGER                     :: i, jbad, k
+  
+  DO k = 1, klev
+     jbad = 0
+     DO i = 1, klon
+        IF (zq(i,k).GT.qmax .OR. zq(i,k).LT.qmin) THEN
+           jbad = jbad + 1
+           jadrs(jbad) = i
+        ENDIF
+     ENDDO
+     IF (jbad.GT.0) THEN
+        WRITE(*,*)comment
+        DO i = 1, jbad
+           WRITE(*,*) "i,k,q=", jadrs(i),k,zq(jadrs(i),k)
+        ENDDO
+     ENDIF
+  ENDDO
+  
+END SUBROUTINE minmaxqfi
Index: DZ4/branches/LMDZ4-dev/libf/phylmd/nflxtr.F
===================================================================
--- /LMDZ4/branches/LMDZ4-dev/libf/phylmd/nflxtr.F	(revision 1190)
+++ 	(revision )
@@ -1,166 +1,0 @@
-!
-! $Header$
-!
-      SUBROUTINE nflxtr(pdtime,pmfu,pmfd,pen_u,pde_u,pen_d,pde_d,
-     .                 pplay,paprs,x,dx) 
-      USE dimphy
-      IMPLICIT NONE 
-c=====================================================================
-c Objet : Melange convectif de traceurs a partir des flux de masse 
-c Date : 13/12/1996 -- 13/01/97
-c Auteur: O. Boucher (LOA) sur inspiration de Z. X. Li (LMD),
-c         Brinkop et Sausen (1996) et Boucher et al. (1996).
-c ATTENTION : meme si cette routine se veut la plus generale possible, 
-c             elle a herite de certaines notations et conventions du 
-c             schema de Tiedtke (1993). 
-c --En particulier, les couches sont numerotees de haut en bas !!!
-c   Ceci est valable pour les flux
-c   mais pas pour les entrees x, pplay, paprs !!!!
-c --pmfu est positif, pmfd est negatif 
-c --Tous les flux d'entrainements et de detrainements sont positifs 
-c   contrairement au schema de Tiedtke d'ou les changements de signe!!!! 
-c=====================================================================
-c
-cym#include "dimensions.h"
-cym#include "dimphy.h"
-#include "YOMCST.h"
-#include "YOECUMF.h" 
-c
-      REAL pdtime
-c--les flux sont definis au 1/2 niveaux
-c--pmfu(klev+1) et pmfd(klev+1) sont implicitement nuls
-      REAL pmfu(klon,klev)  ! flux de masse dans le panache montant 
-      REAL pmfd(klon,klev)  ! flux de masse dans le panache descendant
-      REAL pen_u(klon,klev) ! flux entraine dans le panache montant
-      REAL pde_u(klon,klev) ! flux detraine dans le panache montant
-      REAL pen_d(klon,klev) ! flux entraine dans le panache descendant
-      REAL pde_d(klon,klev) ! flux detraine dans le panache descendant
-
-      REAL pplay(klon,klev)    ! pression aux couches (bas en haut)
-      REAL paprs(klon,klev+1)  ! pression aux 1/2 couches (bas en haut)
-      REAL x(klon,klev)        ! q de traceur (bas en haut) 
-      REAL dx(klon,klev)     ! tendance de traceur  (bas en haut)
-c
-c--flux convectifs mais en variables locales
-      REAL zmfu(klon,klev+1) 
-      REAL zmfd(klon,klev+1) 
-      REAL zen_u(klon,klev) 
-      REAL zde_u(klon,klev)
-      REAL zen_d(klon,klev) 
-      REAL zde_d(klon,klev)
-      real zmfe
-c
-c--variables locales      
-c--les flux de x sont definis aux 1/2 niveaux 
-c--xu et xd sont definis aux niveaux complets
-      REAL xu(klon,klev)        ! q de traceurs dans le panache montant
-      REAL xd(klon,klev)        ! q de traceurs dans le panache descendant
-      REAL zmfux(klon,klev+1)   ! flux de x dans le panache montant
-      REAL zmfdx(klon,klev+1)   ! flux de x dans le panache descendant
-      REAL zmfex(klon,klev+1)   ! flux de x dans l'environnement 
-      INTEGER i, k 
-      REAL zmfmin
-      PARAMETER (zmfmin=1.E-10)
-
-c =========================================
-c
-c
-c   Extension des flux UP et DN sur klev+1 niveaux
-c =========================================
-      do k=1,klev
-         do i=1,klon
-            zmfu(i,k)=pmfu(i,k)
-            zmfd(i,k)=pmfd(i,k)
-         enddo
-      enddo
-      do i=1,klon
-         zmfu(i,klev+1)=0.
-         zmfd(i,klev+1)=0.
-      enddo
-
-c--modif pour diagnostiquer les detrainements
-c =========================================
-c   on privilegie l'ajustement de l'entrainement dans l'ascendance.
-
-      do k=1, klev
-         do i=1, klon
-            zen_d(i,k)=pen_d(i,k)
-            zde_u(i,k)=pde_u(i,k)
-            zde_d(i,k) =-zmfd(i,k+1)+zmfd(i,k)+zen_d(i,k)
-            zen_u(i,k) = zmfu(i,k+1)-zmfu(i,k)+zde_u(i,k)
-         enddo 
-      enddo 
-c
-c--calcul des flux dans le panache montant
-c =========================================
-c
-c Dans la premiere couche, on prend q comme valeur de qu
-c
-      do i=1, klon
-         zmfux(i,1)=0.0 
-      enddo
-c
-c Autres couches
-      do k=1,klev
-         do i=1, klon
-            if ((zmfu(i,k+1)+zde_u(i,k)).lt.zmfmin) THEN
-               xu(i,k)=x(i,k)
-            else
-               xu(i,k)=(zmfux(i,k)+zen_u(i,k)*x(i,k))
-     s               /(zmfu(i,k+1)+zde_u(i,k))
-            endif
-            zmfux(i,k+1)=zmfu(i,k+1)*xu(i,k)
-         enddo
-      enddo
-c
-c--calcul des flux dans le panache descendant
-c =========================================
-c   
-      do i=1, klon
-         zmfdx(i,klev+1)=0.0 
-      enddo
-c
-      do k=klev,1,-1
-         do i=1, klon
-            if ((zde_d(i,k)-zmfd(i,k)).lt.zmfmin) THEN
-               xd(i,k)=x(i,k)
-            else
-               xd(i,k)=(zmfdx(i,k+1)-zen_d(i,k)*x(i,k)) /
-     .               (zmfd(i,k)-zde_d(i,k))
-            endif
-            zmfdx(i,k)=zmfd(i,k)*xd(i,k)
-         enddo
-      enddo
-c
-c--introduction du flux de retour dans l'environnement
-c =========================================
-c
-      do k=2, klev
-         do i=1, klon
-            zmfe=-zmfu(i,k)-zmfd(i,k)
-            if (zmfe.le.0.) then
-               zmfex(i,k)= zmfe*x(i,k)
-            else
-               zmfex(i,k)= zmfe*x(i,k-1)
-            endif
-         enddo
-      enddo
-
-      do i=1, klon
-         zmfex(i,1)=0.
-         zmfex(i,klev+1)=0.
-      enddo
-c
-c--calcul final des tendances
-c
-      do k=1, klev
-         do i=1, klon
-            dx(i,k)=RG/(paprs(i,k)-paprs(i,k+1))*pdtime*
-     .                      ( zmfux(i,k) - zmfux(i,k+1) +
-     .                        zmfdx(i,k) - zmfdx(i,k+1) +
-     .                        zmfex(i,k) - zmfex(i,k+1) )
-         enddo
-      enddo
-c
-      return 
-      end
Index: /LMDZ4/branches/LMDZ4-dev/libf/phylmd/nflxtr.F90
===================================================================
--- /LMDZ4/branches/LMDZ4-dev/libf/phylmd/nflxtr.F90	(revision 1191)
+++ /LMDZ4/branches/LMDZ4-dev/libf/phylmd/nflxtr.F90	(revision 1191)
@@ -0,0 +1,159 @@
+!
+! $Id $
+!
+SUBROUTINE nflxtr(pdtime,pmfu,pmfd,pen_u,pde_u,pen_d,pde_d,pplay,paprs,x,dx) 
+  USE dimphy
+  IMPLICIT NONE 
+!=====================================================================
+! Objet : Melange convectif de traceurs a partir des flux de masse 
+! Date : 13/12/1996 -- 13/01/97
+! Auteur: O. Boucher (LOA) sur inspiration de Z. X. Li (LMD),
+!         Brinkop et Sausen (1996) et Boucher et al. (1996).
+! ATTENTION : meme si cette routine se veut la plus generale possible, 
+!             elle a herite de certaines notations et conventions du 
+!             schema de Tiedtke (1993).
+! 1. En particulier, les couches sont numerotees de haut en bas !!!
+!    Ceci est valable pour les flux
+!    mais pas pour les entrees x, pplay, paprs !!!!
+! 2. pmfu est positif, pmfd est negatif 
+! 3. Tous les flux d'entrainements et de detrainements sont positifs 
+!    contrairement au schema de Tiedtke d'ou les changements de signe!!!! 
+!=====================================================================
+!
+  include "YOMCST.h"
+  include "YOECUMF.h" 
+
+  REAL,INTENT(IN) :: pdtime  ! pdtphys
+!
+! les flux sont definis au 1/2 niveaux 
+! => pmfu(klev+1) et pmfd(klev+1) sont implicitement nuls
+!
+  REAL,DIMENSION(klon,klev),INTENT(IN) :: pmfu  ! flux de masse dans le panache montant 
+  REAL,DIMENSION(klon,klev),INTENT(IN) :: pmfd  ! flux de masse dans le panache descendant
+  REAL,DIMENSION(klon,klev),INTENT(IN) :: pen_u ! flux entraine dans le panache montant
+  REAL,DIMENSION(klon,klev),INTENT(IN) :: pde_u ! flux detraine dans le panache montant
+  REAL,DIMENSION(klon,klev),INTENT(IN) :: pen_d ! flux entraine dans le panache descendant
+  REAL,DIMENSION(klon,klev),INTENT(IN) :: pde_d ! flux detraine dans le panache descendant
+
+  REAL,DIMENSION(klon,klev),INTENT(IN)   :: pplay ! pression aux couches (bas en haut)
+  REAL,DIMENSION(klon,klev+1),INTENT(IN) :: paprs ! pression aux 1/2 couches (bas en haut)
+  REAL,DIMENSION(klon,klev),INTENT(IN)   :: x     ! q de traceur (bas en haut) 
+  REAL,DIMENSION(klon,klev),INTENT(INOUT) :: dx   ! tendance de traceur  (bas en haut)
+
+! flux convectifs mais en variables locales
+  REAL,DIMENSION(klon,klev+1) :: zmfu  ! copie de pmfu avec klev+1 = 0
+  REAL,DIMENSION(klon,klev+1) :: zmfd  ! copie de pmfd avec klev+1 = 0
+  REAL,DIMENSION(klon,klev)   :: zen_u
+  REAL,DIMENSION(klon,klev)   :: zde_u
+  REAL,DIMENSION(klon,klev)   :: zen_d
+  REAL,DIMENSION(klon,klev)   :: zde_d
+  REAL                        :: zmfe
+
+! variables locales      
+! les flux de x sont definis aux 1/2 niveaux 
+! xu et xd sont definis aux niveaux complets
+  REAL,DIMENSION(klon,klev)   :: xu      ! q de traceurs dans le panache montant
+  REAL,DIMENSION(klon,klev)   :: xd      ! q de traceurs dans le panache descendant
+  REAL,DIMENSION(klon,klev+1) :: zmfux   ! flux de x dans le panache montant
+  REAL,DIMENSION(klon,klev+1) :: zmfdx   ! flux de x dans le panache descendant
+  REAL,DIMENSION(klon,klev+1) :: zmfex   ! flux de x dans l'environnement 
+  INTEGER                     :: i, k 
+  REAL,PARAMETER              :: zmfmin=1.E-10
+
+! ==============================================
+! Extension des flux UP et DN sur klev+1 niveaux
+! ==============================================
+  DO k=1,klev
+     DO i=1,klon
+        zmfu(i,k)=pmfu(i,k)
+        zmfd(i,k)=pmfd(i,k)
+     ENDDO
+  ENDDO
+  DO i=1,klon
+     zmfu(i,klev+1)=0.
+     zmfd(i,klev+1)=0.
+  ENDDO
+! ==========================================
+! modif pour diagnostiquer les detrainements
+! ==========================================
+!   on privilegie l'ajustement de l'entrainement dans l'ascendance.
+
+  DO k=1, klev
+     DO i=1, klon
+        zen_d(i,k)=pen_d(i,k)
+        zde_u(i,k)=pde_u(i,k)
+        zde_d(i,k) =-zmfd(i,k+1)+zmfd(i,k)+zen_d(i,k)
+        zen_u(i,k) = zmfu(i,k+1)-zmfu(i,k)+zde_u(i,k)
+     ENDDO
+  ENDDO
+! =========================================
+! calcul des flux dans le panache montant
+! =========================================
+!
+! Dans la premiere couche, on prend q comme valeur de qu
+
+  DO i=1, klon
+     zmfux(i,1)=0.0 
+  ENDDO
+
+! Autres couches
+  DO k=1,klev
+     DO i=1, klon
+        IF ((zmfu(i,k+1)+zde_u(i,k)).lt.zmfmin) THEN
+           xu(i,k)=x(i,k)
+        ELSE
+           xu(i,k)=(zmfux(i,k)+zen_u(i,k)*x(i,k))/(zmfu(i,k+1)+zde_u(i,k))
+        ENDIF
+        zmfux(i,k+1)=zmfu(i,k+1)*xu(i,k)
+     ENDDO
+  ENDDO
+! ==========================================
+! calcul des flux dans le panache descendant
+! ==========================================
+   
+  DO i=1, klon
+     zmfdx(i,klev+1)=0.0 
+  ENDDO
+
+  DO k=klev,1,-1
+     DO i=1, klon
+        IF ((zde_d(i,k)-zmfd(i,k)).lt.zmfmin) THEN
+           xd(i,k)=x(i,k)
+        ELSE
+           xd(i,k)=(zmfdx(i,k+1)-zen_d(i,k)*x(i,k))/(zmfd(i,k)-zde_d(i,k))
+        ENDIF
+        zmfdx(i,k)=zmfd(i,k)*xd(i,k)
+     ENDDO
+  ENDDO
+! ===================================================
+! introduction du flux de retour dans l'environnement
+! ===================================================
+
+  DO k=2, klev
+     DO i=1, klon
+        zmfe=-zmfu(i,k)-zmfd(i,k)
+        IF (zmfe.le.0.) then
+           zmfex(i,k)= zmfe*x(i,k)
+        ELSE
+           zmfex(i,k)= zmfe*x(i,k-1)
+        ENDIF
+     ENDDO
+  ENDDO
+
+  DO i=1, klon
+     zmfex(i,1)=0.
+     zmfex(i,klev+1)=0.
+  ENDDO
+! ==========================
+! calcul final des tendances
+! ==========================
+  DO k=1, klev
+     DO i=1, klon
+        dx(i,k)=RG/(paprs(i,k)-paprs(i,k+1))*pdtime*  &
+             ( zmfux(i,k) - zmfux(i,k+1) +            &
+             zmfdx(i,k) - zmfdx(i,k+1) +              &
+             zmfex(i,k) - zmfex(i,k+1) )
+     ENDDO
+  ENDDO
+  
+END SUBROUTINE nflxtr
Index: /LMDZ4/branches/LMDZ4-dev/libf/phylmd/phyetat0.F
===================================================================
--- /LMDZ4/branches/LMDZ4-dev/libf/phylmd/phyetat0.F	(revision 1190)
+++ /LMDZ4/branches/LMDZ4-dev/libf/phylmd/phyetat0.F	(revision 1191)
@@ -19,4 +19,7 @@
       USE iostart
       USE write_field_phy
+      USE infotrac
+      USE traclmdz_mod,    ONLY : traclmdz_from_restart
+
       IMPLICIT none
 c======================================================================
@@ -48,4 +51,5 @@
       REAL run_off_lic_0(klon)
       REAL fractint(klon)
+      REAL trs(klon,nbtr)
 
       CHARACTER*6 ocean_in
@@ -62,4 +66,5 @@
       INTEGER length
       PARAMETER (length=100)
+      INTEGER it, iiq
       REAL tab_cntrl(length), tabcntr0(length)
       CHARACTER*7 str7
@@ -983,4 +988,27 @@
       PRINT*,'(ecart-type) wake_fip:', xmin, xmax
 c
+c Read and send field trs to traclmdz
+c
+      IF (type_trac == 'lmdz') THEN
+         DO it=1,nbtr
+            iiq=niadv(it+2)
+            CALL get_field("trs_"//tname(iiq),trs(:,it),found)
+            IF (.NOT. found) THEN
+               PRINT*, 
+     $           "phyetat0: Le champ <trs_"//tname(iiq)//"> est absent"
+               PRINT*, "Depart legerement fausse. Mais je continue"
+               trs(:,it) = 0.
+            ENDIF
+            xmin = 1.0E+20
+            xmax = -1.0E+20
+            xmin = MINval(trs(:,it))
+            xmax = MAXval(trs(:,it))
+            PRINT*,"(ecart-type) trs_"//tname(iiq)//" :", xmin, xmax
+
+         END DO
+         
+         CALL traclmdz_from_restart(trs)
+      END IF
+
 
 c on ferme le fichier
@@ -1005,4 +1033,5 @@
       CALL fonte_neige_init(run_off_lic_0)
 
+
       RETURN
       END
Index: /LMDZ4/branches/LMDZ4-dev/libf/phylmd/phyredem.F
===================================================================
--- /LMDZ4/branches/LMDZ4-dev/libf/phylmd/phyredem.F	(revision 1190)
+++ /LMDZ4/branches/LMDZ4-dev/libf/phylmd/phyredem.F	(revision 1191)
@@ -12,4 +12,6 @@
       USE phys_state_var_mod
       USE iostart
+      USE traclmdz_mod, ONLY : traclmdz_to_restart
+      USE infotrac
 
       IMPLICIT none
@@ -42,4 +44,5 @@
       REAL agesno(klon,nbsrf)
       REAL run_off_lic_0(klon)
+      REAL trs(klon,nbtr)
 c
       INTEGER nid, nvarid, idim1, idim2, idim3
@@ -52,5 +55,6 @@
       CHARACTER (len=7) :: str7
       CHARACTER (len=2) :: str2
-
+      INTEGER           :: it, iiq
+      
 c======================================================================
 c 
@@ -311,4 +315,14 @@
       CALL put_field("WAKE_FIP","",wake_fip)
 
+
+! trs from traclmdz_mod
+      IF (type_trac == 'lmdz') THEN
+         CALL traclmdz_to_restart(trs)
+         DO it=1,nbtr
+            iiq=niadv(it+2)
+            CALL put_field("trs_"//tname(iiq),"",trs(:,it))
+         END DO
+      END IF
+
       CALL close_restartphy
 !$OMP BARRIER
Index: /LMDZ4/branches/LMDZ4-dev/libf/phylmd/physiq.F
===================================================================
--- /LMDZ4/branches/LMDZ4-dev/libf/phylmd/physiq.F	(revision 1190)
+++ /LMDZ4/branches/LMDZ4-dev/libf/phylmd/physiq.F	(revision 1191)
@@ -111,6 +111,4 @@
       PARAMETER (ok_stratus=.FALSE.)
 c======================================================================
-      LOGICAL, SAVE :: rnpb=.TRUE.
-c$OMP THREADPRIVATE(rnpb)
       REAL amn, amx
       INTEGER igout
@@ -3133,67 +3131,22 @@
 c====================================================================
 C
-      IF (config_inca /= 'none') rnpb=.FALSE.
-
-      call phytrac (     rnpb,
-     I                   itap, 
-     I                   julien, 
-     I                   gmtime,
-     I                   debut,
-     I                   lafin,
-     I                   nlon,
-     I                   klev,
-     I                   dtime,
-     I                   u,
-     I                   v,
-     I                   t,
-     I                   paprs,
-     I                   pplay,
-     I                   pmfu, 
-     I                   pmfd, 
-     I                   pen_u, 
-     I                   pde_u, 
-     I                   pen_d, 
-     I                   pde_d,
-     I                   cdragh,
-     I                   coefh,
-     I                   fm_therm,
-     I                   entr_therm,
-     I                   u1,
-     I                   v1,
-     I                   ftsol,
-     I                   pctsrf,
-     I                   rlat,
-     I                   frac_impa, 
-     I                   frac_nucl,
-     I                   rlon,
-     I                   presnivs,
-     I                   pphis,
-     I                   pphi,
-     I                   albsol1,
-     I                   qx(1,1,1), 
-     I                   rhcl,
-     I                   cldfra, 
-     I                   rneb, 
-     I                   diafra, 
-     I                   cldliq, 
-     I                   itop_con,
-     I                   ibas_con,
-     I                   pmflxr,
-     I                   pmflxs,
-     I                   prfl,
-     I                   psfl,
-     I                   da,
-     I                   phi,
-     I                   mp,
-     I                   upwd,
-     I                   dnwd,
-     I                   aerosol_couple,
-     I                   flxmass_w,
-     I                   tau_aero, 
-     I                   piz_aero, 
-     I                   cg_aero,
-     I                   ccm,
-     I                   rfname,
-     O                   tr_seri)
+
+      call phytrac (
+     I     itap,     julien,    gmtime,   debut,
+     I     lafin,    dtime,     u, v,     t,
+     I     paprs,    pplay,     pmfu,     pmfd, 
+     I     pen_u,    pde_u,     pen_d,    pde_d,
+     I     cdragh,   coefh,     fm_therm, entr_therm,
+     I     u1,       v1,        ftsol,    pctsrf,
+     I     rlat,     frac_impa, frac_nucl,rlon,
+     I     presnivs, pphis,     pphi,     albsol1,
+     I     qx(:,:,ivap),rhcl,   cldfra,   rneb, 
+     I     diafra,   cldliq,    itop_con, ibas_con,
+     I     pmflxr,   pmflxs,    prfl,     psfl,
+     I     da,       phi,       mp,       upwd,     
+     I     dnwd,     aerosol_couple,      flxmass_w,
+     I     tau_aero, piz_aero,  cg_aero,  ccm,
+     I     rfname,
+     O     tr_seri)
 
       IF (offline) THEN
Index: DZ4/branches/LMDZ4-dev/libf/phylmd/phytrac.F
===================================================================
--- /LMDZ4/branches/LMDZ4-dev/libf/phylmd/phytrac.F	(revision 1190)
+++ 	(revision )
@@ -1,833 +1,0 @@
-!
-c
-c
-      SUBROUTINE phytrac (rnpb,
-     I                    nstep,
-     I                    julien,
-     I                    gmtime,
-     I                    debutphy,
-     I                    lafin,
-     I                    nlon,
-     I                    nlev,
-     I                    pdtphys,
-     I                    u,
-     I                    v,
-     I                    t_seri,
-     I                    paprs,
-     I                    pplay,
-     I                    pmfu,
-     I                    pmfd,
-     I                    pen_u,
-     I                    pde_u,
-     I                    pen_d,
-     I                    pde_d,
-     I                    cdragh,
-     I                    coefh,
-     I                    fm_therm,
-     I                    entr_therm,
-     I                    yu1,
-     I                    yv1,
-     I                    ftsol,
-     I                    pctsrf,
-     I                    xlat,
-     I                    frac_impa,
-     I                    frac_nucl,
-     I                    xlon,
-     I                    presnivs,
-     I                    pphis,
-     I                    pphi,
-     I                    albsol,
-     I                    sh,
-     I                    rh,
-     I                    cldfra,
-     I                    rneb,
-     I                    diafra,
-     I                    cldliq,
-     I                    itop_con,
-     I                    ibas_con,
-     I                    pmflxr,
-     I                    pmflxs,
-     I                    prfl,
-     I                    psfl,
-     I                    da,
-     I                    phi,
-     I                    mp,
-     I                    upwd,
-     I                    dnwd,
-     I                    aerosol_couple,
-     I                    flxmass_w,
-     I                    tau_aero, 
-     I                    piz_aero, 
-     I                    cg_aero,
-     I                    ccm,
-     I                    rfname,
-     O                    tr_seri)
-
-      USE ioipsl
-      USE dimphy
-      USE infotrac
-      USE mod_grid_phy_lmdz
-      USE mod_phys_lmdz_para
-      USE comgeomphy
-      USE iophy
-      USE vampir
-
-      IMPLICIT none
-c======================================================================
-c Auteur(s) FH
-c Objet: Moniteur general des tendances traceurs
-c
-cAA Remarques en vrac:
-cAA--------------------
-cAA 2/ Le choix du radon et du pb se fait juste avec un data 
-cAA    (peu propre). Peut-etre pourrait-on prevoir dans l'avenir 
-cAA    une variable "type de traceur" 
-c======================================================================
-#include "YOMCST.h"
-#include "dimensions.h"
-#include "indicesol.h"
-#include "clesphys.h"
-#include "temps.h"
-#include "paramet.h"
-#include "control.h"
-#include "thermcell.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 nstep  ! appel physique
-      integer julien !jour julien
-      integer itop_con(nlon)
-      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,nbtr) ! traceur  
-      real u(klon,klev)
-      real v(klon,klev)
-      real sh(nlon,nlev)     ! humidite specifique
-      real rh(nlon,nlev)     ! humidite relative
-      real cldliq(nlon,nlev) ! eau liquide nuageuse
-      real cldfra(nlon,nlev) ! fraction nuageuse (tous les nuages)
-      real diafra(nlon,nlev) ! fraction nuageuse (convection ou stratus artificiels)
-      real rneb(nlon,nlev)   ! fraction nuageuse (grande echelle)
-      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,klev) ! geopotentiel
-      real pphis(klon)
-      REAL presnivs(klev)
-      logical debutphy       ! le flag de l'initialisation de la physique
-      logical lafin          ! le flag de la fin de la physique
-c Olivia      
-      integer nsplit
-      REAL pmflxr(klon,klev+1), pmflxs(klon,klev+1)   !--lessivage convection
-      REAL prfl(klon,klev+1),   psfl(klon,klev+1)     !--lessivage large-scale
-      LOGICAL aerosol_couple
-
-      REAL flxmass_w(klon,klev)
-      CHARACTER(len=8) :: solsym(nbtr)
-      integer la
-      REAL              ::    tau_aero(klon,klev,9,2)
-      REAL              ::    piz_aero(klon,klev,9,2)
-      REAL              ::    cg_aero(klon,klev,9,2)
-      character*4       ::    rfname(9) 
-      REAL              ::    ccm(klon,klev,2) 
-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   thermiques:
-c   -----------
-c
-      real fm_therm(klon,klev+1),entr_therm(klon,klev)
-	real fm_therm1(klon,klev)
-c
-      REAL pde_u(nlon,nlev) ! flux detraine dans le panache montant
-      REAL pen_d(nlon,nlev) ! flux entraine dans le panache descendant
-      REAL pde_d(nlon,nlev) ! flux detraine dans le panache descendant
-c KE
-      real da(nlon,nlev),phi(nlon,nlev,nlev),mp(nlon,nlev)
-      REAL upwd(nlon,nlev)      ! saturated updraft mass flux
-      REAL dnwd(nlon,nlev)      ! saturated downdraft mass flux
-
-c
-c   Couche limite:
-c   --------------
-c
-      REAL cdragh(nlon,nlev)! coeff drag pour T et Q
-      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 
-
-c
-c   Lessivage:
-c   ----------
-c
-c pour le ON-LINE
-c
-      REAL frac_impa(nlon,nlev)  ! fraction d'aerosols impactes
-      REAL frac_nucl(nlon,nlev)  ! fraction d'aerosols nuclees
-c
-cAA
-cAA Arguments necessaires pour les sources et puits de traceur:
-cAA ----------------
-cAA
-      real ftsol(nlon,nbsrf)  ! Temperature du sol (surf)(Kelvin)
-      real pctsrf(nlon,nbsrf) ! Pourcentage de sol f(nature du sol)
-c abder
-      real pftsol1(nlon),pftsol2(nlon),pftsol3(nlon),pftsol4(nlon)
-      real ppsrf1(nlon),ppsrf2(nlon),ppsrf3(nlon),ppsrf4(nlon)
-c fin
-cAA ----------------------------
-cAA  VARIABLES LOCALES TRACEURS
-cAA ----------------------------
-cAA
-cAA Sources et puits des traceurs:
-cAA ------------------------------
-cAA
-cAA Pour l'instant seuls les cas du rn et du pb ont ete envisages.
-
-      REAL source(klon,nbtr)       ! a voir lorsque le flux est prescrit 
-cAA 
-cAA Pour la source de radon et son reservoir de sol
-cAA ................................................
- 
-      REAL,save,allocatable :: trs(:,:)    ! Conc. radon ds le sol
-c$OMP THREADPRIVATE(trs)
-      REAL :: trs_tmp(klon_glo)
-      
-      REAL,save,allocatable :: masktr(:,:) ! Masque reservoir de sol traceur
-c                            Masque de l'echange avec la surface
-c                           (1 = reservoir) ou (possible => 1 )
-c$OMP THREADPRIVATE(masktr)
-      REAL,save,allocatable :: fshtr(:,:)  ! Flux surfacique dans le reservoir de sol
-c$OMP THREADPRIVATE(fshtr)
-      REAL,save,allocatable :: hsoltr(:)      ! Epaisseur equivalente du reservoir de sol
-c$OMP THREADPRIVATE(hsoltr)
-      REAL,save,allocatable :: tautr(:)       ! Constante de decroissance radioactive
-c$OMP THREADPRIVATE(tautr)
-      REAL,save,allocatable :: vdeptr(:)      ! Vitesse de depot sec dans la couche Brownienne
-c$OMP THREADPRIVATE(vdeptr)
-      REAL,save,allocatable :: scavtr(:)      ! Coefficient de lessivage
-c$OMP THREADPRIVATE(scavtr)
-cAA
-      CHARACTER*2 itn
-C maf ioipsl
-      CHARACTER*2 str2
-      INTEGER nhori, nvert
-      REAL zsto, zout, zjulian
-      INTEGER, SAVE :: nid_tra
-c$OMP THREADPRIVATE(nid_tra)
-      INTEGER, SAVE :: nid_tra2,nid_tra3
-c$OMP THREADPRIVATE(nid_tra2,nid_tra3)
-      INTEGER ndex(1)
-      INTEGER ndex2d(iim*(jjm+1)),ndex3d(iim*(jjm+1)*klev)
-      REAL zx_tmp_fi2d(klon)      ! variable temporaire grille physique
-      REAL zx_tmp_fi3d(klon,klev) ! variable temporaire pour champs 3D 
-      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
-      logical ok_sync
-      parameter (ok_sync = .true.)
-C
-C nature du traceur
-c
-      logical,save,allocatable :: aerosol(:)  ! Nature du traceur
-c                            ! aerosol(it) = true  => aerosol 
-c                            ! aerosol(it) = false => gaz 
-c                            ! nat_trac(it) = 1. aerosol
-      logical,save,allocatable :: clsol(:)    ! clsol(it) = true => CL sol calculee
-      logical,save,allocatable :: radio(:)    ! radio(it)=true => decroisssance radioactive
-c$OMP THREADPRIVATE(aerosol,clsol,radio)  
-C
-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,nbtr) ! tendance de traceurs  couche limite
-      REAL d_tr_cv(klon,klev,nbtr) ! tendance de traceurs  conv pour chq traceur
-      REAL d_tr_th(klon,klev,nbtr) ! la tendance des thermiques
-      REAL d_tr_dec(klon,klev,nbtr) ! la tendance de la decroissance 
-c                                   ! radioactive du rn - > pb 
-      REAL d_tr_lessi_impa(klon,klev,nbtr) ! la tendance du lessivage 
-c                                          ! par impaction
-      REAL d_tr_lessi_nucl(klon,klev,nbtr) ! la tendance du lessivage 
-c                                          ! par nucleation 
-      REAL fluxrn(klon,klev) 
-      REAL fluxpb(klon,klev) 
-      REAL pbimpa(klon,klev) 
-      REAL pbnucl(klon,klev) 
-      REAL rn(klon,klev) 
-      REAL pb(klon,klev) 
-      REAL flestottr(klon,klev,nbtr) ! flux de lessivage 
-c                                    ! dans chaque couche 
-      real zmasse(klon,klev)
-      real ztra_th(klon,klev)
-
-C
-      character*20 modname
-      character*80 abort_message
-c
-c   Controles
-c-------------
-      logical first,couchelimite,convection,lessivage,sorties,
-     s        rnpb,inirnpb
-      save first,couchelimite,convection,lessivage,
-     s        sorties,inirnpb
-c$OMP THREADPRIVATE(first,couchelimite,convection,lessivage,
-c$OMP+              sorties,inirnpb)
-c      data first,couchelimite,convection,lessivage,sorties
-c     s     /.true.,.true.,.false.,.true.,.true./
-c Olivia
-       data first,couchelimite,convection,lessivage,
-     s      sorties
-     s     /.true.,.true.,.true.,.true.,.true./
-
-! Variables needed for configuration with INCA
-      INTEGER           :: lastgas
-      INTEGER           :: ncsec
-
-      REAL, PARAMETER   :: dry_mass = 28.966
-      REAL, POINTER     :: hbuf(:)
-      REAL, ALLOCATABLE :: obuf(:)
-      REAL              :: calday
-      REAL              :: pdel(klon,klev)
-c
-c======================================================================
-
-         modname='phytrac'
-
-         ps(:)=paprs(:,1)
-
-         if (debutphy) then
-           allocate( trs(klon,nbtr) )
-	   allocate( masktr(klon,nbtr)) 
-           allocate( fshtr(klon,nbtr) )
-           allocate( hsoltr(nbtr))
-           allocate( tautr(nbtr))
-           allocate( vdeptr(nbtr))
-           allocate( scavtr(nbtr))
-           allocate( aerosol(nbtr))
-           allocate( clsol(nbtr))
-           allocate( radio(nbtr))
-
-
-! FH 2008/05/09 correction de la frequence d'ecriture des traceurs
-!         ecrit_tra = FLOAT(NINT(86400./pdtphys *ecritphy)) 
-          print*,'dans phytrac ',pdtphys,ecrit_tra
-
-         inirnpb=rnpb
-         PRINT*, 'La frequence de sortie traceurs est  ', ecrit_tra
-C         
-c=============================================================
-c   Initialisation des sorties
-c=============================================================
-
-#ifdef CPP_IOIPSL
-#include "ini_histrac.h"
-#endif
-
-c======================================================================
-c   Initialisation de certaines variables pour le Rn et le Pb 
-c======================================================================
-
-c Initialisation du traceur dans le sol (couche limite radonique)
-c
-c        print*,'valeur de debut dans phytrac :',debutphy
-         trs(:,:) = 0.
-c$OMP MASTER         
-       if (is_mpi_root) then
-         trs_tmp(:)=0.
-	 open (99,file='starttrac',status='old',
-     .         err=999,form='formatted')
-         read(99,*) (trs_tmp(i),i=1,klon_glo)
-999      close(99)
-       endif
-c$OMP END MASTER
-       call Scatter(trs_tmp,trs(:,1))
-
-c         print*, 'apres starttrac'
-
-c Initialisation de la fraction d'aerosols lessivee
-c
-         d_tr_lessi_impa(:,:,:) = 0.
-         d_tr_lessi_nucl(:,:,:) = 0. 
-c
-c Initialisation de la nature des traceurs
-c
-         DO it = 1, nbtr
-            aerosol(it) = .FALSE.  ! Tous les traceurs sont des gaz par defaut
-            radio(it) = .FALSE.    ! Par defaut pas de passage par radiornpb
-            clsol(it) = .FALSE.  ! Par defaut couche limite avec flux prescrit
-         ENDDO
-c
-      ENDIF  ! fin debutphy 
-c Initialisation du traceur dans le sol (couche limite radonique)
-      if(inirnpb) THEN
-c
-         radio(1)= .true.
-         radio(2)= .true.
-         clsol(1)= .true.
-         clsol(2)= .true.
-         aerosol(2) = .TRUE. ! le Pb est un aerosol 
-c
-         call initrrnpb (ftsol,pctsrf,masktr,fshtr,hsoltr,tautr
-     .                   ,vdeptr,scavtr)
-         inirnpb=.false.
-      endif
-
-
-      IF (config_inca == 'none') THEN
-       DO i=1,nlon
-          pftsol1(i) = ftsol(i,1)
-          pftsol2(i) = ftsol(i,2)
-          pftsol3(i) = ftsol(i,3)
-          pftsol4(i) = ftsol(i,4)
-
-          ppsrf1(i) = pctsrf(i,1)
-          ppsrf2(i) = pctsrf(i,2)
-          ppsrf3(i) = pctsrf(i,3)
-          ppsrf4(i) = pctsrf(i,4)
-
-       ENDDO
-
-      ELSE ! config_inca /=none
-#ifdef INCA
-      call VTe(VTphysiq)
-      call VTb(VTinca)
-
-!======================================================================
-!     Chimie
-!======================================================================
-
-        calday = FLOAT(julien) + gmtime
-        ncsec  = NINT (86400.*gmtime)
-
-        DO k = 1, nlev
-        pdel(:,k) = paprs(:,k) - paprs (:,k+1)
-        END DO
-
-        IF (config_inca == 'aero') THEN
-           CALL aerosolmain (aerosol_couple,
-     $                 tr_seri,
-     $                 pdtphys,
-     $                 pplay,
-     $                 pdel,
-     $                 prfl,
-     $                 pmflxr,
-     $                 psfl,
-     $                 pmflxs,
-     $                 pmfu,
-     $                 itop_con,
-     $                 ibas_con,
-     $                 pphi,
-     $                 airephy, ! paire,
-     $                 nstep,
-     $                 rneb,         ! for chimiaq
-     $                 t_seri,       ! for chimiaq
-     $                 rh,
-     $                 tau_aero,
-     $                 piz_aero,
-     $                 cg_aero,
-     $                 rfname,
-     $                 ccm,
-     $                 lafin)
-        END IF
-
-        CALL chemmain (tr_seri,    !mmr
-     $                 nstep,      !nstep
-     $                 calday,     !calday
-     $                 julien,     !ncdate
-     $                 ncsec,      !ncsec
-     $                 1,          !lat
-     $                 pdtphys,    !delt
-     $                 paprs(1,1), !ps
-     $                 pplay,      !pmid
-     $                 pdel,       !pdel
-     $                 airephy,
-     $                 pctsrf(1,1),!oro
-     $                 ftsol,      !tsurf
-     $                 albsol,     !albs
-     $                 pphi,       !zma
-     $                 pphis,      !phis
-     $                 cldfra,     !cldfr
-     $                 rneb,       !cldfr_st
-     $                 diafra,     !cldfr_cv
-     $                 itop_con,   !cldtop
-     $                 ibas_con,   !cldbot
-     $                 cldliq,     !cwat
-     $                 prfl,       !flxrst
-     $                 pmflxr,     !flxrcv
-     $                 psfl,       !flxsst
-     $                 pmflxs,     !flxscv
-     $                 pmfu,       !flxupd
-     $                 flxmass_w,  !flxmass_w
-     $                 t_seri,     !tfld
-     $                 sh,         !sh
-     $                 rh,         !rh
-     $                 iip1,       !nx
-     $                 jjp1,       !ny
-     $                 source,
-     $                 solsym)
-
-
-      call VTe(VTinca)
-      call VTb(VTphysiq)
-#endif
-      END IF ! config_inca
-c======================================================================
-c   Calcul de l'effet de la convection
-c======================================================================
-c     print*,'Avant convection'
-       do it=1,nbtr
-          WRITE(itn,'(i2)') it
-c        call diagtracphy(tr_seri(:,:,it),paprs,'Avant conv'//itn)
-       enddo
-
-      if (convection) then
-
-c      print*,'Pas de temps dans phytrac : ',pdtphys
-      DO it=1, nbtr
-
-      IF ( config_inca/='none') THEN 
-         IF ( conv_flg(it) == 0 ) CYCLE
-      END IF
-
-      if (iflag_con.lt.2) then
-       d_tr_cv=0.
-      else if (iflag_con.eq.2) then
-c tiedke
-      CALL nflxtr(pdtphys, pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
-     .            pplay, paprs, tr_seri(1,1,it), d_tr_cv(1,1,it))
-      else
-c KE
-      call cvltr(pdtphys, da, phi, mp, paprs,pplay, tr_seri(1,1,it),
-     .           upwd,dnwd,d_tr_cv(1,1,it))
-      endif
-
-       DO k = 1, nlev
-       DO i = 1, klon
-         tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_cv(i,k,it)
-       ENDDO
-       ENDDO
-
-       IF (config_inca == 'none') THEN
-        CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'convection it = '//itn)
-       ELSE
-        CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'convection it = '
-     .         //solsym(it))
-       END IF   
-      
-      ENDDO
-
-      endif ! convection
-c        print*,'Apres convection'
-c      do it=1,nbtr
-c         WRITE(itn,'(i1)') it
-c        call diagtracphy(tr_seri(:,:,it),paprs,'Avant conv'//itn)
-c      enddo
-
-
-c======================================================================
-c   Calcul de l'effet des thermiques
-c======================================================================
-
-      do k=1,klev
-         do i=1,klon
-            zmasse(i,k)=(paprs(i,k)-paprs(i,k+1))/rg
-         enddo
-      enddo
-
-c      print*,'masse dans ph ',zmasse
-      do it=1,nbtr
-         do k=1,klev
-            do i=1,klon
-               d_tr_th(i,k,it)=0.
-               tr_seri(i,k,it)=max(tr_seri(i,k,it),0.)
-               tr_seri(i,k,it)=min(tr_seri(i,k,it),1.e10)
-            enddo
-         enddo
-      enddo
-
-      if (iflag_thermals.gt.0) then
-c        print*,'calcul de leffet des thermiques'
-        nsplit=10
-        DO it=1, nbtr
-c        WRITE(itn,'(i1)') it
-c        CALL minmaxqfi(tr_seri(1,1,it),1.e10,-1.e33,'conv it='//itn)
-c	     print*,'avant dqthermiquesretro'
-c             call dump2d(iim,jjm-1,tr_seri(2,1,1),'TR_SERI      ')
-
-         do isplit=1,nsplit
-c  Abderr 25 11 02
-C Thermiques 
-c	print*,'Avant dans phytrac'
-            call dqthermcell(klon,klev,pdtphys/nsplit
-     .       ,fm_therm,entr_therm,zmasse
-     .       ,tr_seri(1:klon,1:klev,it),d_tr,ztra_th)
-
-            do k=1,klev
-               do i=1,klon
-                  d_tr(i,k)=pdtphys*d_tr(i,k)/nsplit
-                  d_tr_th(i,k,it)=d_tr_th(i,k,it)+d_tr(i,k)
-                  tr_seri(i,k,it)=max(tr_seri(i,k,it)+d_tr(i,k),0.)
-               enddo
-            enddo
-          enddo ! nsplit
-c          print*,'apres thermiques'
-c             call dump2d(iim,jjm-1,d_tr_th(1,1,1),'d_tr_th      ')
-c            do k=1,klev
-c	print*,'d_tr_th(',k,')=',tr_seri(280,k,1)
-c	   enddo
-
-c        WRITE(itn,'(i1)') it
-c        CALL minmaxqfi(tr_seri(1,1,it),1.e10,-1.e33,'therm it='//itn)
-       ENDDO ! it
-       endif ! Thermiques
-c       print*,'ATTENTION: sdans thermniques'
-      
-c======================================================================
-c   Calcul de l'effet de la couche limite
-c======================================================================
-c	print *,'Avant couchelimite'
-c      do it=1,nbtr
-c         WRITE(itn,'(i1)') it
-c        call diagtracphy(tr_seri(:,:,it),paprs,'Avant CL  '//itn)
-c      enddo
-
-      if (couchelimite) then
-
-      DO k = 1, nlev
-      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, nbtr
-
-         IF ( config_inca/='none' ) THEN
-            IF( pbl_flg(it) == 0 ) CYCLE
-         END IF
-
-c     print *,'it',it,clsol(it)
-      if (clsol(it)) then  ! couche limite avec quantite dans le sol calculee
-          CALL cltracrn(it, pdtphys, yu1, yv1,
-     e                    cdragh, coefh,t_seri,ftsol,pctsrf,
-     e                    tr_seri(1,1,it),trs(1,it),
-     e                    paprs, pplay, delp,
-     e                    masktr(1,it),fshtr(1,it),hsoltr(it),
-     e                    tautr(it),vdeptr(it),
-     e                    xlat,
-     s                    d_tr_cl(1,1,it),d_trs)
-          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
-c
-c Traceur ds sol
-c
-          DO i = 1, klon
-            trs(i,it) = trs(i,it) + d_trs(i)
-          END DO
-C
-C maf provisoire suppression des prints
-C         WRITE(itn,'(i1)') it
-C         CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'cltracrn it='//itn)
-      else ! couche limite avec flux prescrit
-
-         IF (config_inca == 'none') THEN
-Cmaf provisoire source / traceur a creer
-            DO i=1, klon
-               source(i,it) = 0.0 ! pas de source, pour l'instant
-            ENDDO
-         END IF
-
-          CALL cltrac(pdtphys, coefh,t_seri,
-     s               tr_seri(1,1,it), source(:,it),
-     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      print*,'Apres couchelimite'
-c      do it=1,nbtr
-c         WRITE(itn,'(i1)') it
-c        call diagtracphy(tr_seri(:,:,it),paprs,'Avant CL  '//itn)
-c      enddo
-
-c======================================================================
-c   Calcul de l'effet du puits radioactif
-c======================================================================
-
-C MAF il faudrait faire une modification pour passer dans radiornpb 
-c si radio=true mais pour l'instant radiornpb propre au cas rnpb
-      if(rnpb) then
-c       print *, 'decroissance radiactive activee'
-        call radiornpb (tr_seri,pdtphys,tautr,d_tr_dec)
-C
-        DO it=1,nbtr
-            if(radio(it)) then
-            DO k = 1, nlev
-               DO i = 1, klon
-                  tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_dec(i,k,it)
-               ENDDO
-            ENDDO
-            WRITE(itn,'(i1)') it
-            CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'puits rn it='//itn)
-            endif
-        ENDDO
-c
-      endif ! rnpb decroissance  radioactive
-C
-c======================================================================
-c   Calcul de l'effet de la precipitation
-c======================================================================
-
-c      print*,'LESSIVAGE =',lessivage
-      IF (lessivage) THEN
-
-c     print*,'avant lessivage'
-
-          d_tr_lessi_nucl(:,:,:) = 0. 
-          d_tr_lessi_impa(:,:,:) = 0. 
-          flestottr(:,:,:) = 0. 
-c
-c tendance des aerosols nuclees et impactes 
-c
-       DO it = 1, nbtr
-         IF (aerosol(it)) THEN
-           DO k = 1, nlev
-              DO i = 1, klon
-               d_tr_lessi_nucl(i,k,it) = d_tr_lessi_nucl(i,k,it) +
-     s                  ( 1 - frac_nucl(i,k) )*tr_seri(i,k,it)
-               d_tr_lessi_impa(i,k,it) = d_tr_lessi_impa(i,k,it) +
-     s                  ( 1 - frac_impa(i,k) )*tr_seri(i,k,it)
-              ENDDO
-           ENDDO
-         ENDIF
-       ENDDO
-c
-c Mises a jour des traceurs + calcul des flux de lessivage 
-c Mise a jour due a l'impaction et a la nucleation
-c
-c      call dump2d(iim,jjm-1,frac_impa(2:klon-1,10),'FRACIMPA')
-c      call dump2d(iim,jjm-1,frac_nucl(2:klon-1,10),'FRACNUCL')
-c      call dump2d(iim,jjm-1,tr_seri(2:klon-1,10,3),'TRACEUR3')
-       DO it = 1, nbtr
-c         print*,'IT=',it,aerosol(it)
-         IF (aerosol(it)) THEN
-c           print*,'IT=',it,' On lessive'
-           DO k = 1, nlev
-              DO i = 1, klon
-               tr_seri(i,k,it)=tr_seri(i,k,it)
-     s         *frac_impa(i,k)*frac_nucl(i,k)
-              ENDDO
-           ENDDO
-         ENDIF
-       ENDDO
-c      call dump2d(iim,jjm-1,tr_seri(2:klon-1,10,3),'TRACEUR3B')
-c
-c Flux lessivage total 
-c
-      DO it = 1, nbtr
-           DO k = 1, nlev
-            DO i = 1, klon
-               flestottr(i,k,it) = flestottr(i,k,it) -
-     s                   ( d_tr_lessi_nucl(i,k,it)   +
-     s                     d_tr_lessi_impa(i,k,it) ) *
-     s                   ( paprs(i,k)-paprs(i,k+1) ) / 
-     s                   (RG * pdtphys)
-            ENDDO
-           ENDDO
-c
-Cmaf        WRITE(itn,'(i1)') it
-Cmaf    CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'tr(lessi) it='//itn)
-      ENDDO
-c
-c     print*,'apres lessivage'
-      ENDIF
-Cc
-      DO k = 1, nlev
-         DO i = 1, klon
-            fluxrn(i,k) = flestottr(i,k,1)
-            fluxpb(i,k) = flestottr(i,k,2)
-            rn(i,k) = tr_seri(i,k,1)
-            pb(i,k) = tr_seri(i,k,2)
-            pbnucl(i,k)=d_tr_lessi_nucl(i,k,2)
-            pbimpa(i,k)=d_tr_lessi_impa(i,k,2)
-         ENDDO
-      ENDDO
-
-c=============================================================
-c   Ecriture des sorties
-c=============================================================
-
-#ifdef CPP_IOIPSL
-#include "write_histrac.h"
-#endif
-
-c=============================================================
-
-      if (lafin) then
-         print*, 'c est la fin de la physique'
-	 call Gather(trs(:,1),trs_tmp)
-c$OMP MASTER	 
-	 if (is_mpi_root) then
-         
-	   open (99,file='restarttrac',  form='formatted')
-           do i=1,klon_glo
-               write(99,*) trs_tmp(i)
-           enddo
-           PRINT*, 'Ecriture du fichier restarttrac'
-           close(99)
-	 endif
-c$OMP END MASTER
-      else
-c         print*, 'physique pas fini'
-      endif
-
-
-      RETURN
-      END
Index: /LMDZ4/branches/LMDZ4-dev/libf/phylmd/phytrac.F90
===================================================================
--- /LMDZ4/branches/LMDZ4-dev/libf/phylmd/phytrac.F90	(revision 1191)
+++ /LMDZ4/branches/LMDZ4-dev/libf/phylmd/phytrac.F90	(revision 1191)
@@ -0,0 +1,410 @@
+!$Id $
+
+SUBROUTINE phytrac(                            &
+     nstep,     julien,   gmtime,   debutphy,  &
+     lafin,     pdtphys,  u, v,     t_seri,     &
+     paprs,     pplay,    pmfu,     pmfd,      &
+     pen_u,     pde_u,    pen_d,    pde_d,     &
+     cdragh,    coefh,    fm_therm, entr_therm,&
+     yu1,       yv1,      ftsol,    pctsrf,    &
+     xlat,      frac_impa,frac_nucl,xlon,      &
+     presnivs,  pphis,    pphi,     albsol,    &
+     sh,        rh,       cldfra,   rneb,      &
+     diafra,    cldliq,   itop_con, ibas_con,  &
+     pmflxr,    pmflxs,   prfl,     psfl,      &
+     da,        phi,      mp,       upwd,      &
+     dnwd,      aerosol_couple,     flxmass_w, &
+     tau_aero,  piz_aero,  cg_aero, ccm,       &
+     rfname,                                   &
+     tr_seri)         
+!
+!======================================================================
+! Auteur(s) FH
+! Objet: Moniteur general des tendances traceurs
+!======================================================================
+
+  USE ioipsl
+  USE dimphy
+  USE infotrac
+  USE mod_grid_phy_lmdz
+  USE mod_phys_lmdz_para
+  USE comgeomphy
+  USE iophy
+  USE traclmdz_mod
+  USE tracinca_mod
+
+
+  IMPLICIT NONE
+
+  INCLUDE "YOMCST.h"
+  INCLUDE "dimensions.h"
+  INCLUDE "indicesol.h"
+  INCLUDE "clesphys.h"
+  INCLUDE "temps.h"
+  INCLUDE "paramet.h"
+  INCLUDE "control.h"
+  INCLUDE "thermcell.h"
+!==========================================================================
+!                   -- ARGUMENT DESCRIPTION --
+!==========================================================================
+
+! Input arguments
+!----------------
+!Configuration grille,temps:
+  INTEGER,INTENT(IN) :: nstep      ! Appel physique
+  INTEGER,INTENT(IN) :: julien     ! Jour julien
+  REAL,INTENT(IN)    :: gmtime
+  REAL,INTENT(IN)    :: pdtphys    ! Pas d'integration pour la physique (seconde)
+  LOGICAL,INTENT(IN) :: debutphy   ! le flag de l'initialisation de la physique
+  LOGICAL,INTENT(IN) :: lafin      ! le flag de la fin de la physique
+  
+  REAL,DIMENSION(klon),INTENT(IN) :: xlat    ! latitudes pour chaque point 
+  REAL,DIMENSION(klon),INTENT(IN) :: xlon    ! longitudes pour chaque point 
+!
+!Physique: 
+!--------
+  REAL,DIMENSION(klon,klev),INTENT(IN)   :: t_seri  ! Temperature
+  REAL,DIMENSION(klon,klev),INTENT(IN)   :: u       ! 
+  REAL,DIMENSION(klon,klev),INTENT(IN)   :: v       ! 
+  REAL,DIMENSION(klon,klev),INTENT(IN)   :: sh      ! humidite specifique
+  REAL,DIMENSION(klon,klev),INTENT(IN)   :: rh      ! humidite relative
+  REAL,DIMENSION(klon,klev+1),INTENT(IN) :: paprs   ! pression pour chaque inter-couche (en Pa)
+  REAL,DIMENSION(klon,klev),INTENT(IN)   :: pplay   ! pression pour le mileu de chaque couche (en Pa)
+  REAL,DIMENSION(klon,klev),INTENT(IN)   :: pphi    ! geopotentiel
+  REAL,DIMENSION(klon),INTENT(IN)        :: pphis
+  REAL,DIMENSION(klev),INTENT(IN)        :: presnivs 
+  REAL,DIMENSION(klon,klev),INTENT(IN)   :: cldliq  ! eau liquide nuageuse
+  REAL,DIMENSION(klon,klev),INTENT(IN)   :: cldfra  ! fraction nuageuse (tous les nuages)
+  REAL,DIMENSION(klon,klev),INTENT(IN)   :: diafra  ! fraction nuageuse (convection ou stratus artificiels)
+  REAL,DIMENSION(klon,klev),INTENT(IN)   :: rneb    ! fraction nuageuse (grande echelle)
+  INTEGER,DIMENSION(klon),INTENT(IN)     :: itop_con
+  INTEGER,DIMENSION(klon),INTENT(IN)     :: ibas_con
+  REAL,DIMENSION(klon),INTENT(IN)        :: albsol  ! albedo surface
+!
+!Convection:
+!----------
+  REAL,DIMENSION(klon,klev),INTENT(IN) :: pmfu  ! flux de masse dans le panache montant
+  REAL,DIMENSION(klon,klev),INTENT(IN) :: pmfd  ! flux de masse dans le panache descendant
+  REAL,DIMENSION(klon,klev),INTENT(IN) :: pen_u ! flux entraine dans le panache montant
+  REAL,DIMENSION(klon,klev),INTENT(IN) :: pde_u ! flux detraine dans le panache montant
+  REAL,DIMENSION(klon,klev),INTENT(IN) :: pen_d ! flux entraine dans le panache descendant
+  REAL,DIMENSION(klon,klev),INTENT(IN) :: pde_d ! flux detraine dans le panache descendant
+
+!...Tiedke     
+  REAL,DIMENSION(klon,klev+1),INTENT(IN)   :: pmflxr, pmflxs ! Flux precipitant de pluie, neige aux interfaces [convection]
+  REAL,DIMENSION(klon,klev+1),INTENT(IN)   :: prfl, psfl ! Flux precipitant de pluie, neige aux interfaces [large-scale]
+
+  LOGICAL,INTENT(IN)                       :: aerosol_couple
+  REAL,DIMENSION(klon,klev),INTENT(IN)     :: flxmass_w
+  REAL,DIMENSION(klon,klev,9,2),INTENT(IN) :: tau_aero
+  REAL,DIMENSION(klon,klev,9,2),INTENT(IN) :: piz_aero
+  REAL,DIMENSION(klon,klev,9,2),INTENT(IN) :: cg_aero
+  CHARACTER(len=4),DIMENSION(9),INTENT(IN) :: rfname 
+  REAL,DIMENSION(klon,klev,2),INTENT(IN)   :: ccm 
+!... K.Emanuel
+  REAL,DIMENSION(klon,klev),INTENT(IN)     :: da
+  REAL,DIMENSION(klon,klev,klev),INTENT(IN):: phi
+  REAL,DIMENSION(klon,klev),INTENT(IN)     :: mp
+  REAL,DIMENSION(klon,klev),INTENT(IN)     :: upwd      ! saturated updraft mass flux
+  REAL,DIMENSION(klon,klev),INTENT(IN)     :: dnwd      ! saturated downdraft mass flux
+!
+!Thermiques:
+!----------
+  REAL,DIMENSION(klon,klev+1),INTENT(IN)   :: fm_therm
+  REAL,DIMENSION(klon,klev),INTENT(IN)     :: entr_therm
+!
+!Couche limite:
+!--------------
+!
+  REAL,DIMENSION(klon,klev),INTENT(IN) :: cdragh ! coeff drag pour T et Q
+  REAL,DIMENSION(klon,klev),INTENT(IN) :: coefh  ! coeff melange CL (m**2/s)
+  REAL,DIMENSION(klon),INTENT(IN)      :: yu1    ! vents au premier niveau
+  REAL,DIMENSION(klon),INTENT(IN)      :: yv1    ! vents au premier niveau
+!
+!Lessivage:
+!----------
+!
+! pour le ON-LINE
+!
+  REAL,DIMENSION(klon,klev),INTENT(IN) :: frac_impa ! fraction d'aerosols non impactes
+  REAL,DIMENSION(klon,klev),INTENT(IN) :: frac_nucl ! fraction d'aerosols non nuclees
+
+! Arguments necessaires pour les sources et puits de traceur:
+  REAL,DIMENSION(klon,nbsrf),INTENT(IN) :: ftsol  ! Temperature du sol (surf)(Kelvin)
+  REAL,DIMENSION(klon,nbsrf),INTENT(IN) :: pctsrf ! Pourcentage de sol (nature du sol)
+
+
+! Output argument
+!----------------
+  REAL,DIMENSION(klon,klev,nbtr),INTENT(INOUT) :: tr_seri ! Concentration Traceur [U/KgA]  
+
+!=======================================================================================
+!                        -- LOCAL VARIABLES --
+!=======================================================================================
+
+  INTEGER :: i, k, it
+  INTEGER :: nsplit
+
+!Sources et Reservoirs de traceurs (ex:Radon):
+!--------------------------------------------
+!
+  REAL,DIMENSION(:,:),ALLOCATABLE,SAVE :: source  ! a voir lorsque le flux de surface est prescrit 
+!$OMP THREADPRIVATE(source)
+
+!
+!Entrees/Sorties: (cf ini_histrac.h et write_histrac.h)  
+!---------------
+  INTEGER                   :: iiq, ierr
+  INTEGER                   :: nhori, nvert
+  REAL                      :: zsto, zout, zjulian
+  INTEGER,SAVE              :: nid_tra     ! pointe vers le fichier histrac.nc         
+!$OMP THREADPRIVATE(nid_tra)
+  REAL,DIMENSION(klon)      :: zx_tmp_fi2d ! variable temporaire grille physique
+  INTEGER                   :: itau_w      ! pas de temps ecriture = nstep + itau_phy
+  LOGICAL,PARAMETER :: ok_sync=.TRUE.
+
+!
+! Nature du traceur
+!------------------
+  LOGICAL,DIMENSION(:),ALLOCATABLE,SAVE :: aerosol  ! aerosol(it) = true  => aerosol => lessivage
+!$OMP THREADPRIVATE(aerosol)                        ! aerosol(it) = false => gaz
+  REAL,DIMENSION(klon,klev)             :: delp     ! epaisseur de couche (Pa)
+!
+! Tendances de traceurs (Td):
+!------------------------
+!
+  REAL,DIMENSION(klon,klev)      :: d_tr     ! Td dans l'atmosphere
+  REAL,DIMENSION(klon,klev,nbtr) :: d_tr_cl  ! Td couche limite/traceur
+  REAL,DIMENSION(klon,klev,nbtr) :: d_tr_cv  ! Td convection/traceur
+  REAL,DIMENSION(klon,klev,nbtr) :: d_tr_th  ! Td thermique
+  REAL,DIMENSION(klon,klev,nbtr) :: d_tr_lessi_impa ! Td du lessivage par impaction
+  REAL,DIMENSION(klon,klev,nbtr) :: d_tr_lessi_nucl ! Td du lessivage par nucleation 
+!
+! Physique
+!----------   
+  REAL,DIMENSION(klon,klev,nbtr) :: flestottr ! flux de lessivage dans chaque couche 
+  REAL,DIMENSION(klon,klev)      :: zmasse    ! densité atmosphérique Kg/m2
+  REAL,DIMENSION(klon,klev)      :: ztra_th
+  
+!Controles:
+!---------
+  LOGICAL,SAVE :: couchelimite=.TRUE.
+  LOGICAL,SAVE :: convection=.TRUE.
+  LOGICAL,SAVE :: lessivage
+!$OMP THREADPRIVATE(couchelimite,convection,lessivage)
+
+  CHARACTER(len=8),DIMENSION(nbtr) :: solsym
+
+
+!######################################################################
+!                    -- INITIALIZATION --
+!######################################################################
+  IF (debutphy) THEN
+     WRITE(*,*) 'FIRST TIME IN PHYTRAC : pdtphys(sec) = ',pdtphys,'ecrit_tra (sec) = ',ecrit_tra
+     ALLOCATE( source(klon,nbtr), stat=ierr)
+     IF (ierr /= 0) CALL abort_gcm('phytrac', 'pb in allocation 1',1)
+     
+     ALLOCATE( aerosol(nbtr), stat=ierr)
+     IF (ierr /= 0) CALL abort_gcm('phytrac', 'pb in allocation 2',1)
+     
+
+     ! Initialize module for specific tracers
+     SELECT CASE(type_trac)
+     CASE('lmdz')
+        CALL traclmdz_init(pctsrf, ftsol, aerosol, lessivage)
+     CASE('inca')
+        source(:,:)=0.
+        CALL tracinca_init(aerosol,lessivage)
+     END SELECT
+!
+! Initialize diagnostic output
+! ----------------------------
+     INCLUDE "ini_histrac.h"
+        
+  END IF
+!############################################ END INITIALIZATION #######
+
+!===============================================================================
+!    -- Do specific treatment according to chemestry model or local LMDZ tracers
+!      
+!===============================================================================
+  SELECT CASE(type_trac)
+  CASE('lmdz')
+     !    -- Traitement des traceurs avec traclmdz
+     
+     CALL traclmdz(&
+          pdtphys,  t_seri,                         &
+          paprs,    pplay,        cdragh,  coefh,   &
+          yu1,      yv1,          ftsol,   pctsrf,  &
+          xlat,     couchelimite,                   &
+          tr_seri,  source,       solsym,  d_tr_cl)
+     
+  CASE('inca')
+     !    -- CHIMIE INCA  config_inca = aero or chem --
+
+     CALL tracinca(&
+          nstep,    julien,   gmtime,         lafin,     &
+          pdtphys,  t_seri,   paprs,          pplay,     &
+          pmfu,     ftsol,    pctsrf,         pphis,     &
+          pphi,     albsol,   sh,             rh,        &
+          cldfra,   rneb,     diafra,         cldliq,    &
+          itop_con, ibas_con, pmflxr,         pmflxs,    &
+          prfl,     psfl,     aerosol_couple, flxmass_w, &
+          tau_aero, piz_aero, cg_aero,        ccm,       &
+          rfname,                                        &
+          tr_seri,  source,   solsym)      
+  END SELECT
+
+!======================================================================
+!       -- Calcul de l'effet de la convection --
+!======================================================================
+  IF (convection) THEN
+     DO it=1, nbtr
+        IF ( conv_flg(it) == 0 ) CYCLE
+        
+        IF (iflag_con.LT.2) THEN
+           d_tr_cv(:,:,:)=0.
+        ELSE IF (iflag_con.EQ.2) THEN
+!..Tiedke
+           CALL nflxtr(pdtphys, pmfu, pmfd, pen_u, pde_u, pen_d, pde_d, &
+                pplay, paprs, tr_seri(:,:,it), d_tr_cv(:,:,it))
+        ELSE
+!..K.Emanuel
+           CALL cvltr(pdtphys, da, phi, mp, paprs,pplay, tr_seri(:,:,it),&
+                upwd,dnwd,d_tr_cv(:,:,it))
+        END IF
+
+        DO k = 1, klev
+           DO i = 1, klon        
+              tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_cv(i,k,it)
+           END DO
+        END DO
+
+        CALL minmaxqfi(tr_seri(:,:,it),0.,1.e33,'convection it = '//solsym(it))
+             
+     END DO ! nbtr
+  END IF ! convection
+
+!======================================================================
+!    -- Calcul de l'effet des thermiques --
+!======================================================================
+
+  DO k=1,klev
+     DO i=1,klon
+        zmasse(i,k)=(paprs(i,k)-paprs(i,k+1))/rg
+     END DO
+  END DO
+
+  DO it=1,nbtr
+     DO k=1,klev
+        DO i=1,klon
+           d_tr_th(i,k,it)=0.
+           tr_seri(i,k,it)=MAX(tr_seri(i,k,it),0.)
+           tr_seri(i,k,it)=MIN(tr_seri(i,k,it),1.e10)
+        END DO
+     END DO
+  END DO
+  
+  IF (iflag_thermals.GT.0) THEN   
+     nsplit=10
+     DO it=1, nbtr
+        DO isplit=1,nsplit
+
+           CALL dqthermcell(klon,klev,pdtphys/nsplit, &
+                fm_therm,entr_therm,zmasse, &
+                tr_seri(1:klon,1:klev,it),d_tr,ztra_th)
+
+           DO k=1,klev
+              DO i=1,klon
+                 d_tr(i,k)=pdtphys*d_tr(i,k)/nsplit
+                 d_tr_th(i,k,it)=d_tr_th(i,k,it)+d_tr(i,k)
+                 tr_seri(i,k,it)=MAX(tr_seri(i,k,it)+d_tr(i,k),0.)
+              END DO
+           END DO
+        END DO ! nsplit
+     END DO ! it
+  END IF ! Thermiques
+
+!======================================================================
+!     -- Calcul de l'effet de la couche limite --
+!======================================================================
+
+  IF (couchelimite) THEN
+
+     DO k = 1, klev
+        DO i = 1, klon
+           delp(i,k) = paprs(i,k)-paprs(i,k+1)
+        END DO
+     END DO
+
+     DO it=1, nbtr
+        
+        IF( pbl_flg(it) /= 0 ) THEN
+        
+           CALL cltrac(pdtphys, coefh,t_seri,       &
+                tr_seri(:,:,it), source(:,it),      &
+                paprs, pplay, delp,                 &
+                d_tr_cl(:,:,it))
+           
+           DO k = 1, klev
+              DO i = 1, klon
+                 tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_cl(i,k,it)
+              END DO
+           END DO
+        END IF
+
+     END DO
+     
+  END IF ! couche limite
+
+
+!======================================================================
+!   Calcul de l'effet de la precipitation
+!======================================================================
+
+  IF (lessivage) THEN
+     
+     d_tr_lessi_nucl(:,:,:) = 0. 
+     d_tr_lessi_impa(:,:,:) = 0.
+     flestottr(:,:,:) = 0. 
+!=========================
+! LESSIVAGE LARGE SCALE : 
+!=========================
+
+! Tendance des aerosols nuclees et impactes 
+! -----------------------------------------
+     DO it = 1, nbtr
+        IF (aerosol(it)) THEN
+           DO k = 1, klev
+              DO i = 1, klon
+                 d_tr_lessi_nucl(i,k,it) = d_tr_lessi_nucl(i,k,it) +    &
+                      ( 1 - frac_nucl(i,k) )*tr_seri(i,k,it)
+                 d_tr_lessi_impa(i,k,it) = d_tr_lessi_impa(i,k,it) +    &
+                      ( 1 - frac_impa(i,k) )*tr_seri(i,k,it)
+
+!
+! Flux lessivage total 
+! ------------------------------------------------------------
+                 flestottr(i,k,it) = flestottr(i,k,it) -   &
+                      ( d_tr_lessi_nucl(i,k,it)   +        &
+                      d_tr_lessi_impa(i,k,it) ) *          &
+                      ( paprs(i,k)-paprs(i,k+1) ) /        &
+                      (RG * pdtphys)
+!
+! Mise a jour des traceurs due a l'impaction,nucleation 
+! ----------------------------------------------------------------------
+                 tr_seri(i,k,it)=tr_seri(i,k,it)*frac_impa(i,k)*frac_nucl(i,k)
+              END DO
+           END DO
+        END IF
+     END DO
+     
+  END IF ! lessivage
+
+!=============================================================
+!   Ecriture des sorties
+!=============================================================
+  INCLUDE "write_histrac.h"
+
+
+END SUBROUTINE phytrac
Index: /LMDZ4/branches/LMDZ4-dev/libf/phylmd/radio_decay.F90
===================================================================
--- /LMDZ4/branches/LMDZ4-dev/libf/phylmd/radio_decay.F90	(revision 1191)
+++ /LMDZ4/branches/LMDZ4-dev/libf/phylmd/radio_decay.F90	(revision 1191)
@@ -0,0 +1,61 @@
+!
+! $Id $
+!
+SUBROUTINE radio_decay(radio,rnpb,dtime,tautr,tr,d_tr) 
+!
+! Caluclate radioactive decay for all tracers with radio(it)=true
+!
+  USE dimphy
+  USE infotrac, ONLY : nbtr
+  IMPLICIT NONE
+!-----------------------------------------------------------------------
+! Auteur(s): AA + CG (LGGE/CNRS) Date 24-06-94
+! Objet: Calcul de la tendance radioactive des traceurs type radioelements
+!CG240694 : Pour un traceur, le radon
+!CG161294 : Plus un 2eme traceur, le 210Pb. Le radon decroit en plomb.
+!-----------------------------------------------------------------------
+!
+! Entrees
+!
+  LOGICAL,DIMENSION(nbtr),INTENT(IN)        :: radio ! .true. = traceur radioactif  
+  LOGICAL,INTENT(IN)                        :: rnpb  ! .true. = decroissance RN = source PB
+  REAL,INTENT(IN)                           :: dtime ! Pas de temps physique (secondes)
+  REAL,DIMENSION(nbtr),INTENT(IN)           :: tautr ! Constante de decroissance radioactive
+  REAL,DIMENSION(klon,klev,nbtr),INTENT(IN) :: tr    ! Concentrations traceurs U/kgA
+!
+! Sortie
+!
+  REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: d_tr  ! Tendance de decroissance radioactive
+!
+! Locales
+!
+  INTEGER  :: i,k,it
+
+
+  DO it = 1,nbtr
+     IF ( radio(it) ) THEN
+        IF (tautr(it) .GT. 0.) THEN
+           DO k = 1,klev
+              DO i = 1,klon
+                 d_tr(i,k,it) = - tr(i,k,it) * dtime / tautr(it)
+              END DO
+           END DO
+        ELSE 
+           d_tr(:,:,it) = 0.
+        END IF
+     ELSE
+        d_tr(:,:,it) = 0.
+     END IF
+  END DO
+!-------------------------------------------------------
+!CG161294 : Cas particulier radon [it=1] => plomb [it=2]
+!-------------------------------------------------------
+  IF ( rnpb ) THEN
+     DO k = 1,klev
+        DO i = 1,klon
+           d_tr(i,k,2) = d_tr(i,k,2) - d_tr(i,k,1)
+        ENDDO
+     ENDDO
+  ENDIF
+
+END SUBROUTINE radio_decay
Index: DZ4/branches/LMDZ4-dev/libf/phylmd/radiornpb.F
===================================================================
--- /LMDZ4/branches/LMDZ4-dev/libf/phylmd/radiornpb.F	(revision 1190)
+++ 	(revision )
@@ -1,55 +1,0 @@
-!
-! $Header$
-!
-      SUBROUTINE radiornpb(tr,dtime,tautr,d_tr) 
-      USE dimphy
-      USE infotrac, ONLY : nbtr
-      IMPLICIT none
-c======================================================================
-c Auteur(s): AA + CG (LGGE/CNRS) Date 24-06-94
-c Objet: Decroissance radioactive d'un traceur dans l'atmosphere
-CG240694 : Pour un traceur, le radon
-CG161294 : Plus un 2eme traceur, le 210Pb. Le radon decroit en plomb.
-c======================================================================
-c Arguments:
-c======================================================================
-cym#include "dimensions.h"
-cym#include "dimphy.h"
-c======================================================================
-C
-      INTEGER i , k , it
-      REAL tr(klon,klev,nbtr) , d_tr(klon,klev,nbtr)
-      REAL dtime
-      REAL tautr(nbtr)
-C
-c      WRITE(*,'(''PASSAGE radiornpb ... '',$)')
-C Attention, pour un pas de temps beaucoup plus petit que la decroissance!!!
-
-      DO it = 1,2
-           IF ( tautr(it) .GT. 0. ) THEN
-          	DO k = 1,klev
-                DO i = 1,klon
-                d_tr(i,k,it) = - tr(i,k,it) * dtime / tautr(it)
-                END DO
-                END DO
-           ELSE
-                DO k = 1,klev
-                DO i = 1,klon
-                d_tr(i,k,it) = 0.
-                END DO
-                END DO
-           END IF
-      END DO
-C
-CG161294 : Cas particulier radon 1 => plomb 2
-c
-      DO k = 1,klev
-        DO i = 1,klon
-          d_tr(i,k,2) = d_tr(i,k,2) - d_tr(i,k,1)
-        ENDDO
-      ENDDO
-c
-c      WRITE(*,*) ' radiornpb OK'
-c
-      RETURN
-      END
Index: /LMDZ4/branches/LMDZ4-dev/libf/phylmd/tracinca_mod.F90
===================================================================
--- /LMDZ4/branches/LMDZ4-dev/libf/phylmd/tracinca_mod.F90	(revision 1191)
+++ /LMDZ4/branches/LMDZ4-dev/libf/phylmd/tracinca_mod.F90	(revision 1191)
@@ -0,0 +1,191 @@
+!$Id $
+!
+MODULE tracinca_mod
+!
+! This module prepares and calls the INCA main subroutines. 
+!
+
+CONTAINS
+
+  SUBROUTINE tracinca_init(aerosol,lessivage)
+    ! This subroutine initialize some control varaibles. 
+
+    USE infotrac
+    IMPLICIT NONE
+    
+    ! Output variables
+    LOGICAL,DIMENSION(nbtr), INTENT(OUT) :: aerosol
+    LOGICAL,INTENT(OUT) :: lessivage
+    
+    
+    ! Initialization
+    lessivage  =.FALSE.
+    aerosol(:) = .FALSE.
+        
+  END SUBROUTINE tracinca_init
+
+  SUBROUTINE tracinca(                                &
+       nstep,    julien,   gmtime,         lafin,     &
+       pdtphys,  t_seri,   paprs,          pplay,     &
+       pmfu,     ftsol,    pctsrf,         pphis,     &
+       pphi,     albsol,   sh,             rh,        &
+       cldfra,   rneb,     diafra,         cldliq,    &
+       itop_con, ibas_con, pmflxr,         pmflxs,    &
+       prfl,     psfl,     aerosol_couple, flxmass_w, &
+       tau_aero, piz_aero, cg_aero,        ccm,       &
+       rfname,                                        &
+       tr_seri,  source,   solsym)      
+
+!========================================================
+!    -- CHIMIE INCA --
+!========================================================
+
+    USE dimphy
+    USE infotrac
+    USE vampir
+    USE comgeomphy
+    
+    IMPLICIT NONE
+    
+    INCLUDE "indicesol.h"
+    INCLUDE "control.h"
+    INCLUDE "dimensions.h"
+    INCLUDE "paramet.h"
+
+!==========================================================================
+!                   -- DESCRIPTION DES ARGUMENTS --
+!==========================================================================
+
+
+! EN ENTREE ...
+!
+!Configuration grille,temps:
+    INTEGER,INTENT(IN) :: nstep      ! Appel physique
+    INTEGER,INTENT(IN) :: julien     ! Jour julien
+    REAL,INTENT(IN)    :: gmtime
+    REAL,INTENT(IN)    :: pdtphys    ! Pas d'integration pour la physique (seconde)
+    LOGICAL,INTENT(IN) :: lafin      ! le flag de la fin de la physique
+    
+
+!Physique: 
+!--------
+    REAL,DIMENSION(klon,klev),INTENT(IN)   :: t_seri  ! Temperature
+    REAL,DIMENSION(klon,klev),INTENT(IN)   :: sh      ! humidite specifique
+    REAL,DIMENSION(klon,klev),INTENT(IN)   :: rh      ! humidite relative
+    REAL,DIMENSION(klon,klev+1),INTENT(IN) :: paprs   ! pression pour chaque inter-couche (en Pa)
+    REAL,DIMENSION(klon,klev),INTENT(IN)   :: pplay   ! pression pour le mileu de chaque couche (en Pa)
+    REAL,DIMENSION(klon,klev),INTENT(IN)   :: pphi    ! geopotentiel
+    REAL,DIMENSION(klon),INTENT(IN)        :: pphis
+    REAL,DIMENSION(klon,klev),INTENT(IN)   :: cldliq  ! eau liquide nuageuse
+    REAL,DIMENSION(klon,klev),INTENT(IN)   :: cldfra  ! fraction nuageuse (tous les nuages)
+    REAL,DIMENSION(klon,klev),INTENT(IN)   :: diafra  ! fraction nuageuse (convection ou stratus artificiels)
+    REAL,DIMENSION(klon,klev),INTENT(IN)   :: rneb    ! fraction nuageuse (grande echelle)
+    INTEGER,DIMENSION(klon),INTENT(IN)     :: itop_con
+    INTEGER,DIMENSION(klon),INTENT(IN)     :: ibas_con
+    REAL,DIMENSION(klon),INTENT(IN)        :: albsol  ! albedo surface
+!
+!Convection:
+!----------
+    REAL,DIMENSION(klon,klev),INTENT(IN) :: pmfu  ! flux de masse dans le panache montant
+
+!...Tiedke     
+    REAL,DIMENSION(klon,klev+1),INTENT(IN)   :: pmflxr, pmflxs ! Flux precipitant de pluie, neige aux interfaces [convection]
+    REAL,DIMENSION(klon,klev+1),INTENT(IN)   :: prfl, psfl ! Flux precipitant de pluie, neige aux interfaces [large-scale]
+
+    LOGICAL,INTENT(IN)                       :: aerosol_couple
+    REAL,DIMENSION(klon,klev),INTENT(IN)     :: flxmass_w
+    REAL,DIMENSION(klon,klev,9,2),INTENT(IN) :: tau_aero
+    REAL,DIMENSION(klon,klev,9,2),INTENT(IN) :: piz_aero
+    REAL,DIMENSION(klon,klev,9,2),INTENT(IN) :: cg_aero
+    CHARACTER(len=4),DIMENSION(9),INTENT(IN) :: rfname 
+    REAL,DIMENSION(klon,klev,2),INTENT(IN)   :: ccm 
+
+! Arguments necessaires pour les sources et puits de traceur:
+    REAL,DIMENSION(klon,nbsrf),INTENT(IN) :: ftsol  ! Temperature du sol (surf)(Kelvin)
+    REAL,DIMENSION(klon,nbsrf),INTENT(IN) :: pctsrf ! Pourcentage de sol f(nature du sol)
+
+
+  ! InOutput argument
+    REAL,DIMENSION(klon,klev,nbtr),INTENT(INOUT) :: tr_seri ! Concentration Traceur [U/KgA]  
+
+  ! Output arguments
+    REAL,DIMENSION(klon,nbtr), INTENT(OUT)        :: source  ! a voir lorsque le flux de surface est prescrit 
+    CHARACTER(len=8),DIMENSION(nbtr), INTENT(OUT) :: solsym
+
+!=======================================================================================
+!                        -- VARIABLES LOCALES TRACEURS --
+!=======================================================================================
+
+    INTEGER :: k
+    REAL,DIMENSION(klon,klev) :: pdel
+    REAL    :: calday
+    INTEGER :: ncsec
+
+    CALL VTe(VTphysiq)
+    CALL VTb(VTinca)
+    
+    calday = FLOAT(julien) + gmtime
+    ncsec  = NINT (86400.*gmtime)
+     
+    DO k = 1, klev
+       pdel(:,k) = paprs(:,k) - paprs (:,k+1)
+    END DO
+    
+    IF (config_inca == 'aero') THEN
+#ifdef INCA
+       CALL aerosolmain(                    &
+            aerosol_couple,tr_seri,pdtphys, &
+            pplay,pdel,prfl,pmflxr,psfl,    &
+            pmflxs,pmfu,itop_con,ibas_con,  &
+            pphi,airephy,nstep,rneb,t_seri, &      
+            rh,tau_aero,piz_aero,cg_aero,   &
+            rfname,ccm,lafin)
+#endif
+    END IF
+
+#ifdef INCA
+    CALL chemmain (tr_seri, &   !mmr
+         nstep,      & !nstep
+         calday,     & !calday
+         julien,     & !ncdate
+         ncsec,      & !ncsec
+         1,          & !lat
+         pdtphys,    & !delt
+         paprs(1,1), & !ps
+         pplay,      & !pmid
+         pdel,       & !pdel
+         airephy,    &
+         pctsrf(1,1),& !oro
+         ftsol,      & !tsurf
+         albsol,     & !albs
+         pphi,       & !zma
+         pphis,      & !phis
+         cldfra,     & !cldfr
+         rneb,       & !cldfr_st
+         diafra,     & !cldfr_cv
+         itop_con,   & !cldtop
+         ibas_con,   & !cldbot
+         cldliq,     & !cwat
+         prfl,       & !flxrst
+         pmflxr,     & !flxrcv
+         psfl,       & !flxsst
+         pmflxs,     & !flxscv
+         pmfu,       & !flxupd
+         flxmass_w,  & !flxmass_w
+         t_seri,     & !tfld
+         sh,         & !sh
+         rh,         & !rh
+         iip1,       & !nx
+         jjp1,       & !ny
+         source,     &
+         solsym)
+#endif
+    
+    CALL VTe(VTinca)
+    CALL VTb(VTphysiq)
+    
+    
+  END SUBROUTINE tracinca
+
+
+END MODULE tracinca_mod
Index: /LMDZ4/branches/LMDZ4-dev/libf/phylmd/traclmdz_mod.F90
===================================================================
--- /LMDZ4/branches/LMDZ4-dev/libf/phylmd/traclmdz_mod.F90	(revision 1191)
+++ /LMDZ4/branches/LMDZ4-dev/libf/phylmd/traclmdz_mod.F90	(revision 1191)
@@ -0,0 +1,329 @@
+!$Id $
+!
+MODULE traclmdz_mod
+! 
+! In this module all tracers specific to LMDZ are treated. This module is used 
+! only if running without any other chemestry model as INCA or REPROBUS.  
+!
+
+  REAL,DIMENSION(:,:),ALLOCATABLE,SAVE :: masktr   ! Masque reservoir de sol traceur
+!$OMP THREADPRIVATE(masktr)                        ! Masque de l'echange avec la surface (1 = reservoir) ou (possible >= 1 )
+  REAL,DIMENSION(:,:),ALLOCATABLE,SAVE :: fshtr    ! Flux surfacique dans le reservoir de sol
+!$OMP THREADPRIVATE(fshtr)
+  REAL,DIMENSION(:),ALLOCATABLE,SAVE   :: hsoltr   ! Epaisseur equivalente du reservoir de sol
+!$OMP THREADPRIVATE(hsoltr)
+!
+!Radioelements:
+!--------------
+!
+  REAL,DIMENSION(:),ALLOCATABLE,SAVE   :: tautr    ! Constante de decroissance radioactive
+!$OMP THREADPRIVATE(tautr)
+  REAL,DIMENSION(:),ALLOCATABLE,SAVE   :: vdeptr   ! Vitesse de depot sec dans la couche Brownienne
+!$OMP THREADPRIVATE(vdeptr)
+  REAL,DIMENSION(:),ALLOCATABLE,SAVE   :: scavtr   ! Coefficient de lessivage
+!$OMP THREADPRIVATE(scavtr)
+  REAL,DIMENSION(:,:),ALLOCATABLE,SAVE :: srcbe    ! Production du beryllium7 dans l atmosphere (U/s/kgA)
+!$OMP THREADPRIVATE(srcbe)
+
+  LOGICAL,DIMENSION(:),ALLOCATABLE,SAVE :: radio    ! radio(it)   = true  => decroisssance radioactive
+!$OMP THREADPRIVATE(radio)  
+
+  REAL,DIMENSION(:,:),ALLOCATABLE,SAVE :: trs     ! Conc. radon ds le sol
+!$OMP THREADPRIVATE(trs)
+
+  INTEGER,SAVE :: id_be       ! Activation et position du traceur Be7 [ id_be=0 -> desactive ] 
+!$OMP THREADPRIVATE(id_be)
+
+  LOGICAL,SAVE :: rnpb=.TRUE. ! Presence du couple Rn222, Pb210
+!$OMP THREADPRIVATE(rnpb)
+
+
+CONTAINS
+
+  SUBROUTINE traclmdz_from_restart(trs_in)
+    ! This subroutine initialize the module saved variable trs with values from restart file (startphy.nc). 
+    ! This subroutine is called from phyetat0 after the field trs_in has been read.
+    
+    USE dimphy
+    USE infotrac
+    IMPLICIT NONE
+    
+    ! Input argument
+    REAL,DIMENSION(klon,nbtr), INTENT(IN) :: trs_in 
+    
+    ! Local variables
+    INTEGER :: ierr
+    
+    ! Allocate restart variables trs
+    ALLOCATE( trs(klon,nbtr), stat=ierr)
+    IF (ierr /= 0) CALL abort_gcm('traclmdz_start', 'pb in allocation 1',1)
+    
+    ! Initialize trs with values read from restart file 
+    trs(:,:) = trs_in(:,:)
+    
+  END SUBROUTINE traclmdz_from_restart
+
+
+  SUBROUTINE traclmdz_init(pctsrf, ftsol, aerosol, lessivage)
+    ! This subroutine allocates and initialize module variables and control variables.
+    USE dimphy
+    USE infotrac
+
+    IMPLICIT NONE
+
+    INCLUDE "indicesol.h"
+
+! Input variables
+    REAL,DIMENSION(klon,nbsrf),INTENT(IN) :: pctsrf ! Pourcentage de sol f(nature du sol)
+    REAL,DIMENSION(klon,nbsrf),INTENT(IN) :: ftsol  ! Temperature du sol (surf)(Kelvin)
+    
+! Output variables
+    LOGICAL,DIMENSION(nbtr), INTENT(OUT) :: aerosol
+    LOGICAL,INTENT(OUT)                  :: lessivage
+        
+! Local variables    
+    INTEGER :: ierr, it, iiq
+    
+! --------------------------------------------
+! Allocation
+! --------------------------------------------
+
+    ALLOCATE( scavtr(nbtr), stat=ierr)
+    IF (ierr /= 0) CALL abort_gcm('phytrac', 'pb in allocation 9',1)
+    scavtr(:)=1.
+    
+    ALLOCATE( radio(nbtr), stat=ierr)
+    IF (ierr /= 0) CALL abort_gcm('phytrac', 'pb in allocation 11',1)
+    radio(:) = .false.    ! Par defaut pas decroissance radioactive
+    
+    ALLOCATE( masktr(klon,nbtr), stat=ierr)
+    IF (ierr /= 0) CALL abort_gcm('phytrac', 'pb in allocation 2',1)
+    
+    ALLOCATE( fshtr(klon,nbtr), stat=ierr)
+    IF (ierr /= 0) CALL abort_gcm('phytrac', 'pb in allocation 3',1)
+    
+    ALLOCATE( hsoltr(nbtr), stat=ierr)
+    IF (ierr /= 0) CALL abort_gcm('phytrac', 'pb in allocation 4',1)
+    
+    ALLOCATE( tautr(nbtr), stat=ierr)
+    IF (ierr /= 0) CALL abort_gcm('phytrac', 'pb in allocation 5',1)
+    tautr(:)  = 0.
+    
+    ALLOCATE( vdeptr(nbtr), stat=ierr)
+    IF (ierr /= 0) CALL abort_gcm('phytrac', 'pb in allocation 6',1)
+    vdeptr(:) = 0.
+
+
+    lessivage  = .TRUE.
+    aerosol(:) = .FALSE.  ! Tous les traceurs sont des gaz par defaut
+    
+!
+! Recherche des traceurs connus : Be7, CO2,...
+! --------------------------------------------
+    DO it=1,nbtr
+       iiq=niadv(it+2)
+       !
+       ! Recherche du Beryllium 7
+       !
+       IF ( tname(iiq) == "BE" .OR. tname(iiq) == "Be" .OR.  &
+            tname(iiq) == "BE7" .OR. tname(iiq) == "Be7" ) THEN  
+          id_be=it
+          ALLOCATE( srcbe(klon,klev) )
+          radio(id_be) = .TRUE.
+          aerosol(id_be) = .TRUE. ! le Be est un aerosol
+          CALL init_be(pctsrf,masktr(:,id_be),tautr(id_be),vdeptr(id_be),scavtr(id_be),srcbe)
+          WRITE(*,*) 'Initialisation srcBe: OK'
+       ELSE
+          id_be=0
+       END IF
+       
+       ! Placer ici la recherche de nouveaux traceurs
+       ! ...
+    END DO
+!
+! Valeurs specifiques pour les traceurs Rn222 et Pb210
+! ----------------------------------------------
+    IF (rnpb) THEN
+        
+       radio(1)= .TRUE.
+       radio(2)= .TRUE.
+       pbl_flg(1) = 0 ! au lieu de clsol=true ! CL au sol calcule
+       pbl_flg(2) = 0 ! au lieu de clsol=true
+       
+       aerosol(2) = .TRUE. ! le Pb est un aerosol
+       
+       CALL initrrnpb (ftsol,pctsrf,masktr,fshtr,hsoltr,tautr,vdeptr,scavtr)
+    END IF
+
+  END SUBROUTINE traclmdz_init
+
+  SUBROUTINE traclmdz(                           &
+       pdtphys,  t_seri,                         &
+       paprs,    pplay,        cdragh,  coefh,   &
+       yu1,      yv1,          ftsol,   pctsrf,  &
+       xlat,     couchelimite,                   &
+       tr_seri,  source,       solsym,  d_tr_cl)
+    
+    USE dimphy
+    USE infotrac
+    
+    IMPLICIT NONE
+    
+    INCLUDE "YOMCST.h"
+    INCLUDE "indicesol.h"
+
+!==========================================================================
+!                   -- DESCRIPTION DES ARGUMENTS --
+!==========================================================================
+
+! Input arguments
+!
+!Configuration grille,temps:
+    REAL,INTENT(IN)    :: pdtphys    ! Pas d'integration pour la physique (seconde)  
+    REAL,DIMENSION(klon),INTENT(IN) :: xlat    ! latitudes pour chaque point 
+
+!
+!Physique: 
+!--------
+    REAL,DIMENSION(klon,klev),INTENT(IN)   :: t_seri  ! Temperature
+    REAL,DIMENSION(klon,klev+1),INTENT(IN) :: paprs   ! pression pour chaque inter-couche (en Pa)
+    REAL,DIMENSION(klon,klev),INTENT(IN)   :: pplay   ! pression pour le mileu de chaque couche (en Pa)
+
+
+!Couche limite:
+!--------------
+!
+    REAL,DIMENSION(klon,klev),INTENT(IN) :: cdragh     ! coeff drag pour T et Q
+    REAL,DIMENSION(klon,klev),INTENT(IN) :: coefh      ! coeff melange CL (m**2/s)
+    REAL,DIMENSION(klon),INTENT(IN)      :: yu1        ! vents au premier niveau
+    REAL,DIMENSION(klon),INTENT(IN)      :: yv1        ! vents au premier niveau
+    LOGICAL,INTENT(IN)                   :: couchelimite
+
+! Arguments necessaires pour les sources et puits de traceur:
+    REAL,DIMENSION(klon,nbsrf),INTENT(IN) :: ftsol  ! Temperature du sol (surf)(Kelvin)
+    REAL,DIMENSION(klon,nbsrf),INTENT(IN) :: pctsrf ! Pourcentage de sol f(nature du sol)
+
+! InOutput argument
+    REAL,DIMENSION(klon,klev,nbtr),INTENT(INOUT)  :: tr_seri ! Concentration Traceur [U/KgA]  
+
+! Output argument
+    CHARACTER(len=8),DIMENSION(nbtr), INTENT(OUT) :: solsym
+    REAL,DIMENSION(klon,nbtr), INTENT(OUT)        :: source  ! a voir lorsque le flux de surface est prescrit 
+    REAL,DIMENSION(klon,klev,nbtr), INTENT(OUT)   :: d_tr_cl ! Td couche limite/traceur
+
+
+!=======================================================================================
+!                        -- VARIABLES LOCALES TRACEURS --
+!=======================================================================================
+
+    INTEGER :: i, k, it
+  
+
+    REAL,DIMENSION(klon)           :: d_trs    ! Td dans le reservoir
+    REAL,DIMENSION(klon,klev)      :: delp     ! epaisseur de couche (Pa)
+    
+    REAL,DIMENSION(klon,klev,nbtr) :: d_tr_dec ! Td radioactive
+    REAL                           :: zrho      ! Masse Volumique de l'air KgA/m3
+
+!
+!
+!=================================================================
+!  Ajout de la production en  Be7 (Beryllium) srcbe U/s/kgA
+!=================================================================
+!
+    IF ( id_be /= 0 ) THEN
+       DO k = 1, klev
+          DO i = 1, klon
+             tr_seri(i,k,id_be) = tr_seri(i,k,id_be)+srcbe(i,k)*pdtphys
+          END DO
+       END DO
+       WRITE(*,*) 'Ajout srcBe dans tr_seri: OK'
+    END IF
+  
+
+    DO it=1,nbtr
+       WRITE(solsym(it),'(i2)') it
+    END DO
+!======================================================================
+!     -- Calcul de l'effet de la couche limite --
+!======================================================================
+
+    IF (couchelimite) THEN             
+       DO it=1, nbtr
+          IF ( it == id_be ) THEN
+             DO i=1, klon
+                zrho = pplay(i,1)/t_seri(i,1)/RD
+                source(i,it) = - vdeptr(it)*tr_seri(i,1,it)*zrho
+             END DO
+          ELSE
+             source(:,it) = 0.0  
+          END IF
+       END DO
+    END IF
+    
+
+    DO k = 1, klev
+       DO i = 1, klon
+          delp(i,k) = paprs(i,k)-paprs(i,k+1)
+       END DO
+    END DO
+    
+    DO it=1, nbtr
+       IF (couchelimite .AND. pbl_flg(it) == 0 ) THEN ! couche limite avec quantite dans le sol calculee
+          CALL cltracrn(it, pdtphys, yu1, yv1,     &
+               cdragh, coefh,t_seri,ftsol,pctsrf,  &
+               tr_seri(:,:,it),trs(:,it),          &
+               paprs, pplay, delp,                 &
+               masktr(:,it),fshtr(:,it),hsoltr(it),&
+               tautr(it),vdeptr(it),               &
+               xlat,d_tr_cl(:,:,it),d_trs)
+          
+          DO k = 1, klev
+             DO i = 1, klon
+                tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_cl(i,k,it)
+             END DO
+          END DO
+        
+          ! Traceur dans le reservoir sol
+          DO i = 1, klon
+             trs(i,it) = trs(i,it) + d_trs(i)
+          END DO
+       END IF
+    END DO
+           
+!======================================================================
+!   Calcul de l'effet du puits radioactif
+!======================================================================
+    CALL radio_decay (radio,rnpb,pdtphys,tautr,tr_seri,d_tr_dec)
+  
+    DO it=1,nbtr
+       IF(radio(it)) then     
+          DO k = 1, klev
+             DO i = 1, klon
+                tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_dec(i,k,it)
+             END DO
+          END DO
+          CALL minmaxqfi(tr_seri(:,:,it),0.,1.e33,'puits rn it='//solsym(it))
+       END IF
+    END DO
+    
+
+  END SUBROUTINE traclmdz
+
+
+  SUBROUTINE traclmdz_to_restart(trs_out)
+    ! This subroutine is called from phyredem.F where the module 
+    ! variable trs is written to restart file (restartphy.nc)
+    USE dimphy
+    USE infotrac
+    
+    IMPLICIT NONE
+    
+    REAL,DIMENSION(klon,nbtr), INTENT(OUT) :: trs_out
+    
+    trs_out(:,:) = trs(:,:)
+    
+  END SUBROUTINE traclmdz_to_restart
+  
+
+END MODULE traclmdz_mod
Index: /LMDZ4/branches/LMDZ4-dev/libf/phylmd/write_histrac.h
===================================================================
--- /LMDZ4/branches/LMDZ4-dev/libf/phylmd/write_histrac.h	(revision 1190)
+++ /LMDZ4/branches/LMDZ4-dev/libf/phylmd/write_histrac.h	(revision 1191)
@@ -1,83 +1,85 @@
+!$Id $
+!***************************************
+!  ECRITURE DU FICHIER :  histrac.nc
+!***************************************
+#ifdef CPP_IOIPSL
+  IF (ecrit_tra > 0. .AND. config_inca == 'none') THEN
+     
+     itau_w = itau_phy + nstep
+     
+     CALL histwrite_phy(nid_tra,"phis",itau_w,pphis)
+     CALL histwrite_phy(nid_tra,"aire",itau_w,airephy)
+
+!TRACEURS
+!----------------
+     DO it=1,nbtr
+        iiq=niadv(it+2)
+
+! CONCENTRATIONS
+        CALL histwrite_phy(nid_tra,tname(iiq),itau_w,tr_seri(:,:,it))
+
+! TD LESSIVAGE       
+        IF (lessivage .AND. aerosol(it)) THEN
+           CALL histwrite_phy(nid_tra,"fl"//tname(iiq),itau_w,flestottr(:,:,it))
+        ENDIF
+
+! TD THERMIQUES
+        IF (iflag_thermals.gt.0) THEN
+           CALL histwrite_phy(nid_tra,"d_tr_th_"//tname(iiq),itau_w,d_tr_th(:,:,it))
+        ENDIF
+
+! TD CONVECTION
+        IF (iflag_con.GE.2) THEN
+           CALL histwrite_phy(nid_tra,"d_tr_cv_"//tname(iiq),itau_w,d_tr_cv(:,:,it))
+        ENDIF
+
+! TD COUCHE-LIMITE
+        CALL histwrite_phy(nid_tra,"d_tr_cl_"//tname(iiq),itau_w,d_tr_cl(:,:,it))
+     ENDDO
+!---------------
 !
-! $Header$
 !
+! VENT (niveau 1)   
+     CALL histwrite_phy(nid_tra,"pyu1",itau_w,yu1)
+     CALL histwrite_phy(nid_tra,"pyv1",itau_w,yv1)
+!
+! TEMPERATURE DU SOL
+     zx_tmp_fi2d(:)=ftsol(:,1)         
+     CALL histwrite_phy(nid_tra,"ftsol1",itau_w,zx_tmp_fi2d)
+     zx_tmp_fi2d(:)=ftsol(:,2)
+     CALL histwrite_phy(nid_tra,"ftsol2",itau_w,zx_tmp_fi2d)
+     zx_tmp_fi2d(:)=ftsol(:,3)
+     CALL histwrite_phy(nid_tra,"ftsol3",itau_w,zx_tmp_fi2d)
+     zx_tmp_fi2d(:)=ftsol(:,4)
+     CALL histwrite_phy(nid_tra,"ftsol4",itau_w,zx_tmp_fi2d)
+!      
+! NATURE DU SOL
+     zx_tmp_fi2d(:)=pctsrf(:,1)
+     CALL histwrite_phy(nid_tra,"psrf1",itau_w,zx_tmp_fi2d)
+     zx_tmp_fi2d(:)=pctsrf(:,2)
+     CALL histwrite_phy(nid_tra,"psrf2",itau_w,zx_tmp_fi2d)
+     zx_tmp_fi2d(:)=pctsrf(:,3)
+     CALL histwrite_phy(nid_tra,"psrf3",itau_w,zx_tmp_fi2d)
+     zx_tmp_fi2d(:)=pctsrf(:,4)
+     CALL histwrite_phy(nid_tra,"psrf4",itau_w,zx_tmp_fi2d)
+ 
+! DIVERS    
+     CALL histwrite_phy(nid_tra,"pplay",itau_w,pplay)     
+     CALL histwrite_phy(nid_tra,"t",itau_w,t_seri)     
+     CALL histwrite_phy(nid_tra,"mfu",itau_w,pmfu)
+     CALL histwrite_phy(nid_tra,"mfd",itau_w,pmfd)
+     CALL histwrite_phy(nid_tra,"en_u",itau_w,pen_u)
+     CALL histwrite_phy(nid_tra,"en_d",itau_w,pen_d)
+     CALL histwrite_phy(nid_tra,"de_d",itau_w,pde_d)
+     CALL histwrite_phy(nid_tra,"de_u",itau_w,pde_u)
+     CALL histwrite_phy(nid_tra,"coefh",itau_w,coefh)
 
-      IF (ecrit_tra>0. .AND. config_inca == 'none') THEN
-      ndex = 0
-      ndex2d = 0
-      ndex3d = 0
-c
-      itau_w = itau_phy + nstep
+     IF (ok_sync) THEN
+!$OMP MASTER
+        CALL histsync(nid_tra)
+!$OMP END MASTER
+     ENDIF
 
-      CALL histwrite_phy(nid_tra,"phis",itau_w,pphis)
-C
-      CALL histwrite_phy(nid_tra,"aire",itau_w,airephy)
+  ENDIF !ecrit_tra>0. .AND. config_inca == 'none'
 
-      DO it=1,nbtr
-       iiq=niadv(it+2)
-
-       CALL histwrite_phy(nid_tra,tname(iiq),itau_w,tr_seri(:,:,it))
-       if (lessivage) THEN
-       CALL histwrite_phy(nid_tra,"fl"//tname(iiq),itau_w,
-     .                                   flestottr(:,:,it))
-      endif
-      
-c----Olivia
-       CALL histwrite_phy(nid_tra,"d_tr_th_"//tname(iiq),itau_w,
-     .                                           d_tr_th(:,:,it))
-
-         if(iflag_con.GE.2) then
-       CALL histwrite_phy(nid_tra,"d_tr_cv_"//tname(iiq),itau_w,
-     .                                           d_tr_cv(:,:,it))
-         endif !(iflag_con.GE.2) then
-       CALL histwrite_phy(nid_tra,"d_tr_cl_"//tname(iiq),itau_w,
-     .                                           d_tr_cl(:,:,it))
-c---fin Olivia      
-      
-      ENDDO
-
-
-C abder
-         CALL histwrite_phy(nid_tra,"pyu1",itau_w,yu1)
-
-         CALL histwrite_phy(nid_tra,"pyv1",itau_w,yv1)
-
-         CALL histwrite_phy(nid_tra,"ftsol1",itau_w,pftsol1)
-
-         CALL histwrite_phy(nid_tra,"ftsol2",itau_w,pftsol2)
-
-         CALL histwrite_phy(nid_tra,"ftsol3",itau_w,pftsol3)
-
-         CALL histwrite_phy(nid_tra,"ftsol4",itau_w,pftsol4)
-
-         CALL histwrite_phy(nid_tra,"psrf1",itau_w,ppsrf1)
-
-         CALL histwrite_phy(nid_tra,"psrf2",itau_w,ppsrf2)
-
-         CALL histwrite_phy(nid_tra,"psrf3",itau_w,ppsrf3)
-
-         CALL histwrite_phy(nid_tra,"psrf4",itau_w,ppsrf4)
-        CALL histwrite_phy(nid_tra,"pplay",itau_w,pplay)
-
-        CALL histwrite_phy(nid_tra,"t",itau_w,t_seri)
-        CALL histwrite_phy(nid_tra,"mfu",itau_w,pmfu)
-        CALL histwrite_phy(nid_tra,"mfd",itau_w,pmfd)
-        CALL histwrite_phy(nid_tra,"en_u",itau_w,pen_u)
-        CALL histwrite_phy(nid_tra,"en_d",itau_w,pen_d)
-        CALL histwrite_phy(nid_tra,"de_d",itau_w,pde_d)
-        CALL histwrite_phy(nid_tra,"de_u",itau_w,pde_u)
-        CALL histwrite_phy(nid_tra,"coefh",itau_w,coefh)
-
-
-c abder
-
-      if (ok_sync) then
-c$OMP MASTER
-        call histsync(nid_tra)
-c$OMP END MASTER
-       endif
-
-       END IF !ecrit_tra>0. .AND. config_inca == 'none'
-
-
-
+#endif
