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