Index: codes/icosagcm/trunk/src/etat0.f90
===================================================================
--- codes/icosagcm/trunk/src/etat0.f90	(revision 324)
+++ codes/icosagcm/trunk/src/etat0.f90	(revision 325)
@@ -25,4 +25,5 @@
     USE etat0_dcmip4_mod, ONLY : etat0_dcmip4=>etat0  
     USE etat0_heldsz_mod, ONLY : etat0_heldsz=>etat0  
+    USE etat0_venus_mod,  ONLY : etat0_venus=>etat0  
     USE etat0_start_file_mod, ONLY : etat0_start_file=>etat0  
 
@@ -72,4 +73,7 @@
        PRINT *,"Held & Suarez (1994) test case"
        CALL etat0_heldsz(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
+    CASE ('venus')
+       CALL etat0_venus(f_ps, f_phis, f_theta_rhodz, f_u, f_q)
+       PRINT *, "Venus (Lebonnois et al., 2012) test case"
     CASE ('dcmip1')
        CALL etat0_dcmip1(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
Index: codes/icosagcm/trunk/src/etat0_venus.f90
===================================================================
--- codes/icosagcm/trunk/src/etat0_venus.f90	(revision 325)
+++ codes/icosagcm/trunk/src/etat0_venus.f90	(revision 325)
@@ -0,0 +1,212 @@
+MODULE etat0_venus_mod
+  USE icosa
+  IMPLICIT NONE
+  PRIVATE
+
+  TYPE(t_field),POINTER :: f_temp_eq( :)
+  TYPE(t_field),POINTER :: f_temp(:) ! buffer used for physics
+
+  REAL(rstd), SAVE :: kfrict
+!$OMP THREADPRIVATE(kfrict)
+
+  PUBLIC :: etat0, init_physics, physics
+  
+CONTAINS
+
+!-------------------------------- "Physics" ----------------------------------------
+
+  SUBROUTINE init_physics
+    USE getin_mod
+    IMPLICIT NONE
+    REAL(rstd),POINTER :: temp(:,:)
+    REAL(rstd) :: friction_time
+    INTEGER :: ind
+
+    friction_time=86400.  !friction 
+    CALL getin('friction_time',friction_time)
+    kfrict=1./friction_time
+
+    CALL allocate_field(f_temp,field_t,type_real,llm) ! Buffer for later use
+
+    PRINT *, 'Initializing Temp_eq (venus)'
+    CALL allocate_field(f_temp_eq,field_t,type_real,llm)
+    DO ind=1,ndomain
+       IF (.NOT. assigned_domain(ind)) CYCLE
+       CALL swap_dimensions(ind)
+       CALL swap_geometry(ind)
+       temp=f_temp_eq(ind)
+       CALL compute_temp_ref(temp, .TRUE.) ! FIXME With meridional gradient
+    ENDDO
+  END SUBROUTINE init_physics
+
+  SUBROUTINE physics(f_ps,f_theta_rhodz,f_u) 
+    USE icosa
+    USE theta2theta_rhodz_mod
+    IMPLICIT NONE
+    TYPE(t_field),POINTER :: f_theta_rhodz(:)
+    TYPE(t_field),POINTER :: f_u(:)
+    TYPE(t_field),POINTER :: f_ps(:)
+    REAL(rstd),POINTER :: temp(:,:)
+    REAL(rstd),POINTER :: temp_eq(:,:)
+    REAL(rstd),POINTER :: u(:,:)
+    INTEGER :: ind
+
+    CALL theta_rhodz2temperature(f_ps,f_theta_rhodz,f_temp)
+    DO ind=1,ndomain
+       IF (.NOT. assigned_domain(ind)) CYCLE
+       CALL swap_dimensions(ind)
+       CALL swap_geometry(ind)
+       u=f_u(ind)
+       temp_eq=f_temp_eq(ind) 
+       temp=f_temp(ind) 
+       CALL compute_physics(temp_eq, temp, u)
+    ENDDO
+    CALL temperature2theta_rhodz(f_ps,f_temp,f_theta_rhodz)
+  END SUBROUTINE physics
+
+  SUBROUTINE compute_physics(temp_eq, temp, u)
+    USE icosa
+    USE theta2theta_rhodz_mod
+    IMPLICIT NONE
+    REAL(rstd),INTENT(IN)    :: temp_eq(iim*jjm,llm) 
+    REAL(rstd),INTENT(INOUT) :: temp(iim*jjm,llm)
+    REAL(rstd),INTENT(INOUT) :: u(3*iim*jjm,llm)
+    REAL(rstd), PARAMETER :: tauCLee=86400*25 ! 25 Earth days, cf Lebonnois 2012
+    INTEGER :: i,j,l,ij
+    DO l=1,llm
+       DO j=jj_begin-1,jj_end+1
+          DO i=ii_begin-1,ii_end+1
+             ij=(j-1)*iim+i
+             temp(ij,l) = temp(ij,l) - (temp(ij,l)-temp_eq(ij,l))*(dt*itau_physics/tauCLee)
+          ENDDO
+       ENDDO
+    ENDDO
+    u(:,1)=u(:,1)*(1.-dt*itau_physics*kfrict)
+  END SUBROUTINE compute_physics
+
+!----------------------------- Initialize to T_eq --------------------------------------
+
+  SUBROUTINE etat0(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
+    USE icosa
+    USE theta2theta_rhodz_mod
+    IMPLICIT NONE
+    TYPE(t_field),POINTER :: f_ps(:)
+    TYPE(t_field),POINTER :: f_phis(:)
+    TYPE(t_field),POINTER :: f_theta_rhodz(:)
+    TYPE(t_field),POINTER :: f_u(:)
+    TYPE(t_field),POINTER :: f_q(:)
+
+    TYPE(t_field),POINTER :: f_temp(:)
+    REAL(rstd),POINTER :: temp(:,:) 
+    REAL(rstd),POINTER :: ps(:)
+    REAL(rstd),POINTER :: phis(:)
+    REAL(rstd),POINTER :: u(:,:)
+    REAL(rstd),POINTER :: q(:,:,:)
+    REAL(rstd) :: lat(iim*jjm)        ! latitude                    
+    REAL(rstd) :: pplay(iim*jjm, llm) ! pressure at full layers
+    INTEGER :: ind
+
+    CALL allocate_field(f_temp,field_t,type_real,llm)
+    DO ind=1,ndomain
+       IF (.NOT. assigned_domain(ind)) CYCLE
+       CALL swap_dimensions(ind)
+       CALL swap_geometry(ind)
+       ps=f_ps(ind)
+       ps(:)=preff 
+       phis=f_phis(ind)
+       phis(:)=0.
+       u=f_u(ind)
+       u(:,:)=0
+       q=f_q(ind)
+       q(:,:,:)=1e2
+       temp=f_temp(ind)
+       CALL compute_temp_ref(temp, .FALSE.) ! Without meridional gradient
+    ENDDO
+    CALL temperature2theta_rhodz(f_ps,f_temp,f_theta_rhodz)
+    CALL deallocate_field(f_temp)
+
+  END SUBROUTINE etat0
+
+!------------------------- Compute reference temperature field ------------------------
+
+  SUBROUTINE compute_temp_ref(theta_eq, gradient)
+    USE icosa
+    USE disvert_mod
+    USE omp_para
+    USE math_const
+    IMPLICIT NONE
+    REAL(rstd), INTENT(OUT) :: theta_eq(iim*jjm,llm) 
+    LOGICAL, INTENT(IN) :: gradient
+    REAL(rstd)  :: clat(iim*jjm)        ! latitude
+    integer :: level
+    REAL ::  pressCLee(31), tempCLee(31), dt_epCLee(31), etaCLee(31)
+    
+    real(rstd) ::  lon,lat, pplay, ztemp,zdt,fact
+    logical, save ::  firstcall
+    integer :: i,j,ij, l,ll
+    
+    data etaCLee / 9.602e-1, 8.679e-1, 7.577e-1, 6.420e-1, 5.299e-1, &
+                      4.273e-1, 3.373e-1, 2.610e-1,1.979e-1,1.472e-1,  &
+                      1.074e-1, 7.672e-2, 5.361e-2,3.657e-2,2.430e-2,  &
+                      1.569e-2, 9.814e-3, 5.929e-3,3.454e-3,1.934e-3,  &
+                      1.043e-3, 5.400e-4, 2.710e-4,1.324e-4,6.355e-5,  &
+                      3.070e-5, 1.525e-5, 7.950e-6,4.500e-6,2.925e-6,  &
+                      2.265e-6/
+
+    DO j=jj_begin-1,jj_end+1
+       DO i=ii_begin-1,ii_end+1
+          ij=(j-1)*iim+i
+          CALL xyz2lonlat(xyz_i(ij,:),lon,lat)
+          clat(ij)=cos(lat) 
+       ENDDO
+    ENDDO
+
+    data   tempCLee/ 728.187, 715.129, 697.876, 677.284, 654.078, 628.885, &
+                     602.225, 574.542, 546.104, 517.339, 488.560, 459.932, &
+                     431.741, 404.202, 377.555, 352.042, 327.887, 305.313, &
+                     284.556, 265.697, 248.844, 233.771, 220.368, 208.247, &
+                     197.127, 187.104, 178.489, 171.800, 167.598, 165.899, &
+                     165.676/
+    data   dt_epCLee/6.101 , 6.136 , 6.176 , 6.410 , 6.634 , 6.678 , &
+                     6.719 , 6.762 , 7.167 , 7.524 , 9.840 ,14.948 ,  &
+                     21.370 ,28.746 ,36.373 ,43.315 ,48.534 ,51.175 ,  &
+                     50.757 ,47.342 ,41.536 ,34.295 ,26.758 ,19.807 ,  &
+                     14.001 , 9.599 , 6.504 , 4.439 , 3.126 , 2.370 ,  &
+                     2.000/
+      
+    pressCLee = etaCLee*9.2e6
+    
+    DO j=jj_begin-1,jj_end+1
+       DO i=ii_begin-1,ii_end+1
+          ij=(j-1)*iim+i
+
+          DO l = 1, llm
+
+             pplay = .5*(ap(l)+ap(l+1)+(bp(l)+bp(l+1))*preff) ! ps=preff          
+             ! look for largest level such that pressCLee(level) > pplay(ij,l))  
+             ! => pressClee(level+1) < pplay(ij,l) < pressClee(level)
+             level = 1
+             DO ll = 1, 30 ! 30 data levels
+                IF(pressCLee(ll) > pplay) THEN
+                   level = ll
+                END IF
+             END DO
+             
+             ! interpolate between level and level+1
+             ! interpolation is linear in log(pressure)
+             fact  = ( log10(pplay)-log10(pressCLee(level))) &
+                    /( log10(pressCLee(level+1))-log10(pressCLee(level)) )
+             ztemp = tempCLee(level)*(1-fact) + tempCLee(level+1)*fact
+             zdt   = dt_epCLee(level)*(1-fact) + dt_epCLee(level+1)*fact
+
+             IF(gradient) THEN
+                theta_eq(ij,l) = ztemp+ zdt*(clat(ij)-Pi/4.)
+             ELSE
+                theta_eq(ij,l) = ztemp
+             END IF
+          END DO
+       END DO
+    END DO
+  END SUBROUTINE compute_temp_ref
+
+END MODULE etat0_venus_mod
Index: codes/icosagcm/trunk/src/physics.f90
===================================================================
--- codes/icosagcm/trunk/src/physics.f90	(revision 324)
+++ codes/icosagcm/trunk/src/physics.f90	(revision 325)
@@ -5,5 +5,5 @@
   PRIVATE
 
-  INTEGER, PARAMETER :: phys_none=0, phys_HS94=1, phys_DCMIP=2
+  INTEGER, PARAMETER :: phys_none=0, phys_HS94=1, phys_DCMIP=2, phys_LB2012=3
 
   INTEGER :: phys_type
@@ -25,4 +25,5 @@
     USE physics_interface_mod
     USE physics_dcmip_mod, ONLY : init_physics_dcmip=>init_physics
+    USE etat0_venus_mod, ONLY : init_phys_venus=>init_physics
     IMPLICIT NONE
 
@@ -36,4 +37,7 @@
     CASE ('held_suarez')
        phys_type = phys_HS94
+    CASE ('Lebonnois2012')
+       phys_type = phys_LB2012
+       CALL init_phys_venus
     CASE ('dcmip')
        CALL allocate_field(f_dulon,field_t,type_real,llm, name='dulon')
@@ -46,5 +50,5 @@
     CASE DEFAULT
        IF(is_mpi_root) PRINT*, 'init_physics : Bad selector for variable physics <',&
-            TRIM(physics_type), '> options are <none>, <held_suarez>, <dcmip>'
+            TRIM(physics_type), '> options are <none>, <held_suarez>, <Lebonnois2012>, <dcmip>'
        STOP
     END SELECT
@@ -58,4 +62,5 @@
     USE physics_dcmip_mod, ONLY : write_physics_dcmip => write_physics
     USE etat0_heldsz_mod
+    USE etat0_venus_mod, ONLY : phys_venus => physics
     IMPLICIT NONE
     INTEGER, INTENT(IN)   :: it
@@ -82,4 +87,6 @@
        CASE(phys_HS94)
           CALL held_suarez(f_ps,f_theta_rhodz,f_ue) 
+       CASE(phys_LB2012)
+          CALL phys_venus(f_ps,f_theta_rhodz,f_ue) 
        CASE DEFAULT
           CALL physics_column(it, f_phis, f_ps, f_theta_rhodz, f_ue, f_q)
