source: trunk/LMDZ.VENUS/libf/phyvenus/new_cloud_sedim.F @ 1521

Last change on this file since 1521 was 1442, checked in by slebonnois, 10 years ago

SL: update of the Venus GCM, + corrections on routines used for newstart/start2archive for Titan and Venus, + some modifications on tools
  • Property svn:executable set to *
File size: 11.7 KB
Line 
1      SUBROUTINE new_cloud_sedim(n_lon,n_lev,ptimestep,
2     &                pmidlay,pbndlay,
3     &                pt,
4     &                pq, pdqsed,pdqs_sed,nq,F_sed)
5
6      USE ioipsl
7      USE infotrac
8      USE control_mod
9      USE dimphy
10      USE comgeomphy
11      USE chemparam_mod
12      IMPLICIT NONE
13
14c=======================================================================
15c     
16c=======================================================================
17
18c-----------------------------------------------------------------------
19c   declarations:
20c   -------------
21#include "YOMCST.h"     
22c#include "dimensions.h"
23c#include "dimphys.h"
24c#include "comcstfi.h"
25c#include "tracer.h"
26c#include "callkeys.h"
27       
28c
29c   arguments:
30c   ----------
31
32      INTEGER n_lon                 ! number of horizontal grid points
33      INTEGER n_lev                 ! number of atmospheric layers
34      REAL ptimestep                ! physics time step (s)
35      REAL pmidlay(n_lon,n_lev)     ! pressure at middle layers (Pa)
36      REAL pt(n_lon,n_lev)          ! temperature at mid-layer (l)
37      REAL pbndlay(n_lon,n_lev+1)   ! pressure at layer boundaries
38
39c    Traceurs :
40      integer nq                    ! number of tracers
41      real pq(n_lon,n_lev,nq)       ! tracers (kg/kg)
42c      real pdqfi(n_lon,n_lev,nq)   ! tendency before sedimentation (kg/kg.s-1)
43      real pdqsed(n_lon,n_lev,2)    ! tendency due to sedimentation (kg/kg)
44      real pdqs_sed(n_lon)          ! surface density (Flux if /ptimestep) at surface due to sedimentation (kg.m-2)
45     
46c   local:
47c   ------
48      integer imode
49      integer ig
50      integer iq
51      integer l
52
53      real zlev(n_lon,n_lev+1)      ! altitude at layer boundaries
54      real zlay(n_lon,n_lev)        ! altitude at the midlle layer
55      real zqi_wv(n_lon,n_lev)      ! to locally store H2O tracer
56      real zqi_sa(n_lon,n_lev)      ! to locally store H2SO4 tracer
57      real m_lay (n_lon,n_lev)      ! Layer Pressure over gravity (Dp/g == kg.m-2)
58      real wq(n_lon,n_lev+1)        ! displaced tracer mass (kg.m-2)
59
60c    Physical constant
61c    ~~~~~~~~~~~~~~~~~
62c     Gas molecular viscosity (N.s.m-2)
63c      real,parameter :: visc=1.e-5       ! CO2
64        REAL :: VISCOSITY_CO2
65c     Effective gas molecular radius (m)
66      real,parameter :: molrad=2.2e-10   ! CO2
67     
68c     Cloud density (kg.m-3)
69c     ~~~~~~~~~~~~~~~~~~~~~~
70c      real, DIMENSION(n_lon,n_lev) ::  rho_droplet
71
72      REAL, DIMENSION(n_lon,n_lev+1) ::
73     + wgt_SA                         ! Fraction of H2SO4 in droplet local
74
75c     Stokes speed and sedimentation flux variable
76c     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
77
78      REAL :: A1,A2,A3,A4,            ! coeff du DL du Flux de sedimentation
79     + D_stokes,                      ! coeff de la vitesse de Stokes
80     + Rp_DL,                         ! "Point" du DL
81     + l_mean,                        ! libre parcours moyen (m)
82     + a,b_exp,c                      ! coeff du calcul du Flux de sedimentation
83      REAL, DIMENSION(n_lon,n_lev+1) ::
84     + F_sed                          ! Flux de sedimentation (kg.m-2.s-1 puis en output kg.m-2)
85     
86     
87      REAL :: R_mode0                 ! Rayon mode 0 (m), rayon le plus frequent
88
89
90
91!      PRINT*,'RHO_DROPLET new_cloud_sedim.F'
92!      PRINT*,'rho_droplet',rho_droplet(16,21)
93!      PRINT*,'T',pt(16,21),'WSA',WH2SO4(16,21)
94
95c-----------------------------------------------------------------------
96c    1. Initialization
97c    -----------------
98   
99c     Updating the droplet mass mixing ratio with the partition H2O/H2SO4
100c     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
101
102      do l=1,n_lev
103         do ig=1,n_lon
104         zqi_wv(ig,l) = pq(ig,l,i_h2oliq)
105         zqi_sa(ig,l) = pq(ig,l,i_h2so4liq)
106         wgt_SA(ig,l) = WH2SO4(ig,l)
107         enddo
108      enddo
109
110c     Init F_sed
111      F_sed(:,:) = 0.0E+0
112
113c     Au niveau top+1 , tout égal a 0     
114      wgt_SA(:,n_lev+1) = 0.0E+0   
115
116c    Computing the different layer properties
117c    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
118c    m_lay (kg.m-2)
119c    Ici g=8.87, conflit pour g entre #include "YOMCST.h"
120c       et #include "comcstfi.h"
121
122      do  l=1,n_lev
123         do ig=1, n_lon
124         m_lay(ig,l)=(pbndlay(ig,l) - pbndlay(ig,l+1)) /8.87E+0
125            IF (m_lay(ig,l).LE.0.0) THEN
126            PRINT*,'!!!! STOP PROBLEME SEDIMENTATION!!!!'
127            PRINT*,'!!!!          m_lay <= 0        !!!!'
128            PRINT*,'!!!! STOP PROBLEME SEDIMENTATION!!!!'
129            ENDIF
130         end do
131      end do
132       
133c         Computing sedimentation for droplet "tracer"
134c         ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
135c           pbndlay(:,51)=0 (en parallèle c'est sûr), ne pas l'utiliser pour Fse
136       
137        DO imode=1, nbr_mode
138         DO l = cloudmin, cloudmax
139            DO ig=1,n_lon
140
141c     RD=1000.*RNAVO*RKBOL/RMD avec RMD=43.44 Masse molaire atm venus en g.mol-1     
142           D_stokes=((rho_droplet(ig,l)-pmidlay(ig,l)/(RD*pt(ig,l))))
143     &  *(2./9.)*(RG/VISCOSITY_CO2(pt(ig,l)))
144     
145           l_mean=(pt(ig,l)/pmidlay(ig,l))*
146     &  (0.707*R/(4.*RPI* molrad*molrad * RNAVO))
147     
148           R_mode0=R_MEDIAN(ig,l,imode)*
149     &     EXP(-LOG(STDDEV(ig,l,imode))**2.)
150              IF ((l_mean/(R_mode0)).GT.10.) THEN
151              Rp_DL=R_MEDIAN(ig,l,imode)*
152     &        EXP(3.*LOG(STDDEV(ig,l,imode))**2.)
153              ELSE
154              Rp_DL=R_MEDIAN(ig,l,imode)*
155     &        EXP(4.*LOG(STDDEV(ig,l,imode))**2.)
156              ENDIF
157               
158           a=1.246*l_mean
159       
160           c=0.87/l_mean
161       
162           b_exp=0.42*l_mean*EXP(-c*Rp_DL)
163       
164           A1=a+b_exp*(1.+c*Rp_DL
165     &  +0.5*(Rp_DL*c)**2
166     &  +1./6.*(Rp_DL*c)**3)
167     
168           A2=1.-b_exp*(c
169     &  +Rp_DL*c**2
170     &  +0.5*(Rp_DL**2)*(c**3))
171       
172           A3=0.5*b_exp*(c**2+Rp_DL*c**3)
173       
174           A4=-b_exp*1./6.*c**3
175
176c     Addition des Flux de tous les modes presents     
177       F_sed(ig,l)=F_sed(ig,l)+(rho_droplet(ig,l)*4./3.*RPI*
178     &  NBRTOT(ig,l,imode)*1.0E6*D_stokes*(
179     &  A1*R_MEDIAN(ig,l,imode)**4
180     &  *EXP(8.0*LOG(STDDEV(ig,l,imode))**2.)
181     &  +A2*R_MEDIAN(ig,l,imode)**5
182     &  *EXP(12.5*LOG(STDDEV(ig,l,imode))**2.)
183     &  +A3*R_MEDIAN(ig,l,imode)**6
184     &  *EXP(18.0*LOG(STDDEV(ig,l,imode))**2.)
185     &  +A4*R_MEDIAN(ig,l,imode)**7
186     &  *EXP(24.5*LOG(STDDEV(ig,l,imode))**2.)))
187     
188c      PRINT*,' APRES dTime: F_sed=',F_sed(ig,l), ig, l
189     
190        IF (F_sed(ig,l).GT.m_lay(ig,l)) THEN
191        PRINT*,'==============================================='
192        PRINT*,'WARNING On a epuise la couche', ig, l
193        PRINT*,'On epuise pas une couche avec une espèce
194     &   minoritaire, c est pas bien maaaaaal'
195            PRINT*,'Water',zqi_wv(ig,l),'Sulfuric Acid',zqi_sa(ig,l)
196        PRINT*,'F_sed:',F_sed(ig,l),'m_lay:',m_lay(ig,l)
197        PRINT*,'F_sed/dtphy',F_sed(ig,l)/ptimestep
198        PRINT*,'Pbnd top',pbndlay(ig,l+1),'Temp',pt(ig,l),'Rho',
199     &  rho_droplet(ig,l)
200                PRINT*,'Ntot',NBRTOT(ig,l,:)
201                PRINT*,'StdDev',STDDEV(ig,l,:),'Rmed',R_MEDIAN(ig,l,:)
202                PRINT*,'K_MASS',K_MASS(ig,l,:)
203                PRINT*,'WSA',WH2SO4(ig,l),'RHO',rho_droplet(ig,l)
204               
205c               ELSE
206c               
207c               PRINT*,'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
208c       PRINT*,'WARNING On a PAS epuise la couche', ig, l
209c       PRINT*,'F_sed:',F_sed(ig,l),'m_lay:',m_lay(ig,l)
210c       PRINT*,'F_sed/dtphy',F_sed(ig,l)/ptimestep
211c       PRINT*,'Pbnd top',pbndlay(ig,l+1),'Temp',pt(ig,l),'Rho',
212c     &         rho_droplet(ig,l)(ig,l)
213c               PRINT*,'Ntot',NBRTOT(ig,l),'Ntot m3',NBRTOT(ig,l)*1.0e6
214c               PRINT*,'StdDev',STDDEV(ig,l),'Rmed',R_MEDIAN(ig,l)
215            STOP               
216              ENDIF
217       
218           IF (F_sed(ig,l).LT.0.0d0) THEN
219              PRINT*,"F_sed est négatif !!!"
220              PRINT*,'F_sed:',F_sed(ig,l),'m_lay:',m_lay(ig,l)
221        PRINT*,'F_sed/dtphy',F_sed(ig,l)/ptimestep
222        PRINT*,'Pbnd top',pbndlay(ig,l+1),'Pmid',pmidlay(ig,l)
223        PRINT*,'Temp',pt(ig,l),'Rho',
224     &  rho_droplet(ig,l)
225                PRINT*,'Ntot',NBRTOT(ig,l,imode),'Ntot m3',
226     &  NBRTOT(ig,l,imode)*1.0e6
227                PRINT*,'StdDev',STDDEV(ig,l,imode),'Rmed',
228     &          R_MEDIAN(ig,l,imode)
229                PRINT*,'A1',A1,'A2',A2
230                PRINT*,'A3',A1,'A4',A2
231                PRINT*,'D_stokes',D_stokes
232                STOP
233           ENDIF
234           
235              ENDDO
236             
237c           ELSE           
238c           F_sed(:,l)=0.0d0           
239c           ENDIF
240           
241         ENDDO
242      ENDDO
243
244c     Passage du Flux au Flux pour un pas de temps (== kg.m-2)     
245      F_sed(:,:)=F_sed(:,:)*ptimestep
246
247
248c       VENUS: le flux à la surface est fixé à 0
249c     les conditions P/T en surface ne permettent pas la condensation
250      DO ig=1,n_lon
251      pdqs_sed(ig) = 0.0d0
252      ENDDO
253     
254c         Compute the final tendency:
255c         ---------------------------
256
257c     Partie H2SO4l
258c     ~~~~~~~~~~~~
259
260      DO l = 1, n_lev
261         DO ig=1,n_lon
262            zqi_sa(ig,l) = zqi_sa(ig,l) + (
263     &                         F_sed(ig,l+1)*wgt_SA(ig,l+1)
264     &                       - F_sed(ig,l)*wgt_SA(ig,l))
265     &                       / m_lay(ig,l)
266c     On peut avoir theoriquement le cas ou on epuise tout le VMR present
267             IF (zqi_sa(ig,l).LT.0.0D0) THEN
268               PRINT*,'STOP sedimentation on epuise tout le VMR present'
269               PRINT*,'couche',ig,'level',l
270c               STOP
271c              Ce n est pas juste mais il faudrait alors adapter les pas
272c              de tps de la phys, microphys et chimie
273c              car dans ce cas, c est comme si on epuisait la couche pour un pdtphys
274c              mais en fait on l epuise pour un pdt<pdtphys
275               zqi_sa(ig,l) = 0.0D0
276             ENDIF
277            pdqsed(ig,l,1) = zqi_sa(ig,l) - pq(ig,l,i_h2so4liq)                       
278         ENDDO
279      ENDDO
280
281c     Partie H2Ol
282c     ~~~~~~~~~~~
283                     
284      DO l = 1, n_lev
285         DO ig=1,n_lon
286            zqi_wv(ig,l) = zqi_wv(ig,l) + (
287     &                         F_sed(ig,l+1)*(1. - wgt_SA(ig,l+1))
288     &                       - F_sed(ig,l)*(1. - wgt_SA(ig,l)))
289     &                       / m_lay(ig,l)
290c     On peut avoir theoriquement le cas ou on epuise tout le VMR present
291             IF (zqi_wv(ig,l).LT.0.0D0) THEN
292               PRINT*,'STOP sedimentation on epuise tout le VMR present'
293               PRINT*,'couche',ig,'level',l
294c               STOP
295c              Ce n est pas juste mais il faudrait alors adapter les pas
296c              de tps de la phys, microphys et chimie
297c              car dans ce cas, c est comme si on epuisait la couche pour un pdtphys
298c              mais en fait on l epuise pour un pdt<pdtphys
299               zqi_wv(ig,l) = 0.0D0
300             ENDIF
301            pdqsed(ig,l,2) = zqi_wv(ig,l) - pq(ig,l,i_h2oliq)                   
302         ENDDO
303      ENDDO
304
305c               Save output file in 1D model
306c               ============================
307               
308c      IF (n_lon .EQ. 1) THEN
309c      PRINT*,'Save output sedim'       
310c      DO l = 1, n_lev
311c       DO ig=1,n_lon
312c       WRITE(77,"(i4,','11(e15.8,','))") l,pdqsed(ig,l),zqi(ig,l),
313c     &         (WH2SO4(ig,l)*pq(ig,l,i_h2so4liq)+
314c     &         (1.-WH2SO4(ig,l))*pq(ig,l,i_h2oliq)),
315c     &         pq(ig,l,i_h2so4liq),pq(ig,l,i_h2oliq)
316c      ENDDO
317c       ENDDO
318c      ENDIF   
319
320      RETURN
321      END
322
323*******************************************************************************
324      REAL FUNCTION VISCOSITY_CO2(temp)
325c       Aurélien Stolzenbach 2015
326c       Calcul de la viscosité dynamique du CO2 80°K -> 300°K
327c       Viscosité dynamique en Pa.s
328c       Source: Johnston & Grilly (1942)
329
330c       température en °K
331        REAL, INTENT(IN) :: temp
332       
333        REAL :: denom, numer
334       
335c       Calcul de la viscosité dynamique grâce à la formule de Jones (Lennard-Jones (1924))
336       
337        numer = 200.**(2.27/4.27)-0.435
338        denom = temp**(2.27/4.27)-0.435
339       
340        VISCOSITY_CO2 = (numer/denom)*1015.*(temp/200.)**(3./2.)
341
342c       convertion de Poises*1e7 -> Pa.s       
343        VISCOSITY_CO2 = VISCOSITY_CO2*1.e-8     
344
345      END FUNCTION VISCOSITY_CO2
346*******************************************************************************
347
348
Note: See TracBrowser for help on using the repository browser.