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