source: palm/trunk/SOURCE/subsidence.f90 @ 1217

 MODULE subsidence_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
!--------------------------------------------------------------------------------!
!
! Current revisions:
! -----------------
!
! Former revisions:
! -----------------
! $Id: subsidence.f90 1037 2012-10-22 14:10:22Z raasch $
!
! 1036 2012-10-22 13:43:42Z raasch
! code put under GPL (PALM 3.9)
!
! 671 2011-01-11 12:04:00Z heinze $
! bugfix: access to ddzu(nzt+2) which is not defined
!
! 667 2010-12-23 12:06:00Z suehring
! nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng
!
! 580 2010-10-05 13:59:11Z heinze
! Renaming of ws_vertical_gradient to subs_vertical_gradient,
! ws_vertical_gradient_level to subs_vertical_gradient_level and
! ws_vertical_gradient_level_ind to subs_vertical_gradient_level_i
!
! Revision 3.7 2009-12-11 14:15:58Z heinze
! Initial revision 
!
! Description:
! ------------
! Impact of large-scale subsidence or ascent as tendency term for use 
! in the prognostic equation of potential temperature. This enables the 
! construction of a constant boundary layer height z_i with time.
!-----------------------------------------------------------------------------!


    IMPLICIT NONE

    PRIVATE
    PUBLIC  init_w_subsidence, subsidence

    INTERFACE init_w_subsidence
       MODULE PROCEDURE init_w_subsidence
    END INTERFACE init_w_subsidence

    INTERFACE subsidence
       MODULE PROCEDURE subsidence
       MODULE PROCEDURE subsidence_ij
    END INTERFACE subsidence

 CONTAINS

    SUBROUTINE init_w_subsidence 

       USE arrays_3d
       USE control_parameters
       USE grid_variables
       USE indices
       USE pegrid
       USE statistics 

       IMPLICIT NONE

       INTEGER :: i, k
       REAL    :: gradient, ws_surface

       IF ( .NOT. ALLOCATED( w_subs )) THEN
          ALLOCATE( w_subs(nzb:nzt+1) )
          w_subs = 0.0
       ENDIF 

      IF ( ocean )  THEN
          message_string = 'Applying large scale vertical motion is not ' // &
                           'allowed for ocean runs'
          CALL message( 'init_w_subsidence', 'PA0324', 2, 2, 0, 6, 0 )
       ENDIF

!
!--   Compute the profile of the subsidence/ascent velocity 
!--   using the given gradients
      i = 1
      gradient = 0.0
      ws_surface = 0.0
      

      subs_vertical_gradient_level_i(1) = 0
      DO  k = 1, nzt+1
         IF ( i < 11 ) THEN
            IF ( subs_vertical_gradient_level(i) < zu(k)  .AND. &
                 subs_vertical_gradient_level(i) >= 0.0 )  THEN
               gradient = subs_vertical_gradient(i) / 100.0
               subs_vertical_gradient_level_i(i) = k - 1
               i = i + 1
            ENDIF
         ENDIF
         IF ( gradient /= 0.0 )  THEN
            IF ( k /= 1 )  THEN
               w_subs(k) = w_subs(k-1) + dzu(k) * gradient
            ELSE
               w_subs(k) = ws_surface   + 0.5 * dzu(k) * gradient
            ENDIF
         ELSE
            w_subs(k) = w_subs(k-1)
         ENDIF
      ENDDO

!
!--   In case of no given gradients for the subsidence/ascent velocity,
!--   choose zero gradient
      IF ( subs_vertical_gradient_level(1) == -9999999.9 )  THEN
         subs_vertical_gradient_level(1) = 0.0
      ENDIF

    END SUBROUTINE init_w_subsidence


    SUBROUTINE subsidence( tendency, var, var_init ) 

       USE arrays_3d
       USE control_parameters
       USE grid_variables
       USE indices
       USE pegrid
       USE statistics 

       IMPLICIT NONE
 
       INTEGER :: i, j, k

       REAL :: tmp_grad
    
       REAL, DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: var, tendency
       REAL, DIMENSION(nzb:nzt+1) :: var_init, var_mod

       var_mod = var_init

!
!--    Influence of w_subsidence on the current tendency term
       DO  i = nxl, nxr
          DO  j = nys, nyn
             DO  k = nzb_s_inner(j,i)+1, nzt 
                IF ( w_subs(k) < 0.0 ) THEN    ! large-scale subsidence
                   tendency(k,j,i) = tendency(k,j,i) - w_subs(k) *  &
                                     ( var(k+1,j,i) - var(k,j,i) ) * ddzu(k+1)
                ELSE                           ! large-scale ascent
                   tendency(k,j,i) = tendency(k,j,i) - w_subs(k) *  &
                                     ( var(k,j,i) - var(k-1,j,i) ) * ddzu(k)
                ENDIF
             ENDDO
          ENDDO
       ENDDO

!
!--    Shifting of the initial profile is especially necessary with Rayleigh 
!--    damping switched on

       DO  k = nzb, nzt
          IF ( w_subs(k) < 0.0 ) THEN      ! large-scale subsidence
             var_mod(k) = var_init(k) - dt_3d * w_subs(k) *  &
                               ( var_init(k+1) - var_init(k) ) * ddzu(k+1)
          ENDIF
       ENDDO
!
!--   At the upper boundary, the initial profile is shifted with aid of 
!--   the gradient tmp_grad. (This is ok if the gradients are linear.)
      IF ( w_subs(nzt) < 0.0 ) THEN
         tmp_grad = ( var_init(nzt+1) - var_init(nzt) ) * ddzu(nzt+1)
         var_mod(nzt+1) = var_init(nzt+1) -  &
                              dt_3d * w_subs(nzt+1) * tmp_grad
      ENDIF
       

      DO  k = nzt+1, nzb+1, -1
         IF ( w_subs(k) >= 0.0 ) THEN  ! large-scale ascent
            var_mod(k) = var_init(k) - dt_3d * w_subs(k) *  &
                               ( var_init(k) - var_init(k-1) ) * ddzu(k) 
         ENDIF
      ENDDO
!
!--   At the lower boundary shifting is not necessary because the 
!--   subsidence velocity w_subs(nzb) vanishes.


      IF ( w_subs(nzb+1) >= 0.0 ) THEN
         var_mod(nzb) = var_init(nzb)
      ENDIF

      var_init = var_mod


 END SUBROUTINE subsidence

 SUBROUTINE subsidence_ij( i, j, tendency, var, var_init ) 

       USE arrays_3d
       USE control_parameters
       USE grid_variables
       USE indices
       USE pegrid
       USE statistics 

       IMPLICIT NONE
 
       INTEGER :: i, j, k

       REAL :: tmp_grad
    
       REAL, DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: var, tendency
       REAL, DIMENSION(nzb:nzt+1) :: var_init, var_mod

       var_mod = var_init

!
!--    Influence of w_subsidence on the current tendency term
       DO  k = nzb_s_inner(j,i)+1, nzt 
          IF ( w_subs(k) < 0.0 ) THEN      ! large-scale subsidence
             tendency(k,j,i) = tendency(k,j,i) - w_subs(k) *  &
                               ( var(k+1,j,i) - var(k,j,i) ) * ddzu(k+1)
          ELSE                             ! large-scale ascent
             tendency(k,j,i) = tendency(k,j,i) - w_subs(k) *  &
                               ( var(k,j,i) - var(k-1,j,i) ) * ddzu(k)
          ENDIF
       ENDDO


!
!--    Shifting of the initial profile is especially necessary with Rayleigh 
!--    damping switched on
       IF ( i == nxl .AND. j == nys ) THEN ! shifting only once per PE 

          DO  k = nzb, nzt
             IF ( w_subs(k) < 0.0 ) THEN      ! large-scale subsidence
                var_mod(k) = var_init(k) - dt_3d * w_subs(k) *  &
                                  ( var_init(k+1) - var_init(k) ) * ddzu(k+1)
             ENDIF
          ENDDO
!
!--       At the upper boundary, the initial profile is shifted with aid of 
!--       the gradient tmp_grad. (This is ok if the gradients are linear.)
          IF ( w_subs(nzt) < 0.0 ) THEN
             tmp_grad = ( var_init(nzt+1) - var_init(nzt) ) * ddzu(nzt+1)
             var_mod(nzt+1) = var_init(nzt+1) -  &
                                  dt_3d * w_subs(nzt+1) * tmp_grad
          ENDIF
       

          DO  k = nzt+1, nzb+1, -1
             IF ( w_subs(k) >= 0.0 ) THEN  ! large-scale ascent
                var_mod(k) = var_init(k) - dt_3d * w_subs(k) *  &
                                   ( var_init(k) - var_init(k-1) ) * ddzu(k)
             ENDIF
          ENDDO
!
!--       At the lower boundary shifting is not necessary because the 
!--       subsidence velocity w_subs(nzb) vanishes.


          IF ( w_subs(nzb+1) >= 0.0 ) THEN
             var_mod(nzb) = var_init(nzb)
          ENDIF

          var_init = var_mod 

       ENDIF

 END SUBROUTINE subsidence_ij


 END MODULE subsidence_mod