| 1 | !======================================================================================================================! |
|---|
| 2 | ! Module: Scheme of co2 cloud formation ===============================================================================! |
|---|
| 3 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 4 | ! Authors: Joaquim Audouard, Constantino Listowski, Anni Määttänen |
|---|
| 5 | ! Date: 09/2016 |
|---|
| 6 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 7 | ! Contains subroutines: |
|---|
| 8 | ! - improvedco2clouds_mod: nucleation |
|---|
| 9 | !======================================================================================================================! |
|---|
| 10 | module improvedco2clouds_mod |
|---|
| 11 | |
|---|
| 12 | implicit none |
|---|
| 13 | |
|---|
| 14 | contains |
|---|
| 15 | !======================================================================================================================! |
|---|
| 16 | ! SUBROUTINE: improvedco2clouds =======================================================================================! |
|---|
| 17 | !======================================================================================================================! |
|---|
| 18 | ! Subject: |
|---|
| 19 | !--------- |
|---|
| 20 | ! This routine is used to form CO2 clouds when a parcel of the GCM is saturated. |
|---|
| 21 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 22 | ! Comments: |
|---|
| 23 | !---------- |
|---|
| 24 | ! Adaptation for CO2 clouds based on improvedclouds_mod.F |
|---|
| 25 | ! |
|---|
| 26 | ! It includes the ability to have supersaturation, a computation of the nucleation rates, growthrates and the |
|---|
| 27 | ! scavenging of dust particles by clouds. It is worth noting that the amount of dust is computed using the dust |
|---|
| 28 | ! optical depth computed in aeropacity.F. |
|---|
| 29 | ! That's why the variable called "tauscaling" is used to convert pq(dust_mass) and pq(dust_number), which are relative |
|---|
| 30 | ! quantities, to absolute and realistic quantities stored in zq. This has to be done to convert the inputs into |
|---|
| 31 | ! absolute values, but also to convert the outputs back into relative values which are then used by the sedimentation |
|---|
| 32 | ! and advection schemes. |
|---|
| 33 | ! CO2 ice particles can nucleate on both dust and water ice particles. When CO2 ice is deposited onto a water ice |
|---|
| 34 | ! particles, the particle is removed from the water tracers. Memory of the origin of the co2 particles is kept and |
|---|
| 35 | ! thus the water cycle shouldn't be modified by this. |
|---|
| 36 | ! There is an energy limit to how much co2 can sublimate/condensate. It is defined by the difference of the GCM |
|---|
| 37 | ! temperature with the co2 condensation temperature. |
|---|
| 38 | ! |
|---|
| 39 | ! /!\ WARNING /!\ |
|---|
| 40 | ! No sedimentation of the water ice origin is performed in the microphysical timestep in co2cloud.F. |
|---|
| 41 | ! |
|---|
| 42 | ! If meteoritic particles are activated and turn into co2 ice particles, then they will be reversed in the dust |
|---|
| 43 | ! tracers if the cloud sublimates. |
|---|
| 44 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 45 | ! Paper: |
|---|
| 46 | !------- |
|---|
| 47 | ! see co2cloud.F90 |
|---|
| 48 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 49 | ! Algorithm: |
|---|
| 50 | !----------- |
|---|
| 51 | ! 0. Firstcall |
|---|
| 52 | ! 0.1. Bonus: meteoritic component, extract data |
|---|
| 53 | ! 1. Initialization |
|---|
| 54 | ! 2. Compute saturation |
|---|
| 55 | ! 3. Bonus: additional meteoritic particles for nucleation |
|---|
| 56 | ! 4. Actual microphysics: Main loop over the GCM's grid |
|---|
| 57 | ! 4.1 Nucleation |
|---|
| 58 | ! 4.2. Ice growth: scheme for radius evolution |
|---|
| 59 | ! 4.3 Dust cores releasing if no more co2 ice |
|---|
| 60 | ! 5. Get cloud tendencies |
|---|
| 61 | !======================================================================================================================! |
|---|
| 62 | subroutine improvedCO2clouds(ngrid, nlay, microtimestep, pplay, pplev, pteff, sum_subpdt, pqeff, sum_subpdq, & |
|---|
| 63 | subpdqcloudco2, subpdtcloudco2, nq, tauscaling, rb_cldco2, sigma_iceco2, dev2) |
|---|
| 64 | |
|---|
| 65 | use comcstfi_h, only: pi, g, cpp |
|---|
| 66 | |
|---|
| 67 | use updaterad, only: updaterice_micro, updaterice_microco2, updaterccnCO2 |
|---|
| 68 | |
|---|
| 69 | use tracer_mod, only: igcm_dust_mass, igcm_dust_number, rho_dust, igcm_h2o_ice, igcm_ccn_mass, igcm_ccn_number, & |
|---|
| 70 | nuice_sed, igcm_co2, igcm_co2_ice, igcm_ccnco2_mass, igcm_ccnco2_number, & |
|---|
| 71 | igcm_ccnco2_h2o_mass_ice, igcm_ccnco2_h2o_mass_ccn, igcm_ccnco2_h2o_number, nuiceco2_sed, & |
|---|
| 72 | nuiceco2_ref, igcm_ccnco2_meteor_mass, igcm_ccnco2_meteor_number |
|---|
| 73 | |
|---|
| 74 | use conc_mod, only: mmean |
|---|
| 75 | use time_phylmdz_mod, only: daysec |
|---|
| 76 | use nucleaco2_mod, only: nucleaco2 |
|---|
| 77 | use datafile_mod, only: datadir |
|---|
| 78 | |
|---|
| 79 | use density_co2_ice_mod, only: density_co2_ice |
|---|
| 80 | |
|---|
| 81 | implicit none |
|---|
| 82 | |
|---|
| 83 | include "callkeys.h" |
|---|
| 84 | include "microphys.h" |
|---|
| 85 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 86 | ! VARIABLES DECLARATION |
|---|
| 87 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 88 | ! Input arguments: |
|---|
| 89 | !----------------- |
|---|
| 90 | integer, intent(in) :: & |
|---|
| 91 | nq, &! number of tracers |
|---|
| 92 | ngrid, &! number of point grid |
|---|
| 93 | nlay ! number of layer |
|---|
| 94 | |
|---|
| 95 | real, intent(in) :: & |
|---|
| 96 | microtimestep, &! physics time step (s) |
|---|
| 97 | pplay(ngrid,nlay), &! mid-layer pressure (Pa) |
|---|
| 98 | pplev(ngrid,nlay+1), &! inter-layer pressure (Pa) |
|---|
| 99 | pteff(ngrid,nlay), &! temperature at the middle of the layers (K) |
|---|
| 100 | sum_subpdt(ngrid,nlay), &! tendency on temperature from previous physical parametrizations |
|---|
| 101 | pqeff(ngrid,nlay,nq), &! tracers (kg/kg) |
|---|
| 102 | tauscaling(ngrid), &! convertion factor for qdust and Ndust |
|---|
| 103 | sum_subpdq(ngrid,nlay,nq) ! tendencies on tracers before condensation (kg/kg.s-1) |
|---|
| 104 | |
|---|
| 105 | real, intent(in) :: & |
|---|
| 106 | sigma_iceco2 ! Variance of the co2 ice and CCN distributions |
|---|
| 107 | |
|---|
| 108 | double precision, intent(in) :: & |
|---|
| 109 | rb_cldco2(nbinco2_cld+1), & ! boundary values of each rad_cldco2 |
|---|
| 110 | dev2 |
|---|
| 111 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 112 | ! Output arguments: |
|---|
| 113 | !------------------ |
|---|
| 114 | real, intent(out) :: & |
|---|
| 115 | subpdtcloudco2(ngrid,nlay), &! tendency on tracers due to CO2 condensation (K/s) |
|---|
| 116 | subpdqcloudco2(ngrid,nlay,nq) ! tendency on tracers due to CO2 condensation (kg/kg.s-1) |
|---|
| 117 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 118 | ! Local: |
|---|
| 119 | !------- |
|---|
| 120 | !----1) Parameters: |
|---|
| 121 | !------------------ |
|---|
| 122 | integer, parameter :: & |
|---|
| 123 | ! ---Meteoritic flux input file |
|---|
| 124 | nbin_meteor = 100, &! Dimension 1 |
|---|
| 125 | nlev_meteor = 130, &! Dimension 2 |
|---|
| 126 | uMeteor = 666, &! File unit |
|---|
| 127 | ! ---Latent heat computation |
|---|
| 128 | l0 = 595594d0, &! coeff from: Azreg-Aïnou (2005) |
|---|
| 129 | l1 = 903.111d0, &! Title: "Low-temperature data for carbon dioxide" |
|---|
| 130 | l2 = -11.5959d0, &! Pulication: eprint arXiv:1403.4403 |
|---|
| 131 | l3 = 0.0528288d0, &! |
|---|
| 132 | l4 = -0.000103183d0 ! |
|---|
| 133 | |
|---|
| 134 | real, parameter :: & |
|---|
| 135 | threshold = 1e-30 ! limit value |
|---|
| 136 | |
|---|
| 137 | character(len=23), parameter:: & |
|---|
| 138 | file_meteoritic_flux = 'Meteo_flux_Plane.dat' |
|---|
| 139 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 140 | !----2) Saved: |
|---|
| 141 | !------------- |
|---|
| 142 | real, save :: & |
|---|
| 143 | sigma_ice ! Variance of the h2o ice and CCN distributions |
|---|
| 144 | |
|---|
| 145 | double precision, save :: & |
|---|
| 146 | pression_meteor(nlev_meteor), &! pressure from meteoritic flux input file |
|---|
| 147 | meteor(nlev_meteor,nbin_meteor), &! Meteoritic flux read from file uMeteor |
|---|
| 148 | dev3 ! 1. / ( sqrt(2.) * sigma_ice ) |
|---|
| 149 | |
|---|
| 150 | logical, save :: & |
|---|
| 151 | firstcall = .true. ! Used to compute saved variables |
|---|
| 152 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 153 | !----3) Variables: |
|---|
| 154 | !----------------- |
|---|
| 155 | integer :: & |
|---|
| 156 | ig, &! loop on ngrid |
|---|
| 157 | l, &! loop on nlay |
|---|
| 158 | i, &! loop on nbinco2 |
|---|
| 159 | ! ---Variables for meteoritic flux |
|---|
| 160 | ibin, &! loop on nbin_meteor |
|---|
| 161 | idx_min,&! index of min(diff_pressure) |
|---|
| 162 | read_ok ! file uMeteor iostat |
|---|
| 163 | |
|---|
| 164 | real :: & |
|---|
| 165 | Nccnco2, & |
|---|
| 166 | Qccnco2, & |
|---|
| 167 | Nccnco2_h2o, & |
|---|
| 168 | Qccnco2_h2o, & |
|---|
| 169 | masse(ngrid,nlay), &! mass layer (kg.m-2) |
|---|
| 170 | rice(ngrid,nlay), &! water ice mass mean radius (m): used for nucleation of CO2 on ice-coated ccns |
|---|
| 171 | zq(ngrid,nlay,nq), &! local value of tracers (kg/kg) |
|---|
| 172 | zq0(ngrid,nlay,nq), &! local init value of tracers (kg/kg) |
|---|
| 173 | zt(ngrid,nlay), &! local value of temperature (K) |
|---|
| 174 | zqsat(ngrid,nlay), &! saturation vapor pressure for CO2 (K) |
|---|
| 175 | tcond(ngrid,nlay), &! condensation temperature of CO2 (K) |
|---|
| 176 | lw, &! Latent heat of sublimation (J.kg-1) |
|---|
| 177 | rdust(ngrid,nlay), &! Dust geometric mean radius (m) |
|---|
| 178 | rhocloud(ngrid,nlay), &! Cloud density (kg.m-3) |
|---|
| 179 | rhocloudco2(ngrid,nlay) ! Cloud density (kg.m-3) |
|---|
| 180 | |
|---|
| 181 | real(kind=8) :: & |
|---|
| 182 | derf ! Error function |
|---|
| 183 | |
|---|
| 184 | double precision :: & |
|---|
| 185 | dMice, &! mass of condensed ice |
|---|
| 186 | facteurmax, &! for energy limit on mass growth |
|---|
| 187 | pco2, &! Co2 vapor partial pressure (Pa) |
|---|
| 188 | satu, &! Co2 vapor saturation ratio over ice |
|---|
| 189 | Mo, &! mass of aerosol particles |
|---|
| 190 | No, &! number of aerosol particles |
|---|
| 191 | Rn, &! logarithm of rdust/rice |
|---|
| 192 | Rm, &! Rn * variance of ice and CCN distribution |
|---|
| 193 | n_derf, &! derf( (rb_cldco2(1)+Rn) *dev3) |
|---|
| 194 | m_derf, &! derf( (rb_cldco2(1)+Rm) *dev2) |
|---|
| 195 | n_aer(nbinco2_cld), &! Radius used by the microphysical scheme (m) |
|---|
| 196 | m_aer(nbinco2_cld), &! number concentration V-1 of particle/each size bin |
|---|
| 197 | n_aer_meteor(nbinco2_cld), &! Radius used by the microphysical scheme (m) |
|---|
| 198 | m_aer_meteor(nbinco2_cld), &! number concentration V-1 of particle/each size bin |
|---|
| 199 | n_aer_h2oice(nbinco2_cld), &! mass mixing ratio of particle/each size bin |
|---|
| 200 | m_aer_h2oice(nbinco2_cld), &! Same - for CO2 nucleation |
|---|
| 201 | Ic_rice, &! Mass transfer rate CO2 ice crystal |
|---|
| 202 | ratioh2o_ccn, &! 1./(zq(ig,l,igcm_h2o_ice) + zq(ig,l,igcm_ccn_mass)*tauscaling(ig)) |
|---|
| 203 | vo2co2, &! volume of co2 ice particle |
|---|
| 204 | dN, &! number of particle of dust used as ccn |
|---|
| 205 | dM, &! mass of dN |
|---|
| 206 | dN_meteor, &! number of particle of dust used as ccn |
|---|
| 207 | dM_meteor, &! mass of dN |
|---|
| 208 | dNh2o, &! number of particle of h2o ice used as ccn |
|---|
| 209 | dMh2o, &! mass of dNh2o |
|---|
| 210 | dNN, &! min(dN,zq(ig,l,igcm_dust_number)) |
|---|
| 211 | dMM, &! min(dM,zq(ig,l,igcm_dust_mass)) |
|---|
| 212 | dNN_meteor, &! number of particle of dust used as ccn |
|---|
| 213 | dMM_meteor, &! mass of dN |
|---|
| 214 | dNNh2o, &! min(dNNh2o,zq(ig,l,igcm_ccn_number)) |
|---|
| 215 | dMh2o_ice, &! min(dMh2o*zq(ig,l,igcm_h2o_ice)*ratioh2o_ccn, zq(ig,l,igcm_h2o_ice)) |
|---|
| 216 | dMh2o_ccn, &! min(dMh2o_ccn,zq(ig,l,igcm_ccn_mass)) |
|---|
| 217 | rate(nbinco2_cld), &! nucleation rate |
|---|
| 218 | rate_meteor(nbinco2_cld), &! nucleation rate |
|---|
| 219 | rateh2o(nbinco2_cld), &! nucleation rate for h2o |
|---|
| 220 | rho_ice_co2T, &! density of co2 ice Temperature-dependent |
|---|
| 221 | riceco2(ngrid,nlay), &! CO2 ice mean radius (m) |
|---|
| 222 | vrat_cld, &! Volume ratio |
|---|
| 223 | Proba, &! 1.0 - exp(-1.*microtimestep*rate(i)) |
|---|
| 224 | Probah2o, &! 1.0 - exp(-1.*microtimestep*rateh2o(i)) |
|---|
| 225 | Proba_meteor, &! 1.0 - exp(-1.*microtimestep*rateh2o(i)) |
|---|
| 226 | diff_pressure(nlev_meteor), &! abs(pression_meteor(:)-pplev(ig,l)) |
|---|
| 227 | meteor_ccn(ngrid,nlay,nbinco2_cld) ! |
|---|
| 228 | |
|---|
| 229 | logical :: & |
|---|
| 230 | file_ok ! test if meteoritic input file exists |
|---|
| 231 | !======================================================================================================================! |
|---|
| 232 | ! BEGIN ===============================================================================================================! |
|---|
| 233 | !======================================================================================================================! |
|---|
| 234 | ! 0. Firstcall |
|---|
| 235 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 236 | if (firstcall) then |
|---|
| 237 | firstcall = .false. |
|---|
| 238 | |
|---|
| 239 | ! Variance of the ice and CCN distributions |
|---|
| 240 | sigma_ice = sqrt(log(1.+nuice_sed)) |
|---|
| 241 | dev3 = 1. / ( sqrt(2.) * sigma_ice ) |
|---|
| 242 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 243 | ! 0.1. Bonus: meteoritic component, extract data |
|---|
| 244 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 245 | if (meteo_flux) then |
|---|
| 246 | ! Check if file exists |
|---|
| 247 | inquire(file=trim(datadir)//'/'//file_meteoritic_flux, exist=file_ok) |
|---|
| 248 | if (.not. file_ok) then |
|---|
| 249 | call abort_physic("CO2clouds", 'file '//file_meteoritic_flux//' should be in'//trim(datadir), 1) |
|---|
| 250 | end if |
|---|
| 251 | |
|---|
| 252 | ! open file |
|---|
| 253 | open(unit=uMeteor,file=trim(datadir)//'/'//file_meteoritic_flux, FORM='formatted') |
|---|
| 254 | |
|---|
| 255 | !skip 1 line |
|---|
| 256 | read(uMeteor,*) |
|---|
| 257 | |
|---|
| 258 | ! extract pressure_meteor |
|---|
| 259 | do i = 1, nlev_meteor |
|---|
| 260 | read(uMeteor,*)pression_meteor(i) |
|---|
| 261 | end do |
|---|
| 262 | |
|---|
| 263 | !skip 1 line |
|---|
| 264 | read(uMeteor,*) |
|---|
| 265 | |
|---|
| 266 | ! extract meteor flux |
|---|
| 267 | do i = 1, nlev_meteor |
|---|
| 268 | ! les mêmes 100 bins size que la distri nuclea : on touche pas |
|---|
| 269 | do ibin = 1, nbin_meteor |
|---|
| 270 | read(uMeteor,'(F12.6)')meteor(i,ibin) |
|---|
| 271 | end do |
|---|
| 272 | end do |
|---|
| 273 | |
|---|
| 274 | ! close file |
|---|
| 275 | close(uMeteor) |
|---|
| 276 | end if ! of if meteo_flux |
|---|
| 277 | end if ! firstcall |
|---|
| 278 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 279 | ! 1. Initialization |
|---|
| 280 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 281 | rdust(:,:) = 0. |
|---|
| 282 | meteor_ccn(:,:,:) = 0d0 |
|---|
| 283 | rice(:,:) = 1.e-8 |
|---|
| 284 | riceco2(:,:) = 1.e-11 |
|---|
| 285 | |
|---|
| 286 | ! Initialize the tendencies |
|---|
| 287 | subpdqcloudco2(:,:,:) = 0. |
|---|
| 288 | subpdtcloudco2(:,:) = 0. |
|---|
| 289 | |
|---|
| 290 | ! pteff temperature layer; sum_subpdt dT.s-1 |
|---|
| 291 | zt(1:ngrid,1:nlay) = 0. |
|---|
| 292 | zt(:,:) = pteff(:,:) + sum_subpdt(:,:) * microtimestep |
|---|
| 293 | |
|---|
| 294 | ! pqeff traceur kg/kg; sum_subpdq tendance idem .s-1 |
|---|
| 295 | zq(:,:,:) = pqeff(:,:,:) + sum_subpdq(:,:,:) * microtimestep |
|---|
| 296 | where( zq(:,:,:) < threshold ) zq(:,:,:) = threshold |
|---|
| 297 | |
|---|
| 298 | zq0(:,:,:) = zq(:,:,:) |
|---|
| 299 | zqsat(:,:) = 0. |
|---|
| 300 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 301 | ! 2. Compute saturation |
|---|
| 302 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 303 | call co2sat(ngrid*nlay,zt,zqsat) |
|---|
| 304 | call tcondco2(ngrid,nlay,pplay, zq(:,:,igcm_co2) + zq(:,:,igcm_co2_ice),tcond) |
|---|
| 305 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 306 | ! 3. Bonus: additional meteoritic particles for nucleation |
|---|
| 307 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 308 | if (meteo_flux) then |
|---|
| 309 | do l = 1, nlay |
|---|
| 310 | do ig = 1, ngrid |
|---|
| 311 | masse(ig,l) = (pplev(ig,l) - pplev(ig,l+1)) / g |
|---|
| 312 | |
|---|
| 313 | diff_pressure(:) = abs(pression_meteor(:)-pplev(ig,l)) |
|---|
| 314 | |
|---|
| 315 | idx_min = minloc(diff_pressure, DIM=1) |
|---|
| 316 | ! Par kg air csi par m carre, x epaisseur/masse pour par kg/air. Check original unit with J. Plane |
|---|
| 317 | meteor_ccn(ig,l,:) = meteor(idx_min,:)/masse(ig,l) |
|---|
| 318 | end do |
|---|
| 319 | end do |
|---|
| 320 | end if |
|---|
| 321 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 322 | ! 4. Actual microphysics: Main loop over the GCM's grid |
|---|
| 323 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 324 | do l = 1, nlay |
|---|
| 325 | do ig = 1, ngrid |
|---|
| 326 | |
|---|
| 327 | ! Get the partial pressure of co2 vapor and its saturation ratio |
|---|
| 328 | pco2 = zq(ig,l,igcm_co2) * (mmean(ig,l)/(mco2*1e3)) * pplay(ig,l) ! mco2 (kg/mol) => mmean (g/mol) |
|---|
| 329 | satu = pco2 / zqsat(ig,l) |
|---|
| 330 | |
|---|
| 331 | !T-dependant CO2 ice density |
|---|
| 332 | call density_co2_ice(dble(zt(ig,l)), rho_ice_co2T) |
|---|
| 333 | |
|---|
| 334 | vo2co2 = m0co2 / rho_ice_co2T |
|---|
| 335 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 336 | ! 4.1 Nucleation |
|---|
| 337 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 338 | ! if there is condensation |
|---|
| 339 | if ( satu >= 1 ) then |
|---|
| 340 | call updaterccnCO2(zq(ig,l,igcm_dust_mass), zq(ig,l,igcm_dust_number), rdust(ig,l), tauscaling(ig)) |
|---|
| 341 | |
|---|
| 342 | ! Expand the dust moments into a binned distribution |
|---|
| 343 | n_aer(:) = 0d0 ! number of aerosol particles |
|---|
| 344 | m_aer(:) = 0d0 ! mass of aerosol particles |
|---|
| 345 | |
|---|
| 346 | No = zq(ig,l,igcm_dust_number) * tauscaling(ig) |
|---|
| 347 | |
|---|
| 348 | Mo = (4./3.) * pi * rho_dust * No * rdust(ig,l)**3 *exp(9.*nuiceco2_ref/2.) |
|---|
| 349 | |
|---|
| 350 | if (No > threshold) then |
|---|
| 351 | Rn = -log(rdust(ig,l)) |
|---|
| 352 | |
|---|
| 353 | Rm = Rn - 3. * sigma_iceco2 * sigma_iceco2 |
|---|
| 354 | |
|---|
| 355 | n_derf = derf( (rb_cldco2(1)+Rn) *dev2) |
|---|
| 356 | m_derf = derf( (rb_cldco2(1)+Rm) *dev2) |
|---|
| 357 | |
|---|
| 358 | do i = 1, nbinco2_cld |
|---|
| 359 | n_aer(i) = -0.5 * No * n_derf |
|---|
| 360 | m_aer(i) = -0.5 * Mo * m_derf |
|---|
| 361 | |
|---|
| 362 | n_derf = derf((rb_cldco2(i+1)+Rn) *dev2) |
|---|
| 363 | m_derf = derf((rb_cldco2(i+1)+Rm) *dev2) |
|---|
| 364 | |
|---|
| 365 | n_aer(i) = n_aer(i) + 0.5 * No * n_derf |
|---|
| 366 | m_aer(i) = m_aer(i) + 0.5 * Mo * m_derf |
|---|
| 367 | end do |
|---|
| 368 | |
|---|
| 369 | ! Call to nucleation routine |
|---|
| 370 | call nucleaco2(dble(pco2), zt(ig,l), dble(satu), n_aer, rate, vo2co2, mtetaco2) |
|---|
| 371 | dN = 0. |
|---|
| 372 | dM = 0. |
|---|
| 373 | |
|---|
| 374 | do i = 1, nbinco2_cld |
|---|
| 375 | Proba = 1.0 - exp(-1.*microtimestep*rate(i)) |
|---|
| 376 | dN = dN + n_aer(i) * Proba |
|---|
| 377 | dM = dM + m_aer(i) * Proba |
|---|
| 378 | end do |
|---|
| 379 | |
|---|
| 380 | ! Now increment CCN tracers and update dust tracers |
|---|
| 381 | dNN = min(dN, zq(ig,l,igcm_dust_number)) ! dNN est devenu DN |
|---|
| 382 | dMM = min(dM, zq(ig,l,igcm_dust_mass)) ! idem dans le min |
|---|
| 383 | |
|---|
| 384 | zq(ig,l,igcm_ccnco2_mass) = zq(ig,l,igcm_ccnco2_mass) + dMM /tauscaling(ig) |
|---|
| 385 | zq(ig,l,igcm_ccnco2_number) = zq(ig,l,igcm_ccnco2_number) + dNN /tauscaling(ig) |
|---|
| 386 | |
|---|
| 387 | zq(ig,l,igcm_dust_mass) = zq(ig,l,igcm_dust_mass) - dMM /tauscaling(ig) |
|---|
| 388 | zq(ig,l,igcm_dust_number) = zq(ig,l,igcm_dust_number) - dNN /tauscaling(ig) |
|---|
| 389 | end if ! of if No > 1e-30 |
|---|
| 390 | |
|---|
| 391 | ! Ajout meteor_ccn particles aux particules de poussière background |
|---|
| 392 | if (meteo_flux) then |
|---|
| 393 | n_aer_meteor(1:nbinco2_cld) = 0d0 |
|---|
| 394 | m_aer_meteor(1:nbinco2_cld) = 0d0 |
|---|
| 395 | do i = 1, nbinco2_cld |
|---|
| 396 | n_aer_meteor(i) = meteor_ccn(ig,l,i) |
|---|
| 397 | m_aer_meteor(i) = (4./3.) * pi * rho_dust * meteor_ccn(ig,l,i) * rad_cldco2(i)**3 |
|---|
| 398 | end do |
|---|
| 399 | ! Call to nucleation routine |
|---|
| 400 | rate_meteor(1:nbinco2_cld) = 0d0 |
|---|
| 401 | call nucleaco2(dble(pco2), zt(ig,l), dble(satu), n_aer_meteor, rate_meteor, vo2co2, mtetaco2) |
|---|
| 402 | |
|---|
| 403 | dN_meteor = 0. |
|---|
| 404 | dM_meteor = 0. |
|---|
| 405 | do i = 1, nbinco2_cld |
|---|
| 406 | Proba_meteor = 1.0 - exp(-1.*microtimestep*rate_meteor(i)) |
|---|
| 407 | dN_meteor = dN_meteor + n_aer_meteor(i) * Proba_meteor |
|---|
| 408 | dM_meteor = dM_meteor + m_aer_meteor(i) * Proba_meteor |
|---|
| 409 | end do |
|---|
| 410 | ! Now increment CCN tracers and update dust tracers |
|---|
| 411 | zq(ig,l,igcm_ccnco2_meteor_mass) = zq(ig,l,igcm_ccnco2_meteor_mass) + dM_meteor |
|---|
| 412 | zq(ig,l,igcm_ccnco2_meteor_number) = zq(ig,l,igcm_ccnco2_meteor_number) + dN_meteor |
|---|
| 413 | |
|---|
| 414 | zq(ig,l,igcm_ccnco2_mass) = zq(ig,l,igcm_ccnco2_mass) + dM_meteor |
|---|
| 415 | zq(ig,l,igcm_ccnco2_number) = zq(ig,l,igcm_ccnco2_number) + dN_meteor |
|---|
| 416 | end if |
|---|
| 417 | |
|---|
| 418 | ! Same but with h2o particles as CCN only if co2useh2o = .true. |
|---|
| 419 | if (co2useh2o) then |
|---|
| 420 | call updaterice_micro(zq(ig,l,igcm_h2o_ice), zq(ig,l,igcm_ccn_mass), zq(ig,l,igcm_ccn_number), & |
|---|
| 421 | tauscaling(ig), rice(ig,l), rhocloud(ig,l)) |
|---|
| 422 | |
|---|
| 423 | ! Total mass of H20 crystals,CCN included |
|---|
| 424 | Mo = zq(ig,l,igcm_h2o_ice) + zq(ig,l,igcm_ccn_mass) * tauscaling(ig) + threshold |
|---|
| 425 | |
|---|
| 426 | No = zq(ig,l,igcm_ccn_number) * tauscaling(ig) + threshold |
|---|
| 427 | if (No > threshold) then |
|---|
| 428 | Rn = -log(rice(ig,l)) |
|---|
| 429 | |
|---|
| 430 | Rm = Rn - 3. * sigma_ice * sigma_ice |
|---|
| 431 | |
|---|
| 432 | n_derf = derf( (rb_cldco2(1)+Rn) *dev3) |
|---|
| 433 | m_derf = derf( (rb_cldco2(1)+Rm) *dev3) |
|---|
| 434 | |
|---|
| 435 | n_aer_h2oice(:)=0. |
|---|
| 436 | m_aer_h2oice(:)=0. |
|---|
| 437 | do i = 1, nbinco2_cld |
|---|
| 438 | n_aer_h2oice(i) = -0.5 * No * n_derf |
|---|
| 439 | m_aer_h2oice(i) = -0.5 * Mo * m_derf |
|---|
| 440 | |
|---|
| 441 | n_derf = derf( (rb_cldco2(i+1)+Rn) *dev3) |
|---|
| 442 | m_derf = derf( (rb_cldco2(i+1)+Rm) *dev3) |
|---|
| 443 | |
|---|
| 444 | n_aer_h2oice(i) = n_aer_h2oice(i) + 0.5 * No * n_derf |
|---|
| 445 | m_aer_h2oice(i) = m_aer_h2oice(i) + 0.5 * Mo * m_derf |
|---|
| 446 | end do |
|---|
| 447 | |
|---|
| 448 | call nucleaco2(dble(pco2), zt(ig,l), dble(satu), n_aer_h2oice, rateh2o, vo2co2, mteta) |
|---|
| 449 | dNh2o = 0. |
|---|
| 450 | dMh2o = 0. |
|---|
| 451 | |
|---|
| 452 | do i = 1, nbinco2_cld |
|---|
| 453 | Probah2o = 1.0 - exp(-1.*microtimestep*rateh2o(i)) |
|---|
| 454 | dNh2o = dNh2o + n_aer_h2oice(i) * Probah2o |
|---|
| 455 | dMh2o = dMh2o + m_aer_h2oice(i) * Probah2o |
|---|
| 456 | end do |
|---|
| 457 | |
|---|
| 458 | ! Update CCN for CO2 nucleating on H2O CCN : Warning: must keep memory of it |
|---|
| 459 | dNNh2o = dNh2o/tauscaling(ig) |
|---|
| 460 | dNNh2o = min(dNNh2o,zq(ig,l,igcm_ccn_number)) |
|---|
| 461 | |
|---|
| 462 | ratioh2o_ccn = 1./(zq(ig,l,igcm_h2o_ice) + zq(ig,l,igcm_ccn_mass)*tauscaling(ig)) |
|---|
| 463 | |
|---|
| 464 | dMh2o_ccn = dMh2o * zq(ig,l,igcm_ccn_mass) * tauscaling(ig) * ratioh2o_ccn |
|---|
| 465 | dMh2o_ccn = dMh2o_ccn/tauscaling(ig) |
|---|
| 466 | dMh2o_ccn = min(dMh2o_ccn,zq(ig,l,igcm_ccn_mass)) |
|---|
| 467 | |
|---|
| 468 | dMh2o_ice = dMh2o * zq(ig,l,igcm_h2o_ice) * ratioh2o_ccn |
|---|
| 469 | dMh2o_ice = min(dMh2o_ice,zq(ig,l,igcm_h2o_ice)) |
|---|
| 470 | |
|---|
| 471 | zq(ig,l,igcm_ccnco2_h2o_mass_ice) = zq(ig,l,igcm_ccnco2_h2o_mass_ice) + dMh2o_ice |
|---|
| 472 | zq(ig,l,igcm_ccnco2_h2o_mass_ccn) = zq(ig,l,igcm_ccnco2_h2o_mass_ccn) + dMh2o_ccn |
|---|
| 473 | zq(ig,l,igcm_ccnco2_h2o_number) = zq(ig,l,igcm_ccnco2_h2o_number) + dNNh2o |
|---|
| 474 | |
|---|
| 475 | zq(ig,l,igcm_ccn_number) = zq(ig,l,igcm_ccn_number) - dNNh2o |
|---|
| 476 | zq(ig,l,igcm_h2o_ice) = zq(ig,l,igcm_h2o_ice) - dMh2o_ice |
|---|
| 477 | zq(ig,l,igcm_ccn_mass) = zq(ig,l,igcm_ccn_mass) - dMh2o_ccn |
|---|
| 478 | |
|---|
| 479 | zq(ig,l,igcm_ccnco2_mass) = zq(ig,l,igcm_ccnco2_mass) + dMh2o_ice + dMh2o_ccn |
|---|
| 480 | zq(ig,l,igcm_ccnco2_number) = zq(ig,l,igcm_ccnco2_number) + dNNh2o |
|---|
| 481 | end if |
|---|
| 482 | end if ! of if co2useh2o |
|---|
| 483 | end if ! of is satu > 1 |
|---|
| 484 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 485 | ! 4.2. Ice growth: scheme for radius evolution |
|---|
| 486 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 487 | ! We trigger crystal growth if and only if there is at least one nuclei (N>1). Indeed, if we are supersaturated |
|---|
| 488 | ! and still don't have at least one nuclei, we should better wait to avoid unrealistic value for nuclei radius |
|---|
| 489 | ! and so on for cases that remain negligible. |
|---|
| 490 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 491 | ! we trigger crystal growth |
|---|
| 492 | if ((zq(ig,l,igcm_ccnco2_number)) * tauscaling(ig) + threshold >= 1) then |
|---|
| 493 | Nccnco2 = dble(zq(ig,l,igcm_ccnco2_number)) |
|---|
| 494 | Qccnco2 = dble(zq(ig,l,igcm_ccnco2_mass)) |
|---|
| 495 | Nccnco2_h2o = 0. |
|---|
| 496 | Qccnco2_h2o = 0. |
|---|
| 497 | |
|---|
| 498 | if (co2useh2o) then |
|---|
| 499 | Nccnco2_h2o = zq(ig,l,igcm_ccnco2_h2o_number) |
|---|
| 500 | Qccnco2_h2o = zq(ig,l,igcm_ccnco2_h2o_mass_ice) + zq(ig,l,igcm_ccnco2_h2o_mass_ccn) |
|---|
| 501 | Nccnco2 = Nccnco2 - Nccnco2_h2o |
|---|
| 502 | Qccnco2 = Qccnco2 - Qccnco2_h2o |
|---|
| 503 | end if |
|---|
| 504 | |
|---|
| 505 | call updaterice_microco2(dble(zq(ig,l,igcm_co2_ice)), dble(Qccnco2), dble(Nccnco2), dble(Qccnco2_h2o), & |
|---|
| 506 | dble(Nccnco2_h2o), zt(ig,l), tauscaling(ig), riceco2(ig,l), rhocloudco2(ig,l)) |
|---|
| 507 | Ic_rice = 0. |
|---|
| 508 | |
|---|
| 509 | ! J.kg-1 |
|---|
| 510 | lw = l0 + l1 * zt(ig,l) + l2 * zt(ig,l)**2 + l3 * zt(ig,l)**3 + l4 * zt(ig,l)**4 |
|---|
| 511 | |
|---|
| 512 | facteurmax = abs(tcond(ig,l)-zt(ig,l)) * (cpp/lw) |
|---|
| 513 | |
|---|
| 514 | ! call scheme of microphys. mass growth for CO2 (evaporation/condensation) |
|---|
| 515 | call massflowrateCO2(pplay(ig,l), zt(ig,l), satu, riceco2(ig,l), mmean(ig,l), Ic_rice) |
|---|
| 516 | |
|---|
| 517 | ! Ic_rice Mass transfer rate (kg/s) for a rice particle > 0 si croissance ! |
|---|
| 518 | if (isnan(Ic_rice) .or. Ic_rice == 0.) then |
|---|
| 519 | Ic_rice = 0. |
|---|
| 520 | subpdtcloudco2(ig,l) = -sum_subpdt(ig,l) |
|---|
| 521 | dMice = 0 |
|---|
| 522 | else |
|---|
| 523 | ! Kg par kg d'air, >0 si croissance ! |
|---|
| 524 | ! kg.s-1 par particule * nb particule par kg air*s = kg par kg air |
|---|
| 525 | dMice = zq(ig,l,igcm_ccnco2_number) * Ic_rice * microtimestep * tauscaling(ig) |
|---|
| 526 | |
|---|
| 527 | ! facteurmax maximum quantity of CO2 that can sublime/condense according to available thermal energy |
|---|
| 528 | ! latent heat release > 0 if growth i.e. if dMice > 0 |
|---|
| 529 | dMice = max(dMice,max(-facteurmax,-zq(ig,l,igcm_co2_ice))) |
|---|
| 530 | dMice = min(dMice,min(facteurmax,zq(ig,l,igcm_co2))) |
|---|
| 531 | |
|---|
| 532 | ! kgco2/kgair* J/kgco2 * 1/(J.kgair-1.K-1)/s = K /s |
|---|
| 533 | subpdtcloudco2(ig,l) = dMice * lw / cpp / microtimestep |
|---|
| 534 | |
|---|
| 535 | !Now update tracers |
|---|
| 536 | zq(ig,l,igcm_co2_ice) = zq(ig,l,igcm_co2_ice) + dMice |
|---|
| 537 | zq(ig,l,igcm_co2) = zq(ig,l,igcm_co2) - dMice |
|---|
| 538 | end if |
|---|
| 539 | end if ! if zq(ccnco2_number) >= 1 |
|---|
| 540 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 541 | ! 4.3 Dust cores releasing if no more co2 ice |
|---|
| 542 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 543 | ! On sublime tout |
|---|
| 544 | if ((zq(ig,l,igcm_co2_ice) <= threshold).or.(zq(ig,l,igcm_ccnco2_number) * tauscaling(ig) < 1.)) then |
|---|
| 545 | zq(ig,l,igcm_dust_mass) = zq(ig,l,igcm_dust_mass) + zq(ig,l,igcm_ccnco2_mass) |
|---|
| 546 | zq(ig,l,igcm_dust_number) = zq(ig,l,igcm_dust_number) + zq(ig,l,igcm_ccnco2_number) |
|---|
| 547 | zq(ig,l,igcm_ccnco2_mass) = 0. |
|---|
| 548 | zq(ig,l,igcm_ccnco2_number) = 0. |
|---|
| 549 | if (meteo_flux) then |
|---|
| 550 | zq(ig,l,igcm_dust_mass) = zq(ig,l,igcm_dust_mass) - zq(ig,l,igcm_ccnco2_meteor_mass) |
|---|
| 551 | zq(ig,l,igcm_dust_number) = zq(ig,l,igcm_dust_number) - zq(ig,l,igcm_ccnco2_meteor_number) |
|---|
| 552 | zq(ig,l,igcm_ccnco2_meteor_mass) = 0. |
|---|
| 553 | zq(ig,l,igcm_ccnco2_meteor_number) = 0. |
|---|
| 554 | end if |
|---|
| 555 | if (co2useh2o) then |
|---|
| 556 | zq(ig,l,igcm_dust_mass) = zq(ig,l,igcm_dust_mass) - zq(ig,l,igcm_ccnco2_h2o_mass_ccn) - & |
|---|
| 557 | zq(ig,l,igcm_ccnco2_h2o_mass_ice) |
|---|
| 558 | zq(ig,l,igcm_dust_number) = zq(ig,l,igcm_dust_number) - zq(ig,l,igcm_ccnco2_h2o_number) |
|---|
| 559 | |
|---|
| 560 | zq(ig,l,igcm_ccn_mass) = zq(ig,l,igcm_ccn_mass) + zq(ig,l,igcm_ccnco2_h2o_mass_ccn) |
|---|
| 561 | zq(ig,l,igcm_h2o_ice) = zq(ig,l,igcm_h2o_ice) + zq(ig,l,igcm_ccnco2_h2o_mass_ice) |
|---|
| 562 | zq(ig,l,igcm_ccn_number) = zq(ig,l,igcm_ccn_number) + zq(ig,l,igcm_ccnco2_h2o_number) |
|---|
| 563 | zq(ig,l,igcm_ccnco2_h2o_number) = 0. |
|---|
| 564 | zq(ig,l,igcm_ccnco2_h2o_mass_ccn) = 0. |
|---|
| 565 | zq(ig,l,igcm_ccnco2_h2o_mass_ice) = 0. |
|---|
| 566 | end if |
|---|
| 567 | zq(ig,l,igcm_co2) = zq(ig,l,igcm_co2) + zq(ig,l,igcm_co2_ice) |
|---|
| 568 | zq(ig,l,igcm_co2_ice) = 0. |
|---|
| 569 | riceco2(ig,l) = 0. |
|---|
| 570 | end if !of if co2_ice < threshold or zq(ccnco2_number) < 1 |
|---|
| 571 | end do ! of ig loop |
|---|
| 572 | end do ! of nlayer loop |
|---|
| 573 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 574 | ! 5. Get cloud tendencies |
|---|
| 575 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 576 | subpdqcloudco2(:,:,igcm_co2) = ( zq(:,:,igcm_co2) - zq0(:,:,igcm_co2) ) / microtimestep |
|---|
| 577 | |
|---|
| 578 | subpdqcloudco2(:,:,igcm_co2_ice) = ( zq(:,:,igcm_co2_ice) - zq0(:,:,igcm_co2_ice) ) / microtimestep |
|---|
| 579 | |
|---|
| 580 | subpdqcloudco2(:,:,igcm_ccnco2_mass) = ( zq(:,:,igcm_ccnco2_mass) - zq0(:,:,igcm_ccnco2_mass) ) / microtimestep |
|---|
| 581 | |
|---|
| 582 | subpdqcloudco2(:,:,igcm_ccnco2_number) = ( zq(:,:,igcm_ccnco2_number) - zq0(:,:,igcm_ccnco2_number))/microtimestep |
|---|
| 583 | |
|---|
| 584 | subpdqcloudco2(:,:,igcm_dust_mass) = ( zq(:,:,igcm_dust_mass) - zq0(:,:,igcm_dust_mass) ) / microtimestep |
|---|
| 585 | |
|---|
| 586 | subpdqcloudco2(:,:,igcm_dust_number) = ( zq(:,:,igcm_dust_number) - zq0(:,:,igcm_dust_number) ) / microtimestep |
|---|
| 587 | |
|---|
| 588 | if (meteo_flux) then |
|---|
| 589 | subpdqcloudco2(:,:,igcm_ccnco2_meteor_mass) = (zq(:,:,igcm_ccnco2_meteor_mass)-zq0(:,:,igcm_ccnco2_meteor_mass) & |
|---|
| 590 | )/microtimestep |
|---|
| 591 | |
|---|
| 592 | subpdqcloudco2(:,:,igcm_ccnco2_meteor_number) = ( zq(:,:,igcm_ccnco2_meteor_number) - & |
|---|
| 593 | zq0(:,:,igcm_ccnco2_meteor_number) )/microtimestep |
|---|
| 594 | end if |
|---|
| 595 | |
|---|
| 596 | if (co2useh2o) then |
|---|
| 597 | subpdqcloudco2(:,:,igcm_h2o_ice) = ( zq(:,:,igcm_h2o_ice) - zq0(:,:,igcm_h2o_ice) ) / microtimestep |
|---|
| 598 | |
|---|
| 599 | subpdqcloudco2(:,:,igcm_ccn_mass) = ( zq(:,:,igcm_ccn_mass) - zq0(:,:,igcm_ccn_mass) ) / microtimestep |
|---|
| 600 | |
|---|
| 601 | subpdqcloudco2(:,:,igcm_ccn_number) = ( zq(:,:,igcm_ccn_number) - zq0(:,:,igcm_ccn_number) ) / microtimestep |
|---|
| 602 | |
|---|
| 603 | subpdqcloudco2(:,:,igcm_ccnco2_h2o_mass_ice) = (zq(:,:,igcm_ccnco2_h2o_mass_ice)-zq0(:,:,igcm_ccnco2_h2o_mass_ice)& |
|---|
| 604 | )/microtimestep |
|---|
| 605 | |
|---|
| 606 | subpdqcloudco2(:,:,igcm_ccnco2_h2o_mass_ccn) = (zq(:,:,igcm_ccnco2_h2o_mass_ccn)-zq0(:,:,igcm_ccnco2_h2o_mass_ccn)& |
|---|
| 607 | )/microtimestep |
|---|
| 608 | |
|---|
| 609 | subpdqcloudco2(:,:,igcm_ccnco2_h2o_number) = ( zq(:,:,igcm_ccnco2_h2o_number) - zq0(:,:,igcm_ccnco2_h2o_number) & |
|---|
| 610 | )/microtimestep |
|---|
| 611 | end if |
|---|
| 612 | !======================================================================================================================! |
|---|
| 613 | ! END =================================================================================================================! |
|---|
| 614 | !======================================================================================================================! |
|---|
| 615 | end subroutine improvedCO2clouds |
|---|
| 616 | |
|---|
| 617 | end module improvedCO2clouds_mod |
|---|
| 618 | |
|---|
