[2175] | 1 | !$Id $ |
---|
| 2 | |
---|
| 3 | SUBROUTINE lsc_scav_orig(pdtime,it,iflag_lscav,oliq,flxr,flxs,rneb,beta_fisrt, & |
---|
| 4 | beta_v1,pplay,paprs,t,tr_seri,d_tr_insc, & |
---|
| 5 | d_tr_bcscav,d_tr_evap,qPrls) |
---|
| 6 | USE ioipsl |
---|
| 7 | USE dimphy |
---|
| 8 | USE mod_grid_phy_lmdz |
---|
| 9 | USE mod_phys_lmdz_para |
---|
| 10 | USE traclmdz_mod |
---|
| 11 | USE infotrac,ONLY : nbtr |
---|
| 12 | USE comgeomphy |
---|
| 13 | USE iophy |
---|
| 14 | |
---|
| 15 | IMPLICIT NONE |
---|
| 16 | !===================================================================== |
---|
| 17 | ! Objet : depot humide (lessivage et evaporation) de traceurs |
---|
| 18 | ! Inspired by routines of Olivier Boucher (mars 1998) |
---|
| 19 | ! author R. Pilon 10 octobre 2012 |
---|
| 20 | ! last modification 16/01/2013 (reformulation partie evaporation) |
---|
| 21 | !===================================================================== |
---|
| 22 | |
---|
| 23 | include "dimensions.h" |
---|
| 24 | include "chem.h" |
---|
| 25 | include "YOMCST.h" |
---|
| 26 | include "YOECUMF.h" |
---|
| 27 | |
---|
| 28 | REAL,INTENT(IN) :: pdtime ! time step (s) |
---|
| 29 | INTEGER,INTENT(IN) :: it ! tracer number |
---|
| 30 | INTEGER,INTENT(IN) :: iflag_lscav ! LS scavenging param: |
---|
| 31 | ! 3=Reddy_Boucher2004, 4=3+RPilon. |
---|
| 32 | REAL,DIMENSION(klon,klev+1),INTENT(IN) :: flxr ! flux precipitant de pluie |
---|
| 33 | REAL,DIMENSION(klon,klev+1),INTENT(IN) :: flxs ! flux precipitant de neige |
---|
| 34 | REAL,INTENT(IN) :: oliq ! contenu en eau liquide dans le nuage (kg/kg) |
---|
| 35 | REAL,DIMENSION(klon,klev),INTENT(IN) :: rneb |
---|
| 36 | REAL,DIMENSION(klon,klev),INTENT(IN) :: pplay ! pression |
---|
| 37 | REAL,DIMENSION(klon,klev+1),INTENT(IN) :: paprs ! pression |
---|
| 38 | REAL,DIMENSION(klon,klev),INTENT(IN) :: t ! temperature |
---|
| 39 | ! tracers |
---|
| 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 | |
---|
| 52 | ! variables locales |
---|
| 53 | LOGICAL,SAVE :: debut=.true. |
---|
| 54 | !$OMP THREADPRIVATE(debut) |
---|
| 55 | ! |
---|
| 56 | REAL,PARAMETER :: henry=1.4 ! constante de Henry en mol/l/atm ~1.4 for gases |
---|
| 57 | REAL :: henry_t ! constante de Henry a T t (mol/l/atm) |
---|
| 58 | REAL,PARAMETER :: kk=2900. ! coefficient de dependence en T (K) |
---|
| 59 | REAL :: f_a ! rapport de la phase aqueuse a la phase gazeuse |
---|
| 60 | REAL :: beta ! taux de conversion de l'eau en pluie |
---|
| 61 | |
---|
| 62 | INTEGER :: i, k |
---|
| 63 | REAL,DIMENSION(klon,klev) :: scav ! water liquid content / fraction aqueuse du constituant |
---|
| 64 | REAL,DIMENSION(klon,klev) :: zrho |
---|
| 65 | REAL,DIMENSION(klon,klev) :: zdz |
---|
| 66 | REAL,DIMENSION(klon,klev) :: zmass ! layer mass |
---|
| 67 | |
---|
| 68 | REAL :: frac_ev ! cste pour la reevaporation : dropplet shrinking |
---|
| 69 | ! frac_ev = frac_gas ou frac_aer |
---|
| 70 | REAL,PARAMETER :: frac_gas=1.0 ! cste pour la reevaporation pour les gaz |
---|
| 71 | REAL :: frac_aer ! cste pour la reevaporation pour les particules |
---|
| 72 | REAL,DIMENSION(klon,klev) :: deltaP ! P(i+1)-P(i) |
---|
| 73 | REAL,DIMENSION(klon,klev) :: beta_ev ! dP/P(i+1) |
---|
| 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 | !!!!!!!!!!!!!!!!!!!! choix lessivage !!!!!!!!!!!!!!!!!!!!!!!! |
---|
| 93 | !! logical,save :: inscav_fisrt |
---|
| 94 | !!! $OMP THREADPRIVATE(inscav_first) |
---|
| 95 | ! |
---|
| 96 | !!!!!!!!!!!!!!!!!!!!!!!!!!! |
---|
| 97 | IF (debut) THEN |
---|
| 98 | ! |
---|
| 99 | ! inscav_fisrt=.true. |
---|
| 100 | ! call getin('inscav_fisrt',inscav_fisrt) |
---|
| 101 | ! if(inscav_fisrt) then |
---|
| 102 | ! print*,'beta from fisrtilp.F90, beta = (z_cond - z_oliq)/z_cond, inscav_fisrt=',inscav_fisrt |
---|
| 103 | ! else |
---|
| 104 | ! print*,'beta from Reddy and Bocuher 2004 (original version), inscav_fisrt=',inscav_fisrt |
---|
| 105 | ! endif |
---|
| 106 | ! |
---|
| 107 | alpha_r=0.001 ! coefficient d'impaction pour la pluie |
---|
| 108 | alpha_s=0.01 ! coefficient d'impaction pour la neige |
---|
| 109 | R_r=0.001 ! mean raindrop radius (m) |
---|
| 110 | R_s=0.001 ! mean snow crystal radius (m) |
---|
| 111 | frac_fine_scav=0.7 |
---|
| 112 | frac_coar_scav=0.7 |
---|
| 113 | ! frac_aer=0.5 ~ droplet size shrinks by evap |
---|
| 114 | frac_aer=0.5 |
---|
| 115 | ! |
---|
| 116 | |
---|
| 117 | !JE to speed up, commented 20140219 |
---|
| 118 | ! |
---|
| 119 | ! OPEN(99,file='lsc_scav_param.data',status='old', & |
---|
| 120 | ! form='formatted',err=9999) |
---|
| 121 | ! READ(99,*,end=9998) alpha_r |
---|
| 122 | ! READ(99,*,end=9998) alpha_s |
---|
| 123 | ! READ(99,*,end=9998) R_r |
---|
| 124 | ! READ(99,*,end=9998) R_s |
---|
| 125 | ! READ(99,*,end=9998) frac_fine_scav |
---|
| 126 | ! READ(99,*,end=9998) frac_coar_scav |
---|
| 127 | ! READ(99,*,end=9998) frac_aer |
---|
| 128 | !9998 Continue |
---|
| 129 | ! CLOSE(99) |
---|
| 130 | !9999 Continue |
---|
| 131 | ! |
---|
| 132 | ! print*,'alpha_r',alpha_r |
---|
| 133 | ! print*,'alpha_s',alpha_s |
---|
| 134 | ! print*,'R_r',R_r |
---|
| 135 | ! print*,'R_s',R_s |
---|
| 136 | ! print*,'frac_fine_scav',frac_fine_scav |
---|
| 137 | ! print*,'frac_coar_scav',frac_coar_scav |
---|
| 138 | ! print*,'frac_aer ev',frac_aer |
---|
| 139 | ! |
---|
| 140 | ! JE endcomment |
---|
| 141 | ! |
---|
| 142 | ENDIF !(debut) |
---|
| 143 | !!!!!!!!!!!!!!!!!!!!!!!!!!! |
---|
| 144 | ! |
---|
| 145 | ! initialization |
---|
| 146 | dxin=0. |
---|
| 147 | dxev=0. |
---|
| 148 | beta_ev=0. |
---|
| 149 | |
---|
| 150 | DO i=1,klon |
---|
| 151 | his_dh(i)=0. |
---|
| 152 | ENDDO |
---|
| 153 | |
---|
| 154 | DO k=1,klev |
---|
| 155 | DO i=1, klon |
---|
| 156 | dxbc(i,k)=0. |
---|
| 157 | beta_v1(i,k)=0. |
---|
| 158 | deltaP(i,k)=0. |
---|
| 159 | ENDDO |
---|
| 160 | ENDDO |
---|
| 161 | |
---|
| 162 | DO k=1,klev |
---|
| 163 | DO i=1, klon |
---|
| 164 | d_tr_insc(i,k,it)=0. |
---|
| 165 | d_tr_bcscav(i,k,it)=0. |
---|
| 166 | d_tr_evap(i,k,it)=0. |
---|
| 167 | ENDDO |
---|
| 168 | ENDDO |
---|
| 169 | |
---|
| 170 | ! pressure and size of the layer |
---|
| 171 | DO k=klev-1, 1, -1 |
---|
| 172 | DO i=1, klon |
---|
| 173 | zrho(i,k)=pplay(i,k)/t(i,k)/RD |
---|
| 174 | zdz(i,k)=(paprs(i,k)-paprs(i,k+1))/zrho(i,k)/RG |
---|
| 175 | zmass(i,k)=(paprs(i,k)-paprs(i,k+1))/RG |
---|
| 176 | ENDDO |
---|
| 177 | ENDDO |
---|
| 178 | |
---|
| 179 | IF (it.gt.1) THEN ! aerosol |
---|
| 180 | frac_ev=frac_aer |
---|
| 181 | ELSE ! gas |
---|
| 182 | frac_ev=frac_gas |
---|
| 183 | ENDIF |
---|
| 184 | |
---|
| 185 | IF(it.gt.1) then ! aerosol |
---|
| 186 | DO k=1, klev |
---|
| 187 | DO i=1, klon |
---|
| 188 | scav(i,k)=frac_fine_scav |
---|
| 189 | ENDDO |
---|
| 190 | ENDDO |
---|
| 191 | ELSE ! gas |
---|
| 192 | DO k=1, klev |
---|
| 193 | DO i=1, klon |
---|
| 194 | henry_t=henry*exp(-kk*(1./298.-1./t(i,k))) ! mol/l/atm |
---|
| 195 | f_a=henry_t/101.325*R*t(i,k)*oliq*zrho(i,k)/rho_water |
---|
| 196 | scav(i,k)=f_a/(1.+f_a) |
---|
| 197 | ENDDO |
---|
| 198 | ENDDO |
---|
| 199 | ENDIF |
---|
| 200 | |
---|
| 201 | DO k=klev-1, 1, -1 |
---|
| 202 | DO i=1, klon |
---|
| 203 | ! incloud scavenging |
---|
| 204 | ! if(inscav_fisrt) then |
---|
| 205 | if (iflag_lscav .eq. 4) then |
---|
| 206 | beta=beta_fisrt(i,k)*rneb(i,k) |
---|
| 207 | else |
---|
| 208 | beta=flxr(i,k)-flxr(i,k+1)+flxs(i,k)-flxs(i,k+1) |
---|
| 209 | ! beta=beta/zdz(i,k)/oliq/zrho(i,k) |
---|
| 210 | beta=beta/zmass(i,k)/oliq |
---|
| 211 | beta=MAX(0.,beta) |
---|
| 212 | endif ! (iflag_lscav .eq. 4) |
---|
| 213 | beta_v1(i,k)=beta !! for output |
---|
| 214 | ! |
---|
| 215 | dxin=tr_seri(i,k,it)*(exp(-scav(i,k)*beta*pdtime)-1.) |
---|
| 216 | ! his_dh(i)=his_dh(i)-dxin*zrho(i,k)*zdz(i,k)/pdtime ! kg/m2/s |
---|
| 217 | his_dh(i)=his_dh(i)-dxin*zmass(i,k)/pdtime ! kg/m2/s |
---|
| 218 | d_tr_insc(i,k,it)=dxin |
---|
| 219 | |
---|
| 220 | ! below-cloud impaction |
---|
| 221 | IF(it.eq.1) then |
---|
| 222 | d_tr_bcscav(i,k,it)=0. |
---|
| 223 | ELSE |
---|
| 224 | pr=0.5*(flxr(i,k)+flxr(i,k+1)) |
---|
| 225 | ps=0.5*(flxs(i,k)+flxs(i,k+1)) |
---|
| 226 | water=pr*alpha_r/R_r/rho_water |
---|
| 227 | ice=ps*alpha_s/R_s/rho_ice |
---|
| 228 | dxbc(i,k)=-3./4.*tr_seri(i,k,it)*pdtime*(water+ice) |
---|
| 229 | ! add tracers from below cloud scav in his_dh |
---|
| 230 | his_dh(i)=his_dh(i)-dxbc(i,k)*zmass(i,k)/pdtime ! kg/m2/s |
---|
| 231 | d_tr_bcscav(i,k,it)=dxbc(i,k) |
---|
| 232 | ENDIF |
---|
| 233 | |
---|
| 234 | ! reevaporation |
---|
| 235 | deltaP(i,k)=flxr(i,k+1)+flxs(i,k+1)-flxr(i,k)-flxs(i,k) |
---|
| 236 | deltaP(i,k)=max(deltaP(i,k),0.) |
---|
| 237 | |
---|
| 238 | if(flxr(i,k+1)+flxs(i,k+1).gt.1.e-16) then |
---|
| 239 | beta_ev(i,k)=deltaP(i,k)/(flxr(i,k+1)+flxs(i,k+1)) |
---|
| 240 | else |
---|
| 241 | beta_ev(i,k)=0. |
---|
| 242 | endif |
---|
| 243 | |
---|
| 244 | beta_ev(i,k)=max(min(1.,beta_ev(i,k)),0.) |
---|
| 245 | |
---|
| 246 | !jyg |
---|
| 247 | |
---|
| 248 | if(abs(1-(1-frac_ev)*beta_ev(i,k)).gt.1.e-16) then |
---|
| 249 | ! remove tracers from precipitation owing to release by evaporation in his_dh |
---|
| 250 | ! dxev=frac_ev*beta_ev(i,k)*his_dh(i) *pdtime/(zrho(i,k)*zdz(i,k)) & |
---|
| 251 | dxev=frac_ev*beta_ev(i,k)*his_dh(i) *pdtime/(zmass(i,k)) & |
---|
| 252 | /(1 -(1-frac_ev)*beta_ev(i,k)) |
---|
| 253 | his_dh(i)=his_dh(i)*(1 - frac_ev*beta_ev(i,k) / (1 -(1-frac_ev)*beta_ev(i,k))) |
---|
| 254 | else |
---|
| 255 | ! dxev=his_dh(i) *pdtime/(zrho(i,k)*zdz(i,k)) |
---|
| 256 | dxev=his_dh(i) *pdtime/(zmass(i,k)) |
---|
| 257 | his_dh(i)=0. |
---|
| 258 | endif |
---|
| 259 | ! print*, k, 'beta_ev',beta_ev |
---|
| 260 | ! remove tracers from precipitation owing to release by evaporation in his_dh |
---|
| 261 | !! dxev=frac_ev*deltaP(i,k)*pdtime * his_dh(i) /(zrho(i,k)*zdz(i,k)) |
---|
| 262 | !rplmd |
---|
| 263 | !! dxev=frac_ev*deltaP(i,k)*his_dh(i) *pdtime/(zrho(i,k)*zdz(i,k)) & |
---|
| 264 | !! /max(flxr(i,k)+flxs(i,k),1.e-16) |
---|
| 265 | |
---|
| 266 | |
---|
| 267 | d_tr_evap(i,k,it)=dxev |
---|
| 268 | !! tendency is further added in phytrac x = x + dx |
---|
| 269 | ENDDO !! do i |
---|
| 270 | ENDDO !! do k |
---|
| 271 | |
---|
| 272 | !jyg (20130114) |
---|
| 273 | DO i = 1,klon |
---|
| 274 | qPrls(i,it) = his_dh(i)/max(flxr(i,1)+flxs(i,1),1.e-16) |
---|
| 275 | ENDDO |
---|
| 276 | !jyg end |
---|
| 277 | |
---|
| 278 | |
---|
| 279 | ! test de conservation |
---|
| 280 | conserv=0. |
---|
| 281 | ! DO k= klev,1,-1 |
---|
| 282 | ! DO i=1, klon |
---|
| 283 | ! conserv=conserv+d_tr_insc(i,k,it)*(paprs(i,k)-paprs(i,k+1))/RG & |
---|
| 284 | ! +d_tr_bcscav(i,k,it)*(paprs(i,k)-paprs(i,k+1))/RG & |
---|
| 285 | ! +d_tr_evap(i,k,it)*(paprs(i,k)-paprs(i,k+1))/RG |
---|
| 286 | ! if(it.eq.3) write(*,'(I2,2X,a,e20.12,2X,a,e20.12,2X,a,e20.12,2X,a,e20.12)'),& |
---|
| 287 | ! 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) |
---|
| 288 | ! ENDDO |
---|
| 289 | ! ENDDO |
---|
| 290 | |
---|
| 291 | END SUBROUTINE lsc_scav_orig |
---|