| 1 | !======================================================================================================================! |
|---|
| 2 | ! Module: CO2 condensation for the CO2 cloud microphysics =============================================================! |
|---|
| 3 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 4 | ! Authors: Christophe Mathé, Anni Määttänen |
|---|
| 5 | ! Date: 16/04/2020 |
|---|
| 6 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 7 | ! Contains subroutines: |
|---|
| 8 | ! - co2condens4micro: condensation/sublimation of CO2 ice on the ground and compute pressure change resulting |
|---|
| 9 | ! |
|---|
| 10 | ! - vl1d: Van-Leer scheme |
|---|
| 11 | !======================================================================================================================! |
|---|
| 12 | module co2condens_mod4micro |
|---|
| 13 | |
|---|
| 14 | implicit none |
|---|
| 15 | |
|---|
| 16 | contains |
|---|
| 17 | !======================================================================================================================! |
|---|
| 18 | ! SUBROUTINE: co2condens4micro ========================================================================================! |
|---|
| 19 | !======================================================================================================================! |
|---|
| 20 | ! Subject: |
|---|
| 21 | !--------- |
|---|
| 22 | ! Condensation/sublimation of CO2 ice on the ground and compute pressure change resulting |
|---|
| 23 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 24 | ! Comments: |
|---|
| 25 | !---------- |
|---|
| 26 | ! Adapted from co2condens_mod.F |
|---|
| 27 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 28 | ! Paper: |
|---|
| 29 | !------- |
|---|
| 30 | ! Forget et al. (2008), "Non condensable gas enrichment and depletion in the Martian polar regions." |
|---|
| 31 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 32 | ! Algorithm: |
|---|
| 33 | !----------- |
|---|
| 34 | ! 1. Initialization |
|---|
| 35 | ! 2. Firstcall |
|---|
| 36 | ! 3. Compute CO2 Volume mixing ratio |
|---|
| 37 | ! 4. Set zcondicea, zfallice from co2clouds condensation rate and set zt |
|---|
| 38 | ! 5. Main co2condens |
|---|
| 39 | ! 5.1. Forecast of ground temperature ztsrf and frost temperature ztcondsol |
|---|
| 40 | ! 5.2. Check if we have condensation/sublimation on the ground |
|---|
| 41 | ! 5.3. Compute zfallheat |
|---|
| 42 | ! 5.4. Compute direct condensation/sublimation of CO2 ice |
|---|
| 43 | ! 5.4.a. If there is not enough CO2 tracer in 1st layer to condense |
|---|
| 44 | ! 5.4.b. If the entire CO2 ice layer sublimes (including what has just condensed in the atmosphere) |
|---|
| 45 | ! 5.5. Changing CO2 ice amount and pressure |
|---|
| 46 | ! 5.6. Surface albedo and emissivity of the surface below the snow (emisref) |
|---|
| 47 | ! 5.6.a. Check that amont of CO2 ice is not problematic |
|---|
| 48 | ! 5.6.b. Set albedo and emissivity of the surface |
|---|
| 49 | ! 5.6.c. Set pemisurf to emissiv when there is bare surface (needed for co2snow) |
|---|
| 50 | ! 5.7. Correction to account for redistribution between sigma or hybrid layers when changing surface pressure (and |
|---|
| 51 | ! warming/cooling of the CO2 currently changing phase). |
|---|
| 52 | ! 5.7.a. Mass fluxes through the sigma levels (kg.m-2.s-1) (>0 when up) |
|---|
| 53 | ! 5.7.b. Mass of each layer at the end of timestep |
|---|
| 54 | ! 5.7.c. Corresponding fluxes for T, U, V and Q (averaging operator for transport) |
|---|
| 55 | ! 5.7.c.i. Value transfert at the surface interface when condensation/sublimation |
|---|
| 56 | ! 5.7.c.ii. Van Leer scheme |
|---|
| 57 | ! 5.7.c.iii Compute tendencies on T, U, V, Q |
|---|
| 58 | ! 6. CO2 snow / clouds scheme |
|---|
| 59 | ! 7. Extra special case for surface temperature tendency pdtsrfc |
|---|
| 60 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 61 | subroutine co2condens4micro(ngrid, nlayer, nq, ptimestep, pcapcal, pplay, pplev, ptsrf, pt, pphi, pdt, pdu, pdv, & |
|---|
| 62 | pdtsrf, pu, pv, pq, pdq, piceco2, psolaralb, pemisurf, pdtc, pdtsrfc, pdpsrf, pduc, & |
|---|
| 63 | pdvc, pdqc, fluxsurf_sw, zls, zdqssed_co2, pcondicea_co2microp, zdtcloudco2) |
|---|
| 64 | |
|---|
| 65 | use tracer_mod, only: noms, igcm_co2, igcm_co2_ice, igcm_h2o_vap, igcm_h2o_ice, igcm_dust_number, igcm_dust_mass, & |
|---|
| 66 | igcm_ccnco2_number, igcm_ccnco2_mass, igcm_ccn_mass |
|---|
| 67 | |
|---|
| 68 | use surfdat_h, only: emissiv, phisfi |
|---|
| 69 | |
|---|
| 70 | use geometry_mod, only: latitude, longitude_deg, latitude_deg |
|---|
| 71 | |
|---|
| 72 | use planete_h, only: obliquit |
|---|
| 73 | |
|---|
| 74 | use comcstfi_h, only: cpp, g, r, pi |
|---|
| 75 | |
|---|
| 76 | #ifndef MESOSCALE |
|---|
| 77 | use vertical_layers_mod, only: ap, bp |
|---|
| 78 | #endif |
|---|
| 79 | |
|---|
| 80 | implicit none |
|---|
| 81 | |
|---|
| 82 | include "callkeys.h" |
|---|
| 83 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 84 | ! VARIABLE DECLARATION |
|---|
| 85 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 86 | ! Input arguments: |
|---|
| 87 | !----------------- |
|---|
| 88 | integer, intent(in) :: & |
|---|
| 89 | nq, &! number of tracers |
|---|
| 90 | ngrid, &! number of atmospheric columns |
|---|
| 91 | nlayer ! number of vertical layers |
|---|
| 92 | |
|---|
| 93 | real, intent(in) :: & |
|---|
| 94 | ptimestep, &! physics timestep (s) |
|---|
| 95 | pcapcal(ngrid), &! surface specific heat |
|---|
| 96 | pplay(ngrid,nlayer), &! mid-layer pressure (Pa) |
|---|
| 97 | pplev(ngrid,nlayer+1), &! inter-layer pressure (Pa) |
|---|
| 98 | ptsrf(ngrid), &! surface temperature (K) |
|---|
| 99 | pphi(ngrid,nlayer), &! geopotential (m2.s-2) |
|---|
| 100 | pt(ngrid,nlayer), &! atmospheric temperature (K) |
|---|
| 101 | pu(ngrid,nlayer), &! zonal wind (m/s) |
|---|
| 102 | pv(ngrid,nlayer), &! meridional wind (m/s) |
|---|
| 103 | pdt(ngrid,nlayer), &! tendency on temperature from previous physical processes (K/s) |
|---|
| 104 | pdu(ngrid,nlayer), &! tendency on zonal wind from previous physical processes(m/s2) |
|---|
| 105 | pdv(ngrid,nlayer), &! tendency on meridional wind from previous physical processes (m/s2) |
|---|
| 106 | pdtsrf(ngrid), &! tendency on surface temperature from previous physical processes (K/s) |
|---|
| 107 | pq(ngrid,nlayer,nq), &! tracers (../kg_air) |
|---|
| 108 | pdq(ngrid,nlayer,nq), &! tendency on tracers from previous physical processes |
|---|
| 109 | zdqssed_co2(ngrid), &! CO2 flux at the surface (kg.m-2.s-1) |
|---|
| 110 | zdtcloudco2(ngrid,nlayer), &! tendency on temperature due to latent heat |
|---|
| 111 | fluxsurf_sw(ngrid,2), &! added to calculate flux dependent albedo |
|---|
| 112 | zls, &! solar longitude (rad) |
|---|
| 113 | pcondicea_co2microp(ngrid,nlayer) ! tendency due to CO2 condensation (kg/kg.s-1) |
|---|
| 114 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 115 | ! In/output arguments: |
|---|
| 116 | !--------------------- |
|---|
| 117 | real, intent(inout) :: & |
|---|
| 118 | piceco2(ngrid), &! CO2 ice on the surface (kg.m-2) |
|---|
| 119 | pemisurf(ngrid), &! emissivity of the surface |
|---|
| 120 | psolaralb(ngrid,2) ! albedo of the surface |
|---|
| 121 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 122 | ! Output arguments: |
|---|
| 123 | !------------------ |
|---|
| 124 | real, intent(out) :: & |
|---|
| 125 | pdtc(ngrid,nlayer), &! tendency on temperature dT/dt due to cond/sub (K/s) |
|---|
| 126 | pdtsrfc(ngrid), &! tendency on surface temperature (K/s) |
|---|
| 127 | pdpsrf(ngrid), &! tendency on surface pressure (Pa/s) |
|---|
| 128 | pduc(ngrid,nlayer), &! tendency on zonal wind (m.s-2) |
|---|
| 129 | pdvc(ngrid,nlayer), &! tendency on meridional wind (m.s-2) |
|---|
| 130 | pdqc(ngrid,nlayer,nq) ! tendency on tracers |
|---|
| 131 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 132 | ! Local: |
|---|
| 133 | !------- |
|---|
| 134 | !----1) Parameters: |
|---|
| 135 | !------------------ |
|---|
| 136 | real, parameter :: & |
|---|
| 137 | latcond = 595594, &! latent heat of solid CO2 ice (J/kg) |
|---|
| 138 | cpice = 1000., &! specific heat of CO2 ice (J.kg-1.K-1) |
|---|
| 139 | tcond1mb = 136.27,&! condensation temperature at 1 mbar (K) |
|---|
| 140 | m_co2 = 44.01E-3, &! CO2 molecular mass (kg/mol) |
|---|
| 141 | m_noco2 = 33.37E-3 ! non condensible molecular mass (kg/mol) |
|---|
| 142 | |
|---|
| 143 | logical, parameter :: & |
|---|
| 144 | improved_ztcond = .true. ! improved_ztcond flag: If set to .true. (AND running with a 'co2' tracer) then |
|---|
| 145 | ! condensation temperature is computed using partial pressure of CO2 |
|---|
| 146 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 147 | !----2) Saved: |
|---|
| 148 | !------------- |
|---|
| 149 | real, save :: & |
|---|
| 150 | A, &! coefficient used to compute mean molecular mass |
|---|
| 151 | B, &! coefficient used to compute mean molecular mass |
|---|
| 152 | acond,&! coefficient used to compute ztcondsol |
|---|
| 153 | bcond,&! coefficient used to compute ztcondsol |
|---|
| 154 | ccond ! coefficient used to compute ztcondsol |
|---|
| 155 | |
|---|
| 156 | logical, save :: & |
|---|
| 157 | firstcall = .true. ! Used to compute saved variables |
|---|
| 158 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 159 | !----3) Variables: |
|---|
| 160 | !----------------- |
|---|
| 161 | integer :: & |
|---|
| 162 | l, &! loop on layers |
|---|
| 163 | ig, &! loop on ngrid points |
|---|
| 164 | iq ! loop on tracer |
|---|
| 165 | |
|---|
| 166 | real :: & |
|---|
| 167 | qco2, &! effective quantity of CO2, used to compute mean molecular mass |
|---|
| 168 | mmean, &! mean molecular mass |
|---|
| 169 | zfallheat, &! aerodynamical friction and energy used to heat the amount of ice |
|---|
| 170 | zmflux(nlayer+1), &! mass fluxes through the sigma levels (kg.m-2.s-1) |
|---|
| 171 | zu(nlayer), &! effective zonal wind |
|---|
| 172 | zv(nlayer), &! effective meridional wind |
|---|
| 173 | zq(nlayer,nq), &! effective tracers quantities |
|---|
| 174 | zq1(nlayer), &! buffer of zq |
|---|
| 175 | ztsrf(ngrid), &! effective temperature at the surface |
|---|
| 176 | ztm(nlayer+1), &! temperature fluxes through the sigma levels |
|---|
| 177 | zum(nlayer+1), &! zonal wind fluxes through the sigma levels |
|---|
| 178 | zvm(nlayer+1), &! meridional wind fluxes through the sigma levels |
|---|
| 179 | zqm(nlayer+1,nq), &! quantity of tracers flux through the sigma levels |
|---|
| 180 | zqm1(nlayer), &! quantity of tracers after Van-Leer scheme |
|---|
| 181 | masse(nlayer), &! mass layer (kg) |
|---|
| 182 | w(nlayer+1), &! total mass fluxes through the sigma levels during ptimestep (kg.m-2) |
|---|
| 183 | availco2, &! available quantity of co2 (kg) |
|---|
| 184 | emisref(ngrid), &! emissivity reference |
|---|
| 185 | vmr_co2(ngrid,nlayer), &! CO2 volume mixing ratio |
|---|
| 186 | zt(nlayer), &! effective temperature in the atmosphere (K) |
|---|
| 187 | ztcond(ngrid,nlayer+1), &! CO2 condensation temperature (atm) |
|---|
| 188 | ztcondsol(ngrid), &! CO2 condensation temperature (surface) |
|---|
| 189 | zdiceco2(ngrid), &! tendency on co2ice surface tracer (kg/m2/s) |
|---|
| 190 | zcondicea(ngrid,nlayer),&! condensation rate in layer l (kg/m2/s) |
|---|
| 191 | zcondices(ngrid), &! condensation rate on the ground (kg/m2/s) |
|---|
| 192 | zfallice(ngrid) ! amount of ice falling from layer l (kg/m2/s) |
|---|
| 193 | |
|---|
| 194 | logical :: & |
|---|
| 195 | condsub(ngrid) ! True if there is condensation/sublimation (used for co2snow) |
|---|
| 196 | !======================================================================================================================! |
|---|
| 197 | ! BEGIN ===============================================================================================================! |
|---|
| 198 | !======================================================================================================================! |
|---|
| 199 | ! 1. Initialization |
|---|
| 200 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 201 | availco2 = 0. |
|---|
| 202 | zfallheat = 0. |
|---|
| 203 | zt(1:nlayer) = 0. |
|---|
| 204 | ztcond(1:ngrid, 1:nlayer+1) = 0. |
|---|
| 205 | ztcondsol(1:ngrid) = 0. |
|---|
| 206 | zmflux(1:nlayer+1) = 0. |
|---|
| 207 | zu(1:nlayer) = 0. |
|---|
| 208 | zv(1:nlayer) = 0. |
|---|
| 209 | zq(1:nlayer, 1:nq) = 0. |
|---|
| 210 | zq1(1:nlayer) = 0. |
|---|
| 211 | ztsrf(1:ngrid) = 0. |
|---|
| 212 | ztm(1:nlayer+1) = 0. |
|---|
| 213 | zum(1:nlayer+1) = 0. |
|---|
| 214 | zvm(1:nlayer+1) = 0. |
|---|
| 215 | zqm(1:nlayer+1, 1:nq) = 0. |
|---|
| 216 | masse(1:nlayer) = 0. |
|---|
| 217 | w(1:nlayer+1) = 0. |
|---|
| 218 | emisref(1:ngrid) = 0. |
|---|
| 219 | vmr_co2(1:ngrid, 1:nlayer) = 0. |
|---|
| 220 | zcondices(1:ngrid) = 0. |
|---|
| 221 | pdtsrfc(1:ngrid) = 0. |
|---|
| 222 | pdpsrf(1:ngrid) = 0. |
|---|
| 223 | zdiceco2(1:ngrid) = 0. |
|---|
| 224 | condsub(1:ngrid) = .false. |
|---|
| 225 | zcondicea(1:ngrid, 1:nlayer) = 0. |
|---|
| 226 | zfallice(1:ngrid) = 0. |
|---|
| 227 | pduc(1:ngrid, 1:nlayer) = 0. |
|---|
| 228 | pdvc(1:ngrid, 1:nlayer) = 0. |
|---|
| 229 | pdqc(1:ngrid, 1:nlayer, 1:nq) = 0. |
|---|
| 230 | pdtc(1:ngrid,1:nlayer) = 0. |
|---|
| 231 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 232 | ! 2. Firstcall |
|---|
| 233 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 234 | ! AS: firstcall OK absolute |
|---|
| 235 | if (firstcall) then |
|---|
| 236 | firstcall = .false. |
|---|
| 237 | |
|---|
| 238 | bcond = 1. / tcond1mb |
|---|
| 239 | ccond = cpp / (g*latcond) |
|---|
| 240 | acond = r / latcond |
|---|
| 241 | |
|---|
| 242 | write(*,*)'CO2condens: improved_ztcond=', improved_ztcond |
|---|
| 243 | write(*,*)'In co2condens:Tcond(P=1mb)=', tcond1mb, ' Lcond=', latcond |
|---|
| 244 | write(*,*)'acond,bcond,ccond', acond, bcond, ccond |
|---|
| 245 | |
|---|
| 246 | ! Prepare Special treatment if one of the tracer is CO2 gas. Compute A and B coefficient use to compute mean molecular |
|---|
| 247 | ! mass Mair defined by: |
|---|
| 248 | ! 1/Mair = q(igcm_co2)/m_co2 + (1-q(igcm_co2))/m_noco2 |
|---|
| 249 | ! 1/Mair = A*q(igcm_co2) + B |
|---|
| 250 | A = (1./m_co2 - 1./m_noco2) |
|---|
| 251 | B = 1./m_noco2 |
|---|
| 252 | end if ! of IF (firstcall) |
|---|
| 253 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 254 | ! 3. Compute CO2 Volume mixing ratio |
|---|
| 255 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 256 | if (improved_ztcond.and.(igcm_co2/=0)) then |
|---|
| 257 | do l = 1, nlayer |
|---|
| 258 | do ig = 1, ngrid |
|---|
| 259 | qco2 = pq(ig,l,igcm_co2) + pdq(ig,l,igcm_co2)*ptimestep |
|---|
| 260 | ! Mean air molecular mass = 1/(q(igcm_co2)/m_co2 + (1-q(igcm_co2))/m_noco2) |
|---|
| 261 | mmean = 1. / (A*qco2 +B) |
|---|
| 262 | vmr_co2(ig,l) = (qco2*mmean) / m_co2 |
|---|
| 263 | end do |
|---|
| 264 | end do |
|---|
| 265 | else |
|---|
| 266 | do l = 1, nlayer |
|---|
| 267 | do ig = 1, ngrid |
|---|
| 268 | vmr_co2(ig,l) = 0.95 |
|---|
| 269 | end do |
|---|
| 270 | end do |
|---|
| 271 | end if |
|---|
| 272 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 273 | ! 4. Set zcondicea, zfallice from co2clouds condensation rate |
|---|
| 274 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 275 | do l = nlayer, 1, -1 |
|---|
| 276 | do ig = 1, ngrid |
|---|
| 277 | zcondicea(ig,l) = pcondicea_co2microp(ig,l) * (pplev(ig,l) - pplev(ig,l+1))/g |
|---|
| 278 | end do |
|---|
| 279 | end do |
|---|
| 280 | |
|---|
| 281 | ! Only sedimentation falls on the ground ! |
|---|
| 282 | do ig = 1, ngrid |
|---|
| 283 | zfallice(ig) = zdqssed_co2(ig) |
|---|
| 284 | piceco2(ig) = piceco2(ig) + zfallice(ig)*ptimestep |
|---|
| 285 | end do |
|---|
| 286 | |
|---|
| 287 | call writediagfi(ngrid, "zcondicea","Condensation rate in layers", " ", 3, zcondicea) |
|---|
| 288 | |
|---|
| 289 | call writediagfi(ngrid,"zfallice", "Sedimentation rate", " ", 2, zfallice) |
|---|
| 290 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 291 | ! 5. Main co2condens |
|---|
| 292 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 293 | do ig = 1, ngrid |
|---|
| 294 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 295 | ! 5.1. Forecast of ground temperature ztsrf and frost temperature ztcondsol |
|---|
| 296 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 297 | ztcondsol(ig) = 1. / (bcond-acond*log(.01*vmr_co2(ig,1)*pplev(ig,1))) |
|---|
| 298 | ztsrf(ig) = ptsrf(ig) + pdtsrf(ig)*ptimestep |
|---|
| 299 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 300 | ! 5.2. Check if we have condensation/sublimation on the ground |
|---|
| 301 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 302 | ! ground condensation || falling snow || ground sublimation |
|---|
| 303 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 304 | if ((ztsrf(ig)<ztcondsol(ig)) .OR. (zfallice(ig)/=0.) .OR. ((ztsrf(ig)>ztcondsol(ig)) .AND. (piceco2(ig)/=0.))) then |
|---|
| 305 | condsub(ig) = .true. |
|---|
| 306 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 307 | ! 5.3. Compute zfallheat |
|---|
| 308 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 309 | if (zfallice(ig)>0) then |
|---|
| 310 | zfallheat = zfallice(ig) * ( pphi(ig,1) - phisfi(ig) + cpice*(ztcond(ig,1) - ztcondsol(ig)) )/latcond |
|---|
| 311 | else |
|---|
| 312 | zfallheat = 0. |
|---|
| 313 | end if |
|---|
| 314 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 315 | ! 5.4. Compute direct condensation/sublimation of CO2 ice |
|---|
| 316 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 317 | zcondices(ig) = pcapcal(ig) * (ztcondsol(ig)-ztsrf(ig)) / (latcond*ptimestep) - zfallheat |
|---|
| 318 | pdtsrfc(ig) = (ztcondsol(ig) - ztsrf(ig)) / ptimestep |
|---|
| 319 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 320 | ! 5.4.a. If there is not enough CO2 tracer in 1st layer to condense |
|---|
| 321 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 322 | if (igcm_co2/=0) then |
|---|
| 323 | ! Available CO2 tracer in layer 1 at end of timestep (kg/m2) |
|---|
| 324 | availco2 = pq(ig,1,igcm_co2) * ( (ap(1)-ap(2)) + (bp(1)-bp(2)) * (pplev(ig,1)/g - zcondices(ig)*ptimestep) ) |
|---|
| 325 | if ( zcondices(ig) * ptimestep>availco2 ) then |
|---|
| 326 | zcondices(ig) = availco2/ptimestep |
|---|
| 327 | zdiceco2(ig) = zcondices(ig) + zfallice(ig) |
|---|
| 328 | pdtsrfc(ig) = (latcond/pcapcal(ig)) * (zcondices(ig)+zfallheat) |
|---|
| 329 | end if |
|---|
| 330 | end if |
|---|
| 331 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 332 | ! 5.4.b. If the entire CO2 ice layer sublimes (including what has just condensed in the atmosphere) |
|---|
| 333 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 334 | if ( (piceco2(ig)/ptimestep) <= -zcondices(ig) ) then |
|---|
| 335 | zcondices(ig) = -piceco2(ig)/ptimestep |
|---|
| 336 | pdtsrfc(ig) = (latcond/pcapcal(ig)) * (zcondices(ig)+zfallheat) |
|---|
| 337 | end if |
|---|
| 338 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 339 | ! 5.5. Changing CO2 ice amount and pressure |
|---|
| 340 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 341 | zdiceco2(ig) = zcondices(ig) + zfallice(ig) + sum(zcondicea(ig,:)) |
|---|
| 342 | piceco2(ig) = piceco2(ig) + zcondices(ig)*ptimestep |
|---|
| 343 | pdpsrf(ig) = -zdiceco2(ig) * g |
|---|
| 344 | |
|---|
| 345 | if (abs(pdpsrf(ig)*ptimestep)>pplev(ig,1)) then |
|---|
| 346 | print *, 'STOP in condens' |
|---|
| 347 | print *, 'condensing more than total mass' |
|---|
| 348 | print *, 'Grid point ', ig |
|---|
| 349 | print *, 'Longitude(degrees): ', longitude_deg(ig) |
|---|
| 350 | print *, 'Latitude(degrees): ', latitude_deg(ig) |
|---|
| 351 | print *, 'Ps = ', pplev(ig,1) |
|---|
| 352 | print *, 'd Ps = ', pdpsrf(ig) |
|---|
| 353 | call abort_physic('co2condens4micro', 'condensing more than total mass', 1) |
|---|
| 354 | end if |
|---|
| 355 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 356 | ! 5.6. Surface albedo and emissivity of the surface below the snow (emisref) |
|---|
| 357 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 358 | ! 5.6.a. Check that amont of CO2 ice is not problematic |
|---|
| 359 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 360 | if(.not.piceco2(ig)>=0.) then |
|---|
| 361 | if(piceco2(ig)<=-5.e-8) then |
|---|
| 362 | write(*,*)'WARNING co2condens piceco2(', ig, ')=', piceco2(ig) |
|---|
| 363 | piceco2(ig) = 0. |
|---|
| 364 | end if |
|---|
| 365 | end if |
|---|
| 366 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 367 | ! 5.6.b. Set albedo and emissivity of the surface |
|---|
| 368 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 369 | call albedocaps(zls,ngrid,piceco2,psolaralb,emisref) |
|---|
| 370 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 371 | ! 5.6.c. Set pemisurf to emissiv when there is bare surface (needed for co2snow) |
|---|
| 372 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 373 | if (piceco2(ig)==0) then |
|---|
| 374 | pemisurf(ig) = emissiv |
|---|
| 375 | end if |
|---|
| 376 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 377 | ! 5.7. Correction to account for redistribution between sigma or hybrid layers when changing surface pressure (and |
|---|
| 378 | ! warming/cooling of the CO2 currently changing phase). |
|---|
| 379 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 380 | do l= 1, nlayer |
|---|
| 381 | zt(l) = pt(ig,l) + pdt(ig,l)*ptimestep |
|---|
| 382 | zu(l) = pu(ig,l) + pdu(ig,l)*ptimestep |
|---|
| 383 | zv(l) = pv(ig,l) + pdv(ig,l)*ptimestep |
|---|
| 384 | do iq=1,nq |
|---|
| 385 | zq(l,iq) = pq(ig,l,iq) + pdq(ig,l,iq)*ptimestep |
|---|
| 386 | end do |
|---|
| 387 | end do |
|---|
| 388 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 389 | ! 5.7.a. Mass fluxes through the sigma levels (kg.m-2.s-1) (>0 when up) |
|---|
| 390 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 391 | zmflux(1) = - zcondices(ig) - zfallice(ig) |
|---|
| 392 | do l = 1, nlayer |
|---|
| 393 | zmflux(l+1) = zmflux(l) - zcondicea(ig,l) & |
|---|
| 394 | #ifndef MESOSCALE |
|---|
| 395 | + (bp(l)-bp(l+1)) * (-pdpsrf(ig)/g) |
|---|
| 396 | ! zmflux set to 0 if very low to avoid: top layer is disappearing in v1ld |
|---|
| 397 | if (abs(zmflux(l+1))<1E-13.OR.bp(l+1)==0.) then |
|---|
| 398 | zmflux(l+1) = 0. |
|---|
| 399 | end if |
|---|
| 400 | #else |
|---|
| 401 | if (abs(zmflux(l+1))<1E-13) then |
|---|
| 402 | zmflux(l+1) = 0. |
|---|
| 403 | end if |
|---|
| 404 | #endif |
|---|
| 405 | end do |
|---|
| 406 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 407 | ! 5.7.b. Mass of each layer at the end of timestep |
|---|
| 408 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 409 | do l = 1, nlayer |
|---|
| 410 | masse(l) = (pplev(ig,l) - pplev(ig,l+1) + (bp(l)-bp(l+1))*pdpsrf(ig)*ptimestep)/g |
|---|
| 411 | end do |
|---|
| 412 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 413 | ! 5.7.c. Corresponding fluxes for T, U, V and Q (averaging operator for transport) |
|---|
| 414 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 415 | ! 5.7.c.i. Value transfert at the surface interface when condensation/sublimation |
|---|
| 416 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 417 | ztm(1) = ztsrf(ig) + pdtsrfc(ig)*ptimestep |
|---|
| 418 | zum(1) = 0 |
|---|
| 419 | zvm(1) = 0 |
|---|
| 420 | |
|---|
| 421 | ! Most tracer do not condense |
|---|
| 422 | do iq = 1, nq |
|---|
| 423 | zqm(1,iq) = 0. |
|---|
| 424 | end do |
|---|
| 425 | |
|---|
| 426 | ! Special case if one of the tracer is CO2 gas |
|---|
| 427 | if (igcm_co2/=0) then |
|---|
| 428 | zqm(1,igcm_co2) = 1. ! flux is 100% CO2 |
|---|
| 429 | end if |
|---|
| 430 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 431 | ! 5.7.c.ii. Van Leer scheme |
|---|
| 432 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 433 | do l=1,nlayer |
|---|
| 434 | w(l)=-zmflux(l)*ptimestep |
|---|
| 435 | end do |
|---|
| 436 | |
|---|
| 437 | call vl1d(nlayer,zt,2.,masse,w,ztm) |
|---|
| 438 | call vl1d(nlayer,zu ,2.,masse,w,zum) |
|---|
| 439 | call vl1d(nlayer,zv ,2.,masse,w,zvm) |
|---|
| 440 | |
|---|
| 441 | do iq=1, nq |
|---|
| 442 | do l=1, nlayer |
|---|
| 443 | zq1(l) = zq(l,iq) |
|---|
| 444 | end do |
|---|
| 445 | |
|---|
| 446 | zqm1(1) = zqm(1,iq) |
|---|
| 447 | zqm1(2:nlayer) = 0. |
|---|
| 448 | |
|---|
| 449 | call vl1d(nlayer,zq1,2.,masse,w,zqm1) |
|---|
| 450 | |
|---|
| 451 | do l = 2, nlayer |
|---|
| 452 | zq(l,iq) = zq1(l) |
|---|
| 453 | zqm(l,iq) = zqm1(l) |
|---|
| 454 | end do |
|---|
| 455 | end do |
|---|
| 456 | |
|---|
| 457 | ! Surface condensation affects low winds |
|---|
| 458 | if (zmflux(1)<0) then |
|---|
| 459 | zum(1) = zu(1) * (w(1)/masse(1)) |
|---|
| 460 | zvm(1) = zv(1) * (w(1)/masse(1)) |
|---|
| 461 | if (w(1)>masse(1)) then ! ensure numerical stability |
|---|
| 462 | zum(1) = ((zu(1)-zum(2))*masse(1)/w(1)) + zum(2) |
|---|
| 463 | zvm(1) = ((zv(1)-zvm(2))*masse(1)/w(1)) + zvm(2) |
|---|
| 464 | end if |
|---|
| 465 | end if |
|---|
| 466 | |
|---|
| 467 | ztm(nlayer+1) = zt(nlayer) ! should not be used, but... |
|---|
| 468 | zum(nlayer+1) = zu(nlayer) ! should not be used, but... |
|---|
| 469 | zvm(nlayer+1) = zv(nlayer) ! should not be used, but... |
|---|
| 470 | |
|---|
| 471 | do iq = 1, nq |
|---|
| 472 | zqm(nlayer+1,iq) = zq(nlayer,iq) |
|---|
| 473 | end do |
|---|
| 474 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 475 | ! 5.7.c.iii Compute tendencies on T, U, V, Q |
|---|
| 476 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 477 | #ifdef MESOSCALE |
|---|
| 478 | ! AS: This part must be commented in the mesoscale model |
|---|
| 479 | ! AS: ... to avoid instabilities. |
|---|
| 480 | ! AS: you have to compile with -DMESOSCALE to do so |
|---|
| 481 | #else |
|---|
| 482 | do l = 1, nlayer |
|---|
| 483 | if (.not. co2clouds) then |
|---|
| 484 | ! Tendencies on T |
|---|
| 485 | pdtc(ig,l) = (1./masse(l)) * ( zmflux(l)*(ztm(l) - zt(l)) - zmflux(l+1)*(ztm(l+1) - zt(l)) & |
|---|
| 486 | + zcondicea(ig,l)*(ztcond(ig,l)-zt(l)) ) |
|---|
| 487 | else |
|---|
| 488 | pdtc(ig,l) = (1./masse(l)) * ( zmflux(l)*(ztm(l) - zt(l)) - zmflux(l+1)*(ztm(l+1) - zt(l)) ) |
|---|
| 489 | end if |
|---|
| 490 | |
|---|
| 491 | ! Tendencies on U |
|---|
| 492 | pduc(ig,l) = (1./masse(l)) * ( zmflux(l)*(zum(l) - zu(l)) - zmflux(l+1)*(zum(l+1) - zu(l)) ) |
|---|
| 493 | |
|---|
| 494 | ! Tendencies on V |
|---|
| 495 | pdvc(ig,l) = (1./masse(l)) * ( zmflux(l)*(zvm(l) - zv(l)) - zmflux(l+1)*(zvm(l+1) - zv(l)) ) |
|---|
| 496 | end do |
|---|
| 497 | #endif |
|---|
| 498 | ! Tendencies on Q |
|---|
| 499 | do iq = 1, nq |
|---|
| 500 | if (iq==igcm_co2) then |
|---|
| 501 | do l = 1, nlayer |
|---|
| 502 | pdqc(ig,l,iq) = (1./masse(l)) * (zmflux(l)*(zqm(l,iq) - zq(l,iq))- zmflux(l+1)*(zqm(l+1,iq) - zq(l,iq))& |
|---|
| 503 | + zcondicea(ig,l)*(zq(l,iq) - 1.)) |
|---|
| 504 | end do |
|---|
| 505 | else |
|---|
| 506 | do l = 1, nlayer |
|---|
| 507 | pdqc(ig,l,iq) = (1./masse(l)) * ( zmflux(l)*(zqm(l,iq)-zq(l,iq)) - zmflux(l+1)*(zqm(l+1,iq)-zq(l,iq))& |
|---|
| 508 | + zcondicea(ig,l)*zq(l,iq) ) |
|---|
| 509 | end do |
|---|
| 510 | end if |
|---|
| 511 | end do |
|---|
| 512 | end if ! if |
|---|
| 513 | end do ! loop on ig |
|---|
| 514 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 515 | ! 6. CO2 snow / clouds scheme |
|---|
| 516 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 517 | call co2snow(ngrid, nlayer, ptimestep, emisref, condsub, pplev, zcondicea, zcondices, zdqssed_co2*ptimestep, & |
|---|
| 518 | pemisurf) |
|---|
| 519 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 520 | ! 7. Extra special case for surface temperature tendency pdtsrfc: |
|---|
| 521 | ! We want to fix the south pole temperature to CO2 condensation temperature. |
|---|
| 522 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 523 | #ifndef MESOSCALE |
|---|
| 524 | if (caps.and.(obliquit<27.)) then |
|---|
| 525 | ! check if last grid point is the south pole |
|---|
| 526 | if (abs(latitude(ngrid)-(-pi/2.))<1.e-5) then |
|---|
| 527 | ! NB: Updated surface pressure, at grid point 'ngrid', is ps(ngrid)=pplev(ngrid,1)+pdpsrf(ngrid)*ptimestep |
|---|
| 528 | ! write(*,*)"co2condens: South pole: latitude(ngrid)=", latitude(ngrid) |
|---|
| 529 | ztcondsol(ngrid) = 1./(bcond-acond*log(.01*vmr_co2(ngrid,1) * (pplev(ngrid,1)+pdpsrf(ngrid)*ptimestep))) |
|---|
| 530 | pdtsrfc(ngrid) = (ztcondsol(ngrid)-ztsrf(ngrid))/ptimestep |
|---|
| 531 | end if |
|---|
| 532 | end if |
|---|
| 533 | #endif |
|---|
| 534 | !======================================================================================================================! |
|---|
| 535 | ! END =================================================================================================================! |
|---|
| 536 | !======================================================================================================================! |
|---|
| 537 | end subroutine co2condens4micro |
|---|
| 538 | |
|---|
| 539 | |
|---|
| 540 | !**********************************************************************************************************************! |
|---|
| 541 | !**********************************************************************************************************************! |
|---|
| 542 | |
|---|
| 543 | |
|---|
| 544 | !======================================================================================================================! |
|---|
| 545 | ! SUBROUTINE: Van-Leer scheme =========================================================================================! |
|---|
| 546 | !======================================================================================================================! |
|---|
| 547 | ! Subject: |
|---|
| 548 | !--------- |
|---|
| 549 | ! Operateur de moyenne inter-couche pour calcul de transport type Van-Leer " pseudo amont " dans la verticale |
|---|
| 550 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 551 | ! Comments: |
|---|
| 552 | !---------- |
|---|
| 553 | ! q,w are input arguments for the s-pg .... |
|---|
| 554 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 555 | ! Paper: |
|---|
| 556 | !------- |
|---|
| 557 | ! Van-Leer (1977), "Towards the Ultimate Conservative Difference Scheme. IV. A New Approach to Numerical Convection" |
|---|
| 558 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 559 | subroutine vl1d(nlayer,q,pente_max,masse,w,qm) |
|---|
| 560 | |
|---|
| 561 | implicit none |
|---|
| 562 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 563 | ! VARIABLE DECLARATION |
|---|
| 564 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 565 | ! Input arguments: |
|---|
| 566 | !----------------- |
|---|
| 567 | integer, intent(in) :: & |
|---|
| 568 | nlayer ! number of layers |
|---|
| 569 | |
|---|
| 570 | real, intent(in) :: & |
|---|
| 571 | pente_max, &! coefficient, pente_max = 2 advised |
|---|
| 572 | masse(nlayer), &! masse layer Dp/g (kg) |
|---|
| 573 | q(nlayer) ! quantity of tracer |
|---|
| 574 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 575 | ! In-Output arguments: |
|---|
| 576 | !--------------------- |
|---|
| 577 | real, intent(inout) :: & |
|---|
| 578 | w(nlayer+1) ! masse d'atm ``transferee'' a chaque pas de temps (kg.m-2) |
|---|
| 579 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 580 | ! Output arguments: |
|---|
| 581 | !------------------ |
|---|
| 582 | real, intent(out) :: & |
|---|
| 583 | qm(nlayer+1) ! quantity of tracer after Van-Leer scheme |
|---|
| 584 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 585 | ! Locals variables: |
|---|
| 586 | !------------------ |
|---|
| 587 | integer :: & |
|---|
| 588 | l, &! loop on nlayer |
|---|
| 589 | m ! index |
|---|
| 590 | |
|---|
| 591 | real :: & |
|---|
| 592 | dzqmax, &! maximum of dzq between two adjacent layers |
|---|
| 593 | sigw, &! |
|---|
| 594 | Mtot, &! |
|---|
| 595 | MQtot, &! |
|---|
| 596 | dzq(nlayer), &! |
|---|
| 597 | dzqw(nlayer),&! |
|---|
| 598 | adzqw(nlayer) ! |
|---|
| 599 | !======================================================================================================================! |
|---|
| 600 | ! BEGIN ===============================================================================================================! |
|---|
| 601 | !======================================================================================================================! |
|---|
| 602 | ! 1. On oriente tout dans le sens de la pression: w > 0 when down |
|---|
| 603 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 604 | do l = 2, nlayer |
|---|
| 605 | dzqw(l) = q(l-1) - q(l) |
|---|
| 606 | adzqw(l) = abs(dzqw(l)) |
|---|
| 607 | end do |
|---|
| 608 | |
|---|
| 609 | do l = 2, nlayer-1 |
|---|
| 610 | if(dzqw(l)*dzqw(l+1)>0.) then |
|---|
| 611 | dzq(l) = 0.5 * (dzqw(l)+dzqw(l+1)) |
|---|
| 612 | else |
|---|
| 613 | dzq(l) = 0. |
|---|
| 614 | end if |
|---|
| 615 | |
|---|
| 616 | dzqmax = pente_max * min(adzqw(l), adzqw(l+1)) |
|---|
| 617 | |
|---|
| 618 | dzq(l) = sign(min(abs(dzq(l)),dzqmax), dzq(l)) |
|---|
| 619 | end do |
|---|
| 620 | |
|---|
| 621 | dzq(1)=0. |
|---|
| 622 | dzq(nlayer)=0. |
|---|
| 623 | |
|---|
| 624 | do l = 1, nlayer-1 |
|---|
| 625 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 626 | ! 2.1. Regular scheme (transfered mass < layer mass) |
|---|
| 627 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 628 | if (w(l+1)>0. .and. w(l+1)<=masse(l+1)) then |
|---|
| 629 | sigw = w(l+1) / masse(l+1) |
|---|
| 630 | qm(l+1) = (q(l+1) + 0.5*(1.-sigw)*dzq(l+1)) |
|---|
| 631 | else if (w(l+1)<=0. .and. -w(l+1)<=masse(l)) then |
|---|
| 632 | sigw = w(l+1) / masse(l) |
|---|
| 633 | qm(l+1) = (q(l) - 0.5*(1.+sigw)*dzq(l)) |
|---|
| 634 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 635 | ! 2.2. Extended scheme (transfered mass > layer mass) |
|---|
| 636 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 637 | else if (w(l+1)>0.) then |
|---|
| 638 | m = l+1 |
|---|
| 639 | Mtot = masse(m) |
|---|
| 640 | MQtot = masse(m)*q(m) |
|---|
| 641 | |
|---|
| 642 | do while ((m<nlayer).and.(w(l+1)>(Mtot+masse(m+1)))) |
|---|
| 643 | m = m+1 |
|---|
| 644 | Mtot = Mtot + masse(m) |
|---|
| 645 | MQtot = MQtot + masse(m)*q(m) |
|---|
| 646 | end do |
|---|
| 647 | |
|---|
| 648 | if (m<nlayer) then |
|---|
| 649 | sigw = (w(l+1)-Mtot) / masse(m+1) |
|---|
| 650 | qm(l+1) = (1/w(l+1))*( MQtot + (w(l+1)-Mtot)* (q(m+1)+0.5*(1.-sigw)*dzq(m+1)) ) |
|---|
| 651 | else |
|---|
| 652 | w(l+1) = Mtot |
|---|
| 653 | qm(l+1) = Mqtot / Mtot |
|---|
| 654 | call abort_physic('co2condens4micro', 'top layer is disapearing !', 1) |
|---|
| 655 | end if |
|---|
| 656 | !----------------------------------------------------------------------------------------------------------------------! |
|---|
| 657 | else ! if(w(l+1).lt.0) |
|---|
| 658 | m = l-1 |
|---|
| 659 | Mtot = masse(m+1) |
|---|
| 660 | MQtot = masse(m+1)*q(m+1) |
|---|
| 661 | if (m>0) then ! because some compilers will have problems evaluating masse(0) |
|---|
| 662 | do while ((m>0).and.(-w(l+1)>(Mtot+masse(m)))) |
|---|
| 663 | m = m-1 |
|---|
| 664 | Mtot = Mtot + masse(m+1) |
|---|
| 665 | MQtot = MQtot + masse(m+1)*q(m+1) |
|---|
| 666 | if (m==0) then |
|---|
| 667 | call abort_physic('co2condens_mod4micro','vl1d',1) |
|---|
| 668 | end if |
|---|
| 669 | end do |
|---|
| 670 | end if |
|---|
| 671 | |
|---|
| 672 | if (m>0) then |
|---|
| 673 | sigw = (w(l+1)+Mtot) / masse(m) |
|---|
| 674 | qm(l+1) = (-1/w(l+1)) * ( MQtot + (-w(l+1)-Mtot) * (q(m)-0.5*(1.+sigw)*dzq(m)) ) |
|---|
| 675 | else |
|---|
| 676 | qm(l+1) = (-1/w(l+1)) * (MQtot + (-w(l+1)-Mtot)*qm(1)) |
|---|
| 677 | end if |
|---|
| 678 | end if |
|---|
| 679 | end do ! l = 1, nlayer-1 |
|---|
| 680 | !======================================================================================================================! |
|---|
| 681 | ! END =================================================================================================================! |
|---|
| 682 | !======================================================================================================================! |
|---|
| 683 | end subroutine vl1d |
|---|
| 684 | |
|---|
| 685 | end module co2condens_mod4micro |
|---|
