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 1242 2013-10-30 11:50:11Z fricke $ ! ! 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