Changeset 1015 for palm/trunk/SOURCE/buoyancy.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/buoyancy.f90 (modified) (4 diffs)

palm/trunk/SOURCE/buoyancy.f90

@@ -4 +4 @@
 ! Currrent revisions:
 ! -----------------
-!
+! accelerator version (*_acc) added
 !
 ! Former revisions:
@@ -55 +55 @@
 
     PRIVATE
-    PUBLIC buoyancy, calc_mean_profile
+    PUBLIC buoyancy, buoyancy_acc, calc_mean_profile
 
     INTERFACE buoyancy
@@ -61 +61 @@
        MODULE PROCEDURE buoyancy_ij
     END INTERFACE buoyancy
+
+    INTERFACE buoyancy_acc
+       MODULE PROCEDURE buoyancy_acc
+    END INTERFACE buoyancy_acc
 
     INTERFACE calc_mean_profile
@@ -164 +168 @@
 
     END SUBROUTINE buoyancy
+
+
+!------------------------------------------------------------------------------!
+! Call for all grid points - accelerator version
+!------------------------------------------------------------------------------!
+    SUBROUTINE buoyancy_acc( var, var_reference, wind_component, pr )
+
+       USE arrays_3d
+       USE control_parameters
+       USE indices
+       USE pegrid
+       USE statistics
+
+       IMPLICIT NONE
+
+       INTEGER ::  i, j, k, pr, wind_component
+       REAL    ::  var_reference
+#if defined( __nopointer )
+       REAL, DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  var
+#else
+       REAL, DIMENSION(:,:,:), POINTER ::  var
+#endif
+
+
+       IF ( .NOT. sloping_surface )  THEN
+!
+!--       Normal case: horizontal surface
+          IF ( use_reference )  THEN
+             DO  i = nxl, nxr
+                DO  j = nys, nyn
+                   DO  k = nzb_s_inner(j,i)+1, nzt-1
+                      tend(k,j,i) = tend(k,j,i) + atmos_ocean_sign * g * 0.5 * &
+                                                            (                  &
+                          ( var(k,j,i)   - hom(k,1,pr,0)   ) / var_reference + &
+                          ( var(k+1,j,i) - hom(k+1,1,pr,0) ) / var_reference   &
+                                                            )
+                   ENDDO
+                ENDDO
+             ENDDO
+          ELSE
+             !$acc kernels present( nzb_s_inner, hom, tend, var )
+             !$acc loop
+             DO  i = nxl, nxr
+                DO  j = nys, nyn
+                   !$acc loop vector(32)
+                   DO  k = 1, nzt-1
+                      IF ( k > nzb_s_inner(j,i) )  THEN
+                         tend(k,j,i) = tend(k,j,i) + atmos_ocean_sign * g * 0.5 * &
+                                                               (                  &
+                             ( var(k,j,i)   - hom(k,1,pr,0)   ) / hom(k,1,pr,0) + &
+                             ( var(k+1,j,i) - hom(k+1,1,pr,0) ) / hom(k+1,1,pr,0) &
+                                                               )
+                      ENDIF
+                   ENDDO
+                ENDDO
+             ENDDO
+             !$acc end kernels
+          ENDIF
+
+       ELSE
+!
+!--       Buoyancy term for a surface with a slope in x-direction. The equations
+!--       for both the u and w velocity-component contain proportionate terms.
+!--       Temperature field at time t=0 serves as environmental temperature.
+!--       Reference temperature (pt_surface) is the one at the lower left corner
+!--       of the total domain.
+          IF ( wind_component == 1 )  THEN
+
+             DO  i = nxlu, nxr
+                DO  j = nys, nyn
+                   DO  k = nzb_s_inner(j,i)+1, nzt-1
+                      tend(k,j,i) = tend(k,j,i) + g * sin_alpha_surface *      &
+                           0.5 * ( ( pt(k,j,i-1)         + pt(k,j,i)         ) &
+                                 - ( pt_slope_ref(k,i-1) + pt_slope_ref(k,i) ) &
+                                 ) / pt_surface
+                   ENDDO
+                ENDDO
+             ENDDO
+
+          ELSEIF ( wind_component == 3 )  THEN
+
+             DO  i = nxl, nxr
+                DO  j = nys, nyn
+                   DO  k = nzb_s_inner(j,i)+1, nzt-1
+                      tend(k,j,i) = tend(k,j,i) + g * cos_alpha_surface *      &
+                           0.5 * ( ( pt(k,j,i)         + pt(k+1,j,i)         ) &
+                                 - ( pt_slope_ref(k,i) + pt_slope_ref(k+1,i) ) &
+                                 ) / pt_surface
+                   ENDDO
+                ENDDO
+            ENDDO
+
+          ELSE
+
+             WRITE( message_string, * ) 'no term for component "',&
+                                       wind_component,'"'
+             CALL message( 'buoyancy', 'PA0159', 1, 2, 0, 6, 0 )
+
+          ENDIF
+
+       ENDIF
+
+    END SUBROUTINE buoyancy_acc