Changeset 915

Timestamp:

06/17/19 18:52:40 (5 years ago)

Author:

dubos

Message:

devel : DYSL for compute_temperature

Location:

codes/icosagcm/devel/src

Files:

• diagnostics/compute_temperature.F90 (modified) (5 diffs)
• kernels_hex/compute_temperature.k90 (added)
• kernels_unst/compute_temperature.k90 (added)

codes/icosagcm/devel/src/diagnostics/compute_temperature.F90

@@ -1 +1 @@
 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 +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 +36 @@
   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 +107 @@
     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 +145 @@
 
   SUBROUTINE Tv2T(f_Tv, f_q, f_T)
+    USE icosa
     TYPE(t_field), POINTER :: f_TV(:)
     TYPE(t_field), POINTER :: f_q(:)