Index: /trunk/LMDZ.VENUS/deftank/context_lmdz_physics-78niv.xml
===================================================================
--- /trunk/LMDZ.VENUS/deftank/context_lmdz_physics-78niv.xml	(revision 2579)
+++ /trunk/LMDZ.VENUS/deftank/context_lmdz_physics-78niv.xml	(revision 2580)
@@ -79,4 +79,106 @@
                    long_name="IR rad. at surface" 
                    unit="W/m2" />
+            <field id="col_co2"
+                   long_name="Column density of CO2" 
+                   unit="kg/m2" />
+            <field id="col_co"
+                   long_name="Column density of CO" 
+                   unit="kg/m2" />
+            <field id="col_h2"
+                   long_name="Column density of H2" 
+                   unit="kg/m2" />
+            <field id="col_h2o"
+                   long_name="Column density of H2O" 
+                   unit="kg/m2" />
+            <field id="col_o1d"
+                   long_name="Column density of O(1D)" 
+                   unit="kg/m2" />
+            <field id="col_o"
+                   long_name="Column density of O" 
+                   unit="kg/m2" />
+            <field id="col_o2"
+                   long_name="Column density of O2" 
+                   unit="kg/m2" />
+            <field id="col_o2dg"
+                   long_name="Column density of O2(DG)" 
+                   unit="kg/m2" />
+            <field id="col_o3"
+                   long_name="Column density of O3" 
+                   unit="kg/m2" />
+            <field id="col_h"
+                   long_name="Column density of H" 
+                   unit="kg/m2" />
+            <field id="col_oh"
+                   long_name="Column density of OH" 
+                   unit="kg/m2" />
+            <field id="col_ho2"
+                   long_name="Column density of HO2" 
+                   unit="kg/m2" />
+            <field id="col_h2o2"
+                   long_name="Column density of H2O2" 
+                   unit="kg/m2" />
+            <field id="col_cl"
+                   long_name="Column density of Cl" 
+                   unit="kg/m2" />
+            <field id="col_clo"
+                   long_name="Column density of ClO" 
+                   unit="kg/m2" />
+            <field id="col_cl2"
+                   long_name="Column density of Cl2" 
+                   unit="kg/m2" />
+            <field id="col_hcl"
+                   long_name="Column density of HCl" 
+                   unit="kg/m2" />
+            <field id="col_hocl"
+                   long_name="Column density of HOCl" 
+                   unit="kg/m2" />
+            <field id="col_clco"
+                   long_name="Column density of ClCO" 
+                   unit="kg/m2" />
+            <field id="col_clco3"
+                   long_name="Column density of ClCO3" 
+                   unit="kg/m2" />
+            <field id="col_cocl2"
+                   long_name="Column density of COCl2" 
+                   unit="kg/m2" />
+            <field id="col_s"
+                   long_name="Column density of S" 
+                   unit="kg/m2" />
+            <field id="col_so"
+                   long_name="Column density of SO" 
+                   unit="kg/m2" />
+            <field id="col_so2"
+                   long_name="Column density of SO2" 
+                   unit="kg/m2" />
+            <field id="col_so3"
+                   long_name="Column density of SO3" 
+                   unit="kg/m2" />
+            <field id="col_s2o2"
+                   long_name="Column density of S2O2" 
+                   unit="kg/m2" />
+            <field id="col_ocs"
+                   long_name="Column density of OCS" 
+                   unit="kg/m2" />
+            <field id="col_hso3"
+                   long_name="Column density of HSO3" 
+                   unit="kg/m2" />
+            <field id="col_h2so4"
+                   long_name="Column density of H2SO4" 
+                   unit="kg/m2" />
+            <field id="col_s2"
+                   long_name="Column density of S2" 
+                   unit="kg/m2" />
+            <field id="col_clso2"
+                   long_name="Column density of ClSO2" 
+                   unit="kg/m2" />
+            <field id="col_oscl"
+                   long_name="Column density of OSCl" 
+                   unit="kg/m2" />
+            <field id="col_n2"
+                   long_name="Column density of N2" 
+                   unit="kg/m2" />
+            <field id="col_he"
+                   long_name="Column density of He" 
+                   unit="kg/m2" />
         </field_group>
 
@@ -101,5 +203,5 @@
                    long_name="Vertical wind" 
                    unit="Pa/s" />
-           <field id="vitwz"
+            <field id="vitwz"
                    long_name="Vertical wind (m/s)" 
                    unit="m/s" />
Index: /trunk/LMDZ.VENUS/deftank/context_lmdz_physics-notrac.xml
===================================================================
--- /trunk/LMDZ.VENUS/deftank/context_lmdz_physics-notrac.xml	(revision 2579)
+++ /trunk/LMDZ.VENUS/deftank/context_lmdz_physics-notrac.xml	(revision 2580)
@@ -79,4 +79,106 @@
                    long_name="IR rad. at surface" 
                    unit="W/m2" />
+            <field id="col_co2"
+                   long_name="Column density of CO2" 
+                   unit="kg/m2" />
+            <field id="col_co"
+                   long_name="Column density of CO" 
+                   unit="kg/m2" />
+            <field id="col_h2"
+                   long_name="Column density of H2" 
+                   unit="kg/m2" />
+            <field id="col_h2o"
+                   long_name="Column density of H2O" 
+                   unit="kg/m2" />
+            <field id="col_o1d"
+                   long_name="Column density of O(1D)" 
+                   unit="kg/m2" />
+            <field id="col_o"
+                   long_name="Column density of O" 
+                   unit="kg/m2" />
+            <field id="col_o2"
+                   long_name="Column density of O2" 
+                   unit="kg/m2" />
+            <field id="col_o2dg"
+                   long_name="Column density of O2(DG)" 
+                   unit="kg/m2" />
+            <field id="col_o3"
+                   long_name="Column density of O3" 
+                   unit="kg/m2" />
+            <field id="col_h"
+                   long_name="Column density of H" 
+                   unit="kg/m2" />
+            <field id="col_oh"
+                   long_name="Column density of OH" 
+                   unit="kg/m2" />
+            <field id="col_ho2"
+                   long_name="Column density of HO2" 
+                   unit="kg/m2" />
+            <field id="col_h2o2"
+                   long_name="Column density of H2O2" 
+                   unit="kg/m2" />
+            <field id="col_cl"
+                   long_name="Column density of Cl" 
+                   unit="kg/m2" />
+            <field id="col_clo"
+                   long_name="Column density of ClO" 
+                   unit="kg/m2" />
+            <field id="col_cl2"
+                   long_name="Column density of Cl2" 
+                   unit="kg/m2" />
+            <field id="col_hcl"
+                   long_name="Column density of HCl" 
+                   unit="kg/m2" />
+            <field id="col_hocl"
+                   long_name="Column density of HOCl" 
+                   unit="kg/m2" />
+            <field id="col_clco"
+                   long_name="Column density of ClCO" 
+                   unit="kg/m2" />
+            <field id="col_clco3"
+                   long_name="Column density of ClCO3" 
+                   unit="kg/m2" />
+            <field id="col_cocl2"
+                   long_name="Column density of COCl2" 
+                   unit="kg/m2" />
+            <field id="col_s"
+                   long_name="Column density of S" 
+                   unit="kg/m2" />
+            <field id="col_so"
+                   long_name="Column density of SO" 
+                   unit="kg/m2" />
+            <field id="col_so2"
+                   long_name="Column density of SO2" 
+                   unit="kg/m2" />
+            <field id="col_so3"
+                   long_name="Column density of SO3" 
+                   unit="kg/m2" />
+            <field id="col_s2o2"
+                   long_name="Column density of S2O2" 
+                   unit="kg/m2" />
+            <field id="col_ocs"
+                   long_name="Column density of OCS" 
+                   unit="kg/m2" />
+            <field id="col_hso3"
+                   long_name="Column density of HSO3" 
+                   unit="kg/m2" />
+            <field id="col_h2so4"
+                   long_name="Column density of H2SO4" 
+                   unit="kg/m2" />
+            <field id="col_s2"
+                   long_name="Column density of S2" 
+                   unit="kg/m2" />
+            <field id="col_clso2"
+                   long_name="Column density of ClSO2" 
+                   unit="kg/m2" />
+            <field id="col_oscl"
+                   long_name="Column density of OSCl" 
+                   unit="kg/m2" />
+            <field id="col_n2"
+                   long_name="Column density of N2" 
+                   unit="kg/m2" />
+            <field id="col_he"
+                   long_name="Column density of He" 
+                   unit="kg/m2" />
         </field_group>
 
@@ -101,5 +203,5 @@
                    long_name="Vertical wind" 
                    unit="Pa/s" />
-           <field id="vitwz"
+            <field id="vitwz"
                    long_name="Vertical wind (m/s)" 
                    unit="m/s" />
Index: /trunk/LMDZ.VENUS/deftank/context_lmdz_physics.xml
===================================================================
--- /trunk/LMDZ.VENUS/deftank/context_lmdz_physics.xml	(revision 2579)
+++ /trunk/LMDZ.VENUS/deftank/context_lmdz_physics.xml	(revision 2580)
@@ -79,4 +79,106 @@
                    long_name="IR rad. at surface" 
                    unit="W/m2" />
+            <field id="col_co2"
+                   long_name="Column density of CO2" 
+                   unit="kg/m2" />
+            <field id="col_co"
+                   long_name="Column density of CO" 
+                   unit="kg/m2" />
+            <field id="col_h2"
+                   long_name="Column density of H2" 
+                   unit="kg/m2" />
+            <field id="col_h2o"
+                   long_name="Column density of H2O" 
+                   unit="kg/m2" />
+            <field id="col_o1d"
+                   long_name="Column density of O(1D)" 
+                   unit="kg/m2" />
+            <field id="col_o"
+                   long_name="Column density of O" 
+                   unit="kg/m2" />
+            <field id="col_o2"
+                   long_name="Column density of O2" 
+                   unit="kg/m2" />
+            <field id="col_o2dg"
+                   long_name="Column density of O2(DG)" 
+                   unit="kg/m2" />
+            <field id="col_o3"
+                   long_name="Column density of O3" 
+                   unit="kg/m2" />
+            <field id="col_h"
+                   long_name="Column density of H" 
+                   unit="kg/m2" />
+            <field id="col_oh"
+                   long_name="Column density of OH" 
+                   unit="kg/m2" />
+            <field id="col_ho2"
+                   long_name="Column density of HO2" 
+                   unit="kg/m2" />
+            <field id="col_h2o2"
+                   long_name="Column density of H2O2" 
+                   unit="kg/m2" />
+            <field id="col_cl"
+                   long_name="Column density of Cl" 
+                   unit="kg/m2" />
+            <field id="col_clo"
+                   long_name="Column density of ClO" 
+                   unit="kg/m2" />
+            <field id="col_cl2"
+                   long_name="Column density of Cl2" 
+                   unit="kg/m2" />
+            <field id="col_hcl"
+                   long_name="Column density of HCl" 
+                   unit="kg/m2" />
+            <field id="col_hocl"
+                   long_name="Column density of HOCl" 
+                   unit="kg/m2" />
+            <field id="col_clco"
+                   long_name="Column density of ClCO" 
+                   unit="kg/m2" />
+            <field id="col_clco3"
+                   long_name="Column density of ClCO3" 
+                   unit="kg/m2" />
+            <field id="col_cocl2"
+                   long_name="Column density of COCl2" 
+                   unit="kg/m2" />
+            <field id="col_s"
+                   long_name="Column density of S" 
+                   unit="kg/m2" />
+            <field id="col_so"
+                   long_name="Column density of SO" 
+                   unit="kg/m2" />
+            <field id="col_so2"
+                   long_name="Column density of SO2" 
+                   unit="kg/m2" />
+            <field id="col_so3"
+                   long_name="Column density of SO3" 
+                   unit="kg/m2" />
+            <field id="col_s2o2"
+                   long_name="Column density of S2O2" 
+                   unit="kg/m2" />
+            <field id="col_ocs"
+                   long_name="Column density of OCS" 
+                   unit="kg/m2" />
+            <field id="col_hso3"
+                   long_name="Column density of HSO3" 
+                   unit="kg/m2" />
+            <field id="col_h2so4"
+                   long_name="Column density of H2SO4" 
+                   unit="kg/m2" />
+            <field id="col_s2"
+                   long_name="Column density of S2" 
+                   unit="kg/m2" />
+            <field id="col_clso2"
+                   long_name="Column density of ClSO2" 
+                   unit="kg/m2" />
+            <field id="col_oscl"
+                   long_name="Column density of OSCl" 
+                   unit="kg/m2" />
+            <field id="col_n2"
+                   long_name="Column density of N2" 
+                   unit="kg/m2" />
+            <field id="col_he"
+                   long_name="Column density of He" 
+                   unit="kg/m2" />
         </field_group>
 
@@ -101,5 +203,5 @@
                    long_name="Vertical wind" 
                    unit="Pa/s" />
-           <field id="vitwz"
+            <field id="vitwz"
                    long_name="Vertical wind (m/s)" 
                    unit="m/s" />
Index: /trunk/LMDZ.VENUS/deftank/physiq-96x96x50-chemistry.def
===================================================================
--- /trunk/LMDZ.VENUS/deftank/physiq-96x96x50-chemistry.def	(revision 2579)
+++ /trunk/LMDZ.VENUS/deftank/physiq-96x96x50-chemistry.def	(revision 2580)
@@ -97,15 +97,18 @@
 ## Thermospheric options 
 ##~~~~~~~~~~~~~~~~~~~~~~
+## Tuning of photochemistry for oxygen production ?
+tuneupperatm = n
+
 #Method to include solar variability?
 #0-> Old method   1-> Variability with E10.7 as observed
 solvarmod=0
 
-## (Solar min=1996.4 ave=1993.4 max=1990.6) ; Only used if solvarmod=0
-solarcondate = 1993.4
+## (Solar min=70 ave=140 max=300)
+fixed_euv_value = 140.
 
 # value for the UV heating efficiency
 ##(experimental values between 0.19 and 0.23, lower values may
 ## be used to compensate for low 15 um cooling)
-euveff = 0.21       
+euveff = 0.20       
 #
 #
Index: /trunk/LMDZ.VENUS/deftank/physiq-96x96x50-notrac.def
===================================================================
--- /trunk/LMDZ.VENUS/deftank/physiq-96x96x50-notrac.def	(revision 2579)
+++ /trunk/LMDZ.VENUS/deftank/physiq-96x96x50-notrac.def	(revision 2580)
@@ -97,15 +97,18 @@
 ## Thermospheric options 
 ##~~~~~~~~~~~~~~~~~~~~~~
+## Tuning of photochemistry for oxygen production ?
+tuneupperatm = n
+
 #Method to include solar variability?
 #0-> Old method   1-> Variability with E10.7 as observed
 solvarmod=0
 
-## (Solar min=1996.4 ave=1993.4 max=1990.6) ; Only used if solvarmod=0
-solarcondate = 1993.4
+## (Solar min=70 ave=140 max=300)
+fixed_euv_value = 140.
 
 # value for the UV heating efficiency
 ##(experimental values between 0.19 and 0.23, lower values may
 ## be used to compensate for low 15 um cooling)
-euveff = 0.21       
+euveff = 0.20       
 #
 #
Index: /trunk/LMDZ.VENUS/deftank/physiq-96x96x78-chemistry.def
===================================================================
--- /trunk/LMDZ.VENUS/deftank/physiq-96x96x78-chemistry.def	(revision 2579)
+++ /trunk/LMDZ.VENUS/deftank/physiq-96x96x78-chemistry.def	(revision 2580)
@@ -91,5 +91,5 @@
 callnirco2=y
 ##CO2 IR Absorption model
-nircorr=1
+nircorr=0
 ##include thermosphere or not
 callthermos=y
@@ -97,9 +97,12 @@
 ## Thermospheric options 
 ##~~~~~~~~~~~~~~~~~~~~~~
+## Tuning of photochemistry for oxygen production ?
+tuneupperatm = y
+
 #Method to include solar variability?
 #0-> Old method   1-> Variability with E10.7 as observed
 solvarmod=0
 
-## (Solar min=80 ave=140 max=300)
+## (Solar min=70 ave=140 max=300)
 fixed_euv_value = 140.
 
Index: /trunk/LMDZ.VENUS/deftank/physiq-96x96x90-chemistry.def
===================================================================
--- /trunk/LMDZ.VENUS/deftank/physiq-96x96x90-chemistry.def	(revision 2579)
+++ /trunk/LMDZ.VENUS/deftank/physiq-96x96x90-chemistry.def	(revision 2580)
@@ -91,5 +91,5 @@
 callnirco2=y
 ##CO2 IR Absorption model
-nircorr=1
+nircorr=0
 ##include thermosphere or not
 callthermos=y
@@ -97,9 +97,12 @@
 ## Thermospheric options 
 ##~~~~~~~~~~~~~~~~~~~~~~
+## Tuning of photochemistry for oxygen production ?
+tuneupperatm = y
+
 #Method to include solar variability?
 #0-> Old method   1-> Variability with E10.7 as observed
 solvarmod=0
 
-## (Solar min=80 ave=140 max=300)
+## (Solar min=70 ave=140 max=300)
 fixed_euv_value = 140.
 
@@ -107,5 +110,5 @@
 ##(experimental values between 0.19 and 0.23, lower values may
 ## be used to compensate for low 15 um cooling)
-euveff = 0.21       
+euveff = 0.20       
 #
 #
Index: /trunk/LMDZ.VENUS/deftank/physiq.def
===================================================================
--- /trunk/LMDZ.VENUS/deftank/physiq.def	(revision 2579)
+++ /trunk/LMDZ.VENUS/deftank/physiq.def	(revision 2580)
@@ -97,15 +97,18 @@
 ## Thermospheric options 
 ##~~~~~~~~~~~~~~~~~~~~~~
+## Tuning of photochemistry for oxygen production ?
+tuneupperatm = n
+
 #Method to include solar variability?
 #0-> Old method   1-> Variability with E10.7 as observed
 solvarmod=0
 
-## (Solar min=1996.4 ave=1993.4 max=1990.6) ; Only used if solvarmod=0
-solarcondate = 1993.4
+## (Solar min=70 ave=140 max=300)
+fixed_euv_value = 140.
 
 # value for the UV heating efficiency
 ##(experimental values between 0.19 and 0.23, lower values may
 ## be used to compensate for low 15 um cooling)
-euveff = 0.21       
+euveff = 0.20       
 #
 #
Index: /trunk/LMDZ.VENUS/libf/phyvenus/blendrad.F
===================================================================
--- /trunk/LMDZ.VENUS/libf/phyvenus/blendrad.F	(revision 2579)
+++ /trunk/LMDZ.VENUS/libf/phyvenus/blendrad.F	(revision 2580)
@@ -36,5 +36,5 @@
       integer l, ig
       real alpha, alpha2
-      real, parameter :: p_lowup = 10.e2
+      real, parameter :: p_lowup = 1.e3
 
 c
Index: /trunk/LMDZ.VENUS/libf/phyvenus/clesphys.h
===================================================================
--- /trunk/LMDZ.VENUS/libf/phyvenus/clesphys.h	(revision 2579)
+++ /trunk/LMDZ.VENUS/libf/phyvenus/clesphys.h	(revision 2580)
@@ -10,5 +10,5 @@
        LOGICAL cycle_diurne,soil_model
        LOGICAL ok_orodr,ok_orolf,ok_gw_nonoro
-       LOGICAL ok_kzmin
+       LOGICAL ok_kzmin,tuneupperatm
        LOGICAL callnlte,callnirco2,callthermos
        LOGICAL ok_cloud, ok_chem, reinit_trac, ok_sedim
@@ -27,5 +27,5 @@
        COMMON/clesphys_l/ cycle_diurne, soil_model,                     &
      &     ok_orodr, ok_orolf, ok_gw_nonoro, ok_kzmin,                  &
-     &     callnlte,callnirco2,callthermos,                             &
+     &     tuneupperatm,callnlte,callnirco2,callthermos,                &
      &     ok_cloud, ok_chem, reinit_trac, ok_sedim,                    &
      &     ok_clmain, physideal, startphy_file
Index: /trunk/LMDZ.VENUS/libf/phyvenus/conf_phys.F90
===================================================================
--- /trunk/LMDZ.VENUS/libf/phyvenus/conf_phys.F90	(revision 2579)
+++ /trunk/LMDZ.VENUS/libf/phyvenus/conf_phys.F90	(revision 2580)
@@ -472,9 +472,17 @@
 !Config Key  = euveff
 !Config Desc =
-!Config Def  = 0.21
-!Config Help =
-!
-  euveff = 0.21
+!Config Def  = 0.22
+!Config Help =
+!
+  euveff = 0.22
   call getin('euveff',euveff)
+
+!Config Key  = tuneupperatm
+!Config Desc =
+!Config Def  = .false.
+!Config Help =
+!
+  tuneupperatm = .false.
+  call getin('tuneupperatm',tuneupperatm)
 
 !
@@ -538,4 +546,5 @@
   write(lunout,*)' fixed_euv_value = ',fixed_euv_value
   write(lunout,*)' euveff = ',euveff
+  write(lunout,*)' tuneupperatm = ',tuneupperatm
 
   end subroutine conf_phys
Index: /trunk/LMDZ.VENUS/libf/phyvenus/nirco2abs.F
===================================================================
--- /trunk/LMDZ.VENUS/libf/phyvenus/nirco2abs.F	(revision 2579)
+++ /trunk/LMDZ.VENUS/libf/phyvenus/nirco2abs.F	(revision 2580)
@@ -88,4 +88,7 @@
       real ztim1,ztim2,ztim3,step
 
+      logical onepeak
+      parameter (onepeak=.false.)
+c     parameter (onepeak=.true.)
 c
 c   local saved variables
@@ -95,6 +98,11 @@
       integer,save :: io=0 ! index of "o" tracer
 
+ccc=================================================
 cccc     parameters for CO2 heating fit
-c
+ccc=================================================
+
+c--------------------------------------------------
+c One-peak martian-type fit => Gabriella (2014+)
+c--------------------------------------------------
 c     n_a  =  heating rate for Venusian day at p0, r0, mu =0 [K day-1]
 c     Here p0 = p_cloud top [Pa]
@@ -104,6 +112,5 @@
       real n_a, n_p0, n_b, p_ctop
 
-    
-cc "Nominal" values used in Gilli+2'17
+cc "Nominal" values used in Gilli+2017
 c       parameter (n_a = 18.13/86400.0)     !c     K/Eday  ---> K/sec   
 c       parameter (p_ctop=13.2e2)
@@ -111,4 +118,5 @@
 
 cc "New" values used to improve SPICAV/SOIR Temperature comparision (D.Quirino)
+cc Gilli+2021
        parameter (n_a = 15.92/86400.0)     !c     K/Eday  ---> K/sec   
        parameter (p_ctop=19.85e2)
@@ -119,5 +127,60 @@
 C       parameter (n_p0=0.01)  
 c       parameter (n_b = 1.3)
-   
+
+c--------------------------------------------------
+c Multi-peaks Roldan-type fit => Laura (2013)
+c New paramaters (Param9*0.5) => Enora (2021)
+c--------------------------------------------------
+c ENORA FINE TUNING used for VCD 1.1 
+c  (fit to fig 12 Roldan-2000)
+      real n_coFB, n_aFB, n_bFB, n_p0FB, n_eFB
+      real n_coISO, n_aISO, n_bISO, n_p0ISO, n_eISO
+      real n_coFH, n_aFH, n_bFH, n_p0FH, n_eFH
+      real n_co43, n_a43, n_b43, n_p043, n_e43
+      real n_co43b, n_a43b, n_b43b, n_p043b, n_e43b
+      real n_conir, n_anir, n_bnir, n_p0nir, n_enir
+
+      parameter (n_coFB=119./86400.0)      !c     K/Eday ---> K/sec
+      parameter (n_aFB=0.185)
+      parameter (n_bFB=3.7)
+      parameter (n_p0FB=2.9e-4)
+      parameter (n_eFB=0.76)
+
+      parameter (n_coISO=265./86400.0)      !c     K/Eday ---> K/sec
+      parameter (n_aISO=0.313)
+      parameter (n_bISO=1.65)
+      parameter (n_p0ISO=0.076)
+      parameter (n_eISO=0.99)
+
+      parameter (n_coFH=2.5/86400.0)      !c     K/Eday ---> K/sec
+      parameter (n_aFH=3.98)
+      parameter (n_bFH=2.9)
+      parameter (n_p0FH=0.17)
+      parameter (n_eFH=2.16)
+
+      parameter (n_co43=55./86400.0)      !c     K/Eday ---> K/sec
+      parameter (n_a43=0.625)
+      parameter (n_b43=2.6)
+      parameter (n_p043=0.043)
+      parameter (n_e43=1.654)
+
+!     parameter (n_co43b=100./86400.0)      !c     K/Eday ---> K/sec
+! => fine tuning: not affected by the *0.5 below (see ENORA FINE TUNING)
+      parameter (n_co43b=200./86400.0)      !c     K/Eday ---> K/sec
+      parameter (n_a43b=5.5)
+      parameter (n_b43b=2.3)
+      parameter (n_p043b=1.)
+      parameter (n_e43b=0.4)
+
+      parameter (n_conir=6.5/86400.0)      !c     K/Eday ---> K/sec
+      parameter (n_anir=35.65)
+      parameter (n_bnir=2.1)
+      parameter (n_p0nir=0.046)
+      parameter (n_enir=0.9)
+
+      real :: picFB(nlon,nlev), picISO(nlon,nlev), picFH(nlon,nlev)
+      real :: pic43(nlon,nlev), pic43b(nlon,nlev), picnir(nlon,nlev)
+
+ccc=================================================
 
 c     Variables added to implement NLTE correction factor (feb 2011)
@@ -126,17 +189,6 @@
       real    p2011,cociente1,merge
       real    cor0,oco2gcm
-!!!!
-c      real :: pic27(nlon,nlev), pic27b(nlon,nlev)
-c      real :: pic43(nlon,nlev), picnir(nlon,nlev)
-
-c     co2heat is the heating by CO2 at p_ctop=13.2e2 for a zero zenithal angle.
-
-      co2heat0=n_a*(0.72/dist_sol)**2     
-
-CCCCCC   TEST: reduce by X% nir Heating
-c      co2heat0  = co2heat0 * 0.8
-
-c----------------------------------------------------------------------
-      
+
+c----------------------------------------------------------------------  
 c     Initialisation
 c     --------------
@@ -159,6 +211,11 @@
         firstcall=.false.
       endif
-
-      
+c     --------------
+c     co2heat0 is correction for dist_sol (is 1 for Venus)
+      co2heat0=(0.7233/dist_sol)**2
+
+      pdtnirco2(:,:)=0.
+c----------------------------------------------------------------------  
+
 c      
 c     Simple calcul for a given sun incident angle (if cycle_diurne=T)
@@ -170,18 +227,16 @@
             zmu(ig)=sqrt(1224.*mu0(ig)*mu0(ig)+1.)/35.
 
-           
+c---------------------------
+           if (onepeak) then
+c---------------------------
             if(nircorr.eq.1) then
                do l=1,nlev
                   pyy(l)=nplay(ig,l)
                enddo
-
                call interpnir(cor1,pyy,nlev,corgcm,pres1d,npres)
                call interpnir(oldoco2,pyy,nlev,oco21d,pres1d,npres)
                call interpnir(alfa2,pyy,nlev,alfa,pres1d,npres)
-               
             endif
-
             do l=1,nlev
-     
 c           Calculations for the O/CO2 correction
                if(nircorr.eq.1) then
@@ -213,12 +268,57 @@
 
               if(fract(ig).gt.0.) pdtnirco2(ig,l)=
-     &             co2heat0*sqrt((p_ctop*zmu(ig))/nplay(ig,l))
+     &             co2heat0*n_a*sqrt((p_ctop*zmu(ig))/nplay(ig,l))
      &             /(1.+n_p0/nplay(ig,l))**n_b
 c           Corrections from tabulation
      $              * cor1(l) * p2011
               
-          enddo
-         enddo
-	 
+            enddo !nlev
+c---------------------------
+           else  ! multipeak
+c---------------------------
+            do l=1,nlev
+               if(fract(ig).gt.0.) then 
+                   picFB(ig,l)=n_coFB
+     &              *((n_aFB/nplay(ig,l))**n_eFB)
+     &		    *zmu(ig)**0.82
+     &             /(1.+n_p0FB/nplay(ig,l))**n_bFB
+
+                   picISO(ig,l)=n_coISO
+     &              *((n_aISO/nplay(ig,l))**n_eISO)
+     &		    *zmu(ig)**0.55
+     &             /(1.+n_p0ISO/nplay(ig,l))**n_bISO
+
+                   picFH(ig,l)=n_coFH
+     &              *((n_aFH/nplay(ig,l))**n_eFH)
+     &		    *zmu(ig)**0.55
+     &             /(1.+n_p0FH/nplay(ig,l))**n_bFH
+
+                   pic43(ig,l)=n_co43
+     &              *((n_a43/nplay(ig,l))**n_e43)
+     &		    *zmu(ig)**0.55
+     &             /(1.+n_p043/nplay(ig,l))**n_b43
+
+                   pic43b(ig,l)=n_co43b
+     &              *((n_a43b/nplay(ig,l))**n_e43b)
+     &		    *zmu(ig)**0.55
+     &             /(1.+n_p043b/nplay(ig,l))**n_b43b
+
+                   picnir(ig,l)=n_conir
+     &              *((n_anir/nplay(ig,l))**n_enir)
+     &		    *zmu(ig)**0.55
+     &             /(1.+n_p0nir/nplay(ig,l))**n_bnir
+
+	           pdtnirco2(ig,l)=co2heat0*
+     &           (picFB(ig,l)+picISO(ig,l)+picFH(ig,l)+pic43(ig,l)
+     &           +pic43b(ig,l)+picnir(ig,l))*0.5  ! *0.5 = ENORA FINE TUNING
+		  
+               endif
+            enddo !nlev
+c---------------------------
+           endif
+c---------------------------
+         enddo  !nlon
+
+
 c     Averaging over diurnal cycle (if diurnal=F)
 c     -------------------------------------------
@@ -237,23 +337,37 @@
      &                  mu0_int,fract_int)
 
-            do ig=1,nlon
+          do ig=1,nlon
                zmu(ig)=sqrt(1224.*mu0_int(ig)*mu0_int(ig)+1.)/35.
 
-               if(nircorr.eq.1) then
-                  do l=1,nlev
-                     pyy(l)=nplay(ig,l)
-                  enddo
-                 call interpnir(cor1,pyy,nlev,corgcm,pres1d,npres)
-                 call interpnir(oldoco2,pyy,nlev,oco21d,pres1d,npres)
-                 call interpnir(alfa2,pyy,nlev,alfa,pres1d,npres)
-               endif
-c
-
+c---------------------------
+           if (onepeak) then
+c---------------------------
+            if(nircorr.eq.1) then
                do l=1,nlev
+                  pyy(l)=nplay(ig,l)
+               enddo
+               call interpnir(cor1,pyy,nlev,corgcm,pres1d,npres)
+               call interpnir(oldoco2,pyy,nlev,oco21d,pres1d,npres)
+               call interpnir(alfa2,pyy,nlev,alfa,pres1d,npres)
+            endif
+            do l=1,nlev
 c           Calculations for the O/CO2 correction
                if(nircorr.eq.1) then
                   cor0=1./(1.+n_p0/nplay(ig,l))**n_b
-                  oco2gcm=pq(ig,l,io)/pq(ig,l,ico2)
+                  if(pq(ig,l,ico2) .gt. 1.e-6) then
+                     oco2gcm=pq(ig,l,io)/pq(ig,l,ico2)
+                     ! handle the rare cases when pq(ig,l,io)<0
+                     if (pq(ig,l,io).lt.0) then
+                       write(*,*) "nirco2abs: warning ig=",ig," l=",l,
+     &                            " pq(ig,l,io)=",pq(ig,l,io)
+                       oco2gcm=1.e6
+                     endif
+                  else
+                     oco2gcm=1.e6
+                  endif
                   cociente1=oco2gcm/oldoco2(l)
+                  
+c                  WRITE(*,*) "nirco2abs line 211", l, cociente1
+
                   merge=alog10(cociente1)*alfa2(l)+alog10(cor0)*
      $                 (1.-alfa2(l))
@@ -266,17 +380,61 @@
                endif
 
-               if(fract_int(ig).gt.0.) pdtnirco2(ig,l)=
+              if(fract(ig).gt.0.) pdtnirco2(ig,l)=
      &              pdtnirco2(ig,l) + (1/float(nstep))*
-     &              co2heat0*sqrt((p_ctop*zmu(ig))/nplay(ig,l))
-     &              /(1.+n_p0/nplay(ig,l))**n_b
-!     Corrections from tabulation
-     $              * cor1(l) * p2011
-
-               enddo
-            enddo
-         end do
+     &             co2heat0*n_a*sqrt((p_ctop*zmu(ig))/nplay(ig,l))
+     &             /(1.+n_p0/nplay(ig,l))**n_b
+c           Corrections from tabulation
+     $              * cor1(l) * p2011 
+              
+            enddo !nlev
+c---------------------------
+           else  ! multipeak
+c---------------------------
+            do l=1,nlev
+               if(fract(ig).gt.0.) then 
+                   picFB(ig,l)=n_coFB
+     &              *((n_aFB/nplay(ig,l))**n_eFB)
+     &		    *zmu(ig)**0.82
+     &             /(1.+n_p0FB/nplay(ig,l))**n_bFB
+
+                   picISO(ig,l)=n_coISO
+     &              *((n_aISO/nplay(ig,l))**n_eISO)
+     &		    *zmu(ig)**0.55
+     &             /(1.+n_p0ISO/nplay(ig,l))**n_bISO
+
+                   picFH(ig,l)=n_coFH
+     &              *((n_aFH/nplay(ig,l))**n_eFH)
+     &		    *zmu(ig)**0.55
+     &             /(1.+n_p0FH/nplay(ig,l))**n_bFH
+
+                   pic43(ig,l)=n_co43
+     &              *((n_a43/nplay(ig,l))**n_e43)
+     &		    *zmu(ig)**0.55
+     &             /(1.+n_p043/nplay(ig,l))**n_b43
+
+                   pic43b(ig,l)=n_co43b
+     &              *((n_a43b/nplay(ig,l))**n_e43b)
+     &		    *zmu(ig)**0.55
+     &             /(1.+n_p043b/nplay(ig,l))**n_b43b
+
+                   picnir(ig,l)=n_conir
+     &              *((n_anir/nplay(ig,l))**n_enir)
+     &		    *zmu(ig)**0.55
+     &             /(1.+n_p0nir/nplay(ig,l))**n_bnir
+
+	           pdtnirco2(ig,l)=
+     &               pdtnirco2(ig,l)+(1/float(nstep))*co2heat0*
+     &           (picFB(ig,l)+picISO(ig,l)+picFH(ig,l)+pic43(ig,l)
+     &           +pic43b(ig,l)+picnir(ig,l))*0.5  ! *0.5 = ENORA FINE TUNING
+		  
+               endif
+            enddo !nlev
+c---------------------------
+           endif
+c---------------------------
+          enddo  !nlon
+         enddo  !nstep
       
-
-      END IF  
+      END IF  ! diurnal cycle
 
       return
Index: /trunk/LMDZ.VENUS/libf/phyvenus/nlte_setup.F
===================================================================
--- /trunk/LMDZ.VENUS/libf/phyvenus/nlte_setup.F	(revision 2579)
+++ /trunk/LMDZ.VENUS/libf/phyvenus/nlte_setup.F	(revision 2580)
@@ -177,9 +177,13 @@
                                 !! k19 & k20
 
-c      k20x = 3.d-12
+c     k20x = 3.d-12
 c  TEST GG: double the values of Kvv as recently found by Sharma et al.2014
-      k20x = 6.d-12   
+c     k20x = 6.d-12   
 c  TEST GG: use the minimum value of the experimental bracket's values [1-6]
 c      k20x = 1.d-12 
+c--------------------
+ccc VCD 1.1 tuning
+      k20x = 5.d-12
+c--------------------
       k20xc = k20x * rf20
       k20xb = 2.d0 * k20xc
Index: /trunk/LMDZ.VENUS/libf/phyvenus/nonoro_gwd_ran_mod.F90
===================================================================
--- /trunk/LMDZ.VENUS/libf/phyvenus/nonoro_gwd_ran_mod.F90	(revision 2580)
+++ /trunk/LMDZ.VENUS/libf/phyvenus/nonoro_gwd_ran_mod.F90	(revision 2580)
@@ -0,0 +1,472 @@
+MODULE nonoro_gwd_ran_mod
+
+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
+      REAL, allocatable, save :: east_gwstress(:, :) ! Eastward stress profile
+      REAL, allocatable, save :: west_gwstress(:, :) ! Westward stress profile
+!$OMP THREADPRIVATE(du_nonoro_gwd,dv_nonoro_gwd,east_gwstress,west_gwstress)
+CONTAINS
+
+      SUBROUTINE nonoro_gwd_ran(ngrid, nlayer, dtime, pp, &
+                   pn2, presnivs, pt, pu, pv, &
+                   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)
+!                      D.BARDET    01/2020  - reproductibility of
+!                                           the launching altitude
+!                                           calculation 
+!                                           - wave characteristic
+!                                           calculation using MOD
+!                                           - adding east_gwstress
+!                                           and west_gwstress variables    
+!
+! ADAPTED FOR GENERIC  D.BARDET    01/2020
+! READAPTED FOR VENUS  S.LEBONNOIS 08/2021
+!---------------------------------------------------------------------------------
+
+     USE ioipsl_getin_p_mod, ONLY : getin_p
+     use geometry_mod, only: cell_area
+!#ifdef CPP_XIOS
+!     use xios_output_mod, only: send_xios_field
+!#endif
+ implicit none
+#include "YOMCST.h"
+
+! 0. DECLARATIONS:
+
+     ! 0.1 INPUTS
+
+     INTEGER, intent(in):: ngrid, nlayer
+     REAL, intent(in):: dtime              ! Time step of the Physics
+     REAL, intent(in):: pp(ngrid, nlayer)  ! Pressure at full levels
+     REAL, intent(in):: pn2(ngrid, nlayer) ! N2 (BV^2) at 1/2 levels
+     REAL, intent(in):: presnivs(nlayer)   ! Approximate pressure for reproductible launching altitude
+     REAL, intent(in):: pt(ngrid, nlayer)  ! Temperature at full levels
+     REAL, intent(in):: pu(ngrid, nlayer)  ! Zonal wind at full levels
+     REAL, intent(in):: pv(ngrid, nlayer)  ! Meridional wind at full levels
+      
+     ! 0.2 OUTPUTS
+
+     REAL, intent(out):: zustr(ngrid)         ! Zonal surface stress
+     REAL, intent(out):: zvstr(ngrid)         ! Meridional surface stress
+     REAL, intent(out):: d_t(ngrid, nlayer)   ! Tendency on temperature
+     REAL, intent(out):: d_u(ngrid, nlayer)   ! Tendency on zonal wind
+     REAL, intent(out):: d_v(ngrid, nlayer)   ! Tendency on meridional wind
+
+     ! 0.3 INTERNAL ARRAYS
+
+     REAL :: tt(ngrid, nlayer) ! Temperature at full levels
+     REAL :: uu(ngrid, nlayer) ! Zonal wind at full levels
+     REAL :: vv(ngrid, nlayer) ! Meridional wind at full levels
+
+     INTEGER ii, jj, ll 
+
+     ! 0.3.0 Time scale of the like cycle of the waves parametrized
+     REAL, parameter:: deltat = 24. * 3600.
+
+     ! 0.3.1 Gravity waves specifications
+     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:: na = 5 ! Number of realizations to get the phase speed
+     INTEGER, parameter:: nw = nk * np *no ! Total numbers of gravity waves
+     INTEGER jk, jp, jo, jw
+     REAL kstar                  ! Control value to constrain the min horizontal wavenumber by the horizontal resolution
+
+!!! TEST Kobs  and BestFit Diogo
+     REAL, parameter:: kmax = 1.e-4 ! Max horizontal wavenumber (lambda min 50 km)
+! generic : kmax=N/u, u(=30) zonal wind velocity at launch
+     REAL, parameter:: kmin = 1.e-5 ! Min horizontal wavenumber (lambda max 628 km)
+! generic : kmin=1/sqrt(DxDy) Dx and Dy horizontal grid
+
+!     REAL, parameter:: cmax = 61.   ! Max horizontal absolute phase velocity  TestGW6
+!----------------------------------------
+! VCD 2.0 tuning
+     REAL, parameter:: cmax = 301.   ! Max horizontal absolute phase velocity  TestGW6
+!----------------------------------------
+
+! generic : cmax=zonal wind at launch
+     REAL, parameter:: cmin = 1.    ! Min horizontal absolute phase velocity
+
+!!! Nominal as Gilli2017
+     REAL, parameter:: epflux_max = 0.005  ! Max EP flux value at launching altitude (previous name: RUWMAX)
+
+     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)
+     INTEGER launch                       ! Launching altitude
+     REAL, parameter:: xlaunch = 5e-3     ! Value for top of cloud convective region
+
+     ! 0.3.2 Parameters of waves dissipations
+     REAL, parameter:: sat   = 0.85   ! saturation parameter
+     REAL, parameter:: rdiss = 0.1    ! coefficient of dissipation
+     REAL, parameter:: zoisec = 1.e-8 ! 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 = 5e5      ! Reference pressure [Pa]   ! VENUS: cloud layer
+     REAL, parameter:: tr = 300.     ! Reference temperature [K] ! VENUS: cloud layer
+     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-8  ! Security to avoid division by 0 pressure
+     REAL bv(ngrid, nlayer + 1)      ! Brunt Vaisala freq. at 1/2 levels
+     REAL, parameter:: bvsec = 1.e-5 ! Security to avoid negative BV
+     REAL href(nlayer + 1)           ! Reference altitude for launching altitude
+
+     REAL ran_num_1, ran_num_2, ran_num_3
+   
+     LOGICAL, save :: firstcall = .true.
+
+!--------------------------------------------------------------------------------------------------  
+! 1. INITIALISATION
+!------------------
+     IF (firstcall) THEN 
+        write(*,*) "nonoro_gwd_ran: FLott non-oro GW scheme is active!"
+        ! Characteristic vertical scale height
+        H0 = RD * tr / RG
+        ! 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
+
+! For VENUS, we use *_seri variables, that already integrate the different previous tendencies
+     tt(:,:) = pt(:,:)
+     uu(:,:) = pu(:,:)
+     vv(:,:) = pv(:,:)
+
+!    print*,'epflux_max just after firstcall:', epflux_max
+
+! 2. EVALUATION OF THE BACKGROUND FLOW AT SEMI-LEVELS
+!----------------------------------------------------
+     ! Pressure and Inv of pressure, temperature at 1/2 level
+     DO ll = 2, nlayer
+        ph(:, ll) = exp((log(pp(:, ll)) + log(pp(:, ll - 1))) / 2.)
+     end DO
+
+     ph(:, nlayer + 1) = 0.
+     ph(:, 1) = 2. * pp(:, 1) - ph(:, 2)
+   
+     ! Launching altitude for reproductible case
+     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
+     end DO
+    
+     ! Log-pressure vertical coordinate
+     DO ll = 1, nlayer + 1
+        zh(:,ll) = H0 * log(pr / (ph(:,ll) + psec))
+     end DO
+     
+     ! Winds and Brunt Vaisala frequency 
+     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 
+! VENUS ATTENTION: CP VARIABLE PSTAB CALCULE EN AMONT DES PARAMETRISATIONS
+        bv(:, ll) = max(bvsec,sqrt(pn2(:,ll)))
+     end DO
+
+     bv(:, 1) = bv(:, 2)
+     uh(:, 1) = 0.
+     vh(:, 1) = 0.
+     bv(:, nlayer + 1) = bv(:, nlayer)
+     uh(:, nlayer + 1) = uu(:, nlayer)
+     vh(:, nlayer + 1) = vv(:, nlayer)
+
+! 3. WAVES CHARACTERISTICS CHOSEN RANDOMLY
+!-----------------------------------------
+! The mod function of here a weird arguments 
+! are used to produce the waves characteristics 
+! in a stochastic way
+     DO jw = 1, nw 
+        DO ii = 1, ngrid
+           
+           ran_num_1 = mod(tt(ii, jw) * 10., 1.)
+           ran_num_2 = mod(tt(ii, jw) * 100., 1.)
+
+!! OPTIONS GENERIC DIFF VENUS !! 
+           ! angle (random)
+           zp(jw, ii) = (sign(1., 0.5 - ran_num_1) &
+                        + 1.) * RPI / 2.
+           ! angle (0 or RPI)
+           !zp(jw, ii) = RPI*mod(jw,2)
+
+           ! horizontal wavenumber amplitude
+!           zk(jw, ii) = kmin + (kmax - kmin) * ran_num_2
+           ! TN+GG April/June 2020
+           ! "Individual waves are not supposed
+           ! to occupy the entire domain, but
+           ! only a faction of it" Lott+2012
+           kstar = RPI / sqrt(cell_area(ii))
+           zk(jw, ii) = max(kmin,kstar) + (kmax - max(kmin,kstar)) * ran_num_2
+
+           ! horizontal phase speed
+! this computation allows to get a gaussian distribution centered on 0 (right ?)
+! then cmin is not useful, and it favors small values.
+! VCD 2.0 tuning
+           cpha = 0. 
+           DO jj = 1, na
+              ran_num_3 = mod(tt(ii, jw + 3 * jj)**2, 1.)
+              cpha = cpha + 2. * cmax *            &
+                     (ran_num_3 - 0.5) *           &
+                     sqrt(3.) / sqrt(na * 1.)
+           end DO
+           IF (cpha < 0.) THEN
+              cpha = - 1. * cpha
+              zp (jw, ii) = zp(jw, ii) + RPI
+           ENDIF
+! otherwise, with the computation below, we get a uniform distribution between cmin and cmax.
+!           ran_num_3 = mod(tt(ii, jw)**2, 1.)
+!           cpha = cmin + (cmax - cmin) * ran_num_3
+
+           ! 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, ii)) * uh(ii, launch) &
+                      + zk(jw, ii) * sin(zp(jw, ii)) * vh(ii, launch)
+
+           ! Momentum flux at launch level
+           epflux_0(jw, ii) = epflux_max                              &
+                        * mod(100. * (uu(ii, jw)**2 + vv(ii, jw)**2), 1.)
+
+        end DO
+     end DO
+!     print*,'epflux_0 just after waves charac. ramdon:', epflux_0
+
+! 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 altutide: 
+        intr_freq_p(jw, :) = zo(jw, :)                                &
+                   - zk(jw, :) * cos(zp(jw, :)) * uh(:, launch)       &
+                   - zk(jw, :) * sin(zp(jw, :)) * vh(:, launch)
+     end DO
+ 
+       ! Vertical velocity at launch level, value to ensure the imposed
+       ! mom flux:
+        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
+     
+     ! 4.2 Initial flux at launching altitude
+     !---------------------------------------
+     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 next done by:
+     !--------------------------------------------------------------------------
+     !     i) conserving the EP flux,
+     !     ii) dissipating a little, 
+     !     iii) testing critical levels, 
+     !     iv) testing the breaking.
+     !--------------------------------------------------------------------------
+     DO ll = launch, nlayer - 1
+        ! Warning! all the physics is here (passage from one level to the next
+        DO jw = 1, nw
+           intr_freq_m(jw, :) = intr_freq_p(jw, :)
+           wwm(jw, :) = wwp(jw, :)
+           ! Intrinsic frequency
+           intr_freq_p(jw, :) = zo(jw, :) - zk(jw, :) * cos(zp(jw, :)) * uh(:, ll + 1)    &
+                      - zk(jw, :) * sin(zp(jw, :)) * vh(:, ll + 1)
+           ! Vertical velocity in flux formulation
+           wwp(jw, :) = min(                                                              &
+                      ! No breaking (Eq. 6):
+                      wwm(jw, :)                                                          & 
+                      ! Dissipation (Eq. 8):
+                      * exp(-rdiss * pr / (ph(:, ll + 1) + ph (:, ll))                    &
+                      * ((bv(:, ll + 1) + bv (:, ll)) / 2.)**3                            &
+                      / 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., intr_freq_p(jw, :) * intr_freq_m(jw, :)))          &
+                      ! Saturation (Eq. 12)
+                      * abs(intr_freq_p(jw, :))**3 / bv(:, ll + 1)                        &
+                      * exp(-zh(:, ll + 1) / H0)                                          &
+                      * sat**2 * kmin**2 / zk(jw, :)**4)
+        end DO
+
+     ! Evaluate EP-flux from Eq. 7 and give the right orientation to the stress
+        DO jw = 1, nw
+           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
+
+        u_epflux_tot(:, ll + 1) = 0. 
+        v_epflux_tot(:, ll + 1) = 0.
+        DO jw = 1, nw
+           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, ll)) / float(nw)
+        end DO
+     end DO ! DO LL = LAUNCH, nlayer - 1
+
+!    print*,'u_epflux_tot just after launching:', u_epflux_tot
+!    print*,'v_epflux_tot just after launching:', v_epflux_tot
+!    print*,'u_epflux_p just after launching:', maxval(u_epflux_p), minval(u_epflux_p)
+!    print*,'v_epflux_p just after launching:', maxval(v_epflux_p), minval(v_epflux_p)
+
+! 5. TENDENCY CALCULATIONS
+!-------------------------
+     ! 5.1 Flow rectification at the top and in the low layers
+     ! --------------------------------------------------------
+     ! Warning, this is the total on all GW...
+
+     u_epflux_tot(:, nlayer + 1) = 0.
+     v_epflux_tot(:, nlayer + 1) = 0.
+
+     ! Here, big change compared to FLott version:
+     ! We compensate (u_epflux_(:, 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 
+
+     ! 5.2 AR-1 RECURSIVE FORMULA (13) IN VERSION 4
+     !---------------------------------------------
+     DO LL = 1, nlayer
+!       d_u(:, LL) = (u_epflux_tot(:, LL + 1) - u_epflux_tot(:, LL))
+       d_u(:, LL) = RG * (u_epflux_tot(:, LL + 1) - u_epflux_tot(:, LL)) &
+                  / (PH(:, LL + 1) - PH(:, LL)) * DTIME
+!       d_v(:, LL) = (v_epflux_tot(:, LL + 1) - v_epflux_tot(:, LL))
+       d_v(:, LL) = RG * (v_epflux_tot(:, LL + 1) - v_epflux_tot(:, LL)) &
+                  / (PH(:, LL + 1) - PH(:, LL)) * DTIME
+     ENDDO
+     d_t(:,:) = 0.
+
+     ! 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(:,:)
+
+!    print*,'u_epflux_tot just after tendency:', u_epflux_tot
+!    print*,'v_epflux_tot just after tendency:', v_epflux_tot
+!    print*,'d_u just after tendency:', maxval(d_u(:,:)), minval(d_u)
+!    print*,'d_v just after tendency:', maxval(d_v(:,:)), minval(d_v)
+
+     ! Cosmetic: evaluation of the cumulated stress
+     ZUSTR(:) = 0.
+     ZVSTR(:) = 0.
+     DO LL = 1, nlayer
+        ZUSTR(:) = ZUSTR(:) + D_U(:, LL) / RG * (PH(:, LL + 1) - PH(:, LL))
+        ZVSTR(:) = ZVSTR(:) + D_V(:, LL) / RG * (PH(:, LL + 1) - PH(:, LL))
+     ENDDO
+
+!#ifdef CPP_XIOS
+!     call send_xios_field("du_nonoro", d_u)
+!     call send_xios_field("dv_nonoro", d_v)
+!     call send_xios_field("east_gwstress", east_gwstress)
+!     call send_xios_field("west_gwstress", west_gwstress)
+!#endif
+
+      END SUBROUTINE nonoro_gwd_ran
+
+! ===================================================================
+! Subroutines used to write variables of memory in start files       
+! ===================================================================
+
+      SUBROUTINE ini_nonoro_gwd_ran(ngrid,nlayer)
+
+      IMPLICIT NONE
+
+      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))
+         allocate(east_gwstress(ngrid, nlayer))
+         allocate(west_gwstress(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)
+         if (allocated(east_gwstress)) deallocate(east_gwstress)
+         if (allocated(west_gwstress)) deallocate(west_gwstress)
+
+      END SUBROUTINE end_nonoro_gwd_ran
+
+END MODULE nonoro_gwd_ran_mod
+
+
Index: /trunk/LMDZ.VENUS/libf/phyvenus/photochemistry_venus.F90
===================================================================
--- /trunk/LMDZ.VENUS/libf/phyvenus/photochemistry_venus.F90	(revision 2579)
+++ /trunk/LMDZ.VENUS/libf/phyvenus/photochemistry_venus.F90	(revision 2580)
@@ -1,3 +1,3 @@
-subroutine photochemistry_venus(nz, n_lon, ptimestep, p, t, tr, mumean, sza_input, lon, lat, nesp, iter)
+subroutine photochemistry_venus(nz, n_lon, ptimestep, p, t, tr, mumean, sza_input, lon, lat, nesp, iter, prod_tr, loss_tr)
 
 use chemparam_mod
@@ -26,4 +26,6 @@
 
 real, dimension(nz,nesp) :: tr    ! tracer mixing ratio
+real, dimension(nz,nesp) :: prod_tr ! production (cm-3.s-1)
+real, dimension(nz,nesp) :: loss_tr ! loss       (cm-3.s-1)
 
 !===================================================================
@@ -86,5 +88,5 @@
 ! production and loss terms (for first-guess solution only)
 
-real, dimension(nesp) :: prod, loss
+real, dimension(nesp) :: prod, loss, lossconc
 
 ! indexes
@@ -170,5 +172,5 @@
 !  first-guess: fill matrix
 
-   call fill_matrix(iz, mat1, prod, loss, c, nesp, nz,                    &
+   call fill_matrix(iz, mat1, prod, loss, lossconc, c, nesp, nz,          &
                     nb_reaction_3_max, nb_reaction_4_max, nb_phot_max,    &
                     v_phot, v_3, v_4)
@@ -223,4 +225,9 @@
 
    tr(iz,:)  = max(c(iz,:)/conc(iz), 1.e-30)
+   
+!  save prod and loss espece
+   
+   prod_tr(iz,:) = prod(:)
+   loss_tr(iz,:) = lossconc(:)
                        
 end do  ! end of loop over vertical levels
@@ -1381,4 +1388,5 @@
 
 implicit none
+#include "clesphys.h"
 
 integer, INTENT(IN) :: nz
@@ -1547,4 +1555,31 @@
    v_phot(:,ij) = j(:,ij)
 end do
+
+!!! TEST: artificial increase of CO2 photodissociation
+if (tuneupperatm) then
+!-- TuneA
+!   v_phot(65:78,4) = v_phot(65:78,4)*10.
+!   v_phot(60:64,4) = v_phot(60:64,4)*3.
+!   v_phot(55:59,3) = v_phot(55:59,3)*2.
+!--
+!-- TuneB
+!   v_phot(65:78,4) = v_phot(65:78,4)*10.
+!   v_phot(55:59,3) = v_phot(55:59,3)*5.
+!--
+!-- TuneC
+! VCD 1.1 tuning
+    v_phot(65:78,4) = v_phot(65:78,4)*10.
+    v_phot(52:59,3) = v_phot(52:59,3)*5.
+!--
+!-- TuneE
+! VCD 2.0 tuning
+!   v_phot(65:90,4) = v_phot(65:90,4)*10.
+!--
+!   v_phot(:,4) = v_phot(:,4)*10.
+!do ij=3,4
+!   v_phot(:,ij) = v_phot(:,ij)*10.
+!end do
+endif
+!!!!!!!!!!!!!!!
 
 !PRINT*,'sza_input: ',sza_input
@@ -2862,5 +2897,5 @@
 !======================================================================
 
- subroutine fill_matrix(ilev, mat, prod, loss, c, nesp, nlayer,            &
+ subroutine fill_matrix(ilev, mat, prod, loss, lossconc, c, nesp, nlayer,  &
                         nb_reaction_3_max, nb_reaction_4_max, nb_phot_max, &
                         v_phot, v_3, v_4)
@@ -2894,5 +2929,5 @@
 
 real, dimension(nesp,nesp), intent(out) :: mat  ! matrix
-real, dimension(nesp), intent(out)      :: prod, loss
+real, dimension(nesp), intent(out)      :: prod, loss, lossconc
 
 ! local
@@ -2911,4 +2946,5 @@
 prod(:)  = 0.
 loss(:)  = 0.
+lossconc(:) = 0.
 
 ! photodissociations
@@ -2926,7 +2962,9 @@
   mat(ind_phot_6,ind_phot_2) = mat(ind_phot_6,ind_phot_2) - indice_phot(iphot)%z5*v_phot(ilev,iphot)
 
-  loss(ind_phot_2) = loss(ind_phot_2) + indice_phot(iphot)%z1*v_phot(ilev,iphot)
-  prod(ind_phot_4) = prod(ind_phot_4) + indice_phot(iphot)%z3*v_phot(ilev,iphot)*c(ilev,ind_phot_2)
-  prod(ind_phot_6) = prod(ind_phot_6) + indice_phot(iphot)%z5*v_phot(ilev,iphot)*c(ilev,ind_phot_2)
+  loss(ind_phot_2)     = loss(ind_phot_2) + indice_phot(iphot)%z1*v_phot(ilev,iphot)
+  lossconc(ind_phot_2) = lossconc(ind_phot_2) + indice_phot(iphot)%z1*v_phot(ilev,iphot)*c(ilev,ind_phot_2)
+  
+  prod(ind_phot_4)     = prod(ind_phot_4) + indice_phot(iphot)%z3*v_phot(ilev,iphot)*c(ilev,ind_phot_2)
+  prod(ind_phot_6)     = prod(ind_phot_6) + indice_phot(iphot)%z5*v_phot(ilev,iphot)*c(ilev,ind_phot_2)
 
 end do
@@ -2944,7 +2982,9 @@
   mat(ind_3_6,ind_3_2) = mat(ind_3_6,ind_3_2) - indice_3(i3)%z5*v_3(ilev,i3)*c(ilev,ind_3_2)
 
-  loss(ind_3_2) = loss(ind_3_2) + indice_3(i3)%z1*v_3(ilev,i3)*c(ilev,ind_3_2)
-  prod(ind_3_4) = prod(ind_3_4) + indice_3(i3)%z3*v_3(ilev,i3)*c(ilev,ind_3_2)*c(ilev,ind_3_2)
-  prod(ind_3_6) = prod(ind_3_6) + indice_3(i3)%z5*v_3(ilev,i3)*c(ilev,ind_3_2)*c(ilev,ind_3_2)
+  loss(ind_3_2)     = loss(ind_3_2) + indice_3(i3)%z1*v_3(ilev,i3)*c(ilev,ind_3_2)
+  lossconc(ind_3_2) = lossconc(ind_3_2) + indice_3(i3)%z1*v_3(ilev,i3)*c(ilev,ind_3_2)*c(ilev,ind_3_2)
+  
+  prod(ind_3_4)     = prod(ind_3_4) + indice_3(i3)%z3*v_3(ilev,i3)*c(ilev,ind_3_2)*c(ilev,ind_3_2)
+  prod(ind_3_6)     = prod(ind_3_6) + indice_3(i3)%z5*v_3(ilev,i3)*c(ilev,ind_3_2)*c(ilev,ind_3_2)
 
 end do
@@ -2974,8 +3014,11 @@
 
 
-  loss(ind_4_2) = loss(ind_4_2) + indice_4(i4)%z1*v_4(ilev,i4)*c(ilev,ind_4_4)
-  loss(ind_4_4) = loss(ind_4_4) + indice_4(i4)%z3*v_4(ilev,i4)*c(ilev,ind_4_2)
-  prod(ind_4_6) = prod(ind_4_6) + indice_4(i4)%z5*v_4(ilev,i4)*c(ilev,ind_4_2)*c(ilev,ind_4_4)
-  prod(ind_4_8) = prod(ind_4_8) + indice_4(i4)%z7*v_4(ilev,i4)*c(ilev,ind_4_2)*c(ilev,ind_4_4)
+  loss(ind_4_2)     = loss(ind_4_2) + indice_4(i4)%z1*v_4(ilev,i4)*c(ilev,ind_4_4)
+  lossconc(ind_4_2) = lossconc(ind_4_2) + indice_4(i4)%z1*v_4(ilev,i4)*c(ilev,ind_4_4)*c(ilev,ind_4_2)
+  loss(ind_4_4)     = loss(ind_4_4) + indice_4(i4)%z3*v_4(ilev,i4)*c(ilev,ind_4_2)
+  lossconc(ind_4_4) = lossconc(ind_4_4) + indice_4(i4)%z3*v_4(ilev,i4)*c(ilev,ind_4_2)*c(ilev,ind_4_4)
+  
+  prod(ind_4_6)     = prod(ind_4_6) + indice_4(i4)%z5*v_4(ilev,i4)*c(ilev,ind_4_2)*c(ilev,ind_4_4)
+  prod(ind_4_8)     = prod(ind_4_8) + indice_4(i4)%z7*v_4(ilev,i4)*c(ilev,ind_4_2)*c(ilev,ind_4_4)
 
 end do
Index: /trunk/LMDZ.VENUS/libf/phyvenus/phyetat0.F90
===================================================================
--- /trunk/LMDZ.VENUS/libf/phyvenus/phyetat0.F90	(revision 2579)
+++ /trunk/LMDZ.VENUS/libf/phyvenus/phyetat0.F90	(revision 2580)
@@ -15,4 +15,6 @@
       USE time_phylmdz_mod, only: itau_phy, raz_date, pdtphys
       USE ioipsl_getin_p_mod, only: getin_p
+  use nonoro_gwd_ran_mod, only: du_nonoro_gwd, dv_nonoro_gwd, &
+                                east_gwstress, west_gwstress
 
 implicit none
@@ -365,4 +367,34 @@
 PRINT*,'Turbulent Kinetic Energy', minval(q2), maxval(q2)
 
+! Non-orographic gravity waves 
+if (startphy_file) then
+   call get_field("du_nonoro_gwd",du_nonoro_gwd,found)
+   if (.not.found) then
+      write(*,*) "phyetat0: <du_nonoro_gwd> not in file"
+      du_nonoro_gwd(:,:)=0.
+   endif
+endif ! of if (startphy_file)
+if (startphy_file) then
+   call get_field("dv_nonoro_gwd",dv_nonoro_gwd,found)
+   if (.not.found) then
+      write(*,*) "phyetat0: <dv_nonoro_gwd> not in file"
+      dv_nonoro_gwd(:,:)=0.
+   endif
+endif ! of if (startphy_file)
+if (startphy_file) then
+   call get_field("east_gwstress",east_gwstress,found)
+   if (.not.found) then
+      write(*,*) "phyetat0: <east_gwstress> not in file"
+      east_gwstress(:,:)=0.
+   endif
+endif ! of if (startphy_file)
+if (startphy_file) then
+   call get_field("west_gwstress",west_gwstress,found)
+   if (.not.found) then
+      write(*,*) "phyetat0: <west_gwstress> not in file"
+      west_gwstress(:,:)=0.
+   endif
+endif ! of if (startphy_file)
+
 ! close file
 IF (startphy_file) call close_startphy
Index: /trunk/LMDZ.VENUS/libf/phyvenus/phyredem.F90
===================================================================
--- /trunk/LMDZ.VENUS/libf/phyvenus/phyredem.F90	(revision 2579)
+++ /trunk/LMDZ.VENUS/libf/phyvenus/phyredem.F90	(revision 2580)
@@ -13,4 +13,6 @@
       use geometry_mod, only: longitude_deg, latitude_deg
       USE time_phylmdz_mod, only: day_end, annee_ref, itau_phy, raz_date
+  use nonoro_gwd_ran_mod, only: du_nonoro_gwd, dv_nonoro_gwd, &
+                                east_gwstress, west_gwstress
 
       implicit none
@@ -99,4 +101,12 @@
       CALL put_field("TANCIEN","T Previous iteration",t_ancien)
 
+! Non-orographic gavity waves
+if (ok_gw_nonoro) then
+   call put_field("du_nonoro_gwd","Zonal wind tendency due to GW",du_nonoro_gwd)
+   call put_field("dv_nonoro_gwd","Meridional wind tendency due to GW",dv_nonoro_gwd)
+   call put_field("east_gwstress","Eastward stress profile due to GW",east_gwstress)
+   call put_field("west_gwstress","Westward stress profile due to GW",west_gwstress)
+endif
+
 ! close file
 
Index: /trunk/LMDZ.VENUS/libf/phyvenus/phys_state_var_mod.F90
===================================================================
--- /trunk/LMDZ.VENUS/libf/phyvenus/phys_state_var_mod.F90	(revision 2579)
+++ /trunk/LMDZ.VENUS/libf/phyvenus/phys_state_var_mod.F90	(revision 2580)
@@ -11,4 +11,5 @@
       USE dimphy
       USE turb_mod
+      use nonoro_gwd_ran_mod, only: ini_nonoro_gwd_ran, end_nonoro_gwd_ran 
 !      INTEGER, SAVE :: radpas
 !!$OMP THREADPRIVATE(radpas)
@@ -173,4 +174,6 @@
       ! allocate turb_mod module variables
       CALL ini_turb_mod(klon,klev)
+      ! allocate arrays in "nonoro_gwd_ran_mod"
+      call ini_nonoro_gwd_ran(klon,klev)
       
 END SUBROUTINE phys_state_var_init
@@ -212,4 +215,6 @@
       ! deallocate turb_mod module variables
       CALL end_turb_mod
+      ! deallocate arrays in "nonoro_gwd_ran_mod"
+      call end_nonoro_gwd_ran
 
 END SUBROUTINE phys_state_var_end
Index: /trunk/LMDZ.VENUS/libf/phyvenus/physiq_mod.F
===================================================================
--- /trunk/LMDZ.VENUS/libf/phyvenus/physiq_mod.F	(revision 2579)
+++ /trunk/LMDZ.VENUS/libf/phyvenus/physiq_mod.F	(revision 2580)
@@ -78,4 +78,5 @@
       use vertical_layers_mod, only: pseudoalt
       use turb_mod, only : sens, turb_resolved
+      use nonoro_gwd_ran_mod, only: nonoro_gwd_ran
       use sed_and_prod_mad, only: aer_sedimentation, drop_sedimentation
 #ifdef CPP_XIOS      
@@ -362,4 +363,6 @@
 c pour sorties
       REAL :: col_dens_tr(klon,nqmax)
+      REAL,allocatable,save :: prod_tr(:,:,:)
+      REAL,allocatable,save :: loss_tr(:,:,:)
 
 c pour ioipsl
@@ -402,4 +405,6 @@
 c-jld ec_conser
 
+c ALBEDO VARIATIONS (VCD)
+      REAL factAlb
 c TEST VENUS...
       REAL mang(klon,klev)    ! moment cinetique
@@ -448,4 +453,6 @@
       IF (debut) THEN
          allocate(source(klon,nqmax))
+         allocate(prod_tr(klon,klev,nqmax))
+         allocate(loss_tr(klon,klev,nqmax))
 
 #ifdef CPP_XIOS
@@ -1067,5 +1074,7 @@
      $                       tr_seri,
      $                       d_tr_chem,
-     $                       iter)
+     $                       iter,
+     $                       prod_tr,
+     $                       loss_tr)
 
          if (ok_sedim) then 
@@ -1514,9 +1523,26 @@
       endif
 
-c albedo variations: test for Yeon Joo Lee
+c ALBEDO VARIATIONS: test for Yeon Joo Lee
 c  increment to increase it for 20 Vd => +80% 
 c       heat(:,:)=heat(:,:)*(1.+0.80*((rjourvrai-356)+gmtime)/20.)
 c  or to decrease it for 20 Vd => 1/1.8 
 c       heat(:,:)=heat(:,:)/(1.+0.80*((rjourvrai-356)+gmtime)/20.)
+
+c ------------ ALBEDO VARIATIONS: scenarios for VCD
+c shape of relative variation from Lee et al 2019 (Fig 13b)
+c between 57 km (4e4 Pa) and 72 km (2.5e3 Pa), peak at 67 km (6e3 Pa)
+c      do j=1,klev
+c       factAlb = 0.
+c       if ((presnivs(j).gt.6.e3).and.(presnivs(j).lt.4.e4)) then
+c        factAlb = (log(presnivs(j))-log(4.e4))/(log(6.e3)-log(4.e4))
+c       elseif ((presnivs(j).lt.6.e3).and.(presnivs(j).gt.2.5e3)) then
+c        factAlb = (log(presnivs(j))-log(2.5e3))/(log(6.e3)-log(2.5e3))
+c       endif
+c Increase by 50% (Minimum albedo)
+c       heat(:,j)=heat(:,j)*(1+factAlb*0.5)
+c Decrease by 30% (Maximum albedo)
+c       heat(:,j)=heat(:,j)*(1-factAlb*0.3)
+c      enddo
+c ------------ END ALBEDO VARIATIONS 
 
 cc---------------------------------------------
@@ -1784,8 +1810,17 @@
        IF(ok_gw_nonoro) then
 
+! Obsolete
+! but used for VCD 1.1
       call flott_gwd_ran(klon,klev,dtime,pplay,zn2,
      e               t_seri, u_seri, v_seri, paprs(klon/2+1,:),
      o               zustrhi,zvstrhi,
      o               d_t_hin, d_u_hin, d_v_hin)
+
+! New routine based on Generic
+! used after VCD 1.1
+!     call nonoro_gwd_ran(klon,klev,dtime,pplay,zn2,presnivs,
+!    e               t_seri, u_seri, v_seri,
+!    o               zustrhi,zvstrhi,
+!    o               d_t_hin, d_u_hin, d_v_hin)
 
 c  ajout des tendances
@@ -2065,4 +2100,23 @@
 
       if (iflag_trac == 1) then
+
+! production and destruction rate, cm-3.s-1
+! Beware of the context*.xml file !!
+!         if ((tr_scheme == 3) .and. (ok_chem)) THEN
+!            do iq = 1,nqmax - nmicro
+!               !if ((iq == i_co) .or. (iq == i_o)) then
+!               !   call send_xios_field('prod_'+tname(iq), prod_tr(:,:,iq))
+!               !   call send_xios_field('loss_'+tname(iq), loss_tr(:,:,iq))
+!               !end if
+!               if (iq == i_o) then
+!                  call send_xios_field('prod_o', prod_tr(:,:,iq))
+!                  call send_xios_field('loss_o', loss_tr(:,:,iq))
+!               end if
+!               if (iq == i_co) then
+!                  call send_xios_field('prod_co', prod_tr(:,:,iq))
+!                  call send_xios_field('loss_co', loss_tr(:,:,iq))
+!               end if
+!            end do
+!         end if
 
 ! tracers in gas phase, volume mixing ratio
Index: /trunk/LMDZ.VENUS/libf/phyvenus/phytrac_chimie.F
===================================================================
--- /trunk/LMDZ.VENUS/libf/phyvenus/phytrac_chimie.F	(revision 2579)
+++ /trunk/LMDZ.VENUS/libf/phyvenus/phytrac_chimie.F	(revision 2580)
@@ -12,8 +12,13 @@
      $                    trac,
      $                    d_tr_chem,
-     $                    iter)
+     $                    iter,
+     $                    prod_tr,
+     $                    loss_tr)
 
       use chemparam_mod
       use conc, only: mmean,rnew
+#ifdef CPP_XIOS      
+      use xios_output_mod, only: send_xios_field
+#endif
 
       implicit none
@@ -43,4 +48,5 @@
       real, dimension(nlon,nlev,nqmax) :: d_tr_chem  ! chemical tendency for each tracer
       integer, dimension(nlon,nlev) :: iter          ! chemical iterations
+      real, dimension(nlon,nlev,nqmax) :: prod_tr, loss_tr  ! produc and loss tracer
 
 !===================================================================
@@ -61,5 +67,5 @@
       real, dimension(nlon,nlev) :: trac_sum
       real, dimension(nlon,nlev,nqmax) :: ztrac  ! local tracer mixing ratio
-     
+      
 !===================================================================
 !     first call : initialisations
@@ -271,10 +277,11 @@
      $                                ztrac(ilon,:,:),
      $                                mmean(ilon,:),
-     $                                sza_local, 
+     $                                sza_local,
      $                                lon(ilon), lat(ilon),
-     $                                nqmax, iter(ilon,:))
-
-         end do
-
+     $                                nqmax, iter(ilon,:),
+     $                                prod_tr(ilon,:,:),
+     $                                loss_tr(ilon,:,:))
+
+         end do         
       end if  ! ok_chem
 
