Ignore:
Timestamp:
Aug 4, 2025, 3:03:07 PM (11 days ago)
Author:
aborella
Message:

Additional diags for contrails + simplified coupling between deep conv and cirrus clouds + small modifsin RRTM for RF of contrails alone

File:
1 edited

Legend:

Unmodified
Added
Removed

  • LMDZ6/branches/contrails/libf/phylmd/lmdz_aviation.f90

    r5790 r5796  
    225225    !--Add a source of contrails from aviation
    226226    IF ( ( potcontfraP(i) .GT. eps ) .AND. ( flight_dist(i) .GT. 1e-20 ) ) THEN
    227       !section_contrails(i) = CONTRAIL_CROSS_SECTION_ONERA()
    228       section_contrails(i) = CONTRAIL_CROSS_SECTION_SCHUMANN( &
     227      section_contrails(i) = CONTRAIL_CROSS_SECTION( &
    229228              dtime, rho(i), flight_dist(i), flight_fuel(i))
    230229      icesat_ratio = qpotcontP / potcontfraP(i) / qsat(i)
    231       !--If Nice init is fixed in the plume
    232       !Nice_per_m_init_contrails = Nice_init_contrails * 1e6 * section_contrails(i)
    233       !--Else if it is parameterized
    234230      CALL CONTRAIL_ICE_NUMBER_CONCENTRATION(temp(i), icesat_ratio, rho(i), &
    235               flight_dist(i), flight_fuel(i), Nice_per_m_init_contrails, &
     231              flight_dist(i), flight_fuel(i), section_contrails(i), &
     232              Nice_per_m_init_contrails, &
    236233              AEI_contrails(i), AEI_surv_contrails(i), fsurv_contrails(i))
    237234
    … …  
    254251
    255252!**********************************************************************************
    256 FUNCTION CONTRAIL_CROSS_SECTION_ONERA()
    257 
    258 USE lmdz_lscp_ini, ONLY: initial_width_contrails, initial_height_contrails
    259 
    260 IMPLICIT NONE
    261 
    262 !
    263 ! Output
    264 !
    265 REAL :: contrail_cross_section_onera ! [m2]
    266 !
    267 ! Local
    268 !
    269 
    270 contrail_cross_section_onera = initial_width_contrails * initial_height_contrails
    271 
    272 END FUNCTION CONTRAIL_CROSS_SECTION_ONERA
    273 
    274 
    275 !**********************************************************************************
    276 !--Based on Schumann (1998)
    277 FUNCTION CONTRAIL_CROSS_SECTION_SCHUMANN(dtime, rho_air, flight_dist, flight_fuel)
    278 
     253FUNCTION CONTRAIL_CROSS_SECTION(dtime, rho_air, flight_dist, flight_fuel)
     254
     255USE lmdz_lscp_ini, ONLY: iflag_cross_sec_contrail
    279256USE lmdz_lscp_ini, ONLY: initial_width_contrails, initial_height_contrails
    280257
    … …  
    291268! Output
    292269!
    293 REAL :: contrail_cross_section_schumann ! [m2]
     270REAL :: contrail_cross_section ! [m2]
    294271!
    295272! Local
    … …  
    297274REAL :: dilution_factor
    298275
    299 !--The contrail is formed on average in the middle of the timestep
    300 dilution_factor = 7000. * (dtime / 2.)**0.8
    301 contrail_cross_section_schumann = flight_fuel / flight_dist * dilution_factor / rho_air
    302 
    303 END FUNCTION CONTRAIL_CROSS_SECTION_SCHUMANN
     276IF ( iflag_cross_sec_contrail .EQ. 0 ) THEN
     277  contrail_cross_section = initial_width_contrails * initial_height_contrails
     278ELSEIF ( iflag_cross_sec_contrail .EQ. 1 ) THEN
     279  !--Based on Schumann (1998)
     280  !--The contrail is formed on average in the middle of the timestep
     281  dilution_factor = 7000. * (dtime / 2.)**0.8
     282  contrail_cross_section = flight_fuel / flight_dist * dilution_factor / rho_air
     283ENDIF
     284
     285END FUNCTION CONTRAIL_CROSS_SECTION
    304286
    305287
    306288!**********************************************************************************
    307289SUBROUTINE CONTRAIL_ICE_NUMBER_CONCENTRATION(temp, icesat_ratio, rho_air, &
    308         flight_dist, flight_fuel, Nice_per_m_init_contrails, &
     290        flight_dist, flight_fuel, cross_section, Nice_per_m_init_contrails, &
    309291        AEI_contrails, AEI_surv_contrails, fsurv_contrails)
    310292
    311 USE lmdz_lscp_ini, ONLY: RPI
     293USE lmdz_lscp_ini, ONLY: RPI, iflag_AEI_contrail, iflag_vortex_loss
     294USE lmdz_lscp_ini, ONLY: Nice_init_contrails, fraction_survival_contrails
    312295USE lmdz_lscp_ini, ONLY: EI_soot_aviation, air_to_fuel_ratio_engine, wingspan
    313296USE lmdz_lscp_ini, ONLY: Naer_amb, raer_amb_mean, raer_amb_std, r_soot_mean, r_soot_std
    … …  
    325308REAL, INTENT(IN)  :: flight_dist  ! flown distance [m/s/m3]
    326309REAL, INTENT(IN)  :: flight_fuel  ! fuel consumed [kg/s/m3]
     310REAL, INTENT(IN)  :: cross_section  ! contrail cross section [m2]
    327311!
    328312! Output
    … …  
    344328REAL :: Nice_per_m, fuel_per_m, frac_surv
    345329
     330! fuel consumption in kg/m flown
     331fuel_per_m = flight_fuel / flight_dist
     332
    346333!------------------------------
    347334!--  INITIAL ICE NUCLEATION  --
    348335!------------------------------
    349336!
    350 !--Karcher et al. (2015), Bier and Burkhardt (2019, 2022)
    351 !log_liqsat_ratio = LOG(liq_satratio)
    352 
    353 !! dry core radius in nm
    354 !! HERE SHOULD IT BE liqsat_ratio OR liqsat_ratio - 1. ?
    355 !rd_act_amb = (4. / 27. / LOG(liqsat_ratio)**2 / 0.5)**(1./3.)
    356 !! Integrate lognormal distribution between rd_act_amb and +inf
    357 !coef_expo = 4. / SQRT(2. * RPI) / LOG(raer_amb_std)
    358 !phi_amb = 1. / (1. + (rd_act_amb / raer_amb_mean)**coef_expo)
    359 !
    360 !! BU22, Eq. A1, dry core radius in nm
    361 !rd_act_soot = 0.96453426 + 1.21152973 / log_liqsat_ratio - 0.00520358 / log_liqsat_ratio**2 &
    362 !        + 2.32286432e-5 / log_liqsat_ratio**3 - 4.36979055e-8 / log_liqsat_ratio**4
    363 !rd_act_soot = MIN(150., MAX(1., rd_act_soot))
    364 !! Integrate lognormal distribution between rd_act_amb and +inf
    365 !coef_expo = 4. / SQRT(2. * RPI) / LOG(r_soot_std)
    366 !phi_amb = 1. / (1. + (rd_act_soot / r_soot_mean)**coef_expo)
    367 !
    368 !dil_coef = (0.01 / t0)**0.9
    369 !
    370 !! BU22, Eq. 5b
    371 !AEI_soot = phi_soot * EI_soot_aviation
    372 !AEI_amb = phi_amb * air_to_fuel_ratio_engine * (1. - dil_coef) / dil_coef &
    373 !        / rho_air * Naer_amb * 1e6
    374 !AEI_contrails = AEI_soot + AEI_amb
    375 AEI_contrails = EI_soot_aviation
    376 
     337IF ( iflag_AEI_contrail .EQ. 0 ) THEN
     338  !--If Nice init is fixed in the plume, in #/cm3
     339  AEI_contrails = Nice_init_contrails * 1e6 * cross_section / fuel_per_m
     340ELSEIF ( iflag_AEI_contrail .EQ. 1 ) THEN
     341  AEI_contrails = EI_soot_aviation
     342!ELSEIF ( iflag_AEI_contrail .EQ. 2 ) THEN
     343!  !--Karcher et al. (2015), Bier and Burkhardt (2019, 2022)
     344!  log_liqsat_ratio = LOG(liq_satratio)
     345
     346!  ! dry core radius in nm
     347!  ! HERE SHOULD IT BE liqsat_ratio OR liqsat_ratio - 1. ?
     348!  rd_act_amb = (4. / 27. / LOG(liqsat_ratio)**2 / 0.5)**(1./3.)
     349!  ! Integrate lognormal distribution between rd_act_amb and +inf
     350!  coef_expo = 4. / SQRT(2. * RPI) / LOG(raer_amb_std)
     351!  phi_amb = 1. / (1. + (rd_act_amb / raer_amb_mean)**coef_expo)
     352! 
     353!  ! BU22, Eq. A1, dry core radius in nm
     354!  rd_act_soot = 0.96453426 + 1.21152973 / log_liqsat_ratio - 0.00520358 / log_liqsat_ratio**2 &
     355!          + 2.32286432e-5 / log_liqsat_ratio**3 - 4.36979055e-8 / log_liqsat_ratio**4
     356!  rd_act_soot = MIN(150., MAX(1., rd_act_soot))
     357!  ! Integrate lognormal distribution between rd_act_amb and +inf
     358!  coef_expo = 4. / SQRT(2. * RPI) / LOG(r_soot_std)
     359!  phi_amb = 1. / (1. + (rd_act_soot / r_soot_mean)**coef_expo)
     360! 
     361!  dil_coef = (0.01 / t0)**0.9
     362! 
     363!  ! BU22, Eq. 5b
     364!  AEI_soot = phi_soot * EI_soot_aviation
     365!  AEI_amb = phi_amb * air_to_fuel_ratio_engine * (1. - dil_coef) / dil_coef &
     366!          / rho_air * Naer_amb * 1e6
     367!  AEI_contrails = AEI_soot + AEI_amb
     368ENDIF
    377369
    378370!----------------------------
    … …  
    383375!--which is an update of Unterstrasser (2016, U16 hereinafter)
    384376
    385 ! fuel consumption in kg/m flown
    386 fuel_per_m = flight_fuel / flight_dist
    387 
    388 ! LU25, Eq. A2
    389 z_atm = 607.46 * (icesat_ratio - 1.)**0.897 * (temp / 205.)**2.225
    390 
    391 ! water vapor emissions [kg/m flown], LU25, Eq. 2
    392 ! U16, Eq. 6
    393 rho_emit = fuel_per_m * EI_H2O_aviation / plume_area_loss
    394 ! LU25, Eq. A3
    395 temp_205 = temp - 205.
    396 z_emit = 1106.6 * (rho_emit * 1e5)**(0.678 + 0.0116 * temp_205) &
    397         * EXP(- (0.0807 + 0.000428 * temp_205) * temp_205)
    398 
    399 ! U16, Eq. 4
    400 z_desc = SQRT(8. * circ_0_loss / RPI / N_Brunt_Vaisala_aviation )
    401 
    402377! initial number of ice crystals [#/m flown], LU25, Eq. 1
    403378Nice_per_m = fuel_per_m * AEI_contrails
    404 ! ice crystals number concentration [#/m3], LU25, Eq. A1
    405 nice_init = Nice_per_m / plume_area_loss
    406 ! LU25, Eq. 9, 13d, 13e, 13f and 13g
    407 z_delta = (nice_init / nice_init_ref_loss)**(-0.16) * (1.27 * z_atm + 0.42 * z_emit) - 0.49 * z_desc
    408 
    409 ! LU25, Eq. 11, 12, 13a, 13b and 13c
    410 fsurv_contrails = 0.42 + 1.31 / RPI * ATAN(-1. + z_delta / 100.)
    411 fsurv_contrails = MIN(1., MAX(0., fsurv_contrails))
     379
     380IF ( iflag_vortex_loss .EQ. 0 ) THEN
     381  fsurv_contrails = fraction_survival_contrails
     382ELSEIF ( iflag_vortex_loss .EQ. 1 ) THEN
     383  ! LU25, Eq. A2
     384  z_atm = 607.46 * (icesat_ratio - 1.)**0.897 * (temp / 205.)**2.225
     385 
     386  ! water vapor emissions [kg/m flown], LU25, Eq. 2
     387  ! U16, Eq. 6
     388  rho_emit = fuel_per_m * EI_H2O_aviation / plume_area_loss
     389  ! LU25, Eq. A3
     390  temp_205 = temp - 205.
     391  z_emit = 1106.6 * (rho_emit * 1e5)**(0.678 + 0.0116 * temp_205) &
     392          * EXP(- (0.0807 + 0.000428 * temp_205) * temp_205)
     393 
     394  ! U16, Eq. 4
     395  z_desc = SQRT(8. * circ_0_loss / RPI / N_Brunt_Vaisala_aviation )
     396 
     397  ! U16, Eq. 3, 10d, 10e, 10f, 10g and 10h
     398  z_delta = (AEI_contrails / 2.8e14)**(-0.18) * (1.7 * z_atm + 1.15 * z_emit) - 0.6 * z_desc
     399 
     400  ! U16, Eq. 9, 10a, 10b and 10c
     401  fsurv_contrails = 0.4 + 1.19 / RPI * ATAN(-1.35 + z_delta / 100.)
     402  fsurv_contrails = MIN(1., MAX(0., fsurv_contrails))
     403ELSEIF ( iflag_vortex_loss .EQ. 2 ) THEN
     404  ! LU25, Eq. A2
     405  z_atm = 607.46 * (icesat_ratio - 1.)**0.897 * (temp / 205.)**2.225
     406 
     407  ! water vapor emissions [kg/m flown], LU25, Eq. 2
     408  ! U16, Eq. 6
     409  rho_emit = fuel_per_m * EI_H2O_aviation / plume_area_loss
     410  ! LU25, Eq. A3
     411  temp_205 = temp - 205.
     412  z_emit = 1106.6 * (rho_emit * 1e5)**(0.678 + 0.0116 * temp_205) &
     413          * EXP(- (0.0807 + 0.000428 * temp_205) * temp_205)
     414 
     415  ! U16, Eq. 4
     416  z_desc = SQRT(8. * circ_0_loss / RPI / N_Brunt_Vaisala_aviation )
     417 
     418  ! ice crystals number concentration [#/m3], LU25, Eq. A1
     419  nice_init = Nice_per_m / plume_area_loss
     420  ! LU25, Eq. 9, 13d, 13e, 13f and 13g
     421  z_delta = (nice_init / nice_init_ref_loss)**(-0.16) * (1.27 * z_atm + 0.42 * z_emit) - 0.49 * z_desc
     422 
     423  ! LU25, Eq. 11, 12, 13a, 13b and 13c
     424  fsurv_contrails = 0.42 + 1.31 / RPI * ATAN(-1. + z_delta / 100.)
     425  fsurv_contrails = MIN(1., MAX(0., fsurv_contrails))
     426ENDIF
    412427
    413428Nice_per_m_init_contrails = Nice_per_m * fsurv_contrails
Note: See TracChangeset for help on using the changeset viewer.