source: palm/trunk/SOURCE/buoyancy.f90 @ 1097

 MODULE buoyancy_mod

!--------------------------------------------------------------------------------!
! This file is part of PALM.
!
! 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-2012  Leibniz University Hannover
!--------------------------------------------------------------------------------!
!
! Currrent revisions:
! -----------------
!
!
! Former revisions:
! -----------------
! $Id: buoyancy.f90 1037 2012-10-22 14:10:22Z raasch $
!
! 1036 2012-10-22 13:43:42Z raasch
! code put under GPL (PALM 3.9)
!
! 1015 2012-09-27 09:23:24Z raasch
! accelerator version (*_acc) added
!
! 1010 2012-09-20 07:59:54Z raasch
! cpp switch __nopointer added for pointer free version
!
! 622 2010-12-10 08:08:13Z raasch
! optional barriers included in order to speed up collective operations
!
! 515 2010-03-18 02:30:38Z raasch
! PGI-compiler creates SIGFPE in routine buoyancy, if opt>1 is used! Therefore,
! opt=1 is forced by compiler-directive.
!
! 247 2009-02-27 14:01:30Z heinze
! Output of messages replaced by message handling routine
!
! 132 2007-11-20 09:46:11Z letzel
! Vertical scalar profiles now based on nzb_s_inner and ngp_2dh_s_inner.
!
! 106 2007-08-16 14:30:26Z raasch
! i loop for u-component (sloping surface) is starting from nxlu (needed for
! non-cyclic boundary conditions)
! 
! 97 2007-06-21 08:23:15Z raasch
! Routine reneralized to be used with temperature AND density:
! argument theta renamed var, new argument var_reference,
! use_pt_reference renamed use_reference,
! calc_mean_pt_profile renamed calc_mean_profile
!
! 57 2007-03-09 12:05:41Z raasch
! Reference temperature pt_reference can be used.
!
! RCS Log replace by Id keyword, revision history cleaned up
!
! Revision 1.19  2006/04/26 12:09:56  raasch
! OpenMP optimization (one dimension added to sums_l)
!
! 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, buoyancy_acc, calc_mean_profile

    INTERFACE buoyancy
       MODULE PROCEDURE buoyancy
       MODULE PROCEDURE buoyancy_ij
    END INTERFACE buoyancy

    INTERFACE buoyancy_acc
       MODULE PROCEDURE buoyancy_acc
    END INTERFACE buoyancy_acc

    INTERFACE calc_mean_profile
       MODULE PROCEDURE calc_mean_profile
    END INTERFACE calc_mean_profile

 CONTAINS


!------------------------------------------------------------------------------!
! Call for all grid points
!------------------------------------------------------------------------------!
    SUBROUTINE buoyancy( 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
             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)   ) / hom(k,1,pr,0) + &
                          ( var(k+1,j,i) - hom(k+1,1,pr,0) ) / hom(k+1,1,pr,0) &
                                                            )
                   ENDDO
                ENDDO
             ENDDO
          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


!------------------------------------------------------------------------------!
! 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


!------------------------------------------------------------------------------!
! 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, 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  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
          ELSE
             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)   ) / hom(k,1,pr,0) + &
                          ( var(k+1,j,i) - hom(k+1,1,pr,0) ) / hom(k+1,1,pr,0) &
                                                                          )
             ENDDO
          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  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

             WRITE( message_string, * ) 'no term for component "',&
                                       wind_component,'"'
             CALL message( 'buoyancy', 'PA0159', 1, 2, 0, 6, 0 )

          ENDIF

       ENDIF

    END SUBROUTINE buoyancy_ij


    SUBROUTINE calc_mean_profile( var, 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
#if defined( __nopointer )
       REAL, DIMENSION(:,:,:) ::  var
#else
       REAL, DIMENSION(:,:,:), POINTER ::  var
#endif

!
!--    Computation of the horizontally averaged profile of variable var, 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.
       IF ( .NOT. flow_statistics_called  .AND. &
            intermediate_timestep_count == 1 )  THEN

!
!--       Horizontal average of variable var
          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_inner(j,i), nzt+1
                   sums_l(k,pr,tn) = sums_l(k,pr,tn) + var(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 )

          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
          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_s_inner(:,0)

       ENDIF

    END SUBROUTINE calc_mean_profile

 END MODULE buoyancy_mod