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