Index: codes/icosagcm/devel/src/initial/etat0_venus.f90
===================================================================
--- codes/icosagcm/devel/src/initial/etat0_venus.f90	(revision 739)
+++ codes/icosagcm/devel/src/initial/etat0_venus.f90	(revision 741)
@@ -3,12 +3,16 @@
   IMPLICIT NONE
   PRIVATE
+  SAVE
 
   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
+  REAL(rstd), ALLOCATABLE :: temp_eq_packed(:,:)
+
+  REAL(rstd) :: kfrict, kv
+!$OMP THREADPRIVATE(kfrict, kv)
+
+  REAL(rstd), PARAMETER :: tauCLee=86400*25 ! 25 Earth days, cf Lebonnois 2012
+
+  PUBLIC :: etat0, init_physics, full_physics
   
 CONTAINS
@@ -16,16 +20,17 @@
 !-------------------------------- "Physics" ----------------------------------------
 
-  SUBROUTINE init_physics
+  SUBROUTINE init_physics_old
     USE getin_mod
-    IMPLICIT NONE
     REAL(rstd),POINTER :: temp(:,:)
     REAL(rstd) :: friction_time
     INTEGER :: ind
 
+    kv=0.15  ! vertical turbulent viscosity
+    CALL getin('venus_diffusion',kv)
     friction_time=86400.  !friction 
-    CALL getin('friction_time',friction_time)
+    CALL getin('venus_friction_time',friction_time)
     kfrict=1./friction_time
 
-    CALL allocate_field(f_temp,field_t,type_real,llm) ! Buffer for later use
+    CALL allocate_field(f_temp,field_t,type_real,llm) ! Buffer for later use by physics
 
     PRINT *, 'Initializing Temp_eq (venus)'
@@ -36,12 +41,10 @@
        CALL swap_geometry(ind)
        temp=f_temp_eq(ind)
-       CALL compute_temp_ref(temp, .TRUE.) ! FIXME With meridional gradient
-    ENDDO
-  END SUBROUTINE init_physics
+       CALL compute_temp_ref(.TRUE.,iim*jjm,lat_i, temp) ! With meridional gradient
+    ENDDO
+  END SUBROUTINE init_physics_old
 
   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(:)
@@ -66,11 +69,7 @@
 
   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
@@ -85,10 +84,128 @@
   END SUBROUTINE compute_physics
 
+!----------------------- Re-implementation using physics_interface_mod -----------------
+
+  SUBROUTINE init_physics
+    USE getin_mod
+    USE physics_interface_mod
+    REAL(rstd),POINTER :: temp(:,:)
+    REAL(rstd) :: friction_time
+    INTEGER :: ngrid
+
+    kv=0.15  ! vertical turbulent viscosity
+    CALL getin('venus_diffusion',kv)
+    friction_time=86400.  !friction 
+    CALL getin('venus_friction_time',friction_time)
+    kfrict=1./friction_time
+
+    ngrid = physics_inout%ngrid
+    ALLOCATE(temp_eq_packed(ngrid,llm))
+    CALL compute_temp_ref(.TRUE.,ngrid,physics_inout%lat, temp_eq_packed) ! With meridional gradient
+  END SUBROUTINE init_physics
+
+  SUBROUTINE full_physics
+    USE physics_interface_mod
+    CALL compute_physics_column(physics_inout%ngrid, physics_inout%dt_phys, &
+         physics_inout%p, physics_inout%geopot, &
+         physics_inout%Temp, physics_inout%ulon, physics_inout%ulat, &
+         physics_inout%dTemp, physics_inout%dulon, physics_inout%dulat)
+  END SUBROUTINE full_physics
+
+  SUBROUTINE compute_physics_column(ngrid,dt_phys,p,geopot,Temp,u,v,dTemp,du,dv)
+    USE earth_const, only : g
+    INTEGER, INTENT(IN)   :: ngrid
+    REAL(rstd),INTENT(IN) :: dt_phys, &
+         p(ngrid,llm+1), geopot(ngrid,llm+1), &
+         Temp(ngrid,llm), u(ngrid,llm), v(ngrid,llm)
+    REAL(rstd),INTENT(OUT) :: dTemp(ngrid,llm), du(ngrid,llm), dv(ngrid,llm)
+    REAL(rstd) :: mass(ngrid,llm), &  ! rho.dz
+                  A(ngrid,llm+1),  &  ! off-diagonal coefficients
+                  B(ngrid,llm),    &  ! diagonal coefficients
+                  Ru(ngrid,llm), Rv(ngrid,llm), & ! right-hand-sides
+                  C(ngrid,llm),    &  ! LU factorization (Thomas algorithm)
+                  xu(ngrid,llm), xv(ngrid, llm)      ! solution
+    REAL(rstd) :: rho, X_ij
+    INTEGER :: l,ij
+    ! Vertical diffusion : rho.du/dt = d/dz (rho*kappa*du/dz) 
+    ! rho.dz = mass.deta
+    ! => mass.du/dt = d/deta ((rho^2 kappa/m)du/deta)
+    ! Backward Euler : (M-S)u_new = M.u_old
+    ! with M.u = mass(l)*u(l)
+    ! and  S.u = A(l+1)*(u(l+1)-u(l)) - A(l)*(u(l)-u(l-1))
+    ! => solve -A(l)u(l-1) + B(l)u(l) - A(l+1)u(l+1) = Ru(l) using Thomas algorithm
+    ! with A=tau*kappa*rho^2/m and B(l) = mass(l)+A(l)+A(l+1)
+    DO l=1,llm
+       !DIR$ SIMD
+       DO ij=1,ngrid
+          mass(ij,l) = (p(ij,l)-p(ij,l+1))*(1./g)
+          rho = (p(ij,l)-p(ij,l+1))/(geopot(ij,l+1)-geopot(ij,l))
+          ! A = kappa.tau.rho^2/m
+          A(ij,l) = (kv*dt_phys)*(rho**2)/mass(ij,l)
+          Ru(ij,l) = mass(ij,l)*u(ij,l)
+          Rv(ij,l) = mass(ij,l)*v(ij,l)
+       ENDDO
+    ENDDO
+    A(:,llm+1)=0.
+    DO l=llm,2,-1
+       !DIR$ SIMD
+       DO ij=1,ngrid
+          A(ij,l) = .5*(A(ij,l)+A(ij,l-1)) ! average A to interfaces
+       ENDDO
+    ENDDO
+    A(:,1)=0.
+    DO l=1,llm
+       !DIR$ SIMD
+       DO ij=1,ngrid
+          B(ij,l) = mass(ij,l)+A(ij,l)+A(ij,l+1)
+       ENDDO
+    ENDDO
+
+    ! Solve -A(l)x(l-1) + B(l)x(l) - A(l+1)x(l+1) = R(l) using Thomas algorithm
+    ! Forward sweep :
+    ! C(0)=0, C(l) = -A(l+1) / (B(l)+A(l)C(l-1)),
+    ! D(0)=0, D(l) = (R(l)+A(l)D(l-1)) / (B(l)+A(l)C(l-1))
+    !DIR$ SIMD
+    DO ij=1,ngrid
+       X_ij = 1./B(ij,1)
+       C(ij,2) = -A(ij,2) * X_ij
+       xu(ij,1) = Ru(ij,1) * X_ij
+       xv(ij,1) = Rv(ij,1) * X_ij
+    END DO
+    DO l = 2,llm
+       !DIR$ SIMD
+       DO ij=1,ngrid
+          X_ij = 1./( B(ij,l) + A(ij,l)*C(ij,l) )
+          C(ij,l+1) = -A(ij,l+1) * X_ij ! zero for l=llm
+          xu(ij,l) = (Ru(ij,l)+A(ij,l)*xu(ij,l-1)) * X_ij
+          xv(ij,l) = (Rv(ij,l)+A(ij,l)*xv(ij,l-1)) * X_ij
+       END DO
+    END DO
+    ! Back substitution :
+    ! x(i) = D(i)-C(i+1)x(i+1), x(llm)=D(llm)
+    !DIR$ SIMD
+    DO ij=1,ngrid
+       ! top layer l=llm
+       du(ij,llm) = (xu(ij,llm)-u(ij,llm))/dt_phys
+       dv(ij,llm) = (xv(ij,llm)-v(ij,llm))/dt_phys
+    END DO
+    ! Back substitution at lower layers
+    DO l = llm-1,1,-1
+       !DIR$ SIMD
+       DO ij=1,ngrid
+          xu(ij,l) = xu(ij,l) - C(ij,l+1)*xu(ij,l+1)
+          xv(ij,l) = xv(ij,l) - C(ij,l+1)*xv(ij,l+1)
+          du(ij,l) = (xu(ij,l)-u(ij,l))/dt_phys
+          dv(ij,l) = (xv(ij,l)-v(ij,l))/dt_phys
+       END DO
+    END DO
+    ! bottom friction + thermal relaxation
+    du(:,1)=du(:,1)-kfrict*u(:,1)
+    dTemp(:,:) = (1./tauCLee)*(temp_eq_packed(:,:)-temp(:,:))
+  END SUBROUTINE compute_physics_column
+
 !----------------------------- 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(:)
@@ -103,9 +220,7 @@
     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)
+    CALL allocate_field(f_temp,field_t,type_real,llm, name='temp_buf_venus')
     DO ind=1,ndomain
        IF (.NOT. assigned_domain(ind)) CYCLE
@@ -121,5 +236,5 @@
        q(:,:,:)=1e2
        temp=f_temp(ind)
-       CALL compute_temp_ref(temp, .FALSE.) ! Without meridional gradient
+       CALL compute_temp_ref(.FALSE., iim*jjm, lat_i, temp) ! Without meridional gradient
     ENDDO
     CALL temperature2theta_rhodz(f_ps,f_temp,f_theta_rhodz)
@@ -130,19 +245,20 @@
 !------------------------- Compute reference temperature field ------------------------
 
-  SUBROUTINE compute_temp_ref(theta_eq, gradient)
-    USE icosa
+  SUBROUTINE compute_temp_ref(gradient,ngrid,lat, temp_eq)
     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)
+    INTEGER, INTENT(IN) :: ngrid
+    REAL(rstd), INTENT(IN) :: lat(ngrid) ! latitude
+    REAL(rstd), INTENT(OUT) :: temp_eq(ngrid,llm) 
+
+    REAL(rstd)  :: clat(ngrid)        
+    INTEGER :: level
+    INTEGER, PARAMETER :: nlev=30
+    REAL ::  pressCLee(nlev+1), tempCLee(nlev+1), dt_epCLee(nlev+1), etaCLee(nlev+1)
     
-    real(rstd) ::  lon,lat, pplay, ztemp,zdt,fact
-    logical, save ::  firstcall
-    integer :: i,j,ij, l,ll
+    REAL(rstd) ::  pplay, ztemp,zdt,fact
+    INTEGER :: ij, l,ll
     
     data etaCLee / 9.602e-1, 8.679e-1, 7.577e-1, 6.420e-1, 5.299e-1, &
@@ -154,11 +270,4 @@
                       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, &
@@ -175,35 +284,32 @@
                      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
+    pressCLee(:) = etaCLee(:)*9.2e6
+    clat(:)=COS(lat(:))
+
+    DO ij=1,ngrid
+       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, nlev ! nlev 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
+             temp_eq(ij,l) = ztemp+ zdt*(clat(ij)-Pi/4.)
+          ELSE
+             temp_eq(ij,l) = ztemp
+          END IF
        END DO
     END DO
Index: codes/icosagcm/devel/src/physics/physics.f90
===================================================================
--- codes/icosagcm/devel/src/physics/physics.f90	(revision 739)
+++ codes/icosagcm/devel/src/physics/physics.f90	(revision 741)
@@ -12,5 +12,4 @@
        phys_lmdz_generic=21, phys_external=22
   INTEGER :: phys_type
-  TYPE(t_field),POINTER,SAVE :: f_extra_physics_2D(:), f_extra_physics_3D(:)
   TYPE(t_field),POINTER,SAVE :: f_dulon(:), f_dulat(:)
   TYPE(t_field),POINTER,SAVE :: f_ulon(:), f_ulat(:)
@@ -50,7 +49,4 @@
     CASE ('held_suarez')
        phys_type = phys_HS94
-    CASE ('Lebonnois2012')
-       phys_type = phys_LB2012
-       CALL init_phys_venus       
     CASE ('phys_lmdz_generic')
        phys_type=phys_lmdz_generic
@@ -79,4 +75,7 @@
           phys_type = phys_DCMIP2016
           CALL init_physics_dcmip2016
+       CASE ('Lebonnois2012')
+          phys_type = phys_LB2012
+          CALL init_phys_venus       
        CASE DEFAULT
           IF(is_mpi_root) PRINT*, 'init_physics : Bad selector for variable physics <',&
@@ -129,5 +128,5 @@
   END SUBROUTINE add_du_phys
 
-  SUBROUTINE physics(it,f_phis, f_ps, f_theta_rhodz, f_ue, f_wflux, f_q)
+  SUBROUTINE physics(it,f_phis, f_geopot, f_ps, f_theta_rhodz, f_ue, f_wflux, f_q)
     USE physics_lmdz_generic_mod, ONLY : physics_lmdz_generic => physics
     USE physics_external_mod, ONLY : physics_external => physics
@@ -135,7 +134,7 @@
     USE physics_dcmip2016_mod, ONLY : write_physics_dcmip2016 => write_physics
     USE etat0_heldsz_mod
-    USE etat0_venus_mod, ONLY : phys_venus => physics
     INTEGER, INTENT(IN)   :: it
     TYPE(t_field),POINTER :: f_phis(:)
+    TYPE(t_field),POINTER :: f_geopot(:)
     TYPE(t_field),POINTER :: f_ps(:)
     TYPE(t_field),POINTER :: f_theta_rhodz(:)
@@ -161,8 +160,6 @@
        CASE (phys_external)
          CALL physics_external(it ,f_phis, f_ps, f_theta_rhodz, f_ue, f_wflux, f_q)
-       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)
+          CALL physics_column(it, f_phis, f_geopot, f_ps, f_theta_rhodz, f_ue, f_q)
        END SELECT
 
@@ -197,11 +194,13 @@
   END SUBROUTINE write_physics_tendencies
     
-  SUBROUTINE physics_column(it, f_phis, f_ps, f_theta_rhodz, f_ue, f_q)
+  SUBROUTINE physics_column(it, f_phis, f_geopot, f_ps, f_theta_rhodz, f_ue, f_q)
     USE physics_dcmip_mod, ONLY : full_physics_dcmip => full_physics
     USE physics_dcmip2016_mod, ONLY : full_physics_dcmip2016 => full_physics
+    USE etat0_venus_mod, ONLY : full_physics_venus=>full_physics
     USE theta2theta_rhodz_mod
     USE mpipara
     USE checksum_mod
     TYPE(t_field),POINTER :: f_phis(:)
+    TYPE(t_field),POINTER :: f_geopot(:)
     TYPE(t_field),POINTER :: f_ps(:)
     TYPE(t_field),POINTER :: f_theta_rhodz(:)
@@ -209,4 +208,5 @@
     TYPE(t_field),POINTER :: f_q(:)
     REAL(rstd),POINTER :: phis(:)
+    REAL(rstd),POINTER :: geopot(:,:)
     REAL(rstd),POINTER :: ps(:)
     REAL(rstd),POINTER :: temp(:,:)
@@ -228,4 +228,5 @@
        CALL swap_geometry(ind)
        phis=f_phis(ind)
+       geopot=f_geopot(ind)
        ps=f_ps(ind)
        temp=f_temp(ind)
@@ -236,5 +237,5 @@
        ulon=f_ulon(ind)
        ulat=f_ulat(ind)
-       CALL pack_physics(pack_info(ind), phis, ps, temp, ue, q, p, pk, ulon, ulat)
+       CALL pack_physics(pack_info(ind), phis, geopot, ps, temp, ue, q, p, pk, ulon, ulat)
     END DO
 
@@ -244,4 +245,6 @@
     CASE (phys_DCMIP2016)
        IF (is_omp_level_master) CALL full_physics_dcmip2016
+    CASE(phys_LB2012)
+       IF (is_omp_level_master) CALL full_physics_venus
     CASE DEFAULT
        IF(is_master) PRINT *,'Internal error : illegal value of phys_type', phys_type
@@ -279,5 +282,5 @@
   END SUBROUTINE physics_column
 
-  SUBROUTINE pack_physics(info, phis, ps, temp, ue, q, p, pk, ulon, ulat )
+  SUBROUTINE pack_physics(info, phis, geopot, ps, temp, ue, q, p, pk, ulon, ulat )
     USE wind_mod
     USE pression_mod
@@ -286,4 +289,5 @@
     TYPE(t_pack_info) :: info
     REAL(rstd) :: phis(iim*jjm)
+    REAL(rstd) :: geopot(iim*jjm,llm+1)
     REAL(rstd) :: ps(iim*jjm)
     REAL(rstd) :: temp(iim*jjm,llm)
@@ -308,4 +312,5 @@
     IF (is_omp_level_master) THEN
       CALL pack_domain(info, phis, physics_inout%phis)
+      CALL pack_domain(info, geopot, physics_inout%geopot)
       CALL pack_domain(info, p, physics_inout%p)
       CALL pack_domain(info, pk, physics_inout%pk)
Index: codes/icosagcm/devel/src/physics/physics_interface.f90
===================================================================
--- codes/icosagcm/devel/src/physics/physics_interface.f90	(revision 739)
+++ codes/icosagcm/devel/src/physics/physics_interface.f90	(revision 741)
@@ -11,5 +11,5 @@
      REAL(rstd), DIMENSION(:), POINTER :: Ai, lon, lat, phis
      ! Input, time-dependent
-     REAL(rstd), DIMENSION(:,:), POINTER :: p, pk, Temp, ulon, ulat
+     REAL(rstd), DIMENSION(:,:), POINTER :: geopot, p, pk, Temp, ulon, ulat
      REAL(rstd), DIMENSION(:,:,:), POINTER :: q
      ! Output arrays
@@ -93,5 +93,5 @@
 
     ngrid=offset
-    ! Input                                                                                                     
+    ! Input
     ALLOCATE(physics_inout%Ai(ngrid))
     ALLOCATE(physics_inout%lon(ngrid))
@@ -99,4 +99,5 @@
     ALLOCATE(physics_inout%phis(ngrid))
     ALLOCATE(physics_inout%p(ngrid,llm+1))
+    ALLOCATE(physics_inout%geopot(ngrid,llm+1))
     ALLOCATE(physics_inout%pk(ngrid,llm))
     ALLOCATE(physics_inout%Temp(ngrid,llm))
Index: codes/icosagcm/devel/src/time/timeloop_gcm.f90
===================================================================
--- codes/icosagcm/devel/src/time/timeloop_gcm.f90	(revision 739)
+++ codes/icosagcm/devel/src/time/timeloop_gcm.f90	(revision 741)
@@ -350,5 +350,5 @@
             f_ps,f_dps,f_u,f_theta_rhodz,f_phis)
           CALL enter_profile(id_phys)
-          CALL physics(it,f_phis, f_ps, f_theta_rhodz, f_u, f_wflux, f_q)        
+          CALL physics(it, f_phis, f_geopot, f_ps, f_theta_rhodz, f_u, f_wflux, f_q)
           CALL exit_profile(id_phys)
           CALL check_conserve_detailed(it, AAM_phys, &
