source: LMDZ6/trunk/libf/phylmd/Dust/lsc_scav_spl.f90 @ 5272

Last change on this file since 5272 was 5271, checked in by abarral, 13 months ago

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