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

Last change on this file since 2203 was 1676, checked in by emillour, 8 years ago

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