source: LMDZ6/branches/LMDZ-INCA-Dyn/libf/phylmd/Dust/lsc_scav_spl.F90 @ 5420

Last change on this file since 5420 was 3786, checked in by asima, 4 years ago

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