! Copyright (2013-2015,2017,2022-2023) Université de Reims Champagne-Ardenne ! Contributors : J. Burgalat (GSMA, URCA), B. de Batz de Trenquelléon (GSMA, URCA) ! email of the author : jeremie.burgalat@univ-reims.fr ! ! This software is a computer program whose purpose is to compute ! microphysics processes using a two-moments scheme. ! ! This library is governed by the CeCILL-B license under French law and ! abiding by the rules of distribution of free software. You can use, ! modify and/ or redistribute the software under the terms of the CeCILL-B ! license as circulated by CEA, CNRS and INRIA at the following URL ! "http://www.cecill.info". ! ! As a counterpart to the access to the source code and rights to copy, ! modify and redistribute granted by the license, users are provided only ! with a limited warranty and the software's author, the holder of the ! economic rights, and the successive licensors have only limited ! liability. ! ! In this respect, the user's attention is drawn to the risks associated ! with loading, using, modifying and/or developing or reproducing the ! software by the user in light of its specific status of free software, ! that may mean that it is complicated to manipulate, and that also ! therefore means that it is reserved for developers and experienced ! professionals having in-depth computer knowledge. Users are therefore ! encouraged to load and test the software's suitability as regards their ! requirements in conditions enabling the security of their systems and/or ! data to be ensured and, more generally, to use and operate it in the ! same conditions as regards security. ! ! The fact that you are presently reading this means that you have had ! knowledge of the CeCILL-B license and that you accept its terms. !! file: mm_clouds.f90 !! summary: Clouds microphysics module !! author: J. Burgalat !! date: 2013-2015,2017,2022-2023 !! modifications: B. de Batz de Trenquelléon MODULE MM_CLOUDS !! Clouds microphysics module. !! !! The module contains all definitions of the microphysics processes related to clouds: !! !! - [nucleation](page/clouds.html#nucleation) !! - [condensation](page/clouds.html#condensation) !! - [sedimentation](page/clouds.html#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. !! !! @note !! Tendencies returned by public methods are always defined from __TOP__ of the atmosphere to the !! __GROUND__. USE MM_MPREC USE MM_GLOBALS USE MM_METHODS IMPLICIT NONE PRIVATE PUBLIC :: mm_cloud_microphysics, mm_cloud_sedimentation, mm_cloud_nucond CONTAINS !============================================================================ ! CLOUDS MICROPHYSICS INTERFACE SUBROUTINE !============================================================================ SUBROUTINE mm_cloud_microphysics(dm0a,dm3a,dm0n,dm3n,dm3i,dgazs) !! Get the evolution of moments tracers through clouds microphysics processes. !! !! The subroutine is a wrapper to the clouds microphysics methods. It computes the tendencies of moments !! tracers for nucleation, condensation and sedimentation processes for the atmospheric column. !! !! @note !! Both __dm3i__ and __dgazs__ are 2D-array with the vertical layers in first dimension and the number !! of ice components in the second. REAL(kind=mm_wp), DIMENSION(:), INTENT(out) :: dm0a !! Tendency of the 0th order moment of the aerosols (fractal mode) (\(m^{-3}\)). REAL(kind=mm_wp), DIMENSION(:), INTENT(out) :: dm3a !! Tendency of the 3rd order moment of the aerosols distribution (fractal mode) (\(m^{3}.m^{-3}\)) . REAL(kind=mm_wp), DIMENSION(:), INTENT(out) :: dm0n !! Tendency of the 0th order moment of the ccn distribution (\(m^{-3}\)). REAL(kind=mm_wp), DIMENSION(:), INTENT(out) :: dm3n !! Tendency of the 3rd order moment of the ccn distribution (\(m^{3}.m^{-3}\)). REAL(kind=mm_wp), DIMENSION(:,:), INTENT(out) :: dm3i !! Tendencies of the 3rd order moments of each ice components (\(m^{3}.m^{-3}\)). REAL(kind=mm_wp), DIMENSION(:,:), INTENT(out) :: dgazs !! Tendencies of each condensible gaz species (\(mol.mol^{-1}\)). REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE :: zdm0n,zdm3n REAL(kind=mm_wp), DIMENSION(:,:), ALLOCATABLE :: zdm3i INTEGER :: i dm0a = 0._mm_wp ; dm3a = 0._mm_wp dm0n = 0._mm_wp ; dm3n = 0._mm_wp dm3i = 0._mm_wp ; dgazs = 0._mm_wp IF (mm_w_cloud_nucond) THEN ! Calls condensation/nucleation (and update saturation ratio diagnostic) call mm_cloud_nucond(dm0a,dm3a,dm0n,dm3n,dm3i,dgazs,mm_gazs_sat) ENDIF IF (mm_w_cloud_sed) THEN ! Calls sedimentation ALLOCATE(zdm0n(mm_nla),zdm3n(mm_nla),zdm3i(mm_nla,mm_nesp)) call mm_cloud_sedimentation(zdm0n,zdm3n,zdm3i) ! Computes settling velocity [m.s-1] of clouds (ccn and ices) mm_ccn_vsed(:) = wsettle(mm_play,mm_temp,mm_zlay,mm_drho,mm_drad) ! Computes flux [kg.m-2.s-1] and precipitation [m.iphysiq] of ccn mm_ccn_flux(:) = get_mass_flux(mm_rhoaer,mm_m3ccn(:)) mm_ccn_prec = SUM(zdm3n*mm_dzlev) ! Computes flux [kg.m-2.s-1] and precipitation [m.iphysiq] of ices DO i = 1, mm_nesp mm_ice_fluxes(:,i) = get_mass_flux(mm_xESPS(i)%rho,(3._mm_wp*mm_m3ice(:,i))/(4._mm_wp*mm_pi)) mm_ice_prec(i) = SUM(zdm3i(:,i)*mm_dzlev) ENDDO ! updates tendencies dm0n = dm0n + zdm0n dm3n = dm3n + zdm3n dm3i = dm3i + zdm3i ENDIF END SUBROUTINE mm_cloud_microphysics !----------------------------------------------------------------------------- ! NUCLEATION/CONDENSATION PROCESS RELATED METHODS !----------------------------------------------------------------------------- SUBROUTINE mm_cloud_nucond(dm0a,dm3a,dm0n,dm3n,dm3i,dgazs,gazsat) !! Get moments tendencies through nucleation/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 only loss of CCN can occur. !! !! 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(IDX)__ !! corresponds to \(CH_{4}\) properties then __m3i(IDX)__ must be the total volume of solid !! \(CH_{4}\) (ice) and __gazs(IDX)__ its vapor mole fraction. REAL(kind=mm_wp), DIMENSION(:), INTENT(out) :: dm0a !! Tendency of the 0th order moment of the aerosols (fractal mode) (\(m^{-3}\)). REAL(kind=mm_wp), DIMENSION(:), INTENT(out) :: dm3a !! Tendency of the 3rd order moment of the aerosols distribution (fractal mode) (\(m^{3}.m^{-3}\)). REAL(kind=mm_wp), DIMENSION(:), INTENT(out) :: dm0n !! Tendency of the 0th order moment of the aerosols distribution (fractal mode) (\(m^{-3}\)). REAL(kind=mm_wp), DIMENSION(:), INTENT(out) :: dm3n !! Tendency of the 3rd order moment of the ccn distribution (\(m^{3}.m^{-3}\)). REAL(kind=mm_wp), DIMENSION(:,:), INTENT(out) :: dm3i !! Tendencies of the 3rd order moments of each ice components (\(m^{3}.m^{-3}\)). REAL(kind=mm_wp), DIMENSION(:,:), INTENT(out) :: dgazs !! Tendencies of each condensible gaz species (\(mol.mol^{-1}\)) . REAL(kind=mm_wp), DIMENSION(:,:), INTENT(out) :: gazsat !! Saturation ratio of each condensible specie. INTEGER :: i,idx REAL(kind=mm_wp), DIMENSION(:,:), ALLOCATABLE :: zdm0a,zdm3a,zdm0n,zdm3n ALLOCATE(zdm0a(mm_nla,mm_nesp),zdm3a(mm_nla,mm_nesp), & zdm0n(mm_nla,mm_nesp),zdm3n(mm_nla,mm_nesp)) zdm0a(:,:) = 0._mm_wp ; zdm3a(:,:) = 0._mm_wp zdm0n(:,:) = 0._mm_wp ; zdm3n(:,:) = 0._mm_wp DO i = 1, mm_nesp call nc_esp(mm_xESPS(i),mm_gazs(:,i),mm_m3ice(:,i),dgazs(:,i),dm3i(:,i), & zdm0a(:,i),zdm3a(:,i),zdm0n(:,i),zdm3n(:,i),gazsat(:,i)) ENDDO DO i=1, mm_nla ! retrieve the index of the maximum tendency of CCN where ice variation is not null. idx = MAXLOC(zdm0n(i,:),DIM=1,MASK=(dm3i(i,:) /= 0._mm_wp .OR. mm_m3ice(i,:)+dm3i(i,:) >= 0._mm_wp)) IF (idx == 0) THEN dm0n(i) = 0._mm_wp dm3n(i) = 0._mm_wp dm0a(i) = 0._mm_wp dm3a(i) = 0._mm_wp ELSE IF (mm_debug .AND. ABS(zdm0n(i,idx)) > 1e3) THEN WRITE(*,'((a),I2.2,(a),ES10.3,(a))') "Z(",i,") = ",mm_play(i)*1e2, & " mbar: Max aer/ccn exchange variation due to specie: "//TRIM(mm_xESPS(idx)%name) ENDIF dm0n(i) = zdm0n(i,idx) dm3n(i) = zdm3n(i,idx) dm0a(i) = zdm0a(i,idx) dm3a(i) = zdm3a(i,idx) ENDIF ENDDO END SUBROUTINE mm_cloud_nucond SUBROUTINE nc_esp(xESP,vapX,m3iX,dvapX,dm3iX,dm0aer,dm3aer,dm0ccn,dm3ccn,Xsat) !! Get moments tendencies through nucleation/condensation/evaporation of a given condensible specie. !! !! The method computes the global tendencies of the aerosols, ccn and "ice" moments through cloud !! microphysics processes (nucleation & condensation). !! !! @warning !! Input quantities __m3iX__,__m3iO__, __m0aer__,__m3aer__, __m0ccn__,__m3ccn__ 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 - \(m^3\) !! for M3) ; __vapX__ must be expressed in term of molar fraction. TYPE(mm_esp), INTENT(in) :: xESP !! Condensate specie properties. REAL(kind=mm_wp),INTENT(in), DIMENSION(:) :: vapX !! Gas specie molar fraction on the vertical grid from __TOP__ to __GROUND__ (\(mol.mol^{-1}\)). REAL(kind=mm_wp),INTENT(in), DIMENSION(:) :: m3iX !! 3rd order moment of the ice component (\(m^{3}.m^{-3}\)). REAL(kind=mm_wp),INTENT(out), DIMENSION(:) :: dvapX !! Tendency of gas specie (\(mol.mol^{-1}\)). REAL(kind=mm_wp),INTENT(out), DIMENSION(:) :: dm3iX !! Tendency of the 3rd order moment of the ice component (\(m^{3}.m^{-3}\)). REAL(kind=mm_wp),INTENT(out), DIMENSION(:) :: dm0aer !! Tendency of the 0th order moment of the fractal mode distribution (\(m^{-3}\)). REAL(kind=mm_wp),INTENT(out), DIMENSION(:) :: dm3aer !! Tendency of the 3rd order moment of the fractal mode distribution (\(m^{3}.m^{-3}\)). REAL(kind=mm_wp),INTENT(out), DIMENSION(:) :: dm0ccn !! Tendency of the 0th order moment of the ccn distribution (\(m^{-3}\)). REAL(kind=mm_wp),INTENT(out), DIMENSION(:) :: dm3ccn !! Tendency of the 3rd order moment of the ccn distribution (\(m^{3}.m^{-3}\)). REAL(kind=mm_wp),INTENT(out), DIMENSION(:) :: Xsat !! Saturation ratio values on the vertical grid (--). INTEGER :: i REAL(kind=mm_wp) :: bef,aft REAL(kind=mm_wp), DIMENSION(SIZE(vapX)) :: sm0a,sm3a,sm0n,sm3n,sm3iX REAL(kind=mm_wp), DIMENSION(SIZE(vapX)) :: zm0a,zm3a,zm0n,zm3n,zm3iX,zvapX REAL(kind=mm_wp), DIMENSION(SIZE(vapX)) :: pX,sig,qsat,seq,up,down,ctot,newvap,nucr,grate,cm0,cm3,drad ! get a copy of drop radius. drad(:) = mm_drad(:) ! Initialization : ! Copy input argument and convert units X.m-3 -> X.kg-1 ! sXXX is the initial converted value saved sm3iX = m3iX/mm_rhoair sm0a = mm_m0aer_f/mm_rhoair ; sm3a = mm_m3aer_f/mm_rhoair sm0n = mm_m0ccn/mm_rhoair ; sm3n = mm_m3ccn/mm_rhoair ! zXXX is our working copy zm3iX = sm3iX ; zm0a = sm0a ; zm3a = sm3a ; zm0n = sm0n ; zm3n = sm3n ! Molar fraction of X specie is set in mass mixing ratio [kg.kg-1] zvapX = vapX * xESP%fmol2fmas ! Surface tension [N.m-1] sig = mm_sigX(mm_temp,xESP) ! X specie mass mixing ratio at saturation [kg.kg-1] qsat = mm_ysatX(mm_temp,mm_play,xESP) * xESP%fmol2fmas ! partial pressure of X specie pX = vapX * mm_play ! Saturation ratio Xsat = zvapX / qsat ! Gets nucleation rate (ccn radius is the monomer !) call nuc_rate((/(mm_rm, i=1,mm_nla)/),mm_temp,xESP,pX,Xsat,nucr) ! 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 specie exchange. ! ! for nucleation we have the following equations: ! dMa(k)/dt = - dMn(k)/dt (conservation of aerosols+ccn) (1) ! dMa(k)/dt = - 4*PI*nucr/rm * Ma(k+3) (2) ! = - 4*PI*nucr/rm * alpha(k+3)/alpha(k) * rc**3 * Ma(k) ! With: ! - Ma(k): k-th order moment of aerosol ! - Mn(k): k-th order moment of ccn ! - nucr : the nucleation rate. ! We solve (implicit scheme) : ! CST_M(k) = 4*PI*nucr/rm * alpha(k+3)/alpha(k)*rc**3 * dt ! Ma(k)[t+dt] = 1/(1+CST_M(k)) * Ma(k)[t] (3) ! Then, from eq. 2: ! Mn(k)[t+dt] = Mn(k)[t] + CST_M(k)/(1+CST_M(k))*Ma(k)[t] (4) cm0 = 4._mm_wp*mm_pi*nucr/mm_rm*mm_alpha_f(3._mm_wp)*mm_rcf**3*mm_dt cm3 = 4._mm_wp*mm_pi*nucr/mm_rm*mm_alpha_f(6._mm_wp)/mm_alpha_f(3._mm_wp)*mm_rcf**3*mm_dt zm0a = 1._mm_wp/(1._mm_wp+cm0) * zm0a zm3a = 1._mm_wp/(1._mm_wp+cm3) * zm3a WHERE (zm0a <= 0._mm_wp .OR. zm3a <= 0._mm_wp) zm0a = 0._mm_wp zm3a = 0._mm_wp zm0n = zm0n + sm0a zm3n = zm3n + sm3a ELSEWHERE zm0n = zm0n + cm0*zm0a zm3n = zm3n + cm3*zm3a ENDWHERE ! update the drop radius (we probably should recompute totally the radius to be in better agreement with the other moments) ! We must manage the case where there is no ices and no ccn ==> drop radius is ZERO, ! but conditions are met to spawn nucleation process: creation of ccn. ! Then 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 (nucr > 0._mm_wp .AND. drad <= mm_drad_min) drad = mm_rm ENDWHERE ! Equilibrium saturation near the drop seq = exp(min(30._mm_wp,2._mm_wp*sig*xESP%masmol/(xESP%rho*mm_rgas*mm_temp*drad))) ! Gets growth rate call growth_rate(mm_temp,mm_play,pX/Xsat,xESP,seq,drad,grate) ctot = zvapx + xESP%rho * zm3iX up = zvapx + mm_dt * grate * 4._mm_wp * mm_pi * xESP%rho * drad * seq * zm0n down = 1._mm_wp + mm_dt * grate * 4._mm_wp * mm_pi * xESP%rho * drad / qsat * zm0n ! gets new vapor X specie 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 grate = grate * (newvap/qsat - seq) ! computes tendencies through condensation DO i=1,mm_nla ! check for the specific case : NO ICE and SUBLIMATION IF (zm3iX(i) <= 0._mm_wp .AND. grate(i) <= 0._mm_wp) THEN zm3iX(i) = 0._mm_wp ELSE ! update ice volume ... zm3iX(i) = zm3iX(i) + mm_dt*grate(i)*4._mm_wp*mm_pi*drad(i)*zm0n(i) ! ... and check if there is ice left in the ccn IF (zm3iX(i) <= 0._mm_wp) THEN zm3iX(i) = 0._mm_wp zm0a(i) = zm0a(i) + zm0n(i) ; zm0n(i) = 0._mm_wp zm3a(i) = zm3a(i) + zm3n(i) ; zm3n(i) = 0._mm_wp ENDIF ENDIF ENDDO ! Computes balance IF (mm_debug) THEN DO i=1,mm_nla bef = sm0a(i) + sm0n(i) aft = zm0a(i) + zm0n(i) IF (ABS(bef-aft)/bef > 1e-10_mm_wp) THEN WRITE(*,'((a),I2.2,(a),ES20.12,(a),ES20.12)') & "[DEBUG] nc_esp("//TRIM(xESP%name)//"): M0 not conserved (z=",i,")",bef," <-> ",aft ENDIF bef = sm3a(i) + sm3n(i) aft = zm3a(i) + zm3n(i) IF (ABS(bef-aft)/bef > 1e-10_mm_wp) THEN WRITE(*,'((a),I2.2,(a),ES20.12,(a),ES20.12)') & "[DEBUG] nc_esp("//TRIM(xESP%name)//"): M3 not conserved (z=",i,")",bef," <-> ",aft ENDIF ENDDO ENDIF ! compute tendencies: ! all of these tendencies are in X.m-3 ! dm0aer = (zm0a - sm0a)*mm_rhoair dm3aer = (zm3a - sm3a)*mm_rhoair dm0ccn = (zm0n - sm0n)*mm_rhoair dm3ccn = (zm3n - sm3n)*mm_rhoair dm3iX = (zm3iX - sm3iX) ! this one in X.kg-1 (temporary) ! dvapX = -xESP%rho * dm3iX / xESP%fmol2fmas ! in order to compute this one in mol.mol-1 dm3iX = dm3iX*mm_rhoair ! update ice tendencies in X.m-3 ! END SUBROUTINE nc_esp SUBROUTINE nuc_rate(rccn,temp,xESP,pvp,sat,rate) !! Get nucleation rate. !! !! The method computes the heterogeneous nucleation rate for the given specie on a fractal particle !! 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). REAL(kind=mm_wp), INTENT(in), DIMENSION(:) :: pvp !! Partial vapor pressure of X specie (Pa). REAL(kind=mm_wp), INTENT(in), DIMENSION(:) :: sat !! Saturation ratio of given specie (--). TYPE(mm_esp), INTENT(in) :: xESP !! X specie properties (--). REAL(kind=mm_wp), INTENT(out), DIMENSION(:) :: rate !! The nucleation rate (\(m^-{2}.s^{-1}\)). INTEGER :: i REAL(kind=mm_wp) :: r,t,s,sig,nX,rstar,gstar,x,zeldov,deltaf,fsh,fstar 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 x = r/rstar fsh = mm_fshape(xESP%mteta,x) fstar = (4._mm_wp*mm_pi/3._mm_wp)*sig*(rstar**2.)*fsh deltaf=MIN(MAX((2.*mm_fdes-mm_fdif-fstar)/(mm_kboltz*t),-100._mm_wp),100._mm_wp) IF (deltaf > -100._mm_wp) THEN gstar = 4._mm_wp*mm_pi*(rstar**3)/(3._mm_wp*xESP%vol) zeldov = dsqrt(fstar/(3._mm_wp*mm_pi*mm_kboltz*t*(gstar**2))) rate(i)= zeldov*mm_kboltz*t*(nX*rstar)**2._mm_wp*dexp(deltaf)/(fsh*mm_nus*xESP%mas) ENDIF ENDIF ENDDO RETURN END SUBROUTINE nuc_rate SUBROUTINE growth_rate(temp,pres,pXsat,xESP,seq,drad,rate) !! Get growth rate through condensation/evaporation process. !! !! The method computes the growth rate a drop through condensation/evaporation processes: !! !! $$ r \times \frac{dr}{dt} = g_{rate} \times (S - S_{eq}) $$ !! !! 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 specie (Pa). TYPE(mm_esp), INTENT(in) :: xESP !! Specie properties. REAL(kind=mm_wp), INTENT(in), DIMENSION(:) :: seq !! 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 (\(m^{2}.s^{-1}\)). REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE :: k,knu,slf,rkc,rdc,l,dv INTEGER :: n n = SIZE(temp) ALLOCATE(k(n),knu(n),slf(n),rkc(n),rdc(n),l(n),dv(n)) ! N2 (air) Thermal conductivity (where does it come from ?) k(:) = ( 2.857e-2_mm_wp*temp-0.5428_mm_wp)*4.184e-3_mm_wp ! Gas mean free path l(:) = mm_lambda_g(temp,pres) ! Diffusion coefficient of X gas 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)) knu(:) = l(:)/drad(:) ! The knudsen number of the drop slf(:) = (1.333_mm_wp+0.71_mm_wp/knu(:))/(1._mm_wp+1._mm_wp/knu(:)) ! Slip flow correction Dv(:) = Dv(:)/(1._mm_wp+slf(:)*knu(:)) ! latent heat resistance coefficient rkc(:) = mm_lheatX(temp(:),xESP)**2 * xESP%rho * xESP%masmol / (k(:)*mm_rgas*temp(:)**2) ! Diffusion resistance coefficient rdc(:) = mm_rgas * temp(:) * xESP%rho / (Dv(:)*pXsat(:)*xESP%masmol) ! Growth rate: rdr/dt = rate * (S-Seq) ; rate is returned rate(:) = 1._mm_wp / (seq(:)*rkc(:)+rdc(:)) 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(module):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. REAL(kind=mm_wp), INTENT(out), DIMENSION(:) :: dm0n !! Tendency of the 0th order moment of the ccn distribution (\(m^{-3}\)). REAL(kind=mm_wp), INTENT(out), DIMENSION(:) :: dm3n !! Tendency of the 3rd order moment of the ccn distribution (\(m^{3}.m^{-3}\)). REAL(kind=mm_wp), INTENT(out), DIMENSION(:,:) :: dm3i !! Tendencies of the 3rd order moment of each ice component of the cloud (\(m^{3}m^{-3}\)). INTEGER :: im,nm REAL(kind=mm_wp), DIMENSION(:,:), ALLOCATABLE :: moms, momsf,chg_matrix nm = 2 + mm_nesp ALLOCATE(moms(mm_nla,nm),momsf(mm_nla,nm),chg_matrix(mm_nla,mm_nla)) ! Initializes moms moms(:,1) = mm_m0ccn * mm_dzlev moms(:,2) = mm_m3ccn * mm_dzlev DO im=1,mm_nesp ; moms(:,2+im) = mm_m3ice(:,im) * mm_dzlev ; ENDDO ! Computes exchange matrix CALL exchange(mm_drad,mm_drho,mm_dt,chg_matrix) ! Computes final moments values DO im=1,nm ; momsf(:,im) = MATMUL(chg_matrix,moms(:,im)) ; ENDDO ! Computes tendencies (converted in X.m-3) dm0n = (momsf(:,1)-moms(:,1))/mm_dzlev dm3n = (momsf(:,2)-moms(:,2))/mm_dzlev DO im=1,mm_nesp ; dm3i(:,im) = (momsf(:,2+im)-moms(:,2+im))/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_clouds(module):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__. REAL(kind=mm_wp), INTENT(in), DIMENSION(:) :: rad !! Cloud drop radius over the atmospheric vertical structure (m). REAL(kind=mm_wp), INTENT(in), DIMENSION(:) :: rhog !! Cloud drop density over the atmospheric vertical structure (\(kg.m^{-3}\)). REAL(kind=mm_wp), INTENT(in) :: dt !! Timestep (s). REAL(kind=mm_wp), INTENT(out) :: matrix(:,:) !! The output _exchange matrix_. 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 matrix = 0._mm_wp ; nz = SIZE(rad) ; ALLOCATE(puit(nz)) zsurf = mm_zlev(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(*,'("[EXCHANGE] speaking: 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 as !! descibed in the following scheme: !! !! ![Cloud sedimentation scheme](|media|/cloud_sed_scheme.svg) !! !! New positions are returned in __zni__ and __zns__ ouptut arguments. !! !! @note !! The method uses directly [[mm_globals(module):mm_play(variable)]], [[mm_globals(module):mm_plev(variable)]], !! [[mm_globals(module):mm_temp(variable)]],[[mm_globals(module):mm_btemp(variable)]], !! [[mm_globals(module):mm_zlay(variable)]] and [[mm_globals(module):mm_zlev(variable)]] and uses __idx__ to !! get the relevant value to use on the vertical grid. INTEGER, INTENT(in) :: ichx !! Velocity extrapolation control flag (0 for linear, 1 for exponential -preferred -). INTEGER, INTENT(in) :: idx !! Initial position of the drop (subscript of vertical layers vectors). REAL(kind=mm_wp), INTENT(in), DIMENSION(:) :: rad !! Cloud drop radius over the atmospheric vertical structure (m). REAL(kind=mm_wp), INTENT(in), DIMENSION(:) :: rho !! Cloud drop density over the atmospheric vertical structure (\(kg.m^{-3}\)). REAL(kind=mm_wp), INTENT(in) :: dt !! Timestep (s). REAL(kind=mm_wp), INTENT(out) :: zni !! Final layer upper boundary position (m). REAL(kind=mm_wp), INTENT(out) :: zns !! Final layer lower boundary position (m). REAL(kind=mm_wp) :: ws,wi,w,zi,zs REAL(kind=mm_wp) :: alpha,argexp,v0,arg1,arg2 INTEGER :: i,nz 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),rho(idx),rad(idx)) IF (idx==nz) THEN ! veloctity center layer wi = wsettle(mm_play(idx),mm_temp(idx),mm_zlay(idx),rho(idx),rad(idx)) ELSEIF(idx 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))') "[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 ! sh** we have to stop IF (i>25) THEN WRITE(*,'(a)')"[getnzs] speaking:" WRITE(*,'(a)') "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 SUBROUTINE getnzs ELEMENTAL FUNCTION wsettle(p,t,z,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__. It accounts for the slip-flow transition (no approximation). 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) :: 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}\)). REAL(kind=mm_wp) :: Us, Fc, kn REAL(kind=mm_wp), PARAMETER :: ra = 1.75e-10_mm_wp ! Knudsen number kn = (mm_kboltz * t) / (p * 4._mm_wp * sqrt(2._mm_wp) * mm_pi * ra**2) / rad ! Computes slip-flow correction : Fc = 1 + (mm_akn * mm_lambda_g(t,p) / rad) Fc = (1._mm_wp + 1.2517_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_g(t)) ! Computes settling velocity (correction factor : x2.0) w = Us * Fc * 2._mm_wp 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 !! [[mm_clouds(module):wsettle(function)]] that compute the _mean_ settling velocity of a !! cloud drop. !! !! @note !! The computed flux is always positive. REAL(kind=mm_wp), INTENT(in) :: rho !! Tracer density (\(kg.m^{-3}\)). REAL(kind=mm_wp), INTENT(in), DIMENSION(:) :: m3 !! Vertical profile of the total volume of tracer (i.e. M3) from __TOP__ to __GROUND__ (\(m^{3}.m^{-3}\)). REAL(kind=mm_wp), DIMENSION(SIZE(m3)) :: flx !! Mass sedimentation fluxes at each layer from __TOP__ to __GROUND__ (\(kg.m^{-2}.s^{-1}\)). REAL(kind=mm_wp), SAVE :: pifac = (4._mm_wp * mm_pi) / 3._mm_wp flx = rho * pifac * m3 * wsettle(mm_play,mm_temp,mm_zlay,mm_drho,mm_drad) RETURN END FUNCTION get_mass_flux END MODULE MM_CLOUDS