- Timestamp:
- Jul 28, 2025, 6:44:28 PM (4 months ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
LMDZ6/branches/contrails/libf/phylmd/lmdz_lscp_condensation.f90
r5779 r5790 97 97 cldfra_in, qvc_in, qliq_in, qice_in, shear, pbl_eps, cell_area, & 98 98 temp, qtot_in, qsat, gamma_cond, ratqs, keepgoing, pt_pron_clds, & 99 dzsed_abv, flsed_abv, cfsed_abv, & 100 dzsed_lincont_abv, flsed_lincont_abv, cfsed_lincont_abv, & 101 dzsed_circont_abv, flsed_circont_abv, cfsed_circont_abv, & 102 dzsed, flsed, cfsed, dzsed_lincont, flsed_lincont, cfsed_lincont, & 103 dzsed_circont, flsed_circont, cfsed_circont, flauto, & 104 cldfra, qincld, qvc, issrfra, qissr, dcf_sub, dcf_con, dcf_mix, dcf_sed, dcf_auto, & 99 dzsed_abv, flsed_abv, cfsed_abv, dzsed, flsed, cfsed, flauto, & 100 cldfra, qincld, qvc, issrfra, qissr, & 101 dcf_sub, dcf_con, dcf_mix, dcf_sed, dcf_auto, & 105 102 dqi_adj, dqi_sub, dqi_con, dqi_mix, dqi_sed, dqi_auto, & 106 103 dqvc_adj, dqvc_sub, dqvc_con, dqvc_mix, dqvc_sed, dqvc_auto, & 107 lincontfra_in, circontfra_in, qtl_in, qtc_in, flight_dist, flight_h2o, &108 lincontfra, circontfra, qlincont, qcircont, &109 Tcritcont, qcritcont, potcontfraP, potcontfraNP, &110 d cfl_ini, dqil_ini, dqtl_ini, dcfl_sub, dqil_sub, dqtl_sub, &111 d cfl_cir, dqtl_cir, dcfl_mix, dqil_mix, dqtl_mix, &112 dcf l_sed, dqil_sed, dqtl_sed, dcfl_auto, dqil_auto, dqtl_auto, &113 dcfc_ sub, dqic_sub, dqtc_sub, dcfc_mix, dqic_mix, dqtc_mix, &114 dcfc_sed, dqic_sed, dqtc_sed, d cfc_auto, dqic_auto, dqtc_auto)104 contfra_in, qtc_in, nic_in, flight_dist, flight_fuel, & 105 contfra, qcont, Ncont, Tcritcont, qcritcont, potcontfraP, potcontfraNP, & 106 AEI_contrails, AEI_surv_contrails, fsurv_contrails, section_contrails, & 107 dzsed_cont_abv, flsed_cont_abv, Nflsed_cont_abv, cfsed_cont_abv, & 108 dzsed_cont, flsed_cont, Nflsed_cont, cfsed_cont, & 109 dcfc_ini, dqic_ini, dqtc_ini, dnic_ini, dcfc_sub, dqic_sub, dqtc_sub, dnic_sub, & 110 dcfc_mix, dqic_mix, dqtc_mix, dnic_mix, dnic_agg, & 111 dcfc_sed, dqic_sed, dqtc_sed, dnic_sed, dcfc_auto, dqic_auto, dqtc_auto, dnic_auto) 115 112 116 113 !---------------------------------------------------------------------- … … 140 137 USE lmdz_lscp_ini, ONLY: chi_sedim 141 138 142 USE lmdz_lscp_ini, ONLY: ok_plane_contrail, aspect_ratio_ lincontrails143 USE lmdz_lscp_ini, ONLY: coef_mixing_ lincontrails, coef_shear_lincontrails144 USE lmdz_lscp_ini, ONLY: chi_mixing_ lincontrails, linear_contrails_lifetime145 USE lmdz_lscp_ ini, ONLY: fallice_linear_contrails, fallice_cirrus_contrails139 USE lmdz_lscp_ini, ONLY: ok_plane_contrail, aspect_ratio_contrails 140 USE lmdz_lscp_ini, ONLY: coef_mixing_contrails, coef_shear_contrails 141 USE lmdz_lscp_ini, ONLY: chi_mixing_contrails, eff2vol_radius_contrails 142 USE lmdz_lscp_tools, ONLY: ICECRYSTAL_VELO 146 143 USE lmdz_aviation, ONLY: contrails_formation 147 144 … … 175 172 REAL, INTENT(IN) , DIMENSION(klon) :: flsed_abv ! sedimentated ice flux FROM THE LAYER ABOVE [kg/s/m2] 176 173 REAL, INTENT(IN) , DIMENSION(klon) :: cfsed_abv ! cloud fraction IN THE LAYER ABOVE [-] 177 REAL, INTENT(IN) , DIMENSION(klon) :: dzsed_lincont_abv ! sedimentated linear contrails height IN THE LAYER ABOVE [m]178 REAL, INTENT(IN) , DIMENSION(klon) :: flsed_lincont_abv ! sedimentated ice flux in linear contrails FROM THE LAYER ABOVE [kg/s/m2]179 REAL, INTENT(IN) , DIMENSION(klon) :: cfsed_lincont_abv ! linear contrails fraction IN THE LAYER ABOVE [-]180 REAL, INTENT(IN) , DIMENSION(klon) :: dzsed_circont_abv ! sedimentated contrails cirrus height IN THE LAYER ABOVE [m]181 REAL, INTENT(IN) , DIMENSION(klon) :: flsed_circont_abv ! sedimentated ice flux in contrails cirrus FROM THE LAYER ABOVE [kg/s/m2]182 REAL, INTENT(IN) , DIMENSION(klon) :: cfsed_circont_abv ! contrails cirrus fraction IN THE LAYER ABOVE [-]183 174 REAL, INTENT(INOUT), DIMENSION(klon) :: dzsed ! sedimentated cloud height [m] 184 175 REAL, INTENT(INOUT), DIMENSION(klon) :: flsed ! sedimentated ice flux [kg/s/m2] 185 176 REAL, INTENT(INOUT), DIMENSION(klon) :: cfsed ! sedimentated cloud fraction [-] 186 REAL, INTENT(INOUT), DIMENSION(klon) :: dzsed_lincont ! sedimentated linear contrails height [m]187 REAL, INTENT(INOUT), DIMENSION(klon) :: flsed_lincont ! sedimentated ice flux in linear contrails [kg/s/m2]188 REAL, INTENT(INOUT), DIMENSION(klon) :: cfsed_lincont ! sedimentated linear contrails fraction [-]189 REAL, INTENT(INOUT), DIMENSION(klon) :: dzsed_circont ! sedimentated contrails cirrus height [m]190 REAL, INTENT(INOUT), DIMENSION(klon) :: flsed_circont ! sedimentated ice flux in contrails cirrus [kg/s/m2]191 REAL, INTENT(INOUT), DIMENSION(klon) :: cfsed_circont ! sedimentated contrails cirrus fraction [-]192 177 REAL, INTENT(INOUT), DIMENSION(klon) :: flauto ! autoconverted ice flux [kg/s/m2] 193 178 ! 194 179 ! Input for aviation 195 180 ! 196 REAL, INTENT(IN) , DIMENSION(klon) :: lincontfra_in ! input linear contrails fraction [-] 197 REAL, INTENT(IN) , DIMENSION(klon) :: circontfra_in ! input contrail cirrus fraction [-] 198 REAL, INTENT(IN) , DIMENSION(klon) :: qtl_in ! input linear contrails total specific humidity [kg/kg] 199 REAL, INTENT(IN) , DIMENSION(klon) :: qtc_in ! input contrail cirrus total specific humidity [kg/kg] 181 REAL, INTENT(IN) , DIMENSION(klon) :: contfra_in ! input contrails fraction [-] 182 REAL, INTENT(IN) , DIMENSION(klon) :: qtc_in ! input contrails total specific humidity [kg/kg] 183 REAL, INTENT(IN) , DIMENSION(klon) :: nic_in ! input contrails ice crystals concentration [#/kg] 200 184 REAL, INTENT(IN) , DIMENSION(klon) :: flight_dist ! aviation distance flown concentration [m/s/m3] 201 REAL, INTENT(IN) , DIMENSION(klon) :: flight_h2o ! aviation emitted H2O concentration [kgH2O/s/m3] 185 REAL, INTENT(IN) , DIMENSION(klon) :: flight_fuel ! aviation fuel consumption concentration [kg/s/m3] 186 REAL, INTENT(IN) , DIMENSION(klon) :: dzsed_cont_abv! sedimentated contrails height IN THE LAYER ABOVE [m] 187 REAL, INTENT(IN) , DIMENSION(klon) :: flsed_cont_abv! sedimentated ice flux in contrails FROM THE LAYER ABOVE [kg/s/m2] 188 REAL, INTENT(IN) , DIMENSION(klon) :: Nflsed_cont_abv! sedimentated ice crystals concentration in contrails FROM THE LAYER ABOVE [#/s/m2] 189 REAL, INTENT(IN) , DIMENSION(klon) :: cfsed_cont_abv! contrails fraction IN THE LAYER ABOVE [-] 190 REAL, INTENT(INOUT), DIMENSION(klon) :: dzsed_cont ! sedimentated contrails height [m] 191 REAL, INTENT(INOUT), DIMENSION(klon) :: flsed_cont ! sedimentated ice flux in contrails [kg/s/m2] 192 REAL, INTENT(INOUT), DIMENSION(klon) :: Nflsed_cont ! sedimentated ice crystals concentration flux in contrails [#/s/m2] 193 REAL, INTENT(INOUT), DIMENSION(klon) :: cfsed_cont ! sedimentated contrails fraction [-] 202 194 ! 203 195 ! Output … … 238 230 ! NB. idem for the INOUT 239 231 ! 240 REAL, INTENT(INOUT), DIMENSION(klon) :: lincontfra ! linear contrail fraction [-] 241 REAL, INTENT(INOUT), DIMENSION(klon) :: circontfra ! contrail cirrus fraction [-] 242 REAL, INTENT(INOUT), DIMENSION(klon) :: qlincont ! linear contrail specific humidity [kg/kg] 243 REAL, INTENT(INOUT), DIMENSION(klon) :: qcircont ! contrail cirrus specific humidity [kg/kg] 232 REAL, INTENT(INOUT), DIMENSION(klon) :: contfra ! contrail fraction [-] 233 REAL, INTENT(INOUT), DIMENSION(klon) :: qcont ! contrail specific humidity [kg/kg] 234 REAL, INTENT(INOUT), DIMENSION(klon) :: Ncont ! contrail ice crystals concentration [#/kg] 244 235 REAL, INTENT(INOUT), DIMENSION(klon) :: Tcritcont ! critical temperature for contrail formation [K] 245 236 REAL, INTENT(INOUT), DIMENSION(klon) :: qcritcont ! critical specific humidity for contrail formation [kg/kg] 246 237 REAL, INTENT(INOUT), DIMENSION(klon) :: potcontfraP ! potential persistent contrail fraction [-] 247 238 REAL, INTENT(INOUT), DIMENSION(klon) :: potcontfraNP ! potential non-persistent contrail fraction [-] 248 REAL, INTENT(INOUT), DIMENSION(klon) :: dcfl_ini ! linear contrails cloud fraction tendency because of initial formation [s-1] 249 REAL, INTENT(INOUT), DIMENSION(klon) :: dqil_ini ! linear contrails ice specific humidity tendency because of initial formation [kg/kg/s] 250 REAL, INTENT(INOUT), DIMENSION(klon) :: dqtl_ini ! linear contrails total specific humidity tendency because of initial formation [kg/kg/s] 251 REAL, INTENT(INOUT), DIMENSION(klon) :: dcfl_sub ! linear contrails cloud fraction tendency because of sublimation [s-1] 252 REAL, INTENT(INOUT), DIMENSION(klon) :: dqil_sub ! linear contrails ice specific humidity tendency because of sublimation [kg/kg/s] 253 REAL, INTENT(INOUT), DIMENSION(klon) :: dqtl_sub ! linear contrails total specific humidity tendency because of sublimation [kg/kg/s] 254 REAL, INTENT(INOUT), DIMENSION(klon) :: dcfl_cir ! linear contrails cloud fraction tendency because of conversion in cirrus [s-1] 255 REAL, INTENT(INOUT), DIMENSION(klon) :: dqtl_cir ! linear contrails total specific humidity tendency because of conversion in cirrus [kg/kg/s] 256 REAL, INTENT(INOUT), DIMENSION(klon) :: dcfl_mix ! linear contrails cloud fraction tendency because of mixing [s-1] 257 REAL, INTENT(INOUT), DIMENSION(klon) :: dqil_mix ! linear contrails ice specific humidity tendency because of mixing [kg/kg/s] 258 REAL, INTENT(INOUT), DIMENSION(klon) :: dqtl_mix ! linear contrails total specific humidity tendency because of mixing [kg/kg/s] 259 REAL, INTENT(INOUT), DIMENSION(klon) :: dcfl_sed ! linear contrails cloud fraction tendency because of ice sedimentation [s-1] 260 REAL, INTENT(INOUT), DIMENSION(klon) :: dqil_sed ! linear contrails ice specific humidity tendency because of ice sedimentation [kg/kg/s] 261 REAL, INTENT(INOUT), DIMENSION(klon) :: dqtl_sed ! linear contrails total specific humidity tendency because of ice sedimentation [kg/kg/s] 262 REAL, INTENT(INOUT), DIMENSION(klon) :: dcfl_auto ! linear contrails cloud fraction tendency because of ice autoconversion [s-1] 263 REAL, INTENT(INOUT), DIMENSION(klon) :: dqil_auto ! linear contrails ice specific humidity tendency because of ice autoconversion [kg/kg/s] 264 REAL, INTENT(INOUT), DIMENSION(klon) :: dqtl_auto ! linear contrails total specific humidity tendency because of ice autoconversion [kg/kg/s] 265 REAL, INTENT(INOUT), DIMENSION(klon) :: dcfc_sub ! contrail cirrus cloud fraction tendency because of sublimation [s-1] 266 REAL, INTENT(INOUT), DIMENSION(klon) :: dqic_sub ! contrail cirrus ice specific humidity tendency because of sublimation [kg/kg/s] 267 REAL, INTENT(INOUT), DIMENSION(klon) :: dqtc_sub ! contrail cirrus total specific humidity tendency because of sublimation [kg/kg/s] 268 REAL, INTENT(INOUT), DIMENSION(klon) :: dcfc_mix ! contrail cirrus cloud fraction tendency because of mixing [s-1] 269 REAL, INTENT(INOUT), DIMENSION(klon) :: dqic_mix ! contrail cirrus ice specific humidity tendency because of mixing [kg/kg/s] 270 REAL, INTENT(INOUT), DIMENSION(klon) :: dqtc_mix ! contrail cirrus total specific humidity tendency because of mixing [kg/kg/s] 271 REAL, INTENT(INOUT), DIMENSION(klon) :: dcfc_sed ! contrail cirrus cloud fraction tendency because of ice sedimentation [s-1] 272 REAL, INTENT(INOUT), DIMENSION(klon) :: dqic_sed ! contrail cirrus ice specific humidity tendency because of ice sedimentation [kg/kg/s] 273 REAL, INTENT(INOUT), DIMENSION(klon) :: dqtc_sed ! contrail cirrus total specific humidity tendency because of ice sedimentation [kg/kg/s] 274 REAL, INTENT(INOUT), DIMENSION(klon) :: dcfc_auto ! contrail cirrus cloud fraction tendency because of ice autoconversion [s-1] 275 REAL, INTENT(INOUT), DIMENSION(klon) :: dqic_auto ! contrail cirrus ice specific humidity tendency because of ice autoconversion [kg/kg/s] 276 REAL, INTENT(INOUT), DIMENSION(klon) :: dqtc_auto ! contrail cirrus total specific humidity tendency because of ice autoconversion [kg/kg/s] 239 REAL, INTENT(INOUT), DIMENSION(klon) :: AEI_contrails ! Apparent emission index contrails [#/kg] 240 REAL, INTENT(INOUT), DIMENSION(klon) :: AEI_surv_contrails ! Apparent emission index contrails after vortex loss [#/kg] 241 REAL, INTENT(INOUT), DIMENSION(klon) :: fsurv_contrails ! Survival fraction after vortex loss [-] 242 REAL, INTENT(INOUT), DIMENSION(klon) :: section_contrails ! Cross section of newly formed contrails [m2] 243 REAL, INTENT(INOUT), DIMENSION(klon) :: dcfc_ini ! contrails cloud fraction tendency because of initial formation [s-1] 244 REAL, INTENT(INOUT), DIMENSION(klon) :: dqic_ini ! contrails ice specific humidity tendency because of initial formation [kg/kg/s] 245 REAL, INTENT(INOUT), DIMENSION(klon) :: dqtc_ini ! contrails total specific humidity tendency because of initial formation [kg/kg/s] 246 REAL, INTENT(INOUT), DIMENSION(klon) :: dnic_ini ! contrails ice crystals concentration tendency because of initial formation [#/kg/s] 247 REAL, INTENT(INOUT), DIMENSION(klon) :: dcfc_sub ! contrails cloud fraction tendency because of sublimation [s-1] 248 REAL, INTENT(INOUT), DIMENSION(klon) :: dqic_sub ! contrails ice specific humidity tendency because of sublimation [kg/kg/s] 249 REAL, INTENT(INOUT), DIMENSION(klon) :: dqtc_sub ! contrails total specific humidity tendency because of sublimation [kg/kg/s] 250 REAL, INTENT(INOUT), DIMENSION(klon) :: dnic_sub ! contrails ice crystals concentration tendency because of sublimation [#/kg/s] 251 REAL, INTENT(INOUT), DIMENSION(klon) :: dcfc_mix ! contrails cloud fraction tendency because of mixing [s-1] 252 REAL, INTENT(INOUT), DIMENSION(klon) :: dqic_mix ! contrails ice specific humidity tendency because of mixing [kg/kg/s] 253 REAL, INTENT(INOUT), DIMENSION(klon) :: dqtc_mix ! contrails total specific humidity tendency because of mixing [kg/kg/s] 254 REAL, INTENT(INOUT), DIMENSION(klon) :: dnic_mix ! contrails ice crystals concentration tendency because of mixing [#/kg/s] 255 REAL, INTENT(INOUT), DIMENSION(klon) :: dnic_agg ! contrails ice crystals concentration tendency because of aggregation [#/kg/s] 256 REAL, INTENT(INOUT), DIMENSION(klon) :: dcfc_sed ! contrails cloud fraction tendency because of ice sedimentation [s-1] 257 REAL, INTENT(INOUT), DIMENSION(klon) :: dqic_sed ! contrails ice specific humidity tendency because of ice sedimentation [kg/kg/s] 258 REAL, INTENT(INOUT), DIMENSION(klon) :: dqtc_sed ! contrails total specific humidity tendency because of ice sedimentation [kg/kg/s] 259 REAL, INTENT(INOUT), DIMENSION(klon) :: dnic_sed ! contrails ice crystals concentration tendency because of ice sedimentation [#/kg/s] 260 REAL, INTENT(INOUT), DIMENSION(klon) :: dcfc_auto ! contrails cloud fraction tendency because of ice autoconversion [s-1] 261 REAL, INTENT(INOUT), DIMENSION(klon) :: dqic_auto ! contrails ice specific humidity tendency because of ice autoconversion [kg/kg/s] 262 REAL, INTENT(INOUT), DIMENSION(klon) :: dqtc_auto ! contrails total specific humidity tendency because of ice autoconversion [kg/kg/s] 263 REAL, INTENT(INOUT), DIMENSION(klon) :: dnic_auto ! contrails ice crystals concentration tendency because of ice autoconversion [#/kg/s] 277 264 ! 278 265 ! Local … … 308 295 ! for sedimentation and autoconversion 309 296 REAL, DIMENSION(klon) :: qised_abv 310 REAL :: iwc, icefall_velo, coef_sed, qice_sedim 297 REAL :: iwc, icefall_velo, coef_sed, qice_sedim, Nice_sedim 311 298 REAL :: sedfra_abv, sedfra1, sedfra2, sedfra3, sedfra_tmp 312 299 REAL :: dcf_sed1, dcf_sed2, dqvc_sed1, dqvc_sed2, dqi_sed1, dqi_sed2 313 300 REAL :: dz_auto, coef_auto, qice_auto 301 ! 302 ! for aggregation 303 REAL :: dei, sticking_coef 304 REAL :: gamma_agg, dispersion_coef 314 305 ! 315 306 ! for mixing … … 324 315 ! 325 316 ! for cell properties 326 REAL :: rho, rhodz, dz 317 REAL, DIMENSION(klon) :: rho 318 REAL :: rhodz, dz 327 319 ! 328 320 ! for contrails 329 REAL :: contrails_conversion_factor 330 REAL, DIMENSION(klon) :: qised_lincont_abv, qised_circont_abv 331 REAL :: dcfl_sed1, dcfl_sed2, dqtl_sed1, dqtl_sed2, dqil_sed1, dqil_sed2 332 REAL :: dcfc_sed1, dcfc_sed2, dqtc_sed1, dqtc_sed2, dqic_sed1, dqic_sed2 333 REAL :: dz_auto_lincont, dz_auto_circont 321 REAL :: mice, Niceincld 322 REAL, DIMENSION(klon) :: qised_cont_abv, Nised_cont_abv 323 REAL :: dcfc_sed1, dcfc_sed2, dqtc_sed1, dnic_sed1, dqtc_sed2, dqic_sed1, dqic_sed2, dnic_sed2 324 REAL :: dz_auto_cont 334 325 335 326 qzero(:) = 0. … … 365 356 ENDIF 366 357 IF ( ok_plane_contrail ) THEN 367 lincontfra(i) = 0. 368 circontfra(i) = 0. 369 qlincont(i) = 0. 370 qcircont(i) = 0. 358 contfra(i) = 0. 359 qcont(i) = 0. 371 360 IF ( ok_ice_sedim ) THEN 372 flsed_lincont(i) = 0. 373 flsed_circont(i) = 0. 374 dzsed_lincont(i) = 0. 375 dzsed_circont(i) = 0. 376 cfsed_lincont(i) = 0. 377 cfsed_circont(i) = 0. 361 flsed_cont(i) = 0. 362 Nflsed_cont(i) = 0. 363 dzsed_cont(i) = 0. 364 cfsed_cont(i) = 0. 378 365 ENDIF 379 366 ENDIF … … 461 448 !--Initialisation of the cell properties 462 449 !--Dry density [kg/m3] 463 rho = pplay(i) / temp(i) / RD450 rho(i) = pplay(i) / temp(i) / RD 464 451 !--Dry air mass [kg/m2] 465 452 rhodz = ( paprsdn(i) - paprsup(i) ) / RG 466 453 !--Cell thickness [m] 467 dz = rhodz / rho 454 dz = rhodz / rho(i) 468 455 469 456 !--If ok_unadjusted_clouds is set to TRUE, then the saturation adjustment … … 537 524 !--NB. the greater kappa_depsub, the more efficient is the 538 525 !--deposition/sublimation process 539 kappa_depsub = 4. * RPI * capa_cond_cirrus * N_ice_volume / rho * corr_incld_depsub &540 / qsat(i) / ( 4. / 3. * RPI * N_ice_volume / rho * rho_ice )**(1./3.) &526 kappa_depsub = 4. * RPI * capa_cond_cirrus * N_ice_volume / rho(i) * corr_incld_depsub & 527 / qsat(i) / ( 4. / 3. * RPI * N_ice_volume / rho(i) * rho_ice )**(1./3.) & 541 528 / ( RV * temp(i) / water_vapor_diff / pres_sat & 542 529 + RLSTT / air_thermal_conduct / temp(i) * ( RLSTT / RV / temp(i) - 1. ) ) … … 544 531 !--If contrails are activated 545 532 IF ( ok_plane_contrail ) THEN 546 lincontfra(i) = MAX(0., lincontfra_in(i))547 circontfra(i) = MAX(0., circontfra_in(i))548 qlincont(i) = MAX(0., qtl_in(i))549 qcircont(i) = MAX(0., qtc_in(i))550 533 !--The following barriers are needed since the advection scheme does not 551 534 !--conserve order relations 552 mixed_fraction = lincontfra(i) + circontfra(i) 553 IF ( mixed_fraction .GT. cldfra(i) ) THEN 554 mixed_fraction = cldfra(i) / mixed_fraction 555 lincontfra(i) = lincontfra(i) * mixed_fraction 556 circontfra(i) = circontfra(i) * mixed_fraction 557 ENDIF 558 mixed_fraction = qlincont(i) + qcircont(i) 559 IF ( mixed_fraction .GT. qcld(i) ) THEN 560 mixed_fraction = qcld(i) / mixed_fraction 561 qlincont(i) = qlincont(i) * mixed_fraction 562 qcircont(i) = qcircont(i) * mixed_fraction 563 ENDIF 564 565 qised_lincont_abv(i) = 0. 566 IF ( dzsed_lincont_abv(i) .GT. eps ) THEN 535 contfra(i) = MAX(0., MIN(cldfra(i), contfra_in(i))) 536 qcont(i) = MAX(0., MIN(qcld(i), qtc_in(i))) 537 Ncont(i) = MAX(0., nic_in(i)) 538 539 qised_cont_abv(i) = 0. 540 Nised_cont_abv(i) = 0. 541 IF ( dzsed_cont_abv(i) .GT. eps ) THEN 567 542 !--If ice sedimentation is activated, the quantity of sedimentated ice was added 568 543 !--to the total water vapor in the precipitation routine. Here we remove it … … 571 546 !--inacurracies in the advection scheme might lead to considering this water 572 547 !--to be cloudy (hence not in the clear sky region) 573 qised_ lincont_abv(i) = MIN(qclr(i), flsed_lincont_abv(i) &548 qised_cont_abv(i) = MIN(qclr(i), flsed_cont_abv(i) & 574 549 / ( paprsdn(i) - paprsup(i) ) * RG * dtime) 575 qclr(i) = qclr(i) - qised_lincont_abv(i) 550 Nised_cont_abv(i) = MIN(qclr(i), Nflsed_cont_abv(i) & 551 / ( paprsdn(i) - paprsup(i) ) * RG * dtime) 552 qclr(i) = qclr(i) - qised_cont_abv(i) 576 553 ENDIF 577 554 578 qised_circont_abv(i) = 0. 579 IF ( dzsed_circont_abv(i) .GT. eps ) THEN 580 !--If ice sedimentation is activated, the quantity of sedimentated ice was added 581 !--to the total water vapor in the precipitation routine. Here we remove it 582 !--(it will be reincluded later) 583 !--The barrier is needed because although the sedimented ice was vaporised, 584 !--inacurracies in the advection scheme might lead to considering this water 585 !--to be cloudy (hence not in the clear sky region) 586 qised_circont_abv(i) = MIN(qclr(i), flsed_circont_abv(i) & 587 / ( paprsdn(i) - paprsup(i) ) * RG * dtime) 588 qclr(i) = qclr(i) - qised_circont_abv(i) 589 ENDIF 590 591 dcfl_ini(i) = 0. 592 dqil_ini(i) = 0. 593 dqtl_ini(i) = 0. 594 dcfl_sub(i) = 0. 595 dqil_sub(i) = 0. 596 dqtl_sub(i) = 0. 597 dcfl_cir(i) = 0. 598 dqtl_cir(i) = 0. 599 dcfl_mix(i) = 0. 600 dqil_mix(i) = 0. 601 dqtl_mix(i) = 0. 602 dcfl_sed(i) = 0. 603 dqil_sed(i) = 0. 604 dqtl_sed(i) = 0. 555 dcfc_ini(i) = 0. 556 dqic_ini(i) = 0. 557 dqtc_ini(i) = 0. 558 dnic_ini(i) = 0. 605 559 dcfc_sub(i) = 0. 606 560 dqic_sub(i) = 0. 607 561 dqtc_sub(i) = 0. 562 dnic_sub(i) = 0. 608 563 dcfc_mix(i) = 0. 609 564 dqic_mix(i) = 0. 610 565 dqtc_mix(i) = 0. 566 dnic_mix(i) = 0. 567 dnic_agg(i) = 0. 611 568 dcfc_sed(i) = 0. 612 569 dqic_sed(i) = 0. 613 570 dqtc_sed(i) = 0. 571 dnic_sed(i) = 0. 572 dcfc_auto(i) = 0. 573 dqic_auto(i) = 0. 574 dqtc_auto(i) = 0. 575 dnic_auto(i) = 0. 614 576 ELSE 615 lincontfra(i) = 0. 616 circontfra(i) = 0. 617 qlincont(i) = 0. 618 qcircont(i) = 0. 577 contfra(i) = 0. 578 qcont(i) = 0. 619 579 ENDIF 620 580 … … 624 584 !---------------------------------------------------------------------- 625 585 626 !--If there is a linearcontrail627 IF ( lincontfra(i) .GT. eps ) THEN586 !--If there is a contrail 587 IF ( contfra(i) .GT. eps ) THEN 628 588 !--We remove contrails from the main class 629 cldfra(i) = MAX(0., cldfra(i) - lincontfra(i))630 qcld(i) = MAX(0., qcld(i) - q lincont(i))631 qvc(i) = MAX(0., MIN(qcld(i), qvc(i) - qsat(i) * lincontfra(i)))589 cldfra(i) = MAX(0., cldfra(i) - contfra(i)) 590 qcld(i) = MAX(0., qcld(i) - qcont(i)) 591 qvc(i) = MAX(0., MIN(qcld(i), qvc(i) - qsat(i) * contfra(i))) 632 592 633 593 !--The contrail is always adjusted to saturation 634 qiceincld = ( qlincont(i) / lincontfra(i) - qsat(i) ) 594 qiceincld = ( qcont(i) / contfra(i) - qsat(i) ) 595 !--The in-cloud ice crystals concentration is conserved 596 Niceincld = Ncont(i) / contfra(i) 635 597 636 598 !--If the ice water content is too low, the cloud is purely sublimated 637 599 IF ( qiceincld .LT. eps ) THEN 638 dcfl_sub(i) = - lincontfra(i) 639 dqil_sub(i) = - qiceincld * lincontfra(i) 640 dqtl_sub(i) = - qlincont(i) 600 dcfc_sub(i) = - contfra(i) 601 dqic_sub(i) = - qiceincld * contfra(i) 602 dqtc_sub(i) = - qcont(i) 603 dnic_sub(i) = - Ncont(i) 641 604 642 605 !--Only a part of the contrail is sublimated … … 649 612 650 613 !--Tendencies and diagnostics 651 dcfl_sub(i) = - lincontfra(i) * pdf_e1 652 dqil_sub(i) = - lincontfra(i) * pdf_e2 / pdf_shape 653 dqtl_sub(i) = dqil_sub(i) + dcfl_sub(i) * qsat(i) 614 dcfc_sub(i) = - contfra(i) * pdf_e1 615 dqic_sub(i) = - contfra(i) * pdf_e2 / pdf_shape 616 dqtc_sub(i) = dqic_sub(i) + dcfc_sub(i) * qsat(i) 617 dnic_sub(i) = dcfc_sub(i) * Niceincld 654 618 ENDIF ! qiceincld .LT. eps 655 619 656 620 !--Add tendencies 657 lincontfra(i) = lincontfra(i) + dcfl_sub(i) 658 qlincont(i) = qlincont(i) + dqtl_sub(i) 659 clrfra(i) = clrfra(i) - dcfl_sub(i) 660 qclr(i) = qclr(i) - dqtl_sub(i) 661 ENDIF ! lincontfra(i) .GT. eps 662 663 !--If there is a contrail cirrus 664 IF ( circontfra(i) .GT. eps ) THEN 665 !--We remove contrails from the main class 666 cldfra(i) = MAX(0., cldfra(i) - circontfra(i)) 667 qcld(i) = MAX(0., qcld(i) - qcircont(i)) 668 qvc(i) = MAX(0., MIN(qcld(i), qvc(i) - qsat(i) * circontfra(i))) 669 670 !--The contrail is always adjusted to saturation 671 qiceincld = ( qcircont(i) / circontfra(i) - qsat(i) ) 672 673 !--If the ice water content is too low, the cloud is purely sublimated 674 IF ( qiceincld .LT. eps ) THEN 675 dcfc_sub(i) = - circontfra(i) 676 dqic_sub(i) = - qiceincld * circontfra(i) 677 dqtc_sub(i) = - qcircont(i) 678 679 !--Only a part of the contrail is sublimated 680 ELSE 681 !--Gamma distribution starting at qvapincld (everything that is below qiceincld_min) 682 pdf_shape = nu_iwc_pdf_lscp / qiceincld 683 pdf_y = pdf_shape * qiceincld_min(i) 684 pdf_e1 = GAMMAINC ( nu_iwc_pdf_lscp , pdf_y ) 685 pdf_e2 = GAMMAINC ( nu_iwc_pdf_lscp + 1. , pdf_y ) 686 687 !--Tendencies and diagnostics 688 dcfc_sub(i) = - circontfra(i) * pdf_e1 689 dqic_sub(i) = - circontfra(i) * pdf_e2 / pdf_shape 690 dqtc_sub(i) = dqic_sub(i) + dcfc_sub(i) * qsat(i) 691 ENDIF ! qiceincld .LT. eps 692 693 !--Add tendencies 694 circontfra(i) = circontfra(i) + dcfc_sub(i) 695 qcircont(i) = qcircont(i) + dqtc_sub(i) 696 clrfra(i) = clrfra(i) - dcfc_sub(i) 697 qclr(i) = qclr(i) - dqtc_sub(i) 698 ENDIF ! circontfra(i) .GT. eps 621 contfra(i) = contfra(i) + dcfc_sub(i) 622 qcont(i) = qcont(i) + dqtc_sub(i) 623 Ncont(i) = Ncont(i) + dnic_sub(i) 624 clrfra(i) = clrfra(i) - dcfc_sub(i) 625 qclr(i) = qclr(i) - dqtc_sub(i) 626 ENDIF ! contfra(i) .GT. eps 699 627 700 628 … … 1118 1046 !-------------------------------------- 1119 1047 1120 IF ( ( lincontfra(i) .GT. eps ) .AND. ( clrfra(i) .GT. eps ) ) THEN1048 IF ( ( contfra(i) .GT. eps ) .AND. ( clrfra(i) .GT. eps ) ) THEN 1121 1049 1122 1050 !-- PART 1 - TURBULENT DIFFUSION … … 1138 1066 !-- clouds_perim = P * N_cld_mix 1139 1067 !-- 1140 bovera = aspect_ratio_ lincontrails1068 bovera = aspect_ratio_contrails 1141 1069 !--P / a is a function of b / a only, that we can calculate 1142 1070 !-- P / a = RPI * ( 3. * ( 1. + b / a ) - SQRT( (3. + b / a) * (1. + 3. * b / a) ) ) … … 1149 1077 !-- cldfra = a * b / a * RPI / (P / a) * coef_mix_lscp * cldfra * ( 1. - cldfra ) 1150 1078 !-- a = (P / a) / ( coef_mix_lscp * RPI * ( 1. - cldfra ) * (b / a) ) 1151 a_mix = Povera / coef_mixing_ lincontrails / RPI / ( 1. - lincontfra(i) ) / bovera1079 a_mix = Povera / coef_mixing_contrails / RPI / ( 1. - contfra(i) ) / bovera 1152 1080 !--and finally, 1153 1081 !-- N_cld_mix = coef_mix_lscp * cldfra * ( 1. - cldfra ) * cell_area / ( P / a * a ) 1154 N_cld_mix = coef_mixing_ lincontrails * lincontfra(i) * ( 1. - lincontfra(i) ) &1082 N_cld_mix = coef_mixing_contrails * contfra(i) * ( 1. - contfra(i) ) & 1155 1083 * cell_area(i) / Povera / a_mix 1156 1084 … … 1182 1110 !--With this, the clouds increase in size along b only, by a factor 1183 1111 !--L_shear * SQRT(2.) / 2. (to account for the 45 degrees orientation of the wind) 1184 L_shear = coef_shear_ lincontrails * shear(i) * dz * dtime1112 L_shear = coef_shear_contrails * shear(i) * dz * dtime 1185 1113 !--therefore, the fraction of clear sky mixed is 1186 1114 !-- N_cld_mix * ( a * (b + L_shear * SQRT(2.) / 2.) - a * b ) * RPI / 2. / cell_area … … 1189 1117 !--(note that they are equal) 1190 1118 shear_fra = RPI * L_shear * a_mix * SQRT(2.) / 2. / 2. * N_cld_mix / cell_area(i) 1119 !shear_fra = RPI * L_shear * a_mix * bovera / 2. * N_cld_mix / cell_area(i) 1191 1120 !--and the environment and cloud mixed fractions are the same, 1192 1121 !--which we add to the previous calculated mixed fractions. … … 1200 1129 1201 1130 clrfra_mix = MIN(clrfra(i), clrfra_mix) 1202 cldfra_mix = MIN( lincontfra(i), cldfra_mix)1131 cldfra_mix = MIN(contfra(i), cldfra_mix) 1203 1132 1204 1133 !--We compute the limit vapor in clear sky where the mixed cloud could not … … 1209 1138 !--cloud's vapor without condensing or sublimating ice crystals 1210 1139 qvapinclr_lim = qsat(i) * ( 1. + cldfra_mix / clrfra_mix ) & 1211 - q lincont(i) / lincontfra(i) * cldfra_mix / clrfra_mix1140 - qcont(i) / contfra(i) * cldfra_mix / clrfra_mix 1212 1141 1213 1142 IF ( qvapinclr_lim .LT. 0. ) THEN 1214 !--Whatever we do, the cloud will increase in size 1215 !--If the linear contrail increases in size, the increment is considered 1216 !--to be a contrail cirrus 1217 dcfl_mix(i) = dcfl_mix(i) + clrfra_mix 1218 dqtl_mix(i) = dqtl_mix(i) + clrfra_mix * qclr(i) / clrfra(i) 1219 dqil_mix(i) = dqil_mix(i) + clrfra_mix * ( qclr(i) / clrfra(i) - qsat(i) ) 1143 !--Whatever we do, contrails will increase in size 1144 dcfc_mix(i) = dcfc_mix(i) + clrfra_mix 1145 dqtc_mix(i) = dqtc_mix(i) + clrfra_mix * qclr(i) / clrfra(i) 1146 dqic_mix(i) = dqic_mix(i) + clrfra_mix * ( qclr(i) / clrfra(i) - qsat(i) ) 1220 1147 ELSE 1221 1148 !--We then calculate the clear sky part where the humidity is lower than … … 1253 1180 !--decreases the size of the clouds 1254 1181 !--For aviation, we increase the chance that the air surrounding contrails 1255 !--is moist. This is quantified with chi_mixing_ lincontrails1182 !--is moist. This is quantified with chi_mixing_contrails 1256 1183 sigma_mix = MAX(1. - pdf_fra_below_lim / clrfra_mix, & 1257 1184 MIN(pdf_fra_above_lim / clrfra_mix, & 1258 chi_mixing_lincontrails * pdf_fra_above_lim & 1259 / ( pdf_fra_below_lim + chi_mixing_lincontrails * pdf_fra_above_lim ))) 1260 1261 IF ( pdf_fra_above_lim .GT. eps ) THEN 1262 dcfl_mix(i) = dcfl_mix(i) + clrfra_mix * sigma_mix 1263 dqtl_mix(i) = dqtl_mix(i) + clrfra_mix * sigma_mix & 1264 * pdf_q_above_lim / pdf_fra_above_lim 1265 dqil_mix(i) = dqil_mix(i) + clrfra_mix * sigma_mix & 1266 * ( pdf_q_above_lim / pdf_fra_above_lim - qsat(i) ) 1267 !--If the linear contrail increases in size, the increment is considered 1268 !--to be a contrail cirrus 1269 !qvapinmix = ( pdf_q_above_lim / pdf_fra_above_lim * clrfra_mix & 1270 ! + qlincont(i) / lincontfra(i) * cldfra_mix ) & 1271 ! / ( clrfra_mix + cldfra_mix ) 1272 !qiceinmix = ( qlincont(i) / lincontfra(i) - qsat(i) ) * cldfra_mix & 1273 ! / ( clrfra_mix + cldfra_mix ) 1274 !dcfc_mix(i) = dcfc_mix(i) + clrfra_mix * sigma_mix 1275 !dqtc_mix(i) = dqtc_mix(i) + clrfra_mix * sigma_mix * qvapinmix 1276 !dqtl_mix(i) = dqtl_mix(i) - cldfra_mix * sigma_mix & 1277 ! * ( qlincont(i) / lincontfra(i) - qvapinmix ) 1278 !dqic_mix(i) = dqic_mix(i) + clrfra_mix * sigma_mix * qiceinmix 1279 !dqil_mix(i) = dqil_mix(i) - cldfra_mix * sigma_mix & 1280 ! * ( qlincont(i) / lincontfra(i) - qsat(i) - qiceinmix ) 1281 ENDIF 1282 1283 IF ( pdf_fra_below_lim .GT. eps ) THEN 1284 dcfl_mix(i) = dcfl_mix(i) - cldfra_mix * ( 1. - sigma_mix ) 1285 dqtl_mix(i) = dqtl_mix(i) - cldfra_mix * ( 1. - sigma_mix ) & 1286 * qlincont(i) / lincontfra(i) 1287 dqil_mix(i) = dqil_mix(i) - cldfra_mix * ( 1. - sigma_mix ) & 1288 * ( qlincont(i) / lincontfra(i) - qsat(i) ) 1289 ENDIF 1290 1291 ENDIF 1292 ENDIF ! ( lincontfra(i) .GT. eps ) .AND. ( clrfra(i) .GT. eps ) 1293 1294 IF ( ( circontfra(i) .GT. eps ) .AND. ( clrfra(i) .GT. eps ) ) THEN 1295 1296 !-- PART 1 - TURBULENT DIFFUSION 1297 1298 !--Clouds within the mesh are assumed to be ellipses. The length of the 1299 !--semi-major axis is a and the length of the semi-minor axis is b. 1300 !--N_cld_mix is the number of clouds in contact with clear sky, and can be non-integer. 1301 !--In particular, it is 0 if cldfra = 1. 1302 !--clouds_perim is the total perimeter of the clouds within the mesh, 1303 !--not considering interfaces with other meshes (only the interfaces with clear 1304 !--sky are taken into account). 1305 !-- 1306 !--The area of each cloud is A = a * b * RPI, 1307 !--and the perimeter of each cloud is 1308 !-- P ~= RPI * ( 3 * (a + b) - SQRT( (3 * a + b) * (a + 3 * b) ) ) 1309 !-- 1310 !--With cell_area the area of the cell, we have: 1311 !-- cldfra = A * N_cld_mix / cell_area 1312 !-- clouds_perim = P * N_cld_mix 1313 !-- 1314 bovera = aspect_ratio_cirrus 1315 !--P / a is a function of b / a only, that we can calculate 1316 !-- P / a = RPI * ( 3. * ( 1. + b / a ) - SQRT( (3. + b / a) * (1. + 3. * b / a) ) ) 1317 Povera = RPI * ( 3. * (1. + bovera) - SQRT( (3. + bovera) * (1. + 3. * bovera) ) ) 1318 1319 !--The clouds perimeter is imposed using the formula from Morcrette 2012, 1320 !--based on observations. 1321 !-- clouds_perim / cell_area = N_cld_mix * ( P / a * a ) / cell_area = coef_mix_lscp * cldfra * ( 1. - cldfra ) 1322 !--With cldfra = a * ( b / a * a ) * RPI * N_cld_mix / cell_area, we have: 1323 !-- cldfra = a * b / a * RPI / (P / a) * coef_mix_lscp * cldfra * ( 1. - cldfra ) 1324 !-- a = (P / a) / ( coef_mix_lscp * RPI * ( 1. - cldfra ) * (b / a) ) 1325 a_mix = Povera / coef_mixing_lscp / RPI / ( 1. - circontfra(i) ) / bovera 1326 !--and finally, 1327 !-- N_cld_mix = coef_mix_lscp * cldfra * ( 1. - cldfra ) * cell_area / ( P / a * a ) 1328 N_cld_mix = coef_mixing_lscp * circontfra(i) * ( 1. - circontfra(i) ) & 1329 * cell_area(i) / Povera / a_mix 1330 1331 !--The time required for turbulent diffusion to homogenize a region of size 1332 !--L_mix is defined as (L_mix**2/tke_dissip)**(1./3.) (Pope, 2000; Field et al., 2014) 1333 !--We compute L_mix and assume that the cloud is mixed over this length 1334 L_mix = SQRT( dtime**3 * pbl_eps(i) ) 1335 !--The mixing length cannot be greater than the semi-minor axis. In this case, 1336 !--the entire cloud is mixed. 1337 L_mix = MIN(L_mix, a_mix * bovera) 1338 1339 !--The fraction of clear sky mixed is 1340 !-- N_cld_mix * ( (a + L_mix) * (b + L_mix) - a * b ) * RPI / cell_area 1341 clrfra_mix = N_cld_mix * RPI / cell_area(i) & 1342 * ( a_mix * ( 1. + bovera ) * L_mix + L_mix**2 ) 1343 !--The fraction of clear sky mixed is 1344 !-- N_cld_mix * ( a * b - (a - L_mix) * (b - L_mix) ) * RPI / cell_area 1345 cldfra_mix = N_cld_mix * RPI / cell_area(i) & 1346 * ( a_mix * ( 1. + bovera ) * L_mix - L_mix**2 ) 1347 1348 1349 !-- PART 2 - SHEARING 1350 1351 !--The clouds are then sheared. We keep the shape and number 1352 !--assumptions from before. The clouds are sheared along their 1353 !--semi-major axis (a_mix), on the entire cell heigh dz. 1354 !--The increase in size is 1355 L_shear = coef_shear_lscp * shear(i) * dz * dtime 1356 !--therefore, the fraction of clear sky mixed is 1357 !-- N_cld_mix * ( (a + L_shear) * b - a * b ) * RPI / 2. / cell_area 1358 !--and the fraction of cloud mixed is 1359 !-- N_cld_mix * ( (a * b) - (a - L_shear) * b ) * RPI / 2. / cell_area 1360 !--(note that they are equal) 1361 shear_fra = RPI * L_shear * a_mix * bovera / 2. * N_cld_mix / cell_area(i) 1362 !--and the environment and cloud mixed fractions are the same, 1363 !--which we add to the previous calculated mixed fractions. 1364 !--We therefore assume that the sheared clouds and the turbulent 1365 !--mixed clouds are different. 1366 clrfra_mix = clrfra_mix + shear_fra 1367 cldfra_mix = cldfra_mix + shear_fra 1368 1369 1370 !-- PART 3 - CALCULATION OF THE MIXING PROPERTIES 1371 1372 clrfra_mix = MIN(clrfra(i), clrfra_mix) 1373 cldfra_mix = MIN(circontfra(i), cldfra_mix) 1374 1375 !--We compute the limit vapor in clear sky where the mixed cloud could not 1376 !--survive if all the ice crystals were sublimated. Note that here we assume, 1377 !--for growth or reduction of the cloud, that the mixed cloud is adjusted 1378 !--to saturation, ie the vapor in the mixed cloud is qsat. This is only a 1379 !--diagnostic, and if the cloud size is increased, we add the new vapor to the 1380 !--cloud's vapor without condensing or sublimating ice crystals 1381 qvapinclr_lim = qsat(i) * ( 1. + cldfra_mix / clrfra_mix ) & 1382 - qcircont(i) / circontfra(i) * cldfra_mix / clrfra_mix 1383 1384 IF ( qvapinclr_lim .LT. 0. ) THEN 1385 !--Whatever we do, the cloud will increase in size 1386 dcfc_mix(i) = dcfc_mix(i) + clrfra_mix 1387 dqtc_mix(i) = dqtc_mix(i) + clrfra_mix * qclr(i) / clrfra(i) 1388 dqic_mix(i) = dqic_mix(i) + clrfra_mix * ( qclr(i) / clrfra(i) - qsat(i) ) 1389 ELSE 1390 !--We then calculate the clear sky part where the humidity is lower than 1391 !--qvapinclr_lim, and the part where it is higher than qvapinclr_lim 1392 !--This is the clear-sky PDF calculated in the condensation section. Note 1393 !--that if we are here, we necessarily went through the condensation part 1394 !--because the clear sky fraction can only be reduced by condensation. 1395 !--Thus the `pdf_xxx` variables are well defined. 1396 1397 rhl_clr = qclr(i) / clrfra(i) / qsatl(i) * 100. 1398 pdf_x = qvapinclr_lim / qsatl(i) * 100. 1399 pdf_loc = rhl_clr - pdf_scale(i) * pdf_gamma(i) 1400 pdf_x = qsat(i) / qsatl(i) * 100. 1401 pdf_y = LOG( MAX( ( pdf_x - pdf_loc ) / pdf_scale(i), eps) ) * pdf_alpha(i) 1402 IF ( pdf_y .GT. 10. ) THEN !--Avoid overflows 1403 pdf_fra_above_lim = 0. 1404 pdf_q_above_lim = 0. 1405 ELSEIF ( pdf_y .LT. -10. ) THEN 1406 pdf_fra_above_lim = clrfra(i) 1407 pdf_q_above_lim = qclr(i) 1408 ELSE 1409 pdf_y = EXP( pdf_y ) 1410 pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha(i) , pdf_y ) 1411 pdf_e3 = pdf_scale(i) * ( 1. - pdf_e3 ) * pdf_gamma(i) 1412 pdf_fra_above_lim = EXP( - pdf_y ) * clrfra(i) 1413 pdf_q_above_lim = ( pdf_e3 * clrfra(i) & 1414 + pdf_loc * pdf_fra_above_lim ) * qsatl(i) / 100. 1415 ENDIF 1416 1417 pdf_fra_below_lim = clrfra(i) - pdf_fra_above_lim 1418 1419 !--sigma_mix is the ratio of the surroundings of the clouds where mixing 1420 !--increases the size of the cloud, to the total surroundings of the clouds. 1421 !--This implies that ( 1. - sigma_mix ) quantifies the ratio where mixing 1422 !--decreases the size of the clouds 1423 sigma_mix = MAX(1. - pdf_fra_below_lim / clrfra_mix, & 1424 MIN(pdf_fra_above_lim / clrfra_mix, & 1425 chi_mixing * pdf_fra_above_lim & 1426 / ( pdf_fra_below_lim + chi_mixing * pdf_fra_above_lim ))) 1185 chi_mixing_contrails * pdf_fra_above_lim & 1186 / ( pdf_fra_below_lim + chi_mixing_contrails * pdf_fra_above_lim ))) 1427 1187 1428 1188 IF ( pdf_fra_above_lim .GT. eps ) THEN … … 1437 1197 dcfc_mix(i) = dcfc_mix(i) - cldfra_mix * ( 1. - sigma_mix ) 1438 1198 dqtc_mix(i) = dqtc_mix(i) - cldfra_mix * ( 1. - sigma_mix ) & 1439 * qc ircont(i) / circontfra(i)1199 * qcont(i) / contfra(i) 1440 1200 dqic_mix(i) = dqic_mix(i) - cldfra_mix * ( 1. - sigma_mix ) & 1441 * ( qcircont(i) / circontfra(i) - qsat(i) ) 1201 * ( qcont(i) / contfra(i) - qsat(i) ) 1202 dnic_mix(i) = dnic_mix(i) - cldfra_mix * ( 1. - sigma_mix ) & 1203 * ( Ncont(i) / contfra(i) ) 1442 1204 ENDIF 1443 1205 1444 1206 ENDIF 1445 ENDIF ! ( c ircontfra(i) .GT. eps ) .AND. ( clrfra(i) .GT. eps )1446 1447 IF ( ( lincontfra(i) + circontfra(i)) .GT. eps ) THEN1207 ENDIF ! ( contfra(i) .GT. eps ) .AND. ( clrfra(i) .GT. eps ) 1208 1209 IF ( contfra(i) .GT. eps ) THEN 1448 1210 !--We balance the mixing tendencies between the different cloud classes 1449 mixed_fraction = dcf_mix(i) + dcf l_mix(i) + dcfc_mix(i)1211 mixed_fraction = dcf_mix(i) + dcfc_mix(i) 1450 1212 IF ( mixed_fraction .GT. clrfra(i) ) THEN 1451 1213 mixed_fraction = clrfra(i) / mixed_fraction … … 1453 1215 dqvc_mix(i) = dqvc_mix(i) * mixed_fraction 1454 1216 dqi_mix(i) = dqi_mix(i) * mixed_fraction 1455 dcfl_mix(i) = dcfl_mix(i) * mixed_fraction1456 dqtl_mix(i) = dqtl_mix(i) * mixed_fraction1457 dqil_mix(i) = dqil_mix(i) * mixed_fraction1458 1217 dcfc_mix(i) = dcfc_mix(i) * mixed_fraction 1459 1218 dqtc_mix(i) = dqtc_mix(i) * mixed_fraction 1460 1219 dqic_mix(i) = dqic_mix(i) * mixed_fraction 1220 dnic_mix(i) = dnic_mix(i) * mixed_fraction 1461 1221 ENDIF 1462 1222 1463 IF ( lincontfra(i) .GT. eps ) THEN 1464 !--Add tendencies 1465 lincontfra(i) = lincontfra(i) + dcfl_mix(i) 1466 qlincont(i) = qlincont(i) + dqtl_mix(i) 1467 clrfra(i) = clrfra(i) - dcfl_mix(i) 1468 qclr(i) = qclr(i) - dqtl_mix(i) 1469 ENDIF 1470 1471 IF ( circontfra(i) .GT. eps ) THEN 1472 !--Add tendencies 1473 circontfra(i) = circontfra(i) + dcfc_mix(i) 1474 qcircont(i) = qcircont(i) + dqtc_mix(i) 1475 clrfra(i) = clrfra(i) - dcfc_mix(i) 1476 qclr(i) = qclr(i) - dqtc_mix(i) 1477 ENDIF 1223 !--Add tendencies 1224 contfra(i) = contfra(i) + dcfc_mix(i) 1225 qcont(i) = qcont(i) + dqtc_mix(i) 1226 Ncont(i) = Ncont(i) + dnic_mix(i) 1227 clrfra(i) = clrfra(i) - dcfc_mix(i) 1228 qclr(i) = qclr(i) - dqtc_mix(i) 1478 1229 ENDIF 1479 1230 … … 1486 1237 1487 1238 1239 !---------------------------------------- 1240 !-- ICE CRYSTALS AGGREGATION -- 1241 !---------------------------------------- 1242 1243 !--Aggregation of ice crystals only occur for 2-moment microphysical clouds, 1244 !--i.e., contrails for now 1245 IF ( ok_plane_contrail ) THEN 1246 IF ( contfra(i) .GT. eps ) THEN 1247 mice = qcont(i) / MAX(eps, Ncont(i)) 1248 icefall_velo = ICECRYSTAL_VELO(mice, temp(i), pplay(i)) 1249 !--Volumic radius (in meters) 1250 dei = (mice / rho_ice * 3. / 4. / RPI)**(1./3.) 1251 !--Effective radius (in meters) 1252 dei = MIN(1e-4, MAX(1e-6, dei / eff2vol_radius_contrails)) 1253 !--Effective radius to effective diameter 1254 dei = 8. / 3. / SQRT(3.) * dei 1255 !--Lin 1983's formula 1256 sticking_coef = EXP( 0.025 * MIN( ( temp(i) - RTT ), 0.) ) 1257 gamma_agg = 1. 1258 dispersion_coef = 0.1 1259 dnic_agg(i) = - 0.5 * gamma_agg * RPI / 6. * sticking_coef * dispersion_coef & 1260 * dei**2 * icefall_velo * ( Ncont(i) / contfra(i) )**2 1261 !--Gridbox average 1262 dnic_agg(i) = dnic_agg(i) * contfra(i) 1263 1264 !--Add tendencies 1265 Ncont(i) = Ncont(i) + dnic_agg(i) 1266 ENDIF 1267 ENDIF 1268 1269 1270 !--------------------------------------- 1271 !-- ICE SEDIMENTATION -- 1272 !--------------------------------------- 1273 1488 1274 IF ( ok_ice_sedim ) THEN 1489 1275 !--Initialisation 1490 1276 dz_auto = 0. 1491 dz_auto_lincont = 0. 1492 dz_auto_circont = 0. 1277 dz_auto_cont = 0. 1493 1278 dcf_sed1 = 0. 1494 1279 dqvc_sed1 = 0. … … 1497 1282 dqvc_sed2 = 0. 1498 1283 dqi_sed2 = 0. 1499 dcfl_sed1 = 0.1500 dqtl_sed1 = 0.1501 dqil_sed1 = 0.1502 dcfl_sed2 = 0.1503 dqtl_sed2 = 0.1504 dqil_sed2 = 0.1505 1284 dcfc_sed1 = 0. 1506 1285 dqtc_sed1 = 0. 1507 1286 dqic_sed1 = 0. 1287 dnic_sed1 = 0. 1508 1288 dcfc_sed2 = 0. 1509 1289 dqtc_sed2 = 0. 1510 1290 dqic_sed2 = 0. 1511 !--------------------------------------- 1512 !-- ICE SEDIMENTATION -- 1513 !--------------------------------------- 1291 dnic_sed2 = 0. 1514 1292 ! 1515 1293 !--First, the current ice is sedimentated into the layer below. The ice fallspeed … … 1519 1297 ! 1520 1298 IF ( cldfra(i) .GT. eps ) THEN 1521 iwc = rho * ( qcld(i) - qvc(i) ) / cldfra(i)1299 iwc = rho(i) * ( qcld(i) - qvc(i) ) / cldfra(i) 1522 1300 icefall_velo = fallice_sedim * cice_velo * MAX(0., iwc)**dice_velo 1523 1301 … … 1555 1333 ENDIF 1556 1334 ! 1557 IF ( lincontfra(i) .GT. eps ) THEN 1558 icefall_velo = fallice_linear_contrails 1335 IF ( contfra(i) .GT. eps ) THEN 1336 mice = qcont(i) / MAX(eps, Ncont(i)) 1337 icefall_velo = ICECRYSTAL_VELO(mice, temp(i), pplay(i)) 1559 1338 1560 1339 !--Sedimentation 1561 1340 coef_sed = MAX(0., MIN(1., 2. - icefall_velo * dtime / dz)) 1562 dzsed_lincont(i) = MIN(dz, icefall_velo * dtime) * coef_sed 1563 cfsed_lincont(i) = lincontfra(i) 1564 qice_sedim = ( qlincont(i) - qsat(i) * lincontfra(i) ) * dzsed_lincont(i) / dz 1341 dzsed_cont(i) = MIN(dz, icefall_velo * dtime) * coef_sed 1342 cfsed_cont(i) = contfra(i) 1343 qice_sedim = ( qcont(i) - qsat(i) * contfra(i) ) * dzsed_cont(i) / dz 1344 Nice_sedim = Ncont(i) * dzsed_cont(i) / dz 1565 1345 !--Tendencies 1566 dcfl_sed1 = - lincontfra(i) * dzsed_lincont(i) / dz 1567 dqil_sed1 = - qice_sedim 1568 dqtl_sed1 = - qlincont(i) * dzsed_lincont(i) / dz 1346 dcfc_sed1 = - contfra(i) * dzsed_cont(i) / dz 1347 dqic_sed1 = - qice_sedim 1348 dqtc_sed1 = - qcont(i) * dzsed_cont(i) / dz 1349 dnic_sed1 = - Nice_sedim 1569 1350 !--Convert to flux 1570 flsed_lincont(i) = qice_sedim * ( paprsdn(i) - paprsup(i) ) / RG / dtime 1351 flsed_cont(i) = qice_sedim * ( paprsdn(i) - paprsup(i) ) / RG / dtime 1352 Nflsed_cont(i) = Nice_sedim * ( paprsdn(i) - paprsup(i) ) / RG / dtime 1571 1353 1572 1354 !--Autoconversion 1573 1355 coef_auto = MAX(0., MIN(1., icefall_velo * dtime / dz - 1.)) 1574 dz_auto_ lincont = MIN(dz, icefall_velo * dtime) * coef_auto1575 qice_auto = ( q lincont(i) - qsat(i) * lincontfra(i) ) * dz_auto_lincont / dz1356 dz_auto_cont = MIN(dz, icefall_velo * dtime) * coef_auto 1357 qice_auto = ( qcont(i) - qsat(i) * contfra(i) ) * dz_auto_cont / dz 1576 1358 !--Tendencies 1577 dcfl_auto(i) = - lincontfra(i) * dz_auto_lincont / dz 1578 dqil_auto(i) = - qice_auto 1579 dqtl_auto(i) = - qlincont(i) * dz_auto_lincont / dz 1359 dcfc_auto(i) = - contfra(i) * dz_auto_cont / dz 1360 dqic_auto(i) = - qice_auto 1361 dqtc_auto(i) = - qcont(i) * dz_auto_cont / dz 1362 dnic_auto(i) = - Ncont(i) * dz_auto_cont / dz 1580 1363 !--Convert to flux 1581 1364 flauto(i) = flauto(i) + qice_auto * ( paprsdn(i) - paprsup(i) ) / RG / dtime 1582 1365 1583 1366 !--Add tendencies 1584 lincontfra(i) = lincontfra(i) + dcfl_sed1 + dcfl_auto(i) 1585 qlincont(i) = qlincont(i) + dqtl_sed1 + dqtl_auto(i) 1586 clrfra(i) = clrfra(i) - dcfl_sed1 - dcfl_auto(i) 1587 qclr(i) = qclr(i) - dqtl_sed1 - dqtl_auto(i) 1588 ENDIF 1589 ! 1590 IF ( circontfra(i) .GT. eps ) THEN 1591 !icefall_velo = fallice_cirrus_contrails 1592 iwc = rho * ( qcircont(i) / circontfra(i) - qsat(i) ) 1593 icefall_velo = fallice_sedim * cice_velo * MAX(0., iwc)**dice_velo 1594 1595 !--Sedimentation 1596 coef_sed = MAX(0., MIN(1., 2. - icefall_velo * dtime / dz)) 1597 dzsed_circont(i) = MIN(dz, icefall_velo) * dtime * coef_sed 1598 cfsed_circont(i) = circontfra(i) 1599 qice_sedim = ( qcircont(i) - qsat(i) * circontfra(i) ) * dzsed_circont(i) / dz 1600 !--Tendencies 1601 dcfc_sed1 = - circontfra(i) * dzsed_circont(i) / dz 1602 dqic_sed1 = - qice_sedim 1603 dqtc_sed1 = - qcircont(i) * dzsed_circont(i) / dz 1604 !--Convert to flux 1605 flsed_circont(i) = qice_sedim * ( paprsdn(i) - paprsup(i) ) / RG / dtime 1606 1607 !--Autoconversion 1608 coef_auto = MAX(0., MIN(1., icefall_velo * dtime / dz - 1.)) 1609 dz_auto_circont = MIN(dz, icefall_velo * dtime) * coef_auto 1610 qice_auto = ( qcircont(i) - qsat(i) * circontfra(i) ) * dz_auto_circont / dz 1611 !--Tendencies 1612 dcfc_auto(i) = - circontfra(i) * dz_auto_circont / dz 1613 dqic_auto(i) = - qice_auto 1614 dqtc_auto(i) = - qcircont(i) * dz_auto_circont / dz 1615 !--Convert to flux 1616 flauto(i) = flauto(i) + qice_auto * ( paprsdn(i) - paprsup(i) ) / RG / dtime 1617 1618 !--Add tendencies 1619 circontfra(i) = circontfra(i) + dcfc_sed1 + dcfc_auto(i) 1620 qcircont(i) = qcircont(i) + dqtc_sed1 + dqtc_auto(i) 1367 contfra(i) = contfra(i) + dcfc_sed1 + dcfc_auto(i) 1368 qcont(i) = qcont(i) + dqtc_sed1 + dqtc_auto(i) 1369 Ncont(i) = Ncont(i) + dnic_sed1 + dnic_auto(i) 1621 1370 clrfra(i) = clrfra(i) - dcfc_sed1 - dcfc_auto(i) 1622 1371 qclr(i) = qclr(i) - dqtc_sed1 - dqtc_auto(i) … … 1707 1456 ENDIF ! clrfra(i) .GT. eps 1708 1457 1709 !--If the sedimented ice falls into a linear contrail, it increases its content 1710 IF ( lincontfra(i) .GT. eps ) THEN 1711 sedfra_tmp = sedfra1 * lincontfra(i) / (totfra_in(i) - cldfra(i)) 1712 dqil_sed2 = dqil_sed2 + sedfra_tmp * qiceincld 1713 ENDIF 1714 1715 !--If the sedimented ice falls into a contrail cirrus, it increases its content 1716 IF ( circontfra(i) .GT. eps ) THEN 1717 sedfra_tmp = sedfra1 * circontfra(i) / (totfra_in(i) - cldfra(i)) 1458 !--If the sedimented ice falls into a contrail, it increases its content 1459 IF ( contfra(i) .GT. eps ) THEN 1460 sedfra_tmp = sedfra1 * contfra(i) / (totfra_in(i) - cldfra(i)) 1718 1461 dqic_sed2 = dqic_sed2 + sedfra_tmp * qiceincld 1719 1462 ENDIF … … 1732 1475 1733 1476 IF ( ok_plane_contrail ) THEN 1734 sedfra_abv = MIN(dzsed_ lincont_abv(i), dz) / dz * cfsed_lincont_abv(i)1477 sedfra_abv = MIN(dzsed_cont_abv(i), dz) / dz * cfsed_cont_abv(i) 1735 1478 IF ( sedfra_abv .GT. eps ) THEN 1736 1479 … … 1744 1487 !--(3) we increase the cloud fraction but use in-cloud water vapor as 1745 1488 !--background vapor 1746 sedfra2 = MIN(cfsed_ lincont(i), cfsed_lincont_abv(i)) &1747 * MAX(0., MIN(dz, dzsed_ lincont_abv(i)) &1748 - dzsed_ lincont(i) - dz_auto_lincont) / dz1749 sedfra3 = MIN(cfsed_ lincont(i), cfsed_lincont_abv(i)) &1750 * MIN(MIN(dz, dzsed_ lincont_abv(i)), &1751 dzsed_ lincont(i) + dz_auto_lincont) / dz1489 sedfra2 = MIN(cfsed_cont(i), cfsed_cont_abv(i)) & 1490 * MAX(0., MIN(dz, dzsed_cont_abv(i)) & 1491 - dzsed_cont(i) - dz_auto_cont) / dz 1492 sedfra3 = MIN(cfsed_cont(i), cfsed_cont_abv(i)) & 1493 * MIN(MIN(dz, dzsed_cont_abv(i)), & 1494 dzsed_cont(i) + dz_auto_cont) / dz 1752 1495 sedfra1 = sedfra_abv - sedfra3 - sedfra2 1753 1496 1754 qiceincld = qised_lincont_abv(i) / sedfra_abv 1497 qiceincld = qised_cont_abv(i) / sedfra_abv 1498 Niceincld = Nised_cont_abv(i) / sedfra_abv 1755 1499 1756 1500 !--For region (1) … … 1787 1531 sedfra1 * chi_sedim * pdf_fra_above_lim & 1788 1532 / ( pdf_fra_below_lim + chi_sedim * pdf_fra_above_lim ))) & 1789 * clrfra(i) / (totfra_in(i) - lincontfra(i)) 1790 1791 IF ( pdf_fra_above_lim .GT. eps ) THEN 1792 dcfl_sed2 = dcfl_sed2 + sedfra_tmp 1793 dqtl_sed2 = dqtl_sed2 + sedfra_tmp & 1794 * (qiceincld + pdf_q_above_lim / pdf_fra_above_lim) 1795 dqil_sed2 = dqil_sed2 + sedfra_tmp & 1796 * (qiceincld + pdf_q_above_lim / pdf_fra_above_lim - qsat(i)) 1797 ENDIF 1798 ENDIF ! clrfra(i) .GT. eps 1799 1800 !--If the sedimented ice falls into a natural cirrus, it increases its content 1801 IF ( cldfra(i) .GT. eps ) THEN 1802 sedfra_tmp = sedfra1 * cldfra(i) / (totfra_in(i) - lincontfra(i)) 1803 dqi_sed2 = dqi_sed2 + sedfra_tmp * qiceincld 1804 ENDIF 1805 1806 !--If the sedimented ice falls into a contrail cirrus, it increases its content 1807 IF ( circontfra(i) .GT. eps ) THEN 1808 sedfra_tmp = sedfra1 * circontfra(i) / (totfra_in(i) - lincontfra(i)) 1809 dqic_sed2 = dqic_sed2 + sedfra_tmp * qiceincld 1810 ENDIF 1811 ENDIF ! sedfra1 .GT. eps 1812 1813 !--For region (2) 1814 dqil_sed2 = dqil_sed2 + sedfra2 * qiceincld 1815 1816 !--For region (3) 1817 IF ( lincontfra(i) .GT. eps ) THEN 1818 dcfl_sed2 = dcfl_sed2 + sedfra3 1819 dqtl_sed2 = dqtl_sed2 + sedfra3 * (qsat(i) + qiceincld) 1820 dqil_sed2 = dqil_sed2 + sedfra3 * qiceincld 1821 ENDIF 1822 ENDIF ! sedfra_abv .GT. eps 1823 1824 ! 1825 sedfra_abv = MIN(dzsed_circont_abv(i), dz) / dz * cfsed_circont_abv(i) 1826 IF ( sedfra_abv .GT. eps ) THEN 1827 1828 !--Three regions to be defined: (1) clear-sky, (2) cloudy, and 1829 !--(3) region where it was previously cloudy but now clear-sky because of 1830 !--ice sedimentation 1831 !--(1) we use the clear-sky PDF to determine the humidity and increase the 1832 !--cloud fraction / sublimate falling ice. NB it can also fall into 1833 !--another cloud class 1834 !--(2) we do not increase cloud fraction and add the falling ice to the cloud 1835 !--(3) we increase the cloud fraction but use in-cloud water vapor as 1836 !--background vapor 1837 sedfra2 = MIN(cfsed_circont(i), cfsed_circont_abv(i)) & 1838 * MAX(0., MIN(dz, dzsed_circont_abv(i)) & 1839 - dzsed_circont(i) - dz_auto_circont) / dz 1840 sedfra3 = MIN(cfsed_circont(i), cfsed_circont_abv(i)) & 1841 * MIN(MIN(dz, dzsed_circont_abv(i)), & 1842 dzsed_circont(i) + dz_auto_circont) / dz 1843 sedfra1 = sedfra_abv - sedfra3 - sedfra2 1844 1845 qiceincld = qised_circont_abv(i) / sedfra_abv 1846 1847 !--For region (1) 1848 IF ( sedfra1 .GT. eps ) THEN 1849 IF ( clrfra(i) .GT. eps ) THEN 1850 qvapinclr_lim = qsat(i) - qiceincld 1851 1852 rhl_clr = qclr(i) / clrfra(i) / qsatl(i) * 100. 1853 pdf_x = qvapinclr_lim / qsatl(i) * 100. 1854 pdf_loc = rhl_clr - pdf_scale(i) * pdf_gamma(i) 1855 pdf_x = qsat(i) / qsatl(i) * 100. 1856 pdf_y = LOG( MAX( ( pdf_x - pdf_loc ) / pdf_scale(i), eps) ) * pdf_alpha(i) 1857 IF ( pdf_y .GT. 10. ) THEN !--Avoid overflows 1858 pdf_fra_above_lim = 0. 1859 pdf_q_above_lim = 0. 1860 ELSEIF ( pdf_y .LT. -10. ) THEN 1861 pdf_fra_above_lim = clrfra(i) 1862 pdf_q_above_lim = qclr(i) 1863 ELSE 1864 pdf_y = EXP( pdf_y ) 1865 pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha(i) , pdf_y ) 1866 pdf_e3 = pdf_scale(i) * ( 1. - pdf_e3 ) * pdf_gamma(i) 1867 pdf_fra_above_lim = EXP( - pdf_y ) * clrfra(i) 1868 pdf_q_above_lim = ( pdf_e3 * clrfra(i) & 1869 + pdf_loc * pdf_fra_above_lim ) * qsatl(i) / 100. 1870 ENDIF 1871 1872 pdf_fra_below_lim = clrfra(i) - pdf_fra_above_lim 1873 1874 !--The first three lines allow to favor ISSR rather than subsaturated sky for 1875 !--the growth. The fourth line indicates that the ice is falling only in 1876 !--the clear-sky area. The rest falls into other cloud classes 1877 sedfra_tmp = MAX(sedfra1 - pdf_fra_below_lim, MIN(pdf_fra_above_lim, & 1878 sedfra1 * chi_sedim * pdf_fra_above_lim & 1879 / ( pdf_fra_below_lim + chi_sedim * pdf_fra_above_lim ))) & 1880 * clrfra(i) / (totfra_in(i) - circontfra(i)) 1533 * clrfra(i) / (totfra_in(i) - contfra(i)) 1881 1534 1882 1535 IF ( pdf_fra_above_lim .GT. eps ) THEN … … 1886 1539 dqic_sed2 = dqic_sed2 + sedfra_tmp & 1887 1540 * (qiceincld + pdf_q_above_lim / pdf_fra_above_lim - qsat(i)) 1541 dnic_sed2 = dnic_sed2 + sedfra_tmp * Niceincld 1888 1542 ENDIF 1889 1543 ENDIF ! clrfra(i) .GT. eps … … 1891 1545 !--If the sedimented ice falls into a natural cirrus, it increases its content 1892 1546 IF ( cldfra(i) .GT. eps ) THEN 1893 sedfra_tmp = sedfra1 * cldfra(i) / (totfra_in(i) - c ircontfra(i))1547 sedfra_tmp = sedfra1 * cldfra(i) / (totfra_in(i) - contfra(i)) 1894 1548 dqi_sed2 = dqi_sed2 + sedfra_tmp * qiceincld 1895 1549 ENDIF 1896 1897 !--If the sedimented ice falls into a contrail cirrus, it increases its content1898 IF ( lincontfra(i) .GT. eps ) THEN1899 sedfra_tmp = sedfra1 * lincontfra(i) / (totfra_in(i) - circontfra(i))1900 dqil_sed2 = dqil_sed2 + sedfra_tmp * qiceincld1901 ENDIF1902 1550 ENDIF ! sedfra1 .GT. eps 1903 1551 1904 1552 !--For region (2) 1905 1553 dqic_sed2 = dqic_sed2 + sedfra2 * qiceincld 1554 dnic_sed2 = dnic_sed2 + sedfra2 * Niceincld 1906 1555 1907 1556 !--For region (3) 1908 IF ( c ircontfra(i) .GT. eps ) THEN1557 IF ( contfra(i) .GT. eps ) THEN 1909 1558 dcfc_sed2 = dcfc_sed2 + sedfra3 1910 1559 dqtc_sed2 = dqtc_sed2 + sedfra3 * (qsat(i) + qiceincld) 1911 1560 dqic_sed2 = dqic_sed2 + sedfra3 * qiceincld 1561 dnic_sed2 = dnic_sed2 + sedfra3 * Niceincld 1912 1562 ENDIF 1913 1563 ENDIF ! sedfra_abv .GT. eps … … 1929 1579 if ( ok_plane_contrail ) THEN 1930 1580 !--Add tendencies 1931 lincontfra(i) = lincontfra(i) + dcfl_sed2 1932 qlincont(i) = qlincont(i) + dqtl_sed2 1933 clrfra(i) = clrfra(i) - dcfl_sed2 1934 qclr(i) = qclr(i) - dqtl_sed2 1581 contfra(i) = contfra(i) + dcfc_sed2 1582 qcont(i) = qcont(i) + dqtc_sed2 1583 Ncont(i) = Ncont(i) + dnic_sed2 1584 clrfra(i) = clrfra(i) - dcfc_sed2 1585 qclr(i) = qclr(i) - dqtc_sed2 1935 1586 !--We re-include sublimated sedimentated ice into clear sky water vapor 1936 qclr(i) = qclr(i) + qised_lincont_abv(i) 1937 !--Diagnostics 1938 dcfl_sed(i) = dcfl_sed1 + dcfl_sed2 1939 dqtl_sed(i) = dqtl_sed1 + dqtl_sed2 1940 dqil_sed(i) = dqil_sed1 + dqil_sed2 1941 1942 !--Add tendencies 1943 circontfra(i) = circontfra(i) + dcfc_sed2 1944 qcircont(i) = qcircont(i) + dqtc_sed2 1945 clrfra(i) = clrfra(i) - dcfc_sed2 1946 qclr(i) = qclr(i) - dqtc_sed2 1947 !--We re-include sublimated sedimentated ice into clear sky water vapor 1948 qclr(i) = qclr(i) + qised_circont_abv(i) 1587 qclr(i) = qclr(i) + qised_cont_abv(i) 1949 1588 !--Diagnostics 1950 1589 dcfc_sed(i) = dcfc_sed1 + dcfc_sed2 1951 1590 dqtc_sed(i) = dqtc_sed1 + dqtc_sed2 1952 1591 dqic_sed(i) = dqic_sed1 + dqic_sed2 1592 dnic_sed(i) = dnic_sed1 + dnic_sed2 1953 1593 ENDIF 1954 1594 … … 1957 1597 1958 1598 !--We put back contrails in the clouds class 1959 IF ( ( lincontfra(i) + circontfra(i)) .GT. 0. ) THEN1960 cldfra(i) = cldfra(i) + lincontfra(i) + circontfra(i)1961 qcld(i) = qcld(i) + q lincont(i) + qcircont(i)1962 qvc(i) = qvc(i) + qsat(i) * ( lincontfra(i) + circontfra(i))1599 IF ( contfra(i) .GT. 0. ) THEN 1600 cldfra(i) = cldfra(i) + contfra(i) 1601 qcld(i) = qcld(i) + qcont(i) 1602 qvc(i) = qvc(i) + qsat(i) * contfra(i) 1963 1603 ENDIF 1964 1604 … … 2042 1682 IF ( ok_plane_contrail ) THEN 2043 1683 1684 rho = pplay / temp / RD 1685 2044 1686 CALL contrails_formation( & 2045 klon, dtime, pplay, temp, qsat, qsatl, gamma_cond, &2046 flight_dist, clrfra, qclr, pdf_scale, pdf_alpha, pdf_gamma, &1687 klon, dtime, pplay, temp, rho, qsat, qsatl, gamma_cond, & 1688 flight_dist, flight_fuel, clrfra, qclr, pdf_scale, pdf_alpha, pdf_gamma, & 2047 1689 keepgoing, pt_pron_clds, & 2048 1690 Tcritcont, qcritcont, potcontfraP, potcontfraNP, & 2049 dcfl_ini, dqil_ini, dqtl_ini) 1691 AEI_contrails, AEI_surv_contrails, fsurv_contrails, section_contrails, & 1692 dcfc_ini, dqic_ini, dqtc_ini, dnic_ini) 2050 1693 2051 1694 DO i = 1, klon 2052 1695 IF ( keepgoing(i) .AND. pt_pron_clds(i) ) THEN 2053 1696 2054 !--Convert existing contrail fraction into "natural" cirrus cloud fraction 2055 IF ( ( clrfra(i) .LE. eps ) .OR. ( lincontfra(i) .LE. eps ) ) THEN 2056 contrails_conversion_factor = 1. 2057 ELSE 2058 contrails_conversion_factor = ( 1. - & 2059 !--First, the linear contrails are continuously degraded in induced cirrus 2060 EXP( - dtime / linear_contrails_lifetime ) & 2061 !--Second, if the cloud fills the entire gridbox, the linear contrails 2062 !--cannot exist. The exponent is set so that this only happens for 2063 !--very cloudy gridboxes 2064 * ( clrfra(i) / totfra_in(i) )**0.1 ) 2065 ENDIF 2066 dcfl_cir(i) = - contrails_conversion_factor * lincontfra(i) 2067 dqtl_cir(i) = - contrails_conversion_factor * qlincont(i) 2068 2069 dcfl_ini(i) = MIN(dcfl_ini(i), issrfra(i)) 2070 dqtl_ini(i) = MIN(dqtl_ini(i), qissr(i)) 1697 dcfc_ini(i) = MIN(dcfc_ini(i), issrfra(i)) 1698 dqtc_ini(i) = MIN(dqtc_ini(i), qissr(i)) 2071 1699 2072 1700 !--Add tendencies 2073 cldfra(i) = cldfra(i) + dcfl_ini(i) 2074 qcld(i) = qcld(i) + dqtl_ini(i) 2075 qvc(i) = qvc(i) + dcfl_ini(i) * qsat(i) 2076 issrfra(i) = issrfra(i) - dcfl_ini(i) 2077 qissr(i) = qissr(i) - dqtl_ini(i) 2078 clrfra(i) = clrfra(i) - dcfl_ini(i) 2079 qclr(i) = qclr(i) - dqtl_ini(i) 2080 2081 lincontfra(i) = lincontfra(i) + dcfl_cir(i) + dcfl_ini(i) 2082 qlincont(i) = qlincont(i) + dqtl_cir(i) + dqtl_ini(i) 2083 circontfra(i) = circontfra(i) - dcfl_cir(i) 2084 qcircont(i) = qcircont(i) - dqtl_cir(i) 1701 cldfra(i) = cldfra(i) + dcfc_ini(i) 1702 qcld(i) = qcld(i) + dqtc_ini(i) 1703 qvc(i) = qvc(i) + dcfc_ini(i) * qsat(i) 1704 issrfra(i) = issrfra(i) - dcfc_ini(i) 1705 qissr(i) = qissr(i) - dqtc_ini(i) 1706 clrfra(i) = clrfra(i) - dcfc_ini(i) 1707 qclr(i) = qclr(i) - dqtc_ini(i) 1708 1709 contfra(i) = contfra(i) + dcfc_ini(i) 1710 qcont(i) = qcont(i) + dqtc_ini(i) 1711 Ncont(i) = Ncont(i) + dnic_ini(i) 2085 1712 2086 1713 … … 2089 1716 !------------------------------------------- 2090 1717 2091 IF ( (lincontfra(i) .LT. eps) .OR. (qlincont(i) .LT. (qsat(i) * lincontfra(i))) ) THEN 2092 cldfra(i) = cldfra(i) - lincontfra(i) 2093 qcld(i) = qcld(i) - qlincont(i) 2094 qvc(i) = qvc(i) - qsat(i) * lincontfra(i) 2095 lincontfra(i) = 0. 2096 qlincont(i) = 0. 2097 ENDIF 2098 2099 IF ( (circontfra(i) .LT. eps) .OR. (qcircont(i) .LT. (qsat(i) * circontfra(i))) ) THEN 2100 cldfra(i) = cldfra(i) - circontfra(i) 2101 qcld(i) = qcld(i) - qcircont(i) 2102 qvc(i) = qvc(i) - qsat(i) * circontfra(i) 2103 circontfra(i) = 0. 2104 qcircont(i) = 0. 1718 IF ( (contfra(i) .LT. eps) .OR. (qcont(i) .LT. (qsat(i) * contfra(i))) ) THEN 1719 cldfra(i) = cldfra(i) - contfra(i) 1720 qcld(i) = qcld(i) - qcont(i) 1721 qvc(i) = qvc(i) - qsat(i) * contfra(i) 1722 contfra(i) = 0. 1723 qcont(i) = 0. 1724 Ncont(i) = 0. 2105 1725 ENDIF 2106 1726 … … 2110 1730 qcld(i) = 0. 2111 1731 qvc(i) = 0. 2112 lincontfra(i) = 0. 2113 qlincont(i) = 0. 2114 circontfra(i) = 0. 2115 qcircont(i) = 0. 1732 contfra(i)= 0. 1733 qcont(i) = 0. 1734 Ncont(i) = 0. 2116 1735 qincld(i) = qsat(i) 2117 1736 ELSE … … 2119 1738 ENDIF 2120 1739 2121 cldfra(i) = MAX(cldfra(i), lincontfra(i) + circontfra(i))2122 qcld(i) = MAX(qcld(i), q lincont(i) + qcircont(i))2123 qvc(i) = MAX(qvc(i), qsat(i) * ( lincontfra(i) + circontfra(i)))1740 cldfra(i) = MAX(cldfra(i), contfra(i)) 1741 qcld(i) = MAX(qcld(i), qcont(i)) 1742 qvc(i) = MAX(qvc(i), qsat(i) * contfra(i)) 2124 1743 2125 1744 !--Diagnostics 2126 dcfl_ini(i) = dcfl_ini(i) / dtime 2127 dqil_ini(i) = dqil_ini(i) / dtime 2128 dqtl_ini(i) = dqtl_ini(i) / dtime 2129 dcfl_sub(i) = dcfl_sub(i) / dtime 2130 dqil_sub(i) = dqil_sub(i) / dtime 2131 dqtl_sub(i) = dqtl_sub(i) / dtime 2132 dcfl_cir(i) = dcfl_cir(i) / dtime 2133 dqtl_cir(i) = dqtl_cir(i) / dtime 2134 dcfl_mix(i) = dcfl_mix(i) / dtime 2135 dqil_mix(i) = dqil_mix(i) / dtime 2136 dqtl_mix(i) = dqtl_mix(i) / dtime 2137 dcfl_sed(i) = dcfl_sed(i) / dtime 2138 dqil_sed(i) = dqil_sed(i) / dtime 2139 dqtl_sed(i) = dqtl_sed(i) / dtime 2140 dcfl_auto(i) = dcfl_auto(i) / dtime 2141 dqil_auto(i) = dqil_auto(i) / dtime 2142 dqtl_auto(i) = dqtl_auto(i) / dtime 1745 dcfc_ini(i) = dcfc_ini(i) / dtime 1746 dqic_ini(i) = dqic_ini(i) / dtime 1747 dqtc_ini(i) = dqtc_ini(i) / dtime 1748 dnic_ini(i) = dnic_ini(i) / dtime 2143 1749 dcfc_sub(i) = dcfc_sub(i) / dtime 2144 1750 dqic_sub(i) = dqic_sub(i) / dtime 2145 1751 dqtc_sub(i) = dqtc_sub(i) / dtime 1752 dnic_sub(i) = dnic_sub(i) / dtime 2146 1753 dcfc_mix(i) = dcfc_mix(i) / dtime 2147 1754 dqic_mix(i) = dqic_mix(i) / dtime 2148 1755 dqtc_mix(i) = dqtc_mix(i) / dtime 1756 dnic_mix(i) = dnic_mix(i) / dtime 1757 dnic_agg(i) = dnic_agg(i) / dtime 2149 1758 dcfc_sed(i) = dcfc_sed(i) / dtime 2150 1759 dqic_sed(i) = dqic_sed(i) / dtime 2151 1760 dqtc_sed(i) = dqtc_sed(i) / dtime 1761 dnic_sed(i) = dnic_sed(i) / dtime 2152 1762 dcfc_auto(i) = dcfc_auto(i) / dtime 2153 1763 dqic_auto(i) = dqic_auto(i) / dtime 2154 1764 dqtc_auto(i) = dqtc_auto(i) / dtime 1765 dnic_auto(i) = dnic_auto(i) / dtime 2155 1766 2156 1767 ENDIF ! keepgoing
Note: See TracChangeset
for help on using the changeset viewer.
