source: dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/lsc_scav.F90 @ 3814

Last change on this file since 3814 was 3814, checked in by ymipsl, 10 years ago

remove all dynamic dependency in LMDZ physics except for the include "dimensions.h"

YM
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_phy,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.