Changeset 604 for codes/icosagcm/trunk/src/kernels

Timestamp:

10/24/17 01:32:57 (7 years ago)

Author:

dubos

Message:

trunk : backported r600-603 from devel

File:

• codes/icosagcm/trunk/src/kernels/energy_fluxes.k90 (modified) (7 diffs)

codes/icosagcm/trunk/src/kernels/energy_fluxes.k90

@@ -11 +11 @@
       END DO
    END DO
+   ! NB : at this point pressure is stored in array pk
+   ! pk then serves as buffer to store temperature
    SELECT CASE(caldyn_thermo)
    CASE(thermo_theta)
@@ -17 +19 @@
             p_ik = pk(ij,l)
             theta_ik = theta_rhodz(ij,l,1)/rhodz(ij,l)
+            theta(ij,l) = theta_ik
             temp_ik = theta_ik*(p_ik/preff)**kappa
             gv = (g*Rd)*temp_ik/p_ik
@@ -29 +32 @@
             theta_ik = theta_rhodz(ij,l,1)/rhodz(ij,l)
             temp_ik = Treff*exp((theta_ik + Rd*log(p_ik/preff))/cpp)
+            theta(ij,l) = Treff*exp(theta_ik/cpp)
             gv = (g*Rd)*temp_ik/p_ik ! specific volume v = Rd*T/p
             pk(ij,l) = temp_ik
@@ -37 +41 @@
    !$OMP BARRIER
    ! Now accumulate energies and energy fluxes
-   ! enthalpy
    ! NB : at this point temperature is stored in array pk
    ! pk then serves as buffer to store energy
+   ! enthalpy
    DO l = ll_begin, ll_end
       !DIR$ SIMD
@@ -73 +77 @@
       END DO
    END DO
+   ! theta
+   DO l = ll_begin, ll_end
+      !DIR$ SIMD
+      DO ij=ij_begin_ext, ij_end_ext
+         energy = theta(ij,l)
+         thetat(ij,l) = thetat(ij,l) + frac*rhodz(ij,l)*energy
+         pk(ij,l) = energy
+      END DO
+   END DO
+   DO l = ll_begin, ll_end
+      !DIR$ SIMD
+      DO ij=ij_begin_ext, ij_end_ext
+         thetat_flux(ij+u_right,l) = thetat_flux(ij+u_right,l) + .5*massflux(ij+u_right,l)*(pk(ij,l)+pk(ij+t_right,l))
+         thetat_flux(ij+u_lup,l) = thetat_flux(ij+u_lup,l) + .5*massflux(ij+u_lup,l)*(pk(ij,l)+pk(ij+t_lup,l))
+         thetat_flux(ij+u_ldown,l) = thetat_flux(ij+u_ldown,l) + .5*massflux(ij+u_ldown,l)*(pk(ij,l)+pk(ij+t_ldown,l))
+      END DO
+   END DO
    ! kinetic energy
    DO l = ll_begin, ll_end
@@ -78 +99 @@
       DO ij=ij_begin_ext, ij_end_ext
          energy=0.d0
-         energy = energy + le(ij+u_rup)*de(ij+u_rup)*u(ij+u_rup,l)**2
-         energy = energy + le(ij+u_lup)*de(ij+u_lup)*u(ij+u_lup,l)**2
-         energy = energy + le(ij+u_left)*de(ij+u_left)*u(ij+u_left,l)**2
-         energy = energy + le(ij+u_ldown)*de(ij+u_ldown)*u(ij+u_ldown,l)**2
-         energy = energy + le(ij+u_rdown)*de(ij+u_rdown)*u(ij+u_rdown,l)**2
-         energy = energy + le(ij+u_right)*de(ij+u_right)*u(ij+u_right,l)**2
+         energy = energy + le(ij+u_rup)*de(ij+u_rup)*ue(ij+u_rup,l)**2
+         energy = energy + le(ij+u_lup)*de(ij+u_lup)*ue(ij+u_lup,l)**2
+         energy = energy + le(ij+u_left)*de(ij+u_left)*ue(ij+u_left,l)**2
+         energy = energy + le(ij+u_ldown)*de(ij+u_ldown)*ue(ij+u_ldown,l)**2
+         energy = energy + le(ij+u_rdown)*de(ij+u_rdown)*ue(ij+u_rdown,l)**2
+         energy = energy + le(ij+u_right)*de(ij+u_right)*ue(ij+u_right,l)**2
          energy = energy * (.25/Ai(ij))
          ekin(ij,l) = ekin(ij,l) + frac*rhodz(ij,l)*energy
@@ -97 +118 @@
       END DO
    END DO
+   ! ulon
+   DO l = ll_begin, ll_end
+      !DIR$ SIMD
+      DO ij=ij_begin_ext, ij_end_ext
+         cx=centroid(ij,1)
+         cy=centroid(ij,2)
+         cz=centroid(ij,3)
+         ux=0. ; uy=0. ; uz=0.
+         ue_le = ne_rup*ue(ij+u_rup,l)*le(ij+u_rup)
+         ux = ux + ue_le*(.5*(xyz_v(ij+z_rup,1)+xyz_v(ij+z_up,1))-cx)
+         uy = uy + ue_le*(.5*(xyz_v(ij+z_rup,2)+xyz_v(ij+z_up,2))-cy)
+         uz = uz + ue_le*(.5*(xyz_v(ij+z_rup,3)+xyz_v(ij+z_up,3))-cz)
+         ue_le = ne_lup*ue(ij+u_lup,l)*le(ij+u_lup)
+         ux = ux + ue_le*(.5*(xyz_v(ij+z_lup,1)+xyz_v(ij+z_up,1))-cx)
+         uy = uy + ue_le*(.5*(xyz_v(ij+z_lup,2)+xyz_v(ij+z_up,2))-cy)
+         uz = uz + ue_le*(.5*(xyz_v(ij+z_lup,3)+xyz_v(ij+z_up,3))-cz)
+         ue_le = ne_left*ue(ij+u_left,l)*le(ij+u_left)
+         ux = ux + ue_le*(.5*(xyz_v(ij+z_lup,1)+xyz_v(ij+z_ldown,1))-cx)
+         uy = uy + ue_le*(.5*(xyz_v(ij+z_lup,2)+xyz_v(ij+z_ldown,2))-cy)
+         uz = uz + ue_le*(.5*(xyz_v(ij+z_lup,3)+xyz_v(ij+z_ldown,3))-cz)
+         ue_le = ne_ldown*ue(ij+u_ldown,l)*le(ij+u_ldown)
+         ux = ux + ue_le*(.5*(xyz_v(ij+z_ldown,1)+xyz_v(ij+z_down,1))-cx)
+         uy = uy + ue_le*(.5*(xyz_v(ij+z_ldown,2)+xyz_v(ij+z_down,2))-cy)
+         uz = uz + ue_le*(.5*(xyz_v(ij+z_ldown,3)+xyz_v(ij+z_down,3))-cz)
+         ue_le = ne_rdown*ue(ij+u_rdown,l)*le(ij+u_rdown)
+         ux = ux + ue_le*(.5*(xyz_v(ij+z_rdown,1)+xyz_v(ij+z_down,1))-cx)
+         uy = uy + ue_le*(.5*(xyz_v(ij+z_rdown,2)+xyz_v(ij+z_down,2))-cy)
+         uz = uz + ue_le*(.5*(xyz_v(ij+z_rdown,3)+xyz_v(ij+z_down,3))-cz)
+         ue_le = ne_right*ue(ij+u_right,l)*le(ij+u_right)
+         ux = ux + ue_le*(.5*(xyz_v(ij+z_rup,1)+xyz_v(ij+z_rdown,1))-cx)
+         uy = uy + ue_le*(.5*(xyz_v(ij+z_rup,2)+xyz_v(ij+z_rdown,2))-cy)
+         uz = uz + ue_le*(.5*(xyz_v(ij+z_rup,3)+xyz_v(ij+z_rdown,3))-cz)
+         ulon_i = ux*elon_i(ij,1) + uy*elon_i(ij,2) + uz*elon_i(ij,3)
+         energy = ulon_i*(1./Ai(ij))
+         ulon(ij,l) = ulon(ij,l) + frac*rhodz(ij,l)*energy
+         pk(ij,l) = energy
+      END DO
+   END DO
+   DO l = ll_begin, ll_end
+      !DIR$ SIMD
+      DO ij=ij_begin_ext, ij_end_ext
+         ulon_flux(ij+u_right,l) = ulon_flux(ij+u_right,l) + .5*massflux(ij+u_right,l)*(pk(ij,l)+pk(ij+t_right,l))
+         ulon_flux(ij+u_lup,l) = ulon_flux(ij+u_lup,l) + .5*massflux(ij+u_lup,l)*(pk(ij,l)+pk(ij+t_lup,l))
+         ulon_flux(ij+u_ldown,l) = ulon_flux(ij+u_ldown,l) + .5*massflux(ij+u_ldown,l)*(pk(ij,l)+pk(ij+t_ldown,l))
+      END DO
+   END DO
    !---------------------------- energy_fluxes ----------------------------------
    !--------------------------------------------------------------------------