Changeset 658 for codes/icosagcm/devel/src/kernels_unst/caldyn_solver.k90

Timestamp:

12/30/17 02:00:38 (7 years ago)

Author:

dubos

Message:

devel/unstructured : updated kernels

File:

• codes/icosagcm/devel/src/kernels_unst/caldyn_solver.k90 (modified) (6 diffs)

codes/icosagcm/devel/src/kernels_unst/caldyn_solver.k90

@@ -1 +1 @@
    !--------------------------------------------------------------------------
    !---------------------------- caldyn_solver ----------------------------------
-   !$OMP DO SCHEDULE(STATIC)
-   DO ij = 1, primal_num
-      l=1
-      m_il(l,ij) = .5*rhodz(l,ij)
-      DO l = 2, llm
-         m_il(l,ij) = .5*(rhodz(l,ij)+rhodz(l-1,ij))
-      END DO
-      l=llm+1
-      m_il(l,ij) = .5*rhodz(l-1,ij)
-   END DO
-   !$OMP END DO
-   !
-   IF(tau>0) THEN ! solve implicit problem for geopotential
-      CALL compute_NH_geopot(tau,PHI_BOT_VAR, rhodz, m_il, theta, W, geopot)
-   END IF
    !
    ! Compute pressure (pres) and Exner function (pk)
@@ -23 +8 @@
    gamma = 1./(1.-kappa)
    vreff = Rd*Treff/preff ! reference specific volume
-   Cvd = cpp-Rd
+   Cvd = 1./(cpp-Rd)
+   Rd_preff = kappa*cpp/preff
    SELECT CASE(caldyn_thermo)
    CASE(thermo_theta)
       !$OMP DO SCHEDULE(STATIC)
       DO ij = 1, primal_num
+         !DIR$ SIMD
          DO l = 1, llm
-            rho_ij = g*rhodz(l,ij)/(geopot(l+1,ij)-geopot(l,ij))
-            X_ij = (cpp/preff)*kappa*theta(l,ij,1)*rho_ij
+            rho_ij = 1./(geopot(l+1,ij)-geopot(l,ij))
+            rho_ij = rho_ij*g*rhodz(l,ij)
+            X_ij = Rd_preff*theta(l,ij,1)*rho_ij
             ! kappa.theta.rho = p/exner
             ! => X = (p/p0)/(exner/Cp)
@@ -44 +32 @@
       !$OMP DO SCHEDULE(STATIC)
       DO ij = 1, primal_num
+         !DIR$ SIMD
          DO l = 1, llm
-            rho_ij = g*rhodz(l,ij)/(geopot(l+1,ij)-geopot(l,ij))
-            T_ij = Treff*exp( (theta(l,ij,1)+Rd*log(vreff*rho_ij))/Cvd )
+            rho_ij = 1./(geopot(l+1,ij)-geopot(l,ij))
+            rho_ij = rho_ij*g*rhodz(l,ij)
+            T_ij = Treff*exp( (theta(l,ij,1)+Rd*log(vreff*rho_ij))*Cvd )
             pres(l,ij) = rho_ij*Rd*T_ij
             pk(l,ij) = T_ij
@@ -64 +54 @@
       W(l,ij) = W(l,ij)+tau*dW(l,ij) ! update W
       dPhi(l,ij) = g*g*W(l,ij)/m_il(l,ij)
+      !DIR$ SIMD
       DO l = 2, llm
          dW(l,ij) = (1./g)*(pres(l-1,ij)-pres(l,ij)) - m_il(l,ij)
@@ -80 +71 @@
    !$OMP DO SCHEDULE(STATIC)
    DO ij = 1, primal_num
+      !DIR$ SIMD
       DO l = 1, llm
          ! compute du = -0.5*g^2.grad(w^2)
@@ -90 +82 @@
       ij_left = left(edge)
       ij_right = right(edge)
+      !DIR$ SIMD
       DO l = 1, llm
          du(l,edge) = 1.*(berni(l,ij_left)-berni(l,ij_right))