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

Last change on this file since 2347 was 2346, checked in by Ehouarn Millour, 9 years ago

Physics/dynamics separation:

  • remove all references to dimensions.h from physics. nbp_lon (==iim) , nbp_lat (==jjm+1) and nbp_lev (==llm) from mod_grid_phy_lmdz should be used instead.
  • added module regular_lonlat_mod in phy_common to store information about the global (lon-lat) grid cell boundaries and centers.

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