Index: codes/icosagcm/devel/src/dynamics/compute_caldyn.f90
===================================================================
--- codes/icosagcm/devel/src/dynamics/compute_caldyn.f90	(revision 1026)
+++ codes/icosagcm/devel/src/dynamics/compute_caldyn.f90	(revision 1027)
@@ -11,60 +11,60 @@
     SUBROUTINE comp_pvort_only(u,rhodz,qu,qv, hv)
       IMPORT
-      REAL(rstd),INTENT(IN)  :: u(iim_jjm_i, llm_)
-      REAL(rstd),INTENT(INOUT) :: rhodz(iim_jjm_i, llm_)
-      REAL(rstd),INTENT(OUT) :: qu(iim_jjm_u, llm_)
-      REAL(rstd),INTENT(OUT) :: qv(iim_jjm_v, llm_)
-      REAL(rstd),INTENT(OUT) :: hv(iim_jjm_v, llm_)
+      REAL(rstd), INTENT(IN)    :: u(:,:)
+      REAL(rstd), INTENT(INOUT) :: rhodz(:,:)
+      REAL(rstd), INTENT(OUT)   :: qu(:,:)
+      REAL(rstd), INTENT(OUT)   :: qv(:,:)
+      REAL(rstd), INTENT(OUT)   :: hv(:,:)
     END SUBROUTINE comp_pvort_only
 
     SUBROUTINE comp_theta(mass_col,theta_rhodz, rhodz,theta)
       IMPORT
-      REAL(rstd),INTENT(IN)    :: mass_col(iim_jjm_i)
-      REAL(rstd),INTENT(IN)    :: theta_rhodz(iim_jjm_i, llm_, nqdyn_)
-      REAL(rstd),INTENT(INOUT) :: rhodz(iim_jjm_i, llm_)
-      REAL(rstd),INTENT(OUT)   :: theta(iim_jjm_i, llm_, nqdyn_)
+      REAL(rstd), INTENT(IN)    :: mass_col(:)
+      REAL(rstd), INTENT(IN)    :: theta_rhodz(:,:,:)
+      REAL(rstd), INTENT(INOUT) :: rhodz(:,:)
+      REAL(rstd), INTENT(OUT)   :: theta(:,:,:)
     END SUBROUTINE comp_theta
 
     SUBROUTINE comp_geopot(rhodz,theta, ps,pk,geopot) 
       IMPORT
-      REAL(rstd),INTENT(IN)    :: rhodz(iim_jjm_i, llm_)         ! rho*dz = mass per unit surface in each full model level
-      REAL(rstd),INTENT(IN)    :: theta(iim_jjm_i, llm_, nqdyn_) ! active scalars : theta/entropy, moisture, ...
-      REAL(rstd),INTENT(INOUT) :: ps(iim_jjm_i)                  ! surface pressure
-      REAL(rstd),INTENT(OUT)   :: pk(iim_jjm_i, llm_)            ! Exner function (compressible) /Lagrange multiplier (Boussinesq)
-      REAL(rstd),INTENT(INOUT) :: geopot(iim_jjm_i, llm1)        ! geopotential
+      REAL(rstd), INTENT(IN)    :: rhodz(:,:)    ! rho*dz = mass per unit surface in each full model level
+      REAL(rstd), INTENT(IN)    :: theta(:,:,:)  ! active scalars : theta/entropy, moisture, ...
+      REAL(rstd), INTENT(INOUT) :: ps(:)         ! surface pressure
+      REAL(rstd), INTENT(OUT)   :: pk(:,:)       ! Exner function (compressible) /Lagrange multiplier (Boussinesq)
+      REAL(rstd), INTENT(INOUT) :: geopot(:,:)   ! geopotential
     END SUBROUTINE comp_geopot
 
     SUBROUTINE comp_caldyn_fast(tau,theta,geopot, pk,berni,du,u)
       IMPORT
-      REAL(rstd),INTENT(IN)    :: tau                ! "solve" u-tau*du/dt = rhs
-      REAL(rstd),INTENT(IN)    :: theta(iim_jjm_i, llm_, nqdyn_)
-      REAL(rstd),INTENT(IN)    :: geopot(iim_jjm_i, llm1)
-      REAL(rstd),INTENT(INOUT) :: pk(iim_jjm_i, llm_)
-      REAL(rstd),INTENT(INOUT) :: berni(iim_jjm_i, llm_)  ! partial Bernoulli function
-      REAL(rstd),INTENT(INOUT) :: du(iim_jjm_u, llm_)
-      REAL(rstd),INTENT(INOUT) :: u(iim_jjm_u, llm_)   ! INOUT if tau>0
+      REAL(rstd), INTENT(IN)    :: tau           ! "solve" u-tau*du/dt = rhs
+      REAL(rstd), INTENT(IN)    :: theta(:,:,:)
+      REAL(rstd), INTENT(IN)    :: geopot(:,:)
+      REAL(rstd), INTENT(INOUT) :: pk(:,:)
+      REAL(rstd), INTENT(INOUT) :: berni(:,:)    ! partial Bernoulli function
+      REAL(rstd), INTENT(INOUT) :: du(:,:)
+      REAL(rstd), INTENT(INOUT) :: u(:,:)        ! INOUT if tau>0
     END SUBROUTINE comp_caldyn_fast
 
     SUBROUTINE comp_caldyn_slow_hydro(zero, u,rhodz,hv,Kv, berni, hflux,du)
       IMPORT
-      LOGICAL, INTENT(IN) :: zero
-      REAL(rstd),INTENT(IN)  :: u(iim_jjm_u, llm_)    ! prognostic "velocity"
-      REAL(rstd),INTENT(IN)  :: rhodz(iim_jjm_i, llm_)
-      REAL(rstd),INTENT(IN)  :: hv(iim_jjm_v, llm_)   ! height/mass averaged to vertices
-      REAL(rstd),INTENT(IN)  :: Kv(iim_jjm_v, llm_)   ! kinetic energy at vertices
-      REAL(rstd), INTENT(OUT) :: berni(iim_jjm_i, llm_)  ! Bernoulli function
-      REAL(rstd),INTENT(OUT) :: hflux(iim_jjm_u, llm_) ! hflux in kg/s 
-      REAL(rstd),INTENT(INOUT) :: du(iim_jjm_u, llm_)
+      LOGICAL,    INTENT(IN)    :: zero
+      REAL(rstd), INTENT(IN)    :: u(:,:)      ! prognostic "velocity"
+      REAL(rstd), INTENT(IN)    :: rhodz(:,:)
+      REAL(rstd), INTENT(IN)    :: hv(:,:)     ! height/mass averaged to vertices
+      REAL(rstd), INTENT(IN)    :: Kv(:,:)     ! kinetic energy at vertices
+      REAL(rstd), INTENT(OUT)   :: berni(:,:)  ! Bernoulli function
+      REAL(rstd), INTENT(OUT)   :: hflux(:,:)  ! hflux in kg/s 
+      REAL(rstd), INTENT(INOUT) :: du(:,:)
     END SUBROUTINE comp_caldyn_slow_hydro
     
     SUBROUTINE comp_caldyn_coriolis(hflux,theta,qu, Ftheta, convm,dtheta_rhodz,du)
       IMPORT
-      REAL(rstd),INTENT(IN)    :: hflux(iim_jjm_u, llm_)  ! hflux in kg/s
-      REAL(rstd),INTENT(IN)    :: theta(iim_jjm_i, llm_, nqdyn_) ! active scalars
-      REAL(rstd),INTENT(IN)    :: qu(iim_jjm_u, llm_)
-      REAL(rstd), INTENT(OUT)  :: Ftheta(iim_jjm_u, llm_)  ! potential temperature flux
-      REAL(rstd),INTENT(OUT)   :: convm(iim_jjm_i, llm_)  ! mass flux convergence
-      REAL(rstd),INTENT(OUT)   :: dtheta_rhodz(iim_jjm_i, llm_, nqdyn_)
-      REAL(rstd),INTENT(INOUT) :: du(iim_jjm_u, llm_)
+      REAL(rstd), INTENT(IN)    :: hflux(:,:)   ! hflux in kg/s
+      REAL(rstd), INTENT(IN)    :: theta(:,:,:) ! active scalars
+      REAL(rstd), INTENT(IN)    :: qu(:,:)
+      REAL(rstd), INTENT(OUT)   :: Ftheta(:,:)  ! potential temperature flux
+      REAL(rstd), INTENT(OUT)   :: convm(:,:)   ! mass flux convergence
+      REAL(rstd), INTENT(OUT)   :: dtheta_rhodz(:,:,:)
+      REAL(rstd), INTENT(INOUT) :: du(:,:)
     END SUBROUTINE comp_caldyn_coriolis
 
Index: codes/icosagcm/devel/src/dynamics/compute_caldyn_Coriolis.F90
===================================================================
--- codes/icosagcm/devel/src/dynamics/compute_caldyn_Coriolis.F90	(revision 1026)
+++ codes/icosagcm/devel/src/dynamics/compute_caldyn_Coriolis.F90	(revision 1027)
@@ -74,5 +74,23 @@
 #endif END_DYSL
 
-  SUBROUTINE compute_caldyn_coriolis_unst(hflux,theta,qu, Ftheta, convm,dtheta_rhodz,du)
+!-------------- Wrappers for F2008 conformity -----------------
+
+    SUBROUTINE compute_caldyn_coriolis_unst(hflux,theta,qu, Ftheta, convm,dtheta_rhodz,du)
+      REAL(rstd), INTENT(IN)    :: hflux(:,:), theta(:,:,:), qu(:,:)
+      REAL(rstd), INTENT(OUT)   :: Ftheta(:,:), convm(:,:), dtheta_rhodz(:,:,:)
+      REAL(rstd), INTENT(INOUT) :: du(:,:)
+      CALL compute_caldyn_coriolis_unst_(hflux,theta,qu, Ftheta, convm,dtheta_rhodz,du)
+    END SUBROUTINE compute_caldyn_coriolis_unst
+
+    SUBROUTINE compute_caldyn_coriolis_hex(hflux,theta,qu, Ftheta, convm,dtheta_rhodz,du)
+      REAL(rstd), INTENT(IN)    :: hflux(:,:), theta(:,:,:), qu(:,:)
+      REAL(rstd), INTENT(OUT)   :: Ftheta(:,:), convm(:,:), dtheta_rhodz(:,:,:)
+      REAL(rstd), INTENT(INOUT) :: du(:,:)
+      CALL compute_caldyn_coriolis_hex_(hflux,theta,qu, Ftheta, convm,dtheta_rhodz,du)
+    END SUBROUTINE compute_caldyn_coriolis_hex
+
+!--------------------------------------------------------------
+
+  SUBROUTINE compute_caldyn_coriolis_unst_(hflux,theta,qu, Ftheta, convm,dtheta_rhodz,du)
     USE ISO_C_BINDING, only : C_DOUBLE, C_FLOAT
     USE data_unstructured_mod, ONLY : enter_trace, exit_trace, &
@@ -89,7 +107,7 @@
 #include "../kernels_unst/coriolis.k90"
     STOP_TRACE
-  END SUBROUTINE compute_caldyn_coriolis_unst
-
-  SUBROUTINE compute_caldyn_Coriolis_hex(hflux,theta,qu, Ftheta, convm,dtheta_rhodz,du)
+  END SUBROUTINE compute_caldyn_coriolis_unst_
+
+  SUBROUTINE compute_caldyn_Coriolis_hex_(hflux,theta,qu, Ftheta, convm,dtheta_rhodz,du)
     USE icosa
     REAL(rstd),INTENT(IN)    :: hflux(3*iim*jjm,llm)  ! hflux in kg/s
@@ -107,5 +125,5 @@
 
     CALL trace_end("compute_caldyn_Coriolis")
-  END SUBROUTINE compute_caldyn_Coriolis_hex
+  END SUBROUTINE compute_caldyn_Coriolis_hex_
 
   SUBROUTINE compute_caldyn_Coriolis_manual(hflux,theta,qu, Ftheta, convm,dtheta_rhodz,du)
Index: codes/icosagcm/devel/src/dynamics/compute_caldyn_fast.F90
===================================================================
--- codes/icosagcm/devel/src/dynamics/compute_caldyn_fast.F90	(revision 1026)
+++ codes/icosagcm/devel/src/dynamics/compute_caldyn_fast.F90	(revision 1027)
@@ -68,5 +68,21 @@
 #endif END_DYSL
 
+!-------------- Wrappers for F2008 conformity -----------------
+
+  SUBROUTINE compute_caldyn_fast_hex(tau,theta,geopot, pk,berni,du,u)
+    REAL(rstd),INTENT(IN)    :: tau, theta(:,:,:), geopot(:,:)
+    REAL(rstd),INTENT(INOUT) :: pk(:,:), berni(:,:), du(:,:), u(:,:)
+    CALL compute_caldyn_fast_hex_(tau,theta,geopot, pk,berni,du,u)
+  END SUBROUTINE compute_caldyn_fast_hex
+
   SUBROUTINE compute_caldyn_fast_unst(tau,theta,geopot, pk,berni,du,u)
+    REAL(rstd),INTENT(IN)    :: tau, theta(:,:,:), geopot(:,:)
+    REAL(rstd),INTENT(INOUT) :: pk(:,:), berni(:,:), du(:,:), u(:,:)
+    CALL compute_caldyn_fast_unst_(tau,theta,geopot, pk,berni,du,u)
+  END SUBROUTINE compute_caldyn_fast_unst
+
+!--------------------------------------------------------------
+
+  SUBROUTINE compute_caldyn_fast_unst_(tau,theta,geopot, pk,berni,du,u)
     USE ISO_C_BINDING, only : C_DOUBLE, C_FLOAT
     USE data_unstructured_mod, ONLY : enter_trace, exit_trace, &
@@ -84,7 +100,7 @@
 #include "../kernels_unst/caldyn_fast.k90"
     STOP_TRACE
-  END SUBROUTINE compute_caldyn_fast_unst
-
-  SUBROUTINE compute_caldyn_fast_hex(tau,theta,geopot, pk,berni,du,u)
+  END SUBROUTINE compute_caldyn_fast_unst_
+
+  SUBROUTINE compute_caldyn_fast_hex_(tau,theta,geopot, pk,berni,du,u)
     USE icosa
     REAL(rstd),INTENT(IN)    :: tau                ! "solve" u-tau*du/dt = rhs
@@ -104,5 +120,5 @@
     CALL trace_end("compute_caldyn_fast")
     
-  END SUBROUTINE compute_caldyn_fast_hex
+  END SUBROUTINE compute_caldyn_fast_hex_
 
   SUBROUTINE compute_caldyn_fast_manual(tau,theta,geopot, pk,berni,du,u)
Index: codes/icosagcm/devel/src/dynamics/compute_caldyn_slow_hydro.F90
===================================================================
--- codes/icosagcm/devel/src/dynamics/compute_caldyn_slow_hydro.F90	(revision 1026)
+++ codes/icosagcm/devel/src/dynamics/compute_caldyn_slow_hydro.F90	(revision 1027)
@@ -53,5 +53,25 @@
 #endif END_DYSL
 
+!-------------- Wrappers for F2008 conformity -----------------
+
+  SUBROUTINE compute_caldyn_slow_hydro_hex(zero, u,rhodz,hv,Kv, berni, hflux,du)
+    LOGICAL, INTENT(IN)      :: zero
+    REAL(rstd),INTENT(IN)    :: u(:,:), rhodz(:,:), hv(:,:), Kv(:,:)
+    REAL(rstd), INTENT(OUT)  :: berni(:,:), hflux(:,:)
+    REAL(rstd),INTENT(INOUT) :: du(:,:)
+    CALL compute_caldyn_slow_hydro_hex_(zero, u,rhodz,hv,Kv, berni, hflux,du)
+  END SUBROUTINE compute_caldyn_slow_hydro_hex
+
   SUBROUTINE compute_caldyn_slow_hydro_unst(zero, u,rhodz,hv,Kv, berni, hflux,du)
+    LOGICAL, INTENT(IN)      :: zero
+    REAL(rstd),INTENT(IN)    :: u(:,:), rhodz(:,:), hv(:,:), Kv(:,:)
+    REAL(rstd), INTENT(OUT)  :: berni(:,:), hflux(:,:)
+    REAL(rstd),INTENT(INOUT) :: du(:,:)
+    CALL compute_caldyn_slow_hydro_unst_(zero, u,rhodz,hv,Kv, berni, hflux,du)
+  END SUBROUTINE compute_caldyn_slow_hydro_unst
+
+!--------------------------------------------------------------
+
+  SUBROUTINE compute_caldyn_slow_hydro_unst_(zero, u,rhodz,hv,Kv, berni, hflux,du)
     USE ISO_C_BINDING, only : C_DOUBLE, C_FLOAT
     USE data_unstructured_mod, ONLY : enter_trace, exit_trace, &
@@ -66,7 +86,7 @@
 #include "../kernels_unst/caldyn_slow_hydro.k90"
     STOP_TRACE
-  END SUBROUTINE compute_caldyn_slow_hydro_unst
-
-  SUBROUTINE compute_caldyn_slow_hydro_hex(zero, u,rhodz,hv,Kv, berni, hflux,du)
+  END SUBROUTINE compute_caldyn_slow_hydro_unst_
+
+  SUBROUTINE compute_caldyn_slow_hydro_hex_(zero, u,rhodz,hv,Kv, berni, hflux,du)
     USE icosa
     USE caldyn_vars_mod
@@ -89,5 +109,5 @@
     CALL trace_end("compute_caldyn_slow_hydro")    
 
-  END SUBROUTINE compute_caldyn_slow_hydro_hex
+  END SUBROUTINE compute_caldyn_slow_hydro_hex_
 
   SUBROUTINE compute_caldyn_slow_hydro_manual(zero, u,rhodz,hv,Kv, berni, hflux,du)
Index: codes/icosagcm/devel/src/dynamics/compute_geopot.F90
===================================================================
--- codes/icosagcm/devel/src/dynamics/compute_geopot.F90	(revision 1026)
+++ codes/icosagcm/devel/src/dynamics/compute_geopot.F90	(revision 1027)
@@ -126,7 +126,24 @@
 #endif END_DYSL
 
+  !-------------- Wrappers for F2008 conformity -----------------
+  !--------------------------------------------------------------
+
+  SUBROUTINE compute_geopot_hex(rhodz,theta, ps,pk,geopot) 
+    REAL(rstd),INTENT(IN)    :: rhodz(:,:), theta(:,:,:) ! active scalars : theta/entropy, moisture, ...
+    REAL(rstd),INTENT(INOUT) :: ps(:), geopot(:,:)       ! geopotential
+    REAL(rstd),INTENT(OUT)   :: pk(:,:)       ! Exner function (compressible) /Lagrange multiplier (Boussinesq)
+    CALL compute_geopot_hex_(rhodz,theta, ps,pk,geopot)
+  END SUBROUTINE compute_geopot_hex
+
+  SUBROUTINE compute_geopot_unst(rhodz,theta, ps,pk,geopot) 
+    REAL(rstd),INTENT(IN)    :: rhodz(:,:), theta(:,:,:) ! active scalars : theta/entropy, moisture, ...
+    REAL(rstd),INTENT(INOUT) :: ps(:), geopot(:,:)       ! geopotential
+    REAL(rstd),INTENT(OUT)   :: pk(:,:)       ! Exner function (compressible) /Lagrange multiplier (Boussinesq)
+    CALL compute_geopot_unst_(rhodz,theta, ps,pk,geopot)
+  END SUBROUTINE compute_geopot_unst
+
   !**************************** Geopotential *****************************
 
-  SUBROUTINE compute_geopot_unst(rhodz,theta,ps,pk,geopot)
+  SUBROUTINE compute_geopot_unst_(rhodz,theta,ps,pk,geopot)
         USE ISO_C_BINDING, only : C_DOUBLE, C_FLOAT
     USE data_unstructured_mod, ONLY : enter_trace, exit_trace, &
@@ -142,7 +159,7 @@
 #include "../kernels_unst/compute_geopot.k90"
     STOP_TRACE
-  END SUBROUTINE compute_geopot_unst
+  END SUBROUTINE compute_geopot_unst_
   
-  SUBROUTINE compute_geopot_hex(rhodz,theta, ps,pk,geopot) 
+  SUBROUTINE compute_geopot_hex_(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, ...
@@ -162,5 +179,5 @@
 
     CALL trace_end("compute_geopot")
-  END SUBROUTINE compute_geopot_hex
+  END SUBROUTINE compute_geopot_hex_
 
   SUBROUTINE compute_geopot_manual(rhodz,theta, ps,pk,geopot) 
Index: codes/icosagcm/devel/src/dynamics/compute_pvort_only.F90
===================================================================
--- codes/icosagcm/devel/src/dynamics/compute_pvort_only.F90	(revision 1026)
+++ codes/icosagcm/devel/src/dynamics/compute_pvort_only.F90	(revision 1027)
@@ -1,4 +1,5 @@
 MODULE compute_pvort_only_mod
   USE grid_param
+  USE prec, ONLY : rstd
   IMPLICIT NONE
   PRIVATE
@@ -36,5 +37,23 @@
 #endif END_DYSL
 
-  SUBROUTINE compute_pvort_only_unst(u,rhodz,qu,qv, hv_)
+!-------------- Wrappers for F2008 conformity -----------------
+
+  SUBROUTINE compute_pvort_only_unst(u,rhodz,qu,qv,hv_)
+    REAL(rstd),INTENT(IN)  :: u(:,:)
+    REAL(rstd),INTENT(INOUT) :: rhodz(:,:)
+    REAL(rstd),INTENT(OUT) :: qu(:,:), qv(:,:), hv_(:,:)
+    CALL compute_pvort_only_unst_(u,rhodz,qu,qv,hv_)
+  END SUBROUTINE compute_pvort_only_unst
+
+  SUBROUTINE compute_pvort_only_hex(u,rhodz,qu,qv,hv_)
+    REAL(rstd),INTENT(IN)  :: u(:,:)
+    REAL(rstd),INTENT(INOUT) :: rhodz(:,:)
+    REAL(rstd),INTENT(OUT) :: qu(:,:), qv(:,:), hv_(:,:)
+    CALL compute_pvort_only_hex_(u,rhodz,qu,qv,hv_)
+  END SUBROUTINE compute_pvort_only_hex
+
+!--------------------------------------------------------------
+
+  SUBROUTINE compute_pvort_only_unst_(u,rhodz,qu,qv, hv_)
     USE ISO_C_BINDING, only : C_DOUBLE, C_FLOAT
     USE geometry, ONLY : Riv2, Av, fv
@@ -51,7 +70,7 @@
 #include "../kernels_unst/pvort_only.k90"
     STOP_TRACE
-  END SUBROUTINE compute_pvort_only_unst
+  END SUBROUTINE compute_pvort_only_unst_
 
-  SUBROUTINE compute_pvort_only_hex(u,rhodz,qu,qv,hv_)
+  SUBROUTINE compute_pvort_only_hex_(u,rhodz,qu,qv,hv_)
     USE icosa
     USE trace, ONLY : trace_start, trace_end
@@ -108,5 +127,5 @@
     CALL trace_end("compute_pvort_only")
 
-  END SUBROUTINE compute_pvort_only_hex
+  END SUBROUTINE compute_pvort_only_hex_
 
 END MODULE compute_pvort_only_mod
Index: codes/icosagcm/devel/src/dynamics/compute_theta.F90
===================================================================
--- codes/icosagcm/devel/src/dynamics/compute_theta.F90	(revision 1026)
+++ codes/icosagcm/devel/src/dynamics/compute_theta.F90	(revision 1027)
@@ -1,3 +1,4 @@
 MODULE compute_theta_mod
+  USE prec, ONLY : rstd
   USE grid_param
   USE disvert_mod, ONLY : mass_dak, mass_dbk, caldyn_eta, eta_mass, ptop
@@ -71,5 +72,23 @@
 #endif END_DYSL
 
-  SUBROUTINE compute_theta_unst(mass_col,theta_rhodz, rhodz,theta)
+!-------------- Wrappers for F2008 conformity -----------------
+
+  SUBROUTINE compute_theta_unst(ps,theta_rhodz, rhodz,theta)
+    REAL(rstd),INTENT(IN)    :: ps(:), theta_rhodz(:,:,:)
+    REAL(rstd),INTENT(INOUT) :: rhodz(:,:)
+    REAL(rstd),INTENT(OUT)   :: theta(:,:,:)
+    CALL compute_theta_unst_(ps,theta_rhodz, rhodz,theta)
+  END SUBROUTINE compute_theta_unst
+
+  SUBROUTINE compute_theta_hex(ps,theta_rhodz, rhodz,theta)
+    REAL(rstd),INTENT(IN)    :: ps(:), theta_rhodz(:,:,:)
+    REAL(rstd),INTENT(INOUT) :: rhodz(:,:)
+    REAL(rstd),INTENT(OUT)   :: theta(:,:,:)
+    CALL compute_theta_hex_(ps,theta_rhodz, rhodz,theta)
+  END SUBROUTINE compute_theta_hex
+
+!--------------------------------------------------------------
+
+  SUBROUTINE compute_theta_unst_(mass_col,theta_rhodz, rhodz,theta)
     USE ISO_C_BINDING, only : C_DOUBLE, C_FLOAT
     USE data_unstructured_mod, ONLY : enter_trace, exit_trace, id_theta
@@ -86,7 +105,7 @@
 #undef MASS_DBK
     STOP_TRACE
-  END SUBROUTINE compute_theta_unst
+  END SUBROUTINE compute_theta_unst_
 
-  SUBROUTINE compute_theta_hex(mass_col,theta_rhodz, rhodz,theta)
+  SUBROUTINE compute_theta_hex_(mass_col,theta_rhodz, rhodz,theta)
     USE icosa
     USE trace, ONLY : trace_start, trace_end
@@ -105,5 +124,5 @@
 #undef MASS_DBK
     CALL trace_end("compute_theta")
-  END SUBROUTINE compute_theta_hex
+  END SUBROUTINE compute_theta_hex_
   
   SUBROUTINE compute_theta_manual(ps,theta_rhodz, rhodz,theta)
