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

Last change on this file since 1661 was 1639, checked in by slebonnois, 8 years ago

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