source: LMDZ6/branches/LMDZ_DECOUPLE/libf/phylmd/lsc_scav.F90 @ 5321

Last change on this file since 5321 was 4815, checked in by nfevrier, 9 months ago

Implémentation de la correction de N. Lebas sur lsc_scav
  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
File size: 11.1 KB
Line 
1!$Id $
2
3SUBROUTINE lsc_scav(pdtime,it,iflag_lscav,  &
4!jyg<
5                    aerosol,  &
6!>jyg
7                    oliq,flxr,flxs,rneb,beta_fisrt,  &
8                    beta_v1,pplay,paprs,t,tr_seri,d_tr_insc,          &
9                    d_tr_bcscav,d_tr_evap,qPrls,beta_st,j)
10  USE ioipsl
11  USE dimphy
12  USE mod_grid_phy_lmdz
13  USE mod_phys_lmdz_para
14  USE traclmdz_mod
15  USE infotrac_phy,ONLY : nbtr,solsym
16  USE iophy
17  IMPLICIT NONE
18!=====================================================================
19! Objet : depot humide (lessivage et evaporation) de traceurs
20! Inspired by routines of Olivier Boucher (mars 1998)
21! author R. Pilon 10 octobre 2012
22! last modification 16/01/2013 (reformulation partie evaporation)
23!=====================================================================
24
25  include "chem.h"
26  include "YOMCST.h"
27  include "YOECUMF.h"
28
29  REAL,INTENT(IN)                        :: pdtime ! time step (s)
30  INTEGER,INTENT(IN)                     :: it     ! tracer number
31  INTEGER,INTENT(IN)                     :: iflag_lscav ! LS scavenging param:
32!                                             3=Reddy_Boucher2004, 4=3+RPilon.
33  REAL,DIMENSION(klon,klev+1),INTENT(IN) :: flxr     ! flux precipitant de pluie
34  REAL,DIMENSION(klon,klev+1),INTENT(IN) :: flxs     ! flux precipitant de neige
35  REAL,INTENT(IN)                        :: oliq ! contenu en eau liquide dans le nuage (kg/kg)
36  REAL,DIMENSION(klon,klev),INTENT(IN)   :: rneb
37  REAL,DIMENSION(klon,klev),INTENT(IN)   :: pplay    ! pression
38  REAL,DIMENSION(klon,klev+1),INTENT(IN) :: paprs    ! pression
39  REAL,DIMENSION(klon,klev),INTENT(IN)   :: t        ! temperature
40! tracers
41  LOGICAL,DIMENSION(nbtr), INTENT(IN)         :: aerosol
42  REAL,DIMENSION(klon,klev,nbtr),INTENT(IN)   :: tr_seri        ! q de traceur 
43  REAL,DIMENSION(klon,klev),INTENT(IN)        :: beta_fisrt     ! taux de conversion de l'eau cond
44  REAL,DIMENSION(klon,klev),INTENT(OUT)       :: beta_v1        ! -- (originale version)
45  REAL,DIMENSION(klon)                        :: his_dh         ! tendance de traceur integre verticalement
46  REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT)  :: d_tr_insc      ! tendance du traceur
47  REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT)  :: d_tr_bcscav  ! tendance de traceur
48  REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT)  :: d_tr_evap
49  ! Nicolas Février
50  REAL,DIMENSION(klon,klev,12),INTENT(OUT)       :: beta_st
51  INTEGER,INTENT(IN)                             :: j           ! indice parmi les 12 traceurs de WETDEP
52
53  REAL,DIMENSION(klon,nbtr),INTENT(OUT)       :: qPrls      !jyg: concentration tra dans pluie LS a la surf.
54  REAL :: dxin,dxev                              ! tendance temporaire de traceur incloud
55  REAL,DIMENSION(klon,klev) :: dxbc       ! tendance temporaire de traceur bc
56
57
58!  variables locales     
59 LOGICAL,SAVE :: debut=.true.
60!$OMP THREADPRIVATE(debut)
61!
62  REAL,PARAMETER :: henry=1.4  ! constante de Henry en mol/l/atm ~1.4 for gases
63  REAL           :: henry_t    !  constante de Henry a T t  (mol/l/atm)
64  REAL,PARAMETER :: kk=2900.   ! coefficient de dependence en T (K)
65  REAL :: f_a     !  rapport de la phase aqueuse a la phase gazeuse
66  REAL :: beta    !  taux de conversion de l'eau en pluie
67
68  INTEGER :: i, k
69  REAL,DIMENSION(klon,klev)    :: scav  !  water liquid content / fraction aqueuse du constituant
70  REAL,DIMENSION(klon,klev)    :: zrho
71  REAL,DIMENSION(klon,klev)    :: zdz
72  REAL,DIMENSION(klon,klev)    :: zmass ! layer mass
73
74  REAL           :: frac_ev       ! cste pour la reevaporation : dropplet shrinking
75!  frac_ev = frac_gas ou frac_aer
76  REAL,PARAMETER :: frac_gas=1.0  ! cste pour la reevaporation pour les gaz
77
78! Nicolas Février
79! Nicolas Lebas : Bugfix in large scale scavenging: variable frac_aer (use for aerosols) was only on the first step and not stored for next timesteps.
80
81  REAL,SAVE                 :: frac_aer      ! cste pour la reevaporation pour les particules
82  REAL,DIMENSION(klon,klev) :: deltaP     ! P(i+1)-P(i)
83  REAL,DIMENSION(klon,klev) :: beta_ev    !  dP/P(i+1)
84
85!$OMP THREADPRIVATE(frac_aer)
86
87!  101.325  m3/l x Pa/atm
88!  R        Pa.m3/mol/K
89!   cste de dissolution pour le depot humide
90  REAL,SAVE :: frac_fine_scav
91  REAL,SAVE :: frac_coar_scav
92!$OMP THREADPRIVATE(frac_fine_scav, frac_coar_scav)
93
94! below-cloud scav variables
95! aerosol : alpha_r=0.001, gas 0.001  (Pruppacher & Klett 1967)
96  REAL,SAVE :: alpha_r  !  coefficient d'impaction pour la pluie
97  REAL,SAVE :: alpha_s  !  coefficient d'impaction pour la neige 
98  REAL,SAVE :: R_r      !  mean raindrop radius (m)
99  REAL,SAVE :: R_s      !  mean snow crystal radius (m)
100!$OMP THREADPRIVATE(alpha_r, alpha_s, R_r, R_s)
101  REAL           :: pr, ps, ice, water
102  real :: conserv
103!
104!!!!!!!!!!!!!!!!!!!! choix lessivage !!!!!!!!!!!!!!!!!!!!!!!!
105!!  logical,save  :: inscav_fisrt
106!!! $OMP THREADPRIVATE(inscav_first)
107!
108!!!!!!!!!!!!!!!!!!!!!!!!!!!
109  IF (debut) THEN
110!
111!  inscav_fisrt=.true.
112!  call getin('inscav_fisrt',inscav_fisrt)
113!  if(inscav_fisrt) then
114!   print*,'beta from fisrtilp.F90, beta = (z_cond - z_oliq)/z_cond, inscav_fisrt=',inscav_fisrt
115!  else
116!   print*,'beta from Reddy and Bocuher 2004 (original version), inscav_fisrt=',inscav_fisrt
117!  endif
118!
119      alpha_r=0.001        !  coefficient d'impaction pour la pluie
120      alpha_s=0.01         !  coefficient d'impaction pour la neige 
121      R_r=0.001            !  mean raindrop radius (m)
122      R_s=0.001            !  mean snow crystal radius (m)
123      !frac_fine_scav=0.7
124      !frac_coar_scav=0.7
125!     frac_aer=0.5 ~ droplet size shrinks by evap
126      frac_aer=0.5
127!
128      OPEN(99,file='lsc_scav_param.data',status='old', &
129                form='formatted',err=9999)
130      READ(99,*,end=9998)  alpha_r
131      READ(99,*,end=9998)  alpha_s
132      READ(99,*,end=9998)  R_r
133      READ(99,*,end=9998)  R_s
134      READ(99,*,end=9998)  frac_fine_scav
135      READ(99,*,end=9998)  frac_coar_scav
136      READ(99,*,end=9998)  frac_aer
1379998  Continue
138      CLOSE(99)
1399999  Continue
140
141   print*,'alpha_r',alpha_r
142   print*,'alpha_s',alpha_s
143   print*,'R_r',R_r
144   print*,'R_s',R_s
145   !print*,'frac_fine_scav',frac_fine_scav
146   !print*,'frac_coar_scav',frac_coar_scav
147   print*,'frac_aer ev',frac_aer
148
149!
150  ENDIF !(debut)
151!!!!!!!!!!!!!!!!!!!!!!!!!!!
152!
153! initialization
154
155  ! Nicolas Février
156  ! Différentiation BC soluble/insoluble
157
158  IF ((TRIM(solsym(it))=='AIBCM') .OR. (TRIM(solsym(it))=='AIPOMM') .OR. (TRIM(solsym(it))=='CIDUSTM')) THEN
159     frac_fine_scav=0.2
160     frac_coar_scav=0.2
161     !PRINT *,'insoluble okidoki'
162  ELSE
163     frac_fine_scav=0.7
164     frac_coar_scav=0.7
165  ENDIF
166
167  dxin=0.
168  dxev=0.
169  beta_ev=0.
170
171  DO i=1,klon
172   his_dh(i)=0.
173  ENDDO
174
175  DO k=1,klev
176   DO i=1, klon
177    dxbc(i,k)=0.
178    beta_v1(i,k)=0.
179    deltaP(i,k)=0.
180   ENDDO
181  ENDDO
182
183  DO k=1,klev
184    DO i=1, klon
185     d_tr_insc(i,k,it)=0.
186     d_tr_bcscav(i,k,it)=0.
187     d_tr_evap(i,k,it)=0.
188    ENDDO
189  ENDDO
190
191!  pressure and size of the layer
192  DO k=klev-1, 1, -1
193   DO i=1, klon
194     zrho(i,k)=pplay(i,k)/t(i,k)/RD   
195     zdz(i,k)=(paprs(i,k)-paprs(i,k+1))/zrho(i,k)/RG
196     zmass(i,k)=(paprs(i,k)-paprs(i,k+1))/RG
197   ENDDO
198  ENDDO
199
200!jyg<
201!!    IF (it.gt.1) THEN                               !  aerosol   
202!! Temporary correction: all non-aerosol tracers are dealt with in the same way.
203!! Should be updated once it has been decided how gases should be dealt with.
204    IF (aerosol(it)) THEN
205!>jyg
206      frac_ev=frac_aer
207    ELSE                                                !  gas
208      frac_ev=frac_gas
209    ENDIF
210
211!jyg<
212!!    IF (it.gt.1) THEN                               !  aerosol   
213    IF (aerosol(it)) THEN
214!>jyg
215     DO k=1, klev
216      DO i=1, klon
217       scav(i,k)=frac_fine_scav
218      ENDDO
219     ENDDO
220    ELSE                  ! gas
221     DO k=1, klev
222      DO i=1, klon
223       henry_t=henry*exp(-kk*(1./298.-1./t(i,k)))    !  mol/l/atm
224       f_a=henry_t/101.325*R*t(i,k)*oliq*zrho(i,k)/rho_water
225       scav(i,k)=f_a/(1.+f_a)
226      ENDDO
227     ENDDO
228    ENDIF
229
230   DO k=klev-1, 1, -1
231    DO i=1, klon
232!  incloud scavenging
233!   if(inscav_fisrt) then
234   if (iflag_lscav .eq. 4) then
235      beta=beta_fisrt(i,k)*rneb(i,k)
236   else
237      beta=flxr(i,k)-flxr(i,k+1)+flxs(i,k)-flxs(i,k+1)
238!      beta=beta/zdz(i,k)/oliq/zrho(i,k)
239      beta=beta/zmass(i,k)/oliq
240      beta=MAX(0.,beta)
241   endif ! (iflag_lscav .eq. 4)
242   beta_v1(i,k)=beta    !! for output
243!
244      dxin=tr_seri(i,k,it)*(exp(-scav(i,k)*beta*pdtime)-1.)
245!      his_dh(i)=his_dh(i)-dxin*zrho(i,k)*zdz(i,k)/pdtime !  kg/m2/s
246      his_dh(i)=his_dh(i)-dxin*zmass(i,k)/pdtime !  kg/m2/s
247      d_tr_insc(i,k,it)=dxin
248      ! Nicolas Février
249      beta_st(i,k,j)=beta_st(i,k,j)-(exp(-scav(i,k)*beta*pdtime)-1.)/pdtime
250
251!  below-cloud impaction
252!jyg<
253!!    IF (it.eq.1) THEN
254    IF (.NOT.aerosol(it)) THEN
255!>jyg
256      d_tr_bcscav(i,k,it)=0.
257    ELSE
258     pr=0.5*(flxr(i,k)+flxr(i,k+1))
259     ps=0.5*(flxs(i,k)+flxs(i,k+1))
260     water=pr*alpha_r/R_r/rho_water
261     ice=ps*alpha_s/R_s/rho_ice
262     dxbc(i,k)=-3./4.*tr_seri(i,k,it)*pdtime*(water+ice)
263!   add tracers from below cloud scav in his_dh
264     his_dh(i)=his_dh(i)-dxbc(i,k)*zmass(i,k)/pdtime !  kg/m2/s
265     d_tr_bcscav(i,k,it)=dxbc(i,k)
266     ! Nicolas Février
267     beta_st(i,k,j)=beta_st(i,k,j)+3./4.*(water+ice)
268    ENDIF
269
270!  reevaporation
271      deltaP(i,k)=flxr(i,k+1)+flxs(i,k+1)-flxr(i,k)-flxs(i,k)
272      deltaP(i,k)=max(deltaP(i,k),0.)
273
274      if(flxr(i,k+1)+flxs(i,k+1).gt.1.e-16) then
275       beta_ev(i,k)=deltaP(i,k)/(flxr(i,k+1)+flxs(i,k+1))
276      else
277       beta_ev(i,k)=0.
278      endif
279
280      beta_ev(i,k)=max(min(1.,beta_ev(i,k)),0.)
281
282!jyg
283     
284      if(abs(1-(1-frac_ev)*beta_ev(i,k)).gt.1.e-16) then
285! remove tracers from precipitation owing to release by evaporation in his_dh
286!      dxev=frac_ev*beta_ev(i,k)*his_dh(i) *pdtime/(zrho(i,k)*zdz(i,k)) &
287      dxev=frac_ev*beta_ev(i,k)*his_dh(i) *pdtime/(zmass(i,k)) &
288                                      /(1 -(1-frac_ev)*beta_ev(i,k))
289       his_dh(i)=his_dh(i)*(1 - frac_ev*beta_ev(i,k) / (1 -(1-frac_ev)*beta_ev(i,k)))
290      else
291!       dxev=his_dh(i) *pdtime/(zrho(i,k)*zdz(i,k))
292       dxev=his_dh(i) *pdtime/(zmass(i,k))
293       his_dh(i)=0.
294      endif
295!      print*,  k, 'beta_ev',beta_ev
296! remove tracers from precipitation owing to release by evaporation in his_dh
297!!      dxev=frac_ev*deltaP(i,k)*pdtime * his_dh(i) /(zrho(i,k)*zdz(i,k))
298!rplmd
299!!      dxev=frac_ev*deltaP(i,k)*his_dh(i) *pdtime/(zrho(i,k)*zdz(i,k)) &
300!!                                      /max(flxr(i,k)+flxs(i,k),1.e-16)
301
302
303      d_tr_evap(i,k,it)=dxev
304!!     tendency is further added in phytrac x = x + dx
305    ENDDO !!  do i
306   ENDDO  !! do k
307
308! Debug beta_st
309    print *,"BETA LSC_SCAV",beta_st(5000,35,7)
310
311
312!jyg (20130114)
313   DO i = 1,klon
314     qPrls(i,it) = his_dh(i)/max(flxr(i,1)+flxs(i,1),1.e-16)
315   ENDDO
316!jyg end
317
318
319! test de conservation
320      conserv=0.
321!      DO k= klev,1,-1
322!        DO i=1, klon
323!         conserv=conserv+d_tr_insc(i,k,it)*(paprs(i,k)-paprs(i,k+1))/RG &
324!                +d_tr_bcscav(i,k,it)*(paprs(i,k)-paprs(i,k+1))/RG  &
325!                +d_tr_evap(i,k,it)*(paprs(i,k)-paprs(i,k+1))/RG
326!      if(it.eq.3) write(*,'(I2,2X,a,e20.12,2X,a,e20.12,2X,a,e20.12,2X,a,e20.12)'),&
327!      k,'lsc conserv ',conserv,'insc',d_tr_insc(i,k,it),'bc',d_tr_bcscav(i,k,it),'ev',d_tr_evap(i,k,it)
328!       ENDDO
329!     ENDDO
330
331!print*, "TAGADATSOINTSOIN"
332
333END SUBROUTINE lsc_scav
Note: See TracBrowser for help on using the repository browser.