Changeset 1015 for palm/trunk/SOURCE/diffusivities.f90

Timestamp:

Sep 27, 2012 9:23:24 AM (12 years ago)

Author:

raasch

Message:

Starting with changes required for GPU optimization. OpenACC statements for using NVIDIA GPUs added.
Adjustment of mixing length to the Prandtl mixing length at first grid point above ground removed.
mask array is set zero for ghost boundaries

File:

• palm/trunk/SOURCE/diffusivities.f90 (modified) (5 diffs)

palm/trunk/SOURCE/diffusivities.f90

@@ -4 +4 @@
 ! Current revisions:
 ! -----------------
+! OpenACC statements added + code changes required for GPU optimization,
+! adjustment of mixing length to the Prandtl mixing length at first grid point
+! above ground removed
 ! 
 ! Former revisions:
@@ -52 +55 @@
     INTEGER ::  i, j, k, omp_get_thread_num, sr, tn
 
-    REAL    ::  dvar_dz, l_stable, var_reference
-
-    REAL, SAVE ::  phi_m = 1.0
+    REAL    ::  dvar_dz, l, ll, l_stable, sqrt_e, var_reference
 
     REAL    ::  var(nzb:nzt+1,nysg:nyng,nxlg:nxrg)
-
-    REAL, DIMENSION(1:nzt) ::  l, ll, sqrt_e
 
 
@@ -71 +70 @@
 !
 !-- Compute the turbulent diffusion coefficient for momentum
-    !$OMP PARALLEL PRIVATE (dvar_dz,i,j,k,l,ll,l_stable,phi_m,sqrt_e,sr,tn)
+    !$OMP PARALLEL PRIVATE (dvar_dz,i,j,k,l,ll,l_stable,sqrt_e,sr,tn)
 !$  tn = omp_get_thread_num()
 
+!
+!-- Data declerations for accelerators
+    !$acc data present( dd2zu, e, km, kh, l_grid, l_wall, nzb_s_inner, rif, var )
+    !$acc kernels
+
+!
+!-- Introduce an optional minimum tke
+    IF ( e_min > 0.0 )  THEN
+       !$OMP DO
+       !$acc loop
+       DO  i = nxlg, nxrg
+          DO  j = nysg, nyng
+             !$acc loop vector( 32 )
+             DO  k = 1, nzt
+                IF ( k > nzb_s_inner(j,i) )  THEN
+                   e(k,j,i) = MAX( e(k,j,i), e_min )
+                ENDIF
+             ENDDO
+          ENDDO
+       ENDDO
+    ENDIF
+
     !$OMP DO
+    !$acc loop
     DO  i = nxlg, nxrg
        DO  j = nysg, nyng
+          !$acc loop vector( 32 )
+          DO  k = 1, nzt
 
+             IF ( k > nzb_s_inner(j,i) )  THEN
+
+                sqrt_e = SQRT( e(k,j,i) )
 !
-!--       Compute the Phi-function for a possible adaption of the mixing length
-!--       to the Prandtl mixing length
-          IF ( adjust_mixing_length  .AND.  prandtl_layer )  THEN
-             IF ( rif(j,i) >= 0.0 )  THEN
-                phi_m = 1.0 + 5.0 * rif(j,i)
-             ELSE
-                phi_m = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rif(j,i) ) )
-             ENDIF
-          ENDIF
-          
+!--             Determine the mixing length
+                dvar_dz = atmos_ocean_sign * &  ! inverse effect of pt/rho gradient
+                          ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
+                IF ( dvar_dz > 0.0 ) THEN
+                   IF ( use_reference )  THEN
+                      l_stable = 0.76 * sqrt_e / &
+                                 SQRT( g / var_reference * dvar_dz ) + 1E-5
+                   ELSE
+                      l_stable = 0.76 * sqrt_e / &
+                                 SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5
+                   ENDIF
+                ELSE
+                   l_stable = l_grid(k)
+                ENDIF
 !
-!--       Introduce an optional minimum tke
-          IF ( e_min > 0.0 )  THEN
-             DO  k = nzb_s_inner(j,i)+1, nzt
-                e(k,j,i) = MAX( e(k,j,i), e_min )
-             ENDDO
-          ENDIF
+!--             Adjustment of the mixing length
+                IF ( wall_adjustment )  THEN
+                   l  = MIN( l_wall(k,j,i), l_grid(k), l_stable )
+                   ll = MIN( l_wall(k,j,i), l_grid(k) )
+                ELSE
+                   l  = MIN( l_grid(k), l_stable )
+                   ll = l_grid(k)
+                ENDIF
 
-!
-!--       Calculate square root of e in a seperate loop, because it is used
-!--       twice in the next loop (better vectorization)
-          DO  k = nzb_s_inner(j,i)+1, nzt
-             sqrt_e(k) = SQRT( e(k,j,i) )
-          ENDDO
+      !
+      !--       Compute diffusion coefficients for momentum and heat
+                km(k,j,i) = 0.1 * l * sqrt_e
+                kh(k,j,i) = ( 1.0 + 2.0 * l / ll ) * km(k,j,i)
 
-!
-!--       Determine the mixing length
-          DO  k = nzb_s_inner(j,i)+1, nzt
-             dvar_dz = atmos_ocean_sign * &  ! inverse effect of pt/rho gradient
-                       ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
-             IF ( dvar_dz > 0.0 ) THEN
-                IF ( use_reference )  THEN
-                   l_stable = 0.76 * sqrt_e(k) / &
-                                     SQRT( g / var_reference * dvar_dz ) + 1E-5
-                ELSE
-                   l_stable = 0.76 * sqrt_e(k) / &
-                                     SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5
-                ENDIF
-             ELSE
-                l_stable = l_grid(k)
-             ENDIF
-!
-!--          Adjustment of the mixing length
-             IF ( wall_adjustment )  THEN
-                l(k)  = MIN( l_wall(k,j,i), l_grid(k), l_stable )
-                ll(k) = MIN( l_wall(k,j,i), l_grid(k) )
-             ELSE
-                l(k)  = MIN( l_grid(k), l_stable )
-                ll(k) = l_grid(k)
-             ENDIF
-             IF ( adjust_mixing_length  .AND.  prandtl_layer )  THEN
-                l(k)  = MIN( l(k),  kappa * &
-                                    ( zu(k) - zw(nzb_s_inner(j,i)) ) / phi_m )
-                ll(k) = MIN( ll(k), kappa * &
-                                    ( zu(k) - zw(nzb_s_inner(j,i)) ) / phi_m )
              ENDIF
 
+#if ! defined( __openacc )
 !
-!--          Compute diffusion coefficients for momentum and heat
-             km(k,j,i) = 0.1 * l(k) * sqrt_e(k)
-             kh(k,j,i) = ( 1.0 + 2.0 * l(k) / ll(k) ) * km(k,j,i)
+!++          Statistics still have to be realized for accelerators
+!--          Summation for averaged profile (cf. flow_statistics)
+             DO  sr = 0, statistic_regions
+                sums_l_l(k,sr,tn) = sums_l_l(k,sr,tn) + l * rmask(j,i,sr)
+             ENDDO
+#endif
 
           ENDDO
-
-!
-!--       Summation for averaged profile (cf. flow_statistics)
-!--       (the IF statement still requires a performance check on NEC machines)
-          DO  sr = 0, statistic_regions
-             IF ( rmask(j,i,sr) /= 0.0 .AND.  &
-                  i >= nxl .AND. i <= nxr .AND. j >= nys .AND. j <= nyn )  THEN
-                DO  k = nzb_s_inner(j,i)+1, nzt
-                   sums_l_l(k,sr,tn) = sums_l_l(k,sr,tn) + l(k)
-                ENDDO
-             ENDIF
-          ENDDO
-
        ENDDO
     ENDDO
 
+#if ! defined( __openacc )
+!
+!++ Statistics still have to be realized for accelerators
     sums_l_l(nzt+1,:,tn) = sums_l_l(nzt,:,tn)   ! quasi boundary-condition for
                                                   ! data output
-
+#endif
     !$OMP END PARALLEL
 
@@ -169 +164 @@
 !-- values of the diffusivities are not needed
     !$OMP PARALLEL DO
+    !$acc loop
     DO  i = nxlg, nxrg
        DO  j = nysg, nyng
@@ -198 +194 @@
     ENDIF
 
+    !$acc end kernels
+    !$acc end data
 
  END SUBROUTINE diffusivities