source: LMDZ5/trunk/libf/phylmd/lsc_scav.F90 @ 2215

Last change on this file since 2215 was 1907, checked in by lguez, 11 years ago

Added a copyright property to every file of the distribution, except
for the fcm files (which have their own copyright). Use svn propget on
a file to see the copyright. For instance:

$ svn propget copyright libf/phylmd/physiq.F90
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

Also added the files defining the CeCILL version 2 license, in French
and English, at the top of the LMDZ tree.
  • 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: 9.7 KB
Line 
1!$Id $
2
3SUBROUTINE lsc_scav(pdtime,it,iflag_lscav,oliq,flxr,flxs,rneb,beta_fisrt,  &
4                    beta_v1,pplay,paprs,t,tr_seri,d_tr_insc,          &
5                    d_tr_bcscav,d_tr_evap,qPrls)
6  USE ioipsl
7  USE dimphy
8  USE mod_grid_phy_lmdz
9  USE mod_phys_lmdz_para
10  USE traclmdz_mod
11  USE infotrac,ONLY : nbtr
12  USE comgeomphy
13  USE iophy
14  IMPLICIT NONE
15!=====================================================================
16! Objet : depot humide (lessivage et evaporation) de traceurs
17! Inspired by routines of Olivier Boucher (mars 1998)
18! author R. Pilon 10 octobre 2012
19! last modification 16/01/2013 (reformulation partie evaporation)
20!=====================================================================
21
22  include "dimensions.h"
23  include "chem.h"
24  include "YOMCST.h"
25  include "YOECUMF.h"
26
27  REAL,INTENT(IN)                        :: pdtime ! time step (s)
28  INTEGER,INTENT(IN)                     :: it     ! tracer number
29  INTEGER,INTENT(IN)                     :: iflag_lscav ! LS scavenging param:
30!                                             3=Reddy_Boucher2004, 4=3+RPilon.
31  REAL,DIMENSION(klon,klev+1),INTENT(IN) :: flxr     ! flux precipitant de pluie
32  REAL,DIMENSION(klon,klev+1),INTENT(IN) :: flxs     ! flux precipitant de neige
33  REAL,INTENT(IN)                        :: oliq ! contenu en eau liquide dans le nuage (kg/kg)
34  REAL,DIMENSION(klon,klev),INTENT(IN)   :: rneb
35  REAL,DIMENSION(klon,klev),INTENT(IN)   :: pplay    ! pression
36  REAL,DIMENSION(klon,klev+1),INTENT(IN) :: paprs    ! pression
37  REAL,DIMENSION(klon,klev),INTENT(IN)   :: t        ! temperature
38! tracers
39  REAL,DIMENSION(klon,klev,nbtr),INTENT(IN)   :: tr_seri        ! q de traceur 
40  REAL,DIMENSION(klon,klev),INTENT(IN)        :: beta_fisrt     ! taux de conversion de l'eau cond
41  REAL,DIMENSION(klon,klev),INTENT(OUT)       :: beta_v1        ! -- (originale version)
42  REAL,DIMENSION(klon)                        :: his_dh         ! tendance de traceur integre verticalement
43  REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT)  :: d_tr_insc      ! tendance du traceur
44  REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT)  :: d_tr_bcscav  ! tendance de traceur
45  REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT)  :: d_tr_evap
46  REAL,DIMENSION(klon,nbtr),INTENT(OUT)       :: qPrls      !jyg: concentration tra dans pluie LS a la surf.
47  REAL :: dxin,dxev                              ! tendance temporaire de traceur incloud
48  REAL,DIMENSION(klon,klev) :: dxbc       ! tendance temporaire de traceur bc
49
50
51!  variables locales     
52 LOGICAL,SAVE :: debut=.true.
53!$OMP THREADPRIVATE(debut)
54!
55  REAL,PARAMETER :: henry=1.4  ! constante de Henry en mol/l/atm ~1.4 for gases
56  REAL           :: henry_t    !  constante de Henry a T t  (mol/l/atm)
57  REAL,PARAMETER :: kk=2900.   ! coefficient de dependence en T (K)
58  REAL :: f_a     !  rapport de la phase aqueuse a la phase gazeuse
59  REAL :: beta    !  taux de conversion de l'eau en pluie
60
61  INTEGER :: i, k
62  REAL,DIMENSION(klon,klev)    :: scav  !  water liquid content / fraction aqueuse du constituant
63  REAL,DIMENSION(klon,klev)    :: zrho
64  REAL,DIMENSION(klon,klev)    :: zdz
65  REAL,DIMENSION(klon,klev)    :: zmass ! layer mass
66
67  REAL           :: frac_ev       ! cste pour la reevaporation : dropplet shrinking
68!  frac_ev = frac_gas ou frac_aer
69  REAL,PARAMETER :: frac_gas=1.0  ! cste pour la reevaporation pour les gaz
70  REAL           :: frac_aer      ! cste pour la reevaporation pour les particules
71  REAL,DIMENSION(klon,klev) :: deltaP     ! P(i+1)-P(i)
72  REAL,DIMENSION(klon,klev) :: beta_ev    !  dP/P(i+1)
73
74!  101.325  m3/l x Pa/atm
75!  R        Pa.m3/mol/K
76!   cste de dissolution pour le depot humide
77  REAL,SAVE :: frac_fine_scav
78  REAL,SAVE :: frac_coar_scav
79!$OMP THREADPRIVATE(frac_fine_scav, frac_coar_scav)
80
81! below-cloud scav variables
82! aerosol : alpha_r=0.001, gas 0.001  (Pruppacher & Klett 1967)
83  REAL,SAVE :: alpha_r  !  coefficient d'impaction pour la pluie
84  REAL,SAVE :: alpha_s  !  coefficient d'impaction pour la neige 
85  REAL,SAVE :: R_r      !  mean raindrop radius (m)
86  REAL,SAVE :: R_s      !  mean snow crystal radius (m)
87!$OMP THREADPRIVATE(alpha_r, alpha_s, R_r, R_s)
88  REAL           :: pr, ps, ice, water
89  real :: conserv
90!
91!!!!!!!!!!!!!!!!!!!! choix lessivage !!!!!!!!!!!!!!!!!!!!!!!!
92!!  logical,save  :: inscav_fisrt
93!!! $OMP THREADPRIVATE(inscav_first)
94!
95!!!!!!!!!!!!!!!!!!!!!!!!!!!
96  IF (debut) THEN
97!
98!  inscav_fisrt=.true.
99!  call getin('inscav_fisrt',inscav_fisrt)
100!  if(inscav_fisrt) then
101!   print*,'beta from fisrtilp.F90, beta = (z_cond - z_oliq)/z_cond, inscav_fisrt=',inscav_fisrt
102!  else
103!   print*,'beta from Reddy and Bocuher 2004 (original version), inscav_fisrt=',inscav_fisrt
104!  endif
105!
106      alpha_r=0.001        !  coefficient d'impaction pour la pluie
107      alpha_s=0.01         !  coefficient d'impaction pour la neige 
108      R_r=0.001            !  mean raindrop radius (m)
109      R_s=0.001            !  mean snow crystal radius (m)
110      frac_fine_scav=0.7
111      frac_coar_scav=0.7
112!     frac_aer=0.5 ~ droplet size shrinks by evap
113      frac_aer=0.5
114!
115      OPEN(99,file='lsc_scav_param.data',status='old', &
116                form='formatted',err=9999)
117      READ(99,*,end=9998)  alpha_r
118      READ(99,*,end=9998)  alpha_s
119      READ(99,*,end=9998)  R_r
120      READ(99,*,end=9998)  R_s
121      READ(99,*,end=9998)  frac_fine_scav
122      READ(99,*,end=9998)  frac_coar_scav
123      READ(99,*,end=9998)  frac_aer
1249998  Continue
125      CLOSE(99)
1269999  Continue
127
128   print*,'alpha_r',alpha_r
129   print*,'alpha_s',alpha_s
130   print*,'R_r',R_r
131   print*,'R_s',R_s
132   print*,'frac_fine_scav',frac_fine_scav
133   print*,'frac_coar_scav',frac_coar_scav
134   print*,'frac_aer ev',frac_aer
135
136!
137  ENDIF !(debut)
138!!!!!!!!!!!!!!!!!!!!!!!!!!!
139!
140! initialization
141  dxin=0.
142  dxev=0.
143  beta_ev=0.
144
145  DO i=1,klon
146   his_dh(i)=0.
147  ENDDO
148
149  DO k=1,klev
150   DO i=1, klon
151    dxbc(i,k)=0.
152    beta_v1(i,k)=0.
153    deltaP(i,k)=0.
154   ENDDO
155  ENDDO
156
157  DO k=1,klev
158    DO i=1, klon
159     d_tr_insc(i,k,it)=0.
160     d_tr_bcscav(i,k,it)=0.
161     d_tr_evap(i,k,it)=0.
162    ENDDO
163  ENDDO
164
165!  pressure and size of the layer
166  DO k=klev-1, 1, -1
167   DO i=1, klon
168     zrho(i,k)=pplay(i,k)/t(i,k)/RD   
169     zdz(i,k)=(paprs(i,k)-paprs(i,k+1))/zrho(i,k)/RG
170     zmass(i,k)=(paprs(i,k)-paprs(i,k+1))/RG
171   ENDDO
172  ENDDO
173
174    IF (it.gt.1) THEN                               !  aerosol   
175      frac_ev=frac_aer
176    ELSE                                                !  gas
177      frac_ev=frac_gas
178    ENDIF
179
180    IF(it.gt.1) then  ! aerosol
181     DO k=1, klev
182      DO i=1, klon
183       scav(i,k)=frac_fine_scav
184      ENDDO
185     ENDDO
186    ELSE                  ! gas
187     DO k=1, klev
188      DO i=1, klon
189       henry_t=henry*exp(-kk*(1./298.-1./t(i,k)))    !  mol/l/atm
190       f_a=henry_t/101.325*R*t(i,k)*oliq*zrho(i,k)/rho_water
191       scav(i,k)=f_a/(1.+f_a)
192      ENDDO
193     ENDDO
194    ENDIF
195
196   DO k=klev-1, 1, -1
197    DO i=1, klon
198!  incloud scavenging
199!   if(inscav_fisrt) then
200   if (iflag_lscav .eq. 4) then
201      beta=beta_fisrt(i,k)*rneb(i,k)
202   else
203      beta=flxr(i,k)-flxr(i,k+1)+flxs(i,k)-flxs(i,k+1)
204!      beta=beta/zdz(i,k)/oliq/zrho(i,k)
205      beta=beta/zmass(i,k)/oliq
206      beta=MAX(0.,beta)
207   endif ! (iflag_lscav .eq. 4)
208   beta_v1(i,k)=beta    !! for output
209!
210      dxin=tr_seri(i,k,it)*(exp(-scav(i,k)*beta*pdtime)-1.)
211!      his_dh(i)=his_dh(i)-dxin*zrho(i,k)*zdz(i,k)/pdtime !  kg/m2/s
212      his_dh(i)=his_dh(i)-dxin*zmass(i,k)/pdtime !  kg/m2/s
213      d_tr_insc(i,k,it)=dxin
214
215!  below-cloud impaction
216    IF(it.eq.1) then
217      d_tr_bcscav(i,k,it)=0.
218    ELSE
219     pr=0.5*(flxr(i,k)+flxr(i,k+1))
220     ps=0.5*(flxs(i,k)+flxs(i,k+1))
221     water=pr*alpha_r/R_r/rho_water
222     ice=ps*alpha_s/R_s/rho_ice
223     dxbc(i,k)=-3./4.*tr_seri(i,k,it)*pdtime*(water+ice)
224!   add tracers from below cloud scav in his_dh
225     his_dh(i)=his_dh(i)-dxbc(i,k)*zmass(i,k)/pdtime !  kg/m2/s
226     d_tr_bcscav(i,k,it)=dxbc(i,k)
227    ENDIF
228
229!  reevaporation
230      deltaP(i,k)=flxr(i,k+1)+flxs(i,k+1)-flxr(i,k)-flxs(i,k)
231      deltaP(i,k)=max(deltaP(i,k),0.)
232
233      if(flxr(i,k+1)+flxs(i,k+1).gt.1.e-16) then
234       beta_ev(i,k)=deltaP(i,k)/(flxr(i,k+1)+flxs(i,k+1))
235      else
236       beta_ev(i,k)=0.
237      endif
238
239      beta_ev(i,k)=max(min(1.,beta_ev(i,k)),0.)
240
241!jyg
242     
243      if(abs(1-(1-frac_ev)*beta_ev(i,k)).gt.1.e-16) then
244! remove tracers from precipitation owing to release by evaporation in his_dh
245!      dxev=frac_ev*beta_ev(i,k)*his_dh(i) *pdtime/(zrho(i,k)*zdz(i,k)) &
246      dxev=frac_ev*beta_ev(i,k)*his_dh(i) *pdtime/(zmass(i,k)) &
247                                      /(1 -(1-frac_ev)*beta_ev(i,k))
248       his_dh(i)=his_dh(i)*(1 - frac_ev*beta_ev(i,k) / (1 -(1-frac_ev)*beta_ev(i,k)))
249      else
250!       dxev=his_dh(i) *pdtime/(zrho(i,k)*zdz(i,k))
251       dxev=his_dh(i) *pdtime/(zmass(i,k))
252       his_dh(i)=0.
253      endif
254!      print*,  k, 'beta_ev',beta_ev
255! remove tracers from precipitation owing to release by evaporation in his_dh
256!!      dxev=frac_ev*deltaP(i,k)*pdtime * his_dh(i) /(zrho(i,k)*zdz(i,k))
257!rplmd
258!!      dxev=frac_ev*deltaP(i,k)*his_dh(i) *pdtime/(zrho(i,k)*zdz(i,k)) &
259!!                                      /max(flxr(i,k)+flxs(i,k),1.e-16)
260
261
262      d_tr_evap(i,k,it)=dxev
263!!     tendency is further added in phytrac x = x + dx
264    ENDDO !!  do i
265   ENDDO  !! do k
266
267!jyg (20130114)
268   DO i = 1,klon
269     qPrls(i,it) = his_dh(i)/max(flxr(i,1)+flxs(i,1),1.e-16)
270   ENDDO
271!jyg end
272
273
274! test de conservation
275      conserv=0.
276!      DO k= klev,1,-1
277!        DO i=1, klon
278!         conserv=conserv+d_tr_insc(i,k,it)*(paprs(i,k)-paprs(i,k+1))/RG &
279!                +d_tr_bcscav(i,k,it)*(paprs(i,k)-paprs(i,k+1))/RG  &
280!                +d_tr_evap(i,k,it)*(paprs(i,k)-paprs(i,k+1))/RG
281!      if(it.eq.3) write(*,'(I2,2X,a,e20.12,2X,a,e20.12,2X,a,e20.12,2X,a,e20.12)'),&
282!      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)
283!       ENDDO
284!     ENDDO
285
286END SUBROUTINE lsc_scav
Note: See TracBrowser for help on using the repository browser.