Index: trunk/LMDZ.GENERIC/libf/phystd/aeropacity.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/aeropacity.F90	(revision 2810)
+++ trunk/LMDZ.GENERIC/libf/phystd/aeropacity.F90	(revision 2831)
@@ -1,4 +1,5 @@
-      Subroutine aeropacity(ngrid,nlayer,nq,pplay,pplev,pq, &
-         aerosol,reffrad,QREFvis3d,QREFir3d,tau_col,cloudfrac,totcloudfrac,clearsky)
+      Subroutine aeropacity(ngrid,nlayer,nq,pplay,pplev,pt, pq, &
+         aerosol,reffrad,nueffrad, QREFvis3d,QREFir3d,tau_col, &
+         cloudfrac,totcloudfrac,clearsky)
 
        use radinc_h, only : L_TAUMAX,naerkind
@@ -6,5 +7,7 @@
                               iaero_aurora, iaero_back2lay, iaero_co2, &
                               iaero_dust, iaero_h2o, iaero_h2so4, &
-                              iaero_nh3, i_rgcs_ice, noaero
+                              iaero_nh3, i_rgcs_ice, noaero, &
+                              iaero_venus1, iaero_venus2, iaero_venus2p, &
+                              iaero_venus3, iaero_venusUV
        USE tracer_h, only: noms,rho_co2,rho_ice,rho_q,mmol
        use comcstfi_mod, only: g, pi, mugaz, avocado
@@ -60,6 +63,8 @@
       REAL,INTENT(IN) :: pplev(ngrid,nlayer+1) ! inter-layer pressure (Pa)
       REAL,INTENT(IN) :: pq(ngrid,nlayer,nq) ! tracers (.../kg_of_air)
+      REAL,INTENT(IN) :: pt(ngrid,nlayer) ! mid-layer temperature (K)
       REAL,INTENT(OUT) :: aerosol(ngrid,nlayer,naerkind) ! aerosol optical depth
       REAL,INTENT(IN) :: reffrad(ngrid,nlayer,naerkind) ! aerosol effective radius
+      REAL,INTENT(IN) :: nueffrad(ngrid,nlayer,naerkind) ! aerosol effective variance
       REAL,INTENT(IN) :: QREFvis3d(ngrid,nlayer,naerkind) ! extinction coefficient in the visible
       REAL,INTENT(IN) :: QREFir3d(ngrid,nlayer,naerkind)
@@ -98,4 +103,7 @@
       logical dummy_bool ! dummy boolean just in case we need one
       ! integer i_rgcs_ice(aerogeneric)
+      !  for venus clouds
+      real      :: p_bot,p_top,h_bot,h_top,mode_dens,h_lay
+
       ! identify tracers
       IF (firstcall) THEN
@@ -148,4 +156,21 @@
           print*,'iaero_aurora= ',iaero_aurora
         endif
+
+        if (iaero_venus1.ne.0) then
+          print*,'iaero_venus1= ',iaero_venus1
+        endif
+        if (iaero_venus2.ne.0) then
+          print*,'iaero_venus2= ',iaero_venus2
+        endif
+        if (iaero_venus2p.ne.0) then
+          print*,'iaero_venus2p= ',iaero_venus2p
+        endif
+        if (iaero_venus3.ne.0) then
+          print*,'iaero_venus3= ',iaero_venus3
+        endif
+        if (iaero_venusUV.ne.0) then
+          print*,'iaero_venusUV= ',iaero_venusUV
+        endif
+
         if (aerogeneric .ne. 0) then 
           print*,"iaero_generic= ",iaero_generic(:)
@@ -679,4 +704,276 @@
       endif !aerogeneric .ne. 0
 
+!==================================================================
+!         Venus clouds (4 modes)
+!   S. Lebonnois (jan 2016)
+!==================================================================
+! distributions from Haus et al, 2013
+! 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=No*(p/p_top)**(h_lay/h_top)      h_lay=RT/g  (in m)
+! p>p_bot: N=No*(p_bot/p)**(h_lay/h_bot)      R=8.31/mu (mu in kg/mol)
+! 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
+          p_bot = 1.e5
+          p_top = 1.e4
+          h_bot = 1.0e3 ! m
+          h_top = 5.0e3
+          
+!       2. Opacity calculation
+
+          DO ig=1,ngrid
+           DO l=1,nlayer-1
+
+             h_lay=8.31*pt(ig,l)/(g*0.044)
+
+             !! 1. below 2e5 Pa: no aerosols
+             IF (pplay(ig,l) .gt. 2.e5) THEN
+               mode_dens = 0.
+
+             !! 2. cloud layer
+             ELSEIF (pplay(ig,l) .le. 2.e5 .and. pplay(ig,l) .gt. p_bot) THEN
+               mode_dens = 1.81e8*(p_bot/pplay(ig,l))**(h_lay/h_bot)
+               
+             ELSEIF (pplay(ig,l) .le. p_bot .and. pplay(ig,l) .gt. p_top) THEN
+               mode_dens = 1.81e8  ! m-3
+               
+             ELSEIF (pplay(ig,l) .le. p_top .and. pplay(ig,l) .gt. 1.e2) THEN
+               mode_dens = 1.81e8*(pplay(ig,l)/p_top)**(h_lay/h_top)
+               
+             !! 3. above 1.e2 Pa: no aerosols
+             ELSEIF (pplay(ig,l) .le. 1.e2) THEN
+               mode_dens = 0.
+             ENDIF
+
+             aerosol(ig,l,iaer) = QREFvis3d(ig,l,iaer)*                       &
+              pi*(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)
+
+           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.1e4
+          p_top = 1.e4
+          h_bot = 3.0e3
+          h_top = 3.5e3
+          
+!       2. Opacity calculation
+
+          DO ig=1,ngrid
+           DO l=1,nlayer-1
+
+             h_lay=8.31*pt(ig,l)/(g*0.044)
+
+             !! 1. below 2e5 Pa: no aerosols
+             IF (pplay(ig,l) .gt. 2.e5) THEN
+               mode_dens = 0.
+
+             !! 2. cloud layer
+             ELSEIF (pplay(ig,l) .le. 2.e5 .and. pplay(ig,l) .gt. p_bot) THEN
+               mode_dens = 1.00e8*(p_bot/pplay(ig,l))**(h_lay/h_bot)
+               
+             ELSEIF (pplay(ig,l) .le. p_bot .and. pplay(ig,l) .gt. p_top) THEN
+               mode_dens = 1.00e8
+               
+             ELSEIF (pplay(ig,l) .le. p_top .and. pplay(ig,l) .gt. 1.e2) THEN
+               mode_dens = 1.00e8*(pplay(ig,l)/p_top)**(h_lay/h_top)
+               
+             !! 3. above 1.e2 Pa: no aerosols
+             ELSEIF (pplay(ig,l) .le. 1.e2) THEN
+               mode_dens = 0.
+             ENDIF
+
+             aerosol(ig,l,iaer) = QREFvis3d(ig,l,iaer)*                       &
+              pi*(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)
+
+           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.e5
+          p_top = 2.3e4
+          h_bot = 0.1e3
+          h_top = 1.0e3
+          
+!       2. Opacity calculation
+
+          DO ig=1,ngrid
+           DO l=1,nlayer-1
+
+             h_lay=8.31*pt(ig,l)/(g*0.044)
+
+             !! 1. below 2e5 Pa: no aerosols
+             IF (pplay(ig,l) .gt. 2.e5) THEN
+               mode_dens = 0.
+
+             !! 2. cloud layer
+             ELSEIF (pplay(ig,l) .le. 2.e5 .and. pplay(ig,l) .gt. p_bot) THEN
+               mode_dens = 5.00e7*(p_bot/pplay(ig,l))**(h_lay/h_bot)
+               
+             ELSEIF (pplay(ig,l) .le. p_bot .and. pplay(ig,l) .gt. p_top) THEN
+               mode_dens = 5.00e7
+               
+             ELSEIF (pplay(ig,l) .le. p_top .and. pplay(ig,l) .gt. 1.e2) THEN
+               mode_dens = 5.00e7*(pplay(ig,l)/p_top)**(h_lay/h_top)
+               
+             !! 3. above 1.e2 Pa: no aerosols
+             ELSEIF (pplay(ig,l) .le. 1.e2) THEN
+               mode_dens = 0.
+             ENDIF
+
+             aerosol(ig,l,iaer) = QREFvis3d(ig,l,iaer)*                       &
+              pi*(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)
+
+           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.e5
+          p_top = 4.e4
+          h_bot = 0.5e3
+          h_top = 1.0e3
+          
+!       2. Opacity calculation
+
+          DO ig=1,ngrid
+           DO l=1,nlayer-1
+ 
+              h_lay=8.31*pt(ig,l)/(g*0.044)
+
+             !! 1. below 2e5 Pa: no aerosols
+             IF (pplay(ig,l) .gt. 2.e5) THEN
+               mode_dens = 0.
+
+             !! 2. cloud layer
+             ELSEIF (pplay(ig,l) .le. 2.e5 .and. pplay(ig,l) .gt. p_bot) THEN
+               mode_dens = 1.40e7*(p_bot/pplay(ig,l))**(h_lay/h_bot)
+               
+             ELSEIF (pplay(ig,l) .le. p_bot .and. pplay(ig,l) .gt. p_top) THEN
+               mode_dens = 1.40e7
+               
+             ELSEIF (pplay(ig,l) .le. p_top .and. pplay(ig,l) .gt. 1.e2) THEN
+               mode_dens = 1.40e7*(pplay(ig,l)/p_top)**(h_lay/h_top)
+               
+             !! 3. above 1.e2 Pa: no aerosols
+             ELSEIF (pplay(ig,l) .le. 1.e2) THEN
+               mode_dens = 0.
+             ENDIF
+
+             aerosol(ig,l,iaer) = QREFvis3d(ig,l,iaer)*                       &
+              pi*(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)
+
+           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
+          p_top = 3.7e3 ! 70 km
+          h_bot = 1.0e3 
+          h_top = 1.0e3
+          
+!       2. Opacity calculation
+
+          DO ig=1,ngrid
+           DO l=1,nlayer-1
+
+             h_lay=8.31*pt(ig,l)/(g*0.044)
+
+             !! 1. below 7e4 Pa: no aerosols
+             IF (pplay(ig,l) .gt. 7.e4) THEN
+               mode_dens = 0.
+
+             !! 2. cloud layer
+             ELSEIF (pplay(ig,l) .le. 7.e4 .and. pplay(ig,l) .gt. p_bot) THEN
+               mode_dens = 1.00e7*(p_bot/pplay(ig,l))**(h_lay/h_bot)
+               
+             ELSEIF (pplay(ig,l) .le. p_bot .and. pplay(ig,l) .gt. p_top) THEN
+               mode_dens = 1.00e7
+               
+             ELSEIF (pplay(ig,l) .le. p_top .and. pplay(ig,l) .gt. 1.e3) THEN
+               mode_dens = 1.00e7*(pplay(ig,l)/p_top)**(h_lay/h_top)
+               
+             !! 3. above 1.e3 Pa: no aerosols
+             ELSEIF (pplay(ig,l) .le. 1.e3) THEN
+               mode_dens = 0.
+             ENDIF
+
+! 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
+               tau_col(ig) = tau_col(ig) &
+                     + aerosol(ig,l,iaer)
+            ENDDO
+         ENDDO
+         DO ig=1,ngrid
+           DO l=1,nlayer-1
+                aerosol(ig,l,iaer)=aerosol(ig,l,iaer)*0.205/tau_col(ig)
+           ENDDO
+         ENDDO
+
+      end if ! UV absorber 
+
+!==================================================================
+!     ig=10
+!      do l=1,nlayer
+!          print*,8.31*pt(ig,l)/(g*0.044),pplay(ig,l),aerosol(ig,l,1),aerosol(ig,l,2),aerosol(ig,l,3),aerosol(ig,l,4)
+!         print*,l,pplay(ig,l),aerosol(ig,l,5)
+!      enddo
+!      stop            
+!==================================================================
+
+
 ! --------------------------------------------------------------------------
 ! Column integrated visible optical depth in each point (used for diagnostic)
@@ -713,5 +1010,5 @@
       !    endif
       ! end do
-      return
+
     end subroutine aeropacity
       
