Index: trunk/LMDZ.GENERIC/README
===================================================================
--- trunk/LMDZ.GENERIC/README	(revision 2822)
+++ trunk/LMDZ.GENERIC/README	(revision 2831)
@@ -1754,2 +1754,12 @@
 was tested although not defined; using tests on "aerogeneric" (number
 of generic aerosols) works better.
+
+== 17/11/2022 == LT
+Removed extrapolations in tpindex. If above Tmax or Pmax use the max value.
+
+== 23/11/2022 == GG
+Add the possibility to include Venus-like aerosols (triggered by option
+aerovenus=.true. in callphys.def); baseline is to use 5 distinct scatterers
+but each may be turned on/off (via aerovenus1, aerovenus2, aerovenus2p,
+aerovenus3, aerovenusUV flags which may be specified in callphys.def).
+
Index: trunk/LMDZ.GENERIC/libf/phystd/aeropacity.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/aeropacity.F90	(revision 2822)
+++ 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
       
Index: trunk/LMDZ.GENERIC/libf/phystd/aerosol_mod.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/aerosol_mod.F90	(revision 2822)
+++ trunk/LMDZ.GENERIC/libf/phystd/aerosol_mod.F90	(revision 2831)
@@ -13,4 +13,5 @@
       integer :: iaero_h2so4 = 0
       logical :: noaero = .false.
+!$OMP THREADPRIVATE(iaero_co2,iaero_h2o,iaero_dust,iaero_h2so4,noaero)
 
 ! two-layer simple aerosol model
@@ -22,8 +23,20 @@
 ! Auroral aerosols
       integer :: iaero_aurora = 0
+!$OMP THREADPRIVATE(iaero_back2lay,iaero_nh3,iaero_nlay,iaero_aurora)
+
 ! Generic aerosols
       integer, dimension(:),allocatable,save :: iaero_generic
       integer,dimension(:),allocatable,save :: i_rgcs_ice
-!$OMP THREADPRIVATE(iaero_co2,iaero_h2o,iaero_dust,iaero_h2so4,noaero,iaero_back2lay,iaero_nh3,iaero_nlay,iaero_aurora,iaero_generic,i_rgcs_ice)
+!$OMP THREADPRIVATE(iaero_generic,i_rgcs_ice)
+
+! Venus clouds
+      integer :: iaero_venus1 = 0
+      integer :: iaero_venus2 = 0
+      integer :: iaero_venus2p = 0
+      integer :: iaero_venus3 = 0
+      integer :: iaero_venusUV = 0
+!$OMP THREADPRIVATE(iaero_venus1,iaero_venus2,iaero_venus2p)
+!$OMP THREADPRIVATE(iaero_venus3,iaero_venusUV)
+
       
 !==================================================================
Index: trunk/LMDZ.GENERIC/libf/phystd/calc_cpp_mugaz.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/calc_cpp_mugaz.F90	(revision 2822)
+++ trunk/LMDZ.GENERIC/libf/phystd/calc_cpp_mugaz.F90	(revision 2831)
@@ -24,4 +24,5 @@
       implicit none
 
+      character(len=20),parameter :: myname="calc_cpp_mugaz"
       real cpp_c   
       real mugaz_c
@@ -64,7 +65,19 @@
                ! C2H2 http://encyclopedia.airliquide.com/Encyclopedia.asp?GasID=1
                mugaz_c = mugaz_c + 26.04*gfrac(igas)
+            ! GG MODIF JAN2019
+            elseif(igas.eq.igas_CO)then
+               mugaz_c = mugaz_c + 28.01*gfrac(igas)
+            elseif(igas.eq.igas_OCS)then
+               ! https://pubchem.ncbi.nlm.nih.gov/compound/carbonyl_sulfide
+               mugaz_c = mugaz_c + 60.0751*gfrac(igas)
+            elseif(igas.eq.igas_HCl)then
+               mugaz_c = mugaz_c + 36.46*gfrac(igas)
+            elseif(igas.eq.igas_HF)then
+               mugaz_c = mugaz_c + 20.01*gfrac(igas)
             else
+              if (is_master) then
                print*,'Error in calc_cpp_mugaz: Gas species not recognised!'
-               call abort
+              endif
+              call abort_physic(myname,'Gas species not recognised!',1)
             endif
          endif
@@ -108,7 +121,18 @@
                ! http://encyclopedia.airliquide.com/Encyclopedia.asp?GasID=1
                cpp_c   = cpp_c   + 1.575*gfrac(igas)*26.04/mugaz_c
+            !!!!! MODIF GG JAN 2019  (check source values !!)
+            elseif(igas.eq.igas_CO)then
+               cpp_c   = cpp_c   + 1.0425*gfrac(igas)*28.01/mugaz_c
+            elseif(igas.eq.igas_OCS)then
+               cpp_c   = cpp_c   + 0.6909*gfrac(igas)*60.07/mugaz_c
+            elseif(igas.eq.igas_HCl)then
+               cpp_c   = cpp_c   + 1.7087*gfrac(igas)*36.46/mugaz_c
+            elseif(igas.eq.igas_HF)then
+               cpp_c   = cpp_c   + 1.4553*gfrac(igas)*20.01/mugaz_c
             else
+              if (is_master) then
                print*,'Error in calc_cpp_mugaz: Gas species not recognised!'
-               call abort
+              endif
+              call abort_physic(myname,'Gas species not recognised!',1)
             endif
          endif
@@ -122,6 +146,8 @@
       cpp_c = 1000.0*cpp_c
 
-      print*,'Cp in calc_cpp_mugaz is ',cpp_c,'J kg^-1 K^-1'
-      print*,'Mg in calc_cpp_mugaz is ',mugaz_c,'amu'
+      if (is_master) then
+        print*,'Cp in calc_cpp_mugaz is ',cpp_c,'J kg^-1 K^-1'
+        print*,'Mg in calc_cpp_mugaz is ',mugaz_c,'amu'
+      endif
 
       if(((abs(1.-cpp/cpp_c).gt.1.e-6) .or.  &
@@ -146,4 +172,3 @@
 
 
-      return
     end subroutine calc_cpp_mugaz
Index: trunk/LMDZ.GENERIC/libf/phystd/callcorrk.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/callcorrk.F90	(revision 2822)
+++ trunk/LMDZ.GENERIC/libf/phystd/callcorrk.F90	(revision 2831)
@@ -28,10 +28,14 @@
       use gases_h, only: ngasmx
       use radii_mod, only : su_aer_radii,co2_reffrad,h2o_reffrad,dust_reffrad,h2so4_reffrad,back2lay_reffrad
-      use aerosol_mod, only : iaero_co2,iaero_h2o,iaero_dust,iaero_h2so4, iaero_back2lay, iaero_aurora
+      use aerosol_mod, only : iaero_co2,iaero_h2o,iaero_dust,iaero_h2so4, &
+                              iaero_back2lay, iaero_aurora,               &
+                              iaero_venus1,iaero_venus2,iaero_venus2p,    &
+                              iaero_venus3,iaero_venusUV
       use tracer_h
       use comcstfi_mod, only: pi, mugaz, cpp
-      use callkeys_mod, only: varactive,diurnal,tracer,water,varfixed,satval,diagdtau,    &
-                              kastprof,strictboundcorrk,specOLR,CLFvarying,               &
-                              tplanckmin,tplanckmax,global1d,generic_condensation
+      use callkeys_mod, only: varactive,diurnal,tracer,water,varfixed,satval, &
+                              diagdtau,kastprof,strictboundcorrk,specOLR, &
+                              CLFvarying,tplanckmin,tplanckmax,global1d, &
+                              generic_condensation, aerovenus
       use optcv_mod, only: optcv
       use optci_mod, only: optci
@@ -294,11 +298,7 @@
 
 #ifndef MESOSCALE
-         call system('rm -f surf_vals_long.out')
+         if (is_master) call system('rm -f surf_vals_long.out')
 #endif
 
-         ! if(naerkind.gt.4)then
-         !    message='Code not general enough to deal with naerkind > 4 yet.'
-         !    call abort_physic(subname,message,1)
-         ! endif
          call su_aer_radii(ngrid,nlayer,reffrad,nueffrad)
          
@@ -498,6 +498,6 @@
 
       ! Get aerosol optical depths.
-      call aeropacity(ngrid,nlayer,nq,pplay,pplev,pq,aerosol,      &
-           reffrad,QREFvis3d,QREFir3d,                             & 
+      call aeropacity(ngrid,nlayer,nq,pplay,pplev, pt,pq,aerosol,      &
+           reffrad,nueffrad,QREFvis3d,QREFir3d,                             & 
            tau_col,cloudfrac,totcloudfrac,clearsky)                
  
@@ -904,9 +904,17 @@
       plevrad(1) = 0.
 !      plevrad(2) = 0.   !! JL18 enabling this line puts the radiative top at p=0 which was the idea before, but does not seem to perform best after all. 
+      if (aerovenus) then
+!!  GG19 modified below after SL routines
+        plevrad(2) = 0.
+      endif
 
       tlevrad(1) = tlevrad(2)
       tlevrad(2*nlayer+1)=tsurf(ig)
       
-      pmid(1) = pplay(ig,nlayer)/scalep
+      pmid(1) = pplay(ig,nlayer)/scalep   
+      if (aerovenus) then
+!! GG19 modified below after SL routines
+        pmid(1) = max(pgasmin,0.0001*plevrad(3))
+      endif
       pmid(2) =  pmid(1)
 
Index: trunk/LMDZ.GENERIC/libf/phystd/callkeys_mod.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/callkeys_mod.F90	(revision 2822)
+++ trunk/LMDZ.GENERIC/libf/phystd/callkeys_mod.F90	(revision 2831)
@@ -48,4 +48,11 @@
       logical,save :: aeronh3, aeronlay, aeroaurora
 !$OMP THREADPRIVATE(aeronh3,aeronlay,aeroaurora)
+
+      logical,save :: aerovenus ! master flag for "Venus-like" aerosol additions
+!$OMP THREADPRIVATE(aerovenus)
+      ! detailed sub-options when with "Venus-like" aerosol additions
+      logical,save :: aerovenus1,aerovenus2,aerovenus2p,aerovenus3,aerovenusUV
+!$OMP THREADPRIVATE(aerovenus1,aerovenus2,aerovenus2p,aerovenus3,aerovenusUV)
+
       logical,save :: aerofixco2,aerofixh2o
 !$OMP THREADPRIVATE(aerofixco2,aerofixh2o)
Index: trunk/LMDZ.GENERIC/libf/phystd/gases_h.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/gases_h.F90	(revision 2822)
+++ trunk/LMDZ.GENERIC/libf/phystd/gases_h.F90	(revision 2831)
@@ -31,7 +31,12 @@
       integer :: igas_C2H2
       integer :: igas_C2H6
+!! GG MODIF Jan 2019
+      integer :: igas_OCS
+      integer :: igas_HCl
+      integer :: igas_HF            
 !!$OMP THREADPRIVATE(ngasmx,vgas,gnom,gfrac,&
 !	!$OMP igas_H2,igas_He,igas_H2O,igas_CO2,igas_CO,igas_N2,&
 !	!$OMP igas_O2,igas_SO2,igas_H2S,igas_CH4,igas_NH3,igas_C2H2,igas_C2H6)
-
+! The saved variables here are not OMP THREADPRIVATE because loaded by
+! the OMP master only (see su_gases.F90) and thus shared by all threads
       end module gases_h
Index: trunk/LMDZ.GENERIC/libf/phystd/iniaerosol.F
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/iniaerosol.F	(revision 2822)
+++ trunk/LMDZ.GENERIC/libf/phystd/iniaerosol.F	(revision 2831)
@@ -2,4 +2,5 @@
 
 
+      use mod_phys_lmdz_para, only : is_master
       use radinc_h, only: naerkind
       use aerosol_mod
@@ -7,5 +8,7 @@
       use callkeys_mod, only: aeroco2,aeroh2o,dusttau,aeroh2so4,
      &		aeroback2lay,aeronh3, nlayaero, aeronlay, aeroaurora,
-     &      aerogeneric
+     &          aerogeneric,
+     &          aerovenus1,aerovenus2,aerovenus2p,aerovenus3,aerovenusUV
+
 
       IMPLICIT NONE
@@ -56,5 +59,5 @@
          iaero_co2=ia
       endif
-      write(*,*) '--- CO2 aerosol = ', iaero_co2
+      if (is_master) write(*,*) '--- CO2 aerosol = ', iaero_co2
  
       if (aeroh2o) then
@@ -62,5 +65,5 @@
          iaero_h2o=ia
       endif
-      write(*,*) '--- H2O aerosol = ', iaero_h2o
+      if (is_master) write(*,*) '--- H2O aerosol = ', iaero_h2o
 
       if (dusttau.gt.0) then
@@ -68,5 +71,5 @@
          iaero_dust=ia
       endif
-      write(*,*) '--- Dust aerosol = ', iaero_dust
+      if (is_master) write(*,*) '--- Dust aerosol = ', iaero_dust
 
       if (aeroh2so4) then
@@ -74,5 +77,5 @@
          iaero_h2so4=ia
       endif
-      write(*,*) '--- H2SO4 aerosol = ', iaero_h2so4
+      if (is_master) write(*,*) '--- H2SO4 aerosol = ', iaero_h2so4
       
       if (aeroback2lay) then
@@ -80,5 +83,6 @@
          iaero_back2lay=ia
       endif
-      write(*,*) '--- Two-layer aerosol = ', iaero_back2lay
+      if (is_master) write(*,*) '--- Two-layer aerosol = ',
+     & iaero_back2lay
 
       if (aeronh3) then
@@ -86,5 +90,5 @@
          iaero_nh3=ia
       endif
-      write(*,*) '--- NH3 Cloud = ', iaero_nh3
+      if (is_master) write(*,*) '--- NH3 Cloud = ', iaero_nh3
 
       if (aeronlay) then
@@ -94,5 +98,5 @@
          enddo
       endif
-      write(*,*) '--- N-layer aerosol = ', iaero_nlay
+      if (is_master) write(*,*) '--- N-layer aerosol = ', iaero_nlay
 
       if (aeroaurora) then
@@ -100,5 +104,41 @@
          iaero_aurora=ia
       endif
-      write(*,*) '--- Auroral aerosols = ', iaero_aurora
+      if (is_master) write(*,*) '--- Auroral aerosols = ', iaero_aurora
+
+      if (aerovenus1) then
+         ia=ia+1
+         iaero_venus1=ia
+      endif
+      if (is_master) write(*,*) '--- Venus cloud, mode 1 aerosol = ',
+     & iaero_venus1
+
+      if (aerovenus2) then
+         ia=ia+1
+         iaero_venus2=ia
+      endif
+      if (is_master) write(*,*) '--- Venus cloud, mode 2 aerosol = ',
+     & iaero_venus2
+
+      if (aerovenus2p) then
+         ia=ia+1
+         iaero_venus2p=ia
+      endif
+      if (is_master) write(*,*) '--- Venus cloud, mode 2p aerosol = ',
+     & iaero_venus2p
+
+      if (aerovenus3) then
+         ia=ia+1
+         iaero_venus3=ia
+      endif
+      if (is_master) write(*,*) '--- Venus cloud, mode 3 aerosol = ',
+     & iaero_venus3
+
+      if (aerovenusUV) then
+         ia=ia+1
+         iaero_venusUV=ia
+      endif
+      if (is_master) write(*,*) '--- Venus cloud, UV absorber = ',
+     & iaero_venusUV
+
 
       if (aerogeneric .ne. 0) then 
@@ -108,7 +148,11 @@
          enddo
       endif 
-      write(*,*)'--- Radiative Generic Condensable Species = '
+      
+      if (is_master) then
+        write(*,*)'--- Radiative Generic Condensable Species = '
      &,iaero_generic
-      write(*,*) '=== Number of aerosols= ', ia
+
+        write(*,*) '=== Number of aerosols= ', ia
+      endif ! of is_master
 
 ! For the zero aerosol case, we currently make a dummy co2 aerosol which is zero everywhere.
@@ -123,9 +167,12 @@
 
       if (ia.ne.naerkind) then
+        if (is_master) then
           print*, 'Aerosols counted not equal to naerkind'
           print*, 'Compile with tag -s',ia,'to run'
           print*, 'or change options in callphys.def'
           print*, 'Abort in iniaerosol.F'
-          call abort
+        endif
+        call abort_physic("iniaerosl",
+     &                    'wrong number of aerosols',1)
       endif
 
Index: trunk/LMDZ.GENERIC/libf/phystd/inifis_mod.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/inifis_mod.F90	(revision 2822)
+++ trunk/LMDZ.GENERIC/libf/phystd/inifis_mod.F90	(revision 2831)
@@ -687,4 +687,42 @@
      endif
 
+     if (is_master) write(*,*)"Radiatively active Venus clouds ?"
+     aerovenus=.false. ! default value
+     call getin_p("aerovenus",aerovenus)
+     if (aerovenus) then
+       aerovenus1=.true.     ! default value
+       aerovenus2=.true.     ! default value
+       aerovenus2p=.true.     ! default value
+       aerovenus3=.true.     ! default value
+       aerovenusUV=.true.     ! default value
+     else
+       aerovenus1=.false.     ! default value
+       aerovenus2=.false.     ! default value
+       aerovenus2p=.false.     ! default value
+       aerovenus3=.false.     ! default value
+       aerovenusUV=.false.     ! default value
+     endif
+     ! in case the user wants to specifically set/unset sub-options
+     call getin_p("aerovenus1",aerovenus1)
+     if (is_master) write(*,*)" aerovenus1 = ",aerovenus1
+     call getin_p("aerovenus2",aerovenus2)
+     if (is_master) write(*,*)" aerovenus2 = ",aerovenus2
+     call getin_p("aerovenus2p",aerovenus2p)
+     if (is_master) write(*,*)" aerovenus2p= ",aerovenus2p
+     call getin_p("aerovenus3",aerovenus3)
+     if (is_master) write(*,*)" aerovenus3 = ",aerovenus3
+     call getin_p("aerovenusUV",aerovenusUV)
+     if (is_master) write(*,*)" aerovenusUV= ",aerovenusUV
+     ! Sanity check: if any of the aerovenus* is set to true
+     ! then aeronovenus should also be set to true
+     if ((.not.aerovenus).and.(aerovenus1.or.aerovenus2.or.aerovenus2p.or.&
+                               aerovenus3.or.aerovenusUV)) then
+      if(is_master) then
+       write(*,*)" Error, if you set some of the aerovenus* to true"
+       write(*,*)" then flag aerovenus should be set to true as well!"
+      endif
+      call abort_physic("inifis"," aerovenus* flags mismatch!",1)
+     endif
+     
      if (is_master) write(*,*)trim(rname)//&
        ": TWOLAY AEROSOL: total optical depth "//&
Index: trunk/LMDZ.GENERIC/libf/phystd/physiq_mod.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/physiq_mod.F90	(revision 2822)
+++ trunk/LMDZ.GENERIC/libf/phystd/physiq_mod.F90	(revision 2831)
@@ -985,4 +985,20 @@
                               clearsky,firstcall,lastcall)     
 
+                !GG (feb2021): Option to "artificially" decrease the raditive time scale in
+                !the deep atmosphere  press > 0.1 bar. Suggested by MT.
+                !! COEFF_RAD to be "tuned" to facilitate convergence of tendency 
+        
+                !coeff_rad=0.   ! 0 values, it doesn't accelerate the convergence
+                !coeff_rad=0.5
+                !coeff_rad=1.                 
+                !do l=1, nlayer
+                !  do ig=1,ngrid
+                !    if(pplay(ig,l).ge.1.d4) then
+                !      zdtsw(ig,l)=zdtsw(ig,l)*(pplay(ig,l)/1.d4)**coeff_rad
+                !      zdtlw(ig,l)=zdtlw(ig,l)*(pplay(ig,l)/1.d4)**coeff_rad
+                !    endif
+                !  enddo
+                !enddo
+
                 ! Option to call scheme once more for clear regions  
                if(CLFvarying)then
@@ -2467,8 +2483,8 @@
 
         ! Temporary inclusions for heating diagnostics.
-        !call writediagfi(ngrid,"zdtsw","SW heating","T s-1",3,zdtsw)
-        !call writediagfi(ngrid,"zdtlw","LW heating","T s-1",3,zdtlw)
+        call writediagfi(ngrid,"zdtsw","SW heating","T s-1",3,zdtsw)
+        call writediagfi(ngrid,"zdtlw","LW heating","T s-1",3,zdtlw)
         call writediagfi(ngrid,"dtrad","radiative heating","K s-1",3,dtrad)
-        ! call writediagfi(ngrid,"zdtdyn","Dyn. heating","T s-1",3,zdtdyn)
+        call writediagfi(ngrid,"zdtdyn","Dyn. heating","T s-1",3,zdtdyn)
         
         ! For Debugging.
Index: trunk/LMDZ.GENERIC/libf/phystd/radii_mod.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/radii_mod.F90	(revision 2822)
+++ trunk/LMDZ.GENERIC/libf/phystd/radii_mod.F90	(revision 2831)
@@ -38,10 +38,14 @@
 !
 !==================================================================
+      use mod_phys_lmdz_para, only : is_master
       use ioipsl_getin_p_mod, only: getin_p
       use radinc_h, only: naerkind
       use aerosol_mod, only: iaero_back2lay, iaero_co2, iaero_dust, &
-                             iaero_h2o, iaero_h2so4, iaero_nh3, iaero_nlay, iaero_aurora, &
-                             iaero_generic, i_rgcs_ice
-      use callkeys_mod, only: size_nh3_cloud, nlayaero, aeronlay_size, aeronlay_nueff,aerogeneric
+                             iaero_h2o, iaero_h2so4, iaero_nh3, iaero_nlay, &
+                             iaero_aurora, iaero_generic, i_rgcs_ice, &
+                             iaero_venus1, iaero_venus2, iaero_venus2p, &
+                             iaero_venus3, iaero_venusUV
+      use callkeys_mod, only: size_nh3_cloud, nlayaero, aeronlay_size, &
+                              aeronlay_nueff,aerogeneric
       use tracer_h, only: radius, nqtot, is_rgcs
       Implicit none
@@ -90,5 +94,5 @@
 
 
-	      if(iaer.eq.iaero_nh3)then ! Nh3 cloud
+	 if(iaer.eq.iaero_nh3)then ! Nh3 cloud
             reffrad(1:ngrid,1:nlayer,iaer) = size_nh3_cloud
             nueffrad(1:ngrid,1:nlayer,iaer) = 0.1 
@@ -102,7 +106,32 @@
          enddo
 
-	      if(iaer.eq.iaero_aurora)then ! Auroral aerosols
+	 if(iaer.eq.iaero_aurora)then ! Auroral aerosols
             reffrad(1:ngrid,1:nlayer,iaer) = 3.e-7
             nueffrad(1:ngrid,1:nlayer,iaer) = 0.1 
+         endif
+
+         if(iaer.eq.iaero_venus1)then ! Venus cloud, mode 1, Haus13 model
+            reffrad(1:ngrid,1:nlayer,iaer)  = 0.49e-6
+            nueffrad(1:ngrid,1:nlayer,iaer) = 0.21
+         endif
+
+         if(iaer.eq.iaero_venus2)then ! Venus cloud, mode 2, Haus13 model
+            reffrad(1:ngrid,1:nlayer,iaer)  = 1.23e-6
+            nueffrad(1:ngrid,1:nlayer,iaer) = 0.067
+         endif
+
+         if(iaer.eq.iaero_venus2p)then ! Venus cloud, mode 2p, Haus13 model
+            reffrad(1:ngrid,1:nlayer,iaer)  = 1.56e-6
+            nueffrad(1:ngrid,1:nlayer,iaer) = 0.044
+         endif
+
+         if(iaer.eq.iaero_venus3)then ! Venus cloud, mode 3, Haus13 model
+            reffrad(1:ngrid,1:nlayer,iaer)  = 4.25e-6
+            nueffrad(1:ngrid,1:nlayer,iaer) = 0.062
+         endif
+
+         if(iaer.eq.iaero_venusUV)then ! Venus cloud, UV abs, 1 val as in table 
+            reffrad(1:ngrid,1:nlayer,iaer)  = 0.5e-6
+            nueffrad(1:ngrid,1:nlayer,iaer) = 0.1
          endif
 
@@ -115,10 +144,4 @@
          enddo  ! generic radiative condensable aerosols
          
-         if(iaer.gt.5)then ! I think this could be commented out by this needs testing. It works with 5 aerosols (LT 2022)
-            print*,'Error in callcorrk, naerkind is too high (>5).'
-            print*,'The code still needs generalisation to arbitrary'
-            print*,'aerosol kinds and number.'
-            call abort
-         endif
       enddo ! iaer=1,naerkind
       
@@ -126,25 +149,25 @@
       if (radfixed) then
 
-         write(*,*)"radius of H2O water particles:"
+         if (is_master) write(*,*)"radius of H2O water particles:"
          rad_h2o=13. ! default value
          call getin_p("rad_h2o",rad_h2o)
-         write(*,*)" rad_h2o = ",rad_h2o
-
-         write(*,*)"radius of H2O ice particles:"
+         if (is_master) write(*,*)" rad_h2o = ",rad_h2o
+
+         if (is_master) write(*,*)"radius of H2O ice particles:"
          rad_h2o_ice=35. ! default value
          call getin_p("rad_h2o_ice",rad_h2o_ice)
-         write(*,*)" rad_h2o_ice = ",rad_h2o_ice
+         if (is_master) write(*,*)" rad_h2o_ice = ",rad_h2o_ice
 
       else
 
-         write(*,*)"Number mixing ratio of H2O water particles:"
+         if (is_master) write(*,*)"Number mixing ratio of H2O water particles:"
          Nmix_h2o=1.e6 ! default value
          call getin_p("Nmix_h2o",Nmix_h2o)
-         write(*,*)" Nmix_h2o = ",Nmix_h2o
-
-         write(*,*)"Number mixing ratio of H2O ice particles:"
+         if (is_master) write(*,*)" Nmix_h2o = ",Nmix_h2o
+
+         if (is_master) write(*,*)"Number mixing ratio of H2O ice particles:"
          Nmix_h2o_ice=Nmix_h2o ! default value
          call getin_p("Nmix_h2o_ice",Nmix_h2o_ice)
-         write(*,*)" Nmix_h2o_ice = ",Nmix_h2o_ice
+         if (is_master) write(*,*)" Nmix_h2o_ice = ",Nmix_h2o_ice
       endif
 
Index: trunk/LMDZ.GENERIC/libf/phystd/setspv.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/setspv.F90	(revision 2822)
+++ trunk/LMDZ.GENERIC/libf/phystd/setspv.F90	(revision 2831)
@@ -130,5 +130,5 @@
          sum         = sum+STELLARF(N)
       end do
-      write(6,'("setspv: Stellar flux at 1 AU = ",f7.2," W m-2")') sum
+      write(6,'("setspv: Stellar flux at 1 AU = ",f9.2," W m-2")') sum
       print*,' '
 
Index: trunk/LMDZ.GENERIC/libf/phystd/su_gases.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/su_gases.F90	(revision 2822)
+++ trunk/LMDZ.GENERIC/libf/phystd/su_gases.F90	(revision 2831)
@@ -109,4 +109,14 @@
            igas_C2H2=igas
            count=count+1
+        !! GG MODIF 11 Jan 2019
+       elseif (trim(gnom(igas)).eq."OCS") then
+           igas_OCS=igas
+           count=count+1
+       elseif (trim(gnom(igas)).eq."HCl") then
+           igas_HCl=igas
+           count=count+1
+       elseif (trim(gnom(igas)).eq."HF") then
+           igas_HF=igas
+           count=count+1
         endif
      enddo
Index: trunk/LMDZ.GENERIC/libf/phystd/suaer_corrk.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/suaer_corrk.F90	(revision 2822)
+++ trunk/LMDZ.GENERIC/libf/phystd/suaer_corrk.F90	(revision 2831)
@@ -232,4 +232,69 @@
 ! added by SG
        endif
+
+! VENUS CLOUDS
+
+       if (iaer.eq.iaero_venus1) then
+         print*, 'naerkind= venus1', iaer
+
+!     visible
+         file_id(iaer,1) = 'optprop_h2so4vis_n50.dat'
+         lamrefvis(iaer)=1.5E-6   ! no idea, must find
+!     infrared
+         file_id(iaer,2) = 'optprop_h2so4ir_n50.dat'
+         lamrefir(iaer)=9.3E-6 ! no idea, must find
+! added by SL
+       endif
+
+       if (iaer.eq.iaero_venus2) then
+         print*, 'naerkind= venus2', iaer
+
+!     visible
+         file_id(iaer,1) = 'optprop_h2so4vis_n50.dat'
+         lamrefvis(iaer)=1.5E-6   ! no idea, must find
+!     infrared
+         file_id(iaer,2) = 'optprop_h2so4ir_n50.dat'
+         lamrefir(iaer)=9.3E-6 ! no idea, must find
+! added by SL
+       endif
+
+       if (iaer.eq.iaero_venus2p) then
+         print*, 'naerkind= venus2p', iaer
+
+!     visible
+         file_id(iaer,1) = 'optprop_h2so4vis_n50.dat'
+         lamrefvis(iaer)=1.5E-6   ! no idea, must find
+!     infrared
+         file_id(iaer,2) = 'optprop_h2so4ir_n50.dat'
+         lamrefir(iaer)=9.3E-6 ! no idea, must find
+! added by SL
+       endif
+
+       if (iaer.eq.iaero_venus3) then
+         print*, 'naerkind= venus3', iaer
+
+!     visible
+         file_id(iaer,1) = 'optprop_h2so4vis_n50.dat'
+         lamrefvis(iaer)=1.5E-6   ! no idea, must find
+!     infrared
+         file_id(iaer,2) = 'optprop_h2so4ir_n50.dat'
+         lamrefir(iaer)=9.3E-6 ! no idea, must find
+! added by SL
+       endif
+
+       if (iaer.eq.iaero_venusUV) then
+         print*, 'naerkind= venusUV', iaer
+
+!     visible
+         file_id(iaer,1) = 'optprop_venusUVvis.dat'
+         lamrefvis(iaer)=3.5E-7   ! Haus et al. 2015
+!     infrared
+         file_id(iaer,2) = 'optprop_venusUVir.dat'
+         lamrefir(iaer)=9.3E-6 ! not used anyway
+! added by SL
+       endif
+
+! END VENUS CLOUDS
+       
 ! the following was added by LT
        do ia=1,aerogeneric ! Read Radiative Generic Condensable Species data
@@ -265,5 +330,5 @@
          endif
        enddo ! ia=1,aerogeneric
-      enddo
+      enddo ! of do iaer=1,naerkind
       
 !------------------------------------------------------------------
