Changeset 919

Timestamp:

06/18/19 22:57:38 (5 years ago)

Author:

dubos

Message:

devel : DYSL for compute_omega

Location:

codes/icosagcm/devel

Files:

• Python/src/hexagonal/macros.jin (modified) (1 diff)
• Python/src/unstructured/macros.jin (modified) (2 diffs)
• make_python (modified) (1 diff)
• src/diagnostics/compute_omega.F90 (moved) (moved from codes/icosagcm/devel/src/diagnostics/omega.f90) (6 diffs)
• src/diagnostics/observable.f90 (modified) (1 diff)

codes/icosagcm/devel/Python/src/hexagonal/macros.jin

@@ -20 +20 @@
 #define WEE wee(ij+{{edge}},{{itrisk}})
 #define DE de(ij+{{edge}})
+#define LE le(ij+{{edge}})
 #define RIV2 Riv2({{vertex}},{{vertex_dir}})
 #define AI Ai(ij)

codes/icosagcm/devel/Python/src/unstructured/macros.jin

@@ -93 +93 @@
 {{ define('EDGE', 'l,edge') }}
 {{ cdef(thecode, 'LE_DE', 'le_de(edge)') }}
+{{ cdef(thecode, 'LE', 'le(edge)') }}
+{{ cdef(thecode, 'DE', 'de(edge)') }}
 {{ cdef(thecode, 'SIGN', '1.') }}
 {{ cdef(thecode, 'CELL1', 'l,ij_left') }}
@@ -285 +287 @@
 {{ define('EDGE', 'l,edge') }}
 {{ cdef(thecode, 'LE_DE', 'le_de(edge)') }}
+{{ cdef(thecode, 'DE', 'de(edge)') }}
+{{ cdef(thecode, 'LE', 'le(edge)') }}
 {{ cdef(thecode, 'SIGN', '1.') }}
 {{ cdef(thecode, 'CELL1', 'l,ij_left') }}

codes/icosagcm/devel/make_python

@@ -57 +57 @@
 # make sure that ./rebuild recompiles dynamics
     cd $DYNAMICO_ROOT
-    touch src/dynamics/*.F90 src/diagnostics/*.F90 src/transport/*.F90 src/unstructured/*.F90
+#    touch src/dynamics/*.F90 src/diagnostics/*.F90 src/transport/*.F90 src/unstructured/*.F90
 }
 

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

@@ -1 @@
-MODULE omega_mod
-
+MODULE compute_omega_mod
   USE icosa
+  IMPLICIT NONE
   PRIVATE
 
@@ -23 +23 @@
   END SUBROUTINE W_omega
 
-
+#ifdef BEGIN_DYSL
+
+KERNEL(compute_omega)
+
+  ! Pressure
+  FORALL_CELLS_EXT()
+    ON_PRIMAL
+      p(CELL) = AP(CELL) + BP(CELL)*ps(HIDX(CELL))
+    END_BLOCK
+  END_BLOCK
+
+  BARRIER
+
+  ! Mass and grad(ps)
+  FORALL_CELLS_EXT()
+    ON_PRIMAL
+      rhodz(CELL) = (p(CELL)-p(UP(CELL)))*(1./g)
+    END_BLOCK
+    ON_EDGES
+      CST_IFTHEN(IS_BOTTOM_LEVEL)
+        gradps(HIDX(EDGE)) = (ps(HIDX(CELL2))-ps(HIDX(CELL1)))*SIGN*LE
+      CST_ENDIF
+    END_BLOCK
+  END_BLOCK
+  
+  ! Mass flux
+  FORALL_CELLS_EXT()
+    ON_EDGES
+      Fe(EDGE)=0.5*(rhodz(CELL1)+rhodz(CELL2)*LE
+    END_BLOCK
+  END_BLOCK
+
+  ! Mass flux divergence
+  ! convm = +div(mass flux), sign convention as in Ringler et al. 2012, eq. 21
+  FORALL_CELLS()
+    ON_PRIMAL
+      divflux=0.
+      FORALL_EDGES
+        divflux = divflux + SIGN*Fe(EDGE)
+      END_BLOCK
+      convm(CELL) = divflux*(1./AI)
+    END_BLOCK
+  END_BLOCK
+
+  ! Barrier needed before and after doing a vertical recurrence
+  BARRIER
+
+  ! vertical integration from up to down
+  SEQUENCE_C1
+    PROLOGUE('llm')
+      convm(CELL)=0.
+    END_BLOCK
+    BODY('llm-1,1,-1') 
+      convm(CELL) = convm(CELL) + convm(UP(CELL))
+    END_BLOCK
+  END_BLOCK
+
+  BARRIER
+
+  ! omega = dp/dt = u.grad p + \pdiff{p}{t}
+  FORALL_CELLS()
+    ON_PRIMAL
+      ugradps=0.
+      FORALL_EDGES
+        ugradps = ugradps + u(EDGE)*gradps(HIDX(EDGE))
+      END_BLOCK
+      ugradps = .5*(BP(CELL)+BP(UP(CELL)))*ugradps/(-4.*AI)  ! sign convention as in Ringler et al. 2010, Eq. 22 p.3072
+      w(CELL) =  ugradps - g*.5*(convm(CELL)+convm(UP(CELL)))
+    END_BLOCK
+  END_BLOCK
+END_BLOCK
+
+#endif END_DYSL
 
   SUBROUTINE compute_omega(ps,u, w)
     USE disvert_mod, ONLY : ap,bp
     USE omp_para
-    IMPLICIT NONE
     REAL(rstd),INTENT(IN) :: u(iim*3*jjm,llm), ps(iim*jjm)
     REAL(rstd),INTENT(OUT):: w(iim*jjm,llm)
     REAL(rstd):: convm(iim*jjm,llm+1)
     REAL(rstd):: p(iim*jjm,llm+1), rhodz(iim*jjm,llm), Fe(iim*3*jjm,llm)
+    REAL(rstd):: gradps(3*iim*jjm)
     REAL(rstd):: ugradps
     
@@ -57 +129 @@
          ENDDO
       ENDDO
+
+      !DIR$ SIMD
+      DO ij=ij_begin_ext, ij_end_ext
+         gradps(ij+u_right) = (ps(ij)-ps(ij+t_right))*ne_right*le(ij+u_right)
+         gradps(ij+u_lup)   = (ps(ij)-ps(ij+t_lup))  *ne_lup  *le(ij+u_lup)
+         gradps(ij+u_ldown) = (ps(ij)-ps(ij+t_ldown))*ne_ldown*le(ij+u_ldown)
+      END DO
     
 !!!  Compute mass flux
@@ -99 +178 @@
       convm(:,llm+1)=0.
 
-!!! Compute dps
-    !  DO j=jj_begin,jj_end
-    !    DO i=ii_begin,ii_end
-    !      ij=(j-1)*iim+i
-    !      ! dps/dt = -int(div flux)dz
-    !      dps(ij)=-convm(ij,1) * g 
-    !      convm(ij,llm+1)=0.
-    !    ENDDO
-    !  ENDDO
-    
-    ! Compute Omega = Dp/Dt
-    ! with p = A(eta)+B(eta)ps
-    !      Dp/Dt = dp/deta.Deta/Dt + B(eta)Dps/Dt
-    !            = -mg.Deta/Dt + B.Dps/Dt
-    ! By definition the mass flux through model levels is W=m.Deta/Dt with m=-1/g dp/deta
-    ! therefore 
-    !  Dp/Dt = -g.W + B.dps/dt + Bu.grad ps
-    !        = B.u.grad ps - g*convm
-    
-!!! Compute vertical flux through model layers
-    !  DO l = 1,llm-1
-    !    DO j=jj_begin,jj_end
-    !      DO i=ii_begin,ii_end
-    !        ij=(j-1)*iim+i
-    !        ! w = int(z,ztop,div(flux)dz) + B(eta)dps/dt
-    !        ! => w>0 for upward transport
-    !        w( ij, l+1 ) = convm( ij, l+1 ) - bp(l+1) * convm( ij, 1 ) ! g.W = g.convm + B dps/dt
-    !      ENDDO
-    !    ENDDO
-    !  ENDDO
-    
-    
 !!! Compute omega
-    ! -grad ps : ( ne(ij,ldown)*ps(ij,l) + ne(ij+t_ldown,rup)*ps(ij+t_ldown,l) ) ) / de(ij+u_ldown)
 
        DO l = 1,llm
@@ -139 +185 @@
              ij=(j-1)*iim+i
              ugradps = &
-                  le(ij+u_right)*u(ij+u_right,l)*( ne(ij,right)*ps(ij) + ne(ij+t_right,left)*ps(ij+t_right) )  &
-                  + le(ij+u_rup)*u(ij+u_rup,l)*( ne(ij,rup)*ps(ij) + ne(ij+t_rup,ldown)*ps(ij+t_rup) )  &
-                  + le(ij+u_lup)*u(ij+u_lup,l)*( ne(ij,lup)*ps(ij) + ne(ij+t_lup,rdown)*ps(ij+t_lup) )         &
-                  + le(ij+u_left)*u(ij+u_left,l)*( ne(ij,left)*ps(ij) + ne(ij+t_left,right)*ps(ij+t_left) )    &
-                  + le(ij+u_ldown)*u(ij+u_ldown,l)*( ne(ij,ldown)*ps(ij) + ne(ij+t_ldown,rup)*ps(ij+t_ldown) ) &
-                  + le(ij+u_rdown)*u(ij+u_rdown,l)*( ne(ij,rdown)*ps(ij) + ne(ij+t_rdown,lup)*ps(ij+t_rdown) )
+                  u(ij+u_rup,l)*gradps(ij+u_rup)       &
+                  + u(ij+u_lup,l)*gradps(ij+u_lup)     &
+                  + u(ij+u_left,l)*gradps(ij+u_left)   &
+                  + u(ij+u_ldown,l)*gradps(ij+u_ldown) &
+                  + u(ij+u_rdown,l)*gradps(ij+u_rdown) &
+                  + u(ij+u_right,l)*gradps(ij+u_right)
+
              ugradps = .5*(bp(l)+bp(l+1)) *ugradps/(-4*Ai(ij))  ! sign convention as in Ringler et al. 2010, Eq. 22 p.3072
              w( ij, l) =  ugradps - g*.5*(convm( ij,l+1)+convm(ij,l)) 
@@ -155 +202 @@
   END SUBROUTINE compute_omega
   
-END MODULE omega_mod
+END MODULE compute_omega_mod

codes/icosagcm/devel/src/diagnostics/observable.f90

@@ -48 +48 @@
     USE vertical_interp_mod
     USE theta2theta_rhodz_mod
-    USE omega_mod
+    USE compute_omega_mod, ONLY : w_omega
     USE kinetic_mod
     USE grid_param