Index: codes/icosagcm/devel/src/diagnostics/compute_temperature.F90
===================================================================
--- codes/icosagcm/devel/src/diagnostics/compute_temperature.F90	(revision 914)
+++ codes/icosagcm/devel/src/diagnostics/compute_temperature.F90	(revision 915)
@@ -1,4 +1,7 @@
 MODULE compute_temperature_mod
-  USE icosa
+  USE earth_const, ONLY : cpp, cppv, kappa, Rd, Rv, preff, Treff, &
+       caldyn_thermo, physics_thermo, &
+       thermo_theta, thermo_entropy, thermo_moist, thermo_fake_moist
+  USE grid_param
   IMPLICIT NONE
   PRIVATE
@@ -7,7 +10,10 @@
   PUBLIC :: temperature, compute_temperature
 
+#include "../unstructured/unstructured.h90"
+
 CONTAINS
   
   SUBROUTINE temperature(f_pmid,f_q,f_temp)
+    USE icosa
     TYPE(t_field), POINTER :: f_pmid(:)         ! IN
     TYPE(t_field), POINTER :: f_q(:)            ! IN
@@ -30,5 +36,70 @@
   END SUBROUTINE temperature
   
-  SUBROUTINE compute_temperature(pmid,q,temp)
+! Macros loop_compute_temperature_fake_moist() and loop_compute_temperature() 
+! are here to inline the thermodynamical formula inside the innermost loop.
+! Tests are made outside the outer loop.
+
+#ifdef BEGIN_DYSL
+
+{%- macro loop_compute_temperature() %}
+{%- set comp_temp=caller() %}
+  FORALL_CELLS()
+    ON_PRIMAL
+      p_ik = pmid(CELL)
+      theta_ik = temp(CELL)
+      qv = q(CELL,1) ! water vapor mixing ratio = mv/md
+      {{ comp_temp }}
+      temp(CELL) = temp_ik
+    END_BLOCK
+  END_BLOCK
+{%- endmacro %}
+
+{%- macro loop_compute_temperature_fake_moist() %}
+{%- set comp_temp=caller() %}
+IF(physics_thermo==thermo_fake_moist) THEN
+   {% call loop_compute_temperature() %}
+     {{ comp_temp }}
+     temp_ik = temp_ik/(1+0.608*qv)
+   {% endcall %}
+ELSE
+   {% call loop_compute_temperature() %}
+     {{ comp_temp }}
+   {% endcall %}
+END IF
+{%- endmacro %}
+
+KERNEL(compute_temperature)
+  SELECT CASE(caldyn_thermo)
+  CASE(thermo_theta)
+     {% call loop_compute_temperature_fake_moist() %}
+        temp_ik = theta_ik*((p_ik/preff)**kappa)
+     {% endcall %}
+  CASE(thermo_entropy)
+     {% call loop_compute_temperature_fake_moist() %}
+        temp_ik = Treff*exp((theta_ik + Rd*log(p_ik/preff))/cpp)
+     {% endcall %}
+  CASE(thermo_moist)
+     {% call loop_compute_temperature() %}
+       Rmix = Rd+qv*Rv
+       chi = ( theta_ik + Rmix*log(p_ik/preff) ) / (cpp + qv*cppv) ! log(T/Treff)
+       temp_ik = Treff*exp(chi)
+     {% endcall %}
+  END SELECT
+END_BLOCK
+
+#endif END_DYSL
+
+  SUBROUTINE compute_temperature_unst(pmid, q, temp)
+    USE prec
+    REAL(rstd),INTENT(IN)    :: pmid(llm, primal_num)
+    REAL(rstd),INTENT(IN)    :: q(llm, primal_num, nqtot)
+    REAL(rstd),INTENT(INOUT) :: temp(llm, primal_num)
+    REAL(rstd) :: p_ik, theta_ik, temp_ik, qv, chi, Rmix  
+    DECLARE_INDICES
+#include "../kernels_unst/compute_temperature.k90"
+  END SUBROUTINE compute_temperature_unst
+
+  SUBROUTINE compute_temperature_hex(pmid,q,temp)
+    USE icosa
     USE omp_para
     REAL(rstd),INTENT(IN)    :: pmid(iim*jjm,llm)
@@ -36,5 +107,17 @@
     REAL(rstd),INTENT(INOUT) :: temp(iim*jjm,llm)
 
-    REAL(rstd) :: Rd, p_ik, theta_ik, temp_ik, qv, chi, Rmix
+    REAL(rstd) :: p_ik, theta_ik, temp_ik, qv, chi, Rmix
+    INTEGER :: ij,l
+#include "../kernels_hex/compute_temperature.k90"
+  END SUBROUTINE compute_temperature_hex
+
+  SUBROUTINE compute_temperature(pmid,q,temp)
+    USE icosa
+    USE omp_para
+    REAL(rstd),INTENT(IN)    :: pmid(iim*jjm,llm)
+    REAL(rstd),INTENT(IN)    :: q(iim*jjm,llm,nqtot)
+    REAL(rstd),INTENT(INOUT) :: temp(iim*jjm,llm)
+
+    REAL(rstd) :: p_ik, theta_ik, temp_ik, qv, chi, Rmix
     INTEGER :: ij,l
 
@@ -62,4 +145,5 @@
 
   SUBROUTINE Tv2T(f_Tv, f_q, f_T)
+    USE icosa
     TYPE(t_field), POINTER :: f_TV(:)
     TYPE(t_field), POINTER :: f_q(:)
Index: codes/icosagcm/devel/src/kernels_hex/compute_temperature.k90
===================================================================
--- codes/icosagcm/devel/src/kernels_hex/compute_temperature.k90	(revision 915)
+++ codes/icosagcm/devel/src/kernels_hex/compute_temperature.k90	(revision 915)
@@ -0,0 +1,69 @@
+   !--------------------------------------------------------------------------
+   !---------------------------- compute_temperature ----------------------------------
+   SELECT CASE(caldyn_thermo)
+   CASE(thermo_theta)
+      IF(physics_thermo==thermo_fake_moist) THEN
+         DO l = ll_begin, ll_end
+            !DIR$ SIMD
+            DO ij=ij_begin, ij_end
+               p_ik = pmid(ij,l)
+               theta_ik = temp(ij,l)
+               qv = q(ij,l,1) ! water vapor mixing ratio = mv/md
+               temp_ik = theta_ik*((p_ik/preff)**kappa)
+               temp_ik = temp_ik/(1+0.608*qv)
+               temp(ij,l) = temp_ik
+            END DO
+         END DO
+      ELSE
+         DO l = ll_begin, ll_end
+            !DIR$ SIMD
+            DO ij=ij_begin, ij_end
+               p_ik = pmid(ij,l)
+               theta_ik = temp(ij,l)
+               qv = q(ij,l,1) ! water vapor mixing ratio = mv/md
+               temp_ik = theta_ik*((p_ik/preff)**kappa)
+               temp(ij,l) = temp_ik
+            END DO
+         END DO
+      END IF
+   CASE(thermo_entropy)
+      IF(physics_thermo==thermo_fake_moist) THEN
+         DO l = ll_begin, ll_end
+            !DIR$ SIMD
+            DO ij=ij_begin, ij_end
+               p_ik = pmid(ij,l)
+               theta_ik = temp(ij,l)
+               qv = q(ij,l,1) ! water vapor mixing ratio = mv/md
+               temp_ik = Treff*exp((theta_ik + Rd*log(p_ik/preff))/cpp)
+               temp_ik = temp_ik/(1+0.608*qv)
+               temp(ij,l) = temp_ik
+            END DO
+         END DO
+      ELSE
+         DO l = ll_begin, ll_end
+            !DIR$ SIMD
+            DO ij=ij_begin, ij_end
+               p_ik = pmid(ij,l)
+               theta_ik = temp(ij,l)
+               qv = q(ij,l,1) ! water vapor mixing ratio = mv/md
+               temp_ik = Treff*exp((theta_ik + Rd*log(p_ik/preff))/cpp)
+               temp(ij,l) = temp_ik
+            END DO
+         END DO
+      END IF
+   CASE(thermo_moist)
+      DO l = ll_begin, ll_end
+         !DIR$ SIMD
+         DO ij=ij_begin, ij_end
+            p_ik = pmid(ij,l)
+            theta_ik = temp(ij,l)
+            qv = q(ij,l,1) ! water vapor mixing ratio = mv/md
+            Rmix = Rd+qv*Rv
+            chi = ( theta_ik + Rmix*log(p_ik/preff) ) / (cpp + qv*cppv) ! log(T/Treff)
+            temp_ik = Treff*exp(chi)
+            temp(ij,l) = temp_ik
+         END DO
+      END DO
+   END SELECT
+   !---------------------------- compute_temperature ----------------------------------
+   !--------------------------------------------------------------------------
Index: codes/icosagcm/devel/src/kernels_unst/compute_temperature.k90
===================================================================
--- codes/icosagcm/devel/src/kernels_unst/compute_temperature.k90	(revision 915)
+++ codes/icosagcm/devel/src/kernels_unst/compute_temperature.k90	(revision 915)
@@ -0,0 +1,79 @@
+   !--------------------------------------------------------------------------
+   !---------------------------- compute_temperature ----------------------------------
+   SELECT CASE(caldyn_thermo)
+   CASE(thermo_theta)
+      IF(physics_thermo==thermo_fake_moist) THEN
+         !$OMP DO SCHEDULE(STATIC)
+         DO ij = 1, primal_num
+            !DIR$ SIMD
+            DO l = 1, llm
+               p_ik = pmid(l,ij)
+               theta_ik = temp(l,ij)
+               qv = q(l,ij,1) ! water vapor mixing ratio = mv/md
+               temp_ik = theta_ik*((p_ik/preff)**kappa)
+               temp_ik = temp_ik/(1+0.608*qv)
+               temp(l,ij) = temp_ik
+            END DO
+         END DO
+         !$OMP END DO
+      ELSE
+         !$OMP DO SCHEDULE(STATIC)
+         DO ij = 1, primal_num
+            !DIR$ SIMD
+            DO l = 1, llm
+               p_ik = pmid(l,ij)
+               theta_ik = temp(l,ij)
+               qv = q(l,ij,1) ! water vapor mixing ratio = mv/md
+               temp_ik = theta_ik*((p_ik/preff)**kappa)
+               temp(l,ij) = temp_ik
+            END DO
+         END DO
+         !$OMP END DO
+      END IF
+   CASE(thermo_entropy)
+      IF(physics_thermo==thermo_fake_moist) THEN
+         !$OMP DO SCHEDULE(STATIC)
+         DO ij = 1, primal_num
+            !DIR$ SIMD
+            DO l = 1, llm
+               p_ik = pmid(l,ij)
+               theta_ik = temp(l,ij)
+               qv = q(l,ij,1) ! water vapor mixing ratio = mv/md
+               temp_ik = Treff*exp((theta_ik + Rd*log(p_ik/preff))/cpp)
+               temp_ik = temp_ik/(1+0.608*qv)
+               temp(l,ij) = temp_ik
+            END DO
+         END DO
+         !$OMP END DO
+      ELSE
+         !$OMP DO SCHEDULE(STATIC)
+         DO ij = 1, primal_num
+            !DIR$ SIMD
+            DO l = 1, llm
+               p_ik = pmid(l,ij)
+               theta_ik = temp(l,ij)
+               qv = q(l,ij,1) ! water vapor mixing ratio = mv/md
+               temp_ik = Treff*exp((theta_ik + Rd*log(p_ik/preff))/cpp)
+               temp(l,ij) = temp_ik
+            END DO
+         END DO
+         !$OMP END DO
+      END IF
+   CASE(thermo_moist)
+      !$OMP DO SCHEDULE(STATIC)
+      DO ij = 1, primal_num
+         !DIR$ SIMD
+         DO l = 1, llm
+            p_ik = pmid(l,ij)
+            theta_ik = temp(l,ij)
+            qv = q(l,ij,1) ! water vapor mixing ratio = mv/md
+            Rmix = Rd+qv*Rv
+            chi = ( theta_ik + Rmix*log(p_ik/preff) ) / (cpp + qv*cppv) ! log(T/Treff)
+            temp_ik = Treff*exp(chi)
+            temp(l,ij) = temp_ik
+         END DO
+      END DO
+      !$OMP END DO
+   END SELECT
+   !---------------------------- compute_temperature ----------------------------------
+   !--------------------------------------------------------------------------
