source: LMDZ6/branches/Amaury_dev/libf/phylmd/lsc_scav.F90 @ 5157

Last change on this file since 5157 was 5144, checked in by abarral, 7 weeks ago

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