| [5099] | 1 | |
|---|
| [3526] | 2 | ! $Id: coagulate.F90 5117 2024-07-24 14:23:34Z abarral $ |
|---|
| [5099] | 3 | |
|---|
| [2690] | 4 | SUBROUTINE COAGULATE(pdtcoag,mdw,tr_seri,t_seri,pplay,dens_aer,is_strato) |
|---|
| 5 | ! ----------------------------------------------------------------------- |
|---|
| [5099] | 6 | |
|---|
| [2690] | 7 | ! Author : Christoph Kleinschmitt (with Olivier Boucher) |
|---|
| 8 | ! ------ |
|---|
| [5099] | 9 | |
|---|
| [2690] | 10 | ! purpose |
|---|
| 11 | ! ------- |
|---|
| [5099] | 12 | |
|---|
| [2690] | 13 | ! interface |
|---|
| 14 | ! --------- |
|---|
| 15 | ! input |
|---|
| 16 | ! pdtphys time step duration [sec] |
|---|
| 17 | ! tr_seri tracer mixing ratios [kg/kg] |
|---|
| 18 | ! mdw # or mass median diameter [m] |
|---|
| [5099] | 19 | |
|---|
| [2690] | 20 | ! method |
|---|
| 21 | ! ------ |
|---|
| [5099] | 22 | |
|---|
| [2690] | 23 | ! ----------------------------------------------------------------------- |
|---|
| 24 | |
|---|
| [5101] | 25 | USE dimphy, ONLY: klon,klev |
|---|
| [2690] | 26 | USE aerophys |
|---|
| [3677] | 27 | USE infotrac_phy |
|---|
| [4950] | 28 | USE phys_local_var_mod, ONLY: DENSO4, DENSO4B, f_r_wet, f_r_wetB |
|---|
| 29 | USE strataer_local_var_mod, ONLY: flag_new_strat_compo |
|---|
| [5098] | 30 | USE lmdz_yomcst |
|---|
| [4950] | 31 | |
|---|
| [2690] | 32 | IMPLICIT NONE |
|---|
| 33 | |
|---|
| 34 | !-------------------------------------------------------- |
|---|
| 35 | |
|---|
| 36 | ! transfer variables when calling this routine |
|---|
| 37 | REAL,INTENT(IN) :: pdtcoag ! Time step in coagulation routine [s] |
|---|
| 38 | REAL,DIMENSION(nbtr_bin),INTENT(IN) :: mdw ! aerosol particle diameter in each bin [m] |
|---|
| 39 | REAL,DIMENSION(klon,klev,nbtr),INTENT(INOUT) :: tr_seri ! Concentration Traceur [U/KgA] |
|---|
| 40 | REAL,DIMENSION(klon,klev),INTENT(IN) :: t_seri ! Temperature |
|---|
| 41 | REAL,DIMENSION(klon,klev),INTENT(IN) :: pplay ! pression pour le mileu de chaque couche (en Pa) |
|---|
| 42 | REAL,DIMENSION(klon,klev) :: dens_aer! density of aerosol [kg/m3 aerosol] with default H2SO4 mass |
|---|
| 43 | LOGICAL,DIMENSION(klon,klev),INTENT(IN) :: is_strato |
|---|
| 44 | |
|---|
| 45 | ! local variables in coagulation routine |
|---|
| 46 | INTEGER :: i,j,k,nb,ilon,ilev |
|---|
| [4950] | 47 | REAL, DIMENSION(nbtr_bin) :: radiusdry ! dry aerosol particle radius in each bin [m] |
|---|
| 48 | REAL, DIMENSION(nbtr_bin) :: radiuswet ! wet aerosol particle radius in each bin [m] |
|---|
| [2690] | 49 | REAL, DIMENSION(klon,klev,nbtr_bin) :: tr_t ! Concentration Traceur at time t [U/KgA] |
|---|
| 50 | REAL, DIMENSION(klon,klev,nbtr_bin) :: tr_tp1 ! Concentration Traceur at time t+1 [U/KgA] |
|---|
| 51 | REAL, DIMENSION(nbtr_bin,nbtr_bin,nbtr_bin) :: ff ! Volume fraction of intermediate particles |
|---|
| [4950] | 52 | REAL, DIMENSION(nbtr_bin) :: Vdry ! Volume dry of bins |
|---|
| 53 | REAL, DIMENSION(nbtr_bin) :: Vwet ! Volume wet of bins |
|---|
| [2690] | 54 | REAL, DIMENSION(nbtr_bin,nbtr_bin) :: Vij ! Volume sum of i and j |
|---|
| 55 | REAL :: eta ! Dynamic viscosity of air |
|---|
| 56 | REAL, PARAMETER :: mair=4.8097E-26 ! Average mass of an air molecule [kg] |
|---|
| 57 | REAL :: zrho ! Density of air |
|---|
| 58 | REAL :: mnfrpth ! Mean free path of air |
|---|
| 59 | REAL, DIMENSION(nbtr_bin) :: Kn ! Knudsen number of particle i |
|---|
| 60 | REAL, DIMENSION(nbtr_bin) :: Di ! Particle diffusion coefficient |
|---|
| 61 | REAL, DIMENSION(nbtr_bin) :: m_par ! Mass of particle i |
|---|
| 62 | REAL, DIMENSION(nbtr_bin) :: thvelpar! Thermal velocity of particle i |
|---|
| 63 | REAL, DIMENSION(nbtr_bin) :: mfppar ! Mean free path of particle i |
|---|
| 64 | REAL, DIMENSION(nbtr_bin) :: delta! delta of particle i (from equation 21) |
|---|
| 65 | REAL, DIMENSION(nbtr_bin,nbtr_bin) :: beta ! Coagulation kernel from Brownian diffusion |
|---|
| 66 | REAL :: beta_const ! Constant coagulation kernel (for comparison) |
|---|
| 67 | REAL :: num |
|---|
| 68 | REAL :: numi |
|---|
| 69 | REAL :: denom |
|---|
| 70 | |
|---|
| [2949] | 71 | ! Additional variables for coagulation enhancement factor due to van der Waals forces |
|---|
| 72 | ! Taken from Chan and Mozurkewich, Measurement of the coagulation rate constant for sulfuric acid |
|---|
| 73 | ! particles as a function of particle size using TDMA, Aerosol Science, 32, 321-339, 2001. |
|---|
| 74 | !--ok_vdw is 0 for no vdW forces, 1 for E(0), 2 for E(infinity) |
|---|
| 75 | INTEGER, PARAMETER :: ok_vdw = 0 |
|---|
| 76 | REAL, PARAMETER :: avdW1 = 0.0757 |
|---|
| 77 | REAL, PARAMETER :: avdW3 = 0.0015 |
|---|
| 78 | REAL, PARAMETER :: bvdW0 = 0.0151 |
|---|
| 79 | REAL, PARAMETER :: bvdW1 = -0.186 |
|---|
| 80 | REAL, PARAMETER :: bvdW3 = -0.0163 |
|---|
| 81 | REAL, PARAMETER :: AvdW = 6.4e-20 !Hamaker constant in J = 1e7 erg |
|---|
| 82 | REAL :: AvdWi |
|---|
| 83 | REAL :: xvdW |
|---|
| 84 | REAL :: EvdW |
|---|
| 85 | |
|---|
| [4950] | 86 | ! ff(i,j,k): Volume fraction of Vi,j that is partitioned to each model bin k |
|---|
| 87 | ! just need to be calculated in model initialization because mdw(:) size is fixed |
|---|
| 88 | ! no need to recalculate radius, Vdry, Vij, and ff every timestep because it is for |
|---|
| 89 | ! dry aerosols |
|---|
| [2690] | 90 | DO i=1, nbtr_bin |
|---|
| [4950] | 91 | radiusdry(i)=mdw(i)/2. |
|---|
| 92 | Vdry(i)=radiusdry(i)**3. !neglecting factor 4*RPI/3 |
|---|
| 93 | Vwet(i)=0.0 |
|---|
| [2690] | 94 | ENDDO |
|---|
| 95 | |
|---|
| 96 | DO j=1, nbtr_bin |
|---|
| [4950] | 97 | DO i=1, nbtr_bin |
|---|
| 98 | Vij(i,j)= Vdry(i)+Vdry(j) |
|---|
| 99 | ENDDO |
|---|
| [2690] | 100 | ENDDO |
|---|
| [4950] | 101 | |
|---|
| [2690] | 102 | !--pre-compute the f(i,j,k) from Jacobson equation 13 |
|---|
| 103 | ff=0.0 |
|---|
| 104 | DO k=1, nbtr_bin |
|---|
| 105 | DO j=1, nbtr_bin |
|---|
| 106 | DO i=1, nbtr_bin |
|---|
| [5082] | 107 | IF (k==1) THEN |
|---|
| [2690] | 108 | ff(i,j,k)= 0.0 |
|---|
| [5082] | 109 | ELSEIF (k>1.AND.Vdry(k-1)<Vij(i,j).AND.Vij(i,j)<Vdry(k)) THEN |
|---|
| [2690] | 110 | ff(i,j,k)= 1.-ff(i,j,k-1) |
|---|
| [5082] | 111 | ELSEIF (k==nbtr_bin) THEN |
|---|
| 112 | IF (Vij(i,j)>=Vdry(k)) THEN |
|---|
| [2690] | 113 | ff(i,j,k)= 1. |
|---|
| 114 | ELSE |
|---|
| 115 | ff(i,j,k)= 0.0 |
|---|
| 116 | ENDIF |
|---|
| [5082] | 117 | ELSEIF (k<=(nbtr_bin-1).AND.Vdry(k)<=Vij(i,j).AND.Vij(i,j)<Vdry(k+1)) THEN |
|---|
| [4950] | 118 | ff(i,j,k)= Vdry(k)/Vij(i,j)*(Vdry(k+1)-Vij(i,j))/(Vdry(k+1)-Vdry(k)) |
|---|
| [2690] | 119 | ENDIF |
|---|
| 120 | ENDDO |
|---|
| 121 | ENDDO |
|---|
| 122 | ENDDO |
|---|
| [4950] | 123 | ! End of just need to be calculated at initialization because mdw(:) size is fixed |
|---|
| 124 | |
|---|
| [2690] | 125 | DO ilon=1, klon |
|---|
| 126 | DO ilev=1, klev |
|---|
| 127 | !only in the stratosphere |
|---|
| 128 | IF (is_strato(ilon,ilev)) THEN |
|---|
| 129 | !compute actual wet particle radius & volume for every grid box |
|---|
| [4950] | 130 | IF(flag_new_strat_compo) THEN |
|---|
| 131 | DO i=1, nbtr_bin |
|---|
| 132 | radiuswet(i)=f_r_wetB(ilon,ilev,i)*mdw(i)/2. |
|---|
| 133 | Vwet(i)= radiuswet(i)**3. !neglecting factor 4*RPI/3 |
|---|
| 134 | !! Vwet(i)= Vdry(i)*(f_r_wetB(ilon,ilev,i)**3) |
|---|
| 135 | ENDDO |
|---|
| 136 | ELSE |
|---|
| 137 | DO i=1, nbtr_bin |
|---|
| 138 | radiuswet(i)=f_r_wet(ilon,ilev)*mdw(i)/2. |
|---|
| 139 | Vwet(i)= radiuswet(i)**3. !neglecting factor 4*RPI/3 |
|---|
| 140 | !! Vwet(i)= Vdry(i)*(f_r_wet(ilon,ilev)**3) |
|---|
| 141 | ENDDO |
|---|
| 142 | ENDIF |
|---|
| 143 | |
|---|
| [2690] | 144 | !--Calculations for the coagulation kernel--------------------------------------------------------- |
|---|
| 145 | |
|---|
| 146 | zrho=pplay(ilon,ilev)/t_seri(ilon,ilev)/RD |
|---|
| 147 | |
|---|
| 148 | !--initialize the tracer at time t and convert from [number/KgA] to [number/m3] |
|---|
| 149 | DO i=1, nbtr_bin |
|---|
| 150 | tr_t(ilon,ilev,i) = tr_seri(ilon,ilev,i+nbtr_sulgas) * zrho |
|---|
| 151 | ENDDO |
|---|
| 152 | |
|---|
| 153 | ! mean free path of air (Pruppacher and Klett, 2010, p.417) [m] |
|---|
| 154 | mnfrpth=6.6E-8*(1.01325E+5/pplay(ilon,ilev))*(t_seri(ilon,ilev)/293.15) |
|---|
| 155 | ! mnfrpth=2.*eta/(zrho*thvelair) |
|---|
| 156 | ! mean free path of air (Prupp. Klett) in [10^-6 m] |
|---|
| 157 | ! ZLAIR = 0.066 *(1.01325E+5/PPLAY)*(T_SERI/293.15)*1.E-06 |
|---|
| 158 | |
|---|
| 159 | ! dynamic viscosity of air (Pruppacher and Klett, 2010, p.417) [kg/(m*s)] |
|---|
| [5082] | 160 | IF (t_seri(ilon,ilev)>=273.15) THEN |
|---|
| [2690] | 161 | eta=(1.718+0.0049*(t_seri(ilon,ilev)-273.15))*1.E-5 |
|---|
| 162 | ELSE |
|---|
| 163 | eta=(1.718+0.0049*(t_seri(ilon,ilev)-273.15)-1.2E-5*(t_seri(ilon,ilev)-273.15)**2)*1.E-5 |
|---|
| 164 | ENDIF |
|---|
| 165 | |
|---|
| 166 | !--pre-compute the particle diffusion coefficient Di(i) from equation 18 |
|---|
| 167 | Di=0.0 |
|---|
| 168 | DO i=1, nbtr_bin |
|---|
| [4950] | 169 | Kn(i)=mnfrpth/radiuswet(i) |
|---|
| 170 | Di(i)=RKBOL*t_seri(ilon,ilev)/(6.*RPI*radiuswet(i)*eta)*(1.+Kn(i)*(1.249+0.42*exp(-0.87/Kn(i)))) |
|---|
| [2690] | 171 | ENDDO |
|---|
| 172 | |
|---|
| 173 | !--pre-compute the thermal velocity of a particle thvelpar(i) from equation 20 |
|---|
| 174 | thvelpar=0.0 |
|---|
| [4950] | 175 | IF(flag_new_strat_compo) THEN |
|---|
| 176 | DO i=1, nbtr_bin |
|---|
| 177 | m_par(i)=4./3.*RPI*radiuswet(i)**3.*DENSO4B(ilon,ilev,i)*1000. |
|---|
| 178 | thvelpar(i)=sqrt(8.*RKBOL*t_seri(ilon,ilev)/(RPI*m_par(i))) |
|---|
| 179 | ENDDO |
|---|
| 180 | ELSE |
|---|
| 181 | DO i=1, nbtr_bin |
|---|
| 182 | m_par(i)=4./3.*RPI*radiuswet(i)**3.*DENSO4(ilon,ilev)*1000. |
|---|
| 183 | thvelpar(i)=sqrt(8.*RKBOL*t_seri(ilon,ilev)/(RPI*m_par(i))) |
|---|
| 184 | ENDDO |
|---|
| 185 | ENDIF |
|---|
| [2690] | 186 | |
|---|
| 187 | !--pre-compute the particle mean free path mfppar(i) from equation 22 |
|---|
| 188 | mfppar=0.0 |
|---|
| 189 | DO i=1, nbtr_bin |
|---|
| 190 | mfppar(i)=8.*Di(i)/(RPI*thvelpar(i)) |
|---|
| 191 | ENDDO |
|---|
| 192 | |
|---|
| 193 | !--pre-compute the mean distance delta(i) from the center of a sphere reached by particles |
|---|
| 194 | !--leaving the surface of the sphere and traveling a distance of particle mfppar(i) from equation 21 |
|---|
| 195 | delta=0.0 |
|---|
| 196 | DO i=1, nbtr_bin |
|---|
| [4950] | 197 | delta(i)=((2.*radiuswet(i)+mfppar(i))**3.-(4.*radiuswet(i)**2.+mfppar(i)**2.)**1.5)/ & |
|---|
| [5087] | 198 | (6.*radiuswet(i)*mfppar(i))-2.*radiuswet(i) |
|---|
| [2690] | 199 | ENDDO |
|---|
| 200 | |
|---|
| [4950] | 201 | ! beta(i,j): coagulation kernel (rate coefficient) of 2 colliding particles i,j |
|---|
| [2949] | 202 | !--pre-compute the beta(i,j) from equation 17 in Jacobson |
|---|
| [2690] | 203 | num=0.0 |
|---|
| 204 | DO j=1, nbtr_bin |
|---|
| 205 | DO i=1, nbtr_bin |
|---|
| [5099] | 206 | |
|---|
| [4950] | 207 | num=4.*RPI*(radiuswet(i)+radiuswet(j))*(Di(i)+Di(j)) |
|---|
| 208 | denom=(radiuswet(i)+radiuswet(j))/(radiuswet(i)+radiuswet(j)+sqrt(delta(i)**2.+delta(j)**2.))+ & |
|---|
| [5087] | 209 | 4.*(Di(i)+Di(j))/(sqrt(thvelpar(i)**2.+thvelpar(j)**2.)*(radiuswet(i)+radiuswet(j))) |
|---|
| [4950] | 210 | beta(i,j)=num/denom |
|---|
| [5099] | 211 | |
|---|
| [2949] | 212 | !--compute enhancement factor due to van der Waals forces |
|---|
| [5082] | 213 | IF (ok_vdw == 0) THEN !--no enhancement factor |
|---|
| [4950] | 214 | Evdw=1.0 |
|---|
| [5082] | 215 | ELSEIF (ok_vdw == 1) THEN !--E(0) case |
|---|
| [4950] | 216 | AvdWi = AvdW/(RKBOL*t_seri(ilon,ilev))*(4.*radiuswet(i)*radiuswet(j))/(radiuswet(i)+radiuswet(j))**2. |
|---|
| 217 | xvdW = LOG(1.+AvdWi) |
|---|
| 218 | EvdW = 1. + avdW1*xvdW + avdW3*xvdW**3 |
|---|
| [5082] | 219 | ELSEIF (ok_vdw == 2) THEN !--E(infinity) case |
|---|
| [4950] | 220 | AvdWi = AvdW/(RKBOL*t_seri(ilon,ilev))*(4.*radiuswet(i)*radiuswet(j))/(radiuswet(i)+radiuswet(j))**2. |
|---|
| 221 | xvdW = LOG(1.+AvdWi) |
|---|
| 222 | EvdW = 1. + SQRT(AvdWi/3.)/(1.+bvdW0*SQRT(AvdWi)) + bvdW1*xvdW + bvdW3*xvdW**3. |
|---|
| [2949] | 223 | ENDIF |
|---|
| [5099] | 224 | |
|---|
| [2949] | 225 | beta(i,j)=beta(i,j)*EvdW |
|---|
| 226 | |
|---|
| [2690] | 227 | ENDDO |
|---|
| 228 | ENDDO |
|---|
| 229 | |
|---|
| 230 | !--external loop for equation 14 |
|---|
| 231 | DO k=1, nbtr_bin |
|---|
| 232 | |
|---|
| 233 | !--calculating denominator sum |
|---|
| 234 | denom=0.0 |
|---|
| 235 | DO j=1, nbtr_bin |
|---|
| [4950] | 236 | ! fraction of coagulation of k and j that is not giving k |
|---|
| 237 | denom=denom+(1.-ff(k,j,k))*beta(k,j)*tr_t(ilon,ilev,j) |
|---|
| [2690] | 238 | ENDDO |
|---|
| 239 | |
|---|
| [5082] | 240 | IF (k==1) THEN |
|---|
| [2690] | 241 | !--calculate new concentration of smallest bin |
|---|
| 242 | tr_tp1(ilon,ilev,k)=tr_t(ilon,ilev,k)/(1.+pdtcoag*denom) |
|---|
| 243 | ELSE |
|---|
| 244 | !--calculating double sum terms in numerator of eq 14 |
|---|
| 245 | num=0.0 |
|---|
| 246 | DO j=1, k |
|---|
| [4950] | 247 | numi=0.0 |
|---|
| 248 | DO i=1, k-1 |
|---|
| [5099] | 249 | |
|---|
| [4950] | 250 | ! see Jacobson: " In order to conserve volume and volume concentration (which |
|---|
| 251 | ! coagulation physically does) while giving up some accuracy in number concentration" |
|---|
| [5099] | 252 | |
|---|
| [4950] | 253 | ! Coagulation of i and j giving k |
|---|
| 254 | ! with V(i) and then V(j) because it considers i,j and j,i with the double loop |
|---|
| [5099] | 255 | |
|---|
| [4950] | 256 | ! BUT WHY WET VOLUME V(i) in old STRATAER? tracers are already dry aerosols and coagulation |
|---|
| 257 | ! kernel beta(i,j) accounts for wet aerosols -> reply below |
|---|
| [5099] | 258 | |
|---|
| [4950] | 259 | ! numi=numi+ff(i,j,k)*beta(i,j)*V(i)*tr_tp1(ilon,ilev,i)*tr_t(ilon,ilev,j) |
|---|
| 260 | numi=numi+ff(i,j,k)*beta(i,j)*Vdry(i)*tr_tp1(ilon,ilev,i)*tr_t(ilon,ilev,j) |
|---|
| 261 | ENDDO |
|---|
| 262 | num=num+numi |
|---|
| [2690] | 263 | ENDDO |
|---|
| 264 | |
|---|
| 265 | !--calculate new concentration of other bins |
|---|
| [4950] | 266 | ! tr_tp1(ilon,ilev,k)=(V(k)*tr_t(ilon,ilev,k)+pdtcoag*num)/( (1.+pdtcoag*denom)*V(k) ) |
|---|
| 267 | tr_tp1(ilon,ilev,k)=(Vdry(k)*tr_t(ilon,ilev,k)+pdtcoag*num)/( (1.+pdtcoag*denom)*Vdry(k) ) |
|---|
| [5099] | 268 | |
|---|
| [4950] | 269 | ! In constant composition (no dependency on aerosol size because no kelvin effect) |
|---|
| 270 | ! V(l)= (f_r_wet(ilon,ilev)**3)*((mdw(l)/2.)**3) = (f_r_wet(ilon,ilev)**3)*Vdry(i) |
|---|
| 271 | ! so numi and num are proportional (f_r_wet(ilon,ilev)**3) |
|---|
| 272 | ! and so |
|---|
| 273 | ! tr_tp1(ilon,ilev,k)=(V(k)*tr_t(ilon,ilev,k)+pdtcoag*num)/( (1.+pdtcoag*denom)*V(k) ) |
|---|
| 274 | ! =(Vdry(k)*tr_t(ilon,ilev,k)+pdtcoag*num_dry)/( (1.+pdtcoag*denom)*Vdry(k) ) |
|---|
| 275 | ! with num_dry=...beta(i,j)*Vdry(i)*.... |
|---|
| [5117] | 276 | ! so in old STRATAER (.NOT.flag_new_strat_compo), it was correct |
|---|
| [2690] | 277 | ENDIF |
|---|
| 278 | |
|---|
| 279 | ENDDO !--end of loop k |
|---|
| 280 | |
|---|
| 281 | !--convert tracer concentration back from [number/m3] to [number/KgA] and write into tr_seri |
|---|
| 282 | DO i=1, nbtr_bin |
|---|
| [4950] | 283 | tr_seri(ilon,ilev,i+nbtr_sulgas) = tr_tp1(ilon,ilev,i) / zrho |
|---|
| [2690] | 284 | ENDDO |
|---|
| 285 | |
|---|
| 286 | ENDIF ! IF IN STRATOSPHERE |
|---|
| 287 | ENDDO !--end of loop klev |
|---|
| 288 | ENDDO !--end of loop klon |
|---|
| [4950] | 289 | ! ********************************************* |
|---|
| [2690] | 290 | |
|---|
| 291 | END SUBROUTINE COAGULATE |
|---|
