Index: trunk/LMDZ.VENUS/libf/phyvenus/YOEGWD.h
===================================================================
--- trunk/LMDZ.VENUS/libf/phyvenus/YOEGWD.h	(revision 175)
+++ trunk/LMDZ.VENUS/libf/phyvenus/YOEGWD.h	(revision 778)
@@ -1,10 +1,10 @@
-C     -----------------------------------------------------------------
-C*    *COMMON* *YOEGWD* - PARAMETERS FOR GRAVITY WAVE DRAG CALCULATIONS
-C     -----------------------------------------------------------------
-C
-      integer NKTOPG,NTOP
-      real GFRCRIT,GKWAKE,GRCRIT,GVCRIT,GKDRAG,GKLIFT
-      real GHMAX,GRAHILO,GSIGCR,GSSEC,GTSEC,GVSEC
-      COMMON/YOEGWD/ GFRCRIT,GKWAKE,GRCRIT,GVCRIT,GKDRAG,GKLIFT
-     *        ,GHMAX,GRAHILO,GSIGCR,NKTOPG,NTOP,GSSEC,GTSEC,GVSEC
+!     -----------------------------------------------------------------
+!*    *COMMON* *YOEGWD* - PARAMETERS FOR GRAVITY WAVE DRAG CALCULATIONS
+!     -----------------------------------------------------------------
+!
+      integer :: NKTOPG,NTOP
+      real    :: GFRCRIT,GKWAKE,GRCRIT,GVCRIT,GKDRAG,GKLIFT
+      real    :: GHMAX,GRAHILO,GSIGCR,GSSEC,GTSEC,GVSEC
+      COMMON/YOEGWD/ GFRCRIT,GKWAKE,GRCRIT,GVCRIT,GKDRAG,GKLIFT  &
+         ,GHMAX,GRAHILO,GSIGCR,NKTOPG,NTOP,GSSEC,GTSEC,GVSEC
 
Index: trunk/LMDZ.VENUS/libf/phyvenus/flott_gwd_ran.F
===================================================================
--- trunk/LMDZ.VENUS/libf/phyvenus/flott_gwd_ran.F	(revision 175)
+++ 	(revision )
@@ -1,407 +1,0 @@
-      SUBROUTINE FLOTT_GWD_RAN(NLON,NLEV,DTIME, pp, pn2,
-     I      tt,uu,vv,
-     O      zustr,zvstr,
-     O      d_t, d_u, d_v)
-
-C ########################################################################
-C Parametrization of the momentum flux deposition due to a discrete
-C numbe of gravity waves. 
-C    F. Lott (2008)
-C                 LMDz model stand alone version
-C ########################################################################
-C 
-C
-      use dimphy
-      implicit none
-
-#include "dimensions.h"
-#include "paramet.h"
-
-#include "YOEGWD.h"
-#include "YOMCST.h"
-
-c VENUS ATTENTION: CP VARIABLE
-      real cp(NLON,NLEV)
-
-C  DECLARATIONS:
-C  INPUTS
-      INTEGER NLON,NLEV 
-      REAL DTIME
-      REAL pp(NLON,NLEV)    ! Pressure at full levels
-c VENUS ATTENTION: CP VARIABLE PN2 CALCULE EN AMONT DES PARAMETRISATIONS
-      REAL pn2(NLON,NLEV)   ! N2 (BV^2) at 1/2 levels
-      REAL TT(NLON,NLEV)                 ! Temp at full levels 
-      REAL UU(NLON,NLEV),VV(NLON,NLEV) ! Temp, and winds at full levels
-
-C  OUTPUTS
-      REAL zustr(NLON),zvstr(NLON)          ! Surface Stresses
-      REAL d_t(NLON,NLEV)                   ! Tendency on Temp.
-      REAL d_u(NLON,NLEV), d_v(NLON,NLEV)   ! Tendencies on winds
-
-c ON CONSERVE LA MEMOIRE un certain temps AVEC UN SAVE ET UN COEF
-      real coef
-      parameter (coef = 0.986)
-      real,save,allocatable :: d_u_sav(:,:),d_v_sav(:,:)
-      LOGICAL firstcall
-      SAVE firstcall
-      DATA firstcall/.true./
-
-C  GRAVITY-WAVES SPECIFICATIONS
-
-      INTEGER JK,NK,JP,NP,JO,NO,JW,NW,II,LL
-      PARAMETER(NK=4,NP=4,NO=4,NW=NK*NP*NO) ! nombres d'ondes envoyées 
-      REAL KMIN,KMAX   !  Min and Max horizontal wavenumbers
-      REAL OMIN,OMAX   ! Max frequency and security for divisions!
-      PARAMETER (KMIN=2.E-5,KMAX=2.E-4,OMIN=2.E-5,OMAX=2.E-3)
-      REAL RUWMAX,SAT  ! MAX EP FLUX AT LAUNCH LEVEL, SATURATION PARAMETER
-      PARAMETER (RUWMAX=1.E-3,SAT=0.1)
-      REAL ZK(NW,KLON),ZP(NW),ZO(NW,KLON)  ! Waves absolute specifications!
-      REAL ZOM(NW,KLON),ZOP(NW,KLON)  ! Intrisic frequency at 1/2 levels
-      REAL PHM(NW,KLON),PHP(NW,KLON)  ! Phi at 1/2 levels
-      REAL RUW0(NW,KLON)              ! Fluxes at lauch level
-      REAL RUWP(NW,KLON),RVWP(NW,KLON)! Fluxes X and Y for each wave
-      INTEGER LAUNCH
-      REAL ALEAS
-      EXTERNAL ALEAS
-
-C  OUTPUT FOR RAW PROFILES:
-
-      INTEGER IU
-
-C  PARAMETERS OF WAVES DISSPATIONS
-      REAL RDISS,ZOISEC  !  COEFF DE DISSIPATION, SECURITY FOR INTRINSIC FREQ
-      PARAMETER (RDISS=10.,ZOISEC=1.E-5)
-
-C  INTERNAL ARRAYS
-      real PR,RH(KLON,KLEV+1) ! Ref Press. intervient dans amplitude ?
-      PARAMETER (PR=9.2e6)    ! VENUS!!
-
-
-      REAL UH(KLON,KLEV+1),VH(KLON,KLEV+1) ! Winds at 1/2 levels
-      REAL PH(KLON,KLEV+1)                 ! Pressure at 1/2 levels
-      REAL PHM1(KLON,KLEV+1)               ! 1/Press  at 1/2 levels
-      REAL BV(KLON,KLEV+1),BVSEC           ! BVF at 1/2 levels
-      REAL PSEC                            ! SECURITY TO AVOID DIVISION
-      PARAMETER (PSEC=1.E-6)               ! BY O PRESSURE VALUE
-      REAL RUW(KLON,KLEV+1)                ! Flux x at semi levels
-      REAL RVW(KLON,KLEV+1)                ! Flux y at semi levels
-
-C  INITIALISATION
-
-      IF (firstcall) THEN
-        allocate(d_u_sav(klon,klev),d_v_sav(klon,klev))
-        firstcall=.false.
-      ENDIF
-
-      BVSEC=1.E-6
-
-c      PRINT *,ALEAS(0.)
-
-C  WAVES CHARACTERISTTICS
-
-      JW=0
-      DO 11 JP=1,NP
-      DO 11 JK=1,NK
-      DO 11 JO=1,NO
-        JW=JW+1
-      ZP(JW)=2.*RPI*FLOAT(JP-1)/FLOAT(NP) ! Angle
-      DO 11 II=1,NLON
-        ZK(JW,II)=KMIN+(KMAX-KMIN)*ALEAS(0.)  ! Hor. Wavenbr
-        ZO(JW,II)=OMIN+(OMAX-OMIN)*ALEAS(0.)  ! Absolute frequency
-        RUW0(JW,II)=RUWMAX/FLOAT(NW)*ALEAS(0.) ! Momentum flux at launch lev     
-11    CONTINUE
-     
-c      print*,"cos(ZP)=",cos(ZP)
-
-c        PRINT *,'Catch 22'
-
-C  2. SETUP
-
-c  VENUS: CP(T) at ph levels
-      do II=1,NLON
-        do LL=2,KLEV
-	  cp(II,LL)=cpdet(0.5*(TT(II,LL)+TT(II,LL-1)))
-        enddo
-      enddo
-C  PRESSURE AT HALPH LEVELS
-
-        II=KLON/2
-
-        DO 20 LL=2,KLEV
-        PH(:,LL)=EXP((ALOG(PP(:,LL))+ALOG(PP(:,LL-1)))/2.)
-        PHM1(:,LL)=1./PH(:,LL)
-20      CONTINUE
-
-        PH(:,KLEV+1)=0.  
-        PHM1(:,KLEV+1)=1./PSEC
-
-        PH(:,1)=2.*PP(:,1)-PH(:,2)
-
-C  LOG-PRESSURE ALTITUDE AND DENSITY
-
-
-        DO 22 LL=1,KLEV+1  
-          RH(:,LL)=PH(:,LL)/PR           !  Reference density
-22      CONTINUE
-
-C   LAUNCHING ALTITUDE
-
-      DO LL=1,NLEV
-      IF(PH(NLON/2,LL)/PH(NLON/2,1).GT.0.8)LAUNCH=LL
-      ENDDO
-c      PRINT *,'LAUNCH:',LAUNCH
-
-C  WINDS AND BV FREQUENCY AT 1/2 LEVELS
-
-        DO 23 LL=2,KLEV
-        UH(:,LL)=0.5*(UU(:,LL)+UU(:,LL-1))
-        VH(:,LL)=0.5*(VV(:,LL)+VV(:,LL-1))
-
-c VENUS ATTENTION: CP VARIABLE PSTAB CALCULE EN AMONT DES PARAMETRISATIONS
-        BV(:,LL)=SQRT(AMAX1(BVSEC,pn2(:,LL)))
-
-23      CONTINUE
-	
-c       print *,'catch 23'
-
-        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)
-
-C  3. COMPUTE THE FLUXES!!
-
-        DO  3 JW=1,NW
-
-        ZOP(JW,:)=ZO(JW,:)-ZK(JW,:)*COS(ZP(JW))*UH(:,LAUNCH) ! Intrinsic
-     C                    -ZK(JW,:)*SIN(ZP(JW))*VH(:,LAUNCH) ! Frequency
-        PHP(JW,:)=
-C       AMIN1(                          ! Phi at launch level
-C  ENSURES THE IMPOSED MOM FLUX:
-     C SQRT(ABS(ZOP(JW,:))*BV(:,LAUNCH)
-     C      *RUW0(JW,:)/RH(:,LAUNCH))/ZK(JW,:)
-        RUWP(JW,:)=COS(ZP(JW))*SIGN(1.,ZOP(JW,:))*RUW0(JW,:)
-        RVWP(JW,:)=SIN(ZP(JW))*SIGN(1.,ZOP(JW,:))*RUW0(JW,:)
-
-3       CONTINUE
-
-        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,:)
-         enddo
-	enddo
-
-        DO 33 LL=LAUNCH,KLEV-1
-
-        DO  34 JW=1,NW
-          ZOM(JW,:)=ZOP(JW,:)
-          PHM(JW,:)=PHP(JW,:)
-
-        ZOP(JW,:)=ZO(JW,:)-ZK(JW,:)*COS(ZP(JW))*UH(:,LL+1) ! Intrinsic
-     C                  -ZK(JW,:)*SIN(ZP(JW))*VH(:,LL+1) ! Frequency
-C
-        PHP(JW,:)=           
-     C             AMIN1(                             
-C  NO BREAKING
-     C               SQRT(PH(:,LL)*PHM1(:,LL+1))*PHM(JW,:)
-     C              *SQRT(BV(:,LL+1)/BV(:,LL)
-     C              *ABS(ZOP(JW,:))/AMAX1(ABS(ZOM(JW,:)),ZOISEC))
-C  CRITICAL LEVELS
-     C           ,
-     C              AMAX1(0.,SIGN(1.,ZOP(JW,:)*ZOM(JW,:)))
-C  SATURATION
-     C              *ZOP(JW,:)**2/ZK(JW,:)**2/SQRT(FLOAT(NK*NO))*SAT
-     C           )
-
-34      CONTINUE
-
-        DO 35 JW=1,NW
-        RUWP(JW,:)=ZOP(JW,:)*COS(ZP(JW))*ZK(JW,:)**2
-     C      *PH(:,LL+1)/PR*PHP(JW,:)**2/BV(:,LL+1)
-     C      /AMAX1(ABS(ZOP(JW,:)),ZOISEC)**2
-        RVWP(JW,:)=ZOP(JW,:)*SIN(ZP(JW))*ZK(JW,:)**2
-     C      *PH(:,LL+1)/PR*PHP(JW,:)**2/BV(:,LL+1)
-     C      /AMAX1(ABS(ZOP(JW,:)),ZOISEC)**2
- 35     CONTINUE
-
-        RUW(:,LL+1)=0.
-        RVW(:,LL+1)=0.
-
-
-        DO 36 JW=1,NW
-        RUW(:,LL+1)=RUW(:,LL+1)+RUWP(JW,:)           
-        RVW(:,LL+1)=RVW(:,LL+1)+RVWP(JW,:)           
- 36     CONTINUE
-
-
-33      CONTINUE
-
-C  4 CALCUL DES TENDANCES:
-
-C  RECTIFICATION AU SOMMET ET DANS LES BASSES COUCHES:
-
-        RUW(:,1)=RUW(:,LAUNCH)
-        RVW(:,1)=RVW(:,LAUNCH)
-        DO 4 LL=2,LAUNCH
-        RUW(:,LL)=RUW(:,LL-1)+(RUW(:,LAUNCH+1)-RUW(:,1))*
-     C            (PH(:,LL)-PH(:,LL-1))/(PH(:,LAUNCH+1)-PH(:,1))
-        RVW(:,LL)=RVW(:,LL-1)+(RVW(:,LAUNCH+1)-RVW(:,1))*
-     C            (PH(:,LL)-PH(:,LL-1))/(PH(:,LAUNCH+1)-PH(:,1))
- 4      CONTINUE
-        
-        DO 41 LL=1,KLEV
-        D_U(:,LL)=RG*(RUW(:,LL+1)-RUW(:,LL))
-     C              /(PH(:,LL+1)-PH(:,LL))*DTIME
-        D_V(:,LL)=RG*(RVW(:,LL+1)-RVW(:,LL))
-     C              /(PH(:,LL+1)-PH(:,LL))*DTIME
- 41     CONTINUE
- 
-c ON CONSERVE LA MEMOIRE un certain temps AVEC UN SAVE ET UN COEF
-        d_u = d_u + coef * d_u_sav
-        d_v = d_v + coef * d_v_sav
-	d_u_sav = d_u
-	d_v_sav = d_v
-
-c        PRINT *,'OUT LOTT, LONG?'
-c        DO JW=1,NW
-c        PRINT *,'ONDE ',COS(ZP(JW))*JW,' KM:',2*RPI/ZK(JW)/1000.
-c     C         ,' HR:',2*RPI/ZO(JW)/3600.,' C:',ZO(JW)/ZK(JW)
-c        ENDDO
-
-        RETURN
-        END
-
-
-c===================================================================
-c===================================================================
-c===================================================================
-c===================================================================
-
-      FUNCTION ALEAS (R)
-C***BEGIN PROLOGUE  ALEAS
-C***PURPOSE  Generate a uniformly distributed random number.
-C***LIBRARY   SLATEC (FNLIB)
-C***CATEGORY  L6A21
-C***TYPE      SINGLE PRECISION (ALEAS-S)
-C***KEYWORDS  FNLIB, ALEAS NUMBER, SPECIAL FUNCTIONS, UNIFORM
-C***AUTHOR  Fullerton, W., (LANL)
-C***DESCRIPTION
-C
-C      This pseudo-random number generator is portable among a wide
-C variety of computers.  RAND(R) undoubtedly is not as good as many
-C readily available installation dependent versions, and so this
-C routine is not recommended for widespread usage.  Its redeeming
-C feature is that the exact same random numbers (to within final round-
-C off error) can be generated from machine to machine.  Thus, programs
-C that make use of random numbers can be easily transported to and
-C checked in a new environment.
-C
-C      The random numbers are generated by the linear congruential
-C method described, e.g., by Knuth in Seminumerical Methods (p.9),
-C Addison-Wesley, 1969.  Given the I-th number of a pseudo-random
-C sequence, the I+1 -st number is generated from
-C             X(I+1) = (A*X(I) + C) MOD M,
-C where here M = 2**22 = 4194304, C = 1731 and several suitable values
-C of the multiplier A are discussed below.  Both the multiplier A and
-C random number X are represented in double precision as two 11-bit
-C words.  The constants are chosen so that the period is the maximum
-C possible, 4194304.
-C
-C      In order that the same numbers be generated from machine to
-C machine, it is necessary that 23-bit integers be reducible modulo
-C 2**11 exactly, that 23-bit integers be added exactly, and that 11-bit
-C integers be multiplied exactly.  Furthermore, if the restart option
-C is used (where R is between 0 and 1), then the product R*2**22 =
-C R*4194304 must be correct to the nearest integer.
-C
-C      The first four random numbers should be .0004127026,
-C .6750836372, .1614754200, and .9086198807.  The tenth random number
-C is .5527787209, and the hundredth is .3600893021 .  The thousandth
-C number should be .2176990509 .
-C
-C      In order to generate several effectively independent sequences
-C with the same generator, it is necessary to know the random number
-C for several widely spaced calls.  The I-th random number times 2**22,
-C where I=K*P/8 and P is the period of the sequence (P = 2**22), is
-C still of the form L*P/8.  In particular we find the I-th random
-C number multiplied by 2**22 is given by
-C I   =  0  1*P/8  2*P/8  3*P/8  4*P/8  5*P/8  6*P/8  7*P/8  8*P/8
-C RAND=  0  5*P/8  2*P/8  7*P/8  4*P/8  1*P/8  6*P/8  3*P/8  0
-C Thus the 4*P/8 = 2097152 random number is 2097152/2**22.
-C
-C      Several multipliers have been subjected to the spectral test
-C (see Knuth, p. 82).  Four suitable multipliers roughly in order of
-C goodness according to the spectral test are
-C    3146757 = 1536*2048 + 1029 = 2**21 + 2**20 + 2**10 + 5
-C    2098181 = 1024*2048 + 1029 = 2**21 + 2**10 + 5
-C    3146245 = 1536*2048 +  517 = 2**21 + 2**20 + 2**9 + 5
-C    2776669 = 1355*2048 + 1629 = 5**9 + 7**7 + 1
-C
-C      In the table below LOG10(NU(I)) gives roughly the number of
-C random decimal digits in the random numbers considered I at a time.
-C C is the primary measure of goodness.  In both cases bigger is better.
-C
-C                   LOG10 NU(I)              C(I)
-C       A       I=2  I=3  I=4  I=5    I=2  I=3  I=4  I=5
-C
-C    3146757    3.3  2.0  1.6  1.3    3.1  1.3  4.6  2.6
-C    2098181    3.3  2.0  1.6  1.2    3.2  1.3  4.6  1.7
-C    3146245    3.3  2.2  1.5  1.1    3.2  4.2  1.1  0.4
-C    2776669    3.3  2.1  1.6  1.3    2.5  2.0  1.9  2.6
-C   Best
-C    Possible   3.3  2.3  1.7  1.4    3.6  5.9  9.7  14.9
-C
-C             Input Argument --
-C R      If R=0., the next random number of the sequence is generated.
-C        If R .LT. 0., the last generated number will be returned for
-C          possible use in a restart procedure.
-C        If R .GT. 0., the sequence of random numbers will start with
-C          the seed R mod 1.  This seed is also returned as the value of
-C          RAND provided the arithmetic is done exactly.
-C
-C             Output Value --
-C RAND   a pseudo-random number between 0. and 1.
-C
-C***REFERENCES  (NONE)
-C***ROUTINES CALLED  (NONE)
-C***REVISION HISTORY  (YYMMDD)
-C   770401  DATE WRITTEN
-C   890531  Changed all specific intrinsics to generic.  (WRB)
-C   890531  REVISION DATE from Version 3.2
-C   891214  Prologue converted to Version 4.0 format.  (BAB)
-C***END PROLOGUE  RAND
-      SAVE IA1, IA0, IA1MA0, IC, IX1, IX0
-      DATA IA1, IA0, IA1MA0 /1536, 1029, 507/
-      DATA IC /1731/
-      DATA IX1, IX0 /0, 0/
-C***FIRST EXECUTABLE STATEMENT  RAND
-      IF (R.LT.0.) GO TO 10
-      IF (R.GT.0.) GO TO 20
-C
-C           A*X = 2**22*IA1*IX1 + 2**11*(IA1*IX1 + (IA1-IA0)*(IX0-IX1)
-C                   + IA0*IX0) + IA0*IX0
-C
-      IY0 = IA0*IX0
-      IY1 = IA1*IX1 + IA1MA0*(IX0-IX1) + IY0
-      IY0 = IY0 + IC
-      IX0 = MOD (IY0, 2048)
-      IY1 = IY1 + (IY0-IX0)/2048
-      IX1 = MOD (IY1, 2048)
-C
- 10   ALEAS = IX1*2048 + IX0
-      ALEAS = ALEAS / 4194304.
-      RETURN
-C
- 20   IX1 = MOD(R,1.)*4194304. + 0.5
-      IX0 = MOD (IX1, 2048)
-      IX1 = (IX1-IX0)/2048
-      GO TO 10
-C
-      END
-
-
Index: trunk/LMDZ.VENUS/libf/phyvenus/flott_gwd_ran.F90
===================================================================
--- trunk/LMDZ.VENUS/libf/phyvenus/flott_gwd_ran.F90	(revision 778)
+++ trunk/LMDZ.VENUS/libf/phyvenus/flott_gwd_ran.F90	(revision 778)
@@ -0,0 +1,374 @@
+      SUBROUTINE FLOTT_GWD_RAN(NLON,NLEV,DTIME, pp, pn2,  &
+                  tt,uu,vv,zustr,zvstr,d_t, d_u, d_v)
+
+    !----------------------------------------------------------------------
+    ! 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
+    ! LMDz model online version      
+    ! ADAPTED FOR VENUS
+    !---------------------------------------------------------------------
+
+      use dimphy
+      implicit none
+
+#include "dimensions.h"
+#include "paramet.h"
+
+#include "YOEGWD.h"
+#include "YOMCST.h"
+
+    ! 0. DECLARATIONS:
+
+    ! 0.1 INPUTS
+    INTEGER, intent(in):: NLON, NLEV 
+    REAL, intent(in):: DTIME ! Time step of the Physics
+    REAL, intent(in):: pp(NLON, NLEV)   ! Pressure at full levels
+! VENUS ATTENTION: CP VARIABLE PN2 CALCULE EN AMONT DES PARAMETRISATIONS
+    REAL, intent(in):: pn2(NLON,NLEV)   ! N2 (BV^2) at 1/2 levels
+    REAL, intent(in):: TT(NLON, NLEV)   ! Temp at full levels 
+
+    REAL, intent(in):: UU(NLON, NLEV) , VV(NLON, NLEV)
+    ! Hor winds at full levels
+
+    ! 0.2 OUTPUTS
+    REAL, intent(out):: zustr(NLON), zvstr(NLON) ! Surface Stresses
+    REAL, intent(inout):: d_t(NLON, NLEV)        ! Tendency on Temp.
+
+    REAL, intent(inout):: d_u(NLON, NLEV), d_v(NLON, NLEV)
+    ! Tendencies on winds
+
+    ! O.3 INTERNAL ARRAYS
+
+    INTEGER II, LL, IEQ
+
+    ! 0.3.0 TIME SCALE OF THE LIFE CYCLE OF THE WAVES PARAMETERIZED
+
+    REAL DELTAT
+
+    ! 0.3.1 GRAVITY-WAVES SPECIFICATIONS
+
+!VENUS INTEGER, PARAMETER:: NK = 4, NP = 4, NO = 4, NW = NK * NP * NO
+    !Online output: change NO
+    INTEGER, PARAMETER:: NK = 1, NP = 2, NO = 10, NW = NK * NP * NO
+    INTEGER JK, JP, JO, JW
+    REAL KMIN, KMAX ! Min and Max horizontal wavenumbers
+    REAL CMIN, CMAX ! Min and Max absolute ph. vel.
+    REAL CPHA ! absolute PHASE VELOCITY frequency
+    REAL ZK(NW, KLON) ! Horizontal wavenumber amplitude
+    REAL ZP(NW)       ! Horizontal wavenumber angle        
+    REAL ZO(NW, KLON) ! Absolute frequency      !
+
+    ! Waves Intr. freq. at the 1/2 lev surrounding the full level
+    REAL ZOM(NW, KLON), ZOP(NW, KLON)
+
+    ! Wave vertical velocities at the 2 1/2 lev surrounding the full level
+    REAL WWM(NW, KLON), WWP(NW, KLON)
+
+    REAL RUW0(NW, KLON) ! Fluxes at launching level
+
+    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
+
+    ! 0.3.2 PARAMETERS OF WAVES DISSIPATIONS
+
+    REAL RDISS, ZOISEC ! COEFF DE DISSIPATION, SECURITY FOR INTRINSIC FREQ
+
+    ! 0.3.3 BACKGROUND FLOW AT 1/2 LEVELS AND VERTICAL COORDINATE
+
+    REAL H0(KLON, KLEV)     ! Characteristic Height of the atmosphere
+    REAL PR ! Reference Pressure 
+
+    REAL ZH(KLON, KLEV + 1) ! Log-pressure altitude
+
+    REAL UH(KLON, KLEV + 1), VH(KLON, KLEV + 1) ! Winds at 1/2 levels
+    REAL PH(KLON, KLEV + 1) ! Pressure at 1/2 levels
+    REAL PSEC ! Security to avoid division by 0 pressure
+    REAL BV(KLON, KLEV + 1) ! Brunt Vaisala freq. (BVF) at 1/2 levels
+    REAL BVSEC ! Security to avoid negative BVF
+
+! COSMETICS TO DIAGNOSE EACH WAVES CONTRIBUTION.
+    logical output
+    data output/.false./
+  ! CAUTION ! IF output is .true. THEN change NO to 10 at least !
+    character*14 outform
+    character*2  str2
+
+! ON CONSERVE LA MEMOIRE un certain temps AVEC UN SAVE
+    real,save,allocatable :: d_u_sav(:,:),d_v_sav(:,:)
+    LOGICAL firstcall
+    SAVE firstcall
+    DATA firstcall/.true./
+
+    !-----------------------------------------------------------------
+    !  1. INITIALISATIONS
+
+      IF (firstcall) THEN
+        allocate(d_u_sav(NLON,NLEV),d_v_sav(NLON,NLEV))
+        firstcall=.false.
+      ENDIF
+    
+    !    1.1 Basic parameter
+
+    !  PARAMETERS CORRESPONDING TO V3:
+    RUWMAX = 0.005      ! Max EP-Flux at Launch altitude
+    SAT    = 0.85       ! Saturation parameter: Sc in (12)
+    RDISS  = 10.        ! 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
+    !Online output: one value only
+    if (output) then
+      KMIN = 1.3E-5
+      KMAX = 1.3E-5
+    endif
+    CMIN = 1.           ! Min phase velocity
+    CMAX = 60.          ! Max phase speed velocity
+    XLAUNCH=0.6         ! Parameter that control launching altitude
+
+    PR = 9.2e6          ! Reference pressure    ! VENUS!!
+
+    BVSEC  = 1.E-5      ! Security to avoid negative BVF  
+    PSEC   = 1.E-6      ! Security to avoid division by 0 pressure
+    ZOISEC = 1.E-6      ! Security FOR 0 INTRINSIC FREQ
+
+    IF(DELTAT.LT.DTIME)THEN
+       PRINT *,'GWD RANDO: DELTAT LT DTIME!'
+       STOP
+    ENDIF
+
+
+    IF (NLEV < NW) THEN
+       PRINT *, 'YOU WILL HAVE PROBLEM WITH RANDOM NUMBERS'
+       PRINT *, 'FLOTT GWD STOP'
+       STOP 1
+    ENDIF
+
+    ! 1.2 WAVES CHARACTERISTICS CHOSEN RANDOMLY
+    !-------------------------------------------
+
+    ! The mod function of here a weird arguments
+    ! are used to produce the waves characteristics
+    ! in a stochastic way
+
+    JW = 0
+    DO JP = 1, NP
+       DO JK = 1, NK
+          DO JO = 1, NO
+             JW = JW + 1
+             !  Angle
+             ZP(JW) = 2. * RPI * REAL(JP - 1) / REAL(NP) 
+             DO II = 1, KLON
+                ! Horizontal wavenumber amplitude
+                ZK(JW, II) = KMIN + (KMAX - KMIN) * MOD(TT(II, JW) * 100., 1.)
+                ! Horizontal phase speed
+                CPHA = CMIN + (CMAX - CMIN) * MOD(TT(II, JW)**2, 1.)
+       !Online output: linear
+                if (output) CPHA = CMIN + (CMAX - CMIN) * JO/NO
+                ! Intrinsic frequency
+                ZO(JW, II) = CPHA * ZK(JW, II)
+                ! Momentum flux at launch lev 
+                ! RUW0(JW, II) = RUWMAX / REAL(NW) &
+                RUW0(JW, II) = RUWMAX &
+                     * MOD(100. * (UU(II, JW)**2 + VV(II, JW)**2), 1.)
+             ENDDO
+          end DO
+       end DO
+    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) 
+    ZH(:,1) = 0.    
+    DO LL = 2, KLEV + 1 
+       H0(:, LL-1) = RD * TT(:, LL-1) / RG 
+       ZH(:, LL) = ZH(:, LL-1) + H0(:, LL-1)*(PH(:, LL-1)-PH(:,LL))/PP(:, LL-1)
+    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) = SQRT(MAX(BVSEC,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
+    !--------------------------
+
+    !  3.1  Vertical velocity at launching altitude to ensure 
+    !       the correct value to the imposed fluxes.
+    !
+    DO JW = 1, NW
+
+       ! Evaluate intrinsic frequency at launching altitude:
+       ZOP(JW, :) = ZO(JW, :) &
+            - 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,:))
+       RUWP(JW, :) = COS(ZP(JW)) * SIGN(1., ZOP(JW, :)) * RUW0(JW, :)
+       RVWP(JW, :) = SIN(ZP(JW)) * SIGN(1., ZOP(JW, :)) * RUW0(JW, :)
+
+    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
+    !      next done by i) conserving the EP flux, ii) dissipating
+    !      a little, iii) testing critical levels, and vi) testing
+    !      the breaking.
+
+    !Online output
+    write(str2,'(i2)') NW+1
+    outform="("//str2//"(E12.4,1X))"
+    if (output) WRITE(11,outform) ZH(IEQ, 1) / 1000., (ZO(JW, IEQ)/ZK(JW, IEQ)*COS(ZP(JW)), JW = 1, NW)
+
+    DO LL = LAUNCH, KLEV - 1
+
+
+       !  W(KB)ARNING: ALL THE PHYSICS IS HERE (PASSAGE FROM ONE LEVEL
+       ! TO THE NEXT)
+       DO JW = 1, NW
+          ZOM(JW, :) = ZOP(JW, :)
+          WWM(JW, :) = WWP(JW, :)
+          ! Intrinsic Frequency
+          ZOP(JW, :) = ZO(JW, :) - ZK(JW, :) * COS(ZP(JW)) * UH(:, LL + 1) &
+               - ZK(JW, :) * SIN(ZP(JW)) * VH(:, LL + 1) 
+
+          WWP(JW, :) = MIN( & 
+               ! No breaking (Eq.6)
+               WWM(JW, :) & 
+               * SQRT(BV(:, LL ) / BV(:, LL+1) &
+               * ABS(ZOP(JW, :)) / MAX(ABS(ZOM(JW, :)), ZOISEC)) &
+               ! 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)
+               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(:,LL)) * SAT)
+       end DO
+
+       ! END OF W(KB)ARNING
+       ! Evaluate EP-flux from Eq. 7 and 
+       ! Give the right orientation to the stress
+
+       DO JW = 1, NW
+          RUWP(JW, :) = ZOP(JW, :)/MAX(ABS(ZOP(JW, :)), ZOISEC)**2 &
+               *BV(:,LL+1)&
+               * COS(ZP(JW)) *  MAX(WWP(JW, :),1e-30)**2
+          RVWP(JW, :) = ZOP(JW, :)/MAX(ABS(ZOP(JW, :)), ZOISEC)**2 &
+               *BV(:,LL+1)&
+               * SIN(ZP(JW)) *  MAX(WWP(JW, :),1e-30)**2
+       end DO
+       !
+       RUW(:, LL + 1) = 0.
+       RVW(:, LL + 1) = 0.
+
+       DO JW = 1, NW
+          RUW(:, LL + 1) = RUW(:, LL + 1) + RUWP(JW, :) 
+          RVW(:, LL + 1) = RVW(:, LL + 1) + RVWP(JW, :) 
+       end DO
+       !Online output
+        if (output)  WRITE(11,outform) ZH(IEQ, LL + 1) / 1000., (RUWP(JW, IEQ), JW = 1, NW)
+
+    end DO
+
+    !Online output
+    if (output) then 
+       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
+    DO LL = 1, KLEV
+       d_u(:, LL) = RG * (RUW(:, LL + 1) - RUW(:, LL)) &
+            / (PH(:, LL + 1) - PH(:, LL)) * DTIME
+       d_v(:, LL) = RG * (RVW(:, LL + 1) - RVW(:, LL)) &
+            / (PH(:, LL + 1) - PH(:, LL)) * DTIME
+    ENDDO
+    ! ON CONSERVE LA MEMOIRE un certain temps AVEC UN SAVE
+        d_u = DTIME/DELTAT/REAL(NW) * d_u + (1.-DTIME/DELTAT) * d_u_sav
+        d_v = DTIME/DELTAT/REAL(NW) * d_v + (1.-DTIME/DELTAT) * d_v_sav
+	d_u_sav = d_u
+	d_v_sav = d_v
+
+    ! Cosmetic: evaluation of the cumulated stress
+
+    ZUSTR(:) = 0.
+    ZVSTR(:) = 0.
+    DO LL = 1, KLEV
+       ZUSTR(:) = ZUSTR(:) + D_U(:, LL) / RG * (PH(:, LL + 1) - PH(:, LL))
+       ZVSTR(:) = ZVSTR(:) + D_V(:, LL) / RG * (PH(:, LL + 1) - PH(:, LL))
+    ENDDO
+
+  END SUBROUTINE FLOTT_GWD_RAN
Index: trunk/LMDZ.VENUS/libf/phyvenus/readstartphy.F
===================================================================
--- trunk/LMDZ.VENUS/libf/phyvenus/readstartphy.F	(revision 175)
+++ trunk/LMDZ.VENUS/libf/phyvenus/readstartphy.F	(revision 778)
@@ -20,16 +20,19 @@
 #include "temps.h"
 c======================================================================
+      integer ngridmx
+      parameter (ngridmx=(2+(jjm-1)*iim - 1/jjm))
+
       CHARACTER*(*) fichnom
-      REAL rlat(klon), rlon(klon)
-      REAL tsol(klon)
-      REAL tsoil(klon,nsoilmx)
-      REAL albe(klon)
+      REAL rlat(ngridmx), rlon(ngridmx)
+      REAL tsol(ngridmx)
+      REAL tsoil(ngridmx,nsoilmx)
+      REAL albe(ngridmx)
 cIM BEG alblw
-      REAL alblw(klon)
+      REAL alblw(ngridmx)
 cIM END alblw
-      REAL radsol(klon)
-      REAL sollw(klon)
-      real solsw(klon)
-      real fder(klon)
+      REAL radsol(ngridmx)
+      REAL sollw(ngridmx)
+      real solsw(ngridmx)
+      real fder(ngridmx)
       INTEGER length
       PARAMETER (length=100)
@@ -114,7 +117,7 @@
          PRINT*, 'phyetat0: Le champ <TS> est present'
 #ifdef NC_DOUBLE
-         ierr = NF_GET_VAR_DOUBLE(nid, nvarid, tsol(1))
-#else
-         ierr = NF_GET_VAR_REAL(nid, nvarid, tsol(1))
+         ierr = NF_GET_VAR_DOUBLE(nid, nvarid, tsol)
+#else
+         ierr = NF_GET_VAR_REAL(nid, nvarid, tsol)
 #endif
          IF (ierr.NE.NF_NOERR) THEN
@@ -124,5 +127,5 @@
          xmin = 1.0E+20
          xmax = -1.0E+20
-         DO i = 1, klon
+         DO i = 1, ngridmx
             xmin = MIN(tsol(i),xmin)
             xmax = MAX(tsol(i),xmax)
@@ -143,5 +146,5 @@
          PRINT*, "phyetat0: Le champ <Tsoil"//str2//"> est absent"
          PRINT*, "          Il prend donc la valeur de surface"
-         DO i=1, klon
+         DO i=1, ngridmx
              tsoil(i,isoil)=tsol(i)
          ENDDO
@@ -158,4 +161,5 @@
       ENDIF
       ENDDO
+
 c
 c Lecture de albedo au sol:
@@ -179,5 +183,5 @@
          xmin = 1.0E+20
          xmax = -1.0E+20
-         DO i = 1, klon
+         DO i = 1, ngridmx
             xmin = MIN(albe(i),xmin)
             xmax = MAX(albe(i),xmax)
@@ -207,5 +211,5 @@
       xmin = 1.0E+20
       xmax = -1.0E+20
-      DO i = 1, klon
+      DO i = 1, ngridmx
          xmin = MIN(solsw(i),xmin)
          xmax = MAX(solsw(i),xmax)
@@ -233,5 +237,5 @@
       xmin = 1.0E+20
       xmax = -1.0E+20
-      DO i = 1, klon
+      DO i = 1, ngridmx
          xmin = MIN(sollw(i),xmin)
          xmax = MAX(sollw(i),xmax)
@@ -260,5 +264,5 @@
       xmin = 1.0E+20
       xmax = -1.0E+20
-      DO i = 1, klon
+      DO i = 1, ngridmx
          xmin = MIN(fder(i),xmin)
          xmax = MAX(fder(i),xmax)
@@ -285,5 +289,5 @@
       xmin = 1.0E+20
       xmax = -1.0E+20
-      DO i = 1, klon
+      DO i = 1, ngridmx
          xmin = MIN(radsol(i),xmin)
          xmax = MAX(radsol(i),xmax)
