MODULE microphysics_mod !------------------------------------------------------------------------------! ! Current revisions: ! ----------------- ! ! Former revisions: ! ----------------- ! $Id$ ! ! 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 REAL :: dqdt_precip 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 REAL :: dqdt_precip 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 REAL :: dqdt_precip 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 REAL :: dqdt_precip 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 REAL :: dqdt_precip 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 REAL :: dqdt_precip 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 REAL :: dqdt_precip 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 REAL :: dr_min = 2.0E-6, dr_max = 1.0E-3 DO k = nzb_2d(j,i)+1, nzt IF ( ( qr(k,j,i) > eps_sb ) .AND. ( nr(k,j,i) > eps_sb ) ) THEN ! !-- Weight averaged diameter of rain drops: dr(k) = ( hyrho(k) * qr(k,j,i) / nr(k,j,i) * & dpirho_l )**( 1.0 / 3.0 ) ! !-- Adjust number of raindrops to avoid nonlinear effects in !-- sedimentation and evaporation of rain drops due to too small or !-- big diameters of rain drops (Stevens and Seifert, 2008). IF ( dr(k) < dr_min ) THEN nr(k,j,i) = qr(k,j,i) * hyrho(k) / dr_min**3 * dpirho_l dr(k) = dr_min ELSEIF ( dr(k) > dr_max ) THEN nr(k,j,i) = qr(k,j,i) * hyrho(k) / dr_max**3 * dpirho_l dr(k) = dr_max ENDIF ! !-- Mean weight of rain drops (Seifert and Beheng, 2006): xr(k) = MIN( MAX( hyrho(k) * qr(k,j,i) / nr(k,j,i), xrmin ), xrmax) ! !-- Shape parameter of gamma distribution (Milbrandt and Yau, 2005; !-- Stevens and Seifert, 2008): IF ( .NOT. mu_constant ) THEN mu_r(k) = 10.0 * ( 1.0 + TANH( 1.2E3 * ( dr(k) - 1.4E-3 ) ) ) ELSE mu_r(k) = mu_constant_value ENDIF ! !-- 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 IMPLICIT NONE INTEGER :: i, j, k REAL :: k_au, autocon, phi_au, tau_cloud, xc, nu_c 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 = 1580.0 * hyrho(k) * ql(k,j,i) - 0.28 ! !-- Mean weight of cloud droplets: xc = hyrho(k) * ql(k,j,i) / nc ! !-- 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 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 ! !-- 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 ) .AND. ( nr(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 + 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 ) .AND. ( nr(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 ) ! !-- Oversaturation: 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, n, n_substep REAL :: sed_nr_tend, sed_qr_tend IF ( intermediate_timestep_count == 1 ) prr(:,j,i) = 0.0 sed_nr = 0.0 sed_qr = 0.0 DO k = nzb_2d(j,i)+1, nzt IF ( ( qr(k,j,i) > eps_sb ) .AND. ( nr(k,j,i) > eps_sb ) ) THEN ! !-- Sedimentation of rain water content and rain drop concentration !-- according to Stevens and Seifert (2008): sed_nr(k) = MIN( 20.0, MAX( 0.1, a_term - b_term * ( 1.0 + & c_term / lambda_r(k) )**( -1.0 * & ( mu_r(k) + 1.0 ) ) ) ) * nr(k,j,i) sed_qr(k) = MIN( 20.0, MAX( 0.1, a_term - b_term * ( 1.0 + & c_term / lambda_r(k) )**( -1.0 * & ( mu_r(k) + 4.0 ) ) ) ) * qr(k,j,i) * hyrho(k) ! !-- Computation of rain rate prr(k,j,i) = prr(k,j,i) + sed_qr(k) / hyrho(k) * & weight_substep(intermediate_timestep_count) ENDIF ENDDO DO k = nzb_2d(j,i)+1, nzt sed_nr_tend = MAX( ( sed_nr(k+1) - sed_nr(k) ) * ddzu(k+1), & -nr(k,j,i) / ( dt_3d * & weight_substep(intermediate_timestep_count) ) ) sed_qr_tend = MAX( ( sed_qr(k+1) - sed_qr(k) ) * ddzu(k+1) / & hyrho(k), & -qr(k,j,i) / ( dt_3d * & weight_substep(intermediate_timestep_count) ) ) tend_nr(k,j,i) = tend_nr(k,j,i) + sed_nr_tend tend_qr(k,j,i) = tend_qr(k,j,i) + sed_qr_tend 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, xx, & ser, tmp, x_gamm, 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