Index: trunk/LMDZ.VENUS/libf/phyvenus/clmain.F
===================================================================
--- trunk/LMDZ.VENUS/libf/phyvenus/clmain.F	(revision 1572)
+++ trunk/LMDZ.VENUS/libf/phyvenus/clmain.F	(revision 1591)
@@ -291,6 +291,6 @@
 c        call coefkzmin(knon,ypaprs,ypplay,yu,yv,yt,ycoefm
 c    .   ,ycoefm0,ycoefh0)
-c      call dump2d(iim,jjm-1,ycoefm(2:klon-1,2), 'KZ         ')
-c      call dump2d(iim,jjm-1,ycoefm0(2:klon-1,2),'KZMIN      ')
+c      call dump2d(nbp_lon,nbp_lat-2,ycoefm(2:klon-1,2), 'KZ         ')
+c      call dump2d(nbp_lon,nbp_lat-2,ycoefm0(2:klon-1,2),'KZMIN      ')
  
          ycoefm0 = 1.e-3
@@ -349,5 +349,5 @@
         ENDIF
 
-c   iflag_pbl peut etre utilise comme longuer de melange
+c   iflag_pbl peut etre utilise comme longueur de melange
 
          if (iflag_pbl.ge.11) then
@@ -558,4 +558,28 @@
       PARAMETER (check=.false.)
 C
+c----------------------
+c ATMOSPHERE PROFONDE
+      real p_lim,p_ref
+      real rsmu_mod(klon,klev),zrsmu_mod(klon,klev)
+      real ratio_mod(klon,klev) 
+! theta_mod = theta * ratio_mod => ratio_mod = mmm0/mmm(p)
+
+      p_lim = 6e6
+      p_ref = 8.9e6
+
+      DO k = 1, klev
+      DO i = 1, klon
+        ratio_mod(i,k) = 1.
+        rsmu_mod(i,k) = RD
+       if (pplay(i,k).gt.p_lim) then
+          ratio_mod(i,k) = RMD /
+     &  (RMD+0.56*(log(pplay(i,k))-log(p_lim))/(log(p_ref)-log(p_lim)))
+           ratio_mod(i,k) = max(ratio_mod(i,k),RMD/(RMD+0.56))
+           rsmu_mod(i,k) = RD*ratio_mod(i,k)
+       endif
+      ENDDO
+      ENDDO
+c----------------------
+
       if (iflag_pbl.eq.1) then
         do k = 3, klev
@@ -584,4 +608,8 @@
       DO i = 1, knon
             zt(i,k)    = (t(i,k)+t(i,k-1)) * 0.5
+c----------------------
+c ATMOSPHERE PROFONDE
+            zrsmu_mod(i,k)    = (rsmu_mod(i,k)+rsmu_mod(i,k-1)) * 0.5
+c----------------------
       ENDDO
       ENDDO
@@ -612,10 +640,13 @@
       DO i = 1, knon
          zx_coef(i,k) = coef(i,k)*RG/(pplay(i,k-1)-pplay(i,k))
-     .                  *(paprs(i,k)/zt(i,k)/RD)**2
+c----------------------
+c ATMOSPHERE PROFONDE
+     .                  *(paprs(i,k)/zt(i,k)/zrsmu_mod(i,k))**2
+c----------------------
          zx_coef(i,k) = zx_coef(i,k) * dtime*RG
       ENDDO
       ENDDO
 c
-c Preparer les flux lies aux contre-gardients
+c Preparer les flux lies aux contre-gardients  OBSOLETE
 c
       DO k = 2, klev
@@ -984,4 +1015,22 @@
       DATA appel1er /.TRUE./
 c
+c----------------------
+c ATMOSPHERE PROFONDE
+      real p_lim,p_ref
+      real modgam(klon,klev)
+
+      p_lim = 6e6
+      p_ref = 8.9e6
+
+      DO k = 1, klev
+      DO i = 1, klon
+        modgam(i,k) = 0.
+       if (pplay(i,k).gt.p_lim) then
+           modgam(i,k) = 0.56/(log(p_ref)-log(p_lim))/R
+       endif
+      ENDDO
+      ENDDO
+c----------------------
+
       isommet=klev
 
@@ -1045,6 +1094,11 @@
             zmgeom(i,k)= zgeop(i,k)-zgeop(i,k-1)
             zdphi      = zmgeom(i,k)/2.
-            ztvd(i,k)  = (t(i,k)   + zdphi/cpdet(zt(i,k)))
-            ztvu(i,k)  = (t(i,k-1) - zdphi/cpdet(zt(i,k)))
+c----------------------
+c ATMOSPHERE PROFONDE
+            ztvd(i,k)  = t(i,k)   + zdphi*
+     &                      (1./cpdet(zt(i,k))+modgam(i,k))
+            ztvu(i,k)  = t(i,k-1) - zdphi*
+     &                      (1./cpdet(zt(i,k))+modgam(i,k))
+c----------------------
             zpk(i,k)   = ppk(i,k)*(paprs(i,k)/pplay(i,k))**RKAPPA
       ENDDO
@@ -1054,5 +1108,9 @@
         zt(i,1)   = ts(i)
         ztvu(i,1) = ts(i) 
-        ztvd(i,1) = t(i,1)+zgeop(i,1)/cpdet(zt(i,1)) 
+c----------------------
+c ATMOSPHERE PROFONDE
+        ztvd(i,1) = t(i,1)+zgeop(i,1)*
+     &                      (1./cpdet(zt(i,1))+modgam(i,1))
+c----------------------
         zpk(i,1)  = ppk(i,1)*(paprs(i,1)/pplay(i,1))**RKAPPA
       ENDDO
Index: trunk/LMDZ.VENUS/libf/phyvenus/concentrations2.F
===================================================================
--- trunk/LMDZ.VENUS/libf/phyvenus/concentrations2.F	(revision 1572)
+++ trunk/LMDZ.VENUS/libf/phyvenus/concentrations2.F	(revision 1591)
@@ -174,10 +174,4 @@
 c         endif
 
-  !	n2 
-c            akin2 = 5.6e-4
-c            cpin2 = 1.034e3
-
-
-         
          ! tell the world about it:
          write(*,*) "concentrations: firstcall, nbq=",nbq
@@ -224,7 +218,6 @@
             end do
             mmean(ig,l) = 1./mmean(ig,l)
-            rnew(ig,l) = 8.314/mmean(ig,l)*1.e3     ! J/kg K		
- 
-c            write(*,*),'Mmean: ',ig, l, mmean(ig,l) 
+            rnew(ig,l) = 8.314/mmean(ig,l)*1.e3     ! J/kg K 
+c            write(*,*),'Mmean in concentration2: ',ig, l, mmean(ig,l) 
          end do
       end do
@@ -244,5 +237,4 @@
 c            write(*,*),'Air density: ',ig, l, rho(0,l)            
 
-!!! --- INSERT N2 values ---- 
 !!  WARNING -> Cp here below doesn't depend on T (cpdet)
 
Index: trunk/LMDZ.VENUS/libf/phyvenus/diffusion.h
===================================================================
--- trunk/LMDZ.VENUS/libf/phyvenus/diffusion.h	(revision 1572)
+++ trunk/LMDZ.VENUS/libf/phyvenus/diffusion.h	(revision 1591)
@@ -9,6 +9,6 @@
       real*8 dzres
 
-      parameter (Pdiff=15.)      ! pressure below which diffusion is computed
-      parameter (tdiffmin=5d0)
-      parameter (dzres=2d0)	 ! grid resolution (km) for diffusion
+      parameter (Pdiff=1.e-1)      ! pressure below which diffusion is computed
+      parameter (tdiffmin=5.d0)
+      parameter (dzres=2.d0)	 ! grid resolution (km) for diffusion
       
Index: trunk/LMDZ.VENUS/libf/phyvenus/flott_gwd_ran.F90
===================================================================
--- trunk/LMDZ.VENUS/libf/phyvenus/flott_gwd_ran.F90	(revision 1572)
+++ trunk/LMDZ.VENUS/libf/phyvenus/flott_gwd_ran.F90	(revision 1591)
@@ -5,11 +5,14 @@
     ! Parametrization of the momentum flux deposition due to a discrete
     ! number of gravity waves. 
-    ! F. Lott (version 9: 16 February, 2012), reproduce v3 but with only
-    !         two waves present at each time step
+    ! F. Lott
+    ! Version 14, Gaussian distribution of the source
     ! LMDz model online version      
-    ! ADAPTED FOR VENUS
+    ! ADAPTED FOR VENUS /  F. LOTT + S. LEBONNOIS
+    ! Version adapted on 03/04/2013:
+    !      - input flux compensated in the deepest layers
     !---------------------------------------------------------------------
 
       use dimphy
+      use moyzon_mod, only: plevmoy
       implicit none
 
@@ -39,5 +42,5 @@
     ! O.3 INTERNAL ARRAYS
 
-    INTEGER II, LL, IEQ
+    INTEGER II, LL
 
     ! 0.3.0 TIME SCALE OF THE LIFE CYCLE OF THE WAVES PARAMETERIZED
@@ -62,21 +65,19 @@
     REAL ZOM(NW, KLON), ZOP(NW, KLON)
 
-    ! Wave vertical velocities at the 2 1/2 lev surrounding the full level
+    ! Wave EP-fluxes at the 2 semi levels surrounding the full level
     REAL WWM(NW, KLON), WWP(NW, KLON)
-
-    REAL RUW0(NW, KLON) ! Fluxes at launching level
-
+    ! Fluxes X and Y for each waves at 1/2 Levels
     REAL RUWP(NW, KLON), RVWP(NW, KLON)
-    ! Fluxes X and Y for each waves at 1/2 Levels
-
-    INTEGER LAUNCH   ! Launching altitude
-
-    REAL RUWMAX,SAT  ! saturation parameter
-    REAL XLAUNCH     ! Controle the launching altitude
     REAL RUW(KLON, KLEV + 1) ! Flux x at semi levels
     REAL RVW(KLON, KLEV + 1) ! Flux y at semi levels
 
+    REAL RUW0(NW, KLON) ! Fluxes at launching level
+    REAL RUWMAX         ! Max value
+    INTEGER LAUNCH      ! Launching altitude
+    REAL XLAUNCH        ! Control the launching altitude
+
     ! 0.3.2 PARAMETERS OF WAVES DISSIPATIONS
 
+    REAL SAT  ! saturation parameter
     REAL RDISS, ZOISEC ! COEFF DE DISSIPATION, SECURITY FOR INTRINSIC FREQ
 
@@ -99,7 +100,9 @@
     logical output
     data output/.false./
-  ! CAUTION ! IF output is .true. THEN change NO to 10 at least !
+!    data output/.true./
+ ! CAUTION ! IF output is .true. THEN change NO to 10 at least !
     character*14 outform
     character*2  str2
+    integer      ieq
 
 ! ON CONSERVE LA MEMOIRE un certain temps AVEC UN SAVE
@@ -127,21 +130,25 @@
     RUWMAX = 0.005      ! Max EP-Flux at Launch altitude
     SAT    = 0.85       ! Saturation parameter: Sc in (12)
-    RDISS  = 10.        ! Diffusion parameter 
+    RDISS  = 0.1        ! Diffusion parameter 
 
     DELTAT=24.*3600.    ! Time scale of the waves (first introduced in 9b)
 
-    KMIN = 1.E-6        ! Min horizontal wavenumber
-    KMAX = 2.E-5        ! Max horizontal wavenumber
+!!!! TEST GG. Values corresponding to min/max horizontal wavel 50-500 km
+!(similar to observations)
+    KMIN = 1.E-5        ! Min horizontal wavenumber
+    KMAX = 1.E-4        ! Max horizontal wavenumber
     !Online output: one value only
     if (output) then
-      KMIN = 6.3E-6
-      KMAX = 6.3E-6
+      KMIN = 3E-5
+      KMAX = 3E-5
     endif
     CMIN = 1.           ! Min phase velocity
     CMAX = 61.          ! Max phase speed velocity
-    XLAUNCH=0.6         ! Parameter that control launching altitude
-
-    PR = 9.2e6          ! Reference pressure    ! VENUS!!
-    TR = 240.           ! Reference Temperature ! VENUS: cloud layer
+!    XLAUNCH=0.6         ! Parameter that control launching altitude
+    XLAUNCH=5e-3        ! Value for top of cloud convective region
+
+!    PR = 9.2e6          ! Reference pressure    ! VENUS!!
+    PR = 5e5            ! Reference pressure    ! VENUS: cloud layer
+    TR = 300.           ! Reference Temperature ! VENUS: cloud layer
     H0 = RD * TR / RG   ! Characteristic vertical scale height
 
@@ -162,5 +169,59 @@
     ENDIF
 
-    ! 1.2 WAVES CHARACTERISTICS CHOSEN RANDOMLY
+    !  2. EVALUATION OF THE BACKGROUND FLOW AT SEMI-LEVELS
+    !-------------------------------------------------------------
+
+    !Online output
+    if (output) OPEN(11,file="impact-gwno.dat")
+
+    ! Pressure and Inv of pressure, Temperature / at 1/2 level
+    DO LL = 2, KLEV
+       PH(:, LL) = EXP((LOG(PP(:, LL)) + LOG(PP(:, LL - 1))) / 2.)
+    end DO
+
+    PH(:, KLEV + 1) = 0. 
+    PH(:, 1) = 2. * PP(:, 1) - PH(:, 2)
+
+    ! Launching altitude
+
+    DO LL = 1, NLEV
+       IF (plevmoy(LL) / plevmoy(1) > XLAUNCH) LAUNCH = LL
+    ENDDO
+! test
+!    print*,"launch=",LAUNCH
+!    print*,"launch p,N2=",plevmoy(LAUNCH),pn2(nlon/2+1,LAUNCH)
+
+    ! Log pressure vert. coordinate
+    DO LL = 1, KLEV + 1 
+       ZH(:, LL) = H0 * LOG(PR / (PH(:, LL) + PSEC))
+    end DO
+
+    if (output) then
+    ! altitude above surface
+       ZHbis(:,1) = 0.    
+       DO LL = 2, KLEV + 1 
+          H0bis(:, LL-1) = RD * TT(:, LL-1) / RG 
+          ZHbis(:, LL) = ZHbis(:, LL-1) &
+           + H0bis(:, LL-1)*(PH(:, LL-1)-PH(:,LL))/PP(:, LL-1)
+       end DO
+    endif
+
+    ! Winds and BV frequency
+    DO LL = 2, KLEV
+       UH(:, LL) = 0.5 * (UU(:, LL) + UU(:, LL - 1)) ! Zonal wind
+       VH(:, LL) = 0.5 * (VV(:, LL) + VV(:, LL - 1)) ! Meridional wind
+       ! BVSEC: BV Frequency
+! VENUS ATTENTION: CP VARIABLE PSTAB CALCULE EN AMONT DES PARAMETRISATIONS
+       BV(:, LL) = MAX(BVSEC,SQRT(pn2(:,LL)))
+    end DO
+    BV(:, 1) = BV(:, 2)
+    UH(:, 1) = 0.
+    VH(:, 1) = 0.
+    BV(:, KLEV + 1) = BV(:, KLEV)
+    UH(:, KLEV + 1) = UU(:, KLEV)
+    VH(:, KLEV + 1) = VV(:, KLEV)
+
+
+    ! 3. WAVES CHARACTERISTICS CHOSEN RANDOMLY
     !-------------------------------------------
 
@@ -168,4 +229,8 @@
     ! are used to produce the waves characteristics
     ! in a stochastic way
+
+!! A REVOIR: 
+!! - utilisation de MOD ou bien de aleas ?
+!! - distribution gaussienne des CPHA ? (avec signe ZP qui est ajuste apres)
 
     JW = 0
@@ -187,4 +252,8 @@
                 ! Intrinsic frequency
                 ZO(JW, II) = CPHA * ZK(JW, II)
+                ! Intrinsic frequency  is imposed
+                    ZO(JW, II) = ZO(JW, II)      &
+                  + ZK(JW, II) * COS(ZP(JW)) * UH(II, LAUNCH) &
+                  + ZK(JW, II) * SIN(ZP(JW)) * VH(II, LAUNCH)
                 ! Momentum flux at launch lev 
                 ! RUW0(JW, II) = RUWMAX / REAL(NW) &
@@ -199,58 +268,8 @@
     end DO
 
-    !  2. EVALUATION OF THE BACKGROUND FLOW AT SEMI-LEVELS
-    !-------------------------------------------------------------
-
-    IEQ = KLON / 2
-    !Online output
-    if (output) OPEN(11,file="impact-gwno.dat")
-
-    ! Pressure and Inv of pressure, Temperature / at 1/2 level
-    DO LL = 2, KLEV
-       PH(:, LL) = EXP((LOG(PP(:, LL)) + LOG(PP(:, LL - 1))) / 2.)
-    end DO
-
-    PH(:, KLEV + 1) = 0. 
-    PH(:, 1) = 2. * PP(:, 1) - PH(:, 2)
-
-    ! Launching altitude
-
-    DO LL = 1, NLEV
-       IF (PH(IEQ, LL) / PH(IEQ, 1) > XLAUNCH) LAUNCH = LL
-    ENDDO
-
-    ! Log pressure vert. coordinate (altitude above surface) 
-    ZHbis(:,1) = 0.    
-    DO LL = 2, KLEV + 1 
-       H0bis(:, LL-1) = RD * TT(:, LL-1) / RG 
-       ZHbis(:, LL) = ZHbis(:, LL-1) &
-           + H0bis(:, LL-1)*(PH(:, LL-1)-PH(:,LL))/PP(:, LL-1)
-    end DO
-    ! Log pressure vert. coordinate
-    DO LL = 1, KLEV + 1 
-       ZH(:, LL) = H0 * LOG(PR / (PH(:, LL) + PSEC))
-    end DO
-
-
-    ! Winds and BV frequency
-    DO LL = 2, KLEV
-       UH(:, LL) = 0.5 * (UU(:, LL) + UU(:, LL - 1)) ! Zonal wind
-       VH(:, LL) = 0.5 * (VV(:, LL) + VV(:, LL - 1)) ! Meridional wind
-       ! BVSEC: BV Frequency
-! VENUS ATTENTION: CP VARIABLE PSTAB CALCULE EN AMONT DES PARAMETRISATIONS
-       BV(:, LL) = MAX(BVSEC,SQRT(pn2(:,LL)))
-    end DO
-    BV(:, 1) = BV(:, 2)
-    UH(:, 1) = 0.
-    VH(:, 1) = 0.
-    BV(:, KLEV + 1) = BV(:, KLEV)
-    UH(:, KLEV + 1) = UU(:, KLEV)
-    VH(:, KLEV + 1) = VV(:, KLEV)
-
-
-    ! 3. COMPUTE THE FLUXES
+    ! 4. COMPUTE THE FLUXES
     !--------------------------
 
-    !  3.1  Vertical velocity at launching altitude to ensure 
+    !  4.1  Vertical velocity at launching altitude to ensure 
     !       the correct value to the imposed fluxes.
     !
@@ -261,8 +280,6 @@
             - ZK(JW, :) * COS(ZP(JW)) * UH(:, LAUNCH) &
             - ZK(JW, :) * SIN(ZP(JW)) * VH(:, LAUNCH) 
-       ! Vertical velocity at launch level, value to ensure the imposed
-       ! mom flux:
-       WWP(JW, :) = SQRT(ABS(ZOP(JW, :)) / MAX(BV(:, LAUNCH),BVSEC) &
-            * RUW0(JW,:))
+       ! WW is directly a flux, here, not vertical velocity anymore
+       WWP(JW, :) = RUW0(JW,:)
        RUWP(JW, :) = COS(ZP(JW)) * SIGN(1., ZOP(JW, :)) * RUW0(JW, :)
        RVWP(JW, :) = SIN(ZP(JW)) * SIGN(1., ZOP(JW, :)) * RUW0(JW, :)
@@ -270,16 +287,14 @@
     end DO
 
-    !  3.2 Uniform values below the launching altitude
-
-    DO LL = 1, LAUNCH
-       RUW(:, LL) = 0
-       RVW(:, LL) = 0
-       DO JW = 1, NW
-          RUW(:, LL) = RUW(:, LL) + RUWP(JW, :)
-          RVW(:, LL) = RVW(:, LL) + RVWP(JW, :)
-       end DO
-    end DO
-
-    !  3.3 Loop over altitudes, with passage from one level to the
+    !  4.2 Initial flux at launching altitude
+
+    RUW(:, LAUNCH) = 0
+    RVW(:, LAUNCH) = 0
+    DO JW = 1, NW
+       RUW(:, LAUNCH) = RUW(:, LAUNCH) + RUWP(JW, :)
+       RVW(:, LAUNCH) = RVW(:, LAUNCH) + RVWP(JW, :)
+    end DO
+
+    !  4.3 Loop over altitudes, with passage from one level to the
     !      next done by i) conserving the EP flux, ii) dissipating
     !      a little, iii) testing critical levels, and vi) testing
@@ -287,8 +302,11 @@
 
     !Online output
-    write(str2,'(i2)') NW+2
-    outform="("//str2//"(E12.4,1X))"
-    if (output) WRITE(11,outform) ZH(IEQ, 1) / 1000., ZHbis(IEQ, 1) / 1000., &
+    if (output) then
+        ieq=nlon/2+1
+        write(str2,'(i2)') NW+2
+        outform="("//str2//"(E12.4,1X))"
+        WRITE(11,outform) ZH(IEQ, 1) / 1000., ZHbis(IEQ, 1) / 1000., &
                (ZO(JW, IEQ)/ZK(JW, IEQ)*COS(ZP(JW)), JW = 1, NW)
+    endif
 
     DO LL = LAUNCH, KLEV - 1
@@ -305,19 +323,16 @@
 
           WWP(JW, :) = MIN( & 
-               ! No breaking (Eq.6)
+       ! No breaking (Eq.6)
                WWM(JW, :) & 
-               * SQRT(BV(:, LL ) / BV(:, LL+1) &
-               * ABS(ZOP(JW, :)) / MAX(ABS(ZOM(JW, :)), ZOISEC)) &
-               ! Dissipation (Eq. 8):
+       ! Dissipation (Eq. 8):
                * EXP(- RDISS * PR / (PH(:, LL + 1) + PH(:, LL)) &
                * ((BV(:, LL + 1) + BV(:, LL)) / 2.)**3 &
                / MAX(ABS(ZOP(JW, :) + ZOM(JW, :)) / 2., ZOISEC)**4 &
                * ZK(JW, :)**3 * (ZH(:, LL + 1) - ZH(:, LL))), &
-               ! Critical levels (forced to zero if intrinsic
-               ! frequency changes sign)
+       ! Critical levels (forced to zero if intrinsic frequency changes sign)
                MAX(0., SIGN(1., ZOP(JW, :) * ZOM(JW, :))) &
-               ! Saturation (Eq. 12)
-               * ZOP(JW, :)**2 / ZK(JW, :)/BV(:, LL+1) & 
-               * EXP(-ZH(:, LL + 1)/2./H0) * SAT)
+       ! Saturation (Eq. 12)
+               * ABS(ZOP(JW, :))**3 /BV(:, LL+1) & 
+               * EXP(-ZH(:, LL + 1)/H0) * SAT**2*KMIN**2/ZK(JW, :)**4)  
        end DO
 
@@ -327,10 +342,6 @@
 
        DO JW = 1, NW
-          RUWP(JW, :) = ZOP(JW, :)/MAX(ABS(ZOP(JW, :)), ZOISEC)**2 &
-               *BV(:,LL+1)&
-               * COS(ZP(JW)) *  WWP(JW, :)**2
-          RVWP(JW, :) = ZOP(JW, :)/MAX(ABS(ZOP(JW, :)), ZOISEC)**2 &
-               *BV(:,LL+1)&
-               * SIN(ZP(JW)) *  WWP(JW, :)**2
+          RUWP(JW, :) = SIGN(1.,ZOP(JW, :))*COS(ZP(JW))*WWP(JW, :)
+          RVWP(JW, :) = SIGN(1.,ZOP(JW, :))*SIN(ZP(JW))*WWP(JW, :)
        end DO
        !
@@ -358,25 +369,40 @@
     end DO
 
+    ! 5 CALCUL DES TENDANCES:
+    !------------------------
+
+    ! 5.1 Rectification des flux au sommet et dans les basses couches:
+! MODIF SL
+ 
+! Attention, ici c'est le total sur toutes les ondes...
+
+    RUW(:, KLEV + 1) = 0.
+    RVW(:, KLEV + 1) = 0.
+
+    ! Here, big change compared to FLott version:
+    ! We compensate (RUW(:, LAUNCH), ie total emitted upward flux
+    !  over the layers max(1,LAUNCH-3) to LAUNCH-1
+    DO LL = 1, max(1,LAUNCH-3)
+      RUW(:, LL) = 0.
+      RVW(:, LL) = 0.
+    end DO
+    DO LL = max(2,LAUNCH-2), LAUNCH-1
+       RUW(:, LL) = RUW(:, LL - 1) + RUW(:, LAUNCH) * &
+            (PH(:,LL)-PH(:,LL-1)) / (PH(:,LAUNCH)-PH(:,max(1,LAUNCH-3)))
+       RVW(:, LL) = RVW(:, LL - 1) + RVW(:, LAUNCH) * &
+            (PH(:,LL)-PH(:,LL-1)) / (PH(:,LAUNCH)-PH(:,max(1,LAUNCH-3)))
+    end DO
+    ! This way, the total flux from GW is zero, but there is a net transport
+    ! (upward) that should be compensated by circulation 
+    ! and induce additional friction at the surface 
+
     !Online output
     if (output) then 
+       DO LL = 1, KLEV - 1
+           WRITE(11,*) ZHbis(IEQ, LL)/1000.,RUW(IEQ,LL)
+       end DO
        CLOSE(11)
        stop
     endif
-
-    ! 4 CALCUL DES TENDANCES:
-    !------------------------
-
-    ! 4.1 Rectification des flux au sommet et dans les basses couches:
-
-    RUW(:, KLEV + 1) = 0.
-    RVW(:, KLEV + 1) = 0.
-    RUW(:, 1) = RUW(:, LAUNCH)
-    RVW(:, 1) = RVW(:, LAUNCH)
-    DO LL = 2, LAUNCH
-       RUW(:, LL) = RUW(:, LL - 1) + (RUW(:, LAUNCH + 1) - RUW(:, 1)) * &
-            (PH(:, LL) - PH(:, LL - 1)) / (PH(:, LAUNCH + 1) - PH(:, 1))
-       RVW(:, LL) = RVW(:, LL - 1) + (RVW(:, LAUNCH + 1) - RVW(:, 1)) * &
-            (PH(:, LL) - PH(:, LL - 1)) / (PH(:, LAUNCH + 1) - PH(:, 1))
-    end DO
 
     ! AR-1 RECURSIVE FORMULA (13) IN VERSION 4
Index: trunk/LMDZ.VENUS/libf/phyvenus/moldiff_red.F90
===================================================================
--- trunk/LMDZ.VENUS/libf/phyvenus/moldiff_red.F90	(revision 1572)
+++ trunk/LMDZ.VENUS/libf/phyvenus/moldiff_red.F90	(revision 1591)
@@ -226,5 +226,5 @@
 !!!
 !! =======> Gab 29 oct 2014
-!!$$$$$$$   THIS UPDATE IS ALREADY DONE in physique (t_seri & tr_seri) $$$$$$
+!!$$$   THIS UPDATE IS ALREADY DONE in physique (t_seri & tr_seri) $$$$$$
 
 ! Update the temperature modified by other processes
@@ -449,10 +449,16 @@
 !        Tdiff=Tdiff/2.
 
-!	if (ig .eq. ij0) then 
-!	print*,'test',ig,Tdiff,tdiffmin,minloc(Tdiffrafmol),minloc(Tdiffrafmol(nlraf,:))
-!	endif
-	if (tdiff .lt. tdiffmin*Mraf(nlraf)) tdiff=tdiffmin*Mraf(nlraf)
-	ntime=anint(dble(ptimestep)/Tdiff)
-!	print*,'ptime',ig,ptimestep,Tdiff,ntime,tdiffmin,Mraf(nlraf)
+        if (tdiff .lt. tdiffmin*Mraf(nlraf)) tdiff=tdiffmin*Mraf(nlraf)
+
+!     if (tdiff .lt. tdiffmin*Mraf(nlraf)) then       
+!        print*, 'Moldiff L454', tdiff, ptimestep, tdiffmin*Mraf(nlraf)
+!        endif
+
+!      JYC + GG aout 2015 add this condition below  
+        if (tdiff .ge. ptimestep) tdiff=ptimestep
+        ! Number of time step
+      ntime=anint(dble(ptimestep)/tdiff)
+      !print*,'ptime',ig,tdiff,ntime,Mraf(nlraf)
+
 ! Adimensionned temperature
 
@@ -592,4 +598,7 @@
 	CALL GCMGRID_P2(Zraf,Praf,Qraf,Traf,Nrafk,Rrafk,qq,qnew,tt,tnew,&
      &   pp,mol_tr,gcmind,nlraf,ncompdiff,nlayer,FacMass,ig)
+
+!!!  Call added by JY+GG Aout 2014
+        CALL MMOY(massemoy,mol_tr,qnew,gcmind,nlayer,ncompdiff)
 
 	CALL RHOTOT(pp,tt,massemoy,qnew,RHOT,RHOK,nlayer,ncompdiff)
@@ -994,12 +1003,13 @@
 	SUBROUTINE ZVERT(P,T,M,Z,nl,ig)
 	IMPLICIT NONE
+#include "YOMCST.h"
 	INTEGER :: nl,l,ig
 	REAL*8,dimension(nl) :: P,T,M,Z,H
 	REAL*8 :: z0
 	REAL*8 :: kbolt,masseU,Konst,g,Hpm
-	masseU=1.660538782d-27
-        kbolt=1.3806504d-23
+	masseU=1.e-3/RNAVO
+        kbolt=RKBOL
 	Konst=kbolt/masseU
-	g=3.72D0
+	g=RG
 
 	z0=0d0
@@ -1025,11 +1035,11 @@
 	SUBROUTINE RHOTOT(P,T,M,qq,rhoN,rhoK,nl,nq)
 	IMPLICIT NONE
-
+#include "YOMCST.h"
 	REAL*8 :: kbolt,masseU,Konst
 	INTEGER :: nl,nq,l,iq
 	REAL*8,DIMENSION(nl) :: P,T,M,RHON
 	REAL*8,DIMENSION(nl,nq) :: RHOK,qq
-	masseU=1.660538782d-27
-        kbolt=1.3806504d-23
+	masseU=1.e-3/RNAVO
+        kbolt=RKBOL
 	Konst=Kbolt/masseU
 
@@ -1049,5 +1059,5 @@
 	use infotrac
 	IMPLICIT NONE
-	
+#include "YOMCST.h"
 	INTEGER :: nl,nq,il,l,i,iq,nlx,iz,ig
 	INTEGER :: gc(nq)
@@ -1060,8 +1070,8 @@
 	REAL*8 :: facZ,dZ,H
 	REAL,DIMENSION(nq) :: mol_tr
-	masseU=1.660538782d-27
-        kbolt=1.3806504d-23
+	masseU=1.e-3/RNAVO
+        kbolt=RKBOL
         Konst=Kbolt/masseU
-	g=3.72d0
+	g=RG
 
 
@@ -1206,12 +1216,13 @@
 	SUBROUTINE VELVERT(nn,T,H,D,Dz,masse,W,nl)
 	IMPLICIT NONE
+#include "YOMCST.h"
 	INTEGER :: l,nl,nn
 	REAL*8,DIMENSION(nl) :: T,H,D,W,DT
 	REAL*8 :: Dz,Hmol,masse
 	REAL*8 :: kbolt,masseU,Konst,g
-	masseU=1.660538782d-27
-        kbolt=1.3806504d-23
+	masseU=1.e-3/RNAVO
+        kbolt=RKBOL
         Konst=Kbolt/masseU
-        g=3.72d0
+	g=RG
 
 	DT(1)=1D0/T(1)*(-3D0*T(1)+4D0*T(2)-T(3))/(2D0*DZ)
@@ -1366,9 +1377,10 @@
 	SUBROUTINE Checkmass2(qold,qnew,P,il,nl,nn,nq)
 	IMPLICIT NONE
+#include "YOMCST.h"
 	INTEGER :: nl,nn,l,nq,il
 	REAL,DIMENSION(nl+1) :: P
 	REAL*8,DIMENSION(nl,nq) :: qold,qnew
 	REAL*8 :: DM,Mold,Mnew,g
-	g=3.72d0
+	g=RG
 	DM=0d0
 	Mold=0d0
@@ -1381,5 +1393,5 @@
 	ENDDO
 	IF (abs(DM/Mold) .gt. 1d-5) THEN
-	Print*,'We dont conserve mas',nn,DM,Mold,Mnew,DM/Mold
+	Print*,'We dont conserve mass',nn,DM,Mold,Mnew,DM/Mold
 	ENDIF
 
@@ -1391,4 +1403,5 @@
 	use infotrac
         IMPLICIT NONE
+#include "YOMCST.h"
 	INTEGER :: nl,nq,nlx,il,nn,iP,ig,compteur
 	INTEGER,DIMENSION(1) :: indP
@@ -1403,8 +1416,8 @@
 	REAL*8 :: kbolt,masseU,Konst,g,Dz,facP,Hi
         REAL*8 :: Znew,Znew2,Pnew,Pnew2
-        masseU=1.660538782d-27
-        kbolt=1.3806504d-23
+	masseU=1.e-3/RNAVO
+        kbolt=RKBOL
         Konst=Kbolt/masseU
-        g=3.72d0
+	g=RG
 	Dz=Z(2)-Z(1)
 	Znew=Z(nl)
@@ -1493,4 +1506,5 @@
 !	 use infotrac
         IMPLICIT NONE
+#include "YOMCST.h"
         INTEGER :: nl,nq,nlx,il,nn,iP,ig,compteur
         INTEGER,DIMENSION(1) :: indP
@@ -1504,8 +1518,8 @@
         REAL*8 :: kbolt,masseU,Konst,g,Dz,facP,Hi
         REAL*8 :: Znew,Znew2,Pnew,Pnew2
-        masseU=1.660538782d-27
-        kbolt=1.3806504d-23
+	masseU=1.e-3/RNAVO
+        kbolt=RKBOL
         Konst=Kbolt/masseU
-        g=3.72d0
+	g=RG
         Dz=Z(2)-Z(1)
         Znew=Z(nl)
Index: trunk/LMDZ.VENUS/libf/phyvenus/nlte_tcool.F
===================================================================
--- trunk/LMDZ.VENUS/libf/phyvenus/nlte_tcool.F	(revision 1572)
+++ trunk/LMDZ.VENUS/libf/phyvenus/nlte_tcool.F	(revision 1591)
@@ -82,7 +82,8 @@
             cpnew_ig(l)=cpnew(ig,l)
 
+c            write (*,*) 'nlte_tcool z_ig', l, z_ig(l)    
          enddo 
 
-                                ! From GCM's grid to NLTE's grid
+c     ! From GCM's grid to NLTE's grid
          call NLTEdlvr11_ZGRID_02 (nlev,
      $        p_ig, t_ig, z_ig,
@@ -302,5 +303,7 @@
       zmin = z_gcm(jlowerboundary)
       zmax = z_gcm(jtopboundary)
-!      print*, zmin, zmax
+
+!      print*, 'nlte boundaries', zmin, zmax
+!      print*, 'Top atmosphere', z_gcm(nlev-1)
 !      stop
       deltaz = (zmax-zmin) / (nl-1)
Index: trunk/LMDZ.VENUS/libf/phyvenus/nlthermeq.F
===================================================================
--- trunk/LMDZ.VENUS/libf/phyvenus/nlthermeq.F	(revision 1572)
+++ trunk/LMDZ.VENUS/libf/phyvenus/nlthermeq.F	(revision 1591)
@@ -64,6 +64,4 @@
       END IF
 
-c    write(*,*) 'LTE rad. calculations up to layer ',  nlaylte
-
 c
  100    return
@@ -73,5 +71,5 @@
    20 format('               half-width (scale heights) ',f6.2)
    30 format('          suggested LTE coverage at least ',f6.2,'Pa')
-   40 format(' nlthermeq: purely NLTE contribution over (nlayer) ',f6.2)  
+   40 format(' nlthermeq: purely NLTE contribution over (nlayer) ',f6.4)
 
       end
Index: trunk/LMDZ.VENUS/libf/phyvenus/radlwsw.F
===================================================================
--- trunk/LMDZ.VENUS/libf/phyvenus/radlwsw.F	(revision 1572)
+++ trunk/LMDZ.VENUS/libf/phyvenus/radlwsw.F	(revision 1591)
@@ -54,5 +54,5 @@
       REAL   PPB(klev+1)
 
-      REAL   zfract, zrmu0
+      REAL   zfract, zrmu0,latdeg
 
       REAL   zheat(klev), zcool(klev)
@@ -273,4 +273,6 @@
 c---------
       znivs=zzlev(j,:)
+      latdeg=abs(latitude_deg(j))
+
 c      CALL SW_venus_ve_1Dglobave(zrmu0, zfract,   ! pour moy globale
 c      CALL SW_venus_ve(zrmu0, zfract,
@@ -279,8 +281,10 @@
 c     S        ztopsw,zsolsw,ZFSNET)
 
-c      CALL SW_venus_cl_1Dglobave(zrmu0, zfract,   ! pour moy globale
-c      CALL SW_venus_cl(zrmu0, zfract,
-c      CALL SW_venus_dc_1Dglobave(zrmu0, zfract,   ! pour moy globale
-      CALL SW_venus_dc(zrmu0, zfract,
+c      CALL SW_venus_cl_1Dglobave(zrmu0,zfract,   ! pour moy globale
+c      CALL SW_venus_cl(zrmu0,zfract,
+c      CALL SW_venus_dc_1Dglobave(zrmu0,zfract,   ! pour moy globale
+c      CALL SW_venus_dc(zrmu0,zfract,
+      CALL SW_venus_rh(zrmu0,zfract,latdeg,
+c      CALL SW_venus_rh_1Dglobave(zrmu0,zfract,   ! pour moy globale
      S        PPB,temp, 
      S        zheat, 
Index: trunk/LMDZ.VENUS/libf/phyvenus/sw_venus_rh.F
===================================================================
--- trunk/LMDZ.VENUS/libf/phyvenus/sw_venus_rh.F	(revision 1591)
+++ trunk/LMDZ.VENUS/libf/phyvenus/sw_venus_rh.F	(revision 1591)
@@ -0,0 +1,221 @@
+      SUBROUTINE SW_venus_rh(PRMU0, PFRAC, latdeg,
+     S              PPB, pt,
+     S              PHEAT, 
+     S              PTOPSW,PSOLSW,ZFSNET)
+      
+      use dimphy
+      use cpdet_mod, only: cpdet
+      IMPLICIT none
+
+#include "YOMCST.h"
+C
+C     ------------------------------------------------------------------
+C
+C     PURPOSE.
+C     --------
+C
+c      this routine loads and interpolates the shortwave radiation
+c     fluxes taken from Rainer Haus calculations for Venus.
+c     Ref: Haus et al. 2016
+C
+C     AUTHOR.
+C     -------
+C        Sebastien Lebonnois
+C
+C     MODIFICATIONS.
+C     --------------
+C        ORIGINAL : 5/2016
+C     ------------------------------------------------------------------
+C
+C* ARGUMENTS:
+C
+c inputs
+
+      REAL   PRMU0  ! COSINE OF ZENITHAL ANGLE
+      REAL   PFRAC  ! fraction de la journee
+      REAL   latdeg ! |latitude| (in degrees)
+      REAL   PPB(klev+1)  ! inter-couches PRESSURE (bar)
+      REAL   pt(klev)     ! mid-layer temperature
+C
+c output
+
+      REAL   PHEAT(klev)  ! SHORTWAVE HEATING (K/s) within each layer
+      REAL   PTOPSW       ! SHORTWAVE FLUX AT T.O.A. (net)
+      REAL   PSOLSW       ! SHORTWAVE FLUX AT SURFACE (net)
+      REAL   ZFSNET(klev+1) ! net solar flux at ppb levels
+
+C
+C* LOCAL VARIABLES:
+C
+      integer nlrh,nszarh,nlatrh
+      parameter (nlrh=118)  ! fichiers Rainer Haus
+      parameter (nszarh=7) ! fichiers Rainer Haus
+      parameter (nlatrh=19) ! fichiers Rainer Haus
+      
+      integer i,j,lat,nsza,nsza0(2),nl0,nlat0
+      real   zsnet(nlrh+1,nszarh+1,nlatrh+1)! net solar flux (W/m**2) (+ vers bas)
+      real   solza(nszarh,nlatrh)       ! solar zenith angles in table
+      real   presrh(nlrh+1)             ! pressure in table (bar)
+      real   altrh(nlrh+1)              ! altitude in table (km)
+      real   latrh(nlatrh)              ! latitude in table (degrees)
+      character*22 nullchar
+      real   sza0,factsza(2),factflux,factlat
+      real   zsnetmoy
+      logical firstcall
+      data    firstcall/.true./
+      save   solza,zsnet,altrh,latrh,presrh
+      save   firstcall
+      
+c ------------------------
+c Loading the file
+c ------------------------
+
+      if (firstcall) then
+
+       zsnet=0.
+
+       open(11,file='SolarNetFlux_RH.dat')
+
+       do i=1,nlrh+1
+          read(11,'(E5.1,4x,F8.2)') altrh(i),presrh(i)
+       enddo
+
+       do lat=1,nlatrh
+         latrh(lat)=5.*(lat-1)
+         read(11,*) nullchar
+         read(11,*) nullchar
+         read(11,'(3x,7(5x,E8.5))') solza(:,lat)
+         read(11,*) nullchar
+
+         do i=1,nlrh+1
+          read(11,'(E6.1,7(2x,F11.5),7x,F11.5)')
+     .          altrh(i),zsnet(i,1:nszarh,lat),zsnetmoy
+         enddo
+         read(11,*) nullchar
+       enddo
+       latrh(nlatrh)=89.
+
+c Correction of factor 2 in the table...
+       zsnet=zsnet*2.
+
+       close(11)
+
+       firstcall=.false.
+      endif
+
+c --------------------------------------
+c Interpolation in the GCM vertical grid
+c --------------------------------------
+
+c Latitude
+c ---------
+      
+      do lat=1,nlatrh
+         if (latrh(lat).le.latdeg) then
+              nlat0 = lat+1
+         endif
+      enddo
+      
+      if (nlat0.ne.nlatrh+1) then
+        factlat = (latdeg-latrh(nlat0-1))/(latrh(nlat0)-latrh(nlat0-1))
+      else
+        factlat = min((latdeg-latrh(nlatrh))/(90.-latrh(nlatrh)), 1.)
+      endif
+
+c Zenith angle
+c ------------
+      
+      sza0 = acos(PRMU0)/3.1416*180.
+c        print*,'Angle Zenithal =',sza0,' PFRAC=',PFRAC
+      nsza0(:)=2
+
+      do nsza=1,nszarh
+         if (solza(nsza,nlat0-1).le.sza0) then
+              nsza0(1) = nsza+1
+         endif
+      enddo
+      
+      if (nsza0(1).ne.nszarh+1) then
+          factsza(1) = (sza0-solza(nsza0(1)-1,nlat0-1))/
+     .        (solza(nsza0(1),nlat0-1)-solza(nsza0(1)-1,nlat0-1))
+      else
+          factsza(1) = min((sza0-solza(nszarh,nlat0-1))/
+     .         (90.-solza(nszarh,nlat0-1)), 1.)
+      endif
+
+      if (nlat0.ne.nlatrh+1) then
+       do nsza=1,nszarh
+         if (solza(nsza,nlat0).le.sza0) then
+              nsza0(2) = nsza+1
+         endif
+       enddo
+      
+       if (nsza0(2).eq.nszarh+1) then
+          factsza(2) = min((sza0-solza(nszarh,nlat0))/
+     .         (90.-solza(nszarh,nlat0)), 1.)
+       elseif ((nsza0(2).eq.2).and.(solza(1,nlat0).gt.sza0)) then
+          factsza(2) = 0.
+       else
+          factsza(2) = (sza0-solza(nsza0(2)-1,nlat0))/
+     .        (solza(nsza0(2),nlat0)-solza(nsza0(2)-1,nlat0))
+       endif
+      else
+        nsza0(2)   = nszarh+1
+        factsza(2) = 1.
+      endif
+
+c Pressure levels
+c ---------------
+
+      do j=1,klev+1
+        nl0 = nlrh
+        do i=nlrh+1,2,-1
+           if (presrh(i).ge.PPB(j)) then
+                nl0 = i-1
+           endif
+        enddo
+
+        factflux = (log10(max(PPB(j),presrh(1)))-log10(presrh(nl0+1)))
+     .            /(log10(presrh(nl0))-log10(presrh(nl0+1)))
+
+        ZFSNET(j) =  factlat*(
+     .      factflux   *  factsza(2)   *zsnet(nl0,nsza0(2),nlat0)
+     . +   factflux   *(1.-factsza(2))*zsnet(nl0,nsza0(2)-1,nlat0)
+     . + (1.-factflux)*  factsza(2)   *zsnet(nl0+1,nsza0(2),nlat0)
+     . + (1.-factflux)*(1.-factsza(2))*zsnet(nl0+1,nsza0(2)-1,nlat0) )
+     .            + (1.-factlat)*(
+     .      factflux   *  factsza(1)   *zsnet(nl0,nsza0(1),nlat0-1)
+     . +   factflux   *(1.-factsza(1))*zsnet(nl0,nsza0(1)-1,nlat0-1)
+     . + (1.-factflux)*  factsza(1)   *zsnet(nl0+1,nsza0(1),nlat0-1)
+     . + (1.-factflux)*(1.-factsza(1))*zsnet(nl0+1,nsza0(1)-1,nlat0-1) )
+
+        ZFSNET(j) = ZFSNET(j)*PFRAC
+        
+      enddo
+
+      PTOPSW = ZFSNET(klev+1)
+      PSOLSW = ZFSNET(1) 
+      
+c Heating rates
+c -------------
+c On utilise le gradient du flux pour calculer le taux de chauffage:
+c   heat(K/s) = d(fluxnet)  (W/m2)
+c              *g           (m/s2)
+c              /(-dp)  (epaisseur couche, en Pa=kg/m/s2)
+c              /cp  (J/kg/K) 
+
+      do j=1,klev
+! ADAPTATION GCM POUR CP(T)
+        PHEAT(j) = (ZFSNET(j+1)-ZFSNET(j))
+     .            *RG/cpdet(pt(j)) / ((PPB(j)-PPB(j+1))*1.e5)
+c-----TEST-------
+c tayloring the solar flux...
+        if ((PPB(j).gt.1.4).and.(PPB(j).le.10.)) then
+          PHEAT(j) = PHEAT(j)*3
+        endif
+c----------------
+      enddo
+
+      return
+      end
+
Index: trunk/LMDZ.VENUS/libf/phyvenus/sw_venus_rh_1Dglobave.F
===================================================================
--- trunk/LMDZ.VENUS/libf/phyvenus/sw_venus_rh_1Dglobave.F	(revision 1591)
+++ trunk/LMDZ.VENUS/libf/phyvenus/sw_venus_rh_1Dglobave.F	(revision 1591)
@@ -0,0 +1,175 @@
+      SUBROUTINE SW_venus_rh_1Dglobave(PRMU0, PFRAC,
+     S              PPB, pt,
+     S              PHEAT, 
+     S              PTOPSW,PSOLSW,ZFSNET)
+      
+      use dimphy
+      use cpdet_mod, only: cpdet
+      IMPLICIT none
+
+#include "YOMCST.h"
+C
+C     ------------------------------------------------------------------
+C
+C     PURPOSE.
+C     --------
+C
+c      this routine loads and interpolates the shortwave radiation
+c     fluxes taken from Rainer Haus calculations for Venus.
+c     Ref: Haus et al. 2016
+C
+C     AUTHOR.
+C     -------
+C        Sebastien Lebonnois
+C
+C     MODIFICATIONS.
+C     --------------
+C        ORIGINAL : 5/2016
+C     ------------------------------------------------------------------
+C
+C* ARGUMENTS:
+C
+c inputs
+
+      REAL   PRMU0  ! COSINE OF ZENITHAL ANGLE
+      REAL   PFRAC  ! fraction de la journee
+      REAL   PPB(klev+1)  ! inter-couches PRESSURE (bar)
+      REAL   pt(klev)     ! mid-layer temperature
+C
+c output
+
+      REAL   PHEAT(klev)  ! SHORTWAVE HEATING (K/s) within each layer
+      REAL   PTOPSW       ! SHORTWAVE FLUX AT T.O.A. (net)
+      REAL   PSOLSW       ! SHORTWAVE FLUX AT SURFACE (net)
+      REAL   ZFSNET(klev+1) ! net solar flux at ppb levels
+
+C
+C* LOCAL VARIABLES:
+C
+      integer nlrh,nszarh,nlatrh
+      parameter (nlrh=118)  ! fichiers Rainer Haus
+      parameter (nszarh=7) ! fichiers Rainer Haus
+      parameter (nlatrh=19) ! fichiers Rainer Haus
+      
+      integer i,j,lat,nsza,nsza0,nl0
+      real   zsnetmoy(nlrh+1,nlatrh)    ! net solar flux (W/m**2) (+ vers bas)
+      real   presrh(nlrh+1)             ! pressure in table (bar)
+      real   altrh(nlrh+1)              ! altitude in table (km)
+      real   latrh(nlatrh)              ! latitude in table (degrees)
+      real   zsolnet(nlrh+1)            ! for mean net solar flux in RH
+      character*22 nullchar
+      real   factflux
+      real   zsnet(nszarh)    ! net solar flux (W/m**2) (+ vers bas)
+      real   deltalat
+      logical firstcall
+      data    firstcall/.true./
+      save   zsolnet,altrh,presrh
+      save   firstcall
+      
+c ------------------------
+c Loading the file
+c ------------------------
+
+      if (firstcall) then
+
+       open(11,file='SolarNetFlux_RH.dat')
+
+       do i=1,nlrh+1
+          read(11,'(E5.1,4x,F8.2)') altrh(i),presrh(i)
+       enddo
+
+       do lat=1,nlatrh
+         latrh(lat)=5.*(lat-1)
+         read(11,*) nullchar
+         read(11,*) nullchar
+         read(11,*) nullchar
+         read(11,*) nullchar
+
+         do i=1,nlrh+1
+          read(11,'(E6.1,7(2x,F11.5),7x,F11.5)')
+     .          altrh(i),zsnet,zsnetmoy(i,lat)
+         enddo
+         read(11,*) nullchar
+       enddo
+       latrh(nlatrh)=89.
+
+       close(11)
+
+c ----------- TEST ------------
+c      Moyenne planetaire
+c -----------------------------
+
+      zsolnet=0.
+      do lat=1,nlatrh-1
+        deltalat=(latrh(lat+1)-latrh(lat))*RPI/180.
+        do j=1,nlrh+1
+        zsolnet(j) = zsolnet(j)+
+     .         (zsnetmoy(j,lat+1)+zsnetmoy(j,lat))/2.*
+     .       deltalat*cos((latrh(lat+1)+latrh(lat))*RPI/360.)
+        enddo
+      enddo
+c -----------------------------
+c --------  FIN TEST ----------
+
+       firstcall=.false.
+      endif
+
+c --------------------------------------
+c Interpolation in the GCM vertical grid
+c --------------------------------------
+
+c Pressure levels
+c ---------------
+
+      do j=1,klev+1
+        nl0 = nlrh
+        do i=nlrh+1,2,-1
+           if (presrh(i).ge.PPB(j)) then
+                nl0 = i-1
+           endif
+        enddo
+        
+        factflux = (log10(max(PPB(j),presrh(1)))-log10(presrh(nl0+1)))
+     .            /(log10(presrh(nl0))-log10(presrh(nl0+1)))
+
+        ZFSNET(j) =  factflux     *zsolnet(nl0)
+     .             + (1.-factflux)*zsolnet(nl0+1)
+        
+c-----TEST-------
+c tayloring the solar flux...
+c        if ((PPB(j).gt.0.236).and.(PPB(j).le.22.52)) then
+c         ZFSNET(j) = ZFSNET(j)+2.5*(1.+cos((log10(PPB(j)/3.5)/
+c     .                                  log10(0.236/3.5))*RPI))
+c        endif
+c----------------
+      enddo
+
+      PTOPSW = ZFSNET(klev+1)
+      PSOLSW = ZFSNET(1) 
+      
+c Heating rates
+c -------------
+c On utilise le gradient du flux pour calculer le taux de chauffage:
+c   heat(K/s) = d(fluxnet)  (W/m2)
+c              *g           (m/s2)
+c              /(-dp)  (epaisseur couche, en Pa=kg/m/s2)
+c              /cp  (J/kg/K) 
+
+      do j=1,klev
+! ADAPTATION GCM POUR CP(T)
+        PHEAT(j) = (ZFSNET(j+1)-ZFSNET(j))
+     .            *RG/cpdet(pt(j)) / ((PPB(j)-PPB(j+1))*1.e5)
+c-----TEST-------
+c tayloring the solar flux...
+        if ((PPB(j).gt.1.4).and.(PPB(j).le.10.)) then
+          PHEAT(j) = PHEAT(j)*3
+c       elseif ((PPB(j).gt.10.).and.(PPB(j).le.30.)) then
+c         PHEAT(j) = PHEAT(j)*1.5
+        endif
+c----------------
+c       print*,PPB(j),ZFSNET(j),PHEAT(j)
+      enddo
+c     stop
+      return
+      end
+
