source: palm/trunk/SOURCE/microphysics.f90 @ 1310

 MODULE microphysics_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-2014 Leibniz Universitaet Hannover
!--------------------------------------------------------------------------------!
!
! Current revisions:
! ------------------
!
!
! Former revisions:
! -----------------
! $Id: microphysics.f90 1310 2014-03-14 08:01:56Z raasch $
!
! 1241 2013-10-30 11:36:58Z heinze
! hyp and rho have to be calculated at each time step if data from external
! file LSF_DATA are used
!
! 1115 2013-03-26 18:16:16Z hoffmann
! microphyical tendencies are calculated in microphysics_control in an optimized
! way; unrealistic values are prevented; bugfix in evaporation; some reformatting
!
! 1106 2013-03-04 05:31:38Z raasch
! small changes in code formatting
!
! 1092 2013-02-02 11:24:22Z raasch
! unused variables removed
! file put under GPL
!
! 1065 2012-11-22 17:42:36Z hoffmann
! Sedimentation process implemented according to Stevens and Seifert (2008).
! Turbulence effects on autoconversion and accretion added (Seifert, Nuijens
! and Stevens, 2010).
!
! 1053 2012-11-13 17:11:03Z hoffmann
! initial revision
! 
! Description:
! ------------
! Calculate cloud microphysics according to the two moment bulk 
! scheme by Seifert and Beheng (2006).
!------------------------------------------------------------------------------!

    PRIVATE
    PUBLIC microphysics_control

    INTERFACE microphysics_control
       MODULE PROCEDURE microphysics_control
       MODULE PROCEDURE microphysics_control_ij
    END INTERFACE microphysics_control

    INTERFACE adjust_cloud
       MODULE PROCEDURE adjust_cloud
       MODULE PROCEDURE adjust_cloud_ij
    END INTERFACE adjust_cloud

    INTERFACE autoconversion
       MODULE PROCEDURE autoconversion
       MODULE PROCEDURE autoconversion_ij
    END INTERFACE autoconversion

    INTERFACE accretion
       MODULE PROCEDURE accretion
       MODULE PROCEDURE accretion_ij
    END INTERFACE accretion

    INTERFACE selfcollection_breakup
       MODULE PROCEDURE selfcollection_breakup
       MODULE PROCEDURE selfcollection_breakup_ij
    END INTERFACE selfcollection_breakup

    INTERFACE evaporation_rain
       MODULE PROCEDURE evaporation_rain
       MODULE PROCEDURE evaporation_rain_ij
    END INTERFACE evaporation_rain

    INTERFACE sedimentation_cloud
       MODULE PROCEDURE sedimentation_cloud
       MODULE PROCEDURE sedimentation_cloud_ij
    END INTERFACE sedimentation_cloud
 
    INTERFACE sedimentation_rain
       MODULE PROCEDURE sedimentation_rain
       MODULE PROCEDURE sedimentation_rain_ij
    END INTERFACE sedimentation_rain

 CONTAINS


!------------------------------------------------------------------------------!
! Call for all grid points
!------------------------------------------------------------------------------!
    SUBROUTINE microphysics_control

       USE arrays_3d
       USE cloud_parameters
       USE control_parameters
       USE grid_variables
       USE indices
       USE statistics

       IMPLICIT NONE

       INTEGER ::  i, j, k

 
       DO  i = nxl, nxr
          DO  j = nys, nyn
             DO  k = nzb_s_inner(j,i)+1, nzt

             ENDDO
          ENDDO
       ENDDO

    END SUBROUTINE microphysics_control

    SUBROUTINE adjust_cloud

       USE arrays_3d
       USE cloud_parameters
       USE indices

       IMPLICIT NONE

       INTEGER ::  i, j, k

 
       DO  i = nxl, nxr
          DO  j = nys, nyn
             DO  k = nzb_s_inner(j,i)+1, nzt

             ENDDO
          ENDDO
       ENDDO

    END SUBROUTINE adjust_cloud


    SUBROUTINE autoconversion

       USE arrays_3d
       USE cloud_parameters
       USE control_parameters
       USE grid_variables
       USE indices

       IMPLICIT NONE

       INTEGER ::  i, j, k

 
       DO  i = nxl, nxr
          DO  j = nys, nyn
             DO  k = nzb_s_inner(j,i)+1, nzt

             ENDDO
          ENDDO
       ENDDO

    END SUBROUTINE autoconversion


    SUBROUTINE accretion

       USE arrays_3d
       USE cloud_parameters
       USE control_parameters
       USE indices

       IMPLICIT NONE

       INTEGER ::  i, j, k

 
       DO  i = nxl, nxr
          DO  j = nys, nyn
             DO  k = nzb_s_inner(j,i)+1, nzt

             ENDDO
          ENDDO
       ENDDO

    END SUBROUTINE accretion


    SUBROUTINE selfcollection_breakup

       USE arrays_3d
       USE cloud_parameters
       USE control_parameters
       USE indices

       IMPLICIT NONE

       INTEGER ::  i, j, k

 
       DO  i = nxl, nxr
          DO  j = nys, nyn
             DO  k = nzb_s_inner(j,i)+1, nzt

             ENDDO
          ENDDO
       ENDDO

    END SUBROUTINE selfcollection_breakup


    SUBROUTINE evaporation_rain

       USE arrays_3d
       USE cloud_parameters
       USE constants
       USE control_parameters
       USE indices

       IMPLICIT NONE

       INTEGER ::  i, j, k

 
       DO  i = nxl, nxr
          DO  j = nys, nyn
             DO  k = nzb_s_inner(j,i)+1, nzt

             ENDDO
          ENDDO
       ENDDO

    END SUBROUTINE evaporation_rain


    SUBROUTINE sedimentation_cloud

       USE arrays_3d
       USE cloud_parameters
       USE constants
       USE control_parameters
       USE indices

       IMPLICIT NONE

       INTEGER ::  i, j, k

 
       DO  i = nxl, nxr
          DO  j = nys, nyn
             DO  k = nzb_s_inner(j,i)+1, nzt

             ENDDO
          ENDDO
       ENDDO

    END SUBROUTINE sedimentation_cloud


    SUBROUTINE sedimentation_rain

       USE arrays_3d
       USE cloud_parameters
       USE constants
       USE control_parameters
       USE indices
       USE statistics

       IMPLICIT NONE

       INTEGER ::  i, j, k

 
       DO  i = nxl, nxr
          DO  j = nys, nyn
             DO  k = nzb_s_inner(j,i)+1, nzt

             ENDDO
          ENDDO
       ENDDO

    END SUBROUTINE sedimentation_rain


!------------------------------------------------------------------------------!
! Call for grid point i,j
!------------------------------------------------------------------------------!

    SUBROUTINE microphysics_control_ij( i, j )

       USE arrays_3d
       USE cloud_parameters
       USE control_parameters
       USE grid_variables
       USE indices
       USE statistics

       IMPLICIT NONE

       INTEGER ::  i, j, k 
       REAL    ::  t_surface

       IF ( large_scale_forcing .AND. lsf_surf ) THEN
!
!--       Calculate:
!--       pt / t : ratio of potential and actual temperature (pt_d_t)
!--       t / pt : ratio of actual and potential temperature (t_d_pt)
!--       p_0(z) : vertical profile of the hydrostatic pressure (hyp)
          t_surface = pt_surface * ( surface_pressure / 1000.0 )**0.286
          DO  k = nzb, nzt+1
             hyp(k)    = surface_pressure * 100.0 * &
                         ( (t_surface - g/cp * zu(k)) / t_surface )**(1.0/0.286)
             pt_d_t(k) = ( 100000.0 / hyp(k) )**0.286
             t_d_pt(k) = 1.0 / pt_d_t(k)
             hyrho(k)  = hyp(k) / ( r_d * t_d_pt(k) * pt_init(k) )       
          ENDDO
!
!--       Compute reference density
          rho_surface = surface_pressure * 100.0 / ( r_d * t_surface )
       ENDIF


       dt_micro = dt_3d * weight_pres(intermediate_timestep_count)
!
!--    Adjust unrealistic values
       IF ( precipitation )  CALL adjust_cloud( i,j ) 
!
!--    Use 1-d arrays
       q_1d(:)  = q(:,j,i)
       pt_1d(:) = pt(:,j,i)
       qc_1d(:) = qc(:,j,i)
       nc_1d(:) = nc_const
       IF ( precipitation )  THEN
          qr_1d(:) = qr(:,j,i)
          nr_1d(:) = nr(:,j,i)
       ENDIF
!
!--    Compute cloud physics
       IF ( precipitation )  THEN
          CALL autoconversion( i,j )
          CALL accretion( i,j )
          CALL selfcollection_breakup( i,j )
          CALL evaporation_rain( i,j )
          CALL sedimentation_rain( i,j )
       ENDIF

       IF ( drizzle )  CALL sedimentation_cloud( i,j )
!
!--    Derive tendencies
       tend_q(:,j,i)  = ( q_1d(:) - q(:,j,i) ) / dt_micro
       tend_pt(:,j,i) = ( pt_1d(:) - pt(:,j,i) ) / dt_micro
       IF ( precipitation )  THEN
          tend_qr(:,j,i) = ( qr_1d(:) - qr(:,j,i) ) / dt_micro
          tend_nr(:,j,i) = ( nr_1d(:) - nr(:,j,i) ) / dt_micro
       ENDIF

    END SUBROUTINE microphysics_control_ij

    SUBROUTINE adjust_cloud_ij( i, j )

       USE arrays_3d
       USE cloud_parameters
       USE indices

       IMPLICIT NONE

       INTEGER ::  i, j, k
!
!--    Adjust number of raindrops to avoid nonlinear effects in 
!--    sedimentation and evaporation of rain drops due to too small or
!--    too big weights of rain drops (Stevens and Seifert, 2008).
!--    The same procedure is applied to cloud droplets if they are determined
!--    prognostically.  
       DO  k = nzb_s_inner(j,i)+1, nzt

          IF ( qr(k,j,i) <= eps_sb )  THEN
             qr(k,j,i) = 0.0
             nr(k,j,i) = 0.0
          ELSE
!
!--          Adjust number of raindrops to avoid nonlinear effects in 
!--          sedimentation and evaporation of rain drops due to too small or
!--          too big weights of rain drops (Stevens and Seifert, 2008).
             IF ( nr(k,j,i) * xrmin > qr(k,j,i) * hyrho(k) )  THEN
                nr(k,j,i) = qr(k,j,i) * hyrho(k) / xrmin
             ELSEIF ( nr(k,j,i) * xrmax < qr(k,j,i) * hyrho(k) )  THEN
                nr(k,j,i) = qr(k,j,i) * hyrho(k) / xrmax
             ENDIF

          ENDIF

       ENDDO

    END SUBROUTINE adjust_cloud_ij


    SUBROUTINE autoconversion_ij( i, j )

       USE arrays_3d
       USE cloud_parameters
       USE control_parameters
       USE grid_variables
       USE indices

       IMPLICIT NONE

       INTEGER ::  i, j, k
       REAL    ::  alpha_cc, autocon, epsilon, k_au, l_mix, nu_c, phi_au,      &
                   r_cc, rc, re_lambda, selfcoll, sigma_cc, tau_cloud, xc                      


       k_au = k_cc / ( 20.0 * x0 )

       DO  k = nzb_s_inner(j,i)+1, nzt

          IF ( qc_1d(k) > eps_sb )  THEN
!
!--          Intern time scale of coagulation (Seifert and Beheng, 2006):
!--          (1.0 - qc(k,j,i) / ( qc(k,j,i) + qr_1d(k) ))
             tau_cloud = 1.0 - qc_1d(k) / ( qr_1d(k) + qc_1d(k) )
!
!--          Universal function for autoconversion process 
!--          (Seifert and Beheng, 2006):
             phi_au    = 600.0 * tau_cloud**0.68 * ( 1.0 - tau_cloud**0.68 )**3
!
!--          Shape parameter of gamma distribution (Geoffroy et al., 2010):
!--          (Use constant nu_c = 1.0 instead?)
             nu_c      = 1.0 !MAX( 0.0, 1580.0 * hyrho(k) * qc(k,j,i) - 0.28 )
!
!--          Mean weight of cloud droplets:
             xc = hyrho(k) * qc_1d(k) / nc_1d(k)
!
!--          Parameterized turbulence effects on autoconversion (Seifert, 
!--          Nuijens and Stevens, 2010)
             IF ( turbulence )  THEN
!
!--             Weight averaged radius of cloud droplets:
                rc = 0.5 * ( xc * dpirho_l )**( 1.0 / 3.0 )

                alpha_cc = ( a_1 + a_2 * nu_c ) / ( 1.0 + a_3 * nu_c )
                r_cc     = ( b_1 + b_2 * nu_c ) / ( 1.0 + b_3 * nu_c )
                sigma_cc = ( c_1 + c_2 * nu_c ) / ( 1.0 + c_3 * nu_c )
!
!--             Mixing length (neglecting distance to ground and stratification)
                l_mix = ( dx * dy * dzu(k) )**( 1.0 / 3.0 )
!
!--             Limit dissipation rate according to Seifert, Nuijens and 
!--             Stevens (2010)
                epsilon = MIN( 0.06, diss(k,j,i) )
!
!--             Compute Taylor-microscale Reynolds number:
                re_lambda = 6.0 / 11.0 * ( l_mix / c_const )**( 2.0 / 3.0 ) *  &
                            SQRT( 15.0 / kin_vis_air ) * epsilon**( 1.0 / 6.0 )
!
!--             The factor of 1.0E4 is needed to convert the dissipation rate
!--             from m2 s-3 to cm2 s-3.
                k_au = k_au * ( 1.0 +                                          &
                       epsilon * 1.0E4 * ( re_lambda * 1.0E-3 )**0.25 *        &
                       ( alpha_cc * EXP( -1.0 * ( ( rc - r_cc ) /              &
                       sigma_cc )**2 ) + beta_cc ) )
             ENDIF
!
!--          Autoconversion rate (Seifert and Beheng, 2006):
             autocon = k_au * ( nu_c + 2.0 ) * ( nu_c + 4.0 ) /    &
                       ( nu_c + 1.0 )**2 * qc_1d(k)**2 * xc**2 *   &
                       ( 1.0 + phi_au / ( 1.0 - tau_cloud )**2 ) * &
                       rho_surface
             autocon = MIN( autocon, qc_1d(k) / dt_micro )

             qr_1d(k) = qr_1d(k) + autocon * dt_micro
             qc_1d(k) = qc_1d(k) - autocon * dt_micro 
             nr_1d(k) = nr_1d(k) + autocon / x0 * hyrho(k) * dt_micro

          ENDIF

       ENDDO

    END SUBROUTINE autoconversion_ij


    SUBROUTINE accretion_ij( i, j )

       USE arrays_3d
       USE cloud_parameters
       USE control_parameters
       USE indices

       IMPLICIT NONE

       INTEGER ::  i, j, k
       REAL    ::  accr, k_cr, phi_ac, tau_cloud, xc

       DO  k = nzb_s_inner(j,i)+1, nzt
          IF ( ( qc_1d(k) > eps_sb )  .AND.  ( qr_1d(k) > eps_sb ) )  THEN
!
!--          Intern time scale of coagulation (Seifert and Beheng, 2006):
             tau_cloud = 1.0 - qc_1d(k) / ( qc_1d(k) + qr_1d(k) ) 
!
!--          Universal function for accretion process 
!--          (Seifert and Beheng, 2001):
             phi_ac = tau_cloud / ( tau_cloud + 5.0E-5 ) 
             phi_ac = ( phi_ac**2 )**2
!
!--          Parameterized turbulence effects on autoconversion (Seifert, 
!--          Nuijens and Stevens, 2010). The factor of 1.0E4 is needed to 
!--          convert the dissipation (diss) from m2 s-3 to cm2 s-3.
             IF ( turbulence )  THEN
                k_cr = k_cr0 * ( 1.0 + 0.05 *                             &
                                 MIN( 600.0, diss(k,j,i) * 1.0E4 )**0.25 )
             ELSE
                k_cr = k_cr0                       
             ENDIF
!
!--          Accretion rate (Seifert and Beheng, 2006):
             accr = k_cr * qc_1d(k) * qr_1d(k) * phi_ac *                 &
                    SQRT( rho_surface * hyrho(k) )
             accr = MIN( accr, qc_1d(k) / dt_micro )

             qr_1d(k) = qr_1d(k) + accr * dt_micro 
             qc_1d(k) = qc_1d(k) - accr * dt_micro

          ENDIF

       ENDDO

    END SUBROUTINE accretion_ij


    SUBROUTINE selfcollection_breakup_ij( i, j )

       USE arrays_3d
       USE cloud_parameters
       USE control_parameters
       USE indices
   
       IMPLICIT NONE

       INTEGER ::  i, j, k
       REAL    ::  breakup, dr, phi_br, selfcoll

       DO  k = nzb_s_inner(j,i)+1, nzt
          IF ( qr_1d(k) > eps_sb )  THEN
!
!--          Selfcollection rate (Seifert and Beheng, 2001):
             selfcoll = k_rr * nr_1d(k) * qr_1d(k) *         &
                        SQRT( hyrho(k) * rho_surface )
!
!--          Weight averaged diameter of rain drops:
             dr = ( hyrho(k) * qr_1d(k) / nr_1d(k) * dpirho_l )**( 1.0 / 3.0 )
!
!--          Collisional breakup rate (Seifert, 2008):
             IF ( dr >= 0.3E-3 )  THEN
                phi_br  = k_br * ( dr - 1.1E-3 )
                breakup = selfcoll * ( phi_br + 1.0 )
             ELSE
                breakup = 0.0
             ENDIF

             selfcoll = MAX( breakup - selfcoll, -nr_1d(k) / dt_micro )
             nr_1d(k) = nr_1d(k) + selfcoll * dt_micro

          ENDIF          
       ENDDO

    END SUBROUTINE selfcollection_breakup_ij


    SUBROUTINE evaporation_rain_ij( i, j )
!
!--    Evaporation of precipitable water. Condensation is neglected for 
!--    precipitable water.

       USE arrays_3d
       USE cloud_parameters
       USE constants
       USE control_parameters
       USE indices

       IMPLICIT NONE

       INTEGER ::  i, j, k
       REAL    ::  alpha, dr, e_s, evap, evap_nr, f_vent, g_evap, lambda_r, &
                   mu_r, mu_r_2, mu_r_5d2, nr_0, q_s, sat, t_l, temp, xr

       DO  k = nzb_s_inner(j,i)+1, nzt
          IF ( qr_1d(k) > eps_sb )  THEN
!
!--          Actual liquid water temperature:
             t_l = t_d_pt(k) * pt_1d(k)
!
!--          Saturation vapor pressure at t_l:
             e_s = 610.78 * EXP( 17.269 * ( t_l - 273.16 ) / ( t_l - 35.86 ) )
!
!--          Computation of saturation humidity:
             q_s = 0.622 * e_s / ( hyp(k) - 0.378 * e_s )
             alpha = 0.622 * l_d_r * l_d_cp / ( t_l * t_l )
             q_s = q_s * ( 1.0 + alpha * q_1d(k) ) / ( 1.0 + alpha * q_s )
!
!--          Supersaturation:
             sat = MIN( 0.0, ( q_1d(k) - qr_1d(k) - qc_1d(k) ) / q_s - 1.0 )
!
!--          Actual temperature:
             temp = t_l + l_d_cp * ( qc_1d(k) + qr_1d(k) )
    
             g_evap = 1.0 / ( ( l_v / ( r_v * temp ) - 1.0 ) * l_v /   &
                      ( thermal_conductivity_l * temp ) + r_v * temp / &
                      ( diff_coeff_l * e_s ) )
!
!--          Mean weight of rain drops
             xr = hyrho(k) * qr_1d(k) / nr_1d(k)
!
!--          Weight averaged diameter of rain drops:
             dr = ( xr * dpirho_l )**( 1.0 / 3.0 )
!
!--          Compute ventilation factor and intercept parameter 
!--          (Seifert and Beheng, 2006; Seifert, 2008):
             IF ( ventilation_effect )  THEN
!
!--             Shape parameter of gamma distribution (Milbrandt and Yau, 2005;
!--             Stevens and Seifert, 2008):
                mu_r = 10.0 * ( 1.0 + TANH( 1.2E3 * ( dr - 1.4E-3 ) ) )
!
!--             Slope parameter of gamma distribution (Seifert, 2008):
                lambda_r = ( ( mu_r + 3.0 ) * ( mu_r + 2.0 ) *       &
                             ( mu_r + 1.0 ) )**( 1.0 / 3.0 ) / dr

                mu_r_2   = mu_r + 2.0
                mu_r_5d2 = mu_r + 2.5 
                f_vent = a_vent * gamm( mu_r_2 ) *                            &
                         lambda_r**( -mu_r_2 ) +                              &
                         b_vent * schmidt_p_1d3 *                             &
                         SQRT( a_term / kin_vis_air ) * gamm( mu_r_5d2 ) *    &
                         lambda_r**( -mu_r_5d2 ) *                            &
                         ( 1.0 - 0.5 * ( b_term / a_term ) *                  &
                         ( lambda_r /                                         &
                         (       c_term + lambda_r ) )**mu_r_5d2 -            &
                                 0.125 * ( b_term / a_term )**2 *             &
                         ( lambda_r /                                         &
                         ( 2.0 * c_term + lambda_r ) )**mu_r_5d2 -            &
                                 0.0625 * ( b_term / a_term )**3 *            &
                         ( lambda_r /                                         &
                         ( 3.0 * c_term + lambda_r ) )**mu_r_5d2 -            &
                                 0.0390625 * ( b_term / a_term )**4 *         &
                         ( lambda_r /                                         &
                         ( 4.0 * c_term + lambda_r ) )**mu_r_5d2 )
                nr_0   = nr_1d(k) * lambda_r**( mu_r + 1.0 ) /                &
                         gamm( mu_r + 1.0 ) 
             ELSE
                f_vent = 1.0
                nr_0   = nr_1d(k) * dr
             ENDIF
!
!--          Evaporation rate of rain water content (Seifert and Beheng, 2006):
             evap = 2.0 * pi * nr_0 * g_evap * f_vent * sat /    &
                    hyrho(k)

             evap    = MAX( evap, -qr_1d(k) / dt_micro )
             evap_nr = MAX( c_evap * evap / xr * hyrho(k), &
                            -nr_1d(k) / dt_micro )

             qr_1d(k) = qr_1d(k) + evap * dt_micro
             nr_1d(k) = nr_1d(k) + evap_nr * dt_micro
          ENDIF          

       ENDDO

    END SUBROUTINE evaporation_rain_ij


    SUBROUTINE sedimentation_cloud_ij( i, j )

       USE arrays_3d
       USE cloud_parameters
       USE constants
       USE control_parameters
       USE indices
       
       IMPLICIT NONE

       INTEGER ::  i, j, k
       REAL    ::  sed_qc_const

       REAL, DIMENSION(nzb:nzt+1) :: sed_qc

!
!--    Sedimentation of cloud droplets (Heus et al., 2010):
       sed_qc_const = k_st * ( 3.0 / ( 4.0 * pi * rho_l ))**( 2.0 / 3.0 ) *   &
                     EXP( 5.0 * LOG( sigma_gc )**2 )

       sed_qc(nzt+1) = 0.0

       DO  k = nzt, nzb_s_inner(j,i)+1, -1
          IF ( qc_1d(k) > eps_sb )  THEN
             sed_qc(k) = sed_qc_const * nc_1d(k)**( -2.0 / 3.0 ) * &
                        ( qc_1d(k) * hyrho(k) )**( 5.0 / 3.0 )
          ELSE
             sed_qc(k) = 0.0
          ENDIF

          sed_qc(k) = MIN( sed_qc(k), hyrho(k) * dzu(k+1) * q_1d(k) /     &
                                      dt_micro + sed_qc(k+1) )

          q_1d(k)  = q_1d(k)  + ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) /  &
                                hyrho(k) * dt_micro
          qc_1d(k) = qc_1d(k) + ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) / & 
                                hyrho(k) * dt_micro
          pt_1d(k) = pt_1d(k) - ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) / &
                                hyrho(k) * l_d_cp * pt_d_t(k) * dt_micro

       ENDDO

    END SUBROUTINE sedimentation_cloud_ij


    SUBROUTINE sedimentation_rain_ij( i, j )

       USE arrays_3d
       USE cloud_parameters
       USE constants
       USE control_parameters
       USE indices
       USE statistics
       
       IMPLICIT NONE

       INTEGER ::  i, j, k, k_run
       REAL    ::  c_run, d_max, d_mean, d_min, dr, dt_sedi, flux, lambda_r,  &
                   mu_r, z_run

       REAL, DIMENSION(nzb:nzt+1) :: c_nr, c_qr, d_nr, d_qr, nr_slope,        &
                                     qr_slope, sed_nr, sed_qr, w_nr, w_qr 
!
!--    Computation of sedimentation flux. Implementation according to Stevens
!--    and Seifert (2008).
       IF ( intermediate_timestep_count == 1 )  prr(:,j,i) = 0.0
!
!--    Compute velocities 
       DO  k = nzb_s_inner(j,i)+1, nzt
          IF ( qr_1d(k) > eps_sb )  THEN
!
!--          Weight averaged diameter of rain drops:
             dr = ( hyrho(k) * qr_1d(k) / nr_1d(k) * dpirho_l )**( 1.0 / 3.0 )
!
!--          Shape parameter of gamma distribution (Milbrandt and Yau, 2005;
!--          Stevens and Seifert, 2008):
             mu_r = 10.0 * ( 1.0 + TANH( 1.2E3 * ( dr - 1.4E-3 ) ) )
!
!--          Slope parameter of gamma distribution (Seifert, 2008):
             lambda_r = ( ( mu_r + 3.0 ) * ( mu_r + 2.0 ) *          &
                        ( mu_r + 1.0 ) )**( 1.0 / 3.0 ) / dr

             w_nr(k) = MAX( 0.1, MIN( 20.0, a_term - b_term * ( 1.0 +          &
                       c_term / lambda_r )**( -1.0 * ( mu_r + 1.0 ) ) ) )
             w_qr(k) = MAX( 0.1, MIN( 20.0, a_term - b_term * ( 1.0 +          &
                       c_term / lambda_r )**( -1.0 * ( mu_r + 4.0 ) ) ) )
          ELSE
             w_nr(k) = 0.0
             w_qr(k) = 0.0
          ENDIF
       ENDDO
!
!--    Adjust boundary values
       w_nr(nzb_s_inner(j,i)) = w_nr(nzb_s_inner(j,i)+1)
       w_qr(nzb_s_inner(j,i)) = w_qr(nzb_s_inner(j,i)+1)
       w_nr(nzt+1) = 0.0
       w_qr(nzt+1) = 0.0
!
!--    Compute Courant number
       DO  k = nzb_s_inner(j,i)+1, nzt
          c_nr(k) = 0.25 * ( w_nr(k-1) + 2.0 * w_nr(k) + w_nr(k+1) ) * &
                    dt_micro * ddzu(k)
          c_qr(k) = 0.25 * ( w_qr(k-1) + 2.0 * w_qr(k) + w_qr(k+1) ) * &
                    dt_micro * ddzu(k)
       ENDDO     
!
!--    Limit slopes with monotonized centered (MC) limiter (van Leer, 1977):
       IF ( limiter_sedimentation )  THEN

          DO k = nzb_s_inner(j,i)+1, nzt
             d_mean = 0.5 * ( qr_1d(k+1) + qr_1d(k-1) )
             d_min  = qr_1d(k) - MIN( qr_1d(k+1), qr_1d(k), qr_1d(k-1) )
             d_max  = MAX( qr_1d(k+1), qr_1d(k), qr_1d(k-1) ) - qr_1d(k)

             qr_slope(k) = SIGN(1.0, d_mean) * MIN ( 2.0 * d_min, 2.0 * d_max, &
                                                     ABS( d_mean ) )

             d_mean = 0.5 * ( nr_1d(k+1) + nr_1d(k-1) )
             d_min  = nr_1d(k) - MIN( nr_1d(k+1), nr_1d(k), nr_1d(k-1) )
             d_max  = MAX( nr_1d(k+1), nr_1d(k), nr_1d(k-1) ) - nr_1d(k)

             nr_slope(k) = SIGN(1.0, d_mean) * MIN ( 2.0 * d_min, 2.0 * d_max, &
                                                     ABS( d_mean ) )
          ENDDO

       ELSE

          nr_slope = 0.0
          qr_slope = 0.0

       ENDIF

       sed_nr(nzt+1) = 0.0
       sed_qr(nzt+1) = 0.0
!
!--    Compute sedimentation flux
       DO  k = nzt, nzb_s_inner(j,i)+1, -1
!
!--       Sum up all rain drop number densities which contribute to the flux 
!--       through k-1/2
          flux  = 0.0
          z_run = 0.0 ! height above z(k)
          k_run = k
          c_run = MIN( 1.0, c_nr(k) )
          DO WHILE ( c_run > 0.0  .AND.  k_run <= nzt )
             flux  = flux + hyrho(k_run) *                                    &
                     ( nr_1d(k_run) + nr_slope(k_run) * ( 1.0 - c_run ) *     &
                     0.5 ) * c_run * dzu(k_run)
             z_run = z_run + dzu(k_run)
             k_run = k_run + 1
             c_run = MIN( 1.0, c_nr(k_run) - z_run * ddzu(k_run) )
          ENDDO
!
!--       It is not allowed to sediment more rain drop number density than 
!--       available
          flux = MIN( flux,                                                   &
                      hyrho(k) * dzu(k+1) * nr_1d(k) + sed_nr(k+1) * dt_micro )

          sed_nr(k) = flux / dt_micro
          nr_1d(k)  = nr_1d(k) + ( sed_nr(k+1) - sed_nr(k) ) * ddzu(k+1) /    &
                                 hyrho(k) * dt_micro
!
!--       Sum up all rain water content which contributes to the flux 
!--       through k-1/2
          flux  = 0.0
          z_run = 0.0 ! height above z(k)
          k_run = k
          c_run = MIN( 1.0, c_qr(k) )

          DO WHILE ( c_run > 0.0  .AND.  k_run <= nzt-1 )

             flux  = flux + hyrho(k_run) *                                    &
                     ( qr_1d(k_run) + qr_slope(k_run) * ( 1.0 - c_run ) *    &
                     0.5 ) * c_run * dzu(k_run)
             z_run = z_run + dzu(k_run)
             k_run = k_run + 1
             c_run = MIN( 1.0, c_qr(k_run) - z_run * ddzu(k_run) )

          ENDDO
!
!--       It is not allowed to sediment more rain water content than available
          flux = MIN( flux,                                                   &
                      hyrho(k) * dzu(k) * qr_1d(k) + sed_qr(k+1) * dt_micro )

          sed_qr(k) = flux / dt_micro

          qr_1d(k) = qr_1d(k) + ( sed_qr(k+1) - sed_qr(k) ) * ddzu(k+1) / &
                                hyrho(k) * dt_micro
          q_1d(k)  = q_1d(k)  + ( sed_qr(k+1) - sed_qr(k) ) * ddzu(k+1) / &
                                hyrho(k) * dt_micro 
          pt_1d(k) = pt_1d(k) - ( sed_qr(k+1) - sed_qr(k) ) * ddzu(k+1) / &
                                hyrho(k) * l_d_cp * pt_d_t(k) * dt_micro
!
!--       Compute the rain rate
          prr(k,j,i) = prr(k,j,i) + sed_qr(k) / hyrho(k) *                    &
                       weight_substep(intermediate_timestep_count)
       ENDDO

!
!--    Precipitation amount
       IF ( intermediate_timestep_count == intermediate_timestep_count_max    &
            .AND.  ( dt_do2d_xy - time_do2d_xy ) <                            &
            precipitation_amount_interval )  THEN

          precipitation_amount(j,i) = precipitation_amount(j,i) +   &
                                      prr(nzb_s_inner(j,i)+1,j,i) *      &
                                      hyrho(nzb_s_inner(j,i)+1) * dt_3d
       ENDIF

    END SUBROUTINE sedimentation_rain_ij


!
!-- This function computes the gamma function (Press et al., 1992). 
!-- The gamma function is needed for the calculation of the evaporation 
!-- of rain drops.
    FUNCTION gamm( xx ) 
       
       USE cloud_parameters
    
       IMPLICIT NONE 
 
       REAL    ::  gamm, ser, tmp, x_gamm, xx, y_gamm
       INTEGER ::  j 

 
       x_gamm = xx 
       y_gamm = x_gamm 
       tmp = x_gamm + 5.5 
       tmp = ( x_gamm + 0.5 ) * LOG( tmp ) - tmp 
       ser = 1.000000000190015 

       DO  j = 1, 6 
          y_gamm = y_gamm + 1.0 
          ser    = ser + cof( j ) / y_gamm 
       ENDDO

! 
!--    Until this point the algorithm computes the logarithm of the gamma 
!--    function. Hence, the exponential function is used.  
!       gamm = EXP( tmp + LOG( stp * ser / x_gamm ) ) 
       gamm = EXP( tmp ) * stp * ser / x_gamm 

       RETURN 
 
    END FUNCTION gamm 

 END MODULE microphysics_mod