source: palm/trunk/SOURCE/disturb_field.f90 @ 1407

 SUBROUTINE disturb_field( nzb_uv_inner, dist1, field )

!--------------------------------------------------------------------------------!
! 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-2014 Leibniz Universitaet Hannover
!--------------------------------------------------------------------------------!
!
! Current revisions:
! ------------------
! 
! 
! Former revisions:
! -----------------
! $Id: disturb_field.f90 1401 2014-05-09 14:15:46Z raasch $
!
! 1400 2014-05-09 14:03:54Z knoop
! Parallel random number generator added
! 
! 1353 2014-04-08 15:21:23Z heinze
! REAL constants provided with KIND-attribute 
!
! 1322 2014-03-20 16:38:49Z raasch
! REAL constants defined as wp-kind
!
! 1320 2014-03-20 08:40:49Z raasch
! ONLY-attribute added to USE-statements,
! kind-parameters added to all INTEGER and REAL declaration statements, 
! kinds are defined in new module kinds, 
! revision history before 2012 removed,
! comment fields (!:) to be used for variable explanations added to
! all variable declaration statements 
!
! 1036 2012-10-22 13:43:42Z raasch
! code put under GPL (PALM 3.9)
!
! Revision 1.1  1998/02/04 15:40:45  raasch
! Initial revision
!
!
! Description:
! ------------
! Imposing a random perturbation on a 3D-array.
! On parallel computers, the random number generator is as well called for all
! gridpoints of the total domain to ensure, regardless of the number of PEs
! used, that the elements of the array have the same values in the same
! order in every case. The perturbation range is steered by dist_range.
!------------------------------------------------------------------------------!

    USE control_parameters,   &
        ONLY:  dist_nxl, dist_nxr, dist_nyn, dist_nys, dist_range,             &
               disturbance_amplitude, disturbance_created,                     &
               disturbance_level_ind_b, disturbance_level_ind_t, iran,         &
               random_generator, topography
                
    USE cpulog,                                                                &
        ONLY:  cpu_log, log_point
        
    USE indices,                                                               &
        ONLY:  nbgp, nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, nzt
        
    USE kinds
    
    USE random_function_mod,                                                   &
        ONLY: random_function
        
    USE random_generator_parallel,                                             &
        ONLY:  random_number_parallel, random_seed_parallel, random_dummy,     &
               id_random_array, seq_random_array

    IMPLICIT NONE

    INTEGER(iwp) ::  i  !:
    INTEGER(iwp) ::  j  !:
    INTEGER(iwp) ::  k  !:
    
    INTEGER(iwp) ::  nzb_uv_inner(nysg:nyng,nxlg:nxrg) !:

    REAL(wp) ::  randomnumber  !:
    
    REAL(wp) ::  dist1(nzb:nzt+1,nysg:nyng,nxlg:nxrg)  !:
    REAL(wp) ::  field(nzb:nzt+1,nysg:nyng,nxlg:nxrg)  !:
    
    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  dist2  !:


    CALL cpu_log( log_point(20), 'disturb_field', 'start' )

!
!-- Create an additional temporary array and initialize the arrays needed
!-- to store the disturbance
    ALLOCATE( dist2(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    dist1 = 0.0_wp
    dist2 = 0.0_wp

!
!-- Create the random perturbation and store it on temporary array
    IF ( random_generator == 'numerical-recipes' )  THEN
       DO  i = dist_nxl(dist_range), dist_nxr(dist_range)
          DO  j = dist_nys(dist_range), dist_nyn(dist_range)
             DO  k = disturbance_level_ind_b, disturbance_level_ind_t
                randomnumber = 3.0_wp * disturbance_amplitude *                &
                               ( random_function( iran ) - 0.5_wp )
                IF ( nxl <= i  .AND.  nxr >= i  .AND.  nys <= j  .AND.         &
                     nyn >= j )                                                &
                THEN
                   dist1(k,j,i) = randomnumber
                ENDIF
             ENDDO
          ENDDO
       ENDDO
    ELSEIF ( random_generator == 'random-parallel' )  THEN
       DO  i = dist_nxl(dist_range), dist_nxr(dist_range)
          DO  j = dist_nys(dist_range), dist_nyn(dist_range)
             CALL random_seed_parallel( put=seq_random_array(:, j, i) )
             DO  k = disturbance_level_ind_b, disturbance_level_ind_t
                CALL random_number_parallel( random_dummy )
                randomnumber = 3.0_wp * disturbance_amplitude *                &
                               ( random_dummy - 0.5_wp )
                IF ( nxl <= i  .AND.  nxr >= i  .AND.  nys <= j  .AND.         &
                     nyn >= j )                                                &
                THEN
                   dist1(k,j,i) = randomnumber
                ENDIF
             ENDDO
!              IF (myid == 0) PRINT*, nxl, i, nxr, i, nys, j, nyn, j
             CALL random_seed_parallel( get=seq_random_array(:, j, i) )
          ENDDO
       ENDDO
    ELSEIF ( random_generator == 'system-specific' )  THEN
       DO  i = dist_nxl(dist_range), dist_nxr(dist_range)
          DO  j = dist_nys(dist_range), dist_nyn(dist_range)
             DO  k = disturbance_level_ind_b, disturbance_level_ind_t
#if defined( __nec )
                randomnumber = 3.0_wp * disturbance_amplitude *                &
                               ( RANDOM( 0 ) - 0.5_wp )
#else
                CALL RANDOM_NUMBER( randomnumber )
                randomnumber = 3.0_wp * disturbance_amplitude *                &
                                ( randomnumber - 0.5_wp )
#endif
                IF ( nxl <= i .AND. nxr >= i .AND. nys <= j .AND. nyn >= j )   &
                THEN
                   dist1(k,j,i) = randomnumber
                ENDIF
             ENDDO
          ENDDO
       ENDDO

    ENDIF

!
!-- Exchange of ghost points for the random perturbation

    CALL exchange_horiz( dist1, nbgp )
!
!-- Applying the Shuman filter in order to smooth the perturbations.
!-- Neighboured grid points in all three directions are used for the
!-- filter operation.
!-- Loop has been splitted to make runs reproducible on HLRN systems using
!-- compiler option -O3
     DO  i = nxl, nxr
        DO  j = nys, nyn
          DO  k = disturbance_level_ind_b-1, disturbance_level_ind_t+1
             dist2(k,j,i) = ( dist1(k,j,i-1) + dist1(k,j,i+1)                  &
                            + dist1(k,j+1,i) + dist1(k+1,j,i)                  &
                            ) / 12.0_wp
          ENDDO
          DO  k = disturbance_level_ind_b-1, disturbance_level_ind_t+1
              dist2(k,j,i) = dist2(k,j,i) + ( dist1(k,j-1,i) + dist1(k-1,j,i)  &
                            + 6.0_wp * dist1(k,j,i)                            &
                            ) / 12.0_wp
          ENDDO
        ENDDO
     ENDDO

!
!-- Exchange of ghost points for the filtered perturbation.
!-- Afterwards, filter operation and exchange of ghost points are repeated.
    CALL exchange_horiz( dist2, nbgp )

    DO  i = nxl, nxr
       DO  j = nys, nyn
          DO  k = disturbance_level_ind_b-2, disturbance_level_ind_t+2
             dist1(k,j,i) = ( dist2(k,j,i-1) + dist2(k,j,i+1) + dist2(k,j-1,i) &
                            + dist2(k,j+1,i) + dist2(k+1,j,i) + dist2(k-1,j,i) &
                            + 6.0_wp * dist2(k,j,i)                            &
                            ) / 12.0_wp
          ENDDO
       ENDDO
    ENDDO

    CALL exchange_horiz( dist1, nbgp )

!
!-- Remove perturbations below topography (including one gridpoint above it
!-- in order to allow for larger timesteps at the beginning of the simulation
!-- (diffusion criterion))
    IF ( TRIM( topography ) /= 'flat' )  THEN
       DO  i = nxlg, nxrg
          DO  j = nysg, nyng
             dist1(nzb:nzb_uv_inner(j,i)+1,j,i) = 0.0_wp
          ENDDO
       ENDDO
    ENDIF

!
!-- Random perturbation is added to the array to be disturbed.
    DO  i = nxlg, nxrg
       DO  j = nysg, nyng
          DO  k = disturbance_level_ind_b-2, disturbance_level_ind_t+2
             field(k,j,i) = field(k,j,i) + dist1(k,j,i)
          ENDDO
       ENDDO
    ENDDO

!
!-- Deallocate the temporary array
    DEALLOCATE( dist2 )

!
!-- Set a flag, which indicates that a random perturbation is imposed
    disturbance_created = .TRUE.


    CALL cpu_log( log_point(20), 'disturb_field', 'stop' )


 END SUBROUTINE disturb_field