| 1 | !======================================================================================================================! |
|---|
| 2 | ! Module: CO2 clouds formation ========================================================================================! |
|---|
| 3 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 4 | ! Authors: Joachim Audouard, Constantino Listowski, Anni Määttänen |
|---|
| 5 | ! Date: 09/2016 |
|---|
| 6 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 7 | ! Contains subroutines: |
|---|
| 8 | ! - co2cloud: of co2 cloud microphysics |
|---|
| 9 | !======================================================================================================================! |
|---|
| 10 | module co2cloud_mod |
|---|
| 11 | |
|---|
| 12 | implicit none |
|---|
| 13 | |
|---|
| 14 | contains |
|---|
| 15 | !======================================================================================================================! |
|---|
| 16 | ! SUBROUTINE: co2cloud ================================================================================================! |
|---|
| 17 | !======================================================================================================================! |
|---|
| 18 | ! Subject: |
|---|
| 19 | !--------- |
|---|
| 20 | ! Main subroutine of co2 cloud microphysics |
|---|
| 21 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 22 | ! Comments: |
|---|
| 23 | !---------- |
|---|
| 24 | ! - Adaptation of the water ice clouds scheme (with specific microphysics) of Montmessin, Navarro et al. |
|---|
| 25 | ! |
|---|
| 26 | ! - Microphysics subroutine is improvedCO2clouds.F |
|---|
| 27 | ! |
|---|
| 28 | ! - There is a time loop specific to cloud formation due to timescales smaller than the GCM integration timestep |
|---|
| 29 | ! |
|---|
| 30 | ! - The microphysics time step is a fraction of the physical one |
|---|
| 31 | ! |
|---|
| 32 | ! - The CO2 clouds tracers (co2_ice, ccn mass and concentration) are sedimented at each microtimestep. pdqs_sedco2 |
|---|
| 33 | ! keeps track of the CO2 flux at the surface |
|---|
| 34 | ! |
|---|
| 35 | ! - The subgrid Temperature distribution is modulated (0 or 1) by Spiga et al. (GRL 2012) |
|---|
| 36 | ! |
|---|
| 37 | ! - Saturation Index to account for GW propagation or dissipation upwards |
|---|
| 38 | ! |
|---|
| 39 | ! - 4D and column opacities are computed using Qext values at 1µm |
|---|
| 40 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 41 | ! Papers: |
|---|
| 42 | !-------- |
|---|
| 43 | ! "Near-pure vapor condensation in the Martian atmosphere: CO2 ice crystal growth", Listowski et al. (2013), JGRE |
|---|
| 44 | ! "Modeling the microphysics of CO2 ice clouds within wave-induced cold pockets in the martian mesosphere", Listowski |
|---|
| 45 | ! et al. (2014), Icarus |
|---|
| 46 | ! "Global climate modeling of the Martian water cycle with improved microphysics and radiatively active water ice |
|---|
| 47 | ! clouds", Navarro et al. (2014), JGRE |
|---|
| 48 | ! "Martian GCM with complete CO2 clouds microphysics", Audouard et al. (2017), EPSC abstract |
|---|
| 49 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 50 | ! Algorithm: |
|---|
| 51 | !----------- |
|---|
| 52 | ! 0. Firstcall |
|---|
| 53 | ! 0.1. Initialization of microtimestep from imicroco2 |
|---|
| 54 | ! 0.2. Compute the radius grid of CO2 ice particles (rb_cldco2) |
|---|
| 55 | ! 0.3. Read file 'optprop_co2ice_1mic.dat' to extract optical properties of CO2 ice at 1 micron (Qext) |
|---|
| 56 | ! 0.4. Interpole the radius grid (rb_cldco2) to get the corresponding exctinction coefficients (Qext1bins) |
|---|
| 57 | ! 0.5. Save the radius grid of CO2 particles (rb_cldco2) |
|---|
| 58 | ! 1. Initialization |
|---|
| 59 | ! 2. Compute mass and thickness layers |
|---|
| 60 | ! 3. Define the sub-grid cloud (CLFvaryingCO2) |
|---|
| 61 | ! 3.1. Representation of sub-grid CO2 ice clouds (CLFvaryingCO2 = True) |
|---|
| 62 | ! 3.1.a. Saturation index CO2 |
|---|
| 63 | ! 3.1.b. Compute tcond |
|---|
| 64 | ! 3.1.c. Compute cloud fraction in cells |
|---|
| 65 | ! 3.2. No sub-grid cloud representation (CLFvaryingCO2 = False) |
|---|
| 66 | ! 4. Microphysics of CO2 cloud formation |
|---|
| 67 | ! 4.1. Effective tracers quantities in the cloud |
|---|
| 68 | ! 4.2. Stepped entry for tendancies: At each micro timestep we add pdt in order to have a stepped entry |
|---|
| 69 | ! 4.3. Main call to the cloud scheme |
|---|
| 70 | ! 4.4. Updating tendencies after cloud scheme |
|---|
| 71 | ! 4.5. Gravitational sedimentation |
|---|
| 72 | ! 4.5.a. Compute cloud density |
|---|
| 73 | ! 4.5.b. Save actualized tracer values to compute sedimentation tendancies |
|---|
| 74 | ! 4.5.c. Sedimentation of co2 ice |
|---|
| 75 | ! 4.5.d. Sedimentation for other tracers |
|---|
| 76 | ! 4.5.e. Compute tendencies due to the sedimation process |
|---|
| 77 | ! 5. Compute final tendencies after time loop |
|---|
| 78 | ! 6. Update clouds physical values in the cloud (for output) |
|---|
| 79 | ! 6.1. Update density of co2 ice, riceco2 and opacity |
|---|
| 80 | ! 6.2. Update rice and rdust |
|---|
| 81 | ! 7. Correction if a lot of subliming CO2 fills the 1st layer |
|---|
| 82 | ! 8. CO2 saturated quantities |
|---|
| 83 | ! 8.1 Compute CO2 saturation in layers |
|---|
| 84 | ! 8.2 Compute CO2 saturated quantities in layers |
|---|
| 85 | ! 9. Everything modified by CO2 microphysics must be wrt co2cloudfrac |
|---|
| 86 | ! 10. Compute opacity at 1 micron |
|---|
| 87 | ! 11. Write outputs in diagfi.nc |
|---|
| 88 | !======================================================================================================================! |
|---|
| 89 | subroutine co2cloud(ngrid, nlay, ptimestep, pplev, pplay, pdpsrf, pzlay, pt, pdt, pq, pdq, pdqcloudco2, pdtcloudco2, & |
|---|
| 90 | nq, tau, tauscaling, rdust, rice, riceco2, nuice, rhocloud, rsedcloudco2, rhocloudco2, pzlev,& |
|---|
| 91 | pdqs_sedco2, pdqs_sedccn, pdu, pu, pcondicea) |
|---|
| 92 | |
|---|
| 93 | use ioipsl_getincom, only: getin |
|---|
| 94 | use dimradmars_mod, only: naerkind |
|---|
| 95 | use comcstfi_h, only: pi, g, cpp |
|---|
| 96 | use updaterad, only: updaterice_microco2, updaterice_micro, updaterdust |
|---|
| 97 | use conc_mod, only: mmean, rnew |
|---|
| 98 | use tracer_mod, only: igcm_co2, igcm_co2_ice, igcm_dust_mass, igcm_dust_number, igcm_h2o_ice, & |
|---|
| 99 | igcm_ccn_mass, igcm_ccn_number, igcm_ccnco2_mass, igcm_ccnco2_number, & |
|---|
| 100 | igcm_ccnco2_h2o_number, igcm_ccnco2_h2o_mass_ice, igcm_ccnco2_h2o_mass_ccn, rho_dust, & |
|---|
| 101 | nuiceco2_sed, nuiceco2_ref, r3n_q, rho_ice, nuice_sed, igcm_ccnco2_meteor_mass, & |
|---|
| 102 | igcm_ccnco2_meteor_number |
|---|
| 103 | |
|---|
| 104 | use newsedim_mod, only: newsedim |
|---|
| 105 | |
|---|
| 106 | use datafile_mod, only: datadir |
|---|
| 107 | use density_co2_ice_mod, only: density_co2_ice |
|---|
| 108 | use tcondco2_mod, only: tcondco2 |
|---|
| 109 | use improvedCO2clouds_mod, only: improvedCO2clouds |
|---|
| 110 | use microphys_h, only: nbinco2_cld, rad_cldco2, mco2 |
|---|
| 111 | #ifndef MESOSCALE |
|---|
| 112 | use write_output_mod, only: write_output |
|---|
| 113 | #endif |
|---|
| 114 | use callkeys_mod, only: sedimentation, CLFvaryingCO2, co2useh2o, & |
|---|
| 115 | satindexco2, meteo_flux |
|---|
| 116 | #ifndef MESOSCALE |
|---|
| 117 | use vertical_layers_mod, only: ap, bp |
|---|
| 118 | #endif |
|---|
| 119 | |
|---|
| 120 | implicit none |
|---|
| 121 | |
|---|
| 122 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 123 | ! VARIABLES DECLARATION |
|---|
| 124 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 125 | ! Input arguments: |
|---|
| 126 | !----------------- |
|---|
| 127 | integer, intent(in) ::& |
|---|
| 128 | ngrid, &! Number of grid points |
|---|
| 129 | nlay, &! Number of layers |
|---|
| 130 | nq ! Number of tracers |
|---|
| 131 | |
|---|
| 132 | real, intent(in) :: & |
|---|
| 133 | ptimestep, &! Physical time step (s) |
|---|
| 134 | pplev(ngrid,nlay+1), &! Inter-layer pressures (Pa) |
|---|
| 135 | pplay(ngrid,nlay), &! Mid-layer pressures (Pa) |
|---|
| 136 | pdpsrf(ngrid), &! Tendency on surface pressure |
|---|
| 137 | pzlay(ngrid,nlay), &! Altitude at the middle of the layers |
|---|
| 138 | pt(ngrid,nlay), &! Temperature at the middle of the layers (K) |
|---|
| 139 | pdt(ngrid,nlay), &! Tendency on temperature from other parametrizations |
|---|
| 140 | pq(ngrid,nlay,nq), &! Tracers (kg/kg) |
|---|
| 141 | pdq(ngrid,nlay,nq), &! Tendencies before condensation (kg/kg.s-1) |
|---|
| 142 | tau(ngrid,naerkind), &! Column dust optical depth at each point |
|---|
| 143 | tauscaling(ngrid), &! Convertion factor for dust amount |
|---|
| 144 | pu(ngrid,nlay), &! Zonal Wind: zu = pu + (pdu * ptimestep) |
|---|
| 145 | pdu(ngrid,nlay), &! Tendency of zonal wind before condensation |
|---|
| 146 | pzlev(ngrid,nlay+1), &! Altitude at the boundaries of the layers |
|---|
| 147 | nuice(ngrid,nlay) ! Estimated effective variance of the size distribution |
|---|
| 148 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 149 | ! Output arguments: |
|---|
| 150 | !------------------ |
|---|
| 151 | real, intent(out) :: & |
|---|
| 152 | rice(ngrid,nlay), & ! Water Ice mass mean radius (m) |
|---|
| 153 | ! rsedcloud(ngrid,nlay), & ! Water Cloud sedimentation radius |
|---|
| 154 | rhocloud(ngrid,nlay), & ! Water Cloud density (kg.m-3) |
|---|
| 155 | pdqs_sedco2(ngrid), & ! CO2 flux at the surface |
|---|
| 156 | pdqs_sedccn(ngrid,nq), & ! CCN flux at the surface |
|---|
| 157 | pdqcloudco2(ngrid,nlay,nq),& ! Tendency due to CO2 condensation (kg/kg.s-1) |
|---|
| 158 | pcondicea(ngrid,nlay), & ! Rate of condensation/sublimation of co2 ice in layers |
|---|
| 159 | pdtcloudco2(ngrid,nlay), & ! Tendency on temperature due to latent heat |
|---|
| 160 | rsedcloudco2(ngrid,nlay) ! Cloud sedimentation radius |
|---|
| 161 | |
|---|
| 162 | real, intent(inout) :: & |
|---|
| 163 | rdust(ngrid,nlay) ! Dust geometric mean radius (m) |
|---|
| 164 | |
|---|
| 165 | double precision, intent(out) :: & |
|---|
| 166 | riceco2(ngrid,nlay) ! Ice mass mean radius (m) r_c in Montmessin et al. (2004) |
|---|
| 167 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 168 | ! Local: |
|---|
| 169 | !------- |
|---|
| 170 | !-----1) Parameters: |
|---|
| 171 | !------------------- |
|---|
| 172 | integer, parameter :: & |
|---|
| 173 | uQext = 555, &! file_qext unit ID |
|---|
| 174 | var_dim_qext = 10000 ! Exact dimension of radv and qextv1mic from file_qext |
|---|
| 175 | |
|---|
| 176 | real, parameter :: & |
|---|
| 177 | mincloud = 0.1, &! Minimum cloud fraction |
|---|
| 178 | beta = 0.85, &! correction for the shape of the particles (see Murphy et al. JGR 1990 vol.95): |
|---|
| 179 | ! beta = 1 for spheres |
|---|
| 180 | ! beta = 0.85 for irregular particles |
|---|
| 181 | ! beta = 0.5 for disk shaped particles |
|---|
| 182 | threshold = 1e-30, & ! limit value considering as zero |
|---|
| 183 | threshold_2 = 1e-13 ! limit value considering the value is physical (below this value => computer noise for dble) |
|---|
| 184 | |
|---|
| 185 | double precision, parameter :: & |
|---|
| 186 | rmin_cld = 1.e-9, &! Minimum radius (m) |
|---|
| 187 | rmax_cld = 5.e-6, &! Maximum radius (m) |
|---|
| 188 | rbmin_cld = 1.e-10,&! Minimum boundary radius (m) |
|---|
| 189 | rbmax_cld = 2.e-4, &! Maximum boundary radius (m) |
|---|
| 190 | Fo = 7.5e-7, &! for sat index (J.m-3) |
|---|
| 191 | lambdaH = 150.e3 ! for sat index (km) |
|---|
| 192 | |
|---|
| 193 | character(len=23), parameter :: & |
|---|
| 194 | file_qext = 'optprop_co2ice_1mic.dat' ! File extinction coefficients of CO2 particles |
|---|
| 195 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 196 | !-----2) Saved: |
|---|
| 197 | !-------------- |
|---|
| 198 | integer, save :: & |
|---|
| 199 | imicroco2 ! Time subsampling for coupled water microphysics sedimentation microtimestep timeloop for microphysics: |
|---|
| 200 | ! if imicroco2 = 1, subpdt is the same as pdt |
|---|
| 201 | real, save :: & |
|---|
| 202 | sigma_iceco2, &! Variance of the ice and CCN distributions |
|---|
| 203 | microtimestep ! Integration timestep for coupled water microphysics & sedimentation |
|---|
| 204 | |
|---|
| 205 | double precision, save :: & |
|---|
| 206 | dev2, &! 1. / ( sqrt(2.) * sigma_iceco2 ) |
|---|
| 207 | Qext1bins(nbinco2_cld), &! Extinction coefficients for rb_cldco2 radius of CO2 ice particles |
|---|
| 208 | Qextv1mic(var_dim_qext), & |
|---|
| 209 | radv(var_dim_qext), & ! radius of CO2 ice at 1 µm (read from file_qext) |
|---|
| 210 | rb_cldco2(nbinco2_cld+1) ! boundary values of each rad_cldco2 bin (m) |
|---|
| 211 | logical, save :: & |
|---|
| 212 | firstcall = .true. ! Used to compute saved variables |
|---|
| 213 | |
|---|
| 214 | !$OMP THREADPRIVATE(imicroco2,sigma_iceco2,microtimestep) |
|---|
| 215 | !$OMP THREADPRIVATE(dev2,Qext1bins,Qextv1mic,radv,rb_cldco2) |
|---|
| 216 | !$OMP THREADPRIVATE(firstcall) |
|---|
| 217 | |
|---|
| 218 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 219 | !-----3) Variables: |
|---|
| 220 | !------------------ |
|---|
| 221 | integer :: & |
|---|
| 222 | iq, &! loop on tracers |
|---|
| 223 | ig, &! loop on grid points |
|---|
| 224 | l, &! loop on layers |
|---|
| 225 | i, &! loop on nbinco2_cld |
|---|
| 226 | nelem, &! number of point between lebon1 and lebon2 => interpolation |
|---|
| 227 | lebon1, &! bound limit for the interpolation |
|---|
| 228 | lebon2, &! bound limit for the interpolation |
|---|
| 229 | microstep ! Time subsampling step variable |
|---|
| 230 | |
|---|
| 231 | real :: & |
|---|
| 232 | ! ---Tendency given by clouds inside the micro loop |
|---|
| 233 | subpdqcloudco2(ngrid,nlay,nq), &! On tracers, cf. pdqcloud |
|---|
| 234 | subpdtcloudco2(ngrid,nlay), &! On temperature, cf. pdtcloud |
|---|
| 235 | ! ---Global tendency (clouds+physics) |
|---|
| 236 | sum_subpdq(ngrid,nlay,nq), &! On tracers, cf. pdqcloud |
|---|
| 237 | sum_subpdt(ngrid,nlay), &! On temperature, cf. pdtcloud |
|---|
| 238 | ! ---Sedimentation |
|---|
| 239 | ztsed(ngrid,nlay), &! Temperature with real-time value in microtimeloop |
|---|
| 240 | zqsed(ngrid,nlay,nq), &! Tracers with real-time value in µloop |
|---|
| 241 | zqsed0(ngrid,nlay,nq), &! For sedimentation tendancy |
|---|
| 242 | subpdqsed(ngrid,nlay,nq), &! Tendancy due to sedimentation |
|---|
| 243 | sum_subpdqs_sedco2(ngrid), &! CO2 flux at the surface |
|---|
| 244 | sum_subpdqs_sedccn(ngrid,nq), &! |
|---|
| 245 | ! ---For sub grid T distribution |
|---|
| 246 | zt(ngrid,nlay), &! Local value of temperature |
|---|
| 247 | zq_co2vap(ngrid, nlay), &! Local value of CO2 vap |
|---|
| 248 | rhocloudco2t(ngrid, nlay), &! Cloud density (kg.m-3) |
|---|
| 249 | ! ---For Saturation Index computation |
|---|
| 250 | zdelt, &! Delta T for the temperature distribution |
|---|
| 251 | co2cloudfrac(ngrid,nlay), &! Cloud fraction used only with CLFvarying is true |
|---|
| 252 | ! ---Misc |
|---|
| 253 | rhocloudco2(ngrid, nlay), &! Cloud density (kg.m-3) |
|---|
| 254 | Nccnco2, &! buffer: number of ccn used for co2 condensation |
|---|
| 255 | Qccnco2, &! buffer: mass of ccn used for co2 condensation |
|---|
| 256 | Nccnco2_h2o, & |
|---|
| 257 | Qccnco2_h2o, & |
|---|
| 258 | Niceco2, &! buffer: mmr co2 ice |
|---|
| 259 | epaisseur(ngrid,nlay), &! Layer thickness (m) |
|---|
| 260 | masse(ngrid,nlay), &! Layer mass (kg.m-2) |
|---|
| 261 | pteff(ngrid, nlay), &! Effective temperature in the cloud |
|---|
| 262 | pqeff(ngrid, nlay, nq), &! Effective tracers quantities in the cloud |
|---|
| 263 | wq(ngrid,nlay+1), &! Displaced tracer mass (kg.m-2) during microtimestep |
|---|
| 264 | satuco2(ngrid,nlay), &! CO2 satu ratio for output diagfi |
|---|
| 265 | zqsatco2(ngrid,nlay), &! Saturation co2 |
|---|
| 266 | myT, &! Temperature scalar for co2 density computation |
|---|
| 267 | availco2,& |
|---|
| 268 | masslayer, & |
|---|
| 269 | tmp, a,b, & |
|---|
| 270 | new_pdq(ngrid,nlay) |
|---|
| 271 | |
|---|
| 272 | double precision :: & |
|---|
| 273 | ! ---Extinction coefficients at 1 micron of CO2 particles |
|---|
| 274 | vrat_cld, &! Volume ratio |
|---|
| 275 | n_derf, &! derf( (rb_cldco2(1)-log(riceco2(ig,l))) *dev2) |
|---|
| 276 | Qtemp, &! mean value in the interval during the interpolation |
|---|
| 277 | ltemp1(var_dim_qext), &! abs(radv(:)-rb_cldco2(i)) |
|---|
| 278 | ltemp2(var_dim_qext), &! abs(radv(:)-rb_cldco2(i+1)) |
|---|
| 279 | n_aer(nbinco2_cld), &! -0.5 * Nccnco2*tauscaling(ig) * n_derf |
|---|
| 280 | tau1mic(ngrid), &! CO2 ice column opacity at 1µm |
|---|
| 281 | Qext1bins2(ngrid,nlay), &! CO2 ice opacities |
|---|
| 282 | ! ---For Saturation Index computation |
|---|
| 283 | rho, &! background density |
|---|
| 284 | zu, &! absolute value of zonal wind field |
|---|
| 285 | NN, &! N^2 static stability |
|---|
| 286 | gradT, &! thermal gradient |
|---|
| 287 | SatIndex(ngrid,nlay), &! Saturation index S in Spiga 2012 paper, assuming like in the paper GW phase speed |
|---|
| 288 | ! (stationary waves): c = 0 m.s-1, lambdaH = 150 km, Fo = 7.5e-7 J.m-3 |
|---|
| 289 | SatIndexmap(ngrid), &! maxval(SatIndex(ig,12:26)) |
|---|
| 290 | ! ---Misc |
|---|
| 291 | tcond(ngrid,nlay) ! CO2 condensation temperature |
|---|
| 292 | |
|---|
| 293 | logical :: & |
|---|
| 294 | file_qext_ok ! Check if file_qext exists |
|---|
| 295 | !======================================================================================================================! |
|---|
| 296 | ! BEGIN ===============================================================================================================! |
|---|
| 297 | !======================================================================================================================! |
|---|
| 298 | ! 0. Firstcall |
|---|
| 299 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 300 | if (firstcall) then |
|---|
| 301 | firstcall=.false. |
|---|
| 302 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 303 | ! 0.1. Initialization of microtimestep from imicroco2 |
|---|
| 304 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 305 | #ifdef MESOSCALE |
|---|
| 306 | imicroco2 = 2 |
|---|
| 307 | #else |
|---|
| 308 | imicroco2 = 30 |
|---|
| 309 | #endif |
|---|
| 310 | call getin("imicroco2", imicroco2) |
|---|
| 311 | microtimestep = ptimestep/real(imicroco2) |
|---|
| 312 | sigma_iceco2 = sqrt(log(1.+nuiceco2_sed)) |
|---|
| 313 | dev2 = 1. / ( sqrt(2.) * sigma_iceco2 ) |
|---|
| 314 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 315 | ! 0.2. Compute the radius grid of CO2 ice particles (rb_cldco2) |
|---|
| 316 | ! > the grid spacing is computed assuming a constant volume ratio between two consecutive bins; i.e. vrat_cld. |
|---|
| 317 | ! - rad_cldco2 is the primary radius grid used for microphysics computation. |
|---|
| 318 | ! - The grid spacing is computed assuming a constant volume ratio between two consecutive bins; i.e. vrat_cld. |
|---|
| 319 | ! - vrat_cld is determined from the boundary values of the size grid: rmin_cld and rmax_cld. |
|---|
| 320 | ! - The rb_cldco2 array contains the boundary values of each rad_cldco2 bin. |
|---|
| 321 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 322 | ! vrat_cld is determined from the boundary values of the size grid: rmin_cld and rmax_cld. |
|---|
| 323 | vrat_cld = exp(log(rmax_cld/rmin_cld) / float(nbinco2_cld-1) * 3.) |
|---|
| 324 | |
|---|
| 325 | ! rad_cldco2 is the primary radius grid used for microphysics computation. |
|---|
| 326 | rad_cldco2(1) = rmin_cld |
|---|
| 327 | do i = 1, nbinco2_cld-1 |
|---|
| 328 | rad_cldco2(i+1) = rad_cldco2(i) * vrat_cld**(1./3.) |
|---|
| 329 | end do |
|---|
| 330 | |
|---|
| 331 | ! rb_cldco2 array contains the boundary values of each rad_cldco2 bin. |
|---|
| 332 | rb_cldco2(1) = rbmin_cld |
|---|
| 333 | do i = 1, nbinco2_cld |
|---|
| 334 | rb_cldco2(i+1) = ( (2.*vrat_cld) / (vrat_cld+1.) )**(1./3.) * rad_cldco2(i) |
|---|
| 335 | end do |
|---|
| 336 | rb_cldco2(nbinco2_cld+1) = rbmax_cld |
|---|
| 337 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 338 | ! 0.3. Read file 'optprop_co2ice_1mic.dat' to extract optical properties of CO2 ice at 1 micron (Qext) |
|---|
| 339 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 340 | ! get information about file_qext |
|---|
| 341 | inquire(file=trim(datadir)//'/'//file_qext, exist=file_qext_ok) |
|---|
| 342 | |
|---|
| 343 | ! if file_qext is missing then stop |
|---|
| 344 | if (.not. file_qext_ok) then |
|---|
| 345 | write(*,*)'file'//file_qext//'should be in ', trim(datadir) |
|---|
| 346 | call abort_physic('co2cloud', 'file missing', 1) |
|---|
| 347 | end if |
|---|
| 348 | |
|---|
| 349 | ! read file_qext |
|---|
| 350 | open(unit=uQext,file=trim(datadir)//'/'//file_qext, form='formatted') |
|---|
| 351 | |
|---|
| 352 | ! skip 1 line |
|---|
| 353 | read(uQext,*) |
|---|
| 354 | |
|---|
| 355 | ! extract radv |
|---|
| 356 | do i = 1, var_dim_qext |
|---|
| 357 | read(uQext,'(E12.5)')radv(i) |
|---|
| 358 | end do |
|---|
| 359 | |
|---|
| 360 | ! skip 1 line |
|---|
| 361 | read(uQext,*) |
|---|
| 362 | |
|---|
| 363 | ! Qextv1mic |
|---|
| 364 | do i = 1 , var_dim_qext |
|---|
| 365 | read(uQext,'(E12.5)')Qextv1mic(i) |
|---|
| 366 | end do |
|---|
| 367 | |
|---|
| 368 | ! close file_qext |
|---|
| 369 | close(uQext) |
|---|
| 370 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 371 | ! 0.4. Interpole the radius grid (rb_cldco2) to get the corresponding exctinction coefficients (Qext1bins), using |
|---|
| 372 | ! file_qext values (radv, Qextv1mic) |
|---|
| 373 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 374 | do i = 1, nbinco2_cld |
|---|
| 375 | ltemp1 = abs(radv(:)-rb_cldco2(i)) |
|---|
| 376 | ltemp2 = abs(radv(:)-rb_cldco2(i+1)) |
|---|
| 377 | lebon1 = minloc(ltemp1,DIM=1) |
|---|
| 378 | lebon2 = min(minloc(ltemp2,DIM=1), var_dim_qext) |
|---|
| 379 | nelem = lebon2 - lebon1 + 1. |
|---|
| 380 | |
|---|
| 381 | ! mean value in the interval |
|---|
| 382 | Qtemp = 0d0 |
|---|
| 383 | do l = 0, nelem |
|---|
| 384 | Qtemp = Qtemp + Qextv1mic(min(lebon1+l, var_dim_qext)) |
|---|
| 385 | end do |
|---|
| 386 | |
|---|
| 387 | Qext1bins(i) = Qtemp / nelem |
|---|
| 388 | end do |
|---|
| 389 | |
|---|
| 390 | Qext1bins(:) = Qext1bins(:) * pi * (rad_cldco2(:)**2) |
|---|
| 391 | |
|---|
| 392 | ! print result of the interpolation |
|---|
| 393 | write(*,*)'--------------------------------------------' |
|---|
| 394 | write(*,*)'Microphysics co2: size bin-Qext information:' |
|---|
| 395 | write(*,*)' i, rad_cldco2(i), Qext1bins(i)' |
|---|
| 396 | do i = 1, nbinco2_cld |
|---|
| 397 | write(*,'(i3,3x,3(e13.6,4x))')i, rad_cldco2(i), Qext1bins(i) |
|---|
| 398 | end do |
|---|
| 399 | write(*,*)'--------------------------------------------' |
|---|
| 400 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 401 | ! 0.5. Save the radius grid of CO2 particles (rb_cldco2) |
|---|
| 402 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 403 | do i = 1, nbinco2_cld+1 |
|---|
| 404 | rb_cldco2(i) = log(rb_cldco2(i)) |
|---|
| 405 | end do |
|---|
| 406 | end if ! of IF (firstcall) |
|---|
| 407 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 408 | ! 1. Initialization |
|---|
| 409 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 410 | sum_subpdq(1:ngrid,1:nlay,1:nq) = 0. |
|---|
| 411 | sum_subpdt(1:ngrid,1:nlay) = 0. |
|---|
| 412 | |
|---|
| 413 | subpdqcloudco2(1:ngrid,1:nlay,1:nq) = 0. |
|---|
| 414 | subpdtcloudco2(1:ngrid,1:nlay) = 0. |
|---|
| 415 | |
|---|
| 416 | pdqcloudco2(1:ngrid,1:nlay,1:nq) = 0. |
|---|
| 417 | pdtcloudco2(1:ngrid,1:nlay) = 0. |
|---|
| 418 | |
|---|
| 419 | ! default value if no ice |
|---|
| 420 | rhocloudco2(1:ngrid,1:nlay) = rho_dust |
|---|
| 421 | rhocloudco2t(1:ngrid,1:nlay) = rho_dust |
|---|
| 422 | epaisseur(1:ngrid,1:nlay) = 0. |
|---|
| 423 | masse(1:ngrid,1:nlay) = 0. |
|---|
| 424 | riceco2(1:ngrid, 1:nlay) = 0. |
|---|
| 425 | zqsed0(1:ngrid,1:nlay,1:nq) = 0. |
|---|
| 426 | sum_subpdqs_sedco2(1:ngrid) = 0. |
|---|
| 427 | sum_subpdqs_sedccn(1:ngrid,1:nq) = 0. |
|---|
| 428 | subpdqsed(1:ngrid,1:nlay,1:nq) = 0. |
|---|
| 429 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 430 | ! 2. Compute mass and thickness layers |
|---|
| 431 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 432 | do l = 1, nlay |
|---|
| 433 | do ig = 1, ngrid |
|---|
| 434 | #ifndef MESOSCALE |
|---|
| 435 | masse(ig,l) = (pplev(ig,l) - pplev(ig,l+1) + (bp(l)-bp(l+1)) ) / g |
|---|
| 436 | #else |
|---|
| 437 | masse(ig,l) = (pplev(ig,l) - pplev(ig,l+1)) / g |
|---|
| 438 | #endif |
|---|
| 439 | epaisseur(ig,l) = pzlev(ig,l+1) - pzlev(ig,l) |
|---|
| 440 | end do |
|---|
| 441 | end do |
|---|
| 442 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 443 | ! 3. Define the sub-grid cloud (CLFvaryingCO2) |
|---|
| 444 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 445 | ! 3.1. Representation of sub-grid CO2 ice clouds (CLFvaryingCO2 = True) |
|---|
| 446 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 447 | if (CLFvaryingCO2) then |
|---|
| 448 | ! effective temperature |
|---|
| 449 | pteff(:,:) = pt(:,:) |
|---|
| 450 | |
|---|
| 451 | ! min co2cloudfrac when there is ice |
|---|
| 452 | co2cloudfrac(:,:) = mincloud |
|---|
| 453 | |
|---|
| 454 | ! temperature |
|---|
| 455 | do l=1,nlay |
|---|
| 456 | do ig=1,ngrid |
|---|
| 457 | zt(ig,l) = pt(ig,l) + pdt(ig,l)*ptimestep |
|---|
| 458 | end do |
|---|
| 459 | end do |
|---|
| 460 | |
|---|
| 461 | ! Quantities of traceurs |
|---|
| 462 | if (igcm_co2 /= 0) then |
|---|
| 463 | do l = 1, nlay |
|---|
| 464 | do ig = 1, ngrid |
|---|
| 465 | zq_co2vap(ig,l) = pq(ig,l,igcm_co2) + pdq(ig,l,igcm_co2)*ptimestep |
|---|
| 466 | end do |
|---|
| 467 | end do |
|---|
| 468 | end if |
|---|
| 469 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 470 | ! 3.1.a. Saturation index CO2 |
|---|
| 471 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 472 | ! if saturation index co2 is true |
|---|
| 473 | if (satindexco2) then |
|---|
| 474 | ! layers 12 --> 26 ~ 12->85 km |
|---|
| 475 | do l = 12, 26 |
|---|
| 476 | do ig = 1, ngrid |
|---|
| 477 | ! compute N^2 static stability |
|---|
| 478 | gradT = (zt(ig,l+1)-zt(ig,l))/(pzlev(ig,l+1)-pzlev(ig,l)) |
|---|
| 479 | NN = sqrt(g/zt(iq,l) * (max(gradT,-g/cpp) + g/cpp)) |
|---|
| 480 | |
|---|
| 481 | ! compute absolute value of zonal wind field |
|---|
| 482 | zu = abs(pu(ig,l) + pdu(ig,l)*ptimestep) |
|---|
| 483 | |
|---|
| 484 | ! compute background density |
|---|
| 485 | rho = pplay(ig,l) / (rnew(ig,l)*zt(ig,l)) |
|---|
| 486 | |
|---|
| 487 | ! saturation index: Modulate the DeltaT by GW propagation index: |
|---|
| 488 | ! -------------------------------------------------------------- |
|---|
| 489 | SatIndex(ig,l) = sqrt(Fo*lambdaH/(2.*pi)*NN / (rho*zu**3) ) |
|---|
| 490 | end do |
|---|
| 491 | end do |
|---|
| 492 | |
|---|
| 493 | ! Then compute Satindex map in layers 12 --> 26 ~ 12->85 km |
|---|
| 494 | do ig = 1, ngrid |
|---|
| 495 | SatIndexmap(ig) = maxval(SatIndex(ig,12:26)) |
|---|
| 496 | end do |
|---|
| 497 | |
|---|
| 498 | #ifndef MESOSCALE |
|---|
| 499 | ! Write outputs in diagfi.nc |
|---|
| 500 | call write_output( "SatIndex", "SatIndex", " ", SatIndex(:,:)) |
|---|
| 501 | |
|---|
| 502 | call write_output( "SatIndexmap", "SatIndexmap", "km", SatIndexmap(:)) |
|---|
| 503 | #endif |
|---|
| 504 | !------------------------------------------------------------------------------------------------------------------! |
|---|
| 505 | ! if saturation index co2 is false, set saturation index to 0.05 |
|---|
| 506 | !------------------------------------------------------------------------------------------------------------------! |
|---|
| 507 | else |
|---|
| 508 | do ig = 1, ngrid |
|---|
| 509 | SatIndexmap(ig)=0.05 |
|---|
| 510 | end do |
|---|
| 511 | end if ! of if (satindexco2) |
|---|
| 512 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 513 | ! 3.1.b. Compute tcond |
|---|
| 514 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 515 | call tcondco2(ngrid,nlay,pplay,zq_co2vap,tcond) |
|---|
| 516 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 517 | ! 3.1.c. Compute cloud fraction in cells |
|---|
| 518 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 519 | do ig = 1, ngrid |
|---|
| 520 | if (SatIndexmap(ig) <= 0.1) then |
|---|
| 521 | do l = 1, nlay-1 |
|---|
| 522 | |
|---|
| 523 | ! The entire fraction is saturated |
|---|
| 524 | if (tcond(ig,l) >= (zt(ig,l)+zdelt) .or. tcond(ig,l) <= 0.) then |
|---|
| 525 | pteff(ig,l) = zt(ig,l) |
|---|
| 526 | co2cloudfrac(ig,l) = 1. |
|---|
| 527 | |
|---|
| 528 | ! No saturation at all |
|---|
| 529 | else if (tcond(ig,l) <= (zt(ig,l)-zdelt)) then |
|---|
| 530 | pteff(ig,l) = zt(ig,l) - zdelt |
|---|
| 531 | co2cloudfrac(ig,l) = mincloud |
|---|
| 532 | |
|---|
| 533 | ! Mean temperature of the cloud fraction |
|---|
| 534 | else |
|---|
| 535 | pteff(ig,l) = (tcond(ig,l)+zt(ig,l)-zdelt) / 2. |
|---|
| 536 | co2cloudfrac(ig,l) = (tcond(ig,l)-zt(ig,l)+zdelt) / (2.0*zdelt) |
|---|
| 537 | end if |
|---|
| 538 | |
|---|
| 539 | pteff(ig,l) = pteff(ig,l) - pdt(ig,l)*ptimestep |
|---|
| 540 | |
|---|
| 541 | ! check boundary values of co2cloudfrac |
|---|
| 542 | if (co2cloudfrac(ig,l) <= mincloud) then |
|---|
| 543 | co2cloudfrac(ig,l) = mincloud |
|---|
| 544 | else if (co2cloudfrac(ig,l)> 1) then |
|---|
| 545 | co2cloudfrac(ig,l) = 1. |
|---|
| 546 | end if |
|---|
| 547 | end do |
|---|
| 548 | |
|---|
| 549 | ! SatIndex not favorable for GW: leave pt untouched |
|---|
| 550 | else |
|---|
| 551 | pteff(ig,l) = pt(ig,l) |
|---|
| 552 | co2cloudfrac(ig,l) = mincloud |
|---|
| 553 | end if ! of if (SatIndexmap <= 0.1) |
|---|
| 554 | end do ! of ngrid |
|---|
| 555 | ! TODO: Totalcloud frac of the column missing here |
|---|
| 556 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 557 | ! 3.2. No sub-grid cloud representation (CLFvaryingCO2 = False) |
|---|
| 558 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 559 | else |
|---|
| 560 | do l = 1, nlay |
|---|
| 561 | do ig = 1, ngrid |
|---|
| 562 | pteff(ig,l) = pt(ig,l) |
|---|
| 563 | end do |
|---|
| 564 | end do |
|---|
| 565 | end if ! end if (CLFvaryingco2) |
|---|
| 566 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 567 | ! 4. Microphysics of CO2 cloud formation |
|---|
| 568 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 569 | pqeff(:,:,:) = pq(:,:,:) |
|---|
| 570 | pteff(:,:) = pt(:,:) |
|---|
| 571 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 572 | ! 4.1. Effective tracers quantities in the cloud |
|---|
| 573 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 574 | if (CLFvaryingCO2) then |
|---|
| 575 | pqeff(:,:,igcm_ccnco2_mass) = pq(:,:,igcm_ccnco2_mass) / co2cloudfrac(:,:) |
|---|
| 576 | |
|---|
| 577 | pqeff(:,:,igcm_ccnco2_number) = pq(:,:,igcm_ccnco2_number) / co2cloudfrac(:,:) |
|---|
| 578 | |
|---|
| 579 | pqeff(:,:,igcm_co2_ice) = pq(:,:,igcm_co2_ice) / co2cloudfrac(:,:) |
|---|
| 580 | end if |
|---|
| 581 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 582 | do microstep = 1, imicroco2 |
|---|
| 583 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 584 | ! 4.2. Stepped entry for tendancies: At each micro timestep we add pdt in order to have a stepped entry |
|---|
| 585 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 586 | do l = 1, nlay |
|---|
| 587 | do ig = 1, ngrid |
|---|
| 588 | ! on temperature |
|---|
| 589 | sum_subpdt(ig,l) = sum_subpdt(ig,l) + pdt(ig,l) |
|---|
| 590 | |
|---|
| 591 | ! on tracers |
|---|
| 592 | sum_subpdq(ig,l,igcm_dust_mass) = sum_subpdq(ig,l,igcm_dust_mass) + pdq(ig,l,igcm_dust_mass) |
|---|
| 593 | |
|---|
| 594 | sum_subpdq(ig,l,igcm_dust_number) = sum_subpdq(ig,l,igcm_dust_number) + pdq(ig,l,igcm_dust_number) |
|---|
| 595 | |
|---|
| 596 | sum_subpdq(ig,l,igcm_ccnco2_mass) = sum_subpdq(ig,l,igcm_ccnco2_mass) + pdq(ig,l,igcm_ccnco2_mass) |
|---|
| 597 | |
|---|
| 598 | sum_subpdq(ig,l,igcm_ccnco2_number) = sum_subpdq(ig,l,igcm_ccnco2_number) + pdq(ig,l,igcm_ccnco2_number) |
|---|
| 599 | |
|---|
| 600 | sum_subpdq(ig,l,igcm_co2_ice) = sum_subpdq(ig,l,igcm_co2_ice) + pdq(ig,l,igcm_co2_ice) |
|---|
| 601 | |
|---|
| 602 | sum_subpdq(ig,l,igcm_co2) = sum_subpdq(ig,l,igcm_co2) + pdq(ig,l,igcm_co2) |
|---|
| 603 | |
|---|
| 604 | if (meteo_flux) then |
|---|
| 605 | sum_subpdq(ig,l,igcm_ccnco2_meteor_number) = sum_subpdq(ig,l,igcm_ccnco2_meteor_number) + & |
|---|
| 606 | pdq(ig,l,igcm_ccnco2_meteor_number) |
|---|
| 607 | |
|---|
| 608 | sum_subpdq(ig,l,igcm_ccnco2_meteor_mass) = sum_subpdq(ig,l,igcm_ccnco2_meteor_mass) + & |
|---|
| 609 | pdq(ig,l,igcm_ccnco2_meteor_mass) |
|---|
| 610 | end if |
|---|
| 611 | if (co2useh2o) then |
|---|
| 612 | sum_subpdq(ig,l,igcm_h2o_ice) = sum_subpdq(ig,l,igcm_h2o_ice) + pdq(ig,l,igcm_h2o_ice) |
|---|
| 613 | |
|---|
| 614 | sum_subpdq(ig,l,igcm_ccn_mass) = sum_subpdq(ig,l,igcm_ccn_mass) + pdq(ig,l,igcm_ccn_mass) |
|---|
| 615 | |
|---|
| 616 | sum_subpdq(ig,l,igcm_ccn_number) = sum_subpdq(ig,l,igcm_ccn_number) + pdq(ig,l,igcm_ccn_number) |
|---|
| 617 | |
|---|
| 618 | sum_subpdq(ig,l,igcm_ccnco2_h2o_number) = sum_subpdq(ig,l,igcm_ccnco2_h2o_number) + & |
|---|
| 619 | pdq(ig,l,igcm_ccnco2_h2o_number) |
|---|
| 620 | |
|---|
| 621 | sum_subpdq(ig,l,igcm_ccnco2_h2o_mass_ice) = sum_subpdq(ig,l,igcm_ccnco2_h2o_mass_ice) + & |
|---|
| 622 | pdq(ig,l,igcm_ccnco2_h2o_mass_ice) |
|---|
| 623 | |
|---|
| 624 | sum_subpdq(ig,l,igcm_ccnco2_h2o_mass_ccn) = sum_subpdq(ig,l,igcm_ccnco2_h2o_mass_ccn) + & |
|---|
| 625 | pdq(ig,l,igcm_ccnco2_h2o_mass_ccn) |
|---|
| 626 | end if |
|---|
| 627 | end do ! ngrid |
|---|
| 628 | end do ! nlay |
|---|
| 629 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 630 | ! 4.3. Main call to the cloud scheme |
|---|
| 631 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 632 | call improvedco2clouds(ngrid, nlay, microtimestep, pplay, pplev, pteff, sum_subpdt, pqeff, sum_subpdq, & |
|---|
| 633 | subpdqcloudco2, subpdtcloudco2, nq, tauscaling, rb_cldco2, sigma_iceco2, dev2) |
|---|
| 634 | |
|---|
| 635 | do l = 1, nlay |
|---|
| 636 | do ig = 1, ngrid |
|---|
| 637 | if(pq(ig,l,igcm_co2_ice) + microtimestep*(sum_subpdq(ig,l,igcm_co2_ice)+subpdqcloudco2(ig,l,igcm_co2_ice)) & |
|---|
| 638 | <= 1.e-12) then |
|---|
| 639 | subpdqcloudco2(ig,l,igcm_co2_ice) = -pq(ig,l,igcm_co2_ice)/microtimestep - sum_subpdq(ig,l,igcm_co2_ice) |
|---|
| 640 | subpdqcloudco2(ig,l,igcm_co2) = -subpdqcloudco2(ig,l,igcm_co2_ice) |
|---|
| 641 | end if |
|---|
| 642 | |
|---|
| 643 | if(pq(ig,l,igcm_co2) + microtimestep*(sum_subpdq(ig,l,igcm_co2)+subpdqcloudco2(ig,l,igcm_co2)) <= 1.e-12) then |
|---|
| 644 | subpdqcloudco2(ig,l,igcm_co2) = - pq(ig,l,igcm_co2)/microtimestep - sum_subpdq(ig,l,igcm_co2) |
|---|
| 645 | subpdqcloudco2(ig,l,igcm_co2_ice) = -subpdqcloudco2(ig,l,igcm_co2) |
|---|
| 646 | end if |
|---|
| 647 | |
|---|
| 648 | ! ccnco2_number and ccnco2_mass |
|---|
| 649 | if (((pq(ig,l,igcm_ccnco2_number)+(sum_subpdq(ig,l,igcm_ccnco2_number)+subpdqcloudco2(ig,l,igcm_ccnco2_number)) & |
|---|
| 650 | *microtimestep)<=1.) .or. & |
|---|
| 651 | (pq(ig,l,igcm_ccnco2_mass)+(sum_subpdq(ig,l,igcm_ccnco2_mass)+subpdqcloudco2(ig,l,igcm_ccnco2_mass)) & |
|---|
| 652 | *microtimestep<=1e-20)) then |
|---|
| 653 | subpdqcloudco2(ig,l,igcm_ccnco2_number) = - pq(ig,l,igcm_ccnco2_number)/microtimestep + 1. & |
|---|
| 654 | - sum_subpdq(ig,l,igcm_ccnco2_number) |
|---|
| 655 | subpdqcloudco2(ig,l,igcm_dust_number) = - subpdqcloudco2(ig,l,igcm_ccnco2_number) |
|---|
| 656 | |
|---|
| 657 | subpdqcloudco2(ig,l,igcm_ccnco2_mass) = - pq(ig,l,igcm_ccnco2_mass)/microtimestep + 1e-20 & |
|---|
| 658 | - sum_subpdq(ig,l,igcm_ccnco2_mass) |
|---|
| 659 | subpdqcloudco2(ig,l,igcm_dust_mass) = - subpdqcloudco2(ig,l,igcm_ccnco2_mass) |
|---|
| 660 | end if |
|---|
| 661 | |
|---|
| 662 | ! ccnco2_meteor_number and ccnco2_meteor_mass |
|---|
| 663 | if (meteo_flux) then |
|---|
| 664 | if (((pq(ig,l,igcm_ccnco2_meteor_number)+(sum_subpdq(ig,l,igcm_ccnco2_meteor_number)+ & |
|---|
| 665 | subpdqcloudco2(ig,l,igcm_ccnco2_meteor_number))*microtimestep)<=1.) .or. & |
|---|
| 666 | (pq(ig,l,igcm_ccnco2_meteor_mass)+(sum_subpdq(ig,l,igcm_ccnco2_meteor_mass)+ & |
|---|
| 667 | subpdqcloudco2(ig,l,igcm_ccnco2_meteor_mass))*microtimestep<=1e-20)) then |
|---|
| 668 | subpdqcloudco2(ig,l,igcm_ccnco2_meteor_number) = - pq(ig,l,igcm_ccnco2_meteor_number)/microtimestep + 1. & |
|---|
| 669 | - sum_subpdq(ig,l,igcm_ccnco2_meteor_number) |
|---|
| 670 | subpdqcloudco2(ig,l,igcm_ccnco2_meteor_mass) = - pq(ig,l,igcm_ccnco2_meteor_mass)/microtimestep + 1e-20 & |
|---|
| 671 | - sum_subpdq(ig,l,igcm_ccnco2_meteor_mass) |
|---|
| 672 | end if |
|---|
| 673 | end if |
|---|
| 674 | |
|---|
| 675 | ! ccnco2_h2o_number and masses |
|---|
| 676 | if (co2useh2o) then |
|---|
| 677 | if (((pq(ig,l,igcm_ccnco2_h2o_number) + (sum_subpdq(ig,l,igcm_ccnco2_h2o_number) + & |
|---|
| 678 | subpdqcloudco2(ig,l,igcm_ccnco2_h2o_number)) & |
|---|
| 679 | *microtimestep)<=1.) .or. & |
|---|
| 680 | (pq(ig,l,igcm_ccnco2_h2o_mass_ice)+pq(ig,l,igcm_ccnco2_h2o_mass_ccn) +& |
|---|
| 681 | (sum_subpdq(ig,l,igcm_ccnco2_h2o_mass_ice)+subpdqcloudco2(ig,l,igcm_ccnco2_h2o_mass_ice) +& |
|---|
| 682 | sum_subpdq(ig,l,igcm_ccnco2_h2o_mass_ccn)+subpdqcloudco2(ig,l,igcm_ccnco2_h2o_mass_ccn)) & |
|---|
| 683 | *microtimestep<=1e-20)) then |
|---|
| 684 | |
|---|
| 685 | subpdqcloudco2(ig,l,igcm_ccnco2_h2o_number) = - pq(ig,l,igcm_ccnco2_h2o_number)/microtimestep + 1. & |
|---|
| 686 | - sum_subpdq(ig,l,igcm_ccnco2_h2o_number) |
|---|
| 687 | subpdqcloudco2(ig,l,igcm_ccn_number) = - subpdqcloudco2(ig,l,igcm_ccnco2_h2o_number) |
|---|
| 688 | |
|---|
| 689 | subpdqcloudco2(ig,l,igcm_ccnco2_h2o_mass_ice) = - pq(ig,l,igcm_ccnco2_h2o_mass_ice)/microtimestep + 1e-20 & |
|---|
| 690 | - sum_subpdq(ig,l,igcm_ccnco2_h2o_mass_ice) |
|---|
| 691 | |
|---|
| 692 | subpdqcloudco2(ig,l,igcm_ccnco2_h2o_mass_ccn) = - pq(ig,l,igcm_ccnco2_h2o_mass_ccn)/microtimestep + 1e-20 & |
|---|
| 693 | - sum_subpdq(ig,l,igcm_ccnco2_h2o_mass_ccn) |
|---|
| 694 | |
|---|
| 695 | subpdqcloudco2(ig,l,igcm_ccn_mass) = - subpdqcloudco2(ig,l,igcm_ccnco2_h2o_mass_ice) & |
|---|
| 696 | - subpdqcloudco2(ig,l,igcm_ccnco2_h2o_mass_ccn) |
|---|
| 697 | end if |
|---|
| 698 | end if |
|---|
| 699 | end do |
|---|
| 700 | end do |
|---|
| 701 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 702 | ! 4.4. Updating tendencies after cloud scheme |
|---|
| 703 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 704 | do l = 1, nlay |
|---|
| 705 | do ig = 1, ngrid |
|---|
| 706 | sum_subpdt(ig,l) = sum_subpdt(ig,l) + subpdtcloudco2(ig,l) |
|---|
| 707 | |
|---|
| 708 | sum_subpdq(ig,l,igcm_dust_mass) = sum_subpdq(ig,l,igcm_dust_mass) + subpdqcloudco2(ig,l,igcm_dust_mass) |
|---|
| 709 | |
|---|
| 710 | sum_subpdq(ig,l,igcm_dust_number) = sum_subpdq(ig,l,igcm_dust_number) + subpdqcloudco2(ig,l,igcm_dust_number) |
|---|
| 711 | |
|---|
| 712 | sum_subpdq(ig,l,igcm_ccnco2_mass) = sum_subpdq(ig,l,igcm_ccnco2_mass) + subpdqcloudco2(ig,l,igcm_ccnco2_mass) |
|---|
| 713 | |
|---|
| 714 | sum_subpdq(ig,l,igcm_ccnco2_number) = sum_subpdq(ig,l,igcm_ccnco2_number) + & |
|---|
| 715 | subpdqcloudco2(ig,l,igcm_ccnco2_number) |
|---|
| 716 | |
|---|
| 717 | sum_subpdq(ig,l,igcm_co2_ice) = sum_subpdq(ig,l,igcm_co2_ice) + subpdqcloudco2(ig,l,igcm_co2_ice) |
|---|
| 718 | |
|---|
| 719 | sum_subpdq(ig,l,igcm_co2) = sum_subpdq(ig,l,igcm_co2) + subpdqcloudco2(ig,l,igcm_co2) |
|---|
| 720 | |
|---|
| 721 | if (meteo_flux) then |
|---|
| 722 | sum_subpdq(ig,l,igcm_ccnco2_meteor_mass) = sum_subpdq(ig,l,igcm_ccnco2_meteor_mass) + & |
|---|
| 723 | subpdqcloudco2(ig,l,igcm_ccnco2_meteor_mass) |
|---|
| 724 | |
|---|
| 725 | sum_subpdq(ig,l,igcm_ccnco2_meteor_number) = sum_subpdq(ig,l,igcm_ccnco2_meteor_number) + & |
|---|
| 726 | subpdqcloudco2(ig,l,igcm_ccnco2_meteor_number) |
|---|
| 727 | end if |
|---|
| 728 | if (co2useh2o) then |
|---|
| 729 | sum_subpdq(ig,l,igcm_h2o_ice) = sum_subpdq(ig,l,igcm_h2o_ice) + subpdqcloudco2(ig,l,igcm_h2o_ice) |
|---|
| 730 | |
|---|
| 731 | sum_subpdq(ig,l,igcm_ccn_mass) = sum_subpdq(ig,l,igcm_ccn_mass) + subpdqcloudco2(ig,l,igcm_ccn_mass) |
|---|
| 732 | |
|---|
| 733 | sum_subpdq(ig,l,igcm_ccn_number) = sum_subpdq(ig,l,igcm_ccn_number) + subpdqcloudco2(ig,l,igcm_ccn_number) |
|---|
| 734 | |
|---|
| 735 | sum_subpdq(ig,l,igcm_ccnco2_h2o_mass_ice) = sum_subpdq(ig,l,igcm_ccnco2_h2o_mass_ice) + & |
|---|
| 736 | subpdqcloudco2(ig,l,igcm_ccnco2_h2o_mass_ice) |
|---|
| 737 | |
|---|
| 738 | sum_subpdq(ig,l,igcm_ccnco2_h2o_mass_ccn) = sum_subpdq(ig,l,igcm_ccnco2_h2o_mass_ccn) + & |
|---|
| 739 | subpdqcloudco2(ig,l,igcm_ccnco2_h2o_mass_ccn) |
|---|
| 740 | |
|---|
| 741 | sum_subpdq(ig,l,igcm_ccnco2_h2o_number) = sum_subpdq(ig,l,igcm_ccnco2_h2o_number) + & |
|---|
| 742 | subpdqcloudco2(ig,l,igcm_ccnco2_number) |
|---|
| 743 | |
|---|
| 744 | end if |
|---|
| 745 | end do ! ngrid |
|---|
| 746 | end do ! nlay |
|---|
| 747 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 748 | ! 4.5. Gravitational sedimentation |
|---|
| 749 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 750 | if (sedimentation) then |
|---|
| 751 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 752 | ! 4.5.a. Compute cloud density |
|---|
| 753 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 754 | do l = 1, nlay |
|---|
| 755 | do ig = 1, ngrid |
|---|
| 756 | ! temperature during the sedimentation process |
|---|
| 757 | ztsed(ig,l) = pteff(ig,l) + sum_subpdt(ig,l) * microtimestep |
|---|
| 758 | |
|---|
| 759 | ! quantities tracers during the sedimentation process |
|---|
| 760 | zqsed(ig,l,:) = pqeff(ig,l,:) + sum_subpdq(ig,l,:) * microtimestep |
|---|
| 761 | |
|---|
| 762 | ! assure positive value of co2_ice mmr, ccnco2 number, ccnco2 mass |
|---|
| 763 | ! meteoritic particle are considered like dust, rho_dust |
|---|
| 764 | Niceco2 = max(zqsed(ig,l,igcm_co2_ice), threshold) |
|---|
| 765 | Nccnco2 = max(zqsed(ig,l,igcm_ccnco2_number), threshold) |
|---|
| 766 | Qccnco2 = max(zqsed(ig,l,igcm_ccnco2_mass), threshold) |
|---|
| 767 | |
|---|
| 768 | Nccnco2_h2o = 0. |
|---|
| 769 | Qccnco2_h2o = 0. |
|---|
| 770 | if (co2useh2o) then |
|---|
| 771 | Nccnco2_h2o = max(zqsed(ig,l,igcm_ccnco2_h2o_number), threshold) |
|---|
| 772 | Qccnco2_h2o = max(zqsed(ig,l,igcm_ccnco2_h2o_mass_ice)+zqsed(ig,l,igcm_ccnco2_h2o_mass_ccn), threshold) |
|---|
| 773 | Nccnco2 = Nccnco2 - Nccnco2_h2o |
|---|
| 774 | Qccnco2 = Qccnco2 - Qccnco2_h2o |
|---|
| 775 | if (Nccnco2 < 0.) then |
|---|
| 776 | Nccnco2 = threshold |
|---|
| 777 | Qccnco2 = threshold |
|---|
| 778 | end if |
|---|
| 779 | end if |
|---|
| 780 | |
|---|
| 781 | ! Get density cloud and co2 ice particle radius |
|---|
| 782 | if (Niceco2/=0d0) then |
|---|
| 783 | call updaterice_microco2(dble(Niceco2), dble(Qccnco2), dble(Nccnco2), dble(Qccnco2_h2o), dble(Nccnco2_h2o),& |
|---|
| 784 | ztsed(ig,l), tauscaling(ig), riceco2(ig,l), rhocloudco2t(ig,l)) |
|---|
| 785 | else |
|---|
| 786 | riceco2(ig,l) = 0. |
|---|
| 787 | rhocloudco2t(ig,l) = 0. |
|---|
| 788 | end if |
|---|
| 789 | end do ! ngrid |
|---|
| 790 | end do ! nlay |
|---|
| 791 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 792 | ! 4.5.b. Save actualized tracer values to compute sedimentation tendancies |
|---|
| 793 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 794 | zqsed0(:,:,igcm_co2_ice) = zqsed(:,:,igcm_co2_ice) |
|---|
| 795 | zqsed0(:,:,igcm_ccnco2_mass) = zqsed(:,:,igcm_ccnco2_mass) |
|---|
| 796 | zqsed0(:,:,igcm_ccnco2_number) = zqsed(:,:,igcm_ccnco2_number) |
|---|
| 797 | |
|---|
| 798 | if (meteo_flux) then |
|---|
| 799 | zqsed0(:,:,igcm_ccnco2_meteor_mass) = zqsed(:,:,igcm_ccnco2_meteor_mass) |
|---|
| 800 | zqsed0(:,:,igcm_ccnco2_meteor_number) = zqsed(:,:,igcm_ccnco2_meteor_number) |
|---|
| 801 | end if |
|---|
| 802 | |
|---|
| 803 | if (co2useh2o) then |
|---|
| 804 | zqsed0(:,:,igcm_ccnco2_h2o_number) = zqsed(:,:,igcm_ccnco2_h2o_number) |
|---|
| 805 | zqsed0(:,:,igcm_ccnco2_h2o_mass_ice) = zqsed(:,:,igcm_ccnco2_h2o_mass_ice) |
|---|
| 806 | zqsed0(:,:,igcm_ccnco2_h2o_mass_ccn) = zqsed(:,:,igcm_ccnco2_h2o_mass_ccn) |
|---|
| 807 | end if |
|---|
| 808 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 809 | ! 4.5.c. Sedimentation of co2 ice |
|---|
| 810 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 811 | do ig = 1, ngrid |
|---|
| 812 | do l = 1, nlay |
|---|
| 813 | rsedcloudco2(ig,l) = max( riceco2(ig,l)*(1.+nuiceco2_sed)*(1.+nuiceco2_sed)*(1.+nuiceco2_sed), & |
|---|
| 814 | rdust(ig,l) ) |
|---|
| 815 | end do |
|---|
| 816 | end do |
|---|
| 817 | |
|---|
| 818 | wq(:,:) = 0. |
|---|
| 819 | call newsedim(ngrid, nlay, ngrid*nlay, ngrid*nlay, microtimestep, pplev, masse, epaisseur, ztsed, & |
|---|
| 820 | rsedcloudco2, rhocloudco2t, zqsed(:,:,igcm_co2_ice), wq, beta) |
|---|
| 821 | |
|---|
| 822 | do ig = 1, ngrid |
|---|
| 823 | sum_subpdqs_sedco2(ig) = sum_subpdqs_sedco2(ig) + wq(ig,1) / microtimestep !wq est en kg.m-2 |
|---|
| 824 | end do |
|---|
| 825 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 826 | ! 4.5.d. Sedimentation for other tracers |
|---|
| 827 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 828 | wq(:,:) = 0. |
|---|
| 829 | ! for ccnco2_mass |
|---|
| 830 | call newsedim(ngrid, nlay, ngrid*nlay, ngrid*nlay, microtimestep, pplev, masse, epaisseur, ztsed, & |
|---|
| 831 | rsedcloudco2, rhocloudco2t, zqsed(:,:,igcm_ccnco2_mass), wq, beta) |
|---|
| 832 | do ig = 1, ngrid |
|---|
| 833 | sum_subpdqs_sedccn(ig,igcm_ccnco2_mass) = sum_subpdqs_sedccn(ig,igcm_ccnco2_mass) + wq(ig,1) / microtimestep |
|---|
| 834 | end do |
|---|
| 835 | |
|---|
| 836 | wq(:,:) = 0. |
|---|
| 837 | ! for ccnco2_number |
|---|
| 838 | call newsedim(ngrid, nlay, ngrid*nlay, ngrid*nlay,microtimestep, pplev, masse, epaisseur, ztsed, & |
|---|
| 839 | rsedcloudco2, rhocloudco2t, zqsed(:,:,igcm_ccnco2_number), wq, beta) |
|---|
| 840 | do ig = 1, ngrid |
|---|
| 841 | sum_subpdqs_sedccn(ig,igcm_ccnco2_number) = sum_subpdqs_sedccn(ig,igcm_ccnco2_number) + wq(ig,1) / microtimestep |
|---|
| 842 | end do |
|---|
| 843 | |
|---|
| 844 | if (meteo_flux) then |
|---|
| 845 | wq(:,:) = 0. |
|---|
| 846 | ! for ccnco2_meteor_mass |
|---|
| 847 | call newsedim(ngrid, nlay, ngrid*nlay, ngrid*nlay, microtimestep, pplev, masse, epaisseur, ztsed, & |
|---|
| 848 | rsedcloudco2, rhocloudco2t, zqsed(:,:,igcm_ccnco2_meteor_mass), wq, beta) |
|---|
| 849 | do ig = 1, ngrid |
|---|
| 850 | sum_subpdqs_sedccn(ig,igcm_ccnco2_meteor_mass) = sum_subpdqs_sedccn(ig,igcm_ccnco2_meteor_mass) + & |
|---|
| 851 | wq(ig,1)/ microtimestep |
|---|
| 852 | end do |
|---|
| 853 | |
|---|
| 854 | wq(:,:) = 0. |
|---|
| 855 | ! for ccnco2_meteor_number |
|---|
| 856 | call newsedim(ngrid, nlay, ngrid*nlay, ngrid*nlay,microtimestep, pplev, masse, epaisseur, ztsed, & |
|---|
| 857 | rsedcloudco2, rhocloudco2t, zqsed(:,:,igcm_ccnco2_meteor_number), wq, beta) |
|---|
| 858 | do ig = 1, ngrid |
|---|
| 859 | sum_subpdqs_sedccn(ig,igcm_ccnco2_meteor_number) = sum_subpdqs_sedccn(ig,igcm_ccnco2_meteor_number) + & |
|---|
| 860 | wq(ig,1) / microtimestep |
|---|
| 861 | end do |
|---|
| 862 | end if |
|---|
| 863 | ! for ccnco2_h2o_mass_ice |
|---|
| 864 | if (co2useh2o) then |
|---|
| 865 | wq(:,:) = 0. |
|---|
| 866 | call newsedim(ngrid, nlay, ngrid*nlay, ngrid*nlay,microtimestep, pplev, masse, epaisseur, ztsed, & |
|---|
| 867 | rsedcloudco2, rhocloudco2t, zqsed(:,:,igcm_ccnco2_h2o_mass_ice), wq, beta) |
|---|
| 868 | do ig = 1, ngrid |
|---|
| 869 | sum_subpdqs_sedccn(ig,igcm_ccnco2_h2o_mass_ice) = sum_subpdqs_sedccn(ig,igcm_ccnco2_h2o_mass_ice) + & |
|---|
| 870 | wq(ig,1) / microtimestep |
|---|
| 871 | end do |
|---|
| 872 | |
|---|
| 873 | wq(:,:) = 0. |
|---|
| 874 | ! for ccnco2_h2o_mass_ccn |
|---|
| 875 | call newsedim(ngrid, nlay, ngrid*nlay, ngrid*nlay,microtimestep, pplev, masse, epaisseur, ztsed, & |
|---|
| 876 | rsedcloudco2, rhocloudco2t, zqsed(:,:,igcm_ccnco2_h2o_mass_ccn), wq, beta) |
|---|
| 877 | do ig = 1, ngrid |
|---|
| 878 | sum_subpdqs_sedccn(ig,igcm_ccnco2_h2o_mass_ccn) = sum_subpdqs_sedccn(ig,igcm_ccnco2_h2o_mass_ccn) + & |
|---|
| 879 | wq(ig,1) / microtimestep |
|---|
| 880 | end do |
|---|
| 881 | |
|---|
| 882 | wq(:,:) = 0. |
|---|
| 883 | ! for ccnco2_h2o_number |
|---|
| 884 | call newsedim(ngrid, nlay, ngrid*nlay, ngrid*nlay,microtimestep, pplev, masse, epaisseur, ztsed, & |
|---|
| 885 | rsedcloudco2, rhocloudco2t, zqsed(:,:,igcm_ccnco2_h2o_number), wq, beta) |
|---|
| 886 | do ig = 1, ngrid |
|---|
| 887 | sum_subpdqs_sedccn(ig,igcm_ccnco2_h2o_number) = sum_subpdqs_sedccn(ig,igcm_ccnco2_h2o_number) + & |
|---|
| 888 | wq(ig,1) / microtimestep |
|---|
| 889 | end do |
|---|
| 890 | end if |
|---|
| 891 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 892 | ! 4.5.e. Compute tendencies due to the sedimation process |
|---|
| 893 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 894 | do l = 1, nlay |
|---|
| 895 | do ig = 1, ngrid |
|---|
| 896 | subpdqsed(ig,l,igcm_ccnco2_mass) = ( zqsed(ig,l,igcm_ccnco2_mass) - zqsed0(ig,l,igcm_ccnco2_mass) ) & |
|---|
| 897 | / microtimestep |
|---|
| 898 | |
|---|
| 899 | subpdqsed(ig,l,igcm_ccnco2_number) = ( zqsed(ig,l,igcm_ccnco2_number) - zqsed0(ig,l,igcm_ccnco2_number) )& |
|---|
| 900 | / microtimestep |
|---|
| 901 | |
|---|
| 902 | subpdqsed(ig,l,igcm_co2_ice) = ( zqsed(ig,l,igcm_co2_ice) - zqsed0(ig,l,igcm_co2_ice) ) / microtimestep |
|---|
| 903 | |
|---|
| 904 | if (meteo_flux) then |
|---|
| 905 | subpdqsed(ig,l,igcm_ccnco2_meteor_mass) = ( zqsed(ig,l,igcm_ccnco2_meteor_mass) - & |
|---|
| 906 | zqsed0(ig,l,igcm_ccnco2_meteor_mass) ) / microtimestep |
|---|
| 907 | |
|---|
| 908 | subpdqsed(ig,l,igcm_ccnco2_meteor_number) = ( zqsed(ig,l,igcm_ccnco2_meteor_number) - & |
|---|
| 909 | zqsed0(ig,l,igcm_ccnco2_meteor_number) ) / microtimestep |
|---|
| 910 | end if |
|---|
| 911 | if (co2useh2o) then |
|---|
| 912 | subpdqsed(ig,l,igcm_ccnco2_h2o_number) = ( zqsed(ig,l,igcm_ccnco2_h2o_number) - & |
|---|
| 913 | zqsed0(ig,l,igcm_ccnco2_h2o_number) ) / microtimestep |
|---|
| 914 | |
|---|
| 915 | subpdqsed(ig,l,igcm_ccnco2_h2o_mass_ice) = ( zqsed(ig,l,igcm_ccnco2_h2o_mass_ice) - & |
|---|
| 916 | zqsed0(ig,l,igcm_ccnco2_h2o_mass_ice) ) / microtimestep |
|---|
| 917 | |
|---|
| 918 | subpdqsed(ig,l,igcm_ccnco2_h2o_mass_ccn) = ( zqsed(ig,l,igcm_ccnco2_h2o_mass_ccn) - & |
|---|
| 919 | zqsed0(ig,l,igcm_ccnco2_h2o_mass_ccn) ) / microtimestep |
|---|
| 920 | end if |
|---|
| 921 | end do |
|---|
| 922 | end do |
|---|
| 923 | ! update subtimestep tendencies with sedimentation input |
|---|
| 924 | do l = 1, nlay |
|---|
| 925 | do ig = 1, ngrid |
|---|
| 926 | sum_subpdq(ig,l,igcm_ccnco2_mass) = sum_subpdq(ig,l,igcm_ccnco2_mass) + subpdqsed(ig,l,igcm_ccnco2_mass) |
|---|
| 927 | |
|---|
| 928 | sum_subpdq(ig,l,igcm_ccnco2_number) = sum_subpdq(ig,l,igcm_ccnco2_number) + subpdqsed(ig,l,igcm_ccnco2_number) |
|---|
| 929 | |
|---|
| 930 | sum_subpdq(ig,l,igcm_co2_ice) = sum_subpdq(ig,l,igcm_co2_ice) + subpdqsed(ig,l,igcm_co2_ice) |
|---|
| 931 | if (meteo_flux) then |
|---|
| 932 | sum_subpdq(ig,l,igcm_ccnco2_meteor_mass) = sum_subpdq(ig,l,igcm_ccnco2_meteor_mass) + & |
|---|
| 933 | subpdqsed(ig,l,igcm_ccnco2_meteor_mass) |
|---|
| 934 | |
|---|
| 935 | sum_subpdq(ig,l,igcm_ccnco2_meteor_number) = sum_subpdq(ig,l,igcm_ccnco2_meteor_number) + & |
|---|
| 936 | subpdqsed(ig,l,igcm_ccnco2_meteor_number) |
|---|
| 937 | end if |
|---|
| 938 | if (co2useh2o) then |
|---|
| 939 | sum_subpdq(ig,l,igcm_ccnco2_h2o_mass_ice) = sum_subpdq(ig,l,igcm_ccnco2_h2o_mass_ice) + & |
|---|
| 940 | subpdqsed(ig,l,igcm_ccnco2_h2o_mass_ice) |
|---|
| 941 | |
|---|
| 942 | sum_subpdq(ig,l,igcm_ccnco2_h2o_mass_ccn) = sum_subpdq(ig,l,igcm_ccnco2_h2o_mass_ccn) + & |
|---|
| 943 | subpdqsed(ig,l,igcm_ccnco2_h2o_mass_ccn) |
|---|
| 944 | |
|---|
| 945 | sum_subpdq(ig,l,igcm_ccnco2_h2o_number) = sum_subpdq(ig,l,igcm_ccnco2_h2o_number) + & |
|---|
| 946 | subpdqsed(ig,l,igcm_ccnco2_h2o_number) |
|---|
| 947 | end if |
|---|
| 948 | end do |
|---|
| 949 | end do |
|---|
| 950 | end if !(end if sedimentation) |
|---|
| 951 | end do ! of do microstep = 1, imicroco2 |
|---|
| 952 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 953 | ! 5. Compute final tendencies after time loop |
|---|
| 954 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 955 | ! condensation/sublimation rate in the atmosphere (kg.kg-1.s-1) |
|---|
| 956 | do l = nlay, 1, -1 |
|---|
| 957 | do ig = 1, ngrid |
|---|
| 958 | pcondicea(ig,l) = sum_subpdq(ig,l,igcm_co2_ice) / real(imicroco2) |
|---|
| 959 | end do |
|---|
| 960 | end do |
|---|
| 961 | |
|---|
| 962 | ! CO2 flux at surface (kg.m-2.s-1) |
|---|
| 963 | do ig = 1, ngrid |
|---|
| 964 | pdqs_sedco2(ig) = sum_subpdqs_sedco2(ig) / real(imicroco2) |
|---|
| 965 | pdqs_sedccn(ig,:) = sum_subpdqs_sedccn(ig,:) / real(imicroco2) |
|---|
| 966 | end do |
|---|
| 967 | |
|---|
| 968 | ! temperature tendency (T.s-1) |
|---|
| 969 | do l = 1, nlay |
|---|
| 970 | do ig = 1, ngrid |
|---|
| 971 | pdtcloudco2(ig,l) = ( sum_subpdt(ig,l)/real(imicroco2) ) - pdt(ig,l) |
|---|
| 972 | if (pdtcloudco2(ig,l) /= pdtcloudco2(ig,l)) then |
|---|
| 973 | write(*,*)ig,l,pdtcloudco2(ig,l), pdt(ig,l), sum_subpdt(ig,l) |
|---|
| 974 | call abort_physic('co2clouds', 'ptdcloudco2 is NaN', 1) |
|---|
| 975 | end if |
|---|
| 976 | end do |
|---|
| 977 | end do |
|---|
| 978 | |
|---|
| 979 | ! tracers tendencies |
|---|
| 980 | do l = 1, nlay |
|---|
| 981 | do ig = 1, ngrid |
|---|
| 982 | pdqcloudco2(ig,l,igcm_co2) = 0. ! here is the trick, this tendency is computed in co2condens |
|---|
| 983 | |
|---|
| 984 | pdqcloudco2(ig,l,igcm_co2_ice) = ( sum_subpdq(ig,l,igcm_co2_ice) / real(imicroco2) ) - pdq(ig,l,igcm_co2_ice) |
|---|
| 985 | |
|---|
| 986 | pdqcloudco2(ig,l,igcm_ccnco2_mass) = ( sum_subpdq(ig,l,igcm_ccnco2_mass)/real(imicroco2) ) - & |
|---|
| 987 | pdq(ig,l,igcm_ccnco2_mass) |
|---|
| 988 | |
|---|
| 989 | pdqcloudco2(ig,l,igcm_ccnco2_number) = ( sum_subpdq(ig,l,igcm_ccnco2_number) / real(imicroco2) ) - & |
|---|
| 990 | pdq(ig,l,igcm_ccnco2_number) |
|---|
| 991 | |
|---|
| 992 | pdqcloudco2(ig,l,igcm_dust_mass) = ( sum_subpdq(ig,l,igcm_dust_mass) / real(imicroco2) ) - & |
|---|
| 993 | pdq(ig,l,igcm_dust_mass) |
|---|
| 994 | |
|---|
| 995 | pdqcloudco2(ig,l,igcm_dust_number) = ( sum_subpdq(ig,l,igcm_dust_number)/real(imicroco2) ) - & |
|---|
| 996 | pdq(ig,l,igcm_dust_number) |
|---|
| 997 | |
|---|
| 998 | if (meteo_flux) then |
|---|
| 999 | pdqcloudco2(ig,l,igcm_ccnco2_meteor_mass) = ( sum_subpdq(ig,l,igcm_ccnco2_meteor_mass)/real(imicroco2) ) & |
|---|
| 1000 | - pdq(ig,l,igcm_ccnco2_meteor_mass) |
|---|
| 1001 | |
|---|
| 1002 | pdqcloudco2(ig,l,igcm_ccnco2_meteor_number) = ( sum_subpdq(ig,l,igcm_ccnco2_meteor_number) / real(imicroco2) ) & |
|---|
| 1003 | - pdq(ig,l,igcm_ccnco2_meteor_number) |
|---|
| 1004 | end if |
|---|
| 1005 | if (co2useh2o) then |
|---|
| 1006 | pdqcloudco2(ig,l,igcm_h2o_ice) = ( sum_subpdq(ig,l,igcm_h2o_ice) / real(imicroco2) ) - & |
|---|
| 1007 | pdq(ig,l,igcm_h2o_ice) |
|---|
| 1008 | |
|---|
| 1009 | pdqcloudco2(ig,l,igcm_ccn_mass) = ( sum_subpdq(ig,l,igcm_ccn_mass) / real(imicroco2) ) - & |
|---|
| 1010 | pdq(ig,l,igcm_ccn_mass) |
|---|
| 1011 | |
|---|
| 1012 | pdqcloudco2(ig,l,igcm_ccn_number) = ( sum_subpdq(ig,l,igcm_ccn_number) / real(imicroco2) ) - & |
|---|
| 1013 | pdq(ig,l,igcm_ccn_number) |
|---|
| 1014 | |
|---|
| 1015 | pdqcloudco2(ig,l,igcm_ccnco2_h2o_number) = ( sum_subpdq(ig,l,igcm_ccnco2_h2o_number) / real(imicroco2) ) - & |
|---|
| 1016 | pdq(ig,l,igcm_ccnco2_h2o_number) |
|---|
| 1017 | |
|---|
| 1018 | pdqcloudco2(ig,l,igcm_ccnco2_h2o_mass_ice) = ( sum_subpdq(ig,l,igcm_ccnco2_h2o_mass_ice) / real(imicroco2) ) - & |
|---|
| 1019 | pdq(ig,l,igcm_ccnco2_h2o_mass_ice) |
|---|
| 1020 | |
|---|
| 1021 | pdqcloudco2(ig,l,igcm_ccnco2_h2o_mass_ccn) = ( sum_subpdq(ig,l,igcm_ccnco2_h2o_mass_ccn) / real(imicroco2) )- & |
|---|
| 1022 | pdq(ig,l,igcm_ccnco2_h2o_mass_ccn) |
|---|
| 1023 | |
|---|
| 1024 | end if |
|---|
| 1025 | end do ! ngrid |
|---|
| 1026 | end do ! nlay |
|---|
| 1027 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 1028 | ! 6. Update clouds physical values in the cloud (for output) |
|---|
| 1029 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 1030 | ! 6.1. Update density of co2 ice, riceco2 and opacity |
|---|
| 1031 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 1032 | do l = 1, nlay |
|---|
| 1033 | do ig = 1, ngrid |
|---|
| 1034 | Niceco2 = pqeff(ig,l,igcm_co2_ice) + ( pdq(ig,l,igcm_co2_ice) + pdqcloudco2(ig,l,igcm_co2_ice) ) * ptimestep |
|---|
| 1035 | Niceco2 = max(Niceco2, threshold) |
|---|
| 1036 | |
|---|
| 1037 | ! meteoritic particles are considered like dust, => rho_dust |
|---|
| 1038 | Nccnco2 = max( (pqeff(ig,l,igcm_ccnco2_number) + (pdq(ig,l,igcm_ccnco2_number) + & |
|---|
| 1039 | pdqcloudco2(ig,l, igcm_ccnco2_number)) * ptimestep) & |
|---|
| 1040 | , threshold) |
|---|
| 1041 | |
|---|
| 1042 | Qccnco2 = max( (pqeff(ig,l,igcm_ccnco2_mass) + (pdq(ig,l,igcm_ccnco2_mass) + & |
|---|
| 1043 | pdqcloudco2(ig,l, igcm_ccnco2_mass)) * ptimestep)& |
|---|
| 1044 | , threshold) |
|---|
| 1045 | |
|---|
| 1046 | myT = pteff(ig,l) + (pdt(ig,l)+pdtcloudco2(ig,l))*ptimestep |
|---|
| 1047 | |
|---|
| 1048 | Nccnco2_h2o = 0. |
|---|
| 1049 | Qccnco2_h2o = 0. |
|---|
| 1050 | if (co2useh2o) then |
|---|
| 1051 | Nccnco2_h2o = max( (pqeff(ig,l,igcm_ccnco2_h2o_number) + (pdq(ig,l,igcm_ccnco2_h2o_number) + & |
|---|
| 1052 | pdqcloudco2(ig,l, igcm_ccnco2_h2o_number)) * ptimestep) & |
|---|
| 1053 | , threshold) |
|---|
| 1054 | |
|---|
| 1055 | Qccnco2_h2o = max( (pqeff(ig,l,igcm_ccnco2_h2o_mass_ice) + pqeff(ig,l,igcm_ccnco2_h2o_mass_ccn) + & |
|---|
| 1056 | (pdq(ig,l,igcm_ccnco2_h2o_mass_ice) + pdq(ig,l,igcm_ccnco2_h2o_mass_ccn) + & |
|---|
| 1057 | pdqcloudco2(ig,l,igcm_ccnco2_h2o_mass_ice) + pdqcloudco2(ig,l,igcm_ccnco2_h2o_mass_ccn)) * & |
|---|
| 1058 | ptimestep)& |
|---|
| 1059 | , threshold) |
|---|
| 1060 | Nccnco2 = Nccnco2 - Nccnco2_h2o |
|---|
| 1061 | Qccnco2 = Qccnco2 - Qccnco2_h2o |
|---|
| 1062 | if (Nccnco2 <= 0) then |
|---|
| 1063 | Nccnco2 = threshold |
|---|
| 1064 | Qccnco2 = threshold |
|---|
| 1065 | end if |
|---|
| 1066 | end if |
|---|
| 1067 | |
|---|
| 1068 | ! Compute particle size |
|---|
| 1069 | call updaterice_microco2(dble(Niceco2), dble(Qccnco2), dble(Nccnco2), dble(Qccnco2_h2o), dble(Nccnco2_h2o), myT, & |
|---|
| 1070 | tauscaling(ig), riceco2(ig,l), rhocloudco2(ig,l)) |
|---|
| 1071 | |
|---|
| 1072 | ! Compute opacities |
|---|
| 1073 | if ( (Niceco2 >= threshold_2 .and. (Nccnco2+Nccnco2_h2o)*tauscaling(ig) >= 1.) ) then |
|---|
| 1074 | n_derf = derf( (rb_cldco2(1)-log(riceco2(ig,l))) *dev2) |
|---|
| 1075 | Qext1bins2(ig,l) = 0. |
|---|
| 1076 | do i = 1, nbinco2_cld |
|---|
| 1077 | n_aer(i) = -0.5 * (Nccnco2+Nccnco2_h2o)*tauscaling(ig) * n_derf |
|---|
| 1078 | |
|---|
| 1079 | n_derf = derf((rb_cldco2(i+1)-log(riceco2(ig,l))) *dev2) |
|---|
| 1080 | n_aer(i) = n_aer(i) + (0.5 * (Nccnco2+Nccnco2_h2o)*tauscaling(ig) * n_derf) |
|---|
| 1081 | |
|---|
| 1082 | Qext1bins2(ig,l) = Qext1bins2(ig,l) + Qext1bins(i) * n_aer(i) |
|---|
| 1083 | end do |
|---|
| 1084 | else |
|---|
| 1085 | riceco2(ig,l) = 0. |
|---|
| 1086 | Qext1bins2(ig,l) = 0. |
|---|
| 1087 | end if |
|---|
| 1088 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 1089 | ! 6.2. Update rice and rdust |
|---|
| 1090 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 1091 | ! update rice water only if co2 use h2o ice as CCN |
|---|
| 1092 | if (co2useh2o) then |
|---|
| 1093 | call updaterice_micro( & |
|---|
| 1094 | pqeff(ig,l,igcm_h2o_ice) + (pdq(ig,l,igcm_h2o_ice) + pdqcloudco2(ig,l,igcm_h2o_ice))*ptimestep, & |
|---|
| 1095 | pqeff(ig,l,igcm_ccn_mass) + (pdq(ig,l,igcm_ccn_mass) + pdqcloudco2(ig,l,igcm_ccn_mass))*ptimestep, & |
|---|
| 1096 | pqeff(ig,l,igcm_ccn_number) + (pdq(ig,l,igcm_ccn_number) + pdqcloudco2(ig,l,igcm_ccn_number))*ptimestep, & |
|---|
| 1097 | tauscaling(ig),rice(ig,l),rhocloud(ig,l)) |
|---|
| 1098 | end if |
|---|
| 1099 | |
|---|
| 1100 | ! update rdust |
|---|
| 1101 | call updaterdust( & |
|---|
| 1102 | pqeff(ig,l,igcm_dust_mass) + (pdq(ig,l,igcm_dust_mass) + pdqcloudco2(ig,l,igcm_dust_mass))*ptimestep, & |
|---|
| 1103 | pqeff(ig,l,igcm_dust_number) + (pdq(ig,l,igcm_dust_number) + pdqcloudco2(ig,l,igcm_dust_number))*ptimestep, & |
|---|
| 1104 | rdust(ig,l)) |
|---|
| 1105 | end do ! ngrid |
|---|
| 1106 | end do ! nlay |
|---|
| 1107 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 1108 | ! 7. A correction if a lot of subliming CO2 fills the 1st layer FF (04/2005). Then that should not affect the ice |
|---|
| 1109 | ! particle radius |
|---|
| 1110 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 1111 | do ig = 1, ngrid |
|---|
| 1112 | if ( pdpsrf(ig)*ptimestep > 0.9*(pplev(ig,1)-pplev(ig,2))) then |
|---|
| 1113 | |
|---|
| 1114 | if ( pdpsrf(ig)*ptimestep > 0.9*(pplev(ig,1)-pplev(ig,3)) ) then |
|---|
| 1115 | riceco2(ig,2) = riceco2(ig,3) |
|---|
| 1116 | end if |
|---|
| 1117 | |
|---|
| 1118 | riceco2(ig,1) = riceco2(ig,2) |
|---|
| 1119 | end if |
|---|
| 1120 | end do |
|---|
| 1121 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 1122 | ! 8. CO2 saturated quantities |
|---|
| 1123 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 1124 | ! 8.1 Compute CO2 saturation in layers |
|---|
| 1125 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 1126 | call co2sat(ngrid*nlay, pteff+(pdt+pdtcloudco2)*ptimestep, zqsatco2) |
|---|
| 1127 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 1128 | ! 8.2 Compute CO2 saturated quantities in layers |
|---|
| 1129 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 1130 | do l = 1, nlay |
|---|
| 1131 | do ig = 1, ngrid |
|---|
| 1132 | satuco2(ig,l) = ( pqeff(ig,l,igcm_co2) + (pdq(ig,l,igcm_co2) + pdqcloudco2(ig,l,igcm_co2))*ptimestep ) * & |
|---|
| 1133 | (mmean(ig,l)/(mco2*1e3)) * pplay(ig,l) / zqsatco2(ig,l) |
|---|
| 1134 | end do |
|---|
| 1135 | end do |
|---|
| 1136 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 1137 | ! 9. Everything modified by CO2 microphysics must be wrt co2cloudfrac |
|---|
| 1138 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 1139 | if (CLFvaryingCO2) then |
|---|
| 1140 | do l = 1, nlay |
|---|
| 1141 | do ig = 1, ngrid |
|---|
| 1142 | pdqcloudco2(ig,l,igcm_ccnco2_mass) = pdqcloudco2(ig,l,igcm_ccnco2_mass) * co2cloudfrac(ig,l) |
|---|
| 1143 | |
|---|
| 1144 | pdqcloudco2(ig,l,igcm_ccnco2_number) = pdqcloudco2(ig,l,igcm_ccnco2_number) * co2cloudfrac(ig,l) |
|---|
| 1145 | |
|---|
| 1146 | pdqcloudco2(ig,l,igcm_dust_mass) = pdqcloudco2(ig,l,igcm_dust_mass) * co2cloudfrac(ig,l) |
|---|
| 1147 | |
|---|
| 1148 | pdqcloudco2(ig,l,igcm_dust_number) = pdqcloudco2(ig,l,igcm_dust_number) * co2cloudfrac(ig,l) |
|---|
| 1149 | |
|---|
| 1150 | pdqcloudco2(ig,l,igcm_co2_ice) = pdqcloudco2(ig,l,igcm_co2_ice) * co2cloudfrac(ig,l) |
|---|
| 1151 | |
|---|
| 1152 | pdqcloudco2(ig,l,igcm_co2) = pdqcloudco2(ig,l,igcm_co2) * co2cloudfrac(ig,l) |
|---|
| 1153 | |
|---|
| 1154 | pdtcloudco2(ig,l) = pdtcloudco2(ig,l) * co2cloudfrac(ig,l) |
|---|
| 1155 | |
|---|
| 1156 | Qext1bins2(ig,l) = Qext1bins2(ig,l) * co2cloudfrac(ig,l) |
|---|
| 1157 | |
|---|
| 1158 | if (meteo_flux) then |
|---|
| 1159 | pdqcloudco2(ig,l,igcm_ccnco2_meteor_mass) = pdqcloudco2(ig,l,igcm_ccnco2_meteor_mass) * co2cloudfrac(ig,l) |
|---|
| 1160 | |
|---|
| 1161 | pdqcloudco2(ig,l,igcm_ccnco2_meteor_number) = pdqcloudco2(ig,l,igcm_ccnco2_meteor_number) * co2cloudfrac(ig,l) |
|---|
| 1162 | end if |
|---|
| 1163 | if (co2useh2o) then |
|---|
| 1164 | pdqcloudco2(ig,l,igcm_h2o_ice) = pdqcloudco2(ig,l,igcm_h2o_ice) * co2cloudfrac(ig,l) |
|---|
| 1165 | |
|---|
| 1166 | pdqcloudco2(ig,l,igcm_ccn_mass) = pdqcloudco2(ig,l,igcm_ccn_mass) * co2cloudfrac(ig,l) |
|---|
| 1167 | |
|---|
| 1168 | pdqcloudco2(ig,l,igcm_ccn_number) = pdqcloudco2(ig,l,igcm_ccn_number) * co2cloudfrac(ig,l) |
|---|
| 1169 | |
|---|
| 1170 | pdqcloudco2(ig,l,igcm_ccnco2_h2o_mass_ice) = pdqcloudco2(ig,l,igcm_ccnco2_h2o_mass_ice) * co2cloudfrac(ig,l) |
|---|
| 1171 | |
|---|
| 1172 | pdqcloudco2(ig,l,igcm_ccnco2_h2o_mass_ccn) = pdqcloudco2(ig,l,igcm_ccnco2_h2o_mass_ccn) * co2cloudfrac(ig,l) |
|---|
| 1173 | |
|---|
| 1174 | pdqcloudco2(ig,l,igcm_ccnco2_h2o_number) = pdqcloudco2(ig,l,igcm_ccnco2_h2o_number) * co2cloudfrac(ig,l) |
|---|
| 1175 | end if |
|---|
| 1176 | end do ! ngrid |
|---|
| 1177 | end do ! nlay |
|---|
| 1178 | end if ! if CLFvaryingCO2 is true |
|---|
| 1179 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 1180 | ! 10. Compute opacity at 1 micron: Opacity in mesh ig is the sum over l of Qext1bins2. Is this true ? |
|---|
| 1181 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 1182 | tau1mic(:)=0. |
|---|
| 1183 | do l = 1, nlay |
|---|
| 1184 | do ig = 1, ngrid |
|---|
| 1185 | tau1mic(ig) = tau1mic(ig) + Qext1bins2(ig,l) |
|---|
| 1186 | end do |
|---|
| 1187 | end do |
|---|
| 1188 | #ifndef MESOSCALE |
|---|
| 1189 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 1190 | ! 11. Write outputs in diagfi.nc |
|---|
| 1191 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 1192 | call write_output( "satuco2", "vap in satu", " ", satuco2(:,:)) |
|---|
| 1193 | |
|---|
| 1194 | call write_output( "precip_co2_ice_rate", "surface deposition rate of co2 ice", "kg.m-2.s-1", pdqs_sedco2(:)) |
|---|
| 1195 | |
|---|
| 1196 | call write_output( "co2ice_cond_rate", "CO2 condensation rate in atm layers", "kg.kg-1.s-1", pcondicea(:,:)) |
|---|
| 1197 | |
|---|
| 1198 | call write_output( "pdtcloudco2", "temperature variation of CO2 latent heat", "K.s-1", pdtcloudco2(:,:)) |
|---|
| 1199 | |
|---|
| 1200 | call write_output( "riceco2", "ice radius", "m", riceco2(:,:)) |
|---|
| 1201 | |
|---|
| 1202 | call write_output( "Tau3D1mic", " co2 ice opacities", " ", Qext1bins2(:,:)) |
|---|
| 1203 | |
|---|
| 1204 | call write_output( "tau1mic", "co2 ice opacity 1 micron", " ", tau1mic(:)) |
|---|
| 1205 | |
|---|
| 1206 | if (CLFvaryingCO2) then |
|---|
| 1207 | call write_output( "co2cloudfrac", "co2 cloud fraction", " ", co2cloudfrac(:,:)) |
|---|
| 1208 | end if |
|---|
| 1209 | #endif |
|---|
| 1210 | !======================================================================================================================! |
|---|
| 1211 | ! END =================================================================================================================! |
|---|
| 1212 | !======================================================================================================================! |
|---|
| 1213 | end subroutine co2cloud |
|---|
| 1214 | |
|---|
| 1215 | end module co2cloud_mod |
|---|
