Index: /codes/icosagcm/trunk/src/caldyn_hevi.f90
===================================================================
--- /codes/icosagcm/trunk/src/caldyn_hevi.f90	(revision 403)
+++ /codes/icosagcm/trunk/src/caldyn_hevi.f90	(revision 404)
@@ -54,5 +54,5 @@
 
 ! temporary shared variable
-    REAL(rstd),POINTER  :: theta(:,:)  
+    REAL(rstd),POINTER  :: theta(:,:,:)  
     REAL(rstd),POINTER  :: pk(:,:)
     REAL(rstd),POINTER  :: geopot(:,:)
@@ -107,5 +107,5 @@
        mass=f_mass(ind)
        theta = f_theta(ind)
-       CALL compute_theta(ps,theta_rhodz(:,:,1), mass,theta)
+       CALL compute_theta(ps,theta_rhodz, mass,theta)
        pk = f_pk(ind)
        geopot = f_geopot(ind)
@@ -113,5 +113,5 @@
        IF(hydrostatic) THEN
           du(:,:)=0.
-          CALL compute_geopot(ps,mass,theta, pk,geopot)
+          CALL compute_geopot(mass,theta, ps,pk,geopot)
        ELSE
           W = f_W(ind)
@@ -163,5 +163,5 @@
           CALL compute_caldyn_slow_NH(u,mass,geopot,W, hflux,du,dPhi,dW)
        END IF
-       CALL compute_caldyn_Coriolis(hflux,theta,qu, convm,dtheta_rhodz(:,:,1),du)
+       CALL compute_caldyn_Coriolis(hflux,theta,qu, convm,dtheta_rhodz,du)
        IF(caldyn_eta==eta_mass) THEN
           wflux=f_wflux(ind)
Index: /codes/icosagcm/trunk/src/caldyn_kernels_base.f90
===================================================================
--- /codes/icosagcm/trunk/src/caldyn_kernels_base.f90	(revision 403)
+++ /codes/icosagcm/trunk/src/caldyn_kernels_base.f90	(revision 404)
@@ -29,8 +29,8 @@
   !**************************** Geopotential *****************************
 
-  SUBROUTINE compute_geopot(ps,rhodz,theta, pk,geopot) 
+  SUBROUTINE compute_geopot(rhodz,theta, ps,pk,geopot) 
+    REAL(rstd),INTENT(IN)    :: rhodz(iim*jjm,llm)
+    REAL(rstd),INTENT(IN)    :: theta(iim*jjm,llm,nqdyn) ! active scalars : theta/entropy, moisture, ...
     REAL(rstd),INTENT(INOUT) :: ps(iim*jjm)
-    REAL(rstd),INTENT(IN)    :: rhodz(iim*jjm,llm)
-    REAL(rstd),INTENT(IN)    :: theta(iim*jjm,llm)    ! potential temperature
     REAL(rstd),INTENT(OUT)   :: pk(iim*jjm,llm)       ! Exner function (compressible) /Lagrange multiplier (Boussinesq)
     REAL(rstd),INTENT(INOUT) :: geopot(iim*jjm,llm+1) ! geopotential
@@ -48,94 +48,77 @@
     Rd = kappa*cpp
 
-    IF(caldyn_eta==eta_mass .AND. .NOT. DEC) THEN
-!!! Compute exner function and geopotential
-       STOP
+    ! Pressure is computed first top-down (temporarily stored in pk)
+    ! Then Exner pressure and geopotential are computed bottom-up
+    ! Works also when caldyn_eta=eta_mass          
+
+    IF(boussinesq) THEN ! compute geopotential and pk=Lagrange multiplier
+       ! specific volume 1 = dphi/g/rhodz
+       !         IF (is_omp_level_master) THEN ! no openMP on vertical due to dependency
        DO l = 1,llm
           !DIR$ SIMD
-          DO ij=ij_omp_begin_ext,ij_omp_end_ext
-             p_ik = ptop + mass_ak(l) + mass_bk(l)*ps(ij) ! FIXME : leave ps for the moment ; change ps to Ms later
-             exner_ik = cpp * (p_ik/preff) ** kappa
-             pk(ij,l) = exner_ik
-             ! specific volume v = kappa*theta*pi/p = dphi/g/rhodz
-             geopot(ij,l+1) = geopot(ij,l) + (g*kappa)*rhodz(ij,l)*theta(ij,l)*exner_ik/p_ik
+          DO ij=ij_omp_begin_ext,ij_omp_end_ext         
+             geopot(ij,l+1) = geopot(ij,l) + g*rhodz(ij,l)
           ENDDO
        ENDDO
-    ELSE 
-       ! We are using DEC or a Lagrangian vertical coordinate
-       ! Pressure is computed first top-down (temporarily stored in pk)
-       ! Then Exner pressure and geopotential are computed bottom-up
-       ! Works also when caldyn_eta=eta_mass          
-
-       IF(boussinesq) THEN ! compute geopotential and pk=Lagrange multiplier
-          ! specific volume 1 = dphi/g/rhodz
-          !         IF (is_omp_level_master) THEN ! no openMP on vertical due to dependency
-          DO l = 1,llm
-             !DIR$ SIMD
-             DO ij=ij_omp_begin_ext,ij_omp_end_ext         
-                geopot(ij,l+1) = geopot(ij,l) + g*rhodz(ij,l)
-             ENDDO
-          ENDDO
-          ! use hydrostatic balance with theta*rhodz to find pk (Lagrange multiplier=pressure) 
-          ! uppermost layer
+       ! use hydrostatic balance with theta*rhodz to find pk (Lagrange multiplier=pressure) 
+       ! uppermost layer
+       !DIR$ SIMD
+       DO ij=ij_begin_ext,ij_end_ext         
+          pk(ij,llm) = ptop + (.5*g)*theta(ij,llm,1)*rhodz(ij,llm)
+       END DO
+       ! other layers
+       DO l = llm-1, 1, -1
+          !          !$OMP DO SCHEDULE(STATIC) 
           !DIR$ SIMD
           DO ij=ij_begin_ext,ij_end_ext         
-             pk(ij,llm) = ptop + (.5*g)*theta(ij,llm)*rhodz(ij,llm)
+             pk(ij,l) = pk(ij,l+1) + (.5*g)*(theta(ij,l,1)*rhodz(ij,l)+theta(ij,l+1,1)*rhodz(ij,l+1))
           END DO
-          ! other layers
-          DO l = llm-1, 1, -1
-             !          !$OMP DO SCHEDULE(STATIC) 
+       END DO
+       ! now pk contains the Lagrange multiplier (pressure)
+    ELSE ! non-Boussinesq, compute pressure, Exner pressure or temperature, then geopotential
+       ! uppermost layer
+       
+       !DIR$ SIMD
+       DO ij=ij_omp_begin_ext,ij_omp_end_ext
+          pk(ij,llm) = ptop + (.5*g)*rhodz(ij,llm)
+       END DO
+       ! other layers
+       DO l = llm-1, 1, -1
+          !DIR$ SIMD
+          DO ij=ij_omp_begin_ext,ij_omp_end_ext         
+             pk(ij,l) = pk(ij,l+1) + (.5*g)*(rhodz(ij,l)+rhodz(ij,l+1))
+          END DO
+       END DO
+       ! surface pressure (for diagnostics)
+       IF(caldyn_eta==eta_lag) THEN
+          DO ij=ij_omp_begin_ext,ij_omp_end_ext         
+             ps(ij) = pk(ij,1) + (.5*g)*rhodz(ij,1)
+          END DO
+       END IF
+       
+       DO l = 1,llm
+          SELECT CASE(caldyn_thermo)
+          CASE(thermo_theta)
              !DIR$ SIMD
-             DO ij=ij_begin_ext,ij_end_ext         
-                pk(ij,l) = pk(ij,l+1) + (.5*g)*(theta(ij,l)*rhodz(ij,l)+theta(ij,l+1)*rhodz(ij,l+1))
-             END DO
-          END DO
-          ! now pk contains the Lagrange multiplier (pressure)
-       ELSE ! non-Boussinesq, compute pressure, Exner pressure or temperature, then geopotential
-          ! uppermost layer
-
-          !DIR$ SIMD
-          DO ij=ij_omp_begin_ext,ij_omp_end_ext
-             pk(ij,llm) = ptop + (.5*g)*rhodz(ij,llm)
-          END DO
-          ! other layers
-          DO l = llm-1, 1, -1
+             DO ij=ij_omp_begin_ext,ij_omp_end_ext
+                p_ik = pk(ij,l)
+                exner_ik = cpp * (p_ik/preff) ** kappa
+                pk(ij,l) = exner_ik
+                ! specific volume v = kappa*theta*pi/p = dphi/g/rhodz
+                geopot(ij,l+1) = geopot(ij,l) + (g*kappa)*rhodz(ij,l)*theta(ij,l,1)*exner_ik/p_ik
+             ENDDO
+          CASE(thermo_entropy) ! theta is in fact entropy = cpp*log(theta/Treff) = cpp*log(T/Treff) - Rd*log(p/preff)
              !DIR$ SIMD
-             DO ij=ij_omp_begin_ext,ij_omp_end_ext         
-                pk(ij,l) = pk(ij,l+1) + (.5*g)*(rhodz(ij,l)+rhodz(ij,l+1))
-             END DO
-          END DO
-          ! surface pressure (for diagnostics)
-          IF(caldyn_eta==eta_lag) THEN
-             DO ij=ij_omp_begin_ext,ij_omp_end_ext         
-                ps(ij) = pk(ij,1) + (.5*g)*rhodz(ij,1)
-             END DO
-          END IF
-
-          DO l = 1,llm
-             SELECT CASE(caldyn_thermo)
-             CASE(thermo_theta)
-                !DIR$ SIMD
-                DO ij=ij_omp_begin_ext,ij_omp_end_ext
-                   p_ik = pk(ij,l)
-                   exner_ik = cpp * (p_ik/preff) ** kappa
-                   pk(ij,l) = exner_ik
-                   ! specific volume v = kappa*theta*pi/p = dphi/g/rhodz
-                   geopot(ij,l+1) = geopot(ij,l) + (g*kappa)*rhodz(ij,l)*theta(ij,l)*exner_ik/p_ik
-                ENDDO
-             CASE(thermo_entropy) ! theta is in fact entropy = cpp*log(theta/Treff) = cpp*log(T/Treff) - Rd*log(p/preff)
-                !DIR$ SIMD
-                DO ij=ij_omp_begin_ext,ij_omp_end_ext
-                   p_ik = pk(ij,l)
-                   temp_ik = Treff*exp((theta(ij,l) + Rd*log(p_ik/preff))/cpp)
-                   pk(ij,l) = temp_ik
-                   ! specific volume v = Rd*T/p = dphi/g/rhodz
-                   geopot(ij,l+1) = geopot(ij,l) + (g*Rd)*rhodz(ij,l)*temp_ik/p_ik
-                ENDDO
-             CASE DEFAULT
-                STOP
-             END SELECT
-          ENDDO
-       END IF
-
+             DO ij=ij_omp_begin_ext,ij_omp_end_ext
+                p_ik = pk(ij,l)
+                temp_ik = Treff*exp((theta(ij,l,1) + Rd*log(p_ik/preff))/cpp)
+                pk(ij,l) = temp_ik
+                ! specific volume v = Rd*T/p = dphi/g/rhodz
+                geopot(ij,l+1) = geopot(ij,l) + (g*Rd)*rhodz(ij,l)*temp_ik/p_ik
+             ENDDO
+          CASE DEFAULT
+             STOP
+          END SELECT
+       ENDDO
     END IF
 
Index: /codes/icosagcm/trunk/src/caldyn_kernels_hevi.f90
===================================================================
--- /codes/icosagcm/trunk/src/caldyn_kernels_hevi.f90	(revision 403)
+++ /codes/icosagcm/trunk/src/caldyn_kernels_hevi.f90	(revision 404)
@@ -20,13 +20,13 @@
 
   SUBROUTINE compute_theta(ps,theta_rhodz, rhodz,theta)
-    REAL(rstd),INTENT(IN)  :: ps(iim*jjm)
-    REAL(rstd),INTENT(IN)  :: theta_rhodz(iim*jjm,llm)
+    REAL(rstd),INTENT(IN)    :: ps(iim*jjm)
+    REAL(rstd),INTENT(IN)    :: theta_rhodz(iim*jjm,llm,nqdyn)
     REAL(rstd),INTENT(INOUT) :: rhodz(iim*jjm,llm)
-    REAL(rstd),INTENT(OUT) :: theta(iim*jjm,llm)
-    INTEGER :: ij,l
+    REAL(rstd),INTENT(OUT)   :: theta(iim*jjm,llm,nqdyn)
+    INTEGER :: ij,l,iq
     REAL(rstd) :: m
-    CALL trace_start("compute_theta")  
-
-    IF(caldyn_eta==eta_mass) THEN ! Compute mass & theta
+    CALL trace_start("compute_theta")
+
+    IF(caldyn_eta==eta_mass) THEN ! Compute mass
        DO l = ll_begin,ll_end
           !DIR$ SIMD
@@ -38,15 +38,16 @@
              END IF
              rhodz(ij,l) = m/g
-             theta(ij,l) = theta_rhodz(ij,l)/rhodz(ij,l)
-          ENDDO
-       ENDDO
-    ELSE ! Compute only theta
-       DO l = ll_begin,ll_end
+          END DO
+       END DO
+    END IF
+
+    DO l = ll_begin,ll_end
+       DO iq=1,nqdyn
           !DIR$ SIMD
           DO ij=ij_begin_ext,ij_end_ext
-             theta(ij,l) = theta_rhodz(ij,l)/rhodz(ij,l)
-          ENDDO
-       ENDDO
-    END IF
+             theta(ij,l,iq) = theta_rhodz(ij,l,iq)/rhodz(ij,l)
+          END DO
+       END DO
+    END DO
 
     CALL trace_end("compute_theta")
@@ -451,13 +452,13 @@
   SUBROUTINE compute_caldyn_Coriolis(hflux,theta,qu, convm,dtheta_rhodz,du)
     REAL(rstd),INTENT(IN)  :: hflux(3*iim*jjm,llm)  ! hflux in kg/s
-    REAL(rstd),INTENT(IN)  :: theta(iim*jjm,llm)  ! potential temperature
+    REAL(rstd),INTENT(IN)  :: theta(iim*jjm,llm,nqdyn) ! active scalars
     REAL(rstd),INTENT(IN)  :: qu(3*iim*jjm,llm)
     REAL(rstd),INTENT(OUT) :: convm(iim*jjm,llm)  ! mass flux convergence
-    REAL(rstd),INTENT(OUT) :: dtheta_rhodz(iim*jjm,llm)
+    REAL(rstd),INTENT(OUT) :: dtheta_rhodz(iim*jjm,llm,nqdyn)
     REAL(rstd),INTENT(INOUT) :: du(3*iim*jjm,llm)
     
     REAL(rstd) :: Ftheta(3*iim*jjm)  ! potential temperature flux
     REAL(rstd) :: uu_right, uu_lup, uu_ldown
-    INTEGER :: ij,l,kdown
+    INTEGER :: ij,iq,l,kdown
 
     CALL trace_start("compute_caldyn_Coriolis")
@@ -466,13 +467,26 @@
        ! compute theta flux
        !DIR$ SIMD
-       DO ij=ij_begin_ext,ij_end_ext      
-          Ftheta(ij+u_right) = 0.5*(theta(ij,l)+theta(ij+t_right,l)) &
+       DO iq=1,nqdyn
+          DO ij=ij_begin_ext,ij_end_ext      
+             Ftheta(ij+u_right) = 0.5*(theta(ij,l,iq)+theta(ij+t_right,l,iq)) &
                                   * hflux(ij+u_right,l)
-          Ftheta(ij+u_lup)   = 0.5*(theta(ij,l)+theta(ij+t_lup,l)) &
-                                  * hflux(ij+u_lup,l)
-          Ftheta(ij+u_ldown) = 0.5*(theta(ij,l)+theta(ij+t_ldown,l)) &
-                                  * hflux(ij+u_ldown,l)
-       ENDDO
-       ! compute horizontal divergence of fluxes
+             Ftheta(ij+u_lup)   = 0.5*(theta(ij,l,iq)+theta(ij+t_lup,l,iq)) &
+                  * hflux(ij+u_lup,l)
+             Ftheta(ij+u_ldown) = 0.5*(theta(ij,l,iq)+theta(ij+t_ldown,l,iq)) &
+                  * hflux(ij+u_ldown,l)
+          END DO
+          ! horizontal divergence of fluxes
+          DO ij=ij_begin,ij_end         
+             ! dtheta_rhodz =  -div(flux.theta)
+             dtheta_rhodz(ij,l,iq)= &
+                  -1./Ai(ij)*(ne_right*Ftheta(ij+u_right)    +  &
+                  ne_rup*Ftheta(ij+u_rup)        +  &  
+                  ne_lup*Ftheta(ij+u_lup)        +  &  
+                  ne_left*Ftheta(ij+u_left)      +  &  
+                  ne_ldown*Ftheta(ij+u_ldown)    +  &
+                  ne_rdown*Ftheta(ij+u_rdown) )
+          END DO
+       END DO
+
        DO ij=ij_begin,ij_end         
           ! convm = -div(mass flux), sign convention as in Ringler et al. 2012, eq. 21
@@ -482,14 +496,7 @@
                ne_left*hflux(ij+u_left,l)     +  &  
                ne_ldown*hflux(ij+u_ldown,l)   +  &  
-               ne_rdown*hflux(ij+u_rdown,l))    
-          ! dtheta_rhodz =  -div(flux.theta)
-          dtheta_rhodz(ij,l)=-1./Ai(ij)*(ne_right*Ftheta(ij+u_right)    +  &
-               ne_rup*Ftheta(ij+u_rup)        +  &  
-               ne_lup*Ftheta(ij+u_lup)        +  &  
-               ne_left*Ftheta(ij+u_left)      +  &  
-               ne_ldown*Ftheta(ij+u_ldown)    +  &  
-               ne_rdown*Ftheta(ij+u_rdown))
-       ENDDO
-    END DO
+               ne_rdown*hflux(ij+u_rdown,l))
+       END DO ! ij
+    END DO ! llm
 
 !!! Compute potential vorticity (Coriolis) contribution to du
Index: /codes/icosagcm/trunk/src/earth_const.f90
===================================================================
--- /codes/icosagcm/trunk/src/earth_const.f90	(revision 403)
+++ /codes/icosagcm/trunk/src/earth_const.f90	(revision 404)
@@ -9,5 +9,5 @@
   REAL(rstd),SAVE :: kappa=0.2857143
   REAL(rstd),SAVE :: cpp=1004.70885
-  REAL(rstd),SAVE :: cppv=1004.70885 ! FIXME
+  REAL(rstd),SAVE :: cppv=1860.
   REAL(rstd),SAVE :: Rv=461.5
   REAL(rstd),SAVE :: Treff=273.
@@ -39,4 +39,6 @@
     CALL getin("kappa",kappa)  
     CALL getin("cpp",cpp)  
+    CALL getin("cppv",cppv)
+    CALL getin("Rv",Rv)
     CALL getin("preff",preff)  
     CALL getin("Treff",Treff)  
Index: /codes/icosagcm/trunk/src/observable.f90
===================================================================
--- /codes/icosagcm/trunk/src/observable.f90	(revision 403)
+++ /codes/icosagcm/trunk/src/observable.f90	(revision 404)
@@ -33,6 +33,6 @@
     CALL allocate_field(f_buf_s,   field_t,type_real, name="buf_s")
 
-    CALL allocate_field(f_theta, field_t,type_real,llm,  name='theta')   ! potential temperature
-    CALL allocate_field(f_pmid,  field_t,type_real,llm,  name='pmid')   ! mid layer pressure
+    CALL allocate_field(f_theta, field_t,type_real,llm,nqdyn,  name='theta') ! potential temperature
+    CALL allocate_field(f_pmid,  field_t,type_real,llm,  name='pmid')       ! mid layer pressure
   END SUBROUTINE init_observable
   
