MODULE MP2M_HAZE !============================================================================ ! ! Purpose ! ------- ! Haze microphysics module. ! ! This module contains all definitions of the microphysics processes related to ! aerosols (production, coagulation, sedimentation). ! ! @Warning ! The tendencies returned by the method are always defined over the vertical grid ! from __TOP__ to __GROUND__. ! ! The module also contains sixteen method: ! - mm_haze_microphysics ! - mm_haze_production ! - mm_haze_sedimentation ! * get_weff ! * let_me_fall_in_peace ! - mm_haze_coagulation ! * g0ssco ! * g0sfco ! * g0ffco ! * g3ssco ! * g3sfco ! * g0ssfm ! * g0sffm ! * g0fffm ! * g3ssfm ! * g3sffm ! ! Authors ! ------- ! B. de Batz de Trenquelléon, J. Burgalat (11/2024) ! !============================================================================ USE MP2M_MPREC USE MP2M_GLOBALS USE MP2M_METHODS IMPLICIT NONE PRIVATE PUBLIC :: mm_haze_microphysics, & mm_haze_production, & mm_haze_sedimentation, & mm_haze_coagulation CONTAINS !============================================================================ ! HAZE MICROPHYSICS INTERFACE SUBROUTINE !============================================================================ SUBROUTINE mm_haze_microphysics(m3a_s_prod,dm0a_s,dm3a_s,dm0a_f,dm3a_f) !! Get the evolution of moments tracers through haze microphysics processes. !! !! The subroutine is a wrapper to the haze microphysics methods. It computes the tendencies !! of moments tracers for production, sedimentation, and coagulation processes for the !! atmospheric column. !! ! Production of the 3rd order moment of the spherical mode distribution from CH4 photolysis (m3.m-3). REAL(kind=mm_wp), INTENT(in), DIMENSION(:) :: m3a_s_prod ! Tendency of the 0th order moment of the spherical mode distribution (m-3). REAL(kind=mm_wp), INTENT(out), DIMENSION(:) :: dm0a_s ! Tendency of the 3rd order moment of the spherical mode distribution (m3.m-3). REAL(kind=mm_wp), INTENT(out), DIMENSION(:) :: dm3a_s ! Tendency of the 0th order moment of the fractal mode distribution (m-3). REAL(kind=mm_wp), INTENT(out), DIMENSION(:) :: dm0a_f ! Tendency of the 3rd order moment of the fractal mode distribution (m3.m-3). REAL(kind=mm_wp), INTENT(out), DIMENSION(:) :: dm3a_f ! Local variables: !~~~~~~~~~~~~~~~~~ REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE :: zdm0as REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE :: zdm3as REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE :: zdm0af REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE :: zdm3af ! Initialization: !~~~~~~~~~~~~~~~~ dm0a_s = 0._mm_wp dm3a_s = 0._mm_wp dm0a_f = 0._mm_wp dm3a_f = 0._mm_wp ALLOCATE(zdm0as(mm_nla),zdm3as(mm_nla),zdm0af(mm_nla),zdm3af(mm_nla)) zdm0as(1:mm_nla) = 0._mm_wp zdm3as(1:mm_nla) = 0._mm_wp zdm0af(1:mm_nla) = 0._mm_wp zdm3af(1:mm_nla) = 0._mm_wp ! Microphysical processes: !~~~~~~~~~~~~~~~~~~~~~~~~~ ! Coagulation IF (mm_w_haze_coag) THEN call mm_haze_coagulation(dm0a_s,dm3a_s,dm0a_f,dm3a_f) ENDIF ! Sedimentation IF (mm_w_haze_sed) THEN call mm_haze_sedimentation(zdm0as,zdm3as,zdm0af,zdm3af) ! Computes precipitations mm_aers_prec = SUM(zdm3as*mm_dzlev*mm_rhoaer) mm_aerf_prec = SUM(zdm3af*mm_dzlev*mm_rhoaer) ! Updates tendencies dm0a_s = dm0a_s + zdm0as dm3a_s = dm3a_s + zdm3as dm0a_f = dm0a_f + zdm0af dm3a_f = dm3a_f + zdm3af ENDIF ! Production IF (mm_w_haze_prod) THEN ! Production by photolysis of CH4 IF (mm_haze_prod_pCH4) THEN dm0a_s = dm0a_s + (m3a_s_prod / (mm_rc_prod**3 * mm_alpha_s(3._mm_wp))) dm3a_s = dm3a_s + m3a_s_prod ELSE call mm_haze_production(zdm0as,zdm3as) dm0a_s = dm0a_s + zdm0as dm3a_s = dm3a_s + zdm3as ENDIF ENDIF RETURN END SUBROUTINE mm_haze_microphysics !============================================================================ ! PRODUCTION PROCESS RELATED METHOD !============================================================================ SUBROUTINE mm_haze_production(dm0s,dm3s) !! Compute the production of aerosols. !! !! @warning !! Spherical aerosols are created at one pressure level. No other source is taken into account. !! This process is controled by two parameters, the pressure level of production and the production rate. !! Then both M0 and M3 of the spherical aerosols distribution are updated in the production zone by addition !! of the production rate along a gaussian shape. !! REAL(kind=mm_wp), INTENT(out), DIMENSION(:) :: dm0s !! Tendency of M0 (m-3). REAL(kind=mm_wp), INTENT(out), DIMENSION(:) :: dm3s !! Tendency of M3 (m3.m-3). ! Local variables: !~~~~~~~~~~~~~~~~~ INTEGER :: i REAL(kind=mm_wp) :: zprod REAL(kind=mm_wp), PARAMETER :: sigz = 20e3_mm_wp, & znorm = dsqrt(2._mm_wp)*sigz ! Locates production altitude !~~~~~~~~~~~~~~~~~~~~~~~~~~~~ zprod = -1._mm_wp DO i=1,mm_nla-1 IF (mm_plev(i) < mm_p_prod.AND.mm_plev(i+1) >= mm_p_prod) THEN zprod = mm_zlay(i) EXIT ENDIF ENDDO ! Sanity check !~~~~~~~~~~~~~ IF (zprod < 0._mm_wp) THEN WRITE(*,'(a)') "cannot find aerosols production altitude" call EXIT(11) ENDIF ! Computes production rate !~~~~~~~~~~~~~~~~~~~~~~~~~ dm3s(:) = mm_tx_prod * (0.75_mm_wp * mm_dt) / (mm_pi * mm_rhoaer) / (2._mm_wp * dsqrt(2._mm_wp) * mm_dzlev(1:mm_nla)) * & (erf((mm_zlev(1:mm_nla)-zprod)/znorm) - erf((mm_zlev(2:mm_nla+1)-zprod)/znorm)) dm0s(:) = dm3s(:) / (mm_rc_prod**3*mm_alpha_s(3._mm_wp)) END SUBROUTINE mm_haze_production !============================================================================ ! SEDIMENTATION PROCESS RELATED METHODS !============================================================================ SUBROUTINE mm_haze_sedimentation(dm0s,dm3s,dm0f,dm3f) !! Interface to sedimentation algorithm. !! !! The subroutine computes the evolution of each moment of the aerosols tracers !! through sedimentation process and returns their tendencies for a timestep. !! ! Tendency of the 0th order moment of the spherical mode (m-3). REAL(kind=mm_wp), INTENT(out), DIMENSION(:) :: dm0s ! Tendency of the 3rd order moment of the spherical mode (m3.m-3). REAL(kind=mm_wp), INTENT(out), DIMENSION(:) :: dm3s ! Tendency of the 0th order moment of the fractal mode (m-3). REAL(kind=mm_wp), INTENT(out), DIMENSION(:) :: dm0f ! Tendency of the 3rd order moment of the fractal mode (m3.m-3). REAL(kind=mm_wp), INTENT(out), DIMENSION(:) :: dm3f ! Local variables: !~~~~~~~~~~~~~~~~~ REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE :: wth ALLOCATE(wth(mm_nle)) ! Force settling velocity to M0 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ IF (mm_wsed_m0) THEN ! Spherical particles call get_weff(0._mm_wp,3._mm_wp,mm_rcs,mm_dt,mm_alpha_s,wth) call let_me_fall_in_peace(mm_m0aer_s,-wth,mm_dt,dm0s) call let_me_fall_in_peace(mm_m3aer_s,-wth,mm_dt,dm3s) ! Get settling velocity and mass flux mm_m0as_vsed(:) = wth(1:mm_nla) mm_m3as_vsed(:) = wth(1:mm_nla) mm_aer_s_flux(:) = (mm_m3aer_s * mm_rhoaer) * wth(2:) ! Fractal particles call get_weff(0._mm_wp,mm_df,mm_rcf,mm_dt,mm_alpha_f,wth) call let_me_fall_in_peace(mm_m0aer_f,-wth,mm_dt,dm0f) call let_me_fall_in_peace(mm_m3aer_f,-wth,mm_dt,dm3f) ! Get settling velocity and mass flux mm_m0af_vsed(:) = wth(1:mm_nla) mm_m3af_vsed(:) = wth(1:mm_nla) mm_aer_f_flux(:) = (mm_m3aer_f * mm_rhoaer) * wth(2:) ! Force settling velocity to M3 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ELSEIF (mm_wsed_m3) THEN ! Spherical particles call get_weff(3._mm_wp,3._mm_wp,mm_rcs,mm_dt,mm_alpha_s,wth) call let_me_fall_in_peace(mm_m3aer_s,-wth,mm_dt,dm3s) call let_me_fall_in_peace(mm_m0aer_s,-wth,mm_dt,dm0s) ! Get settling velocity and mass flux mm_m3as_vsed(:) = wth(1:mm_nla) mm_m0as_vsed(:) = wth(1:mm_nla) mm_aer_s_flux(:) = (mm_m3aer_s * mm_rhoaer) * wth(2:) ! Fractal particles call get_weff(3._mm_wp,mm_df,mm_rcf,mm_dt,mm_alpha_f,wth) call let_me_fall_in_peace(mm_m3aer_f,-wth,mm_dt,dm3f) call let_me_fall_in_peace(mm_m0aer_f,-wth,mm_dt,dm0f) ! Get settling velocity and mass flux mm_m3af_vsed(:) = wth(1:mm_nla) mm_m0af_vsed(:) = wth(1:mm_nla) mm_aer_f_flux(:) = (mm_m3aer_f * mm_rhoaer) * wth(2:) ! No forcing of settling velocity (might be a problem...) !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ELSE ! Spherical particles call get_weff(0._mm_wp,3._mm_wp,mm_rcs,mm_dt,mm_alpha_s,wth) call let_me_fall_in_peace(mm_m0aer_s,-wth,mm_dt,dm0s) ! Get settling velocity mm_m0as_vsed(:) = wth(1:mm_nla) call get_weff(3._mm_wp,3._mm_wp,mm_rcs,mm_dt,mm_alpha_s,wth) call let_me_fall_in_peace(mm_m3aer_s,-wth,mm_dt,dm3s) ! Get settling velocity mm_m3as_vsed(:) = wth(1:mm_nla) ! Get mass flux mm_aer_s_flux(:) = (mm_m3aer_s * mm_rhoaer) * wth(2:) ! Fractal aerosols call get_weff(0._mm_wp,mm_df,mm_rcf,mm_dt,mm_alpha_f,wth) call let_me_fall_in_peace(mm_m0aer_f,-wth,mm_dt,dm0f) ! Get settling velocity mm_m0af_vsed(:) = wth(1:mm_nla) call get_weff(3._mm_wp,mm_df,mm_rcf,mm_dt,mm_alpha_f,wth) call let_me_fall_in_peace(mm_m3aer_f,-wth,mm_dt,dm3f) ! Get settling velocity mm_m3af_vsed(:) = wth(1:mm_nla) ! Get mass flux mm_aer_f_flux(:) = (mm_m3aer_f * mm_rhoaer) * wth(2:) ENDIF DEALLOCATE(wth) RETURN END SUBROUTINE mm_haze_sedimentation SUBROUTINE get_weff(k,df,rc,dt,afun,wth) !! Get the effective settling velocity for aerosols moments. !! !! This method computes the effective settling velocity of the kth order moment of aerosol tracers. !! The basic settling velocity (weff(Mk)) is computed using the following equation: !! Phi_sed (Mk) = \int_0^infty (n(r) . r**k . w(r) . dr) !! = Mk . weff (Mk) !! ! The order of the moment. REAL(kind=mm_wp), INTENT(in) :: k ! Fractal dimension of the aersols. REAL(kind=mm_wp), INTENT(in) :: df ! Characteristic radius associated to the moment at each layer (m). REAL(kind=mm_wp), INTENT(in), DIMENSION(mm_nla) :: rc ! Time step (s). REAL(kind=mm_wp), INTENT(in) :: dt ! Theoretical effective settling velocity at each vertical __levels__ (m.s-1). REAL(kind=mm_wp), INTENT(out), DIMENSION(mm_nle) :: wth ! Inter-moment relation function. INTERFACE FUNCTION afun(order) IMPORT mm_wp REAL(kind=mm_wp), INTENT(in) :: order ! Order of the moment. REAL(kind=mm_wp) :: afun ! Alpha value. END FUNCTION afun END INTERFACE ! Local variables: !~~~~~~~~~~~~~~~~~ INTEGER :: i REAL(kind=mm_wp) :: af1,af2,ar1,ar2 REAL(kind=mm_wp) :: mrcoef, brhoair, thermal_w REAL(kind=mm_wp) :: csto,cslf REAL(kind=mm_wp) :: rb2ra wth(:) = 0._mm_wp ! Molecular's case !~~~~~~~~~~~~~~~~~ mrcoef = 0.74_mm_wp ar1 = (3._mm_wp*df - 6._mm_wp) / df af1 = (df*(k+3._mm_wp) - 6._mm_wp) / df rb2ra = mm_rm**((6._mm_wp-2._mm_wp*df)/df) DO i=1,mm_nle thermal_w = sqrt((8._mm_wp * mm_kboltz * mm_btemp(i)) / (mm_pi * mm_air_mmas)) IF (i == 1) THEN IF (rc(i) <= 0._mm_wp) CYCLE brhoair = mm_rhoair(i) wth(i) = - (mrcoef * mm_effg(mm_zlev(i)) * (mm_rhoaer/brhoair) / thermal_w) * rb2ra * & (afun(af1)/afun(k)) * rc(i)**ar1 ELSEIF (i == mm_nle) THEN IF (rc(i-1) <= 0._mm_wp) CYCLE brhoair = mm_rhoair(mm_nla) + (mm_rhoair(mm_nla) - mm_rhoair(mm_nla-1)) / 2._mm_wp wth(i) = - (mrcoef * mm_effg(mm_zlev(i)) * (mm_rhoaer/brhoair) / thermal_w) * rb2ra * & (afun(af1)/afun(k)) * rc(i-1)**ar1 ELSE IF (rc(i-1) <= 0._mm_wp) CYCLE brhoair = (mm_rhoair(i-1) + mm_rhoair(i)) / 2._mm_wp wth(i) = - (mrcoef * mm_effg(mm_zlev(i)) * (mm_rhoaer/brhoair) / thermal_w) * rb2ra * & (afun(af1)/afun(k)) * rc(i-1)**ar1 ENDIF ENDDO ! Hydrodynamical's case ! In fact: transitory case which is the general case !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ !ar1 = (3._mm_wp*df - 3._mm_wp) / df !ar2 = (3._mm_wp*df - 6._mm_wp) / df ! !af1 = (df*(k+3._mm_wp) - 3._mm_wp) / df !af2 = (df*(k+3._mm_wp) - 6._mm_wp) / df ! !rb2ra = mm_rm**((df-3._mm_wp)/df) ! !DO i=1,mm_nle ! csto = (2._mm_wp*mm_rhoaer*mm_effg(mm_zlev(i))) / (9._mm_wp*mm_eta_air(mm_btemp(i))) ! cslf = mm_akn * mm_lambda_air(mm_btemp(i),mm_plev(i)) / rb2ra ! ! IF (i == 1) THEN ! IF (rc(i) <= 0._mm_wp) CYCLE ! wth(i) = - csto/(rb2ra*afun(k)) * & ! ((afun(af1)*rc(i)**ar1) + (cslf*afun(af2)*rc(i)**ar2)) ! ELSE ! IF (rc(i-1) <= 0._mm_wp) CYCLE ! wth(i) = - csto/(rb2ra*afun(k)) * & ! ((afun(af1)*rc(i-1)**ar1) + (cslf*afun(af2)*rc(i-1)**ar2)) ! ENDIF !ENDDO RETURN END SUBROUTINE get_weff SUBROUTINE let_me_fall_in_peace(mk,ft,dt,dmk) !! Compute the tendency of the moment through sedimentation process. !! !! The method computes the time evolution of the kth order moment through sedimentation: !! dM(k) / dt = Phi(k) / dz !! !! The equation is resolved using a [Crank-Nicolson algorithm](http://en.wikipedia.org/wiki/Crank-Nicolson_method). !! !! Sedimentation algorithm is quite messy. It appeals to the dark side of the Force and uses evil black magic spells !! from ancient times. It is based on \cite{toon1988b,fiadeiro1977,turco1979} and is an update of the algorithm !! originally implemented in the LMDZ-Titan 2D GCM. !! REAL(kind=mm_wp), INTENT(in), DIMENSION(:) :: mk !! kth order moment to sediment (m^k.m^-3). REAL(kind=mm_wp), INTENT(in) :: dt !! Time step (s). REAL(kind=mm_wp), INTENT(in), DIMENSION(:) :: ft !! Downward sedimentation flux (effective velocity of the moment). REAL(kind=mm_wp), INTENT(out), DIMENSION(:) :: dmk !! Tendency of kth order moment (m^k.m^-3). ! Local variables: !~~~~~~~~~~~~~~~~~ INTEGER :: i REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE :: as,bs,cs,mko ALLOCATE(as(mm_nla), bs(mm_nla), cs(mm_nla), mko(mm_nla)) mko(:) = 0._mm_wp cs(1:mm_nla) = ft(2:mm_nla+1) - mm_dzlay(1:mm_nla)/dt IF (ANY(cs > 0._mm_wp)) THEN ! Implicit case as(1:mm_nla) = ft(1:mm_nla) bs(1:mm_nla) = -(ft(2:mm_nle) + mm_dzlay(1:mm_nla)/dt) cs(1:mm_nla) = -mm_dzlay(1:mm_nla) / dt * mk(1:mm_nla) ! (Tri)diagonal matrix inversion mko(1) = cs(1) / bs(1) DO i=2,mm_nla ; mko(i) = (cs(i)-mko(i-1)*as(i))/bs(i) ; ENDDO ELSE ! Explicit case as(1:mm_nla) = -mm_dzlay(1:mm_nla) / dt bs(1:mm_nla) = -ft(1:mm_nla) ! Boundaries mko(1) = cs(1) * mk(1) / as(1) mko(mm_nla) = (bs(mm_nla)*mk(mm_nla-1) + cs(mm_nla)*mk(mm_nla)) / as(mm_nla) ! Interior mko(2:mm_nla-1) = (bs(2:mm_nla-1)*mk(1:mm_nla-2) + cs(2:mm_nla-1)*mk(2:mm_nla-1)) & /as(2:mm_nla-1) ENDIF dmk = mko - mk DEALLOCATE(as,bs,cs,mko) RETURN END SUBROUTINE let_me_fall_in_peace !============================================================================ ! COAGULATION PROCESS RELATED METHODS !============================================================================ SUBROUTINE mm_haze_coagulation(dM0s,dM3s,dM0f,dM3f) !! Get the evolution of the aerosols moments vertical column due to coagulation process. !! !! This is the main method of the coagulation process: !! 1. Computes gamma pre-factor for each parts of the coagulation equation(s). !! 2. Applies the electic correction on the gamma pre-factor. !! 3. Computes the specific flow regime kernels (molecular or hydrodynamic). !! 4. Computes the harmonic mean of the kernels (transitory regime). !! 5. Finally computes the tendencies of the moments. !! !! @Warning !! If the transfert probabilities are set to 1 for the two flow regimes (pco and pfm), !! a floating point exception occured (i.e. a NaN) as we perform a division by zero !! ! Tendency of the 0th order moment of the spherical size-distribution over a time step (m-3). REAL(kind=mm_wp), INTENT(out), DIMENSION(:) :: dM0s ! Tendency of the 3rd order moment of the spherical size-distribution (m3.m-3). REAL(kind=mm_wp), INTENT(out), DIMENSION(:) :: dM3s ! Tendency of the 0th order moment of the fractal size-distribution over a time step (m-3). REAL(kind=mm_wp), INTENT(out), DIMENSION(:) :: dM0f ! Tendency of the 3rd order moment of the fractal size-distribution over a time step (m3.m-3). REAL(kind=mm_wp), INTENT(out), DIMENSION(:) :: dM3f ! Local variables: !~~~~~~~~~~~~~~~~~ INTEGER :: i REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE :: kco,kfm,slf REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE :: g_co,g_fm,pco,pfm,mq,tmp REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE :: a_ss,a_sf,b_ss,b_ff,c_ss,c_sf ! Sanity check: IF (mm_coag_choice < 0 .OR. mm_coag_choice > mm_coag_ss+mm_coag_sf+mm_coag_ff) & STOP "invalid choice for coagulation mode interaction activation" ! Allocates local arrays: ALLOCATE(kco(mm_nla),kfm(mm_nla),slf(mm_nla)) ALLOCATE(g_co(mm_nla),g_fm(mm_nla), & pco(mm_nla),pfm(mm_nla), & mq(mm_nla), tmp(mm_nla)) ALLOCATE(a_ss(mm_nla),a_sf(mm_nla), & b_ss(mm_nla),b_ff(mm_nla), & c_ss(mm_nla),c_sf(mm_nla)) a_ss(:) = 0._mm_wp ; a_sf(:) = 0._mm_wp b_ss(:) = 0._mm_wp ; b_ff(:) = 0._mm_wp c_ss(:) = 0._mm_wp ; c_sf(:) = 0._mm_wp ! 1. Gamma pre-factor: !~~~~~~~~~~~~~~~~~~~~~ ! gets kco, kfm pre-factors kco(:) = mm_get_kco(mm_temp) ; kfm(:) = mm_get_kfm(mm_temp) ! get slf (slip-flow factor) slf(:) = mm_akn * mm_lambda_air(mm_temp,mm_play) ! 2,3,4. Electic correction, Flow regime kernels, Harmonic mean: !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ DO i=1,mm_nla ! SS interactions !~~~~~~~~~~~~~~~~ IF (mm_rcs(i) > mm_rcs_min .AND. IAND(mm_coag_choice,mm_coag_ss) /= 0) THEN ! Probability S --> S for M0/CO and M0/FM pco(i) = mm_ps2s(mm_rcs(i),0,0) pfm(i) = mm_ps2s(mm_rcs(i),0,1) ! Gamma (kernel dependent) g_co(i) = g0ssco(mm_rcs(i),slf(i),kco(i)) g_fm(i) = g0ssfm(mm_rcs(i),kfm(i)) ! Molecular's case !~~~~~~~~~~~~~~~~~ a_ss(i) = (pfm(i) - 2._mm_wp) * g_fm(i) b_ss(i) = (1._mm_wp - pfm(i)) * g_fm(i) ! Hydrodynamical's case ! In fact: transitory case which is the general case !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! (A_SS_CO x A_SS_FM) / (A_SS_CO + A_SS_FM) !IF (g_co(i)*(pco(i)-2._mm_wp)+g_fm(i)*(pfm(i)-2._mm_wp) /= 0._mm_wp) THEN ! a_ss(i) = (g_co(i)*(pco(i)-2._mm_wp)*g_fm(i)*(pfm(i)-2._mm_wp))/(g_co(i)*(pco(i)-2._mm_wp)+g_fm(i)*(pfm(i)-2._mm_wp)) !ENDIF ! (B_SS_CO x B_SS_FM) / (B_SS_CO + B_SS_FM) !IF (g_co(i)*(1._mm_wp-pco(i))+g_fm(i)*(1._mm_wp-pfm(i)) /= 0._mm_wp) THEN ! b_ss(i) = (g_co(i)*(1._mm_wp-pco(i))*g_fm(i)*(1._mm_wp-pfm(i)))/(g_co(i)*(1._mm_wp-pco(i))+g_fm(i)*(1._mm_wp-pfm(i))) !ENDIF ! Eletric charge correction for M0/SS interactions mq(i) = mm_qmean(mm_rcs(i),mm_rcs(i),0,'SS',mm_temp(i)) a_ss(i) = a_ss(i) * mq(i) b_ss(i) = b_ss(i) * mq(i) ! Probability S --> S for M3/CO and M3/FM pco(i) = mm_ps2s(mm_rcs(i),3,0) pfm(i) = mm_ps2s(mm_rcs(i),3,1) ! Gamma (kernel dependent) g_co(i) = g3ssco(mm_rcs(i),slf(i),kco(i)) * (pco(i)-1._mm_wp) g_fm(i) = g3ssfm(mm_rcs(i),kfm(i)) * (pfm(i)-1._mm_wp) ! Molecular's case !~~~~~~~~~~~~~~~~~ c_ss(i) = g_fm(i) ! Hydrodynamical's case ! In fact: transitory case which is the general case !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! (C_SS_CO x C_SS_FM) / (C_SS_CO + C_SS_FM) !IF (g_co(i) + g_fm(i) /= 0._mm_wp) THEN ! c_ss(i) = (g_co(i) * g_fm(i)) / (g_co(i) + g_fm(i)) !ENDIF IF (b_ss(i) <= 0._mm_wp) c_ss(i) = 0._mm_wp ! Eletric charge correction for M3/SS interactions mq(i) = mm_qmean(mm_rcs(i),mm_rcs(i),3,'SS',mm_temp(i)) c_ss(i) = c_ss(i) * mq(i) ENDIF ! end SS interactions ! SF interactions !~~~~~~~~~~~~~~~~ IF (mm_rcs(i) > mm_rcs_min .AND. mm_rcf(i) > mm_rcf_min .AND. IAND(mm_coag_choice,mm_coag_sf) /= 0) THEN ! Gamma (kernel dependent) g_co(i) = g0sfco(mm_rcs(i),mm_rcf(i),slf(i),kco(i)) g_fm(i) = g0sffm(mm_rcs(i),mm_rcf(i),kfm(i)) ! Molecular's case !~~~~~~~~~~~~~~~~~ a_sf(i) = g_fm(i) ! Hydrodynamical's case ! In fact: transitory case which is the general case !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! (A_SF_CO x A_SF_FM) / (A_SF_CO + A_SF_FM) !IF(g_co(i) + g_fm(i) /= 0._mm_wp) THEN ! a_sf(i) = (g_co(i) * g_fm(i)) / (g_co(i) + g_fm(i)) !ENDIF ! Eletric charge correction for M0/SF interactions mq(i) = mm_qmean(mm_rcs(i),mm_rcf(i),0,'SF',mm_temp(i)) a_sf(i) = a_sf(i) * mq(i) ! Gamma (kernel dependent) g_co(i) = g3sfco(mm_rcs(i),mm_rcf(i),slf(i),kco(i)) g_fm(i) = g3sffm(mm_rcs(i),mm_rcf(i),kfm(i)) ! Molecular's case !~~~~~~~~~~~~~~~~~ c_sf(i) = g_fm(i) ! Hydrodynamical's case ! In fact: transitory case which is the general case !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! (C_SF_CO x C_SF_FM) / (C_SF_CO + C_SF_FM) !IF (g_co(i) + g_fm(i) /= 0._mm_wp) THEN ! c_sf(i) = (g_co(i) * g_fm(i)) / (g_co(i) + g_fm(i)) !ENDIF ! Eletric charge correction for M3/SF interactions mq(i) = mm_qmean(mm_rcs(i),mm_rcf(i),3,'SF',mm_temp(i)) c_sf(i) = c_sf(i) * mq(i) ENDIF ! end SF interactions ! FF interactions !~~~~~~~~~~~~~~~~ IF(mm_rcf(i) > mm_rcf_min .AND. IAND(mm_coag_choice,mm_coag_sf) /= 0) THEN ! Gamma (kernel dependent) g_co(i) = g0ffco(mm_rcf(i),slf(i),kco(i)) g_fm(i) = g0fffm(mm_rcf(i),kfm(i)) ! Molecular's case !~~~~~~~~~~~~~~~~~ b_ff(i) = g_fm(i) ! Hydrodynamical's case ! In fact: transitory case which is the general case !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! (B_FF_CO x B_FF_FM) / (B_FF_CO + B_FF_FM) !IF (g_co(i) + g_fm(i) /= 0._mm_wp) THEN ! b_ff(i) = (g_co(i) * g_fm(i)) / (g_co(i) + g_fm(i)) !ENDIF ! Eletric charge correction for M0/FF interactions mq(i) = mm_qmean(mm_rcf(i),mm_rcf(i),0,'FF',mm_temp(i)) b_ff(i) = b_ff(i) * mq(i) ENDIF ! end FF interactions ENDDO ! end n_lay DEALLOCATE(g_co,g_fm,kco,kfm,pco,pfm,slf) ! 5. Tendencies of the moments: !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! Now we will use the kharm two by two to compute : ! dM_0_S/mm_dt = kharm(1) * M_0_S^2 - kharm(2) * M_0_S * m_0_F ! dM_0_F/mm_dt = kharm(3) * M_0_S^2 - kharm(4) * M_0_F^2 ! dM_3_S/mm_dt = kharm(5) * M_3_S^2 - kharm(6) * M_3_S * m_3_F ! dM_3_F/mm_dt = - dM_3_S/mm_dt ! ! We use a (semi) implicit scheme: ! When X appears as square we set one X at t+1, the other a t. ! --- dM0s tmp(:) = mm_dt * (a_ss*mm_m0aer_s - a_sf*mm_m0aer_f) dm0s(:) = mm_m0aer_s * (tmp / (1._mm_wp - tmp)) ! --- dM0f tmp(:) = mm_dt * b_ff * mm_m0aer_f dm0f(:) = (b_ss*mm_dt*mm_m0aer_s**2 - tmp*mm_m0aer_f) / (1._mm_wp + tmp) ! --- dM3s tmp(:) = mm_dt * (c_ss*mm_m3aer_s - c_sf*mm_m3aer_f) dm3s(:) = mm_m3aer_s * (tmp / (1._mm_wp - tmp)) ! --- dM3f dm3f(:) = -dm3s DEALLOCATE(a_ss,a_sf,b_ss,b_ff,c_ss,c_sf,tmp) RETURN END SUBROUTINE mm_haze_coagulation ! ===================== ! Hydrodynamical's case ! ===================== ELEMENTAL FUNCTION g0ssco(rcs,slf,kco) RESULT(res) !! Get gamma pre-factor for the 0th order moment with SS interactions in the Continuous flow regime. !! REAL(kind=mm_wp), INTENT(in) :: rcs !! Characteristic radius of the spherical size-distribution. REAL(kind=mm_wp), INTENT(in) :: slf !! Slip-Flow correction pre-factor. REAL(kind=mm_wp), INTENT(in) :: kco !! Thermodynamic continuous flow regime pre-factor. REAL(kind=mm_wp) :: res !! Gamma coagulation kernel pre-factor. REAL(kind=mm_wp) :: a1, a2, a3 res = 0._mm_wp ; IF (rcs <= 0._mm_wp) RETURN ! Computes mm_alpha coefficients a1 = mm_alpha_s(1._mm_wp) a2 = mm_alpha_s(-1._mm_wp) a3 = mm_alpha_s(-2._mm_wp) ! Computes gamma pre-factor res = kco * (1._mm_wp + a1*a2 + slf/rcs*(a2+a1*a3)) RETURN END FUNCTION g0ssco ELEMENTAL FUNCTION g0sfco(rcs,rcf,slf,kco) RESULT(res) !! Get gamma pre-factor for the 0th order moment with SF interactions in the Continuous flow regime. !! REAL(kind=mm_wp), INTENT(in) :: rcs !! Characteristic radius of the spherical size-distribution. REAL(kind=mm_wp), INTENT(in) :: rcf !! Characteristic radius of the fractal size-distribution. REAL(kind=mm_wp), INTENT(in) :: slf !! Slip-Flow correction pre-factor. REAL(kind=mm_wp), INTENT(in) :: kco !! Thermodynamic continuous flow regime pre-factor. REAL(kind=mm_wp) :: res !! Gamma coagulation kernel pre-factor. REAL(kind=mm_wp) :: a1, a2, a3, a4, a5, a6, e, rcff res = 0._mm_wp ; IF (rcs <= 0._mm_wp .OR. rcf <= 0._mm_wp) RETURN ! Computes mm_alpha and other coefficients e = 3._mm_wp/mm_df rcff = rcf**e * mm_rb2ra a1 = mm_alpha_s(1._mm_wp) ; a2 = mm_alpha_f(-e) ; a3 = mm_alpha_f(e) a4 = mm_alpha_s(-1._mm_wp) ; a5 = mm_alpha_s(-2._mm_wp) ; a6 = mm_alpha_f(-2._mm_wp*e) ! Computes gamma pre-factor res = kco * (2._mm_wp + a1*a2*rcs/rcff + & a4*a3*rcff/rcs + & slf * (a4/rcs + a2/rcff + & a5*a3*rcff/rcs**2 + & a1*a6*rcs/rcff**2)) RETURN END FUNCTION g0sfco ELEMENTAL FUNCTION g0ffco(rcf,slf,kco) RESULT(res) !! Get gamma pre-factor for the 0th order moment with FF interactions in the Continuous flow regime. !! REAL(kind=mm_wp), INTENT(in) :: rcf !! Characteristic radius of the fractal size-distribution. REAL(kind=mm_wp), INTENT(in) :: slf !! Slip-Flow correction pre-factor. REAL(kind=mm_wp), INTENT(in) :: kco !! Thermodynamic continuous flow regime pre-factor. REAL(kind=mm_wp) :: res !! Gamma coagulation kernel pre-factor. REAL(kind=mm_wp) :: a1, a2, a3, e, rcff res = 0._mm_wp ; IF (rcf <= 0._mm_wp) RETURN ! Computes mm_alpha and other coefficients e = 3._mm_wp/mm_df rcff = rcf**e * mm_rb2ra a1 = mm_alpha_f(e) ; a2 = mm_alpha_f(-e) ; a3 = mm_alpha_f(-2._mm_wp*e) ! Computes gamma pre-factor res = kco * (1._mm_wp + a1*a2 + slf/rcff *(a2+a1*a3)) RETURN END FUNCTION g0ffco ELEMENTAL FUNCTION g3ssco(rcs, slf, kco) RESULT(res) !! Get gamma pre-factor for the 3rd order moment with SS interactions in the Continuous flow regime. !! REAL(kind=mm_wp), INTENT(in) :: rcs !! Characteristic radius of the spherical size-distribution. REAL(kind=mm_wp), INTENT(in) :: slf !! Slip-Flow correction pre-factor. REAL(kind=mm_wp), INTENT(in) :: kco !! Thermodynamic continuous flow regime pre-factor. REAL(kind=mm_wp) :: res !! Gamma coagulation kernel pre-factor. REAL(kind=mm_wp) :: a1, a2, a3, a4, a5, a6 res = 0._mm_wp ; IF (rcs <= 0._mm_wp) RETURN ! Computes mm_alpha coefficients a1 = mm_alpha_s(3._mm_wp) ; a2 = mm_alpha_s(2._mm_wp) ; a3 = mm_alpha_s(1._mm_wp) a4 = mm_alpha_s(4._mm_wp) ; a5 = mm_alpha_s(-1._mm_wp) ; a6 = mm_alpha_s(-2._mm_wp) ! Computes gamma pre-factor res = kco/(a1**2*rcs**3) * (2._mm_wp*a1 + a2*a3 + a4*a5 + & slf/rcs*(a3**2 + a4*a6 + a1*a5 + a2)) RETURN END FUNCTION g3ssco ELEMENTAL FUNCTION g3sfco(rcs, rcf, slf, kco) RESULT(res) !! Get gamma pre-factor for the 3rd order moment with SF interactions in the Continuous flow regime. !! REAL(kind=mm_wp), INTENT(in) :: rcs !! Characteristic radius of the spherical size-distribution. REAL(kind=mm_wp), INTENT(in) :: rcf !! Characteristic radius of the fractal size-distribution. REAL(kind=mm_wp), INTENT(in) :: slf !! Slip-Flow correction pre-factor. REAL(kind=mm_wp), INTENT(in) :: kco !! Thermodynamic continuous flow regime pre-factor. REAL(kind=mm_wp) :: res !! Gamma coagulation kernel pre-factor. REAL(kind=mm_wp) :: a1, a2, a3, a4, a5, a6, a7, a8, e, rcff res = 0._mm_wp ; IF (rcs <= 0._mm_wp .OR. rcf <= 0._mm_wp) RETURN ! Computes mm_alpha and other coefficients e = 3._mm_wp/mm_df rcff = rcf**e * mm_rb2ra a1 = mm_alpha_s(3._mm_wp) ; a2 = mm_alpha_s(4._mm_wp) ; a3 = mm_alpha_f(-e) a4 = mm_alpha_s(2._mm_wp) ; a5 = mm_alpha_f(e) ; a6 = mm_alpha_s(1._mm_wp) a7 = mm_alpha_f(-2._mm_wp*e) ; a8 = mm_alpha_f(3._mm_wp) ! Computes gamma pre-factor res = kco/(a1*a8*(rcs*rcf)**3) * (2._mm_wp*a1*rcs**3 + a2*rcs**4*a3/rcff + a4*rcs**2*a5*rcff + & slf * (a4*rcs**2 + a1*rcs**3*a3/rcff + & a6*rcs*a5*rcff + & a2*rcs**4*a7/rcff**2)) RETURN END FUNCTION g3sfco ! ================ ! Molecular's case ! ================ ELEMENTAL FUNCTION g0ssfm(rcs, kfm) RESULT(res) !! Get gamma pre-factor for the 0th order moment with SS interactions in the Free Molecular flow regime. !! REAL(kind=mm_wp), INTENT(in) :: rcs !! Characteristic radius of the spherical size-distribution. REAL(kind=mm_wp), INTENT(in) :: kfm !! Thermodynamic free molecular flow regime pre-factor. REAL(kind=mm_wp) :: res !! Gamma coagulation kernel pre-factor. REAL(kind=mm_wp) :: a1, a2, a3, a4, a5 res = 0._mm_wp ; IF (rcs <= 0._mm_wp) RETURN ! Computes mm_alpha coefficients a1 = mm_alpha_s(0.5_mm_wp) ; a2 = mm_alpha_s(1._mm_wp) ; a3 = mm_alpha_s(-0.5_mm_wp) a4 = mm_alpha_s(2._mm_wp) ; a5 = mm_alpha_s(-1.5_mm_wp) ! Computes gamma pre-factor res = kfm * (rcs**0.5_mm_wp*mm_get_btk(1,0)) * (a1 + 2._mm_wp*a2*a3 + a4*a5) RETURN END FUNCTION g0ssfm ELEMENTAL FUNCTION g0sffm(rcs, rcf, kfm) RESULT(res) !! Get gamma pre-factor for the 0th order moment with SF interactions in the Free Molecular flow regime. !! REAL(kind=mm_wp), INTENT(in) :: rcs !! Characteristic radius of the spherical size-distribution. REAL(kind=mm_wp), INTENT(in) :: rcf !! Characteristic radius of the fractal size-distribution. REAL(kind=mm_wp), INTENT(in) :: kfm !! Thermodynamic free molecular flow regime pre-factor. REAL(kind=mm_wp) :: res !! Gamma coagulation kernel pre-factor. REAL(kind=mm_wp) :: a1, a2, a3, a4, a5, a6, a7, a8, a9, a10 REAL(kind=mm_wp) :: e1, e2, e3, e4 REAL(kind=mm_wp) :: rcff1, rcff2, rcff3, rcff4 res = 0._mm_wp ; IF (rcs <= 0._mm_wp .OR. rcf <= 0._mm_wp) RETURN ! Computes mm_alpha and other coefficients e1 = 3._mm_wp/mm_df e2 = 6._mm_wp/mm_df e3 = (6._mm_wp-3._mm_wp*mm_df)/(2._mm_wp*mm_df) e4 = (12._mm_wp-3._mm_wp*mm_df)/(2._mm_wp*mm_df) rcff1 = rcf**e1 * mm_rb2ra rcff2 = rcff1**2 rcff3 = rcf**e3 * mm_rb2ra rcff4 = rcf**e4 * mm_rb2ra**2 a1 = mm_alpha_s(0.5_mm_wp) ; a2 = mm_alpha_s(-0.5_mm_wp) ; a3 = mm_alpha_f(e1) a4 = mm_alpha_s(-1.5_mm_wp) ; a5 = mm_alpha_f(e2) ; a6 = mm_alpha_s(2._mm_wp) a7 = mm_alpha_f(-1.5_mm_wp) ; a8 = mm_alpha_s(1._mm_wp) ; a9 = mm_alpha_f(e3) a10 = mm_alpha_f(e4) ! Computes gamma pre-factor res = kfm * mm_get_btk(4,0) * (a1*rcs**0.5_mm_wp + 2._mm_wp*a2*a3*rcff1/rcs**0.5_mm_wp + & a4*a5*rcff2/rcs**1.5_mm_wp + & a6*a7*rcs**2/rcf**1.5_mm_wp + & 2._mm_wp*a8*a9*rcs*rcff3 + & a10*rcff4) RETURN END FUNCTION g0sffm ELEMENTAL FUNCTION g0fffm(rcf, kfm) RESULT(res) !! Get gamma pre-factor for the 0th order moment with FF interactions in the Free Molecular flow regime. !! REAL(kind=mm_wp), INTENT(in) :: rcf !! Characteristic radius of the fractal size-distribution. REAL(kind=mm_wp), INTENT(in) :: kfm !! Thermodynamic free molecular flow regime pre-factor. REAL(kind=mm_wp) :: res !! Gamma coagulation kernel pre-factor. REAL(kind=mm_wp) :: a1, a2, a3, a4, a5, e1, e2, e3, rcff res = 0._mm_wp ; IF (rcf <= 0._mm_wp) RETURN ! Computes mm_alpha and other coefficients e1 = 3._mm_wp/mm_df e2 = (6_mm_wp-3._mm_wp*mm_df)/(2._mm_wp*mm_df) e3 = (12._mm_wp-3._mm_wp*mm_df)/(2._mm_wp*mm_df) rcff = rcf**e3 * mm_rb2ra**2 a1 = mm_alpha_f(e3) ; a2 = mm_alpha_f(e1) ; a3 = mm_alpha_f(e2) a4 = mm_alpha_f(-1.5_mm_wp) ; a5 = mm_alpha_f(2._mm_wp*e1) ! Computes gamma pre-factor res = kfm * (rcff*mm_get_btk(3,0)) * (a1 + 2._mm_wp*a2*a3 + a4*a5) RETURN END FUNCTION g0fffm ELEMENTAL FUNCTION g3ssfm(rcs, kfm) RESULT(res) !! Get gamma pre-factor for the 3rd order moment with SS interactions in the Free Molecular flow regime. !! REAL(kind=mm_wp), INTENT(in) :: rcs !! Characteristic radius of the spherical size-distribution. REAL(kind=mm_wp), INTENT(in) :: kfm !! Thermodynamic free molecular flow regime pre-factor. REAL(kind=mm_wp) :: res !! Gamma coagulation kernel pre-factor. REAL(kind=mm_wp) :: a1, a2, a3, a4, a5, a6, a7, a8, a9, a10, a11 res = 0._mm_wp ; IF (rcs <= 0._mm_wp) RETURN ! Computes mm_alpha coefficients a1 = mm_alpha_s(3.5_mm_wp) ; a2 = mm_alpha_s(1._mm_wp) ; a3 = mm_alpha_s(2.5_mm_wp) a4 = mm_alpha_s(2._mm_wp) ; a5 = mm_alpha_s(1.5_mm_wp) ; a6 = mm_alpha_s(5._mm_wp) a7 = mm_alpha_s(-1.5_mm_wp) ; a8 = mm_alpha_s(4._mm_wp) ; a9 = mm_alpha_s(-0.5_mm_wp) a10 = mm_alpha_s(3._mm_wp) ; a11 = mm_alpha_s(0.5_mm_wp) ! Computes gamma pre-factor res = kfm * (mm_get_btk(1,3)/(a10**2*rcs**2.5_mm_wp)) * (a1 + 2._mm_wp*a2*a3 + & a4*a5 + a6*a7 + & 2._mm_wp*a8*a9 + a10*a11) RETURN END FUNCTION g3ssfm ELEMENTAL FUNCTION g3sffm(rcs, rcf, kfm) RESULT(res) !! Get gamma pre-factor for the 3rd order moment with SF interactions in the Free Molecular flow regime. !! REAL(kind=mm_wp), INTENT(in) :: rcs !! Characteristic radius of the spherical size-distribution. REAL(kind=mm_wp), INTENT(in) :: rcf !! Characteristic radius of the fractal size-distribution. REAL(kind=mm_wp), INTENT(in) :: kfm !! Thermodynamic free molecular flow regime pre-factor. REAL(kind=mm_wp) :: res !! Gamma coagulation kernel pre-factor. REAL(kind=mm_wp) :: a1, a2, a3, a4, a5, a6, a7, a8, a9, a10, a11, a12 REAL(kind=mm_wp) :: e1, e2, e3, rcff1, rcff2, rcff3 res = 0._mm_wp ; IF (rcs <= 0._mm_wp .OR. rcf <= 0._mm_wp) RETURN ! Computes mm_alpha and other coefficients e1 = 3._mm_wp/mm_df e2 = (6._mm_wp-3._mm_wp*mm_df)/(2._mm_wp*mm_df) e3 = (12._mm_wp-3._mm_wp*mm_df)/(2._mm_wp*mm_df) rcff1 = rcf**e1 * mm_rb2ra rcff2 = rcf**e2 * mm_rb2ra rcff3 = rcf**e3 * mm_rb2ra**2 a1 = mm_alpha_s(3.5_mm_wp) ; a2 = mm_alpha_s(2.5_mm_wp) ; a3 = mm_alpha_f(e1) a4 = mm_alpha_s(1.5_mm_wp) ; a5 = mm_alpha_f(2._mm_wp*e1) ; a6 = mm_alpha_s(5._mm_wp) a7 = mm_alpha_f(-1.5_mm_wp) ; a8 = mm_alpha_s(4._mm_wp) ; a9 = mm_alpha_f(e2) a10 = mm_alpha_s(3._mm_wp) ; a11 = mm_alpha_f(e3) ; a12 = mm_alpha_f(3._mm_wp) ! Computes gamma pre-factor res = kfm * (mm_get_btk(4,3)/(a10*a12*(rcs*rcf)**3)) * (a1*rcs**3.5_mm_wp + & 2._mm_wp*a2*a3*rcs**2.5_mm_wp*rcff1 + & a4*a5*rcs**1.5_mm_wp*rcff1**2 + & a6*a7*rcs**5/rcf**1.5_mm_wp + & 2._mm_wp*a8*a9*rcs**4*rcff2 + & a10*a11*rcs**3*rcff3) RETURN END FUNCTION g3sffm END MODULE MP2M_HAZE