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

Last change on this file since 3496 was 3496, checked in by debatzbr, 2 weeks ago

Final clean and debug for the microphysical model
File size: 34.4 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 [kg.m-2] of ccn
115      mm_ccn_flux(:) = get_mass_flux(mm_rhoaer,mm_m3ccn(:))
116      mm_ccn_prec = SUM(zdm3n*mm_dzlev*mm_rhoaer)
117
118      ! Computes flux [kg.m-2.s-1] and precipitation [kg.m-2] 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*mm_xESPS(i)%rho)
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   
262    ! Gets nucleation rate (ccn radius is the monomer !)
263    call nuc_rate((/(mm_rm, i=1,mm_nla)/),mm_temp,xESP,pX,Xsat,nucr)
264    ! IMPORTANT: update CCN and aerosols moment from nucleation NOW !
265    ! 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).
266    !
267    ! 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)
268    ! 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.
269    !
270    ! for nucleation we have the following equations:
271    !   dMa(k)/dt = - dMn(k)/dt (conservation of aerosols+ccn)           (1)
272    !   dMa(k)/dt = - 4*PI*nucr/rm * Ma(k+3)                             (2)
273    !             = - 4*PI*nucr/rm * alpha(k+3)/alpha(k) * rc**3 * Ma(k)
274    ! With:
275    !   - Ma(k): k-th order moment of aerosol
276    !   - Mn(k): k-th order moment of ccn
277    !   - nucr : the nucleation rate.
278    ! We solve (implicit scheme) :
279    !   CST_M(k) = 4*PI*nucr/rm * alpha(k+3)/alpha(k)*rc**3 * dt
280    !   Ma(k)[t+dt] = 1/(1+CST_M(k)) * Ma(k)[t]                          (3)
281    ! Then, from eq. 2:
282    ! Mn(k)[t+dt] = Mn(k)[t] + CST_M(k)/(1+CST_M(k))*Ma(k)[t]            (4)
283    cm0 = 4._mm_wp*mm_pi*nucr/mm_rm*mm_alpha_f(3._mm_wp)*mm_rcf**3*mm_dt
284    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
285    zm0a = 1._mm_wp/(1._mm_wp+cm0) * zm0a
286    zm3a = 1._mm_wp/(1._mm_wp+cm3) * zm3a
287    WHERE (zm0a <= 0._mm_wp .OR. zm3a <= 0._mm_wp)
288      zm0a = 0._mm_wp
289      zm3a = 0._mm_wp
290      zm0n = zm0n + sm0a
291      zm3n = zm3n + sm3a
292    ELSEWHERE
293      zm0n = zm0n + cm0*zm0a
294      zm3n = zm3n + cm3*zm3a
295    ENDWHERE
296
297    ! update the drop radius (we probably should recompute totally the radius to be in better agreement with the other moments)
298    ! We must manage the case where there is no ices and no ccn ==> drop radius is ZERO,
299    ! but conditions are met to spawn nucleation process: creation of ccn.
300    ! Then we set the drop radius to the monomer radius.
301    !
302    ! Doing so will prevent a nasty bug to occur later when ice volume is updated !
303    WHERE (nucr > 0._mm_wp .AND. drad <= mm_drad_min)
304      drad = mm_rm
305    ENDWHERE
306
307    ! Equilibrium saturation near the drop
308    seq = exp(min(30._mm_wp,2._mm_wp*sig*xESP%masmol/(xESP%rho*mm_rgas*mm_temp*drad)))
309    ! Gets growth rate
310    call growth_rate(mm_temp,mm_play,pX/Xsat,xESP,seq,drad,grate)
311    ctot = zvapx + xESP%rho * zm3iX
312    up = zvapx + mm_dt * grate * 4._mm_wp * mm_pi * xESP%rho * drad * seq * zm0n
313    down = 1._mm_wp + mm_dt * grate * 4._mm_wp * mm_pi * xESP%rho * drad / qsat * zm0n
314    ! gets new vapor X specie mass mixing ratio : cannot be greater than the
315    ! total gas + ice and lower than nothing :)
316    newvap = max(min(up/down,ctot),0._mm_wp)
317    ! gets "true" growth rate
318    grate = grate * (newvap/qsat - seq)
319
320    ! computes tendencies through condensation
321    DO i=1,mm_nla
322      ! check for the specific case : NO ICE and SUBLIMATION
323      IF (zm3iX(i) <= 0._mm_wp .AND. grate(i) <= 0._mm_wp) THEN
324        zm3iX(i) = 0._mm_wp
325      ELSE
326        ! update ice volume ...
327        zm3iX(i) = zm3iX(i) + mm_dt*grate(i)*4._mm_wp*mm_pi*drad(i)*zm0n(i)
328        ! ... and check if there is ice left in the ccn
329        IF (zm3iX(i) <= 0._mm_wp) THEN
330          zm3iX(i) = 0._mm_wp
331          zm0a(i) = zm0a(i) + zm0n(i) ; zm0n(i) = 0._mm_wp
332          zm3a(i) = zm3a(i) + zm3n(i) ; zm3n(i) = 0._mm_wp
333        ENDIF
334      ENDIF
335    ENDDO
336
337    ! Computes balance
338    IF (mm_debug) THEN
339      DO i=1,mm_nla
340        bef = sm0a(i) + sm0n(i)
341        aft = zm0a(i)  + zm0n(i)
342        IF (ABS(bef-aft)/bef > 1e-10_mm_wp) THEN
343          WRITE(*,'((a),I2.2,(a),ES20.12,(a),ES20.12)') &
344            "[DEBUG] nc_esp("//TRIM(xESP%name)//"): M0 not conserved (z=",i,")",bef," <-> ",aft
345        ENDIF
346        bef = sm3a(i) + sm3n(i)
347        aft = zm3a(i)  + zm3n(i)
348        IF (ABS(bef-aft)/bef > 1e-10_mm_wp) THEN
349          WRITE(*,'((a),I2.2,(a),ES20.12,(a),ES20.12)') &
350            "[DEBUG] nc_esp("//TRIM(xESP%name)//"): M3 not conserved (z=",i,")",bef," <-> ",aft
351        ENDIF
352      ENDDO
353    ENDIF
354
355    ! compute tendencies:
356    ! all of these tendencies are in X.m-3 !
357    dm0aer = (zm0a - sm0a)*mm_rhoair
358    dm3aer = (zm3a - sm3a)*mm_rhoair
359    dm0ccn = (zm0n - sm0n)*mm_rhoair
360    dm3ccn = (zm3n - sm3n)*mm_rhoair
361
362    dm3iX  = (zm3iX - sm3iX)                    ! this one in X.kg-1 (temporary) !
363    dvapX  = -xESP%rho * dm3iX / xESP%fmol2fmas ! in order to compute this one in mol.mol-1
364    dm3iX  = dm3iX*mm_rhoair                    ! update ice tendencies in X.m-3 !
365
366  END SUBROUTINE nc_esp
367
368  SUBROUTINE nuc_rate(rccn,temp,xESP,pvp,sat,rate)
369    !! Get nucleation rate.
370    !!
371    !! The method computes the heterogeneous nucleation rate for the given specie on a fractal particle
372    !! of size __rccn__.
373    !! Except __xESP__, all arguments are vectors of the same size (vertical grid).
374    REAL(kind=mm_wp), INTENT(in), DIMENSION(:)  :: rccn !! Radius of the cloud condensation nuclei (m).
375    REAL(kind=mm_wp), INTENT(in), DIMENSION(:)  :: temp !! Temperature (K).
376    REAL(kind=mm_wp), INTENT(in), DIMENSION(:)  :: pvp  !! Partial vapor pressure of X specie (Pa).
377    REAL(kind=mm_wp), INTENT(in), DIMENSION(:)  :: sat  !! Saturation ratio of given specie (--).
378    TYPE(mm_esp), INTENT(in)                    :: xESP !! X specie properties (--).
379    REAL(kind=mm_wp), INTENT(out), DIMENSION(:) :: rate !! The nucleation rate (\(m^-{2}.s^{-1}\)).
380    INTEGER          :: i
381    REAL(kind=mm_wp) :: r,t,s,sig,nX,rstar,gstar,x,zeldov,deltaf,fsh,fstar
382
383    rate(:) = 0._mm_wp
384    ! Activation condition
385    DO i=1, SIZE(rccn)
386      s = sat(i)
387      IF (s > 1._mm_wp) THEN
388        t = temp(i) ; r = rccn(i)
389        sig = mm_sigX(t,xESP)
390        nX    = pvp(i)/mm_kboltz/t
391        rstar = 2._mm_wp*sig*xESP%vol/(mm_kboltz*t*dlog(s))
392        ! curvature radius
393        x     = r/rstar
394        fsh   = mm_fshape(xESP%mteta,x)
395        fstar = (4._mm_wp*mm_pi/3._mm_wp)*sig*(rstar**2.)*fsh
396        deltaf=MIN(MAX((2.*mm_fdes-mm_fdif-fstar)/(mm_kboltz*t),-100._mm_wp),100._mm_wp)
397        IF (deltaf > -100._mm_wp) THEN
398          gstar  = 4._mm_wp*mm_pi*(rstar**3)/(3._mm_wp*xESP%vol)
399          zeldov = dsqrt(fstar/(3._mm_wp*mm_pi*mm_kboltz*t*(gstar**2)))
400          rate(i)= zeldov*mm_kboltz*t*(nX*rstar)**2._mm_wp*dexp(deltaf)/(fsh*mm_nus*xESP%mas)
401        ENDIF
402      ENDIF
403    ENDDO
404
405    RETURN
406  END SUBROUTINE nuc_rate
407
408  SUBROUTINE growth_rate(temp,pres,pXsat,xESP,seq,drad,rate)
409    !! Get growth rate through condensation/evaporation process.
410    !!
411    !! The method computes the growth rate a drop through condensation/evaporation processes:
412    !!
413    !! $$ r \times \frac{dr}{dt} = g_{rate} \times (S - S_{eq}) $$
414    !!
415    !! Except __xESP__ which is a scalar, all arguments are vectors of the same size (vertical grid).
416    REAL(kind=mm_wp), INTENT(in), DIMENSION(:)  :: temp  !! Temperature (K).
417    REAL(kind=mm_wp), INTENT(in), DIMENSION(:)  :: pres  !! Pressure level (Pa).
418    REAL(kind=mm_wp), INTENT(in), DIMENSION(:)  :: pXsat !! Saturation vapor pressure of specie (Pa).
419    TYPE(mm_esp), INTENT(in)                    :: xESP  !! Specie properties.
420    REAL(kind=mm_wp), INTENT(in), DIMENSION(:)  :: seq   !! Equilibrium saturation near the drop.
421    REAL(kind=mm_wp), INTENT(in), DIMENSION(:)  :: drad  !! Drop radius (m).
422    REAL(kind=mm_wp), INTENT(out), DIMENSION(:) :: rate  !! Growth rate (\(m^{2}.s^{-1}\)).
423    REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE :: k,knu,slf,rkc,rdc,l,dv
424    INTEGER                                     :: n
425    n = SIZE(temp)
426    ALLOCATE(k(n),knu(n),slf(n),rkc(n),rdc(n),l(n),dv(n))
427
428    ! N2 (air) Thermal conductivity (where does it come from ?)
429    k(:) = ( 2.857e-2_mm_wp*temp-0.5428_mm_wp)*4.184e-3_mm_wp
430    ! Gas mean free path
431    l(:) = mm_lambda_g(temp,pres)
432    ! Diffusion coefficient of X gas
433    Dv(:) = 1._mm_wp/3._mm_wp*dsqrt(8._mm_wp*mm_rgas*temp(:)/(mm_pi*xESP%masmol))*mm_kboltz*temp(:) / &
434      (mm_pi*pres(:)*(mm_air_rad+xESP%ray)**2* dsqrt(1._mm_wp+xESP%fmol2fmas))
435    knu(:) = l(:)/drad(:)                                         ! The knudsen number of the drop
436    slf(:) = (1.333_mm_wp+0.71_mm_wp/knu(:))/(1._mm_wp+1._mm_wp/knu(:)) ! Slip flow correction
437    Dv(:)  = Dv(:)/(1._mm_wp+slf(:)*knu(:))
438    ! latent heat resistance coefficient
439    rkc(:) = mm_lheatX(temp(:),xESP)**2 * xESP%rho * xESP%masmol / (k(:)*mm_rgas*temp(:)**2)
440    ! Diffusion resistance coefficient
441    rdc(:) = mm_rgas * temp(:) * xESP%rho / (Dv(:)*pXsat(:)*xESP%masmol)
442    ! Growth rate: rdr/dt = rate * (S-Seq) ; rate is returned
443    rate(:) = 1._mm_wp / (seq(:)*rkc(:)+rdc(:))
444    RETURN
445  END SUBROUTINE growth_rate
446
447
448  !-----------------------------------------------------------------------------
449  ! SEDIMENTATION PROCESS RELATED METHODS
450  !-----------------------------------------------------------------------------
451
452  SUBROUTINE mm_cloud_sedimentation(dm0n,dm3n,dm3i)
453    !! Compute the tendency of _clouds_ related moments through sedimentation process.
454    !!
455    !! The method computes the tendencies of moments related to cloud microphysics through
456    !! sedimentation process. The algorithm used here differs from
457    !! [[mm_haze(module):mm_haze_sedimentation(subroutine)]] as all moments settle with the same
458    !! terminal velocity which is computed with the average drop radius of the size distribution.
459    !! We simply compute an _exchange matrix_ that stores the new positions of each cells through
460    !! sedimentation process and then computes the matrix
461    !! product with input moments values to get final tendencies.
462    REAL(kind=mm_wp), INTENT(out), DIMENSION(:)   :: dm0n
463    !! Tendency of the 0th order moment of the ccn distribution (\(m^{-3}\)).
464    REAL(kind=mm_wp), INTENT(out), DIMENSION(:)   :: dm3n
465    !! Tendency of the 3rd order moment of the ccn distribution (\(m^{3}.m^{-3}\)).
466    REAL(kind=mm_wp), INTENT(out), DIMENSION(:,:) :: dm3i
467    !! Tendencies of the 3rd order moment of each ice component of the cloud (\(m^{3}m^{-3}\)).
468    INTEGER                                   :: im,nm
469    REAL(kind=mm_wp), DIMENSION(:,:), ALLOCATABLE :: moms, momsf,chg_matrix
470
471    nm = 2 + mm_nesp
472    ALLOCATE(moms(mm_nla,nm),momsf(mm_nla,nm),chg_matrix(mm_nla,mm_nla))
473    ! Initializes moms
474    moms(:,1)  = mm_m0ccn * mm_dzlev
475    moms(:,2)  = mm_m3ccn * mm_dzlev
476    DO im=1,mm_nesp ; moms(:,2+im) = mm_m3ice(:,im) * mm_dzlev ; ENDDO
477    ! Computes exchange matrix
478    CALL exchange(mm_drad,mm_drho,mm_dt,chg_matrix)
479    ! Computes final moments values
480    DO im=1,nm ; momsf(:,im) = MATMUL(chg_matrix,moms(:,im)) ; ENDDO
481    ! Computes tendencies (converted in X.m-3)
482    dm0n = (momsf(:,1)-moms(:,1))/mm_dzlev
483    dm3n = (momsf(:,2)-moms(:,2))/mm_dzlev
484    DO im=1,mm_nesp ; dm3i(:,im) = (momsf(:,2+im)-moms(:,2+im))/mm_dzlev ; ENDDO
485    RETURN
486  END SUBROUTINE mm_cloud_sedimentation
487
488  SUBROUTINE exchange(rad,rhog,dt,matrix)
489    !! Compute the exchange matrix.
490    !!
491    !! The subroutine computes the matrix exchange used by
492    !! [[mm_clouds(module):mm_cloud_sedimentation(subroutine)]] to compute moments tendencies
493    !! through sedimentation process. Both __rad__ and __rhog__ must be vector with relevant
494    !! values over the atmospheric vertical structure. __matrix__ is square 2D-array with same
495    !! dimension size than __rad__.
496    REAL(kind=mm_wp), INTENT(in), DIMENSION(:) :: rad
497    !! Cloud drop radius over the atmospheric vertical structure (m).
498    REAL(kind=mm_wp), INTENT(in), DIMENSION(:) :: rhog
499    !! Cloud drop density over the atmospheric vertical structure (\(kg.m^{-3}\)).
500    REAL(kind=mm_wp), INTENT(in)               :: dt
501    !! Timestep (s).
502    REAL(kind=mm_wp), INTENT(out)              :: matrix(:,:)
503    !! The output _exchange matrix_.
504    INTEGER                                     :: nz,i,j,jj,jinf,jsup
505    REAL(kind=mm_wp)                            :: zni,znip1,xf,xft,xcnt
506    REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE :: puit
507    REAL(kind=mm_wp)                            :: cpte,cpte2
508    REAL(kind=mm_wp)                            :: zsurf
509    INTEGER, PARAMETER                          :: ichx  = 1
510    matrix = 0._mm_wp ; nz = SIZE(rad) ; ALLOCATE(puit(nz))
511    zsurf = mm_zlev(nz)
512   
513    ! compute exchange matrix
514    DO i=1,nz
515      puit(i) = 0._mm_wp
516      xcnt = 0._mm_wp
517      ! computes drop move (i.e. its new positions)
518      CALL getnzs(ichx,i,rad,rhog,dt,zni,znip1)
519
520      ! Peculiar case : Ground level precipitation [znip1 < zsurf && (zni < zsurf || zni > zsurf)]
521      ! - complete precipitation [ znip1 <= 0 && zni <= 0 ] :
522      IF(zni <= zsurf .and. znip1 <= zsurf) THEN
523        xft=0._mm_wp
524        xf=1._mm_wp
525        xcnt=xcnt+xf
526        puit(i)=puit(i)+xf
527      ENDIF
528      ! - partial precipitation [ znip1 <= zsurf && zni > zsurf ] :
529      IF (zni > zsurf .and. znip1 <= zsurf) THEN
530        xft=(zni-zsurf)/(zni-znip1)
531        xf=(1.-xft)
532        xcnt=xcnt+xf
533        puit(i)=puit(i)+xf
534      ENDIF
535      ! General case : no ground precipitation [ znip1 > zsurf && zni > zsurf ]
536      IF (zni > zsurf .and. znip1 > zsurf) THEN
537        xft = 1._mm_wp ! on a la totalite de la case
538        xf  = 0._mm_wp
539        xcnt=xcnt+xf
540        puit(i)=puit(i)+xf
541      ENDIF
542      IF (zni < znip1) THEN
543        WRITE(*,'("[EXCHANGES] WARNING, missing this case :",2(2X,ES10.3))') zni, znip1
544      ENDIF
545     
546      ! Fix minimum level to the ground
547      znip1 = MAX(znip1,zsurf)
548      zni   = MAX(zni,zsurf)
549      ! Locate new "drop" position in the verical grid
550      jsup=nz+1
551      jinf=nz+1
552      DO j=1,nz
553        IF (zni<=mm_zlev(j).and.zni>=mm_zlev(j+1)) jsup=j
554        IF (znip1<=mm_zlev(j).and.znip1>=mm_zlev(j+1)) jinf=j
555      ENDDO
556      ! Volume is out of range: (all drops have touched the ground!)
557      ! Note: cannot happen here, it has been treated previously :)
558      IF (jsup>=nz+1.and.jinf==jsup) THEN
559        WRITE(*,'(a)') "[EXCHANGE] speaking: The impossible happened !"
560        call EXIT(666)
561      ENDIF
562      ! Volume inside a single level
563      IF (jsup==jinf.and.jsup<=nz) THEN
564        xf=1._mm_wp
565        xcnt=xcnt+xft*xf
566        matrix(jinf,i)=matrix(jinf,i)+xft*xf
567      ENDIF
568
569      ! Volume over 2 levels
570      IF (jinf==jsup+1) THEN
571        xf=(zni-mm_zlev(jinf))/(zni-znip1)
572        xcnt=xcnt+xf*xft
573        IF(jsup<=nz) THEN
574          matrix(jsup,i)=matrix(jsup,i)+xft*xf
575        ENDIF
576        xf=(mm_zlev(jinf)-znip1)/(zni-znip1)
577        xcnt=xcnt+xf*xft
578        IF (jinf<=nz) THEN
579          matrix(jinf,i)=matrix(jinf,i)+xft*xf
580        ENDIF
581      ENDIF
582
583      ! Volume over 3 or more levels
584      IF (jinf > jsup+1) THEN
585        ! first cell
586        xf=(zni-mm_zlev(jsup+1))/(zni-znip1)
587        xcnt=xcnt+xf*xft
588        matrix(jsup,i)=matrix(jsup,i)+xft*xf
589        ! last cell
590        xf=(mm_zlev(jinf)-znip1)/(zni-znip1)
591        xcnt=xcnt+xf*xft
592        matrix(jinf,i)=matrix(jinf,i)+xft*xf
593        ! other :)
594        DO jj=jsup+1,jinf-1
595          xf=(mm_zlev(jj)-mm_zlev(jj+1))/(zni-znip1)
596          xcnt=xcnt+xf*xft
597          matrix(jj,i)=matrix(jj,i)+xft*xf
598        ENDDO
599      ENDIF
600    ENDDO
601
602    ! checking if everything is alright if debug enabled...
603    IF (mm_debug) THEN
604      cpte=0._mm_wp ; cpte2=0._mm_wp
605      DO j=1,nz
606        DO jj=1,nz
607          cpte=cpte+matrix(jj,j)
608        ENDDO
609      ENDDO
610      cpte2=cpte+sum(puit)
611      IF (abs(cpte2-nz)>1.e-4_mm_wp) THEN
612        WRITE(*,'("[EXCHANGE] speaking: tx expl (/nz):",2(2X,ES10.3))') cpte,cpte2
613      ENDIF
614    ENDIF
615
616    RETURN
617  END SUBROUTINE exchange
618
619  SUBROUTINE getnzs(ichx,idx,rad,rho,dt,zni,zns)
620    !! Compute displacement of a cell under sedimentation process.
621    !!
622    !! The method computes the new position of a _drop cell_ through sedimentation process as
623    !! descibed in the following scheme:
624    !!
625    !! ![Cloud sedimentation scheme](|media|/cloud_sed_scheme.svg)
626    !!
627    !! New positions are returned in __zni__ and __zns__ ouptut arguments.
628    !!
629    !! @note
630    !! The method uses directly [[mm_globals(module):mm_play(variable)]], [[mm_globals(module):mm_plev(variable)]],
631    !! [[mm_globals(module):mm_temp(variable)]],[[mm_globals(module):mm_btemp(variable)]],
632    !! [[mm_globals(module):mm_zlay(variable)]] and [[mm_globals(module):mm_zlev(variable)]] and uses __idx__ to
633    !! get the relevant value to use on the vertical grid.
634    INTEGER, INTENT(in)                        :: ichx
635    !! Velocity extrapolation control flag (0 for linear, 1 for exponential -preferred -).
636    INTEGER, INTENT(in)                        :: idx
637    !! Initial position of the drop (subscript of vertical layers vectors).
638    REAL(kind=mm_wp), INTENT(in), DIMENSION(:) :: rad
639    !! Cloud drop radius over the atmospheric vertical structure (m).
640    REAL(kind=mm_wp), INTENT(in), DIMENSION(:) :: rho
641    !! Cloud drop density over the atmospheric vertical structure (\(kg.m^{-3}\)).
642    REAL(kind=mm_wp), INTENT(in)               :: dt
643    !! Timestep (s).
644    REAL(kind=mm_wp), INTENT(out)              :: zni
645    !! Final layer upper boundary position (m).
646    REAL(kind=mm_wp), INTENT(out)              :: zns
647    !! Final layer lower boundary position (m).
648    REAL(kind=mm_wp)            :: ws,wi,w,zi,zs
649    REAL(kind=mm_wp)            :: alpha,argexp,v0,arg1,arg2
650    INTEGER                 :: i,nz
651    REAL(kind=mm_wp), PARAMETER :: es = 30._mm_wp
652
653    nz = SIZE(rad)
654    ! Linear extrapolation of velocity
655    IF (ichx==0) THEN
656      ! velocity upper interface
657      ws = wsettle(mm_plev(idx),mm_btemp(idx),mm_zlev(idx),rho(idx),rad(idx))
658      IF (idx==nz) THEN
659        ! veloctity center layer
660        wi = wsettle(mm_play(idx),mm_temp(idx),mm_zlay(idx),rho(idx),rad(idx))
661      ELSEIF(idx<nz) THEN
662        ! velocity lower interface
663        wi = wsettle(mm_plev(idx+1),mm_btemp(idx+1),mm_zlev(idx+1), &
664          rho(idx+1),rad(idx+1))
665      ELSE
666        WRITE(*,'(a)') "[getnzs] speaking:"
667        WRITE(*,'(a)') "This is the fatal error..."
668        WRITE(*,'(a)') "index is higher than number of levels"
669        call EXIT(111)
670      ENDIF
671      w   = (ws+wi)/2._mm_wp
672      zni = mm_zlev(idx)-w*dt
673      zns = mm_zlev(idx)-mm_dzlev(idx)-w*dt
674      RETURN
675      ! Exponential extrapolation of velocity
676
677    ELSEIF(ichx==1) THEN
678      zs = mm_zlev(idx)
679      ws = wsettle(mm_plev(idx),mm_btemp(idx),zs,rho(idx),rad(idx))
680      zi=mm_zlay(idx)
681      wi = wsettle(mm_play(idx),mm_temp(idx),zi,rho(idx),rad(idx))
682      ! ws & wi must be different !
683      IF(dabs(wi-ws)/wi <= 1.e-3_mm_wp)  wi=ws/1.001_mm_wp
684      IF (wi /= 0._mm_wp) alpha = dlog(ws/wi)/(zs-zi)  ! alpha < 0 if wi > ws
685      !   -es < argexp < es
686      argexp=MAX(MIN(alpha*zs,es),-es)
687      v0 = ws/dexp(argexp)
688      arg1=1._mm_wp+v0*alpha*dexp(argexp)*dt
689      argexp=MAX(MIN(alpha*(mm_zlev(idx)-mm_dzlev(idx)),es),-es)
690      arg2=1._mm_wp+v0*alpha*dexp(argexp)*dt
691      IF (arg1<=0._mm_wp.OR.arg2<=0._mm_wp) THEN
692        ! correct velocity
693        ! divides the velocity argument in arg1 and arg2 :
694        ! argX=1+alpha*v0*exp(alpha*z)*dt <==> argX-1=alpha*v0*exp(alpha*z)*dt
695        ! ==> argX' = 1 + (argX-1)/2  <==> argX' = (1+argX)/2.
696        DO i=1,25
697          IF (arg1<=0._mm_wp.OR.arg2<=0._mm_wp) THEN
698            IF (mm_debug) &
699              WRITE(*,'((a),I2.2,(a))') "[getnzs] must adjust velocity (iter:",i,"/25)"
700            arg1=(arg1+1._mm_wp)/2._mm_wp ; arg2=(arg2+1._mm_wp)/2._mm_wp
701          ELSE
702            EXIT
703          ENDIF
704        ENDDO
705        ! sh** we have to stop
706        IF (i>25) THEN
707          WRITE(*,'(a)')"[getnzs] speaking:"
708          WRITE(*,'(a)') "Cannot adjust velocities"
709          call EXIT(111)
710        ENDIF
711      ENDIF
712
713      zni = mm_zlev(idx)-dlog(arg1)/alpha
714      zns = mm_zlev(idx)-mm_dzlev(idx)-dlog(arg2)/alpha
715
716      RETURN
717    ENDIF
718  END SUBROUTINE getnzs
719
720  ELEMENTAL FUNCTION wsettle(p,t,z,rho,rad) RESULT(w)
721    !! Compute the settling velocity of a spherical particle.
722    !!
723    !! The method computes the effective settling velocity of spherical particle of
724    !! radius __rad__. It accounts for the slip-flow transition (no approximation).
725    REAL(kind=mm_wp), INTENT(in) :: p   !! The pressure level (Pa).
726    REAL(kind=mm_wp), INTENT(in) :: t   !! The temperature (K).
727    REAL(kind=mm_wp), INTENT(in) :: z   !! The altitude level (m).
728    REAL(kind=mm_wp), INTENT(in) :: rho !! Density of the particle (\(kg.m^{-3}\)).
729    REAL(kind=mm_wp), INTENT(in) :: rad !! Radius of the particle (m).
730    REAL(kind=mm_wp) :: w               !! Settling velocity (\(m.s^{-1}\)).
731    REAL(kind=mm_wp) :: Us, Fc, kn, wtmp, wmax
732    REAL(kind=mm_wp), PARAMETER :: ra = 1.75e-10_mm_wp
733   
734    ! Knudsen number
735    kn = (mm_kboltz * t) / (p * 4._mm_wp * sqrt(2._mm_wp) * mm_pi * ra**2) / rad
736
737    ! Computes slip-flow correction : Fc = 1 + (mm_akn * mm_lambda_g(t,p) / rad)
738    Fc = (1._mm_wp + 1.2517_mm_wp*kn + 0.4_mm_wp*kn*dexp(-1.1_mm_wp/kn))
739   
740    ! Computes Stokes settling velocity
741    Us = (2._mm_wp * rad**2 * rho * mm_effg(z)) / (9._mm_wp * mm_eta_g(t))
742   
743    ! Computes settling velocity (correction factor : x3.0)
744    w = Us * Fc * 3._mm_wp
745
746    ! Imposes a velocity limit
747    wmax = 20._mm_wp ! 20 m/s [Lorenz 1993]
748    wtmp = (1._mm_wp / w) + (1._mm_wp / wmax)
749    w = 1._mm_wp / wtmp
750  END FUNCTION wsettle
751
752  FUNCTION get_mass_flux(rho,m3) RESULT(flx)
753    !> Get the mass flux of (clouds related) moment through sedimention.
754    !!
755    !! @warning
756    !! The method is __only__ valid for cloud moments (i.e. ice or ccn). It calls
757    !! [[mm_clouds(module):wsettle(function)]] that compute the _mean_ settling velocity of a
758    !! cloud drop.
759    !!
760    !! @note
761    !! The computed flux is always positive.
762    REAL(kind=mm_wp), INTENT(in)               :: rho
763    !! Tracer density (\(kg.m^{-3}\)).
764    REAL(kind=mm_wp), INTENT(in), DIMENSION(:) :: m3
765    !! Vertical profile of the total volume of tracer (i.e. M3) from __TOP__ to __GROUND__ (\(m^{3}.m^{-3}\)).
766    REAL(kind=mm_wp), DIMENSION(SIZE(m3))      :: flx
767    !! Mass sedimentation fluxes at each layer from __TOP__ to __GROUND__ (\(kg.m^{-2}.s^{-1}\)).
768    REAL(kind=mm_wp), SAVE :: pifac = (4._mm_wp * mm_pi) / 3._mm_wp
769    flx = rho * pifac * m3 * wsettle(mm_play,mm_temp,mm_zlay,mm_drho,mm_drad)
770    RETURN
771  END FUNCTION get_mass_flux
772
773END MODULE MM_CLOUDS
Note: See TracBrowser for help on using the repository browser.