MODULE MP2M_CLOUDS !============================================================================ ! ! Purpose ! ------- ! Clouds microphysics module. ! ! This module contains all definitions of the microphysics processes related to ! clouds (nucleation, condensation, sedimentation). ! ! The interface methods always use the global variables defined in [[mm_globals(module)]] when values ! (temperature, pressure, moments...) over the vertical grid are required. ! Consequently, all these functions only deal with output arguments which are most of the time the ! tendencies of relevant variables on the atmospheric column. ! ! @Warning ! The tendencies returned by the method are always defined over the vertical grid ! from __TOP__ to __GROUND__. ! ! The module also contains ten methods: ! - mm_cloud_microphysics | Evolution of moments tracers through clouds microphysics processes. ! - mm_cloud_nucond | Get moments tendencies through nucleation/condensation/evaporation. ! - nc_esp | Get moments tendencies through nucleation/condensation/evaporation of a given condensible species. ! - hetnucl_rate_aer | Get heterogeneous nucleation rate on aerosols. ! - growth_rate | Get growth rate through condensation/evaporation process. ! - mm_cloud_sedimentation | Compute the tendency of clouds related moments through sedimentation process. ! - exchange | Compute the exchange matrix. ! - getnzs | Compute displacement of a cell under sedimentation process. ! - wsettle | Compute the settling velocity of a spherical particle. ! - get_mass_flux | Get the mass flux of clouds related moment through sedimention. ! ! Authors ! ------- ! B. de Batz de Trenquelléon (10/2025) ! !============================================================================ USE MP2M_MPREC USE MP2M_GLOBALS USE MP2M_METHODS USE MP2M_CLOUDS_METHODS IMPLICIT NONE PRIVATE PUBLIC :: mm_cloud_microphysics, mm_cloud_nucond, mm_cloud_sedimentation CONTAINS !============================================================================ ! CLOUDS MICROPHYSICS INTERFACE SUBROUTINE !============================================================================ SUBROUTINE mm_cloud_microphysics(Hdm0as,Hdm3as,Hdm0af,Hdm3af,& Cdm0as,Cdm3as,Cdm0af,Cdm3af,& Cdm0ccn,Cdm3ccn,Cdm3ice,Cdmugas) !! Get the evolution of moments tracers through cloud microphysics processes. !! The subroutine is a wrapper to the cloud microphysics methods. It computes the tendencies of moments !! tracers for nucleation, condensation and sedimentation processes for the atmospheric column. !! !! @note !! Both __dm3ice__ and __dmugas__ are 2D-array with the vertical layers in first dimension and the number !! of ice components in the second. !! ! Tendency of the 0th and 3rd order moment of the aerosols (spherical mode) through haze microphysics processes. REAL(kind=mm_wp), DIMENSION(:), INTENT(in) :: Hdm0as ! (m-3) REAL(kind=mm_wp), DIMENSION(:), INTENT(in) :: Hdm3as ! (m3.m-3) ! Tendency of the 0th and 3rd order moment of the aerosols (fractal mode) through haze microphysics processes. REAL(kind=mm_wp), DIMENSION(:), INTENT(in) :: Hdm0af ! (m-3) REAL(kind=mm_wp), DIMENSION(:), INTENT(in) :: Hdm3af ! (m3.m-3) ! Tendency of the 0th and 3rd order moment of the aerosols (spherical mode) through cloud microphysics processes. REAL(kind=mm_wp), DIMENSION(:), INTENT(inout) :: Cdm0as ! (m-3) REAL(kind=mm_wp), DIMENSION(:), INTENT(inout) :: Cdm3as ! (m3.m-3) ! Tendency of the 0th and 3rd order moment of the aerosols (fractal mode) through cloud microphysics processes. REAL(kind=mm_wp), DIMENSION(:), INTENT(inout) :: Cdm0af ! (m-3) REAL(kind=mm_wp), DIMENSION(:), INTENT(inout) :: Cdm3af ! (m3.m-3) ! Tendency of the 0th and 3rd order moment of the ccn distribution through cloud microphysics processes. REAL(kind=mm_wp), DIMENSION(:), INTENT(inout) :: Cdm0ccn ! (m-3) REAL(kind=mm_wp), DIMENSION(:), INTENT(inout) :: Cdm3ccn ! (m3.m-3) ! Tendencies of the 3rd order moments of each ice components through cloud microphysics processes. REAL(kind=mm_wp), DIMENSION(:,:), INTENT(inout) :: Cdm3ice ! (m3.m-3) ! Tendencies of each condensible gaz species through cloud microphysics processes. REAL(kind=mm_wp), DIMENSION(:,:), INTENT(inout) :: Cdmugas ! (mol.mol-1) ! Local variables !~~~~~~~~~~~~~~~~ INTEGER :: i ! Temporary tendencies through sedimentation (usefull for diagnostics) REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE :: zCdm0ccn,zCdm3ccn REAL(kind=mm_wp), DIMENSION(:,:), ALLOCATABLE :: zCdm3ice ! Initialization: !~~~~~~~~~~~~~~~~ ALLOCATE(zCdm0ccn(mm_nla),zCdm3ccn(mm_nla),zCdm3ice(mm_nla,mm_nesp)) Cdm0as(:) = 0._mm_wp ; Cdm3as(:) = 0._mm_wp Cdm0af(:) = 0._mm_wp ; Cdm3af(:) = 0._mm_wp Cdm0ccn(:) = 0._mm_wp ; Cdm3ccn(:) = 0._mm_wp Cdm3ice(:,:) = 0._mm_wp ; Cdmugas(:,:) = 0._mm_wp zCdm0ccn(:) = 0._mm_wp ; zCdm3ccn(:) = 0._mm_wp zCdm3ice(:,:) = 0._mm_wp ! Calls nucleation/condensation: ! And update saturation ratio, nucleation rate and growth rate diagnostic. !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ call mm_cloud_nucond(Hdm0as,Hdm3as,Hdm0af,Hdm3af,Cdm0as,Cdm3as,Cdm0af,Cdm3af,& Cdm0ccn,Cdm3ccn,Cdm3ice,Cdmugas,mm_gas_sat,mm_nrate,mm_grate) ! Calls sedimentation: !~~~~~~~~~~~~~~~~~~~~~ call mm_cloud_sedimentation(zCdm0ccn,zCdm3ccn,zCdm3ice) ! Computes diagnostics: !~~~~~~~~~~~~~~~~~~~~~~ ! Settling velocity [m.s-1] of clouds (ccn and ices) mm_cld_vsed(:) = wsettle(mm_play,mm_temp,mm_zlay,mm_rhoair,mm_drho,mm_drad) ! Flux [kg.m-2.s-1] and precipitation [kg.m-2] of ccn mm_ccn_flux(:) = get_mass_flux(mm_rhoaer,mm_m3ccn(:)) mm_ccn_prec = SUM(zCdm3ccn(:)*mm_dzlev*mm_rhoaer) ! Flux [kg.m-2.s-1] and precipitation [kg.m-2] of ices DO i = 1, mm_nesp mm_ice_fluxes(:,i) = get_mass_flux(mm_xESPS(i)%rho_s,mm_m3ice(:,i)) mm_ice_prec(i) = SUM(zCdm3ice(:,i)*mm_dzlev*mm_xESPS(i)%rho_s) ENDDO ! Updates tendencies: !~~~~~~~~~~~~~~~~~~~~ Cdm0ccn = Cdm0ccn + zCdm0ccn Cdm3ccn = Cdm3ccn + zCdm3ccn Cdm3ice = Cdm3ice + zCdm3ice END SUBROUTINE mm_cloud_microphysics !============================================================================ ! NUCLEATION/CONDENSATION PROCESS RELATED METHODS !============================================================================ SUBROUTINE mm_cloud_nucond(Hdm0as,Hdm3as,Hdm0af,Hdm3af,Cdm0as,Cdm3as,Cdm0af,Cdm3af,& Cdm0ccn,Cdm3ccn,Cdm3ice,Cdmugas,gassat,nrate,grate) !! Get moments tendencies through nucleation and condensation/evaporation. !! !! The method is a wrapper of [[mm_clouds(module):nc_esp(subroutine)]] which computes the !! tendencies of tracers for all the condensible species given in the vector __xESPS__. !! !! Aerosols and CCN distribution evolution depends on the ice components: !! - For nucleation only creation of CCN can occur. !! - For condensation/evaporation only loss of CCN can occur. !! !! @note !! We use the simple following rule to compute the variation of CCN and aerosols: !! The global variation of CCN (and thus aerosols) is determined from the most intense activity !! of the ice components. !! !! @warning !! __xESPS__, __m3i__ and __gazes__ must share the same indexing. !! For example if xESPS(i) corresponds to CH4 properties then m3i(i) must be !! the total volume of CH4 ice and gazs(i) its vapor mole fraction. !! ! Tendency of the 0th and 3rd order moment of the aerosols (spherical mode) through haze microphysics processes. REAL(kind=mm_wp), DIMENSION(:), INTENT(in) :: Hdm0as ! (m-3) REAL(kind=mm_wp), DIMENSION(:), INTENT(in) :: Hdm3as ! (m3.m-3) ! Tendency of the 0th and 3rd order moment of the aerosols (fractal mode) through haze microphysics processes. REAL(kind=mm_wp), DIMENSION(:), INTENT(in) :: Hdm0af ! (m-3) REAL(kind=mm_wp), DIMENSION(:), INTENT(in) :: Hdm3af ! (m3.m-3) ! Tendency of the 0th and 3rd order moment of the aerosols (spherical mode) through cloud microphysics processes. REAL(kind=mm_wp), DIMENSION(:), INTENT(out) :: Cdm0as ! (m-3) REAL(kind=mm_wp), DIMENSION(:), INTENT(out) :: Cdm3as ! (m3.m-3) ! Tendency of the 0th and 3rd order moment of the aerosols (fractal mode) through cloud microphysics processes. REAL(kind=mm_wp), DIMENSION(:), INTENT(out) :: Cdm0af ! (m-3) REAL(kind=mm_wp), DIMENSION(:), INTENT(out) :: Cdm3af ! (m3.m-3) ! Tendency of the 0th and 3rd order moment of the ccn distribution through cloud microphysics processes. REAL(kind=mm_wp), DIMENSION(:), INTENT(out) :: Cdm0ccn ! (m-3) REAL(kind=mm_wp), DIMENSION(:), INTENT(out) :: Cdm3ccn ! (m3.m-3) ! Tendencies of the 3rd order moments of each ice components through cloud microphysics processes. REAL(kind=mm_wp), DIMENSION(:,:), INTENT(out) :: Cdm3ice ! (m3.m-3) ! Tendencies of each condensible gaz species through cloud microphysics processes. REAL(kind=mm_wp), DIMENSION(:,:), INTENT(out) :: Cdmugas ! (mol.mol-1) ! Saturation ratio of each condensible species. REAL(kind=mm_wp), DIMENSION(:,:), INTENT(out) :: gassat ! Nucleation rate values of each condensible species (m-2.s-1). REAL(kind=mm_wp), DIMENSION(:,:),INTENT(out) :: nrate ! Growth rate values of each condensible species (m2.s-1). REAL(kind=mm_wp), DIMENSION(:,:),INTENT(out) :: grate ! Local variables: !~~~~~~~~~~~~~~~~~ INTEGER :: i,idx REAL(kind=mm_wp), DIMENSION(:,:), ALLOCATABLE :: zCdm0as,zCdm3as REAL(kind=mm_wp), DIMENSION(:,:), ALLOCATABLE :: zCdm0af,zCdm3af REAL(kind=mm_wp), DIMENSION(:,:), ALLOCATABLE :: zCdm0ccn,zCdm3ccn ! Initialization: !~~~~~~~~~~~~~~~~ ALLOCATE(zCdm0as(mm_nla,mm_nesp),zCdm3as(mm_nla,mm_nesp)) ALLOCATE(zCdm0af(mm_nla,mm_nesp),zCdm3af(mm_nla,mm_nesp)) ALLOCATE(zCdm0ccn(mm_nla,mm_nesp),zCdm3ccn(mm_nla,mm_nesp)) zCdm0as(:,:) = 0._mm_wp ; zCdm3as(:,:) = 0._mm_wp zCdm0af(:,:) = 0._mm_wp ; zCdm3af(:,:) = 0._mm_wp zCdm0ccn(:,:) = 0._mm_wp ; zCdm3ccn(:,:) = 0._mm_wp ! Calls nucleation/condensation for each species: !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ DO i = 1, mm_nesp CALL nc_esp(mm_xESPS(i),mm_gas(:,i),mm_m3ice(:,i), & Hdm0as(:),Hdm3as(:),Hdm0af(:),Hdm3af(:), & zCdm0as(:,i),zCdm3as(:,i),zCdm0af(:,i),zCdm3af(:,i), & zCdm0ccn(:,i),zCdm3ccn(:,i),Cdm3ice(:,i),Cdmugas(:,i),& gassat(:,i),nrate(:,i),grate(:,i)) ENDDO ! Retrieve the index of the maximum tendency of CCN where ice variation is not null: !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ DO i = 1, mm_nla idx = MAXLOC(zCdm0ccn(i,:),DIM=1,MASK=(Cdm3ice(i,:) /= 0._mm_wp .OR. mm_m3ice(i,:)+Cdm3ice(i,:) >= 0._mm_wp)) IF (idx == 0) THEN Cdm0as(i) = 0._mm_wp Cdm3as(i) = 0._mm_wp Cdm0af(i) = 0._mm_wp Cdm3af(i) = 0._mm_wp Cdm0ccn(i) = 0._mm_wp Cdm3ccn(i) = 0._mm_wp ELSE IF (mm_debug) THEN WRITE(*,'((a),I2.2,(a),ES10.3,(a))') "[MM_DEBUG - mm_cloud_nucond] Z(",i,") = ", mm_play(i)*1e2, & " mbar: Max aer/ccn exchange variation due to species: "//TRIM(mm_xESPS(idx)%name) ENDIF Cdm0as(i) = zCdm0as(i,idx) Cdm3as(i) = zCdm3as(i,idx) Cdm0af(i) = zCdm0af(i,idx) Cdm3af(i) = zCdm3af(i,idx) Cdm0ccn(i) = zCdm0ccn(i,idx) Cdm3ccn(i) = zCdm3ccn(i,idx) ENDIF ENDDO END SUBROUTINE mm_cloud_nucond SUBROUTINE nc_esp(xESP,Xvap,Xm3ice,Hdm0as,Hdm3as,Hdm0af,Hdm3af,& dm0as,dm3as,dm0af,dm3af,dm0ccn,dm3ccn,Xdm3ice,Xdvap,Xsat,Xnrate,Xgrate) !! Get moments tendencies through nucleation/condensation/evaporation of a given condensible species. !! The method computes the global tendencies of the aerosols, ccn and ice moments through cloud !! microphysics processes (nucleation & condensation/evaporation). !! !! @warning !! Input quantities __m0aer__,__m3aer__,__m0ccn__,__m3ccn__,__m3ice__ are assumed to be in (X.m-3), !! where X is the unit of the moment that is, a number for M0 and a volume for M3. !! __Xvap__ must be expressed in term of molar fraction. !! !! Implicit scheme: !! For nucleation we have the following equations: !! (1) dM_aer(k)/dt = - dM_ccn(k)/dt !! (2) dM_aer(k)/dt = - (4 * pi * Jhet) / rm * M_aer(k+3) !! = - (4 * pi * Jhet) / rm * alpha(k+3)/alpha(k) * rc**3 * M_aer(k) !! With: !! CST_M(k) = (4 * pi * Jhet) / rm * alpha(k+3)/alpha(k) * rc**3 * dt !! We solve: !! (3) M_aer(k)[t+dt] = (1 / (1+CST_M(k))) * M_aer(k)[t] !! Then, from (2): !! (4) M_ccn(k)[t+dt] = M_ccn(k)[t] + (CST_M(k)/(1+CST_M(k))) * M_aer(k)[t] !! ! Condensate species properties. TYPE(mm_esp), INTENT(in) :: xESP ! Gas species molar fraction on the vertical grid from __TOP__ to __GROUND__. REAL(kind=mm_wp), DIMENSION(:), INTENT(in) :: Xvap ! (mol.mol-1) ! 3rd order moment of the ice component. REAL(kind=mm_wp), DIMENSION(:), INTENT(in) :: Xm3ice ! (m3.m-3) ! Tendency of the 0th and 3rd order moment of the aerosols (spherical mode) through haze microphysics processes. REAL(kind=mm_wp), DIMENSION(:), INTENT(in) :: Hdm0as ! (m-3) REAL(kind=mm_wp), DIMENSION(:), INTENT(in) :: Hdm3as ! (m3.m-3) ! Tendency of the 0th and 3rd order moment of the aerosols (fractal mode) through haze microphysics processes. REAL(kind=mm_wp), DIMENSION(:), INTENT(in) :: Hdm0af ! (m-3) REAL(kind=mm_wp), DIMENSION(:), INTENT(in) :: Hdm3af ! (m3.m-3) ! Tendency of the 0th and 3rd order moment of the aerosols (spherical mode) through cloud microphysics processes. REAL(kind=mm_wp), DIMENSION(:), INTENT(inout) :: dm0as ! (m-3) REAL(kind=mm_wp), DIMENSION(:), INTENT(inout) :: dm3as ! (m3.m-3) ! Tendency of the 0th and 3rd order moment of the aerosols (fractal mode) through cloud microphysics processes. REAL(kind=mm_wp), DIMENSION(:), INTENT(inout) :: dm0af ! (m-3) REAL(kind=mm_wp), DIMENSION(:), INTENT(inout) :: dm3af ! (m3.m-3) ! Tendency of the 0th and 3rd order moment of the ccn through cloud microphysics processes. REAL(kind=mm_wp), DIMENSION(:), INTENT(inout) :: dm0ccn ! (m-3) REAL(kind=mm_wp), DIMENSION(:), INTENT(inout) :: dm3ccn ! (m3.m-3) ! Tendency of the 3rd order moment of the ice component. REAL(kind=mm_wp), DIMENSION(:), INTENT(inout) :: Xdm3ice ! (m3.m-3) ! Tendency of gas species. REAL(kind=mm_wp), DIMENSION(:), INTENT(inout) :: Xdvap ! (mol.mol-1) ! Saturation ratio values on the vertical grid. REAL(kind=mm_wp), DIMENSION(:), INTENT(inout) :: Xsat ! Nucleation rate values on the vertical grid for the species X. REAL(kind=mm_wp), DIMENSION(:), INTENT(inout) :: Xnrate ! (m-2.s-1) ! Growth rate values on the vertical grid for the species X. REAL(kind=mm_wp), DIMENSION(:), INTENT(inout) :: Xgrate ! (m2.s-1) ! Local variables: !~~~~~~~~~~~~~~~~~ INTEGER :: i REAL(kind=mm_wp) :: bef,aft REAL(kind=mm_wp), DIMENSION(SIZE(Xvap)) :: s_m0as,s_m3as,s_m0af,s_m3af,s_m0ccn,s_m3ccn,s_Xm3ice,s_Xvap REAL(kind=mm_wp), DIMENSION(SIZE(Xvap)) :: z_m0as,z_m3as,z_m0af,z_m3af,z_m0ccn,z_m3ccn,z_Xm3ice,z_Xvap REAL(kind=mm_wp), DIMENSION(SIZE(Xvap)) :: drad,sig,qsat,pX REAL(kind=mm_wp), DIMENSION(SIZE(Xvap)) :: Jhet_aers,cst_m0aers,cst_m3aers REAL(kind=mm_wp), DIMENSION(SIZE(Xvap)) :: Jhet_aerf,cst_m0aerf,cst_m3aerf REAL(kind=mm_wp), DIMENSION(SIZE(Xvap)) :: S_eq,g_rate REAL(kind=mm_wp), DIMENSION(SIZE(Xvap)) :: ctot,up,down,newvap REAL(kind=mm_wp), PARAMETER :: athird = 1._mm_wp / 3._mm_wp REAL(kind=mm_wp), PARAMETER :: pifac = (4._mm_wp * mm_pi) / 3._mm_wp ! Initialization: !~~~~~~~~~~~~~~~~ ! Copy input argument and convert units X.m-3 to X.kg-1 or X.mol-1 to X.kg-1. ! @warning: we must update the aerosol tracers through haze microphysics processes. ! s_XXX is the initial converted value saved: !s_m0as = mm_m0aer_s/mm_rhoair ; s_m3as = mm_m3aer_s/mm_rhoair !s_m0af = mm_m0aer_f/mm_rhoair ; s_m3af = mm_m3aer_f/mm_rhoair s_m0as = (mm_m0aer_s+Hdm0as) / mm_rhoair ; s_m3as = (mm_m3aer_s+Hdm3as) / mm_rhoair s_m0af = (mm_m0aer_f+Hdm0af) / mm_rhoair ; s_m3af = (mm_m3aer_f+Hdm3af) / mm_rhoair s_m0ccn = mm_m0ccn / mm_rhoair s_m3ccn = mm_m3ccn / mm_rhoair s_Xm3ice = Xm3ice / mm_rhoair s_Xvap = Xvap * xESP%fmol2fmas ! z_XXX is our working copy: z_m0as = s_m0as ; z_m3as = s_m3as z_m0af = s_m0af ; z_m3af = s_m3af z_m0ccn = s_m0ccn ; z_m3ccn = s_m3ccn z_Xm3ice = s_Xm3ice z_Xvap = s_Xvap ! Initialize local variables: bef = 0._mm_wp ; aft = 0._mm_wp drad(:) = 0._mm_wp sig(:) = 0._mm_wp ; qsat(:) = 0._mm_wp ; pX(:) = 0._mm_wp Jhet_aers(:) = 0._mm_wp ; cst_m0aers(:) = 0._mm_wp ; cst_m3aers(:) = 0._mm_wp Jhet_aerf(:) = 0._mm_wp ; cst_m0aerf(:) = 0._mm_wp ; cst_m3aerf(:) = 0._mm_wp S_eq(:) = 0._mm_wp ; g_rate(:) = 0._mm_wp ctot(:) = 0._mm_wp ; up(:) = 0._mm_wp ; down(:) = 0._mm_wp newvap(:) = 0._mm_wp ! Get a copy of drop radius: drad(:) = mm_drad(:) ! Surface tension (N.m-1): !~~~~~~~~~~~~~~~~~~~~~~~~~ sig = mm_sigX(mm_temp,xESP) ! Species mass mixing ratio at saturation (kg.kg-1): !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ qsat = mm_qsatX(mm_temp,mm_play,xESP) * xESP%fmol2fmas ! Partial pressure of species: !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ pX = Xvap * mm_play ! Saturation ratio: !~~~~~~~~~~~~~~~~~~ where ((z_Xvap / qsat) > 1e5) Xsat = 1e5 elsewhere Xsat = z_Xvap / qsat endwhere !~~~~~~~~~~~~~~~~~~~~~ ! GETS NUCLEATION RATE !~~~~~~~~~~~~~~~~~~~~~ ! Gets heterogeneous nucleation rate on spherical aerosols: ! Not used yet... !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ !call hetnucl_rate_aer(mm_rcs,mm_temp,xESP,pX,Xsat,Jhet_aers) ! Gets heterogeneous nucleation rate on fractal aerosols (ccn radius is the monomer): !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ call hetnucl_rate_aer((/(mm_rm, i=1,mm_nla)/),mm_temp,xESP,pX,Xsat,Jhet_aerf) ! Diagnostics Xnrate = Jhet_aers + Jhet_aerf ! /!\ IMPORTANT: ! Update CCN and aerosols moment from nucleation now ! ! Doing so should prevent a nasty bug that occurs if we want to generate clouds from scratch ! (i.e. a "dry" atmosphere without any clouds tracers already present). ! ! In such case, we do not produce ice variation on the first call of the method, ! at most only CCN are produced (i.e. dm3i == 0, dm3n != 0) ! But the rules for computing the global tendencies in mm_cloud_nucond state that the global ! variation for CCN is due to the most active species exchange. ! Spherical aerosols: ! Not used yet... !~~~~~~~~~~~~~~~~~~~~ !cst_m0aers = (4._mm_wp*mm_pi*Jhet_aers) / mm_rm * (mm_alpha_s(3._mm_wp)/mm_alpha_s(0._mm_wp)*mm_rcs**3) * mm_dt !cst_m3aers = (4._mm_wp*mm_pi*Jhet_aers) / mm_rm * (mm_alpha_s(6._mm_wp)/mm_alpha_s(3._mm_wp)*mm_rcs**3) * mm_dt !z_m0as = (1._mm_wp/(1._mm_wp+cst_m0aers)) * z_m0as !z_m3as = (1._mm_wp/(1._mm_wp+cst_m3aers)) * z_m3as !where (z_m0as <= 0._mm_wp .OR. z_m3as <= 0._mm_wp) ! z_m0as = 0._mm_wp ! z_m3as = 0._mm_wp ! z_m0ccn = z_m0ccn + s_m0as ! z_m3ccn = z_m3ccn + s_m3as !elsewhere ! z_m0ccn = z_m0ccn + cst_m0aers*z_m0as ! z_m3ccn = z_m3ccn + cst_m3aers*z_m3as !endwhere ! Fractal aerosols !~~~~~~~~~~~~~~~~~ cst_m0aerf = (4._mm_wp*mm_pi*Jhet_aerf) / mm_rm * (mm_alpha_f(3._mm_wp)/mm_alpha_f(0._mm_wp)*mm_rcf**3) * mm_dt cst_m3aerf = (4._mm_wp*mm_pi*Jhet_aerf) / mm_rm * (mm_alpha_f(6._mm_wp)/mm_alpha_f(3._mm_wp)*mm_rcf**3) * mm_dt z_m0af = (1._mm_wp/(1._mm_wp+cst_m0aerf)) * z_m0af z_m3af = (1._mm_wp/(1._mm_wp+cst_m3aerf)) * z_m3af where (z_m0af <= 0._mm_wp .OR. z_m3af <= 0._mm_wp) z_m0af = 0._mm_wp z_m3af = 0._mm_wp z_m0ccn = z_m0ccn + s_m0af z_m3ccn = z_m3ccn + s_m3af elsewhere z_m0ccn = z_m0ccn + cst_m0aerf*z_m0af z_m3ccn = z_m3ccn + cst_m3aerf*z_m3af endwhere ! Update the drop radius: !~~~~~~~~~~~~~~~~~~~~~~~~ ! Heterogeneous nucleation rate on fractal aerosols ==> we set the drop radius to the monomer radius. ! Doing so will prevent a nasty bug to occur later when ice volume is updated ! where (drad <= mm_drad_min .AND. Jhet_aerf > 0._mm_wp) drad = mm_rm endwhere !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! GETS CONDENSATION/EVAPORATION RATE !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! Equilibrium saturation near the drop: !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ S_eq = exp(min((2._mm_wp*sig*xESP%masmol) / (mm_rgas*mm_temp*xESP%rho_s*drad),30._mm_wp)) ! Gets growth rate: !~~~~~~~~~~~~~~~~~~ call growth_rate(mm_temp,mm_play,pX/Xsat,xESP,S_eq,drad,g_rate) ctot = z_Xvap + (z_Xm3ice * xESP%rho_s) up = z_Xvap + mm_dt * g_rate * 4._mm_wp * mm_pi * xESP%rho_s * drad * S_eq * z_m0ccn down = 1._mm_wp + mm_dt * g_rate * 4._mm_wp * mm_pi * xESP%rho_s * drad / qsat * z_m0ccn ! Gets new vapor X species mass mixing ratio: ! Cannot be greater than the total gas + ice and lower than nothing. !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ newvap = max(min(up/down,ctot),0._mm_wp) ! Gets "true" growth rate: !~~~~~~~~~~~~~~~~~~~~~~~~~ g_rate = g_rate * (newvap/qsat - S_eq) ! Diagnostics Xgrate = g_rate ! Update ice volume through condensation/evaporation: !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ DO i = 1, mm_nla ! Check for the specific case: no ice and sublimation IF (z_Xm3ice(i) <= 0._mm_wp .AND. g_rate(i) <= 0._mm_wp) THEN z_Xm3ice(i) = 0._mm_wp ELSE z_Xm3ice(i) = z_Xm3ice(i) + (4._mm_wp * mm_pi * drad(i) * g_rate(i) * z_m0ccn(i) * mm_dt) ! Check if there is ice left in the ccn. ! @note: only fractal aerosols for now... IF (z_Xm3ice(i) <= 0._mm_wp) THEN z_m0af(i) = z_m0af(i) + z_m0ccn(i) z_m3af(i) = z_m3af(i) + z_m3ccn(i) z_m0ccn(i) = 0._mm_wp z_m3ccn(i) = 0._mm_wp z_Xm3ice(i) = 0._mm_wp ENDIF ENDIF ENDDO ! Sanity check: !~~~~~~~~~~~~~~ IF (mm_debug) THEN DO i = 1, mm_nla bef = s_m0as(i) + s_m0af(i) + s_m0ccn(i) aft = z_m0as(i) + z_m0af(i) + z_m0ccn(i) IF (ABS(bef-aft)/bef > 1e-10_mm_wp) THEN WRITE(*,'((a),I2.2,(a),ES20.12,(a),ES20.12)') & "[MM_DEBUG - nc_esp] nc_esp("//TRIM(xESP%name)//"): M0 not conserved (z=",i,")",bef," <-> ",aft ENDIF bef = s_m3as(i) + s_m3af(i) + s_m3ccn(i) aft = z_m3as(i) + z_m3af(i) + z_m3ccn(i) IF (ABS(bef-aft)/bef > 1e-10_mm_wp) THEN WRITE(*,'((a),I2.2,(a),ES20.12,(a),ES20.12)') & "[MM_DEBUG - nc_esp] nc_esp("//TRIM(xESP%name)//"): M3 not conserved (z=",i,")",bef," <-> ",aft ENDIF ENDDO ENDIF ! Compute tendencies (in X.m-3 or X.mol-1): !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ dm0as = (z_m0as - s_m0as) * mm_rhoair ; dm3as = (z_m3as - s_m3as) * mm_rhoair dm0af = (z_m0af - s_m0af) * mm_rhoair ; dm3af = (z_m3af - s_m3af) * mm_rhoair dm0ccn = (z_m0ccn - s_m0ccn) * mm_rhoair ; dm3ccn = (z_m3ccn - s_m3ccn) * mm_rhoair Xdm3ice = (z_Xm3ice - s_Xm3ice) * mm_rhoair Xdvap = - (z_Xm3ice - s_Xm3ice) * xESP%rho_s / xESP%fmol2fmas END SUBROUTINE nc_esp SUBROUTINE hetnucl_rate_aer(rccn,temp,xESP,pvp,sat,rate) !! Get heterogeneous nucleation rate on aerosols. !! The method computes the heterogeneous nucleation rate for the given species on a aerosols of size __rccn__. !! Except __xESP__, all arguments are vectors of the same size (vertical grid). !! REAL(kind=mm_wp), INTENT(in), DIMENSION(:) :: rccn ! Radius of the cloud condensation nuclei (m). REAL(kind=mm_wp), INTENT(in), DIMENSION(:) :: temp ! Temperature (K). TYPE(mm_esp), INTENT(in) :: xESP ! Species properties. REAL(kind=mm_wp), INTENT(in), DIMENSION(:) :: pvp ! Partial vapor pressure of X species (Pa). REAL(kind=mm_wp), INTENT(in), DIMENSION(:) :: sat ! Saturation ratio of given species. REAL(kind=mm_wp), INTENT(out), DIMENSION(:) :: rate ! The nucleation rate (m-2.s-1). ! Local variables: !~~~~~~~~~~~~~~~~~ INTEGER :: i REAL(kind=mm_wp) :: S, T, r REAL(kind=mm_wp) :: sig,nX,rstar REAL(kind=mm_wp) :: x,fsh,deltaFstar REAL(kind=mm_wp) :: deltaF, gstar,zeldov ! Initialization: !~~~~~~~~~~~~~~~~ rate(:) = 0._mm_wp ! Activation condition: !~~~~~~~~~~~~~~~~~~~~~~ DO i = 1, SIZE(rccn) S = sat(i) IF (S > 1._mm_wp) THEN T = temp(i) r = rccn(i) sig = mm_sigX(T,xESP) nX = pvp(i) / (mm_kboltz*T) rstar = (2._mm_wp*sig*xESP%vol) / (mm_kboltz*T*dlog(S)) ! Curvature radius and shape factor x = r / rstar fsh = mm_fshape(xESP%mteta,x) deltaFstar = (4._mm_wp*mm_pi/3._mm_wp) * sig * (rstar**2.) * fsh deltaF = MIN(MAX((2.*xESP%fdes - xESP%fdif - deltaFstar) / (mm_kboltz*T),-100._mm_wp),100._mm_wp) IF (deltaF > -100._mm_wp) THEN gstar = ((4._mm_wp*mm_pi/3._mm_wp) * (rstar**3)) / (xESP%vol) zeldov = dsqrt(deltaFstar / (3._mm_wp*mm_pi*mm_kboltz*T*(gstar**2))) ! Heterogeneous nucleation rate rate(i) = ((zeldov*mm_kboltz*T) / (fsh*xESP%nus*xESP%mas)) * (nX*rstar)**2._mm_wp * dexp(deltaF) ENDIF ! (deltaF > -100._mm_wp) ENDIF ! (S > 1._mm_wp) ENDDO ! SIZE(rccn) RETURN END SUBROUTINE hetnucl_rate_aer SUBROUTINE growth_rate(temp,pres,pXsat,xESP,S_eq,drad,rate) !! Get growth rate through condensation/evaporation process. !! The method computes the growth rate a drop through condensation/evaporation processes: !! Except __xESP__ which is a scalar, all arguments are vectors of the same size (vertical grid). !! REAL(kind=mm_wp), INTENT(in), DIMENSION(:) :: temp ! Temperature (K). REAL(kind=mm_wp), INTENT(in), DIMENSION(:) :: pres ! Pressure level (Pa). REAL(kind=mm_wp), INTENT(in), DIMENSION(:) :: pXsat ! Saturation vapor pressure of species (Pa). TYPE(mm_esp), INTENT(in) :: xESP ! Specie properties. REAL(kind=mm_wp), INTENT(in), DIMENSION(:) :: S_eq ! Equilibrium saturation near the drop. REAL(kind=mm_wp), INTENT(in), DIMENSION(:) :: drad ! Drop radius (m). REAL(kind=mm_wp), INTENT(out), DIMENSION(:) :: rate ! Growth rate (m2.s-1). ! Local variables: !~~~~~~~~~~~~~~~~~ REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE :: k, L, knu, slf REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE :: Dv, Rc, Rd ALLOCATE(k(SIZE(temp)),L(SIZE(temp)),knu(SIZE(temp)),slf(SIZE(temp))) ALLOCATE(Dv(SIZE(temp)),Rc(SIZE(temp)),Rd(SIZE(temp))) ! Initialization: !~~~~~~~~~~~~~~~~ k(:) = 0._mm_wp ; L(:) = 0._mm_wp ; knu(:) = 0._mm_wp ; slf(:) = 0._mm_wp Dv(:) = 0._mm_wp ; Rc(:) = 0._mm_wp ; Rd(:) = 0._mm_wp ! N2 (air) Thermal conductivity: !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ k(:) = (2.857e-2_mm_wp*temp - 0.5428_mm_wp) * 4.184e-3_mm_wp ! Gas mean free path: !~~~~~~~~~~~~~~~~~~~~ L(:) = mm_lambda_air(temp,pres) ! The knudsen number of the drop: !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ knu(:) = L(:) / drad(:) ! Slip flow correction: !~~~~~~~~~~~~~~~~~~~~~~ slf(:) = knu(:) * (1.333_mm_wp + 0.71_mm_wp/knu(:)) / (1._mm_wp + 1._mm_wp/knu(:)) ! Molecular diffusivity coefficient of each species: !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Dv(:) = (1._mm_wp/3._mm_wp) * dsqrt(8._mm_wp*mm_rgas*temp(:) / (mm_pi*xESP%masmol)) * & mm_kboltz*temp(:) / (mm_pi * pres(:) * (mm_air_rad+xESP%ray)**2 * dsqrt(1._mm_wp+xESP%fmol2fmas)) ! Transitional regime: !~~~~~~~~~~~~~~~~~~~~~ Dv(:) = Dv(:) / (1._mm_wp + slf(:)) ! Latent heat resistance coefficient Rc: !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Rc(:) = (mm_LheatX(temp(:),xESP)**2 * xESP%rho_s * xESP%masmol) / (k(:) * mm_rgas * temp(:)**2) ! Molecular diffusion resistance coefficient Rd: !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Rd(:) = (mm_rgas * temp(:) * xESP%rho_s) / (Dv(:) * pXsat(:) * xESP%masmol) ! Growth rate: rdr/dt = rate * (S - S_eq): !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ rate(:) = 1._mm_wp / (S_eq(:)*Rc(:) + Rd(:)) RETURN END SUBROUTINE growth_rate !============================================================================ ! SEDIMENTATION PROCESS RELATED METHODS !============================================================================ SUBROUTINE mm_cloud_sedimentation(dm0n,dm3n,dm3i) !! Compute the tendency of clouds related moments through sedimentation process. !! !! The method computes the tendencies of moments related to cloud microphysics through !! sedimentation process. The algorithm used here differs from !! [[mm_haze_sedimentation(subroutine)]] as all moments settle with the same !! terminal velocity which is computed with the average drop radius of the size distribution. !! We simply compute an exchange matrix that stores the new positions of each cells through !! sedimentation process and then computes the matrix !! product with input moments values to get final tendencies. !! ! Tendency of the 0th order moment of the ccn distribution (m-3). REAL(kind=mm_wp), INTENT(out), DIMENSION(:) :: dm0n ! Tendency of the 3rd order moment of the ccn distribution (m3.m-3). REAL(kind=mm_wp), INTENT(out), DIMENSION(:) :: dm3n ! Tendencies of the 3rd order moment of each ice component of the cloud (m3.m-3). REAL(kind=mm_wp), INTENT(out), DIMENSION(:,:) :: dm3i ! Local variables: !~~~~~~~~~~~~~~~~~ INTEGER :: i, nm REAL(kind=mm_wp), DIMENSION(:,:), ALLOCATABLE :: moms_i, moms_f, chg_matrix ! Initialization: !~~~~~~~~~~~~~~~~ nm = 2 + mm_nesp ALLOCATE(moms_i(mm_nla,nm),moms_f(mm_nla,nm),chg_matrix(mm_nla,mm_nla)) moms_i(:,1) = mm_m0ccn * mm_dzlev moms_i(:,2) = mm_m3ccn * mm_dzlev DO i = 1, mm_nesp moms_i(:,2+i) = mm_m3ice(:,i) * mm_dzlev ENDDO ! Computes exchange matrix: !~~~~~~~~~~~~~~~~~~~~~~~~~~ CALL exchange(mm_drad,mm_drho,mm_dt,chg_matrix) ! Computes final moments values: !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ DO i = 1, nm moms_f(:,i) = MATMUL(chg_matrix,moms_i(:,i)) ENDDO ! Computes tendencies (converted in X.m-3): !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ dm0n = (moms_f(:,1) - moms_i(:,1)) / mm_dzlev dm3n = (moms_f(:,2) - moms_i(:,2)) / mm_dzlev DO i = 1, mm_nesp dm3i(:,i) = (moms_f(:,2+i) - moms_i(:,2+i)) / mm_dzlev ENDDO RETURN END SUBROUTINE mm_cloud_sedimentation SUBROUTINE exchange(rad,rhog,dt,matrix) !! Compute the exchange matrix. !! !! The subroutine computes the matrix exchange used by [[mm_cloud_sedimentation(subroutine)]] !! to compute moments tendencies through sedimentation process. Both rad and rhog must be vector with relevant !! values over the atmospheric vertical structure. !! matrix is square 2D-array with same dimension size than rad. !! ! Cloud drop radius over the atmospheric vertical structure (m). REAL(kind=mm_wp), INTENT(in), DIMENSION(:) :: rad ! Cloud drop density over the atmospheric vertical structure (kg.m-3). REAL(kind=mm_wp), INTENT(in), DIMENSION(:) :: rhog ! Timestep (s). REAL(kind=mm_wp), INTENT(in) :: dt ! The output _exchange matrix. REAL(kind=mm_wp), INTENT(out) :: matrix(:,:) ! Local variables: !~~~~~~~~~~~~~~~~~ INTEGER :: nz,i,j,jj,jinf,jsup REAL(kind=mm_wp) :: zni,znip1,xf,xft,xcnt REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE :: puit REAL(kind=mm_wp) :: cpte,cpte2 REAL(kind=mm_wp) :: zsurf INTEGER, PARAMETER:: ichx = 1 ! Initialization: !~~~~~~~~~~~~~~~~ matrix = 0._mm_wp nz = SIZE(rad) zsurf = mm_zlev(nz) ALLOCATE(puit(nz)) ! Compute exchange matrix: !~~~~~~~~~~~~~~~~~~~~~~~~~ DO i=1,nz puit(i) = 0._mm_wp xcnt = 0._mm_wp ! Computes drop move (i.e. its new positions) CALL getnzs(ichx,i,rad,rhog,dt,zni,znip1) ! Peculiar case : Ground level precipitation [znip1 < zsurf && (zni < zsurf || zni > zsurf)] ! - complete precipitation [ znip1 <= 0 && zni <= 0 ] : IF(zni <= zsurf .and. znip1 <= zsurf) THEN xft=0._mm_wp xf=1._mm_wp xcnt=xcnt+xf puit(i)=puit(i)+xf ENDIF ! - partial precipitation [ znip1 <= zsurf && zni > zsurf ] : IF (zni > zsurf .and. znip1 <= zsurf) THEN xft=(zni-zsurf)/(zni-znip1) xf=(1.-xft) xcnt=xcnt+xf puit(i)=puit(i)+xf ENDIF ! General case : no ground precipitation [ znip1 > zsurf && zni > zsurf ] IF (zni > zsurf .and. znip1 > zsurf) THEN xft = 1._mm_wp ! on a la totalite de la case xf = 0._mm_wp xcnt=xcnt+xf puit(i)=puit(i)+xf ENDIF IF (zni < znip1) THEN WRITE(*,'("[EXCHANGES] WARNING, missing this case :",2(2X,ES10.3))') zni, znip1 ENDIF ! Fix minimum level to the ground znip1 = MAX(znip1,zsurf) zni = MAX(zni,zsurf) ! Locate new "drop" position in the verical grid jsup=nz+1 jinf=nz+1 DO j=1,nz IF (zni<=mm_zlev(j).and.zni>=mm_zlev(j+1)) jsup=j IF (znip1<=mm_zlev(j).and.znip1>=mm_zlev(j+1)) jinf=j ENDDO ! Volume is out of range: (all drops have touched the ground!) ! Note: cannot happen here, it has been treated previously :) IF (jsup>=nz+1.and.jinf==jsup) THEN WRITE(*,'(a)') "[EXCHANGE] speaking: The impossible happened !" call EXIT(666) ENDIF ! Volume inside a single level IF (jsup==jinf.and.jsup<=nz) THEN xf=1._mm_wp xcnt=xcnt+xft*xf matrix(jinf,i)=matrix(jinf,i)+xft*xf ENDIF ! Volume over 2 levels IF (jinf==jsup+1) THEN xf=(zni-mm_zlev(jinf))/(zni-znip1) xcnt=xcnt+xf*xft IF(jsup<=nz) THEN matrix(jsup,i)=matrix(jsup,i)+xft*xf ENDIF xf=(mm_zlev(jinf)-znip1)/(zni-znip1) xcnt=xcnt+xf*xft IF (jinf<=nz) THEN matrix(jinf,i)=matrix(jinf,i)+xft*xf ENDIF ENDIF ! Volume over 3 or more levels IF (jinf > jsup+1) THEN ! first cell xf=(zni-mm_zlev(jsup+1))/(zni-znip1) xcnt=xcnt+xf*xft matrix(jsup,i)=matrix(jsup,i)+xft*xf ! last cell xf=(mm_zlev(jinf)-znip1)/(zni-znip1) xcnt=xcnt+xf*xft matrix(jinf,i)=matrix(jinf,i)+xft*xf ! other :) DO jj=jsup+1,jinf-1 xf=(mm_zlev(jj)-mm_zlev(jj+1))/(zni-znip1) xcnt=xcnt+xf*xft matrix(jj,i)=matrix(jj,i)+xft*xf ENDDO ENDIF ENDDO ! Checking if everything is alright if debug enabled... !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ IF (mm_debug) THEN cpte=0._mm_wp ; cpte2=0._mm_wp DO j=1,nz DO jj=1,nz cpte=cpte+matrix(jj,j) ENDDO ENDDO cpte2=cpte+sum(puit) IF (abs(cpte2-nz)>1.e-4_mm_wp) THEN WRITE(*,'("[MM_DEBUG - exchange] Tx expl (/nz):",2(2X,ES10.3))') cpte,cpte2 ENDIF ENDIF RETURN END SUBROUTINE exchange SUBROUTINE getnzs(ichx,idx,rad,rho,dt,zni,zns) !! Compute displacement of a cell under sedimentation process. !! !! The method computes the new position of a drop cell through sedimentation process. !! New positions are returned in zni and zns ouptut arguments. !! ! Velocity extrapolation control flag (0 for linear, 1 for exponential -preferred -). INTEGER, INTENT(in) :: ichx ! Initial position of the drop (subscript of vertical layers vectors). INTEGER, INTENT(in) :: idx ! Cloud drop radius over the atmospheric vertical structure (m). REAL(kind=mm_wp), INTENT(in), DIMENSION(:) :: rad ! Cloud drop density over the atmospheric vertical structure (kg.m-3). REAL(kind=mm_wp), INTENT(in), DIMENSION(:) :: rho ! Timestep (s). REAL(kind=mm_wp), INTENT(in) :: dt ! Final layer upper boundary position (m). REAL(kind=mm_wp), INTENT(out) :: zni ! Final layer lower boundary position (m). REAL(kind=mm_wp), INTENT(out) :: zns ! Local variables: !~~~~~~~~~~~~~~~~~ INTEGER :: i,nz REAL(kind=mm_wp) :: ws,wi,w,zi,zs,rhoair_mid REAL(kind=mm_wp) :: alpha,argexp,v0,arg1,arg2 REAL(kind=mm_wp), PARAMETER :: es = 30._mm_wp nz = SIZE(rad) ! Linear extrapolation of velocity: !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ IF (ichx==0) THEN ! Velocity upper interface ws = wsettle(mm_plev(idx),mm_btemp(idx),mm_zlev(idx),mm_rhoair(idx),rho(idx),rad(idx)) IF (idx==nz) THEN ! Veloctity center layer rhoair_mid = (mm_rhoair(idx-1) + mm_rhoair(idx)) / 2._mm_wp wi = wsettle(mm_play(idx),mm_temp(idx),mm_zlay(idx),rhoair_mid,rho(idx),rad(idx)) ELSEIF(idx0._mm_wp.AND.(zs-zi)/=0._mm_wp) alpha = dlog(ws/wi)/(zs-zi) ! Alpha < 0 if wi > ws ! -es < argexp < es argexp=MAX(MIN(alpha*zs,es),-es) v0 = ws/dexp(argexp) arg1=1._mm_wp+v0*alpha*dexp(argexp)*dt argexp=MAX(MIN(alpha*(mm_zlev(idx)-mm_dzlev(idx)),es),-es) arg2=1._mm_wp+v0*alpha*dexp(argexp)*dt IF (arg1<=0._mm_wp.OR.arg2<=0._mm_wp) THEN ! correct velocity ! divides the velocity argument in arg1 and arg2 : ! argX=1+alpha*v0*exp(alpha*z)*dt <==> argX-1=alpha*v0*exp(alpha*z)*dt ! ==> argX' = 1 + (argX-1)/2 <==> argX' = (1+argX)/2. DO i=1,25 IF (arg1<=0._mm_wp.OR.arg2<=0._mm_wp) THEN IF (mm_debug) & WRITE(*,'((a),I2.2,(a))') "[MM_DEBUG - getnzs] Must adjust velocity (iter:",i,"/25)" arg1=(arg1+1._mm_wp)/2._mm_wp ; arg2=(arg2+1._mm_wp)/2._mm_wp ELSE EXIT ENDIF ENDDO ! We have to stop IF (i>25) THEN WRITE(*,'(a)')"[ERROR - getnzs] Cannot adjust velocities !" call EXIT(111) ENDIF ENDIF zni = mm_zlev(idx)-dlog(arg1)/alpha zns = mm_zlev(idx)-mm_dzlev(idx)-dlog(arg2)/alpha RETURN ENDIF ! end of ichx END SUBROUTINE getnzs ELEMENTAL FUNCTION wsettle(p,t,z,rhoair,rho,rad) RESULT(w) !! Compute the settling velocity of a spherical particle. !! !! The method computes the effective settling velocity of spherical particle of radius rad. !! REAL(kind=mm_wp), INTENT(in) :: p ! The pressure level (Pa). REAL(kind=mm_wp), INTENT(in) :: t ! The temperature (K). REAL(kind=mm_wp), INTENT(in) :: z ! The altitude level (m). REAL(kind=mm_wp), INTENT(in) :: rhoair ! Density of air (kg.m-3). REAL(kind=mm_wp), INTENT(in) :: rho ! Density of the particle (kg.m-3). REAL(kind=mm_wp), INTENT(in) :: rad ! Radius of the particle (m). REAL(kind=mm_wp) :: w ! Settling velocity (m.s-1). ! Local variables: REAL(kind=mm_wp), PARAMETER :: wmax = 30.0_mm_wp ! Maximal settling velocity (m.s-1) REAL(kind=mm_wp), PARAMETER :: mrcoef = 0.74_mm_wp ! Molecular reflexion coefficient REAL(kind=mm_wp) :: thermal_w REAL(kind=mm_wp) :: kn, Fc, Us ! Molecular's case !~~~~~~~~~~~~~~~~~ ! Thermal velocity thermal_w = sqrt((8._mm_wp * mm_kboltz * t) / (mm_pi * mm_air_mmas)) ! Computes settling velocity w = mrcoef * mm_effg(z) * (mm_rhoaer/rhoair) * rad / thermal_w ! Hydrodynamical's case ! In fact: transitory case which is the general case !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! Knudsen number !kn = (mm_kboltz * t) / (p * 4._mm_wp * sqrt(2._mm_wp) * mm_pi * mm_air_rad**2) / rad ! Computes slip-flow correction !Fc = (1._mm_wp + 1.257_mm_wp*kn + 0.4_mm_wp*kn*dexp(-1.1_mm_wp/kn)) ! Computes Stokes settling velocity !Us = (2._mm_wp * rad**2 * rho * mm_effg(z)) / (9._mm_wp * mm_eta_air(t)) ! Cunningham-Millikan correction !w = Us * Fc ! Velocity limit: drop deformation (Lorenz 1993) !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ w = 1._mm_wp / ((1._mm_wp / w) + (1._mm_wp / wmax)) RETURN END FUNCTION wsettle FUNCTION get_mass_flux(rho,M3) RESULT(flx) !! Get the mass flux of clouds related moment through sedimention. !! !! @warning !! The method is only valid for cloud moments (i.e. ice or ccn). !! It calls [[wsettle(function)]] that compute the mean settling velocity of a cloud drop. !! !! @note !! The computed flux is always positive. !! ! Tracer density (kg.m-3). REAL(kind=mm_wp), INTENT(in) :: rho ! Vertical profile of the total volume of tracer from __TOP__ to __GROUND__ (m3.m-3). REAL(kind=mm_wp), INTENT(in), DIMENSION(:) :: M3 ! Mass sedimentation fluxes at each layer from __TOP__ to __GROUND__ (kg.m-2.s-1). REAL(kind=mm_wp), DIMENSION(SIZE(M3)) :: flx ! Local variables INTEGER :: i REAL(kind=mm_wp), DIMENSION(SIZE(M3)) :: w_cld REAL(kind=mm_wp), SAVE :: pifac = (4._mm_wp * mm_pi) / 3._mm_wp ! Cloud drop settling velocity !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ w_cld = wsettle(mm_play,mm_temp,mm_zlay,mm_rhoair,mm_drho,mm_drad) ! Computes flux through sedimention !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ flx = rho * pifac * M3 * w_cld RETURN END FUNCTION get_mass_flux END MODULE MP2M_CLOUDS