Index: codes/icosagcm/trunk/src/caldyn_kernels.f90
===================================================================
--- codes/icosagcm/trunk/src/caldyn_kernels.f90	(revision 427)
+++ codes/icosagcm/trunk/src/caldyn_kernels.f90	(revision 428)
@@ -33,5 +33,4 @@
     USE icosa
     USE disvert_mod, ONLY : mass_dak, mass_dbk, caldyn_eta, eta_mass
-    USE exner_mod
     USE trace
     USE omp_para
@@ -121,5 +120,4 @@
     USE icosa
     USE disvert_mod
-    USE exner_mod
     USE trace
     USE omp_para
Index: codes/icosagcm/trunk/src/etat0_academic.f90
===================================================================
--- codes/icosagcm/trunk/src/etat0_academic.f90	(revision 427)
+++ codes/icosagcm/trunk/src/etat0_academic.f90	(revision 428)
@@ -1,7 +1,10 @@
 MODULE etat0_academic_mod
-
-
-
-
+  USE icosa
+  IMPLICIT NONE
+
+  PRIVATE
+  
+  PUBLIC :: etat0
+  
 CONTAINS
   
@@ -9,5 +12,4 @@
   USE icosa
   USE kinetic_mod
-  IMPLICIT NONE
     TYPE(t_field),POINTER,SAVE :: f_ps(:)
     TYPE(t_field),POINTER,SAVE :: f_phis(:)
@@ -21,6 +23,5 @@
     REAL(rstd),POINTER :: temp(:)
     INTEGER :: ind
-    
-    
+        
     CALL allocate_field(f_ps,field_t,type_real)
     CALL allocate_field(f_phis,field_t,type_real)
@@ -40,9 +41,7 @@
     
   END SUBROUTINE test_etat0_academic
-   
     
   SUBROUTINE etat0(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
   USE icosa
-  IMPLICIT NONE
     TYPE(t_field),POINTER :: f_ps(:)
     TYPE(t_field),POINTER :: f_phis(:)
@@ -56,5 +55,9 @@
     REAL(rstd),POINTER :: u(:,:)
     INTEGER :: ind
-    
+
+    PRINT *, 'etat0_academic needs an upgrade for 4D theta'
+    PRINT *, 'STOP in etat0_academic.f90/etat0'
+    STOP
+
     DO ind=1,ndomain
       IF (.NOT. assigned_domain(ind)) CYCLE
@@ -77,5 +80,4 @@
   USE geopotential_mod
   USE theta2theta_rhodz_mod
-  IMPLICIT NONE  
   REAL(rstd),INTENT(OUT) :: ps(iim*jjm)
   REAL(rstd),INTENT(OUT) :: phis(iim*jjm)
@@ -99,6 +101,5 @@
   REAL(rstd) :: fact(3*iim*jjm)
   REAL(rstd) :: ut(3*iim*jjm,llm)
-  
-  
+    
     DO l=1,llm
        zsig=ap(l)/preff+bp(l)
@@ -137,7 +138,5 @@
     CALL compute_pression(ps,p,1)     
 !$OMP BARRIER
-    CALL compute_exner(ps,p,pks,pk,1)  
-!$OMP BARRIER
-    CALL compute_geopotential(phis,pks,pk,theta,phi,1)
+    CALL compute_geopotential(phis,ps,theta,phi,1)
 
 
Index: codes/icosagcm/trunk/src/etat0_heldsz.f90
===================================================================
--- codes/icosagcm/trunk/src/etat0_heldsz.f90	(revision 427)
+++ codes/icosagcm/trunk/src/etat0_heldsz.f90	(revision 428)
@@ -23,7 +23,5 @@
 
   SUBROUTINE test_etat0_heldsz
-    USE icosa
     USE kinetic_mod
-    IMPLICIT NONE
     TYPE(t_field),POINTER :: f_ps(:)
     TYPE(t_field),POINTER :: f_phis(:)
@@ -50,5 +48,4 @@
 
   SUBROUTINE etat0(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
-    USE icosa
     USE theta2theta_rhodz_mod
     IMPLICIT NONE
@@ -99,7 +96,5 @@
 
   SUBROUTINE init_Teq
-    USE icosa
     USE disvert_mod, ONLY : ap,bp
-    IMPLICIT NONE
     REAL(rstd),POINTER :: clat(:) 
     REAL(rstd),POINTER :: theta_eq(:,:) 
@@ -159,7 +154,5 @@
 
   SUBROUTINE compute_Teq(lat,theta_eq)
-    USE icosa
     USE disvert_mod
-    IMPLICIT NONE  
     REAL(rstd),INTENT(IN) :: lat(iim*jjm)
     REAL(rstd),INTENT(OUT) :: theta_eq(iim*jjm,llm) 
@@ -184,11 +177,5 @@
 
   SUBROUTINE compute_etat0_heldsz(theta_eq, theta)
-    USE icosa
     USE disvert_mod
-    USE pression_mod
-    USE exner_mod
-    USE geopotential_mod
-    USE theta2theta_rhodz_mod
-    IMPLICIT NONE  
     REAL(rstd),INTENT(IN) :: theta_eq(iim*jjm,llm) 
     REAL(rstd),INTENT(OUT) :: theta(iim*jjm,llm) 
@@ -211,6 +198,4 @@
 
   SUBROUTINE held_suarez(f_ps,f_theta_rhodz,f_u) 
-    USE icosa
-    IMPLICIT NONE
     TYPE(t_field),POINTER :: f_theta_rhodz(:)
     TYPE(t_field),POINTER :: f_u(:)
@@ -238,7 +223,5 @@
 
   SUBROUTINE compute_heldsz(ps,theta_eq,lat, theta_rhodz,u, theta) 
-    USE icosa
     USE theta2theta_rhodz_mod
-    IMPLICIT NONE
     REAL(rstd),INTENT(IN)    :: ps(iim*jjm) 
     REAL(rstd),INTENT(IN)    :: theta_eq(iim*jjm,llm) 
Index: codes/icosagcm/trunk/src/exner.f90
===================================================================
--- codes/icosagcm/trunk/src/exner.f90	(revision 427)
+++ codes/icosagcm/trunk/src/exner.f90	(revision 428)
@@ -11,8 +11,8 @@
   USE icosa
   IMPLICIT NONE
-    TYPE(t_field), POINTER :: f_ps(:)
-    TYPE(t_field), POINTER :: f_p(:)
-    TYPE(t_field), POINTER :: f_pks(:)
-    TYPE(t_field), POINTER :: f_pk(:)
+    TYPE(t_field), POINTER :: f_ps(:)  ! IN
+    TYPE(t_field), POINTER :: f_p(:)   ! IN
+    TYPE(t_field), POINTER :: f_pks(:) ! OUT
+    TYPE(t_field), POINTER :: f_pk(:)  ! OUT
   
     REAL(rstd), POINTER :: ps(:)
@@ -49,82 +49,27 @@
     INTEGER,INTENT(IN) :: offset
     INTEGER :: i,j,ij,l
-    REAL(rstd) :: alpha(iim*jjm,llm),beta(iim*jjm,llm)
-    REAL(rstd) :: delta 
     
-    IF(caldyn_exner == lmdz) THEN          ! LMD-Z style calculation of Exner pressure
-       !! Compute Alpha and Beta
+    ! surface : pks
+    IF (is_omp_level_master) THEN
        
-       IF (is_omp_level_master) THEN
-       ! for llm layer
-         DO j=jj_begin-offset,jj_end+offset
-            DO i=ii_begin-offset,ii_end+offset
-               ij=(j-1)*iim+i
-               alpha(ij,llm) = 0.
-               beta (ij,llm) = 1./ (1+ 2*kappa)
-            ENDDO
-         ENDDO
+       DO j=jj_begin-offset,jj_end+offset
+          DO i=ii_begin-offset,ii_end+offset
+             ij=(j-1)*iim+i
+             pks(ij) = cpp * ( ps(ij)/preff ) ** kappa
+          ENDDO
+       ENDDO
        
-       ! for other layer   
-         DO l = llm-1 , 2 , -1
-            DO j=jj_begin-offset,jj_end+offset
-               DO i=ii_begin-offset,ii_end+offset
-                  ij=(j-1)*iim+i
-                  delta = p(ij,l)* (1+2*kappa) + p(ij,l+1)* ( beta(ij,l+1)- (1+2*kappa) )
-                  alpha(ij,l)  = - p(ij,l+1) / delta * alpha(ij,l+1)
-                  beta (ij,l)  =   p(ij,l  ) / delta   
-               ENDDO
-            ENDDO
-         ENDDO
-       
-         !! Compute pk
-       
-         ! for first layer
-         DO j=jj_begin-offset,jj_end+offset
-            DO i=ii_begin-offset,ii_end+offset
-               ij=(j-1)*iim+i
-               pks(ij) = cpp * ( ps(ij)/preff ) ** kappa
-               pk(ij,1) = ( p(ij,1)*pks(ij) - 0.5*alpha(ij,2)*p(ij,2) )     /    &
-                    (  p(ij,1)* (1.+kappa) + 0.5*( beta(ij,2)-(1.+2*kappa) )* p(ij,2) )
-            ENDDO
-         ENDDO
-       
-       ! for other layers
-       
-         DO l = 2, llm
-            DO j=jj_begin-offset,jj_end+offset
-               DO i=ii_begin-offset,ii_end+offset
-                  ij=(j-1)*iim+i
-                  pk(ij,l) = alpha(ij,l) + beta(ij,l) * pk(ij,l-1)
-               ENDDO
-            ENDDO
-         ENDDO
-
-       ENDIF
-
-    ELSE ! Simple calculation of Exner pressure based on centered average
-       ! surface : pks
-       IF (is_omp_level_master) THEN
-
-         DO j=jj_begin-offset,jj_end+offset
-            DO i=ii_begin-offset,ii_end+offset
-               ij=(j-1)*iim+i
-               pks(ij) = cpp * ( ps(ij)/preff ) ** kappa
-            ENDDO
-         ENDDO
-
-       ENDIF
-
-         ! 3D : pk
-       DO l = 1, llm
-          DO j=jj_begin-offset,jj_end+offset
-             DO i=ii_begin-offset,ii_end+offset
-               ij=(j-1)*iim+i
-               pk(ij,l) = cpp * ((.5/preff)*(p(ij,l)+p(ij,l+1))) ** kappa
+    ENDIF
+    
+    ! 3D : pk
+    DO l = 1, llm
+       DO j=jj_begin-offset,jj_end+offset
+          DO i=ii_begin-offset,ii_end+offset
+             ij=(j-1)*iim+i
+             pk(ij,l) = cpp * ((.5/preff)*(p(ij,l)+p(ij,l+1))) ** kappa
              ENDDO
           ENDDO
        ENDDO
-
-    END IF
-    
+       
   END SUBROUTINE compute_exner
   
Index: codes/icosagcm/trunk/src/geopotential_mod.f90
===================================================================
--- codes/icosagcm/trunk/src/geopotential_mod.f90	(revision 427)
+++ codes/icosagcm/trunk/src/geopotential_mod.f90	(revision 428)
@@ -2,48 +2,49 @@
   IMPLICIT NONE
   PRIVATE
-
+  
   PUBLIC :: compute_geopotential
 CONTAINS
+  
+  SUBROUTINE thetarhodz2geopot(f_ps,f_phis,f_theta_rhodz, f_p,f_theta,f_phi) ! FIXME : never called, dry only
+    USE icosa
+    USE omp_para
+    USE pression_mod
+    USE theta2theta_rhodz_mod
+    TYPE(t_field), POINTER :: f_ps(:), f_phis(:), f_theta_rhodz(:), &  ! IN                                       
+         f_p(:), f_theta(:), f_phi(:)               ! OUT                                                         
+    REAL(rstd),POINTER :: pk(:,:), p(:,:), theta(:,:,:), theta_rhodz(:,:,:), &
+         phi(:,:), phis(:), ps(:)
+    INTEGER :: ind
+    
+    DO ind=1,ndomain
+       IF (.NOT. assigned_domain(ind)) CYCLE
+       CALL swap_dimensions(ind)
+       CALL swap_geometry(ind)
+       ps = f_ps(ind)
+       p  = f_p(ind)
+!$OMP BARRIER                                                                                                     
+       CALL compute_pression(ps,p,0)
+!$OMP BARRIER                                                                                                     
+       theta_rhodz = f_theta_rhodz(ind)
+       theta = f_theta(ind)
+       CALL compute_theta_rhodz2theta(ps, theta_rhodz,theta,0)
+       phis = f_phis(ind)
+       phi = f_phi(ind)
+       CALL compute_geopotential(phis,p,theta,phi,0)
+    END DO
 
-  SUBROUTINE geopotential(f_phis,f_pks,f_pk,f_theta,f_phi) ! ORPHAN
-  USE icosa
-    TYPE(t_field), POINTER :: f_phis(:)
-    TYPE(t_field), POINTER :: f_pks(:)
-    TYPE(t_field), POINTER :: f_pk(:)
-    TYPE(t_field), POINTER :: f_theta(:)
-    TYPE(t_field), POINTER :: f_phi(:)
-  
-    REAL(rstd), POINTER :: phis(:)
-    REAL(rstd), POINTER :: pks(:)
-    REAL(rstd), POINTER :: pk(:,:)
-    REAL(rstd), POINTER :: theta(:,:)
-    REAL(rstd), POINTER :: phi(:,:)
-    INTEGER :: ind
+  END SUBROUTINE thetarhodz2geopot
 
-    DO ind=1,ndomain
-      IF (.NOT. assigned_domain(ind)) CYCLE
-      CALL swap_dimensions(ind)
-      CALL swap_geometry(ind)
-      phis=f_phis(ind)
-      pks=f_pks(ind)
-      pk=f_pk(ind)
-      theta=f_theta(ind)
-      phi=f_phi(ind)
-      CALL compute_geopotential(phis,pks,pk,theta,phi,0)
-    ENDDO
-  
-  END SUBROUTINE geopotential
-  
-  SUBROUTINE compute_geopotential(phis,pks,pk,theta,phi,offset)
-  USE icosa
-  USE omp_para
-  IMPLICIT NONE
+
+  SUBROUTINE compute_geopotential(phis,p,theta,phi,offset)
+    USE icosa
+    USE omp_para
+    REAL(rstd),INTENT(IN) :: p(iim*jjm,llm+1)
     REAL(rstd),INTENT(IN) :: phis(iim*jjm)
-    REAL(rstd),INTENT(IN) :: pks(iim*jjm)
-    REAL(rstd),INTENT(IN) :: pk(iim*jjm,llm)
-    REAL(rstd),INTENT(IN) :: theta(iim*jjm,llm)
+    REAL(rstd),INTENT(IN) :: theta(iim*jjm,llm,nqdyn)
     REAL(rstd),INTENT(OUT) :: phi(iim*jjm,llm)
     INTEGER,INTENT(IN) :: offset
     INTEGER :: i,j,ij,l
+    REAL(rstd) :: mg_ij, p_ij, exner_ij
 
 !!! Compute geopotential
@@ -57,19 +58,22 @@
 ! flush pks, pk thetha, phis
 !$OMP BARRIER
-   IF(is_omp_level_master) THEN
-     DO j=jj_begin-offset,jj_end+offset
-       DO i=ii_begin-offset,ii_end+offset
-         ij=(j-1)*iim+i
-         phi( ij,1 ) = phis( ij ) + theta(ij,1) * ( pks(ij) - pk(ij,1) )
-       ENDDO
-     ENDDO
+    IF(is_omp_level_master) THEN
+       DO j=jj_begin-offset,jj_end+offset
+          DO i=ii_begin-offset,ii_end+offset
+             ij=(j-1)*iim+i
+             phi( ij,1 ) = phis( ij )
+          END DO
+       END DO
           
     ! for other layers
-     DO l = 2, llm
+     DO l = 1, llm-1
        DO j=jj_begin-offset,jj_end+offset
          DO i=ii_begin-offset,ii_end+offset
            ij=(j-1)*iim+i
-           phi(ij,l) = phi(ij,l-1) + 0.5 * ( theta(ij,l)  + theta(ij,l-1) )  & 
-                                         * (  pk(ij,l-1) -  pk(ij,l)    )
+           mg_ij = p(ij,l+1)-p(ij,l)
+           ! v=RT/p=(kappa.cpp).Theta.(p/preff)**kappa /p
+           p_ij=.5*(p(ij,l+1)+p(ij,l))
+           exner_ij = (p_ij/preff)**kappa ! NB : no cpp factor
+           phi(ij,l+1) = phi(ij,l) + (kappa*cpp)*mg_ij*theta(ij,l,1)*exner_ij/p_ij
          ENDDO
        ENDDO
Index: codes/icosagcm/trunk/src/observable.f90
===================================================================
--- codes/icosagcm/trunk/src/observable.f90	(revision 427)
+++ codes/icosagcm/trunk/src/observable.f90	(revision 428)
@@ -158,79 +158,5 @@
   END SUBROUTINE write_output_fields_basic
   
-  SUBROUTINE write_output_fields(f_ps, f_phis, f_dps, f_u, f_theta_rhodz, f_q, &
-       f_buf_i, f_buf_v, f_buf_i3, f_buf1_i, f_buf2_i, f_buf_s, f_buf_p)
-    USE vorticity_mod
-    USE theta2theta_rhodz_mod
-    USE pression_mod
-    USE omega_mod
-    USE write_field_mod
-    USE vertical_interp_mod
-    USE wind_mod
-    TYPE(t_field),POINTER :: f_ps(:), f_phis(:), f_u(:), f_theta_rhodz(:), f_q(:), f_dps(:), &
-         f_buf_i(:), f_buf_v(:), f_buf_i3(:), f_buf1_i(:), f_buf2_i(:), f_buf_s(:), f_buf_p(:)
-    
-    REAL(rstd) :: out_pression_level
-    CHARACTER(LEN=255) :: str_pression
-    CHARACTER(LEN=255) :: physics_type
-    
-    out_pression_level=0.
-    CALL getin("out_pression_level",out_pression_level) 
-    WRITE(str_pression,*) INT(out_pression_level/100)
-    str_pression=ADJUSTL(str_pression)
-    
-    CALL writefield("ps",f_ps)
-    CALL writefield("dps",f_dps)
-    CALL writefield("phis",f_phis)
-    CALL vorticity(f_u,f_buf_v)
-    CALL writefield("vort",f_buf_v)
-
-    CALL w_omega(f_ps, f_u, f_buf_i)
-    CALL writefield('omega', f_buf_i)
-    IF (out_pression_level<=preff .AND. out_pression_level > 0) THEN
-      CALL vertical_interp(f_ps,f_buf_i,f_buf_s,out_pression_level)
-      CALL writefield("omega"//TRIM(str_pression),f_buf_s)
-    ENDIF
-    
-    ! Temperature
-!    CALL theta_rhodz2temperature(f_ps,f_theta_rhodz,f_buf_i) ; ! FIXME
-     
-    CALL getin('physics',physics_type)
-    IF (TRIM(physics_type)=='dcmip') THEN
-      CALL Tv2T(f_buf_i,f_q,f_buf1_i) 
-      CALL writefield("T",f_buf1_i)
-      IF (out_pression_level<=preff .AND. out_pression_level > 0) THEN
-        CALL vertical_interp(f_ps,f_buf1_i,f_buf_s,out_pression_level)
-        CALL writefield("T"//TRIM(str_pression),f_buf_s)
-      ENDIF
-    ELSE
-      CALL writefield("T",f_buf_i)
-      IF (out_pression_level<=preff .AND. out_pression_level > 0) THEN
-        CALL vertical_interp(f_ps,f_buf_i,f_buf_s,out_pression_level)
-        CALL writefield("T"//TRIM(str_pression),f_buf_s)
-      ENDIF
-    ENDIF
-   
-    ! velocity components
-    CALL un2ulonlat(f_u, f_buf1_i, f_buf2_i)
-    CALL writefield("ulon",f_buf1_i)
-    CALL writefield("ulat",f_buf2_i)
-
-    IF (out_pression_level<=preff .AND. out_pression_level > 0) THEN
-      CALL vertical_interp(f_ps,f_buf1_i,f_buf_s,out_pression_level)
-      CALL writefield("ulon"//TRIM(str_pression),f_buf_s)
-      CALL vertical_interp(f_ps,f_buf2_i,f_buf_s,out_pression_level)
-      CALL writefield("ulat"//TRIM(str_pression),f_buf_s)
-    ENDIF
-    
-    ! geopotential ! FIXME
-    CALL thetarhodz2geopot(f_ps,f_phis,f_theta_rhodz, f_buf_s,f_buf_p,f_buf1_i,f_buf2_i,f_buf_i)
-    CALL writefield("p",f_buf_p)
-!    CALL writefield("phi",f_geopot)   ! geopotential
-    CALL writefield("theta",f_buf1_i) ! potential temperature
-    CALL writefield("pk",f_buf2_i)    ! Exner pressure
-  
-  END SUBROUTINE write_output_fields
-
-!------------------- Conversion from prognostic to observable variables ------------------
+ !------------------- Conversion from prognostic to observable variables ------------------
 
   SUBROUTINE progonostic_vel_to_horiz(f_geopot, f_ps, f_rhodz, f_u, f_W, f_uh, f_uz)
@@ -308,39 +234,4 @@
   END SUBROUTINE compute_prognostic_vel_to_horiz
 
-  SUBROUTINE thetarhodz2geopot(f_ps,f_phis,f_theta_rhodz, f_pks,f_p,f_theta,f_pk,f_phi) 
-    USE field_mod
-    USE pression_mod
-    USE exner_mod
-    USE geopotential_mod
-    USE theta2theta_rhodz_mod
-    TYPE(t_field), POINTER :: f_ps(:), f_phis(:), f_theta_rhodz(:), &  ! IN
-         f_pks(:), f_p(:), f_theta(:), f_pk(:), f_phi(:)               ! OUT
-    REAL(rstd),POINTER :: pk(:,:), p(:,:), theta(:,:), theta_rhodz(:,:,:), &
-         phi(:,:), phis(:), ps(:), pks(:)
-    INTEGER :: ind
-    
-    DO ind=1,ndomain
-       IF (.NOT. assigned_domain(ind)) CYCLE
-       CALL swap_dimensions(ind)
-       CALL swap_geometry(ind)
-       ps = f_ps(ind)
-       p  = f_p(ind)
-!$OMP BARRIER
-       CALL compute_pression(ps,p,0)
-       pk = f_pk(ind)
-       pks = f_pks(ind)
-!$OMP BARRIER
-       CALL compute_exner(ps,p,pks,pk,0)
-!$OMP BARRIER
-       theta_rhodz = f_theta_rhodz(ind)
-       theta = f_theta(ind)
-       CALL compute_theta_rhodz2theta(ps, theta_rhodz(:,:,1),theta,0)
-       phis = f_phis(ind)
-       phi = f_phi(ind)
-       CALL compute_geopotential(phis,pks,pk,theta,phi,0)
-    END DO
-
-  END SUBROUTINE thetarhodz2geopot
-  
   SUBROUTINE Tv2T(f_Tv, f_q, f_T)
     TYPE(t_field), POINTER :: f_TV(:)
Index: codes/icosagcm/trunk/src/theta_rhodz.f90
===================================================================
--- codes/icosagcm/trunk/src/theta_rhodz.f90	(revision 427)
+++ codes/icosagcm/trunk/src/theta_rhodz.f90	(revision 428)
@@ -1,10 +1,11 @@
 MODULE theta2theta_rhodz_mod
   USE field_mod
-  
+  PRIVATE
   TYPE(t_field), POINTER, SAVE  :: f_p(:)
-  TYPE(t_field), POINTER, SAVE  :: f_pks(:)
-  TYPE(t_field), POINTER, SAVE  :: f_pk(:)
-
-  PRIVATE :: f_p,f_pk,f_pks 
+
+  PUBLIC :: init_theta2theta_rhodz, theta_rhodz2theta, &
+       theta_rhodz2temperature, temperature2theta_rhodz, &
+       theta2theta_rhodz, &
+       compute_theta2theta_rhodz, compute_theta_rhodz2theta
 
 CONTAINS
@@ -14,10 +15,6 @@
   USE field_mod
   IMPLICIT NONE
-    CALL allocate_field(f_p,field_t,type_real,llm+1,name='p (theta2theta_rhodz_mod)')
-    CALL allocate_field(f_pk,field_t,type_real,llm,name='pk (theta2theta_rhodz_mod)')
-    CALL allocate_field(f_pks,field_t,type_real,name='pks (theta2theta_rhodz_mod)')
-    
+    CALL allocate_field(f_p,field_t,type_real,llm+1,name='p (theta2theta_rhodz_mod)')    
   END SUBROUTINE init_theta2theta_rhodz
-
 
 
@@ -61,6 +58,4 @@
     REAL(rstd), POINTER :: temp(:,:)
     REAL(rstd), POINTER :: p(:,:)
-    REAL(rstd), POINTER :: pk(:,:)
-    REAL(rstd), POINTER :: pks(:)
     INTEGER :: ind
 
@@ -71,6 +66,4 @@
       ps=f_ps(ind)
       p=f_p(ind)
-      pks=f_pks(ind)
-      pk=f_pk(ind)
       theta_rhodz=f_theta_rhodz(ind)
       temp=f_temp(ind)
@@ -79,7 +72,5 @@
       CALL compute_pression(ps,p,0)
 !$OMP BARRIER
-      CALL compute_exner(ps,p,pks,pk,0)
-!$OMP BARRIER
-      CALL compute_theta_rhodz2temperature(p, pk, theta_rhodz(:,:,1),temp,0)
+      CALL compute_theta_rhodz2temperature(p, theta_rhodz(:,:,1),temp,0)
     ENDDO
 !$OMP BARRIER
@@ -100,6 +91,4 @@
     REAL(rstd), POINTER :: temp(:,:)
     REAL(rstd), POINTER :: p(:,:)
-    REAL(rstd), POINTER :: pk(:,:)
-    REAL(rstd), POINTER :: pks(:)
     INTEGER :: ind
 
@@ -110,6 +99,4 @@
       ps=f_ps(ind)
       p=f_p(ind)
-      pks=f_pks(ind)
-      pk=f_pk(ind)
       theta_rhodz=f_theta_rhodz(ind)
       temp=f_temp(ind)
@@ -118,7 +105,5 @@
       CALL compute_pression(ps,p,0)
 !$OMP BARRIER
-      CALL compute_exner(ps,p,pks,pk,0)
-!$OMP BARRIER
-      CALL compute_temperature2theta_rhodz(p, pk, temp, theta_rhodz(:,:,1), 0)
+      CALL compute_temperature2theta_rhodz(p, temp, theta_rhodz(:,:,1), 0)
     ENDDO
 !$OMP BARRIER
@@ -213,5 +198,5 @@
 
 
-  SUBROUTINE compute_theta_rhodz2temperature(p,pk,theta_rhodz,temp,offset)
+  SUBROUTINE compute_theta_rhodz2temperature(p,theta_rhodz,temp,offset)
   USE icosa
   USE pression_mod
@@ -220,8 +205,8 @@
   IMPLICIT NONE
     REAL(rstd),INTENT(IN) :: p(iim*jjm,llm+1)
-    REAL(rstd),INTENT(IN) :: pk(iim*jjm,llm)
     REAL(rstd),INTENT(IN) :: theta_rhodz(iim*jjm,llm)
     REAL(rstd),INTENT(OUT) :: temp(iim*jjm,llm)
     INTEGER,INTENT(IN) :: offset
+    REAL(rstd) :: pk_ij
     INTEGER :: i,j,ij,l
         
@@ -232,5 +217,6 @@
         DO i=ii_begin-offset,ii_end+offset
           ij=(j-1)*iim+i
-          temp(ij,l) = ( theta_rhodz(ij,l) / ((p(ij,l)-p(ij,l+1))/g) ) * pk(ij,l) / cpp
+          pk_ij=((.5/preff)*(p(ij,l)+p(ij,l+1)))**kappa
+          temp(ij,l) = ( theta_rhodz(ij,l) / ((p(ij,l)-p(ij,l+1))/g) ) * pk_ij
         ENDDO
       ENDDO
@@ -241,5 +227,5 @@
   END SUBROUTINE compute_theta_rhodz2temperature
 
-  SUBROUTINE compute_temperature2theta_rhodz(p,pk,temp,theta_rhodz,offset)
+  SUBROUTINE compute_temperature2theta_rhodz(p,temp,theta_rhodz,offset)
   USE icosa
   USE pression_mod
@@ -248,8 +234,8 @@
   IMPLICIT NONE
     REAL(rstd),INTENT(IN)  :: p(iim*jjm,llm+1)
-    REAL(rstd),INTENT(IN)  :: pk(iim*jjm,llm)
     REAL(rstd),INTENT(OUT) :: theta_rhodz(iim*jjm,llm)
     REAL(rstd),INTENT(IN)  :: temp(iim*jjm,llm)
     INTEGER,INTENT(IN) :: offset
+    REAL(rstd) :: pk_ij
     INTEGER :: i,j,ij,l
 
@@ -262,5 +248,6 @@
         DO i=ii_begin-offset,ii_end+offset
           ij=(j-1)*iim+i
-          theta_rhodz(ij,l) = temp(ij,l) * ((p(ij,l)-p(ij,l+1))/g) / (pk(ij,l) / cpp )
+          pk_ij=((.5/preff)*(p(ij,l)+p(ij,l+1)))**kappa
+          theta_rhodz(ij,l) = temp(ij,l) * ((p(ij,l)-p(ij,l+1))/g) / pk_ij
         ENDDO
       ENDDO
