source: trunk/LMDZ.VENUS/libf/phyvenus/cloudvenus/sed_and_prod_mad.F90 @ 3210

Last change on this file since 3210 was 2259, checked in by emillour, 6 years ago

Venus GCM:
Some cleanup to be able to compile with modern gfortran compiler (5.5+) which
complained about the need for an explicit interface in sed_and_prod_mad.F90
for get_weff().
Turned it into a module in the process.
EM
File size: 20.7 KB
Line 
1MODULE SED_AND_PROD_MAD
2
3IMPLICIT NONE
4
5!
6! All parameter variable should be modified regarding the sudied system.
7!
8!
9! All vectors arguments of the subroutines defined here should be defined from TOP of
10! the atmosphere to the Ground.
11!
12! There should be one more vertical level than layers. Layer thickness vector should be
13! defined on the number of layer. The last item dzlay(nla), where nla is the number of layer,
14! is arbitrary and can be set to dzlay(nla-1). Usually dzlay is computed from zlay:
15!
16!    dzlay(1:nla-1) = zlay(1:nla-1)-zlay(2:nla)
17!    dzlay(nla) = dzlay(nla-1)
18!
19! Level thickness (dzlev) is also defined on the number of layer and can be computed as follow:
20!
21!    dzlev(1:nla) = zlev(1:nla)-zlev(2:nle)
22!
23! Where nle is number of vertical level (i.e. nla+1).
24!
25!
26! eta_g and lambda_g functions should be modified regarding the studied atmosphere. The current parameters
27! stands for an N2 atmosphere.
28!
29! eta_g    --> air viscosity at given temperature.
30! lambda_g --> air mean free path at given temperature and pressure.
31!
32! By J.Burgalat
33!
34! ft = Theoretical Settling velocity at each vertical levels with Fiadero correction
35
36CONTAINS
37
38!============================================================================
39! MOMENT RELATED METHODS
40!============================================================================
41FUNCTION get_r0_1(M0,M3) RESULT(ret)
42
43  use free_param, only: sigma1
44
45  IMPLICIT NONE
46
47  !! Compute the mean radius of the size distribution of the first mode.
48  !!
49  !! $$ r_{0,1} = \left(\dfrac{M_{3}}{M_{0}}/ \alpha_{1}(3) \right)^{1/3} $$
50  REAL, INTENT(in) :: M0
51  !! 0th order moment of the size-distribution
52  REAL, INTENT(in) :: M3
53  !! 3rd order moment of the size-distribution
54  REAL             :: ret, alpha_k
55  !! Mean radius of the mode 1
56  ret = (m3/m0 / alpha_k(3.D0,sigma1))**0.333333D0 !SG: ok+
57END FUNCTION get_r0_1
58
59
60!*****************************************************************************
61FUNCTION get_r0_2(M0,M3) RESULT(ret)
62
63  use free_param, only: sigma2
64
65  IMPLICIT NONE
66
67  !! Compute the mean radius of the size distribution of the first mode.
68  !!
69  !! $$ r_{0,2} = \left(\dfrac{M_{3}}{M_{0}}/ \alpha_{2}(3) \right)^{1/3} $$
70  REAL, INTENT(in) :: M0
71  !! 0th order moment of the size-distribution
72  REAL, INTENT(in) :: M3
73  !! 3rd order moment of the size-distribution
74  REAL             :: ret, alpha_k
75  !! Mean radius of the mode 2
76  ret = (m3/m0 / alpha_k(3.D0,sigma2))**0.333333D0
77END FUNCTION get_r0_2 !SG: ok ++
78
79
80!============================================================================
81! ATMOSPHERE RELATED METHODS
82!============================================================================
83FUNCTION effg(z) RESULT(ret)
84  !! Compute effective gravitational acceleration.
85
86  use donnees, only: g0, rpla
87
88  IMPLICIT NONE
89
90  REAL, INTENT(in) :: z !! Altitude in meters
91  REAL :: ret           !! Effective gravitational acceleration in \(m.s^{-2}\)
92  ret = g0 * (rpla/(rpla+z))**2
93  RETURN
94END FUNCTION effg
95
96
97!============================================================================
98!                 PRODUCTION PROCESS RELATED METHOD
99!============================================================================
100SUBROUTINE aer_production(dt,nlay,plev,zlev,zlay,dm0_1,dm3_1)
101
102  use free_param, only: sigz, p_aer, r_aer, rho_aer, sig_aer, tx_prod
103  use donnees, only: pi
104
105  IMPLICIT NONE
106
107  !! Compute the production of aerosols moments.
108  !!
109  !! The method computes the tendencies of M0 and M3 of the CCN precursor.
110  !! Production values are distributed along a normal law in altitude, centered  at
111  !! p_prod pressure level with a fixed sigma of 20km.
112  !!
113  !! First M3 tendency is computed and M0 is retrieved using the inter-moments relation
114  !! of the first mode of drop size-distribution with a fixed mean radius (r_aer).
115  REAL, INTENT(in) :: dt
116  INTEGER, INTENT(in) :: nlay
117  REAL, INTENT(in),  DIMENSION(nlay+1) :: plev    !! Pressure at each vertical level (Pa).
118  REAL, INTENT(in),  DIMENSION(nlay+1) :: zlev    !! Altitude of each vertical level (m).
119  REAL, INTENT(in),  DIMENSION(nlay)   :: zlay    !! Altitude at the center of each vertical layer (m).
120  REAL, INTENT(out), DIMENSION(nlay)   :: dm0_1   !! Tendency of M0 (\(m^{-3}\)).
121  REAL, INTENT(out), DIMENSION(nlay)   :: dm3_1   !! Tendency of M3 (\(m^{3}.m^{-3}\)).
122  INTEGER               :: i,nla,nle
123  REAL, DIMENSION(nlay) :: dzlev
124  REAL                  :: zprod,cprod, alpha_k
125
126  REAL, PARAMETER :: fnorm = 1.D0/(dsqrt(2.D0*pi)*sigz)
127  REAL, PARAMETER :: znorm = dsqrt(2.D0)*sigz
128
129  nle = SIZE(zlev)
130  nla = nlay+1
131  dzlev = zlev(1:nle-1)-zlev(2:nle) !distance between two levels, in meter
132  zprod = -1.0D0 ! initialization
133 
134  ! locate production altitude
135  DO i=1, nla-1
136     IF (plev(i) .lt. p_aer.AND.plev(i+1) .ge. p_aer) THEN
137        zprod = zlay(i) ; EXIT ! the altitude of production is determined here
138     ENDIF
139  ENDDO
140  IF (zprod < 0.D0) THEN ! test
141     WRITE(*,'(a)') "In aer_prod, cannot find production altitude"
142     call EXIT(11)
143  ENDIF
144
145  dm3_1(:)= tx_prod *0.75D0/pi *dt / rho_aer / 2.D0 / dzlev(1:nla) * &
146       (erf((zlev(1:nla)-zprod)/znorm) - &
147       erf((zlev(2:nla+1)-zprod)/znorm))
148
149  dm0_1(:) = dm3_1(:)/(r_aer**3*alpha_k(3.D0,sig_aer))
150
151END SUBROUTINE aer_production
152
153
154!============================================================================
155! SEDIMENTATION PROCESS RELATED METHODS
156!============================================================================
157SUBROUTINE let_me_fall_in_peace(dt,nlay,dzlay,mk,ft,dmk)
158  !! Compute the tendency of the moment through sedimentation process.
159  !!
160  !!
161  !! The method computes the time evolution of the \(k^{th}\) order moment through sedimentation:
162  !!
163  !! $$ \dfrac{dM_{k}}{dt} = \dfrac{\Phi_{k}}{dz} $$
164  !!
165  !! The equation is resolved using a [Crank-Nicolson algorithm](http://en.wikipedia.org/wiki/Crank-Nicolson_method).
166  !!
167  !! Sedimentation algorithm is quite messy. It appeals to the dark side of the Force and uses evil black magic spells
168  !! from ancient times. It is based on \cite{toon1988b,fiadeiro1977,turco1979} and is an update of the algorithm
169  !! originally implemented in the LMDZ-Titan 2D GCM.
170
171  IMPLICIT NONE
172
173  REAL, INTENT(in)  :: dt
174  INTEGER, INTENT(in) :: nlay
175  REAL, INTENT(in), DIMENSION(nlay)  :: dzlay !! Atmospheric thickness between the center of 2 adjacent layers (\(m\))
176  REAL, INTENT(in), DIMENSION(nlay)  :: mk    !! \(k^{th}\) order moment to sediment (in \(m^{k}\)).
177  REAL, INTENT(in), DIMENSION(nlay+1):: ft    !! Downward sedimentation flux  (effective velocity of the moment).
178  REAL, INTENT(out), DIMENSION(nlay) :: dmk   !! Tendency of \(k^{th}\) order moment (in \(m^{k}.m^{-3}\)).
179  INTEGER                            :: i,nla,nle
180  REAL, DIMENSION(nlay) :: as,bs,cs,mko
181
182  nla = SIZE(mk,DIM=1) ; nle = nla+1
183
184  mko(1:nla) = 0.D0 !initialization
185  cs(1:nla) = ft(2:nla+1) - dzlay(1:nla)/dt 
186  IF (ANY(cs > 0.D0)) THEN
187     ! implicit case
188     as(1:nla) = ft(1:nla)
189     bs(1:nla) = -(ft(2:nle)+dzlay(1:nla)/dt)
190     cs(1:nla) = -dzlay(1:nla)/dt*mk(1:nla)
191     ! (Tri)diagonal matrix inversion
192     mko(1) = cs(1)/bs(1) ! at the top
193     DO i=2,nla
194        mko(i) = (cs(i)-mko(i-1)*as(i))/bs(i)
195     ENDDO
196  ELSE
197     ! explicit case
198     as(1:nla)=-dzlay(1:nla)/dt
199     bs(1:nla)=-ft(1:nla)
200     ! boundaries
201     mko(1)=cs(1)*mk(1)/as(1)
202     mko(nla)=(bs(nla)*mk(nla-1)+cs(nla)*mk(nla))/as(nla)
203     ! interior
204     mko(2:nla-1)=(bs(2:nla-1)*mk(1:nla-2) + &
205          cs(2:nla-1)*mk(2:nla-1)   &
206          )/as(2:nla-1)
207  ENDIF
208  dmk = mko - mk
209
210  RETURN
211END SUBROUTINE let_me_fall_in_peace
212
213
214!*****************************************************************************
215SUBROUTINE get_weff(dt,nlay,plev,zlev,dzlay,btemp,mk,k,rc,sig,wth,corf)
216!     call get_weff(plev,zlev,dzlay,btemp,m0ccn,0.D0,r0,sig,wth,fdcor)
217  !! Get the effective settling velocity for aerosols moments.
218  !!
219  !! This method computes the effective settling velocity of the \(k^{th}\) order moment of aerosol
220  !! tracers. The basic settling velocity (\(v^{eff}_{M_{k}}\)) is computed using the following
221  !! equation:
222  !!
223  !! $$
224  !! \begin{eqnarray*}
225  !! \Phi^{sed}_{M_{k}} &=& \int_{0}^{\infty} n(r) r^{k} \times w(r) dr
226  !!                     == M_{k}  \times v^{eff}_{M_{k}} \\
227  !!     v^{eff}_{M_{k} &=& \dfrac{2 \rho g r_{c}^{\dfrac{3D_{f}-3}{D_{f}}}}
228  !!     {r_{m}^{D_{f}-3}/D_{f}} \times \alpha(k)} \times \left( \alpha \right(
229  !!     \frac{D_{f}(k+3)-3}{D_{f}}\left) + \dfrac{A_{kn}\lambda_{g}}{r_{c}^{
230  !!     3/D_{f}}} \alpha \right( \frac{D_{f}(k+3)-6}{D_{f}}\left)\left)
231  !! \end{eqnarray*}
232  !! $$
233  !!
234  !! \(v^{eff}_{M_{k}\) is then corrected to reduce numerical diffusion of the sedimentation algorithm
235  !! as defined in \cite{toon1988b}.
236  !!
237  !! @warning
238  !! Both __df__, __rc__ and __afun__ must be consistent with each other otherwise wrong values will
239  !! be computed.
240
241  use free_param, only: rho_aer
242  use donnees, only: akn, df, fiadero_max, fiadero_min, no_fiadero_w, rmono
243 
244  IMPLICIT NONE
245  REAL, INTENT(in)  :: dt
246  INTEGER, INTENT(in) :: nlay
247  REAL, INTENT(in), DIMENSION(nlay+1)          :: plev
248  !! Pressure at each level (Pa).
249  REAL, INTENT(in), DIMENSION(nlay+1)          :: zlev
250  !! Altitude at each level (m).
251  REAL, INTENT(in), DIMENSION(nlay)            :: dzlay
252  !! Altitude at the center of each layer (Pa).
253  REAL, INTENT(in), DIMENSION(nlay)            :: btemp
254  !! Temperature at each level (T).
255  REAL, INTENT(in), DIMENSION(nlay)            :: mk
256  !! Moment of order __k__ (\(m^{k}.m^{-3}\)).
257  REAL, INTENT(in), DIMENSION(nlay)            :: rc
258  !! Characteristic radius associated to the moment at each vertical __levels__.
259  REAL, INTENT(in)                          :: k, sig
260  !! Fractal dimension of the aersols.
261  REAL, INTENT(out), DIMENSION(nlay+1)           :: wth
262  !! Theoretical Settling velocity at each vertical __levels__ (\( wth \times corf = weff\)).
263  REAL, INTENT(out), DIMENSION(nlay+1), OPTIONAL :: corf
264  !! _Fiadero_ correction factor applied to the theoretical settling velocity at each vertical __levels__.
265  INTEGER                             :: i,nla,nle
266  REAL                       :: af1,af2,ar1,ar2
267  REAL                       :: csto,cslf,ratio,wdt,dzb
268  REAL                       :: rb2ra, VISAIR, FPLAIR, alpha_k
269  REAL, DIMENSION(nlay+1) :: zcorf
270  ! ------------------
271  write(*,*)'yupitururu0'
272  nla = SIZE(mk,DIM=1)
273  nle = nla+1
274
275  write(*,*)'yupitururu1'
276  wth(:) = 0.D0
277  zcorf(:) = 1.D0
278 
279  ar1 = (3.D0*df -3.D0)/df   
280  ar2 = (3.D0*df -6.D0)/df
281  af1 = (df*(k+3.D0)-3.D0)/df
282  af2 = (df*(k+3.D0)-6.D0)/df
283  rb2ra = rmono**((df-3.D0)/df)
284  write(*,*)'yupitururu2'
285  DO i=2,nla
286     IF (rc(i-1) <= 0.D0) CYCLE
287     dzb = (dzlay(i)+dzlay(i-1))/2.D0
288     csto = 2.D0*rho_aer*effg(zlev(i))/(9.D0*VISAIR(btemp(i)))
289  write(*,*)'yupitururu3'
290     ! ***** WARNING ******
291     ! The 1st order slip flow correction is hidden here :
292     !    (rc(i-1)**ar1 * afun(af1) + cslf/rb2ra * rc(i-1)**ar2 * afun(af2))
293     ! wth is simple the classical settling velocity of particle due to gravity (Stokes law).
294     ! It may be enhanced using a Taylor series of the slip flow correction around a given (modal) radius
295     cslf = akn * FPLAIR(btemp(i),plev(i))
296     wth(i) = -csto/(rb2ra*alpha_k(k,sig)) * (rc(i-1)**ar1 * alpha_k(af1,sig) + cslf/rb2ra * rc(i-1)**ar2 * alpha_k(af2,sig))
297     ! now correct velocity to reduce numerical diffusion
298     IF (.NOT.no_fiadero_w) THEN
299        IF (mk(i) <= 0.D0) THEN
300           ratio = fiadero_max
301        ELSE
302           ratio = MAX(MIN(mk(i-1)/mk(i),fiadero_max),fiadero_min)
303        ENDIF
304        ! apply correction
305        IF ((ratio <= 0.9D0 .OR. ratio >= 1.1D0) .AND. wth(i) /= 0.D0) THEN
306           wdt = wth(i)*dt
307           zcorf(i) = dzb/wdt * (exp(-wdt*log(ratio)/dzb)-1.D0) / (1.D0-ratio)
308        ENDIF
309     ENDIF
310  ENDDO
311  ! last value (ground) set to first layer value: arbitrary :)
312  wth(i) = wth(i-1)
313  zcorf(i) = zcorf(i-1)
314  IF (PRESENT(corf)) corf(:) = zcorf(:)
315  RETURN
316END SUBROUTINE get_weff
317
318
319!*****************************************************************************
320SUBROUTINE aer_sedimentation(dt,nlay,m0,m3,plev,zlev,dzlay,btemp,dm0,dm3,aer_flux)
321  !! Interface to sedimentation algorithm.
322  !!
323  !! The subroutine computes the evolution of each moment of the aerosols tracers
324  !! through sedimentation process and returns their tendencies for a timestep.
325  !!
326  !! It is assumed that the aerosols size-distribution is the same as the mode 1 of cloud drops.
327 
328  use free_param, only: rho_aer, sig_aer
329  use donnees, only: pi, wsed_m0, wsed_m3
330 
331  IMPLICIT NONE
332
333  REAL, INTENT(in)  :: dt
334  INTEGER, INTENT(in) :: nlay
335  REAL, INTENT(in), DIMENSION(nlay)   :: m0
336  !! 0th order moment of the CCN precursors (\(m^{-3}\)).
337  REAL, INTENT(in), DIMENSION(nlay)   :: m3
338  !! 3rd order moment of the CCN precursors (\(m^{3}.m^{-3}\))
339  REAL, INTENT(in), DIMENSION(nlay+1) :: plev
340  !! Pressure at each level (Pa).
341  REAL, INTENT(in), DIMENSION(nlay+1) :: zlev
342  !! Altitude at each level (m).
343  REAL, INTENT(in), DIMENSION(nlay)   :: dzlay
344  !! Altitude at the center of each layer (Pa).
345  REAL, INTENT(in), DIMENSION(nlay)   :: btemp
346  !! Temperature at each level (T).
347  REAL, INTENT(out), DIMENSION(nlay)  :: dm0
348  !! Tendency of the 0th order moment of the CCN precursors (\(m^{-3}\)).
349  REAL, INTENT(out), DIMENSION(nlay)  :: dm3
350  !! Tendency of the 3rd order moment of the CCN precursors (\(m^{3}.m^{-3}\)).
351  REAL, INTENT(out), DIMENSION(nlay+1):: aer_flux
352  !! Aerosol mass (downard) flux
353  REAL, DIMENSION(nlay+1):: ft,fdcor,wth
354  REAL, DIMENSION(nlay)  :: m0vsed,m3vsed,r0
355  REAL                   :: m,n,p
356  REAL, PARAMETER        :: fac = 4.D0/3.D0 * pi
357  INTEGER                :: nla,nle,i
358 
359  aer_flux(:) = 0.D0
360  nla = size(m0,dim=1)
361  nle = nla+1
362
363  DO i=1,nlay
364     r0(i) = get_r0_1(m0(i),m3(i))
365  ENDDO
366
367  ! Compute sedimentation process based on control flag
368  ! wsed_m0 = true, force all aerosols moments to fall at M0 settling velocity
369  ! use M0ccn to compute settling velocity
370  IF (wsed_m0) THEN 
371     ! M0
372     call get_weff(dt,nlay,plev,zlev,dzlay,btemp,m0,0.D0,r0,sig_aer,wth,fdcor)
373     ft(:) = wth(:) * fdcor(:)  ; m0vsed(:) = ft(1:nla)
374     call let_me_fall_in_peace(dt,nlay,dzlay,m0,-ft,dm0)
375     ! M3
376     m3vsed(:) = ft(1:nla)
377     call let_me_fall_in_peace(dt,nlay,dzlay,m3,-ft,dm3)
378     ! get sedimentation flux
379     aer_flux(:) = fac * rho_aer * ft(:) * dm3
380     write(*,*)'sedimentation- aerosol flux: ', aer_flux
381
382  ! wsed_m3 = true, force all aerosols moments to fall at M3 settling velocity
383  ! use M3ccn to compute settling velocity
384  ELSEIF (wsed_m3) THEN
385     ! M3
386     call get_weff(dt,nlay,plev,zlev,dzlay,btemp,m3,3.D0,r0,sig_aer,wth,fdcor)
387     ft(:) = wth(:) * fdcor(:)  ; m3vsed(:) = ft(1:nla)
388     call let_me_fall_in_peace(dt,nlay,dzlay,m3,-ft,dm3)
389     ! M0
390     m0vsed(:) = ft(1:nla)
391     call let_me_fall_in_peace(dt,nlay,dzlay,m0,-ft,dm0)
392     ! get sedimentation flux
393     aer_flux(:) = fac * rho_aer * ft(:) * m3
394     write(*,*)'sedimentation- aerosol flux: ', aer_flux
395
396  ! m_wsed_m0 = false and m_wsed_m3 = false, computes the effective settling velocity for each moments
397  ELSE
398     ! M0
399     call get_weff(dt,nlay,plev,zlev,dzlay,btemp,m0,0.D0,r0,sig_aer,wth,fdcor)
400     ft(:) = wth(:) * fdcor(:)  ; m0vsed(:) = ft(1:nla)
401     call let_me_fall_in_peace(dt,nlay,dzlay,m0,-ft,dm0)
402     ! M3
403     call get_weff(dt,nlay,plev,zlev,dzlay,btemp,m3,3.D0,r0,sig_aer,wth,fdcor)
404     ft(:) = wth(:) * fdcor(:)  ; m3vsed(:) = ft(1:nla)
405     call let_me_fall_in_peace(dt,nlay,dzlay,m3,-ft,dm3)
406     ! get sedimentation flux
407     aer_flux(:) = fac * rho_aer * ft(:) * m3
408     write(*,*)'sedimentation- aerosol flux: ', aer_flux
409  ENDIF
410  RETURN
411END SUBROUTINE aer_sedimentation
412
413
414!*****************************************************************************
415SUBROUTINE drop_sedimentation(dt,nlay,M0mode,M3mode,plev,zlev,dzlay,btemp,mode,dm0ccn,dm0,dm3ccn,dm3liq)
416  !! Interface to sedimentation algorithm.
417  !!
418  !! The subroutine computes the evolution of each moment of the cloud drops
419  !! through sedimentation process and returns their tendencies for a timestep.
420  !!
421  !! m0mode,m3ccn should be vector defined on the vertical structure (top to ground).
422  !! m3liq a 2D array with in 1st dimension the vertical distribution (top to ground) and
423  !! in the 2nd the list of condensate.
424  !!
425  !! The global variables wsed_m0 and wsed_m3 control how the settling velocity is
426  !! applied during the sedimentation process:
427  !!   - If both are .false., m0mode and the sum of m3 (ccn+liquid) are treated separatly.
428  !!   - If mm_wsed_m0 is set to .true., all moments fall at m0mode settling velocity.
429  !!   - Otherwise all moments fall at m3sum settling velocity.
430  !!
431  !! note:
432  !!  drop sedimentation flux is not computed here. It can be done as in aer_sedimentation.
433  !!  However the overall density of cloud drop should be used (that is the density of
434  !!  CCN + the density of each liquid condensate)
435  !!     cld_flux = fac *ft * (m3ccn * rho_aer + SUM(rho_liq(:) + m3liq(:)))
436 
437  use free_param, only: sigma1, sigma2
438  use donnees, only: pi, wsed_m0, wsed_m3
439 
440  IMPLICIT NONE
441
442  REAL, INTENT(in)  :: dt
443  INTEGER, INTENT(in) :: nlay
444  REAL, INTENT(in), DIMENSION(nlay,2) :: M0mode
445  !! 0th order moment (\(m^{-3}\)).
446  REAL, INTENT(in), DIMENSION(nlay,3) :: M3mode
447  !! 3rd order moment (\(m^{-3}\)).
448  REAL, INTENT(in), DIMENSION(nlay+1) :: plev
449  !! Pressure at each level (Pa).
450  REAL, INTENT(in), DIMENSION(nlay+1) :: zlev
451  !! Altitude at each level (m).
452  REAL, INTENT(in), DIMENSION(nlay)   :: dzlay
453  !! Altitude at the center of each layer (Pa).
454  REAL, INTENT(in), DIMENSION(nlay)   :: btemp
455  !! Temperature at each level (T).
456  INTEGER, INTENT(in)                 :: mode
457  !! Size-distribution mode.
458  REAL, INTENT(out), DIMENSION(nlay)  :: dm0
459  !! Tendency of the 0th order moment of droplets (\(m^{-3}\)).
460  REAL, INTENT(out), DIMENSION(nlay)  :: dm0ccn
461  !! Tendency of the 0th order moment of the CCN (\(m^{-3}\)).
462  REAL, INTENT(out), DIMENSION(nlay)  :: dm3ccn
463  !! Tendency of the 3rd order moment of the CCN (\(m^{3}.m^{-3}\)).
464  REAL, INTENT(out), DIMENSION(nlay,2):: dm3liq
465  !! Tendency of the 3rd order moments of each condensate (\(m^{3}.m^{-3}\)).
466  REAL, DIMENSION(nlay+1):: ft, fdcor, wth
467  REAL, DIMENSION(nlay)  :: r0, m3sum
468  REAL                   :: m, n, p, sig
469  REAL, PARAMETER        :: fac = 4.D0/3.D0 * pi
470  INTEGER                :: i, nla, nle, nc
471
472  nla = nlay
473  nle = nla + 1 ; nc = nlay
474
475  m3sum(:) = 0.D0
476  DO i=1,nla
477     m3sum(:) = M3mode(i,1) + M3mode(i,2) + M3mode(i,3) ! liq1 + liq2 + ccn
478  ENDDO
479
480  IF (mode == 2) THEN
481     sig = sigma2
482     DO i=1,nla
483        r0(i) = get_r0_2(M0mode(i,1),m3sum(i))
484     ENDDO
485  ELSE
486     sig = sigma1
487     DO i=1,nla
488        r0(i) = get_r0_1(M0mode(i,1),m3sum(i))
489     ENDDO
490  ENDIF
491
492  ! Compute sedimentation process based on control flag
493  ! wsed_m0 = true, force all aerosols moments to fall at M0 settling velocity
494  ! use M0mode to compute settling velocity
495  IF (wsed_m0) THEN
496     ! use M0mode to compute settling velocity
497     call get_weff(dt,nlay,plev,zlev,dzlay,btemp,M0mode(:,1),0.D0,r0,sig,wth,fdcor)
498     ft(:) = wth(:) * fdcor(:)
499
500     ! apply the sedimentation process to all other moments
501     call let_me_fall_in_peace(dt,nlay,dzlay,M0mode(:,1),-ft,dm0)
502     call let_me_fall_in_peace(dt,nlay,dzlay,M0mode(:,2),-ft,dm0ccn)
503     call let_me_fall_in_peace(dt,nlay,dzlay,M3mode(:,3),-ft,dm3ccn)
504     call let_me_fall_in_peace(dt,nlay,dzlay,M3mode(:,2),-ft,dm3liq(:,2))
505     call let_me_fall_in_peace(dt,nlay,dzlay,M3mode(:,1),-ft,dm3liq(:,1))
506
507  ! Compute sedimentation process based on control flag
508  ! wsed_m3 = true, force all aerosols moments to fall at M3 settling velocity
509  ! use M3ccn to compute settling velocity
510  ELSEIF (wsed_m3) THEN
511     ! use M3ccn + M3liq = m3sum to compute settling velocity
512     call get_weff(dt,nlay,plev,zlev,dzlay,btemp,m3sum,3.D0,r0,sig,wth,fdcor)
513     ft(:) = wth(:) * fdcor(:)
514
515     ! apply the sedimentation process to all other moments
516     call let_me_fall_in_peace(dt,nlay,dzlay,M0mode(:,1),-ft,dm0)
517     call let_me_fall_in_peace(dt,nlay,dzlay,M0mode(:,2),-ft,dm0ccn)
518     call let_me_fall_in_peace(dt,nlay,dzlay,M3mode(:,3),-ft,dm3ccn)
519     call let_me_fall_in_peace(dt,nlay,dzlay,M3mode(:,2),-ft,dm3liq(:,2))
520     call let_me_fall_in_peace(dt,nlay,dzlay,M3mode(:,1),-ft,dm3liq(:,1))
521
522  ELSE
523     ! M0mode and M3sum fall independently
524     ! M0
525     call get_weff(dt,nlay,plev,zlev,dzlay,btemp,m0mode,0.D0,r0,sig,wth,fdcor)
526     ft(:) = wth(:) * fdcor(:)
527     call let_me_fall_in_peace(dt,nlay,dzlay,M0mode(:,1),-ft,dm0)
528     call let_me_fall_in_peace(dt,nlay,dzlay,M0mode(:,2),-ft,dm0ccn)
529
530     ! M3
531     call get_weff(dt,nlay,plev,zlev,dzlay,btemp,m3sum,3.D0,r0,sig,wth,fdcor)
532     ft(:) = wth(:) * fdcor(:)
533     call let_me_fall_in_peace(dt,nlay,dzlay,M3mode(:,3),-ft,dm3ccn)
534     call let_me_fall_in_peace(dt,nlay,dzlay,M3mode(:,2),-ft,dm3liq(:,2))
535     call let_me_fall_in_peace(dt,nlay,dzlay,M3mode(:,1),-ft,dm3liq(:,1))
536  ENDIF
537
538  RETURN
539END SUBROUTINE drop_sedimentation
540
541END MODULE SED_AND_PROD_MAD
Note: See TracBrowser for help on using the repository browser.