Changeset 558

Timestamp:

09/21/17 18:18:47 (7 years ago)

Author:

dubos

Message:

devel : Fix OpenMP for NH

Location:

codes/icosagcm/devel/src/dynamics

Files:

• caldyn_gcm.F90 (modified) (2 diffs)
• caldyn_hevi.f90 (modified) (4 diffs)
• caldyn_kernels_base.F90 (modified) (4 diffs)
• caldyn_kernels_hevi.F90 (modified) (9 diffs)

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

@@ -22 +22 @@
     REAL(rstd),POINTER :: planetvel(:)
 
-    #ifdef CPP_DYSL
+#ifdef CPP_DYSL
        IF(is_master) PRINT *,'CPP_DYSL : Using macro-generated compute kernels'
-    #endif
+#endif
 
     hydrostatic=.TRUE.
@@ -122 +122 @@
     CALL allocate_field(f_qv,field_z,type_real,llm) 
     CALL allocate_field(f_pk,    field_t,type_real,llm,  name='pk')
-    CALL allocate_field(f_wwuu,  field_u,type_real,llm+1,name='wwuu')
+    CALL allocate_field(f_wwuu,  field_u,type_real,llm+1,name='wwuu')    
     CALL allocate_field(f_planetvel,  field_u,type_real, name='planetvel') ! planetary velocity at r=a
-
+    IF(.NOT.hydrostatic) THEN
+       CALL allocate_field(f_Fel,      field_u,type_real,llm+1,name='F_el')
+       CALL allocate_field(f_gradPhi2, field_t,type_real,llm+1,name='gradPhi2')
+       CALL allocate_field(f_wil,      field_t,type_real,llm+1,name='w_il')
+       CALL allocate_field(f_Wetadot,  field_t,type_real,llm,name='W_etadot')
+    END IF
   END SUBROUTINE allocate_caldyn
 

codes/icosagcm/devel/src/dynamics/caldyn_hevi.f90

@@ -59 +59 @@
     REAL(rstd),POINTER  :: convm(:,:) 
     REAL(rstd),POINTER  :: wwuu(:,:)
-        
+    REAL(rstd),POINTER  :: F_el(:,:), gradPhi2(:,:), w_il(:,:) , W_etadot(:,:), pres(:,:), m_il(:,:)
     INTEGER :: ind
     LOGICAL,SAVE :: first=.TRUE.
@@ -118 +118 @@
           dW = f_dW_fast(ind)
           dPhi = f_dPhi_fast(ind)
-          CALL compute_caldyn_solver(tau,mass,theta,pk,geopot,W,dPhi,dW,du) ! computes d(Phi,W,du)_fast and updates Phi,W
+          ! reuse buffers
+          m_il = f_wil(ind)
+          pres = f_gradPhi2(ind)
+          CALL compute_caldyn_solver(tau,mass,theta,pk,geopot,W, m_il,pres, dPhi,dW,du) ! computes d(Phi,W,du)_fast and updates Phi,W
        END IF
        u=f_u(ind)
@@ -161 +164 @@
           geopot = f_geopot(ind)
           dPhi = f_dPhi_slow(ind)
-          CALL compute_caldyn_slow_NH(u,mass,geopot,W, hflux,du,dPhi,dW)
+          F_el = f_Fel(ind)
+          gradPhi2 = f_gradPhi2(ind)
+          w_il = f_wil(ind)
+          CALL compute_caldyn_slow_NH(u,mass,geopot,W, F_el,gradPhi2,w_il, hflux,du,dPhi,dW)
        END IF
        CALL compute_caldyn_Coriolis(hflux,theta,qu, convm,dtheta_rhodz,du)
@@ -168 +174 @@
           wwuu=f_wwuu(ind)
           dps=f_dps(ind)
+          W_etadot=f_Wetadot(ind)
           CALL compute_caldyn_vert(u,theta,mass,convm, wflux,wwuu, dps, dtheta_rhodz, du)
           IF(.NOT.hydrostatic) THEN
-             CALL compute_caldyn_vert_NH(mass,geopot,W,wflux, du,dPhi,dW)
+             CALL compute_caldyn_vert_NH(mass,geopot,W,wflux, W_etadot, du,dPhi,dW)
           END IF
        END IF

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

@@ -14 +14 @@
 
   ! temporary shared variables for caldyn
-  TYPE(t_field),POINTER,PUBLIC :: f_pk(:),f_wwuu(:),f_planetvel(:)
+  TYPE(t_field),POINTER,PUBLIC :: f_pk(:),f_wwuu(:),f_planetvel(:), &
+                                  f_Fel(:), f_gradPhi2(:), f_wil(:), f_Wetadot(:)
 
   INTEGER, PUBLIC :: caldyn_conserv
@@ -295 +296 @@
   END SUBROUTINE compute_caldyn_vert
 
-  SUBROUTINE compute_caldyn_vert_NH(mass,geopot,W,wflux, du,dPhi,dW)
+  SUBROUTINE compute_caldyn_vert_NH(mass,geopot,W,wflux, W_etadot, du,dPhi,dW)
     REAL(rstd),INTENT(IN) :: mass(iim*jjm,llm)
     REAL(rstd),INTENT(IN) :: geopot(iim*jjm,llm+1)
@@ -303 +304 @@
     REAL(rstd),INTENT(INOUT) :: dPhi(iim*jjm,llm+1)
     REAL(rstd),INTENT(INOUT) :: dW(iim*jjm,llm+1)
+    REAL(rstd) :: W_etadot(iim*jjm,llm) ! vertical flux of vertical momentum
     ! local arrays
     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
     INTEGER    :: ij, l
@@ -314 +315 @@
 
 #ifdef CPP_DYSL
-!#if 0
+!$OMP BARRIER
 #include "../kernels/caldyn_vert_NH.k90"
+!$OMP BARRIER
 #else
 #define ETA_DOT(ij) eta_dot(ij,1)

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

@@ -125 +125 @@
 
 #ifdef CPP_DYSL
-!#if 0
 #include "../kernels/compute_NH_geopot.k90"
 #else
-
 !    FIXME : vertical OpenMP parallelism will not work
    
@@ -266 +264 @@
   END SUBROUTINE compute_NH_geopot
 
-  SUBROUTINE compute_caldyn_solver(tau,rhodz,theta,pk, geopot,W, dPhi,dW,du)
+  SUBROUTINE compute_caldyn_solver(tau,rhodz,theta,pk, geopot,W, m_il,pres, dPhi,dW,du)
     REAL(rstd),INTENT(IN) :: tau ! "solve" Phi-tau*dPhi/dt = Phi_rhs
     REAL(rstd),INTENT(IN)    :: rhodz(iim*jjm,llm)
@@ -273 +271 @@
     REAL(rstd),INTENT(INOUT) :: geopot(iim*jjm,llm+1)
     REAL(rstd),INTENT(INOUT) :: W(iim*jjm,llm+1) ! OUT if tau>0
+    REAL(rstd),INTENT(OUT)   :: m_il(iim*jjm,llm+1)        ! rhodz averaged to interfaces
+    REAL(rstd),INTENT(OUT)   :: pres(iim*jjm,llm)          ! pressure
     REAL(rstd),INTENT(OUT)   :: dW(iim*jjm,llm+1)
     REAL(rstd),INTENT(OUT)   :: dPhi(iim*jjm,llm+1)
     REAL(rstd),INTENT(OUT)   :: du(3*iim*jjm,llm)
 
-    REAL(rstd) :: m_il(iim*jjm,llm+1)        ! rhodz averaged to interfaces
-    REAL(rstd) :: pres(iim*jjm,llm)          ! pressure
-    REAL(rstd) :: berni(iim*jjm,llm)             ! (W/m_il)^2
+    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
@@ -288 +286 @@
 
 #ifdef CPP_DYSL
-!#if 0
+!$OMP BARRIER
 #include "../kernels/caldyn_solver.k90"
+!$OMP BARRIER
 #else
 #define BERNI(ij) berni(ij,1)
@@ -517 +516 @@
 
 #ifdef CPP_DYSL
-!#if 0
 #include "../kernels/coriolis.k90"
 #else
@@ -675 +673 @@
 
 #ifdef CPP_DYSL
-!#if 0
 #define BERNI(ij,l) berni(ij,l)
 #include "../kernels/caldyn_slow_hydro.k90"
@@ -733 +730 @@
   END SUBROUTINE compute_caldyn_slow_hydro
 
-  SUBROUTINE compute_caldyn_slow_NH(u,rhodz,Phi,W, hflux,du,dPhi,dW)
+  SUBROUTINE compute_caldyn_slow_NH(u,rhodz,Phi,W, F_el,gradPhi2,w_il, hflux,du,dPhi,dW)
     REAL(rstd),INTENT(IN)  :: u(3*iim*jjm,llm)    ! prognostic "velocity"
     REAL(rstd),INTENT(IN)  :: rhodz(iim*jjm,llm)  ! rho*dz
@@ -744 +741 @@
     REAL(rstd),INTENT(OUT) :: dPhi(iim*jjm,llm+1)
 
-    REAL(rstd) :: w_il(3*iim*jjm,llm+1) ! Wil/mil
+    REAL(rstd) :: w_il(iim*jjm,llm+1) ! Wil/mil
     REAL(rstd) :: F_el(3*iim*jjm,llm+1) ! NH mass flux
-    REAL(rstd) :: GradPhi2(3*iim*jjm,llm+1) ! grad_Phi**2
+    REAL(rstd) :: gradPhi2(iim*jjm,llm+1) ! grad_Phi**2
     REAL(rstd) :: DePhil(3*iim*jjm,llm+1) ! grad(Phi)
     
@@ -765 +762 @@
 
 #ifdef CPP_DYSL
+!$OMP BARRIER
 #include "../kernels/caldyn_slow_NH.k90"
+!$OMP BARRIER
 #else
     DO l=ll_begin, ll_endp1 ! compute on l levels (interfaces)