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

Last change on this file since 5182 was 5182, checked in by abarral, 2 months ago

(WIP) Replace REPROBUS CPP KEY by logical
properly name modules
File size: 11.1 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 lmdz_grid_phy
11  USE lmdz_phys_para
12  USE traclmdz_mod
13  USE lmdz_infotrac, ONLY: nbtr
14  USE iophy
15  USE lmdz_yomcst
16  USE lmdz_YOECUMF
17  USE lmdz_dimensions, ONLY: iim, jjm, llm, ndm
18  USE lmdz_chem, ONLY: idms, iso2, iso4, ih2s, idmso, imsa, ih2o2, &
19          n_avogadro, masse_s, masse_so4, rho_water, rho_ice
20
21  IMPLICIT NONE
22  !=====================================================================
23  ! Objet : depot humide (lessivage et evaporation) de traceurs
24  ! Inspired by routines of Olivier Boucher (mars 1998)
25  ! author R. Pilon 10 octobre 2012
26  ! last modification 16/01/2013 (reformulation partie evaporation)
27  !=====================================================================
28  ! SPLA version taken from trunk revision 2041
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    !JE to speed up, commented 20140219
127
128    !      OPEN(99,file='lsc_scav_param.data',status='old', &
129    !                form='formatted',err=9999)
130    !      READ(99,*,end=9998)  alpha_r
131    !      READ(99,*,end=9998)  alpha_s
132    !      READ(99,*,end=9998)  R_r
133    !      READ(99,*,end=9998)  R_s
134    !      READ(99,*,end=9998)  frac_fine_scav
135    !      READ(99,*,end=9998)  frac_coar_scav
136    !      READ(99,*,end=9998)  frac_aer
137    !9998  Continue
138    !      CLOSE(99)
139    !9999  Continue
140
141    !   PRINT*,'JE alpha_r',alpha_r
142    !   PRINT*,'JE alpha_s',alpha_s
143    !   PRINT*,'JE R_r',R_r
144    !   PRINT*,'JE R_s',R_s
145    !   PRINT*,'frac_fine_scav',frac_fine_scav
146    !   PRINT*,'frac_coar_scav',frac_coar_scav
147    !   PRINT*,'frac_aer ev',frac_aer
148
149    ! JE endcomment
150
151  ENDIF !(debut)
152  !!!!!!!!!!!!!!!!!!!!!!!!!!!
153
154  ! initialization
155  dxin = 0.
156  dxev = 0.
157  beta_ev = 0.
158
159  DO i = 1, klon
160    his_dh(i) = 0.
161  ENDDO
162
163  DO k = 1, klev
164    DO i = 1, klon
165      dxbc(i, k) = 0.
166      beta_v1(i, k) = 0.
167      deltaP(i, k) = 0.
168    ENDDO
169  ENDDO
170
171  DO k = 1, klev
172    DO i = 1, klon
173      d_tr_insc(i, k, it) = 0.
174      d_tr_bcscav(i, k, it) = 0.
175      d_tr_evap(i, k, it) = 0.
176    ENDDO
177  ENDDO
178
179  !  pressure and size of the layer
180  DO k = klev, 1, -1
181    DO i = 1, klon
182      zrho(i, k) = pplay(i, k) / t(i, k) / RD
183      zdz(i, k) = (paprs(i, k) - paprs(i, k + 1)) / zrho(i, k) / RG
184      zmass(i, k) = (paprs(i, k) - paprs(i, k + 1)) / RG
185    ENDDO
186  ENDDO
187
188  !JE<<
189  IF (it==id_prec) THEN                               !  gas
190    frac_ev = frac_gas
191  ELSE                                   !aerosol
192    frac_ev = frac_aer
193  ENDIF
194
195  IF (it==id_prec) THEN                               !  gas
196    DO k = 1, klev
197      DO i = 1, klon
198        henry_t = henry(it) * exp(-kk(it) * (1. / 298. - 1. / t(i, k)))    !  mol/l/atm
199        f_a = henry_t / 101.325 * R * t(i, k) * oliq * zrho(i, k) / rho_water
200        scav(i, k) = f_a / (1. + f_a)
201      ENDDO
202    ENDDO
203  ELSE                            !aerosol
204    DO k = 1, klev
205      DO i = 1, klon
206        scav(i, k) = frac_fine_scav
207      ENDDO
208    ENDDO
209  ENDIF
210  !JE>>
211  DO k = klev - 1, 1, -1
212    DO i = 1, klon
213      !  incloud scavenging
214      !   IF(inscav_fisrt) THEN
215      IF (iflag_lscav == 4) THEN
216        beta = beta_fisrt(i, k) * rneb(i, k)
217      else
218        beta = flxr(i, k) - flxr(i, k + 1) + flxs(i, k) - flxs(i, k + 1)
219        !      beta=beta/zdz(i,k)/oliq/zrho(i,k)
220        beta = beta / zmass(i, k) / oliq
221        beta = MAX(0., beta)
222      endif ! (iflag_lscav .EQ. 4)
223      beta_v1(i, k) = beta    !! for output
224
225      dxin = tr_seri(i, k, it) * (exp(-scav(i, k) * beta * pdtime) - 1.)
226      !      his_dh(i)=his_dh(i)-dxin*zrho(i,k)*zdz(i,k)/pdtime !  kg/m2/s
227      his_dh(i) = his_dh(i) - dxin * zmass(i, k) / pdtime !  kg/m2/s
228      d_tr_insc(i, k, it) = dxin
229
230      !  below-cloud impaction
231      IF(it==id_prec) THEN
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(it) / R_r / rho_water
237        ice = ps * alpha_s(it) / 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)>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))>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      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.