Index: palm/trunk/SOURCE/buoyancy.f90
===================================================================
--- palm/trunk/SOURCE/buoyancy.f90	(revision 4901)
+++ palm/trunk/SOURCE/buoyancy.f90	(revision 1)
@@ -1,38 +1,73 @@
-!> @file buoyancy.f90
-!--------------------------------------------------------------------------------------------------!
-! This file is part of the PALM model system.
-!
-! PALM is free software: you can redistribute it and/or modify it under the terms of the GNU General
-! Public License as published by the Free Software Foundation, either version 3 of the License, or
-! (at your option) any later version.
-!
-! PALM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the
-! implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
-! Public License for more details.
-!
-! You should have received a copy of the GNU General Public License along with PALM. If not, see
-! <http://www.gnu.org/licenses/>.
-!
-!
-! Copyright 1997-2021 Leibniz Universitaet Hannover
-!--------------------------------------------------------------------------------------------------!
-!
-! Current revisions:
-! ------------------
+ MODULE buoyancy_mod
+
+!------------------------------------------------------------------------------!
+! Actual revisions:
+! -----------------
 ! 
-! 
+!
 ! Former revisions:
 ! -----------------
-! $Id$
-! file re-formatted to follow the PALM coding standard
-!
-! 4360 2020-01-07 11:25:50Z suehring
-! Corrected "Former revisions" section
-!
-! 3725 2019-02-07 10:11:02Z raasch
-! unused variables removed
-!
-! 3655 2019-01-07 16:51:22Z knoop
-! nopointer option removed
+! $Log: buoyancy.f90,v $
+! Revision 1.19  2006/04/26 12:09:56  raasch
+! OpenMP optimization (one dimension added to sums_l)
+!
+! Revision 1.18  2006/02/23 09:55:33  raasch
+! nzb_2d replaced by nzb_s_inner/outer, nanz_2dh(0) replaced by
+! ngp_2dh_outer(:,0)
+!
+! Revision 1.17  2005/12/06 15:29:16  raasch
+! Mean pt profile is calculated only in case of the first intermediate timestep
+!
+! Revision 1.16  2004/01/30 10:14:31  raasch
+! Scalar lower k index nzb replaced by 2d-array nzb_2d
+!
+! Revision 1.15  2003/10/29 08:41:47  raasch
+! Horizontal mean temperature is now taken from array hom instead of array sums
+!
+! Revision 1.14  2003/03/16 09:27:22  raasch
+! Two underscores (_) are placed in front of all define-strings
+!
+! Revision 1.13  2003/03/12 16:20:47  raasch
+! Full code replaced in the call for all gridpoints instead of calling the
+! _ij version (required by NEC, because otherwise no vectorization)
+!
+! Revision 1.12  2002/12/19 13:47:38  raasch
+! STOP statement replaced by call of subroutine local_stop
+!
+! Revision 1.11  2002/06/11 12:33:07  raasch
+! Former subroutine changed to a module which allows to be called for all grid
+! points of a single vertical column with index i,j or for all grid points by
+! using function overloading.
+! Calculation of the horizontally averaged temperature profile moved to new
+! subroutine calc_mean_pt_profile, which is part of this module.
+!
+! Revision 1.10  2001/03/30 06:54:56  raasch
+! Translation of remaining German identifiers (variables, subroutines, etc.)
+!
+! Revision 1.9  2001/01/22 05:27:03  raasch
+! Module test_variables removed
+!
+! Revision 1.8  2000/07/03 12:55:56  raasch
+! Argument theta declared as pointer
+!
+! Revision 1.7  2000/04/27 07:04:19  raasch
+! in case of a sloping surface the temperature difference is computed relative
+! to the initial (2D) temperature field
+!
+! Revision 1.6  2000/04/13 14:27:34  schroeter
+! considering the influence of humidity to the buoyancy term
+! (only for sloping_surface = false),
+!
+! Revision 1.5  2000/01/21  16:17:23  16:17:23  letzel (Marcus Letzel)
+! All comments translated into English
+! 
+! Revision 1.4  1998/09/22 17:16:22  raasch
+! Auftriebsterm erweitert fuer Rechnungen mit geneigter Ebene
+!
+! Revision 1.3  1998/07/06 12:07:27  raasch
+! + USE test_variables
+!
+! Revision 1.2  1998/03/30 11:33:29  raasch
+! nanz_2dh in nanz_2dh(0) geaendert
 !
 ! Revision 1.1  1997/08/29 08:56:48  raasch
@@ -42,23 +77,10 @@
 ! Description:
 ! ------------
-!> Buoyancy term of the third component of the equation of motion.
-!> @attention Humidity is not regarded when using a sloping surface!
-!--------------------------------------------------------------------------------------------------!
- MODULE buoyancy_mod
-
-    USE basic_constants_and_equations_mod,                                                         &
-        ONLY:  g
-
-    USE arrays_3d,                                                                                 &
-        ONLY:  pt, pt_slope_ref, ref_state, tend
-
-    USE control_parameters,                                                                        &
-        ONLY:  atmos_ocean_sign, cos_alpha_surface, message_string, pt_surface, sin_alpha_surface, &
-               sloping_surface
-
-    USE kinds
+! Buoyancy term of the third component of the equation of motion.
+! WARNING: humidity is not regarded when using a sloping surface!
+!------------------------------------------------------------------------------!
 
     PRIVATE
-    PUBLIC buoyancy
+    PUBLIC buoyancy, calc_mean_pt_profile
 
     INTERFACE buoyancy
@@ -67,26 +89,26 @@
     END INTERFACE buoyancy
 
+    INTERFACE calc_mean_pt_profile
+       MODULE PROCEDURE calc_mean_pt_profile
+    END INTERFACE calc_mean_pt_profile
+
  CONTAINS
 
 
-!--------------------------------------------------------------------------------------------------!
-! Description:
-! ------------
-!> Call for all grid points
-!--------------------------------------------------------------------------------------------------!
-    SUBROUTINE buoyancy( var, wind_component )
-
-       USE indices,                                                                                &
-           ONLY:  nxl, nxlu, nxr, nyn, nys, nzb, nzt
-
+!------------------------------------------------------------------------------!
+! Call for all grid points
+!------------------------------------------------------------------------------!
+    SUBROUTINE buoyancy( theta, wind_component, pr )
+
+       USE arrays_3d
+       USE control_parameters
+       USE indices
+       USE pegrid
+       USE statistics
 
        IMPLICIT NONE
 
-       INTEGER(iwp) ::  i              !<
-       INTEGER(iwp) ::  j              !<
-       INTEGER(iwp) ::  k              !<
-       INTEGER(iwp) ::  wind_component !<
-
-       REAL(wp), DIMENSION(:,:,:), POINTER ::  var
+       INTEGER ::  i, j, k, pr, wind_component
+       REAL, DIMENSION(:,:,:), POINTER  ::  theta
 
 
@@ -94,16 +116,11 @@
 !
 !--       Normal case: horizontal surface
-          !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
-          !$ACC PRESENT(var) &
-          !$ACC PRESENT(ref_state) &
-          !$ACC PRESENT(tend)
           DO  i = nxl, nxr
              DO  j = nys, nyn
-                DO  k = nzb+1, nzt-1
-                   tend(k,j,i) = tend(k,j,i) + atmos_ocean_sign * g * 0.5_wp *                     &
-                                 (                                                                 &
-                                    ( var(k,j,i)   - ref_state(k) )   / ref_state(k) +             &
-                                    ( var(k+1,j,i) - ref_state(k+1) ) / ref_state(k+1)             &
-                                 )
+                DO  k = nzb_s_inner(j,i)+1, nzt-1
+                    tend(k,j,i) = tend(k,j,i) + g * 0.5 * (                    &
+                        ( theta(k,j,i)   - hom(k,1,pr,0)   ) / hom(k,1,pr,0) + &
+                        ( theta(k+1,j,i) - hom(k+1,1,pr,0) ) / hom(k+1,1,pr,0) &
+                                                          )
                 ENDDO
              ENDDO
@@ -112,18 +129,18 @@
        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.
+!--       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.
+!--       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  i = nxl, nxr
                 DO  j = nys, nyn
-                   DO  k = nzb+1, nzt-1
-                      tend(k,j,i) = tend(k,j,i) + g * sin_alpha_surface *                          &
-                                    0.5_wp * ( ( pt(k,j,i-1)         + pt(k,j,i)         )         &
-                                             - ( pt_slope_ref(k,i-1) + pt_slope_ref(k,i) )         &
-                                             ) / pt_surface
+                   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
@@ -134,9 +151,9 @@
              DO  i = nxl, nxr
                 DO  j = nys, nyn
-                   DO  k = nzb+1, nzt-1
-                      tend(k,j,i) = tend(k,j,i) + g * cos_alpha_surface *                          &
-                                    0.5_wp * ( ( pt(k,j,i)         + pt(k+1,j,i)         )         &
-                                             - ( pt_slope_ref(k,i) + pt_slope_ref(k+1,i) )         &
-                                             ) / pt_surface
+                   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
@@ -145,6 +162,7 @@
           ELSE
 
-             WRITE( message_string, * ) 'no term for component "', wind_component, '"'
-             CALL message( 'buoyancy', 'PA0159', 1, 2, 0, 6, 0 )
+             IF ( myid == 0 )  PRINT*, '+++ buoyancy: no term for component "',&
+                                       wind_component,'"'
+             CALL local_stop
 
           ENDIF
@@ -155,73 +173,135 @@
 
 
-!--------------------------------------------------------------------------------------------------!
+!------------------------------------------------------------------------------!
+! Call for grid point i,j
+!------------------------------------------------------------------------------!
+    SUBROUTINE buoyancy_ij( i, j, theta, 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, DIMENSION(:,:,:), POINTER  ::  theta
+
+
+       IF ( .NOT. sloping_surface )  THEN
+!
+!--       Normal case: horizontal surface
+          DO  k = nzb_s_inner(j,i)+1, nzt-1
+              tend(k,j,i) = tend(k,j,i) + g * 0.5 * (                          &
+                        ( theta(k,j,i)   - hom(k,1,pr,0)   ) / hom(k,1,pr,0) + &
+                        ( theta(k+1,j,i) - hom(k+1,1,pr,0) ) / hom(k+1,1,pr,0) &
+                                                    )
+          ENDDO
+
+       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  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
+
+          ELSEIF ( wind_component == 3 )  THEN
+
+             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
+
+          ELSE
+
+             IF ( myid == 0 )  PRINT*, '+++ buoyancy: no term for component "',&
+                                       wind_component,'"'
+             CALL local_stop
+
+          ENDIF
+
+       ENDIF
+
+    END SUBROUTINE buoyancy_ij
+
+
+    SUBROUTINE calc_mean_pt_profile( theta, pr )
+
+!------------------------------------------------------------------------------!
 ! Description:
 ! ------------
-!> Call for grid point i,j
-!> @attention PGI-compiler creates SIGFPE if opt>1 is used! Therefore, opt=1 is forced by
-!>            compiler-directive.
-!--------------------------------------------------------------------------------------------------!
-!pgi$r opt=1
-    SUBROUTINE buoyancy_ij( i, j, var, wind_component )
-
-
-       USE indices,                                                                                &
-           ONLY:  nzb, nzt
+! Calculate the horizontally averaged vertical temperature profile (pr=4 in case
+! of potential temperature and 44 in case of virtual potential temperature).
+!------------------------------------------------------------------------------!
+
+       USE control_parameters
+       USE indices
+       USE pegrid
+       USE statistics
 
        IMPLICIT NONE
 
-       INTEGER(iwp) ::  i              !<
-       INTEGER(iwp) ::  j              !<
-       INTEGER(iwp) ::  k              !<
-       INTEGER(iwp) ::  wind_component !<
-
-       REAL(wp), DIMENSION(:,:,:), POINTER ::  var
-
-
-       IF ( .NOT. sloping_surface )  THEN
-!
-!--       Normal case: horizontal surface
-          DO  k = nzb+1, nzt-1
-              tend(k,j,i) = tend(k,j,i) + atmos_ocean_sign * g * 0.5_wp * (                        &
-                            ( var(k,j,i)   - ref_state(k)   ) / ref_state(k)   +                   &
-                            ( var(k+1,j,i) - ref_state(k+1) ) / ref_state(k+1)                     &
-                                                                          )
+       INTEGER ::  i, j, k, omp_get_thread_num, pr, tn
+       REAL, DIMENSION(:,:,:), POINTER  ::  theta
+
+!
+!--    Computation of the horizontally averaged temperature profile, unless
+!--    already done by the relevant call from flow_statistics. The calculation
+!--    is done only for the first respective intermediate timestep in order to
+!--    spare communication time and to produce identical model results with jobs
+!--    which are calling flow_statistics at different time intervals.
+!--    Although this calculation is not required for model runs with a slope,
+!--    it is nevertheless also computed.
+       IF ( .NOT. flow_statistics_called  .AND. &
+            intermediate_timestep_count == 1 )  THEN
+
+!
+!--       Horizontal average of the potential temperature
+          tn           =   0  ! Default thread number in case of one thread
+          !$OMP PARALLEL PRIVATE( i, j, k, tn )
+!$        tn = omp_get_thread_num()
+          sums_l(:,pr,tn) = 0.0
+          !$OMP DO
+          DO  i = nxl, nxr
+             DO  j =  nys, nyn
+                DO  k = nzb_s_outer(j,i), nzt+1
+                   sums_l(k,pr,tn) = sums_l(k,pr,tn) + theta(k,j,i)
+                ENDDO
+             ENDDO
           ENDDO
-
-       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  k = nzb+1, nzt-1
-                tend(k,j,i) = tend(k,j,i) + g * sin_alpha_surface *                                &
-                              0.5_wp * ( ( pt(k,j,i-1)         + pt(k,j,i)         )               &
-                                       - ( pt_slope_ref(k,i-1) + pt_slope_ref(k,i) )               &
-                                       ) / pt_surface
-             ENDDO
-
-          ELSEIF ( wind_component == 3 )  THEN
-
-             DO  k = nzb+1, nzt-1
-                tend(k,j,i) = tend(k,j,i) + g * cos_alpha_surface *                                &
-                              0.5_wp * ( ( pt(k,j,i)         + pt(k+1,j,i)         )               &
-                                       - ( pt_slope_ref(k,i) + pt_slope_ref(k+1,i) )               &
-                                       ) / pt_surface
-             ENDDO
-
-          ELSE
-
-             WRITE( message_string, * ) 'no term for component "', wind_component, '"'
-             CALL message( 'buoyancy', 'PA0159', 1, 2, 0, 6, 0 )
-
-          ENDIF
+          !$OMP END PARALLEL
+
+          DO  i = 1, threads_per_task-1
+             sums_l(:,pr,0) = sums_l(:,pr,0) + sums_l(:,pr,i)
+          ENDDO
+
+#if defined( __parallel )
+
+          CALL MPI_ALLREDUCE( sums_l(nzb,pr,0), sums(nzb,pr), nzt+2-nzb, &
+                              MPI_REAL, MPI_SUM, comm2d, ierr )
+
+#else
+
+          sums(:,pr) = sums_l(:,pr,0)
+
+#endif
+
+          hom(:,1,pr,0) = sums(:,pr) / ngp_2dh_outer(:,0)
 
        ENDIF
 
-    END SUBROUTINE buoyancy_ij
+    END SUBROUTINE calc_mean_pt_profile
 
  END MODULE buoyancy_mod
