source: LMDZ6/branches/Amaury_dev/libf/phylmd/Dust/lsc_scav_spl.F90 @ 5103

Last change on this file since 5103 was 5103, checked in by abarral, 8 weeks ago

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