Index: /trunk/LMDZ.MARS/README
===================================================================
--- /trunk/LMDZ.MARS/README	(revision 2148)
+++ /trunk/LMDZ.MARS/README	(revision 2149)
@@ -2720,2 +2720,6 @@
 - minor updates: enable output of Ls in diagfi files and make an error
   message more verbose.
+
+== 19/07/2019 == GG+EM
+- Add F.Lott's non-orographic GW parametrization. Disabled by default for now,
+  activated by setting calllott_nonoro=.true. in callphys.def
Index: /trunk/LMDZ.MARS/libf/phymars/callkeys.h
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/callkeys.h	(revision 2148)
+++ /trunk/LMDZ.MARS/libf/phymars/callkeys.h	(revision 2149)
@@ -15,5 +15,5 @@
      &   ,calltherm,callrichsl,callslope,tituscap,callyamada4,co2clouds &
      &   ,co2useh2o,meteo_flux,CLFvaryingCO2,spantCO2,CLFvarying        &
-     &   ,satindexco2,rdstorm
+     &   ,satindexco2,rdstorm,calllott_nonoro
      
       COMMON/callkeys_i/iradia,iaervar,iddist,ilwd,ilwb,ilwn,ncouche    &
@@ -25,5 +25,5 @@
      
       LOGICAL callrad,calldifv,calladj,callcond,callsoil,               &
-     &   season,diurnal,lwrite,calllott                                 &
+     &   season,diurnal,lwrite,calllott,calllott_nonoro                 &
      &   ,callstats,calleofdump                                         &
      &   ,callnirco2,callnlte,callthermos,callconduct,                  &
Index: /trunk/LMDZ.MARS/libf/phymars/conf_phys.F
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/conf_phys.F	(revision 2148)
+++ /trunk/LMDZ.MARS/libf/phymars/conf_phys.F	(revision 2149)
@@ -280,4 +280,10 @@
          write(*,*)" calllott = ",calllott
 
+         write(*,*)"call Lott's non-oro GWs parameterisation ",
+     &             "scheme ?"
+         calllott_nonoro=.false. ! default value
+         call getin("calllott_nonoro",calllott_nonoro)
+         write(*,*)" calllott_nonoro = ",calllott_nonoro
+
 ! rocket dust storm injection scheme
          write(*,*)"call rocket dust storm and slope lifting", 
Index: /trunk/LMDZ.MARS/libf/phymars/nonoro_gwd_ran_mod.F90
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/nonoro_gwd_ran_mod.F90	(revision 2149)
+++ /trunk/LMDZ.MARS/libf/phymars/nonoro_gwd_ran_mod.F90	(revision 2149)
@@ -0,0 +1,621 @@
+MODULE nonoro_gwd_ran_mod
+
+IMPLICIT NONE
+
+CONTAINS
+
+      SUBROUTINE NONORO_GWD_RAN(NLON,NLEV,DTIME, pp,  &
+                  zmax_therm, pt, pu, pv, pdt, pdu, pdv, &
+                  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 14, Gaussian distribution of the source
+    ! LMDz model online version      
+    ! ADAPTED FOR VENUS /  F. LOTT + S. LEBONNOIS
+    ! Version adapted on 03/04/2013:
+    !      - input flux compensated in the deepest layers
+    !                           
+    ! ADAPTED FOR MARS     G.GILLI     02/2016
+    !        Revision with F.Forget    06/2016  Variable EP-flux according to
+    !                                            PBL variation (max velocity thermals)
+    !---------------------------------------------------------------------------------
+
+      use dimphy, only: klon, klev
+!      use moyzon_mod, only: plevmoy
+      use comcstfi_h, only: g, pi, cpp, r
+!      use comgeomfi_h, only: latitude
+      implicit none
+
+      include "dimensions.h"
+      include "paramet.h"
+      include "yoegwd.h" 
+
+    ! 0. DECLARATIONS:
+
+    ! 0.1 INPUTS
+    INTEGER, intent(in):: NLON, NLEV 
+    REAL, intent(in):: DTIME ! Time step of the Physics
+    REAL, intent(in):: zmax_therm(NLON) ! altitude of max velocity thermals (m) 
+    REAL, intent(in):: pp(NLON,NLEV)   ! Pressure at full levels
+    REAL, intent(in):: pt(NLON,NLEV)   ! Temp at full levels 
+    REAL, intent(in):: pu(NLON,NLEV),pv(NLON,NLEV) ! Hor winds at full levels
+    REAL,INTENT(in) :: pdt(nlon,nlev) ! tendency on temperature
+    REAL,INTENT(in) :: pdu(nlon,nlev) ! tendency on zonal wind
+    REAL,INTENT(in) :: pdv(nlon,nlev) ! tendency on meridional wind
+
+    ! 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
+    REAL :: TT(NLON, NLEV)   ! Temp at full levels 
+
+    REAL :: UU(NLON, NLEV) , VV(NLON, NLEV)
+    ! Hor winds at full levels
+
+    INTEGER II, LL
+
+    ! 0.3.0 TIME SCALE OF THE LIFE CYCLE OF THE WAVES PARAMETERIZED
+
+    REAL DELTAT
+    
+    ! 0.3.1 GRAVITY-WAVES SPECIFICATIONS
+
+!!  NO= values of C, NW= total of GW
+    INTEGER, PARAMETER:: NK = 2, NP = 2, NO = 2, NW = NK * NP * NO
+    !Online output: change NO
+!    INTEGER, PARAMETER:: NK = 1, NP = 2, NO = 11, 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 EP-fluxes at the 2 semi levels surrounding the full level
+    REAL WWM(NW, KLON), WWP(NW, KLON)
+    ! Fluxes X and Y for each waves at 1/2 Levels
+    REAL RUWP(NW, KLON), RVWP(NW, KLON)
+    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
+	REAL ZMAXTH_TOP     ! Top of convective layer (approx.)
+
+   ! 0.3.2 PARAMETERS OF WAVES DISSIPATIONS
+
+    REAL SAT  ! saturation parameter
+    REAL RDISS, ZOISEC ! COEFF DE DISSIPATION, SECURITY FOR INTRINSIC FREQ
+
+    ! 0.3.3 BACKGROUND FLOW AT 1/2 LEVELS AND VERTICAL COORDINATE
+
+    REAL H0bis(KLON, KLEV) ! Characteristic Height of the atmosphere
+    REAL H0 ! Characteristic Height of the atmosphere
+    REAL PR, TR ! Reference Pressure and Temperature
+
+    REAL ZH(KLON, KLEV + 1) ! Log-pressure altitude (constant H0)
+    REAL ZHbis(KLON, KLEV + 1) ! Log-pressure altitude (varying H)
+
+    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,save :: output=.false.
+ ! 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
+    real,save,allocatable :: d_u_sav(:,:),d_v_sav(:,:)
+    LOGICAL,SAVE :: firstcall = .true.
+
+!      REAL ALEAS
+!      EXTERNAL ALEAS
+
+   !-----------------------------------------------------------------
+    !  1. INITIALISATIONS
+
+      IF (firstcall) THEN
+        write(*,*) "nonoro_gwd_ran: FLott non-oro GW scheme is active!"
+        allocate(d_u_sav(NLON,NLEV),d_v_sav(NLON,NLEV))
+        d_u_sav(:,:) = 0.
+        d_v_sav(:,:) = 0.
+        firstcall=.false.
+      ENDIF
+
+    ! Compute current values of temperature and winds
+    tt(:,:)=pt(:,:)+dtime*pdt(:,:)
+    uu(:,:)=pu(:,:)+dtime*pdu(:,:)
+    vv(:,:)=pv(:,:)+dtime*pdv(:,:)
+    
+    !    1.1 Basic parameter
+
+    !  PARAMETERS CORRESPONDING TO V3:
+     RUWMAX  = 1.e-4      ! best-case for MCS comparison
+!     RUWMAX  = 1.e-3  
+
+!     SAT    = 0.25       ! Saturation parameter: Sc in (12)        
+      SAT    = 1.         ! GG: This is the bestcase value for MCS comparison, 
+!                         ! but runs for 1 full MY in the thermosphere are instables.                 
+!!!   TESTS                   
+!      SAT    = 0.85
+! 
+!   RDISS  = 40         ! Diffusion parameter     
+    RDISS  =  1
+!    RDISS = 0.5
+
+    DELTAT=24.*3600.    ! Time scale of the waves (first introduced in 9b)
+!    DELTAT=DTIME
+
+    KMIN = 2.E-5        ! Min horizontal wavenumber
+    KMAX = 7.E-4        ! Max horizontal wavenumber
+
+    !Online output: one value only
+    if (output) then
+      KMIN = 6.3E-5
+      KMAX = 6.3E-5
+    endif
+    CMIN = 1.           ! Min phase velocity
+!    CMAX = 10.          ! Max phase speed velocity
+    CMAX = 30.        
+!    CMAX = 60.
+    XLAUNCH=0.4         ! Parameter that control launching altitude P/Ps
+!!! TEST: launching GW near the ground
+!    XLAUNCH =0.8
+
+! MARS above typical convective cell at 250 Pa
+    PR = 250.           ! Reference pressure  [Pa]   ! 
+    TR = 220.           ! Reference Temperature [K]  ! 
+
+    ZMAXTH_TOP  = 8000. ! Reference max altitude of convective layer [m]
+         
+    H0 = r * TR / g   ! Characteristic vertical scale height
+
+    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
+
+    !  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 (PH(NLON / 2, LL) / PH(NLON / 2, 1) > XLAUNCH) LAUNCH = LL
+!       IF (plevmoy(LL) / plevmoy(1) > XLAUNCH) LAUNCH = LL   ! on venus
+    ENDDO
+ 
+    if (output) print*, " WE ARE IN FLOTT GW SCHEME "
+    
+    !print*, "launch=",LAUNCH
+    !print*,"launch p =",PH(NLON/2, LAUNCH) 
+
+    ! Log pressure vert. coordinate
+    DO LL = 1, KLEV + 1 
+       ZH(:, LL) = H0 * LOG(PR / (PH(:, LL) + PSEC))
+     !  print*,'ZH at each level',LL, ZH(NLON/2, LL)
+    end DO
+
+    if (output) then
+    ! altitude above surface
+       ZHbis(:,1) = 0.    
+       DO LL = 2, KLEV + 1 
+          H0bis(:, LL-1) = r * TT(:, LL-1) / g 
+          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
+       ! GG test	
+       !print*, 'TT, UH, VH, ZH at launch', TT(NLON/2,LAUNCH), UH(NLON/2,LAUNCH),VH(NLON/2, LAUNCH), ZH(NLON/2,LAUNCH)	
+       ! BVSEC: BV Frequency
+       BV(:, LL) = 0.5 * (TT(:, LL) + TT(:, LL - 1)) &
+            * r**2 / cpp / H0**2 + (TT(:, LL) &
+            - TT(:, LL - 1)) / (ZH(:, LL) - ZH(:, LL - 1)) * r / H0
+       BV(:,LL) =SQRT(MAX(BVSEC,BV(:,LL)))
+    end DO
+       !GG test
+       !print*, 'BV freq in flott_gwnoro:',LAUNCH,  BV(NLON/2, LAUNCH)  
+
+    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
+    !-------------------------------------------
+
+    ! The mod function of here a weird arguments
+    ! 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
+    DO JP = 1, NP
+       DO JK = 1, NK
+          DO JO = 1, NO
+             JW = JW + 1
+             !  Angle
+             ZP(JW) = 2. * pi * REAL(JP - 1) / REAL(NP) 
+             DO II = 1, KLON
+                ! Horizontal wavenumber amplitude
+!                ZK(JW, II) = KMIN + (KMAX - KMIN) * MOD(TT(II, JW) * 100., 1.)
+                ZK(JW, II) = KMIN + (KMAX - KMIN) * ALEAS(0.)
+                ! Horizontal phase speed
+!                CPHA = CMIN + (CMAX - CMIN) * MOD(TT(II, JW)**2, 1.)
+                CPHA = CMIN + (CMAX - CMIN) * ALEAS(0.)
+       !Online output: linear
+                if (output) CPHA = CMIN + (CMAX - CMIN) * (JO-1)/(NO-1)
+                ! 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 
+!!
+!!! CHANGE on MARS GCM: variable EP-FLUX according to:
+!     TEST1:   PBL diurnal variation
+!!
+!                 RUW0(JW, II) = MAX(RUWMAX / 100.,  &
+!                               zmax_therm(II)/ZMAXTH_TOP * RUWMAX) 
+          
+!!! BELOW: previous version from Flott original version 
+                RUW0(JW, II) = RUWMAX &
+                     * ALEAS(0.)     ! &
+!!! The factor cos**8:Limit the fluxes to tropical regions
+!                    * COS(pi/180.*latitude(II))**8
+       !Online output: fixed to max
+                if (output) RUW0(JW, II) = RUWMAX
+             ENDDO
+          end DO
+       end DO
+    end DO
+
+    ! 4. COMPUTE THE FLUXES
+    !--------------------------
+
+    !  4.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
+
+    !!!  EARTH VERSION  !!!!
+    !  4.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
+
+    !!!  VENUS VERSION  !!!!
+    !  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
+    !      the breaking.
+
+    !Online output
+    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
+
+
+       !  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)
+               * ABS(ZOP(JW, :))**3 /BV(:, LL+1) & 
+               * EXP(-ZH(:, LL + 1)/H0) * SAT**2*KMIN**2/ZK(JW, :)**4)  
+       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)) *  WWP(JW, :)**2
+          RVWP(JW, :) = ZOP(JW, :)/MAX(ABS(ZOP(JW, :)), ZOISEC)**2 &
+               *BV(:,LL+1)&
+               * SIN(ZP(JW)) *  WWP(JW, :)**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) then
+         do JW=1,NW
+            if(RUWP(JW, IEQ).gt.0.) then
+              RUWP(JW, IEQ) = max(RUWP(JW, IEQ), 1.e-99)
+            else
+              RUWP(JW, IEQ) = min(RUWP(JW, IEQ), -1.e-99)
+            endif
+         enddo
+                   WRITE(11,outform) ZH(IEQ, LL+1) / 1000., &
+                                  ZHbis(IEQ, LL+1) / 1000., &
+                                  (RUWP(JW, IEQ), JW = 1, NW)
+       endif
+
+    end DO
+
+    ! 5 CALCUL DES TENDANCES:
+    !------------------------
+
+    ! 5.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) = G * (RUW(:, LL + 1) - RUW(:, LL)) &
+            / (PH(:, LL + 1) - PH(:, LL)) * DTIME
+       d_v(:, LL) = G * (RVW(:, LL + 1) - RVW(:, LL)) &
+            / (PH(:, LL + 1) - PH(:, LL)) * DTIME
+    ENDDO
+    d_t(:,:) = 0.
+    ! 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) / g * (PH(:, LL + 1) - PH(:, LL))
+       ZVSTR(:) = ZVSTR(:) + D_V(:, LL) / g * (PH(:, LL + 1) - PH(:, LL))
+    ENDDO
+
+  END SUBROUTINE NONORO_GWD_RAN
+
+!===================================================================
+!===================================================================
+!===================================================================
+!===================================================================
+
+  REAL FUNCTION ALEAS (R)
+  IMPLICIT NONE
+!***BEGIN PROLOGUE  ALEAS
+!***PURPOSE  Generate a uniformly distributed random number.
+!***LIBRARY   SLATEC (FNLIB)
+!***CATEGORY  L6A21
+!***TYPE      SINGLE PRECISION (ALEAS-S)
+!***KEYWORDS  FNLIB, ALEAS NUMBER, SPECIAL FUNCTIONS, UNIFORM
+!***AUTHOR  Fullerton, W., (LANL)
+!***DESCRIPTION
+!
+!      This pseudo-random number generator is portable among a wide
+! variety of computers.  RAND(R) undoubtedly is not as good as many
+! readily available installation dependent versions, and so this
+! routine is not recommended for widespread usage.  Its redeeming
+! feature is that the exact same random numbers (to within final round-
+! off error) can be generated from machine to machine.  Thus, programs
+! that make use of random numbers can be easily transported to and
+! checked in a new environment.
+!
+!      The random numbers are generated by the linear congruential
+! method described, e.g., by Knuth in Seminumerical Methods (p.9),
+! Addison-Wesley, 1969.  Given the I-th number of a pseudo-random
+! sequence, the I+1 -st number is generated from
+!             X(I+1) = (A*X(I) + C) MOD M,
+! where here M = 2**22 = 4194304, C = 1731 and several suitable values
+! of the multiplier A are discussed below.  Both the multiplier A and
+! random number X are represented in double precision as two 11-bit
+! words.  The constants are chosen so that the period is the maximum
+! possible, 4194304.
+!
+!      In order that the same numbers be generated from machine to
+! machine, it is necessary that 23-bit integers be reducible modulo
+! 2**11 exactly, that 23-bit integers be added exactly, and that 11-bit
+! integers be multiplied exactly.  Furthermore, if the restart option
+! is used (where R is between 0 and 1), then the product R*2**22 =
+! R*4194304 must be correct to the nearest integer.
+!
+!      The first four random numbers should be .0004127026,
+! .6750836372, .1614754200, and .9086198807.  The tenth random number
+! is .5527787209, and the hundredth is .3600893021 .  The thousandth
+! number should be .2176990509 .
+!
+!      In order to generate several effectively independent sequences
+! with the same generator, it is necessary to know the random number
+! for several widely spaced calls.  The I-th random number times 2**22,
+! where I=K*P/8 and P is the period of the sequence (P = 2**22), is
+! still of the form L*P/8.  In particular we find the I-th random
+! number multiplied by 2**22 is given by
+! 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
+! RAND=  0  5*P/8  2*P/8  7*P/8  4*P/8  1*P/8  6*P/8  3*P/8  0
+! Thus the 4*P/8 = 2097152 random number is 2097152/2**22.
+!
+!      Several multipliers have been subjected to the spectral test
+! (see Knuth, p. 82).  Four suitable multipliers roughly in order of
+! goodness according to the spectral test are
+!    3146757 = 1536*2048 + 1029 = 2**21 + 2**20 + 2**10 + 5
+!    2098181 = 1024*2048 + 1029 = 2**21 + 2**10 + 5
+!    3146245 = 1536*2048 +  517 = 2**21 + 2**20 + 2**9 + 5
+!    2776669 = 1355*2048 + 1629 = 5**9 + 7**7 + 1
+!
+!      In the table below LOG10(NU(I)) gives roughly the number of
+! random decimal digits in the random numbers considered I at a time.
+! C is the primary measure of goodness.  In both cases bigger is better.
+!
+!                   LOG10 NU(I)              C(I)
+!       A       I=2  I=3  I=4  I=5    I=2  I=3  I=4  I=5
+!
+!    3146757    3.3  2.0  1.6  1.3    3.1  1.3  4.6  2.6
+!    2098181    3.3  2.0  1.6  1.2    3.2  1.3  4.6  1.7
+!    3146245    3.3  2.2  1.5  1.1    3.2  4.2  1.1  0.4
+!    2776669    3.3  2.1  1.6  1.3    2.5  2.0  1.9  2.6
+!   Best
+!    Possible   3.3  2.3  1.7  1.4    3.6  5.9  9.7  14.9
+!
+!             Input Argument --
+! R      If R=0., the next random number of the sequence is generated.
+!        If R .LT. 0., the last generated number will be returned for
+!          possible use in a restart procedure.
+!        If R .GT. 0., the sequence of random numbers will start with
+!          the seed R mod 1.  This seed is also returned as the value of
+!          RAND provided the arithmetic is done exactly.
+!
+!             Output Value --
+! RAND   a pseudo-random number between 0. and 1.
+!
+!***REFERENCES  (NONE)
+!***ROUTINES CALLED  (NONE)
+!***REVISION HISTORY  (YYMMDD)
+!   770401  DATE WRITTEN
+!   890531  Changed all specific intrinsics to generic.  (WRB)
+!   890531  REVISION DATE from Version 3.2
+!   891214  Prologue converted to Version 4.0 format.  (BAB)
+!***END PROLOGUE  RAND
+      INTEGER IA1, IA0, IA1MA0, IC, IX1, IX0
+      SAVE IA1, IA0, IA1MA0, IC, IX1, IX0
+      DATA IA1, IA0, IA1MA0 /1536, 1029, 507/
+      DATA IC /1731/
+      DATA IX1, IX0 /0, 0/
+      REAL R
+      INTEGER IY0,IY1
+!***FIRST EXECUTABLE STATEMENT  RAND
+!
+!           A*X = 2**22*IA1*IX1 + 2**11*(IA1*IX1 + (IA1-IA0)*(IX0-IX1)
+!                   + IA0*IX0) + IA0*IX0
+!
+      IF (R.EQ.0.) THEN
+       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)
+      ENDIF
+
+      IF (R.GT.0.) THEN
+       IX1 = MOD(R,1.)*4194304. + 0.5
+       IX0 = MOD (IX1, 2048)
+       IX1 = (IX1-IX0)/2048
+      ENDIF
+
+      ALEAS = IX1*2048 + IX0
+      ALEAS = ALEAS / 4194304.
+      RETURN
+
+      END FUNCTION ALEAS
+
+END MODULE nonoro_gwd_ran_mod
Index: /trunk/LMDZ.MARS/libf/phymars/physiq_mod.F
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/physiq_mod.F	(revision 2148)
+++ /trunk/LMDZ.MARS/libf/phymars/physiq_mod.F	(revision 2149)
@@ -62,4 +62,5 @@
       use iono_h, only: allocate_param_iono
       use compute_dtau_mod, only: compute_dtau
+      use nonoro_gwd_ran_mod, only: nonoro_gwd_ran
 #ifdef MESOSCALE
       use comsoil_h, only: mlayer,layer
@@ -414,4 +415,14 @@
       real*8  varerr
 
+C  Non-oro GW drag & Calcul of Brunt-Vaisala freq. (BV2) 
+      REAL ztetalev(ngrid,nlayer)
+      real zdtetalev(ngrid,nlayer), zdzlev(ngrid,nlayer)
+      REAL bv2(ngrid,nlayer)    ! BV2 at zlev   
+c  Non-oro GW tendencies
+      REAL d_u_hin(ngrid,nlayer), d_v_hin(ngrid,nlayer)
+      REAL d_t_hin(ngrid,nlayer)
+c  Diagnostics 2D of gw_nonoro
+      REAL zustrhi(ngrid), zvstrhi(ngrid)
+
 c Variables for PBL
       REAL zz1(ngrid)
@@ -1270,12 +1281,37 @@
       ENDIF ! of IF(calladj)
 
+c-----------------------------------------------------
+c    8. Non orographic Gravity waves :
+c -------------------------------------------------
+
+      IF (calllott_nonoro) THEN
+
+         CALL nonoro_gwd_ran(ngrid,nlayer,ptimestep,zplay,
+     &               zmax_th,                      ! max altitude reached by thermals (m)
+     &               pt, pu, pv,
+     &               pdt, pdu, pdv,
+     &               zustrhi,zvstrhi,
+     &               d_t_hin, d_u_hin, d_v_hin)
+
+!  Update tendencies
+         pdt(1:ngrid,1:nlayer)=pdt(1:ngrid,1:nlayer)
+     &                         +d_t_hin(1:ngrid,1:nlayer)
+!        d_t_hin(:,:)= d_t_hin(:,:)/ptimestep ! K/s (for outputs?)
+         pdu(1:ngrid,1:nlayer)=pdu(1:ngrid,1:nlayer)
+     &                         +d_u_hin(1:ngrid,1:nlayer) 
+!        d_u_hin(:,:)= d_u_hin(:,:)/ptimestep ! (m/s)/s (for outputs?)
+         pdv(1:ngrid,1:nlayer)=pdv(1:ngrid,1:nlayer)
+     &                         +d_v_hin(1:ngrid,1:nlayer)
+!        d_v_hin(:,:)= d_v_hin(:,:)/ptimestep ! (m/s)/s (for outputs?)
+
+      ENDIF ! of IF (calllott_nonoro)
 
 c-----------------------------------------------------------------------
-c   8. Specific parameterizations for tracers 
+c   9. Specific parameterizations for tracers 
 c:   -----------------------------------------
 
       if (tracer) then 
 
-c   8a. Water and ice
+c   9a. Water and ice
 c     ---------------
 
@@ -1357,5 +1393,5 @@
          END IF  ! of IF (water)
 
-c   8a bis. CO2 clouds (CL & JA)
+c   9a bis. CO2 clouds (CL & JA)
 c        ---------------------------------------
 c        CO2 ice cloud condensation in the atmosphere
@@ -1479,5 +1515,5 @@
       END IF                    ! of IF (co2clouds)
 
-c   8b. Aerosol particles
+c   9b. Aerosol particles
 c     -------------------
 c        ----------
@@ -1571,5 +1607,5 @@
       ENDIF
 c
-c   8c. Chemical species
+c   9c. Chemical species
 c     ------------------
 
@@ -1632,5 +1668,5 @@
 #endif
 
-c   8d. Updates
+c   9d. Updates
 c     ---------
 
@@ -1650,5 +1686,5 @@
 #ifndef MESOSCALE
 c-----------------------------------------------------------------------
-c   9. THERMOSPHERE CALCULATION
+c   10. THERMOSPHERE CALCULATION
 c-----------------------------------------------------------------------
 
@@ -1675,5 +1711,5 @@
 #endif
 c-----------------------------------------------------------------------
-c   10. Carbon dioxide condensation-sublimation:
+c   11. Carbon dioxide condensation-sublimation:
 c     (should be the last atmospherical physical process to be computed)
 c   -------------------------------------------
@@ -1748,9 +1784,9 @@
 
 c-----------------------------------------------------------------------
-c   11. Surface  and sub-surface soil temperature
+c   12. Surface  and sub-surface soil temperature
 c-----------------------------------------------------------------------
 c
 c
-c   11.1 Increment Surface temperature:
+c   12.1 Increment Surface temperature:
 c   ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
@@ -1793,5 +1829,5 @@
 
 c
-c   11.2 Compute soil temperatures and subsurface heat flux:
+c   12.2 Compute soil temperatures and subsurface heat flux:
 c   ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
       IF (callsoil) THEN
@@ -1841,5 +1877,5 @@
 
 c-----------------------------------------------------------------------
-c  12. Write output files
+c  13. Write output files
 c  ----------------------
 
@@ -2539,5 +2575,5 @@
          call WRITEDIAGFI(ngrid,"v","Meridional wind","m.s-1",3,zv)
          call WRITEDIAGFI(ngrid,"w","Vertical wind","m.s-1",3,pw)
-         call WRITEDIAGFI(ngrid,"rho","density","none",3,rho)
+         call WRITEDIAGFI(ngrid,"rho","density","kg.m-3",3,rho)
 c        call WRITEDIAGFI(ngrid,"q2","q2","kg.m-3",3,q2)
 c        call WRITEDIAGFI(ngrid,'Teta','T potentielle','K',3,zh)
@@ -2878,4 +2914,14 @@
          end if  ! (tracer.and.(dustbin.ne.0))
 
+c        ----------------------------------------------------------
+c        GW non-oro outputs
+c        ----------------------------------------------------------
+
+         if(calllott_nonoro) then 
+           call WRITEDIAGFI(ngrid,"dugwno","GW non-oro dU","m/s2",
+     $           3,d_u_hin/ptimestep)
+           call WRITEDIAGFI(ngrid,"dvgwno","GW non-oro dV","m/s2",
+     $           3,d_v_hin/ptimestep)
+         endif                  !(calllott_nonoro)
 
 c        ----------------------------------------------------------
Index: /trunk/LMDZ.MARS/libf/phymars/yoegwd.h
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/yoegwd.h	(revision 2148)
+++ /trunk/LMDZ.MARS/libf/phymars/yoegwd.h	(revision 2149)
@@ -1,16 +1,12 @@
-C     -----------------------------------------------------------------
-C*    *COMMON* *YOEGWD* - PARAMETERS FOR GRAVITY WAVE DRAG CALCULATIONS
-C     -----------------------------------------------------------------
-C
-      real GFRCRIT,GRCRIT,GVCRIT,GKDRAG,GKDRAGL,GHMAX
-     *             , GRAHILO,GSIGCR,GSSEC, GTSEC,  GVSEC
-     *             ,GKWAKE
-
-      integer NKTOPG
-
-      COMMON/YOEGWD/ GFRCRIT,GRCRIT,GVCRIT,GKDRAG,GKDRAGL,GHMAX
-     *             , GRAHILO,GSIGCR,NKTOPG,GSSEC, GTSEC,  GVSEC
-     *             ,GKWAKE
-C
+!     -----------------------------------------------------------------
+!*    *COMMON* *YOEGWD* - PARAMETERS FOR GRAVITY WAVE DRAG CALCULATIONS
+!     -----------------------------------------------------------------
+!
+      real    ::  GFRCRIT,GRCRIT,GVCRIT,GKDRAG,GKDRAGL,GHMAX
+      real    ::   GRAHILO,GSIGCR,GSSEC, GTSEC,  GVSEC, GKWAKE
+      integer ::  NKTOPG
+      COMMON/YOEGWD/ GFRCRIT,GRCRIT,GVCRIT,GKDRAG,GKDRAGL,GHMAX         &
+     &             ,GRAHILO,GSIGCR,NKTOPG,GSSEC, GTSEC, GVSEC,GKWAKE
 
 
+
