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