source: trunk/LMDZ.VENUS/libf/phyvenus/cloudvenus/mad_muphy.F90 @ 3094

Last change on this file since 3094 was 1667, checked in by slebonnois, 8 years ago

SL: an other round of bug corrections for Venus clouds
File size: 22.6 KB
RevLine 
[1661]1!-----------------------------------------------------------------------
2! S.GUILBON - LATMOS - 2016.
3!-----------------------------------------------------------------------
4! PROGRAM MAD contains
5!
6!  SUBROUTINE   MAD_MUPHY        Microphysical processes with moments
7!  SUBROUTINE   change_layer     Need if layer changes
8!  SUBROUTINE   change_wsa       Need if wsa changes
9!  SUBROUTINE   change_vapor     Need if number concentration of vapor changes
10!  FUNCTION     logn
11!  FUNCTION     moment
12
13!------------------------------------------------------------------------------------------------
14!------------------------------------------------------------------------------------------------
15SUBROUTINE MAD_MUPHY(dt,TAIR,PAIR,MRWV_loc,MRSA_loc,M0_aeros,M3_aeros,  &
16                           M0_m1drop,M0_m1ccn,M3_m1sa,M3_m1w,M3_m1ccn,  &
17                           M0_m2drop,M0_m2ccn,M3_m2sa,M3_m2w,M3_m2ccn)
18  ! Loop on lat and long outside this subroutine => in phytrac_chimie
19!------------------------------------------------------------------------------------------------
20!------------------------------------------------------------------------------------------------
21
22  use free_param
23  use donnees
24
25  IMPLICIT NONE
26
27  integer :: i, mode
28
29  ! LAYER
30  real, intent(in) :: TAIR, PAIR, dt
31  real, intent(inout) :: MRWV_loc, MRSA_loc
32
33  ! MOMENTS
34  real, intent(inout):: M0_aeros, M3_aeros
35  real, intent(inout):: M0_m1drop,M0_m1ccn,M0_m2drop,M0_m2ccn
36  real, intent(inout):: M3_m1sa,M3_m1w,M3_m1ccn,M3_m2sa,M3_m2w,M3_m2ccn
37
38  ! MODAL TENDENCIES
39  real, dimension(2) :: M0_m1, M0_m2
40  real, dimension(3) :: M3_m1, M3_m2
41  real :: p, k                        ! Modal order
42  real :: sum_M0, sum_M3              ! Sum of moment
43  real :: dM0_merge_m1, dM3_merge_m1  ! Mode merging mode 1
44  real :: dM0_merge_m2, dM3_merge_m2  ! Mode merging mode 2
45  real :: dM0_flux_m1,  dM0_flux_m2   ! Masse flux, only M0, mode 1 and 2
46  real :: dM3_flux_m1,  dM3_flux_m2   ! Masse flux, only M3, mode 1 and 2
47  real :: dM0_aeros, dM3_aeros, ra    ! Heterogeneous nucleation, only aerosol
48  real :: dM0_hom, dM3_hom            ! Homogeneous nucleation, only mode 1
49  real :: dM0_het, dM3_het            ! Heterogeneous nucleation, only mode 1
50  real :: dM0_m1_drop, dM0_m1_ccn, dM3_m1_SA, dM3_m1_WV, dM3_m1_ccn, dM0_aeros_m2
51  real :: dM0_m2_drop, dM0_m2_ccn, dM3_m2_SA, dM3_m2_WV, dM3_m2_ccn, dM3_aeros_m2
52  real :: dM3_m1, dM3_m2, dM3, dM3_SA1, dM3_SA2, dM3_WV1, dM3_WV2, dM3_SA, dM3_WV
53  real :: frac_ccn_m1, frac_ccn_m2
54
55  ! CONSERVATION Variables
56  real :: var_mode1_M0, var_mode1_M3, var_mode2_M0, var_mode2_M3
57  real :: var_MRSA_m1, var_MRWV_m1
58  real :: var_MRSA_m2, var_MRWV_m2
59  real :: var_MR_m1, var_MR_m2, var_MRSA, var_MRWV
60  real :: factor, diff, check, MRTOT
61
62  ! FUNCTIONS
63  real :: moment, D0, alpha_k
64
65  ! OUTPUTS and GRAPHES
66  real, dimension(nbin) :: n_g1
67  integer :: cpt
68
69  !-------------
70
71  ! Initialization of vectors
72  M0_m1(1)=M0_m1drop
73  M0_m1(2)=M0_m1ccn
74  M0_m2(1)=M0_m2drop
75  M0_m2(2)=M0_m2ccn
76  M3_m1(1)=M3_m1sa
77  M3_m1(2)=M3_m1w
78  M3_m1(3)=M3_m1ccn
79  M3_m2(1)=M3_m2sa
80  M3_m2(2)=M3_m2w
81  M3_m2(3)=M3_m2ccn
82 
83  !-------------
84
85  ! Initialization of tendencies   
86  dM0_hom = 0.D0      ;  dM3_hom = 0.D0
87  dM0_het = 0.D0      ;  dM3_het = 0.D0
88  dM0_aeros = 0.D0    ;  dM3_aeros = 0.D0 
89  dM0_flux_m1 = 0.D0  ;  dM0_flux_m2 = 0.D0
90  dM3_flux_m1 = 0.D0  ;  dM3_flux_m2 = 0.D0
91  dM0_merge_m1 = 0.D0 ;  dM3_merge_m1 = 0.D0
92  dM0_merge_m2 = 0.D0 ;  dM3_merge_m2 = 0.D0
93
94  !-------------
95
96  ! To controle the conservation !
97  ! Total mass (vapor + cond) of species (h2o + h2so4)
98  sum_M3 = M3_m1(1) + M3_m1(2) + M3_m1(3) ! mode 1
99  MRTOT = sum_M3/(exp(4.5D0*log(sigma1)))
100  sum_M3 = M3_m2(1) + M3_m2(2) + M3_m2(3) ! mode 2
101  MRTOT = MRTOT + sum_M3/(exp(4.5D0*log(sigma2)))
102  MRTOT = MRTOT*(4.0D0/3.0D0)*PI*RHOSA    ! total condensed mass of SA
103  MRTOT = MRTOT*WSA/(MSA*E) + MRTOT*(1.D0-WSA)/(MWV*E) + MRSA_loc + MRWV_loc ! in vapor quantity
104  check = MRTOT
105
106!-----------------------------------------------------------------------
107!                    INITIALIZATION - THERMO
108!-----------------------------------------------------------------------
109  WSA = 0.01D0    ! Init. of the fraction of sulfuric acid [0.1;1]
110  WSAEQ = 0.0D0   ! WSA at equilibrium
111
112  call change_wsa(TAIR,PAIR,MRSA_loc)            ! Need if wsa changes
113  call change_layer(TAIR,PAIR)                   ! Need if layer changes
114  call change_vapor(TAIR,PAIR,MRSA_loc,MRWV_loc) ! Need if nb vapor concentration changes
115
116  ! Number concentration of H2SO4 (m-3)
117  h2so4_m3 = MRSA_loc*PAIR/(KBZ*TAIR)
118
119!-----------------------------------------------------------------------
120!                           MODE MERGING
121!-----------------------------------------------------------------------
122  IF (SEDIM .and. MERGE) THEN
123     redge = (r1*r2)**0.5D0 ! Edge radius of mode 1
[1667]124!     write(*,*)'The intersection size or virtual bonduary', redge
[1661]125
126     IF (r1 .EQ. redge) THEN
[1667]127!        write(*,*)'Mode 1 merge to mode 2', redge, r1
[1661]128        CALL MERGING(r1,N1,sigma1,sigma2,dM0_merge_m1,dM3_merge_m1, &
129             &   dM0_merge_m2,dM3_merge_m2)
130        ! mode 1
131        M0_m1(1) = dM0_merge_m1
132        M0_m1(2) = dM0_merge_m1*frac_ccn_m1
133        M3_m1(1) = dM3_merge_m1*(1.D0-frac_ccn_m1)*WSA
134        M3_m1(2) = dM3_merge_m1*(1.D0-frac_ccn_m1)*(1.D0-WSA)
135        M3_m1(3) = dM3_merge_m1*frac_ccn_m1
136        ! mode 2
137        M0_m2(1) = M0_m2(1) + dM0_merge_m2
138        M0_m2(2) = M0_m2(2) + dM0_merge_m2*(frac_ccn_m2+frac_ccn_m1)
139        M3_m2(1) = M3_m2(1) + dM3_merge_m2*(1.D0-(frac_ccn_m2+frac_ccn_m1))*WSA
140        M3_m2(2) = M3_m2(2) + dM3_merge_m2*(1.D0-(frac_ccn_m2+frac_ccn_m1))*(1.D0-WSA)
141        M3_m2(3) = M3_m2(3) + dM3_merge_m2*(frac_ccn_m2+frac_ccn_m1)
142     ENDIF
143
144     IF (r2 .EQ. redge) THEN
[1667]145!        write(*,*)'Mode 2 merge to mode 1',redge,r2
[1661]146        CALL MERGING(r2,N2,sigma2,sigma1,dM0_merge_m2,dM3_merge_m2, &
147             &   dM0_merge_m1,dM3_merge_m1)
148        ! mode 1
149        M0_m2(1) = dM0_merge_m2
150        M0_m2(2) = dM0_merge_m2*frac_ccn_m2
151        M3_m2(1) = dM3_merge_m2*(1.D0-frac_ccn_m2)*WSA
152        M3_m2(2) = dM3_merge_m2*(1.D0-frac_ccn_m2)*(1.D0-WSA)
153        M3_m2(3) = dM3_merge_m2*frac_ccn_m2
154        ! mode 2
155        M0_m1(1) = M0_m1(1) + dM0_merge_m1
156        M0_m1(2) = M0_m1(2) + dM0_merge_m1*(frac_ccn_m2+frac_ccn_m1)
157        M3_m1(1) = M3_m1(1) + dM3_merge_m1*(1.D0-(frac_ccn_m2+frac_ccn_m1))*WSA
158        M3_m1(2) = M3_m1(2) + dM3_merge_m1*(1.D0-(frac_ccn_m2+frac_ccn_m1))*(1.D0-WSA)
159        M3_m1(3) = M3_m1(3) + dM3_merge_m1*(frac_ccn_m2+frac_ccn_m1)
160     ENDIF
161  ENDIF
162
163
164!-----------------------------------------------------------------------
165!                           NUCLEATION                       Only mode 1
166!-----------------------------------------------------------------------
167  IF(NUCLEA) then ! HOmogeneous NUcleation with MOment
168     CALL HONUMO(TAIR,dt,M0_m1,M3_m1,dM0_hom,dM3_hom) !#.m(air)-3 and m3.m(air)-3
169  ENDIF
170
171  IF(HETNUCLEA) then ! HETerogeneous NUcleation with MOment
172     ! Calculation of ra (mean radius of aerosols) before heterogeneous nucleation
173     ! Needed ra to calculate the bin idstrib to determine Jhet.
174     IF (M0_aeros .le. ZERO) THEN ! Mean radii MODE AER
175        ra = 0.D0
176     ELSE
177        ra = (1.D0/alpha_k(3,sig_aer) * M3_aeros/M0_aeros)**(1.D0/3.D0)     
178     ENDIF
179
180     CALL HETNUMO(TAIR,dt,ra,M0_aeros,M3_aeros,M0_m1,M3_m1, &
181          & dM0_aeros,dM3_aeros,dM0_het,dM3_het)
182  ENDIF
183
184
185!-----------------------------------------------------------------------
186!                           MASSES FLUX                       Mode 1 & 2
187!-----------------------------------------------------------------------
188  ! We can have condensation/evaporation if we have some particles :)
189  ! Conditions : flag + mode + gas
190  IF (MASSFLUX .and. (M0_m1(1).gt.ZERO)) THEN ! Mode 1
191     CALL FLUX(TAIR,PAIR,dt,sigma1,r1,M3_m1,dM0_flux_m1,dM3_flux_m1)
192  ENDIF
193
194  IF (MASSFLUX .and. (M0_m2(1).gt.ZERO))  THEN ! Mode 2
195     CALL FLUX(TAIR,PAIR,dt,sigma2,r2,M3_m2,dM0_flux_m2,dM3_flux_m2)
196  ENDIF
197
198
199!-----------------------------------------------------------------------
200!                           CONSERVATION                      Mode 1 & 2
201!-----------------------------------------------------------------------
202  ! M3(1) sulfuric acid
203  ! M3(2) water
204  ! M3(3) condensed cloud nuclei
205  !-------------------------------------!
206
207
208  !--------AEROSOLS / HETERONUC---------!
209  ! We just modify the tendencies. They will be add to moments.
210  IF ((M0_aeros + dM0_aeros) .le. ZERO) THEN
211     dM0_aeros = (-1.D0) * M0_aeros
212     dM3_aeros = (-1.D0) * M3_aeros
213     dM0_het = M0_aeros
214     dM3_het = M3_aeros
215  ENDIF
216  !-------------------------------------!
217
218
219  !----------SUM OF TENDENCIES----------!
220  ! +5% of ccn vol to condensed mass to active the mass flux in next time step.
221  ! Mode 1
222  dM0_m1_drop = dM0_hom + dM0_het + dM0_flux_m1
223  dM0_m1_ccn = dM0_het
224  dM3_m1_SA = (dM3_hom + dM3_het + dM3_flux_m1) * WSA
225  dM3_m1_WV = (dM3_hom + dM3_het + dM3_flux_m1) * (1.D0-WSA)
226  dM3_m1_ccn = dM3_het
227  ! Mode 2
228  dM0_m2_drop = dM0_flux_m2
229  dM0_m2_ccn = 0.D0
230  dM3_m2_SA = dM3_flux_m2 * WSA
231  dM3_m2_WV = dM3_flux_m2 * (1.D0-WSA)
232  dM3_m2_ccn = 0.D0
233  ! Aerosols and CCN
234  dM0_aeros_m2 = 0.D0
235  dM3_aeros_m2 = 0.D0
236  !-------------------------------------!
237
238
239  !--------TOO MUCH EVAPORATION---------!
240  ! We just modify the tendencies. They will be add to moments.
241  IF ((M3_m1(1) + dM3_m1_SA) .le. ZERO) THEN
242     dM3_m1_SA = (-1.D0) * M3_m1(1) !-- vol acid
243  ENDIF
244  IF ((M3_m1(2) + dM3_m1_WV) .le. ZERO) THEN
245     dM3_m1_WV = (-1.D0) * M3_m1(2) !-- vol water
246  ENDIF
247  IF ((M0_m1(1) + dM0_m1_drop) .le. ZERO) THEN
248     dM0_m1_drop = (-1.D0)* M0_m1(1) !-- nb droplets
249     dM0_m1_ccn = (-1.D0) * M0_m1(2) !-- nb ccn
250     dM3_m1_SA  = (-1.D0) * M3_m1(1) !-- acid
251     dM3_m1_WV  = (-1.D0) * M3_m1(2) !-- water
252     dM3_m1_ccn = (-1.D0) * M3_m1(3) !-- ccn
253     dM0_aeros = M0_m1(2)        !++ nb ccn
254     dM3_aeros = M3_m1(3)        !++ vol ccn
255  ENDIF
256  ! Second mode of particles !
257  IF ((M3_m2(1) + dM3_m2_SA) .le. ZERO) THEN
258     dM3_m2_SA = (-1.D0) * M3_m2(1) !-- vol acid
259  ENDIF
260  IF ((M3_m2(2) + dM3_m2_WV) .le. ZERO) THEN
261     dM3_m2_WV = (-1.D0) * M3_m2(2) !-- vol water
262  ENDIF
263  IF ((M0_m2(1) + dM0_m2_drop) .le. ZERO) THEN
264     dM0_m2_drop = (-1.D0)* M0_m2(1) !-- nb droplets
265     dM0_m2_ccn = (-1.D0) * M0_m2(2) !-- nb ccn
266     dM3_m2_SA  = (-1.D0) * M3_m2(1) !-- acid
267     dM3_m2_WV  = (-1.D0) * M3_m2(2) !-- water
268     dM3_m2_ccn = (-1.D0) * M3_m2(3) !-- ccn
269     dM0_aeros_m2 = M0_m2(2)         !++ nb ccn
270     dM3_aeros_m2 = M3_m2(3)         !++ vol ccn
271  ENDIF
272  !-------------------------------------!
273
274
275  !--------TOO MUCH CONDENSATION--------!
276  ! SUM with new tendencies
277  dM3_m1 = dM3_m1_SA + dM3_m1_WV
278  dM3_m2 = dM3_m2_SA + dM3_m2_WV
279  ! Vapor mode 1
280  dM3 = dM3_m1 / (exp(4.5D0*log(sigma1)))
281  dM3 = RHOSA * 4.D0/3.D0 * PI * dM3_m1 ! mass
282  dM3_SA1 = (dM3_m1*WSA)/(MSA*E)                     ! Total consumed SA vapor
283  dM3_WV1 = (dM3_m1*(1.D0-WSA))/(MWV*E)              ! Total consumed water vapor
284  ! Vapor mode 2
285  dM3 = dM3_m2/(exp(4.5D0*log(sigma2)))
286  dM3 = RHOSA * 4.D0/3.D0 * PI * dM3_m2 ! mass
287  dM3_SA2 = (dM3_m2*WSA)/(MSA*E)           ! Total consumed SA vapor
288  dM3_WV2 = (dM3_m2*(1.D0-WSA))/(MWV*E)    ! Total consumed water vapor
289  ! Sum with distinguished mode: not same process for mode 1 and mode 2
290  dM3_SA = dM3_SA1 + dM3_SA2 ! h2so4
291  dM3_WV = dM3_WV1 + dM3_WV2 ! h2o
292  ! We just modify the tendencies. They will be add to moments.
293  ! we define a sort of purcentage of participation for each moment.
294  IF (((MRSA_loc + dM3_SA) .le. ZERO).and.(dM3_SA.ne.0.D0).and.(WSA.ne.0.D0)) THEN ! Comsuption of all H2SO4 vapor
295     dM3_m1_SA = dM3_SA1/dM3_SA * MRSA_loc * MSA*E/WSA * 1.D0/(RHOSA*4.D0/3.D0*PI)*alpha_k(3,sigma1)
296     dM3_m2_SA = dM3_SA2/dM3_SA * MRSA_loc * MSA*E/WSA * 1.D0/(RHOSA*4.D0/3.D0*PI)*alpha_k(3,sigma2)
297     dM3_SA = dM3_SA1 + dM3_SA2
298  ENDIF
299  IF (((MRWV_loc + dM3_WV) .le. ZERO).and.(dM3_WV.ne.0.D0)) THEN  ! Comsuption of all H2O vapor
300     dM3_m1_WV = dM3_WV1/dM3_WV * MRWV_loc * MWV*E/(1.D0-WSA) * 1.D0/(RHOSA*4.D0/3.D0*PI)*alpha_k(3,sigma1)
301     dM3_m2_WV = dM3_WV2/dM3_WV * MRWV_loc * MWV*E/(1.D0-WSA) * 1.D0/(RHOSA*4.D0/3.D0*PI)*alpha_k(3,sigma2)
302     dM3_WV = dM3_WV1 + dM3_WV2
303  ENDIF
304  !------------------------------------!
305
306
307  !-------------UPDATES----------------!
308  ! Updates Vapors !
309  MRSA_loc = MRSA_loc + dM3_SA ! acid
310  MRWV_loc = MRWV_loc + dM3_WV ! water
311  ! Updtaes of CCN and Aerosols
312  M0_aeros = M0_aeros + dM0_aeros +  dM0_aeros_m2     ! number of aerosols
313  M3_aeros = M3_aeros + dM3_aeros +  dM3_aeros_m2     ! vol/mass of aerosols 
314  ! Updates moments - Mode 1
315  M0_m1(1) = M0_m1(1) + dM0_m1_drop ! number of droplets
316  M0_m1(2) = M0_m1(2) + dM0_m1_ccn  ! number of CCN in droplets
317  M3_m1(1) = M3_m1(1) + dM3_m1_SA   ! vol/mass of Héso4 in droplets
318  M3_m1(2) = M3_m1(2) + dM3_m1_WV   ! vol/mass of h2o in droplets
319  M3_m1(3) = M3_m1(3) + dM3_m1_ccn  ! vol/mass of CCN in droplets
320  ! Updates moments - Mode 2
321  M0_m2(1) = M0_m2(1) + dM0_m2_drop
322  M0_m2(2) = M0_m2(2) + dM0_m2_ccn
323  M3_m2(1) = M3_m2(1) + dM3_m2_SA
324  M3_m2(2) = M3_m2(2) + dM3_m2_WV
325  M3_m2(3) = M3_m2(3) + dM3_m2_ccn
326  !------------------------------------!
[1667]327!  write(*,*)'l.653: after up', M3_m1(1),M3_m1(2),M3_m1(3),dM3_m1_SA,dM3_m1_WV,dM3_m1_ccn
328!  write(*,*)'l.654: after up', M0_m1(1),M0_m1(2),dM0_m1_drop
[1661]329!-----------------------------------------------------------------------
330  ! To controle/check the mass conservation !
331  MRTOT = (M3_m1(1)+M3_m1(2))/(exp(4.5D0*log(sigma1)))         ! mode 1
332  MRTOT = MRTOT + (M3_m2(1)+M3_m2(2))/(exp(4.5D0*log(sigma2))) ! mode 2
333  MRTOT = MRTOT * (4.0D0/3.0D0)*PI*RHOSA                       ! mass
334  MRTOT = MRTOT*WSA/(MSA*E) + MRTOT*(1.D0-WSA)/(MWV*E)         ! binary solution !     
335  MRTOT = MRTOT + MRSA_loc + MRWV_loc                          ! with vapor gas
336  diff = check - MRTOT                                         ! difference ?
337
338  if (diff .gt. ZERO) then
339     write(*,*) 'MASS CHANGEs !!!', diff
340     stop
341  endif
342
343  ! Updates
344  call change_wsa(TAIR,PAIR,MRSA_loc)
345  call change_vapor(TAIR,PAIR,MRSA_loc,MRWV_loc)
346
347
348!-----------------------------------------------------------------------
349!                           COAGULATION
350!-----------------------------------------------------------------------
351!!$  IF (COAG) THEN
352!!$     CALL drop_coagul(TAIR,PAIR,dt,M0_m1,M3_m1,M0_m2,M3_m2)
353!!$  END IF
354
355
356!-----------------------------------------------------------------------
357!                        R_mean AND N_tot
358!-----------------------------------------------------------------------
359 
360  IF (M0_m1(1) .le. ZERO) THEN ! Mean radii MODE 1
361     r1 = 0.D0
362     N1 = 0.D0
363     M0_m1(1) = 0.D0
364     M0_aeros = M0_aeros + M0_m1(2)
365     M0_m1(2) = 0.D0
366     M3_m1(1) = 0.D0
367     M3_m1(2) = 0.D0
368     M3_aeros = M3_aeros + M3_m1(3)
369     M3_m1(3) = 0.D0
370  ELSE
371     r1 = (1.D0/alpha_k(3,sigma1) * (M3_m1(1)+M3_m1(2)+M3_m1(3))/M0_m1(1))**(1.D0/3.D0)
372     N1 = M0_m1(1)
373     IF ((r1 .lt. ZERO) .or. (N1.lt.ZERO)) THEN
374        r1 = 0.D0
375        N1 = 0.D0
376     ENDIF
[1667]377!     write(*,*)'rmean not ZERO',r1,N1,M0_m1(1),M3_m1(1),WSA
[1661]378  ENDIF
379 
380  IF (M0_m2(1) .le. ZERO) THEN ! Mean radii MODE 2
381     r2 = 0.D0
382     N2 = 0.D0
383     M0_m2(1) = 0.D0
384     M0_aeros = M0_aeros + M0_m2(2)
385     M0_m2(2) = 0.D0
386     M3_m2(1) = 0.D0
387     M3_m2(2) = 0.D0
388     M3_aeros = M3_aeros + M3_m2(3)
389     M3_m2(3) = 0.D0
390  ELSE
391     r2 = (1.D0/alpha_k(3,sigma2) * (M3_m2(1)+M3_m2(2)+M3_m2(3))/M0_m2(1))**(1.D0/3.D0)
392     N2 = M0_m2(1)
393     IF ((r2 .lt. ZERO) .or. (N2.lt.ZERO)) THEN
394        r2 = 0.D0
395        N2 = 0.D0
396     ENDIF
397  ENDIF
398
399  IF (M0_aeros .le. ZERO) THEN ! Mean radii MODE AER
400     ra = 0.D0
401  ELSE
402     ra = (1.D0/alpha_k(3,sig_aer) * M3_aeros/M0_aeros)**(1.D0/3.D0)     
403  ENDIF
404
405  NTOT = N1 + N2
406  ! CCN fraction in distribution => based on mode 1, like WSA ratio !
407  frac_ccn_m1 = M0_m1(2)/(M0_m1(1)+M0_m1(2))
408  frac_ccn_m2 = M0_m2(2)/(M0_m2(1)+M0_m2(2))
409
[1667]410!  write(*,*)'rmean',r1,'NTOT',N1, 'SH2SO4', SH2SO4
[1661]411
412
413!-----------------------------------------------------------------------
414!                           MODE MERGING
415!-----------------------------------------------------------------------
416  IF (MERGE) THEN
417     redge = (r1*r2)**0.5D0 ! Edge radius of mode 1
[1667]418!     write(*,*)'The intersection size or virtual bonduary', redge
[1661]419
420     IF (r1 .EQ. redge) THEN
[1667]421!        write(*,*)'Mode 1 merge to mode 2', redge, r1
[1661]422        CALL MERGING(r1,N1,sigma1,sigma2,dM0_merge_m1,dM3_merge_m1, &
423             &   dM0_merge_m2,dM3_merge_m2)
424        ! mode 1
425        M0_m1(1) = dM0_merge_m1
426        M0_m1(2) = dM0_merge_m1*frac_ccn_m1
427        M3_m1(1) = dM3_merge_m1*(1.D0-frac_ccn_m1)*WSA
428        M3_m1(2) = dM3_merge_m1*(1.D0-frac_ccn_m1)*(1.D0-WSA)
429        M3_m1(3) = dM3_merge_m1*frac_ccn_m1
430        ! mode 2
431        M0_m2(1) = M0_m2(1) + dM0_merge_m2
432        M0_m2(2) = M0_m2(2) + dM0_merge_m2*(frac_ccn_m2+frac_ccn_m1)
433        M3_m2(1) = M3_m2(1) + dM3_merge_m2*(1.D0-(frac_ccn_m2+frac_ccn_m1))*WSA
434        M3_m2(2) = M3_m2(2) + dM3_merge_m2*(1.D0-(frac_ccn_m2+frac_ccn_m1))*(1.D0-WSA)
435        M3_m2(3) = M3_m2(3) + dM3_merge_m2*(frac_ccn_m2+frac_ccn_m1)
436     ENDIF
437
438     IF (r2 .EQ. redge) THEN
[1667]439!        write(*,*)'Mode 2 merge to mode 1',redge,r2
[1661]440        CALL MERGING(r2,N2,sigma2,sigma1,dM0_merge_m2,dM3_merge_m2, &
441             &   dM0_merge_m1,dM3_merge_m1)
[1664]442        ! mode 2
[1661]443        M0_m2(1) = dM0_merge_m2
444        M0_m2(2) = dM0_merge_m2*frac_ccn_m2
445        M3_m2(1) = dM3_merge_m2*(1.D0-frac_ccn_m2)*WSA
446        M3_m2(2) = dM3_merge_m2*(1.D0-frac_ccn_m2)*(1.D0-WSA)
447        M3_m2(3) = dM3_merge_m2*frac_ccn_m2
[1664]448        ! mode 1
[1661]449        M0_m1(1) = M0_m1(1) + dM0_merge_m1
450        M0_m1(2) = M0_m1(2) + dM0_merge_m1*(frac_ccn_m2+frac_ccn_m1)
451        M3_m1(1) = M3_m1(1) + dM3_merge_m1*(1.D0-(frac_ccn_m2+frac_ccn_m1))*WSA
452        M3_m1(2) = M3_m1(2) + dM3_merge_m1*(1.D0-(frac_ccn_m2+frac_ccn_m1))*(1.D0-WSA)
453        M3_m1(3) = M3_m1(3) + dM3_merge_m1*(frac_ccn_m2+frac_ccn_m1)
454     ENDIF
455  ENDIF
456
[1664]457!-----------------------------------------------------------------------
458!                        UPDATE OF OUTPUT VALUES
459!-----------------------------------------------------------------------
[1661]460
[1664]461        ! mode 1
462        M0_m1drop = M0_m1(1)
463        M0_m1ccn  = M0_m1(2)
464        M3_m1sa   = M3_m1(1)
465        M3_m1w    = M3_m1(2)
466        M3_m1ccn  = M3_m1(3)
467        ! mode 2
468        M0_m2drop = M0_m2(1)
469        M0_m2ccn  = M0_m2(2)
470        M3_m2sa   = M3_m2(1)
471        M3_m2w    = M3_m2(2)
472        M3_m2ccn  = M3_m2(3)
473
[1661]474!-----------------------------------------------------------------------
475!                          OUTPUTS => GRAPHES ! :D
476!-----------------------------------------------------------------------
477  if (dt .eq. 1) then
[1667]478!     write(*,*) 'Ntot, rc mode 1: ', N1,r1
[1661]479     ! Log-normal distribution function
480     CALL build_radius_grid(nbin,rmin,rmax,vratio)
481     CALL logdist(sigma1,N1,r1,rad_cld,n_g1)
[1667]482!     DO i = 1, nbin
483!        open(888,file='ND_t3600_nuc-het-flux_dt1_1step_sssat')
484!        write(888,'(I10,2X,3(ES15.7,1X))') i, rad_cld(i), n_g1(i)
485!     ENDDO
[1661]486  endif
487
488
489  RETURN
490END SUBROUTINE MAD_MUPHY
491
492
493!*****************************************************************************
494!                              *************                                 !
495!                               SUBROUTINES                                  !
496!                              *************                                 !
497!*****************************************************************************
498SUBROUTINE change_layer(TAIR,PAIR)
499
500  use free_param
501  use donnees
502
503  IMPLICIT NONE
504
505  real, intent(in) :: TAIR, PAIR
506  real :: FPLAIR, DFWVA, MFPAIR, CDTAIR
507
508  MFPAIR = FPLAIR(TAIR,PAIR) ! Mean free path of air molecules
[1667]509!  write(*,*) 'MFPAIR',MFPAIR, TAIR, PAIR
[1661]510  Kn = MFPAIR/ri1             ! Knudsen number
511
512  D = DFWVA(TAIR,PAIR)      ! Molecular diffusivity of the air (m2/s)
513  Akn = Kn*(1.333D0+0.71D0*(1.0D0/Kn)/(1.0D0+(1.0D0/Kn)))
514  D = D/(1.0D0 + Akn)
515
516  E = PAIR / (KBZ*TAIR)  ! Constant for microphys. processes
517
518  KAIR = CDTAIR(TAIR)        ! Thermal conductivity of air (J/m sec K)
519
520END SUBROUTINE change_layer
521
522
523!*****************************************************************************
524SUBROUTINE change_wsa(TAIR,PAIR,MRSA_loc)
525
526  use free_param
527  use donnees
528
529  IMPLICIT NONE
530
531  real, intent(in) :: TAIR, PAIR, MRSA_loc
532  real :: ROSAS, PSASAS_ZELE,  STSAS
[1667]533!  write(*,*)'checkpoint change_wsa', ST, RHOSA, RHOsasat, PSASAS_ZELE(TAIR,WSA), TAIR, WSA
[1661]534  ST = STSAS(TAIR,WSA)      ! Surface tension of sulfuric acid solution/vapor (N/m)
535
536  RHOSA = ROSAS(TAIR,WSA)          ! Density of liquid sulfuric acid solution (kg/m3)
537  RHOsasat = PSASAS_ZELE(TAIR,WSA) ! Saturated sulfuric acid vapor pressure (Pa)
538
539  SH2SO4 = (PAIR*MRSA_loc)/RHOsasat    ! Saturation Ratio
540
[1667]541!  write(*,*)'checkpoint change_wsa', PSASAS_ZELE(TAIR,WSA)
[1661]542
543END SUBROUTINE change_wsa
544
545
546!*****************************************************************************
547SUBROUTINE change_vapor(TAIR,PAIR,MRSA_loc,MRWV_loc)
548
549  use free_param
550  use donnees
551
552  IMPLICIT NONE
553
554  real, intent(in) :: TAIR, PAIR, MRSA_loc, MRWV_loc
555  real :: waterps
556
557  H2SO4_m3 = MRSA_loc*PAIR/(KBZ*TAIR) ! Number concentration of H2SO4 (m-3)
558
559  PPWV = PAIR*MRWV_loc                   ! Partial pressure of water vapor (Pa or N/m2)
560  PPSA = PAIR*MRSA_loc                   ! Partial Pressure of sulfuric acid (Pa or N/m2)
561
562  RH = PPWV/waterps(TAIR) * 1.0D2    ! Relative humidity (%)
563
564  CHECKWV = PPWV/PAIR
565  CHECKSA = PPSA/PAIR
566
567END SUBROUTINE change_vapor
568
569
570!*****************************************************************************
571FUNCTION logn(x,rm,sig)
572
573  !*     Computes the log-normal distribution at x with parameters rm and sig
574  !*
575  !*     INPUTS
576  !*          x    location
577  !*          rm   mean radius
578  !*          sig  standard deviation
579  !*
580  !*     OUTPUTS
581  !*          The value of the log-normal distribution at location x
582
583  use free_param
584  use donnees
585  IMPLICIT NONE
586
587  real, intent(in) :: sig, rm, x
588  real ::  cste, logn
589
590  cste = 1.D0 / (sqrt(2.D0*PI)*log(sig))
591  logn = cste / x * exp(-0.5D0*((log(x)-log(rm))/log(sig))**2)
592
593  RETURN
594END FUNCTION logn
595
596
597!*****************************************************************************
598FUNCTION moment(k, ntot, rbar, sigma)
599
600  !*     Computes the kth order moment of the log-normale with parameters : ntot,rbar and sigma
601  !*
602  !*     INPUTS
603  !*            k     : the order of the moment to compute
604  !*            ntot  : the total number of particles of the distribution (ie: M0)
605  !*            rbar  : the mean raidus of the ditribution
606  !*            sigma : the standard deviation of the distribution.
607
608  IMPLICIT NONE
609
610  real, intent(in) :: ntot, rbar, sigma, k
611  real :: moment
612
613  moment = ntot * rbar ** k * exp(0.5D0 *(k*log(sigma))**2)
614
615  RETURN
616END FUNCTION moment
617
618
619!*****************************************************************************
620FUNCTION alpha_k(k,sig)
621
622  IMPLICIT NONE
623
624  real :: k, alpha_k
625  real :: sig
626
627  alpha_k = EXP(0.5 * k*k *(LOG(sig))**2)
628
629  return
630end FUNCTION alpha_k
Note: See TracBrowser for help on using the repository browser.