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