Index: LMDZ4/trunk/libf/phylmd/traclmdz_mod.F90
===================================================================
--- LMDZ4/trunk/libf/phylmd/traclmdz_mod.F90	(revision 1409)
+++ LMDZ4/trunk/libf/phylmd/traclmdz_mod.F90	(revision 1421)
@@ -6,5 +6,4 @@
 ! only if running without any other chemestry model as INCA or REPROBUS.  
 !
-
   IMPLICIT NONE
 
@@ -45,17 +44,14 @@
 !$OMP THREADPRIVATE(id_be)
 
-!IM ajout traceurs RR 
-  INTEGER,SAVE :: id_dry !traceur dry intrusions
-!$OMP THREADPRIVATE(id_dry)
   INTEGER,SAVE :: id_pcsat, id_pcocsat, id_pcq ! traceurs pseudo-vapeur CL qsat, qsat_oc, q
 !$OMP THREADPRIVATE(id_pcsat, id_pcocsat, id_pcq)
-  INTEGER,SAVE :: id_pcs0, id_pcos0, id_pcq0 ! traceurs pseudo-vapeur CL qsat, qsat_oc, q
-!                                            ! qui ne sont pas transportes par la convection
+  INTEGER,SAVE :: id_pcs0, id_pcos0, id_pcq0   ! traceurs pseudo-vapeur CL qsat, qsat_oc, q
+!                                              ! qui ne sont pas transportes par la convection
 !$OMP THREADPRIVATE(id_pcs0, id_pcos0, id_pcq0)
 
   INTEGER, SAVE:: id_o3
-  !$OMP THREADPRIVATE(id_o3)
-  ! index of ozone tracer with Cariolle parameterization
-  ! 0 means no ozone tracer
+!$OMP THREADPRIVATE(id_o3)
+! index of ozone tracer with Cariolle parameterization
+! 0 means no ozone tracer
 
   LOGICAL,SAVE :: rnpb=.FALSE. ! Presence du couple Rn222, Pb210
@@ -90,4 +86,5 @@
   SUBROUTINE traclmdz_init(pctsrf, ftsol, tr_seri, t_seri, pplay, sh, aerosol, lessivage)
     ! This subroutine allocates and initialize module variables and control variables.
+    ! Initialization of the tracers should be done here only for those not found in the restart file.
     USE dimphy
     USE infotrac
@@ -95,11 +92,12 @@
     USE press_coefoz_m, ONLY: press_coefoz
     USE carbon_cycle_mod, ONLY : carbon_cycle_init, carbon_cycle_tr, carbon_cycle_cpl
+    USE mod_grid_phy_lmdz
+    USE mod_phys_lmdz_para
 
     INCLUDE "indicesol.h"
-
+    INCLUDE "iniprint.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)
-!IM traceurs RR   REAL,DIMENSION(klon,klev,nbtr),INTENT(IN) :: tr_seri! Concentration Traceur [U/KgA]  
     REAL,DIMENSION(klon,klev,nbtr),INTENT(INOUT) :: tr_seri! Concentration Traceur [U/KgA]  
     REAL,DIMENSION(klon,klev),INTENT(IN)   :: t_seri  ! Temperature
@@ -113,10 +111,11 @@
 ! Local variables    
     INTEGER :: ierr, it, iiq, i, k
-    REAL,DIMENSION(klon,klev)      :: qsat   ! pression de la vapeur a saturation
-    
+    REAL, DIMENSION(klon_glo,klev) :: varglo ! variable temporaire sur la grille global    
+    REAL, DIMENSION(klev)          :: mintmp, maxtmp
+    LOGICAL                        :: zero
+
 ! --------------------------------------------
 ! Allocation
 ! --------------------------------------------
-
     ALLOCATE( scavtr(nbtr), stat=ierr)
     IF (ierr /= 0) CALL abort_gcm('traclmdz_init', 'pb in allocation 9',1)
@@ -152,4 +151,5 @@
 ! --------------------------------------------
     id_rn=0; id_pb=0; id_aga=0; id_be=0; id_o3=0
+    id_pcsat=0; id_pcocsat=0; id_pcq=0; id_pcs0=0; id_pcos0=0; id_pcq0=0
     DO it=1,nbtr
        iiq=niadv(it+2)
@@ -181,5 +181,5 @@
           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'
+          WRITE(lunout,*) 'Initialisation srcBe: OK'
        ELSE IF (tname(iiq)=="O3" .OR. tname(iiq)=="o3") THEN
           ! Recherche de l'ozone : parametrization de la chimie par Cariolle
@@ -187,60 +187,21 @@
           CALL alloc_coefoz   ! allocate ozone coefficients
           CALL press_coefoz   ! read input pressure levels
-       END IF    
-    END DO
-
-    id_dry=0
-    DO it=1,nbtr
-       iiq=niadv(it+2)
-       IF ( tname(iiq) == "dry" .OR. tname(iiq) == "Dry" ) THEN
-        id_dry=it
-       END IF    
-    END DO  
-
-    id_pcsat=0
-    DO it=1,nbtr
-       iiq=niadv(it+2)
-       IF ( tname(iiq) == "pcsat" .OR. tname(iiq) == "Pcsat" ) THEN
-        id_pcsat=it
-      END IF
-    END DO
-
-    id_pcocsat=0
-    DO it=1,nbtr
-       iiq=niadv(it+2)
-       IF ( tname(iiq) == "pcocsat" .OR. tname(iiq) == "Pcocsat" ) THEN
-        id_pcocsat=it
-       END IF
-    END DO
-
-    id_pcq=0
-    DO it=1,nbtr
-       iiq=niadv(it+2)
-       IF ( tname(iiq) == "pcq" .OR. tname(iiq) == "Pcq" ) THEN
-        id_pcq=it
-       END IF
-    END DO
-
-    id_pcs0=0
-    DO it=1,nbtr
-       iiq=niadv(it+2)
-       IF ( tname(iiq) == "pcs0" .OR. tname(iiq) == "Pcs0" ) THEN
-        id_pcs0=it
-      END IF
-    END DO
-
-    id_pcos0=0
-    DO it=1,nbtr
-       iiq=niadv(it+2)
-       IF ( tname(iiq) == "pcos0" .OR. tname(iiq) == "Pcos0" ) THEN
-        id_pcos0=it
-       END IF
-    END DO
-
-    id_pcq0=0
-    DO it=1,nbtr
-       iiq=niadv(it+2)
-       IF ( tname(iiq) == "pcq0" .OR. tname(iiq) == "Pcq0" ) THEN
-        id_pcq0=it
+       ELSE IF ( tname(iiq) == "pcsat" .OR. tname(iiq) == "Pcsat" ) THEN
+          id_pcsat=it
+       ELSE IF ( tname(iiq) == "pcocsat" .OR. tname(iiq) == "Pcocsat" ) THEN
+          id_pcocsat=it
+       ELSE IF ( tname(iiq) == "pcq" .OR. tname(iiq) == "Pcq" ) THEN
+          id_pcq=it
+       ELSE IF ( tname(iiq) == "pcs0" .OR. tname(iiq) == "Pcs0" ) THEN
+          id_pcs0=it
+          conv_flg(it)=0 ! No transport by convection for this tracer
+       ELSE IF ( tname(iiq) == "pcos0" .OR. tname(iiq) == "Pcos0" ) THEN
+          id_pcos0=it
+          conv_flg(it)=0 ! No transport by convection for this tracer
+       ELSE IF ( tname(iiq) == "pcq0" .OR. tname(iiq) == "Pcq0" ) THEN
+          id_pcq0=it
+          conv_flg(it)=0 ! No transport by convection for this tracer
+       ELSE
+          WRITE(lunout,*) 'This is an unknown tracer in LMDZ : ', trim(tname(iiq))
        END IF
     END DO
@@ -268,82 +229,47 @@
     END IF
 
-!IM initialisation traceurs pseudo-vapeurs
-    call q_sat(klon*klev,t_seri,pplay,qsat)
-    IF ( id_pcsat /= 0 ) THEN
-     DO k = 1, klev
-      DO i = 1, klon
-       IF ( pplay(i,k).GE.85000.) THEN
-        tr_seri(i,k,id_pcsat) = qsat(i,k)
-       ELSE
-        tr_seri(i,k,id_pcsat) = 100.
-       END IF
-      END DO
-     END DO
-    END IF
-
-    IF ( id_pcocsat /= 0 ) THEN
-     DO k = 1, klev
-      DO i = 1, klon
-       IF ( pplay(i,k).GE.85000.) THEN
-        IF ( pctsrf (i, is_oce) > 0. ) THEN
-         tr_seri(i,k,id_pcocsat) = qsat(i,k)
-        ELSE
-         tr_seri(i,k,id_pcocsat) = 100.
-        END IF
-       END IF
-      END DO
-     END DO
-    END IF
-
-    IF ( id_pcq /= 0 ) THEN
-     DO k = 1, klev
-      DO i = 1, klon
-       IF ( pplay(i,k).GE.85000.) THEN
-        tr_seri(i,k,id_pcq) = sh(i,k)
-       ELSE
-        tr_seri(i,k,id_pcq) = 100.
-       END IF
-      END DO
-     END DO
-    END IF
-
-    IF ( id_pcs0 /= 0 ) THEN
-     DO k = 1, klev
-      DO i = 1, klon
-       IF ( pplay(i,k).GE.85000.) THEN
-        tr_seri(i,k,id_pcs0) = qsat(i,k)
-       ELSE
-        tr_seri(i,k,id_pcs0) = 100.
-       END IF
-      END DO
-     END DO
-    END IF
-
-    IF ( id_pcos0 /= 0 ) THEN
-     DO k = 1, klev
-      DO i = 1, klon
-       IF ( pplay(i,k).GE.85000.) THEN
-        IF ( pctsrf (i, is_oce) > 0. ) THEN
-         tr_seri(i,k,id_pcos0) = qsat(i,k)
-        ELSE
-         tr_seri(i,k,id_pcos0) = 100.
-        END IF
-       END IF
-      END DO
-     END DO
-    END IF
-
-    IF ( id_pcq0 /= 0 ) THEN
-     DO k = 1, klev
-      DO i = 1, klon
-       IF ( pplay(i,k).GE.85000.) THEN
-        tr_seri(i,k,id_pcq0) = sh(i,k)
-       ELSE
-        tr_seri(i,k,id_pcq0) = 100.
-       END IF
-      END DO
-     END DO
-    END IF
- 
+! Check if all tracers have restart values
+! ----------------------------------------------
+    DO it=1,nbtr
+       iiq=niadv(it+2)
+       ! Test if tracer is zero everywhere. 
+       ! Done by master process MPI and master thread OpenMP
+       CALL gather(tr_seri(:,:,it),varglo)
+       IF (is_mpi_root .AND. is_omp_root) THEN
+          mintmp=MINVAL(varglo,dim=1)
+          maxtmp=MAXVAL(varglo,dim=1)
+          IF (MINVAL(mintmp,dim=1)==0. .AND. MAXVAL(maxtmp,dim=1)==0.) THEN
+             ! Tracer is zero everywhere
+             zero=.TRUE.
+          ELSE
+             zero=.FALSE.
+          END IF
+       END IF
+
+       ! Distribute variable at all processes
+       CALL bcast(zero)
+
+       ! Initalize tracer that was not found in restart file.
+       IF (zero) THEN
+          ! The tracer was not found in restart file or it was equal zero everywhere.
+          WRITE(lunout,*) "The tracer ",trim(tname(iiq))," will be initialized"
+          IF (it==id_pcsat .OR. it==id_pcq .OR. &
+               it==id_pcs0 .OR. it==id_pcq0) THEN
+             tr_seri(:,:,it) = 100.
+          ELSE IF (it==id_pcocsat .OR. it==id_pcos0) THEN
+             DO i = 1, klon
+                IF ( pctsrf (i, is_oce) == 0. ) THEN
+                   tr_seri(i,:,it) = 0.
+                ELSE
+                   tr_seri(i,:,it) = 100.
+                END IF
+             END DO
+          ELSE
+             ! No specific initialization exist for this tracer
+             tr_seri(:,:,it) = 0.
+          END IF
+       END IF
+    END DO
+
   END SUBROUTINE traclmdz_init
 
@@ -391,11 +317,9 @@
     REAL,DIMENSION(klon),INTENT(IN)      :: yv1        ! vents au premier niveau
     LOGICAL,INTENT(IN)                   :: couchelimite
-!IM traceurs RR
-    REAL,DIMENSION(klon,klev),INTENT(IN)   :: sh      ! humidite specifique
+    REAL,DIMENSION(klon,klev),INTENT(IN) :: sh         ! humidite specifique
 
 ! 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
@@ -412,13 +336,11 @@
 
     INTEGER :: i, k, it
-    INTEGER lmt_pas ! number of time steps of "physics" per day
+    INTEGER :: lmt_pas ! number of time steps of "physics" per day
 
     REAL,DIMENSION(klon)           :: d_trs    ! Td dans le reservoir
+    REAL,DIMENSION(klon,klev)      :: qsat     ! pression de la vapeur a saturation
     REAL,DIMENSION(klon,klev,nbtr) :: d_tr_dec ! Td radioactive
-    REAL                           :: zrho      ! Masse Volumique de l'air KgA/m3
-
-!IM traceurs RR
-    REAL,DIMENSION(klon,klev)      :: qsat   ! pression de la vapeur a saturation
-    REAL :: amn, amx
+    REAL                           :: zrho     ! Masse Volumique de l'air KgA/m3
+    REAL                           :: amn, amx
 !
 !=================================================================
@@ -435,13 +357,19 @@
     END IF
   
-!IM ajout traceurs RR
-    call q_sat(klon*klev,t_seri,pplay,qsat)
-    
+
+!=================================================================
+! Update pseudo-vapor tracers 
+!=================================================================
+
+    CALL q_sat(klon*klev,t_seri,pplay,qsat)
+
     IF ( id_pcsat /= 0 ) THEN
      DO k = 1, klev
       DO i = 1, klon
-       IF ( pplay(i,k).GE.85000.) THEN
-        tr_seri(i,k,id_pcsat) = qsat(i,k)
-       END IF
+         IF ( pplay(i,k).GE.85000.) THEN
+            tr_seri(i,k,id_pcsat) = qsat(i,k)
+         ELSE
+            tr_seri(i,k,id_pcsat) = MIN (qsat(i,k), tr_seri(i,k,id_pcsat))            
+         END IF
       END DO 
      END DO
@@ -451,8 +379,12 @@
      DO k = 1, klev
       DO i = 1, klon
-       IF ( pplay(i,k).GE.85000.) THEN
-        IF ( pctsrf (i, is_oce) > 0. ) THEN
-         tr_seri(i,k,id_pcocsat) = qsat(i,k)
-        END IF
+         IF ( pplay(i,k).GE.85000.) THEN
+            IF ( pctsrf (i, is_oce) > 0. ) THEN
+               tr_seri(i,k,id_pcocsat) = qsat(i,k)
+            ELSE
+               tr_seri(i,k,id_pcocsat) = 0.
+          END IF
+       ELSE 
+          tr_seri(i,k,id_pcocsat) = MIN (qsat(i,k), tr_seri(i,k,id_pcocsat))
        END IF
       END DO
@@ -463,46 +395,55 @@
      DO k = 1, klev
       DO i = 1, klon
-       IF ( pplay(i,k).GE.85000.) THEN
-        tr_seri(i,k,id_pcq) = sh(i,k)
-       END IF
+         IF ( pplay(i,k).GE.85000.) THEN
+            tr_seri(i,k,id_pcq) = sh(i,k)
+         ELSE
+            tr_seri(i,k,id_pcq) = MIN (qsat(i,k), tr_seri(i,k,id_pcq))    
+         END IF
       END DO
      END DO
     END IF
+
 
     IF ( id_pcs0 /= 0 ) THEN
      DO k = 1, klev
       DO i = 1, klon
-       IF ( pplay(i,k).GE.85000.) THEN
-        tr_seri(i,k,id_pcs0) = qsat(i,k)
-       END IF
+         IF ( pplay(i,k).GE.85000.) THEN
+            tr_seri(i,k,id_pcs0) = qsat(i,k)
+         ELSE
+            tr_seri(i,k,id_pcs0) = MIN (qsat(i,k), tr_seri(i,k,id_pcs0))    
+         END IF
       END DO
      END DO
     END IF
+
 
     IF ( id_pcos0 /= 0 ) THEN
      DO k = 1, klev
       DO i = 1, klon
-       IF ( pplay(i,k).GE.85000.) THEN
-        IF ( pctsrf (i, is_oce) > 0. ) THEN
-         tr_seri(i,k,id_pcos0) = qsat(i,k)
-        END IF
-       END IF
+         IF ( pplay(i,k).GE.85000.) THEN
+            IF ( pctsrf (i, is_oce) > 0. ) THEN
+               tr_seri(i,k,id_pcos0) = qsat(i,k)
+            ELSE
+               tr_seri(i,k,id_pcos0) = 0.
+            END IF
+         ELSE
+            tr_seri(i,k,id_pcos0) = MIN (qsat(i,k), tr_seri(i,k,id_pcos0))
+         END IF
       END DO
      END DO
     END IF
+
 
     IF ( id_pcq0 /= 0 ) THEN
      DO k = 1, klev
       DO i = 1, klon
-       IF ( pplay(i,k).GE.85000.) THEN
-        tr_seri(i,k,id_pcq0) = sh(i,k)
-       END IF
+         IF ( pplay(i,k).GE.85000.) THEN
+            tr_seri(i,k,id_pcq0) = sh(i,k)
+         ELSE 
+            tr_seri(i,k,id_pcq0) = MIN (qsat(i,k), tr_seri(i,k,id_pcq0))
+         END IF
       END DO
      END DO
     END IF
-
-    DO it=1,nbtr
-       WRITE(solsym(it),'(i2)') it
-    END DO
 
 !=================================================================
@@ -566,69 +507,14 @@
     END DO
           
-!IM traceurs RR
-    IF ( id_pcsat /= 0 ) THEN
-     DO k = 1, klev
-      DO i = 1, klon
-       IF ( pplay(i,k).LT.85000.) THEN
-        tr_seri(i,k,id_pcsat) = MIN (qsat(i,k), tr_seri(i,k,id_pcsat))
-       END IF 
-      END DO
-     END DO
-    END IF
-
-    IF ( id_pcocsat /= 0 ) THEN
-     DO k = 1, klev
-      DO i = 1, klon
-       IF ( pplay(i,k).LT.85000.) THEN
-        tr_seri(i,k,id_pcocsat) = MIN (qsat(i,k), tr_seri(i,k,id_pcocsat))
-       END IF
-      END DO
-     END DO
-    END IF
-
-    IF ( id_pcq /= 0 ) THEN
-     DO k = 1, klev 
-      DO i = 1, klon 
-       IF ( pplay(i,k).LT.85000.) THEN
-        tr_seri(i,k,id_pcq) = MIN (qsat(i,k), tr_seri(i,k,id_pcq))
-       END IF
-      END DO  
-     END DO  
-    END IF  
-
-    IF ( id_pcs0 /= 0 ) THEN
-     DO k = 1, klev
-      DO i = 1, klon
-       IF ( pplay(i,k).LT.85000.) THEN
-        tr_seri(i,k,id_pcs0) = MIN (qsat(i,k), tr_seri(i,k,id_pcs0))
-       END IF
-      END DO
-     END DO
-    END IF
-
-    IF ( id_pcos0 /= 0 ) THEN
-     DO k = 1, klev
-      DO i = 1, klon
-       IF ( pplay(i,k).LT.85000.) THEN
-        tr_seri(i,k,id_pcos0) = MIN (qsat(i,k), tr_seri(i,k,id_pcos0))
-       END IF
-      END DO
-     END DO
-    END IF
-
-    IF ( id_pcq0 /= 0 ) THEN
-     DO k = 1, klev
-      DO i = 1, klon
-       IF ( pplay(i,k).LT.85000.) THEN
-        tr_seri(i,k,id_pcq0) = MIN (qsat(i,k), tr_seri(i,k,id_pcq0))
-       END IF
-      END DO
-     END DO
-    END IF
+
 !======================================================================
 !   Calcul de l'effet du puits radioactif
 !======================================================================
     CALL radio_decay (radio,rnpb,pdtphys,tautr,tr_seri,d_tr_dec)
-  
+
+    DO it=1,nbtr
+       WRITE(solsym(it),'(i2)') it
+    END DO
+
     DO it=1,nbtr
        IF(radio(it)) then     
