Index: /codes/icosagcm/trunk/src/caldyn_gcm.f90
===================================================================
--- /codes/icosagcm/trunk/src/caldyn_gcm.f90	(revision 166)
+++ /codes/icosagcm/trunk/src/caldyn_gcm.f90	(revision 167)
@@ -22,5 +22,5 @@
 
   PUBLIC init_caldyn, caldyn_BC, caldyn, write_output_fields, &
-       req_ps, req_mass, caldyn_eta, eta_mass, eta_lag
+       req_ps, req_mass
 
 CONTAINS
@@ -122,5 +122,5 @@
        f_hflux, f_wflux, f_dps, f_dmass, f_dtheta_rhodz, f_du)
     USE icosa
-    USE disvert_mod
+    USE disvert_mod, ONLY : caldyn_eta, eta_mass
     USE vorticity_mod
     USE kinetic_mod
@@ -285,5 +285,5 @@
 SUBROUTINE compute_pvort(ps,u,theta_rhodz, rhodz,theta,qu,qv)
   USE icosa
-  USE disvert_mod
+  USE disvert_mod, ONLY : mass_dak, mass_dbk, caldyn_eta, eta_mass
   USE exner_mod
   USE trace
@@ -318,6 +318,5 @@
            DO i=ii_begin-1,ii_end+1
               ij=(j-1)*iim+i
-              m = ( ap(l) - ap(l+1) + (bp(l)-bp(l+1))*ps(ij) )/g 
-!              m = ( mass_dak(l)+ps(ij)*mass_dbk(l) )/g
+              m = ( mass_dak(l)+ps(ij)*mass_dbk(l) )/g
               rhodz(ij,l) = m
               theta(ij,l) = theta_rhodz(ij,l)/rhodz(ij,l)
@@ -420,10 +419,9 @@
        ENDDO
 
-    ELSE
-
+    ELSE 
        ! We are using a Lagrangian vertical coordinate
        ! Pressure must be computed first top-down (temporarily stored in pk)
        ! Then Exner pressure and geopotential are computed bottom-up
-       ! Notice that the computation below should work also when caldyn_eta=eta_mass
+       ! Notice that the computation below should work also when caldyn_eta=eta_mass          
 
        ! uppermost layer
@@ -451,17 +449,29 @@
           END DO
        END DO
-       DO l = 1,llm
-          !$OMP DO SCHEDULE(STATIC) 
-          DO j=jj_begin-1,jj_end+1
-             DO i=ii_begin-1,ii_end+1
-                ij=(j-1)*iim+i
-                p_ik = pk(ij,l)
-                exner_ik = cpp * (p_ik/preff) ** kappa
-                ! 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 
-                pk(ij,l) = exner_ik
+
+       IF(boussinesq) THEN ! specific volume 1 = dphi/g/rhodz
+          DO l = 1,llm
+             !$OMP DO SCHEDULE(STATIC) 
+             DO j=jj_begin-1,jj_end+1
+                DO i=ii_begin-1,ii_end+1
+                   ij=(j-1)*iim+i
+                   geopot(ij,l+1) = geopot(ij,l) + g*rhodz(ij,l)
+                ENDDO
              ENDDO
           ENDDO
-       ENDDO
+       ELSE ! specific volume v = kappa*theta*pi/p = dphi/g/rhodz
+          DO l = 1,llm
+             !$OMP DO SCHEDULE(STATIC) 
+             DO j=jj_begin-1,jj_end+1
+                DO i=ii_begin-1,ii_end+1
+                   ij=(j-1)*iim+i
+                   p_ik = pk(ij,l)
+                   exner_ik = cpp * (p_ik/preff) ** kappa
+                   geopot(ij,l+1) = geopot(ij,l) + (g*kappa)*rhodz(ij,l)*theta(ij,l)*exner_ik/p_ik 
+                   pk(ij,l) = exner_ik
+                ENDDO
+             ENDDO
+          ENDDO
+       END IF
 
     END IF
@@ -482,5 +492,5 @@
     REAL(rstd),INTENT(IN)  :: qu(iim*3*jjm,llm)
     REAL(rstd),INTENT(IN)  :: theta(iim*jjm,llm)  ! potential temperature
-    REAL(rstd),INTENT(IN)  :: pk(iim*jjm,llm) ! Exner function
+    REAL(rstd),INTENT(INOUT) :: pk(iim*jjm,llm) ! Exner function
     REAL(rstd),INTENT(IN)  :: geopot(iim*jjm,llm+1)    ! geopotential
 
@@ -647,10 +657,31 @@
   
 !!! Compute bernouilli term = Kinetic Energy + geopotential
-   DO l=ll_begin,ll_end
-    DO j=jj_begin,jj_end
-      DO i=ii_begin,ii_end
-        ij=(j-1)*iim+i
-
-        berni(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1)) &
+    IF(boussinesq) THEN
+
+       DO l=ll_begin,ll_end
+          DO j=jj_begin,jj_end
+             DO i=ii_begin,ii_end
+                ij=(j-1)*iim+i
+                
+                berni(ij,l) = pk(ij,l) + &
+                     + 1/(4*Ai(ij))*(le(ij+u_right)*de(ij+u_right)*u(ij+u_right,l)**2 +    &
+                                     le(ij+u_rup)*de(ij+u_rup)*u(ij+u_rup,l)**2 +          &
+                                     le(ij+u_lup)*de(ij+u_lup)*u(ij+u_lup,l)**2 +          &
+                                     le(ij+u_left)*de(ij+u_left)*u(ij+u_left,l)**2 +       &
+                                     le(ij+u_ldown)*de(ij+u_ldown)*u(ij+u_ldown,l)**2 +    &
+                                     le(ij+u_rdown)*de(ij+u_rdown)*u(ij+u_rdown,l)**2 )  
+                pk(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1))
+             ENDDO
+          ENDDO
+       ENDDO
+
+    ELSE
+
+       DO l=ll_begin,ll_end
+          DO j=jj_begin,jj_end
+             DO i=ii_begin,ii_end
+                ij=(j-1)*iim+i
+                
+                berni(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1)) &
                      + 1/(4*Ai(ij))*(le(ij+u_right)*de(ij+u_right)*u(ij+u_right,l)**2 +    &
                                      le(ij+u_rup)*de(ij+u_rup)*u(ij+u_rup,l)**2 +          &
@@ -660,37 +691,37 @@
                                      le(ij+u_rdown)*de(ij+u_rdown)*u(ij+u_rdown,l)**2 )  
        
-      ENDDO
-    ENDDO
-  ENDDO
-
- 
+             ENDDO
+          ENDDO
+       ENDDO
+
+    END IF ! Boussinesq/compressible
+
 !!! Add gradients of Bernoulli and Exner functions to du
-
-   DO l=ll_begin,ll_end
-    DO j=jj_begin,jj_end
-      DO i=ii_begin,ii_end
-        ij=(j-1)*iim+i
+    DO l=ll_begin,ll_end
+       DO j=jj_begin,jj_end
+          DO i=ii_begin,ii_end
+             ij=(j-1)*iim+i
         
-        du(ij+u_right,l) = du(ij+u_right,l) + 1/de(ij+u_right) * (             &
-                          0.5*(theta(ij,l)+theta(ij+t_right,l))                &
-                             *( ne_right*pk(ij,l)+ne_left*pk(ij+t_right,l))    &
-                        + ne_right*berni(ij,l)+ne_left*berni(ij+t_right,l) )
+             du(ij+u_right,l) = du(ij+u_right,l) + 1/de(ij+u_right) * (             &
+                               0.5*(theta(ij,l)+theta(ij+t_right,l))                &
+                                  *( ne_right*pk(ij,l)+ne_left*pk(ij+t_right,l))    &
+                                   + ne_right*berni(ij,l)+ne_left*berni(ij+t_right,l) )
         
         
-        du(ij+u_lup,l) = du(ij+u_lup,l) +  1/de(ij+u_lup) * (                  &
-                          0.5*(theta(ij,l)+theta(ij+t_lup,l))                  &
-                             *( ne_lup*pk(ij,l)+ne_rdown*pk(ij+t_lup,l))       &
-                        + ne_lup*berni(ij,l)+ne_rdown*berni(ij+t_lup,l) )
+             du(ij+u_lup,l) = du(ij+u_lup,l) +  1/de(ij+u_lup) * (                  &
+                               0.5*(theta(ij,l)+theta(ij+t_lup,l))                  &
+                                  *( ne_lup*pk(ij,l)+ne_rdown*pk(ij+t_lup,l))       &
+                                   + ne_lup*berni(ij,l)+ne_rdown*berni(ij+t_lup,l) )
                 
-        du(ij+u_ldown,l) = du(ij+u_ldown,l) + 1/de(ij+u_ldown) * (             &
-                          0.5*(theta(ij,l)+theta(ij+t_ldown,l))                &
-                             *( ne_ldown*pk(ij,l)+ne_rup*pk(ij+t_ldown,l))     &
-                        + ne_ldown*berni(ij,l)+ne_rup*berni(ij+t_ldown,l) )
-
-      ENDDO
+             du(ij+u_ldown,l) = du(ij+u_ldown,l) + 1/de(ij+u_ldown) * (             &
+                               0.5*(theta(ij,l)+theta(ij+t_ldown,l))                &
+                                  *( ne_ldown*pk(ij,l)+ne_rup*pk(ij+t_ldown,l))     &
+                                   + ne_ldown*berni(ij,l)+ne_rup*berni(ij+t_ldown,l) )
+
+          ENDDO
+       ENDDO
     ENDDO
-  ENDDO
-
-  CALL trace_end("compute_caldyn_horiz")
+ 
+    CALL trace_end("compute_caldyn_horiz")
 
 END SUBROUTINE compute_caldyn_horiz
Index: /codes/icosagcm/trunk/src/dissip_gcm.f90
===================================================================
--- /codes/icosagcm/trunk/src/dissip_gcm.f90	(revision 166)
+++ /codes/icosagcm/trunk/src/dissip_gcm.f90	(revision 167)
@@ -27,5 +27,4 @@
   REAL, save    :: rayleigh_tau
   
-!  INTEGER,SAVE :: itau_dissip
   REAL,SAVE    :: dtdissip
  
@@ -368,5 +367,9 @@
     fact=2
     DO l=1,llm
-       zz= 1. - preff/presnivs(l)
+       IF(ap_bp_present) THEN ! height-dependent dissipation
+          zz= 1. - preff/presnivs(l)
+       ELSE
+          zz = 0.
+       END IF
        zvert=fact-(fact-1)/(1+zz*zz)
        tau_graddiv(l) = tau_graddiv(l)/zvert
@@ -390,5 +393,5 @@
  
   
-  SUBROUTINE dissip(f_ue,f_due,f_ps,f_phis,f_theta_rhodz,f_dtheta_rhodz)
+  SUBROUTINE dissip(f_ue,f_due,f_mass,f_phis,f_theta_rhodz,f_dtheta_rhodz)
   USE icosa
   USE theta2theta_rhodz_mod
@@ -402,5 +405,5 @@
     TYPE(t_field),POINTER :: f_ue(:)
     TYPE(t_field),POINTER :: f_due(:)
-    TYPE(t_field),POINTER :: f_ps(:), f_phis(:)
+    TYPE(t_field),POINTER :: f_mass(:), f_phis(:)
     TYPE(t_field),POINTER :: f_theta_rhodz(:)
     TYPE(t_field),POINTER :: f_dtheta_rhodz(:)
@@ -421,6 +424,6 @@
     CALL gradiv(f_ue,f_due_diss1)
     CALL gradrot(f_ue,f_due_diss2)
-    
-    CALL divgrad_theta_rhodz(f_ps,f_theta_rhodz,f_dtheta_rhodz_diss)
+
+    CALL divgrad_theta_rhodz(f_mass,f_theta_rhodz,f_dtheta_rhodz_diss)
     
 ! later for openmp    
@@ -646,15 +649,14 @@
   END SUBROUTINE divgrad
    
-  SUBROUTINE divgrad_theta_rhodz(f_ps,f_theta_rhodz,f_dtheta_rhodz)
+  SUBROUTINE divgrad_theta_rhodz(f_mass,f_theta_rhodz,f_dtheta_rhodz)
   USE icosa
   USE trace
   USE omp_para
-  USE disvert_mod
   IMPLICIT NONE
-    TYPE(t_field),POINTER :: f_ps(:)
+    TYPE(t_field),POINTER :: f_mass(:)
     TYPE(t_field),POINTER :: f_theta_rhodz(:)
     TYPE(t_field),POINTER :: f_dtheta_rhodz(:)
     
-    REAL(rstd),POINTER :: ps(:)
+    REAL(rstd),POINTER :: mass(:,:)
     REAL(rstd),POINTER :: theta_rhodz(:,:)
     REAL(rstd),POINTER :: dtheta_rhodz(:,:)
@@ -668,5 +670,5 @@
       CALL swap_dimensions(ind)
       CALL swap_geometry(ind)
-      ps=f_ps(ind)
+      mass=f_mass(ind)
       theta_rhodz=f_theta_rhodz(ind)
       dtheta_rhodz=f_dtheta_rhodz(ind) 
@@ -675,5 +677,5 @@
           DO i=ii_begin,ii_end
             ij=(j-1)*iim+i
-            dtheta_rhodz(ij,l) = theta_rhodz(ij,l) / ( ap(l)-ap(l+1) + (bp(l)-bp(l+1))*ps(ij))/g
+            dtheta_rhodz(ij,l) = theta_rhodz(ij,l) / mass(ij,l)
           ENDDO
         ENDDO
@@ -703,5 +705,5 @@
           DO i=ii_begin,ii_end
             ij=(j-1)*iim+i
-            dtheta_rhodz(ij,l) = dtheta_rhodz(ij,l) * ( ap(l)-ap(l+1) + (bp(l)-bp(l+1))*ps(ij))/g
+            dtheta_rhodz(ij,l) = dtheta_rhodz(ij,l) * mass(ij,l)
           ENDDO
         ENDDO
Index: /codes/icosagcm/trunk/src/timeloop_gcm.f90
===================================================================
--- /codes/icosagcm/trunk/src/timeloop_gcm.f90	(revision 166)
+++ /codes/icosagcm/trunk/src/timeloop_gcm.f90	(revision 167)
@@ -96,17 +96,4 @@
 !----------------------------------------------------
 
-
-    def='eta_mass'
-    CALL getin('caldyn_eta',def)
-    SELECT CASE(TRIM(def))
-    CASE('eta_mass')
-       caldyn_eta=eta_mass
-    CASE('eta_lag')
-       caldyn_eta=eta_lag
-    CASE DEFAULT
-       IF (is_mpi_root) PRINT *,'Bad selector for variable caldyn_eta : <', TRIM(def),'> options are <eta_mass>, <eta_lag>'
-       STOP
-    END SELECT
-    IF (is_mpi_root) PRINT *, 'caldyn_eta=',TRIM(def), caldyn_eta, eta_lag, eta_mass
 
 ! Time-independant orography
@@ -284,7 +271,15 @@
 
     IF (MOD(it+1,itau_dissip)==0) THEN
-      CALL dissip(f_u,f_du,f_ps,f_phis, f_theta_rhodz,f_dtheta_rhodz)
-      CALL euler_scheme(.FALSE.)
-    ENDIF
+       IF(caldyn_eta==eta_mass) THEN
+          DO ind=1,ndomain
+             CALL swap_dimensions(ind)
+             CALL swap_geometry(ind)
+             mass=f_mass(ind); ps=f_ps(ind);
+             CALL compute_rhodz(.TRUE., ps, mass)
+          END DO
+       ENDIF
+       CALL dissip(f_u,f_du,f_mass,f_phis, f_theta_rhodz,f_dtheta_rhodz)
+       CALL euler_scheme(.FALSE.)  ! update only u, theta
+    END IF
 
     IF(MOD(it+1,itau_adv)==0) THEN
