source: LMDZ6/trunk/libf/phylmd/lsc_scav.f90 @ 5277

Last change on this file since 5277 was 5274, checked in by abarral, 7 hours ago

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