Changeset 538 for codes/icosagcm/devel/src/dynamics

Timestamp:

06/09/17 16:13:47 (7 years ago)

Author:

dubos

Message:

devel : macro-generated kernels for caldyn_coriolis, compute_NH_geopot, caldyn_solver, caldyn_vert_NH

Location:

codes/icosagcm/devel/src/dynamics

Files:

• caldyn_kernels_base.F90 (modified) (6 diffs)
• caldyn_kernels_hevi.F90 (modified) (13 diffs)

codes/icosagcm/devel/src/dynamics/caldyn_kernels_base.F90

@@ -49 +49 @@
 
 #ifdef CPP_DYSL
+!#if 0
 #include "../kernels/compute_geopot.k90"
 #else
@@ -303 +304 @@
     REAL(rstd),INTENT(INOUT) :: dW(iim*jjm,llm+1)
     ! local arrays
-    REAL(rstd) :: eta_dot(iim*jjm) ! eta_dot in full layers
-    REAL(rstd) :: wcov(iim*jjm) ! covariant vertical momentum in full layers
+    REAL(rstd) :: eta_dot(iim*jjm, llm) ! eta_dot in full layers
+    REAL(rstd) :: wcov(iim*jjm,llm) ! covariant vertical momentum in full layers
     REAL(rstd) :: W_etadot(iim*jjm,llm) ! vertical flux of vertical momentum
     ! indices and temporary values
@@ -312 +313 @@
     CALL trace_start("compute_caldyn_vert_nh")
 
+#ifdef CPP_DYSL
+!#if 0
+#include "../kernels/caldyn_vert_NH.k90"
+#else
+#define ETA_DOT(ij) eta_dot(ij,1)
+#define WCOV(ij) wcov(ij,1)
+
     DO l=ll_begin,ll_end
        ! compute the local arrays
@@ -319 +327 @@
           w_ij = .5*(W(ij,l)+W(ij,l+1))/mass(ij,l)
           W_etadot(ij,l) = wflux_ij*w_ij
-          eta_dot(ij) = wflux_ij / mass(ij,l)
-          wcov(ij) = w_ij*(geopot(ij,l+1)-geopot(ij,l))
+          ETA_DOT(ij) = wflux_ij / mass(ij,l)
+          WCOV(ij) = w_ij*(geopot(ij,l+1)-geopot(ij,l))
        ENDDO
        ! add NH term to du
@@ -326 +334 @@
       DO ij=ij_begin,ij_end
           du(ij+u_right,l) = du(ij+u_right,l) &
-               - .5*(wcov(ij+t_right)+wcov(ij)) &
-               *ne_right*(eta_dot(ij+t_right)-eta_dot(ij))
+               - .5*(WCOV(ij+t_right)+WCOV(ij)) &
+               *ne_right*(ETA_DOT(ij+t_right)-ETA_DOT(ij))
           du(ij+u_lup,l) = du(ij+u_lup,l) &
-               - .5*(wcov(ij+t_lup)+wcov(ij)) &
-               *ne_lup*(eta_dot(ij+t_lup)-eta_dot(ij))
+               - .5*(WCOV(ij+t_lup)+WCOV(ij)) &
+               *ne_lup*(ETA_DOT(ij+t_lup)-ETA_DOT(ij))
           du(ij+u_ldown,l) = du(ij+u_ldown,l) &
-               - .5*(wcov(ij+t_ldown)+wcov(ij)) &
-               *ne_ldown*(eta_dot(ij+t_ldown)-eta_dot(ij))
+               - .5*(WCOV(ij+t_ldown)+WCOV(ij)) &
+               *ne_ldown*(ETA_DOT(ij+t_ldown)-ETA_DOT(ij))
        END DO
     ENDDO
@@ -352 +360 @@
        END DO
     END DO
+
+#undef ETA_DOT
+#undef WCOV
+#endif
+
     CALL trace_end("compute_caldyn_vert_nh")
 

codes/icosagcm/devel/src/dynamics/caldyn_kernels_hevi.F90

@@ -11 +11 @@
   REAL(rstd), PARAMETER :: pbot=1e5, Phi_bot=0., rho_bot=1e6 ! FIXME
 
-  LOGICAL, PARAMETER :: debug_hevi_solver = .FALSE.
+  LOGICAL, SAVE :: debug_hevi_solver = .FALSE.
+  LOGICAL, PARAMETER :: rigid=.TRUE.
 
   PUBLIC :: compute_theta, compute_pvort_only, compute_caldyn_Coriolis, &
@@ -117 +118 @@
     REAL(rstd) :: wil, tau2_g, g2, gm2, ml_g2, c2_mik
 
-    INTEGER    :: iter, ij, l
+    INTEGER    :: iter, ij, l, ij_omp_begin_ext, ij_omp_end_ext
+
+    CALL distrib_level(ij_end_ext-ij_begin_ext+1,ij_omp_begin_ext,ij_omp_end_ext)
+    ij_omp_begin_ext=ij_omp_begin_ext+ij_begin_ext-1
+    ij_omp_end_ext=ij_omp_end_ext+ij_begin_ext-1
+
+#ifdef CPP_DYSL
+!#if 0
+#include "../kernels/compute_NH_geopot.k90"
+#else
 
 !    FIXME : vertical OpenMP parallelism will not work
@@ -251 +261 @@
 
     END DO ! Newton-Raphson
+
+#endif
     
   END SUBROUTINE compute_NH_geopot
@@ -257 +269 @@
     REAL(rstd),INTENT(IN) :: tau ! "solve" Phi-tau*dPhi/dt = Phi_rhs
     REAL(rstd),INTENT(IN)    :: rhodz(iim*jjm,llm)
-    REAL(rstd),INTENT(IN)    :: theta(iim*jjm,llm)
+    REAL(rstd),INTENT(IN)    :: theta(iim*jjm,llm,nqdyn)
     REAL(rstd),INTENT(OUT)   :: pk(iim*jjm,llm)
     REAL(rstd),INTENT(INOUT) :: geopot(iim*jjm,llm+1)
@@ -266 +278 @@
 
     REAL(rstd) :: m_il(iim*jjm,llm+1)        ! rhodz averaged to interfaces
-    REAL(rstd) :: berni(iim*jjm)             ! (W/m_il)^2
-    REAL(rstd) :: gamma, rho_ij, X_ij, Y_ij
+    REAL(rstd) :: pres(iim*jjm,llm)          ! pressure
+    REAL(rstd) :: berni(iim*jjm,llm)             ! (W/m_il)^2
+    REAL(rstd) :: gamma, rho_ij, T_ij, X_ij, Y_ij, vreff, Rd, Cvd
     INTEGER    :: ij, l
 
     CALL trace_start("compute_caldyn_solver")
 
+    Rd=cpp*kappa
+
+#ifdef CPP_DYSL
+!#if 0
+#include "../kernels/caldyn_solver.k90"
+#else
+#define BERNI(ij) berni(ij,1)
     ! FIXME : vertical OpenMP parallelism will not work
 
@@ -303 +323 @@
        DO ij=ij_begin_ext,ij_end_ext
           rho_ij = (g*rhodz(ij,l))/(geopot(ij,l+1)-geopot(ij,l))
-          X_ij = (cpp/preff)*kappa*theta(ij,l)*rho_ij
+          X_ij = (cpp/preff)*kappa*theta(ij,l,1)*rho_ij
           ! kappa.theta.rho = p/exner 
           ! => X = (p/p0)/(exner/Cp)
@@ -342 +362 @@
        DO ij=ij_begin_ext,ij_end_ext
           pk(ij,l) = cpp*((pk(ij,l)/preff)**kappa) ! other formulae possible if exponentiation is slow
-          berni(ij) = (-.25*g*g)*(              &
+          BERNI(ij) = (-.25*g*g)*(              &
                  (W(ij,l)/m_il(ij,l))**2       &
                + (W(ij,l+1)/m_il(ij,l+1))**2 )
        ENDDO
        DO ij=ij_begin,ij_end         
-          du(ij+u_right,l) = ne_right*(berni(ij)-berni(ij+t_right))
-          du(ij+u_lup,l)   = ne_lup  *(berni(ij)-berni(ij+t_lup))
-          du(ij+u_ldown,l) = ne_ldown*(berni(ij)-berni(ij+t_ldown))
+          du(ij+u_right,l) = ne_right*(BERNI(ij)-BERNI(ij+t_right))
+          du(ij+u_lup,l)   = ne_lup  *(BERNI(ij)-BERNI(ij+t_lup))
+          du(ij+u_ldown,l) = ne_ldown*(BERNI(ij)-BERNI(ij+t_ldown))
        ENDDO
     ENDDO
-    
+#undef BERNI
+#endif
+
     CALL trace_end("compute_caldyn_solver")
     
@@ -488 +510 @@
     REAL(rstd),INTENT(INOUT) :: du(3*iim*jjm,llm)
     
-    REAL(rstd) :: Ftheta(3*iim*jjm)  ! potential temperature flux
-    REAL(rstd) :: uu_right, uu_lup, uu_ldown
+    REAL(rstd) :: Ftheta(3*iim*jjm,llm)  ! potential temperature flux
+    REAL(rstd) :: uu_right, uu_lup, uu_ldown, du_trisk, divF
     INTEGER :: ij,iq,l,kdown
 
     CALL trace_start("compute_caldyn_Coriolis")
+
+#ifdef CPP_DYSL
+!#if 0
+#include "../kernels/coriolis.k90"
+#else
+#define FTHETA(ij) Ftheta(ij,1)
 
     DO l=ll_begin, ll_end
@@ -499 +527 @@
        !DIR$ SIMD
           DO ij=ij_begin_ext,ij_end_ext      
-             Ftheta(ij+u_right) = 0.5*(theta(ij,l,iq)+theta(ij+t_right,l,iq)) &
+             FTHETA(ij+u_right) = 0.5*(theta(ij,l,iq)+theta(ij+t_right,l,iq)) &
                                   * hflux(ij+u_right,l)
-             Ftheta(ij+u_lup)   = 0.5*(theta(ij,l,iq)+theta(ij+t_lup,l,iq)) &
+             FTHETA(ij+u_lup)   = 0.5*(theta(ij,l,iq)+theta(ij+t_lup,l,iq)) &
                   * hflux(ij+u_lup,l)
-             Ftheta(ij+u_ldown) = 0.5*(theta(ij,l,iq)+theta(ij+t_ldown,l,iq)) &
+             FTHETA(ij+u_ldown) = 0.5*(theta(ij,l,iq)+theta(ij+t_ldown,l,iq)) &
                   * hflux(ij+u_ldown,l)
           END DO
@@ -511 +539 @@
              ! dtheta_rhodz =  -div(flux.theta)
              dtheta_rhodz(ij,l,iq)= &
-                  -1./Ai(ij)*(ne_right*Ftheta(ij+u_right)    +  &
-                  ne_rup*Ftheta(ij+u_rup)        +  &  
-                  ne_lup*Ftheta(ij+u_lup)        +  &  
-                  ne_left*Ftheta(ij+u_left)      +  &  
-                  ne_ldown*Ftheta(ij+u_ldown)    +  &
-                  ne_rdown*Ftheta(ij+u_rdown) )
+                  -1./Ai(ij)*(ne_right*FTHETA(ij+u_right)    +  &
+                  ne_rup*FTHETA(ij+u_rup)        +  &  
+                  ne_lup*FTHETA(ij+u_lup)        +  &  
+                  ne_left*FTHETA(ij+u_left)      +  &  
+                  ne_ldown*FTHETA(ij+u_ldown)    +  &
+                  ne_rdown*FTHETA(ij+u_rdown) )
           END DO
        END DO
@@ -626 +654 @@
        STOP
     END SELECT
+#undef FTHETA
+#endif
 
     CALL trace_end("compute_caldyn_Coriolis")
@@ -645 +675 @@
 
 #ifdef CPP_DYSL
+!#if 0
 #define BERNI(ij,l) berni(ij,l)
 #include "../kernels/caldyn_slow_hydro.k90"
+#undef BERNI
 #else
 #define BERNI(ij) berni(ij,1)
@@ -696 +728 @@
        END IF
     END DO
+#undef BERNI
 #endif
-#undef BERNI
     CALL trace_end("compute_caldyn_slow_hydro")    
   END SUBROUTINE compute_caldyn_slow_hydro