source: trunk/LMDZ.TITAN/libf/muphytitan/mm_clouds.f90 @ 3094

Last change on this file since 3094 was 3090, checked in by slebonnois, 15 months ago

BdeBatz? : Cleans microphysics and makes few corrections for physics
File size: 34.2 KB
Line 
1! Copyright (2013-2015,2017,2022-2023) Université de Reims Champagne-Ardenne
2! Contributors : J. Burgalat (GSMA, URCA), B. de Batz de Trenquelléon (GSMA, URCA)
3! email of the author : jeremie.burgalat@univ-reims.fr
4!
5! This software is a computer program whose purpose is to compute
6! microphysics processes using a two-moments scheme.
7!
8! This library is governed by the CeCILL-B license under French law and
9! abiding by the rules of distribution of free software.  You can  use,
10! modify and/ or redistribute the software under the terms of the CeCILL-B
11! license as circulated by CEA, CNRS and INRIA at the following URL
12! "http://www.cecill.info".
13!
14! As a counterpart to the access to the source code and  rights to copy,
15! modify and redistribute granted by the license, users are provided only
16! with a limited warranty  and the software's author,  the holder of the
17! economic rights,  and the successive licensors  have only  limited
18! liability.
19!
20! In this respect, the user's attention is drawn to the risks associated
21! with loading,  using,  modifying and/or developing or reproducing the
22! software by the user in light of its specific status of free software,
23! that may mean  that it is complicated to manipulate,  and  that  also
24! therefore means  that it is reserved for developers  and  experienced
25! professionals having in-depth computer knowledge. Users are therefore
26! encouraged to load and test the software's suitability as regards their
27! requirements in conditions enabling the security of their systems and/or
28! data to be ensured and,  more generally, to use and operate it in the
29! same conditions as regards security.
30!
31! The fact that you are presently reading this means that you have had
32! knowledge of the CeCILL-B license and that you accept its terms.
33
34!! file: mm_clouds.f90
35!! summary: Clouds microphysics module
36!! author: J. Burgalat
37!! date: 2013-2015,2017,2022-2023
38!! modifications: B. de Batz de Trenquelléon
39
40MODULE MM_CLOUDS
41  !! Clouds microphysics module.
42  !!
43  !! The module contains all definitions of the microphysics processes related to clouds:
44  !!
45  !! - [nucleation](page/clouds.html#nucleation)
46  !! - [condensation](page/clouds.html#condensation)
47  !! - [sedimentation](page/clouds.html#sedimentation)
48  !!
49  !!
50  !! The interface methods always use the global variables defined in [[mm_globals(module)]] when values
51  !! (temperature, pressure, moments...) over the vertical grid are required.
52  !! Consequently, all these functions only deal with output arguments which are most of the time the
53  !! tendencies of relevant variables on the atmospheric column.
54  !!
55  !! @note
56  !! Tendencies returned by public methods are always defined from  __TOP__ of the atmosphere to the
57  !! __GROUND__.
58  USE MM_MPREC
59  USE MM_GLOBALS
60  USE MM_METHODS
61  IMPLICIT NONE
62
63  PRIVATE
64
65  PUBLIC :: mm_cloud_microphysics, mm_cloud_sedimentation, mm_cloud_nucond
66
67CONTAINS
68
69  !============================================================================
70  ! CLOUDS MICROPHYSICS INTERFACE SUBROUTINE
71  !============================================================================
72
73  SUBROUTINE mm_cloud_microphysics(dm0a,dm3a,dm0n,dm3n,dm3i,dgazs)
74    !! Get the evolution of moments tracers through clouds microphysics processes.
75    !!
76    !! The subroutine is a wrapper to the clouds microphysics methods. It computes the tendencies of moments
77    !! tracers for nucleation, condensation and sedimentation processes for the atmospheric column.
78    !!
79    !! @note
80    !! Both __dm3i__ and __dgazs__ are 2D-array with the vertical layers in first dimension and the number
81    !! of ice components in the second.
82    REAL(kind=mm_wp), DIMENSION(:), INTENT(out)   :: dm0a
83    !! Tendency of the 0th order moment of the aerosols (fractal mode) (\(m^{-3}\)).
84    REAL(kind=mm_wp), DIMENSION(:), INTENT(out)   :: dm3a
85    !! Tendency of the 3rd order moment of the aerosols distribution (fractal mode) (\(m^{3}.m^{-3}\)) .
86    REAL(kind=mm_wp), DIMENSION(:), INTENT(out)   :: dm0n
87    !! Tendency of the 0th order moment of the ccn distribution (\(m^{-3}\)).
88    REAL(kind=mm_wp), DIMENSION(:), INTENT(out)   :: dm3n
89    !! Tendency of the 3rd order moment of the ccn distribution (\(m^{3}.m^{-3}\)).
90    REAL(kind=mm_wp), DIMENSION(:,:), INTENT(out) :: dm3i
91    !! Tendencies of the 3rd order moments of each ice components (\(m^{3}.m^{-3}\)).
92    REAL(kind=mm_wp), DIMENSION(:,:), INTENT(out) :: dgazs
93    !! Tendencies of each condensible gaz species (\(mol.mol^{-1}\)).
94    REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE   :: zdm0n,zdm3n
95    REAL(kind=mm_wp), DIMENSION(:,:), ALLOCATABLE :: zdm3i
96    INTEGER                                    :: i
97    dm0a = 0._mm_wp ; dm3a = 0._mm_wp
98    dm0n = 0._mm_wp ; dm3n = 0._mm_wp
99    dm3i = 0._mm_wp ; dgazs = 0._mm_wp
100
101    IF (mm_w_cloud_nucond) THEN
102      ! Calls condensation/nucleation (and update saturation ratio diagnostic)
103      call mm_cloud_nucond(dm0a,dm3a,dm0n,dm3n,dm3i,dgazs,mm_gazs_sat)
104    ENDIF
105
106    IF (mm_w_cloud_sed) THEN
107      ! Calls sedimentation
108      ALLOCATE(zdm0n(mm_nla),zdm3n(mm_nla),zdm3i(mm_nla,mm_nesp))
109      call mm_cloud_sedimentation(zdm0n,zdm3n,zdm3i)
110
111      ! Computes settling velocity [m.s-1] of clouds (ccn and ices)
112      mm_ccn_vsed(:) = wsettle(mm_play,mm_temp,mm_zlay,mm_drho,mm_drad)
113
114      ! Computes flux [kg.m-2.s-1] and precipitation [m.iphysiq] of ccn
115      mm_ccn_flux(:) = get_mass_flux(mm_rhoaer,mm_m3ccn(:))
116      mm_ccn_prec = SUM(zdm3n*mm_dzlev)
117
118      ! Computes flux [kg.m-2.s-1] and precipitation [m.iphysiq] of ices
119      DO i = 1, mm_nesp
120        mm_ice_fluxes(:,i) = get_mass_flux(mm_xESPS(i)%rho,(3._mm_wp*mm_m3ice(:,i))/(4._mm_wp*mm_pi))
121        mm_ice_prec(i) = SUM(zdm3i(:,i)*mm_dzlev)
122      ENDDO
123
124      ! updates tendencies
125      dm0n = dm0n + zdm0n
126      dm3n = dm3n + zdm3n
127      dm3i = dm3i + zdm3i
128    ENDIF
129
130  END SUBROUTINE mm_cloud_microphysics
131
132  !-----------------------------------------------------------------------------
133  ! NUCLEATION/CONDENSATION PROCESS RELATED METHODS
134  !-----------------------------------------------------------------------------
135
136  SUBROUTINE mm_cloud_nucond(dm0a,dm3a,dm0n,dm3n,dm3i,dgazs,gazsat)
137    !! Get moments tendencies through nucleation/condensation/evaporation.
138    !!
139    !! The method is a wrapper of [[mm_clouds(module):nc_esp(subroutine)]] which computes the
140    !! tendencies of tracers for all the condensible species given in the vector __xESPS__.
141    !!
142    !! Aerosols and CCN distribution evolution depends on the ice components:
143    !!   - For nucleation only creation of CCN can occur.
144    !!   - For condensation only loss of CCN can occur.
145    !!
146    !! We use the simple following rule to compute the variation of CCN and aerosols:
147    !! The global variation of CCN (and thus aerosols) is determined from the most intense activity
148    !! of the ice components.
149    !!
150    !! @warning
151    !! __xESPS__, __m3i__ and __gazes__ must share the same indexing. For example if __xESPS(IDX)__
152    !! corresponds to \(CH_{4}\) properties then __m3i(IDX)__ must be the total volume of solid
153    !! \(CH_{4}\) (ice) and  __gazs(IDX)__ its vapor mole fraction.
154    REAL(kind=mm_wp), DIMENSION(:), INTENT(out)   :: dm0a
155    !! Tendency of the 0th order moment of the aerosols (fractal mode) (\(m^{-3}\)).
156    REAL(kind=mm_wp), DIMENSION(:), INTENT(out)   :: dm3a
157    !! Tendency of the 3rd order moment of the aerosols distribution (fractal mode) (\(m^{3}.m^{-3}\)).
158    REAL(kind=mm_wp), DIMENSION(:), INTENT(out)   :: dm0n
159    !! Tendency of the 0th order moment of the aerosols distribution (fractal mode) (\(m^{-3}\)).
160    REAL(kind=mm_wp), DIMENSION(:), INTENT(out)   :: dm3n
161    !! Tendency of the 3rd order moment of the ccn distribution (\(m^{3}.m^{-3}\)).
162    REAL(kind=mm_wp), DIMENSION(:,:), INTENT(out) :: dm3i
163    !! Tendencies of the 3rd order moments of each ice components (\(m^{3}.m^{-3}\)).
164    REAL(kind=mm_wp), DIMENSION(:,:), INTENT(out) :: dgazs
165    !! Tendencies of each condensible gaz species (\(mol.mol^{-1}\)) .
166    REAL(kind=mm_wp), DIMENSION(:,:), INTENT(out) :: gazsat
167    !! Saturation ratio of each condensible specie.
168    INTEGER                                       :: i,idx
169    REAL(kind=mm_wp), DIMENSION(:,:), ALLOCATABLE :: zdm0a,zdm3a,zdm0n,zdm3n
170
171    ALLOCATE(zdm0a(mm_nla,mm_nesp),zdm3a(mm_nla,mm_nesp), &
172             zdm0n(mm_nla,mm_nesp),zdm3n(mm_nla,mm_nesp))
173    zdm0a(:,:) = 0._mm_wp ; zdm3a(:,:) = 0._mm_wp
174    zdm0n(:,:) = 0._mm_wp ; zdm3n(:,:) = 0._mm_wp
175    DO i = 1, mm_nesp
176      call nc_esp(mm_xESPS(i),mm_gazs(:,i),mm_m3ice(:,i),dgazs(:,i),dm3i(:,i), &
177                  zdm0a(:,i),zdm3a(:,i),zdm0n(:,i),zdm3n(:,i),gazsat(:,i))
178    ENDDO
179
180    DO i=1, mm_nla
181      ! retrieve the index of the maximum tendency of CCN where ice variation is not null.
182      idx = MAXLOC(zdm0n(i,:),DIM=1,MASK=(dm3i(i,:) /= 0._mm_wp .OR. mm_m3ice(i,:)+dm3i(i,:) >= 0._mm_wp))
183      IF (idx == 0) THEN
184        dm0n(i) = 0._mm_wp
185        dm3n(i) = 0._mm_wp
186        dm0a(i) = 0._mm_wp
187        dm3a(i) = 0._mm_wp
188      ELSE
189        IF (mm_debug .AND. ABS(zdm0n(i,idx)) > 1e3) THEN
190          WRITE(*,'((a),I2.2,(a),ES10.3,(a))') "Z(",i,") = ",mm_play(i)*1e2, &
191            " mbar: Max aer/ccn exchange variation due to specie: "//TRIM(mm_xESPS(idx)%name)
192        ENDIF
193        dm0n(i) = zdm0n(i,idx)
194        dm3n(i) = zdm3n(i,idx)
195        dm0a(i) = zdm0a(i,idx)
196        dm3a(i) = zdm3a(i,idx)
197      ENDIF
198    ENDDO
199
200  END SUBROUTINE mm_cloud_nucond
201
202  SUBROUTINE nc_esp(xESP,vapX,m3iX,dvapX,dm3iX,dm0aer,dm3aer,dm0ccn,dm3ccn,Xsat)
203    !! Get moments tendencies through nucleation/condensation/evaporation of a given condensible specie.
204    !!
205    !! The method computes the global tendencies of the aerosols, ccn and "ice" moments through cloud
206    !! microphysics processes (nucleation & condensation).
207    !!
208    !! @warning
209    !! Input quantities __m3iX__,__m3iO__, __m0aer__,__m3aer__, __m0ccn__,__m3ccn__ are assumed to be in
210    !! \(X.m^{-3}\) (where X is the unit of the moment that is, a number for M0 and a volume - \(m^3\)
211    !! for M3) ; __vapX__ must be expressed in term of molar fraction.
212    TYPE(mm_esp), INTENT(in)                   :: xESP
213    !! Condensate specie properties.
214    REAL(kind=mm_wp),INTENT(in), DIMENSION(:)  :: vapX
215    !! Gas specie molar fraction on the vertical grid from __TOP__ to __GROUND__ (\(mol.mol^{-1}\)).
216    REAL(kind=mm_wp),INTENT(in), DIMENSION(:)  :: m3iX
217    !! 3rd order moment of the ice component (\(m^{3}.m^{-3}\)).
218    REAL(kind=mm_wp),INTENT(out), DIMENSION(:) :: dvapX
219    !! Tendency of gas specie (\(mol.mol^{-1}\)).
220    REAL(kind=mm_wp),INTENT(out), DIMENSION(:) :: dm3iX
221    !! Tendency of the 3rd order moment of the ice component (\(m^{3}.m^{-3}\)).
222    REAL(kind=mm_wp),INTENT(out), DIMENSION(:) :: dm0aer
223    !! Tendency of the 0th order moment of the fractal mode distribution (\(m^{-3}\)).
224    REAL(kind=mm_wp),INTENT(out), DIMENSION(:) :: dm3aer
225    !! Tendency of the 3rd order moment of the fractal mode distribution (\(m^{3}.m^{-3}\)).
226    REAL(kind=mm_wp),INTENT(out), DIMENSION(:) :: dm0ccn
227    !! Tendency of the 0th order moment of the ccn distribution (\(m^{-3}\)).
228    REAL(kind=mm_wp),INTENT(out), DIMENSION(:) :: dm3ccn
229    !! Tendency of the 3rd order moment of the ccn distribution (\(m^{3}.m^{-3}\)).
230    REAL(kind=mm_wp),INTENT(out), DIMENSION(:) :: Xsat
231    !! Saturation ratio values on the vertical grid (--).
232    INTEGER                                 :: i
233    REAL(kind=mm_wp)                        :: bef,aft
234    REAL(kind=mm_wp), DIMENSION(SIZE(vapX)) :: sm0a,sm3a,sm0n,sm3n,sm3iX
235    REAL(kind=mm_wp), DIMENSION(SIZE(vapX)) :: zm0a,zm3a,zm0n,zm3n,zm3iX,zvapX
236    REAL(kind=mm_wp), DIMENSION(SIZE(vapX)) :: pX,sig,qsat,seq,up,down,ctot,newvap,nucr,grate,cm0,cm3,drad
237
238    ! get a copy of drop radius.
239    drad(:) = mm_drad(:)
240
241    ! Initialization :
242    ! Copy input argument and convert units X.m-3 -> X.kg-1
243    ! sXXX is the initial converted value saved
244    sm3iX = m3iX/mm_rhoair
245    sm0a = mm_m0aer_f/mm_rhoair ; sm3a = mm_m3aer_f/mm_rhoair
246    sm0n = mm_m0ccn/mm_rhoair ; sm3n = mm_m3ccn/mm_rhoair
247    ! zXXX is our working copy
248    zm3iX = sm3iX ; zm0a = sm0a ; zm3a = sm3a ; zm0n = sm0n ; zm3n = sm3n
249
250    ! Molar fraction of X specie is set in mass mixing ratio [kg.kg-1]
251    zvapX  = vapX  * xESP%fmol2fmas
252    ! Surface tension [N.m-1]
253    sig = mm_sigX(mm_temp,xESP)
254    ! X specie mass mixing ratio at saturation [kg.kg-1]
255    qsat = mm_ysatX(mm_temp,mm_play,xESP) * xESP%fmol2fmas
256    ! partial pressure of X specie
257    pX = vapX * mm_play
258    ! Saturation ratio
259    Xsat = zvapX / qsat
260   
261    ! Gets nucleation rate (ccn radius is the monomer !)
262    call nuc_rate((/(mm_rm, i=1,mm_nla)/),mm_temp,xESP,pX,Xsat,nucr)
263    ! IMPORTANT: update CCN and aerosols moment from nucleation NOW !
264    ! 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).
265    !
266    ! 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)
267    ! 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.
268    !
269    ! for nucleation we have the following equations:
270    !   dMa(k)/dt = - dMn(k)/dt (conservation of aerosols+ccn)           (1)
271    !   dMa(k)/dt = - 4*PI*nucr/rm * Ma(k+3)                             (2)
272    !             = - 4*PI*nucr/rm * alpha(k+3)/alpha(k) * rc**3 * Ma(k)
273    ! With:
274    !   - Ma(k): k-th order moment of aerosol
275    !   - Mn(k): k-th order moment of ccn
276    !   - nucr : the nucleation rate.
277    ! We solve (implicit scheme) :
278    !   CST_M(k) = 4*PI*nucr/rm * alpha(k+3)/alpha(k)*rc**3 * dt
279    !   Ma(k)[t+dt] = 1/(1+CST_M(k)) * Ma(k)[t]                          (3)
280    ! Then, from eq. 2:
281    ! Mn(k)[t+dt] = Mn(k)[t] + CST_M(k)/(1+CST_M(k))*Ma(k)[t]            (4)
282    cm0 = 4._mm_wp*mm_pi*nucr/mm_rm*mm_alpha_f(3._mm_wp)*mm_rcf**3*mm_dt
283    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
284    zm0a = 1._mm_wp/(1._mm_wp+cm0) * zm0a
285    zm3a = 1._mm_wp/(1._mm_wp+cm3) * zm3a
286    WHERE (zm0a <= 0._mm_wp .OR. zm3a <= 0._mm_wp)
287      zm0a = 0._mm_wp
288      zm3a = 0._mm_wp
289      zm0n = zm0n + sm0a
290      zm3n = zm3n + sm3a
291    ELSEWHERE
292      zm0n = zm0n + cm0*zm0a
293      zm3n = zm3n + cm3*zm3a
294    ENDWHERE
295
296    ! update the drop radius (we probably should recompute totally the radius to be in better agreement with the other moments)
297    ! We must manage the case where there is no ices and no ccn ==> drop radius is ZERO,
298    ! but conditions are met to spawn nucleation process: creation of ccn.
299    ! Then we set the drop radius to the monomer radius.
300    !
301    ! Doing so will prevent a nasty bug to occur later when ice volume is updated !
302    WHERE (nucr > 0._mm_wp .AND. drad <= mm_drad_min)
303      drad = mm_rm
304    ENDWHERE
305
306    ! Equilibrium saturation near the drop
307    seq = exp(min(30._mm_wp,2._mm_wp*sig*xESP%masmol/(xESP%rho*mm_rgas*mm_temp*drad)))
308    ! Gets growth rate
309    call growth_rate(mm_temp,mm_play,pX/Xsat,xESP,seq,drad,grate)
310    ctot = zvapx + xESP%rho * zm3iX
311    up = zvapx + mm_dt * grate * 4._mm_wp * mm_pi * xESP%rho * drad * seq * zm0n
312    down = 1._mm_wp + mm_dt * grate * 4._mm_wp * mm_pi * xESP%rho * drad / qsat * zm0n
313    ! gets new vapor X specie mass mixing ratio : cannot be greater than the
314    ! total gas + ice and lower than nothing :)
315    newvap = max(min(up/down,ctot),0._mm_wp)
316    ! gets "true" growth rate
317    grate = grate * (newvap/qsat - seq)
318
319    ! computes tendencies through condensation
320    DO i=1,mm_nla
321      ! check for the specific case : NO ICE and SUBLIMATION
322      IF (zm3iX(i) <= 0._mm_wp .AND. grate(i) <= 0._mm_wp) THEN
323        zm3iX(i) = 0._mm_wp
324      ELSE
325        ! update ice volume ...
326        zm3iX(i) = zm3iX(i) + mm_dt*grate(i)*4._mm_wp*mm_pi*drad(i)*zm0n(i)
327        ! ... and check if there is ice left in the ccn
328        IF (zm3iX(i) <= 0._mm_wp) THEN
329          zm3iX(i) = 0._mm_wp
330          zm0a(i) = zm0a(i) + zm0n(i) ; zm0n(i) = 0._mm_wp
331          zm3a(i) = zm3a(i) + zm3n(i) ; zm3n(i) = 0._mm_wp
332        ENDIF
333      ENDIF
334    ENDDO
335
336    ! Computes balance
337    IF (mm_debug) THEN
338      DO i=1,mm_nla
339        bef = sm0a(i) + sm0n(i)
340        aft = zm0a(i)  + zm0n(i)
341        IF (ABS(bef-aft)/bef > 1e-10_mm_wp) THEN
342          WRITE(*,'((a),I2.2,(a),ES20.12,(a),ES20.12)') &
343            "[DEBUG] nc_esp("//TRIM(xESP%name)//"): M0 not conserved (z=",i,")",bef," <-> ",aft
344        ENDIF
345        bef = sm3a(i) + sm3n(i)
346        aft = zm3a(i)  + zm3n(i)
347        IF (ABS(bef-aft)/bef > 1e-10_mm_wp) THEN
348          WRITE(*,'((a),I2.2,(a),ES20.12,(a),ES20.12)') &
349            "[DEBUG] nc_esp("//TRIM(xESP%name)//"): M3 not conserved (z=",i,")",bef," <-> ",aft
350        ENDIF
351      ENDDO
352    ENDIF
353
354    ! compute tendencies:
355    ! all of these tendencies are in X.m-3 !
356    dm0aer = (zm0a - sm0a)*mm_rhoair
357    dm3aer = (zm3a - sm3a)*mm_rhoair
358    dm0ccn = (zm0n - sm0n)*mm_rhoair
359    dm3ccn = (zm3n - sm3n)*mm_rhoair
360
361    dm3iX  = (zm3iX - sm3iX)                    ! this one in X.kg-1 (temporary) !
362    dvapX  = -xESP%rho * dm3iX / xESP%fmol2fmas ! in order to compute this one in mol.mol-1
363    dm3iX  = dm3iX*mm_rhoair                    ! update ice tendencies in X.m-3 !
364
365  END SUBROUTINE nc_esp
366
367  SUBROUTINE nuc_rate(rccn,temp,xESP,pvp,sat,rate)
368    !! Get nucleation rate.
369    !!
370    !! The method computes the heterogeneous nucleation rate for the given specie on a fractal particle
371    !! of size __rccn__.
372    !! Except __xESP__, all arguments are vectors of the same size (vertical grid).
373    REAL(kind=mm_wp), INTENT(in), DIMENSION(:)  :: rccn !! Radius of the cloud condensation nuclei (m).
374    REAL(kind=mm_wp), INTENT(in), DIMENSION(:)  :: temp !! Temperature (K).
375    REAL(kind=mm_wp), INTENT(in), DIMENSION(:)  :: pvp  !! Partial vapor pressure of X specie (Pa).
376    REAL(kind=mm_wp), INTENT(in), DIMENSION(:)  :: sat  !! Saturation ratio of given specie (--).
377    TYPE(mm_esp), INTENT(in)                    :: xESP !! X specie properties (--).
378    REAL(kind=mm_wp), INTENT(out), DIMENSION(:) :: rate !! The nucleation rate (\(m^-{2}.s^{-1}\)).
379    INTEGER          :: i
380    REAL(kind=mm_wp) :: r,t,s,sig,nX,rstar,gstar,x,zeldov,deltaf,fsh,fstar
381
382    rate(:) = 0._mm_wp
383    ! Activation condition
384    DO i=1, SIZE(rccn)
385      s = sat(i)
386      IF (s > 1._mm_wp) THEN
387        t = temp(i) ; r = rccn(i)
388        sig = mm_sigX(t,xESP)
389        nX    = pvp(i)/mm_kboltz/t
390        rstar = 2._mm_wp*sig*xESP%vol/(mm_kboltz*t*dlog(s))
391        ! curvature radius
392        x     = r/rstar
393        fsh   = mm_fshape(xESP%mteta,x)
394        fstar = (4._mm_wp*mm_pi/3._mm_wp)*sig*(rstar**2.)*fsh
395        deltaf=MIN(MAX((2.*mm_fdes-mm_fdif-fstar)/(mm_kboltz*t),-100._mm_wp),100._mm_wp)
396        IF (deltaf > -100._mm_wp) THEN
397          gstar  = 4._mm_wp*mm_pi*(rstar**3)/(3._mm_wp*xESP%vol)
398          zeldov = dsqrt(fstar/(3._mm_wp*mm_pi*mm_kboltz*t*(gstar**2)))
399          rate(i)= zeldov*mm_kboltz*t*(nX*rstar)**2._mm_wp*dexp(deltaf)/(fsh*mm_nus*xESP%mas)
400        ENDIF
401      ENDIF
402    ENDDO
403
404    RETURN
405  END SUBROUTINE nuc_rate
406
407  SUBROUTINE growth_rate(temp,pres,pXsat,xESP,seq,drad,rate)
408    !! Get growth rate through condensation/evaporation process.
409    !!
410    !! The method computes the growth rate a drop through condensation/evaporation processes:
411    !!
412    !! $$ r \times \frac{dr}{dt} = g_{rate} \times (S - S_{eq}) $$
413    !!
414    !! Except __xESP__ which is a scalar, all arguments are vectors of the same size (vertical grid).
415    REAL(kind=mm_wp), INTENT(in), DIMENSION(:)  :: temp  !! Temperature (K).
416    REAL(kind=mm_wp), INTENT(in), DIMENSION(:)  :: pres  !! Pressure level (Pa).
417    REAL(kind=mm_wp), INTENT(in), DIMENSION(:)  :: pXsat !! Saturation vapor pressure of specie (Pa).
418    TYPE(mm_esp), INTENT(in)                    :: xESP  !! Specie properties.
419    REAL(kind=mm_wp), INTENT(in), DIMENSION(:)  :: seq   !! Equilibrium saturation near the drop.
420    REAL(kind=mm_wp), INTENT(in), DIMENSION(:)  :: drad  !! Drop radius (m).
421    REAL(kind=mm_wp), INTENT(out), DIMENSION(:) :: rate  !! Growth rate (\(m^{2}.s^{-1}\)).
422    REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE :: k,knu,slf,rkc,rdc,l,dv
423    INTEGER                                     :: n
424    n = SIZE(temp)
425    ALLOCATE(k(n),knu(n),slf(n),rkc(n),rdc(n),l(n),dv(n))
426
427    ! N2 (air) Thermal conductivity (where does it come from ?)
428    k(:) = ( 2.857e-2_mm_wp*temp-0.5428_mm_wp)*4.184e-3_mm_wp
429    ! Gas mean free path
430    l(:) = mm_lambda_g(temp,pres)
431    ! Diffusion coefficient of X gas
432    Dv(:) = 1._mm_wp/3._mm_wp*dsqrt(8._mm_wp*mm_rgas*temp(:)/(mm_pi*xESP%masmol))*mm_kboltz*temp(:) / &
433      (mm_pi*pres(:)*(mm_air_rad+xESP%ray)**2* dsqrt(1._mm_wp+xESP%fmol2fmas))
434    knu(:) = l(:)/drad(:)                                         ! The knudsen number of the drop
435    slf(:) = (1.333_mm_wp+0.71_mm_wp/knu(:))/(1._mm_wp+1._mm_wp/knu(:)) ! Slip flow correction
436    Dv(:)  = Dv(:)/(1._mm_wp+slf(:)*knu(:))
437    ! latent heat resistance coefficient
438    rkc(:) = mm_lheatX(temp(:),xESP)**2 * xESP%rho * xESP%masmol / (k(:)*mm_rgas*temp(:)**2)
439    ! Diffusion resistance coefficient
440    rdc(:) = mm_rgas * temp(:) * xESP%rho / (Dv(:)*pXsat(:)*xESP%masmol)
441    ! Growth rate: rdr/dt = rate * (S-Seq) ; rate is returned
442    rate(:) = 1._mm_wp / (seq(:)*rkc(:)+rdc(:))
443    RETURN
444  END SUBROUTINE growth_rate
445
446
447  !-----------------------------------------------------------------------------
448  ! SEDIMENTATION PROCESS RELATED METHODS
449  !-----------------------------------------------------------------------------
450
451  SUBROUTINE mm_cloud_sedimentation(dm0n,dm3n,dm3i)
452    !! Compute the tendency of _clouds_ related moments through sedimentation process.
453    !!
454    !! The method computes the tendencies of moments related to cloud microphysics through
455    !! sedimentation process. The algorithm used here differs from
456    !! [[mm_haze(module):mm_haze_sedimentation(subroutine)]] as all moments settle with the same
457    !! terminal velocity which is computed with the average drop radius of the size distribution.
458    !! We simply compute an _exchange matrix_ that stores the new positions of each cells through
459    !! sedimentation process and then computes the matrix
460    !! product with input moments values to get final tendencies.
461    REAL(kind=mm_wp), INTENT(out), DIMENSION(:)   :: dm0n
462    !! Tendency of the 0th order moment of the ccn distribution (\(m^{-3}\)).
463    REAL(kind=mm_wp), INTENT(out), DIMENSION(:)   :: dm3n
464    !! Tendency of the 3rd order moment of the ccn distribution (\(m^{3}.m^{-3}\)).
465    REAL(kind=mm_wp), INTENT(out), DIMENSION(:,:) :: dm3i
466    !! Tendencies of the 3rd order moment of each ice component of the cloud (\(m^{3}m^{-3}\)).
467    INTEGER                                   :: im,nm
468    REAL(kind=mm_wp), DIMENSION(:,:), ALLOCATABLE :: moms, momsf,chg_matrix
469
470    nm = 2 + mm_nesp
471    ALLOCATE(moms(mm_nla,nm),momsf(mm_nla,nm),chg_matrix(mm_nla,mm_nla))
472    ! Initializes moms
473    moms(:,1)  = mm_m0ccn * mm_dzlev
474    moms(:,2)  = mm_m3ccn * mm_dzlev
475    DO im=1,mm_nesp ; moms(:,2+im) = mm_m3ice(:,im) * mm_dzlev ; ENDDO
476    ! Computes exchange matrix
477    CALL exchange(mm_drad,mm_drho,mm_dt,chg_matrix)
478    ! Computes final moments values
479    DO im=1,nm ; momsf(:,im) = MATMUL(chg_matrix,moms(:,im)) ; ENDDO
480    ! Computes tendencies (converted in X.m-3)
481    dm0n = (momsf(:,1)-moms(:,1))/mm_dzlev
482    dm3n = (momsf(:,2)-moms(:,2))/mm_dzlev
483    DO im=1,mm_nesp ; dm3i(:,im) = (momsf(:,2+im)-moms(:,2+im))/mm_dzlev ; ENDDO
484    RETURN
485  END SUBROUTINE mm_cloud_sedimentation
486
487  SUBROUTINE exchange(rad,rhog,dt,matrix)
488    !! Compute the exchange matrix.
489    !!
490    !! The subroutine computes the matrix exchange used by
491    !! [[mm_clouds(module):mm_cloud_sedimentation(subroutine)]] to compute moments tendencies
492    !! through sedimentation process. Both __rad__ and __rhog__ must be vector with relevant
493    !! values over the atmospheric vertical structure. __matrix__ is square 2D-array with same
494    !! dimension size than __rad__.
495    REAL(kind=mm_wp), INTENT(in), DIMENSION(:) :: rad
496    !! Cloud drop radius over the atmospheric vertical structure (m).
497    REAL(kind=mm_wp), INTENT(in), DIMENSION(:) :: rhog
498    !! Cloud drop density over the atmospheric vertical structure (\(kg.m^{-3}\)).
499    REAL(kind=mm_wp), INTENT(in)               :: dt
500    !! Timestep (s).
501    REAL(kind=mm_wp), INTENT(out)              :: matrix(:,:)
502    !! The output _exchange matrix_.
503    INTEGER                                     :: nz,i,j,jj,jinf,jsup
504    REAL(kind=mm_wp)                            :: zni,znip1,xf,xft,xcnt
505    REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE :: puit
506    REAL(kind=mm_wp)                            :: cpte,cpte2
507    REAL(kind=mm_wp)                            :: zsurf
508    INTEGER, PARAMETER                          :: ichx  = 1
509    matrix = 0._mm_wp ; nz = SIZE(rad) ; ALLOCATE(puit(nz))
510    zsurf = mm_zlev(nz)
511   
512    ! compute exchange matrix
513    DO i=1,nz
514      puit(i) = 0._mm_wp
515      xcnt = 0._mm_wp
516      ! computes drop move (i.e. its new positions)
517      CALL getnzs(ichx,i,rad,rhog,dt,zni,znip1)
518
519      ! Peculiar case : Ground level precipitation [znip1 < zsurf && (zni < zsurf || zni > zsurf)]
520      ! - complete precipitation [ znip1 <= 0 && zni <= 0 ] :
521      IF(zni <= zsurf .and. znip1 <= zsurf) THEN
522        xft=0._mm_wp
523        xf=1._mm_wp
524        xcnt=xcnt+xf
525        puit(i)=puit(i)+xf
526      ENDIF
527      ! - partial precipitation [ znip1 <= zsurf && zni > zsurf ] :
528      IF (zni > zsurf .and. znip1 <= zsurf) THEN
529        xft=(zni-zsurf)/(zni-znip1)
530        xf=(1.-xft)
531        xcnt=xcnt+xf
532        puit(i)=puit(i)+xf
533      ENDIF
534      ! General case : no ground precipitation [ znip1 > zsurf && zni > zsurf ]
535      IF (zni > zsurf .and. znip1 > zsurf) THEN
536        xft = 1._mm_wp ! on a la totalite de la case
537        xf  = 0._mm_wp
538        xcnt=xcnt+xf
539        puit(i)=puit(i)+xf
540      ENDIF
541      IF (zni < znip1) THEN
542        WRITE(*,'("[EXCHANGES] WARNING, missing this case :",2(2X,ES10.3))') zni, znip1
543      ENDIF
544     
545      ! Fix minimum level to the ground
546      znip1 = MAX(znip1,zsurf)
547      zni   = MAX(zni,zsurf)
548      ! Locate new "drop" position in the verical grid
549      jsup=nz+1
550      jinf=nz+1
551      DO j=1,nz
552        IF (zni<=mm_zlev(j).and.zni>=mm_zlev(j+1)) jsup=j
553        IF (znip1<=mm_zlev(j).and.znip1>=mm_zlev(j+1)) jinf=j
554      ENDDO
555      ! Volume is out of range: (all drops have touched the ground!)
556      ! Note: cannot happen here, it has been treated previously :)
557      IF (jsup>=nz+1.and.jinf==jsup) THEN
558        WRITE(*,'(a)') "[EXCHANGE] speaking: The impossible happened !"
559        call EXIT(666)
560      ENDIF
561      ! Volume inside a single level
562      IF (jsup==jinf.and.jsup<=nz) THEN
563        xf=1._mm_wp
564        xcnt=xcnt+xft*xf
565        matrix(jinf,i)=matrix(jinf,i)+xft*xf
566      ENDIF
567
568      ! Volume over 2 levels
569      IF (jinf==jsup+1) THEN
570        xf=(zni-mm_zlev(jinf))/(zni-znip1)
571        xcnt=xcnt+xf*xft
572        IF(jsup<=nz) THEN
573          matrix(jsup,i)=matrix(jsup,i)+xft*xf
574        ENDIF
575        xf=(mm_zlev(jinf)-znip1)/(zni-znip1)
576        xcnt=xcnt+xf*xft
577        IF (jinf<=nz) THEN
578          matrix(jinf,i)=matrix(jinf,i)+xft*xf
579        ENDIF
580      ENDIF
581
582      ! Volume over 3 or more levels
583      IF (jinf > jsup+1) THEN
584        ! first cell
585        xf=(zni-mm_zlev(jsup+1))/(zni-znip1)
586        xcnt=xcnt+xf*xft
587        matrix(jsup,i)=matrix(jsup,i)+xft*xf
588        ! last cell
589        xf=(mm_zlev(jinf)-znip1)/(zni-znip1)
590        xcnt=xcnt+xf*xft
591        matrix(jinf,i)=matrix(jinf,i)+xft*xf
592        ! other :)
593        DO jj=jsup+1,jinf-1
594          xf=(mm_zlev(jj)-mm_zlev(jj+1))/(zni-znip1)
595          xcnt=xcnt+xf*xft
596          matrix(jj,i)=matrix(jj,i)+xft*xf
597        ENDDO
598      ENDIF
599    ENDDO
600
601    ! checking if everything is alright if debug enabled...
602    IF (mm_debug) THEN
603      cpte=0._mm_wp ; cpte2=0._mm_wp
604      DO j=1,nz
605        DO jj=1,nz
606          cpte=cpte+matrix(jj,j)
607        ENDDO
608      ENDDO
609      cpte2=cpte+sum(puit)
610      IF (abs(cpte2-nz)>1.e-4_mm_wp) THEN
611        WRITE(*,'("[EXCHANGE] speaking: tx expl (/nz):",2(2X,ES10.3))') cpte,cpte2
612      ENDIF
613    ENDIF
614
615    RETURN
616  END SUBROUTINE exchange
617
618  SUBROUTINE getnzs(ichx,idx,rad,rho,dt,zni,zns)
619    !! Compute displacement of a cell under sedimentation process.
620    !!
621    !! The method computes the new position of a _drop cell_ through sedimentation process as
622    !! descibed in the following scheme:
623    !!
624    !! ![Cloud sedimentation scheme](|media|/cloud_sed_scheme.svg)
625    !!
626    !! New positions are returned in __zni__ and __zns__ ouptut arguments.
627    !!
628    !! @note
629    !! The method uses directly [[mm_globals(module):mm_play(variable)]], [[mm_globals(module):mm_plev(variable)]],
630    !! [[mm_globals(module):mm_temp(variable)]],[[mm_globals(module):mm_btemp(variable)]],
631    !! [[mm_globals(module):mm_zlay(variable)]] and [[mm_globals(module):mm_zlev(variable)]] and uses __idx__ to
632    !! get the relevant value to use on the vertical grid.
633    INTEGER, INTENT(in)                        :: ichx
634    !! Velocity extrapolation control flag (0 for linear, 1 for exponential -preferred -).
635    INTEGER, INTENT(in)                        :: idx
636    !! Initial position of the drop (subscript of vertical layers vectors).
637    REAL(kind=mm_wp), INTENT(in), DIMENSION(:) :: rad
638    !! Cloud drop radius over the atmospheric vertical structure (m).
639    REAL(kind=mm_wp), INTENT(in), DIMENSION(:) :: rho
640    !! Cloud drop density over the atmospheric vertical structure (\(kg.m^{-3}\)).
641    REAL(kind=mm_wp), INTENT(in)               :: dt
642    !! Timestep (s).
643    REAL(kind=mm_wp), INTENT(out)              :: zni
644    !! Final layer upper boundary position (m).
645    REAL(kind=mm_wp), INTENT(out)              :: zns
646    !! Final layer lower boundary position (m).
647    REAL(kind=mm_wp)            :: ws,wi,w,zi,zs
648    REAL(kind=mm_wp)            :: alpha,argexp,v0,arg1,arg2
649    INTEGER                 :: i,nz
650    REAL(kind=mm_wp), PARAMETER :: es = 30._mm_wp
651
652    nz = SIZE(rad)
653    ! Linear extrapolation of velocity
654    IF (ichx==0) THEN
655      ! velocity upper interface
656      ws = wsettle(mm_plev(idx),mm_btemp(idx),mm_zlev(idx),rho(idx),rad(idx))
657      IF (idx==nz) THEN
658        ! veloctity center layer
659        wi = wsettle(mm_play(idx),mm_temp(idx),mm_zlay(idx),rho(idx),rad(idx))
660      ELSEIF(idx<nz) THEN
661        ! velocity lower interface
662        wi = wsettle(mm_plev(idx+1),mm_btemp(idx+1),mm_zlev(idx+1), &
663          rho(idx+1),rad(idx+1))
664      ELSE
665        WRITE(*,'(a)') "[getnzs] speaking:"
666        WRITE(*,'(a)') "This is the fatal error..."
667        WRITE(*,'(a)') "index is higher than number of levels"
668        call EXIT(111)
669      ENDIF
670      w   = (ws+wi)/2._mm_wp
671      zni = mm_zlev(idx)-w*dt
672      zns = mm_zlev(idx)-mm_dzlev(idx)-w*dt
673      RETURN
674      ! Exponential extrapolation of velocity
675
676    ELSEIF(ichx==1) THEN
677      zs = mm_zlev(idx)
678      ws = wsettle(mm_plev(idx),mm_btemp(idx),zs,rho(idx),rad(idx))
679      zi=mm_zlay(idx)
680      wi = wsettle(mm_play(idx),mm_temp(idx),zi,rho(idx),rad(idx))
681      ! ws & wi must be different !
682      IF(dabs(wi-ws)/wi <= 1.e-3_mm_wp)  wi=ws/1.001_mm_wp
683      IF (wi /= 0._mm_wp) alpha = dlog(ws/wi)/(zs-zi)  ! alpha < 0 if wi > ws
684      !   -es < argexp < es
685      argexp=MAX(MIN(alpha*zs,es),-es)
686      v0 = ws/dexp(argexp)
687      arg1=1._mm_wp+v0*alpha*dexp(argexp)*dt
688      argexp=MAX(MIN(alpha*(mm_zlev(idx)-mm_dzlev(idx)),es),-es)
689      arg2=1._mm_wp+v0*alpha*dexp(argexp)*dt
690      IF (arg1<=0._mm_wp.OR.arg2<=0._mm_wp) THEN
691        ! correct velocity
692        ! divides the velocity argument in arg1 and arg2 :
693        ! argX=1+alpha*v0*exp(alpha*z)*dt <==> argX-1=alpha*v0*exp(alpha*z)*dt
694        ! ==> argX' = 1 + (argX-1)/2  <==> argX' = (1+argX)/2.
695        DO i=1,25
696          IF (arg1<=0._mm_wp.OR.arg2<=0._mm_wp) THEN
697            IF (mm_debug) &
698              WRITE(*,'((a),I2.2,(a))') "[getnzs] must adjust velocity (iter:",i,"/25)"
699            arg1=(arg1+1._mm_wp)/2._mm_wp ; arg2=(arg2+1._mm_wp)/2._mm_wp
700          ELSE
701            EXIT
702          ENDIF
703        ENDDO
704        ! sh** we have to stop
705        IF (i>25) THEN
706          WRITE(*,'(a)')"[getnzs] speaking:"
707          WRITE(*,'(a)') "Cannot adjust velocities"
708          call EXIT(111)
709        ENDIF
710      ENDIF
711
712      zni = mm_zlev(idx)-dlog(arg1)/alpha
713      zns = mm_zlev(idx)-mm_dzlev(idx)-dlog(arg2)/alpha
714
715      RETURN
716    ENDIF
717  END SUBROUTINE getnzs
718
719  ELEMENTAL FUNCTION wsettle(p,t,z,rho,rad) RESULT(w)
720    !! Compute the settling velocity of a spherical particle.
721    !!
722    !! The method computes the effective settling velocity of spherical particle of
723    !! radius __rad__. It accounts for the slip-flow transition (no approximation).
724    REAL(kind=mm_wp), INTENT(in) :: p   !! The pressure level (Pa).
725    REAL(kind=mm_wp), INTENT(in) :: t   !! The temperature (K).
726    REAL(kind=mm_wp), INTENT(in) :: z   !! The altitude level (m).
727    REAL(kind=mm_wp), INTENT(in) :: rho !! Density of the particle (\(kg.m^{-3}\)).
728    REAL(kind=mm_wp), INTENT(in) :: rad !! Radius of the particle (m).
729    REAL(kind=mm_wp) :: w               !! Settling velocity (\(m.s^{-1}\)).
730    REAL(kind=mm_wp) :: Us, Fc, kn
731    REAL(kind=mm_wp), PARAMETER :: ra = 1.75e-10_mm_wp
732   
733    ! Knudsen number
734    kn = (mm_kboltz * t) / (p * 4._mm_wp * sqrt(2._mm_wp) * mm_pi * ra**2) / rad
735
736    ! Computes slip-flow correction : Fc = 1 + (mm_akn * mm_lambda_g(t,p) / rad)
737    Fc = (1._mm_wp + 1.2517_mm_wp*kn + 0.4_mm_wp*kn*dexp(-1.1_mm_wp/kn))
738   
739    ! Computes Stokes settling velocity
740    Us = (2._mm_wp * rad**2 * rho * mm_effg(z)) / (9._mm_wp * mm_eta_g(t))
741   
742    ! Computes settling velocity (correction factor : x2.0)
743    w = Us * Fc * 2._mm_wp
744  END FUNCTION wsettle
745
746  FUNCTION get_mass_flux(rho,m3) RESULT(flx)
747    !> Get the mass flux of (clouds related) moment through sedimention.
748    !!
749    !! @warning
750    !! The method is __only__ valid for cloud moments (i.e. ice or ccn). It calls
751    !! [[mm_clouds(module):wsettle(function)]] that compute the _mean_ settling velocity of a
752    !! cloud drop.
753    !!
754    !! @note
755    !! The computed flux is always positive.
756    REAL(kind=mm_wp), INTENT(in)               :: rho
757    !! Tracer density (\(kg.m^{-3}\)).
758    REAL(kind=mm_wp), INTENT(in), DIMENSION(:) :: m3
759    !! Vertical profile of the total volume of tracer (i.e. M3) from __TOP__ to __GROUND__ (\(m^{3}.m^{-3}\)).
760    REAL(kind=mm_wp), DIMENSION(SIZE(m3))      :: flx
761    !! Mass sedimentation fluxes at each layer from __TOP__ to __GROUND__ (\(kg.m^{-2}.s^{-1}\)).
762    REAL(kind=mm_wp), SAVE :: pifac = (4._mm_wp * mm_pi) / 3._mm_wp
763    flx = rho * pifac * m3 * wsettle(mm_play,mm_temp,mm_zlay,mm_drho,mm_drad)
764    RETURN
765  END FUNCTION get_mass_flux
766
767END MODULE MM_CLOUDS
Note: See TracBrowser for help on using the repository browser.