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 . ! ! Copyright 1997-2012 Leibniz University Hannover !--------------------------------------------------------------------------------! ! ! Current revisions: ! ------------------ ! ! ! Former revisions: ! ----------------- ! $Id: microphysics.f90 1107 2013-03-04 06:23:14Z raasch $ ! ! 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 dsd_properties, autoconversion, accretion, selfcollection_breakup, & evaporation_rain, sedimentation_cloud, sedimentation_rain INTERFACE dsd_properties MODULE PROCEDURE dsd_properties MODULE PROCEDURE dsd_properties_ij END INTERFACE dsd_properties 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 dsd_properties USE arrays_3d USE cloud_parameters USE constants USE indices IMPLICIT NONE INTEGER :: i, j, k DO i = nxl, nxr DO j = nys, nyn DO k = nzb_2d(j,i)+1, nzt ENDDO ENDDO ENDDO END SUBROUTINE dsd_properties SUBROUTINE autoconversion USE arrays_3d USE cloud_parameters USE constants USE indices IMPLICIT NONE INTEGER :: i, j, k DO i = nxl, nxr DO j = nys, nyn DO k = nzb_2d(j,i)+1, nzt ENDDO ENDDO ENDDO END SUBROUTINE autoconversion SUBROUTINE accretion USE arrays_3d USE cloud_parameters USE constants USE indices IMPLICIT NONE INTEGER :: i, j, k DO i = nxl, nxr DO j = nys, nyn DO k = nzb_2d(j,i)+1, nzt ENDDO ENDDO ENDDO END SUBROUTINE accretion SUBROUTINE selfcollection_breakup USE arrays_3d USE cloud_parameters USE constants USE indices IMPLICIT NONE INTEGER :: i, j, k DO i = nxl, nxr DO j = nys, nyn DO k = nzb_2d(j,i)+1, nzt ENDDO ENDDO ENDDO END SUBROUTINE selfcollection_breakup SUBROUTINE evaporation_rain USE arrays_3d USE cloud_parameters USE constants USE indices IMPLICIT NONE INTEGER :: i, j, k DO i = nxl, nxr DO j = nys, nyn DO k = nzb_2d(j,i)+1, nzt ENDDO ENDDO ENDDO END SUBROUTINE evaporation_rain SUBROUTINE sedimentation_cloud USE arrays_3d USE cloud_parameters USE constants USE indices IMPLICIT NONE INTEGER :: i, j, k DO i = nxl, nxr DO j = nys, nyn DO k = nzb_2d(j,i)+1, nzt ENDDO ENDDO ENDDO END SUBROUTINE sedimentation_cloud SUBROUTINE sedimentation_rain USE arrays_3d USE cloud_parameters USE constants USE indices IMPLICIT NONE INTEGER :: i, j, k DO i = nxl, nxr DO j = nys, nyn DO k = nzb_2d(j,i)+1, nzt ENDDO ENDDO ENDDO END SUBROUTINE sedimentation_rain !------------------------------------------------------------------------------! ! Call for grid point i,j !------------------------------------------------------------------------------! SUBROUTINE dsd_properties_ij( i, j ) USE arrays_3d USE cloud_parameters USE constants USE indices USE control_parameters USE user IMPLICIT NONE INTEGER :: i, j, k DO k = nzb_2d(j,i)+1, nzt IF ( qr(k,j,i) <= eps_sb ) THEN qr(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 xr(k) = hyrho(k) * qr(k,j,i) / nr(k,j,i) ! !-- Weight averaged diameter of rain drops: dr(k) = ( hyrho(k) * qr(k,j,i) / nr(k,j,i) * & dpirho_l )**( 1.0 / 3.0 ) ! !-- Shape parameter of gamma distribution (Milbrandt and Yau, 2005; !-- Stevens and Seifert, 2008): mu_r(k) = 10.0 * ( 1.0 + TANH( 1.2E3 * ( dr(k) - 1.4E-3 ) ) ) ! !-- Slope parameter of gamma distribution (Seifert, 2008): lambda_r(k) = ( ( mu_r(k) + 3.0 ) * ( mu_r(k) + 2.0 ) * & ( mu_r(k) + 1.0 ) )**( 1.0 / 3.0 ) / dr(k) ENDIF ENDDO END SUBROUTINE dsd_properties_ij SUBROUTINE autoconversion_ij( i, j ) USE arrays_3d USE cloud_parameters USE constants USE indices USE control_parameters USE statistics USE grid_variables IMPLICIT NONE INTEGER :: i, j, k REAL :: k_au, autocon, phi_au, tau_cloud, xc, nu_c, rc, & l_mix, re_lambda, alpha_cc, r_cc, sigma_cc, epsilon k_au = k_cc / ( 20.0 * x0 ) DO k = nzb_2d(j,i)+1, nzt IF ( ql(k,j,i) > 0.0 ) THEN ! !-- Intern time scale of coagulation (Seifert and Beheng, 2006): !-- (1.0 - ql(k,j,i) / ( ql(k,j,i) + qr(k,j,i) )) tau_cloud = 1.0 - ql(k,j,i) / ( ql(k,j,i) + qr(k,j,i) + 1.0E-20 ) ! !-- 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) * ql(k,j,i) - 0.28 ) ! !-- Mean weight of cloud droplets: xc = hyrho(k) * ql(k,j,i) / nc ! !-- 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 * ql(k,j,i)**2 * xc**2 * & ( 1.0 + phi_au / ( 1.0 - tau_cloud )**2 ) * & rho_surface autocon = MIN( autocon, ql(k,j,i) / ( dt_3d * & weight_substep(intermediate_timestep_count) ) ) ! !-- Tendencies for q, qr, nr, pt: tend_qr(k,j,i) = tend_qr(k,j,i) + autocon tend_q(k,j,i) = tend_q(k,j,i) - autocon tend_nr(k,j,i) = tend_nr(k,j,i) + autocon / x0 * hyrho(k) tend_pt(k,j,i) = tend_pt(k,j,i) + autocon * l_d_cp * pt_d_t(k) ENDIF ENDDO END SUBROUTINE autoconversion_ij SUBROUTINE accretion_ij( i, j ) USE arrays_3d USE cloud_parameters USE constants USE indices USE control_parameters USE statistics IMPLICIT NONE INTEGER :: i, j, k REAL :: accr, phi_ac, tau_cloud, k_cr DO k = nzb_2d(j,i)+1, nzt IF ( ( ql(k,j,i) > 0.0 ) .AND. ( qr(k,j,i) > eps_sb ) ) THEN ! !-- Intern time scale of coagulation (Seifert and Beheng, 2006): tau_cloud = 1.0 - ql(k,j,i) / ( ql(k,j,i) + qr(k,j,i) + 1.0E-20) ! !-- 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 * ql(k,j,i) * qr(k,j,i) * phi_ac * & SQRT( rho_surface * hyrho(k) ) accr = MIN( accr, ql(k,j,i) / ( dt_3d * & weight_substep(intermediate_timestep_count) ) ) ! !-- Tendencies for q, qr, pt: tend_qr(k,j,i) = tend_qr(k,j,i) + accr tend_q(k,j,i) = tend_q(k,j,i) - accr tend_pt(k,j,i) = tend_pt(k,j,i) + accr * l_d_cp * pt_d_t(k) ENDIF ENDDO END SUBROUTINE accretion_ij SUBROUTINE selfcollection_breakup_ij( i, j ) USE arrays_3d USE cloud_parameters USE constants USE indices USE control_parameters USE statistics IMPLICIT NONE INTEGER :: i, j, k REAL :: selfcoll, breakup, phi_br, phi_sc DO k = nzb_2d(j,i)+1, nzt IF ( qr(k,j,i) > eps_sb ) THEN ! !-- Selfcollection rate (Seifert and Beheng, 2006): !-- pirho_l**( 1.0 / 3.0 ) is necessary to convert [lambda_r] = m-1 to !-- kg**( 1.0 / 3.0 ). phi_sc = 1.0 !( 1.0 + kappa_rr / lambda_r(k) * & !pirho_l**( 1.0 / 3.0 ) )**( -9 ) selfcoll = k_rr * nr(k,j,i) * qr(k,j,i) * phi_sc * & SQRT( hyrho(k) * rho_surface ) ! !-- Collisional breakup rate (Seifert, 2008): IF ( dr(k) >= 0.3E-3 ) THEN phi_br = k_br * ( dr(k) - 1.1E-3 ) breakup = selfcoll * ( phi_br + 1.0 ) ELSE breakup = 0.0 ENDIF selfcoll = MAX( breakup - selfcoll, -nr(k,j,i) / ( dt_3d * & weight_substep(intermediate_timestep_count) ) ) ! !-- Tendency for nr: tend_nr(k,j,i) = tend_nr(k,j,i) + selfcoll 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 indices USE control_parameters USE statistics IMPLICIT NONE INTEGER :: i, j, k REAL :: evap, alpha, e_s, q_s, t_l, sat, temp, g_evap, f_vent, & mu_r_2, mu_r_5d2, nr_0 DO k = nzb_2d(j,i)+1, nzt IF ( qr(k,j,i) > eps_sb ) THEN ! !-- Actual liquid water temperature: t_l = t_d_pt(k) * pt(k,j,i) ! !-- 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(k,j,i) ) / ( 1.0 + alpha * q_s ) ! !-- Supersaturation: sat = MIN( 0.0, ( q(k,j,i) - ql(k,j,i) ) / q_s - 1.0 ) ! !-- Actual temperature: temp = t_l + l_d_cp * ql(k,j,i) ! !-- g_evap = ( l_v / ( r_v * temp ) - 1.0 ) * l_v / & ( thermal_conductivity_l * temp ) + rho_l * r_v * temp /& ( diff_coeff_l * e_s ) g_evap = 1.0 / g_evap ! !-- Compute ventilation factor and intercept parameter !-- (Seifert and Beheng, 2006; Seifert, 2008): IF ( ventilation_effect ) THEN mu_r_2 = mu_r(k) + 2.0 mu_r_5d2 = mu_r(k) + 2.5 f_vent = a_vent * gamm( mu_r_2 ) * & lambda_r(k)**( -mu_r_2 ) + & b_vent * schmidt_p_1d3 * & SQRT( a_term / kin_vis_air ) * gamm( mu_r_5d2 ) * & lambda_r(k)**( -mu_r_5d2 ) * & ( 1.0 - 0.5 * ( b_term / a_term ) * & ( lambda_r(k) / & ( c_term + lambda_r(k) ) )**mu_r_5d2 - & 0.125 * ( b_term / a_term )**2 * & ( lambda_r(k) / & ( 2.0 * c_term + lambda_r(k) ) )**mu_r_5d2 - & 0.0625 * ( b_term / a_term )**3 * & ( lambda_r(k) / & ( 3.0 * c_term + lambda_r(k) ) )**mu_r_5d2 - & 0.0390625 * ( b_term / a_term )**4 * & ( lambda_r(k) / & ( 4.0 * c_term + lambda_r(k) ) )**mu_r_5d2 ) nr_0 = nr(k,j,i) * lambda_r(k)**( mu_r(k) + 1.0 ) / & gamm( mu_r(k) + 1.0 ) ELSE f_vent = 1.0 nr_0 = nr(k,j,i) * dr(k) 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(k,j,i) / ( dt_3d * & weight_substep(intermediate_timestep_count) ) ) ! !-- Tendencies for q, qr, nr, pt: tend_qr(k,j,i) = tend_qr(k,j,i) + evap tend_q(k,j,i) = tend_q(k,j,i) - evap tend_nr(k,j,i) = tend_nr(k,j,i) + c_evap * evap / xr(k) * hyrho(k) tend_pt(k,j,i) = tend_pt(k,j,i) + evap * l_d_cp * pt_d_t(k) ENDIF ENDDO END SUBROUTINE evaporation_rain_ij SUBROUTINE sedimentation_cloud_ij( i, j ) USE arrays_3d USE cloud_parameters USE constants USE indices USE control_parameters IMPLICIT NONE INTEGER :: i, j, k REAL :: sed_q_const, sigma_gc = 1.3, k_st = 1.2E8 ! !-- Sedimentation of cloud droplets (Heus et al., 2010): sed_q_const = k_st * ( 3.0 / ( 4.0 * pi * rho_l ))**( 2.0 / 3.0 ) * & EXP( 5.0 * LOG( sigma_gc )**2 ) sed_q = 0.0 DO k = nzb_2d(j,i)+1, nzt IF ( ql(k,j,i) > 0.0 ) THEN sed_q(k) = sed_q_const * nc**( -2.0 / 3.0 ) * & ( ql(k,j,i) * hyrho(k) )**( 5.0 / 3.0 ) ENDIF ENDDO ! !-- Tendency for q, pt: DO k = nzb_2d(j,i)+1, nzt tend_q(k,j,i) = tend_q(k,j,i) + ( sed_q(k+1) - sed_q(k) ) * & ddzu(k+1) / hyrho(k) tend_pt(k,j,i) = tend_pt(k,j,i) - ( sed_q(k+1) - sed_q(k) ) * & ddzu(k+1) / hyrho(k) * l_d_cp * pt_d_t(k) ENDDO END SUBROUTINE sedimentation_cloud_ij SUBROUTINE sedimentation_rain_ij( i, j ) USE arrays_3d USE cloud_parameters USE constants USE indices USE control_parameters USE statistics IMPLICIT NONE INTEGER :: i, j, k, k_run REAL :: c_run, d_max, d_mean, d_min, dt_sedi, flux, z_run REAL, DIMENSION(nzb:nzt) :: c_nr, c_qr, nr_slope, qr_slope, 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 dt_sedi = dt_3d * weight_substep(intermediate_timestep_count) w_nr = 0.0 w_qr = 0.0 ! !-- Compute velocities DO k = nzb_s_inner(j,i)+1, nzt IF ( qr(k,j,i) > eps_sb ) THEN w_nr(k) = MAX( 0.1, MIN( 20.0, a_term - b_term * ( 1.0 + & c_term / lambda_r(k) )**( -1.0 * ( mu_r(k) + 1.0 ) ) ) ) w_qr(k) = MAX( 0.1, MIN( 20.0, a_term - b_term * ( 1.0 + & c_term / lambda_r(k) )**( -1.0 * ( mu_r(k) + 4.0 ) ) ) ) ELSE w_nr(k) = 0.0 w_qr(k) = 0.0 ENDIF ENDDO ! !-- Adjust boundary values w_nr(nzb_2d(j,i)) = w_nr(nzb_2d(j,i)+1) w_qr(nzb_2d(j,i)) = w_qr(nzb_2d(j,i)+1) w_nr(nzt) = w_nr(nzt-1) w_qr(nzt) = w_qr(nzt-1) ! !-- Compute Courant number DO k = nzb_s_inner(j,i)+1, nzt-1 c_nr(k) = 0.25 * ( w_nr(k-1) + 2.0 * w_nr(k) + w_nr(k+1) ) * & dt_sedi * ddzu(k) c_qr(k) = 0.25 * ( w_qr(k-1) + 2.0 * w_qr(k) + w_qr(k+1) ) * & dt_sedi * ddzu(k) ENDDO ! !-- Limit slopes with monotonized centered (MC) limiter (van Leer, 1977): IF ( limiter_sedimentation ) THEN qr(nzb_s_inner(j,i),j,i) = 0.0 nr(nzb_s_inner(j,i),j,i) = 0.0 qr(nzt,j,i) = 0.0 nr(nzt,j,i) = 0.0 DO k = nzb_s_inner(j,i)+1, nzt-1 d_mean = 0.5 * ( qr(k+1,j,i) + qr(k-1,j,i) ) d_min = qr(k,j,i) - MIN( qr(k+1,j,i), qr(k,j,i), qr(k-1,j,i) ) d_max = MAX( qr(k+1,j,i), qr(k,j,i), qr(k-1,j,i) ) - qr(k,j,i) 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(k+1,j,i) + nr(k-1,j,i) ) d_min = nr(k,j,i) - MIN( nr(k+1,j,i), nr(k,j,i), nr(k-1,j,i) ) d_max = MAX( nr(k+1,j,i), nr(k,j,i), nr(k-1,j,i) ) - nr(k,j,i) 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 ! !-- Compute sedimentation flux DO k = nzt-2, 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-1 ) flux = flux + hyrho(k_run) * & ( nr(k_run,j,i) + 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) * nr(k,j,i) + sed_nr(k+1) * dt_sedi ) sed_nr(k) = flux / dt_sedi tend_nr(k,j,i) = tend_nr(k,j,i) + ( sed_nr(k+1) - sed_nr(k) ) * & ddzu(k+1) / hyrho(k) ! !-- 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(k_run,j,i) + 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(k,j,i) + sed_qr(k+1) * dt_sedi ) sed_qr(k) = flux / dt_sedi tend_qr(k,j,i) = tend_qr(k,j,i) + ( sed_qr(k+1) - sed_qr(k) ) * & ddzu(k+1) / hyrho(k) ! !-- 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_2d(j,i)+1,j,i) * & hyrho(nzb_2d(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