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

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

Reorganize geometry and grid modules. Prepare physics for unstructutured grid support. Simplify initialization of physics from dynamic.
Compiled only with dynd3dmem, but not tested for moment.

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