Index: trunk/LMDZ.VENUS/libf/phyvenus/aeropacity.F90
===================================================================
--- trunk/LMDZ.VENUS/libf/phyvenus/aeropacity.F90	(revision 2580)
+++ trunk/LMDZ.VENUS/libf/phyvenus/aeropacity.F90	(revision 2622)
@@ -2,7 +2,7 @@
          aerosol,reffrad,nueffrad,QREFvis3d,tau_col)
 
-       use radinc_h, only : L_TAUMAX,naerkind
+       use radinc_h, only : naerkind
        use aerosol_mod
-                  
+                   
        implicit none
 #include "YOMCST.h"
@@ -49,365 +49,142 @@
       REAL,INTENT(OUT):: tau_col(ngrid) !column integrated visible optical depth
 
-      real aerosol0, obs_tau_col_aurora, pm
-      real pcloud_deck, cloud_slope
-
-      real dp_strato(ngrid)
-      real dp_tropo(ngrid)
-      real dp_layer(ngrid)
-
-      INTEGER l,ig,iq,iaer,ia,ll
+      INTEGER l,ig,iaer,ll
 
       LOGICAL,SAVE :: firstcall=.true.
 !$OMP THREADPRIVATE(firstcall)
-      REAL CBRT
-      EXTERNAL CBRT
-
-      ! for fixed dust profiles
-      real topdust, expfactor, zp
-      REAL taudusttmp(ngrid) ! Temporary dust opacity used before scaling
-      REAL tauh2so4tmp(ngrid) ! Temporary h2so4 opacity used before scaling
-
-      real CLFtot
+      character*80 abort_message
+      character*20 modname
 
       !  for venus clouds
-      real :: p_bot2,p_bot,p_top,p_top2,h_bot2,h_bot,h_top,h_top2
-      real :: mode_dens,h_lay,nmode1,nmode2,nmode2p,nmode3,nmodeuv
-
+      real,ALLOCATABLE,SAVE :: pbot2(:,:),pbot(:,:),ptop(:,:),ptop2(:,:)
+      real,ALLOCATABLE,SAVE :: hbot2(:,:),hbot(:,:),htop(:,:),htop2(:,:)
+      real,ALLOCATABLE,SAVE :: nmode(:,:)
+!$OMP THREADPRIVATE(pbot2,pbot,ptop,ptop2)
+!$OMP THREADPRIVATE(hbot2,hbot,htop,htop2)
+!$OMP THREADPRIVATE(nmode)
+      real :: mode_dens,hlay
+
+! ---------------------------------------------------------
+! ---------------------------------------------------------
 ! First call only
       IF (firstcall) THEN
-
-      ! identify tracers
-        write(*,*) "Active aerosols found in aeropacity:"
-
-        if (iaero_venus1.ne.0) then
+      
+        modname = 'aeropacity'
+
+! verify tracers
+        write(*,*) "Active aerosols found in aeropacity, designed for Venus"
+
+        if (iaero_venus1.eq.1) then
           print*,'iaero_venus1= ',iaero_venus1
-        endif
-        if (iaero_venus2.ne.0) then
+        else
+            abort_message='iaero_venus1 is not 1...'
+            call abort_physic(modname,abort_message,1)
+        endif
+        if (iaero_venus2.eq.2) then
           print*,'iaero_venus2= ',iaero_venus2
-        endif
-        if (iaero_venus2p.ne.0) then
+        else
+            abort_message='iaero_venus2 is not 2...'
+            call abort_physic(modname,abort_message,1)
+        endif
+        if (iaero_venus2p.eq.3) then
           print*,'iaero_venus2p= ',iaero_venus2p
-        endif
-        if (iaero_venus3.ne.0) then
+        else
+            abort_message='iaero_venus2p is not 3...'
+            call abort_physic(modname,abort_message,1)
+        endif
+        if (iaero_venus3.eq.4) then
           print*,'iaero_venus3= ',iaero_venus3
-        endif
-        if (iaero_venusUV.ne.0) then
+        else
+            abort_message='iaero_venus3 is not 4...'
+            call abort_physic(modname,abort_message,1)
+        endif
+        if (iaero_venusUV.eq.5) then
           print*,'iaero_venusUV= ',iaero_venusUV
-        endif
+        else
+            abort_message='iaero_venus1UV is not 5...'
+            call abort_physic(modname,abort_message,1)
+        endif
+
+! cloud model
+        allocate(pbot2(ngrid,naerkind),pbot(ngrid,naerkind))
+        allocate(ptop(ngrid,naerkind),ptop2(ngrid,naerkind))
+        allocate(hbot2(ngrid,naerkind),hbot(ngrid,naerkind))
+        allocate(htop(ngrid,naerkind),htop2(ngrid,naerkind))
+        allocate(nmode(ngrid,naerkind))
+        call cloud_haus16(ngrid,pbot2,pbot,hbot2,hbot,ptop,ptop2,htop,htop2,nmode)
 
         firstcall=.false.
       ENDIF ! of IF (firstcall)
-
-
-!     ---------------------------------------------------------
-
-!==================================================================
-!         Venus clouds (4 modes)
-!   S. Lebonnois (jan 2016)
-!==================================================================
-! distributions from Haus et al, 2016
-! mode             1      2      2p     3
-! r (microns)     0.30   1.05   1.40   3.65
-! sigma           1.56   1.29   1.23   1.28
-! reff (microns)  0.49   1.23   1.56   4.25
-! nueff           0.21   0.067  0.044  0.063
-! (nueff=exp(ln^2 sigma)-1)
-!
-! p_bot <=> zb ; p_top <=> zb+zc ; h_bot <=> Hlo ; h_top <=> Hup
-! p<p_top: N(i+1)=N(i)*(p(i+1)/p(i))**(h_lay(i)/h_top)      
-!    h_lay=R(T(i)+T(i+1))/2/g  (in m)
-!    R=8.31/mu (mu in kg/mol)
-! p>p_bot: N(i-1)=N(i)*(p(i)/p(i-1))**(h_lay(i)/h_bot)      
-! N is in m-3
-!
-! dTau = Qext*[pi*reff**2*exp(-3*ln(1+nueff))]*N*h_lay*(-dp)/p
-
-! Mode 1
-      if (iaero_venus1 .ne.0) then
-          iaer=iaero_venus1
-
-!       1. Initialization
-          aerosol(1:ngrid,1:nlayer,iaer)=0.0
-! two scales heights for below and above max density (cf Haus et al 16)
-          p_bot2 = 2.0e5 ! Pa   2.0e5
-            h_bot2 = 5.0e3 ! m
-          p_bot  = 1.2e5   !    1.2e5
-            h_bot  = 1.0e3
-          nmode1 = 1.935e8 ! m-3
-          p_top  = 1.0e4
-            h_top  = 3.5e3
-          p_top2 = 4.5e2
-            h_top2 = 2.0e3
-
-!       2. Opacity calculation
-
-          DO ig=1,ngrid
+! ---------------------------------------------------------
+! ---------------------------------------------------------
+
+    aerosol(:,:,:) = 0.
+        tau_col(:) = 0.
+
+      do ig=1,ngrid
+        do iaer=1,naerkind
+
+!       Opacity calculation
+
 ! determine the approximate middle of the mode layer
            ll=1
            DO l=1,nlayer-1    ! high p to low p
-             IF (pplay(ig,l) .gt. (p_top+p_bot)/2) ll=l
+             IF (pplay(ig,l) .gt. (ptop(ig,iaer)+pbot(ig,iaer))/2) ll=l
            ENDDO
-! from there go down and up for profile N
-           mode_dens = nmode1  ! m-3
+! from there go DOWN for profile N
+           mode_dens = nmode(ig,iaer)  ! m-3
            DO l=ll,1,-1
-             h_lay=RD*(pt(ig,l+1)+pt(ig,l))/2./RG
-             IF (pplay(ig,l) .le. p_bot) THEN
-               mode_dens = nmode1  ! m-3
-             ELSEIF (pplay(ig,l) .gt. p_bot .and. pplay(ig,l) .le. p_bot2) THEN
-               mode_dens = mode_dens*(pplay(ig,l+1)/pplay(ig,l))**(h_lay/h_bot)
-             ELSEIF (pplay(ig,l) .gt. p_bot2) THEN
-               mode_dens = mode_dens*(pplay(ig,l+1)/pplay(ig,l))**(h_lay/h_bot2)
+             hlay=RD*(pt(ig,l+1)+pt(ig,l))/2./RG
+             IF (pplay(ig,l) .le. pbot(ig,iaer)) THEN
+               mode_dens = nmode(ig,iaer)  ! m-3
+             ELSEIF (pplay(ig,l) .gt. pbot(ig,iaer) .and. pplay(ig,l) .le. pbot2(ig,iaer)) THEN
+               mode_dens = mode_dens*(pplay(ig,l+1)/pplay(ig,l))**(hlay/hbot(ig,iaer))
+             ELSEIF (pplay(ig,l) .gt. pbot2(ig,iaer)) THEN
+               mode_dens = mode_dens*(pplay(ig,l+1)/pplay(ig,l))**(hlay/hbot2(ig,iaer))
              ENDIF
              mode_dens=max(1.e-30,mode_dens)
-             aerosol(ig,l,iaer) = QREFvis3d(ig,l,iaer)*                       &
-              RPI*(reffrad(ig,l,iaer))**2*exp(-3.*log(1+nueffrad(ig,l,iaer)))* & 
-              mode_dens*h_lay*(pplev(ig,l)-pplev(ig,l+1))/pplay(ig,l)
-!             aerosol(ig,l,iaer) = mode_dens
+             if (iaer.lt.5) then
+              aerosol(ig,l,iaer) = QREFvis3d(ig,l,iaer)*                       &
+               RPI*(reffrad(ig,l,iaer))**2*exp(-3.*log(1+nueffrad(ig,l,iaer)))* & 
+               mode_dens*hlay*(pplev(ig,l)-pplev(ig,l+1))/pplay(ig,l)
+             else  ! for UV absorber
+! normalized to 0.35 microns (peak of absorption)
+              aerosol(ig,l,iaer) = QREFvis3d(ig,l,iaer)*mode_dens
+             endif
            ENDDO
-           mode_dens = nmode1  ! m-3
+! ... then UP for profile N
+           mode_dens = nmode(ig,iaer)  ! m-3
            DO l=ll+1,nlayer-1
-             h_lay=RD*(pt(ig,l-1)+pt(ig,l))/2./RG
-             IF (pplay(ig,l) .gt. p_top) THEN
-               mode_dens = nmode1  ! m-3
-             ELSEIF (pplay(ig,l) .le. p_top .and. pplay(ig,l) .gt. p_top2) THEN
-               mode_dens = mode_dens*(pplay(ig,l)/pplay(ig,l-1))**(h_lay/h_top)
-             ELSEIF (pplay(ig,l) .le. p_top2) THEN
-               mode_dens = mode_dens*(pplay(ig,l)/pplay(ig,l-1))**(h_lay/h_top2)
+             hlay=RD*(pt(ig,l-1)+pt(ig,l))/2./RG
+             IF (pplay(ig,l) .gt. ptop(ig,iaer)) THEN
+               mode_dens = nmode(ig,iaer)  ! m-3
+             ELSEIF (pplay(ig,l) .le. ptop(ig,iaer) .and. pplay(ig,l) .gt. ptop2(ig,iaer)) THEN
+               mode_dens = mode_dens*(pplay(ig,l)/pplay(ig,l-1))**(hlay/htop(ig,iaer))
+             ELSEIF (pplay(ig,l) .le. ptop2(ig,iaer)) THEN
+               mode_dens = mode_dens*(pplay(ig,l)/pplay(ig,l-1))**(hlay/htop2(ig,iaer))
              ENDIF
              mode_dens=max(1.e-30,mode_dens)
-             aerosol(ig,l,iaer) = QREFvis3d(ig,l,iaer)*                       &
-              RPI*(reffrad(ig,l,iaer))**2*exp(-3.*log(1+nueffrad(ig,l,iaer)))* & 
-              mode_dens*h_lay*(pplev(ig,l)-pplev(ig,l+1))/pplay(ig,l)
-!             aerosol(ig,l,iaer) = mode_dens
+             if (iaer.lt.5) then
+              aerosol(ig,l,iaer) = QREFvis3d(ig,l,iaer)*                       &
+               RPI*(reffrad(ig,l,iaer))**2*exp(-3.*log(1+nueffrad(ig,l,iaer)))* & 
+               mode_dens*hlay*(pplev(ig,l)-pplev(ig,l+1))/pplay(ig,l)
+             else  ! for UV absorber
+! normalized to 0.35 microns (peak of absorption)
+              aerosol(ig,l,iaer) = QREFvis3d(ig,l,iaer)*mode_dens
+             endif
            ENDDO
-          ENDDO
-
-      end if ! mode 1
-
-! Mode 2
-      if (iaero_venus2 .ne.0) then
-          iaer=iaero_venus2
-
-!       1. Initialization
-          aerosol(1:ngrid,1:nlayer,iaer)=0.0
-          p_bot = 1.0e4
-            h_bot = 3.0e3
-          nmode2 = 1.00e8 ! m-3
-          p_top = 8.0e3
-            h_top = 3.5e3
-          p_top2 = 4.5e2
-            h_top2 = 2.0e3
-
-!       2. Opacity calculation
-
-          DO ig=1,ngrid
-! determine the approximate middle of the mode layer
-           ll=1
-           DO l=1,nlayer-1    ! high p to low p
-             IF (pplay(ig,l) .gt. (p_top+p_bot)/2) ll=l
-           ENDDO
-! from there go down and up for profile N
-           mode_dens = nmode2  ! m-3
-           DO l=ll,1,-1
-             h_lay=RD*(pt(ig,l+1)+pt(ig,l))/2./RG
-             IF (pplay(ig,l) .le. p_bot) THEN
-               mode_dens = nmode2  ! m-3
-             ELSEIF (pplay(ig,l) .gt. p_bot) THEN
-               mode_dens = mode_dens*(pplay(ig,l+1)/pplay(ig,l))**(h_lay/h_bot)
-             ENDIF
-             mode_dens=max(1.e-30,mode_dens)
-             aerosol(ig,l,iaer) = QREFvis3d(ig,l,iaer)*                       &
-              RPI*(reffrad(ig,l,iaer))**2*exp(-3.*log(1+nueffrad(ig,l,iaer)))* & 
-              mode_dens*h_lay*(pplev(ig,l)-pplev(ig,l+1))/pplay(ig,l)
-!             aerosol(ig,l,iaer) = mode_dens
-           ENDDO
-           mode_dens = nmode2  ! m-3
-           DO l=ll+1,nlayer-1
-             h_lay=RD*(pt(ig,l-1)+pt(ig,l))/2./RG
-             IF (pplay(ig,l) .gt. p_top) THEN
-               mode_dens = nmode2  ! m-3
-             ELSEIF (pplay(ig,l) .le. p_top .and. pplay(ig,l) .gt. p_top2) THEN
-               mode_dens = mode_dens*(pplay(ig,l)/pplay(ig,l-1))**(h_lay/h_top)
-             ELSEIF (pplay(ig,l) .le. p_top2) THEN
-               mode_dens = mode_dens*(pplay(ig,l)/pplay(ig,l-1))**(h_lay/h_top2)
-             ENDIF
-             mode_dens=max(1.e-30,mode_dens)
-             aerosol(ig,l,iaer) = QREFvis3d(ig,l,iaer)*                       &
-              RPI*(reffrad(ig,l,iaer))**2*exp(-3.*log(1+nueffrad(ig,l,iaer)))* & 
-              mode_dens*h_lay*(pplev(ig,l)-pplev(ig,l+1))/pplay(ig,l)
-!             aerosol(ig,l,iaer) = mode_dens
-           ENDDO
-          ENDDO
-
-      end if ! mode 2
-
-! Mode 2p
-      if (iaero_venus2p .ne.0) then
-          iaer=iaero_venus2p
-
-!       1. Initialization
-          aerosol(1:ngrid,1:nlayer,iaer)=0.0
-          p_bot = 1.2e5  ! 1.2e5
-            h_bot = 0.1e3
-          nmode2p = 5.0e7 ! m-3
-          p_top = 2.3e4
-            h_top = 1.0e3
-
-!       2. Opacity calculation
-
-          DO ig=1,ngrid
-! determine the approximate middle of the mode layer
-           ll=1
-           DO l=1,nlayer-1    ! high p to low p
-             IF (pplay(ig,l) .gt. (p_top+p_bot)/2) ll=l
-           ENDDO
-! from there go down and up for profile N
-           mode_dens = nmode2p  ! m-3
-           DO l=ll,1,-1
-             h_lay=RD*(pt(ig,l+1)+pt(ig,l))/2./RG
-             IF (pplay(ig,l) .le. p_bot) THEN
-               mode_dens = nmode2p  ! m-3
-             ELSEIF (pplay(ig,l) .gt. p_bot) THEN
-               mode_dens = mode_dens*(pplay(ig,l+1)/pplay(ig,l))**(h_lay/h_bot)
-             ENDIF
-             mode_dens=max(1.e-30,mode_dens)
-             aerosol(ig,l,iaer) = QREFvis3d(ig,l,iaer)*                       &
-              RPI*(reffrad(ig,l,iaer))**2*exp(-3.*log(1+nueffrad(ig,l,iaer)))* & 
-              mode_dens*h_lay*(pplev(ig,l)-pplev(ig,l+1))/pplay(ig,l)
-!             aerosol(ig,l,iaer) = mode_dens
-           ENDDO
-           mode_dens = nmode2p  ! m-3
-           DO l=ll+1,nlayer-1
-             h_lay=RD*(pt(ig,l-1)+pt(ig,l))/2./RG
-             IF (pplay(ig,l) .gt. p_top) THEN
-               mode_dens = nmode2p  ! m-3
-             ELSEIF (pplay(ig,l) .le. p_top) THEN
-               mode_dens = mode_dens*(pplay(ig,l)/pplay(ig,l-1))**(h_lay/h_top)
-             ENDIF
-             mode_dens=max(1.e-30,mode_dens)
-             aerosol(ig,l,iaer) = QREFvis3d(ig,l,iaer)*                       &
-              RPI*(reffrad(ig,l,iaer))**2*exp(-3.*log(1+nueffrad(ig,l,iaer)))* & 
-              mode_dens*h_lay*(pplev(ig,l)-pplev(ig,l+1))/pplay(ig,l)
-!             aerosol(ig,l,iaer) = mode_dens
-           ENDDO
-          ENDDO
-
-      end if ! mode 2p
-
-! Mode 3
-      if (iaero_venus3 .ne.0) then
-          iaer=iaero_venus3
-
-!       1. Initialization
-          aerosol(1:ngrid,1:nlayer,iaer)=0.0
-          p_bot = 1.2e5  ! 1.2e5
-            h_bot = 0.5e3
-          nmode3 = 1.4e7 ! m-3
-          p_top = 3.9e4
-            h_top = 1.0e3
-
-!       2. Opacity calculation
-
-          DO ig=1,ngrid
-! determine the approximate middle of the mode layer
-           ll=1
-           DO l=1,nlayer-1    ! high p to low p
-             IF (pplay(ig,l) .gt. (p_top+p_bot)/2) ll=l
-           ENDDO
-! from there go down and up for profile N
-           mode_dens = nmode3  ! m-3
-           DO l=ll,1,-1
-             h_lay=RD*(pt(ig,l+1)+pt(ig,l))/2./RG
-             IF (pplay(ig,l) .le. p_bot) THEN
-               mode_dens = nmode3  ! m-3
-             ELSEIF (pplay(ig,l) .gt. p_bot) THEN
-               mode_dens = mode_dens*(pplay(ig,l+1)/pplay(ig,l))**(h_lay/h_bot)
-             ENDIF
-             mode_dens=max(1.e-30,mode_dens)
-             aerosol(ig,l,iaer) = QREFvis3d(ig,l,iaer)*                       &
-              RPI*(reffrad(ig,l,iaer))**2*exp(-3.*log(1+nueffrad(ig,l,iaer)))* & 
-              mode_dens*h_lay*(pplev(ig,l)-pplev(ig,l+1))/pplay(ig,l)
-!             aerosol(ig,l,iaer) = mode_dens
-           ENDDO
-           mode_dens = nmode3  ! m-3
-           DO l=ll+1,nlayer-1
-             h_lay=RD*(pt(ig,l-1)+pt(ig,l))/2./RG
-             IF (pplay(ig,l) .gt. p_top) THEN
-               mode_dens = nmode3  ! m-3
-             ELSEIF (pplay(ig,l) .le. p_top) THEN
-               mode_dens = mode_dens*(pplay(ig,l)/pplay(ig,l-1))**(h_lay/h_top)
-             ENDIF
-             mode_dens=max(1.e-30,mode_dens)
-             aerosol(ig,l,iaer) = QREFvis3d(ig,l,iaer)*                       &
-              RPI*(reffrad(ig,l,iaer))**2*exp(-3.*log(1+nueffrad(ig,l,iaer)))* & 
-              mode_dens*h_lay*(pplev(ig,l)-pplev(ig,l+1))/pplay(ig,l)
-!             aerosol(ig,l,iaer) = mode_dens
-           ENDDO
-          ENDDO
-
-      end if ! mode 3 
-
-! UV absorber
-      if (iaero_venusUV .ne.0) then
-          iaer=iaero_venusUV
-
-!       1. Initialization
-          aerosol(1:ngrid,1:nlayer,iaer)=0.0
-          p_bot = 3.3e4  ! 58 km
-            h_bot = 1.0e3
-          nmodeuv = 1.00e7 !m-3
-          p_top = 3.7e3 ! 70 km
-            h_top = 1.0e3
-
-!       2. Opacity calculation
-
-          DO ig=1,ngrid
-! determine the approximate middle of the mode layer
-           ll=1
-           DO l=1,nlayer-1    ! high p to low p
-             IF (pplay(ig,l) .gt. (p_top+p_bot)/2) ll=l
-           ENDDO
-! from there go down and up for profile N
-           mode_dens = nmodeuv ! m-3
-           DO l=ll,1,-1
-             h_lay=RD*(pt(ig,l+1)+pt(ig,l))/2./RG
-             IF (pplay(ig,l) .le. p_bot) THEN
-               mode_dens = nmodeuv ! m-3
-             ELSEIF (pplay(ig,l) .gt. p_bot) THEN
-               mode_dens = mode_dens*(pplay(ig,l+1)/pplay(ig,l))**(h_lay/h_bot)
-             ENDIF
-             mode_dens=max(1.e-30,mode_dens)
-! normalized to 0.35 microns (peak of absorption)
-             aerosol(ig,l,iaer) = QREFvis3d(ig,l,iaer)*mode_dens
-           ENDDO
-           mode_dens = nmodeuv ! m-3
-           DO l=ll+1,nlayer-1
-             h_lay=RD*(pt(ig,l-1)+pt(ig,l))/2./RG
-             IF (pplay(ig,l) .gt. p_top) THEN
-               mode_dens = nmodeuv ! m-3
-             ELSEIF (pplay(ig,l) .le. p_top) THEN
-               mode_dens = mode_dens*(pplay(ig,l)/pplay(ig,l-1))**(h_lay/h_top)
-             ENDIF
-             mode_dens=max(1.e-30,mode_dens)
-! normalized to 0.35 microns (peak of absorption)
-             aerosol(ig,l,iaer) = QREFvis3d(ig,l,iaer)*mode_dens
-           ENDDO
-          ENDDO
-
-!       3. Re-normalize to Haus et al 2015 total column optical depth
-         tau_col(:)=0.0
-         DO l=1,nlayer
-          DO ig=1,ngrid
+
+           if (iaer.eq.5) then  ! for UV absorber
+            DO l=1,nlayer
                tau_col(ig) = tau_col(ig) &
                      + aerosol(ig,l,iaer)
             ENDDO
-         ENDDO
-         DO ig=1,ngrid
-           DO l=1,nlayer-1
+            DO l=1,nlayer
                 aerosol(ig,l,iaer)=aerosol(ig,l,iaer)*0.205/tau_col(ig)
-           ENDDO
-         ENDDO
-
-      end if ! UV absorber 
+            ENDDO
+           endif
+
+        enddo
+      enddo
 
 !==================================================================
@@ -429,7 +206,7 @@
 
       tau_col(:)=0.0
-      do iaer = 1, naerkind
+      do ig=1,ngrid
          do l=1,nlayer
-            do ig=1,ngrid
+            do iaer = 1, naerkind
                tau_col(ig) = tau_col(ig) + aerosol(ig,l,iaer)
             end do
@@ -459,5 +236,149 @@
          endif
       end do
+
       return
     end subroutine aeropacity
+
+! --------------------------------------------------------------------------
+! --------------------------------------------------------------------------
+subroutine cloud_haus16(ngrid,pbot2,pbot,hbot2,hbot,ptop,ptop2,htop,htop2,nmode)
+
+!==================================================================
+!         Venus clouds (4 modes)
+!   S. Lebonnois (jan 2016)
+!==================================================================
+! distributions from Haus et al, 2016
+! ------------------------------------
+!
+! mode             1      2      2p     3
+! r (microns)     0.30   1.05   1.40   3.65
+! sigma           1.56   1.29   1.23   1.28
+! reff (microns)  0.49   1.23   1.56   4.25
+! nueff           0.21   0.067  0.044  0.063
+! (nueff=exp(ln^2 sigma)-1)
+!
+! pbot <=> zb ; ptop <=> zb+zc ; hbot <=> Hlo ; htop <=> Hup
+! p<ptop: N(i+1)=N(i)*(p(i+1)/p(i))**(hlay(i)/htop)      
+!    hlay=R(T(i)+T(i+1))/2/g  (in m)
+!    R=8.31/mu (mu in kg/mol)
+! p>pbot: N(i-1)=N(i)*(p(i)/p(i-1))**(hlay(i)/hbot)      
+! N is in m-3
+!
+! values in each mode below from Table 1 and text (p.186)
+!
+! dTau = Qext*[pi*reff**2*exp(-3*ln(1+nueff))]*N*hlay*(-dp)/p
+!
+! Latitudinal dependence (values from tables 3 and 4):
+! mode factor MF12(lat) for modes 1 and 2
+! mode 2' unchanged
+! mode factor MF3(lat) for mode 3
+! + for mode 2 only : pbot(lat), ptop(lat), htop(lat), htop2(lat)
+!==================================================================
+
+      use radinc_h, only : naerkind
+      USE geometry_mod, ONLY: latitude_deg             
+      implicit none
+       
+      INTEGER,INTENT(IN) :: ngrid  ! number of atmospheric columns
+      real,INTENT(out) :: pbot2(ngrid,naerkind),pbot(ngrid,naerkind)
+      real,INTENT(out) :: ptop(ngrid,naerkind),ptop2(ngrid,naerkind)
+      real,INTENT(out) :: hbot2(ngrid,naerkind),hbot(ngrid,naerkind)
+      real,INTENT(out) :: htop(ngrid,naerkind),htop2(ngrid,naerkind)
+      real,INTENT(out) :: nmode(ngrid,naerkind)
       
+! For latitudinal dependence
+      integer :: i,j,ilat
+      real :: latit,factlat
+      integer,parameter :: nlat_H16 = 16
+      real :: lat_H16(nlat_H16)
+      real :: MF12(nlat_H16),MF3(nlat_H16)
+      real :: pbotM2(nlat_H16),ptopM2(nlat_H16),htopM2(nlat_H16)
+
+!-------------------------
+!!   Equatorial model  
+!-------------------------
+! initialization
+      pbot2(:,:) = 1.e8
+      pbot(:,:)  = 1.e8
+      ptop2(:,:) = 0.
+      ptop(:,:)  = 0.
+      hbot2(:,:) = 1.
+      hbot(:,:)  = 1.
+      htop2(:,:) = 1.
+      htop(:,:)  = 1.
+      nmode(:,:) = 0.
+        
+! table of values(lat,mode)
+
+! Mode 1
+      pbot2(:,1) = 2.0e5 ! Pa
+      hbot2(:,1) = 5.0e3 ! m
+      pbot(:,1)  = 1.2e5
+      hbot(:,1)  = 1.0e3
+      nmode(:,1) = 1.935e8 ! m-3
+      ptop(:,1)  = 1.0e4
+      htop(:,1)  = 3.5e3
+      ptop2(:,1) = 4.5e2
+      htop2(:,1) = 2.0e3
+
+! Mode 2
+      pbot(:,2)  = 1.0e4
+      hbot(:,2)  = 3.0e3
+      nmode(:,2) = 1.00e8 ! m-3
+      ptop(:,2)  = 8.0e3
+      htop(:,2)  = 3.5e3
+      ptop2(:,2) = 4.5e2
+      htop2(:,2) = 2.0e3
+
+! Mode 2p
+      pbot(:,3)  = 1.2e5
+      hbot(:,3)  = 0.1e3
+      nmode(:,3) = 5.0e7 ! m-3
+      ptop(:,3)  = 2.3e4
+      htop(:,3)  = 1.0e3
+
+! Mode 3
+      pbot(:,4)  = 1.2e5
+      hbot(:,4)  = 0.5e3
+      nmode(:,4) = 1.4e7 ! m-3
+      ptop(:,4)  = 3.9e4
+      htop(:,4)  = 1.0e3
+
+! UV absorber
+      pbot(:,5)  = 3.3e4  ! 58 km
+      hbot(:,5)  = 1.0e3
+      nmode(:,5) = 1.00e7 !m-3
+      ptop(:,5)  = 3.7e3  ! 70 km
+      htop(:,5)  = 1.0e3
+
+!----------------------------
+!!   Latitudinal variations  
+!----------------------------
+      lat_H16(:) = (/0.,15.,20.,25.,30.,35.,40.,45.,50.,55.,60.,65.,70.,75.,80.,90./)
+      MF12(:) = (/0.98,0.98,0.99,1.00,0.98,0.94,0.86,0.81,0.73,0.67,0.64,0.61,0.59,0.47,0.36,0.36/)
+      MF3(:)  = (/1.30,1.30,1.26,1.23,1.17,1.13,1.06,1.03,1.04,1.09,1.22,1.51,1.82,2.02,2.09,2.09/)
+      pbotM2(:) = (/1.0e4,1.0e4,1.0e4,1.0e4,1.0e4,1.0e4,1.0e4,1.0e4,9.4e3, &
+                    9.2e3,9.1e3,9.6e3,1.14e4,1.21e4,1.25e4,1.25e4/)
+      ptopM2(:) = (/8.0e3,8.0e3,8.0e3,8.0e3,8.0e3,8.0e3,8.0e3,8.0e3,7.7e3, &
+                    7.5e3,7.5e3,8.0e3,9.2e3,1.0e4,1.06e4,1.06e4/)
+      htopM2(:) = (/3.5e3,3.5e3,3.5e3,3.5e3,3.5e3,3.5e3,3.5e3,3.5e3,3.4e3, &
+                    3.2e3,2.6e3,2.0e3,2.0e3,0.6e3,0.5e3,0.5e3/)
+      
+      do i=1,ngrid
+        latit = abs(latitude_deg(i))
+        ilat=2
+        do j=2,nlat_H16
+           if (latit.gt.lat_H16(j-1)) ilat=j
+        enddo
+        factlat    = (lat_H16(ilat)-latit)/(lat_H16(ilat)-lat_H16(ilat-1))
+        pbot(i,2)  = pbotM2(ilat)*(1.-factlat) + pbotM2(ilat-1)*factlat
+        ptop(i,2)  = ptopM2(ilat)*(1.-factlat) + ptopM2(ilat-1)*factlat
+        htop(i,2)  = htopM2(ilat)*(1.-factlat) + htopM2(ilat-1)*factlat
+        htop2(i,2) = min(htop2(i,2),htop(i,2))
+        nmode(i,1) = nmode(i,1)*(MF12(ilat)*(1.-factlat)+MF12(ilat-1)*factlat)
+        nmode(i,2) = nmode(i,2)*(MF12(ilat)*(1.-factlat)+MF12(ilat-1)*factlat)
+        nmode(i,3) = nmode(i,3)*(MF3(ilat) *(1.-factlat)+MF3(ilat-1) *factlat)
+      enddo
+      
+      return
+end subroutine cloud_haus16
Index: trunk/LMDZ.VENUS/libf/phyvenus/diffusion.h
===================================================================
--- trunk/LMDZ.VENUS/libf/phyvenus/diffusion.h	(revision 2580)
+++ trunk/LMDZ.VENUS/libf/phyvenus/diffusion.h	(revision 2622)
@@ -11,4 +11,4 @@
       parameter (Pdiff=1.e-1)      ! pressure below which diffusion is computed
       parameter (tdiffmin=5.d0)
-      parameter (dzres=0.1d0)	 ! grid resolution (km) for diffusion
+      parameter (dzres=0.2d0)	 ! grid resolution (km) for diffusion
       
Index: trunk/LMDZ.VENUS/libf/phyvenus/hrtherm.F
===================================================================
--- trunk/LMDZ.VENUS/libf/phyvenus/hrtherm.F	(revision 2580)
+++ trunk/LMDZ.VENUS/libf/phyvenus/hrtherm.F	(revision 2622)
@@ -76,5 +76,5 @@
          xabsi(3,i)  = rm(i,i_o)
          xabsi(8,i)  = rm(i,i_n2)
-         xabsi(11,i)  = rm(i,i_co)
+         xabsi(11,i) = rm(i,i_co)
 
 c         xabsi(6,i)  = rm(i,i_h2o2)
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 2622)
@@ -96,10 +96,12 @@
 ! generic : kmin=1/sqrt(DxDy) Dx and Dy horizontal grid
 
-!     REAL, parameter:: cmax = 61.   ! Max horizontal absolute phase velocity  TestGW6
+!----------------------------------------
+! VCD 1.1 tuning
+!     REAL, parameter:: cmax = 61.   ! Max horizontal absolute phase velocity
 !----------------------------------------
 ! VCD 2.0 tuning
-     REAL, parameter:: cmax = 301.   ! Max horizontal absolute phase velocity  TestGW6
-!----------------------------------------
-
+      REAL, parameter:: cmax = 61.   ! Max horizontal absolute phase velocity
+!----------------------------------------
+    
 ! generic : cmax=zonal wind at launch
      REAL, parameter:: cmin = 1.    ! Min horizontal absolute phase velocity
@@ -126,12 +128,27 @@
 
      ! 0.3.2 Parameters of waves dissipations
-     REAL, parameter:: sat   = 0.85   ! saturation parameter
-     REAL, parameter:: rdiss = 0.1    ! coefficient of dissipation
+!----------------------------------------
+! VCD 1.1 tuning
+!    REAL, parameter:: sat   = 0.85   ! saturation parameter
+!    REAL, parameter:: rdiss = 0.1    ! coefficient of dissipation
+!----------------------------------------
+! VCD 2.0 tuning
+     REAL, parameter:: sat   = 0.6    ! saturation parameter
+     REAL, parameter:: rdiss = 8.e-4  ! 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 min_k(ngrid)               ! min(kstar,kmin)
      REAL, save::  H0                ! characteristic height of the atmosphere
-     REAL, parameter:: pr = 5e5      ! Reference pressure [Pa]   ! VENUS: cloud layer
+     REAL, save::  MDR               ! characteristic mass density  
+!----------------------------------------
+! VCD 1.1 tuning
+!    REAL, parameter:: pr = 5e5      ! Reference pressure [Pa]   ! VENUS: cloud layer
+!----------------------------------------
+! VCD 2.0 tuning
+     REAL, parameter:: pr = 5e4      ! 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) 
@@ -156,4 +173,6 @@
         ! Characteristic vertical scale height
         H0 = RD * tr / RG
+        ! Characteristic mass density at launch altitude
+        MDR = pr / ( RD * tr )
         ! Control
         if (deltat .LT. dtime) THEN
@@ -214,4 +233,14 @@
      uh(:, nlayer + 1) = uu(:, nlayer)
      vh(:, nlayer + 1) = vv(:, nlayer)
+
+     ! TN+GG April/June 2020
+     ! "Individual waves are not supposed
+     ! to occupy the entire domain, but
+     ! only a faction of it" Lott+2012
+     ! minimum value of k between kmin and kstar
+     DO ii = 1, ngrid
+        kstar = RPI / sqrt(cell_area(ii)) 
+        min_k(ii) = max(kmin,kstar)
+     end DO
 
 ! 3. WAVES CHARACTERISTICS CHOSEN RANDOMLY
@@ -235,10 +264,5 @@
            ! 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
+           zk(jw, ii) = min_k(ii) + (kmax - min_k(ii)) * ran_num_2
 
            ! horizontal phase speed
@@ -320,5 +344,5 @@
            ! Intrinsic frequency
            intr_freq_p(jw, :) = zo(jw, :) - zk(jw, :) * cos(zp(jw, :)) * uh(:, ll + 1)    &
-                      - zk(jw, :) * sin(zp(jw, :)) * vh(:, ll + 1)
+                                          - zk(jw, :) * sin(zp(jw, :)) * vh(:, ll + 1)
            ! Vertical velocity in flux formulation
            wwp(jw, :) = min(                                                              &
@@ -336,5 +360,7 @@
                       * abs(intr_freq_p(jw, :))**3 / bv(:, ll + 1)                        &
                       * exp(-zh(:, ll + 1) / H0)                                          &
-                      * sat**2 * kmin**2 / zk(jw, :)**4)
+!! Correction for VCD 2.0
+!                      * sat**2 * kmin**2 / zk(jw, :)**4)
+                      * sat**2 * min_k(:)**2 * MDR / zk(jw, :)**4)                       
         end DO
 
Index: trunk/LMDZ.VENUS/libf/phyvenus/photochemistry_venus.F90
===================================================================
--- trunk/LMDZ.VENUS/libf/phyvenus/photochemistry_venus.F90	(revision 2580)
+++ trunk/LMDZ.VENUS/libf/phyvenus/photochemistry_venus.F90	(revision 2622)
@@ -1569,10 +1569,10 @@
 !-- 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.
+!   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(65:90,4) = v_phot(65:90,4)*10.
 !--
 !   v_phot(:,4) = v_phot(:,4)*10.
Index: trunk/LMDZ.VENUS/libf/phyvenus/physiq_mod.F
===================================================================
--- trunk/LMDZ.VENUS/libf/phyvenus/physiq_mod.F	(revision 2580)
+++ trunk/LMDZ.VENUS/libf/phyvenus/physiq_mod.F	(revision 2622)
@@ -1449,7 +1449,8 @@
             mmean(:,:) = 0.
             do iq = 1,nqmax - nmicro
-               mmean(:,:) = mmean(:,:)+tr_seri(:,:,iq)*m_tr(iq)
-               rnew(:,:) = 8.314/mmean(:,:)*1.e3     ! J/kg K 
+               mmean(:,:) = mmean(:,:)+tr_seri(:,:,iq)/m_tr(iq)
             enddo
+            mmean(:,:) = 1./mmean(:,:)
+            rnew(:,:) = 8.314/mmean(:,:)*1.e3     ! J/kg K 
          endif
 
@@ -1812,15 +1813,15 @@
 ! 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,:),
+!     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, for VCD 2.0
+      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)
-
-! 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
Index: trunk/LMDZ.VENUS/libf/phyvenus/phytrac_chimie.F
===================================================================
--- trunk/LMDZ.VENUS/libf/phyvenus/phytrac_chimie.F	(revision 2580)
+++ trunk/LMDZ.VENUS/libf/phyvenus/phytrac_chimie.F	(revision 2622)
@@ -74,9 +74,9 @@
       if (debutphy) then
       
-!!--- Adjustment of Helium amount
-!        if (i_he/=0) then
-!           trac(:,:,i_he)=trac(:,:,i_he)*20.
-!        endif
-!!---
+!--- Adjustment of Helium amount
+!       if (i_he/=0) then
+!          trac(:,:,i_he)=trac(:,:,i_he)/2.
+!       endif
+!---
 
 !-------------------------------------------------------------------
@@ -86,17 +86,28 @@
 
 !!! in this reinitialisation, trac is VOLUME mixing ratio
-            print*, "Tracers are re-initialised"
-            trac(:,:,:) = 1.0e-30
-
-            if ((i_ocs /= 0) .and. (i_co /= 0) .and. (i_hcl /= 0) 
-     $           .and. (i_so2 /= 0) .and. (i_h2o /= 0) .and. (i_n2/ = 0)
-     $           .and. (i_co2 /= 0)) then
-
-               trac(:,1:22,i_ocs) = 3.e-6
-               trac(:,1:22,i_co)  = 25.e-6
-               trac(:,:,i_hcl)    = 0.4e-6
-               trac(:,1:22,i_so2) = 7.e-6
-               trac(:,1:22,i_h2o) = 30.e-6
-               trac(:,:,i_n2)     = 0.35e-1
+! ONLY SO2!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!           convert mass to volume mixing ratio
+            do iq = 1,nqmax - nmicro
+               trac(:,:,iq) = trac(:,:,iq)*mmean(:,:)/m_tr(iq)
+            end do
+            print*, "SO2 is re-initialised"
+            if (i_so2 /= 0) then 
+               trac(:,1:22,i_so2) = 10.e-6
+
+! AL TRACERS!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!           print*, "Tracers are re-initialised"
+!           trac(:,:,:) = 1.0e-30
+
+!           if ((i_ocs /= 0) .and. (i_co /= 0) .and. (i_hcl /= 0) 
+!    $           .and. (i_so2 /= 0) .and. (i_h2o /= 0) .and. (i_n2/ = 0)
+!    $           .and. (i_co2 /= 0)) then
+
+!              trac(:,1:22,i_ocs) = 3.e-6
+!              trac(:,1:22,i_co)  = 25.e-6
+!              trac(:,:,i_hcl)    = 0.4e-6
+!              trac(:,1:22,i_so2) = 7.e-6
+!              trac(:,1:22,i_h2o) = 30.e-6
+!              trac(:,:,i_n2)     = 0.35e-1
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    
 !          ensure that sum of mixing ratios = 1
@@ -124,6 +135,6 @@
             do iq = 1,nqmax - nmicro
                mmean(:,:) = mmean(:,:)+trac(:,:,iq)*m_tr(iq)
-               rnew(:,:) = 8.314/mmean(:,:)*1.e3     ! J/kg K 
             enddo
+            rnew(:,:) = 8.314/mmean(:,:)*1.e3     ! J/kg K 
 
 !           convert volume to mass mixing ratio
@@ -293,8 +304,11 @@
             mmean(:,:) = 0.
             do iq = 1,nqmax - nmicro
-               mmean(:,:) = mmean(:,:)+trac(:,:,iq)*m_tr(iq)
-               rnew(:,:) = 8.314/mmean(:,:)*1.e3     ! J/kg K 
+               mmean(:,:) = mmean(:,:)+ztrac(:,:,iq)*m_tr(iq)
             enddo
-
+            rnew(:,:) = 8.314/mmean(:,:)*1.e3     ! J/kg K 
+
+!===================================================================
+!     convert volume to mass mixing ratio / then tendencies in mmr
+!===================================================================
 !     gas phase
 
