Index: LMDZ5/branches/testing/libf/phylmd/init_be.F90
===================================================================
--- LMDZ5/branches/testing/libf/phylmd/init_be.F90	(revision 1663)
+++ LMDZ5/branches/testing/libf/phylmd/init_be.F90	(revision 1750)
@@ -1,5 +1,6 @@
 !$Id $
 
-SUBROUTINE init_be(pctsrf,masktr,tautr,vdeptr,scavtr,srcbe)
+SUBROUTINE init_be(pctsrf,pplay,masktr,tautr,vdeptr,scavtr,srcbe)
+!!!SUBROUTINE init_be(pctsrf,masktr,tautr,vdeptr,scavtr,srcbe)
 
   USE dimphy
@@ -26,4 +27,5 @@
 !
   REAL,DIMENSION(klon,nbsrf),INTENT(IN) :: pctsrf !Pourcentage de sol (f(nature du sol))
+  REAL,DIMENSION(klon,klev), INTENT(IN) :: pplay  ! Pressions en milieu de couches
 !
 ! Output Arguments
@@ -37,15 +39,32 @@
 ! Local Variables
 !
+!!!  INTEGER              :: iref      ! numero d'un point oceanique donnant la grille de pression de reference
   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
+  INTEGER              :: k,i, kref, k2
+  INTEGER              :: nref
+  PARAMETER (nref=39)
+  REAL,DIMENSION(nref) :: pref      ! grille de pression de reference (bas des couches)
+  DATA pref  /   &
+      101249.99999999994, 100387.17261011522, 99447.35334189111,  98357.43412194174,   &
+      97046.47707771382,  95447.1116450629,   93496.85259615642,  91139.46548240296,   &
+      88326.55568744117,  85019.60710580258,  81192.7404556645,   76836.48366938648,   &
+      71962.81275769137,  66611.56331321516,  60857.914829743604, 54819.84484441629,   &
+      48663.06257114699,  42598.95465845692,  36869.104365898806, 31709.927925633147,  &
+      27296.757208636915, 23682.282929080895, 20766.025578936627, 18336.105961406534,  &
+      16178.04816768436,  14168.286905562818, 12275.719926478887, 10507.798835225762,  &
+      8876.585404909414,  7391.283929569539,  6057.514475749798,  4877.165909157005,   &
+      3848.34936408203,   2965.444753540027,  2219.2391544640013, 1597.15366044666,    &
+      1083.5531161631498, 660.1311067852655,  306.36072267002805 / 
+!$OMP THREADPRIVATE(pref)
 
   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)
-!
+! la source est maintenant définie independemment de la valeur de klev.
+!!! 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
 ! -------------------------
@@ -53,4 +72,18 @@
   vdeptr = 1.E-3 
   scavtr = 0.5 
+!!!!!jyg le 13/03/2013; puis 20/03/2013 : pref est maintenant une table.
+!!!
+!!!   Recherche d'un point rlat=0., rlon=180.
+!!      iref=(klon+1)/2
+!!      DO i = 1,klon
+!!        IF (abs(rlatd(i)) .LT. 0.15 .AND. cos(rlond(i)) .LT. -0.85) iref=i
+!!      ENDDO
+!!!
+!!!   Grille de pression de reference (= approx de sommets de couches)
+!!      pref(1) = pplay(iref,1)+0.5*(pplay(iref,1)-pplay(iref,2))
+!!      DO k = 2,klev
+!!        pref(k) = 0.5*(pplay(iref,k-1)+pplay(iref,k))
+!!      ENDDO
+!!!
 
   WRITE(*,*) '-------------- SOURCE DE BERYLLIUM ------------------- '
@@ -77,10 +110,13 @@
 ! 3-mettre la source de Be ds la bonne unite (en at/kgA/s)
 !
-  glt=78.5*rpi/180.
-  glg=69.0*rpi/180.
+  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)
+!!jyg
+!!     qcos=qcos+cos(glt)*cos(rlatd(i))*cos(rlond(i)+glg)
+     qcos=qcos+cos(glt)*cos(rlatd(i))*cos(rlond(i)-glg)
+!!jyg end
      IF ( qcos .LT. -1.) qcos = -1.
      IF ( qcos .GT. 1.)  qcos = 1.
@@ -88,23 +124,184 @@
   ENDDO
 
-!===========================
-!  Cas 19 niveaux verticaux
-!===========================
-  IF (klev.eq.19) then
+!!!===========================
+!!!  Cas 19 niveaux verticaux
+!!!===========================
+!!  IF (klev.eq.19) then
+!!     DO k = 1,klev
+!!        DO i = 1,klon
+!!!!!jyg le 13/03/2013
+!!!
+!!! k est le niveau dans la grille locale
+!!! Determination du niveau kref dans la grille de refernce
+!!      kref = 1
+!!      DO k2 = 1,klev
+!!        IF (pref(k2) .GT. pplay(i,k)) kref=k2
+!!      ENDDO
+!!!!!
+!!           latgeo=(180./rpi)*abs(rlatgeo(i))
+!!           IF ( kref .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 ( kref .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 ( kref .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 ( kref .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 ( kref .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 ( kref .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 ( kref .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 ( kref .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 ( kref .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 ( kref .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 ( kref .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 ( kref .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 ( kref .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 ( kref .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 ( kref .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 ( kref .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 ( kref .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
+!!
+!!!!!!  IF (klev .ge. 39) then
      DO k = 1,klev
         DO i = 1,klon
+!!!jyg le 13/03/2013
+!
+! k est le niveau dans la grille locale
+! Determination du niveau kref dans la grille de refernce
+      kref = 1
+      DO k2 = 1,nref
+        IF (pref(k2) .GT. pplay(i,k)) kref=k2
+      ENDDO
+!!!
            latgeo=(180./rpi)*abs(rlatgeo(i))
-           IF ( k .EQ. 1 ) THEN
+           IF ( kref .LE. 4 ) THEN
+              IF (latgeo.GE.50.0) srcbe(i,k)=0.07
+           END IF
+           IF ( kref .EQ. 5 ) 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.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. 3 ) THEN
+           IF ( kref .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
@@ -113,151 +310,5 @@
               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 ( kref .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
@@ -267,5 +318,5 @@
               IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=0.09
            END IF
-           IF ( k .EQ. 8 ) THEN
+           IF ( kref .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
@@ -275,5 +326,5 @@
               IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=0.1
            END IF
-           IF ( k .EQ. 9 ) THEN
+           IF ( kref .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
@@ -283,5 +334,5 @@
               IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=0.12
            END IF
-           IF ( k .EQ. 10 ) THEN
+           IF ( kref .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
@@ -291,5 +342,5 @@
               IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=0.14
            END IF
-           IF ( k .EQ. 11 ) THEN
+           IF ( kref .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
@@ -299,5 +350,5 @@
               IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=0.16
            END IF
-           IF ( k .EQ. 12 ) THEN
+           IF ( kref .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
@@ -306,5 +357,5 @@
               IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=0.25
            END IF
-           IF ( k .EQ. 13 ) THEN
+           IF ( kref .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
@@ -314,5 +365,5 @@
               IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=0.35
            END IF
-           IF ( k .EQ. 14 ) THEN
+           IF ( kref .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.
@@ -322,5 +373,5 @@
               IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=0.4
            END IF
-           IF ( k .EQ. 15 ) THEN
+           IF ( kref .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 
@@ -331,5 +382,5 @@
               IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=0.65
            END IF
-           IF ( k .EQ. 16 ) THEN
+           IF ( kref .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
@@ -339,5 +390,5 @@
               IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=0.9
            END IF
-           IF ( k .EQ. 17 ) THEN
+           IF ( kref .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.
@@ -347,5 +398,5 @@
               IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=1.4
            END IF
-           IF ( k .EQ. 18 ) THEN
+           IF ( kref .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.
@@ -355,5 +406,5 @@
               IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=2.
            END IF
-           IF ( k .EQ. 19 ) THEN
+           IF ( kref .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.
@@ -363,5 +414,5 @@
               IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=3.
            END IF
-           IF ( k .EQ. 20 ) THEN
+           IF ( kref .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.
@@ -371,5 +422,5 @@
               IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=3.5
            END IF
-           IF ( k .EQ. 21 ) THEN
+           IF ( kref .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.
@@ -379,5 +430,5 @@
               IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=4.
            END IF
-           IF ( k .EQ. 22 ) THEN
+           IF ( kref .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
@@ -387,5 +438,5 @@
               IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=4.5
            END IF
-           IF ( k .EQ. 23 ) THEN
+           IF ( kref .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.
@@ -395,5 +446,5 @@
               IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=6.
            END IF
-           IF ( k .EQ. 24 ) THEN
+           IF ( kref .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.
@@ -403,5 +454,5 @@
               IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=7.
            END IF
-           IF ( k .EQ. 25 ) THEN
+           IF ( kref .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.
@@ -411,5 +462,5 @@
               IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=8.5
            END IF
-           IF ( k .EQ. 26 ) THEN
+           IF ( kref .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.
@@ -419,7 +470,9 @@
               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
+!!jyg
+!!              IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=6.
+!!jyg end
+           END IF
+           IF ( kref .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.
@@ -430,5 +483,5 @@
               IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=10.
            END IF
-           IF ( k .EQ. 28 ) THEN
+           IF ( kref .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.
@@ -439,5 +492,5 @@
               IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=10.
            END IF
-           IF ( k .EQ. 29 ) THEN
+           IF ( kref .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.
@@ -449,5 +502,5 @@
               IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=10.
            END IF
-           IF ( k .EQ. 30 ) THEN
+           IF ( kref .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.
@@ -459,5 +512,5 @@
               IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=9.
            END IF
-           IF ( k .EQ. 31 ) THEN
+           IF ( kref .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.
@@ -469,5 +522,5 @@
               IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=7.6
            END IF
-           IF ( k .EQ. 32 ) THEN
+           IF ( kref .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.
@@ -479,5 +532,5 @@
               IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=7.4
            END IF
-           IF ( k .GE. 33 ) THEN
+           IF ( kref .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.
@@ -491,5 +544,5 @@
         END DO
      END DO
-  END IF ! fin de 39 niveaux verticaux
+!!!!!!  END IF ! fin de 39 niveaux verticaux
 
 
@@ -503,4 +556,6 @@
        ! 1/(60*1.295) = 0.01287
        srcbe(i,k)=srcbe(i,k)*0.01287
+!!       print *,'  k, srcbe(i,k) ',   &
+!!                  k, srcbe(i,k)
        ! La source est  at/min/m3 -> at/s/m3
        ! srcbe(i,k)=srcbe(i,k)*0.0166667
