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

Last change on this file since 3818 was 3818, checked in by millour, 10 years ago

Some partial cleanup on uses of "dimensions.h" in physics.
At this point 3D gcm compiles and bench seems to run fine :-)
EM

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