Index: codes/icosagcm/trunk/src/dynamics/caldyn_kernels_hevi.F90
===================================================================
--- codes/icosagcm/trunk/src/dynamics/caldyn_kernels_hevi.F90	(revision 561)
+++ codes/icosagcm/trunk/src/dynamics/caldyn_kernels_hevi.F90	(revision 580)
@@ -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,8 @@
     CALL trace_start("compute_pvort_only")  
 !!! Compute shallow-water potential vorticity
-#ifdef CPP_DYSL
+    IF(dysl_pvort_only) THEN
 #include "../kernels/pvort_only.k90"
-#else
+    ELSE
+
     radius_m2=radius**(-2)
     DO l = ll_begin,ll_end
@@ -94,11 +94,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)
@@ -124,7 +126,9 @@
     ij_omp_end_ext=ij_omp_end_ext+ij_begin_ext-1
 
-#ifdef CPP_DYSL
+    IF(dysl) THEN
+#define PHI_BOT(ij) phis(ij)
 #include "../kernels/compute_NH_geopot.k90"
-#else
+#undef PHI_BOT
+    ELSE
 !    FIXME : vertical OpenMP parallelism will not work
    
@@ -164,5 +168,5 @@
           B_il(ij,1) = A_ik(ij,1) + ml_g2 + tau2_g*rho_bot
           R_il(ij,1) = ml_g2*( Phi_il(ij,1)-Phi_star_il(ij,1))  &
-               + tau2_g*( p_ik(ij,1)-pbot+rho_bot*(Phi_il(ij,1)-Phi_bot) )
+               + tau2_g*( p_ik(ij,1)-pbot+rho_bot*(Phi_il(ij,1)-phis(ij)) )
        ENDDO
        ! inner interfaces
@@ -260,10 +264,11 @@
     END DO ! Newton-Raphson
 
-#endif
+    END IF ! dysl
     
   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)
@@ -278,4 +283,5 @@
 
     REAL(rstd) :: berni(iim*jjm,llm)         ! (W/m_il)^2
+    REAL(rstd) :: berni1(iim*jjm)         ! (W/m_il)^2
     REAL(rstd) :: gamma, rho_ij, T_ij, X_ij, Y_ij, vreff, Rd, Cvd
     INTEGER    :: ij, l
@@ -285,10 +291,17 @@
     Rd=cpp*kappa
 
-#ifdef CPP_DYSL
+    IF(dysl) THEN
+
 !$OMP BARRIER
+#define PHI_BOT(ij) phis(ij)
+#define PHI_BOT_VAR phis
 #include "../kernels/caldyn_solver.k90"
+#undef PHI_BOT_VAR
+#undef PHI_BOT
 !$OMP BARRIER
-#else
-#define BERNI(ij) berni(ij,1)
+
+    ELSE
+
+#define BERNI(ij) berni1(ij)
     ! FIXME : vertical OpenMP parallelism will not work
 
@@ -310,5 +323,5 @@
 
     IF(tau>0) THEN ! solve implicit problem for geopotential
-       CALL compute_NH_geopot(tau, rhodz, m_il, theta, W, geopot)
+       CALL compute_NH_geopot(tau, phis, rhodz, m_il, theta, W, geopot)
     END IF
 
@@ -372,5 +385,6 @@
     ENDDO
 #undef BERNI
-#endif
+
+    END IF ! dysl
 
     CALL trace_end("compute_caldyn_solver")
@@ -396,7 +410,8 @@
     Rd=cpp*kappa
 
-#ifdef CPP_DYSL
+    IF(dysl_caldyn_fast) THEN
 #include "../kernels/caldyn_fast.k90"
-#else
+    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,9 @@
     CALL trace_start("compute_caldyn_Coriolis")
 
-#ifdef CPP_DYSL
+    IF(dysl_caldyn_coriolis) THEN
+
 #include "../kernels/coriolis.k90"
-#else
+
+    ELSE
 #define FTHETA(ij) Ftheta(ij,1)
 
@@ -653,5 +671,6 @@
     END SELECT
 #undef FTHETA
-#endif
+
+    END IF ! dysl
 
     CALL trace_end("compute_caldyn_Coriolis")
@@ -667,4 +686,5 @@
     
     REAL(rstd) :: berni(iim*jjm,llm)  ! Bernoulli function
+    REAL(rstd) :: berni1(iim*jjm)  ! Bernoulli function
     REAL(rstd) :: uu_right, uu_lup, uu_ldown, ke, uu
     INTEGER :: ij,l
@@ -672,10 +692,13 @@
     CALL trace_start("compute_caldyn_slow_hydro")
 
-#ifdef CPP_DYSL
+    IF(dysl_slow_hydro) THEN
+
 #define BERNI(ij,l) berni(ij,l)
 #include "../kernels/caldyn_slow_hydro.k90"
 #undef BERNI
-#else
-#define BERNI(ij) berni(ij,1)
+
+     ELSE
+
+#define BERNI(ij) berni1(ij)
 
     DO l = ll_begin, ll_end
@@ -725,6 +748,8 @@
        END IF
     END DO
+
 #undef BERNI
-#endif
+
+    END IF ! dysl
     CALL trace_end("compute_caldyn_slow_hydro")    
   END SUBROUTINE compute_caldyn_slow_hydro
@@ -749,21 +774,26 @@
     REAL(rstd) :: W_el, W2_el, uu_right, uu_lup, uu_ldown, gPhi2, dP, divG, u2, uu
 
-#ifdef CPP_DYSL
     REAL(rstd) :: berni(iim*jjm,llm)  ! Bernoulli function
     REAL(rstd) :: G_el(3*iim*jjm,llm+1) ! horizontal flux of W
     REAL(rstd) :: v_el(3*iim*jjm,llm+1)
-#else
-    REAL(rstd) :: berni(iim*jjm)  ! Bernoulli function
-    REAL(rstd) :: G_el(3*iim*jjm) ! horizontal flux of W
-    REAL(rstd) :: v_el(3*iim*jjm)
-#endif
+
+    REAL(rstd) :: berni1(iim*jjm)  ! Bernoulli function
+    REAL(rstd) :: G_el1(3*iim*jjm) ! horizontal flux of W
+    REAL(rstd) :: v_el1(3*iim*jjm)
 
     CALL trace_start("compute_caldyn_slow_NH")
 
-#ifdef CPP_DYSL
+    IF(dysl) THEN
+
 !$OMP BARRIER
 #include "../kernels/caldyn_slow_NH.k90"
 !$OMP BARRIER
-#else
+     
+     ELSE
+
+#define BERNI(ij) berni1(ij)
+#define G_EL(ij) G_el1(ij)
+#define V_EL(ij) v_el1(ij)
+
     DO l=ll_begin, ll_endp1 ! compute on l levels (interfaces)
        IF(l==1) THEN
@@ -792,6 +822,6 @@
           W2_el = .5*le_de(ij+u_right) * &
                ( W(ij,l)*w_il(ij,l) + W(ij+t_right,l)*w_il(ij+t_right,l) )
-          v_el(ij+u_right) = .5*le_de(ij+u_right)*(u(ij+u_right,kup)+u(ij+u_right,kdown)) ! checked
-          G_el(ij+u_right) = v_el(ij+u_right)*W_el - DePhil(ij+u_right,l)*W2_el
+          V_EL(ij+u_right) = .5*le_de(ij+u_right)*(u(ij+u_right,kup)+u(ij+u_right,kdown)) ! checked
+          G_EL(ij+u_right) = V_EL(ij+u_right)*W_el - DePhil(ij+u_right,l)*W2_el
           ! Compute on edge 'lup'
           W_el  = .5*( W(ij,l)+W(ij+t_lup,l) )
@@ -800,6 +830,6 @@
           W2_el = .5*le_de(ij+u_lup) * &
                ( W(ij,l)*w_il(ij,l) + W(ij+t_lup,l)*w_il(ij+t_lup,l) )
-          v_el(ij+u_lup) = .5*le_de(ij+u_lup)*( u(ij+u_lup,kup) + u(ij+u_lup,kdown)) ! checked
-          G_el(ij+u_lup) = v_el(ij+u_lup)*W_el - DePhil(ij+u_lup,l)*W2_el
+          V_EL(ij+u_lup) = .5*le_de(ij+u_lup)*( u(ij+u_lup,kup) + u(ij+u_lup,kdown)) ! checked
+          G_EL(ij+u_lup) = V_EL(ij+u_lup)*W_el - DePhil(ij+u_lup,l)*W2_el
           ! Compute on edge 'ldown'
           W_el  = .5*( W(ij,l)+W(ij+t_ldown,l) )
@@ -808,6 +838,6 @@
           W2_el = .5*le_de(ij+u_ldown) * &
                ( W(ij,l)*w_il(ij,l) + W(ij+t_ldown,l)*w_il(ij+t_ldown,l) )
-          v_el(ij+u_ldown) = .5*le_de(ij+u_ldown)*( u(ij+u_ldown,kup) + u(ij+u_ldown,kdown)) ! checked
-          G_el(ij+u_ldown) = v_el(ij+u_ldown)*W_el - DePhil(ij+u_ldown,l)*W2_el
+          V_EL(ij+u_ldown) = .5*le_de(ij+u_ldown)*( u(ij+u_ldown,kup) + u(ij+u_ldown,kdown)) ! checked
+          G_EL(ij+u_ldown) = V_EL(ij+u_ldown)*W_el - DePhil(ij+u_ldown,l)*W2_el
        END DO
        ! compute GradPhi2, dPhi, dW
@@ -820,21 +850,20 @@
                le_de(ij+u_ldown)*DePhil(ij+u_ldown,l)**2 +  &
                le_de(ij+u_rdown)*DePhil(ij+u_rdown,l)**2 )
-!          gradPhi2(ij,l) = 0. ! FIXME !!
 
           dPhi(ij,l) = gradPhi2(ij,l)*w_il(ij,l) -1/(2*Ai(ij))* & 
-               ( DePhil(ij+u_right,l)*v_el(ij+u_right) + & ! -v.gradPhi, 
-                 DePhil(ij+u_rup,l)*v_el(ij+u_rup) +     & ! v_el already has le_de
-                 DePhil(ij+u_lup,l)*v_el(ij+u_lup) +     &
-                 DePhil(ij+u_left,l)*v_el(ij+u_left) +   &
-                 DePhil(ij+u_ldown,l)*v_el(ij+u_ldown) + &
-                 DePhil(ij+u_rdown,l)*v_el(ij+u_rdown) )
+               ( DePhil(ij+u_right,l)*V_EL(ij+u_right) + & ! -v.gradPhi, 
+                 DePhil(ij+u_rup,l)*V_EL(ij+u_rup) +     & ! v_el already has le_de
+                 DePhil(ij+u_lup,l)*V_EL(ij+u_lup) +     &
+                 DePhil(ij+u_left,l)*V_EL(ij+u_left) +   &
+                 DePhil(ij+u_ldown,l)*V_EL(ij+u_ldown) + &
+                 DePhil(ij+u_rdown,l)*V_EL(ij+u_rdown) )
 
           dW(ij,l) = -1./Ai(ij)*(           & ! -div(G_el), 
-               ne_right*G_el(ij+u_right) +  & ! G_el already has le_de
-               ne_rup*G_el(ij+u_rup) +      &
-               ne_lup*G_el(ij+u_lup) +      &  
-               ne_left*G_el(ij+u_left) +    &
-               ne_ldown*G_el(ij+u_ldown) +  &
-               ne_rdown*G_el(ij+u_rdown))
+               ne_right*G_EL(ij+u_right) +  & ! G_el already has le_de
+               ne_rup*G_EL(ij+u_rup) +      &
+               ne_lup*G_EL(ij+u_lup) +      &  
+               ne_left*G_EL(ij+u_left) +    &
+               ne_ldown*G_EL(ij+u_ldown) +  &
+               ne_rdown*G_EL(ij+u_rdown))
        END DO
     END DO
@@ -843,5 +872,5 @@
        ! Compute berni at scalar points
        DO ij=ij_begin_ext, ij_end_ext
-          berni(ij) = &
+          BERNI(ij) = &
                1/(4*Ai(ij))*( &
                     le_de(ij+u_right)*u(ij+u_right,l)**2 +    &
@@ -860,18 +889,23 @@
                     -0.5*(F_el(ij+u_right,l)+F_el(ij+u_right,l+1))
           hflux(ij+u_right,l) = uu_right*le_de(ij+u_right)
-          du(ij+u_right,l) = ne_right*(berni(ij)-berni(ij+t_right))
+          du(ij+u_right,l) = ne_right*(BERNI(ij)-BERNI(ij+t_right))
           ! Compute on edge 'lup'
           uu_lup = 0.5*(rhodz(ij,l)+rhodz(ij+t_lup,l))*u(ij+u_lup,l) &
                   -0.5*(F_el(ij+u_lup,l)+F_el(ij+u_lup,l+1))
           hflux(ij+u_lup,l) = uu_lup*le_de(ij+u_lup)
-          du(ij+u_lup,l) = ne_lup*(berni(ij)-berni(ij+t_lup))
+          du(ij+u_lup,l) = ne_lup*(BERNI(ij)-BERNI(ij+t_lup))
           ! Compute on edge 'ldown'
           uu_ldown = 0.5*(rhodz(ij,l)+rhodz(ij+t_ldown,l))*u(ij+u_ldown,l) &
                     -0.5*(F_el(ij+u_ldown,l)+F_el(ij+u_ldown,l+1))
           hflux(ij+u_ldown,l) = uu_ldown*le_de(ij+u_ldown)
-          du(ij+u_ldown,l) = ne_ldown*(berni(ij)-berni(ij+t_ldown))
+          du(ij+u_ldown,l) = ne_ldown*(BERNI(ij)-BERNI(ij+t_ldown))
        END DO
     END DO
-#endif
+
+#undef V_EL
+#undef G_EL
+#undef BERNI
+
+    END IF ! dysl
 
     CALL trace_end("compute_caldyn_slow_NH")
