Index: /trunk/LMDZ.MARS/README
===================================================================
--- /trunk/LMDZ.MARS/README	(revision 2219)
+++ /trunk/LMDZ.MARS/README	(revision 2220)
@@ -2807,4 +2807,10 @@
 Add the observer simulator program simu_MCS_temp.F90 to the LMDZ.MARS/util/ directory, as well as a simu_MCS_temp.def. This program reads a MCS file (binned by Luca Montabone), extracts the daytime and nighttime temperatures from a GCM simulation file, and outputs a netcdf file in the same format as the MCS file.
 More details in simu_MCS_temp.def and in the preface of simu_MCS_temp.F90.
+<<<<<<< .mine
+
+== 22/01/2020 == DB (the return!)
+Update the nonorographic gravity waves drag parametrization (from the r3599 of Earth s model LMDZ6): 1) add east_ and west_gwdstress variables, 2) delete aleas function: random waves are producted by the MOD function, 3) reproductibility concerning the number of procs is now validated, 4) changing name of some variables like RUW, RVW, ZOP, ZOM,..., 5) tendency of winds due to GW drag are now module variables and written in restart files. 
+
+=======
 
 == 07/01/2020 == AD
@@ -2829,2 +2835,3 @@
 Change "latentheat" flag to a more descriptive "latentheat_surfwater" and
 set its default value to .true.
+>>>>>>> .r2219
Index: /trunk/LMDZ.MARS/libf/phymars/callkeys.h
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/callkeys.h	(revision 2219)
+++ /trunk/LMDZ.MARS/libf/phymars/callkeys.h	(revision 2220)
@@ -15,6 +15,11 @@
      &   ,calltherm,callrichsl,callslope,tituscap,callyamada4,co2clouds &
      &   ,co2useh2o,meteo_flux,CLFvaryingCO2,spantCO2,CLFvarying        &
+<<<<<<< .mine
+     &   ,satindexco2,rdstorm,slpwind,calllott_nonoro,latentheat        &
+     &   ,gwd_convective_source
+=======
      &   ,satindexco2,rdstorm,slpwind,calllott_nonoro                   &
      &   ,latentheat_surfwater
+>>>>>>> .r2219
      
       COMMON/callkeys_i/iradia,iaervar,iddist,ilwd,ilwb,ilwn,ncouche    &
@@ -37,4 +42,5 @@
       logical callg2d
       logical linear
+      logical gwd_convective_source
 
       real topdustref
Index: /trunk/LMDZ.MARS/libf/phymars/nonoro_gwd_ran_mod.F90
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/nonoro_gwd_ran_mod.F90	(revision 2219)
+++ /trunk/LMDZ.MARS/libf/phymars/nonoro_gwd_ran_mod.F90	(revision 2220)
@@ -3,9 +3,13 @@
 IMPLICIT NONE
 
+REAL,allocatable,save :: du_nonoro_gwd(:,:) ! Zonal wind tendency due to GWD 
+REAL,allocatable,save :: dv_nonoro_gwd(:,:) ! Meridional wind tendency due to GWD
+
 CONTAINS
 
-      SUBROUTINE NONORO_GWD_RAN(NLON,NLEV,DTIME, pp,  &
+      SUBROUTINE NONORO_GWD_RAN(ngrid,nlayer,DTIME, pp,  &
                   zmax_therm, pt, pu, pv, pdt, pdu, pdv, &
-                  zustr,zvstr,d_t, d_u, d_v)
+                  zustr,zvstr,d_t, d_u, d_v, &
+                  east_gwstress, west_gwstress)
 
     !--------------------------------------------------------------------------------
@@ -21,12 +25,17 @@
     ! ADAPTED FOR MARS     G.GILLI     02/2016
     !        Revision with F.Forget    06/2016  Variable EP-flux according to
-    !                                            PBL variation (max velocity thermals)
+    !                                           PBL variation (max velocity thermals)
+    ! UPDATED              D.BARDET    01/2020  - reproductibility of the
+    !                                           launching altitude calculation 
+    !                                           - wave characteristic
+    !                                           calculation using MOD
+    !                                           - adding east_gwstress and
+    !                                           west_gwstress variables    
     !---------------------------------------------------------------------------------
 
-      use dimphy, only: klon, klev
-!      use moyzon_mod, only: plevmoy
       use comcstfi_h, only: g, pi, cpp, r
-!      use comgeomfi_h, only: latitude
       USE ioipsl_getin_p_mod, ONLY : getin_p
+      use assert_m, only : assert
+      use vertical_layers_mod, only : presnivs
       implicit none
 
@@ -34,87 +43,98 @@
       include "paramet.h"
       include "yoegwd.h" 
+      include "callkeys.h"
+
+      CHARACTER (LEN=20) :: modname='flott_gwd_rando'
+      CHARACTER (LEN=80) :: abort_message
+
 
     ! 0. DECLARATIONS:
 
     ! 0.1 INPUTS
-    INTEGER, intent(in):: NLON, NLEV 
+    INTEGER, intent(in):: ngrid, nlayer 
     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
+    REAL, intent(in):: zmax_therm(ngrid) ! altitude of max velocity thermals (m) 
+    REAL, intent(in):: pp(ngrid,nlayer)   ! Pressure at full levels
+    REAL, intent(in):: pt(ngrid,nlayer)   ! Temp at full levels 
+    REAL, intent(in):: pu(ngrid,nlayer),pv(ngrid,nlayer) ! Hor winds at full levels
+    REAL,INTENT(in) :: pdt(ngrid,nlayer) ! tendency on temperature
+    REAL,INTENT(in) :: pdu(ngrid,nlayer) ! tendency on zonal wind
+    REAL,INTENT(in) :: pdv(ngrid,nlayer) ! 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
+    REAL, intent(out):: zustr(ngrid), zvstr(ngrid) ! Surface Stresses
+    REAL, intent(inout):: d_t(ngrid, nlayer)        ! Tendency on Temp.
+    REAL, intent(inout):: d_u(ngrid, nlayer), d_v(ngrid, nlayer) ! tendency on winds
+    REAL, intent(inout):: east_gwstress(ngrid, nlayer) ! Profile of eastward stress
+    REAL, intent(inout):: west_gwstress(ngrid, nlayer) ! Profile of westward stress
 
     ! 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
+    REAL :: TT(ngrid, nlayer)   ! Temp at full levels 
+    REAL :: UU(ngrid, nlayer) , VV(ngrid, nlayer) ! Hor winds at full levels
+    REAL :: BVLOW(ngrid)
+    REAL :: DZ
+
+    INTEGER II, JJ, LL
 
     ! 0.3.0 TIME SCALE OF THE LIFE CYCLE OF THE WAVES PARAMETERIZED
 
-    REAL DELTAT
+    REAL, parameter:: DELTAT = 24. * 3600.
     
     ! 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, PARAMETER:: NK = 2 ! number of horizontal wavenumbers
+    INTEGER, PARAMETER:: NP = 2 ! directions (eastward and westward) phase speed
+    INTEGER, PARAMETER:: NO = 2 ! absolute values of phase speed
+    INTEGER, PARAMETER:: NW = NK * NP * NO ! Total numbers of gravity waves
     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,SAVE :: RUWMAX         ! Max value
+    INTEGER, PARAMETER:: NA = 5 ! number of realizations to get the phase speed
+    REAL, parameter:: kmax = 7.e-4 ! Max horizontal wavenumber
+    REAL, parameter:: kmin = 2.e-5 ! Min horizontal wavenumber
+    REAL, parameter:: cmax = 30.   ! Max horizontal absolute phase velocity
+    REAL, parameter:: cmin = 1.    ! Min horizontal absolute phase velocity    
+    REAL CPHA                   ! absolute PHASE VELOCITY frequency
+    REAL ZK(NW, ngrid)           ! Horizontal wavenumber amplitude
+    REAL ZP(NW, ngrid)           ! Horizontal wavenumber angle        
+    REAL ZO(NW, ngrid)           ! Absolute frequency
+
+    REAL intr_freq_m(nw, ngrid)          ! Waves Intr. freq. at the 1/2 lev below the full level (previous name: ZOM)
+    REAL intr_freq_p(nw, ngrid)          ! Waves Intr. freq. at the 1/2 lev above the full level (previous name: ZOP)
+    REAL wwm(nw, ngrid)                  ! Wave EP-fluxes at the 1/2 level below the full level
+    REAL wwp(nw, ngrid)                  ! Wave EP-fluxes at the 1/2 level above the full level
+    REAL u_epflux_p(nw, ngrid)           ! Partial zonal flux (=for each wave) at the 1/2 level above the full level (previous name: RUWP)
+    REAL v_epflux_p(nw, ngrid)           ! Partial meridional flux (=for each wave) at the 1/2 level above the full level (previous name: RVWP)
+    REAL u_epflux_tot(ngrid, nlayer + 1) ! Total zonal flux (=for all waves (nw)) at the 1/2 level above the full level (3D) (previous name: RUW)  
+    REAL v_epflux_tot(ngrid, nlayer + 1) ! Total meridional flux (=for all waves (nw)) at the 1/2 level above the full level (3D) (previous name: RVW) 
+    REAL epflux_0(nw, ngrid)             ! Fluxes at launching level (previous name: RUW0)
+    REAL, save :: epflux_max             ! Max EP flux value at launching altitude (previous name: RUWMAX)
     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
+    REAL, parameter:: xlaunch = 0.4      ! Control the launching altitude
+    REAL, parameter:: zmaxth_top = 8000. ! Top of convective layer (approx.)
+
+
+    REAL PREC(ngrid)
+    REAL PRMAX ! Maximum value of PREC, and for which our linear formula
+
+
+    ! 0.3.2 PARAMETERS OF WAVES DISSIPATIONS
+    REAL, parameter:: sat   = 1.     ! saturation parameter
+    REAL, parameter:: rdiss = 1.     ! coefficient of dissipation
+    REAL, parameter:: zoisec = 1.e-6 ! security for intrinsic freguency
+
+    ! 0.3.3 Background flow at 1/2 levels and vertical coordinate
+    REAL H0bis(ngrid, nlayer)          ! Characteristic Height of the atmosphere
+    REAL, save:: H0                 ! Characteristic Height of the atmosphere
+    REAL, parameter:: pr = 250      ! Reference pressure [Pa]
+    REAL, parameter:: tr = 220.     ! Reference temperature [K]
+    REAL ZH(ngrid, nlayer + 1)         ! Log-pressure altitude (constant H0)
+    REAL ZHbis(ngrid, nlayer + 1)      ! Log-pressure altitude (varying H)
+    REAL UH(ngrid, nlayer + 1)         ! zonal wind at 1/2 levels
+    REAL VH(ngrid, nlayer + 1)         ! meridional wind at 1/2 levels
+    REAL PH(ngrid, nlayer + 1)         ! Pressure at 1/2 levels
+    REAL, parameter:: psec = 1.e-6  ! Security to avoid division by 0 pressure
+    REAL BV(ngrid, nlayer + 1)         ! Brunt Vaisala freq. (BVF) at 1/2 levels
+    REAL, parameter:: bvsec = 1.e-5 ! Security to avoid negative BV  
+    REAL HREF(nlayer + 1)             ! Reference altitude for launching alt.
+
 
 ! COSMETICS TO DIAGNOSE EACH WAVES CONTRIBUTION.
@@ -125,24 +145,31 @@
     integer      ieq
 
-! ON CONSERVE LA MEMOIRE un certain temps AVEC UN SAVE
-    real,save,allocatable :: d_u_sav(:,:),d_v_sav(:,:)
+    REAL RAN_NUM_1,RAN_NUM_2,RAN_NUM_3
+
+
     LOGICAL,SAVE :: firstcall = .true.
 
-!      REAL ALEAS
-!      EXTERNAL ALEAS
 
    !-----------------------------------------------------------------
     !  1. INITIALISATIONS
 
-      IF (firstcall) THEN
+     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.
-        RUWMAX=1.E-4
-        CALL getin_p("nonoro_gwd_RUWMAX",RUWMAX)
-        write(*,*) "nonoro_gwd_ran: RUWMAX=",RUWMAX
-        firstcall=.false.
-      ENDIF
+        epflux_max = 7.E-7 ! Mars' value !!
+        call getin_p("nonoro_gwd_epflux_max", epflux_max)
+        write(*,*) "nonoro_gwd_ran: epflux_max=", epflux_max
+        ! Characteristic vertical scale height
+        H0 = r * tr / g
+        ! Control
+        if (deltat .LT. dtime) THEN
+             call abort_physic("nonoro_gwd_ran","gwd random: deltat lower than dtime!",1)
+        endif
+        if (nlayer .LT. nw) THEN
+             call abort_physic("nonoro_gwd_ran","gwd random: nlayer lower than nw!",1)
+        endif
+        firstcall = .false.
+     ENDIF
+
+    gwd_convective_source=.false.
 
     ! Compute current values of temperature and winds
@@ -150,64 +177,6 @@
     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 = 5.e-6
-!     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
@@ -218,27 +187,32 @@
 
     ! Pressure and Inv of pressure, Temperature / at 1/2 level
-    DO LL = 2, KLEV
+    DO LL = 2, nlayer
        PH(:, LL) = EXP((LOG(PP(:, LL)) + LOG(PP(:, LL - 1))) / 2.)
     end DO
 
-    PH(:, KLEV + 1) = 0. 
+    PH(:, nlayer + 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
+    !Pour revenir a la version non reproductible en changeant le nombre de
+    !process
+    ! Reprend la formule qui calcule PH en fonction de PP=play
+    DO LL = 2, nlayer
+       HREF(LL) = EXP((LOG(presnivs(LL))+ LOG(presnivs(LL - 1))) / 2.)
+    end DO
+    HREF(nlayer + 1) = 0.
+    HREF(1) = 2. * presnivs(1) - HREF(2)
+
+    LAUNCH=0
+    DO LL = 1, nlayer
+       IF (HREF(LL) / HREF(1) > XLAUNCH) LAUNCH = LL
     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 
+    DO LL = 1, nlayer + 1 
        ZH(:, LL) = H0 * LOG(PR / (PH(:, LL) + PSEC))
-     !  print*,'ZH at each level',LL, ZH(NLON/2, LL)
     end DO
 
@@ -246,5 +220,5 @@
     ! altitude above surface
        ZHbis(:,1) = 0.    
-       DO LL = 2, KLEV + 1 
+       DO LL = 2, nlayer + 1 
           H0bis(:, LL-1) = r * TT(:, LL-1) / g 
           ZHbis(:, LL) = ZHbis(:, LL-1) &
@@ -254,9 +228,9 @@
 
     ! Winds and BV frequency
-    DO LL = 2, KLEV
+    DO LL = 2, nlayer
        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)	
+       !print*, 'TT, UH, VH, ZH at launch', TT(ngrid/2,LAUNCH), UH(ngrid/2,LAUNCH),VH(ngrid/2, LAUNCH), ZH(ngrid/2,LAUNCH)	
        ! BVSEC: BV Frequency
        BV(:, LL) = 0.5 * (TT(:, LL) + TT(:, LL - 1)) &
@@ -266,12 +240,12 @@
     end DO
        !GG test
-       !print*, 'BV freq in flott_gwnoro:',LAUNCH,  BV(NLON/2, LAUNCH)  
+       !print*, 'BV freq in flott_gwnoro:',LAUNCH,  BV(ngrid/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)
+    BV(:, nlayer + 1) = BV(:, nlayer)
+    UH(:, nlayer + 1) = UU(:, nlayer)
+    VH(:, nlayer + 1) = VV(:, nlayer)
 
 
@@ -283,22 +257,25 @@
     ! 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
+    DO JW = 1, NW
              !  Angle
-             ZP(JW) = 2. * pi * REAL(JP - 1) / REAL(NP) 
-             DO II = 1, KLON
+             DO II = 1, ngrid
+                ! Angle (0 or PI so far)
+                RAN_NUM_1=MOD(TT(II, JW) * 10., 1.)
+                RAN_NUM_2= MOD(TT(II, JW) * 100., 1.)
+                ZP(JW, II) = (SIGN(1., 0.5 - RAN_NUM_1) + 1.) &
+                     * PI / 2.
                 ! Horizontal wavenumber amplitude
-!                ZK(JW, II) = KMIN + (KMAX - KMIN) * MOD(TT(II, JW) * 100., 1.)
-                ZK(JW, II) = KMIN + (KMAX - KMIN) * ALEAS(0.)
+                ZK(JW, II) = KMIN + (KMAX - KMIN) *RAN_NUM_2
                 ! Horizontal phase speed
-!                CPHA = CMIN + (CMAX - CMIN) * MOD(TT(II, JW)**2, 1.)
-                CPHA = CMIN + (CMAX - CMIN) * ALEAS(0.)
+                CPHA = 0.
+                DO JJ = 1, NA
+                    RAN_NUM_3=MOD(TT(II, JW+3*JJ)**2, 1.)
+                    CPHA = CPHA + &
+                    CMAX*2.*(RAN_NUM_3 -0.5)*SQRT(3.)/SQRT(NA*1.)
+                END DO
+                IF (CPHA.LT.0.)  THEN
+                   CPHA = -1.*CPHA
+                   ZP(JW,II) = ZP(JW,II) + PI
+                ENDIF
        !Online output: linear
                 if (output) CPHA = CMIN + (CMAX - CMIN) * (JO-1)/(NO-1)
@@ -307,25 +284,14 @@
                 ! 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)
+                  + ZK(JW, II) * COS(ZP(JW, II)) * UH(II, LAUNCH) &
+                  + ZK(JW, II) * SIN(ZP(JW, II)) * 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
+                ! epflux_0(JW, II) = epflux_max / REAL(NW) &
+                epflux_0(JW, II) = epflux_max &
+                     * MOD(100. * (UU(II, JW)**2 + VV(II, JW)**2), 1.)
        !Online output: fixed to max
-                if (output) RUW0(JW, II) = RUWMAX
+                if (output) epflux_0(JW, II) = epflux_max
              ENDDO
-          end DO
-       end DO
-    end DO
+   end DO
 
     ! 4. COMPUTE THE FLUXES
@@ -336,39 +302,60 @@
     !
     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) 
+       intr_freq_p(JW, :) = ZO(JW, :) &
+            - ZK(JW, :) * COS(ZP(JW, :)) * UH(:, LAUNCH) &
+            - ZK(JW, :) * SIN(ZP(JW, :)) * VH(:, LAUNCH) 
+    end DO
+
+    IF (gwd_convective_source) THEN
+         DO JW = 1, NW
+       ! VERSION WITH CONVECTIVE SOURCE (designed for Earth)
+
+       ! Vertical velocity at launch level, value to ensure the
+       ! imposed mmt flux factor related to the convective forcing:
+       ! precipitations.
+
+       ! tanh limitation to values above prmax:
+!       WWP(JW, :) = epflux_0(JW, :) &
+!            * (r / cpp / H0 * RLVTT * PRMAX * TANH(PREC(:) / PRMAX))**2
+!       Here, we neglected the kinetic energy providing of the thermodynamic
+!       phase change
+
+!
+
+       ! Factor related to the characteristics of the waves:
+            WWP(JW, :) = WWP(JW, :) * ZK(JW, :)**3 / KMIN / BVLOW(:)  &
+                 / MAX(ABS(intr_freq_p(JW, :)), ZOISEC)**3
+
+      ! Moderation by the depth of the source (dz here):
+            WWP(JW, :) = WWP(JW, :) &
+                 * EXP(- BVLOW(:)**2 / MAX(ABS(intr_freq_p(JW, :)), ZOISEC)**2 &
+                 * ZK(JW, :)**2 * DZ**2)
+
+      ! Put the stress in the right direction:
+            u_epflux_p(JW, :) = intr_freq_p(JW, :) / MAX(ABS(intr_freq_p(JW, :)), ZOISEC)**2 &
+                 * BV(:, LAUNCH) * COS(ZP(JW, :)) * WWP(JW, :)**2
+            v_epflux_p(JW, :) = intr_freq_p(JW, :) / MAX(ABS(intr_freq_p(JW, :)), ZOISEC)**2 &
+                 * BV(:, LAUNCH) * SIN(ZP(JW, :)) * WWP(JW, :)**2
+         end DO
+    ELSE ! VERSION WITHOUT CONVECTIVE SOURCE
        ! 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  !!!!
+         DO JW = 1, NW
+       ! WW is directly a flux, here, not vertical velocity anymore
+            WWP(JW, :) = epflux_0(JW,:)
+            u_epflux_p(JW, :) = COS(ZP(JW, :)) * SIGN(1., intr_freq_p(JW, :)) * epflux_0(JW, :)
+            v_epflux_p(JW, :) = SIN(ZP(JW, :)) * SIGN(1., intr_freq_p(JW, :)) * epflux_0(JW, :)
+
+         end DO
+    ENDIF
     !  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
+    u_epflux_tot(:, LAUNCH) = 0
+    v_epflux_tot(:, LAUNCH) = 0
+    DO JW = 1, NW
+       u_epflux_tot(:, LAUNCH) = u_epflux_tot(:, LAUNCH) + u_epflux_p(JW, :)
+       v_epflux_tot(:, LAUNCH) = v_epflux_tot(:, LAUNCH) + v_epflux_p(JW, :)
+    end DO
 
     !  4.3 Loop over altitudes, with passage from one level to the
@@ -379,12 +366,12 @@
     !Online output
     if (output) then
-        ieq=nlon/2+1
+        ieq=ngrid/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)
+               (ZO(JW, IEQ)/ZK(JW, IEQ)*COS(ZP(JW, IEQ)), JW = 1, NW)
     endif
 
-    DO LL = LAUNCH, KLEV - 1
+    DO LL = LAUNCH, nlayer - 1
 
 
@@ -392,24 +379,22 @@
        ! TO THE NEXT)
        DO JW = 1, NW
-          ZOM(JW, :) = ZOP(JW, :)
+          intr_freq_m(JW, :) = intr_freq_p(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) 
+          intr_freq_p(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 &
+               / MAX(ABS(intr_freq_p(JW, :) + intr_freq_m(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, :))) &
+               MAX(0., SIGN(1., intr_freq_p(JW, :) * intr_freq_m(JW, :))) &
        ! Saturation (Eq. 12)
-               * ABS(ZOP(JW, :))**3 /BV(:, LL+1) & 
+               * ABS(intr_freq_p(JW, :))**3 /BV(:, LL+1) & 
                * EXP(-ZH(:, LL + 1)/H0) * SAT**2*KMIN**2/ZK(JW, :)**4)  
        end DO
@@ -419,32 +404,29 @@
        ! 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
+          u_epflux_p(JW, :) = SIGN(1.,intr_freq_p(JW, :)) * COS(ZP(JW, :)) *  WWP(JW, :)
+          v_epflux_p(JW, :) = SIGN(1.,intr_freq_p(JW, :)) * SIN(ZP(JW, :)) *  WWP(JW, :)
        end DO
  
-       !
-       RUW(:, LL + 1) = 0.
-       RVW(:, LL + 1) = 0.
+       u_epflux_tot(:, LL + 1) = 0.
+       v_epflux_tot(:, LL + 1) = 0.
 
        DO JW = 1, NW
-          RUW(:, LL + 1) = RUW(:, LL + 1) + RUWP(JW, :) 
-          RVW(:, LL + 1) = RVW(:, LL + 1) + RVWP(JW, :) 
+          u_epflux_tot(:, LL + 1) = u_epflux_tot(:, LL + 1) + u_epflux_p(JW, :) 
+          v_epflux_tot(:, LL + 1) = v_epflux_tot(:, LL + 1) + v_epflux_p(JW, :) 
+          EAST_GWSTRESS(:, LL)=EAST_GWSTRESS(:, LL)+MAX(0.,u_epflux_p(JW,:))/FLOAT(NW)
+          WEST_GWSTRESS(:, LL)=WEST_GWSTRESS(:, LL)+MIN(0.,u_epflux_p(JW,:))/FLOAT(NW)
        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)
+            if(u_epflux_p(JW, IEQ).gt.0.) then
+              u_epflux_p(JW, IEQ) = max(u_epflux_p(JW, IEQ), 1.e-99)
             else
-              RUWP(JW, IEQ) = min(RUWP(JW, IEQ), -1.e-99)
+              u_epflux_p(JW, IEQ) = min(u_epflux_p(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)
+                                  (u_epflux_p(JW, IEQ), JW = 1, NW)
        endif
 
@@ -456,28 +438,60 @@
     ! 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)) &
+! Attention, ici c'est le total sur toutes les ondes...
+
+    u_epflux_tot(:, nlayer + 1) = 0.
+    v_epflux_tot(:, nlayer + 1) = 0.
+
+    ! Here, big change compared to FLott version:
+    ! We compensate (u_epflux_tot(:, LAUNCH), ie total emitted upward flux
+    !  over the layers max(1,LAUNCH-3) to LAUNCH-1
+    DO LL = 1, max(1,LAUNCH-3)
+      u_epflux_tot(:, LL) = 0.
+      v_epflux_tot(:, LL) = 0.
+    end DO
+    DO LL = max(2,LAUNCH-2), LAUNCH-1
+       u_epflux_tot(:, LL) = u_epflux_tot(:, LL - 1) + u_epflux_tot(:, LAUNCH) * &
+            (PH(:,LL)-PH(:,LL-1)) / (PH(:,LAUNCH)-PH(:,max(1,LAUNCH-3)))
+       v_epflux_tot(:, LL) = v_epflux_tot(:, LL - 1) + v_epflux_tot(:, LAUNCH) * &
+            (PH(:,LL)-PH(:,LL-1)) / (PH(:,LAUNCH)-PH(:,max(1,LAUNCH-3)))
+       EAST_GWSTRESS(:,LL) = EAST_GWSTRESS(:, LL - 1) + &
+            EAST_GWSTRESS(:, LAUNCH) * (PH(:,LL)-PH(:,LL-1))/ &
+            (PH(:,LAUNCH)-PH(:,max(1,LAUNCH-3)))
+       WEST_GWSTRESS(:,LL) = WEST_GWSTRESS(:, LL - 1) + &
+            WEST_GWSTRESS(:, 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, nlayer - 1
+           WRITE(11,*) ZHbis(IEQ, LL)/1000.,u_epflux_tot(IEQ,LL)
+       end DO
+       CLOSE(11)
+       stop
+    endif
+
+
+    ! 5.2 AR-1 RECURSIVE FORMULA (13) IN VERSION 4
+    !---------------------------------------------
+    DO LL = 1, nlayer
+       d_u(:, LL) = G * (u_epflux_tot(:, LL + 1) - u_epflux_tot(:, LL)) &
             / (PH(:, LL + 1) - PH(:, LL)) * DTIME
-       d_v(:, LL) = G * (RVW(:, LL + 1) - RVW(:, LL)) &
+       d_v(:, LL) = G * (v_epflux_tot(:, LL + 1) - v_epflux_tot(:, 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(:,:)
+
+    ! 5.3 Update tendency of wind with the previous (and saved) values
+    !-----------------------------------------------------------------
+    d_u(:,:) = DTIME/DELTAT/REAL(NW) * d_u(:,:)                       &
+             + (1.-DTIME/DELTAT) * du_nonoro_gwd(:,:)
+    d_v(:,:) = DTIME/DELTAT/REAL(NW) * d_v(:,:)                       &
+             + (1.-DTIME/DELTAT) * dv_nonoro_gwd(:,:)
+    du_nonoro_gwd(:,:) = d_u(:,:)
+    dv_nonoro_gwd(:,:) = d_v(:,:)
 
     ! Cosmetic: evaluation of the cumulated stress
@@ -485,5 +499,5 @@
     ZUSTR(:) = 0.
     ZVSTR(:) = 0.
-    DO LL = 1, KLEV
+    DO LL = 1, nlayer
        ZUSTR(:) = ZUSTR(:) + D_U(:, LL) / g * (PH(:, LL + 1) - PH(:, LL))
        ZVSTR(:) = ZVSTR(:) + D_V(:, LL) / g * (PH(:, LL + 1) - PH(:, LL))
@@ -492,135 +506,28 @@
   END SUBROUTINE NONORO_GWD_RAN
 
-!===================================================================
-!===================================================================
-!===================================================================
-!===================================================================
-
-  REAL FUNCTION ALEAS (R)
+! ===================================================================
+! Subroutines used to write variables of memory in start files       
+! ===================================================================
+
+  SUBROUTINE ini_nonoro_gwd_ran(ngrid,nlayer)
+
   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
+
+      INTEGER, INTENT (in) :: ngrid  ! number of atmospheric columns
+      INTEGER, INTENT (in) :: nlayer ! number of atmospheric layers
+
+         allocate(du_nonoro_gwd(ngrid,nlayer))
+         allocate(dv_nonoro_gwd(ngrid,nlayer))
+
+  END SUBROUTINE ini_nonoro_gwd_ran
+! ----------------------------------
+  SUBROUTINE end_nonoro_gwd_ran
+
+  IMPLICIT NONE
+
+         if (allocated(du_nonoro_gwd)) deallocate(du_nonoro_gwd)
+         if (allocated(dv_nonoro_gwd)) deallocate(dv_nonoro_gwd)
+
+  END SUBROUTINE end_nonoro_gwd_ran
 
 END MODULE nonoro_gwd_ran_mod
Index: /trunk/LMDZ.MARS/libf/phymars/phyetat0_mod.F90
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/phyetat0_mod.F90	(revision 2219)
+++ /trunk/LMDZ.MARS/libf/phymars/phyetat0_mod.F90	(revision 2220)
@@ -17,5 +17,5 @@
                      inquire_dimension, inquire_dimension_length
   use ioipsl_getincom, only : getin
-
+  use nonoro_gwd_ran_mod, only: du_nonoro_gwd, dv_nonoro_gwd
   implicit none
   
@@ -415,4 +415,24 @@
 endif
 
+! Non-orographic gravity waves 
+call get_field("du_nonoro_gwd",du_nonoro_gwd,found,indextime)
+if (.not.found) then
+   write(*,*) "phyetat0: <du_nonoro_gwd> not in file"
+   du_nonoro_gwd(:,:)=0.
+else
+   write(*,*) "phyetat0: Memory of zonal wind tendency due to non-orographic GW"
+   write(*,*) " <du_nonoro_gwd> range:", &
+          minval(du_nonoro_gwd), maxval(du_nonoro_gwd)
+endif
+
+call get_field("dv_nonoro_gwd",dv_nonoro_gwd,found,indextime)
+if (.not.found) then
+   write(*,*) "phyetat0: <dv_nonoro_gwd> not in file"
+   dv_nonoro_gwd(:,:)=0.
+else
+   write(*,*) "phyetat0: Memory of meridional wind tendency due to non-orographic GW"
+   write(*,*) " <dv_nonoro_gwd> range:", &
+          minval(dv_nonoro_gwd), maxval(dv_nonoro_gwd)
+endif
 
 ! tracer on surface
Index: /trunk/LMDZ.MARS/libf/phymars/phyredem.F90
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/phyredem.F90	(revision 2219)
+++ /trunk/LMDZ.MARS/libf/phymars/phyredem.F90	(revision 2220)
@@ -24,5 +24,4 @@
   use comcstfi_h, only: g, mugaz, omeg, rad, rcp
   use time_phylmdz_mod, only: daysec
-
   implicit none
  
@@ -156,4 +155,5 @@
                       put_var, put_field
   use tracer_mod, only: noms ! tracer names
+  use nonoro_gwd_ran_mod, only: du_nonoro_gwd, dv_nonoro_gwd
 
   implicit none
@@ -263,5 +263,10 @@
      call put_field("mem_Mh2o_co2","H2O mass integred into CO2 crystal",mem_Mh2o_co2,time)
   endif
-  
+ 
+  ! Non-orographic gavity waves
+  if (calllott_nonoro) then
+     call put_field("du_nonoro_gwd","Zonal wind tendency due to GW",du_nonoro_gwd,time)
+     call put_field("dv_nonoro_gwd","Meridional wind tendency due to GW",dv_nonoro_gwd,time)
+  endif
   ! Close file
   call close_restartphy
Index: /trunk/LMDZ.MARS/libf/phymars/phys_state_var_init_mod.F90
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/phys_state_var_init_mod.F90	(revision 2219)
+++ /trunk/LMDZ.MARS/libf/phymars/phys_state_var_init_mod.F90	(revision 2220)
@@ -54,4 +54,6 @@
       use watercloud_mod, only: ini_watercloud_mod, &
                                 end_watercloud_mod
+      use nonoro_gwd_ran_mod, only: ini_nonoro_gwd_ran, &
+                                    end_nonoro_gwd_ran
 
       IMPLICIT NONE
@@ -135,4 +137,9 @@
       call ini_watercloud_mod(ngrid,nlayer,nq)
 
+      ! allocate arrays in "nonoro_gwd_ran_mod"
+      call end_nonoro_gwd_ran
+      call ini_nonoro_gwd_ran(ngrid,nlayer)
+
+
       END SUBROUTINE phys_state_var_init
 
Index: /trunk/LMDZ.MARS/libf/phymars/physiq_mod.F
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/physiq_mod.F	(revision 2219)
+++ /trunk/LMDZ.MARS/libf/phymars/physiq_mod.F	(revision 2220)
@@ -430,7 +430,8 @@
       REAL d_u_hin(ngrid,nlayer), d_v_hin(ngrid,nlayer)
       REAL d_t_hin(ngrid,nlayer)
+      REAL east_gwstress(ngrid,nlayer)
+      REAL west_gwstress(ngrid,nlayer)
 c  Diagnostics 2D of gw_nonoro
       REAL zustrhi(ngrid), zvstrhi(ngrid)
-
 c Variables for PBL
       REAL zz1(ngrid)
@@ -1381,5 +1382,6 @@
      &               pdt, pdu, pdv,
      &               zustrhi,zvstrhi,
-     &               d_t_hin, d_u_hin, d_v_hin)
+     &               d_t_hin, d_u_hin, d_v_hin, 
+     &               east_gwstress, west_gwstress)
 
 !  Update tendencies
