MODULE buoyancy_mod

!------------------------------------------------------------------------------!
! Actual revisions:
! -----------------
! 
!
! Former revisions:
! -----------------
! $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
! Initial revision
!
!
! Description:
! ------------
! Buoyancy term of the third component of the equation of motion.
! WARNING: humidity is not regarded when using a sloping surface!
!------------------------------------------------------------------------------!

    PRIVATE
    PUBLIC buoyancy, calc_mean_pt_profile

    INTERFACE buoyancy
       MODULE PROCEDURE buoyancy
       MODULE PROCEDURE buoyancy_ij
    END INTERFACE buoyancy

    INTERFACE calc_mean_pt_profile
       MODULE PROCEDURE calc_mean_pt_profile
    END INTERFACE calc_mean_pt_profile

 CONTAINS


!------------------------------------------------------------------------------!
! 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 ::  i, j, k, pr, wind_component
       REAL, DIMENSION(:,:,:), POINTER  ::  theta


       IF ( .NOT. sloping_surface )  THEN
!
!--       Normal case: horizontal surface
          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 * 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
          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  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 * 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

             IF ( myid == 0 )  PRINT*, '+++ buoyancy: no term for component "',&
                                       wind_component,'"'
             CALL local_stop

          ENDIF

       ENDIF

    END SUBROUTINE buoyancy


!------------------------------------------------------------------------------!
! 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:
! ------------
! 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 ::  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
          !$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 calc_mean_pt_profile

 END MODULE buoyancy_mod