Index: /codes/icosagcm/devel/src/dynamics/caldyn.f90
===================================================================
--- /codes/icosagcm/devel/src/dynamics/caldyn.f90	(revision 561)
+++ /codes/icosagcm/devel/src/dynamics/caldyn.f90	(revision 562)
@@ -2,8 +2,9 @@
   USE icosa
   PRIVATE
+  SAVE
   CHARACTER(LEN=255),SAVE :: caldyn_type
 !$OMP THREADPRIVATE(caldyn_type)
-  
-  PUBLIC init_caldyn, caldyn, caldyn_BC
+ 
+  PUBLIC init_caldyn, caldyn, caldyn_BC, dysl
   
 CONTAINS
@@ -28,5 +29,17 @@
     END SELECT
         
-      
+    dysl=.FALSE.
+    CALL getin("dysl",dysl)
+    dysl_geopot=dysl
+    CALL getin("dysl_geopot",dysl_geopot)
+    dysl_caldyn_fast=dysl
+    CALL getin("dysl_caldyn_fast",dysl_caldyn_fast)
+    dysl_slow_hydro=dysl
+    CALL getin("dysl_slow_hydro",dysl_slow_hydro)
+    dysl_pvort_only=dysl
+    CALL getin("dysl_pvort_only",dysl_pvort_only)
+    dysl_caldyn_coriolis=dysl
+    CALL getin("dysl_caldyn_coriolis",dysl_caldyn_coriolis)
+
   END SUBROUTINE init_caldyn
 
Index: /codes/icosagcm/devel/src/dynamics/caldyn_hevi.f90
===================================================================
--- /codes/icosagcm/devel/src/dynamics/caldyn_hevi.f90	(revision 561)
+++ /codes/icosagcm/devel/src/dynamics/caldyn_hevi.f90	(revision 562)
@@ -48,5 +48,5 @@
     TYPE(t_field) :: f_dPhi_fast(:)
     
-    REAL(rstd),POINTER :: ps(:), dps(:)
+    REAL(rstd),POINTER :: ps(:), dps(:), phis(:)
     REAL(rstd),POINTER :: mass(:,:), theta_rhodz(:,:,:), dtheta_rhodz(:,:,:)
     REAL(rstd),POINTER :: du(:,:), dW(:,:), dPhi(:,:), hflux(:,:), wflux(:,:)
@@ -115,4 +115,5 @@
           CALL compute_geopot(mass,theta, ps,pk,geopot)
        ELSE
+          phis = f_phis(ind)
           W = f_W(ind)
           dW = f_dW_fast(ind)
@@ -121,5 +122,5 @@
           m_il = f_wil(ind)
           pres = f_gradPhi2(ind)
-          CALL compute_caldyn_solver(tau,mass,theta,pk,geopot,W, m_il,pres, dPhi,dW,du) ! computes d(Phi,W,du)_fast and updates Phi,W
+          CALL compute_caldyn_solver(tau,phis, mass,theta,pk,geopot,W, m_il,pres, dPhi,dW,du) ! computes d(Phi,W,du)_fast and updates Phi,W
        END IF
        u=f_u(ind)
Index: /codes/icosagcm/devel/src/dynamics/caldyn_kernels_base.F90
===================================================================
--- /codes/icosagcm/devel/src/dynamics/caldyn_kernels_base.F90	(revision 561)
+++ /codes/icosagcm/devel/src/dynamics/caldyn_kernels_base.F90	(revision 562)
@@ -49,9 +49,11 @@
     Rd = kappa*cpp
 
+    IF(dysl_geopot) THEN
 #ifdef CPP_DYSL
-!#if 0
 #include "../kernels/compute_geopot.k90"
 #else
-
+    STOP 'dysl_geopot=.TRUE. but DYNAMICO has been compiled without the -dysl option'
+#endif
+    ELSE
     ! Pressure is computed first top-down (temporarily stored in pk)
     ! Then Exner pressure and geopotential are computed bottom-up
@@ -164,5 +166,5 @@
     END IF
 
-#endif
+    END IF ! dysl
 
     !ym flush geopot
Index: /codes/icosagcm/devel/src/dynamics/caldyn_kernels_hevi.F90
===================================================================
--- /codes/icosagcm/devel/src/dynamics/caldyn_kernels_hevi.F90	(revision 561)
+++ /codes/icosagcm/devel/src/dynamics/caldyn_kernels_hevi.F90	(revision 562)
@@ -9,8 +9,7 @@
   PRIVATE
 
-  REAL(rstd), PARAMETER :: pbot=1e5, Phi_bot=0., rho_bot=1e6 ! FIXME
+  REAL(rstd), PARAMETER :: pbot=1e5, rho_bot=1e6
 
   LOGICAL, SAVE :: debug_hevi_solver = .FALSE.
-  LOGICAL, PARAMETER :: rigid=.TRUE.
 
   PUBLIC :: compute_theta, compute_pvort_only, compute_caldyn_Coriolis, &
@@ -62,7 +61,12 @@
     CALL trace_start("compute_pvort_only")  
 !!! Compute shallow-water potential vorticity
+    IF(dysl_pvort_only) THEN
 #ifdef CPP_DYSL
 #include "../kernels/pvort_only.k90"
 #else
+	STOP 'dysl_pvort_only=.TRUE. but DYNAMICO is compiled without the -dysl option'
+#endif
+    ELSE
+
     radius_m2=radius**(-2)
     DO l = ll_begin,ll_end
@@ -94,11 +98,13 @@
 
     ENDDO
-#endif
+   
+    END IF ! dysl
     CALL trace_end("compute_pvort_only")
 
   END SUBROUTINE compute_pvort_only
 
-  SUBROUTINE compute_NH_geopot(tau, m_ik, m_il, theta, W_il, Phi_il)
+  SUBROUTINE compute_NH_geopot(tau, phis, m_ik, m_il, theta, W_il, Phi_il)
     REAL(rstd),INTENT(IN)    :: tau ! solve Phi-tau*dPhi/dt = Phi_rhs
+    REAL(rstd),INTENT(IN)    :: phis(iim*jjm)
     REAL(rstd),INTENT(IN)    :: m_ik(iim*jjm,llm)
     REAL(rstd),INTENT(IN)    :: m_il(iim*jjm,llm+1)
@@ -125,4 +131,5 @@
 
 #ifdef CPP_DYSL
+#define PHI_BOT(ij) phis(ij)
 #include "../kernels/compute_NH_geopot.k90"
 #else
@@ -264,6 +271,7 @@
   END SUBROUTINE compute_NH_geopot
 
-  SUBROUTINE compute_caldyn_solver(tau,rhodz,theta,pk, geopot,W, m_il,pres, dPhi,dW,du)
+  SUBROUTINE compute_caldyn_solver(tau,phis, rhodz,theta,pk, geopot,W, m_il,pres, dPhi,dW,du)
     REAL(rstd),INTENT(IN) :: tau ! "solve" Phi-tau*dPhi/dt = Phi_rhs
+    REAL(rstd),INTENT(IN)    :: phis(iim*jjm)
     REAL(rstd),INTENT(IN)    :: rhodz(iim*jjm,llm)
     REAL(rstd),INTENT(IN)    :: theta(iim*jjm,llm,nqdyn)
@@ -287,4 +295,6 @@
 #ifdef CPP_DYSL
 !$OMP BARRIER
+#define PHI_BOT(ij) phis(ij)
+#define PHI_BOT_VAR phis
 #include "../kernels/caldyn_solver.k90"
 !$OMP BARRIER
@@ -396,7 +406,12 @@
     Rd=cpp*kappa
 
+    IF(dysl_caldyn_fast) THEN
 #ifdef CPP_DYSL
 #include "../kernels/caldyn_fast.k90"
 #else
+    STOP 'dysl_caldyn_fast=.TRUE. but DYNAMICO is compiled without the -dysl option'
+#endif
+    ELSE
+
     ! Compute Bernoulli term
     IF(boussinesq) THEN
@@ -496,5 +511,6 @@
        END IF
     END DO
-#endif       
+
+    END IF ! dysl
     CALL trace_end("compute_caldyn_fast")
 
@@ -515,7 +531,12 @@
     CALL trace_start("compute_caldyn_Coriolis")
 
+    IF(dysl_caldyn_coriolis) THEN
 #ifdef CPP_DYSL
 #include "../kernels/coriolis.k90"
 #else
+    STOP 'dysl_caldyn_coriolis is .TRUE. but DYNAMICO is compiled without the -dysl option'
+#endif
+
+    ELSE
 #define FTHETA(ij) Ftheta(ij,1)
 
@@ -653,5 +674,6 @@
     END SELECT
 #undef FTHETA
-#endif
+
+    END IF ! dysl
 
     CALL trace_end("compute_caldyn_Coriolis")
@@ -672,4 +694,5 @@
     CALL trace_start("compute_caldyn_slow_hydro")
 
+    IF(dysl_slow_hydro) THEN
 #ifdef CPP_DYSL
 #define BERNI(ij,l) berni(ij,l)
@@ -677,4 +700,9 @@
 #undef BERNI
 #else
+	STOP 'dysl_slow_hydro=.TRUE. but DYNAMICO is compiled without the -dysl option'
+#endif
+
+     ELSE
+
 #define BERNI(ij) berni(ij,1)
 
@@ -726,5 +754,5 @@
     END DO
 #undef BERNI
-#endif
+    END IF ! dysl
     CALL trace_end("compute_caldyn_slow_hydro")    
   END SUBROUTINE compute_caldyn_slow_hydro
