source: trunk/LMDZ.VENUS/libf/phyvenus/cloudvenus/new_cloud_venus.F @ 3556

Last change on this file since 3556 was 3323, checked in by flefevre, 8 months ago

1) revised cloud microphysical parameters (this changes the results)
2) the number of modes is passed as an argument to cloud routines (this does not change the results)
3) some cosmetics to chemparam_mod.F90 (this does not change the results)
  • Property svn:executable set to *
File size: 13.7 KB
Line 
1!******************************************************************************
2!*     venus_SAS_composition SUBROUTINE
3!*     modified from
4!*     PROGRAM PSC_MODEL_E
5!*     by A. Määttänen
6!*     subroutine for LMDZ+photochemistry VENUS
7!*     by A. Stolzenbach
8!*
9!*     Input/Output files:
10!*     -------------------
11!*
12!----------------------------------------------------------------------------
13      SUBROUTINE new_cloud_venus(nbr_mode,
14     + nblev, nblon,
15     + TT,PP,
16     + mrt_wv,mrt_sa,
17     + mr_wv,mr_sa)
18
19      USE chemparam_mod
20      IMPLICIT NONE
21
22#include "YOMCST.h"
23
24      integer :: nbr_mode
25      INTEGER, INTENT(IN) :: nblon  ! nombre de points horizontaux
26      INTEGER, INTENT(IN) :: nblev  ! nombre de couches verticales
27
28!----------------------------------------------------------------------------
29!     Ambient air state variables:
30      REAL, INTENT(IN), DIMENSION(nblon,nblev) :: mrt_wv,mrt_sa,
31     +                                            TT,PP
32      REAL, INTENT(INOUT), DIMENSION(nblon,nblev) :: mr_wv,mr_sa
33!----------------------------------------------------------------------------
34      INTEGER :: ilon, ilev, imode
35!----------------------------------------------------------------------------
36!     Thermodynamic functions:
37      REAL :: ROSAS
38!----------------------------------------------------------------------------
39!     Auxilary variables:
40      REAL :: NH2SO4,NH2O
41      REAL :: H2SO4_liq,H2O_liq
42      REAL :: CONCM
43      REAL :: MCONDTOT
44      REAL :: RMODE
45      REAL :: WSAFLAG
46      REAL :: K_SAV
47!----------------------------------------------------------------------------
48!     Ridder's Method variables:
49      REAL :: WVMIN, WVMAX, WVACC
50
51      INTEGER :: NBROOT
52
53      INTEGER :: MAXITE
54      PARAMETER(MAXITE=20)
55
56      INTEGER :: NBRAC
57      PARAMETER(NBRAC=5)
58
59      INTEGER :: FLAG
60     
61      REAL :: qmass
62c       masse volumique du coeur (kg.m-3)
63c     ATTENTION ! DOIT ETRE COHERENT AVEC new_cloud_sedim !
64      REAL, PARAMETER :: rho_core = 2500.0
65!----------------------------------------------------------------------------
66!     External functions needed:
67      REAL :: IRFRMWV
68!----------------------------------------------------------------------------
69
70
71! >>> Program starts here:
72
73!AM Venus
74! These aerosols will then be given an equilibrium composition for the given size distribution
75
76  ! Hanna Vehkamäki and Markku Kulmala and Ismo Napari
77  ! and Kari E. J. Lehtinen and Claudia Timmreck and Madis Noppel and Ari Laaksonen, 2002,
78  ! An improved parameterization for sulfuric acid/water nucleation rates for tropospheric
79  !and stratospheric conditions, () J. Geophys. Res., 107, PP. 4622-4631
80
81!===========================================
82!     Debut boucle sur niveau et lat,lon
83!===========================================
84!     Init, tous les points=0, cela met les niveaux > cloudmax et < cloudmin a 0
85      NBRTOT(:,:,:)=0.0E+0
86      WH2SO4(:,:)=0.0E+0
87      rho_droplet(:,:)=0.0E+0
88
89      DO ilev=cloudmin, cloudmax
90      DO ilon=1, nblon
91
92!       Boucle sur les modes
93        RMODE=0.0E+0
94        K_SAV = 0.0
95
96        DO imode=1, nbr_mode
97          IF (K_MASS(ilon,ilev,imode).GT.K_SAV) THEN
98!       RMODE est le rayon modal de la distribution en volume du mode le plus
99!       representatif pour la Mtot
100            RMODE=R_MEDIAN(ilon,ilev,imode)*
101     &      EXP(2.*(DLOG(STDDEV(ilon,ilev,imode))**2.))
102            K_SAV=K_MASS(ilon,ilev,imode)
103          ENDIF
104        ENDDO ! FIN boucle imode
105
106!       Initialisation des bornes pour WV
107        WVMIN=1.E-90
108        WVMAX=mrt_wv(ilon,ilev)
109
110!       Accuracy de WVeq
111        WVACC=WVMAX*1.0E-3
112
113!       BRACWV borne la fonction f(WV) - WV = 0
114!       de WV=0 à WV=WVtot on cherche l'intervalle où f(WV) - WV = 0
115!       avec précisément f(WVliq de WSA<=WVinput) + WVinput - WVtot = 0
116!       Elle fait appel à la fct/ssrtine ITERWV()
117
118        CALL BRACWV(TT(ilon,ilev),PP(ilon,ilev),WVMIN,WVMAX,NBRAC,
119     &  RMODE,mrt_wv(ilon,ilev),mrt_sa(ilon,ilev),FLAG,WSAFLAG,NBROOT)
120
121        SELECT CASE(FLAG)
122
123        CASE(1)
124!         Cas NROOT=1 ou NROOT>1 mais dans un intervalle restreint WVTOT (cas courant)
125!       IRFRMWV Ridder's method pour trouver, sur [WVmin,WVmax], WVo tel que f(WVo) - WVo = 0
126!       Elle fait appel ˆ la fct/ssrtine ITERWV()
127
128          WH2SO4(ilon,ilev)=IRFRMWV(TT(ilon,ilev),PP(ilon,ilev),
129     &    WVMIN,WVMAX,WVACC,MAXITE,RMODE,
130     &    mrt_wv(ilon,ilev),mrt_sa(ilon,ilev),NBROOT)
131
132          rho_droplet(ilon,ilev)=ROSAS(TT(ilon,ilev),WH2SO4(ilon,ilev))
133
134!          IF (rho_droplet(ilon,ilev).LT.1100.) THEN
135!            PRINT*,'PROBLEM RHO_DROPLET'
136!            PRINT*,'rho_droplet',rho_droplet(ilon,ilev)
137!            PRINT*,'T',TT(ilon,ilev),'WSA',WH2SO4(ilon,ilev)
138!            PRINT*,'RHODROPLET',ROSAS(TT(ilon,ilev),WH2SO4(ilon,ilev))
139!            PRINT*,'FLAG',FLAG,'NROOT',NBROOT
140!            STOP
141!          ENDIF
142
143          CONCM= PP(ilon,ilev)/(1.3806488E-23*TT(ilon,ilev)) !air number density, molec/m3
144
145              NH2SO4=mrt_sa(ilon,ilev)*CONCM
146              NH2O=mrt_wv(ilon,ilev)*CONCM
147
148          CALL CALCM_SAT(NH2SO4,NH2O,WH2SO4(ilon,ilev),
149     &       rho_droplet(ilon,ilev),TT(ilon,ilev),
150     &       H2SO4_liq,H2O_liq,MCONDTOT)
151
152!       Boucle sur les modes
153! *********************************************
154!       AVEC MODE 3 type J. Cimino 1982, Icarus
155! *********************************************
156!       Mode 1 et 2 avec les Kmass modifies
157          DO imode=1, nbr_mode-1
158!           calcul qmass
159              qmass = (rho_core*qrad**3)/
160     &      (rho_core*qrad**3 + rho_droplet(ilon,ilev)*(1.-qrad**3))
161
162            IF (K_MASS(ilon,ilev,imode).GT.0.) THEN
163              NBRTOT(ilon,ilev,imode)= 1.E-6*3./(4.*RPI)*
164     &        K_MASS(ilon,ilev,imode)/(1.-qmass*K_MASS(ilon,ilev,3))*
165     &        MCONDTOT*EXP(-4.5*DLOG(STDDEV(ilon,ilev,imode))**2.)/
166     &        (R_MEDIAN(ilon,ilev,imode)**3.)
167            ELSE
168              NBRTOT(ilon,ilev,imode)=0.0E+0
169            ENDIF
170          ENDDO
171!       Mode 3 reste identique car on veut le N correspondant aux mesures de Knollenberg
172          IF (K_MASS(ilon,ilev,3).GT.0.) THEN
173            NBRTOT(ilon,ilev,3)= 1.E-6*3./(4.*RPI)*
174     &      K_MASS(ilon,ilev,3)*MCONDTOT*
175     &      EXP(-4.5*DLOG(STDDEV(ilon,ilev,3))**2.)/
176     &      (R_MEDIAN(ilon,ilev,3)**3.)
177          ELSE
178              NBRTOT(ilon,ilev,3)=0.0E+0
179          ENDIF
180
181!       Passage de #/m3 en VMR
182          H2O_liq=H2O_liq/CONCM
183          H2SO4_liq=H2SO4_liq/CONCM
184
185          mr_wv(ilon,ilev)=mrt_wv(ilon,ilev)-H2O_liq
186          mr_sa(ilon,ilev)=mrt_sa(ilon,ilev)-H2SO4_liq
187
188!          Problemes quand on a condense tout, on peut obtenir des -1e-24
189!               aprs la soustraction et conversion de ND ˆ VMR
190          IF (mr_wv(ilon,ilev).LE.0.0) mr_wv(ilon,ilev)=1.0E-30
191          IF (mr_sa(ilon,ilev).LE.0.0) mr_sa(ilon,ilev)=1.0E-30
192
193
194
195        CASE(2)
196!       Cas NROOT=0 mais proche de 0
197
198          WH2SO4(ilon,ilev)=WSAFLAG
199
200          rho_droplet(ilon,ilev)=ROSAS(TT(ilon,ilev),WH2SO4(ilon,ilev))
201
202!     ATTENTION ce IF ne sert a rien en fait, juste a retenir une situation
203!     ubuesque dans mon code ou sans ce IF les valeurs de rho_droplets sont
204!     incoherentes avec TT et WH2SO4 (a priori lorsque NTOT=0)
205!     Juste le fait de METTRE un IF fait que rho_droplet a la bonne valeur
206!     donne par ROSAS (cf test externe en Python), sinon, la valeur est trop
207!     basse (de l'ordre de 1000 kg/m3) et correspond parfois a la valeur avec
208!     WSA=0.1 (pas totalement sur)
209!     En tous cas, incoherent avec ce qui est attendue pour le WSA et T donneeŽ
210!     La version avec le IF (rho<1100 & WSA>0.1) est CORRECTE, rho_droplet a
211!     la bonne valeur (tests externes Python confirment)
212
213          IF ((rho_droplet(ilon,ilev).LT.1100.).AND.
214     &      (WH2SO4(ilon,ilev).GT.0.1))THEN
215            PRINT*,'PROBLEM RHO_DROPLET'
216            PRINT*,'rho_droplet',rho_droplet(ilon,ilev)
217            PRINT*,'T',TT(ilon,ilev),'WSA',WH2SO4(ilon,ilev)
218            PRINT*,'RHODROPLET',ROSAS(TT(ilon,ilev),WH2SO4(ilon,ilev))
219            PRINT*,'FLAG',FLAG,'NROOT',NBROOT
220            STOP
221          ENDIF
222
223
224          CONCM= PP(ilon,ilev)/(1.3806488E-23*TT(ilon,ilev)) !air number density, molec/m3
225
226            NH2SO4=mrt_sa(ilon,ilev)*CONCM
227            NH2O=mrt_wv(ilon,ilev)*CONCM
228
229          CALL CALCM_SAT(NH2SO4,NH2O,WH2SO4(ilon,ilev),
230     &       rho_droplet(ilon,ilev),TT(ilon,ilev),
231     &       H2SO4_liq,H2O_liq,MCONDTOT)
232
233!       Boucle sur les modes
234! *********************************************
235!       AVEC MODE 3 type J. Cimino 1982, Icarus
236! *********************************************
237!       Mode 1 et 2 avec alcul coeff*Kmass
238          DO imode=1, nbr_mode-1
239!           calcul qmass
240              qmass = (rho_core*qrad**3)/
241     &      (rho_core*qrad**3 + rho_droplet(ilon,ilev)*(1.-qrad**3))
242
243            IF (K_MASS(ilon,ilev,imode).GT.0.) THEN
244              NBRTOT(ilon,ilev,imode)= 1.E-6*3./(4.*RPI)*
245     &        K_MASS(ilon,ilev,imode)/(1.-qmass*K_MASS(ilon,ilev,3))*
246     &        MCONDTOT*EXP(-4.5*DLOG(STDDEV(ilon,ilev,imode))**2.)/
247     &        (R_MEDIAN(ilon,ilev,imode)**3.)
248            ELSE
249              NBRTOT(ilon,ilev,imode)=0.0E+0
250            ENDIF
251          ENDDO
252!       Mode 3 reste identique car on veut le N correspondant aux mesures de Knollenberg
253          IF (K_MASS(ilon,ilev,3).GT.0.) THEN
254            NBRTOT(ilon,ilev,3)= 1.E-6*3./(4.*RPI)*
255     &      K_MASS(ilon,ilev,3)*MCONDTOT*
256     &      EXP(-4.5*DLOG(STDDEV(ilon,ilev,3))**2.)/
257     &      (R_MEDIAN(ilon,ilev,3)**3.)
258          ELSE
259              NBRTOT(ilon,ilev,3)=0.0E+0
260          ENDIF
261
262
263!       Passage de #/m3 en VMR
264          H2O_liq=H2O_liq/CONCM
265          H2SO4_liq=H2SO4_liq/CONCM
266
267          mr_wv(ilon,ilev)=mrt_wv(ilon,ilev)-H2O_liq
268          mr_sa(ilon,ilev)=mrt_sa(ilon,ilev)-H2SO4_liq
269
270!          Problmes quand on a condense tout, on peut obtenir des -1e-24
271!               aprs la soustraction et conversion de ND ˆ VMR
272          IF (mr_wv(ilon,ilev).LE.0.0) mr_wv(ilon,ilev)=1.0E-30
273          IF (mr_sa(ilon,ilev).LE.0.0) mr_sa(ilon,ilev)=1.0E-30
274
275        CASE(3)
276!         Cas 0 NROOT
277            mr_wv(ilon,ilev)=mrt_wv(ilon,ilev)
278            mr_sa(ilon,ilev)=mrt_sa(ilon,ilev)
279            rho_droplet(ilon,ilev)=0.0E+0
280            WH2SO4(ilon,ilev)=0.0E+0
281            DO imode=1, nbr_mode
282              NBRTOT(ilon,ilev,imode)=0.0E+0
283            ENDDO
284
285        END SELECT
286      ENDDO   !FIN boucle ilon
287      ENDDO   !FIN boucle ilev
288
289      END SUBROUTINE new_cloud_venus
290
291
292!******************************************************************
293      SUBROUTINE CALCM_SAT(H2SO4,H2O,WSA,DENSO4,
294     + T,H2SO4COND,H2OCOND,RMTOT)
295
296!     DERIVE NO (TOTAL NUMBER OF AEROSOL PARTICLES CONCENTRATION)
297!     FROM TOTAL H2SO4 AND RMOD/SIGMA OF AEROSOL LOG-NORMAL
298!                                       SIZE DISTRIBTUION
299!     ASSUMING ALL THE H2SO4 ABOVE MIXTURE SAT PRESSURE modified by H2SO4 activity IS CONDENSED
300!    ---------------------------------------------------------------
301!     INPUT:
302!     H2SO4: #/m3 of total H2SO4
303!       H2O  : #/m3 of total H2O
304!     WSA: aerosol H2SO4 weight fraction (fraction)
305!     DENSO4: aerosol volumic mass (kg/m3 = aerosol mass/aerosol volume)
306!       for total mass, almost same result with ro=1.67 gr/cm3
307!     RSTDEV: standard deviation of aerosol distribution (no unit)
308!     RADIUS: MEDIAN radius (m)
309!     T: temperature (K)
310!
311!     OUTPUT:
312!     RMTOT: Total condensed "Mass" (M_tot_distrib / rho_droplet), sans dimension
313!            mais rho_droplet et M_tot_distrib doivent etre de meme dimension
314!       H2OCOND
315!       H2SO4COND
316
317
318
319      IMPLICIT NONE
320
321      REAL, INTENT(IN) :: H2SO4, H2O, WSA
322      REAL, INTENT(IN) :: DENSO4, T
323      REAL, INTENT(OUT) :: H2OCOND, H2SO4COND, RMTOT
324!     working variables
325      REAL :: RMH2S4
326      REAL :: DND2,pstand,lpar,acidps
327      REAL :: x1, satpacid
328      REAL , DIMENSION(2):: act
329!
330!     masse of an H2SO4 molecule (kg)
331      RMH2S4=98.078/(6.02214129E+26)
332
333      pstand=1.01325E+5 !Pa  1 atm pressure
334
335        x1=(WSA/98.08)/(WSA/98.08 + ((1.-WSA)/18.0153))
336
337        call zeleznik(x1,t,act)
338
339!pure acid satur vapor pressure
340        lpar= -11.695+DLOG(pstand) ! Zeleznik
341        acidps=1/360.15-1.0/t+0.38/545.
342     + *(1.0+DLOG(360.15/t)-360.15/t)
343        acidps = 10156.0*acidps +lpar
344        acidps = DEXP(acidps)    !Pa
345
346!acid sat.vap.PP over mixture (flat surface):
347        satpacid=act(2)*acidps ! Pa
348
349!       Conversion from Pa to N.D #/m3
350        DND2=satpacid/(1.3806488E-23*T)
351!       Conversion from N.D #/m3 TO #/cm3
352!        DND2=DND2*1.d-6
353
354!       H2SO4COND N.D #/m3 condensee ssi H2SO4>H2SO4sat
355        IF (H2SO4.GE.DND2) THEN
356        H2SO4COND=H2SO4-DND2
357!       calcul de H2O cond correspondant a H2SO4 cond
358        H2OCOND=H2SO4COND*98.078*(1.0-WSA)/(18.0153*WSA)
359
360!     RMTOT: = Mass of H2SO4 satur per m3 of air/ Mass of sulfuric acid part of droplet solution per m3
361!     RMTOT = Mtot/rho_droplet
362!       RMTOT=M_distrib/rho_droplet
363
364        RMTOT=H2SO4COND*RMH2S4/(DENSO4*WSA)
365
366!       Si on a H2SO4<H2SO4sat on ne condense rien et NDTOT=0
367        ELSE
368        H2SO4COND=0.0E+0
369        H2OCOND=0.0E+0
370        RMTOT=0.0E+0
371        END IF
372
373!       Test si H2O en defaut H2Ocond>H2O dispo
374        IF ((H2OCOND.GT.H2O).AND.(H2SO4.GE.DND2)) THEN
375
376!     Si H2O en dŽfaut, on as pas le bon WSA!
377!     En effet, normalement, on a exactement le WSA correspondant a
378!     WVg + WVl = WVtot
379!     Dans les cas o WVtot, SAtot sont trs faibles (Upper Haze) ou
380!     quand T est grand (Lower Haze), le modle reprŽsente mal le WSA
381!     cf carte NCL, avec des max erreur absolue de 0.1 sur le WSA
382
383!      PRINT*,'PROBLEM H2O EN DEFAUT'
384!      PRINT*,'H2OCOND',H2OCOND,'H2O',H2O
385!      PRINT*,'WSA',WSA,'RHO',DENSO4
386!      STOP
387
388
389!       On peut alors condenser tout le H2O dispo
390        H2OCOND=H2O
391!       On met alors egalement a jour le H2SO4 cond correspondant au H2O cond
392        H2SO4COND=H2OCOND*18.0153*WSA/(98.078*(1.0-WSA))
393
394!     RMTOT: = Mass of H2SO4 satur per m3 of air/ Mass of sulfuric acid part of droplet solution per m3
395!       RMTOT=Volume of aerosol m3 /m3 of air
396!       Volume of aerosol/m3 air
397
398        RMTOT=H2SO4COND*RMH2S4/(DENSO4*WSA)
399
400      END IF
401
402      END SUBROUTINE CALCM_SAT
403
Note: See TracBrowser for help on using the repository browser.